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