#include <cassert>
#include <cstdlib>
#include "TSystem.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1D.h"
#include "TF1.h"
#include "Framework/Conventions/GBuild.h"
#include "Framework/Algorithm/AlgFactory.h"
#include "Framework/Conventions/Controls.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/EventGen/GMCJMonitor.h"
#include "Framework/EventGen/EventRecordVisitorI.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/GHEP/GHepRecord.h"
#include "Framework/GHEP/GHepStatus.h"
#include "Framework/Interaction/Interaction.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Ntuple/NtpWriter.h"
#include "Framework/Ntuple/NtpMCFormat.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/Numerical/Spline.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/StringUtils.h"
#include "Framework/Utils/PrintUtils.h"
#include "Framework/Utils/XSecSplineList.h"
#include "Framework/Utils/RunOpt.h"
#include "Framework/Utils/CmdLnArgParser.h"
#include "Physics/HadronTransport/INukeHadroFates.h"
#include "Physics/HadronTransport/INukeUtils.h"

Include dependency graph for gEvGenHadronNucleus.cxx:

Functions
void	GetCommandLineArgs (int argc, char **argv)

const EventRecordVisitorI *	GetIntranuke (void)

double	GenProbeKineticEnergy (void)

EventRecord *	InitializeEvent (void)

void	BuildSpectrum (void)

void	PrintSyntax (void)

int	main (int argc, char **argv)

Variables
int	kDefOptNevents = 10000

Long_t	kDefOptRunNu = 0

string	kDefOptEvFilePrefix = "gntp.inuke"

string	kDefOptMode = "hA"

string	gOptMode

Long_t	gOptRunNu

int	gOptNevents

int	gOptProbePdgCode

int	gOptTgtPdgCode

double	gOptProbeKE

double	gOptProbeKEmin

double	gOptProbeKEmax

string	gOptFlux

string	gOptEvFilePrefix

bool	gOptUsingFlux =false

long int	gOptRanSeed

TH1D *	gSpectrum = 0

Function Documentation

void BuildSpectrum ( void )

Definition at line 323 of file gEvGenHadronNucleus.cxx.

References genie::Spline::Evaluate(), genie::gAbortingInErr, gOptFlux, gOptProbeKEmax, gOptProbeKEmin, gOptUsingFlux, gSpectrum, genie::RandomGen::Instance(), genie::controls::kRjMaxIterations, LOG, pFATAL, pNOTICE, and genie::RandomGen::RndGen().

Referenced by main().

 {
 // Create kinetic energy spectrum from input function
 //
 
   if(!gOptUsingFlux) return;
 
   if(gSpectrum) {
     delete gSpectrum;
     gSpectrum = 0;
   }
 
   LOG("gevgen_hadron", pNOTICE)
       << "Generating a flux histogram ... ";
 
   int    flux_bins    = 300;
   int    flux_entries = 1000000;
   double ke_min       = gOptProbeKEmin;
   double ke_max       = gOptProbeKEmax;
   double dke          = ke_max - ke_min;
   assert(dke>0 && ke_min>=0.);
 
   // kinetic energy spectrum
   gSpectrum  = new TH1D(
     "spectrum","hadron kinetic energy spectrum", flux_bins, ke_min, ke_max);
 
   // check whether the input flux is a file or a functional form
   bool input_is_file = ! gSystem->AccessPathName(gOptFlux.c_str());
   if(input_is_file) {
     Spline * input_flux = new Spline(gOptFlux.c_str());
 
     // generate the flux hisogram from the input flux file
     int    n       = 100;
     double ke_step = (ke_max-ke_min)/(n-1);
     double ymax  = -1, ry = -1, gy = -1, ke = -1;
     for(int i=0; i<n; i++) {
       ke   = ke_min + i*ke_step;
       ymax = TMath::Max(ymax, input_flux->Evaluate(ke));
     }
     ymax *= 1.3;
 
     RandomGen * r = RandomGen::Instance();
 
     for(int ientry=0; ientry<flux_entries; ientry++) {
       bool accept = false;
       unsigned int iter=0;
       while(!accept) {
         iter++;
         if(iter > kRjMaxIterations) {
            LOG("gevgen_hadron", pFATAL)
              << "Couldn't generate a flux histogram";
            gAbortingInErr = true;
            exit(1);
         }
         ke = ke_min + dke * r->RndGen().Rndm();
         gy = ymax * r->RndGen().Rndm();
         ry = input_flux->Evaluate(ke);
         accept = gy < ry;
         if(accept) gSpectrum->Fill(ke);
       }
     }
     delete input_flux;
 
   } else {
     // generate the flux histogram from the input functional form
     TF1 *  input_func = new TF1("input_func", gOptFlux.c_str(), ke_min, ke_max);
     gSpectrum->FillRandom("input_func", flux_entries);
     delete input_func;
   }
   TFile f("./input-hadron-flux.root","recreate");
   gSpectrum->Write();
   f.Close();
 }

