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
106
107 JParser<> zap(
"Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
108
110 zap[
'f'] =
make_field(inputFile,
"input file (output from JMergeCalibrateK40).");
112 zap[
'a'] =
make_field(detectorFile,
"detector file.");
113 zap[
'P'] =
make_field(pmtFile,
"specify PMT file name that can be given to JTriggerEfficiency.") =
"";
115 "Fix time offset(s) of PMT(s) of certain module(s), e.g."
116 "\n-! \"808969848 0 808982077 23\" will fix time offsets of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
117 "\nSame PMT offset can be fixed for all optical modules, e.g."
118 "\n-! \"-1 0 -1 23\" will fix time offsets of PMTs 0 and 23 of all optical modules.") =
JPARSER::initialised();
119 zap[
'r'] =
make_field(reverse,
"reverse TDC constraints due to option -! <TDC>.");
120 zap[
'A'] =
make_field(overwriteDetector,
"overwrite detector file provided through '-a' with fitted time offsets.");
121 zap[
'w'] =
make_field(writeFits,
"write fit results to ROOT file; -ww also write fitted TTS to PMT file.");
122 zap[
'D'] =
make_field(fitAngle,
"fit angular distribution; fix normalisation.");
123 zap[
'B'] =
make_field(fitNoise,
"fit background.");
124 zap[
'M'] =
make_field(fitModel,
"fit angular distribution as well as normalisation; fix PMT QEs = 1.0.");
125 zap[
'Q'] =
make_field(fixQE,
"fix PMT QEs = 1.0.");
129
130 zap(argc, argv);
131 }
132 catch(const exception &error) {
133 FATAL(error.what() << endl);
134 }
135
136
137 if ((fitModel ? 1 : 0) +
138 (fitAngle ? 1 : 0) +
139 (fitNoise ? 1 : 0) +
140 (fixQE ? 1 : 0) > 1) {
141 FATAL(
"Use either option -M, -D, -B or -Q" << endl);
142 }
143
149
150 if (reverse) {
152 }
153
154 for (JTDC_t::const_iterator i = TDC.begin(); i != TDC.end(); ++i) {
155 DEBUG(
"PMT " << setw(10) << i->first <<
' ' << setw(2) << i->second <<
" constrain t0." << endl);
156 }
157
158 try {
160 }
161 catch(const exception &error) {
162 FATAL(error.what() << endl);
163 }
164
166
167 try {
169 }
172 }
173
174
176
177 if (pmtFile != "") {
178 try {
179 parameters.
load(pmtFile.c_str());
180 }
181 catch(const exception& error) {}
182 }
183
185
186
187 TFile* in = TFile::Open(inputFile.c_str(), "exist");
188
189 if (in == NULL || !in->IsOpen()) {
190 FATAL(
"File: " << inputFile <<
" not opened." << endl);
191 }
192
193
195
196 TH1D h0("chi2", NULL, 500, 0.0, 5.0);
197 TH1D hn("hn", NULL, 501, -0.5, 500.0);
198 TH1D hr("rate", NULL, 500, 0.0, 25.0);
199 TH1D h1("p1", NULL, 500, -5.0, +5.0);
200 TH1D h2("p2", NULL, 500, -5.0, +5.0);
201 TH1D h3("p3", NULL, 500, -5.0, +5.0);
202 TH1D h4("p4", NULL, 500, -5.0, +5.0);
203 TH1D hc("cc", NULL, 500, -0.1, +0.1);
204 TH1D hb("bc", NULL, 500, -0.1, +0.1);
205
208
209 TH2D H2("detector", NULL,
210 string.size() + 0, -0.5, string.size() - 0.5,
212
213 for (Int_t i = 1; i <= H2.GetXaxis()->GetNbins(); ++i) {
214 H2.GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(
string.at(i-1)));
215 }
216 for (Int_t i = 1; i <= H2.GetYaxis()->GetNbins(); ++i) {
218 }
219
220 TH2D* HN = (TH2D*) H2.Clone("iterations");
221
223
224 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
225
226 if (module->getFloor() == 0) {
227 continue;
228 }
229
231
232 NOTICE(
"Module " << setw(10) << module->getID() <<
' ' <<
getLabel(module->getLocation()) <<
" !" <<
distance(range.first, range.second) << endl);
233
235
236 if (h2d == NULL || h2d->GetEntries() == 0) {
237
238 NOTICE(
"No data for module " << module->getID() <<
" -> set QEs to 0." << endl);
239
242 }
243
244 continue;
245 }
246
248
250
253
254 for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
255
257
259
260 double V = 0.0;
261 double W = 0.0;
262
263 for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
264
265 const double x = h2d->GetXaxis()->GetBinCenter(ix);
266 const double y = h2d->GetYaxis()->GetBinCenter(iy);
267
268 if (X(x) && Y(y)) {
269
270 double value = h2d->GetBinContent(ix,iy);
271 double error = h2d->GetBinError (ix,iy);
272
273 buffer.push_back(
rate_type(y, value, error));
274
275 double width = h2d->GetYaxis()->GetBinWidth(iy);
276
277 value *= width;
278 error *= width;
279
280 V += value;
281 W += error * error;
282 }
283 }
284
285 W = sqrt(W);
286
287 if (V <= 0.0 - STDEV*W) {
288 count[0][
pair.first] += 1;
289 count[0][
pair.second] += 1;
290 }
291
292 if (V <= MINIMAL_RATE_HZ + STDEV*W) {
293 count[1][
pair.first] += 1;
294 count[1][
pair.second] += 1;
295 }
296 }
297
299
300 if (count[0][pmt] >= MAXIMAL_COUNTS) {
301
302 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
" some rates negative -> fit background" << endl);
303
304 if (fit.value.parameters[pmt].status) {
305 model.parameters[pmt].bg.set();
306 }
307 }
308
309 if (count[1][pmt] == NUMBER_OF_PMTS) {
310
311 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
" all rates to low -> disable" << endl);
312
313 model.parameters[pmt].disable();
314 }
315 }
316
317 DEBUG(
"Start value:" << endl <<
model << endl);
318
319 try {
320
322
324
326
327 ERROR(
"Fit result " << setw(10) << module->getID() <<
" NDF " << setw(5) <<
result.ndf <<
" -> skip" << endl);
328
329 continue;
330 }
331
332 bool refit = false;
333
335
336 if (fit.value.parameters[pmt].status) {
337
338 if (fit.value.parameters[pmt].QE() <= QE_MIN ||
339 fit.value.parameters[pmt].QE() <= 0.0 + STDEV * fit.error.parameters[pmt].QE()) {
340
341 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
' '
342 << "QE = "
343 <<
FIXED(5,3) << fit.value.parameters[pmt].QE() <<
" +/- "
344 <<
FIXED(5,3) << fit.error.parameters[pmt].QE() <<
" "
345 << " -> disable" << (!refit ? " and refit" : "") << endl);
346
347 fit.value.parameters[pmt].disable();
348
349 refit = true;
350 }
351
352 if (fit.value.parameters[pmt].t0.atLimit(T0_NS)) {
353
354 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
' '
355 << "t0 at limit "
356 <<
FIXED(5,3) << fit.value.parameters[pmt].t0() <<
" +/- "
357 <<
FIXED(5,3) << fit.error.parameters[pmt].t0() <<
" "
358 << " -> disable" << (!refit ? " and refit" : "") << endl);
359
360 fit.value.parameters[pmt].disable();
361
362 refit = true;
363 }
364 }
365 }
366
367 if (refit) {
368
370 if (fit.value.parameters[pmt].status) {
372 }
373 }
374
375 refit = false;
377 }
378
379 NOTICE(
"Fit result " << setw(10) << module->getID() <<
" chi2 / NDF " <<
FIXED(10,2) <<
result.chi2 <<
" / " << setw(5) <<
result.ndf <<
' ' << setw(5) << fit.numberOfIterations << endl);
380
383 }
384
386
387 if (writeFits) {
388
390 hn.Fill(fit.numberOfIterations);
391 hr.Fill(fit.value.R );
392 h1.Fill(fit.value.p1);
393 h2.Fill(fit.value.p2);
394 h3.Fill(fit.value.p3);
395 h4.Fill(fit.value.p4);
396 hc.Fill(fit.value.cc);
397 hb.Fill(fit.value.bc);
398
399 TH1D h1t(
MAKE_CSTRING(module->getID() <<
".1t0"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
400 TH1D h1s(
MAKE_CSTRING(module->getID() <<
".1TTS"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
401 TH1D h1q(
MAKE_CSTRING(module->getID() <<
".1QE"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
402
404 h1t.SetBinContent(pmt + 1, fit.value.parameters[pmt].t0 ());
405 h1t.SetBinError (pmt + 1, fit.error.parameters[pmt].t0 () + numeric_limits<double>::epsilon());
406 h1s.SetBinContent(pmt + 1, fit.value.parameters[pmt].TTS());
407 h1s.SetBinError (pmt + 1, fit.error.parameters[pmt].TTS() + numeric_limits<double>::epsilon());
408 h1q.SetBinContent(pmt + 1, fit.value.parameters[pmt].QE ());
409 h1q.SetBinError (pmt + 1, fit.error.parameters[pmt].QE () + numeric_limits<double>::epsilon());
410 }
411
412 out << h1t << h1s << h1q;
413
414 for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
415
417
418 for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
419
420 const double dt_ns = h2d->GetYaxis()->GetBinCenter(iy);
421
422 h2d->SetBinContent(ix, iy, fit.value.getValue(
pair, dt_ns));
423 h2d->SetBinError (ix, iy, 0.0);
424 }
425 }
426
428 h2d->Write();
429
430 const double x =
string.getIndex(module->getString());
431 const double y =
module->getFloor();
432
434 HN->Fill(x, y, fit.numberOfIterations);
435 }
436
437 const double t0 = (fit.value.hasFixedTimeOffset() ? fit.value.getFixedTimeOffset() : 0.0);
438
440
442
444
445 if (P > 0.0)
446 data.QE = fit.value.parameters[pmt].QE / P;
447 else
449
450 if (writeFits > 1) {
451 data.TTS_ns = fit.value.parameters[pmt].TTS();
452 }
453
454 module->getPMT(pmt).addT0(fit.value.parameters[pmt].t0.get() - t0);
455 }
456 }
457 catch(const exception& error) {
458
459 ERROR(
"Module " << setw(10) << module->getID() <<
' ' << error.what() <<
" -> set QEs to 0." << endl);
460
463 }
464 }
465 }
466
467
469
470 {
472 JSTDObjectWriter <JMeta> writer(meta);
473
474 writer << reader;
475 }
476
477 for (vector<JMeta>::const_reverse_iterator i = meta.rbegin(); i != meta.rend(); ++i) {
480 }
481
482 if (overwriteDetector) {
483
484 NOTICE(
"Store calibration data on file " << detectorFile << endl);
485
487 }
488
489 if (pmtFile != "") {
490 parameters.
store(pmtFile.c_str());
491 }
492
493
494 for (vector<JMeta>::const_iterator i = meta.begin(); i != meta.end(); ++i) {
496 }
497
500 }
501
502 if (writeFits) {
503 out << h0 << hn << hr << h1 << h2 << h3 << h4 << hc << hb << H2 << *HN;
504 }
505
506 out.Close();
507
508 return 0;
509}
#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
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)...