Jpp  debug
the software that should make you happy
Functions | Variables
JMakePDF.cc File Reference

Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a muon. More...

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <set>
#include <map>
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JPDF.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/Antares.hh"
#include "JPhysics/KM3NeT.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

double getAbsorptionLength (const double lambda)
 
double getScatteringLength (const double lambda)
 
int main (int argc, char **argv)
 

Variables

double absorptionLengthFactor
 Scaling of absorption and scattering length. More...
 
double scatteringLengthFactor
 

Detailed Description

Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a muon.

The PDFs are tabulated as a function of (R, theta, phi, t), where:

The orientation of the PMT is defined in the coordinate system in which the muon travels along the z-axis and the PMT is located in the x-z plane.

Author
mdejong

Definition in file JMakePDF.cc.

Function Documentation

◆ getAbsorptionLength()

double getAbsorptionLength ( const double  lambda)
inline

Definition at line 26 of file JMakePDF.cc.

27 {
28  return absorptionLengthFactor * NAMESPACE::getAbsorptionLength(lambda);
29 }
getAbsorptionLength
double getAbsorptionLength(const double lambda)
Definition: JMakePDF.cc:26
absorptionLengthFactor
double absorptionLengthFactor
Scaling of absorption and scattering length.
Definition: JMakePDF.cc:22

◆ getScatteringLength()

double getScatteringLength ( const double  lambda)
inline

Definition at line 32 of file JMakePDF.cc.

33 {
34  return scatteringLengthFactor * NAMESPACE::getScatteringLength(lambda);
35 }
getScatteringLength
double getScatteringLength(const double lambda)
Definition: JMakePDF.cc:32
scatteringLengthFactor
double scatteringLengthFactor
Definition: JMakePDF.cc:23

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 52 of file JMakePDF.cc.

