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 
00018 ReactorEvent::ReactorEvent(double newNevents,
00019                            double newRmaxgen,
00020                            double newEmin,
00021                            double newEmax,
00022                            int newXCorder){
00023   Nevents = newNevents;
00024   Rmaxgen = newRmaxgen;
00025   Emin = newEmin;
00026   Emax = newEmax;
00027   XCorder = newXCorder;
00028 
00029   int nmax = k.Ndim;
00030   int nnu  = 0;
00031   int npos = 0;
00032   int nneu = 0;
00033 
00034   double* pnubar = new double[nmax];
00035   double* ppositron = new double[nmax];
00036   double* pneutron = new double[nmax];
00037 
00038   double* pnu  = pnubar;
00039   double* ppos = ppositron;
00040   double* pneu = pneutron;
00041 
00042   double* rnubar = new double[nmax];
00043   double* rpositron = new double[nmax];
00044   double* rgenpos = new double[nmax];
00045   double* rneutron = new double[nmax];
00046   double* rgenneu = new double[nmax];
00047 
00048   double* rnu  = rnubar;
00049   double* rpos = rpositron;
00050   double* rgp = rgenpos;
00051   double* rneu = rneutron;
00052   double* rgn = rgenneu;
00053 
00054   double Enu,z,phi;
00055   double R,zR,phiR;
00056 
00057   // setup the event
00058   fstream fin("src/event_setup.txt", ios::in);
00059   //cout << "setup file: src/event_setup.txt" << endl;
00060 
00061   char paramLine[256];
00062   string dummy;
00063   bool done = false;
00064   bool set = false;
00065 
00066   while(!done){
00067 
00068     // '#' specifies a comment line -- ignore it
00069     if(fin.peek() == '#'){
00070 
00071       fin.getline(paramLine, 256);
00072       continue;
00073     }
00074 
00075     // 'E'ND or 'e'nd specifies end of file -- set done flag to true
00076     if(fin.peek() == 'e' || fin.peek() == 'E'){
00077 
00078       done = true;
00079       continue;
00080     }
00081 
00082     // 'T'YPE or 't'ype specifies the event type
00083     if(fin.peek() == 't' || fin.peek() == 'T'){
00084 
00085       if(set){
00086 
00087         fin.getline( paramLine, 256 );
00088         continue;
00089       }
00090       fin >> dummy;
00091       while(fin.peek() == ' ' || fin.peek() == '=') 
00092         fin.ignore(1);
00093 
00094       fin >> eventType;
00095 
00096       if(eventType == "reactor" || eventType == "cf252" || eventType == "all"){
00097         set = true;
00098       }
00099       else{
00100         cout << "  SETUP ERROR: unrecognized event type!" << endl;
00101       }
00102 
00103       continue;
00104     }
00105 
00106     // 'P'MTRAD or 'p'mtrad turns on/off radioactive decays in the PMTs
00107     if(fin.peek() == 'p' || fin.peek() == 'P'){
00108 
00109       fin >> dummy;
00110       while(fin.peek() == ' ' || fin.peek() == '=') 
00111         fin.ignore(1);
00112 
00113       fin >> pmtRad;
00114       continue;
00115     }
00116 
00117     // ignore extra whitespace
00118     if(fin.peek() == ' ' || fin.peek() == '\n'){
00119 
00120       fin.ignore(1);
00121       continue;
00122     }
00123   }
00124 
00125   // if all events are requested, mix the 'reactor' and 'cf252' 
00126   // events together with 90/10% Cf252/reactor events
00127   if(eventType == "all"){
00128     int len = 1;
00129     float r[len];
00130     ranlux_(r,len);
00131 
00132     if(r[0] < 0.90){
00133       eventType = "cf252";
00134     }
00135     else{
00136       eventType = "reactor";
00137     }
00138   }
00139   //cout << eventType << " type of events" << endl;
00140   //cout << "PMT radiation " << pmtRad << endl;
00141 
00142   if(eventType == "reactor"){
00143 
00144     XCWeight = 1;
00145 
00146     // inverse beta decay reaction
00147     double XCmax = InverseBeta(Emax,-1);
00148     double FluxMax = RMPflux(Emin);
00149 
00150     bool passed=0;
00151 
00152     while(!passed){
00153       int len = 7;
00154       float r[len];
00155       ranlux_(r,len);
00156       Enu = Emin+(Emax-Emin)*r[0];
00157       z = -1+2*r[1];
00158       phi = 2*k.pi*r[2];
00159 
00160       R = Rmaxgen*exp(log(r[3])/3);
00161       zR = -1+2*r[4];
00162       phiR = 2*k.pi*r[5];
00163 
00164       float XCtest = XCmax*FluxMax*r[6];
00165       XCWeight = InverseBeta(Enu,z);
00166       FluxWeight=RMPflux(Enu);
00167       passed= XCWeight*FluxWeight>XCtest;
00168     }
00169   
00170     double E0 = Enu - k.Delta;
00171     double p0 = sqrt(E0*E0-k.Melectron*k.Melectron);
00172     double v0 = p0/E0;
00173     double Ysquared = (k.Delta*k.Delta-k.Melectron*k.Melectron)/2;
00174     double E1 = E0*(1-Enu/k.Mproton*(1-v0*z))-Ysquared/k.Mproton;
00175     double p1 = sqrt(E1*E1-k.Melectron*k.Melectron);
00176 
00177     *pnu = Enu; pnu++;
00178     *pnu = 0; pnu++;
00179     *pnu = 0; pnu++;
00180     *pnu = Enu; pnu++;
00181     *pnu = Enu; pnu++;
00182     *pnu = Enu; pnu++;
00183     *pnu = 1.; pnu++; // inverse beta decay process Id = 1
00184 
00185     *ppos = E1; ppos++;
00186     *ppos = p1*sqrt(1-z*z)*cos(phi); ppos++;
00187     *ppos = p1*sqrt(1-z*z)*sin(phi); ppos++;
00188     *ppos = p1*z; ppos++;
00189     *ppos = E1-k.Melectron; ppos++;
00190     *ppos = p1; ppos++;
00191     *ppos = 1.; ppos++;
00192 
00193     *pneu = Enu+k.Mproton-E1; pneu++;
00194     *pneu = -p1*sqrt(1-z*z)*cos(phi); pneu++;
00195     *pneu = -p1*sqrt(1-z*z)*sin(phi); pneu++;
00196     *pneu = Enu-p1*z; pneu++;
00197     *pneu = Enu-E1-k.Delta; pneu++;
00198     *pneu = sqrt((Enu+k.Mproton-E1)*(Enu+k.Mproton-E1)-
00199                          (k.Mproton+k.Delta)*(k.Mproton+k.Delta)); pneu++;
00200     *pneu = 1.; pneu++;
00201 
00202     *rnu = 0; rnu++;
00203     *rnu = R*sqrt(1-zR*zR)*cos(phiR); rnu++;
00204     *rnu = R*sqrt(1-zR*zR)*sin(phiR); rnu++;
00205     *rnu = R*zR; rnu++;
00206     *rnu = R; rnu++;
00207     *rnu = zR; rnu++;
00208     *rnu = phiR; rnu++;
00209 
00210     *rpos = 0; rpos++;
00211     *rpos = R*sqrt(1-zR*zR)*cos(phiR); rpos++;
00212     *rpos = R*sqrt(1-zR*zR)*sin(phiR); rpos++;
00213     *rpos = R*zR; rpos++;
00214     *rpos = R; rpos++;
00215     *rpos = zR; rpos++;
00216     *rpos = phiR; rpos++;
00217 
00218     *rgp = 0; rgp++;
00219     *rgp = R*sqrt(1-zR*zR)*cos(phiR); rgp++;
00220     *rgp = R*sqrt(1-zR*zR)*sin(phiR); rgp++;
00221     *rgp = R*zR; rgp++;
00222     *rgp = R; rgp++;
00223     *rgp = zR; rgp++;
00224     *rgp = phiR; rgp++;
00225 
00226     *rneu = 0; rneu++;
00227     *rneu = R*sqrt(1-zR*zR)*cos(phiR); rneu++;
00228     *rneu = R*sqrt(1-zR*zR)*sin(phiR); rneu++;
00229     *rneu = R*zR; rneu++;
00230     *rneu = R; rneu++;
00231     *rneu = zR; rneu++;
00232     *rneu = phiR; rneu++;
00233 
00234     *rgn = 0; rgn++;
00235     *rgn = R*sqrt(1-zR*zR)*cos(phiR); rgn++;
00236     *rgn = R*sqrt(1-zR*zR)*sin(phiR); rgn++;
00237     *rgn = R*zR; rgn++;
00238     *rgn = R; rgn++;
00239     *rgn = zR; rgn++;
00240     *rgn = phiR; rgn++;
00241 
00242     // count the nubar, positron, and neutron
00243     nnu++;
00244     npos++;
00245     nneu++;
00246 
00247   } // done with reactor neutrino inverse-beta reaction
00248 
00249   // Cf source reaction
00250   if(eventType == "cf252"){
00251 
00252     int Isotope = 252;
00253     MyCfSource CfSource(Isotope,Rmaxgen);
00254 
00255     int Nsources = CfSource.GetNumCfSources();
00256     //cout << "Number of Cf252 sources = " << Nsources << endl;
00257 
00258     for(int n=0;n<Nsources;n++){
00259       float neutronKE = CfSource.GetCfNeutronEnergy(n);
00260       //cout << "   Cf252 source " << n << " energy = " 
00261       //           << neutronKE << endl;
00262 
00263       float Mneutron = k.Mproton+k.Delta;
00264       float neutronE = neutronKE+Mneutron;
00265 
00266       int len = 2;
00267       float r[len];
00268       ranlux_(r,len);
00269       z = -1+2*r[0];
00270       phi = 2*k.pi*r[1];
00271       double x = sqrt(1-z*z)*cos(phi);
00272       double y = sqrt(1-z*z)*sin(phi);
00273 
00274       *pneu = neutronE; pneu++;
00275       *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*x; pneu++;
00276       *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*y; pneu++;
00277       *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*z; pneu++;
00278       *pneu = neutronKE; pneu++;
00279       *pneu = sqrt(neutronE*neutronE-Mneutron*Mneutron); pneu++;
00280       *pneu = 2.; pneu++; // cf252 source process Id = 2
00281 
00282       R = CfSource.GetCfVertexR(n);
00283       zR = CfSource.GetCfVertexCosTheta(n);
00284       phiR = CfSource.GetCfVertexPhi(n);
00285 
00286       *rneu = 0; rneu++;
00287       *rneu = R*sqrt(1-zR*zR)*cos(phiR); rneu++;
00288       *rneu = R*sqrt(1-zR*zR)*sin(phiR); rneu++;
00289       *rneu = R*zR; rneu++;
00290       *rneu = R; rneu++;
00291       *rneu = zR; rneu++;
00292       *rneu = phiR; rneu++;
00293 
00294       *rgn = 0; rgn++;
00295       *rgn = R*sqrt(1-zR*zR)*cos(phiR); rgn++;
00296       *rgn = R*sqrt(1-zR*zR)*sin(phiR); rgn++;
00297       *rgn = R*zR; rgn++;
00298       *rgn = R; rgn++;
00299       *rgn = zR; rgn++;
00300       *rgn = phiR; rgn++;
00301     }
00302 
00303     // count the neutrons
00304     nneu += Nsources;
00305 
00306   } // done with cf source neutron production
00307 
00308   //
00309   // transfer temp data to storage
00310   //
00311   Nnubar    = nnu;
00312   Npositron = npos;
00313   Nneutron  = nneu;
00314 
00315   Pnubar    = new double[Pdim*nnu];
00316   Ppositron = new double[Pdim*npos];
00317   Pneutron  = new double[Pdim*nneu];
00318 
00319   Rnubar    = new double[Rdim*nnu];
00320   Rpositron = new double[Rdim*npos];
00321   Rgenpos   = new double[Rdim*npos];
00322   Rneutron  = new double[Rdim*nneu];
00323   Rgenneu   = new double[Rdim*nneu];
00324 
00325   double* ptr;
00326 
00327   ptr = Pnubar+Pdim*Nnubar;
00328   while(ptr>Pnubar)*--ptr=*--pnu;
00329 
00330   ptr = Ppositron+Pdim*Npositron;
00331   while(ptr>Ppositron)*--ptr=*--ppos;
00332 
00333   ptr = Pneutron+Pdim*Nneutron;
00334   while(ptr>Pneutron)*--ptr=*--pneu;
00335 
00336   ptr = Rnubar+Rdim*Nnubar;
00337   while(ptr>Rnubar)*--ptr=*--rnu;
00338 
00339   ptr = Rpositron+Rdim*Npositron;
00340   while(ptr>Rpositron)*--ptr=*--rpos;
00341 
00342   ptr = Rgenpos+Rdim*Npositron;
00343   while(ptr>Rgenpos)*--ptr=*--rgp;
00344 
00345   ptr = Rneutron+Rdim*Nneutron;
00346   while(ptr>Rneutron)*--ptr=*--rneu;
00347 
00348   ptr = Rgenneu+Rdim*Nneutron;
00349   while(ptr>Rgenneu)*--ptr=*--rgn;
00350 
00351   /*
00352   cout << Nnubar << " nubar(s):" << endl;
00353   for(int i=0;i<Nnubar*Pdim;i++){
00354     cout << " Pnubar[" << i << "] " << Pnubar[i] << endl;
00355   }
00356   for(int i=0;i<Nnubar*Rdim;i++){
00357     cout << "  Rnubar[" << i << "] " << Rnubar[i] << endl;
00358   }
00359 
00360   cout << Npositron << " positron(s):" << endl;
00361   for(int i=0;i<Npositron*Pdim;i++){
00362     cout << " Ppositron[" << i << "] " << Ppositron[i] << endl;
00363   }
00364   for(int i=0;i<Npositron*Rdim;i++){
00365     cout << "  Rpositron[" << i << "] " << Rpositron[i] << endl;
00366   }
00367   for(int i=0;i<Npositron*Rdim;i++){
00368     cout << "  Rgenpos[" << i << "] " << Rgenpos[i] << endl;
00369   }
00370 
00371   cout << Nneutron << " neutron(s):" << endl;
00372   for(int i=0;i<Nneutron*Pdim;i++){
00373     cout << " Pneutron[" << i << "] " << Pneutron[i] << endl;
00374   }
00375   for(int i=0;i<Nneutron*Rdim;i++){
00376     cout << "  Rneutron[" << i << "] " << Rneutron[i] << endl;
00377   }
00378   for(int i=0;i<Nneutron*Rdim;i++){
00379     cout << "  Rgenneu[" << i << "] " << Rgenneu[i] << endl;
00380   }
00381   */
00382 
00383 }
00384 
00385 ReactorEvent::~ReactorEvent(){
00386   delete[] Pnubar;
00387   delete[] Rnubar;
00388   delete[] Ppositron;
00389   delete[] Rpositron;
00390   delete[] Rgenpos;
00391   delete[] Pneutron;
00392   delete[] Rneutron;
00393   delete[] Rgenneu;
00394 }
00395 
00396 ReactorEvent::ReactorEvent(const ReactorEvent& Event){
00397 
00398   XCorder = Event.XCorder;  
00399   Emin = Event.Emin;
00400   Emax = Event.Emax;
00401   Rmaxgen = Event.Rmaxgen;
00402 
00403   XCWeight = Event.XCWeight;
00404   FluxWeight = Event.FluxWeight;
00405 
00406   Nnubar    = Event.Nnubar;
00407   Npositron = Event.Npositron;
00408   Nneutron  = Event.Nneutron;
00409 
00410   double* pold;
00411   double* pnew;
00412 
00413   Pnubar = new double[Nnubar];
00414   pold = Event.Pnubar+Pdim*Nnubar;
00415   pnew = Pnubar+Pdim*Nnubar;
00416   while(pnew>Pnubar)*--pnew=*--pold;
00417 
00418   Ppositron = new double[Npositron];
00419   pold = Event.Ppositron+Pdim*Npositron;
00420   pnew = Ppositron+Pdim*Npositron;
00421   while(pnew>Ppositron)*--pnew=*--pold;
00422 
00423   Pneutron = new double[Nneutron];
00424   pold = Event.Pneutron+Pdim*Nneutron;
00425   pnew = Pneutron+Pdim*Nneutron;
00426   while(pnew>Pneutron)*--pnew=*--pold;
00427 
00428   double* rold;
00429   double* rnew;
00430 
00431   Rnubar = new double[Nnubar];
00432   rold = Event.Rnubar+Rdim*Nnubar;
00433   rnew = Rnubar+Rdim*Nnubar;
00434   while(rnew>Rnubar)*--rnew=*--rold;
00435 
00436   Rpositron = new double[Npositron];
00437   rold = Event.Rpositron+Rdim*Npositron;
00438   rnew = Rpositron+Rdim*Npositron;
00439   while(rnew>Rpositron)*--rnew=*--rold;
00440 
00441   Rgenpos = new double[Npositron];
00442   rold = Event.Rgenpos+Rdim*Npositron;
00443   rnew = Rgenpos+Rdim*Npositron;
00444   while(rnew>Rgenpos)*--rnew=*--rold;
00445 
00446   Rneutron = new double[Nneutron];
00447   rold = Event.Rneutron+Rdim*Nneutron;
00448   rnew = Rneutron+Rdim*Nneutron;
00449   while(rnew>Rneutron)*--rnew=*--rold;
00450 
00451   Rgenneu = new double[Nneutron];
00452   rold = Event.Rgenneu+Rdim*Nneutron;
00453   rnew = Rgenneu+Rdim*Nneutron;
00454   while(rnew>Rgenneu)*--rnew=*--rold;
00455 
00456 }    
00457 
00458 ReactorEvent& ReactorEvent::operator=(const ReactorEvent& rhs){
00459 
00460   if(this == &rhs)return *this;
00461 
00462   XCorder = rhs.XCorder;  
00463   Emin = rhs.Emin;
00464   Emax = rhs.Emax;
00465   Rmaxgen = rhs.Rmaxgen;
00466 
00467   XCWeight = rhs.XCWeight;
00468   FluxWeight = rhs.FluxWeight;
00469   
00470   Nnubar    = rhs.Nnubar;
00471   Npositron = rhs.Npositron;
00472   Nneutron  = rhs.Nneutron;
00473 
00474   double* pold;
00475   double* pnew;
00476 
00477   delete[] Pnubar;
00478   Pnubar = new double[Nnubar];
00479   pold = rhs.Pnubar+Pdim*Nnubar;
00480   pnew = Pnubar+Pdim*Nnubar;
00481   while(pnew>Pnubar)*--pnew=*--pold;
00482 
00483   delete[] Ppositron;
00484   Ppositron = new double[Npositron];
00485   pold = rhs.Ppositron+Pdim*Npositron;
00486   pnew = Ppositron+Pdim*Npositron;
00487   while(pnew>Ppositron)*--pnew=*--pold;
00488 
00489   delete[] Pneutron;
00490   Pneutron = new double[Nneutron];
00491   pold = rhs.Pneutron+Pdim*Nneutron;
00492   pnew = Pneutron+Pdim*Nneutron;
00493   while(pnew>Pneutron)*--pnew=*--pold;
00494 
00495   double* rold;
00496   double* rnew;
00497 
00498   delete[] Rnubar;
00499   Rnubar = new double[Nnubar];
00500   rold = rhs.Rnubar+Rdim*Nnubar;
00501   rnew = Rnubar+Rdim*Nnubar;
00502   while(rnew>Rnubar)*--rnew=*--rold;
00503 
00504   delete[] Rpositron;
00505   Rpositron = new double[Npositron];
00506   rold = rhs.Rpositron+Rdim*Npositron;
00507   rnew = Rpositron+Rdim*Npositron;
00508   while(rnew>Rpositron)*--rnew=*--rold;
00509 
00510   delete[] Rgenpos;
00511   Rgenpos = new double[Npositron];
00512   rold = rhs.Rgenpos+Rdim*Npositron;
00513   rnew = Rgenpos+Rdim*Npositron;
00514   while(rnew>Rgenpos)*--rnew=*--rold;
00515 
00516   delete[] Rneutron;
00517   Rneutron = new double[Nneutron];
00518   rold = rhs.Rneutron+Rdim*Nneutron;
00519   rnew = Rneutron+Rdim*Nneutron;
00520   while(rnew>Rneutron)*--rnew=*--rold;
00521 
00522   delete[] Rgenneu;
00523   Rgenneu = new double[Nneutron];
00524   rold = rhs.Rgenneu+Rdim*Nneutron;
00525   rnew = Rgenneu+Rdim*Nneutron;
00526   while(rnew>Rgenneu)*--rnew=*--rold;
00527 
00528   return *this;
00529 }