Jpp
Functions | Variables
JMakePDG.cc File Reference
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <set>
#include <map>
#include <cmath>
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JQuadrature.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 shower.

The PDFs are tabulated as a function of (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
mdejong

Definition in file JMakePDG.cc.

Function Documentation

◆ getAbsorptionLength()

double getAbsorptionLength ( const double  lambda)
inline

Definition at line 28 of file JMakePDG.cc.

29 {
30  return absorptionLengthFactor * NAMESPACE::getAbsorptionLength(lambda);
31 }

◆ getScatteringLength()

double getScatteringLength ( const double  lambda)
inline

Definition at line 34 of file JMakePDG.cc.

35 {
36  return scatteringLengthFactor * NAMESPACE::getScatteringLength(lambda);
37 }

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 55 of file JMakePDG.cc.

56 {
57  using namespace std;
58  using namespace JPP;
59 
60  string outputFile;
61  int numberOfPoints;
62  double epsilon;
63  int function;
64  set<double> D; // distance [m]
65  int debug;
66 
67  try {
68 
69  JParser<> zap("Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a shower.");
70 
71  zap['o'] = make_field(outputFile);
72  zap['n'] = make_field(numberOfPoints, "points for integration") = 25;
73  zap['e'] = make_field(epsilon, "precision for integration") = 1.0e-10;
74  zap['A'] = make_field(absorptionLengthFactor, "scaling factor") = 1.0;
75  zap['S'] = make_field(scatteringLengthFactor, "scaling factor") = 1.0;
76  zap['F'] = make_field(function, "PDF type") =
77  SCATTERED_LIGHT_FROM_MUON_5D,
78  DIRECT_LIGHT_FROM_EMSHOWER,
79  SCATTERED_LIGHT_FROM_EMSHOWER;
80  zap['D'] = make_field(D, "distance [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,
97  const double) const;
98 
99 
100  // set global parameters
101 
102  const double P_atm = NAMESPACE::getAmbientPressure();
103  const double wmin = getMinimalWavelength();
104  const double wmax = getMaximalWavelength();
105 
106 
107  const JPDF_C
108  pdf_c(NAMESPACE::getPhotocathodeArea(),
109  NAMESPACE::getQE,
110  NAMESPACE::getAngularAcceptance,
111  getAbsorptionLength,
112  getScatteringLength,
113  NAMESPACE::getScatteringProbability,
114  P_atm,
115  wmin,
116  wmax,
117  numberOfPoints,
118  epsilon);
119 
120 
121  typedef JSplineFunction1D_t JFunction1D_t;
122  typedef JMAPLIST<JPolint1FunctionalMap,
123  JPolint1FunctionalMap,
124  JPolint1FunctionalGridMap,
125  JPolint1FunctionalGridMap>::maplist JMapList_t;
126  typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
127 
128  typedef JPDFTransformer<4, JFunction1D_t::argument_type> JFunction4DTransformer_t;
129  typedef JArray<4, JFunction1D_t::argument_type> JArray_t;
130 
131  JPDF_t pdf;
132 
133 
134  NOTICE("building multi-dimensional function object <" << function << ">... " << flush);
135 
136  const double ng[] = { pdf_c.getIndexOfRefractionGroup(wmax),
137  pdf_c.getIndexOfRefractionGroup(wmin) };
138 
139  map<int, pair<fcn,JFunction4DTransformer_t> > zmap;
140 
141  zmap[SCATTERED_LIGHT_FROM_MUON_5D] = make_pair((fcn) &JPDF::getScatteredLightFromMuon, JFunction4DTransformer_t(21.5, 2, ng[0], 0.0, JGeant(JGeanx(0.33, -9.5)), 6e-4, NAMESPACE::getAngularAcceptance, 0.06));
142  zmap[DIRECT_LIGHT_FROM_EMSHOWER] = make_pair((fcn) &JPDF::getDirectLightFromEMshower, JFunction4DTransformer_t(21.5, 2, ng[0], ng[1], JGeant(JGeanx(0.35, -5.4)), 1e-5, NAMESPACE::getAngularAcceptance, 0.001));
143  zmap[SCATTERED_LIGHT_FROM_EMSHOWER] = make_pair((fcn) &JPDF::getScatteredLightFromEMshower, JFunction4DTransformer_t(21.5, 2, ng[0], 0.0, JGeant(JGeanx(0.55, -4.5)), 1e-2, NAMESPACE::getAngularAcceptance, 0.05));
144 
145  if (zmap.find(function) == zmap.end()) {
146  FATAL("illegal function specifier" << endl);
147  }
148 
149  fcn f = zmap[function].first; // PDF
150  JFunction4DTransformer_t transformer = zmap[function].second; // transformer
151 
152 
153  if (D.empty()) {
154  D.insert( 0.10);
155  D.insert( 0.50);
156  D.insert( 1.00);
157  D.insert( 5.00);
158  D.insert( 10.00);
159  D.insert( 20.00);
160  D.insert( 30.00);
161  D.insert( 40.00);
162  D.insert( 50.00);
163  D.insert( 60.00);
164  D.insert( 70.00);
165  D.insert( 80.00);
166  D.insert( 90.00);
167  D.insert(100.00);
168  D.insert(120.00);
169  D.insert(150.00);
170  D.insert(170.00);
171  D.insert(190.00);
172  D.insert(210.00);
173  D.insert(230.00);
174  D.insert(250.00);
175  D.insert(270.00);
176  D.insert(290.00);
177  D.insert(310.00);
178  }
179 
180  set<double> C; // cosine emission angle
181 
182  JQuadrature qeant(-1.0, +1.0, 60, geanx);
183 
184  for (JQuadrature::const_iterator i = qeant.begin(); i != qeant.end(); ++i)
185  C.insert(i->getX());
186 
187  C.insert(-1.00);
188  C.insert(+1.00);
189 
190 
191  set<double> X; // time [ns]
192 
193  if (function == DIRECT_LIGHT_FROM_EMSHOWER) {
194 
195  for (double buffer[] = { 0.0, 0.005, 0.01, 0.015, -1 }, *x = buffer; *x >= 0; ++x) {
196  X.insert(0.0 + *x);
197  X.insert(1.0 - *x);
198  }
199 
200  for (double x = 0.02; x < 0.99; x += 0.01)
201  X.insert(x);
202 
203  } else {
204 
205  X.insert( 0.00);
206  X.insert( 0.10);
207  X.insert( 0.20);
208  X.insert( 0.30);
209  X.insert( 0.40);
210  X.insert( 0.50);
211  X.insert( 0.60);
212  X.insert( 0.70);
213  X.insert( 0.80);
214  X.insert( 0.90);
215  X.insert( 1.00);
216  X.insert( 1.00);
217  X.insert( 1.10);
218  X.insert( 1.20);
219  X.insert( 1.30);
220  X.insert( 1.40);
221  X.insert( 1.50);
222  X.insert( 1.60);
223  X.insert( 1.70);
224  X.insert( 1.80);
225  X.insert( 1.90);
226  X.insert( 2.00);
227  X.insert( 2.20);
228  X.insert( 2.40);
229  X.insert( 2.60);
230  X.insert( 2.80);
231  X.insert( 3.00);
232  X.insert( 3.25);
233  X.insert( 3.50);
234  X.insert( 3.75);
235  X.insert( 4.00);
236  X.insert( 4.25);
237  X.insert( 4.50);
238  X.insert( 4.75);
239  X.insert( 5.0);
240  X.insert( 6.0);
241  X.insert( 7.0);
242  X.insert( 8.0);
243  X.insert( 9.0);
244  X.insert( 10.0);
245  X.insert( 15.0);
246  X.insert( 20.0);
247  X.insert( 25.0);
248  X.insert( 30.0);
249  X.insert( 40.0);
250  X.insert( 50.0);
251  X.insert( 60.0);
252  X.insert( 70.0);
253  X.insert( 80.0);
254  X.insert( 90.0);
255  X.insert(100.0);
256  X.insert(120.0);
257  X.insert(140.0);
258  X.insert(160.0);
259  X.insert(180.0);
260  X.insert(200.0);
261  X.insert(250.0);
262  X.insert(300.0);
263  X.insert(350.0);
264  X.insert(400.0);
265  X.insert(450.0);
266  X.insert(500.0);
267  X.insert(600.0);
268  X.insert(700.0);
269  X.insert(800.0);
270  X.insert(900.0);
271  X.insert(1200.0);
272  X.insert(1500.0);
273  }
274 
275  const double grid = 7.0; // [deg]
276 
277  const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0)); // azimuth angle unit step size
278 
279 
280  for (set<double>::const_iterator d = D.begin(); d != D.end(); ++d) {
281 
282  const double D_m = *d;
283 
284  for (set<double>::const_iterator c = C.begin(); c != C.end(); ++c) {
285 
286  const double cd = *c;
287 
288  const unsigned int number_of_theta_points = max(2u, (unsigned int) (180.0/(1.4 * grid)));
289 
290  for (double theta = 0.0; theta <= PI + epsilon; theta += PI/number_of_theta_points) {
291 
292  const unsigned int number_of_phi_points = max(2u, (unsigned int) (PI * sin(theta) / alpha));
293 
294  for (double phi = 0.0; phi <= PI + epsilon; phi += PI/number_of_phi_points) {
295 
296  JFunction1D_t& f1 = pdf[D_m][cd][theta][phi];
297 
298  const JArray_t array(D_m, cd, theta, phi);
299 
300  double t_old = transformer.getXn(array, *X.begin());
301  double y_old = 0.0;
302 
303  for (set<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
304 
305  const double t = transformer.getXn(array, *x);
306  const double y = (pdf_c.*f)(D_m, cd, theta, phi, t);
307 
308  if (y != 0.0) {
309 
310  if (*x < 0.0) {
311  WARNING("dt < 0 " << *x << ' ' << D_m << ' ' << t << ' ' << y << endl);
312  }
313 
314  if (y_old == 0.0) {
315  f1[t_old] = y_old;
316  }
317 
318  f1[t] = y;
319 
320  } else {
321 
322  if (y_old != 0.0) {
323  f1[t] = y;
324  }
325  }
326 
327  t_old = t;
328  y_old = y;
329  }
330  }
331  }
332  }
333  }
334 
335  pdf.transform(transformer);
336  pdf.compile();
337 
338  NOTICE("OK" << endl);
339 
340  try {
341 
342  NOTICE("storing output to file " << outputFile << "... " << flush);
343 
344  pdf.store(outputFile.c_str());
345 
346  NOTICE("OK" << endl);
347  }
348  catch(const JException& error) {
349  FATAL(error.what() << endl);
350  }
351 }

