ReactorEvent.cpp Source File

00001 
00007 #include <math.h>
00008 #include <iostream.h>
00009 #include <fstream.h>   // file I/O
00010 #include <iomanip.h>   // format manipulation
00011 #include <string>
00012 #include <stdio.h>
00013 #include "ReactorConstants.hh"
00014 #include "ReactorFortran.hh"
00015 #include "ReactorXC.hh"
00016 #include "ReactorEvent.hh"
00017 #include "MyCfSource.hh"
00018 #include "MuonPropagator.hh"
00019 
00020 
00021 ReactorEvent::ReactorEvent(int newNevents,
00022                            double newRmaxgen,
00023                            double newEmin,
00024                            double newEmax,
00025                            int newXCorder){
00026 
00027   Pdim = 7; 
00028   Rdim = 7;
00029 
00030   Nevents = newNevents;
00031   Rmaxgen_RE = newRmaxgen;
00032   Emin_RE = newEmin;
00033   Emax_RE = newEmax;
00034   XCorder_RE = newXCorder;
00035 
00036   /*
00037   cout << "ReactorEvent::ReactorEvent(" << Nevents << "," 
00038        << Rmaxgen_RE << "," << Emin << "," << Emax << "," 
00039        << XCorder_RE << ")" << endl;
00040   */
00041 
00042   int nmax = Ndim;
00043   int nnu  = 0;
00044   int npos = 0;
00045   int nele = 0;
00046   int nneu = 0;
00047   int ngam = 0;
00048 
00049   static double FracDone;
00050   static double FracDisplay;
00051   static double FracLastDisplay;
00052   static double NeventLastDisplay;
00053   static bool display;
00054 
00055   double* pnubar = new double[nmax*Pdim];
00056   double* pelectron = new double[nmax*Pdim];
00057   double* ppositron = new double[nmax*Pdim];
00058   double* pneutron = new double[nmax*Pdim];
00059   double* pgamma = new double[nmax*Pdim];
00060 
00061   double* pnu  = pnubar;
00062   double* pele = pelectron;
00063   double* ppos = ppositron;
00064   double* pneu = pneutron;
00065   double* pgam = pgamma;
00066 
00067   double* rnubar = new double[nmax*Rdim];
00068   double* relectron = new double[nmax*Rdim];
00069   double* rgenele = new double[nmax*Rdim];
00070   double* rpositron = new double[nmax*Rdim];
00071   double* rgenpos = new double[nmax*Rdim];
00072   double* rneutron = new double[nmax*Rdim];
00073   double* rgenneu = new double[nmax*Rdim];
00074   double* rgamma = new double[nmax*Rdim];
00075 
00076   double* rnu  = rnubar;
00077   double* rele = relectron;
00078   double* rge  = rgenele;
00079   double* rpos = rpositron;
00080   double* rgp  = rgenpos;
00081   double* rneu = rneutron;
00082   double* rgn  = rgenneu;
00083   double* rgam = rgamma;
00084 
00085   double Enu,z,phi;
00086   double R,zR,phiR;
00087 
00088   // setup the event
00089   fstream fin("src/event_setup.txt", ios::in);
00090   //cout << "setup file: src/event_setup.txt" << endl;
00091 
00092   char paramLine[256];
00093   std::string dummy;
00094   bool done = false;
00095   bool set = false;
00096   std::string static eventtype;
00097   std::string static pmtrad;
00098   double static Depth;
00099   bool static first = true;
00100 
00101   if(first){
00102     while(!done){
00103 
00104       // '#' specifies a comment line -- ignore it
00105       if(fin.peek() == '#'){
00106 
00107         fin.getline(paramLine, 256);
00108         continue;
00109       }
00110 
00111       // 'E'ND or 'e'nd specifies end of file -- set done flag to true
00112       if(fin.peek() == 'e' || fin.peek() == 'E'){
00113 
00114         done = true;
00115         continue;
00116       }
00117 
00118       // 'K'IND or 'k'ind specifies the event type
00119       if(fin.peek() == 'k' || fin.peek() == 'K'){
00120         
00121         if(set){
00122 
00123           fin.getline( paramLine, 256 );
00124           continue;
00125         }
00126         fin >> dummy;
00127         while(fin.peek() == ' ' || fin.peek() == '=') 
00128           fin.ignore(1);
00129 
00130         fin >> eventtype;
00131 
00132         if(eventtype == "reactor" || eventtype == "cf252" || 
00133            eventtype == "pmtrad"  || eventtype == "muon" || 
00134            eventtype == "all"     || eventtype == "veto"){
00135           set = true;
00136         }
00137         else{
00138           cout << "SETUP ERROR: unrecognized event type" << eventtype << endl;
00139         }
00140 
00141         continue;
00142       }
00143 
00144       // 'P'MTRAD or 'p'mtrad turns on/off radioactive decays in the PMTs
00145       if(fin.peek() == 'p' || fin.peek() == 'P'){
00146 
00147         fin >> dummy;
00148         while(fin.peek() == ' ' || fin.peek() == '=') 
00149           fin.ignore(1);
00150 
00151         fin >> pmtrad;
00152         continue;
00153       }
00154 
00155       // 'D'EPTH or 'd'epth sets the detector depth in mwe
00156       if(fin.peek() == 'd' || fin.peek() == 'D'){
00157 
00158         fin >> dummy;
00159         while(fin.peek() == ' ' || fin.peek() == '=') 
00160           fin.ignore(1);
00161 
00162         fin >> Depth;
00163         continue;
00164       }
00165 
00166       // skip n, o, r, t, a, b, m, and f
00167       if(fin.peek() == 'N' || fin.peek() == 'n' ||
00168          fin.peek() == 'O' || fin.peek() == 'o' ||
00169          fin.peek() == 'R' || fin.peek() == 'r' ||
00170          fin.peek() == 'T' || fin.peek() == 't' ||
00171          fin.peek() == 'A' || fin.peek() == 'a' ||
00172          fin.peek() == 'B' || fin.peek() == 'b' ||
00173          fin.peek() == 'M' || fin.peek() == 'm' ||
00174          fin.peek() == 'F' || fin.peek() == 'f'){
00175 
00176         fin.getline(paramLine, 256);
00177         continue;
00178       }
00179 
00180       // ignore extra whitespace
00181       if(fin.peek() == ' ' || fin.peek() == '\n'){
00182 
00183         fin.ignore(1);
00184         continue;
00185       }
00186     }
00187 
00188     // log file printout
00189     cout << "===============================" << endl;
00190     cout << "  " << eventtype << " type of events" << endl;
00191     cout << "  PMT radioactivity is " << pmtrad << endl;
00192     cout << "  Detector depth is " << Depth << " mwe" << endl;
00193     cout << "===============================" << endl;
00194 
00195     // zero the event counter
00196     Nevent = 0;
00197 
00198     first = false;
00199   }
00200   eventType = eventtype;
00201   pmtRad = pmtrad;
00202 
00203   // report fraction of processed events
00204   display = 0;
00205   if (Nevent==0) {
00206     // Set up fractions for progress checking
00207     FracDisplay = 0.01;
00208     NeventLastDisplay = 0.;
00209     FracLastDisplay = 0.;
00210     FracDone = 0.;
00211     display = 1;
00212   }
00213   //
00214   // report every 1% of processed events
00215   //
00216   FracDone = (double) Nevent / (double) Nevents;
00217   if (FracDone >= 0.){
00218     if (display || (FracDone-FracLastDisplay) >= FracDisplay){
00219       cout << Nevent << " events (" << 100.*FracDone 
00220            << "%) processed" << endl;
00221       FracLastDisplay += FracDisplay;
00222     }
00223   }
00224   else {
00225     if (display || Nevent-NeventLastDisplay >= FracDisplay*NeventLastDisplay){
00226       cout << Nevent << " events processed" << endl;
00227     }
00228   }
00229   //
00230   // count the event
00231   //
00232   Nevent++;
00233 
00234   //  if(Nevent%1000 == 0) cout << Nevent << " events processed" << endl;
00235   //
00236   // if all events are requested, mix the 'reactor' and 'cf252' 
00237   // events together with 90/10% Cf252/reactor events
00238   //
00239   if(eventType == "all"){
00240     const int len = 1;
00241     float r[len];
00242     ranlux_(r,len);
00243 
00244     if(r[0] < 0.90){
00245       eventType = "cf252";
00246     }
00247     else{
00248       eventType = "reactor";
00249     }
00250   }
00251   //cout << eventType << " type of events" << endl;
00252   //cout << "PMT radiation " << pmtRad << endl;
00253 
00254   if(eventType == "reactor"){
00255 
00256     XCWeight = 1;
00257 
00258     // inverse beta decay reaction
00259     double XCmax = InverseBeta(Emax_RE,-1);
00260     double FluxMax = RMPflux(Emin_RE);
00261 
00262     bool passed=0;
00263 
00264     while(!passed){
00265       const int len = 7;
00266       float r[len];
00267       ranlux_(r,len);
00268       Enu = Emin_RE+(Emax_RE-Emin_RE)*r[0];
00269       z = -1+2*r[1];
00270       phi = 2*PI*r[2];
00271 
00272       R = Rmaxgen_RE*exp(log(r[3])/3);
00273       zR = -1+2*r[4];
00274       phiR = 2*PI*r[5];
00275       //      cout << R << " " << Rmaxgen_RE << " " << r[3] << endl;
00276 
00277       float XCtest = XCmax*FluxMax*r[6];
00278       XCWeight = InverseBeta(Enu,z);
00279       FluxWeight=RMPflux(Enu);
00280       passed= XCWeight*FluxWeight>XCtest;
00281     }
00282     //      cout << "passed " << R << " " << Rmaxgen_RE << " "<< endl;
00283 
00284   
00285     double E0 = Enu - DELTA;
00286     double p0 = sqrt(E0*E0-MELECTRON*MELECTRON);
00287     double v0 = p0/E0;
00288     double Ysquared = (DELTA*DELTA-MELECTRON*MELECTRON)/2;
00289     double E1 = E0*(1-Enu/MPROTON*(1-v0*z))-Ysquared/MPROTON;
00290     double p1 = sqrt(E1*E1-MELECTRON*MELECTRON);
00291 
00292     *pnu = Enu; pnu++;
00293     *pnu = 0; pnu++;
00294     *pnu = 0; pnu++;
00295     *pnu = Enu; pnu++;
00296     *pnu = Enu; pnu++;
00297     *pnu = Enu; pnu++;
00298     *pnu = 1.; pnu++; // inverse beta decay process Id = 1
00299 
00300     *ppos = E1; ppos++;
00301     *ppos = p1*sqrt(1-z*z)*cos(phi); ppos++;
00302     *ppos = p1*sqrt(1-z*z)*sin(phi); ppos++;
00303     *ppos = p1*z; ppos++;
00304     *ppos = E1-MELECTRON; ppos++;
00305     *ppos = p1; ppos++;
00306     *ppos = 1.; ppos++;
00307 
00308     *pneu = Enu+MPROTON-E1; pneu++;
00309     *pneu = -p1*sqrt(1-z*z)*cos(phi); pneu++;
00310     *pneu = -p1*sqrt(1-z*z)*sin(phi); pneu++;
00311     *pneu = Enu-p1*z; pneu++;
00312     *pneu = Enu-E1-DELTA; pneu++;
00313     *pneu = sqrt((Enu+MPROTON-E1)*(Enu+MPROTON-E1)-
00314                          (MPROTON+DELTA)*(MPROTON+DELTA)); pneu++;
00315     *pneu = 1.; pneu++;
00316 
00317     *rnu = 0; rnu++;
00318     *rnu = R*sqrt(1-zR*zR)*cos(phiR); rnu++;
00319     *rnu = R*sqrt(1-zR*zR)*sin(phiR); rnu++;
00320     *rnu = R*zR; rnu++;
00321     *rnu = R; rnu++;
00322     *rnu = zR; rnu++;
00323     *rnu = phiR; rnu++;
00324 
00325     *rpos = 0; rpos++;
00326     *rpos = R*sqrt(1-zR*zR)*cos(phiR); rpos++;
00327     *rpos = R*sqrt(1-zR*zR)*sin(phiR); rpos++;
00328     *rpos = R*zR; rpos++;
00329     *rpos = R; rpos++;
00330     *rpos = zR; rpos++;
00331     *rpos = phiR; rpos++;
00332 
00333     *rgp = 0; rgp++;
00334     *rgp = R*sqrt(1-zR*zR)*cos(phiR); rgp++;
00335     *rgp = R*sqrt(1-zR*zR)*sin(phiR); rgp++;
00336     *rgp = R*zR; rgp++;
00337     *rgp = R; rgp++;
00338     *rgp = zR; rgp++;
00339     *rgp = phiR; rgp++;
00340 
00341     *rneu = 0; rneu++;
00342     *rneu = R*sqrt(1-zR*zR)*cos(phiR); rneu++;
00343     *rneu = R*sqrt(1-zR*zR)*sin(phiR); rneu++;
00344     *rneu = R*zR; rneu++;
00345     *rneu = R; rneu++;
00346     *rneu = zR; rneu++; 
00347     *rneu = phiR; rneu++;
00348 
00349     *rgn = 0; rgn++;
00350     *rgn = R*sqrt(1-zR*zR)*cos(phiR); rgn++;
00351     *rgn = R*sqrt(1-zR*zR)*sin(phiR); rgn++;
00352     *rgn = R*zR; rgn++;
00353     *rgn = R; rgn++;
00354     *rgn = zR; rgn++;
00355     *rgn = phiR; rgn++;
00356 
00357     // count the nubar, positron, and neutron
00358     nnu++;
00359     npos++;
00360     nneu++;
00361 
00362   } // done with reactor neutrino inverse-beta reaction
00363 
00364   // Cf source reaction
00365   else if(eventType == "cf252"){
00366     XCWeight = 1;
00367     FluxWeight = 1;
00368 
00369     int Isotope = 252;
00370     MyCfSource CfSource(Isotope,RMAXGEN);
00371 
00372     int Nsources = CfSource.GetNumCfSources();
00373     //cout << "Number of Cf252 sources = " << Nsources << endl;
00374 
00375     for(int n=0;n<Nsources;n++){
00376 
00377       R = CfSource.GetCfVertexR(n);
00378       zR = CfSource.GetCfVertexCosTheta(n);
00379       phiR = CfSource.GetCfVertexPhi(n);
00380 
00381       // loop over the neutrons
00382       for(int nn=0;nn<CfSource.GetNumNeutron(n);nn++){
00383 
00384         double neutronKE = CfSource.GetCfNeutronEnergy(n,nn);
00385         //cout << "   Cf252 source " << n << ", neutron " << nn 
00386         //     << " energy = " << neutronKE << endl;
00387 
00388         double Tneu = CfSource.GetCfNeutronTime(n,nn);
00389         double Mneutron = MPROTON+DELTA;
00390         double neutronE = neutronKE+Mneutron;
00391 
00392         // pick a direction for neutron momentum
00393         const int len = 2;
00394         float r[len];
00395         ranlux_(r,len);
00396         z = -1+2*r[0];
00397         phi = 2*PI*r[1];
00398         double x = sqrt(1-z*z)*cos(phi);
00399         double y = sqrt(1-z*z)*sin(phi);
00400 
00401         *pneu = neutronE; pneu++;
00402         *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*x; pneu++;
00403         *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*y; pneu++;
00404         *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*z; pneu++;
00405         *pneu = neutronKE; pneu++;
00406         *pneu = sqrt(neutronE*neutronE-Mneutron*Mneutron); pneu++;
00407         *pneu = 2.; pneu++; // cf252 source process Id = 2
00408 
00409         *rneu = Tneu; rneu++;
00410         *rneu = R*sqrt(1-zR*zR)*cos(phiR); rneu++;
00411         *rneu = R*sqrt(1-zR*zR)*sin(phiR); rneu++;
00412         *rneu = R*zR; rneu++;
00413         *rneu = R; rneu++;
00414         *rneu = zR; rneu++;
00415         *rneu = phiR; rneu++;
00416 
00417         *rgn = Tneu; rgn++;
00418         *rgn = R*sqrt(1-zR*zR)*cos(phiR); rgn++;
00419         *rgn = R*sqrt(1-zR*zR)*sin(phiR); rgn++;
00420         *rgn = R*zR; rgn++;
00421         *rgn = R; rgn++;
00422         *rgn = zR; rgn++;
00423         *rgn = phiR; rgn++;
00424       }
00425       // count the neutrons
00426       nneu += CfSource.GetNumNeutron(n);      
00427 
00428       // loop over the prompt gammas (these are real photons and 
00429       // need to be compton scattered in the scint)
00430       //
00431       for(int nn=0;nn<CfSource.GetNumGamma(n);nn++){
00432 
00433         double gammaKE = CfSource.GetCfGammaEnergy(n,nn);
00434         //cout << "   Cf252 source " << n << ", gamma " << nn 
00435         //     << " energy = " << gammaKE << endl;
00436 
00437         double Tgam = CfSource.GetCfGammaTime(n,nn);
00438 
00439         // compton scatter the gammas
00440         //
00441         const int len = 5;
00442         float rv[len];
00443         ranlux_(rv,len);
00444         double cosg = -1+2*rv[0];
00445         double phig = 2*PI*rv[1];
00446         double t[3];
00447         t[0] = cos(phig)*sqrt(1-cosg*cosg);
00448         t[1] = sin(phig)*sqrt(1-cosg*cosg);
00449         t[2] = cosg;
00450         //for(int i=0;i<3;i++) cout << "t["<<i<<"] " << t[i] << " ";
00451         //cout << endl;
00452 
00453         double range;
00454         std::string material;
00455         if(R < Rmaxgen_RE)
00456           material = "scintillator";
00457         else
00458           material = "oil";
00459         range = GammaComptonRange(material,gammaKE);
00460         range*=log(1/rv[2]);
00461         //cout << "range " << range << endl;
00462         //
00463         // get the new x,y,z positions from source R,zR,phiR
00464         //
00465         double xpos = R*sqrt(1-zR*zR)*cos(phiR) + range*t[0];
00466         double ypos = R*sqrt(1-zR*zR)*sin(phiR) + range*t[1];
00467         double zpos = R*zR + range*t[2];
00468         //
00469         // now convert back to R,zR,phiR
00470         //
00471         R = sqrt(xpos*xpos + ypos*ypos + zpos*zpos);
00472         zR = zpos/R;
00473         phiR = atan2(ypos,xpos);
00474         if(phiR < 0) phiR += 2*PI;
00475         //
00476         // pick a direction for gamma momentum
00477         //
00478         float r[len];
00479         ranlux_(r,len);
00480         z = -1+2*rv[3];
00481         phi = 2*PI*rv[4];
00482         double x = sqrt(1-z*z)*cos(phi);
00483         double y = sqrt(1-z*z)*sin(phi);
00484 
00485         *pgam = gammaKE; pgam++;
00486         *pgam = gammaKE*x; pgam++;
00487         *pgam = gammaKE*y; pgam++;
00488         *pgam = gammaKE*z; pgam++;
00489         *pgam = gammaKE; pgam++;
00490         *pgam = gammaKE; pgam++;
00491         *pgam = 2.; pgam++; // cf252 source process Id = 2
00492 
00493         *rgam = Tgam; rgam++;
00494         *rgam = xpos; rgam++;
00495         *rgam = ypos; rgam++;
00496         *rgam = zpos; rgam++;
00497         *rgam = R; rgam++;
00498         *rgam = zR; rgam++;
00499         *rgam = phiR; rgam++;
00500       }
00501       //
00502       ngam += CfSource.GetNumGamma(n);
00503 
00504     } // done looping over sources
00505 
00506   } // done with cf source neutron production
00507 
00508   // PMT radiation
00509   else if(eventType == "pmtrad"){
00510 
00511     XCWeight = 1;
00512     FluxWeight = 1;
00513 
00514     // nothing happens in these events but photons, generated in 
00515     // ReactorDetector::GenerateGammas with pmtRad on
00516     if(pmtRad == "off"){
00517       cout << "SETUP ERROR: PMT radiation set " << pmtRad << " in a " 
00518            << eventType << " type event" << endl;
00519       pmtRad = "on";
00520       cout << "Forcing PMT radiation " << pmtRad << endl;
00521     }
00522   }
00523 
00524   // muons
00525   else if(eventType == "muon"){
00526 
00527     //cout << eventType << endl;
00528 
00529     // create the energy and cosine of zenith angle (theta) arrays
00530     bool static first = true;
00531     double static Energy[Mdim],CosTheta[Mdim];
00532     double dummy;
00533 
00534     if(first){
00535 
00536       // read in the muon energy in GeV and angle
00537       if(Depth == 300.){
00538         fstream fin("data/Totb300a.dat", ios::in);
00539         cout << "setup file: data/Totb300a.dat" << endl;
00540 
00541         for(int i=0;i<Mdim;i++){
00542           fin >> dummy >> Energy[i] >> dummy >> dummy >> dummy >> CosTheta[i];
00543           //cout << i << "  " << Energy[i] << "  " << CosTheta[i] << endl;
00544         }
00545         printf("Read %d input records\n",Mdim);
00546       }
00547       else if(Depth == 450.){
00548         fstream fin("data/Totb450a.dat", ios::in);
00549         cout << "setup file: data/Totb450a.dat" << endl;
00550 
00551         for(int i=0;i<Mdim;i++){
00552           fin >> dummy >> Energy[i] >> dummy >> dummy >> dummy >> CosTheta[i];
00553           //cout << i << "  " << Energy[i] << "  " << CosTheta[i] << endl;
00554         }
00555         printf("Read %d input records\n",Mdim);
00556       }
00557       else if(Depth == 500.){
00558         fstream fin("data/Totb500a.dat", ios::in);
00559         cout << "setup file: data/Totb500a.dat" << endl;
00560 
00561         for(int i=0;i<Mdim;i++){
00562           fin >> dummy >> Energy[i] >> dummy >> dummy >> dummy >> CosTheta[i];
00563           //cout << i << "  " << Energy[i] << "  " << CosTheta[i] << endl;
00564         }
00565         printf("Read %d input records\n",Mdim);
00566       }
00567       else{
00568         cout << "ERROR: no muon data for depth of " 
00569              << Depth << " mwe" << endl;
00570         return;
00571       }
00572   
00573       first = false;
00574     }
00575 
00576     /* get the muon energy and cosine of the zenith angle
00577        the spectra in the Tot dat files are random so start 
00578        with the first entry and cycle through the file
00579     */
00580     static int index = 0;
00581     double energy = 0.;
00582     double cosg = 0.;
00583     //
00584     if(index < Mdim){
00585       //
00586       // muon energy in MeV
00587       energy = Energy[index]*1000.;
00588       cosg = CosTheta[index];
00589       //cout << " muon energy " << energy 
00590       //   << " MeV, cos(theta) " << cosg << endl;
00591       index++;
00592     }
00593     else{
00594       cout << "  WARNING: end of muon input file reached at index " << index 
00595            << ".  Starting over." << endl;
00596       index = 0;
00597     }
00598     //
00599     // call the propagator, which returns particle information
00600     //
00601     MuonPropagator Muon(RMAXGEN,Depth,energy,cosg);
00602 
00603     XCWeight = Muon.GetWeight();
00604     FluxWeight = 1;
00605     //
00606     // loop over neutrons
00607     //
00608     nneu = Muon.GetNumNeutron();
00609     for(int n=0;n<nneu;n++){
00610 
00611       float neutronKE = Muon.GetNeutronKE(n);
00612       //cout << "  neutron " << n << " energy = " << neutronKE << endl;
00613 
00614       float Mneutron = MPROTON+DELTA;
00615       float neutronE = neutronKE+Mneutron;
00616 
00617       // pick a direction for neutron momentum
00618       const int len = 2;
00619       float rv[len];
00620       ranlux_(rv,len);
00621       z = -1+2*rv[0];
00622       phi = 2*PI*rv[1];
00623       double x = sqrt(1-z*z)*cos(phi);
00624       double y = sqrt(1-z*z)*sin(phi);
00625 
00626       *pneu = neutronE; pneu++;
00627       *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*x; pneu++;
00628       *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*y; pneu++;
00629       *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*z; pneu++;
00630       *pneu = neutronKE; pneu++;
00631       *pneu = sqrt(neutronE*neutronE-Mneutron*Mneutron); pneu++;
00632       *pneu = 3.; pneu++; // muon process Id = 3
00633 
00634       double neuxyz[3];
00635       for(int i=0;i<3;i++){
00636         neuxyz[i] = *(Muon.GetNeutronVertexXYZ()+(i+n));
00637         //cout << "  neuxyz[" << i << "] " << neuxyz[i] << endl;
00638       }
00639       double rr = sqrt(neuxyz[0]*neuxyz[0] + 
00640                        neuxyz[1]*neuxyz[1] + 
00641                        neuxyz[2]*neuxyz[2]);
00642       double rphi = atan2(neuxyz[1],neuxyz[0]);
00643       if(rphi<0) rphi += 2*PI;
00644 
00645       double time = Muon.GetNeutronTime(n);
00646       //cout << " n time " << time << endl;
00647 
00648       *rneu = time; rneu++;
00649       *rneu = neuxyz[0]; rneu++;
00650       *rneu = neuxyz[1]; rneu++;
00651       *rneu = neuxyz[2]; rneu++;
00652       *rneu = rr; rneu++;
00653       *rneu = neuxyz[2]/rr; rneu++;
00654       *rneu = rphi; rneu++;
00655 
00656       *rgn = time; rgn++;
00657       *rgn = neuxyz[0]; rgn++;
00658       *rgn = neuxyz[1]; rgn++;
00659       *rgn = neuxyz[2]; rgn++;
00660       *rgn = rr; rgn++;
00661       *rgn = neuxyz[2]/rr; rgn++;
00662       *rgn = rphi; rgn++;
00663     }
00664     //
00665     // loop over electrons
00666     //
00667     nele = Muon.GetNumElectron();
00668     for(int n=0;n<nele;n++){
00669 
00670       float electronKE = Muon.GetElectronKE(n);
00671       //cout << "  electron " << n << " KE = " << electronKE << endl;
00672 
00673       float electronE = electronKE+MELECTRON;
00674       //cout << "  electron " << n << " energy = " << electronE << endl;
00675 
00676       // pick a direction for electron momentum
00677       const int len = 2;
00678       float r[len];
00679       ranlux_(r,len);
00680       z = -1+2*r[0];
00681       phi = 2*PI*r[1];
00682       double x = sqrt(1-z*z)*cos(phi);
00683       double y = sqrt(1-z*z)*sin(phi);
00684 
00685       *pele = electronE; pele++;
00686       *pele = sqrt((electronE*electronE)-(MELECTRON*MELECTRON))*x; pele++;
00687       *pele = sqrt((electronE*electronE)-(MELECTRON*MELECTRON))*y; pele++;
00688       *pele = sqrt((electronE*electronE)-(MELECTRON*MELECTRON))*z; pele++;
00689       *pele = electronKE; pele++;
00690       *pele = sqrt(electronE*electronE-MELECTRON*MELECTRON); pele++;
00691       *pele = 3.; pele++; // muon process Id = 3
00692 
00693       double elexyz[3];
00694       for(int i=0;i<3;i++){
00695         elexyz[i] = *(Muon.GetElectronVertexXYZ()+(i+3*n));
00696         //cout << "  elexyz[" << i << "] " << elexyz[i] << endl;
00697       }
00698       double rr = sqrt(elexyz[0]*elexyz[0] + 
00699                        elexyz[1]*elexyz[1] + 
00700                        elexyz[2]*elexyz[2]);
00701       double rphi = atan2(elexyz[1],elexyz[0]);
00702       if(rphi<0) rphi += 2*PI;
00703 
00704       double time = Muon.GetElectronTime(n);
00705       //cout << " e time " << time << endl;
00706 
00707       *rele = time; rele++;
00708       *rele = elexyz[0]; rele++;
00709       *rele = elexyz[1]; rele++;
00710       *rele = elexyz[2]; rele++;
00711       *rele = rr; rele++;
00712       *rele = elexyz[2]/rr; rele++;
00713       *rele = rphi; rele++;
00714 
00715       *rge = time; rge++;
00716       *rge = elexyz[0]; rge++;
00717       *rge = elexyz[1]; rge++;
00718       *rge = elexyz[2]; rge++;
00719       *rge = rr; rge++;
00720       *rge = elexyz[2]/rr; rge++;
00721       *rge = rphi; rge++;
00722     }
00723     //
00724     // loop over photons (these are 'virtual' to contain the energy
00725     // loss of the muon and therefore are not compton scattered)
00726     //
00727     ngam = Muon.GetNumGamma();
00728     for(int n=0;n<ngam;n++){
00729 
00730       float gammaKE = Muon.GetGammaKE(n);
00731       //cout << "  gamma " << n << " energy = " << gammaKE << endl;
00732 
00733       // pick a direction for gamma momentum
00734       const int len = 2;
00735       float r[len];
00736       ranlux_(r,len);
00737       z = -1+2*r[0];
00738       phi = 2*PI*r[1];
00739       double x = sqrt(1-z*z)*cos(phi);
00740       double y = sqrt(1-z*z)*sin(phi);
00741 
00742       *pgam = gammaKE; pgam++;
00743       *pgam = gammaKE*x; pgam++;
00744       *pgam = gammaKE*y; pgam++;
00745       *pgam = gammaKE*z; pgam++;
00746       *pgam = gammaKE; pgam++;
00747       *pgam = gammaKE; pgam++;
00748       *pgam = 3.; pgam++; // muon process Id = 3
00749 
00750       double gamxyz[3];
00751       for(int i=0;i<3;i++){
00752         gamxyz[i] = *(Muon.GetGammaVertexXYZ()+(i+3*n));
00753         //cout << "  gamxyz[" << i << "] " << gamxyz[i] << endl;
00754       }
00755       double rr = sqrt(gamxyz[0]*gamxyz[0] + 
00756                        gamxyz[1]*gamxyz[1] + 
00757                        gamxyz[2]*gamxyz[2]);
00758       double rphi = atan2(gamxyz[1],gamxyz[0]);
00759       if(rphi<0) rphi += 2*PI;
00760 
00761       double time = Muon.GetGammaTime(n);
00762 
00763       *rgam = time; rgam++;
00764       *rgam = gamxyz[0]; rgam++;
00765       *rgam = gamxyz[1]; rgam++;
00766       *rgam = gamxyz[2]; rgam++;
00767       *rgam = rr; rgam++;
00768       *rgam = gamxyz[2]/rr; rgam++;
00769       *rgam = rphi; rgam++;
00770     }
00771   }
00772 
00773   // veto
00774   else if(eventType == "veto"){
00775 
00776     //cout << eventType << endl;
00777 
00778     // read in the veto file
00779     fstream fveto("data/testveto.ascii", ios::in);
00780     //cout << "setup file: data/testveto.ascii" << endl;
00781 
00782     // convert the event number to a character string
00783     int localeventsize = 0;
00784     if(Nevent < 10)
00785       localeventsize = 1;
00786     else if(Nevent > 9 && Nevent < 100)
00787       localeventsize = 2;
00788     else if(Nevent > 99 && Nevent < 1000)
00789       localeventsize = 3;
00790     else if(Nevent > 999 && Nevent < 10000)
00791       localeventsize = 4;
00792     else if(Nevent > 9999 && Nevent < 100000)
00793       localeventsize = 5;
00794     else if(Nevent > 99999 && Nevent < 1000000)
00795       localeventsize = 6;
00796     else if(Nevent > 999999 && Nevent < 10000000)
00797       localeventsize = 7;
00798     else 
00799       localeventsize = 0;
00800 
00801     char *localevent = new char[localeventsize];
00802 //    char localevent[localeventsize];
00803     sprintf(localevent,"%i",Nevent);
00804     //cout << "localevent " << localevent << endl;
00805 
00806     // counter
00807     int Npart = 0;
00808     // number of variables per particle from veto: id,x,y,z,vx,vy,vz,ke
00809     int Vdim = 8;
00810 
00811     double* VetoData = new double[Ndim*Vdim];
00812 
00813     bool edone = false;
00814     while(!edone){
00815       //
00816       // get the veto event with the matching event number
00817       //
00818       std::string vetoevent;
00819       getline(fveto,vetoevent);
00820       //cout << " vetoevent " << vetoevent << endl;
00821       //
00822       if(vetoevent == localevent){
00823 
00824         bool pdone = false;
00825         while(!pdone){
00826           //
00827           // get the particle(s) information
00828           //
00829           if(fveto.peek() == '0'){
00830             //
00831             // end of particle list
00832             //
00833             edone = true;
00834             pdone = true;
00835           }
00836           else if(fveto.peek() == ' '){
00837             //
00838             // a particle, info is: id,x,y,z,vx,vy,vz,ke
00839             //
00840             for(int i=0;i<Vdim;i++) fveto >> VetoData[i+Vdim*Npart];
00841             //
00842             // close out any remaining whitespace
00843             //
00844             std::string dummy;
00845             getline(fveto,dummy);
00846             //
00847             Npart++;
00848           }
00849         } // close loop over particles in veto event
00850       }
00851     } // close loop over veto events
00852 
00853     /*
00854     cout << Npart << " veto particle(s)" << endl;
00855     for(int n=0;n<Npart*Vdim;n++){
00856       cout << " VetoData[" << n << "] " << VetoData[n] << endl;
00857     }
00858     */
00859 
00860     //
00861     // sort particles into R and P arrays
00862     //
00863     for(int n=0;n<Npart;n++){
00864       //
00865       // PID follows the PDG Monte Carlo numbering scheme
00866       //
00867       int particleid = int(VetoData[n*Vdim]);
00868       //
00869       if(particleid == 11){
00870         //
00871         // electron
00872         //
00873         nele++;
00874 
00875         // electron KE and total E in MeV
00876         double electronKE = VetoData[n*Vdim+7]*1000.;
00877         double electronE = electronKE+MELECTRON;
00878 
00879         *pele = electronE; pele++;
00880         *pele = VetoData[n*Vdim+4]*1000.; pele++;
00881         *pele = VetoData[n*Vdim+5]*1000.; pele++;
00882         *pele = VetoData[n*Vdim+6]*1000.; pele++;
00883         *pele = electronKE; pele++;
00884         *pele = sqrt(electronE*electronE-MELECTRON*MELECTRON); pele++;
00885         *pele = 4.; pele++; // veto process Id = 4
00886 
00887         double elexyz[3];
00888         for(int i=0;i<3;i++){
00889           elexyz[i] = VetoData[n*Vdim+(i+1)];
00890           //cout << "  elexyz[" << i << "] " << elexyz[i] << endl;
00891         }
00892         double rr = sqrt(elexyz[0]*elexyz[0] + 
00893                          elexyz[1]*elexyz[1] + 
00894                          elexyz[2]*elexyz[2]);
00895         double rphi = atan2(elexyz[1],elexyz[0]);
00896         if(rphi<0) rphi += 2*PI;
00897         
00898         double time = 0.;
00899 
00900         *rele = time; rele++;
00901         *rele = elexyz[0]; rele++;
00902         *rele = elexyz[1]; rele++;
00903         *rele = elexyz[2]; rele++;
00904         *rele = rr; rele++;
00905         *rele = elexyz[2]/rr; rele++;
00906         *rele = rphi; rele++;
00907 
00908         *rge = time; rge++;
00909         *rge = elexyz[0]; rge++;
00910         *rge = elexyz[1]; rge++;
00911         *rge = elexyz[2]; rge++;
00912         *rge = rr; rge++;
00913         *rge = elexyz[2]/rr; rge++;
00914         *rge = rphi; rge++;
00915       }
00916       else if(particleid == -11){
00917         //
00918         // positron
00919         //
00920         npos++;
00921 
00922         // positron KE and total E in MeV
00923         double positronKE = VetoData[n*Vdim+7]*1000.;
00924         double positronE = positronKE+MELECTRON;
00925 
00926         *ppos = positronE; ppos++;
00927         *ppos = VetoData[n*Vdim+4]*1000.; ppos++;
00928         *ppos = VetoData[n*Vdim+5]*1000.; ppos++;
00929         *ppos = VetoData[n*Vdim+6]*1000.; ppos++;
00930         *ppos = positronKE; ppos++;
00931         *ppos = sqrt(positronE*positronE-MELECTRON*MELECTRON); ppos++;
00932         *ppos = 4.; ppos++; // veto process Id = 4
00933 
00934         double posxyz[3];
00935         for(int i=0;i<3;i++){
00936           posxyz[i] = VetoData[n*Vdim+(i+1)];
00937           //cout << "  posxyz[" << i << "] " << posxyz[i] << endl;
00938         }
00939         double rr = sqrt(posxyz[0]*posxyz[0] + 
00940                          posxyz[1]*posxyz[1] + 
00941                          posxyz[2]*posxyz[2]);
00942         double rphi = atan2(posxyz[1],posxyz[0]);
00943         if(rphi<0) rphi += 2*PI;
00944         
00945         double time = 0.;
00946 
00947         *rpos = time; rpos++;
00948         *rpos = posxyz[0]; rpos++;
00949         *rpos = posxyz[1]; rpos++;
00950         *rpos = posxyz[2]; rpos++;
00951         *rpos = rr; rpos++;
00952         *rpos = posxyz[2]/rr; rpos++;
00953         *rpos = rphi; rpos++;
00954 
00955         *rgp = time; rgp++;
00956         *rgp = posxyz[0]; rgp++;
00957         *rgp = posxyz[1]; rgp++;
00958         *rgp = posxyz[2]; rgp++;
00959         *rgp = rr; rgp++;
00960         *rgp = posxyz[2]/rr; rgp++;
00961         *rgp = rphi; rgp++;
00962       }
00963       else if(particleid == 13){
00964         //
00965         // muon
00966         //
00967 
00968         // muon KE in MeV
00969         double muonKE = VetoData[n*Vdim+7]*1000.;
00970         double x = VetoData[n*Vdim+1];
00971         double y = VetoData[n*Vdim+2];
00972         double z = VetoData[n*Vdim+3];
00973         double px = VetoData[n*Vdim+4];
00974         double py = VetoData[n*Vdim+5];
00975         double pz = VetoData[n*Vdim+6];
00976 
00977         MuonPropagator Muon(0,Rmaxgen_RE,muonKE,x,y,z,px,py,pz);
00978 
00979         // loop over photons (these are 'virtual' to contain the energy
00980         // loss of the muon and therefore are not compton scattered)
00981         //
00982         ngam = Muon.GetNumGamma();
00983         for(int n=0;n<ngam;n++){
00984 
00985           float gammaKE = Muon.GetGammaKE(n);
00986           //cout << "  gamma " << n << " energy = " << gammaKE << endl;
00987 
00988           // pick a direction for gamma momentum
00989           const int len = 2;
00990           float r[len];
00991           ranlux_(r,len);
00992           z = -1+2*r[0];
00993           phi = 2*PI*r[1];
00994           double x = sqrt(1-z*z)*cos(phi);
00995           double y = sqrt(1-z*z)*sin(phi);
00996 
00997           *pgam = gammaKE; pgam++;
00998           *pgam = gammaKE*x; pgam++;
00999           *pgam = gammaKE*y; pgam++;
01000           *pgam = gammaKE*z; pgam++;
01001           *pgam = gammaKE; pgam++;
01002           *pgam = gammaKE; pgam++;
01003           *pgam = 4.; pgam++; // veto process Id = 4
01004 
01005           double gamxyz[3];
01006           for(int i=0;i<3;i++){
01007             gamxyz[i] = *(Muon.GetGammaVertexXYZ()+(i+3*n));
01008             //cout << "  gamxyz[" << i << "] " << gamxyz[i] << endl;
01009           }
01010           double rr = sqrt(gamxyz[0]*gamxyz[0] + 
01011                            gamxyz[1]*gamxyz[1] + 
01012                            gamxyz[2]*gamxyz[2]);
01013           double rphi = atan2(gamxyz[1],gamxyz[0]);
01014           if(rphi<0) rphi += 2*PI;
01015 
01016           double time = Muon.GetGammaTime(n);
01017 
01018           *rgam = time; rgam++;
01019           *rgam = gamxyz[0]; rgam++;
01020           *rgam = gamxyz[1]; rgam++;
01021           *rgam = gamxyz[2]; rgam++;
01022           *rgam = rr; rgam++;
01023           *rgam = gamxyz[2]/rr; rgam++;
01024           *rgam = rphi; rgam++;
01025         }
01026       }
01027       else if(particleid == 2112){
01028         //
01029         // neutron
01030         //
01031         nneu++;
01032 
01033         // neutron KE and total E in MeV
01034         double neutronKE = VetoData[n*Vdim+7]*1000.;
01035         double Mneutron = MPROTON+DELTA;
01036         double neutronE = neutronKE+Mneutron;
01037 
01038         *pneu = neutronE; pneu++;
01039         *pneu = VetoData[n*Vdim+4]*1000.; pneu++;
01040         *pneu = VetoData[n*Vdim+5]*1000.; pneu++;
01041         *pneu = VetoData[n*Vdim+6]*1000.; pneu++;
01042         *pneu = neutronKE; pneu++;
01043         *pneu = sqrt(neutronE*neutronE-Mneutron*Mneutron); pneu++;
01044         *pneu = 4.; pneu++; // veto process Id = 4
01045 
01046         double neuxyz[3];
01047         for(int i=0;i<3;i++){
01048           neuxyz[i] = VetoData[n*Vdim+(i+1)];
01049           //cout << "  neuxyz[" << i << "] " << neuxyz[i] << endl;
01050         }
01051         double rr = sqrt(neuxyz[0]*neuxyz[0] + 
01052                          neuxyz[1]*neuxyz[1] + 
01053                          neuxyz[2]*neuxyz[2]);
01054         double rphi = atan2(neuxyz[1],neuxyz[0]);
01055         if(rphi<0) rphi += 2*PI;
01056         
01057         double time = 0.;
01058 
01059         *rneu = time; rneu++;
01060         *rneu = neuxyz[0]; rneu++;
01061         *rneu = neuxyz[1]; rneu++;
01062         *rneu = neuxyz[2]; rneu++;
01063         *rneu = rr; rneu++;
01064         *rneu = neuxyz[2]/rr; rneu++;
01065         *rneu = rphi; rneu++;
01066 
01067         *rgn = time; rgn++;
01068         *rgn = neuxyz[0]; rgn++;
01069         *rgn = neuxyz[1]; rgn++;
01070         *rgn = neuxyz[2]; rgn++;
01071         *rgn = rr; rgn++;
01072         *rgn = neuxyz[2]/rr; rgn++;
01073         *rgn = rphi; rgn++;
01074       }
01075       else{
01076         cout << "ERROR: unknown particle ID from veto: " << particleid << endl;
01077       }
01078     } // end loop over particles
01079     delete[] localevent;
01080   }
01081   else{
01082     cout << "ERROR: unknown event type: " << eventType << endl;
01083   }
01084   //
01085   // transfer temp data to storage
01086   //
01087   Nnubar    = nnu;
01088   Nelectron = nele;
01089   Npositron = npos;
01090   Nneutron  = nneu;
01091   Ngamma    = ngam;
01092 
01093   Pnubar    = new double[Pdim*nnu];
01094   Pelectron = new double[Pdim*nele];
01095   Ppositron = new double[Pdim*npos];
01096   Pneutron  = new double[Pdim*nneu];
01097   Pgamma    = new double[Pdim*ngam];
01098 
01099   Rnubar    = new double[Rdim*nnu];
01100   Relectron = new double[Rdim*nele];
01101   Rgenele   = new double[Rdim*nele];
01102   Rpositron = new double[Rdim*npos];
01103   Rgenpos   = new double[Rdim*npos];
01104   Rneutron  = new double[Rdim*nneu];
01105   Rgenneu   = new double[Rdim*nneu];
01106   Rgamma    = new double[Rdim*ngam];
01107 
01108   double* ptr;
01109 
01110   ptr = Pnubar+Pdim*Nnubar;
01111   while(ptr>Pnubar)*--ptr=*--pnu;
01112 
01113   ptr = Pelectron+Pdim*Nelectron;
01114   while(ptr>Pelectron)*--ptr=*--pele;
01115 
01116   ptr = Ppositron+Pdim*Npositron;
01117   while(ptr>Ppositron)*--ptr=*--ppos;
01118 
01119   ptr = Pneutron+Pdim*Nneutron;
01120   while(ptr>Pneutron)*--ptr=*--pneu;
01121 
01122   ptr = Pgamma+Pdim*Ngamma;
01123   while(ptr>Pgamma)*--ptr=*--pgam;
01124 
01125   ptr = Rnubar+Rdim*Nnubar;
01126   while(ptr>Rnubar)*--ptr=*--rnu;
01127 
01128   ptr = Relectron+Rdim*Nelectron;
01129   while(ptr>Relectron)*--ptr=*--rele;
01130 
01131   ptr = Rgenele+Rdim*Nelectron;
01132   while(ptr>Rgenele)*--ptr=*--rge;
01133 
01134   ptr = Rpositron+Rdim*Npositron;
01135   while(ptr>Rpositron)*--ptr=*--rpos;
01136 
01137   ptr = Rgenpos+Rdim*Npositron;
01138   while(ptr>Rgenpos)*--ptr=*--rgp;
01139 
01140   ptr = Rneutron+Rdim*Nneutron;
01141   while(ptr>Rneutron)*--ptr=*--rneu;
01142 
01143   ptr = Rgenneu+Rdim*Nneutron;
01144   while(ptr>Rgenneu)*--ptr=*--rgn;
01145 
01146   ptr = Rgamma+Rdim*Ngamma;
01147   while(ptr>Rgamma)*--ptr=*--rgam;
01148 
01149   /*
01150   cout << Nnubar << " nubar(s):" << endl;
01151   for(int i=0;i<Nnubar*Pdim;i++){
01152     cout << " Pnubar[" << i << "] " << Pnubar[i] << endl;
01153   }
01154   for(int i=0;i<Nnubar*Rdim;i++){
01155     cout << "  Rnubar[" << i << "] " << Rnubar[i] << endl;
01156   }
01157 
01158   cout << Nelectron << " electron(s):" << endl;
01159   for(int i=0;i<Nelectron*Pdim;i++){
01160     cout << " Pelectron[" << i << "] " << Pelectron[i] << endl;
01161   }
01162   for(int i=0;i<Nelectron*Rdim;i++){
01163     cout << "  Relectron[" << i << "] " << Relectron[i] << endl;
01164   }
01165   for(int i=0;i<Nelectron*Rdim;i++){
01166     cout << "  Rgenele[" << i << "] " << Rgenele[i] << endl;
01167   }
01168 
01169   cout << Npositron << " positron(s):" << endl;
01170   for(int i=0;i<Npositron*Pdim;i++){
01171     cout << " Ppositron[" << i << "] " << Ppositron[i] << endl;
01172   }
01173   for(int i=0;i<Npositron*Rdim;i++){
01174     cout << "  Rpositron[" << i << "] " << Rpositron[i] << endl;
01175   }
01176   for(int i=0;i<Npositron*Rdim;i++){
01177     cout << "  Rgenpos[" << i << "] " << Rgenpos[i] << endl;
01178   }
01179 
01180   cout << Nneutron << " neutron(s):" << endl;
01181   for(int i=0;i<Nneutron*Pdim;i++){
01182     cout << " Pneutron[" << i << "] " << Pneutron[i] << endl;
01183   }
01184   for(int i=0;i<Nneutron*Rdim;i++){
01185     cout << "  Rneutron[" << i << "] " << Rneutron[i] << endl;
01186   }
01187   for(int i=0;i<Nneutron*Rdim;i++){
01188     cout << "  Rgenneu[" << i << "] " << Rgenneu[i] << endl;
01189   }
01190 
01191   cout << Ngamma << " photon(s):" << endl;
01192   for(int i=0;i<Ngamma*Pdim;i++){
01193     cout << " Pgamma[" << i << "] " << Pgamma[i] << endl;
01194   }
01195   for(int i=0;i<Ngamma*Rdim;i++){
01196     cout << "  Rgamma[" << i << "] " << Rgamma[i] << endl;
01197   }
01198   */
01199 
01200   // clean up
01201   delete[] pnubar;
01202   delete[] pelectron;
01203   delete[] ppositron;
01204   delete[] pneutron;
01205   delete[] pgamma;
01206   delete[] rnubar;
01207   delete[] relectron;
01208   delete[] rgenele;
01209   delete[] rpositron;
01210   delete[] rgenpos;
01211   delete[] rneutron;
01212   delete[] rgenneu;
01213   delete[] rgamma;
01214 }
01215 
01216 ReactorEvent::~ReactorEvent(){
01217   delete[] Pnubar;
01218   delete[] Rnubar;
01219   delete[] Pelectron;
01220   delete[] Relectron;
01221   delete[] Rgenele;
01222   delete[] Ppositron;
01223   delete[] Rpositron;
01224   delete[] Rgenpos;
01225   delete[] Pneutron;
01226   delete[] Rneutron;
01227   delete[] Rgenneu;
01228   delete[] Pgamma;
01229   delete[] Rgamma;
01230 }
01231 
01232 ReactorEvent::ReactorEvent(const ReactorEvent& Event){
01233 
01234   Pdim = 7; 
01235   Rdim = 7;
01236 
01237   XCorder_RE = Event.XCorder_RE;  
01238   Emin_RE = Event.Emin_RE;
01239   Emax_RE = Event.Emax_RE;
01240   Rmaxgen_RE = Event.Rmaxgen_RE;
01241   Nevents = Event.Nevents;
01242 
01243   XCWeight = Event.XCWeight;
01244   FluxWeight = Event.FluxWeight;
01245 
01246   Nnubar    = Event.Nnubar;
01247   Nelectron = Event.Nelectron;
01248   Npositron = Event.Npositron;
01249   Nneutron  = Event.Nneutron;
01250   Ngamma    = Event.Ngamma;
01251 
01252   double* pold;
01253   double* pnew;
01254 
01255   Pnubar = new double[Nnubar];
01256   pold = Event.Pnubar+Pdim*Nnubar;
01257   pnew = Pnubar+Pdim*Nnubar;
01258   while(pnew>Pnubar)*--pnew=*--pold;
01259 
01260   Pelectron = new double[Nelectron];
01261   pold = Event.Pelectron+Pdim*Nelectron;
01262   pnew = Pelectron+Pdim*Nelectron;
01263   while(pnew>Pelectron)*--pnew=*--pold;
01264 
01265   Ppositron = new double[Npositron];
01266   pold = Event.Ppositron+Pdim*Npositron;
01267   pnew = Ppositron+Pdim*Npositron;
01268   while(pnew>Ppositron)*--pnew=*--pold;
01269 
01270   Pneutron = new double[Nneutron];
01271   pold = Event.Pneutron+Pdim*Nneutron;
01272   pnew = Pneutron+Pdim*Nneutron;
01273   while(pnew>Pneutron)*--pnew=*--pold;
01274 
01275   Pgamma = new double[Ngamma];
01276   pold = Event.Pgamma+Pdim*Ngamma;
01277   pnew = Pgamma+Pdim*Ngamma;
01278   while(pnew>Pgamma)*--pnew=*--pold;
01279 
01280   double* rold;
01281   double* rnew;
01282 
01283   Rnubar = new double[Nnubar];
01284   rold = Event.Rnubar+Rdim*Nnubar;
01285   rnew = Rnubar+Rdim*Nnubar;
01286   while(rnew>Rnubar)*--rnew=*--rold;
01287 
01288   Relectron = new double[Nelectron];
01289   rold = Event.Relectron+Rdim*Nelectron;
01290   rnew = Relectron+Rdim*Nelectron;
01291   while(rnew>Relectron)*--rnew=*--rold;
01292 
01293   Rgenele = new double[Nelectron];
01294   rold = Event.Rgenele+Rdim*Nelectron;
01295   rnew = Rgenele+Rdim*Nelectron;
01296   while(rnew>Rgenele)*--rnew=*--rold;
01297 
01298   Rpositron = new double[Npositron];
01299   rold = Event.Rpositron+Rdim*Npositron;
01300   rnew = Rpositron+Rdim*Npositron;
01301   while(rnew>Rpositron)*--rnew=*--rold;
01302 
01303   Rgenpos = new double[Npositron];
01304   rold = Event.Rgenpos+Rdim*Npositron;
01305   rnew = Rgenpos+Rdim*Npositron;
01306   while(rnew>Rgenpos)*--rnew=*--rold;
01307 
01308   Rneutron = new double[Nneutron];
01309   rold = Event.Rneutron+Rdim*Nneutron;
01310   rnew = Rneutron+Rdim*Nneutron;
01311   while(rnew>Rneutron)*--rnew=*--rold;
01312 
01313   Rgenneu = new double[Nneutron];
01314   rold = Event.Rgenneu+Rdim*Nneutron;
01315   rnew = Rgenneu+Rdim*Nneutron;
01316   while(rnew>Rgenneu)*--rnew=*--rold;
01317 
01318   Rgamma = new double[Ngamma];
01319   rold = Event.Rgamma+Rdim*Ngamma;
01320   rnew = Rgamma+Rdim*Ngamma;
01321   while(rnew>Rgamma)*--rnew=*--rold;
01322 }    
01323 
01324 ReactorEvent& ReactorEvent::operator=(const ReactorEvent& rhs){
01325 
01326   Pdim = 7; 
01327   Rdim = 7;
01328   if(this == &rhs)return *this;
01329 
01330   XCorder_RE = rhs.XCorder_RE;  
01331   Emin_RE = rhs.Emin_RE;
01332   Emax_RE = rhs.Emax_RE;
01333   Rmaxgen_RE = rhs.Rmaxgen_RE;
01334   Nevents = rhs.Nevents;
01335 
01336   XCWeight = rhs.XCWeight;
01337   FluxWeight = rhs.FluxWeight;
01338   
01339   Nnubar    = rhs.Nnubar;
01340   Nelectron = rhs.Nelectron;
01341   Npositron = rhs.Npositron;
01342   Nneutron  = rhs.Nneutron;
01343   Ngamma    = rhs.Ngamma;
01344 
01345   double* pold;
01346   double* pnew;
01347 
01348   delete[] Pnubar;
01349   Pnubar = new double[Nnubar];
01350   pold = rhs.Pnubar+Pdim*Nnubar;
01351   pnew = Pnubar+Pdim*Nnubar;
01352   while(pnew>Pnubar)*--pnew=*--pold;
01353 
01354   delete[] Pelectron;
01355   Pelectron = new double[Nelectron];
01356   pold = rhs.Pelectron+Pdim*Nelectron;
01357   pnew = Pelectron+Pdim*Nelectron;
01358   while(pnew>Pelectron)*--pnew=*--pold;
01359 
01360   delete[] Ppositron;
01361   Ppositron = new double[Npositron];
01362   pold = rhs.Ppositron+Pdim*Npositron;
01363   pnew = Ppositron+Pdim*Npositron;
01364   while(pnew>Ppositron)*--pnew=*--pold;
01365 
01366   delete[] Pneutron;
01367   Pneutron = new double[Nneutron];
01368   pold = rhs.Pneutron+Pdim*Nneutron;
01369   pnew = Pneutron+Pdim*Nneutron;
01370   while(pnew>Pneutron)*--pnew=*--pold;
01371 
01372   delete[] Pgamma;
01373   Pgamma = new double[Ngamma];
01374   pold = rhs.Pgamma+Pdim*Ngamma;
01375   pnew = Pgamma+Pdim*Ngamma;
01376   while(pnew>Pgamma)*--pnew=*--pold;
01377 
01378   double* rold;
01379   double* rnew;
01380 
01381   delete[] Rnubar;
01382   Rnubar = new double[Nnubar];
01383   rold = rhs.Rnubar+Rdim*Nnubar;
01384   rnew = Rnubar+Rdim*Nnubar;
01385   while(rnew>Rnubar)*--rnew=*--rold;
01386 
01387   delete[] Relectron;
01388   Relectron = new double[Nelectron];
01389   rold = rhs.Relectron+Rdim*Nelectron;
01390   rnew = Relectron+Rdim*Nelectron;
01391   while(rnew>Relectron)*--rnew=*--rold;
01392 
01393   delete[] Rgenele;
01394   Rgenele = new double[Nelectron];
01395   rold = rhs.Rgenele+Rdim*Nelectron;
01396   rnew = Rgenele+Rdim*Nelectron;
01397   while(rnew>Rgenele)*--rnew=*--rold;
01398 
01399   delete[] Rpositron;
01400   Rpositron = new double[Npositron];
01401   rold = rhs.Rpositron+Rdim*Npositron;
01402   rnew = Rpositron+Rdim*Npositron;
01403   while(rnew>Rpositron)*--rnew=*--rold;
01404 
01405   delete[] Rgenpos;
01406   Rgenpos = new double[Npositron];
01407   rold = rhs.Rgenpos+Rdim*Npositron;
01408   rnew = Rgenpos+Rdim*Npositron;
01409   while(rnew>Rgenpos)*--rnew=*--rold;
01410 
01411   delete[] Rneutron;
01412   Rneutron = new double[Nneutron];
01413   rold = rhs.Rneutron+Rdim*Nneutron;
01414   rnew = Rneutron+Rdim*Nneutron;
01415   while(rnew>Rneutron)*--rnew=*--rold;
01416 
01417   delete[] Rgenneu;
01418   Rgenneu = new double[Nneutron];
01419   rold = rhs.Rgenneu+Rdim*Nneutron;
01420   rnew = Rgenneu+Rdim*Nneutron;
01421   while(rnew>Rgenneu)*--rnew=*--rold;
01422 
01423   delete[] Rgamma;
01424   Rgamma = new double[Ngamma];
01425   rold = rhs.Rgamma+Rdim*Ngamma;
01426   rnew = Rgamma+Rdim*Ngamma;
01427   while(rnew>Rgamma)*--rnew=*--rold;
01428 
01429   return *this;
01430 }