#include <ReactorEvent.hh>
Public Methods | |
| ReactorEvent (int newNevents, double newRmaxgen=k.Rmaxgen, double newEmin=k.Emin, double newEmax=k.Emax, int newXCorder=k.XCorder) | |
| ReactorEvent constructor. More... | |
| ~ReactorEvent () | |
| ReactorEvent destructor. | |
| ReactorEvent (const ReactorEvent &Event) | |
| ReactorEvent copy constructor. | |
| ReactorEvent & | operator= (const ReactorEvent &rhs) |
| ReactorEvent overloaded = operator. | |
| double | XCReweight () |
| double | XCReweight (double Enu) |
| double | XCReweight (ReactorEvent &Event) |
| std::string | GetPmtRad () const |
| Returns the PMT radiation state: ON/OFF. | |
| int | GetNevent () const |
| Returns the number of events. | |
| int | GetNumNubar () const |
| Returns the number of anti-neutrinos. | |
| int | GetNumPositron () const |
| Returns the number of positrons. | |
| int | GetNumElectron () const |
| Returns the number of electrons. | |
| int | GetNumNeutron () const |
| Returns the number of neutrons. | |
| int | GetNumGamma () const |
| Returns the number of gammas. | |
| double | GetXCWeight () const |
| double | GetFluxWeight () const |
| double | GetEnu (int &n) const |
| Returns the energy of each anti-neutrino produced in the event. More... | |
| double | GetNubarParentProcess (int &n) const |
| Returns the parent process used to generate each anti-neutrino in the event. More... | |
| double | GetElectronKE (int &n) const |
| Returns the energy of each electron produced in the event. More... | |
| double | GetElectronCosTheta (int &n) const |
| double | GetElectronPhi (int &n) const |
| double | GetElectronParentProcess (int &n) const |
| Returns the parent process used to generate each electron in the event. More... | |
| double *const | GetElectronPxyz () |
| double | GetPositronKE (int &n) const |
| Returns the energy of each positron produced in the event. More... | |
| double | GetPositronCosTheta (int &n) const |
| double | GetPositronPhi (int &n) const |
| double | GetPositronParentProcess (int &n) const |
| Returns the parent process used to generate each positron in the event. More... | |
| double *const | GetPositronPxyz () |
| double | GetNeutronKE (int &n) const |
| Returns the energy of each neutron produced in the event. More... | |
| double | GetNeutronCosTheta (int &n) const |
| double | GetNeutronPhi (int &n) const |
| double | GetNeutronParentProcess (int &n) const |
| Returns the parent process used to generate each neutron in the event. More... | |
| double *const | GetNeutronPxyz () |
| double | GetGammaKE (int &n) const |
| Returns the energy of each photon produced in the event. More... | |
| double | GetGammaParentProcess (int &n) const |
| Returns the parent process used to generate each photon in the event. More... | |
| double | GetNubarTime (int &n) const |
| double | GetElectronTime (int &n) const |
| double | GetPositronTime (int &n) const |
| double | GetNeutronTime (int &n) const |
| double | GetGammaTime (int &n) const |
| double | GetNubarVertexR (int &n) const |
| double | GetNubarVertexCosTheta (int &n) const |
| double | GetNubarVertexPhi (int &n) const |
| double *const | GetNubarVertexXYZ () |
| double | GetElectronVertexR (int &n) const |
| double | GetElectronVertexCosTheta (int &n) const |
| double | GetElectronVertexPhi (int &n) const |
| double *const | GetElectronVertexXYZ () |
| double | GetPositronVertexR (int &n) const |
| double | GetPositronVertexCosTheta (int &n) const |
| double | GetPositronVertexPhi (int &n) const |
| double *const | GetPositronVertexXYZ () |
| double | GetNeutronVertexR (int &n) const |
| double | GetNeutronVertexCosTheta (int &n) const |
| double | GetNeutronVertexPhi (int &n) const |
| double *const | GetNeutronVertexXYZ () |
| double | GetGammaVertexR (int &n) const |
| double | GetGammaVertexCosTheta (int &n) const |
| double | GetGammaVertexPhi (int &n) const |
| double *const | GetGammaVertexXYZ () |
| double | GetGeneratedElectronTime (int &n) const |
| double | GetGeneratedElectronVertexR (int &n) const |
| double | GetGeneratedElectronVertexCosTheta (int &n) const |
| double | GetGeneratedElectronVertexPhi (int &n) const |
| double *const | GetGeneratedElectronVertexXYZ () |
| double | GetGeneratedPositronTime (int &n) const |
| double | GetGeneratedPositronVertexR (int &n) const |
| double | GetGeneratedPositronVertexCosTheta (int &n) const |
| double | GetGeneratedPositronVertexPhi (int &n) const |
| double *const | GetGeneratedPositronVertexXYZ () |
| double | GetGeneratedNeutronTime (int &n) const |
| double | GetGeneratedNeutronVertexR (int &n) const |
| double | GetGeneratedNeutronVertexCosTheta (int &n) const |
| double | GetGeneratedNeutronVertexPhi (int &n) const |
| double *const | GetGeneratedNeutronVertexXYZ () |
| void | SetPositronVertex (int &n, double *rpos) |
| void | SetNeutronVertex (int &n, double *rneu) |
Public Attributes | |
| std::string | eventType |
| std::string | pmtRad |
Static Public Attributes | |
| const int | Pdim = 7 |
| const int | Rdim = 7 |
It first sets up the event by reading in the type of event from event_setup.txt, the event setup file. The three types of event currently supported are "reactor," "cf252," and "all." For more on event_setup.txt, see event_readme.txt.
The "reactor" setting generates an anti-neutrino, positron, neutron from the inverse-beta decay reaction: nubar + proton -> positron + neutron. The "cf252" setting generates a neutron from a Cf252 decay. The "all" setting mixes both types of event.
Definition at line 26 of file ReactorEvent.hh.
|
||||||||||||||||||||||||
|
ReactorEvent constructor. Called with the number of generated events set in ReactorFsim.cpp, the maximum allowed radius for generating particles in the simulation, Rmaxgen (default = 200 cm), and the minimum and maximum allowed energies, respectively, for the generated nubar in the inverse-beta decay reaction. The default minimum energy is 1.85; the maximum is 8.15, set in ReactorConstants.hh. Definition at line 20 of file ReactorEvent.cpp. References ReactorConstants::Delta, MyCfSource::GetCfGammaEnergy, MyCfSource::GetCfGammaTime, MyCfSource::GetCfNeutronEnergy, MyCfSource::GetCfNeutronTime, MyCfSource::GetCfVertexCosTheta, MyCfSource::GetCfVertexPhi, MyCfSource::GetCfVertexR, MuonPropagator::GetElectronKE, MuonPropagator::GetElectronTime, MuonPropagator::GetElectronVertexXYZ, MuonPropagator::GetGammaKE, MuonPropagator::GetGammaTime, MuonPropagator::GetGammaVertexXYZ, MuonPropagator::GetNeutronKE, MuonPropagator::GetNeutronTime, MuonPropagator::GetNeutronVertexXYZ, MyCfSource::GetNumCfSources, MuonPropagator::GetNumElectron, MuonPropagator::GetNumGamma, MyCfSource::GetNumGamma, MuonPropagator::GetNumNeutron, MyCfSource::GetNumNeutron, MuonPropagator::GetWeight, ReactorConstants::Melectron, ReactorConstants::Mproton, ReactorConstants::Ndim, and ReactorConstants::pi.
00024 {
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 }
|
|
|
Returns the energy of each electron produced in the event. Called with the integer index for each electron in the Pelectron array, from 0 to the total number of electrons. Definition at line 99 of file ReactorEvent.hh.
00099 {return Pelectron[4+n*Pdim];}
|
|
|
Returns the parent process used to generate each electron in the event. Called with the integer index for each electron in the Pelectron array, from 0 to the total number of electrons. Definition at line 106 of file ReactorEvent.hh.
00106 {return Pelectron[6+n*Pdim];}
|
|
|
Returns the energy of each anti-neutrino produced in the event. Called with the integer index for each anti-neutrino in the Pnubar array, from 0 to the total number of anti-neutrinos. Definition at line 88 of file ReactorEvent.hh.
00088 {return Pnubar[0+n*Pdim];}
|
|
|
Returns the energy of each photon produced in the event. Called with the integer index for each photon in the Pgamma array, from 0 to the total number of photons. Definition at line 141 of file ReactorEvent.hh.
00141 {return Pgamma[0+n*Pdim];}
|
|
|
Returns the parent process used to generate each photon in the event. Called with the integer index for each photon in the Pgamma array, from 0 to the total number of photons. Definition at line 146 of file ReactorEvent.hh.
00146 {return Pgamma[6+n*Pdim];}
|
|
|
Returns the energy of each neutron produced in the event. Called with the integer index for each neutron in the Pneutron array, from 0 to the total number of neutrons. Definition at line 127 of file ReactorEvent.hh.
00127 {return Pneutron[4+n*Pdim];}
|
|
|
Returns the parent process used to generate each neutron in the event. Called with the integer index for each neutron in the Pneutron array, from 0 to the total number of neutrons. Definition at line 134 of file ReactorEvent.hh.
00134 {return Pneutron[6+n*Pdim];}
|
|
|
Returns the parent process used to generate each anti-neutrino in the event. Called with the integer index for each anti-neutrino in the Pnubar array, from 0 to the total number of anti-neutrinos. Definition at line 93 of file ReactorEvent.hh.
00093 {return Pnubar[6+n*Pdim];}
|
|
|
Returns the energy of each positron produced in the event. Called with the integer index for each positron in the Ppositron array, from 0 to the total number of positrons. Definition at line 113 of file ReactorEvent.hh.
00113 {return Ppositron[4+n*Pdim];}
|
|
|
Returns the parent process used to generate each positron in the event. Called with the integer index for each positron in the Ppositron array, from 0 to the total number of positrons. Definition at line 120 of file ReactorEvent.hh.
00120 {return Ppositron[6+n*Pdim];}
|
1.2.14 written by Dimitri van Heesch,
© 1997-2002