Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
JLine3ZRegressor.hh
Go to the documentation of this file.
1#ifndef __JFIT__JLINE3ZREGRESSOR__
2#define __JFIT__JLINE3ZREGRESSOR__
3
4#include <string>
5#include <array>
6#include <memory>
7
10#include "JPhysics/JNPETable.hh"
12#include "JPhysics/JGeane.hh"
16#include "JTools/JRange.hh"
19#include "JMath/JZero.hh"
20#include "JFit/JTimeRange.hh"
21#include "JFit/JLine3Z.hh"
22#include "JFit/JPMTW0.hh"
23#include "JFit/JSimplex.hh"
24#include "JFit/JGandalf.hh"
25#include "JFit/JMEstimator.hh"
26#include "JFit/JRegressor.hh"
27#include "Jeep/JMessage.hh"
28
29
30/**
31 * \file
32 * Data regression method for JFIT::JLine3Z.
33 * \author mdejong
34 */
35namespace JFIT {}
36namespace JPP { using namespace JFIT; }
37
38namespace JFIT {
39
47
48
49 /**
50 * Regressor function object for JLine3Z fit using JSimplex minimiser.
51 */
52 template<>
53 struct JRegressor<JLine3Z, JSimplex> :
54 public JAbstractRegressor<JLine3Z, JSimplex>
55 {
56 using JAbstractRegressor<JLine3Z, JSimplex>::operator();
57
58 /**
59 * Constructor.
60 *
61 * \param sigma time resolution of hit [ns]
62 */
63 JRegressor(double sigma) :
64 estimator(new JMEstimatorNormal())
65 {
66 this->sigma = sigma;
67 }
68
69
70 /**
71 * Fit function.
72 * This method is used to determine the chi2 of given hit with respect to trajectory of muon.
73 *
74 * The template argument <tt>JHit_t</tt> refers to a data structure which should have the following member methods:
75 * - double getX(); // [m]
76 * - double getY(); // [m]
77 * - double getZ(); // [m]
78 * - double getT(); // [ns]
79 *
80 * \param track track
81 * \param hit hit
82 * \return chi2
83 */
84 template<class JHit_t>
85 double operator()(const JLine3Z& track, const JHit_t& hit) const
86 {
87 using namespace JPP;
88
89 JPosition3D D(hit.getX(), hit.getY(), hit.getZ());
90
91 D.sub(track.getPosition());
92
93 const double z = D.getDot(track.getDirection());
94 const double R = sqrt(D.getLengthSquared() - z*z);
95
96 const double t1 = track.getT() + (z + R * getKappaC()) * getInverseSpeedOfLight();
97
98 const double u = (t1 - hit.getT()) / sigma;
99
100 return estimator->getRho(u) * hit.getW();
101 }
102
103 std::shared_ptr<JMEstimator> estimator; //!< M-Estimator function
104 double sigma; //!< Time resolution [ns]
105 };
106
107
108 /**
109 * Template specialisation for storage of PDF tables.
110 */
111 template<>
113 {
119
124
126
127 static const int NUMBER_OF_PDFS = 6; //!< Number of PDFs
128
131
132
133 /**
134 * Default constructor.
135 */
138
139
140 /**
141 * Constructor.
142 *
143 * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
144 * will be replaced by the PDF types listed in JRegressorStorage<JLine3Z, JGandalf>::pdf_t.
145 *
146 * The <tt>TTS</tt> corresponds to the additional time smearing applied to the PDFs.
147 *
148 * \param fileDescriptor PDF file descriptor
149 * \param TTS TTS [ns]
150 * \param numberOfPoints number of points for Gauss-Hermite integration of TTS
151 * \param epsilon precision for Gauss-Hermite integration of TTS
152 */
153 JRegressorStorage(const std::string& fileDescriptor,
154 const double TTS,
155 const int numberOfPoints = 25,
156 const double epsilon = 1.0e-10)
157 {
158 using namespace std;
159 using namespace JPP;
160
161 const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
162
163 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
164
165 const string file_name = getFilename(fileDescriptor, pdf_t[i]);
166
167 _pdf[i].load(file_name.c_str());
168
169 _pdf[i].setExceptionHandler(supervisor);
170 }
171
172 // Add PDFs
173
174 for (int i = 1; i < NUMBER_OF_PDFS; i += 2) {
175
176 _pdf[ i ].add(_pdf[i-1]);
177
178 JPDF_t buffer;
179
180 _pdf[i-1].swap(buffer);
181
182 if (TTS > 0.0) {
183 _pdf[i].blur(TTS, numberOfPoints, epsilon);
184 }
185
186 _npe[ i ] = JNPE_t(_pdf[i]);
187 }
188 }
189
190
191 /**
192 * Get PDFs.
193 *
194 * \return PDFs
195 */
196 const JPDFs_t& getPDF() const
197 {
198 return _pdf;
199 }
200
201
202 /**
203 * Get NPEs.
204 *
205 * \return NPEs
206 */
207 const JNPEs_t& getNPE() const
208 {
209 return _npe;
210 }
211
212
213 /**
214 * Transform PDFs and NPEs.
215 *
216 * \param transformer transformer
217 */
218 void transform(const transformer_type& transformer)
219 {
220 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
221 _pdf[i].transform(transformer);
222 _npe[i].transform(transformer);
223 }
224 }
225
226
227 /**
228 * PDF types.
229 */
230 static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
231
232 private:
233 JPDFs_t _pdf; //!< PDFs
234 JNPEs_t _npe; //!< NPEs
235 };
236
237
238 /**
239 * PDF types.
240 */
241 const JPDFType_t JRegressorStorage<JLine3Z, JGandalf>::pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
242 SCATTERED_LIGHT_FROM_MUON,
243 DIRECT_LIGHT_FROM_DELTARAYS,
244 SCATTERED_LIGHT_FROM_DELTARAYS,
245 DIRECT_LIGHT_FROM_EMSHOWERS,
246 SCATTERED_LIGHT_FROM_EMSHOWERS };
247
248
249 /**
250 * Regressor function object for JLine3Z fit using JGandalf minimiser.
251 */
252 template<>
253 struct JRegressor<JLine3Z, JGandalf> :
254 public JAbstractRegressor<JLine3Z, JGandalf>,
255 public JRegressorStorage <JLine3Z, JGandalf>
256 {
257 using JAbstractRegressor<JLine3Z, JGandalf>::operator();
258
259 typedef JRegressorStorage<JLine3Z, JGandalf> storage_type;
260
261 /**
262 * Default constructor
263 */
264 JRegressor() :
265 storage_type(),
266 pdf(getPDF()),
267 npe(getNPE()),
268 estimator()
269 {}
270
271
272 /**
273 * Constructor.
274 *
275 * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
276 * will be replaced by the PDF types listed in JRegressorStorage<JLine3Z, JGandalf>::pdf_t.
277 *
278 * The <tt>TTS</tt> corresponds to the additional time smearing applied to the PDFs.
279 *
280 * \param fileDescriptor PDF file descriptor
281 * \param TTS TTS [ns]
282 * \param numberOfPoints number of points for Gauss-Hermite integration of TTS
283 * \param epsilon precision for Gauss-Hermite integration of TTS
284 */
285 JRegressor(const std::string& fileDescriptor,
286 const double TTS,
287 const int numberOfPoints = 25,
288 const double epsilon = 1.0e-10) :
289 storage_type(fileDescriptor, TTS, numberOfPoints, epsilon),
290 pdf(getPDF()),
291 npe(getNPE()),
292 estimator(new JMEstimatorNormal())
293 {}
294
295
296 /**
297 * Constructor.
298 *
299 * \param storage PDF storage
300 */
301 JRegressor(const storage_type& storage) :
302 pdf(storage.getPDF()),
303 npe(storage.getNPE()),
304 estimator(new JMEstimatorNormal())
305 {}
306
307
308 /**
309 * Fit function.
310 * This method is used to determine the chi2 and gradient of given hit with respect to trajectory of muon.
311 *
312 * The template argument <tt>JHit_t</tt> refers to a data structure which should have the following member methods:
313 * - double getX(); // [m]
314 * - double getY(); // [m]
315 * - double getZ(); // [m]
316 * - double getDX(); // [u]
317 * - double getDY(); // [u]
318 * - double getDZ(); // [u]
319 * - double getT(); // [ns]
320 *
321 * \param track track
322 * \param hit hit
323 * \return chi2 and gradient
324 */
325 template<class JHit_t>
326 result_type operator()(const JLine3Z& track, const JHit_t& hit) const
327 {
328 using namespace JPP;
329
330 JPosition3D D(hit.getX(), hit.getY(), hit.getZ());
331 JDirection3D U(hit.getDX(), hit.getDY(), hit.getDZ());
332
333 D.sub(track.getPosition());
334
335 const double z = D.getDot(track.getDirection());
336 const double x = D.getX() - z * track.getDX();
337 const double y = D.getY() - z * track.getDY();
338 const double R2 = D.getLengthSquared() - z*z;
339 const double R = (R2 > Rmin_m*Rmin_m ? sqrt(R2) : Rmin_m);
340
341 const double t1 = track.getT() + (z + R * getTanThetaC()) * getInverseSpeedOfLight();
342
343 U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
344
345 const double theta = U.getTheta();
346 const double phi = fabs(U.getPhi());
347
348 const double E = gWater.getE(E_GeV, z);
349 const double dt = T_ns.constrain(hit.getT() - t1);
350
351 JPDF_t::result_type H0 = getH0(hit.getR(), dt);
352 JPDF_t::result_type H1 = getH1(E, R, theta, phi, dt);
353
354 if (H1.V >= Vmax_npe) {
355 H1 *= Vmax_npe / H1.V;
356 }
357
358 H1 += H0; // signal + background
359
360 result_type result;
361
362 result.chi2 = H1.getChi2() - H0.getChi2(); // Likelihood ratio
363
364 const double wc = 1.0 - getTanThetaC() * z / R; // d(ct1)/d(z)
365
366 result.gradient = JLine3Z(JLine1Z(JPosition3D(-getTanThetaC() * D.getX() / R, // d(ct1)/d(x)
367 -getTanThetaC() * D.getY() / R, // d(ct1)/d(y)
368 0.0),
369 getSpeedOfLight()), // d(ct1)/d(t)
370 JVersor3Z(wc * (D.getX() - D.getZ()*track.getDX()/track.getDZ()), // d(ct1)/d(dx)
371 wc * (D.getY() - D.getZ()*track.getDY()/track.getDZ()))); // d(ct1)/d(dy)
372
373 result.gradient.mul(getInverseSpeedOfLight() * (H1.getDerivativeOfChi2() -
374 H0.getDerivativeOfChi2())); // x d(chi2)/d(ct1)
375
376 return result;
377 }
378
379
380 /**
381 * Fit function.
382 * This method is used to determine the chi2 and gradient of given PMT with respect to trajectory of muon.
383 *
384 * \param track track
385 * \param pmt pmt
386 * \return chi2 and gradient
387 */
388 result_type operator()(const JLine3Z& track, const JPMTW0& pmt) const
389 {
390 using namespace JGEOMETRY3D;
391 using namespace JPHYSICS;
392
393 JPosition3D D(pmt.getPosition());
394 JDirection3D U(pmt.getDirection());
395
396 D.sub(track.getPosition());
397
398 const double z = D.getDot(track.getDirection());
399 const double x = D.getX() - z * track.getDX();
400 const double y = D.getY() - z * track.getDY();
401 const double R2 = D.getLengthSquared() - z*z;
402 const double R = (R2 > Rmin_m*Rmin_m ? sqrt(R2) : Rmin_m);
403
404 U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
405
406 const double theta = U.getTheta();
407 const double phi = fabs(U.getPhi());
408
409 const double E = gWater.getE(E_GeV, z);
410
411 JNPE_t::result_type H0 = getH0(pmt.getR());
412 JNPE_t::result_type H1 = getH1(E, R, theta, phi);
413
414 if (H1.f >= Vmax_npe) {
415 H1 *= Vmax_npe / H1.f;
416 }
417
418 H1 += H0;
419
420 const bool hit = pmt.getN() != 0;
421 const double u = H1.getChi2(hit);
422
423 result_type result;
424
425 result.chi2 = estimator->getRho(u);
426
427 result.gradient = JLine3Z(JLine1Z(JPosition3D(-D.getX() / R, // d(R)/d(x)
428 -D.getY() / R, // d(R)/d(y)
429 0.0),
430 0.0),
431 JVersor3Z(-z * D.getX() / R, // d(R)/d(dx)
432 -z * D.getY() / R)); // d(R)/d(dy)
433
434 result.gradient.mul(estimator->getPsi(u));
435 result.gradient.mul(H1.getDerivativeOfChi2(hit)); // x d(chi2)/d(R)
436
437 return result;
438 }
439
440
441 /**
442 * Get background hypothesis value for time differentiated PDF.
443 *
444 * \param R_Hz rate [Hz]
445 * \param t1 time [ns]
446 * \return hypothesis value
447 */
448 JPDF_t::result_type getH0(const double R_Hz,
449 const double t1) const
450 {
451 return JPDF_t::result_type(R_Hz * 1e-9, t1, T_ns);
452 }
453
454
455 /**
456 * Get signal hypothesis value for time differentiated PDF.
457 *
458 * \param E muon energy at minimum distance of approach [GeV]
459 * \param R minimum distance of approach [m]
460 * \param theta PMT zenith angle [rad]
461 * \param phi PMT azimuth angle [rad]
462 * \param t1 arrival time relative to Cherenkov hypothesis [ns]
463 * \return hypothesis value
464 */
465 JPDF_t::result_type getH1(const double E,
466 const double R,
467 const double theta,
468 const double phi,
469 const double t1) const
470 {
471 using namespace std;
472 using namespace JPP;
473
474 JPDF_t::result_type h1 = zero;
475
476 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
477
478 if (!pdf[i].empty() && R <= pdf[i].getXmax()) {
479
480 JPDF_t::result_type y1 = pdf[i](max(R, pdf[i].getXmin()), theta, phi, t1);
481
482 // safety measures
483
484 if (y1.f <= 0.0) {
485 y1.f = 0.0;
486 y1.fp = 0.0;
487 }
488
489 if (y1.v <= 0.0) {
490 y1.v = 0.0;
491 }
492
493 // energy dependence
494
495 if (is_deltarays(pdf_t[i])) {
496 y1 *= getDeltaRaysFromMuon(E);
497 } else if (is_bremsstrahlung(pdf_t[i])) {
498 y1 *= E;
499 }
500
501 h1 += y1;
502 }
503 }
504
505 return h1;
506 }
507
508
509 /**
510 * Get background hypothesis value for time integrated PDF.
511 *
512 * \param R_Hz rate [Hz]
513 * \return hypothesis value
514 */
515 JNPE_t::result_type getH0(const double R_Hz) const
516 {
517 return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength(), 0.0);
518 }
519
520
521 /**
522 * Get signal hypothesis value for time integrated PDF.
523 *
524 * \param E muon energy at minimum distance of approach [GeV]
525 * \param R minimum distance of approach [m]
526 * \param theta PMT zenith angle [rad]
527 * \param phi PMT azimuth angle [rad]
528 * \return hypothesis value
529 */
530 JNPE_t::result_type getH1(const double E,
531 const double R,
532 const double theta,
533 const double phi) const
534 {
535 using namespace std;
536 using namespace JPP;
537
538 JNPE_t::result_type h1 = JMATH::zero;
539
540 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
541
542 if (!npe[i].empty() && R <= npe[i].getXmax()) {
543
544 try {
545
546 JNPE_t::result_type y1 = npe[i](max(R, npe[i].getXmin()), theta, phi);
547
548 // safety measures
549
550 if (y1.f > 0.0) {
551
552 // energy dependence
553
554 if (is_deltarays(pdf_t[i])) {
555 y1 *= getDeltaRaysFromMuon(E);
556 } else if (is_bremsstrahlung(pdf_t[i])) {
557 y1 *= E;
558 }
559
560 h1 += y1;
561 }
562 }
563 catch(JException& error) {
564 ERROR(error << endl);
565 }
566 }
567 }
568
569 return h1;
570 }
571
572
573 /**
574 * Get maximal road width of PDF.
575 *
576 * \return road width [m]
577 */
578 inline double getRmax() const
579 {
580 double xmax = 0.0;
581
582 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
583 if (!pdf[i].empty() && pdf[i].getXmax() > xmax) {
584 xmax = pdf[i].getXmax();
585 }
586 }
587
588 return xmax;
589 }
590
591
592 static JTimeRange T_ns; //!< Time window with respect to Cherenkov hypothesis [ns]
593 static double Vmax_npe; //!< Maximal integral of PDF [npe]
594 static double Rmin_m; //!< Minimal distance of [m]
595
596 const JPDFs_t& pdf; //!< PDF
597 const JNPEs_t& npe; //!< PDF
598
599 double E_GeV = 0.0; //!< Energy of muon at vertex [GeV]
600
601 std::shared_ptr<JMEstimator> estimator; //!< M-Estimator function
602 };
603
604
605 /**
606 * Default values.
607 */
609 double JRegressor<JLine3Z, JGandalf>::Vmax_npe = std::numeric_limits<double>::max();
611}
612
613#endif
Various implementations of functional maps.
Energy loss of muon.
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.
Auxiliary class to define a range between two values.
General purpose data regression method.
int numberOfPoints
Definition JResultPDF.cc:22
Definition of zero value for any class.
Fit method based on the Levenberg-Marquardt method.
Definition JGandalf.hh:87
Data structure for fit of straight line paralel to z-axis.
Definition JLine1Z.hh:29
Data structure for fit of straight line in positive z-direction.
Definition JLine3Z.hh:40
JVersor3D getDirection(const JVector3D &pos) const
Get photon direction of Cherenkov light on PMT.
Definition JLine3Z.hh:255
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition JLine3Z.hh:234
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++,...
Definition JSimplex.hh:44
Data structure for direction in three dimensions.
const JDirection3D & getDirection() const
Get direction.
Data structure for position in three dimensions.
const JPosition3D & getPosition() const
Get position.
Rotation around Z-axis.
Data structure for normalised vector in positive z-direction.
Definition JVersor3Z.hh:41
double getDZ() const
Get z direction.
Definition JVersor3Z.hh:169
double getDY() const
Get y direction.
Definition JVersor3Z.hh:158
double getDX() const
Get x direction.
Definition JVersor3Z.hh:147
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
Interface for weight application and coordinate transformation of function.
const double sigma[]
const double xmax
const double epsilon
double getNPE(const Hit &hit)
Get true charge of hit.
Auxiliary classes and methods for linear and iterative data regression.
Definition JEnergy.hh:15
Auxiliary classes and methods for 3D geometrical objects and operations.
Definition JAngle3D.hh:19
static const JZero zero
Function object to assign zero value.
Definition JZero.hh:105
Auxiliary methods for light properties of deep-sea water.
double getDeltaRaysFromMuon(const double E, const JRange< double > T_GeV=JRange< double >(DELTARAY_TMIN, DELTARAY_TMAX))
Equivalent EM-shower energy due to delta-rays per unit muon track length.
double getKappaC()
Get average R-dependence of arrival time of Cherenkov light (a.k.a.
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition JPDFTypes.hh:151
bool is_bremsstrahlung(const int pdf)
Test if given PDF type corresponds to Cherenkov light from Bremsstrahlung.
Definition JPDFTypes.hh:137
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
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Abstract class for global fit method.
Definition JRegressor.hh:79
Normal M-estimator.
Auxiliary class for handling PMT geometry, rate and response.
Definition JPMTW0.hh:24
int getN() const
Get number of hits.
Definition JPMTW0.hh:67
double getR() const
Get rate.
Definition JPMTW0.hh:56
JPHYSICS::JNPETable< double, double, JNPEMaplist_t > JNPE_t
time integrated PDF
const JNPEs_t & getNPE() const
Get NPEs.
std::array< JPDF_t, NUMBER_OF_PDFS > JPDFs_t
PDFs.
JPHYSICS::JPDFTable< JFunction1D_t, JPDFMaplist_t > JPDF_t
time dependent PDF
JTOOLS::JMAPLIST< JTOOLS::JPolint1FunctionalMap, JTOOLS::JPolint0FunctionalGridMap, JTOOLS::JPolint0FunctionalGridMap >::maplist JPDFMaplist_t
std::array< JNPE_t, NUMBER_OF_PDFS > JNPEs_t
NPEs.
JTOOLS::JMAPLIST< JTOOLS::JPolint1FunctionalMapH, JTOOLS::JPolint1FunctionalGridMap, JTOOLS::JPolint1FunctionalGridMap >::maplist JNPEMaplist_t
JRegressorStorage(const std::string &fileDescriptor, const double TTS, const int numberOfPoints=25, const double epsilon=1.0e-10)
Constructor.
const JPDFs_t & getPDF() const
Get PDFs.
void transform(const transformer_type &transformer)
Transform PDFs and NPEs.
Template data structure for storage of internal data.
Template definition of a data regressor of given model.
Definition JRegressor.hh:70
Auxiliary class to set-up Hit.
Definition JSirene.hh:58
Auxiliary class for recursive map list generation.
Definition JMapList.hh:109
Type definition of a zero degree polynomial interpolation based on a JGridMap implementation.
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 1st degree polynomial interpolation based on a JMap implementation.
Type definition of a spline interpolation method based on a JCollection with JResultPDF result type.