File indexing completed on 2023-10-25 09:59:13
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "DataFormats/Candidate/interface/Candidate.h"
0010 #include "DataFormats/Candidate/interface/CandAssociation.h"
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "DataFormats/DetId/interface/DetId.h"
0013 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0014 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0015 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0016 #include "FWCore/Framework/interface/ConsumesCollector.h"
0017 #include "FWCore/Framework/interface/ESHandle.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/global/EDProducer.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/Utilities/interface/ESGetToken.h"
0023 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0025 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRecHitIsolation.h"
0026 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0027 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0028
0029 class EgammaEcalRecHitIsolationProducer : public edm::global::EDProducer<> {
0030 public:
0031 explicit EgammaEcalRecHitIsolationProducer(const edm::ParameterSet&);
0032
0033 void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0034
0035 private:
0036 const edm::EDGetTokenT<edm::View<reco::Candidate>> emObjectProducer_;
0037 const edm::EDGetTokenT<EcalRecHitCollection> ecalBarrelRecHitCollection_;
0038 const edm::EDGetTokenT<EcalRecHitCollection> ecalEndcapRecHitCollection_;
0039
0040 double egIsoPtMinBarrel_;
0041 double egIsoEMinBarrel_;
0042 double egIsoPtMinEndcap_;
0043 double egIsoEMinEndcap_;
0044 double egIsoConeSizeOut_;
0045 double egIsoConeSizeInBarrel_;
0046 double egIsoConeSizeInEndcap_;
0047 double egIsoJurassicWidth_;
0048
0049 bool useIsolEt_;
0050 bool tryBoth_;
0051 bool subtract_;
0052
0053 bool useNumCrystals_;
0054 bool vetoClustered_;
0055
0056 edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> sevLvToken_;
0057 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometrytoken_;
0058 };
0059
0060 #include "FWCore/Framework/interface/MakerMacros.h"
0061 DEFINE_FWK_MODULE(EgammaEcalRecHitIsolationProducer);
0062
0063 EgammaEcalRecHitIsolationProducer::EgammaEcalRecHitIsolationProducer(const edm::ParameterSet& config)
0064
0065 : emObjectProducer_{consumes(config.getParameter<edm::InputTag>("emObjectProducer"))},
0066 ecalBarrelRecHitCollection_{consumes(config.getParameter<edm::InputTag>("ecalBarrelRecHitCollection"))},
0067 ecalEndcapRecHitCollection_{consumes(config.getParameter<edm::InputTag>("ecalEndcapRecHitCollection"))} {
0068
0069 egIsoPtMinBarrel_ = config.getParameter<double>("etMinBarrel");
0070 egIsoEMinBarrel_ = config.getParameter<double>("eMinBarrel");
0071 egIsoPtMinEndcap_ = config.getParameter<double>("etMinEndcap");
0072 egIsoEMinEndcap_ = config.getParameter<double>("eMinEndcap");
0073 egIsoConeSizeInBarrel_ = config.getParameter<double>("intRadiusBarrel");
0074 egIsoConeSizeInEndcap_ = config.getParameter<double>("intRadiusEndcap");
0075 egIsoConeSizeOut_ = config.getParameter<double>("extRadius");
0076 egIsoJurassicWidth_ = config.getParameter<double>("jurassicWidth");
0077
0078
0079 useIsolEt_ = config.getParameter<bool>("useIsolEt");
0080 tryBoth_ = config.getParameter<bool>("tryBoth");
0081 subtract_ = config.getParameter<bool>("subtract");
0082 useNumCrystals_ = config.getParameter<bool>("useNumCrystals");
0083 vetoClustered_ = config.getParameter<bool>("vetoClustered");
0084
0085
0086 sevLvToken_ = esConsumes();
0087 caloGeometrytoken_ = esConsumes();
0088
0089
0090 produces<edm::ValueMap<double>>();
0091 }
0092
0093
0094 void EgammaEcalRecHitIsolationProducer::produce(edm::StreamID,
0095 edm::Event& iEvent,
0096 const edm::EventSetup& iSetup) const {
0097
0098 auto emObjectHandle = iEvent.getHandle(emObjectProducer_);
0099
0100
0101 auto ecalBarrelRecHitHandle = iEvent.getHandle(ecalBarrelRecHitCollection_);
0102
0103
0104 auto ecalEndcapRecHitHandle = iEvent.getHandle(ecalEndcapRecHitCollection_);
0105
0106 edm::ESHandle<EcalSeverityLevelAlgo> sevlv = iSetup.getHandle(sevLvToken_);
0107 const EcalSeverityLevelAlgo* sevLevel = sevlv.product();
0108
0109
0110 edm::ESHandle<CaloGeometry> pG = iSetup.getHandle(caloGeometrytoken_);
0111 const CaloGeometry* caloGeom = pG.product();
0112
0113
0114 auto isoMap = std::make_unique<edm::ValueMap<double>>();
0115 edm::ValueMap<double>::Filler filler(*isoMap);
0116 std::vector<double> retV(emObjectHandle->size(), 0);
0117
0118 EgammaRecHitIsolation ecalBarrelIsol(egIsoConeSizeOut_,
0119 egIsoConeSizeInBarrel_,
0120 egIsoJurassicWidth_,
0121 egIsoPtMinBarrel_,
0122 egIsoEMinBarrel_,
0123 caloGeom,
0124 *ecalBarrelRecHitHandle,
0125 sevLevel,
0126 DetId::Ecal);
0127 ecalBarrelIsol.setUseNumCrystals(useNumCrystals_);
0128 ecalBarrelIsol.setVetoClustered(vetoClustered_);
0129
0130 EgammaRecHitIsolation ecalEndcapIsol(egIsoConeSizeOut_,
0131 egIsoConeSizeInEndcap_,
0132 egIsoJurassicWidth_,
0133 egIsoPtMinEndcap_,
0134 egIsoEMinEndcap_,
0135 caloGeom,
0136 *ecalEndcapRecHitHandle,
0137 sevLevel,
0138 DetId::Ecal);
0139 ecalEndcapIsol.setUseNumCrystals(useNumCrystals_);
0140 ecalEndcapIsol.setVetoClustered(vetoClustered_);
0141
0142 for (size_t i = 0; i < emObjectHandle->size(); ++i) {
0143
0144
0145
0146
0147 double isoValue = 0.;
0148
0149 reco::SuperClusterRef superClus = emObjectHandle->at(i).get<reco::SuperClusterRef>();
0150
0151 if (tryBoth_) {
0152 if (useIsolEt_)
0153 isoValue =
0154 ecalBarrelIsol.getEtSum(&(emObjectHandle->at(i))) + ecalEndcapIsol.getEtSum(&(emObjectHandle->at(i)));
0155 else
0156 isoValue = ecalBarrelIsol.getEnergySum(&(emObjectHandle->at(i))) +
0157 ecalEndcapIsol.getEnergySum(&(emObjectHandle->at(i)));
0158 } else if (fabs(superClus->eta()) < 1.479) {
0159 if (useIsolEt_)
0160 isoValue = ecalBarrelIsol.getEtSum(&(emObjectHandle->at(i)));
0161 else
0162 isoValue = ecalBarrelIsol.getEnergySum(&(emObjectHandle->at(i)));
0163 } else {
0164 if (useIsolEt_)
0165 isoValue = ecalEndcapIsol.getEtSum(&(emObjectHandle->at(i)));
0166 else
0167 isoValue = ecalEndcapIsol.getEnergySum(&(emObjectHandle->at(i)));
0168 }
0169
0170
0171 double subtractVal = 0;
0172
0173 if (useIsolEt_)
0174 subtractVal = superClus.get()->rawEnergy() * sin(2 * atan(exp(-superClus.get()->eta())));
0175 else
0176 subtractVal = superClus.get()->rawEnergy();
0177
0178 if (subtract_)
0179 isoValue -= subtractVal;
0180
0181 retV[i] = isoValue;
0182
0183
0184 }
0185
0186 filler.insert(emObjectHandle, retV.begin(), retV.end());
0187 filler.fill();
0188
0189 iEvent.put(std::move(isoMap));
0190 }
0191
0192
0193