gAtmoEvGen.cxx File Reference

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

Go to the source code of this file.

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.

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

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

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

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

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

Variable Documentation

string gOptEvFilePrefix

Definition at line 310 of file gAtmoEvGen.cxx.

double gOptEvMax

Definition at line 309 of file gAtmoEvGen.cxx.

double gOptEvMin

Definition at line 308 of file gAtmoEvGen.cxx.

string gOptExtMaxPlXml

Definition at line 304 of file gAtmoEvGen.cxx.

map<int,string> gOptFluxFiles

Definition at line 297 of file gAtmoEvGen.cxx.

string gOptFluxSim

Definition at line 296 of file gAtmoEvGen.cxx.

double gOptGeomDUnits = 0

Definition at line 303 of file gAtmoEvGen.cxx.

double gOptGeomLUnits = 0

Definition at line 302 of file gAtmoEvGen.cxx.

string gOptInpXSecFile

Definition at line 313 of file gAtmoEvGen.cxx.

double gOptKtonYrExposure = -1

Definition at line 306 of file gAtmoEvGen.cxx.

int gOptNev = -1

Definition at line 305 of file gAtmoEvGen.cxx.

long int gOptRanSeed

Definition at line 312 of file gAtmoEvGen.cxx.

double gOptRL = -1

Definition at line 314 of file gAtmoEvGen.cxx.

string gOptRootGeom

Definition at line 300 of file gAtmoEvGen.cxx.

string gOptRootGeomTopVol = ""

Definition at line 301 of file gAtmoEvGen.cxx.

TRotation gOptRot

Definition at line 311 of file gAtmoEvGen.cxx.

double gOptRT = -1

Definition at line 315 of file gAtmoEvGen.cxx.

Long_t gOptRunNu

Definition at line 295 of file gAtmoEvGen.cxx.

double gOptSecExposure = -1

Definition at line 307 of file gAtmoEvGen.cxx.

map<int,double> gOptTgtMix

Definition at line 299 of file gAtmoEvGen.cxx.

bool gOptUsingRootGeom = false

Definition at line 298 of file gAtmoEvGen.cxx.

string kDefOptEvFilePrefix = "gntp"

Definition at line 320 of file gAtmoEvGen.cxx.

double kDefOptEvMax = 50.0

Definition at line 324 of file gAtmoEvGen.cxx.

double kDefOptEvMin = 0.5

Definition at line 323 of file gAtmoEvGen.cxx.

string kDefOptGeomDUnits = "g_cm3"

Definition at line 322 of file gAtmoEvGen.cxx.

string kDefOptGeomLUnits = "mm"

Definition at line 321 of file gAtmoEvGen.cxx.

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 319 of file gAtmoEvGen.cxx.