ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/grimes/L1Menu/src/TriggerRatePlot.cpp

Comparing UserCode/grimes/L1Menu/src/TriggerRatePlot.cpp (file contents):
Revision 1.9 by grimes, Fri Jul 5 13:54:44 2013 UTC vs.
Revision 1.10 by grimes, Wed Jul 24 11:48:55 2013 UTC

#	Line 6 | Line 6
6		#include "l1menu/IEvent.h"
7		#include "l1menu/ISample.h"
8		#include "l1menu/TriggerTable.h"
9	+	#include "l1menu/tools/tools.h"
10	+	#include "l1menu/tools/stringManipulation.h"
11		#include <TH1F.h>
12		#include <sstream>
13		#include <algorithm>
14	<
14	>	#include <stdexcept>
15		#include <iostream>
16
17		l1menu::TriggerRatePlot::TriggerRatePlot( const l1menu::ITrigger& trigger, std::unique_ptr<TH1> pHistogram, const std::string& versusParameter, const std::vector<std::string> scaledParameters )
18		: pHistogram_( std::move(pHistogram) ), versusParameter_(versusParameter), histogramOwnedByMe_(true)
19		{
20	+	initiate( trigger, scaledParameters );
21	+	}
22	+
23	+	l1menu::TriggerRatePlot::TriggerRatePlot( const l1menu::ITrigger& trigger, const std::string& name, size_t numberOfBins, float lowEdge, float highEdge, const std::string& versusParameter, const std::vector<std::string> scaledParameters )
24	+	: pHistogram_( new TH1F( name.c_str(), "This title gets changed later", numberOfBins, lowEdge, highEdge ) ), versusParameter_(versusParameter), histogramOwnedByMe_(true)
25	+	{
26	+	pHistogram_->SetDirectory( NULL ); // Hold in memory only
27	+	initiate( trigger, scaledParameters );
28	+	}
29	+
30	+	l1menu::TriggerRatePlot::TriggerRatePlot( const TH1* pPreExisitingHistogram )
31	+	{
32	+	// All of the information about the trigger is stored in the title, so examine
33	+	// that to see what the trigger is, the version, and all of the parameters.
34	+	// Format is "<name> rate versus <thresholdName> [v<version>,<param1>=<value>,<param2>=<value>,...]"
35	+	// so first split by whitespace.
36	+	std::vector<std::string> titleSplitByWhitespace=l1menu::tools::splitByWhitespace( pPreExisitingHistogram->GetTitle() );
37	+	std::cout << "titleSplitByWhitespace="; for( const auto& element : titleSplitByWhitespace ) std::cout << "\"" << element << "\","; std::cout << std::endl;
38	+
39	+	if( titleSplitByWhitespace.size()!=5 ) throw std::runtime_error( "TriggerRatePlot cannot be loaded from the supplied histogram because the title is not as expected" );
40	+	std::string triggerName=titleSplitByWhitespace[0];
41	+	versusParameter_=titleSplitByWhitespace[3];
42	+
43	+	// Check that the final element begins with "[" and ends with "]"
44	+	if( titleSplitByWhitespace[4].front()!='[' || titleSplitByWhitespace[4].back()!=']' ) throw std::runtime_error( "TriggerRatePlot cannot be loaded from the supplied histogram because the title is not as expected" );
45	+	// If so, chop off those two parts
46	+	std::string parameterList=titleSplitByWhitespace[4].substr( 1, titleSplitByWhitespace[4].size()-2 );
47	+	// and split by commas
48	+	std::vector<std::string> parameterListSplitByCommas=l1menu::tools::splitByDelimeters( parameterList, "," );
49	+	std::cout << "parameterListSplitByCommas="; for( const auto& element : parameterListSplitByCommas ) std::cout << "\"" << element << "\","; std::cout << std::endl;
50	+	// Need to at least have the version listed
51	+	if( parameterListSplitByCommas.empty() ) throw std::runtime_error( "TriggerRatePlot cannot be loaded from the supplied histogram because the title is not as expected" );
52	+
53	+	// Figure out the version number of the trigger. It should be prefixed with a "v".
54	+	if( parameterListSplitByCommas[0].front()!='v' ) throw std::runtime_error( "TriggerRatePlot cannot be loaded from the supplied histogram because the title is not as expected" );
55	+	std::string versionAsString=parameterListSplitByCommas[0].substr( 1 );
56	+	std::cout << "versionAsString=" << versionAsString << std::endl;
57	+	unsigned int versionNumber=l1menu::tools::convertStringToFloat( versionAsString );
58	+
59	+	// Now I have the trigger name and version number, I can create a copy from the TriggerTable.
60	+	std::unique_ptr<l1menu::ITrigger> pTemporaryTriggerCopy=l1menu::TriggerTable::instance().getTrigger( triggerName, versionNumber );
61	+	if( pTemporaryTriggerCopy.get()==NULL ) throw std::runtime_error( "TriggerRatePlot cannot be loaded from the supplied histogram because trigger is not in the TriggerTable" );
62	+
63	+	// Set the main threshold to 1 just as an arbitrary starting point. I might need to scale thresholds
64	+	// off that in a minute.
65	+	pTemporaryTriggerCopy->parameter( versusParameter_ )=1;
66	+
67	+	std::vector<std::string> scaledParameters; // This will build into a list of anything that should scale with the main threshold
68	+
69	+	// Then loop through the other parameters and set them to whatever they're listed as.
70	+	// Start at 1 because I've already done the version.
71	+	for( size_t index=1; index<parameterListSplitByCommas.size(); ++index )
72	+	{
73	+	std::cout << "Checking string " << parameterListSplitByCommas[index] << std::endl;
74	+	// Split by equals sign
75	+	std::vector<std::string> nameAndValue=l1menu::tools::splitByDelimeters( parameterListSplitByCommas[index], "=" );
76	+	if( nameAndValue.size()!=2 ) throw std::runtime_error( "TriggerRatePlot cannot be loaded from the supplied histogram because the title is not as expected" );
77	+
78	+	// The value could be scaled against the main threshold, so I need to check for that.
79	+	std::string valueAsString=nameAndValue[1];
80	+	if( valueAsString.size()>2 )
81	+	{
82	+	if( valueAsString[0]=='x' && valueAsString[1]=='*' )
83	+	{
84	+	// Parameter is scaled with the main threshold. Since I've set that to 1
85	+	// I can just use the number after the "x*" as is, so just chop off the
86	+	// first two characters.
87	+	valueAsString=valueAsString.substr( 2 );
88	+	scaledParameters.push_back( nameAndValue[0] );
89	+	}
90	+	}
91	+
92	+	pTemporaryTriggerCopy->parameter( nameAndValue[0] )=l1menu::tools::convertStringToFloat( valueAsString );
93	+	}
94	+
95	+	// Everything seems to have worked okay, so I can take a copy of the histogram
96	+	// and delegate to the normal initialisation routine.
97	+	pHistogram_.reset( static_cast<TH1*>( pPreExisitingHistogram->Clone() ) );
98	+	pHistogram_->SetDirectory( NULL ); // Hold in memory only
99	+	initiate( *pTemporaryTriggerCopy, scaledParameters );
100	+	}
101	+
102	+	l1menu::TriggerRatePlot::TriggerRatePlot( l1menu::TriggerRatePlot&& otherTriggerRatePlot ) noexcept
103	+	: pTrigger_( std::move(otherTriggerRatePlot.pTrigger_) ), pHistogram_( std::move(otherTriggerRatePlot.pHistogram_) ),
104	+	versusParameter_( std::move(otherTriggerRatePlot.versusParameter_) ), pParameter_(&pTrigger_->parameter(versusParameter_)),
105	+	otherParameterScalings_( std::move(otherTriggerRatePlot.otherParameterScalings_) ), histogramOwnedByMe_(histogramOwnedByMe_)
106	+	{
107	+	// No operation besides the initaliser list
108	+	}
109	+
110	+	l1menu::TriggerRatePlot& l1menu::TriggerRatePlot::operator=( l1menu::TriggerRatePlot&& otherTriggerRatePlot ) noexcept
111	+	{
112	+	// Free up whatever was there before
113	+	pTrigger_.reset();
114	+	if( histogramOwnedByMe_ ) pHistogram_.reset();
115	+	else pHistogram_.release();
116	+
117	+	// Then copy from the other object
118	+	pTrigger_=std::move(otherTriggerRatePlot.pTrigger_);
119	+	pHistogram_=std::move(otherTriggerRatePlot.pHistogram_);
120	+	versusParameter_=std::move(otherTriggerRatePlot.versusParameter_);
121	+	pParameter_=&pTrigger_->parameter(versusParameter_);
122	+	otherParameterScalings_=std::move(otherTriggerRatePlot.otherParameterScalings_);
123	+	histogramOwnedByMe_=otherTriggerRatePlot.histogramOwnedByMe_;
124	+
125	+	return *this;
126	+	}
127	+
128	+	void l1menu::TriggerRatePlot::initiate( const l1menu::ITrigger& trigger, const std::vector<std::string>& scaledParameters )
129	+	{
130	+	// Before making any histograms make sure errors are done properly
131	+	TH1::SetDefaultSumw2();
132	+
133		// Take a copy of the trigger
134		l1menu::TriggerTable& table=l1menu::TriggerTable::instance();
135		pTrigger_=table.copyTrigger( trigger );
#	Line 36 | Line 151 | l1menu::TriggerRatePlot::TriggerRatePlot
151		// I want to make a note of the other parameters set for the trigger. As far as I know TH1
152		// has no way of adding annotations so I'll tag it on the end of the title.
153		std::stringstream description;
154	<	description << pTrigger_->name() << " rate versus " << versusParameter;
154	>	description << pTrigger_->name() << " rate versus " << versusParameter_;
155
156		// Loop over all of the parameters and add their values to the description
157		std::vector<std::string> parameterNames=pTrigger_->parameterNames();
#	Line 44 | Line 159 | l1menu::TriggerRatePlot::TriggerRatePlot
159		if( parameterNames.size()>1 ) description << ",";
160		for( std::vector<std::string>::const_iterator iName=parameterNames.begin(); iName!=parameterNames.end(); ++iName )
161		{
162	<	if( *iName==versusParameter ) continue; // Don't bother adding the parameter I'm plotting against
162	>	if( *iName==versusParameter_ ) continue; // Don't bother adding the parameter I'm plotting against
163
164		// First check to see if this is one of the parameters that are being scaled
165		if( std::find(scaledParameters.begin(),scaledParameters.end(),*iName)==scaledParameters.end() )
#	Line 63 | Line 178 | l1menu::TriggerRatePlot::TriggerRatePlot
178		description << "]";
179
180		pHistogram_->SetTitle( description.str().c_str() );
66	–
67	–	}
68	–
69	–	l1menu::TriggerRatePlot::TriggerRatePlot( l1menu::TriggerRatePlot&& otherTriggerRatePlot ) noexcept
70	–	: pTrigger_( std::move(otherTriggerRatePlot.pTrigger_) ), pHistogram_( std::move(otherTriggerRatePlot.pHistogram_) ),
71	–	versusParameter_( std::move(otherTriggerRatePlot.versusParameter_) ), pParameter_(&pTrigger_->parameter(versusParameter_)),
72	–	otherParameterScalings_( std::move(otherTriggerRatePlot.otherParameterScalings_) ), histogramOwnedByMe_(histogramOwnedByMe_)
73	–	{
74	–	// No operation besides the initaliser list
75	–	}
76	–
77	–	l1menu::TriggerRatePlot& l1menu::TriggerRatePlot::operator=( l1menu::TriggerRatePlot&& otherTriggerRatePlot ) noexcept
78	–	{
79	–	// Free up whatever was there before
80	–	pTrigger_.reset();
81	–	if( histogramOwnedByMe_ ) pHistogram_.reset();
82	–	else pHistogram_.release();
83	–
84	–	// Then copy from the other object
85	–	pTrigger_=std::move(otherTriggerRatePlot.pTrigger_);
86	–	pHistogram_=std::move(otherTriggerRatePlot.pHistogram_);
87	–	versusParameter_=std::move(otherTriggerRatePlot.versusParameter_);
88	–	pParameter_=&pTrigger_->parameter(versusParameter_);
89	–	otherParameterScalings_=std::move(otherTriggerRatePlot.otherParameterScalings_);
90	–	histogramOwnedByMe_=otherTriggerRatePlot.histogramOwnedByMe_;
91	–
92	–	return *this;
181		}
182
183		l1menu::TriggerRatePlot::~TriggerRatePlot()
#	Line 134 | Line 222 | void l1menu::TriggerRatePlot::addSample(
222		void l1menu::TriggerRatePlot::addEvent( const l1menu::IEvent& event, const std::unique_ptr<l1menu::ICachedTrigger>& pCachedTrigger, float weightPerEvent )
223		{
224		//
225	<	// Loop over all of the bins and fill it if that threshold would
226	<	// have passed the trigger.
225	>	// Use bisection to find the bin that passes the trigger and the one
226	>	// immediately after it that fails.
227		//
228	<	for( int binNumber=1; binNumber<pHistogram_->GetNbinsX(); ++binNumber )
229	<	{
230	<	(*pParameter_)=pHistogram_->GetBinLowEdge(binNumber);
228	>	size_t lowBin=1;
229	>	size_t highBin=pHistogram_->GetNbinsX();
230	>
231	>	//
232	>	// First need to perform a check that the first bin passes. If it doesn't then
233	>	// the histogram doesn't need filling at all and I can return.
234	>	//
235	>	(*pParameter_)=pHistogram_->GetBinLowEdge(lowBin);
236	>	// Scale accordingly any other parameters that should be scaled. Remember that
237	>	// in parameterScalingPair, 'first' is a pointer to the threshold to be changed
238	>	// and 'second' is the ratio of the first threshold it should be.
239	>	for( const auto& parameterScalingPair : otherParameterScalings_ ) *(parameterScalingPair.first)=parameterScalingPair.second*(*pParameter_);
240	>	if( !pCachedTrigger->apply(event) ) return;
241
242	<	// Scale accordingly any other parameters that should be scaled. Remember that
243	<	// in parameterScalingPair, 'first' is a pointer to the threshold to be changed
244	<	// and 'second' is the ratio of the first threshold it should be.
245	<	for( const auto& parameterScalingPair : otherParameterScalings_ ) *(parameterScalingPair.first)=parameterScalingPair.second*(*pParameter_);
242	>	//
243	>	// Also check the highest bin. If that passes then I just fill every bin,
244	>	// otherwise I need to find the point at which the trigger fails.
245	>	//
246	>	(*pParameter_)=pHistogram_->GetBinLowEdge(highBin);
247	>	for( const auto& parameterScalingPair : otherParameterScalings_ ) *(parameterScalingPair.first)=parameterScalingPair.second*(*pParameter_);
248
249	<	if( pCachedTrigger->apply(event) ) // If the event passes the trigger
249	>	if( pCachedTrigger->apply(event) ) lowBin=highBin;
250	>	else
251	>	{
252	>	while( highBin-lowBin>1 ) // Loop until I find two bins next to each other
253		{
254	<	pHistogram_->Fill( pHistogram_->GetBinCenter(binNumber), event.weight()*weightPerEvent );
254	>	size_t middleBin=(highBin+lowBin)/2;
255	>
256	>	(*pParameter_)=pHistogram_->GetBinLowEdge(middleBin);
257	>	for( const auto& parameterScalingPair : otherParameterScalings_ ) *(parameterScalingPair.first)=parameterScalingPair.second*(*pParameter_);
258	>
259	>	if( pCachedTrigger->apply(event) ) lowBin=middleBin;
260	>	else highBin=middleBin;
261		}
262	<	else break;
263	<	// Not sure if this "else break" is a good idea. I don't know if triggers
264	<	// will run their thresholds the other way. E.g. a trigger could require
265	<	// a minimum amount of energy in the event, in which case the higher
266	<	// bins will pass.
267	<	} // end of loop over histogram bins
262	>	}
263	>
264	>	//
265	>	// Now I know which bins need filling, loop over them and fill.
266	>	//
267	>	for( size_t binNumber=1; binNumber<=lowBin; ++binNumber )
268	>	{
269	>	pHistogram_->Fill( pHistogram_->GetBinCenter(binNumber), event.weight()*weightPerEvent );
270	>	}
271
272		}
273
#	Line 164 | Line 276 | const l1menu::ITrigger& l1menu::TriggerR
276		return *pTrigger_;
277		}
278
279	+	float l1menu::TriggerRatePlot::findThreshold( float targetRate ) const
280	+	{
281	+	//
282	+	// Loop over all of the bins in the plot and find the first one
283	+	// that is less than the requested rate.
284	+	//
285	+	int binNumber;
286	+	for( binNumber=1; binNumber<pHistogram_->GetNbinsX(); ++binNumber )
287	+	{
288	+	if( pHistogram_->GetBinContent(binNumber)<targetRate ) break;
289	+	}
290	+
291	+	//if( binNumber==pHistogram_->GetNbinsX() ) throw std::runtime_error( "l1menu::TriggerRatePlot::findThreshold() was unable to find a threshold for the given rate" );
292	+
293	+	// I now have the bin number of the bin after the point I'm looking
294	+	// for. Now do a linear fit to interpolate between the bins using the
295	+	// two bins for the point and the two after.
296	+	std::vector<int> binNumberForLinearFit;
297	+	binNumberForLinearFit.push_back( binNumber-2 );
298	+	binNumberForLinearFit.push_back( binNumber-1 );
299	+	binNumberForLinearFit.push_back( binNumber );
300	+	binNumberForLinearFit.push_back( binNumber+1 );
301	+
302	+	std::vector< std::pair<float,float> > dataPoints;
303	+	for( auto& number : binNumberForLinearFit )
304	+	{
305	+	// Make sure all of of the bin numbers are valid
306	+	if( number<1 ) number=1;
307	+	else if( number>=pHistogram_->GetNbinsX() ) number=pHistogram_->GetNbinsX()-1;
308	+	dataPoints.push_back( std::make_pair( pHistogram_->GetBinCenter(number), pHistogram_->GetBinContent(number) ));
309	+	}
310	+
311	+	// Now do a simple linear fit on the data points
312	+	std::pair<float,float> slopeAndIntercept=l1menu::tools::simpleLinearFit( dataPoints );
313	+
314	+	// The straight line formula "y=m*x+c" equates to
315	+	// "targetRate=slopeAndIntercept.first*threshold+slopeAndIntercept.second"
316	+	// so rearrange and return the interpolated threshold.
317	+	return (targetRate-slopeAndIntercept.second)/slopeAndIntercept.first;
318	+	}
319	+
320	+
321		TH1* l1menu::TriggerRatePlot::getPlot()
322		{
323		return pHistogram_.get();

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines