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 "ReactorConstants.hh"
00013 #include "ReactorFortran.hh"
00014 #include "ReactorXC.hh"
00015 #include "ReactorEvent.hh"
00016 #include "MyCfSource.hh"
00017 #include "MuonPropagator.hh"
00018 
00019 
00020 ReactorEvent::ReactorEvent(double newRmaxgen,
00021                            double newEmin,
00022                            double newEmax,
00023                            int newXCorder){
00024 
00025   Rmaxgen = newRmaxgen;
00026   Emin = newEmin;
00027   Emax = newEmax;
00028   XCorder = newXCorder;
00029 
00030   int nmax = k.Ndim;
00031   int nnu  = 0;
00032   int npos = 0;
00033   int nele = 0;
00034   int nneu = 0;
00035   int ngam = 0;
00036 
00037   double* pnubar = new double[nmax*Pdim];
00038   double* pelectron = new double[nmax*Pdim];
00039   double* ppositron = new double[nmax*Pdim];
00040   double* pneutron = new double[nmax*Pdim];
00041   double* pgamma = new double[nmax*Pdim];
00042 
00043   double* pnu  = pnubar;
00044   double* pele = pelectron;
00045   double* ppos = ppositron;
00046   double* pneu = pneutron;
00047   double* pgam = pgamma;
00048 
00049   double* rnubar = new double[nmax*Rdim];
00050   double* relectron = new double[nmax*Rdim];
00051   double* rgenele = new double[nmax*Rdim];
00052   double* rpositron = new double[nmax*Rdim];
00053   double* rgenpos = new double[nmax*Rdim];
00054   double* rneutron = new double[nmax*Rdim];
00055   double* rgenneu = new double[nmax*Rdim];
00056   double* rgamma = new double[nmax*Rdim];
00057 
00058   double* rnu  = rnubar;
00059   double* rele = relectron;
00060   double* rge  = rgenele;
00061   double* rpos = rpositron;
00062   double* rgp  = rgenpos;
00063   double* rneu = rneutron;
00064   double* rgn  = rgenneu;
00065   double* rgam = rgamma;
00066 
00067   double Enu,z,phi;
00068   double R,zR,phiR;
00069 
00070   // setup the event
00071   fstream fin("src/event_setup.txt", ios::in);
00072   //cout << "setup file: src/event_setup.txt" << endl;
00073 
00074   char paramLine[256];
00075   std::string dummy;
00076   bool done = false;
00077   bool set = false;
00078   std::string static eventtype;
00079   std::string static pmtrad;
00080   double static Depth;
00081   bool static first = true;
00082 
00083   if(first){
00084     while(!done){
00085 
00086       // '#' specifies a comment line -- ignore it
00087       if(fin.peek() == '#'){
00088 
00089         fin.getline(paramLine, 256);
00090         continue;
00091       }
00092 
00093       // 'E'ND or 'e'nd specifies end of file -- set done flag to true
00094       if(fin.peek() == 'e' || fin.peek() == 'E'){
00095 
00096         done = true;
00097         continue;
00098       }
00099 
00100       // 'T'YPE or 't'ype specifies the event type
00101       if(fin.peek() == 't' || fin.peek() == 'T'){
00102 
00103         if(set){
00104 
00105           fin.getline( paramLine, 256 );
00106           continue;
00107         }
00108         fin >> dummy;
00109         while(fin.peek() == ' ' || fin.peek() == '=') 
00110           fin.ignore(1);
00111 
00112         fin >> eventtype;
00113 
00114         if(eventtype == "reactor" || eventtype == "cf252" || 
00115            eventtype == "pmtrad"  || eventtype == "muon" || eventtype == "all"){
00116           set = true;
00117         }
00118         else{
00119           cout << "  SETUP ERROR: unrecognized event type!" << endl;
00120         }
00121 
00122         continue;
00123       }
00124 
00125       // 'P'MTRAD or 'p'mtrad turns on/off radioactive decays in the PMTs
00126       if(fin.peek() == 'p' || fin.peek() == 'P'){
00127 
00128         fin >> dummy;
00129         while(fin.peek() == ' ' || fin.peek() == '=') 
00130           fin.ignore(1);
00131 
00132         fin >> pmtrad;
00133         continue;
00134       }
00135 
00136       // 'D'EPTH or 'd'epth sets the detector depth in mwe
00137       if(fin.peek() == 'd' || fin.peek() == 'D'){
00138 
00139         fin >> dummy;
00140         while(fin.peek() == ' ' || fin.peek() == '=') 
00141           fin.ignore(1);
00142 
00143         fin >> Depth;
00144         continue;
00145       }
00146 
00147       // skip n, o, r, and f
00148       if(fin.peek() == 'N' || fin.peek() == 'n' ||
00149          fin.peek() == 'O' || fin.peek() == 'o' ||
00150          fin.peek() == 'R' || fin.peek() == 'r' ||
00151          fin.peek() == 'F' || fin.peek() == 'f'){
00152 
00153         fin.getline(paramLine, 256);
00154         continue;
00155       }
00156 
00157       // ignore extra whitespace
00158       if(fin.peek() == ' ' || fin.peek() == '\n'){
00159 
00160         fin.ignore(1);
00161         continue;
00162       }
00163     }
00164     cout << "  " << eventtype << " type of events" << endl;
00165     cout << "  PMT radioactivity is " << pmtrad << endl;
00166     cout << "===============================" << endl;
00167 
00168     // zero the event counter
00169     Nevent = 0;
00170 
00171     first = false;
00172   }
00173   eventType = eventtype;
00174   pmtRad = pmtrad;
00175   //
00176   // count the event
00177   //
00178   Nevent++;
00179   if(Nevent%100 == 0) cout << Nevent << " events processed" << endl;
00180   //
00181   // if all events are requested, mix the 'reactor' and 'cf252' 
00182   // events together with 90/10% Cf252/reactor events
00183   if(eventType == "all"){
00184     int len = 1;
00185     float r[len];
00186     ranlux_(r,len);
00187 
00188     if(r[0] < 0.90){
00189       eventType = "cf252";
00190     }
00191     else{
00192       eventType = "reactor";
00193     }
00194   }
00195   //cout << eventType << " type of events" << endl;
00196   //cout << "PMT radiation " << pmtRad << endl;
00197 
00198   if(eventType == "reactor"){
00199 
00200     XCWeight = 1;
00201 
00202     // inverse beta decay reaction
00203     double XCmax = InverseBeta(Emax,-1);
00204     double FluxMax = RMPflux(Emin);
00205 
00206     bool passed=0;
00207 
00208     while(!passed){
00209       int len = 7;
00210       float r[len];
00211       ranlux_(r,len);
00212       Enu = Emin+(Emax-Emin)*r[0];
00213       z = -1+2*r[1];
00214       phi = 2*k.pi*r[2];
00215 
00216       R = Rmaxgen*exp(log(r[3])/3);
00217       zR = -1+2*r[4];
00218       phiR = 2*k.pi*r[5];
00219 
00220       float XCtest = XCmax*FluxMax*r[6];
00221       XCWeight = InverseBeta(Enu,z);
00222       FluxWeight=RMPflux(Enu);
00223       passed= XCWeight*FluxWeight>XCtest;
00224     }
00225   
00226     double E0 = Enu - k.Delta;
00227     double p0 = sqrt(E0*E0-k.Melectron*k.Melectron);
00228     double v0 = p0/E0;
00229     double Ysquared = (k.Delta*k.Delta-k.Melectron*k.Melectron)/2;
00230     double E1 = E0*(1-Enu/k.Mproton*(1-v0*z))-Ysquared/k.Mproton;
00231     double p1 = sqrt(E1*E1-k.Melectron*k.Melectron);
00232 
00233     *pnu = Enu; pnu++;
00234     *pnu = 0; pnu++;
00235     *pnu = 0; pnu++;
00236     *pnu = Enu; pnu++;
00237     *pnu = Enu; pnu++;
00238     *pnu = Enu; pnu++;
00239     *pnu = 1.; pnu++; // inverse beta decay process Id = 1
00240 
00241     *ppos = E1; ppos++;
00242     *ppos = p1*sqrt(1-z*z)*cos(phi); ppos++;
00243     *ppos = p1*sqrt(1-z*z)*sin(phi); ppos++;
00244     *ppos = p1*z; ppos++;
00245     *ppos = E1-k.Melectron; ppos++;
00246     *ppos = p1; ppos++;
00247     *ppos = 1.; ppos++;
00248 
00249     *pneu = Enu+k.Mproton-E1; pneu++;
00250     *pneu = -p1*sqrt(1-z*z)*cos(phi); pneu++;
00251     *pneu = -p1*sqrt(1-z*z)*sin(phi); pneu++;
00252     *pneu = Enu-p1*z; pneu++;
00253     *pneu = Enu-E1-k.Delta; pneu++;
00254     *pneu = sqrt((Enu+k.Mproton-E1)*(Enu+k.Mproton-E1)-
00255                          (k.Mproton+k.Delta)*(k.Mproton+k.Delta)); pneu++;
00256     *pneu = 1.; pneu++;
00257 
00258     *rnu = 0; rnu++;
00259     *rnu = R*sqrt(1-zR*zR)*cos(phiR); rnu++;
00260     *rnu = R*sqrt(1-zR*zR)*sin(phiR); rnu++;
00261     *rnu = R*zR; rnu++;
00262     *rnu = R; rnu++;
00263     *rnu = zR; rnu++;
00264     *rnu = phiR; rnu++;
00265 
00266     *rpos = 0; rpos++;
00267     *rpos = R*sqrt(1-zR*zR)*cos(phiR); rpos++;
00268     *rpos = R*sqrt(1-zR*zR)*sin(phiR); rpos++;
00269     *rpos = R*zR; rpos++;
00270     *rpos = R; rpos++;
00271     *rpos = zR; rpos++;
00272     *rpos = phiR; rpos++;
00273 
00274     *rgp = 0; rgp++;
00275     *rgp = R*sqrt(1-zR*zR)*cos(phiR); rgp++;
00276     *rgp = R*sqrt(1-zR*zR)*sin(phiR); rgp++;
00277     *rgp = R*zR; rgp++;
00278     *rgp = R; rgp++;
00279     *rgp = zR; rgp++;
00280     *rgp = phiR; rgp++;
00281 
00282     *rneu = 0; rneu++;
00283     *rneu = R*sqrt(1-zR*zR)*cos(phiR); rneu++;
00284     *rneu = R*sqrt(1-zR*zR)*sin(phiR); rneu++;
00285     *rneu = R*zR; rneu++;
00286     *rneu = R; rneu++;
00287     *rneu = zR; rneu++; 
00288     *rneu = phiR; rneu++;
00289 
00290     *rgn = 0; rgn++;
00291     *rgn = R*sqrt(1-zR*zR)*cos(phiR); rgn++;
00292     *rgn = R*sqrt(1-zR*zR)*sin(phiR); rgn++;
00293     *rgn = R*zR; rgn++;
00294     *rgn = R; rgn++;
00295     *rgn = zR; rgn++;
00296     *rgn = phiR; rgn++;
00297 
00298     // count the nubar, positron, and neutron
00299     nnu++;
00300     npos++;
00301     nneu++;
00302 
00303   } // done with reactor neutrino inverse-beta reaction
00304 
00305   // Cf source reaction
00306   if(eventType == "cf252"){
00307 
00308     XCWeight = 1;
00309 
00310     int Isotope = 252;
00311     MyCfSource CfSource(Isotope,Rmaxgen);
00312 
00313     int Nsources = CfSource.GetNumCfSources();
00314     //cout << "Number of Cf252 sources = " << Nsources << endl;
00315 
00316     for(int n=0;n<Nsources;n++){
00317 
00318       R = CfSource.GetCfVertexR(n);
00319       zR = CfSource.GetCfVertexCosTheta(n);
00320       phiR = CfSource.GetCfVertexPhi(n);
00321 
00322       // loop over the neutrons
00323       for(int nn=0;nn<CfSource.GetNumNeutron(n);nn++){
00324 
00325         double neutronKE = CfSource.GetCfNeutronEnergy(n,nn);
00326         //cout << "   Cf252 source " << n << ", neutron " << nn 
00327         //     << " energy = " << neutronKE << endl;
00328 
00329         double Tneu = CfSource.GetCfNeutronTime(n,nn);
00330         double Mneutron = k.Mproton+k.Delta;
00331         double neutronE = neutronKE+Mneutron;
00332 
00333         // pick a direction for neutron momentum
00334         int len = 2;
00335         float r[len];
00336         ranlux_(r,len);
00337         z = -1+2*r[0];
00338         phi = 2*k.pi*r[1];
00339         double x = sqrt(1-z*z)*cos(phi);
00340         double y = sqrt(1-z*z)*sin(phi);
00341 
00342         *pneu = neutronE; pneu++;
00343         *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*x; pneu++;
00344         *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*y; pneu++;
00345         *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*z; pneu++;
00346         *pneu = neutronKE; pneu++;
00347         *pneu = sqrt(neutronE*neutronE-Mneutron*Mneutron); pneu++;
00348         *pneu = 2.; pneu++; // cf252 source process Id = 2
00349 
00350         *rneu = Tneu; rneu++;
00351         *rneu = R*sqrt(1-zR*zR)*cos(phiR); rneu++;
00352         *rneu = R*sqrt(1-zR*zR)*sin(phiR); rneu++;
00353         *rneu = R*zR; rneu++;
00354         *rneu = R; rneu++;
00355         *rneu = zR; rneu++;
00356         *rneu = phiR; rneu++;
00357 
00358         *rgn = Tneu; rgn++;
00359         *rgn = R*sqrt(1-zR*zR)*cos(phiR); rgn++;
00360         *rgn = R*sqrt(1-zR*zR)*sin(phiR); rgn++;
00361         *rgn = R*zR; rgn++;
00362         *rgn = R; rgn++;
00363         *rgn = zR; rgn++;
00364         *rgn = phiR; rgn++;
00365       }
00366       // count the neutrons
00367       nneu += CfSource.GetNumNeutron(n);      
00368 
00369       // loop over the prompt gammas
00370       for(int nn=0;nn<CfSource.GetNumGamma(n);nn++){
00371 
00372         double gammaKE = CfSource.GetCfGammaEnergy(n,nn);
00373         //cout << "   Cf252 source " << n << ", gamma " << nn 
00374         //     << " energy = " << gammaKE << endl;
00375 
00376         double Tgam = CfSource.GetCfGammaTime(n,nn);
00377 
00378         // pick a direction for gamma momentum
00379         int len = 2;
00380         float r[len];
00381         ranlux_(r,len);
00382         z = -1+2*r[0];
00383         phi = 2*k.pi*r[1];
00384         double x = sqrt(1-z*z)*cos(phi);
00385         double y = sqrt(1-z*z)*sin(phi);
00386 
00387         *pgam = gammaKE; pgam++;
00388         *pgam = gammaKE*x; pgam++;
00389         *pgam = gammaKE*y; pgam++;
00390         *pgam = gammaKE*z; pgam++;
00391         *pgam = gammaKE; pgam++;
00392         *pgam = gammaKE; pgam++;
00393         *pgam = 2.; pgam++; // cf252 source process Id = 2
00394 
00395         *rgam = Tgam; rgam++;
00396         *rgam = R*sqrt(1-zR*zR)*cos(phiR); rgam++;
00397         *rgam = R*sqrt(1-zR*zR)*sin(phiR); rgam++;
00398         *rgam = R*zR; rgam++;
00399         *rgam = R; rgam++;
00400         *rgam = zR; rgam++;
00401         *rgam = phiR; rgam++;
00402       }
00403       //
00404       ngam += CfSource.GetNumGamma(n);
00405 
00406     } // done looping over sources
00407 
00408   } // done with cf source neutron production
00409 
00410   // PMT radiation
00411   if(eventType == "pmtrad"){
00412 
00413     XCWeight = 1;
00414 
00415     // nothing happens in these events but photons, generated in 
00416     // ReactorDetector::GenerateGammas with pmtRad on
00417     if(pmtRad == "off"){
00418       cout << "  SETUP ERROR: PMT radiation set " << pmtRad << " in a " 
00419            << eventType << " type event" << endl;
00420       pmtRad = "on";
00421       cout << "  Forcing PMT radiation " << pmtRad << endl;
00422     }
00423   }
00424 
00425   // muons
00426   if(eventType == "muon"){
00427 
00428     XCWeight = 1;
00429 
00430     //cout << eventType << endl;
00431     MuonPropagator Muon(Rmaxgen,Depth);
00432     //
00433     // loop over neutrons
00434     //
00435     nneu = Muon.GetNumNeutron();
00436     for(int n=0;n<nneu;n++){
00437 
00438       float neutronKE = Muon.GetNeutronKE(n);
00439       //cout << "  neutron " << n << " energy = " << neutronKE << endl;
00440 
00441       float Mneutron = k.Mproton+k.Delta;
00442       float neutronE = neutronKE+Mneutron;
00443 
00444       // pick a direction for neutron momentum
00445       int len = 2;
00446       float rv[len];
00447       ranlux_(rv,len);
00448       z = -1+2*rv[0];
00449       phi = 2*k.pi*rv[1];
00450       double x = sqrt(1-z*z)*cos(phi);
00451       double y = sqrt(1-z*z)*sin(phi);
00452 
00453       *pneu = neutronE; pneu++;
00454       *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*x; pneu++;
00455       *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*y; pneu++;
00456       *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*z; pneu++;
00457       *pneu = neutronKE; pneu++;
00458       *pneu = sqrt(neutronE*neutronE-Mneutron*Mneutron); pneu++;
00459       *pneu = 3.; pneu++; // muon process Id = 3
00460 
00461       double neuxyz[3];
00462       for(int i=0;i<3;i++){
00463         neuxyz[i] = *(Muon.GetNeutronVertexXYZ()+(i+n));
00464         //cout << "  neuxyz[" << i << "] " << neuxyz[i] << endl;
00465       }
00466       double rr = sqrt(neuxyz[0]*neuxyz[0] + 
00467                        neuxyz[1]*neuxyz[1] + 
00468                        neuxyz[2]*neuxyz[2]);
00469       double rphi = atan2(neuxyz[1],neuxyz[0]);
00470       if(rphi<0) rphi += 2*k.pi;
00471 
00472       double time = Muon.GetNeutronTime(n);
00473       //cout << " n time " << time << endl;
00474 
00475       *rneu = time; rneu++;
00476       *rneu = neuxyz[0]; rneu++;
00477       *rneu = neuxyz[1]; rneu++;
00478       *rneu = neuxyz[2]; rneu++;
00479       *rneu = rr; rneu++;
00480       *rneu = neuxyz[2]/rr; rneu++;
00481       *rneu = rphi; rneu++;
00482 
00483       *rgn = time; rgn++;
00484       *rgn = neuxyz[0]; rgn++;
00485       *rgn = neuxyz[1]; rgn++;
00486       *rgn = neuxyz[2]; rgn++;
00487       *rgn = rr; rgn++;
00488       *rgn = neuxyz[2]/rr; rgn++;
00489       *rgn = rphi; rgn++;
00490     }
00491     //
00492     // loop over electrons
00493     //
00494     nele = Muon.GetNumElectron();
00495     for(int n=0;n<nele;n++){
00496 
00497       float electronKE = Muon.GetElectronKE(n);
00498       //cout << "  electron " << n << " KE = " << electronKE << endl;
00499 
00500       float Melectron = k.Melectron;
00501       float electronE = electronKE+Melectron;
00502       //cout << "  electron " << n << " energy = " << electronE << endl;
00503 
00504       // pick a direction for electron momentum
00505       int len = 2;
00506       float r[len];
00507       ranlux_(r,len);
00508       z = -1+2*r[0];
00509       phi = 2*k.pi*r[1];
00510       double x = sqrt(1-z*z)*cos(phi);
00511       double y = sqrt(1-z*z)*sin(phi);
00512 
00513       *pele = electronE; pele++;
00514       *pele = sqrt((electronE*electronE)-(Melectron*Melectron))*x; pele++;
00515       *pele = sqrt((electronE*electronE)-(Melectron*Melectron))*y; pele++;
00516       *pele = sqrt((electronE*electronE)-(Melectron*Melectron))*z; pele++;
00517       *pele = electronKE; pele++;
00518       *pele = sqrt(electronE*electronE-Melectron*Melectron); pele++;
00519       *pele = 3.; pele++; // muon process Id = 3
00520 
00521       double elexyz[3];
00522       for(int i=0;i<3;i++){
00523         elexyz[i] = *(Muon.GetElectronVertexXYZ()+(i+3*n));
00524         //cout << "  elexyz[" << i << "] " << elexyz[i] << endl;
00525       }
00526       double rr = sqrt(elexyz[0]*elexyz[0] + 
00527                        elexyz[1]*elexyz[1] + 
00528                        elexyz[2]*elexyz[2]);
00529       double rphi = atan2(elexyz[1],elexyz[0]);
00530       if(rphi<0) rphi += 2*k.pi;
00531 
00532       double time = Muon.GetElectronTime(n);
00533       //cout << " e time " << time << endl;
00534 
00535       *rele = time; rele++;
00536       *rele = elexyz[0]; rele++;
00537       *rele = elexyz[1]; rele++;
00538       *rele = elexyz[2]; rele++;
00539       *rele = rr; rele++;
00540       *rele = elexyz[2]/rr; rele++;
00541       *rele = rphi; rele++;
00542 
00543       *rge = time; rge++;
00544       *rge = elexyz[0]; rge++;
00545       *rge = elexyz[1]; rge++;
00546       *rge = elexyz[2]; rge++;
00547       *rge = rr; rge++;
00548       *rge = elexyz[2]/rr; rge++;
00549       *rge = rphi; rge++;
00550     }
00551     //
00552     // loop over photons (these contain the energy loss of the muon)
00553     //
00554     ngam = Muon.GetNumGamma();
00555     for(int n=0;n<ngam;n++){
00556 
00557       float gammaKE = Muon.GetGammaKE(n);
00558       //cout << "  gamma " << n << " energy = " << gammaKE << endl;
00559 
00560       // pick a direction for gamma momentum
00561       int len = 2;
00562       float r[len];
00563       ranlux_(r,len);
00564       z = -1+2*r[0];
00565       phi = 2*k.pi*r[1];
00566       double x = sqrt(1-z*z)*cos(phi);
00567       double y = sqrt(1-z*z)*sin(phi);
00568 
00569       *pgam = gammaKE; pgam++;
00570       *pgam = gammaKE*x; pgam++;
00571       *pgam = gammaKE*y; pgam++;
00572       *pgam = gammaKE*z; pgam++;
00573       *pgam = gammaKE; pgam++;
00574       *pgam = gammaKE; pgam++;
00575       *pgam = 3.; pgam++; // muon process Id = 3
00576 
00577       double gamxyz[3];
00578       for(int i=0;i<3;i++){
00579         gamxyz[i] = *(Muon.GetGammaVertexXYZ()+(i+3*n));
00580         //cout << "  gamxyz[" << i << "] " << gamxyz[i] << endl;
00581       }
00582       double rr = sqrt(gamxyz[0]*gamxyz[0] + 
00583                        gamxyz[1]*gamxyz[1] + 
00584                        gamxyz[2]*gamxyz[2]);
00585       double rphi = atan2(gamxyz[1],gamxyz[0]);
00586       if(rphi<0) rphi += 2*k.pi;
00587 
00588       double time = Muon.GetGammaTime(n);
00589 
00590       *rgam = time; rgam++;
00591       *rgam = gamxyz[0]; rgam++;
00592       *rgam = gamxyz[1]; rgam++;
00593       *rgam = gamxyz[2]; rgam++;
00594       *rgam = rr; rgam++;
00595       *rgam = gamxyz[2]/rr; rgam++;
00596       *rgam = rphi; rgam++;
00597     }
00598     XCWeight = Muon.GetWeight();
00599   }
00600   //
00601   // transfer temp data to storage
00602   //
00603   Nnubar    = nnu;
00604   Nelectron = nele;
00605   Npositron = npos;
00606   Nneutron  = nneu;
00607   Ngamma    = ngam;
00608 
00609   Pnubar    = new double[Pdim*nnu];
00610   Pelectron = new double[Pdim*nele];
00611   Ppositron = new double[Pdim*npos];
00612   Pneutron  = new double[Pdim*nneu];
00613   Pgamma    = new double[Pdim*ngam];
00614 
00615   Rnubar    = new double[Rdim*nnu];
00616   Relectron = new double[Rdim*nele];
00617   Rgenele   = new double[Rdim*nele];
00618   Rpositron = new double[Rdim*npos];
00619   Rgenpos   = new double[Rdim*npos];
00620   Rneutron  = new double[Rdim*nneu];
00621   Rgenneu   = new double[Rdim*nneu];
00622   Rgamma    = new double[Rdim*ngam];
00623 
00624   double* ptr;
00625 
00626   ptr = Pnubar+Pdim*Nnubar;
00627   while(ptr>Pnubar)*--ptr=*--pnu;
00628 
00629   ptr = Pelectron+Pdim*Nelectron;
00630   while(ptr>Pelectron)*--ptr=*--pele;
00631 
00632   ptr = Ppositron+Pdim*Npositron;
00633   while(ptr>Ppositron)*--ptr=*--ppos;
00634 
00635   ptr = Pneutron+Pdim*Nneutron;
00636   while(ptr>Pneutron)*--ptr=*--pneu;
00637 
00638   ptr = Pgamma+Pdim*Ngamma;
00639   while(ptr>Pgamma)*--ptr=*--pgam;
00640 
00641   ptr = Rnubar+Rdim*Nnubar;
00642   while(ptr>Rnubar)*--ptr=*--rnu;
00643 
00644   ptr = Relectron+Rdim*Nelectron;
00645   while(ptr>Relectron)*--ptr=*--rele;
00646 
00647   ptr = Rgenele+Rdim*Nelectron;
00648   while(ptr>Rgenele)*--ptr=*--rge;
00649 
00650   ptr = Rpositron+Rdim*Npositron;
00651   while(ptr>Rpositron)*--ptr=*--rpos;
00652 
00653   ptr = Rgenpos+Rdim*Npositron;
00654   while(ptr>Rgenpos)*--ptr=*--rgp;
00655 
00656   ptr = Rneutron+Rdim*Nneutron;
00657   while(ptr>Rneutron)*--ptr=*--rneu;
00658 
00659   ptr = Rgenneu+Rdim*Nneutron;
00660   while(ptr>Rgenneu)*--ptr=*--rgn;
00661 
00662   ptr = Rgamma+Rdim*Ngamma;
00663   while(ptr>Rgamma)*--ptr=*--rgam;
00664 
00665   /*
00666   cout << Nnubar << " nubar(s):" << endl;
00667   for(int i=0;i<Nnubar*Pdim;i++){
00668     cout << " Pnubar[" << i << "] " << Pnubar[i] << endl;
00669   }
00670   for(int i=0;i<Nnubar*Rdim;i++){
00671     cout << "  Rnubar[" << i << "] " << Rnubar[i] << endl;
00672   }
00673 
00674   cout << Nelectron << " electron(s):" << endl;
00675   for(int i=0;i<Nelectron*Pdim;i++){
00676     cout << " Pelectron[" << i << "] " << Pelectron[i] << endl;
00677   }
00678   for(int i=0;i<Nelectron*Rdim;i++){
00679     cout << "  Relectron[" << i << "] " << Relectron[i] << endl;
00680   }
00681   for(int i=0;i<Nelectron*Rdim;i++){
00682     cout << "  Rgenele[" << i << "] " << Rgenele[i] << endl;
00683   }
00684 
00685   cout << Npositron << " positron(s):" << endl;
00686   for(int i=0;i<Npositron*Pdim;i++){
00687     cout << " Ppositron[" << i << "] " << Ppositron[i] << endl;
00688   }
00689   for(int i=0;i<Npositron*Rdim;i++){
00690     cout << "  Rpositron[" << i << "] " << Rpositron[i] << endl;
00691   }
00692   for(int i=0;i<Npositron*Rdim;i++){
00693     cout << "  Rgenpos[" << i << "] " << Rgenpos[i] << endl;
00694   }
00695 
00696   cout << Nneutron << " neutron(s):" << endl;
00697   for(int i=0;i<Nneutron*Pdim;i++){
00698     cout << " Pneutron[" << i << "] " << Pneutron[i] << endl;
00699   }
00700   for(int i=0;i<Nneutron*Rdim;i++){
00701     cout << "  Rneutron[" << i << "] " << Rneutron[i] << endl;
00702   }
00703   for(int i=0;i<Nneutron*Rdim;i++){
00704     cout << "  Rgenneu[" << i << "] " << Rgenneu[i] << endl;
00705   }
00706 
00707   cout << Ngamma << " photon(s):" << endl;
00708   for(int i=0;i<Ngamma*Pdim;i++){
00709     cout << " Pgamma[" << i << "] " << Pgamma[i] << endl;
00710   }
00711   for(int i=0;i<Ngamma*Rdim;i++){
00712     cout << "  Rgamma[" << i << "] " << Rgamma[i] << endl;
00713   }
00714   */
00715 
00716   // clean up
00717   delete[] pnubar;
00718   delete[] pelectron;
00719   delete[] ppositron;
00720   delete[] pneutron;
00721   delete[] pgamma;
00722   delete[] rnubar;
00723   delete[] relectron;
00724   delete[] rgenele;
00725   delete[] rpositron;
00726   delete[] rgenpos;
00727   delete[] rneutron;
00728   delete[] rgenneu;
00729   delete[] rgamma;
00730 }
00731 
00732 ReactorEvent::~ReactorEvent(){
00733   delete[] Pnubar;
00734   delete[] Rnubar;
00735   delete[] Pelectron;
00736   delete[] Relectron;
00737   delete[] Rgenele;
00738   delete[] Ppositron;
00739   delete[] Rpositron;
00740   delete[] Rgenpos;
00741   delete[] Pneutron;
00742   delete[] Rneutron;
00743   delete[] Rgenneu;
00744   delete[] Pgamma;
00745   delete[] Rgamma;
00746 }
00747 
00748 ReactorEvent::ReactorEvent(const ReactorEvent& Event){
00749 
00750   XCorder = Event.XCorder;  
00751   Emin = Event.Emin;
00752   Emax = Event.Emax;
00753   Rmaxgen = Event.Rmaxgen;
00754 
00755   XCWeight = Event.XCWeight;
00756   FluxWeight = Event.FluxWeight;
00757 
00758   Nnubar    = Event.Nnubar;
00759   Nelectron = Event.Nelectron;
00760   Npositron = Event.Npositron;
00761   Nneutron  = Event.Nneutron;
00762   Ngamma    = Event.Ngamma;
00763 
00764   double* pold;
00765   double* pnew;
00766 
00767   Pnubar = new double[Nnubar];
00768   pold = Event.Pnubar+Pdim*Nnubar;
00769   pnew = Pnubar+Pdim*Nnubar;
00770   while(pnew>Pnubar)*--pnew=*--pold;
00771 
00772   Pelectron = new double[Nelectron];
00773   pold = Event.Pelectron+Pdim*Nelectron;
00774   pnew = Pelectron+Pdim*Nelectron;
00775   while(pnew>Pelectron)*--pnew=*--pold;
00776 
00777   Ppositron = new double[Npositron];
00778   pold = Event.Ppositron+Pdim*Npositron;
00779   pnew = Ppositron+Pdim*Npositron;
00780   while(pnew>Ppositron)*--pnew=*--pold;
00781 
00782   Pneutron = new double[Nneutron];
00783   pold = Event.Pneutron+Pdim*Nneutron;
00784   pnew = Pneutron+Pdim*Nneutron;
00785   while(pnew>Pneutron)*--pnew=*--pold;
00786 
00787   Pgamma = new double[Ngamma];
00788   pold = Event.Pgamma+Pdim*Ngamma;
00789   pnew = Pgamma+Pdim*Ngamma;
00790   while(pnew>Pgamma)*--pnew=*--pold;
00791 
00792   double* rold;
00793   double* rnew;
00794 
00795   Rnubar = new double[Nnubar];
00796   rold = Event.Rnubar+Rdim*Nnubar;
00797   rnew = Rnubar+Rdim*Nnubar;
00798   while(rnew>Rnubar)*--rnew=*--rold;
00799 
00800   Relectron = new double[Nelectron];
00801   rold = Event.Relectron+Rdim*Nelectron;
00802   rnew = Relectron+Rdim*Nelectron;
00803   while(rnew>Relectron)*--rnew=*--rold;
00804 
00805   Rgenele = new double[Nelectron];
00806   rold = Event.Rgenele+Rdim*Nelectron;
00807   rnew = Rgenele+Rdim*Nelectron;
00808   while(rnew>Rgenele)*--rnew=*--rold;
00809 
00810   Rpositron = new double[Npositron];
00811   rold = Event.Rpositron+Rdim*Npositron;
00812   rnew = Rpositron+Rdim*Npositron;
00813   while(rnew>Rpositron)*--rnew=*--rold;
00814 
00815   Rgenpos = new double[Npositron];
00816   rold = Event.Rgenpos+Rdim*Npositron;
00817   rnew = Rgenpos+Rdim*Npositron;
00818   while(rnew>Rgenpos)*--rnew=*--rold;
00819 
00820   Rneutron = new double[Nneutron];
00821   rold = Event.Rneutron+Rdim*Nneutron;
00822   rnew = Rneutron+Rdim*Nneutron;
00823   while(rnew>Rneutron)*--rnew=*--rold;
00824 
00825   Rgenneu = new double[Nneutron];
00826   rold = Event.Rgenneu+Rdim*Nneutron;
00827   rnew = Rgenneu+Rdim*Nneutron;
00828   while(rnew>Rgenneu)*--rnew=*--rold;
00829 
00830   Rgamma = new double[Ngamma];
00831   rold = Event.Rgamma+Rdim*Ngamma;
00832   rnew = Rgamma+Rdim*Ngamma;
00833   while(rnew>Rgamma)*--rnew=*--rold;
00834 }    
00835 
00836 ReactorEvent& ReactorEvent::operator=(const ReactorEvent& rhs){
00837 
00838   if(this == &rhs)return *this;
00839 
00840   XCorder = rhs.XCorder;  
00841   Emin = rhs.Emin;
00842   Emax = rhs.Emax;
00843   Rmaxgen = rhs.Rmaxgen;
00844 
00845   XCWeight = rhs.XCWeight;
00846   FluxWeight = rhs.FluxWeight;
00847   
00848   Nnubar    = rhs.Nnubar;
00849   Nelectron = rhs.Nelectron;
00850   Npositron = rhs.Npositron;
00851   Nneutron  = rhs.Nneutron;
00852   Ngamma    = rhs.Ngamma;
00853 
00854   double* pold;
00855   double* pnew;
00856 
00857   delete[] Pnubar;
00858   Pnubar = new double[Nnubar];
00859   pold = rhs.Pnubar+Pdim*Nnubar;
00860   pnew = Pnubar+Pdim*Nnubar;
00861   while(pnew>Pnubar)*--pnew=*--pold;
00862 
00863   delete[] Pelectron;
00864   Pelectron = new double[Nelectron];
00865   pold = rhs.Pelectron+Pdim*Nelectron;
00866   pnew = Pelectron+Pdim*Nelectron;
00867   while(pnew>Pelectron)*--pnew=*--pold;
00868 
00869   delete[] Ppositron;
00870   Ppositron = new double[Npositron];
00871   pold = rhs.Ppositron+Pdim*Npositron;
00872   pnew = Ppositron+Pdim*Npositron;
00873   while(pnew>Ppositron)*--pnew=*--pold;
00874 
00875   delete[] Pneutron;
00876   Pneutron = new double[Nneutron];
00877   pold = rhs.Pneutron+Pdim*Nneutron;
00878   pnew = Pneutron+Pdim*Nneutron;
00879   while(pnew>Pneutron)*--pnew=*--pold;
00880 
00881   delete[] Pgamma;
00882   Pgamma = new double[Ngamma];
00883   pold = rhs.Pgamma+Pdim*Ngamma;
00884   pnew = Pgamma+Pdim*Ngamma;
00885   while(pnew>Pgamma)*--pnew=*--pold;
00886 
00887   double* rold;
00888   double* rnew;
00889 
00890   delete[] Rnubar;
00891   Rnubar = new double[Nnubar];
00892   rold = rhs.Rnubar+Rdim*Nnubar;
00893   rnew = Rnubar+Rdim*Nnubar;
00894   while(rnew>Rnubar)*--rnew=*--rold;
00895 
00896   delete[] Relectron;
00897   Relectron = new double[Nelectron];
00898   rold = rhs.Relectron+Rdim*Nelectron;
00899   rnew = Relectron+Rdim*Nelectron;
00900   while(rnew>Relectron)*--rnew=*--rold;
00901 
00902   delete[] Rgenele;
00903   Rgenele = new double[Nelectron];
00904   rold = rhs.Rgenele+Rdim*Nelectron;
00905   rnew = Rgenele+Rdim*Nelectron;
00906   while(rnew>Rgenele)*--rnew=*--rold;
00907 
00908   delete[] Rpositron;
00909   Rpositron = new double[Npositron];
00910   rold = rhs.Rpositron+Rdim*Npositron;
00911   rnew = Rpositron+Rdim*Npositron;
00912   while(rnew>Rpositron)*--rnew=*--rold;
00913 
00914   delete[] Rgenpos;
00915   Rgenpos = new double[Npositron];
00916   rold = rhs.Rgenpos+Rdim*Npositron;
00917   rnew = Rgenpos+Rdim*Npositron;
00918   while(rnew>Rgenpos)*--rnew=*--rold;
00919 
00920   delete[] Rneutron;
00921   Rneutron = new double[Nneutron];
00922   rold = rhs.Rneutron+Rdim*Nneutron;
00923   rnew = Rneutron+Rdim*Nneutron;
00924   while(rnew>Rneutron)*--rnew=*--rold;
00925 
00926   delete[] Rgenneu;
00927   Rgenneu = new double[Nneutron];
00928   rold = rhs.Rgenneu+Rdim*Nneutron;
00929   rnew = Rgenneu+Rdim*Nneutron;
00930   while(rnew>Rgenneu)*--rnew=*--rold;
00931 
00932   delete[] Rgamma;
00933   Rgamma = new double[Ngamma];
00934   rold = rhs.Rgamma+Rdim*Ngamma;
00935   rnew = Rgamma+Rdim*Ngamma;
00936   while(rnew>Rgamma)*--rnew=*--rold;
00937 
00938   return *this;
00939 }