#include <cassert>
#include <cstdlib>
#include <cctype>
#include <string>
#include <vector>
#include <sstream>
#include <map>
#include <iomanip>
#include <cmath>
#include <TRotation.h>
#include <TMath.h>
#include <TGeoShape.h>
#include <TGeoBBox.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/Ntuple/NtpMCFormat.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Utils/XSecSplineList.h"
#include "Framework/Utils/StringUtils.h"
#include "Framework/Utils/SystemUtils.h"
#include "Framework/Utils/UnitUtils.h"
#include "Framework/Utils/CmdLnArgParser.h"
#include "Framework/Utils/PrintUtils.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"

Include dependency graph for gAtmoEvGen.cxx:

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

void	PrintSyntax (void)

GAtmoFlux *	GetFlux (void)

GeomAnalyzerI *	GetGeometry (void)

int	main (int argc, char **argv)

Variables
Long_t	gOptRunNu

string	gOptFluxSim

map< int, string >	gOptFluxFiles

bool	gOptUsingRootGeom = false

map< int, double >	gOptTgtMix

string	gOptRootGeom

string	gOptRootGeomTopVol = ""

double	gOptGeomLUnits = 0

double	gOptGeomDUnits = 0

string	gOptExtMaxPlXml

int	gOptNev = -1

double	gOptKtonYrExposure = -1

double	gOptSecExposure = -1

double	gOptEvMin

double	gOptEvMax

string	gOptEvFilePrefix

TRotation	gOptRot

long int	gOptRanSeed

string	gOptInpXSecFile

double	gOptRL = -1

double	gOptRT = -1

NtpMCFormat_t	kDefOptNtpFormat = kNFGHEP

string	kDefOptEvFilePrefix = "gntp"

string	kDefOptGeomLUnits = "mm"

string	kDefOptGeomDUnits = "g_cm3"

double	kDefOptEvMin = 0.5

double	kDefOptEvMax = 50.0

Function Documentation

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

Definition at line 563 of file gAtmoEvGen.cxx.

References genie::CmdLnArgParser::ArgAsDouble(), genie::CmdLnArgParser::ArgAsInt(), genie::CmdLnArgParser::ArgAsLong(), genie::CmdLnArgParser::ArgAsString(), genie::utils::system::FileExists(), genie::PDGLibrary::Find(), genie::gAbortingInErr, geom, gOptEvFilePrefix, gOptEvMax, gOptEvMin, gOptExtMaxPlXml, gOptFluxFiles, gOptFluxSim, gOptGeomDUnits, gOptGeomLUnits, gOptInpXSecFile, gOptKtonYrExposure, gOptNev, gOptRanSeed, gOptRL, gOptRootGeom, gOptRootGeomTopVol, gOptRot, gOptRT, gOptRunNu, gOptSecExposure, gOptTgtMix, gOptUsingRootGeom, genie::RunOpt::Instance(), genie::PDGLibrary::Instance(), kDefOptEvFilePrefix, kDefOptEvMax, kDefOptGeomDUnits, kDefOptGeomLUnits, LOG, lunits, genie::CmdLnArgParser::OptionExists(), pDEBUG, pFATAL, pINFO, pNOTICE, genie::utils::print::PrintFramedMesg(), PrintSyntax(), genie::RunOpt::ReadFromCommandLine(), genie::utils::str::Split(), and genie::utils::units::UnitFromString().

