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