Jpp 19.3.0-rc.1
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 range time range [ns]
78 */
79 JRegressorStorage(const std::string& fileDescriptor, const JTimeRange& range = JTimeRange())
80 {
81 using namespace std;
82 using namespace JPP;
83
84 typedef JSplineFunction1D_t JFunction1D_t;
86
87 const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
88
89 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
90
91 JPDF_t pdf;
92
93 const JPDFType_t type = pdf_t[i];
94 const string file_name = getFilename(fileDescriptor, type);
95
96 pdf.load(file_name.c_str());
97
98 pdf.setExceptionHandler(supervisor);
99
100 if (is_bremsstrahlung(type))
101 _YB.push_back(JNPE_t(pdf, range));
102 else if (is_deltarays(type))
103 _YA.push_back(JNPE_t(pdf, range));
104 else
105 _Y1.push_back(JNPE_t(pdf, range));
106 }
107
108 // Add PDFs
109
110 if (_Y1.size() == 2) { _Y1[1].add(_Y1[0]); _Y1.erase(_Y1.begin()); }
111 if (_YA.size() == 2) { _YA[1].add(_YA[0]); _YA.erase(_YA.begin()); }
112 if (_YB.size() == 2) { _YB[1].add(_YB[0]); _YB.erase(_YB.begin()); }
113 }
114
115
116 /**
117 * Get light from muon.
118 *
119 * \return NPEs
120 */
121 const JNPEs_t& getY1() const
122 {
123 return _Y1;
124 }
125
126
127 /**
128 * Get light from delta-rays.
129 *
130 * \return NPEs
131 */
132 const JNPEs_t& getYA() const
133 {
134 return _YA;
135 }
136
137
138 /**
139 * Get light from EM showers.
140 *
141 * \return NPEs
142 */
143 const JNPEs_t& getYB() const
144 {
145 return _YB;
146 }
147
148
149 /**
150 * PDF types.
151 */
152 static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
153
154 private:
155 JNPEs_t _Y1; //!< light from muon
156 JNPEs_t _YA; //!< light from delta-rays
157 JNPEs_t _YB; //!< light from EM showers
158 };
159
160
161 /**
162 * PDF types.
163 */
164 const JPDFType_t JRegressorStorage<JEnergy>::pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
165 SCATTERED_LIGHT_FROM_MUON,
166 DIRECT_LIGHT_FROM_DELTARAYS,
167 SCATTERED_LIGHT_FROM_DELTARAYS,
168 DIRECT_LIGHT_FROM_EMSHOWERS,
169 SCATTERED_LIGHT_FROM_EMSHOWERS };
170
171
172
173 /**
174 * Regressor function object for fit of muon energy.
175 */
176 template<>
177 struct JRegressor<JEnergy> :
178 public JAbstractRegressor<JEnergy>,
179 public JRegressorStorage <JEnergy>
180 {
181 using JAbstractRegressor<JEnergy>::operator();
182
183 typedef JRegressorStorage<JEnergy> storage_type;
184
185 using storage_type::getY1;
186 using storage_type::getYA;
187 using storage_type::getYB;
188
189
190 /**
191 * Constructor.
192 *
193 * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
194 * will be replaced by the PDF types listed in JRegressor<JEnergy>::pdf_t.
195 *
196 * \param fileDescriptor PDF file descriptor
197 * \param range time range [ns]
198 */
199 JRegressor(const std::string& fileDescriptor, const JTimeRange& range = JTimeRange()) :
200 storage_type(fileDescriptor, range),
201 Y1(getY1()),
202 YA(getYA()),
203 YB(getYB()),
204 estimator(new JMEstimatorNormal())
205 {}
206
207
208 /**
209 * Constructor.
210 *
211 * \param storage PDF storage
212 */
213 JRegressor(const storage_type& storage) :
214 Y1(storage.getY1()),
215 YA(storage.getYA()),
216 YB(storage.getYB()),
217 estimator(new JMEstimatorNormal())
218 {}
219
220
221 /**
222 * Fit function.
223 * This method is used to determine chi2 of given number of photo-electrons for given energy of muon.
224 *
225 * \param x energy
226 * \param npe npe
227 * \return chi2
228 */
229 inline double operator()(const JEnergy& x, const JNPEHit& npe) const
230 {
231 const double E = x.getE();
232 const double u = npe.getChi2(E);
233
234 return estimator->getRho(u);
235 }
236
237
238 /**
239 * Get number of photo-electrons from muon.
240 * Note that the PMT geometry should be relative to the muon trajectory,
241 * conform method JGEOMETRY3D::JAxis3D::transform.
242 *
243 * \param axis PMT axis
244 * \return npe
245 */
246 inline double getY1(const JAxis3D& axis) const
247 {
248 return getNPE(Y1, hypot(axis.getX(), axis.getY()), axis.getTheta(), fabs(axis.getPhi()));
249 }
250
251
252 /**
253 * Get number of photo-electrons from delta-rays.
254 * Note that the PMT geometry should be relative to the muon trajectory,
255 * conform method JGEOMETRY3D::JAxis3D::transform.
256 *
257 * \param axis PMT axis
258 * \return npe
259 */
260 inline double getYA(const JAxis3D& axis) const
261 {
262 return getNPE(YA, hypot(axis.getX(), axis.getY()), axis.getTheta(), fabs(axis.getPhi()));
263 }
264
265
266 /**
267 * Get number of photo-electrons from Bremsstrahlung.
268 * Note that the PMT geometry should be relative to the muon trajectory,
269 * conform method JGEOMETRY3D::JAxis3D::transform.
270 *
271 * \param axis PMT axis
272 * \return npe
273 */
274 inline double getYB(const JAxis3D& axis) const
275 {
276 return getNPE(YB, hypot(axis.getX(), axis.getY()), axis.getTheta(), fabs(axis.getPhi()));
277 }
278
279
280 /**
281 * Create data structure for handling light yields for PMT.
282 * Note that the PMT geometry should be relative to the muon trajectory,
283 * conform method JGEOMETRY3D::JAxis3D::transform.
284 *
285 * \param axis PMT axis
286 * \param R_Hz singles rate [Hz]
287 * \return light yields
288 */
289 inline JNPE getNPE(const JAxis3D& axis,
290 const double R_Hz) const
291 {
292 using namespace std;
293 using namespace JPP;
294
295 const double x = axis.getX();
296 const double y = axis.getY();
297 const double R = sqrt(x*x + y*y);
298 const double z = axis.getZ() - R / getTanThetaC();
299
300 const double theta = axis.getTheta();
301 const double phi = fabs(axis.getPhi());
302
303 const double y1 = getNPE(Y1, R, theta, phi);
304 const double yA = getNPE(YA, R, theta, phi);
305 const double yB = getNPE(YB, R, theta, phi);
306
307 return JNPE(getN(T_ns, R_Hz * 1.0e-9), y1, yA, yB, z);
308 }
309
310
311 /**
312 * Get maximal road width of NPE.
313 *
314 * \return road width [m]
315 */
316 inline double getRmax() const
317 {
318 return std::max(getRmax(Y1), std::max(getRmax(YA), getRmax(YB)));
319 }
320
321
322 static JTimeRange T_ns; //!< Time window with respect to Cherenkov hypothesis [ns]
323
324 const JNPEs_t& Y1; //!< PDF
325 const JNPEs_t& YA; //!< PDF
326 const JNPEs_t& YB; //!< PDF
327
328 std::shared_ptr<JMEstimator> estimator; //!< M-Estimator function
329
330 protected:
331 /**
332 * Get maximal road width of PDF.
333 *
334 * \param NPE NPE tables
335 * \return road width [m]
336 */
337 static inline double getRmax(const JNPEs_t& NPE)
338 {
339 double xmax = 0.0;
340
341 for (JNPEs_t::const_iterator i = NPE.begin(); i != NPE.end(); ++i) {
342 if (!i->empty() && i->getXmax() > xmax) {
343 xmax = i->getXmax();
344 }
345 }
346
347 return xmax;
348 }
349
350
351 /**
352 * Get number of photo-electrons.
353 *
354 * \param NPE NPE tables
355 * \param R distance between muon and PMT [m]
356 * \param theta zenith angle orientation PMT [rad]
357 * \param phi azimuth angle orientation PMT [rad]
358 * \return number of photo-electrons
359 */
360 static inline double getNPE(const JNPEs_t& NPE,
361 const double R,
362 const double theta,
363 const double phi)
364 {
365 using JTOOLS::get_value;
366
367 double npe = 0.0;
368
369 for (JNPEs_t::const_iterator i = NPE.begin(); i != NPE.end(); ++i) {
370
371 if (R <= i->getXmax()) {
372
373 try {
374
375 const double y = get_value((*i)(std::max(R, i->getXmin()), theta, phi));
376
377 if (y > 0.0) {
378 npe += y;
379 }
380 }
381 catch(const JLANG::JException& error) {
382 ERROR(error << std::endl);
383 }
384 }
385 }
386
387 return npe;
388 }
389 };
390
391
392 /**
393 * Time range.
394 */
396}
397
398#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
JRegressorStorage(const std::string &fileDescriptor, const JTimeRange &range=JTimeRange())
Constructor.
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
const JNPEs_t & getY1() const
Get light from muon.
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.