GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
genie::flux::GJPARCNuFlux Class Reference

A GENIE flux driver encapsulating the JPARC neutrino flux. It reads-in the official JPARC neutrino flux ntuples. More...

#include <GJPARCNuFlux.h>

Inheritance diagram for genie::flux::GJPARCNuFlux:
Inheritance graph
[legend]
Collaboration diagram for genie::flux::GJPARCNuFlux:
Collaboration graph
[legend]

Public Member Functions

 GJPARCNuFlux ()
 
 ~GJPARCNuFlux ()
 
const PDGCodeListFluxParticles (void)
 declare list of flux neutrinos that can be generated (for init. purposes) More...
 
double MaxEnergy (void)
 declare the max flux neutrino energy that can be generated (for init. purposes) More...
 
bool GenerateNext (void)
 generate the next flux neutrino (return false in err) More...
 
int PdgCode (void)
 returns the flux neutrino pdg code More...
 
double Weight (void)
 returns the flux neutrino weight (if any) More...
 
const TLorentzVector & Momentum (void)
 returns the flux neutrino 4-momentum More...
 
const TLorentzVector & Position (void)
 returns the flux neutrino 4-position (note: expect SI rather than physical units) More...
 
bool End (void)
 true if no more flux nu's can be thrown (eg reaching end of beam sim ntuples) More...
 
long int Index (void)
 returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current entry number) More...
 
void Clear (Option_t *opt)
 reset state variables based on opt More...
 
void GenerateWeighted (bool gen_weighted=true)
 set whether to generate weighted or unweighted neutrinos More...
 
bool LoadBeamSimData (string filename, string det_loc)
 load a jnubeam root flux ntuple More...
 
void SetFluxParticles (const PDGCodeList &particles)
 specify list of flux neutrino species More...
 
void SetMaxEnergy (double Ev)
 specify maximum flx neutrino energy More...
 
void SetFilePOT (double pot)
 flux file norm is in /N POT/det [ND] or /N POT/cm^2 [FD]. Specify N (typically 1E+21) More...
 
void SetUpstreamZ (double z0)
 set flux neutrino initial z position (upstream of the detector) More...
 
void SetNumOfCycles (int n)
 set how many times to cycle through the ntuple (default: 1 / n=0 means 'infinite') More...
 
void DisableOffset (void)
 switch off random offset, must be called before LoadBeamSimData to have any effect More...
 
void RandomOffset (void)
 choose a random offset as starting entry in flux ntuple More...
 
double POT_1cycle (void)
 flux POT per cycle More...
 
double POT_curravg (void)
 current average POT More...
 
long int NFluxNeutrinos (void) const
 number of flux neutrinos looped so far More...
 
double SumWeight (void) const
 intergated weight for flux neutrinos looped so far More...
 
const GJPARCNuFluxPassThroughInfoPassThroughInfo (void)
 GJPARCNuFluxPassThroughInfo. More...
 
- Public Member Functions inherited from genie::GFluxI
virtual ~GFluxI ()
 

Private Member Functions

bool GenerateNext_weighted (void)
 
void Initialize (void)
 
void SetDefaults (void)
 
void CleanUp (void)
 
void ResetCurrent (void)
 
int DLocName2Id (string name)
 

Private Attributes

double fMaxEv
 maximum energy More...
 
PDGCodeListfPdgCList
 list of neutrino pdg-codes More...
 
int fgPdgC
 running generated nu pdg-code More...
 
TLorentzVector fgP4
 running generated nu 4-momentum More...
 
TLorentzVector fgX4
 running generated nu 4-position More...
 
TFile * fNuFluxFile
 input flux file More...
 
TTree * fNuFluxTree
 input flux ntuple More...
 
TChain * fNuFluxChain
 input flux ntuple More...
 
TTree * fNuFluxSumTree
 input summary ntuple for flux file. This tree is only present for later flux versions > 10a More...
 
TChain * fNuFluxSumChain
 input summary ntuple for flux file. This tree is only present for later flux versions > 10a More...
 
bool fNuFluxUsingTree
 are we using a TTree or a TChain to view the input flux file? More...
 
string fDetLoc
 input detector location ('sk','nd1','nd2',...) More...
 
int fDetLocId
 input detector location id (fDetLoc -> jnubeam idfd) More...
 
int fNDetLocIdFound
 per cycle keep track of the number of fDetLocId are found - if this is zero will exit job More...
 
bool fIsFDLoc
 input location is a 'far' detector location? More...
 
bool fIsNDLoc
 input location is a 'near' detector location? More...
 
long int fNEntries
 number of flux ntuple entries More...
 
long int fIEntry
 current flux ntuple entry More...
 
long int fEntriesThisCycle
 keep track of number of entries used so far for this cycle More...
 
long int fOffset
 start looping at entry fOffset More...
 
double fNorm
 current flux ntuple normalisation More...
 
double fMaxWeight
 max flux neutrino weight in input file for the specified detector location More...
 
double fFilePOT
 file POT normalization, typically 1E+21 More...
 
double fZ0
 configurable starting z position for each flux neutrino (in detector coord system) More...
 
int fNCycles
 how many times to cycle through the flux ntuple More...
 
int fICycle
 current cycle More...
 
double fSumWeight
 sum of weights for neutrinos thrown so far More...
 
long int fNNeutrinos
 number of flux neutrinos thrown so far More...
 
double fSumWeightTot1c
 total sum of weights for neutrinos to be thrown / cycle More...
 
long int fNNeutrinosTot1c
 total number of flux neutrinos to be thrown / cycle More...
 
bool fGenerateWeighted
 generate weighted/deweighted flux neutrinos (default is false) More...
 
bool fUseRandomOffset
 whether set random starting point when looping over flux ntuples More...
 
bool fLoadedNeutrino
 set to true when GenerateNext_weighted has been called successfully More...
 
GJPARCNuFluxPassThroughInfofPassThroughInfo
 

Additional Inherited Members

- Protected Member Functions inherited from genie::GFluxI
 GFluxI ()
 

Detailed Description

A GENIE flux driver encapsulating the JPARC neutrino flux. It reads-in the official JPARC neutrino flux ntuples.

References:
See http://jnusrv01.kek.jp/internal/t2k/nubeam/flux/ (Note: T2K internal)
Author
Costas Andreopoulos <c.andreopoulos cern.ch> University of Liverpool
Created:
Feb 04, 2008
License:
Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 50 of file GJPARCNuFlux.h.

Constructor & Destructor Documentation

GJPARCNuFlux::GJPARCNuFlux ( )

Definition at line 45 of file GJPARCNuFlux.cxx.

References Initialize().

46 {
47  this->Initialize();
48 }
GJPARCNuFlux::~GJPARCNuFlux ( )

Definition at line 50 of file GJPARCNuFlux.cxx.

References CleanUp().

51 {
52  this->CleanUp();
53 }

Member Function Documentation

void GJPARCNuFlux::CleanUp ( void  )
private

Definition at line 909 of file GJPARCNuFlux.cxx.

References fNuFluxFile, fPassThroughInfo, fPdgCList, LOG, and pNOTICE.

Referenced by ~GJPARCNuFlux().

910 {
911  LOG("Flux", pNOTICE) << "Cleaning up...";
912 
913  if (fPdgCList) delete fPdgCList;
915 
916  if (fNuFluxFile) {
917  fNuFluxFile->Close();
918  delete fNuFluxFile;
919  }
920 }
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
GJPARCNuFluxPassThroughInfo * fPassThroughInfo
Definition: GJPARCNuFlux.h:143
TFile * fNuFluxFile
input flux file
Definition: GJPARCNuFlux.h:114
PDGCodeList * fPdgCList
list of neutrino pdg-codes
Definition: GJPARCNuFlux.h:108
#define pNOTICE
Definition: Messenger.h:61
void GJPARCNuFlux::Clear ( Option_t *  opt)
virtual

reset state variables based on opt

Implements genie::GFluxI.

Definition at line 806 of file GJPARCNuFlux.cxx.

References fEntriesThisCycle, fICycle, fIEntry, fNNeutrinos, fOffset, fSumWeight, and GenerateWeighted().

807 {
808 // If opt = "CycleHistory" then:
809 // Reset all counters and state variables from any previous cycles. This
810 // should be called if, for instance, before event generation this flux
811 // driver had been used for pre-calculating interaction probabilities. Does
812 // not reset initial state variables such as flux particle types, POTs etc...
813 //
814  if(std::strcmp(opt, "CycleHistory") == 0){
815  // Reset so that generate de-weighted events
816  this->GenerateWeighted(false);
817  // Reset cycle counters
818  fICycle = 0;
819  fSumWeight = 0;
820  fNNeutrinos = 0;
821  fIEntry = fOffset;
822  fEntriesThisCycle = 0;
823  }
824 }
void GenerateWeighted(bool gen_weighted=true)
set whether to generate weighted or unweighted neutrinos
long int fOffset
start looping at entry fOffset
Definition: GJPARCNuFlux.h:128
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GJPARCNuFlux.h:136
long int fIEntry
current flux ntuple entry
Definition: GJPARCNuFlux.h:126
int fICycle
current cycle
Definition: GJPARCNuFlux.h:134
double fSumWeight
sum of weights for neutrinos thrown so far
Definition: GJPARCNuFlux.h:135
long int fEntriesThisCycle
keep track of number of entries used so far for this cycle
Definition: GJPARCNuFlux.h:127
void genie::flux::GJPARCNuFlux::DisableOffset ( void  )
inline

switch off random offset, must be called before LoadBeamSimData to have any effect

Definition at line 83 of file GJPARCNuFlux.h.

References fUseRandomOffset.

Referenced by main().

int GJPARCNuFlux::DLocName2Id ( string  name)
private

Definition at line 922 of file GJPARCNuFlux.cxx.

Referenced by LoadBeamSimData().

923 {
924 // detector location: name -> int id
925 // sk -> -1
926 // nd1 -> +1
927 // nd2 -> +2
928 // ...
929 
930  if(name == "sk" ) return -1;
931 
932  TString temp;
933  for (int i=1; i<51; i++) {
934  temp.Form("nd%d",i);
935  if(name == temp.Data()) return i;
936  }
937 
938  return 0;
939 }
const char * name
bool genie::flux::GJPARCNuFlux::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 66 of file GJPARCNuFlux.h.

References fEntriesThisCycle, fICycle, fNCycles, and fNEntries.

Referenced by GenerateNext(), and main().

66  { return fEntriesThisCycle >= fNEntries
67  && fICycle == fNCycles && fNCycles > 0; }
int fNCycles
how many times to cycle through the flux ntuple
Definition: GJPARCNuFlux.h:133
int fICycle
current cycle
Definition: GJPARCNuFlux.h:134
long int fNEntries
number of flux ntuple entries
Definition: GJPARCNuFlux.h:125
long int fEntriesThisCycle
keep track of number of entries used so far for this cycle
Definition: GJPARCNuFlux.h:127
const PDGCodeList& genie::flux::GJPARCNuFlux::FluxParticles ( void  )
inlinevirtual

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

Implements genie::GFluxI.

