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   IdGen = Ntuple.IdGen;
00020   IdNu  = Ntuple.IdNu;
00021   IdEle = Ntuple.IdEle;
00022   IdPos = Ntuple.IdPos;
00023   IdNeu = Ntuple.IdNeu;
00024   IdDet = Ntuple.IdDet;
00025 
00026   NvarGen = Ntuple.NvarGen;
00027   NvarNu  = Ntuple.NvarNu;
00028   NvarEle = Ntuple.NvarEle;
00029   NvarPos = Ntuple.NvarPos;
00030   NvarNeu = Ntuple.NvarNeu;
00031   NvarDet = Ntuple.NvarDet;
00032 
00033   Offset = Ntuple.Offset;
00034 }
00035 ReactorNtuple& ReactorNtuple::operator=(const ReactorNtuple& rhs){
00036 
00037   if(this == &rhs)return *this;
00038   IdGen = rhs.IdGen;
00039   IdNu  = rhs.IdNu;
00040   IdEle = rhs.IdEle;
00041   IdPos = rhs.IdPos;
00042   IdNeu = rhs.IdNeu;
00043   IdDet = rhs.IdDet;
00044 
00045   NvarGen = rhs.NvarGen;
00046   NvarNu  = rhs.NvarNu;
00047   NvarEle = rhs.NvarEle;
00048   NvarPos = rhs.NvarPos;
00049   NvarNeu = rhs.NvarNeu;
00050   NvarDet = rhs.NvarDet;
00051 
00052   Offset = rhs.Offset;
00053   return *this;
00054 }
00055 
00056 ReactorNtuple::ReactorNtuple(int ntId,
00057                              const char* ntName,
00058                              char* dirName){
00059 
00060   // define ntuple tags here
00061   // to add vars, increment ntVars, and add tag to end of list
00062   // don't change wordsize
00063   //
00064   int wordsize = 8;
00065   int ntVarsMax = 100;
00066   char ntTags[ntVarsMax*wordsize];
00067   //
00068   int ntVars = 0;
00069   int id = 1;
00070   char Name[100];
00071   //
00072   // book the general event veriables, nt=1001
00073   //
00074   strcpy(ntTags+ntVars++*wordsize,"nevent   ");
00075   strcpy(ntTags+ntVars++*wordsize,"wtc      ");
00076   strcpy(ntTags+ntVars++*wordsize,"wtf      ");
00077   //
00078   IdGen = 1000*ntId+id++;
00079   sprintf(Name,"%s%s",ntName,", general");
00080   NvarGen = ntVars;
00081   //
00082   hbookn_(IdGen,Name,NvarGen,dirName,1000,ntTags,
00083           strlen(Name),strlen(dirName),wordsize,NvarGen);
00084   //
00085   // reset
00086   ntVars = 0;
00087   //
00088   // book the nubar event variables, nt=1002
00089   //
00090   // up to 1 nubar per event
00091   strcpy(ntTags+ntVars++*wordsize,"nnu      ");
00092 
00093   strcpy(ntTags+ntVars++*wordsize,"enu1     ");
00094   strcpy(ntTags+ntVars++*wordsize,"nuproc1  ");
00095   strcpy(ntTags+ntVars++*wordsize,"vnu1r    ");
00096   strcpy(ntTags+ntVars++*wordsize,"vnu1cos  ");
00097   strcpy(ntTags+ntVars++*wordsize,"vnu1phi  ");
00098   strcpy(ntTags+ntVars++*wordsize,"vnu1x    ");
00099   strcpy(ntTags+ntVars++*wordsize,"vnu1y    ");
00100   strcpy(ntTags+ntVars++*wordsize,"vnu1z    ");
00101   //
00102   IdNu = 1000*ntId+id++;
00103   sprintf(Name,"%s%s",ntName,", nubar");
00104   NvarNu = ntVars;
00105   //
00106   hbookn_(IdNu,Name,NvarNu,dirName,1000,ntTags,
00107           strlen(Name),strlen(dirName),wordsize,NvarNu);
00108   //
00109   // reset
00110   ntVars = 0;
00111   //
00112   // book the electron event variables, nt=1003
00113   //
00114   // up to 1 electron per event
00115   strcpy(ntTags+ntVars++*wordsize,"nele     ");
00116 
00117   strcpy(ntTags+ntVars++*wordsize,"KEele1   ");
00118   strcpy(ntTags+ntVars++*wordsize,"ele1cos  ");
00119   strcpy(ntTags+ntVars++*wordsize,"eleproc1 ");
00120   strcpy(ntTags+ntVars++*wordsize,"gele1r   ");
00121   strcpy(ntTags+ntVars++*wordsize,"gele1cos ");
00122   strcpy(ntTags+ntVars++*wordsize,"gele1phi ");
00123   strcpy(ntTags+ntVars++*wordsize,"gele1x   ");
00124   strcpy(ntTags+ntVars++*wordsize,"gele1y   ");
00125   strcpy(ntTags+ntVars++*wordsize,"gele1z   ");
00126   strcpy(ntTags+ntVars++*wordsize,"vele1r   ");
00127   strcpy(ntTags+ntVars++*wordsize,"vele1cos ");
00128   strcpy(ntTags+ntVars++*wordsize,"vele1phi ");
00129   strcpy(ntTags+ntVars++*wordsize,"vele1x   ");
00130   strcpy(ntTags+ntVars++*wordsize,"vele1y   ");
00131   strcpy(ntTags+ntVars++*wordsize,"vele1z   ");
00132   strcpy(ntTags+ntVars++*wordsize,"tele1    ");
00133   //
00134   IdEle = 1000*ntId+id++;
00135   sprintf(Name,"%s%s",ntName,", electron");
00136   NvarEle = ntVars;
00137   //
00138   hbookn_(IdEle,Name,NvarEle,dirName,1000,ntTags,
00139           strlen(Name),strlen(dirName),wordsize,NvarEle);
00140   //
00141   // reset
00142   ntVars = 0;
00143   //
00144   // book the positron event variables, nt=1004
00145   //
00146   // up to 1 positron per event
00147   strcpy(ntTags+ntVars++*wordsize,"npos     ");
00148 
00149   strcpy(ntTags+ntVars++*wordsize,"KEpos1   ");
00150   strcpy(ntTags+ntVars++*wordsize,"pos1cos  ");
00151   strcpy(ntTags+ntVars++*wordsize,"posproc1 ");
00152   strcpy(ntTags+ntVars++*wordsize,"gpos1r   ");
00153   strcpy(ntTags+ntVars++*wordsize,"gpos1cos ");
00154   strcpy(ntTags+ntVars++*wordsize,"gpos1phi ");
00155   strcpy(ntTags+ntVars++*wordsize,"gpos1x   ");
00156   strcpy(ntTags+ntVars++*wordsize,"gpos1y   ");
00157   strcpy(ntTags+ntVars++*wordsize,"gpos1z   ");
00158   strcpy(ntTags+ntVars++*wordsize,"vpos1r   ");
00159   strcpy(ntTags+ntVars++*wordsize,"vpos1cos ");
00160   strcpy(ntTags+ntVars++*wordsize,"vpos1phi ");
00161   strcpy(ntTags+ntVars++*wordsize,"vpos1x   ");
00162   strcpy(ntTags+ntVars++*wordsize,"vpos1y   ");
00163   strcpy(ntTags+ntVars++*wordsize,"vpos1z   ");
00164   strcpy(ntTags+ntVars++*wordsize,"tpos1    ");
00165   //
00166   IdPos = 1000*ntId+id++;
00167   sprintf(Name,"%s%s",ntName,", positron");
00168   NvarPos = ntVars;
00169   //
00170   hbookn_(IdPos,Name,NvarPos,dirName,1000,ntTags,
00171           strlen(Name),strlen(dirName),wordsize,NvarPos);
00172   //
00173   // reset
00174   ntVars = 0;
00175   //
00176   // book the neutron event variables, nt=1005
00177   //
00178   // up to 5 neutron per event
00179   strcpy(ntTags+ntVars++*wordsize,"nneu     ");
00180 
00181   strcpy(ntTags+ntVars++*wordsize,"KEneu1   ");
00182   strcpy(ntTags+ntVars++*wordsize,"neu1cos  ");
00183   strcpy(ntTags+ntVars++*wordsize,"neuproc1 ");
00184   strcpy(ntTags+ntVars++*wordsize,"gneu1r   ");
00185   strcpy(ntTags+ntVars++*wordsize,"gneu1cos ");
00186   strcpy(ntTags+ntVars++*wordsize,"gneu1phi ");
00187   strcpy(ntTags+ntVars++*wordsize,"gneu1x   ");
00188   strcpy(ntTags+ntVars++*wordsize,"gneu1y   ");
00189   strcpy(ntTags+ntVars++*wordsize,"gneu1z   ");
00190   strcpy(ntTags+ntVars++*wordsize,"vneu1r   ");
00191   strcpy(ntTags+ntVars++*wordsize,"vneu1cos ");
00192   strcpy(ntTags+ntVars++*wordsize,"vneu1phi ");
00193   strcpy(ntTags+ntVars++*wordsize,"vneu1x   ");
00194   strcpy(ntTags+ntVars++*wordsize,"vneu1y   ");
00195   strcpy(ntTags+ntVars++*wordsize,"vneu1z   ");
00196   strcpy(ntTags+ntVars++*wordsize,"tneu1    ");
00197   strcpy(ntTags+ntVars++*wordsize,"nsteps1  ");
00198   strcpy(ntTags+ntVars++*wordsize,"absorb1  ");
00199 
00200   strcpy(ntTags+ntVars++*wordsize,"KEneu2   ");
00201   strcpy(ntTags+ntVars++*wordsize,"neu2cos  ");
00202   strcpy(ntTags+ntVars++*wordsize,"neuproc2 ");
00203   strcpy(ntTags+ntVars++*wordsize,"gneu2r   ");
00204   strcpy(ntTags+ntVars++*wordsize,"gneu2cos ");
00205   strcpy(ntTags+ntVars++*wordsize,"gneu2phi ");
00206   strcpy(ntTags+ntVars++*wordsize,"gneu2x   ");
00207   strcpy(ntTags+ntVars++*wordsize,"gneu2y   ");
00208   strcpy(ntTags+ntVars++*wordsize,"gneu2z   ");
00209   strcpy(ntTags+ntVars++*wordsize,"vneu2r   ");
00210   strcpy(ntTags+ntVars++*wordsize,"vneu2cos ");
00211   strcpy(ntTags+ntVars++*wordsize,"vneu2phi ");
00212   strcpy(ntTags+ntVars++*wordsize,"vneu2x   ");
00213   strcpy(ntTags+ntVars++*wordsize,"vneu2y   ");
00214   strcpy(ntTags+ntVars++*wordsize,"vneu2z   ");
00215   strcpy(ntTags+ntVars++*wordsize,"tneu2    ");
00216   strcpy(ntTags+ntVars++*wordsize,"nsteps2  ");
00217   strcpy(ntTags+ntVars++*wordsize,"absorb2  ");
00218 
00219   strcpy(ntTags+ntVars++*wordsize,"KEneu3   ");
00220   strcpy(ntTags+ntVars++*wordsize,"neu3cos  ");
00221   strcpy(ntTags+ntVars++*wordsize,"neuproc3 ");
00222   strcpy(ntTags+ntVars++*wordsize,"gneu3r   ");
00223   strcpy(ntTags+ntVars++*wordsize,"gneu3cos ");
00224   strcpy(ntTags+ntVars++*wordsize,"gneu3phi ");
00225   strcpy(ntTags+ntVars++*wordsize,"gneu3x   ");
00226   strcpy(ntTags+ntVars++*wordsize,"gneu3y   ");
00227   strcpy(ntTags+ntVars++*wordsize,"gneu3z   ");
00228   strcpy(ntTags+ntVars++*wordsize,"vneu3r   ");
00229   strcpy(ntTags+ntVars++*wordsize,"vneu3cos ");
00230   strcpy(ntTags+ntVars++*wordsize,"vneu3phi ");
00231   strcpy(ntTags+ntVars++*wordsize,"vneu3x   ");
00232   strcpy(ntTags+ntVars++*wordsize,"vneu3y   ");
00233   strcpy(ntTags+ntVars++*wordsize,"vneu3z   ");
00234   strcpy(ntTags+ntVars++*wordsize,"tneu3    ");
00235   strcpy(ntTags+ntVars++*wordsize,"nsteps3  ");
00236   strcpy(ntTags+ntVars++*wordsize,"absorb3  ");
00237 
00238   strcpy(ntTags+ntVars++*wordsize,"KEneu4   ");
00239   strcpy(ntTags+ntVars++*wordsize,"neu4cos  ");
00240   strcpy(ntTags+ntVars++*wordsize,"neuproc4 ");
00241   strcpy(ntTags+ntVars++*wordsize,"gneu4r   ");
00242   strcpy(ntTags+ntVars++*wordsize,"gneu4cos ");
00243   strcpy(ntTags+ntVars++*wordsize,"gneu4phi ");
00244   strcpy(ntTags+ntVars++*wordsize,"gneu4x   ");
00245   strcpy(ntTags+ntVars++*wordsize,"gneu4y   ");
00246   strcpy(ntTags+ntVars++*wordsize,"gneu4z   ");
00247   strcpy(ntTags+ntVars++*wordsize,"vneu4r   ");
00248   strcpy(ntTags+ntVars++*wordsize,"vneu4cos ");
00249   strcpy(ntTags+ntVars++*wordsize,"vneu4phi ");
00250   strcpy(ntTags+ntVars++*wordsize,"vneu4x   ");
00251   strcpy(ntTags+ntVars++*wordsize,"vneu4y   ");
00252   strcpy(ntTags+ntVars++*wordsize,"vneu4z   ");
00253   strcpy(ntTags+ntVars++*wordsize,"tneu4    ");
00254   strcpy(ntTags+ntVars++*wordsize,"nsteps4  ");
00255   strcpy(ntTags+ntVars++*wordsize,"absorb4  ");
00256 
00257   strcpy(ntTags+ntVars++*wordsize,"KEneu5   ");
00258   strcpy(ntTags+ntVars++*wordsize,"neu5cos  ");
00259   strcpy(ntTags+ntVars++*wordsize,"neuproc5 ");
00260   strcpy(ntTags+ntVars++*wordsize,"gneu5r   ");
00261   strcpy(ntTags+ntVars++*wordsize,"gneu5cos ");
00262   strcpy(ntTags+ntVars++*wordsize,"gneu5phi ");
00263   strcpy(ntTags+ntVars++*wordsize,"gneu5x   ");
00264   strcpy(ntTags+ntVars++*wordsize,"gneu5y   ");
00265   strcpy(ntTags+ntVars++*wordsize,"gneu5z   ");
00266   strcpy(ntTags+ntVars++*wordsize,"vneu5r   ");
00267   strcpy(ntTags+ntVars++*wordsize,"vneu5cos ");
00268   strcpy(ntTags+ntVars++*wordsize,"vneu5phi ");
00269   strcpy(ntTags+ntVars++*wordsize,"vneu5x   ");
00270   strcpy(ntTags+ntVars++*wordsize,"vneu5y   ");
00271   strcpy(ntTags+ntVars++*wordsize,"vneu5z   ");
00272   strcpy(ntTags+ntVars++*wordsize,"tneu5    ");
00273   strcpy(ntTags+ntVars++*wordsize,"nsteps5  ");
00274   strcpy(ntTags+ntVars++*wordsize,"absorb5  ");
00275   //
00276   IdNeu = 1000*ntId+id++;
00277   sprintf(Name,"%s%s",ntName,", neutron");
00278   NvarNeu = ntVars;
00279   //
00280   hbookn_(IdNeu,Name,NvarNeu,dirName,1000,ntTags,
00281           strlen(Name),strlen(dirName),wordsize,NvarNeu);
00282   //
00283   // reset
00284   ntVars = 0;
00285   //
00286   // book the pmt detector variables, nt=1006
00287   //
00288   strcpy(ntTags+ntVars++*wordsize,"ncosbins ");
00289   strcpy(ntTags+ntVars++*wordsize,"nphibins ");
00290 
00291   strcpy(ntTags+ntVars++*wordsize,"phposvis ");
00292   strcpy(ntTags+ntVars++*wordsize,"phposrec ");
00293   strcpy(ntTags+ntVars++*wordsize,"phneuvis ");
00294   strcpy(ntTags+ntVars++*wordsize,"phneurec ");
00295   strcpy(ntTags+ntVars++*wordsize,"phtlvis  ");
00296   strcpy(ntTags+ntVars++*wordsize,"phtlrec  ");
00297   strcpy(ntTags+ntVars++*wordsize,"phbivis  ");
00298   strcpy(ntTags+ntVars++*wordsize,"phbirec  ");
00299   strcpy(ntTags+ntVars++*wordsize,"phcfvis  ");
00300   strcpy(ntTags+ntVars++*wordsize,"phcfrec  ");
00301   strcpy(ntTags+ntVars++*wordsize,"phmuovis ");
00302   strcpy(ntTags+ntVars++*wordsize,"phmuorec ");
00303   strcpy(ntTags+ntVars++*wordsize,"phelevis ");
00304   strcpy(ntTags+ntVars++*wordsize,"phelerec ");
00305 
00306   strcpy(ntTags+ntVars++*wordsize,"npepos   ");
00307   strcpy(ntTags+ntVars++*wordsize,"attenpos ");
00308   strcpy(ntTags+ntVars++*wordsize,"sfactpos ");
00309   strcpy(ntTags+ntVars++*wordsize,"npeneu   ");
00310   strcpy(ntTags+ntVars++*wordsize,"attenneu ");
00311   strcpy(ntTags+ntVars++*wordsize,"sfactneu ");
00312   strcpy(ntTags+ntVars++*wordsize,"npetl    ");
00313   strcpy(ntTags+ntVars++*wordsize,"attentl  ");
00314   strcpy(ntTags+ntVars++*wordsize,"sfacttl  ");
00315   strcpy(ntTags+ntVars++*wordsize,"npebi    ");
00316   strcpy(ntTags+ntVars++*wordsize,"attenbi  ");
00317   strcpy(ntTags+ntVars++*wordsize,"sfactbi  ");
00318   strcpy(ntTags+ntVars++*wordsize,"npecf    ");
00319   strcpy(ntTags+ntVars++*wordsize,"attencf  ");
00320   strcpy(ntTags+ntVars++*wordsize,"sfactcf  ");
00321   strcpy(ntTags+ntVars++*wordsize,"npemuo   ");
00322   strcpy(ntTags+ntVars++*wordsize,"attenmuo ");
00323   strcpy(ntTags+ntVars++*wordsize,"sfactmuo ");
00324   strcpy(ntTags+ntVars++*wordsize,"npeele   ");
00325   strcpy(ntTags+ntVars++*wordsize,"attenele ");
00326   strcpy(ntTags+ntVars++*wordsize,"sfactele ");
00327 
00328   strcpy(ntTags+ntVars++*wordsize,"ngampos  ");
00329   strcpy(ntTags+ntVars++*wordsize,"ngamneu  ");
00330   strcpy(ntTags+ntVars++*wordsize,"ngamtl   ");
00331   strcpy(ntTags+ntVars++*wordsize,"ngambi   ");
00332   strcpy(ntTags+ntVars++*wordsize,"ngamcf   ");
00333   strcpy(ntTags+ntVars++*wordsize,"ngammuo  ");
00334   strcpy(ntTags+ntVars++*wordsize,"ngamele  ");
00335   strcpy(ntTags+ntVars++*wordsize,"fpmtpos  ");
00336   strcpy(ntTags+ntVars++*wordsize,"fpmtneu  ");
00337   strcpy(ntTags+ntVars++*wordsize,"fpmttl   ");
00338   strcpy(ntTags+ntVars++*wordsize,"fpmtbi   ");
00339   strcpy(ntTags+ntVars++*wordsize,"fpmtcf   ");
00340   strcpy(ntTags+ntVars++*wordsize,"fpmtmuo  ");
00341   strcpy(ntTags+ntVars++*wordsize,"fpmtele  ");
00342 
00343   strcpy(ntTags+ntVars++*wordsize,"dposx    ");
00344   strcpy(ntTags+ntVars++*wordsize,"dposy    ");
00345   strcpy(ntTags+ntVars++*wordsize,"dposz    ");
00346   strcpy(ntTags+ntVars++*wordsize,"dneux    ");
00347   strcpy(ntTags+ntVars++*wordsize,"dneuy    ");
00348   strcpy(ntTags+ntVars++*wordsize,"dneuz    ");
00349   //
00350   IdDet = 1000*ntId+id++;
00351   sprintf(Name,"%s%s",ntName,", detector");
00352   NvarDet = ntVars;
00353   //
00354   hbookn_(IdDet,Name,NvarDet,dirName,1000,ntTags,
00355           strlen(Name),strlen(dirName),wordsize,NvarDet);
00356   //
00357   // define some histograms
00358   // Last argument is need for c++ ==> f strings
00359   //
00360   Offset = 10000*(ntId-1);
00361 
00362   hbook1_(Offset+101,"electron PH /PMT/gamma/electron",200,0,0.1,0,
00363           strlen("electron PH /PMT/gamma/electron"));
00364   hbook1_(Offset+102,"positron PH /PMT/gamma/positron",200,0,0.1,0,
00365           strlen("positron PH /PMT/gamma/positron"));
00366   hbook1_(Offset+103,"neutron PH /PMT/gamma/neutron",200,0,0.1,0,
00367           strlen("neutron PH /PMT/gamma/neutron"));
00368   hbook1_(Offset+104,"thallium PH /PMT/gamma/decay",200,0,0.1,0,
00369           strlen("thallium PH /PMT/gamma/decay"));
00370   hbook1_(Offset+105,"bismuth PH /PMT/gamma/decay",200,0,0.1,0,
00371           strlen("bismuth PH /PMT/gamma/decay"));
00372   hbook1_(Offset+106,"cf prompt gamma PH /PMT/gamma/decay",200,0,0.1,0,
00373           strlen("cf prompt gamma PH /PMT/gamma/decay"));
00374   hbook1_(Offset+107,"muon PH /PMT/gamma/muon",200,0,0.1,0,
00375           strlen("muon PH /PMT/gamma/muon"));
00376   //
00377   hbook1_(Offset+201,"electron time /PMT/gamma/electron",250,0,5e8,0,
00378           strlen("electron time /PMT/gamma/electron"));
00379   hbook1_(Offset+202,"positron time /PMT/gamma/positron",250,0,50,0,
00380           strlen("positron time /PMT/gamma/positron"));
00381   hbook1_(Offset+203,"neutron time /PMT/gamma/neutron",250,0,5e5,0,
00382           strlen("neutron time /PMT/gamma/neutron"));
00383   hbook1_(Offset+204,"thallium time /PMT/gamma/decay",250,0,1e6,0,
00384           strlen("thallium time /PMT/gamma/decay"));
00385   hbook1_(Offset+205,"bismuth time /PMT/gamma/decay",250,0,1e6,0,
00386           strlen("bismuth time /PMT/gamma/decay"));
00387   hbook1_(Offset+206,"cf prompt gamma time /PMT/gamma/decay",250,0,50,0,
00388           strlen("cf prompt gamma time /PMT/gamma/decay"));
00389   hbook1_(Offset+207,"muon time /PMT/gamma/muon",250,0,50,0,
00390           strlen("muon time /PMT/gamma/muon"));
00391   //
00392   hbook1_(Offset+301,"number of gamma/electron",10,0,10,0,
00393           strlen("number of gamma/electron"));
00394   hbook1_(Offset+302,"number of gamma/positron",10,0,10,0,
00395           strlen("number of gamma/positron"));
00396   hbook1_(Offset+303,"number of gamma/neutron",20,0,20,0,
00397           strlen("number of gamma/neutron"));
00398   hbook1_(Offset+306,"number of prompt Cf gamma/decay",20,0,20,0,
00399           strlen("number of prompt CF gamma/decay"));
00400   //
00401   hbook1_(Offset+401,"electron gamma energy",100,0,10,0,
00402           strlen("electron gamma energy"));
00403   hbook1_(Offset+402,"positron gamma energy",100,0,10,0,
00404           strlen("positron gamma energy"));
00405   hbook1_(Offset+403,"neutron gamma energy",100,0,10,0,
00406           strlen("neutron gamma energy"));
00407   hbook1_(Offset+406,"prompt Cf gamma energy",100,0,10,0,
00408           strlen("prompt Cf gamma energy"));
00409   //
00410   hbook1_(Offset+501,"electron gamma time",250,0,5e8,0,
00411           strlen("electron gamma time"));
00412   hbook1_(Offset+502,"positron gamma time",250,0,5,0,
00413           strlen("positron gamma time"));
00414   hbook1_(Offset+503,"neutron gamma time",250,0,5e5,0,
00415           strlen("neutron gamma time"));
00416   hbook1_(Offset+506,"prompt Cf gamma time",250,0,5,0,
00417           strlen("prompt Cf gamma time"));
00418   //
00419   hbook1_(Offset+601,"electron gamma radius",300,0,300,0,
00420           strlen("electron gamma radius"));
00421   hbook1_(Offset+602,"positron gamma radius",300,0,300,0,
00422           strlen("positron gamma radius"));
00423   hbook1_(Offset+603,"neutron gamma radius",300,0,300,0,
00424           strlen("neutron gamma radius"));
00425   hbook1_(Offset+606,"prompt Cf gamma radius",300,0,300,0,
00426           strlen("prompt Cf gamma radius"));
00427   //
00428   hbook1_(Offset+701,"electron gamma cos(theta)",20,-1.1,1.1,0,
00429           strlen("electron gamma cos(theta)"));
00430   hbook1_(Offset+702,"positron gamma cos(theta)",20,-1.1,1.1,0,
00431           strlen("positron gamma cos(theta)"));
00432   hbook1_(Offset+703,"neutron gamma cos(theta)",20,-1.1,1.1,0,
00433           strlen("neutron gamma cos(theta)"));
00434   hbook1_(Offset+706,"prompt Cf gamma cos(theta)",20,-1.1,1.1,0,
00435           strlen("prompt Cf gamma cos(theta)"));}
00436 
00437 ReactorNtuple::~ReactorNtuple(){
00438   ;
00439 }
00440 
00441 void ReactorNtuple::Fill(ReactorEvent& Event,ReactorDetector& Detector){
00442 
00443   float* tuple;
00444   int nmax = 0;
00445   int nloop = 0;
00446   //
00447   // fill the general event variables
00448   //
00449   //cout << "NvarGen " << NvarGen << endl;
00450   float TupleGen[NvarGen];
00451   tuple = TupleGen;
00452   //
00453   *tuple++ = Event.GetNevent();
00454   *tuple++ = Event.GetXCWeight();
00455   *tuple++ = Event.GetFluxWeight();
00456   //
00457   hfn_(IdGen,TupleGen);
00458   //
00459   // fill the nubar event variables
00460   //
00461   //cout << "NvarNu " << NvarNu << endl;
00462   float TupleNu[NvarNu];
00463   tuple = TupleNu;
00464   //
00465   // get number of nubar
00466   int nnu  = Event.GetNumNubar();
00467   //cout << nnu << " nubar(s)" << endl;
00468   *tuple++ = nnu;
00469   //
00470   // deal with multiple nubar
00471   double rnu[nnu];
00472   double cnu[nnu];
00473   double fnu[nnu];
00474   double xnu[nnu];
00475   double ynu[nnu];
00476   double znu[nnu];
00477   //
00478   // up to 1 nubar per event, keep the first one
00479   //
00480   nmax = 1;
00481   if(nnu < nmax)
00482     nloop = nnu;
00483   else
00484     nloop = nmax;
00485   //
00486   for(int n=0;n<nloop;n++){
00487     *tuple++ = Event.GetEnu(n);
00488     *tuple++ = Event.GetNubarParentProcess(n);
00489     //
00490     rnu[n]   = Event.GetNubarVertexR(n);
00491     cnu[n]   = Event.GetNubarVertexCosTheta(n);
00492     fnu[n]   = Event.GetNubarVertexPhi(n);
00493     xnu[n]   = rnu[n]*cos(fnu[n])*sqrt(1-cnu[n]*cnu[n]);
00494     ynu[n]   = rnu[n]*sin(fnu[n])*sqrt(1-cnu[n]*cnu[n]);
00495     znu[n]   = rnu[n]*cnu[n];
00496     //
00497     *tuple++ = rnu[n];
00498     *tuple++ = cnu[n];
00499     *tuple++ = fnu[n];
00500     *tuple++ = xnu[n];
00501     *tuple++ = ynu[n];
00502     *tuple++ = znu[n];
00503   }
00504   // zero out remaining tuples
00505   for(int n=nloop;n<nmax;n++){
00506     for(int i=0;i<((NvarNu-1)/nmax);i++) *tuple++ = 0;
00507   }
00508   //
00509   hfn_(IdNu,TupleNu);
00510   //
00511   // fill the electron event variables
00512   //
00513   //cout << "NvarEle " << NvarEle << endl;
00514   float TupleEle[NvarEle];
00515   tuple = TupleEle;
00516   //
00517   // get number of electron
00518   int nele = Event.GetNumElectron();
00519   //cout << nele << " electron(s)" << endl;
00520   *tuple++ = nele;
00521   //
00522   // deal with multiple electrons
00523   double rgenele[nele];
00524   double cgenele[nele];
00525   double fgenele[nele];
00526   double xgenele[nele];
00527   double ygenele[nele];
00528   double zgenele[nele];
00529   //
00530   double rele[nele];
00531   double cele[nele];
00532   double fele[nele];
00533   double xele[nele];
00534   double yele[nele];
00535   double zele[nele];
00536   //
00537   // up to 1 electron per event, keep the first one
00538   //
00539   nmax = 1;
00540   if(nele < nmax)
00541     nloop = nele;
00542   else
00543     nloop = nmax;
00544   //
00545   for(int n=0;n<nloop;n++){
00546     *tuple++ = Event.GetElectronKE(n);
00547     *tuple++ = Event.GetElectronCosTheta(n);
00548     *tuple++ = Event.GetElectronParentProcess(n);
00549     //
00550     rgenele[n]  = Event.GetGeneratedElectronVertexR(n);
00551     cgenele[n]  = Event.GetGeneratedElectronVertexCosTheta(n);
00552     fgenele[n]  = Event.GetGeneratedElectronVertexPhi(n);
00553     xgenele[n]  = rgenele[n]*cos(fgenele[n])*sqrt(1-cgenele[n]*cgenele[n]);
00554     ygenele[n]  = rgenele[n]*sin(fgenele[n])*sqrt(1-cgenele[n]*cgenele[n]);
00555     zgenele[n]  = rgenele[n]*cgenele[n];
00556     //
00557     *tuple++ = rgenele[n];
00558     *tuple++ = cgenele[n];
00559     *tuple++ = fgenele[n];
00560     *tuple++ = xgenele[n];
00561     *tuple++ = ygenele[n];
00562     *tuple++ = zgenele[n];
00563     //
00564     rele[n]  = Event.GetElectronVertexR(n);
00565     cele[n]  = Event.GetElectronVertexCosTheta(n);
00566     fele[n]  = Event.GetElectronVertexPhi(n);
00567     xele[n]  = rele[n]*cos(fele[n])*sqrt(1-cele[n]*cele[n]);
00568     yele[n]  = rele[n]*sin(fele[n])*sqrt(1-cele[n]*cele[n]);
00569     zele[n]  = rele[n]*cele[n];
00570     //
00571     *tuple++ = rele[n];
00572     *tuple++ = cele[n];
00573     *tuple++ = fele[n];
00574     *tuple++ = xele[n];
00575     *tuple++ = yele[n];
00576     *tuple++ = zele[n];
00577     //
00578     *tuple++ = Event.GetElectronTime(n);
00579   }
00580   // zero out remaining tuples
00581   for(int n=nloop;n<nmax;n++){
00582     for(int i=0;i<((NvarEle-1)/nmax);i++) *tuple++ = 0;
00583   }
00584   //
00585   hfn_(IdEle,TupleEle);
00586   //
00587   //
00588   // fill the positron event variables
00589   //
00590   //cout << "NvarPos " << NvarPos << endl;
00591   float TuplePos[NvarPos];
00592   tuple = TuplePos;
00593   //
00594   // get number of positron
00595   int npos = Event.GetNumPositron();
00596   //cout << npos << " positron(s)" << endl;
00597   *tuple++ = npos;
00598   //
00599   // deal with multiple positrons
00600   double rgenpos[npos];
00601   double cgenpos[npos];
00602   double fgenpos[npos];
00603   double xgenpos[npos];
00604   double ygenpos[npos];
00605   double zgenpos[npos];
00606   //
00607   double rpos[npos];
00608   double cpos[npos];
00609   double fpos[npos];
00610   double xpos[npos];
00611   double ypos[npos];
00612   double zpos[npos];
00613   //
00614   // up to 1 positron per event, keep the first one
00615   //
00616   nmax = 1;
00617   if(npos < nmax)
00618     nloop = npos;
00619   else
00620     nloop = nmax;
00621   //
00622   for(int n=0;n<nloop;n++){
00623     *tuple++ = Event.GetPositronKE(n);
00624     *tuple++ = Event.GetPositronCosTheta(n);
00625     *tuple++ = Event.GetPositronParentProcess(n);
00626     //
00627     rgenpos[n]  = Event.GetGeneratedPositronVertexR(n);
00628     cgenpos[n]  = Event.GetGeneratedPositronVertexCosTheta(n);
00629     fgenpos[n]  = Event.GetGeneratedPositronVertexPhi(n);
00630     xgenpos[n]  = rgenpos[n]*cos(fgenpos[n])*sqrt(1-cgenpos[n]*cgenpos[n]);
00631     ygenpos[n]  = rgenpos[n]*sin(fgenpos[n])*sqrt(1-cgenpos[n]*cgenpos[n]);
00632     zgenpos[n]  = rgenpos[n]*cgenpos[n];
00633     //
00634     *tuple++ = rgenpos[n];
00635     *tuple++ = cgenpos[n];
00636     *tuple++ = fgenpos[n];
00637     *tuple++ = xgenpos[n];
00638     *tuple++ = ygenpos[n];
00639     *tuple++ = zgenpos[n];
00640     //
00641     rpos[n]  = Event.GetPositronVertexR(n);
00642     cpos[n]  = Event.GetPositronVertexCosTheta(n);
00643     fpos[n]  = Event.GetPositronVertexPhi(n);
00644     xpos[n]  = rpos[n]*cos(fpos[n])*sqrt(1-cpos[n]*cpos[n]);
00645     ypos[n]  = rpos[n]*sin(fpos[n])*sqrt(1-cpos[n]*cpos[n]);
00646     zpos[n]  = rpos[n]*cpos[n];
00647     //
00648     *tuple++ = rpos[n];
00649     *tuple++ = cpos[n];
00650     *tuple++ = fpos[n];
00651     *tuple++ = xpos[n];
00652     *tuple++ = ypos[n];
00653     *tuple++ = zpos[n];
00654     //
00655     *tuple++ = Event.GetPositronTime(n);
00656   }
00657   // zero out remaining tuples
00658   for(int n=nloop;n<nmax;n++){
00659     for(int i=0;i<((NvarPos-1)/nmax);i++) *tuple++ = 0;
00660   }
00661   //
00662   hfn_(IdPos,TuplePos);
00663   //
00664   //
00665   // fill the neutron event variables
00666   //
00667   //cout << "NvarNeu " << NvarNeu << endl;
00668   float TupleNeu[NvarNeu];
00669   tuple = TupleNeu;
00670   //
00671   // get number of neutron
00672   int nneu = Event.GetNumNeutron();
00673   //cout << nneu << " neutron(s)" << endl;
00674   *tuple++ = nneu;
00675   //
00676   // deal with multiple neutrons
00677   double rgenneu[nneu];
00678   double cgenneu[nneu];
00679   double fgenneu[nneu];
00680   double xgenneu[nneu];
00681   double ygenneu[nneu];
00682   double zgenneu[nneu];
00683   //
00684   double rneu[nneu];
00685   double cneu[nneu];
00686   double fneu[nneu];
00687   double xneu[nneu];
00688   double yneu[nneu];
00689   double zneu[nneu];
00690   //
00691   // up to 5 neutrons per event, keep the first ones
00692   //
00693   nmax = 5;
00694   if(nneu < nmax)
00695     nloop = nneu;
00696   else
00697     nloop = nmax;
00698   //
00699   for(int n=0;n<nloop;n++){
00700     *tuple++ = Event.GetNeutronKE(n);
00701     *tuple++ = Event.GetNeutronCosTheta(n);
00702     *tuple++ = Event.GetNeutronParentProcess(n);
00703     //
00704     rgenneu[n]  = Event.GetGeneratedNeutronVertexR(n);
00705     cgenneu[n]  = Event.GetGeneratedNeutronVertexCosTheta(n);
00706     fgenneu[n]  = Event.GetGeneratedNeutronVertexPhi(n);
00707     xgenneu[n]  = rgenneu[n]*cos(fgenneu[n])*sqrt(1-cgenneu[n]*cgenneu[n]);
00708     ygenneu[n]  = rgenneu[n]*sin(fgenneu[n])*sqrt(1-cgenneu[n]*cgenneu[n]);
00709     zgenneu[n]  = rgenneu[n]*cgenneu[n];
00710     //
00711     *tuple++ = rgenneu[n];
00712     *tuple++ = cgenneu[n]; 
00713     *tuple++ = fgenneu[n];
00714     *tuple++ = xgenneu[n];
00715     *tuple++ = ygenneu[n];
00716     *tuple++ = zgenneu[n];
00717     //
00718     rneu[n]  = Event.GetNeutronVertexR(n);
00719     cneu[n]  = Event.GetNeutronVertexCosTheta(n);
00720     fneu[n]  = Event.GetNeutronVertexPhi(n);
00721     xneu[n]  = rneu[n]*cos(fneu[n])*sqrt(1-cneu[n]*cneu[n]);
00722     yneu[n]  = rneu[n]*sin(fneu[n])*sqrt(1-cneu[n]*cneu[n]);
00723     zneu[n]  = rneu[n]*cneu[n];
00724     //
00725     *tuple++ = rneu[n];
00726     *tuple++ = cneu[n]; 
00727     *tuple++ = fneu[n];
00728     *tuple++ = xneu[n];
00729     *tuple++ = yneu[n];
00730     *tuple++ = zneu[n];
00731     //
00732     *tuple++ = Event.GetNeutronTime(n);
00733     *tuple++ = Detector.GetNeutronSteps(n);
00734     *tuple++ = Detector.GetAbsorber(n);
00735   }
00736   // zero out remaining tuples
00737   for(int n=nloop;n<nmax;n++){
00738     for(int i=0;i<((NvarNeu-1)/nmax);i++) *tuple++ = 0;
00739   }
00740   //
00741   hfn_(IdNeu,TupleNeu);
00742   //
00743   //
00744   // fill the PMT detector variables
00745   //
00746   //cout << "NvarDet " << NvarDet << endl;
00747   float TupleDet[NvarDet];
00748   tuple = TupleDet;
00749   //
00750   *tuple++ = Detector.GetNcosbins();
00751   *tuple++ = Detector.GetNphibins();
00752   //
00753   *tuple++ = Detector.GetPHpositron();
00754   *tuple++ = Detector.GetPHposReco();
00755   *tuple++ = Detector.GetPHneutron();
00756   *tuple++ = Detector.GetPHneuReco();
00757   *tuple++ = Detector.GetPHthallium();
00758   *tuple++ = Detector.GetPHthaReco();
00759   *tuple++ = Detector.GetPHbismuth();
00760   *tuple++ = Detector.GetPHbisReco();
00761   *tuple++ = Detector.GetPHcalifornium();
00762   *tuple++ = Detector.GetPHcalReco();
00763   *tuple++ = Detector.GetPHmuon();
00764   *tuple++ = Detector.GetPHmuoReco();
00765   *tuple++ = Detector.GetPHelectron();
00766   *tuple++ = Detector.GetPHeleReco();
00767   //
00768   *tuple++ = Detector.GetMeanNpePos();
00769   *tuple++ = Detector.GetMeanAttenPos();
00770   *tuple++ = Detector.GetMeanSfactorPos();
00771   *tuple++ = Detector.GetMeanNpeNeu();
00772   *tuple++ = Detector.GetMeanAttenNeu();
00773   *tuple++ = Detector.GetMeanSfactorNeu();
00774   *tuple++ = Detector.GetMeanNpeTl();
00775   *tuple++ = Detector.GetMeanAttenTl();
00776   *tuple++ = Detector.GetMeanSfactorTl();
00777   *tuple++ = Detector.GetMeanNpeBi();
00778   *tuple++ = Detector.GetMeanAttenBi();
00779   *tuple++ = Detector.GetMeanSfactorBi();
00780   *tuple++ = Detector.GetMeanNpeCf();
00781   *tuple++ = Detector.GetMeanAttenCf();
00782   *tuple++ = Detector.GetMeanSfactorCf();
00783   *tuple++ = Detector.GetMeanNpeMuon();
00784   *tuple++ = Detector.GetMeanAttenMuon();
00785   *tuple++ = Detector.GetMeanSfactorMuon();
00786   *tuple++ = Detector.GetMeanNpeEle();
00787   *tuple++ = Detector.GetMeanAttenEle();
00788   *tuple++ = Detector.GetMeanSfactorEle();
00789 
00790   *tuple++ = Detector.GetNgamma("P"); 
00791   *tuple++ = Detector.GetNgamma("N");
00792   *tuple++ = Detector.GetNgamma("Tl");
00793   *tuple++ = Detector.GetNgamma("Bi");
00794   *tuple++ = Detector.GetNgamma("Cf");
00795   *tuple++ = Detector.GetNgamma("muon");
00796   *tuple++ = Detector.GetNgamma("E");
00797   *tuple++ = Detector.GetHitFracPos();  
00798   *tuple++ = Detector.GetHitFracNeu();  
00799   *tuple++ = Detector.GetHitFracTl();  
00800   *tuple++ = Detector.GetHitFracBi();  
00801   *tuple++ = Detector.GetHitFracCf();  
00802   *tuple++ = Detector.GetHitFracMuon();  
00803   *tuple++ = Detector.GetHitFracEle();  
00804   double dippos[3];
00805   double dipneu[3];
00806   for(int i=0;i<3;i++)dippos[i]=*(Detector.GetPositronDipole()+i);
00807   for(int i=0;i<3;i++)dipneu[i]=*(Detector.GetNeutronDipole()+i);
00808   //
00809   *tuple++ = dippos[0];  
00810   *tuple++ = dippos[1];  
00811   *tuple++ = dippos[2];  
00812   *tuple++ = dipneu[0];  
00813   *tuple++ = dipneu[1];  
00814   *tuple++ = dipneu[2];  
00815   //
00816   hfn_(IdDet,TupleDet);
00817   //
00818   //
00819   // do some looping over all the "PMT"
00820   //
00821   int npmt=Detector.GetNpmt();
00822   for(int n=0;n<npmt;n++){
00823 
00824     int nele;
00825     int npos;
00826     int nneu;
00827     int ntha;
00828     int nbis;
00829     int ncal;
00830     int nmuo;
00831     double xyzpmt[3];
00832     for(int i=0;i<3;i++) xyzpmt[i] = 0.;
00833     double phele[100];
00834     double phpos[100];
00835     double phneu[100];
00836     double phtha[100];
00837     double phbis[100];
00838     double phcal[100];
00839     double phmuo[100];
00840     double tele[100];
00841     double tpos[100];
00842     double tneu[100];
00843     double ttha[100];
00844     double tbis[100];
00845     double tcal[100];
00846     double tmuo[100];
00847     for(int i=0;i<100;i++){
00848       phele[i] = 0.;
00849       phpos[i] = 0.;
00850       phneu[i] = 0.;
00851       phtha[i] = 0.;
00852       phbis[i] = 0.;
00853       phcal[i] = 0.;
00854       phmuo[i] = 0.;
00855       tele[i] = 0.;
00856       tpos[i] = 0.;
00857       tneu[i] = 0.;
00858       ttha[i] = 0.;
00859       tbis[i] = 0.;
00860       tcal[i] = 0.;
00861       tmuo[i] = 0.;
00862     }
00863     //
00864     // Get "phototube" information.
00865     //      
00866     Detector.GetPMTdata(n,xyzpmt,nele,phele,tele,
00867                         npos,phpos,tpos,nneu,phneu,tneu,
00868                         ntha,phtha,ttha,nbis,phbis,tbis,
00869                         ncal,phcal,tcal,nmuo,phmuo,tmuo);
00870 
00871     for(int nn=0;nn<nele;nn++){
00872       hfill_(Offset+101,*(phele+nn));
00873       hfill_(Offset+201,*(tele+nn));
00874     }
00875     for(int nn=0;nn<npos;nn++){
00876       hfill_(Offset+102,*(phpos+nn));
00877       hfill_(Offset+202,*(tpos+nn));
00878     }
00879     for(int nn=0;nn<nneu;nn++){
00880       hfill_(Offset+103,*(phneu+nn));
00881       hfill_(Offset+203,*(tneu+nn));
00882     }
00883     for(int nn=0;nn<ntha;nn++){
00884       hfill_(Offset+104,*(phtha+nn));
00885       hfill_(Offset+204,*(ttha+nn));
00886     }
00887     for(int nn=0;nn<nbis;nn++){
00888       hfill_(Offset+105,*(phbis+nn));
00889       hfill_(Offset+205,*(tbis+nn));
00890     }
00891     for(int nn=0;nn<ncal;nn++){
00892       hfill_(Offset+106,*(phcal+nn));
00893       hfill_(Offset+206,*(tcal+nn));
00894     }
00895     for(int nn=0;nn<nmuo;nn++){
00896       hfill_(Offset+107,*(phmuo+nn));
00897       hfill_(Offset+207,*(tmuo+nn));
00898     }
00899   }
00900   //
00901   // Get information about gammas from electron.
00902   //
00903   int nelegamma = Detector.GetNgamma("E"); 
00904   hfill_(Offset+301,float(nelegamma));
00905   double ege[nelegamma];
00906   double tge[nelegamma];
00907   double rge[nelegamma];
00908   double cge[nelegamma];
00909   double fge[nelegamma];
00910   Detector.GetGamma("E",ege,tge,rge,cge,fge);
00911   for (int n=0;n<nelegamma;n++){
00912     hfill_(Offset+401,ege[n]);
00913     hfill_(Offset+501,tge[n]);
00914     hfill_(Offset+601,rge[n]);
00915     hfill_(Offset+701,cge[n]);
00916   }
00917   //
00918   // Get information about gammas from positron.
00919   // Note that the positron dedx is defined as a "gamma"
00920   // originating at the positron range point.
00921   //
00922   int nposgamma = Detector.GetNgamma("P"); 
00923   hfill_(Offset+302,float(nposgamma));
00924   double egp[nposgamma];
00925   double tgp[nposgamma];
00926   double rgp[nposgamma];
00927   double cgp[nposgamma];
00928   double fgp[nposgamma];
00929   Detector.GetGamma("P",egp,tgp,rgp,cgp,fgp);
00930   for (int n=0;n<nposgamma;n++){
00931     hfill_(Offset+402,egp[n]);
00932     hfill_(Offset+502,tgp[n]);
00933     hfill_(Offset+602,rgp[n]);
00934     hfill_(Offset+702,cgp[n]);
00935   }
00936   //
00937   // Get information about gammas from neutron.
00938   //
00939   int nneugamma = Detector.GetNgamma("N"); 
00940   hfill_(Offset+303,float(nneugamma));
00941   double egn[nneugamma];
00942   double tgn[nneugamma];
00943   double rgn[nneugamma];
00944   double cgn[nneugamma];
00945   double fgn[nneugamma];
00946   Detector.GetGamma("N",egn,tgn,rgn,cgn,fgn);
00947   for (int n=0;n<nneugamma;n++){
00948     hfill_(Offset+403,egn[n]);
00949     hfill_(Offset+503,tgn[n]);
00950     hfill_(Offset+603,rgn[n]);
00951     hfill_(Offset+703,cgn[n]);
00952   }
00953   //
00954   // Get information about prompt gammas from 252cf.
00955   //
00956   int ngamma = Detector.GetNgamma("Cf"); 
00957   hfill_(Offset+306,float(ngamma));
00958   double egcf[ngamma];
00959   double tgcf[ngamma];
00960   double rgcf[ngamma];
00961   double cgcf[ngamma];
00962   double fgcf[ngamma];
00963   Detector.GetGamma("Cf",egcf,tgcf,rgcf,cgcf,fgcf);
00964   for (int n=0;n<ngamma;n++){
00965     hfill_(Offset+406,egcf[n]);
00966     hfill_(Offset+506,tgcf[n]);
00967     hfill_(Offset+606,rgcf[n]);
00968     hfill_(Offset+706,cgcf[n]);
00969   }
00970 }
00971