JShowerBjorkenYRegressor.hh

Go to the documentation of this file.

#ifndef __JFIT__JSHOWERBJORKENYREGRESSOR__
#define __JFIT__JSHOWERBJORKENYREGRESSOR__
 
#include <memory>
 
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JNPETable.hh"
 
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JConstants.hh"
#include "JTools/JRange.hh"
#include "JTools/JResult.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JAbstractHistogram.hh"
 
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JVersor3D.hh"
#include "JGeometry3D/JDirection3D.hh"
 
#include "JMath/JZero.hh"
 
#include "JFit/JTimeRange.hh"
#include "JFit/JPMTW0.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JMEstimator.hh"
#include "JFit/JRegressor.hh"
#include "JFit/JShowerEH.hh"
#include "JFit/JFitToolkit.hh"
 
#include "Jeep/JMessage.hh"
 
/**
 * \file
 * Data regression method for JFIT::JShowerEH.
 * \author adomi
 */
 
namespace JFIT {
 
  using JPHYSICS::JPDFType_t;
  using JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER;
  using JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER;
  using JPHYSICS::DIRECT_LIGHT_FROM_BRIGHT_POINT;
  using JPHYSICS::SCATTERED_LIGHT_FROM_BRIGHT_POINT;
  using JTOOLS::get_value;
 
  /**
   * Regressor function object for JShowerEH fit using JSimplex minimiser.
   */
  template<>
  struct JRegressor<JShowerEH, JSimplex> :
    public JAbstractRegressor<JShowerEH, JSimplex>
  {
    using JAbstractRegressor<JShowerEH, JSimplex>::operator();
 
    typedef JTOOLS::JSplineFunction1S_t                                     JFunction1D_t;
    typedef JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                             JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                                              JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap,
                                                               JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap> > > >        JPDFMaplist_t;
    typedef JPHYSICS::JPDFTable<JFunction1D_t, JPDFMaplist_t>                JPDF_t;
    
    typedef JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                             JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                                              JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap,
                                                               JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap> > > >        JNPEMaplist_t;
    typedef JPHYSICS::JNPETable<double, double, JNPEMaplist_t>               JNPE_t;
 
 
    // Part for the Isotropic PDF
    typedef JTOOLS::JMAPLIST<JTOOLS::JPolint1FunctionalMap,
                             JTOOLS::JPolint1FunctionalGridMap>::maplist     JMapList_t2;
    typedef JPHYSICS::JPDFTable<JFunction1D_t, JMapList_t2>                  JPDF_t2;
 
    typedef JTOOLS::JMAPLIST<JTOOLS::JPolint1FunctionalMap,
                             JTOOLS::JPolint1FunctionalGridMap>::maplist     JNPEMapList_t2;
    typedef JPHYSICS::JNPETable<double, double, JNPEMapList_t2>              JNPE_t2;
 
    
    /**
     * Parameterized constructor
     *
     * The PDF file descriptor should contain the wild card character JEEP::FILENAME_WILDCARD which
     * will be replaced by the corresponding PDF types listed in JRegressor<JShower3Z, JGandalf>::pdf_t.
     *
     * \param  fileDescriptor     PDF file descriptor
     */
 
    JRegressor(const std::string& fileDescriptor):
      estimator(new JMEstimatorNull())
    {
      using namespace std;
      using namespace JPP;
 
      const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
      const JPDF_t2::JSupervisor supervisor2(new JPDF_t2::JDefaultResult(JMATH::zero));
      
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
        
        try {
            
          JPDF_t pdf;
          JPDF_t2 pdf2;
 
          const string file_name = getFilename(fileDescriptor, pdf_t[i]);
 
          NOTICE("loading PDF from file " << file_name << "... " << flush);
 
          if(i < 2){
 
            pdf.load(file_name.c_str());
 
            pdf.setExceptionHandler(supervisor);
          
            npe[ i ] = JNPE_t(pdf);
 
          } else {
            
            pdf2.load(file_name.c_str());
          
            NOTICE("OK" << endl);
 
            pdf2.setExceptionHandler(supervisor2);
          
            npe2[ i-2 ] = JNPE_t2(pdf2);
        
          }
        }
        catch(const JException& error) {
          FATAL(error.what() << endl);
        }
      }
 
      // Add PDFs
      for (int i = 1; i < (NUMBER_OF_PDFS-2); i += 2) {
 
        npe[ i ].add(npe[i-1]);
        
        JNPE_t buffer;
        
        npe[i-1].swap(buffer);
          
        npe2[ i ].add(npe2[i-1]);
          
        JNPE_t2 buffer2;
 
        npe2[i-1].swap(buffer2);        
          
      }   
    }
 
