Add shower pdf

CHANGELOG.rst

@@ -3,6 +3,10 @@ Unreleased changes
 
 Version 3
 ---------
+3.2 / 2020-03-02
+~~~~~~~~~~~~~~~~
+* shower PDF functionality added (``jppy.pdf.JShowerPDF``)
+
 3.1 / 2020-02-28
 ~~~~~~~~~~~~~~~~
 * removed ROOT dependency. ``jppy`` now only depends on a few Jpp headers

scripts/get_pdfs.sh

 #!/usr/bin/env bash
+
+export URL="http://pi1139.physik.uni-erlangen.de/data/latest"
 if [ ! -d "pdfs" ]; then
     echo "Retrieving PDFs..."
     mkdir pdfs
     cd pdfs
     for i in {1..6}; do
-        wget "http://pi1139.physik.uni-erlangen.de/data/latest/J${i}p.dat"
+        wget "$URL/J${i}p.dat"
+    done
+    for i in {13..14}; do
+        wget "$URL/J${i}p.dat"
     done
 else
     echo "PDFs already downloaded."

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_<JShowerPDF_t>(m, "JShowerPDF")
+        .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", &JShowerPDF_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"),

tests/test_pdf.py

@@ -11,3 +11,12 @@ class TestMuonPDF(unittest.TestCase):
         self.assertAlmostEqual(-2.22809691e-05, result.fp)
         self.assertAlmostEqual(0.024902763, result.v)
         self.assertAlmostEqual(0.116492968, result.V)
+
+class TestShowerPDF(unittest.TestCase):
+    def test_pdf(self):
+        shower_pdf = jppy.pdf.JShowerPDF(PDFS, 0)
+        result = shower_pdf.calculate(100, 10, 0.1, 0.2, 0.3, 6)
+        self.assertAlmostEqual(0.001612553295068934, result.f)
+        self.assertAlmostEqual(4.000285659029551e-05, result.fp)
+        self.assertAlmostEqual(0.010999553987301543, result.v)
+        self.assertAlmostEqual(0.1527856994106781, result.V)