Jpp  17.3.0-rc.2
the software that should make you happy
 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/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 
202  typename function_type::iterator p = f1.lower_bound(range.getLowerLimit());
203 
204  f1.function_type::container_type::erase(f1.begin(), p);
205 
206  typename function_type::iterator q = f1.lower_bound(range.getUpperLimit());
207 
208  f1.function_type::container_type::erase(++q, f1.end());
209  }
210 
211  this->compile();
212  }
213 
214 
215  /**
216  * Read from input.
217  *
218  * \param in reader
219  * \return reader
220  */
221  virtual JReader& read(JReader& in) override
222  {
223  if (in >> static_cast<transformablemultifunction_type&>(*this)) {
224 
225  // read optional transformer
226 
228 
229  if (buffer.read(in)) {
230 
231  this->transformer.reset(buffer.clone());
232 
233  } else {
234 
235  in.clear();
236 
238  }
239  }
240 
241  this->compile();
242 
243  return in;
244  }
245 
246 
247  /**
248  * Write from input.
249  *
250  * \param out writer
251  * \return writer
252  */
253  virtual JWriter& write(JWriter& out) const override
254  {
255  out << static_cast<const transformablemultifunction_type&>(*this);
256 
257  this->transformer->write(out);
258 
259  return out;
260  }
261  };
262 }
263 
264 #endif
function_type::argument_type argument_type
void insert(const JTransformableMultiFunction< __JFunction_t, __JMaplist_t, __JDistance_t > &input)
Insert multidimensional input.
void compile()
Compilation.
virtual JWriter & write(JWriter &out) const override
Write from input.
Definition: JPDFTable.hh:253
data_type w[N+1][M+1]
Definition: JPolint.hh:778
Interface for binary output.
Q(UTCMax_s-UTCMin_s)-livetime_s
multifunction_type::result_type result_type
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
multimap_type::super_iterator super_iterator
Auxiliary base class for storing and loading a single object to and from a binary file...
Auxiliary methods for mathematics.
transformablemultifunction_type::multimap_type multimap_type
Definition: JPDFTable.hh:53
Transformable multidimensional function.
JMultiMapTransformer< JMapLength< JMaplist_t >::value, argument_type > transformer_type
Template class for distance evaluation.
Definition: JDistance.hh:24
multifunction_type::argument_type argument_type
transformablemultifunction_type::function_type function_type
Definition: JPDFTable.hh:60
Template definition of transformer of the probability density function (PDF) of the time response of ...
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
transformablemultifunction_type::value_type value_type
Definition: JPDFTable.hh:51
JTransformableMultiFunction< JFunction1D_t, JMaplist_t, JDistance_t > transformablemultifunction_type
Definition: JPDFTable.hh:47
JMultiMap< typename JFunction_t::argument_type, JFunction_t, JMaplist_t, JDistance_t > multimap_type
void compress(const JRange< typename function_type::abscissa_type > &range)
Compresses PDF to given abscissa range.
Definition: JPDFTable.hh:196
Forward declaration of binary output.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
function_type::value_type value_type
JArray< NUMBER_OF_DIMENSIONS, argument_type > buffer
Transformable multidimensional histogram.
Multi-dimensional PDF table for arrival time of Cherenkov light.
Definition: JPDFTable.hh:40
JPDFTable()
Default constructor.
Definition: JPDFTable.hh:68
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
Physics constants.
static const double PI
Mathematical constants.
static JMultiMapTransformer * getClone()
Get clone of default transformer.
Interface for binary input.
One dimensional array of template objects with fixed length.
Definition: JArray.hh:40
Range of values.
Definition: JRange.hh:38
JPDFTable(const JTransformableMultiFunction< __JFunction_t, __JMaplist_t, __JDistance_t > &input)
Constructor.
Definition: JPDFTable.hh:79
transformablemultifunction_type::transformer_type transformer_type
Definition: JPDFTable.hh:54
transformablemultifunction_type::argument_type argument_type
Definition: JPDFTable.hh:49
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
virtual JReader & read(JReader &in) override
Read from input.
Definition: JPDFTable.hh:221
Auxiliary class to define a range between two values.
then set_variable NUMBER_OF_TESTS else set_variable NUMBER_OF_TESTS fi function gauss()
int numberOfPoints
Definition: JResultPDF.cc:22
double Gauss(const double x, const double sigma)
Normalised Gauss function.
transformablemultifunction_type::super_const_iterator super_const_iterator
Definition: JPDFTable.hh:58
JPDFTable(const JTransformableMultiHistogram< JHistogram_t, __JMaplist_t, __JDistance_t > &input)
Constructor.
Definition: JPDFTable.hh:90
no fit printf nominal n $STRING awk v X
container_type::const_iterator const_iterator
Definition: JCollection.hh:91
transformablemultifunction_type::super_iterator super_iterator
Definition: JPDFTable.hh:59
int j
Definition: JPolint.hh:703
transformablemultifunction_type::result_type result_type
Definition: JPDFTable.hh:50
multifunction_type::super_iterator super_iterator
data_type v[N+1][M+1]
Definition: JPolint.hh:777
double u[N+1]
Definition: JPolint.hh:776
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:46
JLANG::JSharedPointer< transformer_type > transformer
multifunction_type::super_const_iterator super_const_iterator
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:998
const double epsilon
Definition: JQuadrature.cc:21
double getValue(const double x) const
Function value.
Definition: JPolynome.hh:233