225 #include <TRotation.h> 247 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__ 253 #ifdef __GENIE_GEOM_DRIVERS_ENABLED__ 262 using std::ostringstream;
264 using namespace genie;
309 LOG(
"gmkspl",
pFATAL) <<
" No TuneId in RunOption";
357 LOG(
"gevgen_atmo",
pNOTICE) <<
"Generated event: " << *event;
382 #ifdef __GENIE_GEOM_DRIVERS_ENABLED__ 397 TGeoVolume * topvol = rgeom->
GetGeometry()->GetTopVolume();
399 LOG(
"gevgen_atmo",
pFATAL) <<
" ** Null top ROOT geometry volume!";
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.";
445 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__ 451 atmo_flux_driver =
dynamic_cast<GAtmoFlux *
>(fluka_flux);
455 atmo_flux_driver =
dynamic_cast<GAtmoFlux *
>(bartol_flux);
459 atmo_flux_driver =
dynamic_cast<GAtmoFlux *
>(honda_flux);
470 map<int,string>::const_iterator file_iter =
gOptFluxFiles.begin();
472 int neutrino_code = file_iter->first;
473 string filename = file_iter->second;
474 atmo_flux_driver->
AddFluxFile(neutrino_code, filename);
485 flux_driver =
dynamic_cast<GFluxI *
>(atmo_flux_driver);
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.";
503 LOG(
"gevgen_atmo",
pNOTICE) <<
"Parsing command line arguments";
518 LOG(
"gevgen_atmo",
pDEBUG) <<
"Reading MC run number";
521 LOG(
"gevgen_atmo",
pDEBUG) <<
"Unspecified run number - Using default";
530 bool have_required_statistics =
false;
533 <<
"Reading number of events to generate";
535 have_required_statistics =
true;
539 if(have_required_statistics) {
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";
548 <<
"Reading requested exposure in kton*yrs";
550 have_required_statistics =
true;
552 if(!have_required_statistics) {
554 <<
"You must request exposure either in terms of number of events and kton*yrs" 555 <<
"\nUse any of the -n, -e options";
565 LOG(
"gevgen_atmo",
pDEBUG) <<
"Reading the event filename prefix";
569 <<
"Will set the default event filename prefix";
577 LOG(
"gevgen_atmo",
pINFO) <<
"Reading neutrino energy range";
581 if(nue.find(
",") != string::npos) {
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.);
592 <<
"Invalid energy range. Use `-E emin,emax', eg `-E 0.5,100.";
599 <<
"No -e option. Using default energy range";
611 LOG(
"gevgen_atmo",
pDEBUG) <<
"Getting input flux files";
617 if(jsimend==string::npos) {
619 <<
"You need to specify the flux file source";
632 <<
"The flux file source needs to be one of <FLUKA,BGLRS,HAKKM>";
638 flux.erase(0,jsimend+1);
640 vector<string>::const_iterator fluxiter = fluxv.begin();
641 for( ; fluxiter != fluxv.end(); ++fluxiter) {
642 string filename_and_pdg = *fluxiter;
645 if (open_bracket ==string::npos ||
646 close_bracket==string::npos)
649 <<
"You made an error in specifying the flux info";
658 string flux_filename = filename_and_pdg.substr(ibeg,iend-ibeg);
659 string neutrino_pdg = filename_and_pdg.substr(jbeg,jend-jbeg);
661 map<int,string>::value_type(atoi(neutrino_pdg.c_str()), flux_filename));
665 <<
"You must specify at least one flux file!";
672 LOG(
"gevgen_atmo",
pFATAL) <<
"No flux info was specified! Use the -f option.";
682 string lunits, dunits;
684 LOG(
"gevgen_atmo",
pDEBUG) <<
"Getting input geometry";
688 bool accessible_geom_file =
690 if (accessible_geom_file) {
696 <<
"No geometry option specified - Exiting";
708 <<
"Checking for input geometry length units";
711 LOG(
"gevgen_atmo",
pDEBUG) <<
"Using default geometry length units";
717 <<
"Checking for input geometry density units";
720 LOG(
"gevgen_atmo",
pDEBUG) <<
"Using default geometry density units";
729 LOG(
"gevgen_atmo",
pDEBUG) <<
"Checking for input volume name";
732 LOG(
"gevgen_atmo",
pDEBUG) <<
"Using the <master volume>";
740 <<
"Checking for maximum path lengths XML file";
744 <<
"Will compute the maximum path lengths at job init";
757 if(tgtmix.size()==1) {
758 int pdg = atoi(tgtmix[0].c_str());
760 gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
762 vector<string>::const_iterator tgtmix_iter = tgtmix.begin();
763 for( ; tgtmix_iter != tgtmix.end(); ++tgtmix_iter) {
764 string tgt_with_wgt = *tgtmix_iter;
767 if (open_bracket ==string::npos ||
768 close_bracket==string::npos)
771 <<
"You made an error in specifying the target mix";
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());
783 <<
"Adding to target mix: pdg = " << pdg <<
", wgt = " <<
wgt;
784 gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
798 string convention =
"";
799 if(j==string::npos) { convention =
"X"; }
800 else { convention = rotarg.substr(0,j); }
804 if(euler_angles.size() != 3) {
806 <<
"You didn't specify all 3 Euler angles using the -R option";
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());
815 if(convention.find(
"X")!=string::npos ||
816 convention.find(
"x")!=string::npos)
818 LOG(
"gevgen_atmo",
pNOTICE) <<
"Using X-convention for input Euler angles";
819 gOptRot.SetXEulerAngles(phi,theta,psi);
821 if(convention.find(
"Y")!=string::npos ||
822 convention.find(
"y")!=string::npos)
824 LOG(
"gevgen_atmo",
pNOTICE) <<
"Using Y-convention for input Euler angles";
825 gOptRot.SetYEulerAngles(phi,theta,psi);
828 <<
"Unknown Euler angle convention. Please use the X- or Y-convention";
834 if(convention.find(
"^-1")!=string::npos) {
835 LOG(
"gevgen_atmo",
pNOTICE) <<
"Inverting rotation matrix";
844 LOG(
"gevgen_atmo",
pINFO) <<
"Reading random number seed";
847 LOG(
"gevgen_atmo",
pINFO) <<
"Unspecified random number seed - Using default";
855 LOG(
"gevgen_atmo",
pINFO) <<
"Reading cross-section file";
858 LOG(
"gevgen_atmo",
pINFO) <<
"Unspecified cross-section file";
868 ostringstream gminfo;
870 gminfo <<
"Using ROOT geometry - file: " <<
gOptRootGeom 873 <<
", max{PL} file: " 875 <<
", length units: " << lunits
876 <<
", density units: " << dunits;
878 gminfo <<
"Using target mix - ";
879 map<int,double>::const_iterator iter;
881 int pdg_code = iter->first;
882 double wgt = iter->second;
883 TParticlePDG *
p = pdglib->
Find(pdg_code);
885 string name = p->GetName();
886 gminfo <<
"(" << name <<
") -> " << 100*wgt <<
"% / ";
891 ostringstream fluxinfo;
892 fluxinfo <<
"Using " <<
gOptFluxSim <<
" flux files: ";
893 map<int,string>::const_iterator file_iter =
gOptFluxFiles.begin();
895 int neutrino_code = file_iter->first;
896 string filename = file_iter->second;
897 TParticlePDG *
p = pdglib->
Find(neutrino_code);
899 string name = p->GetName();
900 fluxinfo <<
"(" << name <<
") -> " << filename <<
" / ";
904 ostringstream expinfo;
908 ostringstream rotation;
923 <<
"\n\t" << gminfo.str()
925 <<
"\n\t" << fluxinfo.str()
927 <<
"\n\t" << expinfo.str()
930 <<
"\n @@ Coordinate transformation (Rotation THZ -> User-defined coordinate system)" 931 <<
"\n" << rotation.str()
939 <<
"\n Option to set exposure in terms of kton*yrs not supported just yet!" 940 <<
"\n Try the -n option instead";
951 <<
"\n gevgen_atmo [-h]" 953 <<
"\n -f simulation:flux_file[neutrino_code],..." 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]" 973 <<
" Please also read the detailed documentation at http://www.genie-mc.org"
void RandGen(long int seed)
static void SetPrintLevel(int print_level)
void XSecTable(string inpfile, bool require_table)
double ArgAsDouble(char opt)
map< int, string > gOptFluxFiles
string ArgAsString(char opt)
THE MAIN GENIE PROJECT NAMESPACE
void SetEventGeneratorList(string listname)
bool FileExists(string filename)
void ReadFromCommandLine(int argc, char **argv)
void Update(int iev, const EventRecord *event)
void SetRadii(double Rlongitudinal, double Rtransverse)
map< int, double > gOptTgtMix
int main(int argc, char **argv)
void SetRefreshRate(int rate)
string kDefOptEvFilePrefix
void UseFluxDriver(GFluxI *flux)
Simple class to create & update MC job status files and env. vars. This is used to be able to keep tr...
A GENIE `MC Job Driver'. Can be used for setting up complicated event generation cases involving deta...
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
A flux driver for the FLUKA 3-D Atmospheric Neutrino Flux.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
void ForceMaxEnergy(double emax)
string gOptRootGeomTopVol
void ForceSingleProbScale(void)
double UnitFromString(string u)
void Configure(bool calc_prob_scales=true)
A ROOT/GEANT4 geometry driver.
void Save(void)
get the even tree
void AddEventRecord(int ievent, const EventRecord *ev_rec)
save the event tree
void BuildTune()
build tune and inform XSecSplineList
EventRecord * GenerateEvent(void)
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
void SetUserCoordSystem(TRotation &rotation)
Rotation: Topocentric Horizontal -> User-defined Topocentric Coord System.
void CustomizeFilenamePrefix(string prefix)
void Initialize(void)
add event
static PDGLibrary * Instance(void)
static RunOpt * Instance(void)
vector< string > Split(string input, string delim)
A utility class to facilitate creating the GENIE MC Ntuple from the output GENIE GHEP event records...
Singleton class to load & serve a TDatabasePDG.
TParticlePDG * Find(int pdgc)
void AddFluxFile(int neutrino_pdg, string filename)
assert(nhit_max >=nhit_nbins)
A driver for the HAKKM 3-D atmospheric neutrino flux (commonly known as the `Honda flux') ...
NtpMCFormat_t kDefOptNtpFormat
string PrintFramedMesg(string mesg, unsigned int nl=1, const char f='*')
A vector of EventGeneratorI objects.
void MesgThresholds(string inpfile)
Command line argument parser.
enum genie::ENtpMCFormat NtpMCFormat_t
double gOptKtonYrExposure
void GetCommandLineArgs(int argc, char **argv)
Defines the GENIE Geometry Analyzer Interface.
A base class for the FLUKA, BGLRS and ATMNC atmo. nu. flux drivers. The driver depends on data files ...
void ForceMinEnergy(double emin)
A flux driver for the Bartol Atmospheric Neutrino Flux.
void RecursiveExhaust(TGeoVolume *topvol, string volnames, bool exhaust)
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
void UseGeomAnalyzer(GeomAnalyzerI *geom)
void UseSplines(bool useLogE=true)
The PointGeomAnalyzer class is the simplest implementation of the GeomAnalyserI interface and defines...
bool OptionExists(char opt)
was option set?
void CacheFile(string inpfile)
virtual TGeoManager * GetGeometry(void) const
GENIE Interface for user-defined flux classes.
GeomAnalyzerI * GetGeometry(void)