gtestINukeHadroXSec.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gtestINukeHadroXSec
5 
6 \brief INTRANUKE test program. Reads-in a hadron-nucleus event file in GHEP
7  format and prints-out hadron-nucleus cross sections.
8 
9  Syntax :
10  gtestINukeHadroXSec -f input_ghep_file [-w]
11 
12  -f : Input file
13  -w : If set, writes computed hadron cross sections to a text file
14 
15 \authors Aaron Meyer and Steve Dytman
16 
17 \created July 26, 2010
18 
19 \cpright Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
20  All rights reserved.
21  For the licensing terms see $GENIE/USER_LICENSE.
22 */
23 //____________________________________________________________________________
24 
25 #include <iomanip>
26 #include <fstream>
27 #include <iostream>
28 #include <string>
29 #include <cstdlib>
30 
31 #include <TSystem.h>
32 #include <TFile.h>
33 #include <TTree.h>
34 #include <TMath.h>
35 
49 
50 using std::endl;
51 using std::setw;
52 using std::setfill;
53 
54 using namespace genie;
55 
56 void GetCommandLineArgs (int argc, char ** argv);
57 void PrintSyntax (void);
58 
59 INukeFateHA_t FindhAFate(const GHepRecord * evrec);
60 
61 // command line options
62 string gOptInpFilename = ""; ///< input event file
63 bool gOptWriteOutput = false; ///< write out hadron cross sections
64 string gOptOutputFilename = "gevgen_hadron_xsection.txt";
65 
66 //____________________________________________________________________________
67 int main(int argc, char ** argv)
68 {
69  // parse command line arguments
70  GetCommandLineArgs(argc,argv);
71 
72  //
73  // initialize
74  //
75 
76  const int nfates = 9; // total number of possible fates
77  int countfate [nfates]; // total no. of events with given fate
78  double sigma [nfates]; // cross sections
79  double sigma_err [nfates]; // cross section errors
80  string fatestr [nfates] = " ";
81  INukeFateHA_t fatetype [nfates];
82 
83  fatetype[0] = kIHAFtUndefined;
84  fatetype[1] = kIHAFtNoInteraction;
85  fatetype[2] = kIHAFtCEx;
86  fatetype[3] = kIHAFtElas;
87  fatetype[4] = kIHAFtInelas;
88  fatetype[5] = kIHAFtAbs;
89  fatetype[6] = kIHAFtKo;
90  fatetype[7] = kIHAFtPiProd;
91  fatetype[8] = kIHAFtDCEx;
92 
93  for (int k=0; k<nfates; k++) {
94  countfate[k] = 0;
95  sigma [k] = 0.;
96  sigma_err[k] = 0.;
97  fatestr [k] = INukeHadroFates::AsString(fatetype[k]);
98  }
99 
100  // event sample info (to be extracted from 1st event)
101  int nev = 0;
102  int probe_pdg = 0;
103  int target_pdg = 0;
104  int displayno = 100;
105  double kin_energy = 0.;
106 
107  //
108  // open the input ROOT file and get the event tree
109  //
110 
111  TTree * tree = 0;
112  TTree * ginuke = 0;
113  NtpMCTreeHeader * thdr = 0;
114 
115  TFile file(gOptInpFilename.c_str(),"READ");
116 
117  tree = dynamic_cast <TTree *> ( file.Get("gtree") );
118  ginuke = dynamic_cast <TTree *> ( file.Get("ginuke") );
119  thdr = dynamic_cast <NtpMCTreeHeader *> ( file.Get("header") );
120 
121  /*if(!tree) {
122  LOG("gtestINukeHadroXSec", pERROR)
123  << "No event tree found in input file!";
124  return 1;
125  }*/
126 
127  if (tree) {
128 
129  NtpMCEventRecord * mcrec = 0;
130  tree->SetBranchAddress("gmcrec", &mcrec);
131 
132  // int nev = (int) tree->GetEntries();
133  nev = (int) tree->GetEntries();
134  LOG("gtestINukeHadroXSec", pNOTICE)
135  << "Processing " << nev << " events";
136 
137  //
138  // event loop
139  //
140 
141  for(int ievent = 0; ievent < nev; ievent++) {
142 
143  // get next tree entry
144  tree->GetEntry(ievent);
145 
146  // get the corresponding GENIE event
147  EventRecord & event = *(mcrec->event);
148 
149  // extract info for the event sample
150  if(ievent==0) {
151  kin_energy = event.Particle(0)->KinE();
152  probe_pdg = event.Particle(0)->Pdg();
153  target_pdg = event.Particle(1)->Pdg();
154  }
155 
156  // analyze
157  const GHepRecord * grec = dynamic_cast<const GHepRecord *> (&event);
158  INukeFateHA_t fate = FindhAFate(grec);
159  if(ievent<displayno) {
160  LOG("gtestINukeHadroXSec", pNOTICE)
161  << "fate = " << INukeHadroFates::AsString(fate);
162  }
163 
164  // We don't want the specific fate data, just the main (9) fate types
165  switch (fate){
166 
167  case 0: countfate[0]++; break;
168  case 1: countfate[1]++; break;
169  case 2: countfate[2]++; break;
170  case 3: countfate[3]++; break;
171  case 4: countfate[4]++; break;
172  case 5: countfate[5]++; break;
173  case 6: countfate[6]++; break;
174  case 13: countfate[8]++; break;
175  default:
176  if (7<=fate && fate<=12) countfate[7]++;
177  else {
178  LOG("gtestINukeHadroXSec", pWARN)
179  << "Undefined fate from FindhAFate() : " << fate;
180  }
181  break;
182  }
183 
184  // clear current mc event record
185  mcrec->Clear();
186 
187  } // end event loop
188  } // end if (tree)
189  else if ( ginuke ) {
190  // possibly a ginuke file
191 
192  LOG("gtestINukeHadroXSec", pNOTICE)
193  << "Found ginuke type file";
194 
195  nev = (int) ginuke->GetEntries();
196  LOG("gtestINukeHadroXSec", pNOTICE)
197  << "Processing " << nev << " events";
198 
199  int kmax = 250;
200  int index = 0;
201  int numh = 0;
202  int numpip = 0;
203  int numpi0 = 0;
204  int numpim = 0;
205  int pdg_had[kmax];
206  double E_had [kmax];
207  double energy = 0.0;
208 
209  ginuke->SetBranchAddress("ke", &kin_energy);
210  ginuke->SetBranchAddress("probe",&probe_pdg );
211  ginuke->SetBranchAddress("tgt", &target_pdg);
212  ginuke->SetBranchAddress("nh", &numh );
213  ginuke->SetBranchAddress("npip", &numpip );
214  ginuke->SetBranchAddress("npi0", &numpi0 );
215  ginuke->SetBranchAddress("npim", &numpim );
216  ginuke->SetBranchAddress("pdgh", &pdg_had );
217  ginuke->SetBranchAddress("Eh", &E_had );
218  ginuke->SetBranchAddress("e", &energy );
219 
220 
221  for(int ievent = 0; ievent < nev; ievent++) {
222 
223  // get next tree entry
224  ginuke->GetEntry(ievent);
225 
226  // Determine fates (as defined in Intranuke/INukeUtils.cxx/ utils::intranuke::FindhAFate())
227  if (energy==E_had[0] && numh==1) // No interaction
228  { index=1; }
229  else if (energy!=E_had[0] && numh==1) // Elastic
230  { index=3; }
231  else if ( pdg::IsPion(probe_pdg) && numpip+numpi0+numpim==0) // Absorption
232  { index=5; }
233  else if ( (pdg::IsNucleon(probe_pdg) && numpip+numpi0+numpim==0 && numh>2 )
234  || (probe_pdg==kPdgGamma && energy!=E_had[0] && numpip+numpi0+numpim==0)) // Knock-out
235  { index=6; }
236  else if ( numpip+numpi0+numpim> (pdg::IsPion(probe_pdg) ? 1 : 0) ) // Pion production
237  { index=7; }
238  else if ( numh==2 ) // Inelastic or Charge Exchange
239  {
240  if ( (pdg::IsPion(probe_pdg) && ( probe_pdg==pdg_had[0] || probe_pdg==pdg_had[1] ))
241  || pdg::IsNucleon(probe_pdg) ) index=4;
242  else index=2;
243  }
244  else //Double Charge Exchange or Undefined
245  {
246  bool undef = true;
247  if ( pdg::IsPion(probe_pdg) )
248  {
249  for (int iter = 0; iter < numh; iter++)
250  {
251  if (probe_pdg==211 && pdg_had[iter]==-211) { index=8; undef=false; }
252  else if (probe_pdg==-211 && pdg_had[iter]==211) { index=8; undef=false; }
253  }
254  }
255  if (undef) { index=0; }
256  }
257  countfate[index]++;
258  if (ievent<displayno) {
259  LOG("gtestINukeHadroXSec", pNOTICE)
260  << "fate = " << INukeHadroFates::AsString(fatetype[index]);
261  }
262 
263  }
264  } // end if (ginuke)
265  else {
266  LOG("gtestINukeHadroXSec", pERROR)
267  << "Could not read input file!";
268  return 1;
269  }
270 
271  //
272  // output section
273  //
274 
275  const double fm2tomb = units::fm2 / units::mb;
276  const double dnev = (double) nev;
277  const int NR = 3;
278  const double R0 = 1.4;
279 
280  int A = pdg::IonPdgCodeToA(target_pdg);
281  int Z = pdg::IonPdgCodeToZ(target_pdg);
282  double nuclear_radius = NR * R0 * TMath::Power(A, 1./3.); // fm
283  double area = TMath::Pi() * TMath::Power(nuclear_radius,2);
284 
285  PDGLibrary * pdglib = PDGLibrary::Instance();
286  string probe_name = pdglib->Find(probe_pdg)->GetName();
287  string target_name = pdglib->Find(target_pdg)->GetName();
288 
289  LOG("gtestINukeHadroXSec", pNOTICE)
290  << " Probe = " << probe_name
291  << ", KE = " << kin_energy << " GeV";
292  LOG("gtestINukeHadroXSec", pNOTICE)
293  << " Target = " << target_name
294  << " (Z,A) = (" << Z << ", " << A
295  << "), nuclear radius = " << nuclear_radius
296  << " fm, area = " << area << " fm**2 " << '\n';
297 
298  int cnttot = 0;
299  int nullint = countfate[1]; // no interactions
300  double sigtot = 0;
301  double sigtoterr = 0;
302  double sigtotScat = 0;
303  double sigtotAbs = 0;
304 
305  for(int k=0; k<nfates; k++) {
306  if(k!=1) {
307  cnttot += countfate[k];
308  double ratio = countfate[k]/dnev;
309  sigma[k] = fm2tomb * area * ratio;
310  sigma_err[k] = fm2tomb * area * TMath::Sqrt(ratio*(1-ratio)/dnev);
311  if(sigma_err[k]==0) {
312  sigma_err[k] = fm2tomb * area * TMath::Sqrt(countfate[k])/dnev;
313  }
314  if(countfate[k]>0) {
315  LOG("gtestINukeHadroXSec", pNOTICE)
316  << " --> " << setw(26) << fatestr[k]
317  << ": " << setw(7) << countfate[k] << " events -> "
318  << setw(7) << sigma[k] << " +- " << sigma_err[k] << " (mb)" << '\n';
319  }
320  if(k==5) {
321  sigtotAbs += sigma[k];
322  }
323  else
324  if (k!=1) {
325  sigtotScat += sigma[k];
326  }
327  }//k!=1
328  }//k koop
329 
330  sigtot = fm2tomb * area * cnttot/dnev;
331  sigtoterr = fm2tomb * area * TMath::Sqrt(cnttot)/dnev;
332 
333  double sigtot_noelas = fm2tomb * area * (cnttot-countfate[3])/dnev;
334  double sigtoterr_noelas = fm2tomb * area * TMath::Sqrt(cnttot-countfate[3])/dnev;
335 
336  double ratio_as = (sigtotScat==0) ? 0 : sigtotAbs/(double)sigtotScat;
337 
338  LOG("gtestINukeHadroXSec", pNOTICE)
339  << "\n\n --------------------------------------------------- "
340  << "\n ==> " << setw(28) << " Total: " << setw(7) << cnttot
341  << " events -> " << setw(7) << sigtot << " +- " << sigtoterr << " (mb)"
342  << "\n (-> " << setw(28) << " Hadrons escaped nucleus: "
343  << setw(7) << nullint << " events)"
344  << "\n ==> " << setw(28) << " Ratio (abs/scat) = "
345  << setw(7) << ratio_as
346  << "\n ==> " << setw(28) << " avg. num of int. = "
347  << setw(7) << cnttot/dnev
348  << "\n ==> " << setw(28) << " no interaction = "
349  << setw(7) << (dnev-cnttot)/dnev
350  << "\n ------------------------------------------------------- \n";
351 
352  if(gOptWriteOutput)
353  {
354  ifstream test_file;
355  bool file_exists=false;
356  test_file.open(gOptOutputFilename.c_str(), std::ifstream::in);
357  file_exists=test_file.is_open();
358  test_file.close();
359  ofstream xsec_file;
360  xsec_file.open(gOptOutputFilename.c_str(), std::ios::app);
361  if (!file_exists)
362  {
363  xsec_file << "#KE" << "\t" << "Undef" << "\t"
364  << "sig" << "\t" << "CEx" << "\t"
365  << "sig" << "\t" << "Elas" << "\t"
366  << "sig" << "\t" << "Inelas"<< "\t"
367  << "sig" << "\t" << "Abs" << "\t"
368  << "sig" << "\t" << "KO" << "\t"
369  << "sig" << "\t" << "PiPro" << "\t"
370  << "sig" << "\t" << "DCEx" << "\t"
371  << "sig" << "\t" << "Reac" << "\t"
372  << "sig" << "\t" << "Tot" << "\t" << "sig" << endl;
373  }
374  xsec_file << kin_energy;
375  for(int k=0; k<nfates; k++) {
376  if (k==1) continue;
377  xsec_file << "\t" << sigma[k] << "\t" << sigma_err[k];
378  }
379  xsec_file << "\t" << sigtot_noelas << "\t" << sigtoterr_noelas;
380  xsec_file << "\t" << sigtot << "\t" << sigtoterr << endl;
381  xsec_file.close();
382  }
383 
384  return 0;
385 }
386 //____________________________________________________________________________
388 {
389  // Determine the fate of an hA event
390  // Works for ghAevgen or gntpc
391  // author: S. Dytman -- July 30, 2007
392 
393  double p_KE = evrec->Probe()->KinE();
394  double p_pdg = evrec->Probe()->Pdg();
395 
396  // particle codes
398  // num of particle for numtype
399  int num[] = {0,0,0,0,0,0,0,0,0};
400  int num_t = 0;
401  int num_nu = 0;
402  int num_pi = 0;
403  int num_k = 0;
404  // max KE for numtype
405  double numKE[] = {0,0,0,0,0,0,0,0,0};
406 
408 
409  bool hasBlob = false;
410  int numFsPart = 0;
411 
412  int index = 0;
413  TObjArrayIter piter(evrec);
414  GHepParticle * p = 0;
415  GHepParticle * fs = 0;
416  GHepParticle * probe = evrec->Probe();
417  while((p=(GHepParticle *) piter.Next()))
418  {
419  status=p->Status();
420  if(status==kIStStableFinalState)
421  {
422  switch((int) p->Pdg())
423  {
424  case ((int) kPdgProton) : index = 0; break;
425  case ((int) kPdgNeutron) : index = 1; break;
426  case ((int) kPdgPiP) : index = 2; break;
427  case ((int) kPdgPiM) : index = 3; break;
428  case ((int) kPdgPi0) : index = 4; break;
429  case ((int) kPdgKP) : index = 5; break;
430  case ((int) kPdgKM) : index = 6; break;
431  case ((int) kPdgK0) : index = 7; break;
432  case ((int) kPdgGamma) : index = 8; break;
433  case (2000000002) : index = 9; hasBlob=true; break;
434  default: index = 9; break;
435  }
436 
437  if(index!=9)
438  {
439  if(numFsPart==0) fs=p;
440  numFsPart++;
441  num[index]++;
442  if(p->KinE() > numKE[index]) numKE[index] = p->KinE();
443  }
444  }
445  }
446 
447  if(numFsPart==1)
448  {
449  double dE = TMath::Abs( probe-> E() - fs-> E() );
450  double dPz = TMath::Abs( probe->Pz() - fs->Pz() );
451  double dPy = TMath::Abs( probe->Py() - fs->Py() );
452  double dPx = TMath::Abs( probe->Px() - fs->Px() );
453 
454  if (dE < 1e-15 && dPz < 1e-15 && dPy < 1e-15 && dPx < 1e-15) return kIHAFtNoInteraction;
455  }
456 
457  num_t = num[0]+num[1]+num[2]+num[3]+num[4]+num[5]+num[6]+num[7];
458  num_nu = num[0]+num[1];
459  num_pi = num[2]+num[3]+num[4];
460  num_k = num[5]+num[6]+num[7];
461 
462  if(num_pi>((p_pdg==kPdgPiP || p_pdg==kPdgPiM || p_pdg==kPdgPi0)?(1):(0)))
463  {
464  /* if(num[3]==10 && num[4]==0) return kIHAFtNPip; //fix later
465  else if(num[4]==10) return kIHAFtNPipPi0; //fix later
466  else if(num[4]>0) return kIHAFtInclPi0;
467  else if(num[2]>0) return kIHAFtInclPip;
468  else if(num[3]>0) return kIHAFtInclPim;
469  else */
470  return kIHAFtPiProd;
471  }
472  else if(num_pi<((p_pdg==kPdgPiP || p_pdg==kPdgPiM || p_pdg==kPdgPi0)?(1):(0)))
473  {
474  if (num[0]==1 && num[1]==1) return kIHAFtAbs;
475  else if(num[0]==2 && num[1]==0) return kIHAFtAbs;
476  else if(num[0]==2 && num[1]==1) return kIHAFtAbs;
477  else if(num[0]==1 && num[1]==2) return kIHAFtAbs;
478  else if(num[0]==2 && num[1]==2) return kIHAFtAbs;
479  else if(num[0]==3 && num[1]==2) return kIHAFtAbs;
480  else return kIHAFtAbs;
481  }
482  else if(num_k<((p_pdg==kPdgKP || p_pdg==kPdgKM || p_pdg==kPdgK0)?(1):(0)))
483  {
484  return kIHAFtAbs;
485  }
486  else
487  {
488  if(p_pdg==kPdgPiP || p_pdg==kPdgPiM || p_pdg==kPdgPi0
489  || p_pdg==kPdgKP|| p_pdg==kPdgKM|| p_pdg==kPdgK0)
490  {
491  int fs_pdg, fs_ind;
492  if (num[2]==1) { fs_pdg=kPdgPiP; fs_ind=2; }
493  else if(num[3]==1) { fs_pdg=kPdgPiM; fs_ind=3; }
494  else if(num[4]==1) { fs_pdg=kPdgPi0; fs_ind=4; }
495  else if(num[5]==1) { fs_pdg=kPdgKP; fs_ind=5; }
496  else if(num[6]==1) { fs_pdg=kPdgKM; fs_ind=6; }
497  else { fs_pdg=kPdgK0; fs_ind=7; }
498 
499  if(p_pdg==fs_pdg)
500  {
501  if(num_nu==0) return kIHAFtElas;
502  else return kIHAFtInelas;
503  }
504  else if(((p_pdg==kPdgPiP || p_pdg==kPdgPiM) && fs_ind==4) ||
505  ((fs_ind==2 || fs_ind==3) && p_pdg==kPdgPi0))
506  {
507  return kIHAFtCEx;
508  }
509  else if(((p_pdg==kPdgKP || p_pdg==kPdgKM) && fs_ind==7) ||
510  ((fs_ind==5 || fs_ind==6) && p_pdg==kPdgK0))
511  {
512  return kIHAFtCEx;
513  }
514  else if((p_pdg==kPdgPiP && fs_ind==3) ||
515  (p_pdg==kPdgPiM &&fs_ind==2))
516  {
517  return kIHAFtDCEx;
518  }
519  else if((p_pdg==kPdgKP && fs_ind==6) ||
520  (p_pdg==kPdgKM &&fs_ind==5))
521  {
522  return kIHAFtDCEx;
523  }
524  }
525  else if(p_pdg==kPdgProton || p_pdg==kPdgNeutron)
526  {
527  int fs_ind;
528  if(num[0]>=1) { fs_ind=0; }
529  else { fs_ind=1; }
530 
531  if(num_nu==1)
532  {
533  if(numtype[fs_ind]==p_pdg) return kIHAFtElas;
534  else return kIHAFtUndefined;
535  }
536  else if(num_nu==2)
537  {
538  if(numKE[1]>numKE[0]) { fs_ind=1; }
539 
540  if(numtype[fs_ind]==p_pdg)
541  {
542  //if(numKE[fs_ind]>=(.8*p_KE))
543  //{
544  // if(num[0]==1 && num[1]==1) return kIHAFtKo;
545  // else if(num[0]==2) return kIHAFtKo;
546  // else return kIHAFtKo;
547  //}
548  //else
549  return kIHAFtInelas; //fix later
550  }
551  else
552  {
553  // if(numKE[fs_ind]>=(.8*p_KE)) return kIHAFtInelas;
554  // else
555  // {
556  // if(num[fs_ind]==2)
557  // {
558  // if(num[0]==2) return kIHAFtKo;
559  // else return kIHAFtKo;
560  // }
561  // else return kIHAFtInelas;
562  // }
563  return kIHAFtInelas; //fix later
564  }
565  }
566  else if(num_nu>2)
567  {
568  if (num[0]==2 && num[1]==1) return kIHAFtKo;
569  else if(num[0]==1 && num[1]==2) return kIHAFtKo;
570  else if(num[0]==2 && num[1]==2) return kIHAFtKo;
571  else if(num[0]==3 && num[1]==2) return kIHAFtKo;
572  else return kIHAFtKo;
573  }
574  }
575  else if (p_pdg==kPdgKP || p_pdg==kPdgKM || p_pdg==kPdgK0)
576  {
577  int fs_ind;
578 
579  if (num[5]==1) fs_ind=5;
580  else if (num[6]==1) fs_ind=6;
581  else fs_ind=7; // num[7]==1
582 
583  if(numKE[fs_ind]>=(.8*p_KE)) return kIHAFtElas;
584  else return kIHAFtInelas;
585  }
586  else if (p_pdg==kPdgGamma)
587  {
588  if (num[0]==2 && num[1]==1) return kIHAFtKo;
589  else if(num[0]==1 && num[1]==2) return kIHAFtKo;
590  else if(num[0]==2 && num[1]==2) return kIHAFtKo;
591  else if(num[0]==3 && num[1]==2) return kIHAFtKo;
592  else if(num_nu < 1) return kIHAFtUndefined;
593  else return kIHAFtKo;
594  }
595  }
596 
597  LOG("Intranuke",pWARN) << "---> *** Undefined fate! ***" << "\n" << (*evrec);
598  return kIHAFtUndefined;
599 }
600 //____________________________________________________________________________
601 void GetCommandLineArgs(int argc, char ** argv)
602 {
603  LOG("gtestINukeHadroXSec", pNOTICE) << "Parsing command line arguments";
604 
605  CmdLnArgParser parser(argc,argv);
606 
607  // get input ROOT file (containing a GENIE GHEP event tree)
608  if( parser.OptionExists('f') ) {
609  LOG("gtestINukeHadroXSec", pINFO) << "Reading input filename";
610  gOptInpFilename = parser.ArgAsString('f');
611  } else {
612  LOG("gtestINukeHadroXSec", pFATAL)
613  << "Unspecified input filename - Exiting";
614  PrintSyntax();
615  gAbortingInErr = true;
616  exit(1);
617  }
618  if( parser.OptionExists('o') ) {
619  LOG("gtestINukeHadroXSec", pINFO) << "Reading output filename";
620  gOptOutputFilename = parser.ArgAsString('o');
621  }
622 
623  // write-out events?
624  gOptWriteOutput = parser.OptionExists('w');
625 }
626 //____________________________________________________________________________
627 void PrintSyntax(void)
628 {
629  LOG("gtestINukeHadroXSec", pNOTICE)
630  << "\n\n"
631  << "Syntax:" << "\n"
632  << " gtestINukeHadroXSec -f event_file [-w]"
633  << "\n";
634 }
635 //____________________________________________________________________________
bool IsPion(int pdgc)
Definition: PDGUtils.cxx:289
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:148
void GetCommandLineArgs(int argc, char **argv)
int status
Definition: fabricate.py:1613
string ArgAsString(char opt)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
bool IsNucleon(int pdgc)
Definition: PDGUtils.cxx:309
MINOS-style ntuple record. Each such ntuple record holds a generated EventRecord object. Ntuples of this type are intended for feeding GENIE events into other applications (for example the GEANT4 based MC generation framework of an experiment) if no direct interface exists.
enum genie::EGHepStatus GHepStatus_t
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
TH1 * ratio(TH1 *h1, TH1 *h2)
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:61
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
double dE
GHepStatus_t Status(void) const
Definition: GHepParticle.h:65
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
const int kPdgK0
Definition: PDGCodes.h:151
Float_t Z
Definition: plot.C:38
int Pdg(void) const
Definition: GHepParticle.h:64
const double R0
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void PrintSyntax(void)
const double NR
string gOptInpFilename
input event file
Float_t E
Definition: plot.C:20
bool gOptWriteOutput
write out hadron cross sections
const int kPdgKM
Definition: PDGCodes.h:150
const int kPdgGamma
Definition: PDGCodes.h:166
double energy
Definition: plottest35.C:25
const int kPdgKP
Definition: PDGCodes.h:149
MINOS-style Ntuple Class to hold an output MC Tree Header.
const int kPdgPiP
Definition: PDGCodes.h:135
const int kPdgPi0
Definition: PDGCodes.h:137
int main(int argc, char **argv)
#define pINFO
Definition: Messenger.h:63
double sigma(TH1F *hist, double percentile)
#define pWARN
Definition: Messenger.h:61
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
string gOptOutputFilename
static const double A
Definition: Units.h:82
ifstream in
Definition: comparison.C:7
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
static string AsString(INukeFateHN_t fate)
int num
Definition: f2_nu.C:119
app
Definition: demo.py:189
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
exit(0)
TFile * file
Definition: cellShifts.C:17
const int kPdgPiM
Definition: PDGCodes.h:136
static const double fm2
Definition: Units.h:84
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:53
bool file_exists(std::string const &qualified_filename)
const int kPdgProton
Definition: PDGCodes.h:65
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:62
Float_t e
Definition: plot.C:35
void Clear(Option_t *opt="")
static const double mb
Definition: Units.h:87
bool gAbortingInErr
Definition: Messenger.cxx:56
INukeFateHA_t FindhAFate(const GHepRecord *evrec)
const int kPdgNeutron
Definition: PDGCodes.h:67
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
bool OptionExists(char opt)
was option set?
EventRecord * event
event
enum genie::EINukeFateHA_t INukeFateHA_t
double Py(void) const
Get Py.
Definition: GHepParticle.h:90