SiStripNoises.cc

Go to the documentation of this file.

#include "CondFormats/SiStripObjects/interface/SiStripNoises.h"
#include "FWCore/Utilities/interface/Exception.h"
#include <iostream>
#include <algorithm>
#include <cmath>
#include <iomanip>
#include "CondFormats/SiStripObjects/interface/SiStripDetSummary.h"

SiStripNoises::SiStripNoises(const SiStripNoises& input) {
  v_noises.clear();
  indexes.clear();
  v_noises.insert(v_noises.end(), input.v_noises.begin(), input.v_noises.end());
  indexes.insert(indexes.end(), input.indexes.begin(), input.indexes.end());
}

bool SiStripNoises::put(const uint32_t& DetId, const InputVector& input) {
  std::vector<unsigned char> Vo_CHAR;
  encode(input, Vo_CHAR);

  Registry::iterator p = std::lower_bound(indexes.begin(), indexes.end(), DetId, SiStripNoises::StrictWeakOrdering());
  if (p != indexes.end() && p->detid == DetId)
    return false;

  size_t sd = Vo_CHAR.end() - Vo_CHAR.begin();
  DetRegistry detregistry;
  detregistry.detid = DetId;
  detregistry.ibegin = v_noises.size();
  detregistry.iend = v_noises.size() + sd;
  indexes.insert(p, detregistry);
  v_noises.insert(v_noises.end(), Vo_CHAR.begin(), Vo_CHAR.end());
  return true;
}

const SiStripNoises::Range SiStripNoises::getRange(const uint32_t DetId) const {
  // get SiStripNoises Range of DetId

  RegistryIterator p = std::lower_bound(indexes.begin(), indexes.end(), DetId, SiStripNoises::StrictWeakOrdering());
  if (p == indexes.end() || p->detid != DetId)
    return SiStripNoises::Range(v_noises.end(), v_noises.end());
  else {
    __builtin_prefetch((&v_noises.front()) + p->ibegin);
    __builtin_prefetch((&v_noises.front()) + p->ibegin + 96);
    __builtin_prefetch((&v_noises.front()) + p->iend - 96);
    return SiStripNoises::Range(v_noises.begin() + p->ibegin, v_noises.begin() + p->iend);
  }
}

SiStripNoises::Range SiStripNoises::getRangeByPos(unsigned short pos) const {
  if (pos > indexes.size())
    return Range(v_noises.end(), v_noises.end());
  auto p = indexes.begin() + pos;
  __builtin_prefetch((&v_noises.front()) + p->ibegin);
  __builtin_prefetch((&v_noises.front()) + p->ibegin + 96);
  __builtin_prefetch((&v_noises.front()) + p->iend - 96);
  return Range(v_noises.begin() + p->ibegin, v_noises.begin() + p->iend);
}

void SiStripNoises::getDetIds(std::vector<uint32_t>& DetIds_) const {
  // returns vector of DetIds in map
  SiStripNoises::RegistryIterator begin = indexes.begin();
  SiStripNoises::RegistryIterator end = indexes.end();
  DetIds_.reserve(indexes.size());
  for (SiStripNoises::RegistryIterator p = begin; p != end; ++p) {
    DetIds_.push_back(p->detid);
  }
}

void SiStripNoises::verify(uint16_t strip, const Range& range) {
  if (9 * strip >= (range.second - range.first) * 8)
    throw cms::Exception("CorruptedData")
        << "[SiStripNoises::getNoise] looking for SiStripNoises for a strip out of range: strip " << strip;
}

void SiStripNoises::setData(float noise_, InputVector& v) {
  v.push_back((static_cast<int16_t>(noise_ * 10.0 + 0.5) & 0x01FF));
}

