Go to the documentation of this file.00001
00009 #include "RecoEgamma/EgammaHLTProducers/interface/EgammaHLTGsfTrackVarProducer.h"
00010
00011
00012 #include "FWCore/Framework/interface/Event.h"
00013 #include "FWCore/Framework/interface/EventSetup.h"
00014 #include "DataFormats/Common/interface/Handle.h"
00015 #include "FWCore/Framework/interface/ESHandle.h"
00016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00017 #include "FWCore/Utilities/interface/Exception.h"
00018
00019 #include "DataFormats/EgammaCandidates/interface/Electron.h"
00020 #include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
00021
00022 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
00023 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
00024 #include "DataFormats/Common/interface/RefToBase.h"
00025 #include "DataFormats/Common/interface/Ref.h"
00026 #include "DataFormats/Common/interface/RefProd.h"
00027
00028 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00029 #include "DataFormats/TrackReco/interface/Track.h"
00030 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
00031
00032 #include "DataFormats/TrackReco/interface/Track.h"
00033 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00034 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
00035 #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
00036
00037
00038
00039 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00040 #include "DataFormats/Math/interface/Point3D.h"
00041
00042 #include "RecoEgamma/EgammaTools/interface/ECALPositionCalculator.h"
00043
00044 #include "FWCore/Framework/interface/EventSetup.h"
00045 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00046 #include "MagneticField/Engine/interface/MagneticField.h"
00047 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
00048 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00049 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00050
00051 EgammaHLTGsfTrackVarProducer::EgammaHLTGsfTrackVarProducer(const edm::ParameterSet& config)
00052 {
00053 recoEcalCandTag_ = config.getParameter<edm::InputTag>("recoEcalCandidateProducer");
00054 inputCollectionTag_ = config.getParameter<edm::InputTag>("inputCollection");
00055 beamSpotTag_ = config.getParameter<edm::InputTag>("beamSpotProducer");
00056 upperTrackNrToRemoveCut_ = config.getParameter<int>("upperTrackNrToRemoveCut");
00057 lowerTrackNrToRemoveCut_ = config.getParameter<int>("lowerTrackNrToRemoveCut");
00058
00059
00060
00061
00062 produces < reco::RecoEcalCandidateIsolationMap >( "Deta" ).setBranchAlias( "deta" );
00063 produces < reco::RecoEcalCandidateIsolationMap >( "Dphi" ).setBranchAlias( "dphi" );
00064 }
00065
00066 EgammaHLTGsfTrackVarProducer::~EgammaHLTGsfTrackVarProducer(){}
00067
00068
00069
00070
00071
00072
00073
00074 void
00075 EgammaHLTGsfTrackVarProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00076 {
00077
00078 trackExtrapolator_.setup(iSetup);
00079
00080
00081 edm::Handle<reco::RecoEcalCandidateCollection> recoEcalCandHandle;
00082 iEvent.getByLabel(recoEcalCandTag_,recoEcalCandHandle);
00083
00084 edm::Handle<reco::ElectronCollection> electronHandle;
00085 iEvent.getByLabel(inputCollectionTag_,electronHandle);
00086
00087 edm::Handle<reco::GsfTrackCollection> gsfTracksHandle;
00088 if(!electronHandle.isValid()) iEvent.getByLabel (inputCollectionTag_,gsfTracksHandle);
00089
00090 edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
00091 iEvent.getByLabel(beamSpotTag_,recoBeamSpotHandle);
00092
00093 const reco::BeamSpot& beamSpot = *recoBeamSpotHandle;
00094
00095 edm::ESHandle<MagneticField> theMagField;
00096 iSetup.get<IdealMagneticFieldRecord>().get(theMagField);
00097
00098 reco::RecoEcalCandidateIsolationMap dEtaMap;
00099 reco::RecoEcalCandidateIsolationMap dPhiMap;
00100
00101 for(reco::RecoEcalCandidateCollection::const_iterator iRecoEcalCand = recoEcalCandHandle->begin(); iRecoEcalCand != recoEcalCandHandle->end(); iRecoEcalCand++){
00102 reco::RecoEcalCandidateRef recoEcalCandRef(recoEcalCandHandle,iRecoEcalCand-recoEcalCandHandle->begin());
00103
00104 const reco::SuperClusterRef scRef = recoEcalCandRef->superCluster();
00105
00106
00107 std::vector<const reco::GsfTrack*> gsfTracks;
00108 if(electronHandle.isValid()){
00109 for(reco::ElectronCollection::const_iterator eleIt = electronHandle->begin(); eleIt != electronHandle->end(); eleIt++){
00110 if(eleIt->superCluster()==scRef){
00111 gsfTracks.push_back(&*eleIt->gsfTrack());
00112 }
00113 }
00114 }else{
00115 for(reco::GsfTrackCollection::const_iterator trkIt =gsfTracksHandle->begin();trkIt!=gsfTracksHandle->end();++trkIt){
00116 edm::RefToBase<TrajectorySeed> seed = trkIt->extra()->seedRef() ;
00117 reco::ElectronSeedRef elseed = seed.castTo<reco::ElectronSeedRef>() ;
00118 edm::RefToBase<reco::CaloCluster> caloCluster = elseed->caloCluster() ;
00119 reco::SuperClusterRef scRefFromTrk = caloCluster.castTo<reco::SuperClusterRef>() ;
00120 if(scRefFromTrk==scRef){
00121 gsfTracks.push_back(&*trkIt);
00122 }
00123 }
00124
00125 }
00126 float dEtaInValue=999999;
00127 float dPhiInValue=999999;
00128
00129 if(static_cast<int>(gsfTracks.size())>=upperTrackNrToRemoveCut_){
00130 dEtaInValue=0;
00131 dPhiInValue=0;
00132 }else if(static_cast<int>(gsfTracks.size())<=lowerTrackNrToRemoveCut_){
00133 dEtaInValue=0;
00134 dPhiInValue=0;
00135 }else{
00136 for(size_t trkNr=0;trkNr<gsfTracks.size();trkNr++){
00137
00138 GlobalPoint scPos(scRef->x(),scRef->y(),scRef->z());
00139 GlobalPoint trackExtrapToSC = trackExtrapolator_.extrapolateTrackPosToPoint(*gsfTracks[trkNr],scPos);
00140 EleRelPointPair scAtVtx(scRef->position(),trackExtrapToSC,beamSpot.position());
00141
00142 if(fabs(scAtVtx.dEta())<dEtaInValue) dEtaInValue=fabs(scAtVtx.dEta());
00143 if(fabs(scAtVtx.dPhi())<dPhiInValue) dPhiInValue=fabs(scAtVtx.dPhi());
00144 }
00145 }
00146
00147 dEtaMap.insert(recoEcalCandRef, dEtaInValue);
00148 dPhiMap.insert(recoEcalCandRef, dPhiInValue);
00149 }
00150
00151 std::auto_ptr<reco::RecoEcalCandidateIsolationMap> dEtaMapForEvent(new reco::RecoEcalCandidateIsolationMap(dEtaMap));
00152 std::auto_ptr<reco::RecoEcalCandidateIsolationMap> dPhiMapForEvent(new reco::RecoEcalCandidateIsolationMap(dPhiMap));
00153 iEvent.put(dEtaMapForEvent, "Deta" );
00154 iEvent.put(dPhiMapForEvent, "Dphi" );
00155
00156 }
00157
00158
00159
00160 EgammaHLTGsfTrackVarProducer::TrackExtrapolator::TrackExtrapolator(const EgammaHLTGsfTrackVarProducer::TrackExtrapolator& rhs):
00161 cacheIDTDGeom_(rhs.cacheIDTDGeom_),
00162 cacheIDMagField_(rhs.cacheIDMagField_),
00163 magField_(rhs.magField_),
00164 trackerHandle_(rhs.trackerHandle_),
00165 mtsMode_(rhs.mtsMode_)
00166
00167 {
00168 if(rhs.mtsTransform_) mtsTransform_ = new MultiTrajectoryStateTransform(*rhs.mtsTransform_);
00169 else mtsTransform_ =0;
00170
00171 }
00172
00173 EgammaHLTGsfTrackVarProducer::TrackExtrapolator* EgammaHLTGsfTrackVarProducer::TrackExtrapolator::operator=(const EgammaHLTGsfTrackVarProducer::TrackExtrapolator& rhs)
00174 {
00175 if(this!=&rhs){
00176 cacheIDTDGeom_ = rhs.cacheIDTDGeom_;
00177 cacheIDMagField_ = rhs.cacheIDMagField_;
00178 magField_ = rhs.magField_;
00179 trackerHandle_ = rhs.trackerHandle_;
00180 mtsMode_ = rhs.mtsMode_;
00181
00182 delete mtsTransform_;
00183 if(rhs.mtsTransform_) mtsTransform_ = new MultiTrajectoryStateTransform(*rhs.mtsTransform_);
00184 else mtsTransform_ =0;
00185 }
00186 return this;
00187 }
00188
00189 void EgammaHLTGsfTrackVarProducer::TrackExtrapolator::setup(const edm::EventSetup& iSetup)
00190 {
00191 bool updateField(false);
00192 if (cacheIDMagField_!=iSetup.get<IdealMagneticFieldRecord>().cacheIdentifier()){
00193 updateField = true;
00194 cacheIDMagField_=iSetup.get<IdealMagneticFieldRecord>().cacheIdentifier();
00195 iSetup.get<IdealMagneticFieldRecord>().get(magField_);
00196 }
00197
00198 bool updateGeometry(false);
00199 if (cacheIDTDGeom_!=iSetup.get<TrackerDigiGeometryRecord>().cacheIdentifier()){
00200 updateGeometry = true;
00201 cacheIDTDGeom_=iSetup.get<TrackerDigiGeometryRecord>().cacheIdentifier();
00202 iSetup.get<TrackerDigiGeometryRecord>().get(trackerHandle_);
00203 }
00204
00205 if ( updateField || updateGeometry || !mtsTransform_ ) {
00206 delete mtsTransform_;
00207 mtsTransform_ = new MultiTrajectoryStateTransform(trackerHandle_.product(),magField_.product());
00208 }
00209 }
00210
00211 GlobalPoint EgammaHLTGsfTrackVarProducer::TrackExtrapolator::extrapolateTrackPosToPoint(const reco::GsfTrack& gsfTrack,const GlobalPoint& pointToExtrapTo)
00212 {
00213 TrajectoryStateOnSurface innTSOS = mtsTransform()->innerStateOnSurface(gsfTrack);
00214 TrajectoryStateOnSurface posTSOS = mtsTransform()->extrapolatedState(innTSOS,pointToExtrapTo);
00215 GlobalPoint extrapolatedPos;
00216 mtsMode()->positionFromModeCartesian(posTSOS,extrapolatedPos);
00217 return extrapolatedPos;
00218 }
00219
00220 GlobalVector EgammaHLTGsfTrackVarProducer::TrackExtrapolator::extrapolateTrackMomToPoint(const reco::GsfTrack& gsfTrack,const GlobalPoint& pointToExtrapTo)
00221 {
00222 TrajectoryStateOnSurface innTSOS = mtsTransform()->innerStateOnSurface(gsfTrack);
00223 TrajectoryStateOnSurface posTSOS = mtsTransform()->extrapolatedState(innTSOS,pointToExtrapTo);
00224 GlobalVector extrapolatedMom;
00225 mtsMode()->momentumFromModeCartesian(posTSOS,extrapolatedMom);
00226 return extrapolatedMom;
00227 }
00228
00229
00230