155 GJPARCNuFlux::~GJPARCNuFlux()
168 if(end)
return false;
172 if(!nextok)
continue;
176 <<
"Got flux entry: " << this->
Index()
177 <<
" - Cycle: "<<
fICycle <<
"/ infinite";
180 <<
"Got flux entry: "<< this->
Index()
188 <<
"Curr flux neutrino fractional weight = " <<
f;
191 <<
"** Fractional weight = " << f <<
" > 1 !!";
193 double r = (f < 1.) ? rnd->
RndFlux().Rndm() : 0;
200 <<
"** Rejecting current flux neutrino based on the flux weight only";
216 <<
"The flux driver has not been properly configured";
227 <<
"The input jnubeam flux ntuple contains no entries for detector id " 239 <<
"No more entries in input flux neutrino ntuple";
269 LOG(
"Flux",
pERROR) <<
"Negative flux weight! Will set weight to 0.0";
279 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 281 <<
"Current flux neutrino not at specified detector location";
339 " --> unable to infer neutrino pdg! Aborting job!";
351 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 353 <<
"Current flux neutrino " 354 <<
"not at the list of neutrinos to be considered at this job.";
368 <<
"Flux neutrino energy exceeds declared maximum neutrino energy";
379 fgP4.SetPxPyPzE (pxnu, pynu, pznu, Enu);
393 fgX4.SetXYZT (xnu, ynu, znu, 0.);
395 fgX4.SetXYZT (0.,0.,0.,0.);
399 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 401 <<
"Generated neutrino: " 402 <<
"\n pdg-code: " <<
fgPdgC 418 if (start_pos == std::string::npos) start_pos = 0;
else ++start_pos;
433 <<
"The flux driver has not been properly configured";
454 <<
"The flux driver has not been properly configured";
490 <<
"Loading jnubeam flux tree from ROOT file: " <<
filename;
492 <<
"Detector location: " << detector_location;
500 string fileroot =
"";
501 int firstfile = -1, lastfile = -1;
504 if (filenamev.size() != 3) {
506 <<
"Flux filename should be specfied as either:\n" 507 <<
"\t For a single input file: /dir/root_filename.#.root\n" 508 <<
"\t For multiple input files: /dir/root_filename@#@#";
511 fileroot = filenamev[0];
512 firstfile = atoi(filenamev[1].c_str());
513 lastfile = atoi(filenamev[2].c_str());
515 <<
"Chaining beam simulation output files with stem: " << fileroot
516 <<
" and run numbers in the range: [" << firstfile <<
", " << firstfile <<
"]";
520 bool is_accessible = ! (gSystem->AccessPathName( filename.c_str() ));
521 if (!is_accessible) {
523 <<
"The input jnubeam flux file doesn't exist! Initialization failed!";
528 bool please_exit =
false;
529 for (
int i = firstfile;
i < lastfile+1;
i++) {
530 bool is_accessible = ! (gSystem->AccessPathName( Form(
"%s.%i.root",fileroot.c_str(),
i) ));
531 if (!is_accessible) {
533 <<
"The input jnubeam flux file " << Form(
"%s.%i.root",fileroot.c_str(),
i)
534 <<
"doesn't exist! Initialization failed!";
547 <<
" ** Unknown input detector location: " <<
fDetLoc;
555 string ntuple_name = (
fIsNDLoc) ?
"h3002" :
"h2000";
560 <<
"Could not find tree h3002 in file " << Form(
"%s.%i.root",fileroot.c_str(),firstfile)
561 <<
" Trying tree h3001";
563 ntuple_name =
"h3001";
564 fNuFluxChain =
new TChain(ntuple_name.c_str());
565 result = fNuFluxChain->Add( Form(
"%s.%i.root",fileroot.c_str(),firstfile), 0);
569 <<
"** Couldn't get flux tree: " << ntuple_name;
573 for (
int i = firstfile+1;
i < lastfile+1;
i++) {
574 result =
fNuFluxChain->Add( Form(
"%s.%i.root",fileroot.c_str(),
i), 0 );
577 <<
"** Couldn't get flux tree " << ntuple_name <<
" in file " << Form(
"%s.%i.root",fileroot.c_str(),
i);
586 string ntuple_name = (
fIsNDLoc) ?
"h3002" :
"h2000";
589 ntuple_name =
"h3001";
593 <<
"Getting flux tree: " << ntuple_name;
596 <<
"** Couldn't get flux tree: " << ntuple_name;
608 <<
"Loaded flux tree contains " <<
fNEntries <<
" entries";
611 <<
"Getting tree branches & setting leaf addresses";
615 bool missing_critical =
false;
622 LOG(
"Flux",
pFATAL) <<
"Cannot find flux branch: norm";
623 missing_critical =
true;
628 LOG(
"Flux",
pFATAL) <<
"Cannot find flux branch: Enu";
629 missing_critical =
true;
634 LOG(
"Flux",
pFATAL) <<
"Cannot find flux branch: ppid";
635 missing_critical =
true;
640 LOG(
"Flux",
pFATAL) <<
"Cannot find flux branch: mode";
641 missing_critical =
true;
646 LOG(
"Flux",
pFATAL) <<
"Cannot find flux branch: rnu";
647 missing_critical =
true;
652 LOG(
"Flux",
pFATAL) <<
"Cannot find flux branch: xnu";
653 missing_critical =
true;
658 LOG(
"Flux",
pFATAL) <<
"Cannot find flux branch: ynu";
659 missing_critical =
true;
664 LOG(
"Flux",
pFATAL) <<
"Cannot find flux branch: nnu";
665 missing_critical =
true;
670 LOG(
"Flux",
pFATAL) <<
"Cannot find flux branch: idfd";
671 missing_critical =
true;
674 if(missing_critical){
676 <<
"Unable to find critical information in the flux ntuple! Initialization failed!";
765 for (
int i = firstfile+1;
i < lastfile+1;
i++) {
792 for(
int ientry = 0; ientry <
fNEntries; ientry++) {
799 LOG(
"Flux",
pERROR) <<
"Negative flux weight! Will set weight to 0.0";
814 <<
"The input jnubeam flux ntuple contains no entries for detector id " 822 LOG(
"Flux",
pFATAL) <<
"Non-positive maximum flux weight!";
845 <<
"Declared list of neutrino species: " << *
fPdgCList;
850 fMaxEv = TMath::Max(0.,Ev);
853 <<
"Declared maximum flux neutrino energy: " <<
fMaxEv;
906 LOG(
"Flux",
pERROR) <<
"Setting flux driver to start looping over entries " 907 <<
"with offset of "<<
offset;
919 if(std::strcmp(opt,
"CycleHistory") == 0){
933 LOG(
"Flux",
pNOTICE) <<
"Initializing GJPARCNuFlux driver";
985 LOG(
"Flux",
pNOTICE) <<
"Setting default GJPARCNuFlux driver options";
1008 fgP4.SetPxPyPzE (0.,0.,0.,0.);
1009 fgX4.SetXYZT (0.,0.,0.,0.);
1035 if(name ==
"sk" )
return -1;
1038 for (
int i=1;
i<51;
i++) {
1039 temp.Form(
"nd%d",
i);
1040 if(name == temp.Data())
return i;
1072 for(
int i = 0;
i<3;
i++){
1107 for(
int i = 0;
i < 2;
i++){
1136 for(
int i=0;
i<3;
i++){
1170 for(
int i = 0;
i < 2;
i++){
1178 for(
int i = 0;
i < 3;
i++)
hcur[
i] = -999999.;
1186 stream <<
"\n idfd = " << info.
idfd 1187 <<
"\n norm = " << info.
norm 1190 <<
"\n Enu = " << info.
Enu 1191 <<
"\n geant code = " << info.
ppid 1193 <<
"\n decay mode = " << info.
mode 1194 <<
"\n nvtx0 = " << info.
nvtx0 1195 <<
"\n |momentum| @ decay = " << info.
ppi 1196 <<
"\n position_vector @ decay = (" 1197 << info.
xpi[0] <<
", " 1198 << info.
xpi[1] <<
", " 1199 << info.
xpi[2] <<
")" 1200 <<
"\n direction_vector @ decay = (" 1201 << info.
npi[0] <<
", " 1202 << info.
npi[1] <<
", " 1203 << info.
npi[2] <<
")" 1204 <<
"\n |momentum| @ prod. = " << info.
ppi0 1205 <<
"\n position_vector @ prod. = (" 1206 << info.
xpi0[0] <<
", " 1207 << info.
xpi0[1] <<
", " 1208 << info.
xpi0[2] <<
")" 1209 <<
"\n direction_vector @ prod. = (" 1210 << info.
npi0[0] <<
", " 1211 << info.
npi0[1] <<
", " 1212 << info.
npi0[2] <<
")" 1213 <<
"\n cospibm = " << info.
cospibm 1214 <<
"\n cospi0bm = " << info.
cospi0bm 1215 <<
"\n Plus additional info if flux version is later than 07a"
TTree * fNuFluxTree
input flux ntuple
const XML_Char XML_Encoding * info
double POT_curravg(void)
current average POT
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 fgPdgC
running generated nu pdg-code
string P4AsShortString(const TLorentzVector *p)
double fSumWeightTot1c
total sum of weights for neutrinos to be thrown / cycle
bool LoadBeamSimData(string filename, string det_loc)
load a jnubeam root flux ntuple
THE MAIN GENIE PROJECT NAMESPACE
bool End(void)
true if no more flux nu's can be thrown (eg reaching end of beam sim ntuples)
bool fLoadedNeutrino
set to true when GenerateNext_weighted has been called successfully
void SetUpstreamZ(double z0)
set flux neutrino initial z position (upstream of the detector)
int DLocName2Id(string name)
static RandomGen * Instance()
Access instance.
string fDetLoc
input detector location ('sk','nd1','nd2',...)
TChain * fNuFluxSumChain
input summary ntuple for flux file. This tree is only present for later flux versions > 10a ...
double POT_1cycle(void)
flux POT per cycle
bool GenerateNext_weighted(void)
bool fNuFluxUsingTree
are we using a TTree or a TChain to view the input flux file?
double fFilePOT
file POT normalization, typically 1E+21
TChain * fNuFluxChain
input flux ntuple
double Weight(void)
returns the flux neutrino weight (if any)
bool ExistsInPDGCodeList(int pdg_code) const
void GenerateWeighted(bool gen_weighted=true)
set whether to generate weighted or unweighted neutrinos
long int fOffset
start looping at entry fOffset
long int fNNeutrinos
number of flux neutrinos thrown so far
void SetFluxParticles(const PDGCodeList &particles)
specify list of flux neutrino species
A singleton holding random number generator classes. All random number generation in GENIE should tak...
GJPARCNuFluxPassThroughInfo()
void Initialize(Bool_t useTMVAStyle=kTRUE)
TTree * fNuFluxSumTree
input summary ntuple for flux file. This tree is only present for later flux versions > 10a ...
A GENIE flux driver encapsulating the JPARC neutrino flux. It reads-in the official JPARC neutrino fl...
void SetMaxEnergy(double Ev)
specify maximum flx neutrino energy
int GeantToPdg(int geant_code)
long int fNNeutrinosTot1c
total number of flux neutrinos to be thrown / cycle
TLorentzVector fgP4
running generated nu 4-momentum
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
int fDetLocId
input detector location id (fDetLoc -> jnubeam idfd)
ClassImp(GJPARCNuFluxPassThroughInfo) GJPARCNuFlux
void SetFilePOT(double pot)
flux file norm is in /N POT/det [ND] or /N POT/cm^2 [FD]. Specify N (typically 1E+21) ...
int fNCycles
how many times to cycle through the flux ntuple
static const double kASmallNum
int fNDetLocIdFound
per cycle keep track of the number of fDetLocId are found - if this is zero will exit job ...
long int fIEntry
current flux ntuple entry
bool fIsNDLoc
input location is a 'near' detector location?
GJPARCNuFluxPassThroughInfo * fPassThroughInfo
bool fGenerateWeighted
generate weighted/deweighted flux neutrinos (default is false)
TFile * fNuFluxFile
input flux file
double fMaxWeight
max flux neutrino weight in input file for the specified detector location
void Clear(Option_t *opt)
reset state variables based on opt
vector< string > Split(string input, string delim)
friend ostream & operator<<(ostream &stream, const GJPARCNuFluxPassThroughInfo &info)
long int Index(void)
returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current ent...
bool fIsFDLoc
input location is a 'far' detector location?
double fZ0
configurable starting z position for each flux neutrino (in detector coord system) ...
PDGCodeList * fPdgCList
list of neutrino pdg-codes
fvar< T > floor(const fvar< T > &x)
bool GenerateNext(void)
generate the next flux neutrino (return false in err)
long int fNEntries
number of flux ntuple entries
void SetNumOfCycles(int n)
set how many times to cycle through the ntuple (default: 1 / n=0 means 'infinite') ...
assert(nhit_max >=nhit_nbins)
TLorentzVector fgX4
running generated nu 4-position
bool fUseRandomOffset
whether set random starting point when looping over flux ntuples
void RandomOffset(void)
choose a random offset as starting entry in flux ntuple
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
string X4AsString(const TLorentzVector *x)
double fNorm
current flux ntuple normalisation
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
#define FLUXDRIVERREG4(_nsa, _nsb, _name, _fqname)
void Copy(const PDGCodeList &list)
copy / print
double fSumWeight
sum of weights for neutrinos thrown so far
long int fEntriesThisCycle
keep track of number of entries used so far for this cycle
void push_back(int pdg_code)
double fMaxEv
maximum energy