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

General purpose fit program for 3D 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 "TH3.h"
#include "TProfile.h"
#include "TF3.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 3D 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])*exp(-0.5*z*z/([3]*[3])"
     -@ "p0 = GetMaximum; p1 = 2*GetRMS(1); .."
     -= "p3 = 1"

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

Author
mdejong

Definition in file JFit3D.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 101 of file JFit3D.cc.

102{
103 using namespace std;
104 using namespace JPP;
105
106 typedef JToken<';'> JToken_t;
107 typedef JRange<double> JRange_t;
108
109 vector<JRootObjectID> inputFile;
110 string outputFile;
111 string formula;
112 vector<JToken_t> startValues;
113 vector<JToken_t> startErrors;
114 vector<JToken_t> fixedValues;
115 vector<JToken_t> limitValues;
116 JRange_t X;
117 JRange_t Y;
118 JRange_t Z;
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['z'] = make_field(Z, "abscissa range") = JRange_t();
137 zap['O'] = make_field(option, "Fit option") = "";
138 zap['w'] = make_field(writeFits, "write fit function(s) to ROOT file " << "(\"<name>" << _F3 << "\")");
139 zap['M'] = make_field(minimizer, "ROOT minimizer type and algorithm [debug]") = JMinimizer();
140 zap['d'] = make_field(debug) = 1;
141
142 zap(argc, argv);
143 }
144 catch(const exception &error) {
145 FATAL(error.what() << endl);
146 }
147
148
149 if (option.find('O') == string::npos) { option += "O"; }
150 if (option.find("R") == string::npos) { option += "R"; }
151 if (option.find("S") == string::npos) { option += "S"; }
152 //if (option.find('N') == string::npos) { option += "N"; }
153 if (debug == 0 && option.find('Q') == string::npos) { option += "Q"; }
154
155
156 TFile out(outputFile.c_str(), "recreate");
157
158
159 TF3* fcn = new TF3("user", formula.c_str());
160
161 fcn->SetNpx(300);
162 fcn->SetNpy(300);
163 fcn->SetNpz(300);
164
165 if (fcn->IsZombie()) {
166 FATAL("Function: " << formula << " is zombie." << endl);
167 }
168
169 for (vector<JRootObjectID>::const_iterator input = inputFile.begin(); input != inputFile.end(); ++input) {
170
171 DEBUG("Input: " << *input << endl);
172
173 TDirectory* dir = getDirectory(*input);
174
175 if (dir == NULL) {
176 ERROR("File: " << input->getFullFilename() << " not opened." << endl);
177 continue;
178 }
179
180 const TRegexp regexp(input->getObjectName());
181
182 TIter iter(dir->GetListOfKeys());
183
184 for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
185
186 const TString tag(key->GetName());
187
188 DEBUG("Key: " << tag << " match = " << tag.Contains(regexp) << endl);
189
190 // option match
191
192 if (tag.Contains(regexp) && isTObject(key)) {
193
194 TObject* object = key->ReadObj();
195
196
197 // set fit parameters
198
199 try {
200
201 for (vector<JToken_t>::const_iterator i = startValues.begin(); i != startValues.end(); ++i) {
202 fcn->SetParameter(getParameter(*i), getValue(*i,object));
203 }
204
205 for (Int_t i = 0; i != fcn->GetNpar(); ++i) {
206 fcn->SetParError (i, 0.0);
207 }
208
209 for (vector<JToken_t>::const_iterator i = startErrors.begin(); i != startErrors.end(); ++i) {
210 fcn->SetParError (getParameter(*i), getValue(*i,object));
211 }
212
213 for (vector<JToken_t>::const_iterator i = fixedValues.begin(); i != fixedValues.end(); ++i) {
214 fcn->FixParameter(getParameter(*i), getValue(*i,object));
215 }
216
217 for (vector<JToken_t>::const_iterator i = limitValues.begin(); i != limitValues.end(); ++i) {
218 fcn->SetParLimits(getParameter(*i), getValue(*i,0), getValue(*i,1));
219 }
220
221 }
222 catch(JLANG::JParseError& error) {
223 FATAL(error << endl);
224 }
225
226 DEBUG("Start values " << object->GetName() << endl);
227
228 for (int i = 0; i != fcn->GetNpar(); ++i) {
229 DEBUG(left << setw(12) << fcn->GetParName (i) << ' ' <<
230 SCIENTIFIC(12,5) << fcn->GetParameter(i) << endl);
231 }
232
233 Double_t xmin = numeric_limits<Double_t>::max();
234 Double_t xmax = numeric_limits<Double_t>::lowest();
235 Double_t ymin = numeric_limits<Double_t>::max();
236 Double_t ymax = numeric_limits<Double_t>::lowest();
237 Double_t zmin = numeric_limits<Double_t>::max();
238 Double_t zmax = numeric_limits<Double_t>::lowest();
239
240 {
241 TH3* h3 = dynamic_cast<TH3*>(object);
242
243 if (h3 != NULL) {
244 xmin = min(xmin, h3->GetXaxis()->GetXmin());
245 xmax = max(xmax, h3->GetXaxis()->GetXmax());
246 ymin = min(ymin, h3->GetYaxis()->GetXmin());
247 ymax = max(ymax, h3->GetYaxis()->GetXmax());
248 zmin = min(zmin, h3->GetZaxis()->GetXmin());
249 zmax = max(zmax, h3->GetZaxis()->GetXmax());
250 }
251 }
252
253 if (X != JRange_t()) {
254 xmin = X.getLowerLimit();
255 xmax = X.getUpperLimit();
256 }
257
258 if (Y != JRange_t()) {
259 ymin = Y.getLowerLimit();
260 ymax = Y.getUpperLimit();
261 }
262
263 if (Z != JRange_t()) {
264 zmin = Z.getLowerLimit();
265 zmax = Z.getUpperLimit();
266 }
267
268 fcn->SetRange(xmin, ymin, zmin, xmax, ymax, zmax);
269
270
271 // execute fit
272
273 const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
274
275 JFit fit(*object, fcn, option);
276
277 for_each(fit, JType<TH3>());
278
279 const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
280
281 if (fit.result != -1) {
282
283 // output fit results
284
285 NOTICE("Fit values " << object->GetName() << endl);
286 NOTICE("Fit formula " << formula << endl);
287 NOTICE("chi2/NDF " << FIXED(7,3) << fit.result->Chi2() << '/' << fit.result->Ndf() << ' ' << (fit.result->IsValid() ? "" : "failed") << endl);
288 NOTICE("Number of calls " << fit.result->NCalls() << endl);
289 NOTICE("Elapsed time [us] " << setw(8) << chrono::duration_cast<chrono::microseconds>(t1 - t0).count() << endl);
290
291 for (int j = 0; j != fcn->GetNpar(); ++j) {
292 NOTICE(left << setw(12) << fcn->GetParName (j) << ' ' <<
293 SCIENTIFIC(12,5) << fcn->GetParameter(j) << " +/- " <<
294 SCIENTIFIC(12,5) << fcn->GetParError (j) << endl);
295 }
296
297 } else {
298
299 WARNING("Object: not compatible with ROOT Fit." << endl);
300 }
301
302 out.cd();
303
304 object->Write();
305 fcn ->Write();
306
307 if (writeFits) {
308
309 TObject* p = object->Clone(MAKE_CSTRING(object->GetName() << _F3));
310
311 {
312 TH2* h2 = dynamic_cast<TH2*>(p);
313
314 if (h2 != NULL) {
315 for (Int_t ix = 1; ix <= h2->GetXaxis()->GetNbins(); ++ix) {
316 for (Int_t iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy) {
317 for (Int_t iz = 1; iz <= h2->GetZaxis()->GetNbins(); ++iz) {
318
319 const Double_t x = h2->GetXaxis()->GetBinCenter(ix);
320 const Double_t y = h2->GetYaxis()->GetBinCenter(iy);
321 const Double_t z = h2->GetZaxis()->GetBinCenter(iz);
322
323 h2->SetBinContent(ix, iy, fcn->Eval(x,y,z));
324 h2->SetBinError (ix, iy, 0.0);
325 }
326 }
327 }
328 }
329 }
330 }
331 }
332 }
333
334 dir->Close();
335 }
336
337 out.Write();
338 out.Close();
339}
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::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