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