Main Page   Compound List   File List   Compound Members   File Members  

ReactorNtuple.cpp

Go to the documentation of this file.
00001 
00008 #include <iostream.h>
00009 #include <math.h>
00010 #include <string.h>
00011 #include <cstdio>
00012 #include "ReactorEvent.hh"
00013 #include "ReactorDetector.hh"
00014 #include "ReactorFortran.hh"
00015 #include "ReactorNtuple.hh"
00016 
00017 ReactorNtuple::ReactorNtuple(const ReactorNtuple& Ntuple){
00018 
00019   Id = Ntuple.Id;
00020   Nvar = Ntuple.Nvar;
00021   Offset = Ntuple.Offset;
00022 }
00023 ReactorNtuple& ReactorNtuple::operator=(const ReactorNtuple& rhs){
00024 
00025   if(this == &rhs)return *this;
00026   Id = rhs.Id;
00027   Nvar = rhs.Nvar;
00028   Offset = rhs.Offset;
00029   return *this;
00030 }
00031 
00032 ReactorNtuple::ReactorNtuple(int ntId,
00033                              const char* ntName,
00034                              char* dirName){
00035 
00036   // define ntuple tags here
00037   // to add vars, increment ntVars, and add tag to end of list
00038   // don't change wordsize
00039   //
00040   int wordsize = 10;
00041   int ntVarsMax = 150;
00042   int ntVars = 0;
00043   char *ntTags = new char[ntVarsMax*wordsize];
00044   //
00045   // book the general event veriables
00046   //
00047   strcpy(ntTags+ntVars++*wordsize,"evtnum   ");
00048   strcpy(ntTags+ntVars++*wordsize,"wtc      ");
00049   strcpy(ntTags+ntVars++*wordsize,"wtf      ");
00050   //
00051   // book the nubar event variables
00052   //
00053   strcpy(ntTags+ntVars++*wordsize,"enu      ");
00054   strcpy(ntTags+ntVars++*wordsize,"vnu1r    ");
00055   strcpy(ntTags+ntVars++*wordsize,"vnu1cos  ");
00056   strcpy(ntTags+ntVars++*wordsize,"vnu1phi  ");
00057   strcpy(ntTags+ntVars++*wordsize,"vnu1x    ");
00058   strcpy(ntTags+ntVars++*wordsize,"vnu1y    ");
00059   strcpy(ntTags+ntVars++*wordsize,"vnu1z    ");
00060   //
00061   // book the electron event variables
00062   //
00063   strcpy(ntTags+ntVars++*wordsize,"nele     ");
00064   strcpy(ntTags+ntVars++*wordsize,"KEele1   ");
00065   strcpy(ntTags+ntVars++*wordsize,"ele1cos  ");
00066   strcpy(ntTags+ntVars++*wordsize,"vele1r   ");
00067   strcpy(ntTags+ntVars++*wordsize,"vele1cos ");
00068   strcpy(ntTags+ntVars++*wordsize,"vele1phi ");
00069   strcpy(ntTags+ntVars++*wordsize,"vele1x   ");
00070   strcpy(ntTags+ntVars++*wordsize,"vele1y   ");
00071   strcpy(ntTags+ntVars++*wordsize,"vele1z   ");
00072   strcpy(ntTags+ntVars++*wordsize,"tele1    ");
00073   //
00074   // book the positron event variables
00075   //
00076   strcpy(ntTags+ntVars++*wordsize,"npos     ");
00077   strcpy(ntTags+ntVars++*wordsize,"KEpos1   ");
00078   strcpy(ntTags+ntVars++*wordsize,"pos1cos  ");
00079   strcpy(ntTags+ntVars++*wordsize,"vpos1r   ");
00080   strcpy(ntTags+ntVars++*wordsize,"vpos1cos ");
00081   strcpy(ntTags+ntVars++*wordsize,"vpos1phi ");
00082   strcpy(ntTags+ntVars++*wordsize,"vpos1x   ");
00083   strcpy(ntTags+ntVars++*wordsize,"vpos1y   ");
00084   strcpy(ntTags+ntVars++*wordsize,"vpos1z   ");
00085   strcpy(ntTags+ntVars++*wordsize,"tpos1    ");
00086   //
00087   // book the neutron event variables
00088   //
00089   strcpy(ntTags+ntVars++*wordsize,"nneu     ");
00090   strcpy(ntTags+ntVars++*wordsize,"KEneu1   ");
00091   strcpy(ntTags+ntVars++*wordsize,"neu1cos  ");
00092   strcpy(ntTags+ntVars++*wordsize,"neuproc1 ");
00093   strcpy(ntTags+ntVars++*wordsize,"vneu1r   ");
00094   strcpy(ntTags+ntVars++*wordsize,"vneu1cos ");
00095   strcpy(ntTags+ntVars++*wordsize,"vneu1phi ");
00096   strcpy(ntTags+ntVars++*wordsize,"vneu1x   ");
00097   strcpy(ntTags+ntVars++*wordsize,"vneu1y   ");
00098   strcpy(ntTags+ntVars++*wordsize,"vneu1z   ");
00099   strcpy(ntTags+ntVars++*wordsize,"tneu1    ");
00100   strcpy(ntTags+ntVars++*wordsize,"nsteps1  ");
00101   strcpy(ntTags+ntVars++*wordsize,"absorb1  ");
00102 
00103   strcpy(ntTags+ntVars++*wordsize,"KEneu2   ");
00104   strcpy(ntTags+ntVars++*wordsize,"neu2cos  ");
00105   strcpy(ntTags+ntVars++*wordsize,"neuproc2 ");
00106   strcpy(ntTags+ntVars++*wordsize,"vneu2r   ");
00107   strcpy(ntTags+ntVars++*wordsize,"vneu2cos ");
00108   strcpy(ntTags+ntVars++*wordsize,"vneu2phi ");
00109   strcpy(ntTags+ntVars++*wordsize,"vneu2x   ");
00110   strcpy(ntTags+ntVars++*wordsize,"vneu2y   ");
00111   strcpy(ntTags+ntVars++*wordsize,"vneu2z   ");
00112   strcpy(ntTags+ntVars++*wordsize,"tneu2    ");
00113   strcpy(ntTags+ntVars++*wordsize,"nsteps2  ");
00114   strcpy(ntTags+ntVars++*wordsize,"absorb2  ");
00115 
00116   strcpy(ntTags+ntVars++*wordsize,"KEneu3   ");
00117   strcpy(ntTags+ntVars++*wordsize,"neu3cos  ");
00118   strcpy(ntTags+ntVars++*wordsize,"neuproc3 ");
00119   strcpy(ntTags+ntVars++*wordsize,"vneu3r   ");
00120   strcpy(ntTags+ntVars++*wordsize,"vneu3cos ");
00121   strcpy(ntTags+ntVars++*wordsize,"vneu3phi ");
00122   strcpy(ntTags+ntVars++*wordsize,"vneu3x   ");
00123   strcpy(ntTags+ntVars++*wordsize,"vneu3y   ");
00124   strcpy(ntTags+ntVars++*wordsize,"vneu3z   ");
00125   strcpy(ntTags+ntVars++*wordsize,"tneu3    ");
00126   strcpy(ntTags+ntVars++*wordsize,"nsteps3  ");
00127   strcpy(ntTags+ntVars++*wordsize,"absorb3  ");
00128   //
00129   // book the pmt detector variables
00130   //
00131   strcpy(ntTags+ntVars++*wordsize,"ncosbins ");
00132   strcpy(ntTags+ntVars++*wordsize,"nphibins ");
00133   //
00134   strcpy(ntTags+ntVars++*wordsize,"phposvis ");
00135   strcpy(ntTags+ntVars++*wordsize,"phposrec ");
00136   strcpy(ntTags+ntVars++*wordsize,"phneuvis ");
00137   strcpy(ntTags+ntVars++*wordsize,"phneurec ");
00138   strcpy(ntTags+ntVars++*wordsize,"phtlvis  ");
00139   strcpy(ntTags+ntVars++*wordsize,"phtlrec  ");
00140   strcpy(ntTags+ntVars++*wordsize,"phbivis  ");
00141   strcpy(ntTags+ntVars++*wordsize,"phbirec  ");
00142   strcpy(ntTags+ntVars++*wordsize,"phcfvis  ");
00143   strcpy(ntTags+ntVars++*wordsize,"phcfrec  ");
00144   strcpy(ntTags+ntVars++*wordsize,"phmuovis ");
00145   strcpy(ntTags+ntVars++*wordsize,"phmuorec ");
00146   strcpy(ntTags+ntVars++*wordsize,"phelevis ");
00147   strcpy(ntTags+ntVars++*wordsize,"phelerec ");
00148   //
00149   strcpy(ntTags+ntVars++*wordsize,"npepos   ");
00150   strcpy(ntTags+ntVars++*wordsize,"Qpos     ");
00151   strcpy(ntTags+ntVars++*wordsize,"attenpos ");
00152   strcpy(ntTags+ntVars++*wordsize,"sfactpos ");
00153   strcpy(ntTags+ntVars++*wordsize,"npeneu   ");
00154   strcpy(ntTags+ntVars++*wordsize,"Qneu     ");
00155   strcpy(ntTags+ntVars++*wordsize,"attenneu ");
00156   strcpy(ntTags+ntVars++*wordsize,"sfactneu ");
00157   strcpy(ntTags+ntVars++*wordsize,"npetl    ");
00158   strcpy(ntTags+ntVars++*wordsize,"Qtl      ");
00159   strcpy(ntTags+ntVars++*wordsize,"attentl  ");
00160   strcpy(ntTags+ntVars++*wordsize,"sfacttl  ");
00161   strcpy(ntTags+ntVars++*wordsize,"npebi    ");
00162   strcpy(ntTags+ntVars++*wordsize,"Qbi      ");
00163   strcpy(ntTags+ntVars++*wordsize,"attenbi  ");
00164   strcpy(ntTags+ntVars++*wordsize,"sfactbi  ");
00165   strcpy(ntTags+ntVars++*wordsize,"npecf    ");
00166   strcpy(ntTags+ntVars++*wordsize,"Qcf      ");
00167   strcpy(ntTags+ntVars++*wordsize,"attencf  ");
00168   strcpy(ntTags+ntVars++*wordsize,"sfactcf  ");
00169   strcpy(ntTags+ntVars++*wordsize,"npemuo   ");
00170   strcpy(ntTags+ntVars++*wordsize,"Qmuo     ");
00171   strcpy(ntTags+ntVars++*wordsize,"attenmuo ");
00172   strcpy(ntTags+ntVars++*wordsize,"sfactmuo ");
00173   strcpy(ntTags+ntVars++*wordsize,"npeele   ");
00174   strcpy(ntTags+ntVars++*wordsize,"Qele     ");
00175   strcpy(ntTags+ntVars++*wordsize,"attenele ");
00176   strcpy(ntTags+ntVars++*wordsize,"sfactele ");
00177   //
00178   strcpy(ntTags+ntVars++*wordsize,"ngampos  ");
00179   strcpy(ntTags+ntVars++*wordsize,"ngamneu  ");
00180   strcpy(ntTags+ntVars++*wordsize,"ngamtl   ");
00181   strcpy(ntTags+ntVars++*wordsize,"ngambi   ");
00182   strcpy(ntTags+ntVars++*wordsize,"ngamcf   ");
00183   strcpy(ntTags+ntVars++*wordsize,"ngammuo  ");
00184   strcpy(ntTags+ntVars++*wordsize,"ngamele  ");
00185   //
00186   strcpy(ntTags+ntVars++*wordsize,"fpmtpos  ");
00187   strcpy(ntTags+ntVars++*wordsize,"fpmtneu  ");
00188   strcpy(ntTags+ntVars++*wordsize,"fpmttl   ");
00189   strcpy(ntTags+ntVars++*wordsize,"fpmtbi   ");
00190   strcpy(ntTags+ntVars++*wordsize,"fpmtcf   ");
00191   strcpy(ntTags+ntVars++*wordsize,"fpmtmuo  ");
00192   strcpy(ntTags+ntVars++*wordsize,"fpmtele  ");
00193   //
00194   strcpy(ntTags+ntVars++*wordsize,"hitpmtpo ");
00195   strcpy(ntTags+ntVars++*wordsize,"hitpmtne ");
00196   strcpy(ntTags+ntVars++*wordsize,"hitpmtth ");
00197   strcpy(ntTags+ntVars++*wordsize,"hitpmtbi ");
00198   strcpy(ntTags+ntVars++*wordsize,"hitpmtca ");
00199   strcpy(ntTags+ntVars++*wordsize,"hitpmtmu ");
00200   strcpy(ntTags+ntVars++*wordsize,"hitpmtel ");
00201   //
00202   hbookn_(ntId,ntName,ntVars,dirName,1000,ntTags,
00203           strlen(ntName),strlen(dirName),wordsize,ntVars);
00204   //
00205   // define some histograms
00206   // Last argument is need for c++ ==> f strings
00207   //
00208   Nvar = ntVars;
00209   Id = ntId;
00210   Offset = 10000*(ntId-1);
00211 
00212   hbook1_(Offset+101,"electron PH /PMT/gamma/electron",200,0,0.1,0,
00213           strlen("electron PH /PMT/gamma/electron"));
00214   hbook1_(Offset+102,"positron PH /PMT/gamma/positron",200,0,0.1,0,
00215           strlen("positron PH /PMT/gamma/positron"));
00216   hbook1_(Offset+103,"neutron PH /PMT/gamma/neutron",200,0,0.1,0,
00217           strlen("neutron PH /PMT/gamma/neutron"));
00218   hbook1_(Offset+104,"thallium PH /PMT/gamma/decay",200,0,0.1,0,
00219           strlen("thallium PH /PMT/gamma/decay"));
00220   hbook1_(Offset+105,"bismuth PH /PMT/gamma/decay",200,0,0.1,0,
00221           strlen("bismuth PH /PMT/gamma/decay"));
00222   hbook1_(Offset+106,"cf prompt gamma PH /PMT/gamma/decay",200,0,0.1,0,
00223           strlen("cf prompt gamma PH /PMT/gamma/decay"));
00224   hbook1_(Offset+107,"muon PH /PMT/gamma/muon",200,0,0.1,0,
00225           strlen("muon PH /PMT/gamma/muon"));
00226   //
00227   hbook1_(Offset+201,"electron time /PMT/gamma/electron",250,0,5e8,0,
00228           strlen("electron time /PMT/gamma/electron"));
00229   hbook1_(Offset+202,"positron time /PMT/gamma/positron",250,0,50,0,
00230           strlen("positron time /PMT/gamma/positron"));
00231   hbook1_(Offset+203,"neutron time /PMT/gamma/neutron",250,0,5e5,0,
00232           strlen("neutron time /PMT/gamma/neutron"));
00233   hbook1_(Offset+204,"thallium time /PMT/gamma/decay",250,0,1e6,0,
00234           strlen("thallium time /PMT/gamma/decay"));
00235   hbook1_(Offset+205,"bismuth time /PMT/gamma/decay",250,0,1e6,0,
00236           strlen("bismuth time /PMT/gamma/decay"));
00237   hbook1_(Offset+206,"cf prompt gamma time /PMT/gamma/decay",250,0,50,0,
00238           strlen("cf prompt gamma time /PMT/gamma/decay"));
00239   hbook1_(Offset+207,"muon time /PMT/gamma/muon",300,0,60,0,
00240           strlen("muon time /PMT/gamma/muon"));
00241   //
00242   hbook1_(Offset+301,"number of gamma/electron",10,0,10,0,
00243           strlen("number of gamma/electron"));
00244   hbook1_(Offset+302,"number of gamma/positron",30,0,30,0,
00245           strlen("number of gamma/positron"));
00246   hbook1_(Offset+303,"number of gamma/neutron",20,0,20,0,
00247           strlen("number of gamma/neutron"));
00248   hbook1_(Offset+306,"number of prompt Cf gamma/decay",20,0,20,0,
00249           strlen("number of prompt CF gamma/decay"));
00250   //
00251   hbook1_(Offset+401,"electron gamma energy",100,0,10,0,
00252           strlen("electron gamma energy"));
00253   hbook1_(Offset+402,"positron gamma energy",100,0,10,0,
00254           strlen("positron gamma energy"));
00255   hbook1_(Offset+403,"neutron gamma energy",100,0,10,0,
00256           strlen("neutron gamma energy"));
00257   hbook1_(Offset+406,"prompt Cf gamma energy",100,0,10,0,
00258           strlen("prompt Cf gamma energy"));
00259   //
00260   hbook1_(Offset+501,"electron gamma time",250,0,5e8,0,
00261           strlen("electron gamma time"));
00262   hbook1_(Offset+502,"positron gamma time",250,0,5,0,
00263           strlen("positron gamma time"));
00264   hbook1_(Offset+503,"neutron gamma time",250,0,5e5,0,
00265           strlen("neutron gamma time"));
00266   hbook1_(Offset+506,"prompt Cf gamma time",250,0,5,0,
00267           strlen("prompt Cf gamma time"));
00268   //
00269   hbook1_(Offset+601,"electron gamma radius",400,0,400,0,
00270           strlen("electron gamma radius"));
00271   hbook1_(Offset+602,"positron gamma radius",400,0,400,0,
00272           strlen("positron gamma radius"));
00273   hbook1_(Offset+603,"neutron gamma radius",400,0,400,0,
00274           strlen("neutron gamma radius"));
00275   hbook1_(Offset+606,"prompt Cf gamma radius",400,0,400,0,
00276           strlen("prompt Cf gamma radius"));
00277   //
00278   hbook1_(Offset+701,"electron gamma cos(theta)",20,-1.1,1.1,0,
00279           strlen("electron gamma cos(theta)"));
00280   hbook1_(Offset+702,"positron gamma cos(theta)",20,-1.1,1.1,0,
00281           strlen("positron gamma cos(theta)"));
00282   hbook1_(Offset+703,"neutron gamma cos(theta)",20,-1.1,1.1,0,
00283           strlen("neutron gamma cos(theta)"));
00284   hbook1_(Offset+706,"prompt Cf gamma cos(theta)",20,-1.1,1.1,0,
00285           strlen("prompt Cf gamma cos(theta)"));
00286   //
00287   hbook1_(Offset+1001,"npe in each PMT from each electron",200,-0.5,199.5,0,
00288           strlen("npe in each PMT from each electron"));
00289   hbook1_(Offset+1002,"npe in each PMT from each positron",200,-0.5,199.5,0,
00290           strlen("npe in each PMT from each positron"));
00291   hbook1_(Offset+1003,"npe in each PMT from each neutron",200,-0.5,199.5,0,
00292           strlen("npe in each PMT from each neutron"));
00293   hbook1_(Offset+1004,"npe in each PMT from each thallium",200,-0.5,199.5,0,
00294           strlen("npe in each PMT from each thallium"));
00295   hbook1_(Offset+1005,"npe in each PMT from each bismuth",200,-0.5,199.5,0,
00296           strlen("npe in each PMT from each bismuth"));
00297   hbook1_(Offset+1006,"npe in each PMT from each calif",200,-0.5,199.5,0,
00298           strlen("npe in each PMT from each calif"));
00299   hbook1_(Offset+1007,"npe in each PMT from each muon",200,-0.5,199.5,0,
00300           strlen("npe in each PMT from each muon"));
00301   //
00302   hbook1_(Offset+2001,"total npe in each PMT from electrons",
00303           200,-0.5,199.5,0,
00304           strlen("total npe in each PMT from electrons"));
00305   hbook1_(Offset+2002,"total npe in each PMT from positrons",
00306           200,-0.5,199.5,0,
00307           strlen("total npe in each PMT from positrons"));
00308   hbook1_(Offset+2003,"total npe in each PMT from neutrons",
00309           200,-0.5,199.5,0,
00310           strlen("total npe in each PMT from neutrons"));
00311   hbook1_(Offset+2004,"total npe in each PMT from thallium",
00312           200,-0.5,199.5,0,
00313           strlen("total npe in each PMT from thallium"));
00314   hbook1_(Offset+2005,"total npe in each PMT from bismuth",
00315           200,-0.5,199.5,0,
00316           strlen("total npe in each PMT from bismuth"));
00317   hbook1_(Offset+2006,"total npe in each PMT from calif",
00318           200,-0.5,199.5,0,
00319           strlen("total npe in each PMT from calif"));
00320   hbook1_(Offset+2007,"total npe in each PMT from muons",
00321           800,-0.5,799.5,0,
00322           strlen("total npe in each PMT from muons"));
00323 }
00324 
00325 ReactorNtuple::~ReactorNtuple(){
00326   ;
00327 }
00328 
00329 void ReactorNtuple::Fill(ReactorEvent& Event,ReactorDetector& Detector){
00330 
00331   //cout << "Nvar " << Nvar << endl;
00332   float* Tuple = new float[Nvar+2];
00333   float* tuple = Tuple;
00334   int nmax = 0;
00335   int nloop = 0;
00336   //
00337   // fill the general event variables
00338   //
00339   *tuple++ = Event.GetNevent();
00340   *tuple++ = Event.GetXCWeight();
00341   *tuple++ = Event.GetFluxWeight();
00342   //
00343   // fill the nubar event variables
00344   //
00345   // get number of nubar
00346   int nnu  = Event.GetNumNubar();
00347   //cout << nnu << " nubar(s)" << endl;
00348   //
00349   // deal with multiple nubar
00350   double* rnu = new double[nnu];
00351   double* cnu = new double[nnu];
00352   double* fnu = new double[nnu];
00353   double* xnu = new double[nnu];
00354   double* ynu = new double [nnu];
00355   double* znu = new double[nnu];
00356   //
00357   // up to 1 nubar per event, keep the first one
00358   //
00359   nmax = 1;
00360   if(nnu < nmax)
00361     nloop = nnu;
00362   else
00363     nloop = nmax;
00364   //
00365   for(int n=0;n<nloop;n++){
00366     *tuple++ = Event.GetEnu(n);
00367     //
00368     rnu[n]   = Event.GetNubarVertexR(n);
00369     cnu[n]   = Event.GetNubarVertexCosTheta(n);
00370     fnu[n]   = Event.GetNubarVertexPhi(n);
00371     xnu[n]   = rnu[n]*cos(fnu[n])*sqrt(1-cnu[n]*cnu[n]);
00372     ynu[n]   = rnu[n]*sin(fnu[n])*sqrt(1-cnu[n]*cnu[n]);
00373     znu[n]   = rnu[n]*cnu[n];
00374     //
00375     *tuple++ = rnu[n];
00376     *tuple++ = cnu[n];
00377     *tuple++ = fnu[n];
00378     *tuple++ = xnu[n];
00379     *tuple++ = ynu[n];
00380     *tuple++ = znu[n];
00381   }
00382   delete[] rnu; 
00383   delete[] cnu; 
00384   delete[] fnu; 
00385   delete[] xnu; 
00386   delete[] ynu; 
00387   delete[] znu;
00388   // zero out remaining tuples
00389   for(int n=nloop;n<nmax;n++){
00390     for(int i=0;i<(7/nmax);i++) *tuple++ = 0;
00391   }
00392   //
00393   // fill the electron event variables
00394   //
00395   // get number of electron
00396   int nele = Event.GetNumElectron();
00397   //cout << nele << " electron(s)" << endl;
00398   *tuple++ = nele;
00399   //
00400   // deal with multiple electrons
00401   double* rele = new double[nele];
00402   double* cele = new double[nele];
00403   double* fele = new double[nele];
00404   double* xele = new double[nele];
00405   double* yele = new double[nele];
00406   double* zele = new double[nele];
00407   //
00408   // up to 1 electron per event, keep the first one
00409   //
00410   nmax = 1;
00411   if(nele < nmax)
00412     nloop = nele;
00413   else
00414     nloop = nmax;
00415   //
00416   for(int n=0;n<nloop;n++){
00417     *tuple++ = Event.GetElectronKE(n);
00418     *tuple++ = Event.GetElectronCosTheta(n);
00419     //
00420     rele[n]  = Event.GetElectronVertexR(n);
00421     cele[n]  = Event.GetElectronVertexCosTheta(n);
00422     fele[n]  = Event.GetElectronVertexPhi(n);
00423     xele[n]  = rele[n]*cos(fele[n])*sqrt(1-cele[n]*cele[n]);
00424     yele[n]  = rele[n]*sin(fele[n])*sqrt(1-cele[n]*cele[n]);
00425     zele[n]  = rele[n]*cele[n];
00426     //
00427     *tuple++ = rele[n];
00428     *tuple++ = cele[n];
00429     *tuple++ = fele[n];
00430     *tuple++ = xele[n];
00431     *tuple++ = yele[n];
00432     *tuple++ = zele[n];
00433     //
00434     *tuple++ = Event.GetElectronTime(n);
00435   }
00436   delete[] rele; 
00437   delete[] cele; 
00438   delete[] fele; 
00439   delete[] xele; 
00440   delete[] yele; 
00441   delete[] zele;
00442   // zero out remaining tuples
00443   for(int n=nloop;n<nmax;n++){
00444     for(int i=0;i<((10-1)/nmax);i++) *tuple++ = 0;
00445   }
00446   //
00447   // fill the positron event variables
00448   //
00449   // get number of positron
00450   int npos = Event.GetNumPositron();
00451   //cout << npos << " positron(s)" << endl;
00452   *tuple++ = npos;
00453   //
00454   // deal with multiple positrons
00455   double* rpos = new double[npos];
00456   double* cpos = new double[npos];
00457   double* fpos = new double[npos];
00458   double* xpos = new double[npos];
00459   double* ypos = new double[npos];
00460   double* zpos = new double[npos];
00461   //
00462   // up to 1 positron per event, keep the first one
00463   //
00464   nmax = 1;
00465   if(npos < nmax)
00466     nloop = npos;
00467   else
00468     nloop = nmax;
00469   //
00470   for(int n=0;n<nloop;n++){
00471     *tuple++ = Event.GetPositronKE(n);
00472     *tuple++ = Event.GetPositronCosTheta(n);
00473     //
00474     rpos[n]  = Event.GetPositronVertexR(n);
00475     cpos[n]  = Event.GetPositronVertexCosTheta(n);
00476     fpos[n]  = Event.GetPositronVertexPhi(n);
00477     xpos[n]  = rpos[n]*cos(fpos[n])*sqrt(1-cpos[n]*cpos[n]);
00478     ypos[n]  = rpos[n]*sin(fpos[n])*sqrt(1-cpos[n]*cpos[n]);
00479     zpos[n]  = rpos[n]*cpos[n];
00480     //
00481     *tuple++ = rpos[n];
00482     *tuple++ = cpos[n];
00483     *tuple++ = fpos[n];
00484     *tuple++ = xpos[n];
00485     *tuple++ = ypos[n];
00486     *tuple++ = zpos[n];
00487     //
00488     *tuple++ = Event.GetPositronTime(n);
00489   }
00490   delete[] rpos; 
00491   delete[] cpos; 
00492   delete[] fpos; 
00493   delete[] xpos; 
00494   delete[] ypos;
00495   delete[] zpos;
00496   // zero out remaining tuples
00497   for(int n=nloop;n<nmax;n++){
00498     for(int i=0;i<((10-1)/nmax);i++) *tuple++ = 0;
00499   }
00500   //
00501   // fill the neutron event variables
00502   //
00503   // get number of neutron
00504   int nneu = Event.GetNumNeutron();
00505   //cout << nneu << " neutron(s)" << endl;
00506   *tuple++ = nneu;
00507   //
00508   // deal with multiple neutrons
00509   double* rneu = new double[nneu];
00510   double* cneu = new double[nneu];
00511   double* fneu = new double[nneu];
00512   double* xneu = new double[nneu];
00513   double* yneu = new double[nneu];
00514   double* zneu = new double[nneu];
00515   //
00516   // up to 3 neutrons per event, keep the first ones
00517   //
00518   nmax = 3;
00519   if(nneu < nmax)
00520     nloop = nneu;
00521   else
00522     nloop = nmax;
00523   //
00524   for(int n=0;n<nloop;n++){
00525     *tuple++ = Event.GetNeutronKE(n);
00526     *tuple++ = Event.GetNeutronCosTheta(n);
00527     *tuple++ = Event.GetNeutronParentProcess(n);
00528     //
00529     rneu[n]  = Event.GetNeutronVertexR(n);
00530     cneu[n]  = Event.GetNeutronVertexCosTheta(n);
00531     fneu[n]  = Event.GetNeutronVertexPhi(n);
00532     xneu[n]  = rneu[n]*cos(fneu[n])*sqrt(1-cneu[n]*cneu[n]);
00533     yneu[n]  = rneu[n]*sin(fneu[n])*sqrt(1-cneu[n]*cneu[n]);
00534     zneu[n]  = rneu[n]*cneu[n];
00535     //
00536     *tuple++ = rneu[n];
00537     *tuple++ = cneu[n];
00538     *tuple++ = fneu[n];
00539     *tuple++ = xneu[n];
00540     *tuple++ = yneu[n];
00541     *tuple++ = zneu[n];
00542     //
00543     *tuple++ = Event.GetNeutronTime(n);
00544     *tuple++ = Detector.GetNeutronSteps(n);
00545     *tuple++ = Detector.GetAbsorber(n);
00546   }
00547   delete[] rneu; 
00548   delete[] cneu; 
00549   delete[] fneu; 
00550   delete[] xneu; 
00551   delete[] yneu; 
00552   delete[] zneu;
00553   // zero out remaining tuples
00554   for(int n=nloop;n<nmax;n++){
00555     for(int i=0;i<((37-1)/nmax);i++) *tuple++ = 0;
00556   }
00557   //
00558   // fill the PMT detector variables
00559   //
00560   *tuple++ = Detector.GetNcosbins();
00561   *tuple++ = Detector.GetNphibins();
00562   //
00563   *tuple++ = Detector.GetPHpositron();
00564   *tuple++ = Detector.GetPHposReco();
00565   *tuple++ = Detector.GetPHneutron();
00566   *tuple++ = Detector.GetPHneuReco();
00567   *tuple++ = Detector.GetPHthallium();
00568   *tuple++ = Detector.GetPHthaReco();
00569   *tuple++ = Detector.GetPHbismuth();
00570   *tuple++ = Detector.GetPHbisReco();
00571   *tuple++ = Detector.GetPHcalifornium();
00572   *tuple++ = Detector.GetPHcalReco();
00573   *tuple++ = Detector.GetPHmuon();
00574   *tuple++ = Detector.GetPHmuoReco();
00575   *tuple++ = Detector.GetPHelectron();
00576   *tuple++ = Detector.GetPHeleReco();
00577   //
00578   *tuple++ = Detector.GetMeanNpePos();
00579   *tuple++ = Detector.GetMeanQPos();
00580   *tuple++ = Detector.GetMeanAttenPos();
00581   *tuple++ = Detector.GetMeanSfactorPos();
00582   *tuple++ = Detector.GetMeanNpeNeu();
00583   *tuple++ = Detector.GetMeanQNeu();
00584   *tuple++ = Detector.GetMeanAttenNeu();
00585   *tuple++ = Detector.GetMeanSfactorNeu();
00586   *tuple++ = Detector.GetMeanNpeTl();
00587   *tuple++ = Detector.GetMeanQTl();
00588   *tuple++ = Detector.GetMeanAttenTl();
00589   *tuple++ = Detector.GetMeanSfactorTl();
00590   *tuple++ = Detector.GetMeanNpeBi();
00591   *tuple++ = Detector.GetMeanQBi();
00592   *tuple++ = Detector.GetMeanAttenBi();
00593   *tuple++ = Detector.GetMeanSfactorBi();
00594   *tuple++ = Detector.GetMeanNpeCf();
00595   *tuple++ = Detector.GetMeanQCf();
00596   *tuple++ = Detector.GetMeanAttenCf();
00597   *tuple++ = Detector.GetMeanSfactorCf();
00598   *tuple++ = Detector.GetMeanNpeMuon();
00599   *tuple++ = Detector.GetMeanQMuon();
00600   *tuple++ = Detector.GetMeanAttenMuon();
00601   *tuple++ = Detector.GetMeanSfactorMuon();
00602   *tuple++ = Detector.GetMeanNpeEle();
00603   *tuple++ = Detector.GetMeanQEle();
00604   *tuple++ = Detector.GetMeanAttenEle();
00605   *tuple++ = Detector.GetMeanSfactorEle();
00606   //
00607   *tuple++ = Detector.GetNgamma("P"); 
00608   *tuple++ = Detector.GetNgamma("N");
00609   *tuple++ = Detector.GetNgamma("Tl");
00610   *tuple++ = Detector.GetNgamma("Bi");
00611   *tuple++ = Detector.GetNgamma("Cf");
00612   *tuple++ = Detector.GetNgamma("muon");
00613   *tuple++ = Detector.GetNgamma("E");
00614   //
00615   *tuple++ = Detector.GetHitFracPos();  
00616   *tuple++ = Detector.GetHitFracNeu();  
00617   *tuple++ = Detector.GetHitFracTl();  
00618   *tuple++ = Detector.GetHitFracBi();  
00619   *tuple++ = Detector.GetHitFracCf();  
00620   *tuple++ = Detector.GetHitFracMuon();  
00621   *tuple++ = Detector.GetHitFracEle();  
00622   //
00623   *tuple++ = Detector.GetHitPMTsPos();
00624   *tuple++ = Detector.GetHitPMTsNeu();
00625   *tuple++ = Detector.GetHitPMTsTha();  
00626   *tuple++ = Detector.GetHitPMTsBis();
00627   *tuple++ = Detector.GetHitPMTsCal();  
00628   *tuple++ = Detector.GetHitPMTsMuo();
00629   *tuple++ = Detector.GetHitPMTsEle();
00630   //
00631   hfn_(Id,Tuple);
00632   //
00633   // do some looping over all the "PMT"
00634   //
00635   int npmt=Detector.GetNpmt();
00636   for(int n=0;n<npmt;n++){
00637 
00638     int nele;
00639     int npos;
00640     int nneu;
00641     int ntha;
00642     int nbis;
00643     int ncal;
00644     int nmuo;
00645     double xyzpmt[3];
00646     for(int i=0;i<3;i++) xyzpmt[i] = 0.;
00647     double phele[100];
00648     double phpos[100];
00649     double phneu[100];
00650     double phtha[100];
00651     double phbis[100];
00652     double phcal[100];
00653     double phmuo[100];
00654     double tele[100];
00655     double tpos[100];
00656     double tneu[100];
00657     double ttha[100];
00658     double tbis[100];
00659     double tcal[100];
00660     double tmuo[100];
00661     int npepos[100];
00662     int npeneu[100];
00663     int npetha[100];
00664     int npebis[100];
00665     int npecal[100];
00666     int npemuo[100];
00667     int npeele[100];
00668     int thispmtnpeele = 0;
00669     int thispmtnpepos = 0;
00670     int thispmtnpeneu = 0;
00671     int thispmtnpetha = 0;
00672     int thispmtnpebis = 0;
00673     int thispmtnpecal = 0;
00674     int thispmtnpemuo = 0;
00675 
00676     for(int i=0;i<100;i++){
00677       phele[i] = 0.;
00678       phpos[i] = 0.;
00679       phneu[i] = 0.;
00680       phtha[i] = 0.;
00681       phbis[i] = 0.;
00682       phcal[i] = 0.;
00683       phmuo[i] = 0.;
00684       tele[i] = 0.;
00685       tpos[i] = 0.;
00686       tneu[i] = 0.;
00687       ttha[i] = 0.;
00688       tbis[i] = 0.;
00689       tcal[i] = 0.;
00690       tmuo[i] = 0.;
00691       npepos[i] = 0;
00692       npeneu[i] = 0;
00693       npetha[i] = 0;
00694       npebis[i] = 0;
00695       npecal[i] = 0;
00696       npemuo[i] = 0;
00697       npeele[i] = 0;
00698     }
00699     //
00700     // Get "phototube" information.
00701     //      
00702     Detector.GetPMTdata(n,xyzpmt,nele,phele,tele,
00703                         npos,phpos,tpos,nneu,phneu,tneu,
00704                         ntha,phtha,ttha,nbis,phbis,tbis,
00705                         ncal,phcal,tcal,nmuo,phmuo,tmuo,
00706                         npepos,npeneu,npetha,npebis,npecal,
00707                         npemuo,npeele);
00708 
00709     for(int nn=0;nn<nele;nn++){
00710       hfill_(Offset+101,*(phele+nn));
00711       hfill_(Offset+201,*(tele+nn));
00712       hfill_(Offset+1001,float(*(npeele+nn)));
00713       thispmtnpeele+=*(npeele+nn);
00714     }
00715     hfill_(Offset+2001,float(thispmtnpeele));
00716     for(int nn=0;nn<npos;nn++){
00717       hfill_(Offset+102,*(phpos+nn));
00718       hfill_(Offset+202,*(tpos+nn));
00719       hfill_(Offset+1002,float(*(npepos+nn)));
00720       thispmtnpepos+=*(npepos+nn);
00721     }
00722     hfill_(Offset+2002,float(thispmtnpepos));
00723     for(int nn=0;nn<nneu;nn++){
00724       hfill_(Offset+103,*(phneu+nn));
00725       hfill_(Offset+203,*(tneu+nn));
00726       hfill_(Offset+1003,float(*(npeneu+nn)));
00727       thispmtnpeneu+=*(npeneu+nn);
00728     }
00729     hfill_(Offset+2003,float(thispmtnpeneu));
00730     for(int nn=0;nn<ntha;nn++){
00731       hfill_(Offset+104,*(phtha+nn));
00732       hfill_(Offset+204,*(ttha+nn));
00733       hfill_(Offset+1004,float(*(npetha+nn)));
00734       thispmtnpetha+=*(npetha+nn);
00735     }
00736     hfill_(Offset+2004,float(thispmtnpetha));
00737     for(int nn=0;nn<nbis;nn++){
00738       hfill_(Offset+105,*(phbis+nn));
00739       hfill_(Offset+205,*(tbis+nn));
00740       hfill_(Offset+1005,float(*(npebis+nn)));
00741       thispmtnpebis+=*(npebis+nn);
00742     }
00743     hfill_(Offset+2005,float(thispmtnpebis));
00744     for(int nn=0;nn<ncal;nn++){
00745       hfill_(Offset+106,*(phcal+nn));
00746       hfill_(Offset+206,*(tcal+nn));
00747       hfill_(Offset+1006,float(*(npecal+nn)));
00748       thispmtnpecal+=*(npecal+nn);
00749     }
00750     hfill_(Offset+2006,float(thispmtnpecal));
00751     for(int nn=0;nn<nmuo;nn++){
00752       hfill_(Offset+107,*(phmuo+nn));
00753       hfill_(Offset+207,*(tmuo+nn));
00754       hfill_(Offset+1007,float(*(npemuo+nn)));
00755       thispmtnpemuo+=*(npemuo+nn);
00756     }
00757     hfill_(Offset+2007,float(thispmtnpemuo));
00758   }
00759   //
00760   // Get information about gammas from electron.
00761   //
00762   int nelegamma = Detector.GetNgamma("E"); 
00763   hfill_(Offset+301,float(nelegamma));
00764   double* ege = new double[nelegamma];
00765   double* tge = new double[nelegamma];
00766   double* rge = new double[nelegamma];
00767   double* cge = new double[nelegamma];
00768   double* fge = new double[nelegamma];
00769   Detector.GetGamma("E",ege,tge,rge,cge,fge);
00770   for (int n=0;n<nelegamma;n++){
00771     hfill_(Offset+401,ege[n]);
00772     hfill_(Offset+501,tge[n]);
00773     hfill_(Offset+601,rge[n]);
00774     hfill_(Offset+701,cge[n]);
00775   }
00776   delete[] ege;
00777   delete[] tge;
00778   delete[] rge;
00779   delete[] cge;
00780   delete[] fge;
00781   //
00782   // Get information about gammas from positron.
00783   // Note that the positron dedx is defined as a "gamma"
00784   // originating at the positron range point.
00785   //
00786   int nposgamma = Detector.GetNgamma("P"); 
00787   hfill_(Offset+302,float(nposgamma));
00788   double* egp = new double[nposgamma];
00789   double* tgp = new double[nposgamma];
00790   double* rgp = new double[nposgamma];
00791   double* cgp = new double[nposgamma];
00792   double* fgp = new double[nposgamma];
00793   Detector.GetGamma("P",egp,tgp,rgp,cgp,fgp);
00794   for (int n=0;n<nposgamma;n++){
00795     hfill_(Offset+402,egp[n]);
00796     hfill_(Offset+502,tgp[n]);
00797     hfill_(Offset+602,rgp[n]);
00798     hfill_(Offset+702,cgp[n]);
00799   }
00800   delete[] egp;
00801   delete[] tgp;
00802   delete[] rgp;
00803   delete[] cgp;
00804   delete[] fgp;
00805   //
00806   // Get information about gammas from neutron.
00807   //
00808   int nneugamma = Detector.GetNgamma("N"); 
00809   hfill_(Offset+303,float(nneugamma));
00810   double* egn = new double[nneugamma];
00811   double* tgn = new double[nneugamma];
00812   double* rgn = new double[nneugamma];
00813   double* cgn = new double[nneugamma];
00814   double* fgn = new double[nneugamma];
00815   Detector.GetGamma("N",egn,tgn,rgn,cgn,fgn);
00816   for (int n=0;n<nneugamma;n++){
00817     hfill_(Offset+403,egn[n]);
00818     hfill_(Offset+503,tgn[n]);
00819     hfill_(Offset+603,rgn[n]);
00820     hfill_(Offset+703,cgn[n]);
00821   }
00822   delete[] egn;
00823   delete[] tgn;
00824   delete[] rgn;
00825   delete[] cgn;
00826   delete[] fgn;
00827   //
00828   // Get information about prompt gammas from 252cf.
00829   //
00830   int ngamma = Detector.GetNgamma("Cf"); 
00831   hfill_(Offset+306,float(ngamma));
00832   double* egcf = new double[ngamma];
00833   double* tgcf = new double[ngamma];
00834   double* rgcf = new double[ngamma];
00835   double* cgcf = new double[ngamma];
00836   double* fgcf = new double[ngamma];
00837   Detector.GetGamma("Cf",egcf,tgcf,rgcf,cgcf,fgcf);
00838   for (int n=0;n<ngamma;n++){
00839     hfill_(Offset+406,egcf[n]);
00840     hfill_(Offset+506,tgcf[n]);
00841     hfill_(Offset+606,rgcf[n]);
00842     hfill_(Offset+706,cgcf[n]);
00843   }
00844   delete[] egcf;
00845   delete[] tgcf;
00846   delete[] rgcf;
00847   delete[] cgcf;
00848   delete[] fgcf;
00849 
00850   delete[] Tuple;
00851 }
00852
Generated on Mon Feb 21 16:11:19 2005 for ReactorFsim by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002