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