VHbbAnalysisCode/src/tools.cpp

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

#include <stdexcept>
#include <iostream>

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

trkupgradeanalysis::tools::TFile_auto_ptr::TFile_auto_ptr( const std::string& filename )
{
        this->reset( TFile::Open( filename.c_str() ) );
}

trkupgradeanalysis::tools::TFile_auto_ptr::~TFile_auto_ptr()
{
        // If the pointer is not null, close the file. Let the auto_ptr destructor delete it.
        if( this->get() ) this->get()->Close();
}

bool trkupgradeanalysis::tools::TFile_auto_ptr::operator==( const TFile* pOtherFile )
{
        return this->get()==pOtherFile;
}


trkupgradeanalysis::tools::NTupleRow::NTupleRow( TTree* pTree )
        : pTree_(pTree)
{
        TObjArray* pLeafList=pTree->GetListOfBranches();

        for( int a=0; a<pLeafList->GetEntries(); ++a )
        {
                TBranch* pBranch=(TBranch*)(*pLeafList)[a];
                std::string name=pBranch->GetName();
                std::string title=pBranch->GetTitle();
                std::string typeString=title.substr( name.size() );

                if( typeString=="/D" ) branchAddresses_[name]=0; // Put an entry in the map for this name
        }

        // Now loop over the map names and set the branch addresses. I didn't want to do this earlier
        // because there's the possibility adding new entries to the map might change the addresses.
        for( std::map<std::string,double>::iterator iEntry=branchAddresses_.begin(); iEntry!=branchAddresses_.end(); ++iEntry )
        {
                pTree->SetBranchAddress( iEntry->first.c_str(), &(iEntry->second) );
        }

        maxCandidateNumber_=pTree_->GetEntries();
        candidateNumber_=-1; // Set to minus one, so that the first call to nextCandidate() sets it to the first (0) entry
        pTree_->GetEntry(0); // Load in the first candidate so that the memory contents isn't random. Note the first call to nextCandidate() will keep it on the 0th entry.
}

const double& trkupgradeanalysis::tools::NTupleRow::getDouble( const std::string& name ) const
{
        std::map<std::string,double>::const_iterator iFindResult=branchAddresses_.find(name);
        if( iFindResult!=branchAddresses_.end() ) return iFindResult->second;
        else throw std::runtime_error( "trkupgradeanalysis::tools::NTupleRow::getDouble(..) - TTree doesn't have a branch called "+name );
}

bool trkupgradeanalysis::tools::NTupleRow::nextRow()
{
        ++candidateNumber_;
        if( candidateNumber_>=maxCandidateNumber_ ) return false;

        pTree_->GetEntry(candidateNumber_);
        return true;
}

void trkupgradeanalysis::tools::NTupleRow::returnToStart()
{
        candidateNumber_=-1; // Set to minus one, so that the first call to nextCandidate() sets it to the first (0) entry
}


float trkupgradeanalysis::tools::linearInterpolate( std::pair<float,float> point1, std::pair<float,float> point2, float requiredY )
{
        // Using the equation for a straight line of y=mx+c, work out m and c so
        // that the line passes through point1 and point2
        float m=(point1.second-point2.second)/(point1.first-point2.first);
        float c=point2.second-m*point2.first;

        // Then I can return x=(y-c)/m
        return (requiredY-c)/m;
}


