Jpp 19.3.0-rc.3
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JFit2D.cc File Reference

General purpose fit program for 2D ROOT objects. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include <chrono>
#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH2.h"
#include "TGraph2D.h"
#include "TGraph2DErrors.h"
#include "TF2.h"
#include "TFitResult.h"
#include "TString.h"
#include "TRegexp.h"
#include "JROOT/JMinimizer.hh"
#include "JLang/JType.hh"
#include "JLang/JTypeList.hh"
#include "JLang/JToken.hh"
#include "JTools/JRange.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 for 2D ROOT objects.


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(-0.5*x*x/([1]*[1])*exp(-0.5*y*y/([2]*[2])"
     -@ "p0 = GetMaximum; p1 = 2*GetRMS(1); .."
     -= "p2 = 1"

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

Author
mdejong

Definition in file JFit2D.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 102 of file JFit2D.cc.

103{
104 using namespace std;
105 using namespace JPP;
106
107 typedef JToken<';'> JToken_t;
108 typedef JRange<double> JRange_t;
109
110 vector<JRootObjectID> inputFile;
111 string outputFile;
112 string formula;
113 vector<JToken_t> startValues;
114 vector<JToken_t> startErrors;
115 vector<JToken_t> fixedValues;
116 vector<JToken_t> limitValues;
117 JRange_t X;
118 JRange_t Y;
119 string option;
120 bool writeFits;
121 int debug;
122
123 try {
124
125 JParser<> zap("General purpose fit program for 2D ROOT objects.");
126
127 zap['f'] = make_field(inputFile, "<input file>:<object name>");
128 zap['o'] = make_field(outputFile, "ROOT file with fit results") = "fit.root";
129 zap['F'] = make_field(formula, "fit formula, e.g: \"[0]+[1]*x\"");
130 zap['@'] = make_field(startValues, "start values, e.g: \"p0 = GetMaximum;\"");
131 zap['E'] = make_field(startErrors, "start errors, e.g: \"p0 = 0.01 * GetMaximum;\"") = JPARSER::initialised();
132 zap['='] = make_field(fixedValues, "fixed values, e.g: \"p0 = GetMaximum;\"") = JPARSER::initialised();
133 zap['R'] = make_field(limitValues, "limit values, e.g: \"p0 = <lower limit> <upper limit>;\"") = JPARSER::initialised();
134 zap['x'] = make_field(X, "abscissa range") = JRange_t();
135 zap['y'] = make_field(Y, "abscissa range") = JRange_t();
136 zap['O'] = make_field(option, "Fit option") = "";
137 zap['w'] = make_field(writeFits, "write fit function(s) to ROOT file " << "(\"<name>" << _F2 << "\")");
138 zap['M'] = make_field(minimizer, "ROOT minimizer type and algorithm [debug]") = JMinimizer();
139 zap['d'] = make_field(debug) = 1;
140
141 zap(argc, argv);
142 }
143 catch(const exception &error) {
144 FATAL(error.what() << endl);
145 }
146
147
148 if (option.find('O') == string::npos) { option += "O"; }
149 if (option.find("R") == string::npos) { option += "R"; }
150 if (option.find("S") == string::npos) { option += "S"; }
151 //if (option.find('N') == string::npos) { option += "N"; }
152 if (debug == 0 && option.find('Q') == string::npos) { option += "Q"; }
153
154
155 TFile out(outputFile.c_str(), "recreate");
156
157
158 TF2* fcn = new TF2("user", formula.c_str());
159
160 fcn->SetNpx(1000);
161 fcn->SetNpy(1000);
162
163 if (fcn->IsZombie()) {
164 FATAL("Function: " << formula << " is zombie." << endl);
165 }
166
167 for (vector<JRootObjectID>::const_iterator input = inputFile.begin(); input != inputFile.end(); ++input) {
168
169 DEBUG("Input: " << *input << endl);
170
171 TDirectory* dir = getDirectory(*input);
172
173 if (dir == NULL) {
174 ERROR("File: " << input->getFullFilename() << " not opened." << endl);
175 continue;
176 }
177
178 const TRegexp regexp(input->getObjectName());
179
180 TIter iter(dir->GetListOfKeys());
181
182 for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
183
184 const TString tag(key->GetName());
185
186 DEBUG("Key: " << tag << " match = " << tag.Contains(regexp) << endl);
187
188 // option match
189
190 if (tag.Contains(regexp) && isTObject(key)) {
191
192 TObject* object = key->ReadObj();
193
194
195 // set fit parameters
196
197 try {
198
199 for (vector<JToken_t>::const_iterator i = startValues.begin(); i != startValues.end(); ++i) {
200 fcn->SetParameter(getParameter(*i), getValue(*i,object));
201 }
202
203 for (Int_t i = 0; i != fcn->GetNpar(); ++i) {
204 fcn->SetParError(i, 0.0);
205 }
206
207 for (vector<JToken_t>::const_iterator i = startErrors.begin(); i != startErrors.end(); ++i) {
208 fcn->SetParError (getParameter(*i), getValue(*i,object));
209 }
210
211 for (vector<JToken_t>::const_iterator i = fixedValues.begin(); i != fixedValues.end(); ++i) {
212 fcn->FixParameter(getParameter(*i), getValue(*i,object));
213 }
214
215 for (vector<JToken_t>::const_iterator i = limitValues.begin(); i != limitValues.end(); ++i) {
216 fcn->SetParLimits(getParameter(*i), getValue(*i,0), getValue(*i,1));
217 }
218 }
219 catch(JLANG::JParseError& error) {
220 FATAL(error << endl);
221 }
222
223 DEBUG("Start values " << object->GetName() << endl);
224
225 for (int i = 0; i != fcn->GetNpar(); ++i) {
226 DEBUG(left << setw(12) << fcn->GetParName (i) << ' ' <<
227 SCIENTIFIC(12,5) << fcn->GetParameter(i) << endl);
228 }
229
230 Double_t xmin = numeric_limits<Double_t>::max();
231 Double_t xmax = numeric_limits<Double_t>::lowest();
232 Double_t ymin = numeric_limits<Double_t>::max();
233 Double_t ymax = numeric_limits<Double_t>::lowest();
234
235 {
236 TH2* h2 = dynamic_cast<TH2*>(object);
237
238 if (h2 != NULL) {
239
240 xmin = min(xmin, h2->GetXaxis()->GetXmin());
241 xmax = max(xmax, h2->GetXaxis()->GetXmax());
242 ymin = min(ymin, h2->GetYaxis()->GetXmin());
243 ymax = max(ymax, h2->GetYaxis()->GetXmax());
244
245 if (X != JRange_t()) { h2->GetXaxis()->SetRangeUser(X.getLowerLimit(), X.getUpperLimit()); }
246 if (Y != JRange_t()) { h2->GetYaxis()->SetRangeUser(Y.getLowerLimit(), Y.getUpperLimit()); }
247 }
248 }
249
250 {
251 TGraph2D* g2 = dynamic_cast<TGraph2D*>(object);
252
253 if (g2 != NULL) {
254 for (Int_t i = 0; i != g2->GetN(); ++i) {
255 xmin = min(xmin, g2->GetX()[i]);
256 xmax = max(xmax, g2->GetX()[i]);
257 ymin = min(ymin, g2->GetY()[i]);
258 ymax = max(ymax, g2->GetY()[i]);
259 }
260 }
261 }
262
263 if (X != JRange_t()) {
264 xmin = X.getLowerLimit();
265 xmax = X.getUpperLimit();
266 }
267
268 if (Y != JRange_t()) {
269 ymin = Y.getLowerLimit();
270 ymax = Y.getUpperLimit();
271 }
272
273 fcn->SetRange(xmin, ymin, xmax, ymax);
274
275
276 // execute fit
277
278 const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
279
280 JFit fit(*object, fcn, option);
281
282 for_each(fit, JType<JTYPELIST<TH2, TGraph2D>::typelist>());
283
284 const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
285
286 if (fit.result != -1) {
287
288 // output fit results
289
290 NOTICE("Fit values " << object->GetName() << endl);
291 NOTICE("Fit formula " << formula << endl);
292 NOTICE("chi2/NDF " << FIXED(7,3) << fit.result->Chi2() << '/' << fit.result->Ndf() << ' ' << (fit.result->IsValid() ? "" : "failed") << endl);
293 NOTICE("Number of calls " << fit.result->NCalls() << endl);
294 NOTICE("Elapsed time [us] " << setw(8) << chrono::duration_cast<chrono::microseconds>(t1 - t0).count() << endl);
295
296 for (int j = 0; j != fcn->GetNpar(); ++j) {
297 NOTICE(left << setw(12) << fcn->GetParName (j) << ' ' <<
298 SCIENTIFIC(12,5) << fcn->GetParameter(j) << " +/- " <<
299 SCIENTIFIC(12,5) << fcn->GetParError (j) << endl);
300 }
301
302 } else {
303
304 WARNING("Object: not compatible with ROOT Fit." << endl);
305 }
306
307 out.cd();
308
309 object->Write();
310 fcn ->Write();
311
312 if (writeFits) {
313
314 TObject* p = object->Clone(MAKE_CSTRING(object->GetName() << _F2));
315
316 {
317 TH2* h2 = dynamic_cast<TH2*>(p);
318
319 if (h2 != NULL) {
320 for (Int_t ix = 1; ix <= h2->GetXaxis()->GetNbins(); ++ix) {
321 for (Int_t iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy) {
322
323 const Double_t x = h2->GetXaxis()->GetBinCenter(ix);
324 const Double_t y = h2->GetYaxis()->GetBinCenter(iy);
325
326 h2->SetBinContent(ix, iy, fcn->Eval(x,y));
327 h2->SetBinError (ix, iy, 0.0);
328 }
329 }
330 }
331 }
332
333 {
334 TGraph2D* g2 = dynamic_cast<TGraph2D*>(p);
335
336 if (g2 != NULL) {
337 for (Int_t i = 0; i != g2->GetN(); ++i) {
338
339 const Double_t x = g2->GetX()[i];
340 const Double_t y = g2->GetY()[i];
341
342 g2->GetZ()[i] = fcn->Eval(x,y);
343 }
344 }
345 }
346
347 {
348 TGraph2DErrors* g2 = dynamic_cast<TGraph2DErrors*>(p);
349
350 if (g2 != NULL) {
351 for (Int_t i = 0; i != g2->GetN(); ++i) {
352 g2->SetPointError(i, 0.0, 0.0, 0.0);
353 }
354 }
355 }
356 }
357 }
358 }
359
360 dir->Close();
361 }
362
363 out.Write();
364 out.Close();
365}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
NOTICE
#define NOTICE(A)
Definition JMessage.hh:64
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
WARNING
#define WARNING(A)
Definition JMessage.hh:65
minimizer
static JMinimizer minimizer
ROOT minimizer.
Definition JMinimizer.hh:80
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
JLANG::JParseError
Exception for parsing value.
Definition JException.hh:198
JLANG::JToken
Wrapper class around string.
Definition JToken.hh:26
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JTOOLS::JRange
Range of values.
Definition JRange.hh:42
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_LEGENDRE::xmax
const double xmax
Definition JQuadrature.cc:24
GAUSS_LEGENDRE::xmin
const double xmin
Definition JQuadrature.cc:23
JEEP::getValue
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition JScale.hh:47
JGIZMO::getParameter
int getParameter(const std::string &text)
Get parameter number from text string.
Definition JGizmoToolkit.hh:386
JGIZMO::getDirectory
TDirectory * getDirectory(const JRootObjectID &id)
Get TDirectory pointer.
Definition JGizmoToolkit.hh:144
JGIZMO::isTObject
bool isTObject(const TKey *key)
Check if given key corresponds to a TObject.
Definition JGizmoToolkit.hh:1595
JLANG::for_each
void for_each(JObject_t &object, JType< JTypeList< JHead_t, JTail_t > > typelist)
For each data type method.
Definition JTypeList.hh:414
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:43
JTOOLS::j
int j
Definition JPolint.hh:801
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JAANET::JRange_t
Type definition of range.
Definition JHead.hh:43
JACOUSTICS::JFit
Acoustic single fit.
Definition JAcoustics/JEvt.hh:36
JLANG::JTypeList
Type list.
Definition JTypeList.hh:23
JLANG::JType
Auxiliary class for a type holder.
Definition JType.hh:19
JMinimizer
Auxiliary data structure to define ROOT minimizer.
Definition JMinimizer.hh:14
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488
Generated by doxygen 1.12.0