Jpp  17.1.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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. More...
 
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

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: JDrawPD0.cc:24
ANTARES::getAbsorptionLength
double getAbsorptionLength(const double lambda)
Get absorption length.
Definition: Antares.hh:63
double getScatteringLength ( const double  lambda)
inline

Definition at line 32 of file JMakePDF.cc.

33 {
34  return scatteringLengthFactor * NAMESPACE::getScatteringLength(lambda);
35 }
ANTARES::getScatteringLength
double getScatteringLength(const double lambda)
Get scattering length.
Definition: Antares.hh:148
scatteringLengthFactor
double scatteringLengthFactor
Definition: JDrawPD0.cc:25
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 
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 }
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1517
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
f
o $QUALITY_ROOT d $DEBUG!JPlot1D f
Definition: JDataQuality.sh:66
absorptionLengthFactor
double absorptionLengthFactor
Scaling of absorption and scattering length.
Definition: JDrawPD0.cc:24
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWERS
direct light from EM showers
Definition: JPDFTypes.hh:29
ANTARES::getScatteringProbability
double getScatteringProbability(const double x)
Function to describe light scattering in water.
Definition: Antares.hh:210
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
JPHYSICS::getMinimalWavelength
double getMinimalWavelength()
Get minimal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:26
JTOOLS::r
data_type r[M+1]
Definition: JPolint.hh:779
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JPHYSICS::DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition: JPDFTypes.hh:26
GAUSS_LEGENDRE::f1
const JPolynome f1(1.0, 2.0, 3.0)
Function.
JPHYSICS::SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition: JPDFTypes.hh:27
JPHYSICS::getMaximalWavelength
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:37
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1993
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
ANTARES::getAmbientPressure
double getAmbientPressure()
Get ambient pressure.
Definition: Antares.hh:40
JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition: JPDFTypes.hh:33
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JMATH::PI
static const double PI
Mathematical constants.
Definition: JMath/JConstants.hh:20
ANTARES::getPhotocathodeArea
double getPhotocathodeArea()
Get photo-cathode area of PMT.
Definition: Antares.hh:51
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWERS
scattered light from EM showers
Definition: JPDFTypes.hh:30
R
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:43
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition: JPDFTypes.hh:32
JPARSER::not_initialised
Empty structure for specification of parser element that is not initialised (i.e. does require input)...
Definition: JParser.hh:89
ANTARES::getAbsorptionLength
double getAbsorptionLength(const double lambda)
Get absorption length.
Definition: Antares.hh:63
ANTARES::getScatteringLength
double getScatteringLength(const double lambda)
Get scattering length.
Definition: Antares.hh:148
numberOfPoints
int numberOfPoints
Definition: JResultPDF.cc:22
scatteringLengthFactor
double scatteringLengthFactor
Definition: JDrawPD0.cc:25
X
no fit printf nominal n $STRING awk v X
Definition: detector-XY:fit.sh:102
JTOOLS::u
double u[N+1]
Definition: JPolint.hh:776
GAUSS_LEGENDRE::epsilon
const double epsilon
Definition: JQuadrature.cc:21
debug
int debug
debug level
Definition: archive-put-Detectors.sh:95
getAngularAcceptance
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:68
JPDF
Auxiliary data structure for muon PDF.
Definition: JPDF_t.hh:26

Variable Documentation

double absorptionLengthFactor

Scaling of absorption and scattering length.

Definition at line 22 of file JMakePDF.cc.

double scatteringLengthFactor

Definition at line 23 of file JMakePDF.cc.

Generated by   doxygen 1.8.5