SpikeAndDoubleSpikeCleaner Class Reference

#include <SpikeAndDoubleSpikeCleaner.h>

Inheritance diagram for SpikeAndDoubleSpikeCleaner:

Classes
struct	spike_cleaning

Public Member Functions
void	clean (const edm::Handle< reco::PFRecHitCollection > &input, std::vector< bool > &mask) override
SpikeAndDoubleSpikeCleaner &	operator= (const SpikeAndDoubleSpikeCleaner &)=delete
	SpikeAndDoubleSpikeCleaner (const edm::ParameterSet &conf)
	SpikeAndDoubleSpikeCleaner (const SpikeAndDoubleSpikeCleaner &)=delete
Public Member Functions inherited from RecHitTopologicalCleanerBase
const std::string &	name () const
RecHitTopologicalCleanerBase &	operator= (const RecHitTopologicalCleanerBase &)=delete
	RecHitTopologicalCleanerBase (const edm::ParameterSet &conf)
	RecHitTopologicalCleanerBase (const RecHitTopologicalCleanerBase &)=delete
virtual	~RecHitTopologicalCleanerBase ()=default

Private Attributes
const std::unordered_map< std::string, int >	_layerMap
std::unordered_map< int, spike_cleaning >	_thresholds

Detailed Description

Definition at line 9 of file SpikeAndDoubleSpikeCleaner.h.

Constructor & Destructor Documentation

SpikeAndDoubleSpikeCleaner::SpikeAndDoubleSpikeCleaner

(

const edm::ParameterSet &

conf

)

Definition at line 62 of file SpikeAndDoubleSpikeCleaner.cc.

References SpikeAndDoubleSpikeCleaner::spike_cleaning::_doubleSpikeS6S2, SpikeAndDoubleSpikeCleaner::spike_cleaning::_doubleSpikeThresh, SpikeAndDoubleSpikeCleaner::spike_cleaning::_eneThreshMod, SpikeAndDoubleSpikeCleaner::spike_cleaning::_fracThreshMod, _layerMap, SpikeAndDoubleSpikeCleaner::spike_cleaning::_minS4S1_a, SpikeAndDoubleSpikeCleaner::spike_cleaning::_minS4S1_b, SpikeAndDoubleSpikeCleaner::spike_cleaning::_singleSpikeThresh, _thresholds, clean(), PFLayer::ECAL_BARREL, PFLayer::ECAL_ENDCAP, mps_splice::entry, Exception, PFLayer::HCAL_BARREL1, PFLayer::HCAL_BARREL2, PFLayer::HCAL_ENDCAP, PFLayer::HF_EM, PFLayer::HF_HAD, info(), createfilelist::int, PFLayer::NONE, PFLayer::PS1, PFLayer::PS2, muonDTDigis_cfi::pset, AlCaHLTBitMon_QueryRunRegistry::string, and particleFlowRecHitECAL_cfi::thresholds.

                                                        :
    RecHitTopologicalCleanerBase(conf),
    _layerMap({ {"PS2",(int)PFLayer::PS2},
          {"PS1",(int)PFLayer::PS1},
          {"ECAL_ENDCAP",(int)PFLayer::ECAL_ENDCAP},
          {"ECAL_BARREL",(int)PFLayer::ECAL_BARREL},
          {"NONE",(int)PFLayer::NONE},
          {"HCAL_BARREL1",(int)PFLayer::HCAL_BARREL1},
          {"HCAL_BARREL2_RING0",(int)PFLayer::HCAL_BARREL2},
              // hack to deal with ring1 in HO 
          {"HCAL_BARREL2_RING1",100*(int)PFLayer::HCAL_BARREL2}, 
          {"HCAL_ENDCAP",(int)PFLayer::HCAL_ENDCAP},
          {"HF_EM",(int)PFLayer::HF_EM},
          {"HF_HAD",(int)PFLayer::HF_HAD} }) {
  const std::vector<edm::ParameterSet>& thresholds =
    conf.getParameterSetVector("cleaningByDetector");
  for( const auto& pset : thresholds ) {
    spike_cleaning info;
    const std::string& det = pset.getParameter<std::string>("detector");
    info._minS4S1_a = pset.getParameter<double>("minS4S1_a");    
    info._minS4S1_b = pset.getParameter<double>("minS4S1_b");    
    info._doubleSpikeS6S2 = pset.getParameter<double>("doubleSpikeS6S2");
    info._eneThreshMod = pset.getParameter<double>("energyThresholdModifier");
    info._fracThreshMod = 
      pset.getParameter<double>("fractionThresholdModifier");
    info._doubleSpikeThresh = pset.getParameter<double>("doubleSpikeThresh");
    info._singleSpikeThresh = pset.getParameter<double>("singleSpikeThresh");
    auto entry = _layerMap.find(det);
    if( entry == _layerMap.end() ) {
      throw cms::Exception("InvalidDetectorLayer")
    << "Detector layer : " << det << " is not in the list of recognized"
    << " detector layers!";
    }
    _thresholds.emplace(_layerMap.find(det)->second,info);
  }  
}

