64{
68
70
72
73 string inputFile;
75 string detectorFile;
76 string pmtFile;
78 bool reverse;
79 bool overwriteDetector;
81 bool fitAngle;
82 bool fitNoise;
83 bool fitModel;
84 bool fixQE;
88
89 try {
90
92
105
106 JParser<> zap(
"Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
107
109 zap[
'f'] =
make_field(inputFile,
"input file (output from JMergeCalibrateK40).");
111 zap[
'a'] =
make_field(detectorFile,
"detector file.");
112 zap[
'P'] =
make_field(pmtFile,
"specify PMT file name that can be given to JTriggerEfficiency.") =
"";
114 "Fix time offset(s) of PMT(s) of certain module(s), e.g."
115 "\n-! \"808969848 0 808982077 23\" will fix time offsets of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
116 "\nSame PMT offset can be fixed for all optical modules, e.g."
117 "\n-! \"-1 0 -1 23\" will fix time offsets of PMTs 0 and 23 of all optical modules.") =
JPARSER::initialised();
118 zap[
'r'] =
make_field(reverse,
"reverse TDC constraints due to option -! <TDC>.");
119 zap[
'A'] =
make_field(overwriteDetector,
"overwrite detector file provided through '-a' with fitted time offsets.");
120 zap[
'w'] =
make_field(writeFits,
"write fit results to ROOT file; -ww also write fitted TTS to PMT file.");
121 zap[
'D'] =
make_field(fitAngle,
"fit angular distribution; fix normalisation.");
122 zap[
'B'] =
make_field(fitNoise,
"fit background.");
123 zap[
'M'] =
make_field(fitModel,
"fit angular distribution as well as normalisation; fix PMT QEs = 1.0.");
124 zap[
'Q'] =
make_field(fixQE,
"fix PMT QEs = 1.0.");
128
129 zap(argc, argv);
130 }
131 catch(const exception &error) {
132 FATAL(error.what() << endl);
133 }
134
135
136 if ((fitModel ? 1 : 0) +
137 (fitAngle ? 1 : 0) +
138 (fitNoise ? 1 : 0) +
139 (fixQE ? 1 : 0) > 1) {
140 FATAL(
"Use either option -M, -D, -B or -Q" << endl);
141 }
142
148
149 if (reverse) {
151 }
152
153 for (JTDC_t::const_iterator i = TDC.begin(); i != TDC.end(); ++i) {
154 DEBUG(
"PMT " << setw(10) << i->first <<
' ' << setw(2) << i->second <<
" constrain t0." << endl);
155 }
156
157 try {
159 }
160 catch(const exception &error) {
161 FATAL(error.what() << endl);
162 }
163
165
166 try {
168 }
171 }
172
173
175
176 if (pmtFile != "") {
177 try {
178 parameters.
load(pmtFile.c_str());
179 }
180 catch(const exception& error) {}
181 }
182
184
185
186 TFile* in = TFile::Open(inputFile.c_str(), "exist");
187
188 if (in == NULL || !in->IsOpen()) {
189 FATAL(
"File: " << inputFile <<
" not opened." << endl);
190 }
191
192
194
195 TH1D h0("chi2", NULL, 500, 0.0, 5.0);
196 TH1D hn("hn", NULL, 501, -0.5, 500.0);
197 TH1D hr("rate", NULL, 500, 0.0, 25.0);
198 TH1D h1("p1", NULL, 500, -5.0, +5.0);
199 TH1D h2("p2", NULL, 500, -5.0, +5.0);
200 TH1D h3("p3", NULL, 500, -5.0, +5.0);
201 TH1D h4("p4", NULL, 500, -5.0, +5.0);
202 TH1D hc("cc", NULL, 500, -0.1, +0.1);
203 TH1D hb("bc", NULL, 500, -0.1, +0.1);
204
207
208 TH2D H2("detector", NULL,
209 string.size() + 0, -0.5, string.size() - 0.5,
211
212 for (Int_t i = 1; i <= H2.GetXaxis()->GetNbins(); ++i) {
213 H2.GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(
string.at(i-1)));
214 }
215 for (Int_t i = 1; i <= H2.GetYaxis()->GetNbins(); ++i) {
217 }
218
219 TH2D* HN = (TH2D*) H2.Clone("iterations");
220
222
223 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
224
225 if (module->getFloor() == 0) {
226 continue;
227 }
228
230
231 NOTICE(
"Module " << setw(10) << module->getID() <<
' ' <<
getLabel(module->getLocation()) <<
" !" <<
distance(range.first, range.second) << endl);
232
234
235 if (h2d == NULL || h2d->GetEntries() == 0) {
236
237 NOTICE(
"No data for module " << module->getID() <<
" -> set QEs to 0." << endl);
238
241 }
242
243 continue;
244 }
245
247
249
252
253 for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
254
256
258
259 double V = 0.0;
260 double W = 0.0;
261
262 for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
263
264 const double x = h2d->GetXaxis()->GetBinCenter(ix);
265 const double y = h2d->GetYaxis()->GetBinCenter(iy);
266
267 if (X(x) && Y(y)) {
268
269 double value = h2d->GetBinContent(ix,iy);
270 double error = h2d->GetBinError (ix,iy);
271
272 buffer.push_back(
rate_type(y, value, error));
273
274 double width = h2d->GetYaxis()->GetBinWidth(iy);
275
276 value *= width;
277 error *= width;
278
279 V += value;
280 W += error * error;
281 }
282 }
283
284 W = sqrt(W);
285
286 if (V <= 0.0 - STDEV*W) {
287 count[0][
pair.first] += 1;
288 count[0][
pair.second] += 1;
289 }
290
291 if (V <= MINIMAL_RATE_HZ + STDEV*W) {
292 count[1][
pair.first] += 1;
293 count[1][
pair.second] += 1;
294 }
295 }
296
298
299 if (count[0][pmt] >= MAXIMAL_COUNTS) {
300
301 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
" some rates negative -> fit background" << endl);
302
303 if (fit.value.parameters[pmt].status) {
304 model.parameters[pmt].bg.set();
305 }
306 }
307
308 if (count[1][pmt] == NUMBER_OF_PMTS) {
309
310 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
" all rates to low -> disable" << endl);
311
312 model.parameters[pmt].disable();
313 }
314 }
315
316 DEBUG(
"Start value:" << endl <<
model << endl);
317
318 try {
319
321
323
325
326 ERROR(
"Fit result " << setw(10) << module->getID() <<
" NDF " << setw(5) <<
result.ndf <<
" -> skip" << endl);
327
328 continue;
329 }
330
331 bool refit = false;
332
334
335 if (fit.value.parameters[pmt].status) {
336
337 if (fit.value.parameters[pmt].QE() <= QE_MIN ||
338 fit.value.parameters[pmt].QE() <= 0.0 + STDEV * fit.error.parameters[pmt].QE()) {
339
340 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
' '
341 << "QE = "
342 <<
FIXED(5,3) << fit.value.parameters[pmt].QE() <<
" +/- "
343 <<
FIXED(5,3) << fit.error.parameters[pmt].QE() <<
" "
344 << " -> disable" << (!refit ? " and refit" : "") << endl);
345
346 fit.value.parameters[pmt].disable();
347
348 refit = true;
349 }
350
351 if (fit.value.parameters[pmt].t0.atLimit(T0_NS)) {
352
353 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
' '
354 << "t0 at limit "
355 <<
FIXED(5,3) << fit.value.parameters[pmt].t0() <<
" +/- "
356 <<
FIXED(5,3) << fit.error.parameters[pmt].t0() <<
" "
357 << " -> disable" << (!refit ? " and refit" : "") << endl);
358
359 fit.value.parameters[pmt].disable();
360
361 refit = true;
362 }
363 }
364 }
365
366 if (refit) {
367
369 if (fit.value.parameters[pmt].status) {
371 }
372 }
373
374 refit = false;
376 }
377
378 NOTICE(
"Fit result " << setw(10) << module->getID() <<
" chi2 / NDF " <<
FIXED(10,2) <<
result.chi2 <<
" / " << setw(5) <<
result.ndf <<
' ' << setw(5) << fit.numberOfIterations << endl);
379
382 }
383
385
386 if (writeFits) {
387
389 hn.Fill(fit.numberOfIterations);
390 hr.Fill(fit.value.R );
391 h1.Fill(fit.value.p1);
392 h2.Fill(fit.value.p2);
393 h3.Fill(fit.value.p3);
394 h4.Fill(fit.value.p4);
395 hc.Fill(fit.value.cc);
396 hb.Fill(fit.value.bc);
397
398 TH1D h1t(
MAKE_CSTRING(module->getID() <<
".1t0"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
399 TH1D h1s(
MAKE_CSTRING(module->getID() <<
".1TTS"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
400 TH1D h1q(
MAKE_CSTRING(module->getID() <<
".1QE"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
401
403 h1t.SetBinContent(pmt + 1, fit.value.parameters[pmt].t0 ());
404 h1t.SetBinError (pmt + 1, fit.error.parameters[pmt].t0 () + numeric_limits<double>::epsilon());
405 h1s.SetBinContent(pmt + 1, fit.value.parameters[pmt].TTS());
406 h1s.SetBinError (pmt + 1, fit.error.parameters[pmt].TTS() + numeric_limits<double>::epsilon());
407 h1q.SetBinContent(pmt + 1, fit.value.parameters[pmt].QE ());
408 h1q.SetBinError (pmt + 1, fit.error.parameters[pmt].QE () + numeric_limits<double>::epsilon());
409 }
410
411 out << h1t << h1s << h1q;
412
413 for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
414
416
417 for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
418
419 const double dt_ns = h2d->GetYaxis()->GetBinCenter(iy);
420
421 h2d->SetBinContent(ix, iy, fit.value.getValue(
pair, dt_ns));
422 h2d->SetBinError (ix, iy, 0.0);
423 }
424 }
425
427 h2d->Write();
428
429 const double x =
string.getIndex(module->getString());
430 const double y =
module->getFloor();
431
433 HN->Fill(x, y, fit.numberOfIterations);
434 }
435
436 const double t0 = (fit.value.hasFixedTimeOffset() ? fit.value.getFixedTimeOffset() : 0.0);
437
439
441
443
444 if (P > 0.0)
445 data.QE = fit.value.parameters[pmt].QE / P;
446 else
448
449 if (writeFits > 1) {
450 data.TTS_ns = fit.value.parameters[pmt].TTS();
451 }
452
453 module->getPMT(pmt).addT0(fit.value.parameters[pmt].t0.get() - t0);
454 }
455 }
456 catch(const exception& error) {
457
458 ERROR(
"Module " << setw(10) << module->getID() <<
' ' << error.what() <<
" -> set QEs to 0." << endl);
459
462 }
463 }
464 }
465
466
468
469 {
471 JSTDObjectWriter <JMeta> writer(meta);
472
473 writer << reader;
474 }
475
476 for (vector<JMeta>::const_reverse_iterator i = meta.rbegin(); i != meta.rend(); ++i) {
479 }
480
481 if (overwriteDetector) {
482
483 NOTICE(
"Store calibration data on file " << detectorFile << endl);
484
486 }
487
488 if (pmtFile != "") {
489 parameters.
store(pmtFile.c_str());
490 }
491
492
493 for (vector<JMeta>::const_iterator i = meta.begin(); i != meta.end(); ++i) {
495 }
496
499 }
500
501 if (writeFits) {
502 out << h0 << hn << hr << h1 << h2 << h3 << h4 << hc << hb << H2 << *HN;
503 }
504
505 out.Close();
506
507 return 0;
508}
#define DEBUG(A)
Message macros.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
#define MAKE_CSTRING(A)
Make C-string.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Auxiliary class for map of PMT parameters.
Data structure for PMT parameters.
Utility class to parse parameter values.
Utility class to parse command line options.
Object reading from a list of files.
static double TEROSTAT_R1
scaling factor
static const char *const _2F
Name extension for 2F rate fitted.
@ FIT_PMTS_QE_FIXED_t
fit parameters of PMTs with QE fixed
@ FIT_PMTS_AND_ANGULAR_DEPENDENCE_t
fit parameters of PMTs and angular dependence of K40 rate
@ FIT_MODEL_t
fit parameters of K40 rate and TTSs of PMTs
@ FIT_PMTS_AND_BACKGROUND_t
fit parameters of PMTs and background
@ FIT_PMTS_t
fit parameters of PMTs
static const char *const _2R
Name extension for 2D rate measured.
static double TEROSTAT_DZ
maximal PMT inclination
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.
double getSurvivalProbability(const JPMTParameters ¶meters)
Get model dependent probability that a one photo-electron hit survives the simulation of the PMT assu...
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
KM3NeT DAQ data structures and auxiliaries.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Auxiliary data structure for sequence of same character.
Auxiliary data structure for floating point format specification.
Type definition of range.
Model for fit to acoustics data.
Fit parameters for two-fold coincidence rate due to K40.
static const JK40Parameters & getInstance()
Get default values.
static const JPMTParameters_t & getInstance()
Get default values.
Auxiliary class for TDC constraints.
range_type equal_range(const int id)
Get range of constraints for given module.
void reverse()
Reverse constraints.
bool is_valid(const bool option=false) const
Check validity of TDC constrants.
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Data structure for measured coincidence rate of pair of PMTs.
Router for mapping of string identifier to index.
void store(const char *file_name) const
Store to output file.
void load(const char *file_name)
Load from input file.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...