#include <cassert>
#include <cstdlib>
#include <string>
#include <sstream>
#include <vector>
#include <map>
#include <TSystem.h>
#include <TTree.h>
#include <TFile.h>
#include <TH1D.h>
#include <TMath.h>
#include <TGeoVolume.h>
#include <TGeoShape.h>
#include <TList.h>
#include <TObject.h>
#include "Framework/Conventions/Units.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/EventGen/GFluxI.h"
#include "Framework/EventGen/GMCJDriver.h"
#include "Framework/EventGen/GMCJMonitor.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Ntuple/NtpWriter.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGCodeList.h"
#include "Framework/Ntuple/NtpMCFormat.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/Utils/XSecSplineList.h"
#include "Framework/Utils/StringUtils.h"
#include "Framework/Utils/UnitUtils.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"
#include "Framework/Utils/CmdLnArgParser.h"
#include "Framework/Utils/T2KEvGenMetaData.h"
#include "Framework/Utils/SystemUtils.h"
#include "Framework/Utils/PrintUtils.h"

Include dependency graph for gT2KEvGen.cxx:

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

void	PrintSyntax (void)

int	main (int argc, char **argv)

Variables
string	kDefOptGeomLUnits = "mm"

string	kDefOptGeomDUnits = "g_cm3"

NtpMCFormat_t	kDefOptNtpFormat = kNFGHEP

double	kDefOptFluxNorm = 1E+21

string	kDefOptEvFilePrefix = "gntp"

Long_t	gOptRunNu

bool	gOptUsingRootGeom = false

bool	gOptUsingHistFlux = false

map< int, double >	gOptTgtMix

map< int, TH1D * >	gOptFluxHst

string	gOptRootGeom

string	gOptRootGeomTopVol = ""

double	gOptGeomLUnits = 0

double	gOptGeomDUnits = 0

string	gOptExtMaxPlXml

string	gOptFluxFile

string	gOptDetectorLocation

double	gOptFluxNorm

PDGCodeList	gOptFluxNtpNuList (false)

int	gOptFluxNCycles

int	gOptNev

double	gOptPOT

bool	gOptExitAtEndOfFullFluxCycles

string	gOptEvFilePrefix

bool	gOptUseFluxProbs = false

bool	gOptSaveFluxProbsFile = false

string	gOptFluxProbFileName

string	gOptSaveFluxProbsFileName

bool	gOptRandomFluxOffset = false

long int	gOptRanSeed

string	gOptInpXSecFile

Function Documentation

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

int main	(	int	argc,
		char **	argv
	)

