SiStripObjects/src/SiStripThreshold.cc

0001 #include "CondFormats/SiStripObjects/interface/SiStripThreshold.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 #include <cassert>
0005 #include <algorithm>
0006 #include <cmath>
0007
0008 bool SiStripThreshold::put(const uint32_t& DetId, const InputVector& _vect) {
0009   InputVector vect = _vect;
0010   // put in SiStripThreshold::v_threshold of DetId
0011   Registry::iterator p =
0012       std::lower_bound(indexes.begin(), indexes.end(), DetId, SiStripThreshold::StrictWeakOrdering());
0013   if (p != indexes.end() && p->detid == DetId) {
0014     edm::LogError("SiStripThreshold") << "[" << __PRETTY_FUNCTION__ << "] SiStripThreshold for DetID " << DetId
0015                                       << " is already stored. Skipping this put" << std::endl;
0016     return false;
0017   }
0018
0019   SiStripThreshold::Container::iterator new_end = compact(vect);
0020
0021   size_t sd = new_end - vect.begin();
0022   DetRegistry detregistry;
0023   detregistry.detid = DetId;
0024   detregistry.ibegin = v_threshold.size();
0025   detregistry.iend = v_threshold.size() + sd;
0026   indexes.insert(p, detregistry);
0027
0028   v_threshold.insert(v_threshold.end(), vect.begin(), new_end);
0029
0030   return true;
0031 }
0032
0033 SiStripThreshold::Container::iterator SiStripThreshold::compact(Container& input) {
0034   std::stable_sort(input.begin(), input.end());
0035   return std::unique(input.begin(), input.end());
0036 }
0037
0038 const SiStripThreshold::Range SiStripThreshold::getRange(const uint32_t& DetId) const {
0039   // get SiStripThreshold Range of DetId
0040
0041   RegistryIterator p = std::lower_bound(indexes.begin(), indexes.end(), DetId, SiStripThreshold::StrictWeakOrdering());
0042   if (p == indexes.end() || p->detid != DetId)
0043     return SiStripThreshold::Range(v_threshold.end(), v_threshold.end());
0044   else
0045     return SiStripThreshold::Range(v_threshold.begin() + p->ibegin, v_threshold.begin() + p->iend);
0046 }
0047
0048 void SiStripThreshold::getDetIds(std::vector<uint32_t>& DetIds_) const {
0049   // returns vector of DetIds in map
0050   SiStripThreshold::RegistryIterator begin = indexes.begin();
0051   SiStripThreshold::RegistryIterator end = indexes.end();
0052   for (SiStripThreshold::RegistryIterator p = begin; p != end; ++p) {
0053     DetIds_.push_back(p->detid);
0054   }
0055 }
0056
0057 void SiStripThreshold::setData(const uint16_t& strip, const float& lTh, const float& hTh, Container& vthr) {
0058   Data a;
0059   a.encode(strip, lTh, hTh);
0060   vthr.push_back(a);
0061 }
0062
0063 void SiStripThreshold::setData(
0064     const uint16_t& strip, const float& lTh, const float& hTh, const float& cTh, Container& vthr) {
0065   Data a;
0066   a.encode(strip, lTh, hTh, cTh);
0067   vthr.push_back(a);
0068 }
0069
0070 SiStripThreshold::Data SiStripThreshold::getData(const uint16_t& strip, const Range& range) const {
0071   uint16_t estrip =
0072       (strip & sistrip::FirstThStripMask_) << sistrip::FirstThStripShift_ | (63 & sistrip::HighThStripMask_);
0073   ContainerIterator p = std::upper_bound(range.first, range.second, estrip, SiStripThreshold::dataStrictWeakOrdering());
0074   if (p != range.first) {
0075     return *(--p);
0076   } else {
0077     throw cms::Exception("CorruptedData") << "[SiStripThreshold::getData] asking for data for a strip " << strip
0078                                           << " lower then the first stored strip " << p->getFirstStrip();
0079   }
0080 }
0081
0082 void SiStripThreshold::allThresholds(std::vector<float>& lowThs,
0083                                      std::vector<float>& highThs,
0084                                      const Range& range) const {
0085   ContainerIterator it = range.first;
0086   size_t strips = lowThs.size();
0087   assert(strips == highThs.size());
0088   while (it != range.second) {
0089     size_t firstStrip = it->getFirstStrip();
0090     //std::cout << "First strip is " << firstStrip << std::endl;
0091     float high = it->getHth(), low = it->getLth();
0092     //std::cout << "High is " << high << ", low is " << low << std::endl;
0093     ++it;  // increment the pointer
0094     size_t lastStrip = (it == range.second ? strips : it->getFirstStrip());
0095     //std::cout << "Last strip is " << lastStrip << std::endl;
0096     if (lastStrip > strips) {
0097       it = range.second;   // I should stop here,
0098       lastStrip = strips;  // and fill only 'strips' strips
0099     }
0100     std::fill(&lowThs[firstStrip], &lowThs[lastStrip], low);
0101     std::fill(&highThs[firstStrip], &highThs[lastStrip], high);
0102   }
0103 }
0104
0105 void SiStripThreshold::printDebug(std::stringstream& ss, const TrackerTopology* /*trackerTopo*/) const {
0106   RegistryIterator rit = getRegistryVectorBegin(), erit = getRegistryVectorEnd();
0107   ContainerIterator it, eit;
0108   for (; rit != erit; ++rit) {
0109     it = getDataVectorBegin() + rit->ibegin;
0110     eit = getDataVectorBegin() + rit->iend;
0111     ss << "\ndetid: " << rit->detid << " \t ";
0112     for (; it != eit; ++it) {
0113       ss << "\n \t ";
0114       it->print(ss);
0115     }
0116   }
0117 }
0118
0119 void SiStripThreshold::printSummary(std::stringstream& ss, const TrackerTopology* /*trackerTopo*/) const {
0120   RegistryIterator rit = getRegistryVectorBegin(), erit = getRegistryVectorEnd();
0121   ContainerIterator it, eit, itp;
0122   float meanLth, meanHth, meanCth;  //mean value
0123   float rmsLth, rmsHth, rmsCth;     //rms value
0124   float maxLth, maxHth, maxCth;     //max value
0125   float minLth, minHth, minCth;     //min value
0126   uint16_t n;
0127   uint16_t firstStrip, stripRange;
0128   for (; rit != erit; ++rit) {
0129     it = getDataVectorBegin() + rit->ibegin;
0130     eit = getDataVectorBegin() + rit->iend;
0131     ss << "\ndetid: " << rit->detid << " \t ";
0132
0133     meanLth = 0;
0134     meanHth = 0;
0135     meanCth = 0;  //mean value
0136     rmsLth = 0;
0137     rmsHth = 0;
0138     rmsCth = 0;  //rms value
0139     maxLth = 0;
0140     maxHth = 0;
0141     maxCth = 0;  //max value
0142     minLth = 10000;
0143     minHth = 10000;
0144     minCth = 10000;  //min value
0145     n = 0;
0146     firstStrip = 0;
0147     for (; it != eit; ++it) {
0148       itp = it + 1;
0149       firstStrip = it->getFirstStrip();
0150       if (itp != eit)
0151         stripRange = (itp->getFirstStrip() - firstStrip);
0152       else
0153         stripRange =
0154             firstStrip > 511
0155                 ? 768 - firstStrip
0156                 : 512 -
0157                       firstStrip;  //*FIXME, I dont' know ithis class the strip number of a detector, so I assume wrongly that if the last firstStrip<511 the detector has only 512 strips. Clearly wrong. to be fixed
0158
0159       addToStat(it->getLth(), stripRange, meanLth, rmsLth, minLth, maxLth);
0160       addToStat(it->getHth(), stripRange, meanHth, rmsHth, minHth, maxHth);
0161       addToStat(it->getClusth(), stripRange, meanCth, rmsCth, minCth, maxCth);
0162       n += stripRange;
0163     }
0164     meanLth /= n;
0165     meanHth /= n;
0166     meanCth /= n;
0167     rmsLth = sqrt(rmsLth / n - meanLth * meanLth);
0168     rmsHth = sqrt(rmsHth / n - meanHth * meanHth);
0169     rmsCth = sqrt(rmsCth / n - meanCth * meanCth);
0170     ss << "\nn " << n << " \tmeanLth " << meanLth << " \t rmsLth " << rmsLth << " \t minLth " << minLth << " \t maxLth "
0171        << maxLth;
0172     ss << "\n\tmeanHth " << meanHth << " \t rmsHth " << rmsHth << " \t minHth " << minHth << " \t maxHth " << maxHth;
0173     ss << "\n\tmeanCth " << meanCth << " \t rmsCth " << rmsCth << " \t minCth " << minCth << " \t maxCth " << maxCth;
0174   }
0175 }
0176
0177 void SiStripThreshold::addToStat(float value, uint16_t& range, float& sum, float& sum2, float& min, float& max) const {
0178   sum += value * range;
0179   sum2 += value * value * range;
0180   if (value < min)
0181     min = value;
0182   if (value > max)
0183     max = value;
0184 }