Merge branch 'add_pdf_evaluator' into 'master'

Add geane and pdf evaluator

See merge request !7

src/jppy/__init__.py

+from pkg_resources import get_distribution, DistributionNotFound
+try:
+    version = get_distribution(__name__).version
+except DistributionNotFound:
+    version = "unknown version"
+
+    
+from . import pdf_evaluator
+from . import constants
+from . import geane
+from . import pdf
+from . import npe

src/jppy/constants.cc

+#include <pybind11/pybind11.h>
+
+#include "JPhysics/JConstants.hh"
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(constants, m) {
+  m.doc() = "Jpp constants";
+  m.def("get_speed_of_light", &JPHYSICS::getSpeedOfLight);
+  m.def("get_inverse_speed_of_light", &JPHYSICS::getInverseSpeedOfLight);
+  m.def("get_index_of_refraction", &JPHYSICS::getIndexOfRefraction);
+  m.def("get_index_of_refraction_phase", &JPHYSICS::getIndexOfRefractionPhase);
+  m.def("get_tan_theta_c", &JPHYSICS::getTanThetaC);
+  m.def("get_cos_theta_c", &JPHYSICS::getCosThetaC);
+  m.def("get_sin_theta_c", &JPHYSICS::getSinThetaC);
+  m.def("get_kappa_c", &JPHYSICS::getKappaC);
+}
+

src/jppy/geane.cc

+#include <pybind11/pybind11.h>
+
+#include "JPhysics/JGeane.hh"
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(geane, m) {
+  m.doc() = "Utilities for muon energy losses";
+  m.def("geanc", &JPHYSICS::geanc);
+  py::class_<JPHYSICS::JGeane>(m, "JGeane");
+  py::class_<JPHYSICS::JGeane_t, JPHYSICS::JGeane>(m, "JGeane_t")
+    .def(py::init<const double, const double>(),
+	 py::arg("Energy loss due to ionisation [GeV/m]"),
+	 py::arg("Energy loss due to pair production and Bremsstrahlung [m^-1]"))
+    .def("get_a", &JPHYSICS::JGeane_t::getA)
+    .def("get_b", &JPHYSICS::JGeane_t::getB)
+    .def("get_E", &JPHYSICS::JGeane_t::getE,
+	 py::arg("Energy of muon [GeV]"),
+	 py::arg("Distance traveled [m]"))
+    .def("get_X", &JPHYSICS::JGeane_t::getX,
+	 py::arg("Energy of muon at start [GeV]"),
+	 py::arg("Energy of mun at end [GeV]")
+	 );
+  py::class_<JPHYSICS::JGeaneWater, JPHYSICS::JGeane>(m, "JGeaneWater")
+    .def(py::init<>())
+    .def("get_a", &JPHYSICS::JGeaneWater::getA)
+    .def("get_b", &JPHYSICS::JGeaneWater::getB)
+    .def("get_E", &JPHYSICS::JGeaneWater::getE,
+	 py::arg("Energy of muon [GeV]"),
+	 py::arg("Distance traveled [m]"))
+    .def("get_X", &JPHYSICS::JGeaneWater::getX,
+	 py::arg("Energy of muon at start [GeV]"),
+	 py::arg("Energy of mun at end [GeV]")
+	 );
+}

