GENIEGenerator
Loading...
Searching...
No Matches
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
bool gOptWriteOutput = false
 write out hadron cross sections
string gOptOutputFilename = "gevgen_hadron_xsection.txt"

Function Documentation

◆ FindhAFate()

INukeFateHA_t FindhAFate ( const GHepRecord * evrec)

Definition at line 387 of file gtestINukeHadroXSec.cxx.

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}
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
#define pWARN
Definition Messenger.h:60
STDHEP-like event record entry that can fit a particle or a nucleus.
int Pdg(void) const
double Px(void) const
Get Px.
double Pz(void) const
Get Pz.
double Py(void) const
Get Py.
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
GHepStatus_t Status(void) const
virtual GHepParticle * Probe(void) const
const double e
const int kPdgPiM
Definition PDGCodes.h:159
@ kIStStableFinalState
Definition GHepStatus.h:30
@ kIStUndefined
Definition GHepStatus.h:28
const int kPdgProton
Definition PDGCodes.h:81
const int kPdgPi0
Definition PDGCodes.h:160
const int kPdgKP
Definition PDGCodes.h:172
const int kPdgNeutron
Definition PDGCodes.h:83
enum genie::EGHepStatus GHepStatus_t
@ kIHAFtNoInteraction
@ kIHAFtUndefined
const int kPdgKM
Definition PDGCodes.h:173
const int kPdgPiP
Definition PDGCodes.h:158
const int kPdgGamma
Definition PDGCodes.h:189
const int kPdgK0
Definition PDGCodes.h:174

References e, 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().

◆ GetCommandLineArgs()

void GetCommandLineArgs ( int argc,
char ** argv )

Definition at line 601 of file gtestINukeHadroXSec.cxx.

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}
#define pNOTICE
Definition Messenger.h:61
#define pINFO
Definition Messenger.h:62
#define pFATAL
Definition Messenger.h:56
Command line argument parser.
string gOptInpFilename
Definition gEvDump.cxx:76
string gOptOutputFilename
Definition gXSecComp.cxx:98
bool gOptWriteOutput
write out hadron cross sections
void PrintSyntax(void)
bool gAbortingInErr
Definition Messenger.cxx:34

References genie::CmdLnArgParser::ArgAsString(), genie::gAbortingInErr, gOptInpFilename, gOptOutputFilename, gOptWriteOutput, LOG, genie::CmdLnArgParser::OptionExists(), pFATAL, pINFO, pNOTICE, and PrintSyntax().

Referenced by main().

◆ main()

int main ( int argc,
char ** argv )

Definition at line 67 of file gtestINukeHadroXSec.cxx.

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}
#define pERROR
Definition Messenger.h:59
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition EventRecord.h:37
GENIE's GHEP MC event record.
Definition GHepRecord.h:45
static string AsString(INukeFateHN_t fate)
MINOS-style ntuple record. Each such ntuple record holds a generated EventRecord object....
EventRecord * event
event
void Clear(Option_t *opt="")
MINOS-style Ntuple Class to hold an output MC Tree Header.
Singleton class to load & serve a TDatabasePDG.
Definition PDGLibrary.h:36
static PDGLibrary * Instance(void)
TParticlePDG * Find(int pdgc, bool must_exist=true)
void GetCommandLineArgs(int argc, char **argv)
INukeFateHA_t FindhAFate(const GHepRecord *evrec)
bool file_exists(const std::string &file_name)
Returns true if a given file exists and is accessible, or false otherwise.
bool IsNucleon(int pdgc)
Definition PDGUtils.cxx:346
int IonPdgCodeToZ(int pdgc)
Definition PDGUtils.cxx:55
bool IsPion(int pdgc)
Definition PDGUtils.cxx:326
int IonPdgCodeToA(int pdgc)
Definition PDGUtils.cxx:63
static constexpr double mb
Definition Units.h:79
static constexpr double fm2
Definition Units.h:76
enum genie::EINukeFateHA_t INukeFateHA_t

References 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::kPdgGamma, LOG, genie::units::mb, pERROR, pNOTICE, and pWARN.

◆ PrintSyntax()

void PrintSyntax ( void )

Definition at line 627 of file gtestINukeHadroXSec.cxx.

628{
629 LOG("gtestINukeHadroXSec", pNOTICE)
630 << "\n\n"
631 << "Syntax:" << "\n"
632 << " gtestINukeHadroXSec -f event_file [-w]"
633 << "\n";
634}

References LOG, and pNOTICE.

Referenced by GetCommandLineArgs().

Variable Documentation

◆ gOptInpFilename

string gOptInpFilename = ""

input event file

Definition at line 62 of file gtestINukeHadroXSec.cxx.

◆ gOptOutputFilename

string gOptOutputFilename = "gevgen_hadron_xsection.txt"

Definition at line 64 of file gtestINukeHadroXSec.cxx.

◆ gOptWriteOutput

bool gOptWriteOutput = false

write out hadron cross sections

Definition at line 63 of file gtestINukeHadroXSec.cxx.

Referenced by GetCommandLineArgs(), and main().