Jpp  15.0.1-rc.1-highQE
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 
8 #include "km3net-dataformat/online/JDAQ.hh"
9 #include "km3net-dataformat/online/JDAQHeader.hh"
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"
22 #include "JDetector/JDetectorToolkit.hh"
23 #include "JDetector/JPMTParametersMap.hh"
24 #include "JDetector/JPMTParametersToolkit.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"
31 #include "JSupport/JSingleFileScanner.hh"
32 
33 #include "JCalibrate/JCalibrateK40.hh"
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 (below, 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.")
129  = JPARSER::initialised();
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 
165  JDetector detector;
166 
167  try {
168  load(detectorFile, detector);
169  }
170  catch(const JException& error) {
171  FATAL(error);
172  }
173 
174 
175  JPMTParametersMap parameters;
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  parameters.comment.add(JMeta(argc, argv));
187  detector .comment.add(JMeta(argc, argv));
188 
189 
190  if (option.find('O') == string::npos) { option += '0'; }
191  if (option.find('S') == string::npos) { option += 'S'; }
192  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
193 
194 
195  TFile* in = TFile::Open(inputFile.c_str(), "exist");
196 
197  if (in == NULL || !in->IsOpen()) {
198  FATAL("File: " << inputFile << " not opened." << endl);
199  }
200 
201 
202  TDirectory::AddDirectory(NULL);
203 
204  TFile out(outputFile.c_str(), "recreate");
205 
206  // Histograms with global fit results.
207 
208  TH1D hc("chi2", NULL, 500, 0.0, 5.0);
209 
210  TH1D p1("p1", NULL, 501, -5.0, +5.0);
211  TH1D p2("p2", NULL, 501, -5.0, +5.0);
212  TH1D p3("p3", NULL, 501, -5.0, +5.0);
213  TH1D p4("p4", NULL, 501, -5.0, +5.0);
214  TH1D bg("bg", NULL, 501, -0.1, +0.1);
215  TH1D cc("cc", NULL, 501, -0.1, +0.1);
216 
217 
218  //----------------------------------------------------------
219  // loop over modules in detector file to fit histogram
220  //----------------------------------------------------------
221 
222  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
223 
224  NOTICE("Module " << setw(10) << module->getID() << ' ' << module->getLocation() << endl);
225 
226  if (module->empty()) {
227  continue;
228  }
229 
230  //----------------------------------------------------------
231  // get histogram for this module
232  //----------------------------------------------------------
233 
234  TH2D* h2s = (TH2D*) in->Get(MAKE_CSTRING(module->getID() << _2R));
235 
236  if (h2s == NULL || h2s->GetEntries() == 0) {
237 
238  NOTICE("No data for module " << module->getID() << "; set corresponding QEs to 0." << endl);
239 
240  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
241  parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
242  }
243 
244  continue;
245  }
246 
247  //----------------------------------------------------------
248  // get constraints
249  //----------------------------------------------------------
250 
251  const JTDC_t::range_type range = TDC.equal_range(module->getID());
252 
253  //----------------------------------------------------------
254  // set-up fitting procedure
255  //----------------------------------------------------------
256 
257  JFitK40 fit(*module,
258  h2s->GetXaxis()->GetXmin(),
259  h2s->GetXaxis()->GetXmax(),
260  h2s->GetYaxis()->GetXmin(),
261  h2s->GetYaxis()->GetXmax(),
262  range.first == range.second);
263 
264  fit.setModelParameters(fitk40.getModelParameters());
265 
266  for (JTDC_t::const_iterator i = range.first; i != range.second; ++i) {
267  fixParameter(fit, fit.getModelParameter(i->second, &JPMTParameters_t::t0));
268  }
269 
270  if (fitModel) {
271  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
272  fixParameter(fit, fit.getModelParameter(pmt, &JPMTParameters_t::QE));
273  }
274 
275  } else {
276 
277  fixParameter(fit, fit.getModelParameter(&JFitK40::Rate_Hz));
278 
279  if (!fitAngDist) {
280  fixParameter(fit, fit.getModelParameter(&JFitK40::p1));
281  fixParameter(fit, fit.getModelParameter(&JFitK40::p2));
282  fixParameter(fit, fit.getModelParameter(&JFitK40::p3));
283  fixParameter(fit, fit.getModelParameter(&JFitK40::p4));
284  }
285  }
286 
287  //----------------------------------------------------------
288  // check validity of data
289  //----------------------------------------------------------
290 
291  vector<int> buffer(NUMBER_OF_PMTS, 1);
292 
293  for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
294 
295  Double_t value = 0.0;
296  Double_t error = 0.0;
297 
298  for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
299  value += h2s->GetBinContent(ix,iy);
300  error += h2s->GetBinError(ix,iy) * h2s->GetBinError(ix,iy);
301  }
302 
303  error = sqrt(error);
304 
305  if (value < MINIMAL_RATE_HZ + STDEV*error) {
306 
307  const JFitK40::pair_type pair = fit.getPair(ix-1);
308 
309  buffer[pair.first] += 1;
310  buffer[pair.second] += 1;
311  }
312  }
313 
314  try {
315 
316  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
317 
318  if (buffer[pmt] == NUMBER_OF_PMTS) {
319 
320  WARNING("PMT " << setw(10) << module->getID() << ' ' << setw(2) << pmt << " too low a rate; switch off." << endl);
321 
322  fit.disablePMT(pmt);
323  }
324  }
325  }
326  catch(const JException& error) {
327 
328  NOTICE("Skip fit for module " << module->getID() << "; set corresponding QEs to 0." << endl);
329 
330  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
331  parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
332  }
333 
334  continue;
335  }
336 
337  // constrain fit to specified range(s)
338 
339  if (X != JRange_t()) {
340  h2s->GetXaxis()->SetRangeUser(max(X.getLowerLimit(), h2s->GetXaxis()->GetXmin()),
341  min(X.getUpperLimit(), h2s->GetXaxis()->GetXmax()));
342  }
343  if (Y != JRange_t()) {
344  h2s->GetYaxis()->SetRangeUser(max(Y.getLowerLimit(), h2s->GetYaxis()->GetXmin()),
345  min(Y.getUpperLimit(), h2s->GetYaxis()->GetXmax()));
346  }
347 
348  // 1D-histogram with peak positions
349 
350  TH1D h1(MAKE_CSTRING(module->getID() << _1D), NULL, h2s->GetXaxis()->GetNbins(), h2s->GetXaxis()->GetXmin(), h2s->GetXaxis()->GetXmax());
351 
352  const Double_t sigma[] = {
353 
354  (min(Y.getUpperLimit(), h2s->GetYaxis()->GetXmax()) -
355  max(Y.getLowerLimit(), h2s->GetYaxis()->GetXmin())) / sqrt(12.0),
356 
357  sqrt(JPMTParameters_t().TTS * JPMTParameters_t().TTS * 2.0 + fit.getSigmaK40() * fit.getSigmaK40())
358  };
359 
360 
361  for (Int_t ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
362 
363  Double_t y0 = 0.0;
364  Double_t zmax = 0.0;
365  Double_t value = 0.0;
366  Double_t error = 0.0;
367 
368  for (Int_t iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
369 
370  const Double_t y = h2s->GetYaxis()->GetBinCenter(iy);
371  const Double_t z = h2s->GetBinContent(ix, iy);
372  const Double_t w = h2s->GetBinError (ix, iy);
373 
374  if (Y(y)) {
375 
376  if (z > zmax) {
377  y0 = y;
378  zmax = z;
379  }
380 
381  value += z;
382  error += w*w;
383  }
384  }
385 
386  error = sqrt(error);
387 
388  h1.SetBinContent(ix, y0);
389  h1.SetBinError (ix, value > MINIMAL_RATE_HZ + STDEV * error ? sigma[1] * error / value : sigma[0]);
390  }
391 
392  {
393  JFitK40_t<TF1> f1(fit);
394 
395  TFitResultPtr result = f1(h1, option.c_str());
396 
397  if (result.Get() == NULL) {
398  FATAL("Invalid TFitResultPtr " << h1.GetName() << endl);
399  }
400 
401  if (writeFits) {
402 
403  TH1D h2(MAKE_CSTRING(module->getID() << _1F), NULL, h1.GetXaxis()->GetNbins(), h1.GetXaxis()->GetXmin(), h1.GetXaxis()->GetXmax());
404 
405  for (int ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
406 
407  const Double_t x = h2.GetXaxis()->GetBinCenter(ix);
408 
409  h2.SetBinContent(ix, f1.Eval(x));
410  h2.SetBinError (ix, 0.0);
411  }
412 
413  out << h1 << h2;
414  }
415 
416  if (debug >= JEEP::notice_t) {
417 
418  cout << "Fit result chi2 / NDF " << result->Chi2() << ' ' << result->Ndf() << ' ' << (result->IsValid() ? "" : "failed") << endl;
419 
420  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
421 
422  cout << ' '
423  << setw(2) << pmt << ' '
424  << FIXED(7,2) << f1.getT0 (pmt) << ' ';
425  cout << (result->IsParameterFixed(f1.getModelParameter(pmt, &JPMTParameters_t::t0)) ? " *** fixed t0 *** " : "");
426  cout << (f1.is_disabled (pmt) ? " (disabled)" : "");
427  cout << (f1.is_average_t0(pmt) ? " (averaged)" : "");
428  cout << endl;
429  }
430  }
431 
432  fit.setModelParameters(f1.getModelParameters());
433  }
434 
435 
436  TFitResultPtr result = fit(*h2s, option);
437 
438 
439  bool refit = false;
440 
441  try {
442 
443  const JFitK40Parameters values(fit.GetParameters());
444  const JFitK40Parameters errors(fit.GetParErrors());
445 
446  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
447 
448  if (!result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::QE)) && !is_valid(values.getQE(pmt), errors.getQE(pmt))) {
449 
450  WARNING("PMT " << setw(10) << module->getID() << ' ' << setw(2) << pmt << ' '
451  << "measured QE "
452  << FIXED(5,2) << values.getQE(pmt) << " +/- "
453  << FIXED(5,2) << errors.getQE(pmt)
454  << " compatible with zero -> set QE = 0 and t0 = 0; -> refit." << endl);
455 
456 
457  fit.disablePMT(pmt);
458 
459  refit = true;
460  }
461  }
462  }
463  catch(const JException& error) {
464 
465  NOTICE("Skip fit for module " << module->getID() << "; set corresponding QEs to 0." << endl);
466 
467  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
468  parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
469  }
470 
471  continue;
472  }
473 
474 
475  if (refit) {
476 
477  JFitK40_t<TF1> f1(fit);
478 
479  f1(h1, option);
480 
481  fit.setModelParameters(f1.getModelParameters());
482 
483  result = fit(*h2s, option);
484  }
485 
486 
487  const JFitK40Parameters values(fit.GetParameters());
488  const JFitK40Parameters errors(fit.GetParErrors());
489 
490  hc.Fill(TMath::Log10(result->Chi2() / result->Ndf()));
491 
492  TH2D h2f(MAKE_CSTRING(module->getID() << _2F),
493  NULL,
494  h2s->GetXaxis()->GetNbins(), h2s->GetXaxis()->GetXmin (), h2s->GetXaxis()->GetXmax (),
495  h2s->GetYaxis()->GetNbins(), h2s->GetYaxis()->GetXmin (), h2s->GetYaxis()->GetXmax ());
496 
497  for (int ix = 1; ix <= h2f.GetXaxis()->GetNbins(); ++ix) {
498  for (int iy = 1; iy <= h2f.GetYaxis()->GetNbins(); ++iy) {
499 
500  const Double_t x = h2f.GetXaxis()->GetBinCenter(ix);
501  const Double_t y = h2f.GetYaxis()->GetBinCenter(iy);
502 
503  h2f.SetBinContent(ix, iy, fit.Eval(x,y));
504  h2f.SetBinError (ix, iy, 0.0);
505  }
506  }
507 
508 
509  if (debug >= JEEP::notice_t) {
510 
511  cout << "Fit result chi2 / NDF " << result->Chi2() << ' ' << result->Ndf() << ' ' << (result->IsValid() ? "" : "failed") << endl;
512 
513  cout << "Rate_Hz= " << FIXED(8,4) << values.Rate_Hz << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::Rate_Hz)) ? " *** fixed *** " : "") << endl;
514 
515  cout << "p1= " << FIXED(8,4) << values.p1 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p1)) ? " *** fixed *** " : "") << endl;
516  cout << "p2= " << FIXED(8,4) << values.p2 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p2)) ? " *** fixed *** " : "") << endl;
517  cout << "p3= " << FIXED(8,4) << values.p3 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p2)) ? " *** fixed *** " : "") << endl;
518  cout << "p4= " << FIXED(8,4) << values.p4 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p3)) ? " *** fixed *** " : "") << endl;
519 
520  cout << "Background (correlated): " << FIXED(8,5) << values.bg << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::bg)) ? " *** fixed *** " : "") << endl;
521  cout << "Background (uncorrelated): " << FIXED(8,5) << values.cc << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::cc)) ? " *** fixed *** " : "") << endl;
522 
523  cout << " PMT t0 [ns] TTS [ns] R" << endl;
524 
525  JQuantile Q[3];
526 
527  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
528 
529  cout << ' '
530  << setw(2) << pmt << ' '
531  << FIXED(7,2) << fit.getT0 (pmt) << ' '
532  << FIXED(7,2) << fit.getTTS(pmt) << ' '
533  << FIXED(7,3) << fit.getQE (pmt);
534 
535  cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::QE)) ? " *** fixed QE *** " : "");
536  cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::TTS)) ? " *** fixed TTS *** " : "");
537  cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::t0)) ? " *** fixed t0 *** " : "");
538  cout << (fit.is_disabled (pmt) ? " (disabled)" : "");
539  cout << (fit.is_average_t0(pmt) ? " (averaged)" : "");
540  cout << endl;
541 
542  Q[0].put(fit.getT0 (pmt));
543  Q[1].put(fit.getTTS(pmt));
544  Q[2].put(fit.getQE (pmt));
545  }
546 
547  cout << "<> "
548  << FIXED(7,2) << Q[0].getMean() << ' '
549  << FIXED(7,2) << Q[1].getMean() << ' '
550  << FIXED(7,3) << Q[2].getMean() << endl;
551  cout << "+/- "
552  << FIXED(7,2) << Q[0].getSTDev() << ' '
553  << FIXED(7,2) << Q[1].getSTDev() << ' '
554  << FIXED(7,3) << Q[2].getSTDev() << endl;
555  }
556 
557  h2s->Write();
558  h2f .Write();
559 
560  if (writeFits) {
561 
562  // Histograms with fit results.
563 
564  p1.Fill(values.p1);
565  p2.Fill(values.p2);
566  p3.Fill(values.p3);
567  p4.Fill(values.p4);
568  bg.Fill(values.bg);
569  cc.Fill(values.cc);
570 
571  TH1D h1t(MAKE_CSTRING(module->getID() << ".1t0"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
572  TH1D h1s(MAKE_CSTRING(module->getID() << ".1TTS"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
573  TH1D h1q(MAKE_CSTRING(module->getID() << ".1QE"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
574 
575  h1t.Sumw2();
576  h1s.Sumw2();
577  h1q.Sumw2();
578 
579  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
580 
581  h1t.SetBinContent(pmt + 1, values.getT0 (pmt));
582  h1s.SetBinContent(pmt + 1, values.getTTS(pmt));
583  h1q.SetBinContent(pmt + 1, values.getQE (pmt));
584 
585  h1t.SetBinError (pmt + 1, max(errors.getT0 (pmt), epsilon));
586  h1s.SetBinError (pmt + 1, max(errors.getTTS(pmt), epsilon));
587  h1q.SetBinError (pmt + 1, max(errors.getQE (pmt), epsilon));
588  }
589 
590  out << h1t << h1s << h1q;
591  }
592 
593 
594  // save output
595 
596  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
597 
598  JPMTParameters& data = parameters[JPMTIdentifier(module->getID(), pmt)];
599 
600  // correction to QE
601 
602  const double P = getSurvivalProbability(data);
603 
604  if (P > 0.0)
605  data.QE = values.getQE(pmt) / P;
606  else
607  data.QE = 0.0;
608 
609  if (overwriteDetector) {
610  module->getPMT(pmt).addT0(fit.getT0(pmt));
611  }
612  }
613  }
614 
615  for (JRootFileReader<JDAQHeader> in(inputFile.c_str()); in.hasNext(); ) {
616  putObject(out, *in.next());
617  }
618 
619  if (writeFits) {
620  out << hc << p1 << p2 << p3 << p4 << bg << cc;
621  }
622 
623  out.Close();
624 
625 
626  for (JSingleFileScanner<JMeta> in(inputFile); in.hasNext(); ) {
627 
628  const JMeta meta(*in.next());
629 
630  parameters.comment.add(meta);
631  detector .comment.add(meta);
632  }
633 
634  if (overwriteDetector) {
635 
636  NOTICE("Store calibration data on file " << detectorFile << endl);
637 
638  store(detectorFile, detector);
639  }
640 
641  if (pmtFile != "") {
642  parameters.store(pmtFile.c_str());
643  }
644 
645  return 0;
646 }
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
Q
Q(UTCMax_s-UTCMin_s)-livetime_s
JEEP::debug_t
debug
Definition: JMessage.hh:29
main
int main(int argc, char *argv[])
Definition: Main.cc:15
p1
TPaveText * p1
Definition: JDrawModule3D.cc:35
JTOOLS::JQuantile::getSTDev
double getSTDev() const
Get standard deviation.
Definition: JQuantile.hh:266
gmake_property
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
Definition: JProperties.hh:1291
JTOOLS::JQuantile
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:89
JCALIBRATE::JTDC_t
Auxiliary class for TDC constraints.
Definition: JTDC_t.hh:37
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:586
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:71
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:716
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:33
FIXED
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
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:252
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
JSingleFileScanner.hh
Scanning of objects from a single file according a format that follows from the extension of each fil...
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:123
JDETECTOR::store
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
Definition: JDetectorToolkit.hh:568
JMeta.hh
ROOT I/O of application specific meta data.
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
X
then break fi done getCenter read X Y Z let X
Definition: acoustics-fit.sh:387
JTOOLS::JQuantile::put
void put(const double x, const double w=1.0)
Put value.
Definition: JQuantile.hh:133
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
JROOT::putObject
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
Definition: JRootFileWriter.hh:41
JMessage.hh
General purpose messaging.
JCALIBRATE::JPMTParameters_t
Fit parameters for single PMT.
Definition: JFitK40.hh:46
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
range
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
Definition: module-Z:fit.sh:118
p2
p2
Definition: module-Z:fit.sh:74
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:475
JRange.hh
Auxiliary class to define a range between two values.
JParser.hh
Utility class to parse command line options.
JDETECTOR::getSurvivalProbability
double getSurvivalProbability(const JPMTParameters &parameters)
Get model dependent probability that a one photo-electron hit survives the simulation of the PMT assu...
Definition: JPMTParametersToolkit.hh:32
DEBUG
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
Definition: JCalibratePMT.sh:21
JROOT::fixParameter
bool fixParameter(TF1 &f1, const JFitParameter_t &parameter)
Fix fit parameter.
Definition: JRootToolkit.hh:532
TDC
set_variable TDC
Definition: JPrintTDC.sh:20
JSUPPORT::JSingleFileScanner
Object reading from a list of files.
Definition: JSingleFileScanner.hh:75
values
possible values
Definition: detector-XY:run.sh:43
JDETECTOR::JPMTParameters
Data structure for PMT parameters.
Definition: JPMTParameters.hh:29
detector
do set_variable DETECTOR_TXT $WORKDIR detector
Definition: detector-XY:fit1d.sh:38
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 CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:41
JCALIBRATE::JFitK40Parameters
Fit parameters for two-fold coincidence rate due to K40.
Definition: JFitK40.hh:66
JDETECTOR::JPMTParameters::QE
double QE
relative quantum efficiency
Definition: JPMTParameters.hh:233
JROOT::JRootFileReader
ROOT file reader.
Definition: JRootFileReader.hh:54
p3
p3
Definition: module-Z:fit.sh:74
JCALIBRATE::_1D
static const char *const _1D
Name extension for 1D time offset.
Definition: JCalibrateK40.hh:38
P
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:62
Generated by   doxygen 1.8.5