32 #include <TChainElement.h> 34 #include <TStopwatch.h> 84 GSimpleNtpFlux::~GSimpleNtpFlux()
110 LOG(
"Flux",
pNOTICE) <<
"GenerateNext signaled End() ";
117 if ( ! nextok )
continue;
139 <<
"** Fractional weight = " << f
140 <<
" > 1 !! Bump fMaxWeight estimate to " <<
fMaxWeight 144 double r = (f < 1.) ? rnd->
RndFlux().Rndm() : 0;
145 bool accept = ( r <
f );
148 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 150 <<
"Generated beam neutrino: " 174 <<
"The flux driver has not been properly configured";
183 if ( fIUse < fNUse && fIEntry >= 0 ) {
200 <<
"No more entries in input flux neutrino ntuple, cycle " 212 #ifdef USE_INDEX_FOR_META 213 int nbmeta =
fNuMetaTree->GetEntryWithIndex(metakey);
222 for (
int imeta = 0; imeta < nmeta; ++imeta ) {
229 LOG(
"Flux",
pERROR) <<
"Failed to find right metakey=" << metakey
230 <<
" (was " << oldkey <<
") out of " << nmeta
234 LOG(
"Flux",
pDEBUG) <<
"Get meta " << metakey
235 <<
" (was " << oldkey <<
") " 237 <<
" nb " << nbytes <<
" " << nbmeta;
238 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 242 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 246 <<
" ifile " << ifile <<
" nbytes " << nbytes
279 <<
"Encountered neutrino specie (" << badpdg
280 <<
") that wasn't in SetFluxParticles() list, " 281 <<
"\nDeclared list of neutrino species: " << *
fPdgCList;
297 <<
"Flux neutrino energy exceeds declared maximum neutrino energy" 298 <<
"\nEv = " << Ev <<
"(> Ev{max} = " <<
fMaxEv <<
")";
304 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 306 <<
"Generated neutrino: " <<
fIEntry 318 <<
"Generated neutrino had E_nu = " << Ev <<
" > " <<
fMaxEv 338 double pzusr =
fP4.Pz();
339 if ( TMath::Abs(pzusr) < 1.0
e-30 ) {
342 <<
"MoveToZ0(" << z0usr <<
") not possible due to pz_usr (" << pzusr <<
")";
346 double scale = (z0usr -
fX4.Z()) / pzusr;
378 <<
"The flux driver has not been properly configured";
391 std::vector<int> nfiles_from_pattern;
395 std::set<std::string> fnames;
397 LOG(
"Flux",
pINFO) <<
"LoadBeamSimData was passed " << patterns.size()
400 for (
size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
401 string pattern = patterns[ipatt];
402 nfiles_from_pattern.push_back(0);
404 <<
"Loading flux tree from ROOT file pattern [" 405 << std::setw(3) << ipatt <<
"] \"" << pattern <<
"\"";
408 string dirname = gSystem->UnixPathName(gSystem->WorkingDirectory());
409 size_t slashpos = pattern.find_last_of(
"/");
411 if ( slashpos != std::string::npos ) {
412 dirname = pattern.substr(0,slashpos);
413 LOG(
"Flux",
pDEBUG) <<
"Look for flux using directory " << dirname;
414 fbegin = slashpos + 1;
415 }
else { fbegin = 0; }
417 void* dirp = gSystem->OpenDirectory(gSystem->ExpandPathName(dirname.c_str()));
420 pattern.substr(fbegin,pattern.size()-fbegin);
421 TRegexp re(basename.c_str(),kTRUE);
423 while ( ( onefile = gSystem->GetDirEntry(dirp) ) ) {
424 TString afile = onefile;
425 if ( afile==
"." || afile==
".." )
continue;
426 if ( basename!=afile && afile.Index(re) == kNPOS )
continue;
427 std::string fullname = dirname +
"/" + afile.Data();
428 fnames.insert(fullname);
429 nfiles_from_pattern[ipatt]++;
431 gSystem->FreeDirectory(dirp);
436 std::set<string>::const_iterator sitr = fnames.begin();
437 for ( ; sitr != fnames.end(); ++sitr, ++indx) {
445 if ( ! isok )
continue;
449 TFile*
tf = TFile::Open(filename.c_str(),
"READ");
450 TTree* etree = (TTree*)tf->Get(
"flux");
452 TTree* mtree = (TTree*)tf->Get(
"meta");
454 LOG(
"Flux",
pDEBUG) <<
"AddFile " << filename
455 <<
" etree " << etree <<
" meta " << mtree;
456 this->
AddFile(etree,mtree,filename);
468 <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
470 <<
"Loaded flux tree contains " <<
fNEntries <<
" entries";
472 <<
"Was passed the file patterns: ";
473 for (
size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
474 string pattern = patterns[ipatt];
476 <<
" [" << std::setw(3) << ipatt <<
"] " <<
pattern;
479 <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
482 <<
"Loaded flux tree contains " <<
fNEntries <<
" entries" 483 <<
" from " << fnames.size() <<
" unique files";
484 for (
size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
485 string pattern = patterns[ipatt];
487 <<
" pattern: " << pattern <<
" yielded " 488 << nfiles_from_pattern[ipatt] <<
" files";
492 int sba_status[3] = { -999, -999, -999 };
494 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,26,0) 498 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,26,0) 499 if ( sba_status[0] < 0 ) {
501 <<
"flux chain has no \"entry\" branch " << sba_status[0];
509 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,26,0) 518 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,26,0) 527 <<
" SetBranchAddress status: " 528 <<
" \"entry\"=" << sba_status[0]
529 <<
" \"numi\"=" << sba_status[1]
530 <<
" \"aux\"=" << sba_status[2];
536 <<
"Run ProcessMeta() as part of LoadBeamSimData";
547 if ( config.find(
"no-offset-index") != string::npos ) {
548 LOG(
"Flux",
pINFO) <<
"Config saw \"no-offset-index\"";
559 LOG(
"Flux",
pDEBUG) <<
"about to CalcEffPOTsPerNu";
565 std::vector<std::string>& branchClassNames,
566 std::vector<void**>& branchObjPointers)
574 branchNames.push_back(
"simple");
575 branchClassNames.push_back(
"genie::flux::GSimpleNtpEntry");
576 branchObjPointers.push_back((
void**)&
fCurEntry);
580 branchNames.push_back(
"numi");
581 branchClassNames.push_back(
"genie::flux::GSimpleNtpNuMI");
582 branchObjPointers.push_back((
void**)&
fCurNuMI);
586 branchNames.push_back(
"aux");
587 branchClassNames.push_back(
"genie::flux::GSimpleNtpAux");
588 branchObjPointers.push_back((
void**)&
fCurAux);
599 double minwgt = +1.0e10;
600 double maxwgt = -1.0e10;
607 #ifdef USE_INDEX_FOR_META 609 LOG(
"Flux",
pDEBUG) <<
"ProcessMeta() BuildIndex nindices " << nindices;
612 for (
int imeta = 0; imeta < nmeta; ++imeta ) {
614 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 615 LOG(
"Flux",
pNOTICE) <<
"ProcessMeta() ifile " << imeta
619 minwgt = TMath::Min(minwgt,fCurMeta->minWgt);
620 maxwgt = TMath::Max(maxwgt,fCurMeta->maxWgt);
621 maxenu = TMath::Max(maxenu,fCurMeta->maxEnergy);
623 for (
size_t i = 0;
i < fCurMeta->pdglist.size(); ++
i)
632 LOG(
"Flux",
pFATAL) <<
"ProcessMeta() not all files have metadata";
636 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 638 <<
", energy = " <<
fMaxEv 649 fMaxEv = TMath::Max(0.,Ev);
652 <<
"Declared maximum flux neutrino energy: " <<
fMaxEv;
660 fNUse = TMath::Max(1
L, nuse);
689 LOG(
"Flux",
pWARN) <<
"GSimpleNtpFlux::Clear(" << opt <<
") called";
709 LOG(
"Flux",
pINFO) <<
"Initializing GSimpleNtpFlux driver";
762 LOG(
"Flux",
pINFO) <<
"Setting default GSimpleNtpFlux driver options";
782 LOG(
"Flux",
pINFO) <<
"Cleaning up...";
803 if ( ! fluxtree )
return;
805 ULong64_t
nentries = fluxtree->GetEntries();
808 if ( metatree )
fNuMetaTree->AddFile(fname.c_str());
812 <<
"flux->AddFile() of " << nentries
813 <<
" " << ((metatree)?
"[+meta]":
"[no-meta]")
814 <<
" [status=" << stat <<
"]" 815 <<
" entries in file: " <<
fname;
817 if ( stat != 1 )
return;
830 <<
"no request for \"" << name <<
"\" branch ";
834 if ( (
fNuFluxTree->GetBranch(name.c_str()) ) )
return true;
837 <<
"no \"" << name <<
"\" branch in the \"flux\" tree";
870 tpx = tpy = tpz = 0.;
872 pdpx = pdpy = pdpz = 0.;
873 pppx = pppy = pppz = 0.;
947 if (
pdglist[
i] == nupdg) found =
true;
948 if ( ! found )
pdglist.push_back(nupdg);
978 stream <<
"\nGSimpleNtpEntry " 979 <<
" PDG " << entry.
pdg 980 <<
" wgt " << entry.
wgt 981 <<
" ( metakey " << entry.
metakey <<
" )" 982 <<
"\n vtx [" << entry.
vtxx <<
"," << entry.
vtxy <<
"," 983 << entry.
vtxz <<
"] dist " << entry.
dist 984 <<
"\n p4 [" << entry.
px <<
"," << entry.
py <<
"," 985 << entry.
pz <<
"," << entry.
E <<
"]";
993 stream <<
"\nGSimpleNtpNuMI " 994 <<
"run " << numi.
run 995 <<
" evtno " << numi.
evtno 997 <<
"\n ndecay " << numi.
ndecay <<
" ptype " << numi.
ptype 999 <<
"\n tp[xyz] [" << numi.
tpx <<
"," << numi.
tpy <<
"," << numi.
tpz <<
"]" 1000 <<
"\n v[xyz] [" << numi.
vx <<
"," << numi.
vy <<
"," << numi.
vz <<
"]" 1001 <<
"\n pd[xyz] [" << numi.
pdpx <<
"," << numi.
pdpy <<
"," << numi.
pdpz <<
"]" 1002 <<
"\n pp[xyz] [" << numi.
pppx <<
"," << numi.
pppy <<
"," << numi.
pppz <<
"]" 1010 stream <<
"\nGSimpleNtpAux ";
1011 size_t nInt = aux.
auxint.size();
1012 stream <<
"\n ints: ";
1013 for (
size_t ijInt=0; ijInt < nInt; ++ijInt)
1014 stream <<
" " << aux.
auxint[ijInt];
1015 size_t nDbl = aux.
auxdbl.size();
1016 stream <<
"\n doubles: ";
1017 for (
size_t ijDbl=0; ijDbl < nDbl; ++ijDbl)
1018 stream <<
" " << aux.
auxdbl[ijDbl];
1026 size_t nf = meta.
pdglist.size();
1027 stream <<
"\nGSimpleNtpMeta " << nf <<
" flavors: ";
1028 for (
size_t i=0;
i<nf; ++
i) stream <<
" " << meta.
pdglist[
i];
1034 stream <<
"\n maxEnergy " << meta.
maxEnergy 1036 <<
" protons " << meta.
protons 1037 <<
" metakey " << meta.
metakey 1038 <<
"\n windowBase [" << meta.
windowBase[0] <<
"," 1040 <<
"\n windowDir1 [" << meta.
windowDir1[0] <<
"," 1042 <<
"\n windowDir2 [" << meta.
windowDir2[0] <<
"," 1046 if ( nInt > 0 ) stream <<
"\n aux ints: ";
1047 for (
size_t ijInt=0; ijInt < nInt; ++ijInt)
1051 if ( nDbl > 0 ) stream <<
"\n aux doubles: ";
1052 for (
size_t ijDbl=0; ijDbl < nDbl; ++ijDbl)
1056 stream <<
"\n " <<
nfiles <<
" input files: ";
1057 UInt_t nprint = TMath::Min(UInt_t(nfiles),
1059 for (UInt_t ifiles=0; ifiles < nprint; ++ifiles)
1060 stream <<
"\n " << meta.
infiles[ifiles];
1062 stream <<
"\n input seed: " << meta.
seed;
1076 std::ostringstream
s;
1077 PDGCodeList::const_iterator itr = fPdgCList->begin();
1078 for ( ; itr != fPdgCList->end(); ++itr) s << (*itr) <<
" ";
1080 itr = fPdgCListRej->begin();
1081 for ( ; itr != fPdgCListRej->end(); ++itr) s << (*itr) <<
" ";
1084 std::ostringstream fpattout;
1085 for (
size_t i = 0;
i < fNuFluxFilePatterns.size(); ++
i)
1086 fpattout <<
"\n [" << std::setw(3) <<
i <<
"] " << fNuFluxFilePatterns[
i];
1088 std::ostringstream flistout;
1089 std::vector<std::string>
flist = GetFileList();
1090 for (
size_t i = 0;
i < flist.size(); ++
i)
1091 flistout <<
"\n [" << std::setw(3) <<
i <<
"] " << flist[
i];
1094 <<
"GSimpleNtpFlux Config:" 1095 <<
"\n Enu_max " << fMaxEv
1096 <<
"\n pdg-codes: " << s.str() <<
"\n " 1097 <<
"\"flux\" " << fNEntries <<
" entries, " 1098 <<
"\"meta\" " << fNFiles <<
" entries" 1099 <<
" (FilePOTs " << fFilePOTs <<
") in files:" 1101 <<
"\n from file patterns: " 1103 <<
"\n wgt max=" << fMaxWeight
1104 <<
"\n Z0 pushback " << fZ0
1105 <<
"\n used entry " << fIEntry <<
" " << fIUse <<
"/" << fNUse
1106 <<
" times, in " << fICycle <<
"/" << fNCycles <<
" cycles" 1107 <<
"\n SumWeight " << fSumWeight <<
" for " << fNNeutrinos <<
" neutrinos" 1108 <<
" with " << fNEntriesUsed <<
" entries read" 1109 <<
"\n EffPOTsPerNu " << fEffPOTsPerNu <<
" AccumPOTs " << fAccumPOTs
1110 <<
"\n GenWeighted \"" << (fGenWeighted?
"true":
"false") <<
"\"" 1111 <<
" AlreadyUnwgt \"" << (fAlreadyUnwgt?
"true":
"false") <<
"\"" 1112 <<
" AllFilesMeta \"" << (fAllFilesMeta?
"true":
"false") <<
"\"";
1119 std::vector<std::string>
flist;
1120 TObjArray *fileElements=fNuFluxTree->GetListOfFiles();
1121 TIter
next(fileElements);
1122 TChainElement *chEl=0;
1123 while (( chEl=(TChainElement*)
next() )) {
1124 flist.push_back(chEl->GetTitle());
double UsedPOTs(void) const
of protons-on-target used
void GetFluxWindow(TVector3 &p1, TVector3 &p2, TVector3 &p3) const
3 points define a plane in beam coordinate
Double_t E
energy in lab frame
GSimpleNtpAux * fCurAux
current "aux" branch extra info
A class for generating concrete GFluxI derived classes based on the factory pattern. This code supplies a CPP macro which allows the classes to self-register and thus no modification of this class is needed in order to expand the list of classes it knows about.
Int_t fIFileNumber
which file for the current entry
void Clear(Option_t *opt)
reset state variables based on opt
virtual long int NFluxNeutrinos() const
of rays generated
std::vector< string > fNuFluxFilePatterns
(potentially wildcarded) path(s)
string P4AsShortString(const TLorentzVector *p)
long int fNNeutrinos
number of flux neutrinos thrown so far
void AddFile(TTree *fluxtree, TTree *metatree, string fname)
void PrintCurrent(void)
print current entry from leaves
THE MAIN GENIE PROJECT NAMESPACE
GSimpleNtpMeta * fCurMeta
current meta data
Double_t pdpx
nu parent momentum at time of decay
void SetEntryReuse(long int nuse=1)
of times to use entry before moving to next
long int fIUse
current # of times an entry has been used
Double_t px
x momentum in lab frame
static RandomGen * Instance()
Access instance.
A GENIE flux driver using a simple ntuple format.
double fSumWeight
sum of weights for nus thrown so far
Double_t vtxy
y position in lab frame
Double_t vx
vertex position of hadron/muon decay
void CalcEffPOTsPerNu(void)
bool ExistsInPDGCodeList(int pdg_code) const
Int_t tptype
parent particle type at target exit
ClassImp(EventRecord) namespace genie
bool End(void)
true if no more flux nu's can be thrown (eg reaching end of beam sim ntuples)
std::vector< Double_t > auxdbl
additional doubles associated w/ entry
virtual void SetNumOfCycles(long int ncycle)
limit cycling through input files
A singleton holding random number generator classes. All random number generation in GENIE should tak...
std::vector< Int_t > auxint
additional ints associated w/ entry
void Initialize(Bool_t useTMVAStyle=kTRUE)
virtual TTree * GetMetaDataTree()
void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
void MoveToZ0(double z0)
move ray origin to user coord Z0
virtual void SetUpstreamZ(double z0)
const genie::flux::GSimpleNtpEntry * GetCurrentEntry(void)
GSimpleNtpEntry.
Double_t pppx
nu parent momentum at production point
void ProcessMeta(void)
scan for max flux energy, weight
bool OptionalAttachBranch(std::string bname)
PDGCodeList * fPdgCList
list of neutrino pdg-codes to generate
std::vector< std::string > GetFileList()
list of files currently part of chain
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
TLorentzVector fX4
reconstituted position vector
static constexpr double L
bool fAllFilesMeta
do all files in chain have meta data
Double_t fFilePOTs
of protons-on-target represented by all files
bool fGenWeighted
does GenerateNext() give weights?
Double_t vtxz
z position in lab frame
GENIE interface for uniform flux exposure iterface.
Int_t ppmedium
tracking medium where parent was produced
virtual double GetTotalExposure() const
GFluxExposureI interface.
GSimpleNtpNuMI * fCurNuMICopy
current "numi" branch extra info
long int fNEntriesUsed
number of entries read from files
bool GenerateNext(void)
generate the next flux neutrino (return false in err)
long int fNCycles
times to cycle through the ntuple(s)
void PrintConfig()
print the current configuration
Double_t tpx
parent particle px at target exit
void SetMaxEnergy(double Ev)
specify maximum flx neutrino energy
GSimpleNtpNuMI * fCurNuMI
current "numi" branch extra info
int fNFiles
number of files in chain
virtual void GetBranchInfo(std::vector< std::string > &branchNames, std::vector< std::string > &branchClassNames, std::vector< void ** > &branchObjPointers)
double fMaxWeight
max flux neutrino weight in input file
PDGCodeList * fPdgCListRej
list of nu pdg-codes seen but rejected
double fWeight
current neutrino weight
TChain * fNuMetaTree
TTree // REF ONLY.
GSimpleNtpEntry * fCurEntry
current entry
TLorentzVector fP4
reconstituted p4 vector
bool fAlreadyUnwgt
are input files already unweighted
Long64_t fIEntry
current flux ntuple entry
Long64_t fNEntries
number of flux ntuple entries
Double_t vtxx
x position in lab frame
Double_t pz
z momentum in lab frame
GSimpleNtpAux * fCurAuxCopy
current "aux" branch extra info
long int fNUse
how often to use same entry in a row
assert(nhit_max >=nhit_nbins)
Double_t dist
distance from hadron decay
GSimpleNtpEntry * fCurEntryCopy
current entry
double fEffPOTsPerNu
what a entry is worth ...
double GetDecayDist() const
dist (user units) from dk to current pos
string fNuFluxBranchRequest
list of requested branches "entry,numi,au"
bool fEnd
end condition reached
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
string X4AsString(const TLorentzVector *x)
virtual void LoadBeamSimData(const std::vector< string > &filenames, const std::string &det_loc)
double fMaxEv
maximum energy
UInt_t metakey
key to meta data
void Print(const Option_t *opt="") const
Double_t py
y momentum in lab frame
double Weight(void)
returns the flux neutrino weight (if any)
double fAccumPOTs
POTs used so far.
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
#define FLUXDRIVERREG4(_nsa, _nsb, _name, _fqname)
void Print(const Option_t *opt="") const
void push_back(int pdg_code)
TChain * fNuFluxTree
TTree // REF ONLY.
bool GenerateNext_weighted(void)
void Print(const Option_t *opt="") const
Int_t ptype
parent type (PDG)