Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 ==nullptr) 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 != nullptr);
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 !=nullptr && func->GetNdim() == hdim-1) ndim = hdim-1;
177
178 assert( ndim > 0 );
179 //typedef BinPoint::CoordData CoordData;
180 //CoordData x = CoordData( hfit->GetDimension() );
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 != nullptr) {
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 the number 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 if (fitOpt.fErrors1)
240 dv.Add( x, value );
241 else
242 dv.Add( x, value, error );
243 if (useBinEdges) {
244 dv.AddBinUpEdge( s );
245 }
246 }
247
248
249#ifdef DEBUG
250 std::cout << "bin " << binx << " add point " << x[0] << " " << hfit->GetBinContent(binx) << std::endl;
251#endif
252
253 } // end loop on z bins
254 } // end loop on y bins
255 } // end loop on x axis
256
257
258#ifdef DEBUG
259 std::cout << "THFitInterface::FillData: Hist FitData size is " << dv.Size() << std::endl;
260#endif
261
262}
263
264////////////////////////////////////////////////////////////////////////////////
265/// Compute rough values of parameters for an exponential
266
268{
269 unsigned int n = data.Size();
270 if (n == 0) return;
271
272 // find xmin and xmax of the data
273 double valxmin;
274 double xmin = *(data.GetPoint(0,valxmin));
275 double xmax = xmin;
276 double valxmax = valxmin;
277
278 for (unsigned int i = 1; i < n; ++ i) {
279 double val;
280 double x = *(data.GetPoint(i,val) );
281 if (x < xmin) {
282 xmin = x;
283 valxmin = val;
284 }
285 else if (x > xmax) {
286 xmax = x;
287 valxmax = val;
288 }
289 }
290
291 // avoid negative values of valxmin/valxmax
292 if (valxmin <= 0 && valxmax > 0 ) valxmin = valxmax;
293 else if (valxmax <=0 && valxmin > 0) valxmax = valxmin;
294 else if (valxmin <=0 && valxmax <= 0) { valxmin = 1; valxmax = 1; }
295
296 double slope = std::log( valxmax/valxmin) / (xmax - xmin);
297 double constant = std::log(valxmin) - slope * xmin;
298 f1->SetParameters(constant, slope);
299}
300
301
302////////////////////////////////////////////////////////////////////////////////
303/// Compute Initial values of parameters for a gaussian
304/// derived from function H1InitGaus defined in TH1.cxx
305
307{
308
309 static const double sqrtpi = 2.506628;
310
311 // - Compute mean value and RMS of the data
312 unsigned int n = data.Size();
313 if (n == 0) return;
314 double sumx = 0;
315 double sumx2 = 0;
316 double allcha = 0;
317 double valmax = 0;
318 double rangex = data.Coords(n-1)[0] - data.Coords(0)[0];
319 // to avoid binwidth = 0 set arbitrarly to 1
320 double binwidth = 1;
321 if ( rangex > 0) binwidth = rangex;
322 double x0 = 0;
323 for (unsigned int i = 0; i < n; ++ i) {
324 double val;
325 double x = *(data.GetPoint(i,val) );
326 sumx += val*x;
327 sumx2 += val*x*x;
328 allcha += val;
329 if (val > valmax) valmax = val;
330 if (i > 0) {
331 double dx = x - x0;
332 if (dx < binwidth) binwidth = dx;
333 }
334 x0 = x;
335 }
336
337 if (allcha <= 0) return;
338 double mean = sumx/allcha;
339 double rms = sumx2/allcha - mean*mean;
340
341
342 if (rms > 0)
343 rms = std::sqrt(rms);
344 else
345 rms = binwidth*n/4;
346
347
348 //if the distribution is really gaussian, the best approximation
349 //is binwidx*allcha/(sqrtpi*rms)
350 //However, in case of non-gaussian tails, this underestimates
351 //the normalisation constant. In this case the maximum value
352 //is a better approximation.
353 //We take the average of both quantities
354
355// printf("valmax %f other %f bw %f allcha %f rms %f \n",valmax, binwidth*allcha/(sqrtpi*rms),
356// binwidth, allcha,rms );
357
358 double constant = 0.5*(valmax+ binwidth*allcha/(sqrtpi*rms));
359
360
361 //In case the mean value is outside the histo limits and
362 //the RMS is bigger than the range, we take
363 // mean = center of bins
364 // rms = half range
365// Double_t xmin = curHist->GetXaxis()->GetXmin();
366// Double_t xmax = curHist->GetXaxis()->GetXmax();
367// if ((mean < xmin || mean > xmax) && rms > (xmax-xmin)) {
368// mean = 0.5*(xmax+xmin);
369// rms = 0.5*(xmax-xmin);
370// }
371
372 f1->SetParameter(0,constant);
373 f1->SetParameter(1,mean);
374 f1->SetParameter(2,rms);
375 f1->SetParLimits(2,0,10*rms);
376
377
378#ifdef DEBUG
379 std::cout << "Gaussian initial par values" << constant << " " << mean << " " << rms << std::endl;
380#endif
381
382}
383
384////////////////////////////////////////////////////////////////////////////////
385/// Compute Initial values of parameters for a gaussian
386/// derived from function H1InitGaus defined in TH1.cxx
387
389{
390
391 static const double sqrtpi = 2.506628;
392
393 // - Compute mean value and RMS of the data
394 unsigned int n = data.Size();
395 if (n == 0) return;
396 double sumx = 0, sumy = 0;
397 double sumx2 = 0, sumy2 = 0;
398 double allcha = 0;
399 double valmax = 0;
400 double rangex = data.Coords(n-1)[0] - data.Coords(0)[0];
401 double rangey = data.Coords(n-1)[1] - data.Coords(0)[1];
402 // to avoid binwidthx = 0 set arbitrarly to 1
403 double binwidthx = 1, binwidthy = 1;
404 if ( rangex > 0) binwidthx = rangex;
405 if ( rangey > 0) binwidthy = rangey;
406 double x0 = 0, y0 = 0;
407 for (unsigned int i = 0; i < n; ++i) {
408 double val;
409 const double *coords = data.GetPoint(i,val);
410 double x = coords[0], y = coords[1];
411 sumx += val*x;
412 sumy += val*y;
413 sumx2 += val*x*x;
414 sumy2 += val*y*y;
415 allcha += val;
416 if (val > valmax) valmax = val;
417 if (i > 0) {
418 double dx = x - x0;
419 if (dx < binwidthx) binwidthx = dx;
420 double dy = y - y0;
421 if (dy < binwidthy) binwidthy = dy;
422 }
423 x0 = x;
424 y0 = y;
425 }
426
427 if (allcha <= 0) return;
428 double meanx = sumx/allcha, meany = sumy/allcha;
429 double rmsx = sumx2/allcha - meanx*meanx;
430 double rmsy = sumy2/allcha - meany*meany;
431
432
433 if (rmsx > 0)
434 rmsx = std::sqrt(rmsx);
435 else
436 rmsx = binwidthx*n/4;
437
438 if (rmsy > 0)
439 rmsy = std::sqrt(rmsy);
440 else
441 rmsy = binwidthy*n/4;
442
443
444 //if the distribution is really gaussian, the best approximation
445 //is binwidx*allcha/(sqrtpi*rmsx)
446 //However, in case of non-gaussian tails, this underestimates
447 //the normalisation constant. In this case the maximum value
448 //is a better approximation.
449 //We take the average of both quantities
450
451 double constant = 0.5 * (valmax+ binwidthx*allcha/(sqrtpi*rmsx))*
452 (valmax+ binwidthy*allcha/(sqrtpi*rmsy));
453
454 f1->SetParameter(0,constant);
455 f1->SetParameter(1,meanx);
456 f1->SetParameter(2,rmsx);
457 f1->SetParLimits(2,0,10*rmsx);
458 f1->SetParameter(3,meany);
459 f1->SetParameter(4,rmsy);
460 f1->SetParLimits(4,0,10*rmsy);
461
462#ifdef DEBUG
463 std::cout << "2D Gaussian initial par values"
464 << constant << " "
465 << meanx << " "
466 << rmsx
467 << meany << " "
468 << rmsy
469 << std::endl;
470#endif
471
472}
473
474// filling fit data from TGraph objects
475
477 // get type of data for TGraph objects
478 double *ex = gr->GetEX();
479 double *ey = gr->GetEY();
480 double * eyl = gr->GetEYlow();
481 double * eyh = gr->GetEYhigh();
482
483
484 // default case for graphs (when they have errors)
486 // if all errors are zero set option of using errors to 1
487 if (fitOpt.fErrors1 || ( ey == nullptr && ( eyl == nullptr || eyh == nullptr ) ) ) {
489 }
490 // need to treat case when all errors are zero
491 // note that by default fitOpt.fCoordError is true
492 else if ( ex != nullptr && fitOpt.fCoordErrors) {
493 // check that all errors are not zero
494 int i = 0;
495 while (i < gr->GetN() && type != BinData::kCoordError) {
496 if (ex[i] > 0) type = BinData::kCoordError;
497 ++i;
498 }
499 }
500 // case of asymmetric errors (by default fAsymErrors is true)
501 else if ( ( eyl != nullptr && eyh != nullptr) && fitOpt.fAsymErrors) {
502 // check also if that all errors are non zero's
503 int i = 0;
504 bool zeroErrorX = true;
505 bool zeroErrorY = true;
506 while (i < gr->GetN() && (zeroErrorX || zeroErrorY)) {
507 double e2X = ( gr->GetErrorXlow(i) + gr->GetErrorXhigh(i) );
508 double e2Y = eyl[i] + eyh[i];
509 zeroErrorX &= (e2X <= 0);
510 zeroErrorY &= (e2Y <= 0);
511 ++i;
512 }
513 if (zeroErrorX && zeroErrorY)
515 else if (!zeroErrorX && zeroErrorY)
517 else if (zeroErrorX && !zeroErrorY) {
519 fitOpt.fCoordErrors = false;
520 }
521 else {
523 }
524 }
525
526 // need to look also a case when all errors in y are zero
527 if ( ey != nullptr && type != BinData::kCoordError ) {
528 int i = 0;
529 bool zeroError = true;
530 while (i < gr->GetN() && zeroError) {
531 if (ey[i] > 0) zeroError = false;;
532 ++i;
533 }
534 if (zeroError) type = BinData::kNoError;
535 }
536
537
538#ifdef DEBUG
539 std::cout << "type is " << type << " graph type is " << gr->IsA()->GetName() << std::endl;
540#endif
541
542 return type;
543}
544
546 // get type of data for TGraph2D object
547 double *ex = gr->GetEX();
548 double *ey = gr->GetEY();
549 double *ez = gr->GetEZ();
550
551 // default case for graphs (when they have errors)
553 // if all errors are zero set option of using errors to 1
554 if (fitOpt.fErrors1 || ez == nullptr ) {
556 }
557 else if ( ex != nullptr && ey!=nullptr && fitOpt.fCoordErrors) {
558 // check that all errors are not zero
559 int i = 0;
560 while (i < gr->GetN() && type != BinData::kCoordError) {
561 if (ex[i] > 0 || ey[i] > 0) type = BinData::kCoordError;
562 ++i;
563 }
564 }
565
566
567#ifdef DEBUG
568 std::cout << "type is " << type << " graph2D type is " << gr->IsA()->GetName() << std::endl;
569#endif
570
571 return type;
572}
573
574
575
576void DoFillData ( BinData & dv, const TGraph * gr, BinData::ErrorType type, TF1 * func ) {
577 // internal method to do the actual filling of the data
578 // given a graph and a multigraph
579
580 // get fit option
581 DataOptions & fitOpt = dv.Opt();
582
583 int nPoints = gr->GetN();
584 double *gx = gr->GetX();
585 double *gy = gr->GetY();
586
587 const DataRange & range = dv.Range();
588 bool useRange = ( range.Size(0) > 0);
589 double xmin = 0;
590 double xmax = 0;
591 range.GetRange(xmin,xmax);
592
593 dv.Initialize(nPoints,1, type);
594
595#ifdef DEBUG
596 std::cout << "DoFillData: graph npoints = " << nPoints << " type " << type << std::endl;
597 if (func) {
598 double a1,a2; func->GetRange(a1,a2); std::cout << "func range " << a1 << " " << a2 << std::endl;
599 }
600#endif
601
602 double x[1];
603 for ( int i = 0; i < nPoints; ++i) {
604
605 x[0] = gx[i];
606
607
608 if (useRange && ( x[0] < xmin || x[0] > xmax) ) continue;
609
610 // need to evaluate function to know about rejected points
611 // hugly but no other solutions
612 if (func) {
613 TF1::RejectPoint(false);
614 (*func)( x ); // evaluate using stored function parameters
615 if (TF1::RejectedPoint() ) continue;
616 }
617
618
619 if (fitOpt.fErrors1)
620 dv.Add( gx[i], gy[i] );
621
622 // for the errors use the getters by index to avoid cases when the arrays are zero
623 // (like in a case of a graph)
624 else if (type == BinData::kValueError) {
625 double errorY = gr->GetErrorY(i);
626 // should consider error = 0 as 1 ? Decide to skip points with zero errors
627 // in case want to keep points with error = 0 as errrors=1 need to set the option UseEmpty
628 if (!HFitInterface::AdjustError(fitOpt,errorY) ) continue;
629 dv.Add( gx[i], gy[i], errorY );
630
631#ifdef DEBUG
632 std::cout << "Point " << i << " " << gx[i] << " " << gy[i] << " " << errorY << std::endl;
633#endif
634
635
636 }
637 else { // case use error in x or asym errors
638 double errorX = 0;
639 if (fitOpt.fCoordErrors)
640 // shoulkd take combined average (sqrt(0.5(e1^2+e2^2)) or math average ?
641 // gr->GetErrorX(i) returns combined average
642 // use math average for same behaviour as before
643 errorX = std::max( 0.5 * ( gr->GetErrorXlow(i) + gr->GetErrorXhigh(i) ) , 0. ) ;
644
645
646 // adjust error in y according to option
647 double errorY = std::max(gr->GetErrorY(i), 0.);
648 // we do not check the return value since we check later if error in X and Y is zero for skipping the point
649 HFitInterface::AdjustError(fitOpt, errorY);
650
651 // skip points with total error = 0
652 if ( errorX <=0 && errorY <= 0 ) continue;
653
654
655 if (type == BinData::kAsymError) {
656 // asymmetric errors
657 dv.Add( gx[i], gy[i], errorX, gr->GetErrorYlow(i), gr->GetErrorYhigh(i) );
658 }
659 else {
660 // case symmetric Y errors
661 dv.Add( gx[i], gy[i], errorX, errorY );
662 }
663 }
664
665 }
666
667#ifdef DEBUG
668 std::cout << "TGraphFitInterface::FillData Graph FitData size is " << dv.Size() << std::endl;
669#endif
670
671}
672
673void FillData(SparseData & dv, const TH1 * h1, TF1 * /*func*/)
674{
675 const int dim = h1->GetDimension();
676 std::vector<double> min(dim);
677 std::vector<double> max(dim);
678
679 int ncells = h1->GetNcells();
680 for ( int i = 0; i < ncells; ++i ) {
681// printf("i: %d; OF: %d; UF: %d; C: %f\n"
682// , i
683// , h1->IsBinOverflow(i) , h1->IsBinUnderflow(i)
684// , h1->GetBinContent(i));
685 if ( !( h1->IsBinOverflow(i) || h1->IsBinUnderflow(i) )
686 && h1->GetBinContent(i))
687 {
688 int x,y,z;
689 h1->GetBinXYZ(i, x, y, z);
690
691// std::cout << "FILLDATA: h1(" << i << ")"
692// << "[" << h1->GetXaxis()->GetBinLowEdge(x) << "-" << h1->GetXaxis()->GetBinUpEdge(x) << "]";
693// if ( dim >= 2 )
694// std::cout << "[" << h1->GetYaxis()->GetBinLowEdge(y) << "-" << h1->GetYaxis()->GetBinUpEdge(y) << "]";
695// if ( dim >= 3 )
696// std::cout << "[" << h1->GetZaxis()->GetBinLowEdge(z) << "-" << h1->GetZaxis()->GetBinUpEdge(z) << "]";
697
698// std::cout << h1->GetBinContent(i) << std::endl;
699
700 min[0] = h1->GetXaxis()->GetBinLowEdge(x);
701 max[0] = h1->GetXaxis()->GetBinUpEdge(x);
702 if ( dim >= 2 )
703 {
704 min[1] = h1->GetYaxis()->GetBinLowEdge(y);
705 max[1] = h1->GetYaxis()->GetBinUpEdge(y);
706 }
707 if ( dim >= 3 ) {
708 min[2] = h1->GetZaxis()->GetBinLowEdge(z);
709 max[2] = h1->GetZaxis()->GetBinUpEdge(z);
710 }
711
712 dv.Add(min, max, h1->GetBinContent(i), h1->GetBinError(i));
713 }
714 }
715}
716
717void FillData(SparseData & dv, const THnBase * h1, TF1 * /*func*/)
718{
719 const int dim = h1->GetNdimensions();
720 std::vector<double> min(dim);
721 std::vector<double> max(dim);
722 std::vector<Int_t> coord(dim);
723
724 ULong64_t nEntries = h1->GetNbins();
725 for ( ULong64_t i = 0; i < nEntries; i++ )
726 {
727 double value = h1->GetBinContent( i, &coord[0] );
728 if ( !value ) continue;
729
730// std::cout << "FILLDATA(SparseData): h1(" << i << ")";
731
732 // Exclude underflows and overflows! (defect behaviour with the TH1*)
733 bool insertBox = true;
734 for ( int j = 0; j < dim && insertBox; ++j )
735 {
736 TAxis* axis = h1->GetAxis(j);
737 if ( ( axis->GetBinLowEdge(coord[j]) < axis->GetXmin() ) ||
738 ( axis->GetBinUpEdge(coord[j]) > axis->GetXmax() ) ) {
739 insertBox = false;
740 }
741 min[j] = h1->GetAxis(j)->GetBinLowEdge(coord[j]);
742 max[j] = h1->GetAxis(j)->GetBinUpEdge(coord[j]);
743 }
744 if ( !insertBox ) {
745// std::cout << "NOT INSERTED!"<< std::endl;
746 continue;
747 }
748
749// for ( int j = 0; j < dim; ++j )
750// {
751// std::cout << "[" << h1->GetAxis(j)->GetBinLowEdge(coord[j])
752// << "-" << h1->GetAxis(j)->GetBinUpEdge(coord[j]) << "]";
753// }
754// std::cout << h1->GetBinContent(i) << std::endl;
755
756 dv.Add(min, max, value, h1->GetBinError(i));
757 }
758}
759
760void FillData(BinData & dv, const THnBase * s1, TF1 * func)
761{
762 // Fill the Range of the THnBase
763 unsigned int const ndim = s1->GetNdimensions();
764 std::vector<double> xmin(ndim);
765 std::vector<double> xmax(ndim);
766 for ( unsigned int i = 0; i < ndim; ++i ) {
767 TAxis* axis = s1->GetAxis(i);
768 xmin[i] = axis->GetXmin();
769 xmax[i] = axis->GetXmax();
770 }
771
772 // Put default options, needed for the likelihood fitting of sparse
773 // data.
774 ROOT::Fit::DataOptions& dopt = dv.Opt();
775 //dopt.fUseEmpty = true;
776 // when using sparse data need to set option to use normalized bin volume, because sparse bins are merged together
777 //if (!dopt.fIntegral) dopt.fBinVolume = true;
778 dopt.fBinVolume = true;
779 dopt.fNormBinVolume = true;
780
781 // Get the sparse data
782 ROOT::Fit::SparseData d(ndim, &xmin[0], &xmax[0]);
783 ROOT::Fit::FillData(d, s1, func);
784
785// std::cout << "FillData(BinData & dv, const THnBase * s1, TF1 * func) (1)" << std::endl;
786
787 // Create the bin data from the sparse data
788 d.GetBinDataIntegral(dv);
789
790}
791
792void FillData ( BinData & dv, const TGraph * gr, TF1 * func ) {
793 // fill the data vector from a TGraph. Pass also the TF1 function which is
794 // needed in case to exclude points rejected by the function
795 assert(gr != nullptr);
796
797 // get fit option
798 DataOptions & fitOpt = dv.Opt();
799
801 // adjust option according to type
802 fitOpt.fErrors1 = (type == BinData::kNoError);
803 // set this if we want to have error=1 for points with zero errors (by default they are skipped)
804 // fitOpt.fUseEmpty = true;
805
806 // use coordinate or asym errors in case option is set and type is consistent
808 fitOpt.fAsymErrors &= (type == BinData::kAsymError);
809
810
811 // if data are filled already check if there are consistent - otherwise do nothing
812 if (dv.Size() > 0 && dv.NDim() == 1 ) {
813 // check if size is correct otherwise flag an errors
814 if ( dv.GetErrorType() != type ) {
815 Error("FillData","Inconsistent TGraph with previous data set- skip all graph data");
816 return;
817 }
818 }
819
820 DoFillData(dv, gr, type, func);
821
822}
823
824void FillData ( BinData & dv, const TMultiGraph * mg, TF1 * func ) {
825 // fill the data vector from a TMultiGraph. Pass also the TF1 function which is
826 // needed in case to exclude points rejected by the function
827 assert(mg != nullptr);
828
829 TList * grList = mg->GetListOfGraphs();
830 assert(grList != nullptr);
831
832#ifdef DEBUG
833// grList->Print();
834 TIter itr(grList, kIterBackward);
835 TObject *obj;
836 std::cout << "multi-graph list of graps: " << std::endl;
837 while ((obj = itr())) {
838 std::cout << obj->IsA()->GetName() << std::endl;
839 }
840
841#endif
842
843 // get fit option
844 DataOptions & fitOpt = dv.Opt();
845
846 // loop on the graphs to get the data type (use maximum)
847 TIter next(grList);
848
850 TGraph *gr = nullptr;
851 while ((gr = (TGraph*) next())) {
853 if (t > type ) type = t;
854 }
855 // adjust option according to type
856 fitOpt.fErrors1 = (type == BinData::kNoError);
857 // use coordinate or asym errors in case option is set and type is consistent
859 fitOpt.fAsymErrors &= (type == BinData::kAsymError);
860
861
862#ifdef DEBUG
863 std::cout << "Fitting MultiGraph of type " << type << std::endl;
864#endif
865
866 // fill the data now
867 next = grList;
868 while ((gr = (TGraph*) next())) {
869 DoFillData( dv, gr, type, func);
870 }
871
872#ifdef DEBUG
873 std::cout << "TGraphFitInterface::FillData MultiGraph FitData size is " << dv.Size() << std::endl;
874#endif
875
876}
877
878void FillData ( BinData & dv, const TGraph2D * gr, TF1 * func ) {
879 // fill the data vector from a TGraph2D. Pass also the TF1 function which is
880 // needed in case to exclude points rejected by the function
881 // in case of a pure TGraph
882 assert(gr != nullptr);
883
884 // get fit option
885 DataOptions & fitOpt = dv.Opt();
887 // adjust option according to type
888 fitOpt.fErrors1 = (type == BinData::kNoError);
890 fitOpt.fAsymErrors = false; // a TGraph2D with asymmetric errors does not exist
891
892 int nPoints = gr->GetN();
893 double *gx = gr->GetX();
894 double *gy = gr->GetY();
895 double *gz = gr->GetZ();
896
897 // if all errors are zero set option of using errors to 1
898 if ( gr->GetEZ() == nullptr) fitOpt.fErrors1 = true;
899
900 double x[2];
901 double ex[2];
902
903 // look at data range
904 const DataRange & range = dv.Range();
905 bool useRangeX = ( range.Size(0) > 0);
906 bool useRangeY = ( range.Size(1) > 0);
907 double xmin = 0;
908 double xmax = 0;
909 double ymin = 0;
910 double ymax = 0;
911 range.GetRange(xmin,xmax,ymin,ymax);
912
913 dv.Initialize(nPoints,2, type);
914
915 for ( int i = 0; i < nPoints; ++i) {
916
917 x[0] = gx[i];
918 x[1] = gy[i];
919
920 //if (fitOpt.fUseRange && HFitInterface::IsPointOutOfRange(func, x) ) continue;
921 if (useRangeX && ( x[0] < xmin || x[0] > xmax) ) continue;
922 if (useRangeY && ( x[1] < ymin || x[1] > ymax) ) continue;
923
924 // need to evaluate function to know about rejected points
925 // hugly but no other solutions
926 if (func) {
927 TF1::RejectPoint(false);
928 (*func)( x ); // evaluate using stored function parameters
929 if (TF1::RejectedPoint() ) continue;
930 }
931
932 if (type == BinData::kNoError) {
933 dv.Add( x, gz[i] );
934 continue;
935 }
936
937 double errorZ = gr->GetErrorZ(i);
938 if (!HFitInterface::AdjustError(fitOpt,errorZ) ) continue;
939
940 if (type == BinData::kValueError) {
941 dv.Add( x, gz[i], errorZ );
942 }
943 else if (type == BinData::kCoordError) { // case use error in coordinates (x and y)
944 ex[0] = std::max(gr->GetErrorX(i), 0.);
945 ex[1] = std::max(gr->GetErrorY(i), 0.);
946 dv.Add( x, gz[i], ex, errorZ );
947 }
948 else
949 assert(0); // should not go here
950
951#ifdef DEBUG
952 std::cout << "Point " << i << " " << gx[i] << " " << gy[i] << " " << errorZ << std::endl;
953#endif
954
955 }
956
957#ifdef DEBUG
958 std::cout << "THFitInterface::FillData Graph2D FitData size is " << dv.Size() << std::endl;
959#endif
960
961}
962
963
964// confidence intervals
966 if (h1->GetDimension() != 1) {
967 Error("GetConfidenceIntervals","Invalid object used for storing confidence intervals");
968 return false;
969 }
970 // fill fit data sets with points to estimate cl.
971 BinData d;
972 FillData(d,h1,nullptr);
973 gr->Set(d.NPoints() );
974 double * ci = gr->GetEY(); // make CL values error of the graph
975 result.GetConfidenceIntervals(d,ci,cl);
976 // put function value as abscissa of the graph
977 for (unsigned int ipoint = 0; ipoint < d.NPoints(); ++ipoint) {
978 const double * x = d.Coords(ipoint);
979 const ROOT::Math::IParamMultiFunction * func = result.FittedFunction();
980 gr->SetPoint(ipoint, x[0], (*func)(x) );
981 }
982 return true;
983}
984
985
986} // end namespace Fit
987
988} // end namespace ROOT
#define d(i)
Definition RSha256.hxx:102
#define s1(x)
Definition RSha256.hxx:91
unsigned long long ULong64_t
Definition RtypesCore.h:81
const Bool_t kIterBackward
Definition TCollection.h:43
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:185
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition TError.cxx:229
Option_t Option_t option
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
float xmin
float ymin
float xmax
float ymax
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition BinData.h:52
void AddBinUpEdge(const double *xup)
add the bin width data, a pointer to an array with the bin upper edge information.
Definition BinData.cxx:613
ErrorType GetErrorType() const
retrieve the errortype
Definition BinData.h:562
void Add(double x, double y)
add one dim data with only coordinate and values
Definition BinData.cxx:408
void Initialize(unsigned int newPoints, unsigned int dim=1, ErrorType err=kValueError)
Preallocate a data set with given size, dimension and error type.
Definition BinData.h:122
class describing the range in the coordinates it supports multiple range in a coordinate.
Definition DataRange.h:35
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:71
void GetRange(unsigned int irange, unsigned int icoord, double &xmin, double &xmax) const
get the i-th range for given coordinate.
Definition DataRange.h:104
unsigned int Size() const
return number of fit points
Definition FitData.h:293
unsigned int NDim() const
return coordinate data dimension
Definition FitData.h:301
const DataOptions & Opt() const
access to options
Definition FitData.h:309
const DataRange & Range() const
access to range
Definition FitData.h:321
class containing the result of the fit and all the related information (fitted parameter values,...
Definition FitResult.h:47
SparseData class representing the data of a THNSparse histogram The data needs to be converted to a B...
Definition SparseData.h:35
void Add(std::vector< double > &min, std::vector< double > &max, const double content, const double error=1.0)
Adds a new bin specified by the vectors.
Class to manage histogram axis.
Definition TAxis.h:31
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition TAxis.cxx:478
Double_t GetXmax() const
Definition TAxis.h:140
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
Definition TAxis.cxx:518
virtual Int_t FindFixBin(Double_t x) const
Find bin number corresponding to abscissa x.
Definition TAxis.cxx:419
Int_t GetLast() const
Return last bin on the axis i.e.
Definition TAxis.cxx:469
Double_t GetXmin() const
Definition TAxis.h:139
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width.
Definition TAxis.cxx:540
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition TAxis.cxx:528
Int_t GetFirst() const
Return first bin on the axis i.e.
Definition TAxis.cxx:458
1-Dim function class
Definition TF1.h:233
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:3683
virtual void GetRange(Double_t *xmin, Double_t *xmax) const
Return range of a generic N-D function.
Definition TF1.cxx:2281
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set lower and upper limits for parameter ipar.
Definition TF1.cxx:3507
static Bool_t RejectedPoint()
See TF1::RejectPoint above.
Definition TF1.cxx:3692
virtual void SetParameters(const Double_t *params)
Definition TF1.h:670
virtual void SetParameter(Int_t param, Double_t value)
Definition TF1.h:660
virtual Bool_t IsInside(const Double_t *x) const
return kTRUE if the point is inside the function range
Definition TF1.h:624
virtual Int_t GetNdim() const
Definition TF1.h:511
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition TGraph2D.h:41
A TGraphErrors is a TGraph with error bars.
Double_t GetErrorY(Int_t bin) const override
It returns the error along Y at point i.
Double_t * GetEX() const override
Double_t GetErrorX(Int_t bin) const override
It returns the error along X at point i.
Double_t * GetEY() const override
Double_t GetErrorXhigh(Int_t bin) const override
It returns the error along X at point i.
Double_t GetErrorYlow(Int_t bin) const override
It returns the error along Y at point i.
Double_t GetErrorYhigh(Int_t bin) const override
It returns the error along Y at point i.
TClass * IsA() const override
Double_t GetErrorXlow(Int_t bin) const override
It returns the error along X at point i.
A TGraph is an 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:2319
Double_t * GetY() const
Definition TGraph.h:139
virtual Double_t * GetEYlow() const
Definition TGraph.h:145
Int_t GetN() const
Definition TGraph.h:131
Double_t * GetX() const
Definition TGraph.h:138
virtual Double_t * GetEYhigh() const
Definition TGraph.h:144
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:2254
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
TAxis * GetZaxis()
Definition TH1.h:326
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition TH1.cxx:9027
virtual Int_t GetDimension() const
Definition TH1.h:283
TAxis * GetXaxis()
Definition TH1.h:324
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:4938
virtual Int_t GetNcells() const
Definition TH1.h:300
TAxis * GetYaxis()
Definition TH1.h:325
Bool_t IsBinUnderflow(Int_t bin, Int_t axis=0) const
Return true if the bin is underflow.
Definition TH1.cxx:5178
Bool_t IsBinOverflow(Int_t bin, Int_t axis=0) const
Return true if the bin is overflow.
Definition TH1.cxx:5146
virtual Double_t GetBinLowEdge(Int_t bin) const
Return bin lower edge for 1D histogram.
Definition TH1.cxx:9116
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5025
Multidimensional histogram base.
Definition THnBase.h:43
A doubly linked list.
Definition TList.h:38
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition TMultiGraph.h:34
TList * GetListOfGraphs() const
Definition TMultiGraph.h:68
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
Mother of all ROOT objects.
Definition TObject.h:41
virtual TClass * IsA() const
Definition TObject.h:245
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
Double_t ey[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
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:133
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=nullptr)
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
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
DataOptions : simple structure holding the options on how the data are filled.
Definition DataOptions.h:28
bool fErrors1
use all errors equal to 1, i.e. fit without errors (default is false)
Definition DataOptions.h:52
bool fAsymErrors
use asymmetric errors in the value when available, selecting them according to the on sign of residua...
Definition DataOptions.h:55
bool fNormBinVolume
normalize data by a normalized the bin volume (bin volume divided by a reference value)
Definition DataOptions.h:49
bool fIntegral
use integral of bin content instead of bin center (default is false)
Definition DataOptions.h:47
bool fBinVolume
normalize data by the bin volume (it is used in the Poisson likelihood fits)
Definition DataOptions.h:48
bool fCoordErrors
use errors on the x coordinates when available (default is true)
Definition DataOptions.h:54