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