genie::flux::GNuMIFlux Class Reference

A GENIE flux driver encapsulating the NuMI neutrino flux. It reads-in the official GNUMI neutrino flux ntuples. Supports both geant3 and geant4 formats. More...

#include <GNuMIFlux.h>

Inheritance diagram for genie::flux::GNuMIFlux:

Collaboration diagram for genie::flux::GNuMIFlux:

Public Types
enum	EStdFluxWindow {   kMinosNearDet , kMinosFarDet , kMinosNearRock , kMinosFarRock ,   kMinosNearCenter , kMinosFarCenter }
typedef enum genie::flux::GNuMIFlux::EStdFluxWindow	StdFluxWindow_t

Public Member Functions
	GNuMIFlux ()
	~GNuMIFlux ()
const PDGCodeList &	FluxParticles (void)
	declare list of flux neutrinos that can be generated (for init. purposes)
double	MaxEnergy (void)
	declare the max flux neutrino energy that can be generated (for init. purposes)
bool	GenerateNext (void)
	generate the next flux neutrino (return false in err)
int	PdgCode (void)
	returns the flux neutrino pdg code
double	Weight (void)
	returns the flux neutrino weight (if any)
const TLorentzVector &	Momentum (void)
	returns the flux neutrino 4-momentum
const TLorentzVector &	Position (void)
	returns the flux neutrino 4-position (note: expect SI rather than physical units)
bool	End (void)
	true if no more flux nu's can be thrown (eg reaching end of beam sim ntuples)
long int	Index (void)
	returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current entry number)
void	Clear (Option_t *opt)
	reset state variables based on opt
void	GenerateWeighted (bool gen_weighted)
	set whether to generate weighted or unweighted neutrinos
const GNuMIFluxPassThroughInfo &	PassThroughInfo (void)
	GNuMIFluxPassThroughInfo.
Long64_t	GetEntryNumber ()
	index in chain
double	GetDecayDist () const
	dist (user units) from dk to current pos
void	MoveToZ0 (double z0)
	move ray origin to user coord Z0
virtual double	GetTotalExposure () const
virtual long int	NFluxNeutrinos () const
double	POT_curr (void)
	current average POT (RWH?)
double	UsedPOTs (void) const
double	SumWeight (void) const
	integrated weight for flux neutrinos looped so far
void	PrintCurrent (void)
	print current entry from leaves
void	PrintConfig ()
	print the current configuration
std::vector< std::string >	GetFileList ()
	list of files currently part of chain
virtual void	LoadBeamSimData (const std::vector< std::string > &filenames, const std::string &det_loc)
virtual void	GetBranchInfo (std::vector< std::string > &branchNames, std::vector< std::string > &branchClassNames, std::vector< void ** > &branchObjPointers)
virtual TTree *	GetMetaDataTree ()
bool	LoadConfig (string cfg)
	load a named configuration
void	SetMaxEnergy (double Ev)
	specify maximum flx neutrino energy
void	SetGenWeighted (bool genwgt=false)
	toggle whether GenerateNext() returns weight=1 flux (initial default false)
void	SetEntryReuse (long int nuse=1)
void	SetTreeName (string name)
	set input tree name (default: "h10")
void	ScanForMaxWeight (void)
	scan for max flux weight (before generating unweighted flux neutrinos)
void	SetMaxWgtScan (double fudge=1.05, long int nentries=2500000)
void	SetMaxEFudge (double fudge=1.05)
void	SetApplyWindowTiltWeight (bool apply=true)
void	SetLengthUnits (double user_units)
	Set units assumed by user.
double	LengthUnits (void) const
	Return user units.
void	SetBeamRotation (TRotation beamrot)
	< beam (0,0,0) relative to user frame, beam direction in user frame
void	SetBeamCenter (TVector3 beam0)
TRotation	GetBeamRotation () const
	rotation to apply from beam->user
TVector3	GetBeamCenter () const
	beam origin in user frame
bool	SetFluxWindow (StdFluxWindow_t stdwindow, double padding=0)
	return false if unhandled
void	SetFluxWindow (TVector3 p1, TVector3 p2, TVector3 p3)
	3 points define a plane (by default in user coordinates)
void	GetFluxWindow (TVector3 &p1, TVector3 &p2, TVector3 &p3) const
	3 points define a plane in beam coordinate
void	UseFluxAtNearDetCenter (void)
	force weights at MINOS detector "center" as found in ntuple
void	UseFluxAtFarDetCenter (void)
void	Beam2UserPos (const TLorentzVector &beamxyz, TLorentzVector &usrxyz) const
void	Beam2UserDir (const TLorentzVector &beamdir, TLorentzVector &usrdir) const
void	Beam2UserP4 (const TLorentzVector &beamp4, TLorentzVector &usrp4) const
void	User2BeamPos (const TLorentzVector &usrxyz, TLorentzVector &beamxyz) const
void	User2BeamDir (const TLorentzVector &usrdir, TLorentzVector &beamdir) const
void	User2BeamP4 (const TLorentzVector &usrp4, TLorentzVector &beamp4) const
TVector3	FluxWindowNormal ()
virtual void	LoadBeamSimData (const std::set< std::string > &filenames, const std::string &det_loc)
virtual void	LoadBeamSimData (const std::string &filename, const std::string &det_loc)
Public Member Functions inherited from genie::GFluxI
virtual	~GFluxI ()
Public Member Functions inherited from genie::flux::GFluxExposureI
	GFluxExposureI (genie::flux::Exposure_t etype)
virtual	~GFluxExposureI ()
const char *	GetExposureUnits () const
	what units are returned by GetTotalExposure?
genie::flux::Exposure_t	GetExposureType () const
Public Member Functions inherited from genie::flux::GFluxFileConfigI
	GFluxFileConfigI ()
virtual	~GFluxFileConfigI ()
virtual void	SetXMLFileBase (std::string xmlbasename="")
virtual std::string	GetXMLFileBase () const
virtual void	SetFluxParticles (const PDGCodeList &particles)
	specify list of flux neutrino species
virtual void	SetUpstreamZ (double z0)
virtual void	SetNumOfCycles (long int ncycle)
	limit cycling through input files

Private Member Functions
bool	GenerateNext_weighted (void)
void	Initialize (void)
void	SetDefaults (void)
void	CleanUp (void)
void	ResetCurrent (void)
void	AddFile (TTree *tree, string fname)
void	CalcEffPOTsPerNu (void)

Private Attributes
double	fMaxEv
	maximum energy
bool	fEnd
	end condition reached
std::vector< string >	fNuFluxFilePatterns
	(potentially wildcarded) path(s)
string	fNuFluxTreeName
	Tree name "h10" (g3) or "nudata" (g4)
TChain *	fNuFluxTree
	TTree in g3numi or g4numi // REF ONLY!
string	fNuFluxGen
	"g3numi" "g4numi" or "flugg"
g3numi *	fG3NuMI
	g3numi ntuple
g4numi *	fG4NuMI
	g4numi ntuple
flugg *	fFlugg
	flugg ntuple
int	fNFiles
	number of files in chain
Long64_t	fNEntries
	number of flux ntuple entries
Long64_t	fIEntry
	current flux ntuple entry
Long64_t	fNuTot
	cummulative # of entries (=fNEntries)
Long64_t	fFilePOTs
double	fWeight
	current neutrino weight, =1 if generating unweighted entries
double	fMaxWeight
	max flux neutrino weight in input file
double	fMaxWgtFudge
	fudge factor for estimating max wgt
long int	fMaxWgtEntries
double	fMaxEFudge
	fudge factor for estmating max enu (0=> use fixed 120GeV)
long int	fNUse
	how often to use same entry in a row
long int	fIUse
	current # of times an entry has been used
double	fSumWeight
	sum of weights for nus thrown so far
long int	fNNeutrinos
	number of flux neutrinos thrown so far
double	fEffPOTsPerNu
	what a entry is worth ...
double	fAccumPOTs
	POTs used so far.
bool	fGenWeighted
	does GenerateNext() give weights?
bool	fApplyTiltWeight
	wgt due to window normal not || beam
bool	fDetLocIsSet
	is a flux location (near/far) set?
int	fUseFluxAtDetCenter
	use flux at near (-1) or far (+1) det center from ntuple?
double	fLengthUnits
	units for coord in user exchanges
double	fLengthScaleB2U
	scale factor beam (cm) --> user
double	fLengthScaleU2B
	scale factor beam user --> (cm)
TLorentzVector	fBeamZero
	beam origin in user coords
TLorentzRotation	fBeamRot
	rotation applied beam --> user coord
TLorentzRotation	fBeamRotInv
TVector3	fFluxWindowPtUser [3]
	user points of flux window
TLorentzVector	fFluxWindowBase
	base point for flux window - beam coord
TLorentzVector	fFluxWindowDir1
	extent for flux window (direction 1)
TLorentzVector	fFluxWindowDir2
	extent for flux window (direction 2)
double	fFluxWindowLen1
double	fFluxWindowLen2
TVector3	fWindowNormal
	normal direction for flux window
TLorentzVector	fgX4dkvtx
	decay 4-position beam coord
GNuMIFluxPassThroughInfo *	fCurEntry
	copy of current ntuple entry info (owned structure)

Additional Inherited Members
Static Public Member Functions inherited from genie::flux::GFluxExposureI
static const char *	AsString (genie::flux::Exposure_t etype)
static genie::flux::Exposure_t	StringToEnum (const char *chars, int maxChar=0)
Protected Member Functions inherited from genie::GFluxI
	GFluxI ()
Protected Attributes inherited from genie::flux::GFluxFileConfigI
PDGCodeList *	fPdgCList
	list of neutrino pdg-codes to generate
PDGCodeList *	fPdgCListRej
	list of nu pdg-codes seen but rejected
std::string	fXMLbasename
	XML file that might hold config param_sets.
long int	fNCycles
long int	fICycle
double	fZ0

Detailed Description

A GENIE flux driver encapsulating the NuMI neutrino flux. It reads-in the official GNUMI neutrino flux ntuples. Supports both geant3 and geant4 formats.

GNuMIFlux:

An implementation of the GFluxI interface that provides NuMI flux

References:\n http://www.hep.utexas.edu/~zarko/wwwgnumi/v19/

Author

Robert Hatcher <rhatcher \at fnal.gov> Fermi National Accelerator Laboratory

Created:\n Jun 27, 2008

License:\n Copyright (c) 2003-2025, The GENIE Collaboration

For the full text of the license visit http://copyright.genie-mc.org

Definition at line 217 of file GNuMIFlux.h.

Member Typedef Documentation

StdFluxWindow_t

typedef enum genie::flux::GNuMIFlux::EStdFluxWindow genie::flux::GNuMIFlux::StdFluxWindow_t

Member Enumeration Documentation

EStdFluxWindow

enum genie::flux::GNuMIFlux::EStdFluxWindow

Enumerator
kMinosNearDet
kMinosFarDet
kMinosNearRock
kMinosFarRock
kMinosNearCenter
kMinosFarCenter

Definition at line 330 of file GNuMIFlux.h.

                              {
    kMinosNearDet,      // appropriate for Near Detector
    kMinosFarDet,       // appropriate for Far Detector
    kMinosNearRock,     // appropriate for Near rock generation
    kMinosFarRock,      // appropriate for Far rock generation
    kMinosNearCenter,   // point location mimic near value in ntuple
    kMinosFarCenter     // point location mimic far value in ntuple
  } StdFluxWindow_t;

Constructor & Destructor Documentation

GNuMIFlux()

GNuMIFlux::GNuMIFlux

(

)

Definition at line 113 of file GNuMIFlux.cxx.

                     :
  GFluxExposureI(genie::flux::kPOTs)
{
  this->Initialize();
}

References genie::flux::GFluxExposureI::GFluxExposureI(), Initialize(), and genie::flux::kPOTs.

~GNuMIFlux()

GNuMIFlux::~GNuMIFlux

(

)

Definition at line 119 of file GNuMIFlux.cxx.

{
  this->CleanUp();
}

References CleanUp().

Member Function Documentation

AddFile()

void GNuMIFlux::AddFile ( TTree * tree, string fname )

private

Definition at line 1133 of file GNuMIFlux.cxx.

