63int main(
int argc,
char **argv)
79 bool overwriteDetector;
104 JParser<> zap(
"Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
107 zap[
'f'] =
make_field(inputFile,
"input file (output from JMergeCalibrateK40).");
109 zap[
'a'] =
make_field(detectorFile,
"detector file.");
110 zap[
'P'] =
make_field(pmtFile,
"specify PMT file name that can be given to JTriggerEfficiency.") =
"";
112 "Fix time offset(s) of PMT(s) of certain module(s), e.g."
113 "\n-! \"808969848 0 808982077 23\" will fix time offsets of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
114 "\nSame PMT offset can be fixed for all optical modules, e.g."
115 "\n-! \"-1 0 -1 23\" will fix time offsets of PMTs 0 and 23 of all optical modules.") =
JPARSER::initialised();
116 zap[
'r'] =
make_field(reverse,
"reverse TDC constraints due to option -! <TDC>.");
117 zap[
'A'] =
make_field(overwriteDetector,
"overwrite detector file provided through '-a' with fitted time offsets.");
118 zap[
'w'] =
make_field(writeFits,
"write fit results to ROOT file; -ww also write fitted TTS to PMT file.");
119 zap[
'D'] =
make_field(fitAngle,
"fit angular distribution; fix normalisation.");
120 zap[
'B'] =
make_field(fitNoise,
"fit background.");
121 zap[
'M'] =
make_field(fitModel,
"fit angular distribution as well as normalisation; fix PMT QEs = 1.0.");
122 zap[
'Q'] =
make_field(fixQE,
"fix PMT QEs = 1.0.");
129 catch(
const exception &error) {
130 FATAL(error.what() << endl);
134 if ((fitModel ? 1 : 0) +
137 (fixQE ? 1 : 0) > 1) {
138 FATAL(
"Use either option -M, -D, -B or -Q" << endl);
151 for (JTDC_t::const_iterator i = TDC.begin(); i != TDC.end(); ++i) {
152 DEBUG(
"PMT " << setw(10) << i->first <<
' ' << setw(2) << i->second <<
" constrain t0." << endl);
158 catch(
const exception &error) {
159 FATAL(error.what() << endl);
176 parameters.
load(pmtFile.c_str());
178 catch(
const exception& error) {}
184 TFile* in = TFile::Open(inputFile.c_str(),
"exist");
186 if (in == NULL || !in->IsOpen()) {
187 FATAL(
"File: " << inputFile <<
" not opened." << endl);
193 TH1D h0(
"chi2", NULL, 500, 0.0, 5.0);
194 TH1D hn(
"hn", NULL, 501, -0.5, 500.0);
195 TH1D hr(
"rate", NULL, 500, 0.0, 25.0);
196 TH1D h1(
"p1", NULL, 500, -5.0, +5.0);
197 TH1D h2(
"p2", NULL, 500, -5.0, +5.0);
198 TH1D h3(
"p3", NULL, 500, -5.0, +5.0);
199 TH1D h4(
"p4", NULL, 500, -5.0, +5.0);
200 TH1D hc(
"cc", NULL, 500, -0.1, +0.1);
205 TH2D H2(
"detector", NULL,
206 string.size() + 0, -0.5,
string.size() - 0.5,
209 for (Int_t i = 1; i <= H2.GetXaxis()->GetNbins(); ++i) {
210 H2.GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(
string.at(i-1)));
212 for (Int_t i = 1; i <= H2.GetYaxis()->GetNbins(); ++i) {
216 TH2D* HN = (TH2D*) H2.Clone(
"iterations");
220 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
222 if (module->getFloor() == 0) {
228 NOTICE(
"Module " << setw(10) << module->getID() <<
' ' <<
getLabel(module->getLocation()) <<
" !" <<
distance(range.first, range.second) << endl);
232 if (h2d == NULL || h2d->GetEntries() == 0) {
234 NOTICE(
"No data for module " << module->getID() <<
" -> set QEs to 0." << endl);
250 for (
int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
254 auto& buffer = data[
pair];
259 for (
int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
261 const double x = h2d->GetXaxis()->GetBinCenter(ix);
262 const double y = h2d->GetYaxis()->GetBinCenter(iy);
266 double value = h2d->GetBinContent(ix,iy);
267 double error = h2d->GetBinError (ix,iy);
269 buffer.push_back(
rate_type(y, value, error));
271 double width = h2d->GetYaxis()->GetBinWidth(iy);
283 if (V <= 0.0 - STDEV*W) {
284 count[0][
pair.first] += 1;
285 count[0][
pair.second] += 1;
288 if (V <= MINIMAL_RATE_HZ + STDEV*W) {
289 count[1][
pair.first] += 1;
290 count[1][
pair.second] += 1;
296 if (count[0][pmt] >= MAXIMAL_COUNTS) {
298 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
" some rates negative -> fit background" << endl);
300 if (fit.value.parameters[pmt].status) {
301 model.parameters[pmt].bg.set();
305 if (count[1][pmt] == NUMBER_OF_PMTS) {
307 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
" all rates to low -> disable" << endl);
309 model.parameters[pmt].disable();
313 DEBUG(
"Start value:" << endl <<
model << endl);
323 ERROR(
"Fit result " << setw(10) << module->getID() <<
" NDF " << setw(5) <<
result.ndf <<
" -> skip" << endl);
332 if (fit.value.parameters[pmt].status) {
334 if (fit.value.parameters[pmt].QE() <= QE_MIN ||
335 fit.value.parameters[pmt].QE() <= 0.0 + STDEV * fit.error.parameters[pmt].QE()) {
337 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
' '
339 <<
FIXED(5,3) << fit.value.parameters[pmt].QE() <<
" +/- "
340 <<
FIXED(5,3) << fit.error.parameters[pmt].QE() <<
" "
341 <<
" -> disable" << (!refit ?
" and refit" :
"") << endl);
343 fit.value.parameters[pmt].disable();
348 if (fit.value.parameters[pmt].t0.atLimit(T0_NS)) {
350 WARNING(
"PMT " << setw(10) << module->getID() <<
'.' <<
FILL(2,
'0') << pmt <<
FILL() <<
' '
352 <<
FIXED(5,3) << fit.value.parameters[pmt].t0() <<
" +/- "
353 <<
FIXED(5,3) << fit.error.parameters[pmt].t0() <<
" "
354 <<
" -> disable" << (!refit ?
" and refit" :
"") << endl);
356 fit.value.parameters[pmt].disable();
366 if (fit.value.parameters[pmt].status) {
375 NOTICE(
"Fit result " << setw(10) << module->getID() <<
" chi2 / NDF " <<
FIXED(10,2) <<
result.chi2 <<
" / " << setw(5) <<
result.ndf <<
' ' << setw(5) << fit.numberOfIterations << endl);
382 hn.Fill(fit.numberOfIterations);
383 hr.Fill(fit.value.R );
384 h1.Fill(fit.value.p1);
385 h2.Fill(fit.value.p2);
386 h3.Fill(fit.value.p3);
387 h4.Fill(fit.value.p4);
388 hc.Fill(fit.value.cc);
390 TH1D h1t(
MAKE_CSTRING(module->getID() <<
".1t0"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
391 TH1D h1s(
MAKE_CSTRING(module->getID() <<
".1TTS"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
392 TH1D h1q(
MAKE_CSTRING(module->getID() <<
".1QE"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
395 h1t.SetBinContent(pmt + 1, fit.value.parameters[pmt].t0 ());
396 h1t.SetBinError (pmt + 1, fit.error.parameters[pmt].t0 () + numeric_limits<double>::epsilon());
397 h1s.SetBinContent(pmt + 1, fit.value.parameters[pmt].TTS());
398 h1s.SetBinError (pmt + 1, fit.error.parameters[pmt].TTS() + numeric_limits<double>::epsilon());
399 h1q.SetBinContent(pmt + 1, fit.value.parameters[pmt].QE ());
400 h1q.SetBinError (pmt + 1, fit.error.parameters[pmt].QE () + numeric_limits<double>::epsilon());
403 out << h1t << h1s << h1q;
405 for (
int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
409 for (
int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
411 const double dt_ns = h2d->GetYaxis()->GetBinCenter(iy);
413 h2d->SetBinContent(ix, iy, fit.value.getValue(
pair, dt_ns));
414 h2d->SetBinError (ix, iy, 0.0);
421 const double x =
string.getIndex(module->getString());
422 const double y =
module->getFloor();
425 HN->Fill(x, y, fit.numberOfIterations);
428 const double t0 = (fit.value.hasFixedTimeOffset() ? fit.value.getFixedTimeOffset() : 0.0);
437 data.QE = fit.value.parameters[pmt].QE / P;
442 data.TTS_ns = fit.value.parameters[pmt].TTS();
445 module->getPMT(pmt).addT0(fit.value.parameters[pmt].t0.get() - t0);
448 catch(
const exception& error) {
450 ERROR(
"Module " << setw(10) << module->getID() <<
' ' << error.what() <<
" -> set QEs to 0." << endl);
463 JSTDObjectWriter <JMeta> writer(meta);
468 for (vector<JMeta>::const_reverse_iterator i = meta.rbegin(); i != meta.rend(); ++i) {
473 if (overwriteDetector) {
475 NOTICE(
"Store calibration data on file " << detectorFile << endl);
481 parameters.
store(pmtFile.c_str());
485 for (vector<JMeta>::const_iterator i = meta.begin(); i != meta.end(); ++i) {
494 out << h0 << hn << hr << h1 << h2 << h3 << h4 << hc << H2 << *HN;