Definition at line 59 of file GJPARCNuFlux.h.

References fPdgCList.

59 { return *fPdgCList; }
PDGCodeList * fPdgCList
list of neutrino pdg-codes
Definition: GJPARCNuFlux.h:108
bool GJPARCNuFlux::GenerateNext ( void  )
virtual

generate the next flux neutrino (return false in err)

Implements genie::GFluxI.

Definition at line 55 of file GJPARCNuFlux.cxx.

References End(), fGenerateWeighted, fICycle, fNCycles, GenerateNext_weighted(), Index(), genie::RandomGen::Instance(), genie::controls::kASmallNum, LOG, pERROR, pNOTICE, genie::RandomGen::RndFlux(), and Weight().

56 {
57 // Get next (unweighted) flux ntuple entry on the specified detector location
58 //
60  while(1) {
61  // Check for end of flux ntuple
62  bool end = this->End();
63  if(end) return false;
64 
65  // Get next weighted flux ntuple entry
66  bool nextok = this->GenerateNext_weighted();
67  if(!nextok) continue;
68 
69  if(fNCycles==0) {
70  LOG("Flux", pNOTICE)
71  << "Got flux entry: " << this->Index()
72  << " - Cycle: "<< fICycle << "/ infinite";
73  } else {
74  LOG("Flux", pNOTICE)
75  << "Got flux entry: "<< this->Index()
76  << " - Cycle: "<< fICycle << "/"<< fNCycles;
77  }
78 
79  // If de-weighting get fractional weight & decide whether to accept curr flux neutrino
80  double f = 1.0;
81  if(fGenerateWeighted == false) f = this->Weight();
82  LOG("Flux", pNOTICE)
83  << "Curr flux neutrino fractional weight = " << f;
84  if(f > (1.+controls::kASmallNum)) {
85  LOG("Flux", pERROR)
86  << "** Fractional weight = " << f << " > 1 !!";
87  }
88  double r = (f < 1.) ? rnd->RndFlux().Rndm() : 0;
89  bool accept = (r<f);
90  if(accept) {
91  return true;
92  }
93 
94  LOG("Flux", pNOTICE)
95  << "** Rejecting current flux neutrino based on the flux weight only";
96  }
97  return false;
98 }
bool End(void)
true if no more flux nu&#39;s can be thrown (eg reaching end of beam sim ntuples)
Definition: GJPARCNuFlux.h:66
#define pERROR
Definition: Messenger.h:59
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
double Weight(void)
returns the flux neutrino weight (if any)
Definition: GJPARCNuFlux.h:63
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
int fNCycles
how many times to cycle through the flux ntuple
Definition: GJPARCNuFlux.h:133
static const double kASmallNum
Definition: Controls.h:40
bool fGenerateWeighted
generate weighted/deweighted flux neutrinos (default is false)
Definition: GJPARCNuFlux.h:139
int fICycle
current cycle
Definition: GJPARCNuFlux.h:134
long int Index(void)
returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current ent...
#define pNOTICE
Definition: Messenger.h:61
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:71
bool GJPARCNuFlux::GenerateNext_weighted ( void  )
private

Definition at line 100 of file GJPARCNuFlux.cxx.

References genie::units::cm, genie::flux::GJPARCNuFluxPassThroughInfo::Enu, genie::PDGCodeList::ExistsInPDGCodeList(), fDetLocId, fEntriesThisCycle, fgP4, fgPdgC, fgX4, fICycle, fIEntry, fIsNDLoc, fLoadedNeutrino, genie::flux::GJPARCNuFluxPassThroughInfo::fluxentry, genie::flux::GJPARCNuFluxPassThroughInfo::fluxfilename, fMaxEv, fMaxWeight, fNCycles, fNDetLocIdFound, fNEntries, fNNeutrinos, fNorm, fNuFluxChain, fNuFluxFile, fNuFluxSumChain, fNuFluxSumTree, fNuFluxTree, fNuFluxUsingTree, fOffset, fPassThroughInfo, fPdgCList, fSumWeight, fZ0, genie::flux::GJPARCNuFluxPassThroughInfo::idfd, Index(), genie::controls::kASmallNum, genie::kPdgAntiNuE, genie::kPdgAntiNuMu, genie::kPdgNuE, genie::kPdgNuMu, LOG, genie::units::m, genie::flux::GJPARCNuFluxPassThroughInfo::mode, genie::flux::GJPARCNuFluxPassThroughInfo::nnu, genie::flux::GJPARCNuFluxPassThroughInfo::norm, genie::utils::print::P4AsShortString(), pERROR, pFATAL, pINFO, pNOTICE, pWARN, ResetCurrent(), Weight(), genie::utils::print::X4AsString(), genie::flux::GJPARCNuFluxPassThroughInfo::xnu, and genie::flux::GJPARCNuFluxPassThroughInfo::ynu.

Referenced by GenerateNext().