float trkupgradeanalysis::tools::findDiscriminator( const TH1F* pEventsVersusDiscriminatorHistogram, const float requiredEfficiency, const bool verbose )
{
        if( pEventsVersusDiscriminatorHistogram==NULL ) throw std::runtime_error("findDiscriminator was called with a NULL histogram pointer");

        TAxis* pHistogramXAxis=pEventsVersusDiscriminatorHistogram->GetXaxis();
        float accumulatedEvents=0; // The number of events that would pass the current cut (i.e. number of events in this bin and all higher ones)
        float previousEfficiency=1; // The efficiency on the last run through the loop. Set high for logic to work

        // Loop backwards over each of the bins in the histogram
        for( int bin=pHistogramXAxis->GetNbins(); bin>=1; --bin )
        {
                accumulatedEvents+=pEventsVersusDiscriminatorHistogram->GetBinContent(bin);

                float currentEfficiency=accumulatedEvents/pEventsVersusDiscriminatorHistogram->GetEntries();
                // See if this loop has gone a bin too far. If so interpolate between this and the previous loop.
                if( currentEfficiency>requiredEfficiency )
                {
                        // The working point is somewhere between the previous bin and this bin.
                        // Do a linear interpolation between the bins.
                        std::pair<float,float> point1( pHistogramXAxis->GetBinLowEdge(bin), currentEfficiency );
                        std::pair<float,float> point2( pHistogramXAxis->GetBinLowEdge(bin+1), previousEfficiency );
                        float operatingPoint=trkupgradeanalysis::tools::linearInterpolate( point1, point2, requiredEfficiency );

                        if( verbose )
                        {
                                std::cout << "[findDiscriminator verbose] The operating point of " << requiredEfficiency << " efficiency for \"" << pEventsVersusDiscriminatorHistogram->GetName() << "\" is between:" << "\n"
                                                << "[findDiscriminator verbose] \t" << pHistogramXAxis->GetBinLowEdge(bin) << " -> " << currentEfficiency << " (bin " << bin << ") and " << "\n"
                                                << "[findDiscriminator verbose] \t" << pHistogramXAxis->GetBinLowEdge(bin+1) << " -> " << previousEfficiency << " (bin " << bin+1 << ")." << "\n"
                                                << "[findDiscriminator verbose] Linear interpolation gives:" << "\n"
                                                << "[findDiscriminator verbose] \t" << operatingPoint << std::endl;
                        }

                        return operatingPoint;
                }
                previousEfficiency=currentEfficiency; // Set this ready for the next loop around

        } // end of backwards loop over histogram bins

        // If flow has got this far, something has gone wrong
        throw std::runtime_error("findDiscriminator was unable to find an operating point");

} // end of function findDiscriminator


float trkupgradeanalysis::tools::findEfficiency( const TH1F* pEventsVersusDiscriminatorHistogram, const float requiredDiscriminator, const bool verbose )
{
        if( pEventsVersusDiscriminatorHistogram==NULL ) throw std::runtime_error("findEfficiency was called with a NULL histogram pointer");

        TAxis* pHistogramXAxis=pEventsVersusDiscriminatorHistogram->GetXaxis();
        float accumulatedEvents=0; // The number of events that would pass the current cut (i.e. number of events in this bin and all higher ones)
        float previousEfficiency=1; // The efficiency on the last run through the loop. Set high for logic to work

        // Loop backwards over each of the bins in the histogram
        for( int bin=pHistogramXAxis->GetNbins(); bin>=1; --bin )
        {
                accumulatedEvents+=pEventsVersusDiscriminatorHistogram->GetBinContent(bin);

                float currentEfficiency=accumulatedEvents/pEventsVersusDiscriminatorHistogram->GetEntries();
                // See if this loop step has gone too far, if so then interpolate between this and the previous step.
                if( pHistogramXAxis->GetBinLowEdge(bin)<requiredDiscriminator )
                {
                        // The efficiency is somewhere between the previous bin and this bin.
                        // Do a linear interpolation between the bins. This is similar to findDiscriminator but
                        // I've switched around the x and y coordinates.
                        std::pair<float,float> point1( currentEfficiency, pHistogramXAxis->GetBinLowEdge(bin) );
                        std::pair<float,float> point2( previousEfficiency, pHistogramXAxis->GetBinLowEdge(bin+1) );
                        float efficiency=trkupgradeanalysis::tools::linearInterpolate( point1, point2, requiredDiscriminator );

                        if( verbose )
                        {
                                std::cout << "[findEfficiency verbose] The discriminator of " << requiredDiscriminator << " for \"" << pEventsVersusDiscriminatorHistogram->GetName() << "\" is between:" << "\n"
                                                << "[findEfficiency verbose] \t" << pHistogramXAxis->GetBinLowEdge(bin) << " -> " << currentEfficiency << " (bin " << bin << ") and " << "\n"
                                                << "[findEfficiency verbose] \t" << pHistogramXAxis->GetBinLowEdge(bin+1) << " -> " << previousEfficiency << " (bin " << bin+1 << ")." << "\n"
                                                << "[findEfficiency verbose] Linear interpolation gives an efficiency of:" << "\n"
                                                << "[findEfficiency verbose] \t" << efficiency << std::endl;
                        }

                        return efficiency;
                }
                previousEfficiency=currentEfficiency; // Set this ready for the next loop around

        } // end of backwards loop over histogram bins

        // If flow has got this far, something has gone wrong
        throw std::runtime_error("findEfficiency was unable to find an operating point");

} // end of function findEfficiency


