VHbbAnalysisCode/src/MuonInfoPlotSet.cpp

#include "TrkUpgradeAnalysis/VHbb/interface/MuonInfoPlotSet.h"

#include <stdexcept>

#include <TDirectory.h>
#include <TH1F.h>
#include <TTree.h>

#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"

trkupgradeanalysis::MuonInfoPlotSet::MuonInfoPlotSet( bool createIsolationTree )
        : histogramHaveBeenBooked_(false), createIsolationTree_(createIsolationTree), pIsolationTree_(NULL)
{
        // No operation besides the initialiser list.
}

void trkupgradeanalysis::MuonInfoPlotSet::book( TDirectory* pDirectory )
{
        if( histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::MuonInfoPlotSet::book() - histograms have already been booked" );

        //
        // Note that the root file which TDirectory is part of takes ownership of all
        // of these objects, so I don't need to (and shouldn't) delete them when I'm
        // finished.
        //


        pGlobalChi2_=new TH1F( "globalChi2","Global chi2", 60, 0, 20 );
        pGlobalChi2_->SetDirectory(pDirectory);

        pNumberOfPixelHits_=new TH1F( "numberOfPixelHits", "Number of pixel hits", 11, -0.5, 10.5 );
        pNumberOfPixelHits_->SetDirectory(pDirectory);

        pNumberOfGlobalHits_=new TH1F( "numberOfGlobalHits", "Number of global hits", 31, -0.5, 30.5 );
        pNumberOfGlobalHits_->SetDirectory(pDirectory);

        pNumberOfHits_=new TH1F( "numberOfHits", "Number of hits", 31, -0.5, 30.5 );
        pNumberOfHits_->SetDirectory(pDirectory);

        pCategory_=new TH1F( "category", "Category", 60, 0, 20 );
        pCategory_->SetDirectory(pDirectory);

        pNumberOfMatches_=new TH1F( "numberOfMatches", "Number of matches", 16, -0.5, 15.5 );
        pNumberOfMatches_->SetDirectory(pDirectory);

        pIPDB_=new TH1F( "ipDb", "ipDb, whatever that means", 60, 0, 2 );
        pIPDB_->SetDirectory(pDirectory);

        pEta_=new TH1F( "eta", "Eta", 60, -3, 3 );
        pEta_->SetDirectory(pDirectory);

        pPt_=new TH1F( "pT", "pT", 60, 0, 250 );
        pPt_->SetDirectory(pDirectory);

        pChargedIsolation_=new TH1F( "chargedIsolation", "chargedIsolation", 60, 0, 30 );
        pChargedIsolation_->SetDirectory(pDirectory);

        pPhotonIsolation_=new TH1F( "photonIsolation", "photonIsolation", 60, 0, 30 );
        pPhotonIsolation_->SetDirectory(pDirectory);

        pNeutralIsolation_=new TH1F( "neutralIsolation", "neutralIsolation", 60, 0, 30 );
        pNeutralIsolation_->SetDirectory(pDirectory);

        pPileupIsolation_=new TH1F( "pileupIsolation", "pileupIsolation", 60, 0, 30 );
        pPileupIsolation_->SetDirectory(pDirectory);

        pRelativeChargedIsolation_=new TH1F( "relativeChargedIsolation", "relativeChargedIsolation", 60, 0, 1 );
        pRelativeChargedIsolation_->SetDirectory(pDirectory);

        pRelativePhotonIsolation_=new TH1F( "relativePhotonIsolation", "relativePhotonIsolation", 60, 0, 1 );
        pRelativePhotonIsolation_->SetDirectory(pDirectory);

        pRelativeNeutralIsolation_=new TH1F( "relativeNeutralIsolation", "relativeNeutralIsolation", 60, 0, 1 );
        pRelativeNeutralIsolation_->SetDirectory(pDirectory);

        pRelativePileupIsolation_=new TH1F( "relativePileupIsolation", "relativePileupIsolation", 60, 0, 1 );
        pRelativePileupIsolation_->SetDirectory(pDirectory);

        pRelativeIsolation_=new TH1F( "relativeIsolation", "relativeIsolation", 120, 0, 2 );
        pRelativeIsolation_->SetDirectory(pDirectory);

        pDeltaBetaCorrectedIsolation_=new TH1F( "deltaBetaCorrectedIsolation", "deltaBetaCorrectedIsolation", 120, 0, 2 );
        pDeltaBetaCorrectedIsolation_->SetDirectory(pDirectory);

        pRho25CorrectedIsolation_=new TH1F( "rho25CorrectedIsolation", "rho25CorrectedIsolation", 120, 0, 2 );
        pRho25CorrectedIsolation_->SetDirectory(pDirectory);

        if( createIsolationTree_ )
        {
                pIsolationTree_=new TTree("isolationTree","Data about the muon isolation");
                pIsolationTree_->Branch("numberOfPrimaryVertices",&numberOfPrimaryVertices_branch_,"numberOfPrimaryVertices/I");
                pIsolationTree_->Branch("chargedIsolation",&chargedIsolation_branch_,"chargedIsolation/F");
                pIsolationTree_->Branch("photonIsolation",&photonIsolation_branch_,"photonIsolation/F");
                pIsolationTree_->Branch("neutralIsolation",&neutralIsolation_branch_,"neutralIsolation/F");
                pIsolationTree_->Branch("pileupIsolation",&pileupIsolation_branch_,"pileupIsolation/F");
                pIsolationTree_->Branch("rho25",&rho25_branch_,"rho25/F");
                pIsolationTree_->Branch("pT",&pT_branch_,"pT/F");
                pIsolationTree_->SetDirectory( pDirectory );
        }
        histogramHaveBeenBooked_=true;
}

void trkupgradeanalysis::MuonInfoPlotSet::fill( const VHbbEvent::MuonInfo& muon, const VHbbEventAuxInfo* pAuxInfo )
{
        if( !histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::MuonInfoPlotSet::book() - histograms have not been booked" );

        // Get the additional info about the event if it's available
        int numberOfPrimaryVertices=-1;
        float rho25=-1;
        if( pAuxInfo )
        {
                numberOfPrimaryVertices=pAuxInfo->pvInfo.nVertices;
                rho25=pAuxInfo->puInfo.rho25;
        }

        pGlobalChi2_->Fill( muon.globChi2 );
        pNumberOfPixelHits_->Fill( muon.nPixelHits );
        pNumberOfGlobalHits_->Fill( muon.globNHits );
        pNumberOfHits_->Fill( muon.nHits );
        pCategory_->Fill( muon.cat );
        pNumberOfMatches_->Fill( muon.nMatches );
        pIPDB_->Fill( muon.ipDb );
        pEta_->Fill( muon.p4.Eta() );
        pPt_->Fill( muon.p4.Pt() );
        pChargedIsolation_->Fill( muon.pfChaIso );
        pPhotonIsolation_->Fill( muon.pfPhoIso );
        pNeutralIsolation_->Fill( muon.pfNeuIso );
        pPileupIsolation_->Fill( muon.pfChaPUIso );
        pRelativeChargedIsolation_->Fill( muon.pfChaIso/muon.p4.Pt() );
        pRelativePhotonIsolation_->Fill( muon.pfPhoIso/muon.p4.Pt() );
        pRelativeNeutralIsolation_->Fill( muon.pfNeuIso/muon.p4.Pt() );
        pRelativePileupIsolation_->Fill( muon.pfChaPUIso/muon.p4.Pt() );
        pRelativeIsolation_->Fill( combinedRelativeIsolation(muon) );
        pDeltaBetaCorrectedIsolation_->Fill( deltaBetaCorrectedIsolation(muon) );
        if( pAuxInfo ) pRho25CorrectedIsolation_->Fill( rho25CorrectedIsolation(muon,rho25) ); // Don't bother filling if I couldn't get the rho25 info

        // If the TTree is not null then isolation data has been requested
        if( pIsolationTree_ )
        {
                numberOfPrimaryVertices_branch_=numberOfPrimaryVertices;
                chargedIsolation_branch_=muon.pfChaIso;
                photonIsolation_branch_=muon.pfPhoIso;
                neutralIsolation_branch_=muon.pfNeuIso;
                pileupIsolation_branch_=muon.pfChaPUIso;
                rho25_branch_=rho25;
                pT_branch_=muon.p4.Pt();

                pIsolationTree_->Fill();
        }
}


float trkupgradeanalysis::MuonInfoPlotSet::combinedRelativeIsolation( const VHbbEvent::MuonInfo& muon )
{
        return (muon.pfChaIso+muon.pfPhoIso+muon.pfNeuIso)/muon.p4.Pt();
}

float trkupgradeanalysis::MuonInfoPlotSet::deltaBetaCorrectedIsolation( const VHbbEvent::MuonInfo& muon, float deltaBetaFactor )
{
        // deltaBetaFactor defaults to -0.5
        float correctedNeutralIsolation=muon.pfPhoIso+muon.pfNeuIso+deltaBetaFactor*muon.pfChaPUIso;
        if( correctedNeutralIsolation<0 ) correctedNeutralIsolation=muon.pfPhoIso+muon.pfNeuIso;
        return (muon.pfChaIso+correctedNeutralIsolation)/muon.p4.Pt();
}

float trkupgradeanalysis::MuonInfoPlotSet::rho25CorrectedIsolation( const VHbbEvent::MuonInfo& muon, float rho25 )
{
        float uncorrectedIsolation=combinedRelativeIsolation(muon);
        float correctedIsolation=uncorrectedIsolation-(M_PI*0.4*0.4*rho25)/muon.p4.Pt();
        if( correctedIsolation<0 ) correctedIsolation=0;

        return correctedIsolation;
}
Revision:	1.2
Committed:	Fri Apr 27 13:52:14 2012 UTC (13 years ago) by grimes
Branch:	MAIN
Changes since 1.1:	+106 -3 lines
Log Message:	Added code to analyse muon isolation.
#	Content
1	#include "TrkUpgradeAnalysis/VHbb/interface/MuonInfoPlotSet.h"
2
3	#include <stdexcept>
4
5	#include <TDirectory.h>
6	#include <TH1F.h>
7	#include <TTree.h>
8
9	#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"
10
11	trkupgradeanalysis::MuonInfoPlotSet::MuonInfoPlotSet( bool createIsolationTree )
12	: histogramHaveBeenBooked_(false), createIsolationTree_(createIsolationTree), pIsolationTree_(NULL)
13	{
14	// No operation besides the initialiser list.
15	}
16
17	void trkupgradeanalysis::MuonInfoPlotSet::book( TDirectory* pDirectory )
18	{
19	if( histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::MuonInfoPlotSet::book() - histograms have already been booked" );
20
21	//
22	// Note that the root file which TDirectory is part of takes ownership of all
23	// of these objects, so I don't need to (and shouldn't) delete them when I'm
24	// finished.
25	//
26
27
28	pGlobalChi2_=new TH1F( "globalChi2","Global chi2", 60, 0, 20 );
29	pGlobalChi2_->SetDirectory(pDirectory);
30
31	pNumberOfPixelHits_=new TH1F( "numberOfPixelHits", "Number of pixel hits", 11, -0.5, 10.5 );
32	pNumberOfPixelHits_->SetDirectory(pDirectory);
33
34	pNumberOfGlobalHits_=new TH1F( "numberOfGlobalHits", "Number of global hits", 31, -0.5, 30.5 );
35	pNumberOfGlobalHits_->SetDirectory(pDirectory);
36
37	pNumberOfHits_=new TH1F( "numberOfHits", "Number of hits", 31, -0.5, 30.5 );
38	pNumberOfHits_->SetDirectory(pDirectory);
39
40	pCategory_=new TH1F( "category", "Category", 60, 0, 20 );
41	pCategory_->SetDirectory(pDirectory);
42
43	pNumberOfMatches_=new TH1F( "numberOfMatches", "Number of matches", 16, -0.5, 15.5 );
44	pNumberOfMatches_->SetDirectory(pDirectory);
45
46	pIPDB_=new TH1F( "ipDb", "ipDb, whatever that means", 60, 0, 2 );
47	pIPDB_->SetDirectory(pDirectory);
48
49	pEta_=new TH1F( "eta", "Eta", 60, -3, 3 );
50	pEta_->SetDirectory(pDirectory);
51
52	pPt_=new TH1F( "pT", "pT", 60, 0, 250 );
53	pPt_->SetDirectory(pDirectory);
54
55	pChargedIsolation_=new TH1F( "chargedIsolation", "chargedIsolation", 60, 0, 30 );
56	pChargedIsolation_->SetDirectory(pDirectory);
57
58	pPhotonIsolation_=new TH1F( "photonIsolation", "photonIsolation", 60, 0, 30 );
59	pPhotonIsolation_->SetDirectory(pDirectory);
60
61	pNeutralIsolation_=new TH1F( "neutralIsolation", "neutralIsolation", 60, 0, 30 );
62	pNeutralIsolation_->SetDirectory(pDirectory);
63
64	pPileupIsolation_=new TH1F( "pileupIsolation", "pileupIsolation", 60, 0, 30 );
65	pPileupIsolation_->SetDirectory(pDirectory);
66
67	pRelativeChargedIsolation_=new TH1F( "relativeChargedIsolation", "relativeChargedIsolation", 60, 0, 1 );
68	pRelativeChargedIsolation_->SetDirectory(pDirectory);
69
70	pRelativePhotonIsolation_=new TH1F( "relativePhotonIsolation", "relativePhotonIsolation", 60, 0, 1 );
71	pRelativePhotonIsolation_->SetDirectory(pDirectory);
72
73	pRelativeNeutralIsolation_=new TH1F( "relativeNeutralIsolation", "relativeNeutralIsolation", 60, 0, 1 );
74	pRelativeNeutralIsolation_->SetDirectory(pDirectory);
75
76	pRelativePileupIsolation_=new TH1F( "relativePileupIsolation", "relativePileupIsolation", 60, 0, 1 );
77	pRelativePileupIsolation_->SetDirectory(pDirectory);
78
79	pRelativeIsolation_=new TH1F( "relativeIsolation", "relativeIsolation", 120, 0, 2 );
80	pRelativeIsolation_->SetDirectory(pDirectory);
81
82	pDeltaBetaCorrectedIsolation_=new TH1F( "deltaBetaCorrectedIsolation", "deltaBetaCorrectedIsolation", 120, 0, 2 );
83	pDeltaBetaCorrectedIsolation_->SetDirectory(pDirectory);
84
85	pRho25CorrectedIsolation_=new TH1F( "rho25CorrectedIsolation", "rho25CorrectedIsolation", 120, 0, 2 );
86	pRho25CorrectedIsolation_->SetDirectory(pDirectory);
87
88	if( createIsolationTree_ )
89	{
90	pIsolationTree_=new TTree("isolationTree","Data about the muon isolation");
91	pIsolationTree_->Branch("numberOfPrimaryVertices",&numberOfPrimaryVertices_branch_,"numberOfPrimaryVertices/I");
92	pIsolationTree_->Branch("chargedIsolation",&chargedIsolation_branch_,"chargedIsolation/F");
93	pIsolationTree_->Branch("photonIsolation",&photonIsolation_branch_,"photonIsolation/F");
94	pIsolationTree_->Branch("neutralIsolation",&neutralIsolation_branch_,"neutralIsolation/F");
95	pIsolationTree_->Branch("pileupIsolation",&pileupIsolation_branch_,"pileupIsolation/F");
96	pIsolationTree_->Branch("rho25",&rho25_branch_,"rho25/F");
97	pIsolationTree_->Branch("pT",&pT_branch_,"pT/F");
98	pIsolationTree_->SetDirectory( pDirectory );
99	}
100	histogramHaveBeenBooked_=true;
101	}
102
103	void trkupgradeanalysis::MuonInfoPlotSet::fill( const VHbbEvent::MuonInfo& muon, const VHbbEventAuxInfo* pAuxInfo )
104	{
105	if( !histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::MuonInfoPlotSet::book() - histograms have not been booked" );
106
107	// Get the additional info about the event if it's available
108	int numberOfPrimaryVertices=-1;
109	float rho25=-1;
110	if( pAuxInfo )
111	{
112	numberOfPrimaryVertices=pAuxInfo->pvInfo.nVertices;
113	rho25=pAuxInfo->puInfo.rho25;
114	}
115
116	pGlobalChi2_->Fill( muon.globChi2 );
117	pNumberOfPixelHits_->Fill( muon.nPixelHits );
118	pNumberOfGlobalHits_->Fill( muon.globNHits );
119	pNumberOfHits_->Fill( muon.nHits );
120	pCategory_->Fill( muon.cat );
121	pNumberOfMatches_->Fill( muon.nMatches );
122	pIPDB_->Fill( muon.ipDb );
123	pEta_->Fill( muon.p4.Eta() );
124	pPt_->Fill( muon.p4.Pt() );
125	pChargedIsolation_->Fill( muon.pfChaIso );
126	pPhotonIsolation_->Fill( muon.pfPhoIso );
127	pNeutralIsolation_->Fill( muon.pfNeuIso );
128	pPileupIsolation_->Fill( muon.pfChaPUIso );
129	pRelativeChargedIsolation_->Fill( muon.pfChaIso/muon.p4.Pt() );
130	pRelativePhotonIsolation_->Fill( muon.pfPhoIso/muon.p4.Pt() );
131	pRelativeNeutralIsolation_->Fill( muon.pfNeuIso/muon.p4.Pt() );
132	pRelativePileupIsolation_->Fill( muon.pfChaPUIso/muon.p4.Pt() );
133	pRelativeIsolation_->Fill( combinedRelativeIsolation(muon) );
134	pDeltaBetaCorrectedIsolation_->Fill( deltaBetaCorrectedIsolation(muon) );
135	if( pAuxInfo ) pRho25CorrectedIsolation_->Fill( rho25CorrectedIsolation(muon,rho25) ); // Don't bother filling if I couldn't get the rho25 info
136
137	// If the TTree is not null then isolation data has been requested
138	if( pIsolationTree_ )
139	{
140	numberOfPrimaryVertices_branch_=numberOfPrimaryVertices;
141	chargedIsolation_branch_=muon.pfChaIso;
142	photonIsolation_branch_=muon.pfPhoIso;
143	neutralIsolation_branch_=muon.pfNeuIso;
144	pileupIsolation_branch_=muon.pfChaPUIso;
145	rho25_branch_=rho25;
146	pT_branch_=muon.p4.Pt();
147
148	pIsolationTree_->Fill();
149	}
150	}
151
152
153	float trkupgradeanalysis::MuonInfoPlotSet::combinedRelativeIsolation( const VHbbEvent::MuonInfo& muon )
154	{
155	return (muon.pfChaIso+muon.pfPhoIso+muon.pfNeuIso)/muon.p4.Pt();
156	}
157
158	float trkupgradeanalysis::MuonInfoPlotSet::deltaBetaCorrectedIsolation( const VHbbEvent::MuonInfo& muon, float deltaBetaFactor )
159	{
160	// deltaBetaFactor defaults to -0.5
161	float correctedNeutralIsolation=muon.pfPhoIso+muon.pfNeuIso+deltaBetaFactor*muon.pfChaPUIso;
162	if( correctedNeutralIsolation<0 ) correctedNeutralIsolation=muon.pfPhoIso+muon.pfNeuIso;
163	return (muon.pfChaIso+correctedNeutralIsolation)/muon.p4.Pt();
164	}
165
166	float trkupgradeanalysis::MuonInfoPlotSet::rho25CorrectedIsolation( const VHbbEvent::MuonInfo& muon, float rho25 )
167	{
168	float uncorrectedIsolation=combinedRelativeIsolation(muon);
169	float correctedIsolation=uncorrectedIsolation-(M_PI0.40.4*rho25)/muon.p4.Pt();
170	if( correctedIsolation<0 ) correctedIsolation=0;
171
172	return correctedIsolation;
173	}