Referenced by main().

 {
 // Get the command line arguments
 
   RunOpt::Instance()->ReadFromCommandLine(argc,argv);
 
   LOG("gevgen_atmo", pNOTICE) << "Parsing command line arguments";
 
   CmdLnArgParser parser(argc,argv);
 
   // help?
   bool help = parser.OptionExists('h');
   if(help) {
       PrintSyntax();
       exit(0);
   }
 
   //
   // *** run number
   //
   if( parser.OptionExists('r') ) {
     LOG("gevgen_atmo", pDEBUG) << "Reading MC run number";
     gOptRunNu = parser.ArgAsLong('r');
   } else {
     LOG("gevgen_atmo", pDEBUG) << "Unspecified run number - Using default";
     gOptRunNu = 100000000;
   } //-r
 
   //
   // *** exposure
   //
 
   // in number of events
   bool have_required_statistics = false;
   if( parser.OptionExists('n') ) {
     LOG("gevgen_atmo", pDEBUG)
         << "Reading number of events to generate";
     gOptNev = parser.ArgAsInt('n');
     have_required_statistics = true;
   }//-n?
   // or, in kton*yrs
   if( parser.OptionExists('e') ) {
     if(have_required_statistics) {
       LOG("gevgen_atmo", pFATAL)
          << "Can't request exposure both in terms of number of events and  kton*yrs"
          << "\nUse just one of the -n and -e options";
       PrintSyntax();
       gAbortingInErr = true;
       exit(1);
     }
     LOG("gevgen_atmo", pDEBUG)
         << "Reading requested exposure in kton*yrs";
     gOptKtonYrExposure = parser.ArgAsDouble('e');
     have_required_statistics = true;
   }//-e?
 
   if (parser.OptionExists('T')) {
     if (have_required_statistics) {
       LOG("gevgen_atmo", pFATAL)
          << "Can't request exposure both in terms of number of events or kton*yrs and time"
          << "\nUse just one of the -n, -e, or -T options";
       PrintSyntax();
       gAbortingInErr = true;
       exit(1);
     }
     LOG("gevgen_atmo", pDEBUG)
         << "Reading requested exposure in seconds";
     gOptSecExposure = parser.ArgAsDouble('T');
     have_required_statistics = true;
   }
 
   if(!have_required_statistics) {
     LOG("gevgen_atmo", pFATAL)
        << "You must request exposure either in terms of number of events and  kton*yrs"
        << "\nUse any of the -n, -e options";
     PrintSyntax();
     gAbortingInErr = true;
     exit(1);
   }
 
   //
   // *** event file prefix
   //
   if( parser.OptionExists('o') ) {
     LOG("gevgen_atmo", pDEBUG) << "Reading the event filename prefix";
     gOptEvFilePrefix = parser.ArgAsString('o');
   } else {
     LOG("gevgen_atmo", pDEBUG)
       << "Will set the default event filename prefix";
     gOptEvFilePrefix = kDefOptEvFilePrefix;
   } //-o
 
   //
   // *** neutrino energy range
   //
   if( parser.OptionExists('E') ) {
     LOG("gevgen_atmo", pINFO) << "Reading neutrino energy range";
     string nue = parser.ArgAsString('E');
 
     // must be a comma separated set of values
     if(nue.find(",") != string::npos) {
        // split the comma separated list
        vector<string> nurange = utils::str::Split(nue, ",");
        assert(nurange.size() == 2);
        double emin = atof(nurange[0].c_str());
        double emax = atof(nurange[1].c_str());
        assert(emax>emin && emin>=0.);
        gOptEvMin = emin;
        gOptEvMax = emax;
     } else {
       LOG("gevgen_atmo", pFATAL)
         << "Invalid energy range. Use `-E emin,emax', eg `-E 0.5,100.";
       PrintSyntax();
       gAbortingInErr = true;
       exit(1);
     }
   } else {
      LOG("gevgen_atmo", pNOTICE)
         << "No -e option. Using default energy range";
      gOptEvMin = kDefOptEvMax;
      gOptEvMax = kDefOptEvMax;
   }
 
   //
   // *** flux files
   //
   // syntax:
   // simulation:/path/file.data[neutrino_code],/path/file.data[neutrino_code],...
   //
   if( parser.OptionExists('f') ) {
     LOG("gevgen_atmo", pDEBUG) << "Getting input flux files";
     string flux = parser.ArgAsString('f');
 
     // get flux simulation info (FLUKA,BGLRS) so as to instantiate the
     // appropriate flux driver
     string::size_type jsimend = flux.find_first_of(":",0);
     if(jsimend==string::npos) {
        LOG("gevgen_atmo", pFATAL)
            << "You need to specify the flux file source";
        PrintSyntax();
        gAbortingInErr = true;
        exit(1);
     }
     gOptFluxSim = flux.substr(0,jsimend);
     for(string::size_type i=0; i<gOptFluxSim.size(); i++) {
        gOptFluxSim[i] = toupper(gOptFluxSim[i]);
     }
     if((gOptFluxSim != "FLUKA") &&
        (gOptFluxSim != "BGLRS") &&
        (gOptFluxSim != "HAKKM")) {
         LOG("gevgen_atmo", pFATAL)
              << "The flux file source needs to be one of <FLUKA,BGLRS,HAKKM>";
         PrintSyntax();
         gAbortingInErr = true;
         exit(1);
     }
     // now get the list of input files and the corresponding neutrino codes.
     flux.erase(0,jsimend+1);
     vector<string> fluxv = utils::str::Split(flux,",");
     vector<string>::const_iterator fluxiter = fluxv.begin();
     for( ; fluxiter != fluxv.end(); ++fluxiter) {
        string filename_and_pdg = *fluxiter;
        string::size_type open_bracket  = filename_and_pdg.find("[");
        string::size_type close_bracket = filename_and_pdg.find("]");
        if (open_bracket ==string::npos ||
            close_bracket==string::npos)
        {
            LOG("gevgen_atmo", pFATAL)
               << "You made an error in specifying the flux info";
            PrintSyntax();
            gAbortingInErr = true;
            exit(1);
        }
        string::size_type ibeg = 0;
        string::size_type iend = open_bracket;
        string::size_type jbeg = open_bracket+1;
        string::size_type jend = close_bracket;
        string flux_filename   = filename_and_pdg.substr(ibeg,iend-ibeg);
        string neutrino_pdg    = filename_and_pdg.substr(jbeg,jend-jbeg);
        gOptFluxFiles.insert(
           map<int,string>::value_type(atoi(neutrino_pdg.c_str()), flux_filename));
     }
     if(gOptFluxFiles.size() == 0) {
        LOG("gevgen_atmo", pFATAL)
           << "You must specify at least one flux file!";
        PrintSyntax();
        gAbortingInErr = true;
        exit(1);
     }
 
   } else {
     LOG("gevgen_atmo", pFATAL) << "No flux info was specified! Use the -f option.";
     PrintSyntax();
     gAbortingInErr = true;
     exit(1);
   }
 
   // *** options to fine tune the flux ray generation surface
 
   if( parser.OptionExists("flux-ray-generation-surface-distance") ) {
     LOG("gevgen_atmo", pINFO)
       << "Reading distance of flux ray generation surface";
     gOptRL = parser.ArgAsDouble("flux-ray-generation-surface-distance");
   } else {
     LOG("gevgen_atmo", pINFO)
       << "Unspecified distance of flux ray generation surface - Using default";
   }
 
   if( parser.OptionExists("flux-ray-generation-surface-radius") ) {
     LOG("gevgen_atmo", pINFO)
       << "Reading radius of flux ray generation surface";
     gOptRT = parser.ArgAsDouble("flux-ray-generation-surface-radius");
   } else {
     LOG("gevgen_atmo", pINFO)
       << "Unspecified radius of flux ray generation surface - Using default";
   }
 
   //
   // *** geometry
   //
   string geom = "";
   string lunits, dunits;
   if( parser.OptionExists('g') ) {
     LOG("gevgen_atmo", pDEBUG) << "Getting input geometry";
     geom = parser.ArgAsString('g');
 
     // is it a ROOT file that contains a ROOT geometry?
     bool accessible_geom_file =
         utils::system::FileExists(geom.c_str());
     if (accessible_geom_file) {
       gOptRootGeom      = geom;
       gOptUsingRootGeom = true;
     }
   } else {
       LOG("gevgen_atmo", pFATAL)
         << "No geometry option specified - Exiting";
       PrintSyntax();
       gAbortingInErr = true;
       exit(1);
   } //-g
 
   if(gOptUsingRootGeom) {
      // using a ROOT geometry - get requested geometry units
 
      // legth units:
      if( parser.OptionExists('L') ) {
         LOG("gevgen_atmo", pDEBUG)
            << "Checking for input geometry length units";
         lunits = parser.ArgAsString('L');
      } else {
         LOG("gevgen_atmo", pDEBUG) << "Using default geometry length units";
         lunits = kDefOptGeomLUnits;
      } // -L
      // density units:
      if( parser.OptionExists('D') ) {
         LOG("gevgen_atmo", pDEBUG)
            << "Checking for input geometry density units";
         dunits = parser.ArgAsString('D');
      } else {
         LOG("gevgen_atmo", pDEBUG) << "Using default geometry density units";
         dunits = kDefOptGeomDUnits;
      } // -D
      gOptGeomLUnits = genie::utils::units::UnitFromString(lunits);
      gOptGeomDUnits = genie::utils::units::UnitFromString(dunits);
 
      // check whether an event generation volume name has been
      // specified -- default is the 'top volume'
      if( parser.OptionExists('t') ) {
         LOG("gevgen_atmo", pDEBUG) << "Checking for input volume name";
         gOptRootGeomTopVol = parser.ArgAsString('t');
      } else {
         LOG("gevgen_atmo", pDEBUG) << "Using the <master volume>";
      } // -t
 
      // check whether an XML file with the maximum (density weighted)
      // path lengths for each detector material is specified -
      // otherwise will compute the max path lengths at job init
      if( parser.OptionExists('m') ) {
         LOG("gevgen_atmo", pDEBUG)
               << "Checking for maximum path lengths XML file";
         gOptExtMaxPlXml = parser.ArgAsString('m');
      } else {
         LOG("gevgen_atmo", pDEBUG)
            << "Will compute the maximum path lengths at job init";
         gOptExtMaxPlXml = "";
      } // -m
   } // using root geom?
 
   else {
     // User has specified a target mix.
     // Decode the list of target pdf codes & their corresponding weight fraction
     // (specified as 'pdg_code_1[fraction_1],pdg_code_2[fraction_2],...')
     // See documentation on top section of this file.
     //
     gOptTgtMix.clear();
     vector<string> tgtmix = utils::str::Split(geom,",");
     if(tgtmix.size()==1) {
          int    pdg = atoi(tgtmix[0].c_str());
          double wgt = 1.0;
          gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
     } else {
       vector<string>::const_iterator tgtmix_iter = tgtmix.begin();
       for( ; tgtmix_iter != tgtmix.end(); ++tgtmix_iter) {
          string tgt_with_wgt = *tgtmix_iter;
          string::size_type open_bracket  = tgt_with_wgt.find("[");
          string::size_type close_bracket = tgt_with_wgt.find("]");
          if (open_bracket ==string::npos ||
              close_bracket==string::npos)
          {
              LOG("gevgen_atmo", pFATAL)
                 << "You made an error in specifying the target mix";
              PrintSyntax();
              gAbortingInErr = true;
              exit(1);
          }
          string::size_type ibeg = 0;
          string::size_type iend = open_bracket;
          string::size_type jbeg = open_bracket+1;
          string::size_type jend = close_bracket;
          int    pdg = atoi(tgt_with_wgt.substr(ibeg,iend-ibeg).c_str());
          double wgt = atof(tgt_with_wgt.substr(jbeg,jend-jbeg).c_str());
          LOG("gevgen_atmo", pDEBUG)
             << "Adding to target mix: pdg = " << pdg << ", wgt = " << wgt;
          gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
 
       }// tgtmix_iter
     } // >1 materials in mix
   } // using tgt mix?
 
   //
   // Coordinate rotation matrix
   //
   gOptRot.SetToIdentity();
   if( parser.OptionExists('R') ) {
     string rotarg = parser.ArgAsString('R');
     //get convention
     string::size_type j = rotarg.find_first_of(":",0);
     string convention = "";
     if(j==string::npos) { convention = "X"; }
     else                { convention = rotarg.substr(0,j); }
     //get angles phi,theta,psi
     rotarg.erase(0,j+1);
     vector<string> euler_angles = utils::str::Split(rotarg,",");
     if(euler_angles.size() != 3) {
        LOG("gevgen_atmo", pFATAL)
          << "You didn't specify all 3 Euler angles using the -R option";
        PrintSyntax();
        gAbortingInErr = true;
        exit(1);
     }
     double phi   = atof(euler_angles[0].c_str());
     double theta = atof(euler_angles[1].c_str());
     double psi   = atof(euler_angles[2].c_str());
     //set Euler angles using appropriate convention
     if(convention.find("X")!=string::npos ||
        convention.find("x")!=string::npos)
     {
        LOG("gevgen_atmo", pNOTICE) << "Using X-convention for input Euler angles";
        gOptRot.SetXEulerAngles(phi,theta,psi);
     } else
     if(convention.find("Y")!=string::npos ||
        convention.find("y")!=string::npos)
     {
        LOG("gevgen_atmo", pNOTICE) << "Using Y-convention for input Euler angles";
        gOptRot.SetYEulerAngles(phi,theta,psi);
     } else {
        LOG("gevgen_atmo", pFATAL)
          << "Unknown Euler angle convention. Please use the X- or Y-convention";
        PrintSyntax();
        gAbortingInErr = true;
        exit(1);
     }
     //invert?
     if(convention.find("^-1")!=string::npos) {
        LOG("gevgen_atmo", pNOTICE) << "Inverting rotation matrix";
        gOptRot.Invert();
     }
   }
 
   //
   // *** random number seed
   //
   if( parser.OptionExists("seed") ) {
     LOG("gevgen_atmo", pINFO) << "Reading random number seed";
     gOptRanSeed = parser.ArgAsLong("seed");
   } else {
     LOG("gevgen_atmo", pINFO) << "Unspecified random number seed - Using default";
     gOptRanSeed = -1;
   }
 
   //
   // *** input cross-section file
   //
   if( parser.OptionExists("cross-sections") ) {
     LOG("gevgen_atmo", pINFO) << "Reading cross-section file";
     gOptInpXSecFile = parser.ArgAsString("cross-sections");
   } else {
     LOG("gevgen_atmo", pINFO) << "Unspecified cross-section file";
     gOptInpXSecFile = "";
   }
 
   //
   // print-out summary
   //
 
   PDGLibrary * pdglib = PDGLibrary::Instance();
 
   ostringstream gminfo;
   if (gOptUsingRootGeom) {
     gminfo << "Using ROOT geometry - file: " << gOptRootGeom
            << ", top volume: "
            << ((gOptRootGeomTopVol.size()==0) ? "<master volume>" : gOptRootGeomTopVol)
            << ", max{PL} file: "
            << ((gOptExtMaxPlXml.size()==0) ? "<none>" : gOptExtMaxPlXml)
            << ", length  units: " << lunits
            << ", density units: " << dunits;
   } else {
     gminfo << "Using target mix - ";
     map<int,double>::const_iterator iter;
     for(iter = gOptTgtMix.begin(); iter != gOptTgtMix.end(); ++iter) {
           int    pdg_code = iter->first;
           double wgt      = iter->second;
           TParticlePDG * p = pdglib->Find(pdg_code);
           if(p) {
             string name = p->GetName();
             gminfo << "(" << name << ") -> " << 100*wgt << "% / ";
           }//p?
     }
   }
 
   ostringstream fluxinfo;
   fluxinfo << "Using " << gOptFluxSim << " flux files: ";
   map<int,string>::const_iterator file_iter = gOptFluxFiles.begin();
   for( ; file_iter != gOptFluxFiles.end(); ++file_iter) {
      int neutrino_code = file_iter->first;
      string filename   = file_iter->second;
      TParticlePDG * p = pdglib->Find(neutrino_code);
      if(p) {
         string name = p->GetName();
         fluxinfo << "(" << name << ") -> " << filename << " / ";
      }
   }
   fluxinfo << "Flux ray generation surface - Distance = "
            << gOptRL << " m, Radius = " << gOptRT << " m";
 
   ostringstream expinfo;
   if(gOptNev > 0)            { expinfo << gOptNev            << " events";   }
   if(gOptKtonYrExposure > 0) { expinfo << gOptKtonYrExposure << " kton*yrs"; }
 
   ostringstream rotation;
   rotation << "\t| " <<  gOptRot.XX() << "  " << gOptRot.XY() << "  " << gOptRot.XZ() << " |\n";
   rotation << "\t| " <<  gOptRot.YX() << "  " << gOptRot.YY() << "  " << gOptRot.YZ() << " |\n";
   rotation << "\t| " <<  gOptRot.ZX() << "  " << gOptRot.ZY() << "  " << gOptRot.ZZ() << " |\n";
 
   LOG("gevgen_atmo", pNOTICE)
      << "\n\n"
      << utils::print::PrintFramedMesg("gevgen_atmo job configuration");
 
   LOG("gevgen_atmo", pNOTICE)
      << "\n"
      << "\n @@ Run number: " << gOptRunNu
      << "\n @@ Random number seed: " << gOptRanSeed
      << "\n @@ Using cross-section file: " << gOptInpXSecFile
      << "\n @@ Geometry"
      << "\n\t" << gminfo.str()
      << "\n @@ Flux"
      << "\n\t" << fluxinfo.str()
      << "\n @@ Exposure"
      << "\n\t" << expinfo.str()
      << "\n @@ Cuts"
      << "\n\t Using energy range = (" << gOptEvMin << " GeV, " << gOptEvMax << " GeV)"
      << "\n @@ Coordinate transformation (Rotation THZ -> User-defined coordinate system)"
      << "\n" << rotation.str()
      << "\n\n";
 
   //
   // final checks
   //
   if(gOptKtonYrExposure > 0) {
     LOG("gevgen_atmo", pFATAL)
       << "\n Option to set exposure in terms of kton*yrs not supported just yet!"
       << "\n Try the -n option instead";
     PrintSyntax();
     gAbortingInErr = true;
     exit(1);
   }
 }

