LCFIVertex  0.7.2
DSTAIDAPlotProcessor.cc
1 // If AIDA is being used then this will have been defined
2 #ifndef MARLIN_USE_AIDA
3 
4 #warning "-------------------------------------------------------------------------------"
5 #warning "- DSTAIDAPlotProcessor requires MARLIN_USE_AIDA to be defined. Did you enable -"
6 #warning "- AIDA when compiling Marlin? DSTAIDAPlotProcessor will not be compiled. -"
7 #warning "-------------------------------------------------------------------------------"
8 
9 // Can't do anything else.
10 #else
11 
12 #include "DSTAIDAPlotProcessor.h"
13 #include <iostream>
14 #include <fstream>
15 #include <string>
16 #include <sstream>
17 #include <vector>
18 #include <cmath>
19 #include <set>
20 
21 
22 // AIDA includes...
23 #include <marlin/AIDAProcessor.h>
24 #include <AIDA/IMeasurement.h>
25 #include <AIDA/IHistogramFactory.h>
26 #include <AIDA/IHistogram1D.h>
27 #include <AIDA/IHistogram2D.h>
28 #include <AIDA/IAxis.h>
29 #include <AIDA/ITree.h>
30 #include <AIDA/ITupleFactory.h>
31 #include <AIDA/ICloud2D.h>
32 #include <UTIL/LCRelationNavigator.h>
33 
34 
35 //change USING_RAIDA to USING_JAIDA if you are using JAIDA/AIDAJNI - you will obtain more functionality!
36 #define USING_RAIDA
37 
38 #ifdef USING_RAIDA
39 #warning "USING_RAIDA defined"
40 #else
41 #define USING_JAIDA
42 #warning "USING_JAIDA defined"
43 #endif
44 
45 
46 #ifdef USING_JAIDA//Data point sets aren't implemented in RAIDA - which is a shame as they have functionality not given by histograms
47 //such as the facility to set the error
48 #include <AIDA/IDataPointSet.h>
49 #include <AIDA/IDataPointSetFactory.h>
50 #include <AIDA/IDataPoint.h>
51 #endif
52 
53 
54 
55 #include "EVENT/LCCollection.h"
56 #include "UTIL/PIDHandler.h"
57 #include "IMPL/ReconstructedParticleImpl.h"
58 #include "EVENT/LCParameters.h"
59 #include "EVENT/LCFloatVec.h"
60 #include "EVENT/LCIntVec.h"
61 
62 #include "util/inc/vector3.h"
63 #include "util/inc/util.h"
64 
65 using std::includes;
66 using std::map;
67 using std::string;
68 using std::endl;
69 using std::cout;
70 using std::vector;
71 using std::stringstream;
72 using std::abs;
73 using std::pair;
74 using std::ofstream;
75 using std::cerr;
76 using std::set;
77 using namespace lcio;
78 
79 #include "TypesafeCollection.h"
80 
81 //Needs to be instantiated for Marlin to know about it (I think)
82 DSTAIDAPlotProcessor aDSTAIDAPlotProcessor;
83 
84 DSTAIDAPlotProcessor::DSTAIDAPlotProcessor() : marlin::Processor("DSTAIDAPlotProcessor")
85 {
86  _description = "Plots various outputs from the flavour tag" ;
87 
88  // register steering parameters: name, description, class-variable, default value
89  //The name of the collection of ReconstructedParticles that is the jet
90  registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE,
91  "JetCollectionName" ,
92  "Name of the collection of ReconstructedParticles that is the jet" ,
93  _JetCollectionName ,
94  string("FTSelectedJets") ) ;
95 
96 
97 }
98 
99 DSTAIDAPlotProcessor::~DSTAIDAPlotProcessor()
100 {
101 }
102 
103 void DSTAIDAPlotProcessor::init()
104 {
105  printParameters();
106  cout << _description << endl
107  << "-------------------------------------------------" << endl
108  << endl;
109 
110  _nRun=0;
111 
112 
113  _VertexCatNames.resize(N_VERTEX_CATEGORIES+1);
114  _VertexCatNames[0]="AnyNumberOfVertices";
115  _VertexCatNames[1]="OneVertex";
116  _VertexCatNames[2]="TwoVertices";
117  _VertexCatNames[3]="ThreeOrMoreVertices";
118 
119  _NumVertexCatDir.resize(N_VERTEX_CATEGORIES+1);
120  _NumVertexCatDir[1]="OneVertex";
121  _NumVertexCatDir[2]="TwoVertices";
122  _NumVertexCatDir[3]="ThreeOrMoreVertices";
123  _NumVertexCatDir[0]="AnyNumberOfVertices";
124 
125  _numberOfPoints=100;
126 
127  AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
128  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
129 
130  if( pHistogramFactory!=0 )
131  {
132  if (!(pTree->cd( "/" + name() + "/"))) {
133  pTree->mkdir( "/" + name() + "/" );
134  pTree->cd( "/" + name() + "/");
135  }
136 
137 
138  CreateTagPlots();
139 
140  CreateFlavourTagTuple();
141 
142  } else {
143 
144  std::cerr << "### " << __FILE__ << "(" << __LINE__ << "): Unable to get the histogram factory! No histograms will be made."<< std::endl;
145  }
146 
147 
148 
149 }
150 
151 void DSTAIDAPlotProcessor::processRunHeader( LCRunHeader* pRun )
152 {
153 
154 
155  CreateFlavourTagInputPlots(pRun);
156 
157 }
158 
159 void DSTAIDAPlotProcessor::processEvent( LCEvent* pEvent )
160 {
161  //check that what we need is here
162  //get the jet collection
163  LCCollection* jetCollection;
164  try{
165  jetCollection=pEvent->getCollection( _JetCollectionName );
166  }
167  catch(DataNotAvailableException &e){
168  std::cout << "Collection " <<_JetCollectionName << " is unavailable, can't plot anything " << std::endl;
169  return;
170  }
171 
172  PIDHandler pidh( jetCollection ) ;
173  //get the algorithm id for the MC info
174  try{
175  pidh.getAlgorithmID("MCTruth");
176  }
177  catch(UnknownAlgorithm &e)
178  {
179  std::cout << "MCTruth Algorithm is unavailable, can't plot anything " << std::endl;
180  return;
181  }
182  //get the algorithm id for the MC info
183  try{
184  pidh.getAlgorithmID("LCFIFlavourTag");
185  }
186  catch(UnknownAlgorithm &e){
187  std::cout << "LCFIFlavourTag Algorithm is unavailable, can't plot anything " << std::endl;
188  return;
189  }
190 
191 
192  //apply any cuts on the event here
193  if( _passesEventCuts(pEvent) )
194  {
195 
196  ReconstructedParticle* pJet;
197  //loop over the jets
198  for( int jetNumber=0; jetNumber<jetCollection->getNumberOfElements(); ++jetNumber )
199  {
200  pJet=dynamic_cast<ReconstructedParticle*>( jetCollection->getElementAt(jetNumber) );
201 
202  //only do anything if the jet passes the jet cuts
203  if( _passesJetCuts(pJet) )
204  {
205  FillTagPlots( pEvent, jetNumber );
206  FillInputsPlots( pEvent, jetNumber );
207  }
208  }
209  }
210 
211 
212 
213 }
214 
215 
216 
217 
218 void DSTAIDAPlotProcessor::end()
219 {
220 
221  CalculateIntegralAndBackgroundPlots();
222  CalculateEfficiencyPurityPlots();
223 
224 }
225 
226 // IMPORTANT - If you change the cuts make sure you change the line below to show the changes in the docs
228 bool DSTAIDAPlotProcessor::_passesEventCuts( LCEvent* pEvent )
229 {
230 
231  return true;
232 }
233 
234 // IMPORTANT - If you change the cuts make sure you change the line below to show the changes in the docs
237 bool DSTAIDAPlotProcessor::_passesJetCuts( ReconstructedParticle* pJet )
238 {
239  //
240  // This cut added on the suggestion of Sonja Hillert 12/Jan/07.
241  //
242  // Selects jets for which the cosine of the jet polar
243  // angle theta for all jets is not too large.
244  //
245  // Make something that's easy to search for to track down erroneous cuts:
246  // GREPTAG_CUT : Jet cut on abs(cos(theta)) of jet axis
247  //
248 
249  double CThJ_lower=0; //lower cut on cos(theta) of the jet axis
250  double CThJ_upper=0.95; //upper cut (Sho Ryu Ken!)
251 
252  vertex_lcfi::util::Vector3 zAxis( 0, 0, 1 ); //work out theta from dot product with jet axis
253 
254  const double* mom=pJet->getMomentum();
255  vertex_lcfi::util::Vector3 jetMomentum( mom[0], mom[1], mom[2] );
256 
257  jetMomentum.makeUnit();
258  double cosTheta=jetMomentum.dot( zAxis );
259  if( fabs(cosTheta)<=CThJ_lower || fabs(cosTheta)>=CThJ_upper ) return false;
260 
261 
262  // If control gets to this point then the jet has passed
263  return true;
264 }
265 
266 
267 void DSTAIDAPlotProcessor::_displayCollectionNames( LCEvent* pEvent )
268 {
269  const vector<string>* pCollectionNames=pEvent->getCollectionNames();
270 
271  cout << "The available collections are: (name - type)" << endl;
272  for( vector<string>::const_iterator i=pCollectionNames->begin(); i<pCollectionNames->end(); ++i )
273  {
274  LCCollection* pCollection=pEvent->getCollection( (*i) );
275  const string typeName=pCollection->getTypeName();
276  cout << " " << (*i) << " - " << typeName << endl;
277  }
278  cout << endl;
279 }
280 
281 
282 void DSTAIDAPlotProcessor::CreateTagPlots()
283 {
284  AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
285  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
286 
287 
288 
289  for (unsigned int iVertexCat=0; iVertexCat < N_VERTEX_CATEGORIES+1; ++iVertexCat ){
290  _pLightJetBTag[_VertexCatNames[iVertexCat]]=0;
291  _pLightJetCTag[_VertexCatNames[iVertexCat]]=0;
292  _pLightJetBCTag[_VertexCatNames[iVertexCat]]=0;
293 
294  _pBJetBTag[_VertexCatNames[iVertexCat]]=0;
295  _pBJetCTag[_VertexCatNames[iVertexCat]]=0;
296  _pBJetBCTag[_VertexCatNames[iVertexCat]]=0;
297 
298  _pCJetBTag[_VertexCatNames[iVertexCat]]=0;
299  _pCJetCTag[_VertexCatNames[iVertexCat]]=0;
300  _pCJetBCTag[_VertexCatNames[iVertexCat]]=0;
301 
302  _pBTagBackgroundValues[_VertexCatNames[iVertexCat]]=0;
303  _pCTagBackgroundValues[_VertexCatNames[iVertexCat]]=0;
304  _pBCTagBackgroundValues[_VertexCatNames[iVertexCat]]=0;
305 
306  }
307 
308 
309 
310 
311  for (unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ){
312 
313  std::string nvname = _VertexCatNames[iVertexCat];
314 
315 
316  pTree->cd( "/" + name() + "/" + "/");
317 
318  if (!pTree->cd( _NumVertexCatDir[iVertexCat]+"/")) {
319  pTree->mkdir( _NumVertexCatDir[iVertexCat]+"/");
320  pTree->cd( _NumVertexCatDir[iVertexCat]+"/");
321  }
322 
323  _pLightJetBTag[nvname] = pHistogramFactory->createHistogram1D( "Numbers of light jets by B-tag NN value. ("+ nvname +")",_numberOfPoints , 0, 1.0 );
324  _pLightJetCTag[nvname] = pHistogramFactory->createHistogram1D( "Numbers of light jets by C-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
325  _pLightJetBCTag[nvname] = pHistogramFactory->createHistogram1D( "Numbers of light jets by BC-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
326  _pBJetBTag[nvname] = pHistogramFactory->createHistogram1D( "Numbers of B jets by B-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
327  _pBJetCTag[nvname] = pHistogramFactory->createHistogram1D( "Numbers of B jets by C-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
328  _pBJetBCTag[nvname] = pHistogramFactory->createHistogram1D( "Numbers of B jets by BC-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
329  _pCJetBTag[nvname] = pHistogramFactory->createHistogram1D( "Numbers of C jets by B-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
330  _pCJetCTag[nvname] = pHistogramFactory->createHistogram1D( "Numbers of C jets by C-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
331  _pCJetBCTag[nvname] = pHistogramFactory->createHistogram1D( "Numbers of C jets by BC-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
332 
333  }
334 
335 }
336 
337 void DSTAIDAPlotProcessor::CreateFlavourTagTuple()
338 {
339 
340  //AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
341  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
342 
343 
344 #ifdef USING_JAIDA
345  //something dosen't work for me with the tuples in RAIDA
346  AIDA::ITupleFactory* pTupleFactory=marlin::AIDAProcessor::tupleFactory( this );
347 
348 
349 
350  pTree->cd( "/" + name());
351 
352  //make the ntuple
353  //this breaks the paradigm of reading these in from the flavour tag collections themselves
354  std::string columnNames="int TrueJetFlavour=-1, int NumberOfVertices=-1, int NumberOfTracksInVertices=-1, float DecayLength = -999., float DecayLengthSignificance= -999., float JointProbRPhi= -999., float JointProbZ= -999., float PTCorrectedMass= -999., float RawMomentum= -999., float SecondaryVertexProbability= -999. ";
355 
356 
357  if (!pTree->cd( "/" + name() + "/" )) {
358  pTree->cd( "/" + name() + "/");
359  if (!pTree->cd("/")) {
360  pTree->mkdir( +"/");
361  pTree->cd( "/");
362  }
363  }
364 
365  if (!pTree->cd( "TupleDir/")) {
366  pTree->mkdir( "TupleDir/");
367  pTree->cd( "TupleDir/");
368  }
369 
370  _pMyTuple=pTupleFactory->create( "FlavourTagInputsTuple","FlavourTagInputsTuple", columnNames);
371 
372 
373 
374 #endif
375 
376 
377 
378 }
379 
380 
381 void DSTAIDAPlotProcessor::CreateFlavourTagInputPlots(LCRunHeader* pRun )
382 {
383 
384 
385  AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
386  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
387  std::vector<std::string> VarNames;
388 
389  VarNames.push_back("NumVertices");
390  VarNames.push_back("JointProbRPhi");
391  VarNames.push_back("JointProbZ");
392  VarNames.push_back("NumTracksInVertices");
393  VarNames.push_back("DecayLength");
394  VarNames.push_back("DecayLengthSignificance");
395  VarNames.push_back("RawMomentum");
396  VarNames.push_back("PTCorrectedMass");
397  VarNames.push_back("SecondaryVertexProbability");
398 
399  for (size_t i = 0;i < VarNames.size();++i)
400  {
401 
402  // If there is no histogram for this name then create one
403  if( _inputsHistogramsBJets[VarNames[i]]==0 )
404  {
405  pTree->cd( "/" + name() + "/");
406  if( !pTree->cd("NNInputs" ) )
407  {
408  pTree->mkdir( "NNInputs" ) ;
409  pTree->cd("NNInputs" ) ;
410  }
411 
412  int numberOfPoints=_numberOfPoints/4;
413  double lowerBin=-1;
414  double higerBin=1;
415 
416 
417  if( VarNames[i]=="NumVertices" )
418  {
419  numberOfPoints=5;
420  lowerBin=0.5;
421  higerBin=5.5;
422  }
423  else if( VarNames[i]=="NumTracksInVertices" )
424  {
425  numberOfPoints=16;
426  lowerBin=-0.5;
427  higerBin=15.5;
428  }
429 
430 
431 
432  else if (VarNames[i]=="DecayLengthSignificance")
433  {
434  numberOfPoints=100;
435  lowerBin=0.;
436  higerBin=100.;
437  }
438  else if (VarNames[i]=="DecayLength")
439  {
440  numberOfPoints=100;
441  lowerBin=0.;
442  higerBin=10.;
443  }
444  else if (VarNames[i]=="JointProbRPhi" || VarNames[i]=="JointProbZ"|| VarNames[i]=="SecondaryVertexProbability")
445  {
446  numberOfPoints=100;
447  lowerBin=0.;
448  higerBin=1.0;
449  }
450  else if (VarNames[i]=="RawMomentum" )
451  {
452  numberOfPoints=100;
453  lowerBin=0.;
454  higerBin=50.;
455  }
456  else if (VarNames[i]=="PTCorrectedMass" )
457  {
458  numberOfPoints=100;
459  lowerBin=0.;
460  higerBin=10.;
461  }
462 
463  pTree->cd( "/" + name() + "/" + "NNInputs");
464  if( !pTree->cd( "bJets" ) )
465  {
466  pTree->mkdir( "bJets" ) ;
467  pTree->cd( "bJets" ) ;
468  }
469 
470  _inputsHistogramsBJets[VarNames[i]]=pHistogramFactory->createHistogram1D( VarNames[i], numberOfPoints, lowerBin, higerBin );
471 
472  pTree->cd( "/" + name() + "/" + "NNInputs");
473  if( !pTree->cd( "cJets" ) )
474  {
475  pTree->mkdir( "cJets" ) ;
476  pTree->cd( "cJets" ) ;
477  }
478 
479  _inputsHistogramsCJets[VarNames[i]]=pHistogramFactory->createHistogram1D( VarNames[i], numberOfPoints, lowerBin, higerBin );
480 
481  pTree->cd( "/" + name() + "/" + "NNInputs");
482  if( !pTree->cd( "udsJets/" ) )
483  {
484  pTree->mkdir( "udsJets/" ) ;
485  pTree->cd( "udsJets/" ) ;
486  }
487 
488  _inputsHistogramsUDSJets[VarNames[i]]=pHistogramFactory->createHistogram1D( VarNames[i], numberOfPoints, lowerBin, higerBin );
489 
490  }//end of histogram creation
491  }
492 
493 
494 }
495 
496 
497 
498 
499 
500 void DSTAIDAPlotProcessor::FillTagPlots( LCEvent* pEvent, unsigned int jetNumber)
501 {
502  // int jetType=FindTrueJetType( pEvent, jetNumber );
503 // if( jetType==0 ) return;
504 
505  //get the jet collection
506  LCCollection* jetCollection=pEvent->getCollection( _JetCollectionName );
507  // TypesafeCollection<lcio::ReconstructedParticle> *jetCollection( pEvent, _JetCollectionName );
508 
509  //get the jet
510  ReconstructedParticle* pJet=dynamic_cast<ReconstructedParticle*>( jetCollection->getElementAt(jetNumber) );
511 
512  //get the PIDHandler fot this jet collection
513  PIDHandler pidh( jetCollection ) ;
514 
515 
516  //get the algorithm id for the MC info
517 
518  int mcid = pidh.getAlgorithmID("MCTruth");
519 
520  //get the Particle id object containing the MC info
521  const ParticleID& mcpid = pidh.getParticleID(pJet,mcid);
522 
523  //get the paramters for the MC info
524  FloatVec mcparams = mcpid.getParameters();
525 
526  //get the true jet flavour
527  int jetType = int(mcparams[pidh.getParameterIndex(mcid,"TrueJetFlavour")]);
528 
529  //get the algorithm id for the FT info
530  int alid = pidh.getAlgorithmID("LCFIFlavourTag");
531 
532  //get the Particle id object containing the FT info
533  const ParticleID& pid = pidh.getParticleID(pJet,alid);
534 
535  //get the paramters for the FT info
536  FloatVec params = pid.getParameters();
537 
538  //get the actual flavour tags
539  double bTag= params[pidh.getParameterIndex(alid,"BTag")];
540  double cTag= params[pidh.getParameterIndex(alid,"CTag")];
541  double cTagBBack= params[pidh.getParameterIndex(alid,"BCTag")];
542 
543  unsigned int NumVertices = int(params[pidh.getParameterIndex(alid,"NumVertices")]);
544  //int CQVtx = FindCQVtx(pEvent, jetNumber);
545  //int BQVtx = FindBQVtx(pEvent, jetNumber);
546  //int trueJetCharge = int(FindTrueJetHadronCharge(pEvent,jetNumber));
547 
548  std::string nvname = _VertexCatNames[ (NumVertices>=N_VERTEX_CATEGORIES) ? (N_VERTEX_CATEGORIES) : (NumVertices)];
549 
550 
551 
552  if( jetType==B_JET ) {
553 
554  if( bTag<=1 && bTag>=0 )
555  {
556  _pBJetBTag[nvname]->fill( bTag );
557  }
558  else
559  {
560  _pBJetBTag[nvname]->fill( -0.005 );
561  }
562 
563  if( cTag<=1 && cTag>=0 )
564  {
565  _pBJetCTag[nvname]->fill( cTag );
566  }
567  else
568  {
569  _pBJetCTag[nvname]->fill( -0.005 );
570  }
571  if( cTagBBack<=1 && cTagBBack>=0 )
572  {
573  _pBJetBCTag[nvname]->fill( cTagBBack );
574  }
575  else
576  {
577  _pBJetBCTag[nvname]->fill( -0.005 );
578  }
579 
580  } else if( jetType==C_JET ) {
581 
582  if( bTag<=1 && bTag>=0 )
583  {
584  _pCJetBTag[nvname]->fill( bTag );
585  }
586  else
587  {
588  _pCJetBTag[nvname]->fill( -0.005 );
589  }
590 
591  if( cTag<=1 && cTag>=0 )
592  {
593  _pCJetCTag[nvname]->fill( cTag );
594  }
595  else
596  {
597  _pCJetCTag[nvname]->fill( -0.005 );
598  }
599 
600  if( cTagBBack<=1 && cTagBBack>=0 )
601  {
602  _pCJetBCTag[nvname]->fill( cTagBBack );
603  }
604  else
605  {
606  _pCJetBCTag[nvname]->fill( -0.005 );
607  }
608  } else {
609  if( bTag<=1 && bTag>=0 )
610  {
611  _pLightJetBTag[nvname]->fill( bTag );
612  }
613  else
614  {
615  _pLightJetBTag[nvname]->fill( -0.005 );
616  }
617 
618  if( cTag<=1 && cTag>=0 )
619  {
620  _pLightJetCTag[nvname]->fill( cTag );
621  }
622  else
623  {
624  _pLightJetCTag[nvname]->fill( -0.005 );
625  }
626 
627  if( cTagBBack<=1 && cTagBBack>=0 )
628  {
629  _pLightJetBCTag[nvname]->fill( cTagBBack );
630  }
631  else
632  {
633  _pLightJetBCTag[nvname]->fill( -0.005 );
634  }
635  }
636 
637 }
638 
639 void DSTAIDAPlotProcessor::FillInputsPlots( LCEvent* pEvent, unsigned int jetNumber )
640 {
641 
642 
643 
644 //jet the jet collection
645  LCCollection* jetCollection=pEvent->getCollection( _JetCollectionName );
646  // TypesafeCollection<lcio::ReconstructedParticle> *jetCollection( pEvent, _JetCollectionName );
647 
648  //get the jet
649  ReconstructedParticle* pJet=dynamic_cast<ReconstructedParticle*>( jetCollection->getElementAt(jetNumber) );
650 
651  //get the PIDHandler fot this jet collection
652  PIDHandler pidh( jetCollection ) ;
653 
654 
655  //get the algorithm id for the MC info
656  int mcid = pidh.getAlgorithmID("MCTruth");
657 
658  //get the Particle id object containing the MC info
659  const ParticleID& mcpid = pidh.getParticleID(pJet,mcid);
660 
661  //get the paramters for the MC info
662  FloatVec mcparams = mcpid.getParameters();
663 
664  //get the true jet flavour
665  int jetType = int(mcparams[pidh.getParameterIndex(mcid,"TrueJetFlavour")]);
666 
667 
668  int alid = pidh.getAlgorithmID("LCFIFlavourTag");
669 
670  //get the Particle id object containing the FT info
671  const ParticleID& pid = pidh.getParticleID(pJet,alid);
672 
673  //get the paramters for the FT info
674  FloatVec params = pid.getParameters();
675 
676  //get the actual flavour tags
677 
678 
679  int NumVertices = int(params[pidh.getParameterIndex(alid,"NumVertices")]);
680 
681 
682 
683  float JointProbRPhi=-999.;
684  float JointProbZ=-999.;
685  int NumTracksInVertices=-999;
686  float DecayLength=-999.;
687  float DecayLengthSignificance=-999.;
688  float RawMomentum=-999.;
689  float PTCorrectedMass=-999.;
690  float SecondaryVertexProbability=-999.;
691 
692  if(NumVertices > 1)
693  {
694  JointProbRPhi= params[pidh.getParameterIndex(alid,"JointProbRPhi")];
695  JointProbZ= params[pidh.getParameterIndex(alid,"JointProbZ")];
696  NumTracksInVertices= int(params[pidh.getParameterIndex(alid,"NumTracksInVertices")]);
697  DecayLength=params[pidh.getParameterIndex(alid,"DecayLength")];
698  DecayLengthSignificance= params[pidh.getParameterIndex(alid,"DecayLengthSignificance")];
699  RawMomentum= params[pidh.getParameterIndex(alid,"RawMomentum")];
700  PTCorrectedMass= params[pidh.getParameterIndex(alid,"PTCorrectedMass")];
701  SecondaryVertexProbability=params[pidh.getParameterIndex(alid,"SecondaryVertexProbability")];
702  }
703 
704  if (_pMyTuple) {
705 
706  _pMyTuple->fill( 0, jetType );
707  _pMyTuple->fill( 1, NumVertices );
708  if(NumTracksInVertices > 1)
709  {
710  _pMyTuple->fill( 2, NumTracksInVertices );
711  _pMyTuple->fill( 3, DecayLength);
712  _pMyTuple->fill( 4, DecayLengthSignificance);
713  _pMyTuple->fill( 5, JointProbRPhi);
714  _pMyTuple->fill( 6, JointProbZ);
715  _pMyTuple->fill( 7, PTCorrectedMass);
716  _pMyTuple->fill( 8, RawMomentum);
717  _pMyTuple->fill( 9, SecondaryVertexProbability);
718  }
719  _pMyTuple->addRow();
720 
721  }//endif _pMyTuple
722 
723 
724  if( jetType==B_JET ) _inputsHistogramsBJets["NumVertices"]->fill(NumVertices);
725  else if( jetType==C_JET ) _inputsHistogramsCJets["NumVertices"]->fill(NumVertices);
726  else _inputsHistogramsUDSJets["NumVertices"]->fill(NumVertices);
727 
728 
729  if(NumVertices > 1)
730  {
731  if( jetType==B_JET ) _inputsHistogramsBJets["NumTracksInVertices"]->fill(NumTracksInVertices);
732  else if( jetType==C_JET ) _inputsHistogramsCJets["NumTracksInVertices"]->fill(NumTracksInVertices);
733  else _inputsHistogramsUDSJets["NumTracksInVertices"]->fill(NumTracksInVertices);
734 
735  if( jetType==B_JET ) _inputsHistogramsBJets["DecayLength"]->fill(DecayLength);
736  else if( jetType==C_JET ) _inputsHistogramsCJets["DecayLength"]->fill(DecayLength);
737  else _inputsHistogramsUDSJets["DecayLength"]->fill(DecayLength);
738 
739  if( jetType==B_JET ) _inputsHistogramsBJets["DecayLengthSignificance"]->fill(DecayLengthSignificance);
740  else if( jetType==C_JET ) _inputsHistogramsCJets["DecayLengthSignificance"]->fill(DecayLengthSignificance);
741  else _inputsHistogramsUDSJets["DecayLengthSignificance"]->fill(DecayLengthSignificance);
742 
743  if( jetType==B_JET ) _inputsHistogramsBJets["JointProbRPhi"]->fill(JointProbRPhi);
744  else if( jetType==C_JET ) _inputsHistogramsCJets["JointProbRPhi"]->fill(JointProbRPhi);
745  else _inputsHistogramsUDSJets["JointProbRPhi"]->fill(JointProbRPhi);
746 
747  if( jetType==B_JET ) _inputsHistogramsBJets["JointProbZ"]->fill(JointProbZ);
748  else if( jetType==C_JET ) _inputsHistogramsCJets["JointProbZ"]->fill(JointProbZ);
749  else _inputsHistogramsUDSJets["JointProbZ"]->fill(JointProbZ);
750 
751  if( jetType==B_JET ) _inputsHistogramsBJets["PTCorrectedMass"]->fill(PTCorrectedMass);
752  else if( jetType==C_JET ) _inputsHistogramsCJets["PTCorrectedMass"]->fill(PTCorrectedMass);
753  else _inputsHistogramsUDSJets["PTCorrectedMass"]->fill(PTCorrectedMass);
754 
755  if( jetType==B_JET ) _inputsHistogramsBJets["RawMomentum"]->fill(RawMomentum);
756  else if( jetType==C_JET ) _inputsHistogramsCJets["RawMomentum"]->fill(RawMomentum);
757  else _inputsHistogramsUDSJets["RawMomentum"]->fill(RawMomentum);
758 
759  if( jetType==B_JET ) _inputsHistogramsBJets["SecondaryVertexProbability"]->fill(SecondaryVertexProbability);
760  else if( jetType==C_JET ) _inputsHistogramsCJets["SecondaryVertexProbability"]->fill(SecondaryVertexProbability);
761  else _inputsHistogramsUDSJets["SecondaryVertexProbability"]->fill(SecondaryVertexProbability);
762 
763 
764  }
765 
766 
767 
768 
769 }
770 
771 
772 
773 void DSTAIDAPlotProcessor::CalculateIntegralAndBackgroundPlots() {
774 
775  AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
776  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
777 
778 
779  for (unsigned int iVertexCat=1; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
780  //sum over the different vertex catagories, this information goes into the "AnyNumberOfVertices" directory
781 
782  _pBJetBTag[_VertexCatNames[0]] -> add(*_pBJetBTag[_VertexCatNames[iVertexCat]]);
783  _pBJetCTag[_VertexCatNames[0]] -> add(*_pBJetCTag[_VertexCatNames[iVertexCat]]);
784  _pBJetBCTag[_VertexCatNames[0]] -> add(*_pBJetBCTag[_VertexCatNames[iVertexCat]]);
785  _pCJetBTag[_VertexCatNames[0]] -> add(*_pCJetBTag[_VertexCatNames[iVertexCat]]);
786  _pCJetCTag[_VertexCatNames[0]] -> add(*_pCJetCTag[_VertexCatNames[iVertexCat]]);
787  _pCJetBCTag[_VertexCatNames[0]] -> add(*_pCJetBCTag[_VertexCatNames[iVertexCat]]);
788  _pLightJetBTag[_VertexCatNames[0]] -> add(*_pLightJetBTag[_VertexCatNames[iVertexCat]]);
789  _pLightJetCTag[_VertexCatNames[0]] -> add(*_pLightJetCTag[_VertexCatNames[iVertexCat]]);
790  _pLightJetBCTag[_VertexCatNames[0]] -> add(*_pLightJetBCTag[_VertexCatNames[iVertexCat]]);
791  }
792 
793  for (unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
794  //add up all the background values
795 
796  pTree->cd( "/" + name() + "/" +_NumVertexCatDir[iVertexCat]);
797 
798  _pBTagBackgroundValues[_VertexCatNames[iVertexCat]] = pHistogramFactory->add("Numbers of non-B jets by B-tag NN value. ("+ _VertexCatNames[iVertexCat]+")",*_pLightJetBTag[_VertexCatNames[iVertexCat]],*_pCJetBTag[_VertexCatNames[iVertexCat]]);
799  _pCTagBackgroundValues[_VertexCatNames[iVertexCat]] = pHistogramFactory->add("Numbers of non-C jets by C-tag NN value. ("+ _VertexCatNames[iVertexCat]+")",*_pLightJetCTag[_VertexCatNames[iVertexCat]],*_pBJetCTag[_VertexCatNames[iVertexCat]]);
800  _pBCTagBackgroundValues[_VertexCatNames[iVertexCat]] = pHistogramFactory->add("Numbers of non-C jets by BC-tag NN value. ("+ _VertexCatNames[iVertexCat]+")",*_pLightJetBCTag[_VertexCatNames[iVertexCat]],*_pBJetBCTag[_VertexCatNames[iVertexCat]]);
801  }
802 
803 
804 }
805 
806 void DSTAIDAPlotProcessor::CalculateEfficiencyPurityPlots()
807 {
808  AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
809  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
810 
811 #ifdef USING_JAIDA
812  AIDA::IDataPointSetFactory* pDataPointSetFactory=marlin::AIDAProcessor::dataPointSetFactory(this);
813 #endif
814 
815 
816  for (unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
817 
818  pTree->cd( "/" + name() + "/" +_NumVertexCatDir[iVertexCat] );
819 
820  std::string nvname = _VertexCatNames[iVertexCat];
821 
822 #ifdef USING_JAIDA
823  AIDA::IDataPointSet* _pBJetBTagEfficiency = CreateEfficiencyPlot( _pBJetBTag[nvname] , pDataPointSetFactory->create("B-Tag efficiency ("+ nvname +")",2));
824  AIDA::IDataPointSet* _pCJetCTagEfficiency = CreateEfficiencyPlot( _pCJetCTag[nvname] , pDataPointSetFactory->create("C-Tag efficiency ("+ nvname +")",2));
825  AIDA::IDataPointSet* _pCJetBCTagEfficiency = CreateEfficiencyPlot( _pCJetBCTag[nvname] , pDataPointSetFactory->create("BC-Tag efficiency ("+ nvname +")",2));
826 #endif
827 
828 
829  _pBJetBTagIntegral[nvname] =
830  CreateIntegralHistogram( _pBJetBTag[nvname],
831  pHistogramFactory->createHistogram1D("B-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
832  _pBJetBTag[nvname]->axis().bins(),_pBJetBTag[nvname]->axis().lowerEdge(),_pBJetBTag[nvname]->axis().upperEdge()));
833 
834  _pCJetBTagIntegral[nvname] =
835  CreateIntegralHistogram( _pCJetBTag[nvname],
836  pHistogramFactory->createHistogram1D("C-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
837  _pCJetBTag[nvname]->axis().bins(),_pCJetBTag[nvname]->axis().lowerEdge(),_pCJetBTag[nvname]->axis().upperEdge()));
838 
839  _pLightJetBTagIntegral[nvname] =
840  CreateIntegralHistogram( _pLightJetBTag[nvname],
841  pHistogramFactory->createHistogram1D("Light-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
842  _pLightJetBTag[nvname]->axis().bins(),_pLightJetBTag[nvname]->axis().lowerEdge(),_pLightJetBTag[nvname]->axis().upperEdge()));
843 
844  _pBJetCTagIntegral[nvname] =
845  CreateIntegralHistogram( _pBJetCTag[nvname],
846  pHistogramFactory->createHistogram1D("B-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
847  _pBJetCTag[nvname]->axis().bins(),_pBJetCTag[nvname]->axis().lowerEdge(),_pBJetCTag[nvname]->axis().upperEdge()));
848 
849  _pCJetCTagIntegral[nvname] =
850  CreateIntegralHistogram( _pCJetCTag[nvname],
851  pHistogramFactory->createHistogram1D("C-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
852  _pCJetCTag[nvname]->axis().bins(),_pCJetCTag[nvname]->axis().lowerEdge(),_pCJetCTag[nvname]->axis().upperEdge()));
853 
854  _pLightJetCTagIntegral[nvname] =
855  CreateIntegralHistogram( _pLightJetCTag[nvname],
856  pHistogramFactory->createHistogram1D("Light-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
857  _pLightJetCTag[nvname]->axis().bins(),_pLightJetCTag[nvname]->axis().lowerEdge(),_pLightJetCTag[nvname]->axis().upperEdge()));
858 
859  _pBJetBCTagIntegral[nvname] =
860  CreateIntegralHistogram( _pBJetBCTag[nvname],
861  pHistogramFactory->createHistogram1D("B-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
862  _pBJetBCTag[nvname]->axis().bins(),_pBJetBCTag[nvname]->axis().lowerEdge(),_pBJetBCTag[nvname]->axis().upperEdge()));
863 
864  _pCJetBCTagIntegral[nvname] =
865  CreateIntegralHistogram( _pCJetBCTag[nvname],
866  pHistogramFactory->createHistogram1D("C-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
867  _pCJetBCTag[nvname]->axis().bins(),_pCJetBCTag[nvname]->axis().lowerEdge(),_pCJetBCTag[nvname]->axis().upperEdge()));
868 
869  _pLightJetBCTagIntegral[nvname] =
870  CreateIntegralHistogram( _pLightJetBCTag[nvname],
871  pHistogramFactory->createHistogram1D("Light-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +")",
872  _pLightJetBCTag[nvname]->axis().bins(),_pLightJetBCTag[nvname]->axis().lowerEdge(),_pLightJetBCTag[nvname]->axis().upperEdge()));
873 
874 
875 
876 #ifdef USING_JAIDA
877  AIDA::IDataPointSet* _pBJetBTagPurity = CreatePurityPlot( _pBJetBTag[nvname], _pBTagBackgroundValues[nvname] , pDataPointSetFactory->create("B-Jet purity for B-Tag ("+ nvname +")",2));
878  AIDA::IDataPointSet* _pCJetCTagPurity = CreatePurityPlot( _pCJetCTag[nvname], _pCTagBackgroundValues[nvname] , pDataPointSetFactory->create("C-Jet purity for C-Tag ("+ nvname +")",2));
879  AIDA::IDataPointSet* _pCJetBCTagPurity = CreatePurityPlot( _pCJetBCTag[nvname], _pBJetBCTag[nvname], pDataPointSetFactory->create("C-Jet purity for BC-Tag ("+ nvname +")",2));
880 
881  AIDA::IDataPointSet* _pCJetBTagLeakage = CreateLeakageRatePlot( _pCJetBTag[nvname], pDataPointSetFactory->create("C-Jets: Leakage Rate into B-Tag Sample ("+ nvname +")",2));
882  AIDA::IDataPointSet* _pLightJetBTagLeakage = CreateLeakageRatePlot( _pLightJetBTag[nvname], pDataPointSetFactory->create("Light-Jets: Leakage Rate into B-Tag Sample ("+ nvname +")",2));
883  AIDA::IDataPointSet* _pBJetCTagLeakage = CreateLeakageRatePlot( _pBJetCTag[nvname], pDataPointSetFactory->create("B-Jets: Leakage Rate into C-Tag Sample ("+ nvname +")",2));
884  AIDA::IDataPointSet* _pLightJetCTagLeakage = CreateLeakageRatePlot( _pLightJetCTag[nvname], pDataPointSetFactory->create("Light-Jets: Leakage Rate into C-Tag Sample ("+ nvname +")",2));
885  AIDA::IDataPointSet* _pBJetBCTagLeakage = CreateLeakageRatePlot( _pBJetBCTag[nvname], pDataPointSetFactory->create("B-Jets: Leakage Rate into BC-Tag Sample ("+ nvname +")",2));
886  AIDA::IDataPointSet* _pLightJetBCTagLeakage = CreateLeakageRatePlot( _pLightJetBCTag[nvname], pDataPointSetFactory->create("Light-Jets: Leakage Rate into BC-Tag Sample ("+ nvname +")",2));
887  AIDA::IDataPointSet* _pNonBJetBTagLeakage = CreateLeakageRatePlot( _pBTagBackgroundValues[nvname], pDataPointSetFactory->create("C-Jets: Leakage Rate into B-Tag Sample ("+ nvname +")",2));
888  AIDA::IDataPointSet* _pNonCJetCTagLeakage = CreateLeakageRatePlot( _pCTagBackgroundValues[nvname], pDataPointSetFactory->create("Light-Jets: Leakage Rate into B-Tag Sample ("+ nvname +")",2));
889  AIDA::IDataPointSet* _pNonCJetBCTagLeakage = CreateLeakageRatePlot( _pBCTagBackgroundValues[nvname], pDataPointSetFactory->create("B-Jets: Leakage Rate into C-Tag Sample ("+ nvname +")",2));
890 
891  CreateXYPlot(_pBJetBTagEfficiency, _pBJetBTagPurity, pDataPointSetFactory->create("Purity-Efficiency for B-Tag ("+ nvname +")",2), 1, 1);
892  CreateXYPlot(_pCJetCTagEfficiency, _pCJetCTagPurity, pDataPointSetFactory->create("Purity-Efficiency for C-Tag ("+ nvname +")",2), 1, 1);
893  CreateXYPlot(_pCJetBCTagEfficiency, _pCJetBCTagPurity, pDataPointSetFactory->create("Purity-Efficiency for BC-Tag ("+ nvname +")",2), 1, 1);
894 
895  CreateXYPlot(_pBJetBTagEfficiency, _pNonBJetBTagLeakage, pDataPointSetFactory->create("Leakage Rate-Efficiency for B-Tag ("+ nvname +")",2), 1, 1);
896  CreateXYPlot(_pCJetCTagEfficiency, _pNonCJetCTagLeakage, pDataPointSetFactory->create("Leakage Rate-Efficiency for C-Tag ("+ nvname +")",2), 1, 1);
897  CreateXYPlot(_pCJetBCTagEfficiency,_pNonCJetBCTagLeakage, pDataPointSetFactory->create("Leakage Rate-Efficiency for BC-Tag ("+ nvname +")",2), 1, 1);
898 
899  CreateXYPlot(_pBJetBTagEfficiency, _pCJetBTagLeakage, pDataPointSetFactory->create("Leakage Rate (C into B) vs Efficiency for B-Tag ("+ nvname +")",2), 1, 1);
900  CreateXYPlot(_pBJetBTagEfficiency, _pLightJetBTagLeakage, pDataPointSetFactory->create("Leakage Rate (Light into B) Efficiency for B-Tag ("+ nvname +")",2), 1, 1);
901  CreateXYPlot(_pCJetCTagEfficiency, _pBJetCTagLeakage, pDataPointSetFactory->create("Leakage Rate (B into C) vs Efficiency for C-Tag ("+ nvname +")",2), 1, 1);
902  CreateXYPlot(_pCJetCTagEfficiency, _pLightJetCTagLeakage, pDataPointSetFactory->create("Leakage Rate (Light into C) Efficiency for C-Tag ("+ nvname +")",2), 1, 1);
903  CreateXYPlot(_pCJetBCTagEfficiency, _pBJetBCTagLeakage, pDataPointSetFactory->create("Leakage Rate (B into BC) vs Efficiency for BC-Tag ("+ nvname +")",2), 1, 1);
904  CreateXYPlot(_pCJetCTagEfficiency, _pLightJetBCTagLeakage, pDataPointSetFactory->create("Leakage Rate (Light into BC) Efficiency for BC-Tag ("+ nvname +")",2), 1, 1);
905 #endif
906 
907  }//end of loop of different number of vertices
908 
909 
910 
911 }
912 
913 
914 AIDA::IHistogram1D* DSTAIDAPlotProcessor::CreateIntegralHistogram(const AIDA::IHistogram1D* pNN, AIDA::IHistogram1D* pIntegral)
915 {
916  //the errors on these entries are wrong...
917 
918  const int numberOfBins=pNN->axis().bins();
919 
920  double integral=pNN->binHeight(AIDA::IAxis::OVERFLOW_BIN);
921  pIntegral->fill(pNN->axis().binLowerEdge(AIDA::IAxis::OVERFLOW_BIN)+pNN->axis().binWidth(numberOfBins-1)/2.,integral);
922 
923  for( int binNumber=numberOfBins-1; binNumber>=0; --binNumber )
924  {
925  integral+= pNN->binHeight( binNumber );
926  pIntegral->fill( pNN->axis().binLowerEdge(binNumber)+pNN->axis().binWidth(binNumber)/2.,integral);
927  }
928 
929  integral+= pNN->binHeight(AIDA::IAxis::UNDERFLOW_BIN);
930  pIntegral->fill(pNN->axis().binUpperEdge(AIDA::IAxis::UNDERFLOW_BIN)-pNN->axis().binWidth(0)/2.,integral);
931 
932  return pIntegral;
933 }
934 
935 
936 
937 #ifdef USING_JAIDA
938 AIDA::IDataPointSet* DSTAIDAPlotProcessor::CreateEfficiencyPlot(const AIDA::IHistogram1D* pSignal, AIDA::IDataPointSet* pDataPointSet)
939 {
940  //makes an effiency plot from one histogram - assuming the last bin contains all the events (eg NN distribution)
941  double totalSignal=pSignal->sumBinHeights();
942  double signalPassedCut=0;
943 
944  const int numberOfBins=pSignal->axis().bins();
945  int iPoint=0;
946 
947  for( int binNumber=numberOfBins-1; binNumber>=0; --binNumber, iPoint++ )
948  {
949  signalPassedCut+=pSignal->binHeight( binNumber );
950 
951  double efficiency = signalPassedCut/totalSignal;
952 
953  double efficiencyError = efficiency * (1. - efficiency) / totalSignal;
954  if (efficiencyError>0) efficiencyError = sqrt(efficiencyError);
955 
956 
957  pDataPointSet->addPoint();
958  pDataPointSet->point(iPoint)->coordinate(0)->setValue(pSignal->axis().binLowerEdge(binNumber)+pSignal->axis().binWidth(binNumber)/2.);
959  pDataPointSet->point(iPoint)->coordinate(1)->setValue( efficiency );
960  pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(efficiencyError);
961  pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(efficiencyError);
962  }
963 
964  return pDataPointSet;
965 }
966 
967 AIDA::IDataPointSet* DSTAIDAPlotProcessor::CreatePurityPlot(const AIDA::IHistogram1D* pSignal, const AIDA::IHistogram1D* pBackground, AIDA::IDataPointSet* pDataPointSet)
968 {
969  const int numberOfBins=pSignal->axis().bins();
970  int iPoint=0;
971 
972  double signalPassedCut=0;
973  double backgroundPassedCut=0;
974 
975  for (int binNumber = numberOfBins-1; binNumber >= 0 ; --binNumber, iPoint++ )
976  {
977 
978  signalPassedCut+=pSignal->binHeight( binNumber );
979  backgroundPassedCut+=pBackground->binHeight( binNumber );
980 
981  double purity = signalPassedCut/(signalPassedCut+backgroundPassedCut);
982  double purityError = purity * (1. - purity) / (signalPassedCut+backgroundPassedCut);
983  if (purityError>0) purityError = sqrt(purityError);
984 
985 
986  pDataPointSet->addPoint();
987  pDataPointSet->point(iPoint)->coordinate(0)->setValue(pSignal->axis().binLowerEdge(binNumber)+pSignal->axis().binWidth(binNumber)/2.);
988  pDataPointSet->point(iPoint)->coordinate(1)->setValue(purity);
989  pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(purityError);
990  pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(purityError);
991 
992 
993  }
994 
995  return pDataPointSet;
996 }
997 AIDA::IDataPointSet* DSTAIDAPlotProcessor::CreateLeakageRatePlot(const AIDA::IHistogram1D* pBackground, AIDA::IDataPointSet* pDataPointSet)
998 {
999 
1000  double totalBackground = pBackground->sumBinHeights();
1001  double backgroundPassedCut=0;
1002 
1003  const int numberOfBins=pBackground->axis().bins();
1004  int iPoint=0;
1005 
1006  for( int binNumber=numberOfBins-1; binNumber>=0; --binNumber , iPoint++ )
1007  {
1008 
1009  backgroundPassedCut+=pBackground->binHeight( binNumber );
1010 
1011  double leakageRate = backgroundPassedCut/totalBackground;
1012 
1013  double leakageRateError = leakageRate * (1. - leakageRate) / totalBackground;
1014  if (leakageRateError>0) leakageRateError = sqrt(leakageRateError);
1015 
1016 
1017  pDataPointSet->addPoint();
1018  pDataPointSet->point(iPoint)->coordinate(0)->setValue(pBackground->axis().binLowerEdge(binNumber)+pBackground->axis().binWidth(binNumber)/2.);
1019  pDataPointSet->point(iPoint)->coordinate(1)->setValue(leakageRate);
1020  pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(leakageRateError);
1021  pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(leakageRateError);
1022  }
1023 
1024  return pDataPointSet;
1025 }
1026 
1027 
1028 AIDA::IDataPointSet* DSTAIDAPlotProcessor::CreateXYPlot(const AIDA::IDataPointSet* pDataPointSet0, const AIDA::IDataPointSet* pDataPointSet1, AIDA::IDataPointSet* xyPointSet, const int dim0, const int dim1 )
1029 {
1030 
1031  //need to do some comparision here
1032  if (pDataPointSet0->size() == pDataPointSet1->size()) {
1033 
1034  for (int iPoint = 0 ; iPoint != pDataPointSet1->size(); iPoint++)
1035  {
1036  xyPointSet->addPoint();
1037  xyPointSet->point(iPoint)->coordinate(0)->setValue(pDataPointSet0->point(iPoint)->coordinate(dim0)->value());
1038  xyPointSet->point(iPoint)->coordinate(1)->setValue(pDataPointSet1->point(iPoint)->coordinate(dim1)->value());
1039  xyPointSet->point(iPoint)->coordinate(0)->setErrorPlus(pDataPointSet0->point(iPoint)->coordinate(dim0)->errorPlus());
1040  xyPointSet->point(iPoint)->coordinate(1)->setErrorPlus(pDataPointSet1->point(iPoint)->coordinate(dim1)->errorPlus());
1041  xyPointSet->point(iPoint)->coordinate(0)->setErrorMinus(pDataPointSet0->point(iPoint)->coordinate(dim0)->errorMinus());
1042  xyPointSet->point(iPoint)->coordinate(1)->setErrorMinus(pDataPointSet1->point(iPoint)->coordinate(dim1)->errorMinus());
1043  }
1044  } else {
1045 
1046  //some error message here!
1047  }
1048  return xyPointSet;
1049 }
1050 
1051 #endif
1052 
1053 
1054 #endif // end of "#ifndef MARLIN_USE_AIDA"
1055 
DSTAIDAPlotProcessor::_displayCollectionNames
void _displayCollectionNames(lcio::LCEvent *pEvent)
DSTAIDAPlotProcessor::_passesEventCuts
bool _passesEventCuts(lcio::LCEvent *pEvent)
vertex_lcfi::util::Vector3
Multi-Perpose 3 Vector.
Definition: include/vertex_lcfi/util/inc/vector3.h:33
DSTAIDAPlotProcessor::_passesJetCuts
bool _passesJetCuts(lcio::ReconstructedParticle *pJet)
Generated on Fri Dec 2 2016 12:36:49 for LCFIVertex by   doxygen 1.8.6