ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitProd/TreeFiller/src/FillerElectrons.cc
Revision: 1.44
Committed: Sat May 29 11:58:28 2010 UTC (14 years, 11 months ago) by bendavid
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_014a, Mit_014
Changes since 1.43: +85 -1 lines
Log Message:
Fill electron gsf vertex compatibility variables

File Contents

# User Rev Content
1 bendavid 1.44 // $Id: FillerElectrons.cc,v 1.43 2010/05/06 17:31:24 bendavid 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 bendavid 1.43 #include "DataFormats/VertexReco/interface/VertexFwd.h"
20     #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
21     #include "TrackingTools/TransientTrack/plugins/TransientTrackBuilderESProducer.h"
22 bendavid 1.44 #include "RecoVertex/GaussianSumVertexFit/interface/GsfVertexTrackCompatibilityEstimator.h"
23 bendavid 1.43 #include "TrackingTools/IPTools/interface/IPTools.h"
24     #include "RecoEgamma/EgammaTools/interface/ConversionFinder.h"
25     #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
26     #include "MagneticField/Engine/interface/MagneticField.h"
27 loizides 1.35 #include "MitAna/DataTree/interface/ElectronCol.h"
28     #include "MitAna/DataTree/interface/Names.h"
29     #include "MitAna/DataTree/interface/Track.h"
30 loizides 1.33 #include "MitEdm/DataFormats/interface/RefToBaseToPtr.h"
31 loizides 1.35 #include "MitProd/ObjectService/interface/ObjectService.h"
32 sixie 1.13
33 loizides 1.1 using namespace std;
34     using namespace edm;
35     using namespace mithep;
36    
37 loizides 1.5 //--------------------------------------------------------------------------------------------------
38 loizides 1.28 FillerElectrons::FillerElectrons(const edm::ParameterSet &cfg, const char *name, bool active) :
39     BaseFiller(cfg,name,active),
40 loizides 1.1 edmName_(Conf().getUntrackedParameter<string>("edmName","pixelMatchGsfElectrons")),
41     mitName_(Conf().getUntrackedParameter<string>("mitName",Names::gkElectronBrn)),
42 loizides 1.10 gsfTrackMapName_(Conf().getUntrackedParameter<string>("gsfTrackMapName","")),
43     trackerTrackMapName_(Conf().getUntrackedParameter<string>("trackerTrackMapName","")),
44 sixie 1.12 barrelSuperClusterMapName_(Conf().getUntrackedParameter<string>("barrelSuperClusterMapName","")),
45     endcapSuperClusterMapName_(Conf().getUntrackedParameter<string>("endcapSuperClusterMapName","")),
46 bendavid 1.36 pfSuperClusterMapName_(Conf().getUntrackedParameter<string>("pfSuperClusterMapName","")),
47 loizides 1.16 eIDCutBasedTightName_(Conf().getUntrackedParameter<string>("eIDCutBasedTightName","eidTight")),
48 bendavid 1.43 eIDCutBasedLooseName_(Conf().getUntrackedParameter<string>("eIDCutBasedLooseName","eidLoose")),
49     pvEdmName_(Conf().getUntrackedParameter<string>("pvEdmName","offlinePrimaryVertices")),
50     pvBSEdmName_(Conf().getUntrackedParameter<string>("pvEdmName","offlinePrimaryVerticesWithBS")),
51 loizides 1.8 electrons_(new mithep::ElectronArr(16)),
52 loizides 1.10 gsfTrackMap_(0),
53 sixie 1.11 trackerTrackMap_(0),
54 sixie 1.12 barrelSuperClusterMap_(0),
55     endcapSuperClusterMap_(0)
56 loizides 1.1 {
57     // Constructor.
58     }
59    
60 loizides 1.5 //--------------------------------------------------------------------------------------------------
61 loizides 1.1 FillerElectrons::~FillerElectrons()
62     {
63     // Destructor.
64 loizides 1.6
65     delete electrons_;
66 loizides 1.1 }
67    
68 loizides 1.5 //--------------------------------------------------------------------------------------------------
69 bendavid 1.42 void FillerElectrons::BookDataBlock(TreeWriter &tws)
70 loizides 1.1 {
71 loizides 1.10 // Add electron branch to our tree and get our maps.
72 loizides 1.1
73 loizides 1.29 tws.AddBranch(mitName_,&electrons_);
74     OS()->add<mithep::ElectronArr>(electrons_,mitName_);
75 loizides 1.10
76 loizides 1.29 if (!gsfTrackMapName_.empty()) {
77     gsfTrackMap_ = OS()->get<TrackMap>(gsfTrackMapName_);
78     if (gsfTrackMap_)
79     AddBranchDep(mitName_,gsfTrackMap_->GetBrName());
80     }
81     if (!trackerTrackMapName_.empty()) {
82     trackerTrackMap_ = OS()->get<TrackMap>(trackerTrackMapName_);
83     if (trackerTrackMap_)
84     AddBranchDep(mitName_,trackerTrackMap_->GetBrName());
85     }
86     if (!barrelSuperClusterMapName_.empty()) {
87     barrelSuperClusterMap_ = OS()->get<SuperClusterMap>(barrelSuperClusterMapName_);
88     if (barrelSuperClusterMap_)
89     AddBranchDep(mitName_,barrelSuperClusterMap_->GetBrName());
90     }
91     if (!endcapSuperClusterMapName_.empty()) {
92     endcapSuperClusterMap_ = OS()->get<SuperClusterMap>(endcapSuperClusterMapName_);
93     if (endcapSuperClusterMap_)
94     AddBranchDep(mitName_,endcapSuperClusterMap_->GetBrName());
95     }
96 bendavid 1.36 if (!pfSuperClusterMapName_.empty()) {
97     pfSuperClusterMap_ = OS()->get<SuperClusterMap>(pfSuperClusterMapName_);
98     if (pfSuperClusterMap_)
99     AddBranchDep(mitName_,pfSuperClusterMap_->GetBrName());
100     }
101 loizides 1.1 }
102    
103 loizides 1.5 //--------------------------------------------------------------------------------------------------
104 loizides 1.29 void FillerElectrons::FillDataBlock(const edm::Event &event, const edm::EventSetup &setup)
105 loizides 1.1 {
106     // Fill electrons from edm collection into our collection.
107    
108 bendavid 1.26 electrons_->Delete();
109 loizides 1.1
110 sixie 1.11 Handle<reco::GsfElectronCollection> hElectronProduct;
111 loizides 1.7 GetProduct(edmName_, hElectronProduct, event);
112 loizides 1.1
113 loizides 1.29 // handles to get the electron ID information
114 sixie 1.13 Handle<edm::ValueMap<float> > eidLooseMap;
115     GetProduct(eIDCutBasedLooseName_, eidLooseMap, event);
116     Handle<edm::ValueMap<float> > eidTightMap;
117     GetProduct(eIDCutBasedTightName_, eidTightMap, event);
118 bendavid 1.43
119     edm::Handle<reco::VertexCollection> hVertex;
120     event.getByLabel(pvEdmName_, hVertex);
121     const reco::VertexCollection *pvCol = hVertex.product();
122    
123     edm::Handle<reco::VertexCollection> hVertexBS;
124     event.getByLabel(pvBSEdmName_, hVertexBS);
125     const reco::VertexCollection *pvBSCol = hVertexBS.product();
126    
127     edm::Handle<reco::TrackCollection> hGeneralTracks;
128     event.getByLabel("generalTracks", hGeneralTracks);
129     //const reco::VertexCollection *trackCol = hGeneralTracks.product();
130    
131     edm::ESHandle<TransientTrackBuilder> hTransientTrackBuilder;
132     setup.get<TransientTrackRecord>().get("TransientTrackBuilder",hTransientTrackBuilder);
133     const TransientTrackBuilder *transientTrackBuilder = hTransientTrackBuilder.product();
134 bendavid 1.19
135 bendavid 1.44 GsfVertexTrackCompatibilityEstimator gsfEstimator;
136    
137 bendavid 1.43 //Get Magnetic Field from event setup, taking value at (0,0,0)
138     edm::ESHandle<MagneticField> magneticField;
139     setup.get<IdealMagneticFieldRecord>().get(magneticField);
140     const double bfield = magneticField->inTesla(GlobalPoint(0.,0.,0.)).z();
141    
142 sixie 1.11 const reco::GsfElectronCollection inElectrons = *(hElectronProduct.product());
143 loizides 1.29 // loop over electrons
144 sixie 1.11 for (reco::GsfElectronCollection::const_iterator iM = inElectrons.begin();
145 loizides 1.1 iM != inElectrons.end(); ++iM) {
146 sixie 1.13
147 loizides 1.29 // the index and Ref are needed for the eID association Map
148 sixie 1.13 unsigned int iElectron = iM - inElectrons.begin();
149     reco::GsfElectronRef eRef(hElectronProduct, iElectron);
150    
151 loizides 1.8 mithep::Electron *outElectron = electrons_->AddNew();
152 bendavid 1.31
153     outElectron->SetPtEtaPhi(iM->pt(),iM->eta(),iM->phi());
154 sixie 1.12
155 bendavid 1.40 outElectron->SetCharge(iM->charge());
156     outElectron->SetScPixCharge(iM->scPixCharge());
157    
158 loizides 1.30 outElectron->SetESuperClusterOverP(iM->eSuperClusterOverP());
159     outElectron->SetESeedClusterOverPout(iM->eSeedClusterOverPout());
160     outElectron->SetPIn(iM->trackMomentumAtVtx().R());
161 loizides 1.29 outElectron->SetPOut(iM->trackMomentumOut().R());
162 loizides 1.30 outElectron->SetDeltaEtaSuperClusterTrackAtVtx(iM->deltaEtaSuperClusterTrackAtVtx());
163     outElectron->SetDeltaEtaSeedClusterTrackAtCalo(iM->deltaEtaSeedClusterTrackAtCalo());
164     outElectron->SetDeltaPhiSuperClusterTrackAtVtx(iM->deltaPhiSuperClusterTrackAtVtx());
165     outElectron->SetDeltaPhiSeedClusterTrackAtCalo(iM->deltaPhiSeedClusterTrackAtCalo());
166     outElectron->SetIsEnergyScaleCorrected(iM->isEnergyScaleCorrected());
167     outElectron->SetIsMomentumCorrected(iM->isMomentumCorrected());
168 bendavid 1.36 outElectron->SetNumberOfClusters(iM->basicClustersSize());
169 loizides 1.30 outElectron->SetClassification(iM->classification());
170 bendavid 1.37 outElectron->SetFBrem(iM->fbrem());
171    
172 loizides 1.38 // pflow electron stuff
173 bendavid 1.41 outElectron->SetIsEcalDriven(iM->ecalDrivenSeed());
174     outElectron->SetIsTrackerDriven(iM->trackerDrivenSeed());
175 bendavid 1.37 outElectron->SetMva(iM->mva());
176    
177 loizides 1.29 // shower shape variables
178 bendavid 1.37 outElectron->SetE15(iM->e1x5());
179     outElectron->SetE25Max(iM->e2x5Max());
180     outElectron->SetE55(iM->e5x5());
181     outElectron->SetCovEtaEta(iM->sigmaEtaEta());
182     outElectron->SetCoviEtaiEta(iM->sigmaIetaIeta());
183     outElectron->SetHadronicOverEm(iM->hcalOverEcal());
184     outElectron->SetHcalDepth1OverEcal(iM->hcalDepth1OverEcal());
185     outElectron->SetHcalDepth2OverEcal(iM->hcalDepth2OverEcal());
186 sixie 1.12
187 loizides 1.38 // fill isolation variables for both cone sizes
188 bendavid 1.37 outElectron->SetEcalRecHitIsoDr04(iM->dr04EcalRecHitSumEt());
189     outElectron->SetHcalDepth1TowerSumEtDr04(iM->dr04HcalDepth1TowerSumEt());
190     outElectron->SetHcalDepth2TowerSumEtDr04(iM->dr04HcalDepth2TowerSumEt());
191     outElectron->SetTrackIsolationDr04(iM->dr04TkSumPt());
192     outElectron->SetEcalRecHitIsoDr03(iM->dr03EcalRecHitSumEt());
193     outElectron->SetHcalTowerSumEtDr03(iM->dr03HcalTowerSumEt());
194     outElectron->SetHcalDepth1TowerSumEtDr03(iM->dr03HcalDepth1TowerSumEt());
195     outElectron->SetHcalDepth2TowerSumEtDr03(iM->dr03HcalDepth2TowerSumEt());
196     outElectron->SetTrackIsolationDr03(iM->dr03TkSumPt());
197    
198 loizides 1.38 // fiducial flags
199 bendavid 1.37 outElectron->SetIsEB(iM->isEB());
200     outElectron->SetIsEE(iM->isEE());
201     outElectron->SetIsEBEEGap(iM->isEBEEGap());
202     outElectron->SetIsEBEtaGap(iM->isEBEtaGap());
203     outElectron->SetIsEBPhiGap(iM->isEBPhiGap());
204     outElectron->SetIsEEDeeGap(iM->isEEDeeGap());
205     outElectron->SetIsEERingGap(iM->isEERingGap());
206 sixie 1.12
207 loizides 1.38 // gsf-tracker match quality
208 bendavid 1.37 outElectron->SetFracSharedHits(iM->shFracInnerHits());
209 sixie 1.15
210 loizides 1.29 // make proper links to Tracks and Super Clusters
211 loizides 1.1 if (gsfTrackMap_ && iM->gsfTrack().isNonnull())
212 bendavid 1.18 outElectron->SetGsfTrk(gsfTrackMap_->GetMit(refToPtr(iM->gsfTrack())));
213 loizides 1.29 // make tracker track links, relinking from gsf track associations if configured and
214     // link is otherwise absent
215 bendavid 1.37 if (trackerTrackMap_ && iM->closestCtfTrackRef().isNonnull()) {
216 bendavid 1.36 outElectron->SetTrackerTrk(trackerTrackMap_->GetMit(refToPtr(iM->closestCtfTrackRef())));
217 bendavid 1.19 }
218 loizides 1.38 if (barrelSuperClusterMap_ && endcapSuperClusterMap_ &&
219     pfSuperClusterMap_ && iM->superCluster().isNonnull()) {
220 bendavid 1.37 if(barrelSuperClusterMap_->HasMit(iM->superCluster())) {
221 sixie 1.12 outElectron->SetSuperCluster(barrelSuperClusterMap_->GetMit(iM->superCluster()));
222 bendavid 1.36 }
223     else if (endcapSuperClusterMap_->HasMit(iM->superCluster())) {
224 sixie 1.12 outElectron->SetSuperCluster(endcapSuperClusterMap_->GetMit(iM->superCluster()));
225     }
226 bendavid 1.36 else if (pfSuperClusterMap_->HasMit(iM->superCluster())) {
227     outElectron->SetSuperCluster(pfSuperClusterMap_->GetMit(iM->superCluster()));
228     }
229     else throw edm::Exception(edm::errors::Configuration, "FillerElectrons:FillDataBlock()\n")
230     << "Error! SuperCluster reference in unmapped collection " << edmName_ << endl;
231     }
232 sixie 1.13
233 bendavid 1.43 //compute impact parameter with respect to PV
234     if (iM->gsfTrack().isNonnull()) {
235     const reco::TransientTrack &tt = transientTrackBuilder->build(iM->gsfTrack());
236    
237     const std::pair<bool,Measurement1D> &d0pv = IPTools::absoluteTransverseImpactParameter(tt,pvCol->at(0));
238     if (d0pv.first) {
239     outElectron->SetD0PV(d0pv.second.value());
240     outElectron->SetD0PVErr(d0pv.second.error());
241     }
242 bendavid 1.44 else {
243     outElectron->SetD0PV(-99.0);
244     }
245    
246 bendavid 1.43
247     const std::pair<bool,Measurement1D> &ip3dpv = IPTools::absoluteImpactParameter3D(tt,pvCol->at(0));
248     if (ip3dpv.first) {
249     outElectron->SetIp3dPV(ip3dpv.second.value());
250     outElectron->SetIp3dPVErr(ip3dpv.second.error());
251     }
252 bendavid 1.44 else {
253     outElectron->SetIp3dPV(-99.0);
254     }
255 bendavid 1.43
256     const std::pair<bool,Measurement1D> &d0pvbs = IPTools::absoluteTransverseImpactParameter(tt,pvBSCol->at(0));
257     if (d0pvbs.first) {
258     outElectron->SetD0PVBS(d0pvbs.second.value());
259     outElectron->SetD0PVBSErr(d0pvbs.second.error());
260     }
261 bendavid 1.44 else {
262     outElectron->SetD0PVBS(-99.0);
263     }
264 bendavid 1.43
265     const std::pair<bool,Measurement1D> &ip3dpvbs = IPTools::absoluteImpactParameter3D(tt,pvBSCol->at(0));
266     if (ip3dpvbs.first) {
267     outElectron->SetIp3dPVBS(ip3dpvbs.second.value());
268     outElectron->SetIp3dPVBSErr(ip3dpvbs.second.error());
269     }
270 bendavid 1.44 else {
271     outElectron->SetIp3dPVBS(-99.0);
272     }
273    
274     //compute compatibility with PV using full GSF state mixture (but skip in AOD)
275     if (iM->gsfTrack()->gsfExtra().isAvailable()) {
276    
277     const std::pair<bool,double> &pvGsfCompat = gsfEstimator.estimate(pvCol->at(0),tt);
278     if (pvGsfCompat.first) {
279     outElectron->SetGsfPVCompatibility(pvGsfCompat.second);
280     }
281     else {
282     outElectron->SetGsfPVCompatibility(-99.0);
283     }
284    
285    
286     const std::pair<bool,double> &pvbsGsfCompat = gsfEstimator.estimate(pvBSCol->at(0),tt);
287     if (pvbsGsfCompat.first) {
288     outElectron->SetGsfPVBSCompatibility(pvbsGsfCompat.second);
289     }
290     else {
291     outElectron->SetGsfPVBSCompatibility(-99.0);
292     }
293    
294     //compute signal vertex compatibility with full GSF state mixture excluding matching ckf track
295     //from vertex
296     if (iM->closestCtfTrackRef().isNonnull() && iM->closestCtfTrackRef()->extra().isAvailable()) {
297     const reco::TransientTrack &ttCkf = transientTrackBuilder->build(iM->closestCtfTrackRef());
298    
299     const std::pair<bool,double> &pvGsfCompatMatched = gsfEstimator.estimate(pvCol->at(0),tt, ttCkf);
300     if (pvGsfCompatMatched.first) {
301     outElectron->SetGsfPVCompatibilityMatched(pvGsfCompatMatched.second);
302     }
303     else {
304     outElectron->SetGsfPVCompatibilityMatched(-99.0);
305     }
306    
307    
308     const std::pair<bool,double> &pvbsGsfCompatMatched = gsfEstimator.estimate(pvBSCol->at(0),tt, ttCkf);
309     if (pvbsGsfCompatMatched.first) {
310     outElectron->SetGsfPVBSCompatibilityMatched(pvbsGsfCompatMatched.second);
311     }
312     else {
313     outElectron->SetGsfPVBSCompatibilityMatched(-99.0);
314     }
315    
316     if (verbose_>1) {
317     printf("gsf compat = %5f\n", pvGsfCompat.second);
318     printf("gsf compat matched = %5f\n", pvGsfCompatMatched.second);
319     }
320     }
321     else {
322     //no matching ckf track, so copy existing values
323     if (pvGsfCompat.first) {
324     outElectron->SetGsfPVCompatibilityMatched(pvGsfCompat.second);
325     }
326     else {
327     outElectron->SetGsfPVCompatibilityMatched(-99.0);
328     }
329    
330     if (pvbsGsfCompat.first) {
331     outElectron->SetGsfPVBSCompatibilityMatched(pvbsGsfCompat.second);
332     }
333     else {
334     outElectron->SetGsfPVBSCompatibilityMatched(-99.0);
335     }
336     }
337    
338     }
339 bendavid 1.43
340     if (verbose_>1) {
341     printf("gsf track pt = %5f\n",iM->gsfTrack()->pt());
342     printf("gsf track mode pt = %5f\n",iM->gsfTrack()->ptMode());
343     printf("ttrack pt = %5f\n",tt.initialFreeState().momentum().perp());
344     //printf("ttrackgsf pt = %5f\n",ttgsf.innermostMeasurementState().globalMomentum().perp());
345 bendavid 1.44 printf("ip3dpv reduced chisquared = %5f, probability = %5f\n", ip3dpv.second.value()/ip3dpv.second.error(), TMath::Prob(ip3dpv.second.value()/ip3dpv.second.error(),1));
346     //printf("gsf reduced chisquared = %5f, probability = %5f\n", pvGsfCompat.second/2, TMath::Prob(pvGsfCompat.second,2));
347 bendavid 1.43 }
348    
349     }
350    
351     //fill conversion partner track info
352     ConversionFinder convFinder;
353     ConversionInfo convInfo = convFinder.getConversionInfo(*iM, hGeneralTracks, bfield);
354     outElectron->SetConvPartnerDCotTheta(convInfo.dcot());
355     outElectron->SetConvPartnerDist(convInfo.dist());
356     outElectron->SetConvPartnerRadius(convInfo.radiusOfConversion());
357    
358 loizides 1.29 // fill Electron ID information
359 sixie 1.13 outElectron->SetPassLooseID((*eidLooseMap)[eRef]);
360     outElectron->SetPassTightID((*eidTightMap)[eRef]);
361 bendavid 1.31
362     if (verbose_>1) {
363     double recomass = sqrt(iM->energy()*iM->energy() - iM->p()*iM->p());
364 loizides 1.32 printf(" mithep::Electron, pt=%5f, eta=%5f, phi=%5f, energy=%5f, p=%5f, mass=%5f\n",
365     outElectron->Pt(), outElectron->Eta(), outElectron->Phi(),
366     outElectron->E(), outElectron->P(), outElectron->Mass());
367     printf("reco::GsfElectron , pt=%5f, eta=%5f, phi=%5f, energy=%5f, p=%5f, mass=%5f\n",
368     iM->pt(), iM->eta(), iM->phi(), iM->energy(), iM->p(), recomass);
369 bendavid 1.31 }
370 sixie 1.12 }
371 loizides 1.1 electrons_->Trim();
372     }