ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitProd/TreeFiller/src/FillerElectrons.cc
Revision: 1.32
Committed: Mon Apr 6 19:37:20 2009 UTC (16 years, 1 month ago) by loizides
Content type: text/plain
Branch: MAIN
Changes since 1.31: +6 -4 lines
Log Message:
Coding conv.

File Contents

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