Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JPDFTable.hh
Go to the documentation of this file.
1 #ifndef __JPHYSICS__JPDFTABLE__
2 #define __JPHYSICS__JPDFTABLE__
3 
4 #include <cmath>
5 
6 #include "JIO/JObjectBinaryIO.hh"
8 #include "JTools/JQuantiles.hh"
9 #include "JTools/JSet.hh"
10 #include "JTools/JConstants.hh"
12 
13 
14 /**
15  * \author mdejong
16  */
17 
18 namespace JPHYSICS {}
19 namespace JPP { using namespace JPHYSICS; }
20 
21 namespace JPHYSICS {
22 
23  using JIO::JSerialisable;
24  using JIO::JReader;
25  using JIO::JWriter;
29 
30 
31  /**
32  * Multi-dimensional PDF table for arrival time of Cherenkov light.
33  */
34  template<class JFunction1D_t,
35  class JMaplist_t,
37  class JPDFTable :
38  public JTransformableMultiFunction<JFunction1D_t, JMaplist_t, JDistance_t>,
39  public JSerialisable,
40  public JObjectBinaryIO< JPDFTable<JFunction1D_t, JMaplist_t, JDistance_t> >
41  {
42  public:
43 
45 
49 
52 
54 
58 
60 
61 
62  /**
63  * Default constructor.
64  */
67  {}
68 
69 
70  /**
71  * Constructor.
72  *
73  * \param input multi-dimensional function
74  */
75  template<class JPDF_t, class JPDFMaplist_t, class JPDFDistance_t>
78  {}
79 
80 
81  /**
82  * Constructor.
83  *
84  * \param input multi-dimensional histogram
85  */
86  template<class JHistogram1D_t, class JHistogramMaplist_t, class JHistogramDistance_t>
89  {}
90 
91 
92  /**
93  * Blur PDF.
94  *
95  * The arrival times of Cherenkov light are smeared according to a Gaussian distribution
96  * with the specified width (i.e. TTS) using Gauss-Hermite integration.
97  * An exception is made when the time range according the specified quantile is
98  * smaller than the specified width (TTS) of the Gaussian distribution.
99  * In that case, the resulting PDF is a Gaussian distribution with the specified width (TTS)
100  * and normalisation according to the integral value of the input PDF.
101  * A smooth transition is imposed between the normal regime and this exeption.
102  *
103  * \param TTS TTS [ns]
104  * \param numberOfPoints number of points for Gauss-Hermite integration
105  * \param epsilon precision
106  * \param quantile quantile
107  */
108  void blur(const double TTS,
109  const int numberOfPoints = 25,
110  const double epsilon = 1.0e-10,
111  const double quantile = 0.99)
112  {
113  using namespace std;
114  using namespace JTOOLS;
115 
116  typedef typename transformer_type::array_type array_type;
117 
118  const JGaussHermite engine(numberOfPoints, epsilon);
119 
120  for (super_iterator i = this->super_begin(); i != this->super_end(); ++i) {
121 
122  const array_type array = (*i).getKey();
123  function_type& f1 = (*i).getValue();
124 
125  if (!f1.empty()) {
126 
127  const typename function_type::supervisor_type& supervisor = f1.getSupervisor();
128 
129  const JMultiMapGetTransformer<NUMBER_OF_DIMENSIONS - 1, value_type> get(*(this->transformer), array);
130  const JMultiMapPutTransformer<NUMBER_OF_DIMENSIONS - 1, value_type> put(*(this->transformer), array);
131 
132  f1.transform(get);
133  f1.compile();
134 
135  const JQuantiles Q(f1, quantile);
136 
137  // abscissa
138 
139  JSet<double> X;
140 
141  for (JGaussHermite::const_iterator j = engine.begin(); j != engine.end(); ++j) {
142  X.insert(Q.getX() + TTS*sqrt(2.0)*j->getX());
143  }
144 
145  for (typename function_type::const_iterator j = f1.begin(); j != f1.end(); ++j) {
146 
147  if (j->getX() - TTS < X.getXmin()) {
148  X.insert(j->getX() - TTS);
149  }
150 
151  if (j->getX() + TTS > X.getXmax()) {
152  X.insert(j->getX() + TTS);
153  }
154  }
155 
156 
157  const double W = gauss(Q.getUpperLimit() - Q.getLowerLimit(), TTS);
158 
160 
161  for (JSet<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
162 
163  double y = 0.0;
164 
165  for (JGaussHermite::const_iterator j = engine.begin(); j != engine.end(); ++j) {
166 
167  const double u = j->getX();
168  const double v = j->getY() / sqrt(PI);
169  const double w = get_value(f1(*x + u*TTS*sqrt(2.0)));
170 
171  y += v * w;
172  }
173 
174  buffer[*x] = W * Q.getIntegral() * Gauss(*x - Q.getX(), TTS) + (1.0 - W) * y;
175  }
176 
177  buffer.transform(put);
178  buffer.compile();
179 
180  f1 = buffer;
181 
182  f1.setExceptionHandler(supervisor);
183  }
184  }
185  }
186 
187 
188  /**
189  * Read from input.
190  *
191  * \param in reader
192  * \return reader
193  */
194  virtual JReader& read(JReader& in)
195  {
196  if (in >> static_cast<transformablemultifunction_t&>(*this)) {
197 
198  // read optional transformer
199 
201 
202  if (buffer.read(in)) {
203 
204  this->transformer.reset(buffer.clone());
205 
206  } else {
207 
208  in.clear();
209 
211  }
212  }
213 
214  this->compile();
215 
216  return in;
217  }
218 
219 
220  /**
221  * Write from input.
222  *
223  * \param out writer
224  * \return writer
225  */
226  virtual JWriter& write(JWriter& out) const
227  {
228  out << static_cast<const transformablemultifunction_t&>(*this);
229 
230  this->transformer->write(out);
231 
232  return out;
233  }
234 
235  protected:
236  /**
237  * Gauss function (normalised to 1 at x = 0).
238  *
239  * \param x x
240  * \param sigma sigma
241  * \return function value
242  */
243  static double gauss(const double x, const double sigma)
244  {
245  const double u = x / sigma;
246 
247  if (fabs(u) < 10.0)
248  return exp(-0.5*u*u);
249  else
250  return 0.0;
251  }
252 
253 
254  /**
255  * Normalised Gauss function.
256  *
257  * \param x x
258  * \param sigma sigma
259  * \return function value
260  */
261  static double Gauss(const double x, const double sigma)
262  {
263  return gauss(x, sigma) / sqrt(2.0*JTOOLS::PI) / sigma;
264  }
265 
266 
267  /**
268  * Normalised Gauss function.
269  *
270  * \param x x
271  * \param x0 central value
272  * \param sigma sigma
273  * \return function value
274  */
275  static double Gauss(const double x, const double x0, const double sigma)
276  {
277  return Gauss(x - x0, sigma);
278  }
279  };
280 }
281 
282 #endif
JFunction_t::argument_type argument_type
void compile()
Compilation.
transformablemultifunction_t::value_type value_type
Definition: JPDFTable.hh:48
Interface for binary output.
transformablemultifunction_t::super_iterator super_iterator
Definition: JPDFTable.hh:56
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:180
multimap_type::super_iterator super_iterator
transformablemultifunction_t::result_type result_type
Definition: JPDFTable.hh:47
double getIntegral() const
Get integral of function.
Definition: JQuantiles.hh:346
Auxiliary base class for storing and loading a single object to and from a binary file...
transformablemultifunction_t::multimap_type multimap_type
Definition: JPDFTable.hh:50
Transformable multidimensional function.
Template class for distance evaluation.
Definition: JDistance.hh:24
JFunction_t::value_type value_type
Numerical integrator for W(x) = e^-(x^2).
Definition: JQuadrature.hh:253
Template definition of transformer of the Probability Density Functions of the time response of a PMT...
double getX() const
Get position of maximum.
Definition: JQuantiles.hh:313
static const double PI
Constants.
Definition: JConstants.hh:20
void blur(const double TTS, const int numberOfPoints=25, const double epsilon=1.0e-10, const double quantile=0.99)
Blur PDF.
Definition: JPDFTable.hh:108
transformablemultifunction_t::argument_type argument_type
Definition: JPDFTable.hh:46
virtual JWriter & write(JWriter &out) const
Write from input.
Definition: JPDFTable.hh:226
virtual JReader & read(JReader &in)
Read from input.
Definition: JPDFTable.hh:194
JTransformableMultiFunction< JFunction1D_t, JMaplist_t, JDistance_t > transformablemultifunction_t
Definition: JPDFTable.hh:44
virtual abscissa_type getXmax() const
Get maximal abscissa value.
Definition: JSet.hh:114
transformablemultifunction_t::super_const_iterator super_const_iterator
Definition: JPDFTable.hh:55
Auxiliary class to convert JMultiMapTransformer to JCollectionElementTransformer. ...
Quantile calculator for a given function.
Definition: JQuantiles.hh:106
static double Gauss(const double x, const double sigma)
Normalised Gauss function.
Definition: JPDFTable.hh:261
Forward declaration of binary output.
JPDFTable(const JTransformableMultiFunction< JPDF_t, JPDFMaplist_t, JPDFDistance_t > &input)
Constructor.
Definition: JPDFTable.hh:76
Constants.
JPDFTable(const JTransformableMultiHistogram< JHistogram1D_t, JHistogramMaplist_t, JHistogramDistance_t > &input)
Constructor.
Definition: JPDFTable.hh:87
JArray< NUMBER_OF_DIMENSIONS, argument_type > buffer
Transformable multidimensional histogram.
Multi-dimensional PDF table for arrival time of Cherenkov light.
Definition: JPDFTable.hh:37
JPDFTable()
Default constructor.
Definition: JPDFTable.hh:65
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:191
static double Gauss(const double x, const double x0, const double sigma)
Normalised Gauss function.
Definition: JPDFTable.hh:275
static JMultiMapTransformer * getClone()
Get clone of default transformer.
Interface for binary input.
int numberOfPoints
Definition: JResultPDF.cc:22
Interface for weight application and coordinate transformation of function.
transformablemultifunction_t::transformer_type transformer_type
Definition: JPDFTable.hh:51
container_type::const_iterator const_iterator
Definition: JCollection.hh:89
Auxiliary class to convert JMultiMapTransformer to JCollectionElementTransformer. ...
void insert(const JTransformableMultiFunction< JPDF_t, JPDFMaplist_t, JPDFDistance_t > &input)
Insert multidimensional input.
virtual abscissa_type getXmin() const
Get minimal abscissa value.
Definition: JSet.hh:103
Simple data structure for an abstract collection of non-equidistant abscissa values.
Definition: JSet.hh:29
JLANG::JSharedPointer< transformer_type > transformer
virtual void clear()
Clear status of reader.
JResultEvaluator< JResult_t >::result_type get_value(const JResult_t &value)
Helper method to recursively evaluate a to function value.
Definition: JResult.hh:748
static double gauss(const double x, const double sigma)
Gauss function (normalised to 1 at x = 0).
Definition: JPDFTable.hh:243
transformablemultifunction_t::function_type function_type
Definition: JPDFTable.hh:57
Multidimensional map.
Definition: JMultiMap.hh:46