Variable Documentation

◆ absorptionLengthFactor

double absorptionLengthFactor

Scaling of absorption and scattering length.

Definition at line 24 of file JMakePDG.cc.

◆ scatteringLengthFactor

double scatteringLengthFactor

Definition at line 25 of file JMakePDG.cc.

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.
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
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER
direct light from EM shower
Definition: JPDFTypes.hh:40
scatteringLengthFactor
double scatteringLengthFactor
Definition: JMakePDG.cc:25
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
JTOOLS::C
static const double C
Speed of light in vacuum [m/ns].
Definition: JConstants.hh:22
std::set
Definition: JSTDTypes.hh:13
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
getAbsorptionLength
double getAbsorptionLength(const double lambda)
Definition: JMakePDG.cc:28
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
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER
scattered light from EM shower
Definition: JPDFTypes.hh:41
getScatteringLength
double getScatteringLength(const double lambda)
Definition: JMakePDG.cc:34
std::map
Definition: JSTDTypes.hh:16
JPHYSICS::SCATTERED_LIGHT_FROM_MUON_5D
scattered light from muon
Definition: JPDFTypes.hh:38
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
std
Definition: jaanetDictionary.h:36
JTOOLS::PI
static const double PI
Constants.
Definition: JConstants.hh:20
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::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
absorptionLengthFactor
double absorptionLengthFactor
Scaling of absorption and scattering length.
Definition: JMakePDG.cc:24
Generated by   doxygen 1.8.16