101 {
102 // Get next (weighted) flux ntuple entry on the specified detector location
103 //
104 
105  // Reset previously generated neutrino code / 4-p / 4-x
106  this->ResetCurrent();
107 
108  // Check whether a jnubeam flux ntuple has been loaded
110  LOG("Flux", pFATAL)
111  << "The flux driver has not been properly configured";
112  //return false; // don't do this - creates an infinite loop!
113  exit(1);
114  }
115 
116  // Read next flux ntuple entry. Use fEntriesThisCycle to keep track of when
117  // in new cycle as fIEntry can now have an offset
119  // Exit if have not found neutrino at specified location for whole cycle
120  if(fNDetLocIdFound == 0){
121  LOG("Flux", pFATAL)
122  << "The input jnubeam flux ntuple contains no entries for detector id "
123  << fDetLocId << ". Terminating job!";
124  exit(1);
125  }
126  fNDetLocIdFound = 0; // reset the counter
127  fICycle++;
129  fEntriesThisCycle = 0;
130  // Run out of entries @ the current cycle.
131  // Check whether more (or infinite) number of cycles is requested
132  if(fICycle >= fNCycles && fNCycles != 0){
133  LOG("Flux", pWARN)
134  << "No more entries in input flux neutrino ntuple";
135  return false;
136  }
137  }
138 
139  // In addition to getting info to generate event the following also
140  // updates pass-through info (= info on the flux neutrino parent particle
141  // that may be stored at an extra branch of the output event tree -alongside
142  // with the generated event branch- for use further upstream in the t2k
143  // analysis chain -eg for beam reweighting etc-)
144  bool found_entry;
145  if (fNuFluxUsingTree)
146  found_entry = fNuFluxTree->GetEntry(fIEntry) > 0;
147  else
148  found_entry = fNuFluxChain->GetEntry(fIEntry) > 0;
149  assert(found_entry);
150  fLoadedNeutrino = true;
152  fIEntry = (fIEntry+1) % fNEntries;
153 
154  if (fNuFluxUsingTree) {
155  if(fNuFluxSumTree) fNuFluxSumTree->GetEntry(0); // get entry 0 as only 1 entry in tree
156  }
157  else {
158  // get entry corresponding to current tree number in the chain, as only 1 entry in each tree
159  if(fNuFluxSumChain) fNuFluxSumChain->GetEntry(fNuFluxChain->GetTreeNumber());
160  }
161 
162  // check for negative flux weights
164  LOG("Flux", pERROR) << "Negative flux weight! Will set weight to 0.0";
165  fPassThroughInfo->norm = 0.0;
166  }
167  // remember to update fNorm as no longer set to fNuFluxTree branch address
168  fNorm = (double) fPassThroughInfo->norm;
169 
170  // for 'near detector' flux ntuples make sure that the current entry
171  // corresponds to a flux neutrino at the specified detector location
172  if(fIsNDLoc /* nd */ &&
173  fDetLocId!=fPassThroughInfo->idfd /* doesn't match specified detector location*/) {
174 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
175  LOG("Flux", pNOTICE)
176  << "Current flux neutrino not at specified detector location";
177 #endif
178  return false;
179  }
180 
181 
182  //
183  // Handling neutrinos at specified detector location
184  //
185 
186  // count the number of times we have neutrinos at specified detector location
187  fNDetLocIdFound += 1;
188 
189 
190  // update the sum of weights & number of neutrinos
191  fSumWeight += this->Weight() * fMaxWeight; // Weight returns fNorm/fMaxWeight
192  fNNeutrinos++;
193 
194  // Convert the current parent paticle decay mode into a neutrino pdg code
195  // See:
196  // http://jnusrv01.kek.jp/internal/t2k/nubeam/flux/nemode.h
197  // mode description
198  // 11 numu from pi+
199  // 12 numu from K+
200  // 13 numu from mu-
201  // 21 numu_bar from pi-
202  // 22 numu_bar from K-
203  // 23 numu_bar from mu+
204  // 31 nue from K+ (Ke3)
205  // 32 nue from K0L(Ke3)
206  // 33 nue from Mu+
207  // 41 nue_bar from K- (Ke3)
208  // 42 nue_bar from K0L(Ke3)
209  // 43 nue_bar from Mu-
210  // Since JNuBeam flux version >= 10a the following modes also expected
211  // 14 numu from K+(3)
212  // 15 numu from K0(3)
213  // 24 numu_bar from K-(3)
214  // 25 numu_bar from K0(3)
215  // 34 nue from pi+
216  // 44 nue_bar from pi-
217  // In general expect more modes following the rule:
218  // 11->19 --> numu
219  // 21->29 --> numu_bar
220  // 31->39 --> nue
221  // 41->49 --> nuebar
222  // This is based on example given at:
223  // http://jnusrv01.kek.jp/internal/t2k/nubeam/flux/efill.kumac
224 
225  if(fPassThroughInfo->mode >= 11 && fPassThroughInfo->mode <= 19) fgPdgC = kPdgNuMu;
226  else if(fPassThroughInfo->mode >= 21 && fPassThroughInfo->mode <= 29) fgPdgC = kPdgAntiNuMu;
227  else if(fPassThroughInfo->mode >= 31 && fPassThroughInfo->mode <= 39) fgPdgC = kPdgNuE;
228  else if(fPassThroughInfo->mode >= 41 && fPassThroughInfo->mode <= 49) fgPdgC = kPdgAntiNuE;
229  else {
230  // If here then trying to process a neutrino from an unknown decay mode.
231  // Rather than just skipping this flux neutrino the job is aborted to avoid
232  // unphysical results.
233  LOG("Flux", pFATAL) << "Unexpected decay mode: "<< fPassThroughInfo->mode <<
234  " --> unable to infer neutrino pdg! Aborting job!";
235  exit(1);
236  }
237 
238  // Check neutrino pdg against declared list of neutrino species declared
239  // by the current instance of the JPARC neutrino flux driver.
240  // No undeclared neutrino species will be accepted at this point as GENIE
241  // has already been configured to handle the specified list.
242  // Make sure that the appropriate list of flux neutrino species was set at
243  // initialization via GJPARCNuFlux::SetFluxParticles(const PDGCodeList &)
244 
246 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
247  LOG("Flux", pNOTICE)
248  << "Current flux neutrino "
249  << "not at the list of neutrinos to be considered at this job.";
250 #endif
251  return false;
252  }
253 
254  // Check current neutrino energy against the maximum flux neutrino energy declared
255  // by the current instance of the JPARC neutrino flux driver.
256  // No flux neutrino exceeding that maximum energy will be accepted at this point as
257  // that maximum energy has already been used for normalizing the interaction probabilities.
258  // Make sure that the appropriate maximum flux neutrino energy was set at
259  // initialization via GJPARCNuFlux::SetMaxEnergy(double Ev)
260 
261  if(fPassThroughInfo->Enu > fMaxEv) {
262  LOG("Flux", pWARN)
263  << "Flux neutrino energy exceeds declared maximum neutrino energy";
264  LOG("Flux", pWARN)
265  << "Ev = " << fPassThroughInfo->Enu << "(> Ev{max} = " << fMaxEv << ")";
266  }
267 
268  // Set the current flux neutrino 4-momentum & 4-position
269 
270  double pxnu = fPassThroughInfo->Enu * fPassThroughInfo->nnu[0];
271  double pynu = fPassThroughInfo->Enu * fPassThroughInfo->nnu[1];
272  double pznu = fPassThroughInfo->Enu * fPassThroughInfo->nnu[2];
273  double Enu = fPassThroughInfo->Enu;
274  fgP4.SetPxPyPzE (pxnu, pynu, pznu, Enu);
275 
276  if(fIsNDLoc) {
277  double cm2m = units::cm / units::m;
278  double xnu = cm2m * fPassThroughInfo->xnu;
279  double ynu = cm2m * fPassThroughInfo->ynu;
280  double znu = 0;
281 
282  // projected 4-position (from z=0) back to a configurable plane
283  // position (fZ0) upstream of the detector face.
284  xnu += (fZ0/fPassThroughInfo->nnu[2])*fPassThroughInfo->nnu[0];
285  ynu += (fZ0/fPassThroughInfo->nnu[2])*fPassThroughInfo->nnu[1];
286  znu = fZ0;
287 
288  fgX4.SetXYZT (xnu, ynu, znu, 0.);
289  } else {
290  fgX4.SetXYZT (0.,0.,0.,0.);
291 
292  }
293 
294 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
295  LOG("Flux", pINFO)
296  << "Generated neutrino: "
297  << "\n pdg-code: " << fgPdgC
298  << "\n p4: " << utils::print::P4AsShortString(&fgP4)
299  << "\n x4: " << utils::print::X4AsString(&fgX4);
300 #endif
301  // Update flux pass through info not set as branch addresses of flux ntuples
302  if(fNuFluxUsingTree)
303  fPassThroughInfo->fluxentry = this->Index();
304  else
305  fPassThroughInfo->fluxentry = fNuFluxChain->GetTree()->GetReadEntry();
306 
307  std::string filename;
308  if(fNuFluxUsingTree)
309  filename = fNuFluxFile->GetName();
310  else
311  filename = fNuFluxChain->GetFile()->GetName();
312  std::string::size_type start_pos = filename.rfind("/");
313  if (start_pos == std::string::npos) start_pos = 0; else ++start_pos;
314  std::string basename(filename,start_pos);
315  if(fNuFluxUsingTree)
316  fPassThroughInfo->fluxfilename = basename + ":" + fNuFluxTree->GetName();
317  else
318  fPassThroughInfo->fluxfilename = basename + ":" + fNuFluxChain->GetName();
319  return true;
320 }
static constexpr double cm
Definition: Units.h:68
TTree * fNuFluxTree
input flux ntuple
Definition: GJPARCNuFlux.h:115
int fgPdgC
running generated nu pdg-code
Definition: GJPARCNuFlux.h:110
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:45
const int kPdgNuE
Definition: PDGCodes.h:28
#define pERROR
Definition: Messenger.h:59
bool fLoadedNeutrino
set to true when GenerateNext_weighted has been called successfully
Definition: GJPARCNuFlux.h:141
TChain * fNuFluxSumChain
input summary ntuple for flux file. This tree is only present for later flux versions &gt; 10a ...
Definition: GJPARCNuFlux.h:118
bool fNuFluxUsingTree
are we using a TTree or a TChain to view the input flux file?
Definition: GJPARCNuFlux.h:119
const int kPdgAntiNuE
Definition: PDGCodes.h:29
#define pFATAL
Definition: Messenger.h:56
const int kPdgNuMu
Definition: PDGCodes.h:30
TChain * fNuFluxChain
input flux ntuple
Definition: GJPARCNuFlux.h:116
double Weight(void)
returns the flux neutrino weight (if any)
Definition: GJPARCNuFlux.h:63
bool ExistsInPDGCodeList(int pdg_code) const
long int fOffset
start looping at entry fOffset
Definition: GJPARCNuFlux.h:128
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GJPARCNuFlux.h:136
TTree * fNuFluxSumTree
input summary ntuple for flux file. This tree is only present for later flux versions &gt; 10a ...
Definition: GJPARCNuFlux.h:117
TLorentzVector fgP4
running generated nu 4-momentum
Definition: GJPARCNuFlux.h:111
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
int fDetLocId
input detector location id (fDetLoc -&gt; jnubeam idfd)
Definition: GJPARCNuFlux.h:121
int fNCycles
how many times to cycle through the flux ntuple
Definition: GJPARCNuFlux.h:133
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:62
int fNDetLocIdFound
per cycle keep track of the number of fDetLocId are found - if this is zero will exit job ...
Definition: GJPARCNuFlux.h:122
long int fIEntry
current flux ntuple entry
Definition: GJPARCNuFlux.h:126
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
bool fIsNDLoc
input location is a &#39;near&#39; detector location?
Definition: GJPARCNuFlux.h:124
#define pWARN
Definition: Messenger.h:60
GJPARCNuFluxPassThroughInfo * fPassThroughInfo
Definition: GJPARCNuFlux.h:143
TFile * fNuFluxFile
input flux file
Definition: GJPARCNuFlux.h:114
double fMaxWeight
max flux neutrino weight in input file for the specified detector location
Definition: GJPARCNuFlux.h:130
int fICycle
current cycle
Definition: GJPARCNuFlux.h:134
long int Index(void)
returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current ent...
double fZ0
configurable starting z position for each flux neutrino (in detector coord system) ...
Definition: GJPARCNuFlux.h:132
PDGCodeList * fPdgCList
list of neutrino pdg-codes
Definition: GJPARCNuFlux.h:108
long int fNEntries
number of flux ntuple entries
Definition: GJPARCNuFlux.h:125
TLorentzVector fgX4
running generated nu 4-position
Definition: GJPARCNuFlux.h:112
#define pNOTICE
Definition: Messenger.h:61
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:57
double fNorm
current flux ntuple normalisation
Definition: GJPARCNuFlux.h:129
static constexpr double m
Definition: Units.h:71
double fSumWeight
sum of weights for neutrinos thrown so far
Definition: GJPARCNuFlux.h:135
long int fEntriesThisCycle
keep track of number of entries used so far for this cycle
Definition: GJPARCNuFlux.h:127
double fMaxEv
maximum energy
Definition: GJPARCNuFlux.h:107
void GJPARCNuFlux::GenerateWeighted ( bool  gen_weighted = true)
virtual

set whether to generate weighted or unweighted neutrinos

Implements genie::GFluxI.

Definition at line 781 of file GJPARCNuFlux.cxx.

References fGenerateWeighted.

Referenced by Clear().

782 {
783 // Ignore the flux weights when getting next flux neutrino. Always set to
784 // true for unweighted event generation but need option to switch off when
785 // using pre-calculated interaction probabilities for each flux ray to speed
786 // up event generation.
787  fGenerateWeighted = gen_weighted;
788 }
bool fGenerateWeighted
generate weighted/deweighted flux neutrinos (default is false)
Definition: GJPARCNuFlux.h:139
long int GJPARCNuFlux::Index ( void  )
virtual

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

Implements genie::GFluxI.

Definition at line 363 of file GJPARCNuFlux.cxx.

References fIEntry, fLoadedNeutrino, and fNEntries.

Referenced by GenerateNext(), and GenerateNext_weighted().

364 {
365 // Return the current flux entry index. If GenerateNext has not yet been
366 // called then return -1.
367 //
368  if(fLoadedNeutrino){
369  // subtract 1 as fIEntry was incremented since call to TTree::GetEntry
370  // and deal with special case where fIEntry-1 is last entry in cycle
371  return fIEntry == 0 ? fNEntries - 1 : fIEntry-1;
372  }
373  // return -1 if no neutrino loaded since last call to this->ResetCurrent()
374  return -1;
375 }
bool fLoadedNeutrino
set to true when GenerateNext_weighted has been called successfully
Definition: GJPARCNuFlux.h:141
long int fIEntry
current flux ntuple entry
Definition: GJPARCNuFlux.h:126
long int fNEntries
number of flux ntuple entries
Definition: GJPARCNuFlux.h:125
void GJPARCNuFlux::Initialize ( void  )
private

Definition at line 826 of file GJPARCNuFlux.cxx.

References fDetLoc, fDetLocId, fEntriesThisCycle, fFilePOT, fGenerateWeighted, fICycle, fIEntry, fIsFDLoc, fIsNDLoc, fLoadedNeutrino, fMaxEv, fMaxWeight, fNCycles, fNDetLocIdFound, fNEntries, fNNeutrinos, fNNeutrinosTot1c, fNorm, fNuFluxChain, fNuFluxFile, fNuFluxSumChain, fNuFluxSumTree, fNuFluxTree, fNuFluxUsingTree, fOffset, fPassThroughInfo, fPdgCList, fSumWeight, fSumWeightTot1c, fUseRandomOffset, fZ0, LOG, pNOTICE, ResetCurrent(), and SetDefaults().