{
  // Add a file to the chain
 
  ULong64_t nentries = thetree->GetEntries();
 
  // first/last "evtno" are the proton # of the first/last proton
  // that generated a neutrino ... not necessarily true first/last #
  // estimate we're probably not off by more than 100 ...
  // !Important change
  // Some files (due to some intermediate flugg issues) list evtno==0
  // when it isn't really true, we need to scan nearby values in case the
  // last entry is one of these (otherwise file contributes 0 POTs).
  // Also assume quantization of 500 (instead of 100).
 
  thetree->SetMakeClass(1); // need full ntuple decomposition for
  // the SetBranchAddress to work on g4numi ntuples.  Works fine
  // without it on gnumi (geant3) and flugg ntuples [each with their
  // own slight differences] but shouldn't harm them either.
 
  Int_t evtno = 0;
  TBranch* br_evtno = 0;
  thetree->SetBranchAddress("evtno",&evtno, &br_evtno);
  Int_t evt_1 = 0x7FFFFFFF;
  Int_t evt_N = 1;
#define OLDGUESS
#ifdef OLDGUESS
  for (int j=0; j<50; ++j) {
    thetree->GetEntry(j);
    if (evtno != 0) evt_1 = TMath::Min(evtno,evt_1);
    thetree->GetEntry(nentries-1 -j );
    if (evtno != 0) evt_N = TMath::Max(evtno,evt_N);
  }
  // looks like low counts are due to "evtno" not including
  // protons that miss the actual target (hit baffle, etc)
  // this will vary with beam conditions parameters
  // so we should round way up, those generating flugg files
  // aren't going to quantize less than 1000
  // though 500 would probably be okay, 100 would not.
  const Int_t    nquant = 1000; // 500;  // 100
  const Double_t rquant = nquant;
#else
  for (int j=0; j<50; ++j) {
    thetree->GetEntry(j);
    if (evtno != 0) evt_1 = TMath::Min(evtno,evt_1);
    std::cout << "[" << j << "] evtno=" << evtno << " evt_1=" << evt_1 << std::endl;
  }
  for (int j=0; j<50; ++j) {
    thetree->GetEntry(nentries-1 -j );
    if (evtno != 0) evt_N = TMath::Max(evtno,evt_N);
    std::cout << "[" << (nentries-1-j) << "] evtno=" << evtno << " evt_N=" << evt_N << std::endl;
  }
 
  Int_t    nquant = 1000; // 500;
  Double_t rquant = nquant;
#endif
 
  Int_t est_1 = (TMath::FloorNint(evt_1/rquant))*nquant + 1;
  Int_t est_N = (TMath::FloorNint((evt_N-1)/rquant)+1)*nquant;
  ULong64_t npots = est_N - est_1 + 1;
  LOG("Flux",pNOTICE) //INFO)
    << fNuFluxTreeName << "->AddFile() of " << nentries << " entries ["
    << evt_1 << ":" << evt_N << "%" << nquant
    << "(" <<  est_1 << ":" << est_N << ")="
    << npots <<" POTs] in {" << fNuFluxGen << "} file: " << fname;
  fNuTot    += nentries;
  fFilePOTs += npots;
  fNFiles++;
 
  fNuFluxTree->AddFile(fname.c_str());
}

References fFilePOTs, fNFiles, fNuFluxGen, fNuFluxTree, fNuFluxTreeName, fNuTot, LOG, and pNOTICE.

Referenced by LoadBeamSimData().

Beam2UserDir()

void GNuMIFlux::Beam2UserDir	(	const TLorentzVector &	beamdir,
		TLorentzVector &	usrdir ) const

Definition at line 976 of file GNuMIFlux.cxx.

{
  usrdir = fLengthScaleB2U*(fBeamRot*beamdir);
}

References fBeamRot, and fLengthScaleB2U.

Beam2UserP4()

void GNuMIFlux::Beam2UserP4	(	const TLorentzVector &	beamp4,
		TLorentzVector &	usrp4 ) const

Definition at line 981 of file GNuMIFlux.cxx.

{
  usrp4 = fBeamRot*beamp4;
}

References fBeamRot.

Referenced by GenerateNext_weighted().

Beam2UserPos()

void GNuMIFlux::Beam2UserPos	(	const TLorentzVector &	beamxyz,
		TLorentzVector &	usrxyz ) const

Definition at line 971 of file GNuMIFlux.cxx.

{
  usrxyz = fLengthScaleB2U*(fBeamRot*beamxyz) + fBeamZero;
}

References fBeamRot, fBeamZero, and fLengthScaleB2U.

Referenced by GenerateNext_weighted(), and SetFluxWindow().

CalcEffPOTsPerNu()

void GNuMIFlux::CalcEffPOTsPerNu ( void )

private

Definition at line 439 of file GNuMIFlux.cxx.

{
  // do this if flux window changes or # of files changes
 
  if (!fNuFluxTree) return;  // not yet fully configured
 
  // effpots = mc_pots * (wgtfunction-area) / window-area / wgt-max-est
  //   wgtfunction-area = pi * radius-det-element^2 = pi * (100.cm)^2
 
  // this should match what is used in the CalcEnuWgt()
  const double kRDET = 100.0;   // set to flux per 100 cm radius
  const double kRDET2 = kRDET * kRDET;
  double flux_area = fFluxWindowDir1.Vect().Cross(fFluxWindowDir2.Vect()).Mag();
  LOG("Flux",pNOTICE) << "in CalcEffPOTsPerNu, area = " << flux_area;
 
  if ( flux_area < 1.0e-30 ) {
    LOG("Flux", pWARN)
          << "CalcEffPOTsPerNu called with flux window area effectively zero";
    flux_area = 1;
  }
  double area_ratio = TMath::Pi() * kRDET2 / flux_area;
  fEffPOTsPerNu = area_ratio * ( (double)fFilePOTs / (double)fNEntries );
}

References fEffPOTsPerNu, fFilePOTs, fFluxWindowDir1, fFluxWindowDir2, fNEntries, fNuFluxTree, kRDET, LOG, pNOTICE, and pWARN.

Referenced by LoadBeamSimData(), SetFluxWindow(), and SetFluxWindow().

CleanUp()

void GNuMIFlux::CleanUp ( void )

private

Definition at line 1115 of file GNuMIFlux.cxx.

{
  LOG("Flux", pNOTICE) << "Cleaning up...";
 
  if (fPdgCList)    delete fPdgCList;
  if (fPdgCListRej) delete fPdgCListRej;
  if (fCurEntry)    delete fCurEntry;
 
  if ( fG3NuMI )    delete fG3NuMI;
  if ( fG4NuMI )    delete fG4NuMI;
  if ( fFlugg  )    delete fFlugg;
 
  LOG("Flux", pNOTICE)
    << " flux file cycles: " << fICycle << " of " << fNCycles
    << ", entry " << fIEntry << " use: " << fIUse << " of " << fNUse;
}

References fCurEntry, fFlugg, fG3NuMI, fG4NuMI, genie::flux::GFluxFileConfigI::fICycle, fIEntry, fIUse, genie::flux::GFluxFileConfigI::fNCycles, fNUse, genie::flux::GFluxFileConfigI::fPdgCList, genie::flux::GFluxFileConfigI::fPdgCListRej, LOG, and pNOTICE.

Referenced by ~GNuMIFlux().

Clear()

void GNuMIFlux::Clear ( Option_t * opt )

virtual

reset state variables based on opt

Implements genie::GFluxI.

Definition at line 1009 of file GNuMIFlux.cxx.

{
  // Clear the driver state
  //
  LOG("Flux", pWARN) << "GSimpleNtpFlux::Clear(" << opt << ") called";
  // do it in all cases, but EVGDriver/GMCJDriver will pass "CycleHistory"
 
  fICycle     = 0;
 
  fSumWeight  = 0;
  fNNeutrinos = 0;
  fAccumPOTs  = 0;
}

References fAccumPOTs, genie::flux::GFluxFileConfigI::fICycle, fNNeutrinos, fSumWeight, LOG, and pWARN.

End()

bool genie::flux::GNuMIFlux::End ( void )

inlinevirtual

true if no more flux nu's can be thrown (eg reaching end of beam sim ntuples)

Implements genie::GFluxI.

Definition at line 237 of file GNuMIFlux.h.

{ return  fEnd;                 }

References fEnd.

Referenced by GenerateNext().

FluxParticles()

const PDGCodeList & genie::flux::GNuMIFlux::FluxParticles ( void )

inlinevirtual

declare list of flux neutrinos that can be generated (for init. purposes)

Implements genie::GFluxI.

Definition at line 230 of file GNuMIFlux.h.

{ return *fPdgCList;            }

References genie::flux::GFluxFileConfigI::fPdgCList.

FluxWindowNormal()

TVector3 genie::flux::GNuMIFlux::FluxWindowNormal ( )

inline

Definition at line 368 of file GNuMIFlux.h.

{ return fWindowNormal; }

References fWindowNormal.

Referenced by GenerateNext_weighted().

GenerateNext()

bool GNuMIFlux::GenerateNext ( void )

virtual

generate the next flux neutrino (return false in err)

Implements genie::GFluxI.

Definition at line 139 of file GNuMIFlux.cxx.

{
// Get next (unweighted) flux ntuple entry on the specified detector location
//
  RandomGen* rnd = RandomGen::Instance();
  while ( true ) {
     // Check for end of flux ntuple
     bool end = this->End();
     if ( end ) {
       LOG("Flux", pNOTICE) << "GenerateNext signaled End() ";
       return false;
     }
 
     // Get next weighted flux ntuple entry
     bool nextok = this->GenerateNext_weighted();
     if ( fGenWeighted ) return nextok;
     if ( ! nextok ) continue;
 
     /* RWH - debug purposes
     if ( fNCycles == 0 ) {
       LOG("Flux", pNOTICE)
          << "Got flux entry: " << fIEntry
          << " - Cycle: "<< fICycle << "/ infinite";
     } else {
       LOG("Flux", pNOTICE)
          << "Got flux entry: "<< fIEntry
          << " - Cycle: "<< fICycle << "/"<< fNCycles;
     }
     */
 
     // Get fractional weight & decide whether to accept curr flux neutrino
     double f = this->Weight() / fMaxWeight;
     //LOG("Flux", pNOTICE)
     //   << "Curr flux neutrino fractional weight = " << f;
     if (f > 1.) {
       fMaxWeight = this->Weight() * fMaxWgtFudge; // bump the weight
       LOG("Flux", pERROR)
         << "** Fractional weight = " << f
         << " > 1 !! Bump fMaxWeight estimate to " << fMaxWeight
         << PassThroughInfo();
     }
     double r = (f < 1.) ? rnd->RndFlux().Rndm() : 0;
     bool accept = ( r < f );
     if ( accept ) {
 
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
       LOG("Flux", pNOTICE)
         << "Generated beam neutrino: "
         << "\n pdg-code: " << fCurEntry->fgPdgC
         << "\n p4: " << utils::print::P4AsShortString(&(fCurEntry->fgP4))
         << "\n x4: " << utils::print::X4AsString(&(fCurEntry->fgX4))
         << "\n p4: " << utils::print::P4AsShortString(&(fCurEntry->fgP4User))
         << "\n x4: " << utils::print::X4AsString(&(fCurEntry->fgX4User));
#endif
 
       fWeight = 1.;
       return true;
     }
 
     //LOG("Flux", pNOTICE)
     //  << "** Rejecting current flux neutrino based on the flux weight only";
  }
  return false;
}

References End(), fCurEntry, fGenWeighted, fMaxWeight, fMaxWgtFudge, fWeight, GenerateNext_weighted(), genie::RandomGen::Instance(), LOG, genie::utils::print::P4AsShortString(), PassThroughInfo(), pERROR, pNOTICE, genie::RandomGen::RndFlux(), Weight(), and genie::utils::print::X4AsString().

GenerateNext_weighted()

bool GNuMIFlux::GenerateNext_weighted ( void )

private

user might modify list via SetFluxParticles() in order to reject certain flavors, even if they're found in the file. So don't make a big fuss. Spit out a single message and then stop reporting that flavor as problematic.

Definition at line 204 of file GNuMIFlux.cxx.

