Morgan/src/MuonAnalyzer.cc

#include "../interface/MuonAnalyzer.h"

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

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

MuonAnalyzer::MuonAnalyzer(const edm::ParameterSet& producersNames, const edm::ParameterSet& myConfig, int verbosity):verbosity_(verbosity)
{
        dataType_ = producersNames.getUntrackedParameter<string>("dataType","unknown");
        muonProducer_ = producersNames.getParameter<edm::InputTag>("muonProducer");
        useMC_ = myConfig.getUntrackedParameter<bool>("doMuonMC");
}

MuonAnalyzer::~MuonAnalyzer()
{
}

void MuonAnalyzer::Process(const edm::Event& iEvent, TClonesArray* rootMuons)
{

        Float_t sintheta = 0.;
        unsigned int nMuons=0;

        edm::Handle < std::vector <reco::Muon> > recoMuons;
        if( dataType_=="RECO" || dataType_=="AOD" )
        {
                iEvent.getByLabel(muonProducer_, recoMuons);
                nMuons = recoMuons->size();
        }

        edm::Handle < std::vector <pat::Muon> > patMuons;
        if( dataType_=="PAT" || dataType_=="PATAOD" )
        {
                iEvent.getByLabel(muonProducer_, patMuons);
                nMuons = patMuons->size();
        }

        if(verbosity_>1) std::cout << "   Number of muons = " << nMuons << "   Label: " << muonProducer_.label() << "   Instance: " << muonProducer_.instance() << std::endl;

        for (unsigned int j=0; j<nMuons; j++)
        {
                const reco::Muon* muon = 0;     
                if( dataType_=="RECO" || dataType_=="AOD" ) muon =  &((*recoMuons)[j]);
                if( dataType_=="PAT" || dataType_=="PATAOD" ) muon = (const reco::Muon*) ( & ((*patMuons)[j]) );

                TRootMuon localMuon(
                        muon->px()
                        ,muon->py()
                        ,muon->pz()
                        ,muon->energy()
                        ,muon->vx()
                        ,muon->vy()
                        ,muon->vz()
                        ,muon->pdgId()
                        ,muon->charge()
                ); 

                // Variables from reco::Muon
                sintheta = sin( localMuon.Theta() );

                localMuon.setCaloEnergy(
                        muon->calEnergy().em * sintheta
                        ,muon->calEnergy().emS9 * sintheta
                        ,muon->calEnergy().had * sintheta
                        ,muon->calEnergy().hadS9 * sintheta
                        ,muon->calEnergy().ho * sintheta
                        ,muon->calEnergy().hoS9 * sintheta
                        ,muon->caloCompatibility()
                );

                localMuon.setIsoR03(
                        muon->isolationR03().emEt
                        ,muon->isolationR03().hadEt
                        ,muon->isolationR03().hoEt
                        ,muon->isolationR03().sumPt
                        ,muon->isolationR03().nTracks
                        ,muon->isolationR03().nJets
                );

                localMuon.setIsoR05(
                        muon->isolationR05().emEt
                        ,muon->isolationR05().hadEt
                        ,muon->isolationR05().hoEt
                        ,muon->isolationR05().sumPt
                        ,muon->isolationR05().nTracks
                        ,muon->isolationR05().nJets
                );

                localMuon.setValidity(
                        muon->isEnergyValid()
                        ,muon->isMatchesValid()
                        ,muon->isIsolationValid()
                );

                localMuon.setDirection( muon->time().direction() );
                localMuon.setAlgo( muon->type() );
                localMuon.setID(
                        int(muon->isGood(reco::Muon::TrackerMuonArbitrated))
                        ,int(muon->isGood(reco::Muon::AllArbitrated))
                        ,int(muon->isGood(reco::Muon::GlobalMuonPromptTight))
                        ,int(muon->isGood(reco::Muon::TMLastStationLoose))
                        ,int(muon->isGood(reco::Muon::TMLastStationTight))
                        ,int(muon->isGood(reco::Muon::TM2DCompatibilityLoose))
                        ,int(muon->isGood(reco::Muon::TM2DCompatibilityTight))
                );

                // Variables from reco::Track
                reco::TrackRef track = muon->innerTrack();
                if ( track.isNonnull() )
                {
                        const reco::HitPattern& hit = track->hitPattern();
                        cout << "Muon["<<j<<"] (pt,eta,phi)="<<localMuon.Et()<<","<<localMuon.Eta()<<","<<localMuon.Phi()<<")"
                                << endl << " numberOfValidPixelHits=" << hit.numberOfValidPixelHits()
                                << endl << " numberOfValidPixelBarrelHits=" << hit.numberOfValidPixelBarrelHits()
                                << endl << " numberOfValidPixelEndcapHits=" << hit.numberOfValidPixelEndcapHits()
                                << endl << " numberOfValidTrackerHits=" << hit.numberOfValidTrackerHits()
                                << endl << " numberOfValidStripHits=" << hit.numberOfValidStripHits()
                                << endl << " numberOfValidStripTIBHits=" << hit.numberOfValidStripTIBHits()
                                << endl << " numberOfValidStripTIDHits=" << hit.numberOfValidStripTIDHits()
                                << endl << " numberOfValidStripTOBHits=" << hit.numberOfValidStripTOBHits()
                                << endl << " numberOfValidStripTECHits=" << hit.numberOfValidStripTECHits()
                                << endl << " numberOfValidMuonHits=" << hit.numberOfValidMuonHits()
                                << endl << " numberOfValidMuonDTHits=" << hit.numberOfValidMuonDTHits()
                                << endl << " numberOfValidMuonCSCHits=" << hit.numberOfValidMuonCSCHits()
                                << endl << " numberOfValidMuonRPCHits=" << hit.numberOfValidMuonRPCHits()
                                << endl << " trackerLayersWithMeasurement=" << hit.trackerLayersWithMeasurement()
                                << endl << " pixelLayersWithMeasurement=" << hit.pixelLayersWithMeasurement()
                                << endl << " stripLayersWithMeasurement=" << hit.stripLayersWithMeasurement()
                                << endl << " pixelBarrelLayersWithMeasurement=" << hit.pixelBarrelLayersWithMeasurement()
                                << endl << " pixelEndcapLayersWithMeasurement=" << hit.pixelEndcapLayersWithMeasurement()
                                << endl << " stripTIBLayersWithMeasurement=" << hit.stripTIBLayersWithMeasurement()
                                << endl << " stripTIDLayersWithMeasurement=" << hit.stripTIDLayersWithMeasurement()
                                << endl << " stripTOBLayersWithMeasurement=" << hit.stripTOBLayersWithMeasurement()
                                << endl << " stripTECLayersWithMeasurement=" << hit.stripTECLayersWithMeasurement()
                                << endl;

                        localMuon.setPixelLayersWithMeasurement(hit.pixelLayersWithMeasurement());
                        localMuon.setStripLayersWithMeasurement(hit.stripLayersWithMeasurement());
                        // FIXME - replace with double dxy(const Point& myBeamSpot) method ?
                        localMuon.setD0(track->d0());
                        localMuon.setD0Error(track->d0Error());
                        // FIXME - replace with double dsz(const Point& myBeamSpot) method ?
                        localMuon.setDsz(track->dsz());
                        localMuon.setDszError(track->dszError());
                        localMuon.setNormalizedChi2(track->normalizedChi2());
                        localMuon.setPtError(track->ptError());
                        localMuon.setEtaError(track->etaError());
                        localMuon.setPhiError(track->phiError());
                }


                if( dataType_=="RECO" || dataType_=="AOD" )
                {
                        // Some specific methods requiring  RECO / AOD format
                        // Do association to genParticle ?
                        // Add InnerTrack, OuterTrack, GlobalTrack infos ?
                }

                if( dataType_=="PATAOD" || dataType_=="PAT" )
                {
                        // Some specific methods to pat::Muon
                        const pat::Muon *patMuon = dynamic_cast<const pat::Muon*>(&*muon);
                        // Use existing reference to genParticle [ pat::PATObject::genParticleRef() ] ?
                        // Alternative methode for isolation (isoDeposit) ?
                        //
                        // leptonID apparently not initialised in PAT...
                        // cout << "Valeur pourrie du leptonID=" << patMuon->leptonID() << endl;

                        if(useMC_)
                        {
                                // MC truth associator index
                                if ((patMuon->genParticleRef()).isNonnull()) {
                                        localMuon.setGenParticleIndex((patMuon->genParticleRef()).index());
                                } else {
                                        localMuon.setGenParticleIndex(-1);
                                }
                        }
                }

                new( (*rootMuons)[j] ) TRootMuon(localMuon);
                if(verbosity_>2) cout << "   ["<< setw(3) << j << "] " << localMuon << endl;
        }
}
Revision:	1.8
Committed:	Tue Apr 14 16:50:16 2009 UTC (16 years ago) by lethuill
Content type:	text/plain
Branch:	MAIN
Changes since 1.7:	+27 -2 lines
Log Message:	Replace number of valid pixel and tracker hits with number of pixel and strip layers with measurements
#	User	Rev	Content
1	lethuill	1.2	#include "../interface/MuonAnalyzer.h"
2	lethuill	1.1
3			using namespace std;
4			using namespace reco;
5			using namespace edm;
6
7	lethuill	1.4	MuonAnalyzer::MuonAnalyzer(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	muonProducer_ = producersNames.getParameter<edm::InputTag>("muonProducer");
11			}
12
13	lethuill	1.4	MuonAnalyzer::MuonAnalyzer(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	muonProducer_ = producersNames.getParameter<edm::InputTag>("muonProducer");
17	lethuill	1.4	useMC_ = myConfig.getUntrackedParameter<bool>("doMuonMC");
18	lethuill	1.1	}
19
20			MuonAnalyzer::~MuonAnalyzer()
21			{
22			}
23
24			void MuonAnalyzer::Process(const edm::Event& iEvent, TClonesArray* rootMuons)
25			{
26
27			Float_t sintheta = 0.;
28	lethuill	1.2	unsigned int nMuons=0;
29
30			edm::Handle < std::vector <reco::Muon> > recoMuons;
31			if( dataType_=="RECO" \|\| dataType_=="AOD" )
32			{
33			iEvent.getByLabel(muonProducer_, recoMuons);
34			nMuons = recoMuons->size();
35			}
36
37			edm::Handle < std::vector <pat::Muon> > patMuons;
38			if( dataType_=="PAT" \|\| dataType_=="PATAOD" )
39			{
40			iEvent.getByLabel(muonProducer_, patMuons);
41			nMuons = patMuons->size();
42			}
43
44			if(verbosity_>1) std::cout << " Number of muons = " << nMuons << " Label: " << muonProducer_.label() << " Instance: " << muonProducer_.instance() << std::endl;
45
46			for (unsigned int j=0; j<nMuons; j++)
47	lethuill	1.1	{
48	lethuill	1.3	const reco::Muon* muon = 0;
49			if( dataType_=="RECO" \|\| dataType_=="AOD" ) muon = &((*recoMuons)[j]);
50			if( dataType_=="PAT" \|\| dataType_=="PATAOD" ) muon = (const reco::Muon) ( & ((patMuons)[j]) );
51	lethuill	1.2
52	lethuill	1.1	TRootMuon localMuon(
53	lethuill	1.3	muon->px()
54			,muon->py()
55			,muon->pz()
56			,muon->energy()
57			,muon->vx()
58			,muon->vy()
59			,muon->vz()
60	lethuill	1.6	,muon->pdgId()
61	lethuill	1.3	,muon->charge()
62	lethuill	1.1	);
63
64	lethuill	1.4	// Variables from reco::Muon
65	lethuill	1.1	sintheta = sin( localMuon.Theta() );
66	lethuill	1.2
67	lethuill	1.3	localMuon.setCaloEnergy(
68			muon->calEnergy().em * sintheta
69			,muon->calEnergy().emS9 * sintheta
70			,muon->calEnergy().had * sintheta
71			,muon->calEnergy().hadS9 * sintheta
72			,muon->calEnergy().ho * sintheta
73			,muon->calEnergy().hoS9 * sintheta
74			,muon->caloCompatibility()
75			);
76
77			localMuon.setIsoR03(
78			muon->isolationR03().emEt
79			,muon->isolationR03().hadEt
80			,muon->isolationR03().hoEt
81			,muon->isolationR03().sumPt
82			,muon->isolationR03().nTracks
83			,muon->isolationR03().nJets
84			);
85
86			localMuon.setIsoR05(
87			muon->isolationR05().emEt
88			,muon->isolationR05().hadEt
89			,muon->isolationR05().hoEt
90			,muon->isolationR05().sumPt
91			,muon->isolationR05().nTracks
92			,muon->isolationR05().nJets
93			);
94
95			localMuon.setValidity(
96			muon->isEnergyValid()
97			,muon->isMatchesValid()
98			,muon->isIsolationValid()
99			);
100
101			localMuon.setDirection( muon->time().direction() );
102			localMuon.setAlgo( muon->type() );
103			localMuon.setID(
104			int(muon->isGood(reco::Muon::TrackerMuonArbitrated))
105			,int(muon->isGood(reco::Muon::AllArbitrated))
106			,int(muon->isGood(reco::Muon::GlobalMuonPromptTight))
107			,int(muon->isGood(reco::Muon::TMLastStationLoose))
108			,int(muon->isGood(reco::Muon::TMLastStationTight))
109			,int(muon->isGood(reco::Muon::TM2DCompatibilityLoose))
110			,int(muon->isGood(reco::Muon::TM2DCompatibilityTight))
111			);
112
113	lethuill	1.7	// Variables from reco::Track
114			reco::TrackRef track = muon->innerTrack();
115			if ( track.isNonnull() )
116			{
117			const reco::HitPattern& hit = track->hitPattern();
118	lethuill	1.8	cout << "Muon["<<j<<"] (pt,eta,phi)="<<localMuon.Et()<<","<<localMuon.Eta()<<","<<localMuon.Phi()<<")"
119			<< endl << " numberOfValidPixelHits=" << hit.numberOfValidPixelHits()
120			<< endl << " numberOfValidPixelBarrelHits=" << hit.numberOfValidPixelBarrelHits()
121			<< endl << " numberOfValidPixelEndcapHits=" << hit.numberOfValidPixelEndcapHits()
122			<< endl << " numberOfValidTrackerHits=" << hit.numberOfValidTrackerHits()
123			<< endl << " numberOfValidStripHits=" << hit.numberOfValidStripHits()
124			<< endl << " numberOfValidStripTIBHits=" << hit.numberOfValidStripTIBHits()
125			<< endl << " numberOfValidStripTIDHits=" << hit.numberOfValidStripTIDHits()
126			<< endl << " numberOfValidStripTOBHits=" << hit.numberOfValidStripTOBHits()
127			<< endl << " numberOfValidStripTECHits=" << hit.numberOfValidStripTECHits()
128			<< endl << " numberOfValidMuonHits=" << hit.numberOfValidMuonHits()
129			<< endl << " numberOfValidMuonDTHits=" << hit.numberOfValidMuonDTHits()
130			<< endl << " numberOfValidMuonCSCHits=" << hit.numberOfValidMuonCSCHits()
131			<< endl << " numberOfValidMuonRPCHits=" << hit.numberOfValidMuonRPCHits()
132			<< endl << " trackerLayersWithMeasurement=" << hit.trackerLayersWithMeasurement()
133			<< endl << " pixelLayersWithMeasurement=" << hit.pixelLayersWithMeasurement()
134			<< endl << " stripLayersWithMeasurement=" << hit.stripLayersWithMeasurement()
135			<< endl << " pixelBarrelLayersWithMeasurement=" << hit.pixelBarrelLayersWithMeasurement()
136			<< endl << " pixelEndcapLayersWithMeasurement=" << hit.pixelEndcapLayersWithMeasurement()
137			<< endl << " stripTIBLayersWithMeasurement=" << hit.stripTIBLayersWithMeasurement()
138			<< endl << " stripTIDLayersWithMeasurement=" << hit.stripTIDLayersWithMeasurement()
139			<< endl << " stripTOBLayersWithMeasurement=" << hit.stripTOBLayersWithMeasurement()
140			<< endl << " stripTECLayersWithMeasurement=" << hit.stripTECLayersWithMeasurement()
141			<< endl;
142
143			localMuon.setPixelLayersWithMeasurement(hit.pixelLayersWithMeasurement());
144			localMuon.setStripLayersWithMeasurement(hit.stripLayersWithMeasurement());
145	lethuill	1.7	// FIXME - replace with double dxy(const Point& myBeamSpot) method ?
146			localMuon.setD0(track->d0());
147			localMuon.setD0Error(track->d0Error());
148			// FIXME - replace with double dsz(const Point& myBeamSpot) method ?
149			localMuon.setDsz(track->dsz());
150			localMuon.setDszError(track->dszError());
151			localMuon.setNormalizedChi2(track->normalizedChi2());
152			localMuon.setPtError(track->ptError());
153			localMuon.setEtaError(track->etaError());
154			localMuon.setPhiError(track->phiError());
155			}
156
157
158	lethuill	1.3	if( dataType_=="RECO" \|\| dataType_=="AOD" )
159	lethuill	1.2	{
160	lethuill	1.3	// Some specific methods requiring RECO / AOD format
161			// Do association to genParticle ?
162			// Add InnerTrack, OuterTrack, GlobalTrack infos ?
163	lethuill	1.2	}
164
165	lethuill	1.3	if( dataType_=="PATAOD" \|\| dataType_=="PAT" )
166	lethuill	1.2	{
167			// Some specific methods to pat::Muon
168			const pat::Muon patMuon = dynamic_cast<const pat::Muon>(&*muon);
169	lethuill	1.3	// Use existing reference to genParticle [ pat::PATObject::genParticleRef() ] ?
170			// Alternative methode for isolation (isoDeposit) ?
171			//
172			// leptonID apparently not initialised in PAT...
173			// cout << "Valeur pourrie du leptonID=" << patMuon->leptonID() << endl;
174	lethuill	1.4
175			if(useMC_)
176			{
177			// MC truth associator index
178			if ((patMuon->genParticleRef()).isNonnull()) {
179	lethuill	1.5	localMuon.setGenParticleIndex((patMuon->genParticleRef()).index());
180	lethuill	1.4	} else {
181	lethuill	1.5	localMuon.setGenParticleIndex(-1);
182	lethuill	1.4	}
183			}
184	lethuill	1.2	}
185	lethuill	1.6
186	lethuill	1.1	new( (*rootMuons)[j] ) TRootMuon(localMuon);
187			if(verbosity_>2) cout << " ["<< setw(3) << j << "] " << localMuon << endl;
188			}
189			}