ReactorNtuple.cpp Source 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 = 8;
00041   int ntVarsMax = 127;
00042   int ntVars = 0;
00043   char ntTags[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",10,0,10,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[Nvar];
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[nnu];
00351   double cnu[nnu];
00352   double fnu[nnu];
00353   double xnu[nnu];
00354   double ynu[nnu];
00355   double znu[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   // zero out remaining tuples
00383   for(int n=nloop;n<nmax;n++){
00384     for(int i=0;i<(7/nmax);i++) *tuple++ = 0;
00385   }
00386   //
00387   // fill the electron event variables
00388   //
00389   // get number of electron
00390   int nele = Event.GetNumElectron();
00391   //cout << nele << " electron(s)" << endl;
00392   *tuple++ = nele;
00393   //
00394   // deal with multiple electrons
00395   double rele[nele];
00396   double cele[nele];
00397   double fele[nele];
00398   double xele[nele];
00399   double yele[nele];
00400   double zele[nele];
00401   //
00402   // up to 1 electron per event, keep the first one
00403   //
00404   nmax = 1;
00405   if(nele < nmax)
00406     nloop = nele;
00407   else
00408     nloop = nmax;
00409   //
00410   for(int n=0;n<nloop;n++){
00411     *tuple++ = Event.GetElectronKE(n);
00412     *tuple++ = Event.GetElectronCosTheta(n);
00413     //
00414     rele[n]  = Event.GetElectronVertexR(n);
00415     cele[n]  = Event.GetElectronVertexCosTheta(n);
00416     fele[n]  = Event.GetElectronVertexPhi(n);
00417     xele[n]  = rele[n]*cos(fele[n])*sqrt(1-cele[n]*cele[n]);
00418     yele[n]  = rele[n]*sin(fele[n])*sqrt(1-cele[n]*cele[n]);
00419     zele[n]  = rele[n]*cele[n];
00420     //
00421     *tuple++ = rele[n];
00422     *tuple++ = cele[n];
00423     *tuple++ = fele[n];
00424     *tuple++ = xele[n];
00425     *tuple++ = yele[n];
00426     *tuple++ = zele[n];
00427     //
00428     *tuple++ = Event.GetElectronTime(n);
00429   }
00430   // zero out remaining tuples
00431   for(int n=nloop;n<nmax;n++){
00432     for(int i=0;i<((10-1)/nmax);i++) *tuple++ = 0;
00433   }
00434   //
00435   // fill the positron event variables
00436   //
00437   // get number of positron
00438   int npos = Event.GetNumPositron();
00439   //cout << npos << " positron(s)" << endl;
00440   *tuple++ = npos;
00441   //
00442   // deal with multiple positrons
00443   double rpos[npos];
00444   double cpos[npos];
00445   double fpos[npos];
00446   double xpos[npos];
00447   double ypos[npos];
00448   double zpos[npos];
00449   //
00450   // up to 1 positron per event, keep the first one
00451   //
00452   nmax = 1;
00453   if(npos < nmax)
00454     nloop = npos;
00455   else
00456     nloop = nmax;
00457   //
00458   for(int n=0;n<nloop;n++){
00459     *tuple++ = Event.GetPositronKE(n);
00460     *tuple++ = Event.GetPositronCosTheta(n);
00461     //
00462     rpos[n]  = Event.GetPositronVertexR(n);
00463     cpos[n]  = Event.GetPositronVertexCosTheta(n);
00464     fpos[n]  = Event.GetPositronVertexPhi(n);
00465     xpos[n]  = rpos[n]*cos(fpos[n])*sqrt(1-cpos[n]*cpos[n]);
00466     ypos[n]  = rpos[n]*sin(fpos[n])*sqrt(1-cpos[n]*cpos[n]);
00467     zpos[n]  = rpos[n]*cpos[n];
00468     //
00469     *tuple++ = rpos[n];
00470     *tuple++ = cpos[n];
00471     *tuple++ = fpos[n];
00472     *tuple++ = xpos[n];
00473     *tuple++ = ypos[n];
00474     *tuple++ = zpos[n];
00475     //
00476     *tuple++ = Event.GetPositronTime(n);
00477   }
00478   // zero out remaining tuples
00479   for(int n=nloop;n<nmax;n++){
00480     for(int i=0;i<((10-1)/nmax);i++) *tuple++ = 0;
00481   }
00482   //
00483   // fill the neutron event variables
00484   //
00485   // get number of neutron
00486   int nneu = Event.GetNumNeutron();
00487   //cout << nneu << " neutron(s)" << endl;
00488   *tuple++ = nneu;
00489   //
00490   // deal with multiple neutrons
00491   double rneu[nneu];
00492   double cneu[nneu];
00493   double fneu[nneu];
00494   double xneu[nneu];
00495   double yneu[nneu];
00496   double zneu[nneu];
00497   //
00498   // up to 3 neutrons per event, keep the first ones
00499   //
00500   nmax = 3;
00501   if(nneu < nmax)
00502     nloop = nneu;
00503   else
00504     nloop = nmax;
00505   //
00506   for(int n=0;n<nloop;n++){
00507     *tuple++ = Event.GetNeutronKE(n);
00508     *tuple++ = Event.GetNeutronCosTheta(n);
00509     *tuple++ = Event.GetNeutronParentProcess(n);
00510     //
00511     rneu[n]  = Event.GetNeutronVertexR(n);
00512     cneu[n]  = Event.GetNeutronVertexCosTheta(n);
00513     fneu[n]  = Event.GetNeutronVertexPhi(n);
00514     xneu[n]  = rneu[n]*cos(fneu[n])*sqrt(1-cneu[n]*cneu[n]);
00515     yneu[n]  = rneu[n]*sin(fneu[n])*sqrt(1-cneu[n]*cneu[n]);
00516     zneu[n]  = rneu[n]*cneu[n];
00517     //
00518     *tuple++ = rneu[n];
00519     *tuple++ = cneu[n];
00520     *tuple++ = fneu[n];
00521     *tuple++ = xneu[n];
00522     *tuple++ = yneu[n];
00523     *tuple++ = zneu[n];
00524     //
00525     *tuple++ = Event.GetNeutronTime(n);
00526     *tuple++ = Detector.GetNeutronSteps(n);
00527     *tuple++ = Detector.GetAbsorber(n);
00528   }
00529   // zero out remaining tuples
00530   for(int n=nloop;n<nmax;n++){
00531     for(int i=0;i<((37-1)/nmax);i++) *tuple++ = 0;
00532   }
00533   //
00534   // fill the PMT detector variables
00535   //
00536   *tuple++ = Detector.GetNcosbins();
00537   *tuple++ = Detector.GetNphibins();
00538   //
00539   *tuple++ = Detector.GetPHpositron();
00540   *tuple++ = Detector.GetPHposReco();
00541   *tuple++ = Detector.GetPHneutron();
00542   *tuple++ = Detector.GetPHneuReco();
00543   *tuple++ = Detector.GetPHthallium();
00544   *tuple++ = Detector.GetPHthaReco();
00545   *tuple++ = Detector.GetPHbismuth();
00546   *tuple++ = Detector.GetPHbisReco();
00547   *tuple++ = Detector.GetPHcalifornium();
00548   *tuple++ = Detector.GetPHcalReco();
00549   *tuple++ = Detector.GetPHmuon();
00550   *tuple++ = Detector.GetPHmuoReco();
00551   *tuple++ = Detector.GetPHelectron();
00552   *tuple++ = Detector.GetPHeleReco();
00553   //
00554   *tuple++ = Detector.GetMeanNpePos();
00555   *tuple++ = Detector.GetMeanQPos();
00556   *tuple++ = Detector.GetMeanAttenPos();
00557   *tuple++ = Detector.GetMeanSfactorPos();
00558   *tuple++ = Detector.GetMeanNpeNeu();
00559   *tuple++ = Detector.GetMeanQNeu();
00560   *tuple++ = Detector.GetMeanAttenNeu();
00561   *tuple++ = Detector.GetMeanSfactorNeu();
00562   *tuple++ = Detector.GetMeanNpeTl();
00563   *tuple++ = Detector.GetMeanQTl();
00564   *tuple++ = Detector.GetMeanAttenTl();
00565   *tuple++ = Detector.GetMeanSfactorTl();
00566   *tuple++ = Detector.GetMeanNpeBi();
00567   *tuple++ = Detector.GetMeanQBi();
00568   *tuple++ = Detector.GetMeanAttenBi();
00569   *tuple++ = Detector.GetMeanSfactorBi();
00570   *tuple++ = Detector.GetMeanNpeCf();
00571   *tuple++ = Detector.GetMeanQCf();
00572   *tuple++ = Detector.GetMeanAttenCf();
00573   *tuple++ = Detector.GetMeanSfactorCf();
00574   *tuple++ = Detector.GetMeanNpeMuon();
00575   *tuple++ = Detector.GetMeanQMuon();
00576   *tuple++ = Detector.GetMeanAttenMuon();
00577   *tuple++ = Detector.GetMeanSfactorMuon();
00578   *tuple++ = Detector.GetMeanNpeEle();
00579   *tuple++ = Detector.GetMeanQEle();
00580   *tuple++ = Detector.GetMeanAttenEle();
00581   *tuple++ = Detector.GetMeanSfactorEle();
00582   //
00583   *tuple++ = Detector.GetNgamma("P"); 
00584   *tuple++ = Detector.GetNgamma("N");
00585   *tuple++ = Detector.GetNgamma("Tl");
00586   *tuple++ = Detector.GetNgamma("Bi");
00587   *tuple++ = Detector.GetNgamma("Cf");
00588   *tuple++ = Detector.GetNgamma("muon");
00589   *tuple++ = Detector.GetNgamma("E");
00590   //
00591   *tuple++ = Detector.GetHitFracPos();  
00592   *tuple++ = Detector.GetHitFracNeu();  
00593   *tuple++ = Detector.GetHitFracTl();  
00594   *tuple++ = Detector.GetHitFracBi();  
00595   *tuple++ = Detector.GetHitFracCf();  
00596   *tuple++ = Detector.GetHitFracMuon();  
00597   *tuple++ = Detector.GetHitFracEle();  
00598   //
00599   *tuple++ = Detector.GetHitPMTsPos();
00600   *tuple++ = Detector.GetHitPMTsNeu();
00601   *tuple++ = Detector.GetHitPMTsTha();  
00602   *tuple++ = Detector.GetHitPMTsBis();
00603   *tuple++ = Detector.GetHitPMTsCal();  
00604   *tuple++ = Detector.GetHitPMTsMuo();
00605   *tuple++ = Detector.GetHitPMTsEle();
00606   //
00607   hfn_(Id,Tuple);
00608   //
00609   // do some looping over all the "PMT"
00610   //
00611   int npmt=Detector.GetNpmt();
00612   for(int n=0;n<npmt;n++){
00613 
00614     int nele;
00615     int npos;
00616     int nneu;
00617     int ntha;
00618     int nbis;
00619     int ncal;
00620     int nmuo;
00621     double xyzpmt[3];
00622     for(int i=0;i<3;i++) xyzpmt[i] = 0.;
00623     double phele[100];
00624     double phpos[100];
00625     double phneu[100];
00626     double phtha[100];
00627     double phbis[100];
00628     double phcal[100];
00629     double phmuo[100];
00630     double tele[100];
00631     double tpos[100];
00632     double tneu[100];
00633     double ttha[100];
00634     double tbis[100];
00635     double tcal[100];
00636     double tmuo[100];
00637     int npepos[100];
00638     int npeneu[100];
00639     int npetha[100];
00640     int npebis[100];
00641     int npecal[100];
00642     int npemuo[100];
00643     int npeele[100];
00644     int thispmtnpeele = 0;
00645     int thispmtnpepos = 0;
00646     int thispmtnpeneu = 0;
00647     int thispmtnpetha = 0;
00648     int thispmtnpebis = 0;
00649     int thispmtnpecal = 0;
00650     int thispmtnpemuo = 0;
00651 
00652     for(int i=0;i<100;i++){
00653       phele[i] = 0.;
00654       phpos[i] = 0.;
00655       phneu[i] = 0.;
00656       phtha[i] = 0.;
00657       phbis[i] = 0.;
00658       phcal[i] = 0.;
00659       phmuo[i] = 0.;
00660       tele[i] = 0.;
00661       tpos[i] = 0.;
00662       tneu[i] = 0.;
00663       ttha[i] = 0.;
00664       tbis[i] = 0.;
00665       tcal[i] = 0.;
00666       tmuo[i] = 0.;
00667       npepos[i] = 0;
00668       npeneu[i] = 0;
00669       npetha[i] = 0;
00670       npebis[i] = 0;
00671       npecal[i] = 0;
00672       npemuo[i] = 0;
00673       npeele[i] = 0;
00674     }
00675     //
00676     // Get "phototube" information.
00677     //      
00678     Detector.GetPMTdata(n,xyzpmt,nele,phele,tele,
00679                         npos,phpos,tpos,nneu,phneu,tneu,
00680                         ntha,phtha,ttha,nbis,phbis,tbis,
00681                         ncal,phcal,tcal,nmuo,phmuo,tmuo,
00682                         npepos,npeneu,npetha,npebis,npecal,
00683                         npemuo,npeele);
00684 
00685     for(int nn=0;nn<nele;nn++){
00686       hfill_(Offset+101,*(phele+nn));
00687       hfill_(Offset+201,*(tele+nn));
00688       hfill_(Offset+1001,float(*(npeele+nn)));
00689       thispmtnpeele+=*(npeele+nn);
00690     }
00691     hfill_(Offset+2001,float(thispmtnpeele));
00692     for(int nn=0;nn<npos;nn++){
00693       hfill_(Offset+102,*(phpos+nn));
00694       hfill_(Offset+202,*(tpos+nn));
00695       hfill_(Offset+1002,float(*(npepos+nn)));
00696       thispmtnpepos+=*(npepos+nn);
00697     }
00698     hfill_(Offset+2002,float(thispmtnpepos));
00699     for(int nn=0;nn<nneu;nn++){
00700       hfill_(Offset+103,*(phneu+nn));
00701       hfill_(Offset+203,*(tneu+nn));
00702       hfill_(Offset+1003,float(*(npeneu+nn)));
00703       thispmtnpeneu+=*(npeneu+nn);
00704     }
00705     hfill_(Offset+2003,float(thispmtnpeneu));
00706     for(int nn=0;nn<ntha;nn++){
00707       hfill_(Offset+104,*(phtha+nn));
00708       hfill_(Offset+204,*(ttha+nn));
00709       hfill_(Offset+1004,float(*(npetha+nn)));
00710       thispmtnpetha+=*(npetha+nn);
00711     }
00712     hfill_(Offset+2004,float(thispmtnpetha));
00713     for(int nn=0;nn<nbis;nn++){
00714       hfill_(Offset+105,*(phbis+nn));
00715       hfill_(Offset+205,*(tbis+nn));
00716       hfill_(Offset+1005,float(*(npebis+nn)));
00717       thispmtnpebis+=*(npebis+nn);
00718     }
00719     hfill_(Offset+2005,float(thispmtnpebis));
00720     for(int nn=0;nn<ncal;nn++){
00721       hfill_(Offset+106,*(phcal+nn));
00722       hfill_(Offset+206,*(tcal+nn));
00723       hfill_(Offset+1006,float(*(npecal+nn)));
00724       thispmtnpecal+=*(npecal+nn);
00725     }
00726     hfill_(Offset+2006,float(thispmtnpecal));
00727     for(int nn=0;nn<nmuo;nn++){
00728       hfill_(Offset+107,*(phmuo+nn));
00729       hfill_(Offset+207,*(tmuo+nn));
00730       hfill_(Offset+1007,float(*(npemuo+nn)));
00731       thispmtnpemuo+=*(npemuo+nn);
00732     }
00733     hfill_(Offset+2007,float(thispmtnpemuo));
00734   }
00735   //
00736   // Get information about gammas from electron.
00737   //
00738   int nelegamma = Detector.GetNgamma("E"); 
00739   hfill_(Offset+301,float(nelegamma));
00740   double ege[nelegamma];
00741   double tge[nelegamma];
00742   double rge[nelegamma];
00743   double cge[nelegamma];
00744   double fge[nelegamma];
00745   Detector.GetGamma("E",ege,tge,rge,cge,fge);
00746   for (int n=0;n<nelegamma;n++){
00747     hfill_(Offset+401,ege[n]);
00748     hfill_(Offset+501,tge[n]);
00749     hfill_(Offset+601,rge[n]);
00750     hfill_(Offset+701,cge[n]);
00751   }
00752   //
00753   // Get information about gammas from positron.
00754   // Note that the positron dedx is defined as a "gamma"
00755   // originating at the positron range point.
00756   //
00757   int nposgamma = Detector.GetNgamma("P"); 
00758   hfill_(Offset+302,float(nposgamma));
00759   double egp[nposgamma];
00760   double tgp[nposgamma];
00761   double rgp[nposgamma];
00762   double cgp[nposgamma];
00763   double fgp[nposgamma];
00764   Detector.GetGamma("P",egp,tgp,rgp,cgp,fgp);
00765   for (int n=0;n<nposgamma;n++){
00766     hfill_(Offset+402,egp[n]);
00767     hfill_(Offset+502,tgp[n]);
00768     hfill_(Offset+602,rgp[n]);
00769     hfill_(Offset+702,cgp[n]);
00770   }
00771   //
00772   // Get information about gammas from neutron.
00773   //
00774   int nneugamma = Detector.GetNgamma("N"); 
00775   hfill_(Offset+303,float(nneugamma));
00776   double egn[nneugamma];
00777   double tgn[nneugamma];
00778   double rgn[nneugamma];
00779   double cgn[nneugamma];
00780   double fgn[nneugamma];
00781   Detector.GetGamma("N",egn,tgn,rgn,cgn,fgn);
00782   for (int n=0;n<nneugamma;n++){
00783     hfill_(Offset+403,egn[n]);
00784     hfill_(Offset+503,tgn[n]);
00785     hfill_(Offset+603,rgn[n]);
00786     hfill_(Offset+703,cgn[n]);
00787   }
00788   //
00789   // Get information about prompt gammas from 252cf.
00790   //
00791   int ngamma = Detector.GetNgamma("Cf"); 
00792   hfill_(Offset+306,float(ngamma));
00793   double egcf[ngamma];
00794   double tgcf[ngamma];
00795   double rgcf[ngamma];
00796   double cgcf[ngamma];
00797   double fgcf[ngamma];
00798   Detector.GetGamma("Cf",egcf,tgcf,rgcf,cgcf,fgcf);
00799   for (int n=0;n<ngamma;n++){
00800     hfill_(Offset+406,egcf[n]);
00801     hfill_(Offset+506,tgcf[n]);
00802     hfill_(Offset+606,rgcf[n]);
00803     hfill_(Offset+706,cgcf[n]);
00804   }
00805 }