132 using namespace KM3NETDAQ;
135 typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
138 JParallelFileScanner_t inputFile;
146 size_t numberOfPrefits;
156 JParser<> zap(
"Program to perform fit of muon energy to data.");
161 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
163 zap[
'R'] =
make_field(roadWidth_m) = numeric_limits<double>::max();
168 zap[
'x'] =
make_field(X) = JEnergyRange(1.0, 8.0);
171 zap[
'M'] =
make_field(mestimator) = EM_NORMAL, EM_LORENTZIAN, EM_LINEAR, EM_NULL;
176 catch(
const exception& error) {
177 FATAL(error.what() << endl);
203 JRegressor_t::T_ns = T_ns;
205 const JEnergy ENERGY_RESOLUTION(0.01);
208 JRegressor_t fit(pdfFile);
210 switch (mestimator) {
213 NOTICE(
"Using the normal M-Estimator." << endl);
218 NOTICE(
"Using the Lorentzian M-Estimator." << endl);
223 NOTICE(
"Using the linear M-Estimator." << endl);
228 NOTICE(
"Using the JEnergy likelihood directly." << endl);
233 FATAL(
"Invalid M-Estimator." << endl);
237 if (fit.getRmax() < roadWidth_m) {
238 roadWidth_m = fit.getRmax();
246 while (inputFile.hasNext()) {
248 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
250 multi_pointer_type ps = inputFile.next();
256 JEvt::iterator __end = evt->end();
262 if (evt->begin() != __end) {
264 copy(evt->begin(), __end, back_inserter(out));
266 if (numberOfPrefits > 0) {
267 advance(__end = evt->begin(), min(numberOfPrefits, out.size()));
270 partial_sort(evt->begin(), __end, evt->end(),
qualitySorter);
272 __end = partition(evt->begin(), __end,
JHistory::is_event(evt->begin()->getHistory()));
277 buildL0(*tev, moduleRouter,
true, back_inserter(dataL0));
280 for (JEvt::const_iterator track = evt->begin(); track != __end; ++track) {
289 for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
296 top.insert(i->getPMTIdentifier());
306 for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
310 for (
size_t i = 0; i != module->size(); ++i) {
314 data.push_back(
JNPEHit(fit.getNPE(module->getPMT(i), rates_Hz.getSinglesRate()), count));
320 buffer.push_back(
JNPEHit(fit.getNPE(*module, rates_Hz), total));
328 if (start.is_valid() && !buffer.empty()) {
334 if (track_start != buffer.end()) {
338 __begin = lower_bound(data.begin(), data.end(), track_start->getZ(),
make_comparator(&JNPEHit::getZ));
343 const int NDF = distance(__begin, data.end()) - 1;
354 for (
int i = 0; i != N; ++i) {
355 result[i].x = X.getLowerLimit() + i * (X.getUpperLimit() - X.getLowerLimit()) / (N-1);
362 for (
int i = 0; i != N; ++i) {
365 result[i].chi2 = fit(result[i].x, __begin, data.end());
368 if (result[i].chi2 < result[j].chi2) {
378 result[4] = result[1];
379 result[2] = JResult(0.5 * (result[0].x + result[4].x));
383 result[4] = result[2];
384 result[2] = result[1];
388 result[0] = result[1];
389 result[4] = result[3];
393 result[0] = result[2];
394 result[2] = result[3];
398 result[0] = result[3];
399 result[2] = JResult(0.5 * (result[0].x + result[4].x));
403 result[1] = JResult(0.5 * (result[0].x + result[2].x));
404 result[3] = JResult(0.5 * (result[2].x + result[4].x));
406 }
while (result[4].x - result[0].x > ENERGY_RESOLUTION);
409 if (result[1].chi2 != result[3].chi2) {
411 result[2].x += 0.25 * (result[3].x - result[1].x) * (result[1].chi2 - result[3].chi2) / (result[1].chi2 + result[3].chi2 - 2*result[2].chi2);
412 result[2].chi2 = fit(result[2].x, __begin, data.end());
415 const double chi2 = result[2].chi2;
416 const double E = result[2].x.getE();
421 double noise_likelihood = 0.0;
426 noise_likelihood += log10(dom->getP());
427 nhits += dom->getN();
435 out.rbegin()->setE(correct(E));
Auxiliary class for start or end point evaluation.
Utility class to parse command line options.
number of degrees of freedom from JEnergy.cc
static const JGeane gWater(2.67e-1 *JTOOLS::DENSITY_SEA_WATER, 3.4e-4 *JTOOLS::DENSITY_SEA_WATER)
Function object for Energy loss of muon in sea water.
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
Regressor function object for JEnergy fit.
Template specialisation of L0 builder for JHitL0 data type.
Router for direct addressing of module data in detector data structure.
log likelihood of every hit being K40 from JEnergy.cc
General purpose class for parallel reading of objects from a single file or multiple files...
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
Auxiliary class to test history.
Data structure for track fit results.
Auxiliary class for defining the range of iterations of objects.
Auxiliary class for correction of energy determined by JEnergy.cc.
T & getInstance(const T &object)
Get static instance from temporary object.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Auxiliary class to test history.
JAxis3D getAxis(const Trk &track)
Get axis.
Auxiliary class for simultaneously handling light yields and response of PMT.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
number of hits from JEnergy.cc
Detector subset without binary search functionality.
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
uncorrected energy [GeV] from JEnergy.cc
Data structure for set of track fit results.
JFit & add(const JFitApplication_t &type)
Add event to history.
Data structure for fit of straight line paralel to z-axis.
bool qualitySorter(const JFIT::JFit &first, const JFIT::JFit &second)
Comparison of fit results.
void copy(const Head &from, JHead &to)
Copy header from from to to.
Reduced data structure for L1 hit.
Data structure for fit of energy.
JDirection3D getDirection(const Vec &v)
Get direction.
Object reading from a list of files.
const JLimit & getLimit() const
Get limit.
range of a muon with the reconstructed energy [m] from JEnergy.cc
JPosition3D & rotate(const JRotation3D &R)
Rotate.
#define DEBUG(A)
Message macros.
JPosition3D getPosition(const Vec &v)
Get position.
Auxiliary class for K40 rates.