VHbbAnalysisCode/src/IsolationStudyPlotSet.cpp

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

#include <stdexcept>

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

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

namespace // Use the unnamed namespace for things only used in this file
{
        /** @brief A functor to sort muon collections by pT in the STL sort algorithm.
         * @author Mark Grimes (mark.grimes@bristol.ac.uk)
         * @date 11/Apr/2012
         */
        template <typename T>
        class SortByTransverseMomentum
        {
        public:
                // Return true if the first argument should come before
                // the second, as per normal STL sort comparitors
                bool operator()( const T& first, const T& second )
                {
                        return first.p4.Pt() > second.p4.Pt();
                }
        };

        /** @brief Class to pass as exception when something in the analysis breaks down.
         *
         * Intended to be thrown when the event/particle/whatever cannot be processed further
         * but the analysis as a whole should continue running. A std::runtime_error will be
         * thrown when everything should stop.
         */
        class AnalysisFailureException : public std::runtime_error
        {
        public:
                AnalysisFailureException( const std::string& what ) : std::runtime_error(what) {}
        };
}

trkupgradeanalysis::IsolationStudyPlotSet::IsolationStudyPlotSet()
        : histogramHaveBeenBooked_(false), leastIsolatedDiMuon_(true), leastDeltaBetaCorrectedIsolatedDiMuon_(true)
{
        // No operation besides the initialiser list.
}

void trkupgradeanalysis::IsolationStudyPlotSet::book( TDirectory* pDirectory )
{
        if( histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::IsolationStudyPlotSet::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);

//      pLeastIsolatedDiMuonTree_=new TTree("leastIsolatedDiMuonTree","Data about the least isolated muon in the di-muon pair");
//      pLeastIsolatedDiMuonTree_->Branch("numberOfPrimaryVertices",&numberOfPrimaryVertices_branch_,"numberOfPrimaryVertices/I");
//      pLeastIsolatedDiMuonTree_->Branch("chargedIsolation",&chargedIsolation_branch_,"chargedIsolation/F");
//      pLeastIsolatedDiMuonTree_->Branch("photonIsolation",&photonIsolation_branch_,"photonIsolation/F");
//      pLeastIsolatedDiMuonTree_->Branch("neutralIsolation",&neutralIsolation_branch_,"neutralIsolation/F");
//      pLeastIsolatedDiMuonTree_->Branch("pileupIsolation",&pileupIsolation_branch_,"pileupIsolation/F");
//      pLeastIsolatedDiMuonTree_->Branch("pT",&pT_branch_,"pT/F");
//      pLeastIsolatedDiMuonTree_->SetDirectory( pDirectory );

        TDirectory* pNewDirectory;
        pNewDirectory=pDirectory->mkdir("cleanedMuons");
        cleanedMuons_.book( pNewDirectory );

        pNewDirectory=pDirectory->mkdir("highestPtDiMuon");
        highestPtDiMuon_.book( pNewDirectory );

        pNewDirectory=pDirectory->mkdir("lowestPtDiMuon");
        lowestPtDiMuon_.book( pNewDirectory );

        pNewDirectory=pDirectory->mkdir("mostIsolatedDiMuon");
        mostIsolatedDiMuon_.book( pNewDirectory );

        pNewDirectory=pDirectory->mkdir("leastIsolatedDiMuon");
        leastIsolatedDiMuon_.book( pNewDirectory );

        pNewDirectory=pDirectory->mkdir("leastDeltaBetaCorrectedIsolatedDiMuon");
        leastDeltaBetaCorrectedIsolatedDiMuon_.book( pNewDirectory );

        histogramHaveBeenBooked_=true;
}

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

        try
        {
                // Clean the muons with all the criteria of the VHbb package (apart from isolation)
                std::vector<VHbbEvent::MuonInfo> cleanedMuons=cleanMuons( event.muInfo );
                // Then make sure they're sorted high to low by pT
                std::sort( cleanedMuons.begin(), cleanedMuons.end(), ::SortByTransverseMomentum<VHbbEvent::MuonInfo>() );
                // Then find two with opposite charge
                std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo> diMuons=findOppositelyChargedPair(cleanedMuons);

                int numberOfPrimaryVertices;
                if( pAuxInfo ) numberOfPrimaryVertices=pAuxInfo->pvInfo.nVertices;
                else numberOfPrimaryVertices=-1;

                cleanedMuons_.fill( cleanedMuons, pAuxInfo ); // Fill some plots about the cleaned muons
                highestPtDiMuon_.fill( diMuons.first, pAuxInfo ); // I know the collection is already sorted by pT, so the first has the highest.
                lowestPtDiMuon_.fill( diMuons.second, pAuxInfo );

                if( firstMuonIsMoreIsolated(diMuons) )
                {
                        mostIsolatedDiMuon_.fill( diMuons.first, pAuxInfo );
                        leastIsolatedDiMuon_.fill( diMuons.second, pAuxInfo );
                }
                else
                {
                        mostIsolatedDiMuon_.fill( diMuons.second, pAuxInfo );
                        leastIsolatedDiMuon_.fill( diMuons.first, pAuxInfo );
                }

                if( firstMuonIsMoreDeltaBetaCorrectedIsolated(diMuons) ) leastDeltaBetaCorrectedIsolatedDiMuon_.fill( diMuons.second, pAuxInfo );
                else leastDeltaBetaCorrectedIsolatedDiMuon_.fill( diMuons.first, pAuxInfo );
        }
        catch( ::AnalysisFailureException& exception )
        {
                // Um, not sure what to do. Nothing and carry on I suppose. I'm expecting a fair few cases where two oppositely charged muons
                // can't be found so I won't print a warning.
        }

