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   struct GsfTrackExtrapolations {
0044     GsfTrackExtrapolations() {}
0045     void operator()(const reco::GsfTrack& trk,
0046                     const reco::SuperCluster& sc,
0047                     const MultiTrajectoryStateTransform& mtsTransform);
0048     TrajectoryStateOnSurface innTSOS;
0049     TrajectoryStateOnSurface outTSOS;
0050     TrajectoryStateOnSurface sclTSOS;
0051
0052     GlobalVector innMom, outMom;
0053     GlobalPoint sclPos;
0054   };
0055
0056 public:
0057   explicit EgammaHLTGsfTrackVarProducer(const edm::ParameterSet&);
0058   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0059   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0060   void fillAbsAbleVar(float& existVal, const float newVal) const {
0061     if (std::abs(newVal) < std::abs(existVal)) {
0062       existVal = produceAbsValues_ ? std::abs(newVal) : newVal;
0063     }
0064   }
0065
0066 private:
0067   const edm::EDGetTokenT<reco::RecoEcalCandidateCollection> recoEcalCandToken_;
0068   const edm::EDGetTokenT<reco::ElectronCollection> electronToken_;
0069   const edm::EDGetTokenT<reco::GsfTrackCollection> gsfTrackToken_;
0070   const edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
0071
0072   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0073   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeometryToken_;
0074
0075   const int upperTrackNrToRemoveCut_;
0076   const int lowerTrackNrToRemoveCut_;
0077   const bool useDefaultValuesForBarrel_;
0078   const bool useDefaultValuesForEndcap_;
0079   const bool produceAbsValues_;
0080
0081   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> dEtaMapPutToken_;
0082   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> dEtaSeedMapPutToken_;
0083   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> dPhiMapPutToken_;
0084   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> oneOverESuperMinusOneOverPMapPutToken_;
0085   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> oneOverESeedMinusOneOverPMapPutToken_;
0086   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> missingHitsMapPutToken_;
0087   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> validHitsMapPutToken_;
0088   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> nLayerITMapPutToken_;
0089   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> chi2MapPutToken_;
0090   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> fbremMapPutToken_;
0091 };
0092
0093 namespace {
0094
0095   float calRelDelta(float a, float b, float defaultVal = 0.f) { return a != 0.f ? (a - b) / a : defaultVal; }
0096
0097 }  // namespace
0098
0099 EgammaHLTGsfTrackVarProducer::EgammaHLTGsfTrackVarProducer(const edm::ParameterSet& config)
0100     : recoEcalCandToken_(
0101           consumes<reco::RecoEcalCandidateCollection>(config.getParameter<edm::InputTag>("recoEcalCandidateProducer"))),
0102       electronToken_{consumes<reco::ElectronCollection>(config.getParameter<edm::InputTag>("inputCollection"))},
0103       gsfTrackToken_{consumes<reco::GsfTrackCollection>(config.getParameter<edm::InputTag>("inputCollection"))},
0104       beamSpotToken_{consumes<reco::BeamSpot>(config.getParameter<edm::InputTag>("beamSpotProducer"))},
0105       magneticFieldToken_{esConsumes()},
0106       trackerGeometryToken_{esConsumes()},
0107       upperTrackNrToRemoveCut_{config.getParameter<int>("upperTrackNrToRemoveCut")},
0108       lowerTrackNrToRemoveCut_{config.getParameter<int>("lowerTrackNrToRemoveCut")},
0109       useDefaultValuesForBarrel_{config.getParameter<bool>("useDefaultValuesForBarrel")},
0110       useDefaultValuesForEndcap_{config.getParameter<bool>("useDefaultValuesForEndcap")},
0111       produceAbsValues_{config.getParameter<bool>("produceAbsValues")},
0112       dEtaMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("Deta").setBranchAlias("deta")},
0113       dEtaSeedMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("DetaSeed").setBranchAlias("detaseed")},
0114       dPhiMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("Dphi").setBranchAlias("dphi")},
0115       oneOverESuperMinusOneOverPMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("OneOESuperMinusOneOP")},
0116       oneOverESeedMinusOneOverPMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("OneOESeedMinusOneOP")},
0117       missingHitsMapPutToken_{
0118           produces<reco::RecoEcalCandidateIsolationMap>("MissingHits").setBranchAlias("missinghits")},
0119       validHitsMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("ValidHits").setBranchAlias("validhits")},
0120       nLayerITMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("NLayerIT").setBranchAlias("nlayerit")},
0121       chi2MapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("Chi2").setBranchAlias("chi2")},
0122       fbremMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("fbrem")} {}
0123
0124 void EgammaHLTGsfTrackVarProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0125   edm::ParameterSetDescription desc;
0126   desc.add<edm::InputTag>(("recoEcalCandidateProducer"), edm::InputTag("hltRecoEcalSuperClusterActivityCandidate"));
0127   desc.add<edm::InputTag>(("inputCollection"), edm::InputTag("hltActivityElectronGsfTracks"));
0128   desc.add<edm::InputTag>(("beamSpotProducer"), edm::InputTag("hltOnlineBeamSpot"));
0129   desc.add<int>(("upperTrackNrToRemoveCut"), 9999);
0130   desc.add<int>(("lowerTrackNrToRemoveCut"), -1);
0131   desc.add<bool>(("useDefaultValuesForBarrel"), false);
0132   desc.add<bool>(("useDefaultValuesForEndcap"), false);
0133   desc.add<bool>(("produceAbsValues"), false);
0134
0135   descriptions.add("hltEgammaHLTGsfTrackVarProducer", desc);
0136 }
0137 void EgammaHLTGsfTrackVarProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0138   // Get the HLT filtered objects
0139   auto recoEcalCandHandle = iEvent.getHandle(recoEcalCandToken_);
0140
0141   auto const& beamSpotPosition = iEvent.get(beamSpotToken_).position();
0142   auto const& magneticField = iSetup.getData(magneticFieldToken_);
0143   auto const& trackerGeometry = iSetup.getData(trackerGeometryToken_);
0144
0145   TransverseImpactPointExtrapolator extrapolator{
0146       GsfPropagatorAdapter{AnalyticalPropagator(&magneticField, anyDirection)}};
0147
0148   reco::RecoEcalCandidateIsolationMap dEtaMap(recoEcalCandHandle);
0149   reco::RecoEcalCandidateIsolationMap dEtaSeedMap(recoEcalCandHandle);
0150   reco::RecoEcalCandidateIsolationMap dPhiMap(recoEcalCandHandle);
0151   reco::RecoEcalCandidateIsolationMap oneOverESuperMinusOneOverPMap(recoEcalCandHandle);
0152   reco::RecoEcalCandidateIsolationMap oneOverESeedMinusOneOverPMap(recoEcalCandHandle);
0153   reco::RecoEcalCandidateIsolationMap missingHitsMap(recoEcalCandHandle);
0154   reco::RecoEcalCandidateIsolationMap validHitsMap(recoEcalCandHandle);
0155   reco::RecoEcalCandidateIsolationMap nLayerITMap(recoEcalCandHandle);
0156   reco::RecoEcalCandidateIsolationMap chi2Map(recoEcalCandHandle);
0157   reco::RecoEcalCandidateIsolationMap fbremMap(recoEcalCandHandle);
0158
0159   for (unsigned int iRecoEcalCand = 0; iRecoEcalCand < recoEcalCandHandle->size(); ++iRecoEcalCand) {
0160     reco::RecoEcalCandidateRef recoEcalCandRef(recoEcalCandHandle, iRecoEcalCand);
0161
0162     const reco::SuperClusterRef scRef = recoEcalCandRef->superCluster();
0163     // the idea is that we can take the tracks from properly associated
0164     // electrons or just take all gsf tracks with that sc as a seed
0165     std::vector<const reco::GsfTrack*> gsfTracks;
0166     if (auto electronHandle = iEvent.getHandle(electronToken_)) {
0167       for (auto const& ele : *electronHandle) {
0168         if (ele.superCluster() == scRef) {
0169           gsfTracks.push_back(ele.gsfTrack().get());
0170         }
0171       }
0172     } else {
0173       for (auto const& trk : iEvent.get(gsfTrackToken_)) {
0174         auto elseed = trk.extra()->seedRef().castTo<reco::ElectronSeedRef>();
0175         if (elseed->caloCluster().castTo<reco::SuperClusterRef>() == scRef) {
0176           gsfTracks.push_back(&trk);
0177         }
0178       }
0179     }
0180
0181     int nLayerITValue = -1;
0182     int validHitsValue = 0;
0183     float chi2Value = 9999999.;
0184     float missingHitsValue = 9999999;
0185     float dEtaInValue = 999999;
0186     float dEtaSeedInValue = 999999;
0187     float dPhiInValue = 999999;
0188     float oneOverESuperMinusOneOverPValue = 999999;
0189     float oneOverESeedMinusOneOverPValue = 999999;
0190     float fbrem = 999999;
0191
0192     const int nrTracks = gsfTracks.size();
0193     const bool rmCutsDueToNrTracks = nrTracks <= lowerTrackNrToRemoveCut_ || nrTracks >= upperTrackNrToRemoveCut_;
0194     // to use the default values, we require at least one track...
0195     const bool useDefaultValues = std::abs(recoEcalCandRef->eta()) < 1.479
0196                                       ? useDefaultValuesForBarrel_ && nrTracks >= 1
0197                                       : useDefaultValuesForEndcap_ && nrTracks >= 1;
0198
0199     MultiTrajectoryStateTransform mtsTransform(&trackerGeometry, &magneticField);
0200     GsfTrackExtrapolations gsfTrackExtrapolations;
0201
0202     if (rmCutsDueToNrTracks || useDefaultValues) {
0203       nLayerITValue = 100;
0204       dEtaInValue = 0;
0205       dEtaSeedInValue = 0;
0206       dPhiInValue = 0;
0207       missingHitsValue = 0;
0208       validHitsValue = 100;
0209       chi2Value = 0;
0210       oneOverESuperMinusOneOverPValue = 0;
0211       oneOverESeedMinusOneOverPValue = 0;
0212       fbrem = 0;
0213     } else {
0214       for (size_t trkNr = 0; trkNr < gsfTracks.size(); trkNr++) {
0215         GlobalPoint scPos(scRef->x(), scRef->y(), scRef->z());
0216
0217         gsfTrackExtrapolations(*gsfTracks[trkNr], *scRef, mtsTransform);
0218
0219         EleRelPointPair scAtVtx(scRef->position(), gsfTrackExtrapolations.sclPos, beamSpotPosition);
0220
0221         fbrem = calRelDelta(gsfTrackExtrapolations.innMom.mag(), gsfTrackExtrapolations.outMom.mag(), fbrem);
0222
0223         float trkP = gsfTracks[trkNr]->p();
0224         if (scRef->energy() != 0 && trkP != 0) {
0225           fillAbsAbleVar(oneOverESuperMinusOneOverPValue, 1 / scRef->energy() - 1 / trkP);
0226         }
0227         if (scRef->seed().isNonnull() && scRef->seed()->energy() != 0 && trkP != 0) {
0228           fillAbsAbleVar(oneOverESeedMinusOneOverPValue, 1 / scRef->seed()->energy() - 1 / trkP);
0229         }
0230
0231         if (gsfTracks[trkNr]->missingInnerHits() < missingHitsValue) {
0232           missingHitsValue = gsfTracks[trkNr]->missingInnerHits();
0233         }
0234
0235         // we are saving the best value, and highest value of validHits is the
0236         // best
0237         if (gsfTracks[trkNr]->numberOfValidHits() > validHitsValue) {
0238           validHitsValue = gsfTracks[trkNr]->numberOfValidHits();
0239         }
0240
0241         if (gsfTracks[trkNr]->hitPattern().pixelLayersWithMeasurement() > nLayerITValue) {
0242           nLayerITValue = gsfTracks[trkNr]->hitPattern().pixelLayersWithMeasurement();
0243         }
0244
0245         if (gsfTracks[trkNr]->normalizedChi2() < chi2Value) {
0246           chi2Value = gsfTracks[trkNr]->normalizedChi2();
0247         }
0248
0249         fillAbsAbleVar(dEtaInValue, scAtVtx.dEta());
0250         fillAbsAbleVar(dEtaSeedInValue, scAtVtx.dEta() - scRef->position().eta() + scRef->seed()->position().eta());
0251         fillAbsAbleVar(dPhiInValue, scAtVtx.dPhi());
0252       }
0253     }
0254
0255     dEtaMap.insert(recoEcalCandRef, dEtaInValue);
0256     dEtaSeedMap.insert(recoEcalCandRef, dEtaSeedInValue);
0257     dPhiMap.insert(recoEcalCandRef, dPhiInValue);
0258     oneOverESuperMinusOneOverPMap.insert(recoEcalCandRef, oneOverESuperMinusOneOverPValue);
0259     oneOverESeedMinusOneOverPMap.insert(recoEcalCandRef, oneOverESeedMinusOneOverPValue);
0260     missingHitsMap.insert(recoEcalCandRef, missingHitsValue);
0261     validHitsMap.insert(recoEcalCandRef, validHitsValue);
0262     nLayerITMap.insert(recoEcalCandRef, nLayerITValue);
0263     chi2Map.insert(recoEcalCandRef, chi2Value);
0264     fbremMap.insert(recoEcalCandRef, fbrem);
0265   }
0266
0267   iEvent.emplace(dEtaMapPutToken_, dEtaMap);
0268   iEvent.emplace(dEtaSeedMapPutToken_, dEtaSeedMap);
0269   iEvent.emplace(dPhiMapPutToken_, dPhiMap);
0270   iEvent.emplace(oneOverESuperMinusOneOverPMapPutToken_, oneOverESuperMinusOneOverPMap);
0271   iEvent.emplace(oneOverESeedMinusOneOverPMapPutToken_, oneOverESeedMinusOneOverPMap);
0272   iEvent.emplace(missingHitsMapPutToken_, missingHitsMap);
0273   iEvent.emplace(validHitsMapPutToken_, validHitsMap);
0274   iEvent.emplace(nLayerITMapPutToken_, nLayerITMap);
0275   iEvent.emplace(chi2MapPutToken_, chi2Map);
0276   iEvent.emplace(fbremMapPutToken_, fbremMap);
0277 }
0278
0279 void EgammaHLTGsfTrackVarProducer::GsfTrackExtrapolations::operator()(
0280     const reco::GsfTrack& trk, const reco::SuperCluster& sc, const MultiTrajectoryStateTransform& mtsTransform) {
0281   innTSOS = mtsTransform.innerStateOnSurface(trk);
0282   outTSOS = mtsTransform.outerStateOnSurface(trk);
0283   sclTSOS = mtsTransform.extrapolatedState(innTSOS, GlobalPoint(sc.x(), sc.y(), sc.z()));
0284
0285   multiTrajectoryStateMode::momentumFromModeCartesian(innTSOS, innMom);
0286   multiTrajectoryStateMode::positionFromModeCartesian(sclTSOS, sclPos);
0287   multiTrajectoryStateMode::momentumFromModeCartesian(outTSOS, outMom);
0288 }
0289
0290 #include "FWCore/Framework/interface/MakerMacros.h"
0291 DEFINE_FWK_MODULE(EgammaHLTGsfTrackVarProducer);