ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitProd/TreeFiller/src/FillerElectrons.cc
Revision: 1.61.2.2
Committed: Wed May 23 03:24:07 2012 UTC (12 years, 11 months ago) by paus
Content type: text/plain
Branch: Mit_025c_branch
CVS Tags: Mit_025c_branch0
Changes since 1.61.2.1: +6 -1 lines
Log Message:
Backport.

File Contents

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