Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-01-17 23:31:52

0001 #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripAPVRestorer.h"
0002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/Framework/interface/ConsumesCollector.h"
0005 
0006 #include <cmath>
0007 #include <iostream>
0008 #include <algorithm>
0009 
0010 SiStripAPVRestorer::SiStripAPVRestorer(const edm::ParameterSet& conf, edm::ConsumesCollector iC)
0011     : siStripRawDigiToken_(iC.consumes<edm::DetSetVector<SiStripRawDigi>>(edm::InputTag("siStripDigis", "VirginRaw"))),
0012       siStripProcessedRawDigiToken_(
0013           iC.consumes<edm::DetSetVector<SiStripProcessedRawDigi>>(edm::InputTag("MEANAPVCM"))),
0014       qualityToken_(iC.esConsumes<SiStripQuality, SiStripQualityRcd>()),
0015       noiseToken_(iC.esConsumes<SiStripNoises, SiStripNoisesRcd>()),
0016       pedestalToken_(iC.esConsumes<SiStripPedestals, SiStripPedestalsRcd>()),
0017       forceNoRestore_(conf.getParameter<bool>("ForceNoRestore")),
0018       inspectAlgo_(conf.getParameter<std::string>("APVInspectMode")),
0019       restoreAlgo_(conf.getParameter<std::string>("APVRestoreMode")),
0020       // CM map settings
0021       useRealMeanCM_(conf.getParameter<bool>("useRealMeanCM")),
0022       meanCM_(conf.getParameter<int32_t>("MeanCM")),
0023       deltaCMThreshold_(conf.getParameter<uint32_t>("DeltaCMThreshold")),
0024       // abnormal baseline inspect
0025       fraction_(conf.getParameter<double>("Fraction")),
0026       deviation_(conf.getParameter<uint32_t>("Deviation")),
0027       // nullInspect
0028       restoreThreshold_(conf.getParameter<double>("restoreThreshold")),
0029       // baseline and saturation inspect
0030       nSaturatedStrip_(conf.getParameter<uint32_t>("nSaturatedStrip")),
0031       // FlatRegionsFinder
0032       nSigmaNoiseDerTh_(conf.getParameter<uint32_t>("nSigmaNoiseDerTh")),
0033       consecThreshold_(conf.getParameter<uint32_t>("consecThreshold")),
0034       nSmooth_(conf.getParameter<uint32_t>("nSmooth")),
0035       distortionThreshold_(conf.getParameter<uint32_t>("distortionThreshold")),
0036       applyBaselineCleaner_(conf.getParameter<bool>("ApplyBaselineCleaner")),
0037       cleaningSequence_(conf.getParameter<uint32_t>("CleaningSequence")),
0038       // cleaner_localmin
0039       slopeX_(conf.getParameter<int32_t>("slopeX")),
0040       slopeY_(conf.getParameter<int32_t>("slopeY")),
0041       // cleaner_monotony
0042       hitStripThreshold_(conf.getParameter<uint32_t>("hitStripThreshold")),
0043       // baseline follower (restore)
0044       minStripsToFit_(conf.getParameter<uint32_t>("minStripsToFit")),
0045       applyBaselineRejection_(conf.getParameter<bool>("ApplyBaselineRejection")),
0046       filteredBaselineMax_(conf.getParameter<double>("filteredBaselineMax")),
0047       filteredBaselineDerivativeSumSquare_(conf.getParameter<double>("filteredBaselineDerivativeSumSquare")),
0048       // derivative follower restorer
0049       gradient_threshold_(conf.getParameter<int>("discontinuityThreshold")),
0050       last_gradient_(conf.getParameter<int>("lastGradient")),
0051       size_window_(conf.getParameter<int>("sizeWindow")),
0052       width_cluster_(conf.getParameter<int>("widthCluster")) {
0053   if (restoreAlgo_ == "BaselineFollower" && inspectAlgo_ != "BaselineFollower" && inspectAlgo_ != "Hybrid")
0054     throw cms::Exception("Incompatible Algorithm")
0055         << "The BaselineFollower restore method requires the BaselineFollower (or Hybrid) inspect method";
0056 }
0057 
0058 void SiStripAPVRestorer::init(const edm::EventSetup& es) {
0059   if (noiseWatcher_.check(es)) {
0060     noiseHandle = &es.getData(noiseToken_);
0061   }
0062   if (qualityWatcher_.check(es)) {
0063     qualityHandle = &es.getData(qualityToken_);
0064   }
0065   if (pedestalWatcher_.check(es)) {
0066     pedestalHandle = &es.getData(pedestalToken_);
0067   }
0068 }
0069 
0070 uint16_t SiStripAPVRestorer::inspectAndRestore(uint32_t detId,
0071                                                uint16_t firstAPV,
0072                                                const digivector_t& rawDigisPedSubtracted,
0073                                                digivector_t& processedRawDigi,
0074                                                const medians_t& vmedians) {
0075   auto nAPVFlagged = inspect(detId, firstAPV, rawDigisPedSubtracted, vmedians);
0076   restore(firstAPV, processedRawDigi);
0077   return nAPVFlagged;
0078 }
0079 
0080 uint16_t SiStripAPVRestorer::inspect(uint32_t detId,
0081                                      uint16_t firstAPV,
0082                                      const digivector_t& digis,
0083                                      const medians_t& vmedians) {
0084   detId_ = detId;
0085 
0086   apvFlags_.assign(6, "");
0087   apvFlagsBoolOverride_.assign(6, false);
0088   apvFlagsBool_.clear();
0089   median_.assign(6, -999);
0090   badAPVs_.assign(6, false);
0091   smoothedMaps_.clear();
0092   baselineMap_.clear();
0093   for (const auto& med : vmedians) {
0094     auto iAPV = med.first;
0095     median_[iAPV] = med.second;
0096     badAPVs_[iAPV] = qualityHandle->IsApvBad(detId_, iAPV);
0097   }
0098 
0099   uint16_t nAPVFlagged = 0;
0100   if (inspectAlgo_ == "Hybrid")
0101     nAPVFlagged = hybridFormatInspect(firstAPV, digis);
0102   else if (inspectAlgo_ == "HybridEmulation")
0103     nAPVFlagged = hybridEmulationInspect(firstAPV, digis);
0104   else if (inspectAlgo_ == "BaselineFollower")
0105     nAPVFlagged = baselineFollowerInspect(firstAPV, digis);
0106   else if (inspectAlgo_ == "AbnormalBaseline")
0107     nAPVFlagged = abnormalBaselineInspect(firstAPV, digis);
0108   else if (inspectAlgo_ == "Null")
0109     nAPVFlagged = nullInspect(firstAPV, digis);
0110   else if (inspectAlgo_ == "BaselineAndSaturation")
0111     nAPVFlagged = baselineAndSaturationInspect(firstAPV, digis);
0112   else if (inspectAlgo_ == "DerivativeFollower")
0113     nAPVFlagged = forceRestoreInspect(firstAPV, digis);
0114   else
0115     throw cms::Exception("Unregistered Inspect Algorithm")
0116         << "SiStripAPVRestorer possibilities: (Null), (AbnormalBaseline), (BaselineFollower), (BaselineAndSaturation), "
0117            "(DerivativeFollower), (Hybrid), (HybridEmulation)";
0118 
0119   // update boolean version of APV flags
0120   for (size_t i = 0; i < apvFlags_.size(); ++i) {
0121     apvFlagsBool_.push_back((!apvFlags_[i].empty()) && (!apvFlagsBoolOverride_[i]));
0122   }
0123 
0124   return nAPVFlagged;
0125 }
0126 
0127 void SiStripAPVRestorer::restore(uint16_t firstAPV, digivector_t& digis) {
0128   if (forceNoRestore_)
0129     return;
0130 
0131   for (uint16_t iAPV = firstAPV; iAPV < digis.size() / nTotStripsPerAPV + firstAPV; ++iAPV) {
0132     auto algoToUse = apvFlags_[iAPV];
0133     if (!algoToUse.empty()) {
0134       if (algoToUse == "Flat") {
0135         flatRestore(iAPV, firstAPV, digis);
0136       } else if (algoToUse == "BaselineFollower") {
0137         baselineFollowerRestore(iAPV, firstAPV, median_[iAPV], digis);
0138       } else if (algoToUse == "DerivativeFollower") {
0139         derivativeFollowerRestore(iAPV, firstAPV, digis);
0140       } else {
0141         throw cms::Exception("Unregistered Restore Algorithm")
0142             << "SiStripAPVRestorer possibilities: (Flat), (BaselineFollower), (DerivativeFollower)";
0143       }
0144     }
0145   }
0146 }
0147 
0148 // Inspect method implementations
0149 
0150 inline uint16_t SiStripAPVRestorer::hybridFormatInspect(uint16_t firstAPV, const digivector_t& digis) {
0151   digivector_t singleAPVdigi;
0152   singleAPVdigi.reserve(nTotStripsPerAPV);
0153   uint16_t nAPVflagged = 0;
0154   for (uint16_t iAPV = firstAPV; iAPV < digis.size() / nTotStripsPerAPV + firstAPV; ++iAPV) {
0155     digimap_t smoothedmap;
0156     if (!badAPVs_[iAPV]) {
0157       uint16_t stripsPerAPV = 0;
0158       singleAPVdigi.clear();
0159       for (int16_t strip = (iAPV - firstAPV) * nTotStripsPerAPV; strip < (iAPV - firstAPV + 1) * nTotStripsPerAPV;
0160            ++strip) {
0161         const uint16_t adc = digis[strip];
0162         singleAPVdigi.push_back(adc);
0163         if (adc > 0)
0164           ++stripsPerAPV;
0165       }
0166 
0167       if (stripsPerAPV > 64) {
0168         flatRegionsFinder(singleAPVdigi, smoothedmap, iAPV);
0169         apvFlags_[iAPV] = restoreAlgo_;  //specify any algo to make the restore
0170         nAPVflagged++;
0171       }
0172     }
0173     smoothedMaps_.emplace_hint(smoothedMaps_.end(), iAPV, std::move(smoothedmap));
0174   }
0175   return nAPVflagged;
0176 }
0177 
0178 inline uint16_t SiStripAPVRestorer::baselineFollowerInspect(uint16_t firstAPV, const digivector_t& digis) {
0179   digivector_t singleAPVdigi;
0180   singleAPVdigi.reserve(nTotStripsPerAPV);
0181   uint16_t nAPVflagged = 0;
0182 
0183   auto itCMMap = std::end(meanCMmap_);
0184   if (useRealMeanCM_)
0185     itCMMap = meanCMmap_.find(detId_);
0186 
0187   for (uint16_t iAPV = firstAPV; iAPV < firstAPV + digis.size() / nTotStripsPerAPV; ++iAPV) {
0188     digimap_t smoothedmap;
0189     if (!badAPVs_[iAPV]) {
0190       float MeanAPVCM = meanCM_;
0191       if (useRealMeanCM_ && (std::end(meanCMmap_) != itCMMap))
0192         MeanAPVCM = itCMMap->second[iAPV];
0193 
0194       singleAPVdigi.clear();
0195       std::copy_n(std::begin(digis) + nTotStripsPerAPV * (iAPV - firstAPV),
0196                   nTotStripsPerAPV,
0197                   std::back_inserter(singleAPVdigi));
0198 
0199       const float DeltaCM = median_[iAPV] - MeanAPVCM;
0200       if ((DeltaCM < 0) && (std::abs(DeltaCM) > deltaCMThreshold_)) {
0201         if (!flatRegionsFinder(singleAPVdigi, smoothedmap, iAPV)) {
0202           apvFlags_[iAPV] = restoreAlgo_;  //specify any algo to make the restore
0203           ++nAPVflagged;
0204         }
0205       }
0206     }
0207     smoothedMaps_.emplace_hint(smoothedMaps_.end(), iAPV, std::move(smoothedmap));
0208   }
0209   return nAPVflagged;
0210 }
0211 
0212 inline uint16_t SiStripAPVRestorer::forceRestoreInspect(uint16_t firstAPV, const digivector_t& digis) {
0213   uint16_t nAPVflagged = 0;
0214   for (uint16_t iAPV = firstAPV; iAPV < firstAPV + digis.size() / nTotStripsPerAPV; ++iAPV) {
0215     if (!badAPVs_[iAPV]) {
0216       apvFlags_[iAPV] = restoreAlgo_;  //specify any algo to make the restore
0217       nAPVflagged++;
0218     }
0219   }
0220   return nAPVflagged;
0221 }
0222 
0223 inline uint16_t SiStripAPVRestorer::baselineAndSaturationInspect(uint16_t firstAPV, const digivector_t& digis) {
0224   digivector_t singleAPVdigi;
0225   singleAPVdigi.reserve(nTotStripsPerAPV);
0226   uint16_t nAPVflagged = 0;
0227 
0228   auto itCMMap = std::end(meanCMmap_);
0229   if (useRealMeanCM_)
0230     itCMMap = meanCMmap_.find(detId_);
0231 
0232   for (uint16_t iAPV = firstAPV; iAPV < firstAPV + digis.size() / nTotStripsPerAPV; ++iAPV) {
0233     if (!badAPVs_[iAPV]) {
0234       float MeanAPVCM = meanCM_;
0235       if (useRealMeanCM_ && (std::end(meanCMmap_) != itCMMap))
0236         MeanAPVCM = itCMMap->second[iAPV];
0237 
0238       singleAPVdigi.clear();
0239       uint16_t nSatStrip = 0;
0240       for (int16_t strip = (iAPV - firstAPV) * nTotStripsPerAPV; strip < (iAPV - firstAPV + 1) * nTotStripsPerAPV;
0241            ++strip) {
0242         const uint16_t digi = digis[strip];
0243         singleAPVdigi.push_back(digi);
0244         if (digi >= 1023)
0245           ++nSatStrip;
0246       }
0247 
0248       const float DeltaCM = median_[iAPV] - MeanAPVCM;
0249       if ((DeltaCM < 0) && (std::abs(DeltaCM) > deltaCMThreshold_) && (nSatStrip >= nSaturatedStrip_)) {
0250         apvFlags_[iAPV] = restoreAlgo_;  //specify any algo to make the restore
0251         ++nAPVflagged;
0252       }
0253     }
0254   }
0255   return nAPVflagged;
0256 }
0257 
0258 inline uint16_t SiStripAPVRestorer::abnormalBaselineInspect(uint16_t firstAPV, const digivector_t& digis) {
0259   SiStripQuality::Range detQualityRange = qualityHandle->getRange(detId_);
0260 
0261   uint16_t nAPVflagged = 0;
0262 
0263   auto itCMMap = std::end(meanCMmap_);
0264   if (useRealMeanCM_)
0265     itCMMap = meanCMmap_.find(detId_);
0266 
0267   int devCount = 0, qualityCount = 0, minstrip = 0;
0268   for (uint16_t iAPV = firstAPV; iAPV < firstAPV + digis.size() / nTotStripsPerAPV; ++iAPV) {
0269     if (!badAPVs_[iAPV]) {
0270       float MeanAPVCM = meanCM_;
0271       if (useRealMeanCM_ && (std::end(meanCMmap_) != itCMMap))
0272         MeanAPVCM = itCMMap->second[iAPV];
0273 
0274       for (uint16_t istrip = iAPV * nTotStripsPerAPV; istrip < (iAPV + 1) * nTotStripsPerAPV; ++istrip) {
0275         const auto fs = static_cast<int>(digis[istrip - firstAPV * nTotStripsPerAPV]);
0276         if (!qualityHandle->IsStripBad(detQualityRange, istrip)) {
0277           qualityCount++;
0278           if (std::abs(fs - MeanAPVCM) > static_cast<int>(deviation_)) {
0279             devCount++;
0280             minstrip = std::min(fs, minstrip);
0281           }
0282         }
0283       }
0284       if (devCount > fraction_ * qualityCount) {
0285         apvFlags_[iAPV] = restoreAlgo_;  //specify any algo to make the restore
0286         nAPVflagged++;
0287       }
0288     }
0289   }
0290   return nAPVflagged;
0291 }
0292 
0293 inline uint16_t SiStripAPVRestorer::nullInspect(uint16_t firstAPV, const digivector_t& digis) {
0294   SiStripQuality::Range detQualityRange = qualityHandle->getRange(detId_);
0295 
0296   uint16_t nAPVflagged = 0;
0297   for (uint16_t iAPV = firstAPV; iAPV < firstAPV + digis.size() / nTotStripsPerAPV; ++iAPV) {
0298     if (!badAPVs_[iAPV]) {
0299       int zeroCount = 0, qualityCount = 0;
0300       for (uint16_t istrip = iAPV * nTotStripsPerAPV; istrip < (iAPV + 1) * nTotStripsPerAPV; ++istrip) {
0301         const auto fs = static_cast<int>(digis[istrip - firstAPV * nTotStripsPerAPV]);
0302         if (!qualityHandle->IsStripBad(detQualityRange, istrip)) {
0303           qualityCount++;
0304           if (fs < 1)
0305             zeroCount++;
0306         }
0307       }
0308       if (zeroCount > restoreThreshold_ * qualityCount) {
0309         apvFlags_[iAPV] = restoreAlgo_;  //specify any algo to make the restore
0310         nAPVflagged++;
0311       }
0312     }
0313   }
0314   return nAPVflagged;
0315 }
0316 
0317 inline uint16_t SiStripAPVRestorer::hybridEmulationInspect(uint16_t firstAPV, const digivector_t& digis) {
0318   uint16_t nAPVflagged = 0;
0319 
0320   auto itCMMap = std::end(meanCMmap_);
0321   if (useRealMeanCM_)
0322     itCMMap = meanCMmap_.find(detId_);
0323 
0324   for (uint16_t iAPV = firstAPV; iAPV < firstAPV + digis.size() / nTotStripsPerAPV; ++iAPV) {
0325     if (!badAPVs_[iAPV]) {
0326       float MeanAPVCM = meanCM_;
0327       if (useRealMeanCM_ && (std::end(meanCMmap_) != itCMMap))
0328         MeanAPVCM = itCMMap->second[iAPV];
0329 
0330       const float DeltaCM = median_[iAPV] - (MeanAPVCM + 1024) / 2;
0331       if ((DeltaCM < 0) && (std::abs(DeltaCM) > deltaCMThreshold_ / 2)) {
0332         apvFlags_[iAPV] = "HybridEmulation";
0333         ++nAPVflagged;
0334       }
0335     }
0336   }
0337   return nAPVflagged;
0338 }
0339 
0340 // Restore method implementations
0341 
0342 inline void SiStripAPVRestorer::baselineFollowerRestore(uint16_t apvN,
0343                                                         uint16_t firstAPV,
0344                                                         float median,
0345                                                         digivector_t& digis) {
0346   digivector_t baseline(size_t(nTotStripsPerAPV), 0);
0347 
0348   //============================= Find Flat Regions & Interpolating the baseline & subtracting the baseline  =================
0349   if (!smoothedMaps_.empty()) {
0350     auto itSmootedMap = smoothedMaps_.find(apvN);
0351     baselineFollower(itSmootedMap->second, baseline, median);
0352   } else {
0353     //median=0;
0354     digivector_t singleAPVdigi;
0355     singleAPVdigi.reserve(nTotStripsPerAPV);
0356     std::copy_n(
0357         std::begin(digis) + nTotStripsPerAPV * (apvN - firstAPV), nTotStripsPerAPV, std::back_inserter(singleAPVdigi));
0358     digimap_t smoothedpoints;
0359     flatRegionsFinder(singleAPVdigi, smoothedpoints, apvN);
0360     baselineFollower(smoothedpoints, baseline, median);
0361   }
0362 
0363   if (applyBaselineRejection_) {
0364     if (checkBaseline(baseline))
0365       apvFlagsBoolOverride_[apvN] = true;
0366   }
0367 
0368   //============================= subtracting the baseline =============================================
0369   for (int16_t itStrip = 0; itStrip < nTotStripsPerAPV; ++itStrip) {
0370     digis[nTotStripsPerAPV * (apvN - firstAPV) + itStrip] -= baseline[itStrip] - median;
0371   }
0372 
0373   //============================= storing baseline to the map =============================================
0374   baselineMap_.emplace_hint(baselineMap_.end(), apvN, std::move(baseline));
0375 }
0376 
0377 inline void SiStripAPVRestorer::flatRestore(uint16_t apvN, uint16_t firstAPV, digivector_t& digis) {
0378   digivector_t baseline(size_t(nTotStripsPerAPV), 150);
0379   baseline[0] = 0;
0380   baseline[nTotStripsPerAPV - 1] = 0;
0381 
0382   for (int16_t itStrip = 0; itStrip < nTotStripsPerAPV; ++itStrip) {
0383     digis[nTotStripsPerAPV * (apvN - firstAPV) + itStrip] = baseline[itStrip];
0384   }
0385   baselineMap_.emplace_hint(baselineMap_.end(), apvN, std::move(baseline));
0386 }
0387 
0388 // Baseline calculation implementation =========================================
0389 bool inline SiStripAPVRestorer::flatRegionsFinder(const digivector_t& adcs, digimap_t& smoothedpoints, uint16_t apvN) {
0390   smoothedpoints.clear();
0391 
0392   SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detId_);
0393 
0394   //============================= Height above local minimum ===============================
0395   std::vector<float> adcsLocalMinSubtracted(size_t(nTotStripsPerAPV), 0);
0396   for (uint32_t istrip = 0; istrip < nTotStripsPerAPV; ++istrip) {
0397     float localmin = 999.9;
0398     for (uint16_t jstrip = std::max(0, static_cast<int>(istrip - nSmooth_ / 2));
0399          jstrip < std::min((int)nTotStripsPerAPV, static_cast<int>(istrip + nSmooth_ / 2));
0400          ++jstrip) {
0401       const float nextvalue = adcs[jstrip];
0402       if (nextvalue < localmin)
0403         localmin = nextvalue;
0404     }
0405     adcsLocalMinSubtracted[istrip] = adcs[istrip] - localmin;
0406   }
0407 
0408   //============================= Find regions with stable slopes ========================
0409   digimap_t consecpoints;
0410   std::vector<uint16_t> nConsStrip;
0411   //Creating maps with all the neighborhood strip and putting in a nCosntStip vector how many we have
0412   uint16_t consecStrips = 0;
0413   for (uint32_t istrip = 0; istrip < nTotStripsPerAPV; ++istrip) {
0414     const int16_t adc = adcs[istrip];
0415     //if( adcsLocalMinSubtracted[istrip] < nSigmaNoiseDerTh_ * (float)noiseHandle->getNoise(istrip+apvN*128,detNoiseRange) && (adc - median) < hitStripThreshold_){
0416     if (adcsLocalMinSubtracted[istrip] <
0417         nSigmaNoiseDerTh_ * (float)noiseHandle->getNoiseFast(istrip + apvN * nTotStripsPerAPV, detNoiseRange)) {
0418       consecpoints.emplace_hint(consecpoints.end(), istrip, adc);
0419       ++consecStrips;
0420     } else if (consecStrips > 0) {
0421       nConsStrip.push_back(consecStrips);
0422       consecStrips = 0;
0423     }
0424   }
0425   //to cope with the last flat region of the APV
0426   if (consecStrips > 0)
0427     nConsStrip.push_back(consecStrips);
0428 
0429   //removing from the map the fist and last points in wide flat regions and erasing from the map too small regions
0430   auto itConsecpoints = consecpoints.begin();
0431   float MinSmoothValue = 20000., MaxSmoothValue = 0.;
0432   for (const auto consecStrips : nConsStrip) {
0433     if (consecStrips >= consecThreshold_) {
0434       ++itConsecpoints;  //skipping first point
0435       uint16_t nFirstStrip = itConsecpoints->first;
0436       uint16_t nLastStrip;
0437       float smoothValue = 0.0;
0438       float stripCount = 1;
0439       for (uint16_t n = 0; n < consecStrips - 2; ++n) {
0440         smoothValue += itConsecpoints->second;
0441         if (stripCount == consecThreshold_) {
0442           smoothValue /= (float)stripCount;
0443           nLastStrip = nFirstStrip + stripCount - 1;
0444           smoothedpoints.emplace_hint(smoothedpoints.end(), nFirstStrip, smoothValue);
0445           smoothedpoints.emplace_hint(smoothedpoints.end(), nLastStrip, smoothValue);
0446           if (smoothValue > MaxSmoothValue)
0447             MaxSmoothValue = smoothValue;
0448           if (smoothValue < MinSmoothValue)
0449             MinSmoothValue = smoothValue;
0450           nFirstStrip = nLastStrip + 1;
0451           smoothValue = 0;
0452           stripCount = 0;
0453         }
0454         ++stripCount;
0455         ++itConsecpoints;
0456       }
0457       ++itConsecpoints;  //and putting the pointer to the new seies of point
0458 
0459       if (stripCount > 1) {
0460         //if(smoothValue>0){
0461         --stripCount;
0462         smoothValue /= static_cast<float>(stripCount);
0463         nLastStrip = nFirstStrip + stripCount - 1;
0464         smoothedpoints.emplace_hint(smoothedpoints.end(), nFirstStrip, smoothValue);
0465         smoothedpoints.emplace_hint(smoothedpoints.end(), nLastStrip, smoothValue);
0466         if (smoothValue > MaxSmoothValue)
0467           MaxSmoothValue = smoothValue;
0468         if (smoothValue < MinSmoothValue)
0469           MinSmoothValue = smoothValue;
0470       }
0471     } else {
0472       for (int n = 0; n < consecStrips; ++n)
0473         ++itConsecpoints;
0474     }
0475   }
0476 
0477   if ((MaxSmoothValue - MinSmoothValue) > distortionThreshold_) {
0478     if (applyBaselineCleaner_)
0479       baselineCleaner(adcs, smoothedpoints, apvN);
0480     return false;
0481   }
0482   return true;
0483 }
0484 
0485 void inline SiStripAPVRestorer::baselineCleaner(const digivector_t& adcs, digimap_t& smoothedpoints, uint16_t apvN) {
0486   if (cleaningSequence_ == 0) {  //default sequence used up to now
0487     cleaner_HighSlopeChecker(smoothedpoints);
0488     cleaner_LocalMinimumAdder(adcs, smoothedpoints, apvN);
0489   } else if (cleaningSequence_ == 1) {
0490     cleaner_LocalMinimumAdder(adcs, smoothedpoints, apvN);
0491     cleaner_HighSlopeChecker(smoothedpoints);
0492     cleaner_MonotonyChecker(smoothedpoints);
0493   } else if (cleaningSequence_ == 2) {
0494     cleaner_HighSlopeChecker(smoothedpoints);
0495   } else if (cleaningSequence_ == 3) {
0496     cleaner_LocalMinimumAdder(adcs, smoothedpoints, apvN);
0497     cleaner_HighSlopeChecker(smoothedpoints);
0498   } else {
0499     cleaner_HighSlopeChecker(smoothedpoints);
0500     cleaner_LocalMinimumAdder(adcs, smoothedpoints, apvN);
0501   }
0502 }
0503 
0504 void inline SiStripAPVRestorer::cleaner_MonotonyChecker(digimap_t& smoothedpoints) {
0505   //Removing points without monotony
0506   //--------------------------------------------------------------------------------------------------
0507   if (smoothedpoints.size() < 3)
0508     return;
0509 
0510   auto it = smoothedpoints.begin();
0511   while (smoothedpoints.size() > 3 && it != --(--smoothedpoints.end())) {  //while we are not at the last point
0512     // get info about current and next points
0513     auto it2 = it;
0514     const float adc1 = it2->second;
0515     const float adc2 = (++it2)->second;
0516     const float adc3 = (++it2)->second;
0517 
0518     if ((adc2 - adc1) > hitStripThreshold_ && (adc2 - adc3) > hitStripThreshold_) {
0519       smoothedpoints.erase(--it2);
0520     } else {
0521       ++it;
0522     }
0523   }
0524 }
0525 
0526 void inline SiStripAPVRestorer::cleaner_LocalMinimumAdder(const digivector_t& adcs,
0527                                                           digimap_t& smoothedpoints,
0528                                                           uint16_t apvN) {
0529   SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detId_);
0530   //inserting extra point is case of local minimum
0531   //--------------------------------------------------------------------------------------------------
0532   // these should be reset now for the point-insertion that follows
0533 
0534   if (smoothedpoints.size() >= 2) {
0535     for (auto it = smoothedpoints.begin(); it != --smoothedpoints.end(); ++it) {
0536       auto itNext = it;
0537       ++itNext;
0538       const float strip1 = it->first;
0539       const float strip2 = itNext->first;
0540       const float adc1 = it->second;
0541       const float adc2 = itNext->second;
0542       const float m = (adc2 - adc1) / (strip2 - strip1);
0543 
0544       //2,4
0545       if ((strip2 - strip1) > slopeX_ && std::abs(adc1 - adc2) > slopeY_) {
0546         float itStrip = 1;
0547         float strip = itStrip + strip1;
0548         while (strip < strip2) {
0549           const float adc = adcs[strip];
0550           const auto noise =
0551               static_cast<float>(noiseHandle->getNoiseFast(strip + apvN * nTotStripsPerAPV, detNoiseRange));
0552           if (adc < (adc1 + m * itStrip - 2 * noise)) {
0553             smoothedpoints.emplace_hint(itNext, strip, adc);
0554             ++it;
0555           }
0556           ++itStrip;
0557           ++strip;
0558         }
0559       }
0560     }
0561   }
0562 
0563   const uint16_t firstStripFlat = smoothedpoints.begin()->first;
0564   const uint16_t lastStripFlat = (--smoothedpoints.end())->first;
0565   const int16_t firstStripFlatADC = smoothedpoints.begin()->second;
0566   const int16_t lastStripFlatADC = (--smoothedpoints.end())->second;
0567   if (firstStripFlat > 3) {
0568     auto it = smoothedpoints.begin();
0569     float strip = 0;
0570     while (strip < firstStripFlat) {
0571       const float adc = adcs[strip];
0572       const auto noise = static_cast<float>(noiseHandle->getNoiseFast(strip + apvN * nTotStripsPerAPV, detNoiseRange));
0573       if (adc < (firstStripFlatADC - 2 * noise)) {
0574         smoothedpoints.emplace_hint(it, strip, adc);
0575         ++it;
0576       }
0577       ++strip;
0578     }
0579   }
0580   if (lastStripFlat < (nTotStripsPerAPV - 3)) {
0581     float strip = lastStripFlat + 1;
0582     while (strip < nTotStripsPerAPV) {
0583       const float adc = adcs[strip];
0584       const auto noise = static_cast<float>(noiseHandle->getNoiseFast(strip + apvN * nTotStripsPerAPV, detNoiseRange));
0585       if (adc < (lastStripFlatADC - 2 * noise)) {
0586         smoothedpoints.emplace_hint(smoothedpoints.end(), strip, adc);
0587       }
0588       ++strip;
0589     }
0590   }
0591 }
0592 
0593 void inline SiStripAPVRestorer::cleaner_HighSlopeChecker(digimap_t& smoothedpoints) {
0594   //Removing points in the slope is too high
0595   //----------------------------------------------------------------------------
0596   if (smoothedpoints.size() < 4)
0597     return;
0598 
0599   auto it = smoothedpoints.begin();
0600   while (smoothedpoints.size() > 2 && it != --smoothedpoints.end()) {  //while we are not at the last point
0601     //if(smoothedpoints.size() <2) break;
0602     // get info about current and next points
0603     auto itNext = it;
0604     ++itNext;
0605     const float strip1 = it->first;
0606     const float strip2 = itNext->first;
0607     const float adc1 = it->second;
0608     const float adc2 = itNext->second;
0609     const float m = (adc2 - adc1) / (strip2 - strip1);
0610 
0611     if (m > 2) {  // in case of large positive slope, remove next point and try again from same current point
0612       smoothedpoints.erase(itNext);
0613     } else if (m < -2) {  // in case of large negative slope, remove current point and either...
0614       // move to next point if we have reached the beginning (post-increment to avoid invalidating pointer during erase) or...
0615       if (it == smoothedpoints.begin())
0616         smoothedpoints.erase(it++);
0617       // try again from the previous point if we have not reached the beginning
0618       else
0619         smoothedpoints.erase(it--);
0620     } else {  // in case of a flat enough slope, continue on to the next point
0621       ++it;
0622     }
0623   }
0624 }
0625 
0626 void inline SiStripAPVRestorer::baselineFollower(const digimap_t& smoothedpoints,
0627                                                  digivector_t& baseline,
0628                                                  float median) {
0629   //if not enough points
0630   if (smoothedpoints.size() < minStripsToFit_) {
0631     baseline.insert(baseline.begin(), nTotStripsPerAPV, median);
0632   } else {
0633     baseline.assign(size_t(nTotStripsPerAPV), 0);
0634 
0635     const auto lastIt = --smoothedpoints.end();
0636     const uint16_t firstStripFlat = smoothedpoints.begin()->first;
0637     const uint16_t lastStripFlat = lastIt->first;
0638     const int16_t firstStripFlatADC = smoothedpoints.begin()->second;
0639     const int16_t lastStripFlatADC = lastIt->second;
0640 
0641     //adding here the costant line at the extremities
0642     baseline.erase(baseline.begin(), baseline.begin() + firstStripFlat);
0643     baseline.insert(baseline.begin(), firstStripFlat, firstStripFlatADC);
0644 
0645     baseline.erase(baseline.begin() + lastStripFlat, baseline.end());
0646     baseline.insert(baseline.end(), nTotStripsPerAPV - lastStripFlat, lastStripFlatADC);
0647 
0648     //IMPORTANT: the itSmoothedpointsEnd should be at least smaller than smoothedpoints.end() -1
0649     for (auto itSmoothedpoints = smoothedpoints.begin(); itSmoothedpoints != (--smoothedpoints.end());
0650          ++itSmoothedpoints) {
0651       auto itNextSmoothedpoints = itSmoothedpoints;
0652       ++itNextSmoothedpoints;
0653       const float strip1 = itSmoothedpoints->first;
0654       const float strip2 = itNextSmoothedpoints->first;
0655       const float adc1 = itSmoothedpoints->second;
0656       const float adc2 = itNextSmoothedpoints->second;
0657 
0658       baseline[strip1] = adc1;
0659       baseline[strip2] = adc2;
0660       const float m = (adc2 - adc1) / (strip2 - strip1);
0661       uint16_t itStrip = strip1 + 1;
0662       float stripadc = adc1 + m;
0663       while (itStrip < strip2) {
0664         baseline[itStrip] = stripadc;
0665         ++itStrip;
0666         stripadc += m;
0667       }
0668     }
0669   }
0670 }
0671 
0672 bool SiStripAPVRestorer::checkBaseline(const digivector_t& baseline) const {
0673   // The Savitzky-Golay (S-G) filter of any left length nL, right
0674   // length nR, and order m, and with an optional opt equals the
0675   // derivative order (0 for the magnitude, 1 for the first
0676   // derivative, etc.) can be calculated using the following
0677   // Mathematica code:
0678   //
0679   // SavitzkyGolay[m_?IntegerQ, {nL_?IntegerQ, nR_?IntegerQ},
0680   //   opt___] := Module[
0681   //   {a, d},
0682   //   d = If[opt === Null, 0, If[IntegerQ[opt], opt, 0]];
0683   //   a = Table[
0684   //     If[i == 0 && j == 0, 1, i^j], {i, -nL, nR}, {j, 0,
0685   //      m}]; (Inverse[Transpose[a].a].Transpose[a])[[d + 1]]];
0686   //
0687   // The following coefficients can be then calculated by:
0688   //
0689   // N[Join[Table[SavitzkyGolay[2, {k, 16}], {k, 0, 16}],
0690   //   Table[SavitzkyGolay[2, {16, k}], {k, 15, 0, -1}]]]
0691 
0692   // nLR = max(nL, nR)
0693   static const size_t savitzky_golay_n_l_r = 16;
0694   static const float savitzky_golay_coefficient[2 * savitzky_golay_n_l_r + 1][2 * savitzky_golay_n_l_r + 1] = {
0695       {0.422085,
0696        0.325077,
0697        0.23839,
0698        0.162023,
0699        0.0959752,
0700        0.0402477,
0701        -0.00515996,
0702        -0.0402477,
0703        -0.0650155,
0704        -0.0794634,
0705        -0.0835913,
0706        -0.0773994,
0707        -0.0608875,
0708        -0.0340557,
0709        0.00309598,
0710        0.0505676,
0711        0.108359},
0712       {0.315789,
0713        0.254902,
0714        0.19969,
0715        0.150155,
0716        0.106295,
0717        0.0681115,
0718        0.0356037,
0719        0.00877193,
0720        -0.0123839,
0721        -0.0278638,
0722        -0.0376677,
0723        -0.0417957,
0724        -0.0402477,
0725        -0.0330237,
0726        -0.0201238,
0727        -0.00154799,
0728        0.0227038,
0729        0.0526316},
0730       {0.234586,
0731        0.198496,
0732        0.165207,
0733        0.134719,
0734        0.107032,
0735        0.0821465,
0736        0.0600619,
0737        0.0407784,
0738        0.024296,
0739        0.0106148,
0740        -0.000265369,
0741        -0.00834439,
0742        -0.0136223,
0743        -0.0160991,
0744        -0.0157747,
0745        -0.0126493,
0746        -0.00672269,
0747        0.00200501,
0748        0.0135338},
0749       {0.172078,     0.153076,    0.135099,    0.118148,   0.102221,   0.0873206, 0.073445,
0750        0.0605947,    0.0487697,   0.0379699,   0.0281955,  0.0194463,  0.0117225, 0.00502392,
0751        -0.000649351, -0.00529733, -0.00892003, -0.0115174, -0.0130895, -0.0136364},
0752       {0.123659,  0.116431,   0.109144,    0.101798,    0.0943921,  0.0869268,  0.0794021,
0753        0.0718179, 0.0641743,  0.0564712,   0.0487087,   0.0408868,  0.0330054,  0.0250646,
0754        0.0170644, 0.00900473, 0.000885613, -0.00729294, -0.0155309, -0.0238283, -0.0321852},
0755       {0.0859684, 0.0868154,  0.0869565,   0.0863919,  0.0851214,  0.0831451, 0.080463,  0.0770751,
0756        0.0729814, 0.0681818,  0.0626765,   0.0564653,  0.0495483,  0.0419255, 0.0335968, 0.0245624,
0757        0.0148221, 0.00437606, -0.00677583, -0.0186335, -0.0311971, -0.0444664},
0758       {0.0565217, 0.0628458, 0.0680971,  0.0722756,   0.0753811,  0.0774139,  0.0783738, 0.0782609,
0759        0.0770751, 0.0748165, 0.071485,   0.0670807,   0.0616036,  0.0550536,  0.0474308, 0.0387352,
0760        0.0289667, 0.0181254, 0.00621118, -0.00677583, -0.0208357, -0.0359684, -0.0521739},
0761       {0.0334615, 0.0434281, 0.0521329, 0.0595759,  0.0657571,   0.0706765, 0.0743341,  0.07673,
0762        0.0778641, 0.0777364, 0.0763469, 0.0736957,  0.0697826,   0.0646078, 0.0581712,  0.0504728,
0763        0.0415126, 0.0312907, 0.0198069, 0.00706142, -0.00694588, -0.022215, -0.0387458, -0.0565385},
0764       {0.0153846, 0.0276923, 0.0386622,  0.0482943,  0.0565886,  0.0635452, 0.0691639, 0.0734448, 0.076388,
0765        0.0779933, 0.0782609, 0.0771906,  0.0747826,  0.0710368,  0.0659532, 0.0595318, 0.0517726, 0.0426756,
0766        0.0322408, 0.0204682, 0.00735786, -0.0070903, -0.0228763, -0.04,     -0.0584615},
0767       {0.001221,  0.0149451, 0.027326,  0.0383639,  0.0480586,   0.0564103,  0.0634188,  0.0690842, 0.0734066,
0768        0.0763858, 0.078022,  0.078315,  0.077265,   0.0748718,   0.0711355,  0.0660562,  0.0596337, 0.0518681,
0769        0.0427595, 0.0323077, 0.0205128, 0.00737485, -0.00710623, -0.0229304, -0.0400977, -0.0586081},
0770       {-0.00985222, 0.00463138, 0.0178098, 0.029683,  0.0402509, 0.0495137,   0.0574713,  0.0641236,  0.0694708,
0771        0.0735127,   0.0762494,  0.0776809, 0.0778073, 0.0766284, 0.0741442,   0.0703549,  0.0652604,  0.0588607,
0772        0.0511557,   0.0421456,  0.0318302, 0.0202097, 0.0072839, -0.00694708, -0.0224833, -0.0393247, -0.0574713},
0773       {-0.0184729, -0.00369458, 0.00984169, 0.0221359,   0.0331881,  0.0429982,  0.0515662,
0774        0.0588923,  0.0649762,   0.0698181,  0.073418,    0.0757758,  0.0768915,  0.0767652,
0775        0.0753968,  0.0727864,   0.0689339,  0.0638394,   0.0575028,  0.0499242,  0.0411035,
0776        0.0310408,  0.019736,    0.00718917, -0.00659972, -0.0216307, -0.0379037, -0.0554187},
0777       {-0.025139, -0.0103925,  0.00318873, 0.0156046,  0.0268552, 0.0369405, 0.0458605, 0.0536151,
0778        0.0602045, 0.0656285,   0.0698872,  0.0729806,  0.0749086, 0.0756714, 0.0752688, 0.0737009,
0779        0.0709677, 0.0670692,   0.0620054,  0.0557763,  0.0483818, 0.039822,  0.0300969, 0.0192065,
0780        0.0071508, -0.00607024, -0.0204566, -0.0360083, -0.0527253},
0781       {-0.0302419, -0.0157536, -0.00234785, 0.00997537, 0.021216,   0.0313741, 0.0404497, 0.0484427,
0782        0.0553532,  0.0611811,  0.0659264,   0.0695892,  0.0721695,  0.0736672, 0.0740823, 0.0734149,
0783        0.0716649,  0.0688324,  0.0649174,   0.0599198,  0.0538396,  0.0466769, 0.0384316, 0.0291038,
0784        0.0186934,  0.00720046, -0.00537502, -0.0190331, -0.0337736, -0.0495968},
0785       {-0.0340909, -0.0200147, -0.006937, 0.00514208,  0.0162226,  0.0263045, 0.0353878, 0.0434725,
0786        0.0505587,  0.0566463,  0.0617353, 0.0658257,   0.0689175,  0.0710107, 0.0721054, 0.0722014,
0787        0.0712989,  0.0693978,  0.0664981, 0.0625999,   0.057703,   0.0518076, 0.0449135, 0.0370209,
0788        0.0281297,  0.01824,    0.0073516, -0.00453534, -0.0174209, -0.031305, -0.0461877},
0789       {-0.0369318, -0.0233688, -0.0107221, 0.00100806, 0.0118218,   0.0217192,  0.0307001, 0.0387647,
0790        0.0459128,  0.0521444,  0.0574597,  0.0618585,  0.0653409,   0.0679069,  0.0695565, 0.0702896,
0791        0.0701063,  0.0690066,  0.0669905,  0.0640579,  0.0602089,   0.0554435,  0.0497617, 0.0431635,
0792        0.0356488,  0.0272177,  0.0178702,  0.0076063,  -0.00357405, -0.0156708, -0.028684, -0.0426136},
0793       {-0.038961, -0.025974,  -0.0138249,  -0.00251362, 0.00795978, 0.0175953, 0.026393,  0.0343527, 0.0414747,
0794        0.0477587, 0.0532049,  0.0578132,   0.0615836,   0.0645161,  0.0666108, 0.0678676, 0.0682866, 0.0678676,
0795        0.0666108, 0.0645161,  0.0615836,   0.0578132,   0.0532049,  0.0477587, 0.0414747, 0.0343527, 0.026393,
0796        0.0175953, 0.00795978, -0.00251362, -0.0138249,  -0.025974,  -0.038961},
0797       {-0.0426136, -0.028684, -0.0156708, -0.00357405, 0.0076063,  0.0178702,  0.0272177,  0.0356488,
0798        0.0431635,  0.0497617, 0.0554435,  0.0602089,   0.0640579,  0.0669905,  0.0690066,  0.0701063,
0799        0.0702896,  0.0695565, 0.0679069,  0.0653409,   0.0618585,  0.0574597,  0.0521444,  0.0459128,
0800        0.0387647,  0.0307001, 0.0217192,  0.0118218,   0.00100806, -0.0107221, -0.0233688, -0.0369318},
0801       {-0.0461877, -0.031305, -0.0174209, -0.00453534, 0.0073516, 0.01824,    0.0281297, 0.0370209,
0802        0.0449135,  0.0518076, 0.057703,   0.0625999,   0.0664981, 0.0693978,  0.0712989, 0.0722014,
0803        0.0721054,  0.0710107, 0.0689175,  0.0658257,   0.0617353, 0.0566463,  0.0505587, 0.0434725,
0804        0.0353878,  0.0263045, 0.0162226,  0.00514208,  -0.006937, -0.0200147, -0.0340909},
0805       {-0.0495968, -0.0337736, -0.0190331, -0.00537502, 0.00720046, 0.0186934, 0.0291038, 0.0384316,
0806        0.0466769,  0.0538396,  0.0599198,  0.0649174,   0.0688324,  0.0716649, 0.0734149, 0.0740823,
0807        0.0736672,  0.0721695,  0.0695892,  0.0659264,   0.0611811,  0.0553532, 0.0484427, 0.0404497,
0808        0.0313741,  0.021216,   0.00997537, -0.00234785, -0.0157536, -0.0302419},
0809       {-0.0527253, -0.0360083, -0.0204566, -0.00607024, 0.0071508, 0.0192065, 0.0300969, 0.039822,
0810        0.0483818,  0.0557763,  0.0620054,  0.0670692,   0.0709677, 0.0737009, 0.0752688, 0.0756714,
0811        0.0749086,  0.0729806,  0.0698872,  0.0656285,   0.0602045, 0.0536151, 0.0458605, 0.0369405,
0812        0.0268552,  0.0156046,  0.00318873, -0.0103925,  -0.025139},
0813       {-0.0554187, -0.0379037, -0.0216307, -0.00659972, 0.00718917, 0.019736,    0.0310408,
0814        0.0411035,  0.0499242,  0.0575028,  0.0638394,   0.0689339,  0.0727864,   0.0753968,
0815        0.0767652,  0.0768915,  0.0757758,  0.073418,    0.0698181,  0.0649762,   0.0588923,
0816        0.0515662,  0.0429982,  0.0331881,  0.0221359,   0.00984169, -0.00369458, -0.0184729},
0817       {-0.0574713, -0.0393247, -0.0224833, -0.00694708, 0.0072839, 0.0202097, 0.0318302, 0.0421456,  0.0511557,
0818        0.0588607,  0.0652604,  0.0703549,  0.0741442,   0.0766284, 0.0778073, 0.0776809, 0.0762494,  0.0735127,
0819        0.0694708,  0.0641236,  0.0574713,  0.0495137,   0.0402509, 0.029683,  0.0178098, 0.00463138, -0.00985222},
0820       {-0.0586081, -0.0400977, -0.0229304, -0.00710623, 0.00737485, 0.0205128, 0.0323077, 0.0427595, 0.0518681,
0821        0.0596337,  0.0660562,  0.0711355,  0.0748718,   0.077265,   0.078315,  0.078022,  0.0763858, 0.0734066,
0822        0.0690842,  0.0634188,  0.0564103,  0.0480586,   0.0383639,  0.027326,  0.0149451, 0.001221},
0823       {-0.0584615, -0.04,     -0.0228763, -0.0070903, 0.00735786, 0.0204682, 0.0322408, 0.0426756, 0.0517726,
0824        0.0595318,  0.0659532, 0.0710368,  0.0747826,  0.0771906,  0.0782609, 0.0779933, 0.076388,  0.0734448,
0825        0.0691639,  0.0635452, 0.0565886,  0.0482943,  0.0386622,  0.0276923, 0.0153846},
0826       {-0.0565385, -0.0387458, -0.022215, -0.00694588, 0.00706142, 0.0198069, 0.0312907, 0.0415126,
0827        0.0504728,  0.0581712,  0.0646078, 0.0697826,   0.0736957,  0.0763469, 0.0777364, 0.0778641,
0828        0.07673,    0.0743341,  0.0706765, 0.0657571,   0.0595759,  0.0521329, 0.0434281, 0.0334615},
0829       {-0.0521739, -0.0359684, -0.0208357, -0.00677583, 0.00621118, 0.0181254, 0.0289667, 0.0387352,
0830        0.0474308,  0.0550536,  0.0616036,  0.0670807,   0.071485,   0.0748165, 0.0770751, 0.0782609,
0831        0.0783738,  0.0774139,  0.0753811,  0.0722756,   0.0680971,  0.0628458, 0.0565217},
0832       {-0.0444664, -0.0311971, -0.0186335, -0.00677583, 0.00437606, 0.0148221, 0.0245624, 0.0335968,
0833        0.0419255,  0.0495483,  0.0564653,  0.0626765,   0.0681818,  0.0729814, 0.0770751, 0.080463,
0834        0.0831451,  0.0851214,  0.0863919,  0.0869565,   0.0868154,  0.0859684},
0835       {-0.0321852, -0.0238283, -0.0155309, -0.00729294, 0.000885613, 0.00900473, 0.0170644,
0836        0.0250646,  0.0330054,  0.0408868,  0.0487087,   0.0564712,   0.0641743,  0.0718179,
0837        0.0794021,  0.0869268,  0.0943921,  0.101798,    0.109144,    0.116431,   0.123659},
0838       {-0.0136364, -0.0130895, -0.0115174, -0.00892003, -0.00529733, -0.000649351, 0.00502392,
0839        0.0117225,  0.0194463,  0.0281955,  0.0379699,   0.0487697,   0.0605947,    0.073445,
0840        0.0873206,  0.102221,   0.118148,   0.135099,    0.153076,    0.172078},
0841       {0.0135338,
0842        0.00200501,
0843        -0.00672269,
0844        -0.0126493,
0845        -0.0157747,
0846        -0.0160991,
0847        -0.0136223,
0848        -0.00834439,
0849        -0.000265369,
0850        0.0106148,
0851        0.024296,
0852        0.0407784,
0853        0.0600619,
0854        0.0821465,
0855        0.107032,
0856        0.134719,
0857        0.165207,
0858        0.198496,
0859        0.234586},
0860       {0.0526316,
0861        0.0227038,
0862        -0.00154799,
0863        -0.0201238,
0864        -0.0330237,
0865        -0.0402477,
0866        -0.0417957,
0867        -0.0376677,
0868        -0.0278638,
0869        -0.0123839,
0870        0.00877193,
0871        0.0356037,
0872        0.0681115,
0873        0.106295,
0874        0.150155,
0875        0.19969,
0876        0.254902,
0877        0.315789},
0878       {0.108359,
0879        0.0505676,
0880        0.00309598,
0881        -0.0340557,
0882        -0.0608875,
0883        -0.0773994,
0884        -0.0835913,
0885        -0.0794634,
0886        -0.0650155,
0887        -0.0402477,
0888        -0.00515996,
0889        0.0402477,
0890        0.0959752,
0891        0.162023,
0892        0.23839,
0893        0.325077,
0894        0.422085}};
0895 
0896   float filtered_baseline[nTotStripsPerAPV];
0897   float filtered_baseline_derivative[nTotStripsPerAPV - 1];
0898 
0899   // Zero filtered_baseline
0900   memset(filtered_baseline, 0, nTotStripsPerAPV * sizeof(float));
0901   // Filter the left edge using (nL, nR) = (0, 16) .. (15, 16) S-G
0902   // filters
0903   for (size_t i = 0; i < savitzky_golay_n_l_r; i++) {
0904     for (size_t j = 0; j < savitzky_golay_n_l_r + 1 + i; j++) {
0905       filtered_baseline[i] += savitzky_golay_coefficient[i][j] * baseline[j];
0906     }
0907   }
0908   // Filter the middle section using the (nL, nR) = (16, 16) S-G
0909   // filter, while taking advantage of the symmetry to save 16
0910   // multiplications.
0911   for (size_t i = savitzky_golay_n_l_r; i < nTotStripsPerAPV - savitzky_golay_n_l_r; i++) {
0912     filtered_baseline[i] = savitzky_golay_coefficient[savitzky_golay_n_l_r][savitzky_golay_n_l_r] * baseline[i];
0913     for (size_t j = 0; j < savitzky_golay_n_l_r; j++) {
0914       filtered_baseline[i] += savitzky_golay_coefficient[savitzky_golay_n_l_r][j] *
0915                               (baseline[i + j - savitzky_golay_n_l_r] + baseline[i - j + savitzky_golay_n_l_r]);
0916     }
0917 #if 0
0918   // Test that the indexing above is correct
0919   float test = 0;
0920   for (size_t j = 0; j < 2 * savitzky_golay_n_l_r + 1; j++) {
0921     test +=
0922         savitzky_golay_coefficient[savitzky_golay_n_l_r][j] *
0923         baseline[i + j - savitzky_golay_n_l_r];
0924   }
0925   // test == filtered_baseline[i] should hold now
0926 #endif
0927   }
0928   // Filter the right edge using (nL, nR) = (16, 15) .. (16, 0) S-G
0929   // filters
0930   for (size_t i = nTotStripsPerAPV - savitzky_golay_n_l_r; i < nTotStripsPerAPV; i++) {
0931     for (size_t j = 0; j < nTotStripsPerAPV - i + savitzky_golay_n_l_r; j++) {
0932       filtered_baseline[i] += savitzky_golay_coefficient[2 * savitzky_golay_n_l_r + i + 1 - nTotStripsPerAPV][j] *
0933                               baseline[i + j - savitzky_golay_n_l_r];
0934     }
0935   }
0936   // In lieu of a spearate S-G derivative filter, the finite
0937   // difference is used here (since the output is sufficiently
0938   // smooth).
0939   for (size_t i = 0; i < (nTotStripsPerAPV - 1); i++) {
0940     filtered_baseline_derivative[i] = filtered_baseline[i + 1] - filtered_baseline[i];
0941   }
0942 
0943   // Calculate the maximum deviation between filtered and unfiltered
0944   // baseline, plus the sum square of the derivative.
0945 
0946   double filtered_baseline_max = 0;
0947   double filtered_baseline_derivative_sum_square = 0;
0948 
0949   for (size_t i = 0; i < nTotStripsPerAPV; i++) {
0950     const double d = filtered_baseline[i] - baseline[i];
0951 
0952     filtered_baseline_max = std::max(filtered_baseline_max, static_cast<double>(fabs(d)));
0953   }
0954   for (size_t i = 0; i < (nTotStripsPerAPV - 1); i++) {
0955     filtered_baseline_derivative_sum_square += filtered_baseline_derivative[i] * filtered_baseline_derivative[i];
0956   }
0957 
0958 #if 0
0959   std::cerr << __FILE__ << ':' << __LINE__ << ": "
0960           << filtered_baseline_max << ' '
0961           << filtered_baseline_derivative_sum_square << std::endl;
0962 #endif
0963 
0964   // Apply the cut
0965   return !(filtered_baseline_max >= filteredBaselineMax_ ||
0966            filtered_baseline_derivative_sum_square >= filteredBaselineDerivativeSumSquare_);
0967 }
0968 
0969 //Other methods implementation ==============================================
0970 //==========================================================================
0971 
0972 void SiStripAPVRestorer::loadMeanCMMap(const edm::Event& iEvent) {
0973   if (useRealMeanCM_) {
0974     edm::Handle<edm::DetSetVector<SiStripRawDigi>> input;
0975     iEvent.getByToken(siStripRawDigiToken_, input);
0976     createCMMapRealPed(*input);
0977   } else {
0978     edm::Handle<edm::DetSetVector<SiStripProcessedRawDigi>> inputCM;
0979     iEvent.getByToken(siStripProcessedRawDigiToken_, inputCM);
0980     createCMMapCMstored(*inputCM);
0981   }
0982 }
0983 void SiStripAPVRestorer::createCMMapRealPed(const edm::DetSetVector<SiStripRawDigi>& input) {
0984   meanCMmap_.clear();
0985   for (const auto& rawDigis : input) {
0986     const auto detPedestalRange = pedestalHandle->getRange(rawDigis.id);
0987     const size_t nAPVs = rawDigis.size() / nTotStripsPerAPV;
0988     std::vector<float> meanCMDetSet;
0989     meanCMDetSet.reserve(nAPVs);
0990     for (uint16_t iAPV = 0; iAPV < nAPVs; ++iAPV) {
0991       uint16_t minPed = nTotStripsPerAPV;
0992       for (uint16_t strip = nTotStripsPerAPV * iAPV; strip < nTotStripsPerAPV * (iAPV + 1); ++strip) {
0993         uint16_t ped = static_cast<uint16_t>(pedestalHandle->getPed(strip, detPedestalRange));
0994         if (ped < minPed)
0995           minPed = ped;
0996       }
0997       meanCMDetSet.push_back(minPed);
0998     }
0999     meanCMmap_.emplace(rawDigis.id, std::move(meanCMDetSet));
1000   }
1001 }
1002 void SiStripAPVRestorer::createCMMapCMstored(const edm::DetSetVector<SiStripProcessedRawDigi>& input) {
1003   meanCMmap_.clear();
1004   for (const auto& rawDigis : input) {
1005     std::vector<float> meanCMNValue;
1006     meanCMNValue.reserve(rawDigis.size());
1007     std::transform(
1008         std::begin(rawDigis), std::end(rawDigis), std::back_inserter(meanCMNValue), [](SiStripProcessedRawDigi cm) {
1009           return cm.adc();
1010         });
1011     meanCMmap_.emplace(rawDigis.id, std::move(meanCMNValue));
1012   }
1013 }
1014 
1015 //NEW algorithm designed for implementation in FW=============================
1016 //============================================================================
1017 void SiStripAPVRestorer::derivativeFollowerRestore(uint16_t apvN, uint16_t firstAPV, digivector_t& digis) {
1018   digivector_t singleAPVdigi;
1019   singleAPVdigi.clear();
1020   for (int16_t strip = (apvN - firstAPV) * nTotStripsPerAPV; strip < (apvN - firstAPV + 1) * nTotStripsPerAPV; ++strip)
1021     singleAPVdigi.push_back(digis[strip] + 1024);
1022 
1023   digimap_t discontinuities;  // it will contain the start and the end of each region in which a greadient is present
1024   discontinuities.clear();
1025 
1026   digimap_t::iterator itdiscontinuities;
1027 
1028   //----Variables of the first part---//
1029 
1030   bool isMinimumAndNoMax = false;
1031   bool isFirstStrip = false;
1032 
1033   //---Variables of the second part---//
1034 
1035   int actualStripADC = 0;
1036   int previousStripADC = 0;
1037 
1038   int greadient = 0;
1039   int maximum_value = 0;
1040   const uint32_t n_high_maximum_cluster = 1025 + 1024;
1041   int high_maximum_cluster = n_high_maximum_cluster;
1042   int number_good_minimum = 0;
1043   int first_gradient = 0;
1044   int strip_first_gradient = 0;
1045   int adc_start_point_cluster_pw = 0;
1046   int auxiliary_end_cluster = 0;
1047   int first_start_cluster_strip = 0;
1048   int first_start_cluster_ADC = 0;
1049   bool isAuxiliary_Minimum = false;
1050   bool isPossible_wrong_minimum = false;
1051   bool isMax = false;
1052 
1053   //----------SECOND PART: CLUSTER FINDING--------//
1054 
1055   for (uint16_t strip = 0; strip < singleAPVdigi.size(); ++strip) {
1056     if (strip == 0) {
1057       actualStripADC = singleAPVdigi[strip];
1058       if (std::abs(singleAPVdigi[strip] - singleAPVdigi[strip + 1]) > gradient_threshold_) {
1059         isFirstStrip = true;
1060         isMinimumAndNoMax = true;
1061         discontinuities.insert(discontinuities.end(), std::pair<int, int>(strip, actualStripADC));
1062       } else if (actualStripADC > (gradient_threshold_ + 1024)) {
1063         discontinuities.insert(discontinuities.end(), std::pair<uint16_t, int16_t>(strip, 0 + 1024));
1064         isMinimumAndNoMax = true;
1065         first_start_cluster_strip = strip;
1066         first_start_cluster_ADC = 1024;
1067       }
1068 
1069     } else {
1070       previousStripADC = actualStripADC;
1071       actualStripADC = singleAPVdigi[strip];
1072       greadient = actualStripADC - previousStripADC;
1073 
1074       if (((greadient > gradient_threshold_) && (!isMax) && (!isMinimumAndNoMax))) {
1075         isMax = false;
1076         if (((std::abs(maximum_value - previousStripADC) < (2 * gradient_threshold_)) &&
1077              discontinuities.size() > 1)) {  // add the || with the cluster size
1078           //to verify if the noise do not interfere and is detected fake hits
1079           isPossible_wrong_minimum = true;
1080           itdiscontinuities = --discontinuities.end();
1081           if (discontinuities.size() > 1) {
1082             --itdiscontinuities;
1083           }
1084           strip_first_gradient = itdiscontinuities->first;
1085           adc_start_point_cluster_pw = itdiscontinuities->second;
1086           first_gradient = std::abs(adc_start_point_cluster_pw - singleAPVdigi[strip_first_gradient + 1]);
1087           discontinuities.erase(++itdiscontinuities);
1088           discontinuities.erase(--discontinuities.end());
1089         }
1090 
1091         if ((discontinuities.size() % 2 == 1) && (!discontinuities.empty())) {  //&&(no_minimo == 0)
1092           discontinuities.erase(--discontinuities.end());
1093         }
1094 
1095         discontinuities.insert(discontinuities.end(), std::pair<uint16_t, int16_t>(strip - 1, previousStripADC));
1096         isMinimumAndNoMax = true;
1097         maximum_value = 0;
1098         first_start_cluster_strip = strip - 1;
1099         first_start_cluster_ADC = previousStripADC;
1100 
1101       } else if ((!isMax) && ((actualStripADC - previousStripADC < 0) && (isMinimumAndNoMax))) {
1102         isMax = true;
1103         isMinimumAndNoMax = false;
1104         high_maximum_cluster = n_high_maximum_cluster;
1105         if ((previousStripADC > maximum_value) && (discontinuities.size() % 2 == 1))
1106           maximum_value = previousStripADC;
1107       }
1108 
1109       if ((isMax) && (strip < (nTotStripsPerAPV - 2))) {
1110         if (high_maximum_cluster > (std::abs(singleAPVdigi[strip + 1] - actualStripADC))) {
1111           high_maximum_cluster = singleAPVdigi[strip + 1] - actualStripADC;
1112           auxiliary_end_cluster = strip + 2;
1113         } else {
1114           auxiliary_end_cluster = nTotStripsPerAPV - 1;
1115         }
1116       }
1117 
1118       if ((isMax) && ((actualStripADC - previousStripADC) >= 0) && (size_window_ > 0) &&
1119           (strip <= ((nTotStripsPerAPV - 1) - size_window_ - 1))) {
1120         number_good_minimum = 0;
1121         for (uint16_t wintry = 0; wintry <= size_window_; wintry++) {
1122           if (std::abs(singleAPVdigi[strip + wintry] - singleAPVdigi[strip + wintry + 1]) <= last_gradient_)
1123             ++number_good_minimum;
1124         }
1125         --number_good_minimum;
1126 
1127         if (size_window_ != number_good_minimum) {
1128           isMax = false;  // not Valid end Point
1129           isMinimumAndNoMax = true;
1130         }
1131 
1132       } else if ((isMax) && (strip > ((nTotStripsPerAPV - 1) - size_window_ - 1))) {
1133         isMax = true;
1134         isAuxiliary_Minimum = true;  //for minimums after strip 127-SW-1
1135       }
1136 
1137       if (!discontinuities.empty()) {
1138         itdiscontinuities = --discontinuities.end();
1139       }
1140 
1141       if ((isMax) && (actualStripADC <= first_start_cluster_ADC)) {
1142         if ((std::abs(first_start_cluster_ADC - singleAPVdigi[strip + 2]) > first_gradient) &&
1143             (isPossible_wrong_minimum)) {
1144           discontinuities.erase(itdiscontinuities);
1145           discontinuities.insert(discontinuities.end(),
1146                                  std::pair<int, int>(strip_first_gradient, adc_start_point_cluster_pw));
1147           discontinuities.insert(discontinuities.end(), std::pair<int, int>(strip, adc_start_point_cluster_pw));
1148         } else {
1149           if ((discontinuities.size() % 2 == 0) && (discontinuities.size() > 1)) {
1150             discontinuities.erase(--discontinuities.end());
1151           }
1152           discontinuities.insert(discontinuities.end(), std::pair<int, int>(strip, actualStripADC));
1153         }
1154         isMax = false;
1155         adc_start_point_cluster_pw = 0;
1156         isPossible_wrong_minimum = false;
1157         strip_first_gradient = 0;
1158         first_start_cluster_strip = 0;
1159         first_start_cluster_ADC = 0;
1160       }
1161 
1162       if ((isMax) && ((actualStripADC - previousStripADC) >= 0) &&
1163           (!isAuxiliary_Minimum)) {  //For the end Poit when strip >127-SW-1
1164         if ((std::abs(first_start_cluster_ADC - singleAPVdigi[strip + 1]) > first_gradient) &&
1165             (isPossible_wrong_minimum)) {
1166           discontinuities.erase(itdiscontinuities);
1167           discontinuities.insert(discontinuities.end(),
1168                                  std::pair<int, int>(strip_first_gradient, adc_start_point_cluster_pw));
1169         }
1170 
1171         if ((discontinuities.size() % 2 == 0) && (discontinuities.size() > 1)) {
1172           discontinuities.erase(--discontinuities.end());
1173         }
1174         discontinuities.insert(discontinuities.end(), std::pair<int, int>(strip - 1, previousStripADC));
1175         isMax = false;
1176         adc_start_point_cluster_pw = 0;
1177         isPossible_wrong_minimum = false;
1178         strip_first_gradient = 0;
1179         first_start_cluster_strip = 0;
1180         first_start_cluster_ADC = 0;
1181       }
1182     }
1183   }
1184 
1185   //----------THIRD PART reconstruction of the event without baseline-------//
1186 
1187   if (!discontinuities.empty()) {
1188     if ((first_start_cluster_strip) == (nTotStripsPerAPV - 1) - 1)
1189       discontinuities.insert(discontinuities.end(),
1190                              std::pair<int, int>((nTotStripsPerAPV - 1), first_start_cluster_ADC));
1191     if ((isMax) && (isAuxiliary_Minimum))
1192       discontinuities.insert(discontinuities.end(),
1193                              std::pair<int, int>(auxiliary_end_cluster, first_start_cluster_ADC));
1194   }
1195 
1196   if (isFirstStrip) {
1197     itdiscontinuities = discontinuities.begin();
1198     ++itdiscontinuities;
1199     ++itdiscontinuities;
1200     int firstADC = itdiscontinuities->second;
1201     --itdiscontinuities;
1202     itdiscontinuities->second = firstADC;
1203     --itdiscontinuities;
1204     itdiscontinuities->second = firstADC;
1205   }
1206 
1207   if (!discontinuities.empty()) {
1208     itdiscontinuities = discontinuities.begin();
1209     uint16_t firstStrip = itdiscontinuities->first;
1210     int16_t firstADC = itdiscontinuities->second;
1211     ++itdiscontinuities;
1212     uint16_t lastStrip = itdiscontinuities->first;
1213     int16_t secondADC = itdiscontinuities->second;
1214     ++itdiscontinuities;
1215 
1216     for (uint16_t strip = 0; strip < nTotStripsPerAPV; ++strip) {
1217       if (strip > lastStrip && itdiscontinuities != discontinuities.end()) {
1218         firstStrip = itdiscontinuities->first;
1219         firstADC = itdiscontinuities->second;
1220         ++itdiscontinuities;
1221         lastStrip = itdiscontinuities->first;
1222         secondADC = itdiscontinuities->second;
1223         ++itdiscontinuities;
1224       }
1225 
1226       if ((firstStrip <= strip) && (strip <= lastStrip) &&
1227           (0 < (singleAPVdigi[strip] - firstADC -
1228                 ((secondADC - firstADC) / (lastStrip - firstStrip)) * (strip - firstStrip) - gradient_threshold_))) {
1229         digis[(apvN - firstAPV) * nTotStripsPerAPV + strip] =
1230             singleAPVdigi[strip] - firstADC -
1231             (((secondADC - firstADC) / (lastStrip - firstStrip)) * (strip - firstStrip));
1232         edm::LogWarning("NoBaseline") << "No baseline " << digis[(apvN - firstAPV) * nTotStripsPerAPV + strip] << "\n";
1233       } else {
1234         digis[(apvN - firstAPV) * nTotStripsPerAPV + strip] = 0;
1235       }
1236     }
1237   }
1238 }