SpikeAndDoubleSpikeCleaner Class Reference

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, edm::ConsumesCollector &cc)
	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, edm::ConsumesCollector &cc)
	RecHitTopologicalCleanerBase (const RecHitTopologicalCleanerBase &)=delete
virtual void	update (const edm::EventSetup &)
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 8 of file SpikeAndDoubleSpikeCleaner.cc.

Constructor & Destructor Documentation

SpikeAndDoubleSpikeCleaner() [1/2]

SpikeAndDoubleSpikeCleaner::SpikeAndDoubleSpikeCleaner	(	const edm::ParameterSet &	conf,
		edm::ConsumesCollector &	cc
	)

Definition at line 98 of file SpikeAndDoubleSpikeCleaner.cc.

References PFLayer::ECAL_BARREL, PFLayer::ECAL_ENDCAP, PFLayer::HCAL_BARREL1, PFLayer::HCAL_BARREL2, PFLayer::HCAL_ENDCAP, PFLayer::HF_EM, PFLayer::HF_HAD, createfilelist::int, PFLayer::NONE, PFLayer::PS1, and PFLayer::PS2.

    : RecHitTopologicalCleanerBase(conf, cc),
      _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() [2/2]

SpikeAndDoubleSpikeCleaner::SpikeAndDoubleSpikeCleaner ( const SpikeAndDoubleSpikeCleaner &  )

delete

Member Function Documentation

clean()

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

overridevirtual

Implements RecHitTopologicalCleanerBase.

Definition at line 132 of file SpikeAndDoubleSpikeCleaner.cc.

References _thresholds, funct::abs(), AlCaHLTBitMon_QueryRunRegistry::comp, ALCARECOPPSCalTrackBasedSel_cff::detid, reco::PFRecHit::detId(), PFLayer::ECAL_BARREL, reco::PFRecHit::energy(), RhoEtaPhi::eta(), PVValHelper::eta, PFLayer::HCAL_BARREL2, HcalForward, PFLayer::HF_EM, PFLayer::HF_HAD, hfClusterShapes_cfi::hits, mps_fire::i, heavyIonCSV_trainingSettings::idx, recoMuon::in, input, createfilelist::int, dqmiolumiharvest::j, isotrackApplyRegressor::k, reco::PFRecHit::layer(), ALPAKA_ACCELERATOR_NAMESPACE::pixelClustering::pixelStatus::mask, SiStripPI::min, reco::PFRecHit::neighbours4(), phi, RhoEtaPhi::phi(), reco::PFRecHit::positionREP(), DetId::rawId(), and jetUpdater_cfi::sort.

                                                                                                            {
  //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
}

operator=()

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

delete

Member Data Documentation

_layerMap

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

private

Definition at line 27 of file SpikeAndDoubleSpikeCleaner.cc.

_thresholds

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

private

Definition at line 28 of file SpikeAndDoubleSpikeCleaner.cc.

Referenced by clean().