double GenProbeKineticEnergy ( void )

Definition at line 317 of file gEvGenHadronNucleus.cxx.

References gOptProbeKE, gOptUsingFlux, and gSpectrum.

Referenced by InitializeEvent().

 {
   if(gOptUsingFlux) return gSpectrum->GetRandom(); // spectrum
   else              return gOptProbeKE;            // mono-energetic
 }

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

const EventRecordVisitorI * GetIntranuke ( void )

Definition at line 235 of file gEvGenHadronNucleus.cxx.

References genie::gAbortingInErr, genie::AlgFactory::GetAlgorithm(), gOptMode, genie::AlgFactory::Instance(), LOG, and pFATAL.

Referenced by main().

 {
 // get the requested INTRANUKE module
 
   string sname = "";
   string sconf = "";
 
   if        ( gOptMode.compare("hA") == 0 ) {
      sname = "genie::HAIntranuke";
      sconf = "Default";
   } else if ( gOptMode.compare("hN") == 0 ) {
      sname = "genie::HNIntranuke";
      sconf = "Default";
   } else if ( gOptMode.compare("hA2019") == 0 ) {
      sname = "genie::HAIntranuke2019";
      sconf = "Default";
   } else if ( gOptMode.compare("hN2019") == 0 ) {
      sname = "genie::HNIntranuke2019";
      sconf = "Default";
   } else if ( gOptMode.compare("hA2018") == 0 ) {
      sname = "genie::HAIntranuke2018";
      sconf = "Default";
   } else if ( gOptMode.compare("hN2018") == 0 ) {
      sname = "genie::HNIntranuke2018";
      sconf = "Default";
 #ifdef __GENIE_INCL_ENABLED__
   } else if ( gOptMode.compare("HINCL") == 0 ) {
      sname = "genie::HINCLCascadeIntranuke";
      sconf = "Default";
 #endif
 #ifdef __GENIE_GEANT4_INTERFACE_ENABLED__
   } else if ( gOptMode.compare("HG4BertCasc") == 0 ) {
      sname = "genie::HG4BertCascIntranuke";
      sconf = "Default";
 #endif
   } else {
     LOG("gevgen_hadron", pFATAL) << "Invalid Intranuke mode - Exiting";
     gAbortingInErr = true;
     exit(1);
   }
 
   AlgFactory * algf = AlgFactory::Instance();
   const EventRecordVisitorI * intranuke =
    dynamic_cast<const EventRecordVisitorI *> (algf->GetAlgorithm(sname,sconf));
   assert(intranuke);
 
   return intranuke;
 }

EventRecord * InitializeEvent ( void )

Definition at line 284 of file gEvGenHadronNucleus.cxx.

References genie::GHepRecord::AddParticle(), genie::GHepRecord::AttachSummary(), GenProbeKineticEnergy(), gOptProbePdgCode, gOptTgtPdgCode, genie::PDGLibrary::Instance(), genie::kIStInitialState, and genie::utils::res::Mass().

