Jpp  master_rocky-43-ge265d140c
the software that should make you happy
Functions
JHobbit.cc File Reference

Program to fit time-residuals histogram in output of JCalibrateMuon.cc. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TMath.h"
#include "TFitResult.h"
#include "JROOT/JMinimizer.hh"
#include "JLang/JLangToolkit.hh"
#include "JSupport/JMeta.hh"
#include "JTools/JRange.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTimeRange.hh"
#include "JDetector/JStringRouter.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

Program to fit time-residuals histogram in output of JCalibrateMuon.cc.

Author
mdejong

Definition in file JHobbit.cc.

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 49 of file JHobbit.cc.

50 {
51  using namespace std;
52  using namespace JPP;
53 
54  vector<string> inputFile;
55  string outputFile;
56  string detectorFile;
57  bool overwriteDetector;
58  string function;
59  JTimeRange T_ns;
60  string option;
61  JTimeRange E_ns;
62  double WMin;
63  int ID;
64  int debug;
65 
66  try {
67 
68  JParser<> zap("Program to fit time-residuals histogram in output of JCalibrateMuon.cc.");
69 
70  zap['f'] = make_field(inputFile, "input files (output from JCalibrateMuon).");
71  zap['o'] = make_field(outputFile, "output file.") = "hobbit.root";
72  zap['a'] = make_field(detectorFile, "detector file.");
73  zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with correct time offsets.");
74  zap['F'] = make_field(function, "fit function") = gauss_t, landau_t, emg_t, breitwigner_t;
75  zap['T'] = make_field(T_ns, "ROOT fit range around maximum [ns].") = JTimeRange(-7.5,+5.0);
76  zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "";
77  zap['E'] = make_field(E_ns, "validity range of t0 uncertainty [ns].") = JTimeRange( 0.0, 0.5);
78  zap['W'] = make_field(WMin, "minimal histogram contents.") = 100.0;
79  zap['M'] = make_field(ID, "module identifier") = -1;
80  zap['d'] = make_field(debug) = 1;
81 
82  zap(argc, argv);
83  }
84  catch(const exception &error) {
85  FATAL(error.what() << endl);
86  }
87 
88 
89  if (!T_ns.is_valid()) {
90  FATAL("Invalid fit range [ns] " << T_ns << endl);
91  }
92 
93  JDetector detector;
94 
95  try {
96  load(detectorFile, detector);
97  }
98  catch(const JException& error) {
99  FATAL(error);
100  }
101 
102 
103  if (option.find('S') == string::npos) { option += 'S'; }
104  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
105 
106  TH2D* h2 = NULL;
107 
108  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
109 
110  NOTICE("Processing " << *i << endl);
111 
112  TFile in(i->c_str(), "exist");
113 
114  if (!in.IsOpen()) {
115  FATAL("File " << *i << " not opened." << endl);
116  }
117 
118  TH2D* p = dynamic_cast<TH2D*>(in.Get(h2_t));
119 
120  if (p == NULL) {
121  FATAL("File " << *i << " has no histogram <" << h2_t << ">." << endl);
122  }
123 
124  if (h2 == NULL)
125  h2 = (TH2D*) p->Clone();
126  else
127  h2->Add(p);
128 
129  h2->SetDirectory(0);
130 
131  in.Close();
132  }
133 
134  if (h2 == NULL) {
135  FATAL("No histogram <" << h2_t << ">." << endl);
136  }
137 
138  // auxiliary data structure for fit result
139 
140  struct result_type {
141 
142  result_type() :
143  value(0.0),
144  error(0.0)
145  {}
146 
147  result_type(double value,
148  double error) :
149  value(value),
150  error(error)
151  {}
152 
153  double value;
154  double error;
155  };
156 
157  typedef map<int, result_type> map_type; // module index -> fit result
158 
159  map_type zmap;
160 
161 
162  // histogram fit results
163 
164  const TAxis* x_axis = h2->GetXaxis();
165  const TAxis* y_axis = h2->GetYaxis();
166 
167  TH1D h0("h0", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
168  TH1D hc("hc", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
169  TH1D hq("hq", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
170 
171 
172  TFile out(outputFile.c_str(), "recreate");
173 
174  size_t counts = 0;
175  size_t errors = 0;
176 
177  for (Int_t ix = 1; ix <= x_axis->GetNbins(); ++ix) {
178 
179  const JModule& module = detector[ix - 1];
180 
181  DEBUG("Module " << setw(8) << module.getID() << ' ' << getLabel(module.getLocation()) << endl);
182 
183  if (module.getFloor() == 0) {
184  continue;
185  }
186 
187  if (ID != -1 && ID != module.getID()) {
188  continue;
189  }
190 
191  TH1D h1(MAKE_CSTRING(module.getID() << ".1D"), NULL, y_axis->GetNbins(), y_axis->GetXmin(), y_axis->GetXmax());
192 
193  // start values
194 
195  Double_t ymax = 0.0;
196  Double_t x0 = 0.0; // [ns]
197  Double_t sigma = 4.0; // [ns]
198  Double_t W = 0.0;
199 
200  for (Int_t iy = 1; iy <= y_axis->GetNbins(); ++iy) {
201 
202  h1.SetBinContent(iy, h2->GetBinContent(ix,iy));
203  h1.SetBinError (iy, sqrt(h2->GetBinContent(ix,iy)));
204 
205  const Double_t x = h1.GetBinCenter (iy);
206  const Double_t y = h1.GetBinContent(iy);
207 
208  W += y;
209 
210  if (y > ymax) {
211  ymax = y;
212  x0 = x;
213  }
214  }
215 
216  if (W <= WMin) {
217 
218  WARNING("Not enough entries for slice " << ix << ' ' << W << "; skip fit." << endl);
219 
220  continue;
221  }
222 
223  ymax *= 0.9;
224 
225  const Double_t xmin = x0 + T_ns.getLowerLimit();
226  const Double_t xmax = x0 + T_ns.getUpperLimit();
227 
228 
229  TF1* f1 = NULL;
230 
231  if (function == gauss_t) {
232 
233  f1 = new TF1(function.c_str(), "[0]*TMath::Gaus(x,[1],[2])");
234 
235  f1->SetParameter(0, 0.8*ymax);
236  f1->SetParameter(1, x0);
237  f1->SetParameter(2, sigma);
238 
239  f1->SetParError(0, sqrt(ymax) * 0.1);
240  f1->SetParError(1, 0.01);
241  f1->SetParError(2, 0.01);
242 
243  } else if (function == landau_t) {
244 
245  f1 = new TF1(function.c_str(), "[0]*TMath::Landau(x,[1],[2])");
246 
247  f1->SetParameter(0, 0.8*ymax);
248  f1->SetParameter(1, x0);
249  f1->SetParameter(2, sigma);
250 
251  f1->SetParError(0, sqrt(ymax) * 0.1);
252  f1->SetParError(1, 0.01);
253  f1->SetParError(2, 0.01);
254 
255  } else if (function == emg_t) {
256 
257  f1 = new TF1(function.c_str(), "[0]*TMath::Exp(0.5*[3]*(2.0*[1]+[3]*[2]*[2]-2.0*x))*TMath::Erfc(([1]+[3]*[2]*[2]-x)/(TMath::Sqrt(2.0)*[2]))");
258 
259  f1->SetParameter(0, 0.5*ymax);
260  f1->SetParameter(1, x0 - sigma);
261  f1->SetParameter(2, sigma);
262  f1->SetParameter(3, 0.06);
263 
264  f1->SetParError(0, sqrt(ymax) * 0.1);
265  f1->SetParError(1, 0.01);
266  f1->SetParError(2, 0.01);
267  f1->SetParError(3, 1.0e-4);
268 
269  } else if (function == breitwigner_t) {
270 
271  f1 = new TF1(function.c_str(), "(x <= [1])*[0]*[2]*TMath::BreitWigner(x,[1],[2])+(x > [1])*[0]*[3]*TMath::BreitWigner(x,[1],[3])");
272 
273  f1->SetParameter(0, 0.8*ymax);
274  f1->SetParameter(1, x0);
275  f1->SetParameter(2, 15.0);
276  f1->SetParameter(3, 25.0);
277 
278  f1->SetParError(0, sqrt(ymax) * 0.1);
279  f1->SetParError(1, 0.01);
280  f1->SetParError(2, 0.1);
281  f1->SetParError(3, 0.1);
282 
283  } else {
284 
285  FATAL("Unknown fit function " << function << endl);
286  }
287 
288 
289  DEBUG("Start values:" << endl);
290 
291  for (int i = 0; i != f1->GetNpar(); ++i) {
292  DEBUG(left << setw(12) << f1->GetParName (i) << ' ' <<
293  SCIENTIFIC(12,5) << f1->GetParameter(i) << endl);
294  }
295 
296  // fit
297 
298  TFitResultPtr result = h1.Fit(f1, option.c_str(), "same", xmin, xmax);
299 
300  if (debug >= notice_t || !result->IsValid()) {
301  cout << "Slice: "
302  << setw(4) << ix << ' '
303  << setw(16) << h1.GetName() << ' '
304  << FIXED(7,3) << f1->GetParameter(1) << " +/- "
305  << FIXED(7,3) << f1->GetParError(1) << ' '
306  << FIXED(7,3) << result->Chi2() << '/'
307  << result->Ndf() << ' '
308  << (result->IsValid() ? "" : "failed") << endl;
309  }
310 
311  counts += 1;
312 
313  if (!result->IsValid()) {
314  errors += 1;
315  }
316 
317  if (result->IsValid()) {
318  zmap[ix] = result_type(f1->GetParameter(1), f1->GetParError (1));
319  }
320 
321  if (result->Ndf() > 0) {
322  hc.SetBinContent(ix, result->Chi2() / result->Ndf());
323  }
324 
325  hq.SetBinContent(ix, result->IsValid() ? 1.0 : 0.0);
326 
327  h1.Write();
328 
329  delete f1;
330  }
331 
332 
333  if (!zmap.empty()) {
334 
335  // average time offset
336 
337  double t0 = 0.0;
338 
339  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
340  t0 += i->second.value;
341  }
342 
343  t0 /= zmap.size();
344 
345  NOTICE("Average time offset [ns] " << FIXED(7,2) << t0 << endl);
346  NOTICE("Number of fits passed/failed " << counts << "/" << errors << endl);
347 
348  for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {
349  i->second.value -= t0;
350  }
351 
352  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
353  h0.SetBinContent(i->first, i->second.value);
354  h0.SetBinError (i->first, i->second.error);
355  }
356 
357 
358  const floor_range range = getRangeOfFloors(detector);
359  const JStringRouter string(detector);
360 
361  TH2D hi("hi", NULL,
362  string.size(), -0.5, string.size() - 0.5,
363  range.getUpperLimit(), 1 - 0.5, range.getUpperLimit() + 0.5);
364 
365  for (Int_t i = 1; i <= hi.GetXaxis()->GetNbins(); ++i) {
366  hi.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1)));
367  }
368  for (Int_t i = 1; i <= hi.GetYaxis()->GetNbins(); ++i) {
369  hi.GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i));
370  }
371 
372  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
373  hi.SetBinContent(detector[i->first - 1].getString(),
374  detector[i->first - 1].getFloor(),
375  i->second.value);
376  }
377 
378  hi.Write();
379 
380 
381  if (overwriteDetector) {
382 
383  NOTICE("Store calibration data on file " << detectorFile << endl);
384 
385  detector.comment.add(JMeta(argc, argv));
386 
387  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
388 
389  if (E_ns(i->second.error))
390  detector[i->first - 1].sub(i->second.value);
391  else
392  ERROR("Slice " << setw(4) << i->first << " fit uncertainty " << FIXED(5,2) << i->second.error << " outside specified range (option -E <E_ns>)" << endl);
393  }
394 
395  store(detectorFile, detector);
396  }
397 
398  } else {
399 
400  NOTICE("No calibration results." << endl);
401  }
402 
403  h0 .Write();
404  hc .Write();
405  hq .Write();
406  h2->Write();
407 
408  out.Close();
409 }
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
ERROR
#define ERROR(A)
Definition: JMessage.hh:66
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
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:96
JDETECTOR::JLocation::getFloor
int getFloor() const
Get floor number.
Definition: JLocation.hh:146
JDETECTOR::JLocation::getLocation
const JLocation & getLocation() const
Get location.
Definition: JLocation.hh:70
JDETECTOR::JModule
Data structure for a composite optical module.
Definition: JModule.hh:75
JLANG::JException
General exception.
Definition: JException.hh:24
JLANG::JObjectID::getID
int getID() const
Get identifier.
Definition: JObjectID.hh:50
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1698
JTOOLS::JRange::is_valid
bool is_valid() const
Check validity of range.
Definition: JRange.hh:311
JTOOLS::JRange::getLowerLimit
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
JTOOLS::JRange::getUpperLimit
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
GAUSS_HERMITE::sigma
const double sigma[]
Definition: JQuadrature.cc:74
GAUSS_LEGENDRE::f1
const JPolynome f1(1.0, 2.0, 3.0)
Function.
GAUSS_LEGENDRE::xmax
const double xmax
Definition: JQuadrature.cc:24
GAUSS_LEGENDRE::xmin
const double xmin
Definition: JQuadrature.cc:23
JACOUSTICS::JTimeRange
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [s]).
Definition: JAcoustics/JTimeRange.hh:19
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
JEEP::debug_t
@ debug_t
debug
Definition: JMessage.hh:29
JEEP::notice_t
@ notice_t
notice
Definition: JMessage.hh:32
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
KM3NETDAQ::map_type
std::map< int, range_type > map_type
Definition: JDQSimulator.cc:161
std
Definition: JSTDTypes.hh:14
FIXED
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
JAANET::detector
Detector file.
Definition: JHead.hh:227
JDETECTOR::JStringRouter
Router for mapping of string identifier to index.
Definition: JStringRouter.hh:23
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:72
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:488
Generated by doxygen 1.9.1