827 {
828  LOG("Flux", pNOTICE) << "Initializing GJPARCNuFlux driver";
829 
830  fMaxEv = 0;
831  fPdgCList = new PDGCodeList;
833 
834  fNuFluxFile = 0;
835  fNuFluxTree = 0;
836  fNuFluxChain = 0;
837  fNuFluxSumTree = 0;
838  fNuFluxSumChain = 0;
839  fNuFluxUsingTree = true;
840  fDetLoc = "";
841  fDetLocId = 0;
842  fNDetLocIdFound = 0;
843  fIsFDLoc = false;
844  fIsNDLoc = false;
845 
846  fNEntries = 0;
847  fIEntry = 0;
849  fOffset = 0;
850  fNorm = 0.;
851  fMaxWeight = 0;
852  fFilePOT = 0;
853  fZ0 = 0;
854  fNCycles = 0;
855  fICycle = 0;
856  fSumWeight = 0;
857  fNNeutrinos = 0;
858  fSumWeightTot1c = 0;
859  fNNeutrinosTot1c = 0;
860  fGenerateWeighted= false;
861  fUseRandomOffset = true;
862  fLoadedNeutrino = false;
863 
864  this->SetDefaults();
865  this->ResetCurrent();
866 }
TTree * fNuFluxTree
input flux ntuple
Definition: GJPARCNuFlux.h:115
double fSumWeightTot1c
total sum of weights for neutrinos to be thrown / cycle
Definition: GJPARCNuFlux.h:137
bool fLoadedNeutrino
set to true when GenerateNext_weighted has been called successfully
Definition: GJPARCNuFlux.h:141
string fDetLoc
input detector location (&#39;sk&#39;,&#39;nd1&#39;,&#39;nd2&#39;,...)
Definition: GJPARCNuFlux.h:120
TChain * fNuFluxSumChain
input summary ntuple for flux file. This tree is only present for later flux versions &gt; 10a ...
Definition: GJPARCNuFlux.h:118
bool fNuFluxUsingTree
are we using a TTree or a TChain to view the input flux file?
Definition: GJPARCNuFlux.h:119
double fFilePOT
file POT normalization, typically 1E+21
Definition: GJPARCNuFlux.h:131
TChain * fNuFluxChain
input flux ntuple
Definition: GJPARCNuFlux.h:116
long int fOffset
start looping at entry fOffset
Definition: GJPARCNuFlux.h:128
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GJPARCNuFlux.h:136
TTree * fNuFluxSumTree
input summary ntuple for flux file. This tree is only present for later flux versions &gt; 10a ...
Definition: GJPARCNuFlux.h:117
A list of PDG codes.
Definition: PDGCodeList.h:32
long int fNNeutrinosTot1c
total number of flux neutrinos to be thrown / cycle
Definition: GJPARCNuFlux.h:138
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
int fDetLocId
input detector location id (fDetLoc -&gt; jnubeam idfd)
Definition: GJPARCNuFlux.h:121
int fNCycles
how many times to cycle through the flux ntuple
Definition: GJPARCNuFlux.h:133
int fNDetLocIdFound
per cycle keep track of the number of fDetLocId are found - if this is zero will exit job ...
Definition: GJPARCNuFlux.h:122
long int fIEntry
current flux ntuple entry
Definition: GJPARCNuFlux.h:126
bool fIsNDLoc
input location is a &#39;near&#39; detector location?
Definition: GJPARCNuFlux.h:124
GJPARCNuFluxPassThroughInfo * fPassThroughInfo
Definition: GJPARCNuFlux.h:143
bool fGenerateWeighted
generate weighted/deweighted flux neutrinos (default is false)
Definition: GJPARCNuFlux.h:139
TFile * fNuFluxFile
input flux file
Definition: GJPARCNuFlux.h:114
double fMaxWeight
max flux neutrino weight in input file for the specified detector location
Definition: GJPARCNuFlux.h:130
int fICycle
current cycle
Definition: GJPARCNuFlux.h:134
bool fIsFDLoc
input location is a &#39;far&#39; detector location?
Definition: GJPARCNuFlux.h:123
double fZ0
configurable starting z position for each flux neutrino (in detector coord system) ...
Definition: GJPARCNuFlux.h:132
PDGCodeList * fPdgCList
list of neutrino pdg-codes
Definition: GJPARCNuFlux.h:108
long int fNEntries
number of flux ntuple entries
Definition: GJPARCNuFlux.h:125
bool fUseRandomOffset
whether set random starting point when looping over flux ntuples
Definition: GJPARCNuFlux.h:140
#define pNOTICE
Definition: Messenger.h:61
double fNorm
current flux ntuple normalisation
Definition: GJPARCNuFlux.h:129
double fSumWeight
sum of weights for neutrinos thrown so far
Definition: GJPARCNuFlux.h:135
long int fEntriesThisCycle
keep track of number of entries used so far for this cycle
Definition: GJPARCNuFlux.h:127
double fMaxEv
maximum energy
Definition: GJPARCNuFlux.h:107
bool GJPARCNuFlux::LoadBeamSimData ( string  filename,
string  det_loc 
)

load a jnubeam root flux ntuple

Definition at line 377 of file GJPARCNuFlux.cxx.

References genie::flux::GJPARCNuFluxPassThroughInfo::alpha, genie::flux::GJPARCNuFluxPassThroughInfo::anorm, genie::flux::GJPARCNuFluxPassThroughInfo::bpos, genie::flux::GJPARCNuFluxPassThroughInfo::brms, genie::flux::GJPARCNuFluxPassThroughInfo::btilt, genie::flux::GJPARCNuFluxPassThroughInfo::cospi0bm, genie::flux::GJPARCNuFluxPassThroughInfo::cospibm, DLocName2Id(), genie::flux::GJPARCNuFluxPassThroughInfo::emit, genie::flux::GJPARCNuFluxPassThroughInfo::Enu, genie::flux::GJPARCNuFluxPassThroughInfo::Enusk, fDetLoc, fDetLocId, fICycle, fIsFDLoc, fIsNDLoc, fMaxWeight, fNDetLocIdFound, fNEntries, fNNeutrinosTot1c, fNorm, fNuFluxChain, fNuFluxFile, fNuFluxSumChain, fNuFluxSumTree, fNuFluxTree, fNuFluxUsingTree, fPassThroughInfo, fSumWeightTot1c, fUseRandomOffset, genie::flux::GJPARCNuFluxPassThroughInfo::gamom0, genie::flux::GJPARCNuFluxPassThroughInfo::gcosbm, genie::flux::GJPARCNuFluxPassThroughInfo::gdistal, genie::flux::GJPARCNuFluxPassThroughInfo::gdistc, genie::flux::GJPARCNuFluxPassThroughInfo::gdistfe, genie::flux::GJPARCNuFluxPassThroughInfo::gdistti, genie::flux::GJPARCNuFluxPassThroughInfo::gipart, genie::flux::GJPARCNuFluxPassThroughInfo::gmat, genie::flux::GJPARCNuFluxPassThroughInfo::gmec, genie::flux::GJPARCNuFluxPassThroughInfo::gpid, genie::flux::GJPARCNuFluxPassThroughInfo::gpos0, genie::flux::GJPARCNuFluxPassThroughInfo::gpx, genie::flux::GJPARCNuFluxPassThroughInfo::gpy, genie::flux::GJPARCNuFluxPassThroughInfo::gpz, genie::flux::GJPARCNuFluxPassThroughInfo::gvec0, genie::flux::GJPARCNuFluxPassThroughInfo::gvx, genie::flux::GJPARCNuFluxPassThroughInfo::gvy, genie::flux::GJPARCNuFluxPassThroughInfo::gvz, genie::flux::GJPARCNuFluxPassThroughInfo::hcur, genie::flux::GJPARCNuFluxPassThroughInfo::idfd, genie::controls::kASmallNum, LOG, genie::flux::GJPARCNuFluxPassThroughInfo::mode, genie::flux::GJPARCNuFluxPassThroughInfo::ng, genie::flux::GJPARCNuFluxPassThroughInfo::nnu, genie::flux::GJPARCNuFluxPassThroughInfo::norm, genie::flux::GJPARCNuFluxPassThroughInfo::normsk, genie::flux::GJPARCNuFluxPassThroughInfo::npi, genie::flux::GJPARCNuFluxPassThroughInfo::npi0, genie::flux::GJPARCNuFluxPassThroughInfo::ntrig, genie::flux::GJPARCNuFluxPassThroughInfo::nvtx0, pDEBUG, pERROR, pFATAL, pINFO, genie::flux::GJPARCNuFluxPassThroughInfo::pint, pNOTICE, genie::flux::GJPARCNuFluxPassThroughInfo::ppi, genie::flux::GJPARCNuFluxPassThroughInfo::ppi0, genie::flux::GJPARCNuFluxPassThroughInfo::ppid, genie::flux::GJPARCNuFluxPassThroughInfo::rand, RandomOffset(), genie::flux::GJPARCNuFluxPassThroughInfo::rnu, genie::flux::GJPARCNuFluxPassThroughInfo::rseed, genie::utils::str::Split(), genie::flux::GJPARCNuFluxPassThroughInfo::tuneid, genie::flux::GJPARCNuFluxPassThroughInfo::version, genie::flux::GJPARCNuFluxPassThroughInfo::xnu, genie::flux::GJPARCNuFluxPassThroughInfo::xpi, genie::flux::GJPARCNuFluxPassThroughInfo::xpi0, and genie::flux::GJPARCNuFluxPassThroughInfo::ynu.

Referenced by main().

378 {
379 // Loads in a jnubeam beam simulation root file (converted from hbook format)
380 // into the GJPARCNuFlux driver.
381 // The detector location can be any of:
382 // "sk","nd1" (<-2km),"nd5" (<-nd280),...,"nd10"
383 
384  LOG("Flux", pNOTICE)
385  << "Loading jnubeam flux tree from ROOT file: " << filename;
386  LOG("Flux", pNOTICE)
387  << "Detector location: " << detector_location;
388 
389  // Check to see if its a single flux file (/dir/root_filename.0.root)
390  // or a sequence of files to be tchained (e.g. /dir/root_filename@0@100)
391  fNuFluxUsingTree = true;
392  if (filename.find('@') != string::npos) { fNuFluxUsingTree = false; }
393 
394  vector<string> filenamev = utils::str::Split(filename,"@");
395  string fileroot = "";
396  int firstfile = -1, lastfile = -1;
397 
398  if (!fNuFluxUsingTree) {
399  if (filenamev.size() != 3) {
400  LOG("Flux", pFATAL)
401  << "Flux filename should be specfied as either:\n"
402  << "\t For a single input file: /dir/root_filename.#.root\n"
403  << "\t For multiple input files: /dir/root_filename@#@#";
404  exit(1);
405  }
406  fileroot = filenamev[0];
407  firstfile = atoi(filenamev[1].c_str());
408  lastfile = atoi(filenamev[2].c_str());
409  LOG("Flux", pNOTICE)
410  << "Chaining beam simulation output files with stem: " << fileroot
411  << " and run numbers in the range: [" << firstfile << ", " << firstfile << "]";
412  }
413 
414  if (fNuFluxUsingTree) {
415  bool is_accessible = ! (gSystem->AccessPathName( filename.c_str() ));
416  if (!is_accessible) {
417  LOG("Flux", pFATAL)
418  << "The input jnubeam flux file doesn't exist! Initialization failed!";
419  exit(1);
420  }
421  }
422  else {
423  bool please_exit = false;
424  for (int i = firstfile; i < lastfile+1; i++) {
425  bool is_accessible = ! (gSystem->AccessPathName( Form("%s.%i.root",fileroot.c_str(),i) ));
426  if (!is_accessible) {
427  LOG("Flux", pFATAL)
428  << "The input jnubeam flux file " << Form("%s.%i.root",fileroot.c_str(),i)
429  << "doesn't exist! Initialization failed!";
430  please_exit = true;
431  }
432  }
433  if(please_exit)
434  exit(1);
435  }
436 
437  fDetLoc = detector_location;
438  fDetLocId = this->DLocName2Id(fDetLoc);
439 
440  if(fDetLocId == 0) {
441  LOG("Flux", pERROR)
442  << " ** Unknown input detector location: " << fDetLoc;
443  return false;
444  }
445 
446  fIsFDLoc = (fDetLocId==-1);
447  fIsNDLoc = (fDetLocId>0);
448 
449  if (!fNuFluxUsingTree) {
450  string ntuple_name = (fIsNDLoc) ? "h3002" : "h2000";
451  fNuFluxChain = new TChain(ntuple_name.c_str());
452  int result = fNuFluxChain->Add( Form("%s.%i.root",fileroot.c_str(),firstfile), 0);
453  if (result != 1 && fIsNDLoc) {
454  LOG("Flux", pINFO)
455  << "Could not find tree h3002 in file " << Form("%s.%i.root",fileroot.c_str(),firstfile)
456  << " Trying tree h3001";
457  delete fNuFluxChain;
458  ntuple_name = "h3001";
459  fNuFluxChain = new TChain(ntuple_name.c_str());
460  result = fNuFluxChain->Add( Form("%s.%i.root",fileroot.c_str(),firstfile), 0);
461  }
462  if (result != 1) {
463  LOG("Flux", pERROR)
464  << "** Couldn't get flux tree: " << ntuple_name;
465  return false;
466  }
467 
468  for (int i = firstfile+1; i < lastfile+1; i++) {
469  result = fNuFluxChain->Add( Form("%s.%i.root",fileroot.c_str(),i), 0 );
470  if (result == 0)
471  LOG("Flux", pERROR)
472  << "** Couldn't get flux tree " << ntuple_name << " in file " << Form("%s.%i.root",fileroot.c_str(),i);
473  fNEntries = fNuFluxChain->GetEntries();
474  }
475  }
476 
477  else {
478  fNuFluxFile = new TFile(filename.c_str(), "read");
479  if(fNuFluxFile) {
480  // nd treename can be h3002 or h3001 depending on fluxfile version
481  string ntuple_name = (fIsNDLoc) ? "h3002" : "h2000";
482  fNuFluxTree = (TTree*) fNuFluxFile->Get(ntuple_name.c_str());
483  if(!fNuFluxTree && fIsNDLoc){
484  ntuple_name = "h3001";
485  fNuFluxTree = (TTree*) fNuFluxFile->Get(ntuple_name.c_str());
486  }
487  LOG("Flux", pINFO)
488  << "Getting flux tree: " << ntuple_name;
489  if(!fNuFluxTree) {
490  LOG("Flux", pERROR)
491  << "** Couldn't get flux tree: " << ntuple_name;
492  return false;
493  }
494  } else {
495  LOG("Flux", pERROR) << "** Couldn't open: " << filename;
496  return false;
497  }
498 
499  fNEntries = fNuFluxTree->GetEntries();
500  }
501 
502  LOG("Flux", pNOTICE)
503  << "Loaded flux tree contains " << fNEntries << " entries";
504 
505  LOG("Flux", pDEBUG)
506  << "Getting tree branches & setting leaf addresses";
507 
508  // try to get all the branches that we know about and only set address if
509  // they exist
510  bool missing_critical = false;
511  TBranch * fBr = 0;
513 
514  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("norm")) ) fBr->SetAddress(&info->norm);
515  else if( (fBr = fNuFluxChain->GetBranch("norm")) ) fNuFluxChain->SetBranchAddress("norm",&info->norm);
516  else {
517  LOG("Flux", pFATAL) <<"Cannot find flux branch: norm";
518  missing_critical = true;
519  }
520  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("Enu")) ) fBr->SetAddress(&info->Enu);
521  else if( (fBr = fNuFluxChain->GetBranch("Enu")) ) fNuFluxChain->SetBranchAddress("Enu",&info->Enu);
522  else {
523  LOG("Flux", pFATAL) <<"Cannot find flux branch: Enu";
524  missing_critical = true;
525  }
526  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ppid")) ) fBr->SetAddress(&info->ppid);
527  else if( (fBr = fNuFluxChain->GetBranch("ppid")) ) fNuFluxChain->SetBranchAddress("ppid",&info->ppid);
528  else {
529  LOG("Flux", pFATAL) <<"Cannot find flux branch: ppid";
530  missing_critical = true;
531  }
532  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("mode")) ) fBr->SetAddress(&info->mode);
533  else if( (fBr = fNuFluxChain->GetBranch("mode")) ) fNuFluxChain->SetBranchAddress("mode",&info->mode);
534  else {
535  LOG("Flux", pFATAL) <<"Cannot find flux branch: mode";
536  missing_critical = true;
537  }
538  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("rnu")) ) fBr->SetAddress(&info->rnu);
539  else if( (fBr = fNuFluxChain->GetBranch("rnu")) ) fNuFluxChain->SetBranchAddress("rnu",&info->rnu);
540  else if(fIsNDLoc){ // Only required for ND location
541  LOG("Flux", pFATAL) <<"Cannot find flux branch: rnu";
542  missing_critical = true;
543  }
544  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("xnu")) ) fBr->SetAddress(&info->xnu);
545  else if( (fBr = fNuFluxChain->GetBranch("xnu")) ) fNuFluxChain->SetBranchAddress("xnu",&info->xnu);
546  else if(fIsNDLoc){ // Only required for ND location
547  LOG("Flux", pFATAL) <<"Cannot find flux branch: xnu";
548  missing_critical = true;
549  }
550  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ynu")) ) fBr->SetAddress(&info->ynu);
551  else if( (fBr = fNuFluxChain->GetBranch("ynu")) ) fNuFluxChain->SetBranchAddress("ynu",&info->ynu);
552  else if(fIsNDLoc) { // Only required for ND location
553  LOG("Flux", pFATAL) <<"Cannot find flux branch: ynu";
554  missing_critical = true;
555  }
556  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("nnu")) ) fBr->SetAddress(info->nnu);
557  else if( (fBr = fNuFluxChain->GetBranch("nnu")) ) fNuFluxChain->SetBranchAddress("nnu",info->nnu);
558  else if(fIsNDLoc){ // Only required for ND location
559  LOG("Flux", pFATAL) <<"Cannot find flux branch: nnu";
560  missing_critical = true;
561  }
562  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("idfd")) ) fBr->SetAddress(&info->idfd);
563  else if( (fBr = fNuFluxChain->GetBranch("idfd")) ) fNuFluxChain->SetBranchAddress("idfd",&info->idfd);
564  else if(fIsNDLoc){ // Only required for ND location
565  LOG("Flux", pFATAL) <<"Cannot find flux branch: idfd";
566  missing_critical = true;
567  }
568  // check that have found essential branches
569  if(missing_critical){
570  LOG("Flux", pFATAL)
571  << "Unable to find critical information in the flux ntuple! Initialization failed!";
572  exit(1);
573  }
574  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ppi")) ) fBr->SetAddress(&info->ppi);
575  else if( (fBr = fNuFluxChain->GetBranch("ppi")) ) fNuFluxChain->SetBranchAddress("ppi",&info->ppi);
576  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("xpi")) ) fBr->SetAddress(info->xpi);
577  else if( (fBr = fNuFluxChain->GetBranch("xpi")) ) fNuFluxChain->SetBranchAddress("xpi",info->xpi);
578  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("npi")) ) fBr->SetAddress(info->npi);
579  else if( (fBr = fNuFluxChain->GetBranch("npi")) ) fNuFluxChain->SetBranchAddress("npi",info->npi);
580  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ppi0")) ) fBr->SetAddress(&info->ppi0);
581  else if( (fBr = fNuFluxChain->GetBranch("ppi0")) ) fNuFluxChain->SetBranchAddress("ppi0",&info->ppi0);
582  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("xpi0")) ) fBr->SetAddress(info->xpi0);
583  else if( (fBr = fNuFluxChain->GetBranch("xpi0")) ) fNuFluxChain->SetBranchAddress("xpi0",info->xpi0);
584  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("npi0")) ) fBr->SetAddress(info->npi0);
585  else if( (fBr = fNuFluxChain->GetBranch("npi0")) ) fNuFluxChain->SetBranchAddress("npi0",info->npi0);
586  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("nvtx0")) ) fBr->SetAddress(&info->nvtx0);
587  else if( (fBr = fNuFluxChain->GetBranch("nvtx0")) ) fNuFluxChain->SetBranchAddress("nvtx0",&info->nvtx0);
588  // Following branches only present since flux version 10a
589  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("cospibm")) ) fBr->SetAddress(&info->cospibm);
590  else if( (fBr = fNuFluxChain->GetBranch("cospibm")) ) fNuFluxChain->SetBranchAddress("cospibm",&info->cospibm);
591  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("cospi0bm"))) fBr->SetAddress(&info->cospi0bm);
592  else if( (fBr = fNuFluxChain->GetBranch("cospi0bm"))) fNuFluxChain->SetBranchAddress("cospi0bm",&info->cospi0bm);
593  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gamom0")) ) fBr->SetAddress(&info->gamom0);
594  else if( (fBr = fNuFluxChain->GetBranch("gamom0")) ) fNuFluxChain->SetBranchAddress("gamom0",&info->gamom0);
595  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gipart")) ) fBr->SetAddress(&info->gipart);
596  else if( (fBr = fNuFluxChain->GetBranch("gipart")) ) fNuFluxChain->SetBranchAddress("gipart",&info->gipart);
597  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gvec0")) ) fBr->SetAddress(info->gvec0);
598  else if( (fBr = fNuFluxChain->GetBranch("gvec0")) ) fNuFluxChain->SetBranchAddress("gvec0",info->gvec0);
599  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpos0")) ) fBr->SetAddress(info->gpos0);
600  else if( (fBr = fNuFluxChain->GetBranch("gpos0")) ) fNuFluxChain->SetBranchAddress("gpos0",info->gpos0);
601  // Following branches only present since flux vesion 10d
602  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ng")) ) fBr->SetAddress(&info->ng);
603  else if( (fBr = fNuFluxChain->GetBranch("ng")) ) fNuFluxChain->SetBranchAddress("ng",&info->ng);
604  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpid")) ) fBr->SetAddress(info->gpid);
605  else if( (fBr = fNuFluxChain->GetBranch("gpid")) ) fNuFluxChain->SetBranchAddress("gpid",info->gpid);
606  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gmec")) ) fBr->SetAddress(info->gmec);
607  else if( (fBr = fNuFluxChain->GetBranch("gmec")) ) fNuFluxChain->SetBranchAddress("gmec",info->gmec);
608  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gvx")) ) fBr->SetAddress(info->gvx);
609  else if( (fBr = fNuFluxChain->GetBranch("gvx")) ) fNuFluxChain->SetBranchAddress("gvx",info->gvx);
610  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gvy")) ) fBr->SetAddress(info->gvy);
611  else if( (fBr = fNuFluxChain->GetBranch("gvy")) ) fNuFluxChain->SetBranchAddress("gvy",info->gvy);
612  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gvz")) ) fBr->SetAddress(info->gvz);
613  else if( (fBr = fNuFluxChain->GetBranch("gvz")) ) fNuFluxChain->SetBranchAddress("gvz",info->gvz);
614  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpx")) ) fBr->SetAddress(info->gpx);
615  else if( (fBr = fNuFluxChain->GetBranch("gpx")) ) fNuFluxChain->SetBranchAddress("gpx",info->gpx);
616  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpy")) ) fBr->SetAddress(info->gpy);
617  else if( (fBr = fNuFluxChain->GetBranch("gpy")) ) fNuFluxChain->SetBranchAddress("gpy",info->gpy);
618  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpz")) ) fBr->SetAddress(info->gpz);
619  else if( (fBr = fNuFluxChain->GetBranch("gpz")) ) fNuFluxChain->SetBranchAddress("gpz",info->gpz);
620  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gmat")) ) fBr->SetAddress(info->gmat);
621  else if( (fBr = fNuFluxChain->GetBranch("gmat")) ) fNuFluxChain->SetBranchAddress("gmat",info->gmat);
622  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gdistc")) ) fBr->SetAddress(info->gdistc);
623  else if( (fBr = fNuFluxChain->GetBranch("gdistc")) ) fNuFluxChain->SetBranchAddress("gdistc",info->gdistc);
624  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gdistal")) ) fBr->SetAddress(&info->gdistal);
625  else if( (fBr = fNuFluxChain->GetBranch("gdistal")) ) fNuFluxChain->SetBranchAddress("gdistal",&info->gdistal);
626  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gdistti")) ) fBr->SetAddress(&info->gdistti);
627  else if( (fBr = fNuFluxChain->GetBranch("gdistti")) ) fNuFluxChain->SetBranchAddress("gdistti",&info->gdistti);
628  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gdistfe")) ) fBr->SetAddress(&info->gdistfe);
629  else if( (fBr = fNuFluxChain->GetBranch("gdistfe")) ) fNuFluxChain->SetBranchAddress("gdistfe",&info->gdistfe);
630  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gcosbm")) ) fBr->SetAddress(info->gcosbm);
631  else if( (fBr = fNuFluxChain->GetBranch("gcosbm")) ) fNuFluxChain->SetBranchAddress("gcosbm",info->gcosbm);
632  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("Enusk")) ) fBr->SetAddress(&info->Enusk);
633  else if( (fBr = fNuFluxChain->GetBranch("Enusk")) ) fNuFluxChain->SetBranchAddress("Enusk",&info->Enusk);
634  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("normsk")) ) fBr->SetAddress(&info->normsk);
635  else if( (fBr = fNuFluxChain->GetBranch("normsk")) ) fNuFluxChain->SetBranchAddress("normsk",&info->normsk);
636  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("anorm")) ) fBr->SetAddress(&info->anorm);
637  else if( (fBr = fNuFluxChain->GetBranch("anorm")) ) fNuFluxChain->SetBranchAddress("anorm",&info->anorm);
638 
639  // Look for the flux file summary info tree (only expected for > 10a flux versions)
640  if( fNuFluxUsingTree && (fNuFluxSumTree = (TTree*) fNuFluxFile->Get("h1000")) ){
641  if( (fBr = fNuFluxSumTree->GetBranch("version"))) fBr->SetAddress(&info->version);
642  if( (fBr = fNuFluxSumTree->GetBranch("ntrig")) ) fBr->SetAddress(&info->ntrig);
643  if( (fBr = fNuFluxSumTree->GetBranch("tuneid")) ) fBr->SetAddress(&info->tuneid);
644  if( (fBr = fNuFluxSumTree->GetBranch("pint")) ) fBr->SetAddress(&info->pint);
645  if( (fBr = fNuFluxSumTree->GetBranch("bpos")) ) fBr->SetAddress(info->bpos);
646  if( (fBr = fNuFluxSumTree->GetBranch("btilt")) ) fBr->SetAddress(info->btilt);
647  if( (fBr = fNuFluxSumTree->GetBranch("brms")) ) fBr->SetAddress(info->brms);
648  if( (fBr = fNuFluxSumTree->GetBranch("emit")) ) fBr->SetAddress(info->emit);
649  if( (fBr = fNuFluxSumTree->GetBranch("alpha")) ) fBr->SetAddress(info->alpha);
650  if( (fBr = fNuFluxSumTree->GetBranch("hcur")) ) fBr->SetAddress(info->hcur);
651  if( (fBr = fNuFluxSumTree->GetBranch("rand")) ) fBr->SetAddress(&info->rand);
652  if( (fBr = fNuFluxSumTree->GetBranch("rseed")) ) fBr->SetAddress(info->rseed);
653  }
654 
655  // Look for the flux file summary info tree (only expected for > 10a flux versions)
656  if ( !fNuFluxUsingTree ) {
657  fNuFluxSumChain = new TChain("h1000");
658  int result = fNuFluxSumChain->Add( Form("%s.%i.root",fileroot.c_str(),firstfile), 0 );
659  if (result==1) {
660  for (int i = firstfile+1; i < lastfile+1; i++) {
661  result = fNuFluxSumChain->Add( Form("%s.%i.root",fileroot.c_str(),i), 0 );
662  }
663  if( (fBr = fNuFluxSumChain->GetBranch("version"))) fNuFluxSumChain->SetBranchAddress("version",&info->version);
664  if( (fBr = fNuFluxSumChain->GetBranch("ntrig")) ) fNuFluxSumChain->SetBranchAddress("ntrig",&info->ntrig);
665  if( (fBr = fNuFluxSumChain->GetBranch("tuneid")) ) fNuFluxSumChain->SetBranchAddress("tuneid",&info->tuneid);
666  if( (fBr = fNuFluxSumChain->GetBranch("pint")) ) fNuFluxSumChain->SetBranchAddress("pint",&info->pint);
667  if( (fBr = fNuFluxSumChain->GetBranch("bpos")) ) fNuFluxSumChain->SetBranchAddress("bpos",info->bpos);
668  if( (fBr = fNuFluxSumChain->GetBranch("btilt")) ) fNuFluxSumChain->SetBranchAddress("btilt",info->btilt);
669  if( (fBr = fNuFluxSumChain->GetBranch("brms")) ) fNuFluxSumChain->SetBranchAddress("brms",info->brms);
670  if( (fBr = fNuFluxSumChain->GetBranch("emit")) ) fNuFluxSumChain->SetBranchAddress("emit",info->emit);
671  if( (fBr = fNuFluxSumChain->GetBranch("alpha")) ) fNuFluxSumChain->SetBranchAddress("alpha",info->alpha);
672  if( (fBr = fNuFluxSumChain->GetBranch("hcur")) ) fNuFluxSumChain->SetBranchAddress("hcur",info->hcur);
673  if( (fBr = fNuFluxSumChain->GetBranch("rand")) ) fNuFluxSumChain->SetBranchAddress("rand",&info->rand);
674  if( (fBr = fNuFluxSumChain->GetBranch("rseed")) ) fNuFluxSumChain->SetBranchAddress("rseed",info->rseed);
675  }
676  }
677 
678  // current ntuple cycle # (flux ntuples may be recycled)
679  fICycle = 1;
680 
681  // sum-up weights & number of neutrinos for the specified location
682  // over a complete cycle. Also record the maximum weight as previous
683  // method using TTree::GetV1() seg faulted for more than ~1.5E6 entries
684  fSumWeightTot1c = 0;
685  fNNeutrinosTot1c = 0;
686  fNDetLocIdFound = 0;
687  for(int ientry = 0; ientry < fNEntries; ientry++) {
688  if (fNuFluxUsingTree)
689  fNuFluxTree->GetEntry(ientry);
690  else
691  fNuFluxChain->GetEntry(ientry);
692  // check for negative flux weights
694  LOG("Flux", pERROR) << "Negative flux weight! Will set weight to 0.0";
695  fPassThroughInfo->norm = 0.0;
696  }
697  fNorm = (double) fPassThroughInfo->norm;
698  // update maximum weight
699  fMaxWeight = TMath::Max(fMaxWeight, (double) fPassThroughInfo->norm);
700  // compare detector location (see GenerateNext_weighted() for details)
701  if(fIsNDLoc && fDetLocId!=fPassThroughInfo->idfd) continue;
704  fNDetLocIdFound++;
705  }
706  // Exit if have not found neutrino at specified location for whole cycle
707  if(fNDetLocIdFound == 0){
708  LOG("Flux", pFATAL)
709  << "The input jnubeam flux ntuple contains no entries for detector id "
710  << fDetLocId << ". Terminating job!";
711  exit(1);
712  }
713  fNDetLocIdFound = 0; // reset the counter
714 
715  LOG("Flux", pNOTICE) << "Maximum flux weight = " << fMaxWeight;
716  if(fMaxWeight <=0 ) {
717  LOG("Flux", pFATAL) << "Non-positive maximum flux weight!";
718  exit(1);
719  }
720 
721  LOG("Flux", pINFO)
722  << "Totals / cycle: #neutrinos = " << fNNeutrinosTot1c
723  << ", Sum{Weights} = " << fSumWeightTot1c;
724 
725  if(fUseRandomOffset){
726  this->RandomOffset(); // Random start point when looping over ntuple
727  }
728 
729  return true;
730 }
TTree * fNuFluxTree
input flux ntuple
Definition: GJPARCNuFlux.h:115
double fSumWeightTot1c
total sum of weights for neutrinos to be thrown / cycle
Definition: GJPARCNuFlux.h:137
#define pERROR
Definition: Messenger.h:59
int DLocName2Id(string name)
string fDetLoc
input detector location (&#39;sk&#39;,&#39;nd1&#39;,&#39;nd2&#39;,...)
Definition: GJPARCNuFlux.h:120
TChain * fNuFluxSumChain
input summary ntuple for flux file. This tree is only present for later flux versions &gt; 10a ...
Definition: GJPARCNuFlux.h:118
bool fNuFluxUsingTree
are we using a TTree or a TChain to view the input flux file?
Definition: GJPARCNuFlux.h:119
#define pFATAL
Definition: Messenger.h:56
TChain * fNuFluxChain
input flux ntuple
Definition: GJPARCNuFlux.h:116
TTree * fNuFluxSumTree
input summary ntuple for flux file. This tree is only present for later flux versions &gt; 10a ...
Definition: GJPARCNuFlux.h:117
long int fNNeutrinosTot1c
total number of flux neutrinos to be thrown / cycle
Definition: GJPARCNuFlux.h:138
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
int fDetLocId
input detector location id (fDetLoc -&gt; jnubeam idfd)
Definition: GJPARCNuFlux.h:121
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:62
int fNDetLocIdFound
per cycle keep track of the number of fDetLocId are found - if this is zero will exit job ...
Definition: GJPARCNuFlux.h:122
bool fIsNDLoc
input location is a &#39;near&#39; detector location?
Definition: GJPARCNuFlux.h:124
GJPARCNuFluxPassThroughInfo * fPassThroughInfo
Definition: GJPARCNuFlux.h:143
TFile * fNuFluxFile
input flux file
Definition: GJPARCNuFlux.h:114
double fMaxWeight
max flux neutrino weight in input file for the specified detector location
Definition: GJPARCNuFlux.h:130
int fICycle
current cycle
Definition: GJPARCNuFlux.h:134
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
bool fIsFDLoc
input location is a &#39;far&#39; detector location?
Definition: GJPARCNuFlux.h:123
long int fNEntries
number of flux ntuple entries
Definition: GJPARCNuFlux.h:125
bool fUseRandomOffset
whether set random starting point when looping over flux ntuples
Definition: GJPARCNuFlux.h:140
#define pNOTICE
Definition: Messenger.h:61
void RandomOffset(void)
choose a random offset as starting entry in flux ntuple
double fNorm
current flux ntuple normalisation
Definition: GJPARCNuFlux.h:129
#define pDEBUG
Definition: Messenger.h:63
double genie::flux::GJPARCNuFlux::MaxEnergy ( void  )
inlinevirtual

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

Implements genie::GFluxI.

Definition at line 60 of file GJPARCNuFlux.h.

References fMaxEv.

60 { return fMaxEv; }
double fMaxEv
maximum energy
Definition: GJPARCNuFlux.h:107
const TLorentzVector& genie::flux::GJPARCNuFlux::Momentum ( void  )
inlinevirtual

returns the flux neutrino 4-momentum

Implements genie::GFluxI.

Definition at line 64 of file GJPARCNuFlux.h.

References fgP4.

64 { return fgP4; }
TLorentzVector fgP4
running generated nu 4-momentum
Definition: GJPARCNuFlux.h:111
long int genie::flux::GJPARCNuFlux::NFluxNeutrinos ( void  ) const
inline

number of flux neutrinos looped so far

Definition at line 88 of file GJPARCNuFlux.h.

References fNNeutrinos.

const GJPARCNuFluxPassThroughInfo& genie::flux::GJPARCNuFlux::PassThroughInfo ( void  )
inline

GJPARCNuFluxPassThroughInfo.

Definition at line 92 of file GJPARCNuFlux.h.

References fPassThroughInfo.

Referenced by main().

int genie::flux::GJPARCNuFlux::PdgCode ( void  )
inlinevirtual

returns the flux neutrino pdg code

Implements genie::GFluxI.

Definition at line 62 of file GJPARCNuFlux.h.

References fgPdgC.

62 { return fgPdgC; }
int fgPdgC
running generated nu pdg-code
Definition: GJPARCNuFlux.h:110
const TLorentzVector& genie::flux::GJPARCNuFlux::Position ( void  )
inlinevirtual

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

Implements genie::GFluxI.

Definition at line 65 of file GJPARCNuFlux.h.

References fgX4.

65 { return fgX4; }
TLorentzVector fgX4
running generated nu 4-position
Definition: GJPARCNuFlux.h:112
double GJPARCNuFlux::POT_1cycle ( void  )

flux POT per cycle

Definition at line 322 of file GJPARCNuFlux.cxx.

References fFilePOT, fMaxWeight, fNuFluxChain, fNuFluxTree, fNuFluxUsingTree, LOG, and pWARN.

Referenced by main(), and POT_curravg().

323 {
324 // Compute number of flux POTs / flux ntuple cycle
325 //
327  LOG("Flux", pWARN)
328  << "The flux driver has not been properly configured";
329  return 0;
330  }
331 
332 // double pot = fFilePOT * (fNNeutrinosTot1c/fSumWeightTot1c);
333 //
334 // Use the max weight instead, since flux neutrinos get de-weighted
335 // before thrown to the event generation driver
336 //
337  double pot = fFilePOT / fMaxWeight;
338  return pot;
339 }
TTree * fNuFluxTree
input flux ntuple
Definition: GJPARCNuFlux.h:115
bool fNuFluxUsingTree
are we using a TTree or a TChain to view the input flux file?
Definition: GJPARCNuFlux.h:119
double fFilePOT
file POT normalization, typically 1E+21
Definition: GJPARCNuFlux.h:131
TChain * fNuFluxChain
input flux ntuple
Definition: GJPARCNuFlux.h:116
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pWARN
Definition: Messenger.h:60
double fMaxWeight
max flux neutrino weight in input file for the specified detector location
Definition: GJPARCNuFlux.h:130
double GJPARCNuFlux::POT_curravg ( void  )

