Jpp 20.0.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JGeanz.hh
Go to the documentation of this file.
1#ifndef __JPHYSICS__JGEANZ__
2#define __JPHYSICS__JGEANZ__
3
4#include <istream>
5#include <ostream>
6#include <cmath>
7
9
10/**
11 * \file
12 * Longitudinal emission profile EM-shower.
13 * \author mdejong
14 */
15
16namespace JPHYSICS {}
17namespace JPP { using namespace JPHYSICS; }
18
19namespace JPHYSICS {
20
21 /**
22 * Function object for longitudinal profile of EM-shower.
23 *
24 *
25 * \f[P(z) \propto z^{a-1} \times e^{-z/b}\f]
26 *
27 * where: \f[a = a_{0} + a_{1} \times \ln(E)\f]
28 *
29 * The parametrisation is taken from reference:
30 * C. Kopper, "Performance Studies for the KM3NeT Neutrino Telescope.",
31 * PhD thesis, University of Erlangen.
32 */
33 class JGeanz {
34 public:
35 /**
36 * Constructor.
37 *
38 * \param __a0 power term (constant)
39 * \param __a1 power term (E dependence)
40 * \param __b expontial slope
41 */
42 JGeanz(const double __a0,
43 const double __a1,
44 const double __b) :
45 a0(__a0),
46 a1(__a1),
47 b (__b),
48 Emin(exp(-a0/a1))
49 {}
50
51
52 /**
53 * Probability Density Function
54 *
55 * \param E EM-shower energy [GeV]
56 * \param z z position of light emission point relative to vertex location (z >= 0) [m]
57 * \return dP/dz
58 */
59 double getProbability(const double E,
60 const double z) const
61 {
62 if (E > Emin) {
63
64 const double a = a0 + a1 * log(E);
65 const double y = pow(z,a-1.0) * exp(-z/b) / (pow(b,a) * std::tgamma(a));
66
67 return y;
68 }
69
70 if (z <= getMinimalShowerSize())
71 return 1.0 / getMinimalShowerSize();
72 else
73 return 0.0;
74 }
75
76
77 /**
78 * Probability Density Function
79 *
80 * \param E EM-shower energy [GeV]
81 * \param z z position of light emission point relative to vertex location (z >= 0) [m]
82 * \return dP/dz
83 */
84 double operator()(const double E,
85 const double z) const
86 {
87 return getProbability(E, z);
88 }
89
90
91 /**
92 * Integral of PDF (starting from 0).
93 *
94 * \param E EM-shower energy [GeV]
95 * \param z z position [m] (>= 0)
96 * \return dP
97 */
98 double getIntegral(const double E,
99 const double z) const
100 {
101 if (E > Emin) {
102
103 const double a = a0 + a1 * log(E);
104 const double x = z / b;
105 const double y = JMATH::Gamma(a,x);
106
107 return y;
108 }
109
110 if (z <= getMinimalShowerSize())
111 return z / getMinimalShowerSize();
112 else
113 return 1.0;
114 }
115
116
117 /**
118 * Derivative of PDF.
119 *
120 * \param E EM-shower energy [GeV]
121 * \param z z position [m] (> 0)
122 * \return dP/dz
123 */
124 double getDerivative(const double E,
125 const double z) const
126 {
127 if (E > Emin) {
128
129 const double a = a0 + a1 * log(E);
130
131 return ((a - 1.0)/z - 1.0/b) * getProbability(E, z);
132 }
133
134 return 0.0;
135 }
136
137
138 /**
139 * Get shower length for a given integrated probability.
140 *
141 * \param E EM-shower energy [GeV]
142 * \param P integrated probability [0,1]
143 * \param eps relative precision
144 * \return shower length [m]
145 */
146 double getLength(const double E,
147 const double P,
148 const double eps = 1.0e-5) const
149 {
150 if (E > Emin) {
151
152 double z0 = 0.0;
153 double z1 = getMaximum(E);
154
155 for (int i = 1000; i != 0; --i) {
156
157 const double p = getIntegral(E, z1);
158
159 if (fabs(p-P) < P*eps) {
160 return z1;
161 }
162
163 if (p > P) {
164 z1 = 0.5 * (z0 + z1);
165 } else {
166 z0 = z1;
167 z1 *= 1.5;
168 }
169 }
170
171 return z1;
172
173 } else {
174
175 return 0.0;
176 }
177 }
178
179
180 /**
181 * Get depth of shower maximum
182 *
183 * \param E EM-shower energy[GeV]
184 * \return depth of maximum [m]
185 */
186
187 double getMaximum(const double E) const
188 {
189 const double a = a0 + a1 * log(E);
190
191 return (a-1)*b;
192 }
193
194
195 /**
196 * Get minimal shower size.
197 *
198 * \return size [m]
199 */
200 static double getMinimalShowerSize()
201 {
202 return 1e-6;
203 }
204
205
206 /**
207 * Read longitudinal profile from input.
208 *
209 * \param in input stream
210 * \param object longitudinal profile
211 * \return input stream
212 */
213 friend inline std::istream& operator>>(std::istream& in, JGeanz& object)
214 {
215 in >> object.a0 >> object.a1 >> object.b;
216
217 object.Emin = exp(-object.a0/object.a1);
218
219 return in;
220 }
221
222
223 /**
224 * Write longitidunal profile to output.
225 *
226 * \param out output stream
227 * \param object longitudinal profile
228 * \return output stream
229 */
230 friend inline std::ostream& operator<<(std::ostream& out, const JGeanz& object)
231 {
232 return out << object.a0 << ' ' << object.a1 << ' ' << object.b;
233 }
234
235 protected:
236 double a0;
237 double a1;
238 double b;
239 double Emin;
240 };
241
242
243 /**
244 * Function object for longitudinal EM-shower profile
245 */
246 static const JGeanz geanz(1.85, 0.62, 0.54);
247}
248
249#endif
Auxiliary methods for mathematics.
Function object for longitudinal profile of EM-shower.
Definition JGeanz.hh:33
friend std::ostream & operator<<(std::ostream &out, const JGeanz &object)
Write longitidunal profile to output.
Definition JGeanz.hh:230
double getLength(const double E, const double P, const double eps=1.0e-5) const
Get shower length for a given integrated probability.
Definition JGeanz.hh:146
JGeanz(const double __a0, const double __a1, const double __b)
Constructor.
Definition JGeanz.hh:42
double getIntegral(const double E, const double z) const
Integral of PDF (starting from 0).
Definition JGeanz.hh:98
double getDerivative(const double E, const double z) const
Derivative of PDF.
Definition JGeanz.hh:124
double operator()(const double E, const double z) const
Probability Density Function.
Definition JGeanz.hh:84
friend std::istream & operator>>(std::istream &in, JGeanz &object)
Read longitudinal profile from input.
Definition JGeanz.hh:213
double getProbability(const double E, const double z) const
Probability Density Function.
Definition JGeanz.hh:59
static double getMinimalShowerSize()
Get minimal shower size.
Definition JGeanz.hh:200
double getMaximum(const double E) const
Get depth of shower maximum.
Definition JGeanz.hh:187
double Gamma(const double a, const double x)
Incomplete gamma function.
Auxiliary methods for light properties of deep-sea water.
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).