Jpp  19.1.0-rc.1
the software that should make you happy
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:

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:69
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
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:1714
JPHYSICS::JGeanz::getLength
double getLength(const double E, const double P, const double eps=1.0e-3) const
Get shower length for a given integrated probability.
Definition: JGeanz.hh:122
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::JQuadrature
Type definition for numerical integration.
Definition: JQuadrature.hh:39
GAUSS_LEGENDRE::f1
const JPolynome f1(1.0, 2.0, 3.0)
Function.
GAUSS_LEGENDRE::epsilon
const double epsilon
Definition: JQuadrature.cc:21
JEEP::getFilename
std::string getFilename(const std::string &file_name)
Get file name part, i.e. part after last JEEP::PATHNAME_SEPARATOR if any.
Definition: JeepToolkit.hh:128
JEEP::unit_t
@ unit_t
unit
Definition: JScale.hh:31
JEEP::debug_t
@ debug_t
debug
Definition: JMessage.hh:29
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
JMATH::PI
static const double PI
Mathematical constants.
Definition: JMath/JConstants.hh:20
JPHYSICS::geanx
static const JGeanx geanx(0.35, -5.40)
Function object for the number of photons from EM-shower as a function of emission angle.
JPHYSICS::geanz
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
JPHYSICS::getIndexOfRefraction
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
Definition: JPhysics/JConstants.hh:134
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER
@ SCATTERED_LIGHT_FROM_EMSHOWER
scattered light from EM shower
Definition: JPDFTypes.hh:38
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER
@ DIRECT_LIGHT_FROM_EMSHOWER
direct light from EM shower
Definition: JPDFTypes.hh:37
JPHYSICS::getSpeedOfLight
const double getSpeedOfLight()
Get speed of light.
Definition: JPhysics/JConstants.hh:112
JPHYSICS::C
static const double C
Physics constants.
Definition: JPhysics/JConstants.hh:25
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
JSYSTEM::getMemoryUsage
float getMemoryUsage(JShell &shell, const pid_t pid)
Get memory usage in percent of given process identifier.
Definition: JSystemToolkit.hh:51
std
Definition: JSTDTypes.hh:14
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:84
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.9.1