1 #ifndef __JFIT__JLINE3ZREGRESSOR__
2 #define __JFIT__JLINE3ZREGRESSOR__
35 namespace JPP {
using namespace JFIT; }
83 template<
class JHit_t>
93 const double R = sqrt(
D.getLengthSquared() - z*z);
97 const double u = (t1 - hit.getT()) / sigma;
99 return estimator->getRho(
u) * hit.getW();
150 const double epsilon = 1.0e-10) :
157 const JPDF_t::JSupervisor supervisor(
new JPDF_t::JDefaultResult(
JMATH::zero));
159 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
163 const string file_name =
getFilename(fileDescriptor, pdf_t[i]);
165 NOTICE(
"loading PDF from file " << file_name <<
"... " << flush);
167 pdf[i].load(file_name.c_str());
171 pdf[i].setExceptionHandler(supervisor);
180 for (
int i = 1; i < NUMBER_OF_PDFS; i += 2) {
182 pdf[ i ].add(pdf[i-1]);
186 pdf[i-1].swap(buffer);
190 }
else if (TTS < 0.0) {
191 ERROR(
"Illegal value of TTS [ns]: " << TTS << endl);
194 npe[ i ] =
JNPE_t(pdf[i]);
216 template<
class JHit_t>
222 JDirection3D U(hit.getDX(), hit.getDY(), hit.getDZ());
227 const double x =
D.getX() - z * track.
getDX();
228 const double y =
D.getY() - z * track.
getDY();
229 const double R2 =
D.getLengthSquared() - z*z;
230 const double R = (R2 > Rmin_m*Rmin_m ? sqrt(R2) : Rmin_m);
234 U.rotate(JRotation3Z(-atan2(y,x)));
236 const double theta = U.getTheta();
237 const double phi = fabs(U.getPhi());
240 const double dt = T_ns.constrain(hit.getT() - t1);
245 if (H1.V >= Vmax_npe) {
246 H1 *= Vmax_npe / H1.V;
253 result.chi2 = H1.getChi2() - H0.getChi2();
262 wc * (
D.getY() -
D.getZ()*track.
getDY()/track.
getDZ())));
265 H0.getDerivativeOfChi2()));
281 using namespace JGEOMETRY3D;
282 using namespace JPHYSICS;
290 const double x =
D.getX() - z * track.
getDX();
291 const double y =
D.getY() - z * track.
getDY();
292 const double R2 =
D.getLengthSquared() - z*z;
293 const double R = (R2 > Rmin_m*Rmin_m ? sqrt(R2) : Rmin_m);
295 U.rotate(JRotation3Z(-atan2(y,x)));
297 const double theta = U.getTheta();
298 const double phi = fabs(U.getPhi());
305 if (H1.f >= Vmax_npe) {
306 H1 *= Vmax_npe / H1.f;
311 const bool hit = pmt.
getN() != 0;
312 const double u = H1.getChi2(hit);
316 result.chi2 = estimator->getRho(
u);
325 result.gradient.mul(estimator->getPsi(
u));
326 result.gradient.mul(H1.getDerivativeOfChi2(hit));
340 const double t1)
const
360 const double t1)
const
367 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
369 if (!pdf[i].empty() && R <= pdf[i].getXmax()) {
424 const double phi)
const
431 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
433 if (!npe[i].empty() && R <= npe[i].getXmax()) {
455 ERROR(error << endl);
473 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
474 if (!pdf[i].empty() && pdf[i].getXmax() > xmax) {
475 xmax = pdf[i].getXmax();
490 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
491 if (!pdf[i].empty()) {
502 static const int NUMBER_OF_PDFS = 4;
Template definition of a data regressor of given model.
General purpose data regression method.
do echo Generating $dir eval D
double getR() const
Get rate.
JTOOLS::JSplineFunction1S_t JFunction1D_t
JPDF_t::result_type getH0(const double R_Hz, const double t1) const
Get background hypothesis value for time differentiated PDF.
JPDF_t::result_type getH1(const double E, const double R, const double theta, const double phi, const double t1) const
Get signal hypothesis value for time differentiated PDF.
JNPE_t::result_type getH1(const double E, const double R, const double theta, const double phi) const
Get signal hypothesis value for time integrated PDF.
JRegressor(const std::string &fileDescriptor, const double TTS, const int numberOfPoints=25, const double epsilon=1.0e-10)
Constructor.
Auxiliary class for handling PMT geometry, rate and response.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
const JDirection3D & getDirection() const
Get direction.
JRegressor(double sigma)
Constructor.
direct light from EM showers
void setRmax(const double Rmax)
Set maximal road width of PDF.
void blur(const double TTS, const int numberOfPoints=25, const double epsilon=1.0e-10, const double quantile=0.99)
Blur PDF.
then for HISTOGRAM in h0 h1
Data structure for fit of straight line in positive z-direction.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
static const JZero zero
Function object to assign zero value.
JNPE_t::result_type getH0(const double R_Hz) const
Get background hypothesis value for time integrated PDF.
Definition of zero value for any class.
virtual double getE(const double E, const double dx) const
Get energy of muon after specified distance.
static double Rmin_m
Minimal distance of [m].
double getDY() const
Get y direction.
Various implementations of functional maps.
esac $JPP_DIR examples JDetector JTransitTime o $OUTPUT_FILE n N $NPE T $TTS_NS d $DEBUG for HISTOGRAM in tts tt2 pmt
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
Numbering scheme for PDF types.
result_type operator()(const JLine3Z &track, const JHit_t &hit) const
Fit function.
bool is_bremsstrahlung(const int pdf)
Test if given PDF type corresponds to Cherenkov light from Bremsstrahlung.
JVersor3D getDirection(const JVector3D &pos) const
Get photon direction of Cherenkov light on PMT.
scattered light from muon
minimiser_type::result_type result_type
The template JSharedPointer class can be used to share a pointer to an object.
JTOOLS::JMAPLIST< JTOOLS::JPolint1FunctionalMapH, JTOOLS::JPolint1FunctionalGridMap, JTOOLS::JPolint1FunctionalGridMap >::maplist JNPEMaplist_t
multifunction_t::result_type result_type
double E_GeV
Energy of muon at vertex [GeV].
scattered light from delta-rays
scattered light from EM showers
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
const JPosition3D & getPosition() const
Get position.
double operator()(const JLine3Z &track, const JHit_t &hit) const
Fit function.
result_type operator()(const JLine3Z &track, const JPMTW0 &pmt) const
Fit function.
JPHYSICS::JNPETable< double, double, JNPEMaplist_t > JNPE_t
time integrated PDF
then usage $script[distance] fi case set_variable R
General purpose messaging.
double getRmax() const
Get maximal road width of PDF.
Fit method based on the Levenberg-Marquardt method.
direct light from delta-rays
double sigma
Time resolution [ns].
int getN() const
Get number of hits.
Auxiliary class to define a range between two values.
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++...
static double Vmax_npe
Maximal integral of PDF [npe].
double getDX() const
Get x direction.
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Data structure for fit of straight line paralel to z-axis.
double getDeltaRaysFromMuon(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
std::string getFilename(const std::string &file_name)
Get file name part, i.e. part after last JEEP::PATHNAME_SEPARATOR if any.
JPHYSICS::JPDFTable< JFunction1D_t, JPDFMaplist_t > JPDF_t
time dependent PDF
Data structure for position in three dimensions.
transformablemultifunction_type::result_type result_type
double getDZ() const
Get z direction.
virtual const char * what() const
Get error message.
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
Data structure for normalised vector in positive z-direction.
JRegressor()
Default constructor.
JTOOLS::JMAPLIST< JTOOLS::JPolint1FunctionalMap, JTOOLS::JPolint0FunctionalGridMap, JTOOLS::JPolint0FunctionalGridMap >::maplist JPDFMaplist_t
Abstract class for global fit method.
Maximum likelihood estimator (M-estimators).
then usage $script[input file[working directory[option]]] nWhere option can be E