GJPARCNuFlux.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 - Feb 04, 2008
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Feb 04, 2008 - CA
14  The first implementation of this concrete flux driver was first added in
15  the development version 2.3.1
16  @ Feb 19, 2008 - CA
17  Extended to handle all near detector locations and super-k
18  @ Feb 22, 2008 - CA
19  Added method to report the actual POT.
20  @ Mar 05, 2008 - CA,JD
21  Added method to configure the starting z position (upstream of the detector
22  face, in detector coord system). Added code to project flux neutrinos from
23  z=0 to a configurable z position (somewhere upstream of the detector face)
24  @ Mar 27, 2008 - CA
25  Added option to recycle the flux ntuple
26  @ Mar 31, 2008 - CA
27  Handle the flux ntuple weights. Renamed the old implementation of GFluxI
28  GenerateNext() to GenerateNext_weighted(). Coded-up a new GenerateNext()
29  method using GenerateNext_weighted() + the rejection method to generate
30  unweighted flux neutrinos. Added code in LoadBeamSimData() to scan for the
31  max. flux weight for the input location. The 'actual POT' norm factor is
32  updated after each generated flux neutrino taking into account the flux
33  weight variability. Added NFluxNeutrinos() and SumWeight().
34  @ May 29, 2008 - CA, PG
35  Protect LoadBeamSimData() against non-existent input file
36  @ May 31, 2008 - CA
37  Added option to keep on recyclying the flux ntuple for an 'infinite' number
38  of times (SetNumOfCycles(0)) so that exiting the event generation loop can
39  be controlled by the accumulated POTs or number of events.
40  @ June 1, 2008 - CA
41  At LoadBeamSimData() added code to scan for number of neutrinos and sum of
42  weights for a complete flux ntuple cycle, at the specified detector.
43  Added POT_1cycle() to return the flux POT per flux ntuple cycle.
44  Added POT_curravg() to return the current average number of POTs taking
45  into account the flux ntuple recycling. At complete cycles the reported
46  number is not just the average POT but the exact POT.
47  Removed the older ActualPOT() method.
48  @ June 4, 2008 - CA, AM
49  Small modification at the POT normalization to account for the fact that
50  flux neutrinos get de-weighted before passed to the event generation code.
51  Now pot = file_pot/max_wght rather than file_pot*(nneutrinos/sum_wght)
52  @ June 19, 2008 - CA
53  Removing some LOG() mesgs speeds up GenerateNext() by a factor of 20 (!)
54  @ June 18, 2009 - CA
55  Demote warning mesgs if the current flux neutrino is skipped because it
56  is not in the list of neutrinos to be considered.
57  Now this can be intentional (eg. if generating nu_e only).
58  In GenerateNext_weighted() Moved the code updating the number of neutrinos
59  and sum of weights higher so that force-rejected flux neutrinos still
60  count for normalization purposes.
61  @ March 6, 2010 - JD
62  Made compatible with 10a version of jnubeam flux. Maintained compatibility
63  with 07a version. This involved finding a few more branches - now look for
64  all branches and only use them if they exist. Also added check for required
65  branches. GJPARCNuFluxPassThroughInfo class now passes through flux info in
66  identical format as given in flux file. Any conversions to more usable
67  format happen at later stage (gNtpConv.cxx). In addition also store the flux
68  entry number for current neutrino and the deduced flux version.
69  Changed method to search for maximum flux weight to avoid seg faults when
70  have large number of flux entries in a file (~1.5E6).
71  @ March 8, 2010 - JD
72  Incremented the GJPARCNuFluxPassThroughInfo class def number to reflect
73  changes made in previous commit. Added a failsafe which will terminate
74  the job if go through a whole flux cycle without finding a detector location
75  matching that specified by user. This avoids potential infinite number of
76  cycles.
77  @ Feb 4, 2011 - JD
78  Made compatible with 11a version of jnubeam flux. Still compatable with older
79  versions. Now set the data members of the pass-through class directly as the
80  branch addresses of the input tree to reduce amount of duplicated code.
81  @ Feb 22, 2011 - JD
82  Added functionality to start looping over input flux file from a random
83  offset. This is to avoid any potential biases when processing very large
84  flux files and always starting from the same position. The default is to
85  apply a random offset but this can be switched off using the
86  GJPARCNuFlux::DisableOffset() method.
87  @ Feb 22, 2011 - JD
88  Implemented the new GFluxI::Clear, GFluxI::Index and GFluxI::GenerateWeighted
89  methods needed so that can be used with the new pre-generation of flux
90  interaction probabilities methods added to GMCJDriver.
91  @ Feb 24, 2011 - JD
92  Updated list of expected decay modes for the JNuBeam flux neutrinos for >10a
93  flux mc. The decay mode is used to infer the neutrino pdg and previously we
94  were just skipping them if we didn't recognise the mode - now the job aborts
95  as this can lead to unphysical results.
96  @ Feb 26, 2011 - JD
97  Now check that there is at least one entry with matching flux location (idfd)
98  at the LoadBeamSimData stage. Previously were only checking this after a
99  cycle of calling GenerateNext. This stops case where if only looping over a
100  single cycle the user was not warned that there were no flux location matches.
101  @ Jul 05, 2011 - TD
102  Code used to match nd detector locations 1-10. Now match detector locations
103  up to 51. Change made in preparation for the new sand muon flux (nd13).
104  @ Feb 09, 2012 - TD
105  Added ability to TChain flux files together. This is so when doing vector
106  producion, we can sample all the input flux files, even if the equivalent
107  hadd'ed flux file is too large.
108  @ Mar 14, 2014 - TD
109  Prevent an infinite loop in GenerateNext() when the flux driver has not been
110  properly configured by exiting within GenerateNext_weighted().
111  LoadBeamSimData() now returns bool, so that the user can catch cases when
112  the flux driver has not been properly configured.
113 */
114 //____________________________________________________________________________
115 
116 #include <cstdlib>
117 #include <iostream>
118 #include <cassert>
119 
120 #include <TFile.h>
121 #include <TTree.h>
122 #include <TChain.h>
123 #include <TString.h>
124 #include <TSystem.h>
125 
129 #include "Tools/Flux/GJPARCNuFlux.h"
137 
140 
141 using std::endl;
142 using std::cout;
143 using std::endl;
144 using namespace genie;
145 using namespace genie::flux;
146 
148 
149 //____________________________________________________________________________
151 {
152  this->Initialize();
153 }
154 //___________________________________________________________________________
155 GJPARCNuFlux::~GJPARCNuFlux()
156 {
157  this->CleanUp();
158 }
159 //___________________________________________________________________________
161 {
162 // Get next (unweighted) flux ntuple entry on the specified detector location
163 //
165  while(1) {
166  // Check for end of flux ntuple
167  bool end = this->End();
168  if(end) return false;
169 
170  // Get next weighted flux ntuple entry
171  bool nextok = this->GenerateNext_weighted();
172  if(!nextok) continue;
173 
174  if(fNCycles==0) {
175  LOG("Flux", pNOTICE)
176  << "Got flux entry: " << this->Index()
177  << " - Cycle: "<< fICycle << "/ infinite";
178  } else {
179  LOG("Flux", pNOTICE)
180  << "Got flux entry: "<< this->Index()
181  << " - Cycle: "<< fICycle << "/"<< fNCycles;
182  }
183 
184  // If de-weighting get fractional weight & decide whether to accept curr flux neutrino
185  double f = 1.0;
186  if(fGenerateWeighted == false) f = this->Weight();
187  LOG("Flux", pNOTICE)
188  << "Curr flux neutrino fractional weight = " << f;
189  if(f > (1.+controls::kASmallNum)) {
190  LOG("Flux", pERROR)
191  << "** Fractional weight = " << f << " > 1 !!";
192  }
193  double r = (f < 1.) ? rnd->RndFlux().Rndm() : 0;
194  bool accept = (r<f);
195  if(accept) {
196  return true;
197  }
198 
199  LOG("Flux", pNOTICE)
200  << "** Rejecting current flux neutrino based on the flux weight only";
201  }
202  return false;
203 }
204 //___________________________________________________________________________
206 {
207 // Get next (weighted) flux ntuple entry on the specified detector location
208 //
209 
210  // Reset previously generated neutrino code / 4-p / 4-x
211  this->ResetCurrent();
212 
213  // Check whether a jnubeam flux ntuple has been loaded
215  LOG("Flux", pFATAL)
216  << "The flux driver has not been properly configured";
217  //return false; // don't do this - creates an infinite loop!
218  exit(1);
219  }
220 
221  // Read next flux ntuple entry. Use fEntriesThisCycle to keep track of when
222  // in new cycle as fIEntry can now have an offset
224  // Exit if have not found neutrino at specified location for whole cycle
225  if(fNDetLocIdFound == 0){
226  LOG("Flux", pFATAL)
227  << "The input jnubeam flux ntuple contains no entries for detector id "
228  << fDetLocId << ". Terminating job!";
229  exit(1);
230  }
231  fNDetLocIdFound = 0; // reset the counter
232  fICycle++;
234  fEntriesThisCycle = 0;
235  // Run out of entries @ the current cycle.
236  // Check whether more (or infinite) number of cycles is requested
237  if(fICycle >= fNCycles && fNCycles != 0){
238  LOG("Flux", pWARN)
239  << "No more entries in input flux neutrino ntuple";
240  return false;
241  }
242  }
243 
244  // In addition to getting info to generate event the following also
245  // updates pass-through info (= info on the flux neutrino parent particle
246  // that may be stored at an extra branch of the output event tree -alongside
247  // with the generated event branch- for use further upstream in the t2k
248  // analysis chain -eg for beam reweighting etc-)
249  bool found_entry;
250  if (fNuFluxUsingTree)
251  found_entry = fNuFluxTree->GetEntry(fIEntry) > 0;
252  else
253  found_entry = fNuFluxChain->GetEntry(fIEntry) > 0;
254  assert(found_entry);
255  fLoadedNeutrino = true;
257  fIEntry = (fIEntry+1) % fNEntries;
258 
259  if (fNuFluxUsingTree) {
260  if(fNuFluxSumTree) fNuFluxSumTree->GetEntry(0); // get entry 0 as only 1 entry in tree
261  }
262  else {
263  // get entry corresponding to current tree number in the chain, as only 1 entry in each tree
264  if(fNuFluxSumChain) fNuFluxSumChain->GetEntry(fNuFluxChain->GetTreeNumber());
265  }
266 
267  // check for negative flux weights
269  LOG("Flux", pERROR) << "Negative flux weight! Will set weight to 0.0";
270  fPassThroughInfo->norm = 0.0;
271  }
272  // remember to update fNorm as no longer set to fNuFluxTree branch address
273  fNorm = (double) fPassThroughInfo->norm;
274 
275  // for 'near detector' flux ntuples make sure that the current entry
276  // corresponds to a flux neutrino at the specified detector location
277  if(fIsNDLoc /* nd */ &&
278  fDetLocId!=fPassThroughInfo->idfd /* doesn't match specified detector location*/) {
279 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
280  LOG("Flux", pNOTICE)
281  << "Current flux neutrino not at specified detector location";
282 #endif
283  return false;
284  }
285 
286 
287  //
288  // Handling neutrinos at specified detector location
289  //
290 
291  // count the number of times we have neutrinos at specified detector location
292  fNDetLocIdFound += 1;
293 
294 
295  // update the sum of weights & number of neutrinos
296  fSumWeight += this->Weight() * fMaxWeight; // Weight returns fNorm/fMaxWeight
297  fNNeutrinos++;
298 
299  // Convert the current parent paticle decay mode into a neutrino pdg code
300  // See:
301  // http://jnusrv01.kek.jp/internal/t2k/nubeam/flux/nemode.h
302  // mode description
303  // 11 numu from pi+
304  // 12 numu from K+
305  // 13 numu from mu-
306  // 21 numu_bar from pi-
307  // 22 numu_bar from K-
308  // 23 numu_bar from mu+
309  // 31 nue from K+ (Ke3)
310  // 32 nue from K0L(Ke3)
311  // 33 nue from Mu+
312  // 41 nue_bar from K- (Ke3)
313  // 42 nue_bar from K0L(Ke3)
314  // 43 nue_bar from Mu-
315  // Since JNuBeam flux version >= 10a the following modes also expected
316  // 14 numu from K+(3)
317  // 15 numu from K0(3)
318  // 24 numu_bar from K-(3)
319  // 25 numu_bar from K0(3)
320  // 34 nue from pi+
321  // 44 nue_bar from pi-
322  // In general expect more modes following the rule:
323  // 11->19 --> numu
324  // 21->29 --> numu_bar
325  // 31->39 --> nue
326  // 41->49 --> nuebar
327  // This is based on example given at:
328  // http://jnusrv01.kek.jp/internal/t2k/nubeam/flux/efill.kumac
329 
330  if(fPassThroughInfo->mode >= 11 && fPassThroughInfo->mode <= 19) fgPdgC = kPdgNuMu;
331  else if(fPassThroughInfo->mode >= 21 && fPassThroughInfo->mode <= 29) fgPdgC = kPdgAntiNuMu;
332  else if(fPassThroughInfo->mode >= 31 && fPassThroughInfo->mode <= 39) fgPdgC = kPdgNuE;
333  else if(fPassThroughInfo->mode >= 41 && fPassThroughInfo->mode <= 49) fgPdgC = kPdgAntiNuE;
334  else {
335  // If here then trying to process a neutrino from an unknown decay mode.
336  // Rather than just skipping this flux neutrino the job is aborted to avoid
337  // unphysical results.
338  LOG("Flux", pFATAL) << "Unexpected decay mode: "<< fPassThroughInfo->mode <<
339  " --> unable to infer neutrino pdg! Aborting job!";
340  exit(1);
341  }
342 
343  // Check neutrino pdg against declared list of neutrino species declared
344  // by the current instance of the JPARC neutrino flux driver.
345  // No undeclared neutrino species will be accepted at this point as GENIE
346  // has already been configured to handle the specified list.
347  // Make sure that the appropriate list of flux neutrino species was set at
348  // initialization via GJPARCNuFlux::SetFluxParticles(const PDGCodeList &)
349 
351 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
352  LOG("Flux", pNOTICE)
353  << "Current flux neutrino "
354  << "not at the list of neutrinos to be considered at this job.";
355 #endif
356  return false;
357  }
358 
359  // Check current neutrino energy against the maximum flux neutrino energy declared
360  // by the current instance of the JPARC neutrino flux driver.
361  // No flux neutrino exceeding that maximum energy will be accepted at this point as
362  // that maximum energy has already been used for normalizing the interaction probabilities.
363  // Make sure that the appropriate maximum flux neutrino energy was set at
364  // initialization via GJPARCNuFlux::SetMaxEnergy(double Ev)
365 
366  if(fPassThroughInfo->Enu > fMaxEv) {
367  LOG("Flux", pWARN)
368  << "Flux neutrino energy exceeds declared maximum neutrino energy";
369  LOG("Flux", pWARN)
370  << "Ev = " << fPassThroughInfo->Enu << "(> Ev{max} = " << fMaxEv << ")";
371  }
372 
373  // Set the current flux neutrino 4-momentum & 4-position
374 
375  double pxnu = fPassThroughInfo->Enu * fPassThroughInfo->nnu[0];
376  double pynu = fPassThroughInfo->Enu * fPassThroughInfo->nnu[1];
377  double pznu = fPassThroughInfo->Enu * fPassThroughInfo->nnu[2];
378  double Enu = fPassThroughInfo->Enu;
379  fgP4.SetPxPyPzE (pxnu, pynu, pznu, Enu);
380 
381  if(fIsNDLoc) {
382  double cm2m = units::cm / units::m;
383  double xnu = cm2m * fPassThroughInfo->xnu;
384  double ynu = cm2m * fPassThroughInfo->ynu;
385  double znu = 0;
386 
387  // projected 4-position (from z=0) back to a configurable plane
388  // position (fZ0) upstream of the detector face.
389  xnu += (fZ0/fPassThroughInfo->nnu[2])*fPassThroughInfo->nnu[0];
390  ynu += (fZ0/fPassThroughInfo->nnu[2])*fPassThroughInfo->nnu[1];
391  znu = fZ0;
392 
393  fgX4.SetXYZT (xnu, ynu, znu, 0.);
394  } else {
395  fgX4.SetXYZT (0.,0.,0.,0.);
396 
397  }
398 
399 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
400  LOG("Flux", pINFO)
401  << "Generated neutrino: "
402  << "\n pdg-code: " << fgPdgC
403  << "\n p4: " << utils::print::P4AsShortString(&fgP4)
404  << "\n x4: " << utils::print::X4AsString(&fgX4);
405 #endif
406  // Update flux pass through info not set as branch addresses of flux ntuples
407  if(fNuFluxUsingTree)
408  fPassThroughInfo->fluxentry = this->Index();
409  else
410  fPassThroughInfo->fluxentry = fNuFluxChain->GetTree()->GetReadEntry();
411 
413  if(fNuFluxUsingTree)
414  filename = fNuFluxFile->GetName();
415  else
416  filename = fNuFluxChain->GetFile()->GetName();
417  std::string::size_type start_pos = filename.rfind("/");
418  if (start_pos == std::string::npos) start_pos = 0; else ++start_pos;
419  std::string basename(filename,start_pos);
420  if(fNuFluxUsingTree)
421  fPassThroughInfo->fluxfilename = basename + ":" + fNuFluxTree->GetName();
422  else
423  fPassThroughInfo->fluxfilename = basename + ":" + fNuFluxChain->GetName();
424  return true;
425 }
426 //___________________________________________________________________________
428 {
429 // Compute number of flux POTs / flux ntuple cycle
430 //
432  LOG("Flux", pWARN)
433  << "The flux driver has not been properly configured";
434  return 0;
435  }
436 
437 // double pot = fFilePOT * (fNNeutrinosTot1c/fSumWeightTot1c);
438 //
439 // Use the max weight instead, since flux neutrinos get de-weighted
440 // before thrown to the event generation driver
441 //
442  double pot = fFilePOT / fMaxWeight;
443  return pot;
444 }
445 //___________________________________________________________________________
447 {
448 // Compute current number of flux POTs
449 // On complete cycles, that POT number should be exact.
450 // Within cycles that is only an average number
451 
453  LOG("Flux", pWARN)
454  << "The flux driver has not been properly configured";
455  return 0;
456  }
457 
458 // double pot = fNNeutrinos*fFilePOT/fSumWeightTot1c;
459 //
460 // See also comment at POT_1cycle()
461 //
462  double cnt = (double)fNNeutrinos;
463  double cnt1c = (double)fNNeutrinosTot1c;
464  double pot = (cnt/cnt1c) * this->POT_1cycle();
465  return pot;
466 }
467 //___________________________________________________________________________
468 long int GJPARCNuFlux::Index(void)
469 {
470 // Return the current flux entry index. If GenerateNext has not yet been
471 // called then return -1.
472 //
473  if(fLoadedNeutrino){
474  // subtract 1 as fIEntry was incremented since call to TTree::GetEntry
475  // and deal with special case where fIEntry-1 is last entry in cycle
476  return fIEntry == 0 ? fNEntries - 1 : fIEntry-1;
477  }
478  // return -1 if no neutrino loaded since last call to this->ResetCurrent()
479  return -1;
480 }
481 //___________________________________________________________________________
482 bool GJPARCNuFlux::LoadBeamSimData(string filename, string detector_location)
483 {
484 // Loads in a jnubeam beam simulation root file (converted from hbook format)
485 // into the GJPARCNuFlux driver.
486 // The detector location can be any of:
487 // "sk","nd1" (<-2km),"nd5" (<-nd280),...,"nd10"
488 
489  LOG("Flux", pNOTICE)
490  << "Loading jnubeam flux tree from ROOT file: " << filename;
491  LOG("Flux", pNOTICE)
492  << "Detector location: " << detector_location;
493 
494  // Check to see if its a single flux file (/dir/root_filename.0.root)
495  // or a sequence of files to be tchained (e.g. /dir/root_filename@0@100)
496  fNuFluxUsingTree = true;
497  if (filename.find('@') != string::npos) { fNuFluxUsingTree = false; }
498 
499  vector<string> filenamev = utils::str::Split(filename,"@");
500  string fileroot = "";
501  int firstfile = -1, lastfile = -1;
502 
503  if (!fNuFluxUsingTree) {
504  if (filenamev.size() != 3) {
505  LOG("Flux", pFATAL)
506  << "Flux filename should be specfied as either:\n"
507  << "\t For a single input file: /dir/root_filename.#.root\n"
508  << "\t For multiple input files: /dir/root_filename@#@#";
509  exit(1);
510  }
511  fileroot = filenamev[0];
512  firstfile = atoi(filenamev[1].c_str());
513  lastfile = atoi(filenamev[2].c_str());
514  LOG("Flux", pNOTICE)
515  << "Chaining beam simulation output files with stem: " << fileroot
516  << " and run numbers in the range: [" << firstfile << ", " << firstfile << "]";
517  }
518 
519  if (fNuFluxUsingTree) {
520  bool is_accessible = ! (gSystem->AccessPathName( filename.c_str() ));
521  if (!is_accessible) {
522  LOG("Flux", pFATAL)
523  << "The input jnubeam flux file doesn't exist! Initialization failed!";
524  exit(1);
525  }
526  }
527  else {
528  bool please_exit = false;
529  for (int i = firstfile; i < lastfile+1; i++) {
530  bool is_accessible = ! (gSystem->AccessPathName( Form("%s.%i.root",fileroot.c_str(),i) ));
531  if (!is_accessible) {
532  LOG("Flux", pFATAL)
533  << "The input jnubeam flux file " << Form("%s.%i.root",fileroot.c_str(),i)
534  << "doesn't exist! Initialization failed!";
535  please_exit = true;
536  }
537  }
538  if(please_exit)
539  exit(1);
540  }
541 
542  fDetLoc = detector_location;
543  fDetLocId = this->DLocName2Id(fDetLoc);
544 
545  if(fDetLocId == 0) {
546  LOG("Flux", pERROR)
547  << " ** Unknown input detector location: " << fDetLoc;
548  return false;
549  }
550 
551  fIsFDLoc = (fDetLocId==-1);
552  fIsNDLoc = (fDetLocId>0);
553 
554  if (!fNuFluxUsingTree) {
555  string ntuple_name = (fIsNDLoc) ? "h3002" : "h2000";
556  fNuFluxChain = new TChain(ntuple_name.c_str());
557  int result = fNuFluxChain->Add( Form("%s.%i.root",fileroot.c_str(),firstfile), 0);
558  if (result != 1 && fIsNDLoc) {
559  LOG("Flux", pINFO)
560  << "Could not find tree h3002 in file " << Form("%s.%i.root",fileroot.c_str(),firstfile)
561  << " Trying tree h3001";
562  delete fNuFluxChain;
563  ntuple_name = "h3001";
564  fNuFluxChain = new TChain(ntuple_name.c_str());
565  result = fNuFluxChain->Add( Form("%s.%i.root",fileroot.c_str(),firstfile), 0);
566  }
567  if (result != 1) {
568  LOG("Flux", pERROR)
569  << "** Couldn't get flux tree: " << ntuple_name;
570  return false;
571  }
572 
573  for (int i = firstfile+1; i < lastfile+1; i++) {
574  result = fNuFluxChain->Add( Form("%s.%i.root",fileroot.c_str(),i), 0 );
575  if (result == 0)
576  LOG("Flux", pERROR)
577  << "** Couldn't get flux tree " << ntuple_name << " in file " << Form("%s.%i.root",fileroot.c_str(),i);
578  fNEntries = fNuFluxChain->GetEntries();
579  }
580  }
581 
582  else {
583  fNuFluxFile = new TFile(filename.c_str(), "read");
584  if(fNuFluxFile) {
585  // nd treename can be h3002 or h3001 depending on fluxfile version
586  string ntuple_name = (fIsNDLoc) ? "h3002" : "h2000";
587  fNuFluxTree = (TTree*) fNuFluxFile->Get(ntuple_name.c_str());
588  if(!fNuFluxTree && fIsNDLoc){
589  ntuple_name = "h3001";
590  fNuFluxTree = (TTree*) fNuFluxFile->Get(ntuple_name.c_str());
591  }
592  LOG("Flux", pINFO)
593  << "Getting flux tree: " << ntuple_name;
594  if(!fNuFluxTree) {
595  LOG("Flux", pERROR)
596  << "** Couldn't get flux tree: " << ntuple_name;
597  return false;
598  }
599  } else {
600  LOG("Flux", pERROR) << "** Couldn't open: " << filename;
601  return false;
602  }
603 
604  fNEntries = fNuFluxTree->GetEntries();
605  }
606 
607  LOG("Flux", pNOTICE)
608  << "Loaded flux tree contains " << fNEntries << " entries";
609 
610  LOG("Flux", pDEBUG)
611  << "Getting tree branches & setting leaf addresses";
612 
613  // try to get all the branches that we know about and only set address if
614  // they exist
615  bool missing_critical = false;
616  TBranch * fBr = 0;
618 
619  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("norm")) ) fBr->SetAddress(&info->norm);
620  else if( (fBr = fNuFluxChain->GetBranch("norm")) ) fNuFluxChain->SetBranchAddress("norm",&info->norm);
621  else {
622  LOG("Flux", pFATAL) <<"Cannot find flux branch: norm";
623  missing_critical = true;
624  }
625  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("Enu")) ) fBr->SetAddress(&info->Enu);
626  else if( (fBr = fNuFluxChain->GetBranch("Enu")) ) fNuFluxChain->SetBranchAddress("Enu",&info->Enu);
627  else {
628  LOG("Flux", pFATAL) <<"Cannot find flux branch: Enu";
629  missing_critical = true;
630  }
631  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ppid")) ) fBr->SetAddress(&info->ppid);
632  else if( (fBr = fNuFluxChain->GetBranch("ppid")) ) fNuFluxChain->SetBranchAddress("ppid",&info->ppid);
633  else {
634  LOG("Flux", pFATAL) <<"Cannot find flux branch: ppid";
635  missing_critical = true;
636  }
637  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("mode")) ) fBr->SetAddress(&info->mode);
638  else if( (fBr = fNuFluxChain->GetBranch("mode")) ) fNuFluxChain->SetBranchAddress("mode",&info->mode);
639  else {
640  LOG("Flux", pFATAL) <<"Cannot find flux branch: mode";
641  missing_critical = true;
642  }
643  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("rnu")) ) fBr->SetAddress(&info->rnu);
644  else if( (fBr = fNuFluxChain->GetBranch("rnu")) ) fNuFluxChain->SetBranchAddress("rnu",&info->rnu);
645  else if(fIsNDLoc){ // Only required for ND location
646  LOG("Flux", pFATAL) <<"Cannot find flux branch: rnu";
647  missing_critical = true;
648  }
649  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("xnu")) ) fBr->SetAddress(&info->xnu);
650  else if( (fBr = fNuFluxChain->GetBranch("xnu")) ) fNuFluxChain->SetBranchAddress("xnu",&info->xnu);
651  else if(fIsNDLoc){ // Only required for ND location
652  LOG("Flux", pFATAL) <<"Cannot find flux branch: xnu";
653  missing_critical = true;
654  }
655  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ynu")) ) fBr->SetAddress(&info->ynu);
656  else if( (fBr = fNuFluxChain->GetBranch("ynu")) ) fNuFluxChain->SetBranchAddress("ynu",&info->ynu);
657  else if(fIsNDLoc) { // Only required for ND location
658  LOG("Flux", pFATAL) <<"Cannot find flux branch: ynu";
659  missing_critical = true;
660  }
661  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("nnu")) ) fBr->SetAddress(info->nnu);
662  else if( (fBr = fNuFluxChain->GetBranch("nnu")) ) fNuFluxChain->SetBranchAddress("nnu",info->nnu);
663  else if(fIsNDLoc){ // Only required for ND location
664  LOG("Flux", pFATAL) <<"Cannot find flux branch: nnu";
665  missing_critical = true;
666  }
667  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("idfd")) ) fBr->SetAddress(&info->idfd);
668  else if( (fBr = fNuFluxChain->GetBranch("idfd")) ) fNuFluxChain->SetBranchAddress("idfd",&info->idfd);
669  else if(fIsNDLoc){ // Only required for ND location
670  LOG("Flux", pFATAL) <<"Cannot find flux branch: idfd";
671  missing_critical = true;
672  }
673  // check that have found essential branches
674  if(missing_critical){
675  LOG("Flux", pFATAL)
676  << "Unable to find critical information in the flux ntuple! Initialization failed!";
677  exit(1);
678  }
679  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ppi")) ) fBr->SetAddress(&info->ppi);
680  else if( (fBr = fNuFluxChain->GetBranch("ppi")) ) fNuFluxChain->SetBranchAddress("ppi",&info->ppi);
681  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("xpi")) ) fBr->SetAddress(info->xpi);
682  else if( (fBr = fNuFluxChain->GetBranch("xpi")) ) fNuFluxChain->SetBranchAddress("xpi",info->xpi);
683  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("npi")) ) fBr->SetAddress(info->npi);
684  else if( (fBr = fNuFluxChain->GetBranch("npi")) ) fNuFluxChain->SetBranchAddress("npi",info->npi);
685  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ppi0")) ) fBr->SetAddress(&info->ppi0);
686  else if( (fBr = fNuFluxChain->GetBranch("ppi0")) ) fNuFluxChain->SetBranchAddress("ppi0",&info->ppi0);
687  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("xpi0")) ) fBr->SetAddress(info->xpi0);
688  else if( (fBr = fNuFluxChain->GetBranch("xpi0")) ) fNuFluxChain->SetBranchAddress("xpi0",info->xpi0);
689  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("npi0")) ) fBr->SetAddress(info->npi0);
690  else if( (fBr = fNuFluxChain->GetBranch("npi0")) ) fNuFluxChain->SetBranchAddress("npi0",info->npi0);
691  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("nvtx0")) ) fBr->SetAddress(&info->nvtx0);
692  else if( (fBr = fNuFluxChain->GetBranch("nvtx0")) ) fNuFluxChain->SetBranchAddress("nvtx0",&info->nvtx0);
693  // Following branches only present since flux version 10a
694  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("cospibm")) ) fBr->SetAddress(&info->cospibm);
695  else if( (fBr = fNuFluxChain->GetBranch("cospibm")) ) fNuFluxChain->SetBranchAddress("cospibm",&info->cospibm);
696  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("cospi0bm"))) fBr->SetAddress(&info->cospi0bm);
697  else if( (fBr = fNuFluxChain->GetBranch("cospi0bm"))) fNuFluxChain->SetBranchAddress("cospi0bm",&info->cospi0bm);
698  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gamom0")) ) fBr->SetAddress(&info->gamom0);
699  else if( (fBr = fNuFluxChain->GetBranch("gamom0")) ) fNuFluxChain->SetBranchAddress("gamom0",&info->gamom0);
700  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gipart")) ) fBr->SetAddress(&info->gipart);
701  else if( (fBr = fNuFluxChain->GetBranch("gipart")) ) fNuFluxChain->SetBranchAddress("gipart",&info->gipart);
702  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gvec0")) ) fBr->SetAddress(info->gvec0);
703  else if( (fBr = fNuFluxChain->GetBranch("gvec0")) ) fNuFluxChain->SetBranchAddress("gvec0",info->gvec0);
704  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpos0")) ) fBr->SetAddress(info->gpos0);
705  else if( (fBr = fNuFluxChain->GetBranch("gpos0")) ) fNuFluxChain->SetBranchAddress("gpos0",info->gpos0);
706  // Following branches only present since flux vesion 10d
707  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ng")) ) fBr->SetAddress(&info->ng);
708  else if( (fBr = fNuFluxChain->GetBranch("ng")) ) fNuFluxChain->SetBranchAddress("ng",&info->ng);
709  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpid")) ) fBr->SetAddress(info->gpid);
710  else if( (fBr = fNuFluxChain->GetBranch("gpid")) ) fNuFluxChain->SetBranchAddress("gpid",info->gpid);
711  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gmec")) ) fBr->SetAddress(info->gmec);
712  else if( (fBr = fNuFluxChain->GetBranch("gmec")) ) fNuFluxChain->SetBranchAddress("gmec",info->gmec);
713  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gvx")) ) fBr->SetAddress(info->gvx);
714  else if( (fBr = fNuFluxChain->GetBranch("gvx")) ) fNuFluxChain->SetBranchAddress("gvx",info->gvx);
715  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gvy")) ) fBr->SetAddress(info->gvy);
716  else if( (fBr = fNuFluxChain->GetBranch("gvy")) ) fNuFluxChain->SetBranchAddress("gvy",info->gvy);
717  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gvz")) ) fBr->SetAddress(info->gvz);
718  else if( (fBr = fNuFluxChain->GetBranch("gvz")) ) fNuFluxChain->SetBranchAddress("gvz",info->gvz);
719  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpx")) ) fBr->SetAddress(info->gpx);
720  else if( (fBr = fNuFluxChain->GetBranch("gpx")) ) fNuFluxChain->SetBranchAddress("gpx",info->gpx);
721  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpy")) ) fBr->SetAddress(info->gpy);
722  else if( (fBr = fNuFluxChain->GetBranch("gpy")) ) fNuFluxChain->SetBranchAddress("gpy",info->gpy);
723  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpz")) ) fBr->SetAddress(info->gpz);
724  else if( (fBr = fNuFluxChain->GetBranch("gpz")) ) fNuFluxChain->SetBranchAddress("gpz",info->gpz);
725  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gmat")) ) fBr->SetAddress(info->gmat);
726  else if( (fBr = fNuFluxChain->GetBranch("gmat")) ) fNuFluxChain->SetBranchAddress("gmat",info->gmat);
727  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gdistc")) ) fBr->SetAddress(info->gdistc);
728  else if( (fBr = fNuFluxChain->GetBranch("gdistc")) ) fNuFluxChain->SetBranchAddress("gdistc",info->gdistc);
729  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gdistal")) ) fBr->SetAddress(&info->gdistal);
730  else if( (fBr = fNuFluxChain->GetBranch("gdistal")) ) fNuFluxChain->SetBranchAddress("gdistal",&info->gdistal);
731  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gdistti")) ) fBr->SetAddress(&info->gdistti);
732  else if( (fBr = fNuFluxChain->GetBranch("gdistti")) ) fNuFluxChain->SetBranchAddress("gdistti",&info->gdistti);
733  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gdistfe")) ) fBr->SetAddress(&info->gdistfe);
734  else if( (fBr = fNuFluxChain->GetBranch("gdistfe")) ) fNuFluxChain->SetBranchAddress("gdistfe",&info->gdistfe);
735  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gcosbm")) ) fBr->SetAddress(info->gcosbm);
736  else if( (fBr = fNuFluxChain->GetBranch("gcosbm")) ) fNuFluxChain->SetBranchAddress("gcosbm",info->gcosbm);
737  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("Enusk")) ) fBr->SetAddress(&info->Enusk);
738  else if( (fBr = fNuFluxChain->GetBranch("Enusk")) ) fNuFluxChain->SetBranchAddress("Enusk",&info->Enusk);
739  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("normsk")) ) fBr->SetAddress(&info->normsk);
740  else if( (fBr = fNuFluxChain->GetBranch("normsk")) ) fNuFluxChain->SetBranchAddress("normsk",&info->normsk);
741  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("anorm")) ) fBr->SetAddress(&info->anorm);
742  else if( (fBr = fNuFluxChain->GetBranch("anorm")) ) fNuFluxChain->SetBranchAddress("anorm",&info->anorm);
743 
744  // Look for the flux file summary info tree (only expected for > 10a flux versions)
745  if( fNuFluxUsingTree && (fNuFluxSumTree = (TTree*) fNuFluxFile->Get("h1000")) ){
746  if( (fBr = fNuFluxSumTree->GetBranch("version"))) fBr->SetAddress(&info->version);
747  if( (fBr = fNuFluxSumTree->GetBranch("ntrig")) ) fBr->SetAddress(&info->ntrig);
748  if( (fBr = fNuFluxSumTree->GetBranch("tuneid")) ) fBr->SetAddress(&info->tuneid);
749  if( (fBr = fNuFluxSumTree->GetBranch("pint")) ) fBr->SetAddress(&info->pint);
750  if( (fBr = fNuFluxSumTree->GetBranch("bpos")) ) fBr->SetAddress(info->bpos);
751  if( (fBr = fNuFluxSumTree->GetBranch("btilt")) ) fBr->SetAddress(info->btilt);
752  if( (fBr = fNuFluxSumTree->GetBranch("brms")) ) fBr->SetAddress(info->brms);
753  if( (fBr = fNuFluxSumTree->GetBranch("emit")) ) fBr->SetAddress(info->emit);
754  if( (fBr = fNuFluxSumTree->GetBranch("alpha")) ) fBr->SetAddress(info->alpha);
755  if( (fBr = fNuFluxSumTree->GetBranch("hcur")) ) fBr->SetAddress(info->hcur);
756  if( (fBr = fNuFluxSumTree->GetBranch("rand")) ) fBr->SetAddress(&info->rand);
757  if( (fBr = fNuFluxSumTree->GetBranch("rseed")) ) fBr->SetAddress(info->rseed);
758  }
759 
760  // Look for the flux file summary info tree (only expected for > 10a flux versions)
761  if ( !fNuFluxUsingTree ) {
762  fNuFluxSumChain = new TChain("h1000");
763  int result = fNuFluxSumChain->Add( Form("%s.%i.root",fileroot.c_str(),firstfile), 0 );
764  if (result==1) {
765  for (int i = firstfile+1; i < lastfile+1; i++) {
766  result = fNuFluxSumChain->Add( Form("%s.%i.root",fileroot.c_str(),i), 0 );
767  }
768  if( (fBr = fNuFluxSumChain->GetBranch("version"))) fNuFluxSumChain->SetBranchAddress("version",&info->version);
769  if( (fBr = fNuFluxSumChain->GetBranch("ntrig")) ) fNuFluxSumChain->SetBranchAddress("ntrig",&info->ntrig);
770  if( (fBr = fNuFluxSumChain->GetBranch("tuneid")) ) fNuFluxSumChain->SetBranchAddress("tuneid",&info->tuneid);
771  if( (fBr = fNuFluxSumChain->GetBranch("pint")) ) fNuFluxSumChain->SetBranchAddress("pint",&info->pint);
772  if( (fBr = fNuFluxSumChain->GetBranch("bpos")) ) fNuFluxSumChain->SetBranchAddress("bpos",info->bpos);
773  if( (fBr = fNuFluxSumChain->GetBranch("btilt")) ) fNuFluxSumChain->SetBranchAddress("btilt",info->btilt);
774  if( (fBr = fNuFluxSumChain->GetBranch("brms")) ) fNuFluxSumChain->SetBranchAddress("brms",info->brms);
775  if( (fBr = fNuFluxSumChain->GetBranch("emit")) ) fNuFluxSumChain->SetBranchAddress("emit",info->emit);
776  if( (fBr = fNuFluxSumChain->GetBranch("alpha")) ) fNuFluxSumChain->SetBranchAddress("alpha",info->alpha);
777  if( (fBr = fNuFluxSumChain->GetBranch("hcur")) ) fNuFluxSumChain->SetBranchAddress("hcur",info->hcur);
778  if( (fBr = fNuFluxSumChain->GetBranch("rand")) ) fNuFluxSumChain->SetBranchAddress("rand",&info->rand);
779  if( (fBr = fNuFluxSumChain->GetBranch("rseed")) ) fNuFluxSumChain->SetBranchAddress("rseed",info->rseed);
780  }
781  }
782 
783  // current ntuple cycle # (flux ntuples may be recycled)
784  fICycle = 1;
785 
786  // sum-up weights & number of neutrinos for the specified location
787  // over a complete cycle. Also record the maximum weight as previous
788  // method using TTree::GetV1() seg faulted for more than ~1.5E6 entries
789  fSumWeightTot1c = 0;
790  fNNeutrinosTot1c = 0;
791  fNDetLocIdFound = 0;
792  for(int ientry = 0; ientry < fNEntries; ientry++) {
793  if (fNuFluxUsingTree)
794  fNuFluxTree->GetEntry(ientry);
795  else
796  fNuFluxChain->GetEntry(ientry);
797  // check for negative flux weights
799  LOG("Flux", pERROR) << "Negative flux weight! Will set weight to 0.0";
800  fPassThroughInfo->norm = 0.0;
801  }
802  fNorm = (double) fPassThroughInfo->norm;
803  // update maximum weight
804  fMaxWeight = TMath::Max(fMaxWeight, (double) fPassThroughInfo->norm);
805  // compare detector location (see GenerateNext_weighted() for details)
806  if(fIsNDLoc && fDetLocId!=fPassThroughInfo->idfd) continue;
809  fNDetLocIdFound++;
810  }
811  // Exit if have not found neutrino at specified location for whole cycle
812  if(fNDetLocIdFound == 0){
813  LOG("Flux", pFATAL)
814  << "The input jnubeam flux ntuple contains no entries for detector id "
815  << fDetLocId << ". Terminating job!";
816  exit(1);
817  }
818  fNDetLocIdFound = 0; // reset the counter
819 
820  LOG("Flux", pNOTICE) << "Maximum flux weight = " << fMaxWeight;
821  if(fMaxWeight <=0 ) {
822  LOG("Flux", pFATAL) << "Non-positive maximum flux weight!";
823  exit(1);
824  }
825 
826  LOG("Flux", pINFO)
827  << "Totals / cycle: #neutrinos = " << fNNeutrinosTot1c
828  << ", Sum{Weights} = " << fSumWeightTot1c;
829 
830  if(fUseRandomOffset){
831  this->RandomOffset(); // Random start point when looping over ntuple
832  }
833 
834  return true;
835 }
836 //___________________________________________________________________________
838 {
839  if(!fPdgCList) {
840  fPdgCList = new PDGCodeList;
841  }
842  fPdgCList->Copy(particles);
843 
844  LOG("Flux", pINFO)
845  << "Declared list of neutrino species: " << *fPdgCList;
846 }
847 //___________________________________________________________________________
849 {
850  fMaxEv = TMath::Max(0.,Ev);
851 
852  LOG("Flux", pINFO)
853  << "Declared maximum flux neutrino energy: " << fMaxEv;
854 }
855 //___________________________________________________________________________
857 {
858 // The flux normalization is in /N POT/det [ND] or /N POT/cm^2 [FD]
859 // This method specifies N (typically 1E+21)
860 
861  fFilePOT = pot;
862 }
863 //___________________________________________________________________________
865 {
866 // The flux neutrino position (x,y) is given at the detector coord system
867 // at z=0. This method sets the preferred starting z position upstream of
868 // the upstream detector face. Each flux neutrino will be backtracked from
869 // z=0 to the input z0.
870 
871  fZ0 = z0;
872 }
873 //___________________________________________________________________________
875 {
876 // The flux ntuples can be recycled for a number of times to boost generated
877 // event statistics without requiring enormous beam simulation statistics.
878 // That option determines how many times the driver is going to cycle through
879 // the input flux ntuple.
880 // With n=0 the flux ntuple will be recycled an infinite amount of times so
881 // that the event generation loop can exit only on a POT or event num check.
882 
883  fNCycles = TMath::Max(0, n);
884 }
885 //___________________________________________________________________________
886 void GJPARCNuFlux::GenerateWeighted(bool gen_weighted)
887 {
888 // Ignore the flux weights when getting next flux neutrino. Always set to
889 // true for unweighted event generation but need option to switch off when
890 // using pre-calculated interaction probabilities for each flux ray to speed
891 // up event generation.
892  fGenerateWeighted = gen_weighted;
893 }
894 //___________________________________________________________________________
896 {
897 // Choose a random number between 0-->fNEntries to set as start point for
898 // looping over flux ntuple. May be necessary when looping over very large
899 // flux files as always starting fromthe same point may introduce biases
900 // (inversely proportional to number of cycles). This method resets the
901 // starting value for fIEntry so must be called before any call to GenerateNext
902 // is made.
903 //
904  double ran_frac = RandomGen::Instance()->RndFlux().Rndm();
905  long int offset = (long int) floor(ran_frac * fNEntries);
906  LOG("Flux", pERROR) << "Setting flux driver to start looping over entries "
907  << "with offset of "<< offset;
908  fIEntry = fOffset = offset;
909 }
910 //___________________________________________________________________________
911 void GJPARCNuFlux::Clear(Option_t * opt)
912 {
913 // If opt = "CycleHistory" then:
914 // Reset all counters and state variables from any previous cycles. This
915 // should be called if, for instance, before event generation this flux
916 // driver had been used for pre-calculating interaction probabilities. Does
917 // not reset initial state variables such as flux particle types, POTs etc...
918 //
919  if(std::strcmp(opt, "CycleHistory") == 0){
920  // Reset so that generate de-weighted events
921  this->GenerateWeighted(false);
922  // Reset cycle counters
923  fICycle = 0;
924  fSumWeight = 0;
925  fNNeutrinos = 0;
926  fIEntry = fOffset;
927  fEntriesThisCycle = 0;
928  }
929 }
930 //___________________________________________________________________________
932 {
933  LOG("Flux", pNOTICE) << "Initializing GJPARCNuFlux driver";
934 
935  fMaxEv = 0;
936  fPdgCList = new PDGCodeList;
938 
939  fNuFluxFile = 0;
940  fNuFluxTree = 0;
941  fNuFluxChain = 0;
942  fNuFluxSumTree = 0;
943  fNuFluxSumChain = 0;
944  fNuFluxUsingTree = true;
945  fDetLoc = "";
946  fDetLocId = 0;
947  fNDetLocIdFound = 0;
948  fIsFDLoc = false;
949  fIsNDLoc = false;
950 
951  fNEntries = 0;
952  fIEntry = 0;
954  fOffset = 0;
955  fNorm = 0.;
956  fMaxWeight = 0;
957  fFilePOT = 0;
958  fZ0 = 0;
959  fNCycles = 0;
960  fICycle = 0;
961  fSumWeight = 0;
962  fNNeutrinos = 0;
963  fSumWeightTot1c = 0;
964  fNNeutrinosTot1c = 0;
965  fGenerateWeighted= false;
966  fUseRandomOffset = true;
967  fLoadedNeutrino = false;
968 
969  this->SetDefaults();
970  this->ResetCurrent();
971 }
972 //___________________________________________________________________________
974 {
975 // - Set default neutrino species list (nue, nuebar, numu, numubar) and
976 // maximum energy (25 GeV).
977 // These defaults can be overwritten by user calls (at the driver init) to
978 // GJPARCNuFlux::SetMaxEnergy(double Ev) and
979 // GJPARCNuFlux::SetFluxParticles(const PDGCodeList & particles)
980 // - Set the default file normalization to 1E+21 POT
981 // - Set the default flux neutrino start z position at -5m (z=0 is the
982 // detector centre).
983 // - Set number of cycles to 1
984 
985  LOG("Flux", pNOTICE) << "Setting default GJPARCNuFlux driver options";
986 
987  PDGCodeList particles;
988  particles.push_back(kPdgNuMu);
989  particles.push_back(kPdgAntiNuMu);
990  particles.push_back(kPdgNuE);
991  particles.push_back(kPdgAntiNuE);
992 
993  this->SetFluxParticles(particles);
994  this->SetMaxEnergy(25./*GeV*/);
995  this->SetFilePOT(1E+21);
996  this->SetUpstreamZ(-5.0);
997  this->SetNumOfCycles(1);
998 }
999 //___________________________________________________________________________
1001 {
1002 // reset running values of neutrino pdg-code, 4-position & 4-momentum
1003 // and the input ntuple leaves
1004 
1005  fLoadedNeutrino = false;
1006 
1007  fgPdgC = 0;
1008  fgP4.SetPxPyPzE (0.,0.,0.,0.);
1009  fgX4.SetXYZT (0.,0.,0.,0.);
1010 
1012 }
1013 //___________________________________________________________________________
1015 {
1016  LOG("Flux", pNOTICE) << "Cleaning up...";
1017 
1018  if (fPdgCList) delete fPdgCList;
1019  if (fPassThroughInfo) delete fPassThroughInfo;
1020 
1021  if (fNuFluxFile) {
1022  fNuFluxFile->Close();
1023  delete fNuFluxFile;
1024  }
1025 }
1026 //___________________________________________________________________________
1028 {
1029 // detector location: name -> int id
1030 // sk -> -1
1031 // nd1 -> +1
1032 // nd2 -> +2
1033 // ...
1034 
1035  if(name == "sk" ) return -1;
1036 
1037  TString temp;
1038  for (int i=1; i<51; i++) {
1039  temp.Form("nd%d",i);
1040  if(name == temp.Data()) return i;
1041  }
1042 
1043  return 0;
1044 }
1045 //___________________________________________________________________________
1047 TObject()
1048 {
1049  this->Reset();
1050 }
1051 //___________________________________________________________________________
1054 TObject()
1055 {
1056  fluxentry = info.fluxentry;
1057  fluxfilename = info.fluxfilename;
1058  Enu = info.Enu;
1059  ppid = info.ppid;
1060  mode = info.mode;
1061  ppi = info.ppi;
1062  ppi0 = info.ppi0;
1063  nvtx0 = info.nvtx0;
1064  cospibm = info.cospibm;
1065  cospi0bm = info.cospi0bm;
1066  idfd = info.idfd;
1067  gamom0 = info.gamom0;
1068  gipart = info.gipart;
1069  xnu = info.xnu;
1070  ynu = info.ynu;
1071  rnu = info.rnu;
1072  for(int i = 0; i<3; i++){
1073  nnu[i] = info.nnu[i];
1074  xpi[i] = info.xpi[i];
1075  xpi0[i] = info.xpi0[i];
1076  gpos0[i] = info.gpos0[i];
1077  npi[i] = info.npi[i];
1078  npi0[i] = info.npi0[i];
1079  gvec0[i] = info.gvec0[i];
1080  }
1081  ng = info.ng;
1082  for(int ip = 0; ip<fNgmax; ip++){
1083  gpid[ip] = info.gpid[ip];
1084  gmec[ip] = info.gmec[ip];
1085  gcosbm[ip] = info.gcosbm[ip];
1086  gvx[ip] = info.gvx[ip];
1087  gvy[ip] = info.gvy[ip];
1088  gvz[ip] = info.gvz[ip];
1089  gpx[ip] = info.gpx[ip];
1090  gpy[ip] = info.gpy[ip];
1091  gpz[ip] = info.gpz[ip];
1092  gmat[ip] = info.gmat[ip];
1093  gdistc[ip] = info.gdistc[ip];
1094  gdistal[ip] = info.gdistal[ip];
1095  gdistti[ip] = info.gdistti[ip];
1096  gdistfe[ip] = info.gdistfe[ip];
1097  }
1098  norm = info.norm;
1099  Enusk = info.Enusk;
1100  normsk = info.normsk;
1101  anorm = info.anorm;
1102  version = info.version;
1103  ntrig = info.ntrig;
1104  tuneid = info.tuneid;
1105  pint = info.pint;
1106  rand = info.rand;
1107  for(int i = 0; i < 2; i++){
1108  bpos[i] = info.bpos[i];
1109  btilt[i] = info.btilt[i];
1110  brms[i] = info.brms[i];
1111  emit[i] = info.emit[i];
1112  alpha[i] = info.alpha[i];
1113  rseed[i] = info.rseed[i];
1114  }
1115  for(int i = 0; i < 3; i++) hcur[i] = info.hcur[i];
1116 }
1117 //___________________________________________________________________________
1119  fluxentry= -1;
1120  fluxfilename = "Not-set";
1121  Enu = 0;
1122  ppid = 0;
1123  mode = 0;
1124  ppi = 0.;
1125  norm = 0;
1126  cospibm = 0.;
1127  nvtx0 = 0;
1128  ppi0 = 0.;
1129  idfd = 0;
1130  rnu = -999999.;
1131  xnu = -999999.;
1132  ynu = -999999.;
1133  cospi0bm = -999999.;
1134  gipart = -1;
1135  gamom0 = 0.;
1136  for(int i=0; i<3; i++){
1137  nnu[i] = 0.;
1138  xpi[i] = -999999.;
1139  npi[i] = 0.;
1140  xpi0[i] = -999999.;
1141  npi0[i] = 0.;
1142  gpos0[i] = -999999.;
1143  gvec0[i] = 0.;
1144  }
1145  ng = -1;
1146  for(int ip = 0; ip<fNgmax; ip++){
1147  gpid[ip] = -999999;
1148  gmec[ip] = -999999;
1149  gcosbm[ip] = -999999.;
1150  gvx[ip] = -999999.;
1151  gvy[ip] = -999999.;
1152  gvz[ip] = -999999.;
1153  gpx[ip] = -999999.;
1154  gpy[ip] = -999999.;
1155  gpz[ip] = -999999.;
1156  gmat[ip] = -999999;
1157  gdistc[ip] = -999999.;
1158  gdistal[ip] = -999999.;
1159  gdistti[ip] = -999999.;
1160  gdistfe[ip] = -999999.;
1161  }
1162  Enusk = -999999.;
1163  normsk = -999999.;
1164  anorm = -999999.;
1165  version = -999999.;
1166  ntrig = -999999;
1167  tuneid = -999999;
1168  pint = -999999;
1169  rand = -999999;
1170  for(int i = 0; i < 2; i++){
1171  bpos[i] = -999999.;
1172  btilt[i] = -999999.;
1173  brms[i] = -999999.;
1174  emit[i] = -999999.;
1175  alpha[i] = -999999.;
1176  rseed[i] = -999999;
1177  }
1178  for(int i = 0; i < 3; i++) hcur[i] = -999999.;
1179 }
1180 //___________________________________________________________________________
1181 namespace genie {
1182 namespace flux {
1183  ostream & operator << (
1185  {
1186  stream << "\n idfd = " << info.idfd
1187  << "\n norm = " << info.norm
1188  << "\n flux entry = " << info.fluxentry
1189  << "\n flux file = " << info.fluxfilename
1190  << "\n Enu = " << info.Enu
1191  << "\n geant code = " << info.ppid
1192  << "\n (pdg code) = " << pdg::GeantToPdg(info.ppid)
1193  << "\n decay mode = " << info.mode
1194  << "\n nvtx0 = " << info.nvtx0
1195  << "\n |momentum| @ decay = " << info.ppi
1196  << "\n position_vector @ decay = ("
1197  << info.xpi[0] << ", "
1198  << info.xpi[1] << ", "
1199  << info.xpi[2] << ")"
1200  << "\n direction_vector @ decay = ("
1201  << info.npi[0] << ", "
1202  << info.npi[1] << ", "
1203  << info.npi[2] << ")"
1204  << "\n |momentum| @ prod. = " << info.ppi0
1205  << "\n position_vector @ prod. = ("
1206  << info.xpi0[0] << ", "
1207  << info.xpi0[1] << ", "
1208  << info.xpi0[2] << ")"
1209  << "\n direction_vector @ prod. = ("
1210  << info.npi0[0] << ", "
1211  << info.npi0[1] << ", "
1212  << info.npi0[2] << ")"
1213  << "\n cospibm = " << info.cospibm
1214  << "\n cospi0bm = " << info.cospi0bm
1215  << "\n Plus additional info if flux version is later than 07a"
1216  << endl;
1217 
1218  return stream;
1219  }
1220 }//flux
1221 }//genie
1222 //___________________________________________________________________________
1223 
1224 
TTree * fNuFluxTree
input flux ntuple
Definition: GJPARCNuFlux.h:116
const XML_Char XML_Encoding * info
Definition: expat.h:530
const XML_Char * name
Definition: expat.h:151
static const double m
Definition: Units.h:79
double POT_curravg(void)
current average POT
A class for generating concrete GFluxI derived classes based on the factory pattern. This code supplies a CPP macro which allows the classes to self-register and thus no modification of this class is needed in order to expand the list of classes it knows about.
int fgPdgC
running generated nu pdg-code
Definition: GJPARCNuFlux.h:111
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:52
const int kPdgNuE
Definition: PDGCodes.h:28
double fSumWeightTot1c
total sum of weights for neutrinos to be thrown / cycle
Definition: GJPARCNuFlux.h:138
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
bool fLoadedNeutrino
set to true when GenerateNext_weighted has been called successfully
Definition: GJPARCNuFlux.h:142
void SetUpstreamZ(double z0)
set flux neutrino initial z position (upstream of the detector)
int DLocName2Id(string name)
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
string fDetLoc
input detector location (&#39;sk&#39;,&#39;nd1&#39;,&#39;nd2&#39;,...)
Definition: GJPARCNuFlux.h:121
TChain * fNuFluxSumChain
input summary ntuple for flux file. This tree is only present for later flux versions > 10a ...
Definition: GJPARCNuFlux.h:119
double POT_1cycle(void)
flux POT per cycle
bool fNuFluxUsingTree
are we using a TTree or a TChain to view the input flux file?
Definition: GJPARCNuFlux.h:120
const int kPdgAntiNuE
Definition: PDGCodes.h:29
#define pFATAL
Definition: Messenger.h:57
const int kPdgNuMu
Definition: PDGCodes.h:30
double fFilePOT
file POT normalization, typically 1E+21
Definition: GJPARCNuFlux.h:132
TChain * fNuFluxChain
input flux ntuple
Definition: GJPARCNuFlux.h:117
double Weight(void)
returns the flux neutrino weight (if any)
Definition: GJPARCNuFlux.h:64
bool ExistsInPDGCodeList(int pdg_code) const
void GenerateWeighted(bool gen_weighted=true)
set whether to generate weighted or unweighted neutrinos
TString ip
Definition: loadincs.C:5
string filename
Definition: shutoffs.py:106
long int fOffset
start looping at entry fOffset
Definition: GJPARCNuFlux.h:129
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GJPARCNuFlux.h:137
void SetFluxParticles(const PDGCodeList &particles)
specify list of flux neutrino species
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
void Initialize(Bool_t useTMVAStyle=kTRUE)
Definition: tmvaglob.C:184
Loaders::FluxType flux
TTree * fNuFluxSumTree
input summary ntuple for flux file. This tree is only present for later flux versions > 10a ...
Definition: GJPARCNuFlux.h:118
A list of PDG codes.
Definition: PDGCodeList.h:33
A GENIE flux driver encapsulating the JPARC neutrino flux. It reads-in the official JPARC neutrino fl...
Definition: GJPARCNuFlux.h:51
void SetMaxEnergy(double Ev)
specify maximum flx neutrino energy
int GeantToPdg(int geant_code)
Definition: PDGUtils.cxx:374
long int fNNeutrinosTot1c
total number of flux neutrinos to be thrown / cycle
Definition: GJPARCNuFlux.h:139
TLorentzVector fgP4
running generated nu 4-momentum
Definition: GJPARCNuFlux.h:112
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Definition: typedefs.hpp:11
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
int fDetLocId
input detector location id (fDetLoc -> jnubeam idfd)
Definition: GJPARCNuFlux.h:122
ClassImp(GJPARCNuFluxPassThroughInfo) GJPARCNuFlux
Float_t E
Definition: plot.C:20
void SetFilePOT(double pot)
flux file norm is in /N POT/det [ND] or /N POT/cm^2 [FD]. Specify N (typically 1E+21) ...
const int fNgmax
Definition: GJPARCNuFlux.h:152
#define pot
int fNCycles
how many times to cycle through the flux ntuple
Definition: GJPARCNuFlux.h:134
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:63
int fNDetLocIdFound
per cycle keep track of the number of fDetLocId are found - if this is zero will exit job ...
Definition: GJPARCNuFlux.h:123
long int fIEntry
current flux ntuple entry
Definition: GJPARCNuFlux.h:127
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
bool fIsNDLoc
input location is a &#39;near&#39; detector location?
Definition: GJPARCNuFlux.h:125
#define pWARN
Definition: Messenger.h:61
GJPARCNuFluxPassThroughInfo * fPassThroughInfo
Definition: GJPARCNuFlux.h:144
OStream cout
Definition: OStream.cxx:6
const XML_Char * version
Definition: expat.h:187
bool fGenerateWeighted
generate weighted/deweighted flux neutrinos (default is false)
Definition: GJPARCNuFlux.h:140
TFile * fNuFluxFile
input flux file
Definition: GJPARCNuFlux.h:115
double fMaxWeight
max flux neutrino weight in input file for the specified detector location
Definition: GJPARCNuFlux.h:131
int fICycle
current cycle
Definition: GJPARCNuFlux.h:135
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void Clear(Option_t *opt)
reset state variables based on opt
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:42
friend ostream & operator<<(ostream &stream, const GJPARCNuFluxPassThroughInfo &info)
long int Index(void)
returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current ent...
bool fIsFDLoc
input location is a &#39;far&#39; detector location?
Definition: GJPARCNuFlux.h:124
double fZ0
configurable starting z position for each flux neutrino (in detector coord system) ...
Definition: GJPARCNuFlux.h:133
PDGCodeList * fPdgCList
list of neutrino pdg-codes
Definition: GJPARCNuFlux.h:109
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
bool GenerateNext(void)
generate the next flux neutrino (return false in err)
exit(0)
long int fNEntries
number of flux ntuple entries
Definition: GJPARCNuFlux.h:126
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)
TRandom3 r(0)
TLorentzVector fgX4
running generated nu 4-position
Definition: GJPARCNuFlux.h:113
static const double cm
Definition: Units.h:76
bool fUseRandomOffset
whether set random starting point when looping over flux ntuples
Definition: GJPARCNuFlux.h:141
#define pNOTICE
Definition: Messenger.h:62
void RandomOffset(void)
choose a random offset as starting entry in flux ntuple
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:72
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:64
double fNorm
current flux ntuple normalisation
Definition: GJPARCNuFlux.h:130
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
#define FLUXDRIVERREG4(_nsa, _nsb, _name, _fqname)
void Copy(const PDGCodeList &list)
copy / print
double fSumWeight
sum of weights for neutrinos thrown so far
Definition: GJPARCNuFlux.h:136
long int fEntriesThisCycle
keep track of number of entries used so far for this cycle
Definition: GJPARCNuFlux.h:128
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
double fMaxEv
maximum energy
Definition: GJPARCNuFlux.h:108
#define pDEBUG
Definition: Messenger.h:64