Morgan/src/VertexAnalyzer.cc

#include "../interface/VertexAnalyzer.h"

using namespace std;
using namespace reco;
using namespace edm;

VertexAnalyzer::VertexAnalyzer(const edm::ParameterSet& producersNames):verbosity_(0)
{
        primaryVertexProducer_ = producersNames.getParameter<edm::InputTag>("primaryVertexProducer");
        beamSpotProducer_ = producersNames.getParameter<edm::InputTag>("beamSpotProducer");
        allowMissingCollection_ = producersNames.getUntrackedParameter<bool>("allowMissingCollection", false);
}


VertexAnalyzer::VertexAnalyzer(const edm::ParameterSet& producersNames, int verbosity):verbosity_(verbosity)
{
        primaryVertexProducer_ = producersNames.getParameter<edm::InputTag>("primaryVertexProducer");
        beamSpotProducer_ = producersNames.getParameter<edm::InputTag>("beamSpotProducer");
        allowMissingCollection_ = producersNames.getUntrackedParameter<bool>("allowMissingCollection", false);
}


VertexAnalyzer::~VertexAnalyzer()
{
}


bool VertexAnalyzer::getBeamSpot(const edm::Event& iEvent, TRootBeamSpot* rootBeamSpot)
{
        try
        {
                edm::Handle<reco::BeamSpot> recoBeamSpot;
                iEvent.getByLabel(beamSpotProducer_, recoBeamSpot);
                
                if ( recoBeamSpot.isValid() )
                {
                        if(verbosity_>1) std::cout  << "   Beam Spot - Label: " << beamSpotProducer_.label() << "   Instance: " << beamSpotProducer_.instance() << std::endl;
                        rootBeamSpot->fill(
                        recoBeamSpot->x0()
                        ,recoBeamSpot->y0()
                        ,recoBeamSpot->z0()
                        ,recoBeamSpot->x0Error()
                        ,recoBeamSpot->y0Error()
                        ,recoBeamSpot->z0Error()
                        ,recoBeamSpot->sigmaZ()
                        ,recoBeamSpot->sigmaZ0Error()
                        ,recoBeamSpot->BeamWidthX()
                        ,recoBeamSpot->BeamWidthXError()
                        ,recoBeamSpot->BeamWidthY()
                        ,recoBeamSpot->BeamWidthYError()
                        );
                        if(verbosity_>2) cout << "   "<< *rootBeamSpot << endl << endl;
                }
                else
                {
                        if(verbosity_>1) cout << "  ##### ERROR IN  VertexAnalyzer::getBeamSpot => No beam spot available from EventSetup #####"<<endl;
                        return false;
                }
        }
        catch (cms::Exception& exception)
        {
                if ( !allowMissingCollection_ )
                {
                        cout << "  ##### ERROR IN  VertexAnalyzer::getBeamSpot => No beam spot available from EventSetup #####"<<endl;
                        throw exception;
                }
                if(verbosity_>1) cout <<  "   ===> No beam spot available, skip beam spot info" << endl;
                return false;
        }
        
        return true;
}


bool VertexAnalyzer::getVertices(const edm::Event& iEvent, TClonesArray* rootVertices)
{
        
        edm::Handle< reco::VertexCollection > recoVertices;
        try
        {
                iEvent.getByLabel(primaryVertexProducer_, recoVertices);
                int nVertices = recoVertices->size();
                if(verbosity_>1) std::cout  << "   Number of primary vertices = " << nVertices << "   Label: " << primaryVertexProducer_.label() << "   Instance: " << primaryVertexProducer_.instance() << std::endl;
        }
        catch (cms::Exception& exception)
        {
                if ( !allowMissingCollection_ )
                {
                        cout << "  ##### ERROR IN  VertexAnalyzer::getVertices => Vertex collection is missing #####"<<endl;
                        throw exception;
                }
                if(verbosity_>1) cout <<  "   ===> No primaryVertex collection, skip vertex info" << endl;
                return false;
        }
        
        int iRootVertex = 0;
        for (unsigned int j=0; j<recoVertices->size(); j++)
        {
                const reco::Vertex* vertex = & ((*recoVertices)[j]);
                
                // Put your vertex selection here....
                if (! vertex->isValid() ) continue;
                if ( vertex->isFake() ) continue;
                
                Int_t ntracks = 0;
                Float_t higherPt = 0.;
                Float_t scalarSumPt = 0.;
                Float_t vectorSumPt = 0.;
                math::XYZVector vectorSum(0.,0.,0.);
                
                for( std::vector< reco::TrackBaseRef >::const_iterator it = vertex->tracks_begin(); it != vertex->tracks_end(); it++)
                {
                        scalarSumPt += (**it).pt();
                        vectorSum += (**it).momentum();
                        if( (**it).pt()>higherPt ) higherPt=(**it).pt();
                        ntracks++;
                }
                vectorSumPt = sqrt(vectorSum.Perp2());
                
                // No refitted tracks embeded in reco::Vertex....
                //cout  << "vertex->refittedTracks().size()=" << vertex->refittedTracks().size() << endl;
                
                TRootVertex localVertex(
                vertex->x()
                ,vertex->y()
                ,vertex->z()
                ,vertex->xError()
                ,vertex->yError()
                ,vertex->zError()
                );
                
                localVertex.setAlgoName("RECO");
                localVertex.setChi2( vertex->chi2() );
                localVertex.setNdof( vertex->ndof() );
                localVertex.setNtracks( ntracks );
                localVertex.setHigherTrackPt( higherPt );
                localVertex.setScalarSumPt( scalarSumPt );
                localVertex.setVectorSumPt( vectorSumPt );
                
                new( (*rootVertices)[iRootVertex] ) TRootVertex(localVertex);
                if(verbosity_>2) cout << "   ["<< setw(3) << iRootVertex << "] " << localVertex << endl;
                iRootVertex++;
        }
        
        return true;
}