GFluxI * GetFlux ( void )

Definition at line 505 of file gAtmoEvGen.cxx.

References genie::flux::GAtmoFlux::AddFluxFile(), genie::flux::GAtmoFlux::ForceMaxEnergy(), genie::flux::GAtmoFlux::ForceMinEnergy(), genie::gAbortingInErr, genie::units::GeV, gOptEvMax, gOptEvMin, gOptFluxFiles, gOptFluxSim, gOptRL, gOptRot, gOptRT, genie::flux::GAtmoFlux::LoadFluxData(), LOG, pFATAL, genie::flux::GAtmoFlux::SetRadii(), and genie::flux::GAtmoFlux::SetUserCoordSystem().

Referenced by genie::flux::GAtmoFlux::GenerateNext_1try(), and main().

 {
 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
 
   // Instantiate appropriate concrete flux driver
   GAtmoFlux * atmo_flux_driver = 0;
   if(gOptFluxSim == "FLUKA") {
      GFLUKAAtmoFlux * fluka_flux = new GFLUKAAtmoFlux;
      atmo_flux_driver = dynamic_cast<GAtmoFlux *>(fluka_flux);
   } else
   if(gOptFluxSim == "BGLRS") {
      GBGLRSAtmoFlux * bartol_flux = new GBGLRSAtmoFlux;
      atmo_flux_driver = dynamic_cast<GAtmoFlux *>(bartol_flux);
   } else
   if(gOptFluxSim == "HAKKM") {
      GHAKKMAtmoFlux * honda_flux = new GHAKKMAtmoFlux;
      atmo_flux_driver = dynamic_cast<GAtmoFlux *>(honda_flux);
   } else {
      LOG("gevgen_atmo", pFATAL) << "Unknown flux simulation: " << gOptFluxSim;
      gAbortingInErr = true;
      exit(1);
   }
   // Configure GAtmoFlux options (common to all concrete atmospheric flux drivers)
   // set min/max energy:
   atmo_flux_driver->ForceMinEnergy (gOptEvMin * units::GeV);
   atmo_flux_driver->ForceMaxEnergy (gOptEvMax * units::GeV);
   // set flux files:
   map<int,string>::const_iterator file_iter = gOptFluxFiles.begin();
   for( ; file_iter != gOptFluxFiles.end(); ++file_iter) {
     int neutrino_code = file_iter->first;
     string filename   = file_iter->second;
     atmo_flux_driver->AddFluxFile(neutrino_code, filename);
   }
 
   if (!atmo_flux_driver->LoadFluxData()) {
     LOG("gevgen_atmo", pFATAL) << "Error loading flux data. Quitting...";
     gAbortingInErr = true;
     exit(1);
   }
 
   // configure flux generation surface:
   atmo_flux_driver->SetRadii(gOptRL, gOptRT);
   // set rotation for coordinate tranformation from the topocentric horizontal
   // system to a user-defined coordinate system:
   if(!gOptRot.IsIdentity()) {
      atmo_flux_driver->SetUserCoordSystem(gOptRot);
   }
 
 #else
   LOG("gevgen_atmo", pFATAL) << "You need to enable the GENIE flux drivers first!";
   LOG("gevgen_atmo", pFATAL) << "Use --enable-flux-drivers at the configuration step.";
   gAbortingInErr = true;
   exit(1);
 #endif
 
   return atmo_flux_driver;
 }

