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