dk2nu.h
Go to the documentation of this file.
1 ///----------------------------------------------------------------------------
2 /**
3  * \class bsim::Dk2Nu
4  * \file dk2nu.h
5  *
6  * \brief A class that defines the "dk2nu" object used as the primary
7  * branch for a TTree for the output of neutrino flux simulations
8  * such as g4numi, g4numi_flugg, etc.
9  *
10  * \author (last to touch it) $Author: rhatcher $
11  *
12  * \version $Revision: 1.3 $
13  *
14  * \date $Date: 2012-11-21 04:44:50 $
15  *
16  * Contact: rhatcher@fnal.gov
17  *
18  * $Id: dk2nu.h,v 1.3 2012-11-21 04:44:50 rhatcher Exp $
19  *
20  * Notes tagged with "DK2NU" are questions that should be answered
21  */
22 ///----------------------------------------------------------------------------
23 
24 #ifndef BSIM_DK2NU_H
25 #define BSIM_DK2NU_H
26 
27 #include "TROOT.h"
28 #include "TObject.h"
29 
30 #include <vector>
31 #include <string>
32 
33 #define DK2NUVER 11 // KEEP THIS UP-TO-DATE! increment for each change
34 #define KEEP_ANCESTOR_PPRODPXYZ 1
35 
36 namespace bsim {
37  /**
38  * All the data members are public as these classes are used as
39  * generalized structs. As they will be branches of a TTree no
40  * specialized naming indicators signifying that they are member data
41  * of a class will be used, nor will any fancy capitalization schemes.
42  *
43  * All classes must implement a clear() method that resets their values
44  * to an identifiably invalid state or clears any vectors. Additionally
45  * classes should provide a AsString() method for formatting themselves
46  * for use output.
47  */
48 
49  /**
50  *=======================================================================
51  * Specialized enumerations
52  */
53 
54  /**
55  * Proposed flag bits:
56  */
57  typedef enum flgbitval {
58  kFlgOverflow = 0x00000001,
59  kMaskReserved = 0x0000FFFF,
60  kMaskUser = 0xFFFF0000
61  } flgbitval_t;
62 
63  /**
64  * Enumeration of decay processes, stored in "ndecay"
65  * stored as integer; these are for reference
66  * DK2NU: should there be an associated AsString() method
67  * that returns a text (optionally formatted for latex?)?
68  */
69  typedef enum dkproc {
71  dkp_k0l_nuepimep = 1, ///< k0long => nu_e + pi- + e+
72  dkp_k0l_nuebpipem = 2, ///< k0long => nu_e_bar + p+ + e-
73  dkp_k0l_numupimmup = 3, ///< k0long => nu_mu + pi- + mu+
74  dkp_k0l_numubpipmum = 4, ///< k0long => nu_mu_bar + pi+ + mu-
75  dkp_kp_numumup = 5, ///< k+ => nu_mu + mu+
76  dkp_kp_nuepi0ep = 6, ///< k+ => nu_e + pi0 + e+
77  dkp_kp_numupi0mup = 7, ///< k+ => nu_mu + pi0 + mu+
78  dkp_kp_numubmum = 8, ///< k- => nu_mu_bar + mu-
79  dkp_kp_nuebpi0em = 9, ///< k- => nu_e_bar + pi0 + e-
80  dkp_kp_numubpi0mum = 10, ///< k- => nu_mu_bar + pi0 + mu-
81  dkp_mup_nusep = 11, ///< mu+ => nu_mu_bar + nu_e + e+
82  dkp_mum_nusep = 12, ///< mu- => nu_mu + nu_e_bar + e-
83  dk_pip_numumup = 13, ///< pi+ => nu_mu + mu+
84  dk_pim_numubmum = 14, ///< pi- => nu_mu_bar + mu-
85  dk_mum_capture = 15, ///< mu- + (A,Z) => nu
86  /// (different kinematics from dkp_mum_nusep)
87  dkp_maximum, ///< one-beyond end for iterating
88  dkp_other = 999 ///< flag for unusual cases
89  } dkproc_t;
90 
91  ///---------------------------------------------------------------------------
92  /**
93  *============================================================================
94  * Fixed Decays through a specific point:
95  */
96  class NuRay
97  {
98  public:
99  Double_t px; ///< px for nu at location
100  Double_t py; ///< py for nu at location
101  Double_t pz; ///< pz for nu at location
102  Double_t E; ///< E for nu at location
103  Double_t wgt; ///< weight for nu at location
104 
105  public:
106  NuRay();
107  NuRay(double px, double py, double pz, double E, double wgt);
108  virtual ~NuRay();
109  void clear(const std::string &opt = ""); ///< reset everything
110  std::string AsString(const std::string& opt = "") const;
111 
112  private:
113  ClassDef(bsim::NuRay,DK2NUVER)
114  }; // end-of-class bsim::NuRay
115 
116  ///---------------------------------------------------------------------------
117  /**
118  *============================================================================
119  * Decay Data:
120  * Core information about the neutrino and the decay that gave rise to it.
121  * % = necessary for reweighting
122  */
123  class Decay
124  {
125  public:
126  Int_t norig; ///< not used?
127  Int_t ndecay; ///< decay process (see dkproc_t)
128  Int_t ntype; ///< % neutrino flavor (PDG? code)
129 
130  Double_t vx; ///< % neutrino production vertex x
131  Double_t vy; ///< % neutrino production vertex y
132  Double_t vz; ///< % neutrino production vertex z
133  Double_t pdpx; ///< % px momentum of nu parent at (vx,vy,vz)
134  Double_t pdpy; ///< % py momentum of nu parent at (vx,vy,vz)
135  Double_t pdpz; ///< % pz momentum of nu parent at (vx,vy,vz)
136 
137  /** these are used in muon decay case? */
138  Double_t ppdxdz; ///< % direction of nu parent at its production point
139  Double_t ppdydz; ///< % direction of nu parent at its production point
140  Double_t pppz; ///< % z momentum of nu parent at its production point
141  Double_t ppenergy; ///< % energy of nu parent at its production point
142 
143  Int_t ppmedium; ///< material nu parent was produced in
144  Int_t ptype; ///< % nu parent species (PDG? code)
145 
146  /** momentum and energy of nu grandparent at
147  muons: grandparent decay point
148  hadrons: grandparent production point
149  Huh? this needs better documentation
150  Marco DT says: One should look at parent type ptype.
151  If ptype is a muon, then muparpx,y,z,e
152  are momentum and energy of the neutrino
153  grandparent (muon parent) at its decay point.
154  Otherwise (for all other values of ptype),
155  muparpx,y,z,e refer to neutrino grandparent
156  (a hadron in this case) production point.
157  If the Ancestor List is on , then
158  these variables are superfluous.
159  */
160  Double_t muparpx; ///< %
161  Double_t muparpy; ///< %
162  Double_t muparpz; ///< %
163  Double_t mupare; ///< % energy of nu grandparent
164 
165  Double_t necm; ///< % nu energy in center-of-mass frame
166  Double_t nimpwt; ///< % cumulative importance weight prod to decay
167 
168  public:
169  Decay();
170  virtual ~Decay();
171  void clear(const std::string &opt = ""); ///< reset everything
172  std::string AsString(const std::string& opt = "") const;
173 
174  private:
175  ClassDef(bsim::Decay,DK2NUVER)
176  }; // end-of-class bsim::Decay
177 
178  ///---------------------------------------------------------------------------
179  /**
180  *============================================================================
181  * Ancestor Data:
182  * Information about the chain of particles from the initial proton (indx=0)
183  * to the final neutrino.
184  */
185  class Ancestor
186  {
187  public:
188  Int_t pdg; ///< ancestor species
189 
190  Double_t startx; ///< particle x initial position
191  Double_t starty; ///< particle y initial position
192  Double_t startz; ///< particle z initial position
193  Double_t startt; ///< particle initial time
194 
195  Double_t startpx; ///< particle x initial momentum
196  Double_t startpy; ///< particle y initial momentum
197  Double_t startpz; ///< particle z initial momentum
198 
199  Double_t stoppx; ///< particle x final momentum
200  Double_t stoppy; ///< particle y final momentum
201  Double_t stoppz; ///< particle z final momentum
202 
203  Double_t polx; ///< x component of polarization
204  Double_t poly; ///< y component of polarization
205  Double_t polz; ///< z component of polarization
206 
207  // what are these ... somehow different from stoppx[-1]?
208  // Marco DT says: Yes, they can be different. Nu parent first entry can
209  // be different than [-1]. Exapmle: it could be [-2], then [-1]
210  // contains a parent elastic interaction.
211  // I'm gettint rid of these variables and I'll add a
212  // parIndex instead, see following.
213 #ifdef KEEP_ANCESTOR_PPRODPXYZ
214  Double_t pprodpx; ///< parent x momentum when producing this particle
215  Double_t pprodpy; ///< parent y momentum when producing this particle
216  Double_t pprodpz; ///< parent z momentum when producing this particle
217 #endif
218 
219  Int_t nucleus; ///< nucleus (PDG) type causing the scatter
220 
221  // Marco DT is adding parIndex
222  Int_t parIndex; ///< particle index, from nu (0), parent (1) ... to proton (n)
223 
224  std::string proc; ///< name of the process that creates this particle
225  std::string ivol; ///< name of the volume where the particle starts
226  std::string imat; ///< name of the material where the particle starts
227 
228  public:
229  Ancestor();
230  virtual ~Ancestor();
231  void clear(const std::string &opt = ""); ///< reset everything
232  std::string AsString(const std::string& opt = "") const;
233 
234  /// set triplets
235  void SetStartXYZT(Double_t x, Double_t y, Double_t z, Double_t t);
236  void SetStartP(Double_t px, Double_t py, Double_t pz);
237  void SetStopP(Double_t px, Double_t py, Double_t pz);
238  void SetPProdP(Double_t px, Double_t py, Double_t pz);
239 
240  /// helper functions
241  Double_t r() const;
242  Double_t startpt() const;
243  Double_t startp() const;
244  Double_t stoppt() const;
245  Double_t stopp() const;
246 
247  // helper functions added by Marco
248  bool IsInTarget();
249 
250 #ifdef KEEP_ANCESTOR_PPRODPXYZ
251  Double_t pprodpt() const;
252  Double_t pprodp() const;
253 #endif
254 
255  private:
256  ClassDef(bsim::Ancestor,DK2NUVER)
257  }; // end-of-class bsim::Ancestor
258 
259  ///---------------------------------------------------------------------------
260  /**
261  *============================================================================
262  * these ancestors are possibly, but not necessarily, the direct nu parent
263  * DK2NU: can these be removed in favor of cascade info (ancestor above)?
264  * 2012-11-08: for now keep these
265  */
266  class TgtExit
267  {
268  public:
269  Double_t tvx; ///< x position of nu ancestor as it exits target
270  Double_t tvy; ///< y position of nu ancestor as it exits target
271  Double_t tvz; ///< z position of nu ancestor as it exits target
272  Double_t tpx; ///< x momentum of nu ancestor as it exits target
273  Double_t tpy; ///< y momentum of nu ancestor as it exits target
274  Double_t tpz; ///< z momentum of nu ancestor as it exits target
275  Int_t tptype; ///< species of ancestor exiting the target
276  Int_t tgen; ///< nu parent generation in cascade:
277  ///< 1=primary proton
278  ///< 2=particles produced by proton interaction
279  ///< etc
280  public:
281  TgtExit();
282  virtual ~TgtExit();
283  void clear(const std::string &opt = ""); ///< reset everything
284  std::string AsString(const std::string& opt = "") const;
285 
286  private:
287  ClassDef(bsim::TgtExit,DK2NUVER)
288  }; // end-of-class TgtExit
289 
290  ///---------------------------------------------------------------------------
291  /**
292  *============================================================================
293  * Track points stored at special locations for plotting trajectories
294  */
295  class Traj
296  {
297  public:
298  Double_t trkx;
299  Double_t trky;
300  Double_t trkz;
301  Double_t trkpx;
302  Double_t trkpy;
303  Double_t trkpz;
304 
305  public:
306  Traj();
307  virtual ~Traj();
308  void clear(const std::string &opt = ""); ///< reset everything
309  std::string AsString(const std::string& opt = "") const;
310 
311  private:
312  ClassDef(bsim::Traj,DK2NUVER)
313  }; // end-of-class bsim::Traj
314 
315  ///---------------------------------------------------------------------------
316  /**
317  *============================================================================
318  * This is the structure that is the basis for the flux ntuple
319  */
320  class Dk2Nu
321  {
322  public:
323  Int_t job; ///< identifying job #
324  Int_t potnum; ///< proton # processed by simulation
325  Int_t jobindx; ///< unique # within the job
326  ///< potnum isn't sufficient
327  bsim::Decay decay; ///< basic decay information
328  std::vector<bsim::NuRay> nuray; ///< rays through detector fixed points
329  std::vector<bsim::Ancestor> ancestor; ///< chain from proton to neutrino
330 
331  /**
332  * These are ancestor.vx[size-2] kept, for now, for convenience
333  */
334  Double_t ppvx; ///< production vertex x of nu parent
335  Double_t ppvy; ///< production vertex y of nu parent
336  Double_t ppvz; ///< production vertex z of nu parent
337 
338  bsim::TgtExit tgtexit; ///< info about leaving the target
339  std::vector<bsim::Traj> traj; ///< [optional] trajectory points
340 
341  /**
342  *=======================================================================
343  * Special Info:
344  */
345  Int_t flagbits; ///< bits signify non-std setting such as
346  ///< Geant vs. PDG codes, mm vs. cm, Mev vs. GeV
347  ///< bit00 = ancestore overflow (carry over from g4minerva)
348  std::vector<Int_t> vint; ///< user defined vector of integers
349  std::vector<Double_t> vdbl; ///< user defined vector of doubles
350 
351  public:
352  /**
353  * Public methods for constructing/destruction and resetting the data
354  */
355  Dk2Nu();
356  virtual ~Dk2Nu();
357  void clear(const std::string &opt = ""); ///< reset everything
358  std::string AsString(const std::string& opt = "") const;
359  void Print(Option_t* option = "") const;
360 
361  size_t indxnu() const; ///< ancestor index of nu ancestor.size()-1
362  size_t indxp() const; ///< ancestor index of parent ancestor.size()-2
363  size_t indxgp() const; ///< ancestor index of grandparent ancestor.size()-3
364  bool overflow() const; ///< ancestor list is incomplete (g4 minerva overflow)
365 
366  private:
367  ClassDef(bsim::Dk2Nu,DK2NUVER)
368 
369  };
370 
371 } // end-of-namespace "bsim"
372 
373 // not part of namespace bsim
374 std::ostream& operator<<(std::ostream& os, const bsim::Dk2Nu& dk2nu);
375 std::ostream& operator<<(std::ostream& os, const bsim::NuRay& nuray);
376 std::ostream& operator<<(std::ostream& os, const bsim::Decay& decay);
377 std::ostream& operator<<(std::ostream& os, const bsim::Ancestor& ancestor);
378 std::ostream& operator<<(std::ostream& os, const bsim::TgtExit& tgtexit);
379 std::ostream& operator<<(std::ostream& os, const bsim::Traj& traj);
380 
381 #endif // BSIM_DK2NU_H
k0long => nu_mu_bar + pi+ + mu-
Definition: dk2nu.h:74
Double_t ppdydz
% direction of nu parent at its production point
Definition: dk2nu.h:139
Double_t muparpx
%
Definition: dk2nu.h:160
Int_t job
identifying job #
Definition: dk2nu.h:323
Double_t E
E for nu at location.
Definition: dk2nu.h:102
Double_t polx
x component of polarization
Definition: dk2nu.h:203
Double_t px
px for nu at location
Definition: dk2nu.h:99
Int_t parIndex
particle index, from nu (0), parent (1) ... to proton (n)
Definition: dk2nu.h:222
Double_t pppz
% z momentum of nu parent at its production point
Definition: dk2nu.h:140
Double_t pprodpz
parent z momentum when producing this particle
Definition: dk2nu.h:216
Double_t wgt
weight for nu at location
Definition: dk2nu.h:103
std::string ivol
name of the volume where the particle starts
Definition: dk2nu.h:225
Double_t py
py for nu at location
Definition: dk2nu.h:100
Double_t pdpx
% px momentum of nu parent at (vx,vy,vz)
Definition: dk2nu.h:133
std::vector< bsim::NuRay > nuray
rays through detector fixed points
Definition: dk2nu.h:328
Int_t potnum
proton # processed by simulation
Definition: dk2nu.h:324
std::string proc
name of the process that creates this particle
Definition: dk2nu.h:224
Double_t pprodpx
parent x momentum when producing this particle
Definition: dk2nu.h:214
flgbitval
Definition: dk2nu.h:57
Double_t poly
y component of polarization
Definition: dk2nu.h:204
#define DK2NUVER
Definition: dk2nu.h:33
Int_t tgen
Definition: dk2nu.h:276
void clear(const std::string &opt="")
reset everything
Double_t trkpx
Definition: dk2nu.h:301
Double_t nimpwt
% cumulative importance weight prod to decay
Definition: dk2nu.h:166
Double_t necm
% nu energy in center-of-mass frame
Definition: dk2nu.h:165
bsim::Decay decay
basic decay information
Definition: dk2nu.h:327
Double_t ppvz
production vertex z of nu parent
Definition: dk2nu.h:336
void Decay(const Decayer *decayer, int pdgc, double E, int ndecays)
Definition: gtestDecay.cxx:179
bsim namespace for beam simulation classes and functions
Definition: GDk2NuFlux.h:59
Double_t tvy
y position of nu ancestor as it exits target
Definition: dk2nu.h:270
bsim::Dk2Nu * dk2nu
Double_t ppdxdz
% direction of nu parent at its production point
Definition: dk2nu.h:138
Double_t ppvy
production vertex y of nu parent
Definition: dk2nu.h:335
mu- => nu_mu + nu_e_bar + e-
Definition: dk2nu.h:82
Double_t vz
% neutrino production vertex z
Definition: dk2nu.h:132
Double_t trkz
Definition: dk2nu.h:300
bsim::TgtExit tgtexit
info about leaving the target
Definition: dk2nu.h:338
Double_t startpz
particle z initial momentum
Definition: dk2nu.h:197
Int_t tptype
species of ancestor exiting the target
Definition: dk2nu.h:275
Double_t trkx
Definition: dk2nu.h:298
Int_t pdg
ancestor species
Definition: dk2nu.h:188
virtual ~NuRay()
Double_t pdpz
% pz momentum of nu parent at (vx,vy,vz)
Definition: dk2nu.h:135
Double_t muparpy
%
Definition: dk2nu.h:161
k0long => nu_mu + pi- + mu+
Definition: dk2nu.h:73
std::vector< Double_t > vdbl
user defined vector of doubles
Definition: dk2nu.h:349
k+ => nu_e + pi0 + e+
Definition: dk2nu.h:76
Double_t ppvx
production vertex x of nu parent
Definition: dk2nu.h:334
Double_t muparpz
%
Definition: dk2nu.h:162
Double_t tpz
z momentum of nu ancestor as it exits target
Definition: dk2nu.h:274
Double_t starty
particle y initial position
Definition: dk2nu.h:191
enum bsim::dkproc dkproc_t
Double_t ppenergy
% energy of nu parent at its production point
Definition: dk2nu.h:141
Double_t stoppz
particle z final momentum
Definition: dk2nu.h:201
Double_t tpy
y momentum of nu ancestor as it exits target
Definition: dk2nu.h:273
Int_t ntype
% neutrino flavor (PDG? code)
Definition: dk2nu.h:128
Double_t startpx
particle x initial momentum
Definition: dk2nu.h:195
Double_t stoppx
particle x final momentum
Definition: dk2nu.h:199
enum bsim::flgbitval flgbitval_t
std::string imat
name of the material where the particle starts
Definition: dk2nu.h:226
Double_t trkpy
Definition: dk2nu.h:302
one-beyond end for iterating
Definition: dk2nu.h:87
Int_t nucleus
nucleus (PDG) type causing the scatter
Definition: dk2nu.h:219
Double_t polz
z component of polarization
Definition: dk2nu.h:205
Double_t startx
particle x initial position
Definition: dk2nu.h:190
std::vector< Int_t > vint
user defined vector of integers
Definition: dk2nu.h:348
k+ => nu_mu + pi0 + mu+
Definition: dk2nu.h:77
Double_t startt
particle initial time
Definition: dk2nu.h:193
z
Definition: test.py:28
k- => nu_e_bar + pi0 + e-
Definition: dk2nu.h:79
k0long => nu_e + pi- + e+
Definition: dk2nu.h:71
void Print(std::string prefix, std::string name, std::string suffix="")
Definition: nue_pid_effs.C:68
Double_t trkpz
Definition: dk2nu.h:303
Double_t pz
pz for nu at location
Definition: dk2nu.h:101
k- => nu_mu_bar + mu-
Definition: dk2nu.h:78
dkproc
Definition: dk2nu.h:69
Double_t tvz
z position of nu ancestor as it exits target
Definition: dk2nu.h:271
Double_t stoppy
particle y final momentum
Definition: dk2nu.h:200
Definition: decay.py:1
Double_t startpy
particle y initial momentum
Definition: dk2nu.h:196
std::vector< bsim::Traj > traj
[optional] trajectory points
Definition: dk2nu.h:339
std::string AsString(const std::string &opt="") const
mu+ => nu_mu_bar + nu_e + e+
Definition: dk2nu.h:81
Double_t tpx
x momentum of nu ancestor as it exits target
Definition: dk2nu.h:272
k+ => nu_mu + mu+
Definition: dk2nu.h:75
flag for unusual cases
Definition: dk2nu.h:88
Double_t vx
% neutrino production vertex x
Definition: dk2nu.h:130
TRandom3 r(0)
k- => nu_mu_bar + pi0 + mu-
Definition: dk2nu.h:80
pi- => nu_mu_bar + mu-
Definition: dk2nu.h:84
Int_t ndecay
decay process (see dkproc_t)
Definition: dk2nu.h:127
Int_t jobindx
Definition: dk2nu.h:325
Double_t mupare
% energy of nu grandparent
Definition: dk2nu.h:163
Double_t vy
% neutrino production vertex y
Definition: dk2nu.h:131
Int_t flagbits
Definition: dk2nu.h:345
k0long => nu_e_bar + p+ + e-
Definition: dk2nu.h:72
Int_t norig
not used?
Definition: dk2nu.h:126
Int_t ppmedium
material nu parent was produced in
Definition: dk2nu.h:143
Double_t trky
Definition: dk2nu.h:299
pi+ => nu_mu + mu+
Definition: dk2nu.h:83
Double_t tvx
x position of nu ancestor as it exits target
Definition: dk2nu.h:269
std::vector< bsim::Ancestor > ancestor
chain from proton to neutrino
Definition: dk2nu.h:329
Double_t pdpy
% py momentum of nu parent at (vx,vy,vz)
Definition: dk2nu.h:134
std::ostream & operator<<(std::ostream &os, const bsim::Dk2Nu &dk2nu)
Definition: nuray.py:1
Double_t startz
particle z initial position
Definition: dk2nu.h:192
Int_t ptype
% nu parent species (PDG? code)
Definition: dk2nu.h:144
Double_t pprodpy
parent y momentum when producing this particle
Definition: dk2nu.h:215
enum BeamMode string