JCDFTable1D.hh

Go to the documentation of this file.

#ifndef __JSIRENE__JCDFTABLE1D__
#define __JSIRENE__JCDFTABLE1D__
 
#include <cmath>
 
#include "JPhysics/JConstants.hh"
#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 JPP;
 
      // 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), JOmega3D::range_type(0.0,0.501*PI), 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