Jpp  17.1.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
software/JReconstruction/JEnergyCorrection.cc File Reference

Utility program to determine energy correction. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TMath.h"
#include "TFitResult.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/definitions/reconstruction.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JEnergyCorrection.hh"
#include "JGeometry3D/JTrack3E.hh"
#include "JPhysics/JGeane.hh"
#include "JTools/JQuantile.hh"
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JManager.hh"
#include "JROOT/JGraph.hh"
#include "JROOT/JGraphErrors.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "JLang/JToken.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Utility program to determine energy correction.

Author
mdejong

Definition in file software/JReconstruction/JEnergyCorrection.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 54 of file software/JReconstruction/JEnergyCorrection.cc.

55 {
56  using namespace std;
57  using namespace JPP;
58  using namespace KM3NETDAQ;
59 
60  typedef JTriggeredFileScanner<JEvt> JTriggeredFileScanner_t;
61  typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
62  typedef JToken<';'> JToken_t;
63  typedef JRange<double> JRange_t;
64 
65  JTriggeredFileScanner_t inputFile;
66  JLimit_t& numberOfEvents = inputFile.getLimit();
67  string outputFile;
68  JQualitySorter quality_sorter;
69  string formula;
70  vector<JToken_t> startValues;
71  vector<JToken_t> fixedValues;
72  JRange_t X;
73  string option;
74  string primary;
75  size_t numberOfBins;
76  int debug;
77 
78  try {
79 
80  JParser<> zap("Utility program to determine energy correction.");
81 
82  zap['f'] = make_field(inputFile, "event file(s) or histogram file");
83  zap['o'] = make_field(outputFile) = "energy_correction.root";
84  zap['n'] = make_field(numberOfEvents) = JLimit::max();
85  zap['L'] = make_field(quality_sorter) = JPARSER::initialised();
86  zap['F'] = make_field(formula, "fit formula, e.g: \"[0]+[1]*x\"") = "";
87  zap['@'] = make_field(startValues, "start values, e.g: \"p0 = GetMaximum;\"") = JPARSER::initialised();
88  zap['='] = make_field(fixedValues, "fixed values, e.g: \"p0 = GetMaximum;\"") = JPARSER::initialised();
89  zap['x'] = make_field(X, "abscissa range") = JRange_t();
90  zap['O'] = make_field(option, "Fit option") = "";
91  zap['p'] = make_field(primary) = muon_t, neutrino_t;
92  zap['N'] = make_field(numberOfBins,"number of bins") = 80;
93  zap['d'] = make_field(debug) = 2;
94 
95  zap(argc, argv);
96  }
97  catch(const exception& error) {
98  FATAL(error.what() << endl);
99  }
100 
101 
102  TFile* in = NULL;
103  TH2D* h2 = NULL;
104 
105  if (inputFile.size() == 1) {
106 
107  in = TFile::Open(inputFile[0].c_str(), "READ");
108 
109  if (in != NULL && in->IsOpen()) {
110  h2 = (TH2D*) in->Get("h2");
111  }
112  }
113 
114  if (h2 == NULL) {
115 
116  vector<JHead> enigma;
117 
118  for (JMultipleFileScanner_t::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
119 
120  try {
121 
122  const JHead head = getHeader(*i);
123 
124  vector<JHead>::iterator p = find(enigma.begin(), enigma.end(), head);
125 
126  if (p == enigma.end())
127  enigma.push_back(head);
128  else
129  p->add(head);
130  }
131  catch(const exception& error) {
132  FATAL(error.what());
133  }
134  }
135 
136  double Emin = numeric_limits<double>::max();
137  double Emax = 0.0;
138 
139  for (vector<JHead>::const_iterator i = enigma.begin(); i!= enigma.end(); ++i) {
140 
141  if (!i->is_valid(&JHead::cut_nu)) {
142  FATAL("Missing neutrino parameters." << endl);
143  }
144 
145  if (i->cut_nu.E.getLowerLimit() < Emin) { Emin = i->cut_nu.E.getLowerLimit(); }
146  if (i->cut_nu.E.getUpperLimit() > Emax) { Emax = i->cut_nu.E.getUpperLimit(); }
147 
148  if (i->is_valid(&JHead::genvol)) {
149  cout << "Neutrino energy range: " << i->cut_nu.E << endl;
150  cout << "Number of generated events: " << i->genvol.numberOfEvents << endl;
151  }
152  }
153 
154  const Double_t xmin = log10(Emin);
155  const Double_t xmax = log10(Emax);
156  const Double_t dx = (xmax - xmin) / numberOfBins;
157 
158  h2 = new TH2D("h2", NULL, numberOfBins + 4, xmin - 2*dx, xmax + 2*dx, numberOfBins, xmin, xmax);
159 
160  double W = 1.0;
161  JPosition3D center(0.0, 0.0, 0.0);
162 
163  for (string filename; inputFile.hasNext(); ) {
164 
165  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
166 
167  if (filename != inputFile.getFilename()) {
168 
169  filename = inputFile.getFilename();
170 
171  try {
172 
173  const JHead head = getHeader(filename);
174 
175  center = get<JPosition3D>(head);
176 
177  vector<JHead>::const_iterator p = find(enigma.begin(), enigma.end(), head);
178 
179  if (p->is_valid(&JHead::genvol)) {
180  W = 1.0 / (double) p->genvol.numberOfEvents;
181  }
182 
183  cout << "File: " << filename << " weight " << W << endl;
184  }
185  catch(const exception& error) {
186  FATAL(error.what());
187  }
188  }
189 
190  multi_pointer_type ps = inputFile.next();
191 
192  JEvt* evt = ps;
193  Evt* event = ps;
194 
195  Trk t1;
196 
197  if (primary == muon_t) {
198 
199  for (vector<Trk>::const_iterator i = event->mc_trks.begin(); i != event->mc_trks.end(); ++i) {
200  if (is_muon(*i)) {
201  if (i->E > t1.E) {
202  t1 = *i;
203  }
204  }
205  }
206 
207  if (!is_muon(t1)) {
208  continue;
209  }
210 
211  } else if (primary == neutrino_t) {
212 
213  if (has_neutrino(*event)) {
214  t1 = get_neutrino(*event);
215  }
216 
217  if (!is_neutrino(t1)) {
218  continue;
219  }
220  }
221 
222  JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_event(JMUONENERGY));
223 
224  if (evt->begin() == __end) {
225  continue;
226  }
227 
228  JEvt::iterator p = min_element(evt->begin(), __end, quality_sorter);
229 
230  JTrack3E ta = getTrack(t1);
231  JTrack3E tb = getTrack(*p);
232 
233  ta.add(center);
234 
235  tb.setE(p->getW(JENERGY_ENERGY));
236  tb.move(tb.getIntersection(ta.getPosition()), getSpeedOfLight(), gWater);
237 
238  if (ta.getE() > 0.0 && tb.getE() > 0.0) {
239  h2->Fill(log10(tb.getE()), log10(ta.getE()), (enigma.size() > (size_t) 1 ? W * event->w[2] : 1.0));
240  }
241  }
242  STATUS(endl);
243  }
244 
245  // short hands
246 
247  const Int_t nx = h2->GetXaxis()->GetNbins();
248  const Double_t xmin = h2->GetXaxis()->GetXmin();
249  const Double_t xmax = h2->GetXaxis()->GetXmax();
250 
251  const Int_t ny = h2->GetYaxis()->GetNbins();
252 
253  TH1D h1("h1", NULL, nx, xmin, xmax);
254 
255  JGraph_t g1;
256  JGraph_t g2;
257  JGraphErrors_t g3;
258 
259  JManager<int, TH1D> H0(new TH1D("py[%]", NULL, nx, xmin, xmax));
260 
261  for (Int_t ix = 1; ix <= nx; ++ix) {
262 
263  TH1D* h0 = H0[ix];
264 
265  JQuantile Q("", true); // start values
266 
267  for (int iy = 1; iy <= ny; ++iy) {
268 
269  const Double_t x = h2->GetYaxis()->GetBinCenter(iy);
270  const Double_t y = h2->GetBinContent(ix,iy);
271  const Double_t z = h2->GetBinError (ix,iy);
272 
273  if (y > 0.0) {
274 
275  h0->SetBinContent(iy, y);
276  h0->SetBinError (iy, z);
277 
278  Q.put(x,y);
279  }
280  }
281 
282  if (Q.getCount() >= 3) {
283 
284  if (debug >= debug_t) {
285  Q.print(cout);
286  }
287 
288  TF1 f1("f1", "[0]*TMath::Gaus(x,[1],[2]) + [3]");
289 
290  const double x0 = Q.getQuantile(0.50);
291  const double sigma = 0.5 * (Q.getQuantile(0.66) - Q.getQuantile(0.33));
292 
293  f1.SetParameter(0, Q.getWmax());
294  f1.SetParameter(1, x0);
295  f1.SetParameter(2, sigma);
296  f1.SetParameter(3, 0.0);
297 
298  TFitResultPtr result = h0->Fit(&f1, "SQ", "same");
299 
300  if (result.Get() == NULL) {
301  FATAL("Invalid TFitResultPtr " << h0->GetName() << endl);
302  }
303 
304  if (debug >= debug_t || !result->IsValid()) {
305  cout << "Bin: "
306  << setw(4) << ix << ' '
307  << FIXED(7,3) << h2->GetXaxis()->GetBinCenter(ix) << ' '
308  << FIXED(7,3) << result->Chi2() << '/'
309  << result->Ndf() << ' '
310  << (result->IsValid() ? "" : "failed") << endl;
311  }
312 
313  if (result->IsValid()) {
314 
315  const Double_t x = h2->GetXaxis()->GetBinCenter(ix);
316 
317  h1.SetBinContent(ix, f1.GetParameter(1));
318  h1.SetBinError (ix, f1.GetParError (1));
319 
320  g1.put(x, pow(10.0, x) - pow(10.0, f1.GetParameter(1)));
321  g2.put(x, pow(10.0, abs(f1.GetParameter(2)) - f1.GetParameter(1)));
322  g3.put(x, f1.GetParameter(1), 0.0, f1.GetParameter(2));
323  }
324  }
325  }
326 
327  TF1* f1 = NULL;
328 
329  if (formula != "") { // fit
330 
331  f1 = new TF1("f1", formula.c_str());
332 
333  if (f1->IsZombie()) {
334  FATAL("Function: " << formula << " is zombie." << endl);
335  }
336 
337  try {
338 
339  for (vector<JToken_t>::const_iterator i = startValues.begin(); i != startValues.end(); ++i) {
340  f1->SetParameter(getParameter(*i), getValue(*i,&h1));
341  }
342 
343  for (vector<JToken_t>::const_iterator i = fixedValues.begin(); i != fixedValues.end(); ++i) {
344  f1->FixParameter(getParameter(*i), getValue(*i,&h1));
345  }
346  }
347  catch(JLANG::JParseError& error) {
348  FATAL(error << endl);
349  }
350 
351  if (option.find('S') == string::npos) { option += "S"; }
352 
353  TFitResultPtr result = h1.Fit(f1, option.c_str(), "same", X.getLowerLimit(), X.getUpperLimit());
354 
355  } else { // interpolate
356 
357  Double_t x[2] = { xmin, xmax };
358  Double_t y[2] = { xmin, xmax };
359  Double_t z[2] = { 0.00, 0.00 };
360 
361  ostringstream os;
362 
363  os << " " // extrapolate
364  << "("
365  << "x >= " << FIXED(5,3) << xmax
366  << ")"
367  << " * "
368  << "(x)";
369  os << "\n";
370 
371  for (Int_t ix = nx; ix >= 1; --ix) {
372 
373  x[0] = h1.GetXaxis()->GetBinCenter(ix);
374  y[0] = h1.GetBinContent(ix);
375  z[0] = h1.GetBinError (ix);
376 
377  if (z[0] == 0.0 || !X(x[0])) {
378  continue;
379  }
380 
381  os << " + "
382  << "("
383  << "x >= " << FIXED(5,3) << x[0]
384  << " && "
385  << "x < " << FIXED(5,3) << x[1]
386  << ")"
387  << " * "
388  << "(" << FIXED(5,3) << y[0] << " + " << FIXED(5,3) << (y[1] - y[0]) << " * (x - " << FIXED(5,3) << x[0] << ") / " << FIXED(5,3) << (x[1] - x[0]) << ")";
389  os << "\n";
390 
391  x[1] = x[0];
392  y[1] = y[0];
393  z[1] = z[0];
394  }
395 
396  x[0] = -1;
397  y[0] = -1;
398 
399  os << " + " // extrapolate
400  << "("
401  << "x < " << FIXED(5,3) << x[1]
402  << ")"
403  << " * "
404  << "(" << FIXED(5,3) << y[0] << " + " << FIXED(5,3) << (y[1] - y[0]) << " * (x - " << FIXED(5,3) << x[0] << ") / " << FIXED(5,3) << (x[1] - x[0]) << ")";
405  os << "\n";
406 
407  string buffer = os.str();
408 
409  NOTICE(buffer << endl);
410 
411  for (string::iterator i = buffer.begin(); i != buffer.end(); ++i) {
412  if (*i == '\n') {
413  *i = ' ';
414  }
415  }
416 
417  f1 = new TF1("f1", buffer.c_str(), xmin, xmax);
418 
419  h1.GetListOfFunctions()->Add(f1);
420  }
421 
422 
423  TFile out(outputFile.c_str(), "recreate");
424 
425  out << JMeta(argc, argv);
426 
427  out << *h2 << H0 << h1;
428 
429  out << JGraph (g1, "g1")
430  << JGraph (g2, "g2")
431  << JGraphErrors(g3, "g3");
432 
433  if (f1 != NULL) { // separately store energy correction
434 
435  f1->GetFormula()->SetName(JEnergyCorrection::getName());
436 
437  out << *(f1->GetFormula());
438  }
439 
440  out.Write();
441  out.Close();
442 
443  if (in != NULL) {
444  in->Close();
445  }
446 }
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
GAUSS_LEGENDRE::xmax
const double xmax
Definition: JQuadrature.cc:24
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1500
Q
Q(UTCMax_s-UTCMin_s)-livetime_s
JEEP::getValue
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
JEEP::debug_t
debug
Definition: JMessage.hh:29
JGIZMO::getParameter
int getParameter(const std::string &text)
Get parameter number from text string.
Definition: JGizmoToolkit.hh:332
JAANET::getTrack
JTrack3E getTrack(const Trk &track)
Get track.
Definition: JAAnetToolkit.hh:258
JENERGY_ENERGY
static const int JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
Definition: fitparameters.hh:16
JROOT::JGraph_t
Data structure for graph data.
Definition: JGraph.hh:21
mc_keys::neutrino_t
const char *const neutrino_t
Definition: io_ascii.hh:29
JTOOLS::JQuantile
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JAANET::has_neutrino
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Definition: JAAnetToolkit.hh:435
JGEOMETRY3D::JAxis3D::getIntersection
double getIntersection(const JVector3D &pos) const
Get longitudinal position along axis of position of closest approach with given position.
Definition: JAxis3D.hh:146
JAANET::is_muon
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Definition: JAAnetToolkit.hh:338
JRECONSTRUCTION::JQualitySorter
General purpose sorter of fit results.
Definition: JReconstruction/JEvtToolkit.hh:260
JSUPPORT::JTriggeredFileScanner
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
Definition: JTriggeredFileScanner.hh:39
JROOT::JGraphErrors_t
Data structure for graph data.
Definition: JGraphErrors.hh:21
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
JPHYSICS::gWater
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
FIXED
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
JROOT::JGraph
Auxiliary data structure to build TGraph.
Definition: JGraph.hh:42
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
Trk::E
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
JGEOMETRY3D::JTrack3E
3D track with energy.
Definition: JTrack3E.hh:30
JGEOMETRY3D::JTrack3E::move
void move(const double step, const double velocity, const JGeane &geane)
Move vertex along this track with given velocity.
Definition: JTrack3E.hh:117
GAUSS_HERMITE::sigma
const double sigma[]
Definition: JQuadrature.cc:74
JSUPPORT::getHeader
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Definition: JMonteCarloFileSupportkit.hh:491
JAANET::is_neutrino
bool is_neutrino(const Trk &track)
Test whether given track is a neutrino.
Definition: JAAnetToolkit.hh:320
GAUSS_LEGENDRE::f1
const JPolynome f1(1.0, 2.0, 3.0)
Function.
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
JAANET::JRange_t
Type definition of range.
Definition: JHead.hh:41
JROOT::JManager
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
JGEOMETRY3D::JTrack3E::getE
double getE() const
Get energy.
Definition: JTrack3E.hh:93
JACOUSTICS::JEvt
Acoustic event fit.
Definition: JAcoustics/JEvt.hh:208
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
log10
set_variable E_E log10(E_{fit}/E_{#mu})"
JFIT::JHistory::is_event
Auxiliary class to test history.
Definition: JHistory.hh:105
JTOOLS::result
return result
Definition: JPolint.hh:743
JGEOMETRY3D::JTrack3E::setE
void setE(const double E)
Set energy.
Definition: JTrack3E.hh:104
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JMATH::pow
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:98
debug
int debug
debug level
Definition: JSirene.cc:67
quality_sorter
Reconstruction type dependent comparison of track quality.
Definition: tools/reconstruction.hh:85
JAANET::JHead
Monte Carlo run header.
Definition: JHead.hh:1167
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
GAUSS_LEGENDRE::xmin
const double xmin
Definition: JQuadrature.cc:23
JROOT::JGraphErrors
Auxiliary data structure to build TGraphErrors.
Definition: JGraphErrors.hh:49
JPHYSICS::getSpeedOfLight
const double getSpeedOfLight()
Get speed of light.
Definition: JPhysics/JConstants.hh:106
JLANG::JToken
Wrapper class around string.
Definition: JToken.hh:23
JAANET::primary
Primary particle.
Definition: JHead.hh:1107
JLANG::JParseError
Exception for parsing value.
Definition: JException.hh:180
X
no fit printf nominal n $STRING awk v X
Definition: detector-XY:fit.sh:102
JROOT::getName
const char * getName()
Get ROOT name of given data type.
Definition: JRootToolkit.hh:57
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
JAANET::get_neutrino
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
Definition: JAAnetToolkit.hh:446
in
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 JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:46
Trk
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
numberOfBins
int numberOfBins
number of bins for average CDF integral of optical module
Definition: JSirene.cc:68
JMUONENERGY
static const int JMUONENERGY
Definition: definitions/reconstruction.hh:17
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
g1
Double_t g1(const Double_t x)
Function.
Definition: JQuantiles.cc:25
Generated by   doxygen 1.8.5