Jpp
JHobbit.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <map>
6 
7 #include "TROOT.h"
8 #include "TFile.h"
9 #include "TH1D.h"
10 #include "TH2D.h"
11 #include "TF1.h"
12 #include "TMath.h"
13 #include "TFitResult.h"
14 
15 #include "JLang/JLangToolkit.hh"
16 #include "JSupport/JMeta.hh"
17 
18 #include "JDetector/JDetector.hh"
19 #include "JDetector/JDetectorToolkit.hh"
20 #include "JDetector/JTimeRange.hh"
21 
22 #include "Jeep/JPrint.hh"
23 #include "Jeep/JParser.hh"
24 #include "Jeep/JMessage.hh"
25 
26 namespace {
27 
28  /**
29  * Name of histogram with fit results from JGandalf.
30  */
31  const char* h2_t = "h2";
32 
33  const char* gauss_t = "Gauss"; //!< Gaussian
34  const char* landau_t = "Landau"; //!< Landau
35  const char* emg_t = "EMG"; //!< Exponentially Modified Gaussian
36  const char* breitwigner_t = "BreitWigner"; //!< Breit-Wigner
37 }
38 
39 
40 /**
41  * \file
42  * Fit time-residuals histogram in output of JGandalf.cc in calibration mode.
43  * \author mdejong
44  */
45 int main(int argc, char **argv)
46 {
47  using namespace std;
48  using namespace JPP;
49 
50  vector<string> inputFile;
51  string outputFile;
52  string detectorFile;
53  bool overwriteDetector;
54  string function;
55  JTimeRange T_ns;
56  string option;
57  JTimeRange E_ns;
58  int debug;
59 
60  try {
61 
62  JParser<> zap("Program to fit time-residuals histogram in output of JGandalf.cc in calibration mode.");
63 
64  zap['f'] = make_field(inputFile, "input files (output from JGandalf -c).");
65  zap['o'] = make_field(outputFile, "output file.") = "hobbit.root";
66  zap['a'] = make_field(detectorFile, "detector file.");
67  zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with correct time offsets.");
68  zap['F'] = make_field(function, "fit function") = gauss_t, landau_t, emg_t, breitwigner_t;
69  zap['T'] = make_field(T_ns, "ROOT fit range around maximum [ns].") = JTimeRange(-7.5,+5.0);
70  zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "";
71  zap['E'] = make_field(E_ns, "validity range of t0 uncertainty [ns].") = JTimeRange( 0.0, 0.3);
72  zap['d'] = make_field(debug) = 1;
73 
74  zap(argc, argv);
75  }
76  catch(const exception &error) {
77  FATAL(error.what() << endl);
78  }
79 
80 
81  cout.tie(&cerr);
82 
83  if (!T_ns.is_valid()) {
84  FATAL("Invalid fit range [ns] " << T_ns << endl);
85  }
86 
87  JDetector detector;
88 
89  try {
90  load(detectorFile, detector);
91  }
92  catch(const JException& error) {
93  FATAL(error);
94  }
95 
96  if (option.find('S') == string::npos) { option += 'S'; }
97  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
98 
99 
100  TH2D* h2 = NULL;
101 
102  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
103 
104  NOTICE("Processing " << *i << endl);
105 
106  TFile in(i->c_str(), "exist");
107 
108  if (!in.IsOpen()) {
109  FATAL("File " << *i << " not opened." << endl);
110  }
111 
112  TH2D* p = dynamic_cast<TH2D*>(in.Get(h2_t));
113 
114  if (p == NULL) {
115  FATAL("File " << *i << " has no histogram <" << h2_t << ">." << endl);
116  }
117 
118  if (h2 == NULL)
119  h2 = (TH2D*) p->Clone();
120  else
121  h2->Add(p);
122 
123  h2->SetDirectory(0);
124 
125  in.Close();
126  }
127 
128  if (h2 == NULL) {
129  FATAL("No histogram <" << h2_t << ">." << endl);
130  }
131 
132  // auxiliary data structure for fit result
133 
134  struct result_type {
135 
136  result_type() :
137  value(0.0),
138  error(0.0)
139  {}
140 
141  result_type(double value,
142  double error) :
143  value(value),
144  error(error)
145  {}
146 
147  double value;
148  double error;
149  };
150 
151  typedef map<int, result_type> map_type; // module index -> fit result
152 
153  map_type zmap;
154 
155 
156  // histogram fit results
157 
158  const TAxis* x_axis = h2->GetXaxis();
159  const TAxis* y_axis = h2->GetYaxis();
160 
161  TH1D h0("h0", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
162  TH1D hc("hc", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
163  TH1D hq("hq", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
164 
165 
166  TFile out(outputFile.c_str(), "recreate");
167 
168  for (int ix = 1; ix <= x_axis->GetNbins(); ++ix) {
169 
170  TH1D h1(MAKE_CSTRING((detector.empty() ? ix : detector[ix - 1].getID()) << ".1D"), NULL, y_axis->GetNbins(), y_axis->GetXmin(), y_axis->GetXmax());
171 
172  // start values
173 
174  Double_t ymax = 0.0;
175  Double_t x0 = 0.0; // [ns]
176  Double_t sigma = 4.5; // [ns]
177 
178  for (int iy = 1; iy <= y_axis->GetNbins(); ++iy) {
179 
180  h1.SetBinContent(iy, h2->GetBinContent(ix,iy));
181  h1.SetBinError (iy, sqrt(h2->GetBinContent(ix,iy)));
182 
183  const Double_t x = h1.GetBinCenter (iy);
184  const Double_t y = h1.GetBinContent(iy);
185 
186  if (y > ymax) {
187  ymax = y;
188  x0 = x;
189  }
190  }
191 
192 
193  const Double_t xmin = x0 + T_ns.getLowerLimit();
194  const Double_t xmax = x0 + T_ns.getUpperLimit();
195 
196 
197  TF1* f1 = NULL;
198 
199  if (function == gauss_t) {
200 
201  f1 = new TF1(function.c_str(), "[0]*TMath::Gaus(x,[1],[2])");
202 
203  f1->SetParameter(0, ymax);
204  f1->SetParameter(1, x0);
205  f1->SetParameter(2, sigma);
206 
207  } else if (function == landau_t) {
208 
209  f1 = new TF1(function.c_str(), "[0]*TMath::Landau(x,[1],[2])");
210 
211  f1->SetParameter(0, ymax);
212  f1->SetParameter(1, x0);
213  f1->SetParameter(2, sigma);
214 
215  } else if (function == emg_t) {
216 
217  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]))");
218 
219  f1->SetParameter(0, ymax);
220  f1->SetParameter(1, x0 - sigma);
221  f1->SetParameter(2, sigma);
222  f1->SetParameter(3, 0.04);
223 
224  } else if (function == breitwigner_t) {
225 
226  f1 = new TF1(function.c_str(), "(x <= [1])*[0]*[2]*TMath::BreitWigner(x,[1],[2])+(x > [1])*[0]*[3]*TMath::BreitWigner(x,[1],[3])");
227 
228  f1->SetParameter(0, ymax);
229  f1->SetParameter(1, x0);
230  f1->SetParameter(2, 15.0);
231  f1->SetParameter(3, 25.0);
232 
233  } else {
234 
235  FATAL("Unknown fit function " << function << endl);
236  }
237 
238 
239  DEBUG("Start values:" << endl);
240 
241  for (int i = 0; i != f1->GetNpar(); ++i) {
242  DEBUG(left << setw(12) << f1->GetParName (i) << ' ' <<
243  SCIENTIFIC(12,5) << f1->GetParameter(i) << endl);
244  }
245 
246  // fit
247 
248  TFitResultPtr result = h1.Fit(f1, option.c_str(), "same", xmin, xmax);
249 
250  if (debug >= notice_t || !result->IsValid()) {
251  cout << "Slice: "
252  << setw(4) << ix << ' '
253  << setw(16) << h1.GetName() << ' '
254  << FIXED(7,3) << f1->GetParameter(1) << " +/- "
255  << FIXED(7,3) << f1->GetParError(1) << ' '
256  << FIXED(7,3) << result->Chi2() << '/'
257  << result->Ndf() << ' '
258  << (result->IsValid() ? "" : "failed") << endl;
259  }
260 
261  if (result->IsValid()) {
262  zmap[ix] = result_type(f1->GetParameter(1), f1->GetParError (1));
263  }
264 
265  if (result->Ndf() > 0) {
266  hc.SetBinContent(ix, result->Chi2() / result->Ndf());
267  }
268  hq.SetBinContent(ix, result->IsValid() ? 1.0 : 0.0);
269 
270  h1.Write();
271 
272  delete f1;
273  }
274 
275 
276  if (!zmap.empty()) {
277 
278  // average time offset
279 
280  double t0 = 0.0;
281 
282  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
283  t0 += i->second.value;
284  }
285 
286  t0 /= zmap.size();
287 
288  NOTICE("Average time offset [ns] " << FIXED(7,2) << t0 << endl);
289 
290  for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {
291  i->second.value -= t0;
292  }
293 
294  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
295  h0.SetBinContent(i->first, i->second.value);
296  h0.SetBinError (i->first, i->second.error);
297  }
298 
299  if (!detector.empty()) {
300 
301  const JRange<int> string(detector.begin(), detector.end(), &JModule::getString);
302  const JRange<int> floor (detector.begin(), detector.end(), &JModule::getFloor);
303 
304  TH2D hi("hi", NULL,
305  string.getLength() + 1,
306  string.getLowerLimit() - 0.5,
307  string.getUpperLimit() + 0.5,
308  floor.getLength() + 1,
309  floor.getLowerLimit() - 0.5,
310  floor.getUpperLimit() + 0.5);
311 
312  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
313  hi.SetBinContent(detector[i->first - 1].getString(),
314  detector[i->first - 1].getFloor(),
315  i->second.value);
316  }
317 
318  hi.Write();
319  }
320 
321 
322  if (overwriteDetector) {
323 
324  NOTICE("Store calibration data on file " << detectorFile << endl);
325 
326  detector.comment.add(JMeta(argc, argv));
327 
328  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
329 
330  if (E_ns(i->second.error))
331  detector[i->first - 1].sub(i->second.value);
332  else
333  ERROR("Slice " << setw(4) << i->first << " fit uncertainty " << FIXED(5,2) << i->second.error << " outside specified range (option -E <E_ns>)" << endl);
334  }
335 
336  store(detectorFile, detector);
337  }
338 
339  } else {
340 
341  NOTICE("No calibration results." << endl);
342  }
343 
344  h0 .Write();
345  hc .Write();
346  hq .Write();
347  h2->Write();
348 
349  out.Close();
350 }
JMeta.hh
FIXED
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:481
JMessage.hh
JPrint.hh
JEEP::notice_t
notice
Definition: JMessage.hh:32
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:456
std::vector
Definition: JSTDTypes.hh:12
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JTimeRange.hh
JTOOLS::JTimeRange
JRange< double > JTimeRange
Type definition for time range.
Definition: JTools/JTimeRange.hh:19
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:708
ERROR
#define ERROR(A)
Definition: JMessage.hh:66
debug
int debug
debug level
Definition: JSirene.cc:59
JTOOLS::result
return result
Definition: JPolint.hh:695
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:518
std::map
Definition: JSTDTypes.hh:16
JParser.hh
JDetectorToolkit.hh
JDETECTOR::store
void store(const JString &file_name, const JDetector &detector)
Store detector to output file.
Definition: JDetectorToolkit.hh:532
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std
Definition: jaanetDictionary.h:36
JDetector.hh
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JEEP::debug_t
debug
Definition: JMessage.hh:29
JLangToolkit.hh
main
int main(int argc, char **argv)
Definition: JHobbit.cc:45
Generated by   doxygen 1.8.16