diff --git a/.gitignore b/.gitignore
index 282d5df620c92157dec51f058ba232b07e0735c1..f1eede293f7c629e65692f7ae448d0e05974314f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,12 +1,27 @@
-*.so
+# Byte-compiled / optimized / DLL files
+__pycache__/
 *.py[cod]
-*.egg-info
+*.pyxbldc
+
+# C extensions
+*.so
+
+# Distribution / packaging
+build/
+venv/
+dist/
 _build
 _generate
-build
-venv
-jppy/__init__.py
-jppy/version.py
 pip-wheel-metadata
-scratch
-pdfs
+*.egg-info
+*.egg
+
+# Version info
+src/jppy/version.py
+
+# Temporary folders and meta-data
+tmp/
+
+# Data folders
+pdfs/
+
diff --git a/pyproject.toml b/pyproject.toml
index a0a6409d121e6d9da585868d7f0fed0d3dcff874..fd2ba36859211fca406d55eb4d6d7c59542d549f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,5 +1,6 @@
 [build-system]
-requires = ["setuptools>=42", "wheel", "setuptools_scm[toml]>=3.4"]
+requires = ["setuptools>=42", "wheel", "setuptools_scm[toml]>=3.4", "pybind11>=2.6"]
+build-backend = "setuptools.build_meta"
 
 [tool.setuptools_scm]
-write_to = "jppy/version.py"
+write_to = "src/jppy/version.py"
diff --git a/scripts/get_pdfs.sh b/scripts/get_pdfs.sh
index d1849150f7618383514f8cca70e74ea58780d3de..9dd954d2dc562fa3e986c49b75dbc5a2bc5fd495 100755
--- a/scripts/get_pdfs.sh
+++ b/scripts/get_pdfs.sh
@@ -1,6 +1,6 @@
 #!/usr/bin/env bash
 
-export URL="http://pi1139.physik.uni-erlangen.de/data/latest"
+export URL="http://sftp.km3net.de/data/latest/"
 if [ ! -d "pdfs" ]; then
     echo "Retrieving PDFs..."
     mkdir pdfs
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..2aa5d31fb329d8973b7122f22a5cf620ccaea085
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,91 @@
+[metadata]
+name = jppy
+description = "Jpp Python package"
+long_description = file: README.rst
+long_description_content_type = text/x-rst
+url = https://git.km3net.de/km3py/jppy
+author = Tamas Gal
+author_email = tgal@km3net.de
+maintainer = Tamas Gal
+maintainer_email = tgal@km3net.de
+license = MIT
+license_files = LICENSE
+classifiers =
+    Development Status :: 5 - Production/Stable
+    Intended Audience :: Developers
+    Intended Audience :: Science/Research
+    License :: OSI Approved :: BSD License
+    Operating System :: OS Independent
+    Programming Language :: C++
+    Programming Language :: Python
+    Programming Language :: Python :: 3
+    Programming Language :: Python :: 3 :: Only
+    Programming Language :: Python :: 3.6
+    Programming Language :: Python :: 3.7
+    Programming Language :: Python :: 3.8
+    Programming Language :: Python :: 3.9
+    Programming Language :: Python :: 3.10
+    Topic :: Scientific/Engineering
+keywords =
+    C++11
+    Python bindings
+    KM3NeT
+    muon
+    neutrino
+    shower
+    hadronic
+    electromagnetic
+    Photomultiplier tubes
+    probability density function
+
+[options]
+packages = find:
+install_requires =
+    pybind11>=2.4
+    setuptools>=40.6.2
+    setuptools_scm
+    toml
+python_requires = >=3.6
+include_package_data = True
+package_dir =
+    =src
+
+[build_ext]
+inplace=1
+
+[options.packages.find]
+where = src
+
+[options.extras_require]
+all =
+    ipykernel
+dev =
+    ipykernel
+    pytest
+    pytest-cov
+    pytest-flake8
+    pytest-pylint
+    pytest-watch
+    sphinx-gallery>=0.1.12
+    sphinx-rtd-theme>=0.3
+    sphinx>=1.6.3
+    sphinxcontrib-napoleon>=0.6.1
+    sphinxcontrib-programoutput>=0.11
+    sphinxcontrib-websupport>=1.0.1
+    sphinx-autoapi
+    twine
+    wheel
+
+[bdist_wheel]
+universal = 1
+
+[tool:pytest]
+junit_family = xunit2
+addopts = -vv -rs -Wd
+testpaths =
+    tests
+
+[check-manifest]
+ignore =
+    src/version.py
+
diff --git a/setup.py b/setup.py
index 210b5ee582d720858a86bef0be460951f50604bb..29a355f6f98ac2eebe359ac915510b60c12d7442 100644
--- a/setup.py
+++ b/setup.py
@@ -1,3 +1,4 @@
+#!/usr/bin/env python3
 import os
 from setuptools import setup, Extension
 from setuptools.command.build_ext import build_ext
@@ -19,38 +20,13 @@ class get_pybind_include(object):
         import pybind11
         return pybind11.get_include(self.user)
 
-
+    
 def get_jpp_include():
     jpp_dir = os.getenv("JPP_DIR")
     if jpp_dir:
         return os.path.join(jpp_dir, "software")
