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