BeamFit2D.cpp

Go to the documentation of this file.

#include "BeamFit2D.hpp"
 
Bool_t BeamFitUtils::InEllipseXY(const BeamFit2D &b, Double_t x, Double_t y,
                                 Double_t nx, Double_t ny) {
  Double_t dx = x - b.mu_x;
  Double_t dy = y - b.mu_y;
  Double_t sx = b.sigma_x, sy = b.sigma_y, rho = b.rho;
  if (sx <= 0 || sy <= 0)
    return kFALSE;
  // Correlated 2D Gaussian ellipse (Mahalanobis distance):
  //   χ² = [(dx/σx)² + (dy/σy)² - 2ρ·dx·dy/(σx·σy)] / (1-ρ²)
  // Inside n-sigma: χ² < n²  (nx serves as the n-sigma level).
  Double_t dx_s = dx / sx, dy_s = dy / sy;
  Double_t r2 = rho * rho;
  Double_t chi2 =
      (dx_s * dx_s + dy_s * dy_s - 2.0 * rho * dx_s * dy_s) / (1.0 - r2);
  return chi2 < nx * nx;
}
 
Moments2D BeamFitUtils::ComputeMoments(TH2F *h, Int_t lo_bx, Int_t hi_bx,
                                       Int_t lo_by, Int_t hi_by,
                                       Double_t thresh, Double_t bw_x,
                                       Double_t bw_y) {
  Moments2D m;
  Double_t W = 0, Mx = 0, My = 0, Cxx = 0, Cyy = 0, Cxy = 0;
  for (Int_t ix = lo_bx; ix <= hi_bx; ix++) {
    Double_t x = h->GetXaxis()->GetBinCenter(ix);
    for (Int_t iy = lo_by; iy <= hi_by; iy++) {
      Double_t w = h->GetBinContent(ix, iy);
      if (w < thresh)
        continue;
      Double_t y = h->GetYaxis()->GetBinCenter(iy);
      W += w;
      Mx += w * x;
      My += w * y;
    }
  }
  if (W <= 0)
    return m;
  Mx /= W;
  My /= W;
  for (Int_t ix = lo_bx; ix <= hi_bx; ix++) {
    Double_t x = h->GetXaxis()->GetBinCenter(ix);
    for (Int_t iy = lo_by; iy <= hi_by; iy++) {
      Double_t w = h->GetBinContent(ix, iy);
      if (w < thresh)
        continue;
      Double_t y = h->GetYaxis()->GetBinCenter(iy);
      Double_t dx = x - Mx, dy = y - My;
      Cxx += w * dx * dx;
      Cyy += w * dy * dy;
      Cxy += w * dx * dy;
    }
  }
  Cxx /= W;
  Cyy /= W;
  Cxy /= W;
  m.mu_x = Mx;
  m.mu_y = My;
  m.sigma_x = std::max(std::sqrt(std::max(Cxx, 0.0)), 2.0 * bw_x);
  m.sigma_y = std::max(std::sqrt(std::max(Cyy, 0.0)), 2.0 * bw_y);
  Double_t r = Cxy / (m.sigma_x * m.sigma_y);
  if (r > 0.95)
    r = 0.95;
  if (r < -0.95)
    r = -0.95;
  m.rho = r;
  m.weight = W;
  return m;
}