current average POT

Definition at line 341 of file GJPARCNuFlux.cxx.

References fNNeutrinos, fNNeutrinosTot1c, fNuFluxChain, fNuFluxTree, fNuFluxUsingTree, LOG, POT_1cycle(), and pWARN.

Referenced by main().

342 {
343 // Compute current number of flux POTs
344 // On complete cycles, that POT number should be exact.
345 // Within cycles that is only an average number
346 
348  LOG("Flux", pWARN)
349  << "The flux driver has not been properly configured";
350  return 0;
351  }
352 
353 // double pot = fNNeutrinos*fFilePOT/fSumWeightTot1c;
354 //
355 // See also comment at POT_1cycle()
356 //
357  double cnt = (double)fNNeutrinos;
358  double cnt1c = (double)fNNeutrinosTot1c;
359  double pot = (cnt/cnt1c) * this->POT_1cycle();
360  return pot;
361 }
TTree * fNuFluxTree
input flux ntuple
Definition: GJPARCNuFlux.h:115
double POT_1cycle(void)
flux POT per cycle
bool fNuFluxUsingTree
are we using a TTree or a TChain to view the input flux file?
Definition: GJPARCNuFlux.h:119
TChain * fNuFluxChain
input flux ntuple
Definition: GJPARCNuFlux.h:116
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GJPARCNuFlux.h:136
long int fNNeutrinosTot1c
total number of flux neutrinos to be thrown / cycle
Definition: GJPARCNuFlux.h:138
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pWARN
Definition: Messenger.h:60
void GJPARCNuFlux::RandomOffset ( void  )

