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
#	Content
1	#include "TrkUpgradeAnalysis/VHbb/interface/tools.h"
2
3	#include <stdexcept>
4	#include <iostream>
5
6	#include <TFile.h>
7	#include <TH1F.h>
8	#include <TTree.h>
9	#include <TDirectory.h>
10
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	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
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
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	}