Jpp 19.3.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JPDF_t.hh
Go to the documentation of this file.
1#include "JLang/JException.hh"
3#include "JTools/JMap.hh"
4#include "JTools/JGridMap.hh"
5#include "JTools/JMapList.hh"
6#include "JTools/JSpline.hh"
7#include "JTools/JPolint.hh"
8#include "JTools/JElement.hh"
9#include "JTools/JResult.hh"
10#include "JPhysics/JPDFTable.hh"
12#include "JPhysics/JPDFTypes.hh"
13#include "JMath/JZero.hh"
14
15
16/**
17 * \file
18 *
19 * Auxiliary data structure for muon PDF.
20 * \author mdejong
21 */
22
23/**
24 * Auxiliary data structure for muon PDF.
25 */
26struct JPDF {
27
31
37
38
39 /**
40 * Constructor.
41 *
42 * The <tt>TTS</tt> corresponds to the additional time smearing applied to the PDFs.
43 *
44 * \param file_name file name
45 * \param TTS TTS [ns]
46 * \param numberOfPoints number of points for Gauss-Hermite integration of TTS
47 * \param epsilon precision for Gauss-Hermite integration of TTS
48 */
49 JPDF(const std::string& file_name,
50 const double TTS,
51 const int numberOfPoints = 25,
52 const double epsilon = 1.0e-10)
53 {
54 using namespace std;
55 using namespace JPP;
56
57 const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(zero));
58
59 cout << "loading input from file " << file_name << "... " << flush;
60
61 pdf.load(file_name.c_str());
62
63 pdf.setExceptionHandler(supervisor);
64
65 cout << "OK" << endl;
66
67 type = getPDFType(file_name);
68
69 if (TTS > 0.0) {
70
71 cout << "bluring PDFs... " << flush;
72
73 pdf.blur(TTS, numberOfPoints, epsilon);
74
75 cout << "OK" << endl;
76
77 } else if (TTS < 0.0) {
78
79 THROW(JValueOutOfRange, "Illegal value of TTS [ns]: " << TTS);
80 }
81 }
82
83
84 /**
85 * Get PDF.
86 *
87 * The orientation of the PMT should be defined according this <a href="https://common.pages.km3net.de/jpp/JPDF.PDF">documentation</a>.\n
88 * In this, the zenith and azimuth angles are limited to \f[\left[0, \pi\right]\f].
89 *
90 * \param E muon energy at minimum distance of approach [GeV]
91 * \param R minimum distance of approach [m]
92 * \param theta PMT zenith angle [rad]
93 * \param phi PMT azimuth angle [rad]
94 * \param t1 arrival time relative to Cherenkov hypothesis [ns]
95 * \return hypothesis value
96 */
97 result_type calculate(const double E,
98 const double R,
99 const double theta,
100 const double phi,
101 const double t1) const
102 {
103 using namespace JPP;
104
105 result_type h1 = pdf(R, theta, phi, t1);
106
107 if (is_bremsstrahlung(type)) {
108 h1 *= E;
109 } else if (is_deltarays(type)) {
110 h1 *= getDeltaRaysFromMuon(E);
111 }
112
113 // safety measures
114
115 if (h1.f <= 0.0) {
116 h1.f = 0.0;
117 h1.fp = 0.0;
118 }
119
120 if (h1.v <= 0.0) {
121 h1.v = 0.0;
122 }
123
124 return h1;
125 }
126
127 JPDF_t pdf; //!< PDF
128 int type; //!< PDF type
129};
130
131
132/**
133 * Auxiliary data structure for muon PDF.
134 */
136
140
146
147
148 /**
149 * Constructor.
150 *
151 * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD.\n
152 * The <tt>TTS</tt> corresponds to the additional time smearing applied to the PDFs.
153 *
154 * \param fileDescriptor PDF file descriptor
155 * \param TTS TTS [ns]
156 * \param numberOfPoints number of points for Gauss-Hermite integration of TTS
157 * \param epsilon precision for Gauss-Hermite integration of TTS
158 */
159 JMuonPDF_t(const std::string& fileDescriptor,
160 const double TTS,
161 const int numberOfPoints = 25,
162 const double epsilon = 1.0e-10)
163 {
164 using namespace std;
165 using namespace JPP;
166
167 const JPDFType_t pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
168 SCATTERED_LIGHT_FROM_MUON,
169 DIRECT_LIGHT_FROM_EMSHOWERS,
170 SCATTERED_LIGHT_FROM_EMSHOWERS,
171 DIRECT_LIGHT_FROM_DELTARAYS,
172 SCATTERED_LIGHT_FROM_DELTARAYS };
173
174 const int N = sizeof(pdf_t) / sizeof(pdf_t[0]);
175
176 JPDF_t pdf[N];
177
178 const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(zero));
179
180 for (int i = 0; i != N; ++i) {
181
182 const string file_name = getFilename(fileDescriptor, pdf_t[i]);
183
184 cout << "loading input from file " << file_name << "... " << flush;
185
186 pdf[i].load(file_name.c_str());
187
188 pdf[i].setExceptionHandler(supervisor);
189
190 cout << "OK" << endl;
191 }
192
193 // Add PDFs
194
195 cout << "adding PDFs... " << flush;
196
197 pdfA = pdf[1]; pdfA.add(pdf[0]);
198 pdfB = pdf[3]; pdfB.add(pdf[2]);
199 pdfC = pdf[5]; pdfC.add(pdf[4]);
200
201 cout << "OK" << endl;
202
203 if (TTS > 0.0) {
204
205 cout << "bluring PDFs... " << flush;
206
207 pdfA.blur(TTS, numberOfPoints, epsilon);
208 pdfB.blur(TTS, numberOfPoints, epsilon);
209 pdfC.blur(TTS, numberOfPoints, epsilon);
210
211 cout << "OK" << endl;
212
213 } else if (TTS < 0.0) {
214
215 THROW(JValueOutOfRange, "Illegal value of TTS [ns]: " << TTS);
216 }
217 }
218
219
220 /**
221 * Get PDF.
222 *
223 * The orientation of the PMT should be defined according this <a href="https://common.pages.km3net.de/jpp/JPDF.PDF">documentation</a>.\n
224 * In this, the zenith and azimuth angles are limited to \f[\left[0, \pi\right]\f].
225 *
226 * \param E muon energy at minimum distance of approach [GeV]
227 * \param R minimum distance of approach [m]
228 * \param theta PMT zenith angle [rad]
229 * \param phi PMT azimuth angle [rad]
230 * \param t1 arrival time relative to Cherenkov hypothesis [ns]
231 * \return hypothesis value
232 */
233 result_type calculate(const double E,
234 const double R,
235 const double theta,
236 const double phi,
237 const double t1) const
238 {
239 using namespace JPP;
240
241 result_type h1 = (pdfA(R, theta, phi, t1) +
242 pdfB(R, theta, phi, t1) * E +
243 pdfC(R, theta, phi, t1) * getDeltaRaysFromMuon(E));
244
245 // safety measures
246
247 if (h1.f <= 0.0) {
248 h1.f = 0.0;
249 h1.fp = 0.0;
250 }
251
252 if (h1.v <= 0.0) {
253 h1.v = 0.0;
254 }
255
256 return h1;
257 }
258
259 JPDF_t pdfA; //!< PDF for minimum ionising particle
260 JPDF_t pdfB; //!< PDF for average energy losses
261 JPDF_t pdfC; //!< PDF for delta-rays
262};
263
264
265/**
266 * Auxiliary data structure for shower PDF.
267 */
269
273
280
281
282 /**
283 * Constructor.
284 *
285 * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD.\n
286 * The <tt>TTS</tt> corresponds to the additional time smearing applied to the PDFs.
287 *
288 * \param fileDescriptor PDF file descriptor
289 * \param TTS TTS [ns]
290 * \param numberOfPoints number of points for Gauss-Hermite integration of TTS
291 * \param epsilon precision for Gauss-Hermite integration of TTS
292 */
293 JShowerPDF_t(const std::string& fileDescriptor,
294 const double TTS,
295 const int numberOfPoints = 25,
296 const double epsilon = 1.0e-10)
297 {
298 using namespace std;
299 using namespace JPP;
300
301 const JPDFType_t pdf_t[] = { SCATTERED_LIGHT_FROM_EMSHOWER,
302 DIRECT_LIGHT_FROM_EMSHOWER };
303
304 const int N = sizeof(pdf_t) / sizeof(pdf_t[0]);
305
306 const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(zero));
307
308 for (int i = 0; i != N; ++i) {
309
310 const string file_name = getFilename(fileDescriptor, pdf_t[i]);
311
312 cout << "loading input from file " << file_name << "... " << flush;
313
314 JPDF_t pdf;
315
316 pdf.load(file_name.c_str());
317
318 pdf.setExceptionHandler(supervisor);
319
320 if (pdfA.empty())
321 pdfA = pdf;
322 else
323 pdfA.add(pdf);
324
325 cout << "OK" << endl;
326 }
327
328 if (TTS > 0.0) {
329
330 cout << "bluring PDFs... " << flush;
331
332 pdfA.blur(TTS, numberOfPoints, epsilon);
333
334 cout << "OK" << endl;
335
336 } else if (TTS < 0.0) {
337
338 THROW(JValueOutOfRange, "Illegal value of TTS [ns]: " << TTS);
339 }
340 }
341
342
343 /**
344 * Get PDF.
345 *
346 * The orientation of the PMT should be defined according this <a href="https://common.pages.km3net.de/jpp/JPDF.PDF">documentation</a>.\n
347 * In this, the zenith and azimuth angles are limited to \f[\left[0, \pi\right]\f].
348 *
349 * \param E shower energy [GeV]
350 * \param D distance [m]
351 * \param cd cosine emission angle
352 * \param theta PMT zenith angle [rad]
353 * \param phi PMT azimuth angle [rad]
354 * \param t1 arrival time relative to Cherenkov hypothesis [ns]
355 * \return hypothesis value
356 */
357 result_type calculate(const double E,
358 const double D,
359 const double cd,
360 const double theta,
361 const double phi,
362 const double t1) const
363 {
364 using namespace JPP;
365
366 result_type h1 = pdfA(D, cd, theta, phi, t1) * E;
367
368 // safety measures
369
370 if (h1.f <= 0.0) {
371 h1.f = 0.0;
372 h1.fp = 0.0;
373 }
374
375 if (h1.v <= 0.0) {
376 h1.v = 0.0;
377 }
378
379 return h1;
380 }
381
382 JPDF_t pdfA; //!< PDF for shower
383};
General purpose class for a collection of sorted elements.
The elements in a collection are sorted according to their abscissa values and a given distance opera...
Exceptions.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Auxiliary methods for PDF calculations.
Numbering scheme for PDF types.
int numberOfPoints
Definition JResultPDF.cc:22
This include file containes various data structures that can be used as specific return types for the...
Definition of zero value for any class.
Exception for accessing a value in a collection that is outside of its range.
Multi-dimensional PDF table for arrival time of Cherenkov light.
Definition JPDFTable.hh:44
void blur(const double TTS, const int numberOfPoints=25, const double epsilon=1.0e-10, const double quantile=0.99)
Blur PDF.
Definition JPDFTable.hh:111
General purpose class for collection of elements, see: <a href="JTools.PDF";>Collection of elements....
Definition JSet.hh:22
void setExceptionHandler(const typename function_type::supervisor_type &supervisor)
Set the supervisor for handling of exceptions.
Template class for spline interpolation in 1D.
Definition JSpline.hh:734
function_type::result_type result_type
Definition JSpline.hh:752
void add(const JMultiFunction_t &input)
Add function.
JPDFType_t
PDF types.
Definition JPDFTypes.hh:24
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
void load(const char *file_name)
Load from input file.
Auxiliary data structure for muon PDF.
Definition JPDF_t.hh:135
JPP::JSplineFunction1D< JPP::JSplineElement2S< double, double >, JPP::JCollection, JPP::JResultPDF< double > > JFunction1D_t
Definition JPDF_t.hh:139
JPP::JMAPLIST< JPP::JPolint1FunctionalMap, JPP::JPolint0FunctionalGridMap, JPP::JPolint0FunctionalGridMap >::maplist JPDFMaplist_t
Definition JPDF_t.hh:143
JPP::JPDFTable< JFunction1D_t, JPDFMaplist_t > JPDF_t
Definition JPDF_t.hh:144
JPDF_t pdfB
PDF for average energy losses.
Definition JPDF_t.hh:260
JMuonPDF_t(const std::string &fileDescriptor, const double TTS, const int numberOfPoints=25, const double epsilon=1.0e-10)
Constructor.
Definition JPDF_t.hh:159
JPDF_t pdfC
PDF for delta-rays.
Definition JPDF_t.hh:261
JFunction1D_t::result_type result_type
Definition JPDF_t.hh:145
result_type calculate(const double E, const double R, const double theta, const double phi, const double t1) const
Get PDF.
Definition JPDF_t.hh:233
JPDF_t pdfA
PDF for minimum ionising particle.
Definition JPDF_t.hh:259
Auxiliary data structure for muon PDF.
Definition JPDF_t.hh:26
JPP::JMAPLIST< JPP::JPolint1FunctionalMap, JPP::JPolint0FunctionalGridMap, JPP::JPolint0FunctionalGridMap >::maplist JPDFMaplist_t
Definition JPDF_t.hh:34
JPDF_t pdf
PDF.
Definition JPDF_t.hh:127
JPP::JSplineFunction1D< JPP::JSplineElement2S< double, double >, JPP::JCollection, JPP::JResultPDF< double > > JFunction1D_t
Definition JPDF_t.hh:30
result_type calculate(const double E, const double R, const double theta, const double phi, const double t1) const
Get PDF.
Definition JPDF_t.hh:97
JPDF(const std::string &file_name, const double TTS, const int numberOfPoints=25, const double epsilon=1.0e-10)
Constructor.
Definition JPDF_t.hh:49
int type
PDF type.
Definition JPDF_t.hh:128
JPP::JPDFTable< JFunction1D_t, JPDFMaplist_t > JPDF_t
Definition JPDF_t.hh:35
JFunction1D_t::result_type result_type
Definition JPDF_t.hh:36
Auxiliary data structure for shower PDF.
Definition JPDF_t.hh:268
JShowerPDF_t(const std::string &fileDescriptor, const double TTS, const int numberOfPoints=25, const double epsilon=1.0e-10)
Constructor.
Definition JPDF_t.hh:293
JFunction1D_t::result_type result_type
Definition JPDF_t.hh:279
JPP::JMAPLIST< JPP::JPolint1FunctionalMap, JPP::JPolint1FunctionalMap, JPP::JPolint0FunctionalGridMap, JPP::JPolint0FunctionalGridMap >::maplist JPDFMaplist_t
Definition JPDF_t.hh:277
JPDF_t pdfA
PDF for shower.
Definition JPDF_t.hh:382
result_type calculate(const double E, const double D, const double cd, const double theta, const double phi, const double t1) const
Get PDF.
Definition JPDF_t.hh:357
JPP::JSplineFunction1D< JPP::JSplineElement2S< double, double >, JPP::JCollection, JPP::JResultPDF< double > > JFunction1D_t
Definition JPDF_t.hh:272
JPP::JPDFTable< JFunction1D_t, JPDFMaplist_t > JPDF_t
Definition JPDF_t.hh:278
Auxiliary class for recursive map list generation.
Definition JMapList.hh:109
Type definition of a zero degree polynomial interpolation based on a JGridMap implementation.
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
Data structure for result including value, first derivative and integrals of function.
Definition JResult.hh:339