Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-07-14 00:48:25

0001 // -*- C++ -*-
0002 //
0003 // Package:    HLTrigger/Egamma
0004 // Class:      HLTScoutingEgammaProducer
0005 //
0006 /**\class HLTScoutingEgammaProducer HLTScoutingEgammaProducer.cc HLTrigger/Egamma/src/HLTScoutingEgammaProducer.cc
0007 
0008 Description: Producer for Run3ScoutingElectron and Run3ScoutingPhoton
0009 
0010 */
0011 //
0012 // Original Author:  David G. Sheffield (Rutgers)
0013 //         Created:  Mon, 20 Jul 2015
0014 //
0015 //
0016 
0017 #include "HLTScoutingEgammaProducer.h"
0018 
0019 #include <cstdint>
0020 
0021 // function to find rechhit associated to detid and return energy
0022 float recHitE(const DetId id, const EcalRecHitCollection& recHits) {
0023   if (id == DetId(0)) {
0024     return 0;
0025   } else {
0026     EcalRecHitCollection::const_iterator it = recHits.find(id);
0027     if (it != recHits.end())
0028       return (*it).energy();
0029   }
0030   return 0;
0031 }
0032 
0033 float recHitT(const DetId id, const EcalRecHitCollection& recHits) {
0034   if (id == DetId(0)) {
0035     return 0;
0036   } else {
0037     EcalRecHitCollection::const_iterator it = recHits.find(id);
0038     if (it != recHits.end())
0039       return (*it).time();
0040   }
0041   return 0;
0042 }
0043 
0044 //
0045 // constructors and destructor
0046 //
0047 HLTScoutingEgammaProducer::HLTScoutingEgammaProducer(const edm::ParameterSet& iConfig)
0048     : EgammaCandidateCollection_(
0049           consumes<reco::RecoEcalCandidateCollection>(iConfig.getParameter<edm::InputTag>("EgammaCandidates"))),
0050       EgammaGsfTrackCollection_(
0051           consumes<reco::GsfTrackCollection>(iConfig.getParameter<edm::InputTag>("EgammaGsfTracks"))),
0052       SigmaIEtaIEtaMap_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("SigmaIEtaIEtaMap"))),
0053       R9Map_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("r9Map"))),
0054       HoverEMap_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("HoverEMap"))),
0055       DetaMap_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("DetaMap"))),
0056       DphiMap_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("DphiMap"))),
0057       MissingHitsMap_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("MissingHitsMap"))),
0058       OneOEMinusOneOPMap_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("OneOEMinusOneOPMap"))),
0059       EcalPFClusterIsoMap_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("EcalPFClusterIsoMap"))),
0060       EleGsfTrackIsoMap_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("EleGsfTrackIsoMap"))),
0061       HcalPFClusterIsoMap_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("HcalPFClusterIsoMap"))),
0062       egammaPtCut(iConfig.getParameter<double>("egammaPtCut")),
0063       egammaEtaCut(iConfig.getParameter<double>("egammaEtaCut")),
0064       egammaHoverECut(iConfig.getParameter<double>("egammaHoverECut")),
0065       egammaSigmaIEtaIEtaCut(iConfig.getParameter<std::vector<double>>("egammaSigmaIEtaIEtaCut")),
0066       absEtaBinUpperEdges(iConfig.getParameter<std::vector<double>>("absEtaBinUpperEdges")),
0067       mantissaPrecision(iConfig.getParameter<int>("mantissaPrecision")),
0068       saveRecHitTiming(iConfig.getParameter<bool>("saveRecHitTiming")),
0069       rechitMatrixSize(iConfig.getParameter<int>("rechitMatrixSize")),  //(2n+1)^2
0070       rechitZeroSuppression(iConfig.getParameter<bool>("rechitZeroSuppression")),
0071       ecalRechitEB_(consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("ecalRechitEB"))),
0072       ecalRechitEE_(consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("ecalRechitEE"))) {
0073   // cross-check for compatibility in input vectors
0074   if (absEtaBinUpperEdges.size() != egammaSigmaIEtaIEtaCut.size()) {
0075     throw cms::Exception("IncompatibleVects")
0076         << "size of \"absEtaBinUpperEdges\" (" << absEtaBinUpperEdges.size() << ") and \"egammaSigmaIEtaIEtaCut\" ("
0077         << egammaSigmaIEtaIEtaCut.size() << ") differ";
0078   }
0079 
0080   for (auto aIt = 1u; aIt < absEtaBinUpperEdges.size(); ++aIt) {
0081     if (absEtaBinUpperEdges[aIt - 1] < 0 || absEtaBinUpperEdges[aIt] < 0) {
0082       throw cms::Exception("IncorrectValue") << "absEtaBinUpperEdges entries should be greater than or equal to zero.";
0083     }
0084     if (absEtaBinUpperEdges[aIt - 1] >= absEtaBinUpperEdges[aIt]) {
0085       throw cms::Exception("ImproperBinning") << "absEtaBinUpperEdges entries should be in increasing order.";
0086     }
0087   }
0088 
0089   if (not absEtaBinUpperEdges.empty() and absEtaBinUpperEdges[absEtaBinUpperEdges.size() - 1] < egammaEtaCut) {
0090     throw cms::Exception("IncorrectValue")
0091         << "Last entry in \"absEtaBinUpperEdges\" (" << absEtaBinUpperEdges[absEtaBinUpperEdges.size() - 1]
0092         << ") should have a value larger than \"egammaEtaCut\" (" << egammaEtaCut << ").";
0093   }
0094 
0095   //register products
0096   produces<Run3ScoutingElectronCollection>();
0097   produces<Run3ScoutingPhotonCollection>();
0098   topologyToken_ = esConsumes();
0099 }
0100 
0101 HLTScoutingEgammaProducer::~HLTScoutingEgammaProducer() = default;
0102 
0103 // ------------ method called to produce the data  ------------
0104 void HLTScoutingEgammaProducer::produce(edm::StreamID sid, edm::Event& iEvent, edm::EventSetup const& setup) const {
0105   using namespace edm;
0106 
0107   auto outElectrons = std::make_unique<Run3ScoutingElectronCollection>();
0108   auto outPhotons = std::make_unique<Run3ScoutingPhotonCollection>();
0109 
0110   // Get RecoEcalCandidate
0111   Handle<reco::RecoEcalCandidateCollection> EgammaCandidateCollection;
0112   if (!iEvent.getByToken(EgammaCandidateCollection_, EgammaCandidateCollection)) {
0113     iEvent.put(std::move(outElectrons));
0114     iEvent.put(std::move(outPhotons));
0115     return;
0116   }
0117 
0118   // Get GsfTrack
0119   Handle<reco::GsfTrackCollection> EgammaGsfTrackCollection;
0120   if (!iEvent.getByToken(EgammaGsfTrackCollection_, EgammaGsfTrackCollection)) {
0121     iEvent.put(std::move(outElectrons));
0122     iEvent.put(std::move(outPhotons));
0123     return;
0124   }
0125 
0126   // Get SigmaIEtaIEtaMap
0127   Handle<RecoEcalCandMap> SigmaIEtaIEtaMap;
0128   if (!iEvent.getByToken(SigmaIEtaIEtaMap_, SigmaIEtaIEtaMap)) {
0129     iEvent.put(std::move(outElectrons));
0130     iEvent.put(std::move(outPhotons));
0131     return;
0132   }
0133 
0134   // Get R9Map
0135   Handle<RecoEcalCandMap> R9Map;
0136   if (!iEvent.getByToken(R9Map_, R9Map)) {
0137     iEvent.put(std::move(outElectrons));
0138     iEvent.put(std::move(outPhotons));
0139     return;
0140   }
0141 
0142   // Get HoverEMap
0143   Handle<RecoEcalCandMap> HoverEMap;
0144   if (!iEvent.getByToken(HoverEMap_, HoverEMap)) {
0145     iEvent.put(std::move(outElectrons));
0146     iEvent.put(std::move(outPhotons));
0147     return;
0148   }
0149 
0150   // Get DetaMap
0151   Handle<RecoEcalCandMap> DetaMap;
0152   if (!iEvent.getByToken(DetaMap_, DetaMap)) {
0153     iEvent.put(std::move(outElectrons));
0154     iEvent.put(std::move(outPhotons));
0155     return;
0156   }
0157 
0158   // Get DphiMap
0159   Handle<RecoEcalCandMap> DphiMap;
0160   if (!iEvent.getByToken(DphiMap_, DphiMap)) {
0161     iEvent.put(std::move(outElectrons));
0162     iEvent.put(std::move(outPhotons));
0163     return;
0164   }
0165 
0166   // Get MissingHitsMap
0167   Handle<RecoEcalCandMap> MissingHitsMap;
0168   if (!iEvent.getByToken(MissingHitsMap_, MissingHitsMap)) {
0169     iEvent.put(std::move(outElectrons));
0170     iEvent.put(std::move(outPhotons));
0171     return;
0172   }
0173 
0174   // Get 1/E - 1/p Map
0175   Handle<RecoEcalCandMap> OneOEMinusOneOPMap;
0176   if (!iEvent.getByToken(OneOEMinusOneOPMap_, OneOEMinusOneOPMap)) {
0177     iEvent.put(std::move(outElectrons));
0178     iEvent.put(std::move(outPhotons));
0179     return;
0180   }
0181 
0182   // Get EcalPFClusterIsoMap
0183   Handle<RecoEcalCandMap> EcalPFClusterIsoMap;
0184   if (!iEvent.getByToken(EcalPFClusterIsoMap_, EcalPFClusterIsoMap)) {
0185     iEvent.put(std::move(outElectrons));
0186     iEvent.put(std::move(outPhotons));
0187     return;
0188   }
0189 
0190   // Get EleGsfTrackIsoMap
0191   Handle<RecoEcalCandMap> EleGsfTrackIsoMap;
0192   if (!iEvent.getByToken(EleGsfTrackIsoMap_, EleGsfTrackIsoMap)) {
0193     iEvent.put(std::move(outElectrons));
0194     iEvent.put(std::move(outPhotons));
0195     return;
0196   }
0197 
0198   // Get HcalPFClusterIsoMap
0199   Handle<RecoEcalCandMap> HcalPFClusterIsoMap;
0200   if (!iEvent.getByToken(HcalPFClusterIsoMap_, HcalPFClusterIsoMap)) {
0201     iEvent.put(std::move(outElectrons));
0202     iEvent.put(std::move(outPhotons));
0203     return;
0204   }
0205 
0206   edm::Handle<EcalRecHitCollection> rechitsEB;
0207   edm::Handle<EcalRecHitCollection> rechitsEE;
0208   iEvent.getByToken(ecalRechitEB_, rechitsEB);
0209   iEvent.getByToken(ecalRechitEE_, rechitsEE);
0210 
0211   const CaloTopology* topology = &setup.getData(topologyToken_);
0212 
0213   // Produce electrons and photons
0214   int index = 0;
0215   for (auto& candidate : *EgammaCandidateCollection) {
0216     reco::RecoEcalCandidateRef candidateRef = getRef(EgammaCandidateCollection, index);
0217     ++index;
0218     if (candidateRef.isNull() && !candidateRef.isAvailable())
0219       continue;
0220 
0221     if (candidate.pt() < egammaPtCut)
0222       continue;
0223     if (fabs(candidate.eta()) > egammaEtaCut)
0224       continue;
0225 
0226     reco::SuperClusterRef scRef = candidate.superCluster();
0227     if (scRef.isNull() && !scRef.isAvailable())
0228       continue;
0229 
0230     reco::CaloClusterPtr SCseed = candidate.superCluster()->seed();
0231     const EcalRecHitCollection* rechits = (std::abs(scRef->eta()) < 1.479) ? rechitsEB.product() : rechitsEE.product();
0232     Cluster2ndMoments moments = EcalClusterTools::cluster2ndMoments(*SCseed, *rechits);
0233     float sMin = moments.sMin;
0234     float sMaj = moments.sMaj;
0235 
0236     uint32_t seedId = (*SCseed).seed();
0237 
0238     std::vector<DetId> mDetIds = EcalClusterTools::matrixDetId((topology), (*SCseed).seed(), rechitMatrixSize);
0239 
0240     int detSize = mDetIds.size();
0241     std::vector<uint32_t> mDetIdIds;
0242     std::vector<float> mEnergies;
0243     std::vector<float> mTimes;
0244     mDetIdIds.reserve(detSize);
0245     mEnergies.reserve(detSize);
0246     mTimes.reserve(detSize);
0247 
0248     for (int i = 0; i < detSize; i++) {
0249       auto const recHit_en = recHitE(mDetIds[i], *rechits);
0250       if (not rechitZeroSuppression or recHit_en > 0) {
0251         mDetIdIds.push_back(mDetIds[i]);
0252         mEnergies.push_back(MiniFloatConverter::reduceMantissaToNbitsRounding(recHit_en, mantissaPrecision));
0253         if (saveRecHitTiming) {
0254           mTimes.push_back(
0255               MiniFloatConverter::reduceMantissaToNbitsRounding(recHitT(mDetIds[i], *rechits), mantissaPrecision));
0256         }
0257       }
0258     }
0259 
0260     auto const HoE = candidate.superCluster()->energy() != 0.
0261                          ? ((*HoverEMap)[candidateRef] / candidate.superCluster()->energy())
0262                          : 999.;
0263     if (HoE > egammaHoverECut)
0264       continue;
0265 
0266     if (not absEtaBinUpperEdges.empty()) {
0267       auto const sinin = candidate.superCluster()->energy() != 0. ? (*SigmaIEtaIEtaMap)[candidateRef] : 999.;
0268       auto etaBinIdx = std::distance(
0269           absEtaBinUpperEdges.begin(),
0270           std::lower_bound(absEtaBinUpperEdges.begin(), absEtaBinUpperEdges.end(), std::abs(candidate.eta())));
0271 
0272       if (sinin > egammaSigmaIEtaIEtaCut[etaBinIdx])
0273         continue;
0274     }
0275 
0276     float d0 = 0.0;
0277     float dz = 0.0;
0278     int charge = -999;
0279     for (auto& track : *EgammaGsfTrackCollection) {
0280       RefToBase<TrajectorySeed> seed = track.extra()->seedRef();
0281       reco::ElectronSeedRef elseed = seed.castTo<reco::ElectronSeedRef>();
0282       RefToBase<reco::CaloCluster> caloCluster = elseed->caloCluster();
0283       reco::SuperClusterRef scRefFromTrk = caloCluster.castTo<reco::SuperClusterRef>();
0284       if (scRefFromTrk == scRef) {
0285         d0 = track.d0();
0286         dz = track.dz();
0287         charge = track.charge();
0288       }
0289     }
0290     if (charge == -999) {  // No associated track. Candidate is a scouting photon
0291       outPhotons->emplace_back(candidate.pt(),
0292                                candidate.eta(),
0293                                candidate.phi(),
0294                                candidate.mass(),
0295                                (*SigmaIEtaIEtaMap)[candidateRef],
0296                                HoE,
0297                                (*EcalPFClusterIsoMap)[candidateRef],
0298                                (*HcalPFClusterIsoMap)[candidateRef],
0299                                0.,
0300                                (*R9Map)[candidateRef],
0301                                sMin,
0302                                sMaj,
0303                                seedId,
0304                                mEnergies,
0305                                mDetIdIds,
0306                                mTimes,
0307                                rechitZeroSuppression);  //read for(ieta){for(iphi){}}
0308     } else {                                            // Candidate is a scouting electron
0309       outElectrons->emplace_back(candidate.pt(),
0310                                  candidate.eta(),
0311                                  candidate.phi(),
0312                                  candidate.mass(),
0313                                  d0,
0314                                  dz,
0315                                  (*DetaMap)[candidateRef],
0316                                  (*DphiMap)[candidateRef],
0317                                  (*SigmaIEtaIEtaMap)[candidateRef],
0318                                  HoE,
0319                                  (*OneOEMinusOneOPMap)[candidateRef],
0320                                  (*MissingHitsMap)[candidateRef],
0321                                  charge,
0322                                  (*EcalPFClusterIsoMap)[candidateRef],
0323                                  (*HcalPFClusterIsoMap)[candidateRef],
0324                                  (*EleGsfTrackIsoMap)[candidateRef],
0325                                  (*R9Map)[candidateRef],
0326                                  sMin,
0327                                  sMaj,
0328                                  seedId,
0329                                  mEnergies,
0330                                  mDetIdIds,
0331                                  mTimes,
0332                                  rechitZeroSuppression);  //read for(ieta){for(iphi){}}
0333     }
0334   }
0335 
0336   // Put output
0337   iEvent.put(std::move(outElectrons));
0338   iEvent.put(std::move(outPhotons));
0339 }
0340 
0341 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0342 void HLTScoutingEgammaProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0343   edm::ParameterSetDescription desc;
0344   desc.add<edm::InputTag>("EgammaCandidates", edm::InputTag("hltEgammaCandidates"));
0345   desc.add<edm::InputTag>("EgammaGsfTracks", edm::InputTag("hltEgammaGsfTracks"));
0346   desc.add<edm::InputTag>("SigmaIEtaIEtaMap", edm::InputTag("hltEgammaClusterShape:sigmaIEtaIEta5x5"));
0347   desc.add<edm::InputTag>("r9Map", edm::InputTag("hltEgammaR9ID:r95x5"));
0348   desc.add<edm::InputTag>("HoverEMap", edm::InputTag("hltEgammaHoverE"));
0349   desc.add<edm::InputTag>("DetaMap", edm::InputTag("hltEgammaGsfTrackVars:DetaSeed"));
0350   desc.add<edm::InputTag>("DphiMap", edm::InputTag("hltEgammaGsfTrackVars:Dphi"));
0351   desc.add<edm::InputTag>("MissingHitsMap", edm::InputTag("hltEgammaGsfTrackVars:MissingHits"));
0352   desc.add<edm::InputTag>("OneOEMinusOneOPMap", edm::InputTag("hltEgammaGsfTrackVars:OneOESuperMinusOneOP"));
0353   desc.add<edm::InputTag>("EcalPFClusterIsoMap", edm::InputTag("hltEgammaEcalPFClusterIso"));
0354   desc.add<edm::InputTag>("EleGsfTrackIsoMap", edm::InputTag("hltEgammaEleGsfTrackIso"));
0355   desc.add<edm::InputTag>("HcalPFClusterIsoMap", edm::InputTag("hltEgammaHcalPFClusterIso"));
0356   desc.add<double>("egammaPtCut", 4.0);
0357   desc.add<double>("egammaEtaCut", 2.5);
0358   desc.add<double>("egammaHoverECut", 1.0);
0359   desc.add<std::vector<double>>("egammaSigmaIEtaIEtaCut", {99999.0, 99999.0});
0360   desc.add<std::vector<double>>("absEtaBinUpperEdges", {1.479, 5.0});
0361   desc.add<bool>("saveRecHitTiming", false);
0362   desc.add<int>("mantissaPrecision", 10)->setComment("default float16, change to 23 for float32");
0363   desc.add<int>("rechitMatrixSize", 10);
0364   desc.add<bool>("rechitZeroSuppression", true);
0365   desc.add<edm::InputTag>("ecalRechitEB", edm::InputTag("hltEcalRecHit:EcalRecHitsEB"));
0366   desc.add<edm::InputTag>("ecalRechitEE", edm::InputTag("hltEcalRecHit:EcalRecHitsEE"));
0367   descriptions.add("hltScoutingEgammaProducer", desc);
0368 }
0369 
0370 // declare this class as a framework plugin
0371 #include "FWCore/Framework/interface/MakerMacros.h"
0372 DEFINE_FWK_MODULE(HLTScoutingEgammaProducer);