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()
                        ,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
                        pair < Float_t, Int_t > trackerIso;

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

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

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

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


                        // Electron ID
                        
                        // Old 2.1.X 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"));
                        */

                        // New 2.2.X electron ID
                        // Only Cut Based ID available by default (4 sequential cuts on H/E, DeltaEta, DeltaPhi, SigmaEtaEta)
                        // "Robust" ids (eidRobustLoose, eidRobustTight, eidRobustHighEnergy) corresponds to fixed threshold
                        // eidLoose and eidTight corresponds to the catagory based identification (E/p, fBrem)
                        localElectron.setIDCutBasedFixedThresholdLoose(int(patElectron->electronID("eidRobustLoose")));
                        localElectron.setIDCutBasedFixedThresholdTight(int(patElectron->electronID("eidRobustTight")));
                        localElectron.setIDCutBasedFixedThresholdHighEnergy(int(patElectron->electronID("eidRobustHighEnergy")));
                        localElectron.setIDCutBasedCategorizedLoose(int(patElectron->electronID("eidLoose")));
                        localElectron.setIDCutBasedCategorizedTight(int(patElectron->electronID("eidTight")));
                        //localElectron.setIDLikelihood(patElectron->leptonID("likelihood"));
                        //localElectron.setIDNeuralNet(patElectron->leptonID("neuralnet"));


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

                }

                new( (*rootElectrons)[j] ) TRootElectron(localElectron);
                if(verbosity_>2) cout << "   ["<< setw(3) << j << "] " << localElectron << endl;
        }
}
Revision:	1.7
Committed:	Wed Mar 11 12:44:56 2009 UTC (16 years, 1 month ago) by lethuill
Content type:	text/plain
Branch:	MAIN
CVS Tags:	pat_2_2_5_02, pat_2_2_5_01
Changes since 1.6:	+17 -0 lines
Log Message:	Transition to 2.2.X
#	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	lethuill	1.5	,electron->pdgId()
60	lethuill	1.3	,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			pair < Float_t, Int_t > trackerIso;
159
160			trackerIso = patElectron->trackerIsoDeposit()->depositAndCountWithin(0.1);
161	lethuill	1.6	localElectron.setIsoR01_sumPt(trackerIso.first);
162	lethuill	1.3	localElectron.setIsoR01_nTracks(trackerIso.second);
163
164			trackerIso = patElectron->trackerIsoDeposit()->depositAndCountWithin(0.2);
165	lethuill	1.6	localElectron.setIsoR02_sumPt(trackerIso.first);
166	lethuill	1.3	localElectron.setIsoR02_nTracks(trackerIso.second);
167
168			trackerIso = patElectron->trackerIsoDeposit()->depositAndCountWithin(0.3);
169	lethuill	1.6	localElectron.setIsoR03_sumPt(trackerIso.first);
170	lethuill	1.3	localElectron.setIsoR03_nTracks(trackerIso.second);
171	lethuill	1.6	localElectron.setIsoR03_emEt( patElectron->ecalIsoDeposit()->depositAndCountWithin(0.3).first );
172	lethuill	1.3	localElectron.setIsoR03_hadEt( patElectron->hcalIsoDeposit()->depositAndCountWithin(0.3).first );
173
174			trackerIso = patElectron->trackerIsoDeposit()->depositAndCountWithin(0.5);
175	lethuill	1.6	localElectron.setIsoR05_sumPt(trackerIso.first );
176	lethuill	1.3	localElectron.setIsoR05_nTracks(trackerIso.second);
177	lethuill	1.6	localElectron.setIsoR05_emEt( patElectron->ecalIsoDeposit()->depositAndCountWithin(0.5).first );
178	lethuill	1.3	localElectron.setIsoR05_hadEt( patElectron->hcalIsoDeposit()->depositAndCountWithin(0.5).first );
179
180
181			// Electron ID
182	lethuill	1.7
183			// Old 2.1.X electron ID
184			/*
185	lethuill	1.3	localElectron.setIDPTDRLoose(patElectron->leptonID("ptdrLoose"));
186			localElectron.setIDPTDRMedium(patElectron->leptonID("ptdrMedium"));
187			localElectron.setIDPTDRTight(patElectron->leptonID("ptdrTight"));
188			localElectron.setIDCutBasedLoose(patElectron->leptonID("loose"));
189			localElectron.setIDCutBasedRobust(patElectron->leptonID("robust"));
190			localElectron.setIDCutBasedTight(patElectron->leptonID("tight"));
191			localElectron.setIDLikelihood(patElectron->leptonID("likelihood"));
192			localElectron.setIDNeuralNet(patElectron->leptonID("neuralnet"));
193	lethuill	1.7	*/
194
195			// New 2.2.X electron ID
196			// Only Cut Based ID available by default (4 sequential cuts on H/E, DeltaEta, DeltaPhi, SigmaEtaEta)
197			// "Robust" ids (eidRobustLoose, eidRobustTight, eidRobustHighEnergy) corresponds to fixed threshold
198			// eidLoose and eidTight corresponds to the catagory based identification (E/p, fBrem)
199			localElectron.setIDCutBasedFixedThresholdLoose(int(patElectron->electronID("eidRobustLoose")));
200			localElectron.setIDCutBasedFixedThresholdTight(int(patElectron->electronID("eidRobustTight")));
201			localElectron.setIDCutBasedFixedThresholdHighEnergy(int(patElectron->electronID("eidRobustHighEnergy")));
202			localElectron.setIDCutBasedCategorizedLoose(int(patElectron->electronID("eidLoose")));
203			localElectron.setIDCutBasedCategorizedTight(int(patElectron->electronID("eidTight")));
204			//localElectron.setIDLikelihood(patElectron->leptonID("likelihood"));
205			//localElectron.setIDNeuralNet(patElectron->leptonID("neuralnet"));
206
207	lethuill	1.3
208			// Matched genParticle
209			if(useMC_)
210			{
211			// MC truth associator index
212			if ((patElectron->genParticleRef()).isNonnull()) {
213	lethuill	1.4	localElectron.setGenParticleIndex((patElectron->genParticleRef()).index());
214	lethuill	1.3	} else {
215	lethuill	1.4	localElectron.setGenParticleIndex(-1);
216	lethuill	1.3	}
217			}
218
219			}
220
221			new( (*rootElectrons)[j] ) TRootElectron(localElectron);
222			if(verbosity_>2) cout << " ["<< setw(3) << j << "] " << localElectron << endl;
223			}
224	lethuill	1.1	}