Logo ROOT   6.08/07
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 
266 {
267  unsigned int n = data.Size();
268  if (n == 0) return;
269 
270  // find xmin and xmax of the data
271  double valxmin;
272  double xmin = *(data.GetPoint(0,valxmin));
273  double xmax = xmin;
274  double valxmax = valxmin;
275 
276  for (unsigned int i = 1; i < n; ++ i) {
277  double val;
278  double x = *(data.GetPoint(i,val) );
279  if (x < xmin) {
280  xmin = x;
281  valxmin = val;
282  }
283  else if (x > xmax) {
284  xmax = x;
285  valxmax = val;
286  }
287  }
288 
289  // avoid negative values of valxmin/valxmax
290  if (valxmin <= 0 && valxmax > 0 ) valxmin = valxmax;
291  else if (valxmax <=0 && valxmin > 0) valxmax = valxmin;
292  else if (valxmin <=0 && valxmax <= 0) { valxmin = 1; valxmax = 1; }
293 
294  double slope = std::log( valxmax/valxmin) / (xmax - xmin);
295  double constant = std::log(valxmin) - slope * xmin;
296  f1->SetParameters(constant, slope);
297 }
298 
299 
300 ////////////////////////////////////////////////////////////////////////////////
301 /// -*-*-*-*Compute Initial values of parameters for a gaussian
302 /// derivaed from function H1InitGaus defined in TH1.cxx
303 /// ===================================================
304 
306 {
307 
308  static const double sqrtpi = 2.506628;
309 
310  // - Compute mean value and RMS of the data
311  unsigned int n = data.Size();
312  if (n == 0) return;
313  double sumx = 0;
314  double sumx2 = 0;
315  double allcha = 0;
316  double valmax = 0;
317  double rangex = data.Coords(n-1)[0] - data.Coords(0)[0];
318  // to avoid binwidth = 0 set arbitrarly to 1
319  double binwidth = 1;
320  if ( rangex > 0) binwidth = rangex;
321  double x0 = 0;
322  for (unsigned int i = 0; i < n; ++ i) {
323  double val;
324  double x = *(data.GetPoint(i,val) );
325  sumx += val*x;
326  sumx2 += val*x*x;
327  allcha += val;
328  if (val > valmax) valmax = val;
329  if (i > 0) {
330  double dx = x - x0;
331  if (dx < binwidth) binwidth = dx;
332  }
333  x0 = x;
334  }
335 
336  if (allcha <= 0) return;
337  double mean = sumx/allcha;
338  double rms = sumx2/allcha - mean*mean;
339 
340 
341  if (rms > 0)
342  rms = std::sqrt(rms);
343  else
344  rms = binwidth*n/4;
345 
346 
347  //if the distribution is really gaussian, the best approximation
348  //is binwidx*allcha/(sqrtpi*rms)
349  //However, in case of non-gaussian tails, this underestimates
350  //the normalisation constant. In this case the maximum value
351  //is a better approximation.
352  //We take the average of both quantities
353 
354 // printf("valmax %f other %f bw %f allcha %f rms %f \n",valmax, binwidth*allcha/(sqrtpi*rms),
355 // binwidth, allcha,rms );
356 
357  double constant = 0.5*(valmax+ binwidth*allcha/(sqrtpi*rms));
358 
359 
360  //In case the mean value is outside the histo limits and
361  //the RMS is bigger than the range, we take
362  // mean = center of bins
363  // rms = half range
364 // Double_t xmin = curHist->GetXaxis()->GetXmin();
365 // Double_t xmax = curHist->GetXaxis()->GetXmax();
366 // if ((mean < xmin || mean > xmax) && rms > (xmax-xmin)) {
367 // mean = 0.5*(xmax+xmin);
368 // rms = 0.5*(xmax-xmin);
369 // }
370 
371  f1->SetParameter(0,constant);
372  f1->SetParameter(1,mean);
373  f1->SetParameter(2,rms);
374  f1->SetParLimits(2,0,10*rms);
375 
376 
377 #ifdef DEBUG
378  std::cout << "Gaussian initial par values" << constant << " " << mean << " " << rms << std::endl;
379 #endif
380 
381 }
382 
383 ////////////////////////////////////////////////////////////////////////////////
384 /// -*-*-*-*Compute Initial values of parameters for a gaussian
385 /// derivaed from function H1InitGaus defined in TH1.cxx
386 /// ===================================================
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 == 0 && ( eyl == 0 || eyh == 0 ) ) ) {
488  type = BinData::kNoError;
489  }
490  // need to treat case when all errors are zero
491  // note that by default fitOpt.fCoordError is true
492  else if ( ex != 0 && 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 != 0 && eyh != 0) && 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)
514  type = BinData::kNoError;
515  else if (!zeroErrorX && zeroErrorY)
516  type = BinData::kCoordError;
517  else if (zeroErrorX && !zeroErrorY) {
518  type = BinData::kAsymError;
519  fitOpt.fCoordErrors = false;
520  }
521  else {
522  type = BinData::kAsymError;
523  }
524  }
525 
526  // need to look also a case when all errors in y are zero
527  if ( ey != 0 && 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 == 0 ) {
555  type = BinData::kNoError;
556  }
557  else if ( ex != 0 && ey!=0 && 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 
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 
673 void 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 
717 void 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 oveflows! (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 
760 void 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 
792 void 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 != 0);
796 
797  // get fit option
798  DataOptions & fitOpt = dv.Opt();
799 
800  BinData::ErrorType type = GetDataType(gr,fitOpt);
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
807  fitOpt.fCoordErrors &= (type == BinData::kCoordError) || (type == BinData::kAsymError) ;
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.PointSize() == 2 && type != BinData::kNoError ) {
815  Error("FillData","Inconsistent TGraph with previous data set- skip all graph data");
816  return;
817  }
818  if (dv.PointSize() == 3 && type != BinData::kValueError ) {
819  Error("FillData","Inconsistent TGraph with previous data set- skip all graph data");
820  return;
821  }
822  if (dv.PointSize() == 4 && type != BinData::kCoordError ) {
823  Error("FillData","Inconsistent TGraph with previous data set- skip all graph data");
824  return;
825  }
826  if (dv.PointSize() == 5 && type != BinData::kAsymError ) {
827  Error("FillData","Inconsistent TGraph with previous data set- skip all graph data");
828  return;
829  }
830  }
831 
832  DoFillData(dv, gr, type, func);
833 
834 }
835 
836 void FillData ( BinData & dv, const TMultiGraph * mg, TF1 * func ) {
837  // fill the data vector from a TMultiGraph. Pass also the TF1 function which is
838  // needed in case to exclude points rejected by the function
839  assert(mg != 0);
840 
841  TList * grList = mg->GetListOfGraphs();
842  assert(grList != 0);
843 
844 #ifdef DEBUG
845 // grList->Print();
846  TIter itr(grList, kIterBackward);
847  TObject *obj;
848  std::cout << "multi-graph list of graps: " << std::endl;
849  while ((obj = itr())) {
850  std::cout << obj->IsA()->GetName() << std::endl;
851  }
852 
853 #endif
854 
855  // get fit option
856  DataOptions & fitOpt = dv.Opt();
857 
858  // loop on the graphs to get the data type (use maximum)
859  TIter next(grList);
860 
862  TGraph *gr = 0;
863  while ((gr = (TGraph*) next())) {
864  BinData::ErrorType t = GetDataType(gr,fitOpt);
865  if (t > type ) type = t;
866  }
867  // adjust option according to type
868  fitOpt.fErrors1 = (type == BinData::kNoError);
869  fitOpt.fCoordErrors = (type == BinData::kCoordError);
870  fitOpt.fAsymErrors = (type == BinData::kAsymError);
871 
872 
873 #ifdef DEBUG
874  std::cout << "Fitting MultiGraph of type " << type << std::endl;
875 #endif
876 
877  // fill the data now
878  next = grList;
879  while ((gr = (TGraph*) next())) {
880  DoFillData( dv, gr, type, func);
881  }
882 
883 #ifdef DEBUG
884  std::cout << "TGraphFitInterface::FillData MultiGraph FitData size is " << dv.Size() << std::endl;
885 #endif
886 
887 }
888 
889 void FillData ( BinData & dv, const TGraph2D * gr, TF1 * func ) {
890  // fill the data vector from a TGraph2D. Pass also the TF1 function which is
891  // needed in case to exclude points rejected by the function
892  // in case of a pure TGraph
893  assert(gr != 0);
894 
895  // get fit option
896  DataOptions & fitOpt = dv.Opt();
897  BinData::ErrorType type = GetDataType(gr,fitOpt);
898  // adjust option according to type
899  fitOpt.fErrors1 = (type == BinData::kNoError);
900  fitOpt.fCoordErrors = (type == BinData::kCoordError);
901  fitOpt.fAsymErrors = false; // a TGraph2D with asymmetric errors does not exist
902 
903  int nPoints = gr->GetN();
904  double *gx = gr->GetX();
905  double *gy = gr->GetY();
906  double *gz = gr->GetZ();
907 
908  // if all errors are zero set option of using errors to 1
909  if ( gr->GetEZ() == 0) fitOpt.fErrors1 = true;
910 
911  double x[2];
912  double ex[2];
913 
914  // look at data range
915  const DataRange & range = dv.Range();
916  bool useRangeX = ( range.Size(0) > 0);
917  bool useRangeY = ( range.Size(1) > 0);
918  double xmin = 0;
919  double xmax = 0;
920  double ymin = 0;
921  double ymax = 0;
922  range.GetRange(xmin,xmax,ymin,ymax);
923 
924  dv.Initialize(nPoints,2, type);
925 
926  for ( int i = 0; i < nPoints; ++i) {
927 
928  x[0] = gx[i];
929  x[1] = gy[i];
930 
931  //if (fitOpt.fUseRange && HFitInterface::IsPointOutOfRange(func, x) ) continue;
932  if (useRangeX && ( x[0] < xmin || x[0] > xmax) ) continue;
933  if (useRangeY && ( x[1] < ymin || x[1] > ymax) ) continue;
934 
935  // need to evaluate function to know about rejected points
936  // hugly but no other solutions
937  if (func) {
938  TF1::RejectPoint(false);
939  (*func)( x ); // evaluate using stored function parameters
940  if (TF1::RejectedPoint() ) continue;
941  }
942 
943  if (type == BinData::kNoError) {
944  dv.Add( x, gz[i] );
945  continue;
946  }
947 
948  double errorZ = gr->GetErrorZ(i);
949  if (!HFitInterface::AdjustError(fitOpt,errorZ) ) continue;
950 
951  if (type == BinData::kValueError) {
952  dv.Add( x, gz[i], errorZ );
953  }
954  else if (type == BinData::kCoordError) { // case use error in coordinates (x and y)
955  ex[0] = std::max(gr->GetErrorX(i), 0.);
956  ex[1] = std::max(gr->GetErrorY(i), 0.);
957  dv.Add( x, gz[i], ex, errorZ );
958  }
959  else
960  assert(0); // should not go here
961 
962 #ifdef DEBUG
963  std::cout << "Point " << i << " " << gx[i] << " " << gy[i] << " " << errorZ << std::endl;
964 #endif
965 
966  }
967 
968 #ifdef DEBUG
969  std::cout << "THFitInterface::FillData Graph2D FitData size is " << dv.Size() << std::endl;
970 #endif
971 
972 }
973 
974 
975 // confidence intervals
976 bool GetConfidenceIntervals(const TH1 * h1, const ROOT::Fit::FitResult & result, TGraphErrors * gr, double cl ) {
977  if (h1->GetDimension() != 1) {
978  Error("GetConfidenceIntervals","Invalid object used for storing confidence intervals");
979  return false;
980  }
981  // fill fit data sets with points to estimate cl.
982  BinData d;
983  FillData(d,h1,0);
984  gr->Set(d.NPoints() );
985  double * ci = gr->GetEY(); // make CL values error of the graph
986  result.GetConfidenceIntervals(d,ci,cl);
987  // put function value as abscissa of the graph
988  for (unsigned int ipoint = 0; ipoint < d.NPoints(); ++ipoint) {
989  const double * x = d.Coords(ipoint);
991  gr->SetPoint(ipoint, x[0], (*func)(x) );
992  }
993  return true;
994 }
995 
996 
997 } // end namespace Fit
998 
999 } // end namespace ROOT
1000 
unsigned int Size() const
return number of fit points
Definition: BinData.h:447
Double_t GetBinError(const Int_t *idx) const
Definition: THnBase.h:169
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
virtual Int_t GetNcells() const
Definition: TH1.h:304
virtual Double_t GetErrorYhigh(Int_t bin) const
This function is called by GraphFitChisquare.
Definition: TGraph.cxx:1433
void Initialize(unsigned int maxpoints, unsigned int dim=1, ErrorType err=kValueError)
preallocate a data set with given size , dimension and error type (to get the full point size) If the...
Definition: BinData.cxx:216
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:439
float xmin
Definition: THbookFile.cxx:93
virtual Double_t * GetEX() const
Definition: TGraph.h:142
Bool_t IsBinUnderflow(Int_t bin, Int_t axis=0) const
Return true if the bin is underflow.
Definition: TH1.cxx:4772
Int_t GetFirst() const
Return first bin on the axis i.e.
Definition: TAxis.cxx:444
virtual Double_t GetErrorZ(Int_t bin) const
This function is called by Graph2DFitChisquare.
Definition: TGraph2D.cxx:1019
virtual Double_t GetErrorY(Int_t bin) const
This function is called by GraphFitChisquare.
Definition: TGraph.cxx:1403
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
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:4740
virtual Long64_t GetNbins() const =0
void GetConfidenceIntervals(unsigned int n, unsigned int stride1, unsigned int stride2, const double *x, double *ci, double cl=0.95, bool norm=true) const
get confidence intervals for an array of n points x.
Definition: FitResult.cxx:560
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
Definition: TAxis.cxx:504
float ymin
Definition: THbookFile.cxx:93
virtual Double_t * GetEYlow() const
Definition: TGraph.h:147
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4638
Double_t GetBinContent(const Int_t *idx) const
Definition: THnBase.h:177
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition: TMultiGraph.h:37
virtual Double_t * GetEY() const
Definition: TGraph.h:143
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:4553
unsigned int NPoints() const
return number of fit points
Definition: BinData.h:442
void ExamineRange(const TAxis *axis, std::pair< double, double > range, int &hxfirst, int &hxlast)
const IModelFunction * FittedFunction() const
fitting quantities
Definition: FitResult.h:153
void AddBinUpEdge(const double *xup)
add the bin width data, a pointer to an array with the bin upper edge information.
Definition: BinData.cxx:452
virtual Double_t GetErrorYlow(Int_t bin) const
This function is called by GraphFitChisquare.
Definition: TGraph.cxx:1443
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:514
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:304
virtual Int_t GetDimension() const
Definition: TH1.h:287
double sqrt(double)
Double_t GetXmin() const
Definition: TAxis.h:139
Double_t x[n]
Definition: legend1.C:17
virtual Double_t * GetEY() const
Definition: TGraph2D.h:132
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:464
bool AdjustError(const DataOptions &option, double &error, double value=1)
virtual Int_t GetNdim() const
Definition: TF1.h:350
TH1F * h1
Definition: legend1.C:5
virtual Double_t * GetEZ() const
Definition: TGraph2D.h:133
A doubly linked list.
Definition: TList.h:47
virtual Double_t GetErrorY(Int_t bin) const
This function is called by Graph2DFitChisquare.
Definition: TGraph2D.cxx:1009
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:3223
virtual Double_t GetErrorXlow(Int_t bin) const
This function is called by GraphFitChisquare.
Definition: TGraph.cxx:1423
void FillData(BinData &dv, const TH1 *hist, TF1 *func=0)
fill the data vector from a TH1.
TAxis * GetAxis(Int_t dim) const
Definition: THnBase.h:135
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:455
Class to manage histogram axis.
Definition: TAxis.h:36
IParamFunction interface (abstract class) describing multi-dimensional parameteric functions It is a ...
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
void GetBinDataIntegral(BinData &) const
Definition: SparseData.cxx:320
const double * Coords(unsigned int ipoint) const
return a pointer to the coordinates data for the given fit point
Definition: BinData.h:226
const DataOptions & Opt() const
access to options
Definition: DataVector.h:97
Int_t GetN() const
Definition: TGraph.h:133
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:3407
TAxis * GetYaxis()
Definition: TH1.h:325
Double_t * GetEY() const
Definition: TGraphErrors.h:69
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:61
float xmax
Definition: THbookFile.cxx:93
unsigned int PointSize() const
return the size of a fit point (is the coordinate dimension + 1 for the value and eventually the numb...
Definition: BinData.h:150
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
virtual Double_t * GetEX() const
Definition: TGraph2D.h:131
Double_t * GetY() const
Definition: TGraph2D.h:129
Double_t * GetX() const
Definition: TGraph.h:140
class containg the result of the fit and all the related information (fitted parameter values...
Definition: FitResult.h:52
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:265
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:70
Double_t y[n]
Definition: legend1.C:17
double func(double *x, double *p)
Definition: stressTF1.cxx:213
Double_t ey[n]
Definition: legend1.C:17
void GetRange(unsigned int icoord, double &xmin, double &xmax) const
get the first range for given coordinate.
Definition: DataRange.h:103
Double_t * GetZ() const
Definition: TGraph2D.h:130
The TH1 histogram class.
Definition: TH1.h:80
virtual Double_t GetErrorX(Int_t bin) const
This function is called by Graph2DFitChisquare.
Definition: TGraph2D.cxx:999
Int_t GetNdimensions() const
Definition: THnBase.h:145
TAxis * GetZaxis()
Definition: TH1.h:326
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:405
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
unsigned int NDim() const
return coordinate data dimension
Definition: BinData.h:452
virtual Bool_t IsInside(const Double_t *x) const
return kTRUE if the point is inside the function range
Definition: TF1.h:411
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width.
Definition: TAxis.cxx:526
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2150
TList * GetListOfGraphs() const
Definition: TMultiGraph.h:71
Double_t * GetY() const
Definition: TGraph.h:141
static Bool_t RejectedPoint()
See TF1::RejectPoint above.
Definition: TF1.cxx:3416
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:1989
1-Dim function class
Definition: TF1.h:149
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:53
Int_t GetN() const
Definition: TGraph2D.h:127
TF1 * f1
Definition: legend1.C:11
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:28
Multidimensional histogram base.
Definition: THnBase.h:53
double result[121]
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:431
const DataRange & Range() const
access to range
Definition: DataVector.h:103
virtual Double_t * GetEYhigh() const
Definition: TGraph.h:146
const Bool_t kIterBackward
Definition: TCollection.h:44
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition: TGraph2D.h:50
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:416
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:2099
Double_t GetXmax() const
Definition: TAxis.h:140
Double_t ex[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
virtual Double_t GetErrorXhigh(Int_t bin) const
This function is called by GraphFitChisquare.
Definition: TGraph.cxx:1413
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()
Definition: TH1.h:324
void Error(ErrorHandler_t func, int code, const char *va_(fmt),...)
Write error message and call a handler, if required.
Double_t * GetX() const
Definition: TGraph2D.h:128
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8173