Jpp 20.0.0-27-g39925593c-D
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JHobbit.cc File Reference

Program to fit time-residuals histogram in output of JCalibrateMuon.cc. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#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 "JLang/JLangToolkit.hh"
#include "JSupport/JMeta.hh"
#include "JTools/JRange.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTimeRange.hh"
#include "JDetector/JStringRouter.hh"
#include "JDetector/JLocationRouter.hh"
#include "JCalibrate/JCalibrateMuon.hh"
#include "Jeep/JPrint.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

Program to fit time-residuals histogram in output of JCalibrateMuon.cc.

Author
mdejong

Definition in file JHobbit.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 52 of file JHobbit.cc.

53{
54 using namespace std;
55 using namespace JPP;
56
57 vector<string> inputFile;
58 string outputFile;
59 string detectorFile;
60 bool overwriteDetector;
61 string function;
62 JTimeRange T_ns;
63 string option;
64 JTimeRange E_ns;
65 double WMin;
66 int ID;
67 string target;
68 int debug;
69
70 try {
71
72 JParser<> zap("Program to fit time-residuals histogram in output of JCalibrateMuon.cc.");
73
74 zap['f'] = make_field(inputFile, "input files (output from JCalibrateMuon).");
75 zap['o'] = make_field(outputFile, "output file.") = "hobbit.root";
76 zap['a'] = make_field(detectorFile, "detector file.");
77 zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with correct time offsets.");
78 zap['F'] = make_field(function, "fit function") = gauss_t, landau_t, emg_t, breitwigner_t;
79 zap['T'] = make_field(T_ns, "ROOT fit range around maximum [ns].") = JTimeRange(-7.5,+5.0);
80 zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "";
81 zap['E'] = make_field(E_ns, "validity range of t0 uncertainty [ns].") = JTimeRange( 0.0, 0.5);
82 zap['W'] = make_field(WMin, "minimal histogram contents.") = 100.0;
83 zap['D'] = make_field(ID, "target identifier") = -1;
84 zap['R'] = make_field(target, "option for time calibration.") = module_t, string_t;
85 zap['d'] = make_field(debug) = 1;
86
87 zap(argc, argv);
88 }
89 catch(const exception &error) {
90 FATAL(error.what() << endl);
91 }
92
93
94 if (!T_ns.is_valid()) {
95 FATAL("Invalid fit range [ns] " << T_ns << endl);
96 }
97
98 JDetector detector;
99
100 try {
101 load(detectorFile, detector);
102 }
103 catch(const JException& error) {
104 FATAL(error);
105 }
106
107 const JLocationRouter router(detector);
108
109 const JRouter_t* rpm = NULL;
110
111 if (target == module_t)
112 rpm = new JModuleRouter_t(detector);
113 else if (target == string_t)
114 rpm = new JStringRouter_t(detector);
115 else
116 FATAL("No valid target; check option -R" << endl);
117
118 if (option.find('S') == string::npos) { option += 'S'; }
119 if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
120
121 TH2D* h2 = NULL;
122
123 for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
124
125 NOTICE("Processing " << *i << endl);
126
127 TFile in(i->c_str(), "exist");
128
129 if (!in.IsOpen()) {
130 FATAL("File " << *i << " not opened." << endl);
131 }
132
133 TH2D* p = dynamic_cast<TH2D*>(in.Get(h2_t));
134
135 if (p == NULL) {
136 FATAL("File " << *i << " has no histogram <" << h2_t << ">." << endl);
137 }
138
139 if (h2 == NULL)
140 h2 = (TH2D*) p->Clone();
141 else
142 h2->Add(p);
143
144 h2->SetDirectory(0);
145
146 in.Close();
147 }
148
149 if (h2 == NULL) {
150 FATAL("No histogram <" << h2_t << ">." << endl);
151 }
152
153 // auxiliary data structure for fit result
154
155 struct result_type {
156
157 result_type() :
158 value(0.0),
159 error(0.0)
160 {}
161
162 result_type(double value,
163 double error) :
164 value(value),
165 error(error)
166 {}
167
168 double value;
169 double error;
170 };
171
172 typedef map<int, result_type> map_type; // index -> fit result
173
174 map_type zmap;
175
176
177 // histogram fit results
178
179 const TAxis* x_axis = h2->GetXaxis();
180 const TAxis* y_axis = h2->GetYaxis();
181
182 TH1D h0("h0", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
183 TH1D hc("hc", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
184 TH1D hq("hq", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
185
186 const floor_range floors = getRangeOfFloors(detector);
187 const JStringRouter strings(detector);
188
189 TH2D H2("detector", NULL,
190 strings.size(), -0.5, strings.size() - 0.5,
191 floors.getUpperLimit(), 1 - 0.5, floors.getUpperLimit() + 0.5);
192
193 for (Int_t ix = 1; ix <= H2.GetXaxis()->GetNbins(); ++ix) {
194 H2.GetXaxis()->SetBinLabel(ix, MAKE_CSTRING(strings.at(ix-1)));
195 }
196 for (Int_t iy = 1; iy <= H2.GetYaxis()->GetNbins(); ++iy) {
197 H2.GetYaxis()->SetBinLabel(iy, MAKE_CSTRING(iy));
198 }
199
200 TFile out(outputFile.c_str(), "recreate");
201
202 size_t counts = 0;
203 size_t errors = 0;
204
205 for (Int_t ix = 1; ix <= x_axis->GetNbins(); ++ix) {
206
207 const int index = ix - 1;
208 const int id = rpm->getID(index);
209
210 DEBUG("Bin " << setw(4) << ix << " -> " << setw(8) << id << endl);
211
212 if (ID != -1 && ID != id) {
213 continue;
214 }
215
216 TH1D h1(MAKE_CSTRING(id << ".1D"), NULL, y_axis->GetNbins(), y_axis->GetXmin(), y_axis->GetXmax());
217
218 // start values
219
220 Double_t ymax = 0.0;
221 Double_t x0 = 0.0; // [ns]
222 Double_t sigma = 4.0; // [ns]
223 Double_t W = 0.0;
224
225 for (Int_t iy = 1; iy <= y_axis->GetNbins(); ++iy) {
226
227 h1.SetBinContent(iy, h2->GetBinContent(ix,iy));
228 h1.SetBinError (iy, sqrt(h2->GetBinContent(ix,iy)));
229
230 const Double_t x = h1.GetBinCenter (iy);
231 const Double_t y = h1.GetBinContent(iy);
232
233 W += y;
234
235 if (y > ymax) {
236 ymax = y;
237 x0 = x;
238 }
239 }
240
241 if (W <= WMin) {
242
243 WARNING("Not enough entries for slice " << ix << ' ' << W << "; skip fit." << endl);
244
245 continue;
246 }
247
248 ymax *= 0.9;
249
250 const Double_t xmin = x0 + T_ns.getLowerLimit();
251 const Double_t xmax = x0 + T_ns.getUpperLimit();
252
253
254 TF1* f1 = NULL;
255
256 if (function == gauss_t) {
257
258 f1 = new TF1(function.c_str(), "[0]*TMath::Gaus(x,[1],[2])");
259
260 f1->SetParameter(0, 0.8*ymax);
261 f1->SetParameter(1, x0);
262 f1->SetParameter(2, sigma);
263
264 f1->SetParError(0, sqrt(ymax) * 0.1);
265 f1->SetParError(1, 0.01);
266 f1->SetParError(2, 0.01);
267
268 } else if (function == landau_t) {
269
270 f1 = new TF1(function.c_str(), "[0]*TMath::Landau(x,[1],[2])");
271
272 f1->SetParameter(0, 0.8*ymax);
273 f1->SetParameter(1, x0);
274 f1->SetParameter(2, sigma);
275
276 f1->SetParError(0, sqrt(ymax) * 0.1);
277 f1->SetParError(1, 0.01);
278 f1->SetParError(2, 0.01);
279
280 } else if (function == emg_t) {
281
282 f1 = new TF1(function.c_str(), "[0]*TMath::Exp(0.5*[3]*(2.0*[1]+[3]*[2]*[2]-2.0*x))*TMath::Erfc(([1]+[3]*[2]*[2]-x)/(TMath::Sqrt(2.0)*[2]))");
283
284 f1->SetParameter(0, 0.5*ymax);
285 f1->SetParameter(1, x0 - sigma);
286 f1->SetParameter(2, sigma);
287 f1->SetParameter(3, 0.06);
288
289 f1->SetParError(0, sqrt(ymax) * 0.1);
290 f1->SetParError(1, 0.01);
291 f1->SetParError(2, 0.01);
292 f1->SetParError(3, 1.0e-4);
293
294 } else if (function == breitwigner_t) {
295
296 f1 = new TF1(function.c_str(), "(x <= [1])*[0]*[2]*TMath::BreitWigner(x,[1],[2])+(x > [1])*[0]*[3]*TMath::BreitWigner(x,[1],[3])");
297
298 f1->SetParameter(0, 0.8*ymax);
299 f1->SetParameter(1, x0);
300 f1->SetParameter(2, 15.0);
301 f1->SetParameter(3, 25.0);
302
303 f1->SetParError(0, sqrt(ymax) * 0.1);
304 f1->SetParError(1, 0.01);
305 f1->SetParError(2, 0.1);
306 f1->SetParError(3, 0.1);
307
308 } else {
309
310 FATAL("Unknown fit function " << function << endl);
311 }
312
313
314 DEBUG("Start values:" << endl);
315
316 for (int i = 0; i != f1->GetNpar(); ++i) {
317 DEBUG(left << setw(12) << f1->GetParName (i) << ' ' <<
318 SCIENTIFIC(12,5) << f1->GetParameter(i) << endl);
319 }
320
321 // fit
322
323 TFitResultPtr result = h1.Fit(f1, option.c_str(), "same", xmin, xmax);
324
325 const bool status = result->IsValid() && E_ns(f1->GetParError(1));
326
327 if (debug >= notice_t || !result->IsValid()) {
328 cout << "Slice: "
329 << setw(4) << ix << ' '
330 << setw(16) << h1.GetName() << ' '
331 << FIXED(7,3) << f1->GetParameter(1) << " +/- "
332 << FIXED(7,3) << f1->GetParError(1) << ' '
333 << FIXED(9,2) << result->Chi2() << '/'
334 << result->Ndf() << ' '
335 << result->IsValid()
336 << E_ns(f1->GetParError(1)) << ' '
337 << (status ? "" : "failed") << endl;
338 }
339
340 counts += 1;
341
342 if (!status) {
343 errors += 1;
344 }
345
346 if (status) {
347 zmap[index] = result_type(f1->GetParameter(1), f1->GetParError(1));
348 }
349
350 if (result->Ndf() > 0) {
351 hc.SetBinContent(ix, result->Chi2() / result->Ndf());
352 }
353
354 hq.SetBinContent(ix, status ? 1.0 : 0.0);
355
356 h1.Write();
357
358 delete f1;
359 }
360
361
362 if (!zmap.empty()) {
363
364 // average time offset
365
366 double t0 = 0.0;
367
368 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
369 t0 += i->second.value;
370 }
371
372 t0 /= zmap.size();
373
374 NOTICE("Average time offset [ns] " << FIXED(7,2) << t0 << endl);
375 NOTICE("Number of fits passed/failed " << counts << "/" << errors << endl);
376
377 for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {
378 i->second.value -= t0;
379 }
380
381 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
382 h0.SetBinContent(i->first, i->second.value);
383 h0.SetBinError (i->first, i->second.error);
384 }
385
386 for (Int_t ix = 1; ix <= H2.GetXaxis()->GetNbins(); ++ix) {
387 for (Int_t iy = 1; iy <= H2.GetYaxis()->GetNbins(); ++iy) {
388
389 const JLocation location(strings.at(ix-1), iy);
390
391 if (router.hasLocation(location)) {
392
393 const JModule& module = router.getModule(location);
394 const int index = rpm->getIndex(module.getID());
395
396 if (zmap.count(index) != 0) {
397 H2.SetBinContent(ix, iy, zmap[index].value);
398 }
399 }
400 }
401 }
402
403 if (overwriteDetector) {
404
405 NOTICE("Store calibration data on file " << detectorFile << endl);
406
407 detector.comment.add(JMeta(argc, argv));
408
409 for (size_t string = 0; string != strings.size(); ++string) {
410 for (int floor = 1; floor <= floors.getUpperLimit(); ++floor) {
411
412 const JLocation location(strings.at(string), floor);
413
414 if (router.hasLocation(location)) {
415
416 JModule& module = detector[router.getIndex(location)];
417 const int index = rpm->getIndex(module.getID());
418
419 if (zmap.count(index) != 0) {
420 module.sub(zmap[index].value);
421 }
422 }
423 }
424 }
425
426 store(detectorFile, detector);
427 }
428
429 } else {
430
431 NOTICE("No calibration results." << endl);
432 }
433
434 h0 .Write();
435 hc .Write();
436 hq .Write();
437 h2->Write();
438 H2 .Write();
439
440 out.Close();
441
442 delete rpm;
443}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
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
WARNING
#define WARNING(A)
Definition JMessage.hh:65
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
JDETECTOR::JDetector
Detector data structure.
Definition JDetector.hh:96
JDETECTOR::JLocationRouter
Router for direct addressing of location data in detector data structure.
Definition JLocationRouter.hh:30
JDETECTOR::JLocation
Logical location of module.
Definition JLocation.hh:40
JDETECTOR::JModule
Data structure for a composite optical module.
Definition JModule.hh:76
JLANG::JException
General exception.
Definition JException.hh:24
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JCALIBRATE::module_t
const char *const module_t
routing by module
Definition JCalibrateMuon.hh:172
JCALIBRATE::string_t
const char *const string_t
routing by string
Definition JCalibrateMuon.hh:173
JDETECTOR::getRangeOfFloors
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Definition JDetectorToolkit.hh:295
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:520
JDETECTOR::store
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
Definition JDetectorToolkit.hh:617
JEEP::debug_t
@ debug_t
debug
Definition JMessage.hh:29
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JAANET::detector
Detector file.
Definition JHead.hh:227
JAANET::target
Target.
Definition JHead.hh:300
JCALIBRATE::JModuleRouter_t
Module router.
Definition JCalibrateMuon.hh:65
JCALIBRATE::JRouter_t
Interface for routing module identifier to index and vice versa.
Definition JCalibrateMuon.hh:25
JCALIBRATE::JStringRouter_t
String router.
Definition JCalibrateMuon.hh:120
JDETECTOR::JStringRouter
Router for mapping of string identifier to index.
Definition JStringRouter.hh:23
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488
Generated by doxygen 1.10.0