gEvGenHadronNucleus.cxx

Go to the documentation of this file.

//____________________________________________________________________________
/*!

\program gevgen_hadron

\brief   Generates hadron + nucleus interactions using GENIE's INTRANUKE
         Similar to NEUGEN's pitest (S.Dytman & H.Gallagher)

         Syntax :
           gevgen_hadron [-n nev] -p probe -t tgt [-r run#] -k KE
                         [-f flux] [-o prefix] [-m mode]
                         [--seed random_number_seed]
                         [--message-thresholds xml_file]
                         [--event-record-print-level level]
                         [--mc-job-status-refresh-rate  rate]

         Options :
           [] Denotes an optional argument
           -n
              Specifies the number of events to generate (default: 10000)
           -p
              Specifies the incoming hadron PDG code
           -t
              Specifies the nuclear target PDG code (10LZZZAAAI)
           -r
              Specifies the MC run number (default: 0)
           -k
              Specifies the incoming hadron's kinetic energy (in GeV).
              The same option can be use to specify a kinetic energy range as
              a comma-separated set of numbers (eg 0.1,1.2).
              If no flux is specified then the hadrons will be fired with a
              uniform kinetic energy distribution within the specified range.
              If a kinetic energy spectrum is supplied then the hadrons will be
              fired with a kinetic energy following the input spectrum within
              the specified range.
           -f
              Specifies the incoming hadron's kinetic energy spectrum -
              it can be either a function, eg 'x*x+4*exp(-x)' or a text file
              containing 2 columns corresponding to (kinetic energy {GeV}, 'flux').
           -o
              Output filename prefix
           -m
              INTRANUKE mode <hA, hN> (default: hA)
           --seed
              Random number seed.
           --message-thresholds
              Allows users to customize the message stream thresholds.
              The thresholds are specified using an XML file.
              See $GENIE/config/Messenger.xml for the XML schema.
           --event-record-print-level
              Allows users to set the level of information shown when the event
              record is printed in the screen. See GHepRecord::Print().
           --mc-job-status-refresh-rate
              Allows users to customize the refresh rate of the status file.

         Examples:

         (1) Generate 100k pi^{+}+Fe56 events with a pi^{+} kinetic energy
             of 165 MeV:
             % gevgen_hadron -n 100000 -p 211 -t 1000260560 -k 0.165

         (2) Generate 100k pi^{+}+Fe56 events with the pi^{+} kinetic energy
             distributed uniformly in the [165 MeV, 1200 MeV] range:
             % gevgen_hadron -n 100000 -p 211 -t 1000260560 -k 0.165,1.200

         (3) Generate 100k pi^{+}+Fe56 events with the pi^{+} kinetic energy
             distributed as f(KE) = 1/KE in the [165 MeV, 1200 MeV] range:
             % ghAevgen -n gevgen_hadron -p 211 -t 1000260560 -k 0.165,1.200 -f '1/x'

\authors  Steve Dytman, Minsuk Kim and Aaron Meyer
          University of Pittsburgh

          Costas Andreopoulos,
          University of Liverpool & STFC Rutherford Appleton Laboratory

\version 1.3

\created May 1, 2007

\cpright Copyright (c) 2003-2020, The GENIE Collaboration
         For the full text of the license visit http://copyright.genie-mc.org
         
*/
//____________________________________________________________________________

#include <cassert>
#include <cstdlib>

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

using namespace genie;
using namespace genie::controls;

using namespace genie::utils::intranuke;

// Function prototypes
void                        GetCommandLineArgs    (int argc, char ** argv);
const EventRecordVisitorI * GetIntranuke          (void);
double                      GenProbeKineticEnergy (void);
EventRecord *               InitializeEvent       (void);
void                        BuildSpectrum         (void);
void                        PrintSyntax           (void);

// Default options
int     kDefOptNevents      = 10000;   // n-events to generate
Long_t  kDefOptRunNu        = 0;       // default run number
string  kDefOptEvFilePrefix = "gntp.inuke";  // default output file prefix
string  kDefOptMode         = "hA";    // default mode

// User-specified options:
string   gOptMode;             // mode variable
Long_t   gOptRunNu;            // run number
int      gOptNevents;          // n-events to generate
int      gOptProbePdgCode;     // probe  PDG code
int      gOptTgtPdgCode;       // target PDG code
double   gOptProbeKE;          // incoming hadron kinetic enegy (GeV) - for monoenergetic probes
double   gOptProbeKEmin;       // incoming hadron kinetic enegy (GeV) - if using flux
double   gOptProbeKEmax;       // incoming hadron kinetic enegy (GeV) - if using flux
string   gOptFlux;             // input flux (function or flux file)
string   gOptEvFilePrefix;     // event file prefix
bool     gOptUsingFlux=false;  // using kinetic energy distribution?
long int gOptRanSeed ;         // random number seed

TH1D * gSpectrum  = 0;

//____________________________________________________________________________
int main(int argc, char ** argv)
{
  // 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;
}
//____________________________________________________________________________
const EventRecordVisitorI * GetIntranuke(void)
{
// 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)
{
// 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;
}
//____________________________________________________________________________
double GenProbeKineticEnergy(void)
{
  if(gOptUsingFlux) return gSpectrum->GetRandom(); // spectrum
  else              return gOptProbeKE;            // mono-energetic
}
//____________________________________________________________________________
void BuildSpectrum(void)
{
// 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();
}
//____________________________________________________________________________
void GetCommandLineArgs(int argc, char ** argv)
{
  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();
}
//____________________________________________________________________________
void PrintSyntax(void)
{
  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";
}
//____________________________________________________________________________