ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/TriggerRatePlot.cpp
Revision: 1.10
Committed: Wed Jul 24 11:48:55 2013 UTC (11 years, 9 months ago) by grimes
Branch: MAIN
CVS Tags: HEAD
Changes since 1.9: +201 -47 lines
Log Message:
Big commit of lots of changes before migrating to Git.

File Contents

# User Rev Content
1 grimes 1.1 #include "l1menu/TriggerRatePlot.h"
2    
3     #include "l1menu/ITrigger.h"
4 grimes 1.5 #include "l1menu/ICachedTrigger.h"
5     #include "l1menu/L1TriggerDPGEvent.h"
6 grimes 1.1 #include "l1menu/IEvent.h"
7 grimes 1.8 #include "l1menu/ISample.h"
8 grimes 1.1 #include "l1menu/TriggerTable.h"
9 grimes 1.10 #include "l1menu/tools/tools.h"
10     #include "l1menu/tools/stringManipulation.h"
11 grimes 1.1 #include <TH1F.h>
12     #include <sstream>
13 grimes 1.6 #include <algorithm>
14 grimes 1.10 #include <stdexcept>
15 grimes 1.6 #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 grimes 1.1 : pHistogram_( std::move(pHistogram) ), versusParameter_(versusParameter), histogramOwnedByMe_(true)
19     {
20 grimes 1.10 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 grimes 1.1 // 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 grimes 1.6 // throw an exception. Take a pointer to the parameter so I don't need to keep performing
139     // expensive string comparisons.
140 grimes 1.2 pParameter_=&pTrigger_->parameter(versusParameter_);
141 grimes 1.1
142 grimes 1.6 // 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 grimes 1.1 // 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 grimes 1.10 description << pTrigger_->name() << " rate versus " << versusParameter_;
155 grimes 1.1
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 grimes 1.10 if( *iName==versusParameter_ ) continue; // Don't bother adding the parameter I'm plotting against
163 grimes 1.6
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 grimes 1.1 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,
186     // then I need to tell the unique_ptr not to delete it.
187     if( !histogramOwnedByMe_ )
188     {
189     pHistogram_.release();
190     }
191     }
192    
193 grimes 1.7 void l1menu::TriggerRatePlot::addEvent( const l1menu::IEvent& event )
194 grimes 1.1 {
195 grimes 1.7 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 grimes 1.9 // 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 grimes 1.7 std::unique_ptr<l1menu::ICachedTrigger> pCachedTrigger=sample.createCachedTrigger( *pTrigger_ );
203    
204 grimes 1.9 addEvent( event, pCachedTrigger, weightPerEvent );
205 grimes 1.1 }
206    
207 grimes 1.5 void l1menu::TriggerRatePlot::addSample( const l1menu::ISample& sample )
208     {
209     float weightPerEvent=sample.eventRate()/sample.sumOfWeights();
210    
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 grimes 1.9 addEvent( sample.getEvent(eventNumber), pCachedTrigger, weightPerEvent );
218     } // end of loop over events
219    
220     }
221 grimes 1.5
222 grimes 1.9 void l1menu::TriggerRatePlot::addEvent( const l1menu::IEvent& event, const std::unique_ptr<l1menu::ICachedTrigger>& pCachedTrigger, float weightPerEvent )
223     {
224     //
225 grimes 1.10 // 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 grimes 1.9 //
246 grimes 1.10 (*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 grimes 1.9 {
252 grimes 1.10 while( highBin-lowBin>1 ) // Loop until I find two bins next to each other
253     {
254     size_t middleBin=(highBin+lowBin)/2;
255 grimes 1.6
256 grimes 1.10 (*pParameter_)=pHistogram_->GetBinLowEdge(middleBin);
257     for( const auto& parameterScalingPair : otherParameterScalings_ ) *(parameterScalingPair.first)=parameterScalingPair.second*(*pParameter_);
258 grimes 1.6
259 grimes 1.10 if( pCachedTrigger->apply(event) ) lowBin=middleBin;
260     else highBin=middleBin;
261 grimes 1.9 }
262 grimes 1.10 }
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 grimes 1.6
272 grimes 1.5 }
273    
274 grimes 1.3 const l1menu::ITrigger& l1menu::TriggerRatePlot::getTrigger() const
275     {
276     return *pTrigger_;
277     }
278    
279 grimes 1.10 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 grimes 1.1 TH1* l1menu::TriggerRatePlot::getPlot()
322     {
323     return pHistogram_.get();
324     }
325    
326     TH1* l1menu::TriggerRatePlot::relinquishOwnershipOfPlot()
327     {
328     // Rather than call release() on the unique_ptr, I'll set a flag so that I know
329     // to release() in the destructor instead. This way it's possible to relinquish
330     // ownership but still perform operations on the histograms.
331     histogramOwnedByMe_=false;
332     return pHistogram_.get();
333     }