EgammaHLTProducers/plugins/EgammaHLTPixelMatchVarProducer.cc

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   //first 4 bits are sub detect of each hit (0=barrel, 1 = endcap)
0033   //next 8 bits are layer information (0=no hit, 1 = hit), first 4 are barrel, next 4 are endcap (theres an empty bit here
0034   //next 4 bits are nr of layers info
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 }  // namespace
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   // ----------member data ---------------------------
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   //register your products
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     //note for product names we start from index 1
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   // Get the HLT filtered objects
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       //       std::cout <<"cand "<<candSCRef->energy()<<" E Corr "<<candSCRef->correctedEnergyUncertainty()<<" "<<candSCRef->correctedEnergy()<<" width "<<candSCRef->phiWidth()<<std::endl;
0269       //  float seedClusEFrac = candSCRef->phiWidth();
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);