13 #include "TFitResult.h"
35 const char* h2_t =
"h2";
37 const char* gauss_t =
"Gauss";
38 const char* landau_t =
"Landau";
39 const char* emg_t =
"EMG";
40 const char* breitwigner_t =
"BreitWigner";
49 int main(
int argc,
char **argv)
57 bool overwriteDetector;
68 JParser<> zap(
"Program to fit time-residuals histogram in output of JCalibrateMuon.cc.");
70 zap[
'f'] =
make_field(inputFile,
"input files (output from JCalibrateMuon).");
72 zap[
'a'] =
make_field(detectorFile,
"detector file.");
73 zap[
'A'] =
make_field(overwriteDetector,
"overwrite detector file provided through '-a' with correct time offsets.");
74 zap[
'F'] =
make_field(
function,
"fit function") = gauss_t, landau_t, emg_t, breitwigner_t;
76 zap[
'O'] =
make_field(option,
"ROOT fit option, see TH1::Fit.") =
"";
78 zap[
'W'] =
make_field(WMin,
"minimal histogram contents.") = 100.0;
84 catch(
const exception &error) {
85 FATAL(error.what() << endl);
90 FATAL(
"Invalid fit range [ns] " << T_ns << endl);
103 if (option.find(
'S') == string::npos) { option +=
'S'; }
110 NOTICE(
"Processing " << *i << endl);
112 TFile in(i->c_str(),
"exist");
115 FATAL(
"File " << *i <<
" not opened." << endl);
118 TH2D* p =
dynamic_cast<TH2D*
>(in.Get(h2_t));
121 FATAL(
"File " << *i <<
" has no histogram <" << h2_t <<
">." << endl);
125 h2 = (TH2D*) p->Clone();
135 FATAL(
"No histogram <" << h2_t <<
">." << endl);
147 result_type(
double value,
164 const TAxis* x_axis = h2->GetXaxis();
165 const TAxis* y_axis = h2->GetYaxis();
167 TH1D h0(
"h0", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
168 TH1D hc(
"hc", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
169 TH1D hq(
"hq", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
177 for (Int_t ix = 1; ix <= x_axis->GetNbins(); ++ix) {
191 TH1D h1(
MAKE_CSTRING(module.
getID() <<
".1D"), NULL, y_axis->GetNbins(), y_axis->GetXmin(), y_axis->GetXmax());
197 Double_t
sigma = 4.0;
200 for (Int_t iy = 1; iy <= y_axis->GetNbins(); ++iy) {
202 h1.SetBinContent(iy, h2->GetBinContent(ix,iy));
203 h1.SetBinError (iy, sqrt(h2->GetBinContent(ix,iy)));
205 const Double_t
x = h1.GetBinCenter (iy);
206 const Double_t
y = h1.GetBinContent(iy);
218 WARNING(
"Not enough entries for slice " << ix <<
' ' << W <<
"; skip fit." << endl);
231 if (
function == gauss_t) {
233 f1 =
new TF1(
function.c_str(),
"[0]*TMath::Gaus(x,[1],[2])");
235 f1->SetParameter(0, 0.8*ymax);
236 f1->SetParameter(1, x0);
239 f1->SetParError(0, sqrt(ymax) * 0.1);
240 f1->SetParError(1, 0.01);
241 f1->SetParError(2, 0.01);
243 }
else if (
function == landau_t) {
245 f1 =
new TF1(
function.c_str(),
"[0]*TMath::Landau(x,[1],[2])");
247 f1->SetParameter(0, 0.8*ymax);
248 f1->SetParameter(1, x0);
251 f1->SetParError(0, sqrt(ymax) * 0.1);
252 f1->SetParError(1, 0.01);
253 f1->SetParError(2, 0.01);
255 }
else if (
function == emg_t) {
257 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]))");
259 f1->SetParameter(0, 0.5*ymax);
260 f1->SetParameter(1, x0 -
sigma);
262 f1->SetParameter(3, 0.06);
264 f1->SetParError(0, sqrt(ymax) * 0.1);
265 f1->SetParError(1, 0.01);
266 f1->SetParError(2, 0.01);
267 f1->SetParError(3, 1.0e-4);
269 }
else if (
function == breitwigner_t) {
271 f1 =
new TF1(
function.c_str(),
"(x <= [1])*[0]*[2]*TMath::BreitWigner(x,[1],[2])+(x > [1])*[0]*[3]*TMath::BreitWigner(x,[1],[3])");
273 f1->SetParameter(0, 0.8*ymax);
274 f1->SetParameter(1, x0);
275 f1->SetParameter(2, 15.0);
276 f1->SetParameter(3, 25.0);
278 f1->SetParError(0, sqrt(ymax) * 0.1);
279 f1->SetParError(1, 0.01);
280 f1->SetParError(2, 0.1);
281 f1->SetParError(3, 0.1);
285 FATAL(
"Unknown fit function " <<
function << endl);
289 DEBUG(
"Start values:" << endl);
291 for (
int i = 0; i !=
f1->GetNpar(); ++i) {
292 DEBUG(left << setw(12) <<
f1->GetParName (i) <<
' ' <<
302 << setw(4) << ix <<
' '
303 << setw(16) << h1.GetName() <<
' '
304 <<
FIXED(7,3) <<
f1->GetParameter(1) <<
" +/- "
305 <<
FIXED(7,3) <<
f1->GetParError(1) <<
' '
308 << (
result->IsValid() ?
"" :
"failed") << endl;
318 zmap[ix] = result_type(
f1->GetParameter(1),
f1->GetParError (1));
325 hq.SetBinContent(ix,
result->IsValid() ? 1.0 : 0.0);
339 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
340 t0 += i->second.value;
345 NOTICE(
"Average time offset [ns] " <<
FIXED(7,2) << t0 << endl);
346 NOTICE(
"Number of fits passed/failed " << counts <<
"/" << errors << endl);
348 for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {
349 i->second.value -= t0;
352 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
353 h0.SetBinContent(i->first, i->second.value);
354 h0.SetBinError (i->first, i->second.error);
362 string.size(), -0.5,
string.size() - 0.5,
365 for (Int_t i = 1; i <= hi.GetXaxis()->GetNbins(); ++i) {
366 hi.GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(
string.at(i-1)));
368 for (Int_t i = 1; i <= hi.GetYaxis()->GetNbins(); ++i) {
372 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
373 hi.SetBinContent(
detector[i->first - 1].getString(),
381 if (overwriteDetector) {
383 NOTICE(
"Store calibration data on file " << detectorFile << endl);
387 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
389 if (E_ns(i->second.error))
390 detector[i->first - 1].sub(i->second.value);
392 ERROR(
"Slice " << setw(4) << i->first <<
" fit uncertainty " <<
FIXED(5,2) << i->second.error <<
" outside specified range (option -E <E_ns>)" << endl);
400 NOTICE(
"No calibration results." << endl);
Data structure for detector geometry and calibration.
int main(int argc, char **argv)
General purpose messaging.
#define DEBUG(A)
Message macros.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
Auxiliary class to define a range between two values.
Direct access to string in detector data structure.
int getFloor() const
Get floor number.
const JLocation & getLocation() const
Get location.
Data structure for a composite optical module.
int getID() const
Get identifier.
Utility class to parse command line options.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [s]).
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
std::map< int, range_type > map_type
Auxiliary data structure for floating point format specification.
Router for mapping of string identifier to index.
Auxiliary data structure for floating point format specification.