GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GNuMIFlux.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  Robert Hatcher <rhatcher@fnal.gov>
7  Fermi National Accelerator Laboratory
8 
9  Costas Andreopoulos <c.andreopoulos \at cern.ch>
10  University of Liverpool
11 */
12 //____________________________________________________________________________
13 
14 #include <cstdlib>
15 #include <fstream>
16 #include <vector>
17 #include <sstream>
18 #include <cassert>
19 #include <climits>
20 
21 #include "libxml/xmlmemory.h"
22 #include "libxml/parser.h"
23 
26 
27 #include <TFile.h>
28 #include <TChain.h>
29 #include <TChainElement.h>
30 #include <TSystem.h>
31 #include <TStopwatch.h>
32 
34 #include "Framework/Conventions/GBuild.h"
35 
36 #include "Tools/Flux/GNuMIFlux.h"
38 #include "Tools/Flux/GNuMINtuple/g3numi.C"
40 #include "Tools/Flux/GNuMINtuple/g4numi.C"
42 #include "Tools/Flux/GNuMINtuple/flugg.C"
43 
51 
52 using std::endl;
53 
54 #include <vector>
55 #include <algorithm>
56 #include <iomanip>
57 #include "TRegexp.h"
58 #include "TString.h"
59 
62 
63 #ifdef GNUMI_TEST_XY_WGT
64 static genie::flux::xypartials gpartials; // global one used by CalcEnuWgt()
65 #endif
66 
67 using namespace genie;
68 using namespace genie::flux;
69 
70 // declaration of helper class
71 namespace genie {
72  namespace flux {
74  public:
75  GNuMIFluxXMLHelper(GNuMIFlux* gnumi) : fVerbose(0), fGNuMI(gnumi) { ; }
77  bool LoadConfig(std::string cfg);
78 
79  // these should go in a more general package
80  std::vector<double> GetDoubleVector(std::string str);
81  std::vector<long int> GetIntVector(std::string str);
82 
83  private:
84  bool LoadParamSet(xmlDocPtr&, std::string cfg);
85  void ParseParamSet(xmlDocPtr&, xmlNodePtr&);
86  void ParseBeamDir(xmlDocPtr&, xmlNodePtr&);
87  void ParseBeamPos(std::string);
88  void ParseRotSeries(xmlDocPtr&, xmlNodePtr&);
89  void ParseWindowSeries(xmlDocPtr&, xmlNodePtr&);
90  void ParseEnuMax(std::string);
91  TVector3 AnglesToAxis(double theta, double phi, std::string units = "deg");
92  TVector3 ParseTV3(const std::string& );
93 
94  int fVerbose; ///< how noisy to be when parsing XML
95  // who to apply these changes to
97 
98  // derived offsets/rotations
99  TVector3 fBeamPosXML;
100  TRotation fBeamRotXML;
101  TVector3 fFluxWindowPtXML[3];
102  };
103  }
104 }
105 
107 
108 // some nominal positions used in the original g3 ntuple
109 const TLorentzVector kPosCenterNearBeam(0.,0., 1039.35,0.);
110 const TLorentzVector kPosCenterFarBeam (0.,0.,735340.00,0.);
111 
112 //____________________________________________________________________________
114  GFluxExposureI(genie::flux::kPOTs)
115 {
116  this->Initialize();
117 }
118 //___________________________________________________________________________
119 GNuMIFlux::~GNuMIFlux()
120 {
121  this->CleanUp();
122 }
123 
124 //___________________________________________________________________________
125 double GNuMIFlux::GetTotalExposure() const
126 {
127  // complete the GFluxExposureI interface
128  return UsedPOTs();
129 }
130 
131 //___________________________________________________________________________
132 long int GNuMIFlux::NFluxNeutrinos(void) const
133 {
134  /// number of flux neutrinos ray generated so far
135  return fNNeutrinos;
136 }
137 
138 //___________________________________________________________________________
139 bool GNuMIFlux::GenerateNext(void)
140 {
141 // Get next (unweighted) flux ntuple entry on the specified detector location
142 //
144  while ( true ) {
145  // Check for end of flux ntuple
146  bool end = this->End();
147  if ( end ) {
148  LOG("Flux", pNOTICE) << "GenerateNext signaled End() ";
149  return false;
150  }
151 
152  // Get next weighted flux ntuple entry
153  bool nextok = this->GenerateNext_weighted();
154  if ( fGenWeighted ) return nextok;
155  if ( ! nextok ) continue;
156 
157  /* RWH - debug purposes
158  if ( fNCycles == 0 ) {
159  LOG("Flux", pNOTICE)
160  << "Got flux entry: " << fIEntry
161  << " - Cycle: "<< fICycle << "/ infinite";
162  } else {
163  LOG("Flux", pNOTICE)
164  << "Got flux entry: "<< fIEntry
165  << " - Cycle: "<< fICycle << "/"<< fNCycles;
166  }
167  */
168 
169  // Get fractional weight & decide whether to accept curr flux neutrino
170  double f = this->Weight() / fMaxWeight;
171  //LOG("Flux", pNOTICE)
172  // << "Curr flux neutrino fractional weight = " << f;
173  if (f > 1.) {
174  fMaxWeight = this->Weight() * fMaxWgtFudge; // bump the weight
175  LOG("Flux", pERROR)
176  << "** Fractional weight = " << f
177  << " > 1 !! Bump fMaxWeight estimate to " << fMaxWeight
178  << PassThroughInfo();
179  }
180  double r = (f < 1.) ? rnd->RndFlux().Rndm() : 0;
181  bool accept = ( r < f );
182  if ( accept ) {
183 
184 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
185  LOG("Flux", pNOTICE)
186  << "Generated beam neutrino: "
187  << "\n pdg-code: " << fCurEntry->fgPdgC
188  << "\n p4: " << utils::print::P4AsShortString(&(fCurEntry->fgP4))
189  << "\n x4: " << utils::print::X4AsString(&(fCurEntry->fgX4))
190  << "\n p4: " << utils::print::P4AsShortString(&(fCurEntry->fgP4User))
191  << "\n x4: " << utils::print::X4AsString(&(fCurEntry->fgX4User));
192 #endif
193 
194  fWeight = 1.;
195  return true;
196  }
197 
198  //LOG("Flux", pNOTICE)
199  // << "** Rejecting current flux neutrino based on the flux weight only";
200  }
201  return false;
202 }
203 //___________________________________________________________________________
204 bool GNuMIFlux::GenerateNext_weighted(void)
205 {
206 // Get next (weighted) flux ntuple entry on the specified detector location
207 //
208 
209  // Check whether a flux ntuple has been loaded
210  if ( ! fG3NuMI && ! fG4NuMI && ! fFlugg ) {
211  LOG("Flux", pFATAL)
212  << "The flux driver has not been properly configured";
213  //return false; // don't do this - creates an infinite loop!
214  exit(1);
215  }
216 
217  // Reuse an entry?
218  //std::cout << " ***** iuse " << fIUse << " nuse " << fNUse
219  // << " ientry " << fIEntry << " nentry " << fNEntries
220  // << " icycle " << fICycle << " ncycle " << fNCycles << std::endl;
221  if ( fIUse < fNUse && fIEntry >= 0 ) {
222  // Reuse this entry
223  fIUse++;
224  } else {
225  // Reset previously generated neutrino code / 4-p / 4-x
226  this->ResetCurrent();
227  // Move on, read next flux ntuple entry
228  fIEntry++;
229  if ( fIEntry >= fNEntries ) {
230  // Ran out of entries @ the current cycle of this flux file
231  // Check whether more (or infinite) number of cycles is requested
232  if ( fICycle < fNCycles || fNCycles == 0 ) {
233  fICycle++;
234  fIEntry=0;
235  } else {
236  LOG("Flux", pWARN)
237  << "No more entries in input flux neutrino ntuple, cycle "
238  << fICycle << " of " << fNCycles;
239  fEnd = true;
240  // assert(0);
241  return false;
242  }
243  }
244 
245  if ( fG3NuMI ) {
246  fG3NuMI->GetEntry(fIEntry);
247  fCurEntry->MakeCopy(fG3NuMI);
248  } else if ( fG4NuMI ) {
249  fG4NuMI->GetEntry(fIEntry);
250  fCurEntry->MakeCopy(fG4NuMI);
251  } else if ( fFlugg ) {
252  fFlugg->GetEntry(fIEntry);
253  fCurEntry->MakeCopy(fFlugg);
254  } else {
255  LOG("Flux", pERROR) << "No ntuple configured";
256  fEnd = true;
257  //assert(0);
258  return false;
259  }
260 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
261  LOG("Flux",pDEBUG)
262  << "got " << fNNeutrinos << " new fIEntry " << fIEntry
263  << " evtno " << fCurEntry->evtno;
264 #endif
265 
266  fIUse = 1;
267  fCurEntry->pcodes = 0; // fetched entry has geant codes
268  fCurEntry->units = 0; // fetched entry has original units
269 
270  // Convert the current gnumi neutrino flavor mode into a neutrino pdg code
271  // Also convert other particle codes in GNuMIFluxPassThroughInfo to PDG
272  fCurEntry->ConvertPartCodes();
273  // here we might want to do flavor oscillations or simple mappings
274  fCurEntry->fgPdgC = fCurEntry->ntype;
275  }
276 
277  // Check neutrino pdg against declared list of neutrino species declared
278  // by the current instance of the NuMI neutrino flux driver.
279  // No undeclared neutrino species will be accepted at this point as GENIE
280  // has already been configured to handle the specified list.
281  // Make sure that the appropriate list of flux neutrino species was set at
282  // initialization via GNuMIFlux::SetFluxParticles(const PDGCodeList &)
283 
284  // update the # POTs, number of neutrinos
285  // do this HERE (before rejecting flavors that users might be weeding out)
286  // in order to keep the POT accounting correct. This allows one to get
287  // the right normalization for generating only events from the intrinsic
288  // nu_e entries.
289  fAccumPOTs += fEffPOTsPerNu / fMaxWeight;
290  fNNeutrinos++;
291 
292  if ( ! fPdgCList->ExistsInPDGCodeList(fCurEntry->fgPdgC) ) {
293  /// user might modify list via SetFluxParticles() in order to reject certain
294  /// flavors, even if they're found in the file. So don't make a big fuss.
295  /// Spit out a single message and then stop reporting that flavor as problematic.
296  int badpdg = fCurEntry->fgPdgC;
297  if ( ! fPdgCListRej->ExistsInPDGCodeList(badpdg) ) {
298  fPdgCListRej->push_back(badpdg);
299  LOG("Flux", pWARN)
300  << "Encountered neutrino specie (" << badpdg
301  << " pcodes=" << fCurEntry->pcodes << ")"
302  << " that wasn't in SetFluxParticles() list, "
303  << "\nDeclared list of neutrino species: " << *fPdgCList;
304  }
305  return false;
306  }
307 
308  // Update the curr neutrino weight and energy
309 
310  // Check current neutrino energy against the maximum flux neutrino energy
311  // declared by the current instance of the NuMI neutrino flux driver.
312  // No flux neutrino exceeding that maximum energy will be accepted at this
313  // point as that maximum energy has already been used for normalizing the
314  // interaction probabilities.
315  // Make sure that the appropriate maximum flux neutrino energy was set at
316  // initialization via GNuMIFlux::SetMaxEnergy(double Ev)
317 
318  fCurEntry->fgX4 = fFluxWindowBase;
319 
320  double Ev = 0;
321  double& wgt_xy = fCurEntry->fgXYWgt;
322  switch ( fUseFluxAtDetCenter ) {
323  case -1: // near detector
324  wgt_xy = fCurEntry->nwtnear;
325  Ev = fCurEntry->nenergyn;
326  break;
327  case +1: // far detector
328  wgt_xy = fCurEntry->nwtfar;
329  Ev = fCurEntry->nenergyf;
330  break;
331  default: // recalculate on x-y window
332  RandomGen * rnd = RandomGen::Instance();
333  fCurEntry->fgX4 += ( rnd->RndFlux().Rndm()*fFluxWindowDir1 +
334  rnd->RndFlux().Rndm()*fFluxWindowDir2 );
335  fCurEntry->CalcEnuWgt(fCurEntry->fgX4,Ev,wgt_xy);
336  break;
337  }
338 
339  if (Ev > fMaxEv) {
340  LOG("Flux", pWARN)
341  << "Flux neutrino energy exceeds declared maximum neutrino energy"
342  << "\nEv = " << Ev << "(> Ev{max} = " << fMaxEv << ")";
343  }
344 
345  // Set the current flux neutrino 4-momentum
346  // this is in *beam* coordinates
347  fgX4dkvtx = TLorentzVector( fCurEntry->vx,
348  fCurEntry->vy,
349  fCurEntry->vz, 0.);
350  // don't use TLorentzVector here for Mag() due to - on metric
351  // normalize direction to 1.0
352  TVector3 dirNu = (fCurEntry->fgX4.Vect() - fgX4dkvtx.Vect()).Unit();
353  fCurEntry->fgP4.SetPxPyPzE( Ev*dirNu.X(),
354  Ev*dirNu.Y(),
355  Ev*dirNu.Z(), Ev);
356 
357  // calculate the weight, potentially includes effect from tilted window
358  // must be done *after* neutrino direction is determined
359  fWeight = fCurEntry->nimpwt * fCurEntry->fgXYWgt; // full weight
360  if ( fApplyTiltWeight ) {
361  double tiltwgt = dirNu.Dot( FluxWindowNormal() );
362  fWeight *= TMath::Abs( tiltwgt );
363  }
364 
365  // update sume of weights
366  fSumWeight += this->Weight();
367 
368  // Set the current flux neutrino 4-position, direction in user coord
369  Beam2UserP4(fCurEntry->fgP4,fCurEntry->fgP4User);
370  Beam2UserPos(fCurEntry->fgX4,fCurEntry->fgX4User);
371 
372  // if desired, move to user specified user coord z
373  if ( TMath::Abs(fZ0) < 1.0e30 ) this->MoveToZ0(fZ0);
374 
375 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
376  LOG("Flux", pINFO)
377  << "Generated neutrino: " << fIEntry << " " << fCurEntry->evtno
378  << " nenergyn " << fCurEntry->nenergyn
379  << "\n pdg-code: " << fCurEntry->fgPdgC
380  << "\n p4 beam: " << utils::print::P4AsShortString(&fCurEntry->fgP4)
381  << "\n x4 beam: " << utils::print::X4AsString(&fCurEntry->fgX4)
382  << "\n p4 user: " << utils::print::P4AsShortString(&(fCurEntry->fgP4User))
383  << "\n x4 user: " << utils::print::X4AsString(&(fCurEntry->fgX4User));
384 #endif
385  if ( Ev > fMaxEv ) {
386  LOG("Flux", pFATAL)
387  << "Generated neutrino had E_nu = " << Ev << " > " << fMaxEv
388  << " maximum ";
389  assert(0);
390  }
391 
392  return true;
393 }
394 //___________________________________________________________________________
395 double GNuMIFlux::GetDecayDist() const
396 {
397  // return distance (user units) between dk point and start position
398  // these are in beam units
399  TVector3 x3diff = fCurEntry->fgX4.Vect() - fgX4dkvtx.Vect();
400  return x3diff.Mag() * fLengthScaleB2U;
401 }
402 //___________________________________________________________________________
403 void GNuMIFlux::MoveToZ0(double z0usr)
404 {
405  // move ray origin to specified user z0
406  // move beam coord entry correspondingly
407 
408 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
409  LOG("Flux", pNOTICE)
410  << "MoveToZ0 (z0usr=" << z0usr << ") before:"
411  << "\n p4 user: " << utils::print::P4AsShortString(&(fCurEntry->fgP4User))
412  << "\n x4 user: " << utils::print::X4AsString(&(fCurEntry->fgX4User));
413 #endif
414 
415  double pzusr = fCurEntry->fgP4User.Pz();
416  if ( TMath::Abs(pzusr) < 1.0e-30 ) {
417  // neutrino is moving almost entirely in x-y plane
418  LOG("Flux", pWARN)
419  << "MoveToZ0(" << z0usr << ") not possible due to pz_usr (" << pzusr << ")";
420  return;
421  }
422 
423  double scale = (z0usr - fCurEntry->fgX4User.Z()) / pzusr;
424  fCurEntry->fgX4User += (scale*fCurEntry->fgP4User);
425  fCurEntry->fgX4 += ((fLengthScaleU2B*scale)*fCurEntry->fgP4);
426  // this scaling works for distances, but not the time component
427  fCurEntry->fgX4.SetT(0);
428  fCurEntry->fgX4User.SetT(0);
429 
430 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
431  LOG("Flux", pNOTICE)
432  << "MoveToZ0 (" << z0usr << ") after:"
433  << "\n x4 user: " << utils::print::X4AsString(&(fCurEntry->fgX4User));
434 #endif
435 
436 }
437 
438 //___________________________________________________________________________
439 void GNuMIFlux::CalcEffPOTsPerNu()
440 {
441  // do this if flux window changes or # of files changes
442 
443  if (!fNuFluxTree) return; // not yet fully configured
444 
445  // effpots = mc_pots * (wgtfunction-area) / window-area / wgt-max-est
446  // wgtfunction-area = pi * radius-det-element^2 = pi * (100.cm)^2
447 
448  // this should match what is used in the CalcEnuWgt()
449  const double kRDET = 100.0; // set to flux per 100 cm radius
450  const double kRDET2 = kRDET * kRDET;
451  double flux_area = fFluxWindowDir1.Vect().Cross(fFluxWindowDir2.Vect()).Mag();
452  LOG("Flux",pNOTICE) << "in CalcEffPOTsPerNu, area = " << flux_area;
453 
454  if ( flux_area < 1.0e-30 ) {
455  LOG("Flux", pWARN)
456  << "CalcEffPOTsPerNu called with flux window area effectively zero";
457  flux_area = 1;
458  }
459  double area_ratio = TMath::Pi() * kRDET2 / flux_area;
460  fEffPOTsPerNu = area_ratio * ( (double)fFilePOTs / (double)fNEntries );
461 }
462 
463 //___________________________________________________________________________
464 double GNuMIFlux::UsedPOTs(void) const
465 {
466 // Compute current number of flux POTs
467 
468  if (!fNuFluxTree) {
469  LOG("Flux", pWARN)
470  << "The flux driver has not been properly configured";
471  return 0;
472  }
473  return fAccumPOTs;
474 }
475 
476 //___________________________________________________________________________
477 double GNuMIFlux::POT_curr(void) {
478  // RWH: Not sure what POT_curr is supposed to represent I'll guess for
479  // now that that it means what I mean by UsedPOTs().
480  return UsedPOTs();
481 }
482 
483 //___________________________________________________________________________
484 void GNuMIFlux::LoadBeamSimData(const std::vector<std::string>& patterns,
485  const std::string& config )
486 {
487 // Loads in a gnumi beam simulation root file (converted from hbook format)
488 // into the GNuMIFlux driver.
489 
490  bool found_cfg = this->LoadConfig(config);
491  if ( ! found_cfg ) {
492  LOG("Flux", pFATAL)
493  << "LoadBeamSimData could not find XML config \"" << config << "\"\n";
494  exit(1);
495  }
496 
497  fNuFluxFilePatterns = patterns;
498  std::vector<int> nfiles_from_pattern;
499 
500  // create a (sorted) set of file names
501  // this also helps ensure that the same file isn't listed multiple times
502  std::set<std::string> fnames;
503 
504  LOG("Flux", pINFO) << "LoadBeamSimData was passed " << patterns.size()
505  << " patterns";
506 
507  for (size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
508  string pattern = patterns[ipatt];
509  nfiles_from_pattern.push_back(0);
510 
511  LOG("Flux", pNOTICE)
512  << "Loading gnumi flux tree from ROOT file pattern ["
513  << std::setw(3) << ipatt << "] \"" << pattern << "\"";
514 
515  // !WILDCARD only works for file name ... NOT directory
516  string dirname = gSystem->UnixPathName(gSystem->WorkingDirectory());
517  size_t slashpos = pattern.find_last_of("/");
518  size_t fbegin;
519  if ( slashpos != std::string::npos ) {
520  dirname = pattern.substr(0,slashpos);
521  LOG("Flux", pINFO) << "Look for flux using directory " << dirname;
522  fbegin = slashpos + 1;
523  } else { fbegin = 0; }
524 
525  void* dirp = gSystem->OpenDirectory(gSystem->ExpandPathName(dirname.c_str()));
526  if ( dirp ) {
527  std::string basename =
528  pattern.substr(fbegin,pattern.size()-fbegin);
529  TRegexp re(basename.c_str(),kTRUE);
530  const char* onefile;
531  while ( ( onefile = gSystem->GetDirEntry(dirp) ) ) {
532  TString afile = onefile;
533  if ( afile=="." || afile==".." ) continue;
534  if ( basename!=afile && afile.Index(re) == kNPOS ) continue;
535  std::string fullname = dirname + "/" + afile.Data();
536  fnames.insert(fullname);
537  nfiles_from_pattern[ipatt]++;
538  }
539  gSystem->FreeDirectory(dirp);
540  } // legal directory
541  } // loop over patterns
542 
543  size_t indx = 0;
544  std::set<string>::const_iterator sitr = fnames.begin();
545  for ( ; sitr != fnames.end(); ++sitr, ++indx ) {
546  string filename = *sitr;
547  //std::cout << " [" << std::setw(3) << indx << "] \""
548  // << filename << "\"" << std::endl;
549  bool isok = true;
550  // this next test only works for local files, so we can't do that
551  // if we want to stream via xrootd
552  // ! (gSystem->AccessPathName(filename.c_str()));
553  if ( isok ) {
554  // open the file to see what it contains
555  // h10 => g3numi _or_ flugg
556  // nudata => g4numi
557  // for now distinguish between g3numi/flugg using file name
558  TFile* tf = TFile::Open(filename.c_str(),"READ");
559  int isflugg = ( filename.find("flugg") != string::npos ) ? 1 : 0;
560  const std::string tnames[] = { "h10", "nudata" };
561  const std::string gnames[] = { "g3numi","g4numi","flugg","g4flugg"};
562  for (int j = 0; j < 2 ; ++j ) {
563  TTree* atree = (TTree*)tf->Get(tnames[j].c_str());
564  if ( atree ) {
565  const std::string tname_this = tnames[j];
566  const std::string gname_this = gnames[j+2*isflugg];
567  // create chain if none exists
568  if ( ! fNuFluxTree ) {
569  this->SetTreeName(tname_this);
570  fNuFluxTree = new TChain(fNuFluxTreeName.c_str());
571  fNuFluxGen = gname_this;
572  // here we should scan for estimated POTs/file
573  // also maybe the check on max weight
574  }
575  // sanity check for mixing g3/g4/flugg files
576  if ( fNuFluxTreeName != tname_this ||
577  fNuFluxGen != gname_this ) {
578  LOG("Flux", pFATAL)
579  << "Inconsistent flux file types\n"
580  << "The input gnumi flux file \"" << filename
581  << "\"\ncontains a '" << tname_this << "' " << gname_this
582  << "numi ntuple "
583  << "but a '" << fNuFluxTreeName << "' " << fNuFluxGen
584  << " numi ntuple has alread been seen in the chain";
585  exit(1);
586  } // sanity mix/match g3/g4
587  // add the file to the chain
588  this->AddFile(atree,filename);
589  } // found a tree
590  } // loop over either g3 or g4
591  tf->Close();
592  delete tf;
593  } // loop over tree type
594  } // loop over sorted file names
595 
596  if ( fNuFluxTreeName == "" ) {
597  LOG("Flux", pFATAL)
598  << "The input gnumi flux file doesn't exist! Initialization failed!";
599  exit(1);
600  }
601  if ( fNuFluxGen == "g3numi" ) fG3NuMI = new g3numi(fNuFluxTree);
602  if ( fNuFluxGen == "g4numi" ) fG4NuMI = new g4numi(fNuFluxTree);
603  if ( fNuFluxGen == "flugg" ) fFlugg = new flugg(fNuFluxTree);
604 
605  // this will open all files and read header!!
606  fNEntries = fNuFluxTree->GetEntries();
607 
608  if ( fNEntries == 0 ) {
609  LOG("Flux", pERROR)
610  << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
611  LOG("Flux", pERROR)
612  << "Loaded flux tree contains " << fNEntries << " entries";
613  LOG("Flux", pERROR)
614  << "Was passed the file patterns: ";
615  for (size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
616  string pattern = patterns[ipatt];
617  LOG("Flux", pERROR)
618  << " [" << std::setw(3) << ipatt <<"] " << pattern;
619  }
620  LOG("Flux", pERROR)
621  << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
622  } else {
623  LOG("Flux", pNOTICE)
624  << "Loaded flux tree contains " << fNEntries << " entries"
625  << " from " << fnames.size() << " unique files";
626  for (size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
627  string pattern = patterns[ipatt];
628  LOG("Flux", pINFO)
629  << " pattern: " << pattern << " yielded "
630  << nfiles_from_pattern[ipatt] << " files";
631  }
632  }
633 
634  // we have a file we can work with
635  if (!fDetLocIsSet) {
636  LOG("Flux", pERROR)
637  << "LoadBeamSimData left detector location unset";
638  }
639  if (fMaxWeight<=0) {
640  LOG("Flux", pINFO)
641  << "Run ScanForMaxWeight() as part of LoadBeamSimData";
642  this->ScanForMaxWeight();
643  }
644 
645  // current ntuple cycle # (flux ntuples may be recycled)
646  fICycle = 0;
647  // pick a starting entry index [0:fNEntries-1]
648  // pretend we just used up the the previous one
650  fIUse = 9999999;
651  fIEntry = rnd->RndFlux().Integer(fNEntries) - 1;
652 
653  // don't count things we used to estimate max weight
654  fSumWeight = 0;
655  fNNeutrinos = 0;
656  fAccumPOTs = 0;
657 
658  LOG("Flux",pNOTICE) << "about to CalcEffPOTsPerNu";
659  this->CalcEffPOTsPerNu();
660 
661 }
662 //___________________________________________________________________________
663 void GNuMIFlux::GetBranchInfo(std::vector<std::string>& branchNames,
664  std::vector<std::string>& branchClassNames,
665  std::vector<void**>& branchObjPointers)
666 {
667  // allow flux driver to report back current status and/or ntuple entry
668  // info for possible recording in the output file by supplying
669  // the class name, and a pointer to the object that will be filled
670  // as well as a suggested name for the branch.
671 
672  branchNames.push_back("flux");
673  branchClassNames.push_back("genie::flux::GNuMIFluxPassThroughInfo");
674  branchObjPointers.push_back((void**)&fCurEntry);
675 
676 }
677 TTree* GNuMIFlux::GetMetaDataTree() { return 0; } // there is none
678 
679 //___________________________________________________________________________
680 void GNuMIFlux::ScanForMaxWeight(void)
681 {
682  if (!fDetLocIsSet) {
683  LOG("Flux", pERROR)
684  << "Specify a detector location before scanning for max weight";
685  return;
686  }
687 
688  // scan for the maximum weight
689  int ipos_estimator = fUseFluxAtDetCenter;
690  if ( ipos_estimator == 0 ) {
691  // within 100m of a known point?
692  double zbase = fFluxWindowBase.Z();
693  if ( TMath::Abs(zbase-103648.837) < 10000. ) ipos_estimator = -1; // use NearDet
694  if ( TMath::Abs(zbase-73534000. ) < 10000. ) ipos_estimator = +1; // use FarDet
695  }
696  if ( ipos_estimator != 0 ) {
697 
698  //// one can't really be sure which Nwtfar/Nwtnear this refers to
699  //// some gnumi files have "NOvA" weights
700  const char* ntwgtstrv[4] = { "Nimpwt*Nwtnear",
701  "Nimpwt*Nwtfar",
702  "Nimpwt*NWtNear[0]",
703  "Nimpwt*NWtFar[0]" };
704  int strindx = 0;
705  if ( ipos_estimator > 0 ) strindx = 1;
706  if ( fG4NuMI ) strindx += 2;
707  // set upper limit on how many entries to scan
708  Long64_t nscan = TMath::Min(fNEntries,200000LL);
709 
710  fNuFluxTree->Draw(ntwgtstrv[strindx],"","goff",nscan);
711  //std::cout << " Draw \"" << ntwgtstrv[strindx] << "\"" << std::endl;
712  //std::cout << " SelectedRows " << fNuFluxTree->GetSelectedRows()
713  // << " V1 " << fNuFluxTree->GetV1() << std::endl;
714 
715  Long64_t idx = TMath::LocMax(fNuFluxTree->GetSelectedRows(),
716  fNuFluxTree->GetV1());
717  //std::cout << "idx " << idx << " of " << fNuFluxTree->GetSelectedRows() << std::endl;
718  fMaxWeight = fNuFluxTree->GetV1()[idx];
719  LOG("Flux", pNOTICE) << "Maximum flux weight from Nwt in ntuple = "
720  << fMaxWeight;
721  if ( fMaxWeight <= 0 ) {
722  LOG("Flux", pFATAL) << "Non-positive maximum flux weight!";
723  exit(1);
724  }
725  }
726  // the above works only for things close to the MINOS stored weight
727  // values. otherwise we need to work out our own estimate.
728  double wgtgenmx = 0, enumx = 0;
729  TStopwatch t;
730  t.Start();
731  for (int itry=0; itry < fMaxWgtEntries; ++itry) {
732  this->GenerateNext_weighted();
733  double wgt = this->Weight();
734  if ( wgt > wgtgenmx ) wgtgenmx = wgt;
735  double enu = fCurEntry->fgP4.Energy();
736  if ( enu > enumx ) enumx = enu;
737  }
738  t.Stop();
739  t.Print("u");
740  LOG("Flux", pNOTICE) << "Maximum flux weight for spin = "
741  << wgtgenmx << ", energy = " << enumx
742  << " (" << fMaxWgtEntries << ")";
743 
744  if (wgtgenmx > fMaxWeight ) fMaxWeight = wgtgenmx;
745  // apply a fudge factor to estimated weight
746  fMaxWeight *= fMaxWgtFudge;
747  // adjust max energy?
748  if ( enumx*fMaxEFudge > fMaxEv ) {
749  LOG("Flux", pNOTICE) << "Adjust max: was=" << fMaxEv
750  << " now " << enumx << "*" << fMaxEFudge
751  << " = " << enumx*fMaxEFudge;
752  fMaxEv = enumx * fMaxEFudge;
753  }
754 
755  LOG("Flux", pNOTICE) << "Maximum flux weight = " << fMaxWeight
756  << ", energy = " << fMaxEv;
757 
758 }
759 //___________________________________________________________________________
760 void GNuMIFlux::SetMaxEnergy(double Ev)
761 {
762  fMaxEv = TMath::Max(0.,Ev);
763 
764  LOG("Flux", pINFO)
765  << "Declared maximum flux neutrino energy: " << fMaxEv;
766 }
767 //___________________________________________________________________________
768 void GNuMIFlux::SetEntryReuse(long int nuse)
769 {
770 // With nuse > 1 then the same entry in the file is used "nuse" times
771 // before moving on to the next entry in the ntuple
772 
773  fNUse = TMath::Max(1L, nuse);
774 }
775 //___________________________________________________________________________
776 void GNuMIFlux::SetTreeName(string name)
777 {
778  fNuFluxTreeName = name;
779 }
780 //___________________________________________________________________________
781 void GNuMIFlux::UseFluxAtNearDetCenter(void)
782 {
783  SetLengthUnits(genie::utils::units::UnitFromString("m"));
784  fFluxWindowBase = kPosCenterNearBeam;
785  fFluxWindowDir1 = TLorentzVector(); // no extent
786  fFluxWindowDir2 = TLorentzVector();
787  fUseFluxAtDetCenter = -1;
788  fDetLocIsSet = true;
789 }
790 //___________________________________________________________________________
791 void GNuMIFlux::UseFluxAtFarDetCenter(void)
792 {
793  SetLengthUnits(genie::utils::units::UnitFromString("m"));
794  fFluxWindowBase = kPosCenterFarBeam;
795  fFluxWindowDir1 = TLorentzVector(); // no extent
796  fFluxWindowDir2 = TLorentzVector();
797  fUseFluxAtDetCenter = +1;
798  fDetLocIsSet = true;
799 }
800 //___________________________________________________________________________
801 bool GNuMIFlux::SetFluxWindow(GNuMIFlux::StdFluxWindow_t stdwindow, double padding)
802 {
803  // set some standard flux windows
804  // rwh: should also set detector coord transform
805  // rwh: padding allows add constant padding to pre-existing set
806  double padbeam = padding / fLengthScaleB2U; // user might set different units
807  LOG("Flux",pNOTICE)
808  << "SetBeamFluxWindow " << (int)stdwindow << " padding " << padbeam << " cm";
809 
810 
811  switch ( stdwindow ) {
812 #ifdef THIS_IS_NOT_YET_IMPLEMENTED
813  case kMinosNearDet:
814  SetBeamFluxWindow(103648.837);
815  break;
816  case kMinosFarDear:
817  SetBeamFluxWindow(73534000.);
818  break;
819  case kMinosNearRock:
820  SetBeamFluxWindow();
821  break;
822  case kMinosFarRock:
823  SetBeamFluxWindow();
824  break;
825 #endif
826  case kMinosNearCenter:
827  {
828  SetLengthUnits(genie::utils::units::UnitFromString("m"));
829  fFluxWindowBase = kPosCenterNearBeam;
830  fFluxWindowDir1 = TLorentzVector(); // no extent
831  fFluxWindowDir2 = TLorentzVector();
832  TLorentzVector usrpos;
833  Beam2UserPos(fFluxWindowBase, usrpos);
834  fFluxWindowPtUser[0] = usrpos.Vect();
835  fFluxWindowPtUser[1] = fFluxWindowPtUser[0];
836  fFluxWindowPtUser[2] = fFluxWindowPtUser[0];
837  fFluxWindowLen1 = 0;
838  fFluxWindowLen2 = 0;
839  break;
840  }
841  case kMinosFarCenter:
842  {
843  SetLengthUnits(genie::utils::units::UnitFromString("m"));
844  fFluxWindowBase = kPosCenterFarBeam;
845  fFluxWindowDir1 = TLorentzVector(); // no extent
846  fFluxWindowDir2 = TLorentzVector();
847  TLorentzVector usrpos;
848  Beam2UserPos(fFluxWindowBase, usrpos);
849  fFluxWindowPtUser[0] = usrpos.Vect();
850  fFluxWindowPtUser[1] = fFluxWindowPtUser[0];
851  fFluxWindowPtUser[2] = fFluxWindowPtUser[0];
852  fFluxWindowLen1 = 0;
853  fFluxWindowLen2 = 0;
854  break;
855  }
856  default:
857  LOG("Flux", pERROR)
858  << "SetBeamFluxWindow - StdFluxWindow " << stdwindow
859  << " not yet implemented";
860  return false;
861  }
862  LOG("Flux",pNOTICE) << "about to CalcEffPOTsPerNu";
863  this->CalcEffPOTsPerNu();
864  return true;
865 }
866 //___________________________________________________________________________
867 void GNuMIFlux::SetFluxWindow(TVector3 p0, TVector3 p1, TVector3 p2)
868  // bool inDetCoord) future extension
869 {
870  // set flux window
871  // NOTE: internally these are in "cm", but user might have set a preference
872  fUseFluxAtDetCenter = 0;
873  fDetLocIsSet = true;
874 
875  fFluxWindowPtUser[0] = p0;
876  fFluxWindowPtUser[1] = p1;
877  fFluxWindowPtUser[2] = p2;
878 
879  // convert from user to beam coord and from 3 points to base + 2 directions
880  // apply units conversion
881  TLorentzVector ptbm0, ptbm1, ptbm2;
882  User2BeamPos(TLorentzVector(fFluxWindowPtUser[0],0),ptbm0);
883  User2BeamPos(TLorentzVector(fFluxWindowPtUser[1],0),ptbm1);
884  User2BeamPos(TLorentzVector(fFluxWindowPtUser[2],0),ptbm2);
885 
886  fFluxWindowBase = ptbm0;
887  fFluxWindowDir1 = ptbm1 - ptbm0;
888  fFluxWindowDir2 = ptbm2 - ptbm0;
889 
890  fFluxWindowLen1 = fFluxWindowDir1.Mag();
891  fFluxWindowLen2 = fFluxWindowDir2.Mag();
892  fWindowNormal = fFluxWindowDir1.Vect().Cross(fFluxWindowDir2.Vect()).Unit();
893 
894  double dot = fFluxWindowDir1.Dot(fFluxWindowDir2);
895  if ( TMath::Abs(dot) > 1.0e-8 )
896  LOG("Flux",pWARN) << "Dot product between window direction vectors was "
897  << dot << "; please check for orthoganality";
898 
899  LOG("Flux",pNOTICE) << "about to CalcEffPOTsPerNu";
900  this->CalcEffPOTsPerNu();
901 }
902 
903 //___________________________________________________________________________
904 void GNuMIFlux::GetFluxWindow(TVector3& p0, TVector3& p1, TVector3& p2) const
905 {
906  // return flux window points
907  p0 = fFluxWindowPtUser[0];
908  p1 = fFluxWindowPtUser[1];
909  p2 = fFluxWindowPtUser[2];
910 
911 }
912 //___________________________________________________________________________
913 void GNuMIFlux::SetBeamRotation(TRotation beamrot)
914 {
915  // rotation is really only 3-d vector, but we'll be operating on LorentzV's
916  fBeamRot = TLorentzRotation(beamrot);
917  fBeamRotInv = fBeamRot.Inverse();
918 }
919 
920 void GNuMIFlux::SetBeamCenter(TVector3 beam0)
921 {
922  // set coord transform between detector and beam
923  // NOTE: internally these are in "cm", but user might have set a preference
924  fBeamZero = TLorentzVector(beam0,0); // no time shift
925 }
926 
927 //___________________________________________________________________________
928 TRotation GNuMIFlux::GetBeamRotation() const
929 {
930  // rotation is really only 3-d vector, but we'll be operating on LorentzV's
931  // give people back the original TRotation ... not pretty
932  // ... it think this is right
933  TRotation rot3;
934  const TLorentzRotation& rot4 = fBeamRot;
935  TVector3 newX(rot4.XX(),rot4.XY(),rot4.XZ());
936  TVector3 newY(rot4.YX(),rot4.YY(),rot4.YZ());
937  TVector3 newZ(rot4.ZX(),rot4.ZY(),rot4.ZZ());
938  rot3.RotateAxes(newX,newY,newZ);
939  return rot3.Inverse();
940 }
941 TVector3 GNuMIFlux::GetBeamCenter() const
942 {
943  TVector3 beam0 = fBeamZero.Vect();
944  return beam0;
945 }
946 
947 //___________________________________________________________________________
948 //void GNuMIFlux::SetCoordTransform(TVector3 beam0, TRotation beamrot)
949 //{
950 // // set coord transform between detector and beam
951 // // NOTE: internally these are in "cm", but user might have set a preference
952 //
953 // beam0 *= (1./fLengthScaleB2U);
954 // fDetectorZero = TLorentzVector(beam0,0); // no time shift
955 // fDetectorRot = TLorentzRotation(beamrot);
956 //
957 //}
958 //___________________________________________________________________________
959 //void GNuMIFlux::GetDetectorCoord(TLorentzVector& det0, TLorentzRotation& detrot) const
960 //{
961 // // get coord transform between detector and beam
962 // // NOTE: internally these are in "cm", but user might have set a preference
963 //
964 // det0 = fDetectorZero;
965 // det0 *= fLengthScaleB2U;
966 // detrot = fDetectorRot;
967 //
968 //}
969 //___________________________________________________________________________
970 
971 void GNuMIFlux::Beam2UserPos(const TLorentzVector& beamxyz,
972  TLorentzVector& usrxyz) const
973 {
974  usrxyz = fLengthScaleB2U*(fBeamRot*beamxyz) + fBeamZero;
975 }
976 void GNuMIFlux::Beam2UserDir(const TLorentzVector& beamdir,
977  TLorentzVector& usrdir) const
978 {
979  usrdir = fLengthScaleB2U*(fBeamRot*beamdir);
980 }
981 void GNuMIFlux::Beam2UserP4 (const TLorentzVector& beamp4,
982  TLorentzVector& usrp4 ) const
983 {
984  usrp4 = fBeamRot*beamp4;
985 }
986 
987 void GNuMIFlux::User2BeamPos(const TLorentzVector& usrxyz,
988  TLorentzVector& beamxyz) const
989 {
990  beamxyz = fLengthScaleU2B*(fBeamRotInv*(usrxyz-fBeamZero));
991 }
992 void GNuMIFlux::User2BeamDir(const TLorentzVector& usrdir,
993  TLorentzVector& beamdir) const
994 {
995  beamdir = fLengthScaleU2B*(fBeamRotInv*usrdir);
996 }
997 void GNuMIFlux::User2BeamP4 (const TLorentzVector& usrp4,
998  TLorentzVector& beamp4) const
999 {
1000  beamp4 = fBeamRotInv*usrp4;
1001 }
1002 
1003 //___________________________________________________________________________
1004 void GNuMIFlux::PrintCurrent(void)
1005 {
1006  LOG("Flux", pNOTICE) << "CurrentEntry:" << *fCurEntry;
1007 }
1008 //___________________________________________________________________________
1009 void GNuMIFlux::Clear(Option_t * opt)
1010 {
1011  // Clear the driver state
1012  //
1013  LOG("Flux", pWARN) << "GSimpleNtpFlux::Clear(" << opt << ") called";
1014  // do it in all cases, but EVGDriver/GMCJDriver will pass "CycleHistory"
1015 
1016  fICycle = 0;
1017 
1018  fSumWeight = 0;
1019  fNNeutrinos = 0;
1020  fAccumPOTs = 0;
1021 }
1022 //___________________________________________________________________________
1023 void GNuMIFlux::GenerateWeighted(bool gen_weighted)
1024 {
1025  // Set whether to generate weighted rays
1026  //
1027  fGenWeighted = gen_weighted;
1028 }
1029 //___________________________________________________________________________
1031 {
1032  LOG("Flux", pNOTICE) << "Initializing GNuMIFlux driver";
1033 
1034  fMaxEv = 0;
1035  fEnd = false;
1036 
1037  fCurEntry = new GNuMIFluxPassThroughInfo;
1038 
1039  fNuFluxTree = 0;
1040  fG3NuMI = 0;
1041  fG4NuMI = 0;
1042  fFlugg = 0;
1043  fNuFluxTreeName = "";
1044  fNuFluxGen = "";
1045  fNFiles = 0;
1046 
1047  fNEntries = 0;
1048  fIEntry = -1;
1049  fICycle = 0;
1050  fNUse = 1;
1051  fIUse = 999999;
1052 
1053  fNuTot = 0;
1054  fFilePOTs = 0;
1055 
1056  fMaxWeight = -1;
1057  fMaxWgtFudge = 1.05;
1058  fMaxWgtEntries = 2500000;
1059  fMaxEFudge = 0;
1060 
1061  fSumWeight = 0;
1062  fNNeutrinos = 0;
1063  fEffPOTsPerNu = 0;
1064  fAccumPOTs = 0;
1065 
1066  fGenWeighted = false;
1067  fApplyTiltWeight = true;
1068  fUseFluxAtDetCenter = 0;
1069  fDetLocIsSet = false;
1070  // by default assume user length is cm
1071  SetLengthUnits(genie::utils::units::UnitFromString("m"));
1072 
1073  this->SetDefaults();
1074  this->ResetCurrent();
1075 }
1076 //___________________________________________________________________________
1077 void GNuMIFlux::SetDefaults(void)
1078 {
1079 // - Set default neutrino species list (nue, nuebar, numu, numubar) and
1080 // maximum energy (120 GeV).
1081 // These defaults can be overwritten by user calls (at the driver init) to
1082 // GNuMIlux::SetMaxEnergy(double Ev) and
1083 // GNuMIFlux::SetFluxParticles(const PDGCodeList & particles)
1084 // - Set the default file normalization to 1E+21 POT
1085 // - Set the default flux neutrino start z position at -5m (z=0 is the
1086 // detector centre).
1087 // - Set number of cycles to 1
1088 
1089  LOG("Flux", pNOTICE) << "Setting default GNuMIFlux driver options";
1090 
1091  PDGCodeList particles;
1092  particles.push_back(kPdgNuMu);
1093  particles.push_back(kPdgAntiNuMu);
1094  particles.push_back(kPdgNuE);
1095  particles.push_back(kPdgAntiNuE);
1096 
1097  this->SetFluxParticles (particles);
1098  this->SetMaxEnergy (120./*GeV*/); // was 200, but that would be wasteful
1099  this->SetUpstreamZ (-3.4e38); // way upstream ==> use flux window
1100  this->SetNumOfCycles (0);
1101  this->SetEntryReuse (1);
1102 
1103  this->SetXMLFileBase("GNuMIFlux.xml");
1104 }
1105 //___________________________________________________________________________
1106 void GNuMIFlux::ResetCurrent(void)
1107 {
1108 // reset running values of neutrino pdg-code, 4-position & 4-momentum
1109 // and the input ntuple leaves
1110 
1111  fCurEntry->ResetCurrent();
1112  fCurEntry->ResetCopy();
1113 }
1114 //___________________________________________________________________________
1115 void GNuMIFlux::CleanUp(void)
1116 {
1117  LOG("Flux", pNOTICE) << "Cleaning up...";
1118 
1119  if (fPdgCList) delete fPdgCList;
1120  if (fPdgCListRej) delete fPdgCListRej;
1121  if (fCurEntry) delete fCurEntry;
1122 
1123  if ( fG3NuMI ) delete fG3NuMI;
1124  if ( fG4NuMI ) delete fG4NuMI;
1125  if ( fFlugg ) delete fFlugg;
1126 
1127  LOG("Flux", pNOTICE)
1128  << " flux file cycles: " << fICycle << " of " << fNCycles
1129  << ", entry " << fIEntry << " use: " << fIUse << " of " << fNUse;
1130 }
1131 
1132 //___________________________________________________________________________
1133 void GNuMIFlux::AddFile(TTree* thetree, string fname)
1134 {
1135  // Add a file to the chain
1136 
1137  ULong64_t nentries = thetree->GetEntries();
1138 
1139  // first/last "evtno" are the proton # of the first/last proton
1140  // that generated a neutrino ... not necessarily true first/last #
1141  // estimate we're probably not off by more than 100 ...
1142  // !Important change
1143  // Some files (due to some intermediate flugg issues) list evtno==0
1144  // when it isn't really true, we need to scan nearby values in case the
1145  // last entry is one of these (otherwise file contributes 0 POTs).
1146  // Also assume quantization of 500 (instead of 100).
1147 
1148  thetree->SetMakeClass(1); // need full ntuple decomposition for
1149  // the SetBranchAddress to work on g4numi ntuples. Works fine
1150  // without it on gnumi (geant3) and flugg ntuples [each with their
1151  // own slight differences] but shouldn't harm them either.
1152 
1153  Int_t evtno = 0;
1154  TBranch* br_evtno = 0;
1155  thetree->SetBranchAddress("evtno",&evtno, &br_evtno);
1156  Int_t evt_1 = 0x7FFFFFFF;
1157  Int_t evt_N = 1;
1158 #define OLDGUESS
1159 #ifdef OLDGUESS
1160  for (int j=0; j<50; ++j) {
1161  thetree->GetEntry(j);
1162  if (evtno != 0) evt_1 = TMath::Min(evtno,evt_1);
1163  thetree->GetEntry(nentries-1 -j );
1164  if (evtno != 0) evt_N = TMath::Max(evtno,evt_N);
1165  }
1166  // looks like low counts are due to "evtno" not including
1167  // protons that miss the actual target (hit baffle, etc)
1168  // this will vary with beam conditions parameters
1169  // so we should round way up, those generating flugg files
1170  // aren't going to quantize less than 1000
1171  // though 500 would probably be okay, 100 would not.
1172  const Int_t nquant = 1000; // 500; // 100
1173  const Double_t rquant = nquant;
1174 #else
1175  for (int j=0; j<50; ++j) {
1176  thetree->GetEntry(j);
1177  if (evtno != 0) evt_1 = TMath::Min(evtno,evt_1);
1178  std::cout << "[" << j << "] evtno=" << evtno << " evt_1=" << evt_1 << std::endl;
1179  }
1180  for (int j=0; j<50; ++j) {
1181  thetree->GetEntry(nentries-1 -j );
1182  if (evtno != 0) evt_N = TMath::Max(evtno,evt_N);
1183  std::cout << "[" << (nentries-1-j) << "] evtno=" << evtno << " evt_N=" << evt_N << std::endl;
1184  }
1185 
1186  Int_t nquant = 1000; // 500;
1187  Double_t rquant = nquant;
1188 #endif
1189 
1190  Int_t est_1 = (TMath::FloorNint(evt_1/rquant))*nquant + 1;
1191  Int_t est_N = (TMath::FloorNint((evt_N-1)/rquant)+1)*nquant;
1192  ULong64_t npots = est_N - est_1 + 1;
1193  LOG("Flux",pNOTICE) //INFO)
1194  << fNuFluxTreeName << "->AddFile() of " << nentries << " entries ["
1195  << evt_1 << ":" << evt_N << "%" << nquant
1196  << "(" << est_1 << ":" << est_N << ")="
1197  << npots <<" POTs] in {" << fNuFluxGen << "} file: " << fname;
1198  fNuTot += nentries;
1199  fFilePOTs += npots;
1200  fNFiles++;
1201 
1202  fNuFluxTree->AddFile(fname.c_str());
1203 }
1204 
1205 //___________________________________________________________________________
1206 void GNuMIFlux::SetLengthUnits(double user_units)
1207 {
1208  // Set the scale factor for lengths going from beam (cm) to user coordinates
1209 
1210  // GNuMIFlux uses "cm" as the length unit consistently internally (this is
1211  // the length units used by both the g3 and g4 ntuples). User interactions
1212  // setting the beam-to-detector coordinate transform, flux window, and the
1213  // returned position might need to be in other units. Use:
1214  // double scale = genie::utils::units::UnitFromString("cm");
1215  // ( #include "Utils/UnitUtils.h for declaration )
1216  // to get the correct scale factor to pass in.
1217 
1218  double rescale = fLengthUnits / user_units;
1219  fLengthUnits = user_units;
1220  double cm = genie::utils::units::UnitFromString("cm");
1221  fLengthScaleB2U = cm / user_units;
1222  fLengthScaleU2B = user_units / cm;
1223 
1224  // in case we're changing units without resetting transform/window
1225  // not recommended, but should work
1226  fCurEntry->fgX4User *= rescale;
1227  fBeamZero *= rescale;
1228  fFluxWindowPtUser[0] *= rescale;
1229  fFluxWindowPtUser[1] *= rescale;
1230  fFluxWindowPtUser[2] *= rescale;
1231 
1232  // case GNuMIFlux::kmeter: fLengthScaleB2U = 0.01 ; break;
1233  // case GNuMIFlux::kcm: fLengthScaleB2U = 1. ; break;
1234  // case GNuMIFlux::kmm: fLengthScaleB2U = 10. ; break;
1235  // case GNuMIFlux::kfm: fLengthScaleB2U = 1.e13 ; break; // 10e-2m -> 10e-15m
1236 
1237 }
1238 
1239 //___________________________________________________________________________
1240 double GNuMIFlux::LengthUnits(void) const
1241 {
1242  // Return the scale factor for lengths the user is getting
1243  double cm = genie::utils::units::UnitFromString("cm");
1244  return fLengthScaleB2U * cm ;
1245 }
1246 
1247 //___________________________________________________________________________
1248 GNuMIFluxPassThroughInfo::GNuMIFluxPassThroughInfo()
1249  : TObject()
1250 {
1251  ResetCopy();
1252  ResetCurrent();
1253 }
1254 
1255 //___________________________________________________________________________
1257 {
1258  pcodes = -1;
1259  units = -1;
1260 
1261  run = -1;
1262  evtno = -1;
1263  ndxdz = 0.;
1264  ndydz = 0.;
1265  npz = 0.;
1266  nenergy = 0.;
1267  ndxdznea = 0.;
1268  ndydznea = 0.;
1269  nenergyn = 0.;
1270  nwtnear = 0.;
1271  ndxdzfar = 0.;
1272  ndydzfar = 0.;
1273  nenergyf = 0.;
1274  nwtfar = 0.;
1275  norig = 0;
1276  ndecay = 0;
1277  ntype = 0;
1278  vx = 0.;
1279  vy = 0.;
1280  vz = 0.;
1281  pdpx = 0.;
1282  pdpy = 0.;
1283  pdpz = 0.;
1284  ppdxdz = 0.;
1285  ppdydz = 0.;
1286  pppz = 0.;
1287  ppenergy = 0.;
1288  ppmedium = 0 ;
1289  ptype = 0 ;
1290  ppvx = 0.;
1291  ppvy = 0.;
1292  ppvz = 0.;
1293  muparpx = 0.;
1294  muparpy = 0.;
1295  muparpz = 0.;
1296  mupare = 0.;
1297  necm = 0.;
1298  nimpwt = 0.;
1299  xpoint = 0.;
1300  ypoint = 0.;
1301  zpoint = 0.;
1302  tvx = 0.;
1303  tvy = 0.;
1304  tvz = 0.;
1305  tpx = 0.;
1306  tpy = 0.;
1307  tpz = 0.;
1308  tptype = 0 ;
1309  tgen = 0 ;
1310  tgptype = 0 ;
1311  tgppx = 0.;
1312  tgppy = 0.;
1313  tgppz = 0.;
1314  tprivx = 0.;
1315  tprivy = 0.;
1316  tprivz = 0.;
1317  beamx = 0.;
1318  beamy = 0.;
1319  beamz = 0.;
1320  beampx = 0.;
1321  beampy = 0.;
1322  beampz = 0.;
1323 
1324 #ifndef SKIP_MINERVA_MODS
1325  //=========================================
1326  // The following was inserted by MINERvA
1327  //=========================================
1328  ntrajectory = 0;
1329  overflow = false;
1330 
1331  for ( unsigned int itclear = 0; itclear < MAX_N_TRAJ; itclear++ ) {
1332  pdgcode[itclear] = 0;
1333  trackId[itclear] = 0;
1334  parentId[itclear] = 0;
1335  proc[itclear] = 0;
1336  ivol[itclear] = 0;
1337  fvol[itclear] = 0;
1338  startx[itclear] = 0;
1339  starty[itclear] = 0;
1340  startz[itclear] = 0;
1341  startpx[itclear] = 0;
1342  startpy[itclear] = 0;
1343  startpz[itclear] = 0;
1344  stopx[itclear] = 0;
1345  stopy[itclear] = 0;
1346  stopz[itclear] = 0;
1347  stoppx[itclear] = 0;
1348  stoppy[itclear] = 0;
1349  stoppz[itclear] = 0;
1350  pprodpx[itclear] = 0;
1351  pprodpy[itclear] = 0;
1352  pprodpz[itclear] = 0;
1353  }
1354  //END of minerva additions
1355 #endif
1356 }
1357 
1358 //___________________________________________________________________________
1360 {
1361  // reset the state of the "generated" entry
1362  fgPdgC = 0;
1363  fgXYWgt = 0;
1364  fgP4.SetPxPyPzE(0.,0.,0.,0.);
1365  fgX4.SetXYZT(0.,0.,0.,0.);
1366 }
1367 
1368 //___________________________________________________________________________
1369 /*******
1370 
1371 ALLOW DEFAULT COPY CONSTRUCTOR
1372 GNuMIFluxPassThroughInfo::GNuMIFluxPassThroughInfo(
1373  const GNuMIFluxPassThroughInfo & info) :
1374 TObject(),
1375 pcodes ( info.pcodes ),
1376 units ( info.units ),
1377 run ( info.run ),
1378 evtno ( info.evtno ),
1379 ndxdz ( info.ndxdz ),
1380 ndydz ( info.ndydz ),
1381 npz ( info.npz ),
1382 nenergy ( info.nenergy ),
1383 ndxdznea ( info.ndxdznea ),
1384 ndydznea ( info.ndydznea ),
1385 nenergyn ( info.nenergyn ),
1386 nwtnear ( info.nwtnear ),
1387 ndxdzfar ( info.ndxdzfar ),
1388 ndydzfar ( info.ndydzfar ),
1389 nenergyf ( info.nenergyf ),
1390 nwtfar ( info.nwtfar ),
1391 norig ( info.norig ),
1392 ndecay ( info.ndecay ),
1393 ntype ( info.ntype ),
1394 vx ( info.vx ),
1395 vy ( info.vy ),
1396 vz ( info.vz ),
1397 pdPx ( info.pdPx ),
1398 pdPy ( info.pdPy ),
1399 pdPz ( info.pdPz ),
1400 ppdxdz ( info.ppdxdz ),
1401 ppdydz ( info.ppdydz ),
1402 pppz ( info.pppz ),
1403 ppenergy ( info.ppenergy ),
1404 ppmedium ( info.ppmedium ),
1405 ptype ( info.ptype ),
1406 ppvx ( info.ppvx ),
1407 ppvy ( info.ppvy ),
1408 ppvz ( info.ppvz ),
1409 muparpx ( info.muparpx ),
1410 muparpy ( info.muparpy ),
1411 muparpz ( info.muparpz ),
1412 mupare ( info.mupare ),
1413 necm ( info.necm ),
1414 nimpwt ( info.nimpwt ),
1415 xpoint ( info.xpoint ),
1416 ypoint ( info.ypoint ),
1417 zpoint ( info.zpoint ),
1418 tvx ( info.tvx ),
1419 tvy ( info.tvy ),
1420 tvz ( info.tvz ),
1421 tpx ( info.tpx ),
1422 tpy ( info.tpy ),
1423 tpz ( info.tpz ),
1424 tptype ( info.tptype ),
1425 tgen ( info.tgen ),
1426 tgptype ( info.tgptype ),
1427 tgppx ( info.tgppx ),
1428 tgppy ( info.tgppy ),
1429 tgppz ( info.tgppz ),
1430 tprivx ( info.tprivx ),
1431 tprivy ( info.tprivy ),
1432 tprivz ( info.tprivz ),
1433 beamx ( info.beamx ),
1434 beamy ( info.beamy ),
1435 beamz ( info.beamz ),
1436 beampx ( info.beampx ),
1437 beampy ( info.beampy ),
1438 beampz ( info.beampz )
1439 {
1440 
1441 }
1442 *************/
1443 
1444 //___________________________________________________________________________
1446 {
1447  if ( pcodes == 0 ) {
1448  pcodes = 1; // flag that conversion has been made
1449  switch ( ntype ) {
1450  case 56: ntype = kPdgNuMu; break;
1451  case 55: ntype = kPdgAntiNuMu; break;
1452  case 53: ntype = kPdgNuE; break;
1453  case 52: ntype = kPdgAntiNuE; break;
1454  default:
1455  LOG("Flux", pNOTICE)
1456  << "ConvertPartCodes saw ntype " << ntype << " -- unknown ";
1457  }
1458  if ( ptype != 0 ) ptype = pdg::GeantToPdg(ptype);
1459  if ( tptype != 0 ) tptype = pdg::GeantToPdg(tptype);
1460  if ( tgptype != 0 ) tgptype = pdg::GeantToPdg(tgptype);
1461  } else if ( pcodes != 1 ) {
1462  // flag as unknown state ...
1463  LOG("Flux", pNOTICE)
1464  << "ConvertPartCodes saw pcodes flag as " << pcodes;
1465  }
1466 
1467 }
1468 
1469 //___________________________________________________________________________
1470 void GNuMIFluxPassThroughInfo::Print(const Option_t* /* opt */ ) const
1471 {
1472  std::cout << *this << std::endl;
1473 }
1474 
1475 //___________________________________________________________________________
1477 {
1478  run = g3->run;
1479  evtno = g3->evtno;
1480  ndxdz = g3->Ndxdz;
1481  ndydz = g3->Ndydz;
1482  npz = g3->Npz;
1483  nenergy = g3->Nenergy;
1484  ndxdznea = g3->Ndxdznea;
1485  ndydznea = g3->Ndydznea;
1486  nenergyn = g3->Nenergyn;
1487  nwtnear = g3->Nwtnear;
1488  ndxdzfar = g3->Ndxdzfar;
1489  ndydzfar = g3->Ndydzfar;
1490  nenergyf = g3->Nenergyf;
1491  nwtfar = g3->Nwtfar;
1492  norig = g3->Norig;
1493  ndecay = g3->Ndecay;
1494  ntype = g3->Ntype;
1495  vx = g3->Vx;
1496  vy = g3->Vy;
1497  vz = g3->Vz;
1498  pdpx = g3->pdpx;
1499  pdpy = g3->pdpy;
1500  pdpz = g3->pdpz;
1501  ppdxdz = g3->ppdxdz;
1502  ppdydz = g3->ppdydz;
1503  pppz = g3->pppz;
1504  ppenergy = g3->ppenergy;
1505  ppmedium = g3->ppmedium;
1506  ptype = g3->ptype;
1507  ppvx = g3->ppvx;
1508  ppvy = g3->ppvy;
1509  ppvz = g3->ppvz;
1510  muparpx = g3->muparpx;
1511  muparpy = g3->muparpy;
1512  muparpz = g3->muparpz;
1513  mupare = g3->mupare;
1514 
1515  necm = g3->Necm;
1516  nimpwt = g3->Nimpwt;
1517  xpoint = g3->xpoint;
1518  ypoint = g3->ypoint;
1519  zpoint = g3->zpoint;
1520 
1521  tvx = g3->tvx;
1522  tvy = g3->tvy;
1523  tvz = g3->tvz;
1524  tpx = g3->tpx;
1525  tpy = g3->tpy;
1526  tpz = g3->tpz;
1527  tptype = g3->tptype;
1528  tgen = g3->tgen;
1529  tgptype = g3->tgptype;
1530  tgppx = g3->tgppx;
1531  tgppy = g3->tgppy;
1532  tgppz = g3->tgppz;
1533  tprivx = g3->tprivx;
1534  tprivy = g3->tprivy;
1535  tprivz = g3->tprivz;
1536  beamx = g3->beamx;
1537  beamy = g3->beamy;
1538  beamz = g3->beamz;
1539  beampx = g3->beampx;
1540  beampy = g3->beampy;
1541  beampz = g3->beampz;
1542 
1543 }
1544 
1545 //___________________________________________________________________________
1547 {
1548 
1549  const int kNearIndx = 0;
1550  const int kFarIndx = 0;
1551 
1552  run = g4->run;
1553  evtno = g4->evtno;
1554  ndxdz = g4->Ndxdz;
1555  ndydz = g4->Ndydz;
1556  npz = g4->Npz;
1557  nenergy = g4->Nenergy;
1558  ndxdznea = g4->NdxdzNear[kNearIndx];
1559  ndydznea = g4->NdydzNear[kNearIndx];
1560  nenergyn = g4->NenergyN[kNearIndx];
1561  nwtnear = g4->NWtNear[kNearIndx];
1562  ndxdzfar = g4->NdxdzFar[kFarIndx];
1563  ndydzfar = g4->NdydzFar[kFarIndx];
1564  nenergyf = g4->NenergyF[kFarIndx];
1565  nwtfar = g4->NWtFar[kFarIndx];
1566  norig = g4->Norig;
1567  ndecay = g4->Ndecay;
1568  ntype = g4->Ntype;
1569  vx = g4->Vx;
1570  vy = g4->Vy;
1571  vz = g4->Vz;
1572  pdpx = g4->pdPx;
1573  pdpy = g4->pdPy;
1574  pdpz = g4->pdPz;
1575  ppdxdz = g4->ppdxdz;
1576  ppdydz = g4->ppdydz;
1577  pppz = g4->pppz;
1578  ppenergy = g4->ppenergy;
1579  ppmedium = (Int_t)g4->ppmedium; // int in g3, double in g4!
1580  ptype = g4->ptype;
1581  ppvx = g4->ppvx;
1582  ppvy = g4->ppvy;
1583  ppvz = g4->ppvz;
1584  muparpx = g4->muparpx;
1585  muparpy = g4->muparpy;
1586  muparpz = g4->muparpz;
1587  mupare = g4->mupare;
1588 
1589  necm = g4->Necm;
1590  nimpwt = g4->Nimpwt;
1591  xpoint = g4->xpoint;
1592  ypoint = g4->ypoint;
1593  zpoint = g4->zpoint;
1594 
1595  tvx = g4->tvx;
1596  tvy = g4->tvy;
1597  tvz = g4->tvz;
1598  tpx = g4->tpx;
1599  tpy = g4->tpy;
1600  tpz = g4->tpz;
1601  tptype = g4->tptype;
1602  tgen = g4->tgen;
1603  tgptype = 0 ; // no equivalent in g4
1604  tgppx = 0.;
1605  tgppy = 0.;
1606  tgppz = 0.;
1607  tprivx = 0.;
1608  tprivy = 0.;
1609  tprivz = 0.;
1610  beamx = 0.; // g4->protonX;
1611  beamy = 0.; // g4->protonY;
1612  beamz = 0.; // g4->protonZ;
1613  beampx = 0.; // g4->protonPx;
1614  beampy = 0.; // g4->protonPy;
1615  beampz = 0.; // g4->protonPz;
1616 
1617 #ifndef SKIP_MINERVA_MODS
1618  //=========================================
1619  // The following was inserted by MINERvA
1620  //=========================================
1621  ntrajectory = g4->ntrajectory;
1622  overflow = g4->overflow;
1623  // Limit number of trajectories to MAX_N_TRAJ
1624  Int_t ntraj = ntrajectory;
1625  if ( overflow )
1626  ntraj = MAX_N_TRAJ;
1627 
1628  for ( Int_t ic = 0; ic < ntraj; ++ic ) {
1629  pdgcode[ic] = g4->pdg[ic];
1630  trackId[ic] = g4->trackId[ic];
1631  parentId[ic] = g4->parentId[ic];
1632 
1633  startx[ic] = g4->startx[ic];
1634  starty[ic] = g4->starty[ic];
1635  startz[ic] = g4->startz[ic];
1636  stopx[ic] = g4->stopx[ic];
1637  stopy[ic] = g4->stopy[ic];
1638  stopz[ic] = g4->stopz[ic];
1639  startpx[ic] = g4->startpx[ic];
1640  startpy[ic] = g4->startpy[ic];
1641  startpz[ic] = g4->startpz[ic];
1642  stoppx[ic] = g4->stoppx[ic];
1643  stoppy[ic] = g4->stoppy[ic];
1644  stoppz[ic] = g4->stoppz[ic];
1645  pprodpx[ic] = g4->pprodpx[ic];
1646  pprodpy[ic] = g4->pprodpy[ic];
1647  pprodpz[ic] = g4->pprodpz[ic];
1648 
1649  proc[ic] = getProcessID(g4->proc[ic]);
1650  ivol[ic] = getVolID(g4->ivol[ic]);
1651  fvol[ic] = getVolID(g4->fvol[ic]);
1652  }
1653  //END of minerva additions
1654 #endif
1655 
1656 }
1657 
1658 //___________________________________________________________________________
1660 {
1661  run = f->run;
1662  evtno = f->evtno;
1663  ndxdz = f->Ndxdz;
1664  ndydz = f->Ndydz;
1665  npz = f->Npz;
1666  nenergy = f->Nenergy;
1667  ndxdznea = f->Ndxdznea;
1668  ndydznea = f->Ndydznea;
1669  nenergyn = f->Nenergyn;
1670  nwtnear = f->Nwtnear;
1671  ndxdzfar = f->Ndxdzfar;
1672  ndydzfar = f->Ndydzfar;
1673  nenergyf = f->Nenergyf;
1674  nwtfar = f->Nwtfar;
1675  norig = f->Norig;
1676  ndecay = f->Ndecay;
1677  ntype = f->Ntype;
1678  vx = f->Vx;
1679  vy = f->Vy;
1680  vz = f->Vz;
1681  pdpx = f->pdPx;
1682  pdpy = f->pdPy;
1683  pdpz = f->pdPz;
1684  ppdxdz = f->ppdxdz;
1685  ppdydz = f->ppdydz;
1686  pppz = f->pppz;
1687  ppenergy = f->ppenergy;
1688  ppmedium = f->ppmedium;
1689  ptype = f->ptype;
1690  ppvx = f->ppvx;
1691  ppvy = f->ppvy;
1692  ppvz = f->ppvz;
1693  muparpx = f->muparpx;
1694  muparpy = f->muparpy;
1695  muparpz = f->muparpz;
1696  mupare = f->mupare;
1697 
1698  necm = f->Necm;
1699  nimpwt = f->Nimpwt;
1700  xpoint = f->xpoint;
1701  ypoint = f->ypoint;
1702  zpoint = f->zpoint;
1703 
1704  tvx = f->tvx;
1705  tvy = f->tvy;
1706  tvz = f->tvz;
1707  tpx = f->tpx;
1708  tpy = f->tpy;
1709  tpz = f->tpz;
1710  tptype = f->tptype;
1711  tgen = f->tgen;
1712  tgptype = f->tgptype;
1713  tgppx = f->tgppx;
1714  tgppy = f->tgppy;
1715  tgppz = f->tgppz;
1716  tprivx = f->tprivx;
1717  tprivy = f->tprivy;
1718  tprivz = f->tprivz;
1719  beamx = f->beamx;
1720  beamy = f->beamy;
1721  beamz = f->beamz;
1722  beampx = f->beampx;
1723  beampy = f->beampy;
1724  beampz = f->beampz;
1725 
1726 }
1727 
1728 //___________________________________________________________________________
1729 int GNuMIFluxPassThroughInfo::CalcEnuWgt(const TLorentzVector& xyz,
1730  double& enu, double& wgt_xy) const
1731 {
1732 
1733  // Neutrino Energy and Weigth at arbitrary point
1734  // based on:
1735  // NuMI-NOTE-BEAM-0109 (MINOS DocDB # 109)
1736  // Title: Neutrino Beam Simulation using PAW with Weighted Monte Carlos
1737  // Author: Rick Milburn
1738  // Date: 1995-10-01
1739 
1740  // history:
1741  // jzh 3/21/96 grab R.H.Milburn's weighing routine
1742  // jzh 5/ 9/96 substantially modify the weighting function use dot product
1743  // instead of rotation vecs to get theta get all info except
1744  // det from ADAMO banks neutrino parent is in Particle.inc
1745  // Add weighting factor for polarized muon decay
1746  // jzh 4/17/97 convert more code to double precision because of problems
1747  // with Enu>30 GeV
1748  // rwh 10/ 9/08 transliterate function from f77 to C++
1749 
1750  // original function description:
1751  // Real function for use with PAW Ntuple To transform from destination
1752  // detector geometry to the unit sphere moving with decaying hadron with
1753  // velocity v, BETA=v/c, etc.. For (pseudo)scalar hadrons the decays will
1754  // be isotropic in this sphere so the fractional area (out of 4-pi) is the
1755  // fraction of decays that hit the target. For a given target point and
1756  // area, and given x-y components of decay transverse location and slope,
1757  // and given decay distance from target ans given decay GAMMA and
1758  // rest-frame neutrino energy, the lab energy at the target and the
1759  // fractional solid angle in the rest-frame are determined.
1760  // For muon decays, correction for non-isotropic nature of decay is done.
1761 
1762  // Arguments:
1763  // double x, y, z :: position to evaluate for (enu,wgt_xy)
1764  // in *beam* frame coordinates (cm units)
1765  // double enu, wgt_xy :: resulting energy and weight
1766  // Return:
1767  // int :: error code
1768  // Assumptions:
1769  // Energies given in GeV
1770  // Particle codes have been translated from GEANT into PDG codes
1771 
1772  // for now ... these _should_ come from DB
1773  // but use these hard-coded values to "exactly" reproduce old code
1774  //
1775  const double kPIMASS = 0.13957;
1776  const double kKMASS = 0.49368;
1777  const double kK0MASS = 0.49767;
1778  const double kMUMASS = 0.105658389;
1779  const double kOMEGAMASS = 1.67245;
1780 
1781  const int kpdg_nue = 12; // extended Geant 53
1782  const int kpdg_nuebar = -12; // extended Geant 52
1783  const int kpdg_numu = 14; // extended Geant 56
1784  const int kpdg_numubar = -14; // extended Geant 55
1785 
1786  const int kpdg_muplus = -13; // Geant 5
1787  const int kpdg_muminus = 13; // Geant 6
1788  const int kpdg_pionplus = 211; // Geant 8
1789  const int kpdg_pionminus = -211; // Geant 9
1790  const int kpdg_k0long = 130; // Geant 10 ( K0=311, K0S=310 )
1791  const int kpdg_k0short = 310; // Geant 16
1792  const int kpdg_k0mix = 311;
1793  const int kpdg_kaonplus = 321; // Geant 11
1794  const int kpdg_kaonminus = -321; // Geant 12
1795  const int kpdg_omegaminus = 3334; // Geant 24
1796  const int kpdg_omegaplus = -3334; // Geant 32
1797 
1798  const double kRDET = 100.0; // set to flux per 100 cm radius
1799 
1800  double xpos = xyz.X();
1801  double ypos = xyz.Y();
1802  double zpos = xyz.Z();
1803 
1804  enu = 0.0; // don't know what the final value is
1805  wgt_xy = 0.0; // but set these in case we return early due to error
1806 
1807 
1808  // in principle we should get these from the particle DB
1809  // but for consistency testing use the hardcoded values
1810  double parent_mass = kPIMASS;
1811  switch ( this->ptype ) {
1812  case kpdg_pionplus:
1813  case kpdg_pionminus:
1814  parent_mass = kPIMASS;
1815  break;
1816  case kpdg_kaonplus:
1817  case kpdg_kaonminus:
1818  parent_mass = kKMASS;
1819  break;
1820  case kpdg_k0long:
1821  case kpdg_k0short:
1822  case kpdg_k0mix:
1823  parent_mass = kK0MASS;
1824  break;
1825  case kpdg_muplus:
1826  case kpdg_muminus:
1827  parent_mass = kMUMASS;
1828  break;
1829  case kpdg_omegaminus:
1830  case kpdg_omegaplus:
1831  parent_mass = kOMEGAMASS;
1832  break;
1833  default:
1834  std::cerr << "NU_REWGT unknown particle type " << this->ptype
1835  << std::endl << std::flush;
1836  LOG("Flux",pFATAL) << "NU_REWGT unknown particle type " << this->ptype;
1837  assert(0);
1838  return 1;
1839  }
1840 
1841  double parentp2 = ( this->pdpx*this->pdpx +
1842  this->pdpy*this->pdpy +
1843  this->pdpz*this->pdpz );
1844  double parent_energy = TMath::Sqrt( parentp2 +
1845  parent_mass*parent_mass);
1846  double parentp = TMath::Sqrt( parentp2 );
1847 
1848  double gamma = parent_energy / parent_mass;
1849  double gamma_sqr = gamma * gamma;
1850  double beta_mag = TMath::Sqrt( ( gamma_sqr - 1.0 )/gamma_sqr );
1851 
1852  // Get the neutrino energy in the parent decay CM
1853  double enuzr = this->necm;
1854  // Get angle from parent line of flight to chosen point in beam frame
1855  double rad = TMath::Sqrt( (xpos-this->vx)*(xpos-this->vx) +
1856  (ypos-this->vy)*(ypos-this->vy) +
1857  (zpos-this->vz)*(zpos-this->vz) );
1858 
1859  double emrat = 1.0;
1860  double costh_pardet = -999., theta_pardet = -999.;
1861 
1862  // boost correction, but only if parent hasn't stopped
1863  if ( parentp > 0. ) {
1864  costh_pardet = ( this->pdpx*(xpos-this->vx) +
1865  this->pdpy*(ypos-this->vy) +
1866  this->pdpz*(zpos-this->vz) )
1867  / ( parentp * rad);
1868  if ( costh_pardet > 1.0 ) costh_pardet = 1.0;
1869  if ( costh_pardet < -1.0 ) costh_pardet = -1.0;
1870  theta_pardet = TMath::ACos(costh_pardet);
1871 
1872  // Weighted neutrino energy in beam, approx, good for small theta
1873  emrat = 1.0 / ( gamma * ( 1.0 - beta_mag * costh_pardet ));
1874  }
1875 
1876  enu = emrat * enuzr; // the energy ... normally
1877 
1878  // RWH-debug
1879  bool debug = false;
1880  if (debug) {
1881  std::cout << std::setprecision(15);
1882  std::cout << "xyz (" << xpos << "," << ypos << "," << zpos << ")" << std::endl;
1883  std::cout << "ptype " << this->ptype << " m " << parent_mass
1884  << " p " << parentp << " e " << parent_energy << " gamma " << gamma
1885  << " beta " << beta_mag << std::endl;
1886 
1887  std::cout << " enuzr " << enuzr << " rad " << rad << " costh " << costh_pardet
1888  << " theta " << theta_pardet << " emrat " << emrat
1889  << " enu " << enu << std::endl;
1890  }
1891 
1892 #ifdef GNUMI_TEST_XY_WGT
1893  gpartials.xdet = xpos;
1894  gpartials.ydet = ypos;
1895  gpartials.zdet = zpos;
1896  gpartials.parent_mass = parent_mass;
1897  gpartials.parentp = parentp;
1898  gpartials.parent_energy = parent_energy;
1899  gpartials.gamma = gamma;
1900  gpartials.beta_mag = beta_mag;
1901  gpartials.enuzr = enuzr;
1902  gpartials.rad = rad;
1903  gpartials.costh_pardet = costh_pardet;
1904  gpartials.theta_pardet = theta_pardet;
1905  gpartials.emrat = emrat;
1906  gpartials.eneu = enu;
1907 #endif
1908 
1909  // Get solid angle/4pi for detector element
1910  // small angle approximation, fixed by Alex Radovic
1911  //SAA//double sangdet = ( kRDET*kRDET / ( (zpos-this->vz)*(zpos-this->vz)))/4.0;
1912 
1913  double sanddetcomp = TMath::Sqrt(( (xpos-this->vx)*(xpos-this->vx) ) +
1914  ( (ypos-this->vy)*(ypos-this->vy) ) +
1915  ( (zpos-this->vz)*(zpos-this->vz) ) );
1916  double sangdet = ( 1.0 - TMath::Cos(TMath::ATan( kRDET / sanddetcomp)))/2.0;
1917 
1918  // Weight for solid angle and lorentz boost
1919  wgt_xy = sangdet * ( emrat * emrat ); // ! the weight ... normally
1920 
1921 #ifdef GNUMI_TEST_XY_WGT
1922  gpartials.sangdet = sangdet;
1923  gpartials.wgt = wgt_xy;
1924  gpartials.ptype = this->ptype; // assume already PDG
1925 #endif
1926 
1927  // Done for all except polarized muon decay
1928  // in which case need to modify weight
1929  // (must be done in double precision)
1930  if ( this->ptype == kpdg_muplus || this->ptype == kpdg_muminus) {
1931  double beta[3], p_dcm_nu[4], p_nu[3], p_pcm_mp[3], partial;
1932 
1933  // Boost neu neutrino to mu decay CM
1934  beta[0] = this->pdpx / parent_energy;
1935  beta[1] = this->pdpy / parent_energy;
1936  beta[2] = this->pdpz / parent_energy;
1937  p_nu[0] = (xpos-this->vx)*enu/rad;
1938  p_nu[1] = (ypos-this->vy)*enu/rad;
1939  p_nu[2] = (zpos-this->vz)*enu/rad;
1940  partial = gamma *
1941  (beta[0]*p_nu[0] + beta[1]*p_nu[1] + beta[2]*p_nu[2] );
1942  partial = enu - partial/(gamma+1.0);
1943  // the following calculation is numerically imprecise
1944  // especially p_dcm_nu[2] leads to taking the difference of numbers of order ~10's
1945  // and getting results of order ~0.02's
1946  // for g3numi we're starting with floats (ie. good to ~1 part in 10^7)
1947  p_dcm_nu[0] = p_nu[0] - beta[0]*gamma*partial;
1948  p_dcm_nu[1] = p_nu[1] - beta[1]*gamma*partial;
1949  p_dcm_nu[2] = p_nu[2] - beta[2]*gamma*partial;
1950  p_dcm_nu[3] = TMath::Sqrt( p_dcm_nu[0]*p_dcm_nu[0] +
1951  p_dcm_nu[1]*p_dcm_nu[1] +
1952  p_dcm_nu[2]*p_dcm_nu[2] );
1953 
1954 #ifdef GNUMI_TEST_XY_WGT
1955  gpartials.betanu[0] = beta[0];
1956  gpartials.betanu[1] = beta[1];
1957  gpartials.betanu[2] = beta[2];
1958  gpartials.p_nu[0] = p_nu[0];
1959  gpartials.p_nu[1] = p_nu[1];
1960  gpartials.p_nu[2] = p_nu[2];
1961  gpartials.partial1 = partial;
1962  gpartials.p_dcm_nu[0] = p_dcm_nu[0];
1963  gpartials.p_dcm_nu[1] = p_dcm_nu[1];
1964  gpartials.p_dcm_nu[2] = p_dcm_nu[2];
1965  gpartials.p_dcm_nu[3] = p_dcm_nu[3];
1966 #endif
1967 
1968  // Boost parent of mu to mu production CM
1969  double particle_energy = this->ppenergy;
1970  gamma = particle_energy/parent_mass;
1971  beta[0] = this->ppdxdz * this->pppz / particle_energy;
1972  beta[1] = this->ppdydz * this->pppz / particle_energy;
1973  beta[2] = this->pppz / particle_energy;
1974  partial = gamma * ( beta[0]*this->muparpx +
1975  beta[1]*this->muparpy +
1976  beta[2]*this->muparpz );
1977  partial = this->mupare - partial/(gamma+1.0);
1978  p_pcm_mp[0] = this->muparpx - beta[0]*gamma*partial;
1979  p_pcm_mp[1] = this->muparpy - beta[1]*gamma*partial;
1980  p_pcm_mp[2] = this->muparpz - beta[2]*gamma*partial;
1981  double p_pcm = TMath::Sqrt ( p_pcm_mp[0]*p_pcm_mp[0] +
1982  p_pcm_mp[1]*p_pcm_mp[1] +
1983  p_pcm_mp[2]*p_pcm_mp[2] );
1984 
1985  //std::cout << " muparpxyz " << this->muparpx << " "
1986  // << this->muparpy << " " << this->muparpz << std::endl;
1987  //std::cout << " beta " << beta[0] << " " << beta[1] << " " << beta[2] << std::endl;
1988  //std::cout << " gamma " << gamma << " partial " << partial << std::endl;
1989  //std::cout << " p_pcm_mp " << p_pcm_mp[0] << " " << p_pcm_mp[1] << " "
1990  // << p_pcm_mp[2] << " " << p_pcm << std::endl;
1991 
1992 #ifdef GNUMI_TEST_XY_WGT
1993  gpartials.muparent_px = this->muparpx;
1994  gpartials.muparent_py = this->muparpy;
1995  gpartials.muparent_pz = this->muparpz;
1996  gpartials.gammamp = gamma;
1997  gpartials.betamp[0] = beta[0];
1998  gpartials.betamp[1] = beta[1];
1999  gpartials.betamp[2] = beta[2];
2000  gpartials.partial2 = partial;
2001  gpartials.p_pcm_mp[0] = p_pcm_mp[0];
2002  gpartials.p_pcm_mp[1] = p_pcm_mp[1];
2003  gpartials.p_pcm_mp[2] = p_pcm_mp[2];
2004  gpartials.p_pcm = p_pcm;
2005 #endif
2006 
2007  const double eps = 1.0e-30; // ? what value to use
2008  if ( p_pcm < eps || p_dcm_nu[3] < eps ) {
2009  return 3; // mu missing parent info?
2010  }
2011  // Calc new decay angle w.r.t. (anti)spin direction
2012  double costh = ( p_dcm_nu[0]*p_pcm_mp[0] +
2013  p_dcm_nu[1]*p_pcm_mp[1] +
2014  p_dcm_nu[2]*p_pcm_mp[2] ) /
2015  ( p_dcm_nu[3]*p_pcm );
2016  if ( costh > 1.0 ) costh = 1.0;
2017  if ( costh < -1.0 ) costh = -1.0;
2018  // Calc relative weight due to angle difference
2019  double wgt_ratio = 0.0;
2020  switch ( this->ntype ) {
2021  case kpdg_nue:
2022  case kpdg_nuebar:
2023  wgt_ratio = 1.0 - costh;
2024  break;
2025  case kpdg_numu:
2026  case kpdg_numubar:
2027  {
2028  double xnu = 2.0 * enuzr / kMUMASS;
2029  wgt_ratio = ( (3.0-2.0*xnu ) - (1.0-2.0*xnu)*costh ) / (3.0-2.0*xnu);
2030  break;
2031  }
2032  default:
2033  return 2; // bad neutrino type
2034  }
2035  wgt_xy = wgt_xy * wgt_ratio;
2036 
2037 #ifdef GNUMI_TEST_XY_WGT
2038  gpartials.ntype = this->ntype; // assume converted to PDG
2039  gpartials.costhmu = costh;
2040  gpartials.wgt_ratio = wgt_ratio;
2041 #endif
2042 
2043  } // ptype is muon
2044 
2045  return 0;
2046 }
2047 
2048 //___________________________________________________________________________
2049 
2050 
2051 namespace genie {
2052 namespace flux {
2053  ostream & operator << (
2054  ostream & stream, const genie::flux::GNuMIFluxPassThroughInfo & info)
2055  {
2056  // stream << "\n ndecay = " << info.ndecay << std::endl;
2057  stream << "\nGNuMIFlux run " << info.run << " evtno " << info.evtno
2058  << " (pcodes " << info.pcodes << " units " << info.units << ")"
2059  << "\n random dk: dx/dz " << info.ndxdz
2060  << " dy/dz " << info.ndydz
2061  << " pz " << info.npz << " E " << info.nenergy
2062  << "\n near00 dk: dx/dz " << info.ndxdznea
2063  << " dy/dz " << info.ndydznea
2064  << " E " << info.nenergyn << " wgt " << info.nwtnear
2065  << "\n far00 dk: dx/dz " << info.ndxdzfar
2066  << " dy/dz " << info.ndydzfar
2067  << " E " << info.nenergyf << " wgt " << info.nwtfar
2068  << "\n norig " << info.norig << " ndecay " << info.ndecay
2069  << " ntype " << info.ntype
2070  << "\n had vtx " << info.vx << " " << info.vy << " " << info.vz
2071  << "\n parent p3 @ dk " << info.pdpx << " " << info.pdpy << " " << info.pdpz
2072  << "\n parent prod: dx/dz " << info.ppdxdz
2073  << " dy/dz " << info.ppdydz
2074  << " pz " << info.pppz << " E " << info.ppenergy
2075  << "\n ppmedium " << info.ppmedium << " ptype " << info.ptype
2076  << " ppvtx " << info.ppvx << " " << info.ppvy << " " << info.ppvz
2077  << "\n mu parent p4 " << info.muparpx << " " << info.muparpy
2078  << " " << info.muparpz << " " << info.mupare
2079  << "\n necm " << info.necm << " nimpwt " << info.nimpwt
2080  << "\n point x,y,z " << info.xpoint << " " << info.ypoint
2081  << " " << info.zpoint
2082  << "\n tv x,y,z " << info.tvx << " " << info.tvy << " " << info.tvz
2083  << "\n tptype " << info.tptype << " tgen " << info.tgen
2084  << " tgptype " << info.tgptype
2085  << "\n tgp px,py,pz " << info.tgppx << " " << info.tgppy
2086  << " " << info.tgppz
2087  << "\n tpriv x,y,z " << info.tprivx << " " << info.tprivy
2088  << " " << info.tprivz
2089  << "\n beam x,y,z " << info.beamx << " " << info.beamy
2090  << " " << info.beamz
2091  << "\n beam px,py,pz " << info.beampx << " " << info.beampy
2092  << " " << info.beampz
2093  ;
2094 
2095 #ifndef SKIP_MINERVA_MODS
2096  //=========================================
2097  // The following was inserted by MINERvA
2098  //=========================================
2099  stream << "\nNeutrino History : ntrajectories " << info.ntrajectory
2100  << "\n (trkID, pdg) of nu parents: [";
2101  Int_t ntraj = info.ntrajectory;
2102  if ( info.overflow ) ntraj = GNuMIFluxPassThroughInfo::MAX_N_TRAJ;
2103 
2104  for ( Int_t itt = 0; itt < ntraj; ++itt )
2105  stream << "(" << info.trackId[itt-1] << ", " << info.pdgcode[itt] << ") ";
2106  stream << "]\n";
2107  //END of minerva additions
2108 #endif
2109 
2110  stream << "\nCurrent: pdg " << info.fgPdgC
2111  << " xywgt " << info.fgXYWgt
2112  << "\n p4 (beam): " << utils::print::P4AsShortString(&info.fgP4)
2113  << "\n x4 (beam): " << utils::print::X4AsString(&info.fgX4)
2114  << "\n p4 (user): " << utils::print::P4AsShortString(&info.fgP4User)
2115  << "\n x4 (user): " << utils::print::X4AsString(&info.fgX4User);
2116  ;
2117 #ifdef GNUMI_TEST_XY_WGT
2118  stream << "\n" << xypartials::GetStaticInstance();
2119 #endif
2120 
2121 
2122  /*
2123  //std::cout << "GNuMIFlux::PrintCurrent ....." << std::endl;
2124  //LOG("Flux", pINFO)
2125  LOG("Flux", pNOTICE)
2126  << "Current Leaf Values: "
2127  << " run " << fLf_run << " evtno " << fLf_evtno << "\n"
2128  << " NenergyN " << fLf_NenergyN << " NWtNear " << fLf_NWtNear
2129  << " NenergyF " << fLf_NenergyF << " NWtFar " << fLf_NWtFar << "\n"
2130  << " Norig " << fLf_Norig << " Ndecay " << fLf_Ndecay << " Ntype " << fLf_Ntype << "\n"
2131  << " Vxyz " << fLf_Vx << " " << fLf_Vy << " " << fLf_Vz
2132  << " pdPxyz " << fLf_pdPx << " " << fLf_pdPy << " " << fLf_pdPz << "\n"
2133  << " pp dxdz " << fLf_ppdxdz << " dydz " << fLf_ppdydz << " pz " << fLf_pppz << "\n"
2134  << " pp energy " << fLf_ppenergy << " medium " << fLf_ppmedium
2135  << " ptype " << fLf_ptype
2136  << " ppvxyz " << fLf_ppvx << " " << fLf_ppvy << " " << fLf_ppvz << "\n"
2137  << " muparpxyz " << fLf_muparpx << " " << fLf_muparpy << " " << fLf_muparpz
2138  << " mupare " << fLf_mupare << "\n"
2139  << ;
2140  */
2141 
2142  return stream;
2143  }
2144 }//flux
2145 }//genie
2146 
2147 //___________________________________________________________________________
2148 
2149 #ifdef GNUMI_TEST_XY_WGT
2150 
2151 double fabserr(double a, double b)
2152 { return TMath::Abs(a-b)/TMath::Max(TMath::Abs(b),1.0e-30); }
2153 
2154 int outdiff(double a, double b, double eps, const char* label)
2155 {
2156  double err = fabserr(a,b);
2157  if ( err > eps ) {
2158  std::cout << std::setw(15) << label << " err " << err
2159  << " vals " << a << " " << b << std::endl;
2160  return 1;
2161  }
2162  return 0;
2163 }
2164 
2165 int gnumi2pdg(int igeant)
2166 {
2167  switch ( igeant ) {
2168  case 52: return -12; // nuebar
2169  case 53: return 12; // nue
2170  case 55: return -14; // numubar
2171  case 56: return 14; // numu
2172  case 58: return -16; // nutaubar
2173  case 59: return 16; // nutau
2174  default:
2175  return genie::pdg::GeantToPdg(igeant);
2176  }
2177 }
2178 
2179 void xypartials::ReadStream(std::ifstream& myfile)
2180 {
2181  myfile >> parent_mass >> parentp >> parent_energy;
2182  myfile >> gamma >> beta_mag >> enuzr >> rad;
2183  myfile >> costh_pardet >> theta_pardet >> emrat >> eneu;
2184  myfile >> sangdet >> wgt;
2185  int ptype_g;
2186  myfile >> ptype_g;
2187  ptype = gnumi2pdg(ptype_g);
2188  if ( ptype == 13 || ptype == -13 ) {
2189  //std::cout << "ReadStream saw ptype " << ptype << std::endl;
2190  myfile >> betanu[0] >> betanu[1] >> betanu[2]
2191  >> p_nu[0] >> p_nu[1] >> p_nu[2];
2192  myfile >> partial1
2193  >> p_dcm_nu[0] >> p_dcm_nu[1] >> p_dcm_nu[2] >> p_dcm_nu[3];
2194 
2195  myfile >> muparent_px >> muparent_py >> muparent_pz;
2196  myfile >> gammamp >> betamp[0] >> betamp[1] >> betamp[2];
2197  myfile >> partial2
2198  >> p_pcm_mp[0] >> p_pcm_mp[1] >> p_pcm_mp[2] >> p_pcm;
2199 
2200  if ( p_pcm != 0.0 && p_dcm_nu[3] != 0.0 ) {
2201  int ntype_g;
2202  myfile >> costhmu >> ntype_g >> wgt_ratio;
2203  ntype = gnumi2pdg(ntype_g);
2204  }
2205  }
2206 }
2207 
2208 int xypartials::Compare(const xypartials& other) const
2209 {
2210  double eps1 = 2.5e-5; // 0.003; // 6.0e-4; // 2.5e-4;
2211  double eps2 = 2.5e-5; // 6.0e-4; // 2.5e-4;
2212  double epsX = 2.5e-5; // 2.5e-4;
2213  int np = 0;
2214  np += outdiff(parent_mass ,other.parent_mass ,eps1,"parent_mass");
2215  np += outdiff(parentp ,other.parentp ,eps1,"parentp");
2216  np += outdiff(parent_energy,other.parent_energy,eps1,"parent_energy");
2217  np += outdiff(gamma ,other.gamma ,eps1,"gamma");
2218  np += outdiff(beta_mag ,other.beta_mag ,eps1,"beta_mag");
2219  np += outdiff(enuzr ,other.enuzr ,eps1,"enuzr");
2220  np += outdiff(rad ,other.rad ,eps1,"rad");
2221  np += outdiff(costh_pardet ,other.costh_pardet ,eps1,"costh_pardet");
2222  //np += outdiff(theta_pardet ,other.theta_pardet ,eps1,"theta_pardet");
2223  np += outdiff(emrat ,other.emrat ,eps1,"emrat");
2224  np += outdiff(eneu ,other.eneu ,epsX,"eneu");
2225  np += outdiff(sangdet ,other.sangdet ,eps1,"sangdet");
2226  np += outdiff(wgt ,other.wgt ,epsX,"wgt");
2227  if ( ptype != other.ptype ) {
2228  std::cout << "ptype mismatch " << ptype << " " << other.ptype << std::endl;
2229  np++;
2230  }
2231  if ( TMath::Abs(ptype)==13 || TMath::Abs(other.ptype)==13 ) {
2232  //std::cout << "== ismu " << std::endl;
2233  np += outdiff(betanu[0] ,other.betanu[0] ,eps2,"betanu[0]");
2234  np += outdiff(betanu[1] ,other.betanu[1] ,eps2,"betanu[1]");
2235  np += outdiff(betanu[2] ,other.betanu[2] ,eps2,"betanu[2]");
2236  np += outdiff(p_nu[0] ,other.p_nu[0] ,eps2,"p_nu[0]");
2237  np += outdiff(p_nu[1] ,other.p_nu[1] ,eps2,"p_nu[1]");
2238  np += outdiff(p_nu[2] ,other.p_nu[2] ,eps2,"p_nu[2]");
2239  np += outdiff(partial1 ,other.partial1 ,eps2,"partial1");
2240  np += outdiff(p_dcm_nu[0] ,other.p_dcm_nu[0] ,eps2,"p_dcm_nu[0]");
2241  np += outdiff(p_dcm_nu[1] ,other.p_dcm_nu[1] ,eps2,"p_dcm_nu[1]");
2242  np += outdiff(p_dcm_nu[2] ,other.p_dcm_nu[2] ,eps2,"p_dcm_nu[2]");
2243  np += outdiff(p_dcm_nu[3] ,other.p_dcm_nu[3] ,eps2,"p_dcm_nu[3]");
2244 
2245  np += outdiff(muparent_px ,other.muparent_px ,eps2,"muparent_px");
2246  np += outdiff(muparent_py ,other.muparent_py ,eps2,"muparent_py");
2247  np += outdiff(muparent_pz ,other.muparent_pz ,eps2,"muparent_pz");
2248  np += outdiff(gammamp ,other.gammamp ,eps1,"gammamp");
2249  np += outdiff(betamp[0] ,other.betamp[0] ,eps1,"betamp[0]");
2250  np += outdiff(betamp[1] ,other.betamp[1] ,eps1,"betamp[1]");
2251  np += outdiff(betamp[2] ,other.betamp[2] ,eps1,"betamp[2]");
2252  np += outdiff(partial2 ,other.partial2 ,eps1,"partial2");
2253  np += outdiff(p_pcm_mp[0] ,other.p_pcm_mp[0] ,eps1,"p_pcm_mp[0]");
2254  np += outdiff(p_pcm_mp[1] ,other.p_pcm_mp[1] ,eps1,"p_pcm_mp[1]");
2255  np += outdiff(p_pcm_mp[2] ,other.p_pcm_mp[2] ,eps1,"p_pcm_mp[2]");
2256  np += outdiff(p_pcm ,other.p_pcm ,eps1,"p_pcm");
2257 
2258  if ( ntype != other.ntype ) {
2259  std::cout << "ntype mismatch " << ntype << " " << other.ntype << std::endl;
2260  np++;
2261  }
2262  np += outdiff(costhmu ,other.costhmu ,eps1,"costhmu");
2263  np += outdiff(wgt_ratio ,other.wgt_ratio ,eps1,"wgt_ratio");
2264 
2265  }
2266  return np;
2267 }
2268 
2269 void xypartials::Print(const Option_t* /* opt */) const
2270 {
2271  std::cout << *this << std::endl;
2272 }
2273 
2274 namespace genie {
2275 namespace flux {
2276  ostream & operator << (ostream & stream, const genie::flux::xypartials & xyp )
2277  {
2278  stream << "GNuMIFlux xypartials " << std::endl;
2279  stream << " xyzdet (" << xyp.xdet << "," << xyp.ydet << ","
2280  << xyp.zdet << ")" << std::endl;
2281  stream << " parent: mass=" << xyp.parent_mass << " p=" << xyp.parentp
2282  << " e=" << xyp.parent_energy << " gamma=" << xyp.gamma
2283  << " beta_mag=" << xyp.beta_mag << std::endl;
2284  stream << " enuzr=" << xyp.enuzr << " rad=" << xyp.rad
2285  << " costh_pardet=" << xyp.costh_pardet << std::endl;
2286  stream << " emrat=" << xyp.emrat << " sangdet=" << xyp.sangdet
2287  << " wgt=" << xyp.wgt << std::endl;
2288  stream << " ptype=" << xyp.ptype << " "
2289  << ((TMath::Abs(xyp.ptype) == 13)?"is-muon":"not-muon")
2290  << std::endl;
2291 
2292  if ( TMath::Abs(xyp.ptype)==13 ) {
2293  stream << " betanu: [" << xyp.betanu[0] << "," << xyp.betanu[1]
2294  << "," << xyp.betanu[2] << "]" << std::endl;
2295  stream << " p_nu: [" << xyp.p_nu[0] << "," << xyp.p_nu[1]
2296  << "," << xyp.p_nu[2] << "]" << std::endl;
2297  stream << " partial1=" << xyp.partial1 << std::endl;
2298  stream << " p_dcm_nu: [" << xyp.p_dcm_nu[0] << "," << xyp.p_dcm_nu[1]
2299  << "," << xyp.p_dcm_nu[2]
2300  << "," << xyp.p_dcm_nu[3] << "]" << std::endl;
2301  stream << " muparent_p: [" << xyp.muparent_px << "," << xyp.muparent_py
2302  << "," << xyp.muparent_pz << "]" << std::endl;
2303  stream << " gammamp=" << xyp.gammamp << std::endl;
2304  stream << " betamp: [" << xyp.betamp[0] << "," << xyp.betamp[1] << ","
2305  << xyp.betamp[2] << "]" << std::endl;
2306  stream << " partial2=" << xyp.partial2 << std::endl;
2307  stream << " p_pcm_mp: [" << xyp.p_pcm_mp[0] << "," << xyp.p_pcm_mp[1]
2308  << "," << xyp.p_pcm_mp[2] << "] p_pcm="
2309  << xyp.p_pcm << std::endl;
2310  stream << " ntype=" << xyp.ntype
2311  << " costhmu=" << xyp.costhmu
2312  << " wgt_ratio=" << xyp.wgt_ratio << std::endl;
2313  }
2314  return stream;
2315  }
2316 } // flux
2317 } // genie
2318 
2319 xypartials& xypartials::GetStaticInstance()
2320 { return gpartials; }
2321 #endif
2322 //___________________________________________________________________________
2323 
2324 bool GNuMIFlux::LoadConfig(string cfg)
2325 {
2326  const char* altxml = gSystem->Getenv("GNUMIFLUXXML");
2327  if ( altxml ) {
2328  SetXMLFileBase(altxml);
2329  }
2330  genie::flux::GNuMIFluxXMLHelper helper(this);
2331  return helper.LoadConfig(cfg);
2332 }
2333 
2334 //___________________________________________________________________________
2335 
2337 {
2338 
2339  std::ostringstream s;
2340  PDGCodeList::const_iterator itr = fPdgCList->begin();
2341  for ( ; itr != fPdgCList->end(); ++itr) s << (*itr) << " ";
2342  s << "[rejected: ";
2343  itr = fPdgCListRej->begin();
2344  for ( ; itr != fPdgCListRej->end(); ++itr) s << (*itr) << " ";
2345  s << " ] ";
2346 
2347  std::ostringstream fpattout;
2348  for (size_t i = 0; i < fNuFluxFilePatterns.size(); ++i)
2349  fpattout << "\n [" << std::setw(3) << i << "] " << fNuFluxFilePatterns[i];
2350 
2351  std::ostringstream flistout;
2352  std::vector<std::string> flist = GetFileList();
2353  for (size_t i = 0; i < flist.size(); ++i)
2354  flistout << "\n [" << std::setw(3) << i << "] " << flist[i];
2355 
2356  TLorentzVector usr0(0,0,0,0);
2357  TLorentzVector usr0asbeam;
2358  User2BeamPos(usr0,usr0asbeam);
2359 
2360  const int w=10, p=6;
2361  std::ostringstream beamrot_str, beamrotinv_str;
2362  beamrot_str
2363  << "fBeamRot: " << std::setprecision(p) << "\n"
2364  << " [ "
2365  << std::setw(w) << fBeamRot.XX() << " "
2366  << std::setw(w) << fBeamRot.XY() << " "
2367  << std::setw(w) << fBeamRot.XZ() << " ]\n"
2368  << " [ "
2369  << std::setw(w) << fBeamRot.YX() << " "
2370  << std::setw(w) << fBeamRot.YY() << " "
2371  << std::setw(w) << fBeamRot.YZ() << " ]\n"
2372  << " [ "
2373  << std::setw(w) << fBeamRot.ZX() << " "
2374  << std::setw(w) << fBeamRot.ZY() << " "
2375  << std::setw(w) << fBeamRot.ZZ() << " ]";
2376  beamrotinv_str
2377  << "fBeamRotInv: " << std::setprecision(p) << "\n"
2378  << " [ "
2379  << std::setw(w) << fBeamRotInv.XX() << " "
2380  << std::setw(w) << fBeamRotInv.XY() << " "
2381  << std::setw(w) << fBeamRotInv.XZ() << " ]\n"
2382  << " [ "
2383  << std::setw(w) << fBeamRotInv.YX() << " "
2384  << std::setw(w) << fBeamRotInv.YY() << " "
2385  << std::setw(w) << fBeamRotInv.YZ() << " ]\n"
2386  << " [ "
2387  << std::setw(w) << fBeamRotInv.ZX() << " "
2388  << std::setw(w) << fBeamRotInv.ZY() << " "
2389  << std::setw(w) << fBeamRotInv.ZZ() << " ]";
2390 
2391  LOG("Flux", pNOTICE)
2392  << "GNuMIFlux Config:"
2393  << "\n Enu_max " << fMaxEv
2394  << "\n pdg-codes: " << s.str() << "\n "
2395  << (fG3NuMI?"g3numi":"")
2396  << (fG4NuMI?"g4numi":"")
2397  << (fFlugg?"flugg":"")
2398  << "/" << fNuFluxGen << " "
2399  << "(" << fNuFluxTreeName << "), " << fNEntries << " entries"
2400  << " (FilePOTs " << fFilePOTs << ") "
2401  << "in " << fNFiles << " files: "
2402  << flistout.str()
2403  << "\n from file patterns:"
2404  << fpattout.str()
2405  << "\n wgt max=" << fMaxWeight << " fudge=" << fMaxWgtFudge << " using "
2406  << fMaxWgtEntries << " entries"
2407  << "\n Z0 pushback " << fZ0
2408  << "\n used entry " << fIEntry << " " << fIUse << "/" << fNUse
2409  << " times, in " << fICycle << "/" << fNCycles << " cycles"
2410  << "\n SumWeight " << fSumWeight << " for " << fNNeutrinos << " neutrinos"
2411  << "\n EffPOTsPerNu " << fEffPOTsPerNu << " AccumPOTs " << fAccumPOTs
2412  << "\n GenWeighted \"" << (fGenWeighted?"true":"false") << ", "
2413  << "\", Detector location set \"" << (fDetLocIsSet?"true":"false") << "\""
2414  << "\n LengthUnits " << fLengthUnits << ", scale b2u " << fLengthScaleB2U
2415  << ", scale u2b " << fLengthScaleU2B
2416  << "\n User Flux Window: "
2420  << "\n Flux Window (cm, beam coord): "
2421  << "\n base " << utils::print::X4AsString(&fFluxWindowBase)
2422  << "\n dir1 " << utils::print::X4AsString(&fFluxWindowDir1) << " len " << fFluxWindowLen1
2423  << "\n dir2 " << utils::print::X4AsString(&fFluxWindowDir2) << " len " << fFluxWindowLen2
2424  << "\n normal " << utils::print::Vec3AsString(&(fWindowNormal))
2425  << "\n User Beam Origin: "
2426  << "\n base " << utils::print::X4AsString(&fBeamZero)
2427  << "\n " << beamrot_str.str() << " "
2428  << "\n Detector Origin (beam coord): "
2429  << "\n base " << utils::print::X4AsString(&usr0asbeam)
2430  << "\n " << beamrotinv_str.str() << " "
2431  << "\n UseFluxAtDetCenter " << fUseFluxAtDetCenter;
2432 
2433 }
2434 
2435 //___________________________________________________________________________
2436 std::vector<std::string> GNuMIFlux::GetFileList()
2437 {
2438  std::vector<std::string> flist;
2439  TObjArray *fileElements=fNuFluxTree->GetListOfFiles();
2440  TIter next(fileElements);
2441  TChainElement *chEl=0;
2442  while (( chEl=(TChainElement*)next() )) {
2443  flist.push_back(chEl->GetTitle());
2444  }
2445  return flist;
2446 }
2447 
2448 //___________________________________________________________________________
2449 
2450 std::vector<double> GNuMIFluxXMLHelper::GetDoubleVector(std::string str)
2451 {
2452  // turn string into vector<double>
2453  // be liberal about separators, users might punctuate for clarity
2454  std::vector<std::string> strtokens = genie::utils::str::Split(str," ,;:()[]=");
2455  std::vector<double> vect;
2456  size_t ntok = strtokens.size();
2457 
2458  if ( fVerbose > 2 )
2459  std::cout << "GetDoubleVector \"" << str << "\"" << std::endl;
2460 
2461  for (size_t i=0; i < ntok; ++i) {
2462  std::string trimmed = utils::str::TrimSpaces(strtokens[i]);
2463  if ( " " == trimmed || "" == trimmed ) continue; // skip empty strings
2464  double val = strtod(trimmed.c_str(), (char**)NULL);
2465  if ( fVerbose > 2 )
2466  std::cout << "(" << vect.size() << ") = " << val << std::endl;
2467  vect.push_back(val);
2468  }
2469 
2470  return vect;
2471 }
2472 
2473 std::vector<long int> GNuMIFluxXMLHelper::GetIntVector(std::string str)
2474 {
2475  // turn string into vector<long int>
2476  // be liberal about separators, users might punctuate for clarity
2477  std::vector<std::string> strtokens = genie::utils::str::Split(str," ,;:()[]=");
2478  std::vector<long int> vect;
2479  size_t ntok = strtokens.size();
2480 
2481  if ( fVerbose > 2 )
2482  std::cout << "GetIntVector \"" << str << "\"" << std::endl;
2483 
2484  for (size_t i=0; i < ntok; ++i) {
2485  std::string trimmed = utils::str::TrimSpaces(strtokens[i]);
2486  if ( " " == trimmed || "" == trimmed ) continue; // skip empty strings
2487  long int val = strtol(trimmed.c_str(),(char**)NULL,10);
2488  if ( fVerbose > 2 )
2489  std::cout << "(" << vect.size() << ") = " << val << std::endl;
2490  vect.push_back(val);
2491  }
2492 
2493  return vect;
2494 }
2495 
2497 {
2499 
2500  bool is_accessible = ! (gSystem->AccessPathName(fname.c_str()));
2501  if (!is_accessible) {
2502  SLOG("GNuMIFlux", pERROR)
2503  << "The XML doc doesn't exist! (filename: " << fname << ")";
2504  return false;
2505  }
2506 
2507  xmlDocPtr xml_doc = xmlParseFile( fname.c_str() );
2508  if ( xml_doc == NULL) {
2509  SLOG("GNuMIFlux", pERROR)
2510  << "The XML doc can't be parsed! (filename: " << fname << ")";
2511  return false;
2512  }
2513 
2514  xmlNodePtr xml_root = xmlDocGetRootElement( xml_doc );
2515  if ( xml_root == NULL ) {
2516  SLOG("GNuMIFlux", pERROR)
2517  << "The XML doc is empty! (filename: " << fname << ")";
2518  return false;
2519  }
2520 
2521  string rootele = "gnumi_config";
2522  if ( xmlStrcmp(xml_root->name, (const xmlChar*)rootele.c_str() ) ) {
2523  SLOG("GNuMIFlux", pERROR)
2524  << "The XML doc has invalid root element! (filename: " << fname << ")"
2525  << " expected \"" << rootele << "\", saw \"" << xml_root->name << "\"";
2526  return false;
2527  }
2528 
2529  SLOG("GNuMIFlux", pINFO) << "Attempt to load config \"" << cfg
2530  << "\" from file: " << fname;
2531 
2532  bool found = this->LoadParamSet(xml_doc,cfg);
2533 
2534  xmlFree(xml_doc);
2535  return found;
2536 
2537 }
2538 
2539 bool GNuMIFluxXMLHelper::LoadParamSet(xmlDocPtr& xml_doc, string cfg)
2540 {
2541 
2542  xmlNodePtr xml_root = xmlDocGetRootElement( xml_doc );
2543 
2544  // loop over all xml tree nodes that are children of the root node
2545  // read the entries looking for "param_set" of the right name
2546 
2547  // loop looking for particular config
2548  bool found = false;
2549  xmlNodePtr xml_pset = xml_root->xmlChildrenNode;
2550  for ( ; xml_pset != NULL ; xml_pset = xml_pset->next ) {
2551  if ( ! xmlStrEqual(xml_pset->name, (const xmlChar*)"param_set") ) continue;
2552  // every time there is a 'param_set' tag
2553  string param_set_name =
2555 
2556  if ( param_set_name != cfg ) continue;
2557 
2558  SLOG("GNuMIFlux", pINFO) << "Found config \"" << cfg;
2559 
2560  this->ParseParamSet(xml_doc,xml_pset);
2561  found = true;
2562 
2563  } // loop over elements of root
2564  xmlFree(xml_pset);
2565 
2566  return found;
2567 }
2568 
2569 void GNuMIFluxXMLHelper::ParseParamSet(xmlDocPtr& xml_doc, xmlNodePtr& xml_pset)
2570 {
2571  xmlNodePtr xml_child = xml_pset->xmlChildrenNode;
2572  for ( ; xml_child != NULL ; xml_child = xml_child->next ) {
2573  // handle basic gnumi_config/param_set
2574  // bad cast away const on next line, but function sig requires it
2575  string pname =
2576  utils::xml::TrimSpaces(const_cast<xmlChar*>(xml_child->name));
2577  if ( pname == "text" || pname == "comment" ) continue;
2578  string pval =
2580  xmlNodeListGetString(xml_doc, xml_child->xmlChildrenNode, 1));
2581 
2582  if ( fVerbose > 1 )
2583  SLOG("GNuMIFlux", pINFO)
2584  << " pname \"" << pname << "\", string value \"" << pval << "\"";
2585 
2586  if ( pname == "verbose" ) {
2587  fVerbose = atoi(pval.c_str());
2588 
2589  } else if ( pname == "using_param_set" ) {
2590  SLOG("GNuMIFlux", pWARN) << "start using_param_set: \"" << pval << "\"";
2591  bool found = this->LoadParamSet(xml_doc,pval); // recurse
2592  if ( ! found ) {
2593  SLOG("GNuMIFlux", pFATAL) << "using_param_set: \"" << pval << "\" NOT FOUND";
2594  assert(found);
2595  }
2596  SLOG("GNuMIFlux", pWARN) << "done using_param_set: \"" << pval << "\"";
2597  } else if ( pname == "units" ) {
2598  double scale = genie::utils::units::UnitFromString(pval);
2599  fGNuMI->SetLengthUnits(scale);
2600  SLOG("GNuMIFlux", pINFO) << "set user units to \"" << pval << "\"";
2601 
2602  } else if ( pname == "beamdir" ) {
2603  ParseBeamDir(xml_doc,xml_child);
2605 
2606  } else if ( pname == "beampos" ) {
2607  ParseBeamPos(pval);
2609 
2610  } else if ( pname == "window" ) {
2611  ParseWindowSeries(xml_doc,xml_child);
2612  // RWH !!!! MEMORY LEAK!!!!
2613  //std::cout << " flux window " << std::endl
2614  // << " [0] " << utils::print::X4AsString(new TLorentzVector(fFluxWindowPt[0],0)) << std::endl
2615  // << " [1] " << utils::print::X4AsString(new TLorentzVector(fFluxWindowPt[1],0)) << std::endl
2616  // << " [2] " << utils::print::X4AsString(new TLorentzVector(fFluxWindowPt[2],0)) << std::endl;
2617 
2619  fFluxWindowPtXML[1],
2620  fFluxWindowPtXML[2]);
2621 
2622  } else if ( pname == "enumax" ) {
2623  ParseEnuMax(pval);
2624 
2625  } else if ( pname == "upstreamz" ) {
2626  double z0usr = -3.4e38;
2627  std::vector<double> v = GetDoubleVector(pval);
2628  if ( v.size() > 0 ) z0usr = v[0];
2629  fGNuMI->SetUpstreamZ(z0usr);
2630  SLOG("GNuMIFlux", pINFO) << "set upstreamz = " << z0usr;
2631 
2632  } else if ( pname == "reuse" ) {
2633  long int nreuse = 1;
2634  std::vector<long int> v = GetIntVector(pval);
2635  if ( v.size() > 0 ) nreuse = v[0];
2636  fGNuMI->SetEntryReuse(nreuse);
2637  SLOG("GNuMIFlux", pINFO) << "set entry reuse = " << nreuse;
2638 
2639  } else {
2640  SLOG("GNuMIFlux", pWARN)
2641  << " NOT HANDLED: pname \"" << pname
2642  << "\", string value \"" << pval << "\"";
2643 
2644  }
2645 
2646  } // loop over param_set contents
2647  xmlFree(xml_child);
2648 }
2649 
2650 void GNuMIFluxXMLHelper::ParseBeamDir(xmlDocPtr& xml_doc, xmlNodePtr& xml_beamdir)
2651 {
2652  fBeamRotXML.SetToIdentity(); // start fresh
2653 
2654  string dirtype =
2656  utils::xml::GetAttribute(xml_beamdir,"type"));
2657 
2658  string pval =
2660  xmlNodeListGetString(xml_doc, xml_beamdir->xmlChildrenNode, 1));
2661 
2662  if ( dirtype == "series" ) {
2663  // series of rotations around an axis
2664  ParseRotSeries(xml_doc,xml_beamdir);
2665 
2666  } else if ( dirtype == "thetaphi3") {
2667  // G3 style triplet of (theta,phi) pairs
2668  std::vector<double> thetaphi3 = GetDoubleVector(pval);
2669  string units =
2670  utils::str::TrimSpaces(utils::xml::GetAttribute(xml_beamdir,"units"));
2671  if ( thetaphi3.size() == 6 ) {
2672  TRotation fTempRot;
2673  TVector3 newX = AnglesToAxis(thetaphi3[0],thetaphi3[1],units);
2674  TVector3 newY = AnglesToAxis(thetaphi3[2],thetaphi3[3],units);
2675  TVector3 newZ = AnglesToAxis(thetaphi3[4],thetaphi3[5],units);
2676  fTempRot.RotateAxes(newX,newY,newZ);
2677  fBeamRotXML = fTempRot; //.Inverse();
2678  } else {
2679  SLOG("GNuMIFlux", pWARN)
2680  << " type=\"" << dirtype << "\" within <beamdir> needs 6 values";
2681  }
2682 
2683  } else if ( dirtype == "newxyz" ) {
2684  // G4 style new axis values
2685  std::vector<double> newdir = GetDoubleVector(pval);
2686  if ( newdir.size() == 9 ) {
2687  TRotation fTempRot;
2688  TVector3 newX = TVector3(newdir[0],newdir[1],newdir[2]).Unit();
2689  TVector3 newY = TVector3(newdir[3],newdir[4],newdir[5]).Unit();
2690  TVector3 newZ = TVector3(newdir[6],newdir[7],newdir[8]).Unit();
2691  fTempRot.RotateAxes(newX,newY,newZ);
2692  fBeamRotXML = fTempRot.Inverse(); // weirdly necessary: frame vs. obj rot
2693  } else {
2694  SLOG("GNuMIFlux", pWARN)
2695  << " type=\"" << dirtype << "\" within <beamdir> needs 9 values";
2696  }
2697 
2698  } else {
2699  // yet something else ... what? 3 choices weren't sufficient?
2700  SLOG("GNuMIFlux", pWARN)
2701  << " UNHANDLED type=\"" << dirtype << "\" within <beamdir>";
2702  }
2703 
2704  if ( fVerbose > 1 ) {
2705  int w=10, p=6;
2706  std::cout << " fBeamRotXML: " << std::setprecision(p) << std::endl;
2707  std::cout << " [ "
2708  << std::setw(w) << fBeamRotXML.XX() << " "
2709  << std::setw(w) << fBeamRotXML.XY() << " "
2710  << std::setw(w) << fBeamRotXML.XZ() << endl
2711  << " "
2712  << std::setw(w) << fBeamRotXML.YX() << " "
2713  << std::setw(w) << fBeamRotXML.YY() << " "
2714  << std::setw(w) << fBeamRotXML.YZ() << endl
2715  << " "
2716  << std::setw(w) << fBeamRotXML.ZX() << " "
2717  << std::setw(w) << fBeamRotXML.ZY() << " "
2718  << std::setw(w) << fBeamRotXML.ZZ() << " ] " << std::endl;
2719  std::cout << std::endl;
2720  }
2721 
2722 }
2723 
2725 {
2726  std::vector<double> xyz = GetDoubleVector(str);
2727  if ( xyz.size() == 3 ) {
2728  fBeamPosXML = TVector3(xyz[0],xyz[1],xyz[2]);
2729  } else if ( xyz.size() == 6 ) {
2730  // should check for '=' between triplets but we won't be so pedantic
2731  // ( userx, usery, userz ) = ( beamx, beamy, beamz )
2732  TVector3 userpos(xyz[0],xyz[1],xyz[2]);
2733  TVector3 beampos(xyz[3],xyz[4],xyz[5]);
2734  fBeamPosXML = userpos - fBeamRotXML*beampos;
2735  } else {
2736  SLOG("GNuMIFlux", pWARN)
2737  << "Unable to parse " << xyz.size() << " values in <beampos>";
2738  return;
2739  }
2740  if ( fVerbose > 1 ) {
2741  int w=16, p=10;
2742  std::cout << " fBeamPosXML: [ " << std::setprecision(p)
2743  << std::setw(w) << fBeamPosXML.X() << " , "
2744  << std::setw(w) << fBeamPosXML.Y() << " , "
2745  << std::setw(w) << fBeamPosXML.Z() << " ] "
2746  << std::endl;
2747  }
2748 }
2749 
2750 void GNuMIFluxXMLHelper::ParseEnuMax(std::string str)
2751 {
2752  std::vector<double> v = GetDoubleVector(str);
2753  size_t n = v.size();
2754  if ( n > 0 ) {
2755  fGNuMI->SetMaxEnergy(v[0]);
2756  if ( fVerbose > 1 )
2757  std::cout << "ParseEnuMax SetMaxEnergy(" << v[0] << ") " << std::endl;
2758  }
2759  if ( n > 1 ) {
2760  fGNuMI->SetMaxEFudge(v[1]);
2761  if ( fVerbose > 1 )
2762  std::cout << "ParseEnuMax SetMaxEFudge(" << v[1] << ")" << std::endl;
2763  }
2764  if ( n > 2 ) {
2765  if ( n == 3 ) {
2766  fGNuMI->SetMaxWgtScan(v[2]);
2767  if ( fVerbose > 1 )
2768  std::cout << "ParseEnuMax SetMaxWgtScan(" << v[2] << ")" << std::endl;
2769  } else {
2770  long int nentries = (long int)v[3];
2771  fGNuMI->SetMaxWgtScan(v[2],nentries);
2772  if ( fVerbose > 1 )
2773  std::cout << "ParseEnuMax SetMaxWgtScan(" << v[2] << "," << nentries << ")" << std::endl;
2774  }
2775  }
2776 }
2777 
2778 void GNuMIFluxXMLHelper::ParseRotSeries(xmlDocPtr& xml_doc, xmlNodePtr& xml_pset)
2779 {
2780  TRotation fTempRot; // reset matrix
2781 
2782  xmlNodePtr xml_child = xml_pset->xmlChildrenNode;
2783  for ( ; xml_child != NULL ; xml_child = xml_child->next ) {
2784  // in a <beamdir> of type "series"
2785  // should be a sequence of <rotation> entries
2786  string name =
2787  utils::xml::TrimSpaces(const_cast<xmlChar*>(xml_child->name));
2788  if ( name == "text" || name == "comment" ) continue;
2789 
2790  if ( name == "rotation" ) {
2791  string val = utils::xml::TrimSpaces(
2792  xmlNodeListGetString(xml_doc, xml_child->xmlChildrenNode, 1));
2793  string axis =
2795 
2796  string units =
2798 
2799  double rot = atof(val.c_str());
2800  // assume radians unless given a hint that it's degrees
2801  if ( 'd' == units[0] || 'D' == units[0] ) rot *= TMath::DegToRad();
2802 
2803  if ( fVerbose > 0 )
2804  SLOG("GNuMIFlux", pINFO)
2805  << " rotate " << rot << " radians around " << axis << " axis";
2806 
2807  if ( axis[0] == 'x' || axis[0] == 'X' ) fTempRot.RotateX(rot);
2808  else if ( axis[0] == 'y' || axis[0] == 'Y' ) fTempRot.RotateY(rot);
2809  else if ( axis[0] == 'z' || axis[0] == 'Z' ) fTempRot.RotateZ(rot);
2810  else {
2811  SLOG("GNuMIFlux", pINFO)
2812  << " no " << axis << " to rotate around";
2813  }
2814 
2815  } else {
2816  SLOG("GNuMIFlux", pWARN)
2817  << " found <" << name << "> within <beamdir type=\"series\">";
2818  }
2819  }
2820  // TRotation rotates objects not frames, so we want the inverse
2821  fBeamRotXML = fTempRot.Inverse();
2822  xmlFree(xml_child);
2823 }
2824 
2825 void GNuMIFluxXMLHelper::ParseWindowSeries(xmlDocPtr& xml_doc, xmlNodePtr& xml_pset)
2826 {
2827  int ientry = -1;
2828 
2829  xmlNodePtr xml_child = xml_pset->xmlChildrenNode;
2830  for ( ; xml_child != NULL ; xml_child = xml_child->next ) {
2831  // in a <windowr> element
2832  // should be a sequence of <point> entries
2833  string name =
2834  utils::xml::TrimSpaces(const_cast<xmlChar*>(xml_child->name));
2835  if ( name == "text" || name == "comment" ) continue;
2836 
2837  if ( name == "point" ) {
2838  string val =
2840  xmlNodeListGetString(xml_doc, xml_child->xmlChildrenNode, 1));
2841  string coord =
2843 
2844  std::vector<double> xyz = GetDoubleVector(val);
2845  if ( xyz.size() != 3 || coord != "det" ) {
2846  SLOG("GNuMIFlux", pWARN)
2847  << "parsing <window> found <point> but size=" << xyz.size()
2848  << " (expect 3) and coord=\"" << coord << "\" (expect \"det\")"
2849  << " IGNORE problem";
2850  }
2851  ++ientry;
2852  if ( ientry < 3 && ientry >= 0 ) {
2853  TVector3 pt(xyz[0],xyz[1],xyz[2]);
2854  if ( fVerbose > 0 ) {
2855  int w=16, p=10;
2856  std::cout << " point[" << ientry <<"] = [ " << std::setprecision(p)
2857  << std::setw(w) << pt.X() << " , "
2858  << std::setw(w) << pt.Y() << " , "
2859  << std::setw(w) << pt.Z() << " ] "
2860  << std::endl;
2861  }
2862  fFluxWindowPtXML[ientry] = pt; // save the point
2863  } else {
2864  SLOG("GNuMIFlux", pWARN)
2865  << " <window><point> ientry " << ientry << " out of range (0-2)";
2866  }
2867 
2868  } else {
2869  SLOG("GNuMIFlux", pWARN)
2870  << " found <" << name << "> within <window>";
2871  }
2872  }
2873  xmlFree(xml_child);
2874 }
2875 
2876 TVector3 GNuMIFluxXMLHelper::AnglesToAxis(double theta, double phi, std::string units)
2877 {
2878  double xyz[3];
2879  // assume radians unless given a hint that it's degrees
2880  double scale = ('d'==units[0]||'D'==units[0]) ? TMath::DegToRad() : 1.0 ;
2881 
2882  xyz[0] = TMath::Cos(scale*phi)*TMath::Sin(scale*theta);
2883  xyz[1] = TMath::Sin(scale*phi)*TMath::Sin(scale*theta);
2884  xyz[2] = TMath::Cos(scale*theta);
2885  // condition vector to eliminate most floating point errors
2886  for (int i=0; i<3; ++i) {
2887  const double eps = 1.0e-15;
2888  if (TMath::Abs(xyz[i]) < eps ) xyz[i] = 0;
2889  if (TMath::Abs(xyz[i]-1) < eps ) xyz[i] = 1;
2890  if (TMath::Abs(xyz[i]+1) < eps ) xyz[i] = -1;
2891  }
2892  return TVector3(xyz[0],xyz[1],xyz[2]);
2893 }
2894 
2895 TVector3 GNuMIFluxXMLHelper::ParseTV3(const string& str)
2896 {
2897  std::vector<double> xyz = GetDoubleVector(str);
2898  if ( xyz.size() != 3 ) {
2899  return TVector3();
2900  SLOG("GNuMIFlux", pWARN)
2901  << " ParseTV3 \"" << str << "\" had " << xyz.size() << " elements ";
2902  }
2903  return TVector3(xyz[0],xyz[1],xyz[2]);
2904 
2905 }
2906 //___________________________________________________________________________
2907 
2908 #ifndef SKIP_MINERVA_MODS
2909 //=========================================
2910 // The following was inserted by MINERvA
2911 //=========================================
2913  int ival=0;
2914  if(sval=="AddedLV")ival=1;
2915  else if(sval=="AlTube1LV")ival=2;
2916  else if(sval=="AlTube2LV")ival=3;
2917  else if(sval=="Al_BLK1")ival=4;
2918  else if(sval=="Al_BLK2")ival=5;
2919  else if(sval=="Al_BLK3")ival=6;
2920  else if(sval=="Al_BLK4")ival=7;
2921  else if(sval=="Al_BLK5")ival=8;
2922  else if(sval=="Al_BLK6")ival=9;
2923  else if(sval=="Al_BLK7")ival=10;
2924  else if(sval=="Al_BLK8")ival=11;
2925  else if(sval=="AlholeL")ival=12;
2926  else if(sval=="AlholeR")ival=13;
2927  else if(sval=="BEndLV")ival=14;
2928  else if(sval=="BFrontLV")ival=15;
2929  else if(sval=="BeDWLV")ival=16;
2930  else if(sval=="BeUp1LV")ival=17;
2931  else if(sval=="BeUp2LV")ival=18;
2932  else if(sval=="BeUp3LV")ival=19;
2933  else if(sval=="BodyLV")ival=20;
2934  else if(sval=="BudalMonitor")ival=21;
2935  else if(sval=="CLid1LV")ival=22;
2936  else if(sval=="CLid2LV")ival=23;
2937  else if(sval=="CShld_BLK11")ival=24;
2938  else if(sval=="CShld_BLK12")ival=25;
2939  else if(sval=="CShld_BLK2")ival=26;
2940  else if(sval=="CShld_BLK3")ival=27;
2941  else if(sval=="CShld_BLK4")ival=28;
2942  else if(sval=="CShld_BLK7")ival=29;
2943  else if(sval=="CShld_BLK8")ival=30;
2944  else if(sval=="CShld_stl,BLK")ival=31;
2945  else if(sval=="CerTubeLV")ival=32;
2946  else if(sval=="CeramicRod")ival=33;
2947  else if(sval=="ConcShield")ival=34;
2948  else if(sval=="Concrete Chase Section")ival=35;
2949  else if(sval=="Conn1LV")ival=36;
2950  else if(sval=="Conn2LV")ival=37;
2951  else if(sval=="Conn3LV")ival=38;
2952  else if(sval=="DNWN")ival=39;
2953  else if(sval=="DPIP")ival=40;
2954  else if(sval=="DVOL")ival=41;
2955  else if(sval=="DuratekBlock")ival=42;
2956  else if(sval=="DuratekBlockCovering")ival=43;
2957  else if(sval=="HadCell")ival=44;
2958  else if(sval=="HadronAbsorber")ival=45;
2959  else if(sval=="MuMonAlcvFill_0")ival=46;
2960  else if(sval=="MuMonAlcv_0")ival=47;
2961  else if(sval=="MuMonAlcv_1")ival=48;
2962  else if(sval=="MuMonAlcv_2")ival=49;
2963  else if(sval=="MuMon_0")ival=50;
2964  else if(sval=="MuMon_1")ival=51;
2965  else if(sval=="MuMon_2")ival=52;
2966  else if(sval=="PHorn1CPB1slv")ival=53;
2967  else if(sval=="PHorn1CPB2slv")ival=54;
2968  else if(sval=="PHorn1F")ival=55;
2969  else if(sval=="PHorn1Front")ival=56;
2970  else if(sval=="PHorn1IC")ival=57;
2971  else if(sval=="PHorn1InsRingslv")ival=58;
2972  else if(sval=="PHorn1OC")ival=59;
2973  else if(sval=="PHorn2CPB1slv")ival=60;
2974  else if(sval=="PHorn2CPB2slv")ival=61;
2975  else if(sval=="PHorn2F")ival=62;
2976  else if(sval=="PHorn2Front")ival=63;
2977  else if(sval=="PHorn2IC")ival=64;
2978  else if(sval=="PHorn2InsRingslv")ival=65;
2979  else if(sval=="PHorn2OC")ival=66;
2980  else if(sval=="PVHadMon")ival=67;
2981  else if(sval=="Pipe1")ival=68;
2982  else if(sval=="Pipe1_water")ival=69;
2983  else if(sval=="Pipe2")ival=70;
2984  else if(sval=="Pipe2_water")ival=71;
2985  else if(sval=="Pipe3")ival=72;
2986  else if(sval=="Pipe3_water")ival=73;
2987  else if(sval=="Pipe4")ival=74;
2988  else if(sval=="Pipe4_water")ival=75;
2989  else if(sval=="Pipe5")ival=76;
2990  else if(sval=="Pipe5_water")ival=77;
2991  else if(sval=="Pipe6")ival=78;
2992  else if(sval=="Pipe6_water")ival=79;
2993  else if(sval=="Pipe7")ival=80;
2994  else if(sval=="Pipe8")ival=81;
2995  else if(sval=="Pipe8_water")ival=82;
2996  else if(sval=="Pipe9")ival=83;
2997  else if(sval=="PipeAdapter1")ival=84;
2998  else if(sval=="PipeAdapter1_water")ival=85;
2999  else if(sval=="PipeAdapter2")ival=86;
3000  else if(sval=="PipeAdapter2_water")ival=87;
3001  else if(sval=="PipeBellowB")ival=88;
3002  else if(sval=="PipeBellowB_water")ival=89;
3003  else if(sval=="PipeBellowT")ival=90;
3004  else if(sval=="PipeBellowT_water")ival=91;
3005  else if(sval=="PipeC1")ival=92;
3006  else if(sval=="PipeC1_water")ival=93;
3007  else if(sval=="PipeC2")ival=94;
3008  else if(sval=="PipeC2_water")ival=95;
3009  else if(sval=="ROCK")ival=96;
3010  else if(sval=="Ring1LV")ival=97;
3011  else if(sval=="Ring2LV")ival=98;
3012  else if(sval=="Ring3LV")ival=99;
3013  else if(sval=="Ring4LV")ival=100;
3014  else if(sval=="Ring5LV")ival=101;
3015  else if(sval=="SC01")ival=102;
3016  else if(sval=="SpiderSupport")ival=103;
3017  else if(sval=="Stl_BLK1")ival=104;
3018  else if(sval=="Stl_BLK10")ival=105;
3019  else if(sval=="Stl_BLK2")ival=106;
3020  else if(sval=="Stl_BLK3")ival=107;
3021  else if(sval=="Stl_BLK4")ival=108;
3022  else if(sval=="Stl_BLK5")ival=109;
3023  else if(sval=="Stl_BLK6")ival=110;
3024  else if(sval=="Stl_BLK7")ival=111;
3025  else if(sval=="Stl_BLK8")ival=112;
3026  else if(sval=="Stl_BLK9")ival=113;
3027  else if(sval=="Stlhole")ival=114;
3028  else if(sval=="TGAR")ival=115;
3029  else if(sval=="TGT1")ival=116;
3030  else if(sval=="TGTExitCyl2LV")ival=117;
3031  else if(sval=="TUNE")ival=118;
3032  else if(sval=="Tube1aLV")ival=119;
3033  else if(sval=="Tube1bLV")ival=120;
3034  else if(sval=="UpWn1")ival=121;
3035  else if(sval=="UpWn2")ival=122;
3036  else if(sval=="UpWnAl1SLV")ival=123;
3037  else if(sval=="UpWnAl2SLV")ival=124;
3038  else if(sval=="UpWnAl3SLV")ival=125;
3039  else if(sval=="UpWnFe1SLV")ival=126;
3040  else if(sval=="UpWnFe2SLV")ival=127;
3041  else if(sval=="UpWnPolyCone")ival=128;
3042  else if(sval=="blu_BLK25")ival=129;
3043  else if(sval=="blu_BLK26")ival=130;
3044  else if(sval=="blu_BLK27")ival=131;
3045  else if(sval=="blu_BLK28")ival=132;
3046  else if(sval=="blu_BLK29")ival=133;
3047  else if(sval=="blu_BLK32")ival=134;
3048  else if(sval=="blu_BLK37")ival=135;
3049  else if(sval=="blu_BLK38")ival=136;
3050  else if(sval=="blu_BLK39")ival=137;
3051  else if(sval=="blu_BLK40")ival=138;
3052  else if(sval=="blu_BLK45")ival=139;
3053  else if(sval=="blu_BLK46")ival=140;
3054  else if(sval=="blu_BLK47")ival=141;
3055  else if(sval=="blu_BLK48")ival=142;
3056  else if(sval=="blu_BLK49")ival=143;
3057  else if(sval=="blu_BLK50")ival=144;
3058  else if(sval=="blu_BLK51")ival=145;
3059  else if(sval=="blu_BLK53")ival=146;
3060  else if(sval=="blu_BLK55")ival=147;
3061  else if(sval=="blu_BLK57")ival=148;
3062  else if(sval=="blu_BLK59")ival=149;
3063  else if(sval=="blu_BLK61")ival=150;
3064  else if(sval=="blu_BLK63")ival=151;
3065  else if(sval=="blu_BLK64")ival=152;
3066  else if(sval=="blu_BLK65")ival=153;
3067  else if(sval=="blu_BLK66")ival=154;
3068  else if(sval=="blu_BLK67")ival=155;
3069  else if(sval=="blu_BLK68")ival=156;
3070  else if(sval=="blu_BLK69")ival=157;
3071  else if(sval=="blu_BLK70")ival=158;
3072  else if(sval=="blu_BLK72")ival=159;
3073  else if(sval=="blu_BLK73")ival=160;
3074  else if(sval=="blu_BLK75")ival=161;
3075  else if(sval=="blu_BLK77")ival=162;
3076  else if(sval=="blu_BLK78")ival=163;
3077  else if(sval=="blu_BLK79")ival=164;
3078  else if(sval=="blu_BLK81")ival=165;
3079  else if(sval=="conc_BLK")ival=166;
3080  else if(sval=="pvBaffleMother")ival=167;
3081  else if(sval=="pvDPInnerTrackerTube")ival=168;
3082  else if(sval=="pvMHorn1Mother")ival=169;
3083  else if(sval=="pvMHorn2Mother")ival=170;
3084  else if(sval=="pvTargetMother")ival=171;
3085  else if(sval=="stl_slab1")ival=172;
3086  else if(sval=="stl_slab4")ival=173;
3087  else if(sval=="stl_slab5")ival=174;
3088  else if(sval=="stl_slabL")ival=175;
3089  else if(sval=="stl_slabR")ival=176;
3090  return ival;
3091 }
3092 
3094  int ival=0;
3095  if(sval=="AntiLambdaInelastic")ival=1;
3096  else if(sval=="AntiNeutronInelastic")ival=2;
3097  else if(sval=="AntiOmegaMinusInelastic")ival=3;
3098  else if(sval=="AntiProtonInelastic")ival=4;
3099  else if(sval=="AntiSigmaMinusInelastic")ival=5;
3100  else if(sval=="AntiSigmaPlusInelastic")ival=6;
3101  else if(sval=="AntiXiMinusInelastic")ival=7;
3102  else if(sval=="AntiXiZeroInelastic")ival=8;
3103  else if(sval=="Decay")ival=9;
3104  else if(sval=="KaonMinusInelastic")ival=10;
3105  else if(sval=="KaonPlusInelastic")ival=11;
3106  else if(sval=="KaonZeroLInelastic")ival=12;
3107  else if(sval=="KaonZeroSInelastic")ival=13;
3108  else if(sval=="LambdaInelastic")ival=14;
3109  else if(sval=="NeutronInelastic")ival=15;
3110  else if(sval=="OmegaMinusInelastic")ival=16;
3111  else if(sval=="PionMinusInelastic")ival=17;
3112  else if(sval=="PionPlusInelastic")ival=18;
3113  else if(sval=="Primary")ival=19;
3114  else if(sval=="ProtonInelastic")ival=20;
3115  else if(sval=="SigmaMinusInelastic")ival=21;
3116  else if(sval=="SigmaPlusInelastic")ival=22;
3117  else if(sval=="XiMinusInelastic")ival=23;
3118  else if(sval=="XiZeroInelastic")ival=24;
3119  else if(sval=="hElastic")ival=25;
3120  return ival;
3121 }
3122 //END of minerva additions
3123 #endif
static constexpr double cm
Definition: Units.h:68
Double_t muparpx
Definition: g4numi.h:69
Int_t Ndecay
Definition: g4numi.h:52
TString fvol[10]
Definition: g4numi.h:113
Double_t beamx
Definition: flugg.h:77
Double_t ypoint
Definition: flugg.h:60
Float_t Nwtfar
Definition: g3numi.h:34
Double_t Nwtfar
Definition: flugg.h:34
bool LoadParamSet(xmlDocPtr &, std::string cfg)
Definition: GNuMIFlux.cxx:2539
Float_t muparpx
Definition: g3numi.h:53
double fSumWeight
sum of weights for nus thrown so far
Definition: GNuMIFlux.h:408
TVector3 AnglesToAxis(double theta, double phi, std::string units="deg")
Definition: GNuMIFlux.cxx:2876
Float_t pppz
Definition: g3numi.h:46
Float_t tpx
Definition: g3numi.h:65
Definition: g3numi.h:15
Int_t ppmedium
Definition: flugg.h:48
Int_t tgen
Definition: g4numi.h:85
Float_t Ndxdzfar
Definition: g3numi.h:31
Float_t pdpy
Definition: g3numi.h:42
Float_t zpoint
Definition: g3numi.h:61
std::vector< long int > GetIntVector(std::string str)
Definition: GNuMIFlux.cxx:2473
static constexpr double rad
Definition: Units.h:164
Double_t tpy
Definition: g4numi.h:82
Float_t muparpz
Definition: g3numi.h:55
Double_t stoppy[10]
Definition: g4numi.h:106
Float_t tvx
Definition: g3numi.h:62
TLorentzVector fgP4
generated nu 4-momentum beam coord
Definition: GNuMIFlux.h:102
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:45
Double_t Ndydz
Definition: g4numi.h:40
void SetEntryReuse(long int nuse=1)
of times to use entry before moving to next
Definition: GNuMIFlux.cxx:768
Double_t tvz
Definition: g4numi.h:80
double fLengthUnits
units for coord in user exchanges
Definition: GNuMIFlux.h:418
const int kPdgNuE
Definition: PDGCodes.h:28
Double_t pdPx
Definition: g4numi.h:57
Double_t tvx
Definition: g4numi.h:78
Float_t ppvz
Definition: g3numi.h:52
Double_t Vz
Definition: g4numi.h:56
long int fIUse
current # of times an entry has been used
Definition: GNuMIFlux.h:407
Long64_t fFilePOTs
of protons-on-target represented by all files
Definition: GNuMIFlux.h:398
static const unsigned int MAX_N_TRAJ
Maximum number of trajectories to store.
Definition: GNuMIFlux.h:180
#define pERROR
Definition: Messenger.h:59
Double_t ppenergy
Definition: flugg.h:47
Double_t tpz
Definition: g4numi.h:83
Int_t run
current Tree number in a TChain
Definition: g4numi.h:25
Float_t beampz
Definition: g3numi.h:82
Int_t evtno
Definition: flugg.h:22
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
void SetLengthUnits(double user_units)
Set units assumed by user.
Definition: GNuMIFlux.cxx:1206
Int_t pdg[10]
Definition: g4numi.h:93
Double_t Npz
Definition: flugg.h:25
string TrimSpaces(xmlChar *xmls)
Float_t beamy
Definition: g3numi.h:78
Long64_t fIEntry
current flux ntuple entry
Definition: GNuMIFlux.h:396
bool fGenWeighted
does GenerateNext() give weights?
Definition: GNuMIFlux.h:413
Float_t tgppz
Definition: g3numi.h:73
Int_t tgptype
Definition: flugg.h:70
Double_t xpoint
Definition: g4numi.h:75
const int kPdgAntiNuE
Definition: PDGCodes.h:29
const double kRDET
TVector3 fFluxWindowPtUser[3]
user points of flux window
Definition: GNuMIFlux.h:426
#define pFATAL
Definition: Messenger.h:56
Double_t stoppx[10]
Definition: g4numi.h:105
TLorentzVector fgX4User
generated nu 4-position user coord
Definition: GNuMIFlux.h:105
Double_t Nenergyf
Definition: flugg.h:33
Double_t muparpz
Definition: g4numi.h:71
const int kPdgNuMu
Definition: PDGCodes.h:30
bool LoadConfig(string cfg)
load a named configuration
Definition: GNuMIFlux.cxx:2324
Double_t tpy
Definition: flugg.h:66
static constexpr double s
Definition: Units.h:95
Int_t Norig
Definition: g4numi.h:51
Float_t muparpy
Definition: g3numi.h:54
void ParseParamSet(xmlDocPtr &, xmlNodePtr &)
Definition: GNuMIFlux.cxx:2569
Float_t tprivx
Definition: g3numi.h:74
Double_t ppdxdz
Definition: flugg.h:44
Double_t startpy[10]
Definition: g4numi.h:103
TLorentzVector fFluxWindowDir2
extent for flux window (direction 2)
Definition: GNuMIFlux.h:429
Double_t ppvy
Definition: flugg.h:51
Float_t Ndxdz
Definition: g3numi.h:23
Double_t NdxdzFar[2]
Definition: g4numi.h:47
Float_t Necm
Definition: g3numi.h:57
void ParseBeamDir(xmlDocPtr &, xmlNodePtr &)
Definition: GNuMIFlux.cxx:2650
Double_t Ndydzfar
Definition: flugg.h:32
Double_t muparpx
Definition: flugg.h:53
Definition: g4numi.h:18
int fVerbose
how noisy to be when parsing XML
Definition: GNuMIFlux.cxx:94
Double_t beampz
Definition: flugg.h:82
Double_t NdxdzNear[11]
Definition: g4numi.h:43
Int_t run
current Tree number in a TChain
Definition: flugg.h:21
Long64_t fNEntries
number of flux ntuple entries
Definition: GNuMIFlux.h:395
Float_t ppvx
Definition: g3numi.h:50
Float_t Nenergy
Definition: g3numi.h:26
GNuMIFluxXMLHelper(GNuMIFlux *gnumi)
Definition: GNuMIFlux.cxx:75
Double_t startz[10]
Definition: g4numi.h:98
Double_t mupare
Definition: g4numi.h:72
Float_t Vz
Definition: g3numi.h:40
Double_t tgppy
Definition: flugg.h:72
Double_t pdPz
Definition: flugg.h:43
Double_t pdPx
Definition: flugg.h:41
Double_t beamy
Definition: flugg.h:78
Float_t Nwtnear
Definition: g3numi.h:30
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
Double_t pdPz
Definition: g4numi.h:59
Int_t tgen
Definition: g3numi.h:69
TLorentzVector fBeamZero
beam origin in user coords
Definition: GNuMIFlux.h:422
Float_t tvz
Definition: g3numi.h:64
Double_t Vy
Definition: g4numi.h:55
Int_t evtno
Definition: g3numi.h:22
Double_t Ndxdz
Definition: flugg.h:23
Double_t Necm
Definition: g4numi.h:73
Double_t ypoint
Definition: g4numi.h:76
TString ivol[10]
Definition: g4numi.h:112
void MakeCopy(const g3numi *)
pull in from g3 ntuple
Definition: GNuMIFlux.cxx:1476
Float_t beamz
Definition: g3numi.h:79
Float_t Nenergyf
Definition: g3numi.h:33
A list of PDG codes.
Definition: PDGCodeList.h:32
Double_t tgppz
Definition: flugg.h:73
double fMaxWgtFudge
fudge factor for estimating max wgt
Definition: GNuMIFlux.h:402
virtual void SetXMLFileBase(std::string xmlbasename="")
TChain * fNuFluxTree
TTree in g3numi or g4numi // REF ONLY!
Definition: GNuMIFlux.h:389
Float_t Ndydz
Definition: g3numi.h:24
virtual void SetUpstreamZ(double z0)
static constexpr double b
Definition: Units.h:78
ostream & operator<<(ostream &stream, const genie::flux::GJPARCNuFluxPassThroughInfo &info)
Int_t Ntype
Definition: g3numi.h:37
A GENIE flux driver encapsulating the NuMI neutrino flux. It reads-in the official GNUMI neutrino flu...
Definition: GNuMIFlux.h:217
Double_t pppz
Definition: g4numi.h:62
Double_t tvy
Definition: flugg.h:63
string fNuFluxTreeName
Tree name &quot;h10&quot; (g3) or &quot;nudata&quot; (g4)
Definition: GNuMIFlux.h:388
Double_t startpx[10]
Definition: g4numi.h:102
void User2BeamPos(const TLorentzVector &usrxyz, TLorentzVector &beamxyz) const
Definition: GNuMIFlux.cxx:987
int GeantToPdg(int geant_code)
Definition: PDGUtils.cxx:424
TLorentzVector fgX4
generated nu 4-position beam coord
Definition: GNuMIFlux.h:103
Double_t ppenergy
Definition: g4numi.h:63
Int_t Ntype
Definition: g4numi.h:53
Float_t ppvy
Definition: g3numi.h:51
Double_t zpoint
Definition: flugg.h:61
void ParseBeamPos(std::string)
Definition: GNuMIFlux.cxx:2724
Double_t pprodpz[10]
Definition: g4numi.h:110
g4numi * fG4NuMI
g4numi ntuple
Definition: GNuMIFlux.h:392
std::vector< string > fNuFluxFilePatterns
(potentially wildcarded) path(s)
Definition: GNuMIFlux.h:387
PDGCodeList * fPdgCList
list of neutrino pdg-codes to generate
Double_t Ndxdznea
Definition: flugg.h:27
Double_t pprodpx[10]
Definition: g4numi.h:108
const double e
Int_t tgptype
Definition: g3numi.h:70
std::vector< std::string > GetFileList()
list of files currently part of chain
Definition: GNuMIFlux.cxx:2436
Float_t ppenergy
Definition: g3numi.h:47
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
TLorentzVector fFluxWindowDir1
extent for flux window (direction 1)
Definition: GNuMIFlux.h:428
Double_t startpz[10]
Definition: g4numi.h:104
Double_t NdydzNear[11]
Definition: g4numi.h:44
Int_t run
current Tree number in a TChain
Definition: g3numi.h:21
Double_t starty[10]
Definition: g4numi.h:97
Int_t ptype
Definition: flugg.h:49
Int_t tptype
Definition: g3numi.h:68
Double_t tprivx
Definition: flugg.h:74
const double a
Int_t tptype
Definition: g4numi.h:84
void SetBeamRotation(TRotation beamrot)
&lt; beam (0,0,0) relative to user frame, beam direction in user frame
Definition: GNuMIFlux.cxx:913
Double_t tvx
Definition: flugg.h:62
Double_t Nenergy
Definition: flugg.h:26
Double_t tvy
Definition: g4numi.h:79
Float_t tpz
Definition: g3numi.h:67
Int_t Ndecay
Definition: g3numi.h:36
double UnitFromString(string u)
Definition: UnitUtils.cxx:18
Int_t parentId[10]
Definition: g4numi.h:95
void SetMaxWgtScan(double fudge=1.05, long int nentries=2500000)
Definition: GNuMIFlux.h:298
GENIE interface for uniform flux exposure iterface.
string GetXMLFilePath(string basename)
Double_t startx[10]
Definition: g4numi.h:96
Double_t beamz
Definition: flugg.h:79
void ParseRotSeries(xmlDocPtr &, xmlNodePtr &)
Definition: GNuMIFlux.cxx:2778
Int_t Ndecay
Definition: flugg.h:36
Double_t beampy
Definition: flugg.h:81
Double_t stopz[10]
Definition: g4numi.h:101
Double_t Nwtnear
Definition: flugg.h:30
int CalcEnuWgt(const TLorentzVector &xyz, double &enu, double &wgt_xy) const
Definition: GNuMIFlux.cxx:1729
void SetBeamCenter(TVector3 beam0)
Definition: GNuMIFlux.cxx:920
void Initialize(void)
Double_t stopy[10]
Definition: g4numi.h:100
Double_t NenergyF[2]
Definition: g4numi.h:49
Float_t mupare
Definition: g3numi.h:56
Double_t pppz
Definition: flugg.h:46
Double_t ppvy
Definition: g4numi.h:67
Float_t Vy
Definition: g3numi.h:39
Double_t Ndxdz
Definition: g4numi.h:39
Double_t Vy
Definition: flugg.h:39
Definition: flugg.h:15
Double_t ppvz
Definition: flugg.h:52
#define pINFO
Definition: Messenger.h:62
bool fDetLocIsSet
is a flux location (near/far) set?
Definition: GNuMIFlux.h:415
Double_t muparpz
Definition: flugg.h:55
int fNFiles
number of files in chain
Definition: GNuMIFlux.h:394
Float_t Ndxdznea
Definition: g3numi.h:27
long int fNCycles
times to cycle through the ntuple(s)
bool SetFluxWindow(StdFluxWindow_t stdwindow, double padding=0)
return false if unhandled
Definition: GNuMIFlux.cxx:801
const TLorentzVector kPosCenterNearBeam(0., 0., 1039.35, 0.)
Float_t Nimpwt
Definition: g3numi.h:58
flugg * fFlugg
flugg ntuple
Definition: GNuMIFlux.h:393
Float_t Vx
Definition: g3numi.h:38
double fAccumPOTs
POTs used so far.
Definition: GNuMIFlux.h:411
Double_t ppdydz
Definition: g4numi.h:61
Double_t Ndxdzfar
Definition: flugg.h:31
bool LoadConfig(std::string cfg)
Definition: GNuMIFlux.cxx:2496
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
Double_t Necm
Definition: flugg.h:57
virtual std::string GetXMLFileBase() const
Int_t tptype
Definition: flugg.h:68
Double_t stoppz[10]
Definition: g4numi.h:107
#define pWARN
Definition: Messenger.h:60
Double_t Ndydznea
Definition: flugg.h:28
void Print(const Option_t *opt="") const
Definition: GNuMIFlux.cxx:1470
Float_t Npz
Definition: g3numi.h:25
Double_t NWtFar[2]
Definition: g4numi.h:50
Int_t ntrajectory
Definition: g4numi.h:91
Int_t ppmedium
Definition: g3numi.h:48
string TrimSpaces(string input)
Definition: StringUtils.cxx:18
Double_t tgppx
Definition: flugg.h:71
PDGCodeList * fPdgCListRej
list of nu pdg-codes seen but rejected
Float_t ppdxdz
Definition: g3numi.h:44
TVector3 ParseTV3(const std::string &)
Definition: GNuMIFlux.cxx:2895
Float_t tgppy
Definition: g3numi.h:72
Double_t mupare
Definition: flugg.h:56
Float_t tgppx
Definition: g3numi.h:71
Float_t ypoint
Definition: g3numi.h:60
Float_t xpoint
Definition: g3numi.h:59
Double_t Nenergyn
Definition: flugg.h:29
Double_t NdydzFar[2]
Definition: g4numi.h:48
Int_t Norig
Definition: flugg.h:35
double fMaxWeight
max flux neutrino weight in input file
Definition: GNuMIFlux.h:401
Double_t Vx
Definition: flugg.h:38
const TLorentzVector kPosCenterFarBeam(0., 0., 735340.00, 0.)
Double_t tpx
Definition: flugg.h:65
Double_t ppdxdz
Definition: g4numi.h:60
Double_t Ndydz
Definition: flugg.h:24
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
double fLengthScaleB2U
scale factor beam (cm) –&gt; user
Definition: GNuMIFlux.h:419
long int fMaxWgtEntries
of entries in estimating max wgt
Definition: GNuMIFlux.h:403
Float_t tprivy
Definition: g3numi.h:75
void End(void)
Definition: gXSecComp.cxx:210
int fgPdgC
generated nu pdg-code
Definition: GNuMIFlux.h:98
Double_t tprivz
Definition: flugg.h:76
Float_t Nenergyn
Definition: g3numi.h:29
Double_t zpoint
Definition: g4numi.h:77
Double_t ppdydz
Definition: flugg.h:45
enum genie::flux::GNuMIFlux::EStdFluxWindow StdFluxWindow_t
Double_t pdPy
Definition: g4numi.h:58
Double_t ppvx
Definition: flugg.h:50
g3numi * fG3NuMI
g3numi ntuple
Definition: GNuMIFlux.h:391
void ParseWindowSeries(xmlDocPtr &, xmlNodePtr &)
Definition: GNuMIFlux.cxx:2825
Float_t beamx
Definition: g3numi.h:77
ClassImp(CacheBranchFx)
Float_t Ndydznea
Definition: g3numi.h:28
Int_t ptype
Definition: g3numi.h:49
Double_t muparpy
Definition: flugg.h:54
Float_t tvy
Definition: g3numi.h:63
Double_t tpx
Definition: g4numi.h:81
Float_t ppdydz
Definition: g3numi.h:45
Double_t pdPy
Definition: flugg.h:42
TVector3 fWindowNormal
normal direction for flux window
Definition: GNuMIFlux.h:432
Float_t tprivz
Definition: g3numi.h:76
Int_t trackId[10]
Definition: g4numi.h:94
Bool_t overflow
Definition: g4numi.h:92
Int_t evtno
Definition: g4numi.h:26
Double_t NWtNear[11]
Definition: g4numi.h:46
Float_t beampx
Definition: g3numi.h:80
double fLengthScaleU2B
scale factor beam user –&gt; (cm)
Definition: GNuMIFlux.h:420
Double_t tvz
Definition: flugg.h:64
void SetMaxEFudge(double fudge=1.05)
Definition: GNuMIFlux.h:300
int fUseFluxAtDetCenter
use flux at near (-1) or far (+1) det center from ntuple?
Definition: GNuMIFlux.h:416
Double_t xpoint
Definition: flugg.h:59
double fEffPOTsPerNu
what a entry is worth ...
Definition: GNuMIFlux.h:410
Double_t Vx
Definition: g4numi.h:54
void PrintConfig()
print the current configuration
Definition: GNuMIFlux.cxx:2336
Double_t beampx
Definition: flugg.h:80
void SetMaxEnergy(double Ev)
specify maximum flx neutrino energy
Definition: GNuMIFlux.cxx:760
Int_t tgen
Definition: flugg.h:69
#define pNOTICE
Definition: Messenger.h:61
Double_t Nimpwt
Definition: g4numi.h:74
Double_t muparpy
Definition: g4numi.h:70
TLorentzRotation fBeamRotInv
Definition: GNuMIFlux.h:424
TLorentzVector fgP4User
generated nu 4-momentum user coord
Definition: GNuMIFlux.h:104
std::vector< double > GetDoubleVector(std::string str)
Definition: GNuMIFlux.cxx:2450
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:71
Int_t Ntype
Definition: flugg.h:37
string fNuFluxGen
&quot;g3numi&quot; &quot;g4numi&quot; or &quot;flugg&quot;
Definition: GNuMIFlux.h:390
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:57
Double_t tpz
Definition: flugg.h:67
Float_t beampy
Definition: g3numi.h:81
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
string GetAttribute(xmlNodePtr xml_cur, string attr_name)
Int_t ptype
Definition: g4numi.h:65
TLorentzRotation fBeamRot
rotation applied beam –&gt; user coord
Definition: GNuMIFlux.h:423
Double_t Nenergy
Definition: g4numi.h:42
Double_t Vz
Definition: flugg.h:40
long int fNUse
how often to use same entry in a row
Definition: GNuMIFlux.h:406
Double_t ppvz
Definition: g4numi.h:68
Double_t pprodpy[10]
Definition: g4numi.h:109
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
#define FLUXDRIVERREG4(_nsa, _nsb, _name, _fqname)
Int_t Norig
Definition: g3numi.h:35
string Vec3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:80
double fMaxEv
maximum energy
Definition: GNuMIFlux.h:384
Double_t stopx[10]
Definition: g4numi.h:99
Float_t tpy
Definition: g3numi.h:66
Double_t Npz
Definition: g4numi.h:41
TLorentzVector fFluxWindowBase
base point for flux window - beam coord
Definition: GNuMIFlux.h:427
TString proc[10]
Definition: g4numi.h:111
Float_t Ndydzfar
Definition: g3numi.h:32
Float_t pdpz
Definition: g3numi.h:43
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GNuMIFlux.h:409
Double_t ppvx
Definition: g4numi.h:66
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
Double_t tprivy
Definition: flugg.h:75
Double_t ppmedium
Definition: g4numi.h:64
#define pDEBUG
Definition: Messenger.h:63
Double_t NenergyN[11]
Definition: g4numi.h:45
Float_t pdpx
Definition: g3numi.h:41
Double_t Nimpwt
Definition: flugg.h:58