Jpp
Functions | Variables
JMakePDF.cc File Reference
#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

◆ getAbsorptionLength()

double getAbsorptionLength ( const double  lambda)
inline

Definition at line 26 of file JMakePDF.cc.

27 {
28  return absorptionLengthFactor * NAMESPACE::getAbsorptionLength(lambda);
29 }

◆ getScatteringLength()

double getScatteringLength ( const double  lambda)
inline

Definition at line 32 of file JMakePDF.cc.

33 {
34  return scatteringLengthFactor * NAMESPACE::getScatteringLength(lambda);
35 }

◆ 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 
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 }

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.

JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition: JPDFTypes.hh:35
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
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
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
std::map
Definition: JSTDTypes.hh:16
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWERS
scattered light from EM showers
Definition: JPDFTypes.hh:33
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
JPHYSICS::getMaximalWavelength
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:49
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