Definition at line 524 of file gT2KEvGen.cxx.

 {
   // Parse command line arguments
   GetCommandLineArgs(argc,argv);
 
 
   if ( ! RunOpt::Instance()->Tune() ) {
     LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
     exit(-1);
   }
   RunOpt::Instance()->BuildTune();
 
   // Initialization of random number generators, cross-section table,
   // messenger thresholds, cache file
   utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
   utils::app_init::CacheFile(RunOpt::Instance()->CacheFile());
   utils::app_init::RandGen(gOptRanSeed);
   utils::app_init::XSecTable(gOptInpXSecFile, true);
 
   // Set GHEP print level
   GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
 
   // *************************************************************************
   // * Create / configure the geometry driver
   // *************************************************************************
   GeomAnalyzerI * geom_driver = 0;
   double zmin=0, zmax=0;
 
   if(gOptUsingRootGeom) {
     //
     // *** Using a realistic root-based detector geometry description
     //
 
     // creating & configuring a root geometry driver
     geometry::ROOTGeomAnalyzer * rgeom =
             new geometry::ROOTGeomAnalyzer(gOptRootGeom);
     rgeom -> SetLengthUnits  (gOptGeomLUnits);
     rgeom -> SetDensityUnits (gOptGeomDUnits);
     // getting the bounding box dimensions, before setting topvolume,
     // along z so as to set the appropriate upstream generation surface
     // for the JPARC flux driver
     TGeoVolume * topvol = rgeom->GetGeometry()->GetTopVolume();
     TGeoShape * bounding_box = topvol->GetShape();
     bounding_box->GetAxisRange(3, zmin, zmax);
     zmin *= rgeom->LengthUnits();
     zmax *= rgeom->LengthUnits();
     // now update to the requested topvolume for use in recursive exhaust method
     rgeom -> SetTopVolName   (gOptRootGeomTopVol);
     topvol = rgeom->GetGeometry()->GetTopVolume();
     if(!topvol) {
       LOG("gevgen_t2k", pFATAL) << "Null top ROOT geometry volume!";
       exit(1);
     }
     // switch on/off volumes as requested
     if ( (gOptRootGeomTopVol[0] == '+') || (gOptRootGeomTopVol[0] == '-') ) {
       bool exhaust = (*gOptRootGeomTopVol.c_str() == '+');
       utils::geometry::RecursiveExhaust(topvol, gOptRootGeomTopVol, exhaust);
     }
 
     // casting to the GENIE geometry driver interface
     geom_driver = dynamic_cast<GeomAnalyzerI *> (rgeom);
   }
   else {
     //
     // *** Using a 'point' geometry with the specified target mix
     // *** ( = a list of targets with their corresponding weight fraction)
     //
 
     // creating & configuring a point geometry driver
     geometry::PointGeomAnalyzer * pgeom =
               new geometry::PointGeomAnalyzer(gOptTgtMix);
     // casting to the GENIE geometry driver interface
     geom_driver = dynamic_cast<GeomAnalyzerI *> (pgeom);
   }
 
   // *************************************************************************
   // * Create / configure the flux driver
   // *************************************************************************
   GFluxI * flux_driver = 0;
 
   flux::GJPARCNuFlux *   jparc_flux_driver = 0;
   flux::GCylindTH1Flux * hst_flux_driver   = 0;
 
   if(!gOptUsingHistFlux) {
     //
     // *** Using the detailed JPARC neutrino flux desription by feeding-in
     // *** the JNUBEAM flux simulation ntuples
     //
 
     // creating JPARC neutrino flux driver
     jparc_flux_driver = new flux::GJPARCNuFlux;
     // before loading the beam sim data set whether to use a random offset when looping
     if(gOptRandomFluxOffset == false) jparc_flux_driver->DisableOffset();
     // specify input JNUBEAM file & detector location
     bool beam_sim_data_success = jparc_flux_driver->LoadBeamSimData(gOptFluxFile, gOptDetectorLocation);
     if(!beam_sim_data_success) {
       LOG("gevgen_t2k", pFATAL)
         << "The flux driver has not been properly configured. Exiting";
       exit(1);
     }
     // specify JNUBEAM normalization
     jparc_flux_driver->SetFilePOT(gOptFluxNorm);
     // specify upstream generation surface
     jparc_flux_driver->SetUpstreamZ(zmin);
     // specify which flavours to consider -
     // if no neutrino list was specified then the MC job will consider all flavours
     if( gOptFluxNtpNuList.size() > 0 ) {
        jparc_flux_driver->SetFluxParticles(gOptFluxNtpNuList);
     }
 
     // casting to the GENIE flux driver interface
     flux_driver = dynamic_cast<GFluxI *> (jparc_flux_driver);
   }
   else {
     //
     // *** Using fluxes from histograms (for all specified neutrino species)
     //
 
     // creating & configuring a generic GCylindTH1Flux flux driver
     TVector3 bdir (0,0,1); // dir along +z
     TVector3 bspot(0,0,0);
     hst_flux_driver = new flux::GCylindTH1Flux;
     hst_flux_driver->SetNuDirection      (bdir);
     hst_flux_driver->SetBeamSpot         (bspot);
     hst_flux_driver->SetTransverseRadius (-1);
     map<int,TH1D*>::iterator it = gOptFluxHst.begin();
     for( ; it != gOptFluxHst.end(); ++it) {
         int    pdg_code = it->first;
         TH1D * spectrum = new TH1D(*(it->second));
         hst_flux_driver->AddEnergySpectrum(pdg_code, spectrum);
     }
     // casting to the GENIE flux driver interface
     flux_driver = dynamic_cast<GFluxI *> (hst_flux_driver);
   }
 
   // *************************************************************************
   // * Create/configure the event generation driver
   // *************************************************************************
   GMCJDriver * mcj_driver = new GMCJDriver;
   mcj_driver->SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
   mcj_driver->UseFluxDriver(flux_driver);
   mcj_driver->UseGeomAnalyzer(geom_driver);
   mcj_driver->UseMaxPathLengths(gOptExtMaxPlXml);
   // do not calculate probability scales if using pre-generated flux probs
   bool calc_prob_scales = (gOptSaveFluxProbsFile || gOptUseFluxProbs) ? false : true;
   mcj_driver->Configure(calc_prob_scales);
   mcj_driver->UseSplines();
   mcj_driver->ForceSingleProbScale();
 
   // *************************************************************************
   // * If specified use pre-calculated flux interaction probabilities instead
   // * of preselecting based on max path lengths
   // *************************************************************************
 
   if(gOptUseFluxProbs || gOptSaveFluxProbsFile){
 
     bool success = false;
 
     // set flux probs output file name
     if(gOptSaveFluxProbsFile){
       // default output name is ${FLUFILENAME}.${TOPVOL}.flxprobs.root
       string basename = gOptFluxFile.substr(gOptFluxFile.rfind("/")+1);
       string name = basename.substr(0, basename.rfind("."));
       if(gOptRootGeomTopVol.length()>0)
         name += "."+gOptRootGeomTopVol+".flxprobs.root";
       else
         name += ".master.flxprobs.root";
       // if specified override with cmd line option
       if(gOptSaveFluxProbsFileName.size()>0) name = gOptSaveFluxProbsFileName;
       // Tell the driver save pre-generated probabilities to an output file
       mcj_driver->SaveFluxProbabilities(name);
     }
 
     // Either load pre-generated flux probabilities
     if(gOptFluxProbFileName.size() > 0){
       success = mcj_driver->LoadFluxProbabilities(gOptFluxProbFileName);
     }
     // Or pre-calculate them
     else success = mcj_driver->PreCalcFluxProbabilities();
 
     if(success){
       LOG("gevgen_t2k", pNOTICE)
        << "Successfully calculated/loaded flux interaction probabilities!";
       // Print out a list of expected number of events per POT and per cycle
       // based on the pre-generated flux interaction probabilities
       map<int, double> sum_probs_map = mcj_driver->SumFluxIntProbs();
       map<int, double>::const_iterator sum_probs_it = sum_probs_map.begin();
       double ntot_per_pot = 0.0;
       double ntot_per_cycle = 0.0;
       double pscale = mcj_driver->GlobProbScale();
       double pot_1cycle = jparc_flux_driver->POT_1cycle();
       LOG("T2KProdInfo", pNOTICE) <<
           "Expected event rates based on flux interaction probabilities:";
       for(; sum_probs_it != sum_probs_map.end(); sum_probs_it++){
         double sum_probs = sum_probs_it->second;
         double nevts_per_cycle = sum_probs / pscale;  // take into account rescale
         double nevts_per_pot = sum_probs/pot_1cycle;  // == (sum_probs*pscale)/(pot_1cycle*pscale)
         ntot_per_pot += nevts_per_pot;
         ntot_per_cycle += nevts_per_cycle;
         LOG("T2KProdInfo", pNOTICE) <<
             "     PDG "<< sum_probs_it->first << ": " << nevts_per_cycle <<
             " Events/Cycle, "<< nevts_per_pot << " Events/POT";
       }
       LOG("T2KProdInfo", pNOTICE) << "      -----------";
       LOG("T2KProdInfo", pNOTICE) <<
           "     All neutrino species: " << ntot_per_cycle <<
           " Events/Cycle, "<< ntot_per_pot << " Events/POT";
       LOG("T2KProdInfo", pNOTICE) <<
           "N.B. This assumes input flux file corresponds to "<< pot_1cycle <<
           "POT, ensure this is correct if using these numbers!";
     }
     else {
       LOG("gevgen_t2k", pFATAL)
        << "Failed to calculated/loaded flux interaction probabilities!";
       return 1;
     }
 
     // Exit now if just pre-generating interaction probabilities
     if(gOptSaveFluxProbsFile){
       LOG("gevgen_t2k", pNOTICE)
        << "Will not generate events - just pre-calculating flux interaction"
        << "probabilities";
       return 0;
     }
   } // Pre-calculated flux interaction probabilities
 
   // *************************************************************************
   // * Work out number of cycles for current exposure settings
   // *************************************************************************
 
   if(!gOptUsingHistFlux) {
     // If a number of POT was requested, then work out how many flux ntuple
     // cycles are required for accumulating those statistics
     if(gOptPOT>0) {
       double fpot_1c = jparc_flux_driver->POT_1cycle(); // flux POT / cycle
       double psc     = mcj_driver->GlobProbScale();     // interaction prob. scale
       double pot_1c  = fpot_1c / psc;                   // actual POT / cycle
       int    ncycles = (int) TMath::Max(1., TMath::Ceil(gOptPOT/pot_1c));
 
       LOG("gevgen_t2k", pNOTICE)
          << " *** POT/cycle:  " << pot_1c;
       LOG("gevgen_t2k", pNOTICE)
          << " *** Requested POT will be accumulated in: "
          << ncycles << " flux ntuple cycles";
 
       jparc_flux_driver->SetNumOfCycles(ncycles);
     }
     // If a number of events was requested, then set the number of flux
     // ntuple cycles to 'infinite'
     else if(gOptNev>0) {
        jparc_flux_driver->SetNumOfCycles(0);
     }
     // Just set the number of cycles to the requested value
     else {
        jparc_flux_driver->SetNumOfCycles(gOptFluxNCycles);
     }
   }
 
   // *************************************************************************
   // * Prepare for writing the output event tree & status file
   // *************************************************************************
 
   // Initialize an Ntuple Writer to save GHEP records into a TTree
   NtpWriter ntpw(kDefOptNtpFormat, gOptRunNu, gOptRanSeed);
   ntpw.CustomizeFilenamePrefix(gOptEvFilePrefix);
   ntpw.Initialize();
 
   // Add a custom-branch at the standard GENIE event tree so that
   // info on the flux neutrino parent particle can be passed-through
   flux::GJPARCNuFluxPassThroughInfo * flux_info = 0;
   if(!gOptUsingHistFlux) {
     TBranch * flux = ntpw.EventTree()->Branch("flux",
        "genie::flux::GJPARCNuFluxPassThroughInfo", &flux_info, 32000, 1);
     assert(flux);
     flux->SetAutoDelete(kFALSE);
   }
 
   // Create a MC job monitor for a periodically updated status file
   GMCJMonitor mcjmonitor(gOptRunNu);
   mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());
 
   // *************************************************************************
   // * Event generation loop
   // *************************************************************************
 
   int ievent = 0;
   while (1)
   {
      LOG("gevgen_t2k", pNOTICE)
           << " *** Generating event............ " << ievent;
 
      // In case the required statistics was expressed as 'number of events'
      // then quit if that number has been generated
      if(ievent == gOptNev) break;
 
      // In case the required statistics was expressed as 'number of POT' and
      // the user does not want to wait till the end of the flux cycle to exit
      // the event loop, then quit if the requested POT has been generated.
      // In this case the computed POT may not be as accurate as if the program
      // was waiting for the current flux cycle to be completed.
      if(!gOptExitAtEndOfFullFluxCycles && gOptPOT>0) {
         double fpot = jparc_flux_driver->POT_curravg(); // current POT in flux file
         double psc  = mcj_driver->GlobProbScale();      // interaction prob. scale
         double pot  = fpot / psc;                       // POT for generated sample
         if(pot >= gOptPOT) break;
      }
 
      // Generate a single event using neutrinos coming from the specified flux
      // and hitting the specified geometry or target mix
      EventRecord * event = mcj_driver->GenerateEvent();
 
      // Check whether a null event was returned due to the flux driver reaching
      // the end of the input flux ntuple - exit the event generation loop
      if(!event && jparc_flux_driver->End()) {
         LOG("gevgen_t2k", pNOTICE)
           << "** The JPARC flux driver read all the input flux ntuple entries";
         break;
      }
      if(!event) {
          LOG("gevgen_t2k", pERROR)
              << "Got a null generated neutino event! Retrying ...";
          continue;
      }
      LOG("gevgen_t2k", pINFO)
          << "Generated event: " << *event;
 
      // A valid event was generated: extract flux info (parent decay/prod
      // position/kinematics) for that simulated event so that it can be
      // passed-through.
      // Can only do so if I am generating events using the JNUBEAM flux
      // ntuples, not simple histograms
      if(!gOptUsingHistFlux) {
         flux_info = new flux::GJPARCNuFluxPassThroughInfo(
             jparc_flux_driver->PassThroughInfo());
         LOG("gevgen_t2k", pINFO)
           << "Pass-through flux info associated with generated event: "
           << *flux_info;
      }
 
      // Add event at the output ntuple, refresh the mc job monitor & clean-up
      ntpw.AddEventRecord(ievent, event);
      mcjmonitor.Update(ievent,event);
      delete event;
      if(flux_info) delete flux_info;
      ievent++;
   } //1
 
   LOG("gevgen_t2k", pNOTICE)
     << "The GENIE MC job is done generaing events - Cleaning up & exiting...";
 
   // *************************************************************************
   // * Print job statistics &
   // * calculate normalization factor for the generated sample
   // *************************************************************************
   if(!gOptUsingHistFlux && gOptUsingRootGeom)
   {
     // POT normalization will only be calculated if event generation was based
     // on beam simulation  outputs (not just histograms) & a detailed detector
     // geometry description.
     double fpot = jparc_flux_driver->POT_curravg(); // current POT in flux file
     double psc  = mcj_driver->GlobProbScale();      // interaction prob. scale
     double pot  = fpot / psc;                       // POT for generated sample
     // Get nunber of flux neutrinos read-in by flux friver, number of flux
     // neutrinos actually thrown to the event generation driver and number
     // of neutrino interactions actually generated
     long int nflx_evg = mcj_driver        -> NFluxNeutrinos();
     long int nflx     = jparc_flux_driver -> NFluxNeutrinos();
     long int nev      = ievent;
 
     LOG("gevgen_t2k", pNOTICE)
         << "\n >> Actual JNUBEAM flux file normalization:  " << fpot
             << " POT * " << ((gOptDetectorLocation == "sk") ? "cm^2" : "det")
         << "\n >> Interaction probability scaling factor:  " << psc
         << "\n >> N of flux v read-in by flux driver:      " << nflx
         << "\n >> N of flux v thrown to event gen driver:  " << nflx_evg
         << "\n >> N of generated v interactions:           " << nev
         << "\n ** Normalization for generated sample:      " << pot
             << " POT * " << ((gOptDetectorLocation == "sk") ? "cm^2" : "det");
 
     ntpw.EventTree()->SetWeight(pot); // POT
   }
 
   // *************************************************************************
   // * MC job meta-data
   // *************************************************************************
 
   genie::utils::T2KEvGenMetaData * metadata = new genie::utils::T2KEvGenMetaData;
 
   metadata -> jnubeam_file       = gOptFluxFile;
   metadata -> detector_location  = gOptDetectorLocation;
   metadata -> geom_file          = gOptRootGeom;
   metadata -> geom_top_volume    = gOptRootGeomTopVol;
   metadata -> geom_length_units  = gOptGeomLUnits;
   metadata -> geom_density_units = gOptGeomDUnits;
   metadata -> using_root_geom    = gOptUsingRootGeom;
   metadata -> target_mix         = gOptTgtMix;
   metadata -> using_hist_flux    = gOptUsingHistFlux;
   metadata -> flux_hists         = gOptFluxHst;
 
   ntpw.EventTree()->GetUserInfo()->Add(metadata);
 
   // *************************************************************************
   // * Save & clean-up
   // *************************************************************************
 
   // Save the generated event tree & close the output file
   ntpw.Save();
 
   // Clean-up
   delete geom_driver;
   delete flux_driver;
   delete mcj_driver;
   map<int,TH1D*>::iterator it = gOptFluxHst.begin();
   for( ; it != gOptFluxHst.end(); ++it) {
     TH1D * spectrum = it->second;
     if(spectrum) delete spectrum;
   }
   gOptFluxHst.clear();
 
   LOG("gevgen_t2k", pNOTICE) << "Done!";
 
   return 0;
 }

