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;
159 const JBuildL0 <JHitR0> buildL0;
166 for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
168 if (router.
hasModule(i->getModuleID())) {
170 buffer(*i, router.
getModule(i->getModuleID()));
172 buildL0(buffer, back_inserter(data[i->getModuleID()]));
177 if (!module.empty()) {
197 const double chi2 = std::numeric_limits<double>::max())
208 operator bool()
const
210 return chi2 != std::numeric_limits<double>::max();
243 for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
248 double zmin = numeric_limits<double>::lowest();
256 for (
const auto& module : input.
data) {
269 for (
size_t i = 0; i !=
module->size(); ++i) {
271 if (module.getStatus(i)) {
277 return (hit.
getPMT() == pmt && T_ns(hit.
getT()));
283 } match = { T_ns + t1, i };
285 JPMT pmt =
module->getPMT(i);
289 const JNPEHit hit(this->
getNPE(pmt, module.frame.getRate(i)), count_if(module.begin(), module.end(), match));
291 DEBUG(
"hit: " << setw(8) << module->getID() <<
'.' <<
FILL(2,
'0') << i <<
' '
293 <<
FIXED(7,3) << module.frame.getRate(i) * 1.0e-3 <<
' '
298 << setw(2) << hit.
getN() << endl);
307 const int NDF =
distance(data.begin(), data.end()) - 1;
317 for (
int i = 0; i != N; ++i) {
327 for (
int i = 0; i != N; ++i) {
332 const double chi2 = (*this)(x, data.begin(), data.end());
335 buffer[chi2] = x.getE();
344 for (
int i = 0; i != N; ++i) {
390 result[2].chi2 = (*this)(
result[2].x, data.begin(), data.end());
394 const double chi2 =
result[2].chi2;
395 const double E =
result[2].x.getE();
399 double Emin = numeric_limits<double>::max();
400 double Emax = numeric_limits<double>::lowest();
403 if (i->second < Emin) { Emin = i->second; }
404 if (i->second > Emax) { Emax = i->second; }
407 const double mu_range =
gWater(E);
409 double noise_likelihood = 0.0;
410 int number_of_hits = 0;
412 for (vector<JNPEHit>::const_iterator i = data.begin(); i != data.end(); ++i) {
413 noise_likelihood += log10(
getP(i->getY0(), i->getN()));
414 number_of_hits += i->getN();
419 fit.push_back(event());
429 out.rbegin()->setW(track->getW());
447 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 JDAQModuleIdentifier &module) const
Get 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 resolution
energy resolution [log10(GeV)]
int mestimator
M-estimator.
double EMin_log
minimal energy [log10(GeV)]
double R_Hz
default rate [Hz]
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.