36{
39
42 string detectorFile;
43 bool overwriteDetector;
44 string option;
45 double factor;
47
48 try {
49
50 JParser<> zap(
"Program to fit histogram in output of JMuonCompass.cc.");
51
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;
59
60 zap(argc, argv);
61 }
62 catch(const exception &error) {
63 FATAL(error.what() << endl);
64 }
65
66
67 gErrorIgnoreLevel = kFatal;
68
69
71
72 try {
74 }
77 }
78
80
81 if (option.find('S') == string::npos) { option += 'S'; }
83
84 TH2D* h2 = NULL;
85
86 for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
87
88 NOTICE(
"Processing " << *i << endl);
89
90 TFile in(i->c_str(), "exist");
91
92 if (!in.IsOpen()) {
93 FATAL(
"File " << *i <<
" not opened." << endl);
94 }
95
96 TH2D* p = dynamic_cast<TH2D*>(in.Get(h2_t));
97
98 if (p == NULL) {
99 FATAL(
"File " << *i <<
" has no histogram <" << h2_t <<
">." << endl);
100 }
101
102 if (h2 == NULL)
103 h2 = (TH2D*) p->Clone();
104 else
105 h2->Add(p);
106
107 h2->SetDirectory(0);
108
109 in.Close();
110 }
111
112 if (h2 == NULL) {
113 FATAL(
"No histogram <" << h2_t <<
">." << endl);
114 }
115
116
117
118 struct result_type {
119
120 result_type() :
121 value(0.0),
122 error(0.0)
123 {}
124
125 result_type(double value,
126 double error) :
127 value(value),
128 error(error)
129 {}
130
131 double value;
132 double error;
133 };
134
136
137 map_type zmap;
138
139
140
141
142 const TAxis* x_axis = h2->GetXaxis();
143 const TAxis* y_axis = h2->GetYaxis();
144
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);
148
151
152 TH2D H2("detector", NULL,
153 strings.size(), -0.5, strings.size() - 0.5,
155
156 for (Int_t ix = 1; ix <= H2.GetXaxis()->GetNbins(); ++ix) {
157 H2.GetXaxis()->SetBinLabel(ix,
MAKE_CSTRING(strings.at(ix-1)));
158 }
159 for (Int_t iy = 1; iy <= H2.GetYaxis()->GetNbins(); ++iy) {
161 }
162
164
165 for (Int_t ix = 1; ix <= x_axis->GetNbins(); ++ix) {
166
167 const int index = ix - 1;
168 const int id =
detector[index].getID();
169
170 if (
detector[index].getFloor() == 0) {
171 continue;
172 }
173
174 DEBUG(
"Bin " << setw(4) << ix <<
" -> " << setw(8) <<
id << endl);
175
176 TH1D h1(
MAKE_CSTRING(
id <<
".1D"), NULL, y_axis->GetNbins(), y_axis->GetXmin(), y_axis->GetXmax());
177
178
179
180 Double_t ymin = numeric_limits<double>::max();
181 Double_t ymax = numeric_limits<double>::lowest();
182 Double_t x0 = 0.0;
183
184 for (Int_t iy = 1; iy <= y_axis->GetNbins(); ++iy) {
185
186 h1.SetBinContent(iy, h2->GetBinContent(ix,iy));
187 h1.SetBinError (iy, h2->GetEntries() * 1.0e-7 * factor);
188
189 const Double_t
x = h1.GetBinCenter (iy);
190 const Double_t
y = h1.GetBinContent(iy);
191
192 if (y < ymin) {
194 }
195
196 if (y > ymax) { ymax =
y; }
197 if (y < ymin) { ymin =
y; }
198 }
199
200 Double_t
xmin = x0 - 0.4*PI;
201 Double_t
xmax = x0 + 0.4*PI;
202
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; }
205
206 TF1*
f1 =
new TF1(
"f1",
"[0] - [1]*cos(x - [2])");
207
208 f1->SetParameter(0, 0.5 * (ymin + ymax));
209 f1->SetParameter(1, 0.5 * (ymax - ymin));
210 f1->SetParameter(2, x0);
211
212 f1->SetParError(0, 0.01*(ymin+ymax));
213 f1->SetParError(1, 0.01*(ymin+ymax));
214 f1->SetParError(2, 0.05);
215
216 DEBUG(
"Start values:" << endl);
217
218 for (
int i = 0; i !=
f1->GetNpar(); ++i) {
219 DEBUG(left << setw(12) <<
f1->GetParName (i) <<
' ' <<
221 }
222
223
224
225 TFitResultPtr
result = h1.Fit(f1, option.c_str(),
"same", xmin, xmax);
226
227 const bool status =
result->IsValid();
228
230 cout << "Slice: "
231 << setw(4) << ix << ' '
232 << setw(16) << h1.GetName() << ' '
233 <<
FIXED(7,3) <<
f1->GetParameter(2) <<
" +/- "
234 <<
FIXED(7,3) <<
f1->GetParError(2) <<
' '
238 << (status ? "" : "failed") << endl;
239 }
240
241 if (status) {
242 zmap[index] = result_type(
f1->GetParameter(2),
f1->GetParError(2));
243 }
244
247 }
248
249 hq.SetBinContent(ix, status ? 1.0 : 0.0);
250
251 h1.Write();
252
254 }
255
256 if (zmap.empty()) {
257 NOTICE(
"No calibration results." << endl);
258 }
259
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);
263 }
264
265 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
266
268
269 const int ix = i->first + 1;
270 const int iy = module.getFloor();
271
272 H2.SetBinContent(ix, iy, i->second.value);
273 }
274
275 if (overwriteDetector) {
276
277 NOTICE(
"Store calibration data on file " << detectorFile << endl);
278
280
281 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
282
284 const JVector3D center =
module.getPosition();
285
286 module.sub(center);
287
288 module.rotate(JRotation3Z(i->second.value));
289
290 module.add(center);
291 }
292
294 }
295
296 h0 .Write();
297 hc .Write();
298 hq .Write();
299 h2->Write();
300 H2 .Write();
301
302 out.Close();
303}
#define DEBUG(A)
Message macros.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
#define MAKE_CSTRING(A)
Make C-string.
Router for direct addressing of module data in detector data structure.
Data structure for a composite optical module.
Data structure for vector in three dimensions.
Utility class to parse command line options.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
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).
Auxiliary data structure for floating point format specification.
Router for mapping of string identifier to index.
Auxiliary data structure for floating point format specification.