    /**
     * Fit function.
     * This method is used to determine the chi2 of given PMT with respect to shower hypothesis.
     *
     * \param  shower             shower
     * \param  pmt                pmt
     * \return                    chi2 
     */
    double operator()(const JShowerEH& shower, const JPMTW0& pmt) const
    {
      using namespace JPP;
      using namespace std;
      
      JPosition3D D(pmt.getPosition());
      JDirection3D U(pmt.getDirection());
 
      D.sub(shower.getPosition());
 
      double ct = U.getDot(D) / D.getLength();
      
      JVersor3D shower_dir(shower.getDirection());
 
      const double z  = D.getDot(shower_dir);
      const double x  = D.getX()  -  z * shower.getDX();
      const double y  = D.getY()  -  z * shower.getDY();
      const double cosDelta = z/D.getLength();  // Delta = angle between shower direction and PMT position
 
      U.rotate(JRotation3Z(-atan2(y,x)));                  // rotate PMT axis to x-z plane
 
      const double theta = U.getTheta();
      const double phi   = fabs(U.getPhi());
 
      double H0 = getH0(pmt.getR());  // background
      double H1 = getH1(D.getLength(), ct, cosDelta, theta, phi,
                        shower.getEem(), shower.getEh(), shower.getBy());  // signal
 
      if (H1 >= Vmax_npe) {
        H1 *= Vmax_npe / H1;
      }
      
      H1 += H0; // now H1 is signal + background
 
      const bool hit = pmt.getN() != 0;
      const double u = getChi2(H1, hit); // -log(lik)
      
      return estimator->getRho(u);
    }
 
    /**
     * Get background hypothesis value for time integrated PDF.
     *
     * \param  R_Hz               rate [Hz]
     * \return                    hypothesis value
     */
    double getH0(const double R_Hz) const
    {
      return get_value(JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength()));
    }
 
    /**
     * Get signal hypothesis value for time integrated PDF.
     *
     * \param  D                  PMT distance from shower [m]
     * \param  ct                 angle between shower direction and PMT position
     * \param  cosDelta           angle between shower direction and PMT position
     * \param  theta              PMT zenith  angle [deg]
     * \param  phi                PMT azimuth angle [deg]
     * \param  Eem                EM shower energy [GeV]
     * \param  Eh                 H shower energy [GeV]
     * \param  Y                  Bjorken Y        
     * \return                    hypothesis value
     */
    double getH1(const double D,
                 const double ct,
                 const double cosDelta,
                 const double theta,
                 const double phi,
                 const double Eem,
                 const double Eh,
                 const double Y) const
    {
      
      double h1 = 0;
      
      for (int i = 0; i != (NUMBER_OF_PDFS-1); ++i) {
 
        if (!npe[i].empty() && D <= npe[i].getXmax() && !npe2[i].empty() && D <= npe2[i].getXmax()) {
 
          try {
 
            JNPE_t::result_type P_em;
            JNPE_t2::result_type P_h;
 
            P_em = fabs(Eem) * npe[i](std::max(D, npe[i].getXmin()), cosDelta, theta, phi);
  
            P_h = fabs(Eh) * npe2[i](std::max(D, npe2[i].getXmin()), ct);
 
            double y1 = get_value(P_em) + get_value(P_h);
            
            if(y1 > 0.0){
              h1 += y1;
            }
            
          }
          catch(JLANG::JException& error) {
            ERROR(error << std::endl);
          }
        }
      }
 
      return h1;
    }
 
    static JTimeRange   T_ns;                      //!< Time window with respect to Cherenkov hypothesis [ns]
    static double       Vmax_npe;                  //!< Maximal integral of PDF [npe]
 
    static const int    NUMBER_OF_PDFS = 4;
 
    static const JPDFType_t  pdf_t[NUMBER_OF_PDFS];
 
    JNPE_t              npe[NUMBER_OF_PDFS-2];       //!< PDF
    JNPE_t2             npe2[NUMBER_OF_PDFS-2];      //!< PDF
 
    std::shared_ptr<JMEstimator>         estimator;  //!< M-Estimator function
  };  
 
  /**
   * PDF types.
   */
  const JPDFType_t  JRegressor<JShowerEH, JSimplex>::pdf_t[]  = { DIRECT_LIGHT_FROM_EMSHOWER, 
                                                                  SCATTERED_LIGHT_FROM_EMSHOWER,
                                                                  DIRECT_LIGHT_FROM_BRIGHT_POINT,
                                                                  SCATTERED_LIGHT_FROM_BRIGHT_POINT };
  
  /**
   * Default values.
   */
  JTimeRange JRegressor<JShowerEH, JSimplex>::T_ns;
  double     JRegressor<JShowerEH, JSimplex>::Vmax_npe = std::numeric_limits<double>::max();
 
}
 
#endif