gAtmoEvGen.cxx
Go to the documentation of this file.
1 //________________________________________________________________________________________
2 /*!
3 
4 \program gevgen_atmo
5 
6 \brief A GENIE atmospheric neutrino event generation application.
7 
8  *** Synopsis :
9 
10  gevgen_atmo [-h]
11  [-r run#]
12  -f flux
13  -g geometry
14  [-R coordinate_rotation_matrix]
15  [-t geometry_top_volume_name]
16  [-m max_path_lengths_xml_file]
17  [-L geometry_length_units]
18  [-D geometry_density_units]
19  <-n n_of_events,
20  -e exposure_in_kton_x_yrs >
21  -E min_energy,max_energy
22  [-o output_event_file_prefix]
23  [--seed random_number_seed]
24  [--cross-sections xml_file]
25  [--event-generator-list list_name]
26  [--tune genie_tune]
27  [--message-thresholds xml_file]
28  [--unphysical-event-mask mask]
29  [--event-record-print-level level]
30  [--mc-job-status-refresh-rate rate]
31  [--cache-file root_file]
32 
33  *** Options :
34 
35  [] Denotes an optional argument.
36  <> Denotes a set of arguments out of which only one can be set.
37 
38  -h
39  Prints out the syntax and exits
40  -r
41  Specifies the MC run number
42  [default: 100000000]
43  -f
44  Specifies the input flux files
45  The general syntax is: `-f simulation:/path/file.data[neutrino_code],...'
46  [Notes]
47  - The `simulation' string can be either `FLUKA', `BGLRS' or `HAKKM'.
48  See:
49  - $GENIE/src/Flux/GFLUKAAtmoFlux.h
50  - $GENIE/src/Flux/GBGLRSAtmoFlux.h
51  - $GENIE/src/Flux/GHAKKMAtmoFlux.h
52  - The neutrino codes are the PDG ones.
53  - The /path/file.data,neutrino_code part of the option can be
54  repeated multiple times (separated by commas), once for each
55  flux neutrino species you want to consider,
56  eg. '-f FLUKA:~/data/sdave_numu07.dat[14],~/data/sdave_nue07.dat[12]'
57  eg. '-f BGLRS:~/data/flux10_271003_z.kam_nue[12]'
58  eg. '-f HAKKM:~/data/kam-ally-20-12-solmax.d[12]'
59  -g
60  Input 'geometry'.
61  This option can be used to specify any of:
62  1 > A ROOT file containing a ROOT/GEANT geometry description
63  [Examples]
64  - To use the master volume from the ROOT geometry stored
65  in the nd280-geom.root file, type:
66  '-g /some/path/nd280-geom.root'
67  2 > A mix of target materials, each with its corresponding weight,
68  typed as a comma-separated list of nuclear pdg codes (in the
69  std PDG2006 convention: 10LZZZAAAI) with the weight fractions
70  in brackets, eg code1[fraction1],code2[fraction2],...
71  If that option is used (no detailed input geometry description)
72  then the interaction vertices are distributed in the detector
73  by the detector MC.
74  [Examples]
75  - To use a target mix of 89% O16 and 11% H, type:
76  '-g 1000080160[0.89],1000010010[0.11]'
77  - To use a target which is 100% C12, type:
78  '-g 1000060120'
79  -R
80  Input rotation matrix for transforming the flux neutrino coordinates
81  from the default Topocentric Horizontal (see GENIE manual) coordinate
82  system to the user-defined topocentric coordinate system.
83  The rotation is specified by the 3 Euler angles (phi, theta, psi).
84  The Euler angles are input as a comma separated list as:
85  `-R <convention>:phi,theta,psi',
86  where <convention> is either X (for X-convention), Y (for Y-convention),
87  X^-1 or Y^-1 (as previously, but using the inverse matrix).
88  By default, the X-convention (rotation about Z-axis, then about the
89  new X-axis, then about the Z-axis) is used.
90  Notes:
91  - (Extract from TRotation documentation)
92  "Euler angles usually define the rotation of the new coordinate
93  system with respect to the original system, however, the TRotation
94  class specifies the rotation of the object in the original system
95  (an active rotation). To recover the usual Euler rotations (ie.
96  rotate the system not the object), you must take the inverse of
97  the rotation."
98  Examples:
99  1. To set the Euler angles phi=3.14, theta=1.28, psi=1.0 using the
100  X-convention, type: `-R 3.14,1.28,1.0', or `-R X:3.14,1.28,1.0'
101  2. To set the Euler angles phi=3.14, theta=1.28, psi=1.0 using the
102  Y-convention, type: `-R Y:3.14,1.28,1.0'
103  3. To set the Euler angles phi=3.14, theta=1.28, psi=1.0 using the
104  Y-convention, and then use the inverse rotation matrix, type:
105  `-R Y^-1:3.14,1.28,1.0'
106  -L
107  Input geometry length units, eg 'm', 'cm', 'mm', ...
108  [default: 'mm']
109  -D
110  Input geometry density units, eg 'g_cm3', 'clhep_def_density_unit',...
111  [default: 'g_cm3']
112  -t
113  Input 'top volume' for event generation -
114  can be used to force event generation in given sub-detector
115  [default: the 'master volume' of the input geometry]
116  You can also use the -t option to switch generation on/off at
117  multiple volumes as, for example, in:
118  `-t +Vol1-Vol2+Vol3-Vol4',
119  `-t "+Vol1 -Vol2 +Vol3 -Vol4"',
120  `-t -Vol2-Vol4+Vol1+Vol3',
121  `-t "-Vol2 -Vol4 +Vol1 +Vol3"'
122  where:
123  "+Vol1" and "+Vol3" tells GENIE to `switch on' Vol1 and Vol3, while
124  "-Vol2" and "-Vol4" tells GENIE to `switch off' Vol2 and Vol4.
125  If the very first character is a '+', GENIE will neglect all volumes
126  except the ones explicitly turned on. Vice versa, if the very first
127  character is a `-', GENIE will keep all volumes except the ones
128  explicitly turned off (feature contributed by J.Holeczek).
129  -n
130  Specifies how many events to generate.
131  -e
132  Specifies requested exposure in terms of kton*yrs.
133  -E
134  Specifies the neutrino energy in GeV.
135  Must be a comma-separated pair of numbers, eg `-E 0.3,70'
136  [default: 0.5,50]
137  -o
138  Sets the prefix of the output event file.
139  The output filename is built as:
140  [prefix].[run_number].[event_tree_format].[file_format]
141  The default output filename is:
142  gntp.[run_number].ghep.root
143  This cmd line arguments lets you override 'gntp'
144  --seed
145  Random number seed.
146  --cross-sections
147  Name (incl. full path) of an XML file with pre-computed
148  cross-section values used for constructing splines.
149  --tune
150  Specifies a GENIE comprehensive neutrino interaction model tune.
151  [default: "Default"].
152  --message-thresholds
153  Allows users to customize the message stream thresholds.
154  The thresholds are specified using an XML file.
155  See $GENIE/config/Messenger.xml for the XML schema.
156  Multiple files, delimited with a `:' can be specified.
157  --unphysical-event-mask
158  Allows users to specify a 16-bit mask to allow certain types of
159  unphysical events to be written in the output file.
160  [default: all unphysical events are rejected]
161  --event-record-print-level
162  Allows users to set the level of information shown when the event
163  record is printed in the screen. See GHepRecord::Print().
164  --mc-job-status-refresh-rate
165  Allows users to customize the refresh rate of the status file.
166  --cache-file
167  Allows users to specify a cache file so that the cache can be
168  re-used in subsequent MC jobs.
169 
170  *** Examples:
171 
172  (1) Generate 100k events (run number 999210) in the energy range 1-10 GeV
173  for nu_e and nu_mu only, using the sdave_numu07.dat FLUKA flux file for
174  nu_mu and the sdave_nue07.dat file for nu_e (files in /data/flx/).
175  Use the detector geometry in the /data/geo/SuperK.root file, where the
176  geometry length and density units are m and kgr/m^3. Generate events over
177  the entire geometry volume. Pre-computed cross-section data are loaded
178  from /data/xsec.xml.
179 
180  % gevgen_atmo -r 999210 -n 100000 -E 1,10
181  -f FLUKA:/data/flx/sdave_numu07.dat[14],/data/flx/sdave_nue07.dat[12]
182  -g /data/geo/SuperK.root -L "m" -D "kg_m3"
183  --cross-sections /data/xsec.xml
184 
185  (2) Like above but, instead of generating events in a realistic detector
186  geometry, use a simple target mix (88.79% O16 + 11.21% H, i.e. `water')
187 
188  % gevgen_atmo -r 999210 -n 100000 -E 1,10
189  -f /data/flux/sdave_numu07.dat[14],/data/flux/sdave_nue07.dat[12]
190  -g 1000080160[0.8879],1000010010[0.1121]
191  --cross-sections /data/xsec.xml
192 
193  ... to add more
194 
195  Please read the GENIE User Manual for more information.
196 
197 \created August 20, 2010
198 
199 \author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
200  University of Liverpool & STFC Rutherford Appleton Lab
201 
202  Torben Ferber <torben.ferber \at DESY.DE>
203  DESY
204 
205  Hugh Gallagher <hugh.gallagher \at stfc.ac.uk>
206  Tufts University
207 
208  Tarak Thakore <tarak \at mailhost.tifr.res.in>
209  Tata Institute of Fundamental Research
210 
211 \cpright Copyright (c) 2003-2019, The GENIE Collaboration
212  For the full text of the license visit http://copyright.genie-mc.org
213  or see $GENIE/LICENSE
214 */
215 //_________________________________________________________________________________________
216 
217 #include <cassert>
218 #include <cstdlib>
219 #include <cctype>
220 #include <string>
221 #include <vector>
222 #include <sstream>
223 #include <map>
224 
225 #include <TRotation.h>
226 
244 #include "Framework/Utils/AppInit.h"
245 #include "Framework/Utils/RunOpt.h"
246 
247 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
251 #endif
252 
253 #ifdef __GENIE_GEOM_DRIVERS_ENABLED__
254 #include "Tools/Geometry/GeoUtils.h"
257 #endif
258 
259 using std::string;
260 using std::vector;
261 using std::map;
262 using std::ostringstream;
263 
264 using namespace genie;
265 using namespace genie::flux;
266 
267 void GetCommandLineArgs (int argc, char ** argv);
268 void PrintSyntax (void);
269 GFluxI * GetFlux (void);
270 GeomAnalyzerI * GetGeometry (void);
271 
272 // User-specified options:
273 //
274 Long_t gOptRunNu; // run number
275 string gOptFluxSim; // flux simulation (FLUKA, BGLRS or HAKKM)
276 map<int,string> gOptFluxFiles; // neutrino pdg code -> flux file map
277 bool gOptUsingRootGeom = false; // using root geom or target mix?
278 map<int,double> gOptTgtMix; // target mix (tgt pdg -> wght frac) / if not using detailed root geom
279 string gOptRootGeom; // input ROOT file with realistic detector geometry
280 string gOptRootGeomTopVol = ""; // input geometry top event generation volume
281 double gOptGeomLUnits = 0; // input geometry length units
282 double gOptGeomDUnits = 0; // input geometry density units
283 string gOptExtMaxPlXml; // max path lengths XML file for input geometry
284 int gOptNev = -1; // exposure - in terms of number of events
285 double gOptKtonYrExposure = -1; // exposure - in terms of kton*yrs
286 double gOptEvMin; // minimum neutrino energy
287 double gOptEvMax; // maximum neutrino energy
288 string gOptEvFilePrefix; // event file prefix
289 TRotation gOptRot; // coordinate rotation matrix: topocentric horizontal -> user-defined topocentric system
290 long int gOptRanSeed; // random number seed
291 string gOptInpXSecFile; // cross-section splines
292 
293 // Defaults:
294 //
295 NtpMCFormat_t kDefOptNtpFormat = kNFGHEP; // def event tree format
296 string kDefOptEvFilePrefix = "gntp"; // def output prefix (override with -o)
297 string kDefOptGeomLUnits = "mm"; // def geom length units (override with -L)
298 string kDefOptGeomDUnits = "g_cm3"; // def geom density units (override with -D)
299 double kDefOptEvMin = 0.5; // min neutrino energy (override with -E)
300 double kDefOptEvMax = 50.0; // max neutrino energy (override with -E)
301 
302 //________________________________________________________________________________________
303 int main(int argc, char** argv)
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  }
313 
314  // Iinitialization of random number generators, cross-section table, messenger, cache etc...
315  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
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;
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
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 }
377 //________________________________________________________________________________________
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 }
440 //________________________________________________________________________________________
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 }
496 //________________________________________________________________________________________
497 void GetCommandLineArgs(int argc, char ** argv)
498 {
499 // Get the command line arguments
500 
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 }
946 //________________________________________________________________________________________
947 void PrintSyntax(void)
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 }
976 //________________________________________________________________________________________
int iev
Definition: runWimpSim.h:118
void RandGen(long int seed)
Definition: AppInit.cxx:31
static void SetPrintLevel(int print_level)
Definition: GHepRecord.cxx:992
const XML_Char * name
Definition: expat.h:151
void XSecTable(string inpfile, bool require_table)
Definition: AppInit.cxx:39
long ArgAsLong(char opt)
double ArgAsDouble(char opt)
map< int, string > gOptFluxFiles
Definition: gAtmoEvGen.cxx:276
string gOptEvFilePrefix
Definition: gAtmoEvGen.cxx:288
string ArgAsString(char opt)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
void SetEventGeneratorList(string listname)
Definition: GMCJDriver.cxx:157
bool FileExists(string filename)
Definition: SystemUtils.cxx:90
TRotation gOptRot
Definition: gAtmoEvGen.cxx:289
void ReadFromCommandLine(int argc, char **argv)
Definition: RunOpt.cxx:107
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
double gOptGeomDUnits
Definition: gAtmoEvGen.cxx:282
double kDefOptEvMin
Definition: gAtmoEvGen.cxx:299
void Update(int iev, const EventRecord *event)
Definition: GMCJMonitor.cxx:58
string filename
Definition: shutoffs.py:106
void SetRadii(double Rlongitudinal, double Rtransverse)
Definition: GAtmoFlux.cxx:412
map< int, double > gOptTgtMix
Definition: gAtmoEvGen.cxx:278
double gOptEvMin
Definition: gAtmoEvGen.cxx:286
Loaders::FluxType flux
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:303
GFluxI * GetFlux(void)
Definition: gAtmoEvGen.cxx:441
const double emin
double kDefOptEvMax
Definition: gAtmoEvGen.cxx:300
void SetRefreshRate(int rate)
Definition: GMCJMonitor.cxx:53
string kDefOptEvFilePrefix
Definition: gAtmoEvGen.cxx:296
void UseFluxDriver(GFluxI *flux)
Definition: GMCJDriver.cxx:174
GENIE flux drivers.
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
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Definition: typedefs.hpp:11
A flux driver for the FLUKA 3-D Atmospheric Neutrino Flux.
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
void ForceMaxEnergy(double emax)
Definition: GAtmoFlux.cxx:276
string gOptRootGeomTopVol
Definition: gAtmoEvGen.cxx:280
const double emax
void ForceSingleProbScale(void)
Definition: GMCJDriver.cxx:219
double UnitFromString(string u)
Definition: UnitUtils.cxx:26
int gOptNev
Definition: gAtmoEvGen.cxx:284
void Configure(bool calc_prob_scales=true)
Definition: GMCJDriver.cxx:446
A ROOT/GEANT4 geometry driver.
const double j
Definition: BetheBloch.cxx:29
void Save(void)
get the even tree
Definition: NtpWriter.cxx:214
#define pINFO
Definition: Messenger.h:63
void AddEventRecord(int ievent, const EventRecord *ev_rec)
save the event tree
Definition: NtpWriter.cxx:69
const ana::Var wgt
string gOptRootGeom
Definition: gAtmoEvGen.cxx:279
void BuildTune()
build tune and inform XSecSplineList
Definition: RunOpt.cxx:100
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
void SetUserCoordSystem(TRotation &rotation)
Rotation: Topocentric Horizontal -> User-defined Topocentric Coord System.
Definition: GAtmoFlux.cxx:306
void CustomizeFilenamePrefix(string prefix)
Definition: NtpWriter.cxx:123
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void Initialize(void)
add event
Definition: NtpWriter.cxx:95
string kDefOptGeomLUnits
Definition: gAtmoEvGen.cxx:297
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
static RunOpt * Instance(void)
Definition: RunOpt.cxx:62
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:42
A utility class to facilitate creating the GENIE MC Ntuple from the output GENIE GHEP event records...
Definition: NtpWriter.h:40
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
void AddFluxFile(int neutrino_pdg, string filename)
Definition: GAtmoFlux.cxx:421
exit(0)
assert(nhit_max >=nhit_nbins)
A driver for the HAKKM 3-D atmospheric neutrino flux (commonly known as the `Honda flux&#39;) ...
NtpMCFormat_t kDefOptNtpFormat
Definition: gAtmoEvGen.cxx:295
string PrintFramedMesg(string mesg, unsigned int nl=1, const char f='*')
Definition: PrintUtils.cxx:171
A vector of EventGeneratorI objects.
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:100
Command line argument parser.
enum genie::ENtpMCFormat NtpMCFormat_t
#define pNOTICE
Definition: Messenger.h:62
double gOptKtonYrExposure
Definition: gAtmoEvGen.cxx:285
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:497
Defines the GENIE Geometry Analyzer Interface.
Definition: GeomAnalyzerI.h:30
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
string kDefOptGeomDUnits
Definition: gAtmoEvGen.cxx:298
static const double GeV
Definition: Units.h:29
bool gAbortingInErr
Definition: Messenger.cxx:56
A flux driver for the Bartol Atmospheric Neutrino Flux.
void RecursiveExhaust(TGeoVolume *topvol, string volnames, bool exhaust)
Definition: GeoUtils.cxx:26
double gOptGeomLUnits
Definition: gAtmoEvGen.cxx:281
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
void UseGeomAnalyzer(GeomAnalyzerI *geom)
Definition: GMCJDriver.cxx:179
void UseSplines(bool useLogE=true)
Definition: GMCJDriver.cxx:184
The PointGeomAnalyzer class is the simplest implementation of the GeomAnalyserI interface and defines...
void PrintSyntax(void)
Definition: gAtmoEvGen.cxx:947
bool OptionExists(char opt)
was option set?
void CacheFile(string inpfile)
Definition: AppInit.cxx:118
string gOptExtMaxPlXml
Definition: gAtmoEvGen.cxx:283
bool gOptUsingRootGeom
Definition: gAtmoEvGen.cxx:277
virtual TGeoManager * GetGeometry(void) const
#define pDEBUG
Definition: Messenger.h:64
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