Jpp  test_elongated_shower_pde
the software that should make you happy
 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 "JPhysics/JConstants.hh"
10 #include "JPhysics/JPDFTypes.hh"
11 #include "JPhysics/JGeane.hh"
12 #include "JAAnet/JAAnetToolkit.hh"
14 #include "JGeometry3D/JTime.hh"
15 #include "JGeometry3D/JVersor3Z.hh"
16 #include "JGeometry2D/JVector2D.hh"
17 #include "JLang/JComparator.hh"
18 
19 
20 /**
21  * \author mdejong
22  */
23 
24 
25 /**
26  * Detector simulations.
27  */
28 namespace JSIRENE {}
29 namespace JPP { using namespace JSIRENE; }
30 
31 namespace JSIRENE {
32 
34  using JPHYSICS::MASS_MUON;
36  using JPHYSICS::JGeane;
37  using JPHYSICS::gWater;
40  using JGEOMETRY3D::JTime;
43  using JAANET::JHitType_t;
44 
45 
46  /**
47  * Get number of photo-electrons of a hit given the expectation
48  * values of the number of photo-electrons on a module and PMT.
49  *
50  * The return value is evaluated by pick-and-drop statistics from
51  * the generated number of photo-electrons when the expectation
52  * value of the number of photo-electrons on a module deviates
53  * less than 5 sigmas from 0 (i.e.\ when it is less than 25).
54  * Otherwise, the return value is evaluated by Poisson statistics
55  * from the expectation value of the number of photo-electrons on PMT.
56  *
57  * \param NPE expectation value of npe on module
58  * \param N generated value of npe on module
59  * \param npe expectation value of npe on PMT
60  * \return number of photo-electrons on PMT
61  */
62  inline int getNumberOfPhotoElectrons(const double NPE,
63  const int N,
64  const double npe)
65  {
66  if (NPE < 25.0) {
67 
68  int n = 0;
69 
70  for (int i = N; i != 0; --i) {
71  if (gRandom->Rndm() * NPE < npe) {
72  ++n;
73  }
74  }
75 
76  return n;
77 
78  } else {
79 
80  return gRandom->Poisson(npe);
81  }
82  }
83 
84 
85  /**
86  * Get number of photo-electrons of a hit given number of photo-electrons on PMT.
87  *
88  * The number of photo-electrons of a hit may be larger than unity
89  * to limit the overall number of hits and consequently the memory coonsumption as well as
90  * the number of times the arrival time needs to be evaluated which is CPU intensive.
91  *
92  * \param npe number of photo-electrons on PMT
93  * \return number of photo-electrons of hit
94  */
95  inline int getNumberOfPhotoElectrons(const int npe)
96  {
97  static const int NPE = 20;
98 
99  if (npe > NPE) {
100 
101  const int n = (int) (NPE * log10((double) npe / (double) NPE));
102 
103  if (n == 0)
104  return 1;
105  else
106  return n;
107  }
108 
109  return 1;
110  }
111 
112 
113  /**
114  * Get hit type corresponding to given PDF type.
115  *
116  * \param pdf PDF type
117  * \param shower force origin from neutrino interaction shower
118  * \return hit type
119  */
120  inline JHitType_t getHitType(const JPDFType_t pdf, const bool shower = false)
121  {
122  using namespace JAANET;
123  using namespace JPHYSICS;
124 
125  switch (pdf) {
126 
128  return HIT_TYPE_MUON_DIRECT;
129 
132 
135 
138 
141 
144 
145  default:
146  return HIT_TYPE_UNKNOWN;
147  }
148  }
149 
150 
151  /**
152  * Point along muon trajectory.
153  */
154  struct JPoint :
155  JPosition3D,
156  JTime
157  {
158  /**
159  * Default constructor.
160  */
161  JPoint() :
162  JPosition3D(),
163  JTime(),
164  __E (0.0),
165  __Es(0.0)
166  {}
167 
168 
169  /**
170  * Constructor.
171  *
172  * \param z position [m]
173  * \param t time [ns]
174  * \param E energy [GeV]
175  */
176  JPoint(const double z,
177  const double t,
178  const double E) :
179  JPosition3D(0.0, 0.0, z),
180  JTime(t),
181  __E (E),
182  __Es(0.0)
183  {}
184 
185 
186  /**
187  * Get muon energy.
188  *
189  * \return energy [GeV]
190  */
191  double getE() const
192  {
193  return __E;
194  }
195 
196 
197  /**
198  * Get shower energy.
199  *
200  * \return energy [GeV]
201  */
202  double getEs() const
203  {
204  return __Es;
205  }
206 
207  protected:
208  double __E;
209  double __Es;
210  };
211 
212 
213  /**
214  * Muon trajectory.
215  */
216  struct JTrack :
217  public std::vector<JPoint>
218  {
219  /**
220  * Constructor.
221  *
222  * \param point muon starting point
223  */
224  JTrack(const JPoint& point) :
225  std::vector<JPoint>(1, point)
226  {}
227 
228 
229  /**
230  * Get muon energy at given position along trajectory.
231  *
232  * \param z position [m]
233  * \return energy [GeV]
234  */
235  double getE(const double z) const
236  {
237  if (!this->empty()) {
238 
239  const_iterator p = lower_bound(this->begin(), this->end(), z, JLANG::make_comparator(&JPoint::getZ));
240 
241  if (p != this->end() && p != this->begin()) {
242 
243  --p;
244 
245  return p->getE() - (z - p->getZ()) * gWater.getA();
246  }
247  }
248 
249  return MASS_MUON;
250  }
251 
252 
253  /**
254  * Get muon position at given position along trajectory.
255  *
256  * \param z position [m]
257  * \return position
258  */
259  JVector2D getPosition(const double z) const
260  {
261  using namespace JPP;
262 
263  const double precision = 1.0e-2;
264 
265  if (!this->empty()) {
266 
267  const_iterator p = lower_bound(this->begin(), this->end(), z, JLANG::make_comparator(&JPoint::getZ));
268 
269  if (p == this->end()) {
270  --p;
271  }
272 
273  if (p == this->begin()) {
274 
275  return JVector2D(p->getX(),
276  p->getY());
277 
278  } else {
279 
280  JVector3D pos;
281 
282  pos = p->getPosition();
283 
284  --p;
285 
286  pos -= p->getPosition();
287 
288  const double u = (pos.getZ() > precision ? (z - p->getZ()) / pos.getZ() : 0.0);
289 
290  return JVector2D(p->getX() + u * pos.getX(),
291  p->getY() + u * pos.getY());
292  }
293  }
294 
295  return JVector2D(0.0, 0.0);
296  }
297  };
298 
299 
300  /**
301  * Vertex of energy loss of muon.
302  */
303  struct JVertex :
304  JPoint,
305  JVersor3Z
306  {
307  /**
308  * Default constructor.
309  */
311  JPoint(),
312  JVersor3Z()
313  {}
314 
315 
316  /**
317  * Constructor.
318  *
319  * \param z position [m]
320  * \param t time [ns]
321  * \param E energy [GeV]
322  */
323  JVertex(const double z,
324  const double t,
325  const double E) :
326  JPoint(z, t, E),
327  JVersor3Z()
328  {}
329 
330 
331  /**
332  * Get visible range of muon using default ionisation energy loss.
333  * This method applies only ionisation energy loss.
334  *
335  * \return range [m]
336  */
337  double getRange() const
338  {
339  return getRange(gWater.getA());
340  }
341 
342 
343  /**
344  * Get visible range of muon using given ionisation energy loss.
345  * This method applies only ionisation energy loss.
346  *
347  * \param a ionization parameter [GeV/m]
348  * \return range [m]
349  */
350  double getRange(double a) const
351  {
352  if (__E > MASS_MUON * getIndexOfRefraction())
353  return (__E - MASS_MUON * getIndexOfRefraction()) / a;
354  else
355  return 0.0;
356  }
357 
358 
359  /**
360  * Get range of muon using given energy loss function.
361  *
362  * \param geane energy loss function
363  * \return range [m]
364  */
365  double getRange(const JGeane& geane) const
366  {
367  return geane(__E);
368  }
369 
370 
371  /**
372  * Apply energy loss.
373  *
374  * \param Ts scattering angles
375  * \param Es energy [GeV]
376  */
377  void applyEloss(const JVersor3Z& Ts, const double Es)
378  {
379  static_cast<JVersor3Z&>(*this).add(Ts);
380 
381  this->__E -= Es;
382  this->__Es = Es;
383  }
384 
385 
386  /**
387  * Step using default ionisation energy loss.
388  * This method applies only ionisation energy loss.
389  *
390  * \param ds step [m]
391  * \return this vertex
392  */
393  JVertex& step(const double ds)
394  {
395  return step(gWater.getA(), ds);
396  }
397 
398 
399  /**
400  * Step using given ionisation energy loss.
401  * This method applies only given ionisation energy loss.
402  *
403  * \param a ionization parameter [GeV/m]
404  * \param ds step [m]
405  * \return this vertex
406  */
407  JVertex& step(const double a, const double ds)
408  {
409  __x += this->getDX() * ds;
410  __y += this->getDY() * ds;
411  __z += this->getDZ() * ds;
412  __t += ds / getSpeedOfLight();
413  __E -= ds * a;
414 
415  return *this;
416  }
417 
418 
419  /**
420  * Step using given energy loss function.
421  *
422  * \param geane energy loss function
423  * \param ds step [m]
424  * \return this vertex
425  */
426  JVertex& step(const JGeane& geane, const double ds)
427  {
428  __x += this->getDX() * ds;
429  __y += this->getDY() * ds;
430  __z += this->getDZ() * ds;
431  __t += ds / getSpeedOfLight();
432  __E = geane(__E, ds);
433 
434  return *this;
435  }
436  };
437 
438 
439  /**
440  * Get kinetic energy of particle with given mass.
441  *
442  * \param E energy [GeV]
443  * \param m mass [GeV]
444  * \return energy [GeV]
445  */
446  inline double getKineticEnergy(const double E, const double m)
447  {
448  if (E > m)
449  return sqrt((E - m) * (E + m));
450  else
451  return 0.0;
452  }
453 }
454 
455 #endif
Data structure for vector in two dimensions.
Definition: JVector2D.hh:32
JHitType_t
Enumeration of hit types based on km3 codes.
Scattered light from primary shower.
Point along muon trajectory.
double getE() const
Get muon energy.
JVertex & step(const double ds)
Step using default ionisation energy loss.
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:185
then usage E
Definition: JMuonPostfit.sh:35
double getRange(double a) const
Get visible range of muon using given ionisation energy loss.
Direct light from delta-rays.
Scattered light from muon.
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
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.
static const double MASS_MUON
muon mass [GeV]
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
Scattered light from Bremsstrahlung.
JTrack(const JPoint &point)
Constructor.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:325
Direct light from Bremsstrahlung.
Direct light from primary shower.
JVersor3Z & add(const JVersor3Z &value)
Addition operator.
Definition: JVersor3Z.hh:200
JVertex & step(const double a, const double ds)
Step using given ionisation energy loss.
direct light from muon
Definition: JPDFTypes.hh:29
JPoint()
Default constructor.
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 using given energy loss function.
const int n
Definition: JPolint.hh:676
scattered light from muon
Definition: JPDFTypes.hh:30
set_variable E_E log10(E_{fit}/E_{#mu})"
double getE(const double z) const
Get muon energy at given position along trajectory.
scattered light from delta-rays
Definition: JPDFTypes.hh:36
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
Physics constants.
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:39
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 using given energy loss function.
double getEs() const
Get shower energy.
Vertex of energy loss of muon.
direct light from delta-rays
Definition: JPDFTypes.hh:35
void applyEloss(const JVersor3Z &Ts, const double Es)
Apply energy loss.
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 ...
const double getSpeedOfLight()
Get speed of light.
Direct light from muon.
then JCalibrateToT a
Definition: JTuneHV.sh:116
JVector2D getPosition(const double z) const
Get muon position at given position along trajectory.
Scattered light from delta-rays.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
JPoint(const double z, const double t, const double E)
Constructor.
double u[N+1]
Definition: JPolint.hh:755
double getRange() const
Get visible range of muon using default ionisation energy loss.
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:39
double getZ() const
Get z position.
Definition: JVector3D.hh:115
virtual double getA() const override
Get energy loss constant.
Definition: JGeane.hh:231