Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:18:22 

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       fBremMap_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("fBremMap"))),
0060       EcalPFClusterIsoMap_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("EcalPFClusterIsoMap"))),
0061       EleGsfTrackIsoMap_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("EleGsfTrackIsoMap"))),
0062       HcalPFClusterIsoMap_(consumes<RecoEcalCandMap>(iConfig.getParameter<edm::InputTag>("HcalPFClusterIsoMap"))),
0063       egammaPtCut(iConfig.getParameter<double>("egammaPtCut")),
0064       egammaEtaCut(iConfig.getParameter<double>("egammaEtaCut")),
0065       egammaHoverECut(iConfig.getParameter<double>("egammaHoverECut")),
0066       egammaSigmaIEtaIEtaCut(iConfig.getParameter<std::vector<double>>("egammaSigmaIEtaIEtaCut")),
0067       absEtaBinUpperEdges(iConfig.getParameter<std::vector<double>>("absEtaBinUpperEdges")),
0068       mantissaPrecision(iConfig.getParameter<int>("mantissaPrecision")),
0069       saveRecHitTiming(iConfig.getParameter<bool>("saveRecHitTiming")),
0070       rechitMatrixSize(iConfig.getParameter<int>("rechitMatrixSize")),  //(2n+1)^2
0071       rechitZeroSuppression(iConfig.getParameter<bool>("rechitZeroSuppression")),
0072       ecalRechitEB_(consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("ecalRechitEB"))),
0073       ecalRechitEE_(consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("ecalRechitEE"))) {
0074   // cross-check for compatibility in input vectors
0075   if (absEtaBinUpperEdges.size() != egammaSigmaIEtaIEtaCut.size()) {
0076     throw cms::Exception("IncompatibleVects")
0077         << "size of \"absEtaBinUpperEdges\" (" << absEtaBinUpperEdges.size() << ") and \"egammaSigmaIEtaIEtaCut\" ("
0078         << egammaSigmaIEtaIEtaCut.size() << ") differ";
0079   }
0080 
0081   for (auto aIt = 1u; aIt < absEtaBinUpperEdges.size(); ++aIt) {
0082     if (absEtaBinUpperEdges[aIt - 1] < 0 || absEtaBinUpperEdges[aIt] < 0) {
0083       throw cms::Exception("IncorrectValue") << "absEtaBinUpperEdges entries should be greater than or equal to zero.";
0084     }
0085     if (absEtaBinUpperEdges[aIt - 1] >= absEtaBinUpperEdges[aIt]) {
0086       throw cms::Exception("ImproperBinning") << "absEtaBinUpperEdges entries should be in increasing order.";
0087     }
0088   }
0089 
0090   if (not absEtaBinUpperEdges.empty() and absEtaBinUpperEdges[absEtaBinUpperEdges.size() - 1] < egammaEtaCut) {
0091     throw cms::Exception("IncorrectValue")
0092         << "Last entry in \"absEtaBinUpperEdges\" (" << absEtaBinUpperEdges[absEtaBinUpperEdges.size() - 1]
0093         << ") should have a value larger than \"egammaEtaCut\" (" << egammaEtaCut << ").";
0094   }
0095 
0096   //register products
0097   produces<Run3ScoutingElectronCollection>();
0098   produces<Run3ScoutingPhotonCollection>();
0099   topologyToken_ = esConsumes();
0100 }
0101 
0102 HLTScoutingEgammaProducer::~HLTScoutingEgammaProducer() = default;
0103 
0104 // ------------ method called to produce the data  ------------
0105 void HLTScoutingEgammaProducer::produce(edm::StreamID sid, edm::Event& iEvent, edm::EventSetup const& setup) const {
0106   using namespace edm;
0107 
0108   auto outElectrons = std::make_unique<Run3ScoutingElectronCollection>();
0109   auto outPhotons = std::make_unique<Run3ScoutingPhotonCollection>();
0110 
0111   // Get RecoEcalCandidate
0112   Handle<reco::RecoEcalCandidateCollection> EgammaCandidateCollection;
0113   if (!iEvent.getByToken(EgammaCandidateCollection_, EgammaCandidateCollection)) {
0114     iEvent.put(std::move(outElectrons));
0115     iEvent.put(std::move(outPhotons));
0116     return;
0117   }
0118 
0119   // Get GsfTrack
0120   Handle<reco::GsfTrackCollection> EgammaGsfTrackCollection;
0121   if (!iEvent.getByToken(EgammaGsfTrackCollection_, EgammaGsfTrackCollection)) {
0122     iEvent.put(std::move(outElectrons));
0123     iEvent.put(std::move(outPhotons));
0124     return;
0125   }
0126 
0127   // Get SigmaIEtaIEtaMap
0128   Handle<RecoEcalCandMap> SigmaIEtaIEtaMap;
0129   if (!iEvent.getByToken(SigmaIEtaIEtaMap_, SigmaIEtaIEtaMap)) {
0130     iEvent.put(std::move(outElectrons));
0131     iEvent.put(std::move(outPhotons));
0132     return;
0133   }
0134 
0135   // Get R9Map
0136   Handle<RecoEcalCandMap> R9Map;
0137   if (!iEvent.getByToken(R9Map_, R9Map)) {
0138     iEvent.put(std::move(outElectrons));
0139     iEvent.put(std::move(outPhotons));
0140     return;
0141   }
0142 
0143   // Get HoverEMap
0144   Handle<RecoEcalCandMap> HoverEMap;
0145   if (!iEvent.getByToken(HoverEMap_, HoverEMap)) {
0146     iEvent.put(std::move(outElectrons));
0147     iEvent.put(std::move(outPhotons));
0148     return;
0149   }
0150 
0151   // Get DetaMap
0152   Handle<RecoEcalCandMap> DetaMap;
0153   if (!iEvent.getByToken(DetaMap_, DetaMap)) {
0154     iEvent.put(std::move(outElectrons));
0155     iEvent.put(std::move(outPhotons));
0156     return;
0157   }
0158 
0159   // Get DphiMap
0160   Handle<RecoEcalCandMap> DphiMap;
0161   if (!iEvent.getByToken(DphiMap_, DphiMap)) {
0162     iEvent.put(std::move(outElectrons));
0163     iEvent.put(std::move(outPhotons));
0164     return;
0165   }
0166 
0167   // Get MissingHitsMap
0168   Handle<RecoEcalCandMap> MissingHitsMap;
0169   if (!iEvent.getByToken(MissingHitsMap_, MissingHitsMap)) {
0170     iEvent.put(std::move(outElectrons));
0171     iEvent.put(std::move(outPhotons));
0172     return;
0173   }
0174 
0175   // Get 1/E - 1/p Map
0176   Handle<RecoEcalCandMap> OneOEMinusOneOPMap;
0177   if (!iEvent.getByToken(OneOEMinusOneOPMap_, OneOEMinusOneOPMap)) {
0178     iEvent.put(std::move(outElectrons));
0179     iEvent.put(std::move(outPhotons));
0180     return;
0181   }
0182 
0183   // Get fBrem Map
0184   Handle<RecoEcalCandMap> fBremMap;
0185   if (!iEvent.getByToken(fBremMap_, fBremMap)) {
0186     iEvent.put(std::move(outElectrons));
0187     iEvent.put(std::move(outPhotons));
0188     return;
0189   }
0190 
0191   // Get EcalPFClusterIsoMap
0192   Handle<RecoEcalCandMap> EcalPFClusterIsoMap;
0193   if (!iEvent.getByToken(EcalPFClusterIsoMap_, EcalPFClusterIsoMap)) {
0194     iEvent.put(std::move(outElectrons));
0195     iEvent.put(std::move(outPhotons));
0196     return;
0197   }
0198 
0199   // Get EleGsfTrackIsoMap
0200   Handle<RecoEcalCandMap> EleGsfTrackIsoMap;
0201   if (!iEvent.getByToken(EleGsfTrackIsoMap_, EleGsfTrackIsoMap)) {
0202     iEvent.put(std::move(outElectrons));
0203     iEvent.put(std::move(outPhotons));
0204     return;
0205   }
0206 
0207   // Get HcalPFClusterIsoMap
0208   Handle<RecoEcalCandMap> HcalPFClusterIsoMap;
0209   if (!iEvent.getByToken(HcalPFClusterIsoMap_, HcalPFClusterIsoMap)) {
0210     iEvent.put(std::move(outElectrons));
0211     iEvent.put(std::move(outPhotons));
0212     return;
0213   }
0214 
0215   edm::Handle<EcalRecHitCollection> rechitsEB;
0216   edm::Handle<EcalRecHitCollection> rechitsEE;
0217   iEvent.getByToken(ecalRechitEB_, rechitsEB);
0218   iEvent.getByToken(ecalRechitEE_, rechitsEE);
0219 
0220   const CaloTopology* topology = &setup.getData(topologyToken_);
0221 
0222   // Produce electrons and photons
0223   int index = 0;
0224   for (auto& candidate : *EgammaCandidateCollection) {
0225     reco::RecoEcalCandidateRef candidateRef = getRef(EgammaCandidateCollection, index);
0226     ++index;
0227     if (candidateRef.isNull() && !candidateRef.isAvailable())
0228       continue;
0229 
0230     if (candidate.pt() < egammaPtCut)
0231       continue;
0232     if (fabs(candidate.eta()) > egammaEtaCut)
0233       continue;
0234 
0235     reco::SuperClusterRef scRef = candidate.superCluster();
0236     if (scRef.isNull() && !scRef.isAvailable())
0237       continue;
0238 
0239     reco::CaloClusterPtr SCseed = candidate.superCluster()->seed();
0240     const EcalRecHitCollection* rechits = (std::abs(scRef->eta()) < 1.479) ? rechitsEB.product() : rechitsEE.product();
0241     Cluster2ndMoments moments = EcalClusterTools::cluster2ndMoments(*SCseed, *rechits);
0242     float sMin = moments.sMin;
0243     float sMaj = moments.sMaj;
0244 
0245     uint32_t seedId = (*SCseed).seed();
0246 
0247     std::vector<DetId> mDetIds = EcalClusterTools::matrixDetId((topology), (*SCseed).seed(), rechitMatrixSize);
0248 
0249     int detSize = mDetIds.size();
0250     std::vector<uint32_t> mDetIdIds;
0251     std::vector<float> mEnergies;
0252     std::vector<float> mTimes;
0253     mDetIdIds.reserve(detSize);
0254     mEnergies.reserve(detSize);
0255     mTimes.reserve(detSize);
0256 
0257     for (int i = 0; i < detSize; i++) {
0258       auto const recHit_en = recHitE(mDetIds[i], *rechits);
0259       if (not rechitZeroSuppression or recHit_en > 0) {
0260         mDetIdIds.push_back(mDetIds[i]);
0261         mEnergies.push_back(MiniFloatConverter::reduceMantissaToNbitsRounding(recHit_en, mantissaPrecision));
0262         if (saveRecHitTiming) {
0263           mTimes.push_back(
0264               MiniFloatConverter::reduceMantissaToNbitsRounding(recHitT(mDetIds[i], *rechits), mantissaPrecision));
0265         }
0266       }
0267     }
0268 
0269     auto const HoE = candidate.superCluster()->energy() != 0.
0270                          ? ((*HoverEMap)[candidateRef] / candidate.superCluster()->energy())
0271                          : 999.;
0272     if (HoE > egammaHoverECut)
0273       continue;
0274 
0275     if (not absEtaBinUpperEdges.empty()) {
0276       auto const sinin = candidate.superCluster()->energy() != 0. ? (*SigmaIEtaIEtaMap)[candidateRef] : 999.;
0277       auto etaBinIdx = std::distance(
0278           absEtaBinUpperEdges.begin(),
0279           std::lower_bound(absEtaBinUpperEdges.begin(), absEtaBinUpperEdges.end(), std::abs(candidate.eta())));
0280 
0281       if (sinin > egammaSigmaIEtaIEtaCut[etaBinIdx])
0282         continue;
0283     }
0284 
0285     unsigned int const maxTrkSize = EgammaGsfTrackCollection->size();
0286     std::vector<float> trkd0;
0287     std::vector<float> trkdz;
0288     std::vector<float> trkpt;
0289     std::vector<float> trketa;
0290     std::vector<float> trkphi;
0291     std::vector<float> trkpMode;
0292     std::vector<float> trketaMode;
0293     std::vector<float> trkphiMode;
0294     std::vector<float> trkqoverpModeError;
0295     std::vector<float> trkchi2overndf;
0296     std::vector<int> trkcharge;
0297     trkd0.reserve(maxTrkSize);
0298     trkdz.reserve(maxTrkSize);
0299     trkpt.reserve(maxTrkSize);
0300     trketa.reserve(maxTrkSize);
0301     trkphi.reserve(maxTrkSize);
0302     trkpMode.reserve(maxTrkSize);
0303     trketaMode.reserve(maxTrkSize);
0304     trkphiMode.reserve(maxTrkSize);
0305     trkqoverpModeError.reserve(maxTrkSize);
0306     trkchi2overndf.reserve(maxTrkSize);
0307     trkcharge.reserve(maxTrkSize);
0308 
0309     for (auto& track : *EgammaGsfTrackCollection) {
0310       RefToBase<TrajectorySeed> seed = track.extra()->seedRef();
0311       reco::ElectronSeedRef elseed = seed.castTo<reco::ElectronSeedRef>();
0312       RefToBase<reco::CaloCluster> caloCluster = elseed->caloCluster();
0313       reco::SuperClusterRef scRefFromTrk = caloCluster.castTo<reco::SuperClusterRef>();
0314       if (scRefFromTrk == scRef) {
0315         trkd0.push_back(track.d0());
0316         trkdz.push_back(track.dz());
0317         trkpt.push_back(track.pt());
0318         trketa.push_back(track.eta());
0319         trkphi.push_back(track.phi());
0320         trkpMode.push_back(track.pMode());
0321         trketaMode.push_back(track.etaMode());
0322         trkphiMode.push_back(track.phiMode());
0323         trkqoverpModeError.push_back(track.qoverpModeError());
0324         auto const trackndof = track.ndof();
0325         trkchi2overndf.push_back(((trackndof == 0) ? -1 : (track.chi2() / trackndof)));
0326         trkcharge.push_back(track.charge());
0327       }
0328     }
0329     if (trkcharge.empty()) {  // No associated track. Candidate is a scouting photon
0330       outPhotons->emplace_back(candidate.pt(),
0331                                candidate.eta(),
0332                                candidate.phi(),
0333                                candidate.mass(),
0334                                scRef->rawEnergy(),
0335                                scRef->preshowerEnergy(),
0336                                scRef->correctedEnergyUncertainty(),
0337                                (*SigmaIEtaIEtaMap)[candidateRef],
0338                                HoE,
0339                                (*EcalPFClusterIsoMap)[candidateRef],
0340                                (*HcalPFClusterIsoMap)[candidateRef],
0341                                0.,
0342                                (*R9Map)[candidateRef],
0343                                sMin,
0344                                sMaj,
0345                                seedId,
0346                                scRef->clustersSize(),
0347                                scRef->size(),
0348                                mEnergies,
0349                                mDetIdIds,
0350                                mTimes,
0351                                rechitZeroSuppression);  //read for(ieta){for(iphi){}}
0352     } else {                                            // Candidate is a scouting electron
0353       outElectrons->emplace_back(candidate.pt(),
0354                                  candidate.eta(),
0355                                  candidate.phi(),
0356                                  candidate.mass(),
0357                                  scRef->rawEnergy(),
0358                                  scRef->preshowerEnergy(),
0359                                  scRef->correctedEnergyUncertainty(),
0360                                  trkd0,
0361                                  trkdz,
0362                                  trkpt,
0363                                  trketa,
0364                                  trkphi,
0365                                  trkpMode,
0366                                  trketaMode,
0367                                  trkphiMode,
0368                                  trkqoverpModeError,
0369                                  trkchi2overndf,
0370                                  (*DetaMap)[candidateRef],
0371                                  (*DphiMap)[candidateRef],
0372                                  (*SigmaIEtaIEtaMap)[candidateRef],
0373                                  HoE,
0374                                  (*OneOEMinusOneOPMap)[candidateRef],
0375                                  (*MissingHitsMap)[candidateRef],
0376                                  trkcharge,
0377                                  (*fBremMap)[candidateRef],
0378                                  (*EcalPFClusterIsoMap)[candidateRef],
0379                                  (*HcalPFClusterIsoMap)[candidateRef],
0380                                  (*EleGsfTrackIsoMap)[candidateRef],
0381                                  (*R9Map)[candidateRef],
0382                                  sMin,
0383                                  sMaj,
0384                                  seedId,
0385                                  scRef->clustersSize(),
0386                                  scRef->size(),
0387                                  mEnergies,
0388                                  mDetIdIds,
0389                                  mTimes,
0390                                  rechitZeroSuppression);  //read for(ieta){for(iphi){}}
0391     }
0392   }
0393 
0394   // Put output
0395   iEvent.put(std::move(outElectrons));
0396   iEvent.put(std::move(outPhotons));
0397 }
0398 
0399 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0400 void HLTScoutingEgammaProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0401   edm::ParameterSetDescription desc;
0402   desc.add<edm::InputTag>("EgammaCandidates", edm::InputTag("hltEgammaCandidates"));
0403   desc.add<edm::InputTag>("EgammaGsfTracks", edm::InputTag("hltEgammaGsfTracks"));
0404   desc.add<edm::InputTag>("SigmaIEtaIEtaMap", edm::InputTag("hltEgammaClusterShape:sigmaIEtaIEta5x5"));
0405   desc.add<edm::InputTag>("r9Map", edm::InputTag("hltEgammaR9ID:r95x5"));
0406   desc.add<edm::InputTag>("HoverEMap", edm::InputTag("hltEgammaHoverE"));
0407   desc.add<edm::InputTag>("DetaMap", edm::InputTag("hltEgammaGsfTrackVarsUnseeded:DetaSeed"));
0408   desc.add<edm::InputTag>("DphiMap", edm::InputTag("hltEgammaGsfTrackVarsUnseeded:Dphi"));
0409   desc.add<edm::InputTag>("MissingHitsMap", edm::InputTag("hltEgammaGsfTrackVarsUnseeded:MissingHits"));
0410   desc.add<edm::InputTag>("OneOEMinusOneOPMap", edm::InputTag("hltEgammaGsfTrackVarsUnseeded:OneOESuperMinusOneOP"));
0411   desc.add<edm::InputTag>("fBremMap", edm::InputTag("hltEgammaGsfTrackVarsUnseeded:fbrem"));
0412   desc.add<edm::InputTag>("EcalPFClusterIsoMap", edm::InputTag("hltEgammaEcalPFClusterIso"));
0413   desc.add<edm::InputTag>("EleGsfTrackIsoMap", edm::InputTag("hltEgammaEleGsfTrackIso"));
0414   desc.add<edm::InputTag>("HcalPFClusterIsoMap", edm::InputTag("hltEgammaHcalPFClusterIso"));
0415   desc.add<double>("egammaPtCut", 4.0);
0416   desc.add<double>("egammaEtaCut", 2.5);
0417   desc.add<double>("egammaHoverECut", 1.0);
0418   desc.add<std::vector<double>>("egammaSigmaIEtaIEtaCut", {99999.0, 99999.0});
0419   desc.add<std::vector<double>>("absEtaBinUpperEdges", {1.479, 5.0});
0420   desc.add<bool>("saveRecHitTiming", false);
0421   desc.add<int>("mantissaPrecision", 10)->setComment("default float16, change to 23 for float32");
0422   desc.add<int>("rechitMatrixSize", 10);
0423   desc.add<bool>("rechitZeroSuppression", true);
0424   desc.add<edm::InputTag>("ecalRechitEB", edm::InputTag("hltEcalRecHit:EcalRecHitsEB"));
0425   desc.add<edm::InputTag>("ecalRechitEE", edm::InputTag("hltEcalRecHit:EcalRecHitsEE"));
0426   descriptions.add("hltScoutingEgammaProducer", desc);
0427 }
0428 
0429 // declare this class as a framework plugin
0430 #include "FWCore/Framework/interface/MakerMacros.h"
0431 DEFINE_FWK_MODULE(HLTScoutingEgammaProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!