-    return "jpp"
-
+    return "src/jpp"
 
-ext_modules = [
-    Extension(
-        'jppy.{}'.format(module),
-        ['src/{}.cc'.format(module)],
-        include_dirs=[
-            get_pybind_include(),
-            get_pybind_include(user=True),
-            get_jpp_include()
-        ],
-        language='c++') for module in ['pdf', 'npe']
-]
-
-# Populating the __init__.py with submodule imports, so that one can import
-# the package and use the submodules directly with the dot-sytax.
-with open("jppy/__init__.py", "w") as fobj:
-    fobj.write("""from pkg_resources import get_distribution, DistributionNotFound
-try:
-    version = get_distribution(__name__).version
-except DistributionNotFound:
-    version = "unknown version"
-""")
-    for module in ext_modules:
-        fobj.write("from . import {}\n".format(module.name.split('.')[1]))
-        
 
 # As of Python 3.6, CCompiler has a `has_flag` method.
 # cf http://bugs.python.org/issue26689
@@ -117,18 +93,23 @@ class BuildExt(build_ext):
         build_ext.build_extensions(self)
 
 
-setup(
-    name='jppy',
-    author='Tamas Gal',
-    author_email='tgal@km3net.de',
-    url='https://git.km3net.de/km3py/jppy',
-    description='Jpp Python Package',
-    packages=["jppy"],
-    long_description="jppy - Jpp Python Package",
-    ext_modules=ext_modules,
-    install_requires=['pybind11>=2.4'],
-    setup_requires=['pybind11>=2.4', 'setuptools_scm'],
-    use_scm_version=True,
-    cmdclass={'build_ext': BuildExt},
-    zip_safe=False,
-)
+if __name__ == '__main__':    
+
+    setup_args = dict(
+        ext_modules = [
+            Extension(
+                'jppy.{}'.format(module),
+                ['src/jppy/{}.cc'.format(module)],
+                include_dirs=[
+                    get_pybind_include(),
+                    get_pybind_include(user=True),
+                    get_jpp_include()
+                ],
+                language='c++') for module in ['constants', 'geane', 'pdf', 'npe']
+        ],
+        cmdclass = dict(
+            build_ext = BuildExt
+        )
+    )
+
+    setup(**setup_args)
diff --git a/jpp/JIO/JBufferedIO.hh b/src/jpp/JIO/JBufferedIO.hh
similarity index 100%
rename from jpp/JIO/JBufferedIO.hh
rename to src/jpp/JIO/JBufferedIO.hh
diff --git a/jpp/JIO/JFileStreamIO.hh b/src/jpp/JIO/JFileStreamIO.hh
similarity index 100%
rename from jpp/JIO/JFileStreamIO.hh
rename to src/jpp/JIO/JFileStreamIO.hh
diff --git a/jpp/JIO/JObjectBinaryIO.hh b/src/jpp/JIO/JObjectBinaryIO.hh
similarity index 100%
rename from jpp/JIO/JObjectBinaryIO.hh
rename to src/jpp/JIO/JObjectBinaryIO.hh
diff --git a/jpp/JIO/JSerialisable.hh b/src/jpp/JIO/JSerialisable.hh
similarity index 100%
rename from jpp/JIO/JSerialisable.hh
rename to src/jpp/JIO/JSerialisable.hh
diff --git a/jpp/JIO/JStreamIO.hh b/src/jpp/JIO/JStreamIO.hh
similarity index 100%
rename from jpp/JIO/JStreamIO.hh
rename to src/jpp/JIO/JStreamIO.hh
diff --git a/jpp/JLang/JAbstractObjectStatus.hh b/src/jpp/JLang/JAbstractObjectStatus.hh
similarity index 100%
rename from jpp/JLang/JAbstractObjectStatus.hh
rename to src/jpp/JLang/JAbstractObjectStatus.hh
diff --git a/jpp/JLang/JAbstractPointer.hh b/src/jpp/JLang/JAbstractPointer.hh
similarity index 100%
rename from jpp/JLang/JAbstractPointer.hh
rename to src/jpp/JLang/JAbstractPointer.hh
diff --git a/jpp/JLang/JAssert.hh b/src/jpp/JLang/JAssert.hh
similarity index 100%
rename from jpp/JLang/JAssert.hh
rename to src/jpp/JLang/JAssert.hh
diff --git a/jpp/JLang/JBinaryIO.hh b/src/jpp/JLang/JBinaryIO.hh
similarity index 100%
rename from jpp/JLang/JBinaryIO.hh
rename to src/jpp/JLang/JBinaryIO.hh
diff --git a/jpp/JLang/JBool.hh b/src/jpp/JLang/JBool.hh
similarity index 100%
rename from jpp/JLang/JBool.hh
rename to src/jpp/JLang/JBool.hh
diff --git a/jpp/JLang/JCC.hh b/src/jpp/JLang/JCC.hh
similarity index 100%
rename from jpp/JLang/JCC.hh
rename to src/jpp/JLang/JCC.hh
diff --git a/jpp/JLang/JClass.hh b/src/jpp/JLang/JClass.hh
similarity index 100%
rename from jpp/JLang/JClass.hh
rename to src/jpp/JLang/JClass.hh
diff --git a/jpp/JLang/JClonable.hh b/src/jpp/JLang/JClonable.hh
similarity index 100%
rename from jpp/JLang/JClonable.hh
rename to src/jpp/JLang/JClonable.hh
diff --git a/jpp/JLang/JComparable.hh b/src/jpp/JLang/JComparable.hh
similarity index 100%
rename from jpp/JLang/JComparable.hh
rename to src/jpp/JLang/JComparable.hh
diff --git a/jpp/JLang/JEquals.hh b/src/jpp/JLang/JEquals.hh
similarity index 100%
rename from jpp/JLang/JEquals.hh
rename to src/jpp/JLang/JEquals.hh
diff --git a/jpp/JLang/JException.hh b/src/jpp/JLang/JException.hh
similarity index 100%
rename from jpp/JLang/JException.hh
rename to src/jpp/JLang/JException.hh
diff --git a/jpp/JLang/JForwardIterator.hh b/src/jpp/JLang/JForwardIterator.hh
similarity index 100%
rename from jpp/JLang/JForwardIterator.hh
rename to src/jpp/JLang/JForwardIterator.hh
diff --git a/jpp/JLang/JLangToolkit.hh b/src/jpp/JLang/JLangToolkit.hh
similarity index 100%
rename from jpp/JLang/JLangToolkit.hh
rename to src/jpp/JLang/JLangToolkit.hh
diff --git a/jpp/JLang/JManip.hh b/src/jpp/JLang/JManip.hh
similarity index 100%
rename from jpp/JLang/JManip.hh
rename to src/jpp/JLang/JManip.hh
diff --git a/jpp/JLang/JMemory.hh b/src/jpp/JLang/JMemory.hh
similarity index 100%
rename from jpp/JLang/JMemory.hh
rename to src/jpp/JLang/JMemory.hh
diff --git a/jpp/JLang/JNullType.hh b/src/jpp/JLang/JNullType.hh
similarity index 100%
rename from jpp/JLang/JNullType.hh
rename to src/jpp/JLang/JNullType.hh
diff --git a/jpp/JLang/JObjectID.hh b/src/jpp/JLang/JObjectID.hh
similarity index 100%
rename from jpp/JLang/JObjectID.hh
rename to src/jpp/JLang/JObjectID.hh
diff --git a/jpp/JLang/JObjectIO.hh b/src/jpp/JLang/JObjectIO.hh
similarity index 100%
rename from jpp/JLang/JObjectIO.hh
rename to src/jpp/JLang/JObjectIO.hh
diff --git a/jpp/JLang/JPointer.hh b/src/jpp/JLang/JPointer.hh
similarity index 100%
rename from jpp/JLang/JPointer.hh
rename to src/jpp/JLang/JPointer.hh
diff --git a/jpp/JLang/JSTDTypes.hh b/src/jpp/JLang/JSTDTypes.hh
similarity index 100%
rename from jpp/JLang/JSTDTypes.hh
rename to src/jpp/JLang/JSTDTypes.hh
diff --git a/jpp/JLang/JSharedCounter.hh b/src/jpp/JLang/JSharedCounter.hh
similarity index 100%
rename from jpp/JLang/JSharedCounter.hh
rename to src/jpp/JLang/JSharedCounter.hh
diff --git a/jpp/JLang/JSharedPointer.hh b/src/jpp/JLang/JSharedPointer.hh
similarity index 100%
rename from jpp/JLang/JSharedPointer.hh
rename to src/jpp/JLang/JSharedPointer.hh
diff --git a/jpp/JLang/JSinglePointer.hh b/src/jpp/JLang/JSinglePointer.hh
similarity index 100%
rename from jpp/JLang/JSinglePointer.hh
rename to src/jpp/JLang/JSinglePointer.hh
diff --git a/jpp/JLang/JStorage.hh b/src/jpp/JLang/JStorage.hh
similarity index 100%
rename from jpp/JLang/JStorage.hh
rename to src/jpp/JLang/JStorage.hh
diff --git a/jpp/JLang/JType.hh b/src/jpp/JLang/JType.hh
similarity index 100%
rename from jpp/JLang/JType.hh
rename to src/jpp/JLang/JType.hh
diff --git a/jpp/JLang/JVectorize.hh b/src/jpp/JLang/JVectorize.hh
similarity index 100%
rename from jpp/JLang/JVectorize.hh
rename to src/jpp/JLang/JVectorize.hh
diff --git a/jpp/JLang/JVoid.hh b/src/jpp/JLang/JVoid.hh
similarity index 100%
rename from jpp/JLang/JVoid.hh
rename to src/jpp/JLang/JVoid.hh
diff --git a/jpp/JLang/gzstream.h b/src/jpp/JLang/gzstream.h
similarity index 100%
rename from jpp/JLang/gzstream.h
rename to src/jpp/JLang/gzstream.h
diff --git a/jpp/JMath/JCalculator.hh b/src/jpp/JMath/JCalculator.hh
similarity index 100%
rename from jpp/JMath/JCalculator.hh
rename to src/jpp/JMath/JCalculator.hh
diff --git a/jpp/JMath/JConstants.hh b/src/jpp/JMath/JConstants.hh
similarity index 100%
rename from jpp/JMath/JConstants.hh
rename to src/jpp/JMath/JConstants.hh
diff --git a/jpp/JMath/JLimits.hh b/src/jpp/JMath/JLimits.hh
similarity index 100%
rename from jpp/JMath/JLimits.hh
rename to src/jpp/JMath/JLimits.hh
diff --git a/jpp/JMath/JMath.hh b/src/jpp/JMath/JMath.hh
similarity index 100%
rename from jpp/JMath/JMath.hh
rename to src/jpp/JMath/JMath.hh
diff --git a/jpp/JMath/JMathSupportkit.hh b/src/jpp/JMath/JMathSupportkit.hh
similarity index 100%
rename from jpp/JMath/JMathSupportkit.hh
rename to src/jpp/JMath/JMathSupportkit.hh
diff --git a/jpp/JMath/JZero.hh b/src/jpp/JMath/JZero.hh
similarity index 100%
rename from jpp/JMath/JZero.hh
rename to src/jpp/JMath/JZero.hh
diff --git a/jpp/JPhysics/JConstants.hh b/src/jpp/JPhysics/JConstants.hh
similarity index 100%
rename from jpp/JPhysics/JConstants.hh
rename to src/jpp/JPhysics/JConstants.hh
diff --git a/src/jpp/JPhysics/JGeane.hh b/src/jpp/JPhysics/JGeane.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4c0c3a3d5c771d485d1b35b8e23bff28100e271b
--- /dev/null
+++ b/src/jpp/JPhysics/JGeane.hh
@@ -0,0 +1,334 @@
+#ifndef __JPHYSICS__JGEANE__
+#define __JPHYSICS__JGEANE__
+
+#include <cmath>
+#include <map>
+
+#include "JPhysics/JConstants.hh"
+
+
+/**
+ * \file
+ * Energy loss of muon.
+ * \author mdejong
+ */
+
+namespace JPHYSICS {}
+namespace JPP { using namespace JPHYSICS; }
+
+namespace JPHYSICS {
+
+  /**
+   * Equivalent muon track length per unit shower energy.
+   *
+   * \return        equivalent muon track length [m/Gev]
+   */
+  inline double geanc()
+  {
+    return 4.7;                   // dx/dE [m/GeV]
+  }
+
+
+  /**
+   * Interface for muon energy loss.
+   *
+   * This interface provides for the various function object operators.
+   */
+  class JGeane {
+  public:
+   /**
+     * Get energy loss constant.
+     *
+     * \return         Energy loss due to ionisation [GeV/m]
+     */
+    virtual double getA() const = 0;
+
+
+    /**
+     * Get energy loss constant.
+     *
+     * \return         Energy loss due to pair production and Bremstrahlung [m^-1]
+     */
+    virtual double getB() const = 0;
+
+
+    /**
+     * Get energy of muon after specified distance.
+     *
+     * \param  E       Energy of muon [GeV]
+     * \param  dx      distance traveled [m]
+     * \return         Energy of muon [GeV]
+     */
+    virtual double getE(const double E, const double dx) const = 0;
+
+
+    /**
+     * Get distance traveled by muon.
+     *
+     * \param  E0      Energy of muon at start [GeV]
+     * \param  E1      Energy of muon at end   [GeV]
+     * \return         distance traveled [m]
+     */
+    virtual double getX(const double E0, 
+			const double E1) const = 0;
+
+
+    /**
+     * Energy of muon after specified distance.
+     *
+     * \param  E       Energy of muon [GeV]
+     * \param  dx      distance traveled [m]
+     * \return         Energy of muon [GeV]
+     */
+    double operator()(const double E, const double dx) const
+    {
+      return this->getE(E, dx);
+    }
+    
+
+    /**
+     * Range of muon.
+     *
+     * \param  E       Energy of muon [GeV]
+     * \return         range [m]
+     */
+    double operator()(const double E) const
+    {
+      return this->getX(E, 0.0);
+    }
+
+
+    /**
+     * Equivalent unit track length per unit shower energy and per unit track length.
+     *
+     * \return        equivalent unit track length [Gev^-1]
+     */
+    double operator()() const
+    {
+      return this->getB() * geanc();
+    }
+  };
+
+
+  /**
+   * Function object for the energy loss of the muon.\n
+   * The energy loss can be formulated as:
+   *
+   *     \f[ -\frac{dE}{dx}  =  a + bE\f]
+   *
+   * N.B:
+   * \f$a\f$ and \f$b\f$ are assumed constant (internal units m and GeV, respectively).
+   */
+  class JGeane_t :
+    public JGeane
+  {
+  public:
+    /**
+     * constructor
+     * \param  __a     Energy loss due to ionisation [GeV/m]
+     * \param  __b     Energy loss due to pair production and Bremstrahlung [m^-1]
+     */ 
+    JGeane_t(const double __a,
+	     const double __b) :
+      a(__a),
+      b(__b)
+    {}
+    
+
+    /**
+     * Get energy loss constant.
+     *
+     * \return         Energy loss due to ionisation [GeV/m]
+     */
+    virtual double getA() const override 
+    {
+      return a;
+    }
+
+
+    /**
+     * Get energy loss constant.
+     *
+     * \return         Energy loss due to pair production and Bremstrahlung [m^-1]
+     */
+    virtual double getB() const override 
+    {
+      return b;
+    }
+
+
+    /**
+     * Get energy of muon after specified distance.
+     *
+     * \param  E       Energy of muon [GeV]
+     * \param  dx      distance traveled [m]
+     * \return         Energy of muon [GeV]
+     */
+    virtual double getE(const double E, const double dx) const override 
+    {
+      const double y = (a/b + E) * exp(-b*dx)  -  a/b;
+
+      if (y > 0.0)
+	return y;
+      else
+	return 0.0;
+    }
+
+
+    /**
+     * Get distance traveled by muon.
+     *
+     * \param  E0      Energy of muon at start [GeV]
+     * \param  E1      Energy of muon at end   [GeV]
+     * \return         distance traveled [m]
+     */
+    virtual double getX(const double E0, 
+			const double E1) const override
+    {
+      return -log((a + b*E1) / (a+b*E0)) / b;
+    }
+
+  protected:
+    const double a;
+    const double b;
+  };
+
+
+  /**
+   * Function object for energy dependent energy loss of the muon.
+   *
+   * Approximate values of energy loss parameters taken from reference:
+   * Proceedings of ICRC 2001, "Precise parametrizations of muon energy losses in water",
+   * S. Klimushin, E. Bugaev and I. Sokalski.
+   */
+  class JGeaneWater : 
+    public JGeane,
+    protected std::map<double, JGeane_t>
+  {
+  public:
+    /**
+     * Default constructor.
+     */
+    JGeaneWater()
+    {
+      using namespace std;
+
+      this->insert(make_pair( 0.0e0,  JGeane_t( 2.30e-1 * DENSITY_SEA_WATER, 15.50e-4 * DENSITY_SEA_WATER)));
+      this->insert(make_pair(30.0e0,  JGeane_t( 2.67e-1 * DENSITY_SEA_WATER,  3.40e-4 * DENSITY_SEA_WATER)));
+      this->insert(make_pair(35.3e3,  JGeane_t(-6.50e-1 * DENSITY_SEA_WATER,  3.66e-4 * DENSITY_SEA_WATER)));
+    }
+    
+
+    /**
+     * Get energy loss constant.
+     *
+     * N.B. The return value corresponds to the low-energy regime.
+     *
+     * \return         Energy loss due to ionisation [GeV/m]
+     */
+    virtual double getA() const override 
+    {
+      return 2.30e-1 * DENSITY_SEA_WATER; //This is the value for low energy (<30 GeV), the value used for step(ds) and getRange())
+    }
+
+
+    /**
+     * Get energy loss constant.
+     *
+     * N.B. The return value corresponds to the medium-energy regime.
+     *
+     * \return         Energy loss due to pair production and Bremstrahlung [m^-1]
+     */
+    virtual double getB() const override 
+    {
+      return 3.40e-4 * DENSITY_SEA_WATER;
+    }
+
+
+    /**
+     * Get energy of muon after specified distance.
+     *
+     * \param  E       Energy of muon [GeV]
+     * \param  dx      distance traveled [m]
+     * \return         Energy of muon [GeV]
+     */
+    virtual double getE(const double E, const double dx) const override 
+    {
+      double E1 = E;
+      double x1 = dx;
+
+      if (E1 > MASS_MUON) {
+
+	const_iterator p = this->lower_bound(E1);
+
+	do {
+
+	  --p;
+
+	  const double x2 = p->second.getX(E1, p->first);
+
+	  if (x2 > x1) {
+	    return p->second.getE(E1, x1);
+	  }
+
+	  E1  = p->first;
+	  x1 -= x2;
+
+	} while (p != this->begin());
+      }
+
+      return E1;
+    }
+
+
+    /**
+     * Get distance traveled by muon.
+     *
+     * \param  E0      Energy of muon at start [GeV]
+     * \param  E1      Energy of muon at end   [GeV]
+     * \return         distance traveled [m]
+     */
+    virtual double getX(const double E0, 
+			const double E1) const override
+    {
+      double E  = E0;
+      double dx = 0.0;
+
+      if (E > MASS_MUON) {
+
+	const_iterator p = this->lower_bound(E);
+
+	do {
+
+	  --p;
+
+	  if (E1 > p->first) {
+	    return dx + p->second.getX(E, E1);
+	  }
+
+	  dx += p->second.getX(E, p->first);
+	  E   = p->first;
+
+	} while (p != this->begin());
+      }
+
+      return dx;
+    }
+  };
+
+
+  /**
+   * Function object for energy loss of muon in sea water.
+   */
+  static const JGeaneWater gWater;
+
+
+  /**
+   * Function object for energy loss of muon in rock.
+   */
+  static const JGeane_t    gRock(2.67e-1 * 0.9 * DENSITY_ROCK, 
+				 3.40e-4 * 1.2 * DENSITY_ROCK);
+
+}
+
+#endif
diff --git a/jpp/JPhysics/JGeant_t.hh b/src/jpp/JPhysics/JGeant_t.hh
similarity index 100%
rename from jpp/JPhysics/JGeant_t.hh
rename to src/jpp/JPhysics/JGeant_t.hh
diff --git a/jpp/JPhysics/JGeanz.hh b/src/jpp/JPhysics/JGeanz.hh
similarity index 100%
rename from jpp/JPhysics/JGeanz.hh
rename to src/jpp/JPhysics/JGeanz.hh
diff --git a/jpp/JPhysics/JNPETable.hh b/src/jpp/JPhysics/JNPETable.hh
similarity index 100%
rename from jpp/JPhysics/JNPETable.hh
rename to src/jpp/JPhysics/JNPETable.hh
diff --git a/jpp/JPhysics/JNPE_t.hh b/src/jpp/JPhysics/JNPE_t.hh
similarity index 100%
rename from jpp/JPhysics/JNPE_t.hh
rename to src/jpp/JPhysics/JNPE_t.hh
diff --git a/jpp/JPhysics/JPDFTable.hh b/src/jpp/JPhysics/JPDFTable.hh
similarity index 100%
rename from jpp/JPhysics/JPDFTable.hh
rename to src/jpp/JPhysics/JPDFTable.hh
diff --git a/jpp/JPhysics/JPDFToolkit.hh b/src/jpp/JPhysics/JPDFToolkit.hh
similarity index 100%
rename from jpp/JPhysics/JPDFToolkit.hh
rename to src/jpp/JPhysics/JPDFToolkit.hh
diff --git a/jpp/JPhysics/JPDFTransformer.hh b/src/jpp/JPhysics/JPDFTransformer.hh
similarity index 100%
rename from jpp/JPhysics/JPDFTransformer.hh
rename to src/jpp/JPhysics/JPDFTransformer.hh
diff --git a/jpp/JPhysics/JPDFTypes.hh b/src/jpp/JPhysics/JPDFTypes.hh
similarity index 100%
rename from jpp/JPhysics/JPDFTypes.hh
rename to src/jpp/JPhysics/JPDFTypes.hh
diff --git a/jpp/JPhysics/JPDF_t.hh b/src/jpp/JPhysics/JPDF_t.hh
similarity index 100%
rename from jpp/JPhysics/JPDF_t.hh
rename to src/jpp/JPhysics/JPDF_t.hh
diff --git a/jpp/JTools/JAbstractCollection.hh b/src/jpp/JTools/JAbstractCollection.hh
similarity index 100%
rename from jpp/JTools/JAbstractCollection.hh
rename to src/jpp/JTools/JAbstractCollection.hh
diff --git a/jpp/JTools/JAbstractMultiMap.hh b/src/jpp/JTools/JAbstractMultiMap.hh
similarity index 100%
rename from jpp/JTools/JAbstractMultiMap.hh
rename to src/jpp/JTools/JAbstractMultiMap.hh
diff --git a/jpp/JTools/JArray.hh b/src/jpp/JTools/JArray.hh
similarity index 100%
rename from jpp/JTools/JArray.hh
rename to src/jpp/JTools/JArray.hh
diff --git a/jpp/JTools/JAssembler.hh b/src/jpp/JTools/JAssembler.hh
similarity index 100%
rename from jpp/JTools/JAssembler.hh
rename to src/jpp/JTools/JAssembler.hh
diff --git a/jpp/JTools/JCollection.hh b/src/jpp/JTools/JCollection.hh
similarity index 100%
rename from jpp/JTools/JCollection.hh
rename to src/jpp/JTools/JCollection.hh
diff --git a/jpp/JTools/JConstantFunction1D.hh b/src/jpp/JTools/JConstantFunction1D.hh
similarity index 100%
rename from jpp/JTools/JConstantFunction1D.hh
rename to src/jpp/JTools/JConstantFunction1D.hh
diff --git a/jpp/JTools/JConstants.hh b/src/jpp/JTools/JConstants.hh
similarity index 100%
rename from jpp/JTools/JConstants.hh
rename to src/jpp/JTools/JConstants.hh
diff --git a/jpp/JTools/JDistance.hh b/src/jpp/JTools/JDistance.hh
similarity index 100%
rename from jpp/JTools/JDistance.hh
rename to src/jpp/JTools/JDistance.hh
diff --git a/jpp/JTools/JElement.hh b/src/jpp/JTools/JElement.hh
similarity index 100%
rename from jpp/JTools/JElement.hh
rename to src/jpp/JTools/JElement.hh
diff --git a/jpp/JTools/JFunction1D_t.hh b/src/jpp/JTools/JFunction1D_t.hh
similarity index 100%
rename from jpp/JTools/JFunction1D_t.hh
rename to src/jpp/JTools/JFunction1D_t.hh
diff --git a/jpp/JTools/JFunctional.hh b/src/jpp/JTools/JFunctional.hh
similarity index 100%
rename from jpp/JTools/JFunctional.hh
rename to src/jpp/JTools/JFunctional.hh
diff --git a/jpp/JTools/JFunctionalMap.hh b/src/jpp/JTools/JFunctionalMap.hh
similarity index 100%
rename from jpp/JTools/JFunctionalMap.hh
rename to src/jpp/JTools/JFunctionalMap.hh
diff --git a/jpp/JTools/JFunctionalMap_t.hh b/src/jpp/JTools/JFunctionalMap_t.hh
similarity index 100%
rename from jpp/JTools/JFunctionalMap_t.hh
rename to src/jpp/JTools/JFunctionalMap_t.hh
diff --git a/jpp/JTools/JGarbageCollection.hh b/src/jpp/JTools/JGarbageCollection.hh
similarity index 100%
rename from jpp/JTools/JGarbageCollection.hh
rename to src/jpp/JTools/JGarbageCollection.hh
diff --git a/jpp/JTools/JGrid.hh b/src/jpp/JTools/JGrid.hh
similarity index 100%
rename from jpp/JTools/JGrid.hh
rename to src/jpp/JTools/JGrid.hh
diff --git a/jpp/JTools/JGridCollection.hh b/src/jpp/JTools/JGridCollection.hh
similarity index 100%
rename from jpp/JTools/JGridCollection.hh
rename to src/jpp/JTools/JGridCollection.hh
diff --git a/jpp/JTools/JGridMap.hh b/src/jpp/JTools/JGridMap.hh
similarity index 100%
rename from jpp/JTools/JGridMap.hh
rename to src/jpp/JTools/JGridMap.hh
diff --git a/jpp/JTools/JHermiteSpline.hh b/src/jpp/JTools/JHermiteSpline.hh
similarity index 100%
rename from jpp/JTools/JHermiteSpline.hh
rename to src/jpp/JTools/JHermiteSpline.hh
diff --git a/jpp/JTools/JHistogram.hh b/src/jpp/JTools/JHistogram.hh
similarity index 100%
rename from jpp/JTools/JHistogram.hh
rename to src/jpp/JTools/JHistogram.hh
diff --git a/jpp/JTools/JHistogramMap.hh b/src/jpp/JTools/JHistogramMap.hh
similarity index 100%
rename from jpp/JTools/JHistogramMap.hh
rename to src/jpp/JTools/JHistogramMap.hh
diff --git a/jpp/JTools/JMap.hh b/src/jpp/JTools/JMap.hh
similarity index 100%
rename from jpp/JTools/JMap.hh
rename to src/jpp/JTools/JMap.hh
diff --git a/jpp/JTools/JMapCollection.hh b/src/jpp/JTools/JMapCollection.hh
similarity index 100%
rename from jpp/JTools/JMapCollection.hh
rename to src/jpp/JTools/JMapCollection.hh
diff --git a/jpp/JTools/JMapList.hh b/src/jpp/JTools/JMapList.hh
similarity index 100%
rename from jpp/JTools/JMapList.hh
rename to src/jpp/JTools/JMapList.hh
diff --git a/jpp/JTools/JMappableCollection.hh b/src/jpp/JTools/JMappableCollection.hh
similarity index 100%
rename from jpp/JTools/JMappableCollection.hh
rename to src/jpp/JTools/JMappableCollection.hh
diff --git a/jpp/JTools/JMultiFunction.hh b/src/jpp/JTools/JMultiFunction.hh
similarity index 100%
rename from jpp/JTools/JMultiFunction.hh
rename to src/jpp/JTools/JMultiFunction.hh
diff --git a/jpp/JTools/JMultiHistogram.hh b/src/jpp/JTools/JMultiHistogram.hh
similarity index 100%
rename from jpp/JTools/JMultiHistogram.hh
rename to src/jpp/JTools/JMultiHistogram.hh
diff --git a/jpp/JTools/JMultiKey.hh b/src/jpp/JTools/JMultiKey.hh
similarity index 100%
rename from jpp/JTools/JMultiKey.hh
rename to src/jpp/JTools/JMultiKey.hh
diff --git a/jpp/JTools/JMultiMap.hh b/src/jpp/JTools/JMultiMap.hh
similarity index 100%
rename from jpp/JTools/JMultiMap.hh
rename to src/jpp/JTools/JMultiMap.hh
diff --git a/jpp/JTools/JMultiMapTransformer.hh b/src/jpp/JTools/JMultiMapTransformer.hh
similarity index 100%
rename from jpp/JTools/JMultiMapTransformer.hh
rename to src/jpp/JTools/JMultiMapTransformer.hh
diff --git a/jpp/JTools/JMultiPair.hh b/src/jpp/JTools/JMultiPair.hh
similarity index 100%
rename from jpp/JTools/JMultiPair.hh
rename to src/jpp/JTools/JMultiPair.hh
diff --git a/jpp/JTools/JPair.hh b/src/jpp/JTools/JPair.hh
similarity index 100%
rename from jpp/JTools/JPair.hh
rename to src/jpp/JTools/JPair.hh
diff --git a/jpp/JTools/JPolint.hh b/src/jpp/JTools/JPolint.hh
similarity index 100%
rename from jpp/JTools/JPolint.hh
rename to src/jpp/JTools/JPolint.hh
diff --git a/jpp/JTools/JQuadrature.hh b/src/jpp/JTools/JQuadrature.hh
similarity index 100%
rename from jpp/JTools/JQuadrature.hh
rename to src/jpp/JTools/JQuadrature.hh
diff --git a/jpp/JTools/JQuantiles.hh b/src/jpp/JTools/JQuantiles.hh
similarity index 100%
rename from jpp/JTools/JQuantiles.hh
rename to src/jpp/JTools/JQuantiles.hh
diff --git a/jpp/JTools/JRange.hh b/src/jpp/JTools/JRange.hh
similarity index 100%
rename from jpp/JTools/JRange.hh
rename to src/jpp/JTools/JRange.hh
diff --git a/jpp/JTools/JResult.hh b/src/jpp/JTools/JResult.hh
similarity index 100%
rename from jpp/JTools/JResult.hh
rename to src/jpp/JTools/JResult.hh
diff --git a/jpp/JTools/JResultTransformer.hh b/src/jpp/JTools/JResultTransformer.hh
similarity index 100%
rename from jpp/JTools/JResultTransformer.hh
rename to src/jpp/JTools/JResultTransformer.hh
diff --git a/jpp/JTools/JSet.hh b/src/jpp/JTools/JSet.hh
similarity index 100%
rename from jpp/JTools/JSet.hh
rename to src/jpp/JTools/JSet.hh
diff --git a/jpp/JTools/JSpline.hh b/src/jpp/JTools/JSpline.hh
similarity index 100%
rename from jpp/JTools/JSpline.hh
rename to src/jpp/JTools/JSpline.hh
diff --git a/jpp/JTools/JToolsToolkit.hh b/src/jpp/JTools/JToolsToolkit.hh
similarity index 100%
rename from jpp/JTools/JToolsToolkit.hh
rename to src/jpp/JTools/JToolsToolkit.hh
diff --git a/jpp/JTools/JTransformableMultiFunction.hh b/src/jpp/JTools/JTransformableMultiFunction.hh
similarity index 100%
rename from jpp/JTools/JTransformableMultiFunction.hh
rename to src/jpp/JTools/JTransformableMultiFunction.hh
diff --git a/jpp/JTools/JTransformableMultiHistogram.hh b/src/jpp/JTools/JTransformableMultiHistogram.hh
similarity index 100%
rename from jpp/JTools/JTransformableMultiHistogram.hh
rename to src/jpp/JTools/JTransformableMultiHistogram.hh
diff --git a/jpp/JTools/JTransformer.hh b/src/jpp/JTools/JTransformer.hh
similarity index 100%
rename from jpp/JTools/JTransformer.hh
rename to src/jpp/JTools/JTransformer.hh
diff --git a/jpp/Jeep/JeepToolkit.hh b/src/jpp/Jeep/JeepToolkit.hh
similarity index 100%
rename from jpp/Jeep/JeepToolkit.hh
rename to src/jpp/Jeep/JeepToolkit.hh
diff --git a/jpp/README.md b/src/jpp/README.md
similarity index 100%
rename from jpp/README.md
rename to src/jpp/README.md
diff --git a/jppy/.gitkeep b/src/jppy/.gitkeep
similarity index 100%
rename from jppy/.gitkeep
rename to src/jppy/.gitkeep
diff --git a/src/jppy/__init__.py b/src/jppy/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..23f9374872cbfc8db797a8f2f0e4d68024d7bced
--- /dev/null
+++ b/src/jppy/__init__.py
@@ -0,0 +1,12 @@
+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
diff --git a/src/jppy/constants.cc b/src/jppy/constants.cc
new file mode 100644
index 0000000000000000000000000000000000000000..6f0b808c1eb5ed5f5679c00d407ca2280a2ac941
--- /dev/null
+++ b/src/jppy/constants.cc
@@ -0,0 +1,18 @@
+#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);
+}
+
diff --git a/src/jppy/geane.cc b/src/jppy/geane.cc
new file mode 100644
index 0000000000000000000000000000000000000000..626eb08bac2fbeef276a1bc39fa3aef074300e8f
--- /dev/null
+++ b/src/jppy/geane.cc
@@ -0,0 +1,35 @@
+#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]")
+	 );
+}
diff --git a/src/npe.cc b/src/jppy/npe.cc
similarity index 100%
rename from src/npe.cc
rename to src/jppy/npe.cc
diff --git a/src/pdf.cc b/src/jppy/pdf.cc
similarity index 100%
rename from src/pdf.cc
rename to src/jppy/pdf.cc
diff --git a/src/jppy/pdf_evaluator.py b/src/jppy/pdf_evaluator.py
new file mode 100644
index 0000000000000000000000000000000000000000..3ad447fb7b3813ce7b78900c713779d22a768c1f
--- /dev/null
+++ b/src/jppy/pdf_evaluator.py
@@ -0,0 +1,164 @@
+#!/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)
diff --git a/tests/test_geane.py b/tests/test_geane.py
new file mode 100644
index 0000000000000000000000000000000000000000..484736af0a2726928ad239e2d199783f81560f38
--- /dev/null
+++ b/tests/test_geane.py
@@ -0,0 +1,11 @@
+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)
diff --git a/tests/test_pdf.py b/tests/test_pdf.py
index 60f987a7cfe88144fdfd86d0c60717322272a325..6936c8cfd6898e516efdf1f0e706920250fa4a94 100644
--- a/tests/test_pdf.py
+++ b/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)
diff --git a/tests/test_pdf_evaluator.py b/tests/test_pdf_evaluator.py
new file mode 100644
index 0000000000000000000000000000000000000000..42bd0619d87cb1468c07ad0e912529781f8aece2
--- /dev/null
+++ b/tests/test_pdf_evaluator.py
@@ -0,0 +1,27 @@
+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)
+
+
+