JCDFTable1D.hh

Go to the documentation of this file.

#ifndef __JSIRENE__JCDFTABLE1D__
#define __JSIRENE__JCDFTABLE1D__

#include <cmath>

#include "JPhysics/JCDFTable.hh"
#include "JPhysics/JPDFTransformer.hh"
#include "JTools/JPolint.hh"
#include "JTools/JGridCollection.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JGeometry3D/JOmega3D.hh"
#include "JLang/JSharedPointer.hh"


/**
 * \author mdejong
 */

namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }

namespace JPHYSICS {

  using JTOOLS::JMapList;
  using JTOOLS::JElement2D;
  using JTOOLS::JPolintFunction1D;
  using JTOOLS::JGridCollection;
  using JTOOLS::JPolint1FunctionalGridMap;
  using JTOOLS::JMultiFunction;


  /**
   * Custom class for CDF table in 1 dimension.
   * This class provides an upper limit for the avegare number of photo-electrons
   * as a function of the distance.
   */
  template<class JArgument_t, class JResult_t>
  class JCDFTable1D :
    public JPolintFunction1D<1, JElement2D<JArgument_t,JResult_t>, JGridCollection>
  {
  public:

    typedef JArgument_t                                            argument_type;
    typedef JResult_t                                              result_type;

    typedef JElement2D<argument_type, result_type>                 JElement2D_t;
    typedef JPolintFunction1D<1, JElement2D_t, JGridCollection>    JFunction1D_t;

    typedef JMultiMapTransformer<1, argument_type>                 transformer_type;
  

    /**
     * Default constructor.
     */
    JCDFTable1D()
    {}


    /**
     * Constructor.
     *
     * \param  cdf                 CDF            (4D table)
     * \param  number_of_bins      number of bins (1D table)
     * \param  safety_factor       safety factor
     */
    template<class JCDF_t, class JCDFMaplist_t, class JCDFDistance_t>
    JCDFTable1D(const JCDFTable<JCDF_t, JCDFMaplist_t, JCDFDistance_t>& cdf,
                const int    number_of_bins,
                const double safety_factor = 2.0)
    {
      using namespace JTOOLS;
      using namespace JGEOMETRY3D;

      // extract the weight function

      try {

        typedef JPDFTransformer<3, argument_type> JPDFTransformer_t;

        const JPDFTransformer_t& object = dynamic_cast<const JPDFTransformer_t&>(*(cdf.transformer));

        transformer.reset(object.transformer.clone());
      }
      catch(const std::exception& error) {
        transformer.reset(transformer_type::getClone());
      }


      // build the 1D table

      this->configure(make_grid(number_of_bins, cdf.intensity.begin()->getX(), cdf.intensity.rbegin()->getX()));

      const JOmega3D omega(JDirection3D(0,1,0), JAngleRange(0.0,PI/2), 0.07*PI);

      for (typename JFunction1D_t::iterator j = this->begin(); j != this->end(); ++j) {

        const double R  = j->getX(); 
        const double W  = transformer->getWeight(R);

        double y = 0.0;

        for (JOmega3D::const_iterator dir = omega.begin(); dir != omega.end(); ++dir) {
          y += cdf.getNPE(R, dir->getTheta(), fabs(dir->getPhi()));
        }

        y /= omega.size();

        j->getY() = safety_factor * y / W;
      }

      this->compile();
    }


    /**
     * Get number of photo-electrons.
     *
     * \param  R                   distance between muon and PMT [m]
     * \return                     number of photo-electrons
     */
    double getNPE(const double R) const
    {
      const double buffer[] = { R };

      const double y = this->evaluate(buffer);
      const double W = this->transformer->getWeight(R);

      return y * W;
    }


    JLANG::JSharedPointer<transformer_type> transformer;
  };
}

#endif