void PrintSyntax ( void )

Variable Documentation

string gOptDetectorLocation

Definition at line 507 of file gT2KEvGen.cxx.

string gOptEvFilePrefix

Definition at line 514 of file gT2KEvGen.cxx.

bool gOptExitAtEndOfFullFluxCycles

Definition at line 513 of file gT2KEvGen.cxx.

Referenced by main().

string gOptExtMaxPlXml

Definition at line 505 of file gT2KEvGen.cxx.

string gOptFluxFile

Definition at line 506 of file gT2KEvGen.cxx.

map<int,TH1D*> gOptFluxHst

Definition at line 500 of file gT2KEvGen.cxx.

int gOptFluxNCycles

Definition at line 510 of file gT2KEvGen.cxx.

Referenced by main().

double gOptFluxNorm

Definition at line 508 of file gT2KEvGen.cxx.

Referenced by main().

PDGCodeList gOptFluxNtpNuList(false)

Referenced by main().

string gOptFluxProbFileName

Definition at line 517 of file gT2KEvGen.cxx.

Referenced by main().

double gOptGeomDUnits = 0

Definition at line 504 of file gT2KEvGen.cxx.

double gOptGeomLUnits = 0

Definition at line 503 of file gT2KEvGen.cxx.

string gOptInpXSecFile

