Main Page   Compound List   File List   Compound Members   File Members  

ReactorDetector.cpp

Go to the documentation of this file.
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 "DetectorDefaults.hh"
00018 #include "ReactorXC.hh"
00019 #include "ReactorDetector.hh"
00020 #include "ReactorScint.hh"
00021 #include "MyTlSource.hh"
00022 #include "MyBiSource.hh"
00023 //
00024 // constructor
00025 //
00026 ReactorDetector::ReactorDetector(){
00027   //
00028   // initialize PMT arrays
00029   //
00030   Xpmt = new double[0];
00031   Ypmt = new double[0];
00032   Zpmt = new double[0];
00033   PMTOutputArray = new PMToutput*[0];
00034   PMTdata = new PMT[0];
00035   PMTelectron = new PMT[0];
00036   PMTpositron = new PMT[0];
00037   PMTneutron = new PMT[0];
00038   PMTthallium = new PMT[0];
00039   PMTbismuth = new PMT[0];
00040   PMTcalifornium = new PMT[0];
00041   PMTmuon = new PMT[0];
00042   //
00043   // Set parameters
00044   //
00045   SetNominalParameters();
00046   SetPMT();
00047   SetNeutronPropStuff();
00048   //
00049   // initialize gamma arrays
00050   //
00051   Ngamma = 0;
00052   Rgamma = new double[Ngamma];
00053   Egamma = new double[Ngamma];
00054   Tgamma = new double[Ngamma];
00055   Ogamma = new std::string[Ngamma];
00056   //
00057   // initialize PMT charge spectrum
00058   //
00059   float Qlo = 0.001;
00060   float Qhi = 10.0;
00061   funlxp_(QSpectrum,Qspace,Qlo,Qhi);
00062 
00063   Nneutron = 0;
00064   Absorber = new double[Nneutron];
00065   NeutronSteps = new double[Nneutron];
00066 
00067   pmtRadCount = 0;
00068   MergeTime = 0;
00069 }
00070 //
00071 // destructor
00072 //
00073 ReactorDetector::~ReactorDetector(){
00074   delete[] Xpmt;
00075   delete[] Ypmt;
00076   delete[] Zpmt;
00077   delete[] PMTOutputArray;
00078   delete[] PMTdata;
00079   delete[] PMTelectron;
00080   delete[] PMTpositron;
00081   delete[] PMTneutron;
00082   delete[] PMTthallium;
00083   delete[] PMTbismuth;
00084   delete[] PMTcalifornium;
00085   delete[] PMTmuon;
00086   delete[] Rgamma;
00087   delete[] Egamma;
00088   delete[] Tgamma;
00089   delete[] Ogamma;
00090   delete[] Absorber;
00091   delete[] NeutronSteps;
00092 }
00093 //
00094 // copy constructor
00095 //
00096 ReactorDetector::ReactorDetector(const ReactorDetector& Detector){
00097 
00098   R0=Detector.R0;
00099   R1=Detector.R1;
00100   R2=Detector.R2;
00101   GdConcentration=Detector.GdConcentration;
00102   refracGd=Detector.refracGd;
00103   refracSc=Detector.refracSc;
00104   mfpGd=Detector.mfpGd;
00105   mfpSc=Detector.mfpSc;
00106   MeanNeuDispl = Detector.MeanNeuDispl;
00107   FastNeutronOption = Detector.FastNeutronOption;
00108   PMTDebugOption = Detector.PMTDebugOption;
00109   MergeTime = Detector.MergeTime;
00110   Temperature = Detector.Temperature;
00111   tGd=Detector.tGd;
00112   tSc=Detector.tSc;
00113   GdCaptureFraction=Detector.GdCaptureFraction;
00114   GammasPerGd=Detector.GammasPerGd;
00115   attenlGd=Detector.attenlGd;
00116   attenlSc=Detector.attenlSc;
00117   PhotonsPerMeV=Detector.PhotonsPerMeV;
00118   PMTcoverage=Detector.PMTcoverage;
00119   PMTdiameter=Detector.PMTdiameter;
00120   PMTqe=Detector.PMTqe;
00121   ScintDecayTime=Detector.ScintDecayTime;
00122   for(int i=0;i<4;i++){
00123     a[i] = Detector.a[i];
00124     z[i] = Detector.z[i];
00125     nn[i] = Detector.nn[i];
00126   }
00127 
00128   Npmt = Detector.Npmt;
00129   Ncosbins = Detector.Ncosbins;
00130   Nphibins = Detector.Nphibins;
00131 
00132   Xpmt = new double[Npmt];
00133   Ypmt = new double[Npmt];
00134   Zpmt = new double[Npmt];
00135   PMTOutputArray = new PMToutput*[Npmt];
00136 
00137   double* pold = Detector.Xpmt+Npmt;
00138   double* pnew = Xpmt+Npmt;
00139   while(pnew>Ypmt)*--pnew=*--pold;
00140 
00141   pold = Detector.Ypmt+Npmt;
00142   pnew = Ypmt+Npmt;
00143   while(pnew>Ypmt)*--pnew=*--pold;
00144 
00145   pold = Detector.Zpmt+Npmt;
00146   pnew = Zpmt+Npmt;
00147   while(pnew>Zpmt)*--pnew=*--pold;
00148 
00149   PHelectron = Detector.PHelectron;
00150   PHpositron = Detector.PHpositron;
00151   PHneutron = Detector.PHneutron;
00152   PHthallium = Detector.PHthallium;
00153   PHbismuth = Detector.PHbismuth;
00154   PHcalifornium = Detector.PHcalifornium;
00155   PHmuon = Detector.PHmuon;
00156 
00157   PHeleReco = Detector.PHeleReco;
00158   PHposReco = Detector.PHposReco;
00159   PHneuReco = Detector.PHneuReco;
00160   PHthaReco = Detector.PHthaReco;
00161   PHbisReco = Detector.PHbisReco;
00162   PHcalReco = Detector.PHcalReco;
00163   PHmuoReco = Detector.PHmuoReco;
00164 
00165   TimePositron = Detector.TimePositron;
00166   TimeNeutron = Detector.TimeNeutron;
00167 
00168   for(int i=0;i<3;i++){
00169     DipolePositron[i] = Detector.DipolePositron[i];
00170     DipoleNeutron[i] = Detector.DipoleNeutron[i];
00171   }
00172 
00173   HitPmtsPos = Detector.HitPmtsPos;
00174   HitPmtsNeu = Detector.HitPmtsNeu;
00175   HitPmtsTha = Detector.HitPmtsTha;
00176   HitPmtsBis = Detector.HitPmtsBis;
00177   HitPmtsCal = Detector.HitPmtsCal;
00178   HitPmtsMuo = Detector.HitPmtsMuo;
00179   HitPmtsEle = Detector.HitPmtsEle;
00180 
00181   for(int i=0;i<7;i++){
00182     MeanNpe[i]=Detector.MeanNpe[i];
00183     MeanAtten[i] = Detector.MeanAtten[i];
00184     MeanSfactor[i] = Detector.MeanSfactor[i];
00185   }
00186 
00187   Nneutron = Detector.Nneutron;
00188   Absorber = new double[Nneutron];
00189   pold = Detector.Absorber+Nneutron;
00190   pnew = Absorber+Nneutron;
00191   while(pnew>Absorber)*--pnew=*--pold;
00192 
00193   NeutronSteps = new double[Nneutron];
00194   pold = Detector.NeutronSteps+Nneutron;
00195   pnew = NeutronSteps+Nneutron;
00196   while(pnew>NeutronSteps)*--pnew=*--pold;
00197 
00198   Ngamma = Detector.Ngamma;
00199   Rgamma = new double[Ngamma];
00200   pold = Detector.Rgamma+3*Ngamma;
00201   pnew = Rgamma+3*Ngamma;
00202   while(pnew>Rgamma)*--pnew=*--pold;
00203 
00204   Egamma = new double[Ngamma];
00205   pold = Detector.Egamma+Ngamma;
00206   pnew = Egamma+Ngamma;
00207   while(pnew>Egamma)*--pnew=*--pold;
00208 
00209   Tgamma = new double[Ngamma];
00210   pold = Detector.Tgamma+Ngamma;
00211   pnew = Tgamma+Ngamma;
00212   while(pnew>Tgamma)*--pnew=*--pold;
00213 
00214   Ogamma = new std::string[Ngamma];
00215   std::string* rold = Detector.Ogamma+Ngamma;
00216   std::string* rnew = Ogamma+Ngamma;
00217   while(rnew>Ogamma)*--rnew=*--rold;
00218 
00219   PMTsize = Detector.PMTsize;
00220   PMTlen = Detector.PMTlen;
00221   PMTdata = new PMT[PMTsize];
00222   PMT* qold = Detector.PMTdata+PMTsize;
00223   PMT* qnew = PMTdata+PMTsize;
00224   while(qnew>PMTdata)*--qnew=*--qold;
00225 
00226   PMTelectron = new PMT[Npmt];
00227   qold = Detector.PMTelectron+Npmt;
00228   qnew = PMTelectron+Npmt;
00229   while(qnew>PMTelectron)*--qnew=*--qold;
00230 
00231   PMTpositron = new PMT[Npmt];
00232   qold = Detector.PMTpositron+Npmt;
00233   qnew = PMTpositron+Npmt;
00234   while(qnew>PMTpositron)*--qnew=*--qold;
00235 
00236   PMTneutron = new PMT[Npmt];
00237   qold = Detector.PMTneutron+Npmt;
00238   qnew = PMTneutron+Npmt;
00239   while(qnew>PMTneutron)*--qnew=*--qold;
00240 
00241   PMTthallium = new PMT[Npmt];
00242   qold = Detector.PMTthallium+Npmt;
00243   qnew = PMTthallium+Npmt;
00244   while(qnew>PMTthallium)*--qnew=*--qold;
00245 
00246   PMTbismuth = new PMT[Npmt];
00247   qold = Detector.PMTbismuth+Npmt;
00248   qnew = PMTbismuth+Npmt;
00249   while(qnew>PMTbismuth)*--qnew=*--qold;
00250 
00251   PMTcalifornium = new PMT[Npmt];
00252   qold = Detector.PMTcalifornium+Npmt;
00253   qnew = PMTcalifornium+Npmt;
00254   while(qnew>PMTcalifornium)*--qnew=*--qold;
00255 
00256   PMTmuon = new PMT[Npmt];
00257   qold = Detector.PMTmuon+Npmt;
00258   qnew = PMTmuon+Npmt;
00259   while(qnew>PMTmuon)*--qnew=*--qold;
00260 }    
00261 //
00262 // = operator
00263 //
00264 ReactorDetector& ReactorDetector::operator=(const ReactorDetector& rhs){
00265 
00266   if(this == &rhs)return *this;
00267 
00268   R0=rhs.R0;
00269   R1=rhs.R1;
00270   R2=rhs.R2;
00271   GdConcentration=rhs.GdConcentration;
00272   refracGd=rhs.refracGd;
00273   refracSc=rhs.refracSc;
00274   mfpGd=rhs.mfpGd;
00275   mfpSc=rhs.mfpSc;
00276   MeanNeuDispl = rhs.MeanNeuDispl;
00277   FastNeutronOption = rhs.FastNeutronOption;
00278   PMTDebugOption = rhs.PMTDebugOption;
00279   MergeTime = rhs.MergeTime;
00280   Temperature = rhs.Temperature;
00281   tGd=rhs.tGd;
00282   tSc=rhs.tSc;
00283   GdCaptureFraction=rhs.GdCaptureFraction;
00284   GammasPerGd=rhs.GammasPerGd;
00285   attenlGd=rhs.attenlGd;
00286   attenlSc=rhs.attenlSc;
00287   PhotonsPerMeV=rhs.PhotonsPerMeV;
00288   PMTcoverage=rhs.PMTcoverage;
00289   PMTdiameter=rhs.PMTdiameter;
00290   PMTqe=rhs.PMTqe;
00291   ScintDecayTime=rhs.ScintDecayTime;
00292   for(int i=0;i<4;i++){
00293     a[i] = rhs.a[i];
00294     z[i] = rhs.z[i];
00295     nn[i] = rhs.nn[i];
00296   }
00297 
00298   Npmt = rhs.Npmt;
00299   Ncosbins = rhs.Ncosbins;
00300   Nphibins = rhs.Nphibins;
00301 
00302   delete[] Xpmt;
00303   delete[] Ypmt;
00304   delete[] Zpmt;
00305   delete[] PMTOutputArray;
00306 
00307   Xpmt = new double[Npmt];
00308   Ypmt = new double[Npmt];  
00309   Zpmt = new double[Npmt];
00310   PMTOutputArray = new PMToutput*[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 // PMT sub-class  constructors,destructor,operator
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 // PMT sub-class equal operator
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 // Detector simulation function
00620 //
00621 void ReactorDetector::LightsOut(ReactorEvent& Event, int eventnum){
00622 
00623   PMTlen=0;
00624   delete[] PMTdata;
00625   PMTdata = new PMT[PMTsize];
00626   CleanPMTOutputArray();
00627   //  
00628   // Detector simulation routines.
00629   //
00630   GenerateSecondaries(Event);
00631   GenerateGammas(Event);
00632   GenerateResponse(Event);
00633 
00634   ConsolidatePMT();
00635   if(PMTDebugOption){
00636     if(eventnum == 0){
00637       ofstream file ( "testdump.txt");
00638       if(!file.is_open()){
00639         cout << "File failure, testdump.txt may include old data" << endl;
00640       }
00641       else{
00642         file << "Event \tPMTid \tTime \tEnergy" << endl;
00643         file.close();
00644       }
00645     }
00646     PMTOutputArray[10]->testdump(eventnum);
00647     PMTOutputArray[50]->testdump(eventnum);
00648     PMTOutputArray[100]->testdump(eventnum);
00649     }
00650   /* Reconstruction code.
00651      Position and energy reconstruction only done for positron
00652      and neutron.  Other particle types have ideal position
00653      information from GenerateSecondaries and a simple energy
00654      reconstruction using ideal position information and some
00655      effects such as PMT smearing and attenuation done in
00656      GenerateResponse.
00657   */
00658   CalculateQPosition("P");
00659   CalculateQPosition("N");
00660 
00661   ImproveQPosition("P");
00662   ImproveQPosition("N");
00663 
00664   CalculateEnergy("P");
00665   CalculateEnergy("N");
00666 
00667 }
00668 //
00669 // Access to PMT data
00670 //
00671 void ReactorDetector::GetPMTdata(int npmt,double* XYZpmt,
00672                                  int& nele,double* PHele,double* Tele,
00673                                  int& npos,double* PHpos,double* Tpos,
00674                                  int& nneu,double* PHneu,double* Tneu,
00675                                  int& ntha,double* PHtha,double* Ttha,
00676                                  int& nbis,double* PHbis,double* Tbis,
00677                                  int& ncal,double* PHcal,double* Tcal,
00678                                  int& nmuo,double* PHmuo,double* Tmuo,
00679                                  int* npepos, int* npeneu, int* npetha,
00680                                  int* npebis, int* npecal, int* npemuo,
00681                                  int* npeele){
00682   
00683   *(XYZpmt+0) = *(Xpmt+npmt);
00684   *(XYZpmt+1) = *(Ypmt+npmt);
00685   *(XYZpmt+2) = *(Zpmt+npmt);
00686 
00687   nele=0;
00688   npos=0;
00689   nneu=0;
00690   ntha=0;
00691   nbis=0;
00692   ncal=0;
00693   nmuo=0;
00694 
00695   for(PMT* PMTptr=PMTdata;PMTptr<PMTdata+PMTlen;PMTptr++){
00696     if(npmt==(*PMTptr).npmt){
00697       if((*PMTptr).Parent=="E"){
00698         *(PHele+nele) = (*PMTptr).PHsmr;
00699         *(Tele+nele) = (*PMTptr).TimeSmr;
00700         *(npeele+nele) = (*PMTptr).Npe;
00701         nele++;
00702       }
00703       if((*PMTptr).Parent=="P"){
00704         *(PHpos+npos) = (*PMTptr).PHsmr;
00705         *(Tpos+npos) = (*PMTptr).TimeSmr;
00706         *(npepos+npos) = (*PMTptr).Npe;
00707         npos++;
00708       }
00709       if((*PMTptr).Parent=="N"){
00710         *(PHneu+nneu) = (*PMTptr).PHsmr;
00711         *(Tneu+nneu) = (*PMTptr).TimeSmr;
00712         *(npeneu+nneu) = (*PMTptr).Npe;
00713         nneu++;
00714       }
00715       if((*PMTptr).Parent=="Tl"){
00716         *(PHtha+ntha) = (*PMTptr).PHsmr;
00717         *(Ttha+ntha) = (*PMTptr).TimeSmr;
00718         *(npetha+ntha) = (*PMTptr).Npe;
00719         ntha++;
00720       }
00721       if((*PMTptr).Parent=="Bi"){
00722         *(PHbis+nbis) = (*PMTptr).PHsmr;
00723         *(Tbis+nbis) = (*PMTptr).TimeSmr;
00724         *(npebis+nbis) = (*PMTptr).Npe;
00725         nbis++;
00726       }
00727       if((*PMTptr).Parent=="Cf"){
00728         *(PHcal+ncal) = (*PMTptr).PHsmr;
00729         *(Tcal+ncal) = (*PMTptr).TimeSmr;
00730         *(npecal+ncal) = (*PMTptr).Npe;
00731         ncal++;
00732       }
00733       if((*PMTptr).Parent=="muon"){
00734         *(PHmuo+nmuo) = (*PMTptr).PHsmr;
00735         *(Tmuo+nmuo) = (*PMTptr).TimeSmr;
00736         *(npemuo+nmuo) = (*PMTptr).Npe;
00737         nmuo++;
00738       }
00739     } 
00740   }
00741   return;
00742 }
00743 // Return information about a PMT for a given secondary
00744 void ReactorDetector::GetPMTdata(std::string parent,int pmtindex,
00745                                  double* XYZpmt, int* hits,
00746                                  double* PH, double* Q, double* Time,
00747                                  int* npe){
00748 
00749   *(XYZpmt+0) = *(Xpmt+pmtindex);
00750   *(XYZpmt+1) = *(Ypmt+pmtindex);
00751   *(XYZpmt+2) = *(Zpmt+pmtindex);
00752 
00753 
00754   for(PMT* PMTptr=PMTdata;PMTptr<PMTdata+PMTlen;PMTptr++){
00755     if(pmtindex==(*PMTptr).npmt){
00756       if((*PMTptr).Parent==parent){
00757         *PH += (*PMTptr).PHsmr;
00758         *Time += (*PMTptr).TimeSmr;
00759         *npe += (*PMTptr).Npe;
00760        *Q += (*PMTptr).Q;
00761         ++*hits;
00762       }
00763     }
00764   }
00765   return;
00766 }
00767 
00768 int ReactorDetector::GetHitPMTsPos(){
00769   return HitPmtsPos;
00770 }
00771 int ReactorDetector::GetHitPMTsNeu(){
00772   return HitPmtsNeu;
00773 }
00774 int ReactorDetector::GetHitPMTsTha(){
00775   return HitPmtsTha;
00776 }
00777 int ReactorDetector::GetHitPMTsBis(){
00778   return HitPmtsBis;
00779 }
00780 int ReactorDetector::GetHitPMTsCal(){
00781   return HitPmtsCal;
00782 }
00783 int ReactorDetector::GetHitPMTsMuo(){
00784   return HitPmtsMuo;
00785 }
00786 int ReactorDetector::GetHitPMTsEle(){
00787   return HitPmtsEle;
00788 }
00789 void ReactorDetector::GenerateSecondaries(ReactorEvent& Event){
00790 
00791   //cout << "GenerateSecondaries(Event)" << endl;
00792   //
00793   // Displace the positron vertex by its range along the positron direction.
00794   //
00795   int Npositron = Event.GetNumPositron();
00796   //cout << Npositron << " positron(s):" << endl;
00797 
00798   double* Rpos = new double[Event.Rdim*Npositron];
00799 
00800   // loop over the positrons
00801   for(int n=0;n<Npositron;n++){
00802 
00803     double kepos = Event.GetPositronKE(n);
00804     double tpos = 0.;
00805     double rangepos = 0.;
00806     PositronRange(kepos,rangepos,tpos);
00807     tpos+=Event.GetPositronTime(n);
00808     //
00809     double cospos = Event.GetPositronCosTheta(n);
00810     double phipos = Event.GetPositronPhi(n);
00811 
00812     double rint[3];
00813     for(int i=0;i<3;i++){
00814       rint[i] = *(Event.GetPositronVertexXYZ()+(i+n*Event.Rdim));
00815     }
00816     //
00817     double rpos[3];
00818     rpos[0]=rint[0]+rangepos*sqrt(1-cospos*cospos)*cos(phipos);
00819     rpos[1]=rint[1]+rangepos*sqrt(1-cospos*cospos)*sin(phipos);
00820     rpos[2]=rint[2]+rangepos*cospos;
00821     double rrpos = sqrt(rpos[0]*rpos[0]+rpos[1]*rpos[1]+rpos[2]*rpos[2]);
00822     //
00823     // set positron vertex in Event class
00824     //
00825     Rpos[0+n*Event.Rdim] = tpos;
00826     for(int i=1;i<4;i++)Rpos[i+n*Event.Rdim]=rpos[i-1];
00827     Rpos[4+n*Event.Rdim] = rrpos;
00828     Rpos[5+n*Event.Rdim] = rpos[2]/rrpos;
00829     Rpos[6+n*Event.Rdim] = atan2(rpos[1],rpos[0]);
00830     if(Rpos[6+n*Event.Rdim]<0)Rpos[6+n*Event.Rdim]+=2*PI;
00831     Event.SetPositronVertex(n,Rpos);
00832   }
00833   //  for(int n=0;n<Event.Rdim*Npositron;n++){
00834   //    cout << " Rpos[" << n << "] " << Rpos[n] << endl;
00835   //  }
00836   // done with positrons
00837 
00838   int Nneutron = Event.GetNumNeutron();
00839   //cout << Nneutron << " neutron(s):" << endl;
00840 
00841   delete[] Absorber;
00842   delete[] NeutronSteps;
00843 
00844   Absorber = new double[Nneutron];
00845   NeutronSteps = new double[Nneutron];
00846 
00847   double* Rneu = new double[Event.Rdim*Nneutron];
00848 
00849   // loop over the neutrons
00850   for(int n=0;n<Nneutron;n++){
00851     //
00852     double rneu[3];
00853     double tneu = 0.;
00854     if(FastNeutronOption){
00855       //
00856       // Throw the neutron vertex via a simple Gaussian
00857       //
00858       float rn[3];
00859       int len=3;
00860       rnormx_(rn,len,ranlux_);
00861       //
00862       // Random for throwing neutron time
00863       //
00864       len=1;
00865       float rv[1];
00866       ranlux_(rv,len);
00867       //
00868       //  Simple scheme that is not right at boundaries
00869       //  Neglects discontinuity in mean free path;
00870       //  Time will be taken as proportional to sqrt(path)
00871       //
00872       double rint[3];
00873       for(int i=0;i<3;i++){
00874         rint[i] = *(Event.GetNeutronVertexXYZ()+(i+n*Event.Rdim));
00875       }
00876       //
00877       double rrint = sqrt(rint[0]*rint[0]+rint[1]*rint[1]+rint[2]*rint[2]);
00878       double path = 0.;
00879       if(rrint<=R0){
00880         for (int i=0;i<3;i++)rneu[i]=rint[i]+mfpGd*rn[i]/sqrt(3.);
00881         path = mfpGd*sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
00882         tneu = tGd*sqrt(path/mfpGd)*log(1/rv[0]);
00883       }
00884       else if(rrint<R2){
00885         for (int i=0;i<3;i++)rneu[i]=rint[i]+mfpSc*rn[i]/sqrt(3.);
00886         path = mfpSc*sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
00887         tneu = tSc*sqrt(path/mfpSc)*log(1/rv[0]);
00888       }
00889       else{
00890         for (int i=0;i<3;i++)rneu[i]=999999;
00891         tneu = 999999.;
00892       }
00893       //
00894       // Shift z displacement by mean Chooz value
00895       //
00896       rneu[2]+=MeanNeuDispl;
00897     }
00898     else{
00899       double pneu[3];
00900       for(int i=0;i<3;i++){
00901         pneu[i]=*(Event.GetNeutronPxyz()+(i+n*Event.Pdim));
00902         //cout << "pneu[" << i << "] " << pneu[i] << endl;
00903       }
00904       for(int i=0;i<3;i++){
00905         rneu[i]=*(Event.GetNeutronVertexXYZ()+(i+n*Event.Rdim));
00906         //cout << "rneu[" << i << "] " << rneu[i] << endl;
00907       }
00908       double keneu = Event.GetNeutronKE(n);
00909       NeutronSteps[n]=0.;
00910       Absorber[n]=0.;
00911       NeutronPropagator(R0,R2,keneu,Absorber[n],NeutronSteps[n],
00912                         pneu,rneu,tneu,Temperature,a,z,nn);
00913       //cout << " Absorber[" << n << "] " << Absorber[n] << endl;
00914       //cout << " NeutronSteps[" << n << "] " << NeutronSteps[n] << endl;
00915       tneu+=Event.GetNeutronTime(n);
00916       //cout << " tneu " << tneu << endl;
00917     }
00918     //
00919     // set neutron vertex in Event class
00920     //
00921     Rneu[0+n*Event.Rdim] = tneu;
00922     for(int i=1;i<4;i++)Rneu[i+n*Event.Rdim]=rneu[i-1];
00923     double rrneu = sqrt(rneu[0]*rneu[0]+rneu[1]*rneu[1]+rneu[2]*rneu[2]);
00924     Rneu[4+n*Event.Rdim] = rrneu;
00925     Rneu[5+n*Event.Rdim] = rneu[2]/rrneu;
00926     Rneu[6+n*Event.Rdim] = atan2(rneu[1],rneu[0]);
00927     if(Rneu[6+n*Event.Rdim]<0)Rneu[6+n*Event.Rdim]+=2*PI;
00928     Event.SetNeutronVertex(n,Rneu);
00929   }
00930   //  for(int n=0;n<Event.Rdim*Nneutron;n++){
00931   //    cout << " Rneu[" << n << "] " << Rneu[n] << endl;
00932   //  }
00933   // done with neutrons
00934   delete[] Rpos;
00935   delete[] Rneu;
00936 }
00937 
00938 void ReactorDetector::GenerateResponse(ReactorEvent& Event){
00939 
00940   //cout << "GenerateResponse(Event)" << endl;
00941   //
00942   PHelectron = 0;
00943   PHpositron = 0;
00944   PHneutron = 0;
00945   PHthallium = 0;
00946   PHbismuth = 0;
00947   PHcalifornium = 0;
00948   PHmuon = 0;
00949 
00950   PHeleReco = 0;
00951   PHposReco = 0;
00952   PHneuReco = 0;
00953   PHthaReco = 0;
00954   PHbisReco = 0;
00955   PHcalReco = 0;
00956   PHmuoReco = 0;
00957 
00958   HitPmtsPos = 0;
00959   HitPmtsNeu = 0;
00960   HitPmtsTha = 0;
00961   HitPmtsBis = 0;
00962   HitPmtsCal = 0;
00963   HitPmtsMuo = 0;
00964   HitPmtsEle = 0;
00965 
00966   for(int i=0;i<7;i++){
00967     MeanSfactor[i] = 0;
00968     MeanAtten[i] = 0;
00969     MeanNpe[i] = 0;
00970     MeanQ[i] = 0;
00971   }
00972   //
00973   //
00974   //
00975   int pmtindex=0;
00976   double* xpmt=Xpmt;
00977   double* ypmt=Ypmt;
00978   double* zpmt=Zpmt;
00979   
00980   PMT* PMTptr = PMTdata;
00981 
00982   //cout << Ngamma << " gammas" << endl;
00983   for(int i=0;i<Ncosbins;i++){
00984     for(int j=0;j<Nphibins;j++){
00985       double rpmt[3];
00986       rpmt[0] = *xpmt++;
00987       rpmt[1] = *ypmt++;
00988       rpmt[2] = *zpmt++;
00989       //for (int l=0;l<3;l++) cout << "  rpmt[" << l << "] " << rpmt[l] << endl;
00990       double r[3];
00991       double dr[3];
00992       double tr[3];
00993       //      
00994       double sumele[7]={0,0,0,0,0,0,0};
00995       double sumpos[7]={0,0,0,0,0,0,0};
00996       double sumneu[7]={0,0,0,0,0,0,0};
00997       double sumtha[7]={0,0,0,0,0,0,0};
00998       double sumbis[7]={0,0,0,0,0,0,0};
00999       double sumcal[7]={0,0,0,0,0,0,0};
01000       double summuo[7]={0,0,0,0,0,0,0};
01001       //
01002       double* rgptr=Rgamma;
01003       double* egptr=Egamma;
01004       double* tgptr=Tgamma;
01005       std::string* ogptr=Ogamma;
01006       //
01007       int nelegt0=0;
01008       int nposgt0=0;
01009       int nneugt0=0;
01010       int nthagt0=0;
01011       int nbisgt0=0;
01012       int ncalgt0=0;
01013       int nmuogt0=0;
01014       //
01015       for(int k=0;k<Ngamma;k++){
01016         //
01017         for (int l=0;l<3;l++)r[l]=*rgptr++;
01018         double rdotr = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
01019         //cout << i << " " << j << " " << k << "  r " << sqrt(rdotr) << endl;
01020         //
01021         // nothing to do if r>R1: cycle through useless info
01022         //
01023         if(sqrt(rdotr)>R1){
01024           double dummy = *egptr++;
01025           dummy = *tgptr++;
01026           std::string sdummy = *ogptr++;
01027           continue;
01028         }
01029         //
01030         for (int l=0;l<3;l++)dr[l]=rpmt[l]-r[l];
01031         double ddr = sqrt(dr[0]*dr[0]+dr[1]*dr[1]+dr[2]*dr[2]);
01032         //cout << "  ddr " << ddr;
01033         for (int l=0;l<3;l++)tr[l]=dr[l]/ddr;
01034         double tdotr = r[0]*tr[0]+r[1]*tr[1]+r[2]*tr[2];
01035         //cout << "  tdotr " << tdotr;
01036         //
01037         // solve for intersection with inner radius
01038         //
01039         double sGd;
01040         double sSc;
01041         //
01042         // check for no intersections
01043         //
01044         double test = R0*R0+tdotr*tdotr-rdotr;
01045         bool inter = (test>=0);
01046         //
01047         // A straight line from the interaction vertex to the pmt
01048         // goes through the inner radius at least once.  In this case,
01049         // the light will see Gd-scint. and unloaded scint.
01050         //
01051         //cout << "  " << inter << endl;
01052         if(inter){
01053           double sint1 = -tdotr+sqrt(test);
01054           double sint2 = -tdotr-sqrt(test);
01055           if(sint1>0&&sint2<0){     // 1 intersection at soln 1
01056             sGd = sint1;
01057           }
01058           else if(sint2>0&&sint1<0){     // 1 intersection at soln 2
01059             sGd = sint2;
01060           }
01061           else if(sint2>0&&sint1>0){     // 2 intersections 
01062             sGd = sint1-sint2;
01063             if(sGd<0)sGd*=-1;
01064           }
01065           else{                          // no intersections towards PMT
01066             sGd = 0;
01067           }
01068           sSc = ddr-sGd;
01069         }
01070         //
01071         //  otherwise the light only goes through unloaded scintillator
01072         //
01073         else{
01074           sGd = 0;
01075           sSc = ddr;
01076         }
01077         //cout << "       sGd " << sGd << " sSc " << sSc << endl;
01078         //
01079         //  increment PH-- make atten. relative to center
01080         //
01081         double dsGd = sGd-R0;
01082         double dsSc = sSc-R2+R0;
01083         double atten0 = exp(-sGd/attenlGd-sSc/attenlSc);     
01084         double atten = exp(-dsGd/attenlGd-dsSc/attenlSc);     
01085         //
01086         // solid angle factor
01087         //
01088         double sfactor=0;
01089         for (int l=0;l<3;l++){sfactor+=tr[l]*rpmt[l]/R2;}
01090         sfactor *= R2*R2/ddr/ddr;
01091         //
01092         // photostatistics
01093         //
01094         double phgen = *egptr++;
01095         //cout << "phgen " << phgen << endl;
01096         //
01097         //
01098         // Now throw number of photons as gaussian
01099         // around mean of phgen*PhotonsPerMeV
01100         //
01101         int len = 1;
01102         float photons[1];
01103         rnormx_(photons,len,ranlux_);
01104         photons[0]*=sqrt(phgen*PhotonsPerMeV);
01105         float xpe = (phgen*PhotonsPerMeV) + (double) photons[0];
01106         xpe *= (atten0*sfactor);
01107         xpe *= (PMTcoverage/Npmt*PMTqe);
01108         //
01109         int npe;
01110         int error;
01111         rnpssn_(xpe,npe,error);
01112         //cout << "npe " << npe << endl;
01113         //
01114         float Q = 0.;
01115         Q = GenerateQ(npe);
01116         double phsmr = Q/PMTqe/PMTcoverage/PhotonsPerMeV;   
01117         //cout << "phsmr " << phsmr << endl;
01118         //        
01119         double tgen = *tgptr++;
01120         //cout << "tgen " << tgen << endl;
01121         double tsmr = refracGd*sGd/c+refracSc*sSc/c + tgen; 
01122         //cout << "tsmr " << tsmr << endl;
01123         //
01124         std::string type = *ogptr++;
01125         //cout << "type " << type << endl;
01126         if(type == "E"){
01127           nelegt0++;
01128           sumele[0]+=phsmr;
01129           sumele[1]+=phgen;
01130           sumele[2]+=phgen*atten0*sfactor;
01131           sumele[3]+=npe;
01132           sumele[4]+=sfactor;
01133           sumele[5]+=atten;
01134           sumele[6]+=Q;
01135           //
01136           PHelectron+=phsmr;
01137         }
01138         else if(type == "P"){
01139           nposgt0++;
01140           sumpos[0]+=phsmr;
01141           sumpos[1]+=phgen;
01142           sumpos[2]+=phgen*atten0*sfactor;
01143           sumpos[3]+=npe;
01144           sumpos[4]+=sfactor;
01145           sumpos[5]+=atten;
01146           sumpos[6]+=Q;
01147           //
01148           PHpositron+=phsmr;
01149         }
01150         else if(type == "N"){
01151           nneugt0++;
01152           sumneu[0]+=phsmr;
01153           sumneu[1]+=phgen;
01154           sumneu[2]+=phgen*atten0*sfactor;
01155           sumneu[3]+=npe;
01156           sumneu[4]+=sfactor;
01157           sumneu[5]+=atten;
01158           sumneu[6]+=Q;
01159           //      
01160           PHneutron+=phsmr;
01161         }
01162         else if(type == "Tl"){
01163           nthagt0++;
01164           sumtha[0]+=phsmr;
01165           sumtha[1]+=phgen;
01166           sumtha[2]+=phgen*atten0*sfactor;
01167           sumtha[3]+=npe;
01168           sumtha[4]+=sfactor;
01169           sumtha[5]+=atten;
01170           sumtha[6]+=Q;
01171           //      
01172           PHthallium+=phsmr;
01173         }
01174         else if(type == "Bi"){
01175           nbisgt0++;
01176           sumbis[0]+=phsmr;
01177           sumbis[1]+=phgen;
01178           sumbis[2]+=phgen*atten0*sfactor;
01179           sumbis[3]+=npe;
01180           sumbis[4]+=sfactor;
01181           sumbis[5]+=atten;
01182           sumbis[6]+=Q;
01183           //      
01184           PHbismuth+=phsmr;
01185         }
01186         else if(type == "Cf"){
01187           ncalgt0++;
01188           sumcal[0]+=phsmr;
01189           sumcal[1]+=phgen;
01190           sumcal[2]+=phgen*atten0*sfactor;
01191           sumcal[3]+=npe;
01192           sumcal[4]+=sfactor;
01193           sumcal[5]+=atten;
01194           sumcal[6]+=Q;
01195           //      
01196           PHcalifornium+=phsmr;
01197         }
01198         else if(type == "muon"){
01199           nmuogt0++;
01200           summuo[0]+=phsmr;
01201           summuo[1]+=phgen;
01202           summuo[2]+=phgen*atten0*sfactor;
01203           summuo[3]+=npe;
01204           summuo[4]+=sfactor;
01205           summuo[5]+=atten;
01206           summuo[6]+=Q;
01207           //      
01208           PHmuon+=phsmr;
01209         }
01210         else{
01211           cout << "  ERROR: GenerateResponse has unrecognized type " 
01212                << type << endl;
01213         }
01214         //        
01215         PMT PMTtmp;
01216         PMTtmp.npmt = pmtindex;
01217         for(int l=0;l<3;l++)PMTtmp.xyz[l] = rpmt[l];
01218         PMTtmp.PHgen = phgen;
01219         PMTtmp.PHsmr = phsmr;
01220         PMTtmp.PathGd = sGd;
01221         PMTtmp.PathSc = sSc;
01222         PMTtmp.Atten = atten;
01223         PMTtmp.Sfactor = sfactor;
01224         PMTtmp.TimeGen = tgen;
01225         PMTtmp.TimeSmr = tsmr;
01226         PMTtmp.Xpe = xpe;
01227         PMTtmp.Npe = npe;
01228         PMTtmp.Q = Q;
01229         PMTtmp.Parent = type;
01230         
01231         if(PMTptr<PMTdata+PMTsize){
01232           *PMTptr++=PMTtmp;
01233           PMTlen++;
01234         }
01235         else{
01236           //
01237           // expand the data array first
01238           //
01239           PMTptr = ResizePMT();
01240           *PMTptr++=PMTtmp;
01241           PMTlen++;
01242         }
01243       }
01244       //cout << "PHpositron " << PHpositron << endl;
01245       //cout << "PHneutron  " << PHneutron << endl;
01246       //
01247       //  Crude estimate of reconstructed PH
01248       //  Includes poisson npe smearing and fluctuation in atten length
01249       //
01250       /*  PHposReco, PHneuReco, etc., are recalculated in CalculateEnergy 
01251           in order to use reconstructed position for energy 
01252           reconstruction.  These calculcations are done with ideal
01253           knowledge of position.
01254       */
01255       if(sumpos[2]>0){
01256         PHposReco+=sumpos[0]*sumpos[1]/sumpos[2];
01257         MeanNpe[0]+=sumpos[3]/Npmt;
01258         MeanQ[0]+=sumpos[6]/Npmt;
01259         MeanSfactor[0]+=sumpos[4]/Npmt/nposgt0;
01260         MeanAtten[0]+=sumpos[5]/Npmt/nposgt0;
01261       }
01262       if(sumneu[2]>0){
01263         PHneuReco+=sumneu[0]*sumneu[1]/sumneu[2];
01264         MeanNpe[1]+=sumneu[3]/Npmt;
01265         MeanQ[1]+=sumneu[6]/Npmt;
01266         MeanSfactor[1]+=sumneu[4]/Npmt/nneugt0;
01267         MeanAtten[1]+=sumneu[5]/Npmt/nneugt0;
01268       }
01269       if(sumtha[2]>0){
01270         PHthaReco+=sumtha[0]*sumtha[1]/sumtha[2];
01271         MeanNpe[2]+=sumtha[3]/Npmt;
01272         MeanQ[2]+=sumtha[6]/Npmt;
01273         MeanSfactor[2]+=sumtha[4]/Npmt/nthagt0;
01274         MeanAtten[2]+=sumtha[5]/Npmt/nthagt0;
01275       }
01276       if(sumbis[2]>0){
01277         PHbisReco+=sumbis[0]*sumbis[1]/sumbis[2];
01278         MeanNpe[3]+=sumbis[3]/Npmt;
01279         MeanQ[3]+=sumbis[6]/Npmt;
01280         MeanSfactor[3]+=sumbis[4]/Npmt/nbisgt0;
01281         MeanAtten[3]+=sumbis[5]/Npmt/nbisgt0;
01282       }
01283       if(sumcal[2]>0){
01284         PHcalReco+=sumcal[0]*sumcal[1]/sumcal[2];
01285         MeanNpe[4]+=sumcal[3]/Npmt;
01286         MeanQ[4]+=sumcal[6]/Npmt;
01287         MeanSfactor[4]+=sumcal[4]/Npmt/ncalgt0;
01288         MeanAtten[4]+=sumcal[5]/Npmt/ncalgt0;
01289       }
01290       if(summuo[2]>0){
01291         PHmuoReco+=summuo[0]*summuo[1]/summuo[2];
01292         MeanNpe[5]+=summuo[3]/Npmt;
01293         MeanQ[5]+=summuo[6]/Npmt;
01294         MeanSfactor[5]+=summuo[4]/Npmt/nmuogt0;
01295         MeanAtten[5]+=summuo[5]/Npmt/nmuogt0;
01296       }
01297       if(sumele[2]>0){
01298         PHeleReco+=sumele[0]*sumele[1]/sumele[2];
01299         MeanNpe[6]+=sumele[3]/Npmt;
01300         MeanQ[6]+=sumele[6]/Npmt;
01301         MeanSfactor[6]+=sumele[4]/Npmt/nelegt0;
01302         MeanAtten[6]+=sumele[5]/Npmt/nelegt0;
01303       }
01304       if(sumpos[6]>0){
01305         HitPmtsPos++;
01306       }
01307       if(sumneu[6]>0){
01308         HitPmtsNeu++;
01309       }
01310       if(sumtha[6]>0){
01311         HitPmtsTha++;
01312       }
01313       if(sumbis[6]>0){
01314         HitPmtsBis++;
01315       }
01316       if(sumcal[6]>0){
01317         HitPmtsCal++;
01318       }
01319       if(summuo[6]>0){
01320         HitPmtsMuo++;
01321       }
01322       if(sumele[6]>0){
01323         HitPmtsEle++;
01324       }
01325       //
01326       // increment PMT index
01327       //
01328       pmtindex++;
01329     }
01330   }
01331 }
01332 
01333 void ReactorDetector::GenerateGammas(ReactorEvent& Event){
01334 
01335   //cout << "GenerateGammas(Event)" << endl;
01336   //
01337   // first check for any gammas produced in the event itself
01338   int ngamma = Event.GetNumGamma();
01339   //cout << " " << ngamma << " gammas from the Event" << endl;
01340   //
01341   int maxgamma = 100+ngamma;
01342   double* egamma = new double[maxgamma];
01343   double* tgamma = new double[maxgamma];
01344   double* rgamma = new double[3*maxgamma];
01345   std::string* ogamma = new std::string[maxgamma];
01346   //
01347   double* egptr=egamma;
01348   double* rgptr=rgamma;
01349   double* tgptr=tgamma;
01350   std::string* ogptr=ogamma;
01351   //
01352   // event gammas
01353   //
01354   for(int n=0;n<ngamma;n++){
01355 
01356     *egptr++ = Event.GetGammaKE(n);
01357     //
01358     for(int i=0;i<3;i++) 
01359       *rgptr++ = *(Event.GetGammaVertexXYZ()+(i+n*Event.Rdim));
01360     //
01361     *tgptr++ = Event.GetGammaTime(n);
01362     //
01363     if(Event.GetGammaParentProcess(n) == 1){
01364       cout << "ERROR: gammas are not produced by inverse-beta decay" << endl;
01365     }
01366     else if(Event.GetGammaParentProcess(n) == 2){
01367       *ogptr++ = "Cf";
01368     }
01369     else if(Event.GetGammaParentProcess(n) == 3){
01370       *ogptr++ = "muon";
01371     }
01372     else if(Event.GetGammaParentProcess(n) == 4){
01373       *ogptr++ = "muon";
01374     }
01375     else{
01376       cout << "Process Id " << Event.GetGammaParentProcess(n) << " is unknown" 
01377            << endl;
01378     }
01379   }
01380   //
01381   // radioactive source gammas
01382   //
01383   if(Event.GetPmtRad() == "on"){
01384 
01385     // pick a PMT location for the decay
01386     const int len = 2;
01387     float r[len];
01388     ranlux_(r,len);
01389     int pmtId = int(r[0]*Npmt);
01390 
01391     // pick decay material and isotope 
01392     std::string Stuff;
01393     int Isotope;
01394 
01395     // this assumes even fractions of decays from all isotopes
01396     //
01397     if(r[1] < 0.25){
01398       Stuff = "thallium"; Isotope = 208; // thallium 208
01399     }
01400     else if(0.25 <= r[1] && r[1] < 0.5){
01401       Stuff = "thallium"; Isotope = 210; // thallium 210
01402     }
01403     else if(0.5 <= r[1] && r[1] < 0.75){
01404       Stuff = "bismuth"; Isotope = 212; // bismuth 212
01405     }
01406     else{
01407       Stuff = "bismuth"; Isotope = 214; // bismuth 214
01408     }
01409     //cout << Stuff << "(" << Isotope << ")" << endl;
01410     //
01411     // Thallium
01412     //
01413     if(Stuff == "thallium"){
01414       MyTlSource Thallium(Isotope);
01415 
01416       //cout << " " << Thallium.GetNumTlGammas() << " gammas" << endl;
01417 
01418       // loop over the decay gammas
01419       for(int n=0;n<Thallium.GetNumTlGammas();n++){
01420 
01421         // pick a direction for each gamma
01422         const int len=2;
01423         float rv[len];
01424         ranlux_(rv,len);
01425         double cosg = -1+2*rv[0];
01426         double phig = 2*PI*rv[1];
01427         double t[3];
01428         t[0] = cos(phig)*sqrt(1-cosg*cosg);
01429         t[1] = sin(phig)*sqrt(1-cosg*cosg);
01430         t[2] = cosg;
01431         //for(int i=0;i<3;i++) cout << "t["<<i<<"] " << t[i] << " ";
01432         //cout << endl;
01433 
01434         // get the energy, time, r and p for each Tl gamma
01435         double eg = Thallium.GetTlGammaEnergy(n); 
01436         double tg = 0.;
01437         //cout << " energy " << eg << " time " << tg << endl;
01438         //
01439         double rg[3];
01440         double* rgp = rg;
01441         rg[0] = Xpmt[pmtId];
01442         rg[1] = Ypmt[pmtId];
01443         rg[2] = Zpmt[pmtId];
01444         //for(int i=0;i<3;i++) cout << " rg["<<i<<"] " << rg[i] << " ";
01445         //cout << endl;
01446         //
01447         // move the "PMT" in 30 cm (for R5912 8" PMTs) from the outer vessel
01448         //
01449         rgp = TranslateRadialDistance(Xpmt[pmtId],Ypmt[pmtId],
01450                                       Zpmt[pmtId],-30.);
01451         for(int i=0;i<3;i++) rg[i] = *rgp++;
01452         //for(int i=0;i<3;i++) cout << " rg["<<i<<"] " << rg[i] << " ";
01453         //cout << endl;
01454         //
01455         double pg[3];
01456         for(int i=0;i<3;i++) pg[i] = eg*t[i];
01457         //for(int i=0;i<3;i++) cout << " pg["<<i<<"] " << pg[i] << " ";
01458         //cout << endl;
01459         //
01460         // compton scatter the gammas through the oil
01461         //
01462         ReactorScint Scint(eg,tg,rg,pg,R0,R1,R2);
01463         //
01464         int scatters = Scint.GetPoints();
01465         for(int s=0;s<scatters;s++){
01466           *egptr++=Scint.GetNextEnergy();
01467           *tgptr++=Scint.GetNextTime();
01468           for(int j=0;j<3;j++) *rgptr++=*(Scint.GetNextXYZ()+s);
01469           *ogptr++ = "Tl";
01470           ngamma++;
01471         }
01472       }
01473       bool inscint = false;
01474       double rdotr = 0.;
01475       for(int n=0;n<ngamma;n++){
01476         rdotr = rgamma[3*n]*rgamma[3*n] +
01477           rgamma[3*n+1]*rgamma[3*n+1] + 
01478           rgamma[3*n+2]*rgamma[3*n+2];
01479         if(sqrt(rdotr)<R1) inscint = true;
01480       }
01481       if(inscint){
01482         pmtRadCount++;
01483         cout << "PMTrad Tl " << pmtRadCount << " radius " << sqrt(rdotr) 
01484              << endl;
01485       }
01486     }
01487     //
01488     // Bismuth
01489     //
01490     else if(Stuff == "bismuth"){
01491       MyBiSource Bismuth(Isotope);
01492 
01493       //cout << " " << Bismuth.GetNumBiGammas() << " gammas" << endl;
01494 
01495       // loop over the decay gammas
01496       for(int n=0;n<Bismuth.GetNumBiGammas();n++){
01497 
01498         // pick a direction for each gamma
01499         const int len=2;
01500         float rv[len];
01501         ranlux_(rv,len);
01502         double cosg = -1+2*rv[0];
01503         double phig = 2*PI*rv[1];
01504         double t[3];
01505         t[0] = cos(phig)*sqrt(1-cosg*cosg);
01506         t[1] = sin(phig)*sqrt(1-cosg*cosg);
01507         t[2] = cosg;
01508         //for(int i=0;i<3;i++) cout << "t["<<i<<"] " << t[i] << " ";
01509         //cout << endl;
01510 
01511         // get the energy, time, r and p for each Bi gamma
01512         double eg = Bismuth.GetBiGammaEnergy(n); 
01513         double tg = 0.;
01514         //cout << " energy " << eg << " time " << tg << endl;
01515         //
01516         double rg[3];
01517         double* rgp = rg;
01518         rg[0] = Xpmt[pmtId];
01519         rg[1] = Ypmt[pmtId];
01520         rg[2] = Zpmt[pmtId];
01521         //for(int i=0;i<3;i++) cout << " rg["<<i<<"] " << rg[i] << " ";
01522         //cout << endl;
01523         //
01524         // move the "PMT" in 30 cm (for R5912 8" PMTs) from the outer vessel
01525         //
01526         rgp = TranslateRadialDistance(Xpmt[pmtId],Ypmt[pmtId],
01527                                      Zpmt[pmtId],-30.);
01528         for(int i=0;i<3;i++) rg[i] = *rgp++;
01529         //for(int i=0;i<3;i++) cout << " rg["<<i<<"] " << rg[i] << " ";
01530         //cout << endl;
01531         //
01532         double pg[3];
01533         for(int i=0;i<3;i++) pg[i] = eg*t[i];
01534         //for(int i=0;i<3;i++) cout << " pg["<<i<<"] " << pg[i] << " ";
01535         //cout << endl;
01536         //
01537         // compton scatter the gammas through the oil
01538         //
01539         ReactorScint Scint(eg,tg,rg,pg,R0,R1,R2);
01540         //
01541         int scatters = Scint.GetPoints();
01542         for(int s=0;s<scatters;s++){
01543           *egptr++=Scint.GetNextEnergy();
01544           *tgptr++=Scint.GetNextTime();
01545           for(int j=0;j<3;j++) *rgptr++=*(Scint.GetNextXYZ()+s);
01546           *ogptr++ = "Bi";
01547           ngamma++;
01548         }
01549       }
01550       bool inscint = false;
01551       double rdotr = 0.;
01552       for(int n=0;n<ngamma;n++){
01553         rdotr = rgamma[3*n]*rgamma[3*n] + 
01554           rgamma[3*n+1]*rgamma[3*n+1] + 
01555           rgamma[3*n+2]*rgamma[3*n+2];
01556         if(sqrt(rdotr)<R1) inscint = true;
01557       }
01558       if(inscint){
01559         pmtRadCount++;
01560         cout << "PMTrad Bi " << pmtRadCount << " radius " << sqrt(rdotr) 
01561              << endl;
01562       }
01563     } 
01564   }
01565   //
01566   // electron gammas, needs improvement
01567   //
01568   int Nelectron = Event.GetNumElectron();
01569   //cout << Nelectron << " electron(s):" << endl;
01570 
01571   // loop over the electrons
01572   for(int n=0;n<Nelectron;n++){
01573 
01574     double rele[3];
01575     for(int i=0;i<3;i++){
01576       rele[i] = *(Event.GetElectronVertexXYZ()+(i+n*Event.Rdim));
01577       //cout << " rele[" << i << "] " << rele[i] << endl;
01578     }
01579     double tele = Event.GetElectronTime(n);
01580     //cout << " tele " << tele << endl;
01581     double keele = Event.GetElectronKE(n);
01582     //cout << " keele " << keele << endl;
01583     //
01584     // generate a "fake" gamma with E=electron KE generated at vertex.
01585     //
01586     for(int i=0;i<3;i++)*rgptr++=rele[i];
01587     *tgptr++ = tele;
01588     *egptr++ = keele;
01589     *ogptr++ = "E";
01590     ngamma++;
01591   } // done with electrons
01592   //cout << "  " << ngamma << " gammas, so far" << endl;
01593   //
01594   // positron gammas
01595   //
01596   int Npositron = Event.GetNumPositron();
01597   //cout << Npositron << " positron(s):" << endl;
01598 
01599   // loop over the positrons
01600   for(int n=0;n<Npositron;n++){
01601 
01602     double rpos[3];
01603     for(int i=0;i<3;i++){
01604       rpos[i] = *(Event.GetPositronVertexXYZ()+(i+n*Event.Rdim));
01605       //cout << " rpos[" << i << "] " << rpos[i] << endl;
01606     }
01607     double ppos[3];
01608     for(int i=0;i<3;i++){
01609       ppos[i] = *(Event.GetPositronPxyz()+(i+n*Event.Pdim));
01610       //cout << " ppos[" << i << "] " << ppos[i] << endl;
01611     }
01612     double tpos = Event.GetPositronTime(n);
01613     //cout << " tpos " << tpos << endl;
01614     double kepos = Event.GetPositronKE(n);
01615     //cout << " kepos " << kepos << endl;
01616     //
01617     // the first is a "fake" gamma with E=positron KE generated at positron
01618     // vertex.
01619     //
01620     for(int i=0;i<3;i++)*rgptr++=rpos[i];
01621     *tgptr++ = tpos;
01622     *egptr++ = kepos;
01623     *ogptr++ = "P";
01624     ngamma++;
01625     //
01626     // Now the two annihilation gammas
01627     //
01628     int len=4;
01629     float rv[4];
01630     ranlux_(rv,len);
01631     double cosg = -1+2*rv[0];
01632     double phig = 2*PI*rv[1];
01633     double t[3];
01634     t[0] = cos(phig)*sqrt(1-cosg*cosg);
01635     t[1] = sin(phig)*sqrt(1-cosg*cosg);
01636     t[2] = cosg;
01637     for (int i=0;i<2;i++){
01638       //
01639       // handle multiple compton scatters
01640       //
01641       ReactorScint Scint(MELECTRON,tpos,rpos,ppos,R0,R1,R2);
01642       //
01643       int scatters = Scint.GetPoints();
01644       for(int s=0;s<scatters;s++){
01645         *egptr++=Scint.GetNextEnergy();
01646         *tgptr++=Scint.GetNextTime();
01647         for(int j=0;j<3;j++) *rgptr++=*(Scint.GetNextXYZ()+s);
01648         *ogptr++ = "P";
01649         ngamma++;
01650       }
01651     }
01652   } // done with positrons
01653   //cout << "  " << ngamma << " gammas, so far" << endl;
01654   //
01655   // Generate neutron gammas
01656   //
01657   Nneutron = Event.GetNumNeutron();
01658   //cout << Nneutron << " neutron(s):" << endl;
01659 
01660   // loop over the neutrons
01661   for(int n=0;n<Nneutron;n++){
01662 
01663     double rneu[3];
01664     for(int i=0;i<3;i++){
01665       rneu[i] = *(Event.GetNeutronVertexXYZ()+(i+n*Event.Rdim));
01666       //cout << " rneu[" << i << "] " << rneu[i] << endl;
01667     }
01668     double tneu = Event.GetNeutronTime(n);
01669     //cout << " tneu " << tneu << endl;
01670     //
01671     // pick H or Gd, then which Gd peak
01672     //
01673     double emax=Hpeak;
01674     if(FastNeutronOption){
01675       //cout << "in FastNeutronOption" << endl;
01676       Absorber[n]=1;
01677       double rrneu = sqrt(rneu[0]*rneu[0]+rneu[1]*rneu[1]+rneu[2]*rneu[2]);
01678       if(rrneu<R0){//in Gd region
01679         int len=2;
01680         float rv[2];
01681         ranlux_(rv,len);
01682         if(rv[0]<=GdCaptureFraction){
01683           //
01684           // model Gd spectrum by Gd155 + Gd157
01685           //
01686           if(rv[1]<=Gd155frac){
01687             emax = Gd155peak;
01688             Absorber[n]=155;
01689           }
01690           else{
01691             emax = Gd157peak;
01692             Absorber[n]=157;
01693           }
01694         }
01695       }
01696     }
01697     else{
01698       //cout << "Absorber[" << n << "] " << Absorber[n] << endl;
01699       if(Absorber[n]==1)emax=Hpeak;
01700       if(Absorber[n]==12)emax=Cpeak;
01701       if(Absorber[n]==155)emax=Gd155peak;
01702       if(Absorber[n]==157)emax=Gd157peak;
01703     }
01704     //
01705     // One photon from an H or C capture
01706     //
01707     if(emax==Hpeak || emax==Cpeak){
01708       //cout << " H or C gammas" << endl;
01709       int len=3;
01710       float rv[3];
01711       ranlux_(rv,len);
01712       std::string material = "scintillator";
01713       double range = GammaComptonRange(material,emax);
01714       range*=log(1/rv[0]);
01715       double cosg = -1+2*rv[1];
01716       double phig = 2*PI*rv[2];
01717       double t[3];
01718       t[0] = cos(phig)*sqrt(1-cosg*cosg);
01719       t[1] = sin(phig)*sqrt(1-cosg*cosg);
01720       t[2] = cosg;
01721       for(int j=0;j<3;j++)*rgptr++ = rneu[j]+range*t[j];
01722       *egptr++ = emax;
01723       *tgptr++ = tneu+range/c*refracSc;
01724       *ogptr++ = "N";
01725       ngamma++;
01726       //cout << "  " << ngamma << " gammas" << endl;
01727     }
01728     //
01729     // otherwise model Gd photons
01730     //
01731     else{
01732       //cout << " Gd gammas" << endl;
01733       int ngd=0;
01734       double* pgd = new double[3*maxgamma];
01735       MyGdDecayModel(emax,ngd,pgd);
01736       float* rv2 = new float[ngd];
01737       ranlux_(rv2,ngd);
01738       double* qgd=pgd;
01739       double check[4]={0,0,0,0};
01740       for(int i=0;i<ngd;i++){
01741         double p[3];
01742         for(int j=0;j<3;j++)p[j]=*qgd++;
01743         double e = sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
01744         std::string material = "scintillator";
01745         double range = GammaComptonRange(material,e);
01746         range*=log(1/rv2[i]);
01747         for(int j=0;j<3;j++)*rgptr++ = rneu[j]+range*p[j]/e;
01748         //
01749         for(int j=0;j<3;j++)check[j]+=p[j];
01750         check[3]+=e;
01751         //
01752         *egptr++ = e;
01753         *tgptr++ = tneu+range/c*refracGd;
01754         *ogptr++ = "N";
01755         ngamma++;
01756         //cout << "  " << ngamma << " gammas" << endl;
01757       }
01758         delete[] pgd;
01759         delete[] rv2;
01760     }
01761   } // done with neutrons
01762   //cout << "  " << ngamma << " gammas, total" << endl;
01763   //
01764   // Transfer temp data to storage
01765   //
01766   delete[] Rgamma;
01767   delete[] Egamma;
01768   delete[] Tgamma;
01769   delete[] Ogamma;
01770   
01771   Ngamma = ngamma;
01772   Rgamma = new double[3*ngamma];
01773   Egamma = new double[ngamma];
01774   Tgamma = new double[ngamma];
01775   Ogamma = new std::string[ngamma];
01776 
01777   double* ptr;
01778 
01779   ptr = Rgamma+3*ngamma;
01780   while(ptr>Rgamma)*--ptr=*--rgptr;
01781 
01782   ptr = Egamma+ngamma;
01783   while(ptr>Egamma)*--ptr=*--egptr;
01784 
01785   ptr = Tgamma+ngamma;
01786   while(ptr>Tgamma)*--ptr=*--tgptr;
01787 
01788   std::string* qtr = Ogamma+ngamma;
01789   while(qtr>Ogamma)*--qtr=*--ogptr;
01790 
01791   /*
01792   cout << "  " << Ngamma << " gammas" << endl;
01793   for(int n=0;n<Ngamma;n++){
01794     for(int i=0;i<3;i++) cout << "  Rgamma[" << 3*n+i << "] " 
01795                               << Rgamma[3*n+i];
01796     cout << endl;
01797     cout << "  Egamma[" << n << "] " << Egamma[n] 
01798          << "  Tgamma[" << n << "] " << Tgamma[n]
01799          << "  Ogamma[" << n << "] " << Ogamma[n] << endl;
01800   }
01801   */
01802 
01803   // clean up
01804   delete[] rgamma;
01805   delete[] egamma;
01806   delete[] tgamma;
01807   delete[] ogamma;    
01808 }
01809 
01810 float ReactorDetector::QSpectrum(float const& Q){
01811   float m = 8.1497;
01812   float ng1 = 0.2735E-01;
01813   float ng2 = 0.7944E-03;
01814   float nexp = 0.81228E-01;
01815   float G2 = 4.5; //Relative gain for 2nd Polya for R1408 spectrum
01816   float Q2 = Q/G2;
01817   float Q0 = 0.61;
01818   float gmratio = m/exp(gamma(m));
01819 
01820 
01821   float Qpol1 = gmratio*pow((m*Q),m-1)*exp(-m*Q);
01822   float Qpol2 = gmratio*pow((m*Q2),m-1)*exp(-m*Q2);
01823   return  (ng1*Qpol1 + ng2*Qpol2 + nexp*exp(-Q/Q0));
01824 }
01825 
01826 float ReactorDetector::GenerateQ(int npe){
01827 
01828     float pedwidth = 0.1;
01829     float qtot = 0.;
01830 
01831     if (npe>0) {
01832      float* qran = new float[npe];
01833      float* qsmear= new float[npe];
01834      rnormx_(qsmear,npe,ranlux_);
01835      funlux_(Qspace,qran[0],npe);
01836       for (int ipe=0; ipe<npe; ipe++){
01837             qtot+=(qran[ipe]+pedwidth*qsmear[ipe]);
01838       }
01839         delete[] qran;
01840         delete[] qsmear;
01841     }
01842     else {
01843       float qsmear[1];
01844       int ipe = 1;
01845       rnormx_(qsmear,ipe,ranlux_);
01846       qtot+=pedwidth*qsmear[0];
01847     }
01848     return qtot;
01849  }
01850 
01851 
01852 
01853 ReactorDetector::PMT* ReactorDetector::ResizePMT(){
01854 
01855   int oldsize = PMTsize;
01856   int newsize = 2*oldsize;
01857   PMT* PMThold = new PMT[newsize];
01858   PMT* PMTold = PMTdata+oldsize;
01859   PMT* PMTnew = PMThold+oldsize;
01860   while(PMTnew>PMThold)*--PMTnew=*--PMTold;
01861   delete[] PMTdata;
01862   PMTdata = new PMT[newsize];
01863   PMTold = PMThold+newsize;
01864   PMTnew = PMTdata+newsize;
01865   while(PMTnew>PMTdata)*--PMTnew=*--PMTold;
01866   delete[] PMThold;
01867   PMTsize = newsize;
01868   return PMTdata+oldsize;
01869 }
01870 //
01871 // Fetch number of gammas from positron
01872 //
01873 int ReactorDetector::GetNgamma(std::string parent){
01874   int gamma=0;
01875   for(std::string* optr=Ogamma;optr<Ogamma+Ngamma;optr++){
01876     if (*optr==parent)gamma++;
01877   }
01878   return gamma;
01879 }
01880 //
01881 // Fetch array of positron gammas
01882 //
01883 void ReactorDetector::GetGamma(std::string parent,double* xyz){
01884 
01885   double* Rptr=Rgamma;
01886   double* rptr=xyz;
01887   //
01888   for(std::string* optr=Ogamma;optr<Ogamma+Ngamma;optr++){
01889     if (*optr==parent){
01890       for(int i=0;i<3;i++)*rptr++ = *Rptr++;
01891     }
01892   }
01893 }  
01894 void ReactorDetector::GetGamma(std::string parent,double* e,
01895                                double* t,double* r,
01896                                double* c,double* f){
01897 
01898   //cout << "GetGamma(" << parent << ")" << endl;
01899   double* Eptr=Egamma;
01900   double* Rptr=Rgamma;
01901   double* Tptr=Tgamma;
01902   //double* eptr=e;
01903   //double* rptr=r;
01904   //double* tptr=t;
01905   //double* cptr=c;
01906   //double* fptr=f;
01907   //
01908   for(std::string* optr=Ogamma;optr<Ogamma+Ngamma;optr++){
01909     if (*optr==parent){
01910       *e++=*Eptr;
01911       *t++=*Tptr;
01912       double x=*(Rptr+0);
01913       double y=*(Rptr+1);
01914       double z=*(Rptr+2);
01915       *r++=sqrt(x*x+y*y+z*z);
01916       *c++=z/sqrt(x*x+y*y+z*z);
01917       *f++=atan2(y,x);
01918       if((*f)<0)*f+=2*PI;
01919     }
01920     Eptr++;
01921     Tptr++;
01922     Rptr+=3;
01923   }
01924 }
01925 void ReactorDetector::SetNominalParameters(){
01926   //
01927   // Set nominal detector parameters from ReactorConstants
01928   //
01929   R0 = RD0;
01930   R1 = RD1;
01931   R2 = RD2;
01932   //
01933   GdConcentration=GDCONCENTRATION;
01934   //
01935   refracGd = REFRACGD;
01936   refracSc = REFRACSC;
01937   //
01938   mfpGd = MFPGD;
01939   mfpSc = MFPSC;
01940   //
01941   MeanNeuDispl = MEANNEUDISPL;
01942   //
01943   FastNeutronOption = FASTNEUTRONOPTION;
01944   Temperature = TEMPERATURE;
01945   //
01946   tGd = TGD;
01947   tSc = TSC;
01948   //
01949   GdCaptureFraction = GDCAPTUREFRACTION;
01950   //
01951   GammasPerGd = GAMMASPERGD; 
01952   //
01953   attenlGd = ATTENLGD;
01954   attenlSc = ATTENLSC;
01955   //
01956   PhotonsPerMeV = PHOTONSPERMEV;
01957   //
01958   ScintDecayTime = SCINTDECAYTIME;
01959   //
01960   PMTcoverage = PMTCOVERAGE;
01961   PMTdiameter = PMTDIAMETER;
01962   //
01963   PMTqe = PMTQE;
01964   EScale = ESCALE;
01965 }
01966 void ReactorDetector::SetPMT(){
01967   //
01968   // determine the number of PMT
01969   //
01970   float num_pmt = 16*R2*R2*PMTcoverage/PMTdiameter/PMTdiameter;
01971   Ncosbins = int(sqrt(num_pmt/PI));
01972   Nphibins = int(sqrt(num_pmt*PI));
01973   Npmt = Ncosbins*Nphibins;
01974   //cout << "  Set Npmt " << Npmt << endl;
01975 
01976   delete[] Xpmt;
01977   delete[] Ypmt;
01978   delete[] Zpmt;
01979   delete[] PMTOutputArray;
01980 
01981   Xpmt = new double[Npmt];
01982   Ypmt = new double[Npmt];
01983   Zpmt = new double[Npmt];
01984   PMTOutputArray = new PMToutput*[Npmt];
01985   for(int i=0;i<Npmt;i++){
01986     PMTOutputArray[i]=new PMToutput (i);
01987     }
01988 
01989   //
01990   // Make a guess at PMT size to avoid too much re-sizing
01991   //
01992   delete[] PMTdata;
01993   PMTlen = 0;
01994   PMTsize = 10*Npmt;
01995   PMTdata = new PMT[PMTsize];
01996   //
01997   // loop over cotheta, then phi
01998   //
01999   double* xpmt=Xpmt;
02000   double* ypmt=Ypmt;
02001   double* zpmt=Zpmt;
02002   
02003   double dphi = 2*PI/Nphibins;
02004   double dz = 2.0/Ncosbins;
02005   
02006   for(int i=0;i<Ncosbins;i++){
02007     double z = -1+i*dz+dz/2;
02008     for(int j=0;j<Nphibins;j++){
02009       double phi = j*dphi+dphi/2;
02010       *xpmt++ = R2*sqrt(1-z*z)*cos(phi);
02011       *ypmt++ = R2*sqrt(1-z*z)*sin(phi);
02012       *zpmt++ = R2*z;
02013       //      cout << "z,phi  " << z << " " << phi << endl;
02014     }
02015   }
02016 }
02017 void ReactorDetector::SetParam_GdConcentration(double x){
02018   GdConcentration =x;
02019   double y =(GdConcentration/GdConcentrationRef);
02020   double z = GdCaptureFraction/(1-GdCaptureFraction);
02021   //
02022   GdCaptureFraction = y*z/(1+y*z);
02023   //
02024   mfpGd*=(1+z)/(1+y*z);
02025   tGd*=(1+z)/(1+y*z);
02026   SetNeutronPropStuff();
02027 }
02028 void ReactorDetector::SetNeutronPropStuff(){
02029   //
02030   a[0] = A0;
02031   a[1] = A1;
02032   a[2] = A2;
02033   a[3] = A3;
02034   //
02035   z[0] = Z0;
02036   z[1] = Z1;
02037   z[2] = Z2;
02038   z[3] = Z3;
02039   //
02040   nn[0]=NAVOGADRO*f0*density/A0;
02041   nn[1]=NAVOGADRO*f1*density/A1;
02042   nn[2]=NAVOGADRO*GdConcentration/100*f2*density/A2;
02043   nn[3]=NAVOGADRO*GdConcentration/100*f3*density/A3;
02044   //
02045 }
02046 double ReactorDetector::GetHitFrac(std::string parent){
02047   bool* hit = new bool[Npmt];
02048   for(int i=0;i<Npmt;i++)hit[i]=0;
02049   for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
02050     if(pPMT->Parent==parent&&pPMT->Npe>0)hit[pPMT->npmt]=1;
02051   }
02052   int nhit=0;
02053   for(int i=0;i<Npmt;i++)nhit+=hit[i];
02054   delete[] hit;
02055   return 1.0*nhit/Npmt;
02056 }
02057 
02058 void ReactorDetector::CalculateEnergy(std::string parent){
02059 
02060   // so far this only works for positron and neutron
02061   //
02062   //double rpmt[3];
02063   double revent[3];
02064   if(parent=="P")
02065     for(int i=0;i<3;i++)revent[i]=DipolePositron[i];
02066   else if(parent=="N")
02067     for(int i=0;i<3;i++)revent[i]=DipoleNeutron[i];
02068   else
02069     cout << "ERROR: CalculateEnergy does not work for parent type: "
02070          << parent << endl;
02071 
02072   double rdotr = revent[0]*revent[0]+revent[1]*revent[1]+revent[2]*revent[2];
02073   double* xpmt=Xpmt;
02074   double* ypmt=Ypmt;
02075   double* zpmt=Zpmt;
02076   double dr[3];
02077   double tr[3];
02078   //
02079   //
02080   // Begin by calculating the probability of photon detection
02081   // integrated over the entire PMT sphere.  Event position is
02082   // reconstructed position.  Ultimately, we'd like to include the
02083   // wavelength and angular PMT response here, rather than the flat
02084   // PMTqe.  Most of this was pasted from GenerateResponse.
02085   //
02086   //
02087   double detprob = 0.;
02088   double sGd;
02089   double sSc;
02090   for(int i=0;i<Ncosbins;i++){
02091     for(int j=0;j<Nphibins;j++){
02092       double rpmt[3];
02093       rpmt[0] = *xpmt++;
02094       rpmt[1] = *ypmt++;
02095       rpmt[2] = *zpmt++;
02096       for (int l=0;l<3;l++)dr[l]=rpmt[l]-revent[l];
02097       double ddr = sqrt(dr[0]*dr[0]+dr[1]*dr[1]+dr[2]*dr[2]);
02098       for (int l=0;l<3;l++)tr[l]=dr[l]/ddr;
02099       double tdotr = revent[0]*tr[0]+revent[1]*tr[1]+revent[2]*tr[2];
02100       //
02101       // solve for intersection with inner radius
02102       //
02103       //
02104       // check for no intersections
02105       //
02106       double test = R0*R0+tdotr*tdotr-rdotr;
02107       bool inter = (test>=0);
02108       //
02109       // A straight line from the interaction vertex to the pmt
02110       // goes through the inner radius at least once.  In this case,
02111       // the light will see Gd-scint. and unloaded scint.
02112       //
02113       if(inter){
02114         double sint1 = -tdotr+sqrt(test);
02115         double sint2 = -tdotr-sqrt(test);
02116         if(sint1>0&&sint2<0){     // 1 intersection at soln 1
02117           sGd = sint1;
02118         }
02119         else if(sint2>0&&sint1<0){     // 1 intersection at soln 2
02120           sGd = sint2;
02121         }
02122         else if(sint2>0&&sint1>0){     // 2 intersections
02123           sGd = sint1-sint2;
02124           if(sGd<0)sGd*=-1;
02125         }
02126         else{                          // no intersections towards PMT
02127           sGd = 0;
02128         }
02129         sSc = ddr-sGd;
02130       }
02131       //
02132       //  otherwise the light only goes through unloaded scintillator
02133       //
02134       else{
02135         sGd = 0;
02136         sSc = ddr;
02137       }
02138       // Now we increment detection probability, including
02139       // solid angle factor (sfactor).
02140       double atten = exp(-sGd/attenlGd-sSc/attenlSc);
02141       double sfactor=0;
02142       for (int l=0;l<3;l++){sfactor+=tr[l]*rpmt[l]/R2;}
02143       sfactor *= R2*R2/ddr/ddr;
02144       detprob+=sfactor*atten*PMTqe*(PMTcoverage/Npmt);
02145     }
02146   }
02147   // Now we find out how many photons we actually detected.
02148   double q = 0.;
02149   for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
02150 //    if(pPMT->Parent==parent&&pPMT->Q>0){
02151     if(pPMT->Parent==parent){
02152       q+=pPMT->Q;
02153     }
02154   }
02155 
02156   // Now we can calculate most probable energy, including
02157   // division by energy scale.
02158   double energy = 0.;
02159   energy = q/(detprob*PhotonsPerMeV*EScale);
02160   //cout << " Escale " << EScale << endl;
02161 
02162   if(parent=="E") PHeleReco=energy;
02163   if(parent=="P") PHposReco=energy;
02164   if(parent=="N") PHneuReco=energy;
02165   if(parent=="Tl") PHthaReco=energy;
02166   if(parent=="Bi") PHbisReco=energy;
02167   if(parent=="Cf") PHcalReco=energy;
02168   if(parent=="muon") PHmuoReco=energy;
02169 
02170 }
02171 
02172 void ReactorDetector::CalculateQPosition(std::string parent){
02173 
02174   double d[3] = {0,0,0};
02175   double q=0;
02176   for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
02177     if(pPMT->Parent==parent&&pPMT->Q>0){
02178       q+=pPMT->PHsmr;
02179       for(int i=0;i<3;i++)d[i]+=(pPMT->PHsmr)*(pPMT->xyz[i]);
02180     }
02181   }
02182   if(q>0)for(int i=0;i<3;i++)d[i]/=q;
02183 
02184   if(parent=="P"){
02185     for(int i=0;i<3;i++)DipolePositron[i] = d[i];
02186   }
02187   if(parent=="N"){
02188     for(int i=0;i<3;i++)DipoleNeutron[i] = d[i];
02189   }
02190 }
02191 void ReactorDetector::ImproveQPosition(std::string parent){
02192 
02193   double d[3];
02194   double k = 0.;
02195   PMT* qPMT;
02196   if(parent=="P"){
02197     for(int i=0;i<3;i++)d[i]=DipolePositron[i];
02198     k = PHpositron;
02199     qPMT=PMTpositron;
02200   }
02201   if(parent=="N"){
02202     for(int i=0;i<3;i++)d[i]=DipoleNeutron[i];
02203     k = PHneutron;
02204     qPMT=PMTneutron;
02205   }
02206   //
02207   double C= PMTqe*PMTdiameter*PMTdiameter/16*PhotonsPerMeV;
02208   //
02209   double sumlast = 999999;
02210   double test=999999;
02211   double tol=0.01;
02212   int sign=1;
02213   bool first=1;
02214   double dev=0.1;
02215   while(test>tol){
02216   //
02217     double dlast[3];
02218     for(int i=0;i<3;i++)dlast[i]=d[i];
02219     double klast=k;
02220     double sum[3]={0,0,0};
02221     for(PMT* pPMT=qPMT;pPMT<qPMT+Npmt;pPMT++){
02222       if(pPMT->Parent!=parent)continue;
02223       double r[3];
02224       double rn[3];
02225       double t[3];
02226       for(int i=0;i<3;i++)rn[i]=pPMT->xyz[i];
02227       for(int i=0;i<3;i++)r[i]=rn[i]-d[i];
02228       double rdotr = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
02229       double rdotrn = r[0]*rn[0]+r[1]*rn[1]+r[2]*rn[2];
02230       for(int i=0;i<3;i++)t[i]=r[i]/sqrt(rdotr);
02231       //
02232       double tdotd=d[0]*t[0]+d[1]*t[1]+d[2]*t[2];
02233       double ddotd=d[0]*d[0]+d[1]*d[1]+d[2]*d[2];
02234       double test2=tdotd*tdotd-ddotd+R0*R0;
02235       double sSc;
02236       double sGd;
02237       //
02238       if(test2<0){//no intersections with inner shell
02239         sSc=sqrt(rdotr);
02240         sGd=0;
02241       }
02242       else{
02243         double s1=-tdotd+sqrt(test2);
02244         double s2=-tdotd-sqrt(test2);
02245         if(s1>0&&s2>0){//2 interesections
02246           sSc=sqrt(rdotr)-s1+s2;
02247           sGd=s1-s2;
02248         }
02249         else if(s1>0&&s2<0){//1 interestion
02250           sSc=sqrt(rdotr)-s1;
02251           sGd=s1;
02252         }
02253         else {//no interesctions
02254           sSc=sqrt(rdotr);
02255           sGd=0;
02256         }
02257       }
02258       //
02259       double mn=pPMT->Q;
02260       double w0=C*rdotrn/R2/sqrt(rdotr)/rdotr*
02261         exp(-sGd/attenlGd-sSc/attenlSc);
02262       sum[0]+=mn;
02263       sum[1]+=w0;
02264       sum[2]+=w0*klast-mn*log(w0*klast);
02265     }
02266     if(first){
02267       k = sum[0]/sum[1];
02268       for(int i=0;i<3;i++)d[i]*=(1+dev);
02269       first=0;
02270     }
02271     else{
02272       k = sum[0]/sum[1];
02273       if(sum[2]<sumlast){
02274         for(int i=0;i<3;i++)d[i]*=(1+sign*dev);
02275       }
02276       else{
02277         sign*=-1;
02278         dev/=2;
02279         for(int i=0;i<3;i++)d[i]*=(1+sign*dev);
02280       }
02281     }
02282     test=1-sum[2]/sumlast;
02283     if(test<0)test*=-1;
02284     sumlast=sum[2];
02285   }
02286   if(parent=="P")for(int i=0;i<3;i++)DipolePositron[i]=d[i];
02287   if(parent=="N")for(int i=0;i<3;i++)DipoleNeutron[i]=d[i];
02288 }  
02289 //
02290 //
02291 //
02292 void ReactorDetector::ConsolidatePMT(){
02293   //
02294   // first sum the total xpe
02295   //
02296   double *xpeele = new double[Npmt];
02297   double *Qele = new double[Npmt];
02298   double *xpepos = new double[Npmt];
02299   double *Qpos = new double[Npmt];
02300   double *xpeneu =  new double[Npmt];
02301   double *Qneu = new double[Npmt];
02302   double *xpetha = new double[Npmt];
02303   double *Qtha = new double[Npmt];
02304   double *xpebis = new double[Npmt];
02305   double *Qbis = new double[Npmt];
02306   double *xpecal = new double[Npmt];
02307   double *Qcal = new double[Npmt];
02308   double *xpemuo = new double[Npmt];
02309   double *Qmuo = new double[Npmt];
02310   int* subhitsele = new int[Npmt];
02311   int* subhitspos = new int[Npmt];
02312   int* subhitsneu = new int[Npmt];
02313   int* subhitstha = new int[Npmt];
02314   int* subhitsbis = new int[Npmt];
02315   int* subhitscal = new int[Npmt];
02316   int* subhitsmuo = new int[Npmt];
02317   //
02318   for(int i=0;i<Npmt;i++){
02319     xpeele[i]=0;
02320     xpepos[i]=0;
02321     xpeneu[i]=0;
02322     xpetha[i]=0;
02323     xpebis[i]=0;
02324     xpecal[i]=0;
02325     xpemuo[i]=0;
02326     Qele[i]=0;
02327     Qpos[i]=0;
02328     Qneu[i]=0;
02329     Qtha[i]=0;
02330     Qbis[i]=0;
02331     Qcal[i]=0;
02332     Qmuo[i]=0;
02333     subhitsele[i]=0;
02334     subhitspos[i]=0;
02335     subhitsneu[i]=0;
02336     subhitstha[i]=0;
02337     subhitsbis[i]=0;
02338     subhitscal[i]=0;
02339     subhitsmuo[i]=0;
02340   }
02341   //
02342   for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
02343     int n=pPMT->npmt;
02344     PMTOutputArray[n]->incout(pPMT, MergeTime);
02345     if(pPMT->Parent=="E"){
02346       xpeele[n]+=pPMT->Npe;
02347       Qele[n]+=pPMT->Q;
02348       subhitsele[n]++;
02349     }
02350     if(pPMT->Parent=="P"){
02351       xpepos[n]+=pPMT->Npe;
02352       Qpos[n]+=pPMT->Q;
02353       subhitspos[n]++;
02354     }
02355     if(pPMT->Parent=="N"){
02356       xpeneu[n]+=pPMT->Npe;
02357       Qneu[n]+=pPMT->Q;
02358       subhitsneu[n]++;
02359     }
02360     if(pPMT->Parent=="Tl"){
02361       xpetha[n]+=pPMT->Npe;
02362       Qtha[n]+=pPMT->Q;
02363       subhitstha[n]++;
02364     }
02365     if(pPMT->Parent=="Bi"){
02366       xpebis[n]+=pPMT->Npe;
02367       Qbis[n]+=pPMT->Q;
02368       subhitsbis[n]++;
02369     }
02370     if(pPMT->Parent=="Cf"){
02371       xpecal[n]+=pPMT->Npe;
02372       Qcal[n]+=pPMT->Q;
02373       subhitscal[n]++;
02374     }
02375     if(pPMT->Parent=="muon"){
02376       xpemuo[n]+=pPMT->Npe;
02377       Qmuo[n]+=pPMT->Q;
02378       subhitsmuo[n]++;
02379     }
02380   }
02381   //   
02382   delete[] PMTelectron;
02383   PMTelectron = new PMT[Npmt];
02384   delete[] PMTpositron;
02385   PMTpositron = new PMT[Npmt];
02386   delete[] PMTneutron;
02387   PMTneutron = new PMT[Npmt];
02388   delete[] PMTthallium;
02389   PMTthallium = new PMT[Npmt];
02390   delete[] PMTbismuth;
02391   PMTbismuth = new PMT[Npmt];
02392   delete[] PMTcalifornium;
02393   PMTcalifornium = new PMT[Npmt];
02394   delete[] PMTmuon;
02395   PMTmuon = new PMT[Npmt];
02396   //
02397   for(int i=0;i<Npmt;i++){
02398     (*(PMTelectron+i)).npmt=i;
02399     (*(PMTelectron+i)).Parent="E";
02400     (PMTelectron+i)->SetSubHitsPtr(subhitsele[i]);
02401     //
02402     (*(PMTpositron+i)).npmt=i;
02403     (*(PMTpositron+i)).Parent="P";
02404     (PMTpositron+i)->SetSubHitsPtr(subhitspos[i]);
02405     //
02406     (*(PMTneutron+i)).npmt=i;
02407     (*(PMTneutron+i)).Parent="N";
02408     (PMTneutron+i)->SetSubHitsPtr(subhitsneu[i]);
02409     //
02410     (*(PMTthallium+i)).npmt=i;
02411     (*(PMTthallium+i)).Parent="Tl";
02412     (PMTthallium+i)->SetSubHitsPtr(subhitstha[i]);
02413     //
02414     (*(PMTbismuth+i)).npmt=i;
02415     (*(PMTbismuth+i)).Parent="Bi";
02416     (PMTbismuth+i)->SetSubHitsPtr(subhitsbis[i]);
02417     //
02418     (*(PMTcalifornium+i)).npmt=i;
02419     (*(PMTcalifornium+i)).Parent="Cf";
02420     (PMTcalifornium+i)->SetSubHitsPtr(subhitscal[i]);
02421     //
02422     (*(PMTmuon+i)).npmt=i;
02423     (*(PMTmuon+i)).Parent="muon";
02424     (PMTmuon+i)->SetSubHitsPtr(subhitsmuo[i]);
02425   }
02426   //
02427   for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
02428     int n=pPMT->npmt;
02429     PMT* qPMT;
02430     if(pPMT->Parent=="E"){
02431       qPMT = PMTelectron+n;
02432       *qPMT+=*pPMT;
02433     }
02434     if(pPMT->Parent=="P"){
02435       qPMT = PMTpositron+n;
02436       *qPMT+=*pPMT;
02437     }
02438     if(pPMT->Parent=="N"){
02439       qPMT = PMTneutron+n;
02440       *qPMT+=*pPMT;
02441     }
02442     if(pPMT->Parent=="Tl"){
02443       qPMT = PMTthallium+n;
02444       *qPMT+=*pPMT;
02445     }
02446     if(pPMT->Parent=="Bi"){
02447       qPMT = PMTbismuth+n;
02448       *qPMT+=*pPMT;
02449     }
02450     if(pPMT->Parent=="Cf"){
02451       qPMT = PMTcalifornium+n;
02452       *qPMT+=*pPMT;
02453     }
02454     if(pPMT->Parent=="muon"){
02455       qPMT = PMTmuon+n;
02456       *qPMT+=*pPMT;
02457     }
02458   }
02459   delete[] xpepos;
02460   delete[] Qpos;
02461   delete[] xpeneu;
02462   delete[] Qneu;
02463   delete[] xpeele;
02464   delete[] Qele;
02465   delete[] xpetha;
02466   delete[] Qtha;
02467   delete[] xpebis;
02468   delete[] Qbis;
02469   delete[] Qcal;
02470   delete[] xpecal;
02471   delete[] Qmuo;
02472   delete[] xpemuo;
02473 }    
02474 
02475 double* ReactorDetector::TranslateRadialDistance(double x,double y,
02476                                                  double z,double d){
02477 
02478   // translate a point along the radial axis
02479   // x,y,z cartesian coordinates of point to translate
02480   // d distance to translate along the radial axis (negative -> move inward)
02481 
02482   double p[3];
02483   double* ptr = p;
02484   double r, costheta, phi;
02485 
02486   // convert to spherical coordinates
02487   //
02488   r = sqrt(x*x + y*y + z*z);
02489   if(r != 0.){
02490     costheta = z/r;
02491   }
02492   else{
02493     costheta = 0.;
02494     cout << " ERROR: TranslateRadialDistance initial r = 0!" << endl;
02495   }
02496   phi = atan2(y,x);
02497   if(phi < 0) phi += 2*PI;
02498 
02499   //cout << " " << r << "," << costheta << "," << phi << endl;
02500 
02501   // move a distance d along the radial axis
02502   //
02503   r+=d;
02504 
02505   // return to cartesian
02506   //
02507   p[0] = r*sqrt(1-costheta*costheta)*cos(phi);
02508   p[1] = r*sqrt(1-costheta*costheta)*sin(phi);
02509   p[2] = r*costheta;
02510 
02511   //for(int i=0;i<3;i++) cout << " p["<<i<<"] " << p[i] << " ";
02512   //cout << endl;
02513 
02514   // for some reason this empty cout statement is required to
02515   // make this funtion work and I don't like it AT ALL
02516   cout << "";
02517 
02518   return(ptr);
02519 }
02520 
02521 void ReactorDetector::CleanPMTOutputArray(){
02522   for(int j=0; j<Npmt; j++)
02523     PMTOutputArray[j]->cleanup();
02524   }
02525 
02526 //PMToutput constructor, initializes empty vector and array
02527 PMToutput::PMToutput(int ID){
02528   PMTid=ID;
02529   eventsin=0;
02530   eventsdetected=0;
02531 }
02532 
02533 //PMToutput destructor
02534 PMToutput::~PMToutput(){
02535 }
02536 //PMToutput copy constructor
02537 PMToutput::PMToutput(const PMToutput& PMTout){
02538   PMTid=PMTout.PMTid;
02539   eventsin=PMTout.eventsin;
02540   eventsdetected=PMTout.eventsdetected;
02541   temp=PMTout.temp;
02542   copy(PMTout.events.begin(),PMTout.events.end(),events.begin());
02543   copy(PMTout.output.begin(),PMTout.output.end(),output.begin());
02544   }
02545 
02546 //PMToutput overloaded = operator
02547   PMToutput& PMToutput::operator=(const PMToutput& rhs){
02548   PMTid=rhs.PMTid;
02549   eventsin=rhs.eventsin;
02550   eventsdetected=rhs.eventsdetected;
02551   int i;
02552   for (i=0;i<eventsdetected;i++) output.push_back(rhs.output[i]);
02553   for (i=0;i<eventsin;i++) events.push_back(rhs.events[i]);
02554   return *this;
02555   }
02556 
02557 //function incout - adds a PMT event into the string
02558 void PMToutput::incout(ReactorDetector::PMT* input, double MergeTime){
02559   temp.energy=input->PHsmr;
02560   if(temp.energy==0)return;
02561   temp.time = input->TimeSmr;
02562   events.push_back(input);
02563   eventsin++;
02564   if(eventsdetected==0) {
02565     output.push_back(temp);
02566     eventsdetected++;
02567   }
02568   else{
02569     bool done=false;
02570     int x,y,z;
02571     x=0;
02572     z=output.size()-1;
02573     if(temp.time < output[0].time + MergeTime){
02574       if(abs(temp.time - output[0].time) < MergeTime) output[0]+=temp;
02575       else output.insert(output.begin(),temp);
02576       done=true;
02577       }
02578     else if(temp.time > output[z].time - MergeTime){
02579       if(abs(temp.time - output[z].time) < MergeTime) output[z]+=temp;
02580       else output.push_back(temp);
02581       done=true;
02582       }
02583     while(!done && z-x>1){
02584       y=(x+z)/2;
02585       if(temp.time > output[y].time) x=y;
02586       else z=y;
02587       }
02588     if(!done){
02589       if(abs(temp.time - output[x].time) < MergeTime){
02590         output[x]+=temp;
02591         done=true;
02592         }
02593       else if(abs(temp.time-output[z].time) < MergeTime){
02594         output[z]+=temp;
02595         done=true;
02596         }
02597       }
02598     if(!done){
02599       std::vector<PMTinfo>::iterator it1;
02600       it1=output.begin();
02601       it1+=z;
02602       output.insert(it1,temp);
02603       eventsdetected++;
02604       }
02605 /*  <Algorithm> implementation
02606     std::vector<PMTinfo>::iterator it1;
02607     it2= lower_bound(output.begin(), output.end(), temp, SortByTime);
02608     if ( it1 != output.begin() )
02609       it1--;
02610     if ( temp.time - it1->time < MergeTime ) {
02611       (*it1)+=temp;
02612     }
02613     else if ( it2 != output.end() && it2->time - temp.time < MergeTime ) {
02614       (*it2)+=temp;
02615     }
02616     else {
02617       output.insert(it2, temp);
02618       eventsdetected++;
02619     }
02620 */}
02621 }
02622 
02623 void PMToutput::testdump(int event){
02624   ofstream file ( "testdump.txt", ios::app );
02625   if(!file.is_open()){
02626     cout << "File failure" << endl;
02627     return;
02628     }
02629   else{
02630     std::vector<PMTinfo>::iterator it1;
02631     for (it1 = output.begin(); it1 != output.end(); it1++)
02632       file << event << "\t" << PMTid << "\t" << it1->time << "\t" << it1->energy << endl;
02633     }
02634 }
02635 
02636 
02637 void PMToutput::cleanup(){
02638   output.clear();
02639   events.clear();
02640   eventsin=0;
02641   eventsdetected=0;
02642   }
02643 
02644 PMTinfo&
02645 PMTinfo::operator+=(const PMTinfo& rhs){
02646   energy+=rhs.energy;
02647   if(time>rhs.time) time=rhs.time;
02648   return *this;
02649 }
02650 

Generated on Mon Feb 21 16:11:19 2005 for ReactorFsim by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002