//      pGlobalChi2_->Fill( muon.globChi2 );
}

std::vector<VHbbEvent::MuonInfo> trkupgradeanalysis::IsolationStudyPlotSet::cleanMuons( const std::vector<VHbbEvent::MuonInfo>& muons )
{
        std::vector<VHbbEvent::MuonInfo> returnValue;

        for( std::vector<VHbbEvent::MuonInfo>::const_iterator iMuon=muons.begin(); iMuon!=muons.end(); ++iMuon )
        {
                if( iMuon->globChi2<10 &&
                    iMuon->nPixelHits>= 1 &&
                    iMuon->globNHits != 0 &&
                    iMuon->nHits > 10 &&
                    (iMuon->cat & 0x1) && (iMuon->cat & 0x2) &&
                    iMuon->nMatches >=2 &&
                    iMuon->ipDb<.2 &&
                    fabs(iMuon->p4.Eta())<2.4 &&
                    iMuon->p4.Pt()>20 ) returnValue.push_back(*iMuon);
        }

        return returnValue;
}

std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo> trkupgradeanalysis::IsolationStudyPlotSet::findOppositelyChargedPair( const std::vector<VHbbEvent::MuonInfo>& muons )
{
        // No need to check if there are two or more muons in the collection, because if not the loop will not be
        // entered and the other exception will be triggered.

        size_t otherMuonIndex=0; // Use a default of zero as an invalid value
        for( size_t a=1; a<muons.size(); ++a ) // Run over all muons except the first
        {
                if( muons[0].charge*muons[a].charge < 0 ) // Check to see if this muon has opposite charge to the first
                {
                        otherMuonIndex=a;
                        break;
                }
        }

        if( otherMuonIndex==0 ) throw ::AnalysisFailureException("IsolationStudyPlotSet::findOppositelyChargedPair - Unable to find two oppositely charged muons");

        return std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo>( muons[0], muons[otherMuonIndex] );
}

bool trkupgradeanalysis::IsolationStudyPlotSet::firstMuonIsMoreIsolated( std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo> diMuons )
{
        float firstIsolation=MuonInfoPlotSet::combinedRelativeIsolation(diMuons.first);
        float secondIsolation=MuonInfoPlotSet::combinedRelativeIsolation(diMuons.second);
        // More isolated if the sum of energy deposits around it is less.
        return firstIsolation<secondIsolation;
}

