JPDF_t.hh

Go to the documentation of this file.

#include "JLang/JException.hh"
#include "JTools/JCollection.hh"
#include "JTools/JMap.hh"
#include "JTools/JGridMap.hh"
#include "JTools/JMapList.hh"
#include "JTools/JSpline.hh"
#include "JTools/JPolint.hh"
#include "JTools/JElement.hh"
#include "JTools/JResult.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JMath/JZero.hh"
 
 
/**
 * \file
 *
 * Auxiliary data structure for muon PDF.
 * \author mdejong
 */
 
/**
 * Auxiliary data structure for muon PDF.
 */
struct JPDF {
 
  typedef JPP::JSplineFunction1D<JPP::JSplineElement2S<double, double>, 
                                 JPP::JCollection, 
                                 JPP::JResultPDF<double> >         JFunction1D_t;
 
  typedef JPP::JMAPLIST<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 <tt>TTS</tt> corresponds to the additional time smearing applied to the PDFs.
   *
   * \param  file_name          file name
   * \param  TTS                TTS [ns]
   * \param  numberOfPoints     number of points for Gauss-Hermite integration of TTS
   * \param  epsilon            precision        for Gauss-Hermite integration of TTS
   */
  JPDF(const std::string& file_name,
       const double       TTS,
       const int          numberOfPoints = 25,
       const double       epsilon        = 1.0e-10)
  {
    using namespace std;
    using namespace JPP;
 
    const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(zero));
 
    cout << "loading input from file " << file_name << "... " << flush;
 
    pdf.load(file_name.c_str());
 
    pdf.setExceptionHandler(supervisor);
 
    cout << "OK" << endl;
 
    type = getPDFType(file_name);
 
    if        (TTS > 0.0) {
 
      cout << "bluring PDFs... " << flush;
 
      pdf.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                  muon energy at minimum distance of approach [GeV]
   * \param  R                  minimum distance of approach [m]
   * \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 R,
                        const double theta,
                        const double phi,
                        const double t1) const
  {
    using namespace JPP;
 
    result_type h1 = pdf(R, theta, phi, t1);
 
    if        (is_bremsstrahlung(type)) {
      h1 *= E;
    } else if (is_deltarays(type)) {
      h1 *= getDeltaRaysFromMuon(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 pdf;    //!< PDF
  int    type;   //!< PDF type
};
 
 
/**
 * Auxiliary data structure for muon PDF.
 */
struct JMuonPDF_t {
 
  typedef JPP::JSplineFunction1D<JPP::JSplineElement2S<double, double>, 
                                 JPP::JCollection, 
                                 JPP::JResultPDF<double> >         JFunction1D_t;
 
  typedef JPP::JMAPLIST<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::WILDCARD.\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
   */
  JMuonPDF_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[] = { DIRECT_LIGHT_FROM_MUON,
                                 SCATTERED_LIGHT_FROM_MUON,
                                 DIRECT_LIGHT_FROM_EMSHOWERS,
                                 SCATTERED_LIGHT_FROM_EMSHOWERS,
                                 DIRECT_LIGHT_FROM_DELTARAYS,
                                 SCATTERED_LIGHT_FROM_DELTARAYS };
 
    const  int N = sizeof(pdf_t) / sizeof(pdf_t[0]);
    
    JPDF_t pdf[N];
 
    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;
 
      pdf[i].load(file_name.c_str());
 
      pdf[i].setExceptionHandler(supervisor);
 
      cout << "OK" << endl;
    }
 
    // Add PDFs
 
    cout << "adding PDFs... " << flush;
 
    pdfA = pdf[1];  pdfA.add(pdf[0]);
    pdfB = pdf[3];  pdfB.add(pdf[2]);
    pdfC = pdf[5];  pdfC.add(pdf[4]);
 
    cout << "OK" << endl;
 
    if        (TTS > 0.0) {
 
      cout << "bluring PDFs... " << flush;
 
      pdfA.blur(TTS, numberOfPoints, epsilon);
      pdfB.blur(TTS, numberOfPoints, epsilon);
      pdfC.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                  muon energy at minimum distance of approach [GeV]
   * \param  R                  minimum distance of approach [m]
   * \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 R,
                        const double theta,
                        const double phi,
                        const double t1) const
  {
    using namespace JPP;
 
    result_type h1 = (pdfA(R, theta, phi, t1)                               +
                      pdfB(R, theta, phi, t1) * E                           +
                      pdfC(R, theta, phi, t1) * getDeltaRaysFromMuon(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 minimum ionising particle
  JPDF_t pdfB;   //!< PDF for average energy losses
  JPDF_t pdfC;   //!< PDF for delta-rays
};
 
 
/**
 * 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::WILDCARD.\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 [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
};