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