Classes | Typedefs | Enumerations | Functions | Variables
bsim Namespace Reference

bsim namespace for beam simulation classes and functions More...

Classes

class  Ancestor
 
More...
 
class  Decay
 
More...
 
class  Dk2Nu
 
More...
 
class  DkMeta
 
More...
 
class  Location
 
More...
 
class  NuChoice
 
More...
 
class  NuRay
 
More...
 
class  TgtExit
 
More...
 
class  Traj
 
More...
 

Typedefs

typedef enum bsim::flgbitval flgbitval_t
 
typedef enum bsim::dkproc dkproc_t
 

Enumerations

enum  flgbitval {
  kFlgOverflow = 0x00000001, kMaskReserved = 0x0000FFFF, kMaskUser = 0xFFFF0000, kFlgOverflow = 0x00000001,
  kMaskReserved = 0x0000FFFF, kMaskUser = 0xFFFF0000
}
 
enum  dkproc {
  dkp_unknown = 0, dkp_k0l_nuepimep = 1, dkp_k0l_nuebpipem = 2, dkp_k0l_numupimmup = 3,
  dkp_k0l_numubpipmum = 4, dkp_kp_numumup = 5, dkp_kp_nuepi0ep = 6, dkp_kp_numupi0mup = 7,
  dkp_kp_numubmum = 8, dkp_kp_nuebpi0em = 9, dkp_kp_numubpi0mum = 10, dkp_mup_nusep = 11,
  dkp_mum_nusep = 12, dk_pip_numumup = 13, dk_pim_numubmum = 14, dk_mum_capture = 15,
  dkp_maximum, dkp_other = 999, dkp_unknown = 0, dkp_k0l_nuepimep = 1,
  dkp_k0l_nuebpipem = 2, dkp_k0l_numupimmup = 3, dkp_k0l_numubpipmum = 4, dkp_kp_numumup = 5,
  dkp_kp_nuepi0ep = 6, dkp_kp_numupi0mup = 7, dkp_kp_numubmum = 8, dkp_kp_nuebpi0em = 9,
  dkp_kp_numubpi0mum = 10, dkp_mup_nusep = 11, dkp_mum_nusep = 12, dk_pip_numumup = 13,
  dk_pim_numubmum = 14, dk_mum_capture = 15, dkp_maximum, dkp_other = 999
}
 
enum  flgbitval {
  kFlgOverflow = 0x00000001, kMaskReserved = 0x0000FFFF, kMaskUser = 0xFFFF0000, kFlgOverflow = 0x00000001,
  kMaskReserved = 0x0000FFFF, kMaskUser = 0xFFFF0000
}
 
enum  dkproc {
  dkp_unknown = 0, dkp_k0l_nuepimep = 1, dkp_k0l_nuebpipem = 2, dkp_k0l_numupimmup = 3,
  dkp_k0l_numubpipmum = 4, dkp_kp_numumup = 5, dkp_kp_nuepi0ep = 6, dkp_kp_numupi0mup = 7,
  dkp_kp_numubmum = 8, dkp_kp_nuebpi0em = 9, dkp_kp_numubpi0mum = 10, dkp_mup_nusep = 11,
  dkp_mum_nusep = 12, dk_pip_numumup = 13, dk_pim_numubmum = 14, dk_mum_capture = 15,
  dkp_maximum, dkp_other = 999, dkp_unknown = 0, dkp_k0l_nuepimep = 1,
  dkp_k0l_nuebpipem = 2, dkp_k0l_numupimmup = 3, dkp_k0l_numubpipmum = 4, dkp_kp_numumup = 5,
  dkp_kp_nuepi0ep = 6, dkp_kp_numupi0mup = 7, dkp_kp_numubmum = 8, dkp_kp_nuebpi0em = 9,
  dkp_kp_numubpi0mum = 10, dkp_mup_nusep = 11, dkp_mum_nusep = 12, dk_pip_numumup = 13,
  dk_pim_numubmum = 14, dk_mum_capture = 15, dkp_maximum, dkp_other = 999
}
 

Functions

int calcEnuWgt (const bsim::Decay &decay, const TVector3 &xyz, double &enu, double &wgt_xy)
 workhorse routine More...
 
int calcEnuWgt (const bsim::Dk2Nu *dk2nu, const TVector3 &xyz, double &enu, double &wgt_xy)
 alternate interface More...
 
void calcLocationWeights (const bsim::DkMeta *dkmeta, bsim::Dk2Nu *dk2nu)
 user interface More...
 
Bool_t IsDefault (Float_t intval)
 
Bool_t IsDefault (Int_t intval)
 
Bool_t IsDefault (UInt_t intval)
 
Bool_t IsDefault (Double_t intval)
 
Bool_t IsDefault (Bool_t intval)
 
void readWeightLocations (std::string locfilename, std::vector< bsim::Location > &locations)
 
void readWeightLocations (std::string locfilename, bsim::DkMeta *dkmeta)
 a variant that will fill the dkmeta object More...
 
void printWeightLocations (std::vector< bsim::Location > &locations)
 print the locations More...
 
void printWeightLocations (bsim::DkMeta *dkmeta)
 a variant that prints the locations in the dkmeta object More...
 

Variables

static const Float_t kDfltFloat = -9999.99
 
static const Int_t kDfltInt = -9999
 
static const UInt_t kDfltUInt = 9999
 
static const Double_t kDfltDouble = -9999.99
 
static const Bool_t kDfltBool = false
 
static const Float_t kDfltFloat = -9999.99
 
static const Int_t kDfltInt = -9999
 
static const UInt_t kDfltUInt = 9999
 
