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