CaloClusterCorr Namespace Reference

Classes
class	CaloSwCalibHitsShowerDepth
	Helper to calculate the shower depth as used in the calib hits SW cluster correction. More...
class	DDHelper
struct	et_compare_larem_only
	Helper to compare two cells by Et. More...
class	SamplingHelper
	Sampling calculator helper class. More...
class	SamplingHelper_CaloCellList
	Version of helper class for cells taken from StoreGate. More...
class	SamplingHelper_Cluster
	Version of helper class for cells taken from the cluster itself. More...
class	Segmentation

Functions
float	interpolate (const CaloRec::Array< 2 > &a, float x, unsigned int degree, unsigned int ycol=1, const CaloRec::Array< 1 > &regions=CaloRec::Array< 1 >(), int n_points=-1, bool fixZero=false)
	Polynomial interpolation in a table. More...
void	etaphi_range (const CaloDetDescrManager &dd_man, double eta, double phi, CaloCell_ID::CaloSample sampling, double &deta, double &dphi)
	Return eta/phi ranges encompassing +- 1 cell. More...

Function Documentation

etaphi_range()

void CaloClusterCorr::etaphi_range	(	const CaloDetDescrManager &	dd_man,
		double	eta,
		double	phi,
		CaloCell_ID::CaloSample	sampling,
		double &	deta,
		double &	dphi
	)

Return eta/phi ranges encompassing +- 1 cell.

Parameters

	eta	Central eta value.
	phi	Central phi value.
	sampling	The sampling to use.
[out]	deta	Range in eta.
[out]	dphi	Range in phi.

This can be a little tricky due to misalignments and the fact that cells have different sizes in different regions. Also, CaloLayerCalculator takes only a symmetric eta range. We try to find the neighboring cells by starting from the center cell and looking a little bit more than half its width in either direction, and finding the centers of those cells. Then we use the larger of these for the symmetric range.

Definition at line 64 of file CaloFillRectangularCluster.cxx.

{
  deta = 0;
  dphi = 0;
 
  // Get the DD element for the central cell.
  const CaloDetDescrElement* elt = dd_man.get_element_raw (sampling, eta, phi);
  if (!elt) return;
 
  // Should be smaller than the eta half-width of any cell.
  const double eps = 0.001;
 
  // Now look in the negative eta direction.
  const CaloDetDescrElement* elt_l = dd_man.get_element_raw
    (sampling,
     eta - elt->deta() - eps,
     phi);
  double deta_l = 0; // Eta difference on the low (left) side.
  if (elt_l)
    deta_l = std::abs (eta - elt_l->eta_raw()) + eps;
 
  // Now look in the positive eta direction.
  const CaloDetDescrElement* elt_r = dd_man.get_element_raw
    (sampling,
     eta + elt->deta() + eps,
     phi);
  double deta_r = 0; // Eta difference on the high (right) side.
  if (elt_r)
    deta_r = std::abs (eta - elt_r->eta_raw()) + eps;
 
  // Total deta is twice the maximum.
  deta = 2 * std::max (deta_r, deta_l);
 
  // Now for the phi variation.
  // The phi size can change as a function of eta, but not of phi.
  // Thus we have to look again at the adjacent eta cells, and
  // take the largest variation.
 
  // Now look in the negative eta direction.
  elt_l = dd_man.get_element_raw
    (sampling,
     eta  - elt->deta() - eps,
     CaloPhiRange::fix (phi - elt->dphi() - eps));
  double dphi_l = 0; // Phi difference on the low-eta () side.
  if (elt_l)
    dphi_l = std::abs (CaloPhiRange::fix (phi - elt_l->phi_raw())) + eps;
 
  // Now look in the positive eta direction.
  elt_r = dd_man.get_element_raw
    (sampling,
     eta + elt->deta() + eps,
     CaloPhiRange::fix (phi - elt->dphi() - eps));
  double dphi_r = 0; // Phi difference on the positive (down) side.
  if (elt_r)
    dphi_r = std::abs (CaloPhiRange::fix (phi - elt_r->phi_raw())) + eps;
 
  // Total dphi is twice the maximum.
  dphi = 2 * std::max (dphi_l, dphi_r);
}

