Logo ROOT   6.16/01
Reference Guide
HFitInterface.cxx
Go to the documentation of this file.
1// @(#)root/hist:$Id$
2// Author: L. Moneta Thu Aug 31 10:40:20 2006
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7 * *
8 * *
9 **********************************************************************/
10
11// Implementation file for class TH1Interface
12
13#include "HFitInterface.h"
14
15#include "Fit/BinData.h"
16#include "Fit/SparseData.h"
17#include "Fit/FitResult.h"
18#include "Math/IParamFunction.h"
19
20#include <vector>
21
22#include <cassert>
23#include <cmath>
24
25#include "TH1.h"
26#include "THnBase.h"
27#include "TF1.h"
28#include "TGraph2D.h"
29#include "TGraph.h"
30#include "TGraphErrors.h"
31// #include "TGraphErrors.h"
32// #include "TGraphBentErrors.h"
33// #include "TGraphAsymmErrors.h"
34#include "TMultiGraph.h"
35#include "TList.h"
36#include "TError.h"
37
38
39//#define DEBUG
40#ifdef DEBUG
41#include "TClass.h"
42#include <iostream>
43#endif
44
45
46namespace ROOT {
47
48namespace Fit {
49
50// add a namespace to distinguish from the Graph functions
51namespace HFitInterface {
52
53
54bool IsPointOutOfRange(const TF1 * func, const double * x) {
55 // function to check if a point is outside range
56 if (func ==0) return false;
57 return !func->IsInside(x);
58}
59
60bool AdjustError(const DataOptions & option, double & error, double value = 1) {
61 // adjust the given error according to the option
62 // return false when point must be skipped.
63 // When point error = 0, the point is kept if the option UseEmpty is set or if
64 // fErrors1 is set and the point value is not zero.
65 // The value should be used only for points representing counts (histograms), not for the graph.
66 // In the graph points with zero errors are by default skipped indepentently of the value.
67 // If one wants to keep the points, the option fUseEmpty must be set
68
69 if (error <= 0) {
70 if (option.fUseEmpty || (option.fErrors1 && std::abs(value) > 0 ) )
71 error = 1.; // set error to 1
72 else
73 return false; // skip bins with zero errors or empty
74 } else if (option.fErrors1)
75 error = 1; // set all error to 1 for non-empty bins
76 return true;
77}
78
79void ExamineRange(const TAxis * axis, std::pair<double,double> range,int &hxfirst,int &hxlast) {
80 // examine the range given with the pair on the given histogram axis
81 // correct in case the bin values hxfirst hxlast
82 double xlow = range.first;
83 double xhigh = range.second;
84#ifdef DEBUG
85 std::cout << "xlow " << xlow << " xhigh = " << xhigh << std::endl;
86#endif
87 // ignore ranges specified outside histogram range
88 int ilow = axis->FindFixBin(xlow);
89 int ihigh = axis->FindFixBin(xhigh);
90 if (ilow > hxlast || ihigh < hxfirst) {
91 Warning("ROOT::Fit::FillData","fit range is outside histogram range, no fit data for %s",axis->GetName());
92 }
93 // consider only range defined with-in histogram not oustide. Always exclude underflow/overflow
94 hxfirst = std::min( std::max( ilow, hxfirst), hxlast+1) ;
95 hxlast = std::max( std::min( ihigh, hxlast), hxfirst-1) ;
96 // exclude bins where range coverage is less than half bin width
97 if (hxfirst < hxlast) {
98 if ( axis->GetBinCenter(hxfirst) < xlow) hxfirst++;
99 if ( axis->GetBinCenter(hxlast) > xhigh) hxlast--;
100 }
101}
102
103
104} // end namespace HFitInterface
105
106
107void FillData(BinData & dv, const TH1 * hfit, TF1 * func)
108{
109 // Function to fill the binned Fit data structure from a TH1
110 // The dimension of the data is the same of the histogram dimension
111 // The function pointer is need in case of integral is used and to reject points
112 // rejected in the function
113
114 // the TF1 pointer cannot be constant since EvalPar and InitArgs are not const methods
115
116 // get fit option
117 const DataOptions & fitOpt = dv.Opt();
118
119
120 // store instead of bin center the bin edges
121 bool useBinEdges = fitOpt.fIntegral || fitOpt.fBinVolume;
122
123 assert(hfit != 0);
124
125 //std::cout << "creating Fit Data from histogram " << hfit->GetName() << std::endl;
126
127 int hxfirst = hfit->GetXaxis()->GetFirst();
128 int hxlast = hfit->GetXaxis()->GetLast();
129
130 int hyfirst = hfit->GetYaxis()->GetFirst();
131 int hylast = hfit->GetYaxis()->GetLast();
132
133 int hzfirst = hfit->GetZaxis()->GetFirst();
134 int hzlast = hfit->GetZaxis()->GetLast();
135
136 // function by default has same range (use that one if requested otherwise use data one)
137
138
139 // get the range (add the function range ??)
140 // to check if inclusion/exclusion at end/point
141 const DataRange & range = dv.Range();
142 if (range.Size(0) != 0) {
143 HFitInterface::ExamineRange( hfit->GetXaxis(), range(0), hxfirst, hxlast);
144 if (range.Size(0) > 1 ) {
145 Warning("ROOT::Fit::FillData","support only one range interval for X coordinate");
146 }
147 }
148
149 if (hfit->GetDimension() > 1 && range.Size(1) != 0) {
150 HFitInterface::ExamineRange( hfit->GetYaxis(), range(1), hyfirst, hylast);
151 if (range.Size(1) > 1 )
152 Warning("ROOT::Fit::FillData","support only one range interval for Y coordinate");
153 }
154
155 if (hfit->GetDimension() > 2 && range.Size(2) != 0) {
156 HFitInterface::ExamineRange( hfit->GetZaxis(), range(2), hzfirst, hzlast);
157 if (range.Size(2) > 1 )
158 Warning("ROOT::Fit::FillData","support only one range interval for Z coordinate");
159 }
160
161
162 int n = (hxlast-hxfirst+1)*(hylast-hyfirst+1)*(hzlast-hzfirst+1);
163
164#ifdef DEBUG
165 std::cout << "THFitInterface: ifirst = " << hxfirst << " ilast = " << hxlast
166 << " total bins " << n
167 << std::endl;
168#endif
169
170 // reserve n for more efficient usage
171 //dv.Data().reserve(n);
172
173 int hdim = hfit->GetDimension();
174 int ndim = hdim;
175 // case of function dimension less than histogram
176 if (func !=0 && func->GetNdim() == hdim-1) ndim = hdim-1;
177
178 assert( ndim > 0 );
179 //typedef BinPoint::CoordData CoordData;
180 //CoordData x = CoordData( hfit->GetDimension() );
181 dv.Initialize(n,ndim);
182
183 double x[3];
184 double s[3];
185
186 int binx = 0;
187 int biny = 0;
188 int binz = 0;
189
190 const TAxis *xaxis = hfit->GetXaxis();
191 const TAxis *yaxis = hfit->GetYaxis();
192 const TAxis *zaxis = hfit->GetZaxis();
193
194 for ( binx = hxfirst; binx <= hxlast; ++binx) {
195 if (useBinEdges) {
196 x[0] = xaxis->GetBinLowEdge(binx);
197 s[0] = xaxis->GetBinUpEdge(binx);
198 }
199 else
200 x[0] = xaxis->GetBinCenter(binx);
201
202
203 for ( biny = hyfirst; biny <= hylast; ++biny) {
204 if (useBinEdges) {
205 x[1] = yaxis->GetBinLowEdge(biny);
206 s[1] = yaxis->GetBinUpEdge(biny);
207 }
208 else
209 x[1] = yaxis->GetBinCenter(biny);
210
211 for ( binz = hzfirst; binz <= hzlast; ++binz) {
212 if (useBinEdges) {
213 x[2] = zaxis->GetBinLowEdge(binz);
214 s[2] = zaxis->GetBinUpEdge(binz);
215 }
216 else
217 x[2] = zaxis->GetBinCenter(binz);
218
219 // need to evaluate function to know about rejected points
220 // hugly but no other solutions
221 if (func != 0) {
222 TF1::RejectPoint(false);
223 (*func)( &x[0] ); // evaluate using stored function parameters
224 if (TF1::RejectedPoint() ) continue;
225 }
226
227
228 double value = hfit->GetBinContent(binx, biny, binz);
229 double error = hfit->GetBinError(binx, biny, binz);
230 if (!HFitInterface::AdjustError(fitOpt,error,value) ) continue;
231
232 if (ndim == hdim -1) {
233 // case of fitting a function with dimension -1
234 // point error is bin width y / sqrt(N) where N is nuber of entries in the bin
235 // normalization of error will be wrong - but they will be rescaled in the fit
236 if (hdim == 2) dv.Add( x, x[1], yaxis->GetBinWidth(biny) / error );
237 if (hdim == 3) dv.Add( x, x[2], zaxis->GetBinWidth(binz) / error );
238 } else {
239 dv.Add( x, value, error );
240 if (useBinEdges) {
241 dv.AddBinUpEdge( s );
242 }
243 }
244
245
246#ifdef DEBUG
247 std::cout << "bin " << binx << " add point " << x[0] << " " << hfit->GetBinContent(binx) << std::endl;
248#endif
249
250 } // end loop on z bins
251 } // end loop on y bins
252 } // end loop on x axis
253
254
255#ifdef DEBUG
256 std::cout << "THFitInterface::FillData: Hist FitData size is " << dv.Size() << std::endl;
257#endif
258
259}
260
261////////////////////////////////////////////////////////////////////////////////
262/// Compute rough values of parameters for an exponential
263
265{
266 unsigned int n = data.Size();
267 if (n == 0) return;
268
269 // find xmin and xmax of the data
270 double valxmin;
271 double xmin = *(data.GetPoint(0,valxmin));
272 double xmax = xmin;
273 double valxmax = valxmin;
274
275 for (unsigned int i = 1; i < n; ++ i) {
276 double val;
277 double x = *(data.GetPoint(i,val) );
278 if (x < xmin) {
279 xmin = x;
280 valxmin = val;
281 }
282 else if (x > xmax) {
283 xmax = x;
284 valxmax = val;
285 }
286 }
287
288 // avoid negative values of valxmin/valxmax
289 if (valxmin <= 0 && valxmax > 0 ) valxmin = valxmax;
290 else if (valxmax <=0 && valxmin > 0) valxmax = valxmin;
291 else if (valxmin <=0 && valxmax <= 0) { valxmin = 1; valxmax = 1; }
292
293 double slope = std::log( valxmax/valxmin) / (xmax - xmin);
294 double constant = std::log(valxmin) - slope * xmin;
295 f1->SetParameters(constant, slope);
296}
297
298
299////////////////////////////////////////////////////////////////////////////////
300/// Compute Initial values of parameters for a gaussian
301/// derived from function H1InitGaus defined in TH1.cxx
302
304{
305
306 static const double sqrtpi = 2.506628;
307
308 // - Compute mean value and RMS of the data
309 unsigned int n = data.Size();
310 if (n == 0) return;
311 double sumx = 0;
312 double sumx2 = 0;
313 double allcha = 0;
314 double valmax = 0;
315 double rangex = data.Coords(n-1)[0] - data.Coords(0)[0];
316 // to avoid binwidth = 0 set arbitrarly to 1
317 double binwidth = 1;
318 if ( rangex > 0) binwidth = rangex;
319 double x0 = 0;
320 for (unsigned int i = 0; i < n; ++ i) {
321 double val;
322 double x = *(data.GetPoint(i,val) );
323 sumx += val*x;
324 sumx2 += val*x*x;
325 allcha += val;
326 if (val > valmax) valmax = val;
327 if (i > 0) {
328 double dx = x - x0;
329 if (dx < binwidth) binwidth = dx;
330 }
331 x0 = x;
332 }
333
334 if (allcha <= 0) return;
335 double mean = sumx/allcha;
336 double rms = sumx2/allcha - mean*mean;
337
338
339 if (rms > 0)
340 rms = std::sqrt(rms);
341 else
342 rms = binwidth*n/4;
343
344
345 //if the distribution is really gaussian, the best approximation
346 //is binwidx*allcha/(sqrtpi*rms)
347 //However, in case of non-gaussian tails, this underestimates
348 //the normalisation constant. In this case the maximum value
349 //is a better approximation.
350 //We take the average of both quantities
351
352// printf("valmax %f other %f bw %f allcha %f rms %f \n",valmax, binwidth*allcha/(sqrtpi*rms),
353// binwidth, allcha,rms );
354
355 double constant = 0.5*(valmax+ binwidth*allcha/(sqrtpi*rms));
356
357
358 //In case the mean value is outside the histo limits and
359 //the RMS is bigger than the range, we take
360 // mean = center of bins
361 // rms = half range
362// Double_t xmin = curHist->GetXaxis()->GetXmin();
363// Double_t xmax = curHist->GetXaxis()->GetXmax();
364// if ((mean < xmin || mean > xmax) && rms > (xmax-xmin)) {
365// mean = 0.5*(xmax+xmin);
366// rms = 0.5*(xmax-xmin);
367// }
368
369 f1->SetParameter(0,constant);
370 f1->SetParameter(1,mean);
371 f1->SetParameter(2,rms);
372 f1->SetParLimits(2,0,10*rms);
373
374
375#ifdef DEBUG
376 std::cout << "Gaussian initial par values" << constant << " " << mean << " " << rms << std::endl;
377#endif
378
379}
380
381////////////////////////////////////////////////////////////////////////////////
382/// Compute Initial values of parameters for a gaussian
383/// derived from function H1InitGaus defined in TH1.cxx
384
386{
387
388 static const double sqrtpi = 2.506628;
389
390 // - Compute mean value and RMS of the data
391 unsigned int n = data.Size();
392 if (n == 0) return;
393 double sumx = 0, sumy = 0;
394 double sumx2 = 0, sumy2 = 0;
395 double allcha = 0;
396 double valmax = 0;
397 double rangex = data.Coords(n-1)[0] - data.Coords(0)[0];
398 double rangey = data.Coords(n-1)[1] - data.Coords(0)[1];
399 // to avoid binwidthx = 0 set arbitrarly to 1
400 double binwidthx = 1, binwidthy = 1;
401 if ( rangex > 0) binwidthx = rangex;
402 if ( rangey > 0) binwidthy = rangey;
403 double x0 = 0, y0 = 0;
404 for (unsigned int i = 0; i < n; ++i) {
405 double val;
406 const double *coords = data.GetPoint(i,val);
407 double x = coords[0], y = coords[1];
408 sumx += val*x;
409 sumy += val*y;
410 sumx2 += val*x*x;
411 sumy2 += val*y*y;
412 allcha += val;
413 if (val > valmax) valmax = val;
414 if (i > 0) {
415 double dx = x - x0;
416 if (dx < binwidthx) binwidthx = dx;
417 double dy = y - y0;
418 if (dy < binwidthy) binwidthy = dy;
419 }
420 x0 = x;
421 y0 = y;
422 }
423
424 if (allcha <= 0) return;
425 double meanx = sumx/allcha, meany = sumy/allcha;
426 double rmsx = sumx2/allcha - meanx*meanx;
427 double rmsy = sumy2/allcha - meany*meany;
428
429
430 if (rmsx > 0)
431 rmsx = std::sqrt(rmsx);
432 else
433 rmsx = binwidthx*n/4;
434
435 if (rmsy > 0)
436 rmsy = std::sqrt(rmsy);
437 else
438 rmsy = binwidthy*n/4;
439
440
441 //if the distribution is really gaussian, the best approximation
442 //is binwidx*allcha/(sqrtpi*rmsx)
443 //However, in case of non-gaussian tails, this underestimates
444 //the normalisation constant. In this case the maximum value
445 //is a better approximation.
446 //We take the average of both quantities
447
448 double constant = 0.5 * (valmax+ binwidthx*allcha/(sqrtpi*rmsx))*
449 (valmax+ binwidthy*allcha/(sqrtpi*rmsy));
450
451 f1->SetParameter(0,constant);
452 f1->SetParameter(1,meanx);
453 f1->SetParameter(2,rmsx);
454 f1->SetParLimits(2,0,10*rmsx);
455 f1->SetParameter(3,meany);
456 f1->SetParameter(4,rmsy);
457 f1->SetParLimits(4,0,10*rmsy);
458
459#ifdef DEBUG
460 std::cout << "2D Gaussian initial par values"
461 << constant << " "
462 << meanx << " "
463 << rmsx
464 << meany << " "
465 << rmsy
466 << std::endl;
467#endif
468
469}
470
471// filling fit data from TGraph objects
472
474 // get type of data for TGraph objects
475 double *ex = gr->GetEX();
476 double *ey = gr->GetEY();
477 double * eyl = gr->GetEYlow();
478 double * eyh = gr->GetEYhigh();
479
480
481 // default case for graphs (when they have errors)
483 // if all errors are zero set option of using errors to 1
484 if (fitOpt.fErrors1 || ( ey == 0 && ( eyl == 0 || eyh == 0 ) ) ) {
486 }
487 // need to treat case when all errors are zero
488 // note that by default fitOpt.fCoordError is true
489 else if ( ex != 0 && fitOpt.fCoordErrors) {
490 // check that all errors are not zero
491 int i = 0;
492 while (i < gr->GetN() && type != BinData::kCoordError) {
493 if (ex[i] > 0) type = BinData::kCoordError;
494 ++i;
495 }
496 }
497 // case of asymmetric errors (by default fAsymErrors is true)
498 else if ( ( eyl != 0 && eyh != 0) && fitOpt.fAsymErrors) {
499 // check also if that all errors are non zero's
500 int i = 0;
501 bool zeroErrorX = true;
502 bool zeroErrorY = true;
503 while (i < gr->GetN() && (zeroErrorX || zeroErrorY)) {
504 double e2X = ( gr->GetErrorXlow(i) + gr->GetErrorXhigh(i) );
505 double e2Y = eyl[i] + eyh[i];
506 zeroErrorX &= (e2X <= 0);
507 zeroErrorY &= (e2Y <= 0);
508 ++i;
509 }
510 if (zeroErrorX && zeroErrorY)
512 else if (!zeroErrorX && zeroErrorY)
514 else if (zeroErrorX && !zeroErrorY) {
516 fitOpt.fCoordErrors = false;
517 }
518 else {
520 }
521 }
522
523 // need to look also a case when all errors in y are zero
524 if ( ey != 0 && type != BinData::kCoordError ) {
525 int i = 0;
526 bool zeroError = true;
527 while (i < gr->GetN() && zeroError) {
528 if (ey[i] > 0) zeroError = false;;
529 ++i;
530 }
531 if (zeroError) type = BinData::kNoError;
532 }
533
534
535#ifdef DEBUG
536 std::cout << "type is " << type << " graph type is " << gr->IsA()->GetName() << std::endl;
537#endif
538
539 return type;
540}
541
543 // get type of data for TGraph2D object
544 double *ex = gr->GetEX();
545 double *ey = gr->GetEY();
546 double *ez = gr->GetEZ();
547
548 // default case for graphs (when they have errors)
550 // if all errors are zero set option of using errors to 1
551 if (fitOpt.fErrors1 || ez == 0 ) {
553 }
554 else if ( ex != 0 && ey!=0 && fitOpt.fCoordErrors) {
555 // check that all errors are not zero
556 int i = 0;
557 while (i < gr->GetN() && type != BinData::kCoordError) {
558 if (ex[i] > 0 || ey[i] > 0) type = BinData::kCoordError;
559 ++i;
560 }
561 }
562
563
564#ifdef DEBUG
565 std::cout << "type is " << type << " graph2D type is " << gr->IsA()->GetName() << std::endl;
566#endif
567
568 return type;
569}
570
571
572
573void DoFillData ( BinData & dv, const TGraph * gr, BinData::ErrorType type, TF1 * func ) {
574 // internal method to do the actual filling of the data
575 // given a graph and a multigraph
576
577 // get fit option
578 DataOptions & fitOpt = dv.Opt();
579
580 int nPoints = gr->GetN();
581 double *gx = gr->GetX();
582 double *gy = gr->GetY();
583
584 const DataRange & range = dv.Range();
585 bool useRange = ( range.Size(0) > 0);
586 double xmin = 0;
587 double xmax = 0;
588 range.GetRange(xmin,xmax);
589
590 dv.Initialize(nPoints,1, type);
591
592#ifdef DEBUG
593 std::cout << "DoFillData: graph npoints = " << nPoints << " type " << type << std::endl;
594 if (func) {
595 double a1,a2; func->GetRange(a1,a2); std::cout << "func range " << a1 << " " << a2 << std::endl;
596 }
597#endif
598
599 double x[1];
600 for ( int i = 0; i < nPoints; ++i) {
601
602 x[0] = gx[i];
603
604
605 if (useRange && ( x[0] < xmin || x[0] > xmax) ) continue;
606
607 // need to evaluate function to know about rejected points
608 // hugly but no other solutions
609 if (func) {
610 TF1::RejectPoint(false);
611 (*func)( x ); // evaluate using stored function parameters
612 if (TF1::RejectedPoint() ) continue;
613 }
614
615
616 if (fitOpt.fErrors1)
617 dv.Add( gx[i], gy[i] );
618
619 // for the errors use the getters by index to avoid cases when the arrays are zero
620 // (like in a case of a graph)
621 else if (type == BinData::kValueError) {
622 double errorY = gr->GetErrorY(i);
623 // should consider error = 0 as 1 ? Decide to skip points with zero errors
624 // in case want to keep points with error = 0 as errrors=1 need to set the option UseEmpty
625 if (!HFitInterface::AdjustError(fitOpt,errorY) ) continue;
626 dv.Add( gx[i], gy[i], errorY );
627
628#ifdef DEBUG
629 std::cout << "Point " << i << " " << gx[i] << " " << gy[i] << " " << errorY << std::endl;
630#endif
631
632
633 }
634 else { // case use error in x or asym errors
635 double errorX = 0;
636 if (fitOpt.fCoordErrors)
637 // shoulkd take combined average (sqrt(0.5(e1^2+e2^2)) or math average ?
638 // gr->GetErrorX(i) returns combined average
639 // use math average for same behaviour as before
640 errorX = std::max( 0.5 * ( gr->GetErrorXlow(i) + gr->GetErrorXhigh(i) ) , 0. ) ;
641
642
643 // adjust error in y according to option
644 double errorY = std::max(gr->GetErrorY(i), 0.);
645 // we do not check the return value since we check later if error in X and Y is zero for skipping the point
646 HFitInterface::AdjustError(fitOpt, errorY);
647
648 // skip points with total error = 0
649 if ( errorX <=0 && errorY <= 0 ) continue;
650
651
652 if (type == BinData::kAsymError) {
653 // asymmetric errors
654 dv.Add( gx[i], gy[i], errorX, gr->GetErrorYlow(i), gr->GetErrorYhigh(i) );
655 }
656 else {
657 // case symmetric Y errors
658 dv.Add( gx[i], gy[i], errorX, errorY );
659 }
660 }
661
662 }
663
664#ifdef DEBUG
665 std::cout << "TGraphFitInterface::FillData Graph FitData size is " << dv.Size() << std::endl;
666#endif
667
668}
669
670void FillData(SparseData & dv, const TH1 * h1, TF1 * /*func*/)
671{
672 const int dim = h1->GetDimension();
673 std::vector<double> min(dim);
674 std::vector<double> max(dim);
675
676 int ncells = h1->GetNcells();
677 for ( int i = 0; i < ncells; ++i ) {
678// printf("i: %d; OF: %d; UF: %d; C: %f\n"
679// , i
680// , h1->IsBinOverflow(i) , h1->IsBinUnderflow(i)
681// , h1->GetBinContent(i));
682 if ( !( h1->IsBinOverflow(i) || h1->IsBinUnderflow(i) )
683 && h1->GetBinContent(i))
684 {
685 int x,y,z;
686 h1->GetBinXYZ(i, x, y, z);
687
688// std::cout << "FILLDATA: h1(" << i << ")"
689// << "[" << h1->GetXaxis()->GetBinLowEdge(x) << "-" << h1->GetXaxis()->GetBinUpEdge(x) << "]";
690// if ( dim >= 2 )
691// std::cout << "[" << h1->GetYaxis()->GetBinLowEdge(y) << "-" << h1->GetYaxis()->GetBinUpEdge(y) << "]";
692// if ( dim >= 3 )
693// std::cout << "[" << h1->GetZaxis()->GetBinLowEdge(z) << "-" << h1->GetZaxis()->GetBinUpEdge(z) << "]";
694
695// std::cout << h1->GetBinContent(i) << std::endl;
696
697 min[0] = h1->GetXaxis()->GetBinLowEdge(x);
698 max[0] = h1->GetXaxis()->GetBinUpEdge(x);
699 if ( dim >= 2 )
700 {
701 min[1] = h1->GetYaxis()->GetBinLowEdge(y);
702 max[1] = h1->GetYaxis()->GetBinUpEdge(y);
703 }
704 if ( dim >= 3 ) {
705 min[2] = h1->GetZaxis()->GetBinLowEdge(z);
706 max[2] = h1->GetZaxis()->GetBinUpEdge(z);
707 }
708
709 dv.Add(min, max, h1->GetBinContent(i), h1->GetBinError(i));
710 }
711 }
712}
713
714void FillData(SparseData & dv, const THnBase * h1, TF1 * /*func*/)
715{
716 const int dim = h1->GetNdimensions();
717 std::vector<double> min(dim);
718 std::vector<double> max(dim);
719 std::vector<Int_t> coord(dim);
720
721 ULong64_t nEntries = h1->GetNbins();
722 for ( ULong64_t i = 0; i < nEntries; i++ )
723 {
724 double value = h1->GetBinContent( i, &coord[0] );
725 if ( !value ) continue;
726
727// std::cout << "FILLDATA(SparseData): h1(" << i << ")";
728
729 // Exclude underflows and oveflows! (defect behaviour with the TH1*)
730 bool insertBox = true;
731 for ( int j = 0; j < dim && insertBox; ++j )
732 {
733 TAxis* axis = h1->GetAxis(j);
734 if ( ( axis->GetBinLowEdge(coord[j]) < axis->GetXmin() ) ||
735 ( axis->GetBinUpEdge(coord[j]) > axis->GetXmax() ) ) {
736 insertBox = false;
737 }
738 min[j] = h1->GetAxis(j)->GetBinLowEdge(coord[j]);
739 max[j] = h1->GetAxis(j)->GetBinUpEdge(coord[j]);
740 }
741 if ( !insertBox ) {
742// std::cout << "NOT INSERTED!"<< std::endl;
743 continue;
744 }
745
746// for ( int j = 0; j < dim; ++j )
747// {
748// std::cout << "[" << h1->GetAxis(j)->GetBinLowEdge(coord[j])
749// << "-" << h1->GetAxis(j)->GetBinUpEdge(coord[j]) << "]";
750// }
751// std::cout << h1->GetBinContent(i) << std::endl;
752
753 dv.Add(min, max, value, h1->GetBinError(i));
754 }
755}
756
757void FillData(BinData & dv, const THnBase * s1, TF1 * func)
758{
759 // Fill the Range of the THnBase
760 unsigned int const ndim = s1->GetNdimensions();
761 std::vector<double> xmin(ndim);
762 std::vector<double> xmax(ndim);
763 for ( unsigned int i = 0; i < ndim; ++i ) {
764 TAxis* axis = s1->GetAxis(i);
765 xmin[i] = axis->GetXmin();
766 xmax[i] = axis->GetXmax();
767 }
768
769 // Put default options, needed for the likelihood fitting of sparse
770 // data.
771 ROOT::Fit::DataOptions& dopt = dv.Opt();
772 dopt.fUseEmpty = true;
773 // when using sparse data need to set option to use normalized bin volume, because sparse bins are merged together
774 //if (!dopt.fIntegral) dopt.fBinVolume = true;
775 dopt.fBinVolume = true;
776 dopt.fNormBinVolume = true;
777
778 // Get the sparse data
779 ROOT::Fit::SparseData d(ndim, &xmin[0], &xmax[0]);
780 ROOT::Fit::FillData(d, s1, func);
781
782// std::cout << "FillData(BinData & dv, const THnBase * s1, TF1 * func) (1)" << std::endl;
783
784 // Create the bin data from the sparse data
785 d.GetBinDataIntegral(dv);
786
787}
788
789void FillData ( BinData & dv, const TGraph * gr, TF1 * func ) {
790 // fill the data vector from a TGraph. Pass also the TF1 function which is
791 // needed in case to exclude points rejected by the function
792 assert(gr != 0);
793
794 // get fit option
795 DataOptions & fitOpt = dv.Opt();
796
798 // adjust option according to type
799 fitOpt.fErrors1 = (type == BinData::kNoError);
800 // set this if we want to have error=1 for points with zero errors (by default they are skipped)
801 // fitOpt.fUseEmpty = true;
802
803 // use coordinate or asym errors in case option is set and type is consistent
805 fitOpt.fAsymErrors &= (type == BinData::kAsymError);
806
807
808 // if data are filled already check if there are consistent - otherwise do nothing
809 if (dv.Size() > 0 && dv.NDim() == 1 ) {
810 // check if size is correct otherwise flag an errors
811 if ( dv.GetErrorType() != type ) {
812 Error("FillData","Inconsistent TGraph with previous data set- skip all graph data");
813 return;
814 }
815 }
816
817 DoFillData(dv, gr, type, func);
818
819}
820
821void FillData ( BinData & dv, const TMultiGraph * mg, TF1 * func ) {
822 // fill the data vector from a TMultiGraph. Pass also the TF1 function which is
823 // needed in case to exclude points rejected by the function
824 assert(mg != 0);
825
826 TList * grList = mg->GetListOfGraphs();
827 assert(grList != 0);
828
829#ifdef DEBUG
830// grList->Print();
831 TIter itr(grList, kIterBackward);
832 TObject *obj;
833 std::cout << "multi-graph list of graps: " << std::endl;
834 while ((obj = itr())) {
835 std::cout << obj->IsA()->GetName() << std::endl;
836 }
837
838#endif
839
840 // get fit option
841 DataOptions & fitOpt = dv.Opt();
842
843 // loop on the graphs to get the data type (use maximum)
844 TIter next(grList);
845
847 TGraph *gr = 0;
848 while ((gr = (TGraph*) next())) {
850 if (t > type ) type = t;
851 }
852 // adjust option according to type
853 fitOpt.fErrors1 = (type == BinData::kNoError);
856
857
858#ifdef DEBUG
859 std::cout << "Fitting MultiGraph of type " << type << std::endl;
860#endif
861
862 // fill the data now
863 next = grList;
864 while ((gr = (TGraph*) next())) {
865 DoFillData( dv, gr, type, func);
866 }
867
868#ifdef DEBUG
869 std::cout << "TGraphFitInterface::FillData MultiGraph FitData size is " << dv.Size() << std::endl;
870#endif
871
872}
873
874void FillData ( BinData & dv, const TGraph2D * gr, TF1 * func ) {
875 // fill the data vector from a TGraph2D. Pass also the TF1 function which is
876 // needed in case to exclude points rejected by the function
877 // in case of a pure TGraph
878 assert(gr != 0);
879
880 // get fit option
881 DataOptions & fitOpt = dv.Opt();
883 // adjust option according to type
884 fitOpt.fErrors1 = (type == BinData::kNoError);
886 fitOpt.fAsymErrors = false; // a TGraph2D with asymmetric errors does not exist
887
888 int nPoints = gr->GetN();
889 double *gx = gr->GetX();
890 double *gy = gr->GetY();
891 double *gz = gr->GetZ();
892
893 // if all errors are zero set option of using errors to 1
894 if ( gr->GetEZ() == 0) fitOpt.fErrors1 = true;
895
896 double x[2];
897 double ex[2];
898
899 // look at data range
900 const DataRange & range = dv.Range();
901 bool useRangeX = ( range.Size(0) > 0);
902 bool useRangeY = ( range.Size(1) > 0);
903 double xmin = 0;
904 double xmax = 0;
905 double ymin = 0;
906 double ymax = 0;
907 range.GetRange(xmin,xmax,ymin,ymax);
908
909 dv.Initialize(nPoints,2, type);
910
911 for ( int i = 0; i < nPoints; ++i) {
912
913 x[0] = gx[i];
914 x[1] = gy[i];
915
916 //if (fitOpt.fUseRange && HFitInterface::IsPointOutOfRange(func, x) ) continue;
917 if (useRangeX && ( x[0] < xmin || x[0] > xmax) ) continue;
918 if (useRangeY && ( x[1] < ymin || x[1] > ymax) ) continue;
919
920 // need to evaluate function to know about rejected points
921 // hugly but no other solutions
922 if (func) {
923 TF1::RejectPoint(false);
924 (*func)( x ); // evaluate using stored function parameters
925 if (TF1::RejectedPoint() ) continue;
926 }
927
928 if (type == BinData::kNoError) {
929 dv.Add( x, gz[i] );
930 continue;
931 }
932
933 double errorZ = gr->GetErrorZ(i);
934 if (!HFitInterface::AdjustError(fitOpt,errorZ) ) continue;
935
936 if (type == BinData::kValueError) {
937 dv.Add( x, gz[i], errorZ );
938 }
939 else if (type == BinData::kCoordError) { // case use error in coordinates (x and y)
940 ex[0] = std::max(gr->GetErrorX(i), 0.);
941 ex[1] = std::max(gr->GetErrorY(i), 0.);
942 dv.Add( x, gz[i], ex, errorZ );
943 }
944 else
945 assert(0); // should not go here
946
947#ifdef DEBUG
948 std::cout << "Point " << i << " " << gx[i] << " " << gy[i] << " " << errorZ << std::endl;
949#endif
950
951 }
952
953#ifdef DEBUG
954 std::cout << "THFitInterface::FillData Graph2D FitData size is " << dv.Size() << std::endl;
955#endif
956
957}
958
959
960// confidence intervals
961bool GetConfidenceIntervals(const TH1 * h1, const ROOT::Fit::FitResult & result, TGraphErrors * gr, double cl ) {
962 if (h1->GetDimension() != 1) {
963 Error("GetConfidenceIntervals","Invalid object used for storing confidence intervals");
964 return false;
965 }
966 // fill fit data sets with points to estimate cl.
967 BinData d;
968 FillData(d,h1,0);
969 gr->Set(d.NPoints() );
970 double * ci = gr->GetEY(); // make CL values error of the graph
971 result.GetConfidenceIntervals(d,ci,cl);
972 // put function value as abscissa of the graph
973 for (unsigned int ipoint = 0; ipoint < d.NPoints(); ++ipoint) {
974 const double * x = d.Coords(ipoint);
975 const ROOT::Math::IParamMultiFunction * func = result.FittedFunction();
976 gr->SetPoint(ipoint, x[0], (*func)(x) );
977 }
978 return true;
979}
980
981
982} // end namespace Fit
983
984} // end namespace ROOT
985
#define d(i)
Definition: RSha256.hxx:102
#define s1(x)
Definition: RSha256.hxx:91
unsigned long long ULong64_t
Definition: RtypesCore.h:70
const Bool_t kIterBackward
Definition: TCollection.h:41
void Error(const char *location, const char *msgfmt,...)
void Warning(const char *location, const char *msgfmt,...)
int type
Definition: TGX11.cxx:120
float xmin
Definition: THbookFile.cxx:93
float ymin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
float ymax
Definition: THbookFile.cxx:93
double sqrt(double)
double log(double)
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:53
void AddBinUpEdge(const double *xup)
add the bin width data, a pointer to an array with the bin upper edge information.
Definition: BinData.cxx:627
ErrorType GetErrorType() const
retrieve the errortype
Definition: BinData.h:545
void Add(double x, double y)
add one dim data with only coordinate and values
Definition: BinData.cxx:422
void Initialize(unsigned int newPoints, unsigned int dim=1, ErrorType err=kValueError)
Definition: BinData.cxx:349
class describing the range in the coordinates it supports multiple range in a coordinate.
Definition: DataRange.h:34
unsigned int Size(unsigned int icoord=0) const
return range size for coordinate icoord (starts from zero) Size == 0 indicates no range is present [-...
Definition: DataRange.h:70
void GetRange(unsigned int irange, unsigned int icoord, double &xmin, double &xmax) const
get the i-th range for given coordinate.
Definition: DataRange.h:103
unsigned int Size() const
return number of fit points
Definition: FitData.h:303
unsigned int NDim() const
return coordinate data dimension
Definition: FitData.h:311
const DataOptions & Opt() const
access to options
Definition: FitData.h:319
const DataRange & Range() const
access to range
Definition: FitData.h:331
class containg the result of the fit and all the related information (fitted parameter values,...
Definition: FitResult.h:48
void GetConfidenceIntervals(unsigned int n, unsigned int stride1, unsigned int stride2, const double *x, double *ci, double cl=0.95, bool norm=false) const
get confidence intervals for an array of n points x.
Definition: FitResult.cxx:560
const IModelFunction * FittedFunction() const
fitting quantities
Definition: FitResult.h:149
void Add(std::vector< double > &min, std::vector< double > &max, const double content, const double error=1.0)
Definition: SparseData.cxx:231
Class to manage histogram axis.
Definition: TAxis.h:30
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:464
Double_t GetXmax() const
Definition: TAxis.h:134
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
Definition: TAxis.cxx:504
virtual Int_t FindFixBin(Double_t x) const
Find bin number corresponding to abscissa x.
Definition: TAxis.cxx:405
Int_t GetLast() const
Return last bin on the axis i.e.
Definition: TAxis.cxx:455
Double_t GetXmin() const
Definition: TAxis.h:133
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width.
Definition: TAxis.cxx:526
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition: TAxis.cxx:514
Int_t GetFirst() const
Return first bin on the axis i.e.
Definition: TAxis.cxx:444
1-Dim function class
Definition: TF1.h:211
static void RejectPoint(Bool_t reject=kTRUE)
Static function to set the global flag to reject points the fgRejectPoint global flag is tested by al...
Definition: TF1.cxx:3644
virtual void GetRange(Double_t *xmin, Double_t *xmax) const
Return range of a generic N-D function.
Definition: TF1.cxx:2257
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set limits for parameter ipar.
Definition: TF1.cxx:3491
static Bool_t RejectedPoint()
See TF1::RejectPoint above.
Definition: TF1.cxx:3653
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:628
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:618
virtual Bool_t IsInside(const Double_t *x) const
return kTRUE if the point is inside the function range
Definition: TF1.h:582
virtual Int_t GetNdim() const
Definition: TF1.h:469
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition: TGraph2D.h:40
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:26
Double_t GetErrorXhigh(Int_t bin) const
This function is called by GraphFitChisquare.
Double_t GetErrorXlow(Int_t bin) const
This function is called by GraphFitChisquare.
Double_t GetErrorYlow(Int_t bin) const
This function is called by GraphFitChisquare.
Double_t GetErrorX(Int_t bin) const
This function is called by GraphFitChisquare.
Double_t GetErrorY(Int_t bin) const
This function is called by GraphFitChisquare.
Double_t GetErrorYhigh(Int_t bin) const
This function is called by GraphFitChisquare.
Double_t * GetEY() const
Definition: TGraphErrors.h:67
Double_t * GetEX() const
Definition: TGraphErrors.h:66
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2200
Double_t * GetY() const
Definition: TGraph.h:131
virtual Double_t * GetEYlow() const
Definition: TGraph.h:137
Int_t GetN() const
Definition: TGraph.h:123
Double_t * GetX() const
Definition: TGraph.h:130
virtual Double_t * GetEYhigh() const
Definition: TGraph.h:136
virtual void Set(Int_t n)
Set number of points in the graph Existing coordinates are preserved New coordinates above fNpoints a...
Definition: TGraph.cxx:2135
The TH1 histogram class.
Definition: TH1.h:56
TAxis * GetZaxis()
Definition: TH1.h:318
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8384
virtual Int_t GetDimension() const
Definition: TH1.h:278
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
virtual void GetBinXYZ(Int_t binglobal, Int_t &binx, Int_t &biny, Int_t &binz) const
Return binx, biny, binz corresponding to the global bin number globalbin see TH1::GetBin function abo...
Definition: TH1.cxx:4705
virtual Int_t GetNcells() const
Definition: TH1.h:295
TAxis * GetYaxis()
Definition: TH1.h:317
Bool_t IsBinUnderflow(Int_t bin, Int_t axis=0) const
Return true if the bin is underflow.
Definition: TH1.cxx:4943
Bool_t IsBinOverflow(Int_t bin, Int_t axis=0) const
Return true if the bin is overflow.
Definition: TH1.cxx:4911
virtual Double_t GetBinLowEdge(Int_t bin) const
Return bin lower edge for 1D histogram.
Definition: TH1.cxx:8473
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4790
Multidimensional histogram base.
Definition: THnBase.h:43
A doubly linked list.
Definition: TList.h:44
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition: TMultiGraph.h:35
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Mother of all ROOT objects.
Definition: TObject.h:37
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:357
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Double_t ey[n]
Definition: legend1.C:17
TGraphErrors * gr
Definition: legend1.C:25
Double_t ex[n]
Definition: legend1.C:17
TH1F * h1
Definition: legend1.C:5
TF1 * f1
Definition: legend1.C:11
TFitResultPtr Fit(FitObject *h1, TF1 *f1, Foption_t &option, const ROOT::Math::MinimizerOptions &moption, const char *goption, ROOT::Fit::DataRange &range)
Definition: HFitImpl.cxx:134
void ExamineRange(const TAxis *axis, std::pair< double, double > range, int &hxfirst, int &hxlast)
bool AdjustError(const DataOptions &option, double &error, double value=1)
bool IsPointOutOfRange(const TF1 *func, const double *x)
void Init2DGaus(const ROOT::Fit::BinData &data, TF1 *f1)
compute initial parameter for 2D gaussian function given the fit data Set the sigma limits for zero t...
void FillData(BinData &dv, const TH1 *hist, TF1 *func=0)
fill the data vector from a TH1.
void InitExpo(const ROOT::Fit::BinData &data, TF1 *f1)
compute initial parameter for an exponential function given the fit data Set the constant and slope a...
void InitGaus(const ROOT::Fit::BinData &data, TF1 *f1)
compute initial parameter for gaussian function given the fit data Set the sigma limits for zero top ...
void DoFillData(BinData &dv, const TGraph *gr, BinData::ErrorType type, TF1 *func)
BinData::ErrorType GetDataType(const TGraph *gr, DataOptions &fitOpt)
bool GetConfidenceIntervals(const TH1 *h1, const ROOT::Fit::FitResult &r, TGraphErrors *gr, double cl=0.95)
compute confidence intervals at level cl for a fitted histogram h1 in a TGraphErrors gr
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
static constexpr double mg
static constexpr double s
DataOptions : simple structure holding the options on how the data are filled.
Definition: DataOptions.h:28