Functions
double	Interpolate (const TH1 *histo, const double x)
double	Interpolate (const TH1 *histo, const double x, const double y)
double	Interpolate2D (const TH1 *histo, const double x, const double y, const int xAxis, const int yAxis, const int otherDimBin)
	Helper function for edge cases in 3D interpolation by interpolating in 2D and keeping the third axis axis (edge) fixed The input parameters "xAxis" and "yAxis" indicate which axes should be used for the 2D interpolation: x-axis: 1, y-axis: 2, z-axis: 3 The input parameter "otherDimBin" is the bin on the third axis which will be kept constant in the interpolation.
double	Interpolate (const TH1 *histo, const double x, const double y, const double z)
	Function for 3D interpolation to account for edge cases in TH3::Interpolate.
Int_t	FindBin (const TAxis *axis, const double x)

Function Documentation

◆ FindBin()

Int_t RootHelpers::FindBin	(	const TAxis *	axis,
		const double	x )

Definition at line 14 of file JetCalibTools/Root/RootHelpers.cxx.

{
    // Copied from ROOT directly and trivially modified, all credit to ROOT authors of TAxis FindBin method
    // This is done because I want a const version of bin finding (no expanding on under/overflow)
    const double fXmin = axis->GetXmin();
    const double fXmax = axis->GetXmax();
    const Int_t fNbins = axis->GetNbins();
    const TArrayD* fXbins = axis->GetXbins();
    Int_t bin;
    if (x < fXmin) {              //*-* underflow
       bin = 0;
    } else  if ( !(x < fXmax)) {     //*-* overflow  (note the way to catch NaN
       bin = fNbins+1;
    } else {
       if (!fXbins->fN) {        //*-* fix bins
          bin = 1 + int (fNbins*(x-fXmin)/(fXmax-fXmin) );
       } else {                  //*-* variable bin sizes
          //for (bin =1; x >= fXbins->fArray[bin]; bin++);
          bin = 1 + TMath::BinarySearch(fXbins->fN,fXbins->fArray,x);
       }
    }
    return bin;
}

◆ Interpolate() [1/3]

double RootHelpers::Interpolate	(	const TH1 *	histo,
		const double	x )

Definition at line 39 of file JetCalibTools/Root/RootHelpers.cxx.

{
    // Copied from ROOT directly and trivially modified, all credit to ROOT authors of TH1, TH2, and TH3 Interpolate methods
    // This is done because I want a const version of interpolation, and none of the methods require modification of the histogram
    // Probable reason is that FindBin isn't const, but there should be a const version...
    Int_t xbin = RootHelpers::FindBin(histo->GetXaxis(),x);
    Double_t x0,x1,y0,y1;
    
    if(x<=histo->GetBinCenter(1)) {
       return histo->GetBinContent(1);
    } else if(x>=histo->GetBinCenter(histo->GetNbinsX())) {
       return histo->GetBinContent(histo->GetNbinsX());
    } else {
       if(x<=histo->GetBinCenter(xbin)) {
          y0 = histo->GetBinContent(xbin-1);
          x0 = histo->GetBinCenter(xbin-1);
          y1 = histo->GetBinContent(xbin);
          x1 = histo->GetBinCenter(xbin);
       } else {
          y0 = histo->GetBinContent(xbin);
          x0 = histo->GetBinCenter(xbin);
          y1 = histo->GetBinContent(xbin+1);
          x1 = histo->GetBinCenter(xbin+1);
       }
    return y0 + (x-x0)*((y1-y0)/(x1-x0));
    }
}

◆ Interpolate() [2/3]

double RootHelpers::Interpolate	(	const TH1 *	histo,
		const double	x,
		const double	y )

Definition at line 67 of file JetCalibTools/Root/RootHelpers.cxx.

{
    // Call the unified method for consistency
    return RootHelpers::Interpolate2D(histo,x,y,1,2,-1);
}

◆ Interpolate() [3/3]

double RootHelpers::Interpolate	(	const TH1 *	histo,
		const double	x,
		const double	y,
		const double	z )

Function for 3D interpolation to account for edge cases in TH3::Interpolate.

Definition at line 233 of file JetCalibTools/Root/RootHelpers.cxx.

