StripClusterizerAlgorithm.cc

Go to the documentation of this file.

#include "RecoLocalTracker/SiStripClusterizer/interface/StripClusterizerAlgorithm.h"

#include "FWCore/Framework/interface/EventSetup.h"
#include "CondFormats/DataRecord/interface/SiStripNoisesRcd.h"
#include "CalibTracker/Records/interface/SiStripGainRcd.h"
#include "CalibTracker/Records/interface/SiStripQualityRcd.h"
#include "DataFormats/SiStripDigi/interface/SiStripDigi.h"
#include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
#include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

// #define VIDEBUG

#ifdef VIDEBUG
#define COUT std::cout
#else
#define COUT LogDebug("StripClusterizerAlgorithm")
#endif

#include <string>
#include <algorithm>
#include <cassert>

void StripClusterizerAlgorithm::
initialize(const edm::EventSetup& es) {
  uint32_t n_cache_id = es.get<SiStripNoisesRcd>().cacheIdentifier();
  uint32_t g_cache_id = es.get<SiStripGainRcd>().cacheIdentifier();
  uint32_t q_cache_id = es.get<SiStripQualityRcd>().cacheIdentifier();

  bool mod=false;
  if(g_cache_id != gain_cache_id) {
    es.get<SiStripGainRcd>().get( gainHandle );
    gain_cache_id = g_cache_id;
    mod=true;
  }
  if(n_cache_id != noise_cache_id) {
    es.get<SiStripNoisesRcd>().get( noiseHandle );
    noise_cache_id = n_cache_id;
    mod=true;
  }
  if(q_cache_id != quality_cache_id) {
    es.get<SiStripQualityRcd>().get( qualityLabel, qualityHandle );
    quality_cache_id = q_cache_id;
    mod=true;
  }

  if (mod) { 
    // redo indexing!
    theCabling = qualityHandle->cabling();
    assert(theCabling); 
    auto const & conn = cabling()->connected();
    COUT << "cabling " << conn.size() << std::endl;
    detIds.clear();
    detIds.reserve(conn.size());
    for (auto const & c : conn) { if (!isModuleBad(c.first)) detIds.push_back(c.first);}
    indices.clear();
    indices.resize(detIds.size());
    COUT << "good detIds " << detIds.size() << std::endl;

    if (detIds.empty()) return;

    {
      connections.clear();
      connections.resize(detIds.size());
      //connections (slow, not a big deal)
      auto const & conns = cabling()->getDetCabling();
      for (auto i=0U; i<detIds.size(); ++i) {
    auto c = conns.find(detIds[i]);
    if (c!=conns.end()) connections[i]=(*c).second;
      }
      

    }


    { // quality
      std::vector<uint32_t> dum; qualityHandle->getDetIds(dum); 
      assert(dum.size()<invalidI);
      unsigned short j=0, i=0;
      while (i<dum.size() && j<detIds.size()) {
    if (dum[i]<detIds[j]) ++i;
    else if (detIds[j]<dum[i]) {indices[j].qi=invalidI; ++j;}
    else {
      indices[j].qi=i; ++i; ++j;
    }
      }
      unsigned int  nn=0;
      for(auto k=0U; k<detIds.size();++k) { if (indices[k].qi<invalidI) {++nn; assert(dum[indices[k].qi]==detIds[k]);}}
      assert(nn<=dum.size());
      COUT << "quality " << dum.size() << " " <<nn<< std::endl;
    }
    { //noise
      std::vector<uint32_t> dum; noiseHandle->getDetIds(dum); 
      assert(dum.size()<invalidI);
      unsigned short j=0, i=0;
      while (i<dum.size() && j<detIds.size()) {
    if (dum[i]<detIds[j]) ++i;
    else if (detIds[j]<dum[i]) {indices[j].ni=invalidI; ++j;}
    else {
      indices[j].ni=i; ++i; ++j;
    }
      }
      unsigned int  nn=0;
      for(auto k=0U; k<detIds.size();++k) { if (indices[k].ni<invalidI) {++nn; assert(dum[indices[k].ni]==detIds[k]);}}
      assert(nn<=dum.size());
      COUT << "noise " << dum.size() << " " <<nn<< std::endl;
    }
    { //gain
      std::vector<uint32_t> dum; gainHandle->getDetIds(dum); 
      assert(dum.size()<invalidI);
      unsigned short j=0, i=0;
      while (i<dum.size() && j<detIds.size()) {
    if (dum[i]<detIds[j]) ++i;
    else if (detIds[j]<dum[i]) {indices[j].gi=invalidI; ++j;}
    else {
      indices[j].gi=i; ++i; ++j;
    }
      }
      unsigned int  nn=0;
      for(auto k=0U; k<detIds.size();++k) { if (indices[k].gi<invalidI) {++nn; assert(dum[indices[k].gi]==detIds[k]);}}
      assert(nn<=dum.size());
      COUT << "gain " << dum.size() << " " <<nn<< std::endl;
    }
  }
  
}
  
StripClusterizerAlgorithm::Det
StripClusterizerAlgorithm::
findDetId(const uint32_t id) const {
  auto b = detIds.begin();
  auto e = detIds.end();
  auto p = std::lower_bound(b,e,id);
  if (p==e || id!=(*p)) {
#ifdef NOT_ON_MONTECARLO
    edm::LogWarning("StripClusterizerAlgorithm") 
      <<"id " << id << " not connected. this is impossible on data "
      << "old id " << detId << std::endl;
#endif
    return Det();
  }
  Det det;
  det.ind = p-detIds.begin();

  det.detId = id;
  det.noiseRange = noiseHandle->getRangeByPos(indices[det.ind].ni);
  det.gainRange = gainHandle->getRangeByPos(indices[det.ind].gi);
  det.qualityRange = qualityHandle->getRangeByPos(indices[det.ind].qi);
  det.quality =   qualityHandle.product();

#ifdef EDM_ML_DEBUG
  assert(detIds[det.ind]==det.detId); 
  auto oldg =  gainHandle->getRange(id);
  assert(oldg==det.gainRange);
  auto oldn = noiseHandle->getRange(id);
  assert(oldn==det.noiseRange);
  auto oldq = qualityHandle->getRange(id);
  assert(oldq==det.qualityRange);
#endif
#ifdef EDM_ML_DEBUG
  assert(isModuleUsable( id ));
#endif
  return det;
}

void StripClusterizerAlgorithm::clusterize(const   edm::DetSetVector<SiStripDigi>& input,  output_t& output) const {clusterize_(input, output);}
void StripClusterizerAlgorithm::clusterize(const edmNew::DetSetVector<SiStripDigi>& input, output_t& output) const {clusterize_(input, output);}

StripClusterizerAlgorithm::
InvalidChargeException::InvalidChargeException(const SiStripDigi& digi)
  : cms::Exception("Invalid Charge") {
  std::stringstream s;
  s << "Digi charge of " << digi.adc() << " ADC "
    << "is out of range on strip " << digi.strip() << ".  ";
  this->append(s.str());
}