1616
1717#include "WCSimRootEvent.hh"
1818#include "WCSimRootGeom.hh"
19+ #include "WCSimRootOptions.hh"
1920#include "WCSimEnumerations.hh"
2021
2122using namespace std ;
@@ -30,7 +31,7 @@ TString create_filename(const char * suffix, TString& filename_string, const cha
3031}
3132
3233// Simple example of reading a generated Root file
33- int daq_readfile (const char * filename = NULL ,
34+ int daq_readfilemain (const char * filename = NULL ,
3435 bool verbose = false,
3536 const char * inbranch = "wcsimrootevent" ,
3637 Long64_t max_nevents = 999999999999 ,
@@ -90,6 +91,27 @@ int daq_readfile(const char *filename=NULL,
9091 exit (9 );
9192 }
9293 geotree -> GetEntry (0 );
94+
95+ // Options tree - only need 1 "event"
96+ TTree * opttree = (TTree * )file -> Get ("wcsimRootOptionsT" );
97+ WCSimRootOptions * opt = 0 ;
98+ opttree -> SetBranchAddress ("wcsimrootoptions" , & opt );
99+ if (verbose ) std ::cout << "Options tree has: " << opttree -> GetEntries () << " entries (1 expected)" << std ::endl ;
100+ if (opttree -> GetEntries () == 0 ) {
101+ exit (9 );
102+ }
103+ opttree -> GetEntry (0 );
104+ opt -> Print ();
105+ string detector = opt -> GetDetectorName ();
106+ //detector boundaries, in cm, for determining
107+ //These are setup for HyperK_HybridmPMT_WithOD_Realistic only
108+ const float detIDR = 3240 ;
109+ const float detODRin = 3302.1 ;
110+ const float detODRout = 3402.1 ;
111+ const float detIDhalfZ = 6579.1 / 2 ;
112+ const float detODhalfZin = 6699.3 / 2 ;
113+ const float detODhalfZout = 7099.5 / 2 ;
114+
93115 // start with the main "subevent", as it contains most of the info
94116 // and always exists.
95117 WCSimRootTrigger * wcsimrootevent ;
@@ -133,6 +155,7 @@ int daq_readfile(const char *filename=NULL,
133155 vector < float > tvtrack_energy ;
134156 vector < double > tvtrack_time ;
135157 vector < TVector3 > tvtrack_startpos , tvtrack_stoppos ;
158+ int true_particles_in_od , true_particles_in_id ;
136159 tout -> Branch ("vtx0" , & tvtx0 );
137160 tout -> Branch ("vtx1" , & tvtx1 );
138161 tout -> Branch ("vtx2" , & tvtx2 );
@@ -147,6 +170,8 @@ int daq_readfile(const char *filename=NULL,
147170 tout -> Branch ("track_flag" , & tvtrack_flag );
148171 tout -> Branch ("track_energy" , & tvtrack_energy );
149172 tout -> Branch ("track_time" , & tvtrack_time );
173+ tout -> Branch ("true_particles_in_od_NOT_CHERENKOV_OR_ELECTRON" , & true_particles_in_od );
174+ tout -> Branch ("true_particles_in_id_NOT_CHERENKOV_OR_ELECTRON" , & true_particles_in_id );
150175 //tout->Branch("track_startpos", &tvtrack_startpos);
151176 //tout->Branch("track_stoppos", &tvtrack_stoppos);
152177 //per trigger variables
@@ -194,6 +219,8 @@ int daq_readfile(const char *filename=NULL,
194219 tvtx0 = -999999 ;
195220 tvtx1 = -999999 ;
196221 tvtx2 = -999999 ;
222+ true_particles_in_id = -1 ;
223+ true_particles_in_od = -1 ;
197224 tntriggers = -1 ;
198225 tntracks = -1 ;
199226 tnrawhits = -1 ;
@@ -256,6 +283,10 @@ int daq_readfile(const char *filename=NULL,
256283 tntracks = ntracks ;
257284
258285 // Loop through elements in the TClonesArray of WCSimTracks
286+ if (detector == "HyperK_HybridmPMT_WithOD_Realistic" ) {
287+ true_particles_in_id = 0 ;
288+ true_particles_in_od = 0 ;
289+ }
259290 for (int itrack = 0 ; itrack < ntracks ; itrack ++ )
260291 {
261292 TObject * element = (wcsimrootevent -> GetTracks ())-> At (itrack );
@@ -270,13 +301,35 @@ int daq_readfile(const char *filename=NULL,
270301 tvtrack_stoppos .push_back (TVector3 (wcsimroottrack -> GetStop (0 ),
271302 wcsimroottrack -> GetStop (1 ),
272303 wcsimroottrack -> GetStop (2 )));
304+ if (detector == "HyperK_HybridmPMT_WithOD_Realistic"
305+ && wcsimroottrack -> GetIpnu () != 11 // ignore electrons
306+ && wcsimroottrack -> GetIpnu () != 0 // ignore optical photons
307+ && wcsimroottrack -> GetFlag () != -2 //ignore fake nuclei tracks
308+ ) {
309+ float trueRstart = TMath ::Sqrt (wcsimroottrack -> GetStart (0 ) * wcsimroottrack -> GetStart (0 ) + wcsimroottrack -> GetStart (1 ) * wcsimroottrack -> GetStart (1 ));
310+ float trueZstart = wcsimroottrack -> GetStart (2 );
311+ float trueRstop = TMath ::Sqrt (wcsimroottrack -> GetStop (0 ) * wcsimroottrack -> GetStop (0 ) + wcsimroottrack -> GetStop (1 ) * wcsimroottrack -> GetStop (1 ));
312+ float trueZstop = wcsimroottrack -> GetStop (2 );
313+ if ((trueRstart <= detIDR && abs (trueZstart ) <= detIDhalfZ ) ||
314+ (trueRstop <= detIDR && abs (trueZstop ) <= detIDhalfZ )) {
315+ true_particles_in_id ++ ;
316+ }
317+ if ((trueRstart <= detODRout && abs (trueZstart ) <= detODhalfZout && trueRstart >= detODRin && abs (trueZstart ) >= detODhalfZin ) ||
318+ (trueRstop <= detODRout && abs (trueZstop ) <= detODhalfZout && trueRstop >= detODRin && abs (trueZstop ) >= detODhalfZin ))
319+ true_particles_in_od ++ ;
320+ }
273321 if (verbose ){
274322 printf ("Track ipnu: %d\n" ,wcsimroottrack -> GetIpnu ());
275323 printf ("Track parent ID: %d\n" ,wcsimroottrack -> GetParenttype ());
276324 for (int j = 0 ; j < 3 ; j ++ )
277325 printf ("Track dir: %d %f\n" , j , wcsimroottrack -> GetDir (j ));
278326 }
279327 } //itrack // End of loop over tracks
328+ if (detector == "HyperK_HybridmPMT_WithOD_Realistic" ) {
329+ //treat as a bools
330+ if (true_particles_in_id > 0 ) true_particles_in_id = 1 ;
331+ if (true_particles_in_od > 0 ) true_particles_in_od = 1 ;
332+ }
280333
281334 //get number of hits and digits
282335 const int ncherenkovhits = wcsimrootevent -> GetNcherenkovhits ();
@@ -570,7 +623,7 @@ int main(int argc, char *argv[]){
570623 usage (branches );
571624 return 1 ;
572625 }
573- daq_readfile (argv [1 ], atoi (argv [2 ]), branches [jobtodo ].c_str ());
626+ daq_readfilemain (argv [1 ], atoi (argv [2 ]), branches [jobtodo ].c_str ());
574627
575628 return 0 ;
576629
0 commit comments