131 #if defined(HAVE_FENV_H) && defined(HAVE_FEENABLEEXCEPT) 138 #include <TVector3.h> 166 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__ 167 #ifdef __GENIE_GEOM_DRIVERS_ENABLED__ 168 #define __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__ 178 using std::ostringstream;
180 using namespace genie;
187 #ifdef __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__ 188 void GenerateEventsUsingFluxOrTgtMix();
190 GFluxI * FluxDriver (
void);
191 GFluxI * MonoEnergeticFluxDriver (
void);
192 GFluxI * TH1FluxDriver (
void);
224 #if defined(HAVE_FENV_H) && defined(HAVE_FEENABLEEXCEPT) 225 feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
232 #ifdef __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__ 233 GenerateEventsUsingFluxOrTgtMix();
236 <<
"\n To be able to generate neutrino events from a flux and/or a target mix" 237 <<
"\n you need to add the following config options at your GENIE installation:" 238 <<
"\n --enable-flux-drivers --enable-geom-drivers \n" ;
250 LOG(
"gevgen",
pFATAL) <<
" No TuneId in RunOption";
271 TLorentzVector nu_p4(0.,0.,Ev,Ev);
303 <<
"\n ** Will generate " <<
gOptNevents <<
" events for \n" 304 << init_state <<
" at Ev = " << Ev <<
" GeV";
310 <<
" *** Generating event............ " << ievent;
317 <<
"Last attempt failed. Re-trying....";
322 <<
"Generated Event GHEP Record: " << *event;
326 mcjmonitor.
Update(ievent,event);
336 #ifdef __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__ 338 void GenerateEventsUsingFluxOrTgtMix(
void)
341 GFluxI * flux_driver = FluxDriver();
378 LOG(
"gevgen",
pNOTICE) <<
" *** Generating event............ " << ievent;
383 LOG(
"gevgen",
pNOTICE) <<
"Generated Event GHEP Record: " << *event;
387 mcjmonitor.
Update(ievent,event);
415 else flux_driver = TH1FluxDriver();
420 GFluxI * MonoEnergeticFluxDriver(
void)
430 GFluxI * TH1FluxDriver(
void)
437 int flux_entries = 100000;
445 bool input_is_text_file = ! gSystem->AccessPathName(
gOptFlux.c_str());
446 bool input_is_root_file =
gOptFlux.find(
".root") != string::npos &&
448 if(input_is_text_file) {
454 double estep = (emax-
emin)/(n-1);
455 double ymax = -1, ry = -1, gy = -1,
e = -1;
456 for(
int i=0;
i<
n;
i++) {
458 ymax = TMath::Max(ymax, input_flux->
Evaluate(
e));
463 spectrum =
new TH1D(
"spectrum",
"neutrino flux", 300, emin, emax);
464 spectrum->SetDirectory(0);
466 for(
int ientry=0; ientry<flux_entries; ientry++) {
472 LOG(
"gevgen",
pFATAL) <<
"Couldn't generate a flux histogram";
475 e = emin + de * r->
RndGen().Rndm();
476 gy = ymax * r->
RndGen().Rndm();
479 if(accept) spectrum->Fill(
e);
484 else if(input_is_root_file) {
490 assert( !gSystem->AccessPathName(fv[0].c_str()) );
492 LOG(
"gevgen",
pNOTICE) <<
"Getting input flux from root file: " << fv[0];
493 TFile * flux_file =
new TFile(fv[0].c_str(),
"read");
495 LOG(
"gevgen",
pNOTICE) <<
"Flux name: " << fv[1];
496 TH1D * hst = (TH1D *)flux_file->Get(fv[1].c_str());
502 spectrum = (TH1D*)hst->Clone();
503 spectrum->SetNameTitle(
"spectrum",
"neutrino_flux");
504 spectrum->SetDirectory(0);
506 if(hst->GetBinLowEdge(
ibin) + hst->GetBinWidth(
ibin) > emax ||
508 spectrum->SetBinContent(
ibin, 0);
512 LOG(
"gevgen",
pNOTICE) << spectrum->GetEntries();
517 LOG(
"gevgen",
pNOTICE) << spectrum->GetEntries();
524 spectrum =
new TH1D(
"spectrum",
"neutrino flux", 300, emin, emax);
525 spectrum->SetDirectory(0);
526 spectrum->FillRandom(
"input_func", flux_entries);
531 TFile
f(
"./input-flux.root",
"recreate");
535 TVector3 bdir (0,0,1);
536 TVector3 bspot(0,0,0);
552 LOG(
"gevgen",
pINFO) <<
"Parsing command line arguments";
571 LOG(
"gevgen",
pINFO) <<
"Reading number of events to generate";
575 <<
"Unspecified number of events to generate - Using default";
581 LOG(
"gevgen",
pINFO) <<
"Reading MC run number";
584 LOG(
"gevgen",
pINFO) <<
"Unspecified run number - Using default";
590 LOG(
"gevgen",
pINFO) <<
"Reading output file name";
602 bool using_flux =
false;
604 LOG(
"gevgen",
pINFO) <<
"Reading flux function";
611 <<
"-s option no longer available. Please read the revised code documentation";
622 LOG(
"gevgen",
pINFO) <<
"Reading neutrino energy";
626 if(nue.find(
",") != string::npos) {
629 assert(nurange.size() == 2);
630 double emin = atof(nurange[0].c_str());
631 double emax = atof(nurange[1].c_str());
632 assert(emax>emin && emin>=0);
637 <<
"No flux was specified but an energy range was input!";
639 <<
"Events will be generated at fixed E = " <<
gOptNuEnergy <<
" GeV";
647 LOG(
"gevgen",
pFATAL) <<
"Unspecified neutrino energy - Exiting";
654 LOG(
"gevgen",
pINFO) <<
"Reading neutrino PDG code";
657 LOG(
"gevgen",
pFATAL) <<
"Unspecified neutrino PDG code - Exiting";
663 bool using_tgtmix =
false;
665 LOG(
"gevgen",
pINFO) <<
"Reading target mix";
669 if(tgtmix.size()==1) {
670 int pdg = atoi(tgtmix[0].c_str());
672 gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
675 vector<string>::const_iterator tgtmix_iter = tgtmix.begin();
676 for( ; tgtmix_iter != tgtmix.end(); ++tgtmix_iter) {
677 string tgt_with_wgt = *tgtmix_iter;
684 int pdg = atoi(tgt_with_wgt.substr(ibeg,iend).c_str());
685 double wgt = atof(tgt_with_wgt.substr(jbeg,jend).c_str());
687 <<
"Adding to target mix: pdg = " << pdg <<
", wgt = " <<
wgt;
688 gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
693 LOG(
"gevgen",
pFATAL) <<
"Unspecified target PDG code - Exiting";
702 LOG(
"gevgen",
pINFO) <<
"Reading random number seed";
705 LOG(
"gevgen",
pINFO) <<
"Unspecified random number seed - Using default";
711 LOG(
"gevgen",
pINFO) <<
"Reading cross-section file";
714 LOG(
"gevgen",
pINFO) <<
"Unspecified cross-section file";
731 <<
"Random number seed was not set, using default";
740 <<
"No input cross-section spline file";
748 <<
"Neutrino energy: [" 757 <<
"Target code (PDG) & weight fraction (in case of multiple targets): ";
758 map<int,double>::const_iterator iter;
760 int tgtpdgc = iter->first;
761 double wgt = iter->second;
763 <<
" >> " << tgtpdgc <<
" (weight fraction = " << wgt <<
")";
774 <<
"\n\n" <<
"Syntax:" <<
"\n" 778 <<
"\n -e energy (or energy range) " 779 <<
"\n -p neutrino_pdg" 780 <<
"\n -t target_pdg " 781 <<
"\n [-f flux_description]" 782 <<
"\n [-o outfile_name]" 784 <<
"\n [--seed random_number_seed]" 785 <<
"\n [--cross-sections xml_file]" 786 <<
"\n [--event-generator-list list_name]" 787 <<
"\n [--message-thresholds xml_file]" 788 <<
"\n [--unphysical-event-mask mask]" 789 <<
"\n [--event-record-print-level level]" 790 <<
"\n [--mc-job-status-refresh-rate rate]" 791 <<
"\n [--cache-file root_file]" 792 <<
"\n [--xml-path config_xml_dir]"
void RandGen(long int seed)
static void SetPrintLevel(int print_level)
void XSecTable(string inpfile, bool require_table)
void SetUnphysEventMask(const TBits &mask)
string ArgAsString(char opt)
THE MAIN GENIE PROJECT NAMESPACE
static RandomGen * Instance()
Access instance.
void SetEventGeneratorList(string listname)
void CustomizeFilename(string filename)
void ReadFromCommandLine(int argc, char **argv)
NtpMCFormat_t kDefOptNtpFormat
A numeric analysis tool class for interpolating 1-D functions.
void Update(int iev, const EventRecord *event)
void CustomizeFilename(string filename)
double Evaluate(double x) const
A singleton holding random number generator classes. All random number generation in GENIE should tak...
void SetRefreshRate(int rate)
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
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
void SetEventGeneratorList(string listname)
void ForceSingleProbScale(void)
A generic GENIE flux driver. Generates a 'cylindrical' neutrino beam along the input direction...
A simple GENIE flux driver for monoenergetic neutrinos along the z direction. Can handle a mix of neu...
void Configure(bool calc_prob_scales=true)
void Save(void)
get the even tree
GENIE Event Generation Driver. A minimalist user interface object for generating neutrino interaction...
void AddEventRecord(int ievent, const EventRecord *ev_rec)
save the event tree
void SetTransverseRadius(double Rt)
Misc GENIE control constants.
void SetNuDirection(const TVector3 &direction)
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...
int main(int argc, char **argv)
void Initialize(void)
add event
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...
map< int, double > gOptTgtMix
bool gOptUsingFluxOrTgtMix
void GetCommandLineArgs(int argc, char **argv)
TRandom3 & RndGen(void) const
rnd number generator for generic usage
void Configure(int nu_pdgc, int Z, int A)
assert(nhit_max >=nhit_nbins)
string PrintFramedMesg(string mesg, unsigned int nl=1, const char f='*')
A vector of EventGeneratorI objects.
void GenerateEventsAtFixedInitState(void)
void SetBeamSpot(const TVector3 &spot)
void MesgThresholds(string inpfile)
void AddEnergySpectrum(int nu_pdgc, TH1D *spectrum)
Command line argument parser.
enum genie::ENtpMCFormat NtpMCFormat_t
Defines the GENIE Geometry Analyzer Interface.
static const unsigned int kRjMaxIterations
void SetUnphysEventMask(const TBits &mask)
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)
void EnableBareXSecPreCalc(bool flag)
Initial State information.
GENIE Interface for user-defined flux classes.
EventRecord * GenerateEvent(const TLorentzVector &nu4p)