Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JFitK40.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 #include <map>
6 #include <cmath>
7 
8 #include "TROOT.h"
9 #include "TFile.h"
10 #include "TH1D.h"
11 #include "TH2D.h"
12 #include "TF2.h"
13 #include "TMath.h"
14 #include "TFitResult.h"
15 #include "Math/MinimizerOptions.h"
16 
17 #include "JPhysics/JConstants.hh"
18 #include "km3net-dataformat/online/JDAQ.hh"
19 #include "JDetector/JDetector.hh"
20 #include "JDetector/JDetectorToolkit.hh"
21 #include "JDetector/JPMTParametersMap.hh"
22 
23 #include "JROOT/JRootToolkit.hh"
24 #include "JTools/JRange.hh"
25 #include "JTools/JQuantile.hh"
26 #include "JSupport/JMeta.hh"
27 
28 #include "JCalibrate/JCalibrateK40.hh"
29 #include "JCalibrate/JFitK40.hh"
30 #include "JCalibrate/JTDC_t.hh"
31 
32 #include "Jeep/JProperties.hh"
33 #include "Jeep/JPrint.hh"
34 #include "Jeep/JParser.hh"
35 #include "Jeep/JMessage.hh"
36 
37 namespace {
38 
39  double STDEV = 2.0; //!< number of standard deviations
40  double QE_MIN_VALUE = 0.01; //!< minimal QE value (below, set to zero)
41  double QE_MAX_ERROR = 0.5; //!< maximal QE error (below, set to zero)
42  double MINIMAL_RATE_HZ = 1.0e-2; //!< minimal coincidence rate (bin-by-bin)
43 
44  /**
45  * Check validity of QE measurement.
46  *
47  * The QE measurement is valid when:
48  * - fitted value is at least STDEV standard deviations larger than QE_MIN_VALUE; and
49  * - fitted error is smaller than QE_MAX_ERROR.
50  *
51  * \param qe_value QE value
52  * \param qe_error QE error
53  * \return true if valid; else false
54  */
55  inline bool is_valid(const double qe_value, const double qe_error)
56  {
57  return (qe_value > QE_MIN_VALUE + STDEV * qe_error &&
58  qe_error < QE_MAX_ERROR);
59  }
60 }
61 
62 
63 /**
64  * \file
65  *
66  * Auxiliary program to fit PMT parameters from JMergeCalibrateK40.cc output.\n
67  * To force drawing of fit functions bin-by-bin compatible to the corresponding histogram data, \n
68  * the function values will be stored as histograms, rather than as associated function objects. \n
69  * See JCalibrateK40.hh for the corresponding extensions of histogram names.
70  * \author mdejong
71  */
72 int main(int argc, char **argv)
73 {
74  using namespace std;
75  using namespace JPP;
76  using namespace KM3NETDAQ;
77 
78  typedef JRange<double> JRange_t;
79 
80  const double epsilon = 1.0e-10; // minimal error for ROOT plotting
81  JFitK40Parameters fitk40; // setting of internal fit parameters
82 
83  string inputFile;
84  string outputFile;
85  string detectorFile;
86  string pmtFile;
87  JTDC_t TDC;
88  bool reverse;
89  bool overwriteDetector;
90  bool writeFits;
91  bool fitAngDist;
92  bool fitModel;
93  JRange_t X;
94  JRange_t Y;
95  string option;
96  int debug;
97 
98  try {
99 
100  JProperties properties;
101 
102  properties.insert(gmake_property(STDEV));
103  properties.insert(gmake_property(QE_MIN_VALUE));
104  properties.insert(gmake_property(QE_MAX_ERROR));
105  properties.insert(gmake_property(MINIMAL_RATE_HZ));
106  properties.insert(gmake_property(fitk40.Rate_Hz));
107  properties.insert(gmake_property(fitk40.p1));
108  properties.insert(gmake_property(fitk40.p2));
109  properties.insert(gmake_property(fitk40.p3));
110  properties.insert(gmake_property(fitk40.p4));
111 
112  JParser<> zap("Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
113 
114  zap['@'] = make_field(properties) = JPARSER::initialised();
115  zap['f'] = make_field(inputFile, "input file (output from JMergeCalibrateK40).");
116  zap['o'] = make_field(outputFile, "output file.") = "fit.root";
117  zap['a'] = make_field(detectorFile, "detector file.");
118  zap['P'] = make_field(pmtFile, "specify PMT file name that can be given to JTriggerEfficiency.") = "";
119  zap['!'] = make_field(TDC,
120  "Fix time offset(s) of PMT(s) of certain module(s), e.g."
121  "\n-! \"808969848 0 808982077 23\" will fix offset of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
122  "\nSame PMT offset can be fixed for all optical modules, e.g."
123  "\n-! \"-1 0 -1 23\" will fix offsets of PMTs 0 and 23 for all optical modules.")
124  = JPARSER::initialised();
125  zap['r'] = make_field(reverse, "reverse TDC constraints due to option -! <TDC>.");
126  zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with correct time offsets.");
127  zap['w'] = make_field(writeFits, "write fit results.");
128  zap['D'] = make_field(fitAngDist, "fit angular distribution; fix normalisation.");
129  zap['M'] = make_field(fitModel, "fit angular distribution; fix PMT QEs = 1.0.");
130  zap['x'] = make_field(X, "ROOT fit range (PMT pairs).") = JRange_t();
131  zap['y'] = make_field(Y, "ROOT fit range (time residual).") = JRange_t();
132  zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "";
133  zap['d'] = make_field(debug, "debug.") = 1;
134 
135  zap(argc, argv);
136  }
137  catch(const exception &error) {
138  FATAL(error.what() << endl);
139  }
140 
141 
142  ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls(1000000);
143 
144 
145  if (reverse) {
146  TDC.reverse();
147  }
148 
149  for (JTDC_t::const_iterator i = TDC.begin(); i != TDC.end(); ++i) {
150  DEBUG("PMT " << setw(10) << i->first << ' ' << setw(2) << i->second << " constrain t0." << endl);
151  }
152 
153  TDC.is_valid(true);
154 
155 
156  //----------------------------------------------------------
157  // load detector file and PMT parameters
158  //----------------------------------------------------------
159 
160  JDetector detector;
161 
162  try {
163  load(detectorFile, detector);
164  }
165  catch(const JException& error) {
166  FATAL(error);
167  }
168 
169 
170  JPMTParametersMap parameters;
171 
172  if (pmtFile != "") {
173  try {
174  parameters.load(pmtFile.c_str());
175  }
176  catch(const JException& error) {
177  //ERROR(error.what());
178  }
179  }
180 
181  parameters.comment.add(JMeta(argc, argv));
182 
183 
184  if (option.find('O') == string::npos) { option += '0'; }
185  if (option.find('S') == string::npos) { option += 'S'; }
186  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
187 
188 
189  TFile* in = TFile::Open(inputFile.c_str(), "exist");
190 
191  if (in == NULL || !in->IsOpen()) {
192  FATAL("File: " << inputFile << " not opened." << endl);
193  }
194 
195 
196  TDirectory::AddDirectory(NULL);
197 
198  TFile out(outputFile.c_str(), "recreate");
199 
200  // Histograms with global fit results.
201 
202  TH1D hc("chi2", NULL, 500, 0.0, 5.0);
203 
204  TH1D p1("p1", NULL, 501, -5.0, +5.0);
205  TH1D p2("p2", NULL, 501, -5.0, +5.0);
206  TH1D p3("p3", NULL, 501, -5.0, +5.0);
207  TH1D p4("p4", NULL, 501, -5.0, +5.0);
208  TH1D bg("bg", NULL, 501, -0.1, +0.1);
209  TH1D cc("cc", NULL, 501, -0.1, +0.1);
210 
211 
212  //----------------------------------------------------------
213  // loop over modules in detector file to fit histogram
214  //----------------------------------------------------------
215 
216  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
217 
218  NOTICE("Module " << setw(10) << module->getID() << ' ' << module->getLocation() << endl);
219 
220  if (module->empty()) {
221  continue;
222  }
223 
224  //----------------------------------------------------------
225  // get histogram for this module
226  //----------------------------------------------------------
227 
228  TH2D* h2s = (TH2D*) in->Get(MAKE_CSTRING(module->getID() << _2R));
229 
230  if (h2s == NULL || h2s->GetEntries() == 0) {
231 
232  NOTICE("No data for module " << module->getID() << "; set corresponding QEs to 0." << endl);
233 
234  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
235  parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
236  }
237 
238  continue;
239  }
240 
241  //----------------------------------------------------------
242  // get constraints
243  //----------------------------------------------------------
244 
245  const JTDC_t::range_type range = TDC.equal_range(module->getID());
246 
247  //----------------------------------------------------------
248  // set-up fitting procedure
249  //----------------------------------------------------------
250 
251  JFitK40 fit(*module,
252  h2s->GetXaxis()->GetXmin(),
253  h2s->GetXaxis()->GetXmax(),
254  h2s->GetYaxis()->GetXmin(),
255  h2s->GetYaxis()->GetXmax(),
256  range.first == range.second);
257 
258  fit.setModelParameters(fitk40.getModelParameters());
259 
260  for (JTDC_t::const_iterator i = range.first; i != range.second; ++i) {
261  fixParameter(fit, fit.getModelParameter(i->second, &JPMTParameters_t::t0));
262  }
263 
264  if (fitModel) {
265  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
266  fixParameter(fit, fit.getModelParameter(pmt, &JPMTParameters_t::QE));
267  }
268 
269  } else {
270 
271  fixParameter(fit, fit.getModelParameter(&JFitK40::Rate_Hz));
272 
273  if (!fitAngDist) {
274  fixParameter(fit, fit.getModelParameter(&JFitK40::p1));
275  fixParameter(fit, fit.getModelParameter(&JFitK40::p2));
276  fixParameter(fit, fit.getModelParameter(&JFitK40::p3));
277  fixParameter(fit, fit.getModelParameter(&JFitK40::p4));
278  }
279  }
280 
281  //----------------------------------------------------------
282  // check validity of data
283  //----------------------------------------------------------
284 
285  vector<int> buffer(NUMBER_OF_PMTS, 1);
286 
287  for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
288 
289  Double_t value = 0.0;
290  Double_t error = 0.0;
291 
292  for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
293  value += h2s->GetBinContent(ix,iy);
294  error += h2s->GetBinError(ix,iy) * h2s->GetBinError(ix,iy);
295  }
296 
297  error = sqrt(error);
298 
299  if (value < MINIMAL_RATE_HZ + STDEV*error) {
300 
301  const JFitK40::pair_type pair = fit.getPair(ix-1);
302 
303  buffer[pair.first] += 1;
304  buffer[pair.second] += 1;
305  }
306  }
307 
308  try {
309 
310  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
311 
312  if (buffer[pmt] == NUMBER_OF_PMTS) {
313 
314  WARNING("PMT " << setw(10) << module->getID() << ' ' << setw(2) << pmt << " too low a rate; switch off." << endl);
315 
316  fit.disablePMT(pmt);
317  }
318  }
319  }
320  catch(const JException& error) {
321 
322  NOTICE("Skip fit for module " << module->getID() << "; set corresponding QEs to 0." << endl);
323 
324  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
325  parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
326  }
327 
328  continue;
329  }
330 
331  // constrain fit to specified range(s)
332 
333  if (X != JRange_t()) {
334  h2s->GetXaxis()->SetRangeUser(max(X.getLowerLimit(), h2s->GetXaxis()->GetXmin()),
335  min(X.getUpperLimit(), h2s->GetXaxis()->GetXmax()));
336  }
337  if (Y != JRange_t()) {
338  h2s->GetYaxis()->SetRangeUser(max(Y.getLowerLimit(), h2s->GetYaxis()->GetXmin()),
339  min(Y.getUpperLimit(), h2s->GetYaxis()->GetXmax()));
340  }
341 
342  // 1D-histogram with peak positions
343 
344  TH1D h1(MAKE_CSTRING(module->getID() << _1D), NULL, h2s->GetXaxis()->GetNbins(), h2s->GetXaxis()->GetXmin(), h2s->GetXaxis()->GetXmax());
345 
346  const Double_t sigma[] = {
347 
348  (min(Y.getUpperLimit(), h2s->GetYaxis()->GetXmax()) -
349  max(Y.getLowerLimit(), h2s->GetYaxis()->GetXmin())) / sqrt(12.0),
350 
351  sqrt(JPMTParameters_t().TTS * JPMTParameters_t().TTS * 2.0 + fit.getSigmaK40() * fit.getSigmaK40())
352  };
353 
354 
355  for (Int_t ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
356 
357  Double_t y0 = 0.0;
358  Double_t zmax = 0.0;
359  Double_t value = 0.0;
360  Double_t error = 0.0;
361 
362  for (Int_t iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
363 
364  const Double_t y = h2s->GetYaxis()->GetBinCenter(iy);
365  const Double_t z = h2s->GetBinContent(ix, iy);
366  const Double_t w = h2s->GetBinError (ix, iy);
367 
368  if (Y(y)) {
369 
370  if (z > zmax) {
371  y0 = y;
372  zmax = z;
373  }
374 
375  value += z;
376  error += w*w;
377  }
378  }
379 
380  error = sqrt(error);
381 
382  h1.SetBinContent(ix, y0);
383  h1.SetBinError (ix, value > MINIMAL_RATE_HZ + STDEV * error ? sigma[1] * error / value : sigma[0]);
384  }
385 
386  {
387  JFitK40_t<TF1> f1(fit);
388 
389  TFitResultPtr result = f1(h1, option.c_str());
390 
391  if (writeFits) {
392 
393  TH1D h2(MAKE_CSTRING(module->getID() << _1F), NULL, h1.GetXaxis()->GetNbins(), h1.GetXaxis()->GetXmin(), h1.GetXaxis()->GetXmax());
394 
395  for (int ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
396 
397  const Double_t x = h2.GetXaxis()->GetBinCenter(ix);
398 
399  h2.SetBinContent(ix, f1.Eval(x));
400  h2.SetBinError (ix, 0.0);
401  }
402 
403  out << h1 << h2;
404  }
405 
406  if (debug >= JEEP::notice_t) {
407 
408  cout << "Fit result chi2 / NDF " << result->Chi2() << ' ' << result->Ndf() << ' ' << (result->IsValid() ? "" : "failed") << endl;
409 
410  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
411 
412  cout << ' '
413  << setw(2) << pmt << ' '
414  << FIXED(7,2) << f1.getT0 (pmt) << ' ';
415  cout << (result->IsParameterFixed(f1.getModelParameter(pmt, &JPMTParameters_t::t0)) ? " *** fixed t0 *** " : "");
416  cout << (f1.is_disabled (pmt) ? " (disabled)" : "");
417  cout << (f1.is_average_t0(pmt) ? " (averaged)" : "");
418  cout << endl;
419  }
420  }
421 
422  fit.setModelParameters(f1.getModelParameters());
423  }
424 
425 
426  TFitResultPtr result = fit(*h2s, option);
427 
428 
429  bool refit = false;
430 
431  try {
432 
433  const JFitK40Parameters values(fit.GetParameters());
434  const JFitK40Parameters errors(fit.GetParErrors());
435 
436  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
437 
438  if (!result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::QE)) && !is_valid(values.getQE(pmt), errors.getQE(pmt))) {
439 
440  WARNING("PMT " << setw(10) << module->getID() << ' ' << setw(2) << pmt << ' '
441  << "measured QE "
442  << FIXED(5,2) << values.getQE(pmt) << " +/- "
443  << FIXED(5,2) << errors.getQE(pmt)
444  << " compatible with zero -> set QE = 0 and t0 = 0; -> refit." << endl);
445 
446 
447  fit.disablePMT(pmt);
448 
449  refit = true;
450  }
451  }
452  }
453  catch(const JException& error) {
454 
455  NOTICE("Skip fit for module " << module->getID() << "; set corresponding QEs to 0." << endl);
456 
457  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
458  parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
459  }
460 
461  continue;
462  }
463 
464 
465  if (refit) {
466 
467  JFitK40_t<TF1> f1(fit);
468 
469  f1(h1, option);
470 
471  fit.setModelParameters(f1.getModelParameters());
472 
473  result = fit(*h2s, option);
474  }
475 
476 
477  const JFitK40Parameters values(fit.GetParameters());
478  const JFitK40Parameters errors(fit.GetParErrors());
479 
480  hc.Fill(TMath::Log10(result->Chi2() / result->Ndf()));
481 
482  TH2D h2f(MAKE_CSTRING(module->getID() << _2F),
483  NULL,
484  h2s->GetXaxis()->GetNbins(), h2s->GetXaxis()->GetXmin (), h2s->GetXaxis()->GetXmax (),
485  h2s->GetYaxis()->GetNbins(), h2s->GetYaxis()->GetXmin (), h2s->GetYaxis()->GetXmax ());
486 
487  for (int ix = 1; ix <= h2f.GetXaxis()->GetNbins(); ++ix) {
488  for (int iy = 1; iy <= h2f.GetYaxis()->GetNbins(); ++iy) {
489 
490  const Double_t x = h2f.GetXaxis()->GetBinCenter(ix);
491  const Double_t y = h2f.GetYaxis()->GetBinCenter(iy);
492 
493  h2f.SetBinContent(ix, iy, fit.Eval(x,y));
494  h2f.SetBinError (ix, iy, 0.0);
495  }
496  }
497 
498 
499  if (debug >= JEEP::notice_t) {
500 
501  cout << "Fit result chi2 / NDF " << result->Chi2() << ' ' << result->Ndf() << ' ' << (result->IsValid() ? "" : "failed") << endl;
502 
503  cout << "Rate_Hz= " << FIXED(8,4) << values.Rate_Hz << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::Rate_Hz)) ? " *** fixed *** " : "") << endl;
504 
505  cout << "p1= " << FIXED(8,4) << values.p1 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p1)) ? " *** fixed *** " : "") << endl;
506  cout << "p2= " << FIXED(8,4) << values.p2 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p2)) ? " *** fixed *** " : "") << endl;
507  cout << "p3= " << FIXED(8,4) << values.p3 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p2)) ? " *** fixed *** " : "") << endl;
508  cout << "p4= " << FIXED(8,4) << values.p4 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p3)) ? " *** fixed *** " : "") << endl;
509 
510  cout << "Background (correlated): " << FIXED(8,5) << values.bg << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::bg)) ? " *** fixed *** " : "") << endl;
511  cout << "Background (uncorrelated): " << FIXED(8,5) << values.cc << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::cc)) ? " *** fixed *** " : "") << endl;
512 
513  cout << " PMT t0 [ns] TTS [ns] R" << endl;
514 
515  JQuantile Q[3];
516 
517  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
518 
519  cout << ' '
520  << setw(2) << pmt << ' '
521  << FIXED(7,2) << fit.getT0 (pmt) << ' '
522  << FIXED(7,2) << fit.getTTS(pmt) << ' '
523  << FIXED(7,3) << fit.getQE (pmt);
524 
525  cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::QE)) ? " *** fixed QE *** " : "");
526  cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::TTS)) ? " *** fixed TTS *** " : "");
527  cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::t0)) ? " *** fixed t0 *** " : "");
528  cout << (fit.is_disabled (pmt) ? " (disabled)" : "");
529  cout << (fit.is_average_t0(pmt) ? " (averaged)" : "");
530  cout << endl;
531 
532  Q[0].put(fit.getT0 (pmt));
533  Q[1].put(fit.getTTS(pmt));
534  Q[2].put(fit.getQE (pmt));
535  }
536 
537  cout << "<> "
538  << FIXED(7,2) << Q[0].getMean() << ' '
539  << FIXED(7,2) << Q[1].getMean() << ' '
540  << FIXED(7,3) << Q[2].getMean() << endl;
541  cout << "+/- "
542  << FIXED(7,2) << Q[0].getSTDev() << ' '
543  << FIXED(7,2) << Q[1].getSTDev() << ' '
544  << FIXED(7,3) << Q[2].getSTDev() << endl;
545  }
546 
547  h2s->Write();
548  h2f .Write();
549 
550  if (writeFits) {
551 
552  // Histograms with fit results.
553 
554  p1.Fill(values.p1);
555  p2.Fill(values.p2);
556  p3.Fill(values.p3);
557  p4.Fill(values.p4);
558  bg.Fill(values.bg);
559  cc.Fill(values.cc);
560 
561  TH1D h1t(MAKE_CSTRING(module->getID() << ".1t0"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
562  TH1D h1s(MAKE_CSTRING(module->getID() << ".1TTS"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
563  TH1D h1q(MAKE_CSTRING(module->getID() << ".1QE"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
564 
565  h1t.Sumw2();
566  h1s.Sumw2();
567  h1q.Sumw2();
568 
569  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
570 
571  h1t.SetBinContent(pmt + 1, values.getT0 (pmt));
572  h1s.SetBinContent(pmt + 1, values.getTTS(pmt));
573  h1q.SetBinContent(pmt + 1, values.getQE (pmt));
574 
575  h1t.SetBinError (pmt + 1, max(errors.getT0 (pmt), epsilon));
576  h1s.SetBinError (pmt + 1, max(errors.getTTS(pmt), epsilon));
577  h1q.SetBinError (pmt + 1, max(errors.getQE (pmt), epsilon));
578  }
579 
580  out << h1t << h1s << h1q;
581  }
582 
583 
584  // save output
585 
586  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
587 
588  parameters[JPMTIdentifier(module->getID(), pmt)].QE = values.getQE(pmt);
589 
590  if (overwriteDetector) {
591  module->getPMT(pmt).addT0(fit.getT0(pmt));
592  }
593  }
594  }
595 
596 
597  if (writeFits) {
598  out << hc << p1 << p2 << p3 << p4 << bg << cc;
599  }
600 
601  out.Close();
602 
603 
604  if (overwriteDetector) {
605 
606  NOTICE("Store calibration data on file " << detectorFile << endl);
607 
608  detector.comment.add(JMeta(argc, argv));
609 
610  store(detectorFile, detector);
611  }
612 
613  if (pmtFile != "") {
614  parameters.store(pmtFile.c_str());
615  }
616 
617  return 0;
618 }
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1500
JLANG::JException
General exception.
Definition: JException.hh:23
JCALIBRATE::pair_type
JCombinatorics::pair_type pair_type
Definition: JCalibrateK40.hh:41
JTDC_t.hh
JTOOLS::w
data_type w[N+1][M+1]
Definition: JPolint.hh:741
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
JEEP::debug_t
debug
Definition: JMessage.hh:29
p1
TPaveText * p1
Definition: JDrawModule3D.cc:35
JTOOLS::JQuantile::getSTDev
double getSTDev() const
Get standard deviation.
Definition: JQuantile.hh:310
gmake_property
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
Definition: JProperties.hh:1291
JTOOLS::JQuantile
Quantile calculator.
Definition: JQuantile.hh:88
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:80
JCALIBRATE::JTDC_t
Auxiliary class for TDC constraints.
Definition: JTDC_t.hh:34
JEEP::notice_t
notice
Definition: JMessage.hh:32
JCALIBRATE::JFitK40_t< TF2 >
Template specialisation of two-fold coincidence rate due to K40 and other radioactive decays...
Definition: JFitK40.hh:585
JEEP::JProperties
Utility class to parse parameter values.
Definition: JProperties.hh:496
parameters
*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
Definition: diff-Tuna.sh:38
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:151
h1
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:69
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
JCALIBRATE::JFitK40_t< TF1 >
Template specialisation of two-fold coincidence rate due to K40 and other radioactive decays...
Definition: JFitK40.hh:715
Y
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
Definition: detector-XY:fit1d.sh:32
FIXED
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:445
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JDetector.hh
Data structure for detector geometry and calibration.
JTOOLS::JQuantile::getMean
double getMean() const
Get mean value.
Definition: JQuantile.hh:282
JProperties.hh
Utility class to parse parameter values.
JCALIBRATE::_1F
static const char *const _1F
Name extension for 1D time offset fit.
Definition: JCalibrateK40.hh:39
JCALIBRATE::_2F
static const char *const _2F
Name extension for 2F rate fitted.
Definition: JCalibrateK40.hh:37
JPrint.hh
I/O formatting auxiliaries.
JAANET::detector
Detector file.
Definition: JHead.hh:196
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
JTOOLS::result
return result
Definition: JPolint.hh:727
JDETECTOR::is_valid
bool is_valid(const json &js)
Check validity of JSon data.
Definition: JDetectorCalibration.hh:102
JDETECTOR::store
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
Definition: JDetectorToolkit.hh:516
JMeta.hh
ROOT I/O of application specific meta data.
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JTOOLS::JQuantile::put
void put(const double x, const double w=1.0)
Put value.
Definition: JQuantile.hh:173
JConstants.hh
Physics constants.
JCALIBRATE::_2R
static const char *const _2R
Name extension for 2D rate measured.
Definition: JCalibrateK40.hh:36
JDETECTOR::JPMTParametersMap
Auxiliary class for map of PMT parameters.
Definition: JPMTParametersMap.hh:89
debug
int debug
debug level
Definition: JSirene.cc:63
JMessage.hh
General purpose messaging.
JCALIBRATE::JPMTParameters_t
Fit parameters for single PMT.
Definition: JFitK40.hh:46
range
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi mv $WORKDIR/fit.root $MODULE_ROOT typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
Definition: module-Z:fit.sh:84
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
p2
p2
Definition: module-Z:fit.sh:48
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:440
JRange.hh
Auxiliary class to define a range between two values.
JParser.hh
Utility class to parse command line options.
JROOT::fixParameter
bool fixParameter(TF1 &f1, const JFitParameter_t &parameter)
Fix fit parameter.
Definition: JRootToolkit.hh:491
values
possible values
Definition: detector-XY:run.sh:41
detector
do set_variable DETECTOR_TXT $WORKDIR detector
Definition: detector-XY:fit1d.sh:37
JDAQ.hh
KM3NeT DAQ constants, bit handling, etc.
KM3NETDAQ::NUMBER_OF_PMTS
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
in
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:36
JCALIBRATE::JFitK40Parameters
Fit parameters for two-fold coincidence rate due to K40.
Definition: JFitK40.hh:66
p3
p3
Definition: module-Z:fit.sh:48
JCALIBRATE::_1D
static const char *const _1D
Name extension for 1D time offset.
Definition: JCalibrateK40.hh:38
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
main
int main(int argc, char *argv[])
Definition: Main.cpp:15
Generated by   doxygen 1.8.5