Jpp  18.0.0-rc.4
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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"
17 #include "JPhysics/JConstants.hh"
18 
19 #include "JDetector/JDetector.hh"
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 
40 namespace {
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  */
57 int main(int argc, char **argv)
58 {
59  using namespace std;
60  using namespace JPP;
61  using namespace KM3NETDAQ;
62 
63  typedef JRange<double> JRange_t;
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 
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  ymax,
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 
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1514
General exception.
Definition: JException.hh:23
double gain
gain [unit]
data_type w[N+1][M+1]
Definition: JPolint.hh:778
const Double_t getModelParameter(const int i) const
Get model parameter.
Definition: JFitToT.hh:146
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
debug
Definition: JMessage.hh:29
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
TString replace(const TString &target, const TRegexp &regexp, const T &replacement)
Replace regular expression in input by given replacement.
Definition: JPrintResult.cc:63
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
double gainSpread
gain spread [unit]
static const double FITTOT_GAIN_MAX
Default maximal gain.
Definition: JFitToT.hh:44
Detector data structure.
Definition: JDetector.hh:89
Double_t gainSpread
PMT gain spread.
Definition: JFitToT.hh:169
const double TIME_OVER_THRESHOLD_NS
Specification for time-over-threshold corresponding to a one photo-electron pulse.
static const double FITTOT_GAIN_MIN
Default minimal gain.
Definition: JFitToT.hh:43
Recording of objects on file according a format that follows from the file name extension.
then echo Enter input within $TIMEOUT_S seconds echo n User name
Definition: JCookie.sh:42
Utility class to parse parameter values.
Definition: JProperties.hh:496
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:136
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
string outputFile
Parametrisation of time-over-threshold distribution.
Definition: JFitToT.hh:179
Data structure for detector geometry and calibration.
Fit parameters for two-fold coincidence rate due to K40.
Definition: JFitToT.hh:53
#define MAKE_STRING(A)
Make string.
Definition: JPrint.hh:127
const JPMTAnalogueSignalProcessor & getCPU() const
Access method for the analogue signal processor.
Definition: JFitToT.hh:210
static const std::string FITTOT_SUFFIX
Definition: JFitToT.hh:36
Type list.
Definition: JTypeList.hh:22
Scanning of objects from a single file according a format that follows from the extension of each fil...
const JPolynome f1(1.0, 2.0, 3.0)
Function.
Type definition of range.
Definition: JHead.hh:41
I/O formatting auxiliaries.
double mean_ns
mean time-over-threshold of threshold-band hits [ns]
Detector file.
Definition: JHead.hh:226
double sigma_ns
time-over-threshold standard deviation of threshold-band hits [ns]
static const double FITTOT_GAINSPREAD_MAX
Default maximal gain spread.
Definition: JFitToT.hh:47
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
void setRange(const range_type &range)
Set range.
Definition: JRange.hh:146
ROOT I/O of application specific meta data.
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
#define NOTICE(A)
Definition: JMessage.hh:64
virtual double getNPE(const double tot_ns) const override
Get number of photo-electrons.
Physics constants.
Auxiliary class for map of PMT parameters.
T constrain(argument_type x) const
Constrain value to range.
Definition: JRange.hh:350
Double_t gain
PMT gain.
Definition: JFitToT.hh:168
bool setParLimits(TF1 &f1, const Int_t index, Double_t xmin, Double_t xmax)
Set fit parameter limits.
General purpose messaging.
Auxiliary data structure for sequence of same character.
Definition: JManip.hh:328
#define FATAL(A)
Definition: JMessage.hh:67
void load(const char *file_name)
Load from input file.
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
then $JPP_DIR examples JDetector JToT o $OUTPUT_FILE n N $NPE P gain
Definition: JToT.sh:47
bool is_valid() const
Check validity of range.
Definition: JRange.hh:311
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
Utility class to parse command line options.
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
Definition: JMatrixNZ.sh:106
static const char *const _2SToT
Histogram naming.
PMT analogue signal processor.
static const double FITTOT_GAINSPREAD_MIN
Default minimal gain spread.
Definition: JFitToT.hh:46
Object reading from a list of files.
possible values
Data structure for PMT parameters.
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:46
then echo WARNING
Definition: JTuneHV.sh:91
ROOT file reader.
static const std::string FITTOT_FNAME
Definition: JFitToT.hh:37
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62