void VertexAnalyzer::selectPrimary(TRootEvent* rootEvent, TClonesArray* rootVertices)
{
        
        int ibestVertex = -1;
        float highestScalarSumPt = -1.;
        TRootVertex* localVertex = 0;
        
        for (int ivtx=0; ivtx<rootVertices->GetEntriesFast(); ivtx++)
        {
                localVertex = (TRootVertex*) rootVertices->At(ivtx);
                // Put you primary vertex selection here...
                if ( localVertex->scalarSumPt()>highestScalarSumPt )
                {
                        highestScalarSumPt = localVertex->scalarSumPt();
                        ibestVertex = ivtx;
                }
        }
        
        if ( ibestVertex >= 0)
        {
                localVertex = (TRootVertex*) rootVertices->At(ibestVertex);
                rootEvent->setPrimaryVertexIndex( ibestVertex );
                rootEvent->setPrimaryVertex( localVertex );
        }
        
        if(verbosity_>1) std::cout  << "   Selected primary vertex: TRootVertex[" << ibestVertex << "]" << std::endl;
}
Revision:	1.9
Committed:	Mon Oct 12 08:32:36 2009 UTC (15 years, 6 months ago) by lethuill
Content type:	text/plain
Branch:	MAIN
Changes since 1.8:	+141 -136 lines
Log Message:	Add private primary vertex reconstruction for Zee events
#	User	Rev	Content
1	lethuill	1.2	#include "../interface/VertexAnalyzer.h"
2	lethuill	1.1
3			using namespace std;
4			using namespace reco;
5			using namespace edm;
6
7			VertexAnalyzer::VertexAnalyzer(const edm::ParameterSet& producersNames):verbosity_(0)
8			{
9	lethuill	1.9	primaryVertexProducer_ = producersNames.getParameter<edm::InputTag>("primaryVertexProducer");
10			beamSpotProducer_ = producersNames.getParameter<edm::InputTag>("beamSpotProducer");
11			allowMissingCollection_ = producersNames.getUntrackedParameter<bool>("allowMissingCollection", false);
12	lethuill	1.1	}
13
14	lethuill	1.4
15	lethuill	1.1	VertexAnalyzer::VertexAnalyzer(const edm::ParameterSet& producersNames, int verbosity):verbosity_(verbosity)
16			{
17	lethuill	1.9	primaryVertexProducer_ = producersNames.getParameter<edm::InputTag>("primaryVertexProducer");
18			beamSpotProducer_ = producersNames.getParameter<edm::InputTag>("beamSpotProducer");
19			allowMissingCollection_ = producersNames.getUntrackedParameter<bool>("allowMissingCollection", false);
20	lethuill	1.1	}
21
22	lethuill	1.4
23	lethuill	1.1	VertexAnalyzer::~VertexAnalyzer()
24			{
25			}
26
27	lethuill	1.4
28	lethuill	1.7	bool VertexAnalyzer::getBeamSpot(const edm::Event& iEvent, TRootBeamSpot* rootBeamSpot)
29	lethuill	1.6	{
30	lethuill	1.9	try
31			{
32			edm::Handle<reco::BeamSpot> recoBeamSpot;
33			iEvent.getByLabel(beamSpotProducer_, recoBeamSpot);
34
35			if ( recoBeamSpot.isValid() )
36			{
37			if(verbosity_>1) std::cout << " Beam Spot - Label: " << beamSpotProducer_.label() << " Instance: " << beamSpotProducer_.instance() << std::endl;
38			rootBeamSpot->fill(
39			recoBeamSpot->x0()
40			,recoBeamSpot->y0()
41			,recoBeamSpot->z0()
42			,recoBeamSpot->x0Error()
43			,recoBeamSpot->y0Error()
44			,recoBeamSpot->z0Error()
45			,recoBeamSpot->sigmaZ()
46			,recoBeamSpot->sigmaZ0Error()
47			,recoBeamSpot->BeamWidthX()
48			,recoBeamSpot->BeamWidthXError()
49			,recoBeamSpot->BeamWidthY()
50			,recoBeamSpot->BeamWidthYError()
51			);
52			if(verbosity_>2) cout << " "<< *rootBeamSpot << endl << endl;
53			}
54			else
55			{
56			if(verbosity_>1) cout << " ##### ERROR IN VertexAnalyzer::getBeamSpot => No beam spot available from EventSetup #####"<<endl;
57			return false;
58			}
59			}
60			catch (cms::Exception& exception)
61			{
62			if ( !allowMissingCollection_ )
63			{
64			cout << " ##### ERROR IN VertexAnalyzer::getBeamSpot => No beam spot available from EventSetup #####"<<endl;
65			throw exception;
66			}
67			if(verbosity_>1) cout << " ===> No beam spot available, skip beam spot info" << endl;
68			return false;
69			}
70
71			return true;
72	lethuill	1.6	}
73
74
75	lethuill	1.7	bool VertexAnalyzer::getVertices(const edm::Event& iEvent, TClonesArray* rootVertices)
76	lethuill	1.1	{
77	lethuill	1.9
78			edm::Handle< reco::VertexCollection > recoVertices;
79			try
80			{
81			iEvent.getByLabel(primaryVertexProducer_, recoVertices);
82			int nVertices = recoVertices->size();
83			if(verbosity_>1) std::cout << " Number of primary vertices = " << nVertices << " Label: " << primaryVertexProducer_.label() << " Instance: " << primaryVertexProducer_.instance() << std::endl;
84			}
85			catch (cms::Exception& exception)
86			{
87			if ( !allowMissingCollection_ )
88			{
89			cout << " ##### ERROR IN VertexAnalyzer::getVertices => Vertex collection is missing #####"<<endl;
90			throw exception;
91			}
92			if(verbosity_>1) cout << " ===> No primaryVertex collection, skip vertex info" << endl;
93			return false;
94			}
95
96			int iRootVertex = 0;
97			for (unsigned int j=0; j<recoVertices->size(); j++)
98			{
99			const reco::Vertex* vertex = & ((*recoVertices)[j]);
100
101			// Put your vertex selection here....
102			if (! vertex->isValid() ) continue;
103			if ( vertex->isFake() ) continue;
104
105			Int_t ntracks = 0;
106			Float_t higherPt = 0.;
107			Float_t scalarSumPt = 0.;
108			Float_t vectorSumPt = 0.;
109			math::XYZVector vectorSum(0.,0.,0.);
110
111			for( std::vector< reco::TrackBaseRef >::const_iterator it = vertex->tracks_begin(); it != vertex->tracks_end(); it++)
112			{
113			scalarSumPt += (**it).pt();
114			vectorSum += (**it).momentum();
115			if( (it).pt()>higherPt ) higherPt=(it).pt();
116			ntracks++;
117			}
118			vectorSumPt = sqrt(vectorSum.Perp2());
119
120			// No refitted tracks embeded in reco::Vertex....
121			//cout << "vertex->refittedTracks().size()=" << vertex->refittedTracks().size() << endl;
122
123			TRootVertex localVertex(
124			vertex->x()
125			,vertex->y()
126			,vertex->z()
127			,vertex->xError()
128			,vertex->yError()
129			,vertex->zError()
130			);
131
132			localVertex.setAlgoName("RECO");
133			localVertex.setChi2( vertex->chi2() );
134			localVertex.setNdof( vertex->ndof() );
135			localVertex.setNtracks( ntracks );
136			localVertex.setHigherTrackPt( higherPt );
137			localVertex.setScalarSumPt( scalarSumPt );
138			localVertex.setVectorSumPt( vectorSumPt );
139
140			new( (*rootVertices)[iRootVertex] ) TRootVertex(localVertex);
141			if(verbosity_>2) cout << " ["<< setw(3) << iRootVertex << "] " << localVertex << endl;
142			iRootVertex++;
143			}
144
145			return true;
146	lethuill	1.4	}
147
148
149	lethuill	1.6	void VertexAnalyzer::selectPrimary(TRootEvent* rootEvent, TClonesArray* rootVertices)
150	lethuill	1.4	{
151	lethuill	1.9
152			int ibestVertex = -1;
153			float highestScalarSumPt = -1.;
154			TRootVertex* localVertex = 0;
155
156			for (int ivtx=0; ivtx<rootVertices->GetEntriesFast(); ivtx++)
157			{
158			localVertex = (TRootVertex*) rootVertices->At(ivtx);
159			// Put you primary vertex selection here...
160			if ( localVertex->scalarSumPt()>highestScalarSumPt )
161			{
162			highestScalarSumPt = localVertex->scalarSumPt();
163			ibestVertex = ivtx;
164			}
165			}
166
167			if ( ibestVertex >= 0)
168			{
169			localVertex = (TRootVertex*) rootVertices->At(ibestVertex);
170			rootEvent->setPrimaryVertexIndex( ibestVertex );
171			rootEvent->setPrimaryVertex( localVertex );
172			}
173
174			if(verbosity_>1) std::cout << " Selected primary vertex: TRootVertex[" << ibestVertex << "]" << std::endl;
175	lethuill	1.1	}