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