GeomAnalyzerI * GetGeometry ( void )

Definition at line 433 of file gAtmoEvGen.cxx.

References genie::gAbortingInErr, genie::geometry::ROOTGeomAnalyzer::GetGeometry(), gOptGeomDUnits, gOptGeomLUnits, gOptRL, gOptRootGeom, gOptRootGeomTopVol, gOptRT, gOptTgtMix, gOptUsingRootGeom, genie::geometry::ROOTGeomAnalyzer::LengthUnits(), LOG, pFATAL, pNOTICE, and genie::utils::geometry::RecursiveExhaust().

Referenced by main().

 {
   GeomAnalyzerI * geom_driver = 0;
 
 #ifdef __GENIE_GEOM_DRIVERS_ENABLED__
 
   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);
     rgeom -> SetTopVolName   (gOptRootGeomTopVol);
     // getting the bounding box dimensions along z so as to set the
     // appropriate upstream generation surface for the JPARC flux driver
     TGeoVolume * topvol = rgeom->GetGeometry()->GetTopVolume();
     if(!topvol) {
       LOG("gevgen_atmo", pFATAL) << " ** Null top ROOT geometry volume!";
       gAbortingInErr = true;
       exit(1);
     }
 
     /* If flux generation surface isn't defined, get the bounding box for the
      * geometry and set something appropriate. */
     TGeoShape *bounding_box = topvol->GetShape();
     TGeoBBox *box = (TGeoBBox *) bounding_box;
     double dx = box->GetDX()*rgeom->LengthUnits();
     double dy = box->GetDY()*rgeom->LengthUnits();
     double dz = box->GetDZ()*rgeom->LengthUnits();
 
     if (gOptRL < 0 && gOptRT < 0) {
       gOptRL = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
       gOptRT = gOptRL;
       LOG("gevgen_atmo", pNOTICE) << "Setting flux longitudinal and transverse radius to " << setprecision(2) << gOptRL << " meters based on bounding box of ROOT geometry.";
     }
 
     // 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);
   }
 
 #else
   LOG("gevgen_atmo", pFATAL) << "You need to enable the GENIE geometry drivers first!";
   LOG("gevgen_atmo", pFATAL) << "Use --enable-geom-drivers at the configuration step.";
   gAbortingInErr = true;
   exit(1);
 #endif
 
   return geom_driver;
 }

