ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitProd/TreeFiller/src/FillerElectrons.cc
Revision: 1.66
Committed: Sat May 12 15:55:11 2012 UTC (12 years, 11 months ago) by paus
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_028a, Mit_028
Changes since 1.65: +3 -1 lines
Log Message:
New tag 028 (EcalEnergy added).

File Contents

# User Rev Content
1 paus 1.66 // $Id: FillerElectrons.cc,v 1.65 2012/05/05 16:49:59 paus Exp $
2 paus 1.65
3 pharris 1.59 #include "FWCore/MessageLogger/interface/MessageLogger.h"
4 loizides 1.1 #include "MitProd/TreeFiller/interface/FillerElectrons.h"
5     #include "DataFormats/TrackReco/interface/Track.h"
6     #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
7     #include "DataFormats/TrackReco/interface/TrackFwd.h"
8 sixie 1.11 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
9     #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
10     #include "DataFormats/EgammaReco/interface/ClusterShape.h"
11     #include "DataFormats/EgammaReco/interface/BasicClusterShapeAssociation.h"
12 bendavid 1.18 #include "DataFormats/Common/interface/RefToPtr.h"
13 loizides 1.35 #include "DataFormats/Common/interface/ValueMap.h"
14     #include "AnalysisDataFormats/Egamma/interface/ElectronID.h"
15     #include "AnalysisDataFormats/Egamma/interface/ElectronIDAssociation.h"
16 sixie 1.11 #include "RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h"
17     #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaEcalIsolation.h"
18 sixie 1.14 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTowerIsolation.h"
19 loizides 1.35 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
20 bendavid 1.43 #include "DataFormats/VertexReco/interface/VertexFwd.h"
21     #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
22     #include "TrackingTools/TransientTrack/plugins/TransientTrackBuilderESProducer.h"
23 bendavid 1.44 #include "RecoVertex/GaussianSumVertexFit/interface/GsfVertexTrackCompatibilityEstimator.h"
24 bendavid 1.43 #include "TrackingTools/IPTools/interface/IPTools.h"
25     #include "RecoEgamma/EgammaTools/interface/ConversionFinder.h"
26     #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
27     #include "MagneticField/Engine/interface/MagneticField.h"
28 loizides 1.35 #include "MitAna/DataTree/interface/ElectronCol.h"
29     #include "MitAna/DataTree/interface/Names.h"
30     #include "MitAna/DataTree/interface/Track.h"
31 loizides 1.33 #include "MitEdm/DataFormats/interface/RefToBaseToPtr.h"
32 loizides 1.35 #include "MitProd/ObjectService/interface/ObjectService.h"
33 bendavid 1.45 #include "MitEdm/DataFormats/interface/DecayPart.h"
34     #include "MitEdm/ConversionRejection/interface/ConversionMatcher.h"
35 bendavid 1.51 #include "RecoVertex/VertexTools/interface/LinearizedTrackStateFactory.h"
36     #include "RecoVertex/VertexTools/interface/VertexTrackFactory.h"
37     #include "RecoVertex/VertexPrimitives/interface/VertexTrack.h"
38     #include "RecoVertex/VertexPrimitives/interface/CachingVertex.h"
39     #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexUpdator.h"
40 ksung 1.57 #include "MitEdm/Tools/interface/VertexReProducer.h"
41 paus 1.65 #include "RecoEgamma/EgammaTools/interface/ggPFClusters.h"
42     #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
43     #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
44     #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
45     #include "Geometry/Records/interface/CaloGeometryRecord.h"
46     #include "FWCore/Framework/interface/ESHandle.h"
47     #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
48     #include "DataFormats/Math/interface/deltaR.h"
49 sixie 1.13
50 loizides 1.1 using namespace std;
51     using namespace edm;
52     using namespace mithep;
53    
54 loizides 1.5 //--------------------------------------------------------------------------------------------------
55 loizides 1.28 FillerElectrons::FillerElectrons(const edm::ParameterSet &cfg, const char *name, bool active) :
56 paus 1.65 BaseFiller (cfg,name,active),
57     edmName_ (Conf().getUntrackedParameter<string>("edmName","pixelMatchGsfElectrons")),
58     expectedHitsName_ (Conf().getUntrackedParameter<string>("expectedHitsName","")),
59     mitName_ (Conf().getUntrackedParameter<string>("mitName",Names::gkElectronBrn)),
60     gsfTrackMapName_ (Conf().getUntrackedParameter<string>("gsfTrackMapName","")),
61     trackerTrackMapName_ (Conf().getUntrackedParameter<string>("trackerTrackMapName","")),
62 sixie 1.12 barrelSuperClusterMapName_(Conf().getUntrackedParameter<string>("barrelSuperClusterMapName","")),
63     endcapSuperClusterMapName_(Conf().getUntrackedParameter<string>("endcapSuperClusterMapName","")),
64 paus 1.65 checkClusterActive_ (Conf().getUntrackedParameter<bool>("requireClusterAndGsfMap",true)),
65     pfSuperClusterMapName_ (Conf().getUntrackedParameter<string>("pfSuperClusterMapName","")),
66     pfClusterMapName_ (Conf().getUntrackedParameter<string>("pfClusterMapName","")),
67     electronMapName_ (Conf().getUntrackedParameter<string>("electronMapName","")),
68     eIDCutBasedTightName_ (Conf().getUntrackedParameter<string>("eIDCutBasedTightName","eidTight")),
69     eIDCutBasedLooseName_ (Conf().getUntrackedParameter<string>("eIDCutBasedLooseName","eidLoose")),
70     eIDLikelihoodName_ (Conf().getUntrackedParameter<string>("eIDLikelihoodName","")),
71     pvEdmName_ (Conf().getUntrackedParameter<string>("pvEdmName","offlinePrimaryVertices")),
72     pvBSEdmName_ (Conf().getUntrackedParameter<string>("pvBSEdmName","offlinePrimaryVerticesWithBS")),
73     recomputeConversionInfo_ (Conf().getUntrackedParameter<bool>("recomputeConversionInfo",false)),
74     fitUnbiasedVertex_ (Conf().getUntrackedParameter<bool>("fitUnbiasedVertex",true)),
75     electronMap_ (new mithep::ElectronMap),
76     electrons_ (new mithep::ElectronArr(16)),
77     gsfTrackMap_ (0),
78     trackerTrackMap_ (0),
79     barrelSuperClusterMap_ (0),
80     endcapSuperClusterMap_ (0),
81     pfSuperClusterMap_ (0),
82     pfClusterMap_ (0)
83 loizides 1.1 {
84     // Constructor.
85     }
86    
87 loizides 1.5 //--------------------------------------------------------------------------------------------------
88 loizides 1.1 FillerElectrons::~FillerElectrons()
89     {
90     // Destructor.
91 loizides 1.6
92     delete electrons_;
93 pharris 1.62 delete electronMap_;
94 loizides 1.1 }
95    
96 loizides 1.5 //--------------------------------------------------------------------------------------------------
97 bendavid 1.42 void FillerElectrons::BookDataBlock(TreeWriter &tws)
98 loizides 1.1 {
99 loizides 1.10 // Add electron branch to our tree and get our maps.
100 loizides 1.1
101 loizides 1.29 tws.AddBranch(mitName_,&electrons_);
102     OS()->add<mithep::ElectronArr>(electrons_,mitName_);
103     if (!gsfTrackMapName_.empty()) {
104     gsfTrackMap_ = OS()->get<TrackMap>(gsfTrackMapName_);
105     if (gsfTrackMap_)
106     AddBranchDep(mitName_,gsfTrackMap_->GetBrName());
107     }
108     if (!trackerTrackMapName_.empty()) {
109     trackerTrackMap_ = OS()->get<TrackMap>(trackerTrackMapName_);
110     if (trackerTrackMap_)
111     AddBranchDep(mitName_,trackerTrackMap_->GetBrName());
112     }
113     if (!barrelSuperClusterMapName_.empty()) {
114     barrelSuperClusterMap_ = OS()->get<SuperClusterMap>(barrelSuperClusterMapName_);
115     if (barrelSuperClusterMap_)
116     AddBranchDep(mitName_,barrelSuperClusterMap_->GetBrName());
117     }
118     if (!endcapSuperClusterMapName_.empty()) {
119     endcapSuperClusterMap_ = OS()->get<SuperClusterMap>(endcapSuperClusterMapName_);
120     if (endcapSuperClusterMap_)
121     AddBranchDep(mitName_,endcapSuperClusterMap_->GetBrName());
122     }
123 bendavid 1.36 if (!pfSuperClusterMapName_.empty()) {
124     pfSuperClusterMap_ = OS()->get<SuperClusterMap>(pfSuperClusterMapName_);
125     if (pfSuperClusterMap_)
126     AddBranchDep(mitName_,pfSuperClusterMap_->GetBrName());
127     }
128 paus 1.65 if (!pfClusterMapName_.empty()) {
129     pfClusterMap_ = OS()->get<BasicClusterMap>(pfClusterMapName_);
130     if (pfClusterMap_)
131     AddBranchDep(mitName_,pfClusterMap_->GetBrName());
132     }
133 pharris 1.62 if (!electronMapName_.empty()) {
134     electronMap_->SetBrName(mitName_);
135     OS()->add<ElectronMap>(electronMap_,electronMapName_);
136     }
137 loizides 1.1 }
138    
139 loizides 1.5 //--------------------------------------------------------------------------------------------------
140 loizides 1.29 void FillerElectrons::FillDataBlock(const edm::Event &event, const edm::EventSetup &setup)
141 loizides 1.1 {
142     // Fill electrons from edm collection into our collection.
143    
144 pharris 1.62 electrons_ ->Delete();
145     electronMap_->Reset();
146    
147 sixie 1.11 Handle<reco::GsfElectronCollection> hElectronProduct;
148 loizides 1.7 GetProduct(edmName_, hElectronProduct, event);
149 loizides 1.1
150 loizides 1.29 // handles to get the electron ID information
151 peveraer 1.47 Handle<edm::ValueMap<float> > eidLooseMap;
152     GetProduct(eIDCutBasedLooseName_, eidLooseMap, event);
153     Handle<edm::ValueMap<float> > eidTightMap;
154 peveraer 1.48 GetProduct(eIDCutBasedTightName_, eidTightMap, event);
155 bendavid 1.49 edm::Handle<edm::ValueMap<float> > eidLikelihoodMap;
156     if (!eIDLikelihoodName_.empty()) {
157     GetProduct(eIDLikelihoodName_, eidLikelihoodMap, event);
158     }
159    
160 bendavid 1.43 edm::Handle<reco::VertexCollection> hVertex;
161     event.getByLabel(pvEdmName_, hVertex);
162     const reco::VertexCollection *pvCol = hVertex.product();
163     edm::Handle<reco::VertexCollection> hVertexBS;
164     event.getByLabel(pvBSEdmName_, hVertexBS);
165     const reco::VertexCollection *pvBSCol = hVertexBS.product();
166     edm::Handle<reco::TrackCollection> hGeneralTracks;
167     event.getByLabel("generalTracks", hGeneralTracks);
168     //const reco::VertexCollection *trackCol = hGeneralTracks.product();
169 bendavid 1.45
170 bendavid 1.54 edm::Handle<reco::GsfTrackCollection> hGsfTracks;
171     event.getByLabel("electronGsfTracks", hGsfTracks);
172    
173 bendavid 1.45 edm::Handle<std::vector<mitedm::DecayPart> > hConversions;
174     event.getByLabel("mvfConversionRemoval", hConversions);
175    
176     mitedm::ConversionMatcher convMatcher;
177 bendavid 1.50
178 bendavid 1.43 edm::ESHandle<TransientTrackBuilder> hTransientTrackBuilder;
179     setup.get<TransientTrackRecord>().get("TransientTrackBuilder",hTransientTrackBuilder);
180     const TransientTrackBuilder *transientTrackBuilder = hTransientTrackBuilder.product();
181 bendavid 1.19
182 bendavid 1.44 GsfVertexTrackCompatibilityEstimator gsfEstimator;
183    
184 bendavid 1.51 LinearizedTrackStateFactory lTrackFactory;
185     VertexTrackFactory<5> vTrackFactory;
186     KalmanVertexUpdator<5> updator;
187    
188 paus 1.65 //pf photon stuff
189     edm::Handle< EcalRecHitCollection > pEBRecHits;
190     event.getByLabel("reducedEcalRecHitsEB", pEBRecHits );
191     edm::Handle< EcalRecHitCollection > pEERecHits;
192     event.getByLabel( "reducedEcalRecHitsEE", pEERecHits );
193    
194     edm::ESHandle<CaloGeometry> pGeometry;
195     setup.get<CaloGeometryRecord>().get(pGeometry);
196    
197     const CaloSubdetectorGeometry *geometryEB = pGeometry->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
198     const CaloSubdetectorGeometry *geometryEE = pGeometry->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
199    
200     ggPFClusters pfclusters(pEBRecHits, pEERecHits, geometryEB, geometryEE);
201    
202 pharris 1.61 //Get Magnetic Field from event setup, taking value at (0,0,0)
203 bendavid 1.54 edm::ESHandle<MagneticField> magneticField;
204     setup.get<IdealMagneticFieldRecord>().get(magneticField);
205     const double bfield = magneticField->inTesla(GlobalPoint(0.,0.,0.)).z();
206 bendavid 1.43
207 sixie 1.11 const reco::GsfElectronCollection inElectrons = *(hElectronProduct.product());
208 loizides 1.29 // loop over electrons
209 sixie 1.11 for (reco::GsfElectronCollection::const_iterator iM = inElectrons.begin();
210 loizides 1.1 iM != inElectrons.end(); ++iM) {
211 sixie 1.13
212 loizides 1.29 // the index and Ref are needed for the eID association Map
213 sixie 1.13 unsigned int iElectron = iM - inElectrons.begin();
214     reco::GsfElectronRef eRef(hElectronProduct, iElectron);
215    
216 loizides 1.8 mithep::Electron *outElectron = electrons_->AddNew();
217 bendavid 1.31
218     outElectron->SetPtEtaPhi(iM->pt(),iM->eta(),iM->phi());
219 sixie 1.12
220 bendavid 1.40 outElectron->SetCharge(iM->charge());
221     outElectron->SetScPixCharge(iM->scPixCharge());
222 pharris 1.62 outElectron->SetESuperClusterOverP(iM->eSuperClusterOverP());
223 loizides 1.30 outElectron->SetESeedClusterOverPout(iM->eSeedClusterOverPout());
224 sixie 1.63 outElectron->SetEEleClusterOverPout(iM->eEleClusterOverPout());
225 loizides 1.30 outElectron->SetPIn(iM->trackMomentumAtVtx().R());
226 loizides 1.29 outElectron->SetPOut(iM->trackMomentumOut().R());
227 loizides 1.30 outElectron->SetDeltaEtaSuperClusterTrackAtVtx(iM->deltaEtaSuperClusterTrackAtVtx());
228     outElectron->SetDeltaEtaSeedClusterTrackAtCalo(iM->deltaEtaSeedClusterTrackAtCalo());
229     outElectron->SetDeltaPhiSuperClusterTrackAtVtx(iM->deltaPhiSuperClusterTrackAtVtx());
230     outElectron->SetDeltaPhiSeedClusterTrackAtCalo(iM->deltaPhiSeedClusterTrackAtCalo());
231     outElectron->SetIsEnergyScaleCorrected(iM->isEnergyScaleCorrected());
232 bendavid 1.55 //outElectron->SetIsMomentumCorrected(iM->isMomentumCorrected());
233 bendavid 1.36 outElectron->SetNumberOfClusters(iM->basicClustersSize());
234 loizides 1.30 outElectron->SetClassification(iM->classification());
235 bendavid 1.37 outElectron->SetFBrem(iM->fbrem());
236 paus 1.66 outElectron->SetEcalEnergy(iM->correctedEcalEnergy());
237     outElectron->SetEcalEnergyError(iM->correctedEcalEnergyError());
238 bendavid 1.37
239 loizides 1.38 // pflow electron stuff
240 bendavid 1.41 outElectron->SetIsEcalDriven(iM->ecalDrivenSeed());
241     outElectron->SetIsTrackerDriven(iM->trackerDrivenSeed());
242 bendavid 1.37 outElectron->SetMva(iM->mva());
243    
244 loizides 1.29 // shower shape variables
245 bendavid 1.37 outElectron->SetE15(iM->e1x5());
246     outElectron->SetE25Max(iM->e2x5Max());
247     outElectron->SetE55(iM->e5x5());
248     outElectron->SetCovEtaEta(iM->sigmaEtaEta());
249     outElectron->SetCoviEtaiEta(iM->sigmaIetaIeta());
250     outElectron->SetHadronicOverEm(iM->hcalOverEcal());
251     outElectron->SetHcalDepth1OverEcal(iM->hcalDepth1OverEcal());
252     outElectron->SetHcalDepth2OverEcal(iM->hcalDepth2OverEcal());
253 bendavid 1.64 outElectron->SetHadOverEmTow(iM->hcalOverEcalBc());
254    
255    
256 loizides 1.38 // fill isolation variables for both cone sizes
257 bendavid 1.37 outElectron->SetEcalRecHitIsoDr04(iM->dr04EcalRecHitSumEt());
258     outElectron->SetHcalDepth1TowerSumEtDr04(iM->dr04HcalDepth1TowerSumEt());
259     outElectron->SetHcalDepth2TowerSumEtDr04(iM->dr04HcalDepth2TowerSumEt());
260     outElectron->SetTrackIsolationDr04(iM->dr04TkSumPt());
261 bendavid 1.64 outElectron->SetHCalIsoTowDr04( iM->dr04HcalTowerSumEtBc() );
262 bendavid 1.37 outElectron->SetEcalRecHitIsoDr03(iM->dr03EcalRecHitSumEt());
263     outElectron->SetHcalTowerSumEtDr03(iM->dr03HcalTowerSumEt());
264     outElectron->SetHcalDepth1TowerSumEtDr03(iM->dr03HcalDepth1TowerSumEt());
265     outElectron->SetHcalDepth2TowerSumEtDr03(iM->dr03HcalDepth2TowerSumEt());
266     outElectron->SetTrackIsolationDr03(iM->dr03TkSumPt());
267 bendavid 1.64 outElectron->SetHCalIsoTowDr03( iM->dr03HcalTowerSumEtBc() );
268    
269    
270 bendavid 1.56 //pflow isolation
271     outElectron->SetPFChargedHadronIso(iM->pfIsolationVariables().chargedHadronIso);
272     outElectron->SetPFChargedHadronIso(iM->pfIsolationVariables().neutralHadronIso);
273     outElectron->SetPFChargedHadronIso(iM->pfIsolationVariables().photonIso);
274    
275 loizides 1.38 // fiducial flags
276 bendavid 1.37 outElectron->SetIsEB(iM->isEB());
277     outElectron->SetIsEE(iM->isEE());
278     outElectron->SetIsEBEEGap(iM->isEBEEGap());
279     outElectron->SetIsEBEtaGap(iM->isEBEtaGap());
280     outElectron->SetIsEBPhiGap(iM->isEBPhiGap());
281     outElectron->SetIsEEDeeGap(iM->isEEDeeGap());
282     outElectron->SetIsEERingGap(iM->isEERingGap());
283 sixie 1.12
284 loizides 1.38 // gsf-tracker match quality
285 bendavid 1.37 outElectron->SetFracSharedHits(iM->shFracInnerHits());
286 sixie 1.15
287 loizides 1.29 // make proper links to Tracks and Super Clusters
288 bendavid 1.49 if (gsfTrackMap_ && iM->gsfTrack().isNonnull()) {
289 pharris 1.59 try{outElectron->SetGsfTrk(gsfTrackMap_->GetMit(refToPtr(iM->gsfTrack())));}
290     catch(...) {
291     if(checkClusterActive_) throw edm::Exception(edm::errors::Configuration, "FillerElectrons:FillDataBlock()\n")
292     << "Error! GSF track unmapped collection " << edmName_ << endl;
293     }
294 bendavid 1.49 }
295     // make links to ambigous gsf tracks
296     if (gsfTrackMap_) {
297     for (reco::GsfTrackRefVector::const_iterator agsfi = iM->ambiguousGsfTracksBegin(); agsfi != iM->ambiguousGsfTracksEnd(); ++agsfi) {
298 pharris 1.59 try{outElectron->AddAmbiguousGsfTrack(gsfTrackMap_->GetMit(refToPtr(*agsfi)));}
299     catch(...) {
300     if(checkClusterActive_) throw edm::Exception(edm::errors::Configuration, "FillerElectrons:FillDataBlock()\n")
301     << "Error! GSF track unmapped collection " << edmName_ << endl;
302     }
303 bendavid 1.49 }
304     }
305    
306     // make tracker track links,
307 bendavid 1.37 if (trackerTrackMap_ && iM->closestCtfTrackRef().isNonnull()) {
308 bendavid 1.36 outElectron->SetTrackerTrk(trackerTrackMap_->GetMit(refToPtr(iM->closestCtfTrackRef())));
309 bendavid 1.19 }
310 pharris 1.61
311 loizides 1.38 if (barrelSuperClusterMap_ && endcapSuperClusterMap_ &&
312     pfSuperClusterMap_ && iM->superCluster().isNonnull()) {
313 bendavid 1.37 if(barrelSuperClusterMap_->HasMit(iM->superCluster())) {
314 sixie 1.12 outElectron->SetSuperCluster(barrelSuperClusterMap_->GetMit(iM->superCluster()));
315 bendavid 1.36 }
316     else if (endcapSuperClusterMap_->HasMit(iM->superCluster())) {
317 sixie 1.12 outElectron->SetSuperCluster(endcapSuperClusterMap_->GetMit(iM->superCluster()));
318     }
319 bendavid 1.36 else if (pfSuperClusterMap_->HasMit(iM->superCluster())) {
320     outElectron->SetSuperCluster(pfSuperClusterMap_->GetMit(iM->superCluster()));
321     }
322 pharris 1.59 else if(checkClusterActive_) {
323     throw edm::Exception(edm::errors::Configuration, "FillerElectrons:FillDataBlock()\n")
324 bendavid 1.36 << "Error! SuperCluster reference in unmapped collection " << edmName_ << endl;
325 pharris 1.59 }
326 bendavid 1.36 }
327 paus 1.65
328     if (pfSuperClusterMap_ && iM->pflowSuperCluster().isNonnull()) {
329     outElectron->SetPFSuperCluster(pfSuperClusterMap_->GetMit(iM->pflowSuperCluster()));
330     }
331    
332     //find matching egamma supercluster first by ref, or by geometric matching if only pfsupercluster is linked
333     const reco::SuperCluster *egammasc = 0;
334     double mindr = 999.;
335     if (iM->superCluster().isNonnull() && iM->superCluster() != iM->pflowSuperCluster()) {
336     egammasc = iM->superCluster().get();
337     }
338     else {
339     for (SuperClusterMap::fwdMapType::const_iterator scit = barrelSuperClusterMap_->FwdMap().begin(); scit != barrelSuperClusterMap_->FwdMap().end(); ++scit) {
340     double dr = reco::deltaR(*iM->pflowSuperCluster(),*scit->first);
341     if (dr<0.1 && dr<mindr) {
342     egammasc = scit->first.get();
343     mindr = dr;
344     }
345     }
346     for (SuperClusterMap::fwdMapType::const_iterator scit = endcapSuperClusterMap_->FwdMap().begin(); scit != endcapSuperClusterMap_->FwdMap().end(); ++scit) {
347     double dr = reco::deltaR(*iM->pflowSuperCluster(),*scit->first);
348     if (dr<0.1 && dr<mindr) {
349     egammasc = scit->first.get();
350     mindr = dr;
351     }
352     }
353     }
354    
355     //tag overlapping energy of pflow clusters
356     if (pfClusterMap_ && iM->pflowSuperCluster().isNonnull() && egammasc) {
357     for (reco::CaloCluster_iterator pfcit = iM->pflowSuperCluster()->clustersBegin();
358     pfcit!=iM->pflowSuperCluster()->clustersEnd(); ++pfcit) {
359     float eoverlap = pfclusters.getPFSuperclusterOverlap( **pfcit, *egammasc );
360     const_cast<mithep::BasicCluster*>(pfClusterMap_->GetMit(*pfcit))->SetMatchedEnergy(eoverlap);
361     }
362     }
363 pharris 1.61
364 sixie 1.63 //compute NLayersWithMeasurement for associated ctf track
365     if (iM->closestCtfTrackRef().isNonnull()) {
366     outElectron->SetCTFTrkNLayersWithMeasurement(iM->closestCtfTrackRef()->hitPattern().trackerLayersWithMeasurement());
367     } else {
368     outElectron->SetCTFTrkNLayersWithMeasurement(-1);
369     }
370    
371 bendavid 1.43 //compute impact parameter with respect to PV
372     if (iM->gsfTrack().isNonnull()) {
373     const reco::TransientTrack &tt = transientTrackBuilder->build(iM->gsfTrack());
374 bendavid 1.51
375     reco::TransientTrack ttckf;
376    
377 bendavid 1.58 reco::Vertex thevtx = pvCol->at(0);
378     reco::Vertex thevtxbs = pvBSCol->at(0);
379 bendavid 1.52
380 bendavid 1.58 reco::Vertex thevtxub = pvCol->at(0);
381     reco::Vertex thevtxubbs = pvBSCol->at(0);
382 pharris 1.61
383 bendavid 1.51 //check if closest ctf track is included in PV and if so, remove it before computing impact parameters and uncertainties
384     if (iM->closestCtfTrackRef().isNonnull()) {
385     ttckf = transientTrackBuilder->build(iM->closestCtfTrackRef());
386 ksung 1.57
387     //check if closest ctf track is included in PV and if so, remove it from the collection of
388     //tracks associated with the PV and perform a refit before computing impact parameters and uncertainties
389     reco::TrackCollection newTkCollection;
390     bool foundMatch = false;
391     for(reco::Vertex::trackRef_iterator itk = thevtx.tracks_begin(); itk!=thevtx.tracks_end(); itk++) {
392     if(iM->closestCtfTrack().ctfTrack.isNonnull()) {
393     bool refMatching = (itk->get() == &*(iM->closestCtfTrack().ctfTrack));
394     float shFraction = iM->closestCtfTrack().shFracInnerHits;
395     if(refMatching && shFraction > 0.5) {
396 pharris 1.61 foundMatch = true;
397     continue;
398 ksung 1.57 }
399     }
400     newTkCollection.push_back(*itk->get());
401     }
402 bendavid 1.51
403 pharris 1.61 if(foundMatch && fitUnbiasedVertex_) {
404 ksung 1.57 edm::Handle<reco::BeamSpot> bs;
405     event.getByLabel(edm::InputTag("offlineBeamSpot"),bs);
406 pharris 1.61
407     VertexReProducer revertex(hVertex,event); //Needs to be fixed
408    
409     edm::Handle<reco::BeamSpot> pvbeamspot;
410 ksung 1.57 event.getByLabel(revertex.inputBeamSpot(),pvbeamspot);
411 pharris 1.61 vector<TransientVertex> pvs = revertex.makeVertices(newTkCollection,*pvbeamspot,setup);
412     if(pvs.size()>0) {
413 bendavid 1.58 thevtxub = pvs.front(); // take the first in the list
414     }
415 ksung 1.57 VertexReProducer revertexbs(hVertexBS,event);
416     edm::Handle<reco::BeamSpot> pvbsbeamspot;
417     event.getByLabel(revertexbs.inputBeamSpot(),pvbsbeamspot);
418     vector<TransientVertex> pvbss = revertexbs.makeVertices(newTkCollection,*pvbsbeamspot,setup);
419 pharris 1.61 if(pvbss.size()>0) {
420 bendavid 1.58 thevtxubbs = pvbss.front(); // take the first in the list
421     }
422 ksung 1.57 }
423 bendavid 1.51 }
424     //preserve sign of transverse impact parameter (cross-product definition from track, not lifetime-signing)
425     const double gsfsign = ( (-iM->gsfTrack()->dxy(thevtx.position())) >=0 ) ? 1. : -1.;
426     const double gsfsignbs = ( (-iM->gsfTrack()->dxy(thevtxbs.position())) >=0 ) ? 1. : -1.;
427     const std::pair<bool,Measurement1D> &d0pv = IPTools::absoluteTransverseImpactParameter(tt,thevtx);
428 bendavid 1.43 if (d0pv.first) {
429 bendavid 1.51 outElectron->SetD0PV(gsfsign*d0pv.second.value());
430 bendavid 1.43 outElectron->SetD0PVErr(d0pv.second.error());
431     }
432 bendavid 1.44 else {
433 bendavid 1.51 outElectron->SetD0PV(-999.0);
434 bendavid 1.44 }
435    
436 bendavid 1.51 const std::pair<bool,Measurement1D> &ip3dpv = IPTools::absoluteImpactParameter3D(tt,thevtx);
437 bendavid 1.43 if (ip3dpv.first) {
438 bendavid 1.51 outElectron->SetIp3dPV(gsfsign*ip3dpv.second.value());
439 bendavid 1.43 outElectron->SetIp3dPVErr(ip3dpv.second.error());
440     }
441 bendavid 1.44 else {
442 bendavid 1.51 outElectron->SetIp3dPV(-999.0);
443 bendavid 1.44 }
444 bendavid 1.43
445 bendavid 1.51 const std::pair<bool,Measurement1D> &d0pvbs = IPTools::absoluteTransverseImpactParameter(tt,thevtxbs);
446 bendavid 1.43 if (d0pvbs.first) {
447 bendavid 1.51 outElectron->SetD0PVBS(gsfsignbs*d0pvbs.second.value());
448 bendavid 1.43 outElectron->SetD0PVBSErr(d0pvbs.second.error());
449     }
450 bendavid 1.44 else {
451 bendavid 1.51 outElectron->SetD0PVBS(-999.0);
452 bendavid 1.44 }
453 bendavid 1.43
454 bendavid 1.51 const std::pair<bool,Measurement1D> &ip3dpvbs = IPTools::absoluteImpactParameter3D(tt,thevtxbs);
455 bendavid 1.43 if (ip3dpvbs.first) {
456 bendavid 1.51 outElectron->SetIp3dPVBS(gsfsignbs*ip3dpvbs.second.value());
457 bendavid 1.43 outElectron->SetIp3dPVBSErr(ip3dpvbs.second.error());
458     }
459 bendavid 1.44 else {
460 bendavid 1.51 outElectron->SetIp3dPVBS(-999.0);
461 bendavid 1.44 }
462    
463 bendavid 1.58 const std::pair<bool,Measurement1D> &d0pvub = IPTools::absoluteTransverseImpactParameter(tt,thevtxub);
464     if (d0pvub.first) {
465     outElectron->SetD0PVUB(gsfsign*d0pvub.second.value());
466     outElectron->SetD0PVUBErr(d0pvub.second.error());
467     }
468     else {
469     outElectron->SetD0PVUB(-999.0);
470     }
471    
472    
473     const std::pair<bool,Measurement1D> &ip3dpvub = IPTools::absoluteImpactParameter3D(tt,thevtxub);
474     if (ip3dpvub.first) {
475     outElectron->SetIp3dPVUB(gsfsign*ip3dpvub.second.value());
476     outElectron->SetIp3dPVUBErr(ip3dpvub.second.error());
477     }
478     else {
479     outElectron->SetIp3dPVUB(-999.0);
480     }
481    
482     const std::pair<bool,Measurement1D> &d0pvubbs = IPTools::absoluteTransverseImpactParameter(tt,thevtxubbs);
483     if (d0pvubbs.first) {
484     outElectron->SetD0PVUBBS(gsfsignbs*d0pvubbs.second.value());
485     outElectron->SetD0PVUBBSErr(d0pvubbs.second.error());
486     }
487     else {
488     outElectron->SetD0PVUBBS(-999.0);
489     }
490    
491     const std::pair<bool,Measurement1D> &ip3dpvubbs = IPTools::absoluteImpactParameter3D(tt,thevtxubbs);
492     if (ip3dpvubbs.first) {
493     outElectron->SetIp3dPVUBBS(gsfsignbs*ip3dpvubbs.second.value());
494     outElectron->SetIp3dPVUBBSErr(ip3dpvubbs.second.error());
495     }
496     else {
497     outElectron->SetIp3dPVUBBS(-999.0);
498     }
499    
500 bendavid 1.51 if (iM->closestCtfTrackRef().isNonnull()) {
501    
502     const double ckfsign = ( (-iM->closestCtfTrackRef()->dxy(thevtx.position())) >=0 ) ? 1. : -1.;
503     const double ckfsignbs = ( (-iM->closestCtfTrackRef()->dxy(thevtxbs.position())) >=0 ) ? 1. : -1.;
504    
505     const std::pair<bool,Measurement1D> &d0pvckf = IPTools::absoluteTransverseImpactParameter(ttckf,thevtx);
506     if (d0pvckf.first) {
507     outElectron->SetD0PVCkf(ckfsign*d0pvckf.second.value());
508     outElectron->SetD0PVCkfErr(d0pvckf.second.error());
509     }
510     else {
511     outElectron->SetD0PVCkf(-999.0);
512     }
513    
514    
515     const std::pair<bool,Measurement1D> &ip3dpvckf = IPTools::absoluteImpactParameter3D(ttckf,thevtx);
516     if (ip3dpvckf.first) {
517     outElectron->SetIp3dPVCkf(ckfsign*ip3dpvckf.second.value());
518     outElectron->SetIp3dPVCkfErr(ip3dpvckf.second.error());
519     }
520     else {
521     outElectron->SetIp3dPVCkf(-999.0);
522     }
523    
524     const std::pair<bool,Measurement1D> &d0pvbsckf = IPTools::absoluteTransverseImpactParameter(ttckf,thevtxbs);
525     if (d0pvbsckf.first) {
526     outElectron->SetD0PVBSCkf(ckfsignbs*d0pvbsckf.second.value());
527     outElectron->SetD0PVBSCkfErr(d0pvbsckf.second.error());
528     }
529     else {
530     outElectron->SetD0PVBSCkf(-999.0);
531     }
532    
533     const std::pair<bool,Measurement1D> &ip3dpvbsckf = IPTools::absoluteImpactParameter3D(ttckf,thevtxbs);
534     if (ip3dpvbsckf.first) {
535     outElectron->SetIp3dPVBSCkf(ckfsignbs*ip3dpvbsckf.second.value());
536     outElectron->SetIp3dPVBSCkfErr(ip3dpvbsckf.second.error());
537     }
538     else {
539     outElectron->SetIp3dPVBSCkf(-999.0);
540     }
541 bendavid 1.58
542     const std::pair<bool,Measurement1D> &d0pvubckf = IPTools::absoluteTransverseImpactParameter(ttckf,thevtxub);
543     if (d0pvubckf.first) {
544     outElectron->SetD0PVUBCkf(ckfsign*d0pvubckf.second.value());
545     outElectron->SetD0PVUBCkfErr(d0pvubckf.second.error());
546     }
547     else {
548     outElectron->SetD0PVUBCkf(-999.0);
549     }
550    
551    
552     const std::pair<bool,Measurement1D> &ip3dpvubckf = IPTools::absoluteImpactParameter3D(ttckf,thevtxub);
553     if (ip3dpvubckf.first) {
554     outElectron->SetIp3dPVUBCkf(ckfsign*ip3dpvubckf.second.value());
555     outElectron->SetIp3dPVUBCkfErr(ip3dpvubckf.second.error());
556     }
557     else {
558     outElectron->SetIp3dPVUBCkf(-999.0);
559     }
560    
561     const std::pair<bool,Measurement1D> &d0pvubbsckf = IPTools::absoluteTransverseImpactParameter(ttckf,thevtxubbs);
562     if (d0pvubbsckf.first) {
563     outElectron->SetD0PVUBBSCkf(ckfsignbs*d0pvubbsckf.second.value());
564     outElectron->SetD0PVUBBSCkfErr(d0pvubbsckf.second.error());
565     }
566     else {
567     outElectron->SetD0PVUBBSCkf(-999.0);
568     }
569    
570     const std::pair<bool,Measurement1D> &ip3dpvubbsckf = IPTools::absoluteImpactParameter3D(ttckf,thevtxubbs);
571     if (ip3dpvubbsckf.first) {
572     outElectron->SetIp3dPVUBBSCkf(ckfsignbs*ip3dpvubbsckf.second.value());
573     outElectron->SetIp3dPVUBBSCkfErr(ip3dpvubbsckf.second.error());
574     }
575     else {
576     outElectron->SetIp3dPVUBBSCkf(-999.0);
577     }
578    
579 bendavid 1.51 }
580     else {
581     outElectron->SetD0PVCkf(-999.0);
582     outElectron->SetIp3dPVCkf(-999.0);
583     outElectron->SetD0PVBSCkf(-999.0);
584     outElectron->SetIp3dPVBSCkf(-999.0);
585 bendavid 1.58
586     outElectron->SetD0PVUBCkf(-999.0);
587     outElectron->SetIp3dPVUBCkf(-999.0);
588     outElectron->SetD0PVUBBSCkf(-999.0);
589     outElectron->SetIp3dPVUBBSCkf(-999.0);
590 bendavid 1.51 }
591    
592 bendavid 1.43 if (verbose_>1) {
593     printf("gsf track pt = %5f\n",iM->gsfTrack()->pt());
594     printf("gsf track mode pt = %5f\n",iM->gsfTrack()->ptMode());
595     printf("ttrack pt = %5f\n",tt.initialFreeState().momentum().perp());
596     //printf("ttrackgsf pt = %5f\n",ttgsf.innermostMeasurementState().globalMomentum().perp());
597 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));
598     //printf("gsf reduced chisquared = %5f, probability = %5f\n", pvGsfCompat.second/2, TMath::Prob(pvGsfCompat.second,2));
599 bendavid 1.43 }
600     }
601 bendavid 1.53
602 bendavid 1.43 //fill conversion partner track info
603 bendavid 1.54 if (recomputeConversionInfo_) {
604     ConversionFinder convFinder; outElectron->SetConvPartnerDCotTheta(iM->convDcot());
605     ConversionInfo convInfo = convFinder.getConversionInfo(*iM, hGeneralTracks, hGsfTracks, bfield);
606    
607     outElectron->SetConvFlag(convInfo.flag());
608     outElectron->SetConvPartnerDCotTheta(convInfo.dcot());
609     outElectron->SetConvPartnerDist(convInfo.dist());
610     outElectron->SetConvPartnerRadius(convInfo.radiusOfConversion());
611     reco::TrackRef ckfconvTrackRef = convInfo.conversionPartnerCtfTk();
612     reco::GsfTrackRef gsfconvTrackRef = convInfo.conversionPartnerGsfTk();
613    
614    
615     if ( gsfconvTrackRef.isNonnull() && gsfTrackMap_ ) {
616 pharris 1.59 try {outElectron->SetConvPartnerTrk(gsfTrackMap_->GetMit(edm::refToPtr(gsfconvTrackRef)));}
617     catch(...) {
618     if(checkClusterActive_) throw edm::Exception(edm::errors::Configuration, "FillerElectrons:FillDataBlock()\n")
619     << "Error! GSF track unmapped collection " << edmName_ << endl;
620     }
621    
622 bendavid 1.54 }
623     else if (ckfconvTrackRef.isNonnull() && trackerTrackMap_) {
624     outElectron->SetConvPartnerTrk(trackerTrackMap_->GetMit(edm::refToPtr(ckfconvTrackRef)));
625 bendavid 1.50 }
626 bendavid 1.54 }
627     else {
628     outElectron->SetConvFlag(iM->convFlags());
629     outElectron->SetConvPartnerDCotTheta(iM->convDcot());
630     outElectron->SetConvPartnerDist(iM->convDist());
631     outElectron->SetConvPartnerRadius(iM->convRadius());
632     reco::TrackBaseRef convTrackRef = iM->convPartner();
633     if (convTrackRef.isNonnull()) {
634     if ( dynamic_cast<const reco::GsfTrack*>(convTrackRef.get()) && gsfTrackMap_ ) {
635 pharris 1.59 try{outElectron->SetConvPartnerTrk(gsfTrackMap_->GetMit(mitedm::refToBaseToPtr(convTrackRef)));}
636     catch(...) {
637     if(checkClusterActive_) throw edm::Exception(edm::errors::Configuration, "FillerElectrons:FillDataBlock()\n")
638     << "Error! GSF track unmapped collection " << edmName_ << endl;
639     }
640 bendavid 1.54 }
641     else if (trackerTrackMap_) {
642     outElectron->SetConvPartnerTrk(trackerTrackMap_->GetMit(mitedm::refToBaseToPtr(convTrackRef)));
643     }
644 bendavid 1.50 }
645 bendavid 1.45 }
646 bendavid 1.50
647 loizides 1.29 // fill Electron ID information
648 peveraer 1.47 outElectron->SetPassLooseID((*eidLooseMap)[eRef]);
649     outElectron->SetPassTightID((*eidTightMap)[eRef]);
650 bendavid 1.49 if (!eIDLikelihoodName_.empty()) {
651     outElectron->SetIDLikelihood((*eidLikelihoodMap)[eRef]);
652     }
653     // fill corrected expected inner hits
654 bendavid 1.50 if(iM->gsfTrack().isNonnull()) {
655     outElectron->SetCorrectedNExpectedHitsInner(iM->gsfTrack()->trackerExpectedHitsInner().numberOfHits());
656 bendavid 1.49 }
657 bendavid 1.45 //fill additional conversion flag
658     outElectron->SetMatchesVertexConversion(convMatcher.matchesGoodConversion(*iM,hConversions));
659    
660 pharris 1.62 // add electron to map
661     edm::Ptr<reco::GsfElectron> thePtr(hElectronProduct, iM - inElectrons.begin());
662     electronMap_->Add(thePtr, outElectron);
663    
664 bendavid 1.31 if (verbose_>1) {
665     double recomass = sqrt(iM->energy()*iM->energy() - iM->p()*iM->p());
666 loizides 1.32 printf(" mithep::Electron, pt=%5f, eta=%5f, phi=%5f, energy=%5f, p=%5f, mass=%5f\n",
667     outElectron->Pt(), outElectron->Eta(), outElectron->Phi(),
668     outElectron->E(), outElectron->P(), outElectron->Mass());
669     printf("reco::GsfElectron , pt=%5f, eta=%5f, phi=%5f, energy=%5f, p=%5f, mass=%5f\n",
670     iM->pt(), iM->eta(), iM->phi(), iM->energy(), iM->p(), recomass);
671 bendavid 1.31 }
672 pharris 1.59 }
673 loizides 1.1 electrons_->Trim();
674     }