JLine3ZRegressor.hh

Go to the documentation of this file.

#ifndef __JFIT__JLINE3ZREGRESSOR__
#define __JFIT__JLINE3ZREGRESSOR__
 
#include <string>
#include <array>
#include <memory>
#include <limits>
 
#include "JPhysics/JConstants.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JNPETable.hh"
#include "JPhysics/JDeltaRays.hh"
#include "JPhysics/JGeane.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JRange.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JMath/JZero.hh"
#include "JFit/JTimeRange.hh"
#include "JFit/JLine3Z.hh"
#include "JFit/JPMTW0.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JGandalf.hh"
#include "JFit/JMEstimator.hh"
#include "JFit/JRegressor.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 * Data regression method for JFIT::JLine3Z.
 * \author mdejong
 */
namespace JFIT {}
namespace JPP { using namespace JFIT; }
 
namespace JFIT {
 
  using JPHYSICS::JPDFType_t;
  using JPHYSICS::DIRECT_LIGHT_FROM_MUON;
  using JPHYSICS::SCATTERED_LIGHT_FROM_MUON;
  using JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS;
  using JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS;
  using JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWERS;
  using JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWERS;
  
 
  /**
   * Regressor function object for JLine3Z fit using JSimplex minimiser.
   */
  template<>
  struct JRegressor<JLine3Z, JSimplex> :
    public JAbstractRegressor<JLine3Z, JSimplex>
  {
    using JAbstractRegressor<JLine3Z, JSimplex>::operator();
 
    /**
     * Constructor.
     *
     * \param  sigma              time resolution of hit [ns]
     */
    JRegressor(double sigma) :
      estimator(new JMEstimatorNormal())
    {
      this->sigma = sigma;
    }
 
 
    /**
     * Fit function.
     * This method is used to determine the chi2 of given hit with respect to trajectory of muon.
     *
     * The template argument <tt>JHit_t</tt> refers to a data structure which should have the following member methods:
     *   - double getX();       // [m]
     *   - double getY();       // [m]
     *   - double getZ();       // [m]
     *   - double getT();       // [ns]
     *
     * \param  track              track
     * \param  hit                hit
     * \return                    chi2
     */
    template<class JHit_t>
    double operator()(const JLine3Z& track, const JHit_t& hit) const
    {
      using namespace JPP;
 
      JPosition3D D(hit.getX(), hit.getY(), hit.getZ());
 
      D.sub(track.getPosition());
 
      const double z  = D.getDot(track.getDirection());
      const double R  = sqrt(D.getLengthSquared() - z*z);
 
      const double t1 = track.getT() + (z + R * getKappaC()) * getInverseSpeedOfLight();
 
      const double u  = (t1 - hit.getT()) / sigma;
 
      return estimator->getRho(u) * hit.getW();
    }
 
    std::shared_ptr<JMEstimator>       estimator;  //!< M-Estimator function
    double                             sigma;      //!< Time resolution [ns]
  };
 
 
  /**
   * Template specialisation for storage of PDF tables.
   */
  template<>
  struct JRegressorStorage<JLine3Z, JGandalf>
  {
    typedef JTOOLS::JSplineFunction1S_t                                   JFunction1D_t;
    typedef JTOOLS::JMAPLIST<JTOOLS::JPolint1FunctionalMap,
                             JTOOLS::JPolint0FunctionalGridMap,
                             JTOOLS::JPolint0FunctionalGridMap>::maplist  JPDFMaplist_t;
    typedef JPHYSICS::JPDFTable<JFunction1D_t, JPDFMaplist_t>             JPDF_t;    //!< time dependent PDF
 
    typedef JTOOLS::JMAPLIST<JTOOLS::JPolint1FunctionalMapH,
                             JTOOLS::JPolint1FunctionalGridMap,
                             JTOOLS::JPolint1FunctionalGridMap>::maplist  JNPEMaplist_t;
    typedef JPHYSICS::JNPETable<double, double, JNPEMaplist_t>            JNPE_t;    //!< time integrated PDF
 
    typedef JPDF_t::transformer_type                                      transformer_type;
 
    static const int           NUMBER_OF_PDFS  =  6;                                 //!< Number of PDFs
 
