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