SiStripFedZeroSuppression.cc

Go to the documentation of this file.

#include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripFedZeroSuppression.h"

#include "CondFormats/DataRecord/interface/SiStripNoisesRcd.h"
#include "CondFormats/SiStripObjects/interface/SiStripNoises.h"
#include "CondFormats/DataRecord/interface/SiStripThresholdRcd.h"
#include "CondFormats/SiStripObjects/interface/SiStripThreshold.h"

//#define DEBUG_SiStripZeroSuppression_
//#define ML_DEBUG 
using namespace std;

void SiStripFedZeroSuppression::init(const edm::EventSetup& es){
  uint32_t n_cache_id = es.get<SiStripNoisesRcd>().cacheIdentifier();
  uint32_t t_cache_id = es.get<SiStripThresholdRcd>().cacheIdentifier();

  if(n_cache_id != noise_cache_id) {
    es.get<SiStripNoisesRcd>().get( noiseHandle );
    noise_cache_id = n_cache_id;
  }
  if(t_cache_id != threshold_cache_id) {
    es.get<SiStripThresholdRcd>().get( thresholdHandle );
    threshold_cache_id = t_cache_id;
  }
}

void SiStripFedZeroSuppression::suppress(const std::vector<SiStripDigi>& in, std::vector<SiStripDigi>& selectedSignal,const uint32_t& detID){
  suppress(in, selectedSignal, detID, noiseHandle, thresholdHandle);
}

void SiStripFedZeroSuppression::suppress(const std::vector<SiStripDigi>& in, std::vector<SiStripDigi>& selectedSignal,const uint32_t& detID, 
                            edm::ESHandle<SiStripNoises> & noiseHandle,edm::ESHandle<SiStripThreshold> & thresholdHandle){

  int i;  
  int inSize = in.size();
  SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detID);
  SiStripThreshold::Range detThRange = thresholdHandle->getRange(detID);

  // reserving more than needed, but quicker than one at a time
  selectedSignal.clear();
  selectedSignal.reserve(inSize);
  for (i = 0; i < inSize; i++) {
    //Find adc values for neighbouring strips
    const uint32_t strip = (uint32_t) in[i].strip();
    
    adc   = in[i].adc();

    SiStripThreshold::Data thresholds=thresholdHandle->getData(strip,detThRange);
    theFEDlowThresh  = static_cast<int16_t>(thresholds.getLth()*noiseHandle->getNoiseFast(strip,detNoiseRange)+0.5);
    theFEDhighThresh = static_cast<int16_t>(thresholds.getHth()*noiseHandle->getNoiseFast(strip,detNoiseRange)+0.5);

    adcPrev  = -9999;
    adcNext  = -9999;
    adcPrev2 = -9999;
    adcNext2 = -9999;
    
    /*
      Since we are not running on 
      Raw data we need to initialize
      all the FED threshold
    */
    
    theNextFEDlowThresh  = theFEDlowThresh;
    theNext2FEDlowThresh = theFEDlowThresh;
    thePrevFEDlowThresh  = theFEDlowThresh;
    thePrev2FEDlowThresh = theFEDlowThresh;
    theNeighFEDlowThresh = theFEDlowThresh;

    theNextFEDhighThresh  = theFEDhighThresh;
    thePrevFEDhighThresh  = theFEDhighThresh;
    theNeighFEDhighThresh = theFEDhighThresh;

    if ( ((strip)%128) == 127){ 
      adcNext = 0;
      theNextFEDlowThresh  = 9999;
      theNextFEDhighThresh = 9999;
    }else if (i + 1 < inSize && in[i+1].strip() == strip + 1) {
      adcNext = in[i+1].adc();
      SiStripThreshold::Data thresholds_1=thresholdHandle->getData(strip+1,detThRange);
      theNextFEDlowThresh  = static_cast<int16_t>(thresholds_1.getLth()*noiseHandle->getNoiseFast(strip+1,detNoiseRange)+0.5);
      theNextFEDhighThresh = static_cast<int16_t>(thresholds_1.getHth()*noiseHandle->getNoiseFast(strip+1,detNoiseRange)+0.5);
      if ( ((strip)%128) == 126){ 
    adcNext2 = 0;
    theNext2FEDlowThresh  = 9999;
      }else if (i + 2 < inSize && in[i+2].strip() == strip + 2) {
    adcNext2 = in[i+2].adc();
    theNext2FEDlowThresh  = static_cast<int16_t>(thresholdHandle->getData(strip+2,detThRange).getLth()*noiseHandle->getNoiseFast(strip+2,detNoiseRange)+0.5);
      }
    }

    if ( ((strip)%128) == 0){   
      adcPrev = 0;
      thePrevFEDlowThresh  = 9999;
      thePrevFEDhighThresh = 9999;
    }else if (i - 1 >= 0 && in[i-1].strip() == strip - 1) {
      adcPrev = in[i-1].adc();
      SiStripThreshold::Data thresholds_1=thresholdHandle->getData(strip-1,detThRange);
      thePrevFEDlowThresh  = static_cast<int16_t>(thresholds_1.getLth()*noiseHandle->getNoiseFast(strip-1,detNoiseRange)+0.5);
      thePrevFEDhighThresh = static_cast<int16_t>(thresholds_1.getHth()*noiseHandle->getNoiseFast(strip-1,detNoiseRange)+0.5);
      if ( ((strip)%128) == 1){
    adcPrev2 = 0; 
    thePrev2FEDlowThresh  = 9999;
      }else if (i - 2 >= 0 && in[i-2].strip() == strip - 2) {
    adcPrev2 = in[i-2].adc();
    thePrev2FEDlowThresh  = static_cast<int16_t>(thresholdHandle->getData(strip-2,detThRange).getLth()*noiseHandle->getNoiseFast(strip-2,detNoiseRange)+0.5);
      }
    }

    if ( adcNext <= adcPrev){
      adcMaxNeigh = adcPrev;
      theNeighFEDlowThresh  = thePrevFEDlowThresh;
      theNeighFEDhighThresh = thePrevFEDhighThresh;
    } else {
      adcMaxNeigh = adcNext;
      theNeighFEDlowThresh  = theNextFEDlowThresh;
      theNeighFEDhighThresh = theNextFEDhighThresh;
    }
    
    if (IsAValidDigi()){
      selectedSignal.push_back(SiStripDigi(strip, adc));
    }
  }
}

