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