{
// Get next (weighted) flux ntuple entry on the specified detector location
//
 
  // Check whether a flux ntuple has been loaded
  if ( ! fG3NuMI && ! fG4NuMI && ! fFlugg ) {
     LOG("Flux", pFATAL)
          << "The flux driver has not been properly configured";
     //return false; //  don't do this - creates an infinite loop!
     exit(1);
  }
 
  // Reuse an entry?
  //std::cout << " ***** iuse " << fIUse << " nuse " << fNUse
  //          << " ientry " << fIEntry << " nentry " << fNEntries
  //          << " icycle " << fICycle << " ncycle " << fNCycles << std::endl;
  if ( fIUse < fNUse && fIEntry >= 0 ) {
    // Reuse this entry
    fIUse++;
  } else {
    // Reset previously generated neutrino code / 4-p / 4-x
    this->ResetCurrent();
    // Move on, read next flux ntuple entry
    fIEntry++;
    if ( fIEntry >= fNEntries ) {
      // Ran out of entries @ the current cycle of this flux file
      // Check whether more (or infinite) number of cycles is requested
      if ( fICycle < fNCycles || fNCycles == 0 ) {
        fICycle++;
        fIEntry=0;
      } else {
        LOG("Flux", pWARN)
          << "No more entries in input flux neutrino ntuple, cycle "
          << fICycle << " of " << fNCycles;
        fEnd = true;
        // assert(0);
        return false;
      }
    }
 
    if ( fG3NuMI ) {
      fG3NuMI->GetEntry(fIEntry);
      fCurEntry->MakeCopy(fG3NuMI);
    } else if ( fG4NuMI ) {
      fG4NuMI->GetEntry(fIEntry);
      fCurEntry->MakeCopy(fG4NuMI);
    } else if ( fFlugg ) {
      fFlugg->GetEntry(fIEntry);
      fCurEntry->MakeCopy(fFlugg);
    } else {
      LOG("Flux", pERROR) << "No ntuple configured";
      fEnd = true;
      //assert(0);
      return false;
    }
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("Flux",pDEBUG)
    << "got " << fNNeutrinos << " new fIEntry " << fIEntry
    << " evtno " << fCurEntry->evtno;
#endif
 
    fIUse = 1;
    fCurEntry->pcodes = 0;  // fetched entry has geant codes
    fCurEntry->units  = 0;  // fetched entry has original units
 
    // Convert the current gnumi neutrino flavor mode into a neutrino pdg code
    // Also convert other particle codes in GNuMIFluxPassThroughInfo to PDG
    fCurEntry->ConvertPartCodes();
    // here we might want to do flavor oscillations or simple mappings
    fCurEntry->fgPdgC = fCurEntry->ntype;
  }
 
  // Check neutrino pdg against declared list of neutrino species declared
  // by the current instance of the NuMI neutrino flux driver.
  // No undeclared neutrino species will be accepted at this point as GENIE
  // has already been configured to handle the specified list.
  // Make sure that the appropriate list of flux neutrino species was set at
  // initialization via GNuMIFlux::SetFluxParticles(const PDGCodeList &)
 
  // update the # POTs, number of neutrinos
  // do this HERE (before rejecting flavors that users might be weeding out)
  // in order to keep the POT accounting correct.  This allows one to get
  // the right normalization for generating only events from the intrinsic
  // nu_e entries.
  fAccumPOTs += fEffPOTsPerNu / fMaxWeight;
  fNNeutrinos++;
 
  if ( ! fPdgCList->ExistsInPDGCodeList(fCurEntry->fgPdgC) ) {
     /// user might modify list via SetFluxParticles() in order to reject certain
     /// flavors, even if they're found in the file.  So don't make a big fuss.
     /// Spit out a single message and then stop reporting that flavor as problematic.
     int badpdg = fCurEntry->fgPdgC;
     if ( ! fPdgCListRej->ExistsInPDGCodeList(badpdg) ) {
       fPdgCListRej->push_back(badpdg);
       LOG("Flux", pWARN)
         << "Encountered neutrino specie (" << badpdg
         << " pcodes=" << fCurEntry->pcodes << ")"
         << " that wasn't in SetFluxParticles() list, "
         << "\nDeclared list of neutrino species: " << *fPdgCList;
     }
     return false;
  }
 
  // Update the curr neutrino weight and energy
 
  // Check current neutrino energy against the maximum flux neutrino energy
  // declared by the current instance of the NuMI neutrino flux driver.
  // No flux neutrino exceeding that maximum energy will be accepted at this
  // point as that maximum energy has already been used for normalizing the
  // interaction probabilities.
  // Make sure that the appropriate maximum flux neutrino energy was set at
  // initialization via GNuMIFlux::SetMaxEnergy(double Ev)
 
  fCurEntry->fgX4 = fFluxWindowBase;
 
  double Ev = 0;
  double& wgt_xy = fCurEntry->fgXYWgt;
  switch ( fUseFluxAtDetCenter ) {
  case -1:  // near detector
    wgt_xy   = fCurEntry->nwtnear;
    Ev       = fCurEntry->nenergyn;
    break;
  case +1:  // far detector
    wgt_xy   = fCurEntry->nwtfar;
    Ev       = fCurEntry->nenergyf;
    break;
  default:  // recalculate on x-y window
    RandomGen * rnd = RandomGen::Instance();
    fCurEntry->fgX4 += ( rnd->RndFlux().Rndm()*fFluxWindowDir1 +
                         rnd->RndFlux().Rndm()*fFluxWindowDir2   );
    fCurEntry->CalcEnuWgt(fCurEntry->fgX4,Ev,wgt_xy);
    break;
  }
 
  if (Ev > fMaxEv) {
     LOG("Flux", pWARN)
          << "Flux neutrino energy exceeds declared maximum neutrino energy"
          << "\nEv = " << Ev << "(> Ev{max} = " << fMaxEv << ")";
  }
 
  // Set the current flux neutrino 4-momentum
  // this is in *beam* coordinates
  fgX4dkvtx = TLorentzVector( fCurEntry->vx,
                              fCurEntry->vy,
                              fCurEntry->vz, 0.);
  // don't use TLorentzVector here for Mag() due to - on metric
  // normalize direction to 1.0
  TVector3 dirNu = (fCurEntry->fgX4.Vect() - fgX4dkvtx.Vect()).Unit();
  fCurEntry->fgP4.SetPxPyPzE( Ev*dirNu.X(),
                              Ev*dirNu.Y(),
                              Ev*dirNu.Z(), Ev);
 
  // calculate the weight, potentially includes effect from tilted window
  // must be done *after* neutrino direction is determined
  fWeight = fCurEntry->nimpwt * fCurEntry->fgXYWgt;  // full weight
  if ( fApplyTiltWeight ) {
    double tiltwgt = dirNu.Dot( FluxWindowNormal() );
    fWeight *= TMath::Abs( tiltwgt );
  }
 
  // update sume of weights
  fSumWeight += this->Weight();
 
  // Set the current flux neutrino 4-position, direction in user coord
  Beam2UserP4(fCurEntry->fgP4,fCurEntry->fgP4User);
  Beam2UserPos(fCurEntry->fgX4,fCurEntry->fgX4User);
 
  // if desired, move to user specified user coord z
  if ( TMath::Abs(fZ0) < 1.0e30 ) this->MoveToZ0(fZ0);
 
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("Flux", pINFO)
    << "Generated neutrino: " << fIEntry << " " << fCurEntry->evtno
    << " nenergyn " << fCurEntry->nenergyn
    << "\n pdg-code: " << fCurEntry->fgPdgC
    << "\n p4 beam: " << utils::print::P4AsShortString(&fCurEntry->fgP4)
    << "\n x4 beam: " << utils::print::X4AsString(&fCurEntry->fgX4)
    << "\n p4 user: " << utils::print::P4AsShortString(&(fCurEntry->fgP4User))
    << "\n x4 user: " << utils::print::X4AsString(&(fCurEntry->fgX4User));
#endif
  if ( Ev > fMaxEv ) {
    LOG("Flux", pFATAL)
      << "Generated neutrino had E_nu = " << Ev << " > " << fMaxEv
      << " maximum ";
    assert(0);
  }
 
  return true;
}

References Beam2UserP4(), Beam2UserPos(), fAccumPOTs, fApplyTiltWeight, fCurEntry, fEffPOTsPerNu, fEnd, fFlugg, fFluxWindowBase, fFluxWindowDir1, fFluxWindowDir2, fG3NuMI, fG4NuMI, fgX4dkvtx, genie::flux::GFluxFileConfigI::fICycle, fIEntry, fIUse, FluxWindowNormal(), fMaxEv, fMaxWeight, genie::flux::GFluxFileConfigI::fNCycles, fNEntries, fNNeutrinos, genie::flux::GFluxFileConfigI::fPdgCList, genie::flux::GFluxFileConfigI::fPdgCListRej, fSumWeight, fUseFluxAtDetCenter, fWeight, genie::flux::GFluxFileConfigI::fZ0, genie::RandomGen::Instance(), LOG, MoveToZ0(), genie::utils::print::P4AsShortString(), pDEBUG, pERROR, pFATAL, pINFO, pWARN, ResetCurrent(), genie::RandomGen::RndFlux(), Weight(), and genie::utils::print::X4AsString().

Referenced by GenerateNext(), and ScanForMaxWeight().

GenerateWeighted()

void GNuMIFlux::GenerateWeighted ( bool gen_weighted )

virtual

set whether to generate weighted or unweighted neutrinos

Implements genie::GFluxI.

Definition at line 1023 of file GNuMIFlux.cxx.

{
  // Set whether to generate weighted rays
  //
  fGenWeighted = gen_weighted;
}

References fGenWeighted.

GetBeamCenter()

TVector3 GNuMIFlux::GetBeamCenter

(

)

const

beam origin in user frame

Definition at line 941 of file GNuMIFlux.cxx.

{
  TVector3 beam0 = fBeamZero.Vect();
  return beam0;
}

References fBeamZero.

GetBeamRotation()

TRotation GNuMIFlux::GetBeamRotation

(

)

const

rotation to apply from beam->user

Definition at line 928 of file GNuMIFlux.cxx.

{
  // rotation is really only 3-d vector, but we'll be operating on LorentzV's
  // give people back the original TRotation ... not pretty
  // ... it think this is right
  TRotation rot3;
  const TLorentzRotation& rot4 = fBeamRot;
  TVector3 newX(rot4.XX(),rot4.XY(),rot4.XZ());
  TVector3 newY(rot4.YX(),rot4.YY(),rot4.YZ());
  TVector3 newZ(rot4.ZX(),rot4.ZY(),rot4.ZZ());
  rot3.RotateAxes(newX,newY,newZ);
  return rot3.Inverse();
}

References fBeamRot.

GetBranchInfo()

void GNuMIFlux::GetBranchInfo ( std::vector< std::string > & branchNames, std::vector< std::string > & branchClassNames, std::vector< void ** > & branchObjPointers )

virtual

allow caller to copy current status / ntuple entry info in the output file by providing copies of internal info

Assumes that branch object pointers will not change which may require either a copy be made or, if using the class directly for reading the branch, one must force ROOT to not autodelete: myns::MyClassType* fCurrMyClass = new myns::MyClassType; myTree->SetBranchAddress("bname",&fCurMyClass); //? TBranch* b = myTree->GetBranch("bname"); //? b->SetAutoDelete(false);

ensure vectors are sized sufficiently (or use .push_back()) branchNames[i] = "bname" branchClassNames[i] = "myns::MyClassType" branchObjPointers[i] = (void**)&fCurMyClass;

Reimplemented from genie::flux::GFluxFileConfigI.

Definition at line 663 of file GNuMIFlux.cxx.

{
  // allow flux driver to report back current status and/or ntuple entry
  // info for possible recording in the output file by supplying
  // the class name, and a pointer to the object that will be filled
  // as well as a suggested name for the branch.
 
  branchNames.push_back("flux");
  branchClassNames.push_back("genie::flux::GNuMIFluxPassThroughInfo");
  branchObjPointers.push_back((void**)&fCurEntry);
 
}

References fCurEntry.

GetDecayDist()

double GNuMIFlux::GetDecayDist

(

)

const

dist (user units) from dk to current pos

