Jpp 20.0.0-rc.3
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JMakePDE.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 <cmath>
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JMultiMapTransformer.hh"
#include "JTools/JQuadrature.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JGeanx.hh"
#include "JPhysics/JGeanz.hh"
#include "JSystem/JSystemToolkit.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JConstants.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

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 (E, D, cos(theta), theta, phi, t), where:

  • E is the logarithm of the energy;
  • D is the distance between the vertex and the PMT;
  • cos(theta) the cosine of the photon emission angle;
  • (theta, phi) the orientation of the PMT; and
  • t the arrival time of the light with respect to the Cherenkov hypothesis.

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
jseneca, mdejong

Definition in file JMakePDE.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 42 of file JMakePDE.cc.

43{
44 using namespace std;
45 using namespace JPP;
46
47 string fileDescriptor;
48 string outputFile;
49 int numberOfPoints;
50 size_t elong_sample_count;
51 size_t cosine_angle_count; // Cosine angle is not as sharp when elongated
52 double epsilon;
53 double TTS;
54 set<double> Ds; // distance [m]
55 set<double> Es; // logarithm of energy (GeV)
56 bool option;
57 int debug;
58
59 try {
60
61 JParser<> zap("Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a shower.");
62
63 zap['F'] = make_field(fileDescriptor, "file name descriptor for PDF tables");
64 zap['o'] = make_field(outputFile, "file name output PDF");
65 zap['n'] = make_field(numberOfPoints, "points for integration") = 25;
66 zap['e'] = make_field(epsilon, "precision for integration") = 1.0e-10;
67 zap['D'] = make_field(Ds, "distance [m]") = JPARSER::initialised();
68 zap['E'] = make_field(Es, "energy [log10(GeV)]") = JPARSER::initialised();
69 zap['C'] = make_field(cosine_angle_count, "points for cosine angle") = 60;
70 zap['N'] = make_field(elong_sample_count, "points for elongation") = 20;
71 zap['T'] = make_field(TTS, "sigma of time bluring [ns]") = 2.0;
72 zap['O'] = make_field(option, "optional transformation");
73 zap['d'] = make_field(debug) = 0;
74
75 zap(argc, argv);
76 }
77 catch(const exception &error) {
78 FATAL(error.what() << endl);
79 }
80
81
82 if (cosine_angle_count == 0) {
83 FATAL("Invalid option -C <" << cosine_angle_count << ">" << endl);
84 }
85
86 typedef JPolint1Function1D_t JFunction1D_t;
87
88 typedef JMAPLIST<JPolint1FunctionalMap,
89 JPolint1FunctionalMap,
90 JPolint1FunctionalGridMap,
91 JPolint1FunctionalGridMap>::maplist J4DMapList_t;
92
93 typedef JMAPLIST<JPolint1FunctionalMap,
94 JPolint1FunctionalMap,
95 JPolint1FunctionalMap,
96 JPolint1FunctionalGridMap,
97 JPolint1FunctionalGridMap>::maplist J5DMapList_t;
98
99 typedef JPDFTable<JFunction1D_t, J4DMapList_t> JPDF4d_t;
100 typedef JPDFTable<JFunction1D_t, J5DMapList_t> JPDF5d_t;
101
102 typedef JPDFTransformer<4, JFunction1D_t::argument_type> JFunction4DTransformer_t;
103 typedef JPDFTransformer<5, JFunction1D_t::argument_type> JFunction5DTransformer_t;
104 typedef JArray<5, JFunction1D_t::argument_type> JArray_t;
105
106 JPDF4d_t i_pdf;
107
108 DEBUG("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
109 {
110 const string file_name = getFilename(fileDescriptor, SCATTERED_LIGHT_FROM_EMSHOWER);
111
112 NOTICE("loading input from file " << file_name << "... " << flush);
113
114 i_pdf.load(file_name.c_str());
115
116 NOTICE("OK" << endl);
117 }
118 DEBUG("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
119 {
120 JPDF4d_t pdf;
121
122 const string file_name = getFilename(fileDescriptor, DIRECT_LIGHT_FROM_EMSHOWER);
123
124 NOTICE("loading input from file " << file_name << "... " << flush);
125
126 pdf.load(file_name.c_str());
127
128 NOTICE("OK" << endl);
129 DEBUG("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
130
131 pdf.setExceptionHandler(new JPDF4d_t::JDefaultResult(JMATH::zero));
132
133 NOTICE("Add PDF... " << flush);
134
135 JTimer timer;
136
137 timer.start();
138
139 i_pdf.add(pdf);
140
141 timer.stop();
142
143 if (debug >= debug_t) {
144 timer.print(cout, unit_t);
145 }
146
147 NOTICE("OK" << endl);
148 }
149 DEBUG("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
150
151 i_pdf.setExceptionHandler(new JPDF4d_t::JDefaultResult(JMATH::zero));
152
153 if ( TTS > 0.0 ){
154
155 NOTICE("Smearing combined table with sigma [ns] = " << TTS << " ..." << flush);
156
157 i_pdf.blur( TTS );
158
159 NOTICE("OK" << endl);
160 }
161
162 JPDF5d_t o_pdf;
163
164 const JFunction4DTransformer_t* i_transformer = dynamic_cast<JFunction4DTransformer_t*>(i_pdf.transformer->clone());
165 const JFunction5DTransformer_t* o_transformer = (i_transformer != NULL ? new JFunction5DTransformer_t(*i_transformer) : NULL);
166
167 if (o_transformer != NULL && debug >= debug_t) {
168 o_transformer->print(cout);
169 }
170
171 if (Es.empty()) {
172 Es.insert( 1.0 );
173 Es.insert( 2.0 );
174 Es.insert( 3.0 );
175 Es.insert( 3.5 );
176 Es.insert( 4.0 );
177 Es.insert( 4.5 );
178 Es.insert( 5.0 );
179 Es.insert( 5.5 );
180 Es.insert( 6.0 );
181 Es.insert( 6.5 );
182 Es.insert( 7.0 );
183 Es.insert( 7.5 );
184 Es.insert( 8.0 );
185 };
186
187 if (Ds.empty()) {
188 Ds.insert( 0.10);
189 Ds.insert( 0.50);
190 Ds.insert( 1.00);
191 Ds.insert( 2.00);
192 Ds.insert( 5.00);
193 Ds.insert( 8.00);
194 Ds.insert( 10.00);
195 Ds.insert( 20.00);
196 Ds.insert( 30.00);
197 Ds.insert( 40.00);
198 Ds.insert( 50.00);
199 Ds.insert( 60.00);
200 Ds.insert( 70.00);
201 Ds.insert( 80.00);
202 Ds.insert( 90.00);
203 Ds.insert(100.00);
204 Ds.insert(120.00);
205 Ds.insert(150.00);
206 Ds.insert(170.00);
207 Ds.insert(190.00);
208 Ds.insert(210.00);
209 Ds.insert(230.00);
210 Ds.insert(250.00);
211 Ds.insert(270.00);
212 Ds.insert(290.00);
213 Ds.insert(310.00);
214 Ds.insert(400.00);
215 Ds.insert(800.00); // ARCA is big, we go to 800 m
216 }
217
218 set<double> C; // cosine emission angle
219
220 JQuadrature qeant(-1.0, +1.0, cosine_angle_count, geanx);
221
222 for (JQuadrature::const_iterator i = qeant.begin(); i != qeant.end(); ++i) {
223 C.insert(i->getX());
224 }
225
226 C.insert(-1.00);
227 C.insert(+1.00);
228
229 set<double> X; // time [ns]
230
231 X.insert( 0.00);
232 X.insert( 0.50);
233 X.insert( 1.00);
234 X.insert( 1.50);
235 X.insert( 2.00);
236 X.insert( 2.50);
237 X.insert( 3.00);
238 X.insert( 3.50);
239 X.insert( 4.00);
240 X.insert( 5.0);
241 X.insert( 6.0);
242 X.insert( 8.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( 80.0);
252 X.insert(100.0);
253 X.insert(120.0);
254 X.insert(150.0);
255 X.insert(200.0);
256 X.insert(250.0);
257 X.insert(300.0);
258 X.insert(350.0);
259 X.insert(400.0);
260 X.insert(450.0);
261 X.insert(500.0);
262 X.insert(600.0);
263 X.insert(700.0);
264 X.insert(800.0);
265 X.insert(900.0);
266 X.insert(1200.0);
267 X.insert(1500.0);
268
269 const double grid = 7.0; // [deg]
270 const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0)); // azimuth angle unit step size
271
272 NOTICE("building multi-dimensional function object" << endl);
273
274 for (const double E : Es) {
275
276 DEBUG("E: " << E << endl);
277
278 vector< double > elong_lengths;
279
280 for (size_t i = 0; i != elong_sample_count ; ++i ){
281
282 const double fraction = (double) (i + 0.5) / (double) elong_sample_count;
283 const double len = geanz.getLength(pow(10.0, E), fraction, 1e-6);
284
285 DEBUG("fraction: " << fraction << ", length [m]: "<< len << endl);
286
287 elong_lengths.push_back(len);
288 }
289
290 if (elong_lengths.empty()) {
291 elong_lengths.push_back(0.0);
292 }
293
294
295 for (const double D_m : Ds) {
296
297 DEBUG("D_m: " << D_m << endl);
298
299 for (const double cd : C) {
300
301 NOTICE("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
302
303 const unsigned int number_of_theta_points = max(2u, (unsigned int) (180.0/(1.4 * grid)));
304
305 for (double theta = 0.0; theta <= PI + epsilon; theta += PI/number_of_theta_points) {
306
307 const unsigned int number_of_phi_points = max(2u, (unsigned int) (PI * sin(theta) / alpha));
308
309 for (double phi = 0.0; phi <= PI + epsilon; phi += PI/number_of_phi_points) {
310
311 JFunction1D_t& f1 = o_pdf[E][D_m][cd][theta][phi];
312
313 f1.reserve(X.size());
314
315 const JArray_t array(E, D_m, cd, theta, phi);
316
317 double t_old = (o_transformer != NULL ? o_transformer->getXn(array, *X.begin()) : *X.begin());
318 double y_old = 0.0;
319
320 for (set<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
321
322 const double t = (o_transformer != NULL ? o_transformer->getXn(array, *x) : *x);
323
324 double y = 0.0;
325
326 for (const double Z : elong_lengths) {
327
328 const double Z0 = cd * D_m;
329 const double Z1 = Z0 - Z;
330
331 const double D_m_prime = sqrt(D_m * D_m + Z1 * Z1 - Z0 * Z0);
332 const double cd_prime = Z1 / D_m_prime;
333 const double t_prime = t - (Z) / getSpeedOfLight() - (D_m_prime - D_m) * getIndexOfRefraction() / getSpeedOfLight();
334
335 y += i_pdf(D_m_prime, cd_prime, theta, phi, t_prime);
336 }
337
338 y /= (double) elong_lengths.size();
339
340 if (y != 0.0) {
341
342 if (*x < 0.0) {
343 WARNING("dt < 0 " << *x << ' ' << D_m << ' ' << t << ' ' << y << endl);
344 }
345
346 if (y_old == 0.0) {
347 f1[t_old] = y_old;
348 }
349
350 f1[t] = y;
351
352 } else {
353
354 if (y_old != 0.0) {
355 f1[t] = y;
356 }
357 }
358
359 t_old = t;
360 y_old = y;
361 }
362 }
363 }
364 }
365 }
366 }
367 NOTICE("OK" << endl);
368 NOTICE("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
369
370 o_pdf.compile();
371
372 if (o_transformer != NULL && option) {
373
374 NOTICE("transform... " << flush);
375
376 o_pdf.transform(*o_transformer);
377
378 NOTICE("OK" << endl);
379 }
380
381 try {
382
383 NOTICE("storing output to file " << outputFile << "... " << flush);
384
385 o_pdf.store(outputFile.c_str());
386
387 NOTICE("OK" << endl);
388 }
389 catch(const JException& error) {
390 FATAL(error.what() << endl);
391 }
392}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
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
JEEP::JTimer
Auxiliary class for CPU timing and usage.
Definition JTimer.hh:33
JEEP::JTimer::print
void print(std::ostream &out, const JScale_t scale=milli_t) const
Print timer data.
Definition JTimer.hh:172
JEEP::JTimer::stop
void stop()
Stop timer.
Definition JTimer.hh:127
JEEP::JTimer::start
void start()
Start timer.
Definition JTimer.hh:106
JLANG::JException
General exception.
Definition JException.hh:24
JLANG::JException::what
virtual const char * what() const override
Get error message.
Definition JException.hh:64
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JPHYSICS::JPDFTable
Multi-dimensional PDF table for arrival time of Cherenkov light.
Definition JPDFTable.hh:44
JPHYSICS::JPDFTransformer
Template definition of transformer of the probability density function (PDF) of the time response of ...
Definition JPDFTransformer.hh:620
JTOOLS::JArray
One dimensional array of template objects with fixed length.
Definition JArray.hh:43
JTOOLS::JCollection::insert
pair_type insert(const value_type &element)
Insert element.
Definition JCollection.hh:316
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
JMATH::zero
static const JZero zero
Function object to assign zero value.
Definition JZero.hh:105
JMATH::pow
T pow(const T &x, const double y)
Power .
Definition JMath.hh:97
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
JSYSTEM::getMemoryUsage
float getMemoryUsage(JShell &shell, const pid_t pid)
Get memory usage in percent of given process identifier.
Definition JSystemToolkit.hh:51
JTOOLS::if
if((p==this->begin() &&this->getDistance(x,(p++) ->getX()) > distance_type::precision)||(p==this->end() &&this->getDistance((--p) ->getX(), x) > distance_type::precision))
Template base class for polynomial interpolations with polynomial result.
Definition JPolint.hh:775
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
JTOOLS::JMAPLIST
Auxiliary class for recursive map list generation.
Definition JMapList.hh:109
JTOOLS::JPolint1Function1D_t
Type definition of a 1st degree polynomial interpolation with result type double.
Definition JFunction1D_t.hh:187
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
Generated by doxygen 1.12.0