/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/CondFormats/SiStripObjects/src/SiStripThreshold.cc

Go to the documentation of this file.

#include "CondFormats/SiStripObjects/interface/SiStripThreshold.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"
#include <cassert>
#include <algorithm>
#include <math.h>

bool SiStripThreshold::put(const uint32_t& DetId, InputVector vect) {
  // put in SiStripThreshold::v_threshold of DetId
  Registry::iterator p = std::lower_bound(indexes.begin(),indexes.end(),DetId,SiStripThreshold::StrictWeakOrdering());
  if (p!=indexes.end() && p->detid==DetId){
    edm::LogError("SiStripThreshold") << "[" << __PRETTY_FUNCTION__ << "] SiStripThreshold for DetID " << DetId << " is already stored. Skippig this put" << std::endl;
    return false;
  }
  
  SiStripThreshold::Container::iterator new_end=compact(vect);

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

SiStripThreshold::Container::iterator SiStripThreshold::compact(Container& input) {
  std::stable_sort(input.begin(),input.end());
  return std::unique(input.begin(),input.end());
}

const SiStripThreshold::Range SiStripThreshold::getRange(const uint32_t& DetId) const {
  // get SiStripThreshold Range of DetId
  
  RegistryIterator p = std::lower_bound(indexes.begin(),indexes.end(),DetId,SiStripThreshold::StrictWeakOrdering());
  if (p==indexes.end()|| p->detid!=DetId) 
    return SiStripThreshold::Range(v_threshold.end(),v_threshold.end()); 
  else 
    return SiStripThreshold::Range(v_threshold.begin()+p->ibegin,v_threshold.begin()+p->iend);
}


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

void SiStripThreshold::setData(const uint16_t& strip, const float& lTh,const float& hTh, Container& vthr){
  Data a;
  a.encode(strip,lTh,hTh);
  vthr.push_back(a);
}

void SiStripThreshold::setData(const uint16_t& strip, const float& lTh,const float& hTh, const float& cTh, Container& vthr){
  Data a;
  a.encode(strip,lTh,hTh,cTh);
  vthr.push_back(a);
}

SiStripThreshold::Data SiStripThreshold::getData(const uint16_t& strip, const Range& range) const {
  uint16_t estrip=(strip & sistrip::FirstThStripMask_)<<sistrip::FirstThStripShift_ | (63 & sistrip::HighThStripMask_);
  ContainerIterator p = std::upper_bound(range.first,range.second,estrip,SiStripThreshold::dataStrictWeakOrdering());
  if (p!=range.first){
    return *(--p);
  }
  else{
    throw cms::Exception("CorruptedData")
      << "[SiStripThreshold::getData] asking for data for a strip " << strip << " lower then the first stored strip " << p->getFirstStrip();
  }
}

void SiStripThreshold::allThresholds(std::vector<float> &lowThs, std::vector<float> &highThs, const Range& range) const  {
    ContainerIterator it = range.first;
    size_t strips = lowThs.size(); 
    assert(strips == highThs.size());
    while (it != range.second) {
        size_t firstStrip = it->getFirstStrip();
        //std::cout << "First strip is " << firstStrip << std::endl;
        float high = it->getHth(), low = it->getLth();
        //std::cout << "High is " << high << ", low is " << low << std::endl;
        ++it; // increment the pointer
        size_t lastStrip = (it == range.second ? strips : it->getFirstStrip());
        //std::cout << "Last strip is " << lastStrip << std::endl;
        if (lastStrip > strips) { 
            it = range.second;  // I should stop here,
            lastStrip = strips; // and fill only 'strips' strips
        }
        std::fill( & lowThs[firstStrip] , & lowThs[lastStrip] , low );
        std::fill( & highThs[firstStrip], & highThs[lastStrip], high );
    }
}    

void SiStripThreshold::printDebug(std::stringstream& ss) const{
  RegistryIterator rit=getRegistryVectorBegin(), erit=getRegistryVectorEnd();
  ContainerIterator it,eit;
  for(;rit!=erit;++rit){
    it=getDataVectorBegin()+rit->ibegin;
    eit=getDataVectorBegin()+rit->iend;
    ss << "\ndetid: " << rit->detid << " \t ";
    for(;it!=eit;++it){
      ss << "\n \t ";
      it->print(ss);
    }
  }
}

void SiStripThreshold::printSummary(std::stringstream& ss) const{
  RegistryIterator rit=getRegistryVectorBegin(), erit=getRegistryVectorEnd();
  ContainerIterator it,eit,itp;
  float meanLth, meanHth, meanCth; //mean value 
  float rmsLth, rmsHth, rmsCth; //rms value 
  float maxLth, maxHth, maxCth; //max value 
  float minLth, minHth, minCth; //min value 
  uint16_t n;
  uint16_t firstStrip,stripRange;
  for(;rit!=erit;++rit){
    it=getDataVectorBegin()+rit->ibegin;
    eit=getDataVectorBegin()+rit->iend;
    ss << "\ndetid: " << rit->detid << " \t ";

    meanLth=0; meanHth=0; meanCth=0; //mean value 
    rmsLth=0; rmsHth=0; rmsCth=0; //rms value 
    maxLth=0; maxHth=0; maxCth=0; //max value 
    minLth=10000; minHth=10000; minCth=10000; //min value 
    n=0;
    firstStrip=0;
    for(;it!=eit;++it){
      itp=it+1;
      firstStrip=it->getFirstStrip();
      if(itp!=eit)
        stripRange=(itp->getFirstStrip()-firstStrip);
      else 
        stripRange=firstStrip>511?768-firstStrip:512-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

      addToStat(it->getLth()   ,stripRange,meanLth,rmsLth,minLth,maxLth);
      addToStat(it->getHth()   ,stripRange,meanHth,rmsHth,minHth,maxHth);
      addToStat(it->getClusth(),stripRange,meanCth,rmsCth,minCth,maxCth);
      n+=stripRange;
    }
    meanLth/=n;
    meanHth/=n;
    meanCth/=n;
    rmsLth= sqrt(rmsLth/n-meanLth*meanLth);
    rmsHth= sqrt(rmsHth/n-meanHth*meanHth);
    rmsCth= sqrt(rmsCth/n-meanCth*meanCth);
    ss<< "\nn " << n << " \tmeanLth " << meanLth << " \t rmsLth " << rmsLth << " \t minLth " << minLth << " \t maxLth " << maxLth;  
    ss<< "\n\tmeanHth " << meanHth << " \t rmsHth " << rmsHth << " \t minHth " << minHth << " \t maxHth " << maxHth;  
    ss<< "\n\tmeanCth " << meanCth << " \t rmsCth " << rmsCth << " \t minCth " << minCth << " \t maxCth " << maxCth;  
  }
}

void SiStripThreshold::addToStat(float value, uint16_t& range, float& sum, float& sum2, float& min, float& max) const{
  sum+=value*range;
  sum2+=value*value*range;
  if(value<min)
    min=value;
  if(value>max)
    max=value;
}