Definition at line 395 of file GNuMIFlux.cxx.

{
  // return distance (user units) between dk point and start position
  // these are in beam units
  TVector3 x3diff = fCurEntry->fgX4.Vect() - fgX4dkvtx.Vect();
  return x3diff.Mag() * fLengthScaleB2U;
}

References fCurEntry, fgX4dkvtx, and fLengthScaleB2U.

GetEntryNumber()

Long64_t genie::flux::GNuMIFlux::GetEntryNumber ( )

inline

index in chain

Definition at line 253 of file GNuMIFlux.h.

References fIEntry.

GetFileList()

std::vector< std::string > GNuMIFlux::GetFileList

(

)

list of files currently part of chain

Definition at line 2436 of file GNuMIFlux.cxx.

{
  std::vector<std::string> flist;
  TObjArray *fileElements=fNuFluxTree->GetListOfFiles();
  TIter next(fileElements);
  TChainElement *chEl=0;
  while (( chEl=(TChainElement*)next() )) {
    flist.push_back(chEl->GetTitle());
  }
  return flist;
}

References fNuFluxTree.

Referenced by PrintConfig().

GetFluxWindow()

void GNuMIFlux::GetFluxWindow	(	TVector3 &	p1,
		TVector3 &	p2,
		TVector3 &	p3 ) const

3 points define a plane in beam coordinate

Definition at line 904 of file GNuMIFlux.cxx.

{
  // return flux window points
  p0 = fFluxWindowPtUser[0];
  p1 = fFluxWindowPtUser[1];
  p2 = fFluxWindowPtUser[2];
 
}

References fFluxWindowPtUser.

GetMetaDataTree()

TTree * GNuMIFlux::GetMetaDataTree ( )

virtual

Reimplemented from genie::flux::GFluxFileConfigI.

Definition at line 677 of file GNuMIFlux.cxx.

{ return 0; } // there is none

GetTotalExposure()

double GNuMIFlux::GetTotalExposure ( ) const

virtual

Implements genie::flux::GFluxExposureI.

Definition at line 125 of file GNuMIFlux.cxx.

{
  // complete the GFluxExposureI interface
  return UsedPOTs();
}

References UsedPOTs().

Index()

long int genie::flux::GNuMIFlux::Index ( void )

inlinevirtual

returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current entry number)

Implements genie::GFluxI.

Definition at line 238 of file GNuMIFlux.h.

{ return  fIEntry;              }

References fIEntry.

Initialize()

void GNuMIFlux::Initialize ( void )

private

Definition at line 1030 of file GNuMIFlux.cxx.

{
  LOG("Flux", pNOTICE) << "Initializing GNuMIFlux driver";
 
  fMaxEv           =  0;
  fEnd             =  false;
 
  fCurEntry        = new GNuMIFluxPassThroughInfo;
 
  fNuFluxTree      =  0;
  fG3NuMI          =  0;
  fG4NuMI          =  0;
  fFlugg           =  0;
  fNuFluxTreeName  = "";
  fNuFluxGen       = "";
  fNFiles          =  0;
 
  fNEntries        =  0;
  fIEntry          = -1;
  fICycle          =  0;
  fNUse            =  1;
  fIUse            =  999999;
 
  fNuTot           = 0;
  fFilePOTs        = 0;
 
  fMaxWeight       = -1;
  fMaxWgtFudge     =  1.05;
  fMaxWgtEntries   = 2500000;
  fMaxEFudge       =  0;
 
  fSumWeight       =  0;
  fNNeutrinos      =  0;
  fEffPOTsPerNu    =  0;
  fAccumPOTs       =  0;
 
  fGenWeighted     = false;
  fApplyTiltWeight = true;
  fUseFluxAtDetCenter = 0;
  fDetLocIsSet        = false;
  // by default assume user length is cm
  SetLengthUnits(genie::utils::units::UnitFromString("m"));
 
  this->SetDefaults();
  this->ResetCurrent();
}

References fAccumPOTs, fApplyTiltWeight, fCurEntry, fDetLocIsSet, fEffPOTsPerNu, fEnd, fFilePOTs, fFlugg, fG3NuMI, fG4NuMI, fGenWeighted, genie::flux::GFluxFileConfigI::fICycle, fIEntry, fIUse, fMaxEFudge, fMaxEv, fMaxWeight, fMaxWgtEntries, fMaxWgtFudge, fNEntries, fNFiles, fNNeutrinos, fNuFluxGen, fNuFluxTree, fNuFluxTreeName, fNUse, fNuTot, fSumWeight, fUseFluxAtDetCenter, LOG, pNOTICE, ResetCurrent(), SetDefaults(), SetLengthUnits(), and genie::utils::units::UnitFromString().

Referenced by GNuMIFlux().

LengthUnits()

double GNuMIFlux::LengthUnits

(

void

)

const

Return user units.

Definition at line 1240 of file GNuMIFlux.cxx.

{
  // Return the scale factor for lengths the user is getting
  double cm = genie::utils::units::UnitFromString("cm");
  return fLengthScaleB2U * cm ;
}

References fLengthScaleB2U, and genie::utils::units::UnitFromString().

LoadBeamSimData() [1/3]

void genie::flux::GFluxFileConfigI::LoadBeamSimData ( const std::set< std::string > & filenames, const std::string & det_loc )

virtual

Reimplemented from genie::flux::GFluxFileConfigI.

Definition at line 49 of file GFluxFileConfigI.cxx.

  {
    // Loads a beam simulation root file into the GFluxFileConfig driver.
    // have a set<> want a vector<>
    std::vector<std::string> filevec;
    std::copy(fileset.begin(),fileset.end(),std::back_inserter(filevec));
    LoadBeamSimData(filevec,config); // call the one that takes a vector
  }

LoadBeamSimData() [2/3]

void genie::flux::GFluxFileConfigI::LoadBeamSimData ( const std::string & filename, const std::string & det_loc )

virtual

Reimplemented from genie::flux::GFluxFileConfigI.

Definition at line 52 of file GFluxFileConfigI.cxx.

  {
    // Loads a beam simulation root file into the GFluxFileConfig driver.
    std::vector<std::string> filevec;
    filevec.push_back(filename);
    LoadBeamSimData(filevec,config); // call the one that takes a vector
  }

LoadBeamSimData() [3/3]

void GNuMIFlux::LoadBeamSimData ( const std::vector< std::string > & filenames, const std::string & det_loc )

virtual

first is primary method for loading root flux ntuple files and config others are alternatives that can be overloaded but have sensible defaults to fall back to calling the vector version

Implements genie::flux::GFluxFileConfigI.

Definition at line 484 of file GNuMIFlux.cxx.

{
// Loads in a gnumi beam simulation root file (converted from hbook format)
// into the GNuMIFlux driver.
 
  bool found_cfg = this->LoadConfig(config);
  if ( ! found_cfg ) {
    LOG("Flux", pFATAL)
      << "LoadBeamSimData could not find XML config \"" << config << "\"\n";
    exit(1);
  }
 
  fNuFluxFilePatterns = patterns;
  std::vector<int> nfiles_from_pattern;
 
  // create a (sorted) set of file names
  // this also helps ensure that the same file isn't listed multiple times
  std::set<std::string> fnames;
 
  LOG("Flux", pINFO) << "LoadBeamSimData was passed " << patterns.size()
                     << " patterns";
 
  for (size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
    string pattern = patterns[ipatt];
    nfiles_from_pattern.push_back(0);
 
    LOG("Flux", pNOTICE)
      << "Loading gnumi flux tree from ROOT file pattern ["
      << std::setw(3) << ipatt << "] \"" << pattern << "\"";
 
    // !WILDCARD only works for file name ... NOT directory
    string dirname = gSystem->UnixPathName(gSystem->WorkingDirectory());
    size_t slashpos = pattern.find_last_of("/");
    size_t fbegin;
    if ( slashpos != std::string::npos ) {
      dirname = pattern.substr(0,slashpos);
      LOG("Flux", pINFO) << "Look for flux using directory " << dirname;
      fbegin = slashpos + 1;
    } else { fbegin = 0; }
 
    void* dirp = gSystem->OpenDirectory(gSystem->ExpandPathName(dirname.c_str()));
    if ( dirp ) {
      std::string basename =
      pattern.substr(fbegin,pattern.size()-fbegin);
      TRegexp re(basename.c_str(),kTRUE);
      const char* onefile;
      while ( ( onefile = gSystem->GetDirEntry(dirp) ) ) {
        TString afile = onefile;
        if ( afile=="." || afile==".." ) continue;
        if ( basename!=afile && afile.Index(re) == kNPOS ) continue;
        std::string fullname = dirname + "/" + afile.Data();
        fnames.insert(fullname);
        nfiles_from_pattern[ipatt]++;
      }
      gSystem->FreeDirectory(dirp);
    } // legal directory
  } // loop over patterns
 
  size_t indx = 0;
  std::set<string>::const_iterator sitr = fnames.begin();
  for ( ; sitr != fnames.end(); ++sitr, ++indx ) {
    string filename = *sitr;
    //std::cout << "  [" << std::setw(3) << indx << "]  \""
    //          << filename << "\"" << std::endl;
    bool isok = true;
    // this next test only works for local files, so we can't do that
    // if we want to stream via xrootd
    // ! (gSystem->AccessPathName(filename.c_str()));
    if ( isok ) {
      // open the file to see what it contains
      // h10    => g3numi _or_ flugg
      // nudata => g4numi
      // for now distinguish between g3numi/flugg using file name
      TFile* tf = TFile::Open(filename.c_str(),"READ");
      int isflugg = ( filename.find("flugg") != string::npos ) ? 1 : 0;
      const std::string tnames[] = { "h10", "nudata" };
      const std::string gnames[] = { "g3numi","g4numi","flugg","g4flugg"};
      for (int j = 0; j < 2 ; ++j ) {
        TTree* atree = (TTree*)tf->Get(tnames[j].c_str());
        if ( atree ) {
          const std::string tname_this = tnames[j];
          const std::string gname_this = gnames[j+2*isflugg];
          // create chain if none exists
          if ( ! fNuFluxTree ) {
            this->SetTreeName(tname_this);
            fNuFluxTree = new TChain(fNuFluxTreeName.c_str());
            fNuFluxGen = gname_this;
            // here we should scan for estimated POTs/file
            // also maybe the check on max weight
          }
          // sanity check for mixing g3/g4/flugg files
          if ( fNuFluxTreeName !=  tname_this ||
               fNuFluxGen      !=  gname_this    ) {
            LOG("Flux", pFATAL)
              << "Inconsistent flux file types\n"
              << "The input gnumi flux file \"" << filename
              << "\"\ncontains a '" << tname_this << "' " << gname_this
              << "numi ntuple "
              << "but a '" << fNuFluxTreeName << "' " << fNuFluxGen
              << " numi ntuple has alread been seen in the chain";
            exit(1);
          } // sanity mix/match g3/g4
            // add the file to the chain
          this->AddFile(atree,filename);
        } // found a tree
      } // loop over either g3 or g4
      tf->Close();
      delete tf;
    } // loop over tree type
  } // loop over sorted file names
 
  if ( fNuFluxTreeName == "" ) {
    LOG("Flux", pFATAL)
     << "The input gnumi flux file doesn't exist! Initialization failed!";
    exit(1);
  }
  if ( fNuFluxGen == "g3numi" ) fG3NuMI = new g3numi(fNuFluxTree);
  if ( fNuFluxGen == "g4numi" ) fG4NuMI = new g4numi(fNuFluxTree);
  if ( fNuFluxGen == "flugg"  ) fFlugg  = new flugg(fNuFluxTree);
 
  // this will open all files and read header!!
  fNEntries = fNuFluxTree->GetEntries();
 
  if ( fNEntries == 0 ) {
    LOG("Flux", pERROR)
      << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
    LOG("Flux", pERROR)
      << "Loaded flux tree contains " <<  fNEntries << " entries";
    LOG("Flux", pERROR)
      << "Was passed the file patterns: ";
    for (size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
      string pattern = patterns[ipatt];
      LOG("Flux", pERROR)
        << "  [" << std::setw(3) << ipatt <<"] " << pattern;
    }
    LOG("Flux", pERROR)
      << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
  } else {
    LOG("Flux", pNOTICE)
      << "Loaded flux tree contains " <<  fNEntries << " entries"
      << " from " << fnames.size() << " unique files";
    for (size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
      string pattern = patterns[ipatt];
      LOG("Flux", pINFO)
        << " pattern: " << pattern << " yielded "
        << nfiles_from_pattern[ipatt] << " files";
    }
  }
 
  // we have a file we can work with
  if (!fDetLocIsSet) {
     LOG("Flux", pERROR)
       << "LoadBeamSimData left detector location unset";
  }
  if (fMaxWeight<=0) {
     LOG("Flux", pINFO)
       << "Run ScanForMaxWeight() as part of LoadBeamSimData";
     this->ScanForMaxWeight();
  }
 
  // current ntuple cycle # (flux ntuples may be recycled)
  fICycle =  0;
  // pick a starting entry index [0:fNEntries-1]
  // pretend we just used up the the previous one
  RandomGen* rnd = RandomGen::Instance();
  fIUse   =  9999999;
  fIEntry = rnd->RndFlux().Integer(fNEntries) - 1;
 
  // don't count things we used to estimate max weight
  fSumWeight  = 0;
  fNNeutrinos = 0;
  fAccumPOTs  = 0;
 
  LOG("Flux",pNOTICE) << "about to CalcEffPOTsPerNu";
  this->CalcEffPOTsPerNu();
 
}

