GMCJDriver.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Feb 08, 2008 - CA
14  Modified the global probability scale to be the maximum amongst the maximum
15  interaction probabilities for each neutrino (rather than the sum of maximum
16  probabilities). The modified probability scale still gives unbiased event
17  generation & reduces the 'no-interaction' probability.
18  @ Feb 14, 2008 - CA
19  Significant speed improvements - Most of the rejected flux neutrinos, are
20  rejected before having them propagated through the detector geometry
21  (flux neutrinos are pre-selected using the maximum path-lengths / for most
22  realistic fluxes -high energy tail, low energy peak- most of the selection
23  inefficiency is caused not because the path-lengths are not close the max
24  possible ones, but because the energy is not close to the max possible one).
25  Also some speed improvement was gained by properly using the total cross
26  section splines (before, han repeatedly summing-up the
27  The driver does not assert that _each_ flux neutrino generation & geometry
28  navigation will be succesfull - In the rare event that this may happen, it
29  prints an err mesg and tries again. In next revision I will limit the
30  number of successive trials something may go wrong to prevent the driver
31  from hanging in truly problematic cases.
32  The driver was adapted to handle flux drivers that -at some point- may stop
33  generating more flux neutrinos (eg because they read flux neutrinos by
34  looping over a beam simulation ntuple and they reached its last entry).
35  Code was appropriately restructured and some methods have been renamed.
36  @ Feb 29, 2008 - CA
37  Modified the InteractionProbability() to calculate absolute interaction
38  probabilities. Added NFluxNeutrinos() and GlobProbScale() to get the
39  number of neutrinos thrown by the flux driver towards the geometry and
40  the global interaction probability scale so as to be able to calculate
41  event sample normalization factors.
42  @ Jan 15, 2009 - CA
43  Stopped GMCJDriver from initializing the unphysical event mask so as not
44  to overwrite the values that each GEVGDriver obtains from the environment.
45  @ Jan 16, 2009 - CA
46  Added methods to return pointers to the flux and geometry drivers.
47  @ Mar 11, 2009 - CA
48  In GenerateEvent1Try() handle failure to generate kinematics. Added sanity
49  check on the no interaction probability.
50  @ Mar 04, 2010 - CA
51  Remove unused FilterUnphysical(TBits) method. Now set exclusively via the
52  GUNPHYSMASK env.var.
53  @ Apr 15, 2010 - CA
54  Fix unit error in ComputeEventProbability() - Reported by Corey Reed.
55  The probability stored at the output event was wrong but this doesn't
56  affect any of the existing applications as this number wasn't actually
57  used anywhere.
58  @ Dec 07, 2010 - CA
59  Don't use a fixed bin size in ComputeProbScales() as this was causing
60  errors for low energy applications. Addresses a problem reported by
61  Joachim Kopp.
62  @ Feb 22, 2011 - JD
63  Added a number of new methods to allow pre-calculation of exact flux
64  interaction probabilities for a given set of flux neutrinos from the
65  flux driver. See the comments for the new LoadFluxProbabilities,
66  SaveFluxProbabilities, PreCalcFluxProbabilities and PreSelectEvents
67  methods for details. Using these methods mean that there is no need
68  to generate maximum path lengths as instead use the exact interaction
69  probabilities to pre-select. This can result in very significant speed
70  increases (between factor of 5 and ~300) for event generation over complex
71  detector geometries and with realistic flux drivers. See
72  src/support/t2k/EvGen/gT2KEvGen.cxx for an example of how to use.
73  @ Mar, 7, 2011 - JD
74  Store sum totals of the flux interaction probabilities for various neutrino
75  type in a map relating pdg code to total interaction probability. Also add
76  public getter method so that this can be used in applications to work out
77  expected event rates. See gT2KEvGen.cxx for an example of how to do this.
78  Also save the PDG code for each entry in the flux interaction probabilities
79  tree.
80  @ Mar, 11, 2011 - JD
81  Set the directory of fFluxIntTree to the output file fFluxIntProbFile if
82  saving it later. This is so that it is incrementally saved and fixes bug
83  where getting std::bad_alloc when trying to Write large trees
84  fFluxIntProbFile.
85  @ Jan 31, 2013 - CA
86  Added SetEventGeneratorList(string listname). $GEVGL var no longer in use.
87  @ Feb 01, 2013 - CA
88  The GUNPHYSMASK env. var is no longer used. Added SetUnphysEventMask(const
89  TBits &). Input is propagated accordingly.
90  @ Feb 06, 2013 - CA
91  Fix small problem introduced with recent changes.
92  In PopulateEventGenDriverPool() calls to GEVGDriver::SetEventGeneratorList()
93  and GEVGDriver::Configure() were reversed. Problem reported by W.Huelsnitz.
94  @ July 15, 2014 - HG
95  Incorporated code provided by Jason Koskinen - IceCube
96  Modified ComputeProbScales to evalulate the cross sections at both the high
97  and low edges of the energy bin when calculating the max interaction
98  probability.
99 */
100 //____________________________________________________________________________
101 
102 #include <cassert>
103 
104 #include <TVector3.h>
105 #include <TSystem.h>
106 #include <TStopwatch.h>
107 
129 
130 using namespace genie;
131 using namespace genie::constants;
132 
133 //____________________________________________________________________________
135 {
136  this->InitJob();
137 }
138 //___________________________________________________________________________
140 {
141  if(fUnphysEventMask) delete fUnphysEventMask;
142  if (fGPool) delete fGPool;
143 
144  map<int,TH1D*>::iterator pmax_iter = fPmax.begin();
145  for( ; pmax_iter != fPmax.end(); ++pmax_iter) {
146  TH1D * pmax = pmax_iter->second;
147  if(pmax) {
148  delete pmax; pmax = 0;
149  }
150  }
151  fPmax.clear();
152 
153  if(fFluxIntTree) delete fFluxIntTree;
154  if(fFluxIntProbFile) delete fFluxIntProbFile;
155 }
156 //___________________________________________________________________________
158 {
159  LOG("GMCJDriver", pNOTICE)
160  << "Setting event generator list: " << listname;
161 
162  fEventGenList = listname;
163 }
164 //___________________________________________________________________________
166 {
167  *fUnphysEventMask = mask;
168 
169  LOG("GMCJDriver", pNOTICE)
170  << "Setting unphysical event mask (bits: " << GHepFlags::NFlags() - 1
171  << " -> 0) : " << *fUnphysEventMask;
172 }
173 //___________________________________________________________________________
175 {
176  fFluxDriver = flux_driver;
177 }
178 //___________________________________________________________________________
180 {
181  fGeomAnalyzer = geom_analyzer;
182 }
183 //___________________________________________________________________________
184 void GMCJDriver::UseSplines(bool useLogE)
185 {
186  fUseSplines = true;
187  fUseLogE = useLogE;
188 }
189 //___________________________________________________________________________
190 bool GMCJDriver::UseMaxPathLengths(string xml_filename)
191 {
192 // If you supply the maximum path lengths for all materials in your geometry
193 // (eg for ROOT/GEANT geometries they can be computed running GENIE's gmxpl
194 // application, see $GENIE/src/stdapp/gMaxPathLengths.cxx ) you can speed up
195 // the driver init phase by quite a bit (especially for complex geometries).
196 
197  fMaxPlXmlFilename = xml_filename;
198 
199  bool is_accessible = !(gSystem->AccessPathName(fMaxPlXmlFilename.c_str()));
200 
201  if ( is_accessible ) fUseExtMaxPl = true;
202  else {
203  fUseExtMaxPl = false;
204  LOG("GMCJDriver", pWARN)
205  << "UseMaxPathLengths could not find file: \"" << xml_filename << "\"";
206  }
207  return fUseExtMaxPl;
208 
209 }
210 //___________________________________________________________________________
212 {
213  LOG("GMCJDriver", pNOTICE)
214  << "Keep on throwing flux neutrinos till one interacts? : "
215  << utils::print::BoolAsYNString(keep_on);
216  fKeepThrowingFluxNu = keep_on;
217 }
218 //___________________________________________________________________________
220 {
221 // Use a single probability scale. That generates unweighted events.
222 // (Note that generating unweighted event kinematics is a different thing)
223 //
224  fGenerateUnweighted = true;
225 
226  LOG("GMCJDriver", pNOTICE)
227  << "GMCJDriver will generate un-weighted events. "
228  << "Note: That does not force unweighted event kinematics!";
229 }
230 //___________________________________________________________________________
231 void GMCJDriver::PreSelectEvents(bool preselect)
232 {
233 // Set whether to pre-select events based on a max-path lengths file. This
234 // should be turned off if using pre-generated interaction probabilities
235 // calculated from a given flux file.
236  fPreSelect = preselect;
237 }
238 //___________________________________________________________________________
240 {
241 // Loop over complete set of flux entries satisfying input config options
242 // (such as neutrino type) and save the interaction probability in a tree
243 // relating flux index (entry number in input flux tree) to interaction
244 // probability. If a pre-generated flux interaction probability tree has
245 // already been loaded then just returns true. Also save tree to a TFile
246 // for use in later jobs if flag is set
247 //
248  bool success = true;
249 
250  bool save_to_file = fFluxIntProbFile == 0 && fFluxIntFileName.size()>0;
251 
252  // Clear map storing sum(fBrFluxWeight*fBrFluxIntProb) for each neutrino pdg
253  fSumFluxIntProbs.clear();
254 
255  // check if already loaded flux interaction probs using LoadFluxProbTree
256  if(fFluxIntTree){
257  LOG("GMCJDriver", pNOTICE) <<
258  "Skipping pre-generation of flux interaction probabilities - "<<
259  "using pre-generated file";
260  success = true;
261  }
262  // otherwise create them on the fly now
263  else {
264 
265  if(save_to_file){
266  fFluxIntProbFile = new TFile(fFluxIntFileName.c_str(), "CREATE");
267  if(fFluxIntProbFile->IsZombie()){
268  LOG("GMCJDriver", pFATAL) << "Cannot overwrite an existing file. Exiting!";
269  exit(1);
270  }
271  }
272 
273  // Create the tree to store flux probs
274  fFluxIntTree = new TTree(fFluxIntTreeName.c_str(),
275  "Tree storing pre-calculated flux interaction probs");
276  fFluxIntTree->Branch("FluxIndex", &fBrFluxIndex, "FluxIndex/I");
277  fFluxIntTree->Branch("FluxIntProb", &fBrFluxIntProb, "FluxIntProb/D");
278  fFluxIntTree->Branch("FluxEnu", &fBrFluxEnu, "FluxEnu/D");
279  fFluxIntTree->Branch("FluxWeight", &fBrFluxWeight, "FluxWeight/D");
280  fFluxIntTree->Branch("FluxPDG", &fBrFluxPDG, "FluxPDG/I");
281  // Associate to file otherwise get std::bad_alloc when writing large trees
282  if(save_to_file) fFluxIntTree->SetDirectory(fFluxIntProbFile);
283 
284  fFluxDriver->GenerateWeighted(true);
285 
286  fGlobPmax = 1.0; // Force ComputeInteractionProbabilities to return absolute value
287 
288  // Loop over flux entries and calculate interaction probabilities
289  TStopwatch stopwatch;
290  stopwatch.Start();
291  long int first_index = -1;
292  bool first_loop = true;
293  // loop until at end of flux ntuple
294  while(fFluxDriver->End() == false){
295 
296  // get the next flux neutrino
297  bool gotnext = fFluxDriver->GenerateNext();
298  if(!gotnext){
299  LOG("GMCJDriver", pWARN) << "*** Couldn't generate next flux ray! ";
300  continue;
301  }
302 
303  // stop if completed a full cycle (this check is necessary as fluxdriver
304  // may be set to loop over more than one cycle before reaching end)
305  bool already_been_here = first_loop ? false : first_index == fFluxDriver->Index();
306  if(already_been_here) break;
307 
308  // compute the path lengths for current flux neutrino
309  if(this->ComputePathLengths() == false){ success = false; break;}
310 
311  // compute and store the interaction probability
312  double psum = this->ComputeInteractionProbabilities(false /*Based on actual PLs*/);
313  assert(psum+controls::kASmallNum > 0.);
314  fBrFluxIntProb = psum;
315  fBrFluxIndex = fFluxDriver->Index();
316  fBrFluxEnu = fFluxDriver->Momentum().E();
317  fBrFluxWeight = fFluxDriver->Weight();
318  fBrFluxPDG = fFluxDriver->PdgCode();
319  fFluxIntTree->Fill();
320 
321  // store the first index so know when have cycled exactly once
322  if(first_loop){
323  first_index = fFluxDriver->Index();
324  first_loop = false;
325  }
326  } // flux loop
327  stopwatch.Stop();
328  LOG("GMCJDriver", pNOTICE)
329  << "Finished pre-calculating flux interaction probabilities. "
330  << "Total CPU time to process "<< fFluxIntTree->GetEntries()
331  << " entries: "<< stopwatch.CpuTime();
332 
333  // reset the flux driver so can be used at next stage. N.B. This
334  // should also reset flux driver to throw de-weighted flux neutrinos
335  fFluxDriver->Clear("CycleHistory");
336  }
337 
338  // If successfully calculated/loaded interaction probabilities then set global
339  // probability scale and, if requested, save tree to output file
340  if(success){
341  fGlobPmax = 0.0;
342  double safety_factor = 1.01;
343  for(int i = 0; i< fFluxIntTree->GetEntries(); i++){
344  fFluxIntTree->GetEntry(i);
345  // Check have non-negative probabilities
346  assert(fBrFluxIntProb+controls::kASmallNum > 0.0);
347  assert(fBrFluxWeight+controls::kASmallNum > 0.0);
348  // Update the global maximum
349  fGlobPmax = TMath::Max(fGlobPmax, fBrFluxIntProb*safety_factor);
350  // Update the sum of fBrFluxIntProb*fBrFluxWeight for different species
351  if(fSumFluxIntProbs.find(fBrFluxPDG) == fSumFluxIntProbs.end()){
352  fSumFluxIntProbs[fBrFluxPDG] = 0.0;
353  }
354  fSumFluxIntProbs[fBrFluxPDG] += fBrFluxIntProb * fBrFluxWeight;
355  }
356  LOG("GMCJDriver", pNOTICE) <<
357  "Updated global probability scale to fGlobPmax = "<< fGlobPmax;
358 
359  if(save_to_file){
360  LOG("GMCJDriver", pNOTICE) <<
361  "Saving pre-generated interaction probabilities to file: "<<
362  fFluxIntProbFile->GetName();
363  fFluxIntProbFile->cd();
364  fFluxIntTree->Write();
365  }
366 
367  // Also build index for use later
368  if(fFluxIntTree->BuildIndex("FluxIndex") != fFluxIntTree->GetEntries()){
369  LOG("GMCJDriver", pFATAL) <<
370  "Cannot build index using branch \"FluxIndex\" for flux prob tree!";
371  exit(1);
372  }
373 
374  // Now that have pre-generated flux probabilities need to trun off event
375  // preselection as this is only advantages when using max path lengths
376  this->PreSelectEvents(false);
377 
378  LOG("GMCJDriver", pNOTICE) << "Successfully generated/loaded pre-calculate flux interaction probabilities";
379  }
380  // Otherwise clean up
381  else if(fFluxIntTree){
382  delete fFluxIntTree;
383  fFluxIntTree = 0;
384  }
385 
386  // Return whether have successfully pre-calculated flux interaction probabilities
387  return success;
388 }
389 //___________________________________________________________________________
391 {
392 // Load a pre-generated set of flux interaction probabilities from an external
393 // file. This is recommended when using large flux files (>100k entries) as
394 // for these the time to calculate the interaction probabilities can exceed
395 // ~20 minutes. After loading the input tree we call PreCalcFluxProbabilities
396 // to check that has successfully loaded
397 //
398  if(fFluxIntProbFile){
399  LOG("GMCJDriver", pWARN)
400  << "Can't load flux interaction prob file as one is already loaded";
401  return false;
402  }
403 
404  fFluxIntProbFile = new TFile(filename.c_str(), "OPEN");
405 
406  if(fFluxIntProbFile){
407  fFluxIntTree = dynamic_cast<TTree*>(fFluxIntProbFile->Get(fFluxIntTreeName.c_str()));
408  if(fFluxIntTree){
409  bool set_addresses =
410  fFluxIntTree->SetBranchAddress("FluxIntProb", &fBrFluxIntProb) >= 0 &&
411  fFluxIntTree->SetBranchAddress("FluxIndex", &fBrFluxIndex) >= 0 &&
412  fFluxIntTree->SetBranchAddress("FluxPDG", &fBrFluxPDG) >= 0 &&
413  fFluxIntTree->SetBranchAddress("FluxWeight", &fBrFluxWeight) >= 0 &&
414  fFluxIntTree->SetBranchAddress("FluxEnu", &fBrFluxEnu) >= 0;
415  if(set_addresses){
416  // Finally check that can use them
417  if(this->PreCalcFluxProbabilities()) {
418  LOG("GMCJDriver", pNOTICE)
419  << "Successfully loaded pre-generated flux interaction probabilities";
420  return true;
421  }
422  }
423  // If cannot load then delete tree
424  LOG("GMCJDriver", pERROR) <<
425  "Cannot find expected branches in input flux probability tree!";
426  delete fFluxIntTree; fFluxIntTree = 0;
427  }
428  else LOG("GMCJDriver", pERROR)
429  << "Cannot find tree: "<< fFluxIntTreeName.c_str();
430  }
431 
432  LOG("GMCJDriver", pWARN)
433  << "Unable to load flux interaction probabilities file";
434  return false;
435 }
436 //___________________________________________________________________________
438 {
439 // Configue the flux driver to save the calculated flux interaction
440 // probabilities to the specified output file name for use in later jobs. See
441 // the LoadFluxProbTree method for how they are fed into a later job.
442 //
443  fFluxIntFileName = outfilename;
444 }
445 //___________________________________________________________________________
446 void GMCJDriver::Configure(bool calc_prob_scales)
447 {
448  LOG("GMCJDriver", pNOTICE)
449  << utils::print::PrintFramedMesg("Configuring GMCJDriver");
450 
451  // Get the list of neutrino types from the input flux driver and the list
452  // of target materials from the input geometry driver
453  this->GetParticleLists();
454 
455  // Ask the input GFluxI for the max. neutrino energy (to compute Pmax)
456  this->GetMaxFluxEnergy();
457 
458  // Create all possible initial states and for each one initialize,
459  // configure & store an GEVGDriver event generation driver object.
460  // Once an 'initial state' has been selected from the input flux / geom,
461  // the responsibility for generating the neutrino interaction will be
462  // delegated to one of these drivers.
463  this->PopulateEventGenDriverPool();
464 
465  // If the user wants to use cross section splines in order to speed things
466  // up, then coordinate spline creation from all GEVGDriver objects pushed
467  // into GEVGPool. This will create all xsec splines needed for all (enabled)
468  // processes that can be simulated involving the particles in the input flux
469  // and geometry.
470  // Spline creation will be skipped for every spline that has been pre-loaded
471  // into the the XSecSplineList.
472  // Once more it is noted that computing cross section splines is a huge
473  // overhead. The user is encouraged to generate them in advance and load
474  // them into the XSecSplineList
475  this->BootstrapXSecSplines();
476 
477  // Create cross section splines describing the total interaction xsec
478  // for a given initial state (Create them by summing all xsec splines
479  // for each possible initial state)
480  this->BootstrapXSecSplineSummation();
481 
482  if(calc_prob_scales){
483  // Ask the input geometry driver to compute the max. path length for each
484  // material in the list of target materials (or load a precomputed list)
485  this->GetMaxPathLengthList();
486 
487  // Compute the max. interaction probability to scale all interaction
488  // probabilities to be computed by this driver
489  this->ComputeProbScales();
490  }
491  LOG("GMCJDriver", pNOTICE) << "Finished configuring GMCJDriver\n\n";
492 }
493 //___________________________________________________________________________
495 {
496  fEventGenList = "Default"; // <-- set of event generators to be loaded by this driver
497 
498  fUnphysEventMask = new TBits(GHepFlags::NFlags()); //<-- unphysical event mask
499  //fUnphysEventMask->ResetAllBits(true);
500  for(unsigned int i = 0; i < GHepFlags::NFlags(); i++) {
501  fUnphysEventMask->SetBitNumber(i, true);
502  }
503 
504  fFluxDriver = 0; // <-- flux driver
505  fGeomAnalyzer = 0; // <-- geometry driver
506  fGPool = 0; // <-- pool of GEVGDriver event generation drivers
507  fEmax = 0; // <-- maximum neutrino energy
508  fMaxPlXmlFilename = ""; // <-- XML file with external path lengths
509  fUseExtMaxPl = false;
510  fUseSplines = false;
511  fNFluxNeutrinos = 0; // <-- number of flux neutrinos thrown so far
512 
513  fGlobPmax = 0; // <-- maximum interaction probability (global prob scale)
514  fPmax.clear(); // <-- maximum interaction probability per neutrino & per energy bin
515 
516  fGenerateUnweighted = false; // <-- default opt to generate weighted events
517  fPreSelect = true; // <-- default to use pre-selection based on maximum path lengths
518 
519  fSelTgtPdg = 0;
520  fCurEvt = 0;
521  fCurVtx.SetXYZT(0.,0.,0.,0.);
522 
523  fFluxIntProbFile = 0;
524  fFluxIntTreeName = "gFlxIntProb";
525  fFluxIntFileName = "";
526  fFluxIntTree = 0;
527  fBrFluxIntProb = -1.;
528  fBrFluxIndex = -1;
529  fBrFluxEnu = -1.;
530  fBrFluxWeight = -1.;
531  fBrFluxPDG = 0;
532  fSumFluxIntProbs.clear();
533 
534  // Throw as many flux neutrinos as necessary till one has interacted
535  // so that GenerateEvent() never returns NULL (except when in error)
536  this->KeepOnThrowingFluxNeutrinos(true);
537 
538  // Allow the selected GEVGDriver to go into recursive mode and regenerate
539  // an interaction that turns out to be unphysical.
540  //TBits unphysmask(GHepFlags::NFlags());
541  //unphysmask.ResetAllBits(false);
542  //this->FilterUnphysical(unphysmask);
543 
544  // Force early initialization of singleton objects that are typically
545  // would be initialized at their first use later on.
546  // This is purely cosmetic and I do it to send the banner and some prolific
547  // initialization printout at the top.
550 
551  // Clear the target and neutrino lists
552  fNuList.clear();
553  fTgtList.clear();
554 
555  // Clear the maximum path length list
556  fMaxPathLengths.clear();
557  fCurPathLengths.clear();
558 }
559 //___________________________________________________________________________
561 {
562  // Get the list of flux neutrinos from the flux driver
563  LOG("GMCJDriver", pNOTICE)
564  << "Asking the flux driver for its list of neutrinos";
565  fNuList = fFluxDriver->FluxParticles();
566 
567  LOG("GMCJDriver", pNOTICE) << "Flux particles: " << fNuList;
568 
569  // Get the list of target materials from the geometry driver
570  LOG("GMCJDriver", pNOTICE)
571  << "Asking the geometry driver for its list of targets";
572  fTgtList = fGeomAnalyzer->ListOfTargetNuclei();
573 
574  LOG("GMCJDriver", pNOTICE) << "Target materials: " << fTgtList;
575 }
576 //___________________________________________________________________________
578 {
579  if(fUseExtMaxPl) {
580  LOG("GMCJDriver", pNOTICE)
581  << "Loading external max path-length list for input geometry from "
582  << fMaxPlXmlFilename;
583  fMaxPathLengths.LoadFromXml(fMaxPlXmlFilename);
584 
585  } else {
586  LOG("GMCJDriver", pNOTICE)
587  << "Querying the geometry driver to compute the max path-length list";
588  fMaxPathLengths = fGeomAnalyzer->ComputeMaxPathLengths();
589  }
590  // Print maximum path lengths & neutrino energy
591  LOG("GMCJDriver", pNOTICE)
592  << "Maximum path length list: " << fMaxPathLengths;
593 }
594 //___________________________________________________________________________
596 {
597  LOG("GMCJDriver", pNOTICE)
598  << "Querying the flux driver for the maximum energy of flux neutrinos";
599  fEmax = fFluxDriver->MaxEnergy();
600 
601  LOG("GMCJDriver", pNOTICE)
602  << "Maximum flux neutrino energy = " << fEmax << " GeV";
603 }
604 //___________________________________________________________________________
606 {
607  LOG("GMCJDriver", pDEBUG)
608  << "Creating GEVGPool & adding a GEVGDriver object per init-state";
609 
610  if (fGPool) delete fGPool;
611  fGPool = new GEVGPool;
612 
613  PDGCodeList::const_iterator nuiter;
614  PDGCodeList::const_iterator tgtiter;
615 
616  for(nuiter = fNuList.begin(); nuiter != fNuList.end(); ++nuiter) {
617  for(tgtiter = fTgtList.begin(); tgtiter != fTgtList.end(); ++tgtiter) {
618 
619  int target_pdgc = *tgtiter;
620  int neutrino_pdgc = *nuiter;
621 
622  InitialState init_state(target_pdgc, neutrino_pdgc);
623 
624  LOG("GMCJDriver", pNOTICE)
625  << "\n\n ---- Creating a GEVGDriver object configured for init-state: "
626  << init_state.AsString() << " ----\n\n";
627 
628  GEVGDriver * evgdriver = new GEVGDriver;
629  evgdriver->SetEventGeneratorList(fEventGenList); // specify list of generators
630  evgdriver->Configure(init_state);
631  evgdriver->UseSplines(); // check if all splines needed are loaded
632 
633  LOG("GMCJDriver", pDEBUG) << "Adding new GEVGDriver object to GEVGPool";
634  fGPool->insert( GEVGPool::value_type(init_state.AsString(), evgdriver) );
635  } // targets
636  } // neutrinos
637 
638  LOG("GMCJDriver", pNOTICE)
639  << "All necessary GEVGDriver object were pushed into GEVGPool\n";
640 }
641 //___________________________________________________________________________
643 {
644 // Bootstrap cross section spline generation by the event generation drivers
645 // that handle each initial state.
646 
647  if(!fUseSplines) return;
648 
649  LOG("GMCJDriver", pNOTICE)
650  << "Asking event generation drivers to compute all needed xsec splines";
651 
652  PDGCodeList::const_iterator nuiter;
653  PDGCodeList::const_iterator tgtiter;
654  for(nuiter = fNuList.begin(); nuiter != fNuList.end(); ++nuiter){
655  for(tgtiter = fTgtList.begin(); tgtiter != fTgtList.end(); ++tgtiter) {
656  int target_pdgc = *tgtiter;
657  int neutrino_pdgc = *nuiter;
658  InitialState init_state(target_pdgc, neutrino_pdgc);
659  LOG("GMCJDriver", pINFO)
660  << "Computing all splines needed for init-state: "
661  << init_state.AsString();
662  GEVGDriver * evgdriver = fGPool->FindDriver(init_state);
663  evgdriver->CreateSplines(-1,-1,fUseLogE);
664  } // targets
665  } // neutrinos
666  LOG("GMCJDriver", pINFO) << "Finished creating cross section splines\n";
667 }
668 //___________________________________________________________________________
670 {
671 // Sum-up the cross section splines for all the interaction that can be
672 // simulated for each initial state
673 
674  LOG("GMCJDriver", pNOTICE)
675  << "Summing-up splines to get total cross section for each init state";
676 
677  GEVGPool::iterator diter;
678  for(diter = fGPool->begin(); diter != fGPool->end(); ++diter) {
679  string init_state = diter->first;
680  GEVGDriver * evgdriver = diter->second;
681  assert(evgdriver);
682  LOG("GMCJDriver", pNOTICE)
683  << "**** Summing xsec splines for init-state = " << init_state;
684 
685  Range1D_t rE = evgdriver->ValidEnergyRange();
686  if (fEmax>rE.max || fEmax<rE.min)
687  LOG("GMCJDriver",pFATAL)
688  << " rE (validEnergyRange) [" << rE.min << "," << rE.max << "] "
689  << " fEmax " << fEmax;
690  assert(fEmax<rE.max && fEmax>rE.min);
691 
692  // decide the energy range for the sum spline - extend the spline a little
693  // bit above the maximum beam energy (but below the maximum valid energy)
694  // to avoid the evaluation of the cubic spline around the viscinity of
695  // knots with zero y values (although the GENIE Spline object handles it)
696  double dE = fEmax/10.;
697  double min = rE.min;
698  double max = (fEmax+dE < rE.max) ? fEmax+dE : rE.max;
699  evgdriver->CreateXSecSumSpline(100,min,max,true);
700  }
701  LOG("GMCJDriver", pNOTICE)
702  << "Finished summing all interaction xsec splines per initial state";
703 }
704 //___________________________________________________________________________
706 {
707 // Computing interaction probability scales.
708 // To minimize the numbers or trials before choosing a neutrino+target init
709 // state for generating an event (note: each trial means selecting a flux
710 // neutrino, navigating it through the detector to compute path lengths,
711 // computing interaction probabilities for each material and so on...)
712 // a set of probability scales can be used for different neutrino species
713 // and at different energy bins.
714 // A global probability scale is also being constructed for keeping the correct
715 // proportions between differect flux neutrino species or flux neutrinos of
716 // different energies.
717 
718  LOG("GMCJDriver", pNOTICE)
719  << "Computing the max. interaction probability (probability scale)";
720 
721  // clean up global probability scale and maximum probabilties per neutrino
722  // type & energy bin
723  {
724  fGlobPmax = 0;
725  map<int,TH1D*>::iterator pmax_iter = fPmax.begin();
726  for( ; pmax_iter != fPmax.end(); ++pmax_iter) {
727  TH1D * pmax = pmax_iter->second;
728  if(pmax) {
729  delete pmax; pmax = 0;
730  }
731  }
732  fPmax.clear();
733  }
734 
735  // for maximum interaction probability vs E /for given geometry/ I will
736  // be using 300 bins up to the maximum energy for the input flux
737  // double de = fEmax/300.;//djk june 5, 2013
738  double de = fEmax/300.;//djk june 5, 2013
739  double emin = 0.0;
740  double emax = fEmax + de;
741  int n = 1 + (int) ((emax-emin)/de);
742 
743  PDGCodeList::const_iterator nuiter;
744  PDGCodeList::const_iterator tgtiter;
745 
746  // loop over all neutrino types generated by the flux driver
747  for(nuiter = fNuList.begin(); nuiter != fNuList.end(); ++nuiter) {
748  int neutrino_pdgc = *nuiter;
749  TH1D * pmax_hst = new TH1D("pmax_hst",
750  "max interaction probability vs E | geom",n,emin,emax);
751  pmax_hst->SetDirectory(0);
752 
753  // loop over energy bins
754  for(int ie = 1; ie <= pmax_hst->GetNbinsX(); ie++) {
755  double EvLow = pmax_hst->GetBinCenter(ie) - 0.5*pmax_hst->GetBinWidth(ie);
756  double EvHigh = pmax_hst->GetBinCenter(ie) + 0.5*pmax_hst->GetBinWidth(ie);
757  //double Ev = pmax_hst->GetBinCenter(ie);
758 
759  // loop over targets in input geometry, form initial state and compute
760  // the sum of maximum interaction probabilities at the current energy bin
761  //
762  for(tgtiter = fTgtList.begin(); tgtiter != fTgtList.end(); ++tgtiter) {
763  int target_pdgc = *tgtiter;
764 
765  InitialState init_state(target_pdgc, neutrino_pdgc);
766 
767  LOG("GMCJDriver", pDEBUG)
768  << "Computing Pmax for init-state: " << init_state.AsString() << " E from " << EvLow << "-" << EvHigh;
769 
770  // get the appropriate driver
771  GEVGDriver * evgdriver = fGPool->FindDriver(init_state);
772 
773  // get xsec sum over all modelled processes for given neutrino+target)
774  double sxsecLow = evgdriver->XSecSumSpline()->Evaluate(EvLow);
775  double sxsecHigh = evgdriver->XSecSumSpline()->Evaluate(EvHigh);
776 
777  // get max{path-length x density}
778  double plmax = fMaxPathLengths.PathLength(target_pdgc);
779 
780  // compute/store the max interaction probabiity (for given energy)
781  int A = pdg::IonPdgCodeToA(target_pdgc);
782  double pmaxLow = this->InteractionProbability(sxsecLow, plmax, A);
783  double pmaxHigh = this->InteractionProbability(sxsecHigh, plmax, A);
784 
785  double pmax = pmaxHigh;
786  if ( pmaxLow > pmaxHigh){
787  pmax = pmaxLow;
788  LOG("GMCJDriver", pWARN)
789  << "Lower energy neutrinos have a higher probability of interacting than those at higher energy."
790  << " pmaxLow(E=" << EvLow << ")=" << pmaxLow << " and " << " pmaxHigh(E=" << EvHigh << ")=" << pmaxHigh;
791  }
792 
793  pmax_hst->SetBinContent(ie, pmax_hst->GetBinContent(ie) + pmax);
794 
795  LOG("GMCJDriver", pDEBUG)
796  << "Pmax[" << init_state.AsString() << ", Ev from " << EvLow << "-" << EvHigh << "] = " << pmax;
797  } // targets
798 
799  pmax_hst->SetBinContent(ie, 1.2 * pmax_hst->GetBinContent(ie));
800 
801  LOG("GMCJDriver", pINFO)
802  << "Pmax[nu=" << neutrino_pdgc << ", Ev from " << EvLow << "-" << EvHigh << "] = "
803  << pmax_hst->GetBinContent(ie);
804  } // E
805 
806  fPmax.insert(map<int,TH1D*>::value_type(neutrino_pdgc,pmax_hst));
807  } // nu
808 
809  // Compute global probability scale
810  // Sum Probabilities {
811  // all neutrinos, all targets, @ max path length, @ max energy}
812  //
813  {
814  for(nuiter = fNuList.begin(); nuiter != fNuList.end(); ++nuiter) {
815  int neutrino_pdgc = *nuiter;
816  map<int,TH1D*>::const_iterator pmax_iter = fPmax.find(neutrino_pdgc);
817  assert(pmax_iter != fPmax.end());
818  TH1D * pmax_hst = pmax_iter->second;
819  assert(pmax_hst);
820 // double pmax = pmax_hst->GetBinContent(pmax_hst->FindBin(fEmax));
821  double pmax = pmax_hst->GetMaximum();
822  assert(pmax>0);
823 // fGlobPmax += pmax;
824  fGlobPmax = TMath::Max(pmax, fGlobPmax); // ?;
825  }
826  LOG("GMCJDriver", pNOTICE) << "*** Probability scale = " << fGlobPmax;
827  }
828 }
829 //___________________________________________________________________________
831 {
832  fCurPathLengths.clear();
833  fCurEvt = 0;
834  fSelTgtPdg = 0;
835  fCurVtx.SetXYZT(0.,0.,0.,0.);
836 }
837 //___________________________________________________________________________
839 {
840  LOG("GMCJDriver", pNOTICE) << "Generating next event...";
841 
842  this->InitEventGeneration();
843 
844  while(1) {
845  bool flux_end = fFluxDriver->End();
846  if(flux_end) {
847  LOG("GMCJDriver", pNOTICE)
848  << "No more neutrinos can be thrown by the flux driver";
849  return 0;
850  }
851 
852  EventRecord * event = this->GenerateEvent1Try();
853  if(event) return event;
854 
855  if(fKeepThrowingFluxNu) {
856  LOG("GMCJDriver", pNOTICE)
857  << "Flux neutrino didn't interact - Trying the next one...";
858  continue;
859  }
860  break;
861  } // (w(1)
862 
863  LOG("GMCJDriver", pINFO) << "Returning NULL event!";
864  return 0;
865 }
866 //___________________________________________________________________________
868 {
869 // attempt generating a neutrino interaction by firing a single flux neutrino
870 //
872 
873  double Pno=0, Psum=0;
874  double R = rnd->RndEvg().Rndm();
875  LOG("GMCJDriver", pDEBUG) << "Rndm [0,1] = " << R;
876 
877  // Generate a neutrino using the input GFluxI & get current pdgc/p4/x4
878  bool flux_ok = this->GenerateFluxNeutrino();
879  if(!flux_ok) {
880  LOG("GMCJDriver", pERROR)
881  << "** Rejecting current flux neutrino (flux driver err)";
882  return 0;
883  }
884 
885  // Compute the interaction probabilities assuming max. path lengths
886  // and decide whether the neutrino would interact --
887  // Many flux neutrinos should be rejected here, drastically reducing
888  // the number of neutrinos that I need to propagate through the
889  // actual detector geometry (this is skipped when using
890  // pre-calculated flux interaction probabilities)
891  if(fPreSelect) {
892  LOG("GMCJDriver", pNOTICE)
893  << "Computing interaction probabilities for max. path lengths";
894 
895  Psum = this->ComputeInteractionProbabilities(true /* <- max PL*/);
896  Pno = 1-Psum;
897  LOG("GMCJDriver", pNOTICE)
898  << "The no-interaction probability (max. path lengths) is: "
899  << 100*Pno << " %";
900  if(Pno<0.) {
901  LOG("GMCJDriver", pFATAL)
902  << "Negative no-interaction probability! (P = " << 100*Pno << " %)"
903  << " Particle E=" << fFluxDriver->Momentum().E() << " type=" << fFluxDriver->PdgCode() << "Psum=" << Psum;
904  gAbortingInErr=true;
905  exit(1);
906  }
907  if(R>=1-Pno) {
908  LOG("GMCJDriver", pNOTICE)
909  << "** Rejecting current flux neutrino";
910  return 0;
911  }
912  } // preselect
913 
914  bool pl_ok = false;
915 
916 
917  // If possible use pre-generated flux neutrino interaction probabilities
918  if(fFluxIntTree){
919  Psum = this->PreGenFluxInteractionProbability();
920  }
921  // Else compute them in the usual manner
922  else {
923  // Compute (pathLength x density x weight fraction) for all materials
924  // in the input geometry, for the neutrino generated by the flux driver
925  pl_ok = this->ComputePathLengths();
926  if(!pl_ok) {
927  LOG("GMCJDriver", pERROR)
928  << "** Rejecting current flux neutrino (err computing path-lengths)";
929  return 0;
930  }
931  if(fCurPathLengths.AreAllZero()) {
932  LOG("GMCJDriver", pNOTICE)
933  << "** Rejecting current flux neutrino (misses generation volume)";
934  return 0;
935  }
936  Psum = this->ComputeInteractionProbabilities(false /* <- actual PL */);
937  }
938 
939 
940  if(TMath::Abs(Psum) < controls::kASmallNum){
941  LOG("GMCJDriver", pNOTICE)
942  << "** Rejecting current flux neutrino (has null interaction probability)";
943  return 0;
944  }
945 
946  // Now decide whether the current neutrino interacts
947  Pno = 1-Psum;
948  LOG("GMCJDriver", pNOTICE)
949  << "The actual 'no interaction' probability is: " << 100*Pno << " %";
950  if(Pno<0.) {
951  LOG("GMCJDriver", pFATAL)
952  << "Negative no interactin probability! (P = " << 100*Pno << " %)";
953 
954  // print info about what caused the problem
955  int nupdg = fFluxDriver -> PdgCode ();
956  const TLorentzVector & nup4 = fFluxDriver -> Momentum ();
957  const TLorentzVector & nux4 = fFluxDriver -> Position ();
958 
959  LOG("GMCJDriver", pWARN)
960  << "\n [-] Problematic neutrino: "
961  << "\n |----o PDG-code : " << nupdg
962  << "\n |----o 4-momentum : " << utils::print::P4AsString(&nup4)
963  << "\n |----o 4-position : " << utils::print::X4AsString(&nux4)
964  << "\n Emax : " << fEmax;
965 
966  LOG("GMCJDriver", pWARN)
967  << "\n Problematic path lengths:" << fCurPathLengths;
968 
969  LOG("GMCJDriver", pWARN)
970  << "\n Maximum path lengths:" << fMaxPathLengths;
971 
972  exit(1);
973  }
974  if(R>=1-Pno) {
975  LOG("GMCJDriver", pNOTICE)
976  << "** Rejecting current flux neutrino";
977  return 0;
978  }
979 
980  //
981  // The flux neutrino interacts!
982  //
983 
984  // Calculate path lengths for first time and check potential mismatch if
985  // used pre-generated flux interaction probabilities
986  if(fFluxIntTree){
987  pl_ok = this->ComputePathLengths();
988  if(!pl_ok) {
989  LOG("GMCJDriver", pFATAL) << "** Cannot calculate path lenths!";
990  exit(1);
991  }
992  double Psum_curr = this->ComputeInteractionProbabilities(false /* <- actual PL */);
993  bool mismatch = TMath::Abs(Psum-Psum_curr) > controls::kASmallNum;
994  if(mismatch){
995  LOG("GMCJDriver", pFATAL) <<
996  "** Mismatch between pre-calculated and current interaction "<<
997  "probabilities!";
998  exit(1);
999  }
1000  }
1001 
1002  // Select a target material
1003  fSelTgtPdg = this->SelectTargetMaterial(R);
1004  if(fSelTgtPdg==0) {
1005  LOG("GMCJDriver", pERROR)
1006  << "** Rejecting current flux neutrino (failed to select tgt!)";
1007  return 0;
1008  }
1009 
1010  // Ask the GEVGDriver object to select and generate an interaction and
1011  // its kinematics for the selected initial state & neutrino 4-momentum
1012  this->GenerateEventKinematics();
1013  if(!fCurEvt) {
1014  LOG("GMCJDriver", pWARN)
1015  << "** Couldn't generate kinematics for selected interaction";
1016  return 0;
1017  }
1018 
1019  // Generate an 'interaction position' in the selected material (in the
1020  // detector coord system), along the direction of nup4 & set it
1021  this->GenerateVertexPosition();
1022 
1023  // Set the event probability (probability for this event to happen given
1024  // the detector setup & the selected flux neutrino)
1025  // Note for users:
1026  // The above probability is stored at GHepRecord::Probability()
1027  // For normalization purposes make sure that you take into account the
1028  // GHepRecord::Weight() -if event generation is weighted-, and
1029  // GFluxI::Weight() -if beam simulation is weighted-.
1030  this->ComputeEventProbability();
1031 
1032  return fCurEvt;
1033 }
1034 //___________________________________________________________________________
1036 {
1037 // Ask the neutrino flux driver to generate a flux neutrino and make sure
1038 // that things look ok...
1039 //
1040  LOG("GMCJDriver", pNOTICE) << "Generating a flux neutrino";
1041 
1042  bool ok = fFluxDriver->GenerateNext();
1043  if(!ok) {
1044  LOG("GMCJDriver", pERROR)
1045  << "*** The flux driver couldn't generate a flux neutrino!!";
1046  return false;
1047  }
1048 
1049  fNFluxNeutrinos++;
1050  int nupdg = fFluxDriver -> PdgCode ();
1051  const TLorentzVector & nup4 = fFluxDriver -> Momentum ();
1052  const TLorentzVector & nux4 = fFluxDriver -> Position ();
1053 
1054  LOG("GMCJDriver", pNOTICE)
1055  << "\n [-] Generated flux neutrino: "
1056  << "\n |----o PDG-code : " << nupdg
1057  << "\n |----o 4-momentum : " << utils::print::P4AsString(&nup4)
1058  << "\n |----o 4-position : " << utils::print::X4AsString(&nux4);
1059 
1060  if(nup4.Energy() > fEmax) {
1061  LOG("GMCJDriver", pFATAL)
1062  << "\n *** Flux driver error ***"
1063  << "\n Generated flux v with E = " << nup4.Energy() << " GeV"
1064  << "\n Max v energy (declared by flux driver) = " << fEmax << " GeV"
1065  << "\n My interaction probability scaling is invalidated!!";
1066  return false;
1067  }
1068  if(!fNuList.ExistsInPDGCodeList(nupdg)) {
1069  LOG("GMCJDriver", pFATAL)
1070  << "\n *** Flux driver error ***"
1071  << "\n Generated flux v with pdg = " << nupdg
1072  << "\n It does not belong to the declared list of flux neutrinos"
1073  << "\n I was not configured to handle this!!";
1074  return false;
1075  }
1076  return true;
1077 }
1078 //___________________________________________________________________________
1080 {
1081 // Ask the geometry driver to compute (pathLength x density x weight frac.)
1082 // for all detector materials for the neutrino generated by the flux driver
1083 // and make sure that things look ok...
1084 
1085  fCurPathLengths.clear();
1086 
1087  const TLorentzVector & nup4 = fFluxDriver -> Momentum ();
1088  const TLorentzVector & nux4 = fFluxDriver -> Position ();
1089 
1090  fCurPathLengths = fGeomAnalyzer->ComputePathLengths(nux4, nup4);
1091 
1092  LOG("GMCJDriver", pNOTICE) << fCurPathLengths;
1093 
1094  if(fCurPathLengths.size() == 0) {
1095  LOG("GMCJDriver", pFATAL)
1096  << "\n *** Geometry driver error ***"
1097  << "\n Got an empty PathLengthList - No material found in geometry?";
1098  return false;
1099  }
1100 
1101  if(fCurPathLengths.AreAllZero()) {
1102  LOG("GMCJDriver", pNOTICE)
1103  << "current flux v doesn't cross any geometry material...";
1104  }
1105  return true;
1106 }
1107 //___________________________________________________________________________
1108 double GMCJDriver::ComputeInteractionProbabilities(bool use_max_path_length)
1109 {
1110  LOG("GMCJDriver", pNOTICE)
1111  << "Computing relative interaction probabilities for each material";
1112 
1113  // current flux neutrino code & 4-p
1114  int nupdg = fFluxDriver->PdgCode();
1115  const TLorentzVector & nup4 = fFluxDriver->Momentum();
1116 
1117  fCurCumulProbMap.clear();
1118 
1119  const PathLengthList & path_length_list =
1120  (use_max_path_length) ? fMaxPathLengths : fCurPathLengths;
1121 
1122  double probsum=0;
1123  PathLengthList::const_iterator pliter;
1124 
1125  for(pliter = path_length_list.begin();
1126  pliter != path_length_list.end(); ++pliter) {
1127  int mpdg = pliter->first; // material PDG code
1128  double pl = pliter->second; // density x path-length
1129  int A = pdg::IonPdgCodeToA(mpdg);
1130  double xsec = 0.; // sum of xsecs for all modelled processes for given init state
1131  double prob = 0.; // interaction probability
1132  double probn = 0.; // normalized interaction probability
1133 
1134  // find the GEVGDriver object that is handling the current init state
1135  InitialState init_state(mpdg, nupdg);
1136  GEVGDriver * evgdriver = fGPool->FindDriver(init_state);
1137  if(!evgdriver) {
1138  LOG("GMCJDriver", pFATAL)
1139  << "\n * The MC Job driver isn't properly configured!"
1140  << "\n * No event generation driver could be found for init state: "
1141  << init_state.AsString();
1142  exit(1);
1143  }
1144  // compute the interaction xsec and probability (if path-length>0)
1145  if(pl>0.) {
1146  const Spline * totxsecspl = evgdriver->XSecSumSpline();
1147  if(!totxsecspl) {
1148  LOG("GMCJDriver", pFATAL)
1149  << "\n * The MC Job driver isn't properly configured!"
1150  << "\n * Couldn't retrieve total cross section spline for init state: "
1151  << init_state.AsString();
1152  exit(1);
1153  } else {
1154  xsec = totxsecspl->Evaluate( nup4.Energy() );
1155  }
1156  prob = this->InteractionProbability(xsec,pl,A);
1157  LOG("GMCJDriver", pDEBUG)
1158  << " (xsec, pl, A)=(" << xsec << "," << pl << "," << A << ")";
1159 
1160  // scale the interaction probability to the maximum one so as not
1161  // to have to throw few billions of flux neutrinos before getting
1162  // an interaction...
1163  double pmax = 0;
1164  if(fGenerateUnweighted) pmax = fGlobPmax;
1165  else {
1166  map<int,TH1D*>::const_iterator pmax_iter = fPmax.find(nupdg);
1167  assert(pmax_iter != fPmax.end());
1168  TH1D * pmax_hst = pmax_iter->second;
1169  assert(pmax_hst);
1170  int ie = pmax_hst->FindBin(nup4.Energy());
1171  pmax = pmax_hst->GetBinContent(ie);
1172  }
1173  assert(pmax>0);
1174  LOG("GMCJDriver", pDEBUG)
1175  << "Pmax=" << pmax;
1176  probn = prob/pmax;
1177  }
1178 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
1179  LOG("GMCJDriver", pNOTICE)
1180  << "tgt: " << mpdg << " -> TotXSec = "
1181  << xsec/units::cm2 << " cm^2, Norm.Prob = " << 100*probn << "%";
1182 #endif
1183 
1184  probsum += probn;
1185  fCurCumulProbMap.insert(map<int,double>::value_type(mpdg,probsum));
1186  }
1187  return probsum;
1188 }
1189 //___________________________________________________________________________
1191 {
1192 // Pick a target material using the pre-computed interaction probabilities
1193 // for a flux neutrino that has already been determined that interacts
1194 
1195  LOG("GMCJDriver", pNOTICE) << "Selecting target material";
1196  int tgtpdg = 0;
1197  map<int,double>::const_iterator probiter = fCurCumulProbMap.begin();
1198  for( ; probiter != fCurCumulProbMap.end(); ++probiter) {
1199  double prob = probiter->second;
1200  if(R<prob) {
1201  tgtpdg = probiter->first;
1202  LOG("GMCJDriver", pNOTICE)
1203  << "Selected target material = " << tgtpdg;
1204  return tgtpdg;
1205  }
1206  }
1207  LOG("GMCJDriver", pERROR)
1208  << "Could not select target material for an interacting neutrino";
1209  return 0;
1210 }
1211 //___________________________________________________________________________
1213 {
1214  int nupdg = fFluxDriver->PdgCode();
1215  const TLorentzVector & nup4 = fFluxDriver->Momentum();
1216 
1217  // Find the GEVGDriver object that generates interactions for the
1218  // given initial state (neutrino + target)
1219  InitialState init_state(fSelTgtPdg, nupdg);
1220  GEVGDriver * evgdriver = fGPool->FindDriver(init_state);
1221  if(!evgdriver) {
1222  LOG("GMCJDriver", pFATAL)
1223  << "No GEVGDriver object for init state: " << init_state.AsString();
1224  exit(1);
1225  }
1226 
1227  // propagate current unphysical event mask
1228  evgdriver->SetUnphysEventMask(*fUnphysEventMask);
1229 
1230  // Ask the GEVGDriver object to select and generate an interaction for
1231  // the selected initial state & neutrino 4-momentum
1232  LOG("GMCJDriver", pNOTICE)
1233  << "Asking the selected GEVGDriver object to generate an event";
1234  fCurEvt = evgdriver->GenerateEvent(nup4);
1235 }
1236 //___________________________________________________________________________
1238 {
1239  // Generate an 'interaction position' in the selected material, along
1240  // the direction of nup4
1241  LOG("GMCJDriver", pNOTICE)
1242  << "Asking the geometry analyzer to generate a vertex";
1243 
1244  const TLorentzVector & p4 = fFluxDriver->Momentum ();
1245  const TLorentzVector & x4 = fFluxDriver->Position ();
1246 
1247  const TVector3 & vtx = fGeomAnalyzer->GenerateVertex(x4, p4, fSelTgtPdg);
1248 
1249  TVector3 origin(x4.X(), x4.Y(), x4.Z());
1250  origin-=vtx; // computes vector dr = origin - vtx
1251 
1252  double dL = origin.Mag();
1254  double dt = dL/c;
1255 
1256  LOG("GMCJDriver", pNOTICE)
1257  << "|vtx - origin|: dL = " << dL << " m, dt = " << dt << " sec";
1258 
1259  fCurVtx.SetXYZT(vtx.x(), vtx.y(), vtx.z(), x4.T() + dt);
1260 
1261  fCurEvt->SetVertex(fCurVtx);
1262 }
1263 //___________________________________________________________________________
1265 {
1266 // Compute event probability for the given flux neutrino & detector geometry
1267 
1268  // get interaction cross section
1269  double xsec = fCurEvt->XSec();
1270 
1271  // get path length in detector along v direction for specified target material
1272  PathLengthList::const_iterator pliter = fCurPathLengths.find(fSelTgtPdg);
1273  double path_length = pliter->second;
1274 
1275  // get target material mass number
1276  int A = pdg::IonPdgCodeToA(fSelTgtPdg);
1277 
1278  // calculate interaction probability
1279  double P = this->InteractionProbability(xsec, path_length, A);
1280 
1281  //
1282  // get weight for selected event
1283  //
1284 
1285  GHepParticle * nu = fCurEvt->Probe();
1286  int nu_pdg = nu->Pdg();
1287  double Ev = nu->P4()->Energy();
1288 
1289  double weight = 1.0;
1290  if(!fGenerateUnweighted) {
1291  map<int,TH1D*>::const_iterator pmax_iter = fPmax.find(nu_pdg);
1292  assert(pmax_iter != fPmax.end());
1293  TH1D * pmax_hst = pmax_iter->second;
1294  assert(pmax_hst);
1295  double pmax = pmax_hst->GetBinContent(pmax_hst->FindBin(Ev));
1296  assert(pmax>0);
1297  weight = pmax/fGlobPmax;
1298  }
1299 
1300  // set probability & update weight
1301  fCurEvt->SetProbability(P);
1302  fCurEvt->SetWeight(weight * fCurEvt->Weight());
1303 }
1304 //___________________________________________________________________________
1305 double GMCJDriver::InteractionProbability(double xsec, double pL, int A)
1306 {
1307 // P = Na (Avogadro number, atoms/mole) *
1308 // 1/A (1/mass number, mole/gr) *
1309 // xsec (total interaction cross section, cm^2) *
1310 // pL (density-weighted path-length, gr/cm^2)
1311 //
1312  xsec = xsec / units::cm2;
1314 
1315  return kNA*(xsec*pL)/A;
1316 }
1317 //___________________________________________________________________________
1319 {
1320 // Return the pre-computed interaction probability for the current flux
1321 // neutrino index (entry number in flux file). Exit if not possible as
1322 // using meaningless interaction probability leads to incorrect physics
1323 //
1324  if(!fFluxIntTree){
1325  LOG("GMCJDriver", pERROR) <<
1326  "Cannot get pre-computed flux interaction probability as no tree!";
1327  exit(1);
1328  }
1329 
1330  assert(fFluxDriver->Index() >= 0); // Check trying to find meaningfull index
1331 
1332  // Check if can find relevant entry and no mismatch in energies -->
1333  // using correct pre-gen interaction prob file
1334  bool found_entry = fFluxIntTree->GetEntryWithIndex(fFluxDriver->Index()) > 0;
1335  bool enu_match = false;
1336  if(found_entry){
1337  double rel_err = fBrFluxEnu-fFluxDriver->Momentum().E();
1338  if(fBrFluxEnu > controls::kASmallNum) rel_err /= fBrFluxEnu;
1339  enu_match = TMath::Abs(rel_err)<controls::kASmallNum;
1340  if(enu_match == false){
1341  LOG("GMCJDriver", pERROR) <<
1342  "Mismatch between: Enu_curr = "<< fFluxDriver->Momentum().E() <<
1343  ", Enu_pre_gen = "<< fBrFluxEnu;
1344  }
1345  }
1346  else {
1347  LOG("GMCJDriver", pERROR) << "Cannot find flux entry in interaction prob tree!";
1348  }
1349 
1350  // Exit if not successful
1351  bool success = found_entry && enu_match;
1352  if(!success){
1353  LOG("GMCJDriver", pFATAL) <<
1354  "Cannot find pre-generated interaction probability! Check you "<<
1355  "are using the correct pre-generated interaction prob file " <<
1356  "generated using current flux input file with same input " <<
1357  "config (same geom TopVol, neutrino species list)";
1358  exit(1);
1359  }
1360  assert(fGlobPmax+controls::kASmallNum>0.0);
1361  return fBrFluxIntProb/fGlobPmax;
1362 }
1363 //___________________________________________________________________________
static constexpr double dL
const int pmax
Definition: cellShifts.C:13
static const double second
Definition: Units.h:35
void CreateXSecSumSpline(int nk, double Emin, double Emax, bool inlogE=true)
Definition: GEVGDriver.cxx:441
void SetUnphysEventMask(const TBits &mask)
Definition: GMCJDriver.cxx:165
Basic constants.
void PopulateEventGenDriverPool(void)
Definition: GMCJDriver.cxx:605
bool PreCalcFluxProbabilities(void)
Definition: GMCJDriver.cxx:239
void InitJob(void)
Definition: GMCJDriver.cxx:494
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
void GetMaxPathLengthList(void)
Definition: GMCJDriver.cxx:577
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
void SetEventGeneratorList(string listname)
Definition: GMCJDriver.cxx:157
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
const Var weight
void InitEventGeneration(void)
Definition: GMCJDriver.cxx:830
A simple [min,max] interval for doubles.
Definition: Range1.h:43
string BoolAsYNString(bool b)
Definition: PrintUtils.cxx:115
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:47
#define pFATAL
Definition: Messenger.h:57
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:61
bool GenerateFluxNeutrino(void)
bool ComputePathLengths(void)
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:34
string filename
Definition: shutoffs.py:106
const Spline * XSecSumSpline(void) const
Definition: GEVGDriver.h:88
double Evaluate(double x) const
Definition: Spline.cxx:362
static const double cm2
Definition: Units.h:77
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
void ComputeEventProbability(void)
double dE
const double emin
int Pdg(void) const
Definition: GHepParticle.h:64
Range1D_t ValidEnergyRange(void) const
Definition: GEVGDriver.cxx:669
void UseFluxDriver(GFluxI *flux)
Definition: GMCJDriver.cxx:174
#define P(a, b, c, d, e, x)
Object to be filled with the neutrino path-length, for all detector geometry materials, when starting from a position x and travelling along the direction of the neutrino 4-momentum.
void KeepOnThrowingFluxNeutrinos(bool keep_on)
Definition: GMCJDriver.cxx:211
Definition: Cand.cxx:23
void GetMaxFluxEnergy(void)
Definition: GMCJDriver.cxx:595
string outfilename
knobs that need extra care
static unsigned int NFlags(void)
Definition: GHepFlags.h:77
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void SetEventGeneratorList(string listname)
Definition: GEVGDriver.cxx:349
const double emax
void ForceSingleProbScale(void)
Definition: GMCJDriver.cxx:219
double PreGenFluxInteractionProbability(void)
static Messenger * Instance(void)
Definition: Messenger.cxx:71
bool UseMaxPathLengths(string xml_filename)
Definition: GMCJDriver.cxx:190
void Configure(bool calc_prob_scales=true)
Definition: GMCJDriver.cxx:446
static const double kNA
Definition: Constants.h:50
static const double meter
Definition: Units.h:33
string AsString(void) const
#define R(x)
static const double kilogram
Definition: Units.h:34
static const double kASmallNum
Definition: Controls.h:40
def success(message)
Definition: log.py:5
TRandom3 & RndEvg(void) const
rnd number generator used by the event generation drivers
Definition: RandomGen.h:75
GENIE Event Generation Driver. A minimalist user interface object for generating neutrino interaction...
Definition: GEVGDriver.h:55
#define pINFO
Definition: Messenger.h:63
void BootstrapXSecSplines(void)
Definition: GMCJDriver.cxx:642
double ComputeInteractionProbabilities(bool use_max_path_length)
Double_t xsec[nknots]
Definition: testXsec.C:47
#define pWARN
Definition: Messenger.h:61
static const double kLightSpeed
Definition: Constants.h:32
EventRecord * GenerateEvent(void)
Definition: GMCJDriver.cxx:838
static const double gram
Definition: Units.h:141
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
double InteractionProbability(double xsec, double pl, int A)
static const double A
Definition: Units.h:82
double max
Definition: Range1.h:54
void PreSelectEvents(bool preselect=true)
Definition: GMCJDriver.cxx:231
void GenerateEventKinematics(void)
int SelectTargetMaterial(double R)
void GenerateVertexPosition(void)
static const double m2
Definition: Units.h:80
void ComputeProbScales(void)
Definition: GMCJDriver.cxx:705
exit(0)
void Configure(int nu_pdgc, int Z, int A)
Definition: GEVGDriver.cxx:138
bool LoadFluxProbabilities(string filename)
Definition: GMCJDriver.cxx:390
EventRecord * GenerateEvent1Try(void)
Definition: GMCJDriver.cxx:867
assert(nhit_max >=nhit_nbins)
string PrintFramedMesg(string mesg, unsigned int nl=1, const char f='*')
Definition: PrintUtils.cxx:171
void CreateSplines(int nknots=-1, double emax=-1, bool inLogE=true)
Definition: GEVGDriver.cxx:578
double min
Definition: Range1.h:53
void BootstrapXSecSplineSummation(void)
Definition: GMCJDriver.cxx:669
void UseSplines(void)
Definition: GEVGDriver.cxx:509
#define pNOTICE
Definition: Messenger.h:62
Defines the GENIE Geometry Analyzer Interface.
Definition: GeomAnalyzerI.h:30
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:64
void SetUnphysEventMask(const TBits &mask)
Definition: GEVGDriver.cxx:220
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
bool gAbortingInErr
Definition: Messenger.cxx:56
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
void UseGeomAnalyzer(GeomAnalyzerI *geom)
Definition: GMCJDriver.cxx:179
void UseSplines(bool useLogE=true)
Definition: GMCJDriver.cxx:184
void SaveFluxProbabilities(string outfilename)
Definition: GMCJDriver.cxx:437
void GetParticleLists(void)
Definition: GMCJDriver.cxx:560
static AlgConfigPool * Instance()
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
EventRecord * GenerateEvent(const TLorentzVector &nu4p)
Definition: GEVGDriver.cxx:229
A pool of GEVGDriver objects with an initial state key.
Definition: GEVGPool.h:38