GDk2NuFlux.h
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \class genie::flux::GDk2NuFlux
5 
6 \brief An implementation of the GENIE GFluxI interface ("flux driver")
7  encapsulating reading/processing the "dk2nu" tree structure.
8 
9 \author Robert Hatcher <rhatcher \at fnal.gov>
10  Fermi National Accelerator Laboratory
11 
12 \created Nov 6, 2012
13 
14 \cpright Copyright (c) 2012, GENIE Neutrino MC Generator Collaboration
15  For the full text of the license visit http://copyright.genie-mc.org
16  or see $GENIE/LICENSE
17 */
18 //____________________________________________________________________________
19 
20 #ifndef _GDK2NUFLUX_H_
21 #define _GDK2NUFLUX_H_
22 
23 #include <string>
24 #include <iostream>
25 #include <vector>
26 #include <set>
27 #include <map>
28 
29 #include <TVector3.h>
30 #include <TLorentzVector.h>
31 #include <TLorentzRotation.h>
32 
33 // post-R3 this needs -I ${GENIE}/src/Framework as well as -I ${GENIE}/src
34 // #include "Conventions/GVersion.h"
35 // can't rely on that ... so we'll require -DGENIE_PRE_R3
36 
37 
38 #ifdef GENIE_PRE_R3
39  #include "Conventions/GVersion.h"
40  #include "EVGDrivers/GFluxI.h"
41  #include "PDG/PDGUtils.h"
42  #if __GENIE_RELEASE_CODE__ >= GRELCODE(2,9,0)
43  #include "FluxDrivers/GFluxExposureI.h"
44  #include "FluxDrivers/GFluxFileConfigI.h"
45  #endif
46 #else
52 #endif
53 
54 class TFile;
55 class TChain;
56 class TTree;
57 class TBranch;
58 
59 namespace bsim {
60  class Dk2Nu;
61  class DkMeta;
62  class NuChoice;
63 }
64 
65 namespace genie {
66 namespace flux {
67 
68 // #if __GENIE_RELEASE_CODE__ >= GRELCODE(2,9,0)
69 class GDk2NuFlux : public GFluxI,
70  public GFluxExposureI, public GFluxFileConfigI {
71 // #else
72 // class GDk2NuFlux: public GFluxI {
73 // #endif
74 
75 public :
76  GDk2NuFlux();
77  ~GDk2NuFlux();
78 
79  // Methods implementing the GENIE GFluxI interface, required for integrating
80  // the NuMI neutrino flux simulations with the GENIE event generation drivers
81 
82  const PDGCodeList & FluxParticles (void) { return *fPdgCList; }
83  double MaxEnergy (void) { return fMaxEv; }
84  bool GenerateNext (void);
85  int PdgCode (void);
86  double Weight (void) { return fWeight; }
87  const TLorentzVector & Momentum (void);
88  const TLorentzVector & Position (void);
89  bool End (void) { return fEnd; }
90  long int Index (void) { return fIEntry; }
91  void Clear (Option_t * opt);
92  void GenerateWeighted (bool gen_weighted);
93 
94  // Methods specific to this flux driver,
95  // for configuration/initialization of the flux & event generation drivers
96  // and and for passing-through flux information (e.g. neutrino parent decay
97  // kinematics) not used by the generator but required by analyses/processing
98  // further downstream
99 
100  //
101  // information about or actions on current entry
102  //
103  const bsim::NuChoice & GetNuChoice(void) { return *fCurNuChoice; };
104  const bsim::Dk2Nu & GetDk2Nu(void) { return *fCurDk2Nu; };
105  const bsim::DkMeta & GetDkMeta(void) { LoadDkMeta(); return *fCurDkMeta; };
106 
107  Long64_t GetEntryNumber() { return fIEntry; } ///< index in chain
108 
109  double GetDecayDist() const; ///< dist (user units) from dk to current pos
110  void MoveToZ0(double z0); ///< move ray origin to user coord Z0
111 
112  //
113  // information about the current state
114  //
115 // #if __GENIE_RELEASE_CODE__ >= GRELCODE(2,9,0)
116  virtual double GetTotalExposure() const; // GFluxExposureI interface
117 // #else
118 // // these are elminated with the introduction of GFluxFileConfigI
119 // void SetXMLFile(std::string xmlbasename="GNuMIFlux.xml")
120 // { fXMLbasename = xmlbasename; }
121 // std::string GetXMLFile() const { return fXMLbasename; }
122 // #endif
123 
124  double POT_curr(void); ///< current average POT (RWH?)
125  double UsedPOTs(void) const; ///< # of protons-on-target used
126  long int NFluxNeutrinos(void) const { return fNNeutrinos; } ///< number of flux neutrinos looped so far
127  double SumWeight(void) const { return fSumWeight; } ///< integrated weight for flux neutrinos looped so far
128 
129  void PrintCurrent(void); ///< print current entry from leaves
130  void PrintConfig(); ///< print the current configuration
131 
132  std::vector<std::string> GetFileList(); ///< list of files currently part of chain
133 
134  //
135  // GFluxFileConfigI interface
136  //
137  virtual void LoadBeamSimData(const std::vector<std::string>& filenames,
138  const std::string& det_loc);
139 // #if __GENIE_RELEASE_CODE__ >= GRELCODE(2,9,0)
140  using GFluxFileConfigI::LoadBeamSimData; // inherit the rest
141  virtual void GetBranchInfo(std::vector<std::string>& branchNames,
142  std::vector<std::string>& branchClassNames,
143  std::vector<void**>& branchObjPointers);
144  virtual TTree* GetMetaDataTree();
145 // #else
146 // void LoadBeamSimData(std::set<std::string> filenames, std::string det_loc); ///< load root flux ntuple files and config
147 // void LoadBeamSimData(std::string filename, std::string det_loc); ///< older (obsolete) single file version
148 //
149 // #endif
150 
151  //
152  // configuration of GDk2NuFlux
153  //
154  void SetTreeNames(std::string fname = "dk2nuTree",
155  std::string mname = "dkmetaTree")
156  { fTreeNames[0] = fname; fTreeNames[1] = mname; }
157 
158 
159  bool LoadConfig(std::string cfg); ///< load a named configuration
160 
161  void SetMaxEnergy(double Ev); ///< specify max flx nu energy
162 
163  void SetGenWeighted(bool genwgt=false) { fGenWeighted = genwgt; } ///< toggle whether GenerateNext() returns weight=1 flux (initial default false)
164 
165  void SetEntryReuse(long int nuse=1); ///< # of times to use entry before moving to next
166 
167  void ScanForMaxWeight(void); ///< scan for max flux weight (before generating unweighted flux neutrinos)
168  void SetMaxWgtScan(double fudge = 1.05, long int nentries = 2500000) ///< configuration when estimating max weight
169  { fMaxWgtFudge = fudge; fMaxWgtEntries = nentries; }
170  void SetMaxEFudge(double fudge = 1.05) ///< extra fudge factor in estimating maximum energy
171  { fMaxEFudge = fudge; }
172 
173  void SetMinMaxWeight(double minwgt) ///< user floor on estimating maxweight (to avoid bumping)
174  { fMinMaxWeight = minwgt; }
175  void SetMaxWeightFailModel(int i=0) ///< 0=bump, 1=leave frozen, 2=abort
176  { fMaxWgtFailModel = i; }
177 
178 
179  void SetApplyWindowTiltWeight(bool apply = true) ///< apply wgt due to tilt of flux window relative to beam
180  { fApplyTiltWeight = apply; }
181 
182  // GDk2NuFlux uses "cm" as the length unit consistently internally (this is
183  // the length units used by the dk2nu ntuples). User interactions
184  // setting the beam-to-detector coordinate transform, flux window, and the
185  // returned position might need to be in other units. Use:
186  // double scale = genie::utils::units::UnitFromString("cm");
187  // ( #include "Utils/UnitUtils.h for declaration )
188  // to get the correct scale factor to pass in. This should get set
189  // FIRST before setting detector position/rotation
190 
191  void SetLengthUnits(double user_units); ///< Set units assumed by user
192  double LengthUnits(void) const; ///< Return user units
193 
194  // set relative orientation of user coords vs. beam system, i.e.
195  // x3_user = ( beamrot * x3_beam ) + x0beam_user
196  // p3_user = beamrot * p3_beam
197 
198  ///< beam (0,0,0) relative to user frame, beam direction in user frame
199  void SetBeamRotation(TRotation beamrot);
200  void SetBeamCenter(TVector3 beam0);
201  TRotation GetBeamRotation() const; ///< rotation to apply from beam->user
202  TVector3 GetBeamCenter() const; ///< beam origin in user frame
203 
204  // configure a flux window (or point) where E_nu and weight are evaluated
205 
206  // rwh: potential upgrade: allow flux window set/get in beam coords
207  // as optional flag to *etFluxWindow
208  bool IsFluxSphere() const { return fIsSphere; } ///< flat window or sphere
209  void SetFluxWindow(TVector3 p1, TVector3 p2, TVector3 p3); ///< 3 points define a plane (by default in user coordinates)
210  void GetFluxWindow(TVector3& p1, TVector3& p2, TVector3& p3) const; ///< 3 points define a plane in beam coordinate
211  void SetFluxSphere(TVector3 center, double radius, bool inDetCoord=true); ///< specification of a sphere
212  void GetFluxSphere(TVector3& center, double& radius, bool inDetCoord=true) const; ///< specification of a sphere
213 
214 // #if __GENIE_RELEASE_CODE__ < GRELCODE(2,9,0)
215 // // migrated to GFluxFileConfigI
216 // void SetUpstreamZ(double z0);
217 // void SetNumOfCycles(long int ncycle);
218 // void SetFluxParticles(const PDGCodeList & particles);
219 // #endif
220 
221  //
222  // Actual coordinate transformations b=beam, u=user (e.g. detector)
223  //
224  void Beam2UserPos(const TLorentzVector& beamxyz,
225  TLorentzVector& usrxyz ) const;
226  void Beam2UserDir(const TLorentzVector& beamdir,
227  TLorentzVector& usrdir ) const;
228  void Beam2UserP4 (const TLorentzVector& beamp4,
229  TLorentzVector& usrp4 ) const;
230  void User2BeamPos(const TLorentzVector& usrxyz,
231  TLorentzVector& beamxyz ) const;
232  void User2BeamDir(const TLorentzVector& usrdir,
233  TLorentzVector& beamdir ) const;
234  void User2BeamP4 (const TLorentzVector& usrp4,
235  TLorentzVector& beamp4 ) const;
236 
237  TVector3 FluxWindowNormal() const { return (fIsSphere) ? SphereNormal() : fFluxWindowNormal; }
238  TVector3 SphereNormal() const;
239 
240 private:
241 
242  // Private methods
243  //
244  bool GenerateNext_weighted (void);
245  void Initialize (void);
246  void SetDefaults (void);
247  void CleanUp (void);
248  void ResetCurrent (void);
249  void AddFile (TTree* fluxtree, TTree* metatree, std::string fname);
250  void CalcEffPOTsPerNu (void);
251  void LoadDkMeta (void);
252 
253  // Private data members
254  //
255  double fMaxEv; ///< maximum energy
256 
257  bool fEnd; ///< end condition reached
258 
259 // #if __GENIE_RELEASE_CODE__ < GRELCODE(2,9,0)
260 // // incorporated into GFluxFileConfigI
261 // PDGCodeList * fPdgCList; ///< list of neutrino pdg-codes to generate
262 // PDGCodeList * fPdgCListRej; ///< list of neutrino pdg-codes seen but rejected
263 // std::string fXMLbasename; ///< XML filename for config data
264 // double fZ0; ///< configurable starting z position for each flux neutrino (in detector coord system)
265 // long int fNCycles; ///< # times to cycle through the flux ntuple
266 // long int fICycle; ///< current file cycle
267 // #endif
268 
269  std::vector<std::string> fNuFluxFilePatterns; ///< (potentially wildcarded) path(s)
270 
271  std::string fTreeNames[2]; ///< pair of names "dk2nuTree", "dkmetaTree"
272  TChain* fNuFluxTree; ///< TTree // REF ONLY!
273  TChain* fNuMetaTree; ///< TTree // REF ONLY!
274 
278 
279  int fNFiles; ///< number of files in chain
280  Long64_t fNEntries; ///< number of flux ntuple entries
281  Long64_t fIEntry; ///< current flux ntuple entry
282  Long64_t fNuTot; ///< cummulative # of entries (=fNEntries)
283  Long64_t fFilePOTs; ///< # of protons-on-target represented by all files
284 
285  std:: map<int,int> fJobToMetaIndex; ///< quick lookup from job# to meta chain
286 
287  double fWeight; ///< current neutrino weight, =1 if generating unweighted entries
288  double fMaxWeight; ///< max flux neutrino weight in input file
289  double fMinMaxWeight; ///< user set lower limit on estimate
290  double fMaxWeightScan; ///< initial estimate from scan
291  double fMaxWeightInit; ///< max of scan & minmaxweight
292  double fMaxWeightMax; ///< if "frozen" this is what bump would given
293  double fMaxWgtFudge; ///< fudge factor for estimating max wgt
294  long int fMaxWgtEntries; ///< # of entries in estimating max wgt
295  double fMaxEFudge; ///< fudge factor for estmating max enu (0=> use fixed 120GeV)
296  long int fMaxWgtExceeded; ///< track failures of estimate
297  int fMaxWgtFailModel; ///< what to do ... 0=bump, 1=frozen, 2=abort
298 
299  long int fNUse; ///< how often to use same entry in a row
300  long int fIUse; ///< current # of times an entry has been used
301  double fSumWeight; ///< sum of weights for nus thrown so far
302  long int fNNeutrinos; ///< number of flux neutrinos thrown so far
303  double fEffPOTsPerNu; ///< what a entry is worth ...
304  double fAccumPOTs; ///< POTs used so far
305 
306  bool fGenWeighted; ///< does GenerateNext() give weights?
307  bool fApplyTiltWeight; ///< wgt due to window normal not || beam
308  bool fDetLocIsSet; ///< is a flux location (near/far) set?
309 
310  double fLengthUnits; ///< units for coord in user exchanges
311  double fLengthScaleB2U; ///< scale factor beam (cm) --> user
312  double fLengthScaleU2B; ///< scale factor beam user --> (cm)
313 
314  TLorentzVector fBeamZero; ///< beam origin in user coords
315  TLorentzRotation fBeamRot; ///< rotation applied beam --> user coord
316  TLorentzRotation fBeamRotInv;
317 
318  bool fIsSphere; ///< doing this on a sphere rather than a flat window?
319  TVector3 fFluxWindowPtUser[3]; ///< user points of flux window
320  TLorentzVector fFluxWindowBase; ///< base point for flux window - beam coord
321  TLorentzVector fFluxWindowDir1; ///< extent for flux window (direction 1)
322  TLorentzVector fFluxWindowDir2; ///< extent for flux window (direction 2)
325  TVector3 fFluxWindowNormal; ///< normal direction for flux window -- beam coord
326  TVector3 fFluxSphereCenterUser; ///< center for flux sphere - user coords
327  TVector3 fFluxSphereCenterBeam; ///< center for flux sphere - beam coords
328  double fFluxSphereRadius; ///< radius for flux sphere
329 
330  TLorentzVector fgX4dkvtx; ///< decay 4-position beam coord
331 
332 };
333 
334 } // flux namespace
335 } // genie namespace
336 
337 #endif // _GDK2NUFLUX_H_
bool fApplyTiltWeight
wgt due to window normal not || beam
Definition: GDk2NuFlux.h:307
def apply(command)
Definition: g4zmq.py:38
void SetMaxWgtScan(double fudge=1.05, long int nentries=2500000)
Definition: GDk2NuFlux.h:168
bool IsFluxSphere() const
flat window or sphere
Definition: GDk2NuFlux.h:208
long int NFluxNeutrinos(void) const
number of flux neutrinos looped so far
Definition: GDk2NuFlux.h:126
void SetMaxWeightFailModel(int i=0)
Definition: GDk2NuFlux.h:175
bool fEnd
end condition reached
Definition: GDk2NuFlux.h:257
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
void SetMinMaxWeight(double minwgt)
Definition: GDk2NuFlux.h:173
TLorentzVector fgX4dkvtx
decay 4-position beam coord
Definition: GDk2NuFlux.h:330
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GDk2NuFlux.h:302
bsim::Dk2Nu * fCurDk2Nu
Definition: GDk2NuFlux.h:275
TChain * fNuFluxTree
TTree // REF ONLY!
Definition: GDk2NuFlux.h:272
def AddFile(fout, deets, f)
Long64_t fNuTot
cummulative # of entries (=fNEntries)
Definition: GDk2NuFlux.h:282
double MaxEnergy(void)
declare the max flux neutrino energy that can be generated (for init. purposes)
Definition: GDk2NuFlux.h:83
double fMaxWeightScan
initial estimate from scan
Definition: GDk2NuFlux.h:290
double fLengthScaleB2U
scale factor beam (cm) –> user
Definition: GDk2NuFlux.h:311
double fSumWeight
sum of weights for nus thrown so far
Definition: GDk2NuFlux.h:301
double fWeight
current neutrino weight, =1 if generating unweighted entries
Definition: GDk2NuFlux.h:287
const bsim::NuChoice & GetNuChoice(void)
Definition: GDk2NuFlux.h:103
bsim namespace for beam simulation classes and functions
Definition: GDk2NuFlux.h:59
double fMaxEv
maximum energy
Definition: GDk2NuFlux.h:255
double fMaxWgtFudge
fudge factor for estimating max wgt
Definition: GDk2NuFlux.h:293
bool fIsSphere
doing this on a sphere rather than a flat window?
Definition: GDk2NuFlux.h:318
long int fNUse
how often to use same entry in a row
Definition: GDk2NuFlux.h:299
TLorentzRotation fBeamRotInv
Definition: GDk2NuFlux.h:316
Long64_t fIEntry
current flux ntuple entry
Definition: GDk2NuFlux.h:281
TVector3 fFluxSphereCenterBeam
center for flux sphere - beam coords
Definition: GDk2NuFlux.h:327
double fEffPOTsPerNu
what a entry is worth ...
Definition: GDk2NuFlux.h:303
std::vector< std::string > fNuFluxFilePatterns
(potentially wildcarded) path(s)
Definition: GDk2NuFlux.h:269
double fFluxSphereRadius
radius for flux sphere
Definition: GDk2NuFlux.h:328
Long64_t fFilePOTs
of protons-on-target represented by all files
Definition: GDk2NuFlux.h:283
void Initialize(Bool_t useTMVAStyle=kTRUE)
Definition: tmvaglob.C:184
Loaders::FluxType flux
long int fIUse
current # of times an entry has been used
Definition: GDk2NuFlux.h:300
A list of PDG codes.
Definition: PDGCodeList.h:33
TChain * fNuMetaTree
TTree // REF ONLY!
Definition: GDk2NuFlux.h:273
double fMaxWeight
max flux neutrino weight in input file
Definition: GDk2NuFlux.h:288
void SetTreeNames(std::string fname="dk2nuTree", std::string mname="dkmetaTree")
Definition: GDk2NuFlux.h:154
void SetMaxEFudge(double fudge=1.05)
Definition: GDk2NuFlux.h:170
double fMaxWeightInit
max of scan & minmaxweight
Definition: GDk2NuFlux.h:291
void SetApplyWindowTiltWeight(bool apply=true)
Definition: GDk2NuFlux.h:179
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
bool End(void)
true if no more flux nu&#39;s can be thrown (eg reaching end of beam sim ntuples)
Definition: GDk2NuFlux.h:89
Long64_t nentries
std::map< int, int > fJobToMetaIndex
quick lookup from job# to meta chain
Definition: GDk2NuFlux.h:285
TLorentzVector fBeamZero
beam origin in user coords
Definition: GDk2NuFlux.h:314
double fMaxWeightMax
if "frozen" this is what bump would given
Definition: GDk2NuFlux.h:292
GENIE interface for uniform flux exposure iterface.
double fLengthScaleU2B
scale factor beam user –> (cm)
Definition: GDk2NuFlux.h:312
int fMaxWgtFailModel
what to do ... 0=bump, 1=frozen, 2=abort
Definition: GDk2NuFlux.h:297
Double_t radius
void SetGenWeighted(bool genwgt=false)
toggle whether GenerateNext() returns weight=1 flux (initial default false)
Definition: GDk2NuFlux.h:163
double SumWeight(void) const
integrated weight for flux neutrinos looped so far
Definition: GDk2NuFlux.h:127
Long64_t GetEntryNumber()
index in chain
Definition: GDk2NuFlux.h:107
const PDGCodeList & FluxParticles(void)
declare list of flux neutrinos that can be generated (for init. purposes)
Definition: GDk2NuFlux.h:82
bsim::DkMeta * fCurDkMeta
Definition: GDk2NuFlux.h:276
double fLengthUnits
units for coord in user exchanges
Definition: GDk2NuFlux.h:310
bsim::NuChoice * fCurNuChoice
Definition: GDk2NuFlux.h:277
GENIE interface for uniform flux exposure iterface.
TLorentzVector fFluxWindowDir1
extent for flux window (direction 1)
Definition: GDk2NuFlux.h:321
const bsim::DkMeta & GetDkMeta(void)
Definition: GDk2NuFlux.h:105
int fNFiles
number of files in chain
Definition: GDk2NuFlux.h:279
long int fMaxWgtEntries
of entries in estimating max wgt
Definition: GDk2NuFlux.h:294
double fMaxEFudge
fudge factor for estmating max enu (0=> use fixed 120GeV)
Definition: GDk2NuFlux.h:295
TVector3 fFluxWindowNormal
normal direction for flux window – beam coord
Definition: GDk2NuFlux.h:325
TLorentzVector fFluxWindowBase
base point for flux window - beam coord
Definition: GDk2NuFlux.h:320
double Weight(void)
returns the flux neutrino weight (if any)
Definition: GDk2NuFlux.h:86
An implementation of the GENIE GFluxI interface ("flux driver") encapsulating reading/processing the ...
Definition: GDk2NuFlux.h:69
double fAccumPOTs
POTs used so far.
Definition: GDk2NuFlux.h:304
TVector3 FluxWindowNormal() const
Definition: GDk2NuFlux.h:237
TVector3 fFluxSphereCenterUser
center for flux sphere - user coords
Definition: GDk2NuFlux.h:326
TLorentzVector fFluxWindowDir2
extent for flux window (direction 2)
Definition: GDk2NuFlux.h:322
bool fDetLocIsSet
is a flux location (near/far) set?
Definition: GDk2NuFlux.h:308
bool fGenWeighted
does GenerateNext() give weights?
Definition: GDk2NuFlux.h:306
Long64_t fNEntries
number of flux ntuple entries
Definition: GDk2NuFlux.h:280
long int fMaxWgtExceeded
track failures of estimate
Definition: GDk2NuFlux.h:296
TLorentzRotation fBeamRot
rotation applied beam –> user coord
Definition: GDk2NuFlux.h:315
double fMinMaxWeight
user set lower limit on estimate
Definition: GDk2NuFlux.h:289
const bsim::Dk2Nu & GetDk2Nu(void)
Definition: GDk2NuFlux.h:104
long int Index(void)
returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current ent...
Definition: GDk2NuFlux.h:90
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
enum BeamMode string