References AddFile(), CalcEffPOTsPerNu(), fAccumPOTs, fDetLocIsSet, fFlugg, fG3NuMI, fG4NuMI, genie::flux::GFluxFileConfigI::fICycle, fIEntry, fIUse, fMaxWeight, fNEntries, fNNeutrinos, fNuFluxFilePatterns, fNuFluxGen, fNuFluxTree, fNuFluxTreeName, fSumWeight, genie::RandomGen::Instance(), LoadConfig(), LOG, pERROR, pFATAL, pINFO, pNOTICE, genie::RandomGen::RndFlux(), ScanForMaxWeight(), and SetTreeName().

LoadConfig()

bool GNuMIFlux::LoadConfig

(

string

cfg

)

load a named configuration

Definition at line 2324 of file GNuMIFlux.cxx.

{
  const char* altxml = gSystem->Getenv("GNUMIFLUXXML");
  if ( altxml ) {
    SetXMLFileBase(altxml);
  }
  genie::flux::GNuMIFluxXMLHelper helper(this);
  return helper.LoadConfig(cfg);
}

References genie::flux::GNuMIFluxXMLHelper::LoadConfig(), and genie::flux::GFluxFileConfigI::SetXMLFileBase().

Referenced by LoadBeamSimData().

MaxEnergy()

double genie::flux::GNuMIFlux::MaxEnergy ( void )

inlinevirtual

declare the max flux neutrino energy that can be generated (for init. purposes)

Implements genie::GFluxI.

Definition at line 231 of file GNuMIFlux.h.

{ return  fMaxEv;               }

References fMaxEv.

Momentum()

const TLorentzVector & genie::flux::GNuMIFlux::Momentum ( void )

inlinevirtual

returns the flux neutrino 4-momentum

Implements genie::GFluxI.

Definition at line 235 of file GNuMIFlux.h.

{ return  fCurEntry->fgP4User;  }

References fCurEntry.

MoveToZ0()

void GNuMIFlux::MoveToZ0

(

double

z0

)

move ray origin to user coord Z0

Definition at line 403 of file GNuMIFlux.cxx.

{
  // move ray origin to specified user z0
  // move beam coord entry correspondingly
 
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("Flux", pNOTICE)
    << "MoveToZ0 (z0usr=" << z0usr << ") before:"
    << "\n p4 user: " << utils::print::P4AsShortString(&(fCurEntry->fgP4User))
    << "\n x4 user: " << utils::print::X4AsString(&(fCurEntry->fgX4User));
#endif
 
  double pzusr    = fCurEntry->fgP4User.Pz();
  if ( TMath::Abs(pzusr) < 1.0e-30 ) {
    // neutrino is moving almost entirely in x-y plane
    LOG("Flux", pWARN)
      << "MoveToZ0(" << z0usr << ") not possible due to pz_usr (" << pzusr << ")";
    return;
  }
 
  double scale = (z0usr - fCurEntry->fgX4User.Z()) / pzusr;
  fCurEntry->fgX4User += (scale*fCurEntry->fgP4User);
  fCurEntry->fgX4     += ((fLengthScaleU2B*scale)*fCurEntry->fgP4);
  // this scaling works for distances, but not the time component
  fCurEntry->fgX4.SetT(0);
  fCurEntry->fgX4User.SetT(0);
 
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("Flux", pNOTICE)
    << "MoveToZ0 (" << z0usr << ") after:"
    << "\n x4 user: " << utils::print::X4AsString(&(fCurEntry->fgX4User));
#endif
 
}

References fCurEntry, fLengthScaleU2B, LOG, genie::utils::print::P4AsShortString(), pNOTICE, pWARN, and genie::utils::print::X4AsString().

Referenced by GenerateNext_weighted().

NFluxNeutrinos()

long int GNuMIFlux::NFluxNeutrinos ( void ) const

virtual

of rays generated

number of flux neutrinos ray generated so far

Implements genie::flux::GFluxExposureI.

Definition at line 132 of file GNuMIFlux.cxx.

{
  /// number of flux neutrinos ray generated so far
  return fNNeutrinos;
}

References fNNeutrinos.

PassThroughInfo()

const GNuMIFluxPassThroughInfo & genie::flux::GNuMIFlux::PassThroughInfo ( void )

inline

GNuMIFluxPassThroughInfo.

Definition at line 252 of file GNuMIFlux.h.

References fCurEntry.

Referenced by GenerateNext().

PdgCode()

int genie::flux::GNuMIFlux::PdgCode ( void )

inlinevirtual

returns the flux neutrino pdg code

Implements genie::GFluxI.

Definition at line 233 of file GNuMIFlux.h.

{ return  fCurEntry->fgPdgC;    }

References fCurEntry.

Position()

const TLorentzVector & genie::flux::GNuMIFlux::Position ( void )

inlinevirtual

returns the flux neutrino 4-position (note: expect SI rather than physical units)

Implements genie::GFluxI.

Definition at line 236 of file GNuMIFlux.h.

{ return  fCurEntry->fgX4User;  }

References fCurEntry.

POT_curr()

double GNuMIFlux::POT_curr

(

void

)

current average POT (RWH?)

Definition at line 477 of file GNuMIFlux.cxx.

                               {
  // RWH: Not sure what POT_curr is supposed to represent I'll guess for
  // now that that it means what I mean by UsedPOTs().
  return UsedPOTs();
}

References UsedPOTs().

PrintConfig()

void GNuMIFlux::PrintConfig ( )

virtual

print the current configuration

Implements genie::flux::GFluxFileConfigI.

Definition at line 2336 of file GNuMIFlux.cxx.

{
 
  std::ostringstream s;
  PDGCodeList::const_iterator itr = fPdgCList->begin();
  for ( ; itr != fPdgCList->end(); ++itr) s << (*itr) << " ";
  s << "[rejected: ";
  itr = fPdgCListRej->begin();
  for ( ; itr != fPdgCListRej->end(); ++itr) s << (*itr) << " ";
  s << " ] ";
 
  std::ostringstream fpattout;
  for (size_t i = 0; i < fNuFluxFilePatterns.size(); ++i)
    fpattout << "\n [" << std::setw(3) << i << "] " << fNuFluxFilePatterns[i];
 
  std::ostringstream flistout;
  std::vector<std::string> flist = GetFileList();
  for (size_t i = 0; i < flist.size(); ++i)
    flistout << "\n [" << std::setw(3) << i << "] " << flist[i];
 
  TLorentzVector usr0(0,0,0,0);
  TLorentzVector usr0asbeam;
  User2BeamPos(usr0,usr0asbeam);
 
  const int w=10, p=6;
  std::ostringstream beamrot_str, beamrotinv_str;
  beamrot_str
    << "fBeamRot: " << std::setprecision(p) << "\n"
    << "  [ "
    << std::setw(w) << fBeamRot.XX() << " "
    << std::setw(w) << fBeamRot.XY() << " "
    << std::setw(w) << fBeamRot.XZ() << " ]\n"
    << "  [ "
    << std::setw(w) << fBeamRot.YX() << " "
    << std::setw(w) << fBeamRot.YY() << " "
    << std::setw(w) << fBeamRot.YZ() << " ]\n"
    << "  [ "
    << std::setw(w) << fBeamRot.ZX() << " "
    << std::setw(w) << fBeamRot.ZY() << " "
    << std::setw(w) << fBeamRot.ZZ() << " ]";
  beamrotinv_str
    << "fBeamRotInv: " << std::setprecision(p) << "\n"
    << "  [ "
    << std::setw(w) << fBeamRotInv.XX() << " "
    << std::setw(w) << fBeamRotInv.XY() << " "
    << std::setw(w) << fBeamRotInv.XZ() << " ]\n"
    << "  [ "
    << std::setw(w) << fBeamRotInv.YX() << " "
    << std::setw(w) << fBeamRotInv.YY() << " "
    << std::setw(w) << fBeamRotInv.YZ() << " ]\n"
    << "  [ "
    << std::setw(w) << fBeamRotInv.ZX() << " "
    << std::setw(w) << fBeamRotInv.ZY() << " "
    << std::setw(w) << fBeamRotInv.ZZ() << " ]";
 
  LOG("Flux", pNOTICE)
    << "GNuMIFlux Config:"
    << "\n Enu_max " << fMaxEv
    << "\n pdg-codes: " << s.str() << "\n "
    << (fG3NuMI?"g3numi":"")
    << (fG4NuMI?"g4numi":"")
    << (fFlugg?"flugg":"")
    << "/" << fNuFluxGen << " "
    << "(" << fNuFluxTreeName << "), " << fNEntries << " entries"
    << " (FilePOTs " << fFilePOTs << ") "
    <<  "in " << fNFiles << " files: "
    << flistout.str()
    << "\n from file patterns:"
    << fpattout.str()
    << "\n wgt max=" << fMaxWeight << " fudge=" << fMaxWgtFudge << " using "
    << fMaxWgtEntries << " entries"
    << "\n Z0 pushback " << fZ0
    << "\n used entry " << fIEntry << " " << fIUse << "/" << fNUse
    << " times, in " << fICycle << "/" << fNCycles << " cycles"
    << "\n SumWeight " << fSumWeight << " for " << fNNeutrinos << " neutrinos"
    << "\n EffPOTsPerNu " << fEffPOTsPerNu << " AccumPOTs " << fAccumPOTs
    << "\n GenWeighted \"" << (fGenWeighted?"true":"false") << ", "
    << "\", Detector location set \"" << (fDetLocIsSet?"true":"false") << "\""
    << "\n LengthUnits " << fLengthUnits << ", scale b2u " << fLengthScaleB2U
    << ", scale u2b " << fLengthScaleU2B
    << "\n User Flux Window: "
    << "\n       " << utils::print::Vec3AsString(&(fFluxWindowPtUser[0]))
    << "\n       " << utils::print::Vec3AsString(&(fFluxWindowPtUser[1]))
    << "\n       " << utils::print::Vec3AsString(&(fFluxWindowPtUser[2]))
    << "\n Flux Window (cm, beam coord): "
    << "\n  base " << utils::print::X4AsString(&fFluxWindowBase)
    << "\n  dir1 " << utils::print::X4AsString(&fFluxWindowDir1) << " len " << fFluxWindowLen1
    << "\n  dir2 " << utils::print::X4AsString(&fFluxWindowDir2) << " len " << fFluxWindowLen2
    << "\n  normal " << utils::print::Vec3AsString(&(fWindowNormal))
    << "\n User Beam Origin: "
    << "\n  base " << utils::print::X4AsString(&fBeamZero)
    << "\n " << beamrot_str.str() << " "
    << "\n Detector Origin (beam coord): "
    << "\n  base " << utils::print::X4AsString(&usr0asbeam)
    << "\n " << beamrotinv_str.str() << " "
    << "\n UseFluxAtDetCenter " << fUseFluxAtDetCenter;
 
}

