50{
54
56 typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
57 typedef JToken<
';'> JToken_t;
59
60 JTriggeredFileScanner_t inputFile;
64 string formula;
68 string option;
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");
80 zap[
'F'] =
make_field(formula,
"fit formula, e.g: \"[0]+[1]*x\"") =
"";
84 zap[
'O'] =
make_field(option,
"Fit option") =
"";
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;
96 const double THRESHOLD = 0.05;
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
115
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
126
128
129 for (vector<Trk>::const_iterator i = event->mc_trks.begin(); i != event->mc_trks.end(); ++i) {
132 t1 = *i;
133 }
134 }
135 }
136
138 continue;
139 }
140
142
143 if (evt->begin() == __end) {
144 continue;
145 }
146
147 JEvt::iterator p = min_element(evt->begin(), __end,
quality_sorter);
148
151
154
155 if (ta.
getE() > 0.0 && tb.
getE() > 0.0) {
156 h2->Fill(log10(tb.
getE()), log10(ta.
getE()));
157 }
158 }
160 }
161
162
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
174
175 for (Int_t ix = 1; ix <= nx; ++ix) {
176
177 TH1D* h0 = H0[ix];
178
179
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) {
198 }
199 }
200
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) {
210 }
211 }
212
214
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
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
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) <<
' '
242 <<
FIXED(6,3) <<
f1.GetParameter(2) <<
' '
245 << (
result->IsValid() ?
"" :
"failed") << endl;
246 }
247
249 h1.SetBinContent(ix,
f1.GetParameter(1));
250 h1.SetBinError (ix,
f1.GetParError (1));
251 }
252 }
253 }
254
256
257 if (formula != "") {
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) {
269 }
270
271 for (vector<JToken_t>::const_iterator i = fixedValues.begin(); i != fixedValues.end(); ++i) {
273 }
274 }
276 FATAL(error << endl);
277 }
278
279 if (option.find('S') == string::npos) { option += "S"; }
280
282
283 } else {
284
287 Double_t z[2] = { 0.00, 0.00 };
288
289 ostringstream os;
290
291 os << " "
292 << "("
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
321 z[1] = z[0];
322 }
323
326
327 os << " + "
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
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
352
353 out <<
JMeta(argc, argv);
354
355 out << *h2 << H0 << h1;
356
357 if (f1 != NULL) {
358
360
361 out << *(
f1->GetFormula());
362 }
363
364 out.Write();
365 out.Close();
366
367 if (in != NULL) {
368 in->Close();
369 }
370}
#define DEBUG(A)
Message macros.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
int numberOfBins
number of bins for average CDF integral of optical module
double getIntersection(const JVector3D &pos) const
Get longitudinal position along axis of position of closest approach with given position.
const JPosition3D & getPosition() const
Get position.
void move(const double step, const double velocity, const JGeane &geane)
Move vertex along this track with given velocity.
void setE(const double E)
Set energy.
double getE() const
Get energy.
Exception for parsing value.
Wrapper class around string.
Utility class to parse command line options.
static const char * getName()
Get name of energy correction formula.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
static const int JMUONENERGY
static const int JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
const JPolynome f1(1.0, 2.0, 3.0)
Function.
JTrack3E getTrack(const Trk &track)
Get track.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
int getParameter(const std::string &text)
Get parameter number from text string.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
KM3NeT DAQ data structures and auxiliaries.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Auxiliary data structure for floating point format specification.
Type definition of range.
Auxiliary class to test history.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
General purpose sorter of fit results.
Auxiliary class for defining the range of iterations of objects.
const JLimit & getLimit() const
Get limit.
static counter_type max()
Get maximum counter value.
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
double E
Energy [GeV] (either MC truth or reconstructed)
Reconstruction type dependent comparison of track quality.