JCDFTable.hh

Go to the documentation of this file.

#ifndef __JPHYSICS__JCDFTABLE__
#define __JPHYSICS__JCDFTABLE__

#include <stdarg.h>

#include "JIO/JSerialisable.hh"
#include "JIO/JObjectBinaryIO.hh"
#include "JLang/JException.hh"
#include "JLang/JSharedPointer.hh"
#include "JMath/JZero.hh"
#include "JTools/JArray.hh"
#include "JTools/JMultiKey.hh"
#include "JTools/JTransformableMultiFunction.hh"
#include "JTools/JTransformableMultiHistogram.hh"
#include "JTools/JConstantFunction1D.hh"
#include "JTools/JToolsToolkit.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 JLANG::JException;
  using JTOOLS::JFunctional;
  using JTOOLS::JTransformableMultiFunction;
  using JTOOLS::JTransformableMultiHistogram;
  using JTOOLS::JMultiMap;
  using JTOOLS::JMultiKey;


  /**
   * Multi-dimensional CDF table for arrival time of Cherenkov light.
   * This class can be used to determine the number of photo-electrons as a function of
   * the values of the leading parameter values and to generate random photon arrival times.
   *
   * N.B. The transformation of the PDF is assumed to be linear.
   */
  template<class JFunction1D_t,
           class JMaplist_t,
           class JDistance_t = JTOOLS::JDistance<typename JFunction1D_t::argument_type> >
  class JCDFTable :
    public JSerialisable,
    public JObjectBinaryIO< JCDFTable<JFunction1D_t, JMaplist_t, JDistance_t> >,
    public JFunctional<>
  {
  public:

    typedef typename JFunction1D_t::argument_type                                       argument_type;
    typedef typename JFunction1D_t::result_type                                         result_type;
    typedef typename JFunction1D_t::value_type                                          value_type;

    typedef JTransformableMultiFunction<JFunction1D_t, JMaplist_t, JDistance_t>         transformablemultifunction_t;

    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 JMultiKey<NUMBER_OF_DIMENSIONS - 1, argument_type>                          multikey_t;

    typedef JTOOLS::JConstantFunction1D<double, argument_type>                          JConstantFunction1D_t;
    typedef JTOOLS::JMultiFunction<JConstantFunction1D_t, JMaplist_t, JDistance_t>      JMultiQuantile_t;
    typedef JTOOLS::JMultiFunction<JFunction1D_t,         JMaplist_t, JDistance_t>      JMultiFunction_t;


    /**
     * Default constructor.
     */
    JCDFTable() :
      transformer(transformer_type::getClone())
    {}


    /**
     * Constructor.
     *
     * \param  input              multi-dimensional PDF
     * \param  eps                minimal step size for CDF
     */
    template<class JPDF_t, class JPDFMaplist_t, class JPDFDistance_t>
    JCDFTable(const JTransformableMultiFunction<JPDF_t, JPDFMaplist_t, JPDFDistance_t>& input,
              const typename JPDF_t::ordinate_type                                      eps = JMATH::zero) :
      transformer(transformer_type::getClone())
    {
      this->transformer.reset(input.transformer->clone());

      for (typename JTransformableMultiFunction<JPDF_t, JPDFMaplist_t, JPDFDistance_t>::super_const_iterator i = input.super_begin(); i != input.super_end(); ++i) {
        this->insert((*i).getKey(), (*i).getValue(), eps);
      }

      this->compile();
    }

    
    /**
     * Constructor.
     *
     * \param  input              multi-dimensional histogram
     * \param  eps                minimal step size for CDF
     */
    template<class JHistogram1D_t, class JHistogramMap_t, class JHistogramDistance_t>    
    JCDFTable(const JTransformableMultiHistogram<JHistogram1D_t, JHistogramMap_t, JHistogramDistance_t>& input,
              const typename JHistogram1D_t::ordinate_type                                               eps = JMATH::zero) :
      transformer(transformer_type::getClone())
    {
      this->transformer.reset(input.transformer->clone());
      
      this->insert(JMultiKey<0, typename JHistogram1D_t::abscissa_type>(), input, eps);
      
      this->compile();
    }


    /**
     * Application of weight function.
     *
     * \param  transformer        function transformer
     */
    template<class JFunctionTransformer_t>
    void transform(const JFunctionTransformer_t& transformer)
    {
      for (typename JMultiQuantile_t::super_iterator i = intensity.super_begin(); i != intensity.super_end(); ++i) {

        const typename transformer_type::array_type array    = (*i).getKey();
        JConstantFunction1D_t&                      function = (*i).getValue();

        const double W1 = this->transformer->getWeight(array);
        const double W2 =       transformer .getWeight(array);

        const result_type y  = function(JMATH::zero);

        function = JConstantFunction1D_t(y * W1 / W2);
      }

      this->transformer.reset(transformer.clone());
      this->compile();
    }


    /**
     * Get number of photo-electrons.
     *
     * \param  x                  comma separated argument list
     * \return                    number of photo-electrons
     */
    double getNPE(const argument_type x, ...) const
    {
      va_start(ap, x);

      buffer[0] = x;

      for (int i = 1; i != NUMBER_OF_DIMENSIONS - 1; ++i) {
        buffer[i] = va_arg(ap, argument_type);
      }

      va_end(ap);
 
      const double W   = transformer->getWeight(buffer);
      const double npe = intensity.evaluate(buffer.begin());

      return W * npe;
    }


    /**
     * Generate arrival time.
     *
     * \param  x                  comma seperated argument list (last value is random number between 0 and 1)
     * \return                    arrival time
     */
    double getTime(const argument_type x, ...) const
    {
      va_start(ap, x);

      buffer[0] = x;

      for (int i = 1; i != NUMBER_OF_DIMENSIONS; ++i) {
        buffer[i] = va_arg(ap, argument_type);
      }

      va_end(ap);

      const argument_type y = function.evaluate(buffer.begin());

      return transformer->getXn(buffer, y);
    }


    /**
     * Read CDF from input.
     *
     * \param  in                 reader
     * \return                    reader
     */
    JReader& read(JReader& in)
    {
      in >> intensity;
      in >> function;

      JPDFTransformer<NUMBER_OF_DIMENSIONS - 1, argument_type> buffer;

      if (buffer.read(in))
        transformer.reset(buffer.clone());
      else
        transformer.reset(transformer_type::getClone());

      compile();

      intensity.setExceptionHandler(new typename JMultiQuantile_t::function_type::JDefaultResult(JMATH::zero));
      function .setExceptionHandler(new typename JMultiFunction_t::function_type::JDefaultResult(JMATH::zero));

      return in;
    }


    /**
     * Write CDF to output.
     *
     * \param  out                writer
     * \return                    writer
     */
    virtual JWriter& write(JWriter& out) const
    {
      out << intensity;
      out << function;

      return transformer->write(out);
    }


    JMultiQuantile_t                intensity;     // integrated PDF
    JMultiFunction_t                function;      // normalised CDF
    JLANG::JSharedPointer<transformer_type> transformer;


  protected:
    /**
     * Function compilation.
     */
    virtual void do_compile() 
    {
      intensity.compile();
      function .compile();
    }


    /**
     * Insert value at given multidimensional key.
     *
     * \param  key                multi-dimensional key
     * \param  value              function or histogram
     * \param  eps                minimal step size for CDF
     */
    template<class JPDF_t>
    void insert(const multikey_t&                    key,
                const JPDF_t&                        value,
                const typename JPDF_t::ordinate_type eps)
    {
      try {

        const typename transformer_type::array_type array(key);
        
        JFunction1D_t buffer;
        
        const argument_type z = transformer->getXn(array, 1.0) - transformer->getXn(array, 0.0);
        const result_type   V = JTOOLS::makeCDF(value, buffer, eps);
        
        intensity.insert(key, JConstantFunction1D_t(z*V));
        function .insert(key, buffer);
      } 
      catch(const JException& error) {
        intensity.insert(key, JMATH::zero);
      }
    }


    /**
     * Insert multi-dimensional histogram at multi-dimensional key.
     *
     * \param  key                multidimensional key
     * \param  histogram          multidimensional histogram
     * \param  eps                minimal step size for CDF
     */
    template<unsigned int N, class JKey_t, class JValue_t, class JHistogramMaplist_t, class JHistogramDistance_t>
    void insert(const JMultiKey<N, JKey_t>&                                                   key,
                const JMultiMap<JKey_t, JValue_t, JHistogramMaplist_t, JHistogramDistance_t>& histogram,
                const typename JValue_t::contents_type                                        eps)
    {
      if (histogram.getSize() > 1) {

        for (typename JMultiMap<JKey_t, JValue_t, JHistogramMaplist_t, JHistogramDistance_t>::const_iterator j = histogram.begin(), i = j++; j != histogram.end(); ++i, ++j) {
          
          const JKey_t x = 0.5 * (i->first + j->first);
          
          insert(JMultiKey<N+1, JKey_t>(key, x), i->second, eps);
        } 
      }
    }


    mutable va_list ap;
    mutable JTOOLS::JArray<NUMBER_OF_DIMENSIONS, argument_type> buffer;
  };
}

#endif