neut_code_from_rootracker.C
Go to the documentation of this file.
1 //________________________________________________________________________________________
2 /*!
3 
4 \macro neut_code_from_rootracker.C
5 
6 \brief Macro to read a GENIE event tree in the t2k_rootracker format and calculate the
7  NEUT reaction code. Note that GENIE event files in t2k_rootracker format for
8  versions >= v2.5.1 already include the neut reaction code as a separate tree branch.
9 
10 \usage shell% root
11  root[0] .L neut_code_from_rootracker.C++
12  root[1] neut_code_from_rootracker("./your_rootracker_file.root");
13 
14 \author Costas Andreopoulos <costas.andreopoulos@stfc.ac.uk>
15  University of Liverpool & STFC Rutherford Appleton Lab
16 
17 \created Nov 24, 2008
18 
19 \cpright Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
20  For the full text of the license visit http://copyright.genie-mc.org
21  or see $GENIE/LICENSE
22 */
23 //_________________________________________________________________________________________
24 
25 #include <cassert>
26 #include <iostream>
27 #include <fstream>
28 #include <sstream>
29 #include <string>
30 
31 #include <TFile.h>
32 #include <TTree.h>
33 #include <TBranch.h>
34 #include <TBits.h>
35 #include <TObjString.h>
36 #include <TLorentzVector.h>
37 
38 using std::cout;
39 using std::endl;
40 using std::ostringstream;
41 using std::ofstream;
42 using std::string;
43 
44 void neut_code_from_rootracker(const char * filename)
45 {
46  bool using_new_version = false; // StdHepReScat and G2NeutEvtCode branches available only for versions >= 2.5.1
47  bool event_printout = false;
48 
49  //
50  // constants
51  //
52 
53  // set a max expected number of particles per event
54  const int kNPmax = 100;
55 
56  // status codes
57 //const int kIStUndefined = -1;
58 //const int kIStInitialState = 0;
59  const int kIStStableFinalState = 1;
60 //const int kIStIntermediateState = 2;
61  const int kIStDecayedState = 3;
62 //const int kIStCorrelatedNucleon = 10;
63 //const int kIStNucleonTarget = 11;
64 //const int kIStDISPreFragmHadronicState = 12;
65 //const int kIStPreDecayResonantState = 13;
66  const int kIStHadronInTheNucleus = 14;
67 
68  // define some particle codes needed for converting GENIE -> NEUT reaction codes
69  const int kPdgNuE = 12;
70  const int kPdgAntiNuE = -12;
71  const int kPdgNuMu = 14;
72  const int kPdgAntiNuMu = -14;
73  const int kPdgNuTau = 16;
74  const int kPdgAntiNuTau = -16;
75  const int kPdgGamma = 22; // photon
76  const int kPdgProton = 2212;
77 //const int kPdgAntiProton = -2212;
78  const int kPdgNeutron = 2112;
79 //const int kPdgAntiNeutron = -2112;
80  const int kPdgPiP = 211; // pi+
81  const int kPdgPiM = -211; // pi-
82  const int kPdgPi0 = 111; // pi0
83  const int kPdgEta = 221; // eta
84  const int kPdgKP = 321; // K+
85  const int kPdgKM = -321; // K-
86  const int kPdgK0 = 311; // K0
87  const int kPdgAntiK0 = -311; // \bar{K0}
88  const int kPdgLambda = 3122; // Lambda
89  const int kPdgAntiLambda = -3122; // \bar{Lambda}
90 
91  // stdhep momentum array indices
92  const int kStdHepIdxPx = 0;
93  const int kStdHepIdxPy = 1;
94  const int kStdHepIdxPz = 2;
95  const int kStdHepIdxE = 3;
96 
97  //
98  // get input tree
99  //
100 
101  TFile file(filename, "READ");
102  TTree * tree = (TTree *) file.Get("gRooTracker");
103  assert(tree);
104 
105  //
106  // event info in rootracker files
107  //
108 
109  TBits* EvtFlags = 0; // generator-specific event flags
110  TObjString* EvtCode = 0; // generator-specific string with 'event code'
111  int EvtNum; // event num.
112  double EvtXSec; // cross section for selected event (1E-38 cm2)
113  double EvtDXSec; // cross section for selected event kinematics (1E-38 cm2 /{K^n})
114  double EvtWght; // weight for that event
115  double EvtProb; // probability for that event (given cross section, path lengths, etc)
116  double EvtVtx[4]; // event vertex position in detector coord syst (in geom units)
117  int StdHepN; // number of particles in particle array
118  int StdHepPdg [kNPmax]; // stdhep-like particle array: pdg codes (& generator specific codes for pseudoparticles)
119  int StdHepStatus[kNPmax]; // stdhep-like particle array: generator-specific status code
120  int StdHepRescat[kNPmax]; // stdhep-like particle array: intranuclear rescattering code [ >= v2.5.1 ]
121  double StdHepX4 [kNPmax][4]; // stdhep-like particle array: 4-x (x, y, z, t) of particle in hit nucleus frame (fm)
122  double StdHepP4 [kNPmax][4]; // stdhep-like particle array: 4-p (px,py,pz,E) of particle in LAB frame (GeV)
123  double StdHepPolz [kNPmax][3]; // stdhep-like particle array: polarization vector
124  int StdHepFd [kNPmax]; // stdhep-like particle array: first daughter
125  int StdHepLd [kNPmax]; // stdhep-like particle array: last daughter
126  int StdHepFm [kNPmax]; // stdhep-like particle array: first mother
127  int StdHepLm [kNPmax]; // stdhep-like particle array: last mother
128  int G2NeutEvtCode; // NEUT code for the current GENIE event [ >= v2.5.1 ]
129  int NuParentPdg; // parent hadron pdg code
130  int NuParentDecMode; // parent hadron decay mode
131  double NuParentDecP4 [4]; // parent hadron 4-momentum at decay
132  double NuParentDecX4 [4]; // parent hadron 4-position at decay
133  double NuParentProP4 [4]; // parent hadron 4-momentum at production
134  double NuParentProX4 [4]; // parent hadron 4-position at production
135  int NuParentProNVtx; // parent hadron vtx id
136 
137  //
138  // get tre branches and set branch addresses
139  //
140 
141  TBranch * brEvtFlags = tree -> GetBranch ("EvtFlags");
142  TBranch * brEvtCode = tree -> GetBranch ("EvtCode");
143  TBranch * brEvtNum = tree -> GetBranch ("EvtNum");
144  TBranch * brEvtXSec = tree -> GetBranch ("EvtXSec");
145  TBranch * brEvtDXSec = tree -> GetBranch ("EvtDXSec");
146  TBranch * brEvtWght = tree -> GetBranch ("EvtWght");
147  TBranch * brEvtProb = tree -> GetBranch ("EvtProb");
148  TBranch * brEvtVtx = tree -> GetBranch ("EvtVtx");
149  TBranch * brStdHepN = tree -> GetBranch ("StdHepN");
150  TBranch * brStdHepPdg = tree -> GetBranch ("StdHepPdg");
151  TBranch * brStdHepStatus = tree -> GetBranch ("StdHepStatus");
152  TBranch * brStdHepRescat = (using_new_version) ? tree -> GetBranch ("StdHepRescat") : 0;
153  TBranch * brStdHepX4 = tree -> GetBranch ("StdHepX4");
154  TBranch * brStdHepP4 = tree -> GetBranch ("StdHepP4");
155  TBranch * brStdHepPolz = tree -> GetBranch ("StdHepPolz");
156  TBranch * brStdHepFd = tree -> GetBranch ("StdHepFd");
157  TBranch * brStdHepLd = tree -> GetBranch ("StdHepLd");
158  TBranch * brStdHepFm = tree -> GetBranch ("StdHepFm");
159  TBranch * brStdHepLm = tree -> GetBranch ("StdHepLm");
160  TBranch * brG2NeutEvtCode = (using_new_version) ? tree -> GetBranch ("G2NeutEvtCode") : 0;
161  TBranch * brNuParentPdg = tree -> GetBranch ("NuParentPdg");
162  TBranch * brNuParentDecMode = tree -> GetBranch ("NuParentDecMode");
163  TBranch * brNuParentDecP4 = tree -> GetBranch ("NuParentDecP4");
164  TBranch * brNuParentDecX4 = tree -> GetBranch ("NuParentDecX4");
165  TBranch * brNuParentProP4 = tree -> GetBranch ("NuParentProP4");
166  TBranch * brNuParentProX4 = tree -> GetBranch ("NuParentProX4");
167  TBranch * brNuParentProNVtx = tree -> GetBranch ("NuParentProNVtx");
168 
169  brEvtFlags -> SetAddress ( &EvtFlags );
170  brEvtCode -> SetAddress ( &EvtCode );
171  brEvtNum -> SetAddress ( &EvtNum );
172  brEvtXSec -> SetAddress ( &EvtXSec );
173  brEvtDXSec -> SetAddress ( &EvtDXSec );
174  brEvtWght -> SetAddress ( &EvtWght );
175  brEvtProb -> SetAddress ( &EvtProb );
176  brEvtVtx -> SetAddress ( EvtVtx );
177  brStdHepN -> SetAddress ( &StdHepN );
178  brStdHepPdg -> SetAddress ( StdHepPdg );
179  brStdHepStatus -> SetAddress ( StdHepStatus );
180  if(using_new_version) {
181  brStdHepRescat -> SetAddress ( StdHepRescat );
182  }
183  brStdHepX4 -> SetAddress ( StdHepX4 );
184  brStdHepP4 -> SetAddress ( StdHepP4 );
185  brStdHepPolz -> SetAddress ( StdHepPolz );
186  brStdHepFd -> SetAddress ( StdHepFd );
187  brStdHepLd -> SetAddress ( StdHepLd );
188  brStdHepFm -> SetAddress ( StdHepFm );
189  brStdHepLm -> SetAddress ( StdHepLm );
190  if(using_new_version) {
191  brG2NeutEvtCode -> SetAddress ( &G2NeutEvtCode );
192  }
193  brNuParentPdg -> SetAddress ( &NuParentPdg );
194  brNuParentDecMode -> SetAddress ( &NuParentDecMode );
195  brNuParentDecP4 -> SetAddress ( NuParentDecP4 );
196  brNuParentDecX4 -> SetAddress ( NuParentDecX4 );
197  brNuParentProP4 -> SetAddress ( NuParentProP4 );
198  brNuParentProX4 -> SetAddress ( NuParentProX4 );
199  brNuParentProNVtx -> SetAddress ( &NuParentProNVtx );
200 
201 
202  //
203  // open a text file to save the reaction codes
204  //
205  ostringstream outfilename;
206  outfilename << filename << ".reaction_codes";
207 
208  ofstream outfile(outfilename.str().c_str());
209  outfile << "#" << endl;
210  outfile << "# NEUT reaction code for GENIE file: " << filename << endl;
211  outfile << "#" << endl;
212  outfile << "#" << endl;
213  outfile << "# GENIE evt nu. NEUT code" << endl;
214 
215  //
216  // event loop
217  //
218 
219  int n = tree->GetEntries();
220  printf("Number of entries: %d", n);
221 
222  for(int i=0; i < tree->GetEntries(); i++) {
223 
224  //
225  // get next event
226  //
227 
228  printf("\n\n ** Current entry: %d \n", i);
229  tree->GetEntry(i);
230 
231  //
232  // print event
233  //
234  if(event_printout) {
235  printf("\n -----------------------------------------------------------------------------------------------------------------");
236  printf("\n Event code : %s", EvtCode->String().Data());
237  printf("\n Event x-section : %10.5f * 1E-38* cm^2", EvtXSec);
238  printf("\n Event kinematics x-section : %10.5f * 1E-38 * cm^2/{K^n}", EvtDXSec);
239  printf("\n Event weight : %10.8f", EvtWght);
240  printf("\n Event vertex : x = %8.2f mm, y = %8.2f mm, z = %8.2f mm", EvtVtx[0], EvtVtx[1], EvtVtx[2]);
241  printf("\n * Particle list:");
242  printf("\n --------------------------------------------------------------------------------------------------------------------------");
243  printf("\n | Idx | Ist | PDG | Rescat | Mother | Daughter | Px | Py | Pz | E | x | y | z |");
244  printf("\n | | | | | | |(GeV/c) |(GeV/c) |(GeV/c) | (GeV) | (fm) | (fm) | (fm) |");
245  printf("\n --------------------------------------------------------------------------------------------------------------------------");
246 
247  for(int ip=0; ip<StdHepN; ip++) {
248  printf("\n | %3d | %3d | %10d | %6d | %3d | %3d | %3d | %3d | %6.3f | %6.3f | %6.3f | %6.3f | %6.3f | %6.3f | %6.3f |",
249  ip, StdHepStatus[ip], StdHepPdg[ip], StdHepRescat[ip],
250  StdHepFm[ip], StdHepLm[ip], StdHepFd[ip], StdHepLd[ip],
251  StdHepP4[ip][0], StdHepP4[ip][1], StdHepP4[ip][2], StdHepP4[ip][3],
252  StdHepX4[ip][0], StdHepX4[ip][1], StdHepX4[ip][2]);
253  }
254  printf("\n --------------------------------------------------------------------------------------------------------------------------");
255  printf("\n * Flux Info:");
256  printf("\n Parent hadron pdg code : %d", NuParentPdg);
257  printf("\n Parent hadron decay mode : %d", NuParentDecMode);
258  printf("\n Parent hadron 4p at decay : Px = %6.3f GeV/c, Py = %6.3f GeV/c, Pz = %6.3f GeV/c, E = %6.3f GeV",
259  NuParentDecP4[0], NuParentDecP4[1], NuParentDecP4[2], NuParentDecP4[3]);
260  printf("\n Parent hadron 4p at prod. : Px = %6.3f GeV/c, Py = %6.3f GeV/c, Pz = %6.3f GeV/c, E = %6.3f GeV",
261  NuParentProP4[0], NuParentProP4[1], NuParentProP4[2], NuParentProP4[3]);
262  printf("\n Parent hadron 4x at decay : x = %6.3f m, y = %6.3f m, z = %6.3f m, t = %6.3f s",
263  NuParentDecX4[0], NuParentDecX4[1], NuParentDecX4[2], NuParentDecX4[3]);
264  printf("\n Parent hadron 4x at prod. : x = %6.3f m, y = %6.3f m, z = %6.3f m, t = %6.3f s",
265  NuParentProX4[0], NuParentProX4[1], NuParentProX4[2], NuParentProX4[3]);
266  printf("\n -------------------------------------------------------------------------------------------------------------------------- \n");
267  }
268 
269  //
270  // figure out the NEUT code for the current GENIE event
271  //
272 
273  int evtype = 0;
274 
275  // basic GENIE reaction types
276  string genie_code = EvtCode->GetString().Data();
277  bool is_cc = (genie_code.find("Weak[CC]") != string::npos);
278  bool is_nc = (genie_code.find("Weak[NC]") != string::npos);
279  bool is_charm = (genie_code.find("charm") != string::npos);
280  bool is_qel = (genie_code.find("QES") != string::npos);
281  bool is_dis = (genie_code.find("DIS") != string::npos);
282  bool is_res = (genie_code.find("RES") != string::npos);
283  bool is_cohpi = (genie_code.find("COH") != string::npos);
284  bool is_ve = (genie_code.find("NuEEL") != string::npos);
285  bool is_imd = (genie_code.find("IMD") != string::npos);
286 
287  // basic initial state info
288  int inu_pos = 0;
289  int fsl_pos = StdHepFd[inu_pos]; // final state primary lepton: neutrino daughter
290  assert(fsl_pos>0);
291  int tgt_pos = 1;
292  int hitnuc_pos = -1; // can be at slot 2 (for nuclear targets), slot 1 (for free nuc targets) or undefined (for ve-, IMD, coherent etc)
293 
294  int nu_code = StdHepPdg[inu_pos];
295  bool is_nu = (nu_code == kPdgNuE || nu_code == kPdgNuMu || nu_code == kPdgNuTau );
296  bool is_nubar = (nu_code == kPdgAntiNuE || nu_code == kPdgAntiNuMu || nu_code == kPdgAntiNuTau);
297 
298  int target_code = StdHepPdg[tgt_pos];
299  bool nuclear_target = (target_code > 1000000000);
300 
301  bool hitnuc_set = false;
302  bool is_p = false;
303  bool is_n = false;
304  for(int ip = 1; ip <= 2; ip++)
305  {
306  int ghep_pdgc = StdHepPdg[ip];
307  if(ghep_pdgc == kPdgProton) {
308  hitnuc_pos = ip;
309  hitnuc_set = true;
310  is_p = true;
311  break;
312  }
313  if(ghep_pdgc == kPdgNeutron) {
314  hitnuc_pos = ip;
315  hitnuc_set = true;
316  is_n = true;
317  break;
318  }
319  }
320 
321  // calculate the true (at the hit nucleon rest frame) hadronic invariant mass
322  // calculate only if the hit nucleon is set (calc not needed for coherent events etc)
323 
324  bool W_gt_2 = false;
325  if(hitnuc_pos != -1)
326  {
327  TLorentzVector p4v ( StdHepP4 [inu_pos] [kStdHepIdxPx],
328  StdHepP4 [inu_pos] [kStdHepIdxPy],
329  StdHepP4 [inu_pos] [kStdHepIdxPz],
330  StdHepP4 [inu_pos] [kStdHepIdxE] );
331  TLorentzVector p4l ( StdHepP4 [fsl_pos] [kStdHepIdxPx],
332  StdHepP4 [fsl_pos] [kStdHepIdxPy],
333  StdHepP4 [fsl_pos] [kStdHepIdxPz],
334  StdHepP4 [fsl_pos] [kStdHepIdxE] );
335  TLorentzVector p4n ( StdHepP4 [hitnuc_pos] [kStdHepIdxPx],
336  StdHepP4 [hitnuc_pos] [kStdHepIdxPy],
337  StdHepP4 [hitnuc_pos] [kStdHepIdxPz],
338  StdHepP4 [hitnuc_pos] [kStdHepIdxE] );
339 
340  TLorentzVector q = p4v - p4l;
341  TLorentzVector w = p4n + q;
342 
343  double W = w.Mag();
344  W_gt_2 = (W > 2.0);
345  }
346 
347  // (quasi-)elastic, nc+cc, nu+nubar
348  //
349  if (is_qel && !is_charm && is_cc && is_nu ) evtype = 1;
350  else if (is_qel && !is_charm && is_nc && is_nu && is_p ) evtype = 51;
351  else if (is_qel && !is_charm && is_nc && is_nu && is_n ) evtype = 52;
352  else if (is_qel && !is_charm && is_cc && is_nubar ) evtype = -1;
353  else if (is_qel && !is_charm && is_nc && is_nubar && is_p) evtype = -51;
354  else if (is_qel && !is_charm && is_nc && is_nubar && is_n) evtype = -52;
355 
356  // quasi-elastic charm production
357  //
358  else if (is_qel && is_charm && is_cc && is_nu ) evtype = 25;
359  else if (is_qel && is_charm && is_cc && is_nubar ) evtype = -25;
360 
361  // inverse mu- (tau-) decay and ve- elastic
362  //
363  else if ( is_imd ) evtype = 9;
364  else if ( is_ve ) evtype = 59;
365 
366  // coherent pi, nc+cc, nu+nubar
367  //
368  else if (is_cohpi && is_cc && is_nu ) evtype = 16;
369  else if (is_cohpi && is_cc && is_nubar) evtype = -16;
370  else if (is_cohpi && is_nc && is_nu ) evtype = 36;
371  else if (is_cohpi && is_nc && is_nubar) evtype = -36;
372 
373  // dis, W>2, nc+cc, nu+nubar
374  // (charm DIS not simulated by NEUT, will bundle GENIE charm DIS into this category)
375  //
376  else if (is_dis && W_gt_2 && is_cc && is_nu ) evtype = 26;
377  else if (is_dis && W_gt_2 && is_nc && is_nu ) evtype = 46;
378  else if (is_dis && W_gt_2 && is_cc && is_nubar) evtype = -26;
379  else if (is_dis && W_gt_2 && is_nc && is_nubar) evtype = -46;
380 
381  // resonance or dis with W < 2 GeV
382  //
383  else if ( is_res || (is_dis && !W_gt_2) ) {
384 
385  //cout << "Current event is RES or DIS with W<2" << endl;
386 
387  // check the number of pions and nucleons in the primary hadronic system
388  // (_before_ intranuclear rescattering)
389  //
390  int nn=0, np=0, npi0=0, npip=0, npim=0, nKp=0, nKm=0, nK0=0, neta=0, nlambda=0, ngamma=0;
391 
392  for(int ip = 0; ip < StdHepN; ip++)
393  {
394  int ghep_ist = StdHepStatus[ip];
395  int ghep_pdgc = StdHepPdg[ip];
396  int ghep_fm = StdHepFm[ip];
397  int ghep_fmpdgc = (ghep_fm==-1) ? 0 : StdHepPdg[ghep_fm];
398 
399  // For nuclear targets use hadrons marked as 'hadron in the nucleus'
400  // which are the ones passed in the intranuclear rescattering
401  // For free nucleon targets use particles marked as 'final state'
402  // but make an exception for decayed pi0's,eta's (count them and not their daughters)
403 
404  bool decayed = (ghep_ist==kIStDecayedState && (ghep_pdgc==kPdgPi0 || ghep_pdgc==kPdgEta));
405  bool parent_included = (ghep_fmpdgc==kPdgPi0 || ghep_fmpdgc==kPdgEta);
406 
407  bool count_it =
408  ( nuclear_target && ghep_ist==kIStHadronInTheNucleus) ||
409  (!nuclear_target && decayed) ||
410  (!nuclear_target && ghep_ist==kIStStableFinalState && !parent_included);
411 
412  if(!count_it) continue;
413 
414  if(ghep_pdgc == kPdgProton ) np++; // p
415  if(ghep_pdgc == kPdgNeutron) nn++; // n
416  if(ghep_pdgc == kPdgPiP) npip++; // pi+
417  if(ghep_pdgc == kPdgPiM) npim++; // pi-
418  if(ghep_pdgc == kPdgPi0) npi0++; // pi0
419  if(ghep_pdgc == kPdgEta) neta++; // eta0
420  if(ghep_pdgc == kPdgKP) nKp++; // K+
421  if(ghep_pdgc == kPdgKM) nKm++; // K-
422  if(ghep_pdgc == kPdgK0) nK0++; // K0
423  if(ghep_pdgc == kPdgAntiK0) nK0++; // K0
424  if(ghep_pdgc == kPdgLambda) nlambda++; // Lamda
425  if(ghep_pdgc == kPdgAntiLambda) nlambda++; // Lamda
426  if(ghep_pdgc == kPdgGamma) ngamma++; // photon
427  }
428  if(event_printout) {
429  cout << "Num of primary particles: \n p = " << np << ", n = " << nn
430  << ", pi+ = " << npip << ", pi- = " << npim << ", pi0 = " << npi0
431  << ", eta = " << neta
432  << ", K+ = " << nKp << ", K- = " << nKm << ", K0 = " << nK0
433  << ", Lambda's = " << nlambda
434  << ", gamma's = " << ngamma
435  << endl;
436  }
437  int nnuc = np + nn;
438  int npi = npi0 + npip + npim;
439  int nK = nK0 + nKp + nKm;
440  int neKL = neta + nK + nlambda;
441 
442  bool is_radiative_dec = (nnuc==1) && (npi==0) && (ngamma==1);
443 
444  //
445  // single gamma from resonances
446  //
447 
448  if (is_res && is_nu && is_cc && is_n && is_radiative_dec) evtype = 17;
449  else if (is_res && is_nu && is_nc && is_n && is_radiative_dec) evtype = 38;
450  else if (is_res && is_nu && is_nc && is_p && is_radiative_dec) evtype = 39;
451 
452  else if (is_res && is_nubar && is_cc && is_p && is_radiative_dec) evtype = -17;
453  else if (is_res && is_nubar && is_nc && is_n && is_radiative_dec) evtype = -38;
454  else if (is_res && is_nubar && is_nc && is_p && is_radiative_dec) evtype = -39;
455 
456  //
457  // single pi (res + non-res bkg)
458  //
459 
460  // nu CC
461  else if (is_nu && is_cc && is_p && np==1 && nn==0 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = 11;
462  else if (is_nu && is_cc && is_n && np==1 && nn==0 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = 12;
463  else if (is_nu && is_cc && is_n && np==0 && nn==1 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = 13;
464 
465  // nu NC
466  else if (is_nu && is_nc && is_n && np==0 && nn==1 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = 31;
467  else if (is_nu && is_nc && is_p && np==1 && nn==0 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = 32;
468  else if (is_nu && is_nc && is_n && np==1 && nn==0 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = 33;
469  else if (is_nu && is_nc && is_p && np==0 && nn==1 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = 34;
470 
471  //nubar CC
472  else if (is_nubar && is_cc && is_n && np==0 && nn==1 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = -11;
473  else if (is_nubar && is_cc && is_p && np==0 && nn==1 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = -12;
474  else if (is_nubar && is_cc && is_p && np==1 && nn==0 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = -13;
475 
476  //nubar NC
477  else if (is_nubar && is_nc && is_n && np==0 && nn==1 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = -31;
478  else if (is_nubar && is_nc && is_p && np==1 && nn==0 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = -32;
479  else if (is_nubar && is_nc && is_n && np==1 && nn==0 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = -33;
480  else if (is_nubar && is_nc && is_p && np==0 && nn==1 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = -34;
481 
482  //
483  // single eta from res
484  //
485 
486  else if (is_res && is_nu && is_cc && is_n && np==1 && nn==0 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = 22;
487  else if (is_res && is_nu && is_nc && is_n && np==0 && nn==1 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = 42;
488  else if (is_res && is_nu && is_nc && is_p && np==1 && nn==0 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = 43;
489 
490  else if (is_res && is_nubar && is_cc && is_p && np==0 && nn==1 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = -22;
491  else if (is_res && is_nubar && is_nc && is_n && np==0 && nn==1 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = -42;
492  else if (is_res && is_nubar && is_nc && is_p && np==1 && nn==0 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = -43;
493 
494  //
495  // single K from res
496  //
497 
498  else if (is_res && is_nu && is_cc && is_n && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = 23;
499  else if (is_res && is_nu && is_nc && is_n && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = 44;
500  else if (is_res && is_nu && is_nc && is_p && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = 45;
501 
502  else if (is_res && is_nubar && is_cc && is_p && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = -23;
503  else if (is_res && is_nubar && is_nc && is_n && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = -44;
504  else if (is_res && is_nubar && is_nc && is_p && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = -45;
505 
506  //
507  // multi-pi (res or dis), W<2GeV
508  //
509 
510  else if (is_nu && is_cc && npi>1) evtype = 21;
511  else if (is_nu && is_nc && npi>1) evtype = 41;
512  else if (is_nubar && is_cc && npi>1) evtype = -21;
513  else if (is_nubar && is_nc && npi>1) evtype = -41;
514 
515  //
516  // rare final state for RES or low-W (<2GeV) DIS events
517  // (eg K0\bar{K0} final states, N0(1720) -> Sigma- K+ res decays, etc)
518  // bundled-in with multi-pi
519  //
520  else {
521 
522  cout << "Rare RES/low-W DIS final state: Bundled-in with multi-pi events" << endl;
523 
524  if (is_nu && is_cc) evtype = 21;
525  else if (is_nu && is_nc) evtype = 41;
526  else if (is_nubar && is_cc) evtype = -21;
527  else if (is_nubar && is_nc) evtype = -41;
528  }
529  }
530 
531  cout << " *** GENIE event = " << i << " --> NEUT reaction code = " << evtype << endl;
532  if(using_new_version) {
533  // for validation, use a file generated + converted with the CVS head version of GENIE
534  // where the NeutCode is stored at the t2k_rootracker tree
535  cout << "NEUT reaction code stored at the rootracker file = " << G2NeutEvtCode << endl;
536  assert(evtype == G2NeutEvtCode);
537  }
538 
539  // save at the output file
540  outfile << "\t" << i << "\t" << evtype << endl;
541 
542  } // event loop
543 
544  outfile.close();
545 
546  file.Close();
547 }
genie::kPdgNuE
const int kPdgNuE
Definition: PDGCodes.h:28
GetBranch
ntuple GetBranch("organID") -> SetAddress(&xx)
genie::kPdgLambda
const int kPdgLambda
Definition: PDGCodes.h:69
genie::kPdgAntiNuE
const int kPdgAntiNuE
Definition: PDGCodes.h:29
genie::kPdgNuMu
const int kPdgNuMu
Definition: PDGCodes.h:30
kNPmax
const int kNPmax
Definition: gNtpConv.cxx:228
ip
TString ip
Definition: loadincs.C:5
shutoffs.filename
string filename
Definition: shutoffs.py:106
genie::kPdgK0
const int kPdgK0
Definition: PDGCodes.h:151
make_template_knob_config.outfilename
string outfilename
knobs that need extra care
Definition: make_template_knob_config.py:120
genie::kPdgKM
const int kPdgKM
Definition: PDGCodes.h:150
genie::kPdgGamma
const int kPdgGamma
Definition: PDGCodes.h:166
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:149
printf
printf("%d Experimental points found\n", nlines)
genie::kPdgEta
const int kPdgEta
Definition: PDGCodes.h:138
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:135
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:137
genie::kPdgAntiK0
const int kPdgAntiK0
Definition: PDGCodes.h:152
genie::kPdgAntiNuTau
const int kPdgAntiNuTau
Definition: PDGCodes.h:33
genie::kPdgAntiNuMu
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
om::cout
OStream cout
Definition: OStream.cxx:6
neut_code_from_rootracker
void neut_code_from_rootracker(const char *filename)
Definition: neut_code_from_rootracker.C:44
genie::kPdgNuTau
const int kPdgNuTau
Definition: PDGCodes.h:32
file
TFile * file
Definition: cellShifts.C:17
ana::assert
assert(nhit_max >=nhit_nbins)
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:136
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:65
genie::kPdgAntiLambda
const int kPdgAntiLambda
Definition: PDGCodes.h:70
w
Float_t w
Definition: plot.C:20
W
#define W(x)
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:67
outfile
FILE * outfile
Definition: dump_event.C:13
string
enum BeamMode string
Generated on Sat Feb 27 2021 14:16:12 for NOvA by   doxygen 1.8.11