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