EgammaHLTProducers/plugins/EgammaHLTGsfTrackVarProducer.cc

0001 /** \class EgammaHLTGsfTrackVarProducer
0002  *
0003  *  \author Roberto Covarelli (CERN)
0004  *
0005  * $Id: EgammaHLTGsfTrackVarProducer.cc,v 1.1 2012/01/23 12:56:38 sharper Exp $
0006  *
0007  */
0008
0009 // this class is designed to calculate dEtaIn,dPhiIn gsf track - supercluster
0010 // pairs it can take as input std::vector<Electron> which the gsf track-sc is
0011 // already done or it can run over the std::vector<GsfTrack> directly in which
0012 // case it will pick the smallest dEta,dPhi the dEta, dPhi do not have to be
0013 // from the same track it can optionally set dEta, dPhi to 0 based on the number
0014 // of tracks found
0015
0016 #include "DataFormats/EgammaCandidates/interface/Electron.h"
0017 #include "DataFormats/EgammaCandidates/interface/ElectronIsolationAssociation.h"
0018 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0019 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0020 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0021 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
0022 #include "DataFormats/TrackReco/interface/Track.h"
0023 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026 #include "FWCore/Framework/interface/global/EDProducer.h"
0027 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0029 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0030 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0031 #include "MagneticField/Engine/interface/MagneticField.h"
0032 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0033 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h"
0034 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0035 #include "TrackingTools/GsfTools/interface/GsfPropagatorAdapter.h"
0036 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateMode.h"
0037 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateTransform.h"
0038 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0039 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0040
0041 class EgammaHLTGsfTrackVarProducer : public edm::global::EDProducer<> {
0042 public:
0043   explicit EgammaHLTGsfTrackVarProducer(const edm::ParameterSet&);
0044   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0045   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0046
0047 private:
0048   const edm::EDGetTokenT<reco::RecoEcalCandidateCollection> recoEcalCandToken_;
0049   const edm::EDGetTokenT<reco::ElectronCollection> electronToken_;
0050   const edm::EDGetTokenT<reco::GsfTrackCollection> gsfTrackToken_;
0051   const edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
0052
0053   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0054   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeometryToken_;
0055
0056   const int upperTrackNrToRemoveCut_;
0057   const int lowerTrackNrToRemoveCut_;
0058   const bool useDefaultValuesForBarrel_;
0059   const bool useDefaultValuesForEndcap_;
0060
0061   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> dEtaMapPutToken_;
0062   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> dEtaSeedMapPutToken_;
0063   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> dPhiMapPutToken_;
0064   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> oneOverESuperMinusOneOverPMapPutToken_;
0065   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> oneOverESeedMinusOneOverPMapPutToken_;
0066   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> missingHitsMapPutToken_;
0067   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> validHitsMapPutToken_;
0068   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> nLayerITMapPutToken_;
0069   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> chi2MapPutToken_;
0070 };
0071
0072 EgammaHLTGsfTrackVarProducer::EgammaHLTGsfTrackVarProducer(const edm::ParameterSet& config)
0073     : recoEcalCandToken_(
0074           consumes<reco::RecoEcalCandidateCollection>(config.getParameter<edm::InputTag>("recoEcalCandidateProducer"))),
0075       electronToken_{consumes<reco::ElectronCollection>(config.getParameter<edm::InputTag>("inputCollection"))},
0076       gsfTrackToken_{consumes<reco::GsfTrackCollection>(config.getParameter<edm::InputTag>("inputCollection"))},
0077       beamSpotToken_{consumes<reco::BeamSpot>(config.getParameter<edm::InputTag>("beamSpotProducer"))},
0078       magneticFieldToken_{esConsumes()},
0079       trackerGeometryToken_{esConsumes()},
0080       upperTrackNrToRemoveCut_{config.getParameter<int>("upperTrackNrToRemoveCut")},
0081       lowerTrackNrToRemoveCut_{config.getParameter<int>("lowerTrackNrToRemoveCut")},
0082       useDefaultValuesForBarrel_{config.getParameter<bool>("useDefaultValuesForBarrel")},
0083       useDefaultValuesForEndcap_{config.getParameter<bool>("useDefaultValuesForEndcap")},
0084       dEtaMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("Deta").setBranchAlias("deta")},
0085       dEtaSeedMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("DetaSeed").setBranchAlias("detaseed")},
0086       dPhiMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("Dphi").setBranchAlias("dphi")},
0087       oneOverESuperMinusOneOverPMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("OneOESuperMinusOneOP")},
0088       oneOverESeedMinusOneOverPMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("OneOESeedMinusOneOP")},
0089       missingHitsMapPutToken_{
0090           produces<reco::RecoEcalCandidateIsolationMap>("MissingHits").setBranchAlias("missinghits")},
0091       validHitsMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("ValidHits").setBranchAlias("validhits")},
0092       nLayerITMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("NLayerIT").setBranchAlias("nlayerit")},
0093       chi2MapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("Chi2").setBranchAlias("chi2")} {}
0094
0095 void EgammaHLTGsfTrackVarProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0096   edm::ParameterSetDescription desc;
0097   desc.add<edm::InputTag>(("recoEcalCandidateProducer"), edm::InputTag("hltRecoEcalSuperClusterActivityCandidate"));
0098   desc.add<edm::InputTag>(("inputCollection"), edm::InputTag("hltActivityElectronGsfTracks"));
0099   desc.add<edm::InputTag>(("beamSpotProducer"), edm::InputTag("hltOnlineBeamSpot"));
0100   desc.add<int>(("upperTrackNrToRemoveCut"), 9999);
0101   desc.add<int>(("lowerTrackNrToRemoveCut"), -1);
0102   desc.add<bool>(("useDefaultValuesForBarrel"), false);
0103   desc.add<bool>(("useDefaultValuesForEndcap"), false);
0104
0105   descriptions.add("hltEgammaHLTGsfTrackVarProducer", desc);
0106 }
0107 void EgammaHLTGsfTrackVarProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0108   // Get the HLT filtered objects
0109   auto recoEcalCandHandle = iEvent.getHandle(recoEcalCandToken_);
0110
0111   auto const& beamSpotPosition = iEvent.get(beamSpotToken_).position();
0112   auto const& magneticField = iSetup.getData(magneticFieldToken_);
0113   auto const& trackerGeometry = iSetup.getData(trackerGeometryToken_);
0114
0115   TransverseImpactPointExtrapolator extrapolator{
0116       GsfPropagatorAdapter{AnalyticalPropagator(&magneticField, anyDirection)}};
0117
0118   reco::RecoEcalCandidateIsolationMap dEtaMap(recoEcalCandHandle);
0119   reco::RecoEcalCandidateIsolationMap dEtaSeedMap(recoEcalCandHandle);
0120   reco::RecoEcalCandidateIsolationMap dPhiMap(recoEcalCandHandle);
0121   reco::RecoEcalCandidateIsolationMap oneOverESuperMinusOneOverPMap(recoEcalCandHandle);
0122   reco::RecoEcalCandidateIsolationMap oneOverESeedMinusOneOverPMap(recoEcalCandHandle);
0123   reco::RecoEcalCandidateIsolationMap missingHitsMap(recoEcalCandHandle);
0124   reco::RecoEcalCandidateIsolationMap validHitsMap(recoEcalCandHandle);
0125   reco::RecoEcalCandidateIsolationMap nLayerITMap(recoEcalCandHandle);
0126   reco::RecoEcalCandidateIsolationMap chi2Map(recoEcalCandHandle);
0127
0128   for (unsigned int iRecoEcalCand = 0; iRecoEcalCand < recoEcalCandHandle->size(); ++iRecoEcalCand) {
0129     reco::RecoEcalCandidateRef recoEcalCandRef(recoEcalCandHandle, iRecoEcalCand);
0130
0131     const reco::SuperClusterRef scRef = recoEcalCandRef->superCluster();
0132     // the idea is that we can take the tracks from properly associated
0133     // electrons or just take all gsf tracks with that sc as a seed
0134     std::vector<const reco::GsfTrack*> gsfTracks;
0135     if (auto electronHandle = iEvent.getHandle(electronToken_)) {
0136       for (auto const& ele : *electronHandle) {
0137         if (ele.superCluster() == scRef) {
0138           gsfTracks.push_back(ele.gsfTrack().get());
0139         }
0140       }
0141     } else {
0142       for (auto const& trk : iEvent.get(gsfTrackToken_)) {
0143         auto elseed = trk.extra()->seedRef().castTo<reco::ElectronSeedRef>();
0144         if (elseed->caloCluster().castTo<reco::SuperClusterRef>() == scRef) {
0145           gsfTracks.push_back(&trk);
0146         }
0147       }
0148     }
0149
0150     int nLayerITValue = -1;
0151     int validHitsValue = 0;
0152     float chi2Value = 9999999.;
0153     float missingHitsValue = 9999999;
0154     float dEtaInValue = 999999;
0155     float dEtaSeedInValue = 999999;
0156     float dPhiInValue = 999999;
0157     float oneOverESuperMinusOneOverPValue = 999999;
0158     float oneOverESeedMinusOneOverPValue = 999999;
0159
0160     const int nrTracks = gsfTracks.size();
0161     const bool rmCutsDueToNrTracks = nrTracks <= lowerTrackNrToRemoveCut_ || nrTracks >= upperTrackNrToRemoveCut_;
0162     // to use the default values, we require at least one track...
0163     const bool useDefaultValues = std::abs(recoEcalCandRef->eta()) < 1.479
0164                                       ? useDefaultValuesForBarrel_ && nrTracks >= 1
0165                                       : useDefaultValuesForEndcap_ && nrTracks >= 1;
0166
0167     if (rmCutsDueToNrTracks || useDefaultValues) {
0168       nLayerITValue = 100;
0169       dEtaInValue = 0;
0170       dEtaSeedInValue = 0;
0171       dPhiInValue = 0;
0172       missingHitsValue = 0;
0173       validHitsValue = 100;
0174       chi2Value = 0;
0175       oneOverESuperMinusOneOverPValue = 0;
0176       oneOverESeedMinusOneOverPValue = 0;
0177     } else {
0178       for (size_t trkNr = 0; trkNr < gsfTracks.size(); trkNr++) {
0179         GlobalPoint scPos(scRef->x(), scRef->y(), scRef->z());
0180
0181         GlobalPoint trackExtrapToSC;
0182         {
0183           auto innTSOS =
0184               MultiTrajectoryStateTransform::innerStateOnSurface(*gsfTracks[trkNr], trackerGeometry, &magneticField);
0185           auto posTSOS = extrapolator.extrapolate(innTSOS, scPos);
0186           multiTrajectoryStateMode::positionFromModeCartesian(posTSOS, trackExtrapToSC);
0187         }
0188
0189         EleRelPointPair scAtVtx(scRef->position(), trackExtrapToSC, beamSpotPosition);
0190
0191         float trkP = gsfTracks[trkNr]->p();
0192         if (scRef->energy() != 0 && trkP != 0) {
0193           if (std::abs(1 / scRef->energy() - 1 / trkP) < oneOverESuperMinusOneOverPValue) {
0194             oneOverESuperMinusOneOverPValue = std::abs(1 / scRef->energy() - 1 / trkP);
0195           }
0196         }
0197         if (scRef->seed().isNonnull() && scRef->seed()->energy() != 0 && trkP != 0) {
0198           if (std::abs(1 / scRef->seed()->energy() - 1 / trkP) < oneOverESeedMinusOneOverPValue) {
0199             oneOverESeedMinusOneOverPValue = std::abs(1 / scRef->seed()->energy() - 1 / trkP);
0200           }
0201         }
0202
0203         if (gsfTracks[trkNr]->missingInnerHits() < missingHitsValue) {
0204           missingHitsValue = gsfTracks[trkNr]->missingInnerHits();
0205         }
0206
0207         // we are saving the best value, and highest value of validHits is the
0208         // best
0209         if (gsfTracks[trkNr]->numberOfValidHits() > validHitsValue) {
0210           validHitsValue = gsfTracks[trkNr]->numberOfValidHits();
0211         }
0212
0213         if (gsfTracks[trkNr]->hitPattern().pixelLayersWithMeasurement() > nLayerITValue) {
0214           nLayerITValue = gsfTracks[trkNr]->hitPattern().pixelLayersWithMeasurement();
0215         }
0216
0217         if (gsfTracks[trkNr]->normalizedChi2() < chi2Value) {
0218           chi2Value = gsfTracks[trkNr]->normalizedChi2();
0219         }
0220
0221         if (std::abs(scAtVtx.dEta()) < dEtaInValue) {
0222           // we are allowing them to come from different tracks
0223           dEtaInValue = std::abs(scAtVtx.dEta());
0224         }
0225
0226         if (std::abs(scAtVtx.dEta()) < dEtaSeedInValue) {
0227           dEtaSeedInValue = std::abs(scAtVtx.dEta() - scRef->position().eta() + scRef->seed()->position().eta());
0228         }
0229
0230         if (std::abs(scAtVtx.dPhi()) < dPhiInValue) {
0231           // we are allowing them to come from different tracks
0232           dPhiInValue = std::abs(scAtVtx.dPhi());
0233         }
0234       }
0235     }
0236
0237     dEtaMap.insert(recoEcalCandRef, dEtaInValue);
0238     dEtaSeedMap.insert(recoEcalCandRef, dEtaSeedInValue);
0239     dPhiMap.insert(recoEcalCandRef, dPhiInValue);
0240     oneOverESuperMinusOneOverPMap.insert(recoEcalCandRef, oneOverESuperMinusOneOverPValue);
0241     oneOverESeedMinusOneOverPMap.insert(recoEcalCandRef, oneOverESeedMinusOneOverPValue);
0242     missingHitsMap.insert(recoEcalCandRef, missingHitsValue);
0243     validHitsMap.insert(recoEcalCandRef, validHitsValue);
0244     nLayerITMap.insert(recoEcalCandRef, nLayerITValue);
0245     chi2Map.insert(recoEcalCandRef, chi2Value);
0246   }
0247
0248   iEvent.emplace(dEtaMapPutToken_, dEtaMap);
0249   iEvent.emplace(dEtaSeedMapPutToken_, dEtaSeedMap);
0250   iEvent.emplace(dPhiMapPutToken_, dPhiMap);
0251   iEvent.emplace(oneOverESuperMinusOneOverPMapPutToken_, oneOverESuperMinusOneOverPMap);
0252   iEvent.emplace(oneOverESeedMinusOneOverPMapPutToken_, oneOverESeedMinusOneOverPMap);
0253   iEvent.emplace(missingHitsMapPutToken_, missingHitsMap);
0254   iEvent.emplace(validHitsMapPutToken_, validHitsMap);
0255   iEvent.emplace(nLayerITMapPutToken_, nLayerITMap);
0256   iEvent.emplace(chi2MapPutToken_, chi2Map);
0257 }
0258
0259 #include "FWCore/Framework/interface/MakerMacros.h"
0260 DEFINE_FWK_MODULE(EgammaHLTGsfTrackVarProducer);