Newer
Older
TB_Chris / TbAnalysis / src / .svn / text-base / TbVertexing.cpp.svn-base
  1. // Gaudi
  2. #include "GaudiKernel/PhysicalConstants.h"
  3. #include "GaudiUtils/HistoLabels.h"
  4. #include "GaudiUtils/Aida2ROOT.h"
  5.  
  6. // Tb/TbEvent
  7. #include "Event/TbTrack.h"
  8. #include "Event/TbCluster.h"
  9.  
  10. // Tb/TbKernel
  11. #include "TbKernel/TbConstants.h"
  12. #include "TbKernel/TbModule.h"
  13. #include "TbKernel/TbFunctors.h"
  14.  
  15. // Local
  16. #include "TbVertexing.h"
  17.  
  18. // ROOT
  19. #include "TH1D.h"
  20.  
  21. using namespace Gaudi::Utils::Histos;
  22.  
  23. DECLARE_ALGORITHM_FACTORY(TbVertexing)
  24.  
  25. template <class TYPE>
  26. class itlowerBound {
  27. public:
  28. bool operator()(TYPE lhs, const double t) const { return (*lhs)->htime() < t; }
  29. };
  30.  
  31. //=============================================================================
  32. // Standard constructor
  33. //=============================================================================
  34. TbVertexing::TbVertexing(const std::string& name,
  35. ISvcLocator* pSvcLocator)
  36. : TbAlgorithm(name, pSvcLocator) {
  37. declareProperty("TimeWindow",m_timeWindow);
  38. declareProperty("MaxDoca2",m_maxDoca2);
  39. declareProperty("TrackLocation", m_trackLocation = LHCb::TbTrackLocation::Default);
  40. }
  41.  
  42. //=============================================================================
  43. // Initialization
  44. //=============================================================================
  45. StatusCode TbVertexing::initialize() {
  46.  
  47. // Initialise the base class.
  48. StatusCode sc = TbAlgorithm::initialize();
  49. if (sc.isFailure()) return sc;
  50. m_trackFit = tool<ITbTrackFit>("TbTrackFit", "Fitter", this);
  51. if( m_trackFit == NULL ){
  52. error() << "Failed to initialise trackfit" << endmsg;
  53. return StatusCode::FAILURE;
  54. }
  55. NTupleFilePtr file1(ntupleSvc(), "/NTUPLES/FILE1");
  56. NTuplePtr nt(ntupleSvc(), "/NTUPLES/FILE1/TbTupleWriter/Vertices");
  57. // Check if already booked.
  58. nt = ntupleSvc()->book("/NTUPLES/FILE1/TbTupleWriter/Vertices",
  59. CLID_ColumnWiseTuple, "nTuple of Vertices");
  60. nt->addItem("nPlanes", m_index, 0, 10);
  61. nt->addIndexedItem("X", m_index, m_vtxX);
  62. nt->addIndexedItem("Y", m_index, m_vtxY);
  63. nt->addIndexedItem("Z", m_index, m_vtxZ);
  64. nt->addItem("Matched", m_matched );
  65. // nt->addItem("Tk1Size",m_tkSize1);
  66. // nt->addItem("Tk2Size",m_tkSize2);
  67. return StatusCode::SUCCESS;
  68. }
  69.  
  70. StatusCode TbVertexing::finalize() {
  71.  
  72. return StatusCode::SUCCESS;
  73. }
  74.  
  75. //=============================================================================
  76. // Main execution
  77. //=============================================================================
  78. StatusCode TbVertexing::execute() {
  79.  
  80. // Grab the tracks.
  81.  
  82. const LHCb::TbTracks* tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
  83. if (tracks) {
  84. for( auto& track1 : *tracks ){
  85.  
  86. const double tMin = track1->htime() - m_timeWindow;
  87. const double tMax = track1->htime() + m_timeWindow;
  88. auto track2 = std::lower_bound(tracks->begin(),
  89. tracks->end(),
  90. tMin, lowerBound<const LHCb::TbTrack*>());
  91.  
  92. for(; track2 != tracks->end() && (*track2)->htime() < tMax ; ++track2 ){
  93. if( track1 == *track2 ) continue;
  94. const double DOCA = doca2(track1,*track2);
  95. if( DOCA > m_maxDoca2 ) continue;
  96. /// track separation ///
  97. const double avgX = ( track1->firstState().tx() + (*track2)->firstState().tx() ) /2.;
  98. const double avgY = ( track1->firstState().ty() + (*track2)->firstState().ty() ) /2.;
  99.  
  100. const double dt = pow ( ( track1->firstState().tx() - (*track2)->firstState().tx() ) / avgX, 2.0 ) +
  101. pow ( ( ( track1->firstState().ty() - (*track2)->firstState().ty() ) ) / avgY , 2.0 ) ;
  102. auto chi2_pos = vertexPosition( {track1,*track2} );
  103. plot( sqrt(DOCA), "DOCA", 0, 0.5, 100 ); /// 500 micron
  104. plot( dt, "Angle", -2., 2., 1000);
  105. plot( chi2_pos[0] , "VertexChi2_X",0.,14.,100);
  106. plot( chi2_pos[1] , "VertexChi2_Y",0.,14.,100);
  107. plot( chi2_pos[2] , "VertexChi2_Z",-700.,650.,2500);
  108. m_vtxX[0] = chi2_pos[0];
  109. m_vtxY[0] = chi2_pos[1];
  110. m_vtxZ[0] = chi2_pos[2];
  111.  
  112. ntupleSvc()->writeRecord("/NTUPLES/FILE1/TbTupleWriter/Vertices");
  113. m_index = 10;
  114. m_matched = matched( track1, *track2 );
  115. for( unsigned int plane = 0 ; plane != 9 ; ++plane ){
  116. m_trackFit->maskPlane( plane );
  117. m_trackFit->fit( track1 );
  118. m_trackFit->fit( *track2 );
  119. //m_vtxX[plane+1]
  120. auto up = vertexPosition( {track1,*track2});
  121. m_vtxX[plane+1] = up[0];
  122. m_vtxY[plane+1] = up[1];
  123. m_vtxZ[plane+1] = up[2];
  124. plot( up[2] , "VertexChi2_Z_"+std::to_string(plane),-700.,650.,2500);
  125. m_trackFit->unmaskPlane(plane);
  126. }
  127. m_trackFit->fit( track1 );
  128. m_trackFit->fit( *track2);
  129. }
  130. }
  131.  
  132. }
  133.  
  134. return StatusCode::SUCCESS;
  135. }
  136.