Jpp 20.0.0-27-g39925593c-D
the software that should make you happy
Loading...
Searching...
No Matches
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.
 
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}
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}
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
getAbsorptionLength
double getAbsorptionLength(const double lambda)
Definition JMakePDF.cc:26
getScatteringLength
double getScatteringLength(const double lambda)
Definition JMakePDF.cc:32
NOTICE
#define NOTICE(A)
Definition JMessage.hh:64
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
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
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
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JPHYSICS::JPDF_C
Probability Density Functions of the time response of a PMT with an implementation of the JAbstractPM...
Definition JPDF.hh:2186
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:45
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
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.10.0