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