65{
69
71
73
74 string inputFile;
76 string detectorFile;
77 string pmtFile;
79 bool reverse;
80 bool overwriteDetector;
82 bool fitAngle;
83 bool fitNoise;
84 bool fitModel;
85 bool fixQE;
90
91 try {
92
94
109
110 JParser<> zap(
"Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
111
113 zap[
'f'] =
make_field(inputFile,
"input file (output from JMergeCalibrateK40).");
115 zap[
'a'] =
make_field(detectorFile,
"detector file.");
116 zap[
'P'] =
make_field(pmtFile,
"specify PMT file name that can be given to JTriggerEfficiency.") =
"";
118 "Fix time offset(s) of PMT(s) of certain module(s), e.g."
119 "\n-! \"808969848 0 808982077 23\" will fix time offsets of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
120 "\nSame PMT offset can be fixed for all optical modules, e.g."
121 "\n-! \"-1 0 -1 23\" will fix time offsets of PMTs 0 and 23 of all optical modules.") =
JPARSER::initialised();
122 zap[
'r'] =
make_field(reverse,
"reverse TDC constraints due to option -! <TDC>.");
123 zap[
'A'] =
make_field(overwriteDetector,
"overwrite detector file provided through '-a' with fitted time offsets.");
124 zap[
'w'] =
make_field(writeFits,
"write fit results to ROOT file; -ww also write fitted TTS to PMT file.");
125 zap[
'D'] =
make_field(fitAngle,
"fit angular distribution; fix normalisation.");
126 zap[
'B'] =
make_field(fitNoise,
"fit background.");
127 zap[
'M'] =
make_field(fitModel,
"fit angular distribution as well as normalisation; fix PMT QEs = 1.0.");
128 zap[
'Q'] =
make_field(fixQE,
"fix PMT QEs = 1.0.");
133
134 zap(argc, argv);
135 }
136 catch(const exception &error) {
137 FATAL(error.what() << endl);
138 }
139
140
141 if ((fitModel ? 1 : 0) +
142 (fitAngle ? 1 : 0) +
143 (fitNoise ? 1 : 0) +
144 (fixQE ? 1 : 0) > 1) {
145 FATAL(
"Use either option -M, -D, -B or -Q" << endl);
146 }
147
153
154 if (reverse) {
156 }
157
158 for (JTDC_t::const_iterator i = TDC.begin(); i != TDC.end(); ++i) {
159 DEBUG(
"PMT " << setw(10) << i->first <<
' ' << setw(2) << i->second <<
" constrain t0." << endl);
160 }
161
162 try {
164 }
165 catch(const exception &error) {
166 FATAL(error.what() << endl);
167 }
168
170
171 try {
173 }
176 }
177
178
180
181 if (pmtFile != "") {
182 try {
183 parameters.
load(pmtFile.c_str());
184 }
185 catch(const exception& error) {}
186 }
187
189
190
191 TFile* in = TFile::Open(inputFile.c_str(), "exist");
192
193 if (in == NULL || !in->IsOpen()) {
194 FATAL(
"File: " << inputFile <<
" not opened." << endl);
195 }
196
197
199
200 TH1D h0("chi2", NULL, 500, 0.0, 5.0);
201 TH1D hn("hn", NULL, 501, -0.5, 500.0);
202 TH1D hr("rate", NULL, 500, 0.0, 25.0);
203 TH1D h1("p1", NULL, 500, -5.0, +5.0);
204 TH1D h2("p2", NULL, 500, -5.0, +5.0);
205 TH1D h3("p3", NULL, 500, -5.0, +5.0);
206 TH1D h4("p4", NULL, 500, -5.0, +5.0);
207 TH1D hc("cc", NULL, 500, -0.1, +0.1);
208 TH1D hb("bc", NULL, 500, -0.1, +0.1);
209
212
213 TH2D H2("detector", NULL,
214 string.size() + 0, -0.5, string.size() - 0.5,
216
217 for (Int_t i = 1; i <= H2.GetXaxis()->GetNbins(); ++i) {
218 H2.GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(
string.at(i-1)));
219 }
220 for (Int_t i = 1; i <= H2.GetYaxis()->GetNbins(); ++i) {
222 }
223
224 TH2D* HN = (TH2D*) H2.Clone("iterations");
225
227
228 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
229
230 if (module->getFloor() == 0) {
231 continue;
232 }
233
235
236 NOTICE(
"Module " << setw(10) << module->getID() <<
' ' <<
getLabel(module->getLocation()) <<
" !" <<
distance(range.first, range.second) << endl);
237
239
240 if (h2d == NULL || h2d->GetEntries() == 0) {
241
242 NOTICE(
"No data for module " << module->getID() <<
" -> set QEs to 0." << endl);
243
246 }
247
248 continue;
249 }
250
252
254
258 };
259
260 for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
261
263
265
266 double V = 0.0;
267 double W = 0.0;
268
269 for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
270
271 const double x = h2d->GetXaxis()->GetBinCenter(ix);
272 const double y = h2d->GetYaxis()->GetBinCenter(iy);
273
275
276 double value = h2d->GetBinContent(ix,iy);
277 double error = h2d->GetBinError (ix,iy);
278
279 buffer.push_back(
rate_type(y, value, error));
280
281 double width = h2d->GetYaxis()->GetBinWidth(iy);
282
283 value *= width;
284 error *= width;
285
286 V += value;
287 W += error * error;
288 }
289 }
290
291 W = sqrt(W);
292
293 if (V <= 0.0 - STDEV*W) {
294 count[0][
pair.first] += 1;
295 count[0][
pair.second] += 1;
296 }
297
298 if (V <= MINIMAL_RATE_HZ + STDEV*W) {
299 count[1][
pair.first] += 1;
300 count[1][
pair.second] += 1;
301 }
302 }
303
305
306 if (count[0][pmt] >= MAXIMAL_COUNTS) {
307
308 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
" some rates negative -> fit background" << endl);
309
310 if (fit.value.parameters[pmt].status) {
311 model.parameters[pmt].bg.set();
312 }
313 }
314
315 if (count[1][pmt] == NUMBER_OF_PMTS) {
316
317 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
" all rates to low -> disable" << endl);
318
319 model.parameters[pmt].disable();
320 }
321 }
322
323 DEBUG(
"Start value:" << endl <<
model << endl);
324
325 try {
326
328
330
332
334
335 ERROR(
"Fit result " << setw(10) << module->getID() <<
" NDF " << setw(5) <<
result.ndf <<
" -> skip" << endl);
336
337 continue;
338 }
339
340 bool refit = false;
341
343
344 if (fit.value.parameters[pmt].status) {
345
346 if (fit.value.parameters[pmt].QE() <= QE_MIN ) {
347
348 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
' '
349 << "QE = "
350 <<
FIXED(5,3) << fit.value.parameters[pmt].QE() <<
" +/- "
351 <<
FIXED(5,3) << fit.error.parameters[pmt].QE() <<
" "
352 << " -> disable" << (!refit ? " and refit" : "") << endl);
353
354 fit.value.parameters[pmt].disable();
355
356 refit = true;
357 }
358 }
359 }
360
362
363 if (fit.value.parameters[pmt].status) {
364
365 if (fit.value.parameters[pmt].t0.atLimit(T0_NS)) {
366
367 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
' '
368 << "t0 at limit "
369 <<
FIXED(7,3) << fit.value.parameters[pmt].t0() <<
" +/- "
370 <<
FIXED(7,3) << fit.error.parameters[pmt].t0());
371
372 if (refit == false) {
373
374 WARNING(
" -> disable and refit" << endl);
375
376 fit.value.parameters[pmt].disable();
377
378 refit = true;
379
380 } else {
381
383
384 fit.value.parameters[pmt].t0.set(0.0);
385 }
386 }
387 }
388 }
389
390 if (refit) {
391
393 if (fit.value.parameters[pmt].status) {
395 }
396 }
397
398 refit = false;
400 }
401
402 NOTICE(
"Fit result " << setw(10) << module->getID() <<
" chi2 / NDF " <<
FIXED(10,2) <<
result.chi2 <<
" / " << setw(5) <<
result.ndf <<
' ' << setw(5) << fit.numberOfIterations << endl);
403
405 fit.value.print(cout);
406 }
407
409
410 if (writeFits) {
411
413 hn.Fill(fit.numberOfIterations);
414 hr.Fill(fit.value.R );
415 h1.Fill(fit.value.p1);
416 h2.Fill(fit.value.p2);
417 h3.Fill(fit.value.p3);
418 h4.Fill(fit.value.p4);
419 hc.Fill(fit.value.cc);
420 hb.Fill(fit.value.bc);
421
422 TH1D h1t(
MAKE_CSTRING(module->getID() <<
".1t0"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
423 TH1D h1s(
MAKE_CSTRING(module->getID() <<
".1TTS"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
424 TH1D h1q(
MAKE_CSTRING(module->getID() <<
".1QE"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
425
427 h1t.SetBinContent(pmt + 1, fit.value.parameters[pmt].t0 ());
428 h1t.SetBinError (pmt + 1, fit.error.parameters[pmt].t0 () + numeric_limits<double>::epsilon());
429 h1s.SetBinContent(pmt + 1, fit.value.parameters[pmt].TTS());
430 h1s.SetBinError (pmt + 1, fit.error.parameters[pmt].TTS() + numeric_limits<double>::epsilon());
431 h1q.SetBinContent(pmt + 1, fit.value.parameters[pmt].QE ());
432 h1q.SetBinError (pmt + 1, fit.error.parameters[pmt].QE () + numeric_limits<double>::epsilon());
433 }
434
435 out << h1t << h1s << h1q;
436
437 for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
438
440
441 for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
442
443 const double dt_ns = h2d->GetYaxis()->GetBinCenter(iy);
444
445 h2d->SetBinContent(ix, iy, fit.value.getValue(
pair, dt_ns));
446 h2d->SetBinError (ix, iy, 0.0);
447 }
448 }
449
451 h2d->Write();
452
453 const double x =
string.getIndex(module->getString());
454 const double y =
module->getFloor();
455
457 HN->Fill(x, y, fit.numberOfIterations);
458 }
459
460 const double t0 = (fit.value.hasFixedTimeOffset() ? fit.value.getFixedTimeOffset() : 0.0);
461
463
465
467
468 if (R > 0.0)
469 data.QE = fit.value.parameters[pmt].QE() / R;
470 else
472
473 if (writeFits > 1) {
474 data.TTS_ns = fit.value.parameters[pmt].TTS();
475 }
476
477 module->getPMT(pmt).addT0(fit.value.parameters[pmt].t0() - t0);
478 }
479 }
480 catch(const exception& error) {
481
482 ERROR(
"Module " << setw(10) << module->getID() <<
' ' << error.what() <<
" -> set QEs to 0." << endl);
483
486 }
487 }
488 }
489
490
492
493 {
495 JSTDObjectWriter <JMeta> writer(meta);
496
497 writer << reader;
498 }
499
500 for (vector<JMeta>::const_reverse_iterator i = meta.rbegin(); i != meta.rend(); ++i) {
503 }
504
505 if (overwriteDetector) {
506
507 NOTICE(
"Store calibration data on file " << detectorFile << endl);
508
510 }
511
512 if (pmtFile != "") {
513 parameters.
store(pmtFile.c_str());
514 }
515
516
517 for (vector<JMeta>::const_iterator i = meta.begin(); i != meta.end(); ++i) {
519 }
520
523 }
524
525 if (writeFits) {
526 out << h0 << hn << hr << h1 << h2 << h3 << h4 << hc << hb << H2 << *HN;
527 }
528
529 out.Close();
530
531 return 0;
532}
#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
#define TEST(T)
Test data member.
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
static double BELL_SHAPE
Bell shape.
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)...