read_t2k_rootracker.C
Go to the documentation of this file.
1 //________________________________________________________________________________________
2 /*!
3 
4 \macro read_t2k_rootracker.C
5 
6 \brief Test macro to read a GENIE event tree in the t2k_rootracker format
7  and print-out the neutrino event and flux pass-through info.
8 
9 \usage shell% root
10  root[0] .L read_t2k_rootracker.C++
11  root[1] read_t2k_rootracker("./your_rootracker_file.root");
12 
13 \author Costas Andreopoulos <costas.andreopoulos@stfc.ac.uk>
14  University of Liverpool & STFC Rutherford Appleton Lab
15 
16 \created Nov 24, 2008
17 
18 \cpright Copyright (c) 2003-2017, GENIE Neutrino MC Generator Collaboration
19  For the full text of the license visit http://copyright.genie-mc.org
20  or see $GENIE/LICENSE
21 */
22 //_________________________________________________________________________________________
23 
24 #include <iostream>
25 
26 #include <TFile.h>
27 #include <TTree.h>
28 #include <TBranch.h>
29 #include <TBits.h>
30 #include <TObjString.h>
31 
32 using std::cout;
33 using std::endl;
34 
35 void read_t2k_rootracker(const char * filename)
36 {
37  bool using_new_version = false; // StdHepReScat and G2NeutEvtCode branches available only for versions >= 2.5.1
38 
39  const int kNPmax = 100;
40 
41  TFile file(filename, "READ");
42  TTree * tree = (TTree *) file.Get("gRooTracker");
43  assert(tree);
44 
45  TBits* EvtFlags = 0; // generator-specific event flags
46  TObjString* EvtCode = 0; // generator-specific string with 'event code'
47  int EvtNum; // event num.
48  double EvtXSec; // cross section for selected event (1E-38 cm2)
49  double EvtDXSec; // cross section for selected event kinematics (1E-38 cm2 /{K^n})
50  double EvtWght; // weight for that event
51  double EvtProb; // probability for that event (given cross section, path lengths, etc)
52  double EvtVtx[4]; // event vertex position in detector coord syst (in geom units)
53  int StdHepN; // number of particles in particle array
54  int StdHepPdg [kNPmax]; // stdhep-like particle array: pdg codes (& generator specific codes for pseudoparticles)
55  int StdHepStatus[kNPmax]; // stdhep-like particle array: generator-specific status code
56  int StdHepRescat[kNPmax]; // stdhep-like particle array: intranuclear rescattering code [ >= v2.5.1 ]
57  double StdHepX4 [kNPmax][4]; // stdhep-like particle array: 4-x (x, y, z, t) of particle in hit nucleus frame (fm)
58  double StdHepP4 [kNPmax][4]; // stdhep-like particle array: 4-p (px,py,pz,E) of particle in LAB frame (GeV)
59  double StdHepPolz [kNPmax][3]; // stdhep-like particle array: polarization vector
60  int StdHepFd [kNPmax]; // stdhep-like particle array: first daughter
61  int StdHepLd [kNPmax]; // stdhep-like particle array: last daughter
62  int StdHepFm [kNPmax]; // stdhep-like particle array: first mother
63  int StdHepLm [kNPmax]; // stdhep-like particle array: last mother
64  int G2NeutEvtCode; // NEUT code for the current GENIE event [ >= v2.5.1 ]
65  int NuParentPdg; // parent hadron pdg code
66  int NuParentDecMode; // parent hadron decay mode
67  double NuParentDecP4 [4]; // parent hadron 4-momentum at decay
68  double NuParentDecX4 [4]; // parent hadron 4-position at decay
69  double NuParentProP4 [4]; // parent hadron 4-momentum at production
70  double NuParentProX4 [4]; // parent hadron 4-position at production
71  int NuParentProNVtx; // parent hadron vtx id
72 
73  // get branches
74  TBranch * brEvtFlags = tree -> GetBranch ("EvtFlags");
75  TBranch * brEvtCode = tree -> GetBranch ("EvtCode");
76  TBranch * brEvtNum = tree -> GetBranch ("EvtNum");
77  TBranch * brEvtXSec = tree -> GetBranch ("EvtXSec");
78  TBranch * brEvtDXSec = tree -> GetBranch ("EvtDXSec");
79  TBranch * brEvtWght = tree -> GetBranch ("EvtWght");
80  TBranch * brEvtProb = tree -> GetBranch ("EvtProb");
81  TBranch * brEvtVtx = tree -> GetBranch ("EvtVtx");
82  TBranch * brStdHepN = tree -> GetBranch ("StdHepN");
83  TBranch * brStdHepPdg = tree -> GetBranch ("StdHepPdg");
84  TBranch * brStdHepStatus = tree -> GetBranch ("StdHepStatus");
85  TBranch * brStdHepRescat = (using_new_version) ? tree -> GetBranch ("StdHepRescat") : 0;
86  TBranch * brStdHepX4 = tree -> GetBranch ("StdHepX4");
87  TBranch * brStdHepP4 = tree -> GetBranch ("StdHepP4");
88  TBranch * brStdHepPolz = tree -> GetBranch ("StdHepPolz");
89  TBranch * brStdHepFd = tree -> GetBranch ("StdHepFd");
90  TBranch * brStdHepLd = tree -> GetBranch ("StdHepLd");
91  TBranch * brStdHepFm = tree -> GetBranch ("StdHepFm");
92  TBranch * brStdHepLm = tree -> GetBranch ("StdHepLm");
93  TBranch * brG2NeutEvtCode = (using_new_version) ? tree -> GetBranch ("G2NeutEvtCode") : 0;
94  TBranch * brNuParentPdg = tree -> GetBranch ("NuParentPdg");
95  TBranch * brNuParentDecMode = tree -> GetBranch ("NuParentDecMode");
96  TBranch * brNuParentDecP4 = tree -> GetBranch ("NuParentDecP4");
97  TBranch * brNuParentDecX4 = tree -> GetBranch ("NuParentDecX4");
98  TBranch * brNuParentProP4 = tree -> GetBranch ("NuParentProP4");
99  TBranch * brNuParentProX4 = tree -> GetBranch ("NuParentProX4");
100  TBranch * brNuParentProNVtx = tree -> GetBranch ("NuParentProNVtx");
101 
102  // set address
103  brEvtFlags -> SetAddress ( &EvtFlags );
104  brEvtCode -> SetAddress ( &EvtCode );
105  brEvtNum -> SetAddress ( &EvtNum );
106  brEvtXSec -> SetAddress ( &EvtXSec );
107  brEvtDXSec -> SetAddress ( &EvtDXSec );
108  brEvtWght -> SetAddress ( &EvtWght );
109  brEvtProb -> SetAddress ( &EvtProb );
110  brEvtVtx -> SetAddress ( EvtVtx );
111  brStdHepN -> SetAddress ( &StdHepN );
112  brStdHepPdg -> SetAddress ( StdHepPdg );
113  brStdHepStatus -> SetAddress ( StdHepStatus );
114  if(using_new_version) {
115  brStdHepRescat -> SetAddress ( StdHepRescat );
116  }
117  brStdHepX4 -> SetAddress ( StdHepX4 );
118  brStdHepP4 -> SetAddress ( StdHepP4 );
119  brStdHepPolz -> SetAddress ( StdHepPolz );
120  brStdHepFd -> SetAddress ( StdHepFd );
121  brStdHepLd -> SetAddress ( StdHepLd );
122  brStdHepFm -> SetAddress ( StdHepFm );
123  brStdHepLm -> SetAddress ( StdHepLm );
124  if(using_new_version) {
125  brG2NeutEvtCode -> SetAddress ( &G2NeutEvtCode );
126  }
127  brNuParentPdg -> SetAddress ( &NuParentPdg );
128  brNuParentDecMode -> SetAddress ( &NuParentDecMode );
129  brNuParentDecP4 -> SetAddress ( NuParentDecP4 );
130  brNuParentDecX4 -> SetAddress ( NuParentDecX4 );
131  brNuParentProP4 -> SetAddress ( NuParentProP4 );
132  brNuParentProX4 -> SetAddress ( NuParentProX4 );
133  brNuParentProNVtx -> SetAddress ( &NuParentProNVtx );
134 
135  int n = tree->GetEntries();
136  printf("Number of entries: %d", n);
137 
138  // read tree & print some info for each entry
139 
140  for(int i=0; i < tree->GetEntries(); i++) {
141  printf("\n\n ** Current entry: %d \n", i);
142  tree->GetEntry(i);
143 
144  printf("\n -----------------------------------------------------------------------------------------------------------------");
145  printf("\n Event code : %s", EvtCode->String().Data());
146  printf("\n Event x-section : %10.5f * 1E-38* cm^2", EvtXSec);
147  printf("\n Event kinematics x-section : %10.5f * 1E-38 * cm^2/{K^n}", EvtDXSec);
148  printf("\n Event weight : %10.8f", EvtWght);
149  printf("\n Event vertex : x = %8.2f mm, y = %8.2f mm, z = %8.2f mm", EvtVtx[0], EvtVtx[1], EvtVtx[2]);
150  printf("\n *Particle list:");
151  printf("\n --------------------------------------------------------------------------------------------------------------------------");
152  printf("\n | Idx | Ist | PDG | Rescat | Mother | Daughter | Px | Py | Pz | E | x | y | z |");
153  printf("\n | | | | | | |(GeV/c) |(GeV/c) |(GeV/c) | (GeV) | (fm) | (fm) | (fm) |");
154  printf("\n --------------------------------------------------------------------------------------------------------------------------");
155 
156  for(int ip=0; ip<StdHepN; ip++) {
157  printf("\n | %3d | %3d | %10d | %6d | %3d | %3d | %3d | %3d | %6.3f | %6.3f | %6.3f | %6.3f | %6.3f | %6.3f | %6.3f |",
158  ip, StdHepStatus[ip], StdHepPdg[ip], StdHepRescat[ip],
159  StdHepFm[ip], StdHepLm[ip], StdHepFd[ip], StdHepLd[ip],
160  StdHepP4[ip][0], StdHepP4[ip][1], StdHepP4[ip][2], StdHepP4[ip][3],
161  StdHepX4[ip][0], StdHepX4[ip][1], StdHepX4[ip][2]);
162  }
163  printf("\n --------------------------------------------------------------------------------------------------------------------------");
164  printf("\n *Flux Info:");
165  printf("\n Parent hadron pdg code : %d", NuParentPdg);
166  printf("\n Parent hadron decay mode : %d", NuParentDecMode);
167  printf("\n Parent hadron 4p at decay : Px = %6.3f GeV/c, Py = %6.3f GeV/c, Pz = %6.3f GeV/c, E = %6.3f GeV",
168  NuParentDecP4[0], NuParentDecP4[1], NuParentDecP4[2], NuParentDecP4[3]);
169  printf("\n Parent hadron 4p at prod. : Px = %6.3f GeV/c, Py = %6.3f GeV/c, Pz = %6.3f GeV/c, E = %6.3f GeV",
170  NuParentProP4[0], NuParentProP4[1], NuParentProP4[2], NuParentProP4[3]);
171  printf("\n Parent hadron 4x at decay : x = %6.3f m, y = %6.3f m, z = %6.3f m, t = %6.3f s",
172  NuParentDecX4[0], NuParentDecX4[1], NuParentDecX4[2], NuParentDecX4[3]);
173  printf("\n Parent hadron 4x at prod. : x = %6.3f m, y = %6.3f m, z = %6.3f m, t = %6.3f s",
174  NuParentProX4[0], NuParentProX4[1], NuParentProX4[2], NuParentProX4[3]);
175  printf("\n -------------------------------------------------------------------------------------------------------------------------- \n");
176 
177  } // tree entries
178 
179  file.Close();
180 }
size_t i(0)
const int kNPmax
Definition: gNtpConv.cxx:228
TString ip
Definition: loadincs.C:5
the ParameterSet object passed in for the configuration of a destination should be the only source that can affect the behavior of that destination This is to eliminate the dependencies of configuring a destination from multiple mostly from the defaults It suppresses possible glitches about changing the configuration file somewhere outside of a destination segament might still affect the behavior of that destination In the previous configuration for a specific the value of a certain e may come from following and have been suppressed It the configuring ParameterSet object for each destination will be required to carry a parameter list as complete as possible If a parameter still cannot be found in the ParameterSet the configuration code will go look for a hardwired default directly The model is a great simplicity comparing with the previous especially when looking for default values Another great advantage is most of the parameters now have very limited places that allows to appear Usually they can only appear at one certain level in a configuration file For in the old configuring model or in a default ParameterSet object inside of a or in a category or in a severity object This layout of multiple sources for a single parameter brings great confusion in both writing a configuration and in processing the configuration file Under the new the only allowed place for the parameter limit to appear is inside of a category which is inside of a destination object Other improvements simplify the meaning of a destination name In the old a destination name has multiple folds of meanings the e cout and cerr have the special meaning of logging messages to standard output or standard error the name also serves as the output filename if the destination is a file these names are also references to look up for detailed configurations in configuring the MessageFacility The multi purpose of the destination name might cause some unwanted behavior in either writing or parsing the configuration file To amend in the new model the destination name is now merely a name for a which might represent the literal purpose of this or just an id All other meanings of the destinations names now go into the destination ParameterSet as individual such as the type parameter and filename parameter Following is the deatiled rule for the new configuring Everything that is related with MessageFacility configuration must be wrapped in a single ParameterSet object with the name MessageFacility The MessageFacility ParameterSet object contains a series of top level parameters These parameters can be chosen a vector of string listing the name of debug enabled models Or use *to enable debug messages in all modules a vector of string a vector of string a vector of string a ParameterSet object containing the list of all destinations The destinations ParameterSet object is a combination of ParameterSet objects for individual destinations There are two types of destinations that you can insert in the destinations ParameterSet ordinary including cout
*file AnalysisTree_module cc *brief Module to create a TTree for analysis *authors tjyang fnal sowjanyag phys ksu edu **Taken from uboone code Imported by Karl k warburton sheffield ac uk *with the help of Tingjun Yang **Current implementation with one set of branches for each tracking algorithm *The data structure which hosts the addresses of the tree branches is *dynamically allocated on and it can be optionally destroyed at the *end of each event *The data and it is contained in a C vector of *one per algorithm These structures can also be allocated on demand *Each of these structures is connected to a set of one branch per *data member Data members are vectors of numbers or vectors of fixed size *C arrays The vector index represents the tracks reconstructed by the and each has a fixed size pool for connect to a *ROOT tree(creating the branches they need) and resize.*The AnalysisTreeDataStruct is const ructed with as many tracking algorithms as *there are named in the module configuration(even if they are not backed by *any available tracking data).*By default const ruction
< separator(=)> module_type Type Source location< separator(-)> DummyAnalyzer analyzer< path > DummyAnalyzer_module cc DummyFilter filter< path > DummyFilter_module cc *DummyProducer producer< path > DummyProducer_module cc *DummyProducer producer< path > DummyProducer_module cc< separator(=)> The modules marked *above are degenerate i e specifying the short module_type value leads to an ambiguity In order to use a degenerate in your configuration file
void read_t2k_rootracker(const char *filename)