35int main(
int argc,
char **argv)
43 bool overwriteDetector;
50 JParser<> zap(
"Program to fit histogram in output of JMuonCompass.cc.");
52 zap[
'f'] =
make_field(inputFile,
"input files (output from JMuonCompass).");
54 zap[
'a'] =
make_field(detectorFile,
"detector file.");
55 zap[
'A'] =
make_field(overwriteDetector,
"overwrite detector file provided through '-a' with correct time offsets.");
56 zap[
'O'] =
make_field(option,
"ROOT fit option, see TH1::Fit.") =
"";
57 zap[
'E'] =
make_field(factor,
"scaling factor for error on bin contents") = 1.0;
62 catch(
const exception &error) {
63 FATAL(error.what() << endl);
67 gErrorIgnoreLevel = kFatal;
81 if (option.find(
'S') == string::npos) { option +=
'S'; }
86 for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
88 NOTICE(
"Processing " << *i << endl);
90 TFile in(i->c_str(),
"exist");
93 FATAL(
"File " << *i <<
" not opened." << endl);
96 TH2D* p =
dynamic_cast<TH2D*
>(in.Get(h2_t));
99 FATAL(
"File " << *i <<
" has no histogram <" << h2_t <<
">." << endl);
103 h2 = (TH2D*) p->Clone();
113 FATAL(
"No histogram <" << h2_t <<
">." << endl);
125 result_type(
double value,
142 const TAxis* x_axis = h2->GetXaxis();
143 const TAxis* y_axis = h2->GetYaxis();
145 TH1D h0(
"h0", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
146 TH1D hc(
"hc", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
147 TH1D hq(
"hq", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
152 TH2D H2(
"detector", NULL,
153 strings.size(), -0.5, strings.size() - 0.5,
156 for (Int_t ix = 1; ix <= H2.GetXaxis()->GetNbins(); ++ix) {
157 H2.GetXaxis()->SetBinLabel(ix,
MAKE_CSTRING(strings.at(ix-1)));
159 for (Int_t iy = 1; iy <= H2.GetYaxis()->GetNbins(); ++iy) {
165 for (Int_t ix = 1; ix <= x_axis->GetNbins(); ++ix) {
167 const int index = ix - 1;
168 const int id =
detector[index].getID();
170 if (
detector[index].getFloor() == 0) {
174 DEBUG(
"Bin " << setw(4) << ix <<
" -> " << setw(8) <<
id << endl);
176 TH1D h1(
MAKE_CSTRING(
id <<
".1D"), NULL, y_axis->GetNbins(), y_axis->GetXmin(), y_axis->GetXmax());
180 Double_t ymin = numeric_limits<double>::max();
181 Double_t ymax = numeric_limits<double>::lowest();
184 for (Int_t iy = 1; iy <= y_axis->GetNbins(); ++iy) {
186 h1.SetBinContent(iy, h2->GetBinContent(ix,iy));
187 h1.SetBinError (iy, h2->GetEntries() * 1.0e-7 * factor);
189 const Double_t x = h1.GetBinCenter (iy);
190 const Double_t y = h1.GetBinContent(iy);
196 if (y > ymax) { ymax = y; }
197 if (y < ymin) { ymin = y; }
200 Double_t xmin = x0 - 0.4*PI;
201 Double_t xmax = x0 + 0.4*PI;
203 if (xmin < y_axis->GetXmin()) { xmin = y_axis->GetXmin(); xmax = xmin + 0.8*PI; }
204 if (xmax > y_axis->GetXmax()) { xmax = y_axis->GetXmax(); xmin = xmax - 0.8*PI; }
206 TF1* f1 =
new TF1(
"f1",
"[0] - [1]*cos(x - [2])");
208 f1->SetParameter(0, 0.5 * (ymin + ymax));
209 f1->SetParameter(1, 0.5 * (ymax - ymin));
210 f1->SetParameter(2, x0);
212 f1->SetParError(0, 0.01*(ymin+ymax));
213 f1->SetParError(1, 0.01*(ymin+ymax));
214 f1->SetParError(2, 0.05);
216 DEBUG(
"Start values:" << endl);
218 for (
int i = 0; i != f1->GetNpar(); ++i) {
219 DEBUG(left << setw(12) << f1->GetParName (i) <<
' ' <<
220 SCIENTIFIC(12,5) << f1->GetParameter(i) << endl);
225 TFitResultPtr result = h1.Fit(f1, option.c_str(),
"same", xmin, xmax);
227 const bool status = result->IsValid();
229 if (
debug >= notice_t || !result->IsValid()) {
231 << setw(4) << ix <<
' '
232 << setw(16) << h1.GetName() <<
' '
233 <<
FIXED(7,3) << f1->GetParameter(2) <<
" +/- "
234 <<
FIXED(7,3) << f1->GetParError(2) <<
' '
235 <<
FIXED(9,2) << result->Chi2() <<
'/'
236 << result->Ndf() <<
' '
238 << (status ?
"" :
"failed") << endl;
242 zmap[index] = result_type(f1->GetParameter(2), f1->GetParError(2));
245 if (result->Ndf() > 0) {
246 hc.SetBinContent(ix, result->Chi2() / result->Ndf());
249 hq.SetBinContent(ix, status ? 1.0 : 0.0);
257 NOTICE(
"No calibration results." << endl);
260 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
261 h0.SetBinContent(i->first, i->second.value);
262 h0.SetBinError (i->first, i->second.error);
265 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
269 const int ix = i->first + 1;
270 const int iy =
module.getFloor();
272 H2.SetBinContent(ix, iy, i->second.value);
275 if (overwriteDetector) {
277 NOTICE(
"Store calibration data on file " << detectorFile << endl);
281 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
284 const JVector3D center =
module.getPosition();
288 module.rotate(JRotation3Z(i->second.value));