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_(false), leastDeltaBetaCorrectedIsolatedDiMuon_(false)
{
        // 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.
        //


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

        pNewDirectory=pDirectory->mkdir("cleanedElectrons");
        cleanedElectrons_.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 );
                cleanedMuons_.fill( cleanedMuons, pAuxInfo ); // Fill some plots about the cleaned muons

                // Clean the muons with all the criteria of the VHbb package (apart from isolation)
                std::vector<VHbbEvent::ElectronInfo> cleanedElectrons=cleanElectrons( event.eleInfo );
                cleanedElectrons_.fill( cleanedElectrons, pAuxInfo ); // Fill some plots about the cleaned muons

                // 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;

                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::vector<VHbbEvent::ElectronInfo> trkupgradeanalysis::IsolationStudyPlotSet::cleanElectrons( const std::vector<VHbbEvent::ElectronInfo>& electrons )
{
        std::vector<VHbbEvent::ElectronInfo> returnValue;

        for( std::vector<VHbbEvent::ElectronInfo>::const_iterator iElectron=electrons.begin(); iElectron!=electrons.end(); ++iElectron )
        {
                if( (std::fabs(iElectron->id95-7)<0.1 || std::fabs(iElectron->id95-5)<0.1) &&
                        std::fabs(iElectron->p4.Eta()) < 2.5 &&
                        iElectron->p4.Pt()>20 ) returnValue.push_back(*iElectron);
        }

        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.2
Committed:	Wed Aug 15 22:37:47 2012 UTC (12 years, 8 months ago) by grimes
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.1:	+24 -14 lines
Log Message:	Long overdue commit with several new files
#	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_(false), leastDeltaBetaCorrectedIsolatedDiMuon_(false)
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	TDirectory* pNewDirectory;
59	pNewDirectory=pDirectory->mkdir("cleanedMuons");
60	cleanedMuons_.book( pNewDirectory );
61
62	pNewDirectory=pDirectory->mkdir("cleanedElectrons");
63	cleanedElectrons_.book( pNewDirectory );
64
65	pNewDirectory=pDirectory->mkdir("highestPtDiMuon");
66	highestPtDiMuon_.book( pNewDirectory );
67
68	pNewDirectory=pDirectory->mkdir("lowestPtDiMuon");
69	lowestPtDiMuon_.book( pNewDirectory );
70
71	pNewDirectory=pDirectory->mkdir("mostIsolatedDiMuon");
72	mostIsolatedDiMuon_.book( pNewDirectory );
73
74	pNewDirectory=pDirectory->mkdir("leastIsolatedDiMuon");
75	leastIsolatedDiMuon_.book( pNewDirectory );
76
77	pNewDirectory=pDirectory->mkdir("leastDeltaBetaCorrectedIsolatedDiMuon");
78	leastDeltaBetaCorrectedIsolatedDiMuon_.book( pNewDirectory );
79
80	histogramHaveBeenBooked_=true;
81	}
82
83	void trkupgradeanalysis::IsolationStudyPlotSet::fill( const VHbbEvent& event, const VHbbEventAuxInfo* pAuxInfo )
84	{
85	if( !histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::IsolationStudyPlotSet::book() - histograms have not been booked" );
86
87	try
88	{
89	// Clean the muons with all the criteria of the VHbb package (apart from isolation)
90	std::vector<VHbbEvent::MuonInfo> cleanedMuons=cleanMuons( event.muInfo );
91	cleanedMuons_.fill( cleanedMuons, pAuxInfo ); // Fill some plots about the cleaned muons
92
93	// Clean the muons with all the criteria of the VHbb package (apart from isolation)
94	std::vector<VHbbEvent::ElectronInfo> cleanedElectrons=cleanElectrons( event.eleInfo );
95	cleanedElectrons_.fill( cleanedElectrons, pAuxInfo ); // Fill some plots about the cleaned muons
96
97	// Then make sure they're sorted high to low by pT
98	std::sort( cleanedMuons.begin(), cleanedMuons.end(), ::SortByTransverseMomentum<VHbbEvent::MuonInfo>() );
99	// Then find two with opposite charge
100	std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo> diMuons=findOppositelyChargedPair(cleanedMuons);
101
102	int numberOfPrimaryVertices;
103	if( pAuxInfo ) numberOfPrimaryVertices=pAuxInfo->pvInfo.nVertices;
104	else numberOfPrimaryVertices=-1;
105
106	highestPtDiMuon_.fill( diMuons.first, pAuxInfo ); // I know the collection is already sorted by pT, so the first has the highest.
107	lowestPtDiMuon_.fill( diMuons.second, pAuxInfo );
108
109	if( firstMuonIsMoreIsolated(diMuons) )
110	{
111	mostIsolatedDiMuon_.fill( diMuons.first, pAuxInfo );
112	leastIsolatedDiMuon_.fill( diMuons.second, pAuxInfo );
113	}
114	else
115	{
116	mostIsolatedDiMuon_.fill( diMuons.second, pAuxInfo );
117	leastIsolatedDiMuon_.fill( diMuons.first, pAuxInfo );
118	}
119
120	if( firstMuonIsMoreDeltaBetaCorrectedIsolated(diMuons) ) leastDeltaBetaCorrectedIsolatedDiMuon_.fill( diMuons.second, pAuxInfo );
121	else leastDeltaBetaCorrectedIsolatedDiMuon_.fill( diMuons.first, pAuxInfo );
122	}
123	catch( ::AnalysisFailureException& exception )
124	{
125	// Um, not sure what to do. Nothing and carry on I suppose. I'm expecting a fair few cases where two oppositely charged muons
126	// can't be found so I won't print a warning.
127	}
128
129	// pGlobalChi2_->Fill( muon.globChi2 );
130	}
131
132	std::vector<VHbbEvent::MuonInfo> trkupgradeanalysis::IsolationStudyPlotSet::cleanMuons( const std::vector<VHbbEvent::MuonInfo>& muons )
133	{
134	std::vector<VHbbEvent::MuonInfo> returnValue;
135
136	for( std::vector<VHbbEvent::MuonInfo>::const_iterator iMuon=muons.begin(); iMuon!=muons.end(); ++iMuon )
137	{
138	if( iMuon->globChi2<10 &&
139	iMuon->nPixelHits>= 1 &&
140	iMuon->globNHits != 0 &&
141	iMuon->nHits > 10 &&
142	(iMuon->cat & 0x1) && (iMuon->cat & 0x2) &&
143	iMuon->nMatches >=2 &&
144	iMuon->ipDb<.2 &&
145	fabs(iMuon->p4.Eta())<2.4 &&
146	iMuon->p4.Pt()>20 ) returnValue.push_back(*iMuon);
147	}
148
149	return returnValue;
150	}
151
152	std::vector<VHbbEvent::ElectronInfo> trkupgradeanalysis::IsolationStudyPlotSet::cleanElectrons( const std::vector<VHbbEvent::ElectronInfo>& electrons )
153	{
154	std::vector<VHbbEvent::ElectronInfo> returnValue;
155
156	for( std::vector<VHbbEvent::ElectronInfo>::const_iterator iElectron=electrons.begin(); iElectron!=electrons.end(); ++iElectron )
157	{
158	if( (std::fabs(iElectron->id95-7)<0.1 \|\| std::fabs(iElectron->id95-5)<0.1) &&
159	std::fabs(iElectron->p4.Eta()) < 2.5 &&
160	iElectron->p4.Pt()>20 ) returnValue.push_back(*iElectron);
161	}
162
163	return returnValue;
164	}
165
166	std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo> trkupgradeanalysis::IsolationStudyPlotSet::findOppositelyChargedPair( const std::vector<VHbbEvent::MuonInfo>& muons )
167	{
168	// No need to check if there are two or more muons in the collection, because if not the loop will not be
169	// entered and the other exception will be triggered.
170
171	size_t otherMuonIndex=0; // Use a default of zero as an invalid value
172	for( size_t a=1; a<muons.size(); ++a ) // Run over all muons except the first
173	{
174	if( muons[0].charge*muons[a].charge < 0 ) // Check to see if this muon has opposite charge to the first
175	{
176	otherMuonIndex=a;
177	break;
178	}
179	}
180
181	if( otherMuonIndex==0 ) throw ::AnalysisFailureException("IsolationStudyPlotSet::findOppositelyChargedPair - Unable to find two oppositely charged muons");
182
183	return std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo>( muons[0], muons[otherMuonIndex] );
184	}
185
186	bool trkupgradeanalysis::IsolationStudyPlotSet::firstMuonIsMoreIsolated( std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo> diMuons )
187	{
188	float firstIsolation=MuonInfoPlotSet::combinedRelativeIsolation(diMuons.first);
189	float secondIsolation=MuonInfoPlotSet::combinedRelativeIsolation(diMuons.second);
190	// More isolated if the sum of energy deposits around it is less.
191	return firstIsolation<secondIsolation;
192	}
193
194	bool trkupgradeanalysis::IsolationStudyPlotSet::firstMuonIsMoreDeltaBetaCorrectedIsolated( std::pair<VHbbEvent::MuonInfo,VHbbEvent::MuonInfo> diMuons )
195	{
196	float firstIsolation=MuonInfoPlotSet::deltaBetaCorrectedIsolation(diMuons.first);
197	float secondIsolation=MuonInfoPlotSet::deltaBetaCorrectedIsolation(diMuons.second);
198	// More isolated if the sum of energy deposits around it is less.
199	return firstIsolation<secondIsolation;
200	}