src/jppy/pdf_evaluator.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Auxiliary python wrappers for Jpp PDFs.
+The wrapper classes provide a common interface to evaluate the muon and shower PDF values\n
+on an event-by-event basis.
+"""
+
+from math import sqrt
+
+from abc import ABCMeta, abstractmethod
+
+from jppy import constants as const
+from jppy.geane import JGeaneWater as JGeaneWater
+from jppy.pdf import (
+    JMuonPDF as JMuonPDF,
+    JShowerPDF as JShowerPDF,
+)
+
+class PDF(object, metaclass=ABCMeta):
+    """Abstract base class for event-by-event PDF evaluation"""
+    
+    def __init__(self, energy=0.0, t0=0.0):
+        self._energy = energy
+        self._t0 = t0
+
+    @classmethod
+    def __subclasshook__(cls, subclass):
+        return (hasattr(subclass, 'energy') and
+                callable(subclass.energy) and
+                hasattr(subclass, 't0') and
+                callable(subclass.t0) and
+                hasattr(subclass, 'evaluate') and
+                callable(subclass.evaluate))
+    
+    @property
+    def energy(self):
+        return self._energy
+
+    @energy.setter
+    def energy(self, value):
+        self._energy = float(value)
+
+    @property
+    def t0(self):
+        return self._t0
+
+    @t0.setter
+    def t0(self, value):
+        self._t0 = float(t0)
+
+    @abstractmethod
+    def evaluate(self, D, cd, theta, phi, t_obs):
+        raise NotImplementedError
+
+        
+class MuonPDF(PDF):
+    """Muon PDF evaluator"""
+
+    def __init__(self, PDFS, energy=0.0, t0=0.0, TTS=0.0):
+        """
+        Constructor.
+        
+        Parameters
+        ----------
+        energy : float
+            muon energy at the simulated vertex 
+            (i.e. interaction vertex in the case of a neutrino interaction\n
+             or the vertex corresponding to the can interception for atmospheric muons)
+        t0 : float
+            time corresponding to the simulated vertex [ns]
+            (i.e. interaction vertex in the case of a neutrino interaction\n
+             or the vertex corresponding to the can interception for atmospheric muons)
+        TTS : float
+            transit time spread [ns]
+        """
+
+        super().__init__(energy, t0)
+        
+        self._geane = JGeaneWater()
+        self._pdf = JMuonPDF(PDFS, TTS=TTS)
+
+    def evaluate(self, D, cd, theta, phi, t_obs):
+        """
+        Muon PDF evaluation method.
+
+        Parameters
+        ----------
+        D : array[float], shape=(n,)
+            Hit distances with respect to the simulated vertex\n
+            (i.e. interaction vertex in the case of a neutrino interaction\n
+             or the vertex corresponding to the can interception for atmospheric muons)
+        cd : array[float], shape=(n,)
+            angle between muon direction and PMT position
+        theta : array[float], shape=(n,)
+            PMT longitudinal angle [deg]
+        phi : array[float], shape=(n,)
+            PMT azimuthal angle [deg]
+        t_obs : array[float], shape=(n,)
+            Observed/Simulated hit times
+
+        Returns
+        -------
+        muon pdf values : array[float], shape=(n,)
+        """
+        
+        dz = D * cd
+        R = sqrt((D + dz) * (D - dz))
+
+        E = self._geane.get_E(self.energy, dz)
+        
+        t_exp = self.t0 + (dz + R * const.get_kappa_c()) * const.get_inverse_speed_of_light()
+        dt = t_obs - t_exp
+        
+        return self._pdf.calculate(E, R, theta, phi, dt)
+    
+
+class ShowerPDF(PDF):
+    """Shower PDF evaluator"""
+
+    def __init__(self, PDFS, energy=0.0, t0=0.0, TTS=0.0):
+        """
+        Constructor.
+        
+        Parameters
+        ----------
+        energy : float
+            shower energy [GeV]
+        t0 : float
+            time corresponding to shower vertex [ns]
+        TTS : float
+            transit time spread [ns]
+        """
+
+        super().__init__(energy, t0)
+        
+        self._pdf = JShowerPDF(PDFS, TTS=TTS)
+
+    def evaluate(self, D, cd, theta, phi, t_obs):
+        """
+        Shower PDF evaluation method.
+
+        Parameters
+        ----------
+        D : array[float], shape=(n,)
+            Hit distances with respect to the shower vertex
+        cd : array[float], shape=(n,)
+            angle between shower direction and PMT position
+        theta : array[float], shape=(n,)
+            PMT longitudinal angle [deg]
+        phi : array[float], shape=(n,)
+            PMT azimuthal angle [deg]
+        t_obs : array[float], shape=(n,)
+            Observed/Simulated hit times
+
+        Returns
+        -------
+        shower pdf values : array[float], shape=(n,)
+        """
+        
+        t_exp = self.t0 + D * const.get_inverse_speed_of_light() * const.get_index_of_refraction()
+        dt = t_obs - t_exp
+        
+        return self._pdf.calculate(self.energy, D, cd, theta, phi, dt)

tests/test_geane.py

+import unittest
+import jppy
+
+class TestGeaneWater(unittest.TestCase):
+    def test_geane(self):
+        gwater = jppy.geane.JGeaneWater()
+        density_sea_water = 1.038
+        assert(gwater.get_a() == 2.30e-1 * density_sea_water)
+        assert(gwater.get_b() == 3.40e-4 * density_sea_water)
+        self.assertAlmostEqual(gwater.get_E(4e4, 100), 3.857507637293732e+04)
+        self.assertAlmostEqual(gwater.get_X(4e4, 4e3), 6.069985857980293e+03)

tests/test_pdf.py

@@ -7,16 +7,16 @@ class TestMuonPDF(unittest.TestCase):
     def test_pdf(self):
         muon_pdf = jppy.pdf.JMuonPDF(PDFS, 0)
         result = muon_pdf.calculate(10, 5, 0, 0, 23)
-        self.assertAlmostEqual(0.001535784, result.f)
-        self.assertAlmostEqual(-2.22809691e-05, result.fp)
-        self.assertAlmostEqual(0.024902763, result.v)
-        self.assertAlmostEqual(0.116492968, result.V)
+        self.assertAlmostEqual(0.001545692, result.f)
+        self.assertAlmostEqual(-1.220889709e-05, result.fp)
+        self.assertAlmostEqual(0.022764524, result.v)
+        self.assertAlmostEqual(0.115814468, 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)
+        self.assertAlmostEqual(0.001383343, result.f)
+        self.assertAlmostEqual(5.091930537e-05, result.fp)
+        self.assertAlmostEqual(0.010475308, result.v)
+        self.assertAlmostEqual(0.149635554, result.V)

tests/test_pdf_evaluator.py

+import unittest
+import jppy
+
+PDFS = "pdfs/J%p.dat"
+
+class TestMuonPDFEvaluator(unittest.TestCase):
+    def test_pdf_evaluator(self):
+        E, t0, t_obs, D, cd, theta, phi = [1e3, 56, 292, 50, 0.7, 1.57, 3.14]
+        muon_pdf = jppy.pdf_evaluator.MuonPDF(PDFS, E, t0)
+        result = muon_pdf.evaluate(D, cd, theta, phi, t_obs)
+        self.assertAlmostEqual( 0.003644475, result.f)
+        self.assertAlmostEqual(-0.000715558, result.fp)
+        self.assertAlmostEqual( 0.033748905, result.v)
+        self.assertAlmostEqual( 0.097171157, result.V)
+
+class TestShowerPDFEvaluator(unittest.TestCase):
+    def test_pdf_evaluator(self):
+        E, t0, t_obs, D, cd, theta, phi = [50, 198, 226, 5, 0.6, 0.5, 0.4]
+        shower_pdf = jppy.pdf_evaluator.ShowerPDF(PDFS, E, t0)
+        result = shower_pdf.evaluate(D, cd, theta, phi, t_obs)
+        self.assertAlmostEqual(0.006204617, result.f)
+        self.assertAlmostEqual(0.000641669, result.fp)
+        self.assertAlmostEqual(0.013960066, result.v)
+        self.assertAlmostEqual(0.296589983, result.V)
+
+
+