void SiStripNoises::encode(const InputVector& Vi, std::vector<unsigned char>& Vo) {
  static const uint16_t BITS_PER_STRIP = 9;
  const size_t VoSize = (size_t)((Vi.size() * BITS_PER_STRIP) / 8 + .999);
  Vo.resize(VoSize);
  for (size_t i = 0; i < VoSize; ++i)
    Vo[i] &= 0x00u;

  for (unsigned int stripIndex = 0; stripIndex < Vi.size(); ++stripIndex) {
    unsigned char* data = &Vo[VoSize - 1];
    uint32_t lowBit = stripIndex * BITS_PER_STRIP;
    uint8_t firstByteBit = (lowBit & 0x7);
    uint8_t firstByteNBits = 8 - firstByteBit;
    uint8_t firstByteMask = 0xffu << firstByteBit;
    uint8_t secondByteNbits = (BITS_PER_STRIP - firstByteNBits);
    uint8_t secondByteMask = ~(0xffu << secondByteNbits);

    *(data - lowBit / 8) = (*(data - lowBit / 8) & ~(firstByteMask)) | ((Vi[stripIndex] & 0xffu) << firstByteBit);
    *(data - lowBit / 8 - 1) =
        (*(data - lowBit / 8 - 1) & ~(secondByteMask)) | ((Vi[stripIndex] >> firstByteNBits) & secondByteMask);

    /*
    if(stripIndex   < 25 ){
      std::cout       << "***************ENCODE*********************"<<std::endl
              << "\tdata-lowBit/8     :"<<print_as_binary((*(data-lowBit/8)   & ~(firstByteMask)))
              << "-"<<print_as_binary(((Vi[stripIndex] & 0xffu) <<firstByteBit))
              << "\tdata-lowBit/8-1   :"<<print_as_binary((*(data-lowBit/8-1)   & ~(secondByteMask)))
              << "-"<<print_as_binary((((Vi[stripIndex]>> firstByteNBits) & secondByteMask)))
              << std::endl;
      std::cout       << "strip "<<stripIndex<<"\tvi: " << Vi[stripIndex] <<"\t"
              << print_short_as_binary(Vi[stripIndex])
              << "\tvo1:"<< print_char_as_binary(*(data-lowBit/8))
              << "\tvo2:"<< print_char_as_binary(*(data-lowBit/8-1))
              << "\tlowBit:"<< lowBit
              << "\tfirstByteMask :"<<print_as_binary(firstByteMask)
              << "\tsecondByteMask:"<<print_as_binary(secondByteMask)
              << "\tfirstByteBit:"<<print_as_binary(firstByteBit)
              << std::endl;
    }
    */
  }
}

//============ Methods for bulk-decoding all noises for a module ================

void SiStripNoises::allNoises(std::vector<float>& noises, const Range& range) const {
  size_t mysize = ((range.second - range.first) << 3) / 9;
  size_t size = noises.size();
  if (mysize < size)
    throw cms::Exception("CorruptedData") << "[SiStripNoises::allNoises] Requested noise for " << noises.size()
                                          << " strips, I have it only for " << mysize << " strips\n";
  size_t size8 = size & (~0x7), carry = size & 0x7;  // we have an optimized way of unpacking 8 strips
  const uint8_t* ptr = (&*range.second) - 1;
  std::vector<float>::iterator out = noises.begin(), end8 = noises.begin() + size8;
  // we do it this baroque way instead of just loopin on all the strips because it's faster
  // as the value of 'skip' is a constant, so the compiler can compute the masks directly
  while (out < end8) {
    *out = static_cast<float>(get9bits(ptr, 0) / 10.0f);
    ++out;
    *out = static_cast<float>(get9bits(ptr, 1) / 10.0f);
    ++out;
    *out = static_cast<float>(get9bits(ptr, 2) / 10.0f);
    ++out;
    *out = static_cast<float>(get9bits(ptr, 3) / 10.0f);
    ++out;
    *out = static_cast<float>(get9bits(ptr, 4) / 10.0f);
    ++out;
    *out = static_cast<float>(get9bits(ptr, 5) / 10.0f);
    ++out;
    *out = static_cast<float>(get9bits(ptr, 6) / 10.0f);
    ++out;
    *out = static_cast<float>(get9bits(ptr, 7) / 10.0f);
    ++out;
    --ptr;  // every 8 strips we have to skip one more bit
  }
  for (size_t rem = 0; rem < carry; ++rem) {
    *out = static_cast<float>(get9bits(ptr, rem) / 10.0f);
    ++out;
  }
}

/*
const std::string SiStripNoises::print_as_binary(const uint8_t ch) const
{
  std::string     str;
  int i = CHAR_BIT;
  while (i > 0)
    {
      -- i;
      str.push_back((ch&(1 << i) ? '1' : '0'));
    }
  return str;
}

std::string SiStripNoises::print_char_as_binary(const unsigned char ch) const
{
  std::string     str;
  int i = CHAR_BIT;
  while (i > 0)
    {
      -- i;
      str.push_back((ch&(1 << i) ? '1' : '0'));
    }
  return str;
}

std::string SiStripNoises::print_short_as_binary(const short ch) const
{
  std::string     str;
  int i = CHAR_BIT*2;
  while (i > 0)
    {
      -- i;
      str.push_back((ch&(1 << i) ? '1' : '0'));
    }
  return str;
}
*/

