Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JSireneToolkit.hh
Go to the documentation of this file.
1 #ifndef __JSIRENE__JSIRENETOOLKIT__
2 #define __JSIRENE__JSIRENETOOLKIT__
3 
4 #include <vector>
5 #include <cmath>
6 
7 #include "TRandom3.h"
8 
9 #include "JTools/JConstants.hh"
10 #include "JPhysics/JPDFTypes.hh"
11 #include "JPhysics/JGeane.hh"
12 #include "JAAnet/JAAnetToolkit.hh"
13 #include "JLang/JComparator.hh"
14 
15 
16 /**
17  * \author mdejong
18  */
19 
20 
21 /**
22  * Detector simulations.
23  */
24 namespace JSIRENE {}
25 namespace JPP { using namespace JSIRENE; }
26 
27 namespace JSIRENE {
28 
30  using JTOOLS::MASS_MUON;
32  using JPHYSICS::JGeane;
33  using JPHYSICS::gWater;
35  using JAANET::JHitType_t;
36 
37 
38  /**
39  * Get number of photo-electrons of a hit given the expectation
40  * values of the number of photo-electrons on a module and PMT.
41  *
42  * The return value is evaluated by pick-and-drop statistics from
43  * the generated number of photo-electrons when the expectation
44  * value of the number of photo-electrons on a module deviates
45  * less than 5 sigmas from 0 (i.e.\ when it is less than 25).
46  * Otherwise, the return value is evaluated by Poisson statistics
47  * from the expectation value of the number of photo-electrons on PMT.
48  *
49  * \param NPE expectation value of npe on module
50  * \param N generated value of npe on module
51  * \param npe expectation value of npe on PMT
52  * \return number of photo-electrons on PMT
53  */
54  inline int getNumberOfPhotoElectrons(const double NPE,
55  const int N,
56  const double npe)
57  {
58  if (NPE < 25.0) {
59 
60  int n = 0;
61 
62  for (int i = N; i != 0; --i) {
63  if (gRandom->Rndm() * NPE < npe) {
64  ++n;
65  }
66  }
67 
68  return n;
69 
70  } else {
71 
72  return gRandom->Poisson(npe);
73  }
74  }
75 
76 
77  /**
78  * Get number of photo-electrons of a hit given number of photo-electrons on PMT.
79  *
80  * The number of photo-electrons of a hit may be larger than unity
81  * to limit the overall number of hits and consequently the memory coonsumption as well as
82  * the number of times the arrival time needs to be evaluated which is CPU intensive.
83  *
84  * \param npe number of photo-electrons on PMT
85  * \return number of photo-electrons of hit
86  */
87  inline int getNumberOfPhotoElectrons(const int npe)
88  {
89  static const int NPE = 20;
90 
91  if (npe > NPE) {
92 
93  const int n = (int) (NPE * log10((double) npe / (double) NPE));
94 
95  if (n == 0)
96  return 1;
97  else
98  return n;
99  }
100 
101  return 1;
102  }
103 
104 
105  /**
106  * Get hit type corresponding to given PDF type.
107  *
108  * \param pdf PDF type
109  * \param shower force origin from neutrino interaction shower
110  * \return hit type
111  */
112  inline JHitType_t getHitType(const JPDFType_t pdf, const bool shower = false)
113  {
114  using namespace JAANET;
115  using namespace JPHYSICS;
116 
117  switch (pdf) {
118 
120  return HIT_TYPE_MUON_DIRECT;
121 
124 
127 
130 
133 
136 
137  default:
138  return HIT_TYPE_UNKNOWN;
139  }
140  }
141 
142 
143  /**
144  * Vertex of energy loss of muon.
145  */
146  struct JVertex {
147  /**
148  * Default constructor.
149  */
151  z (0.0),
152  t (0.0),
153  E (0.0),
154  Es(0.0)
155  {}
156 
157 
158  /**
159  * Constructor.
160  *
161  * \param z position [m]
162  * \param t time [ns]
163  * \param E energy [GeV]
164  */
165  JVertex(const double z,
166  const double t,
167  const double E)
168  {
169  this->z = z;
170  this->t = t;
171  this->E = E;
172  this->Es = 0.0;
173  }
174 
175 
176  /**
177  * Get position.
178  *
179  * \return position [m]
180  */
181  double getZ() const
182  {
183  return z;
184  }
185 
186 
187  /**
188  * Get time.
189  *
190  * \return time [ns]
191  */
192  double getT() const
193  {
194  return t;
195  }
196 
197 
198  /**
199  * Get muon energy.
200  *
201  * \return energy [GeV]
202  */
203  double getE() const
204  {
205  return E;
206  }
207 
208 
209  /**
210  * Get shower energy.
211  *
212  * \return energy [GeV]
213  */
214  double getEs() const
215  {
216  return Es;
217  }
218 
219 
220  /**
221  * Get range of muon.
222  * This method applies only ionisation energy loss.
223  *
224  * \return range [m]
225  */
226  double getRange() const
227  {
228  if (E > MASS_MUON * getIndexOfRefraction())
229  return (E - MASS_MUON * getIndexOfRefraction()) / gWater.getA();
230  else
231  return 0.0;
232  }
233 
234 
235  /**
236  * Get range of muon.
237  *
238  * \param geane energy loss
239  * \return range [m]
240  */
241  double getRange(const JGeane& geane) const
242  {
243  return geane(E);
244  }
245 
246 
247  /**
248  * Apply energy loss energy.
249  *
250  * \param Es energy [GeV]
251  */
252  void applyEloss(const double Es)
253  {
254  this->E -= Es;
255  this->Es = Es;
256  }
257 
258 
259  /**
260  * Step.
261  * This method applies only ionisation energy loss.
262  *
263  * \param ds step [m]
264  * \return this vertex
265  */
266  JVertex& step(const double ds)
267  {
268  z += ds;
269  t += ds / getSpeedOfLight();
270  E -= ds * gWater.getA();
271 
272  return *this;
273  }
274 
275 
276  /**
277  * Step.
278  *
279  * \param geane energy loss
280  * \param ds step [m]
281  * \return this vertex
282  */
283  JVertex& step(const JGeane& geane, const double ds)
284  {
285  z += ds;
286  t += ds / getSpeedOfLight();
287  E = geane(E, ds);
288 
289  return *this;
290  }
291 
292  protected:
293  double z;
294  double t;
295  double E;
296  double Es;
297  };
298 
299 
300  /**
301  * Muon trajectory.
302  */
303  struct JTrack :
304  public std::vector<JVertex>
305  {
306  /**
307  * Constructor.
308  *
309  * \param vertex muon starting point
310  */
311  JTrack(const JVertex& vertex) :
312  std::vector<JVertex>(1, vertex)
313  {}
314 
315 
316  /**
317  * Get muon energy at given position.
318  *
319  * \param z position [m]
320  * \return energy [GeV]
321  */
322  double getE(const double z) const
323  {
324  if (!this->empty()) {
325 
326  const_iterator pos = lower_bound(this->begin(), this->end(), z, JLANG::make_comparator(&JVertex::getZ));
327 
328  if (pos != this->end() && pos != this->begin()) {
329 
330  --pos;
331 
332  return pos->getE() - (z - pos->getZ()) * gWater.getA();
333  }
334  }
335 
336  return MASS_MUON;
337  }
338  };
339 
340 
341  /**
342  * Get kinetic energy of particle with given mass.
343  *
344  * \param E energy [GeV]
345  * \param m mass [GeV]
346  * \return energy [GeV]
347  */
348  inline double getKineticEnergy(const double E, const double m)
349  {
350  if (E > m)
351  return sqrt((E - m) * (E + m));
352  else
353  return 0.0;
354  }
355 }
356 
357 #endif
JHitType_t
Enumeration of hit types based on km3 codes.
double getEs() const
Get shower energy.
JTrack(const JVertex &vertex)
Constructor.
Scattered light from primary shower.
double getT() const
Get time.
JVertex & step(const double ds)
Step.
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
Definition: JComparator.hh:143
static const double MASS_MUON
muon mass [GeV]
Definition: JConstants.hh:59
void applyEloss(const double Es)
Apply energy loss energy.
Direct light from delta-rays.
double getIndexOfRefraction()
Get average index of refraction of water.
Definition: JConstants.hh:111
Scattered light from muon.
const double getSpeedOfLight()
Number of bytes in a gigabyte.
Definition: JConstants.hh:89
Energy loss of muon.
scattered light from EM shower
Definition: JPDFTypes.hh:41
JHitType_t getHitType(const JPDFType_t pdf, const bool shower=false)
Get hit type corresponding to given PDF type.
double getE() const
Get muon energy.
Scattered light from Bremsstrahlung.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:328
Direct light from Bremsstrahlung.
Direct light from primary shower.
direct light from muon
Definition: JPDFTypes.hh:29
JVertex()
Default constructor.
JVertex(const double z, const double t, const double E)
Constructor.
Numbering scheme for PDF types.
JVertex & step(const JGeane &geane, const double ds)
Step.
scattered light from muon
Definition: JPDFTypes.hh:30
Constants.
double getZ() const
Get position.
double getE(const double z) const
Get muon energy at given position.
scattered light from delta-rays
Definition: JPDFTypes.hh:36
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
direct light from EM shower
Definition: JPDFTypes.hh:40
Muon trajectory.
JPDFType_t
PDF types.
Definition: JPDFTypes.hh:27
Interface for muon energy loss.
Definition: JGeane.hh:42
double getKineticEnergy(const double E, const double m)
Get kinetic energy of particle with given mass.
double getRange(const JGeane &geane) const
Get range of muon.
Vertex of energy loss of muon.
direct light from delta-rays
Definition: JPDFTypes.hh:35
int getNumberOfPhotoElectrons(const double NPE, const int N, const double npe)
Get number of photo-electrons of a hit given the expectation values of the number of photo-electrons ...
Direct light from muon.
Scattered light from delta-rays.
alias put_queue eval echo n
Definition: qlib.csh:19
double getRange() const
Get range of muon.
then usage $script[input file[working directory[option]]] nWhere option can be N
Definition: JMuonPostfit.sh:37
then usage $script[input file[working directory[option]]] nWhere option can be E
Definition: JMuonPostfit.sh:37
virtual double getA() const
Get energy loss constant.
Definition: JGeane.hh:234