Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Enumerations | Functions
JEnergy.cc File Reference

Program to perform JFIT::JRegressor<JEnergy, JGandalf> fit using I/O of JFIT::JEvt data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <algorithm>
#include <limits>
#include "TMath.h"
#include "evt/Head.hh"
#include "evt/Evt.hh"
#include "JDAQ/JDAQEvent.hh"
#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/JDAQSummaryslice.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorSubset.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JK40Rates.hh"
#include "JTrigger/JHit.hh"
#include "JTrigger/JFrame.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitR1.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JFit/JLine1Z.hh"
#include "JFit/JEnergy.hh"
#include "JFit/JEnergyRegressor.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JFitParameters.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JEnergyCorrection.hh"
#include "JFit/JNPEHit.hh"
#include "JFit/JStart.hh"
#include "JFit/JMEstimator.hh"
#include "JFit/JTimeRange.hh"
#include "JFit/JModel.hh"
#include "JTools/JConstants.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JNPETable.hh"
#include "JLang/JComparator.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Enumerations

enum  JMEstimator_t
 Definition of the various M-Estimators available to use. More...
 

Functions

int main (int argc, char **argv)
 

Detailed Description

Program to perform JFIT::JRegressor<JEnergy, JGandalf> fit using I/O of JFIT::JEvt data.

Author
mdejong

Definition in file JEnergy.cc.

Enumeration Type Documentation

Definition of the various M-Estimators available to use.

Definition at line 113 of file JEnergy.cc.

113  {
114  EM_NORMAL = 0,
115  EM_LORENTZIAN = 1,
116  EM_LINEAR = 2,
117  EM_NULL = 3
118  };

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 128 of file JEnergy.cc.

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;
137 
138  JParallelFileScanner_t inputFile;
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) = 1;
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 
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 
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 
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(JENERGY_ENERGY, E);
440  out.rbegin()->setW(JENERGY_CHI2, chi2);
441  out.rbegin()->setW(JENERGY_MUON_RANGE_METRES, muRange);
442  out.rbegin()->setW(JENERGY_NOISE_LIKELIHOOD, noise_likelihood);
443  out.rbegin()->setW(JENERGY_NDF, NDF);
444  out.rbegin()->setW(JENERGY_NUMBER_OF_HITS, nhits);
445  }
446  }
447 
448  // apply default sorter
449 
450  sort(out.begin(), out.end(), qualitySorter);
451  }
452 
453  outputFile.put(out);
454  }
455  STATUS(endl);
456 
458 
459  io >> outputFile;
460 
461  outputFile.close();
462 }
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
Object writing to file.
Auxiliary class for start or end point evaluation.
Definition: JStart.hh:21
Utility class to parse command line options.
Definition: JParser.hh:1410
General exception.
Definition: JException.hh:40
number of degrees of freedom from JEnergy.cc
chi2 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.
Null M-estimator.
Definition: JMEstimator.hh:53
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
Regressor function object for JEnergy fit.
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:102
#define STATUS(A)
Definition: JMessage.hh:61
Detector data structure.
Definition: JDetector.hh:77
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:108
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.
Definition: JModel.hh:31
Auxiliary class to test history.
Definition: JHistory.hh:149
string outputFile
JRange< double > JTimeRange
Type definition for time range.
Data structure for track fit results.
Definition: JEvt.hh:32
Type list.
Definition: JTypeList.hh:22
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Detector file.
Definition: JHead.hh:126
Auxiliary class for correction of energy determined by JEnergy.cc.
T & getInstance(const T &object)
Get static instance from temporary object.
Definition: JObject.hh:73
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
Auxiliary class to test history.
Definition: JHistory.hh:101
JAxis3D getAxis(const Trk &track)
Get axis.
#define NOTICE(A)
Definition: JMessage.hh:62
Auxiliary class for simultaneously handling light yields and response of PMT.
Definition: JNPEHit.hh:19
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Normal M-estimator.
Definition: JMEstimator.hh:66
int debug
debug level
Definition: JSirene.cc:59
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.
Definition: JCounter.hh:35
uncorrected energy [GeV] from JEnergy.cc
#define FATAL(A)
Definition: JMessage.hh:65
Data structure for set of track fit results.
Definition: JEvt.hh:312
JFit & add(const JFitApplication_t &type)
Add event to history.
Definition: JEvt.hh:136
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
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.
Definition: JHead.cc:40
Reduced data structure for L1 hit.
Definition: JHitR1.hh:31
Data structure for fit of energy.
Definition: JEnergy.hh:24
JDirection3D getDirection(const Vec &v)
Get direction.
Object reading from a list of files.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:72
range of a muon with the reconstructed energy [m] from JEnergy.cc
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:185
Lorentzian M-estimator.
Definition: JMEstimator.hh:79
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
JPosition3D getPosition(const Vec &v)
Get position.
Auxiliary class for K40 rates.
Definition: JK40Rates.hh:41