void SiStripFedZeroSuppression::suppress(const edm::DetSet<SiStripRawDigi>& in, edm::DetSet<SiStripDigi>& out)
{
  const uint32_t detID = out.id;
  SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detID);
  SiStripThreshold::Range detThRange = thresholdHandle->getRange(detID);
#ifdef DEBUG_SiStripZeroSuppression_
  if (edm::isDebugEnabled())
    LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress] Zero suppression on edm::DetSet<SiStripRawDigi>: detID " << detID << " size = " << in.data.size();
#endif
  edm::DetSet<SiStripRawDigi>::const_iterator in_iter=in.data.begin();
  for (;in_iter!=in.data.end();in_iter++){

    const uint32_t strip = (uint32_t) (in_iter-in.data.begin());

#ifdef DEBUG_SiStripZeroSuppression_
    if (edm::isDebugEnabled())
      LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress] detID= " << detID << " strip= " <<  strip << "  adc= " << in_iter->adc();
#endif    
    adc   = in_iter->adc();

    SiStripThreshold::Data thresholds=thresholdHandle->getData(strip,detThRange);
    theFEDlowThresh  = static_cast<int16_t>(thresholds.getLth()*noiseHandle->getNoiseFast(strip,detNoiseRange)+0.5);
    theFEDhighThresh = static_cast<int16_t>(thresholds.getHth()*noiseHandle->getNoiseFast(strip,detNoiseRange)+0.5);

    adcPrev = -9999;
    adcNext = -9999;
    /*
      If a strip is the last one on the chip
      set its next neighbor's thresholds to infinity
      because the FED does not merge clusters across
      chip boundaries right now
    */
    if ( strip%128 == 127 ) { 
      adcNext = 0;      
      theNextFEDlowThresh  = 9999;
      theNextFEDhighThresh = 9999;
    }
    else {
      adcNext = (in_iter+1)->adc();
      SiStripThreshold::Data thresholds_1=thresholdHandle->getData(strip+1,detThRange);
      theNextFEDlowThresh  = static_cast<int16_t>(thresholds_1.getLth()*noiseHandle->getNoiseFast(strip+1,detNoiseRange)+0.5);
      theNextFEDhighThresh = static_cast<int16_t>(thresholds_1.getHth()*noiseHandle->getNoiseFast(strip+1,detNoiseRange)+0.5);
    }
    /*
      Similarily, for the first strip 
      on a chip
    */
    if ( strip%128 == 0 ) {
      adcPrev = 0;
      thePrevFEDlowThresh  = 9999;
      thePrevFEDhighThresh = 9999;   
    }
    else {
      adcPrev = (in_iter-1)->adc();
      SiStripThreshold::Data thresholds_1=thresholdHandle->getData(strip-1,detThRange);
      thePrevFEDlowThresh  = static_cast<int16_t>(thresholds_1.getLth()*noiseHandle->getNoiseFast(strip-1,detNoiseRange)+0.5);
      thePrevFEDhighThresh = static_cast<int16_t>(thresholds_1.getHth()*noiseHandle->getNoiseFast(strip-1,detNoiseRange)+0.5);
    }
    if ( adcNext < adcPrev){
      adcMaxNeigh = adcPrev;
      theNeighFEDlowThresh  = thePrevFEDlowThresh;
      theNeighFEDhighThresh = thePrevFEDhighThresh;
    } else {
      adcMaxNeigh = adcNext;
      theNeighFEDlowThresh  = theNextFEDlowThresh;
      theNeighFEDhighThresh = theNextFEDhighThresh;
    }
    
    //Find adc values for next neighbouring strips
    adcPrev2 = -9999;
    adcNext2 = -9999;
    thePrev2FEDlowThresh  = 1;
    theNext2FEDlowThresh  = 1;
    if ( strip%128 >= 126 ) {
      adcNext2 = 0;
      theNext2FEDlowThresh  = 9999;
    }
    else if ( strip%128 < 126 ) {
      adcNext2 = (in_iter+2)->adc();
     theNext2FEDlowThresh  = static_cast<int16_t>(thresholdHandle->getData(strip+2,detThRange).getLth()*noiseHandle->getNoiseFast(strip+2,detNoiseRange)+0.5);
    }
    if ( strip%128 <= 1 ) {
      adcPrev2 = 0; 
      thePrev2FEDlowThresh  = 9999;
    }
    else if ( strip%128 > 1 ) {
      adcPrev2 = (in_iter-2)->adc();
     thePrev2FEDlowThresh  = static_cast<int16_t>(thresholdHandle->getData(strip-2,detThRange).getLth()*noiseHandle->getNoiseFast(strip-2,detNoiseRange)+0.5);
    }
    //GB 23/6/08: truncation should be done at the very beginning
    if (IsAValidDigi())
      out.data.push_back(SiStripDigi(strip, truncate(in_iter->adc())));
  }
}

