GJPARCNuFlux.h
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \class genie::flux::GJPARCNuFlux
5 
6 \brief A GENIE flux driver encapsulating the JPARC neutrino flux.
7  It reads-in the official JPARC neutrino flux ntuples.
8 
9 \ref See http://jnusrv01.kek.jp/internal/t2k/nubeam/flux/
10  (Note: T2K internal)
11 
12 \author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
13  University of Liverpool & STFC Rutherford Appleton Lab
14 
15 \created Feb 04, 2008
16 
17 \cpright Copyright (c) 2003-2019, The GENIE Collaboration
18  For the full text of the license visit http://copyright.genie-mc.org
19  or see $GENIE/LICENSE
20 */
21 //____________________________________________________________________________
22 
23 #ifndef _GJPARC_NEUTRINO_FLUX_H_
24 #define _GJPARC_NEUTRINO_FLUX_H_
25 
26 #include <string>
27 #include <iostream>
28 
29 #include <TLorentzVector.h>
30 #include <TVector3.h>
31 #include <TH1D.h>
32 
35 
36 class TFile;
37 class TTree;
38 class TChain;
39 class TBranch;
40 
41 using std::string;
42 using std::ostream;
43 
44 namespace genie {
45 namespace flux {
46 
47 class GJPARCNuFluxPassThroughInfo;
48 
49 ostream & operator << (ostream & stream, const GJPARCNuFluxPassThroughInfo & info);
50 
51 class GJPARCNuFlux: public GFluxI {
52 
53 public :
54  GJPARCNuFlux();
55  ~GJPARCNuFlux();
56 
57  // Methods implementing the GENIE GFluxI interface, required for integrating
58  // the JPARC neutrino flux simulations with the GENIE event generation drivers
59 
60  const PDGCodeList & FluxParticles (void) { return *fPdgCList; }
61  double MaxEnergy (void) { return fMaxEv; }
62  bool GenerateNext (void);
63  int PdgCode (void) { return fgPdgC; }
64  double Weight (void) { return fNorm / fMaxWeight; }
65  const TLorentzVector & Momentum (void) { return fgP4; }
66  const TLorentzVector & Position (void) { return fgX4; }
67  bool End (void) { return fEntriesThisCycle >= fNEntries
68  && fICycle == fNCycles && fNCycles > 0; }
69  long int Index (void);
70  void Clear (Option_t * opt);
71  void GenerateWeighted (bool gen_weighted = true);
72 
73  // Methods specific to the JPARC flux driver,
74  // for configuration/initialization of the flux & event generation drivers and
75  // and for passing-through flux information (eg neutrino parent decay kinematics)
76  // not used by the generator but required by analyses/processing further upstream
77 
78  bool LoadBeamSimData (string filename, string det_loc); ///< load a jnubeam root flux ntuple
79  void SetFluxParticles (const PDGCodeList & particles); ///< specify list of flux neutrino species
80  void SetMaxEnergy (double Ev); ///< specify maximum flx neutrino energy
81  void SetFilePOT (double pot); ///< flux file norm is in /N POT/det [ND] or /N POT/cm^2 [FD]. Specify N (typically 1E+21)
82  void SetUpstreamZ (double z0); ///< set flux neutrino initial z position (upstream of the detector)
83  void SetNumOfCycles (int n); ///< set how many times to cycle through the ntuple (default: 1 / n=0 means 'infinite')
84  void DisableOffset (void){fUseRandomOffset = false;} ///< switch off random offset, must be called before LoadBeamSimData to have any effect
85  void RandomOffset (void); ///< choose a random offset as starting entry in flux ntuple
86 
87  double POT_1cycle (void); ///< flux POT per cycle
88  double POT_curravg (void); ///< current average POT
89  long int NFluxNeutrinos (void) const { return fNNeutrinos; } ///< number of flux neutrinos looped so far
90  double SumWeight (void) const { return fSumWeight; } ///< intergated weight for flux neutrinos looped so far
91 
93  PassThroughInfo(void) { return *fPassThroughInfo; } ///< GJPARCNuFluxPassThroughInfo
94 
95 private:
96 
97  // Private methods
98  //
99  bool GenerateNext_weighted (void);
100  void Initialize (void);
101  void SetDefaults (void);
102  void CleanUp (void);
103  void ResetCurrent (void);
104  int DLocName2Id (string name);
105 
106  // Private data members
107  //
108  double fMaxEv; ///< maximum energy
109  PDGCodeList * fPdgCList; ///< list of neutrino pdg-codes
110 
111  int fgPdgC; ///< running generated nu pdg-code
112  TLorentzVector fgP4; ///< running generated nu 4-momentum
113  TLorentzVector fgX4; ///< running generated nu 4-position
114 
115  TFile * fNuFluxFile; ///< input flux file
116  TTree * fNuFluxTree; ///< input flux ntuple
117  TChain * fNuFluxChain; ///< input flux ntuple
118  TTree * fNuFluxSumTree; ///< input summary ntuple for flux file. This tree is only present for later flux versions > 10a
119  TChain * fNuFluxSumChain; ///< input summary ntuple for flux file. This tree is only present for later flux versions > 10a
120  bool fNuFluxUsingTree; ///< are we using a TTree or a TChain to view the input flux file?
121  string fDetLoc; ///< input detector location ('sk','nd1','nd2',...)
122  int fDetLocId; ///< input detector location id (fDetLoc -> jnubeam idfd)
123  int fNDetLocIdFound; ///< per cycle keep track of the number of fDetLocId are found - if this is zero will exit job
124  bool fIsFDLoc; ///< input location is a 'far' detector location?
125  bool fIsNDLoc; ///< input location is a 'near' detector location?
126  long int fNEntries; ///< number of flux ntuple entries
127  long int fIEntry; ///< current flux ntuple entry
128  long int fEntriesThisCycle; ///< keep track of number of entries used so far for this cycle
129  long int fOffset; ///< start looping at entry fOffset
130  double fNorm; ///< current flux ntuple normalisation
131  double fMaxWeight; ///< max flux neutrino weight in input file for the specified detector location
132  double fFilePOT; ///< file POT normalization, typically 1E+21
133  double fZ0; ///< configurable starting z position for each flux neutrino (in detector coord system)
134  int fNCycles; ///< how many times to cycle through the flux ntuple
135  int fICycle; ///< current cycle
136  double fSumWeight; ///< sum of weights for neutrinos thrown so far
137  long int fNNeutrinos; ///< number of flux neutrinos thrown so far
138  double fSumWeightTot1c; ///< total sum of weights for neutrinos to be thrown / cycle
139  long int fNNeutrinosTot1c; ///< total number of flux neutrinos to be thrown / cycle
140  bool fGenerateWeighted; ///< generate weighted/deweighted flux neutrinos (default is false)
141  bool fUseRandomOffset; ///< whether set random starting point when looping over flux ntuples
142  bool fLoadedNeutrino; ///< set to true when GenerateNext_weighted has been called successfully
143 
145 };
146 
147 
148 // A small persistable C-struct - like class that may be stored at an extra branch of
149 // the output event tree -alongside with the generated event branch- for use further
150 // upstream in the t2k analysis chain -eg beam reweighting etc-)
151 //
152 const int fNgmax = 12;
153 class GJPARCNuFluxPassThroughInfo: public TObject {
154 public:
158 
159  void Reset();
160  friend ostream & operator << (ostream & stream, const GJPARCNuFluxPassThroughInfo & info);
161 
162  long fluxentry;
163  string fluxfilename;
164  // Using an instance of this class the following datamembers are set
165  // directly as the branch addresses of jnubeam flux ntuples tree:
166  float Enu; // set to "Enu/F": Nu energy (GeV)
167  int ppid; // set to "ppid/I": Nu parent GEANT particle id
168  int mode; // set to "mode/I": Nu parent decay mode (see http://jnusrv01.kek.jp/internal/t2k/nubeam/flux/nemode.h)
169  float ppi; // set to "ppi/F": Nu parent momentum at its decay point (GeV)
170  float xpi[3]; // set to "xpi[3]/F": Nu parent position vector at decay (cm, in t2k global coord system)
171  float npi[3]; // set to "npi[3]/F": Nu parent direction vector at decay (in t2k global coord system)
172  float norm; // set to "norm/F": Weight to give flux in /N POT/det. [ND] or /N POT/cm^2 [FD], where is N is typically 1E+21
173  int nvtx0; // set to "nvtx0/I": Number of vtx where the nu. parent was produced (made obsolete by nd variable inroduced in 10d flux version)
174  float ppi0; // set to "ppi0/F": Nu parent momentum at its production point (GeV)
175  float xpi0[3]; // set to "xpi0[3]/F": Nu parent position vector at production (cm, in t2k global coord system)
176  float npi0[3]; // set to "npi0[3]/F": Nu parent direction vector at production (in t2k global coord system)
177  float rnu; // set to "rnu/F": Nu radial position (cm, in detector coord system)
178  float xnu; // set to "xnu/F": Nu x position (cm, in detector coord system)
179  float ynu; // set to "ynu/F": Nu y position (cm, in detector coord system)
180  float nnu[3]; // set to "nnu[3]/F": Nu direction (in t2k global coord system)
181  // New since JNuBeam 10a flux version.
182  float cospibm; // set to "cospibm/F": Nu parent direction cosine at decay (with respect to the beam direction)
183  float cospi0bm; // set to "cospi0bm/F": Nu parent direction cosine at production (with respect to the beam direction)
184  int idfd; // set to "idfd/I": Detector ID
185  unsigned char gipart; // set to "gipart/B": Primary particle ID
186  float gpos0[3]; // set to "gpos0[3]/F": Primary particle starting point
187  float gvec0[3]; // set to "gvec0[3]/F": Primary particle direction at the starting point
188  float gamom0; // set to "gamom0/F": Momentum of the primary particle at the starting point
189  // New since JNuBeam 10d and 11a flux version updates
190  int ng; // set to "ng/I": Number of parents (number of generations)
191  float gpx[fNgmax]; // set to "gpx[20]/F": Momentum X of each ancestor particle
192  float gpy[fNgmax]; // set to "gpy[20]/F": Momentum Y of each ancestor particle
193  float gpz[fNgmax]; // set to "gpz[20]/F": Momentum Z of each ancestor particle
194  float gcosbm[fNgmax]; // set to "gcosbm[20]/F": Cosine of the angle between the ancestor particle direction and the beam direction
195  float gvx[fNgmax]; // set to "gvx[20]/F": Vertex X position of each ancestor particle
196  float gvy[fNgmax]; // set to "gvy[20]/F": Vertex Y position of each ancestor particle
197  float gvz[fNgmax]; // set to "gvz[20]/F": Vertex Z position of each ancestor particle
198  int gpid[fNgmax]; // set to "gpid[20]/I": Particle ID of each ancestor particles
199  int gmec[fNgmax]; // set to "gmec[20]/I": Particle production mechanism of each ancestor particle
200  // Next five only present since 11a flux
201  int gmat[fNgmax]; // set to "gmat[fNgmax]/I": Material in which the particle originates
202  float gdistc[fNgmax]; // set to "gdistc[fNgmax]/F": Distance traveled through carbon
203  float gdistal[fNgmax]; // set to "gdista[fNgmax]/F": Distance traveled through aluminum
204  float gdistti[fNgmax];// set to "gdistti[fNgmax]/F": Distance traveled through titanium
205  float gdistfe[fNgmax];// set to "gdistte[fNgmax]/F": Distance traveled through iron
206  float Enusk; // set to "Enusk/F": "Enu" for SK
207  float normsk; // set to "normsk/F": "norm" for SK
208  float anorm; // set to "anorm/F": Norm component from ND acceptance calculation
209  // The following do not change per flux entry as is summary info for the flux
210  // file. For simplicity we just store per flux entry and accept the duplication.
211  float version; // set to "version/F": Jnubeam version
212  int tuneid; // set to "tuneid/I": Parameter set identifier
213  int ntrig; // set to "ntrig/I": Number of Triggers in simulation
214  int pint; // set to "pint/I": Interaction model ID
215  float bpos[2]; // set to "bpos[2]/F": Beam center position
216  float btilt[2]; // set to "btilt[2]/F": Beam Direction
217  float brms[2]; // set to "brms[2]/F": Beam RMS Width
218  float emit[2]; // set to "emit[2]/F": Beam Emittance
219  float alpha[2]; // set to "alpha[2]/F": Beam alpha parameter
220  float hcur[3]; // set to "hcur[3]/F": Horns 1, 2 and 3 Currents
221  int rand; // set to "rand/I": Random seed
222  int rseed[2]; // set to "rseed/I": Random seed
223 
224 ClassDef(GJPARCNuFluxPassThroughInfo,3)
225 };
226 
227 } // flux namespace
228 } // genie namespace
229 
230 #endif // _GJPARC_NEUTRINO_FLUX_H_
TTree * fNuFluxTree
input flux ntuple
Definition: GJPARCNuFlux.h:116
const XML_Char XML_Encoding * info
Definition: expat.h:530
const XML_Char * name
Definition: expat.h:151
double POT_curravg(void)
current average POT
int fgPdgC
running generated nu pdg-code
Definition: GJPARCNuFlux.h:111
double fSumWeightTot1c
total sum of weights for neutrinos to be thrown / cycle
Definition: GJPARCNuFlux.h:138
bool LoadBeamSimData(string filename, string det_loc)
load a jnubeam root flux ntuple
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
bool End(void)
true if no more flux nu&#39;s can be thrown (eg reaching end of beam sim ntuples)
Definition: GJPARCNuFlux.h:67
bool fLoadedNeutrino
set to true when GenerateNext_weighted has been called successfully
Definition: GJPARCNuFlux.h:142
void SetUpstreamZ(double z0)
set flux neutrino initial z position (upstream of the detector)
int DLocName2Id(string name)
string fDetLoc
input detector location (&#39;sk&#39;,&#39;nd1&#39;,&#39;nd2&#39;,...)
Definition: GJPARCNuFlux.h:121
TChain * fNuFluxSumChain
input summary ntuple for flux file. This tree is only present for later flux versions > 10a ...
Definition: GJPARCNuFlux.h:119
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:120
const TLorentzVector & Position(void)
returns the flux neutrino 4-position (note: expect SI rather than physical units) ...
Definition: GJPARCNuFlux.h:66
int PdgCode(void)
returns the flux neutrino pdg code
Definition: GJPARCNuFlux.h:63
double fFilePOT
file POT normalization, typically 1E+21
Definition: GJPARCNuFlux.h:132
TChain * fNuFluxChain
input flux ntuple
Definition: GJPARCNuFlux.h:117
double Weight(void)
returns the flux neutrino weight (if any)
Definition: GJPARCNuFlux.h:64
void GenerateWeighted(bool gen_weighted=true)
set whether to generate weighted or unweighted neutrinos
string filename
Definition: shutoffs.py:106
long int fOffset
start looping at entry fOffset
Definition: GJPARCNuFlux.h:129
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GJPARCNuFlux.h:137
void SetFluxParticles(const PDGCodeList &particles)
specify list of flux neutrino species
Loaders::FluxType flux
TTree * fNuFluxSumTree
input summary ntuple for flux file. This tree is only present for later flux versions > 10a ...
Definition: GJPARCNuFlux.h:118
A list of PDG codes.
Definition: PDGCodeList.h:33
A GENIE flux driver encapsulating the JPARC neutrino flux. It reads-in the official JPARC neutrino fl...
Definition: GJPARCNuFlux.h:51
void SetMaxEnergy(double Ev)
specify maximum flx neutrino energy
ostream & operator<<(ostream &stream, const genie::flux::GJPARCNuFluxPassThroughInfo &info)
long int fNNeutrinosTot1c
total number of flux neutrinos to be thrown / cycle
Definition: GJPARCNuFlux.h:139
long int NFluxNeutrinos(void) const
number of flux neutrinos looped so far
Definition: GJPARCNuFlux.h:89
TLorentzVector fgP4
running generated nu 4-momentum
Definition: GJPARCNuFlux.h:112
double MaxEnergy(void)
declare the max flux neutrino energy that can be generated (for init. purposes)
Definition: GJPARCNuFlux.h:61
int fDetLocId
input detector location id (fDetLoc -> jnubeam idfd)
Definition: GJPARCNuFlux.h:122
double SumWeight(void) const
intergated weight for flux neutrinos looped so far
Definition: GJPARCNuFlux.h:90
void SetFilePOT(double pot)
flux file norm is in /N POT/det [ND] or /N POT/cm^2 [FD]. Specify N (typically 1E+21) ...
void Reset()
const int fNgmax
Definition: GJPARCNuFlux.h:152
#define pot
int fNCycles
how many times to cycle through the flux ntuple
Definition: GJPARCNuFlux.h:134
int fNDetLocIdFound
per cycle keep track of the number of fDetLocId are found - if this is zero will exit job ...
Definition: GJPARCNuFlux.h:123
void DisableOffset(void)
switch off random offset, must be called before LoadBeamSimData to have any effect ...
Definition: GJPARCNuFlux.h:84
long int fIEntry
current flux ntuple entry
Definition: GJPARCNuFlux.h:127
bool fIsNDLoc
input location is a &#39;near&#39; detector location?
Definition: GJPARCNuFlux.h:125
GJPARCNuFluxPassThroughInfo * fPassThroughInfo
Definition: GJPARCNuFlux.h:144
bool fGenerateWeighted
generate weighted/deweighted flux neutrinos (default is false)
Definition: GJPARCNuFlux.h:140
TFile * fNuFluxFile
input flux file
Definition: GJPARCNuFlux.h:115
double fMaxWeight
max flux neutrino weight in input file for the specified detector location
Definition: GJPARCNuFlux.h:131
int fICycle
current cycle
Definition: GJPARCNuFlux.h:135
void Clear(Option_t *opt)
reset state variables based on opt
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 &#39;far&#39; detector location?
Definition: GJPARCNuFlux.h:124
double fZ0
configurable starting z position for each flux neutrino (in detector coord system) ...
Definition: GJPARCNuFlux.h:133
PDGCodeList * fPdgCList
list of neutrino pdg-codes
Definition: GJPARCNuFlux.h:109
bool GenerateNext(void)
generate the next flux neutrino (return false in err)
long int fNEntries
number of flux ntuple entries
Definition: GJPARCNuFlux.h:126
void SetNumOfCycles(int n)
set how many times to cycle through the ntuple (default: 1 / n=0 means &#39;infinite&#39;) ...
TLorentzVector fgX4
running generated nu 4-position
Definition: GJPARCNuFlux.h:113
bool fUseRandomOffset
whether set random starting point when looping over flux ntuples
Definition: GJPARCNuFlux.h:141
void RandomOffset(void)
choose a random offset as starting entry in flux ntuple
const TLorentzVector & Momentum(void)
returns the flux neutrino 4-momentum
Definition: GJPARCNuFlux.h:65
double fNorm
current flux ntuple normalisation
Definition: GJPARCNuFlux.h:130
double fSumWeight
sum of weights for neutrinos thrown so far
Definition: GJPARCNuFlux.h:136
long int fEntriesThisCycle
keep track of number of entries used so far for this cycle
Definition: GJPARCNuFlux.h:128
const GJPARCNuFluxPassThroughInfo & PassThroughInfo(void)
GJPARCNuFluxPassThroughInfo.
Definition: GJPARCNuFlux.h:93
double fMaxEv
maximum energy
Definition: GJPARCNuFlux.h:108
const PDGCodeList & FluxParticles(void)
declare list of flux neutrinos that can be generated (for init. purposes)
Definition: GJPARCNuFlux.h:60
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
enum BeamMode string