1#ifndef __JFIT__JLINE3ZREGRESSOR__
2#define __JFIT__JLINE3ZREGRESSOR__
36namespace JPP {
using namespace JFIT; }
84 template<
class JHit_t>
85 double operator()(
const JLine3Z& track,
const JHit_t& hit)
const
94 const double R = sqrt(D.getLengthSquared() - z*z);
96 const double t1 = track.
getT() + (z + R *
getKappaC()) * getInverseSpeedOfLight();
98 const double u = (t1 - hit.getT()) /
sigma;
100 return estimator->getRho(u) * hit.getW();
103 std::shared_ptr<JMEstimator> estimator;
127 static const int NUMBER_OF_PDFS = 6;
156 const double epsilon = 1.0e-10)
161 const JPDF_t::JSupervisor supervisor(
new JPDF_t::JDefaultResult(
JMATH::zero));
163 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
165 const string file_name = getFilename(fileDescriptor, pdf_t[i]);
167 _pdf[i].load(file_name.c_str());
169 _pdf[i].setExceptionHandler(supervisor);
174 for (
int i = 1; i < NUMBER_OF_PDFS; i += 2) {
176 _pdf[ i ].add(_pdf[i-1]);
180 _pdf[i-1].swap(buffer);
186 _npe[ i ] =
JNPE_t(_pdf[i]);
220 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
221 _pdf[i].transform(transformer);
222 _npe[i].transform(transformer);
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 };
288 const double epsilon = 1.0e-10) :
302 pdf(storage.getPDF()),
325 template<
class JHit_t>
326 result_type operator()(
const JLine3Z& track,
const JHit_t& hit)
const
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);
341 const double t1 = track.
getT() + (z + R * getTanThetaC()) * getInverseSpeedOfLight();
345 const double theta = U.getTheta();
346 const double phi = fabs(U.getPhi());
348 const double E = gWater.getE(E_GeV, z);
349 const double dt = T_ns.constrain(hit.getT() - t1);
351 JPDF_t::result_type H0 = getH0(hit.getR(), dt);
352 JPDF_t::result_type H1 = getH1(E, R, theta, phi, dt);
354 if (H1.V >= Vmax_npe) {
355 H1 *= Vmax_npe / H1.V;
362 result.chi2 = H1.getChi2() - H0.getChi2();
364 const double wc = 1.0 - getTanThetaC() * z / R;
367 -getTanThetaC() * D.getY() / R,
371 wc * (D.getY() - D.getZ()*track.
getDY()/track.
getDZ())));
373 result.gradient.mul(getInverseSpeedOfLight() * (H1.getDerivativeOfChi2() -
374 H0.getDerivativeOfChi2()));
388 result_type operator()(
const JLine3Z& track,
const JPMTW0& pmt)
const
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);
406 const double theta = U.getTheta();
407 const double phi = fabs(U.getPhi());
409 const double E = gWater.getE(E_GeV, z);
411 JNPE_t::result_type H0 = getH0(pmt.
getR());
412 JNPE_t::result_type H1 = getH1(E, R, theta, phi);
414 if (H1.f >= Vmax_npe) {
415 H1 *= Vmax_npe / H1.f;
420 const bool hit = pmt.
getN() != 0;
421 const double u = H1.getChi2(hit);
425 result.chi2 = estimator->getRho(u);
434 result.gradient.mul(estimator->getPsi(u));
435 result.gradient.mul(H1.getDerivativeOfChi2(hit));
448 JPDF_t::result_type getH0(
const double R_Hz,
449 const double t1)
const
451 return JPDF_t::result_type(R_Hz * 1e-9, t1, T_ns);
465 JPDF_t::result_type getH1(
const double E,
469 const double t1)
const
474 JPDF_t::result_type h1 = zero;
476 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
478 if (!pdf[i].empty() && R <= pdf[i].getXmax()) {
480 JPDF_t::result_type y1 = pdf[i](max(R, pdf[i].getXmin()), theta, phi, t1);
515 JNPE_t::result_type getH0(
const double R_Hz)
const
517 return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength(), 0.0);
530 JNPE_t::result_type getH1(
const double E,
533 const double phi)
const
540 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
542 if (!npe[i].empty() && R <= npe[i].getXmax()) {
546 JNPE_t::result_type y1 = npe[i](max(R, npe[i].getXmin()), theta, phi);
564 ERROR(error << endl);
578 inline double getRmax()
const
582 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
583 if (!pdf[i].empty() && pdf[i].getXmax() > xmax) {
584 xmax = pdf[i].getXmax();
593 static double Vmax_npe;
594 static double Rmin_m;
601 std::shared_ptr<JMEstimator> estimator;
Various implementations of functional maps.
Maximum likelihood estimator (M-estimators).
General purpose messaging.
Numbering scheme for PDF types.
Auxiliary class to define a range between two values.
General purpose data regression method.
Definition of zero value for any class.
Fit method based on the Levenberg-Marquardt method.
Data structure for fit of straight line paralel to z-axis.
Data structure for fit of straight line in positive z-direction.
JVersor3D getDirection(const JVector3D &pos) const
Get photon direction of Cherenkov light on PMT.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++,...
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.
Data structure for normalised vector in positive z-direction.
double getDZ() const
Get z direction.
double getDY() const
Get y direction.
double getDX() const
Get x direction.
Custom class for integrated values of the PDF of the arrival time of Cherenkov light.
Multi-dimensional PDF table for arrival time of Cherenkov light.
double getNPE(const Hit &hit)
Get true charge of hit.
Auxiliary classes and methods for linear and iterative data regression.
Auxiliary classes and methods for 3D geometrical objects and operations.
static const JZero zero
Function object to assign zero value.
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.
bool is_bremsstrahlung(const int pdf)
Test if given PDF type corresponds to Cherenkov light from Bremsstrahlung.
@ SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
@ DIRECT_LIGHT_FROM_EMSHOWERS
direct light from EM showers
@ SCATTERED_LIGHT_FROM_EMSHOWERS
scattered light from EM showers
@ SCATTERED_LIGHT_FROM_MUON
scattered light from muon
@ DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
@ DIRECT_LIGHT_FROM_MUON
direct light from muon
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Abstract class for global fit method.
Auxiliary class for handling PMT geometry, rate and response.
int getN() const
Get number of hits.
double getR() const
Get rate.
JPHYSICS::JNPETable< double, double, JNPEMaplist_t > JNPE_t
time integrated PDF
const JNPEs_t & getNPE() const
Get NPEs.
JRegressorStorage()
Default constructor.
JTOOLS::JSplineFunction1S_t JFunction1D_t
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
JPDF_t::transformer_type transformer_type
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.
Auxiliary class to set-up Hit.