LCFIVertex  0.7.2
PlotProcessor.cc
1 
2 #include "PlotProcessor.h"
3 #include <iostream>
4 #include <fstream>
5 #include <string>
6 #include <sstream>
7 #include <vector>
8 #include <cmath>
9 #include <set>
10 #ifdef USEROOT
11 #include "TFile.h"
13 #include "TGraph.h"
14 #include "TH1.h"
15 #include "TCanvas.h"
16 #include "TMultiGraph.h"
18 #endif // of #ifdef USEROOT
19 
20 #include "EVENT/LCCollection.h"
21 #include "IMPL/ReconstructedParticleImpl.h"
22 #include "EVENT/LCParameters.h"
23 #include "EVENT/LCFloatVec.h"
24 #include "EVENT/LCIntVec.h"
25 
26 #include "util/inc/vector3.h"
27 #include "util/inc/util.h"
28 
29 using std::includes;
30 using std::map;
31 using std::string;
32 using std::endl;
33 using std::cout;
34 using std::vector;
35 using std::stringstream;
36 using std::abs;
37 using std::pair;
38 using std::ofstream;
39 using std::cerr;
40 using std::set;
41 using namespace lcio;
42 
43 //Needs to be instantiated for Marlin to know about it (I think)
44 PlotProcessor aPlotProcessor;
45 
46 PlotProcessor::PlotProcessor() : marlin::Processor("Plot")
47 {
48  _description = "Plots various outputs from the flavour tag" ;
49 
50  // register steering parameters: name, description, class-variable, default value
51  //The name of the collection of ReconstructedParticles that is the jet
52  registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE,
53  "JetCollectionName" ,
54  "Name of the collection of ReconstructedParticles that is the jet" ,
55  _JetCollectionName ,
56  string("SGVJets") ) ;
57  vector<string> FlavourTagCollectionNamesDefault;
58  FlavourTagCollectionNamesDefault.push_back("FlavourTag");
59  registerProcessorParameter("FlavourTagCollections" ,
60  "Name of the LCFloatVec Collections that contain the flavour tags (one purity efficiency plot per tag) (in same order as jet collection)" ,
61  _FlavourTagCollectionNames,
62  FlavourTagCollectionNamesDefault) ;
63  registerInputCollection( LCIO::LCFLOATVEC,
64  "TrueJetFlavourCollection" ,
65  "Name of the output collection of LCIntVec (same order as jets)" ,
66  _TrueJetFlavourColName ,
67  std::string("TrueJetFlavour") ) ;
68  //The output filename
69  registerProcessorParameter( "OutputFilename" ,
70  "Filename for the output" ,
71  _OutputFilename ,
72  string("PlotProcessorOutput") ) ;
73 }
74 
75 PlotProcessor::~PlotProcessor()
76 {
77 }
78 
79 void PlotProcessor::init()
80 {
81  printParameters();
82  cout << _description << endl
83  << "-------------------------------------------------" << endl
84  << endl;
85 
86  _nRun=0;
87  _jetEMax=0;
88  //Make as many binning containers as we need
89  for (unsigned int iTag=0; iTag < _FlavourTagCollectionNames.size(); ++iTag)
90  {
94  }
95 }
96 
97 void PlotProcessor::processRunHeader( LCRunHeader* pRun )
98 {
99  //
100  // Perform a check to see if the variable names we need are here
101  //
102  for (unsigned int iTag=0; iTag < _FlavourTagCollectionNames.size(); ++iTag)
103  {
104  vector<string> VarNames;
105  (pRun->parameters()).getStringVals(_FlavourTagCollectionNames[iTag],VarNames);
106  //Fill the map realting names and indexes
107  set<string> AvailableNames;
108  map<string,unsigned int> IndexOf;
109  for (size_t i = 0;i < VarNames.size();++i)
110  {
111  AvailableNames.insert(VarNames[i]);
112  IndexOf[VarNames[i]] = i;
113  }
114 
115  //Add the index to the list
116  _IndexOfForEachTag.push_back(IndexOf);
117 
118  //Check the required information is in the LCFloatVec
119  set<string> RequiredNames;
120  RequiredNames.insert("BTag");
121  RequiredNames.insert("CTag");
122  RequiredNames.insert("BCTag");
123 
124  if (!includes(AvailableNames.begin(),AvailableNames.end(),RequiredNames.begin(),RequiredNames.end()))
125  cerr << _FlavourTagCollectionNames[iTag] << " does not contain information required by PlotProcessor";
126  }
127 }
128 
129 void PlotProcessor::processEvent( LCEvent* pEvent )
130 {
131  //Get the collection of jets. Can't do anything if the collection isn't there
132  //so don't bother catching the exception and terminate.
133  LCCollection* pJetCollection=pEvent->getCollection( _JetCollectionName );
134 
135  //make sure the collection is of the right type
136  if( pJetCollection->getTypeName()!=LCIO::RECONSTRUCTEDPARTICLE )
137  {
138  stringstream message;
139  message << endl
140  << "########################################################################################\n"
141  << "# PlotProcessor - #\n"
142  << "# The jet collection requested (\"" << _JetCollectionName << "\") is not of the type \"" << LCIO::RECONSTRUCTEDPARTICLE << "\" #\n"
143  << "########################################################################################" << endl;
144  throw EventException( message.str() );
145  }
146 
147  //apply any cuts on the event here
148  if( _passesEventCuts(pEvent) )
149  {
150  ReconstructedParticle* pJet;
151  //loop over the jets
152  for( int a=0; a<pJetCollection->getNumberOfElements(); ++a )
153  {
154  //Dynamic casts are not the best programming practice in the world, but I can't see another way of doing this
155  //in the LCIO framework. This cast should be safe though because we've already tested the type.
156  pJet=dynamic_cast<ReconstructedParticle*>( pJetCollection->getElementAt(a) );
157 
158  if( _passesJetCuts(pJet) )
159  {
160  _fillPlots(pEvent, a );
161  _jetEnergy.add(pJet->getEnergy());
162  if (pJet->getEnergy() > _jetEMax) _jetEMax = pJet->getEnergy();
163  }
164  }
165  }
166 }
167 
168 
169 void PlotProcessor::end()
170 {
172 
173  cout << "The largest jet energy was " << _jetEMax << endl;
174 }
175 
176 // IMPORTANT - If you change the cuts make sure you change the line below to show the changes in the docs
178 bool PlotProcessor::_passesEventCuts( LCEvent* pEvent )
179 {
180  //
181  // No event cuts at present
182  //
183 
184  return true;
185 }
186 
187 // IMPORTANT - If you change the cuts make sure you change the line below to show the changes in the docs
190 bool PlotProcessor::_passesJetCuts( ReconstructedParticle* pJet )
191 {
192  //
193  // This cut added on the suggestion of Sonja Hillert 12/Jan/07.
194  //
195  // Selects jets for which the cosine of the jet polar
196  // angle theta for all jets is not too large.
197  //
198  // Make something that's easy to search for to track down erroneous cuts:
199  // GREPTAG_CUT : Jet cut on abs(cos(theta)) of jet axis
200  //
201 
202  double CThJ_lower=0; //lower cut on cos(theta) of the jet axis
203  double CThJ_upper=0.95; //upper cut (Sho Ryu Ken!)
204 
205  vertex_lcfi::util::Vector3 zAxis( 0, 0, 1 ); //work out theta from dot product with jet axis
206 
207  const double* mom=pJet->getMomentum();
208  vertex_lcfi::util::Vector3 jetMomentum( mom[0], mom[1], mom[2] );
209 
210  jetMomentum.makeUnit();
211  double cosTheta=jetMomentum.dot( zAxis );
212  if( fabs(cosTheta)<=CThJ_lower || fabs(cosTheta)>=CThJ_upper ) return false;
213 
214 
215  // If control gets to this point then the jet has passed
216  return true;
217 }
218 
219 void PlotProcessor::_fillPlots( LCEvent* pEvent, unsigned int jet)
220 {
221  LCCollection* pTrueCollection=pEvent->getCollection( _TrueJetFlavourColName );
222  float jetType = *(dynamic_cast<LCFloatVec*>(pTrueCollection->getElementAt(jet))->begin());
223  //LCCollection* pTrueCollection=pEvent->getCollection( "SGVFlavourTagInputs" );
224  //int jetType=static_cast<int>((*dynamic_cast<LCFloatVec*>(pTrueCollection->getElementAt(jet)))[19]);
225  for (unsigned int iTag=0; iTag < _FlavourTagCollectionNames.size(); ++iTag)
226  {
227  LCCollection* pTagCollection=pEvent->getCollection( _FlavourTagCollectionNames[iTag] );
228  double bTag= (*dynamic_cast<LCFloatVec*>(pTagCollection->getElementAt(jet)))[_IndexOfForEachTag[iTag]["BTag"]];
229  double cTag= (*dynamic_cast<LCFloatVec*>(pTagCollection->getElementAt(jet)))[_IndexOfForEachTag[iTag]["CTag"]];
230  double cTagBBack= (*dynamic_cast<LCFloatVec*>(pTagCollection->getElementAt(jet)))[_IndexOfForEachTag[iTag]["BCTag"]];
231 
232  if( jetType==B_JET )
233  {
234  if( bTag<=1 && bTag>=0 ) _BTagEfficiencyPurity[iTag].add_signal( bTag );
235  if( cTag<=1 && cTag>=0 ) _CTagEfficiencyPurity[iTag].add_background( cTag );
236  if( cTagBBack<=1 && cTagBBack>=0 ) _BCTagEfficiencyPurity[iTag].add_background( cTagBBack );
237  }
238  else if( jetType==C_JET )
239  {
240  if( bTag<=1 && bTag>=0 ) _BTagEfficiencyPurity[iTag].add_background( bTag );
241  if( cTag<=1 && cTag>=0 ) _CTagEfficiencyPurity[iTag].add_signal( cTag );
242  if( cTagBBack<=1 && cTagBBack>=0 ) _BCTagEfficiencyPurity[iTag].add_signal( cTagBBack );
243  }
244  else
245  {
246  if( bTag<=1 && bTag>=0 ) _BTagEfficiencyPurity[iTag].add_background( bTag );
247  if( cTag<=1 && cTag>=0 ) _CTagEfficiencyPurity[iTag].add_background( cTag );
248  //don't add background for the bcnet result because we're only considering the b background
249  }
250  }
251 }
252 
253 void PlotProcessor::_outputDataToFile( string filename )
254 {
255 #ifdef USEROOT
256  stringstream filenameStream;
257  filenameStream << filename << ".root";
258  TFile rootFile( filenameStream.str().c_str(), "RECREATE", "Various plots from the plot processor" );
259 
260  //have to put all the results into normal c arrays for root to deal with
261  //there may be an easier way of doing this conversion but I don't know what it is
262 
263  for (unsigned int iTag=0; iTag < _FlavourTagCollectionNames.size(); ++iTag)
264  {
265  int numberOfPoints=100;
266 
267  //get the results in STL form
268  vector< pair<double,double> > BTagefficiencyPurityPairs=_BTagEfficiencyPurity[iTag].eff_pur( numberOfPoints );
269  vector< pair<double,double> > CTagefficiencyPurityPairs=_CTagEfficiencyPurity[iTag].eff_pur( numberOfPoints );
270  vector< pair<double,double> > BCTagefficiencyPurityPairs=_BCTagEfficiencyPurity[iTag].eff_pur( numberOfPoints );
271 
272  //have to put all the results into normal c arrays for root to deal with
273  //there may be an easier way of doing this conversion but I don't know what it is
274  double BNetEfficiency[100];
275  double BNetPurity[100];
276  double CNetEfficiency[100];
277  double CNetPurity[100];
278  double BCNetEfficiency[100];
279  double BCNetPurity[100];
280 
281  for( int a=0; a<numberOfPoints; ++a )
282  {
283  BNetEfficiency[a]=BTagefficiencyPurityPairs[a].first;
284  BNetPurity[a]=BTagefficiencyPurityPairs[a].second;
285  CNetEfficiency[a]=CTagefficiencyPurityPairs[a].first;
286  CNetPurity[a]=CTagefficiencyPurityPairs[a].second;
287  BCNetEfficiency[a]=BCTagefficiencyPurityPairs[a].first;
288  BCNetPurity[a]=BCTagefficiencyPurityPairs[a].second;
289  }
290 
291  TGraph BNetGraph( 99, BNetEfficiency, BNetPurity );
292  BNetGraph.SetName( (_FlavourTagCollectionNames[iTag] + " B Tag").c_str() );
293  TGraph CNetGraph( 99, CNetEfficiency, CNetPurity );
294  CNetGraph.SetName( (_FlavourTagCollectionNames[iTag] + " C Tag").c_str() );
295  TGraph BCNetGraph( 99, BCNetEfficiency, BCNetPurity );
296  BCNetGraph.SetName( (_FlavourTagCollectionNames[iTag] + " BC Tag").c_str() );
297 
298  TMultiGraph TagGraphs;
299  TagGraphs.SetName( (_FlavourTagCollectionNames[iTag]).c_str() );
300  TagGraphs.Add( &BNetGraph );
301  TagGraphs.Add( &CNetGraph );
302  TagGraphs.Add( &BCNetGraph );
303 
304  TAxis* pAxis;
305  TagGraphs.SetTitle( (string("Efficiency-purity for ") + _FlavourTagCollectionNames[iTag]).c_str() );
306 
307  //why the hell doesn't this work?
308  pAxis=TagGraphs.GetXaxis();
309  if(pAxis) pAxis->SetTitle( "Efficiency" );
310 
311  pAxis=TagGraphs.GetYaxis();
312  if(pAxis) pAxis->SetTitle( "Purity" );
313 
314  rootFile.Add( &TagGraphs );
315  rootFile.Write();
316  }
317  // now add in the jet energies
318  TH1F jetEnergyHistogram( "jetEnergyHistogram", "Jet energies", 200, 0, 120 );
319  const vector<double> jetEnergyData=_jetEnergy.sorted_data();
320  for( vector<double>::const_iterator i=jetEnergyData.begin(); i<jetEnergyData.end(); ++i ) jetEnergyHistogram.Fill( (*i) );
321 
322  rootFile.Write();
323  rootFile.Close();
324 
325 #else //of ifdef USEROOT
326 
327  //
328  // Not using root so output everything to comma seperated value files (csv)
329  //
330  ofstream ofile;
331  ofile.open( (filename + ".csv").c_str() );
332  if( ofile.is_open() )
333  {
334  for (unsigned int iTag=0; iTag < _FlavourTagCollectionNames.size(); ++iTag)
335  {
336  string n = _FlavourTagCollectionNames[iTag];
337  ofile << "E," << n << " B Tag,E," << n << " C Tag,E," << n << " BC Tag,";
338  }
339  ofile << endl;
340 
341  vector<vector<pair<double,double> > > BTagefficiencyPurityPairs;
342  vector<vector<pair<double,double> > > CTagefficiencyPurityPairs;
343  vector<vector<pair<double,double> > > BCTagefficiencyPurityPairs;
344 
345  int numberOfPoints=100;
346  for (unsigned int iTag=0; iTag < _FlavourTagCollectionNames.size(); ++iTag)
347  {
348  BTagefficiencyPurityPairs.push_back(_BTagEfficiencyPurity[iTag].eff_pur( numberOfPoints ));
349  CTagefficiencyPurityPairs.push_back(_CTagEfficiencyPurity[iTag].eff_pur( numberOfPoints ));
350  BCTagefficiencyPurityPairs.push_back(_BCTagEfficiencyPurity[iTag].eff_pur( numberOfPoints ));
351  // std::cout << _BTagEfficiencyPurity[iTag].number_of_signal() << " " << _BTagEfficiencyPurity[iTag].number_of_background() << std::endl;
352  // std::cout << _CTagEfficiencyPurity[iTag].number_of_signal() << " " << _CTagEfficiencyPurity[iTag].number_of_background() << std::endl;
353  // std::cout << _BCTagEfficiencyPurity[iTag].number_of_signal() << " " << _BCTagEfficiencyPurity[iTag].number_of_background() << std::endl;
354  }
355  for( int a=0; a<numberOfPoints; ++a )
356  {
357  for (unsigned int iTag=0; iTag < _FlavourTagCollectionNames.size(); ++iTag)
358  {
359  ofile << BTagefficiencyPurityPairs[iTag][a].first << "," << BTagefficiencyPurityPairs[iTag][a].second << ",";
360  ofile << CTagefficiencyPurityPairs[iTag][a].first << "," << CTagefficiencyPurityPairs[iTag][a].second << ",";
361  ofile << BCTagefficiencyPurityPairs[iTag][a].first << "," << BCTagefficiencyPurityPairs[iTag][a].second << ",";
362  }
363  ofile << endl;
364  }
365  ofile.close();
366  }
367  else
368  {
369  cerr << "########################################################################################\n"
370  << "# PlotProcessor - #\n"
371  << "# Unable to open file \"" << filename << "\" for output #\n"
372  << "########################################################################################" << endl;
373  }
374 
375  stringstream filenameStream;
376  filenameStream << filename << "-JetEnergies.csv";
377 
378  ofile.open( filenameStream.str().c_str() );
379  if( ofile.is_open() )
380  {
381  vector< vertex_lcfi::util::bin<double> > binned=_jetEnergy.binned_data( 200, 0, 120 );
382  ofile << "Bin low,Bin High,Frequency" << endl;
383  for( vector< vertex_lcfi::util::bin<double> >::iterator i=binned.begin(); i<binned.end(); ++i )
384  ofile << (*i).region_low() << "," << (*i).region_high() << "," << (*i).contents() << endl;
385  ofile.close();
386  }
387  else
388  {
389  cerr << "########################################################################################\n"
390  << "# PlotProcessor - #\n"
391  << "# Unable to open file \"" << filenameStream.str() << "\" for output #\n"
392  << "########################################################################################" << endl;
393  }
394 #endif //of ifdef USEROOT
395 }
396 
398 {
399  const vector<string>* pCollectionNames=pEvent->getCollectionNames();
400 
401  cout << "The available collections are: (name - type)" << endl;
402  for( vector<string>::const_iterator i=pCollectionNames->begin(); i<pCollectionNames->end(); ++i )
403  {
404  LCCollection* pCollection=pEvent->getCollection( (*i) );
405  const string typeName=pCollection->getTypeName();
406  cout << " " << (*i) << " - " << typeName << endl;
407  }
408  cout << endl;
409 }
std::vector< efficiency_purity< double > > _BCTagEfficiencyPurity
Definition: PlotProcessor.h:86
std::vector< efficiency_purity< double > > _CTagEfficiencyPurity
Definition: PlotProcessor.h:85
std::vector< efficiency_purity< double > > _BTagEfficiencyPurity
Definition: PlotProcessor.h:84
bool _passesJetCuts(lcio::ReconstructedParticle *pJet)
void _fillPlots(LCEvent *pEvent, unsigned int jet)
std::vector< std::string > _FlavourTagCollectionNames
Definition: PlotProcessor.h:76
static const int B_JET
Definition: PlotProcessor.h:90
std::string _TrueJetFlavourColName
Definition: PlotProcessor.h:77
static const int C_JET
Definition: PlotProcessor.h:89
std::string _OutputFilename
Definition: PlotProcessor.h:78
Creates some sample plots from the data calculated by the LCFI vertex package.
Definition: PlotProcessor.h:61
void _displayCollectionNames(lcio::LCEvent *pEvent)
histogram_data< double > _jetEnergy
Definition: PlotProcessor.h:82
std::string _JetCollectionName
Definition: PlotProcessor.h:75
bool _passesEventCuts(lcio::LCEvent *pEvent)
void _outputDataToFile(std::string filename)