Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:26:31 

0001 #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripFedZeroSuppression.h"
0002 
0003 #include "CondFormats/DataRecord/interface/SiStripNoisesRcd.h"
0004 #include "CondFormats/SiStripObjects/interface/SiStripNoises.h"
0005 #include "CondFormats/DataRecord/interface/SiStripThresholdRcd.h"
0006 #include "CondFormats/SiStripObjects/interface/SiStripThreshold.h"
0007 
0008 //#define DEBUG_SiStripZeroSuppression_
0009 //#define ML_DEBUG
0010 using namespace std;
0011 
0012 void SiStripFedZeroSuppression::init(const edm::EventSetup& es) {
0013   if (noiseWatcher_.check(es)) {
0014     noise_ = &es.getData(noiseToken_);
0015   }
0016   if (thresholdWatcher_.check(es)) {
0017     threshold_ = &es.getData(thresholdToken_);
0018   }
0019 }
0020 
0021 namespace {
0022 
0023   constexpr uint8_t zeroTh = 0;
0024   constexpr uint8_t lowTh = 1;
0025   constexpr uint8_t highTh = 2;
0026 
0027   struct Payload {
0028     uint8_t stat;
0029     uint8_t statPrev;
0030     uint8_t statNext;
0031     uint8_t statMaxNeigh;
0032     uint8_t statPrev2;
0033     uint8_t statNext2;
0034   };
0035 
0036   constexpr bool isDigiValid(Payload const& data, uint16_t FEDalgorithm) {
0037     // Decide if this strip should be accepted.
0038     bool accept = false;
0039     switch (FEDalgorithm) {
0040       case 1:
0041         accept = (data.stat >= lowTh);
0042         break;
0043       case 2:
0044         accept = (data.stat >= highTh || (data.stat >= lowTh && data.statMaxNeigh >= lowTh));
0045         break;
0046       case 3:
0047         accept = (data.stat >= highTh || (data.stat >= lowTh && data.statMaxNeigh >= highTh));
0048         break;
0049       case 4:
0050         accept = ((data.stat >= highTh)     //Test for adc>highThresh (same as algorithm 2)
0051                   || ((data.stat >= lowTh)  //Test for adc>lowThresh, with neighbour adc>lowThresh (same as algorithm 2)
0052                       && (data.statMaxNeigh >= lowTh)) ||
0053                   ((data.stat < lowTh)             //Test for adc<lowThresh
0054                    && (((data.statPrev >= highTh)  //with both neighbours>highThresh
0055                         && (data.statNext >= highTh)) ||
0056                        ((data.statPrev >= highTh)    //OR with previous neighbour>highThresh and
0057                         && (data.statNext >= lowTh)  //both the next neighbours>lowThresh
0058                         && (data.statNext2 >= lowTh)) ||
0059                        ((data.statNext >= highTh)    //OR with next neighbour>highThresh and
0060                         && (data.statPrev >= lowTh)  //both the previous neighbours>lowThresh
0061                         && (data.statPrev2 >= lowTh)) ||
0062                        ((data.statNext >= lowTh)      //OR with both next neighbours>lowThresh and
0063                         && (data.statNext2 >= lowTh)  //both the previous neighbours>lowThresh
0064                         && (data.statPrev >= lowTh) && (data.statPrev2 >= lowTh)))));
0065         break;
0066       case 5:
0067         accept = true;  // zero removed in conversion
0068         break;
0069     }
0070     return accept;
0071   }
0072 
0073 }  // namespace
0074 
0075 void SiStripFedZeroSuppression::suppress(const std::vector<SiStripDigi>& in,
0076                                          std::vector<SiStripDigi>& selectedSignal,
0077                                          uint32_t detID) {
0078   suppress(in, selectedSignal, detID, *noise_, *threshold_);
0079 }
0080 
0081 void SiStripFedZeroSuppression::suppress(const std::vector<SiStripDigi>& in,
0082                                          std::vector<SiStripDigi>& selectedSignal,
0083                                          uint32_t detID,
0084                                          const SiStripNoises& noise,
0085                                          const SiStripThreshold& threshold) {
0086   selectedSignal.clear();
0087   size_t inSize = in.size();
0088   if (inSize == 0) {
0089     return;
0090   }
0091 
0092   SiStripNoises::Range detNoiseRange = noise.getRange(detID);
0093   SiStripThreshold::Range detThRange = threshold.getRange(detID);
0094 
0095   // reserving more than needed, but quicker than one at a time
0096   selectedSignal.reserve(inSize);
0097 
0098   // load status
0099   uint8_t stat[inSize];
0100   for (size_t i = 0; i < inSize; i++) {
0101     auto strip = (uint32_t)in[i].strip();
0102 
0103     auto ladc = in[i].adc();
0104     assert(ladc > 0);
0105 
0106     auto thresholds = threshold.getData(strip, detThRange);
0107 
0108     auto highThresh = static_cast<int16_t>(thresholds.getHth() * noise.getNoiseFast(strip, detNoiseRange) + 0.5f);
0109     auto lowThresh = static_cast<int16_t>(thresholds.getLth() * noise.getNoiseFast(strip, detNoiseRange) + 0.5f);
0110 
0111     assert(lowThresh >= 0);
0112     assert(lowThresh <= highThresh);
0113 
0114     stat[i] = zeroTh;
0115     if (ladc >= lowThresh)
0116       stat[i] = lowTh;
0117     if (ladc >= highThresh)
0118       stat[i] = highTh;
0119   }
0120 
0121   for (size_t i = 0; i < inSize; i++) {
0122     auto strip = (uint32_t)in[i].strip();
0123     Payload ldata;
0124 
0125     ldata.stat = stat[i];
0126     ldata.statPrev = zeroTh;
0127     ldata.statNext = zeroTh;
0128     ldata.statPrev2 = zeroTh;
0129     ldata.statNext2 = zeroTh;
0130 
0131     if (((strip) % 128) == 127) {
0132       ldata.statNext = zeroTh;
0133     } else if (i + 1 < inSize && in[i + 1].strip() == strip + 1) {
0134       ldata.statNext = stat[i + 1];
0135       if (((strip) % 128) == 126) {
0136         ldata.statNext2 = zeroTh;
0137       } else if (i + 2 < inSize && in[i + 2].strip() == strip + 2) {
0138         ldata.statNext2 = stat[i + 2];
0139       }
0140     }
0141 
0142     if (((strip) % 128) == 0) {
0143       ldata.statPrev = zeroTh;
0144     } else if (i >= 1 && in[i - 1].strip() == strip - 1) {
0145       ldata.statPrev = stat[i - 1];
0146       if (((strip) % 128) == 1) {
0147         ldata.statPrev2 = zeroTh;
0148       } else if (i >= 2 && in[i - 2].strip() == strip - 2) {
0149         ldata.statPrev2 = stat[i - 2];
0150       }
0151     }
0152 
0153     ldata.statMaxNeigh = std::max(ldata.statPrev, ldata.statNext);
0154 
0155     if (isDigiValid(ldata, theFEDalgorithm)) {
0156       selectedSignal.push_back(SiStripDigi(strip, in[i].adc()));
0157     }
0158   }
0159 }
0160 
0161 void SiStripFedZeroSuppression::suppress(const edm::DetSet<SiStripRawDigi>& in, edm::DetSet<SiStripDigi>& out) {
0162   const uint32_t detID = out.id;
0163   SiStripNoises::Range detNoiseRange = noise_->getRange(detID);
0164   SiStripThreshold::Range detThRange = threshold_->getRange(detID);
0165 #ifdef DEBUG_SiStripZeroSuppression_
0166   if (edm::isDebugEnabled())
0167     LogTrace("SiStripZeroSuppression")
0168         << "[SiStripFedZeroSuppression::suppress] Zero suppression on edm::DetSet<SiStripRawDigi>: detID " << detID
0169         << " size = " << in.data.size();
0170 #endif
0171   edm::DetSet<SiStripRawDigi>::const_iterator in_iter = in.data.begin();
0172   for (; in_iter != in.data.end(); in_iter++) {
0173     const uint32_t strip = (uint32_t)(in_iter - in.data.begin());
0174 
0175 #ifdef DEBUG_SiStripZeroSuppression_
0176     if (edm::isDebugEnabled())
0177       LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress] detID= " << detID
0178                                          << " strip= " << strip << "  adc= " << in_iter->adc();
0179 #endif
0180     adc = in_iter->adc();
0181 
0182     SiStripThreshold::Data thresholds = threshold_->getData(strip, detThRange);
0183     theFEDlowThresh = static_cast<int16_t>(thresholds.getLth() * noise_->getNoiseFast(strip, detNoiseRange) + 0.5);
0184     theFEDhighThresh = static_cast<int16_t>(thresholds.getHth() * noise_->getNoiseFast(strip, detNoiseRange) + 0.5);
0185 
0186     adcPrev = -9999;
0187     adcNext = -9999;
0188     /*
0189       If a strip is the last one on the chip
0190       set its next neighbor's thresholds to infinity
0191       because the FED does not merge clusters across
0192       chip boundaries right now
0193     */
0194     if (strip % 128 == 127) {
0195       adcNext = 0;
0196       theNextFEDlowThresh = 9999;
0197       theNextFEDhighThresh = 9999;
0198     } else {
0199       adcNext = (in_iter + 1)->adc();
0200       SiStripThreshold::Data thresholds_1 = threshold_->getData(strip + 1, detThRange);
0201       theNextFEDlowThresh =
0202           static_cast<int16_t>(thresholds_1.getLth() * noise_->getNoiseFast(strip + 1, detNoiseRange) + 0.5);
0203       theNextFEDhighThresh =
0204           static_cast<int16_t>(thresholds_1.getHth() * noise_->getNoiseFast(strip + 1, detNoiseRange) + 0.5);
0205     }
0206     /*
0207       Similarily, for the first strip 
0208       on a chip
0209     */
0210     if (strip % 128 == 0) {
0211       adcPrev = 0;
0212       thePrevFEDlowThresh = 9999;
0213       thePrevFEDhighThresh = 9999;
0214     } else {
0215       adcPrev = (in_iter - 1)->adc();
0216       SiStripThreshold::Data thresholds_1 = threshold_->getData(strip - 1, detThRange);
0217       thePrevFEDlowThresh =
0218           static_cast<int16_t>(thresholds_1.getLth() * noise_->getNoiseFast(strip - 1, detNoiseRange) + 0.5);
0219       thePrevFEDhighThresh =
0220           static_cast<int16_t>(thresholds_1.getHth() * noise_->getNoiseFast(strip - 1, detNoiseRange) + 0.5);
0221     }
0222     if (adcNext < adcPrev) {
0223       adcMaxNeigh = adcPrev;
0224       theNeighFEDlowThresh = thePrevFEDlowThresh;
0225       theNeighFEDhighThresh = thePrevFEDhighThresh;
0226     } else {
0227       adcMaxNeigh = adcNext;
0228       theNeighFEDlowThresh = theNextFEDlowThresh;
0229       theNeighFEDhighThresh = theNextFEDhighThresh;
0230     }
0231 
0232     //Find adc values for next neighbouring strips
0233     adcPrev2 = -9999;
0234     adcNext2 = -9999;
0235     thePrev2FEDlowThresh = 1;
0236     theNext2FEDlowThresh = 1;
0237     if (strip % 128 >= 126) {
0238       adcNext2 = 0;
0239       theNext2FEDlowThresh = 9999;
0240     } else if (strip % 128 < 126) {
0241       adcNext2 = (in_iter + 2)->adc();
0242       theNext2FEDlowThresh = static_cast<int16_t>(
0243           threshold_->getData(strip + 2, detThRange).getLth() * noise_->getNoiseFast(strip + 2, detNoiseRange) + 0.5);
0244     }
0245     if (strip % 128 <= 1) {
0246       adcPrev2 = 0;
0247       thePrev2FEDlowThresh = 9999;
0248     } else if (strip % 128 > 1) {
0249       adcPrev2 = (in_iter - 2)->adc();
0250       thePrev2FEDlowThresh = static_cast<int16_t>(
0251           threshold_->getData(strip - 2, detThRange).getLth() * noise_->getNoiseFast(strip - 2, detNoiseRange) + 0.5);
0252     }
0253     //GB 23/6/08: truncation should be done at the very beginning
0254     if (isAValidDigi())
0255       out.data.push_back(SiStripDigi(strip, truncate(in_iter->adc())));
0256   }
0257 }
0258 
0259 void SiStripFedZeroSuppression::fillThresholds_(const uint32_t detID, size_t size) {
0260   SiStripNoises::Range detNoiseRange = noise_->getRange(detID);
0261   SiStripThreshold::Range detThRange = threshold_->getRange(detID);
0262 
0263   if (highThr_.size() != size) {
0264     highThr_.resize(size);
0265     lowThr_.resize(size);
0266     noises_.resize(size);
0267     highThrSN_.resize(size);
0268     lowThrSN_.resize(size);
0269   }
0270 
0271   noise_->allNoises(noises_, detNoiseRange);
0272   threshold_->allThresholds(lowThrSN_, highThrSN_, detThRange);  // thresholds as S/N
0273   for (size_t strip = 0; strip < size; ++strip) {
0274     float noise = noises_[strip];
0275     //  uncomment line below to check bluk noise decoding
0276     //assert( noise == noiseHandle->getNoiseFast(strip,detNoiseRange) );
0277     highThr_[strip] = static_cast<int16_t>(highThrSN_[strip] * noise + 0.5 + 1e-6);
0278     lowThr_[strip] = static_cast<int16_t>(lowThrSN_[strip] * noise + 0.5 + 1e-6);
0279     // Note: it's a bit wierd, but there are some cases for which 'highThrSN_[strip]*noise' is an exact integer
0280     //   but due to roundoffs it gets rounded to the integer below if.
0281     //   Apparently the optimized code inlines differently and this changes the roundoff.
0282     //   The +1e-6 fixes the problem.   [GPetruc]
0283   }
0284 }
0285 
0286 void SiStripFedZeroSuppression::suppress(const std::vector<int16_t>& in,
0287                                          uint16_t firstAPV,
0288                                          edm::DetSet<SiStripDigi>& out) {
0289   const uint32_t detID = out.id;
0290   size_t size = in.size();
0291 #ifdef DEBUG_SiStripZeroSuppression_
0292   if (edm::isDebugEnabled())
0293     LogTrace("SiStripZeroSuppression")
0294         << "[SiStripFedZeroSuppression::suppress] Zero suppression on std::vector<int16_t>: detID " << detID
0295         << " size = " << in.size();
0296 #endif
0297 
0298   fillThresholds_(detID, size + firstAPV * 128);  // want to decouple this from the other cost
0299 
0300   std::vector<int16_t>::const_iterator in_iter = in.begin();
0301   uint16_t strip = firstAPV * 128;
0302   for (; strip < size + firstAPV * 128; ++strip, ++in_iter) {
0303     size_t strip_mod_128 = strip & 127;
0304 #ifdef DEBUG_SiStripZeroSuppression_
0305     if (edm::isDebugEnabled())
0306       LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress]  detID= " << detID
0307                                          << " strip= " << strip << "  adc= " << *in_iter;
0308 #endif
0309     adc = *in_iter;
0310 
0311     theFEDlowThresh = lowThr_[strip];
0312     theFEDhighThresh = highThr_[strip];
0313 
0314     //Find adc values for neighbouring strips
0315 
0316     /*
0317       If a strip is the last one on the chip
0318       set its next neighbor's thresholds to infinity
0319       because the FED does not merge clusters across
0320       chip boundaries right now
0321     */
0322 
0323     //adcPrev = -9999;  // useless, they are set
0324     //adcNext = -9999;  // in the next lines in any case
0325     if (strip_mod_128 == 127) {
0326       adcNext = 0;
0327       theNextFEDlowThresh = 9999;
0328       theNextFEDhighThresh = 9999;
0329     } else {
0330       adcNext = *(in_iter + 1);
0331       theNextFEDlowThresh = lowThr_[strip + 1];
0332       theNextFEDhighThresh = highThr_[strip + 1];
0333     }
0334 
0335     /*
0336       Similarily, for the first strip 
0337       on a chip
0338     */
0339     if (strip_mod_128 == 0) {
0340       adcPrev = 0;
0341       thePrevFEDlowThresh = 9999;
0342       thePrevFEDhighThresh = 9999;
0343     } else {
0344       adcPrev = *(in_iter - 1);
0345       thePrevFEDlowThresh = lowThr_[strip - 1];
0346       thePrevFEDhighThresh = highThr_[strip - 1];
0347     }
0348 
0349     if (adcNext < adcPrev) {
0350       adcMaxNeigh = adcPrev;
0351       theNeighFEDlowThresh = thePrevFEDlowThresh;
0352       theNeighFEDhighThresh = thePrevFEDhighThresh;
0353     } else {
0354       adcMaxNeigh = adcNext;
0355       theNeighFEDlowThresh = theNextFEDlowThresh;
0356       theNeighFEDhighThresh = theNextFEDhighThresh;
0357     }
0358 
0359     //Find adc values for next neighbouring strips
0360     //adcPrev2 = -9999;           //
0361     //adcNext2 = -9999;           // useless to set them here
0362     //thePrev2FEDlowThresh  = 1;  // they are overwritten always in the next 8 lines
0363     //theNext2FEDlowThresh  = 1;  //
0364     if (strip_mod_128 >= 126) {
0365       adcNext2 = 0;
0366       theNext2FEDlowThresh = 9999;
0367       //} else if ( strip_mod_128 < 126 ) { // if it's not >= then is <, no need to "if" again
0368     } else {
0369       adcNext2 = *(in_iter + 2);
0370       theNext2FEDlowThresh = lowThr_[strip + 2];
0371     }
0372     if (strip_mod_128 <= 1) {
0373       adcPrev2 = 0;
0374       thePrev2FEDlowThresh = 9999;
0375       //} else if ( strip_mod_128 > 1 ) { // same as above
0376     } else {
0377       adcPrev2 = *(in_iter - 2);
0378       thePrev2FEDlowThresh = lowThr_[strip - 2];
0379       ;
0380     }
0381 
0382     if (isAValidDigi()) {
0383 #ifdef DEBUG_SiStripZeroSuppression_
0384       if (edm::isDebugEnabled())
0385         LogTrace("SiStripZeroSuppression")
0386             << "[SiStripFedZeroSuppression::suppress] DetId " << out.id << " strip " << strip << " adc " << *in_iter
0387             << " digiCollection size " << out.data.size();
0388 #endif
0389       //GB 23/6/08: truncation should be done at the very beginning
0390       out.push_back(SiStripDigi(strip, (*in_iter < 0 ? 0 : truncate(*in_iter))));
0391     }
0392   }
0393 }
0394 
0395 bool SiStripFedZeroSuppression::isAValidDigi() {
0396 #ifdef DEBUG_SiStripZeroSuppression_
0397 
0398   if (edm::isDebugEnabled()) {
0399     LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress] "
0400                                        << "\n\t adc " << adc << "\n\t adcPrev " << adcPrev << "\n\t adcNext " << adcNext
0401                                        << "\n\t adcMaxNeigh " << adcMaxNeigh << "\n\t adcPrev2 " << adcPrev2
0402                                        << "\n\t adcNext2 " << adcNext2 << std::endl;
0403 
0404     LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress] "
0405                                        << "\n\t theFEDlowThresh " << theFEDlowThresh << "\n\t theFEDhighThresh "
0406                                        << theFEDhighThresh << "\n\t thePrevFEDlowThresh " << thePrevFEDlowThresh
0407                                        << "\n\t thePrevFEDhighThresh " << thePrevFEDhighThresh
0408                                        << "\n\t theNextFEDlowThresh " << theNextFEDlowThresh
0409                                        << "\n\t theNextFEDhighThresh " << theNextFEDhighThresh
0410                                        << "\n\t theNeighFEDlowThresh " << theNeighFEDlowThresh
0411                                        << "\n\t theNeighFEDhighThresh " << theNeighFEDhighThresh
0412                                        << "\n\t thePrev2FEDlowThresh " << thePrev2FEDlowThresh
0413                                        << "\n\t theNext2FEDlowThresh " << theNext2FEDlowThresh << std::endl;
0414   }
0415 #endif
0416   // Decide if this strip should be accepted.
0417   bool accept = false;
0418   switch (theFEDalgorithm) {
0419     case 1:
0420       accept = (adc >= theFEDlowThresh);
0421       break;
0422     case 2:
0423       accept = (adc >= theFEDhighThresh || (adc >= theFEDlowThresh && adcMaxNeigh >= theNeighFEDlowThresh));
0424       break;
0425     case 3:
0426       accept = (adc >= theFEDhighThresh || (adc >= theFEDlowThresh && adcMaxNeigh >= theNeighFEDhighThresh));
0427       break;
0428     case 4:
0429       accept =
0430           ((adc >= theFEDhighThresh)     //Test for adc>highThresh (same as algorithm 2)
0431            || ((adc >= theFEDlowThresh)  //Test for adc>lowThresh, with neighbour adc>lowThresh (same as algorithm 2)
0432                && (adcMaxNeigh >= theNeighFEDlowThresh)) ||
0433            ((adc < theFEDlowThresh)                 //Test for adc<lowThresh
0434             && (((adcPrev >= thePrevFEDhighThresh)  //with both neighbours>highThresh
0435                  && (adcNext >= theNextFEDhighThresh)) ||
0436                 ((adcPrev >= thePrevFEDhighThresh)    //OR with previous neighbour>highThresh and
0437                  && (adcNext >= theNextFEDlowThresh)  //both the next neighbours>lowThresh
0438                  && (adcNext2 >= theNext2FEDlowThresh)) ||
0439                 ((adcNext >= theNextFEDhighThresh)    //OR with next neighbour>highThresh and
0440                  && (adcPrev >= thePrevFEDlowThresh)  //both the previous neighbours>lowThresh
0441                  && (adcPrev2 >= thePrev2FEDlowThresh)) ||
0442                 ((adcNext >= theNextFEDlowThresh)       //OR with both next neighbours>lowThresh and
0443                  && (adcNext2 >= theNext2FEDlowThresh)  //both the previous neighbours>lowThresh
0444                  && (adcPrev >= thePrevFEDlowThresh) && (adcPrev2 >= thePrev2FEDlowThresh)))));
0445       break;
0446     case 5:
0447       accept = adc > 0;
0448       break;
0449   }
0450   return accept;
0451 }
[ 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!