GSimpleNtpFlux.h
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \class genie::flux::GSimpleNtpFlux
5 
6 \brief A GENIE flux driver using a simple ntuple format
7 
8 \author Robert Hatcher <rhatcher \at fnal.gov>
9  Fermi National Accelerator Laboratory
10 
11 \created Jan 25, 2010
12 
13 \cpright Copyright (c) 2003-2019, The GENIE Collaboration
14  For the full text of the license visit http://copyright.genie-mc.org
15  or see $GENIE/LICENSE
16 */
17 //____________________________________________________________________________
18 
19 #ifndef _SIMPLE_NTP_FLUX_H_
20 #define _SIMPLE_NTP_FLUX_H_
21 
22 #include <string>
23 #include <iostream>
24 #include <vector>
25 #include <set>
26 
27 #include <TVector3.h>
28 #include <TLorentzVector.h>
29 #include <TLorentzRotation.h>
30 
35 
36 class TFile;
37 class TChain;
38 class TTree;
39 class TBranch;
40 
41 using std::string;
42 using std::ostream;
43 
44 namespace genie {
45 namespace flux {
46 
47 class GSimpleNtpEntry;
48 ostream & operator << (ostream & stream, const GSimpleNtpEntry & info);
49 
50 /// Small persistable C-struct -like classes that makes up the SimpleNtpFlux
51 /// ntuple. This is only valid for a particular flux window (no reweighting,
52 /// no coordinate transformation available).
53 ///
54 /// Order elements from largest to smallest for ROOT alignment purposes
55 
56 /// GSimpleNtpEntry
57 /// =========================
58 /// This is the only required branch ("entry") of the "flux" tree
60  public:
62  /* allow default copy constructor ... for now nothing special
63  GSimpleNtpEntry(const GSimpleNtpEntry & info);
64  */
65  virtual ~GSimpleNtpEntry() { };
66  void Reset();
67  void Print(const Option_t* opt = "") const;
68  friend ostream & operator << (ostream & stream, const GSimpleNtpEntry & info);
69 
70  Double_t wgt; ///< nu weight
71 
72  Double_t vtxx; ///< x position in lab frame
73  Double_t vtxy; ///< y position in lab frame
74  Double_t vtxz; ///< z position in lab frame
75  Double_t dist; ///< distance from hadron decay
76 
77  Double_t px; ///< x momentum in lab frame
78  Double_t py; ///< y momentum in lab frame
79  Double_t pz; ///< z momentum in lab frame
80  Double_t E; ///< energy in lab frame
81 
82  Int_t pdg; ///< nu pdg-code
83  UInt_t metakey; ///< key to meta data
84 
85  ClassDef(GSimpleNtpEntry,1)
86  };
87 
88 
89  class GSimpleNtpNuMI;
90  ostream & operator << (ostream & stream, const GSimpleNtpNuMI & info);
91 
92 /// GSimpleNtpNuMI
93 /// =========================
94 /// Additional elements for NuMI (allow SKZP reweighting and reference
95 /// back to original GNuMI flux entries) as "numi" branch
97  public:
99  virtual ~GSimpleNtpNuMI() { };
100  void Reset();
101  void Print(const Option_t* opt = "") const;
102  friend ostream & operator << (ostream & stream, const GSimpleNtpNuMI & info);
103  Double_t tpx; ///< parent particle px at target exit
104  Double_t tpy;
105  Double_t tpz;
106  Double_t vx; ///< vertex position of hadron/muon decay
107  Double_t vy;
108  Double_t vz;
109  Double_t pdpx; ///< nu parent momentum at time of decay
110  Double_t pdpy;
111  Double_t pdpz;
112  Double_t pppx; ///< nu parent momentum at production point
113  Double_t pppy;
114  Double_t pppz;
115 
116  Int_t ndecay; ///< decay mode
117  Int_t ptype; ///< parent type (PDG)
118  Int_t ppmedium; ///< tracking medium where parent was produced
119  Int_t tptype; ///< parent particle type at target exit
120 
121  Int_t run; ///<
122  Int_t evtno; ///<
123  Int_t entryno; ///<
124 
125  ClassDef(GSimpleNtpNuMI,3)
126  };
127 
128 
129  class GSimpleNtpAux;
130  ostream & operator << (ostream & stream, const GSimpleNtpAux & info);
131 
132 /// GSimpleNtpAux
133 /// =========================
134 /// Additional elements for expansion as "aux" branch
136  public:
137  GSimpleNtpAux();
138  virtual ~GSimpleNtpAux() { };
139  void Reset();
140  void Print(const Option_t* opt = "") const;
141  friend ostream & operator << (ostream & stream, const GSimpleNtpAux & info);
142 
143  std::vector<Int_t> auxint; ///< additional ints associated w/ entry
144  std::vector<Double_t> auxdbl; ///< additional doubles associated w/ entry
145 
146  ClassDef(GSimpleNtpAux,1)
147  };
148 
149 
150  class GSimpleNtpMeta;
151  ostream & operator << (ostream & stream, const GSimpleNtpMeta & info);
152 
153 /// GSimpleNtpMeta
154 /// =========================
155 /// A small persistable C-struct -like class that holds metadata
156 /// about the the SimpleNtpFlux ntple.
157 ///
158  class GSimpleNtpMeta: public TObject {
159  public:
160  GSimpleNtpMeta();
161  /* allow default copy constructor ... for now nothing special
162  GSimpleNtpMeta(const GSimpleNtpMeta & info);
163  */
164  virtual ~GSimpleNtpMeta();
165 
166  void Reset();
167  void AddFlavor(Int_t nupdg);
168  void Print(const Option_t* opt = "") const;
169  friend ostream & operator << (ostream & stream, const GSimpleNtpMeta & info);
170 
171  std::vector<Int_t> pdglist; ///< list of neutrino flavors
172 
173  Double_t maxEnergy; ///< maximum energy
174  Double_t minWgt; ///< minimum weight
175  Double_t maxWgt; ///< maximum weight
176  Double_t protons; ///< represented number of protons-on-target
177 
178  Double_t windowBase[3]; ///< x,y,z position of window base point
179  Double_t windowDir1[3]; ///< dx,dy,dz of window direction 1
180  Double_t windowDir2[3]; ///< dx,dy,dz of window direction 2
181 
182  std::vector<std::string> auxintname; ///< tagname of aux ints associated w/ entry
183  std::vector<std::string> auxdblname; ///< tagname of aux doubles associated w/ entry
184  std::vector<std::string> infiles; ///< list of input files
185 
186  Int_t seed; ///< random seed used in generation
187  UInt_t metakey; ///< index key to tie to individual entries
188 
189  static UInt_t mxfileprint; ///< allow user to limit # of files to print
190 
191  ClassDef(GSimpleNtpMeta,1)
192  };
193 
194 
195 /// GSimpleNtpFlux:
196 /// ==========
197 /// An implementation of the GFluxI interface that provides NuMI flux
198 ///
200  : public genie::GFluxI
203 {
204 
205 public :
206  GSimpleNtpFlux();
207  ~GSimpleNtpFlux();
208 
209  // Methods implementing the GENIE GFluxI interface, required for integrating
210  // the NuMI neutrino flux simulations with the GENIE event generation drivers
211 
212  const PDGCodeList & FluxParticles (void) { return *fPdgCList; }
213  double MaxEnergy (void) { return fMaxEv; }
214  bool GenerateNext (void);
215  int PdgCode (void) { return fCurEntry->pdg; }
216  double Weight (void) { return fWeight; }
217  const TLorentzVector & Momentum (void) { return fP4; }
218  const TLorentzVector & Position (void) { return fX4; }
219  bool End (void) { return fEnd; }
220  long int Index (void) { return fIEntry; }
221  void Clear (Option_t * opt);
222  void GenerateWeighted (bool gen_weighted);
223 
224  // Methods specific to the NuMI flux driver,
225  // for configuration/initialization of the flux & event generation drivers
226  // and and for passing-through flux information (e.g. neutrino parent decay
227  // kinematics) not used by the generator but required by analyses/processing
228  // further downstream
229 
230  //
231  // information about or actions on current entry
232  //
234  GetCurrentEntry(void) { return fCurEntry; } ///< GSimpleNtpEntry
236  GetCurrentNuMI(void) { return fCurNuMI; } ///< GSimpleNtpNuMI
238  GetCurrentAux(void) { return fCurAux; } ///< GSimpleNtpAux
240  GetCurrentMeta(void) { return fCurMeta; } ///< GSimpleNtpMeta
241 
242  // allow access to main tree so we can call Branch() to retrieve extra stuff
243  TChain*
244  GetFluxTChain(void) { return fNuFluxTree; } ///<
245 
246  double GetDecayDist() const; ///< dist (user units) from dk to current pos
247  void MoveToZ0(double z0); ///< move ray origin to user coord Z0
248 
249  //
250  // information about the current state
251  //
252  virtual double GetTotalExposure() const; ///< GFluxExposureI interface
253  virtual long int NFluxNeutrinos() const; ///< # of rays generated
254 
255  double UsedPOTs(void) const; ///< # of protons-on-target used
256 
257  long int NEntriesUsed(void) const { return fNEntriesUsed; } ///< number of entries read from files
258  double SumWeight(void) const { return fSumWeight; } ///< integrated weight for flux neutrinos looped so far
259 
260  void PrintCurrent(void); ///< print current entry from leaves
261  void PrintConfig(); ///< print the current configuration
262 
263  std::vector<std::string> GetFileList(); ///< list of files currently part of chain
264 
265  //
266  // GFluxFileConfigI interface
267  //
268  virtual void LoadBeamSimData(const std::vector<string>& filenames,
269  const std::string& det_loc);
270  using GFluxFileConfigI::LoadBeamSimData; // inherit the rest
271  virtual void GetBranchInfo(std::vector<std::string>& branchNames,
272  std::vector<std::string>& branchClassNames,
273  std::vector<void**>& branchObjPointers);
274  virtual TTree* GetMetaDataTree();
275 
276  //
277  // configuration of GSimpleNtpFlux
278  //
279 
280  void SetRequestedBranchList(string blist="entry,numi,aux") { fNuFluxBranchRequest = blist; }
281 
282  void SetMaxEnergy(double Ev); ///< specify maximum flx neutrino energy
283 
284  void SetGenWeighted(bool genwgt=false) { fGenWeighted = genwgt; } ///< toggle whether GenerateNext() returns weight=1 flux (initial default false)
285 
286  void SetEntryReuse(long int nuse=1); ///< # of times to use entry before moving to next
287 
288  void ProcessMeta(void); ///< scan for max flux energy, weight
289 
290  void GetFluxWindow(TVector3& p1, TVector3& p2, TVector3& p3) const; ///< 3 points define a plane in beam coordinate
291 
292 private:
293 
294  // Private methods
295  //
296  bool GenerateNext_weighted (void);
297  void Initialize (void);
298  void SetDefaults (void);
299  void CleanUp (void);
300  void ResetCurrent (void);
301  void AddFile (TTree* fluxtree, TTree* metatree, string fname);
302  bool OptionalAttachBranch (std::string bname);
303  void CalcEffPOTsPerNu (void);
304  void ScanMeta (void);
305 
306  // Private data members
307  //
308  double fMaxEv; ///< maximum energy
309  bool fEnd; ///< end condition reached
310 
311  std::vector<string> fNuFluxFilePatterns; ///< (potentially wildcarded) path(s)
312  string fNuFluxBranchRequest; ///< list of requested branches "entry,numi,au"
313  TChain* fNuFluxTree; ///< TTree // REF ONLY
314  TChain* fNuMetaTree; ///< TTree // REF ONLY
315 
316  int fNFiles; ///< number of files in chain
317  Long64_t fNEntries; ///< number of flux ntuple entries
318  Long64_t fIEntry; ///< current flux ntuple entry
319  Int_t fIFileNumber; ///< which file for the current entry
320 
321  Double_t fFilePOTs; ///< # of protons-on-target represented by all files
322 
323  double fWeight; ///< current neutrino weight
324  double fMaxWeight; ///< max flux neutrino weight in input file
325 
326  long int fNUse; ///< how often to use same entry in a row
327  long int fIUse; ///< current # of times an entry has been used
328 
329  double fSumWeight; ///< sum of weights for nus thrown so far
330  long int fNNeutrinos; ///< number of flux neutrinos thrown so far
331  long int fNEntriesUsed; ///< number of entries read from files
332  double fEffPOTsPerNu; ///< what a entry is worth ...
333  double fAccumPOTs; ///< POTs used so far
334 
335  bool fGenWeighted; ///< does GenerateNext() give weights?
336  bool fAlreadyUnwgt; ///< are input files already unweighted
337  // i.e. are all entry "wgt" values = 1
338  bool fAllFilesMeta; ///< do all files in chain have meta data
339 
340  GSimpleNtpEntry* fCurEntry; ///< current entry
341  GSimpleNtpNuMI* fCurNuMI; ///< current "numi" branch extra info
342  GSimpleNtpAux* fCurAux; ///< current "aux" branch extra info
343  TLorentzVector fP4; ///< reconstituted p4 vector
344  TLorentzVector fX4; ///< reconstituted position vector
345  GSimpleNtpMeta* fCurMeta; ///< current meta data
346 
347  GSimpleNtpEntry* fCurEntryCopy; ///< current entry
348  GSimpleNtpNuMI* fCurNuMICopy; ///< current "numi" branch extra info
349  GSimpleNtpAux* fCurAuxCopy; ///< current "aux" branch extra info
350 
351 };
352 
353 } // flux namespace
354 } // genie namespace
355 
356 #endif // _SIMPLE_NTP_FLUX_H_
void SetGenWeighted(bool genwgt=false)
toggle whether GenerateNext() returns weight=1 flux (initial default false)
const PDGCodeList & FluxParticles(void)
declare list of flux neutrinos that can be generated (for init. purposes)
const TLorentzVector & Position(void)
returns the flux neutrino 4-position (note: expect SI rather than physical units) ...
Double_t E
energy in lab frame
const XML_Char XML_Encoding * info
Definition: expat.h:530
GSimpleNtpAux * fCurAux
current "aux" branch extra info
Int_t fIFileNumber
which file for the current entry
std::vector< string > fNuFluxFilePatterns
(potentially wildcarded) path(s)
long int fNNeutrinos
number of flux neutrinos thrown so far
virtual void LoadBeamSimData(const std::vector< std::string > &filenames, const std::string &det_loc)=0
int PdgCode(void)
returns the flux neutrino pdg code
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
GSimpleNtpMeta * fCurMeta
current meta data
const TLorentzVector & Momentum(void)
returns the flux neutrino 4-momentum
Double_t pdpx
nu parent momentum at time of decay
long int fIUse
current # of times an entry has been used
Double_t px
x momentum in lab frame
A GENIE flux driver using a simple ntuple format.
double fSumWeight
sum of weights for nus thrown so far
Double_t vtxy
y position in lab frame
def AddFile(fout, deets, f)
friend ostream & operator<<(ostream &stream, const GSimpleNtpEntry &info)
Double_t vx
vertex position of hadron/muon decay
Double_t maxEnergy
maximum energy
std::vector< std::string > auxintname
tagname of aux ints associated w/ entry
double MaxEnergy(void)
declare the max flux neutrino energy that can be generated (for init. purposes)
Int_t tptype
parent particle type at target exit
bool End(void)
true if no more flux nu&#39;s can be thrown (eg reaching end of beam sim ntuples)
std::vector< Double_t > auxdbl
additional doubles associated w/ entry
std::vector< Int_t > auxint
additional ints associated w/ entry
void Initialize(Bool_t useTMVAStyle=kTRUE)
Definition: tmvaglob.C:184
Loaders::FluxType flux
Double_t protons
represented number of protons-on-target
A list of PDG codes.
Definition: PDGCodeList.h:33
Int_t seed
random seed used in generation
ostream & operator<<(ostream &stream, const genie::flux::GJPARCNuFluxPassThroughInfo &info)
const genie::flux::GSimpleNtpEntry * GetCurrentEntry(void)
GSimpleNtpEntry.
Double_t pppx
nu parent momentum at production point
const genie::flux::GSimpleNtpMeta * GetCurrentMeta(void)
GSimpleNtpMeta.
void SetRequestedBranchList(string blist="entry,numi,aux")
std::vector< Int_t > pdglist
list of neutrino flavors
TLorentzVector fX4
reconstituted position vector
bool fAllFilesMeta
do all files in chain have meta data
Double_t fFilePOTs
of protons-on-target represented by all files
bool fGenWeighted
does GenerateNext() give weights?
Double_t vtxz
z position in lab frame
GENIE interface for uniform flux exposure iterface.
Double_t minWgt
minimum weight
Int_t ppmedium
tracking medium where parent was produced
GSimpleNtpNuMI * fCurNuMICopy
current "numi" branch extra info
long int fNEntriesUsed
number of entries read from files
std::vector< std::string > infiles
list of input files
Double_t tpx
parent particle px at target exit
GSimpleNtpNuMI * fCurNuMI
current "numi" branch extra info
double SumWeight(void) const
integrated weight for flux neutrinos looped so far
int fNFiles
number of files in chain
double fMaxWeight
max flux neutrino weight in input file
GENIE interface for uniform flux exposure iterface.
double fWeight
current neutrino weight
TChain * fNuMetaTree
TTree // REF ONLY.
GSimpleNtpEntry * fCurEntry
current entry
const genie::flux::GSimpleNtpAux * GetCurrentAux(void)
GSimpleNtpAux.
Double_t maxWgt
maximum weight
TLorentzVector fP4
reconstituted p4 vector
UInt_t metakey
index key to tie to individual entries
bool fAlreadyUnwgt
are input files already unweighted
Long64_t fIEntry
current flux ntuple entry
Long64_t fNEntries
number of flux ntuple entries
Double_t vtxx
x position in lab frame
Double_t pz
z momentum in lab frame
long int Index(void)
returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current ent...
GSimpleNtpAux * fCurAuxCopy
current "aux" branch extra info
long int fNUse
how often to use same entry in a row
long int NEntriesUsed(void) const
number of entries read from files
Double_t dist
distance from hadron decay
GSimpleNtpEntry * fCurEntryCopy
current entry
double fEffPOTsPerNu
what a entry is worth ...
const genie::flux::GSimpleNtpNuMI * GetCurrentNuMI(void)
GSimpleNtpNuMI.
string fNuFluxBranchRequest
list of requested branches "entry,numi,au"
bool fEnd
end condition reached
double fMaxEv
maximum energy
UInt_t metakey
key to meta data
void Print(const Option_t *opt="") const
Double_t py
y momentum in lab frame
double Weight(void)
returns the flux neutrino weight (if any)
double fAccumPOTs
POTs used so far.
static UInt_t mxfileprint
allow user to limit # of files to print
std::vector< std::string > auxdblname
tagname of aux doubles associated w/ entry
TChain * fNuFluxTree
TTree // REF ONLY.
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
Int_t ptype
parent type (PDG)
enum BeamMode string