GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gMakeSplinesDM.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gmkspl_dm
5 
6 \brief GENIE utility program building XML cross section splines that can
7  be loaded into GENIE to speed-up event generation.
8  The list of neutrino PDG codes is passed from the command line.
9  The list of nuclear target PDG codes is either passed from the
10  command line or extracted from the input ROOT/GEANT geometry.
11 
12  Syntax :
13  gmkspl_dm -m masses
14  <-t target_pdg_codes,
15  -f geometry_file>
16  <-o | --output-cross-sections> output_xml_xsec_file
17  [-g zp_couplings]
18  [-z med_ratios]
19  [-n nknots]
20  [-e max_energy]
21  [--no-copy]
22  [--seed random_number_seed]
23  [--input-cross-sections xml_file]
24  [--event-generator-list list_name]
25  [--tune genie_tune]
26  [--message-thresholds xml_file]
27 
28  Note :
29  [] marks optional arguments.
30  <> marks a list of arguments out of which only one can be
31  selected at any given time.
32 
33  Options :
34  -m
35  A comma separated list of DM masses.
36  -t
37  A comma separated list of tgt PDG codes.
38  PDG code format: 10LZZZAAAI
39  -f
40  A ROOT file containing a ROOT/GEANT geometry description.
41  -o, --output-cross-sections
42  Name of output XML file containing computed cross-section data.
43  Default: `xsec_splines.xml'.
44  -g
45  A comma separated list of Z' coupling constants
46  Default: Value in UserPhysicsOptions.xml
47  -z
48  A comma separated list of ratios of mediator mass to dark
49  matter mass
50  Default: 0.5
51  -n
52  Number of knots per spline.
53  Default: 15 knots per decade of energy range with a minimum
54  of 30 knots totally.
55  -e
56  Maximum energy in spline.
57  Default: The max energy in the validity range of the spline
58  generating thread.
59  --no-copy
60  Does not write out the input cross-sections in the output file
61  --seed
62  Random number seed.
63  --input-cross-sections
64  Name (incl. full path) of an XML file with pre-computed
65  free-nucleon cross-section values. If loaded, it can speed-up
66  cross-section calculation for nuclear targets.
67  --event-generator-list
68  List of event generators to load in event generation drivers.
69  [default: "Default"].
70  --tune
71  Specifies a GENIE comprehensive neutrino interaction model tune.
72  [default: "Default"].
73  --message-thresholds
74  Allows users to customize the message stream thresholds.
75  The thresholds are specified using an XML file.
76  See $GENIE/config/Messenger.xml for the XML schema.
77 
78  *** See the User Manual for more details and examples. ***
79 
80 \author Joshua Berger <jberger \at physics.wisc.edu>
81  University of Wisconsin-Madison
82  Costas Andreopoulos <c.andreopoulos \at cern.ch>
83  University of Liverpool
84 
85 \created September 1, 2017
86 
87 \cpright Copyright (c) 2003-2024, The GENIE Collaboration
88  For the full text of the license visit http://copyright.genie-mc.org
89 
90 */
91 //____________________________________________________________________________
92 
93 #include <cassert>
94 #include <cstdlib>
95 #include <string>
96 #include <vector>
97 
98 #if defined(HAVE_FENV_H) && defined(HAVE_FEENABLEEXCEPT)
99 #include <fenv.h> // for `feenableexcept`
100 #endif
101 
102 #include <TSystem.h>
103 
105 #include "Framework/Conventions/GBuild.h"
113 #include "Framework/Utils/RunOpt.h"
114 #include "Framework/Utils/AppInit.h"
116 //#include "Framework/Utils/SystemUtils.h"
120 
121 #ifdef __GENIE_GEOM_DRIVERS_ENABLED__
123 #endif
124 
125 using std::string;
126 using std::vector;
127 
128 using namespace genie;
129 
130 #ifdef __GENIE_GEOM_DRIVERS_ENABLED__
131 using namespace genie::geometry;
132 #endif
133 
134 // Prototypes:
135 void GetCommandLineArgs (int argc, char ** argv);
136 void PrintSyntax (void);
137 PDGCodeList * GetTargetCodes (void);
138 
139 // User-specified options:
140 string gOptTgtPdgCodeList = "";
141 string gOptGeomFilename = "";
142 int gOptNKnots = -1;
143 double gOptMaxE = -1.;
144 vector<double> gOptDMMasses;
145 vector<double> gOptMedRatios;
146 vector<double> gOptZpCouplings;
147 bool gOptNoCopy = false;
148 long int gOptRanSeed = -1; // random number seed
149 string gOptInpXSecFile = ""; // input cross-section file
150 string gOptOutXSecFile = ""; // output cross-section file
151 
152 //____________________________________________________________________________
153 int main(int argc, char ** argv)
154 {
155  // Parse command line arguments
156  GetCommandLineArgs(argc,argv);
157 
158  if ( ! RunOpt::Instance()->Tune() ) {
159  LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
160  exit(-1);
161  }
163 
164  for (vector<double>::iterator mass = gOptDMMasses.begin(); mass != gOptDMMasses.end(); ++mass) {
165  for (vector<double>::iterator ratio = gOptMedRatios.begin(); ratio != gOptMedRatios.end(); ++ratio) {
166  for (vector<double>::iterator coup = gOptZpCouplings.begin(); coup != gOptZpCouplings.end(); ++coup) {
167  // Add dark matter to the table
169  PDGLibrary::Instance()->AddDarkMatter(*mass,*ratio);
170  if (*coup > 0.) {
171  Registry * r = AlgConfigPool::Instance()->CommonList("Param", "BoostedDarkMatter");
172  r->UnLock();
173  r->Set("ZpCoupling", *coup);
174  r->Lock();
175  }
176 
177 
178  // throw on NaNs and Infs...
179 #if defined(HAVE_FENV_H) && defined(HAVE_FEENABLEEXCEPT)
180  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
181 #endif
182 
183  // Init
184  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
187 
188  // Get list of neutrinos and nuclear targets
189 
190  // PDGCodeList * neutrinos = GetNeutrinoCodes();
191  PDGCodeList * targets = GetTargetCodes();
192 
193  if(!targets || targets->size() == 0 ) {
194  LOG("gmkspl_dm", pFATAL) << "Empty target PDG code list";
195  PrintSyntax();
196  exit(3);
197  }
198 
199  LOG("gmkspl_dm", pINFO) << "Targets: " << *targets;
200 
201  // Loop over all possible input init states and ask the GEVGDriver
202  // to build splines for all the interactions that its loaded list
203  // of event generators can generate.
204 
205  // PDGCodeList::const_iterator nuiter;
206  PDGCodeList::const_iterator tgtiter;
207  // for(nuiter = neutrinos->begin(); nuiter != neutrinos->end(); ++nuiter) {
208  for(tgtiter = targets->begin(); tgtiter != targets->end(); ++tgtiter) {
209  int dmpdgc = kPdgDarkMatter;
210  int tgtpdgc = *tgtiter;
211  InitialState init_state(tgtpdgc, dmpdgc);
212  GEVGDriver driver;
214  driver.Configure(init_state);
216  }
217  delete targets;
218  }
219  }
220  }
221  // }
222 
223  // Save the splines at the requested XML file
225  bool save_init = !gOptNoCopy;
226  xspl->SaveAsXml(gOptOutXSecFile, save_init);
227 
228 
229  return 0;
230 }
231 //____________________________________________________________________________
232 void GetCommandLineArgs(int argc, char ** argv)
233 {
234  LOG("gmkspl_dm", pINFO) << "Parsing command line arguments";
235 
236  // Common run options. Set defaults and read.
239 
240  // Parse run options for this app
241 
242  CmdLnArgParser parser(argc,argv);
243 
244  // output XML file name
245  if( parser.OptionExists('o') ||
246  parser.OptionExists("output-cross-sections") )
247  {
248  LOG("gmkspl_dm", pINFO) << "Reading output filename";
249  if( parser.OptionExists('o') ) {
250  gOptOutXSecFile = parser.ArgAsString('o');
251  }
252  else {
253  gOptOutXSecFile = parser.ArgAsString("output-cross-sections");
254  }
255  } else {
256  LOG("gmkspl_dm", pINFO) << "Unspecified filename - Using default";
257  gOptOutXSecFile = "xsec_splines.xml";
258  }
259 
260  // number of knots
261  if( parser.OptionExists('n') ) {
262  LOG("gmkspl_dm", pINFO) << "Reading number of knots/spline";
263  gOptNKnots = parser.ArgAsInt('n');
264  } else {
265  LOG("gmkspl_dm", pINFO)
266  << "Unspecified number of knots - Using default";
267  gOptNKnots = -1;
268  }
269 
270  // max spline energy (if < max of validity range)
271  if( parser.OptionExists('e') ) {
272  LOG("gmkspl_dm", pINFO) << "Reading maximum spline energy";
273  gOptMaxE = parser.ArgAsDouble('e');
274  } else {
275  LOG("gmkspl_dm", pINFO)
276  << "Unspecified maximum spline energy - Using default";
277  gOptMaxE = -1;
278  }
279 
280  // get the dark matter mass
281  if( parser.OptionExists('m') ) {
282  LOG("gmkspl_dm", pINFO) << "Reading dark matter mass";
283  string dmm = parser.ArgAsString('m');
284 
285  if(dmm.find(",") != string::npos) {
286  vector<string> massrange = utils::str::Split(dmm, ",");
287  assert(massrange.size() > 1);
288  for (vector<string>::iterator m = massrange.begin(); m != massrange.end(); ++m) {
289  gOptDMMasses.push_back(atof((*m).c_str()));
290  }
291  }
292  else {
293  gOptDMMasses.push_back(atof(dmm.c_str()));
294  }
295  } else {
296  LOG("gmkspl_dm", pFATAL)
297  << "Unspecified dark matter masses - Exiting";
298  PrintSyntax();
299  exit(1);
300  }
301 
302  // get the mediator mass ratio
303  if( parser.OptionExists('z') ) {
304  LOG("gmkspl_dm", pINFO) << "Reading mediator mass ratios";
305  string mr = parser.ArgAsString('z');
306 
307  if(mr.find(",") != string::npos) {
308  vector<string> ratiorange = utils::str::Split(mr, ",");
309  assert(ratiorange.size() > 1);
310  for (vector<string>::iterator r = ratiorange.begin(); r != ratiorange.end(); ++r) {
311  gOptMedRatios.push_back(atof((*r).c_str()));
312  }
313  }
314  else {
315  gOptMedRatios.push_back(atof(mr.c_str()));
316  }
317  } else {
318  LOG("gmkspl_dm", pINFO)
319  << "Unspecified mediator mass ratios - Using default";
320  gOptMedRatios.push_back(0.5);
321  }
322 
323  // write out input splines?
324  if( parser.OptionExists("no-copy") ) {
325  LOG("gmkspl_dm", pINFO) << "Not copying input splines to output";
326  gOptNoCopy = true;
327  }
328 
329  // get the mediator coupling
330  if( parser.OptionExists('g') ) {
331  LOG("gmkspl_dm", pINFO) << "Reading mediator couplings";
332  string gzp = parser.ArgAsString('g');
333 
334  if(gzp.find(",") != string::npos) {
335  vector<string> couprange = utils::str::Split(gzp, ",");
336  assert(couprange.size() > 1);
337  for (vector<string>::iterator g = couprange.begin(); g != couprange.end(); ++g) {
338  gOptZpCouplings.push_back(atof((*g).c_str()));
339  }
340  }
341  else {
342  gOptZpCouplings.push_back(atof(gzp.c_str()));
343  }
344  } else {
345  LOG("gmkspl_dm", pNOTICE)
346  << "Unspecified mediator couplings - Using value from config";
347  gOptZpCouplings.push_back(-1.);
348  }
349 
350  // write out input splines?
351  if( parser.OptionExists("no-copy") ) {
352  LOG("gmkspl_dm", pINFO) << "Not copying input splines to output";
353  gOptNoCopy = true;
354  }
355 
356 
357  // comma-separated target PDG code list or input geometry file
358  bool tgt_cmd = true;
359  if( parser.OptionExists('t') ) {
360  LOG("gmkspl_dm", pINFO) << "Reading target nuclei PDG codes";
361  gOptTgtPdgCodeList = parser.ArgAsString('t');
362  } else {
363  LOG("gmkspl_dm", pINFO) << "No code list specified from the command line";
364  tgt_cmd = false;
365  }
366 
367  bool tgt_geom = true;
368  if( parser.OptionExists('f') ) {
369  LOG("gmkspl_dm", pINFO) << "Reading ROOT geometry filename";
370  gOptGeomFilename = parser.ArgAsString('f');
371  } else {
372  LOG("gmkspl_dm", pINFO) << "No geometry file was specified";
373  tgt_cmd = false;
374  }
375 
376  bool both = tgt_geom && tgt_cmd;
377  bool none = !tgt_geom && !tgt_cmd;
378  if(none) {
379  LOG("gmkspl_dm", pFATAL)
380  << "No geom file or cmd line target list was specified - Exiting";
381  PrintSyntax();
382  exit(1);
383  }
384  if(both) {
385  LOG("gmkspl_dm", pFATAL)
386  << "You specified both a geom file and a cmd line target list "
387  << "- Exiting confused";
388  PrintSyntax();
389  exit(1);
390  }
391 
392  // random number seed
393  if( parser.OptionExists("seed") ) {
394  LOG("gmkspl_dm", pINFO) << "Reading random number seed";
395  gOptRanSeed = parser.ArgAsLong("seed");
396  } else {
397  LOG("gmkspl_dm", pINFO) << "Unspecified random number seed - Using default";
398  gOptRanSeed = -1;
399  }
400 
401  // input cross-section file
402  if( parser.OptionExists("input-cross-sections") ) {
403  LOG("gmkspl_dm", pINFO) << "Reading cross-section file";
404  gOptInpXSecFile = parser.ArgAsString("input-cross-sections");
405  } else {
406  LOG("gmkspl_dm", pINFO) << "Unspecified input cross-section file";
407  gOptInpXSecFile = "";
408  }
409 
410  //
411  // print the command-line options
412  //
413  LOG("gmkspl_dm", pNOTICE)
414  << "\n"
415  << utils::print::PrintFramedMesg("gmkspl_dm job configuration")
416  << "\n Target PDG codes : " << gOptTgtPdgCodeList
417  << "\n Input ROOT geometry : " << gOptGeomFilename
418  << "\n Output cross-section file : " << gOptOutXSecFile
419  << "\n Input cross-section file : " << gOptInpXSecFile
420  << "\n Random number seed : " << gOptRanSeed
421  << "\n";
422 
423  // print list of DM masses
424  for (vector<double>::iterator m = gOptDMMasses.begin(); m != gOptDMMasses.end(); ++m) {
425  LOG("gmkspl_dm",pNOTICE)
426  << "Considering DM mass : " << *m;
427  }
428  for (vector<double>::iterator r = gOptMedRatios.begin(); r != gOptMedRatios.end(); ++r) {
429  LOG("gmkspl_dm",pNOTICE)
430  << "Considering mediator mass ratio : " << *r;
431  }
432  for (vector<double>::iterator g = gOptZpCouplings.begin(); g != gOptZpCouplings.end(); ++g) {
433  LOG("gmkspl_dm",pNOTICE)
434  << "Considering Z' couplings : " << *g;
435  }
436  LOG("gmkspl_dm", pNOTICE) << *RunOpt::Instance();
437 }
438 //____________________________________________________________________________
439 void PrintSyntax(void)
440 {
441  LOG("gmkspl_dm", pNOTICE)
442  << "\n\n" << "Syntax:" << "\n"
443  << " gmkspl_dm -m dm_masses "
444  << " <-t tgtpdg, -f geomfile> "
445  << " <-o | --output-cross-section> xsec_xml_file_name"
446  << " [-g zp_couplings] "
447  << " [-z med_ratios] "
448  << " [-n nknots] [-e max_energy] "
449  << " [--seed seed_number]"
450  << " [--input-cross-section xml_file]"
451  << " [--event-generator-list list_name]"
452  << " [--message-thresholds xml_file]\n\n";
453 }
454 //____________________________________________________________________________
456 {
457  bool from_geom_file = ( gOptGeomFilename.size() > 0 );
458  bool from_cmd_line = ( gOptTgtPdgCodeList.size() > 0 );
459 
460  if (from_cmd_line) {
461  // split the comma separated list
462  vector<string> tgtvec = utils::str::Split(gOptTgtPdgCodeList, ",");
463 
464  // fill in the PDG code list
465  PDGCodeList * list = new PDGCodeList;
466  vector<string>::const_iterator iter;
467  for(iter = tgtvec.begin(); iter != tgtvec.end(); ++iter) {
468  list->push_back( atoi(iter->c_str()) );
469  }
470  return list;
471  }
472 
473  if (from_geom_file) {
474 #ifdef __GENIE_GEOM_DRIVERS_ENABLED__
475  // create/configure a geometry driver
476  LOG("gmkspl_dm", pINFO) << "Creating/configuring a ROOT geom. driver";
478 
479  PDGCodeList * list = new PDGCodeList(geom->ListOfTargetNuclei());
480 
481  delete geom;
482  return list;
483 #else
484  LOG("gmkspl_dm", pFATAL)
485  << "To read-in a ROOT geometry you need to enable the geometry drivers!";
486  gAbortingInErr = true;
487  exit(1);
488  return 0;
489 #endif
490 
491  }
492  return 0;
493 }
494 //____________________________________________________________________________
void RandGen(long int seed)
Definition: AppInit.cxx:30
PDGCodeList * GetTargetCodes(void)
void XSecTable(string inpfile, bool require_table)
Definition: AppInit.cxx:38
int gOptNKnots
static constexpr double g
Definition: Units.h:144
void AddDarkMatter(double mass, double med_ratio)
Definition: PDGLibrary.cxx:153
void ReadFromCommandLine(int argc, char **argv)
Definition: RunOpt.cxx:99
#define pFATAL
Definition: Messenger.h:56
const int kPdgDarkMatter
Definition: PDGCodes.h:218
vector< double > gOptDMMasses
static XSecSplineList * Instance()
virtual const PDGCodeList & ListOfTargetNuclei(void)
implement the GeomAnalyzerI interface
Registry * CommonList(const string &file_id, const string &set_name) const
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:327
A list of PDG codes.
Definition: PDGCodeList.h:32
string gOptOutXSecFile
string gOptInpXSecFile
Definition: gAtmoEvGen.cxx:313
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void SetEventGeneratorList(string listname)
Definition: GEVGDriver.cxx:348
vector< double > gOptMedRatios
string geom
string gOptTgtPdgCodeList
vector< double > gOptZpCouplings
A ROOT/GEANT4 geometry driver.
GENIE Event Generation Driver. A minimalist user interface object for generating neutrino interaction...
Definition: GEVGDriver.h:54
#define pINFO
Definition: Messenger.h:62
void Lock(void)
locks the registry
Definition: Registry.cxx:148
void BuildTune()
build tune and inform XSecSplineList
Definition: RunOpt.cxx:92
void SaveAsXml(const string &filename, bool save_init=true) const
void UnLock(void)
unlocks the registry (doesn&#39;t unlock items)
Definition: Registry.cxx:153
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
static RunOpt * Instance(void)
Definition: RunOpt.cxx:54
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
void Configure(int nu_pdgc, int Z, int A)
Definition: GEVGDriver.cxx:137
double gOptMaxE
string PrintFramedMesg(string mesg, unsigned int nl=1, const char f='*')
Definition: PrintUtils.cxx:164
A vector of EventGeneratorI objects.
void CreateSplines(int nknots=-1, double emax=-1, bool inLogE=true)
Definition: GEVGDriver.cxx:577
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:99
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:61
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:563
bool gOptNoCopy
bool gAbortingInErr
Definition: Messenger.cxx:34
List of cross section vs energy splines.
void Set(RgIMapPair entry)
Definition: Registry.cxx:267
static constexpr double m
Definition: Units.h:71
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
void PrintSyntax(void)
void EnableBareXSecPreCalc(bool flag)
Definition: RunOpt.h:62
void ReloadDBase(void)
Definition: PDGLibrary.cxx:224
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
static AlgConfigPool * Instance()
string gOptGeomFilename
Initial State information.
Definition: InitialState.h:48
long int gOptRanSeed
Definition: gAtmoEvGen.cxx:312