SpikeAndDoubleSpikeCleaner::SpikeAndDoubleSpikeCleaner ( const SpikeAndDoubleSpikeCleaner &  )

delete

Member Function Documentation

void SpikeAndDoubleSpikeCleaner::clean ( const edm::Handle< reco::PFRecHitCollection > &  input, std::vector< bool > &  mask  )

overridevirtual

Implements RecHitTopologicalCleanerBase.

Definition at line 101 of file SpikeAndDoubleSpikeCleaner.cc.

References SpikeAndDoubleSpikeCleaner::spike_cleaning::_doubleSpikeS6S2, SpikeAndDoubleSpikeCleaner::spike_cleaning::_doubleSpikeThresh, SpikeAndDoubleSpikeCleaner::spike_cleaning::_eneThreshMod, SpikeAndDoubleSpikeCleaner::spike_cleaning::_fracThreshMod, SpikeAndDoubleSpikeCleaner::spike_cleaning::_minS4S1_a, SpikeAndDoubleSpikeCleaner::spike_cleaning::_minS4S1_b, SpikeAndDoubleSpikeCleaner::spike_cleaning::_singleSpikeThresh, _thresholds, funct::abs(), AlCaHLTBitMon_QueryRunRegistry::comp, reco::PFRecHit::detId(), PFLayer::ECAL_BARREL, reco::PFRecHit::energy(), RhoEtaPhi::eta(), PFLayer::HCAL_BARREL2, HcalForward, PFLayer::HF_EM, PFLayer::HF_HAD, hfClusterShapes_cfi::hits, mps_fire::i, training_settings::idx, HcalDetId::ieta(), recoMuon::in, input, createfilelist::int, HcalDetId::iphi(), gen::k, reco::PFRecHit::layer(), min(), reco::PFRecHit::neighbours4(), RhoEtaPhi::phi(), reco::PFRecHit::positionREP(), and DetId::rawId().