53 {
54  using namespace std;
55  using namespace JPP;
56 
57  string outputFile;
58  int numberOfPoints;
59  double epsilon;
60  int function;
61  set<double> R; // distance [m]
62  int debug;
63 
64  try {
65 
66  JParser<> zap("Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a muon.");
67 
68  zap['o'] = make_field(outputFile);
69  zap['n'] = make_field(numberOfPoints, "points for integration") = 25;
70  zap['e'] = make_field(epsilon, "precision for integration") = 1.0e-10;
71  zap['A'] = make_field(absorptionLengthFactor, "scaling factor") = 1.0;
72  zap['S'] = make_field(scatteringLengthFactor, "scaling factor") = 1.0;
73  zap['F'] = make_field(function, "PDF type") =
74  DIRECT_LIGHT_FROM_MUON,
75  DIRECT_LIGHT_FROM_EMSHOWERS,
76  DIRECT_LIGHT_FROM_DELTARAYS,
77  SCATTERED_LIGHT_FROM_MUON,
78  SCATTERED_LIGHT_FROM_EMSHOWERS,
79  SCATTERED_LIGHT_FROM_DELTARAYS;
80  zap['R'] = make_field(R, "distance of approach [m]") = JPARSER::initialised();
81  zap['d'] = make_field(debug) = 0;
82 
83  zap['F'] = JPARSER::not_initialised();
84 
85  zap(argc, argv);
86  }
87  catch(const exception &error) {
88  FATAL(error.what() << endl);
89  }
90 
91 
92  typedef double (JPDF::*fcn)(const double,
93  const double,
94  const double,
95  const double) const;
96 
97 
98  // set global parameters
99 
100  const double P_atm = NAMESPACE::getAmbientPressure();
101  const double wmin = getMinimalWavelength();
102  const double wmax = getMaximalWavelength();
103 
104 
105  const JPDF_C
106  pdf_c(NAMESPACE::getPhotocathodeArea(),
107  NAMESPACE::getQE,
108  NAMESPACE::getAngularAcceptance,
109  getAbsorptionLength,
110  getScatteringLength,
111  NAMESPACE::getScatteringProbability,
112  P_atm,
113  wmin,
114  wmax,
115  numberOfPoints,
116  epsilon);
117 
118 
119  typedef JSplineFunction1D_t JFunction1D_t;
120  typedef JMAPLIST<JPolint1FunctionalMap,
121  JPolint1FunctionalGridMap,
122  JPolint1FunctionalGridMap>::maplist JMapList_t;
123  typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
124 
125  typedef JPDFTransformer<3, JFunction1D_t::argument_type> JFunction3DTransformer_t;
126  typedef JArray<3, JFunction1D_t::argument_type> JArray_t;
127 
128  JPDF_t pdf;
129 
130 
131  NOTICE("building multi-dimensional function object <" << function << ">... " << flush);
132 
133 
134  const double kmin = pdf_c.getKappa(wmax);
135  const double kmax = pdf_c.getKappa(wmin);
136  const double cmin = pdf_c.getKmin (wmax);
137 
138  map<int, pair<fcn,JFunction3DTransformer_t> > zmap;
139 
140  zmap[DIRECT_LIGHT_FROM_MUON] = make_pair((fcn) &JPDF::getDirectLightFromMuon, JFunction3DTransformer_t(21.5, 2, kmin, kmax, NAMESPACE::getAngularAcceptance, 0.001));
141  zmap[SCATTERED_LIGHT_FROM_MUON] = make_pair((fcn) &JPDF::getScatteredLightFromMuon, JFunction3DTransformer_t(35.0, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
142  zmap[DIRECT_LIGHT_FROM_EMSHOWERS] = make_pair((fcn) &JPDF::getDirectLightFromEMshowers, JFunction3DTransformer_t(21.5, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
143  zmap[SCATTERED_LIGHT_FROM_EMSHOWERS] = make_pair((fcn) &JPDF::getScatteredLightFromEMshowers, JFunction3DTransformer_t(35.0, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
144  zmap[DIRECT_LIGHT_FROM_DELTARAYS] = make_pair((fcn) &JPDF::getDirectLightFromDeltaRays, JFunction3DTransformer_t(21.5, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
145  zmap[SCATTERED_LIGHT_FROM_DELTARAYS] = make_pair((fcn) &JPDF::getScatteredLightFromDeltaRays, JFunction3DTransformer_t(35.0, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
146 
147  if (zmap.find(function) == zmap.end()) {
148  FATAL("illegal function specifier" << endl);
149  }
150 
151  fcn f = zmap[function].first; // PDF
152  JFunction3DTransformer_t transformer = zmap[function].second; // transformer
153 
154 
155  if (R.empty()) {
156  R.insert( 0.10);
157  R.insert( 0.30);
158  R.insert( 0.50);
159  R.insert( 1.00);
160  R.insert( 2.00);
161  R.insert( 3.00);
162  R.insert( 4.00);
163  R.insert( 5.00);
164  R.insert( 6.00);
165  R.insert( 7.00);
166  R.insert( 8.00);
167  R.insert( 9.00);
168  R.insert( 10.00);
169  R.insert( 11.00);
170  R.insert( 12.00);
171  R.insert( 13.00);
172  R.insert( 14.00);
173  R.insert( 15.00);
174  R.insert( 16.00);
175  R.insert( 17.00);
176  R.insert( 18.00);
177  R.insert( 19.00);
178  R.insert( 20.00);
179  R.insert( 22.00);
180  R.insert( 24.00);
181  R.insert( 26.00);
182  R.insert( 28.00);
183  R.insert( 30.00);
184  R.insert( 32.00);
185  R.insert( 34.00);
186  R.insert( 36.00);
187  R.insert( 38.00);
188  R.insert( 40.00);
189  R.insert( 42.00);
190  R.insert( 44.00);
191  R.insert( 46.00);
192  R.insert( 48.00);
193  R.insert( 50.00);
194  R.insert( 55.00);
195  R.insert( 60.00);
196  R.insert( 65.00);
197  R.insert( 70.00);
198  R.insert( 75.00);
199  R.insert( 80.00);
200  R.insert( 85.00);
201  R.insert( 90.00);
202  R.insert( 95.00);
203  R.insert(100.00);
204  R.insert(110.00);
205  R.insert(120.00);
206  R.insert(130.00);
207  R.insert(140.00);
208  R.insert(150.00);
209  R.insert(170.00);
210  R.insert(190.00);
211  R.insert(210.00);
212  R.insert(250.00);
213  }
214 
215  set<double> X;
216 
217  if (function == DIRECT_LIGHT_FROM_MUON) {
218 
219  for (double buffer[] = { -0.01, -0.005, 0.0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, -1.0 }, *x = buffer; *x != -1.0; ++x) {
220  X.insert(0.0 + *x);
221  X.insert(1.0 - *x);
222  }
223 
224  for (double x = 0.01; x < 0.1; x += 0.0025) {
225  X.insert(0.0 + x);
226  X.insert(1.0 - x);
227  }
228 
229  for (double x = 0.10; x < 0.5; x += 0.010) {
230  X.insert(0.0 + x);
231  X.insert(1.0 - x);
232  }
233 
234  } else {
235 
236  X.insert( 0.00);
237  X.insert( 0.01);
238  X.insert( 0.02);
239  X.insert( 0.03);
240  X.insert( 0.04);
241  X.insert( 0.05);
242  X.insert( 0.06);
243  X.insert( 0.07);
244  X.insert( 0.08);
245  X.insert( 0.09);
246  X.insert( 0.10);
247  X.insert( 0.12);
248  X.insert( 0.15);
249  X.insert( 0.20);
250  X.insert( 0.25);
251  X.insert( 0.30);
252  X.insert( 0.40);
253  X.insert( 0.50);
254  X.insert( 0.60);
255  X.insert( 0.70);
256  X.insert( 0.80);
257  X.insert( 0.90);
258  X.insert( 1.00);
259  X.insert( 1.10);
260  X.insert( 1.20);
261  X.insert( 1.30);
262  X.insert( 1.40);
263  X.insert( 1.50);
264  X.insert( 1.60);
265  X.insert( 1.70);
266  X.insert( 1.80);
267  X.insert( 1.90);
268  X.insert( 2.00);
269  X.insert( 2.20);
270  X.insert( 2.40);
271  X.insert( 2.60);
272  X.insert( 2.80);
273  X.insert( 3.00);
274  X.insert( 3.25);
275  X.insert( 3.50);
276  X.insert( 3.75);
277  X.insert( 4.00);
278  X.insert( 4.25);
279  X.insert( 4.50);
280  X.insert( 4.75);
281  X.insert( 5.0);
282  X.insert( 6.0);
283  X.insert( 7.0);
284  X.insert( 8.0);
285  X.insert( 9.0);
286  X.insert( 10.0);
287  X.insert( 12.0);
288  X.insert( 14.0);
289  X.insert( 16.0);
290  X.insert( 18.0);
291  X.insert( 20.0);
292  X.insert( 25.0);
293  X.insert( 30.0);
294  X.insert( 40.0);
295  X.insert( 50.0);
296  X.insert( 60.0);
297  X.insert( 70.0);
298  X.insert( 80.0);
299  X.insert( 90.0);
300  X.insert(100.0);
301  X.insert(120.0);
302  X.insert(140.0);
303  X.insert(160.0);
304  X.insert(180.0);
305  X.insert(200.0);
306  X.insert(250.0);
307  X.insert(300.0);
308  X.insert(350.0);
309  X.insert(400.0);
310  X.insert(450.0);
311  X.insert(500.0);
312  X.insert(600.0);
313  X.insert(700.0);
314  X.insert(800.0);
315  X.insert(900.0);
316  X.insert(1200.0);
317  X.insert(1500.0);
318  }
319 
320 
321  const double grid = 5.0; // [deg]
322 
323  const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0)); // azimuth angle unit step size
324 
325 
326  for (set<double>::const_iterator r = R.begin(); r != R.end(); ++r) {
327 
328  const double R_m = *r;
329 
330  const unsigned int number_of_theta_points = max(2u, (unsigned int) (180.0/(1.4 * grid)));
331 
332  for (double theta = 0.0; theta <= PI + epsilon; theta += PI/number_of_theta_points) {
333 
334  const unsigned int number_of_phi_points = max(2u, (unsigned int) (PI * sin(theta) / alpha));
335 
336  for (double phi = 0.0; phi <= PI + epsilon; phi += PI/number_of_phi_points) {
337 
338  JFunction1D_t& f1 = pdf[R_m][theta][phi];
339 
340  const JArray_t array(R_m, theta, phi);
341 
342  double t_old = transformer.getXn(array, *X.begin());
343  double y_old = 0.0;
344 
345  for (set<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
346 
347  const double t = transformer.getXn(array, *x);
348  const double y = (pdf_c.*f)(R_m, theta, phi, t);
349 
350  if (y != 0.0) {
351 
352  if (*x < 0.0) {
353  WARNING("dt < 0 " << *x << ' ' << R_m << ' ' << t << ' ' << y << endl);
354  }
355 
356  if (y_old == 0.0) {
357  f1[t_old] = y_old;
358  }
359 
360  f1[t] = y;
361 
362  } else {
363 
364  if (y_old != 0.0) {
365  f1[t] = y;
366  }
367  }
368 
369  t_old = t;
370  y_old = y;
371  }
372  }
373  }
374  }
375 
376  pdf.transform(transformer);
377  pdf.compile();
378 
379  NOTICE("OK" << endl);
380 
381  try {
382 
383  NOTICE("storing output to file " << outputFile << "... " << flush);
384 
385  pdf.store(outputFile.c_str());
386 
387  NOTICE("OK" << endl);
388  }
389  catch(const JException& error) {
390  FATAL(error.what() << endl);
391  }
392 }
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
getAngularAcceptance
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:68
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
debug
int debug
debug level
Definition: JSirene.cc:69
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
numberOfPoints
int numberOfPoints
Definition: JResultPDF.cc:22
JLANG::JException
General exception.
Definition: JException.hh:24
JLANG::JException::what
virtual const char * what() const override
Get error message.
Definition: JException.hh:64
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1714
JPHYSICS::JPDFTable
Multi-dimensional PDF table for arrival time of Cherenkov light.
Definition: JPDFTable.hh:44
JPHYSICS::JPDFTransformer
Template definition of transformer of the probability density function (PDF) of the time response of ...
Definition: JPDFTransformer.hh:620
JPHYSICS::JPDF_C
Probability Density Functions of the time response of a PMT with an implementation of the JAbstractPM...
Definition: JPDF.hh:2185
JPHYSICS::JPDF
Low level interface for the calculation of the Probability Density Functions (PDFs) of the arrival ti...
Definition: JPDF.hh:282
JTOOLS::JArray
One dimensional array of template objects with fixed length.
Definition: JArray.hh:43
ANTARES::getAmbientPressure
double getAmbientPressure()
Get ambient pressure.
Definition: Antares.hh:40
ANTARES::getScatteringProbability
double getScatteringProbability(const double x)
Function to describe light scattering in water.
Definition: Antares.hh:210
ANTARES::getPhotocathodeArea
double getPhotocathodeArea()
Get photo-cathode area of PMT.
Definition: Antares.hh:51
GAUSS_LEGENDRE::f1
const JPolynome f1(1.0, 2.0, 3.0)
Function.
GAUSS_LEGENDRE::epsilon
const double epsilon
Definition: JQuadrature.cc:21
JDETECTOR::getQE
double getQE(const double R, const double mu)
Get QE for given ratio of hit probabilities and expectation value of the number of photo-electrons.
Definition: JPMTParametersToolkit.hh:89
JMATH::PI
static const double PI
Mathematical constants.
Definition: JMath/JConstants.hh:20
JPHYSICS::getMinimalWavelength
double getMinimalWavelength()
Get minimal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:35
JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS
@ SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition: JPDFTypes.hh:33
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWERS
@ DIRECT_LIGHT_FROM_EMSHOWERS
direct light from EM showers
Definition: JPDFTypes.hh:29
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWERS
@ SCATTERED_LIGHT_FROM_EMSHOWERS
scattered light from EM showers
Definition: JPDFTypes.hh:30
JPHYSICS::SCATTERED_LIGHT_FROM_MUON
@ SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition: JPDFTypes.hh:27
JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS
@ DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition: JPDFTypes.hh:32
JPHYSICS::DIRECT_LIGHT_FROM_MUON
@ DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition: JPDFTypes.hh:26
JPHYSICS::getMaximalWavelength
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:46
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
JTOOLS::r
data_type r[M+1]
Definition: JPolint.hh:868
std
Definition: JSTDTypes.hh:14
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:84
JPARSER::not_initialised
Empty structure for specification of parser element that is not initialised (i.e. does require input)...
Definition: JParser.hh:90
JTOOLS::JMAPLIST
Auxiliary class for recursive map list generation.
Definition: JMapList.hh:109
JTOOLS::JPolint1FunctionalGridMap
Type definition of a 1st degree polynomial interpolation based on a JGridMap implementation.
Definition: JFunctionalMap_t.hh:93
JTOOLS::JPolint1FunctionalMap
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
Definition: JFunctionalMap_t.hh:57
JTOOLS::JSplineFunction1D_t
Type definition of a spline interpolation method based on a JCollection with double result type.
Definition: JFunction1D_t.hh:37

Variable Documentation

◆ absorptionLengthFactor

double absorptionLengthFactor

Scaling of absorption and scattering length.

Definition at line 22 of file JMakePDF.cc.

◆ scatteringLengthFactor

double scatteringLengthFactor

Definition at line 23 of file JMakePDF.cc.

Generated by doxygen 1.9.1