gEvGenHadronNucleus.cxx File Reference

#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"

Go to the source code of this file.

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.

{
// 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.

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

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

Definition at line 397 of file gEvGenHadronNucleus.cxx.

{
  LOG("gevgen_hadron", pINFO) << "Parsing command line arguments";

  // Common run options.
  RunOpt::Instance()->ReadFromCommandLine(argc,argv);

  // Parse run options for this app

  CmdLnArgParser parser(argc,argv);

  // number of events
  if( parser.OptionExists('n') ) {
    LOG("gevgen_hadron", pINFO) << "Reading number of events to generate";
    gOptNevents = parser.ArgAsInt('n');
  } else {
    LOG("gevgen_hadron", pINFO)
       << "Unspecified number of events to generate - Using default";
    gOptNevents = kDefOptNevents;
  }

  // run number
  if( parser.OptionExists('r') ) {
    LOG("gevgen_hadron", pINFO) << "Reading MC run number";
    gOptRunNu = parser.ArgAsLong('r');
  } else {
    LOG("gevgen_hadron", pINFO) << "Unspecified run number - Using default";
    gOptRunNu = kDefOptRunNu;
  }

  // incoming hadron PDG code
  if( parser.OptionExists('p') ) {
    LOG("gevgen_hadron", pINFO) << "Reading rescattering particle PDG code";
    gOptProbePdgCode = parser.ArgAsInt('p');
  } else {
    LOG("gevgen_hadron", pFATAL) << "Unspecified PDG code - Exiting";
    PrintSyntax();
    gAbortingInErr = true;
    exit(1);
  }

  // target PDG code
  if( parser.OptionExists('t') ) {
    LOG("gevgen_hadron", pINFO) << "Reading target PDG code";
    gOptTgtPdgCode = parser.ArgAsInt('t');
  } else {
    LOG("gevgen_hadron", pFATAL) << "Unspecified target PDG code - Exiting";
    PrintSyntax();
    gAbortingInErr = true;
    exit(1);
  }

  // target PDG code
  if( parser.OptionExists('m') ) {
    LOG("gevgen_hadron", pINFO) << "Reading mode";
    gOptMode = parser.ArgAsString('m');
  } else {
    LOG("gevgen_hadron", pFATAL) << "Unspecified mode - Using default";
    gOptMode = kDefOptMode;
  }

  // flux functional form or flux file
  if( parser.OptionExists('f') ) {
    LOG("gevgen_hadron", pINFO) << "Reading hadron's kinetic energy spectrum";
    gOptFlux = parser.ArgAsString('f');
    gOptUsingFlux = true;
  }

  // incoming hadron kinetic energy (or kinetic energy range, if using flux)
  if( parser.OptionExists('k') ) {
    LOG("gevgen_hadron", pINFO) << "Reading probe kinetic energy";
    string ke = parser.ArgAsString('k');
    // is it just a value or a range (comma separated set of values)
    if(ke.find(",") != string::npos) {
       // split the comma separated list
       vector<string> kerange = utils::str::Split(ke, ",");
       assert(kerange.size() == 2);
       double kemin = atof(kerange[0].c_str());
       double kemax = atof(kerange[1].c_str());
       assert(kemax>kemin && kemin>0);
       gOptProbeKE    = -1;
       gOptProbeKEmin = kemin;
       gOptProbeKEmax = kemax;
       // if no flux was specified, generate uniformly within given range
       if(!gOptUsingFlux) {
          gOptFlux = "1";
          gOptUsingFlux = true;
       }
    } else {
       gOptProbeKE    = atof(ke.c_str());
       gOptProbeKEmin = -1;
       gOptProbeKEmax = -1;
       if(gOptUsingFlux) {
          LOG("gevgen_hadron", pFATAL)
            << "You specified an input flux without a kinetic energy range";
          PrintSyntax();
          gAbortingInErr = true;
          exit(1);
       }
    }
  } else {
    LOG("gevgen_hadron", pFATAL) << "Unspecified kinetic energy - Exiting";
    PrintSyntax();
    gAbortingInErr = true;
    exit(1);
  }

  // event file prefix
  if( parser.OptionExists('o') ) {
    LOG("gevgen_hadron", pINFO) << "Reading the event filename prefix";
    gOptEvFilePrefix = parser.ArgAsString('o');
  } else {
    LOG("gevgen_hadron", pDEBUG)
      << "Will set the default event filename prefix";
    gOptEvFilePrefix = kDefOptEvFilePrefix;
  } //-o

  // INTRANUKE mode
  if( parser.OptionExists('m') ) {
    LOG("gevgen_hadron", pINFO) << "Reading mode";
    gOptMode = parser.ArgAsString('m');
  } else {
    LOG("gevgen_hadron", pDEBUG)
       << "Unspecified mode - Using default";
    gOptMode = kDefOptMode;
  }

  // random number seed
  if( parser.OptionExists("seed") ) {
    LOG("gevgen_hadron", pINFO) << "Reading random number seed";
    gOptRanSeed = parser.ArgAsLong("seed");
  } else {
    LOG("gevgen_hadron", pINFO) << "Unspecified random number seed - Using default";
    gOptRanSeed = -1;
  }


  LOG("gevgen_hadron", pNOTICE)
     << "\n"
     << utils::print::PrintFramedMesg("gevgen_hadron job configuration");

  LOG("gevgen_hadron", pNOTICE) << "MC Run Number      = " << gOptRunNu;
  LOG("gevgen_hadron", pNOTICE) << "Random number seed = " << gOptRanSeed;
  LOG("gevgen_hadron", pNOTICE) << "Mode               = " << gOptMode;
  LOG("gevgen_hadron", pNOTICE) << "Number of events   = " << gOptNevents;
  LOG("gevgen_hadron", pNOTICE) << "Probe PDG code     = " << gOptProbePdgCode;
  LOG("gevgen_hadron", pNOTICE) << "Target PDG code    = " << gOptTgtPdgCode;
  if(gOptProbeKEmin<0 && gOptProbeKEmax<0) {
    LOG("gevgen_hadron", pNOTICE)
        << "Hadron input KE    = " << gOptProbeKE;
  } else {
    LOG("gevgen_hadron", pNOTICE)
        << "Hadron input KE range = ["
        << gOptProbeKEmin << ", " << gOptProbeKEmax << "]";
  }
  if(gOptUsingFlux) {
    LOG("gevgen_hadron", pNOTICE)
        << "Input flux            = "
        << gOptFlux;
  }

  LOG("gevgen_hadron", pNOTICE) << "\n";
  LOG("gevgen_hadron", pNOTICE) << *RunOpt::Instance();
}

const EventRecordVisitorI * GetIntranuke

(

void

)

Definition at line 235 of file gEvGenHadronNucleus.cxx.

{
// 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.

{
// 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.

{
  // 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;
}

void PrintSyntax

(

void

)

Definition at line 562 of file gEvGenHadronNucleus.cxx.

{
  LOG("gevgen_hadron", pNOTICE)
    << "\n\n"
    << "Syntax:" << "\n"
    << "   gevgen_hadron [-r run] [-n nev] -p hadron_pdg -t tgt_pdg -k KE [-m mode] "
    << "                 [-f flux] "
    << "                 [--seed random_number_seed]"
    << "                 [--message-thresholds xml_file]"
    << "                 [--event-record-print-level level]"
    << "                 [--mc-job-status-refresh-rate rate]"
    << "\n";
}

Variable Documentation

string gOptEvFilePrefix

Definition at line 152 of file gEvGenHadronNucleus.cxx.

string gOptFlux

Definition at line 151 of file gEvGenHadronNucleus.cxx.

string gOptMode

Definition at line 143 of file gEvGenHadronNucleus.cxx.

int gOptNevents

Definition at line 145 of file gEvGenHadronNucleus.cxx.

double gOptProbeKE

Definition at line 148 of file gEvGenHadronNucleus.cxx.

double gOptProbeKEmax

Definition at line 150 of file gEvGenHadronNucleus.cxx.

double gOptProbeKEmin

Definition at line 149 of file gEvGenHadronNucleus.cxx.

int gOptProbePdgCode

Definition at line 146 of file gEvGenHadronNucleus.cxx.

long int gOptRanSeed

Definition at line 154 of file gEvGenHadronNucleus.cxx.

Long_t gOptRunNu

Definition at line 144 of file gEvGenHadronNucleus.cxx.

int gOptTgtPdgCode

Definition at line 147 of file gEvGenHadronNucleus.cxx.

bool gOptUsingFlux =false

Definition at line 153 of file gEvGenHadronNucleus.cxx.

TH1D* gSpectrum = 0

Definition at line 156 of file gEvGenHadronNucleus.cxx.

string kDefOptEvFilePrefix = "gntp.inuke"

Definition at line 139 of file gEvGenHadronNucleus.cxx.

string kDefOptMode = "hA"

Definition at line 140 of file gEvGenHadronNucleus.cxx.

int kDefOptNevents = 10000

Definition at line 137 of file gEvGenHadronNucleus.cxx.

Long_t kDefOptRunNu = 0

Definition at line 138 of file gEvGenHadronNucleus.cxx.