GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions | Variables
gtestINukeHadroXSec.cxx File Reference
#include <iomanip>
#include <fstream>
#include <iostream>
#include <string>
#include <cstdlib>
#include <TSystem.h>
#include <TFile.h>
#include <TTree.h>
#include <TMath.h>
#include "Framework/Conventions/Units.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Physics/HadronTransport/INukeHadroFates.h"
#include "Physics/HadronTransport/INukeUtils.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Ntuple/NtpMCFormat.h"
#include "Framework/Ntuple/NtpMCTreeHeader.h"
#include "Framework/Ntuple/NtpMCEventRecord.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Utils/CmdLnArgParser.h"
Include dependency graph for gtestINukeHadroXSec.cxx:

Go to the source code of this file.

Functions

void GetCommandLineArgs (int argc, char **argv)
 
void PrintSyntax (void)
 
INukeFateHA_t FindhAFate (const GHepRecord *evrec)
 
int main (int argc, char **argv)
 

Variables

string gOptInpFilename = ""
 input event file More...
 
bool gOptWriteOutput = false
 write out hadron cross sections More...
 
string gOptOutputFilename = "gevgen_hadron_xsection.txt"
 

Function Documentation

INukeFateHA_t FindhAFate ( const GHepRecord evrec)

Definition at line 387 of file gtestINukeHadroXSec.cxx.

References e, genie::units::fs, genie::kIHAFtAbs, genie::kIHAFtCEx, genie::kIHAFtDCEx, genie::kIHAFtElas, genie::kIHAFtInelas, genie::kIHAFtKo, genie::kIHAFtNoInteraction, genie::kIHAFtPiProd, genie::kIHAFtUndefined, genie::GHepParticle::KinE(), genie::kIStStableFinalState, genie::kIStUndefined, genie::kPdgGamma, genie::kPdgK0, genie::kPdgKM, genie::kPdgKP, genie::kPdgNeutron, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, LOG, genie::GHepParticle::Pdg(), genie::GHepRecord::Probe(), pWARN, genie::GHepParticle::Px(), genie::GHepParticle::Py(), genie::GHepParticle::Pz(), and genie::GHepParticle::Status().

Referenced by main().

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 
407  GHepStatus_t status = kIStUndefined;
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 }
static constexpr double fs
Definition: Units.h:100
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:284
const int kPdgK0
Definition: PDGCodes.h:174
int Pdg(void) const
Definition: GHepParticle.h:63
const double e
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const int kPdgKM
Definition: PDGCodes.h:173
const int kPdgGamma
Definition: PDGCodes.h:189
const int kPdgKP
Definition: PDGCodes.h:172
const int kPdgPiP
Definition: PDGCodes.h:158
const int kPdgPi0
Definition: PDGCodes.h:160
#define pWARN
Definition: Messenger.h:60
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
const int kPdgPiM
Definition: PDGCodes.h:159
const int kPdgProton
Definition: PDGCodes.h:81
const int kPdgNeutron
Definition: PDGCodes.h:83
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
enum genie::EGHepStatus GHepStatus_t
double Py(void) const
Get Py.
Definition: GHepParticle.h:89
void GetCommandLineArgs ( int  argc,
char **  argv 
)
int main ( int  argc,
char **  argv 
)

Definition at line 67 of file gtestINukeHadroXSec.cxx.

References genie::units::A, genie::INukeHadroFates::AsString(), genie::NtpMCEventRecord::Clear(), genie::NtpMCEventRecord::event, genie::PDGLibrary::Find(), FindhAFate(), genie::units::fm2, GetCommandLineArgs(), gOptInpFilename, gOptOutputFilename, gOptWriteOutput, genie::PDGLibrary::Instance(), genie::pdg::IonPdgCodeToA(), genie::pdg::IonPdgCodeToZ(), genie::pdg::IsNucleon(), genie::pdg::IsPion(), genie::kIHAFtAbs, genie::kIHAFtCEx, genie::kIHAFtDCEx, genie::kIHAFtElas, genie::kIHAFtInelas, genie::kIHAFtKo, genie::kIHAFtNoInteraction, genie::kIHAFtPiProd, genie::kIHAFtUndefined, genie::GHepParticle::KinE(), genie::kPdgGamma, LOG, genie::units::mb, genie::GHepRecord::Particle(), pERROR, pNOTICE, and pWARN.

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 }
string gOptOutputFilename
Definition: gXSecComp.cxx:98
bool IsPion(int pdgc)
Definition: PDGUtils.cxx:326
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:104
#define pERROR
Definition: Messenger.h:59
bool IsNucleon(int pdgc)
Definition: PDGUtils.cxx:346
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.
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:63
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double A
Definition: Units.h:74
bool gOptWriteOutput
write out hadron cross sections
static constexpr double mb
Definition: Units.h:79
const int kPdgGamma
Definition: PDGCodes.h:189
static constexpr double fm2
Definition: Units.h:76
MINOS-style Ntuple Class to hold an output MC Tree Header.
#define pWARN
Definition: Messenger.h:60
string gOptInpFilename
Definition: gEvDump.cxx:76
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:37
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
const char * AsString(Resonance_t res)
resonance id -&gt; string
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:55
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
#define pNOTICE
Definition: Messenger.h:61
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:563
void Clear(Option_t *opt="")
INukeFateHA_t FindhAFate(const GHepRecord *evrec)
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
EventRecord * event
event
enum genie::EINukeFateHA_t INukeFateHA_t
void PrintSyntax ( void  )

Variable Documentation

string gOptInpFilename = ""

input event file

Definition at line 62 of file gtestINukeHadroXSec.cxx.

string gOptOutputFilename = "gevgen_hadron_xsection.txt"

Definition at line 64 of file gtestINukeHadroXSec.cxx.

bool gOptWriteOutput = false

write out hadron cross sections

Definition at line 63 of file gtestINukeHadroXSec.cxx.

Referenced by main().