Newer
Older
TB_Chris / TbAlgorithms / src / .svn / text-base / TbCalibration.cpp.svn-base
  1. #include <fstream>
  2.  
  3. // Tb/TbKernel
  4. #include "TbKernel/TbModule.h"
  5.  
  6. // Local
  7. #include "TbCalibration.h"
  8.  
  9. // ROOT
  10. #include "TFitResult.h"
  11. #include "TFitResultPtr.h"
  12. #include "TH1D.h"
  13.  
  14. /// GAUDI
  15. #include "GaudiUtils/Aida2ROOT.h"
  16.  
  17. DECLARE_ALGORITHM_FACTORY(TbCalibration)
  18.  
  19. //=============================================================================
  20. // Standard constructor
  21. //=============================================================================
  22. TbCalibration::TbCalibration(const std::string& name, ISvcLocator* pSvcLocator)
  23. : TbAlgorithm(name, pSvcLocator) {
  24. declareProperty("PixelConfigFile", m_pixelSvcConfig = "PixelConfig.dat");
  25. declareProperty("TimingConfigFile", m_timingSvcConfig = "TimingConfig.dat");
  26. declareProperty("CheckHotPixels", m_checkHotPixels = false);
  27. declareProperty("CheckSynchronisation", m_checkSyncronisation = false);
  28. declareProperty("SyncMethod", m_syncMethod = 1);
  29. declareProperty("CheckColumnOffsets", m_checkColumnOffsets = false);
  30. declareProperty("HitLocation", m_hitLocation = LHCb::TbHitLocation::Default);
  31. declareProperty("DuT", m_dut = 9999);
  32. declareProperty("ClusterLocation",
  33. m_clusterLocation = LHCb::TbClusterLocation::Default);
  34. }
  35.  
  36. //=============================================================================
  37. // Initialisation
  38. //=============================================================================
  39. StatusCode TbCalibration::initialize() {
  40.  
  41. // Initialise the base class.
  42. StatusCode sc = TbAlgorithm::initialize();
  43. if (sc.isFailure()) return sc;
  44.  
  45. if (m_checkColumnOffsets) m_offsets.resize(m_nPlanes, PROFILE1D(128));
  46. if (m_checkHotPixels) m_hitMaps.resize(m_nDevices, PROFILE2D(256, 256));
  47. if (m_checkSyncronisation) m_sync.resize(m_nPlanes);
  48. return sc;
  49. }
  50.  
  51. //=============================================================================
  52. // Main execution
  53. //=============================================================================
  54. StatusCode TbCalibration::execute() {
  55.  
  56. if (m_checkColumnOffsets) {
  57. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  58. LHCb::TbClusters* clusters =
  59. getIfExists<LHCb::TbClusters>(m_clusterLocation + std::to_string(i));
  60. if (!clusters) continue;
  61. columnOffset_execute(clusters);
  62. }
  63. }
  64. if (m_checkHotPixels) {
  65. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  66. LHCb::TbHits* hits =
  67. getIfExists<LHCb::TbHits>(m_hitLocation + std::to_string(i));
  68. if (!hits) continue;
  69. hotPixel_execute(hits);
  70. }
  71. }
  72. if (m_checkSyncronisation) {
  73. if (m_syncMethod == 0) {
  74. std::vector<LHCb::TbClusters*> clusters;
  75. for (unsigned int i = 0; i < m_nPlanes; ++i)
  76. clusters.push_back(getIfExists<LHCb::TbClusters>(m_clusterLocation +
  77. std::to_string(i)));
  78. sync_execute(clusters);
  79. }
  80. }
  81. return StatusCode::SUCCESS;
  82. }
  83.  
  84. //=============================================================================
  85. // Finalisation
  86. //=============================================================================
  87. StatusCode TbCalibration::finalize() {
  88.  
  89. std::ofstream pixelFile(m_pixelSvcConfig.c_str());
  90. std::ofstream timingFile(m_timingSvcConfig.c_str());
  91.  
  92. if (m_checkHotPixels) {
  93. pixelFile << "Mask" << std::endl;
  94. for (unsigned int i = 0; i < m_nDevices; ++i)
  95. hotPixelAnalysis(m_hitMaps[i], geomSvc()->module(i)->id(), pixelFile);
  96. }
  97. if (m_checkSyncronisation) {
  98. if (m_syncMethod == 0)
  99. for (unsigned int i = 0; i < m_nPlanes; ++i)
  100. syncAnalysis(m_sync[i].avg(), geomSvc()->module(i)->id(), timingFile);
  101. else if (m_syncMethod == 1)
  102. syncOffset2(timingFile);
  103. }
  104. if (m_checkColumnOffsets) {
  105. pixelFile << "Offset" << std::endl;
  106. for (unsigned int i = 0; i < m_nPlanes; ++i)
  107. columnOffsetAnalysis(m_offsets[i], geomSvc()->module(i)->id(), pixelFile);
  108. }
  109. pixelFile.close();
  110. timingFile.close();
  111. return StatusCode::SUCCESS;
  112. }
  113.  
  114. //=============================================================================
  115. // Identify hot pixels
  116. //=============================================================================
  117. void TbCalibration::hotPixelAnalysis(const PROFILE2D& hitMap,
  118. const std::string& plane,
  119. std::ostream& os) {
  120.  
  121. // Based on Hella's script for identifying hot pixels.
  122. // Reject pixels which are more than 40x average
  123. for (int col = 0; col < 256; col++) {
  124. for (int row = 0; row < 256; row++) {
  125. double count = (hitMap[col][row]).n();
  126. std::vector<double> nn = hitMap.neighbours(col, row);
  127.  
  128. // Cull outliers if not at edge
  129. if (nn.size() == 8) {
  130. std::sort(nn.begin(), nn.end());
  131. nn.erase(nn.begin(), nn.begin() + 2);
  132. nn.erase(nn.end() - 2, nn.end());
  133. }
  134.  
  135. double total = std::accumulate(nn.begin(), nn.end(), 0);
  136.  
  137. if (count > 10 * total) {
  138. os << plane << " " << std::setw(3) << col << " " << row << std::endl;
  139. }
  140. }
  141. }
  142. }
  143.  
  144. void TbCalibration::syncOffset2(std::ostream& os) {
  145.  
  146. os << "Timing" << std::endl;
  147.  
  148. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  149. std::string title =
  150. "Tb/TbTrackPlots/BiasedResiduals/Time/Plane" + std::to_string(i);
  151. if (i == m_dut)
  152. title = "Tb/TbDUTMonitor/ResidualsTime/Plane" + std::to_string(i);
  153. AIDA::IHistogram1D* hAida = NULL;
  154. StatusCode sc = GaudiAlgorithm::histoSvc()->retrieveObject(title, hAida);
  155. auto hRoot = Gaudi::Utils::Aida2ROOT::aida2root(hAida);
  156. os << geomSvc()->module(i)->id() << std::setw(3) << " "
  157. << -hRoot->GetMean() << std::endl;
  158. }
  159. }
  160.  
  161. //=============================================================================
  162. // Calculate time offsets for each column
  163. //=============================================================================
  164. void TbCalibration::columnOffsetAnalysis(const PROFILE1D& avg_difference,
  165. const std::string& plane,
  166. std::ostream& os) {
  167. const unsigned int nDcols = 128;
  168. std::vector<int> offsets(nDcols, 0);
  169. double sum(0.), total(0.);
  170. for (auto& d : avg_difference) {
  171. sum += d.val();
  172. total += d.n();
  173. }
  174. double avg = sum / total;
  175.  
  176. for (unsigned int dcol = 1; dcol < nDcols; ++dcol) {
  177. const double difference =
  178. avg_difference[dcol].avg() - avg - 25 * offsets[dcol - 1];
  179. if (difference > 10.)
  180. offsets[dcol] = -1;
  181. else if (difference < -10.)
  182. offsets[dcol] = +1;
  183. }
  184. // calculate the average offset
  185.  
  186. // deal with the special case where the first super column is desynchronised,
  187. // in which case, everything will be shifted by 25 ns in this definition.
  188. double avg_offset(0.);
  189. for (const auto& d : offsets) avg_offset += d;
  190. avg_offset /= 128.;
  191.  
  192. if (avg_offset > 0.5)
  193. for (auto& d : offsets) d--;
  194.  
  195. else if (avg_offset < -0.5)
  196. for (auto& d : offsets) d++;
  197.  
  198. for (unsigned int dcol = 0; dcol < nDcols; ++dcol) {
  199. if (offsets[dcol] != 0)
  200. os << plane << " " << std::setw(3) << dcol << " " << offsets[dcol]
  201. << std::endl;
  202. }
  203. }
  204.  
  205. //=============================================================================
  206. // Fill the data for calculating the column time offsets
  207. //=============================================================================
  208. void TbCalibration::columnOffset_execute(const LHCb::TbClusters* clusters) {
  209. if (!clusters || clusters->empty()) return;
  210. LHCb::TbClusters::const_iterator itc;
  211. for (itc = clusters->begin(); itc != clusters->end(); ++itc) {
  212. const unsigned int plane = (*itc)->plane();
  213. if ((*itc)->size() == 1) continue;
  214. auto hits = (*itc)->hits();
  215. for (auto& ih0 : hits) {
  216. for (auto& ih1 : hits) {
  217. const LHCb::TbHit* h0 = ih0;
  218. const LHCb::TbHit* h1 = ih1;
  219. const int col0 = h0->col();
  220. const int col1 = h1->col();
  221. if (abs(col0 - col1) != 1) continue;
  222. if (h0->row() == h1->row() && h0->col() / 2 != h1->col() / 2) {
  223. if (h0->col() > h1->col()) std::swap(h0, h1);
  224. m_offsets[plane][(int)(h1->col() / 2)].add(h1->htime() - h0->htime());
  225. }
  226. }
  227. }
  228. }
  229. }
  230.  
  231. //=============================================================================
  232. // Fill the data for checking the synchronisation.
  233. //=============================================================================
  234. void TbCalibration::sync_execute(
  235. const std::vector<LHCb::TbClusters*>& clusters) {
  236.  
  237. LHCb::TbClusters* plane0_clusters = clusters[0];
  238. for (auto& c : *plane0_clusters) {
  239. for (unsigned int i = 1; i < m_nPlanes; ++i) {
  240. double nearest = nearestHit(c, clusters[i]);
  241. if (nearest != 9999) m_sync[i].add(nearest);
  242. }
  243. }
  244. }
  245.  
  246. //=============================================================================
  247. // Get the smallest time difference of a list of clusters.
  248. //=============================================================================
  249. double TbCalibration::nearestHit(const LHCb::TbCluster* cluster,
  250. const LHCb::TbClusters* clusters) {
  251. if (!clusters || !cluster || clusters->empty()) return 9999;
  252. double minTime = cluster->htime() - 50.;
  253. double maxTime = cluster->htime() + 50.;
  254. LHCb::TbClusters::const_iterator c = std::lower_bound(
  255. clusters->begin(), clusters->end(), minTime, lowerBound());
  256. double nn = 9999;
  257. for (; c != clusters->end() && (*c)->htime() < maxTime; ++c) {
  258. if (std::abs((*c)->htime() - cluster->htime()) < std::abs(nn))
  259. nn = (*c)->htime() - cluster->htime();
  260. }
  261. return nn;
  262. }
  263.  
  264. //=============================================================================
  265. // Fill the hit map for identifying hot pixels
  266. //=============================================================================
  267. void TbCalibration::hotPixel_execute(const LHCb::TbHits* hits) {
  268. if (!hits || hits->empty()) return;
  269. const unsigned int device = (*hits->begin())->device();
  270. for (auto& h : *hits) m_hitMaps[device][h->col()][h->row()].add(1);
  271. }
  272.  
  273. //=============================================================================
  274. // Save synchronisation info to file.
  275. //=============================================================================
  276. void TbCalibration::syncAnalysis(const double& sync, const std::string& plane,
  277. std::ostream& os) {
  278.  
  279. os << plane << std::setw(3) << " " << sync << std::endl;
  280. }