#include <ReactorEvent.hh>
Public Methods | |
| ReactorEvent (double newNevents=0, 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 | GetNevents () const |
| Returns the number of events. | |
| int | GetNumNubar () const |
| Returns the number of anti-neutrinos. | |
| int | GetNumPositron () const |
| Returns the number of positrons. | |
| int | GetNumNeutron () const |
| Returns the number of neutrons. | |
| 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 | 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 | GetNubarTime (int &n) const |
| double | GetPositronTime (int &n) const |
| double | GetNeutronTime (int &n) const |
| double | GetNubarVertexR (int &n) const |
| double | GetNubarVertexCosTheta (int &n) const |
| double | GetNubarVertexPhi (int &n) const |
| double *const | GetNubarVertexXYZ () |
| 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 | 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 18 of file ReactorEvent.cpp. References ReactorConstants::Delta, MyCfSource::GetCfNeutronEnergy, MyCfSource::GetCfVertexCosTheta, MyCfSource::GetCfVertexPhi, MyCfSource::GetCfVertexR, MyCfSource::GetNumCfSources, ReactorConstants::Melectron, ReactorConstants::Mproton, ReactorConstants::Ndim, and ReactorConstants::pi.
00022 {
00023 Nevents = newNevents;
00024 Rmaxgen = newRmaxgen;
00025 Emin = newEmin;
00026 Emax = newEmax;
00027 XCorder = newXCorder;
00028
00029 int nmax = k.Ndim;
00030 int nnu = 0;
00031 int npos = 0;
00032 int nneu = 0;
00033
00034 double* pnubar = new double[nmax];
00035 double* ppositron = new double[nmax];
00036 double* pneutron = new double[nmax];
00037
00038 double* pnu = pnubar;
00039 double* ppos = ppositron;
00040 double* pneu = pneutron;
00041
00042 double* rnubar = new double[nmax];
00043 double* rpositron = new double[nmax];
00044 double* rgenpos = new double[nmax];
00045 double* rneutron = new double[nmax];
00046 double* rgenneu = new double[nmax];
00047
00048 double* rnu = rnubar;
00049 double* rpos = rpositron;
00050 double* rgp = rgenpos;
00051 double* rneu = rneutron;
00052 double* rgn = rgenneu;
00053
00054 double Enu,z,phi;
00055 double R,zR,phiR;
00056
00057 // setup the event
00058 fstream fin("src/event_setup.txt", ios::in);
00059 //cout << "setup file: src/event_setup.txt" << endl;
00060
00061 char paramLine[256];
00062 string dummy;
00063 bool done = false;
00064 bool set = false;
00065
00066 while(!done){
00067
00068 // '#' specifies a comment line -- ignore it
00069 if(fin.peek() == '#'){
00070
00071 fin.getline(paramLine, 256);
00072 continue;
00073 }
00074
00075 // 'E'ND or 'e'nd specifies end of file -- set done flag to true
00076 if(fin.peek() == 'e' || fin.peek() == 'E'){
00077
00078 done = true;
00079 continue;
00080 }
00081
00082 // 'T'YPE or 't'ype specifies the event type
00083 if(fin.peek() == 't' || fin.peek() == 'T'){
00084
00085 if(set){
00086
00087 fin.getline( paramLine, 256 );
00088 continue;
00089 }
00090 fin >> dummy;
00091 while(fin.peek() == ' ' || fin.peek() == '=')
00092 fin.ignore(1);
00093
00094 fin >> eventType;
00095
00096 if(eventType == "reactor" || eventType == "cf252" || eventType == "all"){
00097 set = true;
00098 }
00099 else{
00100 cout << " SETUP ERROR: unrecognized event type!" << endl;
00101 }
00102
00103 continue;
00104 }
00105
00106 // 'P'MTRAD or 'p'mtrad turns on/off radioactive decays in the PMTs
00107 if(fin.peek() == 'p' || fin.peek() == 'P'){
00108
00109 fin >> dummy;
00110 while(fin.peek() == ' ' || fin.peek() == '=')
00111 fin.ignore(1);
00112
00113 fin >> pmtRad;
00114 continue;
00115 }
00116
00117 // ignore extra whitespace
00118 if(fin.peek() == ' ' || fin.peek() == '\n'){
00119
00120 fin.ignore(1);
00121 continue;
00122 }
00123 }
00124
00125 // if all events are requested, mix the 'reactor' and 'cf252'
00126 // events together with 90/10% Cf252/reactor events
00127 if(eventType == "all"){
00128 int len = 1;
00129 float r[len];
00130 ranlux_(r,len);
00131
00132 if(r[0] < 0.90){
00133 eventType = "cf252";
00134 }
00135 else{
00136 eventType = "reactor";
00137 }
00138 }
00139 //cout << eventType << " type of events" << endl;
00140 //cout << "PMT radiation " << pmtRad << endl;
00141
00142 if(eventType == "reactor"){
00143
00144 XCWeight = 1;
00145
00146 // inverse beta decay reaction
00147 double XCmax = InverseBeta(Emax,-1);
00148 double FluxMax = RMPflux(Emin);
00149
00150 bool passed=0;
00151
00152 while(!passed){
00153 int len = 7;
00154 float r[len];
00155 ranlux_(r,len);
00156 Enu = Emin+(Emax-Emin)*r[0];
00157 z = -1+2*r[1];
00158 phi = 2*k.pi*r[2];
00159
00160 R = Rmaxgen*exp(log(r[3])/3);
00161 zR = -1+2*r[4];
00162 phiR = 2*k.pi*r[5];
00163
00164 float XCtest = XCmax*FluxMax*r[6];
00165 XCWeight = InverseBeta(Enu,z);
00166 FluxWeight=RMPflux(Enu);
00167 passed= XCWeight*FluxWeight>XCtest;
00168 }
00169
00170 double E0 = Enu - k.Delta;
00171 double p0 = sqrt(E0*E0-k.Melectron*k.Melectron);
00172 double v0 = p0/E0;
00173 double Ysquared = (k.Delta*k.Delta-k.Melectron*k.Melectron)/2;
00174 double E1 = E0*(1-Enu/k.Mproton*(1-v0*z))-Ysquared/k.Mproton;
00175 double p1 = sqrt(E1*E1-k.Melectron*k.Melectron);
00176
00177 *pnu = Enu; pnu++;
00178 *pnu = 0; pnu++;
00179 *pnu = 0; pnu++;
00180 *pnu = Enu; pnu++;
00181 *pnu = Enu; pnu++;
00182 *pnu = Enu; pnu++;
00183 *pnu = 1.; pnu++; // inverse beta decay process Id = 1
00184
00185 *ppos = E1; ppos++;
00186 *ppos = p1*sqrt(1-z*z)*cos(phi); ppos++;
00187 *ppos = p1*sqrt(1-z*z)*sin(phi); ppos++;
00188 *ppos = p1*z; ppos++;
00189 *ppos = E1-k.Melectron; ppos++;
00190 *ppos = p1; ppos++;
00191 *ppos = 1.; ppos++;
00192
00193 *pneu = Enu+k.Mproton-E1; pneu++;
00194 *pneu = -p1*sqrt(1-z*z)*cos(phi); pneu++;
00195 *pneu = -p1*sqrt(1-z*z)*sin(phi); pneu++;
00196 *pneu = Enu-p1*z; pneu++;
00197 *pneu = Enu-E1-k.Delta; pneu++;
00198 *pneu = sqrt((Enu+k.Mproton-E1)*(Enu+k.Mproton-E1)-
00199 (k.Mproton+k.Delta)*(k.Mproton+k.Delta)); pneu++;
00200 *pneu = 1.; pneu++;
00201
00202 *rnu = 0; rnu++;
00203 *rnu = R*sqrt(1-zR*zR)*cos(phiR); rnu++;
00204 *rnu = R*sqrt(1-zR*zR)*sin(phiR); rnu++;
00205 *rnu = R*zR; rnu++;
00206 *rnu = R; rnu++;
00207 *rnu = zR; rnu++;
00208 *rnu = phiR; rnu++;
00209
00210 *rpos = 0; rpos++;
00211 *rpos = R*sqrt(1-zR*zR)*cos(phiR); rpos++;
00212 *rpos = R*sqrt(1-zR*zR)*sin(phiR); rpos++;
00213 *rpos = R*zR; rpos++;
00214 *rpos = R; rpos++;
00215 *rpos = zR; rpos++;
00216 *rpos = phiR; rpos++;
00217
00218 *rgp = 0; rgp++;
00219 *rgp = R*sqrt(1-zR*zR)*cos(phiR); rgp++;
00220 *rgp = R*sqrt(1-zR*zR)*sin(phiR); rgp++;
00221 *rgp = R*zR; rgp++;
00222 *rgp = R; rgp++;
00223 *rgp = zR; rgp++;
00224 *rgp = phiR; rgp++;
00225
00226 *rneu = 0; rneu++;
00227 *rneu = R*sqrt(1-zR*zR)*cos(phiR); rneu++;
00228 *rneu = R*sqrt(1-zR*zR)*sin(phiR); rneu++;
00229 *rneu = R*zR; rneu++;
00230 *rneu = R; rneu++;
00231 *rneu = zR; rneu++;
00232 *rneu = phiR; rneu++;
00233
00234 *rgn = 0; rgn++;
00235 *rgn = R*sqrt(1-zR*zR)*cos(phiR); rgn++;
00236 *rgn = R*sqrt(1-zR*zR)*sin(phiR); rgn++;
00237 *rgn = R*zR; rgn++;
00238 *rgn = R; rgn++;
00239 *rgn = zR; rgn++;
00240 *rgn = phiR; rgn++;
00241
00242 // count the nubar, positron, and neutron
00243 nnu++;
00244 npos++;
00245 nneu++;
00246
00247 } // done with reactor neutrino inverse-beta reaction
00248
00249 // Cf source reaction
00250 if(eventType == "cf252"){
00251
00252 int Isotope = 252;
00253 MyCfSource CfSource(Isotope,Rmaxgen);
00254
00255 int Nsources = CfSource.GetNumCfSources();
00256 //cout << "Number of Cf252 sources = " << Nsources << endl;
00257
00258 for(int n=0;n<Nsources;n++){
00259 float neutronKE = CfSource.GetCfNeutronEnergy(n);
00260 //cout << " Cf252 source " << n << " energy = "
00261 // << neutronKE << endl;
00262
00263 float Mneutron = k.Mproton+k.Delta;
00264 float neutronE = neutronKE+Mneutron;
00265
00266 int len = 2;
00267 float r[len];
00268 ranlux_(r,len);
00269 z = -1+2*r[0];
00270 phi = 2*k.pi*r[1];
00271 double x = sqrt(1-z*z)*cos(phi);
00272 double y = sqrt(1-z*z)*sin(phi);
00273
00274 *pneu = neutronE; pneu++;
00275 *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*x; pneu++;
00276 *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*y; pneu++;
00277 *pneu = sqrt((neutronE*neutronE)-(Mneutron*Mneutron))*z; pneu++;
00278 *pneu = neutronKE; pneu++;
00279 *pneu = sqrt(neutronE*neutronE-Mneutron*Mneutron); pneu++;
00280 *pneu = 2.; pneu++; // cf252 source process Id = 2
00281
00282 R = CfSource.GetCfVertexR(n);
00283 zR = CfSource.GetCfVertexCosTheta(n);
00284 phiR = CfSource.GetCfVertexPhi(n);
00285
00286 *rneu = 0; rneu++;
00287 *rneu = R*sqrt(1-zR*zR)*cos(phiR); rneu++;
00288 *rneu = R*sqrt(1-zR*zR)*sin(phiR); rneu++;
00289 *rneu = R*zR; rneu++;
00290 *rneu = R; rneu++;
00291 *rneu = zR; rneu++;
00292 *rneu = phiR; rneu++;
00293
00294 *rgn = 0; rgn++;
00295 *rgn = R*sqrt(1-zR*zR)*cos(phiR); rgn++;
00296 *rgn = R*sqrt(1-zR*zR)*sin(phiR); rgn++;
00297 *rgn = R*zR; rgn++;
00298 *rgn = R; rgn++;
00299 *rgn = zR; rgn++;
00300 *rgn = phiR; rgn++;
00301 }
00302
00303 // count the neutrons
00304 nneu += Nsources;
00305
00306 } // done with cf source neutron production
00307
00308 //
00309 // transfer temp data to storage
00310 //
00311 Nnubar = nnu;
00312 Npositron = npos;
00313 Nneutron = nneu;
00314
00315 Pnubar = new double[Pdim*nnu];
00316 Ppositron = new double[Pdim*npos];
00317 Pneutron = new double[Pdim*nneu];
00318
00319 Rnubar = new double[Rdim*nnu];
00320 Rpositron = new double[Rdim*npos];
00321 Rgenpos = new double[Rdim*npos];
00322 Rneutron = new double[Rdim*nneu];
00323 Rgenneu = new double[Rdim*nneu];
00324
00325 double* ptr;
00326
00327 ptr = Pnubar+Pdim*Nnubar;
00328 while(ptr>Pnubar)*--ptr=*--pnu;
00329
00330 ptr = Ppositron+Pdim*Npositron;
00331 while(ptr>Ppositron)*--ptr=*--ppos;
00332
00333 ptr = Pneutron+Pdim*Nneutron;
00334 while(ptr>Pneutron)*--ptr=*--pneu;
00335
00336 ptr = Rnubar+Rdim*Nnubar;
00337 while(ptr>Rnubar)*--ptr=*--rnu;
00338
00339 ptr = Rpositron+Rdim*Npositron;
00340 while(ptr>Rpositron)*--ptr=*--rpos;
00341
00342 ptr = Rgenpos+Rdim*Npositron;
00343 while(ptr>Rgenpos)*--ptr=*--rgp;
00344
00345 ptr = Rneutron+Rdim*Nneutron;
00346 while(ptr>Rneutron)*--ptr=*--rneu;
00347
00348 ptr = Rgenneu+Rdim*Nneutron;
00349 while(ptr>Rgenneu)*--ptr=*--rgn;
00350
00351 /*
00352 cout << Nnubar << " nubar(s):" << endl;
00353 for(int i=0;i<Nnubar*Pdim;i++){
00354 cout << " Pnubar[" << i << "] " << Pnubar[i] << endl;
00355 }
00356 for(int i=0;i<Nnubar*Rdim;i++){
00357 cout << " Rnubar[" << i << "] " << Rnubar[i] << endl;
00358 }
00359
00360 cout << Npositron << " positron(s):" << endl;
00361 for(int i=0;i<Npositron*Pdim;i++){
00362 cout << " Ppositron[" << i << "] " << Ppositron[i] << endl;
00363 }
00364 for(int i=0;i<Npositron*Rdim;i++){
00365 cout << " Rpositron[" << i << "] " << Rpositron[i] << endl;
00366 }
00367 for(int i=0;i<Npositron*Rdim;i++){
00368 cout << " Rgenpos[" << i << "] " << Rgenpos[i] << endl;
00369 }
00370
00371 cout << Nneutron << " neutron(s):" << endl;
00372 for(int i=0;i<Nneutron*Pdim;i++){
00373 cout << " Pneutron[" << i << "] " << Pneutron[i] << endl;
00374 }
00375 for(int i=0;i<Nneutron*Rdim;i++){
00376 cout << " Rneutron[" << i << "] " << Rneutron[i] << endl;
00377 }
00378 for(int i=0;i<Nneutron*Rdim;i++){
00379 cout << " Rgenneu[" << i << "] " << Rgenneu[i] << endl;
00380 }
00381 */
00382
00383 }
|
|
|
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 83 of file ReactorEvent.hh.
00083 {return Pnubar[0+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 108 of file ReactorEvent.hh.
00108 {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 115 of file ReactorEvent.hh.
00115 {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 88 of file ReactorEvent.hh.
00088 {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 94 of file ReactorEvent.hh.
00094 {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 101 of file ReactorEvent.hh.
00101 {return Ppositron[6+n*Pdim];}
|
1.2.14 written by Dimitri van Heesch,
© 1997-2002