static const Double_t kDfltDouble = -9999.99
 
static const Bool_t kDfltBool = false
 

Detailed Description

bsim namespace for beam simulation classes and functions


Typedef Documentation

Enumeration of decay processes, stored in "ndecay" stored as integer; these are for reference DK2NU: should there be an associated AsString() method that returns a text (optionally formatted for latex?)?

All the data members are public as these classes are used as generalized structs. As they will be branches of a TTree no specialized naming indicators signifying that they are member data of a class will be used, nor will any fancy capitalization schemes.

All classes must implement a clear() method that resets their values to an identifiably invalid state or clears any vectors. Additionally classes should provide a AsString() method for formatting themselves

for use output.

Specialized enumerations Proposed flag bits:

Enumeration Type Documentation

Enumeration of decay processes, stored in "ndecay" stored as integer; these are for reference DK2NU: should there be an associated AsString() method that returns a text (optionally formatted for latex?)?

Enumerator
dkp_unknown 
dkp_k0l_nuepimep 

k0long => nu_e + pi- + e+

dkp_k0l_nuebpipem 

k0long => nu_e_bar + p+ + e-

dkp_k0l_numupimmup 

k0long => nu_mu + pi- + mu+

dkp_k0l_numubpipmum 

k0long => nu_mu_bar + pi+ + mu-

dkp_kp_numumup 

k+ => nu_mu + mu+

dkp_kp_nuepi0ep 

k+ => nu_e + pi0 + e+

dkp_kp_numupi0mup 

k+ => nu_mu + pi0 + mu+

dkp_kp_numubmum 

k- => nu_mu_bar + mu-

dkp_kp_nuebpi0em 

k- => nu_e_bar + pi0 + e-

dkp_kp_numubpi0mum 

k- => nu_mu_bar + pi0 + mu-

dkp_mup_nusep 

mu+ => nu_mu_bar + nu_e + e+

dkp_mum_nusep 

mu- => nu_mu + nu_e_bar + e-

dk_pip_numumup 

pi+ => nu_mu + mu+

dk_pim_numubmum 

pi- => nu_mu_bar + mu-

dk_mum_capture 

mu- + (A,Z) => nu (different kinematics from dkp_mum_nusep)

dkp_maximum 

one-beyond end for iterating

dkp_other 

flag for unusual cases

dkp_unknown 
dkp_k0l_nuepimep 

k0long => nu_e + pi- + e+

dkp_k0l_nuebpipem 

k0long => nu_e_bar + p+ + e-

dkp_k0l_numupimmup 

k0long => nu_mu + pi- + mu+

dkp_k0l_numubpipmum 

k0long => nu_mu_bar + pi+ + mu-

dkp_kp_numumup 

k+ => nu_mu + mu+

dkp_kp_nuepi0ep 

k+ => nu_e + pi0 + e+

dkp_kp_numupi0mup 

k+ => nu_mu + pi0 + mu+

dkp_kp_numubmum 

k- => nu_mu_bar + mu-

dkp_kp_nuebpi0em 

k- => nu_e_bar + pi0 + e-

dkp_kp_numubpi0mum 

k- => nu_mu_bar + pi0 + mu-

dkp_mup_nusep 

mu+ => nu_mu_bar + nu_e + e+

dkp_mum_nusep 

mu- => nu_mu + nu_e_bar + e-

dk_pip_numumup 

pi+ => nu_mu + mu+

dk_pim_numubmum 

pi- => nu_mu_bar + mu-

dk_mum_capture 

mu- + (A,Z) => nu (different kinematics from dkp_mum_nusep)

dkp_maximum 

one-beyond end for iterating

dkp_other 

flag for unusual cases

Definition at line 69 of file dk2nu.h.

