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 
00018 ReactorNtuple::ReactorNtuple(const ReactorNtuple& Ntuple){
00019 
00020   IdGen = Ntuple.IdGen;
00021   IdNu  = Ntuple.IdNu;
00022   IdPos = Ntuple.IdPos;
00023   IdNeu = Ntuple.IdNeu;
00024   IdDet = Ntuple.IdDet;
00025 
00026   NvarGen = Ntuple.NvarGen;
00027   NvarNu  = Ntuple.NvarNu;
00028   NvarPos = Ntuple.NvarPos;
00029   NvarNeu = Ntuple.NvarNeu;
00030   NvarDet = Ntuple.NvarDet;
00031 
00032   Offset = Ntuple.Offset;
00033 }
00034 ReactorNtuple& ReactorNtuple::operator=(const ReactorNtuple& rhs){
00035 
00036   if(this == &rhs)return *this;
00037   IdGen = rhs.IdGen;
00038   IdNu  = rhs.IdNu;
00039   IdPos = rhs.IdPos;
00040   IdNeu = rhs.IdNeu;
00041   IdDet = rhs.IdDet;
00042 
00043   NvarGen = rhs.NvarGen;
00044   NvarNu  = rhs.NvarNu;
00045   NvarPos = rhs.NvarPos;
00046   NvarNeu = rhs.NvarNeu;
00047   NvarDet = rhs.NvarDet;
00048 
00049   Offset = rhs.Offset;
00050   return *this;
00051 }
00052 
00053 ReactorNtuple::ReactorNtuple(int ntId,
00054                              const char* ntName,
00055                              char* dirName){
00056 
00057   // define ntuple tags here
00058   // to add vars, increment ntVars, and add tag to end of list
00059   // don't change wordsize
00060   //
00061   int wordsize = 8;
00062   int ntVarsMax = 100;
00063   char ntTags[ntVarsMax*wordsize];
00064   //
00065   int ntVars = 0;
00066   int id = 1;
00067   char Name[100];
00068   //
00069   // book the general event veriables, nt=1001
00070   //
00071   strcpy(ntTags+ntVars++*wordsize,"nevents  ");
00072   strcpy(ntTags+ntVars++*wordsize,"wtc      ");
00073   strcpy(ntTags+ntVars++*wordsize,"wtf      ");
00074   //
00075   IdGen = 1000*ntId+id++;
00076   sprintf(Name,"%s%s",ntName,", general");
00077   NvarGen = ntVars;
00078   //
00079   hbookn_(IdGen,Name,NvarGen,dirName,1000,ntTags,
00080           strlen(Name),strlen(dirName),wordsize,NvarGen);
00081   //
00082   // reset
00083   ntVars = 0;
00084   //
00085   // book the nubar event variables, nt=1002
00086   //
00087   // up to 5 nubar per event
00088   strcpy(ntTags+ntVars++*wordsize,"nnu      ");
00089 
00090   strcpy(ntTags+ntVars++*wordsize,"enu1     ");
00091   strcpy(ntTags+ntVars++*wordsize,"nuproc1  ");
00092   strcpy(ntTags+ntVars++*wordsize,"vnu1r    ");
00093   strcpy(ntTags+ntVars++*wordsize,"vnu1cos  ");
00094   strcpy(ntTags+ntVars++*wordsize,"vnu1phi  ");
00095   strcpy(ntTags+ntVars++*wordsize,"vnu1x    ");
00096   strcpy(ntTags+ntVars++*wordsize,"vnu1y    ");
00097   strcpy(ntTags+ntVars++*wordsize,"vnu1z    ");
00098 
00099   strcpy(ntTags+ntVars++*wordsize,"enu2     ");
00100   strcpy(ntTags+ntVars++*wordsize,"nuproc2  ");
00101   strcpy(ntTags+ntVars++*wordsize,"vnu2r    ");
00102   strcpy(ntTags+ntVars++*wordsize,"vnu2cos  ");
00103   strcpy(ntTags+ntVars++*wordsize,"vnu2phi  ");
00104   strcpy(ntTags+ntVars++*wordsize,"vnu2x    ");
00105   strcpy(ntTags+ntVars++*wordsize,"vnu2y    ");
00106   strcpy(ntTags+ntVars++*wordsize,"vnu2z    ");
00107 
00108   strcpy(ntTags+ntVars++*wordsize,"enu3     ");
00109   strcpy(ntTags+ntVars++*wordsize,"nuproc3  ");
00110   strcpy(ntTags+ntVars++*wordsize,"vnu3r    ");
00111   strcpy(ntTags+ntVars++*wordsize,"vnu3cos  ");
00112   strcpy(ntTags+ntVars++*wordsize,"vnu3phi  ");
00113   strcpy(ntTags+ntVars++*wordsize,"vnu3x    ");
00114   strcpy(ntTags+ntVars++*wordsize,"vnu3y    ");
00115   strcpy(ntTags+ntVars++*wordsize,"vnu3z    ");
00116 
00117   strcpy(ntTags+ntVars++*wordsize,"enu4     ");
00118   strcpy(ntTags+ntVars++*wordsize,"nuproc4  ");
00119   strcpy(ntTags+ntVars++*wordsize,"vnu4r    ");
00120   strcpy(ntTags+ntVars++*wordsize,"vnu4cos  ");
00121   strcpy(ntTags+ntVars++*wordsize,"vnu4phi  ");
00122   strcpy(ntTags+ntVars++*wordsize,"vnu4x    ");
00123   strcpy(ntTags+ntVars++*wordsize,"vnu4y    ");
00124   strcpy(ntTags+ntVars++*wordsize,"vnu4z    ");
00125 
00126   strcpy(ntTags+ntVars++*wordsize,"enu5     ");
00127   strcpy(ntTags+ntVars++*wordsize,"nuproc5  ");
00128   strcpy(ntTags+ntVars++*wordsize,"vnu5r    ");
00129   strcpy(ntTags+ntVars++*wordsize,"vnu5cos  ");
00130   strcpy(ntTags+ntVars++*wordsize,"vnu5phi  ");
00131   strcpy(ntTags+ntVars++*wordsize,"vnu5x    ");
00132   strcpy(ntTags+ntVars++*wordsize,"vnu5y    ");
00133   strcpy(ntTags+ntVars++*wordsize,"vnu5z    ");
00134   //
00135   IdNu = 1000*ntId+id++;
00136   sprintf(Name,"%s%s",ntName,", nubar");
00137   NvarNu = ntVars;
00138   //
00139   hbookn_(IdNu,Name,NvarNu,dirName,1000,ntTags,
00140           strlen(Name),strlen(dirName),wordsize,NvarNu);
00141   //
00142   // reset
00143   ntVars = 0;
00144   //
00145   // book the positron event variables, nt=1003
00146   //
00147   // up to 5 positron per event
00148   strcpy(ntTags+ntVars++*wordsize,"npos     ");
00149 
00150   strcpy(ntTags+ntVars++*wordsize,"KEpos1   ");
00151   strcpy(ntTags+ntVars++*wordsize,"pos1cos  ");
00152   strcpy(ntTags+ntVars++*wordsize,"posproc1 ");
00153   strcpy(ntTags+ntVars++*wordsize,"gpos1r   ");
00154   strcpy(ntTags+ntVars++*wordsize,"gpos1cos ");
00155   strcpy(ntTags+ntVars++*wordsize,"gpos1phi ");
00156   strcpy(ntTags+ntVars++*wordsize,"gpos1x   ");
00157   strcpy(ntTags+ntVars++*wordsize,"gpos1y   ");
00158   strcpy(ntTags+ntVars++*wordsize,"gpos1z   ");
00159   strcpy(ntTags+ntVars++*wordsize,"vpos1r   ");
00160   strcpy(ntTags+ntVars++*wordsize,"vpos1cos ");
00161   strcpy(ntTags+ntVars++*wordsize,"vpos1phi ");
00162   strcpy(ntTags+ntVars++*wordsize,"vpos1x   ");
00163   strcpy(ntTags+ntVars++*wordsize,"vpos1y   ");
00164   strcpy(ntTags+ntVars++*wordsize,"vpos1z   ");
00165   strcpy(ntTags+ntVars++*wordsize,"tpos1    ");
00166 
00167   strcpy(ntTags+ntVars++*wordsize,"KEpos2   ");
00168   strcpy(ntTags+ntVars++*wordsize,"pos2cos  ");
00169   strcpy(ntTags+ntVars++*wordsize,"posproc2 ");
00170   strcpy(ntTags+ntVars++*wordsize,"gpos2r   ");
00171   strcpy(ntTags+ntVars++*wordsize,"gpos2cos ");
00172   strcpy(ntTags+ntVars++*wordsize,"gpos2phi ");
00173   strcpy(ntTags+ntVars++*wordsize,"gpos2x   ");
00174   strcpy(ntTags+ntVars++*wordsize,"gpos2y   ");
00175   strcpy(ntTags+ntVars++*wordsize,"gpos2z   ");
00176   strcpy(ntTags+ntVars++*wordsize,"vpos2r   ");
00177   strcpy(ntTags+ntVars++*wordsize,"vpos2cos ");
00178   strcpy(ntTags+ntVars++*wordsize,"vpos2phi ");
00179   strcpy(ntTags+ntVars++*wordsize,"vpos2x   ");
00180   strcpy(ntTags+ntVars++*wordsize,"vpos2y   ");
00181   strcpy(ntTags+ntVars++*wordsize,"vpos2z   ");
00182   strcpy(ntTags+ntVars++*wordsize,"tpos2    ");
00183 
00184   strcpy(ntTags+ntVars++*wordsize,"KEpos3   ");
00185   strcpy(ntTags+ntVars++*wordsize,"pos3cos  ");
00186   strcpy(ntTags+ntVars++*wordsize,"posproc3 ");
00187   strcpy(ntTags+ntVars++*wordsize,"gpos3r   ");
00188   strcpy(ntTags+ntVars++*wordsize,"gpos3cos ");
00189   strcpy(ntTags+ntVars++*wordsize,"gpos3phi ");
00190   strcpy(ntTags+ntVars++*wordsize,"gpos3x   ");
00191   strcpy(ntTags+ntVars++*wordsize,"gpos3y   ");
00192   strcpy(ntTags+ntVars++*wordsize,"gpos3z   ");
00193   strcpy(ntTags+ntVars++*wordsize,"vpos3r   ");
00194   strcpy(ntTags+ntVars++*wordsize,"vpos3cos ");
00195   strcpy(ntTags+ntVars++*wordsize,"vpos3phi ");
00196   strcpy(ntTags+ntVars++*wordsize,"vpos3x   ");
00197   strcpy(ntTags+ntVars++*wordsize,"vpos3y   ");
00198   strcpy(ntTags+ntVars++*wordsize,"vpos3z   ");
00199   strcpy(ntTags+ntVars++*wordsize,"tpos3    ");
00200 
00201   strcpy(ntTags+ntVars++*wordsize,"KEpos4   ");
00202   strcpy(ntTags+ntVars++*wordsize,"pos4cos  ");
00203   strcpy(ntTags+ntVars++*wordsize,"posproc4 ");
00204   strcpy(ntTags+ntVars++*wordsize,"gpos4r   ");
00205   strcpy(ntTags+ntVars++*wordsize,"gpos4cos ");
00206   strcpy(ntTags+ntVars++*wordsize,"gpos4phi ");
00207   strcpy(ntTags+ntVars++*wordsize,"gpos4x   ");
00208   strcpy(ntTags+ntVars++*wordsize,"gpos4y   ");
00209   strcpy(ntTags+ntVars++*wordsize,"gpos4z   ");
00210   strcpy(ntTags+ntVars++*wordsize,"vpos4r   ");
00211   strcpy(ntTags+ntVars++*wordsize,"vpos4cos ");
00212   strcpy(ntTags+ntVars++*wordsize,"vpos4phi ");
00213   strcpy(ntTags+ntVars++*wordsize,"vpos4x   ");
00214   strcpy(ntTags+ntVars++*wordsize,"vpos4y   ");
00215   strcpy(ntTags+ntVars++*wordsize,"vpos4z   ");
00216   strcpy(ntTags+ntVars++*wordsize,"tpos4    ");
00217 
00218   strcpy(ntTags+ntVars++*wordsize,"KEpos5   ");
00219   strcpy(ntTags+ntVars++*wordsize,"pos5cos  ");
00220   strcpy(ntTags+ntVars++*wordsize,"posproc5 ");
00221   strcpy(ntTags+ntVars++*wordsize,"gpos5r   ");
00222   strcpy(ntTags+ntVars++*wordsize,"gpos5cos ");
00223   strcpy(ntTags+ntVars++*wordsize,"gpos5phi ");
00224   strcpy(ntTags+ntVars++*wordsize,"gpos5x   ");
00225   strcpy(ntTags+ntVars++*wordsize,"gpos5y   ");
00226   strcpy(ntTags+ntVars++*wordsize,"gpos5z   ");
00227   strcpy(ntTags+ntVars++*wordsize,"vpos5r   ");
00228   strcpy(ntTags+ntVars++*wordsize,"vpos5cos ");
00229   strcpy(ntTags+ntVars++*wordsize,"vpos5phi ");
00230   strcpy(ntTags+ntVars++*wordsize,"vpos5x   ");
00231   strcpy(ntTags+ntVars++*wordsize,"vpos5y   ");
00232   strcpy(ntTags+ntVars++*wordsize,"vpos5z   ");
00233   strcpy(ntTags+ntVars++*wordsize,"tpos5    ");
00234   //
00235   IdPos = 1000*ntId+id++;
00236   sprintf(Name,"%s%s",ntName,", positron");
00237   NvarPos = ntVars;
00238   //
00239   hbookn_(IdPos,Name,NvarPos,dirName,1000,ntTags,
00240           strlen(Name),strlen(dirName),wordsize,NvarPos);
00241   //
00242   // reset
00243   ntVars = 0;
00244   //
00245   // book the neutron event variables, nt=1004
00246   //
00247   // up to 5 neutron per event
00248   strcpy(ntTags+ntVars++*wordsize,"nneu     ");
00249 
00250   strcpy(ntTags+ntVars++*wordsize,"KEneu1   ");
00251   strcpy(ntTags+ntVars++*wordsize,"neu1cos  ");
00252   strcpy(ntTags+ntVars++*wordsize,"neuproc1 ");
00253   strcpy(ntTags+ntVars++*wordsize,"gneu1r   ");
00254   strcpy(ntTags+ntVars++*wordsize,"gneu1cos ");
00255   strcpy(ntTags+ntVars++*wordsize,"gneu1phi ");
00256   strcpy(ntTags+ntVars++*wordsize,"gneu1x   ");
00257   strcpy(ntTags+ntVars++*wordsize,"gneu1y   ");
00258   strcpy(ntTags+ntVars++*wordsize,"gneu1z   ");
00259   strcpy(ntTags+ntVars++*wordsize,"vneu1r   ");
00260   strcpy(ntTags+ntVars++*wordsize,"vneu1cos ");
00261   strcpy(ntTags+ntVars++*wordsize,"vneu1phi ");
00262   strcpy(ntTags+ntVars++*wordsize,"vneu1x   ");
00263   strcpy(ntTags+ntVars++*wordsize,"vneu1y   ");
00264   strcpy(ntTags+ntVars++*wordsize,"vneu1z   ");
00265   strcpy(ntTags+ntVars++*wordsize,"tneu1    ");
00266   strcpy(ntTags+ntVars++*wordsize,"nsteps1  ");
00267   strcpy(ntTags+ntVars++*wordsize,"absorb1  ");
00268 
00269   strcpy(ntTags+ntVars++*wordsize,"KEneu2   ");
00270   strcpy(ntTags+ntVars++*wordsize,"neu2cos  ");
00271   strcpy(ntTags+ntVars++*wordsize,"neuproc2 ");
00272   strcpy(ntTags+ntVars++*wordsize,"gneu2r   ");
00273   strcpy(ntTags+ntVars++*wordsize,"gneu2cos ");
00274   strcpy(ntTags+ntVars++*wordsize,"gneu2phi ");
00275   strcpy(ntTags+ntVars++*wordsize,"gneu2x   ");
00276   strcpy(ntTags+ntVars++*wordsize,"gneu2y   ");
00277   strcpy(ntTags+ntVars++*wordsize,"gneu2z   ");
00278   strcpy(ntTags+ntVars++*wordsize,"vneu2r   ");
00279   strcpy(ntTags+ntVars++*wordsize,"vneu2cos ");
00280   strcpy(ntTags+ntVars++*wordsize,"vneu2phi ");
00281   strcpy(ntTags+ntVars++*wordsize,"vneu2x   ");
00282   strcpy(ntTags+ntVars++*wordsize,"vneu2y   ");
00283   strcpy(ntTags+ntVars++*wordsize,"vneu2z   ");
00284   strcpy(ntTags+ntVars++*wordsize,"tneu2    ");
00285   strcpy(ntTags+ntVars++*wordsize,"nsteps2  ");
00286   strcpy(ntTags+ntVars++*wordsize,"absorb2  ");
00287 
00288   strcpy(ntTags+ntVars++*wordsize,"KEneu3   ");
00289   strcpy(ntTags+ntVars++*wordsize,"neu3cos  ");
00290   strcpy(ntTags+ntVars++*wordsize,"neuproc3 ");
00291   strcpy(ntTags+ntVars++*wordsize,"gneu3r   ");
00292   strcpy(ntTags+ntVars++*wordsize,"gneu3cos ");
00293   strcpy(ntTags+ntVars++*wordsize,"gneu3phi ");
00294   strcpy(ntTags+ntVars++*wordsize,"gneu3x   ");
00295   strcpy(ntTags+ntVars++*wordsize,"gneu3y   ");
00296   strcpy(ntTags+ntVars++*wordsize,"gneu3z   ");
00297   strcpy(ntTags+ntVars++*wordsize,"vneu3r   ");
00298   strcpy(ntTags+ntVars++*wordsize,"vneu3cos ");
00299   strcpy(ntTags+ntVars++*wordsize,"vneu3phi ");
00300   strcpy(ntTags+ntVars++*wordsize,"vneu3x   ");
00301   strcpy(ntTags+ntVars++*wordsize,"vneu3y   ");
00302   strcpy(ntTags+ntVars++*wordsize,"vneu3z   ");
00303   strcpy(ntTags+ntVars++*wordsize,"tneu3    ");
00304   strcpy(ntTags+ntVars++*wordsize,"nsteps3  ");
00305   strcpy(ntTags+ntVars++*wordsize,"absorb3  ");
00306 
00307   strcpy(ntTags+ntVars++*wordsize,"KEneu4   ");
00308   strcpy(ntTags+ntVars++*wordsize,"neu4cos  ");
00309   strcpy(ntTags+ntVars++*wordsize,"neuproc4 ");
00310   strcpy(ntTags+ntVars++*wordsize,"gneu4r   ");
00311   strcpy(ntTags+ntVars++*wordsize,"gneu4cos ");
00312   strcpy(ntTags+ntVars++*wordsize,"gneu4phi ");
00313   strcpy(ntTags+ntVars++*wordsize,"gneu4x   ");
00314   strcpy(ntTags+ntVars++*wordsize,"gneu4y   ");
00315   strcpy(ntTags+ntVars++*wordsize,"gneu4z   ");
00316   strcpy(ntTags+ntVars++*wordsize,"vneu4r   ");
00317   strcpy(ntTags+ntVars++*wordsize,"vneu4cos ");
00318   strcpy(ntTags+ntVars++*wordsize,"vneu4phi ");
00319   strcpy(ntTags+ntVars++*wordsize,"vneu4x   ");
00320   strcpy(ntTags+ntVars++*wordsize,"vneu4y   ");
00321   strcpy(ntTags+ntVars++*wordsize,"vneu4z   ");
00322   strcpy(ntTags+ntVars++*wordsize,"tneu4    ");
00323   strcpy(ntTags+ntVars++*wordsize,"nsteps4  ");
00324   strcpy(ntTags+ntVars++*wordsize,"absorb4  ");
00325 
00326   strcpy(ntTags+ntVars++*wordsize,"KEneu5   ");
00327   strcpy(ntTags+ntVars++*wordsize,"neu5cos  ");
00328   strcpy(ntTags+ntVars++*wordsize,"neuproc5 ");
00329   strcpy(ntTags+ntVars++*wordsize,"gneu5r   ");
00330   strcpy(ntTags+ntVars++*wordsize,"gneu5cos ");
00331   strcpy(ntTags+ntVars++*wordsize,"gneu5phi ");
00332   strcpy(ntTags+ntVars++*wordsize,"gneu5x   ");
00333   strcpy(ntTags+ntVars++*wordsize,"gneu5y   ");
00334   strcpy(ntTags+ntVars++*wordsize,"gneu5z   ");
00335   strcpy(ntTags+ntVars++*wordsize,"vneu5r   ");
00336   strcpy(ntTags+ntVars++*wordsize,"vneu5cos ");
00337   strcpy(ntTags+ntVars++*wordsize,"vneu5phi ");
00338   strcpy(ntTags+ntVars++*wordsize,"vneu5x   ");
00339   strcpy(ntTags+ntVars++*wordsize,"vneu5y   ");
00340   strcpy(ntTags+ntVars++*wordsize,"vneu5z   ");
00341   strcpy(ntTags+ntVars++*wordsize,"tneu5    ");
00342   strcpy(ntTags+ntVars++*wordsize,"nsteps5  ");
00343   strcpy(ntTags+ntVars++*wordsize,"absorb5  ");
00344   //
00345   IdNeu = 1000*ntId+id++;
00346   sprintf(Name,"%s%s",ntName,", neutron");
00347   NvarNeu = ntVars;
00348   //
00349   hbookn_(IdNeu,Name,NvarNeu,dirName,1000,ntTags,
00350           strlen(Name),strlen(dirName),wordsize,NvarNeu);
00351   //
00352   // reset
00353   ntVars = 0;
00354   //
00355   // book the pmt detector variables, nt=1005
00356   //
00357   strcpy(ntTags+ntVars++*wordsize,"ncosbins ");
00358   strcpy(ntTags+ntVars++*wordsize,"nphibins ");
00359 
00360   strcpy(ntTags+ntVars++*wordsize,"phposvis ");
00361   strcpy(ntTags+ntVars++*wordsize,"phposrec ");
00362   strcpy(ntTags+ntVars++*wordsize,"phneuvis ");
00363   strcpy(ntTags+ntVars++*wordsize,"phneurec ");
00364 
00365   strcpy(ntTags+ntVars++*wordsize,"npepos   ");
00366   strcpy(ntTags+ntVars++*wordsize,"attenpos ");
00367   strcpy(ntTags+ntVars++*wordsize,"sfactpos ");
00368   strcpy(ntTags+ntVars++*wordsize,"npeneu   ");
00369   strcpy(ntTags+ntVars++*wordsize,"attenneu ");
00370   strcpy(ntTags+ntVars++*wordsize,"sfactneu ");
00371   strcpy(ntTags+ntVars++*wordsize,"ngampos  ");
00372   strcpy(ntTags+ntVars++*wordsize,"ngamneu  ");
00373   strcpy(ntTags+ntVars++*wordsize,"ngamtha  ");
00374   strcpy(ntTags+ntVars++*wordsize,"fpmtpos  ");
00375   strcpy(ntTags+ntVars++*wordsize,"fpmtneu  ");
00376 
00377   strcpy(ntTags+ntVars++*wordsize,"dposx    ");
00378   strcpy(ntTags+ntVars++*wordsize,"dposy    ");
00379   strcpy(ntTags+ntVars++*wordsize,"dposz    ");
00380   strcpy(ntTags+ntVars++*wordsize,"dneux    ");
00381   strcpy(ntTags+ntVars++*wordsize,"dneuy    ");
00382   strcpy(ntTags+ntVars++*wordsize,"dneuz    ");
00383   //
00384   IdDet = 1000*ntId+id++;
00385   sprintf(Name,"%s%s",ntName,", detector");
00386   NvarDet = ntVars;
00387   //
00388   hbookn_(IdDet,Name,NvarDet,dirName,1000,ntTags,
00389           strlen(Name),strlen(dirName),wordsize,NvarDet);
00390   //
00391   // define some histograms
00392   // Last argument is need for c++ ==> f strings
00393   //
00394   Offset = 10000*(ntId-1);
00395   /*
00396   hbook1_(Offset+101,"phpos",200,0,0.1,0,strlen("phpos"));
00397   hbook1_(Offset+102,"phneu",200,0,0.1,0,strlen("phneu"));
00398   //
00399   hbook1_(Offset+201,"tpos",250,0,25,0,strlen("tpos"));
00400   hbook1_(Offset+202,"tneu",250,0,500000,0,strlen("tneu"));
00401   //
00402   hbook1_(Offset+301,"ngpos",10,0,10,0,strlen("ngpos"));
00403   hbook1_(Offset+302,"ngneu",20,0,20,0,strlen("ngneu"));
00404   //
00405   hbook1_(Offset+401,"egpos",200,0,10,0,strlen("egpos"));
00406   hbook1_(Offset+402,"egneu",200,0,20,0,strlen("egneu"));
00407   //
00408   hbook1_(Offset+501,"rgp-rpos",200,0,100,0,strlen("rgp-rpos"));
00409   hbook1_(Offset+502,"rgn-rneu",200,0,100,0,strlen("rgn-rneu"));
00410   //
00411   hbook1_(Offset+601,"rgp-rnu",200,0,100,0,strlen("rgp-rnu"));
00412   hbook1_(Offset+602,"rgn-rnu",200,0,100,0,strlen("rgn-rnu"));
00413   //
00414   hbook1_(Offset+701,"tgp-tpos",250,0,5,0,strlen("tgp-tpos"));
00415   hbook1_(Offset+702,"tgn-tneu",250,0,5,0,strlen("tgn-tneu"));
00416   //
00417   hbook1_(Offset+801,"tgp-tnu",250,0,25,0,strlen("tgp-tnu"));
00418   hbook1_(Offset+802,"tgn-tnu",250,0,1000000,0,strlen("tgn-tnu"));
00419   */
00420 }
00421 
00422 ReactorNtuple::~ReactorNtuple(){
00423   ;
00424 }
00425 
00426 void ReactorNtuple::Fill(ReactorEvent& Event,ReactorDetector& Detector){
00427 
00428   float* tuple;
00429   //
00430   // fill the general event variables
00431   //
00432   //cout << "NvarGen " << NvarGen << endl;
00433   float TupleGen[NvarGen];
00434   tuple = TupleGen;
00435   //
00436   *tuple++ = Event.GetNevents();
00437   *tuple++ = Event.GetXCWeight();
00438   *tuple++ = Event.GetFluxWeight();
00439   //
00440   hfn_(IdGen,TupleGen);
00441   //
00442   // fill the nubar event variables
00443   //
00444   //cout << "NvarNu " << NvarNu << endl;
00445   float TupleNu[NvarNu];
00446   tuple = TupleNu;
00447   //
00448   // get number of nubar
00449   int nnu  = Event.GetNumNubar();
00450   //cout << nnu << " nubar(s)" << endl;
00451   *tuple++ = nnu;
00452   //
00453   // deal with multiple nubar
00454   double rnu[nnu];
00455   double cnu[nnu];
00456   double fnu[nnu];
00457   double xnu[nnu];
00458   double ynu[nnu];
00459   double znu[nnu];
00460   //
00461   for(int n=0;n<nnu;n++){
00462     *tuple++ = Event.GetEnu(n);
00463     *tuple++ = Event.GetNubarParentProcess(n);
00464     //
00465     rnu[n]   = Event.GetNubarVertexR(n);
00466     cnu[n]   = Event.GetNubarVertexCosTheta(n);
00467     fnu[n]   = Event.GetNubarVertexPhi(n);
00468     xnu[n]   = rnu[n]*cos(fnu[n])*sqrt(1-cnu[n]*cnu[n]);
00469     ynu[n]   = rnu[n]*sin(fnu[n])*sqrt(1-cnu[n]*cnu[n]);
00470     znu[n]   = rnu[n]*cnu[n];
00471     //
00472     *tuple++ = rnu[n];
00473     *tuple++ = cnu[n];
00474     *tuple++ = fnu[n];
00475     *tuple++ = xnu[n];
00476     *tuple++ = ynu[n];
00477     *tuple++ = znu[n];
00478   }
00479   // zero out remaining tuples
00480   for(int n=nnu;n<5;n++){
00481     for(int i=0;i<((NvarNu-1)/5);i++) *tuple++ = 0;
00482   }
00483   //
00484   hfn_(IdNu,TupleNu);
00485   //
00486   // fill the positron event variables
00487   //
00488   //cout << "NvarPos " << NvarPos << endl;
00489   float TuplePos[NvarPos];
00490   tuple = TuplePos;
00491   //
00492   // get number of positron
00493   int npos = Event.GetNumPositron();
00494   //cout << npos << " positron(s)" << endl;
00495   *tuple++ = npos;
00496   //
00497   // deal with multiple positrons
00498   double rgenpos[npos];
00499   double cgenpos[npos];
00500   double fgenpos[npos];
00501   double xgenpos[npos];
00502   double ygenpos[npos];
00503   double zgenpos[npos];
00504   //
00505   double rpos[npos];
00506   double cpos[npos];
00507   double fpos[npos];
00508   double xpos[npos];
00509   double ypos[npos];
00510   double zpos[npos];
00511   //
00512   for(int n=0;n<npos;n++){
00513     *tuple++ = Event.GetPositronKE(n);
00514     *tuple++ = Event.GetPositronCosTheta(n);
00515     *tuple++ = Event.GetPositronParentProcess(n);
00516     //
00517     rgenpos[n]  = Event.GetGeneratedPositronVertexR(n);
00518     cgenpos[n]  = Event.GetGeneratedPositronVertexCosTheta(n);
00519     fgenpos[n]  = Event.GetGeneratedPositronVertexPhi(n);
00520     xgenpos[n]  = rgenpos[n]*cos(fgenpos[n])*sqrt(1-cgenpos[n]*cgenpos[n]);
00521     ygenpos[n]  = rgenpos[n]*sin(fgenpos[n])*sqrt(1-cgenpos[n]*cgenpos[n]);
00522     zgenpos[n]  = rgenpos[n]*cgenpos[n];
00523     //
00524     *tuple++ = rgenpos[n];
00525     *tuple++ = cgenpos[n];
00526     *tuple++ = fgenpos[n];
00527     *tuple++ = xgenpos[n];
00528     *tuple++ = ygenpos[n];
00529     *tuple++ = zgenpos[n];
00530     //
00531     rpos[n]  = Event.GetPositronVertexR(n);
00532     cpos[n]  = Event.GetPositronVertexCosTheta(n);
00533     fpos[n]  = Event.GetPositronVertexPhi(n);
00534     xpos[n]  = rpos[n]*cos(fpos[n])*sqrt(1-cpos[n]*cpos[n]);
00535     ypos[n]  = rpos[n]*sin(fpos[n])*sqrt(1-cpos[n]*cpos[n]);
00536     zpos[n]  = rpos[n]*cpos[n];
00537     //
00538     *tuple++ = rpos[n];
00539     *tuple++ = cpos[n];
00540     *tuple++ = fpos[n];
00541     *tuple++ = xpos[n];
00542     *tuple++ = ypos[n];
00543     *tuple++ = zpos[n];
00544     //
00545     *tuple++ = Event.GetPositronTime(n);
00546   }
00547   // zero out remaining tuples
00548   for(int n=npos;n<5;n++){
00549     for(int i=0;i<((NvarPos-1)/5);i++) *tuple++ = 0;
00550   }
00551   //
00552   hfn_(IdPos,TuplePos);
00553   //
00554   //
00555   // fill the neutron event variables
00556   //
00557   //cout << "NvarNeu " << NvarNeu << endl;
00558   float TupleNeu[NvarNeu];
00559   tuple = TupleNeu;
00560   //
00561   // get number of neutron
00562   int nneu = Event.GetNumNeutron();
00563   //cout << nneu << " neutron(s)" << endl;
00564   *tuple++ = nneu;
00565   //
00566   // deal with multiple neutrons
00567   double rgenneu[nneu];
00568   double cgenneu[nneu];
00569   double fgenneu[nneu];
00570   double xgenneu[nneu];
00571   double ygenneu[nneu];
00572   double zgenneu[nneu];
00573   //
00574   double rneu[nneu];
00575   double cneu[nneu];
00576   double fneu[nneu];
00577   double xneu[nneu];
00578   double yneu[nneu];
00579   double zneu[nneu];
00580   //
00581   for(int n=0;n<nneu;n++){
00582     *tuple++ = Event.GetNeutronKE(n);
00583     *tuple++ = Event.GetNeutronCosTheta(n);
00584     *tuple++ = Event.GetNeutronParentProcess(n);
00585     //
00586     rgenneu[n]  = Event.GetGeneratedNeutronVertexR(n);
00587     cgenneu[n]  = Event.GetGeneratedNeutronVertexCosTheta(n);
00588     fgenneu[n]  = Event.GetGeneratedNeutronVertexPhi(n);
00589     xgenneu[n]  = rgenneu[n]*cos(fgenneu[n])*sqrt(1-cgenneu[n]*cgenneu[n]);
00590     ygenneu[n]  = rgenneu[n]*sin(fgenneu[n])*sqrt(1-cgenneu[n]*cgenneu[n]);
00591     zgenneu[n]  = rgenneu[n]*cgenneu[n];
00592     //
00593     *tuple++ = rgenneu[n];
00594     *tuple++ = cgenneu[n]; 
00595     *tuple++ = fgenneu[n];
00596     *tuple++ = xgenneu[n];
00597     *tuple++ = ygenneu[n];
00598     *tuple++ = zgenneu[n];
00599     //
00600     rneu[n]  = Event.GetNeutronVertexR(n);
00601     cneu[n]  = Event.GetNeutronVertexCosTheta(n);
00602     fneu[n]  = Event.GetNeutronVertexPhi(n);
00603     xneu[n]  = rneu[n]*cos(fneu[n])*sqrt(1-cneu[n]*cneu[n]);
00604     yneu[n]  = rneu[n]*sin(fneu[n])*sqrt(1-cneu[n]*cneu[n]);
00605     zneu[n]  = rneu[n]*cneu[n];
00606     //
00607     *tuple++ = rneu[n];
00608     *tuple++ = cneu[n]; 
00609     *tuple++ = fneu[n];
00610     *tuple++ = xneu[n];
00611     *tuple++ = yneu[n];
00612     *tuple++ = zneu[n];
00613     //
00614     *tuple++ = Event.GetNeutronTime(n);
00615     *tuple++ = Detector.GetNeutronSteps(n);
00616     *tuple++ = Detector.GetAbsorber(n);
00617   }
00618   // zero out remaining tuples
00619   for(int n=nneu;n<5;n++){
00620     for(int i=0;i<((NvarNeu-1)/5);i++) *tuple++ = 0;
00621   }
00622   //
00623   hfn_(IdNeu,TupleNeu);
00624   //
00625   //
00626   // fill the PMT detector variables
00627   //
00628   //cout << "NvarDet " << NvarDet << endl;
00629   float TupleDet[NvarDet];
00630   tuple = TupleDet;
00631   //
00632   *tuple++ = Detector.GetNcosbins();
00633   *tuple++ = Detector.GetNphibins();
00634   //
00635   *tuple++ = Detector.GetPHpositron();
00636   *tuple++ = Detector.GetPHposReco();
00637   *tuple++ = Detector.GetPHneutron();
00638   *tuple++ = Detector.GetPHneuReco();
00639    //
00640   *tuple++ = Detector.GetMeanNpePos();
00641   *tuple++ = Detector.GetMeanAttenPos();
00642   *tuple++ = Detector.GetMeanSfactorPos();
00643   *tuple++ = Detector.GetMeanNpeNeu();
00644   *tuple++ = Detector.GetMeanAttenNeu();
00645   *tuple++ = Detector.GetMeanSfactorNeu();
00646   *tuple++ = Detector.GetNgamma('P'); 
00647   *tuple++ = Detector.GetNgamma('N');
00648   *tuple++ = Detector.GetNgamma('T');
00649   *tuple++ = Detector.GetHitFracPos();  
00650   *tuple++ = Detector.GetHitFracNeu();  
00651   double dippos[3];
00652   double dipneu[3];
00653   for(int i=0;i<3;i++)dippos[i]=*(Detector.GetPositronDipole()+i);
00654   for(int i=0;i<3;i++)dipneu[i]=*(Detector.GetNeutronDipole()+i);
00655   //
00656   *tuple++ = dippos[0];  
00657   *tuple++ = dippos[1];  
00658   *tuple++ = dippos[2];  
00659   *tuple++ = dipneu[0];  
00660   *tuple++ = dipneu[1];  
00661   *tuple++ = dipneu[2];  
00662   //
00663   hfn_(IdDet,TupleDet);
00664   //
00665   //
00666   // do some looping over all the "PMT"
00667   //
00668   /*
00669   int npmt=Detector.GetNpmt();
00670   for(int n=0;n<npmt;n++){
00671     int npos;
00672     int nneu;
00673     double xyzpmt[3];
00674     double phpos[100];
00675     double phneu[100];
00676     double tpos[100];
00677     double tneu[100];
00678     //
00679     // Get "phototube" information.
00680     //      
00681     Detector.GetPMTdata(n,xyzpmt,npos,phpos,tpos,nneu,phneu,tneu);
00682     
00683     for(int nn=0;nn<npos;nn++){
00684       hfill_(Offset+101,*(phpos+nn));
00685       hfill_(Offset+201,*(tpos+nn));
00686     }
00687     for(int nn=0;nn<nneu;nn++){
00688       hfill_(Offset+102,*(phneu+nn));
00689       hfill_(Offset+202,*(tneu+nn));
00690     }
00691   }
00692   //
00693   // Get information about gammas from positron.
00694   // Note that the positron dedx is defined as a "gamma"
00695   // originating at the positron range point.
00696   //
00697   int nposgamma = Detector.GetNgamma('P'); 
00698   hfill_(Offset+301,1.0*nposgamma);
00699   double egp[nposgamma];
00700   double tgp[nposgamma];
00701   double rgp[nposgamma];
00702   double cgp[nposgamma];
00703   double fgp[nposgamma];
00704   Detector.GetGamma('P',egp,tgp,rgp,cgp,fgp);
00705   for (int n=0;n<nposgamma;n++){
00706     hfill_(Offset+401,egp[n]);
00707     int n = 0;
00708     hfill_(Offset+701,tgp[n]-Event.GetPositronTime(n));
00709     hfill_(Offset+801,tgp[n]-Event.GetNubarTime());
00710   }
00711   double xyzp[3*nposgamma];
00712   Detector.GetGamma('P',xyzp);
00713   double* pxyz=xyzp;
00714   for (int n=0;n<nposgamma;n++){
00715     double* qxyz=Event.GetPositronVertexXYZ();
00716     double* rxyz=Event.GetNubarVertexXYZ();
00717     double dr1[3];
00718     double dr2[3];
00719     for(int i=0;i<3;i++){
00720       dr1[i] = *pxyz - *qxyz++;
00721       dr2[i] = *pxyz++ - *rxyz++;
00722     }
00723     hfill_(Offset+501,sqrt(dr1[0]*dr1[0]+dr1[1]*dr1[1]+dr1[2]*dr1[2]));
00724     hfill_(Offset+601,sqrt(dr2[0]*dr2[0]+dr2[1]*dr2[1]+dr2[2]*dr2[2]));
00725   }
00726   //
00727   // Get information about gammas from neutron.
00728   //
00729   int nneugamma = Detector.GetNgamma('N'); 
00730   hfill_(Offset+302,1.0*nneugamma);
00731   double egn[nneugamma];
00732   double tgn[nneugamma];
00733   double rgn[nneugamma];
00734   double cgn[nneugamma];
00735   double fgn[nneugamma];
00736   Detector.GetGamma('N',egn,tgn,rgn,cgn,fgn);
00737   for (int n=0;n<nneugamma;n++){
00738     hfill_(Offset+402,egn[n]);
00739     int n = 0;
00740     hfill_(Offset+702,tgn[n]-Event.GetNeutronTime(n));
00741     hfill_(Offset+802,tgn[n]-Event.GetNubarTime());
00742   }
00743   double xyzn[3*nneugamma];
00744   Detector.GetGamma('N',xyzn);
00745   pxyz=xyzn;
00746   for (int n=0;n<nneugamma;n++){
00747     double* qxyz=Event.GetNeutronVertexXYZ();
00748     double* rxyz=Event.GetNubarVertexXYZ();
00749     double dr1[3];
00750     double dr2[3];
00751     for(int i=0;i<3;i++){
00752       dr1[i] = *pxyz - *qxyz++;
00753       dr2[i] = *pxyz++ - *rxyz++;
00754     }
00755     hfill_(Offset+502,sqrt(dr1[0]*dr1[0]+dr1[1]*dr1[1]+dr1[2]*dr1[2]));
00756     hfill_(Offset+602,sqrt(dr2[0]*dr2[0]+dr2[1]*dr2[1]+dr2[2]*dr2[2]));
00757   }
00758   */
00759 }
00760