    typedef std::array<JPDF_t, NUMBER_OF_PDFS>                            JPDFs_t;   //!< PDFs
    typedef std::array<JNPE_t, NUMBER_OF_PDFS>                            JNPEs_t;   //!< NPEs
    
 
    /**
     * Default constructor.
     */
    JRegressorStorage()
    {}
 
 
    /**
     * Constructor.
     *
     * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
     * will be replaced by the PDF types listed in JRegressorStorage<JLine3Z, JGandalf>::pdf_t.
     *
     * The <tt>TTS_ns</tt> corresponds to the additional time smearing applied to the PDFs.
     *
     * \param  fileDescriptor     PDF file descriptor
     * \param  T_ns               time range [ns]
     * \param  TTS_ns             TTS [ns]
     * \param  numberOfPoints     number of points for Gauss-Hermite integration of TTS
     * \param  epsilon            precision        for Gauss-Hermite integration of TTS
     */
    JRegressorStorage(const std::string& fileDescriptor,
                      const JTimeRange&  T_ns, 
                      const double       TTS_ns,
                      const int          numberOfPoints = 25,
                      const double       epsilon        = 1.0e-10) :
      T_ns(T_ns)
    {
      using namespace std;
      using namespace JPP;
 
      const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
 
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
 
        const string file_name = getFilename(fileDescriptor, pdf_t[i]);
 
        _pdf[i].load(file_name.c_str());
 
        _pdf[i].setExceptionHandler(supervisor);
      }
      
      // Add PDFs
 