TDirectory* trkupgradeanalysis::tools::createDirectory( const std::string& fullPath, TDirectory* pParent )
{
        if( pParent==NULL ) throw std::runtime_error( "The parent directory is a Null pointer" );

        TDirectory* pSubDirectory=pParent;
        size_t currentPosition=0;
        size_t nextSlash;
        do
        {
                nextSlash=fullPath.find_first_of('/', currentPosition );
                std::string directoryName=fullPath.substr(currentPosition,nextSlash-currentPosition);
                currentPosition=nextSlash+1;

                TDirectory* pNextSubDirectory=pSubDirectory->GetDirectory( directoryName.c_str() );
                if( pNextSubDirectory==NULL ) pNextSubDirectory=pSubDirectory->mkdir( directoryName.c_str() );
                if( pNextSubDirectory==NULL ) throw std::runtime_error( "Couldn't create the root directory \""+directoryName+"\"" );
                pSubDirectory=pNextSubDirectory;

        } while( nextSlash!=std::string::npos );

        return pSubDirectory;
}
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:	+79 -0 lines
Log Message:	Long overdue commit with several new files
#	User	Rev	Content
1	grimes	1.1	#include "TrkUpgradeAnalysis/VHbb/interface/tools.h"
2
3			#include <stdexcept>
4			#include <iostream>
5
6			#include <TFile.h>
7			#include <TH1F.h>
8	grimes	1.2	#include <TTree.h>
9			#include <TDirectory.h>
10	grimes	1.1
11			trkupgradeanalysis::tools::TFile_auto_ptr::TFile_auto_ptr( const std::string& filename )
12			{
13			this->reset( TFile::Open( filename.c_str() ) );
14			}
15
16			trkupgradeanalysis::tools::TFile_auto_ptr::~TFile_auto_ptr()
17			{
18			// If the pointer is not null, close the file. Let the auto_ptr destructor delete it.
19			if( this->get() ) this->get()->Close();
20			}
21
22			bool trkupgradeanalysis::tools::TFile_auto_ptr::operator==( const TFile* pOtherFile )
23			{
24			return this->get()==pOtherFile;
25			}
26
27
28
29
30
31	grimes	1.2	trkupgradeanalysis::tools::NTupleRow::NTupleRow( TTree* pTree )
32			: pTree_(pTree)
33			{
34			TObjArray* pLeafList=pTree->GetListOfBranches();
35
36			for( int a=0; a<pLeafList->GetEntries(); ++a )
37			{
38			TBranch* pBranch=(TBranch)(pLeafList)[a];
39			std::string name=pBranch->GetName();
40			std::string title=pBranch->GetTitle();
41			std::string typeString=title.substr( name.size() );
42
43			if( typeString=="/D" ) branchAddresses_[name]=0; // Put an entry in the map for this name
44			}
45
46			// Now loop over the map names and set the branch addresses. I didn't want to do this earlier
47			// because there's the possibility adding new entries to the map might change the addresses.
48			for( std::map<std::string,double>::iterator iEntry=branchAddresses_.begin(); iEntry!=branchAddresses_.end(); ++iEntry )
49			{
50			pTree->SetBranchAddress( iEntry->first.c_str(), &(iEntry->second) );
51			}
52
53			maxCandidateNumber_=pTree_->GetEntries();
54			candidateNumber_=-1; // Set to minus one, so that the first call to nextCandidate() sets it to the first (0) entry
55			pTree_->GetEntry(0); // Load in the first candidate so that the memory contents isn't random. Note the first call to nextCandidate() will keep it on the 0th entry.
56			}
57
58			const double& trkupgradeanalysis::tools::NTupleRow::getDouble( const std::string& name ) const
59			{
60			std::map<std::string,double>::const_iterator iFindResult=branchAddresses_.find(name);
61			if( iFindResult!=branchAddresses_.end() ) return iFindResult->second;
62			else throw std::runtime_error( "trkupgradeanalysis::tools::NTupleRow::getDouble(..) - TTree doesn't have a branch called "+name );
63			}
64
65			bool trkupgradeanalysis::tools::NTupleRow::nextRow()
66			{
67			++candidateNumber_;
68			if( candidateNumber_>=maxCandidateNumber_ ) return false;
69
70			pTree_->GetEntry(candidateNumber_);
71			return true;
72			}
73
74			void trkupgradeanalysis::tools::NTupleRow::returnToStart()
75			{
76			candidateNumber_=-1; // Set to minus one, so that the first call to nextCandidate() sets it to the first (0) entry
77			}
78
79
80
81
82
83	grimes	1.1
84			float trkupgradeanalysis::tools::linearInterpolate( std::pair<float,float> point1, std::pair<float,float> point2, float requiredY )
85			{
86			// Using the equation for a straight line of y=mx+c, work out m and c so
87			// that the line passes through point1 and point2
88			float m=(point1.second-point2.second)/(point1.first-point2.first);
89			float c=point2.second-m*point2.first;
90
91			// Then I can return x=(y-c)/m
92			return (requiredY-c)/m;
93			}
94
95
96
97
98
99
100			float trkupgradeanalysis::tools::findDiscriminator( const TH1F* pEventsVersusDiscriminatorHistogram, const float requiredEfficiency, const bool verbose )
101			{
102			if( pEventsVersusDiscriminatorHistogram==NULL ) throw std::runtime_error("findDiscriminator was called with a NULL histogram pointer");
103
104			TAxis* pHistogramXAxis=pEventsVersusDiscriminatorHistogram->GetXaxis();
105			float accumulatedEvents=0; // The number of events that would pass the current cut (i.e. number of events in this bin and all higher ones)
106			float previousEfficiency=1; // The efficiency on the last run through the loop. Set high for logic to work
107
108			// Loop backwards over each of the bins in the histogram
109			for( int bin=pHistogramXAxis->GetNbins(); bin>=1; --bin )
110			{
111			accumulatedEvents+=pEventsVersusDiscriminatorHistogram->GetBinContent(bin);
112
113			float currentEfficiency=accumulatedEvents/pEventsVersusDiscriminatorHistogram->GetEntries();
114			// See if this loop has gone a bin too far. If so interpolate between this and the previous loop.
115			if( currentEfficiency>requiredEfficiency )
116			{
117			// The working point is somewhere between the previous bin and this bin.
118			// Do a linear interpolation between the bins.
119			std::pair<float,float> point1( pHistogramXAxis->GetBinLowEdge(bin), currentEfficiency );
120			std::pair<float,float> point2( pHistogramXAxis->GetBinLowEdge(bin+1), previousEfficiency );
121			float operatingPoint=trkupgradeanalysis::tools::linearInterpolate( point1, point2, requiredEfficiency );
122
123			if( verbose )
124			{
125			std::cout << "[findDiscriminator verbose] The operating point of " << requiredEfficiency << " efficiency for \"" << pEventsVersusDiscriminatorHistogram->GetName() << "\" is between:" << "\n"
126			<< "[findDiscriminator verbose] \t" << pHistogramXAxis->GetBinLowEdge(bin) << " -> " << currentEfficiency << " (bin " << bin << ") and " << "\n"
127			<< "[findDiscriminator verbose] \t" << pHistogramXAxis->GetBinLowEdge(bin+1) << " -> " << previousEfficiency << " (bin " << bin+1 << ")." << "\n"
128			<< "[findDiscriminator verbose] Linear interpolation gives:" << "\n"
129			<< "[findDiscriminator verbose] \t" << operatingPoint << std::endl;
130			}
131
132			return operatingPoint;
133			}
134			previousEfficiency=currentEfficiency; // Set this ready for the next loop around
135
136			} // end of backwards loop over histogram bins
137
138			// If flow has got this far, something has gone wrong
139			throw std::runtime_error("findDiscriminator was unable to find an operating point");
140
141			} // end of function findDiscriminator
142
143
144
145
146
147
148
149			float trkupgradeanalysis::tools::findEfficiency( const TH1F* pEventsVersusDiscriminatorHistogram, const float requiredDiscriminator, const bool verbose )
150			{
151			if( pEventsVersusDiscriminatorHistogram==NULL ) throw std::runtime_error("findEfficiency was called with a NULL histogram pointer");
152
153			TAxis* pHistogramXAxis=pEventsVersusDiscriminatorHistogram->GetXaxis();
154			float accumulatedEvents=0; // The number of events that would pass the current cut (i.e. number of events in this bin and all higher ones)
155			float previousEfficiency=1; // The efficiency on the last run through the loop. Set high for logic to work
156
157			// Loop backwards over each of the bins in the histogram
158			for( int bin=pHistogramXAxis->GetNbins(); bin>=1; --bin )
159			{
160			accumulatedEvents+=pEventsVersusDiscriminatorHistogram->GetBinContent(bin);
161
162			float currentEfficiency=accumulatedEvents/pEventsVersusDiscriminatorHistogram->GetEntries();
163			// See if this loop step has gone too far, if so then interpolate between this and the previous step.
164			if( pHistogramXAxis->GetBinLowEdge(bin)<requiredDiscriminator )
165			{
166			// The efficiency is somewhere between the previous bin and this bin.
167			// Do a linear interpolation between the bins. This is similar to findDiscriminator but
168			// I've switched around the x and y coordinates.
169			std::pair<float,float> point1( currentEfficiency, pHistogramXAxis->GetBinLowEdge(bin) );
170			std::pair<float,float> point2( previousEfficiency, pHistogramXAxis->GetBinLowEdge(bin+1) );
171			float efficiency=trkupgradeanalysis::tools::linearInterpolate( point1, point2, requiredDiscriminator );
172
173			if( verbose )
174			{
175			std::cout << "[findEfficiency verbose] The discriminator of " << requiredDiscriminator << " for \"" << pEventsVersusDiscriminatorHistogram->GetName() << "\" is between:" << "\n"
176			<< "[findEfficiency verbose] \t" << pHistogramXAxis->GetBinLowEdge(bin) << " -> " << currentEfficiency << " (bin " << bin << ") and " << "\n"
177			<< "[findEfficiency verbose] \t" << pHistogramXAxis->GetBinLowEdge(bin+1) << " -> " << previousEfficiency << " (bin " << bin+1 << ")." << "\n"
178			<< "[findEfficiency verbose] Linear interpolation gives an efficiency of:" << "\n"
179			<< "[findEfficiency verbose] \t" << efficiency << std::endl;
180			}
181
182			return efficiency;
183			}
184			previousEfficiency=currentEfficiency; // Set this ready for the next loop around
185
186			} // end of backwards loop over histogram bins
187
188			// If flow has got this far, something has gone wrong
189			throw std::runtime_error("findEfficiency was unable to find an operating point");
190
191			} // end of function findEfficiency
192	grimes	1.2
193
194
195			TDirectory* trkupgradeanalysis::tools::createDirectory( const std::string& fullPath, TDirectory* pParent )
196			{
197			if( pParent==NULL ) throw std::runtime_error( "The parent directory is a Null pointer" );
198
199			TDirectory* pSubDirectory=pParent;
200			size_t currentPosition=0;
201			size_t nextSlash;
202			do
203			{
204			nextSlash=fullPath.find_first_of('/', currentPosition );
205			std::string directoryName=fullPath.substr(currentPosition,nextSlash-currentPosition);
206			currentPosition=nextSlash+1;
207
208			TDirectory* pNextSubDirectory=pSubDirectory->GetDirectory( directoryName.c_str() );
209			if( pNextSubDirectory==NULL ) pNextSubDirectory=pSubDirectory->mkdir( directoryName.c_str() );
210			if( pNextSubDirectory==NULL ) throw std::runtime_error( "Couldn't create the root directory \""+directoryName+"\"" );
211			pSubDirectory=pNextSubDirectory;
212
213			} while( nextSlash!=std::string::npos );
214
215			return pSubDirectory;
216			}