Jpp 20.0.0-195-g190c9e876
the software that should make you happy
Loading...
Searching...
No Matches
Functions
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 "JPhysics/JPDFSupportkit.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

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:

  • R is the minimal distance of approach of the muon to the PMT;
  • (theta, phi) the orientation of the PMT; and
  • t the arrival time of the light with respect to the Cherenkov hypothesis.

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

◆ main()

int main ( int argc,
char ** argv )

Definition at line 36 of file JMakePDF.cc.

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