EgammaHLTProducers/plugins/EgammaHLTPixelMatchElectronAlgo.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    EgammaHLTAlgos
0004 // Class:      EgammaHLTPixelMatchElectronAlgo.
0005 //
0006 /**\class EgammaHLTPixelMatchElectronAlgo EgammaHLTAlgos/EgammaHLTPixelMatchElectronAlgo
0007
0008  Description: top algorithm producing TrackCandidate and Electron objects from supercluster
0009               driven pixel seeded Ckf tracking for HLT
0010 */
0011 //
0012 // Original Author:  Monica Vazquez Acosta (CERN)
0013 //
0014 //
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017
0018 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0019 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0020 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0021 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateClosestToPoint.h"
0022
0023 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0024 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0025 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0026 #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
0027
0028 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
0029 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0030
0031 #include "DataFormats/TrackReco/interface/Track.h"
0032 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0033
0034 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0035 #include "DataFormats/Math/interface/Point3D.h"
0036
0037 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateMode.h"
0038
0039 #include "EgammaHLTPixelMatchElectronAlgo.h"
0040
0041 using namespace edm;
0042 using namespace std;
0043 using namespace reco;
0044
0045 EgammaHLTPixelMatchElectronAlgo::EgammaHLTPixelMatchElectronAlgo(const edm::ParameterSet& conf,
0046                                                                  edm::ConsumesCollector&& iC)
0047     : trackProducer_(iC.consumes(conf.getParameter<edm::InputTag>("TrackProducer"))),
0048       gsfTrackProducer_(iC.consumes(conf.getParameter<edm::InputTag>("GsfTrackProducer"))),
0049       useGsfTracks_(conf.getParameter<bool>("UseGsfTracks")),
0050       bsProducer_(iC.consumes<reco::BeamSpot>(conf.getParameter<edm::InputTag>("BSProducer"))),
0051       magneticFieldToken_(iC.esConsumes()),
0052       trackerGeometryToken_(iC.esConsumes()) {}
0053
0054 void EgammaHLTPixelMatchElectronAlgo::setupES(const edm::EventSetup& eventSetup) {
0055   //services
0056   bool updateField = magneticFieldWatcher_.check(eventSetup);
0057   bool updateGeometry = trackerGeometryWatcher_.check(eventSetup);
0058
0059   if (updateField) {
0060     magField_ = eventSetup.getHandle(magneticFieldToken_);
0061   }
0062
0063   if (useGsfTracks_) {  //only need the geom and mtsTransform if we are doing gsf tracks
0064     if (updateGeometry) {
0065       trackerGeom_ = eventSetup.getHandle(trackerGeometryToken_);
0066     }
0067     if (updateField || updateGeometry || !mtsTransform_) {
0068       mtsTransform_ = std::make_unique<MultiTrajectoryStateTransform>(trackerGeom_.product(), magField_.product());
0069     }
0070   }
0071 }
0072
0073 void EgammaHLTPixelMatchElectronAlgo::run(Event& e, ElectronCollection& outEle) {
0074   // get the input
0075   edm::Handle<TrackCollection> tracksH;
0076   if (!useGsfTracks_)
0077     e.getByToken(trackProducer_, tracksH);
0078
0079   // get the input
0080   edm::Handle<GsfTrackCollection> gsfTracksH;
0081   if (useGsfTracks_)
0082     e.getByToken(gsfTrackProducer_, gsfTracksH);
0083
0084   //Get the Beam Spot position
0085   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0086   e.getByToken(bsProducer_, recoBeamSpotHandle);
0087
0088   // gets its position
0089   const BeamSpot::Point& bsPosition = recoBeamSpotHandle->position();
0090   Global3DPoint bs(bsPosition.x(), bsPosition.y(), 0);
0091   process(tracksH, gsfTracksH, outEle, bs);
0092
0093   return;
0094 }
0095
0096 void EgammaHLTPixelMatchElectronAlgo::process(edm::Handle<TrackCollection> tracksH,
0097                                               edm::Handle<GsfTrackCollection> gsfTracksH,
0098                                               ElectronCollection& outEle,
0099                                               Global3DPoint& bs) {
0100   if (!useGsfTracks_) {
0101     for (unsigned int i = 0; i < tracksH->size(); ++i) {
0102       const TrackRef trackRef = edm::Ref<TrackCollection>(tracksH, i);
0103       edm::RefToBase<TrajectorySeed> seed = trackRef->extra()->seedRef();
0104       ElectronSeedRef elseed = seed.castTo<ElectronSeedRef>();
0105
0106       edm::RefToBase<CaloCluster> caloCluster = elseed->caloCluster();
0107       SuperClusterRef scRef = caloCluster.castTo<SuperClusterRef>();
0108
0109       // Get the momentum at vertex (not at the innermost layer)
0110       TSCPBuilderNoMaterial tscpBuilder;
0111
0112       FreeTrajectoryState fts = trajectoryStateTransform::innerFreeState(*trackRef, magField_.product());
0113       TrajectoryStateClosestToPoint tscp = tscpBuilder(fts, bs);
0114
0115       float scale = scRef->energy() / tscp.momentum().mag();
0116
0117       const math::XYZTLorentzVector momentum(
0118           tscp.momentum().x() * scale, tscp.momentum().y() * scale, tscp.momentum().z() * scale, scRef->energy());
0119
0120       Electron ele(trackRef->charge(), momentum, trackRef->vertex());
0121       ele.setSuperCluster(scRef);
0122       edm::Ref<TrackCollection> myRef(tracksH, i);
0123       ele.setTrack(myRef);
0124       outEle.push_back(ele);
0125     }  // loop over tracks
0126   } else {
0127     // clean gsf tracks
0128     std::vector<unsigned int> flag(gsfTracksH->size(), 0);
0129     if (gsfTracksH->empty())
0130       return;
0131
0132     for (unsigned int i = 0; i < gsfTracksH->size() - 1; ++i) {
0133       const GsfTrackRef trackRef1 = edm::Ref<GsfTrackCollection>(gsfTracksH, i);
0134       ElectronSeedRef elseed1 = trackRef1->extra()->seedRef().castTo<ElectronSeedRef>();
0135       SuperClusterRef scRef1 = elseed1->caloCluster().castTo<SuperClusterRef>();
0136
0137       TrajectoryStateOnSurface inTSOS = mtsTransform_->innerStateOnSurface((*trackRef1));
0138       TrajectoryStateOnSurface fts = mtsTransform_->extrapolatedState(inTSOS, bs);
0139       GlobalVector innMom;
0140       float pin1 = trackRef1->pMode();
0141       if (fts.isValid()) {
0142         multiTrajectoryStateMode::momentumFromModeCartesian(fts, innMom);
0143         pin1 = innMom.mag();
0144       }
0145
0146       for (unsigned int j = i + 1; j < gsfTracksH->size(); ++j) {
0147         const GsfTrackRef trackRef2 = edm::Ref<GsfTrackCollection>(gsfTracksH, j);
0148         ElectronSeedRef elseed2 = trackRef2->extra()->seedRef().castTo<ElectronSeedRef>();
0149         SuperClusterRef scRef2 = elseed2->caloCluster().castTo<SuperClusterRef>();
0150
0151         TrajectoryStateOnSurface inTSOS = mtsTransform_->innerStateOnSurface((*trackRef2));
0152         TrajectoryStateOnSurface fts = mtsTransform_->extrapolatedState(inTSOS, bs);
0153         GlobalVector innMom;
0154         float pin2 = trackRef2->pMode();
0155         if (fts.isValid()) {
0156           multiTrajectoryStateMode::momentumFromModeCartesian(fts, innMom);
0157           pin2 = innMom.mag();
0158         }
0159
0160         if (scRef1 == scRef2) {
0161           bool isSameLayer = false;
0162           bool iGsfInnermostWithLostHits = isInnerMostWithLostHits(trackRef2, trackRef1, isSameLayer);
0163
0164           if (iGsfInnermostWithLostHits) {
0165             flag[j] = 1;
0166           } else if (isSameLayer) {
0167             if (fabs((scRef1->energy() / pin1) - 1) < fabs((scRef2->energy() / pin2) - 1))
0168               flag[j] = 1;
0169           } else {
0170             flag[i] = 1;
0171           }
0172         }
0173       }
0174     }
0175
0176     for (unsigned int i = 0; i < gsfTracksH->size(); ++i) {
0177       if (flag[i] == 1)
0178         continue;
0179
0180       const GsfTrackRef trackRef = edm::Ref<GsfTrackCollection>(gsfTracksH, i);
0181       ElectronSeedRef elseed = trackRef->extra()->seedRef().castTo<ElectronSeedRef>();
0182       SuperClusterRef scRef = elseed->caloCluster().castTo<SuperClusterRef>();
0183
0184       // Get the momentum at vertex (not at the innermost layer)
0185       TrajectoryStateOnSurface inTSOS = mtsTransform_->innerStateOnSurface((*trackRef));
0186       TrajectoryStateOnSurface fts = mtsTransform_->extrapolatedState(inTSOS, bs);
0187       GlobalVector innMom;
0188       multiTrajectoryStateMode::momentumFromModeCartesian(inTSOS, innMom);
0189       if (fts.isValid()) {
0190         multiTrajectoryStateMode::momentumFromModeCartesian(fts, innMom);
0191       }
0192
0193       float scale = scRef->energy() / innMom.mag();
0194       const math::XYZTLorentzVector momentum(
0195           innMom.x() * scale, innMom.y() * scale, innMom.z() * scale, scRef->energy());
0196
0197       Electron ele(trackRef->charge(), momentum, trackRef->vertex());
0198       ele.setSuperCluster(scRef);
0199       edm::Ref<GsfTrackCollection> myRef(gsfTracksH, i);
0200       ele.setGsfTrack(myRef);
0201       outEle.push_back(ele);
0202     }
0203   }
0204 }
0205
0206 bool EgammaHLTPixelMatchElectronAlgo::isInnerMostWithLostHits(const reco::GsfTrackRef& nGsfTrack,
0207                                                               const reco::GsfTrackRef& iGsfTrack,
0208                                                               bool& sameLayer) {
0209   // define closest using the lost hits on the expectedhitsineer
0210   auto nLostHits = nGsfTrack->missingInnerHits();
0211   auto iLostHits = iGsfTrack->missingInnerHits();
0212
0213   if (nLostHits != iLostHits) {
0214     return (nLostHits > iLostHits);
0215   } else {
0216     sameLayer = true;
0217     return false;
0218   }
0219 }