References fAccumPOTs, fBeamRot, fBeamRotInv, fBeamZero, fDetLocIsSet, fEffPOTsPerNu, fFilePOTs, fFlugg, fFluxWindowBase, fFluxWindowDir1, fFluxWindowDir2, fFluxWindowLen1, fFluxWindowLen2, fFluxWindowPtUser, fG3NuMI, fG4NuMI, fGenWeighted, genie::flux::GFluxFileConfigI::fICycle, fIEntry, fIUse, fLengthScaleB2U, fLengthScaleU2B, fLengthUnits, fMaxEv, fMaxWeight, fMaxWgtEntries, fMaxWgtFudge, genie::flux::GFluxFileConfigI::fNCycles, fNEntries, fNFiles, fNNeutrinos, fNuFluxFilePatterns, fNuFluxGen, fNuFluxTreeName, fNUse, genie::flux::GFluxFileConfigI::fPdgCList, genie::flux::GFluxFileConfigI::fPdgCListRej, fSumWeight, fUseFluxAtDetCenter, fWindowNormal, genie::flux::GFluxFileConfigI::fZ0, GetFileList(), LOG, pNOTICE, User2BeamPos(), genie::utils::print::Vec3AsString(), and genie::utils::print::X4AsString().

PrintCurrent()

void GNuMIFlux::PrintCurrent

(

void

)

print current entry from leaves

Definition at line 1004 of file GNuMIFlux.cxx.

{
  LOG("Flux", pNOTICE) << "CurrentEntry:" << *fCurEntry;
}

References fCurEntry, LOG, and pNOTICE.

ResetCurrent()

void GNuMIFlux::ResetCurrent ( void )

private

Definition at line 1106 of file GNuMIFlux.cxx.

{
// reset running values of neutrino pdg-code, 4-position & 4-momentum
// and the input ntuple leaves
 
  fCurEntry->ResetCurrent();
  fCurEntry->ResetCopy();
}

References fCurEntry.

Referenced by GenerateNext_weighted(), and Initialize().

ScanForMaxWeight()

void GNuMIFlux::ScanForMaxWeight

(

void

)

scan for max flux weight (before generating unweighted flux neutrinos)

Definition at line 680 of file GNuMIFlux.cxx.

{
  if (!fDetLocIsSet) {
     LOG("Flux", pERROR)
       << "Specify a detector location before scanning for max weight";
     return;
  }
 
  // scan for the maximum weight
  int ipos_estimator = fUseFluxAtDetCenter;
  if ( ipos_estimator == 0 ) {
    // within 100m of a known point?
    double zbase = fFluxWindowBase.Z();
    if ( TMath::Abs(zbase-103648.837) < 10000. ) ipos_estimator = -1; // use NearDet
    if ( TMath::Abs(zbase-73534000. ) < 10000. ) ipos_estimator = +1; // use FarDet
  }
  if ( ipos_estimator != 0 ) {
 
    //// one can't really be sure which Nwtfar/Nwtnear this refers to
    //// some gnumi files have "NOvA" weights
    const char* ntwgtstrv[4] = { "Nimpwt*Nwtnear",
                                 "Nimpwt*Nwtfar",
                                 "Nimpwt*NWtNear[0]",
                                 "Nimpwt*NWtFar[0]"  };
    int strindx = 0;
    if ( ipos_estimator > 0 ) strindx = 1;
    if ( fG4NuMI ) strindx += 2;
    // set upper limit on how many entries to scan
    Long64_t nscan = TMath::Min(fNEntries,200000LL);
 
    fNuFluxTree->Draw(ntwgtstrv[strindx],"","goff",nscan);
    //std::cout << " Draw \"" << ntwgtstrv[strindx] << "\"" << std::endl;
    //std::cout << " SelectedRows " << fNuFluxTree->GetSelectedRows()
    //          << " V1 " << fNuFluxTree->GetV1() << std::endl;
 
    Long64_t idx = TMath::LocMax(fNuFluxTree->GetSelectedRows(),
                                 fNuFluxTree->GetV1());
    //std::cout << "idx " << idx << " of " << fNuFluxTree->GetSelectedRows() << std::endl;
    fMaxWeight = fNuFluxTree->GetV1()[idx];
    LOG("Flux", pNOTICE) << "Maximum flux weight from Nwt in ntuple = "
                         << fMaxWeight;
    if ( fMaxWeight <= 0 ) {
      LOG("Flux", pFATAL) << "Non-positive maximum flux weight!";
      exit(1);
    }
  }
  // the above works only for things close to the MINOS stored weight
  // values.  otherwise we need to work out our own estimate.
  double wgtgenmx = 0, enumx = 0;
  TStopwatch t;
  t.Start();
  for (int itry=0; itry < fMaxWgtEntries; ++itry) {
    this->GenerateNext_weighted();
    double wgt = this->Weight();
    if ( wgt > wgtgenmx ) wgtgenmx = wgt;
    double enu = fCurEntry->fgP4.Energy();
    if ( enu > enumx ) enumx = enu;
  }
  t.Stop();
  t.Print("u");
  LOG("Flux", pNOTICE) << "Maximum flux weight for spin = "
                       << wgtgenmx << ", energy = " << enumx
                       << " (" << fMaxWgtEntries << ")";
 
  if (wgtgenmx > fMaxWeight ) fMaxWeight = wgtgenmx;
  // apply a fudge factor to estimated weight
  fMaxWeight *= fMaxWgtFudge;
  // adjust max energy?
  if ( enumx*fMaxEFudge > fMaxEv ) {
    LOG("Flux", pNOTICE) << "Adjust max: was=" << fMaxEv
                         << " now " << enumx << "*" << fMaxEFudge
                         << " = " << enumx*fMaxEFudge;
    fMaxEv = enumx * fMaxEFudge;
  }
 
  LOG("Flux", pNOTICE) << "Maximum flux weight = " << fMaxWeight
                       << ", energy = " << fMaxEv;
 
}

References fCurEntry, fDetLocIsSet, fFluxWindowBase, fG4NuMI, fMaxEFudge, fMaxEv, fMaxWeight, fMaxWgtEntries, fMaxWgtFudge, fNEntries, fNuFluxTree, fUseFluxAtDetCenter, GenerateNext_weighted(), LOG, pERROR, pFATAL, pNOTICE, and Weight().

Referenced by LoadBeamSimData().

SetApplyWindowTiltWeight()

void genie::flux::GNuMIFlux::SetApplyWindowTiltWeight ( bool apply = true )

inline

Parameters

apply

apply wgt due to tilt of flux window relative to beam

Definition at line 302 of file GNuMIFlux.h.

            { fApplyTiltWeight = apply; }

References fApplyTiltWeight.

SetBeamCenter()

void GNuMIFlux::SetBeamCenter

(

TVector3

beam0

)

Definition at line 920 of file GNuMIFlux.cxx.

{
  // set coord transform between detector and beam
  // NOTE: internally these are in "cm", but user might have set a preference
  fBeamZero = TLorentzVector(beam0,0);  // no time shift
}

References fBeamZero.

SetBeamRotation()

void GNuMIFlux::SetBeamRotation

(

TRotation

beamrot

)

< beam (0,0,0) relative to user frame, beam direction in user frame

Definition at line 913 of file GNuMIFlux.cxx.

{
  // rotation is really only 3-d vector, but we'll be operating on LorentzV's
  fBeamRot    = TLorentzRotation(beamrot);
  fBeamRotInv = fBeamRot.Inverse();
}

References fBeamRot, and fBeamRotInv.

SetDefaults()

void GNuMIFlux::SetDefaults ( void )

private

Definition at line 1077 of file GNuMIFlux.cxx.

{
// - Set default neutrino species list (nue, nuebar, numu, numubar) and
//   maximum energy (120 GeV).
//   These defaults can be overwritten by user calls (at the driver init) to
//   GNuMIlux::SetMaxEnergy(double Ev) and
//   GNuMIFlux::SetFluxParticles(const PDGCodeList & particles)
// - Set the default file normalization to 1E+21 POT
// - Set the default flux neutrino start z position at -5m (z=0 is the
//   detector centre).
// - Set number of cycles to 1
 
  LOG("Flux", pNOTICE) << "Setting default GNuMIFlux driver options";
 
  PDGCodeList particles;
  particles.push_back(kPdgNuMu);
  particles.push_back(kPdgAntiNuMu);
  particles.push_back(kPdgNuE);
  particles.push_back(kPdgAntiNuE);
 
  this->SetFluxParticles (particles);
  this->SetMaxEnergy     (120./*GeV*/);  // was 200, but that would be wasteful
  this->SetUpstreamZ     (-3.4e38); // way upstream ==> use flux window
  this->SetNumOfCycles   (0);
  this->SetEntryReuse    (1);
 
  this->SetXMLFileBase("GNuMIFlux.xml");
}

References genie::kPdgAntiNuE, genie::kPdgAntiNuMu, genie::kPdgNuE, genie::kPdgNuMu, LOG, pNOTICE, genie::PDGCodeList::push_back(), SetEntryReuse(), genie::flux::GFluxFileConfigI::SetFluxParticles(), SetMaxEnergy(), genie::flux::GFluxFileConfigI::SetNumOfCycles(), genie::flux::GFluxFileConfigI::SetUpstreamZ(), and genie::flux::GFluxFileConfigI::SetXMLFileBase().

Referenced by Initialize().

SetEntryReuse()

void GNuMIFlux::SetEntryReuse

(

long int

nuse = 1

)

of times to use entry before moving to next

Definition at line 768 of file GNuMIFlux.cxx.

{
// With nuse > 1 then the same entry in the file is used "nuse" times
// before moving on to the next entry in the ntuple
 
  fNUse    = TMath::Max(1L, nuse);
}

References fNUse.

Referenced by SetDefaults().

SetFluxWindow() [1/2]

bool GNuMIFlux::SetFluxWindow	(	GNuMIFlux::StdFluxWindow_t	stdwindow,
		double	padding = 0 )

return false if unhandled

Definition at line 801 of file GNuMIFlux.cxx.

{
  // set some standard flux windows
  // rwh:  should also set detector coord transform
  // rwh:  padding allows add constant padding to pre-existing set
  double padbeam = padding / fLengthScaleB2U;  // user might set different units
  LOG("Flux",pNOTICE)
    << "SetBeamFluxWindow " << (int)stdwindow << " padding " << padbeam << " cm";
 
 
  switch ( stdwindow ) {
#ifdef THIS_IS_NOT_YET_IMPLEMENTED
  case kMinosNearDet:
    SetBeamFluxWindow(103648.837);
    break;
  case kMinosFarDear:
    SetBeamFluxWindow(73534000.);
    break;
  case kMinosNearRock:
    SetBeamFluxWindow();
    break;
  case kMinosFarRock:
    SetBeamFluxWindow();
    break;
#endif
  case kMinosNearCenter:
    {
      SetLengthUnits(genie::utils::units::UnitFromString("m"));
      fFluxWindowBase = kPosCenterNearBeam;
      fFluxWindowDir1 = TLorentzVector();  // no extent
      fFluxWindowDir2 = TLorentzVector();
      TLorentzVector usrpos;
      Beam2UserPos(fFluxWindowBase, usrpos);
      fFluxWindowPtUser[0] = usrpos.Vect();
      fFluxWindowPtUser[1] = fFluxWindowPtUser[0];
      fFluxWindowPtUser[2] = fFluxWindowPtUser[0];
      fFluxWindowLen1 = 0;
      fFluxWindowLen2 = 0;
      break;
    }
  case kMinosFarCenter:
    {
      SetLengthUnits(genie::utils::units::UnitFromString("m"));
      fFluxWindowBase = kPosCenterFarBeam;
      fFluxWindowDir1 = TLorentzVector();  // no extent
      fFluxWindowDir2 = TLorentzVector();
      TLorentzVector usrpos;
      Beam2UserPos(fFluxWindowBase, usrpos);
      fFluxWindowPtUser[0] = usrpos.Vect();
      fFluxWindowPtUser[1] = fFluxWindowPtUser[0];
      fFluxWindowPtUser[2] = fFluxWindowPtUser[0];
      fFluxWindowLen1 = 0;
      fFluxWindowLen2 = 0;
      break;
    }
  default:
    LOG("Flux", pERROR)
      << "SetBeamFluxWindow - StdFluxWindow " << stdwindow
      << " not yet implemented";
    return false;
  }
  LOG("Flux",pNOTICE) << "about to CalcEffPOTsPerNu";
  this->CalcEffPOTsPerNu();
  return true;
}