69  {
70  dkp_unknown = 0,
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;
k0long => nu_mu_bar + pi+ + mu-
Definition: dk2nu.h:74
mu- => nu_mu + nu_e_bar + e-
Definition: dk2nu.h:82
k0long => nu_mu + pi- + mu+
Definition: dk2nu.h:73
k+ => nu_e + pi0 + e+
Definition: dk2nu.h:76
enum bsim::dkproc dkproc_t
one-beyond end for iterating
Definition: dk2nu.h:87
k+ => nu_mu + pi0 + mu+
Definition: dk2nu.h:77
k- => nu_e_bar + pi0 + e-
Definition: dk2nu.h:79
k0long => nu_e + pi- + e+
Definition: dk2nu.h:71
k- => nu_mu_bar + mu-
Definition: dk2nu.h:78
mu+ => nu_mu_bar + nu_e + e+
Definition: dk2nu.h:81
k+ => nu_mu + mu+
Definition: dk2nu.h:75
flag for unusual cases
Definition: dk2nu.h:88
k- => nu_mu_bar + pi0 + mu-
Definition: dk2nu.h:80
pi- => nu_mu_bar + mu-
Definition: dk2nu.h:84
k0long => nu_e_bar + p+ + e-
Definition: dk2nu.h:72
pi+ => nu_mu + mu+
Definition: dk2nu.h:83

Enumeration of decay processes, stored in "ndecay" stored as integer; these are for reference DK2NU: should there be an associated AsString() method that returns a text (optionally formatted for latex?)?

Enumerator
dkp_unknown 
dkp_k0l_nuepimep 

k0long => nu_e + pi- + e+

dkp_k0l_nuebpipem 

k0long => nu_e_bar + p+ + e-

dkp_k0l_numupimmup 

k0long => nu_mu + pi- + mu+

dkp_k0l_numubpipmum 

k0long => nu_mu_bar + pi+ + mu-

dkp_kp_numumup 

k+ => nu_mu + mu+

dkp_kp_nuepi0ep 

k+ => nu_e + pi0 + e+

dkp_kp_numupi0mup 

k+ => nu_mu + pi0 + mu+

dkp_kp_numubmum 

k- => nu_mu_bar + mu-

dkp_kp_nuebpi0em 

k- => nu_e_bar + pi0 + e-

dkp_kp_numubpi0mum 

k- => nu_mu_bar + pi0 + mu-

dkp_mup_nusep 

mu+ => nu_mu_bar + nu_e + e+

dkp_mum_nusep 

mu- => nu_mu + nu_e_bar + e-

dk_pip_numumup 

pi+ => nu_mu + mu+

dk_pim_numubmum 

pi- => nu_mu_bar + mu-

dk_mum_capture 

mu- + (A,Z) => nu (different kinematics from dkp_mum_nusep)

dkp_maximum 

one-beyond end for iterating

dkp_other 

flag for unusual cases

dkp_unknown 
dkp_k0l_nuepimep 

k0long => nu_e + pi- + e+

dkp_k0l_nuebpipem 

k0long => nu_e_bar + p+ + e-

dkp_k0l_numupimmup 

k0long => nu_mu + pi- + mu+

dkp_k0l_numubpipmum 

k0long => nu_mu_bar + pi+ + mu-

dkp_kp_numumup 

k+ => nu_mu + mu+

dkp_kp_nuepi0ep 

k+ => nu_e + pi0 + e+

dkp_kp_numupi0mup 

k+ => nu_mu + pi0 + mu+

dkp_kp_numubmum 

k- => nu_mu_bar + mu-

dkp_kp_nuebpi0em 

k- => nu_e_bar + pi0 + e-

dkp_kp_numubpi0mum 

k- => nu_mu_bar + pi0 + mu-

dkp_mup_nusep 

mu+ => nu_mu_bar + nu_e + e+

dkp_mum_nusep 

mu- => nu_mu + nu_e_bar + e-

dk_pip_numumup 

pi+ => nu_mu + mu+

dk_pim_numubmum 

pi- => nu_mu_bar + mu-

dk_mum_capture 

mu- + (A,Z) => nu (different kinematics from dkp_mum_nusep)

dkp_maximum 

one-beyond end for iterating

dkp_other 

flag for unusual cases

Definition at line 69 of file dk2nu.h.

69  {
70  dkp_unknown = 0,
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;
k0long => nu_mu_bar + pi+ + mu-
Definition: dk2nu.h:74
mu- => nu_mu + nu_e_bar + e-
Definition: dk2nu.h:82
k0long => nu_mu + pi- + mu+
Definition: dk2nu.h:73
k+ => nu_e + pi0 + e+
Definition: dk2nu.h:76
enum bsim::dkproc dkproc_t
one-beyond end for iterating
Definition: dk2nu.h:87
k+ => nu_mu + pi0 + mu+
Definition: dk2nu.h:77
k- => nu_e_bar + pi0 + e-
Definition: dk2nu.h:79
k0long => nu_e + pi- + e+
Definition: dk2nu.h:71
k- => nu_mu_bar + mu-
Definition: dk2nu.h:78
mu+ => nu_mu_bar + nu_e + e+
Definition: dk2nu.h:81
k+ => nu_mu + mu+
Definition: dk2nu.h:75
flag for unusual cases
Definition: dk2nu.h:88
k- => nu_mu_bar + pi0 + mu-
Definition: dk2nu.h:80
pi- => nu_mu_bar + mu-
Definition: dk2nu.h:84
k0long => nu_e_bar + p+ + e-
Definition: dk2nu.h:72
pi+ => nu_mu + mu+
Definition: dk2nu.h:83

All the data members are public as these classes are used as generalized structs. As they will be branches of a TTree no specialized naming indicators signifying that they are member data of a class will be used, nor will any fancy capitalization schemes.

All classes must implement a clear() method that resets their values to an identifiably invalid state or clears any vectors. Additionally classes should provide a AsString() method for formatting themselves

for use output.

Specialized enumerations Proposed flag bits:

Enumerator
kFlgOverflow 
kMaskReserved 
kMaskUser 
kFlgOverflow 
kMaskReserved 
kMaskUser 

Definition at line 57 of file dk2nu.h.

57  {
58  kFlgOverflow = 0x00000001,
59  kMaskReserved = 0x0000FFFF,
60  kMaskUser = 0xFFFF0000
61  } flgbitval_t;
enum bsim::flgbitval flgbitval_t

All the data members are public as these classes are used as generalized structs. As they will be branches of a TTree no specialized naming indicators signifying that they are member data of a class will be used, nor will any fancy capitalization schemes.

All classes must implement a clear() method that resets their values to an identifiably invalid state or clears any vectors. Additionally classes should provide a AsString() method for formatting themselves

for use output.

Specialized enumerations Proposed flag bits:

Enumerator
kFlgOverflow 
kMaskReserved 
kMaskUser 
kFlgOverflow 
kMaskReserved 
kMaskUser 

Definition at line 57 of file dk2nu.h.

57  {
58  kFlgOverflow = 0x00000001,
59  kMaskReserved = 0x0000FFFF,
60  kMaskUser = 0xFFFF0000
61  } flgbitval_t;
enum bsim::flgbitval flgbitval_t

Function Documentation

int bsim::calcEnuWgt ( const bsim::Decay decay,
const TVector3 &  xyz,
double &  enu,
double &  wgt_xy 
)

workhorse routine

Definition at line 48 of file calcLocationWeights.cxx.

References beta, om::cerr, dkp_mum_nusep, dkp_mup_nusep, allTimeWatchdog::endl, bsim::Decay::mupare, bsim::Decay::muparpx, bsim::Decay::muparpy, bsim::Decay::muparpz, bsim::Decay::ndecay, bsim::Decay::necm, bsim::Decay::norig, bsim::Decay::ntype, partial, bsim::Decay::pdpx, bsim::Decay::pdpy, bsim::Decay::pdpz, bsim::Decay::ppdxdz, bsim::Decay::ppdydz, bsim::Decay::ppenergy, bsim::Decay::pppz, bsim::Decay::ptype, Munits::rad, ana::Sqrt(), bsim::Decay::vx, bsim::Decay::vy, and bsim::Decay::vz.

Referenced by calcEnuWgt(), calcLocationWeights(), genie::flux::GDk2NuFlux::GenerateNext_weighted(), and genie::flux::GDk2NuFluxXMLHelper::~GDk2NuFluxXMLHelper().

50 {
51  // Neutrino Energy and Weight at arbitrary point
52  // Based on:
53  // NuMI-NOTE-BEAM-0109 (MINOS DocDB # 109)
54  // Title: Neutrino Beam Simulation using PAW with Weighted Monte Carlos
55  // Author: Rick Milburn
56  // Date: 1995-10-01
57 
58  // History:
59  // jzh 3/21/96 grab R.H.Milburn's weighing routine
60  // jzh 5/ 9/96 substantially modify the weighting function use dot product
61  // instead of rotation vecs to get theta get all info except
62  // det from ADAMO banks neutrino parent is in Particle.inc
63  // Add weighting factor for polarized muon decay
64  // jzh 4/17/97 convert more code to double precision because of problems
65  // with Enu>30 GeV
66  // rwh 10/ 9/08 transliterate function from f77 to C++
67 
68  // Original function description:
69  // Real function for use with PAW Ntuple To transform from destination
70  // detector geometry to the unit sphere moving with decaying hadron with
71  // velocity v, BETA=v/c, etc.. For (pseudo)scalar hadrons the decays will
72  // be isotropic in this sphere so the fractional area (out of 4-pi) is the
73  // fraction of decays that hit the target. For a given target point and
74  // area, and given x-y components of decay transverse location and slope,
75  // and given decay distance from target ans given decay GAMMA and
76  // rest-frame neutrino energy, the lab energy at the target and the
77  // fractional solid angle in the rest-frame are determined.
78  // For muon decays, correction for non-isotropic nature of decay is done.
79 
80  // Arguments:
81  // dk2nu :: contains current decay information
82  // xyz :: 3-vector of position to evaluate
83  // in *beam* frame coordinates (cm units)
84  // enu :: resulting energy
85  // wgt_xy :: resulting weight
86  // Return:
87  // (int) :: error code
88  // Assumptions:
89  // Energies given in GeV
90  // Particle codes have been translated from GEANT into PDG codes
91 
92  // for now ... these masses _should_ come from TDatabasePDG
93  // but use these hard-coded values to "exactly" reproduce old code
94  //
95  // old mass values are v01_07_* and before
96  // new mass values (v01_08_* and after) are Geant4 v4.10.3 values
97  //
98 #ifdef HISTORIC_MASS
99  const double kPIMASS = 0.13957;
100  const double kKMASS = 0.49368;
101  const double kK0MASS = 0.49767;
102  const double kMUMASS = 0.105658389;
103  const double kOMEGAMASS = 1.67245;
104 #else
105  const double kPIMASS = 0.1395701; // 0.13957;
106  const double kKMASS = 0.493677; // 0.49368;
107  const double kK0MASS = 0.497614; // 0.49767;
108  const double kMUMASS = 0.1056583715; // 0.105658389;
109  const double kOMEGAMASS = 1.67245; // 1.67245;
110 #endif
111  // from CLHEP/Units/PhysicalConstants.h
112  // used by Geant as CLHEP::neutron_mass_c2
113  const double kNEUTRONMASS = 0.93956536;
114 
115  const int kpdg_nue = 12; // extended Geant 53
116  const int kpdg_nuebar = -12; // extended Geant 52
117  const int kpdg_numu = 14; // extended Geant 56
118  const int kpdg_numubar = -14; // extended Geant 55
119 
120  const int kpdg_muplus = -13; // Geant 5
121  const int kpdg_muminus = 13; // Geant 6
122  const int kpdg_pionplus = 211; // Geant 8
123  const int kpdg_pionminus = -211; // Geant 9
124  const int kpdg_k0long = 130; // Geant 10 ( K0=311, K0S=310 )
125  const int kpdg_k0short = 310; // Geant 16
126  const int kpdg_k0mix = 311;
127  const int kpdg_kaonplus = 321; // Geant 11
128  const int kpdg_kaonminus = -321; // Geant 12
129  const int kpdg_omegaminus = 3334; // Geant 24
130  const int kpdg_omegaplus = -3334; // Geant 32
131  const int kpdg_neutron = 2112;
132  const int kpdg_antineutron = -2112;
133 
134  const double kRDET = 100.0; // set to flux per 100 cm radius
135 
136  double xpos = xyz.X();
137  double ypos = xyz.Y();
138  double zpos = xyz.Z();
139 
140  enu = 0.0; // don't know what the final value is
141  wgt_xy = 0.0; // but set these in case we return early due to error
142 
143 
144  // in principle we should get these from the particle DB
145  // but for consistency testing use the hardcoded values
146  double parent_mass = kPIMASS;
147  switch ( decay.ptype ) {
148  case kpdg_pionplus:
149  case kpdg_pionminus:
150  parent_mass = kPIMASS;
151  break;
152  case kpdg_kaonplus:
153  case kpdg_kaonminus:
154  parent_mass = kKMASS;
155  break;
156  case kpdg_k0long:
157  case kpdg_k0short:
158  case kpdg_k0mix:
159  parent_mass = kK0MASS;
160  break;
161  case kpdg_muplus:
162  case kpdg_muminus:
163  parent_mass = kMUMASS;
164  break;
165  case kpdg_omegaminus:
166  case kpdg_omegaplus:
167  parent_mass = kOMEGAMASS;
168  break;
169  case kpdg_neutron:
170  case kpdg_antineutron:
171  parent_mass = kNEUTRONMASS;
172  break;
173  default:
174  std::cerr << "bsim::calcEnuWgt unknown particle type " << decay.ptype
175  << std::endl << std::flush;
176  enu = 0.0;
177  wgt_xy = 0.0;
178  return 1;
179  }
180 
181  double parentp2 = ( decay.pdpx*decay.pdpx +
182  decay.pdpy*decay.pdpy +
183  decay.pdpz*decay.pdpz );
184  double parent_energy = TMath::Sqrt( parentp2 +
185  parent_mass*parent_mass);
186  double parentp = TMath::Sqrt( parentp2 );
187 
188  double gamma = parent_energy / parent_mass;
189  double gamma_sqr = gamma * gamma;
190  double beta_mag = TMath::Sqrt( ( gamma_sqr - 1.0 )/gamma_sqr );
191 
192  // Get the neutrino energy in the parent decay CM
193  double enuzr = decay.necm;
194  // Get angle from parent line of flight to chosen point in beam frame
195  double rad = TMath::Sqrt( (xpos-decay.vx)*(xpos-decay.vx) +
196  (ypos-decay.vy)*(ypos-decay.vy) +
197  (zpos-decay.vz)*(zpos-decay.vz) );
198 
199  double emrat = 1.0;
200  double costh_pardet = -999., theta_pardet = -999.;
201 
202  // boost correction, but only if parent hasn't stopped
203  if ( parentp > 0. ) {
204  costh_pardet = ( decay.pdpx*(xpos-decay.vx) +
205  decay.pdpy*(ypos-decay.vy) +
206  decay.pdpz*(zpos-decay.vz) )
207  / ( parentp * rad);
208  if ( costh_pardet > 1.0 ) costh_pardet = 1.0;
209  if ( costh_pardet < -1.0 ) costh_pardet = -1.0;
210  theta_pardet = TMath::ACos(costh_pardet);
211 
212  // Weighted neutrino energy in beam, approx, good for small theta
213  emrat = 1.0 / ( gamma * ( 1.0 - beta_mag * costh_pardet ));
214  }
215 
216  enu = emrat * enuzr; // the energy ... normally
217 
218  // Get solid angle/4pi for detector element
219  // small angle approximation, fixed by Alex Radovic
220  //SAA// double sangdet = ( kRDET*kRDET /
221  //SAA// ( (zpos-decay.vz)*(zpos-decay.vz) ) ) / 4.0;
222  double sanddetcomp = TMath::Sqrt( ( (xpos-decay.vx)*(xpos-decay.vx) ) +
223  ( (ypos-decay.vy)*(ypos-decay.vy) ) +
224  ( (zpos-decay.vz)*(zpos-decay.vz) ) );
225  double sangdet = (1.0-TMath::Cos(TMath::ATan( kRDET / sanddetcomp )))/2.0;
226 
227  // Weight for solid angle and lorentz boost
228  wgt_xy = sangdet * ( emrat * emrat ); // ! the weight ... normally
229 
230  // Done for all except polarized muon decay
231  // in which case need to modify weight
232  // (must be done in double precision)
233  // BUT do this only for case of muon decay, not muon capture
234  // until beamline simulation code gets updated these generally show up as
235  // decay.ndecay == 0, but certainly not dkp_mup_nusep or dkp_mum_nusep
236  // so was:
237  // if ( decay.ptype == kpdg_muplus ||
238  // decay.ptype == kpdg_muminus ) {
239  // now:
240  if ( decay.ndecay == bsim::dkp_mup_nusep ||
241  decay.ndecay == bsim::dkp_mum_nusep ) {
242  double beta[3], p_dcm_nu[4], p_nu[3], p_pcm_mp[3], partial;
243 
244  // Boost neu neutrino to mu decay CM
245  beta[0] = decay.pdpx / parent_energy;
246  beta[1] = decay.pdpy / parent_energy;
247  beta[2] = decay.pdpz / parent_energy;
248  p_nu[0] = (xpos-decay.vx)*enu/rad;
249  p_nu[1] = (ypos-decay.vy)*enu/rad;
250  p_nu[2] = (zpos-decay.vz)*enu/rad;
251  partial = gamma *
252  (beta[0]*p_nu[0] + beta[1]*p_nu[1] + beta[2]*p_nu[2] );
253  partial = enu - partial/(gamma+1.0);
254  // the following calculation is numerically imprecise
255  // especially p_dcm_nu[2] leads to taking the difference of numbers
256  // of order ~10's and getting results of order ~0.02's
257  // for g3numi we're starting with floats (ie. good to ~1 part in 10^7)
258  p_dcm_nu[0] = p_nu[0] - beta[0]*gamma*partial;
259  p_dcm_nu[1] = p_nu[1] - beta[1]*gamma*partial;
260  p_dcm_nu[2] = p_nu[2] - beta[2]*gamma*partial;
261  p_dcm_nu[3] = TMath::Sqrt( p_dcm_nu[0]*p_dcm_nu[0] +
262  p_dcm_nu[1]*p_dcm_nu[1] +
263  p_dcm_nu[2]*p_dcm_nu[2] );
264 
265  // Boost parent of mu to mu production CM
266  double particle_energy = decay.ppenergy;
267  gamma = particle_energy/parent_mass;
268  beta[0] = decay.ppdxdz * decay.pppz / particle_energy;
269  beta[1] = decay.ppdydz * decay.pppz / particle_energy;
270  beta[2] = decay.pppz / particle_energy;
271  partial = gamma * ( beta[0]*decay.muparpx +
272  beta[1]*decay.muparpy +
273  beta[2]*decay.muparpz );
274  partial = decay.mupare - partial/(gamma+1.0);
275  p_pcm_mp[0] = decay.muparpx - beta[0]*gamma*partial;
276  p_pcm_mp[1] = decay.muparpy - beta[1]*gamma*partial;
277  p_pcm_mp[2] = decay.muparpz - beta[2]*gamma*partial;
278  double p_pcm = TMath::Sqrt ( p_pcm_mp[0]*p_pcm_mp[0] +
279  p_pcm_mp[1]*p_pcm_mp[1] +
280  p_pcm_mp[2]*p_pcm_mp[2] );
281 
282  const double eps = 1.0e-30; // ? what value to use
283  if ( p_pcm < eps || p_dcm_nu[3] < eps ) {
284  return 3; // mu missing parent info?
285  }
286  // Calc new decay angle w.r.t. (anti)spin direction
287  double costh = ( p_dcm_nu[0]*p_pcm_mp[0] +
288  p_dcm_nu[1]*p_pcm_mp[1] +
289  p_dcm_nu[2]*p_pcm_mp[2] ) /
290  ( p_dcm_nu[3]*p_pcm );
291  // protect against small excursions
292  if ( costh > 1.0 ) costh = 1.0;
293  if ( costh < -1.0 ) costh = -1.0;
294  // Calc relative weight due to angle difference
295  double wgt_ratio = 0.0;
296  switch ( decay.ntype ) {
297  case kpdg_nue:
298  case kpdg_nuebar:
299  wgt_ratio = 1.0 - costh;
300  break;
301  case kpdg_numu:
302  case kpdg_numubar:
303  {
304  double xnu = 2.0 * enuzr / kMUMASS;
305  wgt_ratio = ( (3.0-2.0*xnu ) - (1.0-2.0*xnu)*costh ) / (3.0-2.0*xnu);
306 
307  if ( wgt_ratio < 0.0 ) {
308  std::cerr << "bsim::calcEnuWgt encountered serious problem: "
309  << " wgt_ratio " << wgt_ratio
310  << " enu " << enu << " costh " << costh << " xnu " << xnu
311  << " enuzr=decay.necm " << enuzr << " kMUMASS " << kMUMASS
312  << " norig " << decay.norig
313  << " ndecay " << decay.ndecay
314  << " ntype " << decay.ntype
315  << " ptype " << decay.ptype
316  << std::endl;
317  enu = 0;
318  wgt_xy = 0;
319  return 4; // bad, bad, bad calculation
320  }
321  break;
322  }
323  default:
324  enu = 0.0;
325  wgt_xy = 0.0;
326  return 2; // bad neutrino type for muon decay
327  }
328 
329  wgt_xy = wgt_xy * wgt_ratio;
330 
331  } // ptype is muon
332 
333  return 0;
334 }
Double_t ppdydz
% direction of nu parent at its production point
Definition: dk2nu.h:139
Double_t muparpx
%
Definition: dk2nu.h:160
Double_t pppz
% z momentum of nu parent at its production point
Definition: dk2nu.h:140
Double_t pdpx
% px momentum of nu parent at (vx,vy,vz)
Definition: dk2nu.h:133
Double_t necm
% nu energy in center-of-mass frame
Definition: dk2nu.h:165
OStream cerr
Definition: OStream.cxx:7
Double_t beta
Double_t ppdxdz
% direction of nu parent at its production point
Definition: dk2nu.h:138
mu- => nu_mu + nu_e_bar + e-
Definition: dk2nu.h:82
Definition: lz4.cxx:387
Double_t vz
% neutrino production vertex z
Definition: dk2nu.h:132
double enu
Definition: runWimpSim.h:113
Double_t pdpz
% pz momentum of nu parent at (vx,vy,vz)
Definition: dk2nu.h:135
Double_t muparpy
%
Definition: dk2nu.h:161
Double_t muparpz
%
Definition: dk2nu.h:162
Double_t ppenergy
% energy of nu parent at its production point
Definition: dk2nu.h:141
Int_t ntype
% neutrino flavor (PDG? code)
Definition: dk2nu.h:128
static constexpr Double_t rad
Definition: Munits.h:162
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
mu+ => nu_mu_bar + nu_e + e+
Definition: dk2nu.h:81
Double_t vx
% neutrino production vertex x
Definition: dk2nu.h:130
Int_t ndecay
decay process (see dkproc_t)
Definition: dk2nu.h:127
Double_t mupare
% energy of nu grandparent
Definition: dk2nu.h:163
Double_t vy
% neutrino production vertex y
Definition: dk2nu.h:131
Int_t norig
not used?
Definition: dk2nu.h:126
Double_t pdpy
% py momentum of nu parent at (vx,vy,vz)
Definition: dk2nu.h:134
Int_t ptype
% nu parent species (PDG? code)
Definition: dk2nu.h:144
int bsim::calcEnuWgt ( const bsim::Dk2Nu dk2nu,
const TVector3 &  xyz,
double &  enu,
double &  wgt_xy 
)

alternate interface

Definition at line 337 of file calcLocationWeights.cxx.

References calcEnuWgt(), and bsim::Dk2Nu::decay.

339 {
340  return bsim::calcEnuWgt(dk2nu->decay,xyz,enu,wgt_xy);
341 }
bsim::Decay decay
basic decay information
Definition: dk2nu.h:327
double enu
Definition: runWimpSim.h:113
int calcEnuWgt(const bsim::Decay &decay, const TVector3 &xyz, double &enu, double &wgt_xy)
workhorse routine
void bsim::calcLocationWeights ( const bsim::DkMeta dkmeta,
bsim::Dk2Nu dk2nu 
)

user interface

Definition at line 10 of file calcLocationWeights.cxx.

References calcEnuWgt(), om::cerr, bsim::Dk2Nu::decay, allTimeWatchdog::endl, bsim::DkMeta::location, bsim::Dk2Nu::nuray, make_associated_cosmic_defs::p3, fabricate::status, string, Unit(), bsim::Decay::vx, bsim::Decay::vy, and bsim::Decay::vz.

Referenced by convert_flugg(), convert_g4lbne(), convert_g4minerva(), and test_fill_dk2nu().

11 {
12  size_t nloc = dkmeta->location.size();
13  for (size_t iloc = 0; iloc < nloc; ++iloc ) {
14  // skip calculation for random location ... should already be filled
15  const std::string rkey = "random decay";
16  if ( dkmeta->location[iloc].name == rkey ) {
17  if ( iloc != 0 ) {
18  std::cerr << "calcLocationWeights \"" << rkey << "\""
19  << " isn't the 0-th entry" << std::endl;
20  return; // assert(0);
21  }
22  if ( dk2nu->nuray.size() != 1 ) {
23  std::cerr << "calcLocationWeights \"" << rkey << "\""
24  << " nuenergy[" << iloc << "] not filled" << std::endl;
25  return; //assert(0);
26  }
27  continue;
28  }
29  TVector3 xyzDet(dkmeta->location[iloc].x,
30  dkmeta->location[iloc].y,
31  dkmeta->location[iloc].z); // position to evaluate
32  double enu_xy = 0; // give a default value
33  double wgt_xy = 0; // give a default value
34  int status = bsim::calcEnuWgt(dk2nu,xyzDet,enu_xy,wgt_xy);
35  if ( status != 0 ) {
36  std::cerr << "bsim::calcEnuWgt returned " << status << " for "
37  << dkmeta->location[iloc].name << std::endl;
38  }
39  // with the recalculated energy compute the momentum components
40  TVector3 xyzDk(dk2nu->decay.vx,dk2nu->decay.vy,dk2nu->decay.vz); // origin of decay
41  TVector3 p3 = enu_xy * (xyzDet - xyzDk).Unit();
42  bsim::NuRay anuray(p3.x(), p3.y(), p3.z(), enu_xy, wgt_xy);
43  dk2nu->nuray.push_back(anuray);
44  }
45 }
std::vector< bsim::Location > location
locations
Definition: dkmeta.h:119
int status
Definition: fabricate.py:1613
std::vector< bsim::NuRay > nuray
rays through detector fixed points
Definition: dk2nu.h:328
OStream cerr
Definition: OStream.cxx:7
bsim::Decay decay
basic decay information
Definition: dk2nu.h:327
Double_t vz
% neutrino production vertex z
Definition: dk2nu.h:132
TVector3 Unit() const
int calcEnuWgt(const bsim::Decay &decay, const TVector3 &xyz, double &enu, double &wgt_xy)
workhorse routine
Double_t vx
% neutrino production vertex x
Definition: dk2nu.h:130
Double_t vy
% neutrino production vertex y
Definition: dk2nu.h:131
enum BeamMode string
Bool_t bsim::IsDefault ( Float_t  intval)

Definition at line 4 of file dflt.cxx.

References e, and kDfltFloat.

5 {
6  return (TMath::Abs(fval - kDfltFloat) < 1e-4);
7 }
static const Float_t kDfltFloat
Definition: dflt.h:12
Float_t e
Definition: plot.C:35
Bool_t bsim::IsDefault ( Int_t  intval)

Definition at line 9 of file dflt.cxx.

References kDfltInt.

10 {
11  return (ival == kDfltInt);
12 }
static const Int_t kDfltInt
Definition: dflt.h:13
ival
Definition: test.py:10
Bool_t bsim::IsDefault ( UInt_t  intval)

Definition at line 14 of file dflt.cxx.

References kDfltUInt.

15 {
16  return (uival == kDfltUInt);
17 }
static const UInt_t kDfltUInt
Definition: dflt.h:14
Bool_t bsim::IsDefault ( Double_t  intval)

Definition at line 19 of file dflt.cxx.

References e, and kDfltDouble.

20 {
21  return (TMath::Abs(dval - kDfltDouble) < 1e-4);
22 }
static const Double_t kDfltDouble
Definition: dflt.h:15
Float_t e
Definition: plot.C:35
Bool_t bsim::IsDefault ( Bool_t  intval)

Definition at line 24 of file dflt.cxx.

References kDfltBool.

25 {
26  return (bval == kDfltBool);
27 }
static const Bool_t kDfltBool
Definition: dflt.h:16
void bsim::printWeightLocations ( std::vector< bsim::Location > &  locations)

print the locations

Definition at line 83 of file readWeightLocations.cxx.

References om::cout, submit_hadd::l, and fhicl::detail::nl().

Referenced by ConvertInit(), and printWeightLocations().

84 {
85  size_t nl = locations.size();
86  std::cout << nl << " locations:\n";
87  for ( size_t l = 0; l < nl; ++l ) {
88  std::cout << " [" << std::setw(2) << l << "] " << locations[l] << "\n";
89  }
90 }
OStream cout
Definition: OStream.cxx:6
std::string nl(std::size_t i=1)
void bsim::printWeightLocations ( bsim::DkMeta dkmeta)

a variant that prints the locations in the dkmeta object

Definition at line 92 of file readWeightLocations.cxx.

References bsim::DkMeta::location, and printWeightLocations().

93 {
95 }
std::vector< bsim::Location > location
locations
Definition: dkmeta.h:119
void printWeightLocations(std::vector< bsim::Location > &locations)
print the locations
void bsim::readWeightLocations ( std::string  locfilename,
std::vector< bsim::Location > &  locations 
)

Read a text file that contains a header line followed by quartets of "<xpos> <ypos> <zpos> <text string>" on separate lines. Fill the supplied vectors. Trim off leading/trailing blanks and quotes (single/double) from the string. Convention has it that positions are given in (cm).

Definition at line 17 of file readWeightLocations.cxx.

References plot_validation_datamc::c, MECModelEnuComparisons::i, path, string, tmp, submit_syst::x, submit_syst::y, and test::z.

Referenced by ConvertInit(), readWeightLocations(), test_fill_dk2nu(), and test_read_locations().

19 {
20 
21  const char* path = gSystem->ExpandPathName(locfilename.c_str());
22  std::ifstream locfile(path);
23 
24  int iline=0;
25 
26  char comment_buffer[1000];
27 
28  // read lines
29  char tmp[1001];
30  size_t tmplen = sizeof(tmp);
31  while ( ! locfile.eof() ) {
32  char c;
33  while ( ( c = locfile.get() ) == ' ' ) {}; // eat leading spaces
34  if ( c == '#' ) {
35  // comment_buffer line
36  locfile.getline(comment_buffer,sizeof(comment_buffer));
37  continue;
38  }
39  locfile.putback(c);
40  double x, y, z;
41  locfile >> x >> y >> z;
42  locfile.getline(tmp,tmplen-1);
43  size_t i = locfile.gcount();
44  // make sure the c-string is null terminated
45  size_t inull = i;
46  //if ( inull < 0 ) inull = 0;
47  if ( inull > tmplen-1 ) inull = tmplen-1;
48  tmp[inull] = '\0';
49  std::string name(tmp);
50  // ignore leading & trailing blanks (and any single/double quotes)
51  size_t ilast = name.find_last_not_of(" \t\n'\"");
52  name.erase(ilast+1,std::string::npos); // trim tail
53  size_t ifirst = name.find_first_not_of(" \t\n'\"");
54  name.erase(0,ifirst); // trim head
55 
56  ++iline;
57  if ( ! locfile.good() ) {
58  //if ( verbose)
59  // std::cout << "stopped reading on line " << iline << std::endl;
60  break;
61  }
62  bsim::Location alocation(x,y,z,name);
63  locations.push_back(alocation);
64  }
65 
66 }
const XML_Char * name
Definition: expat.h:151
Float_t tmp
Definition: plot.C:36
z
Definition: test.py:28
const std::string path
Definition: plot_BEN.C:43
enum BeamMode string
void bsim::readWeightLocations ( std::string  locfilename,
bsim::DkMeta dkmeta 
)

a variant that will fill the dkmeta object

read & print the locations where weights are to be calculated

make an entry for the random decay

read and parse the text file for additional positions use the vector version

Definition at line 69 of file readWeightLocations.cxx.

References bsim::DkMeta::location, and readWeightLocations().

70 {
71  /// read & print the locations where weights are to be calculated
72 
73  std::vector<bsim::Location>& locations = dkmeta->location;
74  /// make an entry for the random decay
75  bsim::Location rndmloc(0,0,0,"random decay");
76  locations.push_back(rndmloc);
77 
78  /// read and parse the text file for additional positions
79  /// use the vector version
80  readWeightLocations(locfilename, locations);
81 }
std::vector< bsim::Location > location
locations
Definition: dkmeta.h:119
void readWeightLocations(std::string locfilename, std::vector< bsim::Location > &locations)

Variable Documentation

const Bool_t bsim::kDfltBool = false
static

Definition at line 16 of file dflt.h.

Referenced by IsDefault().

const Bool_t bsim::kDfltBool = false
static

Definition at line 16 of file dflt.h.

const Double_t bsim::kDfltDouble = -9999.99
static

Definition at line 15 of file dflt.h.

const Double_t bsim::kDfltDouble = -9999.99
static

Definition at line 15 of file dflt.h.

Referenced by IsDefault().

const Float_t bsim::kDfltFloat = -9999.99
static

Definition at line 12 of file dflt.h.

Referenced by IsDefault().

const Float_t bsim::kDfltFloat = -9999.99
static

Definition at line 12 of file dflt.h.

const Int_t bsim::kDfltInt = -9999
static

Definition at line 13 of file dflt.h.

Referenced by IsDefault().

const Int_t bsim::kDfltInt = -9999
static

Definition at line 13 of file dflt.h.

const UInt_t bsim::kDfltUInt = 9999
static

Definition at line 14 of file dflt.h.

Referenced by IsDefault().

const UInt_t bsim::kDfltUInt = 9999
static

Definition at line 14 of file dflt.h.