Jpp 19.3.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
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#include <memory>
9
10#include "TROOT.h"
11#include "TFile.h"
12#include "TH1D.h"
13#include "TH2D.h"
14
17
18#include "JTools/JConstants.hh"
23
25#include "JROOT/JRootToolkit.hh"
26
27#include "JTools/JRange.hh"
28
30
31#include "JSupport/JMeta.hh"
33
35#include "JCalibrate/JTDC_t.hh"
36
37#include "JFitK40.hh"
38
39#include "Jeep/JProperties.hh"
40#include "Jeep/JPrint.hh"
41#include "Jeep/JParser.hh"
42#include "Jeep/JMessage.hh"
43
44
45namespace {
46
47 double MINIMAL_RATE_HZ = 1.0e-2; //!< minimal coincidence rate [Hz] (below, disable PMT)
48 double STDEV = 2.0; //!< number of standard deviations
49 size_t MAXIMAL_COUNTS = 10; //!< maximal count of too low rates
50 double QE_MIN = 0.1; //!< minimal QE (below, disable PMT)
51 double T0_NS = 1.0; //!< precision t0 [ns] (at limit, disable PMT)
52}
53
54
55/**
56 * \file
57 *
58 * Auxiliary program to fit PMT parameters from JMergeCalibrateK40.cc output.\n
59 * Note that the contebts of the 2D histograms are converted to standard deviations.
60 *
61 * \author mdejong
62 */
63int main(int argc, char **argv)
64{
65 using namespace std;
66 using namespace JPP;
67 using namespace KM3NETDAQ;
68
70
72
73 string inputFile;
74 string outputFile;
75 string detectorFile;
76 string pmtFile;
77 JTDC_t TDC;
78 bool reverse;
79 bool overwriteDetector;
80 JCounter writeFits;
81 bool fitAngle;
82 bool fitNoise;
83 bool fitModel;
84 bool fixQE;
85 JRange_t X;
86 JRange_t Y;
87 int debug;
88
89 try {
90
91 JProperties properties;
92
93 properties.insert(gmake_property(MINIMAL_RATE_HZ));
94 properties.insert(gmake_property(STDEV));
95 properties.insert(gmake_property(MAXIMAL_COUNTS));
96 properties.insert(gmake_property(QE_MIN));
97 properties.insert(gmake_property(T0_NS));
98 properties.insert(gmake_property(K40.R));
99 properties.insert(gmake_property(K40.p1));
100 properties.insert(gmake_property(K40.p2));
101 properties.insert(gmake_property(K40.p3));
102 properties.insert(gmake_property(K40.p4));
103
104 JParser<> zap("Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
105
106 zap['@'] = make_field(properties) = JPARSER::initialised();
107 zap['f'] = make_field(inputFile, "input file (output from JMergeCalibrateK40).");
108 zap['o'] = make_field(outputFile, "output file.") = "fit.root";
109 zap['a'] = make_field(detectorFile, "detector file.");
110 zap['P'] = make_field(pmtFile, "specify PMT file name that can be given to JTriggerEfficiency.") = "";
111 zap['!'] = make_field(TDC,
112 "Fix time offset(s) of PMT(s) of certain module(s), e.g."
113 "\n-! \"808969848 0 808982077 23\" will fix time offsets of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
114 "\nSame PMT offset can be fixed for all optical modules, e.g."
115 "\n-! \"-1 0 -1 23\" will fix time offsets of PMTs 0 and 23 of all optical modules.") = JPARSER::initialised();
116 zap['r'] = make_field(reverse, "reverse TDC constraints due to option -! <TDC>.");
117 zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with fitted time offsets.");
118 zap['w'] = make_field(writeFits, "write fit results to ROOT file; -ww also write fitted TTS to PMT file.");
119 zap['D'] = make_field(fitAngle, "fit angular distribution; fix normalisation.");
120 zap['B'] = make_field(fitNoise, "fit background.");
121 zap['M'] = make_field(fitModel, "fit angular distribution as well as normalisation; fix PMT QEs = 1.0.");
122 zap['Q'] = make_field(fixQE, "fix PMT QEs = 1.0.");
123 zap['x'] = make_field(X, "fit range (PMT pairs).") = JRange_t();
124 zap['y'] = make_field(Y, "fit range (time residual [ns]).") = JRange_t();
125 zap['d'] = make_field(debug, "debug.") = 1;
126
127 zap(argc, argv);
128 }
129 catch(const exception &error) {
130 FATAL(error.what() << endl);
131 }
132
133
134 if ((fitModel ? 1 : 0) +
135 (fitAngle ? 1 : 0) +
136 (fitNoise ? 1 : 0) +
137 (fixQE ? 1 : 0) > 1) {
138 FATAL("Use either option -M, -D, -B or -Q" << endl);
139 }
140
141 const int option = (fitModel ? FIT_MODEL_t :
143 fitNoise ? FIT_PMTS_AND_BACKGROUND_t :
144 fixQE ? FIT_PMTS_QE_FIXED_t :
145 FIT_PMTS_t);
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 try {
156 TDC.is_valid(true);
157 }
158 catch(const exception &error) {
159 FATAL(error.what() << endl);
160 }
161
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 exception& error) {}
179 }
180
181 parameters.comment.add(JMeta(argc, argv));
182
183
184 TFile* in = TFile::Open(inputFile.c_str(), "exist");
185
186 if (in == NULL || !in->IsOpen()) {
187 FATAL("File: " << inputFile << " not opened." << endl);
188 }
189
190
191 TFile out(outputFile.c_str(), "recreate");
192
193 TH1D h0("chi2", NULL, 500, 0.0, 5.0);
194 TH1D hn("hn", NULL, 501, -0.5, 500.0);
195 TH1D hr("rate", NULL, 500, 0.0, 25.0);
196 TH1D h1("p1", NULL, 500, -5.0, +5.0);
197 TH1D h2("p2", NULL, 500, -5.0, +5.0);
198 TH1D h3("p3", NULL, 500, -5.0, +5.0);
199 TH1D h4("p4", NULL, 500, -5.0, +5.0);
200 TH1D hc("cc", NULL, 500, -0.1, +0.1);
201
203 const JStringRouter string(detector);
204
205 TH2D H2("detector", NULL,
206 string.size() + 0, -0.5, string.size() - 0.5,
207 range.getUpperLimit(), 1 - 0.5, range.getUpperLimit() + 0.5);
208
209 for (Int_t i = 1; i <= H2.GetXaxis()->GetNbins(); ++i) {
210 H2.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1)));
211 }
212 for (Int_t i = 1; i <= H2.GetYaxis()->GetNbins(); ++i) {
213 H2.GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i));
214 }
215
216 TH2D* HN = (TH2D*) H2.Clone("iterations");
217
218 JFit fit(0, debug);
219
220 for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
221
222 if (module->getFloor() == 0) {
223 continue;
224 }
225
226 const JTDC_t::range_type range = TDC.equal_range(module->getID());
227
228 NOTICE("Module " << setw(10) << module->getID() << ' ' << getLabel(module->getLocation()) << " !" << distance(range.first, range.second) << endl);
229
230 TH2D* h2d = (TH2D*) in->Get(MAKE_CSTRING(module->getID() << _2R));
231
232 if (h2d == NULL || h2d->GetEntries() == 0) {
233
234 NOTICE("No data for module " << module->getID() << " -> set QEs to 0." << endl);
235
236 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
237 parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
238 }
239
240 continue;
241 }
242
243 JModel model(*module, K40, range, option);
244
245 data_type data; // input data
246
247 vector<size_t> count[2] = { vector<size_t>(NUMBER_OF_PMTS, 0), // counter (exclude self)
248 vector<size_t>(NUMBER_OF_PMTS, 1) }; // counter (include self)
249
250 for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
251
252 const pair_type pair = model.getPair(ix - 1);
253
254 auto& buffer = data[pair]; // storage
255
256 double V = 0.0; // integrated value
257 double W = 0.0; // integrated error
258
259 for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
260
261 const double x = h2d->GetXaxis()->GetBinCenter(ix);
262 const double y = h2d->GetYaxis()->GetBinCenter(iy);
263
264 if (X(x) && Y(y)) {
265
266 double value = h2d->GetBinContent(ix,iy);
267 double error = h2d->GetBinError (ix,iy);
268
269 buffer.push_back(rate_type(y, value, error));
270
271 double width = h2d->GetYaxis()->GetBinWidth(iy);
272
273 value *= width;
274 error *= width;
275
276 V += value;
277 W += error * error;
278 }
279 }
280
281 W = sqrt(W);
282
283 if (V <= 0.0 - STDEV*W) {
284 count[0][pair.first] += 1;
285 count[0][pair.second] += 1;
286 }
287
288 if (V <= MINIMAL_RATE_HZ + STDEV*W) {
289 count[1][pair.first] += 1;
290 count[1][pair.second] += 1;
291 }
292 }
293
294 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
295
296 if (count[0][pmt] >= MAXIMAL_COUNTS) { // too many paired rates negative
297
298 WARNING("PMT " << setw(10) << module->getID() << '.' << FILL(2,'0') << pmt << FILL() << " some rates negative -> fit background" << endl);
299
300 if (fit.value.parameters[pmt].status) {
301 model.parameters[pmt].bg.set();
302 }
303 }
304
305 if (count[1][pmt] == NUMBER_OF_PMTS) { // all paired rates too low
306
307 WARNING("PMT " << setw(10) << module->getID() << '.' << FILL(2,'0') << pmt << FILL() << " all rates to low -> disable" << endl);
308
309 model.parameters[pmt].disable();
310 }
311 }
312
313 DEBUG("Start value:" << endl << model << endl);
314
315 try {
316
317 fit.value = model; // start value
318
319 auto result = fit(data);
320
321 if (result.ndf <= 0) {
322
323 ERROR("Fit result " << setw(10) << module->getID() << " NDF " << setw(5) << result.ndf << " -> skip" << endl);
324
325 continue;
326 }
327
328 bool refit = false;
329
330 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
331
332 if (fit.value.parameters[pmt].status) {
333
334 if (fit.value.parameters[pmt].QE() <= QE_MIN ||
335 fit.value.parameters[pmt].QE() <= 0.0 + STDEV * fit.error.parameters[pmt].QE()) {
336
337 WARNING("PMT " << setw(10) << module->getID() << '.' << FILL(2,'0') << pmt << FILL() << ' '
338 << "QE = "
339 << FIXED(5,3) << fit.value.parameters[pmt].QE() << " +/- "
340 << FIXED(5,3) << fit.error.parameters[pmt].QE() << " "
341 << " -> disable" << (!refit ? " and refit" : "") << endl);
342
343 fit.value.parameters[pmt].disable();
344
345 refit = true;
346 }
347
348 if (fit.value.parameters[pmt].t0.atLimit(T0_NS)) {
349
350 WARNING("PMT " << setw(10) << module->getID() << '.' << FILL(2,'0') << pmt << FILL() << ' '
351 << "t0 at limit "
352 << FIXED(5,3) << fit.value.parameters[pmt].t0() << " +/- "
353 << FIXED(5,3) << fit.error.parameters[pmt].t0() << " "
354 << " -> disable" << (!refit ? " and refit" : "") << endl);
355
356 fit.value.parameters[pmt].disable();
357
358 refit = true;
359 }
360 }
361 }
362
363 if (refit) {
364
365 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
366 if (fit.value.parameters[pmt].status) {
367 fit.value.parameters[pmt].set(JPMTParameters_t::getInstance());
368 }
369 }
370
371 refit = false;
372 result = fit(data);
373 }
374
375 NOTICE("Fit result " << setw(10) << module->getID() << " chi2 / NDF " << FIXED(10,2) << result.chi2 << " / " << setw(5) << result.ndf << ' ' << setw(5) << fit.numberOfIterations << endl);
376
377 DEBUG(fit.value);
378
379 if (writeFits) {
380
381 h0.Fill(result.chi2 / result.ndf);
382 hn.Fill(fit.numberOfIterations);
383 hr.Fill(fit.value.R );
384 h1.Fill(fit.value.p1);
385 h2.Fill(fit.value.p2);
386 h3.Fill(fit.value.p3);
387 h4.Fill(fit.value.p4);
388 hc.Fill(fit.value.cc);
389
390 TH1D h1t(MAKE_CSTRING(module->getID() << ".1t0"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
391 TH1D h1s(MAKE_CSTRING(module->getID() << ".1TTS"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
392 TH1D h1q(MAKE_CSTRING(module->getID() << ".1QE"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
393
394 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
395 h1t.SetBinContent(pmt + 1, fit.value.parameters[pmt].t0 ());
396 h1t.SetBinError (pmt + 1, fit.error.parameters[pmt].t0 () + numeric_limits<double>::epsilon());
397 h1s.SetBinContent(pmt + 1, fit.value.parameters[pmt].TTS());
398 h1s.SetBinError (pmt + 1, fit.error.parameters[pmt].TTS() + numeric_limits<double>::epsilon());
399 h1q.SetBinContent(pmt + 1, fit.value.parameters[pmt].QE ());
400 h1q.SetBinError (pmt + 1, fit.error.parameters[pmt].QE () + numeric_limits<double>::epsilon());
401 }
402
403 out << h1t << h1s << h1q;
404
405 for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
406
407 const pair_type pair = fit.value.getPair(ix - 1);
408
409 for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
410
411 const double dt_ns = h2d->GetYaxis()->GetBinCenter(iy);
412
413 h2d->SetBinContent(ix, iy, fit.value.getValue(pair, dt_ns));
414 h2d->SetBinError (ix, iy, 0.0);
415 }
416 }
417
418 h2d->SetName(MAKE_CSTRING(module->getID() << _2F));
419 h2d->Write();
420
421 const double x = string.getIndex(module->getString());
422 const double y = module->getFloor();
423
424 H2 .Fill(x, y, result.chi2 / result.ndf);
425 HN->Fill(x, y, fit.numberOfIterations);
426 }
427
428 const double t0 = (fit.value.hasFixedTimeOffset() ? fit.value.getFixedTimeOffset() : 0.0);
429
430 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
431
432 JPMTParameters& data = parameters[JPMTIdentifier(module->getID(), pmt)];
433
434 const double P = getSurvivalProbability(data);
435
436 if (P > 0.0)
437 data.QE = fit.value.parameters[pmt].QE / P;
438 else
439 data.QE = 0.0;
440
441 if (writeFits > 1) {
442 data.TTS_ns = fit.value.parameters[pmt].TTS();
443 }
444
445 module->getPMT(pmt).addT0(fit.value.parameters[pmt].t0.get() - t0);
446 }
447 }
448 catch(const exception& error) {
449
450 ERROR("Module " << setw(10) << module->getID() << ' ' << error.what() << " -> set QEs to 0." << endl);
451
452 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
453 parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
454 }
455 }
456 }
457
458
459 vector<JMeta> meta(1, JMeta(argc, argv));
460
461 {
462 JSingleFileScanner<JMeta> reader(inputFile);
463 JSTDObjectWriter <JMeta> writer(meta);
464
465 writer << reader;
466 }
467
468 for (vector<JMeta>::const_reverse_iterator i = meta.rbegin(); i != meta.rend(); ++i) {
469 parameters.comment.add(*i);
470 detector .comment.add(*i);
471 }
472
473 if (overwriteDetector) {
474
475 NOTICE("Store calibration data on file " << detectorFile << endl);
476
477 store(detectorFile, detector);
478 }
479
480 if (pmtFile != "") {
481 parameters.store(pmtFile.c_str());
482 }
483
484
485 for (vector<JMeta>::const_iterator i = meta.begin(); i != meta.end(); ++i) {
486 putObject(out,*i);
487 }
488
489 for (JRootFileReader<JDAQHeader> in(inputFile.c_str()); in.hasNext(); ) {
490 putObject(out, *in.next());
491 }
492
493 if (writeFits) {
494 out << h0 << hn << hr << h1 << h2 << h3 << h4 << hc << H2 << *HN;
495 }
496
497 out.Close();
498
499 return 0;
500}
string outputFile
KM3NeT DAQ constants, bit handling, etc.
Data structure for detector geometry and calibration.
int main(int argc, char **argv)
Definition JFitK40.cc:63
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define NOTICE(A)
Definition JMessage.hh:64
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
#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.
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
Object reading from a list of files.
Range of values.
Definition JRange.hh:42
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...
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
Definition JGandalf.hh:57
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.
return result
Definition JPolint.hh:862
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
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:678
static const JK40Parameters & getInstance()
Get default values.
Definition JFitK40.hh:695
static const JPMTParameters_t & getInstance()
Get default values.
Definition JFitK40.hh:477
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.