ATLAS Offline Software
RootHelpers.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "JetCalibTools/RootHelpers.h"
6 
7 #include "TH1.h"
8 #include "TMath.h"
9 
10 namespace RootHelpers
11 {
12 
13 
14 Int_t FindBin(const TAxis* axis, const double x)
15 {
16  // Copied from ROOT directly and trivially modified, all credit to ROOT authors of TAxis FindBin method
17  // This is done because I want a const version of bin finding (no expanding on under/overflow)
18  const double fXmin = axis->GetXmin();
19  const double fXmax = axis->GetXmax();
20  const Int_t fNbins = axis->GetNbins();
21  const TArrayD* fXbins = axis->GetXbins();
22  Int_t bin;
23  if (x < fXmin) { //*-* underflow
24  bin = 0;
25  } else if ( !(x < fXmax)) { //*-* overflow (note the way to catch NaN
26  bin = fNbins+1;
27  } else {
28  if (!fXbins->fN) { //*-* fix bins
29  bin = 1 + int (fNbins*(x-fXmin)/(fXmax-fXmin) );
30  } else { //*-* variable bin sizes
31  //for (bin =1; x >= fXbins->fArray[bin]; bin++);
32  bin = 1 + TMath::BinarySearch(fXbins->fN,fXbins->fArray,x);
33  }
34  }
35  return bin;
36 }
37 
38 
39 double Interpolate(const TH1* histo, const double x)
40 {
41  // Copied from ROOT directly and trivially modified, all credit to ROOT authors of TH1, TH2, and TH3 Interpolate methods
42  // This is done because I want a const version of interpolation, and none of the methods require modification of the histogram
43  // Probable reason is that FindBin isn't const, but there should be a const version...
44  Int_t xbin = RootHelpers::FindBin(histo->GetXaxis(),x);
45  Double_t x0,x1,y0,y1;
46 
47  if(x<=histo->GetBinCenter(1)) {
48  return histo->GetBinContent(1);
49  } else if(x>=histo->GetBinCenter(histo->GetNbinsX())) {
50  return histo->GetBinContent(histo->GetNbinsX());
51  } else {
52  if(x<=histo->GetBinCenter(xbin)) {
53  y0 = histo->GetBinContent(xbin-1);
54  x0 = histo->GetBinCenter(xbin-1);
55  y1 = histo->GetBinContent(xbin);
56  x1 = histo->GetBinCenter(xbin);
57  } else {
58  y0 = histo->GetBinContent(xbin);
59  x0 = histo->GetBinCenter(xbin);
60  y1 = histo->GetBinContent(xbin+1);
61  x1 = histo->GetBinCenter(xbin+1);
62  }
63  return y0 + (x-x0)*((y1-y0)/(x1-x0));
64  }
65 }
66 
67 double Interpolate(const TH1* histo, const double x, const double y)
68 {
69  // Call the unified method for consistency
70  return RootHelpers::Interpolate2D(histo,x,y,1,2,-1);
71 }
72 
73 double Interpolate2D(const TH1* histo, const double x, const double y, const int xAxis, const int yAxis, const int otherDimBin)
74 {
75  // Copied from ROOT directly and trivially modified, all credit to ROOT authors of TH1, TH2, and TH3 Interpolate methods
76  // This is done because I want a const version of interpolation, and none of the methods require modification of the histogram
77  // Probable reason is that FindBin isn't const, but there should be a const version...
78  const TAxis* fXaxis = xAxis == 1 ? histo->GetXaxis() : xAxis == 2 ? histo->GetYaxis() : xAxis == 3 ? histo->GetZaxis() : nullptr;
79  const TAxis* fYaxis = yAxis == 1 ? histo->GetXaxis() : yAxis == 2 ? histo->GetYaxis() : yAxis == 3 ? histo->GetZaxis() : nullptr;
80 
81  if (!fXaxis || !fYaxis)
82  {
83  histo->Error("Interpolate2D","Failed to parse axes from inputs");
84  return 0;
85  }
86 
87  Double_t f=0;
88  Double_t x1=0,x2=0,y1=0,y2=0;
89  Double_t dx,dy;
90  Int_t bin_x = RootHelpers::FindBin(fXaxis,x);
91  Int_t bin_y = RootHelpers::FindBin(fYaxis,y);
92  if(bin_x<1 || bin_x>fXaxis->GetNbins() || bin_y<1 || bin_y>fYaxis->GetNbins()) {
93  histo->Error("Interpolate","Cannot interpolate outside histogram domain. (x: %f vs [%f,%f], y: %f vs [%f,%f])",x,fXaxis->GetBinLowEdge(1),fXaxis->GetBinLowEdge(fXaxis->GetNbins()+1),y,fYaxis->GetBinLowEdge(1),fYaxis->GetBinLowEdge(fYaxis->GetNbins()+1));
94  return 0;
95  }
96  Int_t quadrant = 0; // CCW from UR 1,2,3,4
97  // which quadrant of the bin (bin_P) are we in?
98  dx = fXaxis->GetBinUpEdge(bin_x)-x;
99  dy = fYaxis->GetBinUpEdge(bin_y)-y;
100  if (dx<=fXaxis->GetBinWidth(bin_x)/2 && dy<=fYaxis->GetBinWidth(bin_y)/2)
101  quadrant = 1; // upper right
102  if (dx>fXaxis->GetBinWidth(bin_x)/2 && dy<=fYaxis->GetBinWidth(bin_y)/2)
103  quadrant = 2; // upper left
104  if (dx>fXaxis->GetBinWidth(bin_x)/2 && dy>fYaxis->GetBinWidth(bin_y)/2)
105  quadrant = 3; // lower left
106  if (dx<=fXaxis->GetBinWidth(bin_x)/2 && dy>fYaxis->GetBinWidth(bin_y)/2)
107  quadrant = 4; // lower right
108  switch(quadrant) {
109  case 1:
110  x1 = fXaxis->GetBinCenter(bin_x);
111  y1 = fYaxis->GetBinCenter(bin_y);
112  x2 = fXaxis->GetBinCenter(bin_x+1);
113  y2 = fYaxis->GetBinCenter(bin_y+1);
114  break;
115  case 2:
116  x1 = fXaxis->GetBinCenter(bin_x-1);
117  y1 = fYaxis->GetBinCenter(bin_y);
118  x2 = fXaxis->GetBinCenter(bin_x);
119  y2 = fYaxis->GetBinCenter(bin_y+1);
120  break;
121  case 3:
122  x1 = fXaxis->GetBinCenter(bin_x-1);
123  y1 = fYaxis->GetBinCenter(bin_y-1);
124  x2 = fXaxis->GetBinCenter(bin_x);
125  y2 = fYaxis->GetBinCenter(bin_y);
126  break;
127  case 4:
128  x1 = fXaxis->GetBinCenter(bin_x);
129  y1 = fYaxis->GetBinCenter(bin_y-1);
130  x2 = fXaxis->GetBinCenter(bin_x+1);
131  y2 = fYaxis->GetBinCenter(bin_y);
132  break;
133  }
134  Int_t bin_x1 = RootHelpers::FindBin(fXaxis,x1);
135  if(bin_x1<1) bin_x1=1;
136  Int_t bin_x2 = RootHelpers::FindBin(fXaxis,x2);
137  if(bin_x2>fXaxis->GetNbins()) bin_x2=fXaxis->GetNbins();
138  Int_t bin_y1 = RootHelpers::FindBin(fYaxis,y1);
139  if(bin_y1<1) bin_y1=1;
140  Int_t bin_y2 = RootHelpers::FindBin(fYaxis,y2);
141  if(bin_y2>fYaxis->GetNbins()) bin_y2=fYaxis->GetNbins();
142 
143  Double_t q11;
144  Double_t q12;
145  Double_t q21;
146  Double_t q22;
147  if (otherDimBin > 0)
148  {
149  // X,Y variable and Z fixed
150  if (xAxis == 1 && yAxis == 2)
151  {
152  q11 = histo->GetBinContent(histo->GetBin(bin_x1,bin_y1,otherDimBin));
153  q12 = histo->GetBinContent(histo->GetBin(bin_x1,bin_y2,otherDimBin));
154  q21 = histo->GetBinContent(histo->GetBin(bin_x2,bin_y1,otherDimBin));
155  q22 = histo->GetBinContent(histo->GetBin(bin_x2,bin_y2,otherDimBin));
156  }
157  // X,Z variable and Y fixed
158  else if (xAxis == 1 && yAxis == 3)
159  {
160  q11 = histo->GetBinContent(histo->GetBin(bin_x1,otherDimBin,bin_y1));
161  q12 = histo->GetBinContent(histo->GetBin(bin_x1,otherDimBin,bin_y2));
162  q21 = histo->GetBinContent(histo->GetBin(bin_x2,otherDimBin,bin_y1));
163  q22 = histo->GetBinContent(histo->GetBin(bin_x2,otherDimBin,bin_y2));
164  }
165  // Y,Z variable and X fixed
166  else if (xAxis == 2 && yAxis == 3)
167  {
168  q11 = histo->GetBinContent(histo->GetBin(otherDimBin,bin_x1,bin_y1));
169  q12 = histo->GetBinContent(histo->GetBin(otherDimBin,bin_x1,bin_y2));
170  q21 = histo->GetBinContent(histo->GetBin(otherDimBin,bin_x2,bin_y1));
171  q22 = histo->GetBinContent(histo->GetBin(otherDimBin,bin_x2,bin_y2));
172  }
173  // Y,X variable and Z fixed
174  else if (xAxis == 2 && yAxis == 1)
175  {
176  q11 = histo->GetBinContent(histo->GetBin(bin_y1,bin_x1,otherDimBin));
177  q12 = histo->GetBinContent(histo->GetBin(bin_y1,bin_x2,otherDimBin));
178  q21 = histo->GetBinContent(histo->GetBin(bin_y2,bin_x1,otherDimBin));
179  q22 = histo->GetBinContent(histo->GetBin(bin_y2,bin_x2,otherDimBin));
180  }
181  // Z,X variable and Y fixed
182  else if (xAxis == 3 && yAxis == 1)
183  {
184  q11 = histo->GetBinContent(histo->GetBin(bin_y1,otherDimBin,bin_x1));
185  q12 = histo->GetBinContent(histo->GetBin(bin_y1,otherDimBin,bin_x2));
186  q21 = histo->GetBinContent(histo->GetBin(bin_y2,otherDimBin,bin_x1));
187  q22 = histo->GetBinContent(histo->GetBin(bin_y2,otherDimBin,bin_x2));
188  }
189  // Z,Y variable and X fixed
190  else if (xAxis == 3 && yAxis == 2)
191  {
192  q11 = histo->GetBinContent(histo->GetBin(otherDimBin,bin_y1,bin_x1));
193  q12 = histo->GetBinContent(histo->GetBin(otherDimBin,bin_y1,bin_x2));
194  q21 = histo->GetBinContent(histo->GetBin(otherDimBin,bin_y2,bin_x1));
195  q22 = histo->GetBinContent(histo->GetBin(otherDimBin,bin_y2,bin_x2));
196  }
197  else
198  {
199  histo->Error("Interpolate2D","Unsupported axis combination: (x,y)=(%d,%d) with one bin fixed",xAxis,yAxis);
200  return 0;
201  }
202  }
203  else
204  {
205  // X,Y variable, no Z
206  if (xAxis == 1 && yAxis == 2)
207  {
208  q11 = histo->GetBinContent(histo->GetBin(bin_x1,bin_y1));
209  q12 = histo->GetBinContent(histo->GetBin(bin_x1,bin_y2));
210  q21 = histo->GetBinContent(histo->GetBin(bin_x2,bin_y1));
211  q22 = histo->GetBinContent(histo->GetBin(bin_x2,bin_y2));
212  }
213  // Y,X variable, no Z
214  else if (xAxis == 2 && yAxis == 1)
215  {
216  q11 = histo->GetBinContent(histo->GetBin(bin_y1,bin_x1));
217  q12 = histo->GetBinContent(histo->GetBin(bin_y1,bin_x2));
218  q21 = histo->GetBinContent(histo->GetBin(bin_y2,bin_x1));
219  q22 = histo->GetBinContent(histo->GetBin(bin_y2,bin_x2));
220  }
221  else
222  {
223  histo->Error("Interpolate2D","Unsupported axis combination: (x,y)=(%d,%d)",xAxis,yAxis);
224  return 0;
225  }
226  }
227 
228  Double_t d = 1.0*(x2-x1)*(y2-y1);
229  f = 1.0*q11/d*(x2-x)*(y2-y)+1.0*q21/d*(x-x1)*(y2-y)+1.0*q12/d*(x2-x)*(y-y1)+1.0*q22/d*(x-x1)*(y-y1);
230  return f;
231 }
232 
233 double Interpolate(const TH1* histo, const double x, const double y, const double z)
234 {
235  // Copied from ROOT directly and trivially modified, all credit to ROOT authors of TH1, TH2, and TH3 Interpolate methods
236  // This is done because I want a const version of interpolation, and none of the methods require modification of the histogram
237  // Probable reason is that FindBin isn't const, but there should be a const version...
238  const TAxis* fXaxis = histo->GetXaxis();
239  const TAxis* fYaxis = histo->GetYaxis();
240  const TAxis* fZaxis = histo->GetZaxis();
241 
242  // Find the bin by bin edges
243  Int_t ubx = RootHelpers::FindBin(fXaxis,x);
244  Int_t uby = RootHelpers::FindBin(fYaxis,y);
245  Int_t ubz = RootHelpers::FindBin(fZaxis,z);
246 
247  // Check if the value(s) are outside of the bin range(s)
248  if ( ubx < 1 || ubx > histo->GetNbinsX() || uby < 1 || uby > histo->GetNbinsY() || ubz < 1 || ubz > histo->GetNbinsZ() )
249  {
250  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));
251  return 0;
252  }
253 
254  // Now switch from bin edges to bin centres
255  // Note that we want to support edge cases, so it is possible that ub* == ob*
256  // This functionality is not in original ROOT TH3::Interpolate()
257  // This functionality is inspired by TH2::Interpolate()
258  Int_t obx = ubx + 1;
259  Int_t oby = uby + 1;
260  Int_t obz = ubz + 1;
261 
262  // Calculate distance weights before checking under/overflow bins
263  Double_t xw = fXaxis->GetBinCenter(obx) - fXaxis->GetBinCenter(ubx);
264  Double_t yw = fYaxis->GetBinCenter(oby) - fYaxis->GetBinCenter(uby);
265  Double_t zw = fZaxis->GetBinCenter(obz) - fZaxis->GetBinCenter(ubz);
266 
267  if (x < fXaxis->GetBinCenter(ubx)) { ubx -= 1; obx -= 1; }
268  if (ubx < 1) ubx = 1;
269  if (obx > histo->GetNbinsX()) obx = histo->GetNbinsX();
270 
271  if (y < fYaxis->GetBinCenter(uby)) { uby -= 1; oby -= 1; }
272  if (uby < 1) uby = 1;
273  if (oby > histo->GetNbinsY()) oby = histo->GetNbinsY();
274 
275  if (z < fZaxis->GetBinCenter(ubz)) { ubz -= 1; obz -= 1; }
276  if (ubz < 1) ubz = 1;
277  if (obz > histo->GetNbinsZ()) obz = histo->GetNbinsZ();
278 
279  // Edge cases were tried with weights set including the under/overflow bins (to follow what TH2::Interpolate() does for boundaries)
280  // In some cases, it performed quite poorly
281  // Tests of switching to 2D interpolation with the third dimension fixed appeared to work much better
282  // Thus, the below is now a switch to bilinear interpolation when bin(s) are equal in trilinear interpolation
283 
284  // The below code block contains added functionality with respect to standard ROOT TH3::Interpolate()
285  // Fall back on bilinear interpolation
286  {
287  if (ubz == obz)
288  return RootHelpers::Interpolate2D(histo,x,y,1,2,ubz);
289  else if (uby == oby)
290  return RootHelpers::Interpolate2D(histo,x,z,1,3,uby);
291  else if (ubx == obx)
292  return RootHelpers::Interpolate2D(histo,y,z,2,3,ubx);
293 
294  }
295 
296  // Not a boundary case, resume normal ROOT::TH3::Interpolate()
297  Double_t xd = (x - fXaxis->GetBinCenter(ubx)) / xw;
298  Double_t yd = (y - fYaxis->GetBinCenter(uby)) / yw;
299  Double_t zd = (z - fZaxis->GetBinCenter(ubz)) / zw;
300 
301 
302  Double_t v[] = { histo->GetBinContent( ubx, uby, ubz ), histo->GetBinContent( ubx, uby, obz ),
303  histo->GetBinContent( ubx, oby, ubz ), histo->GetBinContent( ubx, oby, obz ),
304  histo->GetBinContent( obx, uby, ubz ), histo->GetBinContent( obx, uby, obz ),
305  histo->GetBinContent( obx, oby, ubz ), histo->GetBinContent( obx, oby, obz ) };
306 
307 
308  Double_t i1 = v[0] * (1 - zd) + v[1] * zd;
309  Double_t i2 = v[2] * (1 - zd) + v[3] * zd;
310  Double_t j1 = v[4] * (1 - zd) + v[5] * zd;
311  Double_t j2 = v[6] * (1 - zd) + v[7] * zd;
312 
313 
314  Double_t w1 = i1 * (1 - yd) + i2 * yd;
315  Double_t w2 = j1 * (1 - yd) + j2 * yd;
316 
317 
318  Double_t result = w1 * (1 - xd) + w2 * xd;
319 
320  return result;
321 }
322 
323 
324 } // end RootHelpers namespace
plotBeamSpotCompare.x1
x1
Definition: plotBeamSpotCompare.py:216
get_generator_info.result
result
Definition: get_generator_info.py:21
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
hist_file_dump.d
d
Definition: hist_file_dump.py:137
RootHelpers::FindBin
Int_t FindBin(const TAxis *axis, const double x)
Definition: RootHelpers.cxx:14
plotBeamSpotCompare.x2
x2
Definition: plotBeamSpotCompare.py:218
yodamerge_tmp.axis
list axis
Definition: yodamerge_tmp.py:241
bin
Definition: BinsDiffFromStripMedian.h:43
x
#define x
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
doubleTestComp.j1
j1
Definition: doubleTestComp.py:21
z
#define z
RootHelpers.h
makeTRTBarrelCans.y2
tuple y2
Definition: makeTRTBarrelCans.py:18
RootHelpers::Interpolate
double Interpolate(const TH1 *histo, const double x)
Definition: RootHelpers.cxx:39
RootHelpers
Definition: RootHelpers.h:9
hist_file_dump.f
f
Definition: hist_file_dump.py:135
plotBeamSpotCompare.xd
xd
Definition: plotBeamSpotCompare.py:220
plotBeamSpotVxVal.bin
int bin
Definition: plotBeamSpotVxVal.py:83
python.PyAthena.v
v
Definition: PyAthena.py:154
makeTRTBarrelCans.dy
tuple dy
Definition: makeTRTBarrelCans.py:21
y
#define y
RootHelpers::Interpolate2D
double Interpolate2D(const TH1 *histo, const double x, const double y, const int xAxis, const int yAxis, const int otherDimBin)
Definition: RootHelpers.cxx:73
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
plotBeamSpotCompare.histo
histo
Definition: plotBeamSpotCompare.py:415
doubleTestComp.j2
j2
Definition: doubleTestComp.py:22
Generated on Thu Nov 7 2024 21:24:53 for ATLAS Offline Software by   doxygen 1.8.18