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;
158 const double epsilon = 1.0e-10) :
164 const JPDF_t::JSupervisor supervisor(
new JPDF_t::JDefaultResult(
JMATH::zero));
166 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
168 const string file_name = getFilename(fileDescriptor, pdf_t[i]);
170 _pdf[i].load(file_name.c_str());
172 _pdf[i].setExceptionHandler(supervisor);
177 for (
int i = 1; i < NUMBER_OF_PDFS; i += 2) {
179 _pdf[ i ].add(_pdf[i-1]);
183 _pdf[i-1].swap(buffer);
189 _npe[ i ] =
JNPE_t(_pdf[i]);
223 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
224 _pdf[i].transform(transformer);
225 _npe[i].transform(transformer);
247 SCATTERED_LIGHT_FROM_MUON,
248 DIRECT_LIGHT_FROM_DELTARAYS,
249 SCATTERED_LIGHT_FROM_DELTARAYS,
250 DIRECT_LIGHT_FROM_EMSHOWERS,
251 SCATTERED_LIGHT_FROM_EMSHOWERS };
295 const double epsilon = 1.0e-10) :
309 pdf(storage.getPDF()),
334 template<
class JHit_t>
335 result_type operator()(
const JLine3Z& track,
const JHit_t& hit)
const
345 const double x = D.getX() - z * track.
getDX();
346 const double y = D.getY() - z * track.
getDY();
347 const double R2 = D.getLengthSquared() - z*z;
348 const double R = (R2 > Rmin_m*Rmin_m ? sqrt(R2) : Rmin_m);
350 const double t1 = track.
getT() + (z + R * getTanThetaC()) * getInverseSpeedOfLight();
354 const double theta = U.getTheta();
355 const double phi = fabs(U.getPhi());
357 const double E = gWater.getE(E_GeV, z);
358 const double dt = T_ns.constrain(hit.getT() - t1);
360 JPDF_t::result_type H0 = getH0(hit.getR(), dt);
361 JPDF_t::result_type H1 = getH1(E, R, theta, phi, dt);
363 if (H1.V >= Vmax_npe) {
364 H1 *= Vmax_npe / H1.V;
371 result.chi2 = H1.getChi2() - H0.getChi2();
373 const double wc = 1.0 - getTanThetaC() * z / R;
376 -getTanThetaC() * D.getY() / R,
380 wc * (D.getY() - D.getZ()*track.
getDY()/track.
getDZ())));
382 result.gradient.mul(getInverseSpeedOfLight() * (H1.getDerivativeOfChi2() -
383 H0.getDerivativeOfChi2()));
397 result_type operator()(
const JLine3Z& track,
const JPMTW0& pmt)
const
408 const double x = D.getX() - z * track.
getDX();
409 const double y = D.getY() - z * track.
getDY();
410 const double R2 = D.getLengthSquared() - z*z;
411 const double R = (R2 > Rmin_m*Rmin_m ? sqrt(R2) : Rmin_m);
415 const double theta = U.getTheta();
416 const double phi = fabs(U.getPhi());
418 const double E = gWater.getE(E_GeV, z);
420 JNPE_t::result_type H0 = getH0(pmt.
getR());
421 JNPE_t::result_type H1 = getH1(E, R, theta, phi);
423 if (H1.f >= Vmax_npe) {
424 H1 *= Vmax_npe / H1.f;
429 const bool hit = pmt.
getN() != 0;
430 const double u = H1.getChi2(hit);
434 result.chi2 = estimator->getRho(u);
443 result.gradient.mul(estimator->getPsi(u));
444 result.gradient.mul(H1.getDerivativeOfChi2(hit));
457 JPDF_t::result_type getH0(
const double R_Hz,
458 const double t1)
const
460 return JPDF_t::result_type(R_Hz * 1e-9, t1, T_ns);
474 JPDF_t::result_type getH1(
const double E,
478 const double t1)
const
483 JPDF_t::result_type h1 = zero;
485 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
487 if (!pdf[i].empty() && R <= pdf[i].getXmax()) {
489 JPDF_t::result_type y1 = pdf[i](max(R, pdf[i].getXmin()), theta, phi, t1);
524 JNPE_t::result_type getH0(
const double R_Hz)
const
526 return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength(), 0.0);
539 JNPE_t::result_type getH1(
const double E,
542 const double phi)
const
549 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
551 if (!npe[i].empty() && R <= npe[i].getXmax()) {
555 JNPE_t::result_type y1 = npe[i](max(R, npe[i].getXmin()), theta, phi);
573 ERROR(error << endl);
587 inline double getRmax()
const
591 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
592 if (!pdf[i].empty() && pdf[i].getXmax() > xmax) {
593 xmax = pdf[i].getXmax();
601 static double Vmax_npe;
602 static double Rmin_m;
609 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
JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
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
const JPDFs_t & getPDF() const
Get PDFs.
JRegressorStorage(const std::string &fileDescriptor, const JTimeRange &T_ns, const double TTS_ns, const int numberOfPoints=25, const double epsilon=1.0e-10)
Constructor.
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.