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.1 by grimes, Tue May 28 23:14:03 2013 UTC vs.
Revision 1.10 by grimes, Wed Jul 24 11:48:55 2013 UTC

#	Line 1 | Line 1
1		#include "l1menu/TriggerRatePlot.h"
2
3		#include "l1menu/ITrigger.h"
4	+	#include "l1menu/ICachedTrigger.h"
5	+	#include "l1menu/L1TriggerDPGEvent.h"
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>
7	–	#include <iostream>
12		#include <sstream>
13	+	#include <algorithm>
14	+	#include <stdexcept>
15	+	#include <iostream>
16
17	<	l1menu::TriggerRatePlot::TriggerRatePlot( const l1menu::ITrigger& trigger, std::unique_ptr<TH1> pHistogram, const std::string versusParameter )
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	<	// Take a copy of the trigger
21	<	l1menu::TriggerTable& table=l1menu::TriggerTable::instance();
15	<	pTrigger_=table.copyTrigger( trigger );
16	<
17	<	// Make sure the versusParameter_ supplied is valid. If it's not then this call will
18	<	// throw an exception.
19	<	pTrigger_->parameter(versusParameter_);
20	>	initiate( trigger, scaledParameters );
21	>	}
22
23	<	// I want to make a note of the other parameters set for the trigger. As far as I know TH1
24	<	// has no way of adding annotations so I'll tag it on the end of the title.
25	<	std::stringstream description;
26	<	description << pTrigger_->name() << " rate versus " << versusParameter;
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	<	// Loop over all of the parameters and add their values to the description
31	<	std::vector<std::string> parameterNames=pTrigger_->parameterNames();
32	<	description << " [v" << pTrigger_->version();
33	<	if( parameterNames.size()>1 ) description << ",";
34	<	for( std::vector<std::string>::const_iterator iName=parameterNames.begin(); iName!=parameterNames.end(); ++iName )
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	<	if( *iName==versusParameter ) continue; // Don't bother adding the parameter I'm plotting against
74	<	description << *iName << "=" << pTrigger_->parameter(*iName);
75	<	if( iName+1!=parameterNames.end() ) description << ","; // Add delimeter between parameter names
76	<	}
77	<	description << "]";
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	<	pHistogram_->SetTitle( description.str().c_str() );
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_) ), histogramOwnedByMe_(histogramOwnedByMe_)
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		}
#	Line 57 | Line 118 | l1menu::TriggerRatePlot& l1menu::Trigger
118		pTrigger_=std::move(otherTriggerRatePlot.pTrigger_);
119		pHistogram_=std::move(otherTriggerRatePlot.pHistogram_);
120		versusParameter_=std::move(otherTriggerRatePlot.versusParameter_);
121	<	histogramOwnedByMe_=histogramOwnedByMe_;
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 );
136	+
137	+	// Make sure the versusParameter_ supplied is valid. If it's not then this call will
138	+	// throw an exception. Take a pointer to the parameter so I don't need to keep performing
139	+	// expensive string comparisons.
140	+	pParameter_=&pTrigger_->parameter(versusParameter_);
141	+
142	+	// If any parameters have been requested to be scaled along with the versusParameter, figure
143	+	// out what the scaling should be and take a note of pointers.
144	+	for( const auto& parameterToScale : scaledParameters )
145	+	{
146	+	if( parameterToScale!=versusParameter_ )
147	+	{
148	+	otherParameterScalings_.push_back( std::make_pair( &pTrigger_->parameter(parameterToScale), pTrigger_->parameter(parameterToScale)/(*pParameter_) ) );
149	+	}
150	+	}
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_;
155	+
156	+	// Loop over all of the parameters and add their values to the description
157	+	std::vector<std::string> parameterNames=pTrigger_->parameterNames();
158	+	description << " [v" << pTrigger_->version();
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
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() )
166	+	{
167	+	// This parameter isn't being scaled, so write the absoulte value in the description
168	+	description << *iName << "=" << pTrigger_->parameter(*iName);
169	+	}
170	+	else
171	+	{
172	+	// This parameter is being scaled, so write what the scaling is in the description
173	+	description << *iName << "=x*" << pTrigger_->parameter(*iName)/(*pParameter_);
174	+	}
175	+
176	+	if( iName+1!=parameterNames.end() ) description << ","; // Add delimeter between parameter names
177	+	}
178	+	description << "]";
179	+
180	+	pHistogram_->SetTitle( description.str().c_str() );
181	+	}
182	+
183		l1menu::TriggerRatePlot::~TriggerRatePlot()
184		{
185		// If the flag has been set telling me that this instance doesn't own the histogram,
#	Line 72 | Line 190 | l1menu::TriggerRatePlot::~TriggerRatePlo
190		}
191		}
192
193	<	void l1menu::TriggerRatePlot::addEvent( const l1menu::IEvent& event ) const
193	>	void l1menu::TriggerRatePlot::addEvent( const l1menu::IEvent& event )
194	>	{
195	>	const l1menu::ISample& sample=event.sample();
196	>	float weightPerEvent=sample.eventRate()/sample.sumOfWeights();
197	>
198	>	// For some implementations of ISample, it is significantly faster to
199	>	// create ICachedTriggers and then loop over those. The addEvent overload
200	>	// I'm about to delegate to loops over the histogram bins, so even for
201	>	// one event this trigger can be called multiple times.
202	>	std::unique_ptr<l1menu::ICachedTrigger> pCachedTrigger=sample.createCachedTrigger( *pTrigger_ );
203	>
204	>	addEvent( event, pCachedTrigger, weightPerEvent );
205	>	}
206	>
207	>	void l1menu::TriggerRatePlot::addSample( const l1menu::ISample& sample )
208		{
209	<	float& threshold1=pTrigger_->parameter(versusParameter_);
209	>	float weightPerEvent=sample.eventRate()/sample.sumOfWeights();
210
211	<	// loop over all of the bins in the histogram, and see if the trigger passes
212	<	// for the value at the center of the bin.
213	<	for( int binNumber=1; binNumber<pHistogram_->GetNbinsX(); ++binNumber )
211	>	// Create a cached trigger, which depending on the concrete type of the ISample
212	>	// may or may not significantly increase the speed at which this next loop happens.
213	>	std::unique_ptr<l1menu::ICachedTrigger> pCachedTrigger=sample.createCachedTrigger( *pTrigger_ );
214	>
215	>	for( size_t eventNumber=0; eventNumber<sample.numberOfEvents(); ++eventNumber )
216		{
217	<	threshold1=pHistogram_->GetBinCenter(binNumber);
218	<	if( pTrigger_->apply( event ) )
217	>	addEvent( sample.getEvent(eventNumber), pCachedTrigger, weightPerEvent );
218	>	} // end of loop over events
219	>
220	>	}
221	>
222	>	void l1menu::TriggerRatePlot::addEvent( const l1menu::IEvent& event, const std::unique_ptr<l1menu::ICachedTrigger>& pCachedTrigger, float weightPerEvent )
223	>	{
224	>	//
225	>	// Use bisection to find the bin that passes the trigger and the one
226	>	// immediately after it that fails.
227	>	//
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	>	//
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) ) lowBin=highBin;
250	>	else
251	>	{
252	>	while( highBin-lowBin>1 ) // Loop until I find two bins next to each other
253		{
254	<	pHistogram_->Fill( threshold1, event.weight() );
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		}
88	–	// could put an "else break" in here, but I don't know if triggers
89	–	// will run their thresholds the other way. E.g. a trigger could require
90	–	// a minimum amount of energy in the event, in which case the higher
91	–	// bins will pass.
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	+
274	+	const l1menu::ITrigger& l1menu::TriggerRatePlot::getTrigger() const
275	+	{
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