Newer
Older
TB_Chris / TbAnalysis / src / .svn / text-base / TbResolutionMonitor.cpp.svn-base
  1. // Gaudi
  2. #include "GaudiKernel/PhysicalConstants.h"
  3.  
  4. // Tb/TbEvent
  5. #include "Event/TbTrack.h"
  6. #include "Event/TbCluster.h"
  7.  
  8. // Tb/TbKernel
  9. #include "TbKernel/TbConstants.h"
  10.  
  11. // Local
  12. #include "TbResolutionMonitor.h"
  13.  
  14. DECLARE_ALGORITHM_FACTORY(TbResolutionMonitor)
  15.  
  16. //=============================================================================
  17. // Standard constructor
  18. //=============================================================================
  19. TbResolutionMonitor::TbResolutionMonitor(const std::string& name,
  20. ISvcLocator* pSvcLocator)
  21. : TbAlgorithm(name, pSvcLocator) {
  22.  
  23. declareProperty("TrackLocation",
  24. m_trackLocation = LHCb::TbTrackLocation::Default);
  25. declareProperty("ClusterLocation",
  26. m_clusterLocation = LHCb::TbClusterLocation::Default);
  27.  
  28. declareProperty("DUTs", m_duts);
  29. }
  30.  
  31. //=============================================================================
  32. // Initialization
  33. //=============================================================================
  34. StatusCode TbResolutionMonitor::initialize() {
  35.  
  36. // Initialise the base class.
  37. StatusCode sc = TbAlgorithm::initialize();
  38. if (sc.isFailure()) return sc;
  39. return StatusCode::SUCCESS;
  40. }
  41.  
  42. //=============================================================================
  43. // Main execution
  44. //=============================================================================
  45. StatusCode TbResolutionMonitor::execute() {
  46.  
  47. // Grab the tracks.
  48. const LHCb::TbTracks* tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
  49. if (!tracks) {
  50. error() << "No tracks in " << m_trackLocation << endmsg;
  51. return StatusCode::FAILURE;
  52. }
  53.  
  54. for (const auto dut : m_duts) {
  55. // Get the clusters for this plane.
  56. const std::string clusterLocation = m_clusterLocation + std::to_string(dut);
  57. const LHCb::TbClusters* clusters = getIfExists<LHCb::TbClusters>(clusterLocation);
  58. for (const LHCb::TbCluster* cluster : *clusters) {
  59. plot(cluster->charge(), "ChargeAll", "ChargeAll", 0., 50000., 200);
  60. }
  61. if (!clusters) continue;
  62. // Loop over the tracks.
  63. for (const LHCb::TbTrack* track : *tracks) {
  64. // Cut on track quality.
  65. if (track->chi2PerNdof() > 20.) continue;
  66. const Gaudi::XYZPoint pGlobal = geomSvc()->intercept(track, dut);
  67. const auto pLocal = geomSvc()->globalToLocal(pGlobal, dut);
  68. if (pLocal.x() < 0. || pLocal.y() < 0.) continue;
  69. const double fcol = pLocal.x() / Tb::PixelPitch;
  70. const double frow = pLocal.y() / Tb::PixelPitch;
  71. const int col = int(fcol);
  72. const int row = int(frow);
  73. if (col > 255 || row > 255) continue;
  74. // Calculate inter-pixel coordinates.
  75. const double xCell = 1.e3 * (fcol - col) * Tb::PixelPitch;
  76. const double yCell = 1.e3 * (frow - row) * Tb::PixelPitch;
  77. const LHCb::TbCluster* match = NULL;
  78. for (const LHCb::TbCluster* cluster : *clusters) {
  79. const double dt = cluster->htime() - track->htime();
  80. if (fabs(dt) > 200.) continue;
  81. const double dx = cluster->xloc() - pLocal.x();
  82. const double dy = cluster->yloc() - pLocal.y();
  83. if (fabs(dx) < 0.15 && fabs(dy) < 0.15) {
  84. match = cluster;
  85. break;
  86. }
  87. }
  88. plot2D(xCell, yCell, "NTracks", "NTracks",
  89. 0., 55., 0., 55., 9, 9);
  90. if (!match) continue;
  91. plot2D(xCell, yCell, "NClusters", "NClusters",
  92. 0., 55., 0., 55., 9, 9);
  93. plot(match->x() - pGlobal.x(), "ResGlobalX", "ResGlobalX", -0.2, 0.2, 100);
  94. plot(match->y() - pGlobal.y(), "ResGlobalY", "ResGlobalY", -0.2, 0.2, 100);
  95. plot(match->xloc() - pLocal.x(), "ResLocalX", "ResLocalX", -0.2, 0.2, 100);
  96. plot(match->yloc() - pLocal.y(), "ResLocalY", "ResLocalX", -0.2, 0.2, 100);
  97. plot(match->htime() - track->htime(), "ResTime", "ResTime", -200., 200., 400);
  98. const unsigned int cls = match->size();
  99. if (cls > 4) continue;
  100. const std::string title = "InterceptCls" + std::to_string(cls);
  101. plot2D(xCell, yCell, title, title,
  102. 0., 55., 0., 55., 50, 50);
  103. }
  104. }
  105. return StatusCode::SUCCESS;
  106. }
  107.