void SiStripFedZeroSuppression::fillThresholds_(const uint32_t detID, size_t size) {
  SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detID);
  SiStripThreshold::Range detThRange = thresholdHandle->getRange(detID);

  if (highThr_.size() != size) { 
    highThr_.resize(size); 
    lowThr_.resize(size);
    noises_.resize(size);
    highThrSN_.resize(size);
    lowThrSN_.resize(size);
  }
  
  noiseHandle->allNoises(noises_, detNoiseRange);
  thresholdHandle->allThresholds(lowThrSN_, highThrSN_, detThRange); // thresholds as S/N
  for (size_t strip = 0; strip < size; ++strip) {
    float noise     = noises_[strip];
    //  uncomment line below to check bluk noise decoding
    //assert( noise == noiseHandle->getNoiseFast(strip,detNoiseRange) ); 
    highThr_[strip] = static_cast<int16_t>(highThrSN_[strip]*noise+0.5+1e-6);
    lowThr_[strip]  = static_cast<int16_t>( lowThrSN_[strip]*noise+0.5+1e-6);
    // Note: it's a bit wierd, but there are some cases for which 'highThrSN_[strip]*noise' is an exact integer
    //   but due to roundoffs it gets rounded to the integer below if. 
    //   Apparently the optimized code inlines differently and this changes the roundoff.
    //   The +1e-6 fixes the problem.   [GPetruc]
  } 
}