Definition at line 521 of file gT2KEvGen.cxx.

int gOptNev

Definition at line 511 of file gT2KEvGen.cxx.

double gOptPOT

Definition at line 512 of file gT2KEvGen.cxx.

bool gOptRandomFluxOffset = false

Definition at line 519 of file gT2KEvGen.cxx.

Referenced by main().

long int gOptRanSeed

Definition at line 520 of file gT2KEvGen.cxx.

string gOptRootGeom

Definition at line 501 of file gT2KEvGen.cxx.

string gOptRootGeomTopVol = ""

Definition at line 502 of file gT2KEvGen.cxx.

Long_t gOptRunNu

Definition at line 496 of file gT2KEvGen.cxx.

bool gOptSaveFluxProbsFile = false

Definition at line 516 of file gT2KEvGen.cxx.

Referenced by main().

string gOptSaveFluxProbsFileName

Definition at line 518 of file gT2KEvGen.cxx.

Referenced by main().

map<int,double> gOptTgtMix

Definition at line 499 of file gT2KEvGen.cxx.

bool gOptUseFluxProbs = false

Definition at line 515 of file gT2KEvGen.cxx.

Referenced by main().

bool gOptUsingHistFlux = false

Definition at line 498 of file gT2KEvGen.cxx.

bool gOptUsingRootGeom = false

Definition at line 497 of file gT2KEvGen.cxx.

string kDefOptEvFilePrefix = "gntp"

Definition at line 492 of file gT2KEvGen.cxx.

double kDefOptFluxNorm = 1E+21

Definition at line 491 of file gT2KEvGen.cxx.

string kDefOptGeomDUnits = "g_cm3"

Definition at line 489 of file gT2KEvGen.cxx.

string kDefOptGeomLUnits = "mm"

Definition at line 488 of file gT2KEvGen.cxx.

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 490 of file gT2KEvGen.cxx.

Functions

Variables

Function Documentation

Variable Documentation