Morgan/src/ElectronAnalyzer.cc

#include "../interface/ElectronAnalyzer.h"

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

ElectronAnalyzer::ElectronAnalyzer(const edm::ParameterSet& producersNames):verbosity_(0),useMC_(false)
{
        dataType_ = producersNames.getUntrackedParameter<string>("dataType","unknown");
        electronProducer_ = producersNames.getParameter<edm::InputTag>("electronProducer");
}

ElectronAnalyzer::ElectronAnalyzer(const edm::ParameterSet& producersNames, const edm::ParameterSet& myConfig, int verbosity):verbosity_(verbosity)
{
        dataType_ = producersNames.getUntrackedParameter<string>("dataType","unknown");
        electronProducer_ = producersNames.getParameter<edm::InputTag>("electronProducer");
        useMC_ = myConfig.getUntrackedParameter<bool>("doElectronMC");
}

ElectronAnalyzer::~ElectronAnalyzer()
{
}

void ElectronAnalyzer::Process(const edm::Event& iEvent, TClonesArray* rootElectrons, EcalClusterLazyTools& lazyTools )
{

        unsigned int nElectrons=0;

        edm::Handle < std::vector <reco::GsfElectron> > recoElectrons;
        if( dataType_=="RECO" || dataType_=="AOD" )
        {
                iEvent.getByLabel(electronProducer_, recoElectrons);
                nElectrons = recoElectrons->size();
        }

        edm::Handle < std::vector <pat::Electron> > patElectrons;
        if( dataType_=="PAT" || dataType_=="PATAOD" )
        {
                iEvent.getByLabel(electronProducer_, patElectrons);
                nElectrons = patElectrons->size();
        }

        if(verbosity_>1) std::cout << "   Number of electrons = " << nElectrons << "   Label: " << electronProducer_.label() << "   Instance: " << electronProducer_.instance() << std::endl;

        for (unsigned int j=0; j<nElectrons; j++)
        {
                const reco::GsfElectron* electron = 0;  
                if( dataType_=="RECO" || dataType_=="AOD" ) electron =  &((*recoElectrons)[j]);
                if( dataType_=="PAT" || dataType_=="PATAOD" ) electron = (const reco::GsfElectron*) ( & ((*patElectrons)[j]) );

                TRootElectron localElectron(
                        electron->px()
                        ,electron->py()
                        ,electron->pz()
                        ,electron->energy()
                        ,electron->vx()
                        ,electron->vy()
                        ,electron->vz()
                        ,abs(electron->pdgId())
                        ,electron->charge()
                ); 

                // Variables from reco::GsfElectron
                localElectron.setClassification(electron->classification());
                localElectron.setCaloEnergy(electron->caloEnergy());
                localElectron.setCaloEnergyError(electron->caloEnergyError());
                // FIXME - Protect against sqrt(0)
                localElectron.setTrackMomentum(sqrt(electron->trackMomentumAtVtx().Mag2()));
                localElectron.setTrackMomentumError(electron->trackMomentumError());
                localElectron.setHadOverEm(electron->hadronicOverEm());
                localElectron.setDeltaEtaIn(electron->deltaEtaSuperClusterTrackAtVtx());
                localElectron.setDeltaPhiIn(electron->deltaPhiSuperClusterTrackAtVtx());
                localElectron.setEnergySuperClusterOverPin(electron->eSuperClusterOverP());
                localElectron.setDeltaEtaOut(electron->deltaEtaSeedClusterTrackAtCalo());
                localElectron.setDeltaPhiOut(electron->deltaPhiSeedClusterTrackAtCalo());
                localElectron.setEnergySeedClusterOverPout(electron->eSeedClusterOverPout());
                localElectron.setEnergyScaleCorrected(electron->isEnergyScaleCorrected());
                localElectron.setMomentumCorrected(electron->isMomentumCorrected());

                // Variables from reco::GsfTrack
                reco::GsfTrackRef gsfTrack = electron->gsfTrack();
                if ( gsfTrack.isNonnull() )
                {
                        localElectron.setD0(gsfTrack->d0());
                        localElectron.setD0Error(gsfTrack->d0Error());
                        localElectron.setDsz(gsfTrack->dsz());
                        localElectron.setDszError(gsfTrack->dszError());
                        localElectron.setNormalizedChi2(gsfTrack->normalizedChi2());
                        localElectron.setPtError(gsfTrack->ptError());
                        localElectron.setEtaError(gsfTrack->etaError());
                        localElectron.setPhiError(gsfTrack->phiError());
                }

                // Variables from reco::SuperCluster
                reco::SuperClusterRef superCluster = electron->superCluster();
                //const reco::SuperCluster* sc = superCluster.get();
                if ( superCluster.isNonnull() )
                {
                        localElectron.setNbClusters(electron->numberOfClusters());
                        localElectron.setSuperClusterRawEnergy(superCluster->rawEnergy());
                        localElectron.setPreshowerEnergy(superCluster->preshowerEnergy());
                        localElectron.setCaloPosition(
                                superCluster->position().X()
                                ,superCluster->position().Y()
                                ,superCluster->position().Z()
                        );

                        // FIXME - associator supercluster <-> electron
                        //localElectron.setSCRef(superCluster->toto());
                }               


                // Cluster Shape variables
                // need reco::SuperCluster and reco::BasicCluster
                if ( superCluster.isNonnull() )
                {
                        reco::BasicClusterRef seedBasicCluster = superCluster->seed();
                        if ( seedBasicCluster.isNonnull() ) localElectron.setClusterAlgo(seedBasicCluster->algo());

                        // dR of the cone centered on the reco::GsfElectron and containing all its basic clusters constituents
                        Float_t caloConeSize = 0;
                        for (reco::basicCluster_iterator basicCluster = (*superCluster).clustersBegin(); basicCluster != (*superCluster).clustersEnd(); ++basicCluster )
                        {
                                Float_t dR = localElectron.DeltaR(TLorentzVector( (*basicCluster)->position().x(), (*basicCluster)->position().y(), (*basicCluster)->position().z(), 0. ) );
                                if (dR > caloConeSize) caloConeSize = dR;
                        }
                        // FIXME - if no BasicCluster collection, init to -999.
                        localElectron.setCaloConeSize(caloConeSize);

                        // need reduced Ecal RecHits Collections for EcalClusterLazyTools
                        // FIXME - check lazyTools is initialized ==> lazyTools has to be a pointer
                        if ( seedBasicCluster.isNonnull() ) 
                        {
                                localElectron.setE2x2(lazyTools.e2x2(*seedBasicCluster));
                                localElectron.setE3x3(lazyTools.e3x3(*seedBasicCluster));
                                localElectron.setE5x5(lazyTools.e5x5(*seedBasicCluster));
                                localElectron.setEMax(lazyTools.eMax(*seedBasicCluster));
                        }
                }


                if( dataType_=="RECO" || dataType_=="AOD" )
                {
                        // Some specific methods requiring  RECO / AOD format
                        // Do association to genParticle ?

                        // Isolation ?
                        // Electron ID ?  RecoEgamma/ElectronIdentification
                }


                if( dataType_=="PATAOD" || dataType_=="PAT" )
                {
                        // Some specific methods to pat::Electron
                        const pat::Electron *patElectron = dynamic_cast<const pat::Electron*>(&*electron);

                        // Isolation
                        // FIXME - how isoDeposit are calculated ? need to substract supercluster energy ?
                        pair < Float_t, Int_t > trackerIso;

                        trackerIso = patElectron->trackerIsoDeposit()->depositAndCountWithin(0.1);
                        localElectron.setIsoR01_sumPt(trackerIso.first - gsfTrack->p() );
                        localElectron.setIsoR01_nTracks(trackerIso.second);

                        trackerIso = patElectron->trackerIsoDeposit()->depositAndCountWithin(0.2);
                        localElectron.setIsoR02_sumPt(trackerIso.first - gsfTrack->p() );
                        localElectron.setIsoR02_nTracks(trackerIso.second);

                        trackerIso = patElectron->trackerIsoDeposit()->depositAndCountWithin(0.3);
                        localElectron.setIsoR03_sumPt(trackerIso.first - gsfTrack->p() );
                        localElectron.setIsoR03_nTracks(trackerIso.second);
                        localElectron.setIsoR03_emEt( patElectron->ecalIsoDeposit()->depositAndCountWithin(0.3).first - superCluster->rawEnergy() );
                        localElectron.setIsoR03_hadEt( patElectron->hcalIsoDeposit()->depositAndCountWithin(0.3).first );

                        trackerIso = patElectron->trackerIsoDeposit()->depositAndCountWithin(0.5);
                        localElectron.setIsoR05_sumPt(trackerIso.first - gsfTrack->p() );
                        localElectron.setIsoR05_nTracks(trackerIso.second);
                        localElectron.setIsoR05_emEt( patElectron->ecalIsoDeposit()->depositAndCountWithin(0.5).first - superCluster->rawEnergy() );
                        localElectron.setIsoR05_hadEt( patElectron->hcalIsoDeposit()->depositAndCountWithin(0.5).first );


                        // Electron ID
                        localElectron.setIDPTDRLoose(patElectron->leptonID("ptdrLoose"));
                        localElectron.setIDPTDRMedium(patElectron->leptonID("ptdrMedium"));
                        localElectron.setIDPTDRTight(patElectron->leptonID("ptdrTight"));
                        localElectron.setIDCutBasedLoose(patElectron->leptonID("loose"));
                        localElectron.setIDCutBasedRobust(patElectron->leptonID("robust"));
                        localElectron.setIDCutBasedTight(patElectron->leptonID("tight"));
                        localElectron.setIDLikelihood(patElectron->leptonID("likelihood"));
                        localElectron.setIDNeuralNet(patElectron->leptonID("neuralnet"));

                        // Matched genParticle
                        if(useMC_)
                        {
                                // MC truth associator index
                                if ((patElectron->genParticleRef()).isNonnull()) {
                                        localElectron.setElectronMCIndex((patElectron->genParticleRef()).index());
                                } else {
                                        localElectron.setElectronMCIndex(-1);
                                }

                                if ( patElectron->genLepton() != NULL )
                                {
                                        localElectron.setMomentumMCElectron(patElectron->genLepton()->px(),patElectron->genLepton()->py(),patElectron->genLepton()->pz(),patElectron->genLepton()->energy());
                                        localElectron.setVertexMCElectron(TVector3(patElectron->genLepton()->vx(),patElectron->genLepton()->vy(),patElectron->genLepton()->vz()));
                                        localElectron.setPdgIdMCElectron(patElectron->genLepton()->pdgId());
                                }
                        }

                }

                new( (*rootElectrons)[j] ) TRootElectron(localElectron);
                if(verbosity_>2) cout << "   ["<< setw(3) << j << "] " << localElectron << endl;
        }
}
Revision:	1.3
Committed:	Mon Dec 15 19:08:52 2008 UTC (16 years, 4 months ago) by lethuill
Content type:	text/plain
Branch:	MAIN
CVS Tags:	pat_2_1_12_01
Changes since 1.2:	+188 -105 lines
Log Message:	MC infos from Silvano
#	User	Rev	Content
1	lethuill	1.2	#include "../interface/ElectronAnalyzer.h"
2	lethuill	1.1
3			using namespace std;
4			using namespace reco;
5			using namespace edm;
6
7	lethuill	1.3	ElectronAnalyzer::ElectronAnalyzer(const edm::ParameterSet& producersNames):verbosity_(0),useMC_(false)
8	lethuill	1.1	{
9	lethuill	1.2	dataType_ = producersNames.getUntrackedParameter<string>("dataType","unknown");
10	lethuill	1.1	electronProducer_ = producersNames.getParameter<edm::InputTag>("electronProducer");
11			}
12
13	lethuill	1.3	ElectronAnalyzer::ElectronAnalyzer(const edm::ParameterSet& producersNames, const edm::ParameterSet& myConfig, int verbosity):verbosity_(verbosity)
14	lethuill	1.1	{
15	lethuill	1.2	dataType_ = producersNames.getUntrackedParameter<string>("dataType","unknown");
16	lethuill	1.1	electronProducer_ = producersNames.getParameter<edm::InputTag>("electronProducer");
17	lethuill	1.3	useMC_ = myConfig.getUntrackedParameter<bool>("doElectronMC");
18	lethuill	1.1	}
19
20			ElectronAnalyzer::~ElectronAnalyzer()
21			{
22			}
23
24	lethuill	1.2	void ElectronAnalyzer::Process(const edm::Event& iEvent, TClonesArray* rootElectrons, EcalClusterLazyTools& lazyTools )
25	lethuill	1.1	{
26
27	lethuill	1.3	unsigned int nElectrons=0;
28
29			edm::Handle < std::vector <reco::GsfElectron> > recoElectrons;
30			if( dataType_=="RECO" \|\| dataType_=="AOD" )
31	lethuill	1.1	{
32	lethuill	1.3	iEvent.getByLabel(electronProducer_, recoElectrons);
33			nElectrons = recoElectrons->size();
34	lethuill	1.1	}
35
36	lethuill	1.3	edm::Handle < std::vector <pat::Electron> > patElectrons;
37			if( dataType_=="PAT" \|\| dataType_=="PATAOD" )
38			{
39			iEvent.getByLabel(electronProducer_, patElectrons);
40			nElectrons = patElectrons->size();
41			}
42
43			if(verbosity_>1) std::cout << " Number of electrons = " << nElectrons << " Label: " << electronProducer_.label() << " Instance: " << electronProducer_.instance() << std::endl;
44
45			for (unsigned int j=0; j<nElectrons; j++)
46			{
47			const reco::GsfElectron* electron = 0;
48			if( dataType_=="RECO" \|\| dataType_=="AOD" ) electron = &((*recoElectrons)[j]);
49			if( dataType_=="PAT" \|\| dataType_=="PATAOD" ) electron = (const reco::GsfElectron) ( & ((patElectrons)[j]) );
50
51			TRootElectron localElectron(
52			electron->px()
53			,electron->py()
54			,electron->pz()
55			,electron->energy()
56			,electron->vx()
57			,electron->vy()
58			,electron->vz()
59			,abs(electron->pdgId())
60			,electron->charge()
61			);
62
63			// Variables from reco::GsfElectron
64			localElectron.setClassification(electron->classification());
65			localElectron.setCaloEnergy(electron->caloEnergy());
66			localElectron.setCaloEnergyError(electron->caloEnergyError());
67			// FIXME - Protect against sqrt(0)
68			localElectron.setTrackMomentum(sqrt(electron->trackMomentumAtVtx().Mag2()));
69			localElectron.setTrackMomentumError(electron->trackMomentumError());
70			localElectron.setHadOverEm(electron->hadronicOverEm());
71			localElectron.setDeltaEtaIn(electron->deltaEtaSuperClusterTrackAtVtx());
72			localElectron.setDeltaPhiIn(electron->deltaPhiSuperClusterTrackAtVtx());
73			localElectron.setEnergySuperClusterOverPin(electron->eSuperClusterOverP());
74			localElectron.setDeltaEtaOut(electron->deltaEtaSeedClusterTrackAtCalo());
75			localElectron.setDeltaPhiOut(electron->deltaPhiSeedClusterTrackAtCalo());
76			localElectron.setEnergySeedClusterOverPout(electron->eSeedClusterOverPout());
77			localElectron.setEnergyScaleCorrected(electron->isEnergyScaleCorrected());
78			localElectron.setMomentumCorrected(electron->isMomentumCorrected());
79
80			// Variables from reco::GsfTrack
81			reco::GsfTrackRef gsfTrack = electron->gsfTrack();
82			if ( gsfTrack.isNonnull() )
83			{
84			localElectron.setD0(gsfTrack->d0());
85			localElectron.setD0Error(gsfTrack->d0Error());
86			localElectron.setDsz(gsfTrack->dsz());
87			localElectron.setDszError(gsfTrack->dszError());
88			localElectron.setNormalizedChi2(gsfTrack->normalizedChi2());
89			localElectron.setPtError(gsfTrack->ptError());
90			localElectron.setEtaError(gsfTrack->etaError());
91			localElectron.setPhiError(gsfTrack->phiError());
92			}
93
94			// Variables from reco::SuperCluster
95			reco::SuperClusterRef superCluster = electron->superCluster();
96			//const reco::SuperCluster* sc = superCluster.get();
97			if ( superCluster.isNonnull() )
98			{
99			localElectron.setNbClusters(electron->numberOfClusters());
100			localElectron.setSuperClusterRawEnergy(superCluster->rawEnergy());
101			localElectron.setPreshowerEnergy(superCluster->preshowerEnergy());
102			localElectron.setCaloPosition(
103			superCluster->position().X()
104			,superCluster->position().Y()
105			,superCluster->position().Z()
106			);
107
108			// FIXME - associator supercluster <-> electron
109			//localElectron.setSCRef(superCluster->toto());
110			}
111
112
113			// Cluster Shape variables
114			// need reco::SuperCluster and reco::BasicCluster
115			if ( superCluster.isNonnull() )
116			{
117			reco::BasicClusterRef seedBasicCluster = superCluster->seed();
118			if ( seedBasicCluster.isNonnull() ) localElectron.setClusterAlgo(seedBasicCluster->algo());
119
120			// dR of the cone centered on the reco::GsfElectron and containing all its basic clusters constituents
121			Float_t caloConeSize = 0;
122			for (reco::basicCluster_iterator basicCluster = (superCluster).clustersBegin(); basicCluster != (superCluster).clustersEnd(); ++basicCluster )
123			{
124			Float_t dR = localElectron.DeltaR(TLorentzVector( (basicCluster)->position().x(), (basicCluster)->position().y(), (*basicCluster)->position().z(), 0. ) );
125			if (dR > caloConeSize) caloConeSize = dR;
126			}
127			// FIXME - if no BasicCluster collection, init to -999.
128			localElectron.setCaloConeSize(caloConeSize);
129
130			// need reduced Ecal RecHits Collections for EcalClusterLazyTools
131			// FIXME - check lazyTools is initialized ==> lazyTools has to be a pointer
132			if ( seedBasicCluster.isNonnull() )
133			{
134			localElectron.setE2x2(lazyTools.e2x2(*seedBasicCluster));
135			localElectron.setE3x3(lazyTools.e3x3(*seedBasicCluster));
136			localElectron.setE5x5(lazyTools.e5x5(*seedBasicCluster));
137			localElectron.setEMax(lazyTools.eMax(*seedBasicCluster));
138			}
139			}
140
141
142			if( dataType_=="RECO" \|\| dataType_=="AOD" )
143			{
144			// Some specific methods requiring RECO / AOD format
145			// Do association to genParticle ?
146
147			// Isolation ?
148			// Electron ID ? RecoEgamma/ElectronIdentification
149			}
150
151
152			if( dataType_=="PATAOD" \|\| dataType_=="PAT" )
153			{
154			// Some specific methods to pat::Electron
155			const pat::Electron patElectron = dynamic_cast<const pat::Electron>(&*electron);
156
157			// Isolation
158			// FIXME - how isoDeposit are calculated ? need to substract supercluster energy ?
159			pair < Float_t, Int_t > trackerIso;
160
161			trackerIso = patElectron->trackerIsoDeposit()->depositAndCountWithin(0.1);
162			localElectron.setIsoR01_sumPt(trackerIso.first - gsfTrack->p() );
163			localElectron.setIsoR01_nTracks(trackerIso.second);
164
165			trackerIso = patElectron->trackerIsoDeposit()->depositAndCountWithin(0.2);
166			localElectron.setIsoR02_sumPt(trackerIso.first - gsfTrack->p() );
167			localElectron.setIsoR02_nTracks(trackerIso.second);
168
169			trackerIso = patElectron->trackerIsoDeposit()->depositAndCountWithin(0.3);
170			localElectron.setIsoR03_sumPt(trackerIso.first - gsfTrack->p() );
171			localElectron.setIsoR03_nTracks(trackerIso.second);
172			localElectron.setIsoR03_emEt( patElectron->ecalIsoDeposit()->depositAndCountWithin(0.3).first - superCluster->rawEnergy() );
173			localElectron.setIsoR03_hadEt( patElectron->hcalIsoDeposit()->depositAndCountWithin(0.3).first );
174
175			trackerIso = patElectron->trackerIsoDeposit()->depositAndCountWithin(0.5);
176			localElectron.setIsoR05_sumPt(trackerIso.first - gsfTrack->p() );
177			localElectron.setIsoR05_nTracks(trackerIso.second);
178			localElectron.setIsoR05_emEt( patElectron->ecalIsoDeposit()->depositAndCountWithin(0.5).first - superCluster->rawEnergy() );
179			localElectron.setIsoR05_hadEt( patElectron->hcalIsoDeposit()->depositAndCountWithin(0.5).first );
180
181
182			// Electron ID
183			localElectron.setIDPTDRLoose(patElectron->leptonID("ptdrLoose"));
184			localElectron.setIDPTDRMedium(patElectron->leptonID("ptdrMedium"));
185			localElectron.setIDPTDRTight(patElectron->leptonID("ptdrTight"));
186			localElectron.setIDCutBasedLoose(patElectron->leptonID("loose"));
187			localElectron.setIDCutBasedRobust(patElectron->leptonID("robust"));
188			localElectron.setIDCutBasedTight(patElectron->leptonID("tight"));
189			localElectron.setIDLikelihood(patElectron->leptonID("likelihood"));
190			localElectron.setIDNeuralNet(patElectron->leptonID("neuralnet"));
191
192			// Matched genParticle
193			if(useMC_)
194			{
195			// MC truth associator index
196			if ((patElectron->genParticleRef()).isNonnull()) {
197			localElectron.setElectronMCIndex((patElectron->genParticleRef()).index());
198			} else {
199			localElectron.setElectronMCIndex(-1);
200			}
201
202			if ( patElectron->genLepton() != NULL )
203			{
204			localElectron.setMomentumMCElectron(patElectron->genLepton()->px(),patElectron->genLepton()->py(),patElectron->genLepton()->pz(),patElectron->genLepton()->energy());
205			localElectron.setVertexMCElectron(TVector3(patElectron->genLepton()->vx(),patElectron->genLepton()->vy(),patElectron->genLepton()->vz()));
206			localElectron.setPdgIdMCElectron(patElectron->genLepton()->pdgId());
207			}
208			}
209
210			}
211
212			new( (*rootElectrons)[j] ) TRootElectron(localElectron);
213			if(verbosity_>2) cout << " ["<< setw(3) << j << "] " << localElectron << endl;
214			}
215	lethuill	1.1	}