void SiStripFedZeroSuppression::suppress(const std::vector<int16_t>& in, const uint16_t& firstAPV,  edm::DetSet<SiStripDigi>& out){

  const uint32_t detID = out.id;
  size_t size = in.size();
#ifdef DEBUG_SiStripZeroSuppression_
  if (edm::isDebugEnabled())
    LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress] Zero suppression on std::vector<int16_t>: detID " << detID << " size = " << in.size();
#endif

  fillThresholds_(detID, size+firstAPV*128); // want to decouple this from the other cost


  std::vector<int16_t>::const_iterator in_iter=in.begin();
  uint16_t strip = firstAPV*128;
  for (; strip < size+firstAPV*128; ++strip, ++in_iter){

    size_t strip_mod_128 = strip & 127;
#ifdef DEBUG_SiStripZeroSuppression_
    if (edm::isDebugEnabled())
      LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress]  detID= " << detID << " strip= " <<  strip << "  adc= " << *in_iter;
#endif    
    adc   = *in_iter;

    theFEDlowThresh  = lowThr_[strip];
    theFEDhighThresh = highThr_[strip];

    //Find adc values for neighbouring strips

     /*
      If a strip is the last one on the chip
      set its next neighbor's thresholds to infinity
      because the FED does not merge clusters across
      chip boundaries right now
    */

    //adcPrev = -9999;  // useless, they are set
    //adcNext = -9999;  // in the next lines in any case
    if ( strip_mod_128 == 127 ) {
      adcNext = 0;      
      theNextFEDlowThresh  = 9999;
      theNextFEDhighThresh = 9999;
    } else {
      adcNext = *(in_iter+1);
      theNextFEDlowThresh  = lowThr_[strip+1];
      theNextFEDhighThresh = highThr_[strip+1];
    }
    
    /*
      Similarily, for the first strip 
      on a chip
    */
    if ( strip_mod_128 == 0 ) {
      adcPrev = 0;
      thePrevFEDlowThresh  = 9999;
      thePrevFEDhighThresh = 9999;   
    } else {
      adcPrev = *(in_iter-1);
      thePrevFEDlowThresh  = lowThr_[strip-1];
      thePrevFEDhighThresh = highThr_[strip-1];
    }

    if ( adcNext < adcPrev){
      adcMaxNeigh = adcPrev;
      theNeighFEDlowThresh  = thePrevFEDlowThresh;
      theNeighFEDhighThresh = thePrevFEDhighThresh;
    } else {
      adcMaxNeigh = adcNext;
      theNeighFEDlowThresh  = theNextFEDlowThresh;
      theNeighFEDhighThresh = theNextFEDhighThresh;
    }
    
    //Find adc values for next neighbouring strips
    //adcPrev2 = -9999;           //
    //adcNext2 = -9999;           // useless to set them here
    //thePrev2FEDlowThresh  = 1;  // they are overwritten always in the next 8 lines
    //theNext2FEDlowThresh  = 1;  //
    if ( strip_mod_128 >=126 ) {
      adcNext2 = 0;
      theNext2FEDlowThresh  = 9999;
    //} else if ( strip_mod_128 < 126 ) { // if it's not >= then is <, no need to "if" again
    } else {
      adcNext2 = *(in_iter+2);
      theNext2FEDlowThresh  = lowThr_[strip+2];
    }
    if ( strip_mod_128 <= 1 ) {
      adcPrev2 = 0; 
      thePrev2FEDlowThresh  = 9999;
    //} else if ( strip_mod_128 > 1 ) { // same as above
    } else {
      adcPrev2 = *(in_iter-2);
      thePrev2FEDlowThresh  = lowThr_[strip-2];;
    }
    
    if (IsAValidDigi()){
#ifdef DEBUG_SiStripZeroSuppression_
      if (edm::isDebugEnabled())
    LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress] DetId " << out.id << " strip " << strip << " adc " << *in_iter << " digiCollection size " << out.data.size() ;
