Logo ROOT   master
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 
46 namespace ROOT {
47 
48 namespace Fit {
49 
50 // add a namespace to distinguish from the Graph functions
51 namespace HFitInterface {
52 
53 
54 bool 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 
60 bool 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 
79 void 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 
107 void 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 
264 void InitExpo(const ROOT::Fit::BinData & data, TF1 * f1)
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 
303 void InitGaus(const ROOT::Fit::BinData & data, TF1 * f1)
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 
385 void Init2DGaus(const ROOT::Fit::BinData & data, TF1 * f1)
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 
573 void 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 
670 void 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 
714 void 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 
757 void 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 
789 void 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 
821 void 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())) {
849  BinData::ErrorType t = GetDataType(gr,fitOpt);
850  if (t > type ) type = t;
851  }
852  // adjust option according to type
853  fitOpt.fErrors1 = (type == BinData::kNoError);
854  // use coordinate or asym errors in case option is set and type is consistent
856  fitOpt.fAsymErrors &= (type == BinData::kAsymError);
857 
858 
859 #ifdef DEBUG
860  std::cout << "Fitting MultiGraph of type " << type << std::endl;
861 #endif
862 
863  // fill the data now
864  next = grList;
865  while ((gr = (TGraph*) next())) {
866  DoFillData( dv, gr, type, func);
867  }
868 
869 #ifdef DEBUG
870  std::cout << "TGraphFitInterface::FillData MultiGraph FitData size is " << dv.Size() << std::endl;
871 #endif
872 
873 }
874 
875 void FillData ( BinData & dv, const TGraph2D * gr, TF1 * func ) {
876  // fill the data vector from a TGraph2D. Pass also the TF1 function which is
877  // needed in case to exclude points rejected by the function
878  // in case of a pure TGraph
879  assert(gr != 0);
880 
881  // get fit option
882  DataOptions & fitOpt = dv.Opt();
884  // adjust option according to type
885  fitOpt.fErrors1 = (type == BinData::kNoError);
886  fitOpt.fCoordErrors = (type == BinData::kCoordError);
887  fitOpt.fAsymErrors = false; // a TGraph2D with asymmetric errors does not exist
888 
889  int nPoints = gr->GetN();
890  double *gx = gr->GetX();
891  double *gy = gr->GetY();
892  double *gz = gr->GetZ();
893 
894  // if all errors are zero set option of using errors to 1
895  if ( gr->GetEZ() == 0) fitOpt.fErrors1 = true;
896 
897  double x[2];
898  double ex[2];
899 
900  // look at data range
901  const DataRange & range = dv.Range();
902  bool useRangeX = ( range.Size(0) > 0);
903  bool useRangeY = ( range.Size(1) > 0);
904  double xmin = 0;
905  double xmax = 0;
906  double ymin = 0;
907  double ymax = 0;
908  range.GetRange(xmin,xmax,ymin,ymax);
909 
910  dv.Initialize(nPoints,2, type);
911 
912  for ( int i = 0; i < nPoints; ++i) {
913 
914  x[0] = gx[i];
915  x[1] = gy[i];
916 
917  //if (fitOpt.fUseRange && HFitInterface::IsPointOutOfRange(func, x) ) continue;
918  if (useRangeX && ( x[0] < xmin || x[0] > xmax) ) continue;
919  if (useRangeY && ( x[1] < ymin || x[1] > ymax) ) continue;
920 
921  // need to evaluate function to know about rejected points
922  // hugly but no other solutions
923  if (func) {
924  TF1::RejectPoint(false);
925  (*func)( x ); // evaluate using stored function parameters
926  if (TF1::RejectedPoint() ) continue;
927  }
928 
929  if (type == BinData::kNoError) {
930  dv.Add( x, gz[i] );
931  continue;
932  }
933 
934  double errorZ = gr->GetErrorZ(i);
935  if (!HFitInterface::AdjustError(fitOpt,errorZ) ) continue;
936 
937  if (type == BinData::kValueError) {
938  dv.Add( x, gz[i], errorZ );
939  }
940  else if (type == BinData::kCoordError) { // case use error in coordinates (x and y)
941  ex[0] = std::max(gr->GetErrorX(i), 0.);
942  ex[1] = std::max(gr->GetErrorY(i), 0.);
943  dv.Add( x, gz[i], ex, errorZ );
944  }
945  else
946  assert(0); // should not go here
947 
948 #ifdef DEBUG
949  std::cout << "Point " << i << " " << gx[i] << " " << gy[i] << " " << errorZ << std::endl;
950 #endif
951 
952  }
953 
954 #ifdef DEBUG
955  std::cout << "THFitInterface::FillData Graph2D FitData size is " << dv.Size() << std::endl;
956 #endif
957 
958 }
959 
960 
961 // confidence intervals
962 bool GetConfidenceIntervals(const TH1 * h1, const ROOT::Fit::FitResult & result, TGraphErrors * gr, double cl ) {
963  if (h1->GetDimension() != 1) {
964  Error("GetConfidenceIntervals","Invalid object used for storing confidence intervals");
965  return false;
966  }
967  // fill fit data sets with points to estimate cl.
968  BinData d;
969  FillData(d,h1,0);
970  gr->Set(d.NPoints() );
971  double * ci = gr->GetEY(); // make CL values error of the graph
972  result.GetConfidenceIntervals(d,ci,cl);
973  // put function value as abscissa of the graph
974  for (unsigned int ipoint = 0; ipoint < d.NPoints(); ++ipoint) {
975  const double * x = d.Coords(ipoint);
976  const ROOT::Math::IParamMultiFunction * func = result.FittedFunction();
977  gr->SetPoint(ipoint, x[0], (*func)(x) );
978  }
979  return true;
980 }
981 
982 
983 } // end namespace Fit
984 
985 } // end namespace ROOT
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual Int_t GetNcells() const
Definition: TH1.h:295
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:638
float xmin
Definition: THbookFile.cxx:93
Double_t * GetEX() const
Definition: TGraphErrors.h:67
Bool_t IsBinUnderflow(Int_t bin, Int_t axis=0) const
Return true if the bin is underflow.
Definition: TH1.cxx:5060
Double_t GetErrorYlow(Int_t bin) const
This function is called by GraphFitChisquare.
Int_t GetFirst() const
Return first bin on the axis i.e.
Definition: TAxis.cxx:455
Returns the available number of logical cores.
Definition: StringConv.hxx:21
Bool_t IsBinOverflow(Int_t bin, Int_t axis=0) const
Return true if the bin is overflow.
Definition: TH1.cxx:5028
ErrorType GetErrorType() const
retrieve the errortype
Definition: BinData.h:551
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
Definition: TAxis.cxx:515
float ymin
Definition: THbookFile.cxx:93
virtual Double_t * GetEYlow() const
Definition: TGraph.h:137
Double_t GetErrorYhigh(Int_t bin) const
This function is called by GraphFitChisquare.
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:545
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4907
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition: TMultiGraph.h:36
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:4822
void ExamineRange(const TAxis *axis, std::pair< double, double > range, int &hxfirst, int &hxlast)
const IModelFunction * FittedFunction() const
fitting quantities
Definition: FitResult.h:148
virtual Double_t GetBinLowEdge(Int_t bin) const
Return bin lower edge for 1D histogram.
Definition: TH1.cxx:8608
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
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 ...
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition: TAxis.cxx:525
const double * Coords(unsigned int ipoint) const
return a pointer to the coordinates data for the given fit point
Definition: FitData.h:246
const double * GetPoint(unsigned int ipoint, double &value) const
retrieve at the same time a pointer to the coordinate data and the fit value More efficient than call...
Definition: BinData.h:370
virtual Int_t GetDimension() const
Definition: TH1.h:278
double sqrt(double)
Double_t GetErrorXlow(Int_t bin) const
This function is called by GraphFitChisquare.
Double_t GetXmin() const
Definition: TAxis.h:133
Double_t x[n]
Definition: legend1.C:17
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
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:475
void Initialize(unsigned int newPoints, unsigned int dim=1, ErrorType err=kValueError)
Definition: BinData.cxx:349
bool AdjustError(const DataOptions &option, double &error, double value=1)
virtual Int_t GetNdim() const
Definition: TF1.h:479
static constexpr double s
TH1F * h1
Definition: legend1.C:5
unsigned int Size() const
return number of fit points
Definition: FitData.h:303
void Error(const char *location, const char *msgfmt,...)
Double_t GetErrorX(Int_t bin) const
This function is called by GraphFitChisquare.
A doubly linked list.
Definition: TList.h:44
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 ...
float ymax
Definition: THbookFile.cxx:93
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set limits for parameter ipar.
Definition: TF1.cxx:3517
void FillData(BinData &dv, const TH1 *hist, TF1 *func=0)
fill the data vector from a TH1.
BinData::ErrorType GetDataType(const TGraph *gr, DataOptions &fitOpt)
DataOptions : simple structure holding the options on how the data are filled.
Definition: DataOptions.h:28
Int_t GetLast() const
Return last bin on the axis i.e.
Definition: TAxis.cxx:466
Class to manage histogram axis.
Definition: TAxis.h:30
void DoFillData(BinData &dv, const TGraph *gr, BinData::ErrorType type, TF1 *func)
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
static constexpr double mg
const DataOptions & Opt() const
access to options
Definition: FitData.h:319
Int_t GetN() const
Definition: TGraph.h:123
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:3670
#define s1(x)
Definition: RSha256.hxx:91
TAxis * GetYaxis()
Definition: TH1.h:317
Double_t * GetEY() const
Definition: TGraphErrors.h:68
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:53
float xmax
Definition: THbookFile.cxx:93
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 Warning(const char *location, const char *msgfmt,...)
TGraphErrors * gr
Definition: legend1.C:25
Double_t * GetX() const
Definition: TGraph.h:130
#define d(i)
Definition: RSha256.hxx:102
class containg the result of the fit and all the related information (fitted parameter values...
Definition: FitResult.h:47
Double_t GetErrorY(Int_t bin) const
This function is called by GraphFitChisquare.
class describing the range in the coordinates it supports multiple range in a coordinate.
Definition: DataRange.h:34
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 Add(double x, double y)
add one dim data with only coordinate and values
Definition: BinData.cxx:422
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
int type
Definition: TGX11.cxx:120
unsigned long long ULong64_t
Definition: RtypesCore.h:72
Double_t y[n]
Definition: legend1.C:17
Double_t ey[n]
Definition: legend1.C:17
The TH1 histogram class.
Definition: TH1.h:56
TAxis * GetZaxis()
Definition: TH1.h:318
Mother of all ROOT objects.
Definition: TObject.h:37
virtual Int_t FindFixBin(Double_t x) const
Find bin number corresponding to abscissa x.
Definition: TAxis.cxx:416
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
virtual Bool_t IsInside(const Double_t *x) const
return kTRUE if the point is inside the function range
Definition: TF1.h:592
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width.
Definition: TAxis.cxx:537
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2261
Double_t * GetY() const
Definition: TGraph.h:131
static Bool_t RejectedPoint()
See TF1::RejectPoint above.
Definition: TF1.cxx:3679
bool IsPointOutOfRange(const TF1 *func, const double *x)
virtual void GetRange(Double_t *xmin, Double_t *xmax) const
Return range of a generic N-D function.
Definition: TF1.cxx:2282
1-Dim function class
Definition: TF1.h:210
A TGraph is an object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
TF1 * f1
Definition: legend1.C:11
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:26
Multidimensional histogram base.
Definition: THnBase.h:43
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:628
unsigned int NDim() const
return coordinate data dimension
Definition: FitData.h:311
const DataRange & Range() const
access to range
Definition: FitData.h:331
virtual Double_t * GetEYhigh() const
Definition: TGraph.h:136
const Bool_t kIterBackward
Definition: TCollection.h:41
Double_t GetErrorXhigh(Int_t bin) const
This function is called by GraphFitChisquare.
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition: TGraph2D.h:41
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:357
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:2196
Double_t GetXmax() const
Definition: TAxis.h:134
Double_t ex[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
void Add(std::vector< double > &min, std::vector< double > &max, const double content, const double error=1.0)
Definition: SparseData.cxx:231
double log(double)
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8519