ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitProd/TreeFiller/src/FillerElectrons.cc
Revision: 1.69
Committed: Wed Jun 26 16:24:05 2013 UTC (11 years, 10 months ago) by pharris
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_032, Mit_031, HEAD
Changes since 1.68: +16 -4 lines
Error occurred while calculating annotation data.
Log Message:
fixed tracker track references in embedded

File Contents

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