Jpp  15.0.1-rc.1-highQE
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 #include "JTools/JRange.hh"
18 
19 #include "JDetector/JDetector.hh"
21 #include "JDetector/JTimeRange.hh"
22 
23 #include "Jeep/JPrint.hh"
24 #include "Jeep/JParser.hh"
25 #include "Jeep/JMessage.hh"
26 
27 namespace {
28 
29  /**
30  * Name of histogram with fit results from JCalibrateMuon.cc.
31  */
32  const char* h2_t = "h2";
33 
34  const char* gauss_t = "Gauss"; //!< Gaussian
35  const char* landau_t = "Landau"; //!< Landau
36  const char* emg_t = "EMG"; //!< Exponentially Modified Gaussian
37  const char* breitwigner_t = "BreitWigner"; //!< Breit-Wigner
38 }
39 
40 
41 /**
42  * \file
43  * Program to fit time-residuals histogram in output of JCalibrateMuon.cc.
44  * \author mdejong
45  */
46 int main(int argc, char **argv)
47 {
48  using namespace std;
49  using namespace JPP;
50 
51  vector<string> inputFile;
52  string outputFile;
53  string detectorFile;
54  bool overwriteDetector;
55  string function;
56  JTimeRange T_ns;
57  string option;
58  JTimeRange E_ns;
59  double WMin;
60  int debug;
61 
62  try {
63 
64  JParser<> zap("Program to fit time-residuals histogram in output of JCalibrateMuon.cc.");
65 
66  zap['f'] = make_field(inputFile, "input files (output from JCalibrateMuon).");
67  zap['o'] = make_field(outputFile, "output file.") = "hobbit.root";
68  zap['a'] = make_field(detectorFile, "detector file.");
69  zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with correct time offsets.");
70  zap['F'] = make_field(function, "fit function") = gauss_t, landau_t, emg_t, breitwigner_t;
71  zap['T'] = make_field(T_ns, "ROOT fit range around maximum [ns].") = JTimeRange(-7.5,+5.0);
72  zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "";
73  zap['E'] = make_field(E_ns, "validity range of t0 uncertainty [ns].") = JTimeRange( 0.0, 0.3);
74  zap['W'] = make_field(WMin, "minimal histogram contents.") = 100.0;
75  zap['d'] = make_field(debug) = 1;
76 
77  zap(argc, argv);
78  }
79  catch(const exception &error) {
80  FATAL(error.what() << endl);
81  }
82 
83 
84  cout.tie(&cerr);
85 
86  if (!T_ns.is_valid()) {
87  FATAL("Invalid fit range [ns] " << T_ns << endl);
88  }
89 
91 
92  try {
93  load(detectorFile, detector);
94  }
95  catch(const JException& error) {
96  FATAL(error);
97  }
98 
99  if (option.find('S') == string::npos) { option += 'S'; }
100  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
101 
102 
103  TH2D* h2 = NULL;
104 
105  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
106 
107  NOTICE("Processing " << *i << endl);
108 
109  TFile in(i->c_str(), "exist");
110 
111  if (!in.IsOpen()) {
112  FATAL("File " << *i << " not opened." << endl);
113  }
114 
115  TH2D* p = dynamic_cast<TH2D*>(in.Get(h2_t));
116 
117  if (p == NULL) {
118  FATAL("File " << *i << " has no histogram <" << h2_t << ">." << endl);
119  }
120 
121  if (h2 == NULL)
122  h2 = (TH2D*) p->Clone();
123  else
124  h2->Add(p);
125 
126  h2->SetDirectory(0);
127 
128  in.Close();
129  }
130 
131  if (h2 == NULL) {
132  FATAL("No histogram <" << h2_t << ">." << endl);
133  }
134 
135  // auxiliary data structure for fit result
136 
137  struct result_type {
138 
139  result_type() :
140  value(0.0),
141  error(0.0)
142  {}
143 
144  result_type(double value,
145  double error) :
146  value(value),
147  error(error)
148  {}
149 
150  double value;
151  double error;
152  };
153 
154  typedef map<int, result_type> map_type; // module index -> fit result
155 
156  map_type zmap;
157 
158 
159  // histogram fit results
160 
161  const TAxis* x_axis = h2->GetXaxis();
162  const TAxis* y_axis = h2->GetYaxis();
163 
164  TH1D h0("h0", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
165  TH1D hc("hc", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
166  TH1D hq("hq", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
167 
168 
169  TFile out(outputFile.c_str(), "recreate");
170 
171  for (int ix = 1; ix <= x_axis->GetNbins(); ++ix) {
172 
173  if (detector[ix - 1].getFloor() == 0) {
174  continue;
175  }
176 
177  TH1D h1(MAKE_CSTRING(detector[ix - 1].getID() << ".1D"), NULL, y_axis->GetNbins(), y_axis->GetXmin(), y_axis->GetXmax());
178 
179  // start values
180 
181  Double_t ymax = 0.0;
182  Double_t x0 = 0.0; // [ns]
183  Double_t sigma = 4.5; // [ns]
184  Double_t W = 0.0;
185 
186  for (int iy = 1; iy <= y_axis->GetNbins(); ++iy) {
187 
188  h1.SetBinContent(iy, h2->GetBinContent(ix,iy));
189  h1.SetBinError (iy, sqrt(h2->GetBinContent(ix,iy)));
190 
191  const Double_t x = h1.GetBinCenter (iy);
192  const Double_t y = h1.GetBinContent(iy);
193 
194  W += y;
195 
196  if (y > ymax) {
197  ymax = y;
198  x0 = x;
199  }
200  }
201 
202  if (W <= WMin) {
203 
204  WARNING("Not enough entries for slice " << ix << ' ' << W << "; skip fit." << endl);
205 
206  continue;
207  }
208 
209 
210  const Double_t xmin = x0 + T_ns.getLowerLimit();
211  const Double_t xmax = x0 + T_ns.getUpperLimit();
212 
213 
214  TF1* f1 = NULL;
215 
216  if (function == gauss_t) {
217 
218  f1 = new TF1(function.c_str(), "[0]*TMath::Gaus(x,[1],[2])");
219 
220  f1->SetParameter(0, ymax);
221  f1->SetParameter(1, x0);
222  f1->SetParameter(2, sigma);
223 
224  } else if (function == landau_t) {
225 
226  f1 = new TF1(function.c_str(), "[0]*TMath::Landau(x,[1],[2])");
227 
228  f1->SetParameter(0, ymax);
229  f1->SetParameter(1, x0);
230  f1->SetParameter(2, sigma);
231 
232  } else if (function == emg_t) {
233 
234  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]))");
235 
236  f1->SetParameter(0, ymax);
237  f1->SetParameter(1, x0 - sigma);
238  f1->SetParameter(2, sigma);
239  f1->SetParameter(3, 0.04);
240 
241  } else if (function == breitwigner_t) {
242 
243  f1 = new TF1(function.c_str(), "(x <= [1])*[0]*[2]*TMath::BreitWigner(x,[1],[2])+(x > [1])*[0]*[3]*TMath::BreitWigner(x,[1],[3])");
244 
245  f1->SetParameter(0, ymax);
246  f1->SetParameter(1, x0);
247  f1->SetParameter(2, 15.0);
248  f1->SetParameter(3, 25.0);
249 
250  } else {
251 
252  FATAL("Unknown fit function " << function << endl);
253  }
254 
255 
256  DEBUG("Start values:" << endl);
257 
258  for (int i = 0; i != f1->GetNpar(); ++i) {
259  DEBUG(left << setw(12) << f1->GetParName (i) << ' ' <<
260  SCIENTIFIC(12,5) << f1->GetParameter(i) << endl);
261  }
262 
263  // fit
264 
265  TFitResultPtr result = h1.Fit(f1, option.c_str(), "same", xmin, xmax);
266 
267  if (debug >= notice_t || !result->IsValid()) {
268  cout << "Slice: "
269  << setw(4) << ix << ' '
270  << setw(16) << h1.GetName() << ' '
271  << FIXED(7,3) << f1->GetParameter(1) << " +/- "
272  << FIXED(7,3) << f1->GetParError(1) << ' '
273  << FIXED(7,3) << result->Chi2() << '/'
274  << result->Ndf() << ' '
275  << (result->IsValid() ? "" : "failed") << endl;
276  }
277 
278  if (result->IsValid()) {
279  zmap[ix] = result_type(f1->GetParameter(1), f1->GetParError (1));
280  }
281 
282  if (result->Ndf() > 0) {
283  hc.SetBinContent(ix, result->Chi2() / result->Ndf());
284  }
285  hq.SetBinContent(ix, result->IsValid() ? 1.0 : 0.0);
286 
287  h1.Write();
288 
289  delete f1;
290  }
291 
292 
293  if (!zmap.empty()) {
294 
295  // average time offset
296 
297  double t0 = 0.0;
298 
299  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
300  t0 += i->second.value;
301  }
302 
303  t0 /= zmap.size();
304 
305  NOTICE("Average time offset [ns] " << FIXED(7,2) << t0 << endl);
306 
307  for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {
308  i->second.value -= t0;
309  }
310 
311  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
312  h0.SetBinContent(i->first, i->second.value);
313  h0.SetBinError (i->first, i->second.error);
314  }
315 
316  if (!detector.empty()) {
317 
318  const JRange<int> string(make_array(detector.begin(), detector.end(), &JModule::getString));
319  const JRange<int> floor (make_array(detector.begin(), detector.end(), &JModule::getFloor));
320 
321  TH2D hi("hi", NULL,
322  string.getLength() + 1,
323  string.getLowerLimit() - 0.5,
324  string.getUpperLimit() + 0.5,
325  floor.getLength() + 1,
326  floor.getLowerLimit() - 0.5,
327  floor.getUpperLimit() + 0.5);
328 
329  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
330  hi.SetBinContent(detector[i->first - 1].getString(),
331  detector[i->first - 1].getFloor(),
332  i->second.value);
333  }
334 
335  hi.Write();
336  }
337 
338 
339  if (overwriteDetector) {
340 
341  NOTICE("Store calibration data on file " << detectorFile << endl);
342 
343  detector.comment.add(JMeta(argc, argv));
344 
345  for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
346 
347  if (E_ns(i->second.error))
348  detector[i->first - 1].sub(i->second.value);
349  else
350  ERROR("Slice " << setw(4) << i->first << " fit uncertainty " << FIXED(5,2) << i->second.error << " outside specified range (option -E <E_ns>)" << endl);
351  }
352 
353  store(detectorFile, detector);
354  }
355 
356  } else {
357 
358  NOTICE("No calibration results." << endl);
359  }
360 
361  h0 .Write();
362  hc .Write();
363  hq .Write();
364  h2->Write();
365 
366  out.Close();
367 }
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
#define WARNING(A)
Definition: JMessage.hh:65
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
Detector data structure.
Definition: JDetector.hh:89
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [ns]).
notice
Definition: JMessage.hh:32
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:151
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
string outputFile
Data structure for detector geometry and calibration.
T getLength() const
Get length (difference between upper and lower limit).
Definition: JRange.hh:289
I/O formatting auxiliaries.
Detector file.
Definition: JHead.hh:196
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
Definition: JVectorize.hh:54
return result
Definition: JPolint.hh:727
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
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
#define ERROR(A)
Definition: JMessage.hh:66
int debug
debug level
Definition: JSirene.cc:63
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
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 usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
do set_variable DETECTOR_TXT $WORKDIR detector
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 source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:41
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:484