int main	(	int	argc,
		char **	argv
	)

Definition at line 327 of file gAtmoEvGen.cxx.

Referenced by genie::AlgConfigPool::LoadMasterConfigs().

 {
   GAtmoFlux* atmo_flux_driver;
   double total_flux, expected_neutrinos;
 
   // Parse command line arguments
   GetCommandLineArgs(argc,argv);
 
   if ( ! RunOpt::Instance()->Tune() ) {
     LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
     exit(-1);
   }
   RunOpt::Instance()->BuildTune();
 
   // Iinitialization of random number generators, cross-section table, messenger, cache etc...
   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);
 
   // get geometry driver
   GeomAnalyzerI * geom_driver = GetGeometry();
 
   if (gOptRT < 0) {
     gOptRT = 1000; // m
     LOG("gevgen_atmo", pWARN) << "Warning! Flux surface transverse radius not specified so using default value of " << gOptRT << " meters!";
   }
 
   if (gOptRL < 0) {
     gOptRL = 1000; // m
     LOG("gevgen_atmo", pWARN) << "Warning! Flux surface longitudinal radius not specified so using default value of " << gOptRL << " meters!";
   }
 
   // get flux driver
   atmo_flux_driver = GetFlux();
 
   // Cast to GFluxI, the generic flux driver interface
   GFluxI *flux_driver = dynamic_cast<GFluxI *>(atmo_flux_driver);
 
   // create the GENIE monte carlo job driver
   GMCJDriver* mcj_driver = new GMCJDriver;
   mcj_driver->SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
   mcj_driver->UseFluxDriver(flux_driver);
   mcj_driver->UseGeomAnalyzer(geom_driver);
   mcj_driver->Configure();
   mcj_driver->UseSplines();
   /* Note: For the method of calculating the total number of events using a
    * livetime we *need* to force a single probability scale. So if you change
    * the next line you need to make sure that the user didn't specify the -T
    * option. */
   mcj_driver->ForceSingleProbScale();
 
   // initialize an ntuple writer
   NtpWriter ntpw(kDefOptNtpFormat, gOptRunNu, gOptRanSeed);
   ntpw.CustomizeFilenamePrefix(gOptEvFilePrefix);
   ntpw.Initialize();
 
   // Create a MC job monitor for a periodically updated status file
   GMCJMonitor mcjmonitor(gOptRunNu);
   mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());
 
   // Set GHEP print level
   GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
 
   total_flux = atmo_flux_driver->GetTotalFluxInEnergyRange();
   LOG("gevgen_atmo", pNOTICE) << "Total atmospheric neutrino flux is " << setprecision(2) << total_flux << " neutrinos per m^2 per second.";
   if (gOptSecExposure > 0) {
     /* Calculate the expected value of the total number of neutrinos we need to
      * throw. We do this by multiplying the total flux by the exposure time in
      * seconds and the area of the flux surface. */
     expected_neutrinos = total_flux*gOptSecExposure*atmo_flux_driver->GetFluxSurfaceArea();
     LOG("gevgen_atmo", pNOTICE) << "Simulating an exposure of " << setprecision(0) << gOptSecExposure << " seconds which corresponds to a total of " << setprecision(0) << expected_neutrinos << " neutrinos.";
   }
 
   // event loop
   for(int iev = 0; gOptNev > 0 ? iev < gOptNev : 1; iev++) {
 
     // generate next event
     EventRecord* event = mcj_driver->GenerateEvent();
 
     // print-out
     LOG("gevgen_atmo", pNOTICE) << "Generated event: " << *event;
 
     // save the event, refresh the mc job monitor
     ntpw.AddEventRecord(iev, event);
     mcjmonitor.Update(iev,event);
 
     // clean-up
     delete event;
 
     if (gOptSecExposure > 0 && mcj_driver->NFluxNeutrinos()/mcj_driver->GlobProbScale() > expected_neutrinos) {
       break;
     }
   }
 
   // save the event file
   ntpw.Save();
 
   // clean-up
   delete geom_driver;
   delete atmo_flux_driver;
   delete mcj_driver;
 
   return 0;
 }

