JRegressorHelper.hh

Go to the documentation of this file.

#ifndef __JFIT_JREGRESSORHELPER__
#define __JFIT_JREGRESSORHELPER__

#include <string>
#include <array>

#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JNPETable.hh"

#include "JTools/JRange.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"

/**
 * \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;

  /**
   * forward declarations
   */
  class JLine3Z;
  template<class JModel_t>
  class JAbstractMinimiser;
  template<class JModel_t>
  class JSimplex;
  template<class JModel_t>
  class JGandalf;


  /**
   * Template data structure for PDF tables.
   */
  template<class JModel_t, template<class> class JMinimiser_t>
  struct JRegressorHelper;


  /**
   * Template specialisation for PDF tables.
   */
  template<>
  struct JRegressorHelper<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
    

    /**
     * PDF types.
     */
    static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
  };


  /**
   * PDF types.
   */
  const JPDFType_t JRegressorHelper<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 };


  
  /**
   * Template data structure for storage for PDF tables.
   */
  template<class JModel_t, template<class> class JMinimiser_t>
  struct JRegressorStorage;


  /**
   * Template specialisation for storage for PDF tables.
   */
  template<>
  struct JRegressorStorage<JLine3Z, JGandalf> :
    public JRegressorHelper<JLine3Z, JGandalf>
  {
    /**
     * Default constructor.
     */
    JRegressorStorage()
    {}


    /**
     * Constructor.
     *
     * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
     * will be replaced by the corresponding PDF types listed in JRegressor<JLine3Z, JGandalf>::pdf_t.
     *
     * 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
     */
    JRegressorStorage(const std::string& fileDescriptor,
                      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(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]);

        JPDF_t buffer;

        _pdf[i-1].swap(buffer);

        if (TTS > 0.0) {
          _pdf[i].blur(TTS, numberOfPoints, epsilon);
        }

        _npe[ i ] = JNPE_t(_pdf[i]);
      }
    }
    

    /**
     * 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);
      }
    }
    

    /**
     * Set maximal road width of PDF.
     *
     * \param  Rmax               road width [m]
     */
    inline void setRmax(const double Rmax)
    {
      const JTOOLS::JRange<typename JFunction1D_t::abscissa_type> range(0.0, Rmax);
      
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
        if (!_pdf[i].empty()) {
          _pdf[i].compress(range);
        }
      }
    }


  private:
    JPDFs_t _pdf;                       //!< PDFs
    JNPEs_t _npe;                       //!< NPEs
  };
}

#endif