00001
00011 #include <math.h>
00012 #include <iostream.h>
00013 #include <string.h>
00014 #include "ReactorEvent.hh"
00015 #include "ReactorFortran.hh"
00016 #include "ReactorConstants.hh"
00017 #include "ReactorXC.hh"
00018 #include "ReactorDetector.hh"
00019 #include "MyTlSource.hh"
00020 #include "MyBiSource.hh"
00021
00022
00023
00024 ReactorDetector::ReactorDetector(){
00025
00026
00027
00028 Xpmt = new double[0];
00029 Ypmt = new double[0];
00030 Zpmt = new double[0];
00031 PMTdata = new PMT[0];
00032 PMTelectron = new PMT[0];
00033 PMTpositron = new PMT[0];
00034 PMTneutron = new PMT[0];
00035 PMTthallium = new PMT[0];
00036 PMTbismuth = new PMT[0];
00037 PMTcalifornium = new PMT[0];
00038 PMTmuon = new PMT[0];
00039
00040
00041
00042 SetNominalParameters();
00043 SetPMT();
00044 SetNeutronPropStuff();
00045
00046
00047
00048 Ngamma = 0;
00049 Rgamma = new double[Ngamma];
00050 Egamma = new double[Ngamma];
00051 Tgamma = new double[Ngamma];
00052 Ogamma = new std::string[Ngamma];
00053
00054
00055
00056 float Qlo = 0.001;
00057 float Qhi = 10.0;
00058 funlxp_(QSpectrum,Qspace,Qlo,Qhi);
00059
00060 Nneutron = 0;
00061 Absorber = new double[Nneutron];
00062 NeutronSteps = new double[Nneutron];
00063
00064 pmtRadCount = 0;
00065 }
00066
00067
00068
00069 ReactorDetector::~ReactorDetector(){
00070 delete[] Xpmt;
00071 delete[] Ypmt;
00072 delete[] Zpmt;
00073 delete[] PMTdata;
00074 delete[] PMTelectron;
00075 delete[] PMTpositron;
00076 delete[] PMTneutron;
00077 delete[] PMTthallium;
00078 delete[] PMTbismuth;
00079 delete[] PMTcalifornium;
00080 delete[] PMTmuon;
00081 delete[] Rgamma;
00082 delete[] Egamma;
00083 delete[] Tgamma;
00084 delete[] Ogamma;
00085 delete[] Absorber;
00086 delete[] NeutronSteps;
00087 }
00088
00089
00090
00091 ReactorDetector::ReactorDetector(const ReactorDetector& Detector){
00092
00093 R0=Detector.R0;
00094 R1=Detector.R1;
00095 R2=Detector.R2;
00096 GdConcentration=Detector.GdConcentration;
00097 refracGd=Detector.refracGd;
00098 refracSc=Detector.refracSc;
00099 mfpGd=Detector.mfpGd;
00100 mfpSc=Detector.mfpSc;
00101 MeanNeuDispl = Detector.MeanNeuDispl;
00102 FastNeutronOption = Detector.FastNeutronOption;
00103 Temperature = Detector.Temperature;
00104 tGd=Detector.tGd;
00105 tSc=Detector.tSc;
00106 GdCaptureFraction=Detector.GdCaptureFraction;
00107 GammasPerGd=Detector.GammasPerGd;
00108 attenlGd=Detector.attenlGd;
00109 attenlSc=Detector.attenlSc;
00110 PhotonsPerMeV=Detector.PhotonsPerMeV;
00111 ShortPosFrac0=Detector.ShortPosFrac0;
00112 ShortPosFrac1=Detector.ShortPosFrac1;
00113 LongPosFrac=Detector.LongPosFrac;
00114 ShortPosTime0=Detector.ShortPosTime0;
00115 ShortPosTime1=Detector.ShortPosTime1;
00116 LongPosTime=Detector.LongPosTime;
00117 PMTcoverage=Detector.PMTcoverage;
00118 PMTdiameter=Detector.PMTdiameter;
00119 PMTqe=Detector.PMTqe;
00120 ScintDecayTime=Detector.ScintDecayTime;
00121 for(int i=0;i<4;i++){
00122 a[i] = Detector.a[i];
00123 z[i] = Detector.z[i];
00124 nn[i] = Detector.nn[i];
00125 }
00126
00127 Npmt = Detector.Npmt;
00128 Ncosbins = Detector.Ncosbins;
00129 Nphibins = Detector.Nphibins;
00130
00131 Xpmt = new double[Npmt];
00132 Ypmt = new double[Npmt];
00133 Zpmt = new double[Npmt];
00134
00135 double* pold = Detector.Xpmt+Npmt;
00136 double* pnew = Xpmt+Npmt;
00137 while(pnew>Ypmt)*--pnew=*--pold;
00138
00139 pold = Detector.Ypmt+Npmt;
00140 pnew = Ypmt+Npmt;
00141 while(pnew>Ypmt)*--pnew=*--pold;
00142
00143 pold = Detector.Zpmt+Npmt;
00144 pnew = Zpmt+Npmt;
00145 while(pnew>Zpmt)*--pnew=*--pold;
00146
00147 PHelectron = Detector.PHelectron;
00148 PHpositron = Detector.PHpositron;
00149 PHneutron = Detector.PHneutron;
00150 PHthallium = Detector.PHthallium;
00151 PHbismuth = Detector.PHbismuth;
00152 PHcalifornium = Detector.PHcalifornium;
00153 PHmuon = Detector.PHmuon;
00154
00155 PHeleReco = Detector.PHeleReco;
00156 PHposReco = Detector.PHposReco;
00157 PHneuReco = Detector.PHneuReco;
00158 PHthaReco = Detector.PHthaReco;
00159 PHbisReco = Detector.PHbisReco;
00160 PHcalReco = Detector.PHcalReco;
00161 PHmuoReco = Detector.PHmuoReco;
00162
00163 TimePositron = Detector.TimePositron;
00164 TimeNeutron = Detector.TimeNeutron;
00165
00166 for(int i=0;i<3;i++){
00167 DipolePositron[i] = Detector.DipolePositron[i];
00168 DipoleNeutron[i] = Detector.DipoleNeutron[i];
00169 }
00170
00171 HitPmtsPos = Detector.HitPmtsPos;
00172 HitPmtsNeu = Detector.HitPmtsNeu;
00173 HitPmtsTha = Detector.HitPmtsTha;
00174 HitPmtsBis = Detector.HitPmtsBis;
00175 HitPmtsCal = Detector.HitPmtsCal;
00176 HitPmtsMuo = Detector.HitPmtsMuo;
00177 HitPmtsEle = Detector.HitPmtsEle;
00178
00179 for(int i=0;i<7;i++){
00180 MeanNpe[i]=Detector.MeanNpe[i];
00181 MeanAtten[i] = Detector.MeanAtten[i];
00182 MeanSfactor[i] = Detector.MeanSfactor[i];
00183 }
00184
00185 Nneutron = Detector.Nneutron;
00186 Absorber = new double[Nneutron];
00187 pold = Detector.Absorber+Nneutron;
00188 pnew = Absorber+Nneutron;
00189 while(pnew>Absorber)*--pnew=*--pold;
00190
00191 NeutronSteps = new double[Nneutron];
00192 pold = Detector.NeutronSteps+Nneutron;
00193 pnew = NeutronSteps+Nneutron;
00194 while(pnew>NeutronSteps)*--pnew=*--pold;
00195
00196 Ngamma = Detector.Ngamma;
00197 Rgamma = new double[Ngamma];
00198 pold = Detector.Rgamma+3*Ngamma;
00199 pnew = Rgamma+3*Ngamma;
00200 while(pnew>Rgamma)*--pnew=*--pold;
00201
00202 Egamma = new double[Ngamma];
00203 pold = Detector.Egamma+Ngamma;
00204 pnew = Egamma+Ngamma;
00205 while(pnew>Egamma)*--pnew=*--pold;
00206
00207 Tgamma = new double[Ngamma];
00208 pold = Detector.Tgamma+Ngamma;
00209 pnew = Tgamma+Ngamma;
00210 while(pnew>Tgamma)*--pnew=*--pold;
00211
00212 Ogamma = new std::string[Ngamma];
00213 std::string* rold = Detector.Ogamma+Ngamma;
00214 std::string* rnew = Ogamma+Ngamma;
00215 while(rnew>Ogamma)*--rnew=*--rold;
00216
00217 PMTsize = Detector.PMTsize;
00218 PMTlen = Detector.PMTlen;
00219 PMTdata = new PMT[PMTsize];
00220 PMT* qold = Detector.PMTdata+PMTsize;
00221 PMT* qnew = PMTdata+PMTsize;
00222 while(qnew>PMTdata)*--qnew=*--qold;
00223
00224 PMTelectron = new PMT[Npmt];
00225 qold = Detector.PMTelectron+Npmt;
00226 qnew = PMTelectron+Npmt;
00227 while(qnew>PMTelectron)*--qnew=*--qold;
00228
00229 PMTpositron = new PMT[Npmt];
00230 qold = Detector.PMTpositron+Npmt;
00231 qnew = PMTpositron+Npmt;
00232 while(qnew>PMTpositron)*--qnew=*--qold;
00233
00234 PMTneutron = new PMT[Npmt];
00235 qold = Detector.PMTneutron+Npmt;
00236 qnew = PMTneutron+Npmt;
00237 while(qnew>PMTneutron)*--qnew=*--qold;
00238
00239 PMTthallium = new PMT[Npmt];
00240 qold = Detector.PMTthallium+Npmt;
00241 qnew = PMTthallium+Npmt;
00242 while(qnew>PMTthallium)*--qnew=*--qold;
00243
00244 PMTbismuth = new PMT[Npmt];
00245 qold = Detector.PMTbismuth+Npmt;
00246 qnew = PMTbismuth+Npmt;
00247 while(qnew>PMTbismuth)*--qnew=*--qold;
00248
00249 PMTcalifornium = new PMT[Npmt];
00250 qold = Detector.PMTcalifornium+Npmt;
00251 qnew = PMTcalifornium+Npmt;
00252 while(qnew>PMTcalifornium)*--qnew=*--qold;
00253
00254 PMTmuon = new PMT[Npmt];
00255 qold = Detector.PMTmuon+Npmt;
00256 qnew = PMTmuon+Npmt;
00257 while(qnew>PMTmuon)*--qnew=*--qold;
00258 }
00259
00260
00261
00262 ReactorDetector& ReactorDetector::operator=(const ReactorDetector& rhs){
00263
00264 if(this == &rhs)return *this;
00265
00266 R0=rhs.R0;
00267 R1=rhs.R1;
00268 R2=rhs.R2;
00269 GdConcentration=rhs.GdConcentration;
00270 refracGd=rhs.refracGd;
00271 refracSc=rhs.refracSc;
00272 mfpGd=rhs.mfpGd;
00273 mfpSc=rhs.mfpSc;
00274 MeanNeuDispl = rhs.MeanNeuDispl;
00275 FastNeutronOption = rhs.FastNeutronOption;
00276 Temperature = rhs.Temperature;
00277 tGd=rhs.tGd;
00278 tSc=rhs.tSc;
00279 GdCaptureFraction=rhs.GdCaptureFraction;
00280 GammasPerGd=rhs.GammasPerGd;
00281 attenlGd=rhs.attenlGd;
00282 attenlSc=rhs.attenlSc;
00283 PhotonsPerMeV=rhs.PhotonsPerMeV;
00284 ShortPosFrac0=rhs.ShortPosFrac0;
00285 ShortPosFrac1=rhs.ShortPosFrac1;
00286 LongPosFrac=rhs.LongPosFrac;
00287 ShortPosTime0=rhs.ShortPosTime0;
00288 ShortPosTime1=rhs.ShortPosTime1;
00289 LongPosTime=rhs.LongPosTime;
00290 PMTcoverage=rhs.PMTcoverage;
00291 PMTdiameter=rhs.PMTdiameter;
00292 PMTqe=rhs.PMTqe;
00293 ScintDecayTime=rhs.ScintDecayTime;
00294 for(int i=0;i<4;i++){
00295 a[i] = rhs.a[i];
00296 z[i] = rhs.z[i];
00297 nn[i] = rhs.nn[i];
00298 }
00299
00300 Npmt = rhs.Npmt;
00301 Ncosbins = rhs.Ncosbins;
00302 Nphibins = rhs.Nphibins;
00303
00304 delete[] Xpmt;
00305 delete[] Ypmt;
00306 delete[] Zpmt;
00307
00308 Xpmt = new double[Npmt];
00309 Ypmt = new double[Npmt];
00310 Zpmt = new double[Npmt];
00311
00312 double* pold = rhs.Xpmt+Npmt;
00313 double* pnew = Xpmt+Npmt;
00314 while(pnew>Ypmt)*--pnew=*--pold;
00315
00316 pold = rhs.Ypmt+Npmt;
00317 pnew = Ypmt+Npmt;
00318 while(pnew>Ypmt)*--pnew=*--pold;
00319
00320 pold = rhs.Zpmt+Npmt;
00321 pnew = Zpmt+Npmt;
00322 while(pnew>Zpmt)*--pnew=*--pold;
00323
00324 PHelectron = rhs.PHelectron;
00325 PHpositron = rhs.PHpositron;
00326 PHneutron = rhs.PHneutron;
00327 PHthallium = rhs.PHthallium;
00328 PHbismuth = rhs.PHbismuth;
00329 PHcalifornium = rhs.PHcalifornium;
00330 PHmuon = rhs.PHmuon;
00331
00332 PHeleReco = rhs.PHeleReco;
00333 PHposReco = rhs.PHposReco;
00334 PHneuReco = rhs.PHneuReco;
00335 PHthaReco = rhs.PHthaReco;
00336 PHbisReco = rhs.PHbisReco;
00337 PHcalReco = rhs.PHcalReco;
00338 PHmuoReco = rhs.PHmuoReco;
00339
00340 TimePositron = rhs.TimePositron;
00341 TimeNeutron = rhs.TimeNeutron;
00342
00343 for(int i=0;i<3;i++){
00344 DipolePositron[i] = rhs.DipolePositron[i];
00345 DipoleNeutron[i] = rhs.DipoleNeutron[i];
00346 }
00347
00348 HitPmtsPos = rhs.HitPmtsPos;
00349 HitPmtsNeu = rhs.HitPmtsNeu;
00350 HitPmtsTha = rhs.HitPmtsTha;
00351 HitPmtsBis = rhs.HitPmtsBis;
00352 HitPmtsCal = rhs.HitPmtsCal;
00353 HitPmtsMuo = rhs.HitPmtsMuo;
00354 HitPmtsEle = rhs.HitPmtsEle;
00355
00356 for (int i=0;i<7;i++){
00357 MeanNpe[i] = rhs.MeanNpe[i];
00358 MeanAtten[i] = rhs.MeanAtten[i];
00359 MeanSfactor[i] = rhs.MeanSfactor[i];
00360 }
00361
00362 Nneutron = rhs.Nneutron;
00363
00364 delete[] Absorber;
00365 Absorber = new double[Nneutron];
00366 pold = rhs.Absorber+Nneutron;
00367 pnew = Absorber+Nneutron;
00368 while(pnew>Absorber)*--pnew=*--pold;
00369
00370 delete[] NeutronSteps;
00371 NeutronSteps = new double[Nneutron];
00372 pold = rhs.NeutronSteps+Nneutron;
00373 pnew = NeutronSteps+Nneutron;
00374 while(pnew>NeutronSteps)*--pnew=*--pold;
00375
00376 Ngamma = rhs.Ngamma;
00377
00378 delete[] Rgamma;
00379 Rgamma = new double[Ngamma];
00380 pold = rhs.Rgamma+3*Ngamma;
00381 pnew = Rgamma+3*Ngamma;
00382 while(pnew>Rgamma)*--pnew=*--pold;
00383
00384 delete[] Egamma;
00385 Egamma = new double[Ngamma];
00386 pold = rhs.Egamma+Ngamma;
00387 pnew = Egamma+Ngamma;
00388 while(pnew>Egamma)*--pnew=*--pold;
00389
00390 delete[] Tgamma;
00391 Tgamma = new double[Ngamma];
00392 pold = rhs.Tgamma+Ngamma;
00393 pnew = Tgamma+Ngamma;
00394 while(pnew>Tgamma)*--pnew=*--pold;
00395
00396 delete[] Ogamma;
00397 Ogamma = new std::string[Ngamma];
00398 std::string* rold = rhs.Ogamma+Ngamma;
00399 std::string* rnew = Ogamma+Ngamma;
00400 while(rnew>Ogamma)*--rnew=*--rold;
00401
00402 PMTsize = rhs.PMTsize;
00403 PMTlen = rhs.PMTlen;
00404 delete[] PMTdata;
00405 PMTdata = new PMT[PMTsize];
00406 PMT* qold = rhs.PMTdata+PMTsize;
00407 PMT* qnew = PMTdata+PMTsize;
00408 while(qnew>PMTdata)*--qnew=*--qold;
00409
00410 delete[] PMTelectron;
00411 PMTelectron = new PMT[Npmt];
00412 qold = rhs.PMTelectron+Npmt;
00413 qnew = PMTelectron+Npmt;
00414 while(qnew>PMTelectron)*--qnew=*--qold;
00415
00416 delete[] PMTpositron;
00417 PMTpositron = new PMT[Npmt];
00418 qold = rhs.PMTpositron+Npmt;
00419 qnew = PMTpositron+Npmt;
00420 while(qnew>PMTpositron)*--qnew=*--qold;
00421
00422 delete[] PMTneutron;
00423 PMTneutron = new PMT[Npmt];
00424 qold = rhs.PMTneutron+Npmt;
00425 qnew = PMTneutron+Npmt;
00426 while(qnew>PMTneutron)*--qnew=*--qold;
00427
00428 delete[] PMTthallium;
00429 PMTthallium = new PMT[Npmt];
00430 qold = rhs.PMTthallium+Npmt;
00431 qnew = PMTthallium+Npmt;
00432 while(qnew>PMTthallium)*--qnew=*--qold;
00433
00434 delete[] PMTbismuth;
00435 PMTbismuth = new PMT[Npmt];
00436 qold = rhs.PMTbismuth+Npmt;
00437 qnew = PMTbismuth+Npmt;
00438 while(qnew>PMTbismuth)*--qnew=*--qold;
00439
00440 delete[] PMTcalifornium;
00441 PMTcalifornium = new PMT[Npmt];
00442 qold = rhs.PMTcalifornium+Npmt;
00443 qnew = PMTcalifornium+Npmt;
00444 while(qnew>PMTcalifornium)*--qnew=*--qold;
00445
00446 delete[] PMTmuon;
00447 PMTmuon = new PMT[Npmt];
00448 qold = rhs.PMTmuon+Npmt;
00449 qnew = PMTmuon+Npmt;
00450 while(qnew>PMTmuon)*--qnew=*--qold;
00451
00452 return *this;
00453
00454 }
00455
00456
00457
00458 ReactorDetector::PMT::PMT(){
00459 npmt=0;
00460 for (int i=0;i<3;i++)xyz[i]=0;
00461 TimeGen=0;
00462 TimeSmr=0;
00463 PHgen=0;
00464 PHsmr=0;
00465 PathGd=0;
00466 PathSc=0;
00467 Sfactor=0;
00468 Atten=0;
00469 Npe=0;
00470 Xpe=0;
00471 Q=0.;
00472 Parent="U";
00473
00474 SubHits=0;
00475 SubHitsPtr = new PMT*[0];
00476 }
00477
00478 ReactorDetector::PMT::~PMT(){
00479 delete[] SubHitsPtr;
00480 }
00481
00482 ReactorDetector::PMT::PMT(const ReactorDetector::PMT& p){
00483
00484 npmt=p.npmt;
00485 for (int i=0;i<3;i++)xyz[i]=p.xyz[i];
00486 TimeGen=p.TimeGen;
00487 TimeSmr=p.TimeSmr;
00488 PHgen=p.PHgen;
00489 PHsmr=p.PHsmr;
00490 PathGd=p.PathGd;
00491 PathSc=p.PathSc;
00492 Sfactor=p.Sfactor;
00493 Atten=p.Atten;
00494 Npe=p.Npe;
00495 Xpe=p.Xpe;
00496 Q=p.Q;
00497 Parent=p.Parent;
00498
00499 SubHits=p.SubHits;
00500 SubHitsPtr = new PMT*[SubHits];
00501 PMT** pold = p.SubHitsPtr+SubHits;
00502 PMT** pnew = SubHitsPtr+SubHits;
00503 while(pnew>SubHitsPtr)*--pnew=*--pold;
00504 }
00505
00506
00507
00508 ReactorDetector::PMT&
00509 ReactorDetector::PMT::operator=(const ReactorDetector::PMT& rhs){
00510
00511 if(this == &rhs)return *this;
00512
00513 npmt=rhs.npmt;
00514 for (int i=0;i<3;i++)xyz[i]=rhs.xyz[i];
00515 TimeGen=rhs.TimeGen;
00516 TimeSmr=rhs.TimeSmr;
00517 PHgen=rhs.PHgen;
00518 PHsmr=rhs.PHsmr;
00519 PathGd=rhs.PathGd;
00520 PathSc=rhs.PathSc;
00521 Sfactor=rhs.Sfactor;
00522 Atten=rhs.Atten;
00523 Npe=rhs.Npe;
00524 Xpe=rhs.Xpe;
00525 Q=rhs.Q;
00526 Parent=rhs.Parent;
00527
00528 SubHits=rhs.SubHits;
00529 delete[] SubHitsPtr;
00530 SubHitsPtr = new PMT*[SubHits];
00531 PMT** pold = rhs.SubHitsPtr+SubHits;
00532 PMT** pnew = SubHitsPtr+SubHits;
00533 while(pnew>SubHitsPtr)*--pnew=*--pold;
00534
00535 return *this;
00536 }
00537 ReactorDetector::PMT
00538 ReactorDetector::PMT::operator+(const ReactorDetector::PMT& rhs){
00539
00540 PMT sum(*this);
00541 if(npmt!=rhs.npmt)return sum;
00542 sum.npmt = npmt;
00543
00544 double w1=Xpe/(Xpe+rhs.Xpe);
00545 double w2=1-w1;
00546
00547 for (int i=0;i<3;i++)sum.xyz[i]=w1*xyz[i]+w2*rhs.xyz[i];
00548 sum.TimeGen=w1*TimeGen+w2*rhs.TimeGen;
00549 sum.TimeSmr=w1*TimeSmr+w2*rhs.TimeSmr;
00550 sum.PHgen=w1*PHgen+w2*rhs.PHgen;
00551 sum.PHsmr=w1*PHsmr+w2*rhs.PHsmr;
00552 sum.PathGd=w1*PathGd+w2*rhs.PathGd;
00553 sum.PathSc=w1*PathSc+w2*rhs.PathSc;
00554 sum.Sfactor=w1*Sfactor+w1*rhs.Sfactor;
00555 sum.Atten=w1*Atten+w2*rhs.Atten;
00556 sum.Npe=Npe+rhs.Npe;
00557 sum.Xpe=Xpe+rhs.Xpe;
00558 if(Parent==rhs.Parent)sum.Parent=rhs.Parent;
00559 else sum.Parent="U";
00560
00561 sum.SubHits = SubHits+rhs.SubHits;
00562 delete[] sum.SubHitsPtr;
00563 sum.SubHitsPtr = new PMT*[sum.SubHits];
00564 PMT** p = sum.SubHitsPtr+sum.SubHits;
00565 PMT** q = SubHitsPtr+SubHits;
00566 PMT** r = rhs.SubHitsPtr+rhs.SubHits;
00567 while(p>sum.SubHitsPtr+SubHits)*--p=*--r;
00568 while(p>sum.SubHitsPtr)*--p=*--q;
00569
00570 return sum;
00571 }
00572 ReactorDetector::PMT&
00573 ReactorDetector::PMT::operator+=(const ReactorDetector::PMT& rhs){
00574
00575 if(npmt!=rhs.npmt)return *this;
00576
00577 double w1=Xpe/(Xpe+rhs.Xpe);
00578 double w2=1-w1;
00579
00580 for (int i=0;i<3;i++)xyz[i]=w1*xyz[i]+w2*rhs.xyz[i];
00581 TimeGen=w1*TimeGen+w2*rhs.TimeGen;
00582 TimeSmr=w1*TimeSmr+w2*rhs.TimeSmr;
00583 PHgen=w1*PHgen+w2*rhs.PHgen;
00584 PHsmr=w1*PHsmr+w2*rhs.PHsmr;
00585 PathGd=w1*PathGd+w2*rhs.PathGd;
00586 PathSc=w1*PathSc+w2*rhs.PathSc;
00587 Sfactor=w1*Sfactor+w1*rhs.Sfactor;
00588 Atten=w1*Atten+w2*rhs.Atten;
00589 Npe+=rhs.Npe;
00590 Xpe+=rhs.Xpe;
00591 Q+=rhs.Q;
00592 if(Parent!=rhs.Parent)Parent="U";
00593
00594 PMT** SubHitsTmp = new PMT*[SubHits];
00595 PMT** p = SubHitsTmp+SubHits;
00596 PMT** q = SubHitsPtr+SubHits;
00597 while(p>SubHitsTmp)*--p=*--q;
00598 delete[] SubHitsPtr;
00599 SubHitsPtr = new PMT*[SubHits+rhs.SubHits];
00600 p = SubHitsPtr+SubHits+rhs.SubHits;
00601 q = SubHitsTmp+SubHits;
00602 PMT** r = rhs.SubHitsPtr+rhs.SubHits;
00603 while(p>SubHitsPtr+SubHits)*--p=*--r;
00604 while(p>SubHitsPtr)*--p=*--q;
00605 SubHits += rhs.SubHits;
00606 delete[] SubHitsTmp;
00607
00608 return *this;
00609 }
00610 void ReactorDetector::PMT::SetSubHitsPtr(int subhits){
00611 delete[] SubHitsPtr;
00612 SubHits = subhits;
00613 SubHitsPtr = new PMT*[subhits];
00614 return;
00615 }
00616
00617
00618
00619
00620
00621 void ReactorDetector::LightsOut(ReactorEvent& Event){
00622
00623 PMTlen=0;
00624 delete[] PMTdata;
00625 PMTdata = new PMT[PMTsize];
00626
00627
00628
00629 GenerateSecondaries(Event);
00630 GenerateGammas(Event);
00631 GenerateResponse(Event);
00632
00633 ConsolidatePMT();
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643 CalculateQPosition("P");
00644 CalculateQPosition("N");
00645
00646 ImproveQPosition("P");
00647 ImproveQPosition("N");
00648
00649 CalculateEnergy("P");
00650 CalculateEnergy("N");
00651
00652 }
00653
00654
00655
00656 void ReactorDetector::GetPMTdata(int npmt,double* XYZpmt,
00657 int& nele,double* PHele,double* Tele,
00658 int& npos,double* PHpos,double* Tpos,
00659 int& nneu,double* PHneu,double* Tneu,
00660 int& ntha,double* PHtha,double* Ttha,
00661 int& nbis,double* PHbis,double* Tbis,
00662 int& ncal,double* PHcal,double* Tcal,
00663 int& nmuo,double* PHmuo,double* Tmuo,
00664 int* npepos, int* npeneu, int* npetha,
00665 int* npebis, int* npecal, int* npemuo,
00666 int* npeele){
00667
00668 *(XYZpmt+0) = *(Xpmt+npmt);
00669 *(XYZpmt+1) = *(Ypmt+npmt);
00670 *(XYZpmt+2) = *(Zpmt+npmt);
00671
00672 nele=0;
00673 npos=0;
00674 nneu=0;
00675 ntha=0;
00676 nbis=0;
00677 ncal=0;
00678 nmuo=0;
00679
00680 for(PMT* PMTptr=PMTdata;PMTptr<PMTdata+PMTlen;PMTptr++){
00681 if(npmt==(*PMTptr).npmt){
00682 if((*PMTptr).Parent=="E"){
00683 *(PHele+nele) = (*PMTptr).PHsmr;
00684 *(Tele+nele) = (*PMTptr).TimeSmr;
00685 *(npeele+nele) = (*PMTptr).Npe;
00686 nele++;
00687 }
00688 if((*PMTptr).Parent=="P"){
00689 *(PHpos+npos) = (*PMTptr).PHsmr;
00690 *(Tpos+npos) = (*PMTptr).TimeSmr;
00691 *(npepos+npos) = (*PMTptr).Npe;
00692 npos++;
00693 }
00694 if((*PMTptr).Parent=="N"){
00695 *(PHneu+nneu) = (*PMTptr).PHsmr;
00696 *(Tneu+nneu) = (*PMTptr).TimeSmr;
00697 *(npeneu+nneu) = (*PMTptr).Npe;
00698 nneu++;
00699 }
00700 if((*PMTptr).Parent=="Tl"){
00701 *(PHtha+ntha) = (*PMTptr).PHsmr;
00702 *(Ttha+ntha) = (*PMTptr).TimeSmr;
00703 *(npetha+ntha) = (*PMTptr).Npe;
00704 ntha++;
00705 }
00706 if((*PMTptr).Parent=="Bi"){
00707 *(PHbis+nbis) = (*PMTptr).PHsmr;
00708 *(Tbis+nbis) = (*PMTptr).TimeSmr;
00709 *(npebis+nbis) = (*PMTptr).Npe;
00710 nbis++;
00711 }
00712 if((*PMTptr).Parent=="Cf"){
00713 *(PHcal+ncal) = (*PMTptr).PHsmr;
00714 *(Tcal+ncal) = (*PMTptr).TimeSmr;
00715 *(npecal+ncal) = (*PMTptr).Npe;
00716 ncal++;
00717 }
00718 if((*PMTptr).Parent=="muon"){
00719 *(PHmuo+nmuo) = (*PMTptr).PHsmr;
00720 *(Tmuo+nmuo) = (*PMTptr).TimeSmr;
00721 *(npemuo+nmuo) = (*PMTptr).Npe;
00722 nmuo++;
00723 }
00724 }
00725 }
00726 return;
00727 }
00728 int ReactorDetector::GetHitPMTsPos(){
00729 return HitPmtsPos;
00730 }
00731 int ReactorDetector::GetHitPMTsNeu(){
00732 return HitPmtsNeu;
00733 }
00734 int ReactorDetector::GetHitPMTsTha(){
00735 return HitPmtsTha;
00736 }
00737 int ReactorDetector::GetHitPMTsBis(){
00738 return HitPmtsBis;
00739 }
00740 int ReactorDetector::GetHitPMTsCal(){
00741 return HitPmtsCal;
00742 }
00743 int ReactorDetector::GetHitPMTsMuo(){
00744 return HitPmtsMuo;
00745 }
00746 int ReactorDetector::GetHitPMTsEle(){
00747 return HitPmtsEle;
00748 }
00749 void ReactorDetector::GenerateSecondaries(ReactorEvent& Event){
00750
00751
00752
00753
00754
00755 int Npositron = Event.GetNumPositron();
00756
00757
00758 double Rpos[Event.Rdim*Npositron];
00759
00760
00761 for(int n=0;n<Npositron;n++){
00762
00763 double kepos = Event.GetPositronKE(n);
00764 double tpos = 0.;
00765 double rangepos = 0.;
00766 PositronRange(kepos,rangepos,tpos);
00767 tpos+=Event.GetPositronTime(n);
00768
00769 double cospos = Event.GetPositronCosTheta(n);
00770 double phipos = Event.GetPositronPhi(n);
00771
00772 double rint[3];
00773 for(int i=0;i<3;i++){
00774 rint[i] = *(Event.GetPositronVertexXYZ()+(i+n*Event.Rdim));
00775 }
00776
00777 double rpos[3];
00778 rpos[0]=rint[0]+rangepos*sqrt(1-cospos*cospos)*cos(phipos);
00779 rpos[1]=rint[1]+rangepos*sqrt(1-cospos*cospos)*sin(phipos);
00780 rpos[2]=rint[2]+rangepos*cospos;
00781 double rrpos = sqrt(rpos[0]*rpos[0]+rpos[1]*rpos[1]+rpos[2]*rpos[2]);
00782
00783
00784
00785 Rpos[0+n*Event.Rdim] = tpos;
00786 for(int i=1;i<4;i++)Rpos[i+n*Event.Rdim]=rpos[i-1];
00787 Rpos[4+n*Event.Rdim] = rrpos;
00788 Rpos[5+n*Event.Rdim] = rpos[2]/rrpos;
00789 Rpos[6+n*Event.Rdim] = atan2(rpos[1],rpos[0]);
00790 if(Rpos[6+n*Event.Rdim]<0)Rpos[6+n*Event.Rdim]+=2*RC.pi;
00791 Event.SetPositronVertex(n,Rpos);
00792 }
00793
00794
00795
00796
00797
00798 int Nneutron = Event.GetNumNeutron();
00799
00800
00801 delete[] Absorber;
00802 delete[] NeutronSteps;
00803
00804 Absorber = new double[Nneutron];
00805 NeutronSteps = new double[Nneutron];
00806
00807 double Rneu[Event.Rdim*Nneutron];
00808
00809
00810 for(int n=0;n<Nneutron;n++){
00811
00812 double rneu[3];
00813 double tneu = 0.;
00814 if(FastNeutronOption){
00815
00816
00817
00818 float rn[3];
00819 int len=3;
00820 rnormx_(rn,len,ranlux_);
00821
00822
00823
00824 len=1;
00825 float rv[1];
00826 ranlux_(rv,len);
00827
00828
00829
00830
00831
00832 double rint[3];
00833 for(int i=0;i<3;i++){
00834 rint[i] = *(Event.GetNeutronVertexXYZ()+(i+n*Event.Rdim));
00835 }
00836
00837 double rrint = sqrt(rint[0]*rint[0]+rint[1]*rint[1]+rint[2]*rint[2]);
00838 double path = 0.;
00839 if(rrint<=R0){
00840 for (int i=0;i<3;i++)rneu[i]=rint[i]+mfpGd*rn[i]/sqrt(3.);
00841 path = mfpGd*sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
00842 tneu = tGd*sqrt(path/mfpGd)*log(1/rv[0]);
00843 }
00844 else if(rrint<R2){
00845 for (int i=0;i<3;i++)rneu[i]=rint[i]+mfpSc*rn[i]/sqrt(3.);
00846 path = mfpSc*sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
00847 tneu = tSc*sqrt(path/mfpSc)*log(1/rv[0]);
00848 }
00849 else{
00850 for (int i=0;i<3;i++)rneu[i]=999999;
00851 tneu = 999999.;
00852 }
00853
00854
00855
00856 rneu[2]+=MeanNeuDispl;
00857 }
00858 else{
00859 double pneu[3];
00860 for(int i=0;i<3;i++){
00861 pneu[i]=*(Event.GetNeutronPxyz()+(i+n*Event.Pdim));
00862
00863 }
00864 for(int i=0;i<3;i++){
00865 rneu[i]=*(Event.GetNeutronVertexXYZ()+(i+n*Event.Rdim));
00866
00867 }
00868 double keneu = Event.GetNeutronKE(n);
00869 NeutronSteps[n]=0.;
00870 Absorber[n]=0.;
00871 NeutronPropagator(R0,R2,keneu,Absorber[n],NeutronSteps[n],
00872 pneu,rneu,tneu,Temperature,a,z,nn);
00873
00874
00875 tneu+=Event.GetNeutronTime(n);
00876
00877 }
00878
00879
00880
00881 Rneu[0+n*Event.Rdim] = tneu;
00882 for(int i=1;i<4;i++)Rneu[i+n*Event.Rdim]=rneu[i-1];
00883 double rrneu = sqrt(rneu[0]*rneu[0]+rneu[1]*rneu[1]+rneu[2]*rneu[2]);
00884 Rneu[4+n*Event.Rdim] = rrneu;
00885 Rneu[5+n*Event.Rdim] = rneu[2]/rrneu;
00886 Rneu[6+n*Event.Rdim] = atan2(rneu[1],rneu[0]);
00887 if(Rneu[6+n*Event.Rdim]<0)Rneu[6+n*Event.Rdim]+=2*RC.pi;
00888 Event.SetNeutronVertex(n,Rneu);
00889 }
00890
00891
00892
00893
00894 }
00895
00896 void ReactorDetector::GenerateResponse(ReactorEvent& Event){
00897
00898
00899
00900 PHelectron = 0;
00901 PHpositron = 0;
00902 PHneutron = 0;
00903 PHthallium = 0;
00904 PHbismuth = 0;
00905 PHcalifornium = 0;
00906 PHmuon = 0;
00907
00908 PHeleReco = 0;
00909 PHposReco = 0;
00910 PHneuReco = 0;
00911 PHthaReco = 0;
00912 PHbisReco = 0;
00913 PHcalReco = 0;
00914 PHmuoReco = 0;
00915
00916 HitPmtsPos = 0;
00917 HitPmtsNeu = 0;
00918 HitPmtsTha = 0;
00919 HitPmtsBis = 0;
00920 HitPmtsCal = 0;
00921 HitPmtsMuo = 0;
00922 HitPmtsEle = 0;
00923
00924 for(int i=0;i<7;i++){
00925 MeanSfactor[i] = 0;
00926 MeanAtten[i] = 0;
00927 MeanNpe[i] = 0;
00928 MeanQ[i] = 0;
00929 }
00930
00931
00932
00933 int pmtindex=0;
00934 double* xpmt=Xpmt;
00935 double* ypmt=Ypmt;
00936 double* zpmt=Zpmt;
00937
00938 PMT* PMTptr = PMTdata;
00939
00940
00941 for(int i=0;i<Ncosbins;i++){
00942 for(int j=0;j<Nphibins;j++){
00943 double rpmt[3];
00944 rpmt[0] = *xpmt++;
00945 rpmt[1] = *ypmt++;
00946 rpmt[2] = *zpmt++;
00947
00948 double r[3];
00949 double dr[3];
00950 double tr[3];
00951
00952 double sumele[7]={0,0,0,0,0,0,0};
00953 double sumpos[7]={0,0,0,0,0,0,0};
00954 double sumneu[7]={0,0,0,0,0,0,0};
00955 double sumtha[7]={0,0,0,0,0,0,0};
00956 double sumbis[7]={0,0,0,0,0,0,0};
00957 double sumcal[7]={0,0,0,0,0,0,0};
00958 double summuo[7]={0,0,0,0,0,0,0};
00959
00960 double* rgptr=Rgamma;
00961 double* egptr=Egamma;
00962 double* tgptr=Tgamma;
00963 std::string* ogptr=Ogamma;
00964
00965 int nelegt0=0;
00966 int nposgt0=0;
00967 int nneugt0=0;
00968 int nthagt0=0;
00969 int nbisgt0=0;
00970 int ncalgt0=0;
00971 int nmuogt0=0;
00972
00973 for(int k=0;k<Ngamma;k++){
00974
00975 for (int l=0;l<3;l++)r[l]=*rgptr++;
00976 double rdotr = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
00977
00978
00979
00980
00981 if(sqrt(rdotr)>R1){
00982 double dummy = *egptr++;
00983 dummy = *tgptr++;
00984 std::string sdummy = *ogptr++;
00985 continue;
00986 }
00987
00988 for (int l=0;l<3;l++)dr[l]=rpmt[l]-r[l];
00989 double ddr = sqrt(dr[0]*dr[0]+dr[1]*dr[1]+dr[2]*dr[2]);
00990
00991 for (int l=0;l<3;l++)tr[l]=dr[l]/ddr;
00992 double tdotr = r[0]*tr[0]+r[1]*tr[1]+r[2]*tr[2];
00993
00994
00995
00996
00997 double sGd;
00998 double sSc;
00999
01000
01001
01002 double test = R0*R0+tdotr*tdotr-rdotr;
01003 bool inter = (test>=0);
01004
01005
01006
01007
01008
01009
01010 if(inter){
01011 double sint1 = -tdotr+sqrt(test);
01012 double sint2 = -tdotr-sqrt(test);
01013 if(sint1>0&&sint2<0){
01014 sGd = sint1;
01015 }
01016 else if(sint2>0&&sint1<0){
01017 sGd = sint2;
01018 }
01019 else if(sint2>0&&sint1>0){
01020 sGd = sint1-sint2;
01021 if(sGd<0)sGd*=-1;
01022 }
01023 else{
01024 sGd = 0;
01025 }
01026 sSc = ddr-sGd;
01027 }
01028
01029
01030
01031 else{
01032 sGd = 0;
01033 sSc = ddr;
01034 }
01035
01036
01037
01038
01039 double dsGd = sGd-R0;
01040 double dsSc = sSc-R2+R0;
01041 double atten0 = exp(-sGd/attenlGd-sSc/attenlSc);
01042 double atten = exp(-dsGd/attenlGd-dsSc/attenlSc);
01043
01044
01045
01046 double sfactor=0;
01047 for (int l=0;l<3;l++){sfactor+=tr[l]*rpmt[l]/R2;}
01048 sfactor *= R2*R2/ddr/ddr;
01049
01050
01051
01052 double phgen = *egptr++;
01053
01054
01055
01056
01057
01058
01059 int len = 1;
01060 float photons[1];
01061 rnormx_(photons,len,ranlux_);
01062 photons[0]*=sqrt(phgen*PhotonsPerMeV);
01063 float xpe = (phgen*PhotonsPerMeV) + (double) photons[0];
01064 xpe *= (atten0*sfactor);
01065 xpe *= (PMTcoverage/Npmt*PMTqe);
01066
01067 int npe;
01068 int error;
01069 rnpssn_(xpe,npe,error);
01070
01071
01072 float Q = 0.;
01073 Q = GenerateQ(npe);
01074 double phsmr = Q/PMTqe/PMTcoverage/PhotonsPerMeV;
01075
01076
01077 double tgen = *tgptr++;
01078
01079 double tsmr = refracGd*sGd/RC.c+refracSc*sSc/RC.c + tgen;
01080
01081
01082 std::string type = *ogptr++;
01083
01084 if(type == "E"){
01085 nelegt0++;
01086 sumele[0]+=phsmr;
01087 sumele[1]+=phgen;
01088 sumele[2]+=phgen*atten0*sfactor;
01089 sumele[3]+=npe;
01090 sumele[4]+=sfactor;
01091 sumele[5]+=atten;
01092 sumele[6]+=Q;
01093
01094 PHelectron+=phsmr;
01095 }
01096 else if(type == "P"){
01097 nposgt0++;
01098 sumpos[0]+=phsmr;
01099 sumpos[1]+=phgen;
01100 sumpos[2]+=phgen*atten0*sfactor;
01101 sumpos[3]+=npe;
01102 sumpos[4]+=sfactor;
01103 sumpos[5]+=atten;
01104 sumpos[6]+=Q;
01105
01106 PHpositron+=phsmr;
01107 }
01108 else if(type == "N"){
01109 nneugt0++;
01110 sumneu[0]+=phsmr;
01111 sumneu[1]+=phgen;
01112 sumneu[2]+=phgen*atten0*sfactor;
01113 sumneu[3]+=npe;
01114 sumneu[4]+=sfactor;
01115 sumneu[5]+=atten;
01116 sumneu[6]+=Q;
01117
01118 PHneutron+=phsmr;
01119 }
01120 else if(type == "Tl"){
01121 nthagt0++;
01122 sumtha[0]+=phsmr;
01123 sumtha[1]+=phgen;
01124 sumtha[2]+=phgen*atten0*sfactor;
01125 sumtha[3]+=npe;
01126 sumtha[4]+=sfactor;
01127 sumtha[5]+=atten;
01128 sumtha[6]+=Q;
01129
01130 PHthallium+=phsmr;
01131 }
01132 else if(type == "Bi"){
01133 nbisgt0++;
01134 sumbis[0]+=phsmr;
01135 sumbis[1]+=phgen;
01136 sumbis[2]+=phgen*atten0*sfactor;
01137 sumbis[3]+=npe;
01138 sumbis[4]+=sfactor;
01139 sumbis[5]+=atten;
01140 sumbis[6]+=Q;
01141
01142 PHbismuth+=phsmr;
01143 }
01144 else if(type == "Cf"){
01145 ncalgt0++;
01146 sumcal[0]+=phsmr;
01147 sumcal[1]+=phgen;
01148 sumcal[2]+=phgen*atten0*sfactor;
01149 sumcal[3]+=npe;
01150 sumcal[4]+=sfactor;
01151 sumcal[5]+=atten;
01152 sumcal[6]+=Q;
01153
01154 PHcalifornium+=phsmr;
01155 }
01156 else if(type == "muon"){
01157 nmuogt0++;
01158 summuo[0]+=phsmr;
01159 summuo[1]+=phgen;
01160 summuo[2]+=phgen*atten0*sfactor;
01161 summuo[3]+=npe;
01162 summuo[4]+=sfactor;
01163 summuo[5]+=atten;
01164 summuo[6]+=Q;
01165
01166 PHmuon+=phsmr;
01167 }
01168 else{
01169 cout << " ERROR: GenerateResponse has unrecognized type "
01170 << type << endl;
01171 }
01172
01173 PMT PMTtmp;
01174 PMTtmp.npmt = pmtindex;
01175 for(int l=0;l<3;l++)PMTtmp.xyz[l] = rpmt[l];
01176 PMTtmp.PHgen = phgen;
01177 PMTtmp.PHsmr = phsmr;
01178 PMTtmp.PathGd = sGd;
01179 PMTtmp.PathSc = sSc;
01180 PMTtmp.Atten = atten;
01181 PMTtmp.Sfactor = sfactor;
01182 PMTtmp.TimeGen = tgen;
01183 PMTtmp.TimeSmr = tsmr;
01184 PMTtmp.Xpe = xpe;
01185 PMTtmp.Npe = npe;
01186 PMTtmp.Q = Q;
01187 PMTtmp.Parent = type;
01188
01189 if(PMTptr<PMTdata+PMTsize){
01190 *PMTptr++=PMTtmp;
01191 PMTlen++;
01192 }
01193 else{
01194
01195
01196
01197 PMTptr = ResizePMT();
01198 *PMTptr++=PMTtmp;
01199 PMTlen++;
01200 }
01201 }
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213 if(sumpos[2]>0){
01214 PHposReco+=sumpos[0]*sumpos[1]/sumpos[2];
01215 MeanNpe[0]+=sumpos[3]/Npmt;
01216 MeanQ[0]+=sumpos[6]/Npmt;
01217 MeanSfactor[0]+=sumpos[4]/Npmt/nposgt0;
01218 MeanAtten[0]+=sumpos[5]/Npmt/nposgt0;
01219 }
01220 if(sumneu[2]>0){
01221 PHneuReco+=sumneu[0]*sumneu[1]/sumneu[2];
01222 MeanNpe[1]+=sumneu[3]/Npmt;
01223 MeanQ[1]+=sumneu[6]/Npmt;
01224 MeanSfactor[1]+=sumneu[4]/Npmt/nneugt0;
01225 MeanAtten[1]+=sumneu[5]/Npmt/nneugt0;
01226 }
01227 if(sumtha[2]>0){
01228 PHthaReco+=sumtha[0]*sumtha[1]/sumtha[2];
01229 MeanNpe[2]+=sumtha[3]/Npmt;
01230 MeanQ[2]+=sumtha[6]/Npmt;
01231 MeanSfactor[2]+=sumtha[4]/Npmt/nthagt0;
01232 MeanAtten[2]+=sumtha[5]/Npmt/nthagt0;
01233 }
01234 if(sumbis[2]>0){
01235 PHbisReco+=sumbis[0]*sumbis[1]/sumbis[2];
01236 MeanNpe[3]+=sumbis[3]/Npmt;
01237 MeanQ[3]+=sumbis[6]/Npmt;
01238 MeanSfactor[3]+=sumbis[4]/Npmt/nbisgt0;
01239 MeanAtten[3]+=sumbis[5]/Npmt/nbisgt0;
01240 }
01241 if(sumcal[2]>0){
01242 PHcalReco+=sumcal[0]*sumcal[1]/sumcal[2];
01243 MeanNpe[4]+=sumcal[3]/Npmt;
01244 MeanQ[4]+=sumcal[6]/Npmt;
01245 MeanSfactor[4]+=sumcal[4]/Npmt/ncalgt0;
01246 MeanAtten[4]+=sumcal[5]/Npmt/ncalgt0;
01247 }
01248 if(summuo[2]>0){
01249 PHmuoReco+=summuo[0]*summuo[1]/summuo[2];
01250 MeanNpe[5]+=summuo[3]/Npmt;
01251 MeanQ[5]+=summuo[6]/Npmt;
01252 MeanSfactor[5]+=summuo[4]/Npmt/nmuogt0;
01253 MeanAtten[5]+=summuo[5]/Npmt/nmuogt0;
01254 }
01255 if(sumele[2]>0){
01256 PHeleReco+=sumele[0]*sumele[1]/sumele[2];
01257 MeanNpe[6]+=sumele[3]/Npmt;
01258 MeanQ[6]+=sumele[6]/Npmt;
01259 MeanSfactor[6]+=sumele[4]/Npmt/nelegt0;
01260 MeanAtten[6]+=sumele[5]/Npmt/nelegt0;
01261 }
01262 if(sumpos[6]>0){
01263 HitPmtsPos++;
01264 }
01265 if(sumneu[6]>0){
01266 HitPmtsNeu++;
01267 }
01268 if(sumtha[6]>0){
01269 HitPmtsTha++;
01270 }
01271 if(sumbis[6]>0){
01272 HitPmtsBis++;
01273 }
01274 if(sumcal[6]>0){
01275 HitPmtsCal++;
01276 }
01277 if(summuo[6]>0){
01278 HitPmtsMuo++;
01279 }
01280 if(sumele[6]>0){
01281 HitPmtsEle++;
01282 }
01283
01284
01285
01286 pmtindex++;
01287 }
01288 }
01289 }
01290
01291 void ReactorDetector::GenerateGammas(ReactorEvent& Event){
01292
01293
01294
01295
01296 int ngamma = Event.GetNumGamma();
01297
01298
01299 int maxgamma = 100+ngamma;
01300 double* egamma = new double[maxgamma];
01301 double* tgamma = new double[maxgamma];
01302 double* rgamma = new double[3*maxgamma];
01303 std::string* ogamma = new std::string[maxgamma];
01304
01305 double* egptr=egamma;
01306 double* rgptr=rgamma;
01307 double* tgptr=tgamma;
01308 std::string* ogptr=ogamma;
01309
01310
01311
01312 for(int n=0;n<ngamma;n++){
01313
01314 *egptr++ = Event.GetGammaKE(n);
01315
01316 for(int i=0;i<3;i++)
01317 *rgptr++ = *(Event.GetGammaVertexXYZ()+(i+n*Event.Rdim));
01318
01319 *tgptr++ = Event.GetGammaTime(n);
01320
01321 if(Event.GetGammaParentProcess(n) == 1){
01322 cout << "ERROR: gammas are not produced by inverse-beta decay" << endl;
01323 }
01324 else if(Event.GetGammaParentProcess(n) == 2){
01325 *ogptr++ = "Cf";
01326 }
01327 else if(Event.GetGammaParentProcess(n) == 3){
01328 *ogptr++ = "muon";
01329 }
01330 else{
01331 cout << "Process Id " << Event.GetGammaParentProcess(n) << " is unknown"
01332 << endl;
01333 }
01334 }
01335
01336
01337
01338 if(Event.GetPmtRad() == "on"){
01339
01340
01341 int len = 2;
01342 float r[len];
01343 ranlux_(r,len);
01344 int pmtId = r[0]*Npmt;
01345
01346
01347 std::string Stuff; int Isotope;
01348
01349
01350
01351 if(r[1] < 0.25){
01352 Stuff = "thallium"; Isotope = 208;
01353 }
01354 else if(0.25 <= r[1] && r[1] < 0.5){
01355 Stuff = "thallium"; Isotope = 210;
01356 }
01357 else if(0.5 <= r[1] && r[1] < 0.75){
01358 Stuff = "bismuth"; Isotope = 212;
01359 }
01360 else{
01361 Stuff = "bismuth"; Isotope = 214;
01362 }
01363
01364
01365
01366
01367 if(Stuff == "thallium"){
01368 MyTlSource Thallium(Isotope);
01369
01370
01371 float PMTWeight[Thallium.GetNumTlGammas()];
01372
01373
01374 for(int n=0;n<Thallium.GetNumTlGammas();n++){
01375
01376
01377 float e = Thallium.GetTlGammaEnergy(n);
01378 *egptr++ = e;
01379
01380
01381 PMTWeight[n] = 1.;
01382 int attempts = 1;
01383 double r[3];
01384 double time;
01385
01386
01387 int len=4;
01388 float rv[len];
01389 ranlux_(rv,len);
01390 double cosg = -1+2*rv[0];
01391 double phig = 2*RC.pi*rv[1];
01392 double t[3];
01393 t[0] = cos(phig)*sqrt(1-cosg*cosg);
01394 t[1] = sin(phig)*sqrt(1-cosg*cosg);
01395 t[2] = cosg;
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405 double range;
01406 std::string material = "oil";
01407 range = GammaComptonRange(material,e);
01408 range*=log(1/rv[2]);
01409
01410
01411 r[0] = Xpmt[pmtId]+range*t[0];
01412 r[1] = Ypmt[pmtId]+range*t[1];
01413 r[2] = Zpmt[pmtId]+range*t[2];
01414
01415 double rdotr = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
01416 if(sqrt(rdotr)<R1){
01417 pmtRadCount++;
01418 cout << "PMTrad Tl " << pmtRadCount << " radius " << sqrt(rdotr)
01419 << ", energy " << e << endl;
01420 }
01421
01422
01423 for(int i=0;i<3;i++) *rgptr++ = r[i];
01424
01425
01426 time = rv[3]*1000000.;
01427
01428 *tgptr++ = time;
01429 *ogptr++ = "Tl";
01430
01431 PMTWeight[n] = 1./float(attempts);
01432
01433
01434 ngamma++;
01435 }
01436 }
01437
01438
01439
01440 else if(Stuff == "bismuth"){
01441 MyBiSource Bismuth(Isotope);
01442
01443
01444 float PMTWeight[Bismuth.GetNumBiGammas()];
01445
01446
01447 for(int n=0;n<Bismuth.GetNumBiGammas();n++){
01448
01449
01450 float e = Bismuth.GetBiGammaEnergy(n);
01451 *egptr++ = e;
01452
01453
01454 PMTWeight[n] = 1.;
01455 int attempts = 1;
01456 double r[3];
01457 double time;
01458
01459
01460 int len=4;
01461 float rv[len];
01462 ranlux_(rv,len);
01463 double cosg = -1+2*rv[0];
01464 double phig = 2*RC.pi*rv[1];
01465 double t[3];
01466 t[0] = cos(phig)*sqrt(1-cosg*cosg);
01467 t[1] = sin(phig)*sqrt(1-cosg*cosg);
01468 t[2] = cosg;
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478 double range;
01479 std::string material = "oil";
01480 range = GammaComptonRange(material,e);
01481 range*=log(1/rv[2]);
01482
01483
01484 r[0] = Xpmt[pmtId]+range*t[0];
01485 r[1] = Ypmt[pmtId]+range*t[1];
01486 r[2] = Zpmt[pmtId]+range*t[2];
01487
01488 double rdotr = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
01489 if(sqrt(rdotr)<R1){
01490 pmtRadCount++;
01491 cout << "PMTrad Bi " << pmtRadCount << " radius " << sqrt(rdotr)
01492 << ", energy " << e << endl;
01493 }
01494
01495
01496 for(int i=0;i<3;i++) *rgptr++ = r[i];
01497
01498
01499 time = rv[3]*1000000.;
01500
01501 *tgptr++ = time;
01502 *ogptr++ = "Bi";
01503
01504 PMTWeight[n] = 1./float(attempts);
01505
01506
01507 ngamma++;
01508 }
01509 }
01510 }
01511
01512
01513
01514 int Nelectron = Event.GetNumElectron();
01515
01516
01517
01518 for(int n=0;n<Nelectron;n++){
01519
01520 double rele[3];
01521 for(int i=0;i<3;i++){
01522 rele[i] = *(Event.GetElectronVertexXYZ()+(i+n*Event.Rdim));
01523
01524 }
01525 double tele = Event.GetElectronTime(n);
01526
01527 double keele = Event.GetElectronKE(n);
01528
01529
01530
01531
01532 for(int i=0;i<3;i++)*rgptr++=rele[i];
01533 *tgptr++ = tele;
01534 *egptr++ = keele;
01535 *ogptr++ = "E";
01536 ngamma++;
01537 }
01538
01539
01540
01541
01542 int Npositron = Event.GetNumPositron();
01543
01544
01545
01546 for(int n=0;n<Npositron;n++){
01547
01548 double rpos[3];
01549 for(int i=0;i<3;i++){
01550 rpos[i] = *(Event.GetPositronVertexXYZ()+(i+n*Event.Rdim));
01551
01552 }
01553 double tpos = Event.GetPositronTime(n);
01554
01555 double kepos = Event.GetPositronKE(n);
01556
01557
01558
01559
01560
01561 for(int i=0;i<3;i++)*rgptr++=rpos[i];
01562 *tgptr++ = tpos;
01563 *egptr++ = kepos;
01564 *ogptr++ = "P";
01565 ngamma++;
01566
01567
01568
01569 int len=4;
01570 float rv[4];
01571 ranlux_(rv,len);
01572 double cosg = -1+2*rv[0];
01573 double phig = 2*RC.pi*rv[1];
01574 double t[3];
01575 t[0] = cos(phig)*sqrt(1-cosg*cosg);
01576 t[1] = sin(phig)*sqrt(1-cosg*cosg);
01577 t[2] = cosg;
01578 double range;
01579 for (int i=0;i<2;i++){
01580 std::string material = "scintillator";
01581 range = GammaComptonRange(material,RC.Melectron);
01582 range*=log(1/rv[i+2]);
01583 if(i==0){
01584 for(int j=0;j<3;j++) *rgptr++ = rpos[j]+range*t[j];
01585 }
01586 if(i==1){
01587 for(int j=0;j<3;j++)*rgptr++ = rpos[j]-range*t[j];
01588 }
01589 *egptr++=RC.Melectron;
01590 *tgptr++=tpos+range/RC.c*refracSc;
01591 *ogptr++ = "P";
01592 ngamma++;
01593 }
01594 }
01595
01596
01597
01598
01599 Nneutron = Event.GetNumNeutron();
01600
01601
01602
01603 for(int n=0;n<Nneutron;n++){
01604
01605 double rneu[3];
01606 for(int i=0;i<3;i++){
01607 rneu[i] = *(Event.GetNeutronVertexXYZ()+(i+n*Event.Rdim));
01608
01609 }
01610 double tneu = Event.GetNeutronTime(n);
01611
01612
01613
01614
01615 double emax=RC.Hpeak;
01616 if(FastNeutronOption){
01617
01618 Absorber[n]=1;
01619 double rrneu = sqrt(rneu[0]*rneu[0]+rneu[1]*rneu[1]+rneu[2]*rneu[2]);
01620 if(rrneu<R0){
01621 int len=2;
01622 float rv[2];
01623 ranlux_(rv,len);
01624 if(rv[0]<=GdCaptureFraction){
01625
01626
01627
01628 if(rv[1]<=RC.Gd155frac){
01629 emax = RC.Gd155peak;
01630 Absorber[n]=155;
01631 }
01632 else{
01633 emax = RC.Gd157peak;
01634 Absorber[n]=157;
01635 }
01636 }
01637 }
01638 }
01639 else{
01640
01641 if(Absorber[n]==1)emax=RC.Hpeak;
01642 if(Absorber[n]==12)emax=RC.Cpeak;
01643 if(Absorber[n]==155)emax=RC.Gd155peak;
01644 if(Absorber[n]==157)emax=RC.Gd157peak;
01645 }
01646
01647
01648
01649 if(emax==RC.Hpeak || emax==RC.Cpeak){
01650
01651 int len=3;
01652 float rv[3];
01653 ranlux_(rv,len);
01654 std::string material = "scintillator";
01655 double range = GammaComptonRange(material,emax);
01656 range*=log(1/rv[0]);
01657 double cosg = -1+2*rv[1];
01658 double phig = 2*RC.pi*rv[2];
01659 double t[3];
01660 t[0] = cos(phig)*sqrt(1-cosg*cosg);
01661 t[1] = sin(phig)*sqrt(1-cosg*cosg);
01662 t[2] = cosg;
01663 for(int j=0;j<3;j++)*rgptr++ = rneu[j]+range*t[j];
01664 *egptr++ = emax;
01665 *tgptr++ = tneu+range/RC.c*refracSc;
01666 *ogptr++ = "N";
01667 ngamma++;
01668
01669 }
01670
01671
01672
01673 else{
01674
01675 int ngd;
01676 double pgd[3*maxgamma];
01677 MyGdDecayModel(emax,ngd,pgd);
01678 float rv2[ngd];
01679 ranlux_(rv2,ngd);
01680 double* qgd=pgd;
01681 double check[4]={0,0,0,0};
01682 for(int i=0;i<ngd;i++){
01683 double p[3];
01684 for(int j=0;j<3;j++)p[j]=*qgd++;
01685 double e = sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
01686 std::string material = "scintillator";
01687 double range = GammaComptonRange(material,e);
01688 range*=log(1/rv2[i]);
01689 for(int j=0;j<3;j++)*rgptr++ = rneu[j]+range*p[j]/e;
01690
01691 for(int j=0;j<3;j++)check[j]+=p[j];
01692 check[3]+=e;
01693
01694 *egptr++ = e;
01695 *tgptr++ = tneu+range/RC.c*refracGd;
01696 *ogptr++ = "N";
01697 ngamma++;
01698
01699 }
01700 }
01701 }
01702
01703
01704
01705
01706 delete[] Rgamma;
01707 delete[] Egamma;
01708 delete[] Tgamma;
01709 delete[] Ogamma;
01710
01711 Ngamma = ngamma;
01712 Rgamma = new double[3*ngamma];
01713 Egamma = new double[ngamma];
01714 Tgamma = new double[ngamma];
01715 Ogamma = new std::string[ngamma];
01716
01717 double* ptr;
01718
01719 ptr = Rgamma+3*ngamma;
01720 while(ptr>Rgamma)*--ptr=*--rgptr;
01721
01722 ptr = Egamma+ngamma;
01723 while(ptr>Egamma)*--ptr=*--egptr;
01724
01725 ptr = Tgamma+ngamma;
01726 while(ptr>Tgamma)*--ptr=*--tgptr;
01727
01728 std::string* qtr = Ogamma+ngamma;
01729 while(qtr>Ogamma)*--qtr=*--ogptr;
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744 delete[] rgamma;
01745 delete[] egamma;
01746 delete[] tgamma;
01747 delete[] ogamma;
01748 }
01749
01750 float ReactorDetector::QSpectrum(const float& Q){
01751 float m = 8.1497;
01752 float ng1 = 0.2735E-01;
01753 float ng2 = 0.7944E-03;
01754 float nexp = 0.81228E-01;
01755 float G2 = 4.5;
01756 float Q2 = Q/G2;
01757 float Q0 = 0.61;
01758 float gmratio = m/exp(gamma(m));
01759
01760
01761 float Qpol1 = gmratio*pow((m*Q),m-1)*exp(-m*Q);
01762 float Qpol2 = gmratio*pow((m*Q2),m-1)*exp(-m*Q2);
01763 return (ng1*Qpol1 + ng2*Qpol2 + nexp*exp(-Q/Q0));
01764 }
01765
01766 float ReactorDetector::GenerateQ(int npe){
01767
01768 float pedwidth = 0.1;
01769 float qtot = 0.;
01770
01771 if (npe>0) {
01772 float qran[npe];
01773 float qsmear[npe];
01774 rnormx_(qsmear,npe,ranlux_);
01775 funlux_(Qspace,qran[0],npe);
01776 for (int ipe=0; ipe<npe; ipe++){
01777 qtot+=(qran[ipe]+pedwidth*qsmear[ipe]);
01778 }
01779 }
01780 else {
01781 float qsmear[1];
01782 int ipe = 1;
01783 rnormx_(qsmear,ipe,ranlux_);
01784 qtot+=pedwidth*qsmear[1];
01785 }
01786 return qtot;
01787 }
01788
01789
01790
01791 ReactorDetector::PMT* ReactorDetector::ResizePMT(){
01792
01793 int oldsize = PMTsize;
01794 int newsize = 2*oldsize;
01795 PMT* PMThold = new PMT[newsize];
01796 PMT* PMTold = PMTdata+oldsize;
01797 PMT* PMTnew = PMThold+oldsize;
01798 while(PMTnew>PMThold)*--PMTnew=*--PMTold;
01799 delete[] PMTdata;
01800 PMTdata = new PMT[newsize];
01801 PMTold = PMThold+newsize;
01802 PMTnew = PMTdata+newsize;
01803 while(PMTnew>PMTdata)*--PMTnew=*--PMTold;
01804 delete[] PMThold;
01805 PMTsize = newsize;
01806 return PMTdata+oldsize;
01807 }
01808
01809
01810
01811 int ReactorDetector::GetNgamma(std::string parent){
01812 int gamma=0;
01813 for(std::string* optr=Ogamma;optr<Ogamma+Ngamma;optr++){
01814 if (*optr==parent)gamma++;
01815 }
01816 return gamma;
01817 }
01818
01819
01820
01821 void ReactorDetector::GetGamma(std::string parent,double* xyz){
01822
01823 double* Rptr=Rgamma;
01824 double* rptr=xyz;
01825
01826 for(std::string* optr=Ogamma;optr<Ogamma+Ngamma;optr++){
01827 if (*optr==parent){
01828 for(int i=0;i<3;i++)*rptr++ = *Rptr++;
01829 }
01830 }
01831 }
01832 void ReactorDetector::GetGamma(std::string parent,double* e,
01833 double* t,double* r,
01834 double* c,double* f){
01835
01836
01837 double* Eptr=Egamma;
01838 double* Rptr=Rgamma;
01839 double* Tptr=Tgamma;
01840 double* eptr=e;
01841 double* rptr=r;
01842 double* tptr=t;
01843 double* cptr=c;
01844 double* fptr=f;
01845
01846 for(std::string* optr=Ogamma;optr<Ogamma+Ngamma;optr++){
01847 if (*optr==parent){
01848 *e++=*Eptr;
01849 *t++=*Tptr;
01850 double x=*(Rptr+0);
01851 double y=*(Rptr+1);
01852 double z=*(Rptr+2);
01853 *r++=sqrt(x*x+y*y+z*z);
01854 *c++=z/sqrt(x*x+y*y+z*z);
01855 *f++=atan2(y,x);
01856 if((*f)<0)*f+=2*RC.pi;
01857 }
01858 Eptr++;
01859 Tptr++;
01860 Rptr+=3;
01861 }
01862 }
01863 void ReactorDetector::SetNominalParameters(){
01864
01865
01866
01867 R0 = RC.R0;
01868 R1 = RC.R1;
01869 R2 = RC.R2;
01870
01871 GdConcentration=RC.GdConcentration;
01872
01873 refracGd = RC.refracGd;
01874 refracSc = RC.refracSc;
01875
01876 mfpGd = RC.mfpGd;
01877 mfpSc = RC.mfpSc;
01878
01879 MeanNeuDispl = RC.MeanNeuDispl;
01880
01881 FastNeutronOption = RC.FastNeutronOption;
01882 Temperature = RC.Temperature;
01883
01884 tGd = RC.tGd;
01885 tSc = RC.tSc;
01886
01887 GdCaptureFraction = RC.GdCaptureFraction;
01888
01889
01890
01891
01892
01893 GammasPerGd = RC.GammasPerGd;
01894
01895 attenlGd = RC.attenlGd;
01896 attenlSc = RC.attenlSc;
01897
01898 PhotonsPerMeV = RC.PhotonsPerMeV;
01899
01900 ShortPosFrac0 = RC.ShortPosFrac0;
01901 ShortPosFrac1 = RC.ShortPosFrac1;
01902 LongPosFrac = RC.LongPosFrac;
01903 ShortPosTime0 = RC.ShortPosTime0;
01904 ShortPosTime1 = RC.ShortPosTime1;
01905 LongPosTime = RC.LongPosTime;
01906
01907 PMTcoverage = RC.PMTcoverage;
01908 PMTdiameter = RC.PMTdiameter;
01909
01910 PMTqe = RC.PMTqe;
01911 EScale = RC.EScale;
01912 }
01913 void ReactorDetector::SetPMT(){
01914
01915
01916
01917
01918 Npmt = 16*R2*R2*PMTcoverage/PMTdiameter/PMTdiameter;
01919 Ncosbins = sqrt(Npmt/RC.pi);
01920 Nphibins = sqrt(Npmt*RC.pi);
01921 Npmt = Ncosbins*Nphibins;
01922
01923
01924 delete[] Xpmt;
01925 delete[] Ypmt;
01926 delete[] Zpmt;
01927
01928 Xpmt = new double[Npmt];
01929 Ypmt = new double[Npmt];
01930 Zpmt = new double[Npmt];
01931
01932
01933
01934
01935 delete[] PMTdata;
01936 PMTlen = 0;
01937 PMTsize = 10*Npmt;
01938 PMTdata = new PMT[PMTsize];
01939
01940
01941
01942 double* xpmt=Xpmt;
01943 double* ypmt=Ypmt;
01944 double* zpmt=Zpmt;
01945
01946 double dphi = 2*RC.pi/Nphibins;
01947 double dz = 2.0/Ncosbins;
01948
01949 for(int i=0;i<Ncosbins;i++){
01950 double z = -1+i*dz+dz/2;
01951 for(int j=0;j<Nphibins;j++){
01952 double phi = j*dphi+dphi/2;
01953 *xpmt++ = R2*sqrt(1-z*z)*cos(phi);
01954 *ypmt++ = R2*sqrt(1-z*z)*sin(phi);
01955 *zpmt++ = R2*z;
01956
01957 }
01958 }
01959 }
01960 void ReactorDetector::SetParam_GdConcentration(double x){
01961 GdConcentration =x;
01962 double y =(GdConcentration/RC.GdConcentrationRef);
01963 double z = GdCaptureFraction/(1-GdCaptureFraction);
01964
01965 GdCaptureFraction = y*z/(1+y*z);
01966
01967 mfpGd*=(1+z)/(1+y*z);
01968 tGd*=(1+z)/(1+y*z);
01969 SetNeutronPropStuff();
01970 }
01971 void ReactorDetector::SetNeutronPropStuff(){
01972
01973 a[0] = RC.A0;
01974 a[1] = RC.A1;
01975 a[2] = RC.A2;
01976 a[3] = RC.A3;
01977
01978 z[0] = RC.Z0;
01979 z[1] = RC.Z1;
01980 z[2] = RC.Z2;
01981 z[3] = RC.Z3;
01982
01983 nn[0]=RC.Navogadro*RC.f0*RC.density/RC.A0;
01984 nn[1]=RC.Navogadro*RC.f1*RC.density/RC.A1;
01985 nn[2]=RC.Navogadro*GdConcentration/100*RC.f2*RC.density/RC.A2;
01986 nn[3]=RC.Navogadro*GdConcentration/100*RC.f3*RC.density/RC.A3;
01987
01988 }
01989 double ReactorDetector::GetHitFrac(std::string parent){
01990 bool hit[Npmt];
01991 for(int i=0;i<Npmt;i++)hit[i]=0;
01992 for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
01993 if(pPMT->Parent==parent&&pPMT->Npe>0)hit[pPMT->npmt]=1;
01994 }
01995 int nhit=0;
01996 for(int i=0;i<Npmt;i++)nhit+=hit[i];
01997 return 1.0*nhit/Npmt;
01998 }
01999
02000 void ReactorDetector::CalculateEnergy(std::string parent){
02001
02002
02003
02004 double rpmt[3];
02005 double revent[3];
02006 if(parent=="P")
02007 for(int i=0;i<3;i++)revent[i]=DipolePositron[i];
02008 else if(parent=="N")
02009 for(int i=0;i<3;i++)revent[i]=DipoleNeutron[i];
02010 else
02011 cout << "ERROR: CalculateEnergy does not work for parent type: "
02012 << parent << endl;
02013
02014 double rdotr = revent[0]*revent[0]+revent[1]*revent[1]+revent[2]*revent[2];
02015 double* xpmt=Xpmt;
02016 double* ypmt=Ypmt;
02017 double* zpmt=Zpmt;
02018 double dr[3];
02019 double tr[3];
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029 double detprob = 0.;
02030 double sGd;
02031 double sSc;
02032 for(int i=0;i<Ncosbins;i++){
02033 for(int j=0;j<Nphibins;j++){
02034 double rpmt[3];
02035 rpmt[0] = *xpmt++;
02036 rpmt[1] = *ypmt++;
02037 rpmt[2] = *zpmt++;
02038 for (int l=0;l<3;l++)dr[l]=rpmt[l]-revent[l];
02039 double ddr = sqrt(dr[0]*dr[0]+dr[1]*dr[1]+dr[2]*dr[2]);
02040 for (int l=0;l<3;l++)tr[l]=dr[l]/ddr;
02041 double tdotr = revent[0]*tr[0]+revent[1]*tr[1]+revent[2]*tr[2];
02042
02043
02044
02045
02046
02047
02048 double test = R0*R0+tdotr*tdotr-rdotr;
02049 bool inter = (test>=0);
02050
02051
02052
02053
02054
02055 if(inter){
02056 double sint1 = -tdotr+sqrt(test);
02057 double sint2 = -tdotr-sqrt(test);
02058 if(sint1>0&&sint2<0){
02059 sGd = sint1;
02060 }
02061 else if(sint2>0&&sint1<0){
02062 sGd = sint2;
02063 }
02064 else if(sint2>0&&sint1>0){
02065 sGd = sint1-sint2;
02066 if(sGd<0)sGd*=-1;
02067 }
02068 else{
02069 sGd = 0;
02070 }
02071 sSc = ddr-sGd;
02072 }
02073
02074
02075
02076 else{
02077 sGd = 0;
02078 sSc = ddr;
02079 }
02080
02081
02082 double atten = exp(-sGd/attenlGd-sSc/attenlSc);
02083 double sfactor=0;
02084 for (int l=0;l<3;l++){sfactor+=tr[l]*rpmt[l]/R2;}
02085 sfactor *= R2*R2/ddr/ddr;
02086 detprob+=sfactor*atten*PMTqe*(PMTcoverage/Npmt);
02087 }
02088 }
02089
02090 double q = 0.;
02091 for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
02092 if(pPMT->Parent==parent&&pPMT->Q>0){
02093 q+=pPMT->Q;
02094 }
02095 }
02096
02097
02098
02099 double energy = 0.;
02100 energy = q/(detprob*PhotonsPerMeV*EScale);
02101
02102 if(parent=="E") PHeleReco=energy;
02103 if(parent=="P") PHposReco=energy;
02104 if(parent=="N") PHneuReco=energy;
02105 if(parent=="Tl") PHthaReco=energy;
02106 if(parent=="Bi") PHbisReco=energy;
02107 if(parent=="Cf") PHcalReco=energy;
02108 if(parent=="muon") PHmuoReco=energy;
02109
02110 }
02111
02112 void ReactorDetector::CalculateQPosition(std::string parent){
02113
02114 double d[3] = {0,0,0};
02115 double q=0;
02116 for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
02117 if(pPMT->Parent==parent&&pPMT->Q>0){
02118 q+=pPMT->PHsmr;
02119 for(int i=0;i<3;i++)d[i]+=(pPMT->PHsmr)*(pPMT->xyz[i]);
02120 }
02121 }
02122 if(q>0)for(int i=0;i<3;i++)d[i]/=q;
02123
02124 if(parent=="P"){
02125 for(int i=0;i<3;i++)DipolePositron[i] = d[i];
02126 }
02127 if(parent=="N"){
02128 for(int i=0;i<3;i++)DipoleNeutron[i] = d[i];
02129 }
02130 }
02131 void ReactorDetector::ImproveQPosition(std::string parent){
02132
02133 double d[3];
02134 double k;
02135 PMT* qPMT;
02136 if(parent=="P"){
02137 for(int i=0;i<3;i++)d[i]=DipolePositron[i];
02138 k = PHpositron;
02139 qPMT=PMTpositron;
02140 }
02141 if(parent=="N"){
02142 for(int i=0;i<3;i++)d[i]=DipoleNeutron[i];
02143 k = PHneutron;
02144 qPMT=PMTneutron;
02145 }
02146
02147 double C= RC.PMTqe*RC.PMTdiameter*RC.PMTdiameter/16*RC.PhotonsPerMeV;
02148
02149 double sumlast = 999999;
02150 double test=999999;
02151 double tol=0.01;
02152 int sign=1;
02153 bool first=1;
02154 double dev=0.1;
02155 while(test>tol){
02156
02157 double dlast[3];
02158 for(int i=0;i<3;i++)dlast[i]=d[i];
02159 double klast=k;
02160 double sum[3]={0,0,0};
02161 for(PMT* pPMT=qPMT;pPMT<qPMT+Npmt;pPMT++){
02162 if(pPMT->Parent!=parent)continue;
02163 double r[3];
02164 double rn[3];
02165 double t[3];
02166 for(int i=0;i<3;i++)rn[i]=pPMT->xyz[i];
02167 for(int i=0;i<3;i++)r[i]=rn[i]-d[i];
02168 double rdotr = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
02169 double rdotrn = r[0]*rn[0]+r[1]*rn[1]+r[2]*rn[2];
02170 for(int i=0;i<3;i++)t[i]=r[i]/sqrt(rdotr);
02171
02172 double tdotd=d[0]*t[0]+d[1]*t[1]+d[2]*t[2];
02173 double ddotd=d[0]*d[0]+d[1]*d[1]+d[2]*d[2];
02174 double test2=tdotd*tdotd-ddotd+R0*R0;
02175 double sSc;
02176 double sGd;
02177
02178 if(test2<0){
02179 sSc=sqrt(rdotr);
02180 sGd=0;
02181 }
02182 else{
02183 double s1=-tdotd+sqrt(test2);
02184 double s2=-tdotd-sqrt(test2);
02185 if(s1>0&&s2>0){
02186 sSc=sqrt(rdotr)-s1+s2;
02187 sGd=s1-s2;
02188 }
02189 else if(s1>0&&s2<0){
02190 sSc=sqrt(rdotr)-s1;
02191 sGd=s1;
02192 }
02193 else {
02194 sSc=sqrt(rdotr);
02195 sGd=0;
02196 }
02197 }
02198
02199 double mn=pPMT->Q;
02200 double w0=C*rdotrn/R2/sqrt(rdotr)/rdotr*
02201 exp(-sGd/RC.attenlGd-sSc/RC.attenlSc);
02202 sum[0]+=mn;
02203 sum[1]+=w0;
02204 sum[2]+=w0*klast-mn*log(w0*klast);
02205 }
02206 if(first){
02207 k = sum[0]/sum[1];
02208 for(int i=0;i<3;i++)d[i]*=(1+dev);
02209 first=0;
02210 }
02211 else{
02212 k = sum[0]/sum[1];
02213 if(sum[2]<sumlast){
02214 for(int i=0;i<3;i++)d[i]*=(1+sign*dev);
02215 }
02216 else{
02217 sign*=-1;
02218 dev/=2;
02219 for(int i=0;i<3;i++)d[i]*=(1+sign*dev);
02220 }
02221 }
02222 test=1-sum[2]/sumlast;
02223 if(test<0)test*=-1;
02224 sumlast=sum[2];
02225 }
02226 if(parent=="P")for(int i=0;i<3;i++)DipolePositron[i]=d[i];
02227 if(parent=="N")for(int i=0;i<3;i++)DipoleNeutron[i]=d[i];
02228 }
02229
02230
02231
02232 void ReactorDetector::ConsolidatePMT(){
02233
02234
02235
02236 double xpeele[Npmt];
02237 double Qele[Npmt];
02238 double xpepos[Npmt];
02239 double Qpos[Npmt];
02240 double xpeneu[Npmt];
02241 double Qneu[Npmt];
02242 double xpetha[Npmt];
02243 double Qtha[Npmt];
02244 double xpebis[Npmt];
02245 double Qbis[Npmt];
02246 double xpecal[Npmt];
02247 double Qcal[Npmt];
02248 double xpemuo[Npmt];
02249 double Qmuo[Npmt];
02250 int subhitsele[Npmt];
02251 int subhitspos[Npmt];
02252 int subhitsneu[Npmt];
02253 int subhitstha[Npmt];
02254 int subhitsbis[Npmt];
02255 int subhitscal[Npmt];
02256 int subhitsmuo[Npmt];
02257
02258 for(int i=0;i<Npmt;i++){
02259 xpeele[i]=0;
02260 xpepos[i]=0;
02261 xpeneu[i]=0;
02262 xpetha[i]=0;
02263 xpebis[i]=0;
02264 xpecal[i]=0;
02265 xpemuo[i]=0;
02266 Qele[i]=0;
02267 Qpos[i]=0;
02268 Qneu[i]=0;
02269 Qtha[i]=0;
02270 Qbis[i]=0;
02271 Qcal[i]=0;
02272 Qmuo[i]=0;
02273 subhitsele[i]=0;
02274 subhitspos[i]=0;
02275 subhitsneu[i]=0;
02276 subhitstha[i]=0;
02277 subhitsbis[i]=0;
02278 subhitscal[i]=0;
02279 subhitsmuo[i]=0;
02280 }
02281
02282 for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
02283 int n=pPMT->npmt;
02284 if(pPMT->Parent=="E"){
02285 xpeele[n]+=pPMT->Npe;
02286 Qele[n]+=pPMT->Q;
02287 subhitsele[n]++;
02288 }
02289 if(pPMT->Parent=="P"){
02290 xpepos[n]+=pPMT->Npe;
02291 Qpos[n]+=pPMT->Q;
02292 subhitspos[n]++;
02293 }
02294 if(pPMT->Parent=="N"){
02295 xpeneu[n]+=pPMT->Npe;
02296 Qneu[n]+=pPMT->Q;
02297 subhitsneu[n]++;
02298 }
02299 if(pPMT->Parent=="Tl"){
02300 xpetha[n]+=pPMT->Npe;
02301 Qtha[n]+=pPMT->Q;
02302 subhitstha[n]++;
02303 }
02304 if(pPMT->Parent=="Bi"){
02305 xpebis[n]+=pPMT->Npe;
02306 Qbis[n]+=pPMT->Q;
02307 subhitsbis[n]++;
02308 }
02309 if(pPMT->Parent=="Cf"){
02310 xpecal[n]+=pPMT->Npe;
02311 Qcal[n]+=pPMT->Q;
02312 subhitscal[n]++;
02313 }
02314 if(pPMT->Parent=="muon"){
02315 xpemuo[n]+=pPMT->Npe;
02316 Qmuo[n]+=pPMT->Q;
02317 subhitsmuo[n]++;
02318 }
02319 }
02320
02321 delete[] PMTelectron;
02322 PMTelectron = new PMT[Npmt];
02323 delete[] PMTpositron;
02324 PMTpositron = new PMT[Npmt];
02325 delete[] PMTneutron;
02326 PMTneutron = new PMT[Npmt];
02327 delete[] PMTthallium;
02328 PMTthallium = new PMT[Npmt];
02329 delete[] PMTbismuth;
02330 PMTbismuth = new PMT[Npmt];
02331 delete[] PMTcalifornium;
02332 PMTcalifornium = new PMT[Npmt];
02333 delete[] PMTmuon;
02334 PMTmuon = new PMT[Npmt];
02335
02336 for(int i=0;i<Npmt;i++){
02337 (*(PMTelectron+i)).npmt=i;
02338 (*(PMTelectron+i)).Parent="E";
02339 (PMTelectron+i)->SetSubHitsPtr(subhitsele[i]);
02340
02341 (*(PMTpositron+i)).npmt=i;
02342 (*(PMTpositron+i)).Parent="P";
02343 (PMTpositron+i)->SetSubHitsPtr(subhitspos[i]);
02344
02345 (*(PMTneutron+i)).npmt=i;
02346 (*(PMTneutron+i)).Parent="N";
02347 (PMTneutron+i)->SetSubHitsPtr(subhitsneu[i]);
02348
02349 (*(PMTthallium+i)).npmt=i;
02350 (*(PMTthallium+i)).Parent="Tl";
02351 (PMTthallium+i)->SetSubHitsPtr(subhitstha[i]);
02352
02353 (*(PMTbismuth+i)).npmt=i;
02354 (*(PMTbismuth+i)).Parent="Tl";
02355 (PMTbismuth+i)->SetSubHitsPtr(subhitsbis[i]);
02356
02357 (*(PMTcalifornium+i)).npmt=i;
02358 (*(PMTcalifornium+i)).Parent="Tl";
02359 (PMTcalifornium+i)->SetSubHitsPtr(subhitscal[i]);
02360
02361 (*(PMTmuon+i)).npmt=i;
02362 (*(PMTmuon+i)).Parent="Tl";
02363 (PMTmuon+i)->SetSubHitsPtr(subhitsmuo[i]);
02364 }
02365
02366 for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
02367 int n=pPMT->npmt;
02368 PMT* qPMT;
02369 if(pPMT->Parent=="E"){
02370 qPMT = PMTelectron+n;
02371 *qPMT+=*pPMT;
02372 }
02373 if(pPMT->Parent=="P"){
02374 qPMT = PMTpositron+n;
02375 *qPMT+=*pPMT;
02376 }
02377 if(pPMT->Parent=="N"){
02378 qPMT = PMTneutron+n;
02379 *qPMT+=*pPMT;
02380 }
02381 if(pPMT->Parent=="Tl"){
02382 qPMT = PMTthallium+n;
02383 *qPMT+=*pPMT;
02384 }
02385 if(pPMT->Parent=="Bi"){
02386 qPMT = PMTbismuth+n;
02387 *qPMT+=*pPMT;
02388 }
02389 if(pPMT->Parent=="Cf"){
02390 qPMT = PMTcalifornium+n;
02391 *qPMT+=*pPMT;
02392 }
02393 if(pPMT->Parent=="muon"){
02394 qPMT = PMTmuon+n;
02395 *qPMT+=*pPMT;
02396 }
02397 }
02398 cout;
02399 }