void SiStripNoises::printDebug(std::stringstream& ss, const TrackerTopology* /*trackerTopo*/) const {
  RegistryIterator rit = getRegistryVectorBegin(), erit = getRegistryVectorEnd();
  uint16_t Nstrips;
  std::vector<float> vstripnoise;

  ss << "detid" << std::setw(15) << "strip" << std::setw(10) << "noise" << std::endl;

  int detId = 0;
  int oldDetId = 0;
  for (; rit != erit; ++rit) {
    Nstrips = (rit->iend - rit->ibegin) * 8 / 9;  //number of strips = number of chars * char size / strip noise size
    vstripnoise.resize(Nstrips);
    allNoises(vstripnoise, make_pair(getDataVectorBegin() + rit->ibegin, getDataVectorBegin() + rit->iend));

    detId = rit->detid;
    if (detId != oldDetId) {
      oldDetId = detId;
      ss << detId;
    } else
      ss << "   ";
    for (size_t i = 0; i < Nstrips; ++i) {
      if (i != 0)
        ss << "   ";
      ss << std::setw(15) << i << std::setw(10) << vstripnoise[i] << std::endl;
    }
  }
}

void SiStripNoises::printSummary(std::stringstream& ss, const TrackerTopology* trackerTopo) const {
  SiStripDetSummary summary{trackerTopo};

  std::stringstream tempss;

  RegistryIterator rit = getRegistryVectorBegin(), erit = getRegistryVectorEnd();
  uint16_t Nstrips;
  std::vector<float> vstripnoise;
  double mean, rms, min, max;
  for (; rit != erit; ++rit) {
    Nstrips = (rit->iend - rit->ibegin) * 8 / 9;  //number of strips = number of chars * char size / strip noise size
    vstripnoise.resize(Nstrips);
    allNoises(vstripnoise, make_pair(getDataVectorBegin() + rit->ibegin, getDataVectorBegin() + rit->iend));
    tempss << "\ndetid: " << rit->detid << " \t ";
    mean = 0;
    rms = 0;
    min = 10000;
    max = 0;

    DetId detId(rit->detid);

    for (size_t i = 0; i < Nstrips; ++i) {
      mean += vstripnoise[i];
      rms += vstripnoise[i] * vstripnoise[i];
      if (vstripnoise[i] < min)
        min = vstripnoise[i];
      if (vstripnoise[i] > max)
        max = vstripnoise[i];

      summary.add(detId, vstripnoise[i]);
    }
    mean /= Nstrips;
    rms = sqrt(rms / Nstrips - mean * mean);

    tempss << "Nstrips " << Nstrips << " \t; mean " << mean << " \t; rms " << rms << " \t; min " << min << " \t; max "
           << max << "\t ";
  }
  ss << std::endl << "Summary:" << std::endl;
  summary.print(ss);
  ss << std::endl;
  ss << tempss.str();
}

std::vector<SiStripNoises::ratioData> SiStripNoises::operator/(const SiStripNoises& d) {
  std::vector<ratioData> result;
  ratioData aData;

  RegistryIterator iter = getRegistryVectorBegin();
  RegistryIterator iterE = getRegistryVectorEnd();

  //Divide result by d
  for (; iter != iterE; ++iter) {
    float value;
    //get noise from d
    aData.detid = iter->detid;
    aData.values.clear();
    Range d_range = d.getRange(iter->detid);
    Range range = Range(v_noises.begin() + iter->ibegin, v_noises.begin() + iter->iend);

    //if denominator is missing, put the ratio value to 0xFFFF (=inf)
    size_t strip = 0, stripE = (range.second - range.first) * 8 / 9;
    for (; strip < stripE; ++strip) {
      if (d_range.first == d_range.second) {
        value = 0xFFFF;
      } else {
        value = getNoise(strip, range) / d.getNoise(strip, d_range);
      }
      aData.values.push_back(value);
    }
    result.push_back(aData);
  }

  iter = d.getRegistryVectorBegin();
  iterE = d.getRegistryVectorEnd();

  //Divide result by d
  for (; iter != iterE; ++iter) {
    float value;
    //get noise from d
    Range range = this->getRange(iter->detid);
    Range d_range = Range(d.v_noises.begin() + iter->ibegin, d.v_noises.begin() + iter->iend);
    if (range.first == range.second) {
      aData.detid = iter->detid;
      aData.values.clear();
      size_t strip = 0, stripE = (d_range.second - d_range.first) * 8 / 9;
      for (; strip < stripE; ++strip) {
        value = 0.;
        aData.values.push_back(value);
      }
      result.push_back(aData);
    }
  }

  return result;
}