void PrintSyntax ( void )

Definition at line 1051 of file gAtmoEvGen.cxx.

References LOG, and pFATAL.

Referenced by DecodeCommandLine(), GetCommandLineArgs(), GetEventRange(), main(), ParseFluxFileConfig(), and ParseFluxHst().

 {
   LOG("gevgen_atmo", pFATAL)
    << "\n **Syntax**"
    << "\n gevgen_atmo [-h]"
    << "\n           [-r run#]"
    << "\n            -f simulation:flux_file[neutrino_code],..."
    << "\n            -g geometry"
    << "\n           [-R coordinate_rotation_matrix]"
    << "\n           [-t geometry_top_volume_name]"
    << "\n           [-m max_path_lengths_xml_file]"
    << "\n           [-L geometry_length_units]"
    << "\n           [-D geometry_density_units]"
    << "\n           <-n n_of_events,"
    << "\n            -e exposure_in_kton_x_yrs"
    << "\n            -T exposure_in_seconds>"
    << "\n            -E min_energy,max_energy"
    << "\n           [-o output_event_file_prefix]"
    << "\n           [--flux-ray-generation-surface-distance]"
    << "\n           [--flux-ray-generation-surface-radius]"
    << "\n           [--seed random_number_seed]"
    << "\n            --cross-sections xml_file"
    << "\n           [--event-generator-list list_name]"
    << "\n           [--message-thresholds xml_file]"
    << "\n           [--unphysical-event-mask mask]"
    << "\n           [--event-record-print-level level]"
    << "\n           [--mc-job-status-refresh-rate  rate]"
    << "\n           [--cache-file root_file]"
    << "\n"
    << " Please also read the detailed documentation at http://www.genie-mc.org"
    << "\n";
 }