Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JEnergyRegressor.hh
Go to the documentation of this file.
1#ifndef __JFIT__JENERGYREGRESSOR__
2#define __JFIT__JENERGYREGRESSOR__
3
4#include <string>
5#include <vector>
6#include <memory>
7
10#include "JPhysics/JNPETable.hh"
17#include "JFit/JRegressor.hh"
18#include "JFit/JEnergy.hh"
19#include "JFit/JNPE.hh"
20#include "JFit/JNPEHit.hh"
21#include "JFit/JMEstimator.hh"
22#include "JFit/JFitToolkit.hh"
23#include "JFit/JTimeRange.hh"
24#include "Jeep/JMessage.hh"
25
26
27/**
28 * \file
29 * Data regression method for JFIT::JEnergy.
30 * \author mdejong
31 */
32namespace JFIT {}
33namespace JPP { using namespace JFIT; }
34
35namespace JFIT {
36
45
46
47 /**
48 * Template specialisation for storage of PDF tables.
49 */
50 template<>
52 {
57
58 static const int NUMBER_OF_PDFS = 6; //!< Number of PDFs
59
60 typedef std::vector<JNPE_t> JNPEs_t; //!< NPEs
61
62
63 /**
64 * Default constructor.
65 */
68
69
70 /**
71 * Constructor.
72 *
73 * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
74 * will be replaced by the PDF types listed in JRegressorStorage<JEnergy>::pdf_t.
75 *
76 * \param fileDescriptor PDF file descriptor
77 * \param T_ns time range [ns]
78 */
79 JRegressorStorage(const std::string& fileDescriptor, const JTimeRange& T_ns = JTimeRange()) :
80 T_ns(T_ns)
81 {
82 using namespace std;
83 using namespace JPP;
84
85 typedef JSplineFunction1D_t JFunction1D_t;
87
88 const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
89
90 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
91
92 JPDF_t pdf;
93
94 const JPDFType_t type = pdf_t[i];
95 const string file_name = getFilename(fileDescriptor, type);
96
97 pdf.load(file_name.c_str());
98
99 pdf.setExceptionHandler(supervisor);
100
101 if (is_bremsstrahlung(type))
102 _YB.push_back(JNPE_t(pdf, T_ns));
103 else if (is_deltarays(type))
104 _YA.push_back(JNPE_t(pdf, T_ns));
105 else
106 _Y1.push_back(JNPE_t(pdf, T_ns));
107 }
108
109 // Add PDFs
110
111 if (_Y1.size() == 2) { _Y1[1].add(_Y1[0]); _Y1.erase(_Y1.begin()); }
112 if (_YA.size() == 2) { _YA[1].add(_YA[0]); _YA.erase(_YA.begin()); }
113 if (_YB.size() == 2) { _YB[1].add(_YB[0]); _YB.erase(_YB.begin()); }
114 }
115
116
117 /**
118 * Get light from muon.
119 *
120 * \return NPEs
121 */
122 const JNPEs_t& getY1() const
123 {
124 return _Y1;
125 }
126
127
128 /**
129 * Get light from delta-rays.
130 *
131 * \return NPEs
132 */
133 const JNPEs_t& getYA() const
134 {
135 return _YA;
136 }
137
138
139 /**
140 * Get light from EM showers.
141 *
142 * \return NPEs
143 */
144 const JNPEs_t& getYB() const
145 {
146 return _YB;
147 }
148
149
150 /**
151 * PDF types.
152 */
153 static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
154
155 JTimeRange T_ns; //!< Time window with respect to Cherenkov hypothesis [ns]
156
157 private:
158 JNPEs_t _Y1; //!< light from muon
159 JNPEs_t _YA; //!< light from delta-rays
160 JNPEs_t _YB; //!< light from EM showers
161 };
162
163
164 /**
165 * PDF types.
166 */
167 const JPDFType_t JRegressorStorage<JEnergy>::pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
168 SCATTERED_LIGHT_FROM_MUON,
169 DIRECT_LIGHT_FROM_DELTARAYS,
170 SCATTERED_LIGHT_FROM_DELTARAYS,
171 DIRECT_LIGHT_FROM_EMSHOWERS,
172 SCATTERED_LIGHT_FROM_EMSHOWERS };
173
174
175
176 /**
177 * Regressor function object for fit of muon energy.
178 */
179 template<>
180 struct JRegressor<JEnergy> :
181 public JAbstractRegressor<JEnergy>,
182 public JRegressorStorage <JEnergy>
183 {
184 using JAbstractRegressor<JEnergy>::operator();
185
186 typedef JRegressorStorage<JEnergy> storage_type;
187
188 using storage_type::getY1;
189 using storage_type::getYA;
190 using storage_type::getYB;
191
192
193 /**
194 * Constructor.
195 *
196 * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
197 * will be replaced by the PDF types listed in JRegressor<JEnergy>::pdf_t.
198 *
199 * \param fileDescriptor PDF file descriptor
200 * \param T_ns time range [ns]
201 */
202 JRegressor(const std::string& fileDescriptor, const JTimeRange& T_ns = JTimeRange()) :
203 storage_type(fileDescriptor, T_ns),
204 Y1(getY1()),
205 YA(getYA()),
206 YB(getYB()),
207 estimator(new JMEstimatorNormal())
208 {}
209
210
211 /**
212 * Constructor.
213 *
214 * \param storage PDF storage
215 */
216 JRegressor(const storage_type& storage) :
217 Y1(storage.getY1()),
218 YA(storage.getYA()),
219 YB(storage.getYB()),
220 estimator(new JMEstimatorNormal())
221 {
222 T_ns = storage.T_ns;
223 }
224
225
226 /**
227 * Fit function.
228 * This method is used to determine chi2 of given number of photo-electrons for given energy of muon.
229 *
230 * \param x energy
231 * \param npe npe
232 * \return chi2
233 */
234 inline double operator()(const JEnergy& x, const JNPEHit& npe) const
235 {
236 const double E = x.getE();
237 const double u = npe.getChi2(E);
238
239 return estimator->getRho(u);
240 }
241
242
243 /**
244 * Get number of photo-electrons from muon.
245 * Note that the PMT geometry should be relative to the muon trajectory,
246 * conform method JGEOMETRY3D::JAxis3D::transform.
247 *
248 * \param axis PMT axis
249 * \return npe
250 */
251 inline double getY1(const JAxis3D& axis) const
252 {
253 return getNPE(Y1, hypot(axis.getX(), axis.getY()), axis.getTheta(), fabs(axis.getPhi()));
254 }
255
256
257 /**
258 * Get number of photo-electrons from delta-rays.
259 * Note that the PMT geometry should be relative to the muon trajectory,
260 * conform method JGEOMETRY3D::JAxis3D::transform.
261 *
262 * \param axis PMT axis
263 * \return npe
264 */
265 inline double getYA(const JAxis3D& axis) const
266 {
267 return getNPE(YA, hypot(axis.getX(), axis.getY()), axis.getTheta(), fabs(axis.getPhi()));
268 }
269
270
271 /**
272 * Get number of photo-electrons from Bremsstrahlung.
273 * Note that the PMT geometry should be relative to the muon trajectory,
274 * conform method JGEOMETRY3D::JAxis3D::transform.
275 *
276 * \param axis PMT axis
277 * \return npe
278 */
279 inline double getYB(const JAxis3D& axis) const
280 {
281 return getNPE(YB, hypot(axis.getX(), axis.getY()), axis.getTheta(), fabs(axis.getPhi()));
282 }
283
284
285 /**
286 * Create data structure for handling light yields for PMT.
287 * Note that the PMT geometry should be relative to the muon trajectory,
288 * conform method JGEOMETRY3D::JAxis3D::transform.
289 *
290 * \param axis PMT axis
291 * \param R_Hz singles rate [Hz]
292 * \return light yields
293 */
294 inline JNPE getNPE(const JAxis3D& axis,
295 const double R_Hz) const
296 {
297 using namespace std;
298 using namespace JPP;
299
300 const double x = axis.getX();
301 const double y = axis.getY();
302 const double R = sqrt(x*x + y*y);
303 const double z = axis.getZ() - R / getTanThetaC();
304
305 const double theta = axis.getTheta();
306 const double phi = fabs(axis.getPhi());
307
308 const double y1 = getNPE(Y1, R, theta, phi);
309 const double yA = getNPE(YA, R, theta, phi);
310 const double yB = getNPE(YB, R, theta, phi);
311
312 return JNPE(getN(T_ns, R_Hz * 1.0e-9), y1, yA, yB, z);
313 }
314
315
316 /**
317 * Get maximal road width of NPE.
318 *
319 * \return road width [m]
320 */
321 inline double getRmax() const
322 {
323 return std::max(getRmax(Y1), std::max(getRmax(YA), getRmax(YB)));
324 }
325
326
327 const JNPEs_t& Y1; //!< PDF
328 const JNPEs_t& YA; //!< PDF
329 const JNPEs_t& YB; //!< PDF
330
331 std::shared_ptr<JMEstimator> estimator; //!< M-Estimator function
332
333 protected:
334 /**
335 * Get maximal road width of PDF.
336 *
337 * \param NPE NPE tables
338 * \return road width [m]
339 */
340 static inline double getRmax(const JNPEs_t& NPE)
341 {
342 double xmax = 0.0;
343
344 for (JNPEs_t::const_iterator i = NPE.begin(); i != NPE.end(); ++i) {
345 if (!i->empty() && i->getXmax() > xmax) {
346 xmax = i->getXmax();
347 }
348 }
349
350 return xmax;
351 }
352
353
354 /**
355 * Get number of photo-electrons.
356 *
357 * \param NPE NPE tables
358 * \param R distance between muon and PMT [m]
359 * \param theta zenith angle orientation PMT [rad]
360 * \param phi azimuth angle orientation PMT [rad]
361 * \return number of photo-electrons
362 */
363 static inline double getNPE(const JNPEs_t& NPE,
364 const double R,
365 const double theta,
366 const double phi)
367 {
368 using JTOOLS::get_value;
369
370 double npe = 0.0;
371
372 for (JNPEs_t::const_iterator i = NPE.begin(); i != NPE.end(); ++i) {
373
374 if (R <= i->getXmax()) {
375
376 try {
377
378 const double y = get_value((*i)(std::max(R, i->getXmin()), theta, phi));
379
380 if (y > 0.0) {
381 npe += y;
382 }
383 }
384 catch(const JLANG::JException& error) {
385 ERROR(error << std::endl);
386 }
387 }
388 }
389
390 return npe;
391 }
392 };
393}
394
395#endif
Auxiliary methods to evaluate Poisson probabilities and chi2.
Various implementations of functional maps.
Maximum likelihood estimator (M-estimators).
General purpose messaging.
#define ERROR(A)
Definition JMessage.hh:66
Auxiliary methods for PDF calculations.
Numbering scheme for PDF types.
Physics constants.
General purpose data regression method.
Data structure for fit of energy.
Definition JEnergy.hh:31
Axis object.
Definition JAxis3D.hh:41
double getY() const
Get y position.
Definition JVector3D.hh:104
double getZ() const
Get z position.
Definition JVector3D.hh:115
double getX() const
Get x position.
Definition JVector3D.hh:94
double getTheta() const
Get theta angle.
Definition JVersor3D.hh:128
double getPhi() const
Get phi angle.
Definition JVersor3D.hh:144
General exception.
Definition JException.hh:24
Custom class for integrated values of the PDF of the arrival time of Cherenkov light.
Definition JNPETable.hh:46
Multi-dimensional PDF table for arrival time of Cherenkov light.
Definition JPDFTable.hh:44
const double xmax
double getNPE(const Hit &hit)
Get true charge of hit.
Auxiliary classes and methods for linear and iterative data regression.
Definition JEnergy.hh:15
static const JZero zero
Function object to assign zero value.
Definition JZero.hh:105
JPDFType_t
PDF types.
Definition JPDFTypes.hh:24
@ SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition JPDFTypes.hh:33
@ DIRECT_LIGHT_FROM_EMSHOWERS
direct light from EM showers
Definition JPDFTypes.hh:29
@ SCATTERED_LIGHT_FROM_EMSHOWERS
scattered light from EM showers
Definition JPDFTypes.hh:30
@ SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition JPDFTypes.hh:27
@ DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition JPDFTypes.hh:32
@ DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition JPDFTypes.hh:26
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [ns]).
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
double getN(const JRange< T > &range, const double R)
Get expected number of occurrences due to given rate within specified interval.
Definition JRange.hh:704
JResultEvaluator< JResult_t >::result_type get_value(const JResult_t &value)
Helper method to recursively evaluate a to function value.
Definition JResult.hh:998
Abstract class for global fit method.
Definition JRegressor.hh:79
Normal M-estimator.
Auxiliary class for simultaneously handling light yields and response of PMT.
Definition JNPEHit.hh:21
double getChi2(const double E_GeV) const
Get chi2.
Definition JNPEHit.hh:87
Auxiliary class for handling various light yields.
Definition JNPE.hh:32
JRegressorStorage()
Default constructor.
JNPEs_t _YB
light from EM showers
JNPEs_t _YA
light from delta-rays
const JNPEs_t & getYB() const
Get light from EM showers.
std::vector< JNPE_t > JNPEs_t
NPEs.
JTOOLS::JMAPLIST< JTOOLS::JPolint1FunctionalMap, JTOOLS::JPolint1FunctionalGridMap, JTOOLS::JPolint1FunctionalGridMap >::maplist JNPEMaplist_t
JRegressorStorage(const std::string &fileDescriptor, const JTimeRange &T_ns=JTimeRange())
Constructor.
const JNPEs_t & getY1() const
Get light from muon.
JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
const JNPEs_t & getYA() const
Get light from delta-rays.
JPHYSICS::JNPETable< double, double, JNPEMaplist_t > JNPE_t
time integrated PDF
Template data structure for storage of internal data.
Template definition of a data regressor of given model.
Definition JRegressor.hh:70
Auxiliary class for recursive map list generation.
Definition JMapList.hh:109
Type definition of a 1st degree polynomial interpolation based on a JGridMap implementation.
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
Type definition of a spline interpolation method based on a JCollection with double result type.