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