interpolate()

float CaloClusterCorr::interpolate	(	const CaloRec::Array< 2 > &	a,
		float	x,
		unsigned int	degree,
		unsigned int	ycol = 1,
		const CaloRec::Array< 1 > &	regions = CaloRec::Array<1>(),
		int	n_points = -1,
		bool	fixZero = false
	)

Polynomial interpolation in a table.

Parameters

a	Interpolation table. The x values are in the first column (column 0), and by default, the y values are in the second column (column 1). The x values must be in ascending order, with no duplicates.
x	The value to interpolate.
degree	The degree of the interpolating polynomial.
ycol	The column number of the y values. (0-based.)
regions	Sometimes, you want to divide the range being interpolated into several distinct regions, where the interpolated function may be discontinuous at the boundaries between regions. To do this, supply this argument, which should contain a list (in ascending order) of the x-coordinates of the region boundaries. When the interpolation runs against a boundary, the algorithm will add a copy of the last point before the boundary, positioned at the boundary. This helps to control runaway extrapolation leading up to the boundary.
n_points	The number of interpolation points in the table to use. If -1, then the entire table is used.
fixZero	If true, remove zeros among y-values by averaging adjacent points.

Returns

The interpolated value.

The method used is Newtonian interpolation. Based on the cernlib routine divdif.

Parameters

a	Interpolation table. The x values are in the first column (column 0), and by default, the y values are in the second column (column 1).
x	The value to interpolate.
degree	The degree of the interpolating polynomial.
ycol	The column number of the y values. (0-based.)
regions	Sometimes, you want to divide the range being interpolated into several distinct regions, where the interpolated function may be discontinuous at the boundaries between regions. To do this, supply this argument, which should contain a list (in ascending order) of the x-coordinates of the region boundaries. When the interpolation runs against a boundary, the algorithm will add a copy of the last point before the boundary, positioned at the boundary. This helps to control runaway extrapolation leading up to the boundary.
n_points	The number of interpolation points in the table to use. If -1, then the entire table is used.
fixZero	If true, remove zeros among y-values by averaging adjacent points.

Returns

The interpolated value.

The method used is Newtonian interpolation. Based on the cernlib routine divdif.

Definition at line 75 of file interpolate.cxx.

