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