72 using namespace KM3NETDAQ;
78 const double epsilon = 1.0e-10;
87 bool overwriteDetector;
111 JParser<> zap(
"Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
114 zap[
'f'] =
make_field(inputFile,
"input file (output from JMergeCalibrateK40).");
116 zap[
'a'] =
make_field(detectorFile,
"detector file.");
117 zap[
'P'] =
make_field(pmtFile,
"specify PMT file name that can be given to JTriggerEfficiency.") =
"";
119 "Fix time offset(s) of PMT(s) of certain module(s), e.g."
120 "\n-! \"808969848 0 808982077 23\" will fix offset of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
121 "\nSame PMT offset can be fixed for all optical modules, e.g."
122 "\n-! \"-1 0 -1 23\" will fix offsets of PMTs 0 and 23 for all optical modules.")
124 zap[
'r'] =
make_field(reverse,
"reverse TDC constraints due to option -! <TDC>.");
125 zap[
'A'] =
make_field(overwriteDetector,
"overwrite detector file provided through '-a' with correct time offsets.");
126 zap[
'w'] =
make_field(writeFits,
"write fit results.");
127 zap[
'D'] =
make_field(fitAngDist,
"fit angular distribution; fix normalisation.");
128 zap[
'M'] =
make_field(fitModel,
"fit angular distribution; fix PMT QEs = 1.0.");
129 zap[
'E'] =
make_field(mu,
"expectation value for npe given two-fold coincidence") = 0.0;
132 zap[
'O'] =
make_field(option,
"ROOT fit option, see TH1::Fit.") =
"";
137 catch(
const exception &error) {
138 FATAL(error.what() << endl);
142 ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls(1000000);
176 const range_type range = TDC.equal_range(-1);
178 if (range.first != range.second) {
180 const JTDC_t buffer(range.first, range.second);
182 TDC.erase(range.first, range.second);
184 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
186 for (JTDC_t::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
192 if (i->second == -1) {
195 TDC.insert(make_pair(module->getID(),
pmt));
200 TDC.insert(make_pair(module->getID(), i->second));
215 for (JTDC_t::const_iterator __p = TDC.begin(); __p != TDC.end(); ) {
217 JTDC_t::const_iterator __q = __p;
219 for ( ; __q != TDC.end() && __q->first == __p->first; ++__q) {}
223 JTDC_t::const_iterator __i = __p;
225 for ( ; __i != __q && __i->second != i; ++__i) {}
228 buffer.insert(std::make_pair(__p->first,i));
238 for (JTDC_t::const_iterator tdc = TDC.begin(); tdc != TDC.end(); ++tdc) {
239 DEBUG(
"PMT " << setw(10) << tdc->first <<
' ' << setw(2) << tdc->second <<
" constrain t0." << endl);
242 for (JTDC_t::const_iterator tdc = TDC.begin(); tdc != TDC.end(); ++tdc) {
243 if (tdc->first < 0) {
244 FATAL(
"Illegal module identifier: " << tdc->first << endl);
247 FATAL(
"Illegal TDC (= readout channel) identifier: " << tdc->second <<
" {0, .., " <<
NUMBER_OF_PMTS - 1 <<
"}" << endl);
252 if (option.find(
'O') == string::npos) { option +=
'0'; }
253 if (option.find(
'S') == string::npos) { option +=
'S'; }
257 TFile*
in = TFile::Open(inputFile.c_str(),
"exist");
259 if (
in == NULL || !
in->IsOpen()) {
260 FATAL(
"File: " << inputFile <<
" not opened." << endl);
265 TH1D hc(
"chi2", NULL, 500, 0.0, 5.0);
267 TH1D
p1(
"p1", NULL, 501, -5.0, +5.0);
268 TH1D p2(
"p2", NULL, 501, -5.0, +5.0);
269 TH1D p3(
"p3", NULL, 501, -5.0, +5.0);
270 TH1D p4(
"p4", NULL, 501, -5.0, +5.0);
271 TH1D bg(
"bg", NULL, 501, -0.1, +0.1);
272 TH1D cc(
"cc", NULL, 501, -0.1, +0.1);
281 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
283 NOTICE(
"Module " << setw(10) << module->getID() << endl);
291 if (h2s == NULL || h2s->GetEntries() == 0) {
293 NOTICE(
"No data for module " << module->getID() <<
"; set corresponding QEs to 0." << endl);
306 const range_type range = TDC.equal_range(module->getID());
313 h2s->GetXaxis()->GetXmin(),
314 h2s->GetXaxis()->GetXmax(),
315 h2s->GetYaxis()->GetXmin(),
316 h2s->GetYaxis()->GetXmax(),
317 range.first == range.second);
321 for (JTDC_t::const_iterator i = range.first; i != range.second; ++i) {
322 fixParameter(fit, fit.getModelParameter(i->second, &JPMTParameters_t::t0));
332 fixParameter(fit, fit.getModelParameter(&JFitK40::Rate_Hz));
348 for (
int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
350 Double_t value = 0.0;
351 Double_t error = 0.0;
353 for (
int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
354 value += h2s->GetBinContent(ix,iy);
355 error += h2s->GetBinError(ix,iy) * h2s->GetBinError(ix,iy);
360 if (value < MINIMAL_RATE_HZ + STDEV*error) {
364 buffer[pair.first] += 1;
365 buffer[pair.second] += 1;
375 WARNING(
"PMT " << setw(10) << module->getID() <<
' ' << setw(2) <<
pmt <<
" too low a rate; switch off." << endl);
383 NOTICE(
"Skip fit for module " << module->getID() <<
"; set corresponding QEs to 0." << endl);
395 h2s->GetXaxis()->SetRangeUser(max(X.getLowerLimit(), h2s->GetXaxis()->GetXmin()),
396 min(X.getUpperLimit(), h2s->GetXaxis()->GetXmax()));
399 h2s->GetYaxis()->SetRangeUser(max(Y.getLowerLimit(), h2s->GetYaxis()->GetXmin()),
400 min(Y.getUpperLimit(), h2s->GetYaxis()->GetXmax()));
404 TFitResultPtr
result = fit(*h2s, option);
416 if (!result->IsParameterFixed(fit.getModelParameter(
pmt, &JPMTParameters_t::QE)) && !
is_valid(values.getQE(
pmt), errors.getQE(
pmt))) {
418 WARNING(
"PMT " << setw(10) << module->getID() <<
' ' << setw(2) <<
pmt <<
' '
420 <<
FIXED(5,2) << values.getQE(
pmt) <<
" +/- "
422 <<
" compatible with zero -> set QE = 0 and t0 = 0; -> refit." << endl);
433 NOTICE(
"Skip fit for module " << module->getID() <<
"; set corresponding QEs to 0." << endl);
444 result = fit(*h2s, option);
451 hc.Fill(TMath::Log10(result->Chi2() / result->Ndf()));
458 h2s->GetXaxis()->GetNbins(), h2s->GetXaxis()->GetXmin (), h2s->GetXaxis()->GetXmax (),
459 h2s->GetYaxis()->GetNbins(), h2s->GetYaxis()->GetXmin (), h2s->GetYaxis()->GetXmax ());
461 for (
int ix = 1; ix <= h2f.GetXaxis()->GetNbins(); ++ix) {
462 for (
int iy = 1; iy <= h2f.GetYaxis()->GetNbins(); ++iy) {
464 const Double_t x = h2f.GetXaxis()->GetBinCenter(ix);
465 const Double_t y = h2f.GetYaxis()->GetBinCenter(iy);
467 h2f.SetBinContent(ix, iy, fit.Eval(x,y));
468 h2f.SetBinError (ix, iy, 0.0);
475 cout <<
"Fit result chi2 / NDF " << result->Chi2() <<
' ' << result->Ndf() <<
' ' << (result->IsValid() ?
"" :
"failed") << endl;
477 cout <<
"Rate_Hz= " <<
FIXED(8,4) << values.Rate_Hz << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::Rate_Hz)) ?
" *** fixed *** " :
"") << endl;
479 cout <<
"p1= " <<
FIXED(8,4) << values.p1 << (result->IsParameterFixed(fit.getModelParameter(&
JFitK40::p1)) ?
" *** fixed *** " :
"") << endl;
480 cout <<
"p2= " <<
FIXED(8,4) << values.p2 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p2)) ?
" *** fixed *** " :
"") << endl;
481 cout <<
"p3= " <<
FIXED(8,4) << values.p3 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p2)) ?
" *** fixed *** " :
"") << endl;
482 cout <<
"p4= " <<
FIXED(8,4) << values.p4 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p3)) ?
" *** fixed *** " :
"") << endl;
484 cout <<
"Background (correlated): " <<
FIXED(8,5) << values.bg << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::bg)) ?
" *** fixed *** " :
"") << endl;
485 cout <<
"Background (uncorrelated): " <<
FIXED(8,5) << values.cc << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::cc)) ?
" *** fixed *** " :
"") << endl;
487 cout <<
" PMT t0 [ns] TTS [ns] R" << endl;
494 << setw(2) <<
pmt <<
' '
495 <<
FIXED(7,2) << fit.getT0 (
pmt) <<
' '
496 <<
FIXED(7,2) << fit.getTTS(
pmt) <<
' '
499 cout << (result->IsParameterFixed(fit.getModelParameter(
pmt, &JPMTParameters_t::QE)) ?
" *** fixed QE *** " :
"");
500 cout << (result->IsParameterFixed(fit.getModelParameter(
pmt, &JPMTParameters_t::TTS)) ?
" *** fixed TTS *** " :
"");
501 cout << (result->IsParameterFixed(fit.getModelParameter(
pmt, &JPMTParameters_t::t0)) ?
" *** fixed t0 *** " :
"");
502 cout << (fit.is_disabled (
pmt) ?
" (disabled)" :
"");
503 cout << (fit.is_average_t0(
pmt) ?
" (averaged)" :
"");
506 Q[0].
put(fit.getT0 (
pmt));
507 Q[1].
put(fit.getTTS(
pmt));
508 Q[2].
put(fit.getQE (
pmt));
545 h1t.SetBinContent(
pmt + 1, values.getT0 (
pmt));
546 h1s.SetBinContent(
pmt + 1, values.getTTS(
pmt));
547 h1q.SetBinContent(
pmt + 1, values.getQE (
pmt));
549 h1t.SetBinError (
pmt + 1, max(errors.getT0 (
pmt), epsilon));
550 h1s.SetBinError (
pmt + 1, max(errors.getTTS(
pmt), epsilon));
551 h1q.SetBinError (
pmt + 1, max(errors.getQE (
pmt), epsilon));
554 out << h1t << h1s << h1q;
564 if (overwriteDetector) {
565 module->getPMT(
pmt).addT0(fit.getT0(
pmt));
572 out << hc <<
p1 << p2 << p3 << p4 << bg << cc;
578 if (overwriteDetector) {
580 NOTICE(
"Store calibration data on file " << detectorFile << endl);
600 ofstream out(pmtFile.c_str());
Utility class to parse command line options.
Parametrisation of two-fold coincidence rate due to K40 and other radioactive decays.
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
esac print_variable DETECTOR INPUT_FILE OUTPUT_FILE CDF for TYPE in
Utility class to parse parameter values.
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
#define MAKE_CSTRING(A)
Make C-string.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Auxiliary data structure for floating point format specification.
esac $JPP_DIR examples JDetector JTransitTime o $OUTPUT_FILE n N $NPE T $TTS_NS d $DEBUG for HISTOGRAM in tts tt2 pmt
static const char *const _2F
Name extension for 2F rate.
void setModelParameters(const Double_t *data)
Set model parameters.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
bool is_valid(const json &js)
Check validity of JSon data.
static const char *const _2R
Name extension for 2D rate.
Auxiliary class for map of PMT parameters.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
JRange< Double_t > JRange_t
bool fixParameter(TF1 &f1, const JFitParameter_t ¶meter)
Fix fit parameter.
void store(const JString &file_name, const JDetector &detector)
Store detector to output file.
virtual const char * what() const
Get error message.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Fit parameters for two-fold coincidence rate due to K40.
#define DEBUG(A)
Message macros.