EgammaElectronProducers/plugins/ElectronSeedProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    ElectronProducers
0004 // Class:      ElectronSeedProducer
0005 //
0006 /**\class ElectronSeedProducer RecoEgamma/ElectronProducers/src/ElectronSeedProducer.cc
0007
0008  Description: EDProducer of ElectronSeed objects
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Ursula Berthon, Claude Charlot
0015 //         Created:  Mon Mar 27 13:22:06 CEST 2006
0016 //
0017 //
0018
0019 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0020 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronSeedGenerator.h"
0021 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronHcalHelper.h"
0022 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0026 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0027 #include "FWCore/Utilities/interface/transform.h"
0028 #include "FWCore/Framework/interface/stream/EDProducer.h"
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "RecoLocalCalo/HGCalRecAlgos/interface/ClusterTools.h"
0031 #include "RecoEcal/EgammaCoreTools/interface/EcalTools.h"
0032 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0033 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0034 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0035
0036 class ElectronSeedProducer : public edm::stream::EDProducer<> {
0037 public:
0038   explicit ElectronSeedProducer(const edm::ParameterSet&);
0039
0040   void produce(edm::Event&, const edm::EventSetup&) final;
0041
0042   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0043
0044 private:
0045   reco::SuperClusterRefVector filterClusters(math::XYZPoint const& beamSpotPosition,
0046                                              const edm::Handle<reco::SuperClusterCollection>& superClusters,
0047                                              HcalPFCuts const* hcalCuts) const;
0048
0049   edm::EDGetTokenT<reco::SuperClusterCollection> superClusters_[2];
0050   std::vector<edm::EDGetTokenT<TrajectorySeedCollection>> initialSeeds_;
0051   edm::EDGetTokenT<reco::BeamSpot> beamSpotTag_;
0052
0053   std::unique_ptr<ElectronSeedGenerator> matcher_;
0054
0055   bool applyHOverECut_ = true;
0056   std::unique_ptr<ElectronHcalHelper> hcalHelper_ = nullptr;
0057   double maxHOverEBarrel_;
0058   double maxHOverEEndcaps_;
0059   double SCEtCut_;
0060
0061   bool allowHGCal_;
0062   std::unique_ptr<hgcal::ClusterTools> hgcClusterTools_;
0063
0064   edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> hcalCutsToken_;
0065   bool cutsFromDB_;
0066 };
0067
0068 using namespace reco;
0069
0070 ElectronSeedProducer::ElectronSeedProducer(const edm::ParameterSet& conf)
0071     :
0072
0073       initialSeeds_{
0074           edm::vector_transform(conf.getParameter<std::vector<edm::InputTag>>("initialSeedsVector"),
0075                                 [this](edm::InputTag const& tag) { return consumes<TrajectorySeedCollection>(tag); })}
0076
0077 {
0078   SCEtCut_ = conf.getParameter<double>("SCEtCut");
0079   auto theconsumes = consumesCollector();
0080
0081   // new beamSpot tag
0082   beamSpotTag_ = consumes(conf.getParameter<edm::InputTag>("beamSpot"));
0083
0084   // for H/E
0085   applyHOverECut_ = conf.getParameter<bool>("applyHOverECut");
0086   if (applyHOverECut_) {
0087     ElectronHcalHelper::Configuration hcalCfg{};
0088     hcalCfg.hOverEConeSize = conf.getParameter<double>("hOverEConeSize");
0089     if (hcalCfg.hOverEConeSize > 0) {
0090       hcalCfg.onlyBehindCluster = false;
0091       hcalCfg.checkHcalStatus = false;
0092
0093       hcalCfg.hbheRecHits = consumes<HBHERecHitCollection>(conf.getParameter<edm::InputTag>("hbheRecHits"));
0094
0095       hcalCfg.eThresHB = conf.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
0096       hcalCfg.maxSeverityHB = conf.getParameter<int>("maxHcalRecHitSeverity");
0097       hcalCfg.eThresHE = conf.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
0098       hcalCfg.maxSeverityHE = hcalCfg.maxSeverityHB;
0099     }
0100     hcalHelper_ = std::make_unique<ElectronHcalHelper>(hcalCfg, consumesCollector());
0101
0102     allowHGCal_ = conf.getParameter<bool>("allowHGCal");
0103     if (allowHGCal_) {
0104       const edm::ParameterSet& hgcCfg = conf.getParameterSet("HGCalConfig");
0105       hgcClusterTools_ = std::make_unique<hgcal::ClusterTools>(hgcCfg, theconsumes);
0106     }
0107
0108     maxHOverEBarrel_ = conf.getParameter<double>("maxHOverEBarrel");
0109     maxHOverEEndcaps_ = conf.getParameter<double>("maxHOverEEndcaps");
0110   }
0111
0112   //Retrieve HCAL PF thresholds - from config or from DB
0113   cutsFromDB_ = conf.getParameter<bool>("usePFThresholdsFromDB");
0114   if (cutsFromDB_) {
0115     hcalCutsToken_ = esConsumes<HcalPFCuts, HcalPFCutsRcd>(edm::ESInputTag("", "withTopo"));
0116   }
0117
0118   ElectronSeedGenerator::Tokens esg_tokens;
0119   esg_tokens.token_bs = beamSpotTag_;
0120   esg_tokens.token_vtx = mayConsume<reco::VertexCollection>(conf.getParameter<edm::InputTag>("vertices"));
0121
0122   matcher_ = std::make_unique<ElectronSeedGenerator>(conf, esg_tokens, consumesCollector());
0123
0124   superClusters_[0] = consumes(conf.getParameter<edm::InputTag>("barrelSuperClusters"));
0125   superClusters_[1] = consumes(conf.getParameter<edm::InputTag>("endcapSuperClusters"));
0126
0127   //register your products
0128   produces<ElectronSeedCollection>();
0129 }
0130
0131 void ElectronSeedProducer::produce(edm::Event& e, const edm::EventSetup& iSetup) {
0132   LogDebug("ElectronSeedProducer") << "[ElectronSeedProducer::produce] entering ";
0133
0134   HcalPFCuts const* hcalCuts = nullptr;
0135   if (cutsFromDB_) {
0136     hcalCuts = &iSetup.getData(hcalCutsToken_);
0137   }
0138
0139   std::vector<TrajectorySeedCollection const*> initialSeedCollections;
0140   std::unique_ptr<TrajectorySeedCollection> initialSeedCollectionPtr = nullptr;  //created on the fly
0141
0142   if (hcalHelper_) {
0143     hcalHelper_->beginEvent(e, iSetup);
0144     if (allowHGCal_) {
0145       hgcClusterTools_->getEventSetup(iSetup);
0146       hgcClusterTools_->getEvent(e);
0147     }
0148   }
0149
0150   matcher_->setupES(iSetup);
0151
0152   // get initial TrajectorySeeds
0153   initialSeedCollections.clear();
0154   for (auto const& seeds : initialSeeds_) {
0155     initialSeedCollections.push_back(&e.get(seeds));
0156   }
0157
0158   auto seeds = std::make_unique<ElectronSeedCollection>();
0159   auto const& beamSpotPosition = e.get(beamSpotTag_).position();
0160
0161   // loop over barrel + endcap
0162   for (unsigned int i = 0; i < 2; i++) {
0163     auto clusterRefs = filterClusters(beamSpotPosition, e.getHandle(superClusters_[i]), hcalCuts);
0164     matcher_->run(e, clusterRefs, initialSeedCollections, *seeds);
0165   }
0166
0167   // store the accumulated result
0168 #ifdef EDM_ML_DEBUG
0169   for (auto const& seed : *seeds) {
0170     SuperClusterRef superCluster = seed.caloCluster().castTo<SuperClusterRef>();
0171     LogDebug("ElectronSeedProducer") << "new seed with " << seed.nHits() << " hits"
0172                                      << ", charge " << seed.getCharge() << " and cluster energy "
0173                                      << superCluster->energy() << " PID " << superCluster.id();
0174   }
0175 #endif
0176   e.put(std::move(seeds));
0177 }
0178
0179 //===============================
0180 // Filter the superclusters
0181 // - with EtCut
0182 // - with HoE using calo cone
0183 //===============================
0184
0185 SuperClusterRefVector ElectronSeedProducer::filterClusters(
0186     math::XYZPoint const& beamSpotPosition,
0187     const edm::Handle<reco::SuperClusterCollection>& superClusters,
0188     HcalPFCuts const* hcalCuts) const {
0189   SuperClusterRefVector sclRefs;
0190
0191   for (unsigned int i = 0; i < superClusters->size(); ++i) {
0192     auto const& scl = (*superClusters)[i];
0193     double sclEta = EleRelPoint(scl.position(), beamSpotPosition).eta();
0194     if (scl.energy() / cosh(sclEta) > SCEtCut_) {
0195       if (applyHOverECut_) {
0196         bool hoeVeto = false;
0197         double had = hcalHelper_->hcalESum(scl, 0, hcalCuts);
0198         double scle = scl.energy();
0199         int det_group = scl.seed()->hitsAndFractions()[0].first.det();
0200         int detector = scl.seed()->hitsAndFractions()[0].first.subdetId();
0201         if (detector == EcalBarrel && had / scle < maxHOverEBarrel_)
0202           hoeVeto = true;
0203         else if (!allowHGCal_ && detector == EcalEndcap && had / scle < maxHOverEEndcaps_)
0204           hoeVeto = true;
0205         else if (allowHGCal_ && EcalTools::isHGCalDet((DetId::Detector)det_group)) {
0206           float had_fraction = hgcClusterTools_->getClusterHadronFraction(*(scl.seed()));
0207           hoeVeto = (had_fraction >= 0.f && had_fraction < maxHOverEEndcaps_);
0208         }
0209         if (hoeVeto) {
0210           sclRefs.push_back({superClusters, i});
0211         }
0212       } else {
0213         sclRefs.push_back({superClusters, i});
0214       }
0215     }
0216   }
0217   LogDebug("ElectronSeedProducer") << "Filtered out " << sclRefs.size() << " superclusters from "
0218                                    << superClusters->size();
0219
0220   return sclRefs;
0221 }
0222
0223 void ElectronSeedProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0224   edm::ParameterSetDescription desc;
0225
0226   // steering
0227   desc.add<std::vector<edm::InputTag>>("initialSeedsVector", {});
0228   desc.add<bool>("useRecoVertex", false);
0229   desc.add<edm::InputTag>("vertices", {"offlinePrimaryVerticesWithBS"});
0230   desc.add<edm::InputTag>("beamSpot", {"offlineBeamSpot"});
0231   desc.add<bool>("dynamicPhiRoad", true);
0232
0233   // SC filtering
0234   desc.add<double>("SCEtCut", 0.0);
0235
0236   // H/E
0237   desc.add<bool>("applyHOverECut", true);
0238   desc.add<double>("hOverEConeSize", 0.15);
0239   desc.add<double>("maxHOverEBarrel", 0.15);
0240   desc.add<double>("maxHOverEEndcaps", 0.15);
0241   desc.add<edm::InputTag>("hbheRecHits", {"hbhereco"});
0242   desc.add<std::vector<double>>("recHitEThresholdHB", {0., 0., 0., 0.});
0243   desc.add<std::vector<double>>("recHitEThresholdHE", {0., 0., 0., 0., 0., 0., 0.});
0244   desc.add<int>("maxHcalRecHitSeverity", 999999);
0245   desc.add<bool>("usePFThresholdsFromDB", false);
0246
0247   // H/E equivalent for HGCal
0248   desc.add<bool>("allowHGCal", false);
0249   edm::ParameterSetDescription psd4;
0250   psd4.add<edm::InputTag>("HGCEEInput", {"HGCalRecHit", "HGCEERecHits"});
0251   psd4.add<edm::InputTag>("HGCFHInput", {"HGCalRecHit", "HGCHEFRecHits"});
0252   psd4.add<edm::InputTag>("HGCBHInput", {"HGCalRecHit", "HGCHEBRecHits"});
0253   psd4.add<edm::InputTag>("hgcalHitMap", {"recHitMapProducer", "hgcalRecHitMap"});
0254   desc.add<edm::ParameterSetDescription>("HGCalConfig", psd4);
0255
0256   // r/z windows
0257   desc.add<double>("nSigmasDeltaZ1", 5.0);      // in case beam spot is used for the matching
0258   desc.add<double>("deltaZ1WithVertex", 25.0);  // in case reco vertex is used for the matching
0259   desc.add<double>("z2MaxB", 0.09);
0260   desc.add<double>("r2MaxF", 0.15);
0261   desc.add<double>("rMaxI", 0.2);  // intermediate region SC in EB and 2nd hits in PXF
0262
0263   // phi windows (dynamic)
0264   desc.add<double>("LowPtThreshold", 5.0);
0265   desc.add<double>("HighPtThreshold", 35.0);
0266   desc.add<double>("SizeWindowENeg", 0.675);
0267   desc.add<double>("DeltaPhi1Low", 0.23);
0268   desc.add<double>("DeltaPhi1High", 0.08);
0269   desc.add<double>("DeltaPhi2B", 0.008);
0270   desc.add<double>("DeltaPhi2F", 0.012);
0271
0272   // phi windows (non dynamic, overwritten in case dynamic is selected)
0273   desc.add<double>("ePhiMin1", -0.125);
0274   desc.add<double>("ePhiMax1", 0.075);
0275   desc.add<double>("PhiMax2B", 0.002);
0276   desc.add<double>("PhiMax2F", 0.003);
0277
0278   desc.add<edm::InputTag>("barrelSuperClusters",
0279                           {"particleFlowSuperClusterECAL", "particleFlowSuperClusterECALBarrel"});
0280   desc.add<edm::InputTag>("endcapSuperClusters",
0281                           {"particleFlowSuperClusterECAL", "particleFlowSuperClusterECALEndcapWithPreshower"});
0282
0283   descriptions.add("ecalDrivenElectronSeedsDefault", desc);
0284 }
0285
0286 #include "FWCore/Framework/interface/MakerMacros.h"
0287 DEFINE_FWK_MODULE(ElectronSeedProducer);