Newer
Older
TB_Chris / TbAnalysis / src / .svn / text-base / TbEfficiency.cpp.svn-base
  1. // Gaudi
  2. #include "GaudiKernel/PhysicalConstants.h"
  3. // TODO (hschindl)
  4. #include "GaudiKernel/IEventProcessor.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.  
  14. // Local
  15. #include "TbEfficiency.h"
  16.  
  17. DECLARE_ALGORITHM_FACTORY(TbEfficiency)
  18.  
  19. //=============================================================================
  20. // Standard constructor
  21. //=============================================================================
  22. TbEfficiency::TbEfficiency(const std::string& name,
  23. ISvcLocator* pSvcLocator)
  24. : TbAlgorithm(name, pSvcLocator),
  25. m_nTracksConsidered(0),
  26. m_nTracksAssociated(0),
  27. m_event(0),
  28. m_pitch(0.055) {
  29.  
  30. declareProperty("TrackLocation",
  31. m_trackLocation = LHCb::TbTrackLocation::Default);
  32. declareProperty("ClusterLocation",
  33. m_clusterLocation = LHCb::TbClusterLocation::Default);
  34. declareProperty("DUT", m_dut = 4);
  35. declareProperty("CheckHitDUT", m_checkHitDUT = true);
  36. declareProperty("CheckHitAlivePixel", m_checkHitAlivePixel = true);
  37. declareProperty("PointingResAllowance", m_pointingResAllowance = 0.01);
  38. declareProperty("PointingResAllowanceDeadPixel", m_pointingResAllowance_deadPixels = 1);
  39. declareProperty("TakeDeadPixelsFromFile", m_takeDeadPixelsFromFile = true);
  40. declareProperty("nTotalTracks", m_nTotalTracks = 0);
  41. declareProperty("MaxChi", m_maxChi = 0);
  42. }
  43.  
  44. //=============================================================================
  45. // Finalization
  46. //=============================================================================
  47. StatusCode TbEfficiency::finalize() {
  48.  
  49. info() <<"DUT Efficiency: "<<m_nTracksAssociated/(1.*m_nTracksConsidered)<<endmsg;
  50. info()<<"nTracksAssociated: "<<m_nTracksAssociated<<endmsg;
  51. info()<<"nTracksConsidered: "<<m_nTracksConsidered<<endmsg;
  52. info()<<"cf: "<<m_eff->GetEfficiency(1)<<"\t ±"<<m_eff->GetEfficiencyErrorLow(1)<<endmsg;
  53.  
  54. TFile * f = new TFile("EfficiencyResults.root", "RECREATE");
  55. h_row->Write();
  56. h_col->Write();
  57. m_eff->Write();
  58. h_hitmap->Write();
  59. h_pixel->Write();
  60. h_pixel2x2->Write();
  61.  
  62. m_eff->CreateGraph()->Write();
  63. h_hitmap->CreateHistogram()->Write();
  64. h_col->CreateGraph()->Write();
  65. h_row->CreateGraph()->Write();
  66. h_pixel->CreateHistogram()->Write();
  67. h_pixel2x2->CreateHistogram()->Write();
  68. f->Close();
  69.  
  70. return TbAlgorithm::finalize();
  71. }
  72.  
  73.  
  74. //=============================================================================
  75. // Initialization
  76. //=============================================================================
  77. StatusCode TbEfficiency::initialize() {
  78. // Initialise the base class.
  79. StatusCode sc = TbAlgorithm::initialize();
  80. if (sc.isFailure()) return sc;
  81.  
  82. uint nBins = 80;
  83. h_row = new TEfficiency("row", "row", 256, 0, 256);
  84. h_col = new TEfficiency("col", "col", 3*256+2, 0, 3*256+2);
  85. m_eff = new TEfficiency("eff", "eff", 1, 0, 1);
  86. h_hitmap = new TEfficiency("hitmap", "hitmap", 3*256+2, 0, 3*256+2, 256, 0, 256);
  87. h_pixel = new TEfficiency("pixel", "pixel", nBins, 0, m_pitch, nBins, 0, m_pitch);
  88. h_pixel2x2 = new TEfficiency("pixel2x2", "pixel2x2", 2*nBins, 0, 2*m_pitch, 2*nBins, 0, 2*m_pitch);
  89.  
  90. if (m_takeDeadPixelsFromFile) {
  91. TFile * fIn = new TFile("telPlaneRef.root", "READ");
  92. const std::string name = "Tb/TbHitMonitor/HitMap/Plane" + std::to_string(m_dut);
  93. m_deadPixelMap = (TH2F*)fIn->Get(name.c_str());
  94. info() <<"DeadPixelMap name: "<<m_deadPixelMap->GetName()<<endmsg;
  95. }
  96.  
  97. return StatusCode::SUCCESS;
  98. }
  99.  
  100. //=============================================================================
  101. // Main execution
  102. //=============================================================================
  103. StatusCode TbEfficiency::execute() {
  104. // Grab the tracks.
  105. const LHCb::TbTracks* tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
  106. if (!tracks) {
  107. return Error("No tracks in " + m_trackLocation);
  108. }
  109.  
  110. // Loop over the tracks.
  111. for (LHCb::TbTrack* track : *tracks) {
  112. if (m_nTotalTracks != 0 && m_nTracksConsidered >= m_nTotalTracks) {
  113. SmartIF<IEventProcessor> app(serviceLocator()->service("ApplicationMgr"));
  114. if (app) app->stopRun();
  115. return StatusCode::SUCCESS;
  116. }
  117.  
  118. // Find the local position of the track intercept with the DUT.
  119. const Gaudi::XYZPoint interceptUG = geomSvc()->intercept(track, m_dut);
  120. const auto interceptUL = geomSvc()->globalToLocal(interceptUG, m_dut);
  121.  
  122. if (!passSelectionCriteria(track, interceptUL)) continue;
  123. m_nTracksConsidered++;
  124. bool pass = false;
  125. if (track->associatedClusters().size() > 0) pass = true;
  126.  
  127. int row = interceptUL.y()/m_pitch - 0.5;
  128. int col = interceptUL.x()/m_pitch - 0.5;
  129.  
  130. h_pixel->Fill(pass, fmod(interceptUL.x(), m_pitch), fmod(interceptUL.y(), m_pitch));
  131. h_pixel2x2->Fill(pass, fmod(interceptUL.x(), 2*m_pitch), fmod(interceptUL.y(), 2*m_pitch));
  132. h_row->Fill(pass, row);
  133. h_col->Fill(pass, col);
  134. h_hitmap->Fill(pass, col, row);
  135. m_eff->Fill(pass, 0);
  136. if (pass) m_nTracksAssociated++;
  137. if (!pass && m_event == 1050) info()<<track->htime()<<endmsg;
  138. }
  139.  
  140. m_event++;
  141. return StatusCode::SUCCESS;
  142. }
  143.  
  144.  
  145. //=============================================================================
  146.  
  147. bool TbEfficiency::passSelectionCriteria(LHCb::TbTrack * track, Gaudi::XYZPoint interceptUL)
  148. {
  149. if (m_maxChi > 0 && track->chi2PerNdof() > m_maxChi) return false;
  150. if (m_checkHitDUT && !passedThroughDUT(interceptUL)) {
  151. return false;
  152. }
  153. if (m_checkHitAlivePixel && !passedThroughAlivePixel(interceptUL)) return false;
  154.  
  155. return true;
  156. }
  157.  
  158.  
  159. //=============================================================================
  160.  
  161. bool TbEfficiency::passedThroughDUT(Gaudi::XYZPoint interceptUL)
  162. {
  163. double xLow = 0.0 + m_pointingResAllowance;
  164. double yLow = 0.0 + m_pointingResAllowance;
  165. double yUp = 14.08 - m_pointingResAllowance;
  166. double xUp = 14.08;
  167.  
  168. if (geomSvc()->modules().at(m_dut)->nChips() == 1)
  169. xUp = 14.08 - m_pointingResAllowance;
  170. else if (geomSvc()->modules().at(m_dut)->nChips() == 3)
  171. xUp = 42.35 - m_pointingResAllowance;
  172. else error()<<"Unknown device!"<<endmsg;
  173.  
  174. if (interceptUL.x() < xLow ||
  175. interceptUL.x() > xUp ||
  176. interceptUL.y() < yLow ||
  177. interceptUL.y() > yUp) {
  178. counter("nTracksNotStrikingDUT") += 1;
  179. return false;
  180. }
  181. return true;
  182. }
  183.  
  184.  
  185. //=============================================================================
  186.  
  187. bool TbEfficiency::passedThroughAlivePixel(Gaudi::XYZPoint interceptUL)
  188. {
  189. int row = interceptUL.y()/0.055 - 0.5;
  190. int col = interceptUL.x()/0.055 - 0.5;
  191. int one = 1;
  192. for (int iRow = row - m_pointingResAllowance_deadPixels;
  193. iRow != row+m_pointingResAllowance_deadPixels+one; iRow++) {
  194. for (int iCol = col - m_pointingResAllowance_deadPixels;
  195. iCol != col + m_pointingResAllowance_deadPixels+one; iCol++) {
  196.  
  197. int colLim;
  198. if (geomSvc()->modules().at(m_dut)->nChips() == 3) colLim = 3*256+2;
  199. else colLim = 256;
  200. if (iRow < 0 || iCol < 0 || iRow >= 256 || iCol >= colLim) continue;
  201. if (pixelSvc()->isMasked(pixelSvc()->address(iCol, iRow), m_dut)) {
  202. counter("nTracksStrikingMaskedPixels") += 1;
  203. return false;
  204. }
  205. else if (m_takeDeadPixelsFromFile && m_deadPixelMap->GetBinContent(iCol, iRow) == 0) {
  206. counter("nTracksStrikingDeadPixels") += 1;
  207. return false;
  208. }
  209.  
  210. }
  211. }
  212. return true;
  213. }
  214.  
  215.  
  216. //=============================================================================