{
  const int xcol = 0;
 
  // Check arguments.
  if (n_points < 0)
    n_points = static_cast<int> (a.size());
  if (n_points < 2 || degree < 1) {
    std::abort();
  }
  degree = std::min (degree, static_cast<unsigned int> (n_points) - 1);
 
  // Find subscripts of the input value in the input arrays.
  unsigned int ix = std::lower_bound (a.begin(), a.begin()+n_points, x,
                                      xcompare()) -
                      a.begin();
  unsigned int ir = std::lower_bound (regions.begin(), regions.end(), x) -
                      regions.begin();
 
  // Number of points to try for.
  // Either degree+1 or degree+2, whichever is even,
  // to give the same number of points on each side.
  // If we run up against an edge or a boundary, we'll
  // fall back to using just degree+1 points (or fewer if we can't
  // even get that many).
  // If we end up using degree+2 points, we'll do two interpolations
  // of degree degree and average them.
  unsigned int npts = degree + 2 - (degree%2);
 
  // If we run up against the edge of a region boundary,
  // we'll want to add a pseudopoint right at the boundary
  // (copying the point closest to the boundary) instead of the
  // point farthest away from it.
  bool extralo = false;
  bool extrahi = false;
 
  // Starting point index, not considering edges or boundaries.
  int ilo = ix - npts/2;
 
  // Make sure this point is within the array range and has not
  // crossed a region boundary.
  if (ilo < 0) {
    ilo = 0;
    npts = degree+1;
  }
  while (ilo < n_points &&
         ir > 0 && a[ilo][xcol] < regions[ir-1])
  {
    ++ilo;
    npts = degree+1;
    extralo = true;
  }
 
  // Same deal for the right hand edge.
  // ihi is one past the last point to use.
  bool himoved = false;
  int ihi = ilo + npts;
  if (ihi > n_points) {
    ihi = n_points;
    npts = degree+1;
    himoved = true;
  }
  while (ihi > 0 && ir < regions.size() && a[ihi-1][xcol] >= regions[ir]) {
    --ihi;
    npts = degree+1;
    himoved = true;
    extrahi = true;
  }
 
  // The adjustment due to ihi may have moved the low point of the range
  // back over a boundary.  If so, we lose points from the fit.
  bool lomoved = false;
  ilo = ihi - npts;
  if (ilo < 0) {
    ilo = 0;
    lomoved = true;
  }
  while (ilo < n_points &&
         ir > 0 && a[ilo][xcol] < regions[ir-1])
  {
    ++ilo;
    extralo = true;
    lomoved = true;
  }
 
  // Copy the points we're going to use to arrays t and d.
  // t gets the x coordinates, d gets the y coordinates.
  // Note that the order doesn't matter.
  // Reserve two extra points in case we end up adding some.
  npts = ihi - ilo;
  std::vector<float> t;
  t.reserve (npts+2);
  std::vector<float> d;
  d.reserve (npts+2);
 
  // Add the pseudopoints, if needed, removing a point from the other end.
  // Be careful not to duplicate a point if there's already one
  // right at the boundary.
  if (extralo && is_different (a[ilo][xcol], regions[ir-1])) {
    if (!himoved)
      --ihi;
    else
      ++npts;
    t.push_back (regions[ir-1]);
    d.push_back (a[ilo][ycol]);
  }
  if (extrahi && is_different (a[ihi-1][xcol], regions[ir])) {
    if (!lomoved)
      ++ilo;
    else
      ++npts;
    t.push_back (regions[ir]);
    d.push_back (a[ihi-1][ycol]);
  }
 
  // Add the rest if the points.
  for (; ilo < ihi; ++ilo) {
    t.push_back (a[ilo][xcol]);
    d.push_back (a[ilo][ycol]);
  }
 
  // Now figure out the interpolation degree we're really going to use.
  assert (t.size() == npts);
  assert (d.size() == npts);
  degree = std::min (degree, npts-1);
 
  // Option to remove zeros in the interpolation table, by averaging
  // adjacent points.  Used to handle a couple bad points in the rfac-v5
  // table.
  if (fixZero) {
    for (size_t i = 1; i < npts-1; i++) {
      if (d[i] == 0) {
        d[i] = (d[i-1] + d[i+1])/2;
      }
    }
  }
 
  // True if we're averaging together two interpolations
  // (to get a symmetric range).
  bool extra = npts != degree+1;
 
  // If we're averaging two interpolations, we need to be sure the
  // two extreme points are at the end of the table.
  // (If extra is true, extrahi and extralo must be false.)
  if (extra) {
    std::swap (t[0], t[npts-2]);
    std::swap (d[0], d[npts-2]);
  }
 
  // Here we start the actual interpolation.
  // Replace d by the leading diagonal of a divided-difference table,
  // supplemented by an extra line if extra is true.
  for (unsigned int l=0; l<degree; l++) {
    if (extra)
      d[degree+1] = (d[degree+1]-d[degree-1])/(t[degree+1]-t[degree-1-l]);
    for (unsigned int i = degree; i > l; --i)
      d[i] = (d[i]-d[i-1])/(t[i]-t[i-1-l]);
  }
 
  // Evaluate the Newton interpolation formula at x, averaging
  // two values of the last difference if extra is true.
  float sum = d[degree];
  if (extra)
    sum = 0.5*(sum+d[degree+1]);
  for (int j = degree-1; j >= 0; --j)
    sum = d[j] + (x - t[j]) * sum;
 
  return sum;
}