Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:45:03

0001 
0002 #include "Alignment/LaserAlignment/interface/LASPeakFinder.h"
0003 
0004 ///
0005 ///
0006 ///
0007 LASPeakFinder::LASPeakFinder() {
0008   amplitudeThreshold = 0.;  // default
0009 }
0010 
0011 ///
0012 /// set a profile to work on and start peak finding;
0013 /// the pair<> will return mean/meanError (in strips);
0014 /// offset is necessary for tob modules which are hit off-center
0015 ///
0016 bool LASPeakFinder::FindPeakIn(const LASModuleProfile& aProfile,
0017                                std::pair<double, double>& result,
0018                                TH1D* histogram,
0019                                const double offset) {
0020   // this function will be propagated to the histo output file,
0021   // therefore we do some cosmetics already here
0022   TF1 fitFunction("fitFunction", "gaus");
0023   fitFunction.SetLineColor(2);
0024   fitFunction.SetLineWidth(2);
0025 
0026   std::pair<int, double> largestAmplitude(0, 0.);  // strip, amplitude
0027   double anAmplitude = 0.;
0028 
0029   // need approximate position -> cast to int
0030   const int approxOffset = static_cast<int>(offset);
0031 
0032   // expected beam position (in strips)
0033   const unsigned int meanPosition = 256 + approxOffset;
0034 
0035   // backplane "alignment hole" approx. half size (in strips)
0036   const unsigned int halfWindowSize = 33;
0037 
0038   // loop over the strips in the "alignment hole"
0039   // to fill the histogram
0040   // and determine fit parameter estimates
0041   double sum0 = 0.;  // this is the sum of all amplitudes
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   // loop outside the "alignment hole"
0053   // to determine the noise level = sqrt(variance)
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   // noise as sqrt of the amplitude variance
0067   const double noise = sqrt(1. / (nStrips - 1) * (sum2 - pow(sum1, 2) / nStrips));
0068 
0069   // empty profile?
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   // no reasonable peak?
0076   if (largestAmplitude.second < amplitudeThreshold * noise) {
0077     std::cerr << " [LASPeakFinder::FindPeakIn] ** WARNING: No reasonably large peak." << std::endl;
0078     return false;
0079   }
0080 
0081   // prepare fit function: starting values..
0082   fitFunction.SetParameter(0, largestAmplitude.second);  // amp
0083   fitFunction.SetParameter(1, largestAmplitude.first);   // mean
0084   fitFunction.SetParameter(2, 3.);                       // width
0085 
0086   // ..and parameter limits
0087   fitFunction.SetParLimits(
0088       0, largestAmplitude.second * 0.3, largestAmplitude.second * 1.8);  // amp of the order of the peak height
0089   fitFunction.SetParLimits(
0090       1, largestAmplitude.first - 12, largestAmplitude.first + 12);  // mean around the peak maximum
0091   fitFunction.SetParLimits(2, 0.5, 8.);                              // reasonable width
0092 
0093   // fit options: Quiet / all Weights=1 / better Errors / enable fix&Bounds on parameters
0094   histogram->Fit(&fitFunction, "QWEB", "", largestAmplitude.first - 12, largestAmplitude.first + 12);
0095 
0096   // prepare output
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; }