choose a random offset as starting entry in flux ntuple

Definition at line 790 of file GJPARCNuFlux.cxx.

References fIEntry, fNEntries, fOffset, genie::RandomGen::Instance(), LOG, pERROR, and genie::RandomGen::RndFlux().

Referenced by LoadBeamSimData().

791 {
792 // Choose a random number between 0-->fNEntries to set as start point for
793 // looping over flux ntuple. May be necessary when looping over very large
794 // flux files as always starting fromthe same point may introduce biases
795 // (inversely proportional to number of cycles). This method resets the
796 // starting value for fIEntry so must be called before any call to GenerateNext
797 // is made.
798 //
799  double ran_frac = RandomGen::Instance()->RndFlux().Rndm();
800  long int offset = (long int) floor(ran_frac * fNEntries);
801  LOG("Flux", pERROR) << "Setting flux driver to start looping over entries "
802  << "with offset of "<< offset;
803  fIEntry = fOffset = offset;
804 }
#define pERROR
Definition: Messenger.h:59
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
long int fOffset
start looping at entry fOffset
Definition: GJPARCNuFlux.h:128
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
long int fIEntry
current flux ntuple entry
Definition: GJPARCNuFlux.h:126
long int fNEntries
number of flux ntuple entries
Definition: GJPARCNuFlux.h:125
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:71
void GJPARCNuFlux::ResetCurrent ( void  )
private

