PFBadHcalPseudoClusterProducer.cc

Go to the documentation of this file.

// system include files
#include <iostream>
#include <memory>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"

#include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
#include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"

#include "DataFormats/HcalDetId/interface/HcalDetId.h"

#include "CondFormats/HcalObjects/interface/HcalChannelStatus.h"
#include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
#include "CondFormats/HcalObjects/interface/HcalCondObjectContainer.h"
#include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloTopology/interface/HcalTopology.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"

#include "FWCore/Framework/interface/ESHandle.h"

using namespace std;
using namespace edm;

//
// class declaration
//

class PFBadHcalPseudoClusterProducer : public edm::stream::EDProducer<> {
public:
  explicit PFBadHcalPseudoClusterProducer(const edm::ParameterSet&);
  ~PFBadHcalPseudoClusterProducer() override;

  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
  virtual void init(const EventSetup& c);
  void produce(edm::Event&, const edm::EventSetup&) override;

  bool enabled_;
  bool debug_;

  edm::ESHandle<HcalChannelQuality> hQuality_;
  edm::ESHandle<CaloGeometry> hGeom_;
  unsigned long long cacheId_quality_, cacheId_geom_;
  std::vector<reco::PFCluster> badAreasC_;
  std::vector<reco::PFRecHit> badAreasRH_;
};

PFBadHcalPseudoClusterProducer::PFBadHcalPseudoClusterProducer(const edm::ParameterSet& ps)
    : enabled_(ps.getParameter<bool>("enable")),
      debug_(ps.getUntrackedParameter<bool>("debug", false)),
      cacheId_quality_(0),
      cacheId_geom_(0) {
  produces<std::vector<reco::PFCluster>>();
  produces<std::vector<reco::PFRecHit>>("hits");
}

PFBadHcalPseudoClusterProducer::~PFBadHcalPseudoClusterProducer() {}

void PFBadHcalPseudoClusterProducer::init(const EventSetup& iSetup) {
  badAreasC_.clear();
  badAreasRH_.clear();

  edm::ESHandle<HcalChannelQuality> hQuality_;
  iSetup.get<HcalChannelQualityRcd>().get("withTopo", hQuality_);
  const HcalChannelQuality& chanquality = *hQuality_;

  edm::ESHandle<CaloGeometry> hGeom_;
  iSetup.get<CaloGeometryRecord>().get(hGeom_);
  const CaloGeometry& caloGeom = *hGeom_;
  const CaloSubdetectorGeometry* hbGeom = caloGeom.getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
  const CaloSubdetectorGeometry* heGeom = caloGeom.getSubdetectorGeometry(DetId::Hcal, HcalEndcap);

  int statusMask = ((1 << HcalChannelStatus::HcalCellOff) | (1 << HcalChannelStatus::HcalCellMask) |
                    (1 << HcalChannelStatus::HcalCellDead));
  // histogram the number of bad depths at each ieta, iphi
  std::map<std::pair<int, int>, int> good, bads;
  std::map<std::pair<int, int>, std::pair<int, HcalSubdetector>> minDepths;
  for (const DetId& i : chanquality.getAllChannels()) {
    if (i.det() != DetId::Hcal)
      continue;  // not an hcal cell
    HcalDetId id = HcalDetId(i);
    if (id.subdet() != HcalBarrel && id.subdet() != HcalEndcap)
      continue;  // we don't deal with HO and HF here
    int status = chanquality.getValues(id)->getValue();
    auto tower = std::make_pair(id.ieta(), id.iphi());
    if (status & statusMask) {
      bads[tower]++;
      if (debug_)
        std::cout << "Channel " << i() << " (subdet " << id.subdet() << ", zside " << id.zside() << ", ieta "
                  << id.ieta() << ", iphi " << id.iphi() << " depth " << id.depth() << " has status " << status
                  << "  masked " << (status & statusMask) << std::endl;
    } else {
      good[tower]++;
    }
    auto& minD = minDepths[tower];
    if (minD.second == HcalEmpty || minD.first > id.depth()) {
      minD.first = id.depth();
      minD.second = id.subdet();
    }
  }

  const float dummyEnergy = 1e-5;  // non-zero, but small (even if it shouldn't ever be used)
  for (const auto& rec : bads) {
    int ieta = rec.first.first, iphi = rec.first.second, nbad = rec.second, ngood = good[rec.first];
    auto minDepth = minDepths[rec.first];
    bool barrel = minDepth.second == HcalBarrel;
    HcalDetId id(minDepth.second, ieta, iphi, minDepth.first);
    bool isBad = (nbad > 0 && nbad >= ngood);
    if (debug_)
      std::cout << "At ieta " << id.ieta() << ", iphi " << id.iphi() << " I have " << nbad << " bad depths, " << ngood
                << " good depths. First depth is in " << (barrel ? "HB" : "HE") << " depth " << minDepth.first << "; "
                << (isBad ? " MARK BAD" : " ignore") << std::endl;
    if (!isBad)
      continue;
    PFLayer::Layer layer = (barrel ? PFLayer::HCAL_BARREL1 : PFLayer::HCAL_ENDCAP);
    // make a PF RecHit
    std::shared_ptr<const CaloCellGeometry> thisCell = (barrel ? hbGeom : heGeom)->getGeometry(id);
    const GlobalPoint& pos = thisCell->getPosition();
    badAreasRH_.emplace_back(thisCell, id(), layer, dummyEnergy);
    // make a PF cluster (but can't add the rechit, as for that I need an edm::Ref)
    badAreasC_.emplace_back(layer, dummyEnergy, pos.x(), pos.y(), pos.z());
    badAreasC_.back().setFlags(reco::CaloCluster::badHcalMarker);
  }

  cacheId_quality_ = iSetup.get<HcalChannelQualityRcd>().cacheIdentifier();
  cacheId_geom_ = iSetup.get<CaloGeometryRecord>().cacheIdentifier();
}

// ------------ method called to produce the data  ------------
void PFBadHcalPseudoClusterProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  if (enabled_) {
    if (cacheId_quality_ != iSetup.get<HcalChannelQualityRcd>().cacheIdentifier() ||
        cacheId_geom_ != iSetup.get<CaloGeometryRecord>().cacheIdentifier()) {
      init(iSetup);
    }
  }

  auto outRH = std::make_unique<std::vector<reco::PFRecHit>>(badAreasRH_);
  auto rhHandle = iEvent.put(std::move(outRH), "hits");

  auto outC = std::make_unique<std::vector<reco::PFCluster>>(badAreasC_);

  // now we go set references
  for (unsigned int i = 0, n = rhHandle->size(); i < n; ++i) {
    (*outC)[i].addRecHitFraction(reco::PFRecHitFraction(reco::PFRecHitRef(rhHandle, i), 1.0));
  }

  iEvent.put(std::move(outC));
}

void PFBadHcalPseudoClusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.add<bool>("enable", false)
      ->setComment("activate the module (if false, it doesn't check the DB and produces an empty collection)");
  desc.addUntracked<bool>("debug", false);
  descriptions.add("particleFlowBadHcalPseudoCluster", desc);
}

//define this as a plug-in
DEFINE_FWK_MODULE(PFBadHcalPseudoClusterProducer);