gT2KEvGen.cxx
Go to the documentation of this file.
1 //________________________________________________________________________________________
2 /*!
3 
4 \program gevgen_t2k
5 
6 \brief A GENIE event generation driver 'customized' for T2K.
7 
8  This driver can handle the JPARC neutrino flux files generated by JNUBEAM and the
9  realistic detector geometries / target mix of T2K detectors.
10 
11  T2K users should note that the generic GENIE event generation driver (gevgen) may
12  still be a more appropriate tool to use for the simpler event generation cases
13  required for many 4-vector level / systematic studies.
14  Please see the GENIE documentation (http://www.genie-mc.org) and contact me
15  <costas.andreopoulos \at stfc.ac.uk> if in doubt.
16 
17  *** Synopsis :
18 
19  gevgen_t2k [-h]
20  [-r run#]
21  -f flux
22  -g geometry
23  [-p pot_normalization_of_flux_file]
24  [-t top_volume_name_at_geom || -t +Vol1-Vol2...]
25  [-P pre_gen_prob_file_name]
26  [-S] [output_name]
27  [-m max_path_lengths_xml_file]
28  [-L length_units_at_geom]
29  [-D density_units_at_geom]
30  [-n n_of_events]
31  [-c flux_cycles]
32  [-e, -E exposure_in_POTs]
33  [-o output_event_file_prefix]
34  [-R]
35  [--seed random_number_seed]
36  --cross-sections xml_file
37  [--tune genie_tune]
38  [--message-thresholds xml_file]
39  [--unphysical-event-mask mask]
40  [--event-record-print-level level]
41  [--mc-job-status-refresh-rate rate]
42  [--cache-file root_file]
43 
44  *** Options :
45 
46  [] Denotes an optional argument
47 
48  -h
49  Prints out the gevgen_t2k syntax and exits.
50  -r
51  Specifies the MC run number [default: 1000].
52  -g
53  Input 'geometry'.
54  This option can be used to specify any of:
55  1 > A ROOT file containing a ROOT/GEANT geometry description
56  [Note]
57  - This is the standard option for generating events in the
58  nd280, 2km and INGRID detectors.
59  [Examples]
60  - To use the master volume from the ROOT geometry stored
61  in the nd280-geom.root file, type:
62  '-g /some/path/nd280-geom.root'
63  2 > A mix of target materials, each with its corresponding weight,
64  typed as a comma-separated list of nuclear PDG codes (in the
65  std PDG2006 convention: 10LZZZAAAI) with the weight fractions
66  in brackets, eg code1[fraction1],code2[fraction2],...
67  If that option is used (no detailed input geometry description)
68  then the interaction vertices are distributed in the detector
69  by the detector MC.
70  [Note]
71  - This is the standard option for generating events in the
72  SuperK detector.
73  [Examples]
74  - To use a target mix of 95% O16 and 5% H type:
75  '-g 1000080160[0.95],1000010010[0.05]'
76  - To use a target which is 100% C12, type:
77  '-g 1000060120'
78  -t
79  Input 'top volume' for event generation -
80  can be used to force event generation in given sub-detector.
81  [default: the 'master volume' of the input geometry]
82  You can also use the -t option to switch generation on/off at
83  multiple volumes as, for example, in:
84  `-t +Vol1-Vol2+Vol3-Vol4',
85  `-t "+Vol1 -Vol2 +Vol3 -Vol4"',
86  `-t -Vol2-Vol4+Vol1+Vol3',
87  `-t "-Vol2 -Vol4 +Vol1 +Vol3"'m
88  where:
89  "+Vol1" and "+Vol3" tells GENIE to `switch on' Vol1 and Vol3, while
90  "-Vol2" and "-Vol4" tells GENIE to `switch off' Vol2 and Vol4.
91  If the very first character is a '+', GENIE will neglect all volumes
92  except the ones explicitly turned on. Vice versa, if the very first
93  character is a `-', GENIE will keep all volumes except the ones
94  explicitly turned off (feature contributed by J.Holeczek).
95  -P
96  Use exact interaction probabilities (stored in the file specified
97  via this option) corresponding to the flux file input in this job.
98  For complex geometries this will dramatically speed up event generation!
99  If this option is chosen then no max path length file needs to be provided.
100  Instead of calculating the interaction probabilities on the fly
101  per job you can pre-generate them using a dedicated job (see the
102  -S option) and tell the event generator to use them via the -P option.
103  This is advised when processing flux files with more than ~100k entries
104  as the time to pre-calculate the interaction probabilities becomes
105  comparable to the event generation time. For smaller flux files there is
106  less book-keeping if you just calculate them per job and on the fly.
107  -S [output_name]
108  Pre-generate flux interaction probabilities and save to root
109  output file for use with future event generation jobs. With this
110  option no events are generated, it is just a preparatory step
111  before actual event generation. When using this option care
112  should be taken to run with same arguments used for the actual
113  event generation job (same ROOT geomtery, flux file, probe
114  particle types, etc...).
115  The output root file is specified as the input for event
116  generation using the [pre_gen_prob_file] optional value of the
117  -P option. The default output interaction probabilities file
118  name is constructed as: [FLUXFILENAME].[TOPVOL].flxprobs.root.
119  Specifying [output_name] will override this.
120  Introducing multiple functionality to the executable is not
121  desirable but is less error prone than duplicating a lot of the
122  functionality in a separate application.
123  -m
124  An XML file (generated by gmxpl) with the max (density weighted)
125  path-lengths for each target material in the input ROOT geometry.
126  If no file is input, then the geometry will be scanned at MC job
127  initialization to determine those max path lengths.
128  Supplying this file can speed-up the MC job initialization.
129  -L
130  Input geometry length units, eg 'm', 'cm', 'mm', ...
131  [default: 'mm']
132  Note that typically:
133  - nd280m uses: 'mm'
134  - ...
135  -D
136  Input geometry density units, eg 'g_cm3', 'clhep_def_density_unit',...
137  [default: 'g_cm3']
138  Note that typically:
139  - nd280m uses: 'clhep_def_density_unit'
140  - ...
141  -f
142  Input 'neutrino flux'.
143  This option can be used to specify any of:
144  1 > A JNUBEAM beam simulation output file and the detector location.
145  The general sytax is:
146  -f /full/path/flux_file.root,detector_location(,list_of_neutrino_codes)
147  [Notes]
148  - For more information on the flux ntuples, see (T2K internal):
149  http://jnusrv01.kek.jp/internal/t2k/nubeam/flux/
150  - The original HBOOK JNUBAM ntuples need to be converted to a ROOT
151  format using the h2root ROOT utility.
152  - The detector location can be any of 'sk' or the near detector
153  positions 'nd1',...,'nd6','nd13' simulated in JNUBEAM.
154  See the above JNUBEAM web page for more info.
155  - The neutrino codes are the PDG ones.
156  - The (comma separated) list of neutrino codes is optional.
157  It may be used for considering only certain neutrino flavours
158  (eg. nu_e only).
159  If no neutrino list is specified then GENIE will consider all
160  possible flavours (nu_e, nu_e_bar, nu_mu, nu_mu_bar).
161  See examples below.
162  - The JNUBEAM flux ntuples are read via GENIE's GJPARCNuFlux
163  driver. This customized GENIE event generation driver
164  passes-through the complete JPARC input flux information
165  (eg parent decay kinematics / position etc) for each neutrino
166  event it generates (an additional 'flux' branch is added at
167  the output event tree).
168  [Examples]
169  - To use the SuperK flux ntuple from the flux.root file,
170  type:
171  '-f /path/flux.root,sk'
172  - To do the same as above, but considering only nu_mu and
173  nu_mu_bar, type:
174  '-f /path/flux.root,sk,14,-14'
175  - To use the 2km flux ntuple [near detector position 'nd1'
176  in the JNUBEAM flux simulation] from the flux.root file,
177  type:
178  '-f /path/flux.root,nd1'
179  - To use the nd280 flux ntuple [near detector position 'nd5'
180  in the JNUBEAM flux simulation] from the flux.root file,
181  type:
182  '-f /path/flux.root,nd5'
183  - To do the same as above, but considering only nu_e
184  type:
185  '-f /path/flux.root,nd5,12'
186  2 > A list of JNUBEAM beam simulation output files and the detector location.
187  The general sytax is:
188  -f /full/path/flux_file_prefix@first_file_number@last_file_number,detector_location(,list_of_neutrino_codes)
189  [Notes]
190  - The ".root" is assumed.
191  - All the files in the series between flux_file_prefix.first_file_number.root
192  and flux_file_prefix.last_file_number.root should exist.
193  - It is important that you set the -p option correctly. See note below.
194  - Also see notes from option 1.
195  [Examples]
196  - To use the SuperK flux ntuples from the files: flux.0.root, flux.1.root, flux.2.root, flux.3.root
197  type:
198  '-f /path/flux.@0@3,sk'
199  - To use the nd280 flux ntuple [near detector position 'nd5'
200  in the JNUBEAM flux simulation] from the files in the series
201  flux.0.root file to flux.100.root, considering only nu_e, type:
202  '-f /path/flux.@0@100,nd5,12'
203  3 > A set of histograms stored in a ROOT file.
204  The general syntax is:
205  -f /path/histogram_file.root,neutrino_code[histo_name],...
206  [Notes]
207  - The neutrino codes are the PDG ones.
208  - The 'neutrino_code[histogram_name]' part of the option can be
209  repeated multiple times (separated by commas), once for each
210  flux neutrino species you want to consider, eg
211  '-f somefile.root,12[nuehst],-12[nuebarhst],14[numuhst]'
212  - When using flux from histograms then there is no point in using
213  a 'detailed detector geometry description' as your flux input
214  contains no directional information for those flux neutrinos.
215  The neutrino direction is conventionally set to be +z {x=0,y=0}.
216  So, when using this option you must be using a simple 'target mix'
217  See the -g option for possible geometry settings.
218  If you want to use the detailed detector geometry description
219  then you should be feeding this driver with the JNUBEAM flux
220  simulation outputs.
221  - When using flux from histograms there is no branch with neutrino
222  parent information (code, decay mode, 4-momentum at prod & decay)
223  added in the output event tree as your flux input contains no
224  such information. If you want to be getting the neutrino parent
225  information written out as well then you should be feeding this
226  driver with the JNUBEAM flux simulation outputs.
227  - Note that the relative normalization of the flux histograms is
228  taken into account and is reflected in the relative frequency
229  of flux neutrinos thrown by the flux driver
230  [Examples]
231  - To use the histogram 'h100' (representing the nu_mu flux) and
232  the histogram 'h300' (representing the nu_e flux) and the
233  histogram 'h301' (representing the nu_e_bar flux) from the
234  flux.root file in /path/
235  type:
236  '-f /path/flux.root,14[h100],12[h300],-12[h301]
237  -p
238  POT normalization of the input flux file.
239  [default: The 'standard' JNUBEAM flux ntuple normalization of
240  1E+21 POT/detector for the near detectors and
241  1E+21 POT/cm2 for the far detector]
242  That will be used to interpret the flux weights & calculate the actual
243  POT normalization for the generated event sample.
244  [Note]
245  - If you are using the multiple JNUBEAM flux file entry method it is
246  very important that you set this. It should be set to the total POT
247  of all input flux files.
248  [Examples]
249  - If you have 10 standard JNUBEAM files use '-p 10E+21'
250  -c
251  Specifies how many times to cycle a JNUBEAM flux ntuple.
252  Due to the large rejection factor when generating unweighted events
253  in the full energy range (approximately ~500 flux neutrinos will be
254  rejected before getting an interaction in nd280), an option is
255  provided to recycle the flux ntuples for a number of times.
256  That option can be used to boost the generated statistics without
257  requiring enormous flux files.
258  See also 'Note on exposure / statistics' below.
259  -e
260  Specifies how many POTs to generate.
261  If that option is set, gevgen_t2k will work out how many times it has
262  to cycle through the input flux ntuple in order to accumulate the
263  requested statistics. The program will stop at the earliest complete
264  flux ntuple cycle after accumulating the required statistics, so the
265  actual statistics will 'slightly' overshoot that number.
266  See also 'Note on exposure / statistics' below.
267  -E
268  Specifies how many POTs to generate.
269  That option is similar to -e but the program will stop immediatelly
270  after the requested POT has been accumulated. That reduces the
271  generated POT overshoot of the requested POT, but the POT calculation
272  may not be as exact as with -e.
273  See also 'Note on exposure / statistics' below.
274  -n
275  Specifies how many events to generate.
276 
277  --------------------------
278  [Note on setting the exposure / statistics]
279  All -c, -e (-E) and -n options can be used to set the exposure.
280  - If the input flux is a JNUBEAM ntuple then any of these options can
281  be used (only one at a time).
282  If no option is set, then the program will automatically set '-c 1'
283  - If the input flux is described with histograms then only the -n
284  option is available.
285  --------------------------
286 
287  -o
288  Sets the prefix of the output event file.
289  The output filename is built as:
290  [prefix].[run_number].[event_tree_format].[file_format]
291  The default output filename is:
292  gntp.[run_number].ghep.root
293  This cmd line arguments lets you override 'gntp'
294  -R
295  Tell the flux driver to start looping over the flux ntuples with a
296  random offset. May be necessary to avoid biases introduced by always
297  starting at the same point when using very large flux input files.
298  --seed
299  Random number seed.
300  --cross-sections
301  Name (incl. full path) of an XML file with pre-computed
302  cross-section values used for constructing splines.
303  --tune
304  Specifies a GENIE comprehensive neutrino interaction model tune.
305  [default: "Default"].
306  --message-thresholds
307  Allows users to customize the message stream thresholds.
308  The thresholds are specified using an XML file.
309  See $GENIE/config/Messenger.xml for the XML schema.
310  --unphysical-event-mask
311  Allows users to specify a 16-bit mask to allow certain types of
312  unphysical events to be written in the output file.
313  [default: all unphysical events are rejected]
314  --event-record-print-level
315  Allows users to set the level of information shown when the event
316  record is printed in the screen. See GHepRecord::Print().
317  --mc-job-status-refresh-rate
318  Allows users to customize the refresh rate of the status file.
319  --cache-file
320  Allows users to specify a cache file so that the cache can be
321  re-used in subsequent MC jobs.
322 
323  *** Examples:
324 
325  (1) shell% gevgen_t2k
326  -r 1001
327  -f /data/t2k/flux/07a/jnubeam001.root,nd5
328  -g /data/t2k/geom/nd280.root
329  -L mm -D clhep_def_density_unit
330  -e 5E+17
331  --cross-sections /data/t2k/xsec/xsec.xml
332 
333  Generate events (run number 1001) using the JNUBEAM flux ntuple in
334  /data/t2k/flux/07a/jnubeam001.root & picking up the flux entries for
335  the detector location nd5 (:nd280m). The job will load the nd280
336  detector geometry description from /data/t2k/geom/nd280.root and
337  use it thinking that the geometry length unit is 'mm' and the density
338  unit is 'clhep_def_density_unit' (g_cm3 / 0.62415185185E+19)
339  The job will stop on the first complete flux ntuple cycle after
340  generating 5E+17 POT. Pre-computed cross-sections for all relevant
341  initial states are loaded from /data/t2k/xsec/xsec.xml.
342 
343  (2) shell% gevgen_t2k
344  -r 1001
345  -f /data/t2k/flux/07a/jnubeam001.root,nd5
346  -g /data/t2k/geom/nd280.root
347  -L mm -D clhep_def_density_unit
348  -c 100
349  --cross-sections /data/t2k/xsec/xsec.xml
350 
351  As before, but now the job will stop after 100 flux ntuple cycles -
352  whatever POT & number of events that may correspond to.
353 
354  (3) shell% gevgen_t2k
355  -r 1001
356  -f /data/t2k/flux/07a/jnubeam001.root,nd5
357  -g /data/t2k/geom/nd280.root
358  -L mm -D clhep_def_density_unit
359  -n 100000
360  --cross-sections /data/t2k/xsec/xsec.xml
361 
362  As before, but now the job will stop after generating 100000 events -
363  whatever POT & number of flux ntuple cycles that may correspond to.
364 
365  (4) shell% gevgen_t2k
366  -r 1001
367  -f /data/t2k/flux/07a/jnubeam001.root,nd5,-12,12
368  -g /data/t2k/geom/nd280.root
369  -L mm -D clhep_def_density_unit
370  -n 100000
371  --cross-sections /data/t2k/xsec/xsec.xml
372 
373  As before, but now the job will consider flux nu_e and nu_e_bar only!
374 
375  (5) shell% gevgen_t2k
376  -r 1001
377  -f /data/t2k/flux/07a/jnubeam001.root,sk
378  -g 1000080160[0.95],1000010010[0.05]
379  -n 50000
380  --cross-sections /data/t2k/xsec/xsec.xml
381 
382  Generate events (run number 1001) using the JNUBEAM flux ntuple in
383  /data/t2k/flux/07a/jnubeam001.root & picking up the flux entries for
384  the SuperK detector location. This time, the job will not use any
385  detailed detector geometry description but just (95% O16 + 5% H)
386  target mix. The job will stop after generating 50000 events.
387 
388  (6) shell% gevgen_t2k
389  -r 1001
390  -f /data/t2k/flux/hst/flux.root,12[h100],-12[h101],14[h200]
391  -g 1000080160[0.95],1000010010[0.05]
392  -n 50000
393  --cross-sections /data/t2k/xsec/xsec.xml
394 
395  As before, but in this case the flux description is not based on a JNUBEAM
396  ntuple but a set of histograms at the /data/t2k/flux/hst/flux.root file:
397  The histogram named 'h100' will be used for the nu_e flux, 'h101' will
398  will be used for the nu_e_bar flux, and 'h200' for the nu_mu flux.
399 
400 
401  Please read the GENIE User Manual for more information.
402 
403 \author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
404  University of Liverpool & STFC Rutherford Appleton Lab
405 
406 \created February 05, 2008
407 
408 \cpright Copyright (c) 2003-2019, The GENIE Collaboration
409  For the full text of the license visit http://copyright.genie-mc.org
410  or see $GENIE/LICENSE
411 */
412 //_________________________________________________________________________________________
413 
414 #include <cassert>
415 #include <cstdlib>
416 #include <string>
417 #include <sstream>
418 #include <vector>
419 #include <map>
420 
421 #include <TSystem.h>
422 #include <TTree.h>
423 #include <TFile.h>
424 #include <TH1D.h>
425 #include <TMath.h>
426 #include <TGeoVolume.h>
427 #include <TGeoShape.h>
428 #include <TList.h>
429 #include <TObject.h>
430 
446 #include "Framework/Utils/AppInit.h"
447 #include "Framework/Utils/RunOpt.h"
452 
453 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
454 #include "Tools/Flux/GJPARCNuFlux.h"
456 #endif
457 
458 #ifdef __GENIE_GEOM_DRIVERS_ENABLED__
459 #include "Tools/Geometry/GeoUtils.h"
462 #endif
463 
464 using std::string;
465 using std::vector;
466 using std::map;
467 using std::ostringstream;
468 
469 using namespace genie;
470 
471 void GetCommandLineArgs (int argc, char ** argv);
472 void PrintSyntax (void);
473 
474 // Default options (override them using the command line arguments):
475 //
476 string kDefOptGeomLUnits = "mm"; // default geometry length units
477 string kDefOptGeomDUnits = "g_cm3"; // default geometry density units
478 NtpMCFormat_t kDefOptNtpFormat = kNFGHEP; // default event tree format
479 double kDefOptFluxNorm = 1E+21; // std JNUBEAM flux ntuple norm. (POT*detector [NDs] or POT*cm^2 [SK])
480 string kDefOptEvFilePrefix = "gntp";
481 
482 // User-specified options:
483 //
484 Long_t gOptRunNu; // run number
485 bool gOptUsingRootGeom = false; // using root geom or target mix?
486 bool gOptUsingHistFlux = false; // using JNUBEAM flux ntuples or flux from histograms?
487 map<int,double> gOptTgtMix; // target mix (tgt pdg -> wght frac) / if not using detailed root geom
488 map<int,TH1D*> gOptFluxHst; // flux histos (nu pdg -> spectrum) / if not using beam sim ntuples
489 string gOptRootGeom; // input ROOT file with realistic detector geometry
490 string gOptRootGeomTopVol = ""; // input geometry top event generation volume
491 double gOptGeomLUnits = 0; // input geometry length units
492 double gOptGeomDUnits = 0; // input geometry density units
493 string gOptExtMaxPlXml; // max path lengths XML file for input geometry
494 string gOptFluxFile; // ROOT file with JNUBEAM flux ntuple
495 string gOptDetectorLocation; // detector location ('sk','nd1','nd2',...)
496 double gOptFluxNorm; // JNUBEAM flux ntuple normalization
497 PDGCodeList gOptFluxNtpNuList(false); // JNUBEAM flux ntuple neutrinos to consider (can be used for considering certain flavours only)
498 int gOptFluxNCycles; // number of flux ntuple cycles
499 int gOptNev; // number of events to generate
500 double gOptPOT; // exposure (in POT)
501 bool gOptExitAtEndOfFullFluxCycles; // once POT >= requested_POT, stop at once or at the end of the flux cycle?
502 string gOptEvFilePrefix; // event file prefix
503 bool gOptUseFluxProbs = false; // use pre-calculated flux interaction probs instead of estimating them using the max paths
504 bool gOptSaveFluxProbsFile = false; // special mode: no events generated, calculate and save flux interaction probs to root file
505 string gOptFluxProbFileName; // filename for file containg flux probs
506 string gOptSaveFluxProbsFileName; // output filename for pre-generated flux probabilities
507 bool gOptRandomFluxOffset = false; // start looping over flux file from random start entry
508 long int gOptRanSeed; // random number seed
509 string gOptInpXSecFile; // cross-section splines
510 
511 //____________________________________________________________________________
512 int main(int argc, char ** argv)
513 {
514  // Parse command line arguments
515  GetCommandLineArgs(argc,argv);
516 
517 
518  if ( ! RunOpt::Instance()->Tune() ) {
519  LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
520  exit(-1);
521  }
523 
524  // Initialization of random number generators, cross-section table,
525  // messenger thresholds, cache file
526  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
530 
531  // Set GHEP print level
532  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
533 
534  // *************************************************************************
535  // * Create / configure the geometry driver
536  // *************************************************************************
537  GeomAnalyzerI * geom_driver = 0;
538  double zmin=0, zmax=0;
539 
540  if(gOptUsingRootGeom) {
541  //
542  // *** Using a realistic root-based detector geometry description
543  //
544 
545  // creating & configuring a root geometry driver
548  rgeom -> SetLengthUnits (gOptGeomLUnits);
549  rgeom -> SetDensityUnits (gOptGeomDUnits);
550  // getting the bounding box dimensions, before setting topvolume,
551  // along z so as to set the appropriate upstream generation surface
552  // for the JPARC flux driver
553  TGeoVolume * topvol = rgeom->GetGeometry()->GetTopVolume();
554  TGeoShape * bounding_box = topvol->GetShape();
555  bounding_box->GetAxisRange(3, zmin, zmax);
556  zmin *= rgeom->LengthUnits();
557  zmax *= rgeom->LengthUnits();
558  // now update to the requested topvolume for use in recursive exhaust method
559  rgeom -> SetTopVolName (gOptRootGeomTopVol);
560  topvol = rgeom->GetGeometry()->GetTopVolume();
561  if(!topvol) {
562  LOG("gevgen_t2k", pFATAL) << "Null top ROOT geometry volume!";
563  exit(1);
564  }
565  // switch on/off volumes as requested
566  if ( (gOptRootGeomTopVol[0] == '+') || (gOptRootGeomTopVol[0] == '-') ) {
567  bool exhaust = (*gOptRootGeomTopVol.c_str() == '+');
569  }
570 
571  // casting to the GENIE geometry driver interface
572  geom_driver = dynamic_cast<GeomAnalyzerI *> (rgeom);
573  }
574  else {
575  //
576  // *** Using a 'point' geometry with the specified target mix
577  // *** ( = a list of targets with their corresponding weight fraction)
578  //
579 
580  // creating & configuring a point geometry driver
583  // casting to the GENIE geometry driver interface
584  geom_driver = dynamic_cast<GeomAnalyzerI *> (pgeom);
585  }
586 
587  // *************************************************************************
588  // * Create / configure the flux driver
589  // *************************************************************************
590  GFluxI * flux_driver = 0;
591 
592  flux::GJPARCNuFlux * jparc_flux_driver = 0;
593  flux::GCylindTH1Flux * hst_flux_driver = 0;
594 
595  if(!gOptUsingHistFlux) {
596  //
597  // *** Using the detailed JPARC neutrino flux desription by feeding-in
598  // *** the JNUBEAM flux simulation ntuples
599  //
600 
601  // creating JPARC neutrino flux driver
602  jparc_flux_driver = new flux::GJPARCNuFlux;
603  // before loading the beam sim data set whether to use a random offset when looping
604  if(gOptRandomFluxOffset == false) jparc_flux_driver->DisableOffset();
605  // specify input JNUBEAM file & detector location
606  bool beam_sim_data_success = jparc_flux_driver->LoadBeamSimData(gOptFluxFile, gOptDetectorLocation);
607  if(!beam_sim_data_success) {
608  LOG("gevgen_t2k", pFATAL)
609  << "The flux driver has not been properly configured. Exiting";
610  exit(1);
611  }
612  // specify JNUBEAM normalization
613  jparc_flux_driver->SetFilePOT(gOptFluxNorm);
614  // specify upstream generation surface
615  jparc_flux_driver->SetUpstreamZ(zmin);
616  // specify which flavours to consider -
617  // if no neutrino list was specified then the MC job will consider all flavours
618  if( gOptFluxNtpNuList.size() > 0 ) {
619  jparc_flux_driver->SetFluxParticles(gOptFluxNtpNuList);
620  }
621 
622  // casting to the GENIE flux driver interface
623  flux_driver = dynamic_cast<GFluxI *> (jparc_flux_driver);
624  }
625  else {
626  //
627  // *** Using fluxes from histograms (for all specified neutrino species)
628  //
629 
630  // creating & configuring a generic GCylindTH1Flux flux driver
631  TVector3 bdir (0,0,1); // dir along +z
632  TVector3 bspot(0,0,0);
633  hst_flux_driver = new flux::GCylindTH1Flux;
634  hst_flux_driver->SetNuDirection (bdir);
635  hst_flux_driver->SetBeamSpot (bspot);
636  hst_flux_driver->SetTransverseRadius (-1);
637  map<int,TH1D*>::iterator it = gOptFluxHst.begin();
638  for( ; it != gOptFluxHst.end(); ++it) {
639  int pdg_code = it->first;
640  TH1D * spectrum = new TH1D(*(it->second));
641  hst_flux_driver->AddEnergySpectrum(pdg_code, spectrum);
642  }
643  // casting to the GENIE flux driver interface
644  flux_driver = dynamic_cast<GFluxI *> (hst_flux_driver);
645  }
646 
647  // *************************************************************************
648  // * Create/configure the event generation driver
649  // *************************************************************************
650  GMCJDriver * mcj_driver = new GMCJDriver;
652  mcj_driver->UseFluxDriver(flux_driver);
653  mcj_driver->UseGeomAnalyzer(geom_driver);
654  mcj_driver->UseMaxPathLengths(gOptExtMaxPlXml);
655  // do not calculate probability scales if using pre-generated flux probs
656  bool calc_prob_scales = (gOptSaveFluxProbsFile || gOptUseFluxProbs) ? false : true;
657  mcj_driver->Configure(calc_prob_scales);
658  mcj_driver->UseSplines();
659  mcj_driver->ForceSingleProbScale();
660 
661  // *************************************************************************
662  // * If specified use pre-calculated flux interaction probabilities instead
663  // * of preselecting based on max path lengths
664  // *************************************************************************
665 
667 
668  bool success = false;
669 
670  // set flux probs output file name
672  // default output name is ${FLUFILENAME}.${TOPVOL}.flxprobs.root
673  string basename = gOptFluxFile.substr(gOptFluxFile.rfind("/")+1);
674  string name = basename.substr(0, basename.rfind("."));
675  if(gOptRootGeomTopVol.length()>0)
676  name += "."+gOptRootGeomTopVol+".flxprobs.root";
677  else
678  name += ".master.flxprobs.root";
679  // if specified override with cmd line option
681  // Tell the driver save pre-generated probabilities to an output file
682  mcj_driver->SaveFluxProbabilities(name);
683  }
684 
685  // Either load pre-generated flux probabilities
686  if(gOptFluxProbFileName.size() > 0){
687  success = mcj_driver->LoadFluxProbabilities(gOptFluxProbFileName);
688  }
689  // Or pre-calculate them
690  else success = mcj_driver->PreCalcFluxProbabilities();
691 
692  if(success){
693  LOG("gevgen_t2k", pNOTICE)
694  << "Successfully calculated/loaded flux interaction probabilities!";
695  // Print out a list of expected number of events per POT and per cycle
696  // based on the pre-generated flux interaction probabilities
697  map<int, double> sum_probs_map = mcj_driver->SumFluxIntProbs();
698  map<int, double>::const_iterator sum_probs_it = sum_probs_map.begin();
699  double ntot_per_pot = 0.0;
700  double ntot_per_cycle = 0.0;
701  double pscale = mcj_driver->GlobProbScale();
702  double pot_1cycle = jparc_flux_driver->POT_1cycle();
703  LOG("T2KProdInfo", pNOTICE) <<
704  "Expected event rates based on flux interaction probabilities:";
705  for(; sum_probs_it != sum_probs_map.end(); sum_probs_it++){
706  double sum_probs = sum_probs_it->second;
707  double nevts_per_cycle = sum_probs / pscale; // take into account rescale
708  double nevts_per_pot = sum_probs/pot_1cycle; // == (sum_probs*pscale)/(pot_1cycle*pscale)
709  ntot_per_pot += nevts_per_pot;
710  ntot_per_cycle += nevts_per_cycle;
711  LOG("T2KProdInfo", pNOTICE) <<
712  " PDG "<< sum_probs_it->first << ": " << nevts_per_cycle <<
713  " Events/Cycle, "<< nevts_per_pot << " Events/POT";
714  }
715  LOG("T2KProdInfo", pNOTICE) << " -----------";
716  LOG("T2KProdInfo", pNOTICE) <<
717  " All neutrino species: " << ntot_per_cycle <<
718  " Events/Cycle, "<< ntot_per_pot << " Events/POT";
719  LOG("T2KProdInfo", pNOTICE) <<
720  "N.B. This assumes input flux file corresponds to "<< pot_1cycle <<
721  "POT, ensure this is correct if using these numbers!";
722  }
723  else {
724  LOG("gevgen_t2k", pFATAL)
725  << "Failed to calculated/loaded flux interaction probabilities!";
726  return 1;
727  }
728 
729  // Exit now if just pre-generating interaction probabilities
731  LOG("gevgen_t2k", pNOTICE)
732  << "Will not generate events - just pre-calculating flux interaction"
733  << "probabilities";
734  return 0;
735  }
736  } // Pre-calculated flux interaction probabilities
737 
738  // *************************************************************************
739  // * Work out number of cycles for current exposure settings
740  // *************************************************************************
741 
742  if(!gOptUsingHistFlux) {
743  // If a number of POT was requested, then work out how many flux ntuple
744  // cycles are required for accumulating those statistics
745  if(gOptPOT>0) {
746  double fpot_1c = jparc_flux_driver->POT_1cycle(); // flux POT / cycle
747  double psc = mcj_driver->GlobProbScale(); // interaction prob. scale
748  double pot_1c = fpot_1c / psc; // actual POT / cycle
749  int ncycles = (int) TMath::Max(1., TMath::Ceil(gOptPOT/pot_1c));
750 
751  LOG("gevgen_t2k", pNOTICE)
752  << " *** POT/cycle: " << pot_1c;
753  LOG("gevgen_t2k", pNOTICE)
754  << " *** Requested POT will be accumulated in: "
755  << ncycles << " flux ntuple cycles";
756 
757  jparc_flux_driver->SetNumOfCycles(ncycles);
758  }
759  // If a number of events was requested, then set the number of flux
760  // ntuple cycles to 'infinite'
761  else if(gOptNev>0) {
762  jparc_flux_driver->SetNumOfCycles(0);
763  }
764  // Just set the number of cycles to the requested value
765  else {
766  jparc_flux_driver->SetNumOfCycles(gOptFluxNCycles);
767  }
768  }
769 
770  // *************************************************************************
771  // * Prepare for writing the output event tree & status file
772  // *************************************************************************
773 
774  // Initialize an Ntuple Writer to save GHEP records into a TTree
777  ntpw.Initialize();
778 
779  // Add a custom-branch at the standard GENIE event tree so that
780  // info on the flux neutrino parent particle can be passed-through
781  flux::GJPARCNuFluxPassThroughInfo * flux_info = 0;
782  if(!gOptUsingHistFlux) {
783  TBranch * flux = ntpw.EventTree()->Branch("flux",
784  "genie::flux::GJPARCNuFluxPassThroughInfo", &flux_info, 32000, 1);
785  assert(flux);
786  flux->SetAutoDelete(kFALSE);
787  }
788 
789  // Create a MC job monitor for a periodically updated status file
790  GMCJMonitor mcjmonitor(gOptRunNu);
791  mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());
792 
793  // *************************************************************************
794  // * Event generation loop
795  // *************************************************************************
796 
797  int ievent = 0;
798  while (1)
799  {
800  LOG("gevgen_t2k", pNOTICE)
801  << " *** Generating event............ " << ievent;
802 
803  // In case the required statistics was expressed as 'number of events'
804  // then quit if that number has been generated
805  if(ievent == gOptNev) break;
806 
807  // In case the required statistics was expressed as 'number of POT' and
808  // the user does not want to wait till the end of the flux cycle to exit
809  // the event loop, then quit if the requested POT has been generated.
810  // In this case the computed POT may not be as accurate as if the program
811  // was waiting for the current flux cycle to be completed.
813  double fpot = jparc_flux_driver->POT_curravg(); // current POT in flux file
814  double psc = mcj_driver->GlobProbScale(); // interaction prob. scale
815  double pot = fpot / psc; // POT for generated sample
816  if(pot >= gOptPOT) break;
817  }
818 
819  // Generate a single event using neutrinos coming from the specified flux
820  // and hitting the specified geometry or target mix
821  EventRecord * event = mcj_driver->GenerateEvent();
822 
823  // Check whether a null event was returned due to the flux driver reaching
824  // the end of the input flux ntuple - exit the event generation loop
825  if(!event && jparc_flux_driver->End()) {
826  LOG("gevgen_t2k", pNOTICE)
827  << "** The JPARC flux driver read all the input flux ntuple entries";
828  break;
829  }
830  if(!event) {
831  LOG("gevgen_t2k", pERROR)
832  << "Got a null generated neutino event! Retrying ...";
833  continue;
834  }
835  LOG("gevgen_t2k", pINFO)
836  << "Generated event: " << *event;
837 
838  // A valid event was generated: extract flux info (parent decay/prod
839  // position/kinematics) for that simulated event so that it can be
840  // passed-through.
841  // Can only do so if I am generating events using the JNUBEAM flux
842  // ntuples, not simple histograms
843  if(!gOptUsingHistFlux) {
844  flux_info = new flux::GJPARCNuFluxPassThroughInfo(
845  jparc_flux_driver->PassThroughInfo());
846  LOG("gevgen_t2k", pINFO)
847  << "Pass-through flux info associated with generated event: "
848  << *flux_info;
849  }
850 
851  // Add event at the output ntuple, refresh the mc job monitor & clean-up
852  ntpw.AddEventRecord(ievent, event);
853  mcjmonitor.Update(ievent,event);
854  delete event;
855  if(flux_info) delete flux_info;
856  ievent++;
857  } //1
858 
859  LOG("gevgen_t2k", pNOTICE)
860  << "The GENIE MC job is done generaing events - Cleaning up & exiting...";
861 
862  // *************************************************************************
863  // * Print job statistics &
864  // * calculate normalization factor for the generated sample
865  // *************************************************************************
867  {
868  // POT normalization will only be calculated if event generation was based
869  // on beam simulation outputs (not just histograms) & a detailed detector
870  // geometry description.
871  double fpot = jparc_flux_driver->POT_curravg(); // current POT in flux file
872  double psc = mcj_driver->GlobProbScale(); // interaction prob. scale
873  double pot = fpot / psc; // POT for generated sample
874  // Get nunber of flux neutrinos read-in by flux friver, number of flux
875  // neutrinos actually thrown to the event generation driver and number
876  // of neutrino interactions actually generated
877  long int nflx_evg = mcj_driver -> NFluxNeutrinos();
878  long int nflx = jparc_flux_driver -> NFluxNeutrinos();
879  long int nev = ievent;
880 
881  LOG("gevgen_t2k", pNOTICE)
882  << "\n >> Actual JNUBEAM flux file normalization: " << fpot
883  << " POT * " << ((gOptDetectorLocation == "sk") ? "cm^2" : "det")
884  << "\n >> Interaction probability scaling factor: " << psc
885  << "\n >> N of flux v read-in by flux driver: " << nflx
886  << "\n >> N of flux v thrown to event gen driver: " << nflx_evg
887  << "\n >> N of generated v interactions: " << nev
888  << "\n ** Normalization for generated sample: " << pot
889  << " POT * " << ((gOptDetectorLocation == "sk") ? "cm^2" : "det");
890 
891  ntpw.EventTree()->SetWeight(pot); // POT
892  }
893 
894  // *************************************************************************
895  // * MC job meta-data
896  // *************************************************************************
897 
899 
900  metadata -> jnubeam_file = gOptFluxFile;
901  metadata -> detector_location = gOptDetectorLocation;
902  metadata -> geom_file = gOptRootGeom;
903  metadata -> geom_top_volume = gOptRootGeomTopVol;
904  metadata -> geom_length_units = gOptGeomLUnits;
905  metadata -> geom_density_units = gOptGeomDUnits;
906  metadata -> using_root_geom = gOptUsingRootGeom;
907  metadata -> target_mix = gOptTgtMix;
908  metadata -> using_hist_flux = gOptUsingHistFlux;
909  metadata -> flux_hists = gOptFluxHst;
910 
911  ntpw.EventTree()->GetUserInfo()->Add(metadata);
912 
913  // *************************************************************************
914  // * Save & clean-up
915  // *************************************************************************
916 
917  // Save the generated event tree & close the output file
918  ntpw.Save();
919 
920  // Clean-up
921  delete geom_driver;
922  delete flux_driver;
923  delete mcj_driver;
924  map<int,TH1D*>::iterator it = gOptFluxHst.begin();
925  for( ; it != gOptFluxHst.end(); ++it) {
926  TH1D * spectrum = it->second;
927  if(spectrum) delete spectrum;
928  }
929  gOptFluxHst.clear();
930 
931  LOG("gevgen_t2k", pNOTICE) << "Done!";
932 
933  return 0;
934 }
935 //____________________________________________________________________________
936 void GetCommandLineArgs(int argc, char ** argv)
937 {
938  LOG("gevgen_t2k", pINFO) << "Parsing command line arguments";
939 
940  // Common run options. Set defaults and read.
943 
944  // Parse run options for this app
945 
946  CmdLnArgParser parser(argc,argv);
947 
948  // help?
949  bool help = parser.OptionExists('h');
950  if(help) {
951  PrintSyntax();
952  exit(0);
953  }
954 
955  // run number:
956  if( parser.OptionExists('r') ) {
957  LOG("gevgen_t2k", pDEBUG) << "Reading MC run number";
958  gOptRunNu = parser.ArgAsLong('r');
959  } else {
960  LOG("gevgen_t2k", pDEBUG) << "Unspecified run number - Using default";
961  gOptRunNu = 0;
962  } //-r
963 
964  //
965  // *** geometry
966  //
967 
968  string geom = "";
969  string lunits, dunits;
970  if( parser.OptionExists('g') ) {
971  LOG("gevgen_t2k", pDEBUG) << "Getting input geometry";
972  geom = parser.ArgAsString('g');
973 
974  // is it a ROOT file that contains a ROOT geometry?
975  bool accessible_geom_file =
976  ! (gSystem->AccessPathName(geom.c_str()));
977  if (accessible_geom_file) {
978  gOptRootGeom = geom;
979  gOptUsingRootGeom = true;
980  }
981  } else {
982  LOG("gevgen_t2k", pFATAL)
983  << "No geometry option specified - Exiting";
984  PrintSyntax();
985  exit(1);
986  } //-g
987 
988  if(gOptUsingRootGeom) {
989  // using a ROOT geometry - get requested geometry units
990 
991  // legth units:
992  if( parser.OptionExists('L') ) {
993  LOG("gevgen_t2k", pDEBUG)
994  << "Checking for input geometry length units";
995  lunits = parser.ArgAsString('L');
996  } else {
997  LOG("gevgen_t2k", pDEBUG) << "Using default geometry length units";
998  lunits = kDefOptGeomLUnits;
999  } // -L
1000  // density units:
1001  if( parser.OptionExists('D') ) {
1002  LOG("gevgen_t2k", pDEBUG)
1003  << "Checking for input geometry density units";
1004  dunits = parser.ArgAsString('D');
1005  } else {
1006  LOG("gevgen_t2k", pDEBUG) << "Using default geometry density units";
1007  dunits = kDefOptGeomDUnits;
1008  } // -D
1011 
1012  // check whether an event generation volume name has been
1013  // specified -- default is the 'top volume'
1014  if( parser.OptionExists('t') ) {
1015  LOG("gevgen_t2k", pDEBUG) << "Checking for input volume name";
1016  gOptRootGeomTopVol = parser.ArgAsString('t');
1017  } else {
1018  LOG("gevgen_t2k", pDEBUG) << "Using the <master volume>";
1019  } // -t
1020 
1021  // check whether an XML file with the maximum (density weighted)
1022  // path lengths for each detector material is specified -
1023  // otherwise will compute the max path lengths at job init
1024  if( parser.OptionExists('m') ) {
1025  LOG("gevgen_t2k", pDEBUG)
1026  << "Checking for maximum path lengths XML file";
1027  gOptExtMaxPlXml = parser.ArgAsString('m');
1028  } else {
1029  LOG("gevgen_t2k", pDEBUG)
1030  << "Will compute the maximum path lengths at job init";
1031  gOptExtMaxPlXml = "";
1032  } // -m
1033  } // using root geom?
1034 
1035  else {
1036  // User has specified a target mix.
1037  // Decode the list of target pdf codes & their corresponding weight fraction
1038  // (specified as 'pdg_code_1[fraction_1],pdg_code_2[fraction_2],...')
1039  // See documentation on top section of this file.
1040  //
1041  gOptTgtMix.clear();
1042  vector<string> tgtmix = utils::str::Split(geom,",");
1043  if(tgtmix.size()==1) {
1044  int pdg = atoi(tgtmix[0].c_str());
1045  double wgt = 1.0;
1046  gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
1047  } else {
1048  vector<string>::const_iterator tgtmix_iter = tgtmix.begin();
1049  for( ; tgtmix_iter != tgtmix.end(); ++tgtmix_iter) {
1050  string tgt_with_wgt = *tgtmix_iter;
1051  string::size_type open_bracket = tgt_with_wgt.find("[");
1052  string::size_type close_bracket = tgt_with_wgt.find("]");
1053  if (open_bracket ==string::npos ||
1054  close_bracket==string::npos)
1055  {
1056  LOG("gevgen_t2k", pFATAL)
1057  << "You made an error in specifying the target mix";
1058  PrintSyntax();
1059  exit(1);
1060  }
1061  string::size_type ibeg = 0;
1062  string::size_type iend = open_bracket;
1063  string::size_type jbeg = open_bracket+1;
1064  string::size_type jend = close_bracket;
1065  int pdg = atoi(tgt_with_wgt.substr(ibeg,iend-ibeg).c_str());
1066  double wgt = atof(tgt_with_wgt.substr(jbeg,jend-jbeg).c_str());
1067  LOG("gevgen_t2k", pDEBUG)
1068  << "Adding to target mix: pdg = " << pdg << ", wgt = " << wgt;
1069  gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
1070 
1071  }// tgtmix_iter
1072  } // >1 materials in mix
1073  } // using tgt mix?
1074 
1075  //
1076  // *** flux
1077  //
1078 
1079  if( parser.OptionExists('f') ) {
1080  LOG("gevgen_t2k", pDEBUG) << "Getting input flux";
1081  string flux = parser.ArgAsString('f');
1082  gOptUsingHistFlux = (flux.find("[") != string::npos);
1083 
1084  if(!gOptUsingHistFlux) {
1085  // Using JNUBEAM flux ntuples
1086  // Extract JNUBEAM flux (root) file name & detector location.
1087  // Also extract the list of JNUBEAM neutrinos to consider (if none
1088  // is specified here then I will consider all flavours)
1089  //
1090  vector<string> fluxv = utils::str::Split(flux,",");
1091  if(fluxv.size()<2) {
1092  LOG("gevgen_t2k", pFATAL)
1093  << "You need to specify both a flux ntuple ROOT file "
1094  << " _AND_ a detector location";
1095  PrintSyntax();
1096  exit(1);
1097  }
1098  gOptFluxFile = fluxv[0];
1099  gOptDetectorLocation = fluxv[1];
1100 
1101  // Extract the list of neutrinos to consider (if any).
1102  //
1103  for(unsigned int inu = 2; inu < fluxv.size(); inu++)
1104  {
1105  gOptFluxNtpNuList.push_back( atoi(fluxv[inu].c_str()) );
1106  }
1107 
1108  } else {
1109  // Using flux from histograms
1110  // Extract the root file name & the list of histogram names & neutrino
1111  // species (specified as 'filename,histo1[species1],histo2[species2],...')
1112  // See documentation on top section of this file.
1113  //
1114  vector<string> fluxv = utils::str::Split(flux,",");
1115  if(fluxv.size()<2) {
1116  LOG("gevgen_t2k", pFATAL)
1117  << "You need to specify both a flux ntuple ROOT file "
1118  << " _AND_ a detector location";
1119  PrintSyntax();
1120  exit(1);
1121  }
1122  gOptFluxFile = fluxv[0];
1123  bool accessible_flux_file =
1124  !(gSystem->AccessPathName(gOptFluxFile.c_str()));
1125  if (!accessible_flux_file) {
1126  LOG("gevgen_t2k", pFATAL)
1127  << "Can not access flux file: " << gOptFluxFile;
1128  PrintSyntax();
1129  exit(1);
1130  }
1131  // Extract energy spectra for all specified neutrino species
1132  TFile flux_file(gOptFluxFile.c_str(), "read");
1133  for(unsigned int inu=1; inu<fluxv.size(); inu++) {
1134  string nutype_and_histo = fluxv[inu];
1135  string::size_type open_bracket = nutype_and_histo.find("[");
1136  string::size_type close_bracket = nutype_and_histo.find("]");
1137  if (open_bracket ==string::npos ||
1138  close_bracket==string::npos)
1139  {
1140  LOG("gevgen_t2k", pFATAL)
1141  << "You made an error in specifying the flux histograms";
1142  PrintSyntax();
1143  exit(1);
1144  }
1145  string::size_type ibeg = 0;
1146  string::size_type iend = open_bracket;
1147  string::size_type jbeg = open_bracket+1;
1148  string::size_type jend = close_bracket;
1149  string nutype = nutype_and_histo.substr(ibeg,iend-ibeg);
1150  string histo = nutype_and_histo.substr(jbeg,jend-jbeg);
1151  // access specified histogram from the input root file
1152  TH1D * ihst = (TH1D*) flux_file.Get(histo.c_str());
1153  if(!ihst) {
1154  LOG("gevgen_t2k", pFATAL)
1155  << "Can not find histogram: " << histo
1156  << " in flux file: " << gOptFluxFile;
1157  PrintSyntax();
1158  exit(1);
1159  }
1160  // create a local copy of the input histogram
1161  TString origname = ihst->GetName();
1162  TString tmpname; tmpname.Form("%s_", origname.Data());
1163  TH1D * spectrum = new TH1D(
1164  tmpname.Data(), ihst->GetName(), ihst->GetNbinsX(),
1165  ihst->GetXaxis()->GetXmin(), ihst->GetXaxis()->GetXmax());
1166  spectrum->SetDirectory(0);
1167  for(int ibin = 1; ibin <= ihst->GetNbinsX(); ibin++) {
1168  spectrum->SetBinContent(ibin, ihst->GetBinContent(ibin));
1169  }
1170  // get rid of original
1171  delete ihst;
1172  // rename copy
1173  spectrum->SetName(origname.Data());
1174 
1175  // convert neutrino name -> pdg code
1176  int pdg = atoi(nutype.c_str());
1177  if(!pdg::IsNeutrino(pdg) && !pdg::IsAntiNeutrino(pdg)) {
1178  LOG("gevgen_t2k", pFATAL)
1179  << "Unknown neutrino type: " << nutype;
1180  PrintSyntax();
1181  exit(1);
1182  }
1183  // store flux neutrino code / energy spectrum
1184  LOG("gevgen_t2k", pDEBUG)
1185  << "Adding energy spectrum for flux neutrino: pdg = " << pdg;
1186  gOptFluxHst.insert(map<int, TH1D*>::value_type(pdg, spectrum));
1187  }//inu
1188  if(gOptFluxHst.size()<1) {
1189  LOG("gevgen_t2k", pFATAL)
1190  << "You have not specified any flux histogram!";
1191  PrintSyntax();
1192  exit(1);
1193  }
1194  flux_file.Close();
1195  } // flux from histograms or from JNUBEAM ntuples?
1196 
1197  } else {
1198  LOG("gevgen_t2k", pFATAL) << "No flux info was specified - Exiting";
1199  PrintSyntax();
1200  exit(1);
1201  }
1202 
1203  // Use a random offset when looping over flux entries
1204  if(parser.OptionExists('R')) gOptRandomFluxOffset = true;
1205 
1206  //
1207  // *** pre-calculated flux interaction probabilities
1208  //
1209 
1210  // using pre-calculated flux interaction probabilities
1211  if( parser.OptionExists('P') ){
1212  gOptFluxProbFileName = parser.ArgAsString('P');
1213  if(gOptFluxProbFileName.length() > 0){
1214  gOptUseFluxProbs = true;
1215  bool accessible =
1216  !(gSystem->AccessPathName(gOptFluxProbFileName.c_str()));
1217  if(!accessible){
1218  LOG("gevgen_t2k", pFATAL)
1219  << "Can not access pre-calculated flux probabilities file: " << gOptFluxProbFileName;
1220  PrintSyntax();
1221  exit(1);
1222  }
1223  }
1224  else {
1225  LOG("gevgen_t2k", pFATAL)
1226  << "No flux interaction probabilites were specified - exiting";
1227  PrintSyntax();
1228  exit(1);
1229  }
1230  }
1231 
1232  // pre-generating interaction probs and saving to output file
1233  if( parser.OptionExists('S') ){
1234  gOptSaveFluxProbsFile = true;
1235  gOptSaveFluxProbsFileName = parser.ArgAsString('S');
1236  }
1237 
1238  // cannot save and run at the same time
1240  LOG("gevgen_t2k", pFATAL)
1241  << "Cannot specify both the -P and -S options at the same time!";
1242  exit(1);
1243  }
1244 
1245  // only makes sense to be setting these options for a realistic flux
1247  LOG("gevgen_t2k", pFATAL)
1248  << "Using pre-calculated flux interaction probabilities only makes "
1249  << "sense when using JNUBEAM flux option!";
1250  exit(1);
1251  }
1252 
1253  // flux file POT normalization
1254  // only relevant when using the JNUBEAM flux ntuples
1255  if( parser.OptionExists('p') ) {
1256  LOG("gevgen_t2k", pDEBUG) << "Reading flux file normalization";
1257  gOptFluxNorm = parser.ArgAsDouble('p');
1258  } else {
1259  LOG("gevgen_t2k", pDEBUG)
1260  << "Setting standard normalization for JNUBEAM flux ntuples";
1262  } //-p
1263 
1264  // number of times to cycle through the JNUBEAM flux ntuple contents
1265  if( parser.OptionExists('c') ) {
1266  LOG("gevgen_t2k", pDEBUG) << "Reading number of flux ntuple cycles";
1267  gOptFluxNCycles = parser.ArgAsInt('c');
1268  } else {
1269  LOG("gevgen_t2k", pDEBUG)
1270  << "Setting standard number of cycles for JNUBEAM flux ntuples";
1271  gOptFluxNCycles = -1;
1272  } //-c
1273 
1274  // limit on max number of events that can be generated
1275  if( parser.OptionExists('n') ) {
1276  LOG("gevgen_t2k", pDEBUG)
1277  << "Reading limit on number of events to generate";
1278  gOptNev = parser.ArgAsInt('n');
1279  } else {
1280  LOG("gevgen_t2k", pDEBUG)
1281  << "Will keep on generating events till the flux driver stops";
1282  gOptNev = -1;
1283  } //-n
1284 
1285  // exposure (in POT)
1286  bool uppc_e = parser.OptionExists('E');
1287  char pot_args = 'e';
1288  bool pot_exit = true;
1289  if(uppc_e) {
1290  pot_args = 'E';
1291  pot_exit = false;
1292  }
1293  gOptExitAtEndOfFullFluxCycles = pot_exit;
1294  if( parser.OptionExists(pot_args) ) {
1295  LOG("gevgen_t2k", pDEBUG) << "Reading requested exposure in POT";
1296  gOptPOT = parser.ArgAsDouble(pot_args);
1297  } else {
1298  LOG("gevgen_t2k", pDEBUG) << "No POT exposure was requested";
1299  gOptPOT = -1;
1300  } //-e, -E
1301 
1302  // event file prefix
1303  if( parser.OptionExists('o') ) {
1304  LOG("gevgen_t2k", pDEBUG) << "Reading the event filename prefix";
1305  gOptEvFilePrefix = parser.ArgAsString('o');
1306  } else {
1307  LOG("gevgen_t2k", pDEBUG)
1308  << "Will set the default event filename prefix";
1310  } //-o
1311 
1312 
1313  // random number seed
1314  if( parser.OptionExists("seed") ) {
1315  LOG("gevgen_t2k", pINFO) << "Reading random number seed";
1316  gOptRanSeed = parser.ArgAsLong("seed");
1317  } else {
1318  LOG("gevgen_t2k", pINFO) << "Unspecified random number seed - Using default";
1319  gOptRanSeed = -1;
1320  }
1321 
1322  // input cross-section file
1323  if( parser.OptionExists("cross-sections") ) {
1324  LOG("gevgen_t2k", pINFO) << "Reading cross-section file";
1325  gOptInpXSecFile = parser.ArgAsString("cross-sections");
1326  } else {
1327  LOG("gevgen_t2k", pINFO) << "Unspecified cross-section file";
1328  gOptInpXSecFile = "";
1329  }
1330 
1331  //
1332  // >>> perform 'sanity' checks on command line arguments
1333  //
1334 
1335  // If we use a JNUBEAM flux ntuple, the 'exposure' may be set
1336  // either as:
1337  // - a number of POTs (whichever number of flux ntuple cycles that corresponds to)
1338  // - a number of generated events (whichever number of POTs that corresponds to)
1339  // - a number of flux ntuple cycles (whichever number of POTs that corresponds to)
1340  // Only one of those options can be set.
1341  if(!gOptUsingHistFlux) {
1342  int nset=0;
1343  if(gOptPOT > 0) nset++;
1344  if(gOptFluxNCycles > 0) nset++;
1345  if(gOptNev > 0) nset++;
1346  if(nset==0) {
1347  LOG("gevgen_t2k", pWARN)
1348  << "** To use a JNUBEAM flux ntuple you need to specify an exposure, "
1349  << "either via the -c, -e or -n options";
1350  LOG("gevgen_t2k", pWARN)
1351  << "** gevgen_t2k automatically sets the exposure via '-c 1'";
1352  gOptFluxNCycles = 1;
1353  }
1354  if(nset>1) {
1355  LOG("gevgen_t2k", pFATAL)
1356  << "You can not specify more than one of the -c, -e or -n options";
1357  PrintSyntax();
1358  exit(1);
1359  }
1360  }
1361  // If we use a flux histograms (not JNUBEAM flux ntuples) then -currently- the
1362  // only way to control exposure is via a number of events
1363  if(gOptUsingHistFlux) {
1364  if(gOptNev < 0) {
1365  LOG("gevgen_t2k", pFATAL)
1366  << "If you're using flux from histograms you need to specify the -n option";
1367  PrintSyntax();
1368  exit(1);
1369  }
1370  }
1371  // If we don't use a detailed ROOT detector geometry (but just a target mix) then
1372  // don't accept POT as a way to control job statistics (not enough info is passed
1373  // in the target mix to compute POT & the calculation can be easily done offline)
1374  if(!gOptUsingRootGeom) {
1375  if(gOptPOT > 0) {
1376  LOG("gevgen_t2k", pFATAL)
1377  << "You may not use the -e, -E options "
1378  << "without a detailed detector geometry description input";
1379  exit(1);
1380  }
1381  }
1382 
1383  //
1384  // >>> print the command line options
1385  //
1386  PDGLibrary * pdglib = PDGLibrary::Instance();
1387 
1388  ostringstream gminfo;
1389  if (gOptUsingRootGeom) {
1390  gminfo << "Using ROOT geometry - file: " << gOptRootGeom
1391  << ", top volume: "
1392  << ((gOptRootGeomTopVol.size()==0) ? "<master volume>" : gOptRootGeomTopVol)
1393  << ", max{PL} file: "
1394  << ((gOptExtMaxPlXml.size()==0) ? "<none>" : gOptExtMaxPlXml)
1395  << ", length units: " << lunits
1396  << ", density units: " << dunits;
1397  } else {
1398  gminfo << "Using target mix - ";
1399  map<int,double>::const_iterator iter;
1400  for(iter = gOptTgtMix.begin(); iter != gOptTgtMix.end(); ++iter) {
1401  int pdg_code = iter->first;
1402  double wgt = iter->second;
1403  TParticlePDG * p = pdglib->Find(pdg_code);
1404  if(p) {
1405  string name = p->GetName();
1406  gminfo << "(" << name << ") -> " << 100*wgt << "% / ";
1407  }//p?
1408  }
1409  }
1410 
1411  ostringstream fluxinfo;
1412  if(gOptUsingHistFlux) {
1413  fluxinfo << "Using flux histograms - ";
1414  map<int,TH1D*>::const_iterator iter;
1415  for(iter = gOptFluxHst.begin(); iter != gOptFluxHst.end(); ++iter) {
1416  int pdg_code = iter->first;
1417  TH1D * spectrum = iter->second;
1418  TParticlePDG * p = pdglib->Find(pdg_code);
1419  if(p) {
1420  string name = p->GetName();
1421  fluxinfo << "(" << name << ") -> " << spectrum->GetName() << " / ";
1422  }//p?
1423  }
1424  } else {
1425  fluxinfo << "Using JNUBEAM flux ntuple - "
1426  << "file: " << gOptFluxFile
1427  << ", location: " << gOptDetectorLocation
1428  << ", pot norm: " << gOptFluxNorm;
1429 
1430  if( gOptFluxNtpNuList.size() > 0 ) {
1431  fluxinfo << ", ** this job is considering only: ";
1432  for(unsigned int inu = 0; inu < gOptFluxNtpNuList.size(); inu++) {
1433  fluxinfo << gOptFluxNtpNuList[inu];
1434  if(inu < gOptFluxNtpNuList.size()-1) fluxinfo << ",";
1435  }
1436  }
1437 
1438  }
1439 
1440  ostringstream exposure;
1441  if(gOptPOT > 0)
1442  exposure << "Number of POTs = " << gOptPOT;
1443  if(gOptFluxNCycles > 0)
1444  exposure << "Number of flux cycles = " << gOptFluxNCycles;
1445  if(gOptNev > 0)
1446  exposure << "Number of events = " << gOptNev;
1447 
1448  LOG("gevgen_t2k", pNOTICE)
1449  << "\n\n"
1450  << utils::print::PrintFramedMesg("T2K event generation job configuration");
1451 
1452  LOG("gevgen_t2k", pNOTICE)
1453  << "\n - Run number: " << gOptRunNu
1454  << "\n - Random number seed: " << gOptRanSeed
1455  << "\n - Using cross-section file: " << gOptInpXSecFile
1456  << "\n - Flux @ " << fluxinfo.str()
1457  << "\n - Geometry @ " << gminfo.str()
1458  << "\n - Exposure @ " << exposure.str();
1459 
1460  LOG("gevgen_t2k", pNOTICE) << *RunOpt::Instance();
1461 }
1462 //____________________________________________________________________________
1463 void PrintSyntax(void)
1464 {
1465  LOG("gevgen_t2k", pFATAL)
1466  << "\n **Syntax**"
1467  << "\n gevgen_t2k [-h] "
1468  << "\n [-r run#]"
1469  << "\n -f flux"
1470  << "\n -g geometry"
1471  << "\n [-p pot_normalization_of_flux_file]"
1472  << "\n [-t top_volume_name_at_geom]"
1473  << "\n [-P pre_gen_prob_file]"
1474  << "\n [-S] [output_name]"
1475  << "\n [-m max_path_lengths_xml_file]"
1476  << "\n [-L length_units_at_geom]"
1477  << "\n [-D density_units_at_geom]"
1478  << "\n [-n n_of_events]"
1479  << "\n [-c flux_cycles]"
1480  << "\n [-e, -E exposure_in_POTs]"
1481  << "\n [-o output_event_file_prefix]"
1482  << "\n [-R]"
1483  << "\n [--seed random_number_seed]"
1484  << "\n --cross-sections xml_file"
1485  << "\n [--event-generator-list list_name]"
1486  << "\n [--message-thresholds xml_file]"
1487  << "\n [--unphysical-event-mask mask]"
1488  << "\n [--event-record-print-level level]"
1489  << "\n [--mc-job-status-refresh-rate rate]"
1490  << "\n [--cache-file root_file]"
1491  << "\n"
1492  << " Please also read the detailed documentation at http://www.genie-mc.org"
1493  << " or look at the source code: $GENIE/src/Apps/gT2KEvGen.cxx"
1494  << "\n";
1495 }
1496 //____________________________________________________________________________
map< int, TH1D * > gOptFluxHst
Definition: gT2KEvGen.cxx:488
string gOptRootGeomTopVol
Definition: gT2KEvGen.cxx:490
string gOptRootGeom
Definition: gT2KEvGen.cxx:489
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
double gOptFluxNorm
Definition: gT2KEvGen.cxx:496
double POT_curravg(void)
current average POT
void XSecTable(string inpfile, bool require_table)
Definition: AppInit.cxx:39
long ArgAsLong(char opt)
double ArgAsDouble(char opt)
bool PreCalcFluxProbabilities(void)
Definition: GMCJDriver.cxx:239
map< int, double > SumFluxIntProbs(void) const
Definition: GMCJDriver.h:74
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
string ArgAsString(char opt)
set< int >::iterator it
bool LoadBeamSimData(string filename, string det_loc)
load a jnubeam root flux ntuple
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
bool End(void)
true if no more flux nu&#39;s can be thrown (eg reaching end of beam sim ntuples)
Definition: GJPARCNuFlux.h:67
#define pERROR
Definition: Messenger.h:60
void SetUpstreamZ(double z0)
set flux neutrino initial z position (upstream of the detector)
void SetEventGeneratorList(string listname)
Definition: GMCJDriver.cxx:157
string kDefOptGeomDUnits
Definition: gT2KEvGen.cxx:477
double POT_1cycle(void)
flux POT per cycle
void ReadFromCommandLine(int argc, char **argv)
Definition: RunOpt.cxx:107
map< int, double > gOptTgtMix
Definition: gT2KEvGen.cxx:487
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
bool gOptUsingHistFlux
Definition: gT2KEvGen.cxx:486
string gOptDetectorLocation
Definition: gT2KEvGen.cxx:495
NtpMCFormat_t kDefOptNtpFormat
Definition: gT2KEvGen.cxx:478
void Update(int iev, const EventRecord *event)
Definition: GMCJMonitor.cxx:58
double gOptPOT
Definition: gT2KEvGen.cxx:500
void SetFluxParticles(const PDGCodeList &particles)
specify list of flux neutrino species
Long_t gOptRunNu
Definition: gT2KEvGen.cxx:484
Loaders::FluxType flux
string gOptFluxProbFileName
Definition: gT2KEvGen.cxx:505
string gOptFluxFile
Definition: gT2KEvGen.cxx:494
A list of PDG codes.
Definition: PDGCodeList.h:33
int gOptNev
Definition: gT2KEvGen.cxx:499
string gOptInpXSecFile
Definition: gT2KEvGen.cxx:509
A GENIE flux driver encapsulating the JPARC neutrino flux. It reads-in the official JPARC neutrino fl...
Definition: GJPARCNuFlux.h:51
void SetRefreshRate(int rate)
Definition: GMCJMonitor.cxx:53
void UseFluxDriver(GFluxI *flux)
Definition: GMCJDriver.cxx:174
bool gOptRandomFluxOffset
Definition: gT2KEvGen.cxx:507
unsigned int nutype(unsigned int u)
Definition: runWimpSim.h:122
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
double GlobProbScale(void) const
Definition: GMCJDriver.h:72
string kDefOptGeomLUnits
Definition: gT2KEvGen.cxx:476
string gOptExtMaxPlXml
Definition: gT2KEvGen.cxx:493
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Definition: typedefs.hpp:11
void PrintSyntax(void)
Definition: gT2KEvGen.cxx:1463
TTree * EventTree(void)
Definition: NtpWriter.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Float_t E
Definition: plot.C:20
void ForceSingleProbScale(void)
Definition: GMCJDriver.cxx:219
void SetFilePOT(double pot)
flux file norm is in /N POT/det [ND] or /N POT/cm^2 [FD]. Specify N (typically 1E+21) ...
virtual double LengthUnits(void) const
A generic GENIE flux driver. Generates a &#39;cylindrical&#39; neutrino beam along the input direction...
bool gOptUseFluxProbs
Definition: gT2KEvGen.cxx:503
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
bool UseMaxPathLengths(string xml_filename)
Definition: GMCJDriver.cxx:190
double UnitFromString(string u)
Definition: UnitUtils.cxx:26
#define pot
int main(int argc, char **argv)
Definition: gT2KEvGen.cxx:512
void Configure(bool calc_prob_scales=true)
Definition: GMCJDriver.cxx:446
string kDefOptEvFilePrefix
Definition: gT2KEvGen.cxx:480
A ROOT/GEANT4 geometry driver.
double kDefOptFluxNorm
Definition: gT2KEvGen.cxx:479
def success(message)
Definition: log.py:5
void Save(void)
get the even tree
Definition: NtpWriter.cxx:214
TH2D * histo
#define pINFO
Definition: Messenger.h:63
void GetCommandLineArgs(int argc, char **argv)
Definition: gT2KEvGen.cxx:936
void AddEventRecord(int ievent, const EventRecord *ev_rec)
save the event tree
Definition: NtpWriter.cxx:69
void SetTransverseRadius(double Rt)
void DisableOffset(void)
switch off random offset, must be called before LoadBeamSimData to have any effect ...
Definition: GJPARCNuFlux.h:84
const ana::Var wgt
void SetNuDirection(const TVector3 &direction)
void BuildTune()
build tune and inform XSecSplineList
Definition: RunOpt.cxx:100
bool gOptUsingRootGeom
Definition: gT2KEvGen.cxx:485
#define pWARN
Definition: Messenger.h:61
PDGCodeList gOptFluxNtpNuList(false)
string gOptSaveFluxProbsFileName
Definition: gT2KEvGen.cxx:506
EventRecord * GenerateEvent(void)
Definition: GMCJDriver.cxx:838
long int gOptRanSeed
Definition: gT2KEvGen.cxx:508
double gOptGeomDUnits
Definition: gT2KEvGen.cxx:492
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
void CustomizeFilenamePrefix(string prefix)
Definition: NtpWriter.cxx:123
bool gOptSaveFluxProbsFile
Definition: gT2KEvGen.cxx:504
void Initialize(void)
add event
Definition: NtpWriter.cxx:95
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
string gOptEvFilePrefix
Definition: gT2KEvGen.cxx:502
int gOptFluxNCycles
Definition: gT2KEvGen.cxx:498
double gOptGeomLUnits
Definition: gT2KEvGen.cxx:491
void geom(int which=0)
Definition: geom.C:163
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
exit(0)
bool LoadFluxProbabilities(string filename)
Definition: GMCJDriver.cxx:390
void SetNumOfCycles(int n)
set how many times to cycle through the ntuple (default: 1 / n=0 means &#39;infinite&#39;) ...
assert(nhit_max >=nhit_nbins)
string PrintFramedMesg(string mesg, unsigned int nl=1, const char f='*')
Definition: PrintUtils.cxx:171
A vector of EventGeneratorI objects.
void SetBeamSpot(const TVector3 &spot)
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:100
void AddEnergySpectrum(int nu_pdgc, TH1D *spectrum)
Utility class to store MC job meta-data.
Command line argument parser.
enum genie::ENtpMCFormat NtpMCFormat_t
#define pNOTICE
Definition: Messenger.h:62
bool gOptExitAtEndOfFullFluxCycles
Definition: gT2KEvGen.cxx:501
Defines the GENIE Geometry Analyzer Interface.
Definition: GeomAnalyzerI.h:30
void RecursiveExhaust(TGeoVolume *topvol, string volnames, bool exhaust)
Definition: GeoUtils.cxx:26
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...
bool OptionExists(char opt)
was option set?
void CacheFile(string inpfile)
Definition: AppInit.cxx:118
void SaveFluxProbabilities(string outfilename)
Definition: GMCJDriver.cxx:437
void EnableBareXSecPreCalc(bool flag)
Definition: RunOpt.h:61
const GJPARCNuFluxPassThroughInfo & PassThroughInfo(void)
GJPARCNuFluxPassThroughInfo.
Definition: GJPARCNuFlux.h:93
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
virtual TGeoManager * GetGeometry(void) const
#define pDEBUG
Definition: Messenger.h:64
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
enum BeamMode string