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
#	User	Rev	Content
1	grimes	1.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	grimes	1.2	: histogramHaveBeenBooked_(false), leastIsolatedDiMuon_(false), leastDeltaBetaCorrectedIsolatedDiMuon_(false)
43	grimes	1.1	{
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	grimes	1.2	pNewDirectory=pDirectory->mkdir("cleanedElectrons");
63			cleanedElectrons_.book( pNewDirectory );
64
65	grimes	1.1	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	grimes	1.2	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	grimes	1.1	// 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	grimes	1.2	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	grimes	1.1	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			}