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
0045       int layerOrDiskNr = seed.layerOrDiskNr(hitNr);
0046       if (layerOrDiskNr == 0 || layerOrDiskNr > 4) {
0047         if (layerOrDiskNr == std::numeric_limits<int>::max()) {
0048           throw cms::Exception("LogicError")
0049               << "The layerOrDiskNr of hitnr " << hitNr << " is " << layerOrDiskNr
0050               << " which is equal to numeric_limits::max which implies it was not filled correctly.\nWhile this is "
0051                  "valid for the matching variables as it can signal the failure to find a postive or negative seed "
0052                  "trajectory,\nit is not valid for the layer number as that is always known.\nThus there is a logic "
0053                  "error in the code or possible corrupt data";
0054         } else {
0055           throw cms::Exception("InvalidData") << "The layerOrDiskNr of hitnr " << hitNr << " is " << layerOrDiskNr
0056                                               << " and is not in the range of 1<=x<=4 which implies a new pixel "
0057                                                  "detector\nthis code does not support as the current one has only 4 "
0058                                                  "layers numbered 1-4 or corrupt data\n";
0059         }
0060       }
0061       int layerBit = 0x1 << layerOffset << layerOrDiskNr;
0062       info |= layerBit;
0063
0064       int nrLayersAlongTrajShifted = seed.nrLayersAlongTraj() << 12;
0065       info |= nrLayersAlongTrajShifted;
0066     }
0067     return info;
0068   }
0069
0070 }  // namespace
0071
0072 struct PixelData {
0073 public:
0074   PixelData(std::string name,
0075             size_t hitNr,
0076             float (reco::ElectronSeed::*func)(size_t) const,
0077             const edm::Handle<reco::RecoEcalCandidateCollection>& candHandle)
0078       : name_(std::move(name)), hitNr_(hitNr), func_(func), val_(std::numeric_limits<float>::max()), valInfo_(0) {
0079     valMap_ = std::make_unique<reco::RecoEcalCandidateIsolationMap>(candHandle);
0080     valInfoMap_ = std::make_unique<reco::RecoEcalCandidateIsolationMap>(candHandle);
0081   }
0082   PixelData(PixelData&& rhs) = default;
0083
0084   void resetVal() {
0085     val_ = std::numeric_limits<float>::max();
0086     valInfo_ = 0;
0087   }
0088   void fill(const reco::ElectronSeed& seed) {
0089     if (hitNr_ < seed.hitInfo().size()) {
0090       float seedVal = (seed.*func_)(hitNr_);
0091       if (std::abs(seedVal) < std::abs(val_)) {
0092         val_ = seedVal;
0093         valInfo_ = makeSeedInfo(seed);
0094       }
0095     }
0096   }
0097   void fill(const reco::RecoEcalCandidateRef& candRef) {
0098     valMap_->insert(candRef, val_);
0099     valInfoMap_->insert(candRef, valInfo_);
0100     val_ = std::numeric_limits<float>::max();
0101     valInfo_ = 0;
0102   }
0103
0104   void putInto(edm::Event& event) {
0105     event.put(std::move(valMap_), name_ + std::to_string(hitNr_ + 1));
0106     event.put(std::move(valInfoMap_), name_ + std::to_string(hitNr_ + 1) + "Info");
0107   }
0108
0109 private:
0110   std::unique_ptr<reco::RecoEcalCandidateIsolationMap> valMap_;
0111   std::unique_ptr<reco::RecoEcalCandidateIsolationMap> valInfoMap_;
0112   std::string name_;
0113   size_t hitNr_;
0114   float (reco::ElectronSeed::*func_)(size_t) const;
0115   float val_;
0116   float valInfo_;
0117 };
0118
0119 class EgammaHLTPixelMatchVarProducer : public edm::stream::EDProducer<> {
0120 public:
0121   explicit EgammaHLTPixelMatchVarProducer(const edm::ParameterSet&);
0122   ~EgammaHLTPixelMatchVarProducer() override;
0123
0124   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0125   void produce(edm::Event&, const edm::EventSetup&) override;
0126   std::array<float, 4> calS2(const reco::ElectronSeed& seed, int charge) const;
0127
0128 private:
0129   // ----------member data ---------------------------
0130
0131   const edm::EDGetTokenT<reco::RecoEcalCandidateCollection> recoEcalCandidateToken_;
0132   const edm::EDGetTokenT<reco::ElectronSeedCollection> pixelSeedsToken_;
0133
0134   egPM::Param<reco::ElectronSeed> dPhi1Para_;
0135   egPM::Param<reco::ElectronSeed> dPhi2Para_;
0136   egPM::Param<reco::ElectronSeed> dRZ2Para_;
0137
0138   int productsToWrite_;
0139   size_t nrHits_;
0140 };
0141
0142 EgammaHLTPixelMatchVarProducer::EgammaHLTPixelMatchVarProducer(const edm::ParameterSet& config)
0143     : recoEcalCandidateToken_(
0144           consumes<reco::RecoEcalCandidateCollection>(config.getParameter<edm::InputTag>("recoEcalCandidateProducer"))),
0145       pixelSeedsToken_(
0146           consumes<reco::ElectronSeedCollection>(config.getParameter<edm::InputTag>("pixelSeedsProducer"))),
0147       dPhi1Para_(config.getParameter<edm::ParameterSet>("dPhi1SParams")),
0148       dPhi2Para_(config.getParameter<edm::ParameterSet>("dPhi2SParams")),
0149       dRZ2Para_(config.getParameter<edm::ParameterSet>("dRZ2SParams")),
0150       productsToWrite_(config.getParameter<int>("productsToWrite")),
0151       nrHits_(4)
0152
0153 {
0154   //register your products
0155   produces<reco::RecoEcalCandidateIsolationMap>("s2");
0156   if (productsToWrite_ >= 1) {
0157     produces<reco::RecoEcalCandidateIsolationMap>("dPhi1BestS2");
0158     produces<reco::RecoEcalCandidateIsolationMap>("dPhi2BestS2");
0159     produces<reco::RecoEcalCandidateIsolationMap>("dzBestS2");
0160   }
0161   if (productsToWrite_ >= 2) {
0162     //note for product names we start from index 1
0163     for (size_t hitNr = 1; hitNr <= nrHits_; hitNr++) {
0164       produces<reco::RecoEcalCandidateIsolationMap>("dPhi" + std::to_string(hitNr));
0165       produces<reco::RecoEcalCandidateIsolationMap>("dPhi" + std::to_string(hitNr) + "Info");
0166       produces<reco::RecoEcalCandidateIsolationMap>("dRZ" + std::to_string(hitNr));
0167       produces<reco::RecoEcalCandidateIsolationMap>("dRZ" + std::to_string(hitNr) + "Info");
0168     }
0169     produces<reco::RecoEcalCandidateIsolationMap>("nrClus");
0170     produces<reco::RecoEcalCandidateIsolationMap>("seedClusEFrac");
0171     produces<reco::RecoEcalCandidateIsolationMap>("phiWidth");
0172     produces<reco::RecoEcalCandidateIsolationMap>("etaWidth");
0173   }
0174 }
0175
0176 EgammaHLTPixelMatchVarProducer::~EgammaHLTPixelMatchVarProducer() {}
0177
0178 void EgammaHLTPixelMatchVarProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0179   edm::ParameterSetDescription desc;
0180   desc.add<edm::InputTag>(("recoEcalCandidateProducer"), edm::InputTag("hltL1SeededRecoEcalCandidate"));
0181   desc.add<edm::InputTag>(("pixelSeedsProducer"), edm::InputTag("electronPixelSeeds"));
0182
0183   edm::ParameterSetDescription varParamDesc;
0184   edm::ParameterSetDescription binParamDesc;
0185
0186   auto binDescCases = "AbsEtaClus" >> (edm::ParameterDescription<double>("xMin", 0.0, true) and
0187                                        edm::ParameterDescription<double>("xMax", 3.0, true) and
0188                                        edm::ParameterDescription<int>("yMin", 0, true) and
0189                                        edm::ParameterDescription<int>("yMax", 99999, true) and
0190                                        edm::ParameterDescription<std::string>("funcType", "pol0", true) and
0191                                        edm::ParameterDescription<std::vector<double>>("funcParams", {0.}, true)) or
0192                       "AbsEtaClusPhi" >> (edm::ParameterDescription<double>("xMin", 0.0, true) and
0193                                           edm::ParameterDescription<double>("xMax", 3.0, true) and
0194                                           edm::ParameterDescription<int>("yMin", 0, true) and
0195                                           edm::ParameterDescription<int>("yMax", 99999, true) and
0196                                           edm::ParameterDescription<std::string>("funcType", "pol0", true) and
0197                                           edm::ParameterDescription<std::vector<double>>("funcParams", {0.}, true)) or
0198                       "AbsEtaClusEt" >> (edm::ParameterDescription<double>("xMin", 0.0, true) and
0199                                          edm::ParameterDescription<double>("xMax", 3.0, true) and
0200                                          edm::ParameterDescription<int>("yMin", 0, true) and
0201                                          edm::ParameterDescription<int>("yMax", 99999, true) and
0202                                          edm::ParameterDescription<std::string>("funcType", "pol0", true) and
0203                                          edm::ParameterDescription<std::vector<double>>("funcParams", {0.}, true));
0204
0205   binParamDesc.ifValue(edm::ParameterDescription<std::string>("binType", "AbsEtaClus", true), std::move(binDescCases));
0206
0207   varParamDesc.addVPSet("bins", binParamDesc);
0208   desc.add("dPhi1SParams", varParamDesc);
0209   desc.add("dPhi2SParams", varParamDesc);
0210   desc.add("dRZ2SParams", varParamDesc);
0211   desc.add<int>("productsToWrite", 0);
0212   descriptions.add(("hltEgammaHLTPixelMatchVarProducer"), desc);
0213 }
0214
0215 void EgammaHLTPixelMatchVarProducer::produce(edm::Event& iEvent, const edm::EventSetup&) {
0216   // Get the HLT filtered objects
0217   auto recoEcalCandHandle = iEvent.getHandle(recoEcalCandidateToken_);
0218
0219   auto pixelSeedsHandle = iEvent.getHandle(pixelSeedsToken_);
0220
0221   if (!recoEcalCandHandle.isValid())
0222     return;
0223   else if (!pixelSeedsHandle.isValid()) {
0224     auto s2Map = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0225     for (unsigned int candNr = 0; candNr < recoEcalCandHandle->size(); candNr++) {
0226       reco::RecoEcalCandidateRef candRef(recoEcalCandHandle, candNr);
0227       s2Map->insert(candRef, 0);
0228     }
0229     iEvent.put(std::move(s2Map), "s2");
0230     return;
0231   }
0232
0233   auto dPhi1BestS2Map = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0234   auto dPhi2BestS2Map = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0235   auto dzBestS2Map = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0236   auto s2Map = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0237
0238   auto nrClusMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0239   auto seedClusEFracMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0240   auto phiWidthMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0241   auto etaWidthMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
0242
0243   std::vector<PixelData> pixelData;
0244   for (size_t hitNr = 0; hitNr < nrHits_; hitNr++) {
0245     pixelData.emplace_back(PixelData("dPhi", hitNr, &reco::ElectronSeed::dPhiBest, recoEcalCandHandle));
0246     pixelData.emplace_back(PixelData("dRZ", hitNr, &reco::ElectronSeed::dRZBest, recoEcalCandHandle));
0247   }
0248
0249   for (unsigned int candNr = 0; candNr < recoEcalCandHandle->size(); candNr++) {
0250     reco::RecoEcalCandidateRef candRef(recoEcalCandHandle, candNr);
0251     reco::SuperClusterRef candSCRef = candRef->superCluster();
0252
0253     std::array<float, 4> bestS2{{std::numeric_limits<float>::max(),
0254                                  std::numeric_limits<float>::max(),
0255                                  std::numeric_limits<float>::max(),
0256                                  std::numeric_limits<float>::max()}};
0257     for (auto& seed : *pixelSeedsHandle) {
0258       const edm::RefToBase<reco::CaloCluster>& pixelClusterRef = seed.caloCluster();
0259       reco::SuperClusterRef pixelSCRef = pixelClusterRef.castTo<reco::SuperClusterRef>();
0260       if (&(*candSCRef) == &(*pixelSCRef)) {
0261         std::array<float, 4> s2Data = calS2(seed, -1);
0262         std::array<float, 4> s2DataPos = calS2(seed, +1);
0263         if (s2Data[0] < bestS2[0])
0264           bestS2 = s2Data;
0265         if (s2DataPos[0] < bestS2[0])
0266           bestS2 = s2DataPos;
0267
0268         if (productsToWrite_ >= 2) {
0269           for (auto& pixelDatum : pixelData) {
0270             pixelDatum.fill(seed);
0271           }
0272         }
0273       }
0274     }
0275
0276     s2Map->insert(candRef, bestS2[0]);
0277     if (productsToWrite_ >= 1) {
0278       dPhi1BestS2Map->insert(candRef, bestS2[1]);
0279       dPhi2BestS2Map->insert(candRef, bestS2[2]);
0280       dzBestS2Map->insert(candRef, bestS2[3]);
0281     }
0282     if (productsToWrite_ >= 2) {
0283       nrClusMap->insert(candRef, candSCRef->clustersSize());
0284       float seedClusEFrac = candSCRef->rawEnergy() > 0 ? candSCRef->seed()->energy() / candSCRef->rawEnergy() : 0.;
0285       //       std::cout <<"cand "<<candSCRef->energy()<<" E Corr "<<candSCRef->correctedEnergyUncertainty()<<" "<<candSCRef->correctedEnergy()<<" width "<<candSCRef->phiWidth()<<std::endl;
0286       //  float seedClusEFrac = candSCRef->phiWidth();
0287       seedClusEFracMap->insert(candRef, seedClusEFrac);
0288       float phiWidth = candSCRef->phiWidth();
0289       float etaWidth = candSCRef->etaWidth();
0290       phiWidthMap->insert(candRef, phiWidth);
0291       etaWidthMap->insert(candRef, etaWidth);
0292
0293       for (auto& pixelDatum : pixelData) {
0294         pixelDatum.fill(candRef);
0295       }
0296     }
0297   }
0298
0299   iEvent.put(std::move(s2Map), "s2");
0300   if (productsToWrite_ >= 1) {
0301     iEvent.put(std::move(dPhi1BestS2Map), "dPhi1BestS2");
0302     iEvent.put(std::move(dPhi2BestS2Map), "dPhi2BestS2");
0303     iEvent.put(std::move(dzBestS2Map), "dzBestS2");
0304   }
0305   if (productsToWrite_ >= 2) {
0306     for (auto& pixelDatum : pixelData) {
0307       pixelDatum.putInto(iEvent);
0308     }
0309     iEvent.put(std::move(nrClusMap), "nrClus");
0310     iEvent.put(std::move(seedClusEFracMap), "seedClusEFrac");
0311     iEvent.put(std::move(phiWidthMap), "phiWidth");
0312     iEvent.put(std::move(etaWidthMap), "etaWidth");
0313   }
0314 }
0315
0316 std::array<float, 4> EgammaHLTPixelMatchVarProducer::calS2(const reco::ElectronSeed& seed, int charge) const {
0317   const float dPhi1Const = dPhi1Para_(seed);
0318   const float dPhi2Const = dPhi2Para_(seed);
0319   const float dRZ2Const = dRZ2Para_(seed);
0320
0321   float dPhi1 = (charge < 0 ? seed.dPhiNeg(0) : seed.dPhiPos(0)) / dPhi1Const;
0322   float dPhi2 = (charge < 0 ? seed.dPhiNeg(1) : seed.dPhiPos(1)) / dPhi2Const;
0323   float dRz2 = (charge < 0 ? seed.dRZNeg(1) : seed.dRZPos(1)) / dRZ2Const;
0324
0325   float s2 = dPhi1 * dPhi1 + dPhi2 * dPhi2 + dRz2 * dRz2;
0326   return std::array<float, 4>{{s2, dPhi1, dPhi2, dRz2}};
0327 }
0328
0329 DEFINE_FWK_MODULE(EgammaHLTPixelMatchVarProducer);