Jpp
JEnergy.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 #include <vector>
6 #include <set>
7 #include <algorithm>
8 #include <limits>
9 
10 #include "TMath.h"
11 
12 #include "evt/Head.hh"
13 #include "evt/Evt.hh"
14 #include "JDAQ/JDAQEvent.hh"
15 #include "JDAQ/JDAQTimeslice.hh"
16 #include "JDAQ/JDAQSummaryslice.hh"
17 
18 #include "JDetector/JDetector.hh"
19 #include "JDetector/JDetectorSubset.hh"
20 #include "JDetector/JDetectorToolkit.hh"
21 #include "JDetector/JModuleRouter.hh"
22 #include "JDetector/JK40Rates.hh"
23 
24 #include "JTrigger/JHit.hh"
25 #include "JTrigger/JFrame.hh"
26 #include "JTrigger/JTimeslice.hh"
27 #include "JTrigger/JHitL0.hh"
28 #include "JTrigger/JHitR1.hh"
29 #include "JTrigger/JBuildL0.hh"
30 #include "JTrigger/JTriggerParameters.hh"
31 
32 #include "JSupport/JMultipleFileScanner.hh"
33 #include "JSupport/JParallelFileScanner.hh"
34 #include "JSupport/JFileRecorder.hh"
35 #include "JSupport/JSupport.hh"
36 #include "JSupport/JMeta.hh"
37 
38 #include "JFit/JLine1Z.hh"
39 #include "JFit/JEnergy.hh"
40 #include "JFit/JEnergyRegressor.hh"
41 #include "JFit/JFitToolkit.hh"
42 #include "JFit/JFitParameters.hh"
43 #include "JFit/JEvt.hh"
44 #include "JFit/JEvtToolkit.hh"
45 #include "JFit/JEnergyCorrection.hh"
46 #include "JFit/JNPEHit.hh"
47 #include "JFit/JStart.hh"
48 #include "JFit/JMEstimator.hh"
49 #include "JFit/JTimeRange.hh"
50 #include "JFit/JModel.hh"
51 
52 #include "JTools/JConstants.hh"
53 #include "JTools/JFunction1D_t.hh"
54 #include "JTools/JFunctionalMap_t.hh"
55 
56 #include "JPhysics/JPDFTable.hh"
57 #include "JPhysics/JPDFToolkit.hh"
58 #include "JPhysics/JNPETable.hh"
59 
60 #include "JLang/JComparator.hh"
61 
62 #include "Jeep/JParser.hh"
63 #include "Jeep/JMessage.hh"
64 
65 
66 namespace {
67 
68  using namespace JPP;
69 
70 
71  /**
72  * Auxiliary class for energy estimation.
73  */
74  struct JResult {
75  /**
76  * Constructor.
77  *
78  * \param x Energy [log(E/GeV)]
79  * \param chi2 chi2
80  */
81  JResult(const JEnergy& x = 0.0,
82  const double chi2 = std::numeric_limits<double>::max())
83  {
84  this->x = x;
85  this->chi2 = chi2;
86  }
87 
88 
89  /**
90  * Type conversion.
91  *
92  * \return true if valid chi2; else false
93  */
94  operator bool() const
95  {
96  return chi2 != std::numeric_limits<double>::max();
97  }
98 
99  JEnergy x; //!< Energy
100  double chi2; //!< chi2
101  };
102 
103 
104  /**
105  * Type definition of energy range.
106  */
107  typedef JTOOLS::JRange<JEnergy> JEnergyRange;
108 
109 
110  /**
111  * Definition of the various M-Estimators available to use
112  */
113  enum JMEstimator_t {
114  EM_NORMAL = 0,
115  EM_LORENTZIAN = 1,
116  EM_LINEAR = 2,
117  EM_NULL = 3
118  };
119 }
120 
121 
122 /**
123  * \file
124  *
125  * Program to perform JFIT::JRegressor<JEnergy, JGandalf> fit using I/O of JFIT::JEvt data.
126  * \author mdejong
127  */
128 int main(int argc, char **argv)
129 {
130  using namespace std;
131  using namespace JPP;
132  using namespace KM3NETDAQ;
133 
134  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> > JParallelFileScanner_t;
135  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
136  typedef JTYPELIST<JAllTypes_t, JEvt>::typelist typelist;
137 
138  JParallelFileScanner_t inputFile;
139  JFileRecorder<typelist> outputFile;
140  JLimit_t& numberOfEvents = inputFile.getLimit();
141  string detectorFile;
142  string pdfFile;
143  double roadWidth_m;
144  JTimeRange T_ns;
145  JK40Rates rates_Hz;
146  size_t numberOfPrefits;
147  JEnergyCorrection correct;
148  JEnergyRange X;
149  JStart start;
150  bool reprocess;
151  int mestimator;
152  int debug;
153 
154  try {
155 
156  JParser<> zap("Program to perform fit of muon energy to data.");
157 
158  zap['f'] = make_field(inputFile);
159  zap['o'] = make_field(outputFile) = "energy.root";
160  zap['a'] = make_field(detectorFile);
161  zap['n'] = make_field(numberOfEvents) = JLimit::max();
162  zap['P'] = make_field(pdfFile);
163  zap['R'] = make_field(roadWidth_m) = numeric_limits<double>::max();
164  zap['T'] = make_field(T_ns) = JTimeRange(-50.0, +450.0); // ns
165  zap['B'] = make_field(rates_Hz) = JK40Rates::getInstance();
166  zap['N'] = make_field(numberOfPrefits) = 0;
167  zap['E'] = make_field(correct) = JEnergyCorrection(); // off
168  zap['x'] = make_field(X) = JEnergyRange(1.0, 8.0); // log10(E/GeV)
169  zap['s'] = make_field(start) = JStart(); // off
170  zap['r'] = make_field(reprocess);
171  zap['M'] = make_field(mestimator) = EM_NORMAL, EM_LORENTZIAN, EM_LINEAR, EM_NULL;
172  zap['d'] = make_field(debug) = 1;
173 
174  zap(argc, argv);
175  }
176  catch(const exception& error) {
177  FATAL(error.what() << endl);
178  }
179 
180 
181 
182  cout.tie(&cerr);
183 
184 
185  JDetector detector;
186 
187  try {
188  load(detectorFile, detector);
189  }
190  catch(const JException& error) {
191  FATAL(error);
192  }
193 
194  const JModuleRouter moduleRouter(detector);
195 
196  typedef vector<JHitL0> JDataL0_t;
197  const JBuildL0<JHitL0> buildL0;
198 
199 
200  typedef JRegressor<JEnergy> JRegressor_t;
201 
202  JRegressor_t::debug = debug;
203  JRegressor_t::T_ns = T_ns;
204 
205  const JEnergy ENERGY_RESOLUTION(0.01); // Energy resolution log10(E/GeV)
206 
207 
208  JRegressor_t fit(pdfFile);
209 
210  switch (mestimator) {
211 
212  case EM_NORMAL:
213  NOTICE("Using the normal M-Estimator." << endl);
214  fit.estimator.reset(new JMEstimatorNormal());
215  break;
216 
217  case EM_LORENTZIAN:
218  NOTICE("Using the Lorentzian M-Estimator." << endl);
219  fit.estimator.reset(new JMEstimatorLorentzian());
220  break;
221 
222  case EM_LINEAR:
223  NOTICE("Using the linear M-Estimator." << endl);
224  fit.estimator.reset(new JMEstimatorNormal());
225  break;
226 
227  case EM_NULL:
228  NOTICE("Using the JEnergy likelihood directly." << endl);
229  fit.estimator.reset(new JMEstimatorNull());
230  break;
231 
232  default:
233  FATAL("Invalid M-Estimator." << endl);
234  }
235 
236 
237  if (fit.getRmax() < roadWidth_m) {
238  roadWidth_m = fit.getRmax();
239  }
240 
241 
242  outputFile.open();
243 
244  outputFile.put(JMeta(argc, argv));
245 
246  while (inputFile.hasNext()) {
247 
248  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
249 
250  multi_pointer_type ps = inputFile.next();
251 
252  JDAQEvent* tev = ps;
253  JEvt* evt = ps;
254  JEvt out;
255 
256  JEvt::iterator __end = evt->end();
257 
258  if (reprocess) {
259  __end = partition(evt->begin(), __end, JHistory::is_not_event(JMUONENERGY));
260  }
261 
262  if (evt->begin() != __end) {
263 
264  copy(evt->begin(), __end, back_inserter(out));
265 
266  if (numberOfPrefits > 0) {
267  advance(__end = evt->begin(), min(numberOfPrefits, out.size()));
268  }
269 
270  partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
271 
272  __end = partition(evt->begin(), __end, JHistory::is_event(evt->begin()->getHistory()));
273 
274 
275  JDataL0_t dataL0;
276 
277  buildL0(*tev, moduleRouter, true, back_inserter(dataL0));
278 
279 
280  for (JEvt::const_iterator track = evt->begin(); track != __end; ++track) {
281 
282  const JRotation3D R (getDirection(*track));
283  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
284  const JModel<JLine1Z> match(tz, roadWidth_m, JRegressor_t::T_ns);
285 
286 
287  multiset<JDAQPMTIdentifier> top;
288 
289  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
290 
291  JHitR1 hit(*i);
292 
293  hit.rotate(R);
294 
295  if (match(hit)) {
296  top.insert(i->getPMTIdentifier());
297  }
298  }
299 
300 
301  vector<JNPEHit> data;
302  vector<JNPEHit> buffer;
303 
304  const JDetectorSubset_t subdetector(detector, getAxis(*track), roadWidth_m);
305 
306  for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
307 
308  size_t total = 0;
309 
310  for (size_t i = 0; i != module->size(); ++i) {
311 
312  size_t count = top.count(JDAQPMTIdentifier(module->getID(), i));
313 
314  data.push_back(JNPEHit(fit.getNPE(module->getPMT(i), rates_Hz.getSinglesRate()), count));
315 
316  total += count;
317  }
318 
319  if (total != 0) {
320  buffer.push_back(JNPEHit(fit.getNPE(*module, rates_Hz), total));
321  }
322  }
323 
324  // start position of track
325 
326  vector<JNPEHit>::iterator __begin = data.begin();
327 
328  if (start.is_valid() && !buffer.empty()) {
329 
330  sort(buffer.begin(), buffer.end(), make_comparator(&JNPEHit::getZ));
331 
332  vector<JNPEHit>::const_iterator track_start = start.find(buffer.begin(), buffer.end());
333 
334  if (track_start != buffer.end()) {
335 
336  sort(data.begin(), data.end(), make_comparator(&JNPEHit::getZ));
337 
338  __begin = lower_bound(data.begin(), data.end(), track_start->getZ(), make_comparator(&JNPEHit::getZ));
339  }
340  }
341 
342 
343  const int NDF = distance(__begin, data.end()) - 1;
344 
345  if (NDF >= 0) {
346 
347 
348  // 5-point search between given limits
349 
350  const int N = 5;
351 
352  JResult result[N];
353 
354  for (int i = 0; i != N; ++i) {
355  result[i].x = X.getLowerLimit() + i * (X.getUpperLimit() - X.getLowerLimit()) / (N-1);
356  }
357 
358  do {
359 
360  int j = 0;
361 
362  for (int i = 0; i != N; ++i) {
363 
364  if (!result[i]) {
365  result[i].chi2 = fit(result[i].x, __begin, data.end());
366  }
367 
368  if (result[i].chi2 < result[j].chi2) {
369  j = i;
370  }
371  }
372 
373  // squeeze range
374 
375  switch (j) {
376 
377  case 0:
378  result[4] = result[1];
379  result[2] = JResult(0.5 * (result[0].x + result[4].x));
380  break;
381 
382  case 1:
383  result[4] = result[2];
384  result[2] = result[1];
385  break;
386 
387  case 2:
388  result[0] = result[1];
389  result[4] = result[3];
390  break;
391 
392  case 3:
393  result[0] = result[2];
394  result[2] = result[3];
395  break;
396 
397  case 4:
398  result[0] = result[3];
399  result[2] = JResult(0.5 * (result[0].x + result[4].x));
400  break;
401  }
402 
403  result[1] = JResult(0.5 * (result[0].x + result[2].x));
404  result[3] = JResult(0.5 * (result[2].x + result[4].x));
405 
406  } while (result[4].x - result[0].x > ENERGY_RESOLUTION);
407 
408 
409  if (result[1].chi2 != result[3].chi2) {
410 
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());
413  }
414 
415  const double chi2 = result[2].chi2;
416  const double E = result[2].x.getE();
417 
418  // Calculate additional variables
419  const double muRange = JPHYSICS::gWater(E); // The range of a muon with the reconstructed energy
420 
421  double noise_likelihood = 0.0; // The log likelihood of every hit being from K40
422  int nhits = 0; // The number of hits selected for JEnergy
423 
424  for (vector<JNPEHit>::const_iterator dom = data.begin(); dom != data.end(); ++dom) {
425 
426  noise_likelihood += log10(dom->getP()); // The probability of getting the observed multiplicity with K40
427  nhits += dom->getN(); // DOM multiplicity
428  }
429 
430 
431  out.push_back(JFit(*track).add(JMUONENERGY));
432 
433  // set corrected energy
434 
435  out.rbegin()->setE(correct(E));
436 
437  // set additional values
438 
439  out.rbegin()->setW(track->getW());
440  out.rbegin()->setW(JENERGY_ENERGY, E);
441  out.rbegin()->setW(JENERGY_CHI2, chi2);
442  out.rbegin()->setW(JENERGY_MUON_RANGE_METRES, muRange);
443  out.rbegin()->setW(JENERGY_NOISE_LIKELIHOOD, noise_likelihood);
444  out.rbegin()->setW(JENERGY_NDF, NDF);
445  out.rbegin()->setW(JENERGY_NUMBER_OF_HITS, nhits);
446  }
447  }
448 
449  // apply default sorter
450 
451  sort(out.begin(), out.end(), qualitySorter);
452  }
453 
454  outputFile.put(out);
455  }
456  STATUS(endl);
457 
458  JSingleFileScanner<JRemove<typelist, JEvt>::typelist> io(inputFile);
459 
460  io >> outputFile;
461 
462  outputFile.close();
463 }
JFIT::JMEstimatorLorentzian
Lorentzian M-estimator.
Definition: JMEstimator.hh:79
JEnergyRegressor.hh
JMeta.hh
JFIT::JMUONENERGY
JEnergy.cc.
Definition: JFitApplications.hh:26
KM3NETDAQ::JDAQEvent
DAQ Event.
Definition: JDAQEvent.hh:34
JTriggerParameters.hh
JFIT::JENERGY_NDF
number of degrees of freedom from JEnergy.cc
Definition: JFitParameters.hh:32
JFIT::JEnergy
Data structure for fit of energy.
Definition: JEnergy.hh:24
JFileRecorder.hh
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
JFIT::JENERGY_CHI2
chi2 from JEnergy.cc
Definition: JFitParameters.hh:22
JMessage.hh
JROOT::advance
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
Definition: JCounter.hh:35
JPHYSICS::gWater
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:328
JAANET::getAxis
JAxis3D getAxis(const Trk &track)
Get axis.
Definition: JAAnetToolkit.hh:244
JLANG::getInstance
T & getInstance(const T &object)
Get static instance from temporary object.
Definition: JObject.hh:73
JFitToolkit.hh
qualitySorter
bool qualitySorter(const JFIT::JFit &first, const JFIT::JFit &second)
Comparison of fit results.
Definition: quality_sorter.cc:10
JFIT::JStart
Auxiliary class for start or end point evaluation.
Definition: JStart.hh:21
JTimeslice.hh
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:456
JComparator.hh
JPDFToolkit.hh
std::vector< JHitL0 >
JTOOLS::JRange
Range of values.
Definition: JRange.hh:34
JFIT::JNPEHit
Auxiliary class for simultaneously handling light yields and response of PMT.
Definition: JNPEHit.hh:19
JEvtToolkit.hh
JFIT::JENERGY_NOISE_LIKELIHOOD
log likelihood of every hit being K40 from JEnergy.cc
Definition: JFitParameters.hh:31
JTOOLS::j
int j
Definition: JPolint.hh:634
JK40Rates.hh
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition: PhysicsEvent.hh:434
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JDAQTimeslice.hh
JFIT::JEvt
Data structure for set of track fit results.
Definition: JEvt.hh:292
JFIT::JStart::is_valid
bool is_valid() const
Check validity of start or end point conditions.
Definition: JStart.hh:50
JEnergyCorrection.hh
JAANET::copy
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:152
JNPEHit.hh
JDAQSummaryslice.hh
JGEOMETRY3D::JPosition3D::rotate
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:185
JTOOLS::JTimeRange
JRange< double > JTimeRange
Type definition for time range.
Definition: JTools/JTimeRange.hh:19
JSupport.hh
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JDETECTOR::JK40Rates
Auxiliary class for K40 rates.
Definition: JK40Rates.hh:41
JFIT::JRegressor< JEnergy >
Regressor function object for JEnergy fit.
Definition: JEnergyRegressor.hh:49
JEvt.hh
JFunction1D_t.hh
JFIT::JMEstimatorNull
Null M-estimator.
Definition: JMEstimator.hh:53
JStart.hh
JBuildL0.hh
JDETECTOR::JK40Rates::getSinglesRate
double getSinglesRate() const
Get singles rate.
Definition: JK40Rates.hh:97
debug
int debug
debug level
Definition: JSirene.cc:59
JConstants.hh
JTOOLS::result
return result
Definition: JPolint.hh:695
JFunctionalMap_t.hh
JFIT::JModel< JLine1Z >
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
Definition: JModel.hh:34
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
std::multiset
Definition: JSTDTypes.hh:14
JHit.hh
JMEstimator_t
JMEstimator_t
Definition of the various M-Estimators available to use.
Definition: JEnergy.cc:113
JModuleRouter.hh
JFIT::JHistory::is_event
Auxiliary class to test history.
Definition: JHistory.hh:101
JNPETable.hh
JParallelFileScanner.hh
main
int main(int argc, char **argv)
Definition: JEnergy.cc:128
JHitR1.hh
JFIT::JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
Definition: JFitParameters.hh:21
JMultipleFileScanner.hh
JTimeRange.hh
JLANG::JTypeList
Type list.
Definition: JTypeList.hh:22
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JFIT::JMEstimatorNormal
Normal M-estimator.
Definition: JMEstimator.hh:66
JParser.hh
JDetectorToolkit.hh
JFIT::JLine1Z
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
JFIT::JHistory::is_not_event
Auxiliary class to test history.
Definition: JHistory.hh:149
JDETECTOR::JModuleRouter
Router for direct addressing of module data in detector data structure.
Definition: JModuleRouter.hh:34
JModel.hh
JTRIGGER::JHitR1
Reduced data structure for L1 hit.
Definition: JHitR1.hh:31
JAANET::getDirection
JDirection3D getDirection(const Vec &v)
Get direction.
Definition: JAAnetToolkit.hh:221
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
JFitParameters.hh
JFIT::JENERGY_NUMBER_OF_HITS
number of hits from JEnergy.cc
Definition: JFitParameters.hh:33
JSUPPORT::JParallelFileScanner
General purpose class for parallel reading of objects from a single file or multiple files.
Definition: JParallelFileScanner.hh:31
KM3NETDAQ::JDAQPMTIdentifier
PMT identifier.
Definition: JDAQPMTIdentifier.hh:25
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:80
JAANET::getPosition
JPosition3D getPosition(const Vec &v)
Get position.
Definition: JAAnetToolkit.hh:197
JAANET::detector
Detector file.
Definition: JHead.hh:130
JMEstimator.hh
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std
Definition: jaanetDictionary.h:36
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
JEnergy.hh
JDAQEvent.hh
JFIT::JENERGY_MUON_RANGE_METRES
range of a muon with the reconstructed energy [m] from JEnergy.cc
Definition: JFitParameters.hh:30
JTRIGGER::JBuildL0< JHitL0 >
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:101
JDetector.hh
JLANG::make_comparator
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
Definition: JComparator.hh:143
JDetectorSubset.hh
JFIT::JFit::add
JFit & add(const JFitApplication_t &type)
Add event to history.
Definition: JEvt.hh:135
JHitL0.hh
JFIT::JStart::find
T find(T __begin, T __end) const
Get start or end point of muon trajectory.
Definition: JStart.hh:108
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JSUPPORT::JFileRecorder
Object writing to file.
Definition: JFileRecorder.hh:41
JLANG::JException
General exception.
Definition: JException.hh:40
JDETECTOR::JDetectorSubset_t
Detector subset without binary search functionality.
Definition: JDetectorSubset.hh:37
JPDFTable.hh
JGEOMETRY3D::JRotation3D
Rotation matrix.
Definition: JRotation3D.hh:111
JSUPPORT::JSingleFileScanner
Object reading from a list of files.
Definition: JSingleFileScanner.hh:75
JFrame.hh
JFIT::JEnergyCorrection
Auxiliary class for correction of energy determined by JEnergy.cc.
Definition: JEnergyCorrection.hh:39
JLine1Z.hh
JFIT::JFit
Data structure for track fit results.
Definition: JEvt.hh:31
Generated by   doxygen 1.8.16