File indexing completed on 2023-03-17 11:17:38
0001
0002
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/stream/EDProducer.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008
0009 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
0010 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0011 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0013
0014 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0015 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
0016
0017 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0018 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0019 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0020 #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
0021
0022 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0023 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateFwd.h"
0024
0025 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0026 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0027 #include "FWCore/Framework/interface/ESHandle.h"
0028
0029 #include "EgammaHLTPixelMatchParamObjects.h"
0030
0031 namespace {
0032
0033
0034
0035 int makeSeedInfo(const reco::ElectronSeed& seed) {
0036 int info = 0;
0037 for (size_t hitNr = 0; hitNr < seed.hitInfo().size(); hitNr++) {
0038 int subDetBit = 0x1 << hitNr;
0039 if (seed.subDet(hitNr) == PixelSubdetector::PixelEndcap)
0040 info |= subDetBit;
0041 int layerOffset = 3;
0042 if (seed.subDet(hitNr) == PixelSubdetector::PixelEndcap)
0043 layerOffset += 4;
0044 int layerBit = 0x1 << layerOffset << seed.layerOrDiskNr(hitNr);
0045 info |= layerBit;
0046
0047 int nrLayersAlongTrajShifted = seed.nrLayersAlongTraj() << 12;
0048 info |= nrLayersAlongTrajShifted;
0049 }
0050 return info;
0051 }
0052
0053 }
0054
0055 struct PixelData {
0056 public:
0057 PixelData(std::string name,
0058 size_t hitNr,
0059 float (reco::ElectronSeed::*func)(size_t) const,
0060 const edm::Handle<reco::RecoEcalCandidateCollection>& candHandle)
0061 : name_(std::move(name)), hitNr_(hitNr), func_(func), val_(std::numeric_limits<float>::max()), valInfo_(0) {
0062 valMap_ = std::make_unique<reco::RecoEcalCandidateIsolationMap>(candHandle);
0063 valInfoMap_ = std::make_unique<reco::RecoEcalCandidateIsolationMap>(candHandle);
0064 }
0065 PixelData(PixelData&& rhs) = default;
0066
0067 void resetVal() {
0068 val_ = std::numeric_limits<float>::max();
0069 valInfo_ = 0;
0070 }
0071 void fill(const reco::ElectronSeed& seed) {
0072 if (hitNr_ < seed.hitInfo().size()) {
0073 float seedVal = (seed.*func_)(hitNr_);
0074 if (std::abs(seedVal) < std::abs(val_)) {
0075 val_ = seedVal;
0076 valInfo_ = makeSeedInfo(seed);
0077 }
0078 }
0079 }
0080 void fill(const reco::RecoEcalCandidateRef& candRef) {
0081 valMap_->insert(candRef, val_);
0082 valInfoMap_->insert(candRef, valInfo_);
0083 val_ = std::numeric_limits<float>::max();
0084 valInfo_ = 0;
0085 }
0086
0087 void putInto(edm::Event& event) {
0088 event.put(std::move(valMap_), name_ + std::to_string(hitNr_ + 1));
0089 event.put(std::move(valInfoMap_), name_ + std::to_string(hitNr_ + 1) + "Info");
0090 }
0091
0092 private:
0093 std::unique_ptr<reco::RecoEcalCandidateIsolationMap> valMap_;
0094 std::unique_ptr<reco::RecoEcalCandidateIsolationMap> valInfoMap_;
0095 std::string name_;
0096 size_t hitNr_;
0097 float (reco::ElectronSeed::*func_)(size_t) const;
0098 float val_;
0099 float valInfo_;
0100 };
0101
0102 class EgammaHLTPixelMatchVarProducer : public edm::stream::EDProducer<> {
0103 public:
0104 explicit EgammaHLTPixelMatchVarProducer(const edm::ParameterSet&);
0105 ~EgammaHLTPixelMatchVarProducer() override;
0106
0107 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0108 void produce(edm::Event&, const edm::EventSetup&) override;
0109 std::array<float, 4> calS2(const reco::ElectronSeed& seed, int charge) const;
0110
0111 private:
0112
0113
0114 const edm::EDGetTokenT<reco::RecoEcalCandidateCollection> recoEcalCandidateToken_;
0115 const edm::EDGetTokenT<reco::ElectronSeedCollection> pixelSeedsToken_;
0116
0117 egPM::Param<reco::ElectronSeed> dPhi1Para_;
0118 egPM::Param<reco::ElectronSeed> dPhi2Para_;
0119 egPM::Param<reco::ElectronSeed> dRZ2Para_;
0120
0121 int productsToWrite_;
0122 size_t nrHits_;
0123 };
0124
0125 EgammaHLTPixelMatchVarProducer::EgammaHLTPixelMatchVarProducer(const edm::ParameterSet& config)
0126 : recoEcalCandidateToken_(
0127 consumes<reco::RecoEcalCandidateCollection>(config.getParameter<edm::InputTag>("recoEcalCandidateProducer"))),
0128 pixelSeedsToken_(
0129 consumes<reco::ElectronSeedCollection>(config.getParameter<edm::InputTag>("pixelSeedsProducer"))),
0130 dPhi1Para_(config.getParameter<edm::ParameterSet>("dPhi1SParams")),
0131 dPhi2Para_(config.getParameter<edm::ParameterSet>("dPhi2SParams")),
0132 dRZ2Para_(config.getParameter<edm::ParameterSet>("dRZ2SParams")),
0133 productsToWrite_(config.getParameter<int>("productsToWrite")),
0134 nrHits_(4)
0135
0136 {
0137
0138 produces<reco::RecoEcalCandidateIsolationMap>("s2");
0139 if (productsToWrite_ >= 1) {
0140 produces<reco::RecoEcalCandidateIsolationMap>("dPhi1BestS2");
0141 produces<reco::RecoEcalCandidateIsolationMap>("dPhi2BestS2");
0142 produces<reco::RecoEcalCandidateIsolationMap>("dzBestS2");
0143 }
0144 if (productsToWrite_ >= 2) {
0145
0146 for (size_t hitNr = 1; hitNr <= nrHits_; hitNr++) {
0147 produces<reco::RecoEcalCandidateIsolationMap>("dPhi" + std::to_string(hitNr));
0148 produces<reco::RecoEcalCandidateIsolationMap>("dPhi" + std::to_string(hitNr) + "Info");
0149 produces<reco::RecoEcalCandidateIsolationMap>("dRZ" + std::to_string(hitNr));
0150 produces<reco::RecoEcalCandidateIsolationMap>("dRZ" + std::to_string(hitNr) + "Info");
0151 }
0152 produces<reco::RecoEcalCandidateIsolationMap>("nrClus");
0153 produces<reco::RecoEcalCandidateIsolationMap>("seedClusEFrac");
0154 produces<reco::RecoEcalCandidateIsolationMap>("phiWidth");
0155 produces<reco::RecoEcalCandidateIsolationMap>("etaWidth");
0156 }
0157 }
0158
0159 EgammaHLTPixelMatchVarProducer::~EgammaHLTPixelMatchVarProducer() {}
0160
0161 void EgammaHLTPixelMatchVarProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0162 edm::ParameterSetDescription desc;
0163 desc.add<edm::InputTag>(("recoEcalCandidateProducer"), edm::InputTag("hltL1SeededRecoEcalCandidate"));
0164 desc.add<edm::InputTag>(("pixelSeedsProducer"), edm::InputTag("electronPixelSeeds"));
0165
0166 edm::ParameterSetDescription varParamDesc;
0167 edm::ParameterSetDescription binParamDesc;
0168
0169 auto binDescCases = "AbsEtaClus" >> (edm::ParameterDescription<double>("xMin", 0.0, true) and
0170 edm::ParameterDescription<double>("xMax", 3.0, true) and
0171 edm::ParameterDescription<int>("yMin", 0, true) and
0172 edm::ParameterDescription<int>("yMax", 99999, true) and
0173 edm::ParameterDescription<std::string>("funcType", "pol0", true) and
0174 edm::ParameterDescription<std::vector<double>>("funcParams", {0.}, true)) or
0175 "AbsEtaClusPhi" >> (edm::ParameterDescription<double>("xMin", 0.0, true) and
0176 edm::ParameterDescription<double>("xMax", 3.0, true) and
0177 edm::ParameterDescription<int>("yMin", 0, true) and
0178 edm::ParameterDescription<int>("yMax", 99999, true) and
0179 edm::ParameterDescription<std::string>("funcType", "pol0", true) and
0180 edm::ParameterDescription<std::vector<double>>("funcParams", {0.}, true)) or
0181 "AbsEtaClusEt" >> (edm::ParameterDescription<double>("xMin", 0.0, true) and
0182 edm::ParameterDescription<double>("xMax", 3.0, true) and
0183 edm::ParameterDescription<int>("yMin", 0, true) and
0184 edm::ParameterDescription<int>("yMax", 99999, true) and
0185 edm::ParameterDescription<std::string>("funcType", "pol0", true) and
0186 edm::ParameterDescription<std::vector<double>>("funcParams", {0.}, true));
0187
0188 binParamDesc.ifValue(edm::ParameterDescription<std::string>("binType", "AbsEtaClus", true), std::move(binDescCases));
0189
0190 varParamDesc.addVPSet("bins", binParamDesc);
0191 desc.add("dPhi1SParams", varParamDesc);
0192 desc.add("dPhi2SParams", varParamDesc);
0193 desc.add("dRZ2SParams", varParamDesc);
0194 desc.add<int>("productsToWrite", 0);
0195 descriptions.add(("hltEgammaHLTPixelMatchVarProducer"), desc);
0196 }
0197
0198 void EgammaHLTPixelMatchVarProducer::produce(edm::Event& iEvent, const edm::EventSetup&) {
0199
0200 auto recoEcalCandHandle = iEvent.getHandle(recoEcalCandidateToken_);
0201
0202 auto pixelSeedsHandle = iEvent.getHandle(pixelSeedsToken_);
0203
0204 if (!recoEcalCandHandle.isValid())
0205 return;
0206 else if (!pixelSeedsHandle.isValid()) {
0207 auto s2Map = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0208 for (unsigned int candNr = 0; candNr < recoEcalCandHandle->size(); candNr++) {
0209 reco::RecoEcalCandidateRef candRef(recoEcalCandHandle, candNr);
0210 s2Map->insert(candRef, 0);
0211 }
0212 iEvent.put(std::move(s2Map), "s2");
0213 return;
0214 }
0215
0216 auto dPhi1BestS2Map = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0217 auto dPhi2BestS2Map = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0218 auto dzBestS2Map = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0219 auto s2Map = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0220
0221 auto nrClusMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0222 auto seedClusEFracMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0223 auto phiWidthMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0224 auto etaWidthMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0225
0226 std::vector<PixelData> pixelData;
0227 for (size_t hitNr = 0; hitNr < nrHits_; hitNr++) {
0228 pixelData.emplace_back(PixelData("dPhi", hitNr, &reco::ElectronSeed::dPhiBest, recoEcalCandHandle));
0229 pixelData.emplace_back(PixelData("dRZ", hitNr, &reco::ElectronSeed::dRZBest, recoEcalCandHandle));
0230 }
0231
0232 for (unsigned int candNr = 0; candNr < recoEcalCandHandle->size(); candNr++) {
0233 reco::RecoEcalCandidateRef candRef(recoEcalCandHandle, candNr);
0234 reco::SuperClusterRef candSCRef = candRef->superCluster();
0235
0236 std::array<float, 4> bestS2{{std::numeric_limits<float>::max(),
0237 std::numeric_limits<float>::max(),
0238 std::numeric_limits<float>::max(),
0239 std::numeric_limits<float>::max()}};
0240 for (auto& seed : *pixelSeedsHandle) {
0241 const edm::RefToBase<reco::CaloCluster>& pixelClusterRef = seed.caloCluster();
0242 reco::SuperClusterRef pixelSCRef = pixelClusterRef.castTo<reco::SuperClusterRef>();
0243 if (&(*candSCRef) == &(*pixelSCRef)) {
0244 std::array<float, 4> s2Data = calS2(seed, -1);
0245 std::array<float, 4> s2DataPos = calS2(seed, +1);
0246 if (s2Data[0] < bestS2[0])
0247 bestS2 = s2Data;
0248 if (s2DataPos[0] < bestS2[0])
0249 bestS2 = s2DataPos;
0250
0251 if (productsToWrite_ >= 2) {
0252 for (auto& pixelDatum : pixelData) {
0253 pixelDatum.fill(seed);
0254 }
0255 }
0256 }
0257 }
0258
0259 s2Map->insert(candRef, bestS2[0]);
0260 if (productsToWrite_ >= 1) {
0261 dPhi1BestS2Map->insert(candRef, bestS2[1]);
0262 dPhi2BestS2Map->insert(candRef, bestS2[2]);
0263 dzBestS2Map->insert(candRef, bestS2[3]);
0264 }
0265 if (productsToWrite_ >= 2) {
0266 nrClusMap->insert(candRef, candSCRef->clustersSize());
0267 float seedClusEFrac = candSCRef->rawEnergy() > 0 ? candSCRef->seed()->energy() / candSCRef->rawEnergy() : 0.;
0268
0269
0270 seedClusEFracMap->insert(candRef, seedClusEFrac);
0271 float phiWidth = candSCRef->phiWidth();
0272 float etaWidth = candSCRef->etaWidth();
0273 phiWidthMap->insert(candRef, phiWidth);
0274 etaWidthMap->insert(candRef, etaWidth);
0275
0276 for (auto& pixelDatum : pixelData) {
0277 pixelDatum.fill(candRef);
0278 }
0279 }
0280 }
0281
0282 iEvent.put(std::move(s2Map), "s2");
0283 if (productsToWrite_ >= 1) {
0284 iEvent.put(std::move(dPhi1BestS2Map), "dPhi1BestS2");
0285 iEvent.put(std::move(dPhi2BestS2Map), "dPhi2BestS2");
0286 iEvent.put(std::move(dzBestS2Map), "dzBestS2");
0287 }
0288 if (productsToWrite_ >= 2) {
0289 for (auto& pixelDatum : pixelData) {
0290 pixelDatum.putInto(iEvent);
0291 }
0292 iEvent.put(std::move(nrClusMap), "nrClus");
0293 iEvent.put(std::move(seedClusEFracMap), "seedClusEFrac");
0294 iEvent.put(std::move(phiWidthMap), "phiWidth");
0295 iEvent.put(std::move(etaWidthMap), "etaWidth");
0296 }
0297 }
0298
0299 std::array<float, 4> EgammaHLTPixelMatchVarProducer::calS2(const reco::ElectronSeed& seed, int charge) const {
0300 const float dPhi1Const = dPhi1Para_(seed);
0301 const float dPhi2Const = dPhi2Para_(seed);
0302 const float dRZ2Const = dRZ2Para_(seed);
0303
0304 float dPhi1 = (charge < 0 ? seed.dPhiNeg(0) : seed.dPhiPos(0)) / dPhi1Const;
0305 float dPhi2 = (charge < 0 ? seed.dPhiNeg(1) : seed.dPhiPos(1)) / dPhi2Const;
0306 float dRz2 = (charge < 0 ? seed.dRZNeg(1) : seed.dRZPos(1)) / dRZ2Const;
0307
0308 float s2 = dPhi1 * dPhi1 + dPhi2 * dPhi2 + dRz2 * dRz2;
0309 return std::array<float, 4>{{s2, dPhi1, dPhi2, dRz2}};
0310 }
0311
0312 DEFINE_FWK_MODULE(EgammaHLTPixelMatchVarProducer);