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