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 #include "JDetector/JPMTParametersToolkit.hh"
23 
24 #include "JROOT/JRootToolkit.hh"
25 #include "JTools/JRange.hh"
26 #include "JTools/JQuantile.hh"
27 #include "JSupport/JMeta.hh"
28 #include "JSupport/JSingleFileScanner.hh"
29 
30 #include "JCalibrate/JCalibrateK40.hh"
31 #include "JCalibrate/JFitK40.hh"
32 #include "JCalibrate/JTDC_t.hh"
33 
34 #include "Jeep/JProperties.hh"
35 #include "Jeep/JPrint.hh"
36 #include "Jeep/JParser.hh"
37 #include "Jeep/JMessage.hh"
38 
39 namespace {
40 
41  double STDEV = 2.0; //!< number of standard deviations
42  double QE_MIN_VALUE = 0.01; //!< minimal QE value (below, set to zero)
43  double QE_MAX_ERROR = 0.5; //!< maximal QE error (below, set to zero)
44  double MINIMAL_RATE_HZ = 1.0e-2; //!< minimal coincidence rate (bin-by-bin)
45 
46  /**
47  * Check validity of QE measurement.
48  *
49  * The QE measurement is valid when:
50  * - fitted value is at least STDEV standard deviations larger than QE_MIN_VALUE; and
51  * - fitted error is smaller than QE_MAX_ERROR.
52  *
53  * \param qe_value QE value
54  * \param qe_error QE error
55  * \return true if valid; else false
56  */
57  inline bool is_valid(const double qe_value, const double qe_error)
58  {
59  return (qe_value > QE_MIN_VALUE + STDEV * qe_error &&
60  qe_error < QE_MAX_ERROR);
61  }
62 }
63 
64 
65 /**
66  * \file
67  *
68  * Auxiliary program to fit PMT parameters from JMergeCalibrateK40.cc output.\n
69  * To force drawing of fit functions bin-by-bin compatible to the corresponding histogram data, \n
70  * the function values will be stored as histograms, rather than as associated function objects. \n
71  * See JCalibrateK40.hh for the corresponding extensions of histogram names.
72  * \author mdejong
73  */
74 int main(int argc, char **argv)
75 {
76  using namespace std;
77  using namespace JPP;
78  using namespace KM3NETDAQ;
79 
80  typedef JRange<double> JRange_t;
81 
82  const double epsilon = 1.0e-10; // minimal error for ROOT plotting
83  JFitK40Parameters fitk40; // setting of internal fit parameters
84 
85  string inputFile;
86  string outputFile;
87  string detectorFile;
88  string pmtFile;
89  JTDC_t TDC;
90  bool reverse;
91  bool overwriteDetector;
92  bool writeFits;
93  bool fitAngDist;
94  bool fitModel;
95  JRange_t X;
96  JRange_t Y;
97  string option;
98  int debug;
99 
100  try {
101 
102  JProperties properties;
103 
104  properties.insert(gmake_property(STDEV));
105  properties.insert(gmake_property(QE_MIN_VALUE));
106  properties.insert(gmake_property(QE_MAX_ERROR));
107  properties.insert(gmake_property(MINIMAL_RATE_HZ));
108  properties.insert(gmake_property(fitk40.Rate_Hz));
109  properties.insert(gmake_property(fitk40.p1));
110  properties.insert(gmake_property(fitk40.p2));
111  properties.insert(gmake_property(fitk40.p3));
112  properties.insert(gmake_property(fitk40.p4));
113 
114  JParser<> zap("Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
115 
116  zap['@'] = make_field(properties) = JPARSER::initialised();
117  zap['f'] = make_field(inputFile, "input file (output from JMergeCalibrateK40).");
118  zap['o'] = make_field(outputFile, "output file.") = "fit.root";
119  zap['a'] = make_field(detectorFile, "detector file.");
120  zap['P'] = make_field(pmtFile, "specify PMT file name that can be given to JTriggerEfficiency.") = "";
121  zap['!'] = make_field(TDC,
122  "Fix time offset(s) of PMT(s) of certain module(s), e.g."
123  "\n-! \"808969848 0 808982077 23\" will fix offset of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
124  "\nSame PMT offset can be fixed for all optical modules, e.g."
125  "\n-! \"-1 0 -1 23\" will fix offsets of PMTs 0 and 23 for all optical modules.")
126  = JPARSER::initialised();
127  zap['r'] = make_field(reverse, "reverse TDC constraints due to option -! <TDC>.");
128  zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with correct time offsets.");
129  zap['w'] = make_field(writeFits, "write fit results.");
130  zap['D'] = make_field(fitAngDist, "fit angular distribution; fix normalisation.");
131  zap['M'] = make_field(fitModel, "fit angular distribution; fix PMT QEs = 1.0.");
132  zap['x'] = make_field(X, "ROOT fit range (PMT pairs).") = JRange_t();
133  zap['y'] = make_field(Y, "ROOT fit range (time residual).") = JRange_t();
134  zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "";
135  zap['d'] = make_field(debug, "debug.") = 1;
136 
137  zap(argc, argv);
138  }
139  catch(const exception &error) {
140  FATAL(error.what() << endl);
141  }
142 
143 
144  ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls(1000000);
145 
146 
147  if (reverse) {
148  TDC.reverse();
149  }
150 
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);
153  }
154 
155  TDC.is_valid(true);
156 
157 
158  //----------------------------------------------------------
159  // load detector file and PMT parameters
160  //----------------------------------------------------------
161 
162  JDetector detector;
163 
164  try {
165  load(detectorFile, detector);
166  }
167  catch(const JException& error) {
168  FATAL(error);
169  }
170 
171 
172  JPMTParametersMap parameters;
173 
174  if (pmtFile != "") {
175  try {
176  parameters.load(pmtFile.c_str());
177  }
178  catch(const JException& error) {
179  //ERROR(error.what());
180  }
181  }
182 
183  parameters.comment.add(JMeta(argc, argv));
184  detector .comment.add(JMeta(argc, argv));
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 (writeFits) {
395 
396  TH1D h2(MAKE_CSTRING(module->getID() << _1F), NULL, h1.GetXaxis()->GetNbins(), h1.GetXaxis()->GetXmin(), h1.GetXaxis()->GetXmax());
397 
398  for (int ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
399 
400  const Double_t x = h2.GetXaxis()->GetBinCenter(ix);
401 
402  h2.SetBinContent(ix, f1.Eval(x));
403  h2.SetBinError (ix, 0.0);
404  }
405 
406  out << h1 << h2;
407  }
408 
409  if (debug >= JEEP::notice_t) {
410 
411  cout << "Fit result chi2 / NDF " << result->Chi2() << ' ' << result->Ndf() << ' ' << (result->IsValid() ? "" : "failed") << endl;
412 
413  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
414 
415  cout << ' '
416  << setw(2) << pmt << ' '
417  << FIXED(7,2) << f1.getT0 (pmt) << ' ';
418  cout << (result->IsParameterFixed(f1.getModelParameter(pmt, &JPMTParameters_t::t0)) ? " *** fixed t0 *** " : "");
419  cout << (f1.is_disabled (pmt) ? " (disabled)" : "");
420  cout << (f1.is_average_t0(pmt) ? " (averaged)" : "");
421  cout << endl;
422  }
423  }
424 
425  fit.setModelParameters(f1.getModelParameters());
426  }
427 
428 
429  TFitResultPtr result = fit(*h2s, option);
430 
431 
432  bool refit = false;
433 
434  try {
435 
436  const JFitK40Parameters values(fit.GetParameters());
437  const JFitK40Parameters errors(fit.GetParErrors());
438 
439  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
440 
441  if (!result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::QE)) && !is_valid(values.getQE(pmt), errors.getQE(pmt))) {
442 
443  WARNING("PMT " << setw(10) << module->getID() << ' ' << setw(2) << pmt << ' '
444  << "measured QE "
445  << FIXED(5,2) << values.getQE(pmt) << " +/- "
446  << FIXED(5,2) << errors.getQE(pmt)
447  << " compatible with zero -> set QE = 0 and t0 = 0; -> refit." << endl);
448 
449 
450  fit.disablePMT(pmt);
451 
452  refit = true;
453  }
454  }
455  }
456  catch(const JException& error) {
457 
458  NOTICE("Skip fit for module " << module->getID() << "; set corresponding QEs to 0." << endl);
459 
460  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
461  parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
462  }
463 
464  continue;
465  }
466 
467 
468  if (refit) {
469 
470  JFitK40_t<TF1> f1(fit);
471 
472  f1(h1, option);
473 
474  fit.setModelParameters(f1.getModelParameters());
475 
476  result = fit(*h2s, option);
477  }
478 
479 
480  const JFitK40Parameters values(fit.GetParameters());
481  const JFitK40Parameters errors(fit.GetParErrors());
482 
483  hc.Fill(TMath::Log10(result->Chi2() / result->Ndf()));
484 
485  TH2D h2f(MAKE_CSTRING(module->getID() << _2F),
486  NULL,
487  h2s->GetXaxis()->GetNbins(), h2s->GetXaxis()->GetXmin (), h2s->GetXaxis()->GetXmax (),
488  h2s->GetYaxis()->GetNbins(), h2s->GetYaxis()->GetXmin (), h2s->GetYaxis()->GetXmax ());
489 
490  for (int ix = 1; ix <= h2f.GetXaxis()->GetNbins(); ++ix) {
491  for (int iy = 1; iy <= h2f.GetYaxis()->GetNbins(); ++iy) {
492 
493  const Double_t x = h2f.GetXaxis()->GetBinCenter(ix);
494  const Double_t y = h2f.GetYaxis()->GetBinCenter(iy);
495 
496  h2f.SetBinContent(ix, iy, fit.Eval(x,y));
497  h2f.SetBinError (ix, iy, 0.0);
498  }
499  }
500 
501 
502  if (debug >= JEEP::notice_t) {
503 
504  cout << "Fit result chi2 / NDF " << result->Chi2() << ' ' << result->Ndf() << ' ' << (result->IsValid() ? "" : "failed") << endl;
505 
506  cout << "Rate_Hz= " << FIXED(8,4) << values.Rate_Hz << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::Rate_Hz)) ? " *** fixed *** " : "") << endl;
507 
508  cout << "p1= " << FIXED(8,4) << values.p1 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p1)) ? " *** fixed *** " : "") << endl;
509  cout << "p2= " << FIXED(8,4) << values.p2 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p2)) ? " *** fixed *** " : "") << endl;
510  cout << "p3= " << FIXED(8,4) << values.p3 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p2)) ? " *** fixed *** " : "") << endl;
511  cout << "p4= " << FIXED(8,4) << values.p4 << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p3)) ? " *** fixed *** " : "") << endl;
512 
513  cout << "Background (correlated): " << FIXED(8,5) << values.bg << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::bg)) ? " *** fixed *** " : "") << endl;
514  cout << "Background (uncorrelated): " << FIXED(8,5) << values.cc << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::cc)) ? " *** fixed *** " : "") << endl;
515 
516  cout << " PMT t0 [ns] TTS [ns] R" << endl;
517 
518  JQuantile Q[3];
519 
520  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
521 
522  cout << ' '
523  << setw(2) << pmt << ' '
524  << FIXED(7,2) << fit.getT0 (pmt) << ' '
525  << FIXED(7,2) << fit.getTTS(pmt) << ' '
526  << FIXED(7,3) << fit.getQE (pmt);
527 
528  cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::QE)) ? " *** fixed QE *** " : "");
529  cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::TTS)) ? " *** fixed TTS *** " : "");
530  cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::t0)) ? " *** fixed t0 *** " : "");
531  cout << (fit.is_disabled (pmt) ? " (disabled)" : "");
532  cout << (fit.is_average_t0(pmt) ? " (averaged)" : "");
533  cout << endl;
534 
535  Q[0].put(fit.getT0 (pmt));
536  Q[1].put(fit.getTTS(pmt));
537  Q[2].put(fit.getQE (pmt));
538  }
539 
540  cout << "<> "
541  << FIXED(7,2) << Q[0].getMean() << ' '
542  << FIXED(7,2) << Q[1].getMean() << ' '
543  << FIXED(7,3) << Q[2].getMean() << endl;
544  cout << "+/- "
545  << FIXED(7,2) << Q[0].getSTDev() << ' '
546  << FIXED(7,2) << Q[1].getSTDev() << ' '
547  << FIXED(7,3) << Q[2].getSTDev() << endl;
548  }
549 
550  h2s->Write();
551  h2f .Write();
552 
553  if (writeFits) {
554 
555  // Histograms with fit results.
556 
557  p1.Fill(values.p1);
558  p2.Fill(values.p2);
559  p3.Fill(values.p3);
560  p4.Fill(values.p4);
561  bg.Fill(values.bg);
562  cc.Fill(values.cc);
563 
564  TH1D h1t(MAKE_CSTRING(module->getID() << ".1t0"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
565  TH1D h1s(MAKE_CSTRING(module->getID() << ".1TTS"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
566  TH1D h1q(MAKE_CSTRING(module->getID() << ".1QE"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
567 
568  h1t.Sumw2();
569  h1s.Sumw2();
570  h1q.Sumw2();
571 
572  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
573 
574  h1t.SetBinContent(pmt + 1, values.getT0 (pmt));
575  h1s.SetBinContent(pmt + 1, values.getTTS(pmt));
576  h1q.SetBinContent(pmt + 1, values.getQE (pmt));
577 
578  h1t.SetBinError (pmt + 1, max(errors.getT0 (pmt), epsilon));
579  h1s.SetBinError (pmt + 1, max(errors.getTTS(pmt), epsilon));
580  h1q.SetBinError (pmt + 1, max(errors.getQE (pmt), epsilon));
581  }
582 
583  out << h1t << h1s << h1q;
584  }
585 
586 
587  // save output
588 
589  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
590 
591  JPMTParameters& data = parameters[JPMTIdentifier(module->getID(), pmt)];
592 
593  // correction to QE
594 
595  const double P = getSurvivalProbability(data);
596 
597  if (P > 0.0)
598  data.QE = values.getQE(pmt) / P;
599  else
600  data.QE = 0.0;
601 
602  if (overwriteDetector) {
603  module->getPMT(pmt).addT0(fit.getT0(pmt));
604  }
605  }
606  }
607 
608 
609  if (writeFits) {
610  out << hc << p1 << p2 << p3 << p4 << bg << cc;
611  }
612 
613  out.Close();
614 
615 
616  for (JSingleFileScanner<JMeta> in(inputFile); in.hasNext(); ) {
617 
618  const JMeta meta(*in.next());
619 
620  parameters.comment.add(meta);
621  detector .comment.add(meta);
622  }
623 
624  if (overwriteDetector) {
625 
626  NOTICE("Store calibration data on file " << detectorFile << endl);
627 
628  store(detectorFile, detector);
629  }
630 
631  if (pmtFile != "") {
632  parameters.store(pmtFile.c_str());
633  }
634 
635  return 0;
636 }
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
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: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
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:103
JDETECTOR::store
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
Definition: JDetectorToolkit.hh:521
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:445
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
JROOT::fixParameter
bool fixParameter(TF1 &f1, const JFitParameter_t &parameter)
Fix fit parameter.
Definition: JRootToolkit.hh:538
JSUPPORT::JSingleFileScanner
Object reading from a list of files.
Definition: JSingleFileScanner.hh:75
values
possible values
Definition: detector-XY:run.sh:41
JDETECTOR::JPMTParameters
Data structure for PMT parameters.
Definition: JPMTParameters.hh:29
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
JDETECTOR::JPMTParameters::QE
double QE
relative quantum efficiency
Definition: JPMTParameters.hh:233
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
P
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:60
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Generated by   doxygen 1.8.5