ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitProd/TreeFiller/src/FillerElectrons.cc
Revision: 1.36
Committed: Thu Jun 18 23:01:46 2009 UTC (15 years, 10 months ago) by bendavid
Content type: text/plain
Branch: MAIN
Changes since 1.35: +40 -25 lines
Log Message:
First attempt at updated electron id and isolation for 31X, updated supercluster link filling

File Contents

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