Jpp 19.3.0-rc.5
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 63 of file JFitK40.cc.

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