PFClusterProducer/plugins/PFBadHcalPseudoClusterProducer.cc

0001 // system include files
0002 #include <iostream>
0003 #include <memory>
0004
0005 // user include files
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/stream/EDProducer.h"
0008
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0016
0017 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0018 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
0019 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0020
0021 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0022
0023 #include "CondFormats/HcalObjects/interface/HcalChannelStatus.h"
0024 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
0025 #include "CondFormats/HcalObjects/interface/HcalCondObjectContainer.h"
0026 #include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
0027 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0028 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0029 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0030 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0031
0032 #include "FWCore/Framework/interface/ESHandle.h"
0033
0034 using namespace std;
0035 using namespace edm;
0036
0037 //
0038 // class declaration
0039 //
0040
0041 class PFBadHcalPseudoClusterProducer : public edm::stream::EDProducer<> {
0042 public:
0043   explicit PFBadHcalPseudoClusterProducer(const edm::ParameterSet&);
0044   ~PFBadHcalPseudoClusterProducer() override;
0045
0046   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0047
0048 private:
0049   virtual void init(const EventSetup& c);
0050   void produce(edm::Event&, const edm::EventSetup&) override;
0051
0052   bool enabled_;
0053   bool debug_;
0054
0055   edm::ESHandle<HcalChannelQuality> hQuality_;
0056   edm::ESHandle<CaloGeometry> hGeom_;
0057   unsigned long long cacheId_quality_, cacheId_geom_;
0058   std::vector<reco::PFCluster> badAreasC_;
0059   std::vector<reco::PFRecHit> badAreasRH_;
0060
0061   edm::ESGetToken<HcalChannelQuality, HcalChannelQualityRcd> hQualityToken_;
0062   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> hGeomToken_;
0063 };
0064
0065 PFBadHcalPseudoClusterProducer::PFBadHcalPseudoClusterProducer(const edm::ParameterSet& ps)
0066     : enabled_(ps.getParameter<bool>("enable")),
0067       debug_(ps.getUntrackedParameter<bool>("debug", false)),
0068       cacheId_quality_(0),
0069       cacheId_geom_(0),
0070       hQualityToken_(esConsumes(edm::ESInputTag("", "withTopo"))),
0071       hGeomToken_(esConsumes()) {
0072   produces<std::vector<reco::PFCluster>>();
0073   produces<std::vector<reco::PFRecHit>>("hits");
0074 }
0075
0076 PFBadHcalPseudoClusterProducer::~PFBadHcalPseudoClusterProducer() {}
0077
0078 void PFBadHcalPseudoClusterProducer::init(const EventSetup& iSetup) {
0079   badAreasC_.clear();
0080   badAreasRH_.clear();
0081
0082   hQuality_ = iSetup.getHandle(hQualityToken_);
0083   const HcalChannelQuality& chanquality = *hQuality_;
0084
0085   hGeom_ = iSetup.getHandle(hGeomToken_);
0086   const CaloGeometry& caloGeom = *hGeom_;
0087   const CaloSubdetectorGeometry* hbGeom = caloGeom.getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0088   const CaloSubdetectorGeometry* heGeom = caloGeom.getSubdetectorGeometry(DetId::Hcal, HcalEndcap);
0089
0090   int statusMask = ((1 << HcalChannelStatus::HcalCellOff) | (1 << HcalChannelStatus::HcalCellMask) |
0091                     (1 << HcalChannelStatus::HcalCellDead));
0092   // histogram the number of bad depths at each ieta, iphi
0093   std::map<std::pair<int, int>, int> good, bads;
0094   std::map<std::pair<int, int>, std::pair<int, HcalSubdetector>> minDepths;
0095   for (const DetId& i : chanquality.getAllChannels()) {
0096     if (i.det() != DetId::Hcal)
0097       continue;  // not an hcal cell
0098     HcalDetId id = HcalDetId(i);
0099     if (id.subdet() != HcalBarrel && id.subdet() != HcalEndcap)
0100       continue;  // we don't deal with HO and HF here
0101     int status = chanquality.getValues(id)->getValue();
0102     auto tower = std::make_pair(id.ieta(), id.iphi());
0103     if (status & statusMask) {
0104       bads[tower]++;
0105       if (debug_)
0106         std::cout << "Channel " << i() << " (subdet " << id.subdet() << ", zside " << id.zside() << ", ieta "
0107                   << id.ieta() << ", iphi " << id.iphi() << " depth " << id.depth() << " has status " << status
0108                   << "  masked " << (status & statusMask) << std::endl;
0109     } else {
0110       good[tower]++;
0111     }
0112     auto& minD = minDepths[tower];
0113     if (minD.second == HcalEmpty || minD.first > id.depth()) {
0114       minD.first = id.depth();
0115       minD.second = id.subdet();
0116     }
0117   }
0118
0119   const float dummyEnergy = 1e-5;  // non-zero, but small (even if it shouldn't ever be used)
0120   for (const auto& rec : bads) {
0121     int ieta = rec.first.first, iphi = rec.first.second, nbad = rec.second, ngood = good[rec.first];
0122     auto minDepth = minDepths[rec.first];
0123     bool barrel = minDepth.second == HcalBarrel;
0124     HcalDetId id(minDepth.second, ieta, iphi, minDepth.first);
0125     bool isBad = (nbad > 0 && nbad >= ngood);
0126     if (debug_)
0127       std::cout << "At ieta " << id.ieta() << ", iphi " << id.iphi() << " I have " << nbad << " bad depths, " << ngood
0128                 << " good depths. First depth is in " << (barrel ? "HB" : "HE") << " depth " << minDepth.first << "; "
0129                 << (isBad ? " MARK BAD" : " ignore") << std::endl;
0130     if (!isBad)
0131       continue;
0132     PFLayer::Layer layer = (barrel ? PFLayer::HCAL_BARREL1 : PFLayer::HCAL_ENDCAP);
0133     // make a PF RecHit
0134     std::shared_ptr<const CaloCellGeometry> thisCell = (barrel ? hbGeom : heGeom)->getGeometry(id);
0135     const GlobalPoint& pos = thisCell->getPosition();
0136     badAreasRH_.emplace_back(thisCell, id(), layer, dummyEnergy);
0137     // make a PF cluster (but can't add the rechit, as for that I need an edm::Ref)
0138     badAreasC_.emplace_back(layer, dummyEnergy, pos.x(), pos.y(), pos.z());
0139     badAreasC_.back().setFlags(reco::CaloCluster::badHcalMarker);
0140   }
0141
0142   cacheId_quality_ = iSetup.get<HcalChannelQualityRcd>().cacheIdentifier();
0143   cacheId_geom_ = iSetup.get<CaloGeometryRecord>().cacheIdentifier();
0144 }
0145
0146 // ------------ method called to produce the data  ------------
0147 void PFBadHcalPseudoClusterProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0148   if (enabled_) {
0149     if (cacheId_quality_ != iSetup.get<HcalChannelQualityRcd>().cacheIdentifier() ||
0150         cacheId_geom_ != iSetup.get<CaloGeometryRecord>().cacheIdentifier()) {
0151       init(iSetup);
0152     }
0153   }
0154
0155   auto outRH = std::make_unique<std::vector<reco::PFRecHit>>(badAreasRH_);
0156   auto rhHandle = iEvent.put(std::move(outRH), "hits");
0157
0158   auto outC = std::make_unique<std::vector<reco::PFCluster>>(badAreasC_);
0159
0160   // now we go set references
0161   for (unsigned int i = 0, n = rhHandle->size(); i < n; ++i) {
0162     (*outC)[i].addRecHitFraction(reco::PFRecHitFraction(reco::PFRecHitRef(rhHandle, i), 1.0));
0163   }
0164
0165   iEvent.put(std::move(outC));
0166 }
0167
0168 void PFBadHcalPseudoClusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0169   edm::ParameterSetDescription desc;
0170   desc.add<bool>("enable", true)
0171       ->setComment("activate the module (if false, it doesn't check the DB and produces an empty collection)");
0172   desc.addUntracked<bool>("debug", false);
0173   descriptions.add("particleFlowBadHcalPseudoCluster", desc);
0174 }
0175
0176 //define this as a plug-in
0177 DEFINE_FWK_MODULE(PFBadHcalPseudoClusterProducer);