GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BLI2D.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 
13 #include <cassert>
14 #include <limits>
15 
18 
19 using namespace genie;
20 
22 
23 //___________________________________________________________________________
25 {
26 
27 }
28 //___________________________________________________________________________
30 {
31  if (fX) { delete [] fX; }
32  if (fY) { delete [] fY; }
33  if (fZ) { delete [] fZ; }
34 }
35 //___________________________________________________________________________
36 int BLI2DGrid::IdxZ(int ix, int iy) const
37 {
38  return ix*fNY+iy;
39 }
40 //___________________________________________________________________________
41 //___________________________________________________________________________
42 //___________________________________________________________________________
44 
46 {
47  this->Init();
48 }
49 //___________________________________________________________________________
51  int nx, double xmin, double xmax, int ny, double ymin, double ymax)
52 {
53  this->Init(nx, xmin, xmax, ny, ymin, ymax);
54 }
55 //___________________________________________________________________________
57  int nx, int ny, double *x, double *y, double *z)
58 {
59  double xmin = x[0];
60  double xmax = x[nx-1];
61  double ymin = y[0];
62  double ymax = y[nx-1];
63 
64  this->Init(nx, xmin, xmax, ny, ymin, ymax);
65 
66  for(int ix=0; ix<nx; ix++) {
67  for(int iy=0; iy<ny; iy++) {
68  this->AddPoint(x[ix], y[iy], z[this->IdxZ(ix,iy)]);
69  }
70  }
71 }
72 //___________________________________________________________________________
73 bool BLI2DUnifGrid::AddPoint(double x, double y, double z)
74 {
75  int ix = TMath::FloorNint( (x - fXmin + fDX/2) / fDX );
76  int iy = TMath::FloorNint( (y - fYmin + fDY/2) / fDY );
77  int iz = this->IdxZ(ix,iy);
78 
79  fZ[iz] = z;
80 
81  fZmin = TMath::Min(z, fZmin);
82  fZmax = TMath::Max(z, fZmax);
83 
84  LOG("BLI2DUnifGrid", pDEBUG)
85  << "Added x = " << x << " (ix = " << ix << ")"
86  << " y = " << y << " (iy = " << iy << ") -> "
87  << " z = " << z << " (iz = " << iz << ")";
88 
89  return true;
90 }
91 //___________________________________________________________________________
92 double BLI2DUnifGrid::Evaluate(double x, double y) const
93 {
94  if(x < fXmin || x > fXmax) return 0.;
95  if(y < fYmin || y > fYmax) return 0.;
96 
97  int ix_lo = TMath::FloorNint( (x - fXmin) / fDX );
98  int iy_lo = TMath::FloorNint( (y - fYmin) / fDY );
99  int ix_hi = ix_lo + 1;
100  int iy_hi = iy_lo + 1;
101 
102  double x1 = fX[ix_lo];
103  double x2 = fX[ix_hi];
104  double y1 = fY[iy_lo];
105  double y2 = fY[iy_hi];
106 
107  double z11 = fZ[ this->IdxZ(ix_lo,iy_lo) ];
108  double z21 = fZ[ this->IdxZ(ix_hi,iy_lo) ];
109  double z12 = fZ[ this->IdxZ(ix_lo,iy_hi) ];
110  double z22 = fZ[ this->IdxZ(ix_hi,iy_hi) ];
111 
112  double z1 = z11 * (x2-x)/(x2-x1) + z21 * (x-x1)/(x2-x1);
113  double z2 = z12 * (x2-x)/(x2-x1) + z22 * (x-x1)/(x2-x1);
114  double z = z1 * (y2-y)/(y2-y1) + z2 * (y-y1)/(y2-y1);
115 
116 /*
117  LOG("BLI2DUnifGrid", pDEBUG) << "x = " << x << " -> nearby nodes: " << x1 << ", " << x2;
118  LOG("BLI2DUnifGrid", pDEBUG) << "y = " << y << " -> nearby nodes: " << y1 << ", " << y2;
119  LOG("BLI2DUnifGrid", pDEBUG) << "z11 := z(" << this->IdxZ(ix_lo,iy_lo) << ") = " << z11;
120  LOG("BLI2DUnifGrid", pDEBUG) << "z21 := z(" << this->IdxZ(ix_hi,iy_lo) << ") = " << z21;
121  LOG("BLI2DUnifGrid", pDEBUG) << "z12 := z(" << this->IdxZ(ix_lo,iy_hi) << ") = " << z12;
122  LOG("BLI2DUnifGrid", pDEBUG) << "z22 := z(" << this->IdxZ(ix_hi,iy_hi) << ") = " << z22;
123  LOG("BLI2DUnifGrid", pDEBUG) << "z1 = " << z1 << ", z2 = " << z2;
124  LOG("BLI2DUnifGrid", pDEBUG) << "interpolated z(x,y) = " << z;
125 */
126 
127  return z;
128 }
129 //___________________________________________________________________________
131  int nx, double xmin, double xmax, int ny, double ymin, double ymax)
132 {
133  fNX = 0;
134  fNY = 0;
135  fNZ = 0;
136  fXmin = 0.;
137  fXmax = 0.;
138  fYmin = 0.;
139  fYmax = 0.;
140  fZmin = std::numeric_limits<double>::max();
141  fZmax = std::numeric_limits<double>::min();
142  fDX = 0.;
143  fDY = 0.;
144  fX = 0;
145  fY = 0;
146  fZ = 0;
147 
148  if(nx>1 && ny>1) {
149  fNX = nx;
150  fNY = ny;
151  fNZ = nx * ny;
152 
153  fXmin = xmin;
154  fXmax = xmax;
155  fYmin = ymin;
156  fYmax = ymax;
157  fZmin = std::numeric_limits<double>::max();
158  fZmax = std::numeric_limits<double>::min();
159 
160  fDX = (xmax-xmin)/(nx-1);
161  fDY = (ymax-ymin)/(ny-1);
162 
163  fX = new double[fNX];
164  fY = new double[fNY];
165  fZ = new double[fNZ];
166 
167  for(int i=0; i<fNX; i++) { fX[i] = xmin + i*fDX; }
168  for(int i=0; i<fNY; i++) { fY[i] = ymin + i*fDY; }
169  for(int i=0; i<fNZ; i++) { fZ[i] = 0.; }
170  }
171 }
172 //___________________________________________________________________________
173 //___________________________________________________________________________
174 //___________________________________________________________________________
176 
178 {
179  this->Init();
180 }
181 //___________________________________________________________________________
183  int nx, double xmin, double xmax, int ny, double ymin, double ymax)
184 {
185  this->Init(nx, xmin, xmax, ny, ymin, ymax);
186 }
187 //___________________________________________________________________________
189  int nx, int ny, double *x, double *y, double *z)
190 {
191  double xmin = x[0];
192  double xmax = x[nx-1];
193  double ymin = y[0];
194  double ymax = y[ny-1];
195 
196  this->Init(nx, xmin, xmax, ny, ymin, ymax);
197 
198  for(int ix=0; ix<nx; ix++) {
199  for(int iy=0; iy<ny; iy++) {
200  this->AddPoint(x[ix], y[iy], z[this->IdxZ(ix,iy)]);
201  }
202  }
203 }
204 //___________________________________________________________________________
205 bool BLI2DNonUnifGrid::AddPoint(double x, double y, double z)
206 {
207 
208  // check the x,y values' existence before moving anything
209  // if they do, use them
210  // if they don't, add them in the appropriate place
211  //
212  int xidx = -1;
213  int yidx = -1;
214  bool changex = true;
215  bool changey = true;
216  // first, check and see if the x,y values exist
217  // NOTE: == should not be used with double values:
218  // instead, a tolerance of ~5 significant digits is checked
219  for(int i=0;i<fNFillX;i++)
220  {
221  if (!(TMath::Abs(x-fX[i])<=.5*.0000001*(TMath::Abs(x)+TMath::Abs(fX[i]))) && x<fX[i])
222  {
223  // missed one, make sure there aren't too many x and move everything up
224  if (fNFillX==fNX)
225  {
226  LOG("BLI2DNonUnifGrid", pWARN) << "Warning: too many x values, cannot add x = "<<x;
227  return false;
228  }
229  else
230  {
231  xidx=i;
232  changex=true;
233  break;
234  }
235  }
236  else
237  {
238  if (TMath::Abs(x-fX[i])<=.5*.0000001*(TMath::Abs(x)+TMath::Abs(fX[i]))) {
239  xidx=i;
240  LOG("BLI2DNonUnifGrid", pDEBUG) << "x value found at index "<<i;
241  changex=false; break;}
242  changex=true;
243  }
244  }
245  if (changex && xidx<0) xidx=fNFillX;
246  if (changex) LOG("BLI2DNonUnifGrid", pDEBUG) << "Adding x value at index "<<xidx;
247 
248  for(int i=0;i<fNFillY;i++)
249  {
250  if (!(TMath::Abs(y-fY[i])<=.5*.0000001*(TMath::Abs(y)+TMath::Abs(fY[i]))) && y<fY[i])
251  {
252  if (fNFillY==fNY)
253  {
254  LOG("BLI2DNonUnifGrid", pWARN) << "Warning: too many y values, cannot add y = "<<y;
255  return false;
256  }
257  else
258  {
259  yidx=i;
260  changey=true;
261  break;
262  }
263  }
264  else
265  {
266  if (TMath::Abs(y-fY[i])<=.5*.0000001*(TMath::Abs(y)+TMath::Abs(fY[i]))) {
267  yidx=i;
268  LOG("BLI2DNonUnifGrid", pDEBUG) << "y value found at index "<<i;
269  changey=false; break;}
270  changey=true;
271  }
272  }
273  if (changey && yidx<0) yidx=fNFillY;
274  if (changey) LOG("BLI2DNonUnifGrid", pDEBUG) << "Adding y value at index "<<yidx;
275 
276  // make new entries if needed
277  if (changex && xidx>=0)
278  {
279  for (int j=0;j<fNFillX-xidx;j++)
280  {
281  fX[fNFillX-j]=fX[fNFillX-j-1];
282  for (int k=0;k<fNFillY;k++)
283  {
284  fZ[ this->IdxZ(fNFillX-j,k) ]
285  = fZ[ this->IdxZ(fNFillX-j-1,k) ];
286  }
287  }
288  fX[xidx]=x;
289  fNFillX++;
290  fXmin = TMath::Min(x,fXmin);
291  fXmax = TMath::Max(x,fXmax);
292  }
293  if (changey && yidx>=0)
294  {
295  for (int j=0;j<fNFillY-yidx;j++)
296  {
297  fY[fNFillY-j]=fY[fNFillY-j-1];
298  for (int k=0;k<fNFillX;k++)
299  {
300  fZ[ this->IdxZ(k,fNFillY-j) ]
301  = fZ[ this->IdxZ(k,fNFillY-j-1) ];
302  }
303  }
304  fY[yidx]=y;
305  fNFillY++;
306  fYmin = TMath::Min(y,fYmin);
307  fYmax = TMath::Max(y,fYmax);
308  }
309 
310  int iz = this->IdxZ(xidx,yidx);
311 
312  fZ[iz] = z;
313 
314  fZmin = TMath::Min(z, fZmin);
315  fZmax = TMath::Max(z, fZmax);
316 
317  LOG("BLI2DNonUnifGrid", pDEBUG)
318  << "Added x = " << x << " (ix = " << xidx << ")"
319  << " y = " << y << " (iy = " << yidx << ") -> "
320  << " z = " << z << " (iz = " << iz << ")";
321 
322  return true;
323 }
324 //___________________________________________________________________________
325 double BLI2DNonUnifGrid::Evaluate(double x, double y) const
326 {
327  double evalx=TMath::Min(x,fXmax);
328  evalx=TMath::Max(evalx,fXmin);
329  double evaly=TMath::Min(y,fYmax);
330  evaly=TMath::Max(evaly,fYmin);
331 
332  int ix_lo = -2;
333  int iy_lo = -2;
334  for (int i=0;i<fNFillX;i++)
335  {
336  if (evalx<=fX[fNFillX-1-i]) ix_lo=fNFillX-2-i;
337  else break;
338  }
339  for (int i=0;i<fNFillY;i++)
340  {
341  if (evaly<=fY[fNFillY-1-i]) iy_lo=fNFillY-2-i;
342  else break;
343  }
344  int ix_hi = ix_lo + 1;
345  int iy_hi = iy_lo + 1;
346 
347  // in case x = xmin
348  if (ix_lo==-1) {ix_lo++; ix_hi++;}
349  if (iy_lo==-1) {iy_lo++; iy_hi++;}
350  // in case x = xmax
351  if (ix_lo==-2) {ix_lo=fNFillX-2; ix_hi=fNFillX-1;}
352  if (iy_lo==-2) {iy_lo=fNFillY-2; iy_hi=fNFillY-1;}
353  // if an error occurs
354  if (ix_lo<0 || iy_lo<0 ) return 0.;
355  if (ix_hi>fNFillX || iy_hi>fNFillY) return 0.;
356 
357  double x1 = fX[ix_lo];
358  double x2 = fX[ix_hi];
359  double y1 = fY[iy_lo];
360  double y2 = fY[iy_hi];
361 
362  double z11 = fZ[ this->IdxZ(ix_lo,iy_lo) ];
363  double z21 = fZ[ this->IdxZ(ix_hi,iy_lo) ];
364  double z12 = fZ[ this->IdxZ(ix_lo,iy_hi) ];
365  double z22 = fZ[ this->IdxZ(ix_hi,iy_hi) ];
366 
367  double z1 = z11 * (x2-evalx)/(x2-x1) + z21 * (evalx-x1)/(x2-x1);
368  double z2 = z12 * (x2-evalx)/(x2-x1) + z22 * (evalx-x1)/(x2-x1);
369  double z = z1 * (y2-evaly)/(y2-y1) + z2 * (evaly-y1)/(y2-y1);
370 
371  /*
372  LOG("BLI2DNonUnifGrid", pINFO) << "x = " << evalx << " -> nearby nodes: " << x1 << ", " << x2;
373  LOG("BLI2DNonUnifGrid", pINFO) << "y = " << evaly << " -> nearby nodes: " << y1 << ", " << y2;
374  LOG("BLI2DNonUnifGrid", pINFO) << "xmin = " << fXmin << ", xmax = " << fXmax;
375  LOG("BLI2DNonUnifGrid", pINFO) << "ymin = " << fYmin << ", ymax = " << fYmax;
376  LOG("BLI2DNonUnifGrid", pINFO) << "z11 := z(" << this->IdxZ(ix_lo,iy_lo) << ") = " << z11;
377  LOG("BLI2DNonUnifGrid", pINFO) << "z21 := z(" << this->IdxZ(ix_hi,iy_lo) << ") = " << z21;
378  LOG("BLI2DNonUnifGrid", pINFO) << "z12 := z(" << this->IdxZ(ix_lo,iy_hi) << ") = " << z12;
379  LOG("BLI2DNonUnifGrid", pINFO) << "z22 := z(" << this->IdxZ(ix_hi,iy_hi) << ") = " << z22;
380  LOG("BLI2DNonUnifGrid", pINFO) << "z1 = " << z1 << ", z2 = " << z2;
381  LOG("BLI2DNonUnifGrid", pINFO) << "interpolated z(x,y) = " << z;
382  */
383 
384  return z;
385 }
386 //___________________________________________________________________________
388  int nx, double xmin, double xmax, int ny, double ymin, double ymax)
389 {
390  fNX = 0;
391  fNY = 0;
392  fNZ = 0;
393  fNFillX= 0;
394  fNFillY= 0;
395  fXmin = 0.;
396  fXmax = 0.;
397  fYmin = 0.;
398  fYmax = 0.;
399  fZmin = std::numeric_limits<double>::max();
400  fZmax = std::numeric_limits<double>::min();
401  fX = 0;
402  fY = 0;
403  fZ = 0;
404 
405  if(nx>1 && ny>1) {
406  fNX = nx;
407  fNY = ny;
408  fNZ = nx * ny;
409 
410  fXmin = xmin;
411  fXmax = xmax;
412  fYmin = ymin;
413  fYmax = ymax;
414  fZmin = std::numeric_limits<double>::max();
415  fZmax = std::numeric_limits<double>::min();
416 
417  fX = new double[fNX];
418  fY = new double[fNY];
419  fZ = new double[fNZ];
420 
421  for(int i=0; i<fNX; i++) { fX[i] = 0.; }
422  for(int i=0; i<fNY; i++) { fY[i] = 0.; }
423  for(int i=0; i<fNZ; i++) { fZ[i] = 0.; }
424  }
425 }
426 //___________________________________________________________________________
double fXmin
Definition: BLI2D.h:62
double fYmax
Definition: BLI2D.h:65
double * fY
Definition: BLI2D.h:58
Bilinear interpolation of 2D functions on a regular grid.
Definition: BLI2D.h:75
bool AddPoint(double x, double y, double z)
Definition: BLI2D.cxx:205
double fZmax
Definition: BLI2D.h:67
double * fX
Definition: BLI2D.h:57
double Evaluate(double x, double y) const
Definition: BLI2D.cxx:325
double fDX
Definition: BLI2D.h:60
void Init(int nx=0, double xmin=0, double xmax=0, int ny=0, double ymin=0, double ymax=0)
Definition: BLI2D.cxx:130
bool AddPoint(double x, double y, double z)
Definition: BLI2D.cxx:73
double fDY
Definition: BLI2D.h:61
void Init(int nx=0, double xmin=0, double xmax=0, int ny=0, double ymin=0, double ymax=0)
Definition: BLI2D.cxx:387
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void Init(void)
Definition: gXSecComp.cxx:138
#define pWARN
Definition: Messenger.h:60
double Evaluate(double x, double y) const
Definition: BLI2D.cxx:92
int IdxZ(int ix, int iy) const
Definition: BLI2D.cxx:36
ClassImp(CacheBranchFx)
double fYmin
Definition: BLI2D.h:64
virtual ~BLI2DGrid()
Definition: BLI2D.cxx:29
double fXmax
Definition: BLI2D.h:63
double fZmin
Definition: BLI2D.h:66
double * fZ
Definition: BLI2D.h:59
#define pDEBUG
Definition: Messenger.h:63