Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 10:00:46 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!