diff --git a/hist/hist/src/HFitInterface.cxx b/hist/hist/src/HFitInterface.cxx index dc287f91df8ff..eb5b35ee43267 100644 --- a/hist/hist/src/HFitInterface.cxx +++ b/hist/hist/src/HFitInterface.cxx @@ -45,941 +45,940 @@ namespace ROOT { -namespace Fit { - -// add a namespace to distinguish from the Graph functions -namespace HFitInterface { - - -bool IsPointOutOfRange(const TF1 * func, const double * x) { - // function to check if a point is outside range - if (func ==0) return false; - return !func->IsInside(x); -} - -bool AdjustError(const DataOptions & option, double & error, double value = 1) { - // adjust the given error according to the option - // return false when point must be skipped. - // When point error = 0, the point is kept if the option UseEmpty is set or if - // fErrors1 is set and the point value is not zero. - // The value should be used only for points representing counts (histograms), not for the graph. - // In the graph points with zero errors are by default skipped indepentently of the value. - // If one wants to keep the points, the option fUseEmpty must be set - - if (error <= 0) { - if (option.fUseEmpty || (option.fErrors1 && std::abs(value) > 0 ) ) - error = 1.; // set error to 1 - else - return false; // skip bins with zero errors or empty - } else if (option.fErrors1) - error = 1; // set all error to 1 for non-empty bins - return true; -} - -void ExamineRange(const TAxis * axis, std::pair range,int &hxfirst,int &hxlast) { - // examine the range given with the pair on the given histogram axis - // correct in case the bin values hxfirst hxlast - double xlow = range.first; - double xhigh = range.second; -#ifdef DEBUG - std::cout << "xlow " << xlow << " xhigh = " << xhigh << std::endl; -#endif - // ignore ranges specified outside histogram range - int ilow = axis->FindFixBin(xlow); - int ihigh = axis->FindFixBin(xhigh); - if (ilow > hxlast || ihigh < hxfirst) { - Warning("ROOT::Fit::FillData","fit range is outside histogram range, no fit data for %s",axis->GetName()); - } - // consider only range defined with-in histogram not oustide. Always exclude underflow/overflow - hxfirst = std::min( std::max( ilow, hxfirst), hxlast+1) ; - hxlast = std::max( std::min( ihigh, hxlast), hxfirst-1) ; - // exclude bins where range coverage is less than half bin width - if (hxfirst < hxlast) { - if ( axis->GetBinCenter(hxfirst) < xlow) hxfirst++; - if ( axis->GetBinCenter(hxlast) > xhigh) hxlast--; - } -} + namespace Fit { + // add a namespace to distinguish from the Graph functions + namespace HFitInterface { -} // end namespace HFitInterface + bool IsPointOutOfRange(const TF1 * func, const double * x) { + // function to check if a point is outside range + if (func ==0) return false; + return !func->IsInside(x); + } -void FillData(BinData & dv, const TH1 * hfit, TF1 * func) -{ - // Function to fill the binned Fit data structure from a TH1 - // The dimension of the data is the same of the histogram dimension - // The function pointer is need in case of integral is used and to reject points - // rejected in the function + bool AdjustError(const DataOptions & option, double & error, double value = 1) { + // adjust the given error according to the option + // return false when point must be skipped. + // When point error = 0, the point is kept if the option UseEmpty is set or if + // fErrors1 is set and the point value is not zero. + // The value should be used only for points representing counts (histograms), not for the graph. + // In the graph points with zero errors are by default skipped indepentently of the value. + // If one wants to keep the points, the option fUseEmpty must be set + + if (error <= 0) { + if (option.fUseEmpty || (option.fErrors1 && std::abs(value) > 0 ) ) + error = 1.; // set error to 1 + else + return false; // skip bins with zero errors or empty + } else if (option.fErrors1) + error = 1; // set all error to 1 for non-empty bins + return true; + } - // the TF1 pointer cannot be constant since EvalPar and InitArgs are not const methods + void ExamineRange(const TAxis * axis, std::pair range,int &hxfirst,int &hxlast) { + // examine the range given with the pair on the given histogram axis + // correct in case the bin values hxfirst hxlast + double xlow = range.first; + double xhigh = range.second; +#ifdef DEBUG + std::cout << "xlow " << xlow << " xhigh = " << xhigh << std::endl; +#endif + // ignore ranges specified outside histogram range + int ilow = axis->FindFixBin(xlow); + int ihigh = axis->FindFixBin(xhigh); + if (ilow > hxlast || ihigh < hxfirst) { + Warning("ROOT::Fit::FillData","fit range is outside histogram range, no fit data for %s",axis->GetName()); + } + // consider only range defined with-in histogram not oustide. Always exclude underflow/overflow + hxfirst = std::min( std::max( ilow, hxfirst), hxlast+1) ; + hxlast = std::max( std::min( ihigh, hxlast), hxfirst-1) ; + // exclude bins where range coverage is less than half bin width + if (hxfirst < hxlast) { + if ( axis->GetBinCenter(hxfirst) < xlow) hxfirst++; + if ( axis->GetBinCenter(hxlast) > xhigh) hxlast--; + } + } - // get fit option - const DataOptions & fitOpt = dv.Opt(); + } // end namespace HFitInterface - // store instead of bin center the bin edges - bool useBinEdges = fitOpt.fIntegral || fitOpt.fBinVolume; - assert(hfit != 0); + void FillData(BinData & dv, const TH1 * hfit, TF1 * func) + { + // Function to fill the binned Fit data structure from a TH1 + // The dimension of the data is the same of the histogram dimension + // The function pointer is need in case of integral is used and to reject points + // rejected in the function - //std::cout << "creating Fit Data from histogram " << hfit->GetName() << std::endl; + // the TF1 pointer cannot be constant since EvalPar and InitArgs are not const methods - int hxfirst = hfit->GetXaxis()->GetFirst(); - int hxlast = hfit->GetXaxis()->GetLast(); + // get fit option + const DataOptions & fitOpt = dv.Opt(); - int hyfirst = hfit->GetYaxis()->GetFirst(); - int hylast = hfit->GetYaxis()->GetLast(); - int hzfirst = hfit->GetZaxis()->GetFirst(); - int hzlast = hfit->GetZaxis()->GetLast(); + // store instead of bin center the bin edges + bool useBinEdges = fitOpt.fIntegral || fitOpt.fBinVolume; - // function by default has same range (use that one if requested otherwise use data one) + assert(hfit != 0); + //std::cout << "creating Fit Data from histogram " << hfit->GetName() << std::endl; - // get the range (add the function range ??) - // to check if inclusion/exclusion at end/point - const DataRange & range = dv.Range(); - if (range.Size(0) != 0) { - HFitInterface::ExamineRange( hfit->GetXaxis(), range(0), hxfirst, hxlast); - if (range.Size(0) > 1 ) { - Warning("ROOT::Fit::FillData","support only one range interval for X coordinate"); - } - } + int hxfirst = hfit->GetXaxis()->GetFirst(); + int hxlast = hfit->GetXaxis()->GetLast(); - if (hfit->GetDimension() > 1 && range.Size(1) != 0) { - HFitInterface::ExamineRange( hfit->GetYaxis(), range(1), hyfirst, hylast); - if (range.Size(1) > 1 ) - Warning("ROOT::Fit::FillData","support only one range interval for Y coordinate"); - } + int hyfirst = hfit->GetYaxis()->GetFirst(); + int hylast = hfit->GetYaxis()->GetLast(); - if (hfit->GetDimension() > 2 && range.Size(2) != 0) { - HFitInterface::ExamineRange( hfit->GetZaxis(), range(2), hzfirst, hzlast); - if (range.Size(2) > 1 ) - Warning("ROOT::Fit::FillData","support only one range interval for Z coordinate"); - } + int hzfirst = hfit->GetZaxis()->GetFirst(); + int hzlast = hfit->GetZaxis()->GetLast(); + // function by default has same range (use that one if requested otherwise use data one) - int n = (hxlast-hxfirst+1)*(hylast-hyfirst+1)*(hzlast-hzfirst+1); + const DataRange & range = dv.Range(); #ifdef DEBUG - std::cout << "THFitInterface: ifirst = " << hxfirst << " ilast = " << hxlast - << " total bins " << n - << std::endl; + std::cout << "THFitInterface: ifirst = " << hxfirst << " ilast = " << hxlast + << " total bins " << n + << std::endl; #endif - // reserve n for more efficient usage - //dv.Data().reserve(n); - - int hdim = hfit->GetDimension(); - int ndim = hdim; - // case of function dimension less than histogram - if (func !=0 && func->GetNdim() == hdim-1) ndim = hdim-1; - - assert( ndim > 0 ); - //typedef BinPoint::CoordData CoordData; - //CoordData x = CoordData( hfit->GetDimension() ); - dv.Initialize(n,ndim); - - double x[3]; - double s[3]; - - int binx = 0; - int biny = 0; - int binz = 0; - - const TAxis *xaxis = hfit->GetXaxis(); - const TAxis *yaxis = hfit->GetYaxis(); - const TAxis *zaxis = hfit->GetZaxis(); - - for ( binx = hxfirst; binx <= hxlast; ++binx) { - if (useBinEdges) { - x[0] = xaxis->GetBinLowEdge(binx); - s[0] = xaxis->GetBinUpEdge(binx); - } - else - x[0] = xaxis->GetBinCenter(binx); - - - for ( biny = hyfirst; biny <= hylast; ++biny) { - if (useBinEdges) { - x[1] = yaxis->GetBinLowEdge(biny); - s[1] = yaxis->GetBinUpEdge(biny); - } - else - x[1] = yaxis->GetBinCenter(biny); - - for ( binz = hzfirst; binz <= hzlast; ++binz) { - if (useBinEdges) { - x[2] = zaxis->GetBinLowEdge(binz); - s[2] = zaxis->GetBinUpEdge(binz); - } - else - x[2] = zaxis->GetBinCenter(binz); - - // need to evaluate function to know about rejected points - // hugly but no other solutions - if (func != 0) { - TF1::RejectPoint(false); - (*func)( &x[0] ); // evaluate using stored function parameters - if (TF1::RejectedPoint() ) continue; + // reserve n for more efficient usage + //dv.Data().reserve(n); + + int hdim = hfit->GetDimension(); + int ndim = hdim; + // case of function dimension less than histogram + if (func !=0 && func->GetNdim() == hdim-1) ndim = hdim-1; + + assert( ndim > 0 ); + //typedef BinPoint::CoordData CoordData; + //CoordData x = CoordData( hfit->GetDimension() ); + + double x[3]; + double s[3]; + + int binx = 0; + int biny = 0; + int binz = 0; + + const TAxis *xaxis = hfit->GetXaxis(); + const TAxis *yaxis = hfit->GetYaxis(); + const TAxis *zaxis = hfit->GetZaxis(); + + for (int ix=0; ix<(int)range.Size(0); ++ix) { + hxfirst = xaxis->GetFirst(); + hxlast = xaxis->GetLast(); + HFitInterface::ExamineRange( xaxis, range(0, ix), hxfirst, hxlast ); //set hxfirst and hxlast for this ix range + + for (int iy=0; iy<=(int)range.Size(1); ++iy) { + if (iy > 0 && iy == (int)range.Size(1)) { continue; } + hyfirst = yaxis->GetFirst(); + hylast = yaxis->GetLast(); + if (ndim > 1 && range.Size(1) != 0) { + HFitInterface::ExamineRange( yaxis, range(1, iy), hyfirst, hylast ); //set hyfirst and hylast for this iy range + } + + for (int iz=0; iz<=(int)range.Size(2); ++iz) { + if (iz > 0 && iz == (int)range.Size(2)) { continue; } + hzfirst = zaxis->GetFirst(); + hzlast = zaxis->GetLast(); + if (ndim > 1 && range.Size(2) != 0) { + HFitInterface::ExamineRange( zaxis, range(2, iz), hzfirst, hzlast ); //set hzfirst and hzlast for this iz range } - - double value = hfit->GetBinContent(binx, biny, binz); - double error = hfit->GetBinError(binx, biny, binz); - if (!HFitInterface::AdjustError(fitOpt,error,value) ) continue; - - if (ndim == hdim -1) { - // case of fitting a function with dimension -1 - // point error is bin width y / sqrt(N) where N is nuber of entries in the bin - // normalization of error will be wrong - but they will be rescaled in the fit - if (hdim == 2) dv.Add( x, x[1], yaxis->GetBinWidth(biny) / error ); - if (hdim == 3) dv.Add( x, x[2], zaxis->GetBinWidth(binz) / error ); - } else { - dv.Add( x, value, error ); - if (useBinEdges) { - dv.AddBinUpEdge( s ); - } - } + int nPoints = (hxlast - hxfirst + 1) * (hylast - hyfirst + 1) * (hzlast - hzfirst + 1); + dv.Append(nPoints, ndim); + for ( binx = hxfirst; binx <= hxlast; ++binx) { + if (useBinEdges) { + x[0] = xaxis->GetBinLowEdge(binx); + s[0] = xaxis->GetBinUpEdge(binx); + } + else + x[0] = xaxis->GetBinCenter(binx); + + for ( biny = hyfirst; biny <= hylast; ++biny) { + if (useBinEdges) { + x[1] = yaxis->GetBinLowEdge(biny); + s[1] = yaxis->GetBinUpEdge(biny); + } + else + x[1] = yaxis->GetBinCenter(biny); + + for ( binz = hzfirst; binz <= hzlast; ++binz) { + if (useBinEdges) { + x[2] = zaxis->GetBinLowEdge(binz); + s[2] = zaxis->GetBinUpEdge(binz); + } + else + x[2] = zaxis->GetBinCenter(binz); + + // need to evaluate function to know about rejected points + // hugly but no other solutions + if (func != 0) { + TF1::RejectPoint(false); + (*func)( &x[0] ); // evaluate using stored function parameters + if (TF1::RejectedPoint() ) continue; + } + + double value = hfit->GetBinContent(binx, biny, binz); + double error = hfit->GetBinError(binx, biny, binz); + if (!HFitInterface::AdjustError(fitOpt,error,value) ) continue; + + if (ndim == hdim -1) { + // case of fitting a function with dimension -1 + // point error is bin width y / sqrt(N) where N is the number of entries in the bin + // normalization of error will be wrong - but they will be rescaled in the fit + if (hdim == 2) dv.Add( x, x[1], yaxis->GetBinWidth(biny) / error ); + if (hdim == 3) dv.Add( x, x[2], zaxis->GetBinWidth(binz) / error ); + } else { + dv.Add( x, value, error ); + if (useBinEdges) { + dv.AddBinUpEdge( s ); + } + } #ifdef DEBUG - std::cout << "bin " << binx << " add point " << x[0] << " " << hfit->GetBinContent(binx) << std::endl; + std::cout << "bin " << binx << " add point " << x[0] << " " << hfit->GetBinContent(binx) << std::endl; #endif - } // end loop on z bins - } // end loop on y bins - } // end loop on x axis + } // end loop on z bins + } // end loop on y bins + } // end loop on x axis + } // end loop on iz + } // end loop on iy + } //end loop on ix #ifdef DEBUG - std::cout << "THFitInterface::FillData: Hist FitData size is " << dv.Size() << std::endl; + std::cout << "THFitInterface::FillData: Hist FitData size is " << dv.Size() << std::endl; #endif -} - -//////////////////////////////////////////////////////////////////////////////// -/// Compute rough values of parameters for an exponential - -void InitExpo(const ROOT::Fit::BinData & data, TF1 * f1) -{ - unsigned int n = data.Size(); - if (n == 0) return; - - // find xmin and xmax of the data - double valxmin; - double xmin = *(data.GetPoint(0,valxmin)); - double xmax = xmin; - double valxmax = valxmin; - - for (unsigned int i = 1; i < n; ++ i) { - double val; - double x = *(data.GetPoint(i,val) ); - if (x < xmin) { - xmin = x; - valxmin = val; + } + + //////////////////////////////////////////////////////////////////////////////// + /// Compute rough values of parameters for an exponential + + void InitExpo(const ROOT::Fit::BinData & data, TF1 * f1) + { + unsigned int n = data.Size(); + if (n == 0) return; + + // find xmin and xmax of the data + double valxmin; + double xmin = *(data.GetPoint(0,valxmin)); + double xmax = xmin; + double valxmax = valxmin; + + for (unsigned int i = 1; i < n; ++ i) { + double val; + double x = *(data.GetPoint(i,val) ); + if (x < xmin) { + xmin = x; + valxmin = val; + } + else if (x > xmax) { + xmax = x; + valxmax = val; + } } - else if (x > xmax) { - xmax = x; - valxmax = val; - } - } - - // avoid negative values of valxmin/valxmax - if (valxmin <= 0 && valxmax > 0 ) valxmin = valxmax; - else if (valxmax <=0 && valxmin > 0) valxmax = valxmin; - else if (valxmin <=0 && valxmax <= 0) { valxmin = 1; valxmax = 1; } - - double slope = std::log( valxmax/valxmin) / (xmax - xmin); - double constant = std::log(valxmin) - slope * xmin; - f1->SetParameters(constant, slope); -} - - -//////////////////////////////////////////////////////////////////////////////// -/// Compute Initial values of parameters for a gaussian -/// derived from function H1InitGaus defined in TH1.cxx - -void InitGaus(const ROOT::Fit::BinData & data, TF1 * f1) -{ - - static const double sqrtpi = 2.506628; - - // - Compute mean value and RMS of the data - unsigned int n = data.Size(); - if (n == 0) return; - double sumx = 0; - double sumx2 = 0; - double allcha = 0; - double valmax = 0; - double rangex = data.Coords(n-1)[0] - data.Coords(0)[0]; - // to avoid binwidth = 0 set arbitrarly to 1 - double binwidth = 1; - if ( rangex > 0) binwidth = rangex; - double x0 = 0; - for (unsigned int i = 0; i < n; ++ i) { - double val; - double x = *(data.GetPoint(i,val) ); - sumx += val*x; - sumx2 += val*x*x; - allcha += val; - if (val > valmax) valmax = val; - if (i > 0) { - double dx = x - x0; - if (dx < binwidth) binwidth = dx; + + // avoid negative values of valxmin/valxmax + if (valxmin <= 0 && valxmax > 0 ) valxmin = valxmax; + else if (valxmax <=0 && valxmin > 0) valxmax = valxmin; + else if (valxmin <=0 && valxmax <= 0) { valxmin = 1; valxmax = 1; } + + double slope = std::log( valxmax/valxmin) / (xmax - xmin); + double constant = std::log(valxmin) - slope * xmin; + f1->SetParameters(constant, slope); + } + + + //////////////////////////////////////////////////////////////////////////////// + /// Compute Initial values of parameters for a gaussian + /// derived from function H1InitGaus defined in TH1.cxx + + void InitGaus(const ROOT::Fit::BinData & data, TF1 * f1) + { + + static const double sqrtpi = 2.506628; + + // - Compute mean value and RMS of the data + unsigned int n = data.Size(); + if (n == 0) return; + double sumx = 0; + double sumx2 = 0; + double allcha = 0; + double valmax = 0; + double rangex = data.Coords(n-1)[0] - data.Coords(0)[0]; + // to avoid binwidth = 0 set arbitrarly to 1 + double binwidth = 1; + if ( rangex > 0) binwidth = rangex; + double x0 = 0; + for (unsigned int i = 0; i < n; ++ i) { + double val; + double x = *(data.GetPoint(i,val) ); + sumx += val*x; + sumx2 += val*x*x; + allcha += val; + if (val > valmax) valmax = val; + if (i > 0) { + double dx = x - x0; + if (dx < binwidth) binwidth = dx; + } + x0 = x; } - x0 = x; - } - if (allcha <= 0) return; - double mean = sumx/allcha; - double rms = sumx2/allcha - mean*mean; + if (allcha <= 0) return; + double mean = sumx/allcha; + double rms = sumx2/allcha - mean*mean; - if (rms > 0) - rms = std::sqrt(rms); - else - rms = binwidth*n/4; + if (rms > 0) + rms = std::sqrt(rms); + else + rms = binwidth*n/4; - //if the distribution is really gaussian, the best approximation - //is binwidx*allcha/(sqrtpi*rms) - //However, in case of non-gaussian tails, this underestimates - //the normalisation constant. In this case the maximum value - //is a better approximation. - //We take the average of both quantities + //if the distribution is really gaussian, the best approximation + //is binwidx*allcha/(sqrtpi*rms) + //However, in case of non-gaussian tails, this underestimates + //the normalisation constant. In this case the maximum value + //is a better approximation. + //We take the average of both quantities -// printf("valmax %f other %f bw %f allcha %f rms %f \n",valmax, binwidth*allcha/(sqrtpi*rms), -// binwidth, allcha,rms ); + // printf("valmax %f other %f bw %f allcha %f rms %f \n",valmax, binwidth*allcha/(sqrtpi*rms), + // binwidth, allcha,rms ); - double constant = 0.5*(valmax+ binwidth*allcha/(sqrtpi*rms)); + double constant = 0.5*(valmax+ binwidth*allcha/(sqrtpi*rms)); - //In case the mean value is outside the histo limits and - //the RMS is bigger than the range, we take - // mean = center of bins - // rms = half range -// Double_t xmin = curHist->GetXaxis()->GetXmin(); -// Double_t xmax = curHist->GetXaxis()->GetXmax(); -// if ((mean < xmin || mean > xmax) && rms > (xmax-xmin)) { -// mean = 0.5*(xmax+xmin); -// rms = 0.5*(xmax-xmin); -// } + //In case the mean value is outside the histo limits and + //the RMS is bigger than the range, we take + // mean = center of bins + // rms = half range + // Double_t xmin = curHist->GetXaxis()->GetXmin(); + // Double_t xmax = curHist->GetXaxis()->GetXmax(); + // if ((mean < xmin || mean > xmax) && rms > (xmax-xmin)) { + // mean = 0.5*(xmax+xmin); + // rms = 0.5*(xmax-xmin); + // } - f1->SetParameter(0,constant); - f1->SetParameter(1,mean); - f1->SetParameter(2,rms); - f1->SetParLimits(2,0,10*rms); + f1->SetParameter(0,constant); + f1->SetParameter(1,mean); + f1->SetParameter(2,rms); + f1->SetParLimits(2,0,10*rms); #ifdef DEBUG - std::cout << "Gaussian initial par values" << constant << " " << mean << " " << rms << std::endl; + std::cout << "Gaussian initial par values" << constant << " " << mean << " " << rms << std::endl; #endif -} - -//////////////////////////////////////////////////////////////////////////////// -/// Compute Initial values of parameters for a gaussian -/// derived from function H1InitGaus defined in TH1.cxx - -void Init2DGaus(const ROOT::Fit::BinData & data, TF1 * f1) -{ - - static const double sqrtpi = 2.506628; - - // - Compute mean value and RMS of the data - unsigned int n = data.Size(); - if (n == 0) return; - double sumx = 0, sumy = 0; - double sumx2 = 0, sumy2 = 0; - double allcha = 0; - double valmax = 0; - double rangex = data.Coords(n-1)[0] - data.Coords(0)[0]; - double rangey = data.Coords(n-1)[1] - data.Coords(0)[1]; - // to avoid binwidthx = 0 set arbitrarly to 1 - double binwidthx = 1, binwidthy = 1; - if ( rangex > 0) binwidthx = rangex; - if ( rangey > 0) binwidthy = rangey; - double x0 = 0, y0 = 0; - for (unsigned int i = 0; i < n; ++i) { - double val; - const double *coords = data.GetPoint(i,val); - double x = coords[0], y = coords[1]; - sumx += val*x; - sumy += val*y; - sumx2 += val*x*x; - sumy2 += val*y*y; - allcha += val; - if (val > valmax) valmax = val; - if (i > 0) { - double dx = x - x0; - if (dx < binwidthx) binwidthx = dx; - double dy = y - y0; - if (dy < binwidthy) binwidthy = dy; + } + + //////////////////////////////////////////////////////////////////////////////// + /// Compute Initial values of parameters for a gaussian + /// derived from function H1InitGaus defined in TH1.cxx + + void Init2DGaus(const ROOT::Fit::BinData & data, TF1 * f1) + { + + static const double sqrtpi = 2.506628; + + // - Compute mean value and RMS of the data + unsigned int n = data.Size(); + if (n == 0) return; + double sumx = 0, sumy = 0; + double sumx2 = 0, sumy2 = 0; + double allcha = 0; + double valmax = 0; + double rangex = data.Coords(n-1)[0] - data.Coords(0)[0]; + double rangey = data.Coords(n-1)[1] - data.Coords(0)[1]; + // to avoid binwidthx = 0 set arbitrarly to 1 + double binwidthx = 1, binwidthy = 1; + if ( rangex > 0) binwidthx = rangex; + if ( rangey > 0) binwidthy = rangey; + double x0 = 0, y0 = 0; + for (unsigned int i = 0; i < n; ++i) { + double val; + const double *coords = data.GetPoint(i,val); + double x = coords[0], y = coords[1]; + sumx += val*x; + sumy += val*y; + sumx2 += val*x*x; + sumy2 += val*y*y; + allcha += val; + if (val > valmax) valmax = val; + if (i > 0) { + double dx = x - x0; + if (dx < binwidthx) binwidthx = dx; + double dy = y - y0; + if (dy < binwidthy) binwidthy = dy; + } + x0 = x; + y0 = y; } - x0 = x; - y0 = y; - } - if (allcha <= 0) return; - double meanx = sumx/allcha, meany = sumy/allcha; - double rmsx = sumx2/allcha - meanx*meanx; - double rmsy = sumy2/allcha - meany*meany; + if (allcha <= 0) return; + double meanx = sumx/allcha, meany = sumy/allcha; + double rmsx = sumx2/allcha - meanx*meanx; + double rmsy = sumy2/allcha - meany*meany; - if (rmsx > 0) - rmsx = std::sqrt(rmsx); - else - rmsx = binwidthx*n/4; + if (rmsx > 0) + rmsx = std::sqrt(rmsx); + else + rmsx = binwidthx*n/4; - if (rmsy > 0) - rmsy = std::sqrt(rmsy); - else - rmsy = binwidthy*n/4; + if (rmsy > 0) + rmsy = std::sqrt(rmsy); + else + rmsy = binwidthy*n/4; - //if the distribution is really gaussian, the best approximation - //is binwidx*allcha/(sqrtpi*rmsx) - //However, in case of non-gaussian tails, this underestimates - //the normalisation constant. In this case the maximum value - //is a better approximation. - //We take the average of both quantities + //if the distribution is really gaussian, the best approximation + //is binwidx*allcha/(sqrtpi*rmsx) + //However, in case of non-gaussian tails, this underestimates + //the normalisation constant. In this case the maximum value + //is a better approximation. + //We take the average of both quantities - double constant = 0.5 * (valmax+ binwidthx*allcha/(sqrtpi*rmsx))* - (valmax+ binwidthy*allcha/(sqrtpi*rmsy)); + double constant = 0.5 * (valmax+ binwidthx*allcha/(sqrtpi*rmsx))* + (valmax+ binwidthy*allcha/(sqrtpi*rmsy)); - f1->SetParameter(0,constant); - f1->SetParameter(1,meanx); - f1->SetParameter(2,rmsx); - f1->SetParLimits(2,0,10*rmsx); - f1->SetParameter(3,meany); - f1->SetParameter(4,rmsy); - f1->SetParLimits(4,0,10*rmsy); + f1->SetParameter(0,constant); + f1->SetParameter(1,meanx); + f1->SetParameter(2,rmsx); + f1->SetParLimits(2,0,10*rmsx); + f1->SetParameter(3,meany); + f1->SetParameter(4,rmsy); + f1->SetParLimits(4,0,10*rmsy); #ifdef DEBUG - std::cout << "2D Gaussian initial par values" - << constant << " " - << meanx << " " - << rmsx - << meany << " " - << rmsy - << std::endl; + std::cout << "2D Gaussian initial par values" + << constant << " " + << meanx << " " + << rmsx + << meany << " " + << rmsy + << std::endl; #endif -} - -// filling fit data from TGraph objects - -BinData::ErrorType GetDataType(const TGraph * gr, DataOptions & fitOpt) { - // get type of data for TGraph objects - double *ex = gr->GetEX(); - double *ey = gr->GetEY(); - double * eyl = gr->GetEYlow(); - double * eyh = gr->GetEYhigh(); - - - // default case for graphs (when they have errors) - BinData::ErrorType type = BinData::kValueError; - // if all errors are zero set option of using errors to 1 - if (fitOpt.fErrors1 || ( ey == 0 && ( eyl == 0 || eyh == 0 ) ) ) { - type = BinData::kNoError; - } - // need to treat case when all errors are zero - // note that by default fitOpt.fCoordError is true - else if ( ex != 0 && fitOpt.fCoordErrors) { - // check that all errors are not zero - int i = 0; - while (i < gr->GetN() && type != BinData::kCoordError) { - if (ex[i] > 0) type = BinData::kCoordError; - ++i; - } - } - // case of asymmetric errors (by default fAsymErrors is true) - else if ( ( eyl != 0 && eyh != 0) && fitOpt.fAsymErrors) { - // check also if that all errors are non zero's - int i = 0; - bool zeroErrorX = true; - bool zeroErrorY = true; - while (i < gr->GetN() && (zeroErrorX || zeroErrorY)) { - double e2X = ( gr->GetErrorXlow(i) + gr->GetErrorXhigh(i) ); - double e2Y = eyl[i] + eyh[i]; - zeroErrorX &= (e2X <= 0); - zeroErrorY &= (e2Y <= 0); - ++i; + } + + // filling fit data from TGraph objects + + BinData::ErrorType GetDataType(const TGraph * gr, DataOptions & fitOpt) { + // get type of data for TGraph objects + double *ex = gr->GetEX(); + double *ey = gr->GetEY(); + double * eyl = gr->GetEYlow(); + double * eyh = gr->GetEYhigh(); + + + // default case for graphs (when they have errors) + BinData::ErrorType type = BinData::kValueError; + // if all errors are zero set option of using errors to 1 + if (fitOpt.fErrors1 || ( ey == 0 && ( eyl == 0 || eyh == 0 ) ) ) { + type = BinData::kNoError; } - if (zeroErrorX && zeroErrorY) - type = BinData::kNoError; - else if (!zeroErrorX && zeroErrorY) - type = BinData::kCoordError; - else if (zeroErrorX && !zeroErrorY) { - type = BinData::kAsymError; - fitOpt.fCoordErrors = false; + // need to treat case when all errors are zero + // note that by default fitOpt.fCoordError is true + else if ( ex != 0 && fitOpt.fCoordErrors) { + // check that all errors are not zero + int i = 0; + while (i < gr->GetN() && type != BinData::kCoordError) { + if (ex[i] > 0) type = BinData::kCoordError; + ++i; + } } - else { - type = BinData::kAsymError; + // case of asymmetric errors (by default fAsymErrors is true) + else if ( ( eyl != 0 && eyh != 0) && fitOpt.fAsymErrors) { + // check also if that all errors are non zero's + int i = 0; + bool zeroErrorX = true; + bool zeroErrorY = true; + while (i < gr->GetN() && (zeroErrorX || zeroErrorY)) { + double e2X = ( gr->GetErrorXlow(i) + gr->GetErrorXhigh(i) ); + double e2Y = eyl[i] + eyh[i]; + zeroErrorX &= (e2X <= 0); + zeroErrorY &= (e2Y <= 0); + ++i; + } + if (zeroErrorX && zeroErrorY) + type = BinData::kNoError; + else if (!zeroErrorX && zeroErrorY) + type = BinData::kCoordError; + else if (zeroErrorX && !zeroErrorY) { + type = BinData::kAsymError; + fitOpt.fCoordErrors = false; + } + else { + type = BinData::kAsymError; + } } - } - - // need to look also a case when all errors in y are zero - if ( ey != 0 && type != BinData::kCoordError ) { - int i = 0; - bool zeroError = true; - while (i < gr->GetN() && zeroError) { - if (ey[i] > 0) zeroError = false;; - ++i; + + // need to look also a case when all errors in y are zero + if ( ey != 0 && type != BinData::kCoordError ) { + int i = 0; + bool zeroError = true; + while (i < gr->GetN() && zeroError) { + if (ey[i] > 0) zeroError = false;; + ++i; + } + if (zeroError) type = BinData::kNoError; } - if (zeroError) type = BinData::kNoError; - } #ifdef DEBUG - std::cout << "type is " << type << " graph type is " << gr->IsA()->GetName() << std::endl; + std::cout << "type is " << type << " graph type is " << gr->IsA()->GetName() << std::endl; #endif - return type; -} - -BinData::ErrorType GetDataType(const TGraph2D * gr, const DataOptions & fitOpt) { - // get type of data for TGraph2D object - double *ex = gr->GetEX(); - double *ey = gr->GetEY(); - double *ez = gr->GetEZ(); - - // default case for graphs (when they have errors) - BinData::ErrorType type = BinData::kValueError; - // if all errors are zero set option of using errors to 1 - if (fitOpt.fErrors1 || ez == 0 ) { - type = BinData::kNoError; - } - else if ( ex != 0 && ey!=0 && fitOpt.fCoordErrors) { - // check that all errors are not zero - int i = 0; - while (i < gr->GetN() && type != BinData::kCoordError) { - if (ex[i] > 0 || ey[i] > 0) type = BinData::kCoordError; - ++i; + return type; + } + + BinData::ErrorType GetDataType(const TGraph2D * gr, const DataOptions & fitOpt) { + // get type of data for TGraph2D object + double *ex = gr->GetEX(); + double *ey = gr->GetEY(); + double *ez = gr->GetEZ(); + + // default case for graphs (when they have errors) + BinData::ErrorType type = BinData::kValueError; + // if all errors are zero set option of using errors to 1 + if (fitOpt.fErrors1 || ez == 0 ) { + type = BinData::kNoError; + } + else if ( ex != 0 && ey!=0 && fitOpt.fCoordErrors) { + // check that all errors are not zero + int i = 0; + while (i < gr->GetN() && type != BinData::kCoordError) { + if (ex[i] > 0 || ey[i] > 0) type = BinData::kCoordError; + ++i; + } } - } #ifdef DEBUG - std::cout << "type is " << type << " graph2D type is " << gr->IsA()->GetName() << std::endl; + std::cout << "type is " << type << " graph2D type is " << gr->IsA()->GetName() << std::endl; #endif - return type; -} + return type; + } -void DoFillData ( BinData & dv, const TGraph * gr, BinData::ErrorType type, TF1 * func ) { - // internal method to do the actual filling of the data - // given a graph and a multigraph + void DoFillData ( BinData & dv, const TGraph * gr, BinData::ErrorType type, TF1 * func ) { + // internal method to do the actual filling of the data + // given a graph and a multigraph - // get fit option - DataOptions & fitOpt = dv.Opt(); + // get fit option + DataOptions & fitOpt = dv.Opt(); - int nPoints = gr->GetN(); - double *gx = gr->GetX(); - double *gy = gr->GetY(); + int nPoints = gr->GetN(); + double *gx = gr->GetX(); + double *gy = gr->GetY(); - const DataRange & range = dv.Range(); - bool useRange = ( range.Size(0) > 0); - double xmin = 0; - double xmax = 0; - range.GetRange(xmin,xmax); + const DataRange & range = dv.Range(); + bool useRange = ( range.Size(0) > 0); + double xmin = 0; + double xmax = 0; + range.GetRange(xmin,xmax); - dv.Initialize(nPoints,1, type); + dv.Initialize(nPoints,1, type); #ifdef DEBUG - std::cout << "DoFillData: graph npoints = " << nPoints << " type " << type << std::endl; - if (func) { - double a1,a2; func->GetRange(a1,a2); std::cout << "func range " << a1 << " " << a2 << std::endl; - } + std::cout << "DoFillData: graph npoints = " << nPoints << " type " << type << std::endl; + if (func) { + double a1,a2; func->GetRange(a1,a2); std::cout << "func range " << a1 << " " << a2 << std::endl; + } #endif - double x[1]; - for ( int i = 0; i < nPoints; ++i) { + double x[1]; + for ( int i = 0; i < nPoints; ++i) { - x[0] = gx[i]; + x[0] = gx[i]; - if (useRange && ( x[0] < xmin || x[0] > xmax) ) continue; + if (useRange && ( x[0] < xmin || x[0] > xmax) ) continue; - // need to evaluate function to know about rejected points - // hugly but no other solutions - if (func) { - TF1::RejectPoint(false); - (*func)( x ); // evaluate using stored function parameters - if (TF1::RejectedPoint() ) continue; - } + // need to evaluate function to know about rejected points + // hugly but no other solutions + if (func) { + TF1::RejectPoint(false); + (*func)( x ); // evaluate using stored function parameters + if (TF1::RejectedPoint() ) continue; + } - if (fitOpt.fErrors1) - dv.Add( gx[i], gy[i] ); + if (fitOpt.fErrors1) + dv.Add( gx[i], gy[i] ); - // for the errors use the getters by index to avoid cases when the arrays are zero - // (like in a case of a graph) - else if (type == BinData::kValueError) { - double errorY = gr->GetErrorY(i); - // should consider error = 0 as 1 ? Decide to skip points with zero errors - // in case want to keep points with error = 0 as errrors=1 need to set the option UseEmpty - if (!HFitInterface::AdjustError(fitOpt,errorY) ) continue; - dv.Add( gx[i], gy[i], errorY ); + // for the errors use the getters by index to avoid cases when the arrays are zero + // (like in a case of a graph) + else if (type == BinData::kValueError) { + double errorY = gr->GetErrorY(i); + // should consider error = 0 as 1 ? Decide to skip points with zero errors + // in case want to keep points with error = 0 as errrors=1 need to set the option UseEmpty + if (!HFitInterface::AdjustError(fitOpt,errorY) ) continue; + dv.Add( gx[i], gy[i], errorY ); #ifdef DEBUG - std::cout << "Point " << i << " " << gx[i] << " " << gy[i] << " " << errorY << std::endl; + std::cout << "Point " << i << " " << gx[i] << " " << gy[i] << " " << errorY << std::endl; #endif - } - else { // case use error in x or asym errors - double errorX = 0; - if (fitOpt.fCoordErrors) + } + else { // case use error in x or asym errors + double errorX = 0; + if (fitOpt.fCoordErrors) // shoulkd take combined average (sqrt(0.5(e1^2+e2^2)) or math average ? // gr->GetErrorX(i) returns combined average // use math average for same behaviour as before errorX = std::max( 0.5 * ( gr->GetErrorXlow(i) + gr->GetErrorXhigh(i) ) , 0. ) ; - // adjust error in y according to option - double errorY = std::max(gr->GetErrorY(i), 0.); - // we do not check the return value since we check later if error in X and Y is zero for skipping the point - HFitInterface::AdjustError(fitOpt, errorY); + // adjust error in y according to option + double errorY = std::max(gr->GetErrorY(i), 0.); + // we do not check the return value since we check later if error in X and Y is zero for skipping the point + HFitInterface::AdjustError(fitOpt, errorY); - // skip points with total error = 0 - if ( errorX <=0 && errorY <= 0 ) continue; + // skip points with total error = 0 + if ( errorX <=0 && errorY <= 0 ) continue; - if (type == BinData::kAsymError) { + if (type == BinData::kAsymError) { // asymmetric errors dv.Add( gx[i], gy[i], errorX, gr->GetErrorYlow(i), gr->GetErrorYhigh(i) ); - } - else { + } + else { // case symmetric Y errors dv.Add( gx[i], gy[i], errorX, errorY ); - } - } + } + } - } + } #ifdef DEBUG - std::cout << "TGraphFitInterface::FillData Graph FitData size is " << dv.Size() << std::endl; + std::cout << "TGraphFitInterface::FillData Graph FitData size is " << dv.Size() << std::endl; #endif -} - -void FillData(SparseData & dv, const TH1 * h1, TF1 * /*func*/) -{ - const int dim = h1->GetDimension(); - std::vector min(dim); - std::vector max(dim); - - int ncells = h1->GetNcells(); - for ( int i = 0; i < ncells; ++i ) { -// printf("i: %d; OF: %d; UF: %d; C: %f\n" -// , i -// , h1->IsBinOverflow(i) , h1->IsBinUnderflow(i) -// , h1->GetBinContent(i)); - if ( !( h1->IsBinOverflow(i) || h1->IsBinUnderflow(i) ) - && h1->GetBinContent(i)) - { - int x,y,z; - h1->GetBinXYZ(i, x, y, z); - -// std::cout << "FILLDATA: h1(" << i << ")" -// << "[" << h1->GetXaxis()->GetBinLowEdge(x) << "-" << h1->GetXaxis()->GetBinUpEdge(x) << "]"; -// if ( dim >= 2 ) -// std::cout << "[" << h1->GetYaxis()->GetBinLowEdge(y) << "-" << h1->GetYaxis()->GetBinUpEdge(y) << "]"; -// if ( dim >= 3 ) -// std::cout << "[" << h1->GetZaxis()->GetBinLowEdge(z) << "-" << h1->GetZaxis()->GetBinUpEdge(z) << "]"; - -// std::cout << h1->GetBinContent(i) << std::endl; - - min[0] = h1->GetXaxis()->GetBinLowEdge(x); - max[0] = h1->GetXaxis()->GetBinUpEdge(x); - if ( dim >= 2 ) - { - min[1] = h1->GetYaxis()->GetBinLowEdge(y); - max[1] = h1->GetYaxis()->GetBinUpEdge(y); - } - if ( dim >= 3 ) { - min[2] = h1->GetZaxis()->GetBinLowEdge(z); - max[2] = h1->GetZaxis()->GetBinUpEdge(z); - } - - dv.Add(min, max, h1->GetBinContent(i), h1->GetBinError(i)); - } - } -} - -void FillData(SparseData & dv, const THnBase * h1, TF1 * /*func*/) -{ - const int dim = h1->GetNdimensions(); - std::vector min(dim); - std::vector max(dim); - std::vector coord(dim); - - ULong64_t nEntries = h1->GetNbins(); - for ( ULong64_t i = 0; i < nEntries; i++ ) - { - double value = h1->GetBinContent( i, &coord[0] ); - if ( !value ) continue; - -// std::cout << "FILLDATA(SparseData): h1(" << i << ")"; - - // Exclude underflows and oveflows! (defect behaviour with the TH1*) - bool insertBox = true; - for ( int j = 0; j < dim && insertBox; ++j ) - { - TAxis* axis = h1->GetAxis(j); - if ( ( axis->GetBinLowEdge(coord[j]) < axis->GetXmin() ) || - ( axis->GetBinUpEdge(coord[j]) > axis->GetXmax() ) ) { - insertBox = false; - } - min[j] = h1->GetAxis(j)->GetBinLowEdge(coord[j]); - max[j] = h1->GetAxis(j)->GetBinUpEdge(coord[j]); - } - if ( !insertBox ) { -// std::cout << "NOT INSERTED!"<< std::endl; - continue; - } + } + + void FillData(SparseData & dv, const TH1 * h1, TF1 * /*func*/) + { + const int dim = h1->GetDimension(); + std::vector min(dim); + std::vector max(dim); + + int ncells = h1->GetNcells(); + for ( int i = 0; i < ncells; ++i ) { + // printf("i: %d; OF: %d; UF: %d; C: %f\n" + // , i + // , h1->IsBinOverflow(i) , h1->IsBinUnderflow(i) + // , h1->GetBinContent(i)); + if ( !( h1->IsBinOverflow(i) || h1->IsBinUnderflow(i) ) + && h1->GetBinContent(i)) + { + int x,y,z; + h1->GetBinXYZ(i, x, y, z); + + // std::cout << "FILLDATA: h1(" << i << ")" + // << "[" << h1->GetXaxis()->GetBinLowEdge(x) << "-" << h1->GetXaxis()->GetBinUpEdge(x) << "]"; + // if ( dim >= 2 ) + // std::cout << "[" << h1->GetYaxis()->GetBinLowEdge(y) << "-" << h1->GetYaxis()->GetBinUpEdge(y) << "]"; + // if ( dim >= 3 ) + // std::cout << "[" << h1->GetZaxis()->GetBinLowEdge(z) << "-" << h1->GetZaxis()->GetBinUpEdge(z) << "]"; + + // std::cout << h1->GetBinContent(i) << std::endl; + + min[0] = h1->GetXaxis()->GetBinLowEdge(x); + max[0] = h1->GetXaxis()->GetBinUpEdge(x); + if ( dim >= 2 ) + { + min[1] = h1->GetYaxis()->GetBinLowEdge(y); + max[1] = h1->GetYaxis()->GetBinUpEdge(y); + } + if ( dim >= 3 ) { + min[2] = h1->GetZaxis()->GetBinLowEdge(z); + max[2] = h1->GetZaxis()->GetBinUpEdge(z); + } -// for ( int j = 0; j < dim; ++j ) -// { -// std::cout << "[" << h1->GetAxis(j)->GetBinLowEdge(coord[j]) -// << "-" << h1->GetAxis(j)->GetBinUpEdge(coord[j]) << "]"; -// } -// std::cout << h1->GetBinContent(i) << std::endl; - - dv.Add(min, max, value, h1->GetBinError(i)); - } -} - -void FillData(BinData & dv, const THnBase * s1, TF1 * func) -{ - // Fill the Range of the THnBase - unsigned int const ndim = s1->GetNdimensions(); - std::vector xmin(ndim); - std::vector xmax(ndim); - for ( unsigned int i = 0; i < ndim; ++i ) { - TAxis* axis = s1->GetAxis(i); - xmin[i] = axis->GetXmin(); - xmax[i] = axis->GetXmax(); - } - - // Put default options, needed for the likelihood fitting of sparse - // data. - ROOT::Fit::DataOptions& dopt = dv.Opt(); - dopt.fUseEmpty = true; - // when using sparse data need to set option to use normalized bin volume, because sparse bins are merged together - //if (!dopt.fIntegral) dopt.fBinVolume = true; - dopt.fBinVolume = true; - dopt.fNormBinVolume = true; - - // Get the sparse data - ROOT::Fit::SparseData d(ndim, &xmin[0], &xmax[0]); - ROOT::Fit::FillData(d, s1, func); - -// std::cout << "FillData(BinData & dv, const THnBase * s1, TF1 * func) (1)" << std::endl; - - // Create the bin data from the sparse data - d.GetBinDataIntegral(dv); - -} - -void FillData ( BinData & dv, const TGraph * gr, TF1 * func ) { - // fill the data vector from a TGraph. Pass also the TF1 function which is - // needed in case to exclude points rejected by the function - assert(gr != 0); - - // get fit option - DataOptions & fitOpt = dv.Opt(); - - BinData::ErrorType type = GetDataType(gr,fitOpt); - // adjust option according to type - fitOpt.fErrors1 = (type == BinData::kNoError); - // set this if we want to have error=1 for points with zero errors (by default they are skipped) - // fitOpt.fUseEmpty = true; - - // use coordinate or asym errors in case option is set and type is consistent - fitOpt.fCoordErrors &= (type == BinData::kCoordError) || (type == BinData::kAsymError) ; - fitOpt.fAsymErrors &= (type == BinData::kAsymError); - - - // if data are filled already check if there are consistent - otherwise do nothing - if (dv.Size() > 0 && dv.NDim() == 1 ) { - // check if size is correct otherwise flag an errors - if ( dv.GetErrorType() != type ) { - Error("FillData","Inconsistent TGraph with previous data set- skip all graph data"); - return; + dv.Add(min, max, h1->GetBinContent(i), h1->GetBinError(i)); + } + } + } + + void FillData(SparseData & dv, const THnBase * h1, TF1 * /*func*/) + { + const int dim = h1->GetNdimensions(); + std::vector min(dim); + std::vector max(dim); + std::vector coord(dim); + + ULong64_t nEntries = h1->GetNbins(); + for ( ULong64_t i = 0; i < nEntries; i++ ) + { + double value = h1->GetBinContent( i, &coord[0] ); + if ( !value ) continue; + + // std::cout << "FILLDATA(SparseData): h1(" << i << ")"; + + // Exclude underflows and oveflows! (defect behaviour with the TH1*) + bool insertBox = true; + for ( int j = 0; j < dim && insertBox; ++j ) + { + TAxis* axis = h1->GetAxis(j); + if ( ( axis->GetBinLowEdge(coord[j]) < axis->GetXmin() ) || + ( axis->GetBinUpEdge(coord[j]) > axis->GetXmax() ) ) { + insertBox = false; + } + min[j] = h1->GetAxis(j)->GetBinLowEdge(coord[j]); + max[j] = h1->GetAxis(j)->GetBinUpEdge(coord[j]); + } + if ( !insertBox ) { + // std::cout << "NOT INSERTED!"<< std::endl; + continue; + } + + // for ( int j = 0; j < dim; ++j ) + // { + // std::cout << "[" << h1->GetAxis(j)->GetBinLowEdge(coord[j]) + // << "-" << h1->GetAxis(j)->GetBinUpEdge(coord[j]) << "]"; + // } + // std::cout << h1->GetBinContent(i) << std::endl; + + dv.Add(min, max, value, h1->GetBinError(i)); + } + } + + void FillData(BinData & dv, const THnBase * s1, TF1 * func) + { + // Fill the Range of the THnBase + unsigned int const ndim = s1->GetNdimensions(); + std::vector xmin(ndim); + std::vector xmax(ndim); + for ( unsigned int i = 0; i < ndim; ++i ) { + TAxis* axis = s1->GetAxis(i); + xmin[i] = axis->GetXmin(); + xmax[i] = axis->GetXmax(); } - } - - DoFillData(dv, gr, type, func); - -} - -void FillData ( BinData & dv, const TMultiGraph * mg, TF1 * func ) { - // fill the data vector from a TMultiGraph. Pass also the TF1 function which is - // needed in case to exclude points rejected by the function - assert(mg != 0); - TList * grList = mg->GetListOfGraphs(); - assert(grList != 0); + // Put default options, needed for the likelihood fitting of sparse + // data. + ROOT::Fit::DataOptions& dopt = dv.Opt(); + dopt.fUseEmpty = true; + // when using sparse data need to set option to use normalized bin volume, because sparse bins are merged together + //if (!dopt.fIntegral) dopt.fBinVolume = true; + dopt.fBinVolume = true; + dopt.fNormBinVolume = true; -#ifdef DEBUG -// grList->Print(); - TIter itr(grList, kIterBackward); - TObject *obj; - std::cout << "multi-graph list of graps: " << std::endl; - while ((obj = itr())) { - std::cout << obj->IsA()->GetName() << std::endl; - } + // Get the sparse data + ROOT::Fit::SparseData d(ndim, &xmin[0], &xmax[0]); + ROOT::Fit::FillData(d, s1, func); -#endif + // std::cout << "FillData(BinData & dv, const THnBase * s1, TF1 * func) (1)" << std::endl; - // get fit option - DataOptions & fitOpt = dv.Opt(); + // Create the bin data from the sparse data + d.GetBinDataIntegral(dv); - // loop on the graphs to get the data type (use maximum) - TIter next(grList); + } - BinData::ErrorType type = BinData::kNoError; - TGraph *gr = 0; - while ((gr = (TGraph*) next())) { - BinData::ErrorType t = GetDataType(gr,fitOpt); - if (t > type ) type = t; - } - // adjust option according to type - fitOpt.fErrors1 = (type == BinData::kNoError); - // use coordinate or asym errors in case option is set and type is consistent - fitOpt.fCoordErrors &= (type == BinData::kCoordError) || (type == BinData::kAsymError); - fitOpt.fAsymErrors &= (type == BinData::kAsymError); + void FillData ( BinData & dv, const TGraph * gr, TF1 * func ) { + // fill the data vector from a TGraph. Pass also the TF1 function which is + // needed in case to exclude points rejected by the function + assert(gr != 0); + // get fit option + DataOptions & fitOpt = dv.Opt(); -#ifdef DEBUG - std::cout << "Fitting MultiGraph of type " << type << std::endl; -#endif + BinData::ErrorType type = GetDataType(gr,fitOpt); + // adjust option according to type + fitOpt.fErrors1 = (type == BinData::kNoError); + // set this if we want to have error=1 for points with zero errors (by default they are skipped) + // fitOpt.fUseEmpty = true; - // fill the data now - next = grList; - while ((gr = (TGraph*) next())) { - DoFillData( dv, gr, type, func); - } + // use coordinate or asym errors in case option is set and type is consistent + fitOpt.fCoordErrors &= (type == BinData::kCoordError) || (type == BinData::kAsymError) ; + fitOpt.fAsymErrors &= (type == BinData::kAsymError); -#ifdef DEBUG - std::cout << "TGraphFitInterface::FillData MultiGraph FitData size is " << dv.Size() << std::endl; -#endif -} + // if data are filled already check if there are consistent - otherwise do nothing + if (dv.Size() > 0 && dv.NDim() == 1 ) { + // check if size is correct otherwise flag an errors + if ( dv.GetErrorType() != type ) { + Error("FillData","Inconsistent TGraph with previous data set- skip all graph data"); + return; + } + } -void FillData ( BinData & dv, const TGraph2D * gr, TF1 * func ) { - // fill the data vector from a TGraph2D. Pass also the TF1 function which is - // needed in case to exclude points rejected by the function - // in case of a pure TGraph - assert(gr != 0); + DoFillData(dv, gr, type, func); - // get fit option - DataOptions & fitOpt = dv.Opt(); - BinData::ErrorType type = GetDataType(gr,fitOpt); - // adjust option according to type - fitOpt.fErrors1 = (type == BinData::kNoError); - fitOpt.fCoordErrors = (type == BinData::kCoordError); - fitOpt.fAsymErrors = false; // a TGraph2D with asymmetric errors does not exist + } - int nPoints = gr->GetN(); - double *gx = gr->GetX(); - double *gy = gr->GetY(); - double *gz = gr->GetZ(); + void FillData ( BinData & dv, const TMultiGraph * mg, TF1 * func ) { + // fill the data vector from a TMultiGraph. Pass also the TF1 function which is + // needed in case to exclude points rejected by the function + assert(mg != 0); - // if all errors are zero set option of using errors to 1 - if ( gr->GetEZ() == 0) fitOpt.fErrors1 = true; + TList * grList = mg->GetListOfGraphs(); + assert(grList != 0); - double x[2]; - double ex[2]; +#ifdef DEBUG + // grList->Print(); + TIter itr(grList, kIterBackward); + TObject *obj; + std::cout << "multi-graph list of graps: " << std::endl; + while ((obj = itr())) { + std::cout << obj->IsA()->GetName() << std::endl; + } - // look at data range - const DataRange & range = dv.Range(); - bool useRangeX = ( range.Size(0) > 0); - bool useRangeY = ( range.Size(1) > 0); - double xmin = 0; - double xmax = 0; - double ymin = 0; - double ymax = 0; - range.GetRange(xmin,xmax,ymin,ymax); +#endif - dv.Initialize(nPoints,2, type); + // get fit option + DataOptions & fitOpt = dv.Opt(); - for ( int i = 0; i < nPoints; ++i) { + // loop on the graphs to get the data type (use maximum) + TIter next(grList); - x[0] = gx[i]; - x[1] = gy[i]; + BinData::ErrorType type = BinData::kNoError; + TGraph *gr = 0; + while ((gr = (TGraph*) next())) { + BinData::ErrorType t = GetDataType(gr,fitOpt); + if (t > type ) type = t; + } + // adjust option according to type + fitOpt.fErrors1 = (type == BinData::kNoError); + // use coordinate or asym errors in case option is set and type is consistent + fitOpt.fCoordErrors &= (type == BinData::kCoordError) || (type == BinData::kAsymError); + fitOpt.fAsymErrors &= (type == BinData::kAsymError); - //if (fitOpt.fUseRange && HFitInterface::IsPointOutOfRange(func, x) ) continue; - if (useRangeX && ( x[0] < xmin || x[0] > xmax) ) continue; - if (useRangeY && ( x[1] < ymin || x[1] > ymax) ) continue; - // need to evaluate function to know about rejected points - // hugly but no other solutions - if (func) { - TF1::RejectPoint(false); - (*func)( x ); // evaluate using stored function parameters - if (TF1::RejectedPoint() ) continue; - } +#ifdef DEBUG + std::cout << "Fitting MultiGraph of type " << type << std::endl; +#endif - if (type == BinData::kNoError) { - dv.Add( x, gz[i] ); - continue; + // fill the data now + next = grList; + while ((gr = (TGraph*) next())) { + DoFillData( dv, gr, type, func); } - double errorZ = gr->GetErrorZ(i); - if (!HFitInterface::AdjustError(fitOpt,errorZ) ) continue; +#ifdef DEBUG + std::cout << "TGraphFitInterface::FillData MultiGraph FitData size is " << dv.Size() << std::endl; +#endif - if (type == BinData::kValueError) { - dv.Add( x, gz[i], errorZ ); - } - else if (type == BinData::kCoordError) { // case use error in coordinates (x and y) - ex[0] = std::max(gr->GetErrorX(i), 0.); - ex[1] = std::max(gr->GetErrorY(i), 0.); - dv.Add( x, gz[i], ex, errorZ ); - } - else - assert(0); // should not go here + } + + void FillData ( BinData & dv, const TGraph2D * gr, TF1 * func ) { + // fill the data vector from a TGraph2D. Pass also the TF1 function which is + // needed in case to exclude points rejected by the function + // in case of a pure TGraph + assert(gr != 0); + + // get fit option + DataOptions & fitOpt = dv.Opt(); + BinData::ErrorType type = GetDataType(gr,fitOpt); + // adjust option according to type + fitOpt.fErrors1 = (type == BinData::kNoError); + fitOpt.fCoordErrors = (type == BinData::kCoordError); + fitOpt.fAsymErrors = false; // a TGraph2D with asymmetric errors does not exist + + int nPoints = gr->GetN(); + double *gx = gr->GetX(); + double *gy = gr->GetY(); + double *gz = gr->GetZ(); + + // if all errors are zero set option of using errors to 1 + if ( gr->GetEZ() == 0) fitOpt.fErrors1 = true; + + double x[2]; + double ex[2]; + + // look at data range + const DataRange & range = dv.Range(); + bool useRangeX = ( range.Size(0) > 0); + bool useRangeY = ( range.Size(1) > 0); + double xmin = 0; + double xmax = 0; + double ymin = 0; + double ymax = 0; + range.GetRange(xmin,xmax,ymin,ymax); + + dv.Initialize(nPoints,2, type); + + for ( int i = 0; i < nPoints; ++i) { + + x[0] = gx[i]; + x[1] = gy[i]; + + //if (fitOpt.fUseRange && HFitInterface::IsPointOutOfRange(func, x) ) continue; + if (useRangeX && ( x[0] < xmin || x[0] > xmax) ) continue; + if (useRangeY && ( x[1] < ymin || x[1] > ymax) ) continue; + + // need to evaluate function to know about rejected points + // hugly but no other solutions + if (func) { + TF1::RejectPoint(false); + (*func)( x ); // evaluate using stored function parameters + if (TF1::RejectedPoint() ) continue; + } + + if (type == BinData::kNoError) { + dv.Add( x, gz[i] ); + continue; + } + + double errorZ = gr->GetErrorZ(i); + if (!HFitInterface::AdjustError(fitOpt,errorZ) ) continue; + + if (type == BinData::kValueError) { + dv.Add( x, gz[i], errorZ ); + } + else if (type == BinData::kCoordError) { // case use error in coordinates (x and y) + ex[0] = std::max(gr->GetErrorX(i), 0.); + ex[1] = std::max(gr->GetErrorY(i), 0.); + dv.Add( x, gz[i], ex, errorZ ); + } + else + assert(0); // should not go here #ifdef DEBUG - std::cout << "Point " << i << " " << gx[i] << " " << gy[i] << " " << errorZ << std::endl; + std::cout << "Point " << i << " " << gx[i] << " " << gy[i] << " " << errorZ << std::endl; #endif - } + } #ifdef DEBUG - std::cout << "THFitInterface::FillData Graph2D FitData size is " << dv.Size() << std::endl; + std::cout << "THFitInterface::FillData Graph2D FitData size is " << dv.Size() << std::endl; #endif -} - - -// confidence intervals -bool GetConfidenceIntervals(const TH1 * h1, const ROOT::Fit::FitResult & result, TGraphErrors * gr, double cl ) { - if (h1->GetDimension() != 1) { - Error("GetConfidenceIntervals","Invalid object used for storing confidence intervals"); - return false; - } - // fill fit data sets with points to estimate cl. - BinData d; - FillData(d,h1,0); - gr->Set(d.NPoints() ); - double * ci = gr->GetEY(); // make CL values error of the graph - result.GetConfidenceIntervals(d,ci,cl); - // put function value as abscissa of the graph - for (unsigned int ipoint = 0; ipoint < d.NPoints(); ++ipoint) { - const double * x = d.Coords(ipoint); - const ROOT::Math::IParamMultiFunction * func = result.FittedFunction(); - gr->SetPoint(ipoint, x[0], (*func)(x) ); - } - return true; -} - - -} // end namespace Fit + } + + + // confidence intervals + bool GetConfidenceIntervals(const TH1 * h1, const ROOT::Fit::FitResult & result, TGraphErrors * gr, double cl ) { + if (h1->GetDimension() != 1) { + Error("GetConfidenceIntervals","Invalid object used for storing confidence intervals"); + return false; + } + // fill fit data sets with points to estimate cl. + BinData d; + FillData(d,h1,0); + gr->Set(d.NPoints() ); + double * ci = gr->GetEY(); // make CL values error of the graph + result.GetConfidenceIntervals(d,ci,cl); + // put function value as abscissa of the graph + for (unsigned int ipoint = 0; ipoint < d.NPoints(); ++ipoint) { + const double * x = d.Coords(ipoint); + const ROOT::Math::IParamMultiFunction * func = result.FittedFunction(); + gr->SetPoint(ipoint, x[0], (*func)(x) ); + } + return true; + } + + + } // end namespace Fit } // end namespace ROOT diff --git a/math/mathcore/inc/Fit/DataRange.h b/math/mathcore/inc/Fit/DataRange.h index 19d0a387802aa..0110eb41f0776 100644 --- a/math/mathcore/inc/Fit/DataRange.h +++ b/math/mathcore/inc/Fit/DataRange.h @@ -222,7 +222,7 @@ class DataRange { internal function to remove all the existing ranges between xmin and xmax called when a new range is inserted */ - void CleanRangeSet(unsigned int icoord, double xmin, double xmax); + void CleanRangeSet(unsigned int icoord, double &xmin, double &xmax); // get the full range (-inf, +inf) static void GetInfRange(double &x1, double &x2); diff --git a/math/mathcore/src/DataRange.cxx b/math/mathcore/src/DataRange.cxx index 2b4ad06b2d263..feeb07813985a 100644 --- a/math/mathcore/src/DataRange.cxx +++ b/math/mathcore/src/DataRange.cxx @@ -113,8 +113,8 @@ void DataRange::AddRange(unsigned int icoord , double xmin, double xmax ) { // case of an already existing range // need to establish a policy (use OR or AND ) - CleanRangeSet(icoord,xmin,xmax); - // add the new one + CleanRangeSet(icoord, xmin, xmax); + // add the new one rs.push_back(std::make_pair(xmin,xmax) ); // sort range in increasing values of xmax std::sort( rs.begin(), rs.end() , lessRange); @@ -161,24 +161,38 @@ void DataRange::Clear(unsigned int icoord ) { fRanges[icoord].clear(); } - -void DataRange::CleanRangeSet(unsigned int icoord, double xmin, double xmax) { - // remove all the existing ranges between xmin and xmax - // called when a new range is inserted - - // loop on existing ranges - RangeSet & ranges = fRanges[icoord]; - for (RangeSet::iterator itr = ranges.begin(); itr != ranges.end(); ++itr) { - // delete included ranges - if ( itr->first >= xmin && itr->second <= xmax) { - itr = ranges.erase(itr); - // itr goes to next element, so go back before adding - --itr; - } - } - +void DataRange::CleanRangeSet(unsigned int icoord, double &xmin, double &xmax) { + // loop on existing ranges + RangeSet & ranges = fRanges[icoord]; + for (RangeSet::iterator itr = ranges.begin(); itr != ranges.end(); ++itr) { + // delete included ranges + if ( itr->first >= xmin && itr->second <= xmax) { + itr = ranges.erase(itr); + // itr goes to next element, so go back before adding + --itr; + continue; + } + + if (xmin >= itr->first && xmax <= itr->second) { + xmin = itr->first; + xmax = itr->second; + itr = ranges.erase(itr); + // itr goes to next element, so go back before adding + --itr; + continue; + } + if ( (itr->first <= xmin && xmin <= itr->second) + || (itr->first <= xmax && xmax <= itr->second) ) { + xmin = std::min(itr->first, xmin); + xmax = std::max(itr->second, xmax); + itr = ranges.erase(itr); + --itr; + } + } } + + void DataRange::GetInfRange(double &xmin, double &xmax) { // get the full range [-inf, +inf] for xmin and xmax xmin = -std::numeric_limits::infinity(); diff --git a/math/mathcore/test/CMakeLists.txt b/math/mathcore/test/CMakeLists.txt index 95af3ed733073..97e2ea860c2b1 100644 --- a/math/mathcore/test/CMakeLists.txt +++ b/math/mathcore/test/CMakeLists.txt @@ -84,7 +84,9 @@ ROOT_ADD_GTEST(MulmodUnitNoInt128 mulmod_noint128.cxx) ROOT_ADD_GTEST(RanluxLCGUnit ranlux_lcg.cxx) ROOT_ADD_GTEST(RanluxppEngineTests RanluxppEngine.cxx - LIBRARIES Core MathCore) + LIBRARIES Core MathCore) + +ROOT_ADD_GTEST(testRange test_Range.cxx LIBRARIES MathCore) if(veccore AND vc) ROOT_ADD_GTEST(VectorizedTMathUnit testVectorizedTMath.cxx diff --git a/math/mathcore/test/test_Range.cxx b/math/mathcore/test/test_Range.cxx new file mode 100644 index 0000000000000..ac97e55cfdf80 --- /dev/null +++ b/math/mathcore/test/test_Range.cxx @@ -0,0 +1,57 @@ +#include "gtest/gtest.h" + +#include "Fit/DataRange.h" + +//Overlap Range +TEST(Range, Overlap) +{ + ROOT::Fit::DataRange range; + range.AddRange(0,5); + range.AddRange(2,3); + + EXPECT_EQ(range.Size(), 1); + + EXPECT_EQ(range(0,0).first, 0); + EXPECT_EQ(range(0,0).second, 5); + + range.AddRange(-1,6); + EXPECT_EQ(range.Size(), 1); + + EXPECT_EQ(range(0,0).first, -1); + EXPECT_EQ(range(0,0).second, 6); + + range.AddRange(-2,4); + EXPECT_EQ(range.Size(), 1); + + EXPECT_EQ(range(0,0).first, -2); + EXPECT_EQ(range(0,0).second, 6); + + range.AddRange(5,7); + EXPECT_EQ(range.Size(), 1); + + EXPECT_EQ(range(0,0).first, -2); + EXPECT_EQ(range(0,0).second, 7); + + range.AddRange(20,25); + EXPECT_EQ(range.Size(), 2); + + EXPECT_EQ(range(0,0).first, -2); + EXPECT_EQ(range(0,0).second, 7); + + EXPECT_EQ(range(0,1).first, 20); + EXPECT_EQ(range(0,1).second, 25); + + range.AddRange(24,26); + EXPECT_EQ(range.Size(), 2); + EXPECT_EQ(range(0,1).first, 20); + EXPECT_EQ(range(0,1).second, 26); + + range.AddRange(19,20); + EXPECT_EQ(range(0,1).first, 19); + EXPECT_EQ(range(0,1).second, 26); + + range.AddRange(6,20); + EXPECT_EQ(range.Size(), 1); + EXPECT_EQ(range(0,0).first, -2); + EXPECT_EQ(range(0,0).second, 26); +}