ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Morgan/src/ZeeVertexAnalyzer.cc
Revision: 1.3
Committed: Fri Oct 23 14:25:45 2009 UTC (15 years, 6 months ago) by lethuill
Content type: text/plain
Branch: MAIN
CVS Tags: all_3_3_2_01, all_3_2_5_02, HEAD
Changes since 1.2: +66 -4 lines
Error occurred while calculating annotation data.
Log Message:
Stock infos on two electron tracks in TRootBardak

File Contents

# Content
1 #include "../interface/ZeeVertexAnalyzer.h"
2
3 using namespace std;
4 using namespace reco;
5 using namespace edm;
6
7 ZeeVertexAnalyzer::ZeeVertexAnalyzer(const edm::ParameterSet& iConfig, const edm::ParameterSet& producersNames, int verbosity):verbosity_(verbosity), config_(iConfig)
8 {
9 allowMissingCollection_ = producersNames.getUntrackedParameter<bool>("allowMissingCollection", false);
10 primaryVertexProducer_ = producersNames.getParameter<edm::InputTag>("primaryVertexProducer");
11 beamSpotProducer_ = producersNames.getParameter<edm::InputTag>("beamSpotProducer");
12 trackProducer_ = producersNames.getParameter<edm::InputTag>("trackProducer");
13 electronProducer_ = producersNames.getParameter<edm::InputTag>("electronProducer");
14 }
15
16
17 ZeeVertexAnalyzer::~ZeeVertexAnalyzer()
18 {
19 }
20
21
22 bool ZeeVertexAnalyzer::getVertices(const edm::Event& iEvent, const edm::EventSetup& iSetup, TClonesArray* rootVertices, TRootBardak* rootBardak)
23 {
24 using namespace edm;
25 int iRootVertex = 0;
26
27 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
28 // Very simple Zee selection: at least 2 electrons pT>10 GeV, |eta|<2.5, 60 < Mee < 120
29 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
30
31 edm::Handle< reco::GsfElectronCollection > electrons;
32 iEvent.getByLabel(electronProducer_, electrons);
33 int nElectrons = electrons->size();
34
35 if(verbosity_>4) std::cout << std::endl << " Number of electrons = " << nElectrons << " (" << electronProducer_.label() << " producer)" << std::endl;
36 if(nElectrons<2) return false;
37
38 // Highest ET electron
39 Float_t maxPt = 0;
40 unsigned int iElectron1 = 999999;
41 for (unsigned int j=0; j<electrons->size(); j++)
42 {
43 const reco::GsfElectron* electron = & (electrons->at(j));
44 if ( electron->pt() > maxPt && abs(electron->eta()) < 2.5 )
45 {
46 maxPt = electron->pt();
47 iElectron1 = j;
48 }
49 }
50 if(iElectron1==999999) return false;
51
52
53 // Second highest ET electron
54 maxPt = 0;
55 unsigned int iElectron2 = 999999;
56 for (unsigned int j=0; j<electrons->size(); j++)
57 {
58 const reco::GsfElectron* electron = & (electrons->at(j));
59 if ( electron->pt() > maxPt && abs(electron->eta()) < 2.5 && j!=iElectron1 )
60 {
61 maxPt = electron->pt();
62 iElectron2 = j;
63 }
64 }
65 if(iElectron2==999999) return false;
66
67
68 // Z mass
69 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > p4Electron1= (electrons->at(iElectron1)).p4();
70 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > p4Electron2= (electrons->at(iElectron2)).p4();
71 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > p4Z = p4Electron1 + p4Electron2;
72
73 if(verbosity_>4) std::cout << " Electron 1 - (Et,eta,phi)=(" << electrons->at(iElectron1).pt() << "," << electrons->at(iElectron1).eta() << "," << electrons->at(iElectron1).phi() << ")" << std::endl;
74 if(verbosity_>4) std::cout << " Electron 2 - (Et,eta,phi)=(" << electrons->at(iElectron2).pt() << "," << electrons->at(iElectron2).eta() << "," << electrons->at(iElectron2).phi() << ")" << std::endl;
75 if(verbosity_>4) std::cout << " Mee=" << p4Z.mass() << std::endl;
76
77 //if( p4Z.mass()<60.0 || p4Z.mass()>120.0) return false;
78
79 // Find CTF track best matching the GSF track assoicated to the electron
80 // We will remove these tracks in the track collection for the primary vertex reconstruction
81 reco::TrackRef tk1 = (electrons->at(iElectron1)).closestCtfTrackRef();
82 reco::TrackRef tk2 = (electrons->at(iElectron2)).closestCtfTrackRef();
83
84
85 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
86 // Private Primary Vertex reconstruction with all CTF tracks
87 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
88
89 // Get BeamSpot
90 reco::BeamSpot vertexBeamSpot;
91 edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
92 iEvent.getByLabel(beamSpotProducer_, recoBeamSpotHandle);
93 if (recoBeamSpotHandle.isValid()) vertexBeamSpot = *recoBeamSpotHandle;
94
95 // Get Track Collection
96 edm::Handle<reco::TrackCollection> recoTracks;
97 iEvent.getByLabel(trackProducer_, recoTracks);
98
99 // Get Transient Track Collection
100 edm::ESHandle<TransientTrackBuilder> trackBuilder;
101 iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",trackBuilder);
102 vector<reco::TransientTrack> transientTracks = (*trackBuilder).build(recoTracks, vertexBeamSpot);
103
104 // Call Primary Vertex Reconstruction using all tracks in the event
105 PrimaryVertexProducerAlgorithm theVertexProducer(config_);
106 vector<TransientVertex> transientZeeVertices = theVertexProducer.vertices(transientTracks, vertexBeamSpot);
107
108 reco::VertexCollection zeeVertices;
109 for (vector<TransientVertex>::const_iterator vertex = transientZeeVertices.begin(); vertex != transientZeeVertices.end(); vertex++)
110 {
111 reco::Vertex v = *vertex;
112 zeeVertices.push_back(v);
113 }
114
115 if(verbosity_>4) std::cout << std::endl << " Private Primary Vertex reconstruction with all tracks - Number of vertices = " << zeeVertices.size() << std::endl;
116
117 for (unsigned int j=0; j<zeeVertices.size(); j++)
118 {
119 reco::Vertex* vertex = & (zeeVertices.at(j));
120
121 // Put your vertex selection here....
122 if (! vertex->isValid() ) continue;
123 if ( vertex->isFake() ) continue;
124
125 Int_t ntracks = 0;
126 Float_t higherPt = 0.;
127 Float_t scalarSumPt = 0.;
128 Float_t vectorSumPt = 0.;
129 math::XYZVector vectorSum(0.,0.,0.);
130
131 for( std::vector< reco::TrackBaseRef >::const_iterator it = vertex->tracks_begin(); it != vertex->tracks_end(); it++)
132 {
133 scalarSumPt += (**it).pt();
134 vectorSum += (**it).momentum();
135 if( (**it).pt()>higherPt ) higherPt=(**it).pt();
136 ntracks++;
137 }
138 vectorSumPt = sqrt(vectorSum.Perp2());
139
140 // No refitted tracks embeded in reco::Vertex....
141 //cout << "vertex->refittedTracks().size()=" << vertex->refittedTracks().size() << endl;
142
143 TRootVertex localVertex(
144 vertex->x()
145 ,vertex->y()
146 ,vertex->z()
147 ,vertex->xError()
148 ,vertex->yError()
149 ,vertex->zError()
150 );
151
152 localVertex.setAlgoName("AllTracks");
153 localVertex.setChi2( vertex->chi2() );
154 localVertex.setNdof( vertex->ndof() );
155 localVertex.setNtracks( ntracks );
156 localVertex.setHigherTrackPt( higherPt );
157 localVertex.setScalarSumPt( scalarSumPt );
158 localVertex.setVectorSumPt( vectorSumPt );
159
160 new( (*rootVertices)[iRootVertex] ) TRootVertex(localVertex);
161 if(verbosity_>2) cout << " ["<< setw(3) << iRootVertex << "] " << localVertex << endl;
162 iRootVertex++;
163 }
164
165
166 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
167 // Private Primary Vertex reconstruction with all CTF tracks except the 2 two tracks matching selected electrons
168 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
169
170 // Remove the two CTF tracks matching the two electron GSF tracks from the transientTracks collection
171 // Keep the two CTF electron tracks in a separate vector (eeTransientTracks)
172 unsigned int itk1 = 999999;
173 unsigned int itk2 = 999999;
174
175 for (unsigned int j=0; j<transientTracks.size(); j++)
176 {
177 const reco::Track* tk = & (transientTracks.at(j).track());
178 if(verbosity_>5) std::cout << " Tk[" << j << "] (Et,eta,phi)=(" << tk->pt() << "," << tk->eta() << "," << tk->phi() << ")" << std::endl;
179 if ( tk->pt()==tk1->pt() && tk->eta()==tk1->eta() && tk->phi()==tk1->phi() )
180 {
181 itk1 = j;
182 if(verbosity_>5) std::cout << " ===> this is tk1" << std::endl;
183 }
184 if ( tk->pt()==tk2->pt() && tk->eta()==tk2->eta() && tk->phi()==tk2->phi() )
185 {
186 itk2 = j;
187 if(verbosity_>5) std::cout << " ===> this is tk2" << std::endl;
188 }
189 }
190 if (itk1==999999 || itk2==999999)
191 {
192 std::cout << " ***** NO CTF TRACKS MATCHING THE ELECTRON *****" << std::endl;
193 return false;
194 }
195
196 if(verbosity_>4)
197 {
198 std::cout << " CTF Track 1 - (Et,eta,phi)=(" << transientTracks.at(itk1).track().pt() << "," << transientTracks.at(itk1).track().eta() << "," << transientTracks.at(itk1).track().phi() << " - Fraction of common hits with GSF =" << (electrons->at(iElectron1)).shFracInnerHits();
199 if (!transientTracks.at(itk1).stateAtBeamLine().isValid()) std::cout << " IPSig=-1";
200 else std::cout << " IPSig=" << transientTracks.at(itk1).stateAtBeamLine().transverseImpactParameter().significance();
201 std::cout << " normChi2=" << transientTracks.at(itk1).normalizedChi2() << " nSiHits=" << transientTracks.at(itk1).hitPattern().numberOfValidHits() << " nPxlHits=" << transientTracks.at(itk1).hitPattern().numberOfValidPixelHits() << std::endl;
202
203 std::cout << " CTF Track 2 - (Et,eta,phi)=(" << transientTracks.at(itk2).track().pt() << "," << transientTracks.at(itk2).track().eta() << "," << transientTracks.at(itk2).track().phi() << " - Fraction of common hits with GSF =" << (electrons->at(iElectron2)).shFracInnerHits();
204 if (!transientTracks.at(itk2).stateAtBeamLine().isValid()) std::cout << " IPSig=-1";
205 else std::cout << " IPSig=" << transientTracks.at(itk2).stateAtBeamLine().transverseImpactParameter().significance();
206 std::cout << " normChi2=" << transientTracks.at(itk2).normalizedChi2() << " nSiHits=" << transientTracks.at(itk2).hitPattern().numberOfValidHits() << " nPxlHits=" << transientTracks.at(itk2).hitPattern().numberOfValidPixelHits()<< std::endl;
207 }
208
209
210 rootBardak->setEle1(
211 TRootTrack(
212 transientTracks.at(itk1).track().px()
213 ,transientTracks.at(itk1).track().py()
214 ,transientTracks.at(itk1).track().pz()
215 ,transientTracks.at(itk1).track().p()
216 ,transientTracks.at(itk1).track().vx()
217 ,transientTracks.at(itk1).track().vy()
218 ,transientTracks.at(itk1).track().vz()
219 ,0
220 ,transientTracks.at(itk1).track().charge()
221 ,transientTracks.at(itk1).hitPattern().pixelLayersWithMeasurement()
222 ,transientTracks.at(itk1).hitPattern().stripLayersWithMeasurement()
223 ,transientTracks.at(itk1).track().chi2()
224 ,transientTracks.at(itk1).track().d0()
225 ,transientTracks.at(itk1).track().d0Error()
226 ,transientTracks.at(itk1).track().dz()
227 ,transientTracks.at(itk1).track().dzError()
228 )
229 ,transientTracks.at(itk1).hitPattern().numberOfValidHits()
230 ,transientTracks.at(itk1).hitPattern().numberOfValidPixelHits()
231 ,transientTracks.at(itk1).stateAtBeamLine().transverseImpactParameter().significance()
232 ,transientTracks.at(itk1).normalizedChi2()
233 );
234
235 rootBardak->setEle2(
236 TRootTrack(
237 transientTracks.at(itk2).track().px()
238 ,transientTracks.at(itk2).track().py()
239 ,transientTracks.at(itk2).track().pz()
240 ,transientTracks.at(itk2).track().p()
241 ,transientTracks.at(itk2).track().vx()
242 ,transientTracks.at(itk2).track().vy()
243 ,transientTracks.at(itk2).track().vz()
244 ,0
245 ,transientTracks.at(itk2).track().charge()
246 ,transientTracks.at(itk2).hitPattern().pixelLayersWithMeasurement()
247 ,transientTracks.at(itk2).hitPattern().stripLayersWithMeasurement()
248 ,transientTracks.at(itk2).track().chi2()
249 ,transientTracks.at(itk2).track().d0()
250 ,transientTracks.at(itk2).track().d0Error()
251 ,transientTracks.at(itk2).track().dz()
252 ,transientTracks.at(itk2).track().dzError()
253 )
254 ,transientTracks.at(itk2).hitPattern().numberOfValidHits()
255 ,transientTracks.at(itk2).hitPattern().numberOfValidPixelHits()
256 ,transientTracks.at(itk2).stateAtBeamLine().transverseImpactParameter().significance()
257 ,transientTracks.at(itk2).normalizedChi2()
258 );
259
260 vector<reco::TransientTrack> eeTransientTracks;
261 eeTransientTracks.push_back(transientTracks.at(itk1));
262 eeTransientTracks.push_back(transientTracks.at(itk2));
263 if (itk1<itk2) itk2--;
264 transientTracks.erase (transientTracks.begin()+itk1);
265 transientTracks.erase (transientTracks.begin()+itk2);
266
267 if(verbosity_>5)
268 {
269 for (unsigned int j=0; j<transientTracks.size(); j++)
270 {
271 const reco::Track* tk = & (transientTracks.at(j).track());
272 std::cout << " New Tk[" << j << "] (Et,eta,phi)=(" << tk->pt() << "," << tk->eta() << "," << tk->phi() << ")" << std::endl;
273 }
274 }
275
276 // Call Primary Vertex Reconstruction using all tracks in the event except the two electrons
277 vector<TransientVertex> transientNoEEVertices = theVertexProducer.vertices(transientTracks, vertexBeamSpot);
278
279 reco::VertexCollection noEEVertices;
280 for (vector<TransientVertex>::const_iterator vertex = transientNoEEVertices.begin(); vertex != transientNoEEVertices.end(); vertex++)
281 {
282 reco::Vertex v = *vertex;
283 noEEVertices.push_back(v);
284 }
285
286 if(verbosity_>4) std::cout << std::endl << " Private Primary Vertex reconstruction with all tracks except the two electrons - Number of vertices = " << noEEVertices.size() << std::endl;
287
288 for (unsigned int j=0; j<noEEVertices.size(); j++)
289 {
290 reco::Vertex* vertex = & (noEEVertices.at(j));
291
292 // Put your vertex selection here....
293 if (! vertex->isValid() ) continue;
294 if ( vertex->isFake() ) continue;
295
296 Int_t ntracks = 0;
297 Float_t higherPt = 0.;
298 Float_t scalarSumPt = 0.;
299 Float_t vectorSumPt = 0.;
300 math::XYZVector vectorSum(0.,0.,0.);
301
302 for( std::vector< reco::TrackBaseRef >::const_iterator it = vertex->tracks_begin(); it != vertex->tracks_end(); it++)
303 {
304 scalarSumPt += (**it).pt();
305 vectorSum += (**it).momentum();
306 if( (**it).pt()>higherPt ) higherPt=(**it).pt();
307 ntracks++;
308 }
309 vectorSumPt = sqrt(vectorSum.Perp2());
310
311 // No refitted tracks embeded in reco::Vertex....
312 //cout << "vertex->refittedTracks().size()=" << vertex->refittedTracks().size() << endl;
313
314 TRootVertex localVertex(
315 vertex->x()
316 ,vertex->y()
317 ,vertex->z()
318 ,vertex->xError()
319 ,vertex->yError()
320 ,vertex->zError()
321 );
322
323 localVertex.setAlgoName("NoEE");
324 localVertex.setChi2( vertex->chi2() );
325 localVertex.setNdof( vertex->ndof() );
326 localVertex.setNtracks( ntracks );
327 localVertex.setHigherTrackPt( higherPt );
328 localVertex.setScalarSumPt( scalarSumPt );
329 localVertex.setVectorSumPt( vectorSumPt );
330
331 new( (*rootVertices)[iRootVertex] ) TRootVertex(localVertex);
332 if(verbosity_>2) cout << " ["<< setw(3) << iRootVertex << "] " << localVertex << endl;
333 iRootVertex++;
334 }
335
336
337 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
338 // Private Primary Vertex reconstruction with only the two CTF tracks matched to the selected electrons
339 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
340
341 vector<TransientVertex> transientOnlyEEVertices = theVertexProducer.vertices(eeTransientTracks, vertexBeamSpot);
342 reco::VertexCollection onlyEEVertices;
343
344 for (vector<TransientVertex>::const_iterator vertex = transientOnlyEEVertices.begin(); vertex != transientOnlyEEVertices.end(); vertex++)
345 {
346 reco::Vertex v = *vertex;
347 onlyEEVertices.push_back(v);
348 }
349
350 if(verbosity_>4) std::cout << std::endl << " Private Primary Vertex reconstruction using only the two electron CTF tracks - Number of vertices = " << onlyEEVertices.size() << std::endl;
351
352 for (unsigned int j=0; j<onlyEEVertices.size(); j++)
353 {
354 reco::Vertex* vertex = & (onlyEEVertices.at(j));
355
356 // Put your vertex selection here....
357 if (! vertex->isValid() ) continue;
358 if ( vertex->isFake() ) continue;
359
360 Int_t ntracks = 0;
361 Float_t higherPt = 0.;
362 Float_t scalarSumPt = 0.;
363 Float_t vectorSumPt = 0.;
364 math::XYZVector vectorSum(0.,0.,0.);
365
366 for( std::vector< reco::TrackBaseRef >::const_iterator it = vertex->tracks_begin(); it != vertex->tracks_end(); it++)
367 {
368 scalarSumPt += (**it).pt();
369 vectorSum += (**it).momentum();
370 if( (**it).pt()>higherPt ) higherPt=(**it).pt();
371 ntracks++;
372 }
373 vectorSumPt = sqrt(vectorSum.Perp2());
374
375 // No refitted tracks embeded in reco::Vertex....
376 //cout << "vertex->refittedTracks().size()=" << vertex->refittedTracks().size() << endl;
377
378 TRootVertex localVertex(
379 vertex->x()
380 ,vertex->y()
381 ,vertex->z()
382 ,vertex->xError()
383 ,vertex->yError()
384 ,vertex->zError()
385 );
386
387 localVertex.setAlgoName("OnlyEE");
388 localVertex.setChi2( vertex->chi2() );
389 localVertex.setNdof( vertex->ndof() );
390 localVertex.setNtracks( ntracks );
391 localVertex.setHigherTrackPt( higherPt );
392 localVertex.setScalarSumPt( scalarSumPt );
393 localVertex.setVectorSumPt( vectorSumPt );
394
395 new( (*rootVertices)[iRootVertex] ) TRootVertex(localVertex);
396 if(verbosity_>2) cout << " ["<< setw(3) << iRootVertex << "] " << localVertex << endl;
397 iRootVertex++;
398 }
399
400 return true;
401 }