Referenced by main().

 {
 // Initialize event record. Inserting the probe and target particles.
 
   EventRecord * evrec = new EventRecord();
   Interaction * interaction = new Interaction;
   evrec->AttachSummary(interaction);
 
   // dummy vertex position
   TLorentzVector x4null(0.,0.,0.,0.);
 
   // incident hadron & target nucleon masses
   PDGLibrary * pdglib = PDGLibrary::Instance();
   double mh  = pdglib -> Find (gOptProbePdgCode) -> Mass();
   double M   = pdglib -> Find (gOptTgtPdgCode  ) -> Mass();
 
   // incident hadron kinetic energy
   double ke = GenProbeKineticEnergy();
 
   // form  incident hadron and target 4-momenta
   double Eh  = mh + ke;
   double pzh = TMath::Sqrt(TMath::Max(0.,Eh*Eh-mh*mh));
   TLorentzVector p4h   (0.,0.,pzh,Eh);
   TLorentzVector p4tgt (0.,0.,0., M);
 
   // insert probe and target entries
   GHepStatus_t ist = kIStInitialState;
   evrec->AddParticle(gOptProbePdgCode, ist, -1,-1,-1,-1, p4h,   x4null);
   evrec->AddParticle(gOptTgtPdgCode,   ist, -1,-1,-1,-1, p4tgt, x4null);
 
   return evrec;
 }

int main	(	int	argc,
		char **	argv
	)

Definition at line 159 of file gEvGenHadronNucleus.cxx.

References genie::NtpWriter::AddEventRecord(), BuildSpectrum(), genie::RunOpt::BuildTune(), genie::NtpWriter::CustomizeFilenamePrefix(), GetCommandLineArgs(), GetIntranuke(), gOptEvFilePrefix, gOptNevents, gOptRanSeed, gOptRunNu, gSpectrum, genie::NtpWriter::Initialize(), InitializeEvent(), genie::RunOpt::Instance(), genie::kNFGHEP, LOG, genie::utils::app_init::MesgThresholds(), pFATAL, pNOTICE, genie::EventRecordVisitorI::ProcessEventRecord(), genie::utils::app_init::RandGen(), genie::NtpWriter::Save(), genie::GHepRecord::SetPrintLevel(), and genie::GMCJMonitor::Update().

 {
   // Parse command line arguments
   GetCommandLineArgs(argc,argv);
 
   if ( ! RunOpt::Instance()->Tune() ) {
     LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
     exit(-1);
   }
   RunOpt::Instance()->BuildTune();
 
   // Init random number generator generator with user-specified seed number,
   // set user-specified mesg thresholds, set user-specified GHEP print-level
   utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
   utils::app_init::RandGen(gOptRanSeed);
   GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
 
   // Build the incident hadron kinetic energy spectrum, if required
   BuildSpectrum();
 
   LOG("gevgen_hadron",pNOTICE) << "finish setup";
 
   // Get the specified INTRANUKE model
   const EventRecordVisitorI * intranuke = GetIntranuke();
 
   // Initialize an Ntuple Writer to save GHEP records into a ROOT tree
   NtpWriter ntpw(kNFGHEP, gOptRunNu, gOptRanSeed);
   ntpw.CustomizeFilenamePrefix(gOptEvFilePrefix);
   ntpw.Initialize();
 
   // Create an MC job monitor
   GMCJMonitor mcjmonitor(gOptRunNu);
 
   LOG("gevgen_hadron",pNOTICE) << "ready to generate events";
 
   //
   // Generate events
   //
 
   int ievent = 0;
   while (ievent < gOptNevents) {
       LOG("gevgen_hadron", pNOTICE)
          << " *** Generating event............ " << ievent;
 
       // initialize
       EventRecord * evrec = InitializeEvent();
 
       // generate full h+A event
       intranuke->ProcessEventRecord(evrec);
 
       // print generated events
       LOG("gevgen_hadron", pNOTICE ) << *evrec;
 
       // add event at the output ntuple
       ntpw.AddEventRecord(ievent, evrec);
 
       // refresh the mc job monitor
       mcjmonitor.Update(ievent,evrec);
 
       ievent++;
       delete evrec;
 
   } // end loop events
 
   // Save the generated MC events
   ntpw.Save();
 
   // Clean-up
   if(gSpectrum) {
     delete gSpectrum;
     gSpectrum = 0;
   }
 
   return 0;
 }