Jpp  19.1.0-rc.1
the software that should make you happy
Functions
JFitK40.cc File Reference

Auxiliary program to fit PMT parameters from JMergeCalibrateK40.cc output. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <set>
#include <vector>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "km3net-dataformat/online/JDAQ.hh"
#include "km3net-dataformat/online/JDAQHeader.hh"
#include "JTools/JConstants.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JStringRouter.hh"
#include "JROOT/JRootFileWriter.hh"
#include "JROOT/JRootToolkit.hh"
#include "JTools/JRange.hh"
#include "JLang/JSTDObjectWriter.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JSingleFileScanner.hh"
#include "JCalibrate/JCalibrateK40.hh"
#include "JCalibrate/JTDC_t.hh"
#include "JFitK40.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Auxiliary program to fit PMT parameters from JMergeCalibrateK40.cc output.


Note that the contebts of the 2D histograms are converted to standard deviations.

Author
mdejong

Definition in file JFitK40.cc.

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 61 of file JFitK40.cc.

62 {
63  using namespace std;
64  using namespace JPP;
65  using namespace KM3NETDAQ;
66 
67  typedef JRange<double> JRange_t;
68 
69  JK40Parameters K40 = JK40Parameters::getInstance();
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 :
139  fitAngle ? FIT_PMTS_AND_ANGULAR_DEPENDENCE_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 
159  JDetector detector;
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 }
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
ERROR
#define ERROR(A)
Definition: JMessage.hh:66
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
debug
int debug
debug level
Definition: JSirene.cc:69
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:136
gmake_property
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Definition: JProperties.hh:1293
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition: PhysicsEvent.hh:434
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:96
JDETECTOR::JPMTParametersMap
Auxiliary class for map of PMT parameters.
Definition: JPMTParametersMap.hh:93
JDETECTOR::JPMTParameters
Data structure for PMT parameters.
Definition: JPMTParameters.hh:29
JEEP::JProperties
Utility class to parse parameter values.
Definition: JProperties.hh:501
JLANG::JException
General exception.
Definition: JException.hh:24
JLANG::JSTDObjectWriter
Implementation of object output from STD container.
Definition: JSTDObjectWriter.hh:27
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1714
JROOT::JRootFileReader
ROOT file reader.
Definition: JRootFileReader.hh:54
JSUPPORT::JSingleFileScanner
Object reading from a list of files.
Definition: JSingleFileScanner.hh:172
JTOOLS::JRange::getUpperLimit
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
JCALIBRATE::_2F
static const char *const _2F
Name extension for 2F rate fitted.
Definition: JCalibrateK40.hh:37
JCALIBRATE::FIT_PMTS_QE_FIXED_t
@ FIT_PMTS_QE_FIXED_t
fit parameters of PMTs with QE fixed
Definition: JFitK40.hh:56
JCALIBRATE::FIT_PMTS_AND_ANGULAR_DEPENDENCE_t
@ FIT_PMTS_AND_ANGULAR_DEPENDENCE_t
fit parameters of PMTs and angular dependence of K40 rate
Definition: JFitK40.hh:54
JCALIBRATE::FIT_MODEL_t
@ FIT_MODEL_t
fit parameters of K40 rate and TTSs of PMTs
Definition: JFitK40.hh:57
JCALIBRATE::FIT_PMTS_AND_BACKGROUND_t
@ FIT_PMTS_AND_BACKGROUND_t
fit parameters of PMTs and background
Definition: JFitK40.hh:55
JCALIBRATE::FIT_PMTS_t
@ FIT_PMTS_t
fit parameters of PMTs
Definition: JFitK40.hh:53
JCALIBRATE::_2R
static const char *const _2R
Name extension for 2D rate measured.
Definition: JCalibrateK40.hh:36
JDETECTOR::getLabel
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
Definition: JLocation.hh:246
JDETECTOR::getRangeOfFloors
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Definition: JDetectorToolkit.hh:246
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:484
JDETECTOR::store
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
Definition: JDetectorToolkit.hh:581
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
JLANG::getIndex
int getIndex()
Get index for user I/O manipulation.
Definition: JManip.hh:26
JLANG::getInstance
T & getInstance(const T &object)
Get static instance from temporary object.
Definition: JObject.hh:75
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
JROOT::putObject
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
Definition: JRootFileWriter.hh:41
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
KM3NETDAQ::NUMBER_OF_PMTS
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
std
Definition: JSTDTypes.hh:14
FILL
Auxiliary data structure for sequence of same character.
Definition: JManip.hh:330
FIXED
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
JAANET::JRange_t
Type definition of range.
Definition: JHead.hh:43
JAANET::K40
Livetime of noise data.
Definition: JHead.hh:1062
JAANET::detector
Detector file.
Definition: JHead.hh:227
JACOUSTICS::JFit
Acoustic single fit.
Definition: JAcoustics/JEvt.hh:37
JACOUSTICS::JModel
Model for fit to acoustics data.
Definition: JAcoustics/JModel.hh:480
JCALIBRATE::JK40Parameters
Fit parameters for two-fold coincidence rate due to K40.
Definition: JFitK40.hh:610
JCALIBRATE::JTDC_t
Auxiliary class for TDC constraints.
Definition: JTDC_t.hh:39
JCALIBRATE::JTDC_t::equal_range
range_type equal_range(const int id)
Get range of constraints for given module.
Definition: JTDC_t.hh:101
JCALIBRATE::JTDC_t::reverse
void reverse()
Reverse constraints.
Definition: JTDC_t.hh:137
JCALIBRATE::JTDC_t::is_valid
bool is_valid(const bool option=false) const
Check validity of TDC constrants.
Definition: JTDC_t.hh:171
JCALIBRATE::data_type
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition: JFitK40.hh:103
JCALIBRATE::rate_type
Data structure for measured coincidence rate of pair of PMTs.
Definition: JFitK40.hh:66
JDETECTOR::JStringRouter
Router for mapping of string identifier to index.
Definition: JStringRouter.hh:23
JEEP::JComment::add
JComment & add(const std::string &comment)
Add comment.
Definition: JComment.hh:100
JLANG::JObjectStreamIO::store
void store(const char *file_name) const
Store to output file.
Definition: JObjectStreamIO.hh:41
JLANG::JObjectStreamIO::load
void load(const char *file_name)
Load from input file.
Definition: JObjectStreamIO.hh:30
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:84
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:72
JTOOLS::JCombinatorics::pair_type
Data structure for a pair of indices.
Definition: JCombinatorics.hh:29
Generated by doxygen 1.9.1