Jpp
Functions
JFit.cc File Reference
#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH1.h"
#include "TGraph.h"
#include "TProfile.h"
#include "TF1.h"
#include "TString.h"
#include "TRegexp.h"
#include "JLang/JType.hh"
#include "JLang/JTypeList.hh"
#include "JLang/JToken.hh"
#include "JTools/JRange.hh"
#include "JTools/JTuple.hh"
#include "JGizmo/JRootObjectID.hh"
#include "JGizmo/JGizmoToolkit.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

General purpose fit program using ROOT. The option -f corresponds to <file name>:<object name>.

The expressions for the fit function, start and fixed values should comply with ROOT class TFormula.

In the expressions of the start and fixed values, names of member methods of corresponding class of the fit object may appear, such as TH1::GetMaximum, TH1::GetRMS, etc., e.g:

     -F "[0]*exp(x/[1])+[2]"
     -@ "p0 = GetMaximum; p1 = 2*GetRMS"
     -= "p2 = 0"

The result of the formulas for the start and fixed values will be evaluated for each histogram separately.

Author
mdejong

Definition in file JFit.cc.

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 92 of file JFit.cc.

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.Contains(regexp) << endl);
178 
179  // option match
180 
181  if (tag.Contains(regexp)) {
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 }
TObject
Definition: JRoot.hh:19
JGIZMO::getParameter
int getParameter(const std::string &text)
Get parameter number from text string.
Definition: JGizmoToolkit.hh:288
JLANG::JType
Auxiliary class for a type holder.
Definition: JType.hh:19
g1
Double_t g1(const Double_t x)
Function.
Definition: JQuantiles.cc:25
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e.
Definition: JParser.hh:63
std::vector
Definition: JSTDTypes.hh:12
JTOOLS::JRange< double >
JTOOLS::j
int j
Definition: JPolint.hh:634
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
ERROR
#define ERROR(A)
Definition: JMessage.hh:66
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
debug
int debug
debug level
Definition: JSirene.cc:59
JGIZMO::getDirectory
TDirectory * getDirectory(const JRootObjectID &id)
Get TDirectory pointer.
Definition: JGizmoToolkit.hh:121
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:518
JLANG::for_each
JObject_t & for_each(JObject_t &object, JType< JTypeList< JHead_t, JTail_t > > typelist)
For each data type method.
Definition: JTypeList.hh:572
JLANG::JTypeList
Type list.
Definition: JTypeList.hh:22
JLANG::JParseError
Exception for parsing value.
Definition: JException.hh:180
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
JLANG::JToken
Wrapper class around string.
Definition: JToken.hh:23
std
Definition: jaanetDictionary.h:36
JEEP::getValue
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JTOOLS::make_tuple
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
JFIT::JFit
Data structure for track fit results.
Definition: JEvt.hh:31