55 int main(
int argc,
char **argv)
69 JParser<> zap(
"Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a shower.");
73 zap[
'e'] =
make_field(epsilon,
"precision for integration") = 1.0e-10;
87 catch(
const exception &error) {
88 FATAL(error.what() << endl);
93 typedef double (JPDF::*fcn)(
const double,
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;
128 typedef JPDFTransformer<4, JFunction1D_t::argument_type> JFunction4DTransformer_t;
129 typedef JArray<4, JFunction1D_t::argument_type> JArray_t;
134 NOTICE(
"building multi-dimensional function object <" <<
function <<
">... " << flush);
136 const double ng[] = { pdf_c.getIndexOfRefractionGroup(wmax),
137 pdf_c.getIndexOfRefractionGroup(wmin) };
145 if (zmap.find(
function) == zmap.end()) {
146 FATAL(
"illegal function specifier" << endl);
149 fcn f = zmap[
function].first;
150 JFunction4DTransformer_t transformer = zmap[
function].second;
182 JQuadrature qeant(-1.0, +1.0, 60,
geanx);
184 for (JQuadrature::const_iterator i = qeant.begin(); i != qeant.end(); ++i)
195 for (
double buffer[] = { 0.0, 0.005, 0.01, 0.015, -1 }, *x = buffer; *x >= 0; ++x) {
200 for (
double x = 0.02; x < 0.99; x += 0.01)
275 const double grid = 7.0;
277 const double alpha = 2.0 * sqrt(1.0 - cos(grid *
PI / 180.0));
282 const double D_m = *d;
286 const double cd = *c;
288 const unsigned int number_of_theta_points = max(2
u, (
unsigned int) (180.0/(1.4 * grid)));
290 for (
double theta = 0.0; theta <=
PI + epsilon; theta +=
PI/number_of_theta_points) {
292 const unsigned int number_of_phi_points = max(2
u, (
unsigned int) (
PI * sin(theta) / alpha));
294 for (
double phi = 0.0; phi <=
PI + epsilon; phi +=
PI/number_of_phi_points) {
296 JFunction1D_t& f1 = pdf[D_m][cd][theta][phi];
298 const JArray_t array(D_m, cd, theta, phi);
300 double t_old = transformer.getXn(array, *X.begin());
305 const double t = transformer.getXn(array, *x);
306 const double y = (pdf_c.*f)(D_m, cd, theta, phi, t);
311 WARNING(
"dt < 0 " << *x <<
' ' << D_m <<
' ' << t <<
' ' << y << endl);
335 pdf.transform(transformer);
348 catch(
const JException& error) {
349 FATAL(error.what() << endl);