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