Definition at line 895 of file GJPARCNuFlux.cxx.

References fgP4, fgPdgC, fgX4, fLoadedNeutrino, fPassThroughInfo, and genie::flux::GJPARCNuFluxPassThroughInfo::Reset().

Referenced by GenerateNext_weighted(), and Initialize().

896 {
897 // reset running values of neutrino pdg-code, 4-position & 4-momentum
898 // and the input ntuple leaves
899 
900  fLoadedNeutrino = false;
901 
902  fgPdgC = 0;
903  fgP4.SetPxPyPzE (0.,0.,0.,0.);
904  fgX4.SetXYZT (0.,0.,0.,0.);
905 
907 }
int fgPdgC
running generated nu pdg-code
Definition: GJPARCNuFlux.h:110
bool fLoadedNeutrino
set to true when GenerateNext_weighted has been called successfully
Definition: GJPARCNuFlux.h:141
TLorentzVector fgP4
running generated nu 4-momentum
Definition: GJPARCNuFlux.h:111
GJPARCNuFluxPassThroughInfo * fPassThroughInfo
Definition: GJPARCNuFlux.h:143
TLorentzVector fgX4
running generated nu 4-position
Definition: GJPARCNuFlux.h:112
void GJPARCNuFlux::SetDefaults ( void  )
private

Definition at line 868 of file GJPARCNuFlux.cxx.

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

Referenced by Initialize().

