#ifndef __JPHYSICS__JCDFTABLE__
#define __JPHYSICS__JCDFTABLE__
#include <memory>
#include "JIO/JSerialisable.hh"
#include "JIO/JObjectBinaryIO.hh"
#include "JLang/JException.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 "JTools/JTransformableMultiHistogram.hh"
#include "JTools/JHistogramMap.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;
using JTOOLS::JArray;
using JTOOLS::JHistogramMap;
/**
* 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_type;
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 __JFunction_t, class __JMaplist_t, class __JDistance_t>
JCDFTable(const JTransformableMultiFunction<__JFunction_t, __JMaplist_t, __JDistance_t>& input,
const typename __JFunction_t::ordinate_type eps = JMATH::zero) :
transformer(transformer_type::getClone())
{
this->transformer.reset(input.transformer->clone());
for (auto 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 JHistogram_t, class __JMaplist_t, class __JDistance_t>
JCDFTable(const JTransformableMultiHistogram<JHistogram_t, __JMaplist_t, __JDistance_t>& input,
const typename JHistogram_t::ordinate_type eps = JMATH::zero) :
transformer(transformer_type::getClone())
{
this->transformer.reset(input.transformer->clone());
this->insert(JMultiKey<0, typename JHistogram_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 args comma separated argument list
* \return number of photo-electrons
*/
template<class ...Args>
double getNPE(const Args& ...args) const
{
const JArray<NUMBER_OF_DIMENSIONS - 1, argument_type> buffer(args...);
const double W = transformer->getWeight(buffer);
const double npe = intensity.evaluate(buffer.data());
return W * npe;
}
/**
* Generate arrival time.
*
* \param args comma seperated argument list (last value is random number between 0 and 1)
* \return arrival time
*/
template<class ...Args>
double getTime(const Args& ...args) const
{
const JArray<NUMBER_OF_DIMENSIONS, argument_type> buffer(args...);
const argument_type y = function.evaluate(buffer.data());
return transformer->getXn(buffer, y);
}
/**
* Read CDF from input.
*
* \param in reader
* \return reader
*/
virtual JReader& read(JReader& in) override
{
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 override
{
out << intensity;
out << function;
return transformer->write(out);
}
JMultiQuantile_t intensity; // integrated PDF
JMultiFunction_t function; // normalised CDF
std::shared_ptr<transformer_type> transformer;
protected:
/**
* Function compilation.
*/
virtual void do_compile() override
{
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 method JTOOLS::makeCDF
*/
template<class JValue_t>
void insert(const multikey_type& key,
const JValue_t& value,
const typename JValue_t::ordinate_type eps)
{
using namespace JPP;
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 = 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 value multidimensional histogram
* \param eps minimal step size for CDF
*/
template<unsigned int N,
class __JAbscissa_t,
class __JContents_t,
template<class, class, class> class __JMap_t,
class __JDistance_t>
void insert(const JMultiKey<N, argument_type>& key,
const JHistogramMap<__JAbscissa_t, __JContents_t, __JMap_t, __JDistance_t>& value,
const __JContents_t eps)
{
if (value.getSize() > 1) {
for (auto j = value.begin(), i = j++; j != value.end(); ++i, ++j) {
const __JAbscissa_t x = 0.5 * (i->first + j->first);
insert(JMultiKey<N+1, __JAbscissa_t>(key, x), i->second, eps);
}
}
}
};
}
#endif