bool trkupgradeanalysis::IsolationStudyPlotSet::firstMuonIsMoreDeltaBetaCorrectedIsolated( std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo> diMuons )
{
        float firstIsolation=MuonInfoPlotSet::deltaBetaCorrectedIsolation(diMuons.first);
        float secondIsolation=MuonInfoPlotSet::deltaBetaCorrectedIsolation(diMuons.second);
        // More isolated if the sum of energy deposits around it is less.
        return firstIsolation<secondIsolation;
}
Revision:	1.1
Committed:	Fri Apr 27 13:52:14 2012 UTC (13 years ago) by grimes
Branch:	MAIN
Log Message:	Added code to analyse muon isolation.
#	Content
1	#include "TrkUpgradeAnalysis/VHbb/interface/IsolationStudyPlotSet.h"
2
3	#include <stdexcept>
4
5	#include <TDirectory.h>
6	#include <TH1F.h>
7
8	#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"
9
10	namespace // Use the unnamed namespace for things only used in this file
11	{
12	/** @brief A functor to sort muon collections by pT in the STL sort algorithm.
13	* @author Mark Grimes (mark.grimes@bristol.ac.uk)
14	* @date 11/Apr/2012
15	*/
16	template <typename T>
17	class SortByTransverseMomentum
18	{
19	public:
20	// Return true if the first argument should come before
21	// the second, as per normal STL sort comparitors
22	bool operator()( const T& first, const T& second )
23	{
24	return first.p4.Pt() > second.p4.Pt();
25	}
26	};
27
28	/** @brief Class to pass as exception when something in the analysis breaks down.
29	*
30	* Intended to be thrown when the event/particle/whatever cannot be processed further
31	* but the analysis as a whole should continue running. A std::runtime_error will be
32	* thrown when everything should stop.
33	*/
34	class AnalysisFailureException : public std::runtime_error
35	{
36	public:
37	AnalysisFailureException( const std::string& what ) : std::runtime_error(what) {}
38	};
39	}
40
41	trkupgradeanalysis::IsolationStudyPlotSet::IsolationStudyPlotSet()
42	: histogramHaveBeenBooked_(false), leastIsolatedDiMuon_(true), leastDeltaBetaCorrectedIsolatedDiMuon_(true)
43	{
44	// No operation besides the initialiser list.
45	}
46
47	void trkupgradeanalysis::IsolationStudyPlotSet::book( TDirectory* pDirectory )
48	{
49	if( histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::IsolationStudyPlotSet::book() - histograms have already been booked" );
50
51	//
52	// Note that the root file which TDirectory is part of takes ownership of all
53	// of these objects, so I don't need to (and shouldn't) delete them when I'm
54	// finished.
55	//
56
57
58	// pGlobalChi2_=new TH1F( "globalChi2","Global chi2", 60, 0, 20 );
59	// pGlobalChi2_->SetDirectory(pDirectory);
60
61	// pLeastIsolatedDiMuonTree_=new TTree("leastIsolatedDiMuonTree","Data about the least isolated muon in the di-muon pair");
62	// pLeastIsolatedDiMuonTree_->Branch("numberOfPrimaryVertices",&numberOfPrimaryVertices_branch_,"numberOfPrimaryVertices/I");
63	// pLeastIsolatedDiMuonTree_->Branch("chargedIsolation",&chargedIsolation_branch_,"chargedIsolation/F");
64	// pLeastIsolatedDiMuonTree_->Branch("photonIsolation",&photonIsolation_branch_,"photonIsolation/F");
65	// pLeastIsolatedDiMuonTree_->Branch("neutralIsolation",&neutralIsolation_branch_,"neutralIsolation/F");
66	// pLeastIsolatedDiMuonTree_->Branch("pileupIsolation",&pileupIsolation_branch_,"pileupIsolation/F");
67	// pLeastIsolatedDiMuonTree_->Branch("pT",&pT_branch_,"pT/F");
68	// pLeastIsolatedDiMuonTree_->SetDirectory( pDirectory );
69
70	TDirectory* pNewDirectory;
71	pNewDirectory=pDirectory->mkdir("cleanedMuons");
72	cleanedMuons_.book( pNewDirectory );
73
74	pNewDirectory=pDirectory->mkdir("highestPtDiMuon");
75	highestPtDiMuon_.book( pNewDirectory );
76
77	pNewDirectory=pDirectory->mkdir("lowestPtDiMuon");
78	lowestPtDiMuon_.book( pNewDirectory );
79
80	pNewDirectory=pDirectory->mkdir("mostIsolatedDiMuon");
81	mostIsolatedDiMuon_.book( pNewDirectory );
82
83	pNewDirectory=pDirectory->mkdir("leastIsolatedDiMuon");
84	leastIsolatedDiMuon_.book( pNewDirectory );
85
86	pNewDirectory=pDirectory->mkdir("leastDeltaBetaCorrectedIsolatedDiMuon");
87	leastDeltaBetaCorrectedIsolatedDiMuon_.book( pNewDirectory );
88
89	histogramHaveBeenBooked_=true;
90	}
91
92	void trkupgradeanalysis::IsolationStudyPlotSet::fill( const VHbbEvent& event, const VHbbEventAuxInfo* pAuxInfo )
93	{
94	if( !histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::IsolationStudyPlotSet::book() - histograms have not been booked" );
95
96	try
97	{
98	// Clean the muons with all the criteria of the VHbb package (apart from isolation)
99	std::vector<VHbbEvent::MuonInfo> cleanedMuons=cleanMuons( event.muInfo );
100	// Then make sure they're sorted high to low by pT
101	std::sort( cleanedMuons.begin(), cleanedMuons.end(), ::SortByTransverseMomentum<VHbbEvent::MuonInfo>() );
102	// Then find two with opposite charge
103	std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo> diMuons=findOppositelyChargedPair(cleanedMuons);
104
105	int numberOfPrimaryVertices;
106	if( pAuxInfo ) numberOfPrimaryVertices=pAuxInfo->pvInfo.nVertices;
107	else numberOfPrimaryVertices=-1;
108
109	cleanedMuons_.fill( cleanedMuons, pAuxInfo ); // Fill some plots about the cleaned muons
110	highestPtDiMuon_.fill( diMuons.first, pAuxInfo ); // I know the collection is already sorted by pT, so the first has the highest.
111	lowestPtDiMuon_.fill( diMuons.second, pAuxInfo );
112
113	if( firstMuonIsMoreIsolated(diMuons) )
114	{
115	mostIsolatedDiMuon_.fill( diMuons.first, pAuxInfo );
116	leastIsolatedDiMuon_.fill( diMuons.second, pAuxInfo );
117	}
118	else
119	{
120	mostIsolatedDiMuon_.fill( diMuons.second, pAuxInfo );
121	leastIsolatedDiMuon_.fill( diMuons.first, pAuxInfo );
122	}
123
124	if( firstMuonIsMoreDeltaBetaCorrectedIsolated(diMuons) ) leastDeltaBetaCorrectedIsolatedDiMuon_.fill( diMuons.second, pAuxInfo );
125	else leastDeltaBetaCorrectedIsolatedDiMuon_.fill( diMuons.first, pAuxInfo );
126	}
127	catch( ::AnalysisFailureException& exception )
128	{
129	// Um, not sure what to do. Nothing and carry on I suppose. I'm expecting a fair few cases where two oppositely charged muons
130	// can't be found so I won't print a warning.
131	}
132
133	// pGlobalChi2_->Fill( muon.globChi2 );
134	}
135
136	std::vector<VHbbEvent::MuonInfo> trkupgradeanalysis::IsolationStudyPlotSet::cleanMuons( const std::vector<VHbbEvent::MuonInfo>& muons )
137	{
138	std::vector<VHbbEvent::MuonInfo> returnValue;
139
140	for( std::vector<VHbbEvent::MuonInfo>::const_iterator iMuon=muons.begin(); iMuon!=muons.end(); ++iMuon )
141	{
142	if( iMuon->globChi2<10 &&
143	iMuon->nPixelHits>= 1 &&
144	iMuon->globNHits != 0 &&
145	iMuon->nHits > 10 &&
146	(iMuon->cat & 0x1) && (iMuon->cat & 0x2) &&
147	iMuon->nMatches >=2 &&
148	iMuon->ipDb<.2 &&
149	fabs(iMuon->p4.Eta())<2.4 &&
150	iMuon->p4.Pt()>20 ) returnValue.push_back(*iMuon);
151	}
152
153	return returnValue;
154	}
155
156	std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo> trkupgradeanalysis::IsolationStudyPlotSet::findOppositelyChargedPair( const std::vector<VHbbEvent::MuonInfo>& muons )
157	{
158	// No need to check if there are two or more muons in the collection, because if not the loop will not be
159	// entered and the other exception will be triggered.
160
161	size_t otherMuonIndex=0; // Use a default of zero as an invalid value
162	for( size_t a=1; a<muons.size(); ++a ) // Run over all muons except the first
163	{
164	if( muons[0].charge*muons[a].charge < 0 ) // Check to see if this muon has opposite charge to the first
165	{
166	otherMuonIndex=a;
167	break;
168	}
169	}
170
171	if( otherMuonIndex==0 ) throw ::AnalysisFailureException("IsolationStudyPlotSet::findOppositelyChargedPair - Unable to find two oppositely charged muons");
172
173	return std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo>( muons[0], muons[otherMuonIndex] );
174	}
175
176	bool trkupgradeanalysis::IsolationStudyPlotSet::firstMuonIsMoreIsolated( std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo> diMuons )
177	{
178	float firstIsolation=MuonInfoPlotSet::combinedRelativeIsolation(diMuons.first);
179	float secondIsolation=MuonInfoPlotSet::combinedRelativeIsolation(diMuons.second);
180	// More isolated if the sum of energy deposits around it is less.
181	return firstIsolation<secondIsolation;
182	}
183
184	bool trkupgradeanalysis::IsolationStudyPlotSet::firstMuonIsMoreDeltaBetaCorrectedIsolated( std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo> diMuons )
185	{
186	float firstIsolation=MuonInfoPlotSet::deltaBetaCorrectedIsolation(diMuons.first);
187	float secondIsolation=MuonInfoPlotSet::deltaBetaCorrectedIsolation(diMuons.second);
188	// More isolated if the sum of energy deposits around it is less.
189	return firstIsolation<secondIsolation;
190	}