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
0022 useRealMeanCM_(conf.getParameter<bool>("useRealMeanCM")),
0023 meanCM_(conf.getParameter<int32_t>("MeanCM")),
0024 deltaCMThreshold_(conf.getParameter<uint32_t>("DeltaCMThreshold")),
0025
0026 fraction_(conf.getParameter<double>("Fraction")),
0027 deviation_(conf.getParameter<uint32_t>("Deviation")),
0028
0029 restoreThreshold_(conf.getParameter<double>("restoreThreshold")),
0030
0031 nSaturatedStrip_(conf.getParameter<uint32_t>("nSaturatedStrip")),
0032
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
0040 slopeX_(conf.getParameter<int32_t>("slopeX")),
0041 slopeY_(conf.getParameter<int32_t>("slopeY")),
0042
0043 hitStripThreshold_(conf.getParameter<uint32_t>("hitStripThreshold")),
0044
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
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
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
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_;
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_;
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_;
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_;
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_;
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_;
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
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
0350 if (!smoothedMaps_.empty()) {
0351 auto itSmootedMap = smoothedMaps_.find(apvN);
0352 baselineFollower(itSmootedMap->second, baseline, median);
0353 } else {
0354
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
0370 for (int16_t itStrip = 0; itStrip < nTotStripsPerAPV; ++itStrip) {
0371 digis[nTotStripsPerAPV * (apvN - firstAPV) + itStrip] -= baseline[itStrip] - median;
0372 }
0373
0374
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
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
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
0410 digimap_t consecpoints;
0411 std::vector<uint16_t> nConsStrip;
0412
0413 uint16_t consecStrips = 0;
0414 for (uint32_t istrip = 0; istrip < nTotStripsPerAPV; ++istrip) {
0415 const int16_t adc = adcs[istrip];
0416
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
0427 if (consecStrips > 0)
0428 nConsStrip.push_back(consecStrips);
0429
0430
0431 auto itConsecpoints = consecpoints.begin();
0432 float MinSmoothValue = 20000., MaxSmoothValue = 0.;
0433 for (const auto consecStrips : nConsStrip) {
0434 if (consecStrips >= consecThreshold_) {
0435 ++itConsecpoints;
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;
0459
0460 if (stripCount > 1) {
0461
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) {
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
0507
0508 if (smoothedpoints.size() < 3)
0509 return;
0510
0511 auto it = smoothedpoints.begin();
0512 while (smoothedpoints.size() > 3 && it != --(--smoothedpoints.end())) {
0513
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
0532
0533
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
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
0596
0597 if (smoothedpoints.size() < 4)
0598 return;
0599
0600 auto it = smoothedpoints.begin();
0601 while (smoothedpoints.size() > 2 && it != --smoothedpoints.end()) {
0602
0603
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) {
0613 smoothedpoints.erase(itNext);
0614 } else if (m < -2) {
0615
0616 if (it == smoothedpoints.begin())
0617 smoothedpoints.erase(it++);
0618
0619 else
0620 smoothedpoints.erase(it--);
0621 } else {
0622 ++it;
0623 }
0624 }
0625 }
0626
0627 void inline SiStripAPVRestorer::baselineFollower(const digimap_t& smoothedpoints,
0628 digivector_t& baseline,
0629 float median) {
0630
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
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
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
0675
0676
0677
0678
0679
0680
0681
0682
0683
0684
0685
0686
0687
0688
0689
0690
0691
0692
0693
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
0901 memset(filtered_baseline, 0, nTotStripsPerAPV * sizeof(float));
0902
0903
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
0910
0911
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
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
0927 #endif
0928 }
0929
0930
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
0938
0939
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
0945
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
0966 return !(filtered_baseline_max >= filteredBaselineMax_ ||
0967 filtered_baseline_derivative_sum_square >= filteredBaselineDerivativeSumSquare_);
0968 }
0969
0970
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
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;
1025 discontinuities.clear();
1026
1027 digimap_t::iterator itdiscontinuities;
1028
1029
1030
1031 bool isMinimumAndNoMax = false;
1032 bool isFirstStrip = false;
1033
1034
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
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)) {
1079
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())) {
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;
1130 isMinimumAndNoMax = true;
1131 }
1132
1133 } else if ((isMax) && (strip > ((nTotStripsPerAPV - 1) - size_window_ - 1))) {
1134 isMax = true;
1135 isAuxiliary_Minimum = true;
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)) {
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
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 }