Functions | Variables
gAtmoEvGen.cxx File Reference
#include <cassert>
#include <cstdlib>
#include <cctype>
#include <string>
#include <vector>
#include <sstream>
#include <map>
#include <TRotation.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)
 
GFluxIGetFlux (void)
 
GeomAnalyzerIGetGeometry (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 gOptEvMin
 
double gOptEvMax
 
string gOptEvFilePrefix
 
TRotation gOptRot
 
long int gOptRanSeed
 
string gOptInpXSecFile
 
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 497 of file gAtmoEvGen.cxx.

References genie::CmdLnArgParser::ArgAsDouble(), genie::CmdLnArgParser::ArgAsInt(), genie::CmdLnArgParser::ArgAsLong(), genie::CmdLnArgParser::ArgAsString(), ana::assert(), emax, emin, exit(), genie::utils::system::FileExists(), shutoffs::filename, genie::PDGLibrary::Find(), flux, genie::gAbortingInErr, geom(), gOptEvFilePrefix, gOptEvMax, gOptEvMin, gOptExtMaxPlXml, gOptFluxFiles, gOptFluxSim, gOptGeomDUnits, gOptGeomLUnits, gOptInpXSecFile, gOptKtonYrExposure, gOptNev, gOptRanSeed, gOptRootGeom, gOptRootGeomTopVol, gOptRot, gOptRunNu, gOptTgtMix, gOptUsingRootGeom, plot_validation_datamc::help, MECModelEnuComparisons::i, genie::PDGLibrary::Instance(), genie::RunOpt::Instance(), calib::j, kDefOptEvFilePrefix, kDefOptEvMax, kDefOptGeomDUnits, kDefOptGeomLUnits, LOG, genie::CmdLnArgParser::OptionExists(), plot_validation_datamc::parser, pDEBUG, make_root_from_grid_output::pdg, pFATAL, pINFO, pNOTICE, genie::utils::print::PrintFramedMesg(), PrintSyntax(), genie::RunOpt::ReadFromCommandLine(), genie::utils::str::Split(), chisquared::theta, genie::utils::units::UnitFromString(), and wgt.

Referenced by main().

498 {
499 // Get the command line arguments
500 
501  RunOpt::Instance()->ReadFromCommandLine(argc,argv);
502 
503  LOG("gevgen_atmo", pNOTICE) << "Parsing command line arguments";
504 
505  CmdLnArgParser parser(argc,argv);
506 
507  // help?
508  bool help = parser.OptionExists('h');
509  if(help) {
510  PrintSyntax();
511  exit(0);
512  }
513 
514  //
515  // *** run number
516  //
517  if( parser.OptionExists('r') ) {
518  LOG("gevgen_atmo", pDEBUG) << "Reading MC run number";
519  gOptRunNu = parser.ArgAsLong('r');
520  } else {
521  LOG("gevgen_atmo", pDEBUG) << "Unspecified run number - Using default";
522  gOptRunNu = 100000000;
523  } //-r
524 
525  //
526  // *** exposure
527  //
528 
529  // in number of events
530  bool have_required_statistics = false;
531  if( parser.OptionExists('n') ) {
532  LOG("gevgen_atmo", pDEBUG)
533  << "Reading number of events to generate";
534  gOptNev = parser.ArgAsInt('n');
535  have_required_statistics = true;
536  }//-n?
537  // or, in kton*yrs
538  if( parser.OptionExists('e') ) {
539  if(have_required_statistics) {
540  LOG("gevgen_atmo", pFATAL)
541  << "Can't request exposure both in terms of number of events and kton*yrs"
542  << "\nUse just one of the -n and -e options";
543  PrintSyntax();
544  gAbortingInErr = true;
545  exit(1);
546  }
547  LOG("gevgen_atmo", pDEBUG)
548  << "Reading requested exposure in kton*yrs";
549  gOptKtonYrExposure = parser.ArgAsDouble('e');
550  have_required_statistics = true;
551  }//-e?
552  if(!have_required_statistics) {
553  LOG("gevgen_atmo", pFATAL)
554  << "You must request exposure either in terms of number of events and kton*yrs"
555  << "\nUse any of the -n, -e options";
556  PrintSyntax();
557  gAbortingInErr = true;
558  exit(1);
559  }
560 
561  //
562  // *** event file prefix
563  //
564  if( parser.OptionExists('o') ) {
565  LOG("gevgen_atmo", pDEBUG) << "Reading the event filename prefix";
566  gOptEvFilePrefix = parser.ArgAsString('o');
567  } else {
568  LOG("gevgen_atmo", pDEBUG)
569  << "Will set the default event filename prefix";
571  } //-o
572 
573  //
574  // *** neutrino energy range
575  //
576  if( parser.OptionExists('E') ) {
577  LOG("gevgen_atmo", pINFO) << "Reading neutrino energy range";
578  string nue = parser.ArgAsString('E');
579 
580  // must be a comma separated set of values
581  if(nue.find(",") != string::npos) {
582  // split the comma separated list
583  vector<string> nurange = utils::str::Split(nue, ",");
584  assert(nurange.size() == 2);
585  double emin = atof(nurange[0].c_str());
586  double emax = atof(nurange[1].c_str());
587  assert(emax>emin && emin>=0.);
588  gOptEvMin = emin;
589  gOptEvMax = emax;
590  } else {
591  LOG("gevgen_atmo", pFATAL)
592  << "Invalid energy range. Use `-E emin,emax', eg `-E 0.5,100.";
593  PrintSyntax();
594  gAbortingInErr = true;
595  exit(1);
596  }
597  } else {
598  LOG("gevgen_atmo", pNOTICE)
599  << "No -e option. Using default energy range";
602  }
603 
604  //
605  // *** flux files
606  //
607  // syntax:
608  // simulation:/path/file.data[neutrino_code],/path/file.data[neutrino_code],...
609  //
610  if( parser.OptionExists('f') ) {
611  LOG("gevgen_atmo", pDEBUG) << "Getting input flux files";
612  string flux = parser.ArgAsString('f');
613 
614  // get flux simulation info (FLUKA,BGLRS) so as to instantiate the
615  // appropriate flux driver
616  string::size_type jsimend = flux.find_first_of(":",0);
617  if(jsimend==string::npos) {
618  LOG("gevgen_atmo", pFATAL)
619  << "You need to specify the flux file source";
620  PrintSyntax();
621  gAbortingInErr = true;
622  exit(1);
623  }
624  gOptFluxSim = flux.substr(0,jsimend);
625  for(string::size_type i=0; i<gOptFluxSim.size(); i++) {
626  gOptFluxSim[i] = toupper(gOptFluxSim[i]);
627  }
628  if((gOptFluxSim != "FLUKA") &&
629  (gOptFluxSim != "BGLRS") &&
630  (gOptFluxSim != "HAKKM")) {
631  LOG("gevgen_atmo", pFATAL)
632  << "The flux file source needs to be one of <FLUKA,BGLRS,HAKKM>";
633  PrintSyntax();
634  gAbortingInErr = true;
635  exit(1);
636  }
637  // now get the list of input files and the corresponding neutrino codes.
638  flux.erase(0,jsimend+1);
639  vector<string> fluxv = utils::str::Split(flux,",");
640  vector<string>::const_iterator fluxiter = fluxv.begin();
641  for( ; fluxiter != fluxv.end(); ++fluxiter) {
642  string filename_and_pdg = *fluxiter;
643  string::size_type open_bracket = filename_and_pdg.find("[");
644  string::size_type close_bracket = filename_and_pdg.find("]");
645  if (open_bracket ==string::npos ||
646  close_bracket==string::npos)
647  {
648  LOG("gevgen_atmo", pFATAL)
649  << "You made an error in specifying the flux info";
650  PrintSyntax();
651  gAbortingInErr = true;
652  exit(1);
653  }
654  string::size_type ibeg = 0;
655  string::size_type iend = open_bracket;
656  string::size_type jbeg = open_bracket+1;
657  string::size_type jend = close_bracket;
658  string flux_filename = filename_and_pdg.substr(ibeg,iend-ibeg);
659  string neutrino_pdg = filename_and_pdg.substr(jbeg,jend-jbeg);
660  gOptFluxFiles.insert(
661  map<int,string>::value_type(atoi(neutrino_pdg.c_str()), flux_filename));
662  }
663  if(gOptFluxFiles.size() == 0) {
664  LOG("gevgen_atmo", pFATAL)
665  << "You must specify at least one flux file!";
666  PrintSyntax();
667  gAbortingInErr = true;
668  exit(1);
669  }
670 
671  } else {
672  LOG("gevgen_atmo", pFATAL) << "No flux info was specified! Use the -f option.";
673  PrintSyntax();
674  gAbortingInErr = true;
675  exit(1);
676  }
677 
678  //
679  // *** geometry
680  //
681  string geom = "";
682  string lunits, dunits;
683  if( parser.OptionExists('g') ) {
684  LOG("gevgen_atmo", pDEBUG) << "Getting input geometry";
685  geom = parser.ArgAsString('g');
686 
687  // is it a ROOT file that contains a ROOT geometry?
688  bool accessible_geom_file =
689  utils::system::FileExists(geom.c_str());
690  if (accessible_geom_file) {
691  gOptRootGeom = geom;
692  gOptUsingRootGeom = true;
693  }
694  } else {
695  LOG("gevgen_atmo", pFATAL)
696  << "No geometry option specified - Exiting";
697  PrintSyntax();
698  gAbortingInErr = true;
699  exit(1);
700  } //-g
701 
702  if(gOptUsingRootGeom) {
703  // using a ROOT geometry - get requested geometry units
704 
705  // legth units:
706  if( parser.OptionExists('L') ) {
707  LOG("gevgen_atmo", pDEBUG)
708  << "Checking for input geometry length units";
709  lunits = parser.ArgAsString('L');
710  } else {
711  LOG("gevgen_atmo", pDEBUG) << "Using default geometry length units";
712  lunits = kDefOptGeomLUnits;
713  } // -L
714  // density units:
715  if( parser.OptionExists('D') ) {
716  LOG("gevgen_atmo", pDEBUG)
717  << "Checking for input geometry density units";
718  dunits = parser.ArgAsString('D');
719  } else {
720  LOG("gevgen_atmo", pDEBUG) << "Using default geometry density units";
721  dunits = kDefOptGeomDUnits;
722  } // -D
725 
726  // check whether an event generation volume name has been
727  // specified -- default is the 'top volume'
728  if( parser.OptionExists('t') ) {
729  LOG("gevgen_atmo", pDEBUG) << "Checking for input volume name";
730  gOptRootGeomTopVol = parser.ArgAsString('t');
731  } else {
732  LOG("gevgen_atmo", pDEBUG) << "Using the <master volume>";
733  } // -t
734 
735  // check whether an XML file with the maximum (density weighted)
736  // path lengths for each detector material is specified -
737  // otherwise will compute the max path lengths at job init
738  if( parser.OptionExists('m') ) {
739  LOG("gevgen_atmo", pDEBUG)
740  << "Checking for maximum path lengths XML file";
741  gOptExtMaxPlXml = parser.ArgAsString('m');
742  } else {
743  LOG("gevgen_atmo", pDEBUG)
744  << "Will compute the maximum path lengths at job init";
745  gOptExtMaxPlXml = "";
746  } // -m
747  } // using root geom?
748 
749  else {
750  // User has specified a target mix.
751  // Decode the list of target pdf codes & their corresponding weight fraction
752  // (specified as 'pdg_code_1[fraction_1],pdg_code_2[fraction_2],...')
753  // See documentation on top section of this file.
754  //
755  gOptTgtMix.clear();
756  vector<string> tgtmix = utils::str::Split(geom,",");
757  if(tgtmix.size()==1) {
758  int pdg = atoi(tgtmix[0].c_str());
759  double wgt = 1.0;
760  gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
761  } else {
762  vector<string>::const_iterator tgtmix_iter = tgtmix.begin();
763  for( ; tgtmix_iter != tgtmix.end(); ++tgtmix_iter) {
764  string tgt_with_wgt = *tgtmix_iter;
765  string::size_type open_bracket = tgt_with_wgt.find("[");
766  string::size_type close_bracket = tgt_with_wgt.find("]");
767  if (open_bracket ==string::npos ||
768  close_bracket==string::npos)
769  {
770  LOG("gevgen_atmo", pFATAL)
771  << "You made an error in specifying the target mix";
772  PrintSyntax();
773  gAbortingInErr = true;
774  exit(1);
775  }
776  string::size_type ibeg = 0;
777  string::size_type iend = open_bracket;
778  string::size_type jbeg = open_bracket+1;
779  string::size_type jend = close_bracket;
780  int pdg = atoi(tgt_with_wgt.substr(ibeg,iend-ibeg).c_str());
781  double wgt = atof(tgt_with_wgt.substr(jbeg,jend-jbeg).c_str());
782  LOG("gevgen_atmo", pDEBUG)
783  << "Adding to target mix: pdg = " << pdg << ", wgt = " << wgt;
784  gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
785 
786  }// tgtmix_iter
787  } // >1 materials in mix
788  } // using tgt mix?
789 
790  //
791  // Coordinate rotation matrix
792  //
793  gOptRot.SetToIdentity();
794  if( parser.OptionExists('R') ) {
795  string rotarg = parser.ArgAsString('R');
796  //get convention
797  string::size_type j = rotarg.find_first_of(":",0);
798  string convention = "";
799  if(j==string::npos) { convention = "X"; }
800  else { convention = rotarg.substr(0,j); }
801  //get angles phi,theta,psi
802  rotarg.erase(0,j+1);
803  vector<string> euler_angles = utils::str::Split(rotarg,",");
804  if(euler_angles.size() != 3) {
805  LOG("gevgen_atmo", pFATAL)
806  << "You didn't specify all 3 Euler angles using the -R option";
807  PrintSyntax();
808  gAbortingInErr = true;
809  exit(1);
810  }
811  double phi = atof(euler_angles[0].c_str());
812  double theta = atof(euler_angles[1].c_str());
813  double psi = atof(euler_angles[2].c_str());
814  //set Euler angles using appropriate convention
815  if(convention.find("X")!=string::npos ||
816  convention.find("x")!=string::npos)
817  {
818  LOG("gevgen_atmo", pNOTICE) << "Using X-convention for input Euler angles";
819  gOptRot.SetXEulerAngles(phi,theta,psi);
820  } else
821  if(convention.find("Y")!=string::npos ||
822  convention.find("y")!=string::npos)
823  {
824  LOG("gevgen_atmo", pNOTICE) << "Using Y-convention for input Euler angles";
825  gOptRot.SetYEulerAngles(phi,theta,psi);
826  } else {
827  LOG("gevgen_atmo", pFATAL)
828  << "Unknown Euler angle convention. Please use the X- or Y-convention";
829  PrintSyntax();
830  gAbortingInErr = true;
831  exit(1);
832  }
833  //invert?
834  if(convention.find("^-1")!=string::npos) {
835  LOG("gevgen_atmo", pNOTICE) << "Inverting rotation matrix";
836  gOptRot.Invert();
837  }
838  }
839 
840  //
841  // *** random number seed
842  //
843  if( parser.OptionExists("seed") ) {
844  LOG("gevgen_atmo", pINFO) << "Reading random number seed";
845  gOptRanSeed = parser.ArgAsLong("seed");
846  } else {
847  LOG("gevgen_atmo", pINFO) << "Unspecified random number seed - Using default";
848  gOptRanSeed = -1;
849  }
850 
851  //
852  // *** input cross-section file
853  //
854  if( parser.OptionExists("cross-sections") ) {
855  LOG("gevgen_atmo", pINFO) << "Reading cross-section file";
856  gOptInpXSecFile = parser.ArgAsString("cross-sections");
857  } else {
858  LOG("gevgen_atmo", pINFO) << "Unspecified cross-section file";
859  gOptInpXSecFile = "";
860  }
861 
862  //
863  // print-out summary
864  //
865 
866  PDGLibrary * pdglib = PDGLibrary::Instance();
867 
868  ostringstream gminfo;
869  if (gOptUsingRootGeom) {
870  gminfo << "Using ROOT geometry - file: " << gOptRootGeom
871  << ", top volume: "
872  << ((gOptRootGeomTopVol.size()==0) ? "<master volume>" : gOptRootGeomTopVol)
873  << ", max{PL} file: "
874  << ((gOptExtMaxPlXml.size()==0) ? "<none>" : gOptExtMaxPlXml)
875  << ", length units: " << lunits
876  << ", density units: " << dunits;
877  } else {
878  gminfo << "Using target mix - ";
879  map<int,double>::const_iterator iter;
880  for(iter = gOptTgtMix.begin(); iter != gOptTgtMix.end(); ++iter) {
881  int pdg_code = iter->first;
882  double wgt = iter->second;
883  TParticlePDG * p = pdglib->Find(pdg_code);
884  if(p) {
885  string name = p->GetName();
886  gminfo << "(" << name << ") -> " << 100*wgt << "% / ";
887  }//p?
888  }
889  }
890 
891  ostringstream fluxinfo;
892  fluxinfo << "Using " << gOptFluxSim << " flux files: ";
893  map<int,string>::const_iterator file_iter = gOptFluxFiles.begin();
894  for( ; file_iter != gOptFluxFiles.end(); ++file_iter) {
895  int neutrino_code = file_iter->first;
896  string filename = file_iter->second;
897  TParticlePDG * p = pdglib->Find(neutrino_code);
898  if(p) {
899  string name = p->GetName();
900  fluxinfo << "(" << name << ") -> " << filename << " / ";
901  }
902  }
903 
904  ostringstream expinfo;
905  if(gOptNev > 0) { expinfo << gOptNev << " events"; }
906  if(gOptKtonYrExposure > 0) { expinfo << gOptKtonYrExposure << " kton*yrs"; }
907 
908  ostringstream rotation;
909  rotation << "\t| " << gOptRot.XX() << " " << gOptRot.XY() << " " << gOptRot.XZ() << " |\n";
910  rotation << "\t| " << gOptRot.YX() << " " << gOptRot.YY() << " " << gOptRot.YZ() << " |\n";
911  rotation << "\t| " << gOptRot.ZX() << " " << gOptRot.ZY() << " " << gOptRot.ZZ() << " |\n";
912 
913  LOG("gevgen_atmo", pNOTICE)
914  << "\n\n"
915  << utils::print::PrintFramedMesg("gevgen_atmo job configuration");
916 
917  LOG("gevgen_atmo", pNOTICE)
918  << "\n"
919  << "\n @@ Run number: " << gOptRunNu
920  << "\n @@ Random number seed: " << gOptRanSeed
921  << "\n @@ Using cross-section file: " << gOptInpXSecFile
922  << "\n @@ Geometry"
923  << "\n\t" << gminfo.str()
924  << "\n @@ Flux"
925  << "\n\t" << fluxinfo.str()
926  << "\n @@ Exposure"
927  << "\n\t" << expinfo.str()
928  << "\n @@ Cuts"
929  << "\n\t Using energy range = (" << gOptEvMin << " GeV, " << gOptEvMax << " GeV)"
930  << "\n @@ Coordinate transformation (Rotation THZ -> User-defined coordinate system)"
931  << "\n" << rotation.str()
932  << "\n\n";
933 
934  //
935  // final checks
936  //
937  if(gOptKtonYrExposure > 0) {
938  LOG("gevgen_atmo", pFATAL)
939  << "\n Option to set exposure in terms of kton*yrs not supported just yet!"
940  << "\n Try the -n option instead";
941  PrintSyntax();
942  gAbortingInErr = true;
943  exit(1);
944  }
945 }
const XML_Char * name
Definition: expat.h:151
map< int, string > gOptFluxFiles
Definition: gAtmoEvGen.cxx:276
string gOptEvFilePrefix
Definition: gAtmoEvGen.cxx:288
bool FileExists(string filename)
Definition: SystemUtils.cxx:90
TRotation gOptRot
Definition: gAtmoEvGen.cxx:289
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
double gOptGeomDUnits
Definition: gAtmoEvGen.cxx:282
string filename
Definition: shutoffs.py:106
map< int, double > gOptTgtMix
Definition: gAtmoEvGen.cxx:278
double gOptEvMin
Definition: gAtmoEvGen.cxx:286
Loaders::FluxType flux
const double emin
double kDefOptEvMax
Definition: gAtmoEvGen.cxx:300
string kDefOptEvFilePrefix
Definition: gAtmoEvGen.cxx:296
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Definition: typedefs.hpp:11
string gOptFluxSim
Definition: gAtmoEvGen.cxx:275
string gOptInpXSecFile
Definition: gAtmoEvGen.cxx:291
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
string gOptRootGeomTopVol
Definition: gAtmoEvGen.cxx:280
const double emax
double UnitFromString(string u)
Definition: UnitUtils.cxx:26
int gOptNev
Definition: gAtmoEvGen.cxx:284
const double j
Definition: BetheBloch.cxx:29
#define pINFO
Definition: Messenger.h:63
const ana::Var wgt
string gOptRootGeom
Definition: gAtmoEvGen.cxx:279
Long_t gOptRunNu
Definition: gAtmoEvGen.cxx:274
string kDefOptGeomLUnits
Definition: gAtmoEvGen.cxx:297
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:42
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
double gOptEvMax
Definition: gAtmoEvGen.cxx:287
void geom(int which=0)
Definition: geom.C:163
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
exit(0)
assert(nhit_max >=nhit_nbins)
string PrintFramedMesg(string mesg, unsigned int nl=1, const char f='*')
Definition: PrintUtils.cxx:171
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:62
double gOptKtonYrExposure
Definition: gAtmoEvGen.cxx:285
string kDefOptGeomDUnits
Definition: gAtmoEvGen.cxx:298
bool gAbortingInErr
Definition: Messenger.cxx:56
double gOptGeomLUnits
Definition: gAtmoEvGen.cxx:281
void PrintSyntax(void)
Definition: gAtmoEvGen.cxx:947
string gOptExtMaxPlXml
Definition: gAtmoEvGen.cxx:283
bool gOptUsingRootGeom
Definition: gAtmoEvGen.cxx:277
#define pDEBUG
Definition: Messenger.h:64
long int gOptRanSeed
Definition: gAtmoEvGen.cxx:290
GFluxI * GetFlux ( void  )

Definition at line 441 of file gAtmoEvGen.cxx.

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

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

442 {
443  GFluxI * flux_driver = 0;
444 
445 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
446 
447  // Instantiate appropriate concrete flux driver
448  GAtmoFlux * atmo_flux_driver = 0;
449  if(gOptFluxSim == "FLUKA") {
450  GFLUKAAtmoFlux * fluka_flux = new GFLUKAAtmoFlux;
451  atmo_flux_driver = dynamic_cast<GAtmoFlux *>(fluka_flux);
452  } else
453  if(gOptFluxSim == "BGLRS") {
454  GBGLRSAtmoFlux * bartol_flux = new GBGLRSAtmoFlux;
455  atmo_flux_driver = dynamic_cast<GAtmoFlux *>(bartol_flux);
456  } else
457  if(gOptFluxSim == "HAKKM") {
458  GHAKKMAtmoFlux * honda_flux = new GHAKKMAtmoFlux;
459  atmo_flux_driver = dynamic_cast<GAtmoFlux *>(honda_flux);
460  } else {
461  LOG("gevgen_atmo", pFATAL) << "Uknonwn flux simulation: " << gOptFluxSim;
462  gAbortingInErr = true;
463  exit(1);
464  }
465  // Configure GAtmoFlux options (common to all concrete atmospheric flux drivers)
466  // set min/max energy:
467  atmo_flux_driver->ForceMinEnergy (gOptEvMin * units::GeV);
468  atmo_flux_driver->ForceMaxEnergy (gOptEvMax * units::GeV);
469  // set flux files:
470  map<int,string>::const_iterator file_iter = gOptFluxFiles.begin();
471  for( ; file_iter != gOptFluxFiles.end(); ++file_iter) {
472  int neutrino_code = file_iter->first;
473  string filename = file_iter->second;
474  atmo_flux_driver->AddFluxFile(neutrino_code, filename);
475  }
476  atmo_flux_driver->LoadFluxData();
477  // configure flux generation surface:
478  atmo_flux_driver->SetRadii(1, 1);
479  // set rotation for coordinate tranformation from the topocentric horizontal
480  // system to a user-defined coordinate system:
481  if(!gOptRot.IsIdentity()) {
482  atmo_flux_driver->SetUserCoordSystem(gOptRot);
483  }
484  // Cast to GFluxI, the generic flux driver interface
485  flux_driver = dynamic_cast<GFluxI *>(atmo_flux_driver);
486 
487 #else
488  LOG("gevgen_atmo", pFATAL) << "You need to enable the GENIE flux drivers first!";
489  LOG("gevgen_atmo", pFATAL) << "Use --enable-flux-drivers at the configuration step.";
490  gAbortingInErr = true;
491  exit(1);
492 #endif
493 
494  return flux_driver;
495 }
map< int, string > gOptFluxFiles
Definition: gAtmoEvGen.cxx:276
TRotation gOptRot
Definition: gAtmoEvGen.cxx:289
#define pFATAL
Definition: Messenger.h:57
string filename
Definition: shutoffs.py:106
void SetRadii(double Rlongitudinal, double Rtransverse)
Definition: GAtmoFlux.cxx:412
double gOptEvMin
Definition: gAtmoEvGen.cxx:286
A flux driver for the FLUKA 3-D Atmospheric Neutrino Flux.
string gOptFluxSim
Definition: gAtmoEvGen.cxx:275
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void ForceMaxEnergy(double emax)
Definition: GAtmoFlux.cxx:276
static constexpr Double_t GeV
Definition: Munits.h:291
void SetUserCoordSystem(TRotation &rotation)
Rotation: Topocentric Horizontal -> User-defined Topocentric Coord System.
Definition: GAtmoFlux.cxx:306
double gOptEvMax
Definition: gAtmoEvGen.cxx:287
void AddFluxFile(int neutrino_pdg, string filename)
Definition: GAtmoFlux.cxx:421
exit(0)
A driver for the HAKKM 3-D atmospheric neutrino flux (commonly known as the `Honda flux&#39;) ...
A base class for the FLUKA, BGLRS and ATMNC atmo. nu. flux drivers. The driver depends on data files ...
Definition: GAtmoFlux.h:61
void ForceMinEnergy(double emin)
Definition: GAtmoFlux.cxx:270
bool gAbortingInErr
Definition: Messenger.cxx:56
A flux driver for the Bartol Atmospheric Neutrino Flux.
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
GeomAnalyzerI * GetGeometry ( void  )

Definition at line 378 of file gAtmoEvGen.cxx.

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

Referenced by main(), and rawp::CheckDAQChannelMap::Reco().

379 {
380  GeomAnalyzerI * geom_driver = 0;
381 
382 #ifdef __GENIE_GEOM_DRIVERS_ENABLED__
383 
384  if(gOptUsingRootGeom) {
385  //
386  // *** Using a realistic root-based detector geometry description
387  //
388 
389  // creating & configuring a root geometry driver
392  rgeom -> SetLengthUnits (gOptGeomLUnits);
393  rgeom -> SetDensityUnits (gOptGeomDUnits);
394  rgeom -> SetTopVolName (gOptRootGeomTopVol);
395  // getting the bounding box dimensions along z so as to set the
396  // appropriate upstream generation surface for the JPARC flux driver
397  TGeoVolume * topvol = rgeom->GetGeometry()->GetTopVolume();
398  if(!topvol) {
399  LOG("gevgen_atmo", pFATAL) << " ** Null top ROOT geometry volume!";
400  gAbortingInErr = true;
401  exit(1);
402  }
403 /*
404  TGeoShape * bounding_box = topvol->GetShape();
405  bounding_box->GetAxisRange(3, zmin, zmax);
406  zmin *= rgeom->LengthUnits();
407  zmax *= rgeom->LengthUnits();
408 */
409  // switch on/off volumes as requested
410  if ( (gOptRootGeomTopVol[0] == '+') || (gOptRootGeomTopVol[0] == '-') ) {
411  bool exhaust = (*gOptRootGeomTopVol.c_str() == '+');
413  }
414 
415  // casting to the GENIE geometry driver interface
416  geom_driver = dynamic_cast<GeomAnalyzerI *> (rgeom);
417  }
418  else {
419  //
420  // *** Using a 'point' geometry with the specified target mix
421  // *** ( = a list of targets with their corresponding weight fraction)
422  //
423 
424  // creating & configuring a point geometry driver
427  // casting to the GENIE geometry driver interface
428  geom_driver = dynamic_cast<GeomAnalyzerI *> (pgeom);
429  }
430 
431 #else
432  LOG("gevgen_atmo", pFATAL) << "You need to enable the GENIE geometry drivers first!";
433  LOG("gevgen_atmo", pFATAL) << "Use --enable-geom-drivers at the configuration step.";
434  gAbortingInErr = true;
435  exit(1);
436 #endif
437 
438  return geom_driver;
439 }
#define pFATAL
Definition: Messenger.h:57
double gOptGeomDUnits
Definition: gAtmoEvGen.cxx:282
map< int, double > gOptTgtMix
Definition: gAtmoEvGen.cxx:278
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
string gOptRootGeomTopVol
Definition: gAtmoEvGen.cxx:280
A ROOT/GEANT4 geometry driver.
string gOptRootGeom
Definition: gAtmoEvGen.cxx:279
exit(0)
Defines the GENIE Geometry Analyzer Interface.
Definition: GeomAnalyzerI.h:30
bool gAbortingInErr
Definition: Messenger.cxx:56
void RecursiveExhaust(TGeoVolume *topvol, string volnames, bool exhaust)
Definition: GeoUtils.cxx:26
double gOptGeomLUnits
Definition: gAtmoEvGen.cxx:281
The PointGeomAnalyzer class is the simplest implementation of the GeomAnalyserI interface and defines...
bool gOptUsingRootGeom
Definition: gAtmoEvGen.cxx:277
virtual TGeoManager * GetGeometry(void) const
int main ( int  argc,
char **  argv 
)

Definition at line 303 of file gAtmoEvGen.cxx.

References genie::NtpWriter::AddEventRecord(), genie::RunOpt::BuildTune(), genie::utils::app_init::CacheFile(), genie::GMCJDriver::Configure(), genie::NtpWriter::CustomizeFilenamePrefix(), exit(), genie::GMCJDriver::ForceSingleProbScale(), genie::GMCJDriver::GenerateEvent(), GetCommandLineArgs(), GetFlux(), GetGeometry(), gOptEvFilePrefix, gOptInpXSecFile, gOptNev, gOptRanSeed, gOptRunNu, iev, genie::NtpWriter::Initialize(), genie::RunOpt::Instance(), kDefOptNtpFormat, LOG, genie::utils::app_init::MesgThresholds(), pFATAL, pNOTICE, genie::utils::app_init::RandGen(), genie::NtpWriter::Save(), genie::GMCJDriver::SetEventGeneratorList(), genie::GHepRecord::SetPrintLevel(), genie::GMCJMonitor::SetRefreshRate(), genie::GMCJMonitor::Update(), genie::GMCJDriver::UseFluxDriver(), genie::GMCJDriver::UseGeomAnalyzer(), genie::GMCJDriver::UseSplines(), and genie::utils::app_init::XSecTable().

304 {
305  // Parse command line arguments
306  GetCommandLineArgs(argc,argv);
307 
308  if ( ! RunOpt::Instance()->Tune() ) {
309  LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
310  exit(-1);
311  }
312  RunOpt::Instance()->BuildTune();
313 
314  // Iinitialization of random number generators, cross-section table, messenger, cache etc...
315  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
316  utils::app_init::CacheFile(RunOpt::Instance()->CacheFile());
319 
320  // get flux driver
321  GFluxI * flux_driver = GetFlux();
322 
323  // get geometry driver
324  GeomAnalyzerI * geom_driver = GetGeometry();
325 
326  // create the GENIE monte carlo job driver
327  GMCJDriver* mcj_driver = new GMCJDriver;
328  mcj_driver->SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
329  mcj_driver->UseFluxDriver(flux_driver);
330  mcj_driver->UseGeomAnalyzer(geom_driver);
331  mcj_driver->Configure();
332  mcj_driver->UseSplines();
333  mcj_driver->ForceSingleProbScale();
334 
335  // initialize an ntuple writer
337  ntpw.CustomizeFilenamePrefix(gOptEvFilePrefix);
338  ntpw.Initialize();
339 
340  // Create a MC job monitor for a periodically updated status file
341  GMCJMonitor mcjmonitor(gOptRunNu);
342  mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());
343 
344  // Set GHEP print level
345  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
346 
347  // event loop
348  for(int iev = 0; iev < gOptNev; iev++) {
349 
350  // generate next event
351  EventRecord* event = mcj_driver->GenerateEvent();
352 
353  // set weight (if using a weighted flux)
354  //event->SetWeight(event->Weight()*flux_driver->Weight());
355 
356  // print-out
357  LOG("gevgen_atmo", pNOTICE) << "Generated event: " << *event;
358 
359  // save the event, refresh the mc job monitor
360  ntpw.AddEventRecord(iev, event);
361  mcjmonitor.Update(iev,event);
362 
363  // clean-up
364  delete event;
365  }
366 
367  // save the event file
368  ntpw.Save();
369 
370  // clean-up
371  delete geom_driver;
372  delete flux_driver;
373  delete mcj_driver;
374 
375  return 0;
376 }
int iev
Definition: runWimpSim.h:118
void RandGen(long int seed)
Definition: AppInit.cxx:31
void XSecTable(string inpfile, bool require_table)
Definition: AppInit.cxx:39
string gOptEvFilePrefix
Definition: gAtmoEvGen.cxx:288
void SetEventGeneratorList(string listname)
Definition: GMCJDriver.cxx:157
#define pFATAL
Definition: Messenger.h:57
GFluxI * GetFlux(void)
Definition: gAtmoEvGen.cxx:441
void UseFluxDriver(GFluxI *flux)
Definition: GMCJDriver.cxx:174
Simple class to create & update MC job status files and env. vars. This is used to be able to keep tr...
Definition: GMCJMonitor.h:30
A GENIE `MC Job Driver&#39;. Can be used for setting up complicated event generation cases involving deta...
Definition: GMCJDriver.h:47
string gOptInpXSecFile
Definition: gAtmoEvGen.cxx:291
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void ForceSingleProbScale(void)
Definition: GMCJDriver.cxx:219
int gOptNev
Definition: gAtmoEvGen.cxx:284
void Configure(bool calc_prob_scales=true)
Definition: GMCJDriver.cxx:446
EventRecord * GenerateEvent(void)
Definition: GMCJDriver.cxx:838
Long_t gOptRunNu
Definition: gAtmoEvGen.cxx:274
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
A utility class to facilitate creating the GENIE MC Ntuple from the output GENIE GHEP event records...
Definition: NtpWriter.h:40
exit(0)
NtpMCFormat_t kDefOptNtpFormat
Definition: gAtmoEvGen.cxx:295
A vector of EventGeneratorI objects.
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:100
#define pNOTICE
Definition: Messenger.h:62
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:497
Defines the GENIE Geometry Analyzer Interface.
Definition: GeomAnalyzerI.h:30
void UseGeomAnalyzer(GeomAnalyzerI *geom)
Definition: GMCJDriver.cxx:179
void UseSplines(bool useLogE=true)
Definition: GMCJDriver.cxx:184
void CacheFile(string inpfile)
Definition: AppInit.cxx:118
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
GeomAnalyzerI * GetGeometry(void)
Definition: gAtmoEvGen.cxx:378
long int gOptRanSeed
Definition: gAtmoEvGen.cxx:290
void PrintSyntax ( void  )

Definition at line 947 of file gAtmoEvGen.cxx.

References LOG, and pFATAL.

Referenced by GetCommandLineArgs().

948 {
949  LOG("gevgen_atmo", pFATAL)
950  << "\n **Syntax**"
951  << "\n gevgen_atmo [-h]"
952  << "\n [-r run#]"
953  << "\n -f simulation:flux_file[neutrino_code],..."
954  << "\n -g geometry"
955  << "\n [-R coordinate_rotation_matrix]"
956  << "\n [-t geometry_top_volume_name]"
957  << "\n [-m max_path_lengths_xml_file]"
958  << "\n [-L geometry_length_units]"
959  << "\n [-D geometry_density_units]"
960  << "\n <-n n_of_events,"
961  << "\n -e exposure_in_kton_x_yrs>"
962  << "\n -E min_energy,max_energy"
963  << "\n [-o output_event_file_prefix]"
964  << "\n [--seed random_number_seed]"
965  << "\n --cross-sections xml_file"
966  << "\n [--event-generator-list list_name]"
967  << "\n [--message-thresholds xml_file]"
968  << "\n [--unphysical-event-mask mask]"
969  << "\n [--event-record-print-level level]"
970  << "\n [--mc-job-status-refresh-rate rate]"
971  << "\n [--cache-file root_file]"
972  << "\n"
973  << " Please also read the detailed documentation at http://www.genie-mc.org"
974  << "\n";
975 }
#define pFATAL
Definition: Messenger.h:57
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97

Variable Documentation

string gOptEvFilePrefix

Definition at line 288 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

double gOptEvMax

Definition at line 287 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and GetFlux().

double gOptEvMin

Definition at line 286 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and GetFlux().

string gOptExtMaxPlXml

Definition at line 283 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs().

map<int,string> gOptFluxFiles

Definition at line 276 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and GetFlux().

string gOptFluxSim

Definition at line 275 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and GetFlux().

double gOptGeomDUnits = 0

Definition at line 282 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and GetGeometry().

double gOptGeomLUnits = 0

Definition at line 281 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and GetGeometry().

string gOptInpXSecFile

Definition at line 291 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

double gOptKtonYrExposure = -1

Definition at line 285 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs().

int gOptNev = -1

Definition at line 284 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

long int gOptRanSeed

Definition at line 290 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

string gOptRootGeom

Definition at line 279 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and GetGeometry().

string gOptRootGeomTopVol = ""

Definition at line 280 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and GetGeometry().

TRotation gOptRot

Definition at line 289 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and GetFlux().

Long_t gOptRunNu

Definition at line 274 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

map<int,double> gOptTgtMix

Definition at line 278 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and GetGeometry().

bool gOptUsingRootGeom = false

Definition at line 277 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs(), and GetGeometry().

string kDefOptEvFilePrefix = "gntp"

Definition at line 296 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs().

double kDefOptEvMax = 50.0

Definition at line 300 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs().

double kDefOptEvMin = 0.5

Definition at line 299 of file gAtmoEvGen.cxx.

string kDefOptGeomDUnits = "g_cm3"

Definition at line 298 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs().

string kDefOptGeomLUnits = "mm"

Definition at line 297 of file gAtmoEvGen.cxx.

Referenced by GetCommandLineArgs().

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 295 of file gAtmoEvGen.cxx.

Referenced by main().