Jpp  18.3.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMakePD0.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 #include <cmath>
8 
11 #include "JTools/JQuadrature.hh"
12 #include "JPhysics/JPDF.hh"
13 #include "JPhysics/JPDFTable.hh"
14 #include "JPhysics/Antares.hh"
15 #include "JPhysics/KM3NeT.hh"
16 
17 #include "Jeep/JParser.hh"
18 #include "Jeep/JMessage.hh"
19 
20 
21 /**
22  * Scaling of absorption and scattering length.
23  */
26 
27 
28 inline double getAbsorptionLength(const double lambda)
29 {
31 }
32 
33 
34 inline double getScatteringLength(const double lambda)
35 {
37 }
38 
39 
40 /**
41  * \file
42  *
43  * Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a bright point.
44  *
45  * The PDFs are tabulated as a function of <tt>(D, cos(theta), t)</tt>, where:
46  * - <tt>D</tt> is the distance between the bright point and the PMT;
47  * - <tt>cos(theta)</tt> the cosine of the PMT angle;
48  * - <tt>t</tt> the arrival time of the light with respect to the Cherenkov hypothesis.
49  *
50  * The orientation of the PMT is defined in the coordinate system in which
51  * the position of the bright point and that of the PMT are along the x-axis and the PMT is oriented in the x-y plane.
52  * \author mdejong
53  */
54 int main(int argc, char **argv)
55 {
56  using namespace std;
57  using namespace JPP;
58 
59  string outputFile;
60  int numberOfPoints;
61  double epsilon;
62  int function;
63  set<double> D; // distance [m]
64  int debug;
65 
66  try {
67 
68  JParser<> zap("Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a bright point.");
69 
70  zap['o'] = make_field(outputFile);
71  zap['n'] = make_field(numberOfPoints, "points for integration") = 25;
72  zap['e'] = make_field(epsilon, "precision for integration") = 1.0e-10;
73  zap['A'] = make_field(absorptionLengthFactor, "scaling factor") = 1.0;
74  zap['S'] = make_field(scatteringLengthFactor, "scaling factor") = 1.0;
75  zap['F'] = make_field(function, "PDF type") =
78  zap['D'] = make_field(D, "distance [m]") = JPARSER::initialised();
79  zap['d'] = make_field(debug) = 0;
80 
81  zap['F'] = JPARSER::not_initialised();
82 
83  zap(argc, argv);
84  }
85  catch(const exception &error) {
86  FATAL(error.what() << endl);
87  }
88 
89 
90  typedef double (JPDF::*fcn)(const double,
91  const double,
92  const double) const;
93 
94 
95  // set global parameters
96 
97  const double P_atm = NAMESPACE::getAmbientPressure();
98  const double wmin = getMinimalWavelength();
99  const double wmax = getMaximalWavelength();
100 
101 
102  const JPDF_C
109  P_atm,
110  wmin,
111  wmax,
113  epsilon);
114 
115 
116  typedef JSplineFunction1D_t JFunction1D_t;
117  typedef JMAPLIST<JPolint1FunctionalMap,
118  JPolint1FunctionalGridMap>::maplist JMapList_t;
119  typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
120 
121  typedef JPDFTransformer<2, JFunction1D_t::argument_type> JFunction2DTransformer_t;
122  typedef JArray<2, JFunction1D_t::argument_type> JArray_t;
123 
124  JPDF_t pdf;
125 
126 
127  NOTICE("building multi-dimensional function object <" << function << ">... " << flush);
128 
129  const double ng[] = { pdf_c.getIndexOfRefractionGroup(wmax),
130  pdf_c.getIndexOfRefractionGroup(wmin) };
131 
133 
134  zmap[DIRECT_LIGHT_FROM_BRIGHT_POINT] = make_pair((fcn) &JPDF::getDirectLightFromBrightPoint, JFunction2DTransformer_t(21.5, 2, ng[0], ng[1]));
135  zmap[SCATTERED_LIGHT_FROM_BRIGHT_POINT] = make_pair((fcn) &JPDF::getScatteredLightFromBrightPoint, JFunction2DTransformer_t(21.5, 2, ng[0], 0.0));
136 
137  if (zmap.find(function) == zmap.end()) {
138  FATAL("illegal function specifier" << endl);
139  }
140 
141  fcn f = zmap[function].first; // PDF
142  JFunction2DTransformer_t transformer = zmap[function].second; // transformer
143 
144 
145  if (D.empty()) {
146  D.insert( 0.10);
147  D.insert( 0.50);
148  D.insert( 1.00);
149  D.insert( 5.00);
150  D.insert( 10.00);
151  D.insert( 20.00);
152  D.insert( 30.00);
153  D.insert( 40.00);
154  D.insert( 50.00);
155  D.insert( 60.00);
156  D.insert( 70.00);
157  D.insert( 80.00);
158  D.insert( 90.00);
159  D.insert(100.00);
160  }
161 
162  set<double> X; // time [ns]
163 
164  if (function == DIRECT_LIGHT_FROM_BRIGHT_POINT) {
165 
166  for (double buffer[] = { 0.0, 0.005, 0.01, 0.015, -1 }, *x = buffer; *x >= 0; ++x) {
167  X.insert(0.0 + *x);
168  X.insert(1.0 - *x);
169  }
170 
171  for (double x = 0.02; x < 0.99; x += 0.01)
172  X.insert(x);
173 
174  } else {
175 
176  X.insert( 0.00);
177  X.insert( 0.10);
178  X.insert( 0.20);
179  X.insert( 0.30);
180  X.insert( 0.40);
181  X.insert( 0.50);
182  X.insert( 0.60);
183  X.insert( 0.70);
184  X.insert( 0.80);
185  X.insert( 0.90);
186  X.insert( 1.00);
187  X.insert( 1.00);
188  X.insert( 1.10);
189  X.insert( 1.20);
190  X.insert( 1.30);
191  X.insert( 1.40);
192  X.insert( 1.50);
193  X.insert( 1.60);
194  X.insert( 1.70);
195  X.insert( 1.80);
196  X.insert( 1.90);
197  X.insert( 2.00);
198  X.insert( 2.20);
199  X.insert( 2.40);
200  X.insert( 2.60);
201  X.insert( 2.80);
202  X.insert( 3.00);
203  X.insert( 3.25);
204  X.insert( 3.50);
205  X.insert( 3.75);
206  X.insert( 4.00);
207  X.insert( 4.25);
208  X.insert( 4.50);
209  X.insert( 4.75);
210  X.insert( 5.0);
211  X.insert( 6.0);
212  X.insert( 7.0);
213  X.insert( 8.0);
214  X.insert( 9.0);
215  X.insert( 10.0);
216  X.insert( 15.0);
217  X.insert( 20.0);
218  X.insert( 25.0);
219  X.insert( 30.0);
220  X.insert( 40.0);
221  X.insert( 50.0);
222  X.insert( 60.0);
223  X.insert( 70.0);
224  X.insert( 80.0);
225  X.insert( 90.0);
226  X.insert(100.0);
227  X.insert(120.0);
228  X.insert(140.0);
229  X.insert(160.0);
230  X.insert(180.0);
231  X.insert(200.0);
232  }
233 
234 
235  for (set<double>::const_iterator d = D.begin(); d != D.end(); ++d) {
236 
237  const double D_m = *d;
238 
239  for (double dc = 0.1, ct = -1.0; ct < +1.0 + 0.5*dc; ct += dc) {
240 
241  JFunction1D_t& f1 = pdf[D_m][ct];
242 
243  const JArray_t array(D_m, ct);
244 
245  double t_old = transformer.getXn(array, *X.begin());
246  double y_old = 0.0;
247 
248  for (set<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
249 
250  const double t = transformer.getXn(array, *x);
251  const double y = (pdf_c.*f)(D_m, ct, t);
252 
253  if (y != 0.0) {
254 
255  if (*x < 0.0) {
256  WARNING("dt < 0 " << *x << ' ' << D_m << ' ' << t << ' ' << y << endl);
257  }
258 
259  if (y_old == 0.0) {
260  f1[t_old] = y_old;
261  }
262 
263  f1[t] = y;
264 
265  } else {
266 
267  if (y_old != 0.0) {
268  f1[t] = y;
269  }
270  }
271 
272  t_old = t;
273  y_old = y;
274  }
275  }
276  }
277 
278  pdf.transform(transformer);
279  pdf.compile();
280 
281  NOTICE("OK" << endl);
282 
283  try {
284 
285  NOTICE("storing output to file " << outputFile << "... " << flush);
286 
287  pdf.store(outputFile.c_str());
288 
289  NOTICE("OK" << endl);
290  }
291  catch(const JException& error) {
292  FATAL(error.what() << endl);
293  }
294 }
Utility class to parse command line options.
Definition: JParser.hh:1514
#define WARNING(A)
Definition: JMessage.hh:65
int main(int argc, char *argv[])
Definition: Main.cc:15
o $QUALITY_ROOT d $DEBUG!CHECK_EXIT_CODE JPlot1D f
Definition: JDataQuality.sh:76
Properties of Antares PMT and deep-sea water.
double absorptionLengthFactor
Scaling of absorption and scattering length.
Definition: JDrawPD0.cc:24
double getScatteringProbability(const double x)
Function to describe light scattering in water.
Definition: Antares.hh:210
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
double getMinimalWavelength()
Get minimal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:26
direct light from bright point
Definition: JPDFTypes.hh:42
string outputFile
Various implementations of functional maps.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
Properties of KM3NeT PMT and deep-sea water.
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:37
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
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...
double getAmbientPressure()
Get ambient pressure.
Definition: Antares.hh:40
#define NOTICE(A)
Definition: JMessage.hh:64
double getPhotocathodeArea()
Get photo-cathode area of PMT.
Definition: Antares.hh:51
General purpose messaging.
scattered light from bright point
Definition: JPDFTypes.hh:43
#define FATAL(A)
Definition: JMessage.hh:67
Empty structure for specification of parser element that is not initialised (i.e. does require input)...
Definition: JParser.hh:89
Auxiliary classes for numerical integration.
double getAbsorptionLength(const double lambda)
Get absorption length.
Definition: Antares.hh:63
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
Definition: JMuonPath.sh:47
Utility class to parse command line options.
double getScatteringLength(const double lambda)
Get scattering length.
Definition: Antares.hh:148
int numberOfPoints
Definition: JResultPDF.cc:22
double scatteringLengthFactor
Definition: JDrawPD0.cc:25
no fit printf nominal n $STRING awk v X
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
const double epsilon
Definition: JQuadrature.cc:21
int debug
debug level
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:68
Auxiliary data structure for muon PDF.
Definition: JPDF_t.hh:26