Jpp master_rocky-44-g75b7c4f75
the software that should make you happy
Loading...
Searching...
No Matches
Functions
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 "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 61 of file JFitK40.cc.

62{
63 using namespace std;
64 using namespace JPP;
65 using namespace KM3NETDAQ;
66
67 typedef JRange<double> JRange_t;
68
69 JK40Parameters K40 = JK40Parameters::getInstance();
70
71 string inputFile;
72 string outputFile;
73 string detectorFile;
74 string pmtFile;
75 JTDC_t TDC;
76 bool reverse;
77 bool overwriteDetector;
78 JCounter writeFits;
79 bool fitAngle;
80 bool fitNoise;
81 bool fitModel;
82 bool fixQE;
83 JRange_t X;
84 JRange_t Y;
85 int debug;
86
87 try {
88
89 JProperties properties;
90
91 properties.insert(gmake_property(MINIMAL_RATE_HZ));
92 properties.insert(gmake_property(STDEV));
93 properties.insert(gmake_property(MAXIMAL_COUNTS));
94 properties.insert(gmake_property(QE_MIN));
95 properties.insert(gmake_property(K40.R));
96 properties.insert(gmake_property(K40.p1));
97 properties.insert(gmake_property(K40.p2));
98 properties.insert(gmake_property(K40.p3));
99 properties.insert(gmake_property(K40.p4));
100
101 JParser<> zap("Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
102
103 zap['@'] = make_field(properties) = JPARSER::initialised();
104 zap['f'] = make_field(inputFile, "input file (output from JMergeCalibrateK40).");
105 zap['o'] = make_field(outputFile, "output file.") = "fit.root";
106 zap['a'] = make_field(detectorFile, "detector file.");
107 zap['P'] = make_field(pmtFile, "specify PMT file name that can be given to JTriggerEfficiency.") = "";
108 zap['!'] = make_field(TDC,
109 "Fix time offset(s) of PMT(s) of certain module(s), e.g."
110 "\n-! \"808969848 0 808982077 23\" will fix time offsets of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
111 "\nSame PMT offset can be fixed for all optical modules, e.g."
112 "\n-! \"-1 0 -1 23\" will fix time offsets of PMTs 0 and 23 of all optical modules.") = JPARSER::initialised();
113 zap['r'] = make_field(reverse, "reverse TDC constraints due to option -! <TDC>.");
114 zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with fitted time offsets.");
115 zap['w'] = make_field(writeFits, "write fit results to ROOT file; -ww also write fitted TTS to PMT file.");
116 zap['D'] = make_field(fitAngle, "fit angular distribution; fix normalisation.");
117 zap['B'] = make_field(fitNoise, "fit background.");
118 zap['M'] = make_field(fitModel, "fit angular distribution as well as normalisation; fix PMT QEs = 1.0.");
119 zap['Q'] = make_field(fixQE, "fix PMT QEs = 1.0.");
120 zap['x'] = make_field(X, "fit range (PMT pairs).") = JRange_t();
121 zap['y'] = make_field(Y, "fit range (time residual [ns]).") = JRange_t();
122 zap['d'] = make_field(debug, "debug.") = 1;
123
124 zap(argc, argv);
125 }
126 catch(const exception &error) {
127 FATAL(error.what() << endl);
128 }
129
130
131 if ((fitModel ? 1 : 0) +
132 (fitAngle ? 1 : 0) +
133 (fitNoise ? 1 : 0) +
134 (fixQE ? 1 : 0) > 1) {
135 FATAL("Use either option -M, -D, -B or -Q" << endl);
136 }
137
138 const int option = (fitModel ? FIT_MODEL_t :
139 fitAngle ? FIT_PMTS_AND_ANGULAR_DEPENDENCE_t :
140 fitNoise ? FIT_PMTS_AND_BACKGROUND_t :
141 fixQE ? FIT_PMTS_QE_FIXED_t :
142 FIT_PMTS_t);
143
144 if (reverse) {
145 TDC.reverse();
146 }
147
148 for (JTDC_t::const_iterator i = TDC.begin(); i != TDC.end(); ++i) {
149 DEBUG("PMT " << setw(10) << i->first << ' ' << setw(2) << i->second << " constrain t0." << endl);
150 }
151
152 try {
153 TDC.is_valid(true);
154 }
155 catch(const exception &error) {
156 FATAL(error.what() << endl);
157 }
158
159 JDetector detector;
160
161 try {
162 load(detectorFile, detector);
163 }
164 catch(const JException& error) {
165 FATAL(error);
166 }
167
168
169 JPMTParametersMap parameters;
170
171 if (pmtFile != "") {
172 try {
173 parameters.load(pmtFile.c_str());
174 }
175 catch(const exception& error) {}
176 }
177
178 parameters.comment.add(JMeta(argc, argv));
179
180
181 TFile* in = TFile::Open(inputFile.c_str(), "exist");
182
183 if (in == NULL || !in->IsOpen()) {
184 FATAL("File: " << inputFile << " not opened." << endl);
185 }
186
187
188 TFile out(outputFile.c_str(), "recreate");
189
190 TH1D h0("chi2", NULL, 500, 0.0, 5.0);
191 TH1D hn("hn", NULL, 501, -0.5, 500.0);
192 TH1D hr("rate", NULL, 500, 0.0, 25.0);
193 TH1D h1("p1", NULL, 500, -5.0, +5.0);
194 TH1D h2("p2", NULL, 500, -5.0, +5.0);
195 TH1D h3("p3", NULL, 500, -5.0, +5.0);
196 TH1D h4("p4", NULL, 500, -5.0, +5.0);
197 TH1D hc("cc", NULL, 500, -0.1, +0.1);
198
199 const floor_range range = getRangeOfFloors(detector);
200 const JStringRouter string(detector);
201
202 TH2D H2("detector", NULL,
203 string.size() + 0, -0.5, string.size() - 0.5,
204 range.getUpperLimit(), 1 - 0.5, range.getUpperLimit() + 0.5);
205
206 for (Int_t i = 1; i <= H2.GetXaxis()->GetNbins(); ++i) {
207 H2.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1)));
208 }
209 for (Int_t i = 1; i <= H2.GetYaxis()->GetNbins(); ++i) {
210 H2.GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i));
211 }
212
213
214 JFit fit(0, debug);
215
216 for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
217
218 if (module->getFloor() == 0) {
219 continue;
220 }
221
222 const JTDC_t::range_type range = TDC.equal_range(module->getID());
223
224 NOTICE("Module " << setw(10) << module->getID() << ' ' << getLabel(module->getLocation()) << " !" << distance(range.first, range.second) << endl);
225
226 TH2D* h2d = (TH2D*) in->Get(MAKE_CSTRING(module->getID() << _2R));
227
228 if (h2d == NULL || h2d->GetEntries() == 0) {
229
230 NOTICE("No data for module " << module->getID() << " -> set QEs to 0." << endl);
231
232 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
233 parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
234 }
235
236 continue;
237 }
238
239 JModel model(*module, K40, range, option);
240
241 data_type data; // input data
242
243 vector<size_t> count[2] = { vector<size_t>(NUMBER_OF_PMTS, 0), // counter (exclude self)
244 vector<size_t>(NUMBER_OF_PMTS, 1) }; // counter (include self)
245
246 for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
247
248 const pair_type pair = model.getPair(ix - 1);
249
250 auto& buffer = data[pair]; // storage
251
252 double V = 0.0; // integrated value
253 double W = 0.0; // integrated error
254
255 for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
256
257 const double x = h2d->GetXaxis()->GetBinCenter(ix);
258 const double y = h2d->GetYaxis()->GetBinCenter(iy);
259
260 if (X(x) && Y(y)) {
261
262 double value = h2d->GetBinContent(ix,iy);
263 double error = h2d->GetBinError (ix,iy);
264
265 buffer.push_back(rate_type(y, value, error));
266
267 double width = h2d->GetYaxis()->GetBinWidth(iy);
268
269 value *= width;
270 error *= width;
271
272 V += value;
273 W += error * error;
274 }
275 }
276
277 W = sqrt(W);
278
279 if (V <= 0.0 - STDEV*W) {
280 count[0][pair.first] += 1;
281 count[0][pair.second] += 1;
282 }
283
284 if (V <= MINIMAL_RATE_HZ + STDEV*W) {
285 count[1][pair.first] += 1;
286 count[1][pair.second] += 1;
287 }
288 }
289
290 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
291
292 if (count[0][pmt] >= MAXIMAL_COUNTS) { // too many paired rates negative
293
294 WARNING("PMT " << setw(10) << module->getID() << '.' << FILL(2,'0') << pmt << FILL() << " some rates negative -> fit background" << endl);
295
296 if (fit.value.parameters[pmt].status) {
297 model.parameters[pmt].bg.set();
298 }
299 }
300
301 if (count[1][pmt] == NUMBER_OF_PMTS) { // all paired rates too low
302
303 WARNING("PMT " << setw(10) << module->getID() << '.' << FILL(2,'0') << pmt << FILL() << " all rates to low -> disable" << endl);
304
305 model.parameters[pmt].disable();
306 }
307 }
308
309 DEBUG("Start value:" << endl << model << endl);
310
311 try {
312
313 fit.value = model; // start value
314
315 auto result = fit(data);
316
317 bool refit = false;
318
319 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
320
321 if (fit.value.parameters[pmt].QE() < QE_MIN && fit.value.parameters[pmt].status) {
322
323 WARNING("PMT " << setw(10) << module->getID() << '.' << FILL(2,'0') << pmt << FILL() << ' ' << FIXED(5,3) << fit.value.parameters[pmt].QE() << " < " << FIXED(5,3) << QE_MIN << " -> disable" << (!refit ? " and refit" : "") << endl);
324
325 fit.value.parameters[pmt].disable();
326
327 refit = true;
328 }
329 }
330
331 if (refit) {
332 refit = false;
333 result = fit(data);
334 }
335
336 NOTICE("Fit result " << setw(10) << module->getID() << " chi2 / NDF " << FIXED(10,2) << result.chi2 << " / " << setw(5) << result.ndf << ' ' << setw(5) << fit.numberOfIterations << endl);
337
338 DEBUG(fit.value);
339
340 if (writeFits) {
341
342 h0.Fill(result.chi2 / result.ndf);
343 hn.Fill(fit.numberOfIterations);
344 hr.Fill(fit.value.R );
345 h1.Fill(fit.value.p1);
346 h2.Fill(fit.value.p2);
347 h3.Fill(fit.value.p3);
348 h4.Fill(fit.value.p4);
349 hc.Fill(fit.value.cc);
350
351 TH1D h1t(MAKE_CSTRING(module->getID() << ".1t0"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
352 TH1D h1s(MAKE_CSTRING(module->getID() << ".1TTS"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
353 TH1D h1q(MAKE_CSTRING(module->getID() << ".1QE"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
354
355 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
356 h1t.SetBinContent(pmt + 1, fit.value.parameters[pmt].t0 ());
357 h1s.SetBinContent(pmt + 1, fit.value.parameters[pmt].TTS());
358 h1q.SetBinContent(pmt + 1, fit.value.parameters[pmt].QE ());
359 }
360
361 out << h1t << h1s << h1q;
362
363 for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
364
365 const pair_type pair = fit.value.getPair(ix - 1);
366
367 for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
368
369 const double dt_ns = h2d->GetYaxis()->GetBinCenter(iy);
370
371 h2d->SetBinContent(ix, iy, fit.value.getValue(pair, dt_ns));
372 h2d->SetBinError (ix, iy, 0.0);
373 }
374 }
375
376 h2d->SetName(MAKE_CSTRING(module->getID() << _2F));
377 h2d->Write();
378
379 H2.Fill((double) string.getIndex(module->getString()), (double) module->getFloor(), result.chi2 / result.ndf);
380 }
381
382 const double t0 = (fit.value.hasFixedTimeOffset() ? fit.value.getFixedTimeOffset() : 0.0);
383
384 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
385
386 JPMTParameters& data = parameters[JPMTIdentifier(module->getID(), pmt)];
387
388 const double P = getSurvivalProbability(data);
389
390 if (P > 0.0)
391 data.QE = fit.value.parameters[pmt].QE / P;
392 else
393 data.QE = 0.0;
394
395 if (writeFits > 1) {
396 data.TTS_ns = fit.value.parameters[pmt].TTS();
397 }
398
399 module->getPMT(pmt).addT0(fit.value.parameters[pmt].t0.get() - t0);
400 }
401 }
402 catch(const exception& error) {
403
404 ERROR("Module " << setw(10) << module->getID() << ' ' << error.what() << " -> set QEs to 0." << endl);
405
406 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
407 parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
408 }
409 }
410 }
411
412
413 vector<JMeta> meta(1, JMeta(argc, argv));
414
415 {
416 JSingleFileScanner<JMeta> reader(inputFile);
417 JSTDObjectWriter <JMeta> writer(meta);
418
419 writer << reader;
420 }
421
422 for (vector<JMeta>::const_reverse_iterator i = meta.rbegin(); i != meta.rend(); ++i) {
423 parameters.comment.add(*i);
424 detector .comment.add(*i);
425 }
426
427 if (overwriteDetector) {
428
429 NOTICE("Store calibration data on file " << detectorFile << endl);
430
431 store(detectorFile, detector);
432 }
433
434 if (pmtFile != "") {
435 parameters.store(pmtFile.c_str());
436 }
437
438
439 for (vector<JMeta>::const_iterator i = meta.begin(); i != meta.end(); ++i) {
440 putObject(out,*i);
441 }
442
443 for (JRootFileReader<JDAQHeader> in(inputFile.c_str()); in.hasNext(); ) {
444 putObject(out, *in.next());
445 }
446
447 if (writeFits) {
448 out << h0 << hn << hr << h1 << h2 << h3 << h4 << hc << H2;
449 }
450
451 out.Close();
452
453 return 0;
454}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
NOTICE
#define NOTICE(A)
Definition JMessage.hh:64
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:69
WARNING
#define WARNING(A)
Definition JMessage.hh:65
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
gmake_property
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Definition JProperties.hh:1293
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition PhysicsEvent.hh:434
JDETECTOR::JDetector
Detector data structure.
Definition JDetector.hh:96
JDETECTOR::JPMTParametersMap
Auxiliary class for map of PMT parameters.
Definition JPMTParametersMap.hh:93
JDETECTOR::JPMTParameters
Data structure for PMT parameters.
Definition JPMTParameters.hh:29
JEEP::JProperties
Utility class to parse parameter values.
Definition JProperties.hh:501
JLANG::JException
General exception.
Definition JException.hh:24
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JSUPPORT::JSingleFileScanner
Object reading from a list of files.
Definition JSingleFileScanner.hh:172
JTOOLS::JRange
Range of values.
Definition JRange.hh:42
JTOOLS::JRange::getUpperLimit
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
JCALIBRATE::_2F
static const char *const _2F
Name extension for 2F rate fitted.
Definition JCalibrateK40.hh:37
JCALIBRATE::FIT_PMTS_QE_FIXED_t
@ FIT_PMTS_QE_FIXED_t
fit parameters of PMTs with QE fixed
Definition JFitK40.hh:56
JCALIBRATE::FIT_PMTS_AND_ANGULAR_DEPENDENCE_t
@ FIT_PMTS_AND_ANGULAR_DEPENDENCE_t
fit parameters of PMTs and angular dependence of K40 rate
Definition JFitK40.hh:54
JCALIBRATE::FIT_MODEL_t
@ FIT_MODEL_t
fit parameters of K40 rate and TTSs of PMTs
Definition JFitK40.hh:57
JCALIBRATE::FIT_PMTS_AND_BACKGROUND_t
@ FIT_PMTS_AND_BACKGROUND_t
fit parameters of PMTs and background
Definition JFitK40.hh:55
JCALIBRATE::FIT_PMTS_t
@ FIT_PMTS_t
fit parameters of PMTs
Definition JFitK40.hh:53
JCALIBRATE::_2R
static const char *const _2R
Name extension for 2D rate measured.
Definition JCalibrateK40.hh:36
JDETECTOR::getLabel
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
Definition JLocation.hh:247
JDETECTOR::getRangeOfFloors
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Definition JDetectorToolkit.hh:267
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:505
JDETECTOR::store
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
Definition JDetectorToolkit.hh:602
JDETECTOR::getSurvivalProbability
double getSurvivalProbability(const JPMTParameters &parameters)
Get model dependent probability that a one photo-electron hit survives the simulation of the PMT assu...
Definition JPMTParametersToolkit.hh:32
JFIT::model
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
Definition JGandalf.hh:56
JLANG::getIndex
int getIndex()
Get index for user I/O manipulation.
Definition JManip.hh:26
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:43
JROOT::putObject
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
Definition JRootFileWriter.hh:41
JTOOLS::result
return result
Definition JPolint.hh:853
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
KM3NETDAQ::NUMBER_OF_PMTS
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition JDAQ.hh:26
FILL
Auxiliary data structure for sequence of same character.
Definition JManip.hh:330
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JAANET::JRange_t
Type definition of range.
Definition JHead.hh:43
JAANET::K40
Livetime of noise data.
Definition JHead.hh:1062
JAANET::detector
Detector file.
Definition JHead.hh:227
JACOUSTICS::JFit
Acoustic single fit.
Definition JAcoustics/JEvt.hh:37
JACOUSTICS::JModel
Model for fit to acoustics data.
Definition JAcoustics/JModel.hh:480
JCALIBRATE::JK40Parameters
Fit parameters for two-fold coincidence rate due to K40.
Definition JFitK40.hh:610
JCALIBRATE::JK40Parameters::getInstance
static const JK40Parameters & getInstance()
Get default values.
Definition JFitK40.hh:627
JCALIBRATE::JTDC_t
Auxiliary class for TDC constraints.
Definition JTDC_t.hh:39
JCALIBRATE::JTDC_t::equal_range
range_type equal_range(const int id)
Get range of constraints for given module.
Definition JTDC_t.hh:101
JCALIBRATE::JTDC_t::reverse
void reverse()
Reverse constraints.
Definition JTDC_t.hh:137
JCALIBRATE::JTDC_t::is_valid
bool is_valid(const bool option=false) const
Check validity of TDC constrants.
Definition JTDC_t.hh:171
JCALIBRATE::data_type
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition JFitK40.hh:103
JCALIBRATE::rate_type
Data structure for measured coincidence rate of pair of PMTs.
Definition JFitK40.hh:66
JDETECTOR::JStringRouter
Router for mapping of string identifier to index.
Definition JStringRouter.hh:23
JEEP::JComment::add
JComment & add(const std::string &comment)
Add comment.
Definition JComment.hh:100
JLANG::JObjectStreamIO::store
void store(const char *file_name) const
Store to output file.
Definition JObjectStreamIO.hh:41
JLANG::JObjectStreamIO::load
void load(const char *file_name)
Load from input file.
Definition JObjectStreamIO.hh:30
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72
JTOOLS::JCombinatorics::pair_type
Data structure for a pair of indices.
Definition JCombinatorics.hh:29
Generated by doxygen 1.12.0