Jpp  test_elongated_shower_pde
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JPDF_t.hh
Go to the documentation of this file.
1 #include "JLang/JException.hh"
2 #include "JTools/JCollection.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"
11 #include "JPhysics/JPDFToolkit.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  */
26 struct JPDF {
27 
28  typedef JPP::JSplineFunction1D<JPP::JSplineElement2S<double, double>,
29  JPP::JCollection,
30  JPP::JResultPDF<double> > JFunction1D_t;
31 
32  typedef JPP::JMAPLIST<JPP::JPolint1FunctionalMap,
33  JPP::JPolint0FunctionalGridMap,
34  JPP::JPolint0FunctionalGridMap>::maplist JPDFMaplist_t;
35  typedef JPP::JPDFTable<JFunction1D_t, JPDFMaplist_t> JPDF_t;
36  typedef JFunction1D_t::result_type result_type;
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  */
135 struct JMuonPDF_t {
136 
137  typedef JPP::JSplineFunction1D<JPP::JSplineElement2S<double, double>,
138  JPP::JCollection,
139  JPP::JResultPDF<double> > JFunction1D_t;
140 
141  typedef JPP::JMAPLIST<JPP::JPolint1FunctionalMap,
142  JPP::JPolint0FunctionalGridMap,
143  JPP::JPolint0FunctionalGridMap>::maplist JPDFMaplist_t;
144  typedef JPP::JPDFTable<JFunction1D_t, JPDFMaplist_t> JPDF_t;
145  typedef JFunction1D_t::result_type result_type;
146 
147 
148  /**
149  * Constructor.
150  *
151  * The PDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD.\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  */
268 struct JShowerPDF_t {
269 
270  typedef JPP::JSplineFunction1D<JPP::JSplineElement2S<double, double>,
271  JPP::JCollection,
272  JPP::JResultPDF<double> > JFunction1D_t;
273 
274  typedef JPP::JMAPLIST<JPP::JPolint1FunctionalMap,
275  JPP::JPolint1FunctionalMap,
276  JPP::JPolint0FunctionalGridMap,
277  JPP::JPolint0FunctionalGridMap>::maplist JPDFMaplist_t;
278  typedef JPP::JPDFTable<JFunction1D_t, JPDFMaplist_t> JPDF_t;
279  typedef JFunction1D_t::result_type result_type;
280 
281 
282  /**
283  * Constructor.
284  *
285  * The PDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD.\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 };
JIO::JObjectBinaryIO::load
void load(const char *file_name)
Load from input file.
Definition: JObjectBinaryIO.hh:29
JException.hh
Exceptions.
JMuonPDF_t::JPDF_t
JPP::JPDFTable< JFunction1D_t, JPDFMaplist_t > JPDF_t
Definition: JPDF_t.hh:144
JShowerPDF_t::JPDFMaplist_t
JPP::JMAPLIST< JPP::JPolint1FunctionalMap, JPP::JPolint1FunctionalMap, JPP::JPolint0FunctionalGridMap, JPP::JPolint0FunctionalGridMap >::maplist JPDFMaplist_t
Definition: JPDF_t.hh:277
JPDFToolkit.hh
Auxiliary methods for PDF calculations.
JTOOLS::JCollection
General purpose class for collection of elements, see: &lt;a href=&quot;JTools.PDF&quot;;&gt;Collection of elements...
Definition: JCollection.hh:73
E
then usage E
Definition: JMuonPostfit.sh:35
JElement.hh
The elements in a collection are sorted according to their abscissa values and a given distance opera...
JPDF::JPDF
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
JResult.hh
This include file containes various data structures that can be used as specific return types for the...
JPDF::pdf
JPDF_t pdf
PDF.
Definition: JPDF_t.hh:127
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER
scattered light from EM shower
Definition: JPDFTypes.hh:41
JMuonPDF_t::JFunction1D_t
JPP::JSplineFunction1D< JPP::JSplineElement2S< double, double >, JPP::JCollection, JPP::JResultPDF< double > > JFunction1D_t
Definition: JPDF_t.hh:139
THROW
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:696
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWERS
direct light from EM showers
Definition: JPDFTypes.hh:32
N
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
Definition: JShowerPostfit.sh:95
JPHYSICS::JPDFTable::blur
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
JShowerPDF_t::pdfA
JPDF_t pdfA
PDF for shower.
Definition: JPDF_t.hh:382
h1
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
JMATH::zero
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
JPDF::JPDFMaplist_t
JPP::JMAPLIST< JPP::JPolint1FunctionalMap, JPP::JPolint0FunctionalGridMap, JPP::JPolint0FunctionalGridMap >::maplist JPDFMaplist_t
Definition: JPDF_t.hh:34
JZero.hh
Definition of zero value for any class.
JShowerPDF_t
Auxiliary data structure for shower PDF.
Definition: JPDF_t.hh:268
JMuonPDF_t::JPDFMaplist_t
JPP::JMAPLIST< JPP::JPolint1FunctionalMap, JPP::JPolint0FunctionalGridMap, JPP::JPolint0FunctionalGridMap >::maplist JPDFMaplist_t
Definition: JPDF_t.hh:143
JPHYSICS::DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition: JPDFTypes.hh:29
JMuonPDF_t::JMuonPDF_t
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
JPDFTypes.hh
Numbering scheme for PDF types.
JTOOLS::JSplineFunction1D
Template class for spline interpolation in 1D.
Definition: JSpline.hh:666
JMuonPDF_t::calculate
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
JPHYSICS::is_bremsstrahlung
bool is_bremsstrahlung(const int pdf)
Test if given PDF type corresponds to Cherenkov light from Bremsstrahlung.
Definition: JPDFTypes.hh:165
JPDF::result_type
JFunction1D_t::result_type result_type
Definition: JPDF_t.hh:36
JPHYSICS::SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition: JPDFTypes.hh:30
JShowerPDF_t::result_type
JFunction1D_t::result_type result_type
Definition: JPDF_t.hh:279
JMuonPDF_t::result_type
JFunction1D_t::result_type result_type
Definition: JPDF_t.hh:145
JTOOLS::JPolint1FunctionalMap
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
Definition: JFunctionalMap_t.hh:55
JMuonPDF_t::pdfB
JPDF_t pdfB
PDF for average energy losses.
Definition: JPDF_t.hh:260
JTOOLS::JResultPDF
Data structure for result including value, first derivative and integrals of function.
Definition: JResult.hh:337
JTOOLS::JTransformableMultiFunction::add
void add(const JMultiFunction_t &input)
Add function.
Definition: JTransformableMultiFunction.hh:155
JTOOLS::JSplineFunction1D::result_type
function_type::result_type result_type
Definition: JSpline.hh:687
JShowerPDF_t::calculate
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
JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition: JPDFTypes.hh:36
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER
direct light from EM shower
Definition: JPDFTypes.hh:40
JPDF::JFunction1D_t
JPP::JSplineFunction1D< JPP::JSplineElement2S< double, double >, JPP::JCollection, JPP::JResultPDF< double > > JFunction1D_t
Definition: JPDF_t.hh:30
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWERS
scattered light from EM showers
Definition: JPDFTypes.hh:33
JCollection.hh
General purpose class for a collection of sorted elements.
JPHYSICS::JPDFType_t
JPDFType_t
PDF types.
Definition: JPDFTypes.hh:27
JTOOLS::JMAPLIST
Auxiliary class for recursive map list generation.
Definition: JMapList.hh:108
R
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:43
JPDF::calculate
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
JMuonPDF_t
Auxiliary data structure for muon PDF.
Definition: JPDF_t.hh:135
JMuonPDF_t::pdfC
JPDF_t pdfC
PDF for delta-rays.
Definition: JPDF_t.hh:261
JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition: JPDFTypes.hh:35
JPHYSICS::getPDFType
int getPDFType(const std::string &file_name)
Get PDF type.
Definition: JPDFTypes.hh:83
JTOOLS::JMultiFunction::setExceptionHandler
void setExceptionHandler(const typename function_type::supervisor_type &supervisor)
Set the supervisor for handling of exceptions.
Definition: JMultiFunction.hh:157
JPDF::JPDF_t
JPP::JPDFTable< JFunction1D_t, JPDFMaplist_t > JPDF_t
Definition: JPDF_t.hh:35
JTOOLS::JPolint0FunctionalGridMap
Type definition of a zero degree polynomial interpolation based on a JGridMap implementation.
Definition: JFunctionalMap_t.hh:82
numberOfPoints
int numberOfPoints
Definition: JResultPDF.cc:22
JPHYSICS::is_deltarays
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition: JPDFTypes.hh:179
JPHYSICS::getDeltaRaysFromMuon
double getDeltaRaysFromMuon(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:67
JEEP::getFilename
std::string getFilename(const std::string &file_name)
Get file name part, i.e. part after last JEEP::PATHNAME_SEPARATOR if any.
Definition: JeepToolkit.hh:88
JPDF::type
int type
PDF type.
Definition: JPDF_t.hh:128
JShowerPDF_t::JShowerPDF_t
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
JMuonPDF_t::pdfA
JPDF_t pdfA
PDF for minimum ionising particle.
Definition: JPDF_t.hh:259
JShowerPDF_t::JFunction1D_t
JPP::JSplineFunction1D< JPP::JSplineElement2S< double, double >, JPP::JCollection, JPP::JResultPDF< double > > JFunction1D_t
Definition: JPDF_t.hh:272
D
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
JShowerPDF_t::JPDF_t
JPP::JPDFTable< JFunction1D_t, JPDFMaplist_t > JPDF_t
Definition: JPDF_t.hh:278
JPDF
Auxiliary data structure for muon PDF.
Definition: JPDF_t.hh:26
Generated by   doxygen 1.8.5