37 typedef JAbstractHistogram<double> JHistogram_t;
50 JParser<> zap(
"Program to plot PDF of Cherenkov light from EM-shower using interpolation tables.");
54 zap[
'D'] =
make_field(dir,
"(theta, phi) of PMT [rad]");
55 zap[
'x'] =
make_field(x,
"histogram x-binning") = JHistogram_t();
56 zap[
'y'] =
make_field(y,
"histogram y-binning") = JHistogram_t();
57 zap[
'E'] =
make_field(E,
"Energy [GeV]") = 0.0;
63 catch(
const exception &error) {
64 FATAL(error.what() << endl);
77 typedef JPolint0Function1D_t JFunction1D_t;
79 JMapList<JPolint1FunctionalMap,
80 JMapList<JPolint1FunctionalMap,
81 JMapList<JPolint1FunctionalGridMap,
82 JMapList<JPolint1FunctionalGridMap> > > > JMapList_t;
83 typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
84 typedef JNPETable<double, double, JMapList_t> JNPE_t;
86 JDistance<double>::precision = 1.0e-10;
88 const int N = inputFile.size();
95 for (
int i = 0; i != N; ++i) {
97 NOTICE(
"loading input from file " << inputFile[i] <<
"... " << flush);
99 pdf[i].load(inputFile[i].c_str());
101 pdf[i].setExceptionHandler(
new JFunction1D_t::JDefaultResult(
JMATH::zero));
103 npe[i] = JNPE_t(pdf[i]);
108 catch(
const JException& error) {
109 FATAL(error.what() << endl);
115 if (!x.is_valid()) { x = JHistogram_t(150, 0.0, 150.0); }
116 if (!y.is_valid()) { y = JHistogram_t(200, -1.0, +1.0); }
118 TH2D h0(
"h0",
"PDF Projection; D [m]; cos #theta_{0}",
119 x.getNumberOfBins(), x.getLowerLimit(), x.getUpperLimit(),
120 y.getNumberOfBins(), y.getLowerLimit(), y.getUpperLimit());
122 for (
int i = 1; i <= h0.GetNbinsX(); ++i) {
123 for (
int j = 1; j <= h0.GetNbinsY(); ++j) {
125 const double cd = h0.GetYaxis()->GetBinCenter(j);
126 const double D = h0.GetXaxis()->GetBinCenter(i);
132 const double W = 1.0 / (double) Z.size();
136 const double __D = sqrt(D*D - 2.0*(D*cd)*(*z) + (*z)*(*z));
137 const double __cd = (D * cd - (*z)) / __D;
139 for (
int k = 0; k != N; ++k) {
141 Y += W * npe[k](__D, __cd, dir.getTheta(), dir.getPhi());
143 catch(
const exception& error) {}
149 for (
int k = 0; k != N; ++k) {
150 Y += npe[k](D, cd, dir.getTheta(), dir.getPhi());
154 h0.SetBinContent(i, j, Y);
Utility class to parse command line options.
static const JZero zero
Function object to assign zero value.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
double getLength(const double E, const double P, const double eps=1.0e-3) const
Get shower length for a given integrated probability.