Jpp  19.0.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JVisibleEnergyToolkit.hh
Go to the documentation of this file.
1 #ifndef __JSIRENE__JVISIBLEENERGYTOOLKIT__
2 #define __JSIRENE__JVISIBLEENERGYTOOLKIT__
3 
8 
10 #include "JGeometry2D/JCircle2D.hh"
11 
14 
15 #include "JPhysics/JConstants.hh"
16 #include "JPhysics/JPDFToolkit.hh"
18 #include "JPhysics/JGeane.hh"
19 
20 #include "JAAnet/JAAnetToolkit.hh"
21 #include "JAAnet/JParticleTypes.hh"
22 #include "JAAnet/JPDB.hh"
23 
24 #include "JSirene/pythia.hh"
25 
26 
27 /**
28  * \file
29  *
30  * Auxiliary methods for evaluating visible energies.
31  * \author bjung
32  */
33 namespace JSIRENE {}
34 namespace JPP { using namespace JSIRENE; }
35 
36 namespace JSIRENE {
37 
39 
40 
41  /**
42  * Auxiliary function to retrieve the maximum cylindrical containment volume.
43  *
44  * \return maximum cylindrical containment volume
45  */
47  {
48  using namespace JPP;
49 
50  const double R_Earth = R_EARTH_KM * 1e3; // m
51 
52  const JVector2D center(0.0, 0.0);
53  const JCircle2D circle(center, R_Earth);
54 
55  return JCylinder3D(circle, -R_Earth, R_Earth);
56  }
57 
58 
59  /**
60  * Get the visible energy of a track.\n
61  * This method accounts for muon radiative energy losses.
62  *
63  * \param track track
64  * \param can detector can
65  * \return visible energy [GeV]
66  */
67  double getVisibleEnergy(const Trk& track,
69  {
70  using namespace std;
71  using namespace JPP;
72 
73  double Evis = 0.0;
74 
75  if (is_finalstate(track)) {
76 
77  if (is_muon(track)) {
78 
79  // Determine muon pathlength inside detector [m]
80 
81  const JCylinder3D::intersection_type& intersection = can.getIntersection(getAxis(track));
82 
83  const double Lmuon = gWater.getX(track.E, MASS_MUON / getSinThetaC());
84  const double Leff = (min(Lmuon, max(intersection.second, 0.0)) -
85  min(Lmuon, max(intersection.first, 0.0)));
86 
87 
88  // Determine visible energy deposition [GeV]
89 
90  const double Emidpoint = gWater.getE(track.E, Lmuon/2.0);
91 
92  const double dEb = gWater.getEb(track.E, Leff);
93  const double dEc = Leff / geanc();
94  const double dEd = Leff * getDeltaRaysFromMuon(Emidpoint);
95 
96  Evis = dEb + dEc + dEd;
97 
98  } else if (!is_neutrino(track) && JPDB::getInstance().hasPDG(track.type)) {
99 
100  Evis = pythia(track.type, getKineticEnergy(track));
101  }
102  }
103 
104  return Evis;
105  }
106 
107 
108  /**
109  * Get the visible energy vector of a track.\n
110  * This method accounts for muon radiative energy losses.
111  *
112  * \param track track
113  * \param can detector can
114  * \return visible energy vector [GeV]
115  */
116  inline Vec getVisibleEnergyVector(const Trk& track,
117  const JCylinder3D& can = getMaximumContainmentVolume()) {
118  return getVisibleEnergy(track, can) * track.dir;
119  }
120 
121 
122  /**
123  * Get the visible energy of a given range of tracks.\n
124  * This method accounts for muon radiative energy losses.
125  *
126  * \param __begin start of track data
127  * \param __end end of track data
128  * \param can detector can
129  * \return visible energy [GeV]
130  */
134  {
135  using namespace std;
136 
137  double Evis = 0.0;
138 
139  for (vector<Trk>::const_iterator track = __begin; track != __end; ++track) {
140  Evis += getVisibleEnergy(*track, can);
141  }
142 
143  return Evis;
144  }
145 
146 
147  /**
148  * Get the visible energy vector of a given range of tracks.\n
149  * This method accounts for muon radiative energy losses.
150  *
151  * \param __begin start of track data
152  * \param __end end of track data
153  * \param can detector can
154  * \return visible energy vector [GeV]
155  */
159  {
160  using namespace std;
161 
162  Vec Evis(0.0, 0.0, 0.0);
163 
164  for (vector<Trk>::const_iterator track = __begin; track != __end; ++track) {
165  Evis += getVisibleEnergyVector(*track, can);
166  }
167 
168  return Evis;
169  }
170 
171 
172  /**
173  * Get the visible energy vector of an event.\n
174  * This method accounts for muon radiative energy losses.
175  *
176  * \param evt event
177  * \param can detector can
178  * \return visible energy [GeV]
179  */
180  inline double getVisibleEnergy(const Evt& evt,
182  {
183  return getVisibleEnergy(evt.mc_trks.begin(), evt.mc_trks.end(), can);
184  }
185 
186 
187  /**
188  * Get the visible energy vector of an event.\n
189  * This method accounts for muon radiative energy losses.
190  *
191  * \param evt event
192  * \param can detector can
193  * \return visible energy vector [GeV]
194  */
195  inline Vec getVisibleEnergyVector(const Evt& evt,
197  {
198  return getVisibleEnergyVector(evt.mc_trks.begin(), evt.mc_trks.end(), can);
199  }
200 }
201 
202 #endif
Data structure for vector in two dimensions.
Definition: JVector2D.hh:32
Auxiliary methods for physics calculations.
static const double R_EARTH_KM
Geophysics constants.
Auxiliary methods for PDF calculations.
double geanc()
Equivalent muon track length per unit shower energy.
Definition: JGeane.hh:28
Energy loss of muon.
This file contains converted Fortran code from km3.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Vec dir
track direction
Definition: Trk.hh:18
static const double MASS_MUON
muon mass [GeV]
Vec getVisibleEnergyVector(const Trk &track, const JCylinder3D &can=getMaximumContainmentVolume())
Get the visible energy vector of a track.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
virtual double getX(const double E0, const double E1) const override
Get distance traveled by muon.
Definition: JGeane.hh:349
static const JPythia pythia
Function object for relative light yield as a function of GEANT particle code.
Definition: JPythia.hh:96
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
double getEb(const double E, const double dx) const
Get energy loss due to pair production and bremsstrahlung.
Definition: JGeane.hh:334
static const JPDB & getInstance()
Get particle data book.
Definition: JPDB.hh:131
double getDeltaRaysFromMuon(const double E, const JRange< double > T_GeV=JRange< double >(DELTARAY_TMIN, DELTARAY_TMAX))
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:260
bool is_neutrino(const Trk &track)
Test whether given track is a neutrino.
bool is_finalstate(const Trk &track)
Test whether given track corresponds to a final state particle.
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition: Vec.hh:12
Cylinder object.
Definition: JCylinder3D.hh:39
JAxis3D getAxis(const Trk &track)
Get axis.
const JCylinder3D getMaximumContainmentVolume()
Auxiliary function to retrieve the maximum cylindrical containment volume.
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
double getVisibleEnergy(const Trk &track, const JCylinder3D &can=getMaximumContainmentVolume())
Get the visible energy of a track.
Physics constants.
double getKineticEnergy(const Trk &trk)
Get track kinetic energy.
int type
MC: particle type in PDG encoding.
Definition: Trk.hh:24
virtual double getE(const double E, const double dx) const override
Get energy of muon after specified distance.
Definition: JGeane.hh:257
Definition of particle types.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
std::vector< Trk > mc_trks
MC: list of MC truth tracks.
Definition: Evt.hh:49
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20
double getSinThetaC()
Get average sine of Cherenkov angle of water corresponding to group velocity.