Jpp
JMakePDF.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <fstream>
4 #include <iomanip>
5 #include <set>
6 #include <map>
7 
8 #include "JTools/JFunction1D_t.hh"
9 #include "JTools/JFunctionalMap_t.hh"
10 #include "JPhysics/JPDF.hh"
11 #include "JPhysics/JPDFTable.hh"
12 #include "JPhysics/Antares.hh"
13 #include "JPhysics/KM3NeT.hh"
14 
15 #include "Jeep/JParser.hh"
16 #include "Jeep/JMessage.hh"
17 
18 
19 /**
20  * Scaling of absorption and scattering length.
21  */
22 double absorptionLengthFactor;
23 double scatteringLengthFactor;
24 
25 
26 inline double getAbsorptionLength(const double lambda)
27 {
28  return absorptionLengthFactor * NAMESPACE::getAbsorptionLength(lambda);
29 }
30 
31 
32 inline double getScatteringLength(const double lambda)
33 {
34  return scatteringLengthFactor * NAMESPACE::getScatteringLength(lambda);
35 }
36 
37 
38 /**
39  * \file
40  *
41  * Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a muon.
42  *
43  * The PDFs are tabulated as a function of <tt>(R, theta, phi, t)</tt>, where:
44  * - <tt>R</tt> is the minimal distance of approach of the muon to the PMT;
45  * - <tt>(theta, phi)</tt> the orientation of the PMT; and
46  * - <tt>t</tt> the arrival time of the light with respect to the Cherenkov hypothesis.
47  *
48  * The orientation of the PMT is defined in the coordinate system in which
49  * the muon travels along the z-axis and the PMT is located in the x-z plane.
50  * \author mdejong
51  */
52 int main(int argc, char **argv)
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 
93  typedef double (JPDF::*fcn)(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  JPolint1FunctionalGridMap,
123  JPolint1FunctionalGridMap>::maplist JMapList_t;
124  typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
125 
126  typedef JPDFTransformer<3, JFunction1D_t::argument_type> JFunction3DTransformer_t;
127  typedef JArray<3, JFunction1D_t::argument_type> JArray_t;
128 
129  JPDF_t pdf;
130 
131 
132  NOTICE("building multi-dimensional function object <" << function << ">... " << flush);
133 
134 
135  const double kmin = pdf_c.getKappa(wmax);
136  const double kmax = pdf_c.getKappa(wmin);
137  const double cmin = pdf_c.getKmin (wmax);
138 
139  map<int, pair<fcn,JFunction3DTransformer_t> > zmap;
140 
141  zmap[DIRECT_LIGHT_FROM_MUON] = make_pair((fcn) &JPDF::getDirectLightFromMuon, JFunction3DTransformer_t(21.5, 2, kmin, kmax, NAMESPACE::getAngularAcceptance, 0.001));
142  zmap[SCATTERED_LIGHT_FROM_MUON] = make_pair((fcn) &JPDF::getScatteredLightFromMuon, JFunction3DTransformer_t(35.0, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
143  zmap[DIRECT_LIGHT_FROM_EMSHOWERS] = make_pair((fcn) &JPDF::getDirectLightFromEMshowers, JFunction3DTransformer_t(21.5, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
144  zmap[SCATTERED_LIGHT_FROM_EMSHOWERS] = make_pair((fcn) &JPDF::getScatteredLightFromEMshowers, JFunction3DTransformer_t(35.0, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
145  zmap[DIRECT_LIGHT_FROM_DELTARAYS] = make_pair((fcn) &JPDF::getDirectLightFromDeltaRays, JFunction3DTransformer_t(21.5, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
146  zmap[SCATTERED_LIGHT_FROM_DELTARAYS] = make_pair((fcn) &JPDF::getScatteredLightFromDeltaRays, JFunction3DTransformer_t(35.0, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
147 
148  if (zmap.find(function) == zmap.end()) {
149  FATAL("illegal function specifier" << endl);
150  }
151 
152  fcn f = zmap[function].first; // PDF
153  JFunction3DTransformer_t transformer = zmap[function].second; // transformer
154 
155 
156  if (R.empty()) {
157  R.insert( 0.10);
158  R.insert( 0.30);
159  R.insert( 0.50);
160  R.insert( 1.00);
161  R.insert( 2.00);
162  R.insert( 3.00);
163  R.insert( 4.00);
164  R.insert( 5.00);
165  R.insert( 6.00);
166  R.insert( 7.00);
167  R.insert( 8.00);
168  R.insert( 9.00);
169  R.insert( 10.00);
170  R.insert( 11.00);
171  R.insert( 12.00);
172  R.insert( 13.00);
173  R.insert( 14.00);
174  R.insert( 15.00);
175  R.insert( 16.00);
176  R.insert( 17.00);
177  R.insert( 18.00);
178  R.insert( 19.00);
179  R.insert( 20.00);
180  R.insert( 22.00);
181  R.insert( 24.00);
182  R.insert( 26.00);
183  R.insert( 28.00);
184  R.insert( 30.00);
185  R.insert( 32.00);
186  R.insert( 34.00);
187  R.insert( 36.00);
188  R.insert( 38.00);
189  R.insert( 40.00);
190  R.insert( 42.00);
191  R.insert( 44.00);
192  R.insert( 46.00);
193  R.insert( 48.00);
194  R.insert( 50.00);
195  R.insert( 55.00);
196  R.insert( 60.00);
197  R.insert( 65.00);
198  R.insert( 70.00);
199  R.insert( 75.00);
200  R.insert( 80.00);
201  R.insert( 85.00);
202  R.insert( 90.00);
203  R.insert( 95.00);
204  R.insert(100.00);
205  R.insert(110.00);
206  R.insert(120.00);
207  R.insert(130.00);
208  R.insert(140.00);
209  R.insert(150.00);
210  R.insert(170.00);
211  R.insert(190.00);
212  R.insert(210.00);
213  R.insert(250.00);
214  }
215 
216  set<double> X;
217 
218  if (function == DIRECT_LIGHT_FROM_MUON) {
219 
220  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) {
221  X.insert(0.0 + *x);
222  X.insert(1.0 - *x);
223  }
224 
225  for (double x = 0.01; x < 0.1; x += 0.0025) {
226  X.insert(0.0 + x);
227  X.insert(1.0 - x);
228  }
229 
230  for (double x = 0.10; x < 0.5; x += 0.010) {
231  X.insert(0.0 + x);
232  X.insert(1.0 - x);
233  }
234 
235  } else {
236 
237  X.insert( 0.00);
238  X.insert( 0.01);
239  X.insert( 0.02);
240  X.insert( 0.03);
241  X.insert( 0.04);
242  X.insert( 0.05);
243  X.insert( 0.06);
244  X.insert( 0.07);
245  X.insert( 0.08);
246  X.insert( 0.09);
247  X.insert( 0.10);
248  X.insert( 0.12);
249  X.insert( 0.15);
250  X.insert( 0.20);
251  X.insert( 0.25);
252  X.insert( 0.30);
253  X.insert( 0.40);
254  X.insert( 0.50);
255  X.insert( 0.60);
256  X.insert( 0.70);
257  X.insert( 0.80);
258  X.insert( 0.90);
259  X.insert( 1.00);
260  X.insert( 1.10);
261  X.insert( 1.20);
262  X.insert( 1.30);
263  X.insert( 1.40);
264  X.insert( 1.50);
265  X.insert( 1.60);
266  X.insert( 1.70);
267  X.insert( 1.80);
268  X.insert( 1.90);
269  X.insert( 2.00);
270  X.insert( 2.20);
271  X.insert( 2.40);
272  X.insert( 2.60);
273  X.insert( 2.80);
274  X.insert( 3.00);
275  X.insert( 3.25);
276  X.insert( 3.50);
277  X.insert( 3.75);
278  X.insert( 4.00);
279  X.insert( 4.25);
280  X.insert( 4.50);
281  X.insert( 4.75);
282  X.insert( 5.0);
283  X.insert( 6.0);
284  X.insert( 7.0);
285  X.insert( 8.0);
286  X.insert( 9.0);
287  X.insert( 10.0);
288  X.insert( 12.0);
289  X.insert( 14.0);
290  X.insert( 16.0);
291  X.insert( 18.0);
292  X.insert( 20.0);
293  X.insert( 25.0);
294  X.insert( 30.0);
295  X.insert( 40.0);
296  X.insert( 50.0);
297  X.insert( 60.0);
298  X.insert( 70.0);
299  X.insert( 80.0);
300  X.insert( 90.0);
301  X.insert(100.0);
302  X.insert(120.0);
303  X.insert(140.0);
304  X.insert(160.0);
305  X.insert(180.0);
306  X.insert(200.0);
307  X.insert(250.0);
308  X.insert(300.0);
309  X.insert(350.0);
310  X.insert(400.0);
311  X.insert(450.0);
312  X.insert(500.0);
313  X.insert(600.0);
314  X.insert(700.0);
315  X.insert(800.0);
316  X.insert(900.0);
317  X.insert(1200.0);
318  X.insert(1500.0);
319  }
320 
321 
322  const double grid = 5.0; // [deg]
323 
324  const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0)); // azimuth angle unit step size
325 
326 
327  for (set<double>::const_iterator r = R.begin(); r != R.end(); ++r) {
328 
329  const double R_m = *r;
330 
331  const unsigned int number_of_theta_points = max(2u, (unsigned int) (180.0/(1.4 * grid)));
332 
333  for (double theta = 0.0; theta <= PI + epsilon; theta += PI/number_of_theta_points) {
334 
335  const unsigned int number_of_phi_points = max(2u, (unsigned int) (PI * sin(theta) / alpha));
336 
337  for (double phi = 0.0; phi <= PI + epsilon; phi += PI/number_of_phi_points) {
338 
339  JFunction1D_t& f1 = pdf[R_m][theta][phi];
340 
341  const JArray_t array(R_m, theta, phi);
342 
343  double t_old = transformer.getXn(array, *X.begin());
344  double y_old = 0.0;
345 
346  for (set<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
347 
348  const double t = transformer.getXn(array, *x);
349  const double y = (pdf_c.*f)(R_m, theta, phi, t);
350 
351  if (y != 0.0) {
352 
353  if (*x < 0.0) {
354  WARNING("dt < 0 " << *x << ' ' << R_m << ' ' << t << ' ' << y << endl);
355  }
356 
357  if (y_old == 0.0) {
358  f1[t_old] = y_old;
359  }
360 
361  f1[t] = y;
362 
363  } else {
364 
365  if (y_old != 0.0) {
366  f1[t] = y;
367  }
368  }
369 
370  t_old = t;
371  y_old = y;
372  }
373  }
374  }
375  }
376 
377  pdf.transform(transformer);
378  pdf.compile();
379 
380  NOTICE("OK" << endl);
381 
382  try {
383 
384  NOTICE("storing output to file " << outputFile << "... " << flush);
385 
386  pdf.store(outputFile.c_str());
387 
388  NOTICE("OK" << endl);
389  }
390  catch(const JException& error) {
391  FATAL(error.what() << endl);
392  }
393 }
JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition: JPDFTypes.hh:35
Antares.hh
JMessage.hh
JPDF.hh
getAbsorptionLength
double getAbsorptionLength(const double lambda)
Definition: JMakePDF.cc:26
numberOfPoints
int numberOfPoints
Definition: JResultPDF.cc:22
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e.
Definition: JParser.hh:63
ANTARES::getScatteringProbability
double getScatteringProbability(const double x)
Function to describe light scattering in water.
Definition: Antares.hh:292
JTOOLS::u
double u[N+1]
Definition: JPolint.hh:706
getScatteringLength
double getScatteringLength(const double lambda)
Definition: JMakePDF.cc:32
scatteringLengthFactor
double scatteringLengthFactor
Definition: JMakePDF.cc:23
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
main
int main(int argc, char **argv)
Definition: JMakePDF.cc:52
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
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
std::set
Definition: JSTDTypes.hh:13
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JFunction1D_t.hh
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
debug
int debug
debug level
Definition: JSirene.cc:59
ANTARES::getAmbientPressure
double getAmbientPressure()
Get ambient pressure.
Definition: Antares.hh:30
JPHYSICS::getMinimalWavelength
double getMinimalWavelength()
Get minimal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:38
JFunctionalMap_t.hh
std::map
Definition: JSTDTypes.hh:16
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWERS
scattered light from EM showers
Definition: JPDFTypes.hh:33
JParser.hh
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
absorptionLengthFactor
double absorptionLengthFactor
Scaling of absorption and scattering length.
Definition: JMakePDF.cc:22
std
Definition: jaanetDictionary.h:36
JTOOLS::PI
static const double PI
Constants.
Definition: JConstants.hh:20
JTOOLS::r
data_type r[M+1]
Definition: JPolint.hh:709
JPHYSICS::DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition: JPDFTypes.hh:29
getAngularAcceptance
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:84
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWERS
direct light from EM showers
Definition: JPDFTypes.hh:32
JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition: JPDFTypes.hh:36
KM3NeT.hh
JPHYSICS::getMaximalWavelength
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:49
JPDFTable.hh
JPARSER::not_initialised
Empty structure for specification of parser element that is not initialised (i.e.
Definition: JParser.hh:69
ANTARES::getPhotocathodeArea
const double getPhotocathodeArea()
Photo-cathode area 10 inch PMT.
Definition: Antares.hh:41
JPHYSICS::SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition: JPDFTypes.hh:30
Generated by   doxygen 1.8.16