SiStripClusterizer/src/SiStripApvShotCleaner.cc

0001 #include "RecoLocalTracker/SiStripClusterizer/interface/SiStripApvShotCleaner.h"
0002 #include <algorithm>
0003 #include <memory>
0004
0005 //Uncomment the following #define to have print debug
0006 //#define DEBUGME
0007
0008 SiStripApvShotCleaner::SiStripApvShotCleaner()
0009     : maxNumOfApvs(6),  //FED Default: 6 (i.e. max num apvs )
0010       stripsPerApv(128),
0011       stripsForMedian(64) {}
0012
0013 bool SiStripApvShotCleaner::clean(const edm::DetSet<SiStripDigi>& in,
0014                                   edm::DetSet<SiStripDigi>::const_iterator& scan,
0015                                   edm::DetSet<SiStripDigi>::const_iterator& end) {
0016   if (in.size() < 64)
0017     return false;
0018
0019   if (loop(in)) {
0020     reset(scan, end);
0021     return true;
0022   }
0023   return false;
0024 }
0025
0026 bool SiStripApvShotCleaner::loop(const edm::DetSet<SiStripDigi>& in) {
0027 #ifdef DEBUGME
0028   std::stringstream ss;
0029   ss << __func__ << " working on detid " << in.detId() << " for a digi.size=" << in.size();
0030 #endif
0031
0032   shots_ = false;
0033   for (auto& val : shotApv_)
0034     val = false;
0035
0036   cacheDetId = in.detId();
0037
0038   //Find the position in the DetSet where the first strip of an apv should be inserted
0039   // needed to deduce if at least stripsForMedian strips per apv have been fired
0040   for (size_t i = 0; i <= maxNumOfApvs; ++i) {
0041     SiStripDigi d(i * stripsPerApv, 0);  //Fake digi, at the edge of the apv
0042     pFirstDigiOfApv[i] = std::lower_bound(in.begin(), in.end(), d);
0043
0044     //if satisfied it means that the number of digis in the apv i-1 is above stripsForMedia -> apvShot
0045     if (i > 0 && pFirstDigiOfApv[i] - pFirstDigiOfApv[i - 1] > stripsForMedian) {
0046       shots_ = true;
0047       shotApv_[i - 1] = true;
0048 #ifdef DEBUGME
0049       ss << " found an apv shot of " << pFirstDigiOfApv[i] - pFirstDigiOfApv[i - 1] << " digis in detid " << in.detId()
0050          << " apv " << i << std::endl;
0051 #endif
0052     }
0053
0054     //---------------------
0055     //Just for debug REMOVE
0056     /*
0057     if(i>0){
0058       ss << "detid " << in.detId() << " apv " << i-1 << " number digis " << pFirstDigiOfApv[i]-pFirstDigiOfApv[i-1] << " \t shot " << shotApv_[i-1] << std::endl;
0059       if(pFirstDigiOfApv[i]-pFirstDigiOfApv[i-1]>stripsForMedian-2){
0060     edm::DetSet<SiStripDigi>::const_iterator dig=pFirstDigiOfApv[i-1];
0061     while(dig!=pFirstDigiOfApv[i]){
0062       ss << "\t strip " << dig->strip() << " dig.adc " << dig->adc();
0063       dig++;
0064     }
0065     ss << std::endl;
0066       }
0067     }
0068     */
0069     //-------------------------------
0070   }
0071
0072 #ifdef DEBUGME
0073   edm::LogInfo("ApvShot") << ss.str();
0074 #endif
0075
0076   if (!shots_)
0077     return false;
0078
0079   dumpInVector(pFirstDigiOfApv, maxNumOfApvs);
0080   return true;
0081 }
0082
0083 void SiStripApvShotCleaner::dumpInVector(edm::DetSet<SiStripDigi>::const_iterator* pFirstDigiOfApv,
0084                                          size_t maxNumOfApvs) {
0085   vdigis.clear();
0086   //loop on Apvs and remove shots. if an apv doesn't have shots, copy it
0087   for (size_t i = 0; i < maxNumOfApvs; ++i) {
0088     apvDigis.clear();
0089
0090     if (shotApv_[i]) {
0091       apvDigis.insert(apvDigis.end(), pFirstDigiOfApv[i], pFirstDigiOfApv[i + 1]);
0092       subtractCM();
0093       std::stable_sort(apvDigis.begin(), apvDigis.end());
0094       vdigis.insert(vdigis.end(), apvDigis.begin(), apvDigis.end());
0095     } else {
0096       vdigis.insert(vdigis.end(), pFirstDigiOfApv[i], pFirstDigiOfApv[i + 1]);
0097     }
0098   }
0099
0100 #ifdef DEBUGME
0101   std::stringstream ss;
0102   ss << "detid " << cacheDetId << " new digi.size " << vdigis.size() << "\n";
0103   for (size_t i = 0; i < vdigis.size(); ++i)
0104     ss << "\t " << i << " strip " << vdigis[i].strip() << " adc " << vdigis[i].adc();
0105   edm::LogInfo("ApvShot") << ss.str() << std::endl;
0106 #endif
0107 }
0108
0109 void SiStripApvShotCleaner::subtractCM() {
0110   //order by charge
0111   std::stable_sort(
0112       apvDigis.begin(), apvDigis.end(), [](SiStripDigi const& a, SiStripDigi const& b) { return a.adc() > b.adc(); });
0113
0114   //ignore case where 64th strip is 0ADC
0115   if (apvDigis[stripsForMedian].adc() == 0) {
0116 #ifdef DEBUGME
0117     std::stringstream ss;
0118     ss << "case with strip64=0 --> detid= " << cacheDetId << "\n";
0119     edm::LogInfo("ApvShot") << ss.str();
0120 #endif
0121     return;
0122   }
0123
0124   //Find the Median
0125   float CM = 0.5f * (apvDigis[stripsForMedian].adc() + apvDigis[stripsForMedian - 1].adc());
0126
0127   if (CM <= 0)
0128     return;
0129
0130   //Subtract the median
0131   const bool is10bit = apvDigis[0].adc() > 255;  // approximation; definitely 10bit in this case
0132   size_t i = 0;
0133   for (; i < stripsForMedian && apvDigis[i].adc() > CM; ++i) {
0134     const uint16_t adc =
0135         ((apvDigis[i].adc() > 253) && !is10bit) ? apvDigis[i].adc() : (uint16_t)(apvDigis[i].adc() - CM);
0136     apvDigis[i] = SiStripDigi(apvDigis[i].strip(), adc);
0137   }
0138   apvDigis.resize(i);
0139
0140 #ifdef DEBUGME
0141   std::stringstream ss;
0142   ss << "[subtractCM] detid " << cacheDetId << "  CM is " << CM << " the remaining strips after CM subtraction are  "
0143      << i;
0144   edm::LogInfo("ApvShot") << ss.str();
0145 #endif
0146 }
0147
0148 void SiStripApvShotCleaner::reset(edm::DetSet<SiStripDigi>::const_iterator& a,
0149                                   edm::DetSet<SiStripDigi>::const_iterator& b) {
0150   pDetSet = std::make_unique<edm::DetSet<SiStripDigi>>(cacheDetId);
0151   pDetSet->data.swap(vdigis);
0152   a = pDetSet->begin();
0153   b = pDetSet->end();
0154 }