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.
#	User	Rev	Content
1	grimes	1.1	#include "TrkUpgradeAnalysis/VHbb/interface/MuonInfoPlotSet.h"
2
3			#include <stdexcept>
4
5			#include <TDirectory.h>
6			#include <TH1F.h>
7	grimes	1.2	#include <TTree.h>
8	grimes	1.1
9	grimes	1.2	#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"
10	grimes	1.1
11	grimes	1.2	trkupgradeanalysis::MuonInfoPlotSet::MuonInfoPlotSet( bool createIsolationTree )
12			: histogramHaveBeenBooked_(false), createIsolationTree_(createIsolationTree), pIsolationTree_(NULL)
13	grimes	1.1	{
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	grimes	1.2	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	grimes	1.1	histogramHaveBeenBooked_=true;
101			}
102
103	grimes	1.2	void trkupgradeanalysis::MuonInfoPlotSet::fill( const VHbbEvent::MuonInfo& muon, const VHbbEventAuxInfo* pAuxInfo )
104	grimes	1.1	{
105			if( !histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::MuonInfoPlotSet::book() - histograms have not been booked" );
106
107	grimes	1.2	// 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	grimes	1.1	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	grimes	1.2	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	grimes	1.1	}
151
152	grimes	1.2
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			}