gNtpConv.cxx
Go to the documentation of this file.
1 //_____________________________________________________________________________________________
2 /*!
3 
4 \program gntpc
5 
6 \brief Converts a native GENIE (GHEP/ROOT) event tree file to a host of
7  plain text, XML or bare-ROOT formats.
8 
9  Syntax:
10  gntpc -i input_file [-o output_file] -f format [-n nev] [-v vrs] [-c]
11  [--seed random_number_seed]
12  [--message-thresholds xml_file]
13  [--event-record-print-level level]
14 
15 
16  Options :
17 
18  [] denotes an optional argument
19 
20  -n
21  Number of events to convert
22  (optional, default: convert all events)
23  -v
24  Output format version, if multiple versions are supported
25  (optional, default: use latest version of each format)
26  -c
27  Copy MC job metadata (gconfig and genv TFolders) from the input GHEP file.
28  -f
29  A string that specifies the output file format.
30  >>
31  >> Generic formats:
32  >>
33  * `gst':
34  The 'definite' GENIE summary tree format (gst).
35  * `gxml':
36  GENIE XML event format
37  * `ghep_mock_data':
38  Output file has the same format as the input file (GHEP) but
39  all information other than final state particles is hidden
40  * `rootracker':
41  A bare-ROOT STDHEP-like GENIE event tree.
42  * `rootracker_mock_data':
43  As the `rootracker' format but hiddes all information
44  except the final state particles.
45  >>
46  >> Experiment-specific formats:
47  >>
48  * `t2k_rootracker':
49  A variance of the `rootracker' format used by the nd280, INGRID and 2km.
50  Includes, in addition, JPARC flux pass-through info.
51  * `numi_rootracker':
52  A variance of the `rootracker' format for the NuMI expts.
53  Includes, in addition, NuMI flux pass-through info.
54  * `t2k_tracker':
55  A tracker-type format with tweaks required by the SuperK MC (SKDETSIM):
56  - Converting K0, \bar{K0} -> KO_{long}, K0_{short}
57  - Emulating 'NEUT' reaction codes
58  - Appropriate $track ordering for SKDETSIM
59  - Passing detailed GENIE MC truth and JPARC flux info
60  using the tracker $info lines. This information,
61  propaged by SKDETSIM to the DSTs, is identical with the
62  one used at the near detectors and can be used for
63  global systematic studies.
64  >>
65  >> GENIE test / cross-generator comparison formats:
66  >>
67  * `ghad':
68  NEUGEN-style text-based format for hadronization studies
69  * `ginuke':
70  A summary ntuple for intranuclear-rescattering studies using simulated
71  hadron-nucleus samples
72  >>
73  >> Other (depreciated) formats:
74  >>
75  * `nuance_tracker':
76  NUANCE-style tracker text-based format
77  -o
78  Specifies the output filename.
79  If not specified a the default filename is constructed by the
80  input base name and an extension depending on the file format:
81  `gst' -> *.gst.root
82  `gxml' -> *.gxml
83  `ghep_mock_data' -> *.mockd.ghep.root
84  `rootracker' -> *.gtrac.root
85  `rootracker_mock_data' -> *.mockd.gtrac.root
86  `t2k_rootracker' -> *.gtrac.root
87  `numi_rootracker' -> *.gtrac.root
88  `t2k_tracker' -> *.gtrac.dat
89  `nuance_tracker' -> *.gtrac_legacy.dat
90  `ghad' -> *.ghad.dat
91  `ginuke' -> *.ginuke.root
92  --seed
93  Random number seed.
94  --message-thresholds
95  Allows users to customize the message stream thresholds.
96  The thresholds are specified using an XML file.
97  See $GENIE/config/Messenger.xml for the XML schema.
98  --event-record-print-level
99  Allows users to set the level of information shown when the event
100  record is printed in the screen. See GHepRecord::Print().
101 
102  Examples:
103  (1) shell% gntpc -i myfile.ghep.root -f t2k_rootracker
104 
105  Converts all events in the GHEP file myfile.ghep.root into the
106  t2k_rootracker format.
107  The output file is named myfile.gtrac.root
108 
109 \author Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
110  University of Liverpool & STFC Rutherford Appleton Laboratory
111 
112 \created September 23, 2005
113 
114 \cpright Copyright (c) 2003-2020, The GENIE Collaboration
115  For the full text of the license visit http://copyright.genie-mc.org
116 
117 */
118 //_____________________________________________________________________________________________
119 
120 #include <cassert>
121 #include <string>
122 #include <sstream>
123 #include <fstream>
124 #include <iomanip>
125 #include <vector>
126 #include <algorithm>
127 
128 #include "libxml/parser.h"
129 #include "libxml/xmlmemory.h"
130 
131 #include <TSystem.h>
132 #include <TFile.h>
133 #include <TTree.h>
134 #include <TFolder.h>
135 #include <TBits.h>
136 #include <TObjString.h>
137 #include <TMath.h>
138 #include "Framework/ParticleData/BaryonResonance.h"
139 #include "Framework/ParticleData/BaryonResUtils.h"
140 #include "Framework/Conventions/GBuild.h"
141 #include "Framework/Conventions/Constants.h"
142 #include "Framework/Conventions/Units.h"
143 #include "Framework/EventGen/EventRecord.h"
144 #include "Framework/GHEP/GHepStatus.h"
145 #include "Framework/GHEP/GHepParticle.h"
146 #include "Framework/GHEP/GHepUtils.h"
147 #include "Framework/Ntuple/NtpMCFormat.h"
148 #include "Framework/Ntuple/NtpMCTreeHeader.h"
149 #include "Framework/Ntuple/NtpMCEventRecord.h"
150 #include "Framework/Ntuple/NtpWriter.h"
151 #include "Framework/Numerical/RandomGen.h"
152 #include "Framework/Messenger/Messenger.h"
153 #include "Framework/ParticleData/PDGCodes.h"
154 #include "Framework/ParticleData/PDGUtils.h"
155 #include "Framework/ParticleData/PDGLibrary.h"
156 #include "Framework/Utils/AppInit.h"
157 #include "Framework/Utils/RunOpt.h"
158 #include "Framework/Utils/CmdLnArgParser.h"
159 #include "Framework/Utils/SystemUtils.h"
160 #include "Framework/Utils/T2KEvGenMetaData.h"
161 
162 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
163 #include "Tools/Flux/GJPARCNuFlux.h"
164 #include "Tools/Flux/GNuMIFlux.h"
165 #endif
166 
167 //define __GHAD_NTP__
168 
169 using std::string;
170 using std::ostringstream;
171 using std::ofstream;
172 using std::endl;
173 using std::setw;
174 using std::setprecision;
175 using std::setfill;
176 using std::ios;
177 using std::setiosflags;
178 using std::vector;
179 
180 using namespace genie;
181 using namespace genie::constants;
182 
183 //func prototypes
184 void ConvertToGST (void);
185 void ConvertToGXML (void);
186 void ConvertToGHepMock (void);
187 void ConvertToGTracker (void);
188 void ConvertToGRooTracker (void);
189 void ConvertToGHad (void);
190 void ConvertToGINuke (void);
191 void GetCommandLineArgs (int argc, char ** argv);
192 void PrintSyntax (void);
193 string DefaultOutputFile (void);
194 int LatestFormatVersionNumber (void);
195 bool CheckRootFilename (string filename);
196 int HAProbeFSI (int, int, int, double [], int [], int, int, int); //Test code
197 //format enum
198 typedef enum EGNtpcFmt {
199  kConvFmt_undef = 0,
200  kConvFmt_gst,
201  kConvFmt_gxml,
202  kConvFmt_ghep_mock_data,
203  kConvFmt_rootracker,
204  kConvFmt_rootracker_mock_data,
205  kConvFmt_t2k_rootracker,
206  kConvFmt_numi_rootracker,
207  kConvFmt_t2k_tracker,
208  kConvFmt_nuance_tracker,
209  kConvFmt_ghad,
210  kConvFmt_ginuke
211 } GNtpcFmt_t;
212 
213 //input options (from command line arguments):
214 string gOptInpFileName; ///< input file name
215 string gOptOutFileName; ///< output file name
216 GNtpcFmt_t gOptOutFileFormat; ///< output file format id
217 int gOptVersion; ///< output file format version
218 Long64_t gOptN; ///< number of events to process
219 bool gOptCopyJobMeta = false; ///< copy MC job metadata (gconfig, genv TFolders)
220 long int gOptRanSeed; ///< random number seed
221 
222 //genie version used to generate the input event file
223 int gFileMajorVrs = -1;
224 int gFileMinorVrs = -1;
225 int gFileRevisVrs = -1;
226 
227 //consts
228 const int kNPmax = 250;
229 //____________________________________________________________________________________
230 int main(int argc, char ** argv)
231 {
232  GetCommandLineArgs(argc, argv);
233 
234  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
235  utils::app_init::RandGen(gOptRanSeed);
236 
237  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
238 
239  PDGLibrary::Instance()->AddDarkMatter( 1.0, 0.5 ) ;
240 
241  // Call the appropriate conversion function
242  switch(gOptOutFileFormat) {
243 
244  case (kConvFmt_gst) :
245 
246  ConvertToGST();
247  break;
248 
249  case (kConvFmt_gxml) :
250 
251  ConvertToGXML();
252  break;
253 
254  case (kConvFmt_ghep_mock_data) :
255 
256  ConvertToGHepMock();
257  break;
258 
259  case (kConvFmt_rootracker ) :
260  case (kConvFmt_rootracker_mock_data) :
261  case (kConvFmt_t2k_rootracker ) :
262  case (kConvFmt_numi_rootracker ) :
263 
264  ConvertToGRooTracker();
265  break;
266 
267  case (kConvFmt_t2k_tracker ) :
268  case (kConvFmt_nuance_tracker) :
269 
270  ConvertToGTracker();
271  break;
272 
273  case (kConvFmt_ghad) :
274 
275  ConvertToGHad();
276  break;
277 
278  case (kConvFmt_ginuke) :
279 
280  ConvertToGINuke();
281  break;
282 
283  default:
284  LOG("gntpc", pFATAL)
285  << "Invalid output format [" << gOptOutFileFormat << "]";
286  PrintSyntax();
287  gAbortingInErr = true;
288  exit(3);
289  }
290  return 0;
291 }
292 //____________________________________________________________________________________
293 // GENIE GHEP EVENT TREE FORMAT -> GENIE SUMMARY NTUPLE
294 //____________________________________________________________________________________
295 void ConvertToGST(void)
296 {
297  // Some constants
298  const double e_h = 1.3; // typical e/h ratio used for computing mean `calorimetric response'
299 
300  // Define branch variables
301  //
302  int brIev = 0; // Event number
303  int brNeutrino = 0; // Neutrino pdg code
304  int brFSPrimLept = 0; // Final state primary lepton pdg code
305  int brTarget = 0; // Nuclear target pdg code (10LZZZAAAI)
306  int brTargetZ = 0; // Nuclear target Z (extracted from pdg code above)
307  int brTargetA = 0; // Nuclear target A (extracted from pdg code above)
308  int brHitNuc = 0; // Hit nucleon pdg code (not set for COH,IMD and NuEL events)
309  int brHitQrk = 0; // Hit quark pdg code (set for DIS events only)
310  bool brFromSea = false; // Hit quark is from sea (set for DIS events only)
311  int brResId = 0; // Produced baryon resonance (set for resonance events only)
312  bool brIsQel = false; // Is QEL?
313  bool brIsRes = false; // Is RES?
314  bool brIsDis = false; // Is DIS?
315  bool brIsCoh = false; // Is Coherent?
316  bool brIsMec = false; // Is MEC?
317  bool brIsDfr = false; // Is Diffractive?
318  bool brIsImd = false; // Is IMD?
319  bool brIsSingleK = false; // Is single kaon?
320  bool brIsImdAnh = false; // Is IMD annihilation?
321  bool brIsNuEL = false; // Is ve elastic?
322  bool brIsEM = false; // Is EM process?
323  bool brIsCC = false; // Is Weak CC process?
324  bool brIsNC = false; // Is Weak NC process?
325  bool brIsCharmPro = false; // Produces charm?
326  bool brIsAMNuGamma = false; // is anomaly mediated nu gamma
327  int brCodeNeut = 0; // The equivalent NEUT reaction code (if any)
328  int brCodeNuance = 0; // The equivalent NUANCE reaction code (if any)
329  double brWeight = 0; // Event weight
330  double brKineXs = 0; // Bjorken x as was generated during kinematical selection; takes fermi momentum / off-shellness into account
331  double brKineYs = 0; // Inelasticity y as was generated during kinematical selection; takes fermi momentum / off-shellness into account
332  double brKineTs = 0; // Energy transfer to nucleus at COH events as was generated during kinematical selection
333  double brKineQ2s = 0; // Momentum transfer Q^2 as was generated during kinematical selection; takes fermi momentum / off-shellness into account
334  double brKineWs = 0; // Hadronic invariant mass W as was generated during kinematical selection; takes fermi momentum / off-shellness into account
335  double brKineX = 0; // Experimental-like Bjorken x; neglects fermi momentum / off-shellness
336  double brKineY = 0; // Experimental-like inelasticity y; neglects fermi momentum / off-shellness
337  double brKineT = 0; // Experimental-like energy transfer to nucleus at COH events
338  double brKineQ2 = 0; // Experimental-like momentum transfer Q^2; neglects fermi momentum / off-shellness
339  double brKineW = 0; // Experimental-like hadronic invariant mass W; neglects fermi momentum / off-shellness
340  double brEvRF = 0; // Neutrino energy @ the rest-frame of the hit-object (eg nucleon for CCQE, e- for ve- elastic,...)
341  double brEv = 0; // Neutrino energy @ LAB
342  double brPxv = 0; // Neutrino px @ LAB
343  double brPyv = 0; // Neutrino py @ LAB
344  double brPzv = 0; // Neutrino pz @ LAB
345  double brEn = 0; // Initial state hit nucleon energy @ LAB
346  double brPxn = 0; // Initial state hit nucleon px @ LAB
347  double brPyn = 0; // Initial state hit nucleon py @ LAB
348  double brPzn = 0; // Initial state hit nucleon pz @ LAB
349  double brEl = 0; // Final state primary lepton energy @ LAB
350  double brPxl = 0; // Final state primary lepton px @ LAB
351  double brPyl = 0; // Final state primary lepton py @ LAB
352  double brPzl = 0; // Final state primary lepton pz @ LAB
353  double brPl = 0; // Final state primary lepton p @ LAB
354  double brCosthl = 0; // Final state primary lepton cos(theta) wrt to neutrino direction
355  int brNfP = 0; // Nu. of final state p's + \bar{p}'s (after intranuclear rescattering)
356  int brNfN = 0; // Nu. of final state n's + \bar{n}'s
357  int brNfPip = 0; // Nu. of final state pi+'s
358  int brNfPim = 0; // Nu. of final state pi-'s
359  int brNfPi0 = 0; // Nu. of final state pi0's (
360  int brNfKp = 0; // Nu. of final state K+'s
361  int brNfKm = 0; // Nu. of final state K-'s
362  int brNfK0 = 0; // Nu. of final state K0's + \bar{K0}'s
363  int brNfEM = 0; // Nu. of final state gammas and e-/e+
364  int brNfOther = 0; // Nu. of heavier final state hadrons (D+/-,D0,Ds+/-,Lamda,Sigma,Lamda_c,Sigma_c,...)
365  int brNiP = 0; // Nu. of `primary' (: before intranuclear rescattering) p's + \bar{p}'s
366  int brNiN = 0; // Nu. of `primary' n's + \bar{n}'s
367  int brNiPip = 0; // Nu. of `primary' pi+'s
368  int brNiPim = 0; // Nu. of `primary' pi-'s
369  int brNiPi0 = 0; // Nu. of `primary' pi0's
370  int brNiKp = 0; // Nu. of `primary' K+'s
371  int brNiKm = 0; // Nu. of `primary' K-'s
372  int brNiK0 = 0; // Nu. of `primary' K0's + \bar{K0}'s
373  int brNiEM = 0; // Nu. of `primary' gammas and e-/e+
374  int brNiOther = 0; // Nu. of other `primary' hadron shower particles
375  int brNf = 0; // Nu. of final state particles in hadronic system
376  int brPdgf [kNPmax]; // Pdg code of k^th final state particle in hadronic system
377  double brEf [kNPmax]; // Energy of k^th final state particle in hadronic system @ LAB
378  double brPxf [kNPmax]; // Px of k^th final state particle in hadronic system @ LAB
379  double brPyf [kNPmax]; // Py of k^th final state particle in hadronic system @ LAB
380  double brPzf [kNPmax]; // Pz of k^th final state particle in hadronic system @ LAB
381  double brPf [kNPmax]; // P of k^th final state particle in hadronic system @ LAB
382  double brCosthf[kNPmax]; // cos(theta) of k^th final state particle in hadronic system @ LAB wrt to neutrino direction
383  int brNi = 0; // Nu. of particles in 'primary' hadronic system (before intranuclear rescattering)
384  int brPdgi[kNPmax]; // Pdg code of k^th particle in 'primary' hadronic system
385  int brResc[kNPmax]; // FSI code of k^th particle in 'primary' hadronic system
386  double brEi [kNPmax]; // Energy of k^th particle in 'primary' hadronic system @ LAB
387  double brPxi [kNPmax]; // Px of k^th particle in 'primary' hadronic system @ LAB
388  double brPyi [kNPmax]; // Py of k^th particle in 'primary' hadronic system @ LAB
389  double brPzi [kNPmax]; // Pz of k^th particle in 'primary' hadronic system @ LAB
390  double brVtxX; // Vertex x in detector coord system (SI)
391  double brVtxY; // Vertex y in detector coord system (SI)
392  double brVtxZ; // Vertex z in detector coord system (SI)
393  double brVtxT; // Vertex t in detector coord system (SI)
394  double brSumKEf; // Sum of kinetic energies of all final state particles
395  double brCalResp0; // Approximate calorimetric response to the hadronic system computed as sum of
396  // - (kinetic energy) for pi+, pi-, p, n
397  // - (energy + 2*mass) for antiproton, antineutron
398  // - ((e/h) * energy) for pi0, gamma, e-, e+, where e/h is set to 1.3
399  // - (kinetic energy) for other particles
400 
401  Double_t brXSec; // the event cross section in 1E-38cm^2
402  Double_t brDXSec; // is the differential cross section for the selected in 1E-38cm^2/{K^n}
403  UInt_t brKPS; // phase space that the xsec has been evaluated into
404 
405  // Open output file & create output summary tree & create the tree branches
406  //
407  LOG("gntpc", pNOTICE)
408  << "*** Saving summary tree to: " << gOptOutFileName;
409  TFile fout(gOptOutFileName.c_str(),"recreate");
410 
411  TTree * s_tree = new TTree("gst","GENIE Summary Event Tree");
412 
413  // Create tree branches
414  //
415  s_tree->Branch("iev", &brIev, "iev/I" );
416  s_tree->Branch("neu", &brNeutrino, "neu/I" );
417  s_tree->Branch("fspl", &brFSPrimLept, "fspl/I" );
418  s_tree->Branch("tgt", &brTarget, "tgt/I" );
419  s_tree->Branch("Z", &brTargetZ, "Z/I" );
420  s_tree->Branch("A", &brTargetA, "A/I" );
421  s_tree->Branch("hitnuc", &brHitNuc, "hitnuc/I" );
422  s_tree->Branch("hitqrk", &brHitQrk, "hitqrk/I" );
423  s_tree->Branch("resid", &brResId, "resid/I" );
424  s_tree->Branch("sea", &brFromSea, "sea/O" );
425  s_tree->Branch("qel", &brIsQel, "qel/O" );
426  s_tree->Branch("mec", &brIsMec, "mec/O" );
427  s_tree->Branch("res", &brIsRes, "res/O" );
428  s_tree->Branch("dis", &brIsDis, "dis/O" );
429  s_tree->Branch("coh", &brIsCoh, "coh/O" );
430  s_tree->Branch("dfr", &brIsDfr, "dfr/O" );
431  s_tree->Branch("imd", &brIsImd, "imd/O" );
432  s_tree->Branch("imdanh", &brIsImdAnh, "imdanh/O" );
433  s_tree->Branch("singlek", &brIsSingleK, "singlek/O" );
434  s_tree->Branch("nuel", &brIsNuEL, "nuel/O" );
435  s_tree->Branch("em", &brIsEM, "em/O" );
436  s_tree->Branch("cc", &brIsCC, "cc/O" );
437  s_tree->Branch("nc", &brIsNC, "nc/O" );
438  s_tree->Branch("charm", &brIsCharmPro, "charm/O" );
439  s_tree->Branch("amnugamma", &brIsAMNuGamma, "amnugamma/O" );
440  s_tree->Branch("neut_code", &brCodeNeut, "neut_code/I" );
441  s_tree->Branch("nuance_code", &brCodeNuance, "nuance_code/I" );
442  s_tree->Branch("wght", &brWeight, "wght/D" );
443  s_tree->Branch("xs", &brKineXs, "xs/D" );
444  s_tree->Branch("ys", &brKineYs, "ys/D" );
445  s_tree->Branch("ts", &brKineTs, "ts/D" );
446  s_tree->Branch("Q2s", &brKineQ2s, "Q2s/D" );
447  s_tree->Branch("Ws", &brKineWs, "Ws/D" );
448  s_tree->Branch("x", &brKineX, "x/D" );
449  s_tree->Branch("y", &brKineY, "y/D" );
450  s_tree->Branch("t", &brKineT, "t/D" );
451  s_tree->Branch("Q2", &brKineQ2, "Q2/D" );
452  s_tree->Branch("W", &brKineW, "W/D" );
453  s_tree->Branch("EvRF", &brEvRF, "EvRF/D" );
454  s_tree->Branch("Ev", &brEv, "Ev/D" );
455  s_tree->Branch("pxv", &brPxv, "pxv/D" );
456  s_tree->Branch("pyv", &brPyv, "pyv/D" );
457  s_tree->Branch("pzv", &brPzv, "pzv/D" );
458  s_tree->Branch("En", &brEn, "En/D" );
459  s_tree->Branch("pxn", &brPxn, "pxn/D" );
460  s_tree->Branch("pyn", &brPyn, "pyn/D" );
461  s_tree->Branch("pzn", &brPzn, "pzn/D" );
462  s_tree->Branch("El", &brEl, "El/D" );
463  s_tree->Branch("pxl", &brPxl, "pxl/D" );
464  s_tree->Branch("pyl", &brPyl, "pyl/D" );
465  s_tree->Branch("pzl", &brPzl, "pzl/D" );
466  s_tree->Branch("pl", &brPl, "pl/D" );
467  s_tree->Branch("cthl", &brCosthl, "cthl/D" );
468  s_tree->Branch("nfp", &brNfP, "nfp/I" );
469  s_tree->Branch("nfn", &brNfN, "nfn/I" );
470  s_tree->Branch("nfpip", &brNfPip, "nfpip/I" );
471  s_tree->Branch("nfpim", &brNfPim, "nfpim/I" );
472  s_tree->Branch("nfpi0", &brNfPi0, "nfpi0/I" );
473  s_tree->Branch("nfkp", &brNfKp, "nfkp/I" );
474  s_tree->Branch("nfkm", &brNfKm, "nfkm/I" );
475  s_tree->Branch("nfk0", &brNfK0, "nfk0/I" );
476  s_tree->Branch("nfem", &brNfEM, "nfem/I" );
477  s_tree->Branch("nfother", &brNfOther, "nfother/I" );
478  s_tree->Branch("nip", &brNiP, "nip/I" );
479  s_tree->Branch("nin", &brNiN, "nin/I" );
480  s_tree->Branch("nipip", &brNiPip, "nipip/I" );
481  s_tree->Branch("nipim", &brNiPim, "nipim/I" );
482  s_tree->Branch("nipi0", &brNiPi0, "nipi0/I" );
483  s_tree->Branch("nikp", &brNiKp, "nikp/I" );
484  s_tree->Branch("nikm", &brNiKm, "nikm/I" );
485  s_tree->Branch("nik0", &brNiK0, "nik0/I" );
486  s_tree->Branch("niem", &brNiEM, "niem/I" );
487  s_tree->Branch("niother", &brNiOther, "niother/I" );
488  s_tree->Branch("ni", &brNi, "ni/I" );
489  s_tree->Branch("pdgi", brPdgi, "pdgi[ni]/I" );
490  s_tree->Branch("resc", brResc, "resc[ni]/I" );
491  s_tree->Branch("Ei", brEi, "Ei[ni]/D" );
492  s_tree->Branch("pxi", brPxi, "pxi[ni]/D" );
493  s_tree->Branch("pyi", brPyi, "pyi[ni]/D" );
494  s_tree->Branch("pzi", brPzi, "pzi[ni]/D" );
495  s_tree->Branch("nf", &brNf, "nf/I" );
496  s_tree->Branch("pdgf", brPdgf, "pdgf[nf]/I" );
497  s_tree->Branch("Ef", brEf, "Ef[nf]/D" );
498  s_tree->Branch("pxf", brPxf, "pxf[nf]/D" );
499  s_tree->Branch("pyf", brPyf, "pyf[nf]/D" );
500  s_tree->Branch("pzf", brPzf, "pzf[nf]/D" );
501  s_tree->Branch("pf", brPf, "pf[nf]/D" );
502  s_tree->Branch("cthf", brCosthf, "cthf[nf]/D" );
503  s_tree->Branch("vtxx", &brVtxX, "vtxx/D" );
504  s_tree->Branch("vtxy", &brVtxY, "vtxy/D" );
505  s_tree->Branch("vtxz", &brVtxZ, "vtxz/D" );
506  s_tree->Branch("vtxt", &brVtxT, "vtxt/D" );
507  s_tree->Branch("sumKEf", &brSumKEf, "sumKEf/D" );
508  s_tree->Branch("calresp0", &brCalResp0, "calresp0/D" );
509  s_tree->Branch("XSec", &brXSec, "XSec/D" );
510  s_tree->Branch("DXSec", &brDXSec, "DXSec/D" );
511  s_tree->Branch("KPS", &brKPS, "KPS/i" );
512 
513  // Open the ROOT file and get the TTree & its header
514  TFile fin(gOptInpFileName.c_str(),"READ");
515  TTree * er_tree = 0;
516  NtpMCTreeHeader * thdr = 0;
517  er_tree = dynamic_cast <TTree *> ( fin.Get("gtree") );
518  thdr = dynamic_cast <NtpMCTreeHeader *> ( fin.Get("header") );
519  if (!er_tree) {
520  LOG("gntpc", pERROR) << "Null input GHEP event tree";
521  return;
522  }
523  LOG("gntpc", pINFO) << "Input tree header: " << *thdr;
524 
525  // Get the mc record
526  NtpMCEventRecord * mcrec = 0;
527  er_tree->SetBranchAddress("gmcrec", &mcrec);
528  if (!mcrec) {
529  LOG("gntpc", pERROR) << "Null MC record";
530  return;
531  }
532 
533  // Figure out how many events to analyze
534  Long64_t nmax = (gOptN<0) ?
535  er_tree->GetEntries() : TMath::Min( er_tree->GetEntries(), gOptN );
536  if (nmax<0) {
537  LOG("gntpc", pERROR) << "Number of events = 0";
538  return;
539  }
540 
541  LOG("gntpc", pNOTICE) << "*** Analyzing: " << nmax << " events";
542 
543  TLorentzVector pdummy(0,0,0,0);
544 
545  // Event loop
546  for(Long64_t iev = 0; iev < nmax; iev++) {
547  er_tree->GetEntry(iev);
548 
549  NtpMCRecHeader rec_header = mcrec->hdr;
550  EventRecord & event = *(mcrec->event);
551 
552  LOG("gntpc", pINFO) << rec_header;
553  LOG("gntpc", pINFO) << event;
554 
555  // Go further only if the event is physical
556  bool is_unphysical = event.IsUnphysical();
557  if(is_unphysical) {
558  LOG("gntpc", pINFO) << "Skipping unphysical event";
559  mcrec->Clear();
560  continue;
561  }
562 
563  // Clean-up arrays
564  //
565  for(int j=0; j<kNPmax; j++) {
566  brPdgi [j] = 0;
567  brResc [j] = -1;
568  brEi [j] = 0;
569  brPxi [j] = 0;
570  brPyi [j] = 0;
571  brPzi [j] = 0;
572  brPdgf [j] = 0;
573  brEf [j] = 0;
574  brPxf [j] = 0;
575  brPyf [j] = 0;
576  brPzf [j] = 0;
577  brPf [j] = 0;
578  brCosthf [j] = 0;
579  }
580 
581  // Computing event characteristics
582  //
583 
584  //input particles
585  GHepParticle * neutrino = event.Probe();
586  GHepParticle * target = event.Particle(1);
587  assert(target);
588  GHepParticle * fsl = event.FinalStatePrimaryLepton();
589  GHepParticle * hitnucl = event.HitNucleon();
590 
591  int tgtZ = 0;
592  int tgtA = 0;
593  if(pdg::IsIon(target->Pdg())) {
594  tgtZ = pdg::IonPdgCodeToZ(target->Pdg());
595  tgtA = pdg::IonPdgCodeToA(target->Pdg());
596  }
597  if(target->Pdg() == kPdgProton ) { tgtZ = 1; tgtA = 1; }
598  if(target->Pdg() == kPdgNeutron ) { tgtZ = 0; tgtA = 1; }
599 
600  // Summary info
601  const Interaction * interaction = event.Summary();
602  const InitialState & init_state = interaction->InitState();
603  const ProcessInfo & proc_info = interaction->ProcInfo();
604  const Kinematics & kine = interaction->Kine();
605  const XclsTag & xcls = interaction->ExclTag();
606  const Target & tgt = init_state.Tgt();
607 
608  // Vertex in detector coord system
609  TLorentzVector * vtx = event.Vertex();
610 
611  // Process id
612  bool is_qel = proc_info.IsQuasiElastic();
613  bool is_res = proc_info.IsResonant();
614  bool is_dis = proc_info.IsDeepInelastic();
615  bool is_coh = proc_info.IsCoherentProduction();
616  bool is_coh_el = proc_info.IsCoherentElastic();
617  bool is_dfr = proc_info.IsDiffractive();
618  bool is_imd = proc_info.IsInverseMuDecay();
619  bool is_imdanh = proc_info.IsIMDAnnihilation();
620  bool is_singlek = proc_info.IsSingleKaon();
621  bool is_nuel = proc_info.IsNuElectronElastic();
622  bool is_em = proc_info.IsEM();
623  bool is_weakcc = proc_info.IsWeakCC();
624  bool is_weaknc = proc_info.IsWeakNC();
625  bool is_mec = proc_info.IsMEC();
626  bool is_amnugamma = proc_info.IsAMNuGamma();
627 
628  if (!hitnucl && neutrino) {
629  assert(is_coh || is_imd || is_imdanh || is_nuel | is_amnugamma || is_coh_el);
630  }
631 
632  // Hit quark - set only for DIS events
633  int qrk = (is_dis) ? tgt.HitQrkPdg() : 0;
634  bool seaq = (is_dis) ? tgt.HitSeaQrk() : false;
635 
636  // Resonance id ($GENIE/src/BaryonResonance/BaryonResonance.h) -
637  // set only for resonance neutrinoproduction
638  int resid = (is_res) ? EResonance(xcls.Resonance()) : -99;
639 
640  // (qel or dis) charm production?
641  bool charm = xcls.IsCharmEvent();
642 
643  // Get NEUT and NUANCE equivalent reaction codes (if any)
644  brCodeNeut = utils::ghep::NeutReactionCode(&event);
645  brCodeNuance = utils::ghep::NuanceReactionCode(&event);
646 
647  // Get event weight
648  double weight = event.Weight();
649 
650  // Access kinematical params _exactly_ as they were selected internally
651  // (at the hit nucleon rest frame;
652  // for bound nucleons: taking into account fermi momentum and off-shell kinematics)
653  //
654  bool get_selected = true;
655  double xs = kine.x (get_selected);
656  double ys = kine.y (get_selected);
657  double ts = (is_coh || is_dfr) ? kine.t (get_selected) : -1;
658  double Q2s = kine.Q2(get_selected);
659  double Ws = kine.W (get_selected);
660 
661  LOG("gntpc", pDEBUG)
662  << "[Select] Q2 = " << Q2s << ", W = " << Ws
663  << ", x = " << xs << ", y = " << ys << ", t = " << ts;
664 
665  // Calculate the same kinematical params but now as an experimentalist would
666  // measure them by neglecting the fermi momentum and off-shellness of bound nucleons
667  //
668 
669  const TLorentzVector & k1 = (neutrino) ? *(neutrino->P4()) : pdummy; // v 4-p (k1)
670  const TLorentzVector & k2 = (fsl) ? *(fsl->P4()) : pdummy; // l 4-p (k2)
671  const TLorentzVector & p1 = (hitnucl) ? *(hitnucl->P4()) : pdummy; // N 4-p (p1)
672 
673  double M = kNucleonMass;
674  TLorentzVector q = k1-k2; // q=k1-k2, 4-p transfer
675  double Q2 = -1 * q.M2(); // momemtum transfer
676 
677  double v = (hitnucl) ? q.Energy() : -1; // v (E transfer to the nucleus)
678  double x, y, W2, W;
679  if(!is_coh){
680 
681  x = (hitnucl) ? 0.5*Q2/(M*v) : -1; // Bjorken x
682  y = (hitnucl) ? v/k1.Energy() : -1; // Inelasticity, y = q*P1/k1*P1
683 
684  W2 = (hitnucl) ? M*M + 2*M*v - Q2 : -1; // Hadronic Invariant mass ^ 2
685  W = (hitnucl) ? TMath::Sqrt(W2) : -1;
686  } else{
687 
688  v = q.Energy();
689  x = 0.5*Q2/(M*v); // Bjorken x
690  y = v/k1.Energy(); // Inelasticity, y = q*P1/k1*P1
691 
692  W2 = M*M + 2*M*v - Q2; // Hadronic Invariant mass ^ 2
693  W = TMath::Sqrt(W2);
694 
695  }
696 
697  double t = (is_coh || is_dfr) ? kine.t (get_selected) : -1;
698 
699  // Get v 4-p at hit nucleon rest-frame
700  TLorentzVector k1_rf = k1;
701  if(hitnucl) {
702  k1_rf.Boost(-1.*p1.BoostVector());
703  }
704 
705 // if(is_mec){
706 // v = q.Energy();
707 // x = 0.5*Q2/(M*v);
708 // y = v/k1.Energy();
709 // W2 = M*M + 2*M*v - Q2;
710 // W = TMath::Sqrt(W2);
711 // }
712 
713  LOG("gntpc", pDEBUG)
714  << "[Calc] Q2 = " << Q2 << ", W = " << W
715  << ", x = " << x << ", y = " << y << ", t = " << t;
716 
717  // Extract more info on the hadronic system
718  // Only for QEL/RES/DIS/COH/MEC events
719  //
720  bool study_hadsyst = (is_qel || is_res || is_dis || is_coh || is_dfr || is_mec || is_singlek);
721 
722  //
723  TObjArrayIter piter(&event);
724  GHepParticle * p = 0;
725  int ip=-1;
726 
727  //
728  // Extract the final state system originating from the hadronic vertex
729  // (after the intranuclear rescattering step)
730  //
731 
732  LOG("gntpc", pDEBUG) << "Extracting final state hadronic system";
733 
734  vector<int> final_had_syst;
735  while( (p = (GHepParticle *) piter.Next()) && study_hadsyst)
736  {
737  ip++;
738  // don't count final state lepton as part hadronic system
739  //if(!is_coh && event.Particle(ip)->FirstMother()==0) continue;
740  if(event.Particle(ip)->FirstMother()==0) continue;
741  if(pdg::IsPseudoParticle(p->Pdg())) continue;
742  int pdgc = p->Pdg();
743  int ist = p->Status();
744  if(ist==kIStStableFinalState) {
745  if (pdgc == kPdgGamma || pdgc == kPdgElectron || pdgc == kPdgPositron) {
746  int igmom = p->FirstMother();
747  if(igmom!=-1) {
748  // only count e+'s e-'s or gammas not from decay of pi0
749  if(event.Particle(igmom)->Pdg() != kPdgPi0) { final_had_syst.push_back(ip); }
750  }
751  } else {
752  final_had_syst.push_back(ip);
753  }
754  }
755  // now add pi0's that were decayed as short lived particles
756  else if(pdgc == kPdgPi0){
757  int ifd = p->FirstDaughter();
758  if ( ifd != -1 ) {
759  int fd_pdgc = event.Particle(ifd)->Pdg();
760  // just require that first daughter is one of gamma, e+ or e-
761  if(fd_pdgc == kPdgGamma || fd_pdgc == kPdgElectron || fd_pdgc == kPdgPositron){
762  final_had_syst.push_back(ip);
763  }
764  }
765  }
766  }//particle-loop
767 
768  if( count(final_had_syst.begin(), final_had_syst.end(), -1) > 0) {
769  mcrec->Clear();
770  continue;
771  }
772 
773  //
774  // Extract info on the primary hadronic system (before any intranuclear rescattering)
775  // looking for particles with status_code == kIStHadronInTheNucleus
776  // An exception is the coherent production and scattering off free nucleon targets
777  // (no intranuclear rescattering) in which case primary hadronic system is set to be
778  // 'identical' with the final state hadronic system
779  //
780 
781  LOG("gntpc", pDEBUG) << "Extracting primary hadronic system";
782 
783  ip = -1;
784  TObjArrayIter piter_prim(&event);
785 
786  vector<int> prim_had_syst;
787  if(study_hadsyst) {
788  // if coherent or free nucleon target set primary states equal to final states
789 
790  if(!pdg::IsIon(target->Pdg()) || (is_coh)) {
791 
792  for( vector<int>::const_iterator hiter = final_had_syst.begin();
793  hiter != final_had_syst.end(); ++hiter) {
794 
795  prim_had_syst.push_back(*hiter);
796  }
797  }
798 
799  else {
800 
801  // otherwise loop over all particles and store indices of those which are hadrons
802  // created within the nucleus
803  /* else {
804  while( (p = (GHepParticle *) piter_prim.Next()) ){
805  ip++;
806  int ist_comp = p->Status();
807  if(ist_comp==kIStHadronInTheNucleus) {
808  prim_had_syst.push_back(ip);
809  }
810  }//particle-loop */
811  //
812 
813 
814  //to find the true particles emitted from the principal vertex,
815  // looping over all Ist=14 particles ok for hA, but doesn't
816  // work for hN. We must now look specifically for these particles.
817  int ist_store = -10;
818  if(is_res){
819  while( (p = (GHepParticle *) piter_prim.Next()) ){
820  ip++;
821  int ist_comp = p->Status();
822  if(ist_comp==kIStDecayedState) {
823  ist_store = ip; //store this mother
824  continue;
825  }
826  // LOG("gntpc",pNOTICE) << p->FirstMother()<< " "<<ist_store;
827  if(p->FirstMother()==ist_store) {
828  prim_had_syst.push_back(ip);
829  }
830  }
831  }
832  if(is_dis){
833  while( (p = (GHepParticle *) piter_prim.Next()) ){
834  ip++;
835  int ist_comp = p->Status();
836  if(ist_comp==kIStDISPreFragmHadronicState) {
837  ist_store = ip; //store this mother
838  continue;
839  }
840  if(p->FirstMother()==ist_store) {
841  prim_had_syst.push_back(ip);
842  }
843  }
844  }
845  if(is_qel){
846  while( (p = (GHepParticle *) piter_prim.Next()) ){
847  ip++;
848  int ist_comp = p->Status();
849  if(ist_comp==kIStNucleonTarget) {
850  ist_store = ip; //store this mother
851  continue;
852  }
853  // LOG("gntpc",pNOTICE) << p->FirstMother()<< " "<<ist_store;
854  if(p->FirstMother()==ist_store) {
855  prim_had_syst.push_back(ip);
856  }
857  }
858  }
859  if(is_mec){
860  while( (p = (GHepParticle *) piter_prim.Next()) ){
861  ip++;
862  int ist_comp = p->Status();
863  if(ist_comp==kIStDecayedState) {
864  ist_store = ip; //store this mother
865  continue;
866  }
867  // LOG("gntpc",pNOTICE) << "MEC: " << p->FirstMother()<< " "<<ist_store;
868  if(p->FirstMother()==ist_store) {
869  prim_had_syst.push_back(ip);
870  }
871  }
872  }
873 
874 
875  // also include gammas from nuclear de-excitations (appearing in the daughter list of the
876  // hit nucleus, earlier than the primary hadronic system extracted above)
877  for(int i = target->FirstDaughter(); i <= target->LastDaughter(); i++) {
878  if(i<0) continue;
879  if(event.Particle(i)->Status()==kIStStableFinalState) { prim_had_syst.push_back(i); }
880  }
881 
882 
883  } // else from ( not ion or coherent )
884 
885  }//study_hadsystem?
886 
887  if( count(prim_had_syst.begin(), prim_had_syst.end(), -1) > 0) {
888  mcrec->Clear();
889  continue;
890  }
891 
892  //
893  // Al information has been assembled -- Start filling up the tree branches
894  //
895  brIev = (int) iev;
896  brNeutrino = (neutrino) ? neutrino->Pdg() : 0;
897  brFSPrimLept = (fsl) ? fsl->Pdg() : 0;
898  brTarget = target->Pdg();
899  brTargetZ = tgtZ;
900  brTargetA = tgtA;
901  brHitNuc = (hitnucl) ? hitnucl->Pdg() : 0;
902  brHitQrk = qrk;
903  brFromSea = seaq;
904  brResId = resid;
905  brIsQel = is_qel;
906  brIsRes = is_res;
907  brIsDis = is_dis;
908  brIsCoh = is_coh;
909  brIsDfr = is_dfr;
910  brIsImd = is_imd;
911  brIsSingleK = is_singlek;
912  brIsNuEL = is_nuel;
913  brIsEM = is_em;
914  brIsMec = is_mec;
915  brIsCC = is_weakcc;
916  brIsNC = is_weaknc;
917  brIsCharmPro = charm;
918  brIsAMNuGamma= is_amnugamma;
919  brWeight = weight;
920  brKineXs = xs;
921  brKineYs = ys;
922  brKineTs = ts;
923  brKineQ2s = Q2s;
924  brKineWs = Ws;
925  brKineX = x;
926  brKineY = y;
927  brKineT = t;
928  brKineQ2 = Q2;
929  brKineW = W;
930  brEvRF = k1_rf.Energy();
931  brEv = k1.Energy();
932  brPxv = k1.Px();
933  brPyv = k1.Py();
934  brPzv = k1.Pz();
935  brEn = (hitnucl) ? p1.Energy() : 0;
936  brPxn = (hitnucl) ? p1.Px() : 0;
937  brPyn = (hitnucl) ? p1.Py() : 0;
938  brPzn = (hitnucl) ? p1.Pz() : 0;
939  brEl = k2.Energy();
940  brPxl = k2.Px();
941  brPyl = k2.Py();
942  brPzl = k2.Pz();
943  brPl = k2.P();
944  brCosthl = TMath::Cos( k2.Vect().Angle(k1.Vect()) );
945 
946  // XSec Info
947 
948  brXSec = event.XSec()*(1E+38/units::cm2);
949  brDXSec = event.DiffXSec()*(1E+38/units::cm2);
950  brKPS = event.DiffXSecVars();
951 
952  // Primary hadronic system (from primary neutrino interaction, before FSI)
953  brNiP = 0;
954  brNiN = 0;
955  brNiPip = 0;
956  brNiPim = 0;
957  brNiPi0 = 0;
958  brNiKp = 0;
959  brNiKm = 0;
960  brNiK0 = 0;
961  brNiEM = 0;
962  brNiOther = 0;
963  brNi = prim_had_syst.size();
964  for(int j=0; j<brNi; j++) {
965  p = event.Particle(prim_had_syst[j]);
966  assert(p);
967  brPdgi[j] = p->Pdg();
968  brResc[j] = p->RescatterCode();
969  brEi [j] = p->Energy();
970  brPxi [j] = p->Px();
971  brPyi [j] = p->Py();
972  brPzi [j] = p->Pz();
973 
974  if (p->Pdg() == kPdgProton || p->Pdg() == kPdgAntiProton) brNiP++;
975  else if (p->Pdg() == kPdgNeutron || p->Pdg() == kPdgAntiNeutron) brNiN++;
976  else if (p->Pdg() == kPdgPiP) brNiPip++;
977  else if (p->Pdg() == kPdgPiM) brNiPim++;
978  else if (p->Pdg() == kPdgPi0) brNiPi0++;
979  else if (p->Pdg() == kPdgKP) brNiKp++;
980  else if (p->Pdg() == kPdgKM) brNiKm++;
981  else if (p->Pdg() == kPdgK0 || p->Pdg() == kPdgAntiK0) brNiK0++;
982  else if (p->Pdg() == kPdgGamma || p->Pdg() == kPdgElectron || p->Pdg() == kPdgPositron) brNiEM++;
983  else brNiOther++;
984 
985  LOG("gntpc", pINFO)
986  << "Counting in primary hadronic system: idx = " << prim_had_syst[j]
987  << " -> " << p->Name();
988  }
989 
990  LOG("gntpc", pINFO)
991  << "N(p):" << brNiP
992  << ", N(n):" << brNiN
993  << ", N(pi+):" << brNiPip
994  << ", N(pi-):" << brNiPim
995  << ", N(pi0):" << brNiPi0
996  << ", N(K+,K-,K0):" << brNiKp+brNiKm+brNiK0
997  << ", N(gamma,e-,e+):" << brNiEM
998  << ", N(etc):" << brNiOther << "\n";
999 
1000  // Final state (visible) hadronic system
1001  brNfP = 0;
1002  brNfN = 0;
1003  brNfPip = 0;
1004  brNfPim = 0;
1005  brNfPi0 = 0;
1006  brNfKp = 0;
1007  brNfKm = 0;
1008  brNfK0 = 0;
1009  brNfEM = 0;
1010  brNfOther = 0;
1011 
1012  brSumKEf = (fsl) ? fsl->KinE() : 0;
1013  brCalResp0 = 0;
1014 
1015  brNf = final_had_syst.size();
1016  for(int j=0; j<brNf; j++) {
1017  p = event.Particle(final_had_syst[j]);
1018  assert(p);
1019 
1020  int hpdg = p->Pdg();
1021  double hE = p->Energy();
1022  double hKE = p->KinE();
1023  double hpx = p->Px();
1024  double hpy = p->Py();
1025  double hpz = p->Pz();
1026  double hp = TMath::Sqrt(hpx*hpx + hpy*hpy + hpz*hpz);
1027  double hm = p->Mass();
1028  double hcth = TMath::Cos( p->P4()->Vect().Angle(k1.Vect()) );
1029 
1030  brPdgf [j] = hpdg;
1031  brEf [j] = hE;
1032  brPxf [j] = hpx;
1033  brPyf [j] = hpy;
1034  brPzf [j] = hpz;
1035  brPf [j] = hp;
1036  brCosthf[j] = hcth;
1037 
1038  brSumKEf += hKE;
1039 
1040  if ( hpdg == kPdgProton ) { brNfP++; brCalResp0 += hKE; }
1041  else if ( hpdg == kPdgAntiProton ) { brNfP++; brCalResp0 += (hE + 2*hm);}
1042  else if ( hpdg == kPdgNeutron ) { brNfN++; brCalResp0 += hKE; }
1043  else if ( hpdg == kPdgAntiNeutron ) { brNfN++; brCalResp0 += (hE + 2*hm);}
1044  else if ( hpdg == kPdgPiP ) { brNfPip++; brCalResp0 += hKE; }
1045  else if ( hpdg == kPdgPiM ) { brNfPim++; brCalResp0 += hKE; }
1046  else if ( hpdg == kPdgPi0 ) { brNfPi0++; brCalResp0 += (e_h * hE); }
1047  else if ( hpdg == kPdgKP ) { brNfKp++; brCalResp0 += hKE; }
1048  else if ( hpdg == kPdgKM ) { brNfKm++; brCalResp0 += hKE; }
1049  else if ( hpdg == kPdgK0 ) { brNfK0++; brCalResp0 += hKE; }
1050  else if ( hpdg == kPdgAntiK0 ) { brNfK0++; brCalResp0 += hKE; }
1051  else if ( hpdg == kPdgGamma ) { brNfEM++; brCalResp0 += (e_h * hE); }
1052  else if ( hpdg == kPdgElectron ) { brNfEM++; brCalResp0 += (e_h * hE); }
1053  else if ( hpdg == kPdgPositron ) { brNfEM++; brCalResp0 += (e_h * hE); }
1054  else { brNfOther++; brCalResp0 += hKE; }
1055 
1056  LOG("gntpc", pINFO)
1057  << "Counting in f/s system from hadronic vtx: idx = " << final_had_syst[j]
1058  << " -> " << p->Name();
1059  }
1060 
1061  LOG("gntpc", pINFO)
1062  << "N(p):" << brNfP
1063  << ", N(n):" << brNfN
1064  << ", N(pi+):" << brNfPip
1065  << ", N(pi-):" << brNfPim
1066  << ", N(pi0):" << brNfPi0
1067  << ", N(K+,K-,K0):" << brNfKp+brNfKm+brNfK0
1068  << ", N(gamma,e-,e+):" << brNfEM
1069  << ", N(etc):" << brNfOther << "\n";
1070 
1071  brVtxX = vtx->X();
1072  brVtxY = vtx->Y();
1073  brVtxZ = vtx->Z();
1074  brVtxT = vtx->T();
1075 
1076  s_tree->Fill();
1077 
1078  mcrec->Clear();
1079 
1080  } // event loop
1081 
1082 
1083  // Copy MC job metadata (gconfig and genv TFolders)
1084  if(gOptCopyJobMeta) {
1085  TFolder * genv = (TFolder*) fin.Get("genv");
1086  TFolder * gconfig = (TFolder*) fin.Get("gconfig");
1087  fout.cd();
1088  genv -> Write("genv");
1089  gconfig -> Write("gconfig");
1090  }
1091 
1092  fin.Close();
1093 
1094  fout.Write();
1095  fout.Close();
1096 }
1097 //____________________________________________________________________________________
1098 // GENIE GHEP EVENT TREE FORMAT -> GENIE XML EVENT FILE FORMAT
1099 //____________________________________________________________________________________
1100 void ConvertToGXML(void)
1101 {
1102  //-- open the ROOT file and get the TTree & its header
1103  TFile fin(gOptInpFileName.c_str(),"READ");
1104  TTree * tree = 0;
1105  NtpMCTreeHeader * thdr = 0;
1106  tree = dynamic_cast <TTree *> ( fin.Get("gtree") );
1107  thdr = dynamic_cast <NtpMCTreeHeader *> ( fin.Get("header") );
1108 
1109  LOG("gntpc", pINFO) << "Input tree header: " << *thdr;
1110 
1111  //-- get mc record
1112  NtpMCEventRecord * mcrec = 0;
1113  tree->SetBranchAddress("gmcrec", &mcrec);
1114 
1115  //-- open the output stream
1116  ofstream output(gOptOutFileName.c_str(), ios::out);
1117 
1118  //-- add required header
1119  output << "<?xml version=\"1.0\" encoding=\"ISO-8859-1\"?>";
1120  output << endl << endl;
1121  output << "<!-- generated by GENIE gntpc utility -->";
1122  output << endl << endl;
1123  output << "<genie_event_list version=\"1.00\">" << endl;
1124 
1125  //-- figure out how many events to analyze
1126  Long64_t nmax = (gOptN<0) ?
1127  tree->GetEntries() : TMath::Min(tree->GetEntries(), gOptN);
1128  if (nmax<0) {
1129  LOG("gntpc", pERROR) << "Number of events = 0";
1130  return;
1131  }
1132  LOG("gntpc", pNOTICE) << "*** Analyzing: " << nmax << " events";
1133 
1134  //-- event loop
1135  for(Long64_t iev = 0; iev < nmax; iev++) {
1136  tree->GetEntry(iev);
1137  NtpMCRecHeader rec_header = mcrec->hdr;
1138  EventRecord & event = *(mcrec->event);
1139 
1140  LOG("gntpc", pINFO) << rec_header;
1141  LOG("gntpc", pINFO) << event;
1142 
1143  //
1144  // convert the current event
1145  //
1146 
1147  output << endl << endl;
1148  output << " <!-- GENIE GHEP event -->" << endl;
1149  output << " <ghep np=\"" << event.GetEntries()
1150  << "\" unphysical=\""
1151  << (event.IsUnphysical() ? "true" : "false") << "\">" << endl;
1152  output << setiosflags(ios::scientific);
1153 
1154  // write-out the event-wide properties
1155  output << " ";
1156  output << " <!-- event weight -->";
1157  output << " <wgt> " << event.Weight() << " </wgt>";
1158  output << endl;
1159  output << " ";
1160  output << " <!-- cross sections -->";
1161  output << " <xsec_evnt> " << event.XSec() << " </xsec_evnt>";
1162  output << " <xsec_kine> " << event.DiffXSec() << " </xsec_kine>";
1163  output << endl;
1164  output << " ";
1165  output << " <!-- event vertex -->";
1166  output << " <vx> " << event.Vertex()->X() << " </vx>";
1167  output << " <vy> " << event.Vertex()->Y() << " </vy>";
1168  output << " <vz> " << event.Vertex()->Z() << " </vz>";
1169  output << " <vt> " << event.Vertex()->T() << " </vt>";
1170  output << endl;
1171 
1172  // write-out the generated particle list
1173  output << " <!-- particle list -->" << endl;
1174  unsigned int i=0;
1175  GHepParticle * p = 0;
1176  TIter event_iter(&event);
1177  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
1178  string type = "U";
1179  if (pdg::IsPseudoParticle(p->Pdg())) type = "F";
1180  else if (pdg::IsParticle (p->Pdg())) type = "P";
1181  else if (pdg::IsIon (p->Pdg())) type = "N";
1182 
1183  output << " <p idx=\"" << i << "\" type=\"" << type << "\">" << endl;
1184  output << " ";
1185  output << " <pdg> " << p->Pdg() << " </pdg>";
1186  output << " <ist> " << p->Status() << " </ist>";
1187  output << endl;
1188  output << " ";
1189  output << " <mother> "
1190  << " <fst> " << setfill(' ') << setw(3) << p->FirstMother() << " </fst> "
1191  << " <lst> " << setfill(' ') << setw(3) << p->LastMother() << " </lst> "
1192  << " </mother>";
1193  output << endl;
1194  output << " ";
1195  output << " <daughter> "
1196  << " <fst> " << setfill(' ') << setw(3) << p->FirstDaughter() << " </fst> "
1197  << " <lst> " << setfill(' ') << setw(3) << p->LastDaughter() << " </lst> "
1198  << " </daughter>";
1199  output << endl;
1200  output << " ";
1201  output << " <px> " << setfill(' ') << setw(20) << p->Px() << " </px>";
1202  output << " <py> " << setfill(' ') << setw(20) << p->Py() << " </py>";
1203  output << " <pz> " << setfill(' ') << setw(20) << p->Pz() << " </pz>";
1204  output << " <E> " << setfill(' ') << setw(20) << p->E() << " </E> ";
1205  output << endl;
1206  output << " ";
1207  output << " <x> " << setfill(' ') << setw(20) << p->Vx() << " </x> ";
1208  output << " <y> " << setfill(' ') << setw(20) << p->Vy() << " </y> ";
1209  output << " <z> " << setfill(' ') << setw(20) << p->Vz() << " </z> ";
1210  output << " <t> " << setfill(' ') << setw(20) << p->Vt() << " </t> ";
1211  output << endl;
1212 
1213  if(p->PolzIsSet()) {
1214  output << " ";
1215  output << " <ppolar> " << p->PolzPolarAngle() << " </ppolar>";
1216  output << " <pazmth> " << p->PolzAzimuthAngle() << " </pazmth>";
1217  output << endl;
1218  }
1219 
1220  if(p->RescatterCode() != -1) {
1221  output << " ";
1222  output << " <rescatter> " << p->RescatterCode() << " </rescatter>";
1223  output << endl;
1224  }
1225 
1226  output << " </p>" << endl;
1227  i++;
1228  }
1229  output << " </ghep>" << endl;
1230 
1231  mcrec->Clear();
1232  } // event loop
1233 
1234  //-- add required footer
1235  output << endl << endl;
1236  output << "<genie_event_list version=\"1.00\">";
1237 
1238  output.close();
1239  fin.Close();
1240 
1241  LOG("gntpc", pINFO) << "\nDone converting GENIE's GHEP ntuple";
1242 }
1243 //____________________________________________________________________________________
1244 // GENIE GHEP FORMAT -> GHEP MOCK DATA FORMAT
1245 //____________________________________________________________________________________
1246 void ConvertToGHepMock(void)
1247 {
1248  //-- open the ROOT file and get the TTree & its header
1249  TFile fin(gOptInpFileName.c_str(),"READ");
1250  TTree * tree = 0;
1251  NtpMCTreeHeader * thdr = 0;
1252  tree = dynamic_cast <TTree *> ( fin.Get("gtree") );
1253  thdr = dynamic_cast <NtpMCTreeHeader *> ( fin.Get("header") );
1254 
1255  LOG("gntpc", pINFO) << "Input tree header: " << *thdr;
1256 
1257  //-- get mc record
1258  NtpMCEventRecord * mcrec = 0;
1259  tree->SetBranchAddress("gmcrec", &mcrec);
1260 
1261  //-- figure out how many events to analyze
1262  Long64_t nmax = (gOptN<0) ?
1263  tree->GetEntries() : TMath::Min(tree->GetEntries(), gOptN);
1264  if (nmax<0) {
1265  LOG("gntpc", pERROR) << "Number of events = 0";
1266  return;
1267  }
1268  LOG("gntpc", pNOTICE) << "*** Analyzing: " << nmax << " events";
1269 
1270  //-- initialize an Ntuple Writer
1271  NtpWriter ntpw(kNFGHEP, thdr->runnu);
1272  ntpw.CustomizeFilename(gOptOutFileName);
1273  ntpw.Initialize();
1274 
1275  //-- event loop
1276  for(Long64_t iev = 0; iev < nmax; iev++) {
1277  tree->GetEntry(iev);
1278  NtpMCRecHeader rec_header = mcrec->hdr;
1279  EventRecord & event = *(mcrec->event);
1280 
1281  LOG("gntpc", pINFO) << rec_header;
1282  LOG("gntpc", pINFO) << event;
1283 
1284  EventRecord * stripped_event = new EventRecord;
1285  Interaction * nullint = new Interaction;
1286 
1287  stripped_event -> AttachSummary (nullint);
1288  stripped_event -> SetWeight (event.Weight());
1289  stripped_event -> SetVertex (*event.Vertex());
1290 
1291  GHepParticle * p = 0;
1292  TIter iter(&event);
1293  while( (p = (GHepParticle *)iter.Next()) ) {
1294  if(!p) continue;
1295  GHepStatus_t ist = p->Status();
1296  if(ist!=kIStStableFinalState) continue;
1297  stripped_event->AddParticle(
1298  p->Pdg(), ist, -1,-1,-1,-1, *p->P4(), *p->X4());
1299  }//p
1300 
1301  ntpw.AddEventRecord(iev,stripped_event);
1302 
1303  mcrec->Clear();
1304  } // event loop
1305 
1306  //-- save the generated MC events
1307  ntpw.Save();
1308 
1309  //-- rename the output file
1310 
1311  fin.Close();
1312 
1313  LOG("gntpc", pINFO) << "\nDone converting GENIE's GHEP ntuple";
1314 }
1315 //____________________________________________________________________________________
1316 // GENIE GHEP EVENT TREE FORMAT -> TRACKER FORMATS
1317 //____________________________________________________________________________________
1318 void ConvertToGTracker(void)
1319 {
1320  //-- open the ROOT file and get the TTree & its header
1321  TFile fin(gOptInpFileName.c_str(),"READ");
1322  TTree * tree = 0;
1323  NtpMCTreeHeader * thdr = 0;
1324  tree = dynamic_cast <TTree *> ( fin.Get("gtree") );
1325  thdr = dynamic_cast <NtpMCTreeHeader *> ( fin.Get("header") );
1326 
1327  LOG("gntpc", pINFO) << "Input tree header: " << *thdr;
1328 
1329  gFileMajorVrs = utils::system::GenieMajorVrsNum(thdr->cvstag.GetString().Data());
1330  gFileMinorVrs = utils::system::GenieMinorVrsNum(thdr->cvstag.GetString().Data());
1331  gFileRevisVrs = utils::system::GenieRevisVrsNum(thdr->cvstag.GetString().Data());
1332 
1333  //-- get mc record
1334  NtpMCEventRecord * mcrec = 0;
1335  tree->SetBranchAddress("gmcrec", &mcrec);
1336 
1337 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
1338  flux::GJPARCNuFluxPassThroughInfo * flux_info = 0;
1339  tree->SetBranchAddress("flux", &flux_info);
1340 #else
1341  LOG("gntpc", pWARN)
1342  << "\n Flux drivers are not enabled."
1343  << "\n No flux pass-through information will be written-out in the rootracker file"
1344  << "\n If this isn't what you are supposed to be doing then build GENIE by adding "
1345  << "--with-flux-drivers in the configuration step.";
1346 #endif
1347 
1348  //-- open the output stream
1349  ofstream output(gOptOutFileName.c_str(), ios::out);
1350 
1351  //-- figure out how many events to analyze
1352  Long64_t nmax = (gOptN<0) ?
1353  tree->GetEntries() : TMath::Min(tree->GetEntries(), gOptN);
1354  if (nmax<0) {
1355  LOG("gntpc", pERROR) << "Number of events = 0";
1356  return;
1357  }
1358  LOG("gntpc", pNOTICE) << "*** Analyzing: " << nmax << " events";
1359 
1360  //-- event loop
1361  for(Long64_t iev = 0; iev < nmax; iev++) {
1362  tree->GetEntry(iev);
1363  NtpMCRecHeader rec_header = mcrec->hdr;
1364  EventRecord & event = *(mcrec->event);
1365  Interaction * interaction = event.Summary();
1366 
1367  LOG("gntpc", pINFO) << rec_header;
1368  LOG("gntpc", pINFO) << event;
1369 
1370  GHepParticle * p = 0;
1371  TIter event_iter(&event);
1372  int iparticle = -1;
1373 
1374  // **** Convert the current event:
1375 
1376  //
1377  // -- Add tracker begin tag
1378  //
1379  output << "$ begin" << endl;
1380 
1381  //
1382  // -- Add the appropriate reaction code
1383  //
1384 
1385  // add 'NEUT'-like event type
1386  if(gOptOutFileFormat == kConvFmt_t2k_tracker) {
1387  int evtype = utils::ghep::NeutReactionCode(&event);
1388  LOG("gntpc", pNOTICE) << "NEUT-like event type = " << evtype;
1389  output << "$ genie " << evtype << endl;
1390  } //neut code
1391 
1392  // add 'NUANCE'-like event type
1393  else if(gOptOutFileFormat == kConvFmt_nuance_tracker) {
1394  int evtype = utils::ghep::NuanceReactionCode(&event);
1395  LOG("gntpc", pNOTICE) << "NUANCE-like event type = " << evtype;
1396  output << "$ nuance " << evtype << endl;
1397  } // nuance code
1398 
1399  else {
1400  gAbortingInErr = true;
1401  exit(1);
1402  }
1403 
1404  //
1405  // -- Add '$vertex' line
1406  //
1407  output << "$ vertex "
1408  << event.Vertex()->X() << " "
1409  << event.Vertex()->Y() << " "
1410  << event.Vertex()->Z() << " "
1411  << event.Vertex()->T() << endl;
1412 
1413  //
1414  // -- Add '$track' lines
1415  //
1416 
1417  // Loop over the generated GHEP particles and decide which ones
1418  // to write-out in $track lines
1419  vector<int> tracks;
1420 
1421  event_iter.Reset();
1422  iparticle = -1;
1423  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) )
1424  {
1425  iparticle++;
1426 
1427  int ghep_pdgc = p->Pdg();
1428  GHepStatus_t ghep_ist = (GHepStatus_t) p->Status();
1429 
1430  // Neglect all GENIE pseudo-particles
1431  if(pdg::IsPseudoParticle(ghep_pdgc)) continue;
1432 
1433  //
1434  // Keep 'initial state', 'nucleon target', 'hadron in the nucleus' and 'final state' particles.
1435  // Neglect pi0 decays if they were performed within GENIE (write out the decayed pi0 and neglect
1436  // the {gamma + gamma} or {gamma + e- + e+} final state
1437  //
1438 
1439  // is pi0 decay?
1440  bool is_pi0_dec = false;
1441  if(ghep_ist == kIStDecayedState && ghep_pdgc == kPdgPi0) {
1442  vector<int> pi0dv; // daughters vector
1443  int ghep_fd = p->FirstDaughter();
1444  int ghep_ld = p->LastDaughter();
1445  for(int jd = ghep_fd; jd <= ghep_ld; jd++) {
1446  if(jd!=-1) {
1447  pi0dv.push_back(event.Particle(jd)->Pdg());
1448  }
1449  }
1450  sort(pi0dv.begin(), pi0dv.end());
1451  is_pi0_dec = (pi0dv.size()==2 && pi0dv[0]==kPdgGamma && pi0dv[1]==kPdgGamma) ||
1452  (pi0dv.size()==3 && pi0dv[0]==kPdgPositron && pi0dv[1]==kPdgElectron && pi0dv[2]==kPdgGamma);
1453  }
1454 
1455  // is pi0 decay product?
1456  int ghep_fm = p->FirstMother();
1457  int ghep_fmpdgc = (ghep_fm==-1) ? 0 : event.Particle(ghep_fm)->Pdg();
1458  bool is_pi0_dpro = (ghep_pdgc == kPdgGamma && ghep_fmpdgc == kPdgPi0) ||
1459  (ghep_pdgc == kPdgElectron && ghep_fmpdgc == kPdgPi0) ||
1460  (ghep_pdgc == kPdgPositron && ghep_fmpdgc == kPdgPi0);
1461 
1462  bool keep = (ghep_ist == kIStInitialState) ||
1463  (ghep_ist == kIStNucleonTarget) ||
1464  (ghep_ist == kIStHadronInTheNucleus) ||
1465  (ghep_ist == kIStDecayedState && is_pi0_dec ) ||
1466  (ghep_ist == kIStStableFinalState && !is_pi0_dpro);
1467  if(!keep) continue;
1468 
1469  // Apparently SKDETSIM chokes with O16 - Neglect the nuclear target in this case
1470  //
1471  if (gOptOutFileFormat == kConvFmt_t2k_tracker && pdg::IsIon(p->Pdg())) continue;
1472 
1473  tracks.push_back(iparticle);
1474  }
1475 
1476  //bool info_added = false;
1477 
1478  // Looping twice to ensure that all final state particle are grouped together.
1479  // On the second loop add only f/s particles. On the first loop add all but f/s particles
1480  for(int iloop=0; iloop<=1; iloop++)
1481  {
1482  for(vector<int>::const_iterator ip = tracks.begin(); ip != tracks.end(); ++ip)
1483  {
1484  iparticle = *ip;
1485  p = event.Particle(iparticle);
1486 
1487  int ghep_pdgc = p->Pdg();
1488  GHepStatus_t ghep_ist = (GHepStatus_t) p->Status();
1489 
1490  bool fs = (ghep_ist==kIStStableFinalState) ||
1491  (ghep_ist==kIStDecayedState && ghep_pdgc==kPdgPi0);
1492 
1493  if(iloop==0 && fs) continue;
1494  if(iloop==1 && !fs) continue;
1495 
1496  // Convert GENIE's GHEP pdgc & status to NUANCE's equivalent
1497  //
1498  int ist;
1499  switch (ghep_ist) {
1500  case kIStInitialState: ist = -1; break;
1501  case kIStStableFinalState: ist = 0; break;
1502  case kIStIntermediateState: ist = -2; break;
1503  case kIStDecayedState: ist = (ghep_pdgc==kPdgPi0) ? 0 : -2; break;
1504  case kIStNucleonTarget: ist = -1; break;
1505  case kIStDISPreFragmHadronicState: ist = -999; break;
1506  case kIStPreDecayResonantState: ist = -999; break;
1507  case kIStHadronInTheNucleus: ist = -2; break;
1508  case kIStUndefined: ist = -999; break;
1509  default: ist = -999; break;
1510  }
1511  // Convert GENIE pdg code -> nuance PDG code
1512  // For most particles both generators use the standard PDG codes.
1513  // For nuclei GENIE follows the PDG-convention: 10LZZZAAAI
1514  // NUANCE is using: ZZZAAA
1515  int pdgc = ghep_pdgc;
1516  if ( pdg::IsIon(p->Pdg()) ) {
1517  int Z = pdg::IonPdgCodeToZ(ghep_pdgc);
1518  int A = pdg::IonPdgCodeToA(ghep_pdgc);
1519  pdgc = 1000*Z + A;
1520  }
1521 
1522  // The SK detector MC expects K0_Long, K0_Short - not K0, \bar{K0}
1523  // Do the conversion here:
1524  if(gOptOutFileFormat == kConvFmt_t2k_tracker) {
1525  if(pdgc==kPdgK0 || pdgc==kPdgAntiK0) {
1526  RandomGen * rnd = RandomGen::Instance();
1527  double R = rnd->RndGen().Rndm();
1528  if(R>0.5) pdgc = kPdgK0L;
1529  else pdgc = kPdgK0S;
1530  }
1531  }
1532  // Get particle's energy & momentum
1533  const TLorentzVector * p4 = p->P4();
1534  double E = p4->Energy() / units::MeV;
1535  double Px = p4->Px() / units::MeV;
1536  double Py = p4->Py() / units::MeV;
1537  double Pz = p4->Pz() / units::MeV;
1538  double P = p4->P() / units::MeV;
1539  // Compute direction cosines
1540  double dcosx = (P>0) ? Px/P : -999;
1541  double dcosy = (P>0) ? Py/P : -999;
1542  double dcosz = (P>0) ? Pz/P : -999;
1543 
1544 // <obsolte/>
1545 // GHepStatus_t gist = (GHepStatus_t) p->Status();
1546 // bool is_init =
1547 // (gist == kIStInitialState || gist == kIStNucleonTarget);
1548 //
1549 // if(!is_init && !info_added) {
1550 // // Add nuance obsolete and flux info (not filled in by
1551 // // GENIE here). Add it once after the initial state particles
1552 // output << "$ info 2 949000 0.0000E+00" << endl;
1553 // info_added = true;
1554 // }
1555 // </obsolte>
1556 
1557  LOG("gntpc", pNOTICE)
1558  << "Adding $track corrsponding to GHEP particle at position: " << iparticle
1559  << " (tracker status code: " << ist << ")";
1560 
1561  output << "$ track " << pdgc << " " << E << " "
1562  << dcosx << " " << dcosy << " " << dcosz << " "
1563  << ist << endl;
1564 
1565  }//tracks
1566  }//iloop
1567 
1568  //
1569  // -- Add $info lines as necessary
1570  //
1571 
1572  if(gOptOutFileFormat == kConvFmt_t2k_tracker) {
1573  //
1574  // Writing $info lines with information identical to the one saved at the rootracker-format
1575  // files for the nd280MC. SKDETSIM can propagate all that complete MC truth information into
1576  // friend event trees that can be 'linked' with the SK DSTs.
1577  // Having identical generator info for both SK and nd280 will enable global studies
1578  //
1579  // The $info lines are formatted as follows:
1580  //
1581  // version 1:
1582  //
1583  // $ info event_num err_flag string_event_code
1584  // $ info xsec_event diff_xsec_kinematics weight prob
1585  // $ info vtxx vtxy vtxz vtxt
1586  // $ info nparticles
1587  // $ info 0 pdg_code status_code first_daughter last_daughter first_mother last_mother px py pz E x y z t polx poly polz
1588  // $ info 1 pdg_code status_code first_daughter last_daughter first_mother last_mother px py pz E x y z t polx poly polz
1589  // ... ... ...
1590  // $ info n pdg_code status_code first_daughter last_daughter first_mother last_mother px py pz E x y z t polx poly polz
1591  //
1592  // version 2:
1593  //
1594  // $ info event_num err_flag string_event_code
1595  // $ info xsec_event diff_xsec_kinematics weight prob
1596  // $ info vtxx vtxy vtxz vtxt
1597  // $ info etc
1598  // $ info nparticles
1599  // $ info 0 pdg_code status_code first_daughter last_daughter first_mother last_mother px py pz E x y z t polx poly polz rescatter_code
1600  // $ info 1 pdg_code status_code first_daughter last_daughter first_mother last_mother px py pz E x y z t polx poly polz rescatter_code
1601  // ... ... ...
1602  // $ info n pdg_code status_code first_daughter last_daughter first_mother last_mother px py pz E x y z t polx poly polz rescatter_code
1603  //
1604  // Comments:
1605  // - The err_flag is a bit field (16 bits)
1606  // - The string_event_code is a rather long string which encapsulates lot of summary info on the event
1607  // (neutrino/nuclear target/hit nucleon/hit quark(if any)/process type/...).
1608  // Information on how to parse that string code is available at the T2K event reweighting package.
1609  // - event_xsec is the event cross section in 1E-38cm^2
1610  // - diff_event_xsec is the cross section for the selected in 1E-38cm^2/{K^n}
1611  // - weight is the event weight (1 for unweighted MC)
1612  // - prob is the event probability (given cross sectios and density-weighted path-length)
1613  // - vtxx,y,z,t is the vertex position/time in SI units
1614  // - etc (added in format vrs >= 2) is used to pass any additional information with event-scope.
1615  // For the time being it is being used to pass the hit quark id (for DIS events) that was lost before
1616  // as SKDETSIM doesn't read the string_event_code where this info is nominally contained.
1617  // The quark id is set as (quark_pdg_code) x 10 + i, where i=0 for valence and i=1 for sea quarks. Set to -1 for non-DIS events.
1618  // - nparticles is the number of particles in the GHEP record (number of $info lines to follow before the start of the JNUBEAM block)
1619  // - first_/last_daughter first_/last_mother indicate the particle
1620  // - px,py,pz,E is the particle 4-momentum at the LAB frame (in GeV)
1621  // - x,y,z,t is the particle 4-position at the hit nucleus coordinate system (in fm, t is not set)
1622  // - polx,y,z is the particle polarization vector
1623  // - rescatter_code (added in format vrs >= 2) is a model-dependent intranuclear rescattering code
1624  // added to simplify the event analysis (although, in principle, it is recoverable from the particle record).
1625  // See $GENIE/src/HadronTransport/INukeHadroFates.h for the meaning of various codes when INTRANUKE is in use.
1626  // The rescattering code is stored at the GHEP event record for files generated with GENIE vrs >= 2.5.1.
1627  // See also ConvertToGRooTracker() for further descriptions of the variables stored at
1628  // the rootracker files.
1629  //
1630  // event info
1631  //
1632  output << "$ info " << (int) iev << " " << *(event.EventFlags()) << " " << interaction->AsString() << endl;
1633  output << "$ info " << (1E+38/units::cm2) * event.XSec() << " "
1634  << (1E+38/units::cm2) * event.DiffXSec() << " "
1635  << event.Weight() << " "
1636  << event.Probability()
1637  << endl;
1638  output << "$ info " << event.Vertex()->X() << " "
1639  << event.Vertex()->Y() << " "
1640  << event.Vertex()->Z() << " "
1641  << event.Vertex()->T()
1642  << endl;
1643 
1644  // insert etc info line for format versions >= 2
1645  if(gOptVersion >= 2) {
1646  int quark_id = -1;
1647  if( interaction->ProcInfo().IsDeepInelastic() && interaction->InitState().Tgt().HitQrkIsSet() ) {
1648  int quark_pdg = interaction->InitState().Tgt().HitQrkPdg();
1649  int sorv = ( interaction->InitState().Tgt().HitSeaQrk() ) ? 1 : 0; // sea q: 1, valence q: 0
1650  quark_id = 10 * quark_pdg + sorv;
1651  }
1652  output << "$ info " << quark_id << endl;
1653  }
1654 
1655  //
1656  // copy stdhep-like particle list
1657  //
1658  iparticle = 0;
1659  event_iter.Reset();
1660  output << "$ info " << event.GetEntries() << endl;
1661  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) )
1662  {
1663  assert(p);
1664  output << "$ info "
1665  << iparticle << " "
1666  << p->Pdg() << " " << (int) p->Status() << " "
1667  << p->FirstDaughter() << " " << p->LastDaughter() << " "
1668  << p->FirstMother() << " " << p->LastMother() << " "
1669  << p->X4()->X() << " " << p->X4()->Y() << " " << p->X4()->Z() << " " << p->X4()->T() << " "
1670  << p->P4()->Px() << " " << p->P4()->Py() << " " << p->P4()->Pz() << " " << p->P4()->E() << " ";
1671  if(p->PolzIsSet()) {
1672  output << TMath::Sin(p->PolzPolarAngle()) * TMath::Cos(p->PolzAzimuthAngle()) << " "
1673  << TMath::Sin(p->PolzPolarAngle()) * TMath::Sin(p->PolzAzimuthAngle()) << " "
1674  << TMath::Cos(p->PolzPolarAngle());
1675  } else {
1676  output << "0. 0. 0.";
1677  }
1678 
1679  // append rescattering code for format versions >= 2
1680  if(gOptVersion >= 2) {
1681  int rescat_code = -1;
1682  bool have_rescat_code = false;
1683  if(gFileMajorVrs >= 2) {
1684  if(gFileMinorVrs >= 5) {
1685  if(gFileRevisVrs >= 1) {
1686  have_rescat_code = true;
1687  }
1688  }
1689  }
1690  if(have_rescat_code) {
1691  rescat_code = p->RescatterCode();
1692  }
1693  output << " ";
1694  output << rescat_code;
1695  }
1696 
1697  output << endl;
1698  iparticle++;
1699  }
1700  //
1701  // JNUBEAM flux info - this info will only be available if events were generated
1702  // by gT2Kevgen using JNUBEAM flux ntuples as inputs
1703  //
1704 /*
1705 The T2K/SK collaboration produces MC based on JNUBEAM flux histograms, not flux ntuples.
1706 Therefore JNUBEAM flux pass-through info is never available for generated events.
1707 Commented-out the following info so as not to maintain/support unused code.
1708 If this section is ever re-instated the JNUBEAM passed-through info needs to be matched
1709 to the latest version of JNUBEAM and an appropriate updated t2k_tracker format needs to
1710 be agreed with the SKDETSIM maintainers.
1711 
1712 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
1713  PDGLibrary * pdglib = PDGLibrary::Instance();
1714  if(flux_info) {
1715  // parent hadron pdg code and decay mode
1716  output << "$ info " << pdg::GeantToPdg(flux_info->ppid) << " " << flux_info->mode << endl;
1717  // parent hadron px,py,pz,E at decay
1718  output << "$ info " << flux_info->ppi * flux_info->npi[0] << " "
1719  << flux_info->ppi * flux_info->npi[1] << " "
1720  << flux_info->ppi * flux_info->npi[2] << " "
1721  << TMath::Sqrt(
1722  TMath::Power(pdglib->Find(pdg::GeantToPdg(flux_info->ppid))->Mass(), 2.)
1723  + TMath::Power(flux_info->ppi, 2.)
1724  ) << endl;
1725  // parent hadron x,y,z,t at decay
1726  output << "$ info " << flux_info->xpi[0] << " "
1727  << flux_info->xpi[1] << " "
1728  << flux_info->xpi[2] << " "
1729  << "0."
1730  << endl;
1731  // parent hadron px,py,pz,E at production
1732  output << "$ info " << flux_info->ppi0 * flux_info->npi0[0] << " "
1733  << flux_info->ppi0 * flux_info->npi0[1] << " "
1734  << flux_info->ppi0 * flux_info->npi0[2] << " "
1735  << TMath::Sqrt(
1736  TMath::Power(pdglib->Find(pdg::GeantToPdg(flux_info->ppid))->Mass(), 2.)
1737  + TMath::Power(flux_info->ppi0, 2.)
1738  ) << endl;
1739  // parent hadron x,y,z,t at production
1740  output << "$ info " << flux_info->xpi0[0] << " "
1741  << flux_info->xpi0[1] << " "
1742  << flux_info->xpi0[2] << " "
1743  << "0."
1744  << endl;
1745  // nvtx
1746  output << "$ info " << output << "$info " << endl;
1747  }
1748 #endif
1749 */
1750  }//fmt==kConvFmt_t2k_tracker
1751 
1752  //
1753  // -- Add tracker end tag
1754  //
1755  output << "$ end" << endl;
1756 
1757  mcrec->Clear();
1758  } // event loop
1759 
1760  // add tracker end-of-file tag
1761  output << "$ stop" << endl;
1762 
1763  output.close();
1764  fin.Close();
1765 
1766  LOG("gntpc", pINFO) << "\nDone converting GENIE's GHEP ntuple";
1767 }
1768 //____________________________________________________________________________________
1769 // GENIE GHEP EVENT TREE FORMAT -> ROOTRACKER FORMATS
1770 //____________________________________________________________________________________
1771 void ConvertToGRooTracker(void)
1772 {
1773  //-- define the output rootracker tree branches
1774 
1775  // event info
1776 
1777  TBits* brEvtFlags = 0; // Generator-specific event flags
1778  TObjString* brEvtCode = 0; // Generator-specific string with 'event code'
1779  int brEvtNum; // Event num.
1780  double brEvtXSec; // Cross section for selected event (1E-38 cm2)
1781  double brEvtDXSec; // Cross section for selected event kinematics (1E-38 cm2 /{K^n})
1782  double brEvtWght; // Weight for that event
1783  double brEvtProb; // Probability for that event (given cross section, path lengths, etc)
1784  double brEvtVtx[4]; // Event vertex position in detector coord syst (SI)
1785  int brStdHepN; // Number of particles in particle array
1786  // stdhep-like particle array:
1787  int brStdHepPdg [kNPmax]; // Pdg codes (& generator specific codes for pseudoparticles)
1788  int brStdHepStatus[kNPmax]; // Generator-specific status code
1789  int brStdHepRescat[kNPmax]; // Hadron transport model - specific rescattering code
1790  double brStdHepX4 [kNPmax][4]; // 4-x (x, y, z, t) of particle in hit nucleus frame (fm)
1791  double brStdHepP4 [kNPmax][4]; // 4-p (px,py,pz,E) of particle in LAB frame (GeV)
1792  double brStdHepPolz [kNPmax][3]; // Polarization vector
1793  int brStdHepFd [kNPmax]; // First daughter
1794  int brStdHepLd [kNPmax]; // Last daughter
1795  int brStdHepFm [kNPmax]; // First mother
1796  int brStdHepLm [kNPmax]; // Last mother
1797 
1798  //
1799  // >> info available at the t2k rootracker variance only
1800  //
1801  TObjString* brNuFileName = 0; // flux file name
1802  long brNuFluxEntry; // entry number from flux file
1803 
1804  // neutrino parent info (passed-through from the beam-line MC / quantities in 'jnubeam' units)
1805  int brNuParentPdg; // parent hadron pdg code
1806  int brNuParentDecMode; // parent hadron decay mode
1807  double brNuParentDecP4 [4]; // parent hadron 4-momentum at decay
1808  double brNuParentDecX4 [4]; // parent hadron 4-position at decay
1809  double brNuParentProP4 [4]; // parent hadron 4-momentum at production
1810  double brNuParentProX4 [4]; // parent hadron 4-position at production
1811  int brNuParentProNVtx; // parent hadron vtx id
1812  // variables added since 10a flux compatibility changes
1813  int brNuIdfd; // detector location id
1814  float brNuCospibm; // cosine of the angle between the parent particle direction and the beam direction
1815  float brNuCospi0bm; // same as above except at the production of the parent particle
1816  int brNuGipart; // primary particle ID
1817  float brNuGpos0[3]; // primary particle starting point
1818  float brNuGvec0[3]; // primary particle direction at the starting point
1819  float brNuGamom0; // momentum of the primary particle at the starting point
1820  // variables added since 10d and 11a flux compatibility changes
1821  float brNuRnu; // neutrino r position at ND5/6 plane
1822  float brNuXnu[2]; // neutrino (x,y) position at ND5/6 plane
1823  // interation history information
1824  int brNuNg; // number of parents (number of generations)
1825  int brNuGpid[flux::fNgmax]; // particle ID of each ancestor particles
1826  int brNuGmec[flux::fNgmax]; // particle production mechanism of each ancestor particle
1827  float brNuGcosbm[flux::fNgmax]; // ancestor particle cos(theta) relative to beam
1828  float brNuGv[flux::fNgmax][3]; // X,Y and Z vertex position of each ancestor particle
1829  float brNuGp[flux::fNgmax][3]; // Px,Px and Pz directional momentum of each ancestor particle
1830  // out-of-target secondary interactions
1831  int brNuGmat[flux::fNgmax]; // material in which the particle originates
1832  float brNuGdistc[flux::fNgmax]; // distance traveled through carbon
1833  float brNuGdistal[flux::fNgmax]; // distance traveled through aluminum
1834  float brNuGdistti[flux::fNgmax]; // distance traveled through titanium
1835  float brNuGdistfe[flux::fNgmax]; // distance traveled through iron
1836 
1837  float brNuNorm; // normalisation weight (makes no sense to apply this when generating unweighted events)
1838  float brNuEnusk; // "Enu" for SK
1839  float brNuNormsk; // "norm" for SK
1840  float brNuAnorm; // Norm component from ND acceptance calculation
1841  float brNuVersion; // Jnubeam version
1842  int brNuNtrig; // Number of Triggers in simulation
1843  int brNuTuneid; // Parameter set identifier
1844  int brNuPint; // Interaction model ID
1845  float brNuBpos[2]; // Beam center position
1846  float brNuBtilt[2]; // Beam Direction
1847  float brNuBrms[2]; // Beam RMS Width
1848  float brNuEmit[2]; // Beam Emittance
1849  float brNuAlpha[2]; // Beam alpha parameter
1850  float brNuHcur[3]; // Horns 1, 2 and 3 Currents
1851  int brNuRand; // Random seed
1852  // codes for T2K cross-generator comparisons
1853  int brNeutCode; // NEUT-like reaction code for the GENIE event
1854 
1855  //
1856  // >> info available at the numi rootracker variance only
1857  //
1858 
1859  // neutrino parent info (GNuMI passed-through info)
1860  // see http://www.hep.utexas.edu/~zarko/wwwgnumi/v19/[/v19/output_gnumi.html]
1861  int brNumiFluxRun; // Run number
1862  int brNumiFluxEvtno; // Event number (proton on target)
1863  double brNumiFluxNdxdz; // Neutrino direction slope (dx/dz) for a random decay
1864  double brNumiFluxNdydz; // Neutrino direction slope (dy/dz) for a random decay
1865  double brNumiFluxNpz; // Neutrino momentum (GeV/c) along z direction (beam axis)
1866  double brNumiFluxNenergy; // Neutrino energy (GeV/c) for a random decay
1867  double brNumiFluxNdxdznea; // Neutrino direction slope (dx/dz) for a decay forced at center of near detector
1868  double brNumiFluxNdydznea; // Neutrino direction slope (dy/dz) for a decay forced at center of near detector
1869  double brNumiFluxNenergyn; // Neutrino energy for a decay forced at center of near detector
1870  double brNumiFluxNwtnear; // Neutrino weight for a decay forced at center of near detector
1871  double brNumiFluxNdxdzfar; // Neutrino direction slope (dx/dz) for a decay forced at center of far detector
1872  double brNumiFluxNdydzfar; // Neutrino direction slope (dy/dz) for a decay forced at center of far detector
1873  double brNumiFluxNenergyf; // Neutrino energy for a decay forced at center of far detector
1874  double brNumiFluxNwtfar; // Neutrino weight for a decay forced at center of far detector
1875  int brNumiFluxNorig; // Obsolete
1876  int brNumiFluxNdecay; // Decay mode that produced neutrino:
1877  // - 1 K0L -> nue pi- e+
1878  // - 2 K0L -> nuebar pi+ e-
1879  // - 3 K0L -> numu pi- mu+
1880  // - 4 K0L -> numubar pi+ mu-
1881  // - 5 K+ -> numu mu+
1882  // - 6 K+ -> nue pi0 e+
1883  // - 7 K+ -> numu pi0 mu+
1884  // - 8 K- -> numubar mu-
1885  // - 9 K- -> nuebar pi0 e-
1886  // - 10 K- -> numubar pi0 mu-
1887  // - 11 mu+ -> numubar nue e+
1888  // - 12 mu- -> numu nuebar e-
1889  // - 13 pi+ -> numu mu+
1890  // - 14 pi- -> numubar mu-
1891  int brNumiFluxNtype; // Neutrino flavor
1892  double brNumiFluxVx; // Position of hadron/muon decay, X coordinate
1893  double brNumiFluxVy; // Position of hadron/muon decay, Y coordinate
1894  double brNumiFluxVz; // Position of hadron/muon decay, Z coordinate
1895  double brNumiFluxPdpx; // Parent momentum at decay point, X - component
1896  double brNumiFluxPdpy; // Parent momentum at decay point, Y - component
1897  double brNumiFluxPdpz; // Parent momentum at decay point, Z - component
1898  double brNumiFluxPpdxdz; // Parent dx/dz direction at production
1899  double brNumiFluxPpdydz; // Parent dy/dz direction at production
1900  double brNumiFluxPppz; // Parent Z momentum at production
1901  double brNumiFluxPpenergy; // Parent energy at production
1902  int brNumiFluxPpmedium; // Tracking medium number where parent was produced
1903  int brNumiFluxPtype; // Parent particle ID (PDG)
1904  double brNumiFluxPpvx; // Parent production vertex, X coordinate (cm)
1905  double brNumiFluxPpvy; // Parent production vertex, Y coordinate (cm)
1906  double brNumiFluxPpvz; // Parent production vertex, Z coordinate (cm)
1907  double brNumiFluxMuparpx; // Repeat of information above, but for muon neutrino parents
1908  double brNumiFluxMuparpy; // ...
1909  double brNumiFluxMuparpz; // ...
1910  double brNumiFluxMupare; // ...
1911  double brNumiFluxNecm; // Neutrino energy in COM frame
1912  double brNumiFluxNimpwt; // Weight of neutrino parent
1913  double brNumiFluxXpoint; // Unused
1914  double brNumiFluxYpoint; // Unused
1915  double brNumiFluxZpoint; // Unused
1916  double brNumiFluxTvx; // Exit point of parent particle at the target, X coordinate
1917  double brNumiFluxTvy; // Exit point of parent particle at the target, Y coordinate
1918  double brNumiFluxTvz; // Exit point of parent particle at the target, Z coordinate
1919  double brNumiFluxTpx; // Parent momentum exiting the target, X - component
1920  double brNumiFluxTpy; // Parent momentum exiting the target, Y - component
1921  double brNumiFluxTpz; // Parent momentum exiting the target, Z - component
1922  double brNumiFluxTptype; // Parent particle ID exiting the target
1923  double brNumiFluxTgen; // Parent generation in cascade
1924  // - 1 primary proton
1925  // - 2 particles produced by proton interaction
1926  // - 3 particles produced by interactions of the 2's, ...
1927  double brNumiFluxTgptype; // Type of particle that created a particle flying of the target
1928  double brNumiFluxTgppx; // Momentum of a particle, that created a particle that flies off
1929  // the target (at the interaction point), X - component
1930  double brNumiFluxTgppy; // Momentum of a particle, that created a particle that flies off
1931  // the target (at the interaction point), Y - component
1932  double brNumiFluxTgppz; // Momentum of a particle, that created a particle that flies off
1933  // the target (at the interaction point), Z - component
1934  double brNumiFluxTprivx; // Primary particle interaction vertex, X coordinate
1935  double brNumiFluxTprivy; // Primary particle interaction vertex, Y coordinate
1936  double brNumiFluxTprivz; // Primary particle interaction vertex, Z coordinate
1937  double brNumiFluxBeamx; // Primary proton origin, X coordinate
1938  double brNumiFluxBeamy; // Primary proton origin, Y coordinate
1939  double brNumiFluxBeamz; // Primary proton origin, Z coordinate
1940  double brNumiFluxBeampx; // Primary proton momentum, X - component
1941  double brNumiFluxBeampy; // Primary proton momentum, Y - component
1942  double brNumiFluxBeampz; // Primary proton momentum, Z - component
1943 
1944  //-- open the output ROOT file
1945  TFile fout(gOptOutFileName.c_str(), "RECREATE");
1946 
1947  //-- create the output ROOT tree
1948  TTree * rootracker_tree = new TTree("gRooTracker","GENIE event tree rootracker format");
1949 
1950  //-- is it a `mock data' variance?
1951  bool hide_truth = (gOptOutFileFormat == kConvFmt_rootracker_mock_data);
1952 
1953  //-- create the output ROOT tree branches
1954 
1955  // branches common to all rootracker(_mock_data) formats
1956  if(!hide_truth) {
1957  // full version
1958  rootracker_tree->Branch("EvtFlags", "TBits", &brEvtFlags, 32000, 1);
1959  rootracker_tree->Branch("EvtCode", "TObjString", &brEvtCode, 32000, 1);
1960  rootracker_tree->Branch("EvtNum", &brEvtNum, "EvtNum/I");
1961  rootracker_tree->Branch("EvtXSec", &brEvtXSec, "EvtXSec/D");
1962  rootracker_tree->Branch("EvtDXSec", &brEvtDXSec, "EvtDXSec/D");
1963  rootracker_tree->Branch("EvtWght", &brEvtWght, "EvtWght/D");
1964  rootracker_tree->Branch("EvtProb", &brEvtProb, "EvtProb/D");
1965  rootracker_tree->Branch("EvtVtx", brEvtVtx, "EvtVtx[4]/D");
1966  rootracker_tree->Branch("StdHepN", &brStdHepN, "StdHepN/I");
1967  rootracker_tree->Branch("StdHepPdg", brStdHepPdg, "StdHepPdg[StdHepN]/I");
1968  rootracker_tree->Branch("StdHepStatus", brStdHepStatus, "StdHepStatus[StdHepN]/I");
1969  rootracker_tree->Branch("StdHepRescat", brStdHepRescat, "StdHepRescat[StdHepN]/I");
1970  rootracker_tree->Branch("StdHepX4", brStdHepX4, "StdHepX4[StdHepN][4]/D");
1971  rootracker_tree->Branch("StdHepP4", brStdHepP4, "StdHepP4[StdHepN][4]/D");
1972  rootracker_tree->Branch("StdHepPolz", brStdHepPolz, "StdHepPolz[StdHepN][3]/D");
1973  rootracker_tree->Branch("StdHepFd", brStdHepFd, "StdHepFd[StdHepN]/I");
1974  rootracker_tree->Branch("StdHepLd", brStdHepLd, "StdHepLd[StdHepN]/I");
1975  rootracker_tree->Branch("StdHepFm", brStdHepFm, "StdHepFm[StdHepN]/I");
1976  rootracker_tree->Branch("StdHepLm", brStdHepLm, "StdHepLm[StdHepN]/I");
1977  } else {
1978  // for mock_data variances
1979  rootracker_tree->Branch("EvtNum", &brEvtNum, "EvtNum/I");
1980  rootracker_tree->Branch("EvtWght", &brEvtWght, "EvtWght/D");
1981  rootracker_tree->Branch("EvtVtx", brEvtVtx, "EvtVtx[4]/D");
1982  rootracker_tree->Branch("StdHepN", &brStdHepN, "StdHepN/I");
1983  rootracker_tree->Branch("StdHepPdg", brStdHepPdg, "StdHepPdg[StdHepN]/I");
1984  rootracker_tree->Branch("StdHepX4", brStdHepX4, "StdHepX4[StdHepN][4]/D");
1985  rootracker_tree->Branch("StdHepP4", brStdHepP4, "StdHepP4[StdHepN][4]/D");
1986  }
1987 
1988  // extra branches of the t2k rootracker variance
1989  if(gOptOutFileFormat == kConvFmt_t2k_rootracker)
1990  {
1991  // NEUT-like reaction code
1992  rootracker_tree->Branch("G2NeutEvtCode", &brNeutCode, "G2NeutEvtCode/I");
1993  // JNUBEAM pass-through info
1994  rootracker_tree->Branch("NuFileName", "TObjString", &brNuFileName, 32000, 1);
1995  rootracker_tree->Branch("NuParentPdg", &brNuParentPdg, "NuParentPdg/I");
1996  rootracker_tree->Branch("NuParentDecMode", &brNuParentDecMode, "NuParentDecMode/I");
1997  rootracker_tree->Branch("NuParentDecP4", brNuParentDecP4, "NuParentDecP4[4]/D");
1998  rootracker_tree->Branch("NuParentDecX4", brNuParentDecX4, "NuParentDecX4[4]/D");
1999  rootracker_tree->Branch("NuParentProP4", brNuParentProP4, "NuParentProP4[4]/D");
2000  rootracker_tree->Branch("NuParentProX4", brNuParentProX4, "NuParentProX4[4]/D");
2001  rootracker_tree->Branch("NuParentProNVtx", &brNuParentProNVtx, "NuParentProNVtx/I");
2002  // Branches added since JNUBEAM '10a' compatibility changes
2003  rootracker_tree->Branch("NuFluxEntry", &brNuFluxEntry, "NuFluxEntry/L");
2004  rootracker_tree->Branch("NuIdfd", &brNuIdfd, "NuIdfd/I");
2005  rootracker_tree->Branch("NuCospibm", &brNuCospibm, "NuCospibm/F");
2006  rootracker_tree->Branch("NuCospi0bm", &brNuCospi0bm, "NuCospi0bm/F");
2007  rootracker_tree->Branch("NuGipart", &brNuGipart, "NuGipart/I");
2008  rootracker_tree->Branch("NuGpos0", brNuGpos0, "NuGpos0[3]/F");
2009  rootracker_tree->Branch("NuGvec0", brNuGvec0, "NuGvec0[3]/F");
2010  rootracker_tree->Branch("NuGamom0", &brNuGamom0, "NuGamom0/F");
2011  // Branches added since JNUBEAM '10d' compatibility changes
2012  rootracker_tree->Branch("NuXnu", brNuXnu, "NuXnu[2]/F");
2013  rootracker_tree->Branch("NuRnu", &brNuRnu, "NuRnu/F");
2014  rootracker_tree->Branch("NuNg", &brNuNg, "NuNg/I");
2015  rootracker_tree->Branch("NuGpid", brNuGpid, "NuGpid[NuNg]/I");
2016  rootracker_tree->Branch("NuGmec", brNuGmec, "NuGmec[NuNg]/I");
2017  rootracker_tree->Branch("NuGv", brNuGv, "NuGv[NuNg][3]/F");
2018  rootracker_tree->Branch("NuGp", brNuGp, "NuGp[NuNg][3]/F");
2019  rootracker_tree->Branch("NuGcosbm", brNuGcosbm, "NuGcosbm[NuNg]/F");
2020  rootracker_tree->Branch("NuGmat", brNuGmat, "NuGmat[NuNg]/I");
2021  rootracker_tree->Branch("NuGdistc", brNuGdistc, "NuGdistc[NuNg]/F");
2022  rootracker_tree->Branch("NuGdistal", brNuGdistal, "NuGdistal[NuNg]/F");
2023  rootracker_tree->Branch("NuGdistti", brNuGdistti, "NuGdistti[NuNg]/F");
2024  rootracker_tree->Branch("NuGdistfe", brNuGdistfe, "NuGdistfe[NuNg]/F");
2025  rootracker_tree->Branch("NuNorm", &brNuNorm, "NuNorm/F");
2026  rootracker_tree->Branch("NuEnusk", &brNuEnusk, "NuEnusk/F");
2027  rootracker_tree->Branch("NuNormsk", &brNuNormsk, "NuNormsk/F");
2028  rootracker_tree->Branch("NuAnorm", &brNuAnorm, "NuAnorm/F");
2029  rootracker_tree->Branch("NuVersion", &brNuVersion, "NuVersion/F");
2030  rootracker_tree->Branch("NuNtrig", &brNuNtrig, "NuNtrig/I");
2031  rootracker_tree->Branch("NuTuneid", &brNuTuneid, "NuTuneid/I");
2032  rootracker_tree->Branch("NuPint", &brNuPint, "NuPint/I");
2033  rootracker_tree->Branch("NuBpos", brNuBpos, "NuBpos[2]/F");
2034  rootracker_tree->Branch("NuBtilt", brNuBtilt, "NuBtilt[2]/F");
2035  rootracker_tree->Branch("NuBrms", brNuBrms, "NuBrms[2]/F");
2036  rootracker_tree->Branch("NuEmit", brNuEmit, "NuEmit[2]/F");
2037  rootracker_tree->Branch("NuAlpha", brNuAlpha, "NuAlpha[2]/F");
2038  rootracker_tree->Branch("NuHcur", brNuHcur, "NuHcur[3]/F");
2039  rootracker_tree->Branch("NuRand", &brNuRand, "NuRand/I");
2040 
2041  }
2042 
2043  // extra branches of the numi rootracker variance
2044  if(gOptOutFileFormat == kConvFmt_numi_rootracker)
2045  {
2046  // GNuMI pass-through info
2047  rootracker_tree->Branch("NumiFluxRun", &brNumiFluxRun, "NumiFluxRun/I");
2048  rootracker_tree->Branch("NumiFluxEvtno", &brNumiFluxEvtno, "NumiFluxEvtno/I");
2049  rootracker_tree->Branch("NumiFluxNdxdz", &brNumiFluxNdxdz, "NumiFluxNdxdz/D");
2050  rootracker_tree->Branch("NumiFluxNdydz", &brNumiFluxNdydz, "NumiFluxNdydz/D");
2051  rootracker_tree->Branch("NumiFluxNpz", &brNumiFluxNpz, "NumiFluxNpz/D");
2052  rootracker_tree->Branch("NumiFluxNenergy", &brNumiFluxNenergy, "NumiFluxNenergy/D");
2053  rootracker_tree->Branch("NumiFluxNdxdznea", &brNumiFluxNdxdznea, "NumiFluxNdxdznea/D");
2054  rootracker_tree->Branch("NumiFluxNdydznea", &brNumiFluxNdydznea, "NumiFluxNdydznea/D");
2055  rootracker_tree->Branch("NumiFluxNenergyn", &brNumiFluxNenergyn, "NumiFluxNenergyn/D");
2056  rootracker_tree->Branch("NumiFluxNwtnear", &brNumiFluxNwtnear, "NumiFluxNwtnear/D");
2057  rootracker_tree->Branch("NumiFluxNdxdzfar", &brNumiFluxNdxdzfar, "NumiFluxNdxdzfar/D");
2058  rootracker_tree->Branch("NumiFluxNdydzfar", &brNumiFluxNdydzfar, "NumiFluxNdydzfar/D");
2059  rootracker_tree->Branch("NumiFluxNenergyf", &brNumiFluxNenergyf, "NumiFluxNenergyf/D");
2060  rootracker_tree->Branch("NumiFluxNwtfar", &brNumiFluxNwtfar, "NumiFluxNwtfar/D");
2061  rootracker_tree->Branch("NumiFluxNorig", &brNumiFluxNorig, "NumiFluxNorig/I");
2062  rootracker_tree->Branch("NumiFluxNdecay", &brNumiFluxNdecay, "NumiFluxNdecay/I");
2063  rootracker_tree->Branch("NumiFluxNtype", &brNumiFluxNtype, "NumiFluxNtype/I");
2064  rootracker_tree->Branch("NumiFluxVx", &brNumiFluxVx, "NumiFluxVx/D");
2065  rootracker_tree->Branch("NumiFluxVy", &brNumiFluxVy, "NumiFluxVy/D");
2066  rootracker_tree->Branch("NumiFluxVz", &brNumiFluxVz, "NumiFluxVz/D");
2067  rootracker_tree->Branch("NumiFluxPdpx", &brNumiFluxPdpx, "NumiFluxPdpx/D");
2068  rootracker_tree->Branch("NumiFluxPdpy", &brNumiFluxPdpy, "NumiFluxPdpy/D");
2069  rootracker_tree->Branch("NumiFluxPdpz", &brNumiFluxPdpz, "NumiFluxPdpz/D");
2070  rootracker_tree->Branch("NumiFluxPpdxdz", &brNumiFluxPpdxdz, "NumiFluxPpdxdz/D");
2071  rootracker_tree->Branch("NumiFluxPpdydz", &brNumiFluxPpdydz, "NumiFluxPpdydz/D");
2072  rootracker_tree->Branch("NumiFluxPppz", &brNumiFluxPppz, "NumiFluxPppz/D");
2073  rootracker_tree->Branch("NumiFluxPpenergy", &brNumiFluxPpenergy, "NumiFluxPpenergy/D");
2074  rootracker_tree->Branch("NumiFluxPpmedium", &brNumiFluxPpmedium, "NumiFluxPpmedium/I");
2075  rootracker_tree->Branch("NumiFluxPtype", &brNumiFluxPtype, "NumiFluxPtype/I");
2076  rootracker_tree->Branch("NumiFluxPpvx", &brNumiFluxPpvx, "NumiFluxPpvx/D");
2077  rootracker_tree->Branch("NumiFluxPpvy", &brNumiFluxPpvy, "NumiFluxPpvy/D");
2078  rootracker_tree->Branch("NumiFluxPpvz", &brNumiFluxPpvz, "NumiFluxPpvz/D");
2079  rootracker_tree->Branch("NumiFluxMuparpx", &brNumiFluxMuparpx, "NumiFluxMuparpx/D");
2080  rootracker_tree->Branch("NumiFluxMuparpy", &brNumiFluxMuparpy, "NumiFluxMuparpy/D");
2081  rootracker_tree->Branch("NumiFluxMuparpz", &brNumiFluxMuparpz, "NumiFluxMuparpz/D");
2082  rootracker_tree->Branch("NumiFluxMupare", &brNumiFluxMupare, "NumiFluxMupare/D");
2083  rootracker_tree->Branch("NumiFluxNecm", &brNumiFluxNecm, "NumiFluxNecm/D");
2084  rootracker_tree->Branch("NumiFluxNimpwt", &brNumiFluxNimpwt, "NumiFluxNimpwt/D");
2085  rootracker_tree->Branch("NumiFluxXpoint", &brNumiFluxXpoint, "NumiFluxXpoint/D");
2086  rootracker_tree->Branch("NumiFluxYpoint", &brNumiFluxYpoint, "NumiFluxYpoint/D");
2087  rootracker_tree->Branch("NumiFluxZpoint", &brNumiFluxZpoint, "NumiFluxZpoint/D");
2088  rootracker_tree->Branch("NumiFluxTvx", &brNumiFluxTvx, "NumiFluxTvx/D");
2089  rootracker_tree->Branch("NumiFluxTvy", &brNumiFluxTvy, "NumiFluxTvy/D");
2090  rootracker_tree->Branch("NumiFluxTvz", &brNumiFluxTvz, "NumiFluxTvz/D");
2091  rootracker_tree->Branch("NumiFluxTpx", &brNumiFluxTpx, "NumiFluxTpx/D");
2092  rootracker_tree->Branch("NumiFluxTpy", &brNumiFluxTpy, "NumiFluxTpy/D");
2093  rootracker_tree->Branch("NumiFluxTpz", &brNumiFluxTpz, "NumiFluxTpz/D");
2094  rootracker_tree->Branch("NumiFluxTptype", &brNumiFluxTptype, "NumiFluxTptype/I");
2095  rootracker_tree->Branch("NumiFluxTgen", &brNumiFluxTgen, "NumiFluxTgen/I");
2096  rootracker_tree->Branch("NumiFluxTgptype", &brNumiFluxTgptype, "NumiFluxTgptype/I");
2097  rootracker_tree->Branch("NumiFluxTgppx", &brNumiFluxTgppx, "NumiFluxTgppx/D");
2098  rootracker_tree->Branch("NumiFluxTgppy", &brNumiFluxTgppy, "NumiFluxTgppy/D");
2099  rootracker_tree->Branch("NumiFluxTgppz", &brNumiFluxTgppz, "NumiFluxTgppz/D");
2100  rootracker_tree->Branch("NumiFluxTprivx", &brNumiFluxTprivx, "NumiFluxTprivx/D");
2101  rootracker_tree->Branch("NumiFluxTprivy", &brNumiFluxTprivy, "NumiFluxTprivy/D");
2102  rootracker_tree->Branch("NumiFluxTprivz", &brNumiFluxTprivz, "NumiFluxTprivz/D");
2103  rootracker_tree->Branch("NumiFluxBeamx", &brNumiFluxBeamx, "NumiFluxBeamx/D");
2104  rootracker_tree->Branch("NumiFluxBeamy", &brNumiFluxBeamy, "NumiFluxBeamy/D");
2105  rootracker_tree->Branch("NumiFluxBeamz", &brNumiFluxBeamz, "NumiFluxBeamz/D");
2106  rootracker_tree->Branch("NumiFluxBeampx", &brNumiFluxBeampx, "NumiFluxBeampx/D");
2107  rootracker_tree->Branch("NumiFluxBeampy", &brNumiFluxBeampy, "NumiFluxBeampy/D");
2108  rootracker_tree->Branch("NumiFluxBeampz", &brNumiFluxBeampz, "NumiFluxBeampz/D");
2109  }
2110 
2111  //-- open the input GENIE ROOT file and get the TTree & its header
2112  TFile fin(gOptInpFileName.c_str(),"READ");
2113  TTree * gtree = 0;
2114  NtpMCTreeHeader * thdr = 0;
2115  gtree = dynamic_cast <TTree *> ( fin.Get("gtree") );
2116  thdr = dynamic_cast <NtpMCTreeHeader *> ( fin.Get("header") );
2117 
2118  LOG("gntpc", pINFO) << "Input tree header: " << *thdr;
2119 
2120  //-- get mc record
2121  NtpMCEventRecord * mcrec = 0;
2122  gtree->SetBranchAddress("gmcrec", &mcrec);
2123 
2124  //-- print-out metadata associated with the input event file in case the
2125  // event file was generated using the gT2Kevgen driver
2126  // (assuming this is the case if the requested output format is the t2k_rootracker format)
2127  if(gOptOutFileFormat == kConvFmt_t2k_rootracker)
2128  {
2129  // Check can find the MetaData
2130  genie::utils::T2KEvGenMetaData * metadata = NULL;
2131  metadata = (genie::utils::T2KEvGenMetaData *) gtree->GetUserInfo()->At(0);
2132  if(metadata){
2133  LOG("gntpc", pINFO) << "Found T2KMetaData!";
2134  LOG("gntpc", pINFO) << *metadata;
2135  }
2136  else {
2137  LOG("gntpc", pWARN)
2138  << "Could not find T2KMetaData attached to the event tree!";
2139  }
2140  }
2141 
2142 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
2143  flux::GJPARCNuFluxPassThroughInfo * jnubeam_flux_info = 0;
2144  if(gOptOutFileFormat == kConvFmt_t2k_rootracker) {
2145  gtree->SetBranchAddress("flux", &jnubeam_flux_info);
2146  }
2147  flux::GNuMIFluxPassThroughInfo * gnumi_flux_info = 0;
2148  if(gOptOutFileFormat == kConvFmt_numi_rootracker) {
2149  gtree->SetBranchAddress("flux", &gnumi_flux_info);
2150  }
2151 #else
2152  LOG("gntpc", pWARN)
2153  << "\n Flux drivers are not enabled."
2154  << "\n No flux pass-through information will be written-out in the rootracker file"
2155  << "\n If this isn't what you are supposed to be doing then build GENIE by adding "
2156  << "--with-flux-drivers in the configuration step.";
2157 #endif
2158 
2159  //-- figure out how many events to analyze
2160  Long64_t nmax = (gOptN<0) ?
2161  gtree->GetEntries() : TMath::Min(gtree->GetEntries(), gOptN);
2162  if (nmax<0) {
2163  LOG("gntpc", pERROR) << "Number of events = 0";
2164  return;
2165  }
2166  LOG("gntpc", pNOTICE) << "*** Analyzing: " << nmax << " events";
2167 
2168  //-- event loop
2169  for(Long64_t iev = 0; iev < nmax; iev++) {
2170  gtree->GetEntry(iev);
2171 
2172  NtpMCRecHeader rec_header = mcrec->hdr;
2173  EventRecord & event = *(mcrec->event);
2174  Interaction * interaction = event.Summary();
2175 
2176  LOG("gntpc", pINFO) << rec_header;
2177  LOG("gntpc", pINFO) << event;
2178  LOG("gntpc", pINFO) << *interaction;
2179 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
2180  if(gOptOutFileFormat == kConvFmt_t2k_rootracker) {
2181  if(jnubeam_flux_info) {
2182  LOG("gntpc", pINFO) << *jnubeam_flux_info;
2183  } else {
2184  LOG("gntpc", pINFO) << "No JNUBEAM flux info associated with this event";
2185  }
2186  }
2187 #endif
2188 
2189  //
2190  // clear output tree branches
2191  //
2192  if(brEvtFlags) delete brEvtFlags;
2193  brEvtFlags = 0;
2194  if(brEvtCode) delete brEvtCode;
2195  brEvtCode = 0;
2196  brEvtNum = 0;
2197  brEvtXSec = 0;
2198  brEvtDXSec = 0;
2199  brEvtWght = 0;
2200  brEvtProb = 0;
2201  for(int k=0; k<4; k++) {
2202  brEvtVtx[k] = 0;
2203  }
2204  brStdHepN = 0;
2205  for(int i=0; i<kNPmax; i++) {
2206  brStdHepPdg [i] = 0;
2207  brStdHepStatus[i] = -1;
2208  brStdHepRescat[i] = -1;
2209  for(int k=0; k<4; k++) {
2210  brStdHepX4 [i][k] = 0;
2211  brStdHepP4 [i][k] = 0;
2212  }
2213  for(int k=0; k<3; k++) {
2214  brStdHepPolz [i][k] = 0;
2215  }
2216  brStdHepFd [i] = 0;
2217  brStdHepLd [i] = 0;
2218  brStdHepFm [i] = 0;
2219  brStdHepLm [i] = 0;
2220  }
2221  brNuParentPdg = 0;
2222  brNuParentDecMode = 0;
2223  for(int k=0; k<4; k++) {
2224  brNuParentDecP4 [k] = 0;
2225  brNuParentDecX4 [k] = 0;
2226  brNuParentProP4 [k] = 0;
2227  brNuParentProX4 [k] = 0;
2228  }
2229  brNuParentProNVtx = 0;
2230  brNeutCode = 0;
2231  brNuFluxEntry = -1;
2232  brNuIdfd = -999999;
2233  brNuCospibm = -999999.;
2234  brNuCospi0bm = -999999.;
2235  brNuGipart = -1;
2236  brNuGamom0 = -999999.;
2237  for(int k=0; k< 3; k++){
2238  brNuGvec0[k] = -999999.;
2239  brNuGpos0[k] = -999999.;
2240  }
2241  // variables added since 10d flux compatibility changes
2242  for(int k=0; k<2; k++) {
2243  brNuXnu[k] = brNuBpos[k] = brNuBtilt[k] = brNuBrms[k] = brNuEmit[k] = brNuAlpha[k] = -999999.;
2244  }
2245  for(int k=0; k<3; k++) brNuHcur[k] = -999999.;
2246  for(int np = 0; np < flux::fNgmax; np++){
2247  for(int k=0; k<3; k++){
2248  brNuGv[np][k] = -999999.;
2249  brNuGp[np][k] = -999999.;
2250  }
2251  brNuGpid[np] = -999999;
2252  brNuGmec[np] = -999999;
2253  brNuGmat[np] = -999999;
2254  brNuGcosbm[np] = -999999.;
2255  brNuGdistc[np] = -999999.;
2256  brNuGdistal[np] = -999999.;
2257  brNuGdistti[np] = -999999.;
2258  brNuGdistfe[np] = -999999.;
2259  }
2260  brNuNg = -999999;
2261  brNuRnu = -999999.;
2262  brNuNorm = -999999.;
2263  brNuEnusk = -999999.;
2264  brNuNormsk = -999999.;
2265  brNuAnorm = -999999.;
2266  brNuVersion= -999999.;
2267  brNuNtrig = -999999;
2268  brNuTuneid = -999999;
2269  brNuPint = -999999;
2270  brNuRand = -999999;
2271  if(brNuFileName) delete brNuFileName;
2272  brNuFileName = 0;
2273 
2274  //
2275  // copy current event info to output tree
2276  //
2277 
2278  brEvtFlags = new TBits(*event.EventFlags());
2279  brEvtCode = new TObjString(event.Summary()->AsString().c_str());
2280  brEvtNum = (int) iev;
2281  brEvtXSec = (1E+38/units::cm2) * event.XSec();
2282  brEvtDXSec = (1E+38/units::cm2) * event.DiffXSec();
2283  brEvtWght = event.Weight();
2284  brEvtProb = event.Probability();
2285  brEvtVtx[0] = event.Vertex()->X();
2286  brEvtVtx[1] = event.Vertex()->Y();
2287  brEvtVtx[2] = event.Vertex()->Z();
2288  brEvtVtx[3] = event.Vertex()->T();
2289 
2290  int iparticle=0;
2291  GHepParticle * p = 0;
2292  TIter event_iter(&event);
2293  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
2294  assert(p);
2295 
2296  // for mock_data variances write out only stable final state particles
2297  if(hide_truth && p->Status() != kIStStableFinalState) continue;
2298 
2299  brStdHepPdg [iparticle] = p->Pdg();
2300  brStdHepStatus[iparticle] = (int) p->Status();
2301  brStdHepRescat[iparticle] = p->RescatterCode();
2302  brStdHepX4 [iparticle][0] = p->X4()->X();
2303  brStdHepX4 [iparticle][1] = p->X4()->Y();
2304  brStdHepX4 [iparticle][2] = p->X4()->Z();
2305  brStdHepX4 [iparticle][3] = p->X4()->T();
2306  brStdHepP4 [iparticle][0] = p->P4()->Px();
2307  brStdHepP4 [iparticle][1] = p->P4()->Py();
2308  brStdHepP4 [iparticle][2] = p->P4()->Pz();
2309  brStdHepP4 [iparticle][3] = p->P4()->E();
2310  if(p->PolzIsSet()) {
2311  brStdHepPolz [iparticle][0] = TMath::Sin(p->PolzPolarAngle()) * TMath::Cos(p->PolzAzimuthAngle());
2312  brStdHepPolz [iparticle][1] = TMath::Sin(p->PolzPolarAngle()) * TMath::Sin(p->PolzAzimuthAngle());
2313  brStdHepPolz [iparticle][2] = TMath::Cos(p->PolzPolarAngle());
2314  }
2315  brStdHepFd [iparticle] = p->FirstDaughter();
2316  brStdHepLd [iparticle] = p->LastDaughter();
2317  brStdHepFm [iparticle] = p->FirstMother();
2318  brStdHepLm [iparticle] = p->LastMother();
2319  iparticle++;
2320  }
2321  brStdHepN = iparticle;
2322 
2323  //
2324  // fill in additional info for the t2k_rootracker format
2325  //
2326  if(gOptOutFileFormat == kConvFmt_t2k_rootracker) {
2327 
2328  // map GENIE event to NEUT reaction codes
2329  brNeutCode = utils::ghep::NeutReactionCode(&event);
2330 
2331 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
2332  // Copy flux info if this is the t2k rootracker variance.
2333  // The flux may not be available, eg if events were generated using plain flux
2334  // histograms and not the JNUBEAM simulation's output flux ntuples.
2335  PDGLibrary * pdglib = PDGLibrary::Instance();
2336  if(jnubeam_flux_info) {
2337  brNuParentPdg = pdg::GeantToPdg(jnubeam_flux_info->ppid);
2338  brNuParentDecMode = jnubeam_flux_info->mode;
2339 
2340  brNuParentDecP4 [0] = jnubeam_flux_info->ppi * jnubeam_flux_info->npi[0]; // px
2341  brNuParentDecP4 [1] = jnubeam_flux_info->ppi * jnubeam_flux_info->npi[1]; // py
2342  brNuParentDecP4 [2] = jnubeam_flux_info->ppi * jnubeam_flux_info->npi[2]; // px
2343  brNuParentDecP4 [3] = TMath::Sqrt(
2344  TMath::Power(pdglib->Find(brNuParentPdg)->Mass(), 2.)
2345  + TMath::Power(jnubeam_flux_info->ppi, 2.)
2346  ); // E
2347  brNuParentDecX4 [0] = jnubeam_flux_info->xpi[0]; // x
2348  brNuParentDecX4 [1] = jnubeam_flux_info->xpi[1]; // y
2349  brNuParentDecX4 [2] = jnubeam_flux_info->xpi[2]; // x
2350  brNuParentDecX4 [3] = 0; // t
2351 
2352  brNuParentProP4 [0] = jnubeam_flux_info->ppi0 * jnubeam_flux_info->npi0[0]; // px
2353  brNuParentProP4 [1] = jnubeam_flux_info->ppi0 * jnubeam_flux_info->npi0[1]; // py
2354  brNuParentProP4 [2] = jnubeam_flux_info->ppi0 * jnubeam_flux_info->npi0[2]; // px
2355  brNuParentProP4 [3] = TMath::Sqrt(
2356  TMath::Power(pdglib->Find(brNuParentPdg)->Mass(), 2.)
2357  + TMath::Power(jnubeam_flux_info->ppi0, 2.)
2358  ); // E
2359  brNuParentProX4 [0] = jnubeam_flux_info->xpi0[0]; // x
2360  brNuParentProX4 [1] = jnubeam_flux_info->xpi0[1]; // y
2361  brNuParentProX4 [2] = jnubeam_flux_info->xpi0[2]; // x
2362  brNuParentProX4 [3] = 0; // t
2363 
2364  brNuParentProNVtx = jnubeam_flux_info->nvtx0;
2365 
2366  // Copy info added post JNUBEAM '10a' compatibility changes
2367  brNuFluxEntry = jnubeam_flux_info->fluxentry;
2368  brNuIdfd = jnubeam_flux_info->idfd;
2369  brNuCospibm = jnubeam_flux_info->cospibm;
2370  brNuCospi0bm = jnubeam_flux_info->cospi0bm;
2371  brNuGipart = jnubeam_flux_info->gipart;
2372  brNuGamom0 = jnubeam_flux_info->gamom0;
2373  for(int k=0; k<3; k++){
2374  brNuGpos0[k] = (double) jnubeam_flux_info->gpos0[k];
2375  brNuGvec0[k] = (double) jnubeam_flux_info->gvec0[k];
2376  }
2377  // Copy info added post JNUBEAM '10d' compatibility changes
2378  brNuXnu[0] = (double) jnubeam_flux_info->xnu;
2379  brNuXnu[1] = (double) jnubeam_flux_info->ynu;
2380  brNuRnu = (double) jnubeam_flux_info->rnu;
2381  for(int k=0; k<2; k++){
2382  brNuBpos[k] = (double) jnubeam_flux_info->bpos[k];
2383  brNuBtilt[k] = (double) jnubeam_flux_info->btilt[k];
2384  brNuBrms[k] = (double) jnubeam_flux_info->brms[k];
2385  brNuEmit[k] = (double) jnubeam_flux_info->emit[k];
2386  brNuAlpha[k] = (double) jnubeam_flux_info->alpha[k];
2387  }
2388  for(int k=0; k<3; k++) brNuHcur[k] = jnubeam_flux_info->hcur[k];
2389  for(int np = 0; np < flux::fNgmax; np++){
2390  brNuGv[np][0] = jnubeam_flux_info->gvx[np];
2391  brNuGv[np][1] = jnubeam_flux_info->gvy[np];
2392  brNuGv[np][2] = jnubeam_flux_info->gvz[np];
2393  brNuGp[np][0] = jnubeam_flux_info->gpx[np];
2394  brNuGp[np][1] = jnubeam_flux_info->gpy[np];
2395  brNuGp[np][2] = jnubeam_flux_info->gpz[np];
2396  brNuGpid[np] = jnubeam_flux_info->gpid[np];
2397  brNuGmec[np] = jnubeam_flux_info->gmec[np];
2398  brNuGcosbm[np] = jnubeam_flux_info->gcosbm[np];
2399  brNuGmat[np] = jnubeam_flux_info->gmat[np];
2400  brNuGdistc[np] = jnubeam_flux_info->gdistc[np];
2401  brNuGdistal[np] = jnubeam_flux_info->gdistal[np];
2402  brNuGdistti[np] = jnubeam_flux_info->gdistti[np];
2403  brNuGdistfe[np] = jnubeam_flux_info->gdistfe[np];
2404  }
2405  brNuNg = jnubeam_flux_info->ng;
2406  brNuNorm = jnubeam_flux_info->norm;
2407  brNuEnusk = jnubeam_flux_info->Enusk;
2408  brNuNormsk = jnubeam_flux_info->normsk;
2409  brNuAnorm = jnubeam_flux_info->anorm;
2410  brNuVersion= jnubeam_flux_info->version;
2411  brNuNtrig = jnubeam_flux_info->ntrig;
2412  brNuTuneid = jnubeam_flux_info->tuneid;
2413  brNuPint = jnubeam_flux_info->pint;
2414  brNuRand = jnubeam_flux_info->rand;
2415  brNuFileName = new TObjString(jnubeam_flux_info->fluxfilename.c_str());
2416  }//jnubeam_flux_info
2417 #endif
2418  }//kConvFmt_t2k_rootracker
2419 
2420  //
2421  // fill in additional info for the numi_rootracker format
2422  //
2423  if(gOptOutFileFormat == kConvFmt_numi_rootracker) {
2424 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
2425  // Copy flux info if this is the numi rootracker variance.
2426  if(gnumi_flux_info) {
2427  brNumiFluxRun = gnumi_flux_info->run;
2428  brNumiFluxEvtno = gnumi_flux_info->evtno;
2429  brNumiFluxNdxdz = gnumi_flux_info->ndxdz;
2430  brNumiFluxNdydz = gnumi_flux_info->ndydz;
2431  brNumiFluxNpz = gnumi_flux_info->npz;
2432  brNumiFluxNenergy = gnumi_flux_info->nenergy;
2433  brNumiFluxNdxdznea = gnumi_flux_info->ndxdznea;
2434  brNumiFluxNdydznea = gnumi_flux_info->ndydznea;
2435  brNumiFluxNenergyn = gnumi_flux_info->nenergyn;
2436  brNumiFluxNwtnear = gnumi_flux_info->nwtnear;
2437  brNumiFluxNdxdzfar = gnumi_flux_info->ndxdzfar;
2438  brNumiFluxNdydzfar = gnumi_flux_info->ndydzfar;
2439  brNumiFluxNenergyf = gnumi_flux_info->nenergyf;
2440  brNumiFluxNwtfar = gnumi_flux_info->nwtfar;
2441  brNumiFluxNorig = gnumi_flux_info->norig;
2442  brNumiFluxNdecay = gnumi_flux_info->ndecay;
2443  brNumiFluxNtype = gnumi_flux_info->ntype;
2444  brNumiFluxVx = gnumi_flux_info->vx;
2445  brNumiFluxVy = gnumi_flux_info->vy;
2446  brNumiFluxVz = gnumi_flux_info->vz;
2447  brNumiFluxPdpx = gnumi_flux_info->pdpx;
2448  brNumiFluxPdpy = gnumi_flux_info->pdpy;
2449  brNumiFluxPdpz = gnumi_flux_info->pdpz;
2450  brNumiFluxPpdxdz = gnumi_flux_info->ppdxdz;
2451  brNumiFluxPpdydz = gnumi_flux_info->ppdydz;
2452  brNumiFluxPppz = gnumi_flux_info->pppz;
2453  brNumiFluxPpenergy = gnumi_flux_info->ppenergy;
2454  brNumiFluxPpmedium = gnumi_flux_info->ppmedium;
2455  brNumiFluxPtype = gnumi_flux_info->ptype;
2456  brNumiFluxPpvx = gnumi_flux_info->ppvx;
2457  brNumiFluxPpvy = gnumi_flux_info->ppvy;
2458  brNumiFluxPpvz = gnumi_flux_info->ppvz;
2459  brNumiFluxMuparpx = gnumi_flux_info->muparpx;
2460  brNumiFluxMuparpy = gnumi_flux_info->muparpy;
2461  brNumiFluxMuparpz = gnumi_flux_info->muparpz;
2462  brNumiFluxMupare = gnumi_flux_info->mupare;
2463  brNumiFluxNecm = gnumi_flux_info->necm;
2464  brNumiFluxNimpwt = gnumi_flux_info->nimpwt;
2465  brNumiFluxXpoint = gnumi_flux_info->xpoint;
2466  brNumiFluxYpoint = gnumi_flux_info->ypoint;
2467  brNumiFluxZpoint = gnumi_flux_info->zpoint;
2468  brNumiFluxTvx = gnumi_flux_info->tvx;
2469  brNumiFluxTvy = gnumi_flux_info->tvy;
2470  brNumiFluxTvz = gnumi_flux_info->tvz;
2471  brNumiFluxTpx = gnumi_flux_info->tpx;
2472  brNumiFluxTpy = gnumi_flux_info->tpy;
2473  brNumiFluxTpz = gnumi_flux_info->tpz;
2474  brNumiFluxTptype = gnumi_flux_info->tptype;
2475  brNumiFluxTgen = gnumi_flux_info->tgen;
2476  brNumiFluxTgptype = gnumi_flux_info->tgptype;
2477  brNumiFluxTgppx = gnumi_flux_info->tgppx;
2478  brNumiFluxTgppy = gnumi_flux_info->tgppy;
2479  brNumiFluxTgppz = gnumi_flux_info->tgppz;
2480  brNumiFluxTprivx = gnumi_flux_info->tprivx;
2481  brNumiFluxTprivy = gnumi_flux_info->tprivy;
2482  brNumiFluxTprivz = gnumi_flux_info->tprivz;
2483  brNumiFluxBeamx = gnumi_flux_info->beamx;
2484  brNumiFluxBeamy = gnumi_flux_info->beamy;
2485  brNumiFluxBeamz = gnumi_flux_info->beamz;
2486  brNumiFluxBeampx = gnumi_flux_info->beampx;
2487  brNumiFluxBeampy = gnumi_flux_info->beampy;
2488  brNumiFluxBeampz = gnumi_flux_info->beampz;
2489  } // gnumi_flux_info
2490 #endif
2491  } // kConvFmt_numi_rootracker
2492 
2493  // fill tree
2494  rootracker_tree->Fill();
2495  mcrec->Clear();
2496 
2497  } // event loop
2498 
2499  // Copy POT normalization for the generated sample
2500  double pot = gtree->GetWeight();
2501  rootracker_tree->SetWeight(pot);
2502 
2503  // Copy MC job metadata (gconfig and genv TFolders)
2504  if(gOptCopyJobMeta) {
2505  TFolder * genv = (TFolder*) fin.Get("genv");
2506  TFolder * gconfig = (TFolder*) fin.Get("gconfig");
2507  fout.cd();
2508  genv -> Write("genv");
2509  gconfig -> Write("gconfig");
2510  }
2511 
2512  fin.Close();
2513 
2514  fout.Write();
2515  fout.Close();
2516 
2517  LOG("gntpc", pINFO) << "\nDone converting GENIE's GHEP ntuple";
2518 }
2519 //____________________________________________________________________________________
2520 // GENIE GHEP EVENT TREE -> NEUGEN-style format for AGKY studies
2521 //____________________________________________________________________________________
2522 void ConvertToGHad(void)
2523 {
2524 // Neugen-style text format for the AGKY hadronization model studies
2525 // Format:
2526 // (blank line)
2527 // event number, neutrino particle code, CCNC, IM, A, Z
2528 // int_type, x, y, w, ihadmod
2529 // neutrino particle code, 5 vec
2530 // lepton particle code, 5-vec
2531 // outgoing hadronic system, 5-vec
2532 // number of stable daughters of hadronic system
2533 // ... then for each stable daughter
2534 // particle id, 5 vec
2535 
2536  //-- open the ROOT file and get the TTree & its header
2537  TFile fin(gOptInpFileName.c_str(),"READ");
2538  TTree * tree = 0;
2539  NtpMCTreeHeader * thdr = 0;
2540  tree = dynamic_cast <TTree *> ( fin.Get("gtree") );
2541  thdr = dynamic_cast <NtpMCTreeHeader *> ( fin.Get("header") );
2542 
2543  LOG("gntpc", pINFO) << "Input tree header: " << *thdr;
2544 
2545  //-- get mc record
2546  NtpMCEventRecord * mcrec = 0;
2547  tree->SetBranchAddress("gmcrec", &mcrec);
2548 
2549  //-- open the output stream
2550  ofstream output(gOptOutFileName.c_str(), ios::out);
2551 
2552  //-- open output root file and create ntuple -- if required
2553 #ifdef __GHAD_NTP__
2554  TFile fout("ghad.root","recreate");
2555  TTree * ghad = new TTree("ghad","");
2556  ghad->Branch("i", &brIev, "i/I" );
2557  ghad->Branch("W", &brW, "W/D" );
2558  ghad->Branch("n", &brN, "n/I" );
2559  ghad->Branch("pdg", brPdg, "pdg[n]/I" );
2560  ghad->Branch("E", brE, "E[n]/D" );
2561  ghad->Branch("px", brPx, "px[n]/D" );
2562  ghad->Branch("py", brPy, "py[n]/D" );
2563  ghad->Branch("pz", brPz, "pz[n]/D" );
2564 #endif
2565 
2566  //-- figure out how many events to analyze
2567  Long64_t nmax = (gOptN<0) ?
2568  tree->GetEntries() : TMath::Min(tree->GetEntries(), gOptN);
2569  if (nmax<0) {
2570  LOG("gntpc", pERROR) << "Number of events = 0";
2571  return;
2572  }
2573  LOG("gntpc", pNOTICE) << "*** Analyzing: " << nmax << " events";
2574 
2575  //-- event loop
2576  for(Long64_t iev = 0; iev < nmax; iev++) {
2577  tree->GetEntry(iev);
2578  NtpMCRecHeader rec_header = mcrec->hdr;
2579  EventRecord & event = *(mcrec->event);
2580 
2581  LOG("gntpc", pINFO) << rec_header;
2582  LOG("gntpc", pINFO) << event;
2583 
2584 #ifdef __GHAD_NTP__
2585  brN = 0;
2586  for(int k=0; k<kNPmax; k++) {
2587  brPdg[k]=0;
2588  brE [k]=0;
2589  brPx [k]=0;
2590  brPy [k]=0;
2591  brPz [k]=0;
2592  }
2593 #endif
2594 
2595  //
2596  // convert the current event
2597  //
2598  const Interaction * interaction = event.Summary();
2599  const ProcessInfo & proc_info = interaction->ProcInfo();
2600  const InitialState & init_state = interaction->InitState();
2601 
2602  bool is_dis = proc_info.IsDeepInelastic();
2603  bool is_res = proc_info.IsResonant();
2604  bool is_cc = proc_info.IsWeakCC();
2605 
2606  bool pass = is_cc && (is_dis || is_res);
2607  if(!pass) {
2608  mcrec->Clear();
2609  continue;
2610  }
2611 
2612  int ccnc = is_cc ? 1 : 0;
2613  int inttyp = 3;
2614 
2615  int im = -1;
2616  if (init_state.IsNuP ()) im = 1;
2617  else if (init_state.IsNuN ()) im = 2;
2618  else if (init_state.IsNuBarP ()) im = 3;
2619  else if (init_state.IsNuBarN ()) im = 4;
2620  else return;
2621 
2622  GHepParticle * neutrino = event.Probe();
2623  assert(neutrino);
2624  GHepParticle * target = event.Particle(1);
2625  assert(target);
2626  GHepParticle * fsl = event.FinalStatePrimaryLepton();
2627  assert(fsl);
2628  GHepParticle * hitnucl = event.HitNucleon();
2629  assert(hitnucl);
2630 
2631  int nupdg = neutrino->Pdg();
2632  int fslpdg = fsl->Pdg();
2633  int A = target->A();
2634  int Z = target->Z();
2635 
2636  const TLorentzVector & k1 = *(neutrino->P4()); // v 4-p (k1)
2637  const TLorentzVector & k2 = *(fsl->P4()); // l 4-p (k2)
2638 // const TLorentzVector & p1 = *(hitnucl->P4()); // N 4-p (p1)
2639 // const TLorentzVector & ph = *(hadsyst->P4()); // had-syst 4-p
2640 
2641  TLorentzVector ph;
2642  if(is_dis) {
2643  GHepParticle * hadsyst = event.FinalStateHadronicSystem();
2644  assert(hadsyst);
2645  ph = *(hadsyst->P4());
2646  }
2647  if(is_res) {
2648  GHepParticle * hadres = event.Particle(hitnucl->FirstDaughter());
2649  ph = *(hadres->P4());
2650  }
2651 
2652  const Kinematics & kine = interaction->Kine();
2653  bool get_selected = true;
2654  double x = kine.x (get_selected);
2655  double y = kine.y (get_selected);
2656  double W = kine.W (get_selected);
2657 
2658  int hadmod = -1;
2659  int ihadmom = -1;
2660  TIter event_iter(&event);
2661  GHepParticle * p = 0;
2662  int i=-1;
2663  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
2664  i++;
2665  int pdg = p->Pdg();
2666  if (pdg == kPdgHadronicSyst ) { hadmod= 2; ihadmom=i; }
2667  if (pdg == kPdgString ) { hadmod=11; ihadmom=i; }
2668  if (pdg == kPdgCluster ) { hadmod=12; ihadmom=i; }
2669  if (pdg == kPdgIndep ) { hadmod=13; ihadmom=i; }
2670  }
2671 
2672  output << endl;
2673  output << iev << "\t"
2674  << nupdg << "\t" << ccnc << "\t" << im << "\t"
2675  << A << "\t" << Z << endl;
2676  output << inttyp << "\t" << x << "\t" << y << "\t" << W << "\t"
2677  << hadmod << endl;
2678  output << nupdg << "\t"
2679  << k1.Px() << "\t" << k1.Py() << "\t" << k1.Pz() << "\t"
2680  << k1.Energy() << "\t" << k1.M() << endl;
2681  output << fslpdg << "\t"
2682  << k2.Px() << "\t" << k2.Py() << "\t" << k2.Pz() << "\t"
2683  << k2.Energy() << "\t" << k2.M() << endl;
2684  output << 111111 << "\t"
2685  << ph.Px() << "\t" << ph.Py() << "\t" << ph.Pz() << "\t"
2686  << ph.Energy() << "\t" << ph.M() << endl;
2687 
2688  vector<int> hadv;
2689 
2690  event_iter.Reset();
2691  i=-1;
2692  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
2693  i++;
2694  if(i<ihadmom) continue;
2695 
2696  GHepStatus_t ist = p->Status();
2697  int pdg = p->Pdg();
2698 
2699  if(ist == kIStDISPreFragmHadronicState) continue;
2700 
2701  if(ist == kIStStableFinalState) {
2702  GHepParticle * mom = event.Particle(p->FirstMother());
2703  GHepStatus_t mom_ist = mom->Status();
2704  int mom_pdg = mom->Pdg();
2705  bool skip = (mom_pdg == kPdgPi0 && mom_ist== kIStDecayedState);
2706  if(!skip) { hadv.push_back(i); }
2707  }
2708 
2709  if(pdg==kPdgPi0 && ist==kIStDecayedState) { hadv.push_back(i); }
2710  }
2711 
2712  output << hadv.size() << endl;
2713 
2714 #ifdef __GHAD_NTP__
2715  brIev = (int) iev;
2716  brW = W;
2717  brN = hadv.size();
2718  int k=0;
2719 #endif
2720 
2721  vector<int>::const_iterator hiter = hadv.begin();
2722  for( ; hiter != hadv.end(); ++hiter) {
2723  int id = *hiter;
2724  GHepParticle * particle = event.Particle(id);
2725  int pdg = particle->Pdg();
2726  double px = particle->P4()->Px();
2727  double py = particle->P4()->Py();
2728  double pz = particle->P4()->Pz();
2729  double E = particle->P4()->Energy();
2730  double m = particle->P4()->M();
2731  output << pdg << "\t"
2732  << px << "\t" << py << "\t" << pz << "\t"
2733  << E << "\t" << m << endl;
2734 
2735 #ifdef __GHAD_NTP__
2736  brPx[k] = px;
2737  brPy[k] = py;
2738  brPz[k] = pz;
2739  brE[k] = E;
2740  brPdg[k] = pdg;
2741  k++;
2742 #endif
2743  }
2744 
2745 #ifdef __GHAD_NTP__
2746  ghad->Fill();
2747 #endif
2748 
2749  mcrec->Clear();
2750 
2751  } // event loop
2752 
2753  output.close();
2754  fin.Close();
2755 
2756 #ifdef __GHAD_NTP__
2757  ghad->Write("ghad");
2758  fout.Write();
2759  fout.Close();
2760 #endif
2761 
2762  LOG("gntpc", pINFO) << "\nDone converting GENIE's GHEP ntuple";
2763 }
2764 //____________________________________________________________________________________
2765 // GENIE GHEP EVENT TREE -> Summary tree for INTRANUKE studies
2766 //____________________________________________________________________________________
2767 void ConvertToGINuke(void)
2768 {
2769  //-- output tree branch variables
2770  //
2771  int brIEv = 0; // Event number
2772  int brProbe = 0; // Incident hadron code
2773  int brTarget = 0; // Nuclear target pdg code (10LZZZAAAI)
2774  double brKE = 0; // Probe kinetic energy
2775  double brE = 0; // Probe energy
2776  double brP = 0; // Probe momentum
2777  int brTgtA = 0; // Target A (mass number)
2778  int brTgtZ = 0; // Target Z (atomic number)
2779  double brVtxX = 0; // "Vertex x" (initial placement of h /in h+A events/ on the nuclear boundary)
2780  double brVtxY = 0; // "Vertex y"
2781  double brVtxZ = 0; // "Vertex z"
2782  int brProbeFSI = 0; // Rescattering code for incident hadron
2783  double brDist = 0; // Distance travelled by h before interacting (if at all before escaping)
2784  int brNh = 0; // Number of final state hadrons
2785  int brPdgh [kNPmax]; // Pdg code of i^th final state hadron
2786  double brEh [kNPmax]; // Energy of i^th final state hadron
2787  double brPh [kNPmax]; // P of i^th final state hadron
2788  double brPxh [kNPmax]; // Px of i^th final state hadron
2789  double brPyh [kNPmax]; // Py of i^th final state hadron
2790  double brPzh [kNPmax]; // Pz of i^th final state hadron
2791  double brCosth [kNPmax]; // Cos(th) of i^th final state hadron
2792  double brMh [kNPmax]; // Mass of i^th final state hadron
2793  int brNp = 0; // Number of final state p
2794  int brNn = 0; // Number of final state n
2795  int brNpip = 0; // Number of final state pi+
2796  int brNpim = 0; // Number of final state pi-
2797  int brNpi0 = 0; // Number of final state pi0
2798 
2799  //-- open output file & create output summary tree & create the tree branches
2800  //
2801  LOG("gntpc", pNOTICE)
2802  << "*** Saving summary tree to: " << gOptOutFileName;
2803  TFile fout(gOptOutFileName.c_str(),"recreate");
2804 
2805 TTree * tEvtTree = new TTree("ginuke","GENIE INuke Summary Tree");
2806  assert(tEvtTree);
2807 
2808  //-- create tree branches
2809  //
2810  tEvtTree->Branch("iev", &brIEv, "iev/I" );
2811  tEvtTree->Branch("probe", &brProbe, "probe/I" );
2812  tEvtTree->Branch("tgt" , &brTarget, "tgt/I" );
2813  tEvtTree->Branch("ke", &brKE, "ke/D" );
2814  tEvtTree->Branch("e", &brE, "e/D" );
2815  tEvtTree->Branch("p", &brP, "p/D" );
2816  tEvtTree->Branch("A", &brTgtA, "A/I" );
2817  tEvtTree->Branch("Z", &brTgtZ, "Z/I" );
2818  tEvtTree->Branch("vtxx", &brVtxX, "vtxx/D" );
2819  tEvtTree->Branch("vtxy", &brVtxY, "vtxy/D" );
2820  tEvtTree->Branch("vtxz", &brVtxZ, "vtxz/D" );
2821  tEvtTree->Branch("probe_fsi", &brProbeFSI, "probe_fsi/I" );
2822  tEvtTree->Branch("dist", &brDist, "dist/D" );
2823  tEvtTree->Branch("nh", &brNh, "nh/I" );
2824  tEvtTree->Branch("pdgh", brPdgh, "pdgh[nh]/I" );
2825  tEvtTree->Branch("Eh", brEh, "Eh[nh]/D" );
2826  tEvtTree->Branch("ph", brPh, "ph[nh]/D" );
2827  tEvtTree->Branch("pxh", brPxh, "pxh[nh]/D" );
2828  tEvtTree->Branch("pyh", brPyh, "pyh[nh]/D" );
2829  tEvtTree->Branch("pzh", brPzh, "pzh[nh]/D" );
2830  tEvtTree->Branch("cth", brCosth, "cth[nh]/D" );
2831  tEvtTree->Branch("mh", brMh, "mh[nh]/D" );
2832  tEvtTree->Branch("np", &brNp, "np/I" );
2833  tEvtTree->Branch("nn", &brNn, "nn/I" );
2834  tEvtTree->Branch("npip", &brNpip, "npip/I" );
2835  tEvtTree->Branch("npim", &brNpim, "npim/I" );
2836  tEvtTree->Branch("npi0", &brNpi0, "npi0/I" );
2837 
2838  //-- open the ROOT file and get the TTree & its header
2839  TFile fin(gOptInpFileName.c_str(),"READ");
2840  TTree * er_tree = 0;
2841  NtpMCTreeHeader * thdr = 0;
2842  er_tree = dynamic_cast <TTree *> ( fin.Get("gtree") );
2843  thdr = dynamic_cast <NtpMCTreeHeader *> ( fin.Get("header") );
2844  if (!er_tree) {
2845  LOG("gntpc", pERROR) << "Null input tree";
2846  return;
2847  }
2848  LOG("gntpc", pINFO) << "Input tree header: " << *thdr;
2849 
2850  //-- get the mc record
2851  NtpMCEventRecord * mcrec = 0;
2852  er_tree->SetBranchAddress("gmcrec", &mcrec);
2853  if (!mcrec) {
2854  LOG("gntpc", pERROR) << "Null MC record";
2855  return;
2856  }
2857 
2858  //-- figure out how many events to analyze
2859  Long64_t nmax = (gOptN<0) ?
2860  er_tree->GetEntries() : TMath::Min( er_tree->GetEntries(), gOptN );
2861  if (nmax<0) {
2862  LOG("gntpc", pERROR) << "Number of events = 0";
2863  return;
2864  }
2865  LOG("gntpc", pNOTICE) << "*** Analyzing: " << nmax << " events";
2866 
2867  for(Long64_t iev = 0; iev < nmax; iev++) {
2868  brIEv = iev;
2869  er_tree->GetEntry(iev);
2870  NtpMCRecHeader rec_header = mcrec->hdr;
2871  EventRecord & event = *(mcrec->event);
2872 
2873  LOG("gntpc", pINFO) << rec_header;
2874  LOG("gntpc", pINFO) << event;
2875 
2876  // analyze current event and fill the summary ntuple
2877 
2878  // clean-up arrays
2879  //
2880  for(int j=0; j<kNPmax; j++) {
2881  brPdgh[j] = 0;
2882  brEh [j] = 0;
2883  brPxh [j] = 0;
2884  brPyh [j] = 0;
2885  brPzh [j] = 0;
2886  brMh [j] = 0;
2887  }
2888 
2889  //
2890  // convert the current event
2891  //
2892 
2893  GHepParticle * probe = event.Particle(0);
2894  GHepParticle * target = event.Particle(1);
2895  assert(probe && target);
2896 
2897  brProbe = probe -> Pdg();
2898  brTarget = target -> Pdg();
2899  brKE = probe -> KinE();
2900  brE = probe -> E();
2901  brP = probe -> P4()->Vect().Mag();
2902  brTgtA = pdg::IonPdgCodeToA(target->Pdg());
2903  brTgtZ = pdg::IonPdgCodeToZ(target->Pdg());
2904  brVtxX = probe -> Vx();
2905  brVtxY = probe -> Vy();
2906  brVtxZ = probe -> Vz();
2907  brProbeFSI = probe -> RescatterCode();
2908  GHepParticle * rescattered_hadron = event.Particle(probe->FirstDaughter());
2909  assert(rescattered_hadron);
2910  if(rescattered_hadron->Status() == kIStStableFinalState) {
2911  brDist = -1; // hadron escaped nucleus before interacting;
2912  }
2913  else {
2914  double x = rescattered_hadron->Vx();
2915  double y = rescattered_hadron->Vy();
2916  double z = rescattered_hadron->Vz();
2917  double d2 = TMath::Power(brVtxX-x,2) +
2918  TMath::Power(brVtxY-y,2) +
2919  TMath::Power(brVtxZ-z,2);
2920  brDist = TMath::Sqrt(d2);
2921  }
2922 
2923  brNp = 0;
2924  brNn = 0;
2925  brNpip = 0;
2926  brNpim = 0;
2927  brNpi0 = 0;
2928 
2929  int i=0;
2930  GHepParticle * p = 0;
2931  TIter event_iter(&event);
2932  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
2933  if(pdg::IsPseudoParticle(p->Pdg())) continue;
2934  if(p->Status() != kIStStableFinalState) continue;
2935 
2936  brPdgh[i] = p->Pdg();
2937  brEh [i] = p->E();
2938  brPxh [i] = p->Px();
2939  brPyh [i] = p->Py();
2940  brPzh [i] = p->Pz();
2941  brPh [i] =
2942  TMath::Sqrt(brPxh[i]*brPxh[i]+brPyh[i]*brPyh[i]
2943  +brPzh[i]*brPzh[i]);
2944  brCosth[i] = brPzh[i]/brPh[i];
2945  brMh [i] = p->Mass();
2946 
2947  if ( p->Pdg() == kPdgProton ) brNp++;
2948  if ( p->Pdg() == kPdgNeutron ) brNn++;
2949  if ( p->Pdg() == kPdgPiP ) brNpip++;
2950  if ( p->Pdg() == kPdgPiM ) brNpim++;
2951  if ( p->Pdg() == kPdgPi0 ) brNpi0++;
2952 
2953  i++;
2954  }
2955  brNh = i;
2956 
2957  ///////////////Test Code///////////////////////
2958  int tempProbeFSI = brProbeFSI;
2959  brProbeFSI = HAProbeFSI(tempProbeFSI, brProbe, brNh, brEh, brPdgh, brNpip, brNpim, brNpi0);
2960  //////////////End Test/////////////////////////
2961 
2962 
2963  // fill the summary tree
2964  tEvtTree->Fill();
2965 
2966  mcrec->Clear();
2967 
2968  } // event loop
2969 
2970  fin.Close();
2971 
2972  fout.Write();
2973  fout.Close();
2974 
2975  LOG("gntpc", pINFO) << "\nDone converting GENIE's GHEP ntuple";
2976 }
2977 //____________________________________________________________________________________
2978 // FUNCTIONS FOR PARSING CMD-LINE ARGUMENTS
2979 //____________________________________________________________________________________
2980 void GetCommandLineArgs(int argc, char ** argv)
2981 {
2982  // Common run options.
2983  RunOpt::Instance()->ReadFromCommandLine(argc,argv);
2984 
2985  // Parse run options for this app
2986 
2987  CmdLnArgParser parser(argc,argv);
2988 
2989  // get input ROOT file (containing a GENIE GHEP event tree)
2990  if( parser.OptionExists('i') ) {
2991  LOG("gntpc", pINFO) << "Reading input filename";
2992  gOptInpFileName = parser.ArgAsString('i');
2993  } else {
2994  LOG("gntpc", pFATAL)
2995  << "Unspecified input filename - Exiting";
2996  PrintSyntax();
2997  gAbortingInErr = true;
2998  exit(1);
2999  }
3000 
3001  // check input GENIE ROOT file
3002  bool inpok = !(gSystem->AccessPathName(gOptInpFileName.c_str()));
3003  if (!inpok) {
3004  LOG("gntpc", pFATAL)
3005  << "The input ROOT file ["
3006  << gOptInpFileName << "] is not accessible";
3007  gAbortingInErr = true;
3008  exit(2);
3009  }
3010 
3011  // get output file format
3012  if( parser.OptionExists('f') ) {
3013  LOG("gntpc", pINFO) << "Reading output file format";
3014  string fmt = parser.ArgAsString('f');
3015 
3016  if (fmt == "gst") { gOptOutFileFormat = kConvFmt_gst; }
3017  else if (fmt == "gxml") { gOptOutFileFormat = kConvFmt_gxml; }
3018  else if (fmt == "ghep_mock_data") { gOptOutFileFormat = kConvFmt_ghep_mock_data; }
3019  else if (fmt == "rootracker") { gOptOutFileFormat = kConvFmt_rootracker; }
3020  else if (fmt == "rootracker_mock_data") { gOptOutFileFormat = kConvFmt_rootracker_mock_data; }
3021  else if (fmt == "t2k_rootracker") { gOptOutFileFormat = kConvFmt_t2k_rootracker; }
3022  else if (fmt == "numi_rootracker") { gOptOutFileFormat = kConvFmt_numi_rootracker; }
3023  else if (fmt == "t2k_tracker") { gOptOutFileFormat = kConvFmt_t2k_tracker; }
3024  else if (fmt == "nuance_tracker" ) { gOptOutFileFormat = kConvFmt_nuance_tracker; }
3025  else if (fmt == "ghad") { gOptOutFileFormat = kConvFmt_ghad; }
3026  else if (fmt == "ginuke") { gOptOutFileFormat = kConvFmt_ginuke; }
3027  else { gOptOutFileFormat = kConvFmt_undef; }
3028 
3029  if(gOptOutFileFormat == kConvFmt_undef) {
3030  LOG("gntpc", pFATAL) << "Unknown output file format (" << fmt << ")";
3031  gAbortingInErr = true;
3032  exit(3);
3033  }
3034 
3035  } else {
3036  LOG("gntpc", pFATAL) << "Unspecified output file format";
3037  gAbortingInErr = true;
3038  exit(4);
3039  }
3040 
3041  // get output file name
3042  if( parser.OptionExists('o') ) {
3043  LOG("gntpc", pINFO) << "Reading output filename";
3044  gOptOutFileName = parser.ArgAsString('o');
3045  } else {
3046  LOG("gntpc", pINFO)
3047  << "Unspecified output filename - Using default";
3048  gOptOutFileName = DefaultOutputFile();
3049  }
3050 
3051  // get number of events to convert
3052  if( parser.OptionExists('n') ) {
3053  LOG("gntpc", pINFO) << "Reading number of events to analyze";
3054  gOptN = parser.ArgAsInt('n');
3055  } else {
3056  LOG("gntpc", pINFO)
3057  << "Unspecified number of events to analyze - Use all";
3058  gOptN = -1;
3059  }
3060 
3061  // get format version number
3062  if( parser.OptionExists('v') ) {
3063  LOG("gntpc", pINFO) << "Reading format version number";
3064  gOptVersion = parser.ArgAsInt('v');
3065  LOG("gntpc", pINFO)
3066  << "Using version number: " << gOptVersion;
3067  } else {
3068  LOG("gntpc", pINFO)
3069  << "Unspecified version number - Use latest";
3070  gOptVersion = LatestFormatVersionNumber();
3071  LOG("gntpc", pINFO)
3072  << "Latest version number: " << gOptVersion;
3073  }
3074 
3075  // check whether to copy MC job metadata (only if output file is in ROOT format)
3076  gOptCopyJobMeta = parser.OptionExists('c');
3077 
3078  // random number seed
3079  if( parser.OptionExists("seed") ) {
3080  LOG("gntpc", pINFO) << "Reading random number seed";
3081  gOptRanSeed = parser.ArgAsLong("seed");
3082  } else {
3083  LOG("gntpc", pINFO) << "Unspecified random number seed - Using default";
3084  gOptRanSeed = -1;
3085  }
3086 
3087  LOG("gntpc", pNOTICE) << "Input filename = " << gOptInpFileName;
3088  LOG("gntpc", pNOTICE) << "Output filename = " << gOptOutFileName;
3089  LOG("gntpc", pNOTICE) << "Conversion to format = " << gOptRanSeed
3090  << ", vrs = " << gOptVersion;
3091  LOG("gntpc", pNOTICE) << "Number of events to be converted = " << gOptN;
3092  LOG("gntpc", pNOTICE) << "Copy metadata? = " << ((gOptCopyJobMeta) ? "Yes" : "No");
3093  LOG("gntpc", pNOTICE) << "Random number seed = " << gOptRanSeed;
3094 
3095  LOG("gntpc", pNOTICE) << *RunOpt::Instance();
3096 }
3097 //____________________________________________________________________________________
3098 string DefaultOutputFile(void)
3099 {
3100  // filename extension - depending on file format
3101  string ext="";
3102  if (gOptOutFileFormat == kConvFmt_gst ) { ext = "gst.root"; }
3103  else if (gOptOutFileFormat == kConvFmt_gxml ) { ext = "gxml"; }
3104  else if (gOptOutFileFormat == kConvFmt_ghep_mock_data ) { ext = "mockd.ghep.root"; }
3105  else if (gOptOutFileFormat == kConvFmt_rootracker ) { ext = "gtrac.root"; }
3106  else if (gOptOutFileFormat == kConvFmt_rootracker_mock_data ) { ext = "mockd.gtrac.root"; }
3107  else if (gOptOutFileFormat == kConvFmt_t2k_rootracker ) { ext = "gtrac.root"; }
3108  else if (gOptOutFileFormat == kConvFmt_numi_rootracker ) { ext = "gtrac.root"; }
3109  else if (gOptOutFileFormat == kConvFmt_t2k_tracker ) { ext = "gtrac.dat"; }
3110  else if (gOptOutFileFormat == kConvFmt_nuance_tracker ) { ext = "gtrac_legacy.dat"; }
3111  else if (gOptOutFileFormat == kConvFmt_ghad ) { ext = "ghad.dat"; }
3112  else if (gOptOutFileFormat == kConvFmt_ginuke ) { ext = "ginuke.root"; }
3113 
3114  string inpname = gOptInpFileName;
3115  unsigned int L = inpname.length();
3116 
3117  // if the last 4 characters are "root" (ROOT file extension) then
3118  // remove them
3119  if(inpname.substr(L-4, L).find("root") != string::npos) {
3120  inpname.erase(L-4, L);
3121  }
3122 
3123  // remove ghep.
3124  size_t pos = inpname.find("ghep.");
3125  if(pos != string::npos) {
3126  inpname.erase(pos, pos+4);
3127  }
3128 
3129  ostringstream name;
3130  name << inpname << ext;
3131 
3132  return gSystem->BaseName(name.str().c_str());
3133 }
3134 //____________________________________________________________________________________
3135 int LatestFormatVersionNumber(void)
3136 {
3137  if (gOptOutFileFormat == kConvFmt_gst ) return 1;
3138  else if (gOptOutFileFormat == kConvFmt_gxml ) return 1;
3139  else if (gOptOutFileFormat == kConvFmt_ghep_mock_data ) return 1;
3140  else if (gOptOutFileFormat == kConvFmt_rootracker ) return 1;
3141  else if (gOptOutFileFormat == kConvFmt_rootracker_mock_data ) return 1;
3142  else if (gOptOutFileFormat == kConvFmt_t2k_rootracker ) return 1;
3143  else if (gOptOutFileFormat == kConvFmt_numi_rootracker ) return 1;
3144  else if (gOptOutFileFormat == kConvFmt_t2k_tracker ) return 2;
3145  else if (gOptOutFileFormat == kConvFmt_nuance_tracker ) return 1;
3146  else if (gOptOutFileFormat == kConvFmt_ghad ) return 1;
3147  else if (gOptOutFileFormat == kConvFmt_ginuke ) return 1;
3148 
3149  return -1;
3150 }
3151 //____________________________________________________________________________________
3152 void PrintSyntax(void)
3153 {
3154  string basedir = string( gSystem->Getenv("GENIE") );
3155  string thisfile = basedir + string("/src/Apps/gNtpConv.cxx");
3156  string cmd = "less " + thisfile;
3157 
3158  gSystem->Exec(cmd.c_str());
3159 }
3160 //____________________________________________________________________________________
3161 /* Converting HN probe_fsi to HA probe_fsi */
3162 int HAProbeFSI(int probe_fsi, int probe_pdg, int numh, double E_had[], int pdg_had[], int numpip, int numpim, int numpi0)
3163 {
3164  int index = -1;
3165  double energy = 0;
3166 
3167  for(int i=0; i<numh; i++)
3168  { energy += E_had[i]; }
3169 
3170 
3171 // Determine fates (as defined in Intranuke/INukeUtils.cxx/ utils::intranuke::FindhAFate())
3172  if (probe_fsi==3 && numh==1) // Elastic
3173  { index=3; }
3174  else if (energy==E_had[0] && numh==1) // No interaction
3175  { index=1; }
3176  else if ( pdg::IsPion(probe_pdg) && numpip+numpi0+numpim==0) // Absorption
3177  { index=5; }
3178  else if ( (pdg::IsNucleon(probe_pdg) && numpip+numpi0+numpim==0 && numh>2 )
3179  || (probe_pdg==kPdgGamma && energy!=E_had[0] && numpip+numpi0+numpim==0)) // Knock-out
3180  { index=6; }
3181  else if ( numpip+numpi0+numpim > (pdg::IsPion(probe_pdg) ? 1 : 0) ) // Pion production
3182  { index=7; }
3183  else if ( numh>=2 ) // Inelastic or Charge Exchange
3184  {
3185  for(int i = 0; i < numh; i++)
3186  {
3187  if ( (pdg::IsPion(probe_pdg) && ( probe_pdg==pdg_had[i] ))
3188  || pdg::IsNucleon(probe_pdg) )
3189  {index=4;}
3190  if(index!=4)
3191  {index=2;}
3192  }
3193  }
3194  else //Double Charge Exchange or Undefined
3195  {
3196  bool undef = true;
3197  if ( pdg::IsPion(probe_pdg) )
3198  {
3199  for (int iter = 0; iter < numh; iter++)
3200  {
3201  if (probe_pdg==211 && pdg_had[iter]==-211) { index=8; undef=false; }
3202  else if (probe_pdg==-211 && pdg_had[iter]==211) { index=8; undef=false; }
3203  }
3204  }
3205  if (undef) { index=0; }
3206  }
3207 
3208  return index;
3209 }
name
static QCString name
Definition: declinfo.cpp:673
genie::ProcessInfo::IsResonant
bool IsResonant(void) const
Definition: ProcessInfo.cxx:94
genie::utils::app_init::RandGen
void RandGen(long int seed)
Definition: AppInit.cxx:30
genie::GHepRecord::SetPrintLevel
static void SetPrintLevel(int print_level)
Definition: GHepRecord.cxx:948
ConvertToGHad
void ConvertToGHad(void)
Definition: gNtpConv.cxx:2522
main
int main(int argc, char **argv)
Definition: gNtpConv.cxx:230
genie::utils::ghep::NeutReactionCode
int NeutReactionCode(const GHepRecord *evrec)
Definition: GHepUtils.cxx:22
genie::Kinematics::W
double W(bool selected=false) const
Definition: Kinematics.cxx:157
genie::GHepParticle::Z
int Z(void) const
Definition: GHepParticle.cxx:174
genie::pdg::IsPion
bool IsPion(int pdgc)
Definition: PDGUtils.cxx:323
genie::CmdLnArgParser::ArgAsLong
long ArgAsLong(char opt)
Definition: CmdLnArgParser.cxx:165
genie::constants
Basic constants.
genie::GHepParticle::RescatterCode
int RescatterCode(void) const
Definition: GHepParticle.h:65
reco_momentum_tuples.py
py
Definition: reco_momentum_tuples.py:219
genie::utils::system::GenieMajorVrsNum
int GenieMajorVrsNum(string tag)
Definition: SystemUtils.cxx:59
genie::pdg::IsParticle
bool IsParticle(int pdgc)
not ion or pseudo-particle
Definition: PDGUtils.cxx:44
genie::Target::HitSeaQrk
bool HitSeaQrk(void) const
Definition: Target.cxx:299
genie::ProcessInfo::IsWeakCC
bool IsWeakCC(void) const
Definition: ProcessInfo.cxx:183
genie::InitialState::IsNuBarN
bool IsNuBarN(void) const
is anti-neutrino + neutron?
Definition: InitialState.cxx:226
generate_CCQE_events.energy
string energy
Definition: generate_CCQE_events.py:15
make_norm_csv.output
list output
Definition: make_norm_csv.py:48
genie::NtpMCRecordI::hdr
NtpMCRecHeader hdr
record header
Definition: NtpMCRecordI.h:38
genie::CmdLnArgParser::ArgAsString
string ArgAsString(char opt)
Definition: CmdLnArgParser.cxx:106
genie
THE MAIN GENIE PROJECT NAMESPACE
Definition: AlgCmp.h:25
pERROR
#define pERROR
Definition: Messenger.h:59
genie::GHepParticle::E
double E(void) const
Get energy.
Definition: GHepParticle.h:91
genie::PDGLibrary::AddDarkMatter
void AddDarkMatter(double mass, double med_ratio)
Definition: PDGLibrary.cxx:142
genie::constants::kNucleonMass
static const double kNucleonMass
Definition: Constants.h:77
genie::utils::kinematics::Q2
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
genie::NtpWriter::CustomizeFilename
void CustomizeFilename(string filename)
Definition: NtpWriter.cxx:128
genie::GHepParticle::P4
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
string
std::string string
Definition: nybbler.cc:12
genie::GHepParticle::FirstDaughter
int FirstDaughter(void) const
Definition: GHepParticle.h:68
genie::pdg::IsNucleon
bool IsNucleon(int pdgc)
Definition: PDGUtils.cxx:343
ConvertToGXML
void ConvertToGXML(void)
Definition: gNtpConv.cxx:1100
genie::NtpMCEventRecord
MINOS-style ntuple record. Each such ntuple record holds a generated EventRecord object. Ntuples of this type are intended for feeding GENIE events into other applications (for example the GEANT4 based MC generation framework of an experiment) if no direct interface exists.
Definition: NtpMCEventRecord.h:37
genie::GHepParticle::PolzPolarAngle
double PolzPolarAngle(void) const
Definition: GHepParticle.h:119
genie::Target::HitQrkPdg
int HitQrkPdg(void) const
Definition: Target.cxx:242
genie::ProcessInfo::IsInverseMuDecay
bool IsInverseMuDecay(void) const
Definition: ProcessInfo.cxx:121
genie::RunOpt::ReadFromCommandLine
void ReadFromCommandLine(int argc, char **argv)
Definition: RunOpt.cxx:99
genie::ProcessInfo::IsQuasiElastic
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
genie::InitialState::IsNuN
bool IsNuN(void) const
is neutrino + neutron?
Definition: InitialState.cxx:208
genie::utils::ghep::NuanceReactionCode
int NuanceReactionCode(const GHepRecord *evrec)
Definition: GHepUtils.cxx:284
pFATAL
#define pFATAL
Definition: Messenger.h:56
tracks.hpdg
hpdg
Definition: tracks.py:120
gOptCopyJobMeta
bool gOptCopyJobMeta
copy MC job metadata (gconfig, genv TFolders)
Definition: gNtpConv.cxx:219
genie::pdg::IonPdgCodeToA
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:60
genie::Kinematics
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
vector
struct vector vector
genie::units::fs
static constexpr double fs
Definition: Units.h:100
genie::Kinematics::x
double x(bool selected=false) const
Definition: Kinematics.cxx:99
kNPmax
const int kNPmax
Definition: gNtpConv.cxx:228
gar::P
std::pair< float, std::string > P
Definition: anatree_module.cc:75
ConvertToGST
void ConvertToGST(void)
Definition: gNtpConv.cxx:295
run_cnn_1class.tracks
tracks
Definition: run_cnn_1class.py:54
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
genie::ProcessInfo::IsDiffractive
bool IsDiffractive(void) const
Definition: ProcessInfo.cxx:161
genie::ProcessInfo::IsCoherentProduction
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:99
genie::GHepParticle::Mass
double Mass(void) const
Mass that corresponds to the PDG code.
Definition: GHepParticle.cxx:143
genie::ProcessInfo::IsIMDAnnihilation
bool IsIMDAnnihilation(void) const
Definition: ProcessInfo.cxx:126
reco_momentum_tuples.px
px
Definition: reco_momentum_tuples.py:218
genie::kPdgElectron
const int kPdgElectron
Definition: PDGCodes.h:35
genie::units::MeV
static constexpr double MeV
Definition: Units.h:129
train.filename
string filename
Definition: train.py:213
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
gOptInpFileName
string gOptInpFileName
input file name
Definition: gNtpConv.cxx:214
genie::GHepParticle::Pz
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
generate_CCQE_events.target
target
Definition: generate_CCQE_events.py:44
genie::Interaction::AsString
string AsString(void) const
Definition: Interaction.cxx:249
genie::GHepParticle::Status
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
test.weight
weight
Definition: test.py:257
setprecision
Q_EXPORT QTSManip setprecision(int p)
Definition: qtextstream.h:343
genie::GHepParticle::Energy
double Energy(void) const
Get energy.
Definition: GHepParticle.h:92
freeze_graph.argv
argv
Definition: freeze_graph.py:218
genie::XclsTag
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
ConvertToGINuke
void ConvertToGINuke(void)
Definition: gNtpConv.cxx:2767
genie::GHepParticle::Px
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
genie::kPdgK0
const int kPdgK0
Definition: PDGCodes.h:174
genie::Kinematics::y
double y(bool selected=false) const
Definition: Kinematics.cxx:112
genie::GHepParticle::LastMother
int LastMother(void) const
Definition: GHepParticle.h:67
genie::XclsTag::IsCharmEvent
bool IsCharmEvent(void) const
Definition: XclsTag.h:50
genie::GHepParticle::Vt
double Vt(void) const
Get production time.
Definition: GHepParticle.h:97
CheckRootFilename
bool CheckRootFilename(string filename)
Definition: gEvComp.cxx:2033
genie::ProcessInfo::IsSingleKaon
bool IsSingleKaon(void) const
Definition: ProcessInfo.cxx:79
genie::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:63
genie::GHepParticle::FirstMother
int FirstMother(void) const
Definition: GHepParticle.h:66
GNtpcFmt_t
enum EGNtpcFmt GNtpcFmt_t
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:135
DefaultOutputFile
string DefaultOutputFile(void)
Definition: gNtpConv.cxx:3098
genie::utils::system::GenieMinorVrsNum
int GenieMinorVrsNum(string tag)
Definition: SystemUtils.cxx:66
genie::Interaction
Summary information for an interaction.
Definition: Interaction.h:56
gFileMinorVrs
int gFileMinorVrs
Definition: gNtpConv.cxx:224
genie::pdg::GeantToPdg
int GeantToPdg(int geant_code)
Definition: PDGUtils.cxx:416
genie::GHepParticle::LastDaughter
int LastDaughter(void) const
Definition: GHepParticle.h:69
genie::ProcessInfo::IsWeakNC
bool IsWeakNC(void) const
Definition: ProcessInfo.cxx:188
genie::InitialState::IsNuP
bool IsNuP(void) const
is neutrino + proton?
Definition: InitialState.cxx:199
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::ProcessInfo::IsCoherentElastic
bool IsCoherentElastic(void) const
Definition: ProcessInfo.cxx:104
ConvertToGRooTracker
void ConvertToGRooTracker(void)
Definition: gNtpConv.cxx:1771
genie::ProcessInfo::IsNuElectronElastic
bool IsNuElectronElastic(void) const
Definition: ProcessInfo.cxx:116
genie::units::cm2
static constexpr double cm2
Definition: Units.h:69
pass
HLTPathStatus const pass
Definition: EventSelExc_t.cpp:48
genie::kPdgKM
const int kPdgKM
Definition: PDGCodes.h:173
genie::ProcessInfo
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
genie::ProcessInfo::IsAMNuGamma
bool IsAMNuGamma(void) const
Definition: ProcessInfo.cxx:151
reco_momentum_tuples.t
t
Definition: reco_momentum_tuples.py:25
genie::kPdgGamma
const int kPdgGamma
Definition: PDGCodes.h:189
dump_to_simple_cpp.W
W
Definition: dump_to_simple_cpp.py:43
gOptN
Long64_t gOptN
number of events to process
Definition: gNtpConv.cxx:218
genie::Interaction::Kine
const Kinematics & Kine(void) const
Definition: Interaction.h:71
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
vector::size
size_t size
Definition: lodepng.cpp:55
genie::kPdgIndep
const int kPdgIndep
Definition: PDGCodes.h:231
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:172
genie::flux::fNgmax
const int fNgmax
Definition: GJPARCNuFlux.h:151
genie::NtpMCTreeHeader
MINOS-style Ntuple Class to hold an output MC Tree Header.
Definition: NtpMCTreeHeader.h:37
test.p
p
Definition: test.py:223
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
genie::kPdgString
const int kPdgString
Definition: PDGCodes.h:230
genie::NtpWriter::Save
void Save(void)
get the even tree
Definition: NtpWriter.cxx:225
pINFO
#define pINFO
Definition: Messenger.h:62
genie::kPdgAntiK0
const int kPdgAntiK0
Definition: PDGCodes.h:175
genie::kPdgK0L
const int kPdgK0L
Definition: PDGCodes.h:176
genie::NtpWriter::AddEventRecord
void AddEventRecord(int ievent, const EventRecord *ev_rec)
save the event tree
Definition: NtpWriter.cxx:57
genie::XclsTag::Resonance
Resonance_t Resonance(void) const
Definition: XclsTag.h:69
HAProbeFSI
int HAProbeFSI(int, int, int, double[], int[], int, int, int)
Definition: gNtpConv.cxx:3162
gOptOutFileFormat
GNtpcFmt_t gOptOutFileFormat
output file format id
Definition: gNtpConv.cxx:216
genie::GHepParticle::PolzIsSet
bool PolzIsSet(void) const
Definition: GHepParticle.cxx:308
pWARN
#define pWARN
Definition: Messenger.h:60
setw
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
genie::ProcessInfo::IsMEC
bool IsMEC(void) const
Definition: ProcessInfo.cxx:156
genie::ProcessInfo::IsEM
bool IsEM(void) const
Definition: ProcessInfo.cxx:173
cache_state.m
m
Definition: cache_state.py:415
genie::EventRecord
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:37
genie::ProcessInfo::IsDeepInelastic
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:84
reco_momentum_tuples.pz
pz
Definition: reco_momentum_tuples.py:220
genie::NtpWriter::Initialize
void Initialize(void)
add event
Definition: NtpWriter.cxx:83
genie::GHepParticle::KinE
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
Definition: GHepParticle.cxx:159
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
genie::RunOpt::Instance
static RunOpt * Instance(void)
Definition: RunOpt.cxx:54
genie::utils::system::GenieRevisVrsNum
int GenieRevisVrsNum(string tag)
Definition: SystemUtils.cxx:73
gOptVersion
int gOptVersion
output file format version
Definition: gNtpConv.cxx:217
genie::NtpWriter
A utility class to facilitate creating the GENIE MC Ntuple from the output GENIE GHEP event records...
Definition: NtpWriter.h:39
genie::pdg::IsIon
bool IsIon(int pdgc)
Definition: PDGUtils.cxx:39
genie::PDGLibrary
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:32
genie::Target::HitQrkIsSet
bool HitQrkIsSet(void) const
Definition: Target.cxx:292
genie::GHepParticle::Vz
double Vz(void) const
Get production z.
Definition: GHepParticle.h:96
genie::GHepParticle::X4
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
genie::kPdgAntiNeutron
const int kPdgAntiNeutron
Definition: PDGCodes.h:84
ConvertToGHepMock
void ConvertToGHepMock(void)
Definition: gNtpConv.cxx:1246
genie::pdg::IsPseudoParticle
bool IsPseudoParticle(int pdgc)
Definition: PDGUtils.cxx:24
genie::Interaction::ExclTag
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
E
E
Definition: 018_def.c:13
ConvertToGTracker
void ConvertToGTracker(void)
Definition: gNtpConv.cxx:1318
genie::InitialState::IsNuBarP
bool IsNuBarP(void) const
is anti-neutrino + proton?
Definition: InitialState.cxx:217
genie::RandomGen::RndGen
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:80
genie::kPdgAntiProton
const int kPdgAntiProton
Definition: PDGCodes.h:82
A
#define A
Definition: memgrp.cpp:38
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:159
gFileMajorVrs
int gFileMajorVrs
Definition: gNtpConv.cxx:223
run_cnn_1class.pdg
pdg
Definition: run_cnn_1class.py:54
reco_momentum_tuples.count
int count
Definition: reco_momentum_tuples.py:202
genie::GHepRecord::AddParticle
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:491
genie::Interaction::InitState
const InitialState & InitState(void) const
Definition: Interaction.h:69
LatestFormatVersionNumber
int LatestFormatVersionNumber(void)
Definition: gNtpConv.cxx:3135
genie::Interaction::ProcInfo
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
train.x
list x
Definition: train.py:276
genie::Kinematics::t
double t(bool selected=false) const
Definition: Kinematics.cxx:170
genie::pdg::IonPdgCodeToZ
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:52
gOptRanSeed
long int gOptRanSeed
random number seed
Definition: gNtpConv.cxx:220
genie::utils::app_init::MesgThresholds
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:99
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
EGNtpcFmt
EGNtpcFmt
Definition: gNtpConvDUNE.cxx:81
genie::utils::T2KEvGenMetaData
Utility class to store MC job meta-data.
Definition: T2KEvGenMetaData.h:38
Write
hadnt Write("hadnt")
genie::NtpMCRecHeader
MINOS-style Ntuple Class to hold an MC Event Record Header.
Definition: NtpMCRecHeader.h:32
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
genie::GHepParticle::Vy
double Vy(void) const
Get production y.
Definition: GHepParticle.h:95
genie::kPdgCluster
const int kPdgCluster
Definition: PDGCodes.h:229
genie::CmdLnArgParser
Command line argument parser.
Definition: CmdLnArgParser.h:29
genie::GHepParticle::A
int A(void) const
Definition: GHepParticle.cxx:184
genie::Kinematics::Q2
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
getreco.cmd
list cmd
Definition: getreco.py:22
genie::InitialState::Tgt
const Target & Tgt(void) const
Definition: InitialState.h:66
reco_momentum_tuples.fout
fout
Definition: reco_momentum_tuples.py:65
genie::GHepParticle::PolzAzimuthAngle
double PolzAzimuthAngle(void) const
Definition: GHepParticle.h:120
genie::NtpMCEventRecord::Clear
void Clear(Option_t *opt="")
Definition: NtpMCEventRecord.cxx:71
genie::kPdgK0S
const int kPdgK0S
Definition: PDGCodes.h:177
genie::units::ys
static constexpr double ys
Definition: Units.h:103
GNtpcFmt_t
enum EGNtpcFmt GNtpcFmt_t
genie::kPdgPositron
const int kPdgPositron
Definition: PDGCodes.h:36
genie::gAbortingInErr
bool gAbortingInErr
Definition: Messenger.cxx:34
genie::kPdgHadronicSyst
const int kPdgHadronicSyst
Definition: PDGCodes.h:210
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
setfill
Q_EXPORT QTSManip setfill(int f)
Definition: qtextstream.h:337
genie::GHepParticle
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
genie::CmdLnArgParser::OptionExists
bool OptionExists(char opt)
was option set?
Definition: CmdLnArgParser.cxx:88
gFileRevisVrs
int gFileRevisVrs
Definition: gNtpConv.cxx:225
genie::NtpMCEventRecord::event
EventRecord * event
event
Definition: NtpMCEventRecord.h:55
make_resolution.tree
tree
Definition: make_resolution.py:22
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
event
Event finding and building.
Definition: EventCheater_module.cc:32
genie::GHepStatus_t
enum genie::EGHepStatus GHepStatus_t
GetCommandLineArgs
void GetCommandLineArgs(int argc, char **argv)
Definition: gNtpConv.cxx:2980
genie::GHepParticle::Vx
double Vx(void) const
Get production x.
Definition: GHepParticle.h:94
PrintSyntax
void PrintSyntax(void)
Definition: gNtpConv.cxx:3152
genie::InitialState
Initial State information.
Definition: InitialState.h:48
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
genie::GHepParticle::Py
double Py(void) const
Get Py.
Definition: GHepParticle.h:89
gOptOutFileName
string gOptOutFileName
output file name
Definition: gNtpConv.cxx:215
Generated by   doxygen 1.8.11