869 {
870 // - Set default neutrino species list (nue, nuebar, numu, numubar) and
871 // maximum energy (25 GeV).
872 // These defaults can be overwritten by user calls (at the driver init) to
873 // GJPARCNuFlux::SetMaxEnergy(double Ev) and
874 // GJPARCNuFlux::SetFluxParticles(const PDGCodeList & particles)
875 // - Set the default file normalization to 1E+21 POT
876 // - Set the default flux neutrino start z position at -5m (z=0 is the
877 // detector centre).
878 // - Set number of cycles to 1
879 
880  LOG("Flux", pNOTICE) << "Setting default GJPARCNuFlux driver options";
881 
882  PDGCodeList particles;
883  particles.push_back(kPdgNuMu);
884  particles.push_back(kPdgAntiNuMu);
885  particles.push_back(kPdgNuE);
886  particles.push_back(kPdgAntiNuE);
887 
888  this->SetFluxParticles(particles);
889  this->SetMaxEnergy(25./*GeV*/);
890  this->SetFilePOT(1E+21);
891  this->SetUpstreamZ(-5.0);
892  this->SetNumOfCycles(1);
893 }
const int kPdgNuE
Definition: PDGCodes.h:28
void SetUpstreamZ(double z0)
set flux neutrino initial z position (upstream of the detector)
const int kPdgAntiNuE
Definition: PDGCodes.h:29
const int kPdgNuMu
Definition: PDGCodes.h:30
void SetFluxParticles(const PDGCodeList &particles)
specify list of flux neutrino species
A list of PDG codes.
Definition: PDGCodeList.h:32
void SetMaxEnergy(double Ev)
specify maximum flx neutrino energy
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void SetFilePOT(double pot)
flux file norm is in /N POT/det [ND] or /N POT/cm^2 [FD]. Specify N (typically 1E+21) ...
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
void SetNumOfCycles(int n)
set how many times to cycle through the ntuple (default: 1 / n=0 means &#39;infinite&#39;) ...
#define pNOTICE
Definition: Messenger.h:61
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
void GJPARCNuFlux::SetFilePOT ( double  pot)

flux file norm is in /N POT/det [ND] or /N POT/cm^2 [FD]. Specify N (typically 1E+21)

Definition at line 751 of file GJPARCNuFlux.cxx.

References fFilePOT.

Referenced by main(), and SetDefaults().

752 {
753 // The flux normalization is in /N POT/det [ND] or /N POT/cm^2 [FD]
754 // This method specifies N (typically 1E+21)
755 
756  fFilePOT = pot;
757 }
double fFilePOT
file POT normalization, typically 1E+21
Definition: GJPARCNuFlux.h:131
void GJPARCNuFlux::SetFluxParticles ( const PDGCodeList particles)

specify list of flux neutrino species

Definition at line 732 of file GJPARCNuFlux.cxx.

References genie::PDGCodeList::Copy(), fPdgCList, LOG, and pINFO.

Referenced by main(), and SetDefaults().

733 {
734  if(!fPdgCList) {
735  fPdgCList = new PDGCodeList;
736  }
737  fPdgCList->Copy(particles);
738 
739  LOG("Flux", pINFO)
740  << "Declared list of neutrino species: " << *fPdgCList;
741 }
A list of PDG codes.
Definition: PDGCodeList.h:32
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pINFO
Definition: Messenger.h:62
PDGCodeList * fPdgCList
list of neutrino pdg-codes
Definition: GJPARCNuFlux.h:108
void Copy(const PDGCodeList &list)
copy / print
void GJPARCNuFlux::SetMaxEnergy ( double  Ev)

specify maximum flx neutrino energy

Definition at line 743 of file GJPARCNuFlux.cxx.

References fMaxEv, LOG, and pINFO.

Referenced by SetDefaults().

744 {
745  fMaxEv = TMath::Max(0.,Ev);
746 
747  LOG("Flux", pINFO)
748  << "Declared maximum flux neutrino energy: " << fMaxEv;
749 }
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pINFO
Definition: Messenger.h:62
double fMaxEv
maximum energy
Definition: GJPARCNuFlux.h:107
void GJPARCNuFlux::SetNumOfCycles ( int  n)

set how many times to cycle through the ntuple (default: 1 / n=0 means 'infinite')

Definition at line 769 of file GJPARCNuFlux.cxx.

References fNCycles.

Referenced by main(), and SetDefaults().

770 {
771 // The flux ntuples can be recycled for a number of times to boost generated
772 // event statistics without requiring enormous beam simulation statistics.
773 // That option determines how many times the driver is going to cycle through
774 // the input flux ntuple.
775 // With n=0 the flux ntuple will be recycled an infinite amount of times so
776 // that the event generation loop can exit only on a POT or event num check.
777 
778  fNCycles = TMath::Max(0, n);
779 }
int fNCycles
how many times to cycle through the flux ntuple
Definition: GJPARCNuFlux.h:133
void GJPARCNuFlux::SetUpstreamZ ( double  z0)

set flux neutrino initial z position (upstream of the detector)

Definition at line 759 of file GJPARCNuFlux.cxx.

References fZ0.

Referenced by main(), and SetDefaults().

760 {
761 // The flux neutrino position (x,y) is given at the detector coord system
762 // at z=0. This method sets the preferred starting z position upstream of
763 // the upstream detector face. Each flux neutrino will be backtracked from
764 // z=0 to the input z0.
765 
766  fZ0 = z0;
767 }
double fZ0
configurable starting z position for each flux neutrino (in detector coord system) ...
Definition: GJPARCNuFlux.h:132
double genie::flux::GJPARCNuFlux::SumWeight ( void  ) const
inline

intergated weight for flux neutrinos looped so far

Definition at line 89 of file GJPARCNuFlux.h.

References fSumWeight.

double genie::flux::GJPARCNuFlux::Weight ( void  )
inlinevirtual

returns the flux neutrino weight (if any)

Implements genie::GFluxI.

Definition at line 63 of file GJPARCNuFlux.h.

References fMaxWeight, and fNorm.

Referenced by GenerateNext(), and GenerateNext_weighted().

63 { return fNorm / fMaxWeight; }
double fMaxWeight
max flux neutrino weight in input file for the specified detector location
Definition: GJPARCNuFlux.h:130
double fNorm
current flux ntuple normalisation
Definition: GJPARCNuFlux.h:129

Member Data Documentation

string genie::flux::GJPARCNuFlux::fDetLoc
private

input detector location ('sk','nd1','nd2',...)

Definition at line 120 of file GJPARCNuFlux.h.

Referenced by Initialize(), and LoadBeamSimData().

int genie::flux::GJPARCNuFlux::fDetLocId
private

input detector location id (fDetLoc -> jnubeam idfd)

Definition at line 121 of file GJPARCNuFlux.h.

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

long int genie::flux::GJPARCNuFlux::fEntriesThisCycle
private

keep track of number of entries used so far for this cycle

Definition at line 127 of file GJPARCNuFlux.h.

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

double genie::flux::GJPARCNuFlux::fFilePOT
private

file POT normalization, typically 1E+21

Definition at line 131 of file GJPARCNuFlux.h.

Referenced by Initialize(), POT_1cycle(), and SetFilePOT().

bool genie::flux::GJPARCNuFlux::fGenerateWeighted
private

generate weighted/deweighted flux neutrinos (default is false)

Definition at line 139 of file GJPARCNuFlux.h.

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

TLorentzVector genie::flux::GJPARCNuFlux::fgP4
private

running generated nu 4-momentum

Definition at line 111 of file GJPARCNuFlux.h.

Referenced by GenerateNext_weighted(), Momentum(), and ResetCurrent().

int genie::flux::GJPARCNuFlux::fgPdgC
private

running generated nu pdg-code

Definition at line 110 of file GJPARCNuFlux.h.

Referenced by GenerateNext_weighted(), PdgCode(), and ResetCurrent().

TLorentzVector genie::flux::GJPARCNuFlux::fgX4
private

running generated nu 4-position

Definition at line 112 of file GJPARCNuFlux.h.

Referenced by GenerateNext_weighted(), Position(), and ResetCurrent().

int genie::flux::GJPARCNuFlux::fICycle
private

current cycle

Definition at line 134 of file GJPARCNuFlux.h.

Referenced by Clear(), End(), GenerateNext(), GenerateNext_weighted(), Initialize(), and LoadBeamSimData().

long int genie::flux::GJPARCNuFlux::fIEntry
private

current flux ntuple entry

Definition at line 126 of file GJPARCNuFlux.h.

Referenced by Clear(), GenerateNext_weighted(), Index(), Initialize(), and RandomOffset().

bool genie::flux::GJPARCNuFlux::fIsFDLoc
private

input location is a 'far' detector location?

Definition at line 123 of file GJPARCNuFlux.h.

Referenced by Initialize(), and LoadBeamSimData().

bool genie::flux::GJPARCNuFlux::fIsNDLoc
private

input location is a 'near' detector location?

Definition at line 124 of file GJPARCNuFlux.h.

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

bool genie::flux::GJPARCNuFlux::fLoadedNeutrino
private

set to true when GenerateNext_weighted has been called successfully

Definition at line 141 of file GJPARCNuFlux.h.

Referenced by GenerateNext_weighted(), Index(), Initialize(), and ResetCurrent().

double genie::flux::GJPARCNuFlux::fMaxEv
private

maximum energy

Definition at line 107 of file GJPARCNuFlux.h.

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

double genie::flux::GJPARCNuFlux::fMaxWeight
private

max flux neutrino weight in input file for the specified detector location

Definition at line 130 of file GJPARCNuFlux.h.

Referenced by GenerateNext_weighted(), Initialize(), LoadBeamSimData(), POT_1cycle(), and Weight().

int genie::flux::GJPARCNuFlux::fNCycles
private

how many times to cycle through the flux ntuple

Definition at line 133 of file GJPARCNuFlux.h.

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

int genie::flux::GJPARCNuFlux::fNDetLocIdFound
private

per cycle keep track of the number of fDetLocId are found - if this is zero will exit job

Definition at line 122 of file GJPARCNuFlux.h.

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

long int genie::flux::GJPARCNuFlux::fNEntries
private

number of flux ntuple entries

Definition at line 125 of file GJPARCNuFlux.h.

Referenced by End(), GenerateNext_weighted(), Index(), Initialize(), LoadBeamSimData(), and RandomOffset().

long int genie::flux::GJPARCNuFlux::fNNeutrinos
private

number of flux neutrinos thrown so far

Definition at line 136 of file GJPARCNuFlux.h.

Referenced by Clear(), GenerateNext_weighted(), Initialize(), NFluxNeutrinos(), and POT_curravg().

long int genie::flux::GJPARCNuFlux::fNNeutrinosTot1c
private

total number of flux neutrinos to be thrown / cycle

Definition at line 138 of file GJPARCNuFlux.h.

Referenced by Initialize(), LoadBeamSimData(), and POT_curravg().

double genie::flux::GJPARCNuFlux::fNorm
private

current flux ntuple normalisation

Definition at line 129 of file GJPARCNuFlux.h.

Referenced by GenerateNext_weighted(), Initialize(), LoadBeamSimData(), and Weight().

TChain* genie::flux::GJPARCNuFlux::fNuFluxChain
private

input flux ntuple

Definition at line 116 of file GJPARCNuFlux.h.

Referenced by GenerateNext_weighted(), Initialize(), LoadBeamSimData(), POT_1cycle(), and POT_curravg().

TFile* genie::flux::GJPARCNuFlux::fNuFluxFile
private

input flux file

Definition at line 114 of file GJPARCNuFlux.h.

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

TChain* genie::flux::GJPARCNuFlux::fNuFluxSumChain
private

input summary ntuple for flux file. This tree is only present for later flux versions > 10a

Definition at line 118 of file GJPARCNuFlux.h.

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

TTree* genie::flux::GJPARCNuFlux::fNuFluxSumTree
private

input summary ntuple for flux file. This tree is only present for later flux versions > 10a

Definition at line 117 of file GJPARCNuFlux.h.

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

TTree* genie::flux::GJPARCNuFlux::fNuFluxTree
private

input flux ntuple

Definition at line 115 of file GJPARCNuFlux.h.

Referenced by GenerateNext_weighted(), Initialize(), LoadBeamSimData(), POT_1cycle(), and POT_curravg().

bool genie::flux::GJPARCNuFlux::fNuFluxUsingTree
private

are we using a TTree or a TChain to view the input flux file?

Definition at line 119 of file GJPARCNuFlux.h.

Referenced by GenerateNext_weighted(), Initialize(), LoadBeamSimData(), POT_1cycle(), and POT_curravg().

long int genie::flux::GJPARCNuFlux::fOffset
private

start looping at entry fOffset

Definition at line 128 of file GJPARCNuFlux.h.

Referenced by Clear(), GenerateNext_weighted(), Initialize(), and RandomOffset().

GJPARCNuFluxPassThroughInfo* genie::flux::GJPARCNuFlux::fPassThroughInfo
private
PDGCodeList* genie::flux::GJPARCNuFlux::fPdgCList
private

list of neutrino pdg-codes

Definition at line 108 of file GJPARCNuFlux.h.

Referenced by CleanUp(), FluxParticles(), GenerateNext_weighted(), Initialize(), and SetFluxParticles().

double genie::flux::GJPARCNuFlux::fSumWeight
private

sum of weights for neutrinos thrown so far

Definition at line 135 of file GJPARCNuFlux.h.

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

double genie::flux::GJPARCNuFlux::fSumWeightTot1c
private

total sum of weights for neutrinos to be thrown / cycle

Definition at line 137 of file GJPARCNuFlux.h.

Referenced by Initialize(), and LoadBeamSimData().

bool genie::flux::GJPARCNuFlux::fUseRandomOffset
private

whether set random starting point when looping over flux ntuples

Definition at line 140 of file GJPARCNuFlux.h.

Referenced by DisableOffset(), Initialize(), and LoadBeamSimData().

double genie::flux::GJPARCNuFlux::fZ0
private

configurable starting z position for each flux neutrino (in detector coord system)

Definition at line 132 of file GJPARCNuFlux.h.

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


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