Add shower pdf

src/pdf.cc

@@ -148,6 +148,126 @@ struct JMuonPDF_t {
 };
 
 
+/**
+ * Auxiliary data structure for shower PDF.
+ */
+struct JShowerPDF_t {
+
+  typedef JPP::JSplineFunction1D<JPP::JSplineElement2S<double, double>,
+				 JPP::JCollection,
+				 JPP::JResultPDF<double> >         JFunction1D_t;
+
+  typedef JPP::JMAPLIST<JPP::JPolint1FunctionalMap,
+			JPP::JPolint1FunctionalMap,
+			JPP::JPolint0FunctionalGridMap,
+			JPP::JPolint0FunctionalGridMap>::maplist   JPDFMaplist_t;
+  typedef JPP::JPDFTable<JFunction1D_t, JPDFMaplist_t>             JPDF_t;
+  typedef JFunction1D_t::result_type                               result_type;
+
+
+  /**
+   * Constructor.
+   *
+   * The PDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD.\n
+   * The <tt>TTS</tt> corresponds to the additional time smearing applied to the PDFs.
+   *
+   * \param  fileDescriptor     PDF file descriptor
+   * \param  TTS                TTS [ns]
+   * \param  numberOfPoints     number of points for Gauss-Hermite integration of TTS
+   * \param  epsilon            precision        for Gauss-Hermite integration of TTS
+   */
+  JShowerPDF_t(const std::string& fileDescriptor,
+	       const double       TTS,
+	       const int          numberOfPoints = 25,
+	       const double       epsilon        = 1.0e-10)
+  {
+    using namespace std;
+    using namespace JPP;
+
+    const JPDFType_t pdf_t[] = { SCATTERED_LIGHT_FROM_EMSHOWER,
+				 DIRECT_LIGHT_FROM_EMSHOWER };
+
+    const  int N = sizeof(pdf_t) / sizeof(pdf_t[0]);
+
+    const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(zero));
+
+    for (int i = 0; i != N; ++i) {
+
+      const string file_name = getFilename(fileDescriptor, pdf_t[i]);
+
+      cout << "loading input from file " << file_name << "... " << flush;
+
+      JPDF_t pdf;
+
+      pdf.load(file_name.c_str());
+
+      pdf.setExceptionHandler(supervisor);
+
+      if (pdfA.empty())
+	pdfA = pdf;
+      else
+	pdfA.add(pdf);
+
+      cout << "OK" << endl;
+    }
+
+    if        (TTS > 0.0) {
+
+      cout << "bluring PDFs... " << flush;
+
+      pdfA.blur(TTS, numberOfPoints, epsilon);
+
+      cout << "OK" << endl;
+
+    } else if (TTS < 0.0) {
+
+      THROW(JValueOutOfRange, "Illegal value of TTS [ns]: " << TTS);
+    }
+  }
+
+
+  /**
+   * Get PDF.
+   *
+   * The orientation of the PMT should be defined according this <a href="https://common.pages.km3net.de/jpp/JPDF.PDF">documentation</a>.\n
+   * In this, the zenith and azimuth angles are limited to \f[\left[0, \pi\right]\f].
+   *
+   * \param  E                  shower energy at minimum distance of approach [GeV]
+   * \param  D                  distance [m]
+   * \param  cd                 cosine emission angle
+   * \param  theta              PMT zenith  angle [rad]
+   * \param  phi                PMT azimuth angle [rad]
+   * \param  t1                 arrival time relative to Cherenkov hypothesis [ns]
+   * \return                    hypothesis value
+   */
+  result_type calculate()(const double E,
+			 const double D,
+			 const double cd,
+			 const double theta,
+			 const double phi,
+			 const double t1) const
+  {
+    using namespace JPP;
+
+    result_type h1 = pdfA(D, cd, theta, phi, t1) * E;
+
+    // safety measures
+
+    if (h1.f <= 0.0) {
+      h1.f  = 0.0;
+      h1.fp = 0.0;
+    }
+
+    if (h1.v <= 0.0) {
+      h1.v  = 0.0;
+    }
+
+    return h1;
+  }
+
+  JPDF_t pdfA;   //!< PDF for shower
+};
+
 /**
    Python bindings.
  */
@@ -169,6 +289,20 @@ PYBIND11_MODULE(pdf, m) {
              py::arg("phi"),
              py::arg("t1")
             ),
+    py::class_<JShoweverPDF_t>(m, "JShoweerPDF")
+        .def(py::init<const std::string &, double, int, double>(),
+             py::arg("file_descriptor"),
+             py::arg("TTS"),
+             py::arg("number_of_points") = 25,
+             py::arg("epsilon") = 1e-10)
+        .def("calculate", &JMuonPDF_t::calculate,
+             py::arg("E"),
+             py::arg("D"),
+             py::arg("cd"),
+             py::arg("theta"),
+             py::arg("phi"),
+             py::arg("t1")
+            ),
     py::class_<JTOOLS::JResultPDF<double>>(m, "JResultPDF")
         .def(py::init<double, double, double, double>(),
              py::arg("f"),