#endif            
      //GB 23/6/08: truncation should be done at the very beginning
      out.push_back(SiStripDigi(strip, (*in_iter<0 ? 0 : truncate( *in_iter ) )));
    }
  }
}


bool SiStripFedZeroSuppression::IsAValidDigi()
{

#ifdef DEBUG_SiStripZeroSuppression_


  if (edm::isDebugEnabled()){

    LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress] "
                       << "\n\t adc " << adc 
                       << "\n\t adcPrev " << adcPrev
                       << "\n\t adcNext " << adcNext 
                       << "\n\t adcMaxNeigh " << adcMaxNeigh 
                       << "\n\t adcPrev2 " << adcPrev2 
                       << "\n\t adcNext2 " << adcNext2 
                       <<std::endl;
  
    LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress] "
                       << "\n\t theFEDlowThresh " <<  theFEDlowThresh 
                       << "\n\t theFEDhighThresh " << theFEDhighThresh 
                       << "\n\t thePrevFEDlowThresh " <<  thePrevFEDlowThresh 
                       << "\n\t thePrevFEDhighThresh " << thePrevFEDhighThresh 
                       << "\n\t theNextFEDlowThresh " <<  theNextFEDlowThresh 
                       << "\n\t theNextFEDhighThresh " << theNextFEDhighThresh 
                       << "\n\t theNeighFEDlowThresh " <<  theNeighFEDlowThresh 
                       << "\n\t theNeighFEDhighThresh " << theNeighFEDhighThresh 
                       << "\n\t thePrev2FEDlowThresh " <<  thePrev2FEDlowThresh 
                       << "\n\t theNext2FEDlowThresh " <<  theNext2FEDlowThresh
                       <<std::endl;
  }
#endif  
  // Decide if this strip should be accepted.
  bool accept = false;
  switch (theFEDalgorithm) {
  case 1:
    accept = (adc >= theFEDlowThresh);
    break;
  case 2:
    accept = (adc >= theFEDhighThresh || (adc >= theFEDlowThresh &&
                      adcMaxNeigh >= theNeighFEDlowThresh));
    break;
  case 3:
    accept = (adc >= theFEDhighThresh || (adc >= theFEDlowThresh &&
                      adcMaxNeigh >= theNeighFEDhighThresh));
    break;
  case 4:
    accept = (
          (adc >= theFEDhighThresh)            //Test for adc>highThresh (same as algorithm 2)
          ||
          (
           (adc >= theFEDlowThresh)            //Test for adc>lowThresh, with neighbour adc>lowThresh (same as algorithm 2)
           &&
           (adcMaxNeigh >= theNeighFEDlowThresh)
           ) 
          ||
          (
           (adc < theFEDlowThresh)        //Test for adc<lowThresh
           &&     
           (
        (
         (adcPrev  >= thePrevFEDhighThresh)    //with both neighbours>highThresh
         &&
         (adcNext  >= theNextFEDhighThresh)
         ) 
        ||
        (
         (adcPrev  >= thePrevFEDhighThresh)    //OR with previous neighbour>highThresh and
         &&
         (adcNext  >= theNextFEDlowThresh)     //both the next neighbours>lowThresh
         &&
         (adcNext2 >= theNext2FEDlowThresh)
         )  
        ||
        (
         (adcNext  >= theNextFEDhighThresh)    //OR with next neighbour>highThresh and
         &&
         (adcPrev  >= thePrevFEDlowThresh)     //both the previous neighbours>lowThresh
         &&
         (adcPrev2 >= thePrev2FEDlowThresh)
         )  
        ||
        (
         (adcNext  >= theNextFEDlowThresh)     //OR with both next neighbours>lowThresh and
         &&
         (adcNext2 >= theNext2FEDlowThresh)   //both the previous neighbours>lowThresh
         &&
         (adcPrev  >= thePrevFEDlowThresh)  
         &&
         (adcPrev2 >= thePrev2FEDlowThresh)
         )
        )
           )
          );
    break;
  case 5:
    accept = adc > 0;
    break;

  }
  return accept;
}