Jpp 20.0.0-27-g39925593c-D
the software that should make you happy
Loading...
Searching...
No Matches
Functions | Variables
JMakePDG.cc File Reference

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

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <set>
#include <map>
#include <cmath>
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JQuadrature.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 shower.

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

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

Author
mdejong

Definition in file JMakePDG.cc.

Function Documentation

◆ getAbsorptionLength()

double getAbsorptionLength ( const double lambda)
inline

Definition at line 28 of file JMakePDG.cc.

29{
30 return absorptionLengthFactor * NAMESPACE::getAbsorptionLength(lambda);
31}
absorptionLengthFactor
double absorptionLengthFactor
Scaling of absorption and scattering length.
Definition JMakePDG.cc:24

◆ getScatteringLength()

double getScatteringLength ( const double lambda)
inline

Definition at line 34 of file JMakePDG.cc.

35{
36 return scatteringLengthFactor * NAMESPACE::getScatteringLength(lambda);
37}
scatteringLengthFactor
double scatteringLengthFactor
Definition JMakePDG.cc:25

◆ main()

int main ( int argc,
char ** argv )

Definition at line 55 of file JMakePDG.cc.

56{
57 using namespace std;
58 using namespace JPP;
59
60 string outputFile;
61 int numberOfPoints;
62 double epsilon;
63 int function;
64 set<double> D; // distance [m]
65 int debug;
66
67 try {
68
69 JParser<> zap("Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a shower.");
70
71 zap['o'] = make_field(outputFile);
72 zap['n'] = make_field(numberOfPoints, "points for integration") = 25;
73 zap['e'] = make_field(epsilon, "precision for integration") = 1.0e-10;
74 zap['A'] = make_field(absorptionLengthFactor, "scaling factor") = 1.0;
75 zap['S'] = make_field(scatteringLengthFactor, "scaling factor") = 1.0;
76 zap['F'] = make_field(function, "PDF type") =
77 SCATTERED_LIGHT_FROM_MUON_5D,
78 DIRECT_LIGHT_FROM_EMSHOWER,
79 SCATTERED_LIGHT_FROM_EMSHOWER;
80 zap['D'] = make_field(D, "distance [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,
96 const double) const;
97
98
99 // set global parameters
100
101 const double P_atm = NAMESPACE::getAmbientPressure();
102 const double wmin = getMinimalWavelength();
103 const double wmax = getMaximalWavelength();
104
105
106 const JPDF_C
107 pdf_c(NAMESPACE::getPhotocathodeArea(),
108 NAMESPACE::getQE,
109 NAMESPACE::getAngularAcceptance,
110 getAbsorptionLength,
111 getScatteringLength,
112 NAMESPACE::getScatteringProbability,
113 P_atm,
114 wmin,
115 wmax,
116 numberOfPoints,
117 epsilon);
118
119
120 typedef JSplineFunction1D_t JFunction1D_t;
121 typedef JMAPLIST<JPolint1FunctionalMap,
122 JPolint1FunctionalMap,
123 JPolint1FunctionalGridMap,
124 JPolint1FunctionalGridMap>::maplist JMapList_t;
125 typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
126
127 typedef JPDFTransformer<4, JFunction1D_t::argument_type> JFunction4DTransformer_t;
128 typedef JArray<4, JFunction1D_t::argument_type> JArray_t;
129
130 JPDF_t pdf;
131
132
133 NOTICE("building multi-dimensional function object <" << function << ">... " << flush);
134
135 const double ng[] = { pdf_c.getIndexOfRefractionGroup(wmax),
136 pdf_c.getIndexOfRefractionGroup(wmin) };
137
138 map<int, pair<fcn,JFunction4DTransformer_t> > zmap;
139
140 zmap[SCATTERED_LIGHT_FROM_MUON_5D] = make_pair((fcn) &JPDF::getScatteredLightFromMuon, JFunction4DTransformer_t(21.5, 2, ng[0], 0.0, JGeant(JGeanx(0.33, -9.5)), 6e-4, NAMESPACE::getAngularAcceptance, 0.06));
141 zmap[DIRECT_LIGHT_FROM_EMSHOWER] = make_pair((fcn) &JPDF::getDirectLightFromEMshower, JFunction4DTransformer_t(21.5, 2, ng[0], ng[1], JGeant(JGeanx(0.35, -5.4)), 1e-5, NAMESPACE::getAngularAcceptance, 0.001));
142 zmap[SCATTERED_LIGHT_FROM_EMSHOWER] = make_pair((fcn) &JPDF::getScatteredLightFromEMshower, JFunction4DTransformer_t(21.5, 2, ng[0], 0.0, JGeant(JGeanx(0.55, -4.5)), 1e-2, NAMESPACE::getAngularAcceptance, 0.05));
143
144 if (zmap.find(function) == zmap.end()) {
145 FATAL("illegal function specifier" << endl);
146 }
147
148 fcn f = zmap[function].first; // PDF
149 JFunction4DTransformer_t transformer = zmap[function].second; // transformer
150
151
152 if (D.empty()) {
153 D.insert( 0.10);
154 D.insert( 0.50);
155 D.insert( 1.00);
156 D.insert( 5.00);
157 D.insert( 10.00);
158 D.insert( 20.00);
159 D.insert( 30.00);
160 D.insert( 40.00);
161 D.insert( 50.00);
162 D.insert( 60.00);
163 D.insert( 70.00);
164 D.insert( 80.00);
165 D.insert( 90.00);
166 D.insert(100.00);
167 D.insert(120.00);
168 D.insert(150.00);
169 D.insert(170.00);
170 D.insert(190.00);
171 D.insert(210.00);
172 D.insert(230.00);
173 D.insert(250.00);
174 D.insert(270.00);
175 D.insert(290.00);
176 D.insert(310.00);
177 }
178
179 set<double> C; // cosine emission angle
180
181 JQuadrature qeant(-1.0, +1.0, 60, geanx);
182
183 for (JQuadrature::const_iterator i = qeant.begin(); i != qeant.end(); ++i)
184 C.insert(i->getX());
185
186 C.insert(-1.00);
187 C.insert(+1.00);
188
189
190 set<double> X; // time [ns]
191
192 if (function == DIRECT_LIGHT_FROM_EMSHOWER) {
193
194 for (double buffer[] = { 0.0, 0.005, 0.01, 0.015, -1 }, *x = buffer; *x >= 0; ++x) {
195 X.insert(0.0 + *x);
196 X.insert(1.0 - *x);
197 }
198
199 for (double x = 0.02; x < 0.99; x += 0.01)
200 X.insert(x);
201
202 } else {
203
204 X.insert( 0.00);
205 X.insert( 0.10);
206 X.insert( 0.20);
207 X.insert( 0.30);
208 X.insert( 0.40);
209 X.insert( 0.50);
210 X.insert( 0.60);
211 X.insert( 0.70);
212 X.insert( 0.80);
213 X.insert( 0.90);
214 X.insert( 1.00);
215 X.insert( 1.00);
216 X.insert( 1.10);
217 X.insert( 1.20);
218 X.insert( 1.30);
219 X.insert( 1.40);
220 X.insert( 1.50);
221 X.insert( 1.60);
222 X.insert( 1.70);
223 X.insert( 1.80);
224 X.insert( 1.90);
225 X.insert( 2.00);
226 X.insert( 2.20);
227 X.insert( 2.40);
228 X.insert( 2.60);
229 X.insert( 2.80);
230 X.insert( 3.00);
231 X.insert( 3.25);
232 X.insert( 3.50);
233 X.insert( 3.75);
234 X.insert( 4.00);
235 X.insert( 4.25);
236 X.insert( 4.50);
237 X.insert( 4.75);
238 X.insert( 5.0);
239 X.insert( 6.0);
240 X.insert( 7.0);
241 X.insert( 8.0);
242 X.insert( 9.0);
243 X.insert( 10.0);
244 X.insert( 15.0);
245 X.insert( 20.0);
246 X.insert( 25.0);
247 X.insert( 30.0);
248 X.insert( 40.0);
249 X.insert( 50.0);
250 X.insert( 60.0);
251 X.insert( 70.0);
252 X.insert( 80.0);
253 X.insert( 90.0);
254 X.insert(100.0);
255 X.insert(120.0);
256 X.insert(140.0);
257 X.insert(160.0);
258 X.insert(180.0);
259 X.insert(200.0);
260 X.insert(250.0);
261 X.insert(300.0);
262 X.insert(350.0);
263 X.insert(400.0);
264 X.insert(450.0);
265 X.insert(500.0);
266 X.insert(600.0);
267 X.insert(700.0);
268 X.insert(800.0);
269 X.insert(900.0);
270 X.insert(1200.0);
271 X.insert(1500.0);
272 }
273
274 const double grid = 7.0; // [deg]
275
276 const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0)); // azimuth angle unit step size
277
278
279 for (set<double>::const_iterator d = D.begin(); d != D.end(); ++d) {
280
281 const double D_m = *d;
282
283 for (set<double>::const_iterator c = C.begin(); c != C.end(); ++c) {
284
285 const double cd = *c;
286
287 const unsigned int number_of_theta_points = max(2u, (unsigned int) (180.0/(1.4 * grid)));
288
289 for (double theta = 0.0; theta <= PI + epsilon; theta += PI/number_of_theta_points) {
290
291 const unsigned int number_of_phi_points = max(2u, (unsigned int) (PI * sin(theta) / alpha));
292
293 for (double phi = 0.0; phi <= PI + epsilon; phi += PI/number_of_phi_points) {
294
295 JFunction1D_t& f1 = pdf[D_m][cd][theta][phi];
296
297 const JArray_t array(D_m, cd, theta, phi);
298
299 double t_old = transformer.getXn(array, *X.begin());
300 double y_old = 0.0;
301
302 for (set<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
303
304 const double t = transformer.getXn(array, *x);
305 const double y = (pdf_c.*f)(D_m, cd, theta, phi, t);
306
307 if (y != 0.0) {
308
309 if (*x < 0.0) {
310 WARNING("dt < 0 " << *x << ' ' << D_m << ' ' << t << ' ' << y << endl);
311 }
312
313 if (y_old == 0.0) {
314 f1[t_old] = y_old;
315 }
316
317 f1[t] = y;
318
319 } else {
320
321 if (y_old != 0.0) {
322 f1[t] = y;
323 }
324 }
325
326 t_old = t;
327 y_old = y;
328 }
329 }
330 }
331 }
332 }
333
334 pdf.transform(transformer);
335 pdf.compile();
336
337 NOTICE("OK" << endl);
338
339 try {
340
341 NOTICE("storing output to file " << outputFile << "... " << flush);
342
343 pdf.store(outputFile.c_str());
344
345 NOTICE("OK" << endl);
346 }
347 catch(const JException& error) {
348 FATAL(error.what() << endl);
349 }
350}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
getAbsorptionLength
double getAbsorptionLength(const double lambda)
Definition JMakePDG.cc:28
getScatteringLength
double getScatteringLength(const double lambda)
Definition JMakePDG.cc:34
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::JGeant
Function object for the probability density function of photon emission from EM-shower as a function ...
Definition JGeant.hh:32
JPHYSICS::JGeanx
Probability density function of photon emission from EM-shower as a function of cosine of the emissio...
Definition JGeanx.hh:32
JPHYSICS::JPDF_C
Probability Density Functions of the time response of a PMT with an implementation of the JAbstractPM...
Definition JPDF.hh:2186
JTOOLS::JCollection< JElement2D_t >::const_iterator
container_type::const_iterator const_iterator
Definition JCollection.hh:91
JTOOLS::JQuadrature
Type definition for numerical integration.
Definition JQuadrature.hh:39
JPHYSICS::getMinimalWavelength
double getMinimalWavelength()
Get minimal wavelength for PDF evaluations.
Definition JPDFToolkit.hh:27
JPHYSICS::SCATTERED_LIGHT_FROM_MUON_5D
@ SCATTERED_LIGHT_FROM_MUON_5D
scattered light from muon
Definition JPDFTypes.hh:35
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 24 of file JMakePDG.cc.

◆ scatteringLengthFactor

double scatteringLengthFactor

Definition at line 25 of file JMakePDG.cc.

Generated by doxygen 1.10.0