Jpp master_rocky-44-g75b7c4f75
the software that should make you happy
Loading...
Searching...
No Matches
Functions
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 "JROOT/JMinimizer.hh"
#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

◆ main()

int main ( int argc,
char ** argv )

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

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