Referenced by SpikeAndDoubleSpikeCleaner().

                              {
  //need to run over energy sorted rechits
  auto const & hits = *input;
  std::vector<unsigned > ordered_hits(hits.size());
  for( unsigned i = 0; i < hits.size(); ++i ) ordered_hits[i]=i;
  std::sort(ordered_hits.begin(),ordered_hits.end(),[&](unsigned i, unsigned j) { return hits[i].energy()>hits[j].energy();});
    
  for( const auto& idx : ordered_hits ) {
    const unsigned i = idx;
    if( !mask[i] ) continue; // don't need to re-mask things :-)
    const reco::PFRecHit& rechit = hits[i];
    int hitlayer = (int)rechit.layer();
    if( hitlayer == PFLayer::HCAL_BARREL2 && 
    std::abs(rechit.positionREP().eta()) > 0.34 ) {
      hitlayer *= 100;
    }    
    const spike_cleaning& clean = _thresholds.find(hitlayer)->second;    
    if( rechit.energy() < clean._singleSpikeThresh ) continue;

    //Fix needed for HF.  Here, we find the (up to) five companion rechits 
    //to work in conjunction with the neighbours4() implementation below for the full HF surrounding energy
    float compsumE = 0.0;
    if ((hitlayer==PFLayer::HF_EM || hitlayer==PFLayer::HF_HAD)) {
      int comp = 1;
      if (hitlayer==PFLayer::HF_EM) comp = 2;
      const HcalDetId& detid = (HcalDetId)rechit.detId();
      int heta = detid.ieta();
      int hphi = detid.iphi();
      
      //At eta>39, phi segmentation changes 
      int predphi =2;
      if (std::abs(heta)>39) predphi=4;

      int curphiL = hphi-predphi;
      int curphiH = hphi+predphi;

      //HcalDetId valid phi range (0-72)
      while (curphiL<0) curphiL+=72;
      while (curphiH>72) curphiH-=72;

      std::pair<std::vector<int>, std::vector<int>>  phietas({heta,heta+1,heta-1,heta,heta},{hphi,hphi,hphi,curphiL,curphiH});

      std::vector<uint32_t> rawDetIds;
      for(unsigned in=0;in<phietas.first.size();in++) {
        HcalDetId tempID (HcalForward, phietas.first[in], phietas.second[in], comp);
        rawDetIds.push_back(tempID.rawId());
      }

      for( const auto& jdx : ordered_hits ) {
        const unsigned j = jdx;
        const reco::PFRecHit& matchrechit = hits[j];
        for( const auto& iID : rawDetIds ) if (iID==matchrechit.detId())compsumE+=matchrechit.energy();  
      }
    }
    //End of fix needed for HF

    const double rhenergy = rechit.energy();
    // single spike cleaning
    auto const & neighbours4 = rechit.neighbours4();
    double surroundingEnergy = compsumE;
    for( auto k : neighbours4 ) {
      if( !mask[k] ) continue;
      auto const & neighbour = hits[k];
      const double sum = neighbour.energy(); //energyUp is just rechit energy?
      surroundingEnergy += sum;
    }    
    //   wannaBeSeed.energyUp()/wannaBeSeed.energy() : 1.;
    // Fraction 1 is the balance between the hit and its neighbours 
    // from both layers
    const double fraction1 = surroundingEnergy/rhenergy;    
    // removed spurious comments from old pfcluster algo... 
    // look there if you want more history
    const double f1Cut = ( clean._minS4S1_a*std::log10(rechit.energy()) + 
               clean._minS4S1_b );
    if( fraction1 < f1Cut ) { 
      const double eta = rechit.positionREP().eta();
      const double aeta = std::abs(eta);
      const double phi = rechit.positionREP().phi();
      std::pair<double,double> dcr = dCrack(phi,eta);
      const double dcrmin = ( rechit.layer() ==  PFLayer::ECAL_BARREL ? 
                  std::min(dcr.first,dcr.second):
                  dcr.second );
      if( aeta < 5.0 && 
      ( (aeta < 2.85 && dcrmin > 1.0) || 
        (rhenergy > clean._eneThreshMod*clean._singleSpikeThresh && 
         fraction1 < f1Cut/clean._fracThreshMod ) ) ) { 
    mask[i] = false;
      }
    }//if initial fraction cut (single spike)
    // double spike removal
    if( mask[i] && rhenergy > clean._doubleSpikeThresh ) {
      //Determine energy surrounding the seed and the most energetic neighbour
      double surroundingEnergyi = 0.0;      
      double enmax = -999.0;
      unsigned int  mostEnergeticNeighbour=0;
      auto const & neighbours4i = rechit.neighbours4();
      for( auto k : neighbours4i ) {
    if( !mask[k] ) continue;
        auto const & neighbour = hits[k];
    const double nenergy = neighbour.energy();
    surroundingEnergyi += nenergy;
    if( nenergy > enmax ) {
      enmax = nenergy;
      mostEnergeticNeighbour = k;
    }
      }
      // is there an energetic neighbour
      if( enmax > 0.0 ) {
    double surroundingEnergyj = 0.0;
    auto const & neighbours4j = 
      hits[mostEnergeticNeighbour].neighbours4();
    for(auto k : neighbours4j ) {
      //if( !mask[k] &&  k != i) continue; // leave out?
      surroundingEnergyj += hits[k].energy();
    }
    // The energy surrounding the double spike candidate 
    const double surroundingEnergyFraction = 
      (surroundingEnergyi+surroundingEnergyj) / 
      (rechit.energy()+hits[mostEnergeticNeighbour].energy()) - 1.;
    if ( surroundingEnergyFraction < clean._doubleSpikeS6S2 ) { 
      const double eta = rechit.positionREP().eta();
      const double aeta = std::abs(eta);
      const double phi = rechit.positionREP().phi();
      std::pair<double,double> dcr = dCrack(phi,eta);
      const double dcrmin = ( rechit.layer() == PFLayer::ECAL_BARREL ? 
                  std::min(dcr.first,dcr.second):
                  dcr.second );
      if( aeta < 5.0 && 
          ( (aeta < 2.85 && dcrmin > 1.0) || 
        (rhenergy > clean._eneThreshMod*clean._doubleSpikeThresh && 
         surroundingEnergyFraction < clean._doubleSpikeS6S2/clean._fracThreshMod ) 
        ) ) {       
        mask[i] = false;
        mask[mostEnergeticNeighbour] = false;
      }
    }
      } // was there an energetic neighbour ? 
    }// if double spike thresh
  } // rechit loop
}

SpikeAndDoubleSpikeCleaner& SpikeAndDoubleSpikeCleaner::operator= ( const SpikeAndDoubleSpikeCleaner &  )

delete

Member Data Documentation

const std::unordered_map<std::string,int> SpikeAndDoubleSpikeCleaner::_layerMap

private

Definition at line 30 of file SpikeAndDoubleSpikeCleaner.h.

Referenced by SpikeAndDoubleSpikeCleaner().

std::unordered_map<int,spike_cleaning> SpikeAndDoubleSpikeCleaner::_thresholds

private

Definition at line 31 of file SpikeAndDoubleSpikeCleaner.h.

Referenced by clean(), and SpikeAndDoubleSpikeCleaner().