Jpp - the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMuonEnergy.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <set>
6 #include <algorithm>
7 #include <limits>
8 #include <map>
9 #include <memory>
10 
11 #include "TMath.h"
12 
16 #include "JDAQ/JDAQEventIO.hh"
17 #include "JDAQ/JDAQTimesliceIO.hh"
19 
20 #include "JDetector/JDetector.hh"
24 
25 #include "JDynamics/JDynamics.hh"
26 
27 #include "JTrigger/JHit.hh"
28 #include "JTrigger/JFrame.hh"
29 #include "JTrigger/JTimeslice.hh"
30 #include "JTrigger/JHitL0.hh"
31 #include "JTrigger/JHitR1.hh"
32 #include "JTrigger/JBuildL0.hh"
34 
38 #include "JSupport/JSupport.hh"
39 #include "JSupport/JMeta.hh"
41 
42 #include "JFit/JLine1Z.hh"
43 #include "JFit/JEnergy.hh"
44 #include "JFit/JEnergyRegressor.hh"
45 #include "JFit/JFitToolkit.hh"
46 #include "JFit/JNPEHit.hh"
47 #include "JFit/JMEstimator.hh"
48 #include "JFit/JTimeRange.hh"
49 #include "JFit/JModel.hh"
50 
51 #include "JReconstruction/JEvt.hh"
55 
56 #include "JPhysics/JConstants.hh"
57 #include "JTools/JRange.hh"
58 #include "JTools/JFunction1D_t.hh"
60 
61 #include "JPhysics/JPDFTable.hh"
62 #include "JPhysics/JPDFToolkit.hh"
63 #include "JPhysics/JNPETable.hh"
64 
65 #include "JLang/JComparator.hh"
66 
67 #include "Jeep/JParser.hh"
68 #include "Jeep/JMessage.hh"
69 
70 
71 namespace {
72 
73  using namespace JPP;
74 
75 
76  /**
77  * Auxiliary class for energy estimation.
78  */
79  struct JResult {
80  /**
81  * Constructor.
82  *
83  * \param x Energy [log(E/GeV)]
84  * \param chi2 chi2
85  */
86  JResult(const JEnergy& x = 0.0,
87  const double chi2 = std::numeric_limits<double>::max())
88  {
89  this->x = x;
90  this->chi2 = chi2;
91  }
92 
93 
94  /**
95  * Type conversion.
96  *
97  * \return true if valid chi2; else false
98  */
99  operator bool() const
100  {
101  return chi2 != std::numeric_limits<double>::max();
102  }
103 
104  JEnergy x; //!< Energy
105  double chi2; //!< chi2
106  };
107 
108 
109  /**
110  * Type definition of energy range.
111  */
112  typedef JTOOLS::JRange<JEnergy> JEnergyRange;
113 }
114 
115 
116 /**
117  * \file
118  *
119  * Program to perform JFIT::JRegressor<JEnergy, JGandalf> fit using I/O of JFIT::JEvt data.
120  * \author mdejong
121  */
122 int main(int argc, char **argv)
123 {
124  using namespace std;
125  using namespace JPP;
126  using namespace KM3NETDAQ;
127 
129  typedef JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> > JParallelFileScanner_t;
130  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
132  JACOUSTICS::JEvt>::typelist calibration_types;
133  typedef JMultipleFileScanner<calibration_types> JCalibration_t;
134 
135  JParallelFileScanner_t inputFile;
137  JLimit_t& numberOfEvents = inputFile.getLimit();
138  string detectorFile;
139  JCalibration_t calibrationFile;
140  double Tmax_s;
141  string pdfFile;
142  JEnergyCorrection correct;
144  int debug;
145 
146  try {
147 
148  JParser<> zap("Program to perform fit of muon energy to data.");
149 
150  zap['f'] = make_field(inputFile);
151  zap['o'] = make_field(outputFile) = "energy.root";
152  zap['a'] = make_field(detectorFile);
153  zap['+'] = make_field(calibrationFile) = JPARSER::initialised();
154  zap['T'] = make_field(Tmax_s) = 100.0;
155  zap['n'] = make_field(numberOfEvents) = JLimit::max();
156  zap['P'] = make_field(pdfFile);
157  zap['E'] = make_field(correct) = JEnergyCorrection(); // off
159  zap['d'] = make_field(debug) = 1;
160 
161  zap(argc, argv);
162  }
163  catch(const exception& error) {
164  FATAL(error.what() << endl);
165  }
166 
167 
168  cout.tie(&cerr);
169 
170 
172 
173  try {
174  load(detectorFile, detector);
175  }
176  catch(const JException& error) {
177  FATAL(error);
178  }
179 
180  getMechanics.load(detector.getID());
181 
182  unique_ptr<JDynamics> dynamics;
183 
184  try {
185 
186  dynamics.reset(new JDynamics(detector, Tmax_s));
187 
188  dynamics->load(calibrationFile);
189  }
190  catch(const exception& error) {
191  if (!calibrationFile.empty()) {
192  FATAL(error.what());
193  }
194  }
195 
196  const JModuleRouter router(dynamics ? dynamics->getDetector() : detector);
197 
198  JSummaryFileRouter summary(inputFile, parameters.R_Hz);
199 
200 
201  typedef vector<JHitL0> JDataL0_t;
202  const JBuildL0<JHitL0> buildL0;
203 
204  typedef JRegressor<JEnergy> JRegressor_t;
205 
207  JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
208 
209  JRegressor_t fit(pdfFile);
210 
211  fit.estimator.reset(getMEstimator(parameters.mestimator));
212 
213  if (!fit.estimator.is_valid()) {
214  FATAL("Invalid M-Estimator." << endl);
215  }
216 
217  if (fit.getRmax() < parameters.roadWidth_m) {
218  parameters.roadWidth_m = fit.getRmax();
219  }
220 
221 
222  outputFile.open();
223 
224  outputFile.put(JMeta(argc, argv));
225 
226  while (inputFile.hasNext()) {
227 
228  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
229 
230  multi_pointer_type ps = inputFile.next();
231 
232  const JDAQEvent* tev = ps;
233  const JFIT::JEvt* in = ps;
234 
235  summary.update(*tev);
236 
237  if (dynamics) {
238  dynamics->update(*tev);
239  }
240 
241  // select start values
242 
243  JFIT::JEvt cp = *in;
244  JFIT::JEvt out;
245 
246  if (parameters.reprocess) {
248  }
249 
250  cp.select(parameters.numberOfPrefits, qualitySorter);
251 
252  if (!cp.empty()) {
253  cp.select(JHistory::is_event(cp.begin()->getHistory()));
254  }
255 
256 
257  JDataL0_t dataL0;
258 
259  buildL0(*tev, router, true, back_inserter(dataL0));
260 
261 
262  for (JFIT::JEvt::const_iterator track = cp.begin(); track != cp.end(); ++track) {
263 
264  const JRotation3D R (getDirection(*track));
265  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
266  const JFIT::JModel<JLine1Z> match(tz, parameters.roadWidth_m, JRegressor_t::T_ns);
267 
268 
270 
271  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
272 
273  JHitR1 hit(*i);
274 
275  hit.rotate(R);
276 
277  if (match(hit)) {
278  top.insert(i->getPMTIdentifier());
279  }
280  }
281 
282 
283  double EMin_log = parameters.EMin_log;
284  double EMax_log = parameters.EMax_log;
285  JRange<double> Z_m;
286 
287  if (track->hasW(JSTART_LENGTH_METRES) && track->getW(JSTART_LENGTH_METRES) > 0.0) {
288 
289  Z_m.setLowerLimit(parameters.ZMin_m);
290 
291  //EMin_log = log10(track->getW(JSTART_LENGTH_METRES) * gWater.getA());
292  }
293 
294 
295  vector<JNPEHit> data;
296 
297  const JDetectorSubset_t subdetector(detector, getAxis(*track), parameters.roadWidth_m, Z_m);
298 
299  for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
300 
301  if (summary.hasSummaryFrame(module->getID())) {
302 
303  const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
304 
305  for (size_t i = 0; i != module->size(); ++i) {
306 
307  if (getDAQStatus(frame, *module, i) &&
308  getPMTStatus(frame, *module, i) &&
309  !module->getPMT(i).has(PMT_DISABLE)) {
310 
311  const double rate_Hz = frame.getRate(i);
312  const size_t count = top.count(JDAQPMTIdentifier(module->getID(), i));
313 
314  data.push_back(JNPEHit(fit.getNPE(module->getPMT(i), rate_Hz), count));
315  }
316  }
317  }
318  }
319 
320  const int NDF = distance(data.begin(), data.end()) - 1;
321 
322  if (NDF >= 0) {
323 
324 
325  // 5-point search between given limits
326 
327  const int N = 5;
328 
329  JResult result[N];
330 
331  for (int i = 0; i != N; ++i) {
332  result[i].x = EMin_log + i * (EMax_log - EMin_log) / (N-1);
333  }
334 
335  map<double, double> buffer;
336 
337  do {
338 
339  int j = 0;
340 
341  for (int i = 0; i != N; ++i) {
342 
343  if (!result[i]) {
344 
345  const JEnergy x = result[i].x;
346  const double chi2 = fit(x, data.begin(), data.end());
347 
348  result[i].chi2 = chi2;
349  buffer[chi2] = x.getE();
350  }
351 
352  if (result[i].chi2 < result[j].chi2) {
353  j = i;
354  }
355  }
356 
357  for (int i = 0; i != N; ++i) {
358  DEBUG(' ' << FIXED(5,2) << result[i].x << ' ' << FIXED(9,3) << result[i].chi2);
359  }
360  DEBUG(endl);
361 
362 
363  // squeeze range
364 
365  switch (j) {
366 
367  case 0:
368  result[4] = result[1];
369  result[2] = result[0];
370  result[0] = JResult(2*result[2].x - result[4].x);
371  break;
372 
373  case 1:
374  result[4] = result[2];
375  result[2] = result[1];
376  break;
377 
378  case 2:
379  result[0] = result[1];
380  result[4] = result[3];
381  break;
382 
383  case 3:
384  result[0] = result[2];
385  result[2] = result[3];
386  break;
387 
388  case 4:
389  result[0] = result[3];
390  result[2] = result[4];
391  result[4] = JResult(2*result[2].x - result[0].x);
392  break;
393  }
394 
395  result[1] = JResult(0.5 * (result[0].x + result[2].x));
396  result[3] = JResult(0.5 * (result[2].x + result[4].x));
397 
398  } while (result[4].x - result[0].x > parameters.resolution);
399 
400 
401  if (result[1].chi2 != result[3].chi2) {
402 
403  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);
404  result[2].chi2 = fit(result[2].x, data.begin(), data.end());
405  }
406 
407  const double chi2 = result[2].chi2;
408  const double E = result[2].x.getE();
409 
410  // calculate additional variables
411 
412  double Emin = numeric_limits<double>::max();
413  double Emax = numeric_limits<double>::lowest();
414 
415  for (map<double, double>::const_iterator i = buffer.begin(); i != buffer.end() && i->first <= chi2 + 1.0; ++i) {
416  if (i->second < Emin) { Emin = i->second; }
417  if (i->second > Emax) { Emax = i->second; }
418  }
419 
420  const double mu_range = gWater(E); // range of a muon with the reconstructed energy
421 
422  double noise_likelihood = 0.0; // log-likelihood of every hit being from K40
423  int number_of_hits = 0; // number of hits selected for JEnergy
424 
425  for (vector<JNPEHit>::const_iterator i = data.begin(); i != data.end(); ++i) {
426  noise_likelihood += log10(i->getP()); // probability of getting the observed multiplicity with K40
427  number_of_hits += i->getN(); // hit multiplicity
428  }
429 
430 
431  out.push_back(JFIT::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, mu_range);
443  out.rbegin()->setW(JENERGY_NOISE_LIKELIHOOD, noise_likelihood);
444  out.rbegin()->setW(JENERGY_NDF, NDF);
445  out.rbegin()->setW(JENERGY_NUMBER_OF_HITS, number_of_hits);
446  out.rbegin()->setW(JENERGY_MINIMAL_ENERGY, Emin);
447  out.rbegin()->setW(JENERGY_MAXIMAL_ENERGY, Emax);
448  }
449  }
450 
451  if (dynamics) {
452 
453  const JDynamics::coverage_type coverage = dynamics->getCoverage();
454 
455  for (JFIT::JEvt::iterator i = out.begin(); i != out.end(); ++i) {
456  i->setW(JPP_COVERAGE_ORIENTATION, coverage.orientation);
457  i->setW(JPP_COVERAGE_POSITION, coverage.position);
458  }
459  }
460 
461  // apply default sorter
462 
463  sort(out.begin(), out.end(), qualitySorter);
464 
465  copy(in->begin(), in->end(), back_inserter(out));
466 
467  outputFile.put(out);
468  }
469  STATUS(endl);
470 
472 
473  io >> outputFile;
474 
475  outputFile.close();
476 }
Auxiliary methods to evaluate Poisson probabilities and chi2.
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
double getE() const
Get energy.
Definition: JEnergy.hh:170
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
double getRate(const int tdc, const double factor=1.0) const
Get count rate.
Auxiliary methods for PDF calculations.
double position
coverage of detector by available position calibration [0,1]
Definition: JDynamics.hh:536
static JDetectorMechanics getMechanics
Function object to get string mechanics.
Definition: JMechanics.hh:243
static const int JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Regressor function object for JEnergy fit.
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:102
static const int JENERGY_MAXIMAL_ENERGY
maximal energy [GeV] from JEnergy.cc
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:80
Recording of objects on file according a format that follows from the file name extension.
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
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: JFit/JModel.hh:34
Auxiliary class to test history.
Definition: JHistory.hh:152
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:323
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
Basic data structure for time and time over threshold information of hit.
JFit & add(const int type)
Add event to history.
string outputFile
static const int JENERGY_CHI2
chi2 from JEnergy.cc
Data structure for detector geometry and calibration.
Various implementations of functional maps.
static const int JENERGY_NOISE_LIKELIHOOD
log likelihood of every hit being K40 from JEnergy.cc
Basic data structure for L0 hit.
Data structure for track fit results.
double orientation
coverage of detector by available orientation calibration [0,1]
Definition: JDynamics.hh:535
static const int JENERGY_NDF
number of degrees of freedom from JEnergy.cc
static const int JENERGY_MINIMAL_ENERGY
minimal energy [GeV] from JEnergy.cc
Type list.
Definition: JTypeList.hh:22
Scanning of objects from a single file according a format that follows from the extension of each fil...
Auxiliary class to extract a subset of optical modules from a detector.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Detector file.
Definition: JHead.hh:196
JDirection3D getDirection(const Vec &dir)
Get direction.
Auxiliary class for correction of energy determined by JEnergy.cc.
Acoustic event fit.
Auxiliary class for recursive type list generation.
Definition: JTypeList.hh:351
static const int JENERGY_NUMBER_OF_HITS
number of hits from JEnergy.cc
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
Data storage class for rate measurements of all PMTs in one module.
Auxiliary class to test history.
Definition: JHistory.hh:104
JAxis3D getAxis(const Trk &track)
Get axis.
return result
Definition: JPolint.hh:727
Parallel scanning of objects from a single file or multiple files according a format that follows fro...
static const int PMT_DISABLE
KM3NeT Data Definitions v2.0.0-15-g59d2e2b https://git.km3net.de/common/km3net-dataformat.
Definition: pmt_status.hh:12
ROOT I/O of application specific meta data.
JPosition3D getPosition(const Vec &pos)
Get position.
Physics constants.
Dynamic detector calibration.
File router for fast addressing of summary data.
Auxiliary class for simultaneously handling light yields and response of PMT.
Definition: JNPEHit.hh:19
void load(const std::string &file_name)
Load mechanical model parameters from file.
Definition: JMechanics.hh:142
int debug
debug level
Definition: JSirene.cc:63
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:40
Range of values.
Definition: JRange.hh:38
General purpose messaging.
Data structure for coverage of dynamic calibrations.
Definition: JDynamics.hh:534
Detector subset without binary search functionality.
Data structure for fit parameters.
#define FATAL(A)
Definition: JMessage.hh:67
Direct access to module in detector data structure.
static const int JPP_COVERAGE_ORIENTATION
coverage of dynamic orientation calibration from any Jpp application
Reduced data structure for L1 hit.
Dynamic detector calibration.
Definition: JDynamics.hh:61
bool getPMTStatus(const JStatus &status)
Test status of PMT.
static const int JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Data structure for set of track fit results.
std::vector< int > count
Definition: JAlgorithm.hh:184
Auxiliary class to define a range between two values.
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:139
then cp
Reduced data structure for L1 hit.
Definition: JHitR1.hh:31
Data structure for fit of energy.
Definition: JEnergy.hh:28
Orientation of module.
Object reading from a list of files.
int j
Definition: JPolint.hh:666
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
Data regression method for JFIT::JEnergy.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:166
do set_variable DETECTOR_TXT $WORKDIR detector
bool getDAQStatus(const JDAQFrameStatus &frame, const JStatus &status)
Test status of DAQ.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:38
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
void select(const JSelector_t &selector)
Select fits.
then usage $script[input file[working directory[option]]] nWhere option can be N
Definition: JMuonPostfit.sh:37
static const int JENERGY_MUON_RANGE_METRES
range of a muon with the reconstructed energy [m] from JEnergy.cc
void setLowerLimit(argument_type x)
Set lower limit.
Definition: JRange.hh:224
Maximum likelihood estimator (M-estimators).
bool qualitySorter(const JRECONSTRUCTION::JFit &first, const JRECONSTRUCTION::JFit &second)
Comparison of fit results.
static const int JMUONENERGY
then usage $script[input file[working directory[option]]] nWhere option can be E
Definition: JMuonPostfit.sh:37
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
static const int JPP_COVERAGE_POSITION
coverage of dynamic position calibration from any Jpp application