1#ifndef __JRECONSTRUCTION__JMUONENERGY__
2#define __JRECONSTRUCTION__JMUONENERGY__
75 using JRegressor_t::operator();
119 const storage_type& storage,
121 const int debug = 0):
132 JRegressor_t::debug =
debug;
160 const JBuildL0 <JHitR0> buildL0;
167 for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
169 if (router.
hasModule(i->getModuleID())) {
171 buffer(*i, router.
getModule(i->getModuleID()));
173 buildL0(buffer, back_inserter(data[i->getModuleID()]));
178 if (!module.empty()) {
198 const double chi2 = std::numeric_limits<double>::max())
209 operator bool()
const
211 return chi2 != std::numeric_limits<double>::max();
244 for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
249 double zmin = numeric_limits<double>::lowest();
257 for (
const auto& module : input.
data) {
270 for (
size_t i = 0; i !=
module->size(); ++i) {
272 if (module.getStatus(i)) {
278 return (hit.
getPMT() == pmt && T_ns(hit.
getT()));
284 } match = { JRegressor_t::T_ns + t1, i };
286 JPMT pmt =
module->getPMT(i);
290 const JNPEHit hit(this->
getNPE(pmt, module.frame.getRate(i)), count_if(module.begin(), module.end(), match));
292 DEBUG(
"hit: " << setw(8) << module->getID() <<
'.' <<
FILL(2,
'0') << i <<
' '
294 <<
FIXED(7,3) << module.frame.getRate(i) * 1.0e-3 <<
' '
299 << setw(2) << hit.
getN() << endl);
308 const int NDF =
distance(data.begin(), data.end()) - 1;
318 for (
int i = 0; i != N; ++i) {
328 for (
int i = 0; i != N; ++i) {
333 const double chi2 = (*this)(x, data.begin(), data.end());
336 buffer[chi2] = x.getE();
345 for (
int i = 0; i != N; ++i) {
391 result[2].chi2 = (*this)(
result[2].x, data.begin(), data.end());
395 const double chi2 =
result[2].chi2;
396 const double E =
result[2].x.getE();
400 double Emin = numeric_limits<double>::max();
401 double Emax = numeric_limits<double>::lowest();
404 if (i->second < Emin) { Emin = i->second; }
405 if (i->second > Emax) { Emax = i->second; }
408 const double mu_range =
gWater(E);
410 double noise_likelihood = 0.0;
411 int number_of_hits = 0;
413 for (vector<JNPEHit>::const_iterator i = data.begin(); i != data.end(); ++i) {
414 noise_likelihood += log10(
getP(i->getY0(), i->getN()));
415 number_of_hits += i->getN();
420 fit.push_back(event());
430 out.rbegin()->setW(track->getW());
448 copy(input.
in.begin(), input.
in.end(), back_inserter(out));
Coverage of dynamical detector calibration.
Auxiliary class to extract a subset of optical modules from a detector.
Data structure for detector geometry and calibration.
Data regression method for JFIT::JEnergy.
Basic data structure for L0 hit.
#define DEBUG(A)
Message macros.
Direct access to module in detector data structure.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Router for direct addressing of module data in detector data structure.
bool hasModule(const JObjectID &id) const
Has module.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Data structure for a composite optical module.
Data structure for PMT geometry, calibration and status.
Data structure for fit of energy.
Data structure for set of track fit results.
void select(const JSelector_t &selector)
Select fits.
Data structure for track fit results with history and optional associated values.
void setE(const double E)
Set energy.
Data structure for fit of straight line paralel to z-axis.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
void transform(const JAxis3D &axis)
Transform axis to reference frame of given axis.
Data structure for position in three dimensions.
void transform(const JRotation3D &R, const JVector3D &pos)
Transform position.
JPosition3D & rotate(const JRotation3D &R)
Rotate.
const JPosition3D & getPosition() const
Get position.
double getY() const
Get y position.
double getZ() const
Get z position.
double getX() const
Get x position.
const JClass_t & getReference() const
Get reference to object.
Auxiliary class for correction of energy determined by JEnergy.cc.
Auxiliary class to to determine muon energy.
JMuonEnergy(const JMuonEnergyParameters_t ¶meters, const storage_type &storage, const JEnergyCorrection &correct, const int debug=0)
Constructor.
JEvt operator()(const input_type &input)
Fit function.
JRegressor< JEnergy > JRegressor_t
std::vector< module_type > detector_type
input_type getInput(const JModuleRouter &router, const JSummaryRouter &summary, const JDAQEvent &event, const JEvt &in, const coverage_type &coverage) const
Get input data.
JEnergyCorrection correct
Router for fast addressing of summary data in KM3NETDAQ::JDAQSummaryslice data structure as a functio...
const JDAQSummaryFrame & getSummaryFrame() const
Get default summary frame.
Reduced data structure for L0 hit.
JPMT_t getPMT() const
Get PMT.
double getT() const
Get calibrated time of hit.
2-dimensional frame with time calibrated data from one optical module.
Data storage class for rate measurements of all PMTs in one module.
static const int JMUONENERGY
static const int JENERGY_NDF
number of degrees of freedom from JEnergy.cc
static const int JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
static const int JENERGY_NOISE_LIKELIHOOD
log likelihood of every hit being K40 from JEnergy.cc
static const int JENERGY_CHI2
chi2 from JEnergy.cc
static const int JENERGY_MUON_RANGE_METRES
range of a muon with the reconstructed energy [m] from JEnergy.cc
static const int JPP_COVERAGE_POSITION
coverage of dynamic position calibration from any Jpp application
static const int JENERGY_NUMBER_OF_HITS
number of hits from JEnergy.cc
static const int JSTART_LENGTH_METRES
distance between projected positions on the track of optical modules for which the response does not ...
static const int JPP_COVERAGE_ORIENTATION
coverage of dynamic orientation calibration from any Jpp application
static const int JENERGY_MINIMAL_ENERGY
minimal energy [GeV] from JEnergy.cc
static const int JENERGY_MAXIMAL_ENERGY
maximal energy [GeV] from JEnergy.cc
void copy(const Head &from, JHead &to)
Copy header from from to to.
double getNPE(const Hit &hit)
Get true charge of hit.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
double getP(const double expval, bool hit)
Get Poisson probability to observe a hit or not for given expectation value for the number of hits.
double getKappaC()
Get average R-dependence of arrival time of Cherenkov light (a.k.a.
const double getInverseSpeedOfLight()
Get inverse speed of light.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
JPosition3D getPosition(const JFit &fit)
Get position.
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
JDirection3D getDirection(const JFit &fit)
Get direction.
Auxiliary classes and methods for triggering.
Auxiliary data structure for sequence of same character.
Auxiliary data structure for floating point format specification.
Data structure for coverage of detector by dynamical calibrations.
double position
coverage of detector by available position calibration [0,1]
double orientation
coverage of detector by available orientation calibration [0,1]
Auxiliary class for historical event.
Auxiliary class to test history.
double getY0() const
Get light yield due to random background.
Auxiliary class for simultaneously handling light yields and response of PMT.
size_t getN() const
Get number of hits.
double getYB() const
Get light yield due to bremsstrahlung.
double getY1() const
Get light yield due to minimum ionizing particle.
double getYA() const
Get light yield due to delta-rays.
Template definition of a data regressor of given model.
Auxiliary class for handling module response.
Data structure for fit parameters.
size_t numberOfPrefits
number of prefits
double roadWidth_m
road width [m]
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double resolution
energy resolution [log10(GeV)]
int mestimator
M-estimator.
double EMin_log
minimal energy [log10(GeV)]
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
double ZMin_m
minimal z-position [m]
double EMax_log
maximal energy [log10(GeV)]
Auxiliary class for energy estimation.
JResult(const JEnergy &x=0.0, const double chi2=std::numeric_limits< double >::max())
Constructor.
Auxiliary data structure for floating point format specification.