Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-05-12 02:41:30 

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