References Beam2UserPos(), CalcEffPOTsPerNu(), fFluxWindowBase, fFluxWindowDir1, fFluxWindowDir2, fFluxWindowLen1, fFluxWindowLen2, fFluxWindowPtUser, fLengthScaleB2U, kMinosFarCenter, kMinosFarRock, kMinosNearCenter, kMinosNearDet, kMinosNearRock, kPosCenterFarBeam(), LOG, pERROR, pNOTICE, SetLengthUnits(), and genie::utils::units::UnitFromString().

SetFluxWindow() [2/2]

void GNuMIFlux::SetFluxWindow	(	TVector3	p1,
		TVector3	p2,
		TVector3	p3 )

3 points define a plane (by default in user coordinates)

Definition at line 867 of file GNuMIFlux.cxx.

{
  // set flux window
  // NOTE: internally these are in "cm", but user might have set a preference
  fUseFluxAtDetCenter  = 0;
  fDetLocIsSet         = true;
 
  fFluxWindowPtUser[0] = p0;
  fFluxWindowPtUser[1] = p1;
  fFluxWindowPtUser[2] = p2;
 
  // convert from user to beam coord and from 3 points to base + 2 directions
  // apply units conversion
  TLorentzVector ptbm0, ptbm1, ptbm2;
  User2BeamPos(TLorentzVector(fFluxWindowPtUser[0],0),ptbm0);
  User2BeamPos(TLorentzVector(fFluxWindowPtUser[1],0),ptbm1);
  User2BeamPos(TLorentzVector(fFluxWindowPtUser[2],0),ptbm2);
 
  fFluxWindowBase = ptbm0;
  fFluxWindowDir1 = ptbm1 - ptbm0;
  fFluxWindowDir2 = ptbm2 - ptbm0;
 
  fFluxWindowLen1 = fFluxWindowDir1.Mag();
  fFluxWindowLen2 = fFluxWindowDir2.Mag();
  fWindowNormal = fFluxWindowDir1.Vect().Cross(fFluxWindowDir2.Vect()).Unit();
 
  double dot = fFluxWindowDir1.Dot(fFluxWindowDir2);
  if ( TMath::Abs(dot) > 1.0e-8 )
    LOG("Flux",pWARN) << "Dot product between window direction vectors was "
                      << dot << "; please check for orthoganality";
 
  LOG("Flux",pNOTICE) << "about to CalcEffPOTsPerNu";
  this->CalcEffPOTsPerNu();
}

References CalcEffPOTsPerNu(), fDetLocIsSet, fFluxWindowBase, fFluxWindowDir1, fFluxWindowDir2, fFluxWindowLen1, fFluxWindowLen2, fFluxWindowPtUser, fUseFluxAtDetCenter, fWindowNormal, LOG, pNOTICE, pWARN, and User2BeamPos().

SetGenWeighted()

void genie::flux::GNuMIFlux::SetGenWeighted ( bool genwgt = false )

inline

toggle whether GenerateNext() returns weight=1 flux (initial default false)

Definition at line 292 of file GNuMIFlux.h.

References fGenWeighted.

SetLengthUnits()

void GNuMIFlux::SetLengthUnits

(

double

user_units

)

Set units assumed by user.

Definition at line 1206 of file GNuMIFlux.cxx.

{
  // Set the scale factor for lengths going from beam (cm) to user coordinates
 
  // GNuMIFlux uses "cm" as the length unit consistently internally (this is
  // the length units used by both the g3 and g4 ntuples).  User interactions
  // setting the beam-to-detector coordinate transform, flux window, and the
  // returned position might need to be in other units.  Use:
  //     double scale = genie::utils::units::UnitFromString("cm");
  // ( #include "Utils/UnitUtils.h for declaration )
  // to get the correct scale factor to pass in.
 
  double rescale = fLengthUnits / user_units;
  fLengthUnits = user_units;
  double cm = genie::utils::units::UnitFromString("cm");
  fLengthScaleB2U = cm / user_units;
  fLengthScaleU2B = user_units / cm;
 
  // in case we're changing units without resetting transform/window
  // not recommended, but should work
  fCurEntry->fgX4User  *= rescale;
  fBeamZero            *= rescale;
  fFluxWindowPtUser[0] *= rescale;
  fFluxWindowPtUser[1] *= rescale;
  fFluxWindowPtUser[2] *= rescale;
 
  // case GNuMIFlux::kmeter:  fLengthScaleB2U =   0.01  ; break;
  // case GNuMIFlux::kcm:     fLengthScaleB2U =   1.    ; break;
  // case GNuMIFlux::kmm:     fLengthScaleB2U =  10.    ; break;
  // case GNuMIFlux::kfm:     fLengthScaleB2U =   1.e13 ; break;  // 10e-2m -> 10e-15m
 
}

References fBeamZero, fCurEntry, fFluxWindowPtUser, fLengthScaleB2U, fLengthScaleU2B, fLengthUnits, and genie::utils::units::UnitFromString().

Referenced by Initialize(), SetFluxWindow(), UseFluxAtFarDetCenter(), and UseFluxAtNearDetCenter().

SetMaxEFudge()

void genie::flux::GNuMIFlux::SetMaxEFudge ( double fudge = 1.05 )

inline

Parameters

fudge

extra fudge factor in estimating maximum energy

Definition at line 300 of file GNuMIFlux.h.

            { fMaxEFudge = fudge; }

References fMaxEFudge.

SetMaxEnergy()

void GNuMIFlux::SetMaxEnergy

(

double

Ev

)

specify maximum flx neutrino energy

Definition at line 760 of file GNuMIFlux.cxx.

{
  fMaxEv = TMath::Max(0.,Ev);
 
  LOG("Flux", pINFO)
    << "Declared maximum flux neutrino energy: " << fMaxEv;
}

References fMaxEv, LOG, and pINFO.

Referenced by SetDefaults().

SetMaxWgtScan()

void genie::flux::GNuMIFlux::SetMaxWgtScan ( double fudge = 1.05, long int nentries = 2500000 )

inline

Parameters

nentries

configuration when estimating max weight

Definition at line 298 of file GNuMIFlux.h.

            { fMaxWgtFudge = fudge; fMaxWgtEntries = nentries; }

References fMaxWgtEntries, and fMaxWgtFudge.

SetTreeName()

void GNuMIFlux::SetTreeName

(

string

name

)

set input tree name (default: "h10")

Definition at line 776 of file GNuMIFlux.cxx.

{
  fNuFluxTreeName = name;
}

References fNuFluxTreeName.

Referenced by LoadBeamSimData().

SumWeight()

double genie::flux::GNuMIFlux::SumWeight ( void ) const

inline

integrated weight for flux neutrinos looped so far

Definition at line 267 of file GNuMIFlux.h.

References fSumWeight.

UsedPOTs()

double GNuMIFlux::UsedPOTs

(

void

)

const

of protons-on-target used

Definition at line 464 of file GNuMIFlux.cxx.

{
// Compute current number of flux POTs
 
  if (!fNuFluxTree) {
     LOG("Flux", pWARN)
          << "The flux driver has not been properly configured";
     return 0;
  }
  return fAccumPOTs;
}

References fAccumPOTs, fNuFluxTree, LOG, and pWARN.

Referenced by GetTotalExposure(), and POT_curr().

UseFluxAtFarDetCenter()

void GNuMIFlux::UseFluxAtFarDetCenter

(

void

)

Definition at line 791 of file GNuMIFlux.cxx.

{
  SetLengthUnits(genie::utils::units::UnitFromString("m"));
  fFluxWindowBase      = kPosCenterFarBeam;
  fFluxWindowDir1      = TLorentzVector();  // no extent
  fFluxWindowDir2      = TLorentzVector();
  fUseFluxAtDetCenter  = +1;
  fDetLocIsSet         = true;
}

References fDetLocIsSet, fFluxWindowBase, fFluxWindowDir1, fFluxWindowDir2, fUseFluxAtDetCenter, kPosCenterFarBeam(), SetLengthUnits(), and genie::utils::units::UnitFromString().

UseFluxAtNearDetCenter()

void GNuMIFlux::UseFluxAtNearDetCenter

(

void

)

force weights at MINOS detector "center" as found in ntuple

Definition at line 781 of file GNuMIFlux.cxx.

{
  SetLengthUnits(genie::utils::units::UnitFromString("m"));
  fFluxWindowBase      = kPosCenterNearBeam;
  fFluxWindowDir1      = TLorentzVector();   // no extent
  fFluxWindowDir2      = TLorentzVector();
  fUseFluxAtDetCenter  = -1;
  fDetLocIsSet         = true;
}

References fDetLocIsSet, fFluxWindowBase, fFluxWindowDir1, fFluxWindowDir2, fUseFluxAtDetCenter, SetLengthUnits(), and genie::utils::units::UnitFromString().

User2BeamDir()

void GNuMIFlux::User2BeamDir	(	const TLorentzVector &	usrdir,
		TLorentzVector &	beamdir ) const

Definition at line 992 of file GNuMIFlux.cxx.

{
  beamdir = fLengthScaleU2B*(fBeamRotInv*usrdir);
}

References fBeamRotInv, and fLengthScaleU2B.

User2BeamP4()

void GNuMIFlux::User2BeamP4	(	const TLorentzVector &	usrp4,
		TLorentzVector &	beamp4 ) const

Definition at line 997 of file GNuMIFlux.cxx.

{
  beamp4 = fBeamRotInv*usrp4;
}

References fBeamRotInv.

User2BeamPos()

void GNuMIFlux::User2BeamPos	(	const TLorentzVector &	usrxyz,
		TLorentzVector &	beamxyz ) const

Definition at line 987 of file GNuMIFlux.cxx.

{
  beamxyz = fLengthScaleU2B*(fBeamRotInv*(usrxyz-fBeamZero));
}

References fBeamRotInv, fBeamZero, and fLengthScaleU2B.

Referenced by PrintConfig(), and SetFluxWindow().

Weight()

double genie::flux::GNuMIFlux::Weight ( void )

inlinevirtual

returns the flux neutrino weight (if any)

Implements genie::GFluxI.

Definition at line 234 of file GNuMIFlux.h.

{ return  fWeight;              }

References fWeight.

Referenced by GenerateNext(), GenerateNext_weighted(), and ScanForMaxWeight().

Member Data Documentation

fAccumPOTs

double genie::flux::GNuMIFlux::fAccumPOTs

private

POTs used so far.

Definition at line 411 of file GNuMIFlux.h.

Referenced by Clear(), GenerateNext_weighted(), Initialize(), LoadBeamSimData(), PrintConfig(), and UsedPOTs().

fApplyTiltWeight

bool genie::flux::GNuMIFlux::fApplyTiltWeight

private

wgt due to window normal not || beam

Definition at line 414 of file GNuMIFlux.h.

Referenced by GenerateNext_weighted(), Initialize(), and SetApplyWindowTiltWeight().

fBeamRot

TLorentzRotation genie::flux::GNuMIFlux::fBeamRot

private

rotation applied beam --> user coord

Definition at line 423 of file GNuMIFlux.h.

Referenced by Beam2UserDir(), Beam2UserP4(), Beam2UserPos(), GetBeamRotation(), PrintConfig(), and SetBeamRotation().

fBeamRotInv

TLorentzRotation genie::flux::GNuMIFlux::fBeamRotInv

private

Definition at line 424 of file GNuMIFlux.h.

Referenced by PrintConfig(), SetBeamRotation(), User2BeamDir(), User2BeamP4(), and User2BeamPos().

fBeamZero

TLorentzVector genie::flux::GNuMIFlux::fBeamZero

private

beam origin in user coords

Definition at line 422 of file GNuMIFlux.h.

Referenced by Beam2UserPos(), GetBeamCenter(), PrintConfig(), SetBeamCenter(), SetLengthUnits(), and User2BeamPos().

fCurEntry

