LCFIVertex  0.7.2
LCFIAIDAPlotProcessor.cc
1 // First of all, make sure that AIDA was enabled by setting the environment
2 // variable MARLIN_USE_AIDA when Marlin was compiled. The makefile will then
3 // have done the setup and defined this macro.
4 
5 #ifndef MARLIN_USE_AIDA
6 
7 #warning "--------------------------------------------------------------------------------"
8 #warning "- LCFIAIDAPlotProcessor requires MARLIN_USE_AIDA to be defined. Did you enable -"
9 #warning "- AIDA when compiling Marlin? LCFIAIDAPlotProcessor will not be compiled. -"
10 #warning "--------------------------------------------------------------------------------"
11 
12 // Can't do anything else.
13 #else
14 
15 #include "LCFIAIDAPlotProcessor.h"
16 
17 #ifdef __APPLE__
18 #define uint uint32_t
19 #endif
20 
21 // standard library includes
22 #include <iostream>
23 #include <iomanip>
24 #include <string>
25 #include <vector>
26 #include <set>
27 //#include <math.h>
28 #include <cmath>
29 #include <cstdlib>
30 #include <algorithm>
31 
32 // LCIO includes...
33 #include "EVENT/LCCollection.h"
34 #include "EVENT/LCParameters.h"
35 #include "EVENT/LCIntVec.h"
36 #include "EVENT/LCFloatVec.h"
37 #include "EVENT/Vertex.h"
38 #include "EVENT/MCParticle.h"
39 #include "EVENT/Track.h"
40 
41 
42 // Marlin includes
43 #include <marlin/Exceptions.h>
44 
45 // AIDA includes...
46 #include <marlin/AIDAProcessor.h>
47 #include <AIDA/IMeasurement.h>
48 #include <AIDA/IHistogramFactory.h>
49 #include <AIDA/IHistogram1D.h>
50 #include <AIDA/IHistogram2D.h>
51 #include <AIDA/IAxis.h>
52 #include <AIDA/ITree.h>
53 #include <AIDA/ITupleFactory.h>
54 #include <AIDA/ICloud2D.h>
55 #include <UTIL/LCRelationNavigator.h>
56 
57 //change USING_RAIDA to USING_JAIDA if you are using JAIDA/AIDAJNI - you will obtain more functionality!
58 #define USING_RAIDA
59 
60 #ifdef USING_RAIDA
61 #warning "USING_RAIDA defined"
62 #else
63 #define USING_JAIDA
64 #warning "USING_JAIDA defined"
65 #endif
66 
67 
68 #ifdef USING_JAIDA//Data point sets aren't implemented in RAIDA - which is a shame as they have functionality not given by histograms
69 //such as the facility to set the error
70 #include <AIDA/IDataPointSet.h>
71 #include <AIDA/IDataPointSetFactory.h>
72 #include <AIDA/IDataPoint.h>
73 #endif
74 
75 #include "TypesafeCollection.h"
76 
77 // There needs to be at least one instantiation for the base constructor to register the processor with
78 // the Marlin processor manager. This is it.
79 LCFIAIDAPlotProcessor aLCFIAIDAPlotProcessor;
80 
81 LCFIAIDAPlotProcessor::LCFIAIDAPlotProcessor() : marlin::Processor( "LCFIAIDAPlotProcessor" )
82 {
83 
84  _description="Creates an AIDA plot of the LCFIVertex tagging efficiency-purity values and various other things. Make sure that MarlinAIDAProcessor is run before this.";
85 
86  registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE,
87  "JetCollectionName" ,
88  "Name of the collection of ReconstructedParticles that is the jet" ,
89  _JetCollectionName ,
90  std::string("FTSelectedJets") );
91 
92  std::vector<std::string> FlavourTagCollectionNamesDefault;
93  FlavourTagCollectionNamesDefault.push_back("FlavourTag");
94  registerProcessorParameter("FlavourTagCollections" ,
95  "Names of the LCFloatVec Collections that contain the flavour tags (one purity efficiency plot per tag) (in same order as jet collection)" ,
96  _FlavourTagCollectionNames,
97  FlavourTagCollectionNamesDefault) ;
98 
99  registerOptionalParameter( "TrueJetFlavourCollection" ,
100  "Name of the LCFloatVec collection containing the true flavour of the jets (same order as jets)" ,
101  _TrueJetFlavourColName ,
102  std::string("TrueJetFlavour") ) ;
103 
104  registerInputCollection( lcio::LCIO::MCPARTICLE,
105  "MCParticleCollection" ,
106  "Name of the collection that holds all MC particles. " ,
107  _MCParticleColName ,
108  std::string("MCParticle") ) ;
109 
110  registerOptionalParameter("VertexCollection",
111  "Name of the collection that holds the Vertices",
112  _VertexColName,
113  std::string("ZVRESVertices") ) ;
114 
115  registerOptionalParameter("CVertexChargeCollection",
116  "Name of collection containing the vertex charge of the jets, assuming they are C-jets",
117  _CVertexChargeCollection,
118  std::string("CCharge") );
119 
120  registerOptionalParameter( "BVertexChargeCollection",
121  "Name of collection containing the vertex charge of the jets, assuming they are B-jets",
122  _BVertexChargeCollection,
123  std::string("BCharge") ) ;
124 
125  registerOptionalParameter("TrueTracksToMCPCollection",
126  "Name of collection linking the tracks and the Monte Carlo Particles",
127  _TrueTracksToMCPCollection,
128  std::string("LDCTracksMCP") ) ;
129 
130  registerOptionalParameter( "ZVRESDecayChainCollection" ,
131  "Name of the ZVRES DecayChain collection" ,
132  _ZVRESDecayChainCollection ,
133  std::string("ZVRESDecayChains") );
134 
135  registerOptionalParameter( "ZVRESSelectedJetsCollection" ,
136  "Name of the ZVRES Selected Jets collection" ,
137  _ZVRESSelectedJetsCollection ,
138  std::string("ZVRESSelectedJets") );
139 
140  registerOptionalParameter("ZVRESDecayChainTrackCollection" ,
141  "Name of the ZVRES Decay Chain Tracks Collection" ,
142  _ZVRESDecayChainRPTracksCollection ,
143  std::string("ZVRESDecayChainRPTracks") );
144 
145  FlavourTagCollectionNamesDefault.clear();
146  FlavourTagCollectionNamesDefault.push_back("FlavourTagInputs");
147  registerProcessorParameter("TagInputsCollections" ,
148  "Names of the LCFloatVec Collections that contain the flavour tag inputs (in same order as jet collection)" ,
149  _FlavourTagInputsCollectionNames,
150  FlavourTagCollectionNamesDefault) ;
151 
152  registerOptionalParameter( "CosThetaJetMax",
153  "Cut determining the maximum cos(theta) of the jet. Default: |cos(theta)|<0.9" ,
154  _CosThetaJetMax,
155  double(0.9)) ;
156 
157  registerOptionalParameter( "CosThetaJetMin",
158  "Cut determining the minimum cos(theta) of the jet. Default: no lower cut." ,
159  _CosThetaJetMin,
160  double(0.0)) ;
161 
162  registerOptionalParameter("PJetMax",
163  "Cut determining the maximum momentum of the jet. Default: 10000 GeV/c" ,
164  _PJetMax,
165  double(10000.)) ;
166 
167  registerOptionalParameter( "PJetMin",
168  "Cut determining the minimum momentum of the jet. Default: no lower cut." ,
169  _PJetMin,
170  double(0.0)) ;
171 
172  registerOptionalParameter( "PrintTrackVertexOutput",
173  "Set true if you want a print-out of the track-vertex association purity",
174  _PrintTrackVertexOutput,
175  bool(false));
176 
177  registerOptionalParameter( "MakePurityEfficiencyPlots",
178  "Set true if you want to make the purity-efficiency plots, and leakage rates plots for the various flavour tags",
179  _MakePurityEfficiencyPlots,
180  bool(true));
181 
182  registerOptionalParameter( "PrintPurityEfficiencyValues",
183  "Set true if you want a print-out of the purity-efficiency for the various flavour tags",
184  _PrintPurityEfficiencyValues,
185  bool(true));
186 
187  registerOptionalParameter( "MakeAdditionalPlots",
188  "Set true if you want to make all other plots (i.e. the non purity-efficiency plots) and the other functionality provided by LCFIAIDAPlotProcessor",
189  _MakeAdditionalPlots,
190  bool(false));
191 
192  registerOptionalParameter( "PurityEfficiencyOutputFile" ,
193  "Output filename for the Purity-Efficiency values. Only used if PrintPurityEfficiencyValues parameter is true. If left blank, output will be directed to standard out.",
194  _PurityEfficiencyOutputFile,
195  std::string("PurityEfficiencyOutput.txt") ) ;
196 
197  registerOptionalParameter( "TrackVertexOutputFile" ,
198  "Output filename for the table of the Track-Vertex association. Only used if PrintTrackVertexOutput parameter is true. If left blank, output will be directed to standard out.",
199  _TrackVertexOutputFile,
200  std::string("TrackVertexOutput.txt") ) ;
201 
202  registerOptionalParameter( "MakeTuple",
203  "Set true to make a tuple of the flavour tag input variables. Default is false (only works with jaida).",
204  _MakeTuple,
205  bool(false));
206 
207  registerOptionalParameter( "CTagNNCut",
208  "Cut determining the Neural Net cut used to select C-Jets",
209  _CTagNNCut,
210  double(0.7));
211 
212  registerOptionalParameter( "BTagNNCut",
213  "Cut determining the Neural Net cut used to select B-Jets",
214  _BTagNNCut,
215  double(0.7));
216 
217  registerOptionalParameter( "UseFlavourTagCollectionForVertexCharge",
218  "Integer parameter determing which FlavourTag Collection to use the determine C-Jets and B-Jets in Vertex Charge Plots",
219  _iVertexChargeTagCollection,
220  int(0));
221 
222 }
223 
224 LCFIAIDAPlotProcessor::~LCFIAIDAPlotProcessor()
225 {
226 }
227 
228 void LCFIAIDAPlotProcessor::init()
229 {
230 
231  if ((_iVertexChargeTagCollection >= int(_FlavourTagCollectionNames.size()) || _iVertexChargeTagCollection < 0) && _FlavourTagCollectionNames.size()!=0) {
232  std::cerr << " In " << __FILE__ << "(" << __LINE__ << "): Invalid parameter for UseFlavourTagCollectionForVertexCharge. Setting to 0." << std::endl;
233  _myVertexChargeTagCollection = 0;
234  } else if (_FlavourTagCollectionNames.size()==0) {
235  _myVertexChargeTagCollection = 0;
236  } else {
237  _myVertexChargeTagCollection = unsigned(_iVertexChargeTagCollection);
238  }
239 
240  _ZoomedVarNames.push_back("D0Significance1");
241  _ZoomedVarNames.push_back("D0Significance2");
242  _ZoomedVarNames.push_back("Z0Significance1");
243  _ZoomedVarNames.push_back("Z0Significance2");
244 
245  _VertexCatNames.resize(N_VERTEX_CATEGORIES+1);
246  _VertexCatNames[0]="AnyNumberOfVertices";
247  _VertexCatNames[1]="OneVertex";
248  _VertexCatNames[2]="TwoVertices";
249  _VertexCatNames[3]="ThreeOrMoreVertices";
250 
251 
252  _NumVertexCatDir.resize(N_VERTEX_CATEGORIES+1);
253  _NumVertexCatDir[1]="OneVertex";
254  _NumVertexCatDir[2]="TwoVertices";
255  _NumVertexCatDir[3]="ThreeOrMoreVertices";
256  _NumVertexCatDir[0]="AnyNumberOfVertices";
257 
258 
259  _numberOfPoints=100;
260 
261 
262 
263  AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
264  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
265 
266  if( pHistogramFactory!=0 )
267  {
268  if (!(pTree->cd( "/" + name() + "/"))) {
269  pTree->mkdir( "/" + name() + "/" );
270  pTree->cd( "/" + name() + "/");
271  }
272 
273 
274  if (_MakePurityEfficiencyPlots) CreateTagPlots();
275  if (_MakeAdditionalPlots) CreateAdditionalPlots();
276  if (_MakeTuple) CreateFlavourTagTuple();
277 
278  } else {
279 
280  std::cerr << "### " << __FILE__ << "(" << __LINE__ << "): Unable to get the histogram factory! No histograms will be made."<< std::endl;
281  }
282 
283  _lastRunHeaderProcessed=-1;
284  _suppressOutputForRun=-1;
285 
286 
287  InternalVectorInitialisation();
288 
289 }
290 
291 void LCFIAIDAPlotProcessor::processRunHeader( LCRunHeader* pRun )
292 {
293 
294  // Marlin doesn't necessarily process the run header, e.g. if you use the "SkipNEvents"
295  // parameter in the steering file. The flavour tag variable/tag value names are held in
296  // the run header though, so this processor has to have access to it. Set this variable
297  // so that "processEvent" can tell if "processRunHeader" has been called for the run
298  // it's in.
299  _lastRunHeaderProcessed=pRun->getRunNumber();
300 
301  InitialiseFlavourTagInputs(pRun);
302 
303 
304  //
305  // Perform a check to see if the variable names we need are here
306  //
307  for (unsigned int iTag=0; iTag < _FlavourTagCollectionNames.size(); ++iTag) // Loop over the different tag collection names given in the steering
308  {
309  std::vector<std::string> VarNames;
310  (pRun->parameters()).getStringVals(_FlavourTagCollectionNames[iTag],VarNames);
311 
312 
313  //Fill a map so that we can get the array index from just the string
314  std::set<std::string> AvailableNames;
315  std::map<std::string,unsigned int> IndexOf;
316 
317  for (size_t i = 0;i < VarNames.size();++i)
318  {
319  AvailableNames.insert(VarNames[i]);
320  IndexOf[VarNames[i]] = i;
321  }
322 
323  //Add the index to the list
324  _IndexOfForEachTag.push_back(IndexOf);
325 
326  //Check the required information is in the LCFloatVec
327  std::set<std::string> RequiredNames;
328  RequiredNames.insert("BTag");
329  RequiredNames.insert("CTag");
330  RequiredNames.insert("BCTag");
331 
332  if (!std::includes(AvailableNames.begin(),AvailableNames.end(),RequiredNames.begin(),RequiredNames.end()))
333  {
334  std::cerr << __FILE__ << "(" << __LINE__ << "): The collection \"" << _FlavourTagCollectionNames[iTag]
335  << "\" (if it exists) does not contain the tag values required by " << type() << "." << std::endl;
336  std::cerr << __FILE__ << "(" << __LINE__ << "): The collection \"" << _FlavourTagCollectionNames[iTag]
337  << "\" (if it exists) does not contain the tag values required by " << type() << "." << std::endl;
338  }
339  }
340 
341  if (_MakeTuple) CreateFlavourTagInputPlots(pRun);
342 
343 }
344 
345 void LCFIAIDAPlotProcessor::CreateFlavourTagTuple()
346 {
347 
348  //AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
349  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
350 
351 
352 #ifdef USING_JAIDA
353  //something dosen't work for me with the tuples in RAIDA
354  AIDA::ITupleFactory* pTupleFactory=marlin::AIDAProcessor::tupleFactory( this );
355 
356 
357  if (_MakeTuple) {
358  pTree->cd( "/" + name());
359 
360  //make the ntuple
361  //this breaks the paradigm of reading these in from the flavour tag collections themselves
362  std::string columnNames="int TrueJetFlavour=-1, int NumberOfVertices=-1, int NumberOfTracksInVertices=-1, float D0Significance1 = -999., float D0Significance2 = -999., float DecayLength = -999., float DecayLength_SeedToIP= -999., float DecayLengthSignificance= -999., float JointProbRPhi= -999., float JointProbZ= -999., float Momentum1= -999.,float Momentum2= -999., float PTCorrectedMass= -999., float RawMomentum= -999., float SecondaryVertexProbability= -999., float Z0Significance1= -999., float Z0Significance2= -999., int BQVtx=-10, int CQVtx=-10";
363 
364 
365  if (!pTree->cd( "/" + name() + "/" +_FlavourTagCollectionNames[_myVertexChargeTagCollection] + "/")) {
366  pTree->cd( "/" + name() + "/");
367  if (!pTree->cd( _FlavourTagCollectionNames[_myVertexChargeTagCollection])) {
368  pTree->mkdir( _FlavourTagCollectionNames[_myVertexChargeTagCollection] + "/");
369  pTree->cd( _FlavourTagCollectionNames[_myVertexChargeTagCollection]+ "/");
370  }
371  }
372 
373  if (!pTree->cd( "TupleDir/")) {
374  pTree->mkdir( "TupleDir/");
375  pTree->cd( "TupleDir/");
376  }
377 
378  _pMyTuple=pTupleFactory->create( "FlavourTagInputsTuple","FlavourTagInputsTuple", columnNames);
379  }
380 #endif
381 
382 
383 }
384 
385 void LCFIAIDAPlotProcessor::InitialiseFlavourTagInputs(LCRunHeader* pRun )
386 {
387  for (unsigned int iInputCollection=0; iInputCollection < _FlavourTagInputsCollectionNames.size(); ++iInputCollection)
388  {
389 
390  std::vector<std::string> VarNames;
391  (pRun->parameters()).getStringVals(_FlavourTagInputsCollectionNames[iInputCollection],VarNames);
392  //Fill the map relating names and indexes
393  std::map<std::string,unsigned int> IndexOf;
394  for (size_t i = 0;i < VarNames.size();++i)
395  {
396  IndexOf[VarNames[i]] = i;
397  }
398 
399  _InputsIndex.push_back(IndexOf);
400  }
401 
402  std::vector<std::string> trueJetFlavourVarNames;
403  (pRun->parameters()).getStringVals(_TrueJetFlavourColName,trueJetFlavourVarNames);
404 
405  for (size_t i = 0;i < trueJetFlavourVarNames.size();++i)
406  {
407  _FlavourIndex[trueJetFlavourVarNames[i]] = i;
408  }
409 
410 }
411 
412 
413 
414 void LCFIAIDAPlotProcessor::CreateFlavourTagInputPlots(LCRunHeader* pRun )
415 {
416 
417  _inputsHistogramsBJets.resize( _FlavourTagInputsCollectionNames.size() );
418  _inputsHistogramsCJets.resize( _FlavourTagInputsCollectionNames.size() );
419  _inputsHistogramsUDSJets.resize( _FlavourTagInputsCollectionNames.size() );
420 
421  _zoomedInputsHistogramsBJets.resize( _FlavourTagInputsCollectionNames.size() );
422  _zoomedInputsHistogramsCJets.resize( _FlavourTagInputsCollectionNames.size() );
423  _zoomedInputsHistogramsUDSJets.resize( _FlavourTagInputsCollectionNames.size() );
424 
425  AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
426  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
427 
428  for (unsigned int iInputCollection=0; iInputCollection < _FlavourTagInputsCollectionNames.size(); ++iInputCollection)
429  {
430 
431  std::vector<std::string> VarNames;
432  (pRun->parameters()).getStringVals(_FlavourTagInputsCollectionNames[iInputCollection],VarNames);
433 
434  //Fill the map relating names and indexes
435  for (size_t i = 0;i < VarNames.size();++i)
436  {
437 
438  // If there is no histogram for this name then create one
439  if( _inputsHistogramsBJets[iInputCollection][VarNames[i]]==0 )
440  {
441  pTree->cd( "/" + name() + "/");
442  if( !pTree->cd(_FlavourTagInputsCollectionNames[iInputCollection] ) )
443  {
444  pTree->mkdir( _FlavourTagInputsCollectionNames[iInputCollection] ) ;
445  pTree->cd(_FlavourTagInputsCollectionNames[iInputCollection] ) ;
446  }
447 
448  int numberOfPoints=_numberOfPoints/4;
449  double lowerBin=-1;
450  double higerBin=1;
451 
452  //binning variables: if the name is not listed here it will use the default above
453  if( VarNames[i]=="BQVtx" || VarNames[i]=="CQVtx" )
454  {
455  numberOfPoints=9;
456  lowerBin=-4.5;
457  higerBin=4.5;
458  }
459  else if( VarNames[i]=="NumVertices" )
460  {
461  numberOfPoints=5;
462  lowerBin=0.5;
463  higerBin=5.5;
464  }
465  else if( VarNames[i]=="NumTracksInVertices" )
466  {
467  numberOfPoints=16;
468  lowerBin=-0.5;
469  higerBin=15.5;
470  }
471  else if ( VarNames[i]=="D0Significance1" || VarNames[i]=="D0Significance2" )
472  {
473  numberOfPoints=120;
474  lowerBin=-20.;
475  higerBin=100.;
476  }
477  else if ( VarNames[i]=="Z0Significance1" || VarNames[i]=="Z0Significance2")
478  {
479  numberOfPoints=100;
480  lowerBin=-50.;
481  higerBin=50.;
482  }
483  else if (VarNames[i]=="DecayLengthSignificance")
484  {
485  numberOfPoints=100;
486  lowerBin=0.;
487  higerBin=100.;
488  }
489  else if (VarNames[i]=="DecayLength" || VarNames[i]=="DecayLength(SeedToIP)" )
490  {
491  numberOfPoints=100;
492  lowerBin=0.;
493  higerBin=10.;
494  }
495  else if (VarNames[i]=="JointProbRPhi" || VarNames[i]=="JointProbZ"|| VarNames[i]=="SecondaryVertexProbability")
496  {
497  numberOfPoints=100;
498  lowerBin=0.;
499  higerBin=1.0;
500  }
501  else if (VarNames[i]=="Momentum1" || VarNames[i]=="Momentum2" || VarNames[i]=="RawMomentum" )
502  {
503  numberOfPoints=100;
504  lowerBin=0.;
505  higerBin=50.;
506  }
507  else if (VarNames[i]=="PTCorrectedMass" )
508  {
509  numberOfPoints=100;
510  lowerBin=0.;
511  higerBin=10.;
512  }
513 
514  pTree->cd( "/" + name() + "/" + _FlavourTagInputsCollectionNames[iInputCollection]);
515  if( !pTree->cd( "bJets" ) )
516  {
517  pTree->mkdir( "bJets" ) ;
518  pTree->cd( "bJets" ) ;
519  }
520 
521  _inputsHistogramsBJets[iInputCollection][VarNames[i]]=pHistogramFactory->createHistogram1D( VarNames[i], numberOfPoints, lowerBin, higerBin );
522 
523  pTree->cd( "/" + name() + "/" + _FlavourTagInputsCollectionNames[iInputCollection]);
524  if( !pTree->cd( "cJets" ) )
525  {
526  pTree->mkdir( "cJets" ) ;
527  pTree->cd( "cJets" ) ;
528  }
529 
530  _inputsHistogramsCJets[iInputCollection][VarNames[i]]=pHistogramFactory->createHistogram1D( VarNames[i], numberOfPoints, lowerBin, higerBin );
531 
532  pTree->cd( "/" + name() + "/" + _FlavourTagInputsCollectionNames[iInputCollection]);
533  if( !pTree->cd( "udsJets/" ) )
534  {
535  pTree->mkdir( "udsJets/" ) ;
536  pTree->cd( "udsJets/" ) ;
537  }
538 
539  _inputsHistogramsUDSJets[iInputCollection][VarNames[i]]=pHistogramFactory->createHistogram1D( VarNames[i], numberOfPoints, lowerBin, higerBin );
540 
541  }//end of histogram creation
542  }
543 
544  if (isFirstEvent()) {
545  //We'd like to make zoomed histograms of some of the flavour tag inputs too
546  for (size_t i = 0;i < _ZoomedVarNames.size();++i) {
547 
548  std::string zoomed_name = _ZoomedVarNames[i] + " (zoomed)";
549 
550  pTree->cd( "/" + name() + "/" + _FlavourTagInputsCollectionNames[iInputCollection]);
551  if( !pTree->cd( "bJets" ) )
552  {
553  pTree->mkdir( "bJets" ) ;
554  pTree->cd( "bJets" ) ;
555  }
556 
557  _zoomedInputsHistogramsBJets[iInputCollection][zoomed_name] = pHistogramFactory->createHistogram1D( zoomed_name, 100, -10., 20.);
558 
559  pTree->cd( "/" + name() + "/" + _FlavourTagInputsCollectionNames[iInputCollection]);
560  if( !pTree->cd( "cJets" ) )
561  {
562  pTree->mkdir( "cJets" ) ;
563  pTree->cd( "cJets" ) ;
564  }
565 
566  _zoomedInputsHistogramsCJets[iInputCollection][zoomed_name] = pHistogramFactory->createHistogram1D( zoomed_name, 100, -10., 20.);
567 
568 
569  pTree->cd( "/" + name() + "/" + _FlavourTagInputsCollectionNames[iInputCollection]);
570  if( !pTree->cd( "udsJets" ) )
571  {
572  pTree->mkdir( "udsJets" ) ;
573  pTree->cd( "udsJets" ) ;
574  }
575 
576  _zoomedInputsHistogramsUDSJets[iInputCollection][zoomed_name] = pHistogramFactory->createHistogram1D( zoomed_name, 100, -10., 20.);
577  }
578  }
579  }
580 }
581 
582 
583 void LCFIAIDAPlotProcessor::processEvent( LCEvent* pEvent )
584 {
585 
586  // Make sure that "processRunHeader" has been called for this run (see the comment in that method).
587  if( (_lastRunHeaderProcessed != pEvent->getRunNumber()) && (_suppressOutputForRun != pEvent->getRunNumber()) )
588  {
589  std::cerr << __FILE__ << "(" << __LINE__ << "): processRunHeader() was not called for run " << pEvent->getRunNumber()
590  << " (did you use \"SkipNEvents\"?). The order of the information in the flavour tag collection(s) is going to be guessed." << std::endl;
591 
592  //Only want to do this once for this run, so set a marker that this run has been done
593  _suppressOutputForRun=pEvent->getRunNumber();
594 
595  // Just assume that the elements are in the order "BTag", "CTag", "BCTag"
596  std::map<std::string,unsigned int> guessedOrder;
597  guessedOrder["BTag"]=0; guessedOrder["CTag"]=1; guessedOrder["BCTag"]=2;
598  _IndexOfForEachTag.clear();
599  for (unsigned int iTag=0; iTag < _FlavourTagCollectionNames.size(); ++iTag) _IndexOfForEachTag.push_back( guessedOrder );
600  }
601 
602 
603  // Try and get the jet collection. If unable, show why. No need to worry about quitting
604  // because getNumberOfElements will return zero, so flow will never go into the loop.
605  // TypesafeCollection is just a wrapper around LCCollection with more error checking and
606  // things. Just makes the code a bit easier to read (in my opinion).
607 
608  TypesafeCollection<lcio::ReconstructedParticle> jetCollection( pEvent, _JetCollectionName );
609 
610 
611 
612  //apply any cuts on the event here
613  if( PassesEventCuts(pEvent) )
614  {
615 
616  ReconstructedParticle* pJet;
617  //loop over the jets
618  for( int jetNumber=0; jetNumber<jetCollection.getNumberOfElements(); ++jetNumber )
619  {
620  pJet=jetCollection.getElementAt(jetNumber);
621 
622  //only do anything if the jet passes the jet cuts
623  if( PassesJetCuts(pJet) )
624  {
625 
626  if (_MakePurityEfficiencyPlots || _PrintPurityEfficiencyValues) FillTagPlots( pEvent, jetNumber );
627  if (_MakeAdditionalPlots) FillVertexChargePlots( pEvent, jetNumber );
628  if (_MakeTuple) FillInputsPlots( pEvent, jetNumber );
629  if (_MakeAdditionalPlots) FillVertexPlots( pEvent, jetNumber );
630  //Look at the track-vertex matching
631  if (_PrintTrackVertexOutput) FillZVRESTable(pEvent);
632 
633  }
634  }
635  }
636 }
637 
638 
639 
640 void LCFIAIDAPlotProcessor::check( LCEvent* pEvent )
641 {
642 }
643 
644 
645 
646 void LCFIAIDAPlotProcessor::end()
647 {
648  if (_MakePurityEfficiencyPlots||_PrintPurityEfficiencyValues) CalculateIntegralAndBackgroundPlots();
649  if (_MakePurityEfficiencyPlots) CalculateEfficiencyPurityPlots();
650  if (_MakeAdditionalPlots) CalculateAdditionalPlots();
651  if (_PrintPurityEfficiencyValues) PrintNNOutput();
652  if (_PrintTrackVertexOutput) PrintZVRESTable();
653 
654 }
655 
656 // IMPORTANT - If you change the cuts make sure you change the line below to show the changes in the docs
658 bool LCFIAIDAPlotProcessor::PassesEventCuts( LCEvent* pEvent )
659 {
660  //
661  // No event cuts at present
662  //
663 
664  return true;
665 
666 }
667 // IMPORTANT - If you change the cuts make sure you change the line below to show the changes in the docs
668 
669 bool LCFIAIDAPlotProcessor::PassesJetCuts( ReconstructedParticle* pJet )
670 {
671  //
672  // This cut added on the suggestion of Sonja Hillert 12/Jan/07.
673  //
674  // Selects jets for which the cosine of the jet polar
675  // angle theta for all jets is not too large.
676  //
677  // Make something that's easy to search for to track down erroneous cuts:
678  // (too many bad experiences of long forgotten about cuts hiding somewhere)
679  // GREPTAG_CUT : Jet cut on abs(cos(theta)) of jet axis
680  //
681 
682 
683  const double* jetMomentum=pJet->getMomentum();
684 
685  double momentumMagnitude = sqrt(pow(jetMomentum[0],2)+pow(jetMomentum[1],2)+pow(jetMomentum[2],2));
686 
687  double cosTheta = jetMomentum[2] / momentumMagnitude;
688  if( fabs(cosTheta) < _CosThetaJetMin || fabs(cosTheta) > _CosThetaJetMax ) return false;
689 
690  if (momentumMagnitude > _PJetMax || momentumMagnitude < _PJetMin) return false;
691 
692 
693  // If control gets to this point then the jet has passed
694  return true;
695 }
696 
697 
698 void LCFIAIDAPlotProcessor::CalculateIntegralAndBackgroundPlots() {
699 
700  AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
701  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
702 
703  _pBJetBTagIntegral.resize( _FlavourTagCollectionNames.size() );
704  _pCJetBTagIntegral.resize( _FlavourTagCollectionNames.size() );
705  _pLightJetBTagIntegral.resize( _FlavourTagCollectionNames.size() );
706  _pBJetCTagIntegral.resize( _FlavourTagCollectionNames.size() );
707  _pCJetCTagIntegral.resize( _FlavourTagCollectionNames.size() );
708  _pLightJetCTagIntegral.resize( _FlavourTagCollectionNames.size() );
709  _pBJetBCTagIntegral.resize( _FlavourTagCollectionNames.size() );
710  _pCJetBCTagIntegral.resize( _FlavourTagCollectionNames.size() );
711  _pLightJetBCTagIntegral.resize( _FlavourTagCollectionNames.size() );
712 
713 
714  for (unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection )
715  {
716 
717  for (unsigned int iVertexCat=1; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
718  //sum over the different vertex catagories, this information goes into the "AnyNumberOfVertices" directory
719 
720  _pBJetBTag[iTagCollection][_VertexCatNames[0]] -> add(*_pBJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]);
721  _pBJetCTag[iTagCollection][_VertexCatNames[0]] -> add(*_pBJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
722  _pBJetBCTag[iTagCollection][_VertexCatNames[0]] -> add(*_pBJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
723  _pCJetBTag[iTagCollection][_VertexCatNames[0]] -> add(*_pCJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]);
724  _pCJetCTag[iTagCollection][_VertexCatNames[0]] -> add(*_pCJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
725  _pCJetBCTag[iTagCollection][_VertexCatNames[0]] -> add(*_pCJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
726  _pLightJetBTag[iTagCollection][_VertexCatNames[0]] -> add(*_pLightJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]);
727  _pLightJetCTag[iTagCollection][_VertexCatNames[0]] -> add(*_pLightJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
728  _pLightJetBCTag[iTagCollection][_VertexCatNames[0]] -> add(*_pLightJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
729  }
730 
731  for (unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
732  //add up all the background values
733 
734  pTree->cd( "/" + name() + "/" + _FlavourTagCollectionNames[iTagCollection]+ "/" +_NumVertexCatDir[iVertexCat]);
735 
736  _pBTagBackgroundValues[iTagCollection][_VertexCatNames[iVertexCat]] = pHistogramFactory->add("Numbers of non-B jets by B-tag NN value. ("+ _VertexCatNames[iVertexCat]+")",*_pLightJetBTag[iTagCollection][_VertexCatNames[iVertexCat]],*_pCJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]);
737  _pCTagBackgroundValues[iTagCollection][_VertexCatNames[iVertexCat]] = pHistogramFactory->add("Numbers of non-C jets by C-tag NN value. ("+ _VertexCatNames[iVertexCat]+")",*_pLightJetCTag[iTagCollection][_VertexCatNames[iVertexCat]],*_pBJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
738  _pBCTagBackgroundValues[iTagCollection][_VertexCatNames[iVertexCat]] = pHistogramFactory->add("Numbers of non-C jets by BC-tag NN value. ("+ _VertexCatNames[iVertexCat]+")",*_pLightJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]],*_pBJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
739 
740  }
741 
742  }
743 }
744 
745 
746 
747 void LCFIAIDAPlotProcessor::CalculateEfficiencyPurityPlots()
748 {
749  AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
750  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
751 
752 #ifdef USING_JAIDA
753  AIDA::IDataPointSetFactory* pDataPointSetFactory=marlin::AIDAProcessor::dataPointSetFactory(this);
754 #endif
755 
756 
757  //now calculate the efficiencies, leakage rate and purity
758  for (unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection )
759  {
760 
761  for (unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
762 
763  pTree->cd( "/" + name() + "/" + _FlavourTagCollectionNames[iTagCollection]+ "/" +_NumVertexCatDir[iVertexCat] );
764 
765  std::string nvname = _VertexCatNames[iVertexCat];
766 
767 #ifdef USING_JAIDA
768  AIDA::IDataPointSet* _pBJetBTagEfficiency = CreateEfficiencyPlot( _pBJetBTag[iTagCollection][nvname] , pDataPointSetFactory->create("B-Tag efficiency ("+ nvname +")",2));
769  AIDA::IDataPointSet* _pCJetCTagEfficiency = CreateEfficiencyPlot( _pCJetCTag[iTagCollection][nvname] , pDataPointSetFactory->create("C-Tag efficiency ("+ nvname +")",2));
770  AIDA::IDataPointSet* _pCJetBCTagEfficiency = CreateEfficiencyPlot( _pCJetBCTag[iTagCollection][nvname] , pDataPointSetFactory->create("BC-Tag efficiency ("+ nvname +")",2));
771 #endif
772 
773 
774  _pBJetBTagIntegral[iTagCollection][nvname] =
775  CreateIntegralHistogram( _pBJetBTag[iTagCollection][nvname],
776  pHistogramFactory->createHistogram1D("B-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
777  _pBJetBTag[iTagCollection][nvname]->axis().bins(),_pBJetBTag[iTagCollection][nvname]->axis().lowerEdge(),_pBJetBTag[iTagCollection][nvname]->axis().upperEdge()));
778 
779  _pCJetBTagIntegral[iTagCollection][nvname] =
780  CreateIntegralHistogram( _pCJetBTag[iTagCollection][nvname],
781  pHistogramFactory->createHistogram1D("C-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
782  _pCJetBTag[iTagCollection][nvname]->axis().bins(),_pCJetBTag[iTagCollection][nvname]->axis().lowerEdge(),_pCJetBTag[iTagCollection][nvname]->axis().upperEdge()));
783 
784  _pLightJetBTagIntegral[iTagCollection][nvname] =
785  CreateIntegralHistogram( _pLightJetBTag[iTagCollection][nvname],
786  pHistogramFactory->createHistogram1D("Light-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
787  _pLightJetBTag[iTagCollection][nvname]->axis().bins(),_pLightJetBTag[iTagCollection][nvname]->axis().lowerEdge(),_pLightJetBTag[iTagCollection][nvname]->axis().upperEdge()));
788 
789  _pBJetCTagIntegral[iTagCollection][nvname] =
790  CreateIntegralHistogram( _pBJetCTag[iTagCollection][nvname],
791  pHistogramFactory->createHistogram1D("B-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
792  _pBJetCTag[iTagCollection][nvname]->axis().bins(),_pBJetCTag[iTagCollection][nvname]->axis().lowerEdge(),_pBJetCTag[iTagCollection][nvname]->axis().upperEdge()));
793 
794  _pCJetCTagIntegral[iTagCollection][nvname] =
795  CreateIntegralHistogram( _pCJetCTag[iTagCollection][nvname],
796  pHistogramFactory->createHistogram1D("C-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
797  _pCJetCTag[iTagCollection][nvname]->axis().bins(),_pCJetCTag[iTagCollection][nvname]->axis().lowerEdge(),_pCJetCTag[iTagCollection][nvname]->axis().upperEdge()));
798 
799  _pLightJetCTagIntegral[iTagCollection][nvname] =
800  CreateIntegralHistogram( _pLightJetCTag[iTagCollection][nvname],
801  pHistogramFactory->createHistogram1D("Light-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
802  _pLightJetCTag[iTagCollection][nvname]->axis().bins(),_pLightJetCTag[iTagCollection][nvname]->axis().lowerEdge(),_pLightJetCTag[iTagCollection][nvname]->axis().upperEdge()));
803 
804  _pBJetBCTagIntegral[iTagCollection][nvname] =
805  CreateIntegralHistogram( _pBJetBCTag[iTagCollection][nvname],
806  pHistogramFactory->createHistogram1D("B-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
807  _pBJetBCTag[iTagCollection][nvname]->axis().bins(),_pBJetBCTag[iTagCollection][nvname]->axis().lowerEdge(),_pBJetBCTag[iTagCollection][nvname]->axis().upperEdge()));
808 
809  _pCJetBCTagIntegral[iTagCollection][nvname] =
810  CreateIntegralHistogram( _pCJetBCTag[iTagCollection][nvname],
811  pHistogramFactory->createHistogram1D("C-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +") (DON'T TRUST ERRORS!)",
812  _pCJetBCTag[iTagCollection][nvname]->axis().bins(),_pCJetBCTag[iTagCollection][nvname]->axis().lowerEdge(),_pCJetBCTag[iTagCollection][nvname]->axis().upperEdge()));
813 
814  _pLightJetBCTagIntegral[iTagCollection][nvname] =
815  CreateIntegralHistogram( _pLightJetBCTag[iTagCollection][nvname],
816  pHistogramFactory->createHistogram1D("Light-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +")",
817  _pLightJetBCTag[iTagCollection][nvname]->axis().bins(),_pLightJetBCTag[iTagCollection][nvname]->axis().lowerEdge(),_pLightJetBCTag[iTagCollection][nvname]->axis().upperEdge()));
818 
819  //Examples of the integral plots - instead of histograms - the histogram calculate the errors wrongly
820 
821  //integralplots _pBJetBTagIntegral[iTagCollection][nvname] =
822  //integralplots// CreateIntegralPlot( _pBJetBTag[iTagCollection][nvname], pDataPointSetFactory->create("B-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +")",2));
823  //integralplots// _pCJetBTagIntegral[iTagCollection][nvname] =
824  //integralplots// CreateIntegralPlot( _pCJetBTag[iTagCollection][nvname], pDataPointSetFactory->create("C-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +")",2));
825  //integralplots// _pLightJetBTagIntegral[iTagCollection][nvname] =
826  //integralplots// CreateIntegralPlot( _pLightJetBTag[iTagCollection][nvname], pDataPointSetFactory->create("Light-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +")",2));
827  //integralplots// _pBJetCTagIntegral[iTagCollection][nvname] =
828  //integralplots// CreateIntegralPlot( _pBJetCTag[iTagCollection][nvname], pDataPointSetFactory->create("B-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +")",2));
829  //integralplots// _pCJetCTagIntegral[iTagCollection][nvname] =
830  //integralplots// CreateIntegralPlot( _pCJetCTag[iTagCollection][nvname], pDataPointSetFactory->create("C-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +")",2));
831  //integralplots// _pLightJetCTagIntegral[iTagCollection][nvname] =
832  //integralplots// CreateIntegralPlot( _pLightJetCTag[iTagCollection][nvname], pDataPointSetFactory->create("Light-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +")",2));
833  //integralplots// _pBJetBCTagIntegral[iTagCollection][nvname] =
834  //integralplots// CreateIntegralPlot( _pBJetBCTag[iTagCollection][nvname], pDataPointSetFactory->create("B-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +")",2));
835  //integralplots// _pCJetBCTagIntegral[iTagCollection][nvname] =
836  //integralplots// CreateIntegralPlot( _pCJetBCTag[iTagCollection][nvname], pDataPointSetFactory->create("C-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +")",2));
837  //integralplots// _pLightJetBCTagIntegral[iTagCollection][nvname] =
838  //integralplots// CreateIntegralPlot( _pLightJetBCTag[iTagCollection][nvname], pDataPointSetFactory->create("Light-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +")",2));
839 
840 #ifdef USING_JAIDA
841  AIDA::IDataPointSet* _pBJetBTagPurity = CreatePurityPlot( _pBJetBTag[iTagCollection][nvname], _pBTagBackgroundValues[iTagCollection][nvname] , pDataPointSetFactory->create("B-Jet purity for B-Tag ("+ nvname +")",2));
842  AIDA::IDataPointSet* _pCJetCTagPurity = CreatePurityPlot( _pCJetCTag[iTagCollection][nvname], _pCTagBackgroundValues[iTagCollection][nvname] , pDataPointSetFactory->create("C-Jet purity for C-Tag ("+ nvname +")",2));
843  AIDA::IDataPointSet* _pCJetBCTagPurity = CreatePurityPlot( _pCJetBCTag[iTagCollection][nvname], _pBJetBCTag[iTagCollection][nvname], pDataPointSetFactory->create("C-Jet purity for BC-Tag ("+ nvname +")",2));
844 
845  AIDA::IDataPointSet* _pCJetBTagLeakage = CreateLeakageRatePlot( _pCJetBTag[iTagCollection][nvname], pDataPointSetFactory->create("C-Jets: Leakage Rate into B-Tag Sample ("+ nvname +")",2));
846  AIDA::IDataPointSet* _pLightJetBTagLeakage = CreateLeakageRatePlot( _pLightJetBTag[iTagCollection][nvname], pDataPointSetFactory->create("Light-Jets: Leakage Rate into B-Tag Sample ("+ nvname +")",2));
847  AIDA::IDataPointSet* _pBJetCTagLeakage = CreateLeakageRatePlot( _pBJetCTag[iTagCollection][nvname], pDataPointSetFactory->create("B-Jets: Leakage Rate into C-Tag Sample ("+ nvname +")",2));
848  AIDA::IDataPointSet* _pLightJetCTagLeakage = CreateLeakageRatePlot( _pLightJetCTag[iTagCollection][nvname], pDataPointSetFactory->create("Light-Jets: Leakage Rate into C-Tag Sample ("+ nvname +")",2));
849  AIDA::IDataPointSet* _pBJetBCTagLeakage = CreateLeakageRatePlot( _pBJetBCTag[iTagCollection][nvname], pDataPointSetFactory->create("B-Jets: Leakage Rate into BC-Tag Sample ("+ nvname +")",2));
850  AIDA::IDataPointSet* _pLightJetBCTagLeakage = CreateLeakageRatePlot( _pLightJetBCTag[iTagCollection][nvname], pDataPointSetFactory->create("Light-Jets: Leakage Rate into BC-Tag Sample ("+ nvname +")",2));
851  AIDA::IDataPointSet* _pNonBJetBTagLeakage = CreateLeakageRatePlot( _pBTagBackgroundValues[iTagCollection][nvname], pDataPointSetFactory->create("C-Jets: Leakage Rate into B-Tag Sample ("+ nvname +")",2));
852  AIDA::IDataPointSet* _pNonCJetCTagLeakage = CreateLeakageRatePlot( _pCTagBackgroundValues[iTagCollection][nvname], pDataPointSetFactory->create("Light-Jets: Leakage Rate into B-Tag Sample ("+ nvname +")",2));
853  AIDA::IDataPointSet* _pNonCJetBCTagLeakage = CreateLeakageRatePlot( _pBCTagBackgroundValues[iTagCollection][nvname], pDataPointSetFactory->create("B-Jets: Leakage Rate into C-Tag Sample ("+ nvname +")",2));
854 
855  CreateXYPlot(_pBJetBTagEfficiency, _pBJetBTagPurity, pDataPointSetFactory->create("Purity-Efficiency for B-Tag ("+ nvname +")",2), 1, 1);
856  CreateXYPlot(_pCJetCTagEfficiency, _pCJetCTagPurity, pDataPointSetFactory->create("Purity-Efficiency for C-Tag ("+ nvname +")",2), 1, 1);
857  CreateXYPlot(_pCJetBCTagEfficiency, _pCJetBCTagPurity, pDataPointSetFactory->create("Purity-Efficiency for BC-Tag ("+ nvname +")",2), 1, 1);
858 
859  CreateXYPlot(_pBJetBTagEfficiency, _pNonBJetBTagLeakage, pDataPointSetFactory->create("Leakage Rate-Efficiency for B-Tag ("+ nvname +")",2), 1, 1);
860  CreateXYPlot(_pCJetCTagEfficiency, _pNonCJetCTagLeakage, pDataPointSetFactory->create("Leakage Rate-Efficiency for C-Tag ("+ nvname +")",2), 1, 1);
861  CreateXYPlot(_pCJetBCTagEfficiency,_pNonCJetBCTagLeakage, pDataPointSetFactory->create("Leakage Rate-Efficiency for BC-Tag ("+ nvname +")",2), 1, 1);
862 
863  CreateXYPlot(_pBJetBTagEfficiency, _pCJetBTagLeakage, pDataPointSetFactory->create("Leakage Rate (C into B) vs Efficiency for B-Tag ("+ nvname +")",2), 1, 1);
864  CreateXYPlot(_pBJetBTagEfficiency, _pLightJetBTagLeakage, pDataPointSetFactory->create("Leakage Rate (Light into B) Efficiency for B-Tag ("+ nvname +")",2), 1, 1);
865  CreateXYPlot(_pCJetCTagEfficiency, _pBJetCTagLeakage, pDataPointSetFactory->create("Leakage Rate (B into C) vs Efficiency for C-Tag ("+ nvname +")",2), 1, 1);
866  CreateXYPlot(_pCJetCTagEfficiency, _pLightJetCTagLeakage, pDataPointSetFactory->create("Leakage Rate (Light into C) Efficiency for C-Tag ("+ nvname +")",2), 1, 1);
867  CreateXYPlot(_pCJetBCTagEfficiency, _pBJetBCTagLeakage, pDataPointSetFactory->create("Leakage Rate (B into BC) vs Efficiency for BC-Tag ("+ nvname +")",2), 1, 1);
868  CreateXYPlot(_pCJetCTagEfficiency, _pLightJetBCTagLeakage, pDataPointSetFactory->create("Leakage Rate (Light into BC) Efficiency for BC-Tag ("+ nvname +")",2), 1, 1);
869 #endif
870 
871  }//end of loop of different number of vertices
872 
873  }//end loop over iTagCollection
874 
875 }
876 
877 
878 
879 void LCFIAIDAPlotProcessor::CalculateAdditionalPlots()
880 {
881 
882  //AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
883  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
884 
885 #ifdef USING_JAIDA
886  AIDA::IDataPointSetFactory* pDataPointSetFactory=marlin::AIDAProcessor::dataPointSetFactory(this);
887 #endif
888 
889 
890  //now calculate the efficiencies, leakage rate and purity
891  for (unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection )
892  {
893 
894  for (unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
895 
896  pTree->cd( "/" + name() + "/" + _FlavourTagCollectionNames[iTagCollection]+ "/" +_NumVertexCatDir[iVertexCat] );
897 
898  std::string nvname = _VertexCatNames[iVertexCat];
899 
900 
901  }//end of loop of different number of vertices
902 
903  pTree->cd( "/" + name() + "/" + _FlavourTagCollectionNames[iTagCollection]+ "/");
904 
905 #ifdef USING_JAIDA
906  CreateEfficiencyPlot2( _pCDecayLengthAll[iTagCollection], _pCDecayLengthTwoVertices[iTagCollection], pDataPointSetFactory->create("Efficiency to reconstruct secondary vertex vs true decay length for true C-events",2));
907  CreateEfficiencyPlot2( _pBDecayLengthAll[iTagCollection], _pBDecayLengthTwoVertices[iTagCollection], pDataPointSetFactory->create("Efficiency to reconstruct secondary vertex vs true decay length for true B-events",2));
908 #endif
909 
910  }//end loop over iTagCollection
911 
912  //create vertex charge leakage rate plots
913  pTree->cd( "/" + name() + "/VertexChargePlots/");
914 
915 
916 #ifdef USING_JAIDA
917  AIDA::IDataPointSet* pBJetVtxChargeDPS = pDataPointSetFactory->create("B-Jets: Vertex Charge Leakage",2);
918  AIDA::IDataPointSet* pCJetVtxChargeDPS = pDataPointSetFactory->create("C-Jets: Vertex Charge Leakage",2);
919  CreateVertexChargeLeakagePlot(pBJetVtxChargeDPS, pCJetVtxChargeDPS);
920 #endif
921 
922 }
923 
924 
925 void LCFIAIDAPlotProcessor::FillInputsPlots( LCEvent* pEvent, unsigned int jetNumber )
926 {
927  int jetType=FindTrueJetType( pEvent, jetNumber );
928  if( jetType==0 ) return;
929 
930 
931 
932  for (unsigned int iInputsCollection=0; iInputsCollection < _FlavourTagInputsCollectionNames.size(); ++iInputsCollection )
933  {
934  TypesafeCollection<lcio::LCFloatVec> inputsCollection( pEvent, _FlavourTagInputsCollectionNames[iInputsCollection] );
935  if( !inputsCollection.is_valid() )
936  {
937  std::cerr << _FlavourTagCollectionNames[iInputsCollection] << "is not valid, not filling flavour tag input plots" << std::endl;
938  std::cerr << "set MakeAdditionalPlots to false to remove this warning" << std::endl;
939  }
940  else
941  {
942  //Do stuff...
943  lcio::LCFloatVec* pInputs=inputsCollection.getElementAt( jetNumber );
944  if( !pInputs )
945  {
946  }
947  else
948  {
949  //ok everything is okay with the data
950 #ifdef USING_JAIDA
951  //I have a problem with tuples in RAIDA
952 
953  int CQVtx = FindCQVtx(pEvent, jetNumber);
954  int BQVtx = FindBQVtx(pEvent, jetNumber);
955 
956  if (_MakeTuple && iInputsCollection==0) {
957 
958  //this could probably be done automatically
959 
960  int NumVertices = int((*pInputs)[_InputsIndex[iInputsCollection]["NumVertices"]]);
961  int NumTracksInVertices = int((*pInputs)[_InputsIndex[iInputsCollection] ["NumTracksInVertices"]]);
962  float D0Significance1=(*pInputs)[_InputsIndex[iInputsCollection] ["D0Significance1"]];
963  float D0Significance2=(*pInputs)[_InputsIndex[iInputsCollection] ["D0Significance2"]];
964  float DecayLength=(*pInputs)[_InputsIndex[iInputsCollection] ["DecayLength"]];
965  float DecayLength_SeedToIP=(*pInputs)[_InputsIndex[iInputsCollection]["DecayLength(SeedToIP)"]];
966  float DecayLengthSignificance=(*pInputs)[_InputsIndex[iInputsCollection] ["DecayLengthSignificance"]];
967  float JointProbRPhi=(*pInputs)[_InputsIndex[iInputsCollection] ["JointProbRPhi"]];
968  float JointProbZ=(*pInputs)[_InputsIndex[iInputsCollection] ["JointProbZ"]];
969  float Momentum1=(*pInputs)[_InputsIndex[iInputsCollection] ["Momentum1"]];
970  float Momentum2=(*pInputs)[_InputsIndex[iInputsCollection] ["Momentum2"]];
971  float PTCorrectedMass=(*pInputs)[_InputsIndex[iInputsCollection] ["PTCorrectedMass"]];
972  float RawMomentum=(*pInputs)[_InputsIndex[iInputsCollection] ["RawMomentum"]];
973  float SecondaryVertexProbability=(*pInputs)[_InputsIndex[iInputsCollection]["SecondaryVertexProbability"]];
974  float Z0Significance1=(*pInputs)[_InputsIndex[iInputsCollection] ["Z0Significance1"]];
975  float Z0Significance2=(*pInputs)[_InputsIndex[iInputsCollection] ["Z0Significance2"]];
976 
977  if (_pMyTuple) {
978 
979  _pMyTuple->fill( 0, jetType );
980  _pMyTuple->fill( 1, NumVertices );
981  _pMyTuple->fill( 2, NumTracksInVertices );
982  _pMyTuple->fill( 3, D0Significance1);
983  _pMyTuple->fill( 4, D0Significance2);
984  _pMyTuple->fill( 5, DecayLength);
985  _pMyTuple->fill( 6, DecayLength_SeedToIP);
986  _pMyTuple->fill( 7, DecayLengthSignificance);
987  _pMyTuple->fill( 8, JointProbRPhi);
988  _pMyTuple->fill( 9, JointProbZ);
989  _pMyTuple->fill( 10, Momentum1);
990  _pMyTuple->fill( 11, Momentum2);
991  _pMyTuple->fill( 12, PTCorrectedMass);
992  _pMyTuple->fill( 13, RawMomentum);
993  _pMyTuple->fill( 14, SecondaryVertexProbability);
994  _pMyTuple->fill( 15, Z0Significance1);
995  _pMyTuple->fill( 16 ,Z0Significance2);
996 
997  _pMyTuple->fill( 17, BQVtx );
998  _pMyTuple->fill( 18, CQVtx );
999 
1000  _pMyTuple->addRow();
1001 
1002  }//endif _pMyTuple
1003  }
1004 #endif
1005 
1006  for( std::map<std::string,unsigned int>::iterator iTagNames=_InputsIndex[iInputsCollection].begin();
1007  iTagNames!=_InputsIndex[iInputsCollection].end(); ++iTagNames ) {
1008 
1009  double input=(*pInputs)[(*iTagNames).second];
1010 
1011  //if the quantity relates to the second vertex, and there is no second vertex, then don't plot it
1012  if (! ((*pInputs)[_InputsIndex[iInputsCollection]["NumVertices"]] < 2 &&
1013  ((*iTagNames).first == "DecayLength" || (*iTagNames).first == "RawMomentum" ||
1014  (*iTagNames).first == "SecondaryVertexProbability" || (*iTagNames).first == "PTCorrectedMass" ||
1015  (*iTagNames).first == "DecayLength(SeedToIP)" || (*iTagNames).first == "DecayLengthSignificance") )) {
1016 
1017  if( jetType==B_JET ) _inputsHistogramsBJets[iInputsCollection][(*iTagNames).first]->fill(input);
1018  else if( jetType==C_JET ) _inputsHistogramsCJets[iInputsCollection][(*iTagNames).first]->fill(input);
1019  else _inputsHistogramsUDSJets[iInputsCollection][(*iTagNames).first]->fill(input);
1020  }
1021 
1022  //fill a few extra histograms created by hand
1023  for (std::vector<std::string>::const_iterator iter = _ZoomedVarNames.begin() ; iter != _ZoomedVarNames.end(); iter++){
1024  if ((*iTagNames).first == *iter) {
1025  std::string zoomed_name = (*iTagNames).first + " (zoomed)";
1026  if( jetType==B_JET ) _zoomedInputsHistogramsBJets[iInputsCollection][zoomed_name]->fill(input);
1027  else if( jetType==C_JET ) _zoomedInputsHistogramsCJets[iInputsCollection][zoomed_name]->fill(input);
1028  else _zoomedInputsHistogramsUDSJets[iInputsCollection][zoomed_name]->fill(input);
1029 
1030  }
1031  }
1032 
1033  }
1034 
1035  }
1036  }
1037  }
1038 }
1039 
1040 
1041 
1042 void LCFIAIDAPlotProcessor::FillTagPlots( LCEvent* pEvent, unsigned int jetNumber)
1043 {
1044  int jetType=FindTrueJetType( pEvent, jetNumber );
1045  if( jetType==0 ) return;
1046 
1047  //needs to be tidied up
1048  TypesafeCollection<lcio::ReconstructedParticle> jetCollection( pEvent, _JetCollectionName );
1049  ReconstructedParticle* pJet;
1050  pJet=jetCollection.getElementAt(jetNumber);
1051 
1052  for (unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection)
1053  {
1054  TypesafeCollection<lcio::LCFloatVec> tagCollection( pEvent, _FlavourTagCollectionNames[iTagCollection] );
1055  if( !tagCollection.is_valid() )
1056  {
1057  std::cerr << _FlavourTagCollectionNames[iTagCollection] << "is not valid, not filling tag plots" << std::endl;
1058  std::cerr << "set MakeAdditionalPlots to false to remove this warning" << std::endl;
1059 
1060  }
1061  else
1062  {
1063  lcio::LCFloatVec* pJetFlavourTags=tagCollection.getElementAt( jetNumber );
1064  if( !pJetFlavourTags )
1065  {
1066  }
1067  else
1068  {
1069 
1070  //const double* jetMomentum = pJet->getMomentum();
1071  //double cosTheta = jetMomentum[2] / sqrt(pow(jetMomentum[0],2)+pow(jetMomentum[1],2)+pow(jetMomentum[2],2));
1072 
1073  double bTag= (*pJetFlavourTags)[_IndexOfForEachTag[iTagCollection]["BTag"]];
1074  double cTag= (*pJetFlavourTags)[_IndexOfForEachTag[iTagCollection]["CTag"]];
1075  double cTagBBack= (*pJetFlavourTags)[_IndexOfForEachTag[iTagCollection]["BCTag"]];
1076  unsigned int NumVertices = FindNumVertex(pEvent, jetNumber, iTagCollection);
1077  //int CQVtx = FindCQVtx(pEvent, jetNumber);
1078  //int BQVtx = FindBQVtx(pEvent, jetNumber);
1079  //int trueJetCharge = int(FindTrueJetHadronCharge(pEvent,jetNumber));
1080 
1081  std::string nvname = _VertexCatNames[ (NumVertices>=N_VERTEX_CATEGORIES) ? (N_VERTEX_CATEGORIES) : (NumVertices)];
1082 
1083  if( jetType==B_JET ) {
1084 
1085  if( bTag<=1 && bTag>=0 )
1086  {
1087  _pBJetBTag[iTagCollection][nvname]->fill( bTag );
1088  }
1089  else
1090  {
1091  _pBJetBTag[iTagCollection][nvname]->fill( -0.005 );
1092  }
1093 
1094  if( cTag<=1 && cTag>=0 )
1095  {
1096  _pBJetCTag[iTagCollection][nvname]->fill( cTag );
1097  }
1098  else
1099  {
1100  _pBJetCTag[iTagCollection][nvname]->fill( -0.005 );
1101  }
1102  if( cTagBBack<=1 && cTagBBack>=0 )
1103  {
1104  _pBJetBCTag[iTagCollection][nvname]->fill( cTagBBack );
1105  }
1106  else
1107  {
1108  _pBJetBCTag[iTagCollection][nvname]->fill( -0.005 );
1109  }
1110 
1111  } else if( jetType==C_JET ) {
1112 
1113  if( bTag<=1 && bTag>=0 )
1114  {
1115  _pCJetBTag[iTagCollection][nvname]->fill( bTag );
1116  }
1117  else
1118  {
1119  _pCJetBTag[iTagCollection][nvname]->fill( -0.005 );
1120  }
1121 
1122  if( cTag<=1 && cTag>=0 )
1123  {
1124  _pCJetCTag[iTagCollection][nvname]->fill( cTag );
1125  }
1126  else
1127  {
1128  _pCJetCTag[iTagCollection][nvname]->fill( -0.005 );
1129  }
1130 
1131  if( cTagBBack<=1 && cTagBBack>=0 )
1132  {
1133  _pCJetBCTag[iTagCollection][nvname]->fill( cTagBBack );
1134  }
1135  else
1136  {
1137  _pCJetBCTag[iTagCollection][nvname]->fill( -0.005 );
1138  }
1139  } else {
1140  if( bTag<=1 && bTag>=0 )
1141  {
1142  _pLightJetBTag[iTagCollection][nvname]->fill( bTag );
1143  }
1144  else
1145  {
1146  _pLightJetBTag[iTagCollection][nvname]->fill( -0.005 );
1147  }
1148 
1149  if( cTag<=1 && cTag>=0 )
1150  {
1151  _pLightJetCTag[iTagCollection][nvname]->fill( cTag );
1152  }
1153  else
1154  {
1155  _pLightJetCTag[iTagCollection][nvname]->fill( -0.005 );
1156  }
1157 
1158  if( cTagBBack<=1 && cTagBBack>=0 )
1159  {
1160  _pLightJetBCTag[iTagCollection][nvname]->fill( cTagBBack );
1161  }
1162  else
1163  {
1164  _pLightJetBCTag[iTagCollection][nvname]->fill( -0.005 );
1165  }
1166  }
1167  }
1168  }
1169  }
1170 }
1171 
1172 void LCFIAIDAPlotProcessor::CreateAdditionalPlots()
1173 {
1174  _pBJetCharge.resize( _FlavourTagCollectionNames.size() );
1175  _pCJetCharge.resize( _FlavourTagCollectionNames.size() );
1176  _pCDecayLengthAll.resize( _FlavourTagCollectionNames.size() );
1177  _pBDecayLengthAll.resize( _FlavourTagCollectionNames.size() );
1178  _pCDecayLengthTwoVertices.resize( _FlavourTagCollectionNames.size() );
1179  _pBDecayLengthTwoVertices.resize( _FlavourTagCollectionNames.size() );
1180 
1181 
1182  AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
1183  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
1184 
1185  if (!pTree->cd("/" + name() + "/ZVTOPPlots/")) {
1186  pTree->cd( "/" + name() + "/");
1187  pTree->mkdir("ZVTOPPlots/");
1188  pTree->cd( "ZVTOPPlots/");
1189  }
1190 
1191 
1192 
1193  _decayLengthBJet2D = pHistogramFactory->createHistogram2D( "B jets: Reconstructed secondary decay length vs MC B decay length (mm)",50,0.,20.,50,0.,20.);
1194  _decayLengthCJet2D = pHistogramFactory->createHistogram2D( "B jets: Reconstructed sec-ter decay length vs MC D decay length (mm)",50,0.,20.,50,0.,20.);
1195  _decayLengthBJetCloud2D = pHistogramFactory->createCloud2D( "B jets: Reconstructed secondary decay length vs MC B decay length (mm) (cloud)");
1196  _decayLengthCJetCloud2D = pHistogramFactory->createCloud2D( "B jets: Reconstructed sec-ter decay length vs MC D decay length (mm) (cloud)");
1197 
1198  _reconstructedSecondaryDecayLength = pHistogramFactory->createHistogram1D( "All jets: Reconstructed secondary decay length (mm)",100,0.,20.);
1199  _reconstructedSecTerDecayLength = pHistogramFactory->createHistogram1D( "All jets: Reconstructed secondary-tertiary decay length",100,0.,20.);
1200 
1201  _recoDecayLengthBJet = pHistogramFactory->createHistogram1D( "B jets: Reconstructed secondary decay length (mm)",50,0.,20.);
1202  _recoDecayLengthCJet = pHistogramFactory->createHistogram1D( "C jets: Reconstructed secondary decay length (mm)",50,0.,20.);
1203  _recoDecayLengthBCJet = pHistogramFactory->createHistogram1D( "B jets: Reconstructed secondary-tertiary decay length (mm)",50,0.,20.);
1204  _recoDecayLengthLightJet = pHistogramFactory->createHistogram1D( "Light jets: Reconstructed secondary decay length (mm)",50,0.,10.);
1205 
1206  _nVerticesBJet = pHistogramFactory->createHistogram1D("B jets: Number of reconstructed vertices",10,-0.5,9.5);
1207  _nVerticesCJet = pHistogramFactory->createHistogram1D( "C jets: Number of reconstructed vertices",10,-0.5,9.5);
1208  _nVerticesLightJet = pHistogramFactory->createHistogram1D( "Light jets: Number of reconstructed vertices",10,-0.5,9.5);
1209 
1210  _decayLengthBJetTrue = pHistogramFactory->createHistogram1D( "B jets: MC secondary decay length (mm)",50,0.,20.);;
1211  _decayLengthBCJetTrue = pHistogramFactory->createHistogram1D( "B jets: MC secondary-tertiary decay length (mm)",50,0.,20.);;
1212  _decayLengthBJetTrue = pHistogramFactory->createHistogram1D( "C jets: MC secondary decay length (mm)",50,0.,20.);;
1213 
1214 
1215  //some plots of vertex charge
1216  if (!pTree->cd("/" + name() + "/VertexChargePlots/")) {
1217  pTree->cd( "/" + name() + "/");
1218  pTree->mkdir("VertexChargePlots/");
1219  pTree->cd( "VertexChargePlots/");
1220  }
1221 
1222  _pBJetCharge2D = pHistogramFactory->createHistogram2D( "B Jets: Reconstructed Vertex Charge vs True Jet Charge",7,-3.5,+3.5,7,-3.5,+3.5);
1223  _pCJetCharge2D = pHistogramFactory->createHistogram2D( "C Jets: Reconstructed Vertex Charge vs True Jet Charge",7,-3.5,+3.5,7,-3.5,+3.5);
1224 
1225  _pBJetVertexCharge = pHistogramFactory->createHistogram1D( "B Jets: Reconstructed Vertex Charge",9,-4.5,+4.5);
1226  _pCJetVertexCharge = pHistogramFactory->createHistogram1D( "C Jets: Reconstructed Vertex Charge",9,-4.5,+4.5);
1227 
1228  _pCJetLeakageRate = pHistogramFactory->createHistogram1D("C Jets: Charged Leakage Rate (DON'T TRUST ERRORS)", N_JETANGLE_BINS,0.,1.);
1229  _pBJetLeakageRate = pHistogramFactory->createHistogram1D("B Jets: Charged Leakage Rate (DON'T TRUST ERRORS)", N_JETANGLE_BINS,0.,1.);
1230 
1231 
1232 
1233  for (unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection )
1234  {
1235 
1236  if (!pTree->cd( "/" + name() + "/" +_FlavourTagCollectionNames[iTagCollection] + "/")) {
1237  pTree->cd( "/" + name() + "/");
1238  if (!pTree->cd( _FlavourTagCollectionNames[iTagCollection])) {
1239  pTree->mkdir( _FlavourTagCollectionNames[iTagCollection] + "/");
1240  pTree->cd( _FlavourTagCollectionNames[iTagCollection]+ "/");
1241  }
1242  }
1243 
1244  if (!pTree->cd( "VertexChargePlots/")) {
1245  pTree->mkdir( "VertexChargePlots/");
1246  pTree->cd( "VertexChargePlots/");
1247  }
1248 
1249  _pCDecayLengthAll[iTagCollection] = pHistogramFactory->createHistogram1D( "True MC decay length of all reconstructed C-jets", 25, 0, 10.0 );
1250  _pBDecayLengthAll[iTagCollection] = pHistogramFactory->createHistogram1D( "True MC decay length of all reconstructed B-jets", 25, 0, 10.5 );
1251  _pCDecayLengthTwoVertices[iTagCollection] = pHistogramFactory->createHistogram1D( "True MC decay length of C-jets with more than one vertex", 25, 0, 10.0 );
1252  _pBDecayLengthTwoVertices[iTagCollection] = pHistogramFactory->createHistogram1D( "True MC decay length of B-jets with more than one vertex", 25, 0, 10.5 );
1253 
1254  _pBJetCharge[iTagCollection] = pHistogramFactory->createHistogram2D( "B Jets: Reconstructed Jet Charge vs True Jet Charge",7,-3.5,+3.5,7,-3.5,+3.5);
1255  _pCJetCharge[iTagCollection] = pHistogramFactory->createHistogram2D( "C Jets: Reconstructed Jet Charge vs True Jet Charge",7,-3.5,+3.5,7,-3.5,+3.5);
1256 
1257  }
1258 
1259  pTree->cd( "/" + name());
1260 
1261  if (!pTree->cd( "VertexPlots/")) {
1262  pTree->mkdir( "VertexPlots/");
1263  pTree->cd( "VertexPlots/");
1264  }
1265 
1266  _pVertexDistanceFromIP = pHistogramFactory->createHistogram1D( "Reconstructed Vertex distance from IP",100, 0., 10.);
1267  _pVertexPositionX = pHistogramFactory->createHistogram1D( "Non-primary vertex: x-position", 100, -10., 10.) ;
1268  _pVertexPositionY = pHistogramFactory->createHistogram1D( "Non-primary vertex: y-position", 100, -10., 10.);
1269  _pVertexPositionZ = pHistogramFactory->createHistogram1D( "Non-primary vertex: z-position", 100, -10., 10.);
1270 
1271  _pPrimaryVertexPullX = pHistogramFactory->createHistogram1D( "Non-primary vertex: x-pull", 100, -10., 10.);
1272  _pPrimaryVertexPullY = pHistogramFactory->createHistogram1D( "Non-primary vertex: y-pull", 100, -10., 10.);
1273  _pPrimaryVertexPullZ = pHistogramFactory->createHistogram1D( "Non-primary vertex: z-pull", 100, -10., 10.);
1274  _pPrimaryVertexPositionX = pHistogramFactory->createHistogram1D( "Primary vertex: x-position", 100, -10., 10.);
1275  _pPrimaryVertexPositionY = pHistogramFactory->createHistogram1D( "Primary vertex: y-position", 100, -10., 10.);
1276  _pPrimaryVertexPositionZ = pHistogramFactory->createHistogram1D( "Primary vertex: z-position", 100, -10., 10.);
1277 
1278 
1279 }
1280 
1281 
1282 void LCFIAIDAPlotProcessor::CreateTagPlots()
1283 {
1284  AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory( this );
1285  AIDA::ITree* pTree=marlin::AIDAProcessor::tree( this );
1286 
1287  _pBJetBTag.resize( _FlavourTagCollectionNames.size() );
1288  _pBJetCTag.resize( _FlavourTagCollectionNames.size() );
1289  _pBJetBCTag.resize( _FlavourTagCollectionNames.size() );
1290 
1291  _pCJetCTag.resize( _FlavourTagCollectionNames.size() );
1292  _pCJetBTag.resize( _FlavourTagCollectionNames.size() );
1293  _pCJetBCTag.resize( _FlavourTagCollectionNames.size() );
1294 
1295  _pLightJetBTag.resize( _FlavourTagCollectionNames.size() );
1296  _pLightJetCTag.resize( _FlavourTagCollectionNames.size() );
1297  _pLightJetBCTag.resize( _FlavourTagCollectionNames.size() );
1298 
1299  _pBTagBackgroundValues.resize( _FlavourTagCollectionNames.size() );
1300  _pCTagBackgroundValues.resize( _FlavourTagCollectionNames.size() );
1301  _pBCTagBackgroundValues.resize( _FlavourTagCollectionNames.size() );
1302 
1303 
1304  for (unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection )
1305  {
1306  for (unsigned int iVertexCat=0; iVertexCat < N_VERTEX_CATEGORIES+1; ++iVertexCat ){
1307  _pLightJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1308  _pLightJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1309  _pLightJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1310 
1311  _pBJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1312  _pBJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1313  _pBJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1314 
1315  _pCJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1316  _pCJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1317  _pCJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1318 
1319  _pBTagBackgroundValues[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1320  _pCTagBackgroundValues[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1321  _pBCTagBackgroundValues[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1322  }
1323  }
1324 
1325 
1326  for (unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection )
1327  {
1328 
1329  if (!pTree->cd( "/" + name() + "/" +_FlavourTagCollectionNames[iTagCollection] + "/")) {
1330  pTree->cd( "/" + name() + "/");
1331  if (!pTree->cd( _FlavourTagCollectionNames[iTagCollection])) {
1332  pTree->mkdir( _FlavourTagCollectionNames[iTagCollection] + "/");
1333  pTree->cd( _FlavourTagCollectionNames[iTagCollection]+ "/");
1334  }
1335  }
1336 
1337  for (unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ){
1338 
1339  std::string nvname = _VertexCatNames[iVertexCat];
1340 
1341 
1342  pTree->cd( "/" + name() + "/" +_FlavourTagCollectionNames[iTagCollection] + "/");
1343 
1344  if (!pTree->cd( _NumVertexCatDir[iVertexCat]+"/")) {
1345  pTree->mkdir( _NumVertexCatDir[iVertexCat]+"/");
1346  pTree->cd( _NumVertexCatDir[iVertexCat]+"/");
1347  }
1348 
1349  _pLightJetBTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D( "Numbers of light jets by B-tag NN value. ("+ nvname +")",_numberOfPoints , 0, 1.0 );
1350  _pLightJetCTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D( "Numbers of light jets by C-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
1351  _pLightJetBCTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D( "Numbers of light jets by BC-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
1352  _pBJetBTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D( "Numbers of B jets by B-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
1353  _pBJetCTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D( "Numbers of B jets by C-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
1354  _pBJetBCTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D( "Numbers of B jets by BC-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
1355  _pCJetBTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D( "Numbers of C jets by B-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
1356  _pCJetCTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D( "Numbers of C jets by C-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
1357  _pCJetBCTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D( "Numbers of C jets by BC-tag NN value. ("+ nvname +")", _numberOfPoints, 0, 1.0 );
1358 
1359  }
1360  }
1361 }
1362 
1363 void LCFIAIDAPlotProcessor::FillVertexChargePlots(LCEvent* pEvent, unsigned int jetNumber)
1364 {
1365  int jetType=FindTrueJetType( pEvent, jetNumber );
1366  if( jetType==0 ) return;
1367 
1368  //needs to be tidied up
1369  TypesafeCollection<lcio::ReconstructedParticle> jetCollection( pEvent, _JetCollectionName );
1370  ReconstructedParticle* pJet;
1371  pJet=jetCollection.getElementAt(jetNumber);
1372 
1373  for (unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection)
1374  {
1375  TypesafeCollection<lcio::LCFloatVec> tagCollection( pEvent, _FlavourTagCollectionNames[iTagCollection] );
1376  if( !tagCollection.is_valid() )
1377  {
1378  std::cerr << _FlavourTagCollectionNames[iTagCollection] << "is not valid, not filling vertex charge plots" << std::endl;
1379  std::cerr << "set MakeAdditionalPlots to false to remove this warning" << std::endl;
1380  }
1381  else
1382  {
1383  lcio::LCFloatVec* pJetFlavourTags=tagCollection.getElementAt( jetNumber );
1384  if( !pJetFlavourTags )
1385  {
1386  }
1387  else
1388  {
1389  const double* jetMomentum = pJet->getMomentum();
1390  double cosTheta = jetMomentum[2] / sqrt(pow(jetMomentum[0],2)+pow(jetMomentum[1],2)+pow(jetMomentum[2],2));
1391 
1392 
1393  if (iTagCollection == _myVertexChargeTagCollection) {
1394 
1395 
1396  double bTag= (*pJetFlavourTags)[_IndexOfForEachTag[iTagCollection]["BTag"]];
1397  double cTag= (*pJetFlavourTags)[_IndexOfForEachTag[iTagCollection]["CTag"]];
1398  //double cTagBBack= (*pJetFlavourTags)[_IndexOfForEachTag[iTagCollection]["BCTag"]];
1399  unsigned int NumVertices = FindNumVertex(pEvent, jetNumber, iTagCollection);
1400  int CQVtx = FindCQVtx(pEvent, jetNumber);
1401  int BQVtx = FindBQVtx(pEvent, jetNumber);
1402  int trueJetCharge = int(FindTrueJetHadronCharge(pEvent,jetNumber));
1403 
1404  std::string nvname = _VertexCatNames[ (NumVertices>=N_VERTEX_CATEGORIES) ? (N_VERTEX_CATEGORIES) : (NumVertices)];
1405 
1406  //vertex charge plots
1407  if( jetType==C_JET && cTag > _CTagNNCut) {
1408 
1409  int bin = _pCJetLeakageRate->coordToIndex(fabs(cosTheta));
1410 
1411  if (trueJetCharge==+2) _cJet_truePlus2++;
1412  if (trueJetCharge==+1) _cJet_truePlus++;
1413  if (trueJetCharge==0) _cJet_trueNeut++;
1414  if (trueJetCharge==-1) _cJet_trueMinus++;
1415  if (trueJetCharge==-2) _cJet_trueMinus2++;
1416 
1417  if (trueJetCharge==+2){
1418  if(CQVtx>0) _cJet_truePlus2_recoPlus++;
1419  if(CQVtx==0) _cJet_truePlus2_recoNeut++;
1420  if(CQVtx<0) _cJet_truePlus2_recoMinus++;
1421  }
1422  if (trueJetCharge==+1){
1423  if(CQVtx>0) _cJet_truePlus_recoPlus++;
1424  if(CQVtx==0) _cJet_truePlus_recoNeut++;
1425  if(CQVtx<0) _cJet_truePlus_recoMinus++;
1426  }
1427  if (trueJetCharge==0) {
1428  if(CQVtx>0) _cJet_trueNeut_recoPlus++;
1429  if(CQVtx==0) _cJet_trueNeut_recoNeut++;
1430  if(CQVtx<0) _cJet_trueNeut_recoMinus++;
1431  }
1432  if (trueJetCharge==-1) {
1433  if(CQVtx>0) _cJet_trueMinus_recoPlus++;
1434  if(CQVtx==0) _cJet_trueMinus_recoNeut++;
1435  if(CQVtx<0) _cJet_trueMinus_recoMinus++;
1436  }
1437  if (trueJetCharge==-2) {
1438  if(CQVtx>0) _cJet_trueMinus2_recoPlus++;
1439  if(CQVtx==0) _cJet_trueMinus2_recoNeut++;
1440  if(CQVtx<0) _cJet_trueMinus2_recoMinus++;
1441  }
1442 
1443  _pCJetVertexCharge->fill(CQVtx);
1444  _pCJetCharge2D->fill(trueJetCharge,CQVtx);
1445 
1446  if (trueJetCharge==+2) _cJet_truePlus2_angle[bin]++;
1447  if (trueJetCharge==+1) _cJet_truePlus_angle[bin]++;
1448  if (trueJetCharge==0) _cJet_trueNeut_angle[bin]++;
1449  if (trueJetCharge==-1) _cJet_trueMinus_angle[bin]++;
1450  if (trueJetCharge==-2) _cJet_trueMinus2_angle[bin]++;
1451 
1452  if (trueJetCharge==+2){
1453  if(CQVtx>0) _cJet_truePlus2_recoPlus_angle[bin]++;
1454  if(CQVtx==0) _cJet_truePlus2_recoNeut_angle[bin]++;
1455  if(CQVtx<0) _cJet_truePlus2_recoMinus_angle[bin]++;
1456  }
1457  if (trueJetCharge==+1){
1458  if(CQVtx>0) _cJet_truePlus_recoPlus_angle[bin]++;
1459  if(CQVtx==0) _cJet_truePlus_recoNeut_angle[bin]++;
1460  if(CQVtx<0) _cJet_truePlus_recoMinus_angle[bin]++;
1461  }
1462  if (trueJetCharge==0) {
1463  if(CQVtx>0) _cJet_trueNeut_recoPlus_angle[bin]++;
1464  if(CQVtx==0) _cJet_trueNeut_recoNeut_angle[bin]++;
1465  if(CQVtx<0) _cJet_trueNeut_recoMinus_angle[bin]++;
1466  }
1467  if (trueJetCharge==-1) {
1468  if(CQVtx>0) _cJet_trueMinus_recoPlus_angle[bin]++;
1469  if(CQVtx==0) _cJet_trueMinus_recoNeut_angle[bin]++;
1470  if(CQVtx<0) _cJet_trueMinus_recoMinus_angle[bin]++;
1471  }
1472  if (trueJetCharge==-2) {
1473  if(CQVtx>0) _cJet_trueMinus2_recoPlus_angle[bin]++;
1474  if(CQVtx==0) _cJet_trueMinus2_recoNeut_angle[bin]++;
1475  if(CQVtx<0) _cJet_trueMinus2_recoMinus_angle[bin]++;
1476  }
1477 
1478  } else if ( jetType==B_JET && bTag > _BTagNNCut) {
1479 
1480  int bin = _pBJetLeakageRate->coordToIndex(fabs(cosTheta));
1481 
1482  if (trueJetCharge==+2) _bJet_truePlus2++;
1483  if (trueJetCharge==+1) _bJet_truePlus++;
1484  if (trueJetCharge==0) _bJet_trueNeut++;
1485  if (trueJetCharge==-1) _bJet_trueMinus++;
1486  if (trueJetCharge==-2) _bJet_trueMinus2++;
1487 
1488  if (trueJetCharge==+2){
1489  if(BQVtx>0) _bJet_truePlus2_recoPlus++;
1490  if(BQVtx==0) _bJet_truePlus2_recoNeut++;
1491  if(BQVtx<0) _bJet_truePlus2_recoMinus++;
1492  }
1493  if (trueJetCharge==+1){
1494  if(BQVtx>0) _bJet_truePlus_recoPlus++;
1495  if(BQVtx==0) _bJet_truePlus_recoNeut++;
1496  if(BQVtx<0) _bJet_truePlus_recoMinus++;
1497  }
1498  if (trueJetCharge==0) {
1499  if(BQVtx>0) _bJet_trueNeut_recoPlus++;
1500  if(BQVtx==0) _bJet_trueNeut_recoNeut++;
1501  if(BQVtx<0) _bJet_trueNeut_recoMinus++;
1502  }
1503  if (trueJetCharge==-1) {
1504  if(BQVtx>0) _bJet_trueMinus_recoPlus++;
1505  if(BQVtx==0) _bJet_trueMinus_recoNeut++;
1506  if(BQVtx<0) _bJet_trueMinus_recoMinus++;
1507  }
1508  if (trueJetCharge==-2) {
1509  if(BQVtx>0) _bJet_trueMinus2_recoPlus++;
1510  if(BQVtx==0) _bJet_trueMinus2_recoNeut++;
1511  if(BQVtx<0) _bJet_trueMinus2_recoMinus++;
1512  }
1513 
1514 
1515  if (trueJetCharge==+2) _bJet_truePlus2_angle[bin]++;
1516  if (trueJetCharge==+1) _bJet_truePlus_angle[bin]++;
1517  if (trueJetCharge==0) _bJet_trueNeut_angle[bin]++;
1518  if (trueJetCharge==-1) _bJet_trueMinus_angle[bin]++;
1519  if (trueJetCharge==-2) _bJet_trueMinus2_angle[bin]++;
1520 
1521  if (trueJetCharge==+2){
1522  if(BQVtx>0) _bJet_truePlus2_recoPlus_angle[bin]++;
1523  if(BQVtx==0) _bJet_truePlus2_recoNeut_angle[bin]++;
1524  if(BQVtx<0) _bJet_truePlus2_recoMinus_angle[bin]++;
1525  }
1526  if (trueJetCharge==+1){
1527  if(BQVtx>0) _bJet_truePlus_recoPlus_angle[bin]++;
1528  if(BQVtx==0) _bJet_truePlus_recoNeut_angle[bin]++;
1529  if(BQVtx<0) _bJet_truePlus_recoMinus_angle[bin]++;
1530  }
1531  if (trueJetCharge==0) {
1532  if(BQVtx>0) _bJet_trueNeut_recoPlus_angle[bin]++;
1533  if(BQVtx==0) _bJet_trueNeut_recoNeut_angle[bin]++;
1534  if(BQVtx<0) _bJet_trueNeut_recoMinus_angle[bin]++;
1535  }
1536  if (trueJetCharge==-1) {
1537  if(BQVtx>0) _bJet_trueMinus_recoPlus_angle[bin]++;
1538  if(BQVtx==0) _bJet_trueMinus_recoNeut_angle[bin]++;
1539  if(BQVtx<0) _bJet_trueMinus_recoMinus_angle[bin]++;
1540  }
1541  if (trueJetCharge==-2) {
1542  if(BQVtx>0) _bJet_trueMinus2_recoPlus_angle[bin]++;
1543  if(BQVtx==0) _bJet_trueMinus2_recoNeut_angle[bin]++;
1544  if(BQVtx<0) _bJet_trueMinus2_recoMinus_angle[bin]++;
1545  }
1546 
1547  _pBJetVertexCharge->fill(BQVtx);
1548  _pBJetCharge2D->fill(trueJetCharge,BQVtx);
1549  }
1550  }
1551  }
1552  }
1553  }
1554 }
1555 
1556 
1557 int LCFIAIDAPlotProcessor::FindTrueJetPDGCode( LCEvent* pEvent, unsigned int jetNumber )
1558 {
1559  float pdgCode=0;
1560 
1561  TypesafeCollection<lcio::LCFloatVec> trueJetFlavourCollection( pEvent, _TrueJetFlavourColName);
1562  if( !trueJetFlavourCollection.is_valid() )
1563  {
1564  std::cerr << " In " << __FILE__ << "(" << __LINE__ << "): Collection " << _TrueJetFlavourColName << " is not valid " << std::endl;
1565  return 0; //can't do anything without this collection
1566  }
1567 
1568  lcio::LCFloatVec* pJetFlavour = trueJetFlavourCollection.getElementAt( jetNumber );
1569 
1570  if( !pJetFlavour )
1571  {
1572  std::cerr << __FILE__ << "(" << __LINE__ << "): Unable to get the LCFloatVec for jet " << jetNumber << " from the collection " << _TrueJetFlavourColName
1573  << " for event " << pEvent->getEventNumber() << " in run " << pEvent->getRunNumber() << "." << std::endl;
1574  }
1575  else
1576  {
1577  pdgCode = (*pJetFlavour)[_FlavourIndex["TrueJetPDGFlavour"]];
1578  }
1579 
1580  return int(pdgCode+0.01);//just to be safe
1581 }
1582 
1583 
1584 
1585 
1586 
1587 float LCFIAIDAPlotProcessor::FindTrueJetPartonCharge( LCEvent* pEvent, unsigned int jetNumber )
1588 {
1589  TypesafeCollection<lcio::LCFloatVec> trueJetFlavourCollection( pEvent, _TrueJetFlavourColName);
1590  if( !trueJetFlavourCollection.is_valid() )
1591  {
1592  std::cerr << " In " << __FILE__ << "(" << __LINE__ << "): Collection " << _TrueJetFlavourColName << " is not valid " << std::endl;
1593  return 0; //can't do anything without this collection
1594  }
1595 
1596  lcio::LCFloatVec* pJetFlavour = trueJetFlavourCollection.getElementAt( jetNumber );
1597  if( !pJetFlavour )
1598  {
1599  std::cerr << __FILE__ << "(" << __LINE__ << "): Unable to get the LCFloatVec for jet " << jetNumber << " from the collection " << _TrueJetFlavourColName
1600  << " for event " << pEvent->getEventNumber() << " in run " << pEvent->getRunNumber() << "." << std::endl;
1601  }
1602  else
1603  {
1604  return (*pJetFlavour)[_FlavourIndex["TruePartonCharge"]];
1605  }
1606  return 0.;
1607 }
1608 
1609 
1610 int LCFIAIDAPlotProcessor::FindTrueJetFlavour( LCEvent* pEvent, unsigned int jetNumber )
1611 {
1612  float jetFlavour=0.;
1613 
1614  TypesafeCollection<lcio::LCFloatVec> trueJetFlavourCollection( pEvent, _TrueJetFlavourColName);
1615  if( !trueJetFlavourCollection.is_valid() )
1616  {
1617  std::cerr << " In " << __FILE__ << "(" << __LINE__ << "): Collection " << _TrueJetFlavourColName << " is not valid " << std::endl;
1618  return 0; //can't do anything without this collection
1619  }
1620 
1621  lcio::LCFloatVec* pJetFlavour = trueJetFlavourCollection.getElementAt( jetNumber );
1622  if( !pJetFlavour )
1623  {
1624  std::cerr << __FILE__ << "(" << __LINE__ << "): Unable to get the LCFloatVec for jet " << jetNumber << " from the collection " << _TrueJetFlavourColName
1625  << " for event " << pEvent->getEventNumber() << " in run " << pEvent->getRunNumber() << "." << std::endl;
1626  }
1627  else
1628  {
1629  jetFlavour=(*pJetFlavour)[_FlavourIndex["TrueJetFlavour"]];
1630  }
1631  return int(jetFlavour+0.001);
1632 }
1633 
1634 float LCFIAIDAPlotProcessor::FindTrueJetHadronCharge( LCEvent* pEvent, unsigned int jetNumber )
1635 {
1636  TypesafeCollection<lcio::LCFloatVec> trueJetFlavourCollection( pEvent, _TrueJetFlavourColName);
1637  if( !trueJetFlavourCollection.is_valid() )
1638  {
1639  std::cerr << " In " << __FILE__ << "(" << __LINE__ << "): Collection " << _TrueJetFlavourColName << " is not valid " << std::endl;
1640  return 0; //can't do anything without this collection
1641  }
1642 
1643  lcio::LCFloatVec* pJetFlavour = trueJetFlavourCollection.getElementAt( jetNumber );
1644  if( !pJetFlavour )
1645  {
1646  std::cerr << __FILE__ << "(" << __LINE__ << "): Unable to get the LCFloatVec for jet " << jetNumber << " from the collection " << _TrueJetFlavourColName
1647  << " for event " << pEvent->getEventNumber() << " in run " << pEvent->getRunNumber() << "." << std::endl;
1648  }
1649  else
1650  {
1651  return (*pJetFlavour)[_FlavourIndex["TrueJetHadronCharge"]];
1652  }
1653  return 0.;
1654 }
1655 
1656 
1657 
1658 
1659 int LCFIAIDAPlotProcessor::FindTrueJetType( LCEvent* pEvent, unsigned int jetNumber )
1660 {
1661  return FindTrueJetFlavour(pEvent, jetNumber);
1662 }
1663 
1664 
1665 
1666 int LCFIAIDAPlotProcessor::FindNumVertex( LCEvent* pEvent, unsigned int jetNumber, unsigned int iInputsCollection)
1667 {
1668  TypesafeCollection<lcio::LCFloatVec> inputsCollection( pEvent, _FlavourTagInputsCollectionNames[iInputsCollection] );
1669 
1670  if( !inputsCollection.is_valid() ) {
1671 
1672  std::cerr << __FILE__ << "(" << __LINE__ << "): Unable to get the LCFloatVec from the collection " << _FlavourTagInputsCollectionNames[iInputsCollection]
1673  << " for event " << pEvent->getEventNumber() << " in run " << pEvent->getRunNumber() << "." << std::endl;
1674 
1675  return 0; //can't do anything without this collection
1676  }
1677  else
1678  {
1679  //Do stuff...
1680  lcio::LCFloatVec* pInputs=inputsCollection.getElementAt( jetNumber );
1681 
1682  if( !pInputs )
1683  {
1684  std::cerr << __FILE__ << "(" << __LINE__ << "): Unable to get the LCFloatVec for jet " << jetNumber << " from the collection " << _FlavourTagInputsCollectionNames[iInputsCollection]
1685  << " for event " << pEvent->getEventNumber() << " in run " << pEvent->getRunNumber() << "." << std::endl;
1686  }
1687  else
1688  {
1689  return int((*pInputs)[_InputsIndex[iInputsCollection]["NumVertices"]]);
1690  }
1691 
1692  }
1693  return 0;
1694 }
1695 
1696 int LCFIAIDAPlotProcessor::FindBQVtx( LCEvent* pEvent, unsigned int jetNumber)
1697 {
1698 
1699  TypesafeCollection<lcio::LCFloatVec> inputsCollection( pEvent, _BVertexChargeCollection);
1700 
1701  if( !inputsCollection.is_valid() ) {
1702 
1703  std::cerr << "In " << __FILE__ << "(" << __LINE__ << "): Cannot find collection " << _BVertexChargeCollection << " for event " << pEvent->getEventNumber() << " in run " << pEvent->getRunNumber() << " BQVtx will be invalid" << std::endl;
1704  return -99;
1705 
1706  } else { //inputsCollection.is_valid()
1707 
1708  float bqvtx;
1709  lcio::LCFloatVec* pBVtxChargeVector =inputsCollection.getElementAt( jetNumber );
1710 
1711  //bool evaluation is done left to right...
1712  if( pBVtxChargeVector && pBVtxChargeVector->size() == 1) {
1713 
1714  bqvtx = pBVtxChargeVector->back();
1715 
1716  } else {
1717 
1718  std::cerr << "In " << __FILE__ << "(" << __LINE__ << "): Cannot find collection element in " << _BVertexChargeCollection << " for event " << pEvent->getEventNumber() << " in run " << pEvent->getRunNumber() << " corresponding to jet number: " << jetNumber << " BQVtx will be invalid" << std::endl;
1719  return -99;
1720  }
1721 
1722  return int(bqvtx);
1723  }
1724 
1725  //should never get here
1726  return -99;
1727 }
1728 
1729 
1730 int LCFIAIDAPlotProcessor::FindCQVtx( LCEvent* pEvent, unsigned int jetNumber)
1731 {
1732 
1733  TypesafeCollection<lcio::LCFloatVec> inputsCollection( pEvent, _CVertexChargeCollection);
1734 
1735  if( !inputsCollection.is_valid() ) {
1736 
1737  std::cerr << "In " << __FILE__ << "(" << __LINE__ << "): Cannot find collection " << _CVertexChargeCollection << " for event " << pEvent->getEventNumber() << " in run " << pEvent->getRunNumber() << " CQVtx will be invalid" << std::endl;
1738  return -99;
1739 
1740  } else { //inputsCollection.is_valid()
1741 
1742  float cqvtx;
1743  lcio::LCFloatVec* pCVtxChargeVector =inputsCollection.getElementAt( jetNumber );
1744 
1745  //bool evaluation is done left to right...
1746  if( pCVtxChargeVector && pCVtxChargeVector->size() == 1) {
1747 
1748  cqvtx = pCVtxChargeVector->back();
1749 
1750  } else {
1751 
1752  std::cerr << "In " << __FILE__ << "(" << __LINE__ << "): Cannot find collection element in " << _CVertexChargeCollection << " for event " << pEvent->getEventNumber() << " in run " << pEvent->getRunNumber() << " corresponding to jet number: " << jetNumber << " CQVtx will be invalid" << std::endl;
1753  return -99;
1754  }
1755 
1756  return int(cqvtx);
1757  }
1758  //should never get here
1759  return -99;
1760 }
1761 
1762 
1763 AIDA::IHistogram1D* LCFIAIDAPlotProcessor::CreateIntegralHistogram(const AIDA::IHistogram1D* pNN, AIDA::IHistogram1D* pIntegral)
1764 {
1765  //the errors on these entries are wrong...
1766 
1767  const int numberOfBins=pNN->axis().bins();
1768 
1769  double integral=pNN->binHeight(AIDA::IAxis::OVERFLOW_BIN);
1770  pIntegral->fill(pNN->axis().binLowerEdge(AIDA::IAxis::OVERFLOW_BIN)+pNN->axis().binWidth(numberOfBins-1)/2.,integral);
1771 
1772  for( int binNumber=numberOfBins-1; binNumber>=0; --binNumber )
1773  {
1774  integral+= pNN->binHeight( binNumber );
1775  pIntegral->fill( pNN->axis().binLowerEdge(binNumber)+pNN->axis().binWidth(binNumber)/2.,integral);
1776  }
1777 
1778  integral+= pNN->binHeight(AIDA::IAxis::UNDERFLOW_BIN);
1779  pIntegral->fill(pNN->axis().binUpperEdge(AIDA::IAxis::UNDERFLOW_BIN)-pNN->axis().binWidth(0)/2.,integral);
1780 
1781  return pIntegral;
1782 }
1783 
1784 
1785 
1786 void LCFIAIDAPlotProcessor::CalculateTagEfficiency(const AIDA::IHistogram1D* pSignal, std::vector<double>& tagEfficiency, std::vector<double>& tagEfficiencyError)
1787 {
1788  //makes an effiency plot from one histogram - assuming the last bin contains all the events (eg NN distribution)
1789  double totalSignal=pSignal->sumBinHeights();
1790  double signalPassedCut=0;
1791  const int numberOfBins=pSignal->axis().bins();
1792 
1793  tagEfficiency.clear();
1794  tagEfficiencyError.clear();
1795 
1796 
1797  for( int binNumber=numberOfBins-1; binNumber>=0; --binNumber)
1798  {
1799  signalPassedCut+=pSignal->binHeight( binNumber );
1800 
1801  double efficiency = signalPassedCut/totalSignal;
1802 
1803  double efficiencyError = efficiency * (1. - efficiency) / totalSignal;
1804  if (efficiencyError>0) efficiencyError = sqrt(efficiencyError);
1805 
1806  tagEfficiency.push_back(efficiency);
1807  tagEfficiencyError.push_back(efficiencyError);
1808 
1809  }
1810  return;
1811 }
1812 
1813 void LCFIAIDAPlotProcessor::CalculateTagPurity(const AIDA::IHistogram1D* pSignal, const AIDA::IHistogram1D* pBackground, std::vector<double>& tagPurity, std::vector<double>& tagPurityError)
1814 {
1815  const int numberOfBins=pSignal->axis().bins();
1816 
1817  double signalPassedCut=0;
1818  double backgroundPassedCut=0;
1819 
1820  tagPurity.clear();
1821  tagPurityError.clear();
1822 
1823 
1824  for (int binNumber = numberOfBins-1; binNumber >= 0 ; --binNumber)
1825  {
1826 
1827  signalPassedCut+=pSignal->binHeight( binNumber );
1828  backgroundPassedCut+=pBackground->binHeight( binNumber );
1829 
1830  double purity = ((signalPassedCut+backgroundPassedCut)>0) ? signalPassedCut/(signalPassedCut+backgroundPassedCut) : 1.;
1831  double purityError = ((signalPassedCut+backgroundPassedCut)>0) ? purity * (1. - purity) / (signalPassedCut+backgroundPassedCut) : 0.;
1832  if (purityError>0) purityError = sqrt(purityError);
1833  else purityError = -sqrt(-purityError);
1834 
1835  tagPurity.push_back(purity);
1836  tagPurityError.push_back(purityError);
1837 
1838  }
1839 
1840  return;
1841 }
1842 
1843 
1844 
1845 #ifdef USING_JAIDA
1846 AIDA::IDataPointSet* LCFIAIDAPlotProcessor::CreateEfficiencyPlot(const AIDA::IHistogram1D* pSignal, AIDA::IDataPointSet* pDataPointSet)
1847 {
1848  //makes an effiency plot from one histogram - assuming the last bin contains all the events (eg NN distribution)
1849  double totalSignal=pSignal->sumBinHeights();
1850  double signalPassedCut=0;
1851 
1852  const int numberOfBins=pSignal->axis().bins();
1853  int iPoint=0;
1854 
1855  for( int binNumber=numberOfBins-1; binNumber>=0; --binNumber, iPoint++ )
1856  {
1857  signalPassedCut+=pSignal->binHeight( binNumber );
1858 
1859  double efficiency = signalPassedCut/totalSignal;
1860 
1861  double efficiencyError = efficiency * (1. - efficiency) / totalSignal;
1862  if (efficiencyError>0) efficiencyError = sqrt(efficiencyError);
1863 
1864 
1865  pDataPointSet->addPoint();
1866  pDataPointSet->point(iPoint)->coordinate(0)->setValue(pSignal->axis().binLowerEdge(binNumber)+pSignal->axis().binWidth(binNumber)/2.);
1867  pDataPointSet->point(iPoint)->coordinate(1)->setValue( efficiency );
1868  pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(efficiencyError);
1869  pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(efficiencyError);
1870  }
1871 
1872  return pDataPointSet;
1873 }
1874 
1875 
1876 AIDA::IDataPointSet* LCFIAIDAPlotProcessor::CreateEfficiencyPlot2(const AIDA::IHistogram1D* pAllEvents, const AIDA::IHistogram1D *pPassEvents, AIDA::IDataPointSet* pDataPointSet)
1877 {
1878  //I'm returning a datapoint set as then I can set the errors to the correct values - not possible for a histogram
1879  //make a "traditional" efficiency plot, divides the 2nd histogram by the 1st one
1880 
1881  const int numberOfBins=pAllEvents->axis().bins();
1882  if (pPassEvents->axis().bins() != numberOfBins) return pDataPointSet;
1883 
1884  int iPoint=0;
1885 
1886  for( int binNumber=numberOfBins-1; binNumber>=0; --binNumber, iPoint++ )
1887  {
1888  double all = pAllEvents->binHeight(binNumber);
1889  double pass = pPassEvents->binHeight(binNumber);
1890 
1891  double efficiency = (all>0) ? pass/all : 0.;
1892 
1893  double efficiencyError = (all>0) ? efficiency * (1. - efficiency) / all : 0.;
1894  if (efficiencyError>0.) efficiencyError = sqrt(efficiencyError);
1895 
1896  pDataPointSet->addPoint();
1897  pDataPointSet->point(iPoint)->coordinate(0)->setValue(pAllEvents->axis().binLowerEdge(binNumber)+pAllEvents->axis().binWidth(binNumber)/2.);
1898  pDataPointSet->point(iPoint)->coordinate(1)->setValue( efficiency );
1899  pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(efficiencyError);
1900  pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(efficiencyError);
1901  }
1902 
1903  return pDataPointSet;
1904 }
1905 
1906 
1907 
1908 AIDA::IDataPointSet* LCFIAIDAPlotProcessor::CreateIntegralPlot(const AIDA::IHistogram1D* pNN, AIDA::IDataPointSet* pDataPointSet )
1909 {
1910  //it might make more sense to make this a histogram, but then the error on the entries are wrong
1911 
1912  const int numberOfBins=pNN->axis().bins();
1913 
1914  double integral=0;
1915 
1916  for (int binNumber = 0; binNumber < numberOfBins ; binNumber++ )
1917  {
1918 
1919  integral+= pNN->binHeight( binNumber );
1920  pDataPointSet->addPoint();
1921  pDataPointSet->point(binNumber)->coordinate(0)->setValue(pNN->axis().binLowerEdge(binNumber)+pNN->axis().binWidth(binNumber)/2.);
1922  pDataPointSet->point(binNumber)->coordinate(1)->setValue(integral);
1923  pDataPointSet->point(binNumber)->coordinate(1)->setErrorPlus(sqrt(integral));
1924  pDataPointSet->point(binNumber)->coordinate(1)->setErrorMinus(sqrt(integral));
1925  }
1926  return pDataPointSet;
1927 }
1928 
1929 AIDA::IDataPointSet* LCFIAIDAPlotProcessor::CreatePurityPlot(const AIDA::IHistogram1D* pSignal, const AIDA::IHistogram1D* pBackground, AIDA::IDataPointSet* pDataPointSet)
1930 {
1931  const int numberOfBins=pSignal->axis().bins();
1932  int iPoint=0;
1933 
1934  double signalPassedCut=0;
1935  double backgroundPassedCut=0;
1936 
1937  for (int binNumber = numberOfBins-1; binNumber >= 0 ; --binNumber, iPoint++ )
1938  {
1939 
1940  signalPassedCut+=pSignal->binHeight( binNumber );
1941  backgroundPassedCut+=pBackground->binHeight( binNumber );
1942 
1943  double purity = signalPassedCut/(signalPassedCut+backgroundPassedCut);
1944  double purityError = purity * (1. - purity) / (signalPassedCut+backgroundPassedCut);
1945  if (purityError>0) purityError = sqrt(purityError);
1946 
1947 
1948  pDataPointSet->addPoint();
1949  pDataPointSet->point(iPoint)->coordinate(0)->setValue(pSignal->axis().binLowerEdge(binNumber)+pSignal->axis().binWidth(binNumber)/2.);
1950  pDataPointSet->point(iPoint)->coordinate(1)->setValue(purity);
1951  pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(purityError);
1952  pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(purityError);
1953 
1954 
1955  }
1956 
1957  return pDataPointSet;
1958 }
1959 
1960 AIDA::IDataPointSet* LCFIAIDAPlotProcessor::CreateLeakageRatePlot(const AIDA::IHistogram1D* pBackground, AIDA::IDataPointSet* pDataPointSet)
1961 {
1962 
1963  double totalBackground = pBackground->sumBinHeights();
1964  double backgroundPassedCut=0;
1965 
1966  const int numberOfBins=pBackground->axis().bins();
1967  int iPoint=0;
1968 
1969  for( int binNumber=numberOfBins-1; binNumber>=0; --binNumber , iPoint++ )
1970  {
1971 
1972  backgroundPassedCut+=pBackground->binHeight( binNumber );
1973 
1974  double leakageRate = backgroundPassedCut/totalBackground;
1975 
1976  double leakageRateError = leakageRate * (1. - leakageRate) / totalBackground;
1977  if (leakageRateError>0) leakageRateError = sqrt(leakageRateError);
1978 
1979 
1980  pDataPointSet->addPoint();
1981  pDataPointSet->point(iPoint)->coordinate(0)->setValue(pBackground->axis().binLowerEdge(binNumber)+pBackground->axis().binWidth(binNumber)/2.);
1982  pDataPointSet->point(iPoint)->coordinate(1)->setValue(leakageRate);
1983  pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(leakageRateError);
1984  pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(leakageRateError);
1985  }
1986 
1987  return pDataPointSet;
1988 }
1989 
1990 
1991 
1992 AIDA::IDataPointSet* LCFIAIDAPlotProcessor::CreateXYPlot(const AIDA::IDataPointSet* pDataPointSet0, const AIDA::IDataPointSet* pDataPointSet1, AIDA::IDataPointSet* xyPointSet, const int dim0, const int dim1 )
1993 {
1994 
1995  //need to do some comparision here
1996  if (pDataPointSet0->size() == pDataPointSet1->size()) {
1997 
1998  for (int iPoint = 0 ; iPoint != pDataPointSet1->size(); iPoint++)
1999  {
2000  xyPointSet->addPoint();
2001  xyPointSet->point(iPoint)->coordinate(0)->setValue(pDataPointSet0->point(iPoint)->coordinate(dim0)->value());
2002  xyPointSet->point(iPoint)->coordinate(1)->setValue(pDataPointSet1->point(iPoint)->coordinate(dim1)->value());
2003  xyPointSet->point(iPoint)->coordinate(0)->setErrorPlus(pDataPointSet0->point(iPoint)->coordinate(dim0)->errorPlus());
2004  xyPointSet->point(iPoint)->coordinate(1)->setErrorPlus(pDataPointSet1->point(iPoint)->coordinate(dim1)->errorPlus());
2005  xyPointSet->point(iPoint)->coordinate(0)->setErrorMinus(pDataPointSet0->point(iPoint)->coordinate(dim0)->errorMinus());
2006  xyPointSet->point(iPoint)->coordinate(1)->setErrorMinus(pDataPointSet1->point(iPoint)->coordinate(dim1)->errorMinus());
2007  }
2008  } else {
2009 
2010  //some error message here!
2011  }
2012  return xyPointSet;
2013 }
2014 
2015 void LCFIAIDAPlotProcessor::CreateVertexChargeLeakagePlot(AIDA::IDataPointSet* pBJetVtxChargeDPS, AIDA::IDataPointSet* pCJetVtxChargeDPS)
2016 {
2017  for (int j = 0 ; j < N_JETANGLE_BINS ; j++) {
2018 
2019  double c_numerator = _cJet_truePlus2_recoPlus_angle[j]+_cJet_truePlus_recoPlus_angle[j]
2020  +_cJet_trueMinus2_recoMinus_angle[j]+_cJet_trueMinus_recoMinus_angle[j];
2021  double c_domininator = _cJet_truePlus2_angle[j] +_cJet_truePlus_angle[j]
2022  +_cJet_trueMinus2_angle[j] +_cJet_trueMinus_angle[j];
2023 
2024  double b_numerator = _bJet_truePlus2_recoPlus_angle[j]+_bJet_truePlus_recoPlus_angle[j]
2025  +_bJet_trueMinus2_recoMinus_angle[j]+_bJet_trueMinus_recoMinus_angle[j];
2026  double b_domininator = _bJet_truePlus2_angle[j] +_bJet_truePlus_angle[j]
2027  +_bJet_trueMinus2_angle[j] +_bJet_trueMinus_angle[j];
2028 
2029  double b_leakage = 1. - b_numerator/b_domininator;
2030  double c_leakage = 1. - c_numerator/c_domininator;
2031 
2032  double b_leakage_error = sqrt( b_leakage * (1. - b_leakage) / b_domininator);
2033  double c_leakage_error = sqrt( c_leakage * (1. - c_leakage) / c_domininator);
2034 
2035 
2036  pBJetVtxChargeDPS->addPoint();
2037  pBJetVtxChargeDPS->point(j)->coordinate(0)->setValue(_pBJetLeakageRate->axis().binLowerEdge(j)+_pBJetLeakageRate->axis().binWidth(j)/2.);
2038  pBJetVtxChargeDPS->point(j)->coordinate(1)->setValue(b_leakage);
2039  pBJetVtxChargeDPS->point(j)->coordinate(1)->setErrorPlus(b_leakage_error);
2040  pBJetVtxChargeDPS->point(j)->coordinate(1)->setErrorMinus(b_leakage_error);
2041 
2042 
2043  pCJetVtxChargeDPS->addPoint();
2044  pCJetVtxChargeDPS->point(j)->coordinate(0)->setValue(_pCJetLeakageRate->axis().binLowerEdge(j)+_pCJetLeakageRate->axis().binWidth(j)/2.);
2045  pCJetVtxChargeDPS->point(j)->coordinate(1)->setValue(c_leakage);
2046  pCJetVtxChargeDPS->point(j)->coordinate(1)->setErrorPlus(c_leakage_error);
2047  pCJetVtxChargeDPS->point(j)->coordinate(1)->setErrorMinus(c_leakage_error);
2048 
2049  _pCJetLeakageRate->fill(_pCJetLeakageRate->axis().binLowerEdge(j)+_pCJetLeakageRate->axis().binWidth(j)/2.,c_leakage);
2050  _pBJetLeakageRate->fill(_pBJetLeakageRate->axis().binLowerEdge(j)+_pBJetLeakageRate->axis().binWidth(j)/2.,b_leakage);
2051  }
2052 }
2053 #endif
2054 
2055 
2056 
2057 void LCFIAIDAPlotProcessor::CreateVertexChargeLeakagePlot()
2058 {
2059 
2060  for (int j = 0 ; j < N_JETANGLE_BINS ; j++) {
2061 
2062  double c_numerator = _cJet_truePlus2_recoPlus_angle[j]+_cJet_truePlus_recoPlus_angle[j]
2063  +_cJet_trueMinus2_recoMinus_angle[j]+_cJet_trueMinus_recoMinus_angle[j];
2064  double c_domininator = _cJet_truePlus2_angle[j] +_cJet_truePlus_angle[j]
2065  +_cJet_trueMinus2_angle[j] +_cJet_trueMinus_angle[j];
2066 
2067  double b_numerator = _bJet_truePlus2_recoPlus_angle[j]+_bJet_truePlus_recoPlus_angle[j]
2068  +_bJet_trueMinus2_recoMinus_angle[j]+_bJet_trueMinus_recoMinus_angle[j];
2069  double b_domininator = _bJet_truePlus2_angle[j] +_bJet_truePlus_angle[j]
2070  +_bJet_trueMinus2_angle[j] +_bJet_trueMinus_angle[j];
2071 
2072  double b_leakage = 1. - b_numerator/b_domininator;
2073  double c_leakage = 1. - c_numerator/c_domininator;
2074 
2075  //double b_leakage_error = sqrt( b_leakage * (1. - b_leakage) / b_domininator);
2076  //double c_leakage_error = sqrt( c_leakage * (1. - c_leakage) / c_domininator);
2077 
2078  _pCJetLeakageRate->fill(_pCJetLeakageRate->axis().binLowerEdge(j)+_pCJetLeakageRate->axis().binWidth(j)/2.,c_leakage);
2079  _pBJetLeakageRate->fill(_pBJetLeakageRate->axis().binLowerEdge(j)+_pBJetLeakageRate->axis().binWidth(j)/2.,b_leakage);
2080 
2081  }
2082 }
2083 
2084 
2085 void LCFIAIDAPlotProcessor::FillVertexPlots(LCEvent* pEvent, unsigned int jetNumber)
2086 {
2087  int jetType=FindTrueJetType( pEvent, jetNumber );
2088 
2089  if( jetType==0 ) return;
2090 
2091  std::vector<double> mcDecayLengthVector;
2092  std::vector<double> bMCDecayLengthVector;
2093  std::vector<double> cMCDecayLengthVector;
2094  double b_mc_decay_length, c_mc_decay_length;
2095 
2096  FindTrueJetDecayLength(pEvent, jetNumber, mcDecayLengthVector, bMCDecayLengthVector, cMCDecayLengthVector);
2097  FindTrueJetDecayLength2(pEvent, jetNumber, b_mc_decay_length, c_mc_decay_length);
2098 
2099 
2100  for (unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection) {
2101 
2102  unsigned int NumVertices = FindNumVertex(pEvent, jetNumber, iTagCollection);
2103 
2104  if( jetType==B_JET ) {
2105  _pBDecayLengthAll[iTagCollection]->fill(b_mc_decay_length);
2106  if (NumVertices > 1) _pBDecayLengthTwoVertices[iTagCollection]->fill(b_mc_decay_length);
2107  } else if( jetType==C_JET ) {
2108  _pCDecayLengthAll[iTagCollection]->fill(c_mc_decay_length);
2109  if (NumVertices > 1) _pCDecayLengthTwoVertices[iTagCollection]->fill(c_mc_decay_length);
2110  }
2111 
2112  }
2113 
2114 
2115  TypesafeCollection<lcio::ReconstructedParticle> zvresjetCollection( pEvent, _ZVRESSelectedJetsCollection );
2116  TypesafeCollection<lcio::ReconstructedParticle> zvresdcCollection( pEvent, _ZVRESDecayChainCollection );
2117 
2118  if (!zvresjetCollection.is_valid() || !zvresdcCollection.is_valid() ) {
2119 
2120  std::cerr << "Input Collections " << _ZVRESSelectedJetsCollection << " and/or " << _ZVRESDecayChainCollection << " not valid. Will not fill vertex plots." <<std::endl;
2121  std::cerr << "Set MakeAdditionalPlots to false to remove this warning" << std::endl;
2122  }
2123  else
2124  {
2125  // zvresjetCollection and zvresdcCollection is valid
2126 
2127 
2128 
2129 
2130  for (int ijet = 0 ; ijet < zvresdcCollection.getNumberOfElements(); ijet++) {
2131 
2132 
2133  int jetType = FindTrueJetType( pEvent, ijet );
2134 
2135  lcio::ReconstructedParticle* pJet=zvresdcCollection.getElementAt(ijet);
2136 
2137  std::vector<Vertex*> myVertexVector;
2138  std::pair<Vertex*,Vertex*> myVertexPair;
2139  std::vector< std::pair<Vertex*,Vertex*> > myVertexPairVector;
2140 
2141  std::vector<lcio::ReconstructedParticle*> decayChainPartColl = pJet->getParticles();
2142  for (unsigned int ipart = 0 ; ipart < decayChainPartColl.size(); ipart++) {
2143 
2144  lcio::ReconstructedParticle* decayChainPart = dynamic_cast< lcio::ReconstructedParticle*> (decayChainPartColl[ipart]);
2145 
2146  //want to make an exclusive list of vertices
2147 
2148  bool alreadyContainsStart = false;
2149  bool alreadyContainsEnd = false;
2150 
2151  for (std::vector<Vertex*>::const_iterator iVert=myVertexVector.begin(); iVert!=myVertexVector.end(); iVert++) {
2152  if (decayChainPart->getStartVertex() == *iVert) alreadyContainsStart = true;
2153  if (decayChainPart->getEndVertex() == *iVert) alreadyContainsEnd = true;
2154  }
2155 
2156  if (!alreadyContainsStart && decayChainPart->getStartVertex()!=0) myVertexVector.push_back(decayChainPart->getStartVertex());
2157  if (!alreadyContainsEnd && decayChainPart->getEndVertex()!=0) myVertexVector.push_back(decayChainPart->getEndVertex());
2158 
2159 
2160  myVertexPair.first = decayChainPart->getStartVertex();
2161  myVertexPair.second = decayChainPart->getEndVertex();
2162  bool pairAlreadyContained = false;
2163 
2164  for (std::vector< std::pair<Vertex*,Vertex*> >::const_iterator iter = myVertexPairVector.begin(); iter != myVertexPairVector.end(); iter++) {
2165  if ((*iter).first == myVertexPair.first && (*iter).second == myVertexPair.second) pairAlreadyContained = true;
2166  }
2167 
2168  if (myVertexPair.first == 0 || myVertexPair.second == 0) pairAlreadyContained = true;
2169 
2170  if (!pairAlreadyContained) myVertexPairVector.push_back(myVertexPair);
2171  }
2172 
2173  //want to find distance between Primary Vertex and next vertex (aka seconday vertex); then the secondary vertex and the tertiary vertex
2174 
2175  float reconstructedSecondaryDecayLength = -1.;
2176  float reconstructedSecTerDecayLength = - 1.;
2177  std::vector< Vertex* > secondaryVertexVector;
2178  std::vector< float > secondaryDecayLengthVector;
2179 
2180  for (std::vector< std::pair<Vertex*,Vertex*> >::const_iterator iter = myVertexPairVector.begin();
2181  iter != myVertexPairVector.end(); iter++) {
2182 
2183  const Vertex* startV = (*iter).first;
2184 
2185  if (startV->isPrimary()) {
2186 
2187  const Vertex* endV = (*iter).second;
2188 
2189  const float* startPos = startV->getPosition();
2190  const float* endPos = endV->getPosition();
2191 
2192  secondaryVertexVector.push_back((*iter).second);
2193 
2194  reconstructedSecondaryDecayLength = CalculateDistance(startPos,endPos);
2195 
2196  _reconstructedSecondaryDecayLength->fill(reconstructedSecondaryDecayLength);
2197 
2198 
2199  if( jetType==B_JET ) {
2200  _recoDecayLengthBJet->fill(reconstructedSecondaryDecayLength);
2201 
2202  } else if( jetType==C_JET ) {
2203  _recoDecayLengthCJet->fill(reconstructedSecondaryDecayLength);
2204 
2205  } else {
2206  _recoDecayLengthLightJet->fill(reconstructedSecondaryDecayLength);
2207 
2208  }
2209  }
2210 
2211  }
2212 
2213 
2214  for (std::vector< std::pair<Vertex*,Vertex*> >::const_iterator iter = myVertexPairVector.begin();
2215  iter != myVertexPairVector.end(); iter++) {
2216 
2217  const Vertex* startV = (*iter).first;
2218 
2219  bool isSecondary=false;
2220  for (std::vector< Vertex* >::const_iterator iter2 = secondaryVertexVector.begin(); iter2 !=secondaryVertexVector.end(); iter2++)
2221  {
2222  if (startV == (*iter2)) isSecondary=true;
2223  }
2224 
2225  if (isSecondary) {
2226  const Vertex* endV = (*iter).second;
2227 
2228  const float* startPos = startV->getPosition();
2229  const float* endPos = endV->getPosition();
2230  reconstructedSecTerDecayLength = CalculateDistance(startPos,endPos);
2231  secondaryDecayLengthVector.push_back(reconstructedSecTerDecayLength);
2232 
2233  _reconstructedSecTerDecayLength->fill(reconstructedSecTerDecayLength);
2234 
2235  if (jetType==B_JET) _recoDecayLengthBCJet->fill(reconstructedSecTerDecayLength);
2236 
2237  }
2238  }
2239 
2240 
2241  int nRecoVertices = myVertexVector.size();
2242 
2243  if( jetType==B_JET ) {
2244  _nVerticesBJet->fill(nRecoVertices);
2245  _decayLengthBJet2D->fill(b_mc_decay_length,reconstructedSecondaryDecayLength);
2246  _decayLengthBJetCloud2D->fill(b_mc_decay_length,reconstructedSecondaryDecayLength);
2247 
2248  _decayLengthCJet2D->fill(c_mc_decay_length,reconstructedSecTerDecayLength);
2249  _decayLengthCJetCloud2D->fill(c_mc_decay_length,reconstructedSecTerDecayLength);
2250 
2251  _decayLengthBJetTrue->fill(b_mc_decay_length);
2252  _decayLengthBCJetTrue->fill(c_mc_decay_length);
2253 
2254  } else if( jetType==C_JET ) {
2255  _nVerticesCJet->fill(nRecoVertices);
2256  _decayLengthBJetTrue->fill(c_mc_decay_length);
2257 
2258  } else {
2259  _nVerticesLightJet->fill(nRecoVertices);
2260  }
2261 
2262 
2263  }
2264  }
2265 
2266  // Plots of the vertex quantities, if there is a vertex collection
2267  TypesafeCollection<lcio::Vertex> vertexCol(pEvent, _VertexColName);
2268 
2269  if (vertexCol.is_valid()) {
2270 
2271  float primaryVertexPostion[] = {0.,0.,0.};
2272 
2273  for (int ii=0 ; ii < vertexCol.getNumberOfElements() ; ii++){
2274 
2275  Vertex* myVertex = dynamic_cast<Vertex*>(vertexCol.getElementAt(ii));
2276 
2277  const float* vertexPosition = myVertex->getPosition();
2278  if (myVertex->isPrimary()) {
2279  primaryVertexPostion[0] = vertexPosition[0];
2280  primaryVertexPostion[1] = vertexPosition[1];
2281  primaryVertexPostion[2] = vertexPosition[2];
2282  }
2283  }
2284 
2285  for (int ii=0 ; ii < vertexCol.getNumberOfElements() ; ii++){
2286 
2287  Vertex* myVertex = dynamic_cast<Vertex*>(vertexCol.getElementAt(ii));
2288  const float* vertexPosition = myVertex->getPosition();
2289 
2290  double px = double(vertexPosition[0]);
2291  double py = double(vertexPosition[1]);
2292  double pz = double(vertexPosition[2]);
2293  double ex = sqrt((double)myVertex->getCovMatrix()[0]);
2294  double ey = sqrt((double)myVertex->getCovMatrix()[2]);
2295  double ez = sqrt((double)myVertex->getCovMatrix()[5]);
2296 
2297 
2298  if (! myVertex->isPrimary() ) {
2299 
2300  double distanceIP = sqrt(pow(px-primaryVertexPostion[0],2)+pow(py-primaryVertexPostion[1],2)+pow(pz-primaryVertexPostion[1],2));
2301 
2302  _pVertexDistanceFromIP->fill(distanceIP);
2303  _pVertexPositionX->fill(px);
2304  _pVertexPositionY->fill(py);
2305  _pVertexPositionZ->fill(pz);
2306 
2307  } else { //it is the primary vertex
2308 
2309  _pPrimaryVertexPositionX->fill(px);
2310  _pPrimaryVertexPositionY->fill(py);
2311  _pPrimaryVertexPositionZ->fill(pz);
2312 
2313  _pPrimaryVertexPullX->fill(px/ex);
2314  _pPrimaryVertexPullY->fill(py/ey);
2315  _pPrimaryVertexPullZ->fill(pz/ez);
2316  }
2317  }
2318  } else {
2319 
2320  std::cerr << "Input Collections " << _VertexColName << " is not valid. Not filling all vertex plots." <<std::endl;
2321  std::cerr << "Set MakeAdditionalPlots to false to remove this warning" << std::endl;
2322 
2323  }
2324 
2325 
2326 
2327 }
2328 
2329 void LCFIAIDAPlotProcessor::FillZVRESTable(LCEvent* pEvent)
2330 {
2331  TypesafeCollection<lcio::ReconstructedParticle> ftjetCollection( pEvent, _JetCollectionName );
2332  TypesafeCollection<lcio::ReconstructedParticle> zvresjetCollection( pEvent, _ZVRESSelectedJetsCollection );
2333  TypesafeCollection<lcio::ReconstructedParticle> zvresdcrptrackCollection( pEvent, _ZVRESDecayChainRPTracksCollection);
2334  TypesafeCollection<lcio::ReconstructedParticle> zvresdcCollection( pEvent, _ZVRESDecayChainCollection );
2335 
2336  LCCollection* relationTrackCollection = pEvent->getCollection(_TrueTracksToMCPCollection);
2337 
2338  UTIL::LCRelationNavigator* MCRelationNavigator = (relationTrackCollection!=NULL) ? new LCRelationNavigator(relationTrackCollection) : NULL;
2339 
2340 
2341  if (MCRelationNavigator==NULL)
2342  std::cerr << __FILE__ << "(" << __LINE__ << "), run number: " << pEvent->getRunNumber() << ", event number: " << pEvent->getEventNumber()<< " LCCollection " << _TrueTracksToMCPCollection << " is not valid." << std::endl;
2343 
2344  if (!ftjetCollection.is_valid())
2345  std::cerr << __FILE__ << "(" << __LINE__ << "), run number: " << pEvent->getRunNumber() << ", event number: " << pEvent->getEventNumber()<< " ReconstructedParticle " << _JetCollectionName << " is not valid." << std::endl;
2346  if (!zvresjetCollection.is_valid())
2347  std::cerr << __FILE__ << "(" << __LINE__ << "), run number: " << pEvent->getRunNumber() << ", event number: " << pEvent->getEventNumber()<< " ReconstructedParticle " << _ZVRESSelectedJetsCollection << " is not valid." << std::endl;
2348  if (!zvresdcrptrackCollection.is_valid())
2349  std::cerr << __FILE__ << "(" << __LINE__ << "), run number: " << pEvent->getRunNumber() << ", event number: " << pEvent->getEventNumber()<< " ReconstructedParticle " << _ZVRESDecayChainRPTracksCollection << " is not valid." << std::endl;
2350  if (!zvresdcCollection.is_valid())
2351  std::cerr << __FILE__ << "(" << __LINE__ << "), run number: " << pEvent->getRunNumber() << ", event number: " << pEvent->getEventNumber()<< " ReconstructedParticle " << _ZVRESDecayChainCollection << " is not valid." << std::endl;
2352 
2353 
2354  if (!zvresdcCollection.is_valid() || !ftjetCollection.is_valid() || !zvresjetCollection.is_valid() || !zvresdcrptrackCollection.is_valid()) {
2355 
2356  std::cerr << "Vertex-Track Association cannot be analysed. Set MakeAdditionalPlots to false to remove this warning" << std::endl;
2357 
2358  } else {
2359 
2360 
2361  //these are to store the primary, secondary and tertiary vertices in the jet
2362  std::vector<Vertex*> primaryVertices;
2363  std::vector<Vertex*> secondaryVertices;
2364  std::vector<Vertex*> tertiaryVertices;
2365  std::vector<Vertex*> allVertices;
2366  //to store pairs of vertices with a track in between
2367  std::vector< std::pair<Vertex*,Vertex*> > myVertexPairVector;
2368  primaryVertices.clear();
2369  secondaryVertices.clear();
2370  tertiaryVertices.clear();
2371  allVertices.clear();
2372  myVertexPairVector.clear();
2373 
2374 
2375  //loop over all the jets to find the verticies. Only ZVRESDecayChains know about their vertices.
2376 
2377  for (int iDecayChain = 0 ; iDecayChain < zvresdcCollection.getNumberOfElements(); iDecayChain++) {
2378 
2379  lcio::ReconstructedParticle* pDecayChain=zvresdcCollection.getElementAt(iDecayChain);
2380 
2381  //to store pairs of vertices with a track in between
2382  std::pair<Vertex*,Vertex*> myVertexPair;
2383  std::vector< std::pair<Vertex*,Vertex*> > myVertexPairVector;
2384 
2385  myVertexPairVector.clear();
2386 
2387  //these are all the particles within the jet - these guys know about their vertices
2388  std::vector<lcio::ReconstructedParticle*> decayChainPartColl = pDecayChain->getParticles();
2389 
2390  for (unsigned int iPart = 0 ; iPart < decayChainPartColl.size(); iPart++) {
2391 
2392  lcio::ReconstructedParticle* decayChainPart = dynamic_cast< lcio::ReconstructedParticle*> (decayChainPartColl[iPart]);
2393 
2394  myVertexPair.first = decayChainPart->getStartVertex();
2395  myVertexPair.second = decayChainPart->getEndVertex();
2396 
2397  bool pairAlreadyContained = false;
2398 
2399  for (std::vector< std::pair<Vertex*,Vertex*> >::const_iterator iter = myVertexPairVector.begin(); iter != myVertexPairVector.end(); iter++) {
2400  if ((*iter).first == myVertexPair.first && (*iter).second == myVertexPair.second) pairAlreadyContained = true;
2401  }
2402 
2403  if (!pairAlreadyContained) myVertexPairVector.push_back(myVertexPair);
2404  }
2405 
2406 
2407 
2408  //first find the primary, and secondary vertices, and count up the total number of vertices
2409  for (std::vector< std::pair<Vertex*,Vertex*> >::const_iterator iter = myVertexPairVector.begin(); iter != myVertexPairVector.end(); iter++) {
2410 
2411  //look to see if we have a primary vertex. It will always be in the first element, due to the way we set up the pair
2412  if ((*iter).first && (*iter).first->isPrimary())
2413  {
2414  bool primaryAlreadyContained = false;
2415  for (std::vector<Vertex*>::const_iterator piter=primaryVertices.begin(); piter != primaryVertices.end(); piter++)
2416  {
2417  if (*piter == (*iter).first) primaryAlreadyContained = true;
2418  }
2419  if (!primaryAlreadyContained) primaryVertices.push_back((*iter).first);
2420 
2421 
2422  // if the first vertex in the pair is a primary, the second must be a secondary
2423  if ((*iter).second != 0)
2424  {
2425  bool secondaryAlreadyContained = false;
2426  for (std::vector<Vertex*>::const_iterator siter=secondaryVertices.begin(); siter != secondaryVertices.end(); siter++)
2427  {
2428  if (*siter == (*iter).second) secondaryAlreadyContained = true;
2429  }
2430  if (!secondaryAlreadyContained) secondaryVertices.push_back((*iter).second);
2431  }
2432  }
2433 
2434 
2435  // make an exclusive list of vertices
2436  bool firstAlreadyContained=false, secondAlreadyContained=false;
2437  for (std::vector<Vertex*>::const_iterator aiter = allVertices.begin(); aiter != allVertices.end(); aiter++)
2438  {
2439  if ( (*iter).first && *aiter == (*iter).first) firstAlreadyContained=true;
2440  if ( (*iter).second && *aiter == (*iter).second ) secondAlreadyContained=true;
2441 
2442  }
2443  if ( !firstAlreadyContained && (*iter).first) allVertices.push_back((*iter).first);
2444  if (!secondAlreadyContained && (*iter).second) allVertices.push_back((*iter).second);
2445  }
2446 
2447  //next find the tertiary vertices
2448  for (std::vector< std::pair<Vertex*,Vertex*> >::const_iterator iter = myVertexPairVector.begin(); iter != myVertexPairVector.end(); iter++)
2449  {
2450  for (std::vector<Vertex*>::const_iterator siter = secondaryVertices.begin(); siter != secondaryVertices.end(); siter++)
2451  {
2452  if ((*iter).first && (*iter).first == (*siter) && (*iter).second)
2453  {
2454  bool tertiaryAlreadyContained = false;
2455  for (std::vector<Vertex*>::const_iterator titer=tertiaryVertices.begin(); titer != tertiaryVertices.end(); titer++)
2456  {
2457  if (*titer == (*iter).second) tertiaryAlreadyContained=false;
2458  }
2459  if (!tertiaryAlreadyContained) tertiaryVertices.push_back((*iter).second);
2460  }
2461  }
2462 
2463  }
2464 
2465 
2466  }
2467 
2468 
2469  //ok now loop over all the ZVRESSelectedJets
2470  for (int iJet = 0 ; iJet < zvresjetCollection.getNumberOfElements(); iJet++) {
2471  lcio::ReconstructedParticle* pJet=zvresjetCollection.getElementAt(iJet);
2472 
2473  //get the particle within the jets
2474  std::vector<lcio::ReconstructedParticle*> jetParticles = pJet->getParticles();
2475 
2476  //get the true jet flavour
2477  int jetType = FindTrueJetType( pEvent, iJet );
2478 
2479  if (jetType == B_JET || jetType == C_JET) {
2480  //loop over the particles within the jet
2481  for (unsigned int iPart = 0 ; iPart < jetParticles.size(); iPart++)
2482  {
2483  lcio::ReconstructedParticle* myJetParticle = dynamic_cast< lcio::ReconstructedParticle*> (jetParticles[iPart]);
2484 
2485  //do we find a match in the decay chain?
2486  bool decayChainMatch = false;
2487  //do we find a match in the decay chain RPs?
2488  bool decayChainRPMatch = false;
2489 
2490  //look for particle in ZVRESDecayChain collection
2491  lcio::ReconstructedParticle* myDecayChain = 0;
2492  //look for particle in ZVRESDecayChainsRPTrack collection
2493  lcio::ReconstructedParticle* myDecayChainRP = 0;
2494 
2495  //count the number of vertices in the DecayChain
2496  std::vector<Vertex*> myVertexVector;
2497  myVertexVector.clear();
2498 
2499  for (int iDCJet = 0 ; iDCJet < zvresdcrptrackCollection.getNumberOfElements(); iDCJet++) {
2500 
2501  lcio::ReconstructedParticle* myDCJetRP = zvresdcrptrackCollection.getElementAt(iDCJet);
2502 
2503  lcio::ReconstructedParticle* myDCJetRPTrack = myDCJetRP->getParticles()[0];
2504 
2505  if (myDCJetRPTrack==myJetParticle) {
2506  decayChainRPMatch = true;
2507  myDecayChainRP = myDCJetRP;
2508  }
2509  }
2510 
2511  //look for a matching decay chain - if we have found a matching RP - use that information
2512  //otherwise look directly in ZVRESDecayChain and assume they match directly to the ZVRESSelectedJet
2513 
2514  if (decayChainRPMatch) {
2515 
2516  //need to find the decay chain the RP belongs to, in order to find out the number of vertices
2517  for (int iDecayChain = 0 ; iDecayChain < zvresdcCollection.getNumberOfElements(); iDecayChain++) {
2518  lcio::ReconstructedParticle* pDecayChain=zvresdcCollection.getElementAt(iDecayChain);
2519  std::vector<lcio::ReconstructedParticle*> decayChainParticles = pDecayChain->getParticles();
2520  for (unsigned int iDCPart = 0 ; iDCPart < decayChainParticles.size(); iDCPart++) {
2521  lcio::ReconstructedParticle* myDecayChainParticle = dynamic_cast< lcio::ReconstructedParticle*> (decayChainParticles[iDCPart]);
2522  if (myDecayChainParticle == myDecayChainRP) {
2523  decayChainMatch = true;
2524  myDecayChain = pDecayChain;
2525  }
2526  }
2527  }
2528  } else {
2529 
2530  //we are looking at the iJet-th ZVRESSelectedJet
2531  //assume that the iJet-th ZVRESDecayChain refer to the same jet
2532 
2533  lcio::ReconstructedParticle* pDecayChain=zvresdcCollection.getElementAt(iJet);
2534  myDecayChain = pDecayChain;
2535  if (pDecayChain) decayChainMatch = true;
2536  }
2537 
2538  if (decayChainMatch) {
2539 
2540  std::vector<lcio::ReconstructedParticle*> myDecayChainParticles = myDecayChain->getParticles();
2541  for (unsigned int iDCPart = 0 ; iDCPart < myDecayChainParticles.size(); iDCPart++) {
2542 
2543  lcio::ReconstructedParticle* myDecayChainParticle = dynamic_cast< lcio::ReconstructedParticle*> (myDecayChainParticles[iDCPart]);
2544 
2545  bool startAlreadyContained = false;
2546  bool endAlreadyContained = false;
2547 
2548  for (std::vector<Vertex*>::const_iterator iter=myVertexVector.begin(); iter!=myVertexVector.end(); iter++)
2549  {
2550  if (*iter == myDecayChainParticle->getStartVertex()) startAlreadyContained = true;
2551  }
2552  if (!startAlreadyContained && myDecayChainParticle->getStartVertex()) myVertexVector.push_back(myDecayChainParticle->getStartVertex());
2553 
2554 
2555  for (std::vector<Vertex*>::const_iterator iter=myVertexVector.begin(); iter!=myVertexVector.end(); iter++)
2556  {
2557  if (*iter == myDecayChainParticle->getEndVertex()) endAlreadyContained = true;
2558  }
2559  if (!endAlreadyContained && myDecayChainParticle->getEndVertex()) myVertexVector.push_back(myDecayChainParticle->getEndVertex());
2560 
2561  }
2562  }
2563 
2564  int numVertex=myVertexVector.size();
2565 
2566  bool trackFromPrimaryLight = false, trackFromSecondaryC = false, trackFromSecondaryB = false, trackFromTertiaryC = false;
2567  bool trackHasNoAssociatedMCP = false;
2568 
2569  //find the origin of these tracks through the MCRelationNavigator.
2570  std::vector<Track*> myTracks;
2571 
2572  if (decayChainRPMatch) myTracks = myDecayChainRP->getParticles()[0]->getTracks();
2573  else {
2574  myTracks = myJetParticle->getTracks();
2575 
2576  }
2577 
2578  if (myTracks.size()==0)
2579  {
2580 
2581  //look at the matching MC version info
2582 
2583  trackHasNoAssociatedMCP = true;
2584 
2585  }
2586 
2587  else if (myTracks.size()>0)
2588  {
2589 
2590  lcio::Track* Track = myTracks[0];
2591 
2592  std::vector<lcio::LCObject*> RelatedMCParticles = MCRelationNavigator->getRelatedToObjects(Track);
2593 
2594  //VJM 14/05/08
2595  //If more than one related MC particle, seems (on a brief visual inspection) that they just all comes from the same original hadron
2596  // - just from different parts of the decay chain.
2597  //So at the moment I'll just use the first MCP in the list. It will (should) give the same answer as any other.
2598 
2599  lcio::MCParticle* MCP = dynamic_cast<lcio::MCParticle*>(RelatedMCParticles[0]);
2600  std::vector<lcio::MCParticle*> Parents = dynamic_cast<lcio::MCParticle*>(RelatedMCParticles[0])->getParents();
2601 
2602  if (Parents.empty())
2603  std::cerr << " In " << __FILE__ << "(" << __LINE__ << "), run number: " << pEvent->getRunNumber() << ", event number: " << pEvent->getEventNumber()<< " : MCP has no parents" << std::endl;
2604  //ERROR NO PARENT;
2605  else
2606  {
2607  int truePDGCode=MCP->getPDG();
2608  int truePDGFlavour = GetPDGFlavour(truePDGCode);
2609  std::vector<int> pdgCodeVector;
2610  std::vector<int> pdgFlavourVector;
2611  std::vector<float> pdgDecayLengths;
2612 
2613  pdgCodeVector.push_back(truePDGCode);
2614  pdgFlavourVector.push_back(truePDGFlavour);
2615 
2616  std::vector<lcio::MCParticle*> ParentVec0 = MCP->getParents();
2617 
2618  //If there are two parents, then they point to the same original flavour
2619  //if (ParentVec0.size()>1) {
2620  // std::cerr << " In " << __FILE__ << "(" << __LINE__ << "), run number: " << pEvent->getRunNumber() << ", event number: ";
2621  // std::cerr << pEvent->getEventNumber()<< " ; Size of ParentVec0: " << ParentVec0.size() << std::endl;
2622  //}
2623  if (ParentVec0.size()>0) {
2624  lcio::MCParticle* Parent0 = ParentVec0[0];
2625 
2626  if (Parent0) {
2627 
2628  pdgCodeVector.push_back(Parent0->getPDG());
2629  pdgFlavourVector.push_back(GetPDGFlavour(Parent0->getPDG()));
2630  std::vector<lcio::MCParticle*> ParentVec1 = Parent0->getParents();
2631  pdgDecayLengths.push_back(CalculateDistance(Parent0->getVertex(),Parent0->getEndpoint()));
2632 
2633  if (ParentVec1.size()>0) {
2634  lcio::MCParticle* Parent1 = ParentVec1[0];
2635 
2636  if (Parent1) {
2637  pdgCodeVector.push_back(Parent1->getPDG());
2638  pdgFlavourVector.push_back(GetPDGFlavour(Parent1->getPDG()));
2639  std::vector<lcio::MCParticle*> ParentVec2 = Parent1->getParents();
2640  pdgDecayLengths.push_back(CalculateDistance(Parent1->getVertex(),Parent1->getEndpoint()));
2641 
2642  if (ParentVec2.size()>0) {
2643  lcio::MCParticle* Parent2 = ParentVec2[0];
2644 
2645  if (Parent2) {
2646  pdgCodeVector.push_back(Parent2->getPDG());
2647  pdgFlavourVector.push_back(GetPDGFlavour(Parent2->getPDG()));
2648  std::vector<lcio::MCParticle*> ParentVec3 = Parent2->getParents();
2649  pdgDecayLengths.push_back(CalculateDistance(Parent2->getVertex(),Parent2->getEndpoint()));
2650 
2651  if (ParentVec3.size()>0) {
2652  lcio::MCParticle* Parent3 = ParentVec3[0];
2653 
2654  if (Parent3) {
2655  pdgCodeVector.push_back(Parent3->getPDG());
2656  pdgFlavourVector.push_back(GetPDGFlavour(Parent3->getPDG()));
2657  std::vector<lcio::MCParticle*> ParentVec4 = Parent3->getParents();
2658  pdgDecayLengths.push_back(CalculateDistance(Parent3->getVertex(),Parent3->getEndpoint()));
2659 
2660  if (ParentVec4.size()>0) {
2661  lcio::MCParticle* Parent4 = ParentVec4[0];
2662 
2663  if (Parent4) {
2664  pdgCodeVector.push_back(Parent4->getPDG());
2665  pdgFlavourVector.push_back(GetPDGFlavour(Parent4->getPDG()));
2666  std::vector<lcio::MCParticle*> ParentVec5 = Parent4->getParents();
2667  pdgDecayLengths.push_back(CalculateDistance(Parent4->getVertex(),Parent4->getEndpoint()));
2668 
2669  if (ParentVec5.size()>0) {
2670  lcio::MCParticle* Parent5 = ParentVec5[0];
2671 
2672  if (Parent5) {
2673  pdgCodeVector.push_back(Parent5->getPDG());
2674  pdgFlavourVector.push_back(GetPDGFlavour(Parent5->getPDG()));
2675  std::vector<lcio::MCParticle*> ParentVec6 = Parent5->getParents();
2676  pdgDecayLengths.push_back(CalculateDistance(Parent5->getVertex(),Parent5->getEndpoint()));
2677 
2678  if (ParentVec6.size()>0) {
2679  lcio::MCParticle* Parent6 = ParentVec6[0];
2680 
2681  if (Parent6) {
2682  pdgCodeVector.push_back(Parent6->getPDG());
2683  pdgFlavourVector.push_back(GetPDGFlavour(Parent6->getPDG()));
2684  pdgDecayLengths.push_back(CalculateDistance(Parent6->getVertex(),Parent6->getEndpoint()));
2685 
2686  } //eh hope that's enough...
2687  }
2688  }
2689  }
2690  }
2691  }
2692  }
2693  }
2694  }
2695  }
2696  }
2697  }
2698  }
2699  }
2700 
2701 
2702  //assign a flavour to the jet.
2703  bool isLight = false;
2704  bool isC = false;
2705  bool isB = false;
2706 
2707 
2708  for (std::vector<int>::const_iterator iter = pdgFlavourVector.begin();
2709  iter != pdgFlavourVector.end(); iter++) {
2710  if ((*iter) == 1) isLight = true;
2711  if ((*iter) == C_JET) isC = true;
2712  if ((*iter) == B_JET) isB = true;
2713  }
2714 
2715 
2716  //if it's all "1"s --> then primary
2717  //if there's a "4" and no "5" --> then D (secondary)
2718  //if there's a "4" AND a "5" --> then B and D (secondary and tertiary)
2719 
2720  if (isLight && !isC && !isB) trackFromPrimaryLight = true;
2721  if (isC && !isB) trackFromSecondaryC = true;
2722  if (!isC && isB) trackFromSecondaryB = true;
2723  if (isC && isB) trackFromTertiaryC = true;
2724 
2725  }
2726  }
2727 
2728  //find the vertex depth - only possible if we've found the RP corresponding to the track
2729 
2730  bool isPrimary=false, isSecondary=false, isTertiary=false;
2731 
2732  if (decayChainRPMatch) {
2733 
2734  lcio::Vertex* partVertex = myDecayChainRP->getStartVertex();
2735 
2736  for (std::vector<Vertex*>::const_iterator piter = primaryVertices.begin(); piter != primaryVertices.end(); piter++)
2737  {
2738  if ((*piter) == partVertex) isPrimary=true;
2739  }
2740  for (std::vector<Vertex*>::const_iterator siter = secondaryVertices.begin(); siter != secondaryVertices.end(); siter++)
2741  {
2742  if ((*siter) == partVertex) isSecondary=true;
2743  }
2744  for (std::vector<Vertex*>::const_iterator titer = tertiaryVertices.begin(); titer != tertiaryVertices.end(); titer++)
2745  {
2746  if ((*titer) == partVertex) isTertiary=true;
2747  }
2748 
2749  }
2750 
2751  if (numVertex == 2)
2752  {
2753  if (trackHasNoAssociatedMCP)
2754  {
2755  if (isPrimary) _nb_twoVertex_Primary_noMCP++;
2756  else if (isSecondary) _nb_twoVertex_Secondary_noMCP++;
2757  else _nb_twoVertex_Isolated_noMCP++;
2758  }
2759 
2760  if (trackFromSecondaryB && jetType == B_JET )
2761  {
2762  if (isPrimary) _nb_twoVertex_bTrack_Primary++;
2763  else if (isSecondary) _nb_twoVertex_bTrack_Secondary++;
2764  //else if (isTertiary) _nb_twoVertex_bTrack_Tertiary++;
2765  else _nb_twoVertex_bTrack_Isolated++;
2766  }
2767  else if ((trackFromSecondaryC || trackFromTertiaryC) && jetType == B_JET )
2768  {
2769  if (isPrimary) _nb_twoVertex_cTrack_Primary++;
2770  else if (isSecondary) _nb_twoVertex_cTrack_Secondary++;
2771  //else if (isTertiary) _nb_twoVertex_cTrack_Tertiary++;
2772  else _nb_twoVertex_cTrack_Isolated++;
2773  }
2774  else if (trackFromPrimaryLight && jetType == B_JET)
2775  {
2776  if (isPrimary) _nb_twoVertex_lTrack_Primary++;
2777  else if (isSecondary) _nb_twoVertex_lTrack_Secondary++;
2778  //else if (isTertiary) _nb_twoVertex_lTrack_Tertiary++;
2779  else _nb_twoVertex_lTrack_Isolated++;
2780  }
2781 
2782  if (trackHasNoAssociatedMCP && jetType == C_JET) {
2783 
2784  if (isPrimary) _nc_twoVertex_Primary_noMCP++;
2785  else if (isSecondary) _nc_twoVertex_Secondary_noMCP++;
2786  else _nc_twoVertex_Isolated_noMCP++;
2787  }
2788 
2789  if (trackFromSecondaryB && jetType == C_JET )
2790  {
2791  if (isPrimary) _nc_twoVertex_bTrack_Primary++;
2792  else if (isSecondary) _nc_twoVertex_bTrack_Secondary++;
2793  //else if (isTertiary) _nc_twoVertex_bTrack_Tertiary++;
2794  else _nc_twoVertex_bTrack_Isolated++;
2795  }
2796  else if ((trackFromSecondaryC || trackFromTertiaryC) && jetType == C_JET)
2797  {
2798  if (isPrimary) _nc_twoVertex_cTrack_Primary++;
2799  else if (isSecondary) _nc_twoVertex_cTrack_Secondary++;
2800  //else if (isTertiary) _nc_twoVertex_cTrack_Tertiary++;
2801  else _nc_twoVertex_cTrack_Isolated++;
2802  }
2803  else if (trackFromPrimaryLight && jetType == C_JET)
2804  {
2805  if (isPrimary) _nc_twoVertex_lTrack_Primary++;
2806  else if (isSecondary) _nc_twoVertex_lTrack_Secondary++;
2807  //else if (isTertiary) _nc_twoVertex_lTrack_Tertiary++;
2808  else _nc_twoVertex_lTrack_Isolated++;
2809  }
2810  }
2811  else if (numVertex == 3 )
2812  {
2813  if (trackHasNoAssociatedMCP && jetType == B_JET)
2814  {
2815  if (isPrimary) _nb_twoVertex_Primary_noMCP++;
2816  else if (isSecondary) _nb_twoVertex_Secondary_noMCP++;
2817  else if (isTertiary) _nb_twoVertex_Tertiary_noMCP++;
2818  else _nb_twoVertex_Isolated_noMCP++;
2819  }
2820 
2821  if (trackFromSecondaryB && jetType == B_JET)
2822  {
2823  if (isPrimary) _nb_threeVertex_bTrack_Primary++;
2824  else if (isSecondary) _nb_threeVertex_bTrack_Secondary++;
2825  else if (isTertiary) _nb_threeVertex_bTrack_Tertiary++;
2826  else _nb_threeVertex_bTrack_Isolated++;
2827  }
2828  else if ((trackFromSecondaryC || trackFromTertiaryC) && jetType == B_JET)
2829  {
2830  if (isPrimary) _nb_threeVertex_cTrack_Primary++;
2831  else if (isSecondary) _nb_threeVertex_cTrack_Secondary++;
2832  else if (isTertiary) _nb_threeVertex_cTrack_Tertiary++;
2833  else _nb_threeVertex_cTrack_Isolated++;
2834  }
2835  else if (trackFromPrimaryLight && jetType == B_JET)
2836  {
2837  if (isPrimary) _nb_threeVertex_lTrack_Primary++;
2838  else if (isSecondary) _nb_threeVertex_lTrack_Secondary++;
2839  else if (isTertiary) _nb_threeVertex_lTrack_Tertiary++;
2840  else _nb_threeVertex_lTrack_Isolated++;
2841  }
2842  if (trackHasNoAssociatedMCP && jetType == C_JET)
2843  {
2844  if (isPrimary) _nc_twoVertex_Primary_noMCP++;
2845  else if (isSecondary) _nc_twoVertex_Secondary_noMCP++;
2846  else if (isTertiary) _nc_twoVertex_Tertiary_noMCP++;
2847  else _nc_twoVertex_Isolated_noMCP++;
2848  }
2849  else if (trackFromSecondaryB && jetType == C_JET)
2850  {
2851  if (isPrimary) _nc_threeVertex_bTrack_Primary++;
2852  else if (isSecondary) _nc_threeVertex_bTrack_Secondary++;
2853  else if (isTertiary) _nc_threeVertex_bTrack_Tertiary++;
2854  else _nc_threeVertex_bTrack_Isolated++;
2855  }
2856  else if ((trackFromSecondaryC || trackFromTertiaryC) && jetType == C_JET)
2857  {
2858  if (isPrimary) _nc_threeVertex_cTrack_Primary++;
2859  else if (isSecondary) _nc_threeVertex_cTrack_Secondary++;
2860  else if (isTertiary) _nc_threeVertex_cTrack_Tertiary++;
2861  else _nc_threeVertex_cTrack_Isolated++;
2862  }
2863  else if (trackFromPrimaryLight && jetType == C_JET)
2864  {
2865  if (isPrimary) _nc_threeVertex_lTrack_Primary++;
2866  else if (isSecondary) _nc_threeVertex_lTrack_Secondary++;
2867  else if (isTertiary) _nc_threeVertex_lTrack_Tertiary++;
2868  else _nc_threeVertex_lTrack_Isolated++;
2869  }
2870  }
2871  }
2872  } //jet type is B_JET or C_JET
2873  }
2874  }//collections are valid
2875 }
2876 
2877 float LCFIAIDAPlotProcessor::CalculateDistance(const float* pos1, const float* pos2){
2878  return sqrt(pow((pos1[0]-pos2[0]),2)+pow((pos1[1]-pos2[1]),2)+pow((pos1[2]-pos2[2]),2));
2879 }
2880 
2881 double LCFIAIDAPlotProcessor::CalculateDistance(const double* pos1, const double* pos2){
2882  return sqrt(pow((pos1[0]-pos2[0]),2)+pow((pos1[1]-pos2[1]),2)+pow((pos1[2]-pos2[2]),2));
2883 }
2884 
2885 void LCFIAIDAPlotProcessor::PrintZVRESTable()
2886 {
2887 
2888  //if there is a _TrackVertexOutputFile string defined use that as the output stream, if not use std::cout
2889  std::filebuf* fb = new std::filebuf;
2890 
2891  std::ostream outputFile( (!_TrackVertexOutputFile.empty()) ?
2892  fb->open(_TrackVertexOutputFile.c_str(),
2893  std::ios_base::out|std::ios_base::trunc):
2894  std::cout.rdbuf());
2895 
2896 
2897  //no longer required
2898  //float twoVertex_bTrack = _nb_twoVertex_bTrack_Primary+_nb_twoVertex_bTrack_Secondary+_nb_twoVertex_bTrack_Isolated;
2899  //float twoVertex_cTrack = _nb_twoVertex_cTrack_Primary+_nb_twoVertex_cTrack_Secondary+_nb_twoVertex_cTrack_Isolated;
2900  //float twoVertex_lTrack = _nb_twoVertex_lTrack_Primary+_nb_twoVertex_lTrack_Secondary+_nb_twoVertex_lTrack_Isolated;
2901 
2902  //float threeVertex_bTrack = _nb_threeVertex_bTrack_Primary+_nb_threeVertex_bTrack_Secondary+_nb_threeVertex_bTrack_Tertiary+_nb_threeVertex_bTrack_Isolated;
2903  //float threeVertex_cTrack = _nb_threeVertex_cTrack_Primary+_nb_threeVertex_cTrack_Secondary+_nb_threeVertex_cTrack_Tertiary+_nb_threeVertex_cTrack_Isolated;
2904  //float threeVertex_lTrack = _nb_threeVertex_lTrack_Primary+_nb_threeVertex_lTrack_Secondary+_nb_threeVertex_lTrack_Tertiary+_nb_threeVertex_lTrack_Isolated;
2905 
2906 
2907  float b_twoVertex_Primary = _nb_twoVertex_bTrack_Primary+_nb_twoVertex_cTrack_Primary+_nb_twoVertex_lTrack_Primary;
2908  float b_twoVertex_Secondary = _nb_twoVertex_bTrack_Secondary+_nb_twoVertex_cTrack_Secondary+_nb_twoVertex_lTrack_Secondary;
2909  float b_twoVertex_Isolated = _nb_twoVertex_bTrack_Isolated+_nb_twoVertex_cTrack_Isolated+_nb_twoVertex_lTrack_Isolated;
2910  float b_twoVertex = b_twoVertex_Primary+b_twoVertex_Secondary+b_twoVertex_Isolated;
2911 
2912  float b_threeVertex_Primary = _nb_threeVertex_bTrack_Primary+_nb_threeVertex_cTrack_Primary+_nb_threeVertex_lTrack_Primary;
2913  float b_threeVertex_Secondary = _nb_threeVertex_bTrack_Secondary+_nb_threeVertex_cTrack_Secondary+_nb_threeVertex_lTrack_Secondary;
2914  float b_threeVertex_Tertiary = _nb_threeVertex_bTrack_Tertiary+_nb_threeVertex_cTrack_Tertiary+_nb_threeVertex_lTrack_Tertiary;
2915  float b_threeVertex_Isolated = _nb_threeVertex_bTrack_Isolated+_nb_threeVertex_cTrack_Isolated+_nb_threeVertex_lTrack_Isolated;
2916  float b_threeVertex = b_threeVertex_Primary+b_threeVertex_Secondary+b_threeVertex_Tertiary+b_threeVertex_Isolated;
2917 
2918 
2919 
2920  float frac_b_twoVertex_bTrack_Primary = 100.*float(_nb_twoVertex_bTrack_Primary) / float(b_twoVertex_Primary);
2921  float frac_b_twoVertex_bTrack_Secondary = 100.*float(_nb_twoVertex_bTrack_Secondary) / float(b_twoVertex_Secondary);
2922  //float frac_b_twoVertex_bTrack_Tertiary = 100.*float(_nb_twoVertex_bTrack_Tertiary) / float(b_twoVertex_Tertiary);
2923  float frac_b_twoVertex_bTrack_Isolated = 100.*float(_nb_twoVertex_bTrack_Isolated) / float(b_twoVertex_Isolated);
2924 
2925  float frac_b_twoVertex_cTrack_Primary = 100.*float(_nb_twoVertex_cTrack_Primary) / float(b_twoVertex_Primary);
2926  float frac_b_twoVertex_cTrack_Secondary = 100.*float(_nb_twoVertex_cTrack_Secondary) / float(b_twoVertex_Secondary);
2927  //float frac_b_twoVertex_cTrack_Tertiary = 100.*float(_nb_twoVertex_cTrack_Tertiary) / float(b_twoVertex_Tertiary);
2928  float frac_b_twoVertex_cTrack_Isolated = 100.*float(_nb_twoVertex_cTrack_Isolated) / float(b_twoVertex_Isolated);
2929 
2930  float frac_b_twoVertex_lTrack_Primary = 100.*float(_nb_twoVertex_lTrack_Primary) / float(b_twoVertex_Primary);
2931  float frac_b_twoVertex_lTrack_Secondary = 100.*float(_nb_twoVertex_lTrack_Secondary) / float(b_twoVertex_Secondary);
2932  //float frac_b_twoVertex_lTrack_Tertiary = 100.*float(_nb_twoVertex_lTrack_Tertiary) / float(b_twoVertex_Tertiary);
2933  float frac_b_twoVertex_lTrack_Isolated = 100.*float(_nb_twoVertex_lTrack_Isolated) / float(b_twoVertex_Isolated);
2934 
2935  float frac_b_threeVertex_bTrack_Primary = 100.*float(_nb_threeVertex_bTrack_Primary) / float(b_threeVertex_Primary);
2936  float frac_b_threeVertex_bTrack_Secondary = 100.*float(_nb_threeVertex_bTrack_Secondary) / float(b_threeVertex_Secondary);
2937  float frac_b_threeVertex_bTrack_Tertiary = 100.*float(_nb_threeVertex_bTrack_Tertiary) / float(b_threeVertex_Tertiary);
2938  float frac_b_threeVertex_bTrack_Isolated = 100.*float(_nb_threeVertex_bTrack_Isolated) / float(b_threeVertex_Isolated);
2939 
2940  float frac_b_threeVertex_cTrack_Primary = 100.*float(_nb_threeVertex_cTrack_Primary) / float(b_threeVertex_Primary);
2941  float frac_b_threeVertex_cTrack_Secondary = 100.*float(_nb_threeVertex_cTrack_Secondary) / float(b_threeVertex_Secondary);
2942  float frac_b_threeVertex_cTrack_Tertiary = 100.*float(_nb_threeVertex_cTrack_Tertiary) / float(b_threeVertex_Tertiary);
2943  float frac_b_threeVertex_cTrack_Isolated = 100.*float(_nb_threeVertex_cTrack_Isolated) / float(b_threeVertex_Isolated);
2944 
2945  float frac_b_threeVertex_lTrack_Primary = 100.*float(_nb_threeVertex_lTrack_Primary) / float(b_threeVertex_Primary);
2946  float frac_b_threeVertex_lTrack_Secondary = 100.*float(_nb_threeVertex_lTrack_Secondary) / float(b_threeVertex_Secondary);
2947  float frac_b_threeVertex_lTrack_Tertiary = 100.*float(_nb_threeVertex_lTrack_Tertiary) / float(b_threeVertex_Tertiary);
2948  float frac_b_threeVertex_lTrack_Isolated = 100.*float(_nb_threeVertex_lTrack_Isolated) / float(b_threeVertex_Isolated);
2949 
2950  float frac_b_twoVertex_Primary = 100.*float(b_twoVertex_Primary) / float(b_twoVertex);
2951  float frac_b_twoVertex_Secondary = 100.*float(b_twoVertex_Secondary) / float(b_twoVertex);
2952  float frac_b_twoVertex_Isolated = 100.*float(b_twoVertex_Isolated) / float(b_twoVertex);
2953 
2954  float frac_b_threeVertex_Primary = 100.*float(b_threeVertex_Primary) / float(b_threeVertex);
2955  float frac_b_threeVertex_Secondary = 100.*float(b_threeVertex_Secondary) / float(b_threeVertex);
2956  float frac_b_threeVertex_Tertiary = 100.*float(b_threeVertex_Tertiary) / float(b_threeVertex);
2957  float frac_b_threeVertex_Isolated = 100.*float(b_threeVertex_Isolated) / float(b_threeVertex);
2958 
2959 
2960 
2961 //check// //just for checks...
2962 //check// frac_b_twoVertex_bTrack_Primary = float(_nb_twoVertex_bTrack_Primary);
2963 //check// frac_b_twoVertex_bTrack_Secondary = float(_nb_twoVertex_bTrack_Secondary);
2964 //check// //frac_b_twoVertex_bTrack_Tertiary = float(_nb_twoVertex_bTrack_Tertiary);
2965 //check// frac_b_twoVertex_bTrack_Isolated = float(_nb_twoVertex_bTrack_Isolated);
2966 //check//
2967 //check// frac_b_twoVertex_cTrack_Primary = float(_nb_twoVertex_cTrack_Primary);
2968 //check// frac_b_twoVertex_cTrack_Secondary = float(_nb_twoVertex_cTrack_Secondary);
2969 //check// //frac_b_twoVertex_cTrack_Tertiary = float(_nb_twoVertex_cTrack_Tertiary);
2970 //check// frac_b_twoVertex_cTrack_Isolated = float(_nb_twoVertex_cTrack_Isolated);
2971 //check//
2972 //check// frac_b_twoVertex_lTrack_Primary = float(_nb_twoVertex_lTrack_Primary);
2973 //check// frac_b_twoVertex_lTrack_Secondary = float(_nb_twoVertex_lTrack_Secondary);
2974 //check// //frac_b_twoVertex_lTrack_Tertiary = float(_nb_twoVertex_lTrack_Tertiary);
2975 //check// frac_b_twoVertex_lTrack_Isolated = float(_nb_twoVertex_lTrack_Isolated);
2976 //check//
2977 //check// frac_b_threeVertex_bTrack_Primary = float(_nb_threeVertex_bTrack_Primary);
2978 //check// frac_b_threeVertex_bTrack_Secondary = float(_nb_threeVertex_bTrack_Secondary);
2979 //check// frac_b_threeVertex_bTrack_Tertiary = float(_nb_threeVertex_bTrack_Tertiary);
2980 //check// frac_b_threeVertex_bTrack_Isolated = float(_nb_threeVertex_bTrack_Isolated);
2981 //check//
2982 //check// frac_b_threeVertex_cTrack_Primary = float(_nb_threeVertex_cTrack_Primary);
2983 //check// frac_b_threeVertex_cTrack_Secondary = float(_nb_threeVertex_cTrack_Secondary);
2984 //check// frac_b_threeVertex_cTrack_Tertiary = float(_nb_threeVertex_cTrack_Tertiary);
2985 //check// frac_b_threeVertex_cTrack_Isolated = float(_nb_threeVertex_cTrack_Isolated);
2986 //check//
2987 //check// frac_b_threeVertex_lTrack_Primary = float(_nb_threeVertex_lTrack_Primary);
2988 //check// frac_b_threeVertex_lTrack_Secondary = float(_nb_threeVertex_lTrack_Secondary);
2989 //check// frac_b_threeVertex_lTrack_Tertiary = float(_nb_threeVertex_lTrack_Tertiary);
2990 //check// frac_b_threeVertex_lTrack_Isolated = float(_nb_threeVertex_lTrack_Isolated);
2991 //check// //end of checks
2992 
2993 
2994 
2995  float c_twoVertex_Primary = _nc_twoVertex_bTrack_Primary+_nc_twoVertex_cTrack_Primary+_nc_twoVertex_lTrack_Primary;
2996  float c_twoVertex_Secondary = _nc_twoVertex_bTrack_Secondary+_nc_twoVertex_cTrack_Secondary+_nc_twoVertex_lTrack_Secondary;
2997  float c_twoVertex_Isolated = _nc_twoVertex_bTrack_Isolated+_nc_twoVertex_cTrack_Isolated+_nc_twoVertex_lTrack_Isolated;
2998  float c_twoVertex = c_twoVertex_Primary+c_twoVertex_Secondary+c_twoVertex_Isolated;
2999 
3000  float c_threeVertex_Primary = _nc_threeVertex_bTrack_Primary+_nc_threeVertex_cTrack_Primary+_nc_threeVertex_lTrack_Primary;
3001  float c_threeVertex_Secondary = _nc_threeVertex_bTrack_Secondary+_nc_threeVertex_cTrack_Secondary+_nc_threeVertex_lTrack_Secondary;
3002  float c_threeVertex_Tertiary = _nc_threeVertex_bTrack_Tertiary+_nc_threeVertex_cTrack_Tertiary+_nc_threeVertex_lTrack_Tertiary;
3003  float c_threeVertex_Isolated = _nc_threeVertex_bTrack_Isolated+_nc_threeVertex_cTrack_Isolated+_nc_threeVertex_lTrack_Isolated;
3004  float c_threeVertex = c_threeVertex_Primary+c_threeVertex_Secondary+c_threeVertex_Tertiary+c_threeVertex_Isolated;
3005 
3006 
3007 
3008  // float frac_c_twoVertex_bTrack_Primary = 100.*float(_nc_twoVertex_bTrack_Primary) / float(c_twoVertex_Primary);
3009  //float frac_c_twoVertex_bTrack_Secondary = 100.*float(_nc_twoVertex_bTrack_Secondary) / float(c_twoVertex_Secondary);
3010  //float frac_c_twoVertex_bTrack_Tertiary = 100.*float(_nc_twoVertex_bTrack_Tertiary) / float(c_twoVertex_Tertiary);
3011  //float frac_c_twoVertex_bTrack_Isolated = 100.*float(_nc_twoVertex_bTrack_Isolated) / float(c_twoVertex_Isolated);
3012 
3013  float frac_c_twoVertex_cTrack_Primary = 100.*float(_nc_twoVertex_cTrack_Primary) / float(c_twoVertex_Primary);
3014  float frac_c_twoVertex_cTrack_Secondary = 100.*float(_nc_twoVertex_cTrack_Secondary) / float(c_twoVertex_Secondary);
3015  //float frac_c_twoVertex_cTrack_Tertiary = 100.*float(_nc_twoVertex_cTrack_Tertiary) / float(c_twoVertex_Tertiary);
3016  float frac_c_twoVertex_cTrack_Isolated = 100.*float(_nc_twoVertex_cTrack_Isolated) / float(c_twoVertex_Isolated);
3017 
3018  float frac_c_twoVertex_lTrack_Primary = 100.*float(_nc_twoVertex_lTrack_Primary) / float(c_twoVertex_Primary);
3019  float frac_c_twoVertex_lTrack_Secondary = 100.*float(_nc_twoVertex_lTrack_Secondary) / float(c_twoVertex_Secondary);
3020  //float frac_c_twoVertex_lTrack_Tertiary = 100.*float(_nc_twoVertex_lTrack_Tertiary) / float(c_twoVertex_Tertiary);
3021  float frac_c_twoVertex_lTrack_Isolated = 100.*float(_nc_twoVertex_lTrack_Isolated) / float(c_twoVertex_Isolated);
3022 
3023  //float frac_c_threeVertex_bTrack_Primary = 100.*float(_nc_threeVertex_bTrack_Primary) / float(c_threeVertex_Primary);
3024  //float frac_c_threeVertex_bTrack_Secondary = 100.*float(_nc_threeVertex_bTrack_Secondary) / float(c_threeVertex_Secondary);
3025  //float frac_c_threeVertex_bTrack_Tertiary = 100.*float(_nc_threeVertex_bTrack_Tertiary) / float(c_threeVertex_Tertiary);
3026  //float frac_c_threeVertex_bTrack_Isolated = 100.*float(_nc_threeVertex_bTrack_Isolated) / float(c_threeVertex_Isolated);
3027 
3028  float frac_c_threeVertex_cTrack_Primary = 100.*float(_nc_threeVertex_cTrack_Primary) / float(c_threeVertex_Primary);
3029  float frac_c_threeVertex_cTrack_Secondary = 100.*float(_nc_threeVertex_cTrack_Secondary) / float(c_threeVertex_Secondary);
3030  float frac_c_threeVertex_cTrack_Tertiary = 100.*float(_nc_threeVertex_cTrack_Tertiary) / float(c_threeVertex_Tertiary);
3031  float frac_c_threeVertex_cTrack_Isolated = 100.*float(_nc_threeVertex_cTrack_Isolated) / float(c_threeVertex_Isolated);
3032 
3033  float frac_c_threeVertex_lTrack_Primary = 100.*float(_nc_threeVertex_lTrack_Primary) / float(c_threeVertex_Primary);
3034  float frac_c_threeVertex_lTrack_Secondary = 100.*float(_nc_threeVertex_lTrack_Secondary) / float(c_threeVertex_Secondary);
3035  float frac_c_threeVertex_lTrack_Tertiary = 100.*float(_nc_threeVertex_lTrack_Tertiary) / float(c_threeVertex_Tertiary);
3036  float frac_c_threeVertex_lTrack_Isolated = 100.*float(_nc_threeVertex_lTrack_Isolated) / float(c_threeVertex_Isolated);
3037 
3038  float frac_c_twoVertex_Primary = 100.*float(c_twoVertex_Primary) / float(c_twoVertex);
3039  float frac_c_twoVertex_Secondary = 100.*float(c_twoVertex_Secondary) / float(c_twoVertex);
3040  float frac_c_twoVertex_Isolated = 100.*float(c_twoVertex_Isolated) / float(c_twoVertex);
3041 
3042  float frac_c_threeVertex_Primary = 100.*float(c_threeVertex_Primary) / float(c_threeVertex);
3043  float frac_c_threeVertex_Secondary = 100.*float(c_threeVertex_Secondary) / float(c_threeVertex);
3044  float frac_c_threeVertex_Tertiary = 100.*float(c_threeVertex_Tertiary) / float(c_threeVertex);
3045  float frac_c_threeVertex_Isolated = 100.*float(c_threeVertex_Isolated) / float(c_threeVertex);
3046 
3047 
3048 
3049 //check// //just for checks...
3050 //check// frac_c_twoVertex_bTrack_Primary = float(_nc_twoVertex_bTrack_Primary);
3051 //check// frac_c_twoVertex_bTrack_Secondary = float(_nc_twoVertex_bTrack_Secondary);
3052 //check// //frac_c_twoVertex_bTrack_Tertiary = float(_nc_twoVertex_bTrack_Tertiary);
3053 //check// frac_c_twoVertex_bTrack_Isolated = float(_nc_twoVertex_bTrack_Isolated);
3054 //check//
3055 //check// frac_c_twoVertex_cTrack_Primary = float(_nc_twoVertex_cTrack_Primary);
3056 //check// frac_c_twoVertex_cTrack_Secondary = float(_nc_twoVertex_cTrack_Secondary);
3057 //check// //frac_c_twoVertex_cTrack_Tertiary = float(_nc_twoVertex_cTrack_Tertiary);
3058 //check// frac_c_twoVertex_cTrack_Isolated = float(_nc_twoVertex_cTrack_Isolated);
3059 //check//
3060 //check// frac_c_twoVertex_lTrack_Primary = float(_nc_twoVertex_lTrack_Primary);
3061 //check// frac_c_twoVertex_lTrack_Secondary = float(_nc_twoVertex_lTrack_Secondary);
3062 //check// //frac_c_twoVertex_lTrack_Tertiary = float(_nc_twoVertex_lTrack_Tertiary);
3063 //check// frac_c_twoVertex_lTrack_Isolated = float(_nc_twoVertex_lTrack_Isolated);
3064 //check//
3065 //check// frac_c_threeVertex_bTrack_Primary = float(_nc_threeVertex_bTrack_Primary);
3066 //check// frac_c_threeVertex_bTrack_Secondary = float(_nc_threeVertex_bTrack_Secondary);
3067 //check// frac_c_threeVertex_bTrack_Tertiary = float(_nc_threeVertex_bTrack_Tertiary);
3068 //check// frac_c_threeVertex_bTrack_Isolated = float(_nc_threeVertex_bTrack_Isolated);
3069 //check//
3070 //check// frac_c_threeVertex_cTrack_Primary = float(_nc_threeVertex_cTrack_Primary);
3071 //check// frac_c_threeVertex_cTrack_Secondary = float(_nc_threeVertex_cTrack_Secondary);
3072 //check// frac_c_threeVertex_cTrack_Tertiary = float(_nc_threeVertex_cTrack_Tertiary);
3073 //check// frac_c_threeVertex_cTrack_Isolated = float(_nc_threeVertex_cTrack_Isolated);
3074 //check//
3075 //check// frac_c_threeVertex_lTrack_Primary = float(_nc_threeVertex_lTrack_Primary);
3076 //check// frac_c_threeVertex_lTrack_Secondary = float(_nc_threeVertex_lTrack_Secondary);
3077 //check// frac_c_threeVertex_lTrack_Tertiary = float(_nc_threeVertex_lTrack_Tertiary);
3078 //check// frac_c_threeVertex_lTrack_Isolated = float(_nc_threeVertex_lTrack_Isolated);
3079 //check// //end of checks
3080 
3081 
3082 
3083  outputFile << std::endl;
3084  outputFile << " ---------------------------------------------------------------------------" << std::endl;
3085  outputFile << " B-JETS " << std::endl;
3086  outputFile << " ---------------------------------------------------------------------------" << std::endl;
3087  outputFile << " Purity of reconstructed track-vertex association (%) " << std::endl;
3088  outputFile << " ---------------------------------------------------------------------------" << std::endl;
3089  outputFile << " Monte Carlo Reconstructed track-vertex association " << std::endl;
3090  outputFile << " track origin Two-vertex case Three-vertex case " << std::endl;
3091  outputFile << " ---------------------------------------------------------------------------" << std::endl;
3092  outputFile << " Pri. Sec. Iso. Pri. Sec. Ter. Iso. " << std::endl;
3093  outputFile << " ---------------------------------------------------------------------------" << std::endl;
3094  outputFile << " Primary ";
3095  outputFile.width(7);
3096  outputFile.precision(3);
3097  outputFile << frac_b_twoVertex_lTrack_Primary << " " ;
3098  outputFile.width(7);
3099  outputFile.precision(3);
3100  outputFile << frac_b_twoVertex_lTrack_Secondary << " " ;
3101  outputFile.width(7);
3102  outputFile.precision(3);
3103  outputFile << frac_b_twoVertex_lTrack_Isolated << " " ;
3104  outputFile.width(7);
3105  outputFile.precision(3);
3106  outputFile << frac_b_threeVertex_lTrack_Primary << " " ;
3107  outputFile.width(7);
3108  outputFile.precision(3);
3109  outputFile << frac_b_threeVertex_lTrack_Secondary << " " ;
3110  outputFile.width(7);
3111  outputFile.precision(3);
3112  outputFile << frac_b_threeVertex_lTrack_Tertiary << " " ;
3113  outputFile.width(7);
3114  outputFile.precision(3);
3115  outputFile << frac_b_threeVertex_lTrack_Isolated << " " ;
3116  outputFile << std::endl;
3117 
3118  outputFile << " B decay ";
3119  outputFile.width(7);
3120  outputFile.precision(3);
3121  outputFile << frac_b_twoVertex_bTrack_Primary << " " ;
3122  outputFile.width(7);
3123  outputFile.precision(3);
3124  outputFile << frac_b_twoVertex_bTrack_Secondary << " " ;
3125  outputFile.width(7);
3126  outputFile.precision(3);
3127  outputFile << frac_b_twoVertex_bTrack_Isolated << " " ;
3128  outputFile.width(7);
3129  outputFile.precision(3);
3130  outputFile << frac_b_threeVertex_bTrack_Primary << " " ;
3131  outputFile.width(7);
3132  outputFile.precision(3);
3133  outputFile << frac_b_threeVertex_bTrack_Secondary << " " ;
3134  outputFile.width(7);
3135  outputFile.precision(3);
3136  outputFile << frac_b_threeVertex_bTrack_Tertiary << " " ;
3137  outputFile.width(7);
3138  outputFile.precision(3);
3139  outputFile << frac_b_threeVertex_bTrack_Isolated << " " ;
3140  outputFile << std::endl;
3141 
3142  outputFile << " D decay ";
3143  outputFile.width(7);
3144  outputFile.precision(3);
3145  outputFile << frac_b_twoVertex_cTrack_Primary << " " ;
3146  outputFile.width(7);
3147  outputFile.precision(3);
3148  outputFile << frac_b_twoVertex_cTrack_Secondary << " " ;
3149  outputFile.width(7);
3150  outputFile.precision(3);
3151  outputFile << frac_b_twoVertex_cTrack_Isolated << " " ;
3152  outputFile.width(7);
3153  outputFile.precision(3);
3154  outputFile << frac_b_threeVertex_cTrack_Primary << " " ;
3155  outputFile.width(7);
3156  outputFile.precision(3);
3157  outputFile << frac_b_threeVertex_cTrack_Secondary << " " ;
3158  outputFile.width(7);
3159  outputFile.precision(3);
3160  outputFile << frac_b_threeVertex_cTrack_Tertiary << " " ;
3161  outputFile.width(7);
3162  outputFile.precision(3);
3163  outputFile << frac_b_threeVertex_cTrack_Isolated << " " ;
3164  outputFile << std::endl;
3165  outputFile << " ----------------------------------------------------------------------------" << std::endl;
3166  outputFile << " All above ";
3167  outputFile.width(7);
3168  outputFile.precision(3);
3169  outputFile << frac_b_twoVertex_Primary << " " ;
3170  outputFile.width(7);
3171  outputFile.precision(3);
3172  outputFile << frac_b_twoVertex_Secondary << " " ;
3173  outputFile.width(7);
3174  outputFile.precision(3);
3175  outputFile << frac_b_twoVertex_Isolated << " " ;
3176  outputFile.width(7);
3177  outputFile.precision(3);
3178  outputFile << frac_b_threeVertex_Primary << " " ;
3179  outputFile.width(7);
3180  outputFile.precision(3);
3181  outputFile << frac_b_threeVertex_Secondary << " " ;
3182  outputFile.width(7);
3183  outputFile.precision(3);
3184  outputFile << frac_b_threeVertex_Tertiary << " " ;
3185  outputFile.width(7);
3186  outputFile.precision(3);
3187  outputFile << frac_b_threeVertex_Isolated << " " ;
3188  outputFile << std::endl;
3189  outputFile << " ----------------------------------------------------------------------------" << std::endl;
3190  outputFile << std::endl;
3191  outputFile << std::endl;
3192 
3193 
3194  outputFile << std::endl;
3195  outputFile << " ---------------------------------------------------------------------------" << std::endl;
3196  outputFile << " C-JETS " << std::endl;
3197  outputFile << " ---------------------------------------------------------------------------" << std::endl;
3198  outputFile << " Purity of reconstructed track-vertex association (%) " << std::endl;
3199  outputFile << " ---------------------------------------------------------------------------" << std::endl;
3200  outputFile << " Monte Carlo Reconstructed track-vertex association " << std::endl;
3201  outputFile << " track origin Two-vertex case Three-vertex case " << std::endl;
3202  outputFile << " ---------------------------------------------------------------------------" << std::endl;
3203  outputFile << " Pri. Sec. Iso. Pri. Sec. Ter. Iso. " << std::endl;
3204  outputFile << " ---------------------------------------------------------------------------" << std::endl;
3205  outputFile << " Primary ";
3206  outputFile.width(7);
3207  outputFile.precision(3);
3208  outputFile << frac_c_twoVertex_lTrack_Primary << " " ;
3209  outputFile.width(7);
3210  outputFile.precision(3);
3211  outputFile << frac_c_twoVertex_lTrack_Secondary << " " ;
3212  outputFile.width(7);
3213  outputFile.precision(3);
3214  outputFile << frac_c_twoVertex_lTrack_Isolated << " " ;
3215  outputFile.width(7);
3216  outputFile.precision(3);
3217  outputFile << frac_c_threeVertex_lTrack_Primary << " " ;
3218  outputFile.width(7);
3219  outputFile.precision(3);
3220  outputFile << frac_c_threeVertex_lTrack_Secondary << " " ;
3221  outputFile.width(7);
3222  outputFile.precision(3);
3223  outputFile << frac_c_threeVertex_lTrack_Tertiary << " " ;
3224  outputFile.width(7);
3225  outputFile.precision(3);
3226  outputFile << frac_c_threeVertex_lTrack_Isolated << " " ;
3227  outputFile << std::endl;
3228 
3229 //check// outputFile << " B decay ";
3230 //check// outputFile.width(7);
3231 //check// outputFile.precision(3);
3232 //check// outputFile << frac_c_twoVertex_bTrack_Primary << " " ;
3233 //check// outputFile.width(7);
3234 //check// outputFile.precision(3);
3235 //check// outputFile << frac_c_twoVertex_bTrack_Secondary << " " ;
3236 //check// outputFile.width(7);
3237 //check// outputFile.precision(3);
3238 //check// outputFile << frac_c_twoVertex_bTrack_Isolated << " " ;
3239 //check// outputFile.width(7);
3240 //check// outputFile.precision(3);
3241 //check// outputFile << frac_c_threeVertex_bTrack_Primary << " " ;
3242 //check// outputFile.width(7);
3243 //check// outputFile.precision(3);
3244 //check// outputFile << frac_c_threeVertex_bTrack_Secondary << " " ;
3245 //check// outputFile.width(7);
3246 //check// outputFile.precision(3);
3247 //check// outputFile << frac_c_threeVertex_bTrack_Tertiary << " " ;
3248 //check// outputFile.width(7);
3249 //check// outputFile.precision(3);
3250 //check// outputFile << frac_c_threeVertex_bTrack_Isolated << " " ;
3251 //check// outputFile << std::endl;
3252 
3253  outputFile << " D decay ";
3254  outputFile.width(7);
3255  outputFile.precision(3);
3256  outputFile << frac_c_twoVertex_cTrack_Primary << " " ;
3257  outputFile.width(7);
3258  outputFile.precision(3);
3259  outputFile << frac_c_twoVertex_cTrack_Secondary << " " ;
3260  outputFile.width(7);
3261  outputFile.precision(3);
3262  outputFile << frac_c_twoVertex_cTrack_Isolated << " " ;
3263  outputFile.width(7);
3264  outputFile.precision(3);
3265  outputFile << frac_c_threeVertex_cTrack_Primary << " " ;
3266  outputFile.width(7);
3267  outputFile.precision(3);
3268  outputFile << frac_c_threeVertex_cTrack_Secondary << " " ;
3269  outputFile.width(7);
3270  outputFile.precision(3);
3271  outputFile << frac_c_threeVertex_cTrack_Tertiary << " " ;
3272  outputFile.width(7);
3273  outputFile.precision(3);
3274  outputFile << frac_c_threeVertex_cTrack_Isolated << " " ;
3275  outputFile << std::endl;
3276  outputFile << " ----------------------------------------------------------------------------" << std::endl;
3277  outputFile << " All above ";
3278  outputFile.width(7);
3279  outputFile.precision(3);
3280  outputFile << frac_c_twoVertex_Primary << " " ;
3281  outputFile.width(7);
3282  outputFile.precision(3);
3283  outputFile << frac_c_twoVertex_Secondary << " " ;
3284  outputFile.width(7);
3285  outputFile.precision(3);
3286  outputFile << frac_c_twoVertex_Isolated << " " ;
3287  outputFile.width(7);
3288  outputFile.precision(3);
3289  outputFile << frac_c_threeVertex_Primary << " " ;
3290  outputFile.width(7);
3291  outputFile.precision(3);
3292  outputFile << frac_c_threeVertex_Secondary << " " ;
3293  outputFile.width(7);
3294  outputFile.precision(3);
3295  outputFile << frac_c_threeVertex_Tertiary << " " ;
3296  outputFile.width(7);
3297  outputFile.precision(3);
3298  outputFile << frac_c_threeVertex_Isolated << " " ;
3299  outputFile << std::endl;
3300  outputFile << " ----------------------------------------------------------------------------" << std::endl;
3301  outputFile << std::endl;
3302  outputFile << std::endl;
3303 
3304 
3305 
3306 }
3307 
3308 void LCFIAIDAPlotProcessor::PrintNNOutput(){
3309 
3310  //if there is a _PurityEfficiencyOutputFile string defined use that as the output stream, if not use std::cout
3311  std::filebuf* fb = new std::filebuf;
3312 
3313  std::ostream outputFile( (!_PurityEfficiencyOutputFile.empty()) ?
3314  fb->open(_PurityEfficiencyOutputFile.c_str(),
3315  std::ios_base::out|std::ios_base::trunc):
3316  std::cout.rdbuf());
3317 
3318  if (outputFile.rdbuf() != fb)
3319  {
3320  delete fb;
3321  std::cerr << " In " << __FILE__ << "(" << __LINE__ << "): Unable to open output file " << _PurityEfficiencyOutputFile << "! Redirecting output to standard out." << std::endl;
3322  outputFile << std::endl;
3323  }
3324 
3325  for (unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection)
3326  {
3327  outputFile << "\n\nRESULTS for " << iTagCollection << "-th Flavour Tag Collection named " << _FlavourTagCollectionNames[iTagCollection] << std::endl;
3328  outputFile << " any number of ZVTOP vertices found " ;
3329  outputFile << " N(b) = " ;
3330  outputFile.width(10);
3331  outputFile << _pBJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->sumAllBinHeights();
3332  outputFile << " N(c) = " ;
3333  outputFile.width(10);
3334  outputFile << _pCJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->sumAllBinHeights();
3335  outputFile << " N(light) = ";
3336  outputFile.width(10);
3337  outputFile << _pLightJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->sumAllBinHeights();
3338  outputFile << std::endl << std::endl;
3339 
3340  for (unsigned int iVertexCat=1; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
3341 
3342  std::string nvname = _VertexCatNames[iVertexCat];
3343 
3344 
3345 
3346  outputFile << " ";
3347  outputFile.width(2);
3348  outputFile << iVertexCat;
3349  outputFile << (
3350  (iVertexCat == 1) ? " ZVTOP vertex found " : " ZVTOP vertices found ") ;
3351 
3352  outputFile << " N(b) = " ;
3353  outputFile.width(10);
3354  outputFile << _pBJetBTagIntegral[iTagCollection][nvname]->sumAllBinHeights();
3355  outputFile << " N(c) = " ;
3356  outputFile.width(10);
3357  outputFile << _pCJetBTagIntegral[iTagCollection][nvname]->sumAllBinHeights();
3358  outputFile << " N(light) = ";
3359  outputFile.width(10);
3360  outputFile << _pLightJetBTagIntegral[iTagCollection][nvname]->sumAllBinHeights();
3361  outputFile << std::endl << std::endl;
3362  }
3363 
3364 
3365 
3366 
3367  outputFile.precision(5);
3368  outputFile.setf(std::ios::fixed,std::ios::floatfield);
3369 
3370 
3371  for (unsigned int iVertexCat=1; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
3372 
3373  std::string nvname = _VertexCatNames[iVertexCat];
3374 
3375  outputFile << "numbers of jets in cuts for ";
3376  outputFile.width(2);
3377  outputFile << iVertexCat;
3378  outputFile << " ZVTOP vertices found" << std::endl;
3379  outputFile << "cut b-tag b b-tag other c-tag c c-tag other c-tagbb c c-tagbb other" << std::endl;
3380 
3381  int numberOfBins=_pBJetBTagIntegral[iTagCollection][nvname]->axis().bins();
3382 
3383  for (int binNumber = 0; binNumber < numberOfBins ; binNumber++ ) {
3384 
3385  outputFile.width(5);
3386  outputFile.precision(3);
3387  outputFile << _pBJetBTagIntegral[iTagCollection][nvname]->axis().binUpperEdge(binNumber) << " ";
3388  outputFile.width(10);
3389  outputFile << int(_pBJetBTagIntegral[iTagCollection][nvname]->binHeight(binNumber)) << " ";
3390  outputFile.width(10);
3391  outputFile << int(_pCJetBTagIntegral[iTagCollection][nvname]->binHeight(binNumber)+_pLightJetBTagIntegral[iTagCollection][nvname]->binHeight(binNumber))<< " ";
3392  outputFile.width(10);
3393  outputFile << int(_pCJetCTagIntegral[iTagCollection][nvname]->binHeight(binNumber)) << " ";
3394  outputFile.width(10);
3395  outputFile << int(_pBJetCTagIntegral[iTagCollection][nvname]->binHeight(binNumber)+_pLightJetCTagIntegral[iTagCollection][nvname]->binHeight(binNumber))<< " ";
3396  outputFile.width(10);
3397  outputFile << int(_pCJetBCTagIntegral[iTagCollection][nvname]->binHeight(binNumber)) << " ";
3398  outputFile.width(10);
3399  outputFile << int(_pBJetBCTagIntegral[iTagCollection][nvname]->binHeight(binNumber)+_pLightJetBCTagIntegral[iTagCollection][nvname]->binHeight(binNumber))<< " ";
3400  outputFile << std::endl;
3401 
3402  }
3403  outputFile << std::endl;
3404  }
3405 
3406  outputFile << "numbers of jets in cuts summed (for any number of vertices)" << std::endl;
3407  outputFile << "cut b-tag b b-tag other c-tag c c-tag other c-tagbb c c-tagbb other" << std::endl;
3408 
3409  int numberOfBins=_pBJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->axis().bins();
3410 
3411  for (int binNumber = 0; binNumber < numberOfBins ; binNumber++ ) {
3412 
3413  outputFile.width(5);
3414  outputFile.precision(3);
3415  outputFile << _pBJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->axis().binUpperEdge(binNumber) << " ";
3416  outputFile.width(10);
3417  outputFile << int(_pBJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber)) << " ";
3418  outputFile.width(10);
3419  outputFile << int(_pCJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber)+_pLightJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber))<< " ";
3420  outputFile.width(10);
3421  outputFile << int(_pCJetCTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber)) << " ";
3422  outputFile.width(10);
3423  outputFile << int(_pBJetCTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber)+_pLightJetCTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber))<< " ";
3424  outputFile.width(10);
3425  outputFile << int(_pCJetBCTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber)) << " ";
3426  outputFile.width(10);
3427  outputFile << int(_pBJetBCTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber)+_pLightJetBCTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber))<< " ";
3428  outputFile << std::endl;
3429  }
3430  outputFile << std::endl;
3431  outputFile << std::endl;
3432 
3433  }
3434 
3435  for (unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection)
3436  {
3437 
3438  int numberOfBins=_pBJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->axis().bins();
3439 
3440  outputFile << "\n\nRESULTS for " << iTagCollection << "-th Flavour Tag Collection named " << _FlavourTagCollectionNames[iTagCollection] << std::endl;
3441 
3442  for (unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
3443 
3444  std::string nvname = _VertexCatNames[iVertexCat];
3445 
3446  std::vector<double> bJetBTagEfficiency;
3447  std::vector<double> bJetBTagEfficiencyError;
3448  std::vector<double> cJetCTagEfficiency;
3449  std::vector<double> cJetCTagEfficiencyError;
3450  std::vector<double> cJetBCTagEfficiency;
3451  std::vector<double> cJetBCTagEfficiencyError;
3452 
3453  CalculateTagEfficiency( _pBJetBTag[iTagCollection][nvname], bJetBTagEfficiency, bJetBTagEfficiencyError);
3454  CalculateTagEfficiency( _pCJetCTag[iTagCollection][nvname], cJetCTagEfficiency, cJetCTagEfficiencyError);
3455  CalculateTagEfficiency( _pCJetBCTag[iTagCollection][nvname], cJetBCTagEfficiency, cJetBCTagEfficiencyError);
3456 
3457  std::vector<double> bJetBTagPurity;
3458  std::vector<double> bJetBTagPurityError;
3459  std::vector<double> cJetCTagPurity;
3460  std::vector<double> cJetCTagPurityError;
3461  std::vector<double> cJetBCTagPurity;
3462  std::vector<double> cJetBCTagPurityError;
3463 
3464  CalculateTagPurity( _pBJetBTag[iTagCollection][nvname], _pBTagBackgroundValues[iTagCollection][nvname], bJetBTagPurity, bJetBTagPurityError);
3465  CalculateTagPurity( _pCJetCTag[iTagCollection][nvname], _pCTagBackgroundValues[iTagCollection][nvname], cJetCTagPurity, cJetCTagPurityError);
3466  CalculateTagPurity( _pCJetBCTag[iTagCollection][nvname], _pBJetBCTag[iTagCollection][nvname], cJetBCTagPurity, cJetBCTagPurityError);
3467 
3468  outputFile << "\n Purity-Efficiencies Values for ";
3469  if (iVertexCat==0) {
3470  outputFile << " any number of ZVTOP vertices found " << std::endl;
3471  } else {
3472  outputFile << iVertexCat;
3473  outputFile << (
3474  (iVertexCat == 1) ? " ZVTOP vertex found " : " ZVTOP vertices found ") ;
3475  }
3476  outputFile << "\n-----------------------------------------------------------------------------------------------" << std::endl;
3477  outputFile << "cut value efficiency(b) purity(b) efficiency(c) purity(c) efficiency(bc) purity(bc) " << std::endl;
3478 
3479  for (int binNumber = 0; binNumber < numberOfBins ; binNumber++ ) {
3480  outputFile.width(5);
3481  outputFile.precision(3);
3482  outputFile << _pBJetBTagIntegral[iTagCollection][nvname]->axis().binUpperEdge(binNumber) << " ";
3483 
3484  outputFile << bJetBTagEfficiency[binNumber];
3485  outputFile << " +/- ";
3486  outputFile << bJetBTagEfficiencyError[binNumber];
3487  outputFile << " ";
3488  outputFile << bJetBTagPurity[binNumber];
3489  outputFile << " +/- ";
3490  outputFile << bJetBTagPurityError[binNumber];
3491  outputFile << " ";
3492 
3493  outputFile << cJetCTagEfficiency[binNumber];
3494  outputFile << " +/- ";
3495  outputFile << cJetCTagEfficiencyError[binNumber];
3496  outputFile << " ";
3497  outputFile << cJetCTagPurity[binNumber];
3498  outputFile << " +/- ";
3499  outputFile << cJetCTagPurityError[binNumber];
3500  outputFile << " ";
3501 
3502  outputFile << cJetBCTagEfficiency[binNumber];
3503  outputFile << " +/- ";
3504  outputFile << cJetBCTagEfficiencyError[binNumber];
3505  outputFile << " ";
3506  outputFile << cJetBCTagPurity[binNumber];
3507  outputFile << " +/- ";
3508  outputFile << cJetBCTagPurityError[binNumber];
3509  outputFile << " ";
3510 
3511  outputFile << std::endl;
3512  }
3513  }
3514 
3515  outputFile << std::endl;
3516  outputFile << std::endl;
3517 
3518  }
3519 
3520 }
3521 
3522 void LCFIAIDAPlotProcessor::InternalVectorInitialisation()
3523 {
3524  //sets the size of some vectors and initialises some counters to zero
3525 
3526  _cJet_truePlus2_angle.resize(N_JETANGLE_BINS);
3527  _cJet_truePlus_angle.resize(N_JETANGLE_BINS);
3528  _cJet_trueNeut_angle.resize(N_JETANGLE_BINS);
3529  _cJet_trueMinus_angle.resize(N_JETANGLE_BINS);
3530  _cJet_trueMinus2_angle.resize(N_JETANGLE_BINS);
3531 
3532  _cJet_truePlus2_recoPlus_angle.resize(N_JETANGLE_BINS);
3533  _cJet_truePlus2_recoNeut_angle.resize(N_JETANGLE_BINS);
3534  _cJet_truePlus2_recoMinus_angle.resize(N_JETANGLE_BINS);
3535  _cJet_truePlus_recoPlus_angle.resize(N_JETANGLE_BINS);
3536  _cJet_truePlus_recoNeut_angle.resize(N_JETANGLE_BINS);
3537  _cJet_truePlus_recoMinus_angle.resize(N_JETANGLE_BINS);
3538  _cJet_trueNeut_recoPlus_angle.resize(N_JETANGLE_BINS);
3539  _cJet_trueNeut_recoNeut_angle.resize(N_JETANGLE_BINS);
3540  _cJet_trueNeut_recoMinus_angle.resize(N_JETANGLE_BINS);
3541  _cJet_trueMinus_recoPlus_angle.resize(N_JETANGLE_BINS);
3542  _cJet_trueMinus_recoNeut_angle.resize(N_JETANGLE_BINS);
3543  _cJet_trueMinus_recoMinus_angle.resize(N_JETANGLE_BINS);
3544  _cJet_trueMinus2_recoPlus_angle.resize(N_JETANGLE_BINS);
3545  _cJet_trueMinus2_recoNeut_angle.resize(N_JETANGLE_BINS);
3546  _cJet_trueMinus2_recoMinus_angle.resize(N_JETANGLE_BINS);
3547 
3548  _bJet_truePlus2_angle.resize(N_JETANGLE_BINS);
3549  _bJet_truePlus_angle.resize(N_JETANGLE_BINS);
3550  _bJet_trueNeut_angle.resize(N_JETANGLE_BINS);
3551  _bJet_trueMinus_angle.resize(N_JETANGLE_BINS);
3552  _bJet_trueMinus2_angle.resize(N_JETANGLE_BINS);
3553  _bJet_truePlus2_recoPlus_angle.resize(N_JETANGLE_BINS);
3554  _bJet_truePlus2_recoNeut_angle.resize(N_JETANGLE_BINS);
3555  _bJet_truePlus2_recoMinus_angle.resize(N_JETANGLE_BINS);
3556  _bJet_truePlus_recoPlus_angle.resize(N_JETANGLE_BINS);
3557  _bJet_truePlus_recoNeut_angle.resize(N_JETANGLE_BINS);
3558  _bJet_truePlus_recoMinus_angle.resize(N_JETANGLE_BINS);
3559  _bJet_trueNeut_recoPlus_angle.resize(N_JETANGLE_BINS);
3560  _bJet_trueNeut_recoNeut_angle.resize(N_JETANGLE_BINS);
3561  _bJet_trueNeut_recoMinus_angle.resize(N_JETANGLE_BINS);
3562  _bJet_trueMinus_recoPlus_angle.resize(N_JETANGLE_BINS);
3563  _bJet_trueMinus_recoNeut_angle.resize(N_JETANGLE_BINS);
3564  _bJet_trueMinus_recoMinus_angle.resize(N_JETANGLE_BINS);
3565  _bJet_trueMinus2_recoPlus_angle.resize(N_JETANGLE_BINS);
3566  _bJet_trueMinus2_recoNeut_angle.resize(N_JETANGLE_BINS);
3567  _bJet_trueMinus2_recoMinus_angle.resize(N_JETANGLE_BINS);
3568 
3569 
3570  for (int j = 0 ; j < N_JETANGLE_BINS ; j++) {
3571  _cJet_truePlus2_recoPlus_angle[j]=0;
3572  _cJet_truePlus2_recoNeut_angle[j]=0;
3573  _cJet_truePlus2_recoMinus_angle[j]=0;
3574  _cJet_truePlus_recoPlus_angle[j]=0;
3575  _cJet_truePlus_recoNeut_angle[j]=0;
3576  _cJet_truePlus_recoMinus_angle[j]=0;
3577  _cJet_trueNeut_recoPlus_angle[j]=0;
3578  _cJet_trueNeut_recoNeut_angle[j]=0;
3579  _cJet_trueNeut_recoMinus_angle[j]=0;
3580  _cJet_trueMinus_recoPlus_angle[j]=0;
3581  _cJet_trueMinus_recoNeut_angle[j]=0;
3582  _cJet_trueMinus_recoMinus_angle[j]=0;
3583  _cJet_trueMinus2_recoPlus_angle[j]=0;
3584  _cJet_trueMinus2_recoNeut_angle[j]=0;
3585  _cJet_trueMinus2_recoMinus_angle[j]=0;
3586 
3587  _bJet_truePlus2_angle[j]=0;
3588  _bJet_truePlus_angle[j]=0;
3589  _bJet_trueNeut_angle[j]=0;
3590  _bJet_trueMinus_angle[j]=0;
3591  _bJet_trueMinus2_angle[j]=0;
3592  _bJet_truePlus2_recoPlus_angle[j]=0;
3593  _bJet_truePlus2_recoNeut_angle[j]=0;
3594  _bJet_truePlus2_recoMinus_angle[j]=0;
3595  _bJet_truePlus_recoPlus_angle[j]=0;
3596  _bJet_truePlus_recoNeut_angle[j]=0;
3597  _bJet_truePlus_recoMinus_angle[j]=0;
3598  _bJet_trueNeut_recoPlus_angle[j]=0;
3599  _bJet_trueNeut_recoNeut_angle[j]=0;
3600  _bJet_trueNeut_recoMinus_angle[j]=0;
3601  _bJet_trueMinus_recoPlus_angle[j]=0;
3602  _bJet_trueMinus_recoNeut_angle[j]=0;
3603  _bJet_trueMinus_recoMinus_angle[j]=0;
3604  _bJet_trueMinus2_recoPlus_angle[j]=0;
3605  _bJet_trueMinus2_recoNeut_angle[j]=0;
3606  _bJet_trueMinus2_recoMinus_angle[j]=0;
3607  }
3608 
3609 
3610  _cJet_truePlus2=0;
3611  _cJet_truePlus=0;
3612  _cJet_trueNeut=0;
3613  _cJet_trueMinus=0;
3614  _cJet_trueMinus2=0;
3615  _cJet_truePlus2_recoPlus=0;
3616  _cJet_truePlus2_recoNeut=0;
3617  _cJet_truePlus2_recoMinus=0;
3618  _cJet_truePlus_recoPlus=0;
3619  _cJet_truePlus_recoNeut=0;
3620  _cJet_truePlus_recoMinus=0;
3621  _cJet_trueNeut_recoPlus=0;
3622  _cJet_trueNeut_recoNeut=0;
3623  _cJet_trueNeut_recoMinus=0;
3624  _cJet_trueMinus_recoPlus=0;
3625  _cJet_trueMinus_recoNeut=0;
3626  _cJet_trueMinus_recoMinus=0;
3627  _cJet_trueMinus2_recoPlus=0;
3628  _cJet_trueMinus2_recoNeut=0;
3629  _cJet_trueMinus2_recoMinus=0;
3630 
3631  _bJet_truePlus2=0;
3632  _bJet_truePlus=0;
3633  _bJet_trueNeut=0;
3634  _bJet_trueMinus=0;
3635  _bJet_trueMinus2=0;
3636  _bJet_truePlus2_recoPlus=0;
3637  _bJet_truePlus2_recoNeut=0;
3638  _bJet_truePlus2_recoMinus=0;
3639  _bJet_truePlus_recoPlus=0;
3640  _bJet_truePlus_recoNeut=0;
3641  _bJet_truePlus_recoMinus=0;
3642  _bJet_trueNeut_recoPlus=0;
3643  _bJet_trueNeut_recoNeut=0;
3644  _bJet_trueNeut_recoMinus=0;
3645  _bJet_trueMinus_recoPlus=0;
3646  _bJet_trueMinus_recoNeut=0;
3647  _bJet_trueMinus_recoMinus=0;
3648  _bJet_trueMinus2_recoPlus=0;
3649  _bJet_trueMinus2_recoNeut=0;
3650  _bJet_trueMinus2_recoMinus=0;
3651 
3652  _nb_twoVertex_bTrack_Primary=0;
3653  _nb_twoVertex_bTrack_Secondary=0;
3654  _nb_twoVertex_bTrack_Tertiary=0;
3655  _nb_twoVertex_bTrack_Isolated=0;
3656  _nb_twoVertex_cTrack_Primary=0;
3657  _nb_twoVertex_cTrack_Secondary=0;
3658  _nb_twoVertex_cTrack_Tertiary=0;
3659  _nb_twoVertex_cTrack_Isolated=0;
3660  _nb_twoVertex_lTrack_Primary=0;
3661  _nb_twoVertex_lTrack_Secondary=0;
3662  _nb_twoVertex_lTrack_Tertiary=0;
3663  _nb_twoVertex_lTrack_Isolated=0;
3664  _nb_threeVertex_bTrack_Primary=0;
3665  _nb_threeVertex_bTrack_Secondary=0;
3666  _nb_threeVertex_bTrack_Tertiary=0;
3667  _nb_threeVertex_bTrack_Isolated=0;
3668  _nb_threeVertex_cTrack_Primary=0;
3669  _nb_threeVertex_cTrack_Secondary=0;
3670  _nb_threeVertex_cTrack_Tertiary=0;
3671  _nb_threeVertex_cTrack_Isolated=0;
3672  _nb_threeVertex_lTrack_Primary=0;
3673  _nb_threeVertex_lTrack_Secondary=0;
3674  _nb_threeVertex_lTrack_Tertiary=0;
3675  _nb_threeVertex_lTrack_Isolated=0;
3676 
3677  _nb_threeVertex_Primary_noMCP=0;
3678  _nb_threeVertex_Secondary_noMCP=0;
3679  _nb_threeVertex_Tertiary_noMCP=0;
3680  _nb_threeVertex_Isolated_noMCP=0;
3681  _nb_twoVertex_Primary_noMCP=0;
3682  _nb_twoVertex_Secondary_noMCP=0;
3683  _nb_twoVertex_Tertiary_noMCP=0;
3684  _nb_twoVertex_Isolated_noMCP=0;
3685 
3686  _nc_twoVertex_bTrack_Primary=0;
3687  _nc_twoVertex_bTrack_Secondary=0;
3688  _nc_twoVertex_bTrack_Tertiary=0;
3689  _nc_twoVertex_bTrack_Isolated=0;
3690  _nc_twoVertex_cTrack_Primary=0;
3691  _nc_twoVertex_cTrack_Secondary=0;
3692  _nc_twoVertex_cTrack_Tertiary=0;
3693  _nc_twoVertex_cTrack_Isolated=0;
3694  _nc_twoVertex_lTrack_Primary=0;
3695  _nc_twoVertex_lTrack_Secondary=0;
3696  _nc_twoVertex_lTrack_Tertiary=0;
3697  _nc_twoVertex_lTrack_Isolated=0;
3698  _nc_threeVertex_bTrack_Primary=0;
3699  _nc_threeVertex_bTrack_Secondary=0;
3700  _nc_threeVertex_bTrack_Tertiary=0;
3701  _nc_threeVertex_bTrack_Isolated=0;
3702  _nc_threeVertex_cTrack_Primary=0;
3703  _nc_threeVertex_cTrack_Secondary=0;
3704  _nc_threeVertex_cTrack_Tertiary=0;
3705  _nc_threeVertex_cTrack_Isolated=0;
3706  _nc_threeVertex_lTrack_Primary=0;
3707  _nc_threeVertex_lTrack_Secondary=0;
3708  _nc_threeVertex_lTrack_Tertiary=0;
3709  _nc_threeVertex_lTrack_Isolated=0;
3710 
3711  _nc_threeVertex_Primary_noMCP=0;
3712  _nc_threeVertex_Secondary_noMCP=0;
3713  _nc_threeVertex_Tertiary_noMCP=0;
3714  _nc_threeVertex_Isolated_noMCP=0;
3715  _nc_twoVertex_Primary_noMCP=0;
3716  _nc_twoVertex_Secondary_noMCP=0;
3717  _nc_twoVertex_Tertiary_noMCP=0;
3718  _nc_twoVertex_Isolated_noMCP=0;
3719 
3720 
3721 
3722 }
3723 
3724 
3725 void LCFIAIDAPlotProcessor::FindTrueJetDecayLength( LCEvent* pEvent, unsigned int jetNumber, std::vector<double>& decaylengthvector,
3726  std::vector<double>& bjetdecaylengthvector, std::vector<double>& cjetdecaylengthvector)
3727 {
3728 
3729  //looks in the MC to find the longest decay length for the b and c jets
3730  //returns -99 for light jets
3731 
3732  int jettype = FindTrueJetType( pEvent, jetNumber );
3733 
3734  if (abs(jettype)!= C_JET && abs(jettype)!=B_JET)
3735  {
3736  return;
3737  } else {
3738 
3739  int pdgcode = FindTrueJetPDGCode( pEvent, jetNumber );
3740 
3741  TypesafeCollection<lcio::MCParticle> mcParticleCollection( pEvent, _MCParticleColName );
3742 
3743  if (!mcParticleCollection.is_valid())
3744  {
3745  std::cerr << "In " << __FILE__ << "(" << __LINE__ << "): for event " << pEvent->getEventNumber() << " in run " << pEvent->getRunNumber() << "mcParticleCollection is not valid" << std::endl;
3746  return;
3747  }
3748  else
3749  {
3750  for (int imc = 0 ; imc < mcParticleCollection.getNumberOfElements(); imc++ ) {
3751 
3752  lcio::MCParticle* myMCParticle = mcParticleCollection.getElementAt( imc );
3753 
3754  if (myMCParticle->getPDG() == pdgcode) {
3755 
3756  //v./ std::cout << "1: " << pdgcode << " " << GetPDGFlavour(pdgcode) << std::endl;
3757 
3758  //get decay length
3759  const double* vertex = myMCParticle->getVertex();
3760  const double* endpoint = myMCParticle->getEndpoint();
3761  double decay_length = CalculateDistance(vertex, endpoint);
3762 
3763  //if (decay_length>0) {
3764 
3765  decaylengthvector.push_back(decay_length);
3766  if (GetPDGFlavour(pdgcode) == C_JET ) cjetdecaylengthvector.push_back(decay_length);
3767  if (GetPDGFlavour(pdgcode) == B_JET ) bjetdecaylengthvector.push_back(decay_length);
3768  //}
3769  //get the daughters
3770  const MCParticleVec& daughters = myMCParticle->getDaughters();
3771 
3772  for ( MCParticleVec::const_iterator idaughter = daughters.begin(); idaughter!=daughters.end(); idaughter++)
3773  {
3774  int code = (*idaughter)->getPDG();
3775 
3776  //std::cout << "PDG code of daughter 1: " << code << std::endl;
3777  //check if they are still bs or ds
3778  int flavour = GetPDGFlavour(code);
3779  //std::cout << "Daughter 1 flavour is: " << flavour << std::endl;
3780  if (abs(flavour) == C_JET || abs(flavour) == B_JET) {
3781 
3782  //v./ std::cout << "2: " << code << " " << flavour << std::endl;
3783 
3784 
3785  //get decay length
3786  const double* vertex = (*idaughter)->getVertex();
3787  const double* endpoint = (*idaughter)->getEndpoint();
3788  double decay_length = CalculateDistance(vertex, endpoint);
3789 
3790  //if (decay_length>0) {
3791  decaylengthvector.push_back(decay_length);
3792  if (abs(flavour) == C_JET ) cjetdecaylengthvector.push_back(decay_length);
3793  if (abs(flavour) == B_JET ) bjetdecaylengthvector.push_back(decay_length);
3794 
3795 
3796  // }
3797 
3798 
3799  // std::cout << "Decay length 2 is: " << decay_length << std::endl;
3800 
3801  //look for daughters again - just in case...
3802 
3803  const MCParticleVec& daughters2 = (*idaughter)->getDaughters();
3804 
3805  for ( MCParticleVec::const_iterator idaughter2 = daughters2.begin(); idaughter2!=daughters2.end(); idaughter2++)
3806  {
3807  int code2 = (*idaughter2)->getPDG();
3808 
3809  //check if they are still bs or ds
3810  int flavour2 = GetPDGFlavour(code2);
3811 
3812  if (abs(flavour2) == C_JET || abs (flavour2) == B_JET) {
3813 
3814  //v./ std::cout << "3: " << code2 << " " << flavour2 << std::endl;
3815 
3816  //get decay length
3817  const double* vertex = (*idaughter2)->getVertex();
3818  const double* endpoint = (*idaughter2)->getEndpoint();
3819  double decay_length = CalculateDistance(vertex, endpoint);
3820 
3821  //if (decay_length>0) {
3822  decaylengthvector.push_back(decay_length);
3823  if (abs(flavour2) == C_JET ) cjetdecaylengthvector.push_back(decay_length);
3824  if (abs(flavour2) == B_JET ) bjetdecaylengthvector.push_back(decay_length);
3825  //}
3826 
3827  const MCParticleVec& daughters3 = (*idaughter2)->getDaughters();
3828 
3829  for ( MCParticleVec::const_iterator idaughter3 = daughters3.begin(); idaughter3!=daughters3.end(); idaughter3++)
3830  {
3831  int code3 = (*idaughter3)->getPDG();
3832 
3833  //check if they are still bs or ds
3834  int flavour3 = GetPDGFlavour(code3);
3835 
3836  if (abs(flavour3)==C_JET || abs(flavour3)==B_JET) {
3837 
3838  //v./ std::cout << "4: " << code3 << " " << flavour3 << std::endl;
3839 
3840  //get decay length
3841  const double* vertex = (*idaughter3)->getVertex();
3842  const double* endpoint = (*idaughter3)->getEndpoint();
3843  double decay_length = CalculateDistance(vertex, endpoint);
3844 
3845  //if (decay_length>0) {
3846  decaylengthvector.push_back(decay_length);
3847  if (abs(flavour3) == C_JET ) cjetdecaylengthvector.push_back(decay_length);
3848  if (abs(flavour3) == B_JET ) bjetdecaylengthvector.push_back(decay_length);
3849  //}
3850 
3851  const MCParticleVec& daughters4 = (*idaughter3)->getDaughters();
3852 
3853  for ( MCParticleVec::const_iterator idaughter4 = daughters4.begin(); idaughter4!=daughters4.end(); idaughter4++)
3854  {
3855  int code4 = (*idaughter4)->getPDG();
3856 
3857  //check if they are still bs or ds
3858  int flavour4 = GetPDGFlavour(code4);
3859 
3860  if (abs(flavour4)==C_JET || abs(flavour4)==B_JET) {
3861 
3862  //v./ std::cout << "5: " << code4 << " " << flavour4 << " " << std::endl;
3863 
3864  //get decay length
3865  const double* vertex = (*idaughter4)->getVertex();
3866  const double* endpoint = (*idaughter4)->getEndpoint();
3867  double decay_length = CalculateDistance(vertex, endpoint);
3868  //if (decay_length>0) {
3869  decaylengthvector.push_back(decay_length);
3870  if (abs(flavour4) == C_JET ) cjetdecaylengthvector.push_back(decay_length);
3871  if (abs(flavour4) == B_JET ) bjetdecaylengthvector.push_back(decay_length);
3872  //}
3873  }
3874  }
3875 
3876  }
3877  }
3878 
3879 
3880 
3881  }//b or c
3882 
3883  }//loop over daughter's daugthers
3884 
3885  }//daughter is b or c
3886  }//loop over daughters
3887  }
3888  }
3889  }
3890  }
3891 
3892  return;
3893 }
3894 
3895 void LCFIAIDAPlotProcessor::FindTrueJetDecayLength2( LCEvent* pEvent, unsigned int jetNumber, double& bjetdecaylength, double& cjetdecaylength)
3896 {
3897  double b_decay_length0 = -1.;
3898  double c_decay_length0 = -1.;
3899 
3900  //looks in the MC to find the longest decay length for the b and c jets
3901 
3902  int jettype = FindTrueJetType( pEvent, jetNumber );
3903 
3904  if (abs(jettype)!=B_JET && abs(jettype)!=C_JET )
3905  {
3906  return;
3907  } else {
3908 
3909  int pdgcode = FindTrueJetPDGCode( pEvent, jetNumber );
3910 
3911  TypesafeCollection<lcio::MCParticle> mcParticleCollection( pEvent, _MCParticleColName );
3912 
3913  if (!mcParticleCollection.is_valid())
3914  {
3915  std::cerr << __FILE__ << "(" << __LINE__ << "): MCParticle collection " << _MCParticleColName << " is not valid"
3916  << " for event " << pEvent->getEventNumber() << " in run " << pEvent->getRunNumber() << std::endl;
3917  return;
3918  }
3919  else
3920  {
3921  for (int imc = 0 ; imc < mcParticleCollection.getNumberOfElements(); imc++ ) {
3922 
3923  lcio::MCParticle* myMCParticle = mcParticleCollection.getElementAt( imc );
3924  if (myMCParticle->getPDG() == pdgcode) {
3925 
3926 
3927  std::vector<double> BJetDecayLengthVector;
3928  std::vector<double> CJetDecayLengthVector;
3929  std::vector<int> BJetDecayVector;
3930  std::vector<int> CJetDecayVector;
3931  //std::cout << "1: " << pdgcode << " " << GetPDGFlavour(pdgcode) << std::endl;
3932 
3933  int code0 = pdgcode;
3934  int flavour0 = GetPDGFlavour(code0);
3935 
3936 
3937  //get start vertex
3938  const double* vertex0 = myMCParticle->getVertex();
3939 
3940  //get the daughters
3941  const MCParticleVec& daughters = myMCParticle->getDaughters();
3942 
3943  for ( MCParticleVec::const_iterator idaughter = daughters.begin(); idaughter!=daughters.end(); idaughter++)
3944  {
3945  int code1 = (*idaughter)->getPDG();
3946  int flavour1 = GetPDGFlavour(code1);
3947 
3948  if (abs(flavour1) != B_JET && abs(flavour0) == B_JET) {
3949 
3950  const double* vertex1 = (*idaughter)->getVertex();
3951  double decay_length = CalculateDistance(vertex0, vertex1);
3952  BJetDecayLengthVector.push_back(decay_length);
3953  }
3954 
3955  if (abs(flavour1) != C_JET && abs(flavour0) == C_JET) {
3956  const double* vertex1 = (*idaughter)->getVertex();
3957  double decay_length = CalculateDistance(vertex0, vertex1);
3958  CJetDecayLengthVector.push_back(decay_length);
3959  }
3960 
3961  //look for daughters again - just in case...
3962 
3963  const MCParticleVec& daughters2 = (*idaughter)->getDaughters();
3964 
3965  for ( MCParticleVec::const_iterator idaughter2 = daughters2.begin(); idaughter2!=daughters2.end(); idaughter2++)
3966  {
3967  int code2 = (*idaughter2)->getPDG();
3968  int flavour2 = GetPDGFlavour(code2);
3969 
3970  if (abs(flavour2) != B_JET && abs(flavour1) == B_JET) {
3971  const double* vertex2 = (*idaughter2)->getVertex();
3972  double decay_length = CalculateDistance(vertex0, vertex2);
3973  BJetDecayLengthVector.push_back(decay_length);
3974  }
3975 
3976  if (abs(flavour2) != C_JET && abs(flavour1) == C_JET) {
3977  const double* vertex2 = (*idaughter2)->getVertex();
3978  double decay_length = CalculateDistance(vertex0, vertex2);
3979  CJetDecayLengthVector.push_back(decay_length);
3980  }
3981 
3982  const MCParticleVec& daughters3 = (*idaughter2)->getDaughters();
3983 
3984  for ( MCParticleVec::const_iterator idaughter3 = daughters3.begin(); idaughter3!=daughters3.end(); idaughter3++)
3985  {
3986  int code3 = (*idaughter3)->getPDG();
3987  int flavour3 = GetPDGFlavour(code3);
3988 
3989  if (abs(flavour3) != B_JET && abs(flavour2) == B_JET) {
3990  const double* vertex3 = (*idaughter3)->getVertex();
3991  double decay_length = CalculateDistance(vertex0, vertex3);
3992  BJetDecayLengthVector.push_back(decay_length);
3993  }
3994 
3995  if (abs(flavour3) != C_JET && abs(flavour2) == C_JET) {
3996  const double* vertex3 = (*idaughter3)->getVertex();
3997  double decay_length = CalculateDistance(vertex0, vertex3);
3998  CJetDecayLengthVector.push_back(decay_length);
3999  }
4000 
4001  const MCParticleVec& daughters4 = (*idaughter3)->getDaughters();
4002 
4003  for ( MCParticleVec::const_iterator idaughter4 = daughters4.begin(); idaughter4!=daughters4.end(); idaughter4++)
4004  {
4005  int code4 = (*idaughter4)->getPDG();
4006  int flavour4 = GetPDGFlavour(code4);
4007 
4008  if (abs(flavour4) != B_JET && abs(flavour3) == B_JET) {
4009  const double* vertex4 = (*idaughter4)->getVertex();
4010  double decay_length = CalculateDistance(vertex0, vertex4);
4011  BJetDecayLengthVector.push_back(decay_length);
4012  }
4013  if (abs(flavour4) != C_JET && abs(flavour3) == C_JET) {
4014  const double* vertex4 = (*idaughter4)->getVertex();
4015  double decay_length = CalculateDistance(vertex0, vertex4);
4016  CJetDecayLengthVector.push_back(decay_length);
4017  }
4018 
4019 
4020  const MCParticleVec& daughters5 = (*idaughter4)->getDaughters();
4021 
4022  for ( MCParticleVec::const_iterator idaughter5 = daughters5.begin(); idaughter5!=daughters5.end(); idaughter5++)
4023  {
4024  int code5 = (*idaughter5)->getPDG();
4025  int flavour5 = GetPDGFlavour(code5);
4026 
4027  if (abs(flavour5) != B_JET && abs(flavour4) == B_JET) {
4028  const double* vertex5 = (*idaughter5)->getVertex();
4029  double decay_length = CalculateDistance(vertex0, vertex5);
4030  BJetDecayLengthVector.push_back(decay_length);
4031  }
4032  if (abs(flavour5) != C_JET && abs(flavour4) == C_JET) {
4033  const double* vertex5 = (*idaughter5)->getVertex();
4034  double decay_length = CalculateDistance(vertex0, vertex5);
4035  CJetDecayLengthVector.push_back(decay_length);
4036  }
4037 
4038  //std::cout << code0 << " " << code1 << " " << code2 << " " << code3 << " " << code4 << " " << code5 << std::endl;
4039 
4040  }
4041  }
4042  }
4043  }
4044  }
4045  //most of the time we will only have one entry per vector
4046  for (std::vector<double>::const_iterator iter = BJetDecayLengthVector.begin(); iter != BJetDecayLengthVector.end() ; iter++)
4047  {
4048  double decay_length = *iter;
4049  if (decay_length > b_decay_length0 ) b_decay_length0 = decay_length;
4050  }
4051 
4052  for (std::vector<double>::const_iterator iter = CJetDecayLengthVector.begin(); iter != CJetDecayLengthVector.end() ; iter++)
4053  {
4054  double decay_length = *iter;
4055  if (decay_length > c_decay_length0 ) c_decay_length0 = decay_length;
4056  }
4057 
4058  }
4059  }
4060  }
4061  }
4062 
4063  bjetdecaylength = b_decay_length0;
4064  cjetdecaylength = c_decay_length0;
4065 
4066  return;
4067 }
4068 
4069 
4070 int LCFIAIDAPlotProcessor::GetPDGFlavour(int code){
4071  //this little addition will take care of wierder mesons
4072 
4073  code = abs(code);
4074 
4075 
4076  //only really want the weak decaying particles
4077 
4078  if (code==411) return C_JET;
4079  if (code==421) return C_JET;
4080  if (code==431) return C_JET;
4081 
4082  if (code==511) return B_JET;
4083  if (code==521) return B_JET;
4084  if (code==531) return B_JET;
4085  if (code==541) return B_JET;
4086 
4087  if (code==4122) return C_JET;
4088  if (code==4222) return C_JET;
4089  if (code==4112) return C_JET;
4090  if (code==4212) return C_JET;
4091  if (code==4232) return C_JET;
4092  if (code==4132) return C_JET;
4093 
4094 
4095  if (code==5122) return B_JET;
4096  if (code==5222) return B_JET;
4097  if (code==5112) return B_JET;
4098  if (code==5212) return B_JET;
4099  if (code==5132) return B_JET;
4100  if (code==5232) return B_JET;
4101  if (code==5332) return B_JET;
4102 
4103  return 1;
4104 
4105 
4106  if( code > 10000 )
4107  {
4108  code = code%1000;
4109  }
4110 
4111  if( code > 1000 )
4112  {
4113 
4114  if( code/ 1000 == C_JET )
4115  {
4116  return C_JET;
4117  }
4118 
4119  if( code/ 1000 == B_JET )
4120  {
4121  return B_JET;
4122  }
4123  }
4124 
4125  if( code > 100 )
4126  {
4127  if( code / 100 == C_JET )
4128  {
4129  return C_JET;
4130 
4131  }
4132 
4133  if( code/100 == B_JET )
4134  {
4135  return B_JET;
4136  }
4137  }
4138 
4139  return 1;
4140 }
4141 
4142 
4143 
4144 #endif // endif of the check to see if MARLIN_USE_AIDA has been defined
LCFIAIDAPlotProcessor::FindCQVtx
int FindCQVtx(LCEvent *pEvent, unsigned int jetNumber)
Finds the vertex charge of the jet - using cuts tuned to find vertex charge for C-jets (from CVertexC...
LCFIAIDAPlotProcessor::CreateIntegralPlot
AIDA::IDataPointSet * CreateIntegralPlot(const AIDA::IHistogram1D *pNN, AIDA::IDataPointSet *pIntegral)
Makes a DataPointSet integrating a histogram from the first bin to the last bin – NOT USED...
LCFIAIDAPlotProcessor::FindTrueJetDecayLength
void FindTrueJetDecayLength(LCEvent *pEvent, unsigned int jetNumber, std::vector< double > &decaylengthvector, std::vector< double > &bjetdecaylengthvector, std::vector< double > &cjetdecaylengthvector)
Finds the true decay length of the longest b- or c- hadron in a jet.
LCFIAIDAPlotProcessor::CreateXYPlot
AIDA::IDataPointSet * CreateXYPlot(const AIDA::IDataPointSet *pDataPointSet0, const AIDA::IDataPointSet *pDataPointSet1, AIDA::IDataPointSet *xyPointSet, const int dim0=0, const int dim1=0)
Plots two DataPointSets against each other.
LCFIAIDAPlotProcessor::FindTrueJetFlavour
int FindTrueJetFlavour(LCEvent *pEvent, unsigned int jetNumber)
Finds the true flavour of the jet (uses TrueJetFlavourCollection)
LCFIAIDAPlotProcessor::CreateEfficiencyPlot
AIDA::IDataPointSet * CreateEfficiencyPlot(const AIDA::IHistogram1D *pSignal, AIDA::IDataPointSet *pDataPointSet)
Makes a DataPointSet of the tag efficiency e.g number of B-jets passing a given B-tag NN cut...
LCFIAIDAPlotProcessor::FindNumVertex
int FindNumVertex(LCEvent *pEvent, unsigned int jetNumber, unsigned int iInputsCollection)
Finds the number of vertices in an event (from the flavour tag inputs)
LCFIAIDAPlotProcessor::CreatePurityPlot
AIDA::IDataPointSet * CreatePurityPlot(const AIDA::IHistogram1D *pSignal, const AIDA::IHistogram1D *pBackground, AIDA::IDataPointSet *pDataPointSet)
Makes a DataPointSet of the tag purity e.g. N(B-jets passing NN cut)/N(all-jets passing NN cut) for a...
LCFIAIDAPlotProcessor::FindBQVtx
int FindBQVtx(LCEvent *pEvent, unsigned int jetNumber)
Finds the vertex charge of the jet - using cuts tuned to find vertex charge for B-jets (from BVertexC...
LCFIAIDAPlotProcessor::FindTrueJetPartonCharge
float FindTrueJetPartonCharge(LCEvent *pEvent, unsigned int jetNumber)
Finds the true charge of the parton producing a jet (uses TrueJetFlavourCollection) ...
LCFIAIDAPlotProcessor::CreateEfficiencyPlot2
AIDA::IDataPointSet * CreateEfficiencyPlot2(const AIDA::IHistogram1D *pAllEvents, const AIDA::IHistogram1D *pPassEvents, AIDA::IDataPointSet *pDataPointSet)
Makes a DataPointSet of histogram 1 divide by histogram 2 - this is an IDataPointSet as a histrogram ...
LCFIAIDAPlotProcessor::CreateLeakageRatePlot
AIDA::IDataPointSet * CreateLeakageRatePlot(const AIDA::IHistogram1D *pBackground, AIDA::IDataPointSet *pDataPointSet)
Makes a DataPointSet showing the tagging leakage e.g. the number of non-B-jets passing a given B-tag ...
LCFIAIDAPlotProcessor::CreateIntegralHistogram
AIDA::IHistogram1D * CreateIntegralHistogram(const AIDA::IHistogram1D *pNN, AIDA::IHistogram1D *pIntegral)
Makes a histogram integrating a histogram from the first bin to the last bin - THE ERRORS RETURNED AR...
LCFIAIDAPlotProcessor::FindTrueJetPDGCode
int FindTrueJetPDGCode(LCEvent *pEvent, unsigned int jetNumber)
Finds the PDG code of the hadron producing a jet (uses TrueJetFlavourCollection)
LCFIAIDAPlotProcessor::FindTrueJetType
int FindTrueJetType(LCEvent *pEvent, unsigned int jetNumber)
Finds the true flavour of a jet (uses TrueJetFlavourCollection)
LCFIAIDAPlotProcessor
LCFIAIDAPlotProcessor Class - make plots of the LCFI flavour tag and vertex charge code...
Definition: LCFIAIDAPlotProcessor.h:91
LCFIAIDAPlotProcessor::CreateVertexChargeLeakagePlot
void CreateVertexChargeLeakagePlot(AIDA::IDataPointSet *pBJetVtxChargeDPS, AIDA::IDataPointSet *pCJetVtxChargeDPS)
Makes DataPointSets for the number of.
LCFIAIDAPlotProcessor::FindTrueJetHadronCharge
float FindTrueJetHadronCharge(LCEvent *pEvent, unsigned int jetNumber)
Finds the true charge of the hadron producing a jet (uses TrueJetFlavourCollection) ...
Generated on Fri Dec 2 2016 12:36:49 for LCFIVertex by   doxygen 1.8.6