Jpp 21.0.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
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 <memory>
#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 64 of file JFitK40.cc.

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
#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
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2140
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:57
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
#define TEST(T)
Test data member.
Definition JRootClass.cc:35
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.