Newer
Older
TB_Chris / TbAnalysis / src / .svn / text-base / TbVertexing.h.svn-base
  1. #pragma once
  2.  
  3. // AIDA
  4. #include "AIDA/IProfile1D.h"
  5. #include "AIDA/IHistogram1D.h"
  6. #include "AIDA/IHistogram2D.h"
  7. // Tb/TbKernel
  8. #include "TbKernel/TbAlgorithm.h"
  9. #include "TbKernel/ITbTrackFit.h"
  10. #include "Event/TbTrack.h"
  11.  
  12. #include "TVector3.h"
  13. #include "TMatrixD.h"
  14. #include "TVectorD.h"
  15.  
  16. class TbVertexing : public TbAlgorithm {
  17. public:
  18. /// Constructor
  19. TbVertexing(const std::string& name, ISvcLocator* pSvcLocator);
  20. /// Destructor
  21. virtual ~TbVertexing() {}
  22. virtual StatusCode finalize();
  23. virtual StatusCode initialize(); ///< Algorithm initialization
  24. virtual StatusCode execute(); ///< Algorithm execution
  25.  
  26. private:
  27. double m_timeWindow;
  28. double m_maxDoca2;
  29. NTuple::Item<unsigned int> m_index;
  30. NTuple::Item<unsigned int> m_matched;
  31. NTuple::Item<unsigned int> m_common;
  32. NTuple::Array<double> m_vtxX;
  33. NTuple::Array<double> m_vtxY;
  34. NTuple::Array<double> m_vtxZ;
  35. ITbTrackFit* m_trackFit;
  36.  
  37. std::string m_trackLocation;
  38.  
  39. double doca2( const LHCb::TbTrack* track1, const LHCb::TbTrack* track2 ) const {
  40. const TVector3 r1( track1->firstState().x(), track1->firstState().y(), 0. );
  41. const TVector3 r2( track2->firstState().y(), track2->firstState().y(), 0. );
  42. const TVector3 t1( track1->firstState().tx(), track1->firstState().ty(), 1.);
  43. const TVector3 t2( track2->firstState().tx(), track2->firstState().ty(), 1.);
  44. const double z0 = - (t1 -t2).Dot( r1 - r2 ) / (t1-t2).Mag2();
  45. return (r1-r2).Mag2() + 2*z0*(t1-t2).Dot(r1-r2) + z0*z0*(t1-t2).Mag2();
  46. };
  47.  
  48. /// checks if has > 3 vertices in the forward direction, with unmatched clusters ///
  49. int matched( const LHCb::TbTrack* track1, const LHCb::TbTrack* track2 ) const {
  50. auto clusters1 = track1->clusters();
  51. auto clusters2 = track2->clusters();
  52. int nCommon = 0;
  53. int nMatched = 0;
  54. for( auto& cluster1 : clusters1 ){
  55. for( auto& cluster2 : clusters2 ){
  56. if( cluster1->plane() == cluster2->plane() ){
  57. if( cluster1 != cluster2 ) nMatched++;
  58. else nCommon++;
  59. break;
  60. }
  61. }
  62. }
  63. return nMatched;
  64. }
  65.  
  66. TVectorD vertexPosition( const std::vector<const LHCb::TbTrack*>& tracks ) const {
  67. TMatrixD cov(3,3);
  68. TVectorD res(3);
  69. for( auto& track : tracks ){
  70. TVectorD t(3);
  71. t[0] = track->firstState().tx();
  72. t[1] = track->firstState().ty();
  73. t[2] = 1;
  74. const double norm = t.Norm2Sqr();
  75. TVectorD r(3);
  76. r[0] = track->firstState().x();
  77. r[1] = track->firstState().y();
  78. r[2] = 0;
  79. double ip = t[0]*r[0] + t[1]*r[1] + t[2]*r[2];
  80.  
  81. res += r - (ip /norm ) * t;
  82.  
  83. for( unsigned int i=0;i!=3;++i){
  84. for( unsigned int j=0; j !=3 ;++j){
  85. if(i==j) cov[i][j] += 1 - t[i]*t[j] / norm;
  86. else cov[i][j] += - t[i]*t[j] / norm;
  87. }
  88. }
  89. }
  90. TMatrixD inv = cov.Invert();
  91. return inv * res ;
  92. }
  93.  
  94. };