Jpp 20.0.0-27-g39925593c-D
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/JStats.hh"
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JManager.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 49 of file software/JReconstruction/JEnergyCorrection.cc.

50{
51 using namespace std;
52 using namespace JPP;
53 using namespace KM3NETDAQ;
54
55 typedef JTriggeredFileScanner<JEvt> JTriggeredFileScanner_t;
56 typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
57 typedef JToken<';'> JToken_t;
58 typedef JRange<double> JRange_t;
59
60 JTriggeredFileScanner_t inputFile;
61 JLimit_t& numberOfEvents = inputFile.getLimit();
62 string outputFile;
63 JQualitySorter quality_sorter;
64 string formula;
65 vector<JToken_t> startValues;
66 vector<JToken_t> fixedValues;
67 JRange_t X;
68 string option;
69 size_t numberOfBins;
70 int debug;
71
72 try {
73
74 JParser<> zap("Utility program to determine energy correction.");
75
76 zap['f'] = make_field(inputFile, "event file(s) or histogram file");
77 zap['o'] = make_field(outputFile) = "energy_correction.root";
78 zap['n'] = make_field(numberOfEvents) = JLimit::max();
79 zap['L'] = make_field(quality_sorter) = JPARSER::initialised();
80 zap['F'] = make_field(formula, "fit formula, e.g: \"[0]+[1]*x\"") = "";
81 zap['@'] = make_field(startValues, "start values, e.g: \"p0 = GetMaximum;\"") = JPARSER::initialised();
82 zap['='] = make_field(fixedValues, "fixed values, e.g: \"p0 = GetMaximum;\"") = JPARSER::initialised();
83 zap['x'] = make_field(X, "abscissa range") = JRange_t();
84 zap['O'] = make_field(option, "Fit option") = "";
85 zap['N'] = make_field(numberOfBins, "number of bins") = 80;
86 zap['d'] = make_field(debug) = 2;
87
88 zap(argc, argv);
89 }
90 catch(const exception& error) {
91 FATAL(error.what() << endl);
92 }
93
94
95 const double NUMBER_OF_ENTRIES = 100.0; // minimum number of histogram entries
96 const double THRESHOLD = 0.05; // minimal bin content relative to maximal bin content
97
98 TFile* in = NULL;
99 TH2D* h2 = NULL;
100
101 if (inputFile.size() == 1) {
102
103 in = TFile::Open(inputFile[0].c_str(), "READ");
104
105 if (in != NULL && in->IsOpen()) {
106 h2 = (TH2D*) in->Get("h2");
107 }
108 }
109
110 if (h2 == NULL) {
111
112 const Double_t xmin = X.getLowerLimit();
113 const Double_t xmax = X.getUpperLimit();
114 const Double_t dx = (xmax - xmin) / numberOfBins;
115
116 h2 = new TH2D("h2", NULL, numberOfBins + 4, xmin - 2*dx, xmax + 2*dx, 3*numberOfBins, xmin, xmax);
117
118 while (inputFile.hasNext()) {
119
120 STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
121
122 multi_pointer_type ps = inputFile.next();
123
124 JEvt* evt = ps;
125 Evt* event = ps;
126
127 Trk t1;
128
129 for (vector<Trk>::const_iterator i = event->mc_trks.begin(); i != event->mc_trks.end(); ++i) {
130 if (is_muon(*i)) {
131 if (i->E > t1.E) {
132 t1 = *i;
133 }
134 }
135 }
136
137 if (!is_muon(t1)) {
138 continue;
139 }
140
141 JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_application(JMUONENERGY));
142
143 if (evt->begin() == __end) {
144 continue;
145 }
146
147 JEvt::iterator p = min_element(evt->begin(), __end, quality_sorter);
148
149 JTrack3E ta = getTrack(t1);
150 JTrack3E tb = getTrack(*p);
151
152 tb.setE(p->getW(JENERGY_ENERGY));
153 tb.move(tb.getIntersection(ta.getPosition()), getSpeedOfLight(), gWater);
154
155 if (ta.getE() > 0.0 && tb.getE() > 0.0) {
156 h2->Fill(log10(tb.getE()), log10(ta.getE()));
157 }
158 }
159 STATUS(endl);
160 }
161
162 // short hands
163
164 const Int_t nx = h2->GetXaxis()->GetNbins();
165 const Double_t xmin = h2->GetXaxis()->GetXmin();
166 const Double_t xmax = h2->GetXaxis()->GetXmax();
167 const Int_t ny = h2->GetYaxis()->GetNbins();
168 const Double_t ymin = h2->GetYaxis()->GetXmin();
169 const Double_t ymax = h2->GetYaxis()->GetXmax();
170
171 TH1D h1("h1", NULL, nx, xmin, xmax);
172
173 JManager<int, TH1D> H0(new TH1D("py[%]", NULL, ny, ymin, ymax));
174
175 for (Int_t ix = 1; ix <= nx; ++ix) {
176
177 TH1D* h0 = H0[ix];
178
179 // start values
180
181 double x0 = 0.0;
182 double y0 = 0.0;
183
184 for (int iy = 1; iy <= ny; ++iy) {
185
186 const Double_t x = h2->GetYaxis()->GetBinCenter(iy);
187 const Double_t y = h2->GetBinContent(ix,iy);
188 const Double_t z = h2->GetBinError (ix,iy);
189
190 if (y > 0.0) {
191 h0->SetBinContent(iy, y);
192 h0->SetBinError (iy, z);
193 }
194
195 if (y >= y0) {
196 x0 = x;
197 y0 = y;
198 }
199 }
200
201 JStats Q;
202
203 for (int iy = 1 ; iy <= ny; ++iy) {
204
205 const Double_t x = h0->GetXaxis()->GetBinCenter(iy);
206 const Double_t y = h0->GetBinContent(iy);
207
208 if (y > THRESHOLD * y0) {
209 Q.put(x,y);
210 }
211 }
212
213 if (Q.getCount() >= 5 && Q.getTotal() >= NUMBER_OF_ENTRIES) {
214
215 if (debug >= debug_t) {
216 Q.print(cout);
217 }
218
219 TF1 f1("f1", "(x < [1]) * [0]*TMath::Gaus(x,[1],[2]) + (x >= [1]) * ([0] * (1.0 - [3])*TMath::Gaus(x,[1],[2]) + [3]*[0])");
220
221 const double sigma = 0.7 * Q.getSTDev();
222
223 f1.SetParameter(0, y0);
224 f1.SetParameter(1, x0);
225 f1.SetParameter(2, sigma);
226 f1.SetParameter(3, 0.05);
227 f1.SetParLimits(3, 0.0, 1.0);
228
229 TFitResultPtr result = h0->Fit(&f1, "SQL", "same", x0 - 2.5*sigma, x0 + 2.5*sigma);
230
231 if (result.Get() == NULL) {
232 FATAL("Invalid TFitResultPtr " << h0->GetName() << endl);
233 }
234
235 if (debug >= status_t || !result->IsValid()) {
236 cout << "Bin: "
237 << setw(3) << ix << ' '
238 << FIXED(6,3) << h2->GetXaxis()->GetBinCenter(ix) << ' '
239 << FIXED(6,3) << x0 << " -> "
240 << FIXED(6,3) << f1.GetParameter(1) << ' '
241 << FIXED(6,3) << sigma << " -> "
242 << FIXED(6,3) << f1.GetParameter(2) << ' '
243 << FIXED(7,3) << result->Chi2() << '/'
244 << result->Ndf() << ' '
245 << (result->IsValid() ? "" : "failed") << endl;
246 }
247
248 if (result->IsValid()) {
249 h1.SetBinContent(ix, f1.GetParameter(1));
250 h1.SetBinError (ix, f1.GetParError (1));
251 }
252 }
253 }
254
255 TF1* f1 = NULL;
256
257 if (formula != "") { // fit
258
259 f1 = new TF1("f1", formula.c_str());
260
261 if (f1->IsZombie()) {
262 FATAL("Function: " << formula << " is zombie." << endl);
263 }
264
265 try {
266
267 for (vector<JToken_t>::const_iterator i = startValues.begin(); i != startValues.end(); ++i) {
268 f1->SetParameter(getParameter(*i), getValue(*i,&h1));
269 }
270
271 for (vector<JToken_t>::const_iterator i = fixedValues.begin(); i != fixedValues.end(); ++i) {
272 f1->FixParameter(getParameter(*i), getValue(*i,&h1));
273 }
274 }
275 catch(JLANG::JParseError& error) {
276 FATAL(error << endl);
277 }
278
279 if (option.find('S') == string::npos) { option += "S"; }
280
281 TFitResultPtr result = h1.Fit(f1, option.c_str(), "same", X.getLowerLimit(), X.getUpperLimit());
282
283 } else { // interpolate
284
285 Double_t x[2] = { xmin, xmax };
286 Double_t y[2] = { xmin, xmax };
287 Double_t z[2] = { 0.00, 0.00 };
288
289 ostringstream os;
290
291 os << " " // extrapolate
292 << "("
293 << "x >= " << FIXED(5,3) << xmax
294 << ")"
295 << " * "
296 << "(x)";
297 os << "\n";
298
299 for (Int_t ix = nx; ix >= 1; --ix) {
300
301 x[0] = h1.GetXaxis()->GetBinCenter(ix);
302 y[0] = h1.GetBinContent(ix);
303 z[0] = h1.GetBinError (ix);
304
305 if (z[0] == 0.0 || !X(x[0])) {
306 continue;
307 }
308
309 os << " + "
310 << "("
311 << "x >= " << FIXED(5,3) << x[0]
312 << " && "
313 << "x < " << FIXED(5,3) << x[1]
314 << ")"
315 << " * "
316 << "(" << 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]) << ")";
317 os << "\n";
318
319 x[1] = x[0];
320 y[1] = y[0];
321 z[1] = z[0];
322 }
323
324 x[0] = -1;
325 y[0] = -1;
326
327 os << " + " // extrapolate
328 << "("
329 << "x < " << FIXED(5,3) << x[1]
330 << ")"
331 << " * "
332 << "(" << 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]) << ")";
333 os << "\n";
334
335 string buffer = os.str();
336
337 NOTICE(buffer << endl);
338
339 for (string::iterator i = buffer.begin(); i != buffer.end(); ++i) {
340 if (*i == '\n') {
341 *i = ' ';
342 }
343 }
344
345 f1 = new TF1("f1", buffer.c_str(), xmin, xmax);
346
347 h1.GetListOfFunctions()->Add(f1);
348 }
349
350
351 TFile out(outputFile.c_str(), "recreate");
352
353 out << JMeta(argc, argv);
354
355 out << *h2 << H0 << h1;
356
357 if (f1 != NULL) { // separately store energy correction
358
359 f1->GetFormula()->SetName(JEnergyCorrection::getName());
360
361 out << *(f1->GetFormula());
362 }
363
364 out.Write();
365 out.Close();
366
367 if (in != NULL) {
368 in->Close();
369 }
370}
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:72
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
numberOfBins
int numberOfBins
number of bins for average CDF integral of optical module
Definition JSirene.cc:73
JGEOMETRY3D::JTrack3E
3D track with energy.
Definition JTrack3E.hh:34
JLANG::JParseError
Exception for parsing value.
Definition JException.hh:198
JLANG::JToken
Wrapper class around string.
Definition JToken.hh:26
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JRECONSTRUCTION::JEnergyCorrection::getName
static const char * getName()
Get name of energy correction formula.
Definition JEnergyCorrection.hh:144
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
JENERGY_ENERGY
static const int JENERGY_ENERGY
uncorrected energy [GeV] see JRECONSTRUCTION::JMuonEnergy
Definition fitparameters.hh:16
JAANET::getTrack
JTrack3E getTrack(const Trk &track)
Get track.
Definition JAAnetToolkit.hh:278
JAANET::is_muon
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Definition JAAnetToolkit.hh:366
JEEP::getValue
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition JScale.hh:47
JEEP::status_t
@ status_t
status
Definition JMessage.hh:30
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:387
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
JTOOLS::result
return result
Definition JPolint.hh:862
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
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
JACOUSTICS::JEvt
Acoustic event fit.
Definition JAcoustics/JEvt.hh:310
JFIT::JHistory::is_application
Auxiliary class to test history.
Definition JHistory.hh:188
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:426
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
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
JTOOLS::JStats
Auxiliary data structure for running average and standard deviation.
Definition JStats.hh:44
JTOOLS::JStats::getCount
long long int getCount() const
Get total count.
Definition JStats.hh:168
JTOOLS::JStats::getTotal
double getTotal() const
Get total weight.
Definition JStats.hh:179
JTOOLS::JStats::print
std::ostream & print(std::ostream &out, bool lpr=true) const
Print stats.
Definition JStats.hh:320
JTOOLS::JStats::put
void put(const double x, const double w=1.0)
Put value.
Definition JStats.hh:119
JTOOLS::JStats::getSTDev
double getSTDev() const
Get standard deviation.
Definition JStats.hh:263
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
quality_sorter
Reconstruction type dependent comparison of track quality.
Definition tools/reconstruction.hh:85
Generated by doxygen 1.10.0