Jpp  master_rocky
the software that should make you happy
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/JRange.hh"
11 #include "JMath/JMathSupportkit.hh"
12 #include "JPhysics/JConstants.hh"
14 
15 
16 /**
17  * \author mdejong
18  */
19 
20 namespace JPHYSICS {}
21 namespace JPP { using namespace JPHYSICS; }
22 
23 namespace JPHYSICS {
24 
25  using JIO::JSerialisable;
26  using JIO::JReader;
27  using JIO::JWriter;
31  using JTOOLS::JRange;
32 
33 
34  /**
35  * Multi-dimensional PDF table for arrival time of Cherenkov light.
36  */
37  template<class JFunction1D_t,
38  class JMaplist_t,
40  class JPDFTable :
41  public JTransformableMultiFunction<JFunction1D_t, JMaplist_t, JDistance_t>,
42  public JSerialisable,
43  public JObjectBinaryIO< JPDFTable<JFunction1D_t, JMaplist_t, JDistance_t> >
44  {
45  public:
46 
48 
52 
55 
57 
61 
63 
64 
65  /**
66  * Default constructor.
67  */
70  {}
71 
72 
73  /**
74  * Constructor.
75  *
76  * \param input multi-dimensional function
77  */
78  template<class __JFunction_t, class __JMaplist_t, class __JDistance_t>
81  {}
82 
83 
84  /**
85  * Constructor.
86  *
87  * \param input multi-dimensional histogram
88  */
89  template<class JHistogram_t, class __JMaplist_t, class __JDistance_t>
92  {}
93 
94 
95  /**
96  * Blur PDF.
97  *
98  * The arrival times of Cherenkov light are smeared according to a Gaussian distribution
99  * with the specified width (i.e.\ TTS) using Gauss-Hermite integration.
100  * An exception is made when the time range according the specified quantile is
101  * smaller than the specified width (TTS) of the Gaussian distribution.
102  * In that case, the resulting PDF is a Gaussian distribution with the specified width (TTS)
103  * and normalisation according to the integral value of the input PDF.
104  * A smooth transition is imposed between the normal regime and this exeption.
105  *
106  * \param TTS TTS [ns]
107  * \param numberOfPoints number of points for Gauss-Hermite integration
108  * \param epsilon precision
109  * \param quantile quantile
110  */
111  void blur(const double TTS,
112  const int numberOfPoints = 25,
113  const double epsilon = 1.0e-10,
114  const double quantile = 0.99)
115  {
116  using namespace std;
117  using namespace JPP;
118 
119  typedef typename transformer_type::array_type array_type;
120 
121  const JGaussHermite engine(numberOfPoints, epsilon);
122 
123  for (super_iterator i = this->super_begin(); i != this->super_end(); ++i) {
124 
125  const array_type array = (*i).getKey();
126  function_type& f1 = (*i).getValue();
127 
128  if (!f1.empty()) {
129 
130  const typename function_type::supervisor_type& supervisor = f1.getSupervisor();
131 
132  const JMultiMapGetTransformer<NUMBER_OF_DIMENSIONS - 1, value_type> get(*(this->transformer), array);
133  const JMultiMapPutTransformer<NUMBER_OF_DIMENSIONS - 1, value_type> put(*(this->transformer), array);
134 
135  f1.transform(get);
136  f1.compile();
137 
138  const JQuantiles Q(f1, quantile);
139 
140  // abscissa
141 
142  JSet<double> X;
143 
144  for (JGaussHermite::const_iterator j = engine.begin(); j != engine.end(); ++j) {
145  X.insert(Q.getX() + TTS*sqrt(2.0)*j->getX());
146  }
147 
148  for (typename function_type::const_iterator j = f1.begin(); j != f1.end(); ++j) {
149 
150  if (j->getX() - TTS < X.getXmin()) {
151  X.insert(j->getX() - TTS);
152  }
153 
154  if (j->getX() + TTS > X.getXmax()) {
155  X.insert(j->getX() + TTS);
156  }
157  }
158 
159 
160  const double W = gauss(Q.getUpperLimit() - Q.getLowerLimit(), TTS);
161 
163 
164  for (JSet<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
165 
166  double y = 0.0;
167 
168  for (JGaussHermite::const_iterator j = engine.begin(); j != engine.end(); ++j) {
169 
170  const double u = j->getX();
171  const double v = j->getY() / sqrt(PI);
172  const double w = get_value(f1(*x + u*TTS*sqrt(2.0)));
173 
174  y += v * w;
175  }
176 
177  buffer[*x] = W * Q.getIntegral() * Gauss(*x - Q.getX(), TTS) + (1.0 - W) * y;
178  }
179 
180  buffer.transform(put);
181  buffer.compile();
182 
183  f1 = buffer;
184 
185  f1.setExceptionHandler(supervisor);
186  }
187  }
188  }
189 
190 
191  /**
192  * Compresses PDF to given abscissa range.
193  *
194  * \param range abscissa range
195  */
197  {
198  for (super_iterator i = this->super_begin(); i != this->super_end(); ++i) {
199 
200  function_type& f1 = i.getValue();
201  function_type f2;
202 
203  typename function_type::iterator p = f1.lower_bound(range.getLowerLimit());
204  typename function_type::iterator q = f1.lower_bound(range.getUpperLimit());
205 
206  for (++q; p != q; ++p) {
207  f2.insert(*p);
208  }
209 
210  f2.swap(f1);
211  }
212 
213  this->compile();
214  }
215 
216 
217  /**
218  * Read from input.
219  *
220  * \param in reader
221  * \return reader
222  */
223  virtual JReader& read(JReader& in) override
224  {
225  if (in >> static_cast<transformablemultifunction_type&>(*this)) {
226 
227  // read optional transformer
228 
230 
231  if (buffer.read(in)) {
232 
233  this->transformer.reset(buffer.clone());
234 
235  } else {
236 
237  in.clear();
238 
240  }
241  }
242 
243  this->compile();
244 
245  return in;
246  }
247 
248 
249  /**
250  * Write from input.
251  *
252  * \param out writer
253  * \return writer
254  */
255  virtual JWriter& write(JWriter& out) const override
256  {
257  out << static_cast<const transformablemultifunction_type&>(*this);
258 
259  this->transformer->write(out);
260 
261  return out;
262  }
263  };
264 }
265 
266 #endif
Auxiliary methods for mathematics.
Physics constants.
Auxiliary class to define a range between two values.
int numberOfPoints
Definition: JResultPDF.cc:22
Interface for binary input.
virtual void clear()
Clear status of reader.
Forward declaration of binary output.
Interface for binary output.
Multi-dimensional PDF table for arrival time of Cherenkov light.
Definition: JPDFTable.hh:44
transformablemultifunction_type::argument_type argument_type
Definition: JPDFTable.hh:49
void compress(const JRange< typename function_type::abscissa_type > &range)
Compresses PDF to given abscissa range.
Definition: JPDFTable.hh:196
transformablemultifunction_type::result_type result_type
Definition: JPDFTable.hh:50
transformablemultifunction_type::transformer_type transformer_type
Definition: JPDFTable.hh:54
JPDFTable()
Default constructor.
Definition: JPDFTable.hh:68
JPDFTable(const JTransformableMultiFunction< __JFunction_t, __JMaplist_t, __JDistance_t > &input)
Constructor.
Definition: JPDFTable.hh:79
transformablemultifunction_type::value_type value_type
Definition: JPDFTable.hh:51
transformablemultifunction_type::super_iterator super_iterator
Definition: JPDFTable.hh:59
virtual JWriter & write(JWriter &out) const override
Write from input.
Definition: JPDFTable.hh:255
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:111
JTransformableMultiFunction< JFunction1D_t, JMaplist_t, JDistance_t > transformablemultifunction_type
Definition: JPDFTable.hh:47
JPDFTable(const JTransformableMultiHistogram< JHistogram_t, __JMaplist_t, __JDistance_t > &input)
Constructor.
Definition: JPDFTable.hh:90
virtual JReader & read(JReader &in) override
Read from input.
Definition: JPDFTable.hh:223
transformablemultifunction_type::multimap_type multimap_type
Definition: JPDFTable.hh:53
transformablemultifunction_type::function_type function_type
Definition: JPDFTable.hh:60
transformablemultifunction_type::super_const_iterator super_const_iterator
Definition: JPDFTable.hh:58
Template definition of transformer of the probability density function (PDF) of the time response of ...
One dimensional array of template objects with fixed length.
Definition: JArray.hh:43
container_type::const_iterator const_iterator
Definition: JCollection.hh:91
Numerical integrator for .
Definition: JQuadrature.hh:251
Auxiliary class to convert JMultiMapTransformer to JCollectionElementTransformer.
Auxiliary class to convert JMultiMapTransformer to JCollectionElementTransformer.
Interface for weight application and coordinate transformation of function.
static JMultiMapTransformer * getClone()
Get clone of default transformer.
Multidimensional map.
Definition: JMultiMap.hh:52
Quantile calculator for a given function.
Definition: JQuantiles.hh:108
double getIntegral() const
Get integral of function.
Definition: JQuantiles.hh:346
double getX() const
Get position of maximum.
Definition: JQuantiles.hh:313
Range of values.
Definition: JRange.hh:42
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
Transformable multidimensional function.
JLANG::JSharedPointer< transformer_type > transformer
multifunction_type::result_type result_type
multifunction_type::super_const_iterator super_const_iterator
multifunction_type::super_iterator super_iterator
multifunction_type::argument_type argument_type
void insert(const JTransformableMultiFunction< __JFunction_t, __JMaplist_t, __JDistance_t > &input)
Insert multidimensional input.
Transformable multidimensional histogram.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
const double epsilon
Definition: JQuadrature.cc:21
double Gauss(const double x, const double sigma)
Normalised Gauss function.
static const double PI
Mathematical constants.
double gauss(const double x, const double sigma)
Gauss function (normalised to 1 at x = 0).
Auxiliary methods for light properties of deep-sea water.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
data_type w[N+1][M+1]
Definition: JPolint.hh:867
double u[N+1]
Definition: JPolint.hh:865
JResultEvaluator< JResult_t >::result_type get_value(const JResult_t &value)
Helper method to recursively evaluate a to function value.
Definition: JResult.hh:998
data_type v[N+1][M+1]
Definition: JPolint.hh:866
int j
Definition: JPolint.hh:792
Definition: JSTDTypes.hh:14
Auxiliary base class for storing and loading a single object to and from a binary file,...
double getValue(const double x) const
Function value.
Definition: JMathlib.hh:1421
Template class for distance evaluation.
Definition: JDistance.hh:24
Simple data structure for an abstract collection of non-equidistant abscissa values.
Definition: JSet.hh:34
virtual abscissa_type getXmin() const override
Get minimal abscissa value.
Definition: JSet.hh:105
virtual abscissa_type getXmax() const override
Get maximal abscissa value.
Definition: JSet.hh:116