ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Morgan/src/ZeeVertexAnalyzer.cc
Revision: 1.2
Committed: Tue Oct 13 14:07:25 2009 UTC (15 years, 6 months ago) by lethuill
Content type: text/plain
Branch: MAIN
Changes since 1.1: +77 -9 lines
Log Message:
Add Primary Vertex reconstruction using two electron tracks only

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     bool ZeeVertexAnalyzer::getVertices(const edm::Event& iEvent, const edm::EventSetup& iSetup, TClonesArray* rootVertices)
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     if(verbosity_>4) std::cout << " CTF Track 1 - (Et,eta,phi)=(" << tk1->pt() << "," << tk1->eta() << "," << tk1->phi() << ") - Fraction of common hits with GSF =" << (electrons->at(iElectron1)).shFracInnerHits() << std::endl;
84     if(verbosity_>4) std::cout << " CTF Track 2 - (Et,eta,phi)=(" << tk2->pt() << "," << tk2->eta() << "," << tk2->phi() << ") - Fraction of common hits with GSF =" << (electrons->at(iElectron2)).shFracInnerHits() << std::endl;
85    
86    
87     ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
88     // Private Primary Vertex reconstruction with all CTF tracks
89     ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
90    
91     // Get BeamSpot
92     reco::BeamSpot vertexBeamSpot;
93     edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
94     iEvent.getByLabel(beamSpotProducer_, recoBeamSpotHandle);
95     if (recoBeamSpotHandle.isValid()) vertexBeamSpot = *recoBeamSpotHandle;
96    
97     // Get Track Collection
98     edm::Handle<reco::TrackCollection> recoTracks;
99     iEvent.getByLabel(trackProducer_, recoTracks);
100    
101     // Get Transient Track Collection
102     edm::ESHandle<TransientTrackBuilder> trackBuilder;
103     iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",trackBuilder);
104     vector<reco::TransientTrack> transientTracks = (*trackBuilder).build(recoTracks, vertexBeamSpot);
105    
106     // Call Primary Vertex Reconstruction using all tracks in the event
107     PrimaryVertexProducerAlgorithm theVertexProducer(config_);
108     vector<TransientVertex> transientZeeVertices = theVertexProducer.vertices(transientTracks, vertexBeamSpot);
109    
110     reco::VertexCollection zeeVertices;
111     for (vector<TransientVertex>::const_iterator vertex = transientZeeVertices.begin(); vertex != transientZeeVertices.end(); vertex++)
112     {
113     reco::Vertex v = *vertex;
114     zeeVertices.push_back(v);
115     }
116    
117     if(verbosity_>4) std::cout << std::endl << " Private Primary Vertex reconstruction with all tracks - Number of vertices = " << zeeVertices.size() << std::endl;
118 lethuill 1.2
119 lethuill 1.1 for (unsigned int j=0; j<zeeVertices.size(); j++)
120     {
121     reco::Vertex* vertex = & (zeeVertices.at(j));
122    
123     // Put your vertex selection here....
124     if (! vertex->isValid() ) continue;
125     if ( vertex->isFake() ) continue;
126    
127     Int_t ntracks = 0;
128     Float_t higherPt = 0.;
129     Float_t scalarSumPt = 0.;
130     Float_t vectorSumPt = 0.;
131     math::XYZVector vectorSum(0.,0.,0.);
132    
133     for( std::vector< reco::TrackBaseRef >::const_iterator it = vertex->tracks_begin(); it != vertex->tracks_end(); it++)
134     {
135     scalarSumPt += (**it).pt();
136     vectorSum += (**it).momentum();
137     if( (**it).pt()>higherPt ) higherPt=(**it).pt();
138     ntracks++;
139     }
140     vectorSumPt = sqrt(vectorSum.Perp2());
141    
142     // No refitted tracks embeded in reco::Vertex....
143     //cout << "vertex->refittedTracks().size()=" << vertex->refittedTracks().size() << endl;
144    
145     TRootVertex localVertex(
146     vertex->x()
147     ,vertex->y()
148     ,vertex->z()
149     ,vertex->xError()
150     ,vertex->yError()
151     ,vertex->zError()
152     );
153    
154     localVertex.setAlgoName("AllTracks");
155     localVertex.setChi2( vertex->chi2() );
156     localVertex.setNdof( vertex->ndof() );
157     localVertex.setNtracks( ntracks );
158     localVertex.setHigherTrackPt( higherPt );
159     localVertex.setScalarSumPt( scalarSumPt );
160     localVertex.setVectorSumPt( vectorSumPt );
161    
162     new( (*rootVertices)[iRootVertex] ) TRootVertex(localVertex);
163     if(verbosity_>2) cout << " ["<< setw(3) << iRootVertex << "] " << localVertex << endl;
164     iRootVertex++;
165     }
166    
167    
168     ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
169     // Private Primary Vertex reconstruction with all CTF tracks except the 2 two tracks matching selected electrons
170     ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
171    
172 lethuill 1.2 // Remove the two CTF tracks matching the two electron GSF tracks from the transientTracks collection
173     // Keep the two CTF electron tracks in a separate vector (eeTransientTracks)
174 lethuill 1.1 unsigned int itk1 = 999999;
175     unsigned int itk2 = 999999;
176 lethuill 1.2
177 lethuill 1.1 for (unsigned int j=0; j<transientTracks.size(); j++)
178     {
179     const reco::Track* tk = & (transientTracks.at(j).track());
180     if(verbosity_>5) std::cout << " Tk[" << j << "] (Et,eta,phi)=(" << tk->pt() << "," << tk->eta() << "," << tk->phi() << ")" << std::endl;
181     if ( tk->pt()==tk1->pt() && tk->eta()==tk1->eta() && tk->phi()==tk1->phi() )
182     {
183     itk1 = j;
184     if(verbosity_>5) std::cout << " ===> this is tk1" << std::endl;
185     }
186     if ( tk->pt()==tk2->pt() && tk->eta()==tk2->eta() && tk->phi()==tk2->phi() )
187     {
188     itk2 = j;
189     if(verbosity_>5) std::cout << " ===> this is tk2" << std::endl;
190     }
191     }
192     if (itk1==999999 || itk2==999999)
193     {
194     std::cout << " ***** NO CTF TRACKS MATCHING THE ELECTRON *****" << std::endl;
195     return false;
196     }
197    
198 lethuill 1.2 vector<reco::TransientTrack> eeTransientTracks;
199     eeTransientTracks.push_back(transientTracks.at(itk1));
200     eeTransientTracks.push_back(transientTracks.at(itk2));
201 lethuill 1.1 if (itk1<itk2) itk2--;
202     transientTracks.erase (transientTracks.begin()+itk1);
203     transientTracks.erase (transientTracks.begin()+itk2);
204    
205     if(verbosity_>5)
206     {
207     for (unsigned int j=0; j<transientTracks.size(); j++)
208     {
209     const reco::Track* tk = & (transientTracks.at(j).track());
210     std::cout << " New Tk[" << j << "] (Et,eta,phi)=(" << tk->pt() << "," << tk->eta() << "," << tk->phi() << ")" << std::endl;
211     }
212     }
213    
214     // Call Primary Vertex Reconstruction using all tracks in the event except the two electrons
215     vector<TransientVertex> transientNoEEVertices = theVertexProducer.vertices(transientTracks, vertexBeamSpot);
216    
217     reco::VertexCollection noEEVertices;
218     for (vector<TransientVertex>::const_iterator vertex = transientNoEEVertices.begin(); vertex != transientNoEEVertices.end(); vertex++)
219     {
220     reco::Vertex v = *vertex;
221     noEEVertices.push_back(v);
222     }
223    
224     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;
225    
226     for (unsigned int j=0; j<noEEVertices.size(); j++)
227     {
228     reco::Vertex* vertex = & (noEEVertices.at(j));
229    
230     // Put your vertex selection here....
231     if (! vertex->isValid() ) continue;
232     if ( vertex->isFake() ) continue;
233 lethuill 1.2
234 lethuill 1.1 Int_t ntracks = 0;
235     Float_t higherPt = 0.;
236     Float_t scalarSumPt = 0.;
237     Float_t vectorSumPt = 0.;
238     math::XYZVector vectorSum(0.,0.,0.);
239 lethuill 1.2
240 lethuill 1.1 for( std::vector< reco::TrackBaseRef >::const_iterator it = vertex->tracks_begin(); it != vertex->tracks_end(); it++)
241     {
242     scalarSumPt += (**it).pt();
243     vectorSum += (**it).momentum();
244     if( (**it).pt()>higherPt ) higherPt=(**it).pt();
245     ntracks++;
246     }
247     vectorSumPt = sqrt(vectorSum.Perp2());
248 lethuill 1.2
249 lethuill 1.1 // No refitted tracks embeded in reco::Vertex....
250     //cout << "vertex->refittedTracks().size()=" << vertex->refittedTracks().size() << endl;
251 lethuill 1.2
252 lethuill 1.1 TRootVertex localVertex(
253     vertex->x()
254     ,vertex->y()
255     ,vertex->z()
256     ,vertex->xError()
257     ,vertex->yError()
258     ,vertex->zError()
259     );
260 lethuill 1.2
261 lethuill 1.1 localVertex.setAlgoName("NoEE");
262     localVertex.setChi2( vertex->chi2() );
263     localVertex.setNdof( vertex->ndof() );
264     localVertex.setNtracks( ntracks );
265     localVertex.setHigherTrackPt( higherPt );
266     localVertex.setScalarSumPt( scalarSumPt );
267     localVertex.setVectorSumPt( vectorSumPt );
268 lethuill 1.2
269     new( (*rootVertices)[iRootVertex] ) TRootVertex(localVertex);
270     if(verbosity_>2) cout << " ["<< setw(3) << iRootVertex << "] " << localVertex << endl;
271     iRootVertex++;
272     }
273 lethuill 1.1
274 lethuill 1.2
275     ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
276     // Private Primary Vertex reconstruction with only the two CTF tracks matched to the selected electrons
277     ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
278    
279     vector<TransientVertex> transientOnlyEEVertices = theVertexProducer.vertices(eeTransientTracks, vertexBeamSpot);
280     reco::VertexCollection onlyEEVertices;
281    
282     for (vector<TransientVertex>::const_iterator vertex = transientOnlyEEVertices.begin(); vertex != transientOnlyEEVertices.end(); vertex++)
283     {
284     reco::Vertex v = *vertex;
285     onlyEEVertices.push_back(v);
286     }
287    
288     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;
289    
290     for (unsigned int j=0; j<onlyEEVertices.size(); j++)
291     {
292     reco::Vertex* vertex = & (onlyEEVertices.at(j));
293    
294     // Put your vertex selection here....
295     if (! vertex->isValid() ) continue;
296     if ( vertex->isFake() ) continue;
297    
298     Int_t ntracks = 0;
299     Float_t higherPt = 0.;
300     Float_t scalarSumPt = 0.;
301     Float_t vectorSumPt = 0.;
302     math::XYZVector vectorSum(0.,0.,0.);
303    
304     for( std::vector< reco::TrackBaseRef >::const_iterator it = vertex->tracks_begin(); it != vertex->tracks_end(); it++)
305     {
306     scalarSumPt += (**it).pt();
307     vectorSum += (**it).momentum();
308     if( (**it).pt()>higherPt ) higherPt=(**it).pt();
309     ntracks++;
310     }
311     vectorSumPt = sqrt(vectorSum.Perp2());
312    
313     // No refitted tracks embeded in reco::Vertex....
314     //cout << "vertex->refittedTracks().size()=" << vertex->refittedTracks().size() << endl;
315    
316     TRootVertex localVertex(
317     vertex->x()
318     ,vertex->y()
319     ,vertex->z()
320     ,vertex->xError()
321     ,vertex->yError()
322     ,vertex->zError()
323     );
324    
325     localVertex.setAlgoName("OnlyEE");
326     localVertex.setChi2( vertex->chi2() );
327     localVertex.setNdof( vertex->ndof() );
328     localVertex.setNtracks( ntracks );
329     localVertex.setHigherTrackPt( higherPt );
330     localVertex.setScalarSumPt( scalarSumPt );
331     localVertex.setVectorSumPt( vectorSumPt );
332    
333 lethuill 1.1 new( (*rootVertices)[iRootVertex] ) TRootVertex(localVertex);
334     if(verbosity_>2) cout << " ["<< setw(3) << iRootVertex << "] " << localVertex << endl;
335     iRootVertex++;
336     }
337    
338 lethuill 1.2 return true;
339 lethuill 1.1 }