File indexing completed on 2024-04-06 11:56:26
0001
0002 #include "Alignment/LaserAlignment/interface/LASPeakFinder.h"
0003
0004
0005
0006
0007 LASPeakFinder::LASPeakFinder() {
0008 amplitudeThreshold = 0.;
0009 }
0010
0011
0012
0013
0014
0015
0016 bool LASPeakFinder::FindPeakIn(const LASModuleProfile& aProfile,
0017 std::pair<double, double>& result,
0018 TH1D* histogram,
0019 const double offset) {
0020
0021
0022 TF1 fitFunction("fitFunction", "gaus");
0023 fitFunction.SetLineColor(2);
0024 fitFunction.SetLineWidth(2);
0025
0026 std::pair<int, double> largestAmplitude(0, 0.);
0027 double anAmplitude = 0.;
0028
0029
0030 const int approxOffset = static_cast<int>(offset);
0031
0032
0033 const unsigned int meanPosition = 256 + approxOffset;
0034
0035
0036 const unsigned int halfWindowSize = 33;
0037
0038
0039
0040
0041 double sum0 = 0.;
0042 for (unsigned int strip = meanPosition - halfWindowSize; strip < meanPosition + halfWindowSize; ++strip) {
0043 anAmplitude = aProfile.GetValue(strip);
0044 sum0 += anAmplitude;
0045 histogram->SetBinContent(1 + strip, anAmplitude);
0046 if (anAmplitude > largestAmplitude.second) {
0047 largestAmplitude.first = strip;
0048 largestAmplitude.second = anAmplitude;
0049 }
0050 }
0051
0052
0053
0054 double sum1 = 0., sum2 = 0.;
0055 int nStrips = 0;
0056 for (unsigned int strip = 0; strip < 512; ++strip) {
0057 if (strip < meanPosition - halfWindowSize || strip > meanPosition + halfWindowSize) {
0058 anAmplitude = aProfile.GetValue(strip);
0059 sum0 += anAmplitude;
0060 sum1 += anAmplitude;
0061 sum2 += pow(anAmplitude, 2);
0062 nStrips++;
0063 }
0064 }
0065
0066
0067 const double noise = sqrt(1. / (nStrips - 1) * (sum2 - pow(sum1, 2) / nStrips));
0068
0069
0070 if (fabs(sum0) < 1.e-3) {
0071 std::cerr << " [LASPeakFinder::FindPeakIn] ** WARNING: Empty profile (sum=" << sum0 << ")." << std::endl;
0072 return false;
0073 }
0074
0075
0076 if (largestAmplitude.second < amplitudeThreshold * noise) {
0077 std::cerr << " [LASPeakFinder::FindPeakIn] ** WARNING: No reasonably large peak." << std::endl;
0078 return false;
0079 }
0080
0081
0082 fitFunction.SetParameter(0, largestAmplitude.second);
0083 fitFunction.SetParameter(1, largestAmplitude.first);
0084 fitFunction.SetParameter(2, 3.);
0085
0086
0087 fitFunction.SetParLimits(
0088 0, largestAmplitude.second * 0.3, largestAmplitude.second * 1.8);
0089 fitFunction.SetParLimits(
0090 1, largestAmplitude.first - 12, largestAmplitude.first + 12);
0091 fitFunction.SetParLimits(2, 0.5, 8.);
0092
0093
0094 histogram->Fit(&fitFunction, "QWEB", "", largestAmplitude.first - 12, largestAmplitude.first + 12);
0095
0096
0097 result.first = fitFunction.GetParameter(1);
0098 result.second = fitFunction.GetParError(1);
0099
0100 return true;
0101 }
0102
0103
0104
0105
0106 void LASPeakFinder::SetAmplitudeThreshold(double aThreshold) { amplitudeThreshold = aThreshold; }