GNuMIFluxPassThroughInfo* genie::flux::GNuMIFlux::fCurEntry

private

copy of current ntuple entry info (owned structure)

Definition at line 436 of file GNuMIFlux.h.

Referenced by CleanUp(), GenerateNext(), GenerateNext_weighted(), GetBranchInfo(), GetDecayDist(), Initialize(), Momentum(), MoveToZ0(), PassThroughInfo(), PdgCode(), Position(), PrintCurrent(), ResetCurrent(), ScanForMaxWeight(), and SetLengthUnits().

fDetLocIsSet

bool genie::flux::GNuMIFlux::fDetLocIsSet

private

is a flux location (near/far) set?

Definition at line 415 of file GNuMIFlux.h.

Referenced by Initialize(), LoadBeamSimData(), PrintConfig(), ScanForMaxWeight(), SetFluxWindow(), UseFluxAtFarDetCenter(), and UseFluxAtNearDetCenter().

fEffPOTsPerNu

double genie::flux::GNuMIFlux::fEffPOTsPerNu

private

what a entry is worth ...

Definition at line 410 of file GNuMIFlux.h.

Referenced by CalcEffPOTsPerNu(), GenerateNext_weighted(), Initialize(), and PrintConfig().

fEnd

bool genie::flux::GNuMIFlux::fEnd

private

end condition reached

Definition at line 385 of file GNuMIFlux.h.

Referenced by End(), GenerateNext_weighted(), and Initialize().

fFilePOTs

Long64_t genie::flux::GNuMIFlux::fFilePOTs

private

of protons-on-target represented by all files

Definition at line 398 of file GNuMIFlux.h.

Referenced by AddFile(), CalcEffPOTsPerNu(), Initialize(), and PrintConfig().

fFlugg

flugg* genie::flux::GNuMIFlux::fFlugg

private

flugg ntuple

Definition at line 393 of file GNuMIFlux.h.

Referenced by CleanUp(), GenerateNext_weighted(), Initialize(), LoadBeamSimData(), and PrintConfig().

fFluxWindowBase

TLorentzVector genie::flux::GNuMIFlux::fFluxWindowBase

private

base point for flux window - beam coord

Definition at line 427 of file GNuMIFlux.h.

Referenced by GenerateNext_weighted(), PrintConfig(), ScanForMaxWeight(), SetFluxWindow(), SetFluxWindow(), UseFluxAtFarDetCenter(), and UseFluxAtNearDetCenter().

fFluxWindowDir1

TLorentzVector genie::flux::GNuMIFlux::fFluxWindowDir1

private

extent for flux window (direction 1)

Definition at line 428 of file GNuMIFlux.h.

Referenced by CalcEffPOTsPerNu(), GenerateNext_weighted(), PrintConfig(), SetFluxWindow(), SetFluxWindow(), UseFluxAtFarDetCenter(), and UseFluxAtNearDetCenter().

fFluxWindowDir2

TLorentzVector genie::flux::GNuMIFlux::fFluxWindowDir2

private

extent for flux window (direction 2)

Definition at line 429 of file GNuMIFlux.h.

Referenced by CalcEffPOTsPerNu(), GenerateNext_weighted(), PrintConfig(), SetFluxWindow(), SetFluxWindow(), UseFluxAtFarDetCenter(), and UseFluxAtNearDetCenter().

fFluxWindowLen1

double genie::flux::GNuMIFlux::fFluxWindowLen1

private

Definition at line 430 of file GNuMIFlux.h.

Referenced by PrintConfig(), SetFluxWindow(), and SetFluxWindow().

fFluxWindowLen2

double genie::flux::GNuMIFlux::fFluxWindowLen2

private

Definition at line 431 of file GNuMIFlux.h.

Referenced by PrintConfig(), SetFluxWindow(), and SetFluxWindow().

fFluxWindowPtUser

TVector3 genie::flux::GNuMIFlux::fFluxWindowPtUser[3]

private

user points of flux window

Definition at line 426 of file GNuMIFlux.h.

Referenced by GetFluxWindow(), PrintConfig(), SetFluxWindow(), SetFluxWindow(), and SetLengthUnits().

fG3NuMI

g3numi* genie::flux::GNuMIFlux::fG3NuMI

private

g3numi ntuple

Definition at line 391 of file GNuMIFlux.h.

Referenced by CleanUp(), GenerateNext_weighted(), Initialize(), LoadBeamSimData(), and PrintConfig().

fG4NuMI

g4numi* genie::flux::GNuMIFlux::fG4NuMI

private

g4numi ntuple

Definition at line 392 of file GNuMIFlux.h.

Referenced by CleanUp(), GenerateNext_weighted(), Initialize(), LoadBeamSimData(), PrintConfig(), and ScanForMaxWeight().

fGenWeighted

bool genie::flux::GNuMIFlux::fGenWeighted

private

does GenerateNext() give weights?

Definition at line 413 of file GNuMIFlux.h.

Referenced by GenerateNext(), GenerateWeighted(), Initialize(), PrintConfig(), and SetGenWeighted().

fgX4dkvtx

TLorentzVector genie::flux::GNuMIFlux::fgX4dkvtx

private

decay 4-position beam coord

Definition at line 434 of file GNuMIFlux.h.

Referenced by GenerateNext_weighted(), and GetDecayDist().

fIEntry

Long64_t genie::flux::GNuMIFlux::fIEntry

private

current flux ntuple entry

Definition at line 396 of file GNuMIFlux.h.

Referenced by CleanUp(), GenerateNext_weighted(), GetEntryNumber(), Index(), Initialize(), LoadBeamSimData(), and PrintConfig().

fIUse

long int genie::flux::GNuMIFlux::fIUse

private

current # of times an entry has been used

Definition at line 407 of file GNuMIFlux.h.

Referenced by CleanUp(), GenerateNext_weighted(), Initialize(), LoadBeamSimData(), and PrintConfig().

fLengthScaleB2U

double genie::flux::GNuMIFlux::fLengthScaleB2U

private

scale factor beam (cm) --> user

Definition at line 419 of file GNuMIFlux.h.

Referenced by Beam2UserDir(), Beam2UserPos(), GetDecayDist(), LengthUnits(), PrintConfig(), SetFluxWindow(), and SetLengthUnits().

fLengthScaleU2B

double genie::flux::GNuMIFlux::fLengthScaleU2B

private

scale factor beam user --> (cm)

Definition at line 420 of file GNuMIFlux.h.

Referenced by MoveToZ0(), PrintConfig(), SetLengthUnits(), User2BeamDir(), and User2BeamPos().

fLengthUnits

double genie::flux::GNuMIFlux::fLengthUnits

private

units for coord in user exchanges

Definition at line 418 of file GNuMIFlux.h.

Referenced by PrintConfig(), and SetLengthUnits().

fMaxEFudge

double genie::flux::GNuMIFlux::fMaxEFudge

private

fudge factor for estmating max enu (0=> use fixed 120GeV)

Definition at line 404 of file GNuMIFlux.h.

Referenced by Initialize(), ScanForMaxWeight(), and SetMaxEFudge().

fMaxEv

double genie::flux::GNuMIFlux::fMaxEv

private

maximum energy

Definition at line 384 of file GNuMIFlux.h.

Referenced by GenerateNext_weighted(), Initialize(), MaxEnergy(), PrintConfig(), ScanForMaxWeight(), and SetMaxEnergy().

fMaxWeight

double genie::flux::GNuMIFlux::fMaxWeight

private

max flux neutrino weight in input file

Definition at line 401 of file GNuMIFlux.h.

Referenced by GenerateNext(), GenerateNext_weighted(), Initialize(), LoadBeamSimData(), PrintConfig(), and ScanForMaxWeight().

fMaxWgtEntries

long int genie::flux::GNuMIFlux::fMaxWgtEntries

private

of entries in estimating max wgt

Definition at line 403 of file GNuMIFlux.h.

Referenced by Initialize(), PrintConfig(), ScanForMaxWeight(), and SetMaxWgtScan().

fMaxWgtFudge

double genie::flux::GNuMIFlux::fMaxWgtFudge

private

fudge factor for estimating max wgt

Definition at line 402 of file GNuMIFlux.h.

Referenced by GenerateNext(), Initialize(), PrintConfig(), ScanForMaxWeight(), and SetMaxWgtScan().

fNEntries

Long64_t genie::flux::GNuMIFlux::fNEntries

private

number of flux ntuple entries

Definition at line 395 of file GNuMIFlux.h.

Referenced by CalcEffPOTsPerNu(), GenerateNext_weighted(), Initialize(), LoadBeamSimData(), PrintConfig(), and ScanForMaxWeight().

fNFiles

int genie::flux::GNuMIFlux::fNFiles

private

number of files in chain

Definition at line 394 of file GNuMIFlux.h.

Referenced by AddFile(), Initialize(), and PrintConfig().

fNNeutrinos

long int genie::flux::GNuMIFlux::fNNeutrinos

private

number of flux neutrinos thrown so far

Definition at line 409 of file GNuMIFlux.h.

Referenced by Clear(), GenerateNext_weighted(), Initialize(), LoadBeamSimData(), NFluxNeutrinos(), and PrintConfig().

fNuFluxFilePatterns

std::vector<string> genie::flux::GNuMIFlux::fNuFluxFilePatterns

private

(potentially wildcarded) path(s)

Definition at line 387 of file GNuMIFlux.h.

Referenced by LoadBeamSimData(), and PrintConfig().

fNuFluxGen

string genie::flux::GNuMIFlux::fNuFluxGen

private

"g3numi" "g4numi" or "flugg"

Definition at line 390 of file GNuMIFlux.h.

Referenced by AddFile(), Initialize(), LoadBeamSimData(), and PrintConfig().

fNuFluxTree

TChain* genie::flux::GNuMIFlux::fNuFluxTree

private

TTree in g3numi or g4numi // REF ONLY!

Definition at line 389 of file GNuMIFlux.h.

Referenced by AddFile(), CalcEffPOTsPerNu(), GetFileList(), Initialize(), LoadBeamSimData(), ScanForMaxWeight(), and UsedPOTs().

fNuFluxTreeName

string genie::flux::GNuMIFlux::fNuFluxTreeName

private

Tree name "h10" (g3) or "nudata" (g4)

Definition at line 388 of file GNuMIFlux.h.

Referenced by AddFile(), Initialize(), LoadBeamSimData(), PrintConfig(), and SetTreeName().

fNUse

long int genie::flux::GNuMIFlux::fNUse

private

how often to use same entry in a row

Definition at line 406 of file GNuMIFlux.h.

Referenced by CleanUp(), Initialize(), PrintConfig(), and SetEntryReuse().

fNuTot

Long64_t genie::flux::GNuMIFlux::fNuTot

private

cummulative # of entries (=fNEntries)

Definition at line 397 of file GNuMIFlux.h.

Referenced by AddFile(), and Initialize().

fSumWeight

double genie::flux::GNuMIFlux::fSumWeight

private

sum of weights for nus thrown so far

Definition at line 408 of file GNuMIFlux.h.

Referenced by Clear(), GenerateNext_weighted(), Initialize(), LoadBeamSimData(), PrintConfig(), and SumWeight().

fUseFluxAtDetCenter

int genie::flux::GNuMIFlux::fUseFluxAtDetCenter

private

use flux at near (-1) or far (+1) det center from ntuple?

Definition at line 416 of file GNuMIFlux.h.

Referenced by GenerateNext_weighted(), Initialize(), PrintConfig(), ScanForMaxWeight(), SetFluxWindow(), UseFluxAtFarDetCenter(), and UseFluxAtNearDetCenter().

fWeight

double genie::flux::GNuMIFlux::fWeight

private

current neutrino weight, =1 if generating unweighted entries

Definition at line 400 of file GNuMIFlux.h.

Referenced by GenerateNext(), GenerateNext_weighted(), and Weight().

fWindowNormal

TVector3 genie::flux::GNuMIFlux::fWindowNormal

private

normal direction for flux window

Definition at line 432 of file GNuMIFlux.h.

Referenced by FluxWindowNormal(), PrintConfig(), and SetFluxWindow().

---

The documentation for this class was generated from the following files:

• GNuMIFlux.h
• GNuMIFlux.cxx