Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JFit.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <algorithm>
6 
7 #include "TROOT.h"
8 #include "TFile.h"
9 #include "TKey.h"
10 #include "TH1.h"
11 #include "TGraph.h"
12 #include "TProfile.h"
13 #include "TF1.h"
14 #include "TString.h"
15 #include "TRegexp.h"
16 
17 #include "JLang/JType.hh"
18 #include "JLang/JTypeList.hh"
19 #include "JLang/JToken.hh"
20 #include "JTools/JRange.hh"
21 #include "JTools/JTuple.hh"
22 
23 #include "JGizmo/JRootObjectID.hh"
24 #include "JGizmo/JGizmoToolkit.hh"
25 
26 #include "Jeep/JPrint.hh"
27 #include "Jeep/JParser.hh"
28 #include "Jeep/JMessage.hh"
29 
30 
31 namespace {
32 
33  using namespace std;
34  using namespace JPP;
35 
36  /**
37  * Auxiliary class for ROOT fit.
38  */
39  struct JFit {
40  /**
41  * Constructor.
42  *
43  * \param object object
44  */
45  JFit(TObject& object) :
46  object(object)
47  {}
48 
49 
50  /**
51  * Attempt fit for given data type.
52  *
53  * \param type data type
54  * \param tuple argument data
55  */
56  template<class T>
57  void operator()(const JType<T>& type, const JTuple<JTYPELIST<TF1*, string, double, double>::typelist>& tuple)
58  {
59  try {
60  result = dynamic_cast<T&>(object).Fit(tuple.first, tuple.second.first.c_str(), "same", tuple.second.second.first, tuple.second.second.second);
61  }
62  catch(const std::exception&) {}
63  }
64 
65 
66  TObject& object;
67  TFitResultPtr result;
68  };
69 }
70 
71 
72 /**
73  * \file
74  * General purpose fit program using ROOT.
75  * The option <tt>-f</tt> corresponds to <tt><file name>:<object name></tt>.
76  *
77  * The expressions for the fit function, start and fixed values should comply
78  * with ROOT class TFormula.
79  *
80  * In the expressions of the start and fixed values, names of member methods
81  * of corresponding class of the fit object may appear,
82  * such as TH1::GetMaximum, TH1::GetRMS, etc., e.g:
83  * <pre>
84  * -F "[0]*exp(x/[1])+[2]"
85  * -@ "p0 = GetMaximum; p1 = 2*GetRMS"
86  * -= "p2 = 0"
87  * </pre>
88  * The result of the formulas for the start and fixed values will be evaluated
89  * for each histogram separately.
90  * \author mdejong
91  */
92 int main(int argc, char **argv)
93 {
94  using namespace std;
95  using namespace JPP;
96 
97  typedef JLANG::JToken<';'> JToken_t;
98  typedef JRange<double> JRange_t;
99 
100  vector<JRootObjectID> inputFile;
101  string outputFile;
102  string formula;
103  vector<JToken_t> startValues;
104  vector<JToken_t> fixedValues;
105  JRange_t X;
106  JRange_t Y;
107  char project;
108  string option;
109  int debug;
110 
111  try {
112 
113  JParser<> zap("General purpose fit program using ROOT.");
114 
115  zap['f'] = make_field(inputFile, "<input file>:<object name>");
116  zap['o'] = make_field(outputFile, "ROOT file with fit results") = "fit.root";
117  zap['F'] = make_field(formula, "fit formula, e.g: \"[0]+[1]*x\"");
118  zap['@'] = make_field(startValues, "start values, e.g: \"p0 = GetMaximum;\"");
119  zap['='] = make_field(fixedValues, "fixed values, e.g: \"p0 = GetMaximum;\"");
120  zap['x'] = make_field(X, "abscissa range") = JRange_t();
121  zap['y'] = make_field(Y, "ordinate range") = JRange_t();
122  zap['P'] = make_field(project, "projection") = '\0', 'x', 'X', 'y', 'Y';
123  zap['O'] = make_field(option, "Fit option") = "";
124  zap['d'] = make_field(debug) = 1;
125 
126  zap['='] = JPARSER::initialised();
127 
128  zap(argc, argv);
129  }
130  catch(const exception &error) {
131  FATAL(error.what() << endl);
132  }
133 
134 
135  if (option.find('O') == string::npos) { option += "O"; }
136  //if (option.find('N') == string::npos) { option += "N"; }
137  if (debug == 0 && option.find('Q') == string::npos) { option += "Q"; }
138 
139  bool px = (project == 'x' || project == 'X'); // projection on x-axis
140  bool py = (project == 'y' || project == 'Y'); // projection on y-axis
141 
142  if (py) {
143  swap(Y, X); // Y becomes x-axis range and X becomes y-axis range
144  }
145 
146 
147  TFile out(outputFile.c_str(), "recreate");
148 
149 
150  TF1* fcn = new TF1("user", formula.c_str());
151 
152  fcn->SetNpx(1000);
153 
154  if (fcn->IsZombie()) {
155  FATAL("Function: " << formula << " is zombie." << endl);
156  }
157 
158  for (vector<JRootObjectID>::const_iterator input = inputFile.begin(); input != inputFile.end(); ++input) {
159 
160  DEBUG("Input: " << *input << endl);
161 
162  TDirectory* dir = getDirectory(*input);
163 
164  if (dir == NULL) {
165  ERROR("File: " << input->getFullFilename() << " not opened." << endl);
166  continue;
167  }
168 
169  const TRegexp regexp(input->getObjectName());
170 
171  TIter iter(dir->GetListOfKeys());
172 
173  for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
174 
175  const TString tag(key->GetName());
176 
177  DEBUG("Key: " << tag << " match = " << tag.Index(regexp) << " (-1 /= OK)" << endl);
178 
179  // option match
180 
181  if (tag.Index(regexp) != -1) {
182 
183  TObject* object = key->ReadObj();
184 
185  try {
186 
187  TProfile& h1 = dynamic_cast<TProfile&>(*object);
188 
189  object = h1.ProjectionX();
190  }
191  catch(exception&) {}
192 
193  try {
194 
195  TH2& h2 = dynamic_cast<TH2&>(*object);
196 
197  if (px) {
198 
199  object = h2.ProjectionX("_px",
200  h2.GetYaxis()->FindBin(Y.getLowerLimit()),
201  h2.GetYaxis()->FindBin(Y.getUpperLimit()));
202 
203  } else if (py) {
204 
205  object = h2.ProjectionY("_py",
206  h2.GetXaxis()->FindBin(Y.getLowerLimit()),
207  h2.GetXaxis()->FindBin(Y.getUpperLimit()));
208 
209  } else {
210 
211  ERROR("For 2D histograms, use option option -P for projections." << endl);
212 
213  continue;
214  }
215  }
216  catch(exception&) {}
217 
218 
219  // set fit parameters
220 
221  try {
222 
223  for (vector<JToken_t>::const_iterator j = startValues.begin(); j != startValues.end(); ++j) {
224  fcn->SetParameter(getParameter(*j), getValue(*j,object));
225  }
226 
227  for (vector<JToken_t>::const_iterator j = fixedValues.begin(); j != fixedValues.end(); ++j) {
228  fcn->FixParameter(getParameter(*j), getValue(*j,object));
229  }
230  }
231  catch(JLANG::JParseError& error) {
232  FATAL(error << endl);
233  }
234 
235  DEBUG("Start values " << object->GetName() << endl);
236 
237  for (int j = 0; j != fcn->GetNpar(); ++j) {
238  DEBUG(left << setw(12) << fcn->GetParName (j) << ' ' <<
239  SCIENTIFIC(12,5) << fcn->GetParameter(j) << endl);
240  }
241 
242  Double_t xmin = +numeric_limits<Double_t>::max();
243  Double_t xmax = -numeric_limits<Double_t>::max();
244 
245  {
246  TH1* h1 = dynamic_cast<TH1*>(object);
247 
248  if (h1 != NULL) {
249  xmin = min(xmin, h1->GetXaxis()->GetXmin());
250  xmax = max(xmax, h1->GetXaxis()->GetXmax());
251  }
252  }
253 
254  {
255  TGraph* g1 = dynamic_cast<TGraph*>(object);
256 
257  if (g1 != NULL) {
258  for (Int_t i = 0; i != g1->GetN(); ++i) {
259  xmin = min(xmin, g1->GetX()[i]);
260  xmax = max(xmax, g1->GetX()[i]);
261  }
262  }
263  }
264 
265  if (X != JRange_t()) {
266  xmin = X.getLowerLimit();
267  xmax = X.getUpperLimit();
268  }
269 
270 
271  // execute fit
272 
273  JFit fit(*object);
274 
275  for_each(fit, JType<JTYPELIST<TH1, TGraph>::typelist>(), JPP::make_tuple(fcn, option, xmin, xmax));
276 
277 
278  if (fit.result != -1) {
279 
280  // output fit results
281 
282  NOTICE("Fit values " << object->GetName() << endl);
283  NOTICE("Fit formula " << formula << endl);
284 
285  for (int j = 0; j != fcn->GetNpar(); ++j) {
286  NOTICE(left << setw(12) << fcn->GetParName (j) << ' ' <<
287  SCIENTIFIC(12,5) << fcn->GetParameter(j) << " +/- " <<
288  SCIENTIFIC(12,5) << fcn->GetParError (j) << endl);
289  }
290 
291  } else {
292 
293  WARNING("Object: not compatible with ROOT Fit." << endl);
294  }
295 
296  out.cd();
297 
298  object->Write();
299  }
300  }
301 
302  dir->Close();
303  }
304 
305  out.Write();
306  out.Close();
307 }
Utility class to parse command line options.
Definition: JParser.hh:1410
Template data structure.
Definition: JTuple.hh:45
#define WARNING(A)
Definition: JMessage.hh:63
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
int getParameter(const std::string &text)
Get parameter number from text string.
Definition: JRoot.hh:19
Empty structure for specification of parser element that is initialised (i.e.
Definition: JParser.hh:64
Auxiliary class for a type holder.
Definition: JType.hh:19
string outputFile
Data structure for track fit results.
Definition: JEvt.hh:32
Type list.
Definition: JTypeList.hh:22
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
#define NOTICE(A)
Definition: JMessage.hh:62
#define ERROR(A)
Definition: JMessage.hh:64
int debug
debug level
Definition: JSirene.cc:59
General purpose messaging.
JObject_t & for_each(JObject_t &object, JType< JTypeList< JHead_t, JTail_t > > typelist)
For each data type method.
Definition: JTypeList.hh:572
#define FATAL(A)
Definition: JMessage.hh:65
Auxiliary class to define a range between two values.
Utility class to parse command line options.
Wrapper class around string.
Definition: JToken.hh:23
Data structure based on type list.
Exception for parsing value.
Definition: JException.hh:162
TDirectory * getDirectory(const JRootObjectID &id)
Get TDirectory pointer.
JTuple< JTypeList< JNullType, JNullType > > make_tuple(const JNullType &b, const JNullType &c, const JNullType &d, const JNullType &e, const JNullType &f, const JNullType &g, const JNullType &h, const JNullType &i, const JNullType &j, const JNullType &k, const JNullType &l, const JNullType &m, const JNullType &n, const JNullType &o, const JNullType &p, const JNullType &q, const JNullType &r, const JNullType &s, const JNullType &t, const JNullType &u, const JNullType &v, const JNullType &w, const JNullType &x, const JNullType &y, const JNullType &z)
Helper method for tuple.
Definition: JTuple.hh:724
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:498
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
Double_t g1(const Double_t x)
Function.
Definition: JQuantiles.cc:25
int main(int argc, char *argv[])
Definition: Main.cpp:15