Functions
genie::utils::ghep Namespace Reference

GHEP event record utilities. More...

Functions

int NeutReactionCode (const GHepRecord *evrec)
 
int NuanceReactionCode (const GHepRecord *evrec)
 

Detailed Description

GHEP event record utilities.

Author
Costas Andreopoulos <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

Nov 30, 2008

Copyright (c) 2003-2020, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Function Documentation

int genie::utils::ghep::NeutReactionCode ( const GHepRecord *  evrec)

Definition at line 22 of file GHepUtils.cxx.

23 {
24 // Ryan Terri, Yoshinari Hayato, Costas Andreopoulos
25 //
26 // A description of NEUT event types can be seen here:
27 // http://t2k.phy.duke.edu/bin/view/Main/NeutModes
28 // Any extension used here has been agreed with SK (Hayato et al)
29 //
30 // Updated for NEUT 5.4.0 by Christophe Bronner
31 
32 
33  if(!event) {
34  LOG("GHepUtils", pWARN) << "Null event!";
35  return 0;
36  }
37 
38  int evtype = 0;
39 
40  Interaction * interaction = event->Summary();
41 
42  const ProcessInfo & proc = interaction->ProcInfo();
43  const InitialState & init = interaction->InitState();
44  const XclsTag & xcls = interaction->ExclTag();
45  const Kinematics & kine = interaction->Kine();
46  const Target & tgt = init.Tgt();
47 
48  bool is_cc = proc.IsWeakCC();
49  bool is_nc = proc.IsWeakNC();
50  bool is_charm = xcls.IsCharmEvent();
51  bool is_qel = proc.IsQuasiElastic();
52  bool is_dis = proc.IsDeepInelastic();
53  bool is_res = proc.IsResonant();
54  bool is_coh_pr = proc.IsCoherentProduction();
55  bool is_ve = proc.IsNuElectronElastic();
56  bool is_mec = proc.IsMEC();
57  bool is_imd = proc.IsInverseMuDecay();
58  bool is_ask = proc.IsSingleKaon();
59  bool is_diff = proc.IsDiffractive();
60  bool is_p = tgt.HitNucIsSet() ? tgt.HitNucPdg()==kPdgProton : false;
61  bool is_n = tgt.HitNucIsSet() ? tgt.HitNucPdg()==kPdgNeutron : false;
62  bool is_nu = pdg::IsNeutrino (init.ProbePdg());
63  bool is_nubar = pdg::IsAntiNeutrino(init.ProbePdg());
64  bool W_gt_2 = kine.KVSet(kKVW) ? (kine.W() > 2.0) : false;
65 
66  // (quasi-)elastic, nc+cc, nu+nubar
67  //
68  if (is_qel && !is_charm && is_cc && is_nu ) evtype = 1;
69  else if (is_qel && !is_charm && is_nc && is_nu && is_p ) evtype = 51;
70  else if (is_qel && !is_charm && is_nc && is_nu && is_n ) evtype = 52;
71  else if (is_qel && !is_charm && is_cc && is_nubar ) evtype = -1;
72  else if (is_qel && !is_charm && is_nc && is_nubar && is_p) evtype = -51;
73  else if (is_qel && !is_charm && is_nc && is_nubar && is_n) evtype = -52;
74 
75  // MEC - only CC implemented in NEUT
76  else if (is_mec && !is_charm && is_cc && is_nu ) evtype = 2;
77  else if (is_mec && !is_charm && is_cc && is_nubar ) evtype = -2;
78 
79  // quasi-elastic charm production
80  // Part of the DIS W>2GeV mode in NEUT - CB
81  else if (is_qel && is_charm && is_cc && is_nu ) evtype = 26;
82  else if (is_qel && is_charm && is_cc && is_nubar ) evtype = -26;
83 
84  // inverse mu- (tau-) decay and ve- elastic
85  //Those modes don't actually exist in NEUT, 9 and 59 used as place holders
86  else if ( is_imd ) evtype = 9;
87  else if ( is_ve ) evtype = 59;
88 
89  // coherent pi, nc+cc, nu+nubar
90  //
91  else if (is_coh_pr && is_cc && is_nu ) evtype = 16;
92  else if (is_coh_pr && is_cc && is_nubar) evtype = -16;
93  else if (is_coh_pr && is_nc && is_nu ) evtype = 36;
94  else if (is_coh_pr && is_nc && is_nubar) evtype = -36;
95 
96  // dis, W>2, nc+cc, nu+nubar
97  // (charm DIS not simulated by NEUT, will bundle GENIE charm DIS into this category)
98  //
99  else if (is_dis && W_gt_2 && is_cc && is_nu ) evtype = 26;
100  else if (is_dis && W_gt_2 && is_nc && is_nu ) evtype = 46;
101  else if (is_dis && W_gt_2 && is_cc && is_nubar) evtype = -26;
102  else if (is_dis && W_gt_2 && is_nc && is_nubar) evtype = -46;
103 
104  // resonance or dis with W < 2 GeV or single kaon
105  //
106  else if ( is_res || (is_dis && !W_gt_2) || is_ask ) {
107 
108  LOG("GHepUtils", pNOTICE) << "Current event is RES or DIS with W<2";
109 
110  // check the number of pions and nucleons in the primary hadronic system
111  // (_before_ intranuclear rescattering)
112  //
113  int nn=0, np=0, npi0=0, npip=0, npim=0, nKp=0, nKm=0, nK0=0, neta=0, nlambda=0, ngamma=0;
114  bool nuclear_target = init.Tgt().IsNucleus();
115 
116  TIter event_iter(event);
117  GHepParticle * p = 0;
118 
119  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) )
120  {
121  GHepStatus_t ghep_ist = (GHepStatus_t) p->Status();
122  int ghep_pdgc = p->Pdg();
123  int ghep_fm = p->FirstMother();
124  int ghep_fmpdgc = (ghep_fm==-1) ? 0 : event->Particle(ghep_fm)->Pdg();
125 
126  // For nuclear targets use hadrons marked as 'hadron in the nucleus'
127  // which are the ones passed in the intranuclear rescattering
128  // For free nucleon targets use particles marked as 'final state'
129  // but make an exception for decayed pi0's,eta's (count them and not their daughters)
130 
131  bool decayed = (ghep_ist==kIStDecayedState && (ghep_pdgc==kPdgPi0 || ghep_pdgc==kPdgEta));
132  bool parent_included = (ghep_fmpdgc==kPdgPi0 || ghep_fmpdgc==kPdgEta);
133 
134  bool count_it =
135  ( nuclear_target && ghep_ist==kIStHadronInTheNucleus) ||
136  (!nuclear_target && decayed) ||
137  (!nuclear_target && ghep_ist==kIStStableFinalState && !parent_included);
138 
139  if(!count_it) continue;
140 
141  if(ghep_pdgc == kPdgProton ) np++; // p
142  if(ghep_pdgc == kPdgNeutron) nn++; // n
143  if(ghep_pdgc == kPdgPiP) npip++; // pi+
144  if(ghep_pdgc == kPdgPiM) npim++; // pi-
145  if(ghep_pdgc == kPdgPi0) npi0++; // pi0
146  if(ghep_pdgc == kPdgEta) neta++; // eta0
147  if(ghep_pdgc == kPdgKP) nKp++; // K+
148  if(ghep_pdgc == kPdgKM) nKm++; // K-
149  if(ghep_pdgc == kPdgK0) nK0++; // K0
150  if(ghep_pdgc == kPdgAntiK0) nK0++; // K0
151  if(ghep_pdgc == kPdgLambda) nlambda++; // Lamda
152  if(ghep_pdgc == kPdgAntiLambda) nlambda++; // Lamda
153  if(ghep_pdgc == kPdgGamma) ngamma++; // photon
154  }
155  LOG("GHepUtils", pNOTICE)
156  << "Num of primary particles: \n p = " << np << ", n = " << nn
157  << ", pi+ = " << npip << ", pi- = " << npim << ", pi0 = " << npi0
158  << ", eta = " << neta
159  << ", K+ = " << nKp << ", K- = " << nKm << ", K0 = " << nK0
160  << ", Lambda's = " << nlambda
161  << ", gamma's = " << ngamma;
162 
163  int nnuc = np + nn;
164  int npi = npi0 + npip + npim;
165  int nK = nK0 + nKp + nKm;
166  int neKL = neta + nK + nlambda;
167 
168  bool is_radiative_dec = (nnuc==1) && (npi==0) && (ngamma==1);
169 
170  //
171  // single gamma from resonances
172  //
173 
174  if (is_res && is_nu && is_cc && is_n && is_radiative_dec) evtype = 17;
175  else if (is_res && is_nu && is_nc && is_n && is_radiative_dec) evtype = 38;
176  else if (is_res && is_nu && is_nc && is_p && is_radiative_dec) evtype = 39;
177 
178  else if (is_res && is_nubar && is_cc && is_p && is_radiative_dec) evtype = -17;
179  else if (is_res && is_nubar && is_nc && is_n && is_radiative_dec) evtype = -38;
180  else if (is_res && is_nubar && is_nc && is_p && is_radiative_dec) evtype = -39;
181 
182  //
183  // single pi (res + non-res bkg)
184  //
185 
186  // nu CC
187  else if (is_nu && is_cc && is_p && np==1 && nn==0 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = 11;
188  else if (is_nu && is_cc && is_n && np==1 && nn==0 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = 12;
189  else if (is_nu && is_cc && is_n && np==0 && nn==1 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = 13;
190 
191  // nu NC
192  else if (is_nu && is_nc && is_n && np==0 && nn==1 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = 31;
193  else if (is_nu && is_nc && is_p && np==1 && nn==0 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = 32;
194  else if (is_nu && is_nc && is_n && np==1 && nn==0 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = 33;
195  else if (is_nu && is_nc && is_p && np==0 && nn==1 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = 34;
196 
197  //nubar CC
198  else if (is_nubar && is_cc && is_n && np==0 && nn==1 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = -11;
199  else if (is_nubar && is_cc && is_p && np==0 && nn==1 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = -12;
200  else if (is_nubar && is_cc && is_p && np==1 && nn==0 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = -13;
201 
202  //nubar NC
203  else if (is_nubar && is_nc && is_n && np==0 && nn==1 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = -31;
204  else if (is_nubar && is_nc && is_p && np==1 && nn==0 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = -32;
205  else if (is_nubar && is_nc && is_n && np==1 && nn==0 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = -33;
206  else if (is_nubar && is_nc && is_p && np==0 && nn==1 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = -34;
207 
208  //
209  // single eta from res
210  //
211 
212  else if (is_res && is_nu && is_cc && is_n && np==1 && nn==0 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = 22;
213  else if (is_res && is_nu && is_nc && is_n && np==0 && nn==1 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = 42;
214  else if (is_res && is_nu && is_nc && is_p && np==1 && nn==0 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = 43;
215 
216  else if (is_res && is_nubar && is_cc && is_p && np==0 && nn==1 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = -22;
217  else if (is_res && is_nubar && is_nc && is_n && np==0 && nn==1 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = -42;
218  else if (is_res && is_nubar && is_nc && is_p && np==1 && nn==0 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = -43;
219 
220  //
221  // single K from res (dS=0)
222  //
223 
224  else if (is_res && is_nu && is_cc && is_n && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = 23;
225  else if (is_res && is_nu && is_nc && is_n && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = 44;
226  else if (is_res && is_nu && is_nc && is_p && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = 45;
227 
228  else if (is_res && is_nubar && is_cc && is_p && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = -23;
229  else if (is_res && is_nubar && is_nc && is_n && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = -44;
230  else if (is_res && is_nubar && is_nc && is_p && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = -45;
231 
232  //
233  // single K from AtharSingleKaon (dS=1)
234  //
235  //Those modes are assigned but not used (xsec=0) in NEUT
236  else if (is_ask && is_nu && is_cc && is_n && nn==1 && np==0 && nKp==1 && neKL==1) evtype = 18;
237  else if (is_ask && is_nu && is_cc && is_n && nn==0 && np==1 && nK0==1 && neKL==1) evtype = 19;
238  else if (is_ask && is_nu && is_cc && is_p && nn==0 && np==1 && nKp==1 && neKL==1) evtype = 20;
239 
240 
241  // antineutrino modes not yet implemented
242  //else if (is_ask && is_nubar && is_cc && is_n && nn==1 && np==0 && nKp==1 && neKL==1) evtype = -18;
243  //else if (is_ask && is_nubar && is_cc && is_n && nn==0 && np==1 && nK0==1 && neKL==1) evtype = -19;
244  //else if (is_ask && is_nubar && is_cc && is_p && nn==0 && np==1 && nKp==1 && neKL==1) evtype = -20;
245 
246  //
247  // multi-pi (res or dis (W<2GeV)
248  //
249 
250  else if (is_nu && is_cc && npi>1) evtype = 21;
251  else if (is_nu && is_nc && npi>1) evtype = 41;
252  else if (is_nubar && is_cc && npi>1) evtype = -21;
253  else if (is_nubar && is_nc && npi>1) evtype = -41;
254 
255  //
256  // rare final state for RES or low-W (<2GeV) DIS events
257  // (eg K0\bar{K0} final states, N0(1720) -> Sigma- K+ res decays, etc)
258  // bundled-in with multi-pi
259  //
260  else {
261  LOG("GHepUtils", pWARN)
262  << "Rare RES/low-W DIS final state: Bundled-in with multi-pi events";
263 
264  if (is_nu && is_cc) evtype = 21;
265  else if (is_nu && is_nc) evtype = 41;
266  else if (is_nubar && is_cc) evtype = -21;
267  else if (is_nubar && is_nc) evtype = -41;
268  }
269  }
270 
271  // Weak diffractive processes
272  else if ( is_diff && is_cc ) {
273  if ( is_nu ) evtype = 15;
274  else if ( is_nubar ) evtype = -15;
275  }
276  else if ( is_diff && is_nc ) {
277  if ( is_nu ) evtype = 35;
278  else if ( is_nubar ) evtype = -35;
279  }
280 
281  return evtype;
282 }
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
genie::kPdgLambda
const int kPdgLambda
Definition: PDGCodes.h:85
train.init
init
Definition: train.py:42
genie::kPdgK0
const int kPdgK0
Definition: PDGCodes.h:174
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::kPdgKM
const int kPdgKM
Definition: PDGCodes.h:173
genie::kPdgGamma
const int kPdgGamma
Definition: PDGCodes.h:189
genie::pdg::IsAntiNeutrino
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:115
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:172
genie::kPdgEta
const int kPdgEta
Definition: PDGCodes.h:161
test.p
p
Definition: test.py:223
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
genie::kPdgAntiK0
const int kPdgAntiK0
Definition: PDGCodes.h:175
pWARN
#define pWARN
Definition: Messenger.h:60
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:159
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::kPdgAntiLambda
const int kPdgAntiLambda
Definition: PDGCodes.h:86
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
event
Event finding and building.
Definition: EventCheater_module.cc:32
genie::GHepStatus_t
enum genie::EGHepStatus GHepStatus_t
int genie::utils::ghep::NuanceReactionCode ( const GHepRecord *  evrec)

Definition at line 284 of file GHepUtils.cxx.

285 {
286 // Josh Spitz, Costas Andreopoulos
287 //
288  if(!event) {
289  LOG("GHepUtils", pWARN) << "Null event!";
290  return 0;
291  }
292 
293  int evtype = 0;
294 
295  Interaction * interaction = event->Summary();
296 
297  const ProcessInfo & proc = interaction->ProcInfo();
298  const InitialState & init = interaction->InitState();
299  if (proc.IsQuasiElastic() && proc.IsWeakCC()) evtype = 1;
300  else if (proc.IsQuasiElastic() && proc.IsWeakNC()) evtype = 2;
301  else if (proc.IsDeepInelastic() && proc.IsWeakCC()) evtype = 91;
302  else if (proc.IsDeepInelastic() && proc.IsWeakNC()) evtype = 92;
303  else if (proc.IsCoherentProduction() && proc.IsWeakNC()) evtype = 96;
304  else if (proc.IsCoherentProduction() && proc.IsWeakCC()) evtype = 97;
305  else if (proc.IsNuElectronElastic()) evtype = 98;
306  else if (proc.IsInverseMuDecay()) evtype = 99;
307  else if (proc.IsResonant()) {
308  int nn=0, np=0, npi0=0, npip=0, npim=0;
309  bool nuclear_target = init.Tgt().IsNucleus();
310  GHepStatus_t matched_ist = (nuclear_target) ?
311  kIStHadronInTheNucleus : kIStStableFinalState;
312 
313  TIter event_iter(event);
314  GHepParticle * p = 0;
315 
316  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) )
317  {
318  GHepStatus_t ghep_ist = (GHepStatus_t) p->Status();
319  if(ghep_ist != matched_ist) continue;
320 
321  int ghep_pdgc = p->Pdg();
322  if(ghep_pdgc == kPdgProton ) np++;
323  if(ghep_pdgc == kPdgNeutron) nn++;
324  if(ghep_pdgc == kPdgPi0) npi0++;
325  if(ghep_pdgc == kPdgPiP) npip++;
326  if(ghep_pdgc == kPdgPiM) npim++;
327  }
328  if(proc.IsWeakCC() && init.IsNuP()) {
329  // v p -> l- p pi+
330  if(np==1 && nn==0 && npip==1 && npi0==0 && npim==0) evtype = 3;
331  }
332  if(proc.IsWeakCC() && init.IsNuN()) {
333  // v n -> l- p pi0
334  if(np==1 && nn==0 && npip==0 && npi0==1 && npim==0) evtype = 4;
335  // v n -> l- n pi+
336  if(np==0 && nn==1 && npip==1 && npi0==0 && npim==0) evtype = 5;
337  }
338  if(proc.IsWeakNC() && init.IsNuP()) {
339  // v p -> v p pi0
340  if(np==1 && nn==0 && npip==0 && npi0==1 && npim==0) evtype = 6;
341  // v p -> v n pi+
342  if(np==0 && nn==1 && npip==1 && npi0==0 && npim==0) evtype = 7;
343  }
344  if(proc.IsWeakNC() && init.IsNuN()) {
345  // v n -> v n pi0
346  if(np==0 && nn==1 && npip==0 && npi0==1 && npim==0) evtype = 8;
347  // v n -> v p pi-
348  if(np==1 && nn==0 && npip==0 && npi0==0 && npim==1) evtype = 9;
349  }
350  if(proc.IsWeakCC() && init.IsNuBarN()) {
351  // \bar{v} n -> l+ n pi-
352  if(np==1 && nn==0 && npip==1 && npi0==0 && npim==0) evtype = 10;
353  }
354  if(proc.IsWeakCC() && init.IsNuBarP()) {
355  // \bar{v} p -> l+ n pi0
356  if(np==1 && nn==0 && npip==0 && npi0==1 && npim==0) evtype = 11;
357  // \bar{v} p -> l+ p pi-
358  if(np==0 && nn==1 && npip==1 && npi0==0 && npim==0) evtype = 12;
359  }
360  if(proc.IsWeakNC() && init.IsNuBarP()) {
361  // \bar{v} p -> \bar{v} p pi0
362  if(np==1 && nn==0 && npip==0 && npi0==1 && npim==0) evtype = 13;
363  // \bar{v} p -> \bar{v} n pi+
364  if(np==0 && nn==1 && npip==1 && npi0==0 && npim==0) evtype = 14;
365  }
366  if(proc.IsWeakNC() && init.IsNuBarN()) {
367  // \bar{v} n -> \bar{v} n pi0
368  if(np==0 && nn==1 && npip==0 && npi0==1 && npim==0) evtype = 15;
369  // \bar{v} n -> \bar{v} p pi-
370  if(np==1 && nn==0 && npip==0 && npi0==0 && npim==1) evtype = 16;
371  }
372  }
373 
374  return evtype;
375 }
train.init
init
Definition: train.py:42
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
test.p
p
Definition: test.py:223
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
pWARN
#define pWARN
Definition: Messenger.h:60
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:159
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
event
Event finding and building.
Definition: EventCheater_module.cc:32
genie::GHepStatus_t
enum genie::EGHepStatus GHepStatus_t
Generated by   doxygen 1.8.11