JCDFTable2D.hh

Go to the documentation of this file.

#ifndef __JSIRENE__JCDFTABLE2D__
#define __JSIRENE__JCDFTABLE2D__
 
#include <memory>
#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 "JTools/JMultiGrid.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 2 dimensions.
   * This class provides an upper limit for the average number of photo-electrons
   * as a function of the distance and the cosine of the emission angle.
   */
  template<class JArgument_t, class JResult_t>
  class JCDFTable2D :
    public JMultiFunction<JPolintFunction1D<1, JElement2D<JArgument_t,JResult_t>, JGridCollection>, 
                          JMapList<JPolint1FunctionalGridMap> >
  {
  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 JMapList<JPolint1FunctionalGridMap>                    JMaplist_t;
    typedef JMultiFunction<JFunction1D_t, JMaplist_t>              JFunction2D_t;
 
    typedef JMultiMapTransformer<2, argument_type>                 transformer_type;
  
 
    /**
     * Default constructor.
     */
    JCDFTable2D()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  cdf                 CDF            (5D table)
     * \param  number_of_bins      number of bins (2D table)
     * \param  safety_factor       safety factor
     */
    template<class JCDF_t, class JCDFMaplist_t, class JCDFDistance_t>
    JCDFTable2D(const JCDFTable<JCDF_t, JCDFMaplist_t, JCDFDistance_t>& cdf,
                const int    number_of_bins,
                const double safety_factor = 2.0) :
      JFunction2D_t()
    { 
      using namespace JPP;
 
      // extract the weight function
 
      try {
 
        typedef JPDFTransformer<4, 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 2D 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 JFunction2D_t::iterator i = this->begin(); i != this->end(); ++i) {
 
        i->getY().configure(make_grid(number_of_bins, -1.0, +1.0));
 
        for (typename JFunction1D_t::iterator j = i->getY().begin(); j != i->getY().end(); ++j) {
 
          const double D  = i->getX(); 
          const double cd = j->getX();
          const double W  = transformer->getWeight(D, cd);
 
          double y = 0.0;
 
          for (JOmega3D::const_iterator dir = omega.begin(); dir != omega.end(); ++dir) {
            y += cdf.getNPE(D, cd, dir->getTheta(), fabs(dir->getPhi()));
          }
 
          y /= omega.size();
 
          j->getY() = safety_factor * y / W;
        }
      }
 
      this->compile();
    }
 
 
    /**
     * Get number of photo-electrons.
     *
     * \param  D                   distance between EM shower and PMT [m]
     * \param  cd                  cosine angle EM shower direction and EM shower - PMT position
     * \return                     number of photo-electrons
     */
    double getNPE(const double D,
                  const double cd) const
    {
      const double buffer[] = { D, cd };
 
      const double y = this->evaluate(buffer);
      const double W = this->transformer->getWeight(D, cd);
 
      return y * W;
    }
 
 
    std::shared_ptr<transformer_type> transformer;
  };
}
 
#endif