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