JPDFTable.hh

Go to the documentation of this file.

#ifndef __JPHYSICS__JPDFTABLE__
#define __JPHYSICS__JPDFTABLE__
 
#include <cmath>
 
#include "JIO/JObjectBinaryIO.hh"
#include "JTools/JTransformableMultiFunction.hh"
#include "JTools/JQuantiles.hh"
#include "JTools/JSet.hh"
#include "JTools/JRange.hh"
#include "JTools/JConstants.hh"
#include "JPhysics/JPDFTransformer.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }
 
namespace JPHYSICS {
 
  using JIO::JSerialisable;
  using JIO::JReader;
  using JIO::JWriter;
  using JIO::JObjectBinaryIO;
  using JTOOLS::JTransformableMultiFunction;
  using JTOOLS::JTransformableMultiHistogram;
  using JTOOLS::JRange;
  
 
  /**
   * Multi-dimensional PDF table for arrival time of Cherenkov light.
   */
  template<class JFunction1D_t,
           class JMaplist_t,
           class JDistance_t = JTOOLS::JDistance<typename JFunction1D_t::argument_type> >
  class JPDFTable :
    public JTransformableMultiFunction<JFunction1D_t, JMaplist_t, JDistance_t>,
    public JSerialisable,
    public JObjectBinaryIO< JPDFTable<JFunction1D_t, JMaplist_t, JDistance_t> >
  {
  public:
 
    typedef JTransformableMultiFunction<JFunction1D_t, JMaplist_t, JDistance_t>      transformablemultifunction_t;
 
    typedef typename transformablemultifunction_t::argument_type                     argument_type;
    typedef typename transformablemultifunction_t::result_type                       result_type;
    typedef typename transformablemultifunction_t::value_type                        value_type;
 
    typedef typename transformablemultifunction_t::multimap_type                     multimap_type;
    typedef typename transformablemultifunction_t::transformer_type                  transformer_type;
 
    enum { NUMBER_OF_DIMENSIONS = transformablemultifunction_t::NUMBER_OF_DIMENSIONS };
 
    typedef typename transformablemultifunction_t::super_const_iterator              super_const_iterator;
    typedef typename transformablemultifunction_t::super_iterator                    super_iterator;
    typedef typename transformablemultifunction_t::function_type                     function_type;
 
    using transformablemultifunction_t::insert;
 
 
    /**
     * Default constructor.
     */
    JPDFTable() :
      transformablemultifunction_t()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  input                multi-dimensional function
     */
    template<class JPDF_t, class JPDFMaplist_t, class JPDFDistance_t>    
    JPDFTable(const JTransformableMultiFunction<JPDF_t, JPDFMaplist_t, JPDFDistance_t>& input) :
      transformablemultifunction_t(input)
    {}
 
    
    /**
     * Constructor.
     *
     * \param  input                multi-dimensional histogram
     */
    template<class JHistogram1D_t, class JHistogramMaplist_t, class JHistogramDistance_t>
    JPDFTable(const JTransformableMultiHistogram<JHistogram1D_t, JHistogramMaplist_t, JHistogramDistance_t>& input) :
      transformablemultifunction_t(input)
    {}
 
 
    /**
     * Blur PDF.
     *
     * The arrival times of Cherenkov light are smeared according to a Gaussian distribution
     * with the specified width (i.e. TTS) using Gauss-Hermite integration.
     * An exception is made when the time range according the specified quantile is 
     * smaller than the specified width (TTS) of the Gaussian distribution.
     * In that case, the resulting PDF is a Gaussian distribution with the specified width (TTS)
     * and normalisation according to the integral value of the input PDF.
     * A smooth transition is imposed between the normal regime and this exeption.
     *
     * \param  TTS                  TTS [ns]
     * \param  numberOfPoints       number of points for Gauss-Hermite integration
     * \param  epsilon              precision
     * \param  quantile             quantile
     */
    void blur(const double TTS, 
              const int    numberOfPoints = 25,
              const double epsilon        = 1.0e-10,
              const double quantile       = 0.99)
    {
      using namespace std;
      using namespace JTOOLS;
 
      typedef typename transformer_type::array_type  array_type;
 
      const JGaussHermite engine(numberOfPoints, epsilon);
 
      for (super_iterator i = this->super_begin(); i != this->super_end(); ++i) {
        
        const array_type array = (*i).getKey();
        function_type&   f1    = (*i).getValue();
 
        if (!f1.empty()) {
 
          const typename function_type::supervisor_type& supervisor = f1.getSupervisor();
 
          const JMultiMapGetTransformer<NUMBER_OF_DIMENSIONS - 1, value_type> get(*(this->transformer), array);
          const JMultiMapPutTransformer<NUMBER_OF_DIMENSIONS - 1, value_type> put(*(this->transformer), array);
 
          f1.transform(get);
          f1.compile();
 
          const JQuantiles Q(f1, quantile);
 
          // abscissa 
 
          JSet<double> X;
 
          for (JGaussHermite::const_iterator j = engine.begin(); j != engine.end(); ++j) {
            X.insert(Q.getX() + TTS*sqrt(2.0)*j->getX());
          }
 
          for (typename function_type::const_iterator j = f1.begin(); j != f1.end(); ++j) {
 
            if (j->getX() - TTS < X.getXmin()) {
              X.insert(j->getX() - TTS);
            }
            
            if (j->getX() + TTS > X.getXmax()) {
              X.insert(j->getX() + TTS);
            }
          }
 
 
          const double W = gauss(Q.getUpperLimit() - Q.getLowerLimit(), TTS);
 
          function_type buffer;
 
          for (JSet<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
 
            double y = 0.0;
 
            for (JGaussHermite::const_iterator j = engine.begin(); j != engine.end(); ++j) {
 
              const double u = j->getX();
              const double v = j->getY() / sqrt(PI);
              const double w = get_value(f1(*x + u*TTS*sqrt(2.0)));
        
              y += v * w;
            }
 
            buffer[*x]  =  W * Q.getIntegral() * Gauss(*x - Q.getX(), TTS)  +  (1.0 - W) * y;
          }
 
          buffer.transform(put);
          buffer.compile();
 
          f1 = buffer;
 
          f1.setExceptionHandler(supervisor);
        }
      }
    }
 
    /**
     * Compresses PDF to abscissa range specified by Crange 
     *
     * \param  Crange               Compression Range [PDF abscissa_type]
     */
    void compress(const JRange<typename function_type::abscissa_type>& Crange)
    {  
      
      for (super_iterator pdf_it=this->super_begin(); pdf_it!=this->super_end(); ++pdf_it) {
        function_type& f1 = pdf_it.getValue();
        
        typename function_type::iterator f1_newB = f1.lower_bound(Crange.getLowerLimit());
        f1.function_type::container_type::erase(f1.begin(),f1_newB);
        typename function_type::iterator f1_newE = f1.lower_bound(Crange.getUpperLimit());
        f1.function_type::container_type::erase(++f1_newE,f1.end());
      }
      
    }
 
    /**
     * Read from input.
     *
     * \param  in                   reader
     * \return                      reader
     */
    virtual JReader& read(JReader& in)
    {
      if (in >> static_cast<transformablemultifunction_t&>(*this)) {
 
        // read optional transformer
 
        JPDFTransformer<NUMBER_OF_DIMENSIONS - 1, argument_type> buffer;
 
        if (buffer.read(in)) {
 
          this->transformer.reset(buffer.clone());
 
        } else {
 
          in.clear();
          
          this->transformer.reset(transformer_type::getClone());
        }
      }
      
      this->compile();
      
      return in;
    }
    
    
    /**
     * Write from input.
     *
     * \param  out                  writer
     * \return                      writer
     */
    virtual JWriter& write(JWriter& out) const
    {
      out << static_cast<const transformablemultifunction_t&>(*this);
 
      this->transformer->write(out);
 
      return out;
    }
 
  protected:
    /**
     * Gauss function (normalised to 1 at x = 0).
     *
     * \param  x                    x
     * \param  sigma                sigma
     * \return                      function value
     */
    static double gauss(const double x, const double sigma)
    {
      const double u = x / sigma;
 
      if (fabs(u) < 10.0)
        return exp(-0.5*u*u);
      else
        return 0.0;
    }
 
 
    /**
     * Normalised Gauss function.
     *
     * \param  x                    x
     * \param  sigma                sigma
     * \return                      function value
     */
    static double Gauss(const double x, const double sigma)
    {
      return gauss(x, sigma) / sqrt(2.0*JTOOLS::PI) / sigma;
    }
 
 
    /**
     * Normalised Gauss function.
     *
     * \param  x                    x
     * \param  x0                   central value
     * \param  sigma                sigma
     * \return                      function value
     */
    static double Gauss(const double x, const double x0, const double sigma)
    {
      return Gauss(x - x0, sigma);
    }
  };
}
 
#endif