CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoEgamma/EgammaHLTProducers/src/EgammaHLTGsfTrackVarProducer.cc

Go to the documentation of this file.
00001 
00009 #include "RecoEgamma/EgammaHLTProducers/interface/EgammaHLTGsfTrackVarProducer.h"
00010 
00011 // Framework
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 //#include "DataFormats/EgammaCandidates/interface/ElectronIsolationAssociation.h"
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 //#include "CondFormats/DataRecord/interface/BeamSpotObjectsRcd.h"//needed?
00037 //#include "CondFormats/BeamSpotObjects/interface/BeamSpotObjects.h"//needed?
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"); //zeros out dEtaIn,dPhiIn if nrTracks>= this
00057   lowerTrackNrToRemoveCut_  = config.getParameter<int>("lowerTrackNrToRemoveCut"); //zeros out dEtaIn,dPhiIn if nrTracks<= this
00058  
00059   
00060   
00061   //register your products
00062   produces < reco::RecoEcalCandidateIsolationMap >( "Deta" ).setBranchAlias( "deta" );
00063   produces < reco::RecoEcalCandidateIsolationMap >( "Dphi" ).setBranchAlias( "dphi" ); 
00064 }
00065 
00066 EgammaHLTGsfTrackVarProducer::~EgammaHLTGsfTrackVarProducer(){}
00067 
00068 
00069 //
00070 // member functions
00071 //
00072 
00073 // ------------ method called to produce the data  ------------
00074 void
00075 EgammaHLTGsfTrackVarProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00076 {
00077 
00078   trackExtrapolator_.setup(iSetup);
00079 
00080   // Get the HLT filtered objects
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   // gets its position
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     //the idea is that we can take the tracks from properly associated electrons or just take all gsf tracks with that sc as a seed
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()); //we are allowing them to come from different tracks
00143         if(fabs(scAtVtx.dPhi())<dPhiInValue) dPhiInValue=fabs(scAtVtx.dPhi());//we are allowing them to come from different tracks
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){ //just to ensure we're not copying ourselves
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 //define this as a plug-in
00230 //DEFINE_FWK_MODULE(EgammaHLTTrackIsolationProducers);