plugins/importers/SuperClusterImporter.cc

0001 #include "RecoParticleFlow/PFProducer/interface/BlockElementImporterBase.h"
0002 #include "RecoParticleFlow/PFProducer/interface/PhotonSelectorAlgo.h"
0003 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0004 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementSuperCluster.h"
0005 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0006 #include "RecoParticleFlow/PFProducer/interface/PFBlockElementSCEqual.h"
0007 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0008 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHcalIsolation.h"
0009 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0010 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0011 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0012 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0013 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0014 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
0015 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
0016 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
0017 #include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
0018
0019 //quick pT for superclusters
0020 inline double ptFast(const double energy, const math::XYZPoint& position, const math::XYZPoint& origin) {
0021   const auto v = position - origin;
0022   return energy * std::sqrt(v.perp2() / v.mag2());
0023 }
0024
0025 #include <memory>
0026
0027 #include <unordered_map>
0028
0029 class SuperClusterImporter : public BlockElementImporterBase {
0030 public:
0031   SuperClusterImporter(const edm::ParameterSet&, edm::ConsumesCollector&);
0032
0033   void updateEventSetup(const edm::EventSetup& es) override;
0034
0035   void importToBlock(const edm::Event&, ElementList&) const override;
0036
0037 private:
0038   edm::EDGetTokenT<reco::SuperClusterCollection> _srcEB, _srcEE;
0039   const double _maxHoverE, _pTbyPass, _minSCPt;
0040   const edm::EDGetTokenT<HBHERecHitCollection> hbheRecHitsTag_;
0041   const int maxSeverityHB_;
0042   const int maxSeverityHE_;
0043   bool cutsFromDB;
0044   CaloTowerConstituentsMap const* towerMap_;
0045   CaloGeometry const* caloGeom_;
0046   HcalTopology const* hcalTopo_;
0047   HcalChannelQuality const* hcalChannelQual_;
0048   HcalSeverityLevelComputer const* hcalSev_;
0049   bool _superClustersArePF;
0050   static const math::XYZPoint _zero;
0051
0052   const edm::ESGetToken<CaloTowerConstituentsMap, CaloGeometryRecord> _ctmapToken;
0053   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0054   const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> hcalTopologyToken_;
0055   const edm::ESGetToken<HcalChannelQuality, HcalChannelQualityRcd> hcalChannelQualityToken_;
0056   const edm::ESGetToken<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd> hcalSevLvlComputerToken_;
0057   edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> hcalCutsToken_;
0058   HcalPFCuts const* hcalCuts = nullptr;
0059 };
0060
0061 const math::XYZPoint SuperClusterImporter::_zero = math::XYZPoint(0, 0, 0);
0062
0063 DEFINE_EDM_PLUGIN(BlockElementImporterFactory, SuperClusterImporter, "SuperClusterImporter");
0064
0065 SuperClusterImporter::SuperClusterImporter(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0066     : BlockElementImporterBase(conf, cc),
0067       _srcEB(cc.consumes<reco::SuperClusterCollection>(conf.getParameter<edm::InputTag>("source_eb"))),
0068       _srcEE(cc.consumes<reco::SuperClusterCollection>(conf.getParameter<edm::InputTag>("source_ee"))),
0069       _maxHoverE(conf.getParameter<double>("maximumHoverE")),
0070       _pTbyPass(conf.getParameter<double>("minPTforBypass")),
0071       _minSCPt(conf.getParameter<double>("minSuperClusterPt")),
0072       hbheRecHitsTag_(cc.consumes(conf.getParameter<edm::InputTag>("hbheRecHitsTag"))),
0073       maxSeverityHB_(conf.getParameter<int>("maxSeverityHB")),
0074       maxSeverityHE_(conf.getParameter<int>("maxSeverityHE")),
0075       cutsFromDB(conf.getParameter<bool>("usePFThresholdsFromDB")),
0076       _superClustersArePF(conf.getParameter<bool>("superClustersArePF")),
0077       _ctmapToken(cc.esConsumes<edm::Transition::BeginLuminosityBlock>()),
0078       caloGeometryToken_{cc.esConsumes<edm::Transition::BeginLuminosityBlock>()},
0079       hcalTopologyToken_{cc.esConsumes<edm::Transition::BeginLuminosityBlock>()},
0080       hcalChannelQualityToken_{cc.esConsumes<edm::Transition::BeginLuminosityBlock>(edm::ESInputTag("", "withTopo"))},
0081       hcalSevLvlComputerToken_{cc.esConsumes<edm::Transition::BeginLuminosityBlock>()} {
0082   if (cutsFromDB) {
0083     hcalCutsToken_ = cc.esConsumes<HcalPFCuts, HcalPFCutsRcd, edm::Transition::BeginLuminosityBlock>(
0084         edm::ESInputTag("", "withTopo"));
0085   }
0086 }
0087
0088 void SuperClusterImporter::updateEventSetup(const edm::EventSetup& es) {
0089   towerMap_ = &es.getData(_ctmapToken);
0090   if (cutsFromDB) {
0091     hcalCuts = &es.getData(hcalCutsToken_);
0092   }
0093   caloGeom_ = &es.getData(caloGeometryToken_);
0094   hcalTopo_ = &es.getData(hcalTopologyToken_);
0095   hcalChannelQual_ = &es.getData(hcalChannelQualityToken_);
0096   hcalSev_ = &es.getData(hcalSevLvlComputerToken_);
0097 }
0098
0099 void SuperClusterImporter::importToBlock(const edm::Event& e, BlockElementImporterBase::ElementList& elems) const {
0100   auto eb_scs = e.getHandle(_srcEB);
0101   auto ee_scs = e.getHandle(_srcEE);
0102   elems.reserve(elems.size() + eb_scs->size() + ee_scs->size());
0103   // setup our elements so that all the SCs are grouped together
0104   auto SCs_end =
0105       std::partition(elems.begin(), elems.end(), [](auto const& a) { return a->type() == reco::PFBlockElement::SC; });
0106   // add eb superclusters
0107   auto bsc = eb_scs->cbegin();
0108   auto esc = eb_scs->cend();
0109   reco::PFBlockElementSuperCluster* scbe = nullptr;
0110   reco::SuperClusterRef scref;
0111
0112   EgammaHcalIsolation thisHcalVar_ = EgammaHcalIsolation(EgammaHcalIsolation::InclusionRule::isBehindClusterSeed,
0113                                                          0,  //outercone
0114                                                          EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
0115                                                          0,  //innercone
0116                                                          {{0, 0, 0, 0}},
0117                                                          {{0, 0, 0, 0}},
0118                                                          maxSeverityHB_,
0119                                                          {{0, 0, 0, 0, 0, 0, 0}},
0120                                                          {{0, 0, 0, 0, 0, 0, 0}},
0121                                                          maxSeverityHE_,
0122                                                          e.get(hbheRecHitsTag_),
0123                                                          caloGeom_,
0124                                                          hcalTopo_,
0125                                                          hcalChannelQual_,
0126                                                          hcalSev_,
0127                                                          towerMap_);
0128
0129   for (auto sc = bsc; sc != esc; ++sc) {
0130     scref = reco::SuperClusterRef(eb_scs, std::distance(bsc, sc));
0131     PFBlockElementSCEqual myEqual(scref);
0132     auto sc_elem = std::find_if(elems.begin(), SCs_end, myEqual);
0133     const double scpT = ptFast(sc->energy(), sc->position(), _zero);
0134     const double H_tower = thisHcalVar_.getHcalESumBc(scref.get(), 0, hcalCuts);
0135     const double HoverE = H_tower / sc->energy();
0136     if (sc_elem == SCs_end && scpT > _minSCPt && (scpT > _pTbyPass || HoverE < _maxHoverE)) {
0137       scbe = new reco::PFBlockElementSuperCluster(scref);
0138       scbe->setFromPFSuperCluster(_superClustersArePF);
0139       SCs_end = elems.emplace(SCs_end, scbe);
0140       ++SCs_end;  // point to element *after* the new one
0141     }
0142   }  // loop on eb superclusters
0143   // add ee superclusters
0144   bsc = ee_scs->cbegin();
0145   esc = ee_scs->cend();
0146   for (auto sc = bsc; sc != esc; ++sc) {
0147     scref = reco::SuperClusterRef(ee_scs, std::distance(bsc, sc));
0148     PFBlockElementSCEqual myEqual(scref);
0149     auto sc_elem = std::find_if(elems.begin(), SCs_end, myEqual);
0150     const double scpT = ptFast(sc->energy(), sc->position(), _zero);
0151     const double H_tower = thisHcalVar_.getHcalESumBc(scref.get(), 0, hcalCuts);
0152     const double HoverE = H_tower / sc->energy();
0153     if (sc_elem == SCs_end && scpT > _minSCPt && (scpT > _pTbyPass || HoverE < _maxHoverE)) {
0154       scbe = new reco::PFBlockElementSuperCluster(scref);
0155       scbe->setFromPFSuperCluster(_superClustersArePF);
0156       SCs_end = elems.emplace(SCs_end, scbe);
0157       ++SCs_end;  // point to element *after* the new one
0158     }
0159   }  // loop on ee superclusters
0160   elems.shrink_to_fit();
0161 }