Jpp  master_rocky-40-g5f0272dcd
the software that should make you happy
JFitK40.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <cmath>
5 #include <set>
6 #include <vector>
7 #include <algorithm>
8 
9 #include "TROOT.h"
10 #include "TFile.h"
11 #include "TH1D.h"
12 #include "TH2D.h"
13 
16 
17 #include "JTools/JConstants.hh"
18 #include "JDetector/JDetector.hh"
22 
23 #include "JROOT/JRootFileWriter.hh"
24 #include "JROOT/JRootToolkit.hh"
25 
26 #include "JTools/JRange.hh"
27 
29 
30 #include "JSupport/JMeta.hh"
32 
34 #include "JCalibrate/JTDC_t.hh"
35 
36 #include "JFitK40.hh"
37 
38 #include "Jeep/JProperties.hh"
39 #include "Jeep/JPrint.hh"
40 #include "Jeep/JParser.hh"
41 #include "Jeep/JMessage.hh"
42 
43 
44 namespace {
45 
46  double MINIMAL_RATE_HZ = 1.0e-2; //!< minimal coincidence rate [Hz] (below, disable PMT)
47  double STDEV = 2.0; //!< number of standard deviations
48  size_t MAXIMAL_COUNTS = 10; //!< maximal count of too low rates
49  double QE_MIN = 0.1; //!< minimal QE (below, disable PMT)
50 }
51 
52 
53 /**
54  * \file
55  *
56  * Auxiliary program to fit PMT parameters from JMergeCalibrateK40.cc output.\n
57  * Note that the contebts of the 2D histograms are converted to standard deviations.
58  *
59  * \author mdejong
60  */
61 int main(int argc, char **argv)
62 {
63  using namespace std;
64  using namespace JPP;
65  using namespace KM3NETDAQ;
66 
67  typedef JRange<double> JRange_t;
68 
70 
71  string inputFile;
72  string outputFile;
73  string detectorFile;
74  string pmtFile;
75  JTDC_t TDC;
76  bool reverse;
77  bool overwriteDetector;
78  JCounter writeFits;
79  bool fitAngle;
80  bool fitNoise;
81  bool fitModel;
82  bool fixQE;
83  JRange_t X;
84  JRange_t Y;
85  int debug;
86 
87  try {
88 
89  JProperties properties;
90 
91  properties.insert(gmake_property(MINIMAL_RATE_HZ));
92  properties.insert(gmake_property(STDEV));
93  properties.insert(gmake_property(MAXIMAL_COUNTS));
94  properties.insert(gmake_property(QE_MIN));
95  properties.insert(gmake_property(K40.R));
96  properties.insert(gmake_property(K40.p1));
97  properties.insert(gmake_property(K40.p2));
98  properties.insert(gmake_property(K40.p3));
99  properties.insert(gmake_property(K40.p4));
100 
101  JParser<> zap("Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
102 
103  zap['@'] = make_field(properties) = JPARSER::initialised();
104  zap['f'] = make_field(inputFile, "input file (output from JMergeCalibrateK40).");
105  zap['o'] = make_field(outputFile, "output file.") = "fit.root";
106  zap['a'] = make_field(detectorFile, "detector file.");
107  zap['P'] = make_field(pmtFile, "specify PMT file name that can be given to JTriggerEfficiency.") = "";
108  zap['!'] = make_field(TDC,
109  "Fix time offset(s) of PMT(s) of certain module(s), e.g."
110  "\n-! \"808969848 0 808982077 23\" will fix time offsets of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
111  "\nSame PMT offset can be fixed for all optical modules, e.g."
112  "\n-! \"-1 0 -1 23\" will fix time offsets of PMTs 0 and 23 of all optical modules.") = JPARSER::initialised();
113  zap['r'] = make_field(reverse, "reverse TDC constraints due to option -! <TDC>.");
114  zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with fitted time offsets.");
115  zap['w'] = make_field(writeFits, "write fit results to ROOT file; -ww also write fitted TTS to PMT file.");
116  zap['D'] = make_field(fitAngle, "fit angular distribution; fix normalisation.");
117  zap['B'] = make_field(fitNoise, "fit background.");
118  zap['M'] = make_field(fitModel, "fit angular distribution as well as normalisation; fix PMT QEs = 1.0.");
119  zap['Q'] = make_field(fixQE, "fix PMT QEs = 1.0.");
120  zap['x'] = make_field(X, "fit range (PMT pairs).") = JRange_t();
121  zap['y'] = make_field(Y, "fit range (time residual [ns]).") = JRange_t();
122  zap['d'] = make_field(debug, "debug.") = 1;
123 
124  zap(argc, argv);
125  }
126  catch(const exception &error) {
127  FATAL(error.what() << endl);
128  }
129 
130 
131  if ((fitModel ? 1 : 0) +
132  (fitAngle ? 1 : 0) +
133  (fitNoise ? 1 : 0) +
134  (fixQE ? 1 : 0) > 1) {
135  FATAL("Use either option -M, -D, -B or -Q" << endl);
136  }
137 
138  const int option = (fitModel ? FIT_MODEL_t :
140  fitNoise ? FIT_PMTS_AND_BACKGROUND_t :
141  fixQE ? FIT_PMTS_QE_FIXED_t :
142  FIT_PMTS_t);
143 
144  if (reverse) {
145  TDC.reverse();
146  }
147 
148  for (JTDC_t::const_iterator i = TDC.begin(); i != TDC.end(); ++i) {
149  DEBUG("PMT " << setw(10) << i->first << ' ' << setw(2) << i->second << " constrain t0." << endl);
150  }
151 
152  try {
153  TDC.is_valid(true);
154  }
155  catch(const exception &error) {
156  FATAL(error.what() << endl);
157  }
158 
160 
161  try {
162  load(detectorFile, detector);
163  }
164  catch(const JException& error) {
165  FATAL(error);
166  }
167 
168 
169  JPMTParametersMap parameters;
170 
171  if (pmtFile != "") {
172  try {
173  parameters.load(pmtFile.c_str());
174  }
175  catch(const exception& error) {}
176  }
177 
178  parameters.comment.add(JMeta(argc, argv));
179 
180 
181  TFile* in = TFile::Open(inputFile.c_str(), "exist");
182 
183  if (in == NULL || !in->IsOpen()) {
184  FATAL("File: " << inputFile << " not opened." << endl);
185  }
186 
187 
188  TFile out(outputFile.c_str(), "recreate");
189 
190  TH1D h0("chi2", NULL, 500, 0.0, 5.0);
191  TH1D hn("hn", NULL, 501, -0.5, 500.0);
192  TH1D hr("rate", NULL, 500, 0.0, 25.0);
193  TH1D h1("p1", NULL, 500, -5.0, +5.0);
194  TH1D h2("p2", NULL, 500, -5.0, +5.0);
195  TH1D h3("p3", NULL, 500, -5.0, +5.0);
196  TH1D h4("p4", NULL, 500, -5.0, +5.0);
197  TH1D hc("cc", NULL, 500, -0.1, +0.1);
198 
199  const floor_range range = getRangeOfFloors(detector);
200  const JStringRouter string(detector);
201 
202  TH2D H2("detector", NULL,
203  string.size() + 0, -0.5, string.size() - 0.5,
204  range.getUpperLimit(), 1 - 0.5, range.getUpperLimit() + 0.5);
205 
206  for (Int_t i = 1; i <= H2.GetXaxis()->GetNbins(); ++i) {
207  H2.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1)));
208  }
209  for (Int_t i = 1; i <= H2.GetYaxis()->GetNbins(); ++i) {
210  H2.GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i));
211  }
212 
213 
214  JFit fit(0, debug);
215 
216  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
217 
218  if (module->getFloor() == 0) {
219  continue;
220  }
221 
222  const JTDC_t::range_type range = TDC.equal_range(module->getID());
223 
224  NOTICE("Module " << setw(10) << module->getID() << ' ' << getLabel(module->getLocation()) << " !" << distance(range.first, range.second) << endl);
225 
226  TH2D* h2d = (TH2D*) in->Get(MAKE_CSTRING(module->getID() << _2R));
227 
228  if (h2d == NULL || h2d->GetEntries() == 0) {
229 
230  NOTICE("No data for module " << module->getID() << " -> set QEs to 0." << endl);
231 
232  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
233  parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
234  }
235 
236  continue;
237  }
238 
239  JModel model(*module, K40, range, option);
240 
241  data_type data; // input data
242 
243  vector<size_t> count[2] = { vector<size_t>(NUMBER_OF_PMTS, 0), // counter (exclude self)
244  vector<size_t>(NUMBER_OF_PMTS, 1) }; // counter (include self)
245 
246  for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
247 
248  const pair_type pair = model.getPair(ix - 1);
249 
250  auto& buffer = data[pair]; // storage
251 
252  double V = 0.0; // integrated value
253  double W = 0.0; // integrated error
254 
255  for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
256 
257  const double x = h2d->GetXaxis()->GetBinCenter(ix);
258  const double y = h2d->GetYaxis()->GetBinCenter(iy);
259 
260  if (X(x) && Y(y)) {
261 
262  double value = h2d->GetBinContent(ix,iy);
263  double error = h2d->GetBinError (ix,iy);
264 
265  buffer.push_back(rate_type(y, value, error));
266 
267  double width = h2d->GetYaxis()->GetBinWidth(iy);
268 
269  value *= width;
270  error *= width;
271 
272  V += value;
273  W += error * error;
274  }
275  }
276 
277  W = sqrt(W);
278 
279  if (V <= 0.0 - STDEV*W) {
280  count[0][pair.first] += 1;
281  count[0][pair.second] += 1;
282  }
283 
284  if (V <= MINIMAL_RATE_HZ + STDEV*W) {
285  count[1][pair.first] += 1;
286  count[1][pair.second] += 1;
287  }
288  }
289 
290  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
291 
292  if (count[0][pmt] >= MAXIMAL_COUNTS) { // too many paired rates negative
293 
294  WARNING("PMT " << setw(10) << module->getID() << '.' << FILL(2,'0') << pmt << FILL() << " some rates negative -> fit background" << endl);
295 
296  if (fit.value.parameters[pmt].status) {
297  model.parameters[pmt].bg.set();
298  }
299  }
300 
301  if (count[1][pmt] == NUMBER_OF_PMTS) { // all paired rates too low
302 
303  WARNING("PMT " << setw(10) << module->getID() << '.' << FILL(2,'0') << pmt << FILL() << " all rates to low -> disable" << endl);
304 
305  model.parameters[pmt].disable();
306  }
307  }
308 
309  DEBUG("Start value:" << endl << model << endl);
310 
311  try {
312 
313  fit.value = model; // start value
314 
315  auto result = fit(data);
316 
317  bool refit = false;
318 
319  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
320 
321  if (fit.value.parameters[pmt].QE() < QE_MIN && fit.value.parameters[pmt].status) {
322 
323  WARNING("PMT " << setw(10) << module->getID() << '.' << FILL(2,'0') << pmt << FILL() << ' ' << FIXED(5,3) << fit.value.parameters[pmt].QE() << " < " << FIXED(5,3) << QE_MIN << " -> disable" << (!refit ? " and refit" : "") << endl);
324 
325  fit.value.parameters[pmt].disable();
326 
327  refit = true;
328  }
329  }
330 
331  if (refit) {
332  refit = false;
333  result = fit(data);
334  }
335 
336  NOTICE("Fit result " << setw(10) << module->getID() << " chi2 / NDF " << FIXED(10,2) << result.chi2 << " / " << setw(5) << result.ndf << ' ' << setw(5) << fit.numberOfIterations << endl);
337 
338  DEBUG(fit.value);
339 
340  if (writeFits) {
341 
342  h0.Fill(result.chi2 / result.ndf);
343  hn.Fill(fit.numberOfIterations);
344  hr.Fill(fit.value.R );
345  h1.Fill(fit.value.p1);
346  h2.Fill(fit.value.p2);
347  h3.Fill(fit.value.p3);
348  h4.Fill(fit.value.p4);
349  hc.Fill(fit.value.cc);
350 
351  TH1D h1t(MAKE_CSTRING(module->getID() << ".1t0"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
352  TH1D h1s(MAKE_CSTRING(module->getID() << ".1TTS"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
353  TH1D h1q(MAKE_CSTRING(module->getID() << ".1QE"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
354 
355  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
356  h1t.SetBinContent(pmt + 1, fit.value.parameters[pmt].t0 ());
357  h1s.SetBinContent(pmt + 1, fit.value.parameters[pmt].TTS());
358  h1q.SetBinContent(pmt + 1, fit.value.parameters[pmt].QE ());
359  }
360 
361  out << h1t << h1s << h1q;
362 
363  for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
364 
365  const pair_type pair = fit.value.getPair(ix - 1);
366 
367  for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
368 
369  const double dt_ns = h2d->GetYaxis()->GetBinCenter(iy);
370 
371  h2d->SetBinContent(ix, iy, fit.value.getValue(pair, dt_ns));
372  h2d->SetBinError (ix, iy, 0.0);
373  }
374  }
375 
376  h2d->SetName(MAKE_CSTRING(module->getID() << _2F));
377  h2d->Write();
378 
379  H2.Fill((double) string.getIndex(module->getString()), (double) module->getFloor(), result.chi2 / result.ndf);
380  }
381 
382  const double t0 = (fit.value.hasFixedTimeOffset() ? fit.value.getFixedTimeOffset() : 0.0);
383 
384  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
385 
386  JPMTParameters& data = parameters[JPMTIdentifier(module->getID(), pmt)];
387 
388  const double P = getSurvivalProbability(data);
389 
390  if (P > 0.0)
391  data.QE = fit.value.parameters[pmt].QE / P;
392  else
393  data.QE = 0.0;
394 
395  if (writeFits > 1) {
396  data.TTS_ns = fit.value.parameters[pmt].TTS();
397  }
398 
399  module->getPMT(pmt).addT0(fit.value.parameters[pmt].t0.get() - t0);
400  }
401  }
402  catch(const exception& error) {
403 
404  ERROR("Module " << setw(10) << module->getID() << ' ' << error.what() << " -> set QEs to 0." << endl);
405 
406  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
407  parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
408  }
409  }
410  }
411 
412 
413  vector<JMeta> meta(1, JMeta(argc, argv));
414 
415  {
416  JSingleFileScanner<JMeta> reader(inputFile);
417  JSTDObjectWriter <JMeta> writer(meta);
418 
419  writer << reader;
420  }
421 
422  for (vector<JMeta>::const_reverse_iterator i = meta.rbegin(); i != meta.rend(); ++i) {
423  parameters.comment.add(*i);
424  detector .comment.add(*i);
425  }
426 
427  if (overwriteDetector) {
428 
429  NOTICE("Store calibration data on file " << detectorFile << endl);
430 
431  store(detectorFile, detector);
432  }
433 
434  if (pmtFile != "") {
435  parameters.store(pmtFile.c_str());
436  }
437 
438 
439  for (vector<JMeta>::const_iterator i = meta.begin(); i != meta.end(); ++i) {
440  putObject(out,*i);
441  }
442 
443  for (JRootFileReader<JDAQHeader> in(inputFile.c_str()); in.hasNext(); ) {
444  putObject(out, *in.next());
445  }
446 
447  if (writeFits) {
448  out << h0 << hn << hr << h1 << h2 << h3 << h4 << hc << H2;
449  }
450 
451  out.Close();
452 
453  return 0;
454 }
string outputFile
KM3NeT DAQ constants, bit handling, etc.
Data structure for detector geometry and calibration.
int main(int argc, char **argv)
Definition: JFitK40.cc:61
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
#define ERROR(A)
Definition: JMessage.hh:66
#define NOTICE(A)
Definition: JMessage.hh:64
#define FATAL(A)
Definition: JMessage.hh:67
int debug
debug level
Definition: JSirene.cc:69
#define WARNING(A)
Definition: JMessage.hh:65
ROOT I/O of application specific meta data.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2142
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:72
Utility class to parse parameter values.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Auxiliary class to define a range between two values.
Scanning of objects from a single file according a format that follows from the extension of each fil...
Direct access to string in detector data structure.
Constants.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Detector data structure.
Definition: JDetector.hh:96
Auxiliary class for map of PMT parameters.
Data structure for PMT parameters.
Utility class to parse parameter values.
Definition: JProperties.hh:501
General exception.
Definition: JException.hh:24
Implementation of object output from STD container.
Utility class to parse command line options.
Definition: JParser.hh:1698
ROOT file reader.
Object reading from a list of files.
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
static const char *const _2F
Name extension for 2F rate fitted.
@ FIT_PMTS_QE_FIXED_t
fit parameters of PMTs with QE fixed
Definition: JFitK40.hh:56
@ FIT_PMTS_AND_ANGULAR_DEPENDENCE_t
fit parameters of PMTs and angular dependence of K40 rate
Definition: JFitK40.hh:54
@ FIT_MODEL_t
fit parameters of K40 rate and TTSs of PMTs
Definition: JFitK40.hh:57
@ FIT_PMTS_AND_BACKGROUND_t
fit parameters of PMTs and background
Definition: JFitK40.hh:55
@ FIT_PMTS_t
fit parameters of PMTs
Definition: JFitK40.hh:53
static const char *const _2R
Name extension for 2D rate measured.
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
Definition: JLocation.hh:247
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
double getSurvivalProbability(const JPMTParameters &parameters)
Get model dependent probability that a one photo-electron hit survives the simulation of the PMT assu...
int getIndex()
Get index for user I/O manipulation.
Definition: JManip.hh:26
T & getInstance(const T &object)
Get static instance from temporary object.
Definition: JObject.hh:75
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
Definition: JSTDTypes.hh:14
Auxiliary data structure for sequence of same character.
Definition: JManip.hh:330
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
Type definition of range.
Definition: JHead.hh:43
Livetime of noise data.
Definition: JHead.hh:1062
Detector file.
Definition: JHead.hh:227
Acoustic single fit.
Model for fit to acoustics data.
Fit parameters for two-fold coincidence rate due to K40.
Definition: JFitK40.hh:610
Auxiliary class for TDC constraints.
Definition: JTDC_t.hh:39
range_type equal_range(const int id)
Get range of constraints for given module.
Definition: JTDC_t.hh:101
void reverse()
Reverse constraints.
Definition: JTDC_t.hh:137
bool is_valid(const bool option=false) const
Check validity of TDC constrants.
Definition: JTDC_t.hh:171
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition: JFitK40.hh:103
Data structure for measured coincidence rate of pair of PMTs.
Definition: JFitK40.hh:66
Router for mapping of string identifier to index.
JComment & add(const std::string &comment)
Add comment.
Definition: JComment.hh:100
void store(const char *file_name) const
Store to output file.
void load(const char *file_name)
Load from input file.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:68
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:72
Data structure for a pair of indices.