ATLAS Offline Software
Loading...
Searching...
No Matches
Functions
JetHelpers Namespace Reference

Functions

double Interpolate (const TH1 *histo, const double x)
double Interpolate (const TH1 *histo, const double x, const double y)
double Interpolate (const TH1 *histo, const double x, const double y, const double z)
double Interpolate2D (const TH1 *histo, const double x, const double y, const int xAxis=1, const int yAxis=2, const int otherDimBin=-1)

Function Documentation

◆ Interpolate() [1/3]

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

Definition at line 16 of file JetHelpers.cxx.

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}
x
#define x

◆ Interpolate() [2/3]

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

Definition at line 44 of file JetHelpers.cxx.

45{
46 // Call the unified method for consistency
47 return Interpolate2D(histo,x,y);
48}
y
#define y
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

◆ Interpolate() [3/3]

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

Definition at line 225 of file JetHelpers.cxx.

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}
z
#define z
get_generator_info.result
result
Definition get_generator_info.py:21

◆ Interpolate2D()

double JetHelpers::Interpolate2D ( const TH1 * histo,
const double x,
const double y,
const int xAxis = 1,
const int yAxis = 2,
const int otherDimBin = -1 )

Definition at line 50 of file JetHelpers.cxx.

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}
hist_file_dump.d
d
Definition hist_file_dump.py:142
hist_file_dump.f
f
Definition hist_file_dump.py:140
Generated on for ATLAS Offline Software by doxygen 1.14.0