      for (int i = 1; i < NUMBER_OF_PDFS; i += 2) {
 
        _pdf[ i ].add(_pdf[i-1]);
        _pdf[ i ].compress(numeric_limits<double>::max(), T_ns);
 
        JPDF_t buffer;
 
        _pdf[i-1].swap(buffer);
 
        if (TTS_ns > 0.0) {
          _pdf[i].blur(TTS_ns, numberOfPoints, epsilon);
        }
 
        _npe[ i ] = JNPE_t(_pdf[i], T_ns);
      }
    }
    
 
    /**
     * Get PDFs.
     *
     * \return                    PDFs
     */
    const JPDFs_t& getPDF() const
    {
      return _pdf;
    }
    
 
    /**
     * Get NPEs.
     *
     * \return                     NPEs
     */
    const JNPEs_t& getNPE() const
    {
      return _npe;
    }
 
 
    /**
     * Transform PDFs and NPEs.
     *
     * \param  transformer        transformer
     */
    void transform(const transformer_type& transformer)
    {
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
        _pdf[i].transform(transformer);
        _npe[i].transform(transformer);
      }
    }
 
 
    /**
     * PDF types.
     */
    static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
 
    JTimeRange T_ns;                    //!< Time window with respect to Cherenkov hypothesis [ns]
 
  private:
    JPDFs_t _pdf;                       //!< PDFs
    JNPEs_t _npe;                       //!< NPEs
  };
 
 
  /**
   * PDF types.
   */
  const JPDFType_t JRegressorStorage<JLine3Z, JGandalf>::pdf_t[] = {
    DIRECT_LIGHT_FROM_MUON,
    SCATTERED_LIGHT_FROM_MUON,
    DIRECT_LIGHT_FROM_DELTARAYS,
    SCATTERED_LIGHT_FROM_DELTARAYS,
    DIRECT_LIGHT_FROM_EMSHOWERS,
    SCATTERED_LIGHT_FROM_EMSHOWERS
  };
 
 
  /**
   * Regressor function object for JLine3Z fit using JGandalf minimiser.
   */
  template<>
  struct JRegressor<JLine3Z, JGandalf> :
    public JAbstractRegressor<JLine3Z, JGandalf>,
    public JRegressorStorage <JLine3Z, JGandalf>
  {
    using JAbstractRegressor<JLine3Z, JGandalf>::operator();
 
    typedef JRegressorStorage<JLine3Z, JGandalf>                          storage_type;
 
    /**
     * Default constructor
     */
    JRegressor() :
      storage_type(),
      pdf(getPDF()),
      npe(getNPE()),
      estimator()
    {}
 
 
    /**
     * Constructor.
     *
     * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
     * will be replaced by the PDF types listed in JRegressorStorage<JLine3Z, JGandalf>::pdf_t.
     *
     * The <tt>TTS_ns</tt> corresponds to the additional Gaussian time smearing applied to the PDFs.
     *
     * \param  fileDescriptor     PDF file descriptor
     * \param  T_ns               time range [ns]
     * \param  TTS_ns             TTS [ns]
     * \param  numberOfPoints     number of points for Gauss-Hermite integration of TTS
     * \param  epsilon            precision        for Gauss-Hermite integration of TTS
     */
    JRegressor(const std::string& fileDescriptor,
               const JTimeRange&  T_ns, 
               const double       TTS_ns,
               const int          numberOfPoints = 25,
               const double       epsilon        = 1.0e-10) :
      storage_type(fileDescriptor, T_ns, TTS_ns, numberOfPoints, epsilon),
      pdf(getPDF()),
      npe(getNPE()),
      estimator(new JMEstimatorNormal())
    {}
 
 
    /**
     * Constructor.
     *
     * \param  storage            PDF storage
     */
    JRegressor(const storage_type& storage) :
      pdf(storage.getPDF()),
      npe(storage.getNPE()),
      estimator(new JMEstimatorNormal())
    {
      T_ns = storage.T_ns;
    }
 
 
    /**
     * Fit function.
     * This method is used to determine the chi2 and gradient of given hit with respect to trajectory of muon.
     *
     * The template argument <tt>JHit_t</tt> refers to a data structure which should have the following member methods:
     *   - double getX();       // [m]
     *   - double getY();       // [m]
     *   - double getZ();       // [m]
     *   - double getDX();      // [u]
     *   - double getDY();      // [u]
     *   - double getDZ();      // [u]
     *   - double getT();       // [ns]
     *   - double getQE();      // relative quantum efficiency
     *   - double getR();       // rate [Hz]
     *
     * \param  track              track
     * \param  hit                hit
     * \return                    chi2 and gradient
     */
    template<class JHit_t>
    result_type operator()(const JLine3Z& track, const JHit_t& hit) const
    {
      using namespace JPP;
      
      JPosition3D  D(hit.getX(),  hit.getY(),  hit.getZ());
      JDirection3D U(hit.getDX(), hit.getDY(), hit.getDZ());
 
      D.sub(track.getPosition());
 
      const double z  = D.getDot(track.getDirection());
      const double x  = D.getX()  -  z * track.getDX();
      const double y  = D.getY()  -  z * track.getDY();
      const double R2 = D.getLengthSquared() - z*z;
      const double R  = sqrt(R2);
 
      const double t1 = track.getT() + (z + R * getTanThetaC()) * getInverseSpeedOfLight();
 
      U.rotate(JRotation3Z(-atan2(y,x)));                  // rotate PMT axis to x-z plane
 
      const double theta = U.getTheta();
      const double phi   = fabs(U.getPhi());
 
      const double E  = gWater.getE(E_GeV, z);
      const double dt = T_ns.constrain(hit.getT()  -  t1);
 
      JPDF_t::result_type H0 = getH0(hit.getR(), dt);
      JPDF_t::result_type H1 = getH1(E, R, theta, phi, dt);
 
      H1 *= hit.getQE();
 
      if (H1.V >= Vmax_npe) {
        H1 *= Vmax_npe / H1.V;
      }
 
      H1 += H0;                                                                                      // signal + background
 
      result_type result;
 
      result.chi2     = H1.getChi2() - H0.getChi2();                                                 // Likelihood ratio
 
      const double wc = 1.0  -  getTanThetaC() * z / R;                                              // d(ct1)/d(z)
 
      result.gradient = JLine3Z(JLine1Z(JPosition3D(-getTanThetaC() * D.getX() / R,                  // d(ct1)/d(x)
                                                    -getTanThetaC() * D.getY() / R,                  // d(ct1)/d(y)
                                                    0.0),
                                        getSpeedOfLight()),                                          // d(ct1)/d(t)
                                JVersor3Z(wc * (D.getX() - D.getZ()*track.getDX()/track.getDZ()),    // d(ct1)/d(dx)
                                          wc * (D.getY() - D.getZ()*track.getDY()/track.getDZ())));  // d(ct1)/d(dy)
      
      result.gradient.mul(getInverseSpeedOfLight() * (H1.getDerivativeOfChi2() -
                                                      H0.getDerivativeOfChi2()));                    // x d(chi2)/d(ct1)
 
      return result;
    }
 
 
    /**
     * Fit function.
     * This method is used to determine the chi2 and gradient of given PMT with respect to trajectory of muon.
     *
     * \param  track              track
     * \param  pmt                pmt
     * \return                    chi2 and gradient
     */
    result_type operator()(const JLine3Z& track, const JPMTW0& pmt) const
    {
      using namespace JGEOMETRY3D;
      using namespace JPHYSICS;
      
      JPosition3D  D(pmt.getPosition());
      JDirection3D U(pmt.getDirection());
 
      D.sub(track.getPosition());
 
      const double z  = D.getDot(track.getDirection());
      const double x  = D.getX()  -  z * track.getDX();
      const double y  = D.getY()  -  z * track.getDY();
      const double R2 = D.getLengthSquared() - z*z;
      const double R  = sqrt(R2);
 
      U.rotate(JRotation3Z(-atan2(y,x)));                  // rotate PMT axis to x-z plane
 
      const double theta = U.getTheta();
      const double phi   = fabs(U.getPhi());
 
      const double E  = gWater.getE(E_GeV, z);
 
      JNPE_t::result_type H0 = getH0(pmt.getR());
      JNPE_t::result_type H1 = getH1(E, R, theta, phi);
 
      H1 *= pmt.getQE();
 
      if (H1.f >= Vmax_npe) {
        H1 *= Vmax_npe / H1.f;
      }
 
      H1 += H0;
 
      const bool   hit = pmt.getN() != 0;
      const double u   = H1.getChi2(hit);
 
      result_type result;
 
      result.chi2     = estimator->getRho(u);
 
      result.gradient = JLine3Z(JLine1Z(JPosition3D(-D.getX() / R,   // d(R)/d(x)
                                                    -D.getY() / R,   // d(R)/d(y)
                                                    0.0),
                                        0.0),
                                JVersor3Z(-z * D.getX() / R,         // d(R)/d(dx)
                                          -z * D.getY() / R));       // d(R)/d(dy)
 
      result.gradient.mul(estimator->getPsi(u));
      result.gradient.mul(H1.getDerivativeOfChi2(hit));              // x d(chi2)/d(R)
 
      return result;
    }
 
 
    /**
     * Get background hypothesis value for time differentiated PDF.
     *
     * \param  R_Hz               rate [Hz]
     * \param  t1                 time [ns]
     * \return                    hypothesis value
     */
    JPDF_t::result_type getH0(const double R_Hz,
                              const double t1) const
    {
      return JPDF_t::result_type(R_Hz * 1e-9, t1, T_ns);
    }
 
 
    /**
     * Get signal hypothesis value for time differentiated PDF.
     *
     * \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
     */
    JPDF_t::result_type getH1(const double E,
                              const double R,
                              const double theta,
                              const double phi,
                              const double t1) const
    {
      using namespace std;
      using namespace JPP;
 
      JPDF_t::result_type h1 = zero;
 
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
 
        if (!pdf[i].empty() && R <= pdf[i].getXmax()) {
 
          JPDF_t::result_type y1 = pdf[i](max(R, pdf[i].getXmin()), theta, phi, t1);
          
          // safety measures
          
          if (y1.f <= 0.0) {
            y1.f  = 0.0;
            y1.fp = 0.0;
          }
          
          if (y1.v <= 0.0) {
            y1.v  = 0.0;
          }
          
          // energy dependence
          
          if        (is_deltarays(pdf_t[i])) {
            y1 *= JDeltaRays::getEnergyLossFromMuon(E);
          } else if (is_bremsstrahlung(pdf_t[i])) {
            y1 *= E;
          }
          
          h1 += y1;
        }
      }
 
      return h1;
    }
 
 
    /**
     * Get background hypothesis value for time integrated PDF.
     *
     * \param  R_Hz               rate [Hz]
     * \return                    hypothesis value
     */
    JNPE_t::result_type getH0(const double R_Hz) const
    {
      return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength(), 0.0);
    }
 
 
    /**
     * Get signal hypothesis value for time integrated PDF.
     *
     * \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]
     * \return                    hypothesis value
     */
    JNPE_t::result_type getH1(const double E,
                              const double R,
                              const double theta,
                              const double phi) const
    {
      using namespace std;
      using namespace JPP;
 
      JNPE_t::result_type h1 = JMATH::zero;
        
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
 
        if (!npe[i].empty() && R <= npe[i].getXmax()) {
 
          try {
 
            JNPE_t::result_type y1 = npe[i](max(R, npe[i].getXmin()), theta, phi);
          
            // safety measures
            
            if (y1.f > 0.0) {
            
              // energy dependence
            
              if        (is_deltarays(pdf_t[i])) {
                y1 *= JDeltaRays::getEnergyLossFromMuon(E);
              } else if (is_bremsstrahlung(pdf_t[i])) {
                y1 *= E;
              }
          
              h1 += y1;
            }
          }
          catch(JException& error) {
            ERROR(error << endl);
          }
        }
      }
 
      return h1;
    }
    
 
    /**
     * Get maximal road width of PDF.
     *
     * \return                    road width [m]
     */
    inline double getRmax() const
    {
      double xmax = 0.0;
 
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
        if (!pdf[i].empty() && pdf[i].getXmax() > xmax) {
          xmax = pdf[i].getXmax();
        }
      }
 
      return xmax;
    }
 
 
    static double       Vmax_npe;                  //!< Maximal integral of PDF [npe]
 
    const JPDFs_t&      pdf;                       //!< PDF
    const JNPEs_t&      npe;                       //!< PDF
 
    double              E_GeV = 0.0;               //!< Energy of muon at vertex [GeV]
 
    std::shared_ptr<JMEstimator> estimator;        //!< M-Estimator function
  };
 
 
  /**
   * Default values.
   */
  double     JRegressor<JLine3Z, JGandalf>::Vmax_npe = std::numeric_limits<double>::max();
}
 
#endif