{
    // Copied from ROOT directly and trivially modified, all credit to ROOT authors of TH1, TH2, and TH3 Interpolate methods
    // This is done because I want a const version of interpolation, and none of the methods require modification of the histogram
    // Probable reason is that FindBin isn't const, but there should be a const version...
    const TAxis* fXaxis = histo->GetXaxis();
    const TAxis* fYaxis = histo->GetYaxis();
    const TAxis* fZaxis = histo->GetZaxis();
 
    // Find the bin by bin edges
    Int_t ubx = RootHelpers::FindBin(fXaxis,x);
    Int_t uby = RootHelpers::FindBin(fYaxis,y);
    Int_t ubz = RootHelpers::FindBin(fZaxis,z);
    
    // Check if the value(s) are outside of the bin range(s)
    if ( ubx < 1 || ubx > histo->GetNbinsX() || uby < 1 || uby > histo->GetNbinsY() || ubz < 1 || ubz > histo->GetNbinsZ() )
    {
       histo->Error("Interpolate","Cannot interpolate outside histogram domain. (x: %f vs [%f,%f], y: %f vs [%f,%f], z: %f vs [%f,%f])",x,fXaxis->GetBinLowEdge(1),fXaxis->GetBinLowEdge(histo->GetNbinsX()+1),y,fYaxis->GetBinLowEdge(1),fYaxis->GetBinLowEdge(histo->GetNbinsY()+1),z,fZaxis->GetBinLowEdge(1),fZaxis->GetBinLowEdge(histo->GetNbinsZ()+1));
       return 0;
    }
 
    // Now switch from bin edges to bin centres
    // Note that we want to support edge cases, so it is possible that ub* == ob*
    // This functionality is not in original ROOT TH3::Interpolate()
    // This functionality is inspired by TH2::Interpolate()
    Int_t obx = ubx + 1;
    Int_t oby = uby + 1;
    Int_t obz = ubz + 1;
    
    // Calculate distance weights before checking under/overflow bins
    Double_t xw = fXaxis->GetBinCenter(obx) - fXaxis->GetBinCenter(ubx);
    Double_t yw = fYaxis->GetBinCenter(oby) - fYaxis->GetBinCenter(uby);
    Double_t zw = fZaxis->GetBinCenter(obz) - fZaxis->GetBinCenter(ubz);
    
    if (x < fXaxis->GetBinCenter(ubx)) { ubx -= 1; obx -= 1; }
    if (ubx < 1) ubx = 1;
    if (obx > histo->GetNbinsX()) obx = histo->GetNbinsX();
 
    if (y < fYaxis->GetBinCenter(uby)) { uby -= 1; oby -= 1; }
    if (uby < 1) uby = 1;
    if (oby > histo->GetNbinsY()) oby = histo->GetNbinsY();
 
    if (z < fZaxis->GetBinCenter(ubz)) { ubz -= 1; obz -= 1; }
    if (ubz < 1) ubz = 1;
    if (obz > histo->GetNbinsZ()) obz = histo->GetNbinsZ();
 
    // Edge cases were tried with weights set including the under/overflow bins (to follow what TH2::Interpolate() does for boundaries)
    // In some cases, it performed quite poorly
    // Tests of switching to 2D interpolation with the third dimension fixed appeared to work much better
    // Thus, the below is now a switch to bilinear interpolation when bin(s) are equal in trilinear interpolation
    
    // The below code block contains added functionality with respect to standard ROOT TH3::Interpolate()
    // Fall back on bilinear interpolation
    {
        if (ubz == obz)
            return RootHelpers::Interpolate2D(histo,x,y,1,2,ubz);
        else if (uby == oby)
            return RootHelpers::Interpolate2D(histo,x,z,1,3,uby);
        else if (ubx == obx)
            return RootHelpers::Interpolate2D(histo,y,z,2,3,ubx);
        
    }
    
    // Not a boundary case, resume normal ROOT::TH3::Interpolate()
    Double_t xd = (x - fXaxis->GetBinCenter(ubx)) / xw;
    Double_t yd = (y - fYaxis->GetBinCenter(uby)) / yw;
    Double_t zd = (z - fZaxis->GetBinCenter(ubz)) / zw;
    
    
    Double_t v[] = { histo->GetBinContent( ubx, uby, ubz ), histo->GetBinContent( ubx, uby, obz ),
                     histo->GetBinContent( ubx, oby, ubz ), histo->GetBinContent( ubx, oby, obz ),
                     histo->GetBinContent( obx, uby, ubz ), histo->GetBinContent( obx, uby, obz ),
                     histo->GetBinContent( obx, oby, ubz ), histo->GetBinContent( obx, oby, obz ) };
    
    
    Double_t i1 = v[0] * (1 - zd) + v[1] * zd;
    Double_t i2 = v[2] * (1 - zd) + v[3] * zd;
    Double_t j1 = v[4] * (1 - zd) + v[5] * zd;
    Double_t j2 = v[6] * (1 - zd) + v[7] * zd;
    
    
    Double_t w1 = i1 * (1 - yd) + i2 * yd;
    Double_t w2 = j1 * (1 - yd) + j2 * yd;
    
    
    Double_t result = w1 * (1 - xd) + w2 * xd;
    
    return result;
}

◆ Interpolate2D()

double RootHelpers::Interpolate2D	(	const TH1 *	histo,
		const double	x,
		const double	y,
		const int	xAxis,
		const int	yAxis,
		const int	otherDimBin )

Helper function for edge cases in 3D interpolation by interpolating in 2D and keeping the third axis axis (edge) fixed The input parameters "xAxis" and "yAxis" indicate which axes should be used for the 2D interpolation: x-axis: 1, y-axis: 2, z-axis: 3 The input parameter "otherDimBin" is the bin on the third axis which will be kept constant in the interpolation.

Definition at line 73 of file JetCalibTools/Root/RootHelpers.cxx.