LASPeakFinder.cc

Go to the documentation of this file.

#include "Alignment/LaserAlignment/interface/LASPeakFinder.h"

LASPeakFinder::LASPeakFinder() {
  amplitudeThreshold = 0.;  // default
}

bool LASPeakFinder::FindPeakIn(const LASModuleProfile& aProfile,
                               std::pair<double, double>& result,
                               TH1D* histogram,
                               const double offset) {
  // this function will be propagated to the histo output file,
  // therefore we do some cosmetics already here
  TF1 fitFunction("fitFunction", "gaus");
  fitFunction.SetLineColor(2);
  fitFunction.SetLineWidth(2);

  std::pair<int, double> largestAmplitude(0, 0.);  // strip, amplitude
  double anAmplitude = 0.;

  // need approximate position -> cast to int
  const int approxOffset = static_cast<int>(offset);

  // expected beam position (in strips)
  const unsigned int meanPosition = 256 + approxOffset;

  // backplane "alignment hole" approx. half size (in strips)
  const unsigned int halfWindowSize = 33;

  // loop over the strips in the "alignment hole"
  // to fill the histogram
  // and determine fit parameter estimates
  double sum0 = 0.;  // this is the sum of all amplitudes
  for (unsigned int strip = meanPosition - halfWindowSize; strip < meanPosition + halfWindowSize; ++strip) {
    anAmplitude = aProfile.GetValue(strip);
    sum0 += anAmplitude;
    histogram->SetBinContent(1 + strip, anAmplitude);
    if (anAmplitude > largestAmplitude.second) {
      largestAmplitude.first = strip;
      largestAmplitude.second = anAmplitude;
    }
  }

  // loop outside the "alignment hole"
  // to determine the noise level = sqrt(variance)
  double sum1 = 0., sum2 = 0.;
  int nStrips = 0;
  for (unsigned int strip = 0; strip < 512; ++strip) {
    if (strip < meanPosition - halfWindowSize || strip > meanPosition + halfWindowSize) {
      anAmplitude = aProfile.GetValue(strip);
      sum0 += anAmplitude;
      sum1 += anAmplitude;
      sum2 += pow(anAmplitude, 2);
      nStrips++;
    }
  }

  // noise as sqrt of the amplitude variance
  const double noise = sqrt(1. / (nStrips - 1) * (sum2 - pow(sum1, 2) / nStrips));

  // empty profile?
  if (fabs(sum0) < 1.e-3) {
    std::cerr << " [LASPeakFinder::FindPeakIn] ** WARNING: Empty profile (sum=" << sum0 << ")." << std::endl;
    return false;
  }

  // no reasonable peak?
  if (largestAmplitude.second < amplitudeThreshold * noise) {
    std::cerr << " [LASPeakFinder::FindPeakIn] ** WARNING: No reasonably large peak." << std::endl;
    return false;
  }

  // prepare fit function: starting values..
  fitFunction.SetParameter(0, largestAmplitude.second);  // amp
  fitFunction.SetParameter(1, largestAmplitude.first);   // mean
  fitFunction.SetParameter(2, 3.);                       // width

  // ..and parameter limits
  fitFunction.SetParLimits(
      0, largestAmplitude.second * 0.3, largestAmplitude.second * 1.8);  // amp of the order of the peak height
  fitFunction.SetParLimits(
      1, largestAmplitude.first - 12, largestAmplitude.first + 12);  // mean around the peak maximum
  fitFunction.SetParLimits(2, 0.5, 8.);                              // reasonable width

  // fit options: Quiet / all Weights=1 / better Errors / enable fix&Bounds on parameters
  histogram->Fit(&fitFunction, "QWEB", "", largestAmplitude.first - 12, largestAmplitude.first + 12);

  // prepare output
  result.first = fitFunction.GetParameter(1);
  result.second = fitFunction.GetParError(1);

  return true;
}

void LASPeakFinder::SetAmplitudeThreshold(double aThreshold) { amplitudeThreshold = aThreshold; }