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
Log Message:
Stock infos on two electron tracks in TRootBardak

File Contents

# User Rev Content
1 lethuill 1.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 lethuill 1.3 bool ZeeVertexAnalyzer::getVertices(const edm::Event& iEvent, const edm::EventSetup& iSetup, TClonesArray* rootVertices, TRootBardak* rootBardak)
23 lethuill 1.1 {
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 lethuill 1.3 //if( p4Z.mass()<60.0 || p4Z.mass()>120.0) return false;
78 lethuill 1.1
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 lethuill 1.2
117 lethuill 1.1 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 lethuill 1.2 // 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 lethuill 1.1 unsigned int itk1 = 999999;
173     unsigned int itk2 = 999999;
174 lethuill 1.2
175 lethuill 1.1 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 lethuill 1.3
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 lethuill 1.1
260 lethuill 1.2 vector<reco::TransientTrack> eeTransientTracks;
261     eeTransientTracks.push_back(transientTracks.at(itk1));
262     eeTransientTracks.push_back(transientTracks.at(itk2));
263 lethuill 1.1 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 lethuill 1.2
296 lethuill 1.1 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 lethuill 1.2
302 lethuill 1.1 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 lethuill 1.2
311 lethuill 1.1 // No refitted tracks embeded in reco::Vertex....
312     //cout << "vertex->refittedTracks().size()=" << vertex->refittedTracks().size() << endl;
313 lethuill 1.2
314 lethuill 1.1 TRootVertex localVertex(
315     vertex->x()
316     ,vertex->y()
317     ,vertex->z()
318     ,vertex->xError()
319     ,vertex->yError()
320     ,vertex->zError()
321     );
322 lethuill 1.2
323 lethuill 1.1 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 lethuill 1.2
331     new( (*rootVertices)[iRootVertex] ) TRootVertex(localVertex);
332     if(verbosity_>2) cout << " ["<< setw(3) << iRootVertex << "] " << localVertex << endl;
333     iRootVertex++;
334     }
335 lethuill 1.1
336 lethuill 1.2
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 lethuill 1.1 new( (*rootVertices)[iRootVertex] ) TRootVertex(localVertex);
396     if(verbosity_>2) cout << " ["<< setw(3) << iRootVertex << "] " << localVertex << endl;
397     iRootVertex++;
398     }
399    
400 lethuill 1.2 return true;
401 lethuill 1.1 }