ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitProd/TreeFiller/src/FillerElectrons.cc
Revision: 1.21
Committed: Thu Nov 6 13:09:23 2008 UTC (16 years, 5 months ago) by sixie
Content type: text/plain
Branch: MAIN
Changes since 1.20: +3 -2 lines
Log Message:
Update tower isolation to work with new tag from Egamma group

File Contents

# User Rev Content
1 sixie 1.21 // $Id: FillerElectrons.cc,v 1.20 2008/11/05 10:47:34 bendavid Exp $
2 loizides 1.1
3     #include "MitProd/TreeFiller/interface/FillerElectrons.h"
4     #include "FWCore/MessageLogger/interface/MessageLogger.h"
5     #include "DataFormats/Common/interface/Handle.h"
6     #include "DataFormats/TrackReco/interface/Track.h"
7     #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
8     #include "DataFormats/TrackReco/interface/TrackFwd.h"
9 sixie 1.11 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
10     #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
11     #include "DataFormats/EgammaReco/interface/ClusterShape.h"
12     #include "DataFormats/EgammaReco/interface/BasicClusterShapeAssociation.h"
13 bendavid 1.18 #include "DataFormats/Common/interface/RefToPtr.h"
14 sixie 1.11 #include "RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h"
15     #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaEcalIsolation.h"
16 sixie 1.14 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTowerIsolation.h"
17 bendavid 1.19 #include "MitEdm/Producers/interface/RefToBaseToPtr.h"
18 loizides 1.1 #include "MitAna/DataTree/interface/Track.h"
19     #include "MitAna/DataTree/interface/Names.h"
20 sixie 1.12 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
21     #include "AnalysisDataFormats/Egamma/interface/ElectronID.h"
22     #include "AnalysisDataFormats/Egamma/interface/ElectronIDAssociation.h"
23 loizides 1.1
24 sixie 1.13 #include "DataFormats/Common/interface/ValueMap.h"
25    
26 loizides 1.1 using namespace std;
27     using namespace edm;
28     using namespace mithep;
29    
30 loizides 1.5 //--------------------------------------------------------------------------------------------------
31 loizides 1.10 FillerElectrons::FillerElectrons(const edm::ParameterSet &cfg, bool active) :
32 loizides 1.8 BaseFiller(cfg,"Electrons",active),
33 loizides 1.1 edmName_(Conf().getUntrackedParameter<string>("edmName","pixelMatchGsfElectrons")),
34     mitName_(Conf().getUntrackedParameter<string>("mitName",Names::gkElectronBrn)),
35 loizides 1.10 gsfTrackMapName_(Conf().getUntrackedParameter<string>("gsfTrackMapName","")),
36     trackerTrackMapName_(Conf().getUntrackedParameter<string>("trackerTrackMapName","")),
37 sixie 1.12 barrelEcalRecHitName_(Conf().getUntrackedParameter<string>("barrelEcalRecHitName","")),
38     endcapEcalRecHitName_(Conf().getUntrackedParameter<string>("endcapEcalRecHitName","")),
39     barrelSuperClusterName_(Conf().getUntrackedParameter<string>("barrelSuperClusterName", "")),
40     endcapSuperClusterName_(Conf().getUntrackedParameter<string>("endcapSuperClusterName", "")),
41     barrelBasicClusterName_(Conf().getUntrackedParameter<string>("barrelBasicClusterName", "")),
42     endcapBasicClusterName_(Conf().getUntrackedParameter<string>("endcapBasicClusterName", "")),
43     barrelSuperClusterMapName_(Conf().getUntrackedParameter<string>("barrelSuperClusterMapName","")),
44     endcapSuperClusterMapName_(Conf().getUntrackedParameter<string>("endcapSuperClusterMapName","")),
45 loizides 1.16 eIDCutBasedTightName_(Conf().getUntrackedParameter<string>("eIDCutBasedTightName","eidTight")),
46 sixie 1.13 eIDCutBasedLooseName_(Conf().getUntrackedParameter<string>("eIDCutBasedLooseName","eidLoose")),
47 sixie 1.12 eIDLikelihoodName_(Conf().getUntrackedParameter<string>("eIDLikelihood","eidLikelihood")),
48     eIDNeuralNetName_(Conf().getUntrackedParameter<string>("eIDNeuralNet","eidNeuralNet")),
49 loizides 1.16 isoTrackColName_(Conf().getUntrackedParameter<string>
50     ("IsolationTrackCollectionName","generalTracks")),
51     isoCaloTowerColName_(Conf().getUntrackedParameter<string>
52     ("IsolationCaloTowerCollectionName","towerMaker")),
53     ecalJurassicIsoName_(Conf().getUntrackedParameter<string>
54     ("EcalJurassicIsolationName","eleIsoFromDepsEcalFromHits")),
55     hcalJurassicIsoName_(Conf().getUntrackedParameter<string>
56     ("HcalJurassicIsolationName","eleIsoFromDepsHcalFromHits")),
57 peveraer 1.17 trackerIsoName_(Conf().getUntrackedParameter<string>
58     ("TrackerIsolationName","eleIsoFromDepsTk")),
59 bendavid 1.19 gsfTrackAssocName_(Conf().getUntrackedParameter<string>("gsfTrackAssocName","")),
60 loizides 1.8 electrons_(new mithep::ElectronArr(16)),
61 loizides 1.10 gsfTrackMap_(0),
62 sixie 1.11 trackerTrackMap_(0),
63 sixie 1.12 barrelSuperClusterMap_(0),
64     endcapSuperClusterMap_(0)
65 loizides 1.1 {
66     // Constructor.
67     }
68    
69 loizides 1.5 //--------------------------------------------------------------------------------------------------
70 loizides 1.1 FillerElectrons::~FillerElectrons()
71     {
72     // Destructor.
73 loizides 1.6
74     delete electrons_;
75 loizides 1.1 }
76    
77 loizides 1.5 //--------------------------------------------------------------------------------------------------
78 loizides 1.1 void FillerElectrons::BookDataBlock(TreeWriter &tws)
79     {
80 loizides 1.10 // Add electron branch to our tree and get our maps.
81 loizides 1.1
82     tws.AddBranch(mitName_.c_str(),&electrons_);
83 loizides 1.10
84     if (!gsfTrackMapName_.empty())
85 bendavid 1.18 gsfTrackMap_ = OS()->get<TrackMap>(gsfTrackMapName_.c_str());
86 loizides 1.10 if (!trackerTrackMapName_.empty())
87 sixie 1.11 trackerTrackMap_ = OS()->get<TrackMap>(trackerTrackMapName_.c_str());
88 sixie 1.12 if (!barrelSuperClusterMapName_.empty())
89     barrelSuperClusterMap_ = OS()->get<SuperClusterMap>(barrelSuperClusterMapName_.c_str());
90     if (!endcapSuperClusterMapName_.empty())
91     endcapSuperClusterMap_ = OS()->get<SuperClusterMap>(endcapSuperClusterMapName_.c_str());
92 loizides 1.1 }
93    
94 loizides 1.5 //--------------------------------------------------------------------------------------------------
95 sixie 1.11 void FillerElectrons::FillDataBlock(const edm::Event &event, const edm::EventSetup &setup)
96 loizides 1.1 {
97     // Fill electrons from edm collection into our collection.
98    
99     electrons_->Reset();
100    
101 sixie 1.11 Handle<reco::GsfElectronCollection> hElectronProduct;
102 loizides 1.7 GetProduct(edmName_, hElectronProduct, event);
103 loizides 1.1
104 sixie 1.13 //Handles to get the electron ID information
105     Handle<edm::ValueMap<float> > eidLooseMap;
106     GetProduct(eIDCutBasedLooseName_, eidLooseMap, event);
107     Handle<edm::ValueMap<float> > eidTightMap;
108     GetProduct(eIDCutBasedTightName_, eidTightMap, event);
109     Handle<edm::ValueMap<float> > eidLikelihoodMap;
110     GetProduct(eIDLikelihoodName_, eidLikelihoodMap, event);
111 bendavid 1.19
112     //get gsf track association map if needed
113     mitedm::TrackAssociation gsfAssociation;
114     if (trackerTrackMap_ && !gsfTrackAssocName_.empty()) {
115     Handle<mitedm::TrackAssociation> gsfAssociationProduct;
116     GetProduct(gsfTrackAssocName_, gsfAssociationProduct, event);
117     gsfAssociation = *(gsfAssociationProduct.product());
118     }
119    
120 sixie 1.13
121 sixie 1.11 const reco::GsfElectronCollection inElectrons = *(hElectronProduct.product());
122     //loop over electrons
123 sixie 1.13
124 sixie 1.11 for (reco::GsfElectronCollection::const_iterator iM = inElectrons.begin();
125 loizides 1.1 iM != inElectrons.end(); ++iM) {
126 sixie 1.13
127     //the index and Ref are needed for the eID association Map
128     unsigned int iElectron = iM - inElectrons.begin();
129     reco::GsfElectronRef eRef(hElectronProduct, iElectron);
130    
131 loizides 1.8 mithep::Electron *outElectron = electrons_->AddNew();
132 sixie 1.12
133 sixie 1.11 outElectron->SetESuperClusterOverP( iM->eSuperClusterOverP() ) ;
134     outElectron->SetESeedClusterOverPout( iM->eSeedClusterOverPout() ) ;
135 sixie 1.13 outElectron->SetPIn( iM->trackMomentumAtVtx().R() ) ;
136     outElectron->SetPOut( iM->trackMomentumOut().R() );
137 sixie 1.11 outElectron->SetDeltaEtaSuperClusterTrackAtVtx( iM->deltaEtaSuperClusterTrackAtVtx() ) ;
138     outElectron->SetDeltaEtaSeedClusterTrackAtCalo( iM->deltaEtaSeedClusterTrackAtCalo() ) ;
139     outElectron->SetDeltaPhiSuperClusterTrackAtVtx( iM->deltaPhiSuperClusterTrackAtVtx() ) ;
140     outElectron->SetDeltaPhiSeedClusterTrackAtCalo( iM->deltaPhiSeedClusterTrackAtCalo() ) ;
141     outElectron->SetHadronicOverEm( iM->hadronicOverEm() ) ;
142     outElectron->SetIsEnergyScaleCorrected( iM->isEnergyScaleCorrected() ) ;
143     outElectron->SetIsMomentumCorrected( iM->isMomentumCorrected() ) ;
144     outElectron->SetNumberOfClusters( iM->numberOfClusters() ) ;
145     outElectron->SetClassification( iM->classification() ) ;
146    
147 loizides 1.16 //shower shape variables
148     EcalClusterLazyTools lazyTools(event, setup, edm::InputTag(barrelEcalRecHitName_),
149     edm::InputTag(endcapEcalRecHitName_));
150     outElectron->SetE33(lazyTools.e3x3(*(iM->superCluster()->seed())));
151     outElectron->SetE55(lazyTools.e5x5(*(iM->superCluster()->seed())));
152     std::vector<float> vCov = lazyTools.covariances(*(iM->superCluster()->seed()));
153 sixie 1.12 outElectron->SetCovEtaEta(vCov[0]);
154     outElectron->SetCovEtaPhi(vCov[1]);
155     outElectron->SetCovPhiPhi(vCov[2]);
156    
157 loizides 1.16 //compute isolations
158     //get the barrel BasicClusters
159 sixie 1.12 edm::Handle<reco::BasicClusterCollection> barrelBasicClusterHandle;
160     GetProduct(barrelBasicClusterName_, barrelBasicClusterHandle, event);
161 sixie 1.11 const reco::BasicClusterCollection* barrelBasicClusters = barrelBasicClusterHandle.product();
162    
163 loizides 1.16 //get the endcap BasicClusters
164 sixie 1.11 edm::Handle<reco::BasicClusterCollection> endcapBasicClusterHandle;
165 sixie 1.12 GetProduct(endcapBasicClusterName_, endcapBasicClusterHandle, event);
166 sixie 1.11 const reco::BasicClusterCollection* endcapBasicClusters = endcapBasicClusterHandle.product();
167    
168 loizides 1.16 //get the barrel SuperClusters
169 sixie 1.12 edm::Handle<reco::SuperClusterCollection> barrelSuperClusterHandle;
170     GetProduct(barrelSuperClusterName_, barrelSuperClusterHandle, event);
171     const reco::SuperClusterCollection* barrelSuperClusters = barrelSuperClusterHandle.product();
172    
173 loizides 1.16 //Get the endcap SuperClusters
174 sixie 1.12 edm::Handle<reco::SuperClusterCollection> endcapSuperClusterHandle;
175     GetProduct(endcapSuperClusterName_, endcapSuperClusterHandle, event);
176     const reco::SuperClusterCollection* endcapSuperClusters = endcapSuperClusterHandle.product();
177    
178     //find out whether this electron super cluster is in the barrel or endcap
179     bool isBarrel=false ;
180     if(barrelSuperClusterMap_->HasMit(iM->superCluster()))
181     isBarrel = true;
182    
183 sixie 1.11 //compute ECAL isolation
184     double extRadius = 0.3;
185     double etLow = 0.0;
186 loizides 1.16 EgammaEcalIsolation *myEcalIsolation = 0;
187 sixie 1.12 if (!isBarrel) {
188     myEcalIsolation =
189     new EgammaEcalIsolation(extRadius,etLow,endcapBasicClusters,endcapSuperClusters);
190 loizides 1.16 } else {
191     myEcalIsolation =
192     new EgammaEcalIsolation(extRadius,etLow,barrelBasicClusters,barrelSuperClusters);
193 sixie 1.12 }
194 loizides 1.16
195 sixie 1.12 double ecalIsoValue = myEcalIsolation->getEcalEtSum(&(*iM));
196    
197 loizides 1.16 //compute CaloTower isolation
198 sixie 1.14 edm::Handle<CaloTowerCollection> caloTowers;
199 loizides 1.16 GetProduct(isoCaloTowerColName_, caloTowers, event);
200 sixie 1.14 extRadius = 0.3;
201 loizides 1.16 etLow = 0.0;
202 sixie 1.14 double intRadius = 0.02;
203 sixie 1.21 int hcalDepth = -1; //-1 means we take all depths.
204 sixie 1.14 EgammaTowerIsolation *myTowerIsolation =
205 sixie 1.21 new EgammaTowerIsolation (extRadius, intRadius, etLow, hcalDepth, caloTowers.product());
206 sixie 1.14 double towerIsoValue = myTowerIsolation->getTowerEtSum(&(*iM));
207     outElectron->SetCaloTowerIsolation( towerIsoValue );
208 peveraer 1.17
209 loizides 1.16 //fill the isolation values
210 sixie 1.11 outElectron->SetCaloIsolation( ecalIsoValue ) ;
211 peveraer 1.17
212     //get and fill Track isolation
213     Handle<edm::ValueMap<double> > eleIsoFromDepsTkValueMap;
214     GetProduct(trackerIsoName_, eleIsoFromDepsTkValueMap, event);
215     outElectron->SetTrackIsolation((*eleIsoFromDepsTkValueMap)[eRef]);
216 sixie 1.11
217 loizides 1.16 //get and fill Jurassic isolation values
218 sixie 1.15 Handle<edm::ValueMap<double> > eleIsoFromDepsEcalFromHitsValueMap;
219 loizides 1.16 GetProduct(ecalJurassicIsoName_, eleIsoFromDepsEcalFromHitsValueMap, event);
220 sixie 1.15 Handle<edm::ValueMap<double> > eleIsoFromDepsHcalFromHitsValueMap;
221 loizides 1.16 GetProduct(hcalJurassicIsoName_, eleIsoFromDepsHcalFromHitsValueMap, event);
222 sixie 1.15
223     outElectron->SetEcalJurassicIsolation((*eleIsoFromDepsEcalFromHitsValueMap)[eRef]);
224     outElectron->SetHcalJurassicIsolation((*eleIsoFromDepsHcalFromHitsValueMap)[eRef]);
225    
226 loizides 1.16 //make proper links to Tracks and Super Clusters
227 loizides 1.1 if (gsfTrackMap_ && iM->gsfTrack().isNonnull())
228 bendavid 1.18 outElectron->SetGsfTrk(gsfTrackMap_->GetMit(refToPtr(iM->gsfTrack())));
229 bendavid 1.19 //make tracker track links, relinking from gsf track associations if configured and
230     //link is otherwise absent
231     if (trackerTrackMap_) {
232     if (iM->track().isNonnull())
233     outElectron->SetTrackerTrk(trackerTrackMap_->GetMit(refToPtr(iM->track())));
234     else if (!gsfTrackAssocName_.empty() && iM->gsfTrack().isNonnull()) {
235     reco::TrackBaseRef gsfRef(iM->gsfTrack());
236     std::vector<std::pair<reco::TrackBaseRef, double> > matchedTracks;
237     try {
238     matchedTracks = gsfAssociation[gsfRef];
239     }
240     catch (edm::Exception &ex) {
241     }
242     //take the best match, but only if more than 50% of the hits came from the original
243     //gsf track
244     reco::TrackBaseRef trackerRef;
245     double rMax = 0.0;
246 bendavid 1.20 for (uint imatch=0; imatch<matchedTracks.size(); ++imatch) {
247 bendavid 1.19 std::pair<reco::TrackBaseRef, double> &match = matchedTracks.at(imatch);
248     double r = match.second;
249     if ( r>rMax && r>0.5 ) {
250     rMax = r;
251     trackerRef = match.first;
252     }
253     }
254     if (trackerRef.isNonnull())
255     outElectron->SetTrackerTrk(trackerTrackMap_->GetMit(mitedm::refToBaseToPtr(trackerRef)));
256     }
257     }
258 sixie 1.12 if (barrelSuperClusterMap_ && endcapSuperClusterMap_ && iM->superCluster().isNonnull())
259     if(isBarrel) {
260     outElectron->SetSuperCluster(barrelSuperClusterMap_->GetMit(iM->superCluster()));
261     } else {
262     outElectron->SetSuperCluster(endcapSuperClusterMap_->GetMit(iM->superCluster()));
263     }
264 sixie 1.13
265 loizides 1.16 //fill Electron ID information
266 sixie 1.13 outElectron->SetPassLooseID((*eidLooseMap)[eRef]);
267     outElectron->SetPassTightID((*eidTightMap)[eRef]);
268     outElectron->SetIDLikelihood((*eidLikelihoodMap)[eRef]);
269 sixie 1.12 }
270 loizides 1.1 electrons_->Trim();
271     }