JCDFTable2D.hh

Go to the documentation of this file.

#ifndef __JSIRENE__JCDFTABLE2D__
#define __JSIRENE__JCDFTABLE2D__

#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;
    }


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

#endif