SiStripObjects/src/SiStripNoises.cc

0001 #include "CondFormats/SiStripObjects/interface/SiStripNoises.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003 #include <iostream>
0004 #include <algorithm>
0005 #include <cmath>
0006 #include <iomanip>
0007 #include "CondFormats/SiStripObjects/interface/SiStripDetSummary.h"
0008
0009 SiStripNoises::SiStripNoises(const SiStripNoises& input) {
0010   v_noises.clear();
0011   indexes.clear();
0012   v_noises.insert(v_noises.end(), input.v_noises.begin(), input.v_noises.end());
0013   indexes.insert(indexes.end(), input.indexes.begin(), input.indexes.end());
0014 }
0015
0016 bool SiStripNoises::put(const uint32_t& DetId, const InputVector& input) {
0017   std::vector<unsigned char> Vo_CHAR;
0018   encode(input, Vo_CHAR);
0019
0020   Registry::iterator p = std::lower_bound(indexes.begin(), indexes.end(), DetId, SiStripNoises::StrictWeakOrdering());
0021   if (p != indexes.end() && p->detid == DetId)
0022     return false;
0023
0024   size_t sd = Vo_CHAR.end() - Vo_CHAR.begin();
0025   DetRegistry detregistry;
0026   detregistry.detid = DetId;
0027   detregistry.ibegin = v_noises.size();
0028   detregistry.iend = v_noises.size() + sd;
0029   indexes.insert(p, detregistry);
0030   v_noises.insert(v_noises.end(), Vo_CHAR.begin(), Vo_CHAR.end());
0031   return true;
0032 }
0033
0034 const SiStripNoises::Range SiStripNoises::getRange(const uint32_t DetId) const {
0035   // get SiStripNoises Range of DetId
0036
0037   RegistryIterator p = std::lower_bound(indexes.begin(), indexes.end(), DetId, SiStripNoises::StrictWeakOrdering());
0038   if (p == indexes.end() || p->detid != DetId)
0039     return SiStripNoises::Range(v_noises.end(), v_noises.end());
0040   else {
0041     __builtin_prefetch((&v_noises.front()) + p->ibegin);
0042     __builtin_prefetch((&v_noises.front()) + p->ibegin + 96);
0043     __builtin_prefetch((&v_noises.front()) + p->iend - 96);
0044     return SiStripNoises::Range(v_noises.begin() + p->ibegin, v_noises.begin() + p->iend);
0045   }
0046 }
0047
0048 SiStripNoises::Range SiStripNoises::getRangeByPos(unsigned short pos) const {
0049   if (pos > indexes.size())
0050     return Range(v_noises.end(), v_noises.end());
0051   auto p = indexes.begin() + pos;
0052   __builtin_prefetch((&v_noises.front()) + p->ibegin);
0053   __builtin_prefetch((&v_noises.front()) + p->ibegin + 96);
0054   __builtin_prefetch((&v_noises.front()) + p->iend - 96);
0055   return Range(v_noises.begin() + p->ibegin, v_noises.begin() + p->iend);
0056 }
0057
0058 void SiStripNoises::getDetIds(std::vector<uint32_t>& DetIds_) const {
0059   // returns vector of DetIds in map
0060   SiStripNoises::RegistryIterator begin = indexes.begin();
0061   SiStripNoises::RegistryIterator end = indexes.end();
0062   DetIds_.reserve(indexes.size());
0063   for (SiStripNoises::RegistryIterator p = begin; p != end; ++p) {
0064     DetIds_.push_back(p->detid);
0065   }
0066 }
0067
0068 void SiStripNoises::verify(uint16_t strip, const Range& range) {
0069   if (9 * strip >= (range.second - range.first) * 8)
0070     throw cms::Exception("CorruptedData")
0071         << "[SiStripNoises::getNoise] looking for SiStripNoises for a strip out of range: strip " << strip;
0072 }
0073
0074 void SiStripNoises::setData(float noise_, InputVector& v) {
0075   v.push_back((static_cast<int16_t>(noise_ * 10.0 + 0.5) & 0x01FF));
0076 }
0077
0078 void SiStripNoises::encode(const InputVector& Vi, std::vector<unsigned char>& Vo) {
0079   static const uint16_t BITS_PER_STRIP = 9;
0080   const size_t VoSize = (size_t)((Vi.size() * BITS_PER_STRIP) / 8 + .999);
0081   Vo.resize(VoSize);
0082   for (size_t i = 0; i < VoSize; ++i)
0083     Vo[i] &= 0x00u;
0084
0085   for (unsigned int stripIndex = 0; stripIndex < Vi.size(); ++stripIndex) {
0086     unsigned char* data = &Vo[VoSize - 1];
0087     uint32_t lowBit = stripIndex * BITS_PER_STRIP;
0088     uint8_t firstByteBit = (lowBit & 0x7);
0089     uint8_t firstByteNBits = 8 - firstByteBit;
0090     uint8_t firstByteMask = 0xffu << firstByteBit;
0091     uint8_t secondByteNbits = (BITS_PER_STRIP - firstByteNBits);
0092     uint8_t secondByteMask = ~(0xffu << secondByteNbits);
0093
0094     *(data - lowBit / 8) = (*(data - lowBit / 8) & ~(firstByteMask)) | ((Vi[stripIndex] & 0xffu) << firstByteBit);
0095     *(data - lowBit / 8 - 1) =
0096         (*(data - lowBit / 8 - 1) & ~(secondByteMask)) | ((Vi[stripIndex] >> firstByteNBits) & secondByteMask);
0097
0098     /*
0099     if(stripIndex   < 25 ){
0100       std::cout       << "***************ENCODE*********************"<<std::endl
0101               << "\tdata-lowBit/8     :"<<print_as_binary((*(data-lowBit/8)   & ~(firstByteMask)))
0102               << "-"<<print_as_binary(((Vi[stripIndex] & 0xffu) <<firstByteBit))
0103               << "\tdata-lowBit/8-1   :"<<print_as_binary((*(data-lowBit/8-1)   & ~(secondByteMask)))
0104               << "-"<<print_as_binary((((Vi[stripIndex]>> firstByteNBits) & secondByteMask)))
0105               << std::endl;
0106       std::cout       << "strip "<<stripIndex<<"\tvi: " << Vi[stripIndex] <<"\t"
0107               << print_short_as_binary(Vi[stripIndex])
0108               << "\tvo1:"<< print_char_as_binary(*(data-lowBit/8))
0109               << "\tvo2:"<< print_char_as_binary(*(data-lowBit/8-1))
0110               << "\tlowBit:"<< lowBit
0111               << "\tfirstByteMask :"<<print_as_binary(firstByteMask)
0112               << "\tsecondByteMask:"<<print_as_binary(secondByteMask)
0113               << "\tfirstByteBit:"<<print_as_binary(firstByteBit)
0114               << std::endl;
0115     }
0116     */
0117   }
0118 }
0119
0120 //============ Methods for bulk-decoding all noises for a module ================
0121
0122 void SiStripNoises::allNoises(std::vector<float>& noises, const Range& range) const {
0123   size_t mysize = ((range.second - range.first) << 3) / 9;
0124   size_t size = noises.size();
0125   if (mysize < size)
0126     throw cms::Exception("CorruptedData") << "[SiStripNoises::allNoises] Requested noise for " << noises.size()
0127                                           << " strips, I have it only for " << mysize << " strips\n";
0128   size_t size8 = size & (~0x7), carry = size & 0x7;  // we have an optimized way of unpacking 8 strips
0129   const uint8_t* ptr = (&*range.second) - 1;
0130   std::vector<float>::iterator out = noises.begin(), end8 = noises.begin() + size8;
0131   // we do it this baroque way instead of just loopin on all the strips because it's faster
0132   // as the value of 'skip' is a constant, so the compiler can compute the masks directly
0133   while (out < end8) {
0134     *out = static_cast<float>(get9bits(ptr, 0) / 10.0f);
0135     ++out;
0136     *out = static_cast<float>(get9bits(ptr, 1) / 10.0f);
0137     ++out;
0138     *out = static_cast<float>(get9bits(ptr, 2) / 10.0f);
0139     ++out;
0140     *out = static_cast<float>(get9bits(ptr, 3) / 10.0f);
0141     ++out;
0142     *out = static_cast<float>(get9bits(ptr, 4) / 10.0f);
0143     ++out;
0144     *out = static_cast<float>(get9bits(ptr, 5) / 10.0f);
0145     ++out;
0146     *out = static_cast<float>(get9bits(ptr, 6) / 10.0f);
0147     ++out;
0148     *out = static_cast<float>(get9bits(ptr, 7) / 10.0f);
0149     ++out;
0150     --ptr;  // every 8 strips we have to skip one more bit
0151   }
0152   for (size_t rem = 0; rem < carry; ++rem) {
0153     *out = static_cast<float>(get9bits(ptr, rem) / 10.0f);
0154     ++out;
0155   }
0156 }
0157
0158 /*
0159 const std::string SiStripNoises::print_as_binary(const uint8_t ch) const
0160 {
0161   std::string     str;
0162   int i = CHAR_BIT;
0163   while (i > 0)
0164     {
0165       -- i;
0166       str.push_back((ch&(1 << i) ? '1' : '0'));
0167     }
0168   return str;
0169 }
0170
0171 std::string SiStripNoises::print_char_as_binary(const unsigned char ch) const
0172 {
0173   std::string     str;
0174   int i = CHAR_BIT;
0175   while (i > 0)
0176     {
0177       -- i;
0178       str.push_back((ch&(1 << i) ? '1' : '0'));
0179     }
0180   return str;
0181 }
0182
0183 std::string SiStripNoises::print_short_as_binary(const short ch) const
0184 {
0185   std::string     str;
0186   int i = CHAR_BIT*2;
0187   while (i > 0)
0188     {
0189       -- i;
0190       str.push_back((ch&(1 << i) ? '1' : '0'));
0191     }
0192   return str;
0193 }
0194 */
0195
0196 void SiStripNoises::printDebug(std::stringstream& ss, const TrackerTopology* /*trackerTopo*/) const {
0197   RegistryIterator rit = getRegistryVectorBegin(), erit = getRegistryVectorEnd();
0198   uint16_t Nstrips;
0199   std::vector<float> vstripnoise;
0200
0201   ss << "detid" << std::setw(15) << "strip" << std::setw(10) << "noise" << std::endl;
0202
0203   int detId = 0;
0204   int oldDetId = 0;
0205   for (; rit != erit; ++rit) {
0206     Nstrips = (rit->iend - rit->ibegin) * 8 / 9;  //number of strips = number of chars * char size / strip noise size
0207     vstripnoise.resize(Nstrips);
0208     allNoises(vstripnoise, make_pair(getDataVectorBegin() + rit->ibegin, getDataVectorBegin() + rit->iend));
0209
0210     detId = rit->detid;
0211     if (detId != oldDetId) {
0212       oldDetId = detId;
0213       ss << detId;
0214     } else
0215       ss << "   ";
0216     for (size_t i = 0; i < Nstrips; ++i) {
0217       if (i != 0)
0218         ss << "   ";
0219       ss << std::setw(15) << i << std::setw(10) << vstripnoise[i] << std::endl;
0220     }
0221   }
0222 }
0223
0224 void SiStripNoises::printSummary(std::stringstream& ss, const TrackerTopology* trackerTopo) const {
0225   SiStripDetSummary summary{trackerTopo};
0226
0227   std::stringstream tempss;
0228
0229   RegistryIterator rit = getRegistryVectorBegin(), erit = getRegistryVectorEnd();
0230   uint16_t Nstrips;
0231   std::vector<float> vstripnoise;
0232   double mean, rms, min, max;
0233   for (; rit != erit; ++rit) {
0234     Nstrips = (rit->iend - rit->ibegin) * 8 / 9;  //number of strips = number of chars * char size / strip noise size
0235     vstripnoise.resize(Nstrips);
0236     allNoises(vstripnoise, make_pair(getDataVectorBegin() + rit->ibegin, getDataVectorBegin() + rit->iend));
0237     tempss << "\ndetid: " << rit->detid << " \t ";
0238     mean = 0;
0239     rms = 0;
0240     min = 10000;
0241     max = 0;
0242
0243     DetId detId(rit->detid);
0244
0245     for (size_t i = 0; i < Nstrips; ++i) {
0246       mean += vstripnoise[i];
0247       rms += vstripnoise[i] * vstripnoise[i];
0248       if (vstripnoise[i] < min)
0249         min = vstripnoise[i];
0250       if (vstripnoise[i] > max)
0251         max = vstripnoise[i];
0252
0253       summary.add(detId, vstripnoise[i]);
0254     }
0255     mean /= Nstrips;
0256     rms = sqrt(rms / Nstrips - mean * mean);
0257
0258     tempss << "Nstrips " << Nstrips << " \t; mean " << mean << " \t; rms " << rms << " \t; min " << min << " \t; max "
0259            << max << "\t ";
0260   }
0261   ss << std::endl << "Summary:" << std::endl;
0262   summary.print(ss);
0263   ss << std::endl;
0264   ss << tempss.str();
0265 }
0266
0267 std::vector<SiStripNoises::ratioData> SiStripNoises::operator/(const SiStripNoises& d) {
0268   std::vector<ratioData> result;
0269   ratioData aData;
0270
0271   RegistryIterator iter = getRegistryVectorBegin();
0272   RegistryIterator iterE = getRegistryVectorEnd();
0273
0274   //Divide result by d
0275   for (; iter != iterE; ++iter) {
0276     float value;
0277     //get noise from d
0278     aData.detid = iter->detid;
0279     aData.values.clear();
0280     Range d_range = d.getRange(iter->detid);
0281     Range range = Range(v_noises.begin() + iter->ibegin, v_noises.begin() + iter->iend);
0282
0283     //if denominator is missing, put the ratio value to 0xFFFF (=inf)
0284     size_t strip = 0, stripE = (range.second - range.first) * 8 / 9;
0285     for (; strip < stripE; ++strip) {
0286       if (d_range.first == d_range.second) {
0287         value = 0xFFFF;
0288       } else {
0289         value = getNoise(strip, range) / d.getNoise(strip, d_range);
0290       }
0291       aData.values.push_back(value);
0292     }
0293     result.push_back(aData);
0294   }
0295
0296   iter = d.getRegistryVectorBegin();
0297   iterE = d.getRegistryVectorEnd();
0298
0299   //Divide result by d
0300   for (; iter != iterE; ++iter) {
0301     float value;
0302     //get noise from d
0303     Range range = this->getRange(iter->detid);
0304     Range d_range = Range(d.v_noises.begin() + iter->ibegin, d.v_noises.begin() + iter->iend);
0305     if (range.first == range.second) {
0306       aData.detid = iter->detid;
0307       aData.values.clear();
0308       size_t strip = 0, stripE = (d_range.second - d_range.first) * 8 / 9;
0309       for (; strip < stripE; ++strip) {
0310         value = 0.;
0311         aData.values.push_back(value);
0312       }
0313       result.push_back(aData);
0314     }
0315   }
0316
0317   return result;
0318 }