Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
JFitToT.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4
5#include "TROOT.h"
6#include "TFile.h"
7#include "TH1D.h"
8#include "TH2D.h"
9#include "TF1.h"
10#include "TFitResult.h"
11#include "TMath.h"
12
15
16#include "JLang/JLangToolkit.hh"
18
23
24#include "JROOT/JRootToolkit.hh"
25
26#include "JTools/JRange.hh"
27
28#include "JSupport/JMeta.hh"
29#include "JSupport/JSupport.hh"
32
34#include "JCalibrate/JFitToT.hh"
35
36#include "Jeep/JPrint.hh"
37#include "Jeep/JParser.hh"
38#include "Jeep/JMessage.hh"
39
40namespace {
41
42 /**
43 * Wild card for histogram naming.\n
44 * The wild card will replaced by the module identifier.
45 */
46 static const char* const WILDCARD = "%";
47}
48
49
50/**
51 * \file
52 *
53 * Auxiliary program to fit time-over-threshold distributions.
54 * The fitted parameters are the PMT gain and gain spread.
55 * \author mkarel
56 */
57int main(int argc, char **argv)
58{
59 using namespace std;
60 using namespace JPP;
61 using namespace KM3NETDAQ;
62
65
66 string inputFile;
68 string detectorFile;
69 string pmtFile;
70 bool writeFits;
71
72 double Wmin = 5e3; // Minimal histogram weight in fit range
73 double LeftMargin = 10.0; // Fit range left of time-over-threshold peak
74 double RightMargin = 10.0; // Fit range right of time-over-threshold peak
75 double gradientThreshold = -0.005; // Minimal right-to-left gradient value in time-over-threshold peak-search,\n
76 // normalized with respect to the total number of histogram entries
77 JRange_t gainLimits (FITTOT_GAIN_MIN,
78 FITTOT_GAIN_MAX); // gain limits
79 JRange_t gainSpreadLimits(FITTOT_GAINSPREAD_MIN,
80 FITTOT_GAINSPREAD_MAX); // gainspread limits
81
82 JRange_t ToTfitRange;
83 string option;
84
85 string regexp;
86 int debug;
87
88 try {
89
90 JProperties properties;
91
92 properties.insert(gmake_property(Wmin));
93 properties.insert(gmake_property(LeftMargin));
94 properties.insert(gmake_property(RightMargin));
95 properties.insert(gmake_property(gradientThreshold));
96 properties.insert(gmake_property(gainLimits));
97 properties.insert(gmake_property(gainSpreadLimits));
98
99 JParser<> zap("Auxiliary program to fit time-over-threshold distributions.");
100
101 zap['f'] = make_field(inputFile, "input file (output from JCalibrateToT).");
102 zap['o'] = make_field(outputFile, "output file.") = "fit.root";
103 zap['a'] = make_field(detectorFile, "detector file.");
104 zap['P'] = make_field(pmtFile, "specify PMT file name that can be given to JTriggerEfficiency.") = "";
105 zap['w'] = make_field(writeFits, "write fit results.");
106 zap['O'] = make_field(option, "ROOT fit options, see TH1::Fit") = "";
107 zap['@'] = make_field(properties, "fit properties") = JPARSER::initialised();
108 zap['x'] = make_field(ToTfitRange, "ROOT fit range (time-over-threshold).") = JRange_t(0.0, 35.0);
109 zap['R'] = make_field(regexp, "regular expression for histogram name.") = MAKE_CSTRING(WILDCARD << _2SToT);
110 zap['d'] = make_field(debug, "debug.") = 1;
111
112 zap(argc, argv);
113 }
114 catch(const exception &error) {
115 FATAL(error.what() << endl);
116 }
117
118
119 //----------------------------------------------------------
120 // load detector file and PMT parameters
121 //----------------------------------------------------------
122
124
125 try {
126 load(detectorFile, detector);
127 }
128 catch(const JException& error) {
129 FATAL(error);
130 }
131
132
133 JPMTParametersMap parametersMap;
134
135 if (!pmtFile.empty()) {
136
137 try {
138 parametersMap.load(pmtFile.c_str());
139 }
140 catch(const JException& error) {}
141 }
142
143 // Set ROOT fit options
144 if (option.find('R') == string::npos) { option += 'R'; }
145 if (option.find('S') == string::npos) { option += 'S'; }
146 if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
147
148
149 TFile* in = TFile::Open(inputFile.c_str(), "exist");
150
151 if (in == NULL || !in->IsOpen()) {
152 FATAL("File: " << inputFile << " not opened." << endl);
153 }
154
155
156 TH2D gain("gain", NULL,
157 100, 0.0, 2.0,
158 100, 0.0, 1.0);
159 TH1D chi2("chi2", NULL, 100, 0.0, 5.0);
160
161 outputFile.open();
162
163 if (!outputFile.is_open()) {
164 FATAL("Error opening file " << outputFile << endl);
165 }
166
167 outputFile.put(JMeta(argc, argv));
168
169 //----------------------------------------------------------
170 // loop over modules in detector file to fit histogram
171 //----------------------------------------------------------
172
173 for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
174
175 NOTICE("Module " << setw(10) << module->getID() << ' ' << module->getLocation() << endl);
176
177 if (module->empty()) {
178 continue;
179 }
180
181 //----------------------------------------------------------
182 // get histogram for this module
183 //----------------------------------------------------------
184
185 TH2D* h2s = (TH2D*) in->Get(replace(regexp, WILDCARD, MAKE_STRING(module->getID())).c_str());
186
187 if (h2s == NULL) {
188
189 WARNING("No histogram " << module->getID() << "; skip fit." << endl);
190
191 continue;
192 }
193
194 const double ny = h2s->GetYaxis()->GetNbins();
195 const double ymin = h2s->GetYaxis()->GetXmin();
196 const double ymax = h2s->GetYaxis()->GetXmax();
197
198 TH1D chi2_1d (MAKE_CSTRING(module->getID() << ".1chi2"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
199 TH1D gain_1d (MAKE_CSTRING(module->getID() << ".1gain"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
200 TH1D gainspread_1d(MAKE_CSTRING(module->getID() << ".1gainspread"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
201
202 for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
203
204 const JPMTIdentifier id(module->getID(), ix-1);
205
206 const string name = MAKE_STRING(module->getID() << '.' << FILL(2,'0') <<
207 id.getPMTAddress() << FITTOT_SUFFIX);
208
209 DEBUG("Processing histogram " << name << endl);
210
211 TH1D h1(name.c_str(), NULL, ny, ymin, ymax);
212
213 h1.Sumw2();
214
215 for (int iy = h2s->GetNbinsY(); iy >= 1; --iy) {
216
217 const double w = h1.GetBinWidth (iy);
218 const double y = h2s->GetBinContent(ix, iy);
219
220 h1.SetBinContent(iy, y/w);
221 h1.SetBinError (iy, sqrt(y)/w);
222 }
223
224 JRange_t range(ToTfitRange);
225
226 const double weight = h1.Integral(h1.FindBin(range.getLowerLimit()),
227 h1.FindBin(range.getUpperLimit()));
228
229 if (weight <= Wmin) {
230
231 WARNING("Insufficient histogram entries for PMT " << id <<
232 "(" << weight << " < " << Wmin << "); skip fit." << endl);
233
234 if (writeFits) {
235 outputFile.put(h1);
236 }
237
238 continue;
239 }
240
241 //----------------------------------------------------------
242 // find the mode of the time-over-threshold distribution
243 //----------------------------------------------------------
244
245 double max = 0.0;
246 double ToTmode = JDETECTOR::TIME_OVER_THRESHOLD_NS;
247
248 JPMTParameters& parameters = parametersMap[id];
249
250 DEBUG(RIGHT(10) << "ToT" << RIGHT(10) << "Counts" << RIGHT(10) << "gradient" <<
251 RIGHT(10) << "threshold" << RIGHT(10) << "Mode" << RIGHT(10) << "Max" << endl);
252
253 for (int i = h1.GetBinCenter(ny-1);
254 i > h1.GetBinCenter(h1.FindBin(parameters.mean_ns + parameters.sigma_ns));
255 --i) {
256
257 const double x = h1.GetBinCenter (i);
258 const double y = h1.GetBinContent(i);
259
260 const double gradient = ( (h1.GetBinContent(i-1) - h1.GetBinContent(i+1)) /
261 (h1.GetBinCenter (i+1) - h1.GetBinCenter (i-1)) );
262
263 DEBUG(FIXED(10,1) << x << FIXED(10,1) << y << FIXED(10,1) << gradient <<
264 FIXED(10,1) << -fabs(gradientThreshold) * h1.Integral() << FIXED(10,1) << ToTmode << FIXED(10,1) << max << endl);
265
266
267 if (y > max) {
268
269 ToTmode = x;
270 max = y;
271
272 } else if (gradient < -fabs(gradientThreshold) * h1.Integral()) {
273
274 break;
275 }
276 }
277
278 //----------------------------------------------------------
279 // set fit range
280 //----------------------------------------------------------
281
282 if (!range.is_valid()) {
283 range.setRange(ToTmode - LeftMargin,
284 ToTmode + RightMargin);
285 }
286
287 range.setRange(std::max(h1.GetXaxis()->GetXmin(), range.getLowerLimit()),
288 std::min(h1.GetXaxis()->GetXmax(), range.getUpperLimit()));
289
290 DEBUG(LEFT(10) << "Fit-range: " << FIXED(20,1) << range << endl);
291
292 //----------------------------------------------------------
293 // perform minimum chi-square fit
294 //----------------------------------------------------------
295
296 JFitToT fit(parameters, range);
297
298 setParLimits(fit, fit.getModelParameter(&JFitToT::JFitToTParameters::gain),
299 gainLimits.getLowerLimit(), gainLimits.getUpperLimit());
300 setParLimits(fit, fit.getModelParameter(&JFitToT::JFitToTParameters::gainSpread),
301 gainSpreadLimits.getLowerLimit(), gainSpreadLimits.getUpperLimit());
302
303 const double initGain = fit.getCPU().getNPE(ToTmode);
304 const double initGainSpread = initGain / 3.0;
305
306 if (!gainLimits(initGain)) { // Invalid initial gain value
307
308 WARNING("Invalid initial gain for PMT " << id << "; set default values." << endl);
309
310 parameters.gain = gainLimits .constrain(initGain);
311 parameters.gainSpread = gainSpreadLimits.constrain(initGainSpread);
312
313 if (writeFits) {
314 outputFile.put(h1);
315 }
316
317 continue;
318 }
319
320 fit.gain = initGain;
321 fit.gainSpread = initGainSpread;
322
323 const TFitResultPtr result = fit(h1, option);
324
325 if (result.Get() == NULL) {
326 FATAL("Invalid TFitResultPtr " << h1.GetName() << endl);
327 }
328
329 if (int(result) < 0 || !result->IsValid()) { // Fit failed
330
331 WARNING("Fit failure for PMT " << id << "; set initialization values." << endl);
332
333 h1.GetListOfFunctions()->Clear();
334
335 parameters.gain = initGain;
336 parameters.gainSpread = initGainSpread;
337
338 if (writeFits) {
339 outputFile.put(h1);
340 }
341
342 continue;
343 }
344
345 if ((ymin < range.getLowerLimit() || ymax > range.getUpperLimit()) &&
346 (option.find("0") == string::npos && option.find("N") == string::npos)) {
347
348 // add function with full range
349
350 TF1* f1 = new TF1(MAKE_CSTRING(FITTOT_FNAME << "_fullRange"),
351 &fit,
353 ymin,
354 round(parameters.saturation) - 0.5,
355 JFitToT::getNumberOfModelParameters());
356 f1->SetParameters(fit.GetParameters());
357 f1->SetLineStyle(kDotted);
358 f1->SetNpx(1000);
359
360 h1.GetListOfFunctions()->Add(f1);
361 }
362
363 // store fit results
364
365 const double reducedChi2 = result->Chi2() / result->Ndf();
366
367 const JFitToTParameters values(fit.GetParameters());
368 const JFitToTParameters errors(fit.GetParErrors());
369
370 chi2_1d. SetBinContent(ix, reducedChi2);
371 gain_1d. SetBinContent(ix, values.gain);
372 gain_1d. SetBinError (ix, errors.gain);
373 gainspread_1d.SetBinContent(ix, values.gainSpread);
374 gainspread_1d.SetBinError (ix, errors.gainSpread);
375
376 gain.Fill(fit.gain, fit.gainSpread);
377 chi2.Fill(TMath::Log10(reducedChi2));
378
379 NOTICE("PMT " << id << " chi2 / NDF " << result->Chi2() << ' ' << result->Ndf() << endl);
380
381 parameters.gain = fit.gain;
382 parameters.gainSpread = fit.gainSpread;
383
384 if (writeFits) {
385 outputFile.put(h1);
386 }
387 }
388
389 if (writeFits) {
390 outputFile.put(chi2_1d);
391 outputFile.put(gain_1d);
392 outputFile.put(gainspread_1d);
393 }
394 }
395
396 {
397 vector<JMeta> buffer(1, JMeta(argc, argv));
398
399 for (JSingleFileScanner<JMeta> in(inputFile); in.hasNext(); ) {
400
401 const JMeta* meta = in.next();
402
403 buffer.push_back(*meta);
404
405 outputFile.put(*meta);
406 }
407
408 for (vector<JMeta>::const_reverse_iterator i = buffer.rbegin(); i != buffer.rend(); ++i) {
409 parametersMap.comment.add(*i);
410 }
411 }
412
413 for (JRootFileReader<JDAQHeader> in(inputFile.c_str()); in.hasNext(); ) {
414 outputFile.put(*in.next());
415 }
416
417
418 if (pmtFile != "") {
419 parametersMap.store(pmtFile.c_str());
420 }
421
422 outputFile.put(gain);
423 outputFile.put(chi2);
424
425 outputFile.close();
426}
427
428
string outputFile
KM3NeT DAQ constants, bit handling, etc.
Data structure for detector geometry and calibration.
Recording of objects on file according a format that follows from the file name extension.
int main(int argc, char **argv)
Definition JFitToT.cc:57
General purpose messaging.
#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
ROOT I/O of application specific meta data.
PMT analogue signal processor.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Physics constants.
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
#define MAKE_STRING(A)
Make string.
Definition JPrint.hh:63
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Auxiliary class to define a range between two values.
Scanning of objects from a single file according a format that follows from the extension of each fil...
ROOT TTree parameter settings of various packages.
Detector data structure.
Definition JDetector.hh:96
Auxiliary class for map of PMT parameters.
Data structure for PMT parameters.
double sigma_ns
time-over-threshold standard deviation of threshold-band hits [ns]
double gainSpread
gain spread [unit]
double mean_ns
mean time-over-threshold of threshold-band hits [ns]
double saturation
saturation [ns]
Utility class to parse parameter values.
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
Object writing to file.
Object reading from a list of files.
virtual bool hasNext() override
Check availability of next element.
Range of values.
Definition JRange.hh:42
void setRange(const range_type &range)
Set range.
Definition JRange.hh:146
bool is_valid() const
Check validity of range.
Definition JRange.hh:311
T constrain(argument_type x) const
Constrain value to range.
Definition JRange.hh:350
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
static const char *const _2SToT
Histogram naming.
static const double FITTOT_GAINSPREAD_MAX
Default maximal gain spread.
Definition JFitToT.hh:48
static const std::string FITTOT_FNAME
Definition JFitToT.hh:38
static const std::string FITTOT_SUFFIX
Definition JFitToT.hh:37
static const double FITTOT_GAIN_MAX
Default maximal gain.
Definition JFitToT.hh:45
static const double FITTOT_GAIN_MIN
Default minimal gain.
Definition JFitToT.hh:44
static const double FITTOT_GAINSPREAD_MIN
Default minimal gain spread.
Definition JFitToT.hh:47
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
const double TIME_OVER_THRESHOLD_NS
Specification for time-over-threshold corresponding to a one photo-electron pulse.
@ debug_t
debug
Definition JMessage.hh:29
std::string replace(const std::string &input, const std::string &target, const std::string &replacement)
Replace tokens in string.
@ LEFT
Definition JTwosome.hh:18
@ RIGHT
Definition JTwosome.hh:18
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool setParLimits(TF1 &f1, const Int_t index, Double_t xmin, Double_t xmax)
Set fit parameter limits.
return result
Definition JPolint.hh:862
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
static const char WILDCARD
Definition JDAQTags.hh:56
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
Detector file.
Definition JHead.hh:227
Fit parameters for two-fold coincidence rate due to K40.
Definition JFitToT.hh:54
Double_t gainSpread
PMT gain spread.
Definition JFitToT.hh:170
const Double_t getModelParameter(const int i) const
Get model parameter.
Definition JFitToT.hh:147
Double_t gain
PMT gain.
Definition JFitToT.hh:169
Parametrisation of time-over-threshold distribution.
Definition JFitToT.hh:183
Double_t getValue(const Double_t *x, const Double_t *par)
Get rate as a function of the fit parameters.
Definition JFitToT.hh:270
const JPMTAnalogueSignalProcessor & getCPU() const
Access method for the analogue signal processor.
Definition JFitToT.hh:211
virtual double getNPE(const double tot_ns) const override
Get number of photo-electrons.
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.
Type list.
Definition JTypeList.hh:23
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