Newer
Older
TB_Chris / TbAlgorithms / src / .svn / text-base / TbClustering.cpp.svn-base
  1. // Gaudi
  2. #include "GaudiKernel/PhysicalConstants.h"
  3.  
  4. // Tb/TbKernel
  5. #include "TbKernel/TbFunctors.h"
  6. #include "TbKernel/TbConstants.h"
  7. #include "TbKernel/TbModule.h"
  8. #include "TbKernel/TbCondFile.h"
  9.  
  10. // Local
  11. #include "TbClustering.h"
  12.  
  13. DECLARE_ALGORITHM_FACTORY(TbClustering)
  14.  
  15. //=============================================================================
  16. // Standard constructor
  17. //=============================================================================
  18. TbClustering::TbClustering(const std::string& name, ISvcLocator* pSvcLocator)
  19. : TbAlgorithm(name, pSvcLocator) {
  20.  
  21. declareProperty("HitLocation", m_hitLocation = LHCb::TbHitLocation::Default);
  22. declareProperty("ClusterLocation",
  23. m_clusterLocation = LHCb::TbClusterLocation::Default);
  24. declareProperty("TimeWindow", m_twindow = 100. * Gaudi::Units::ns);
  25. declareProperty("SearchDist", m_searchDist = 1);
  26. declareProperty("ClusterErrorMethod", m_clusterErrorMethod = 0);
  27. }
  28.  
  29. //=============================================================================
  30. // Initialization
  31. //=============================================================================
  32. StatusCode TbClustering::initialize() {
  33.  
  34. // Initialise the base class.
  35. StatusCode sc = TbAlgorithm::initialize();
  36. if (sc.isFailure()) return sc;
  37.  
  38. m_eta.resize(m_nPlanes);
  39. for (const auto& filename : dataSvc()->getEtaConfig()) {
  40. info() << "Importing eta corrections from " << filename << endmsg;
  41. readEta(filename);
  42. }
  43. return StatusCode::SUCCESS;
  44. }
  45.  
  46. //=============================================================================
  47. // Main execution
  48. //=============================================================================
  49. StatusCode TbClustering::execute() {
  50.  
  51. // Get the hits to be clustered - loop over planes.
  52. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  53. const std::string ext = std::to_string(i);
  54. const std::string hitLocation = m_hitLocation + ext;
  55. const LHCb::TbHits* hits = getIfExists<LHCb::TbHits>(hitLocation);
  56. if (!hits) return Error("No hits in " + hitLocation);
  57. // Create a cluster container and transfer its ownership to the TES.
  58. LHCb::TbClusters* clusters = new LHCb::TbClusters();
  59. put(clusters, m_clusterLocation + ext);
  60. // Skip masked planes.
  61. if (masked(i)) continue;
  62. // Keep track of which hits have been clustered.
  63. std::vector<bool> used(hits->size(), false);
  64. std::vector<const LHCb::TbHit*> pixels;
  65. pixels.reserve(100);
  66. // Cycle over the (time ordered) hits.
  67. for (LHCb::TbHits::const_iterator it = hits->begin(); it != hits->end();
  68. ++it) {
  69. // Skip hits which are already part of a cluster.
  70. if (used[(*it)->key()]) continue;
  71. // Start a cluster from this seed pixel and tag the seed as used.
  72. used[(*it)->key()] = true;
  73. pixels.clear();
  74. pixels.push_back(*it);
  75. // Get the search range.
  76. const double tMax = (*it)->htime() + m_twindow;
  77. LHCb::TbHits::const_iterator end = std::upper_bound(
  78. it, hits->end(), tMax, lowerBound<const LHCb::TbHit*>());
  79. // Add neighbouring hits in the cluster time window.
  80. addNeighbouringHits(pixels, it, end, used);
  81. // Sort the hits by time.
  82. std::sort(pixels.begin(), pixels.end(),
  83. TbFunctors::LessByTime<const LHCb::TbHit*>());
  84. // Finally, set the remaining cluster attributes according to its hits.
  85. completeCluster(i, pixels, clusters);
  86. }
  87. // Sort the clusters by time.
  88. std::sort(clusters->begin(), clusters->end(),
  89. TbFunctors::LessByTime<LHCb::TbCluster*>());
  90. // Fill the counter.
  91. counter("NumberOfClusters" + ext) += clusters->size();
  92. }
  93.  
  94. return StatusCode::SUCCESS;
  95. }
  96.  
  97. //=============================================================================
  98. // Finding touching hits
  99. //=============================================================================
  100. void TbClustering::addNeighbouringHits(std::vector<const LHCb::TbHit*>& pixels,
  101. LHCb::TbHits::const_iterator begin,
  102. LHCb::TbHits::const_iterator end,
  103. std::vector<bool>& used) {
  104.  
  105. // Keep track of the cluster's bounding box.
  106. int scolmin = pixels.front()->scol();
  107. int scolmax = scolmin;
  108. int rowmin = pixels.front()->row();
  109. int rowmax = rowmin;
  110. // Try adding hits to the cluster.
  111. bool hitAdded = true;
  112. while (hitAdded) {
  113. hitAdded = false;
  114. unsigned int nUnused = 0;
  115. for (LHCb::TbHits::const_iterator it = begin; it != end; ++it) {
  116. // Skip hits which are already part of a cluster.
  117. if (used[(*it)->key()]) continue;
  118. ++nUnused;
  119. const int scol = (*it)->scol();
  120. const int row = (*it)->row();
  121. if (scol < scolmin - m_searchDist) continue;
  122. if (scol > scolmax + m_searchDist) continue;
  123. if (row < rowmin - m_searchDist) continue;
  124. if (row > rowmax + m_searchDist) continue;
  125. // Ask if the hit is touching the cluster (with its current set of hits).
  126. if (hitTouchesCluster((*it)->scol(), (*it)->row(), pixels)) {
  127. // Add the hit to the cluster.
  128. pixels.push_back(*it);
  129. used[(*it)->key()] = true;
  130. --nUnused;
  131. hitAdded = true;
  132. // Update the bounding box.
  133. if (scol < scolmin)
  134. scolmin = scol;
  135. else if (scol > scolmax)
  136. scolmax = scol;
  137. if (row < rowmin)
  138. rowmin = row;
  139. else if (row > rowmax)
  140. rowmax = row;
  141. }
  142. }
  143. if (nUnused == 0) break;
  144. }
  145. }
  146.  
  147. //=============================================================================
  148. // Complete remaining cluster attributes
  149. //=============================================================================
  150. void TbClustering::completeCluster(
  151. const unsigned int plane, const std::vector<const LHCb::TbHit*>& pixels,
  152. LHCb::TbClusters* clusters) {
  153.  
  154. // Create a new cluster object.
  155. LHCb::TbCluster* cluster = new LHCb::TbCluster();
  156. cluster->setPlane(plane);
  157. // Set cluster width along the column and row direction.
  158. auto cols = std::minmax_element(pixels.cbegin(), pixels.cend(),
  159. [](const LHCb::TbHit* h1, const LHCb::TbHit* h2) {
  160. return h1->scol() < h2->scol();
  161. });
  162. auto rows = std::minmax_element(pixels.cbegin(), pixels.cend(),
  163. [](const LHCb::TbHit* h1, const LHCb::TbHit* h2) {
  164. return h1->row() < h2->row();
  165. });
  166. const unsigned int nCols =
  167. 1 + (*cols.second)->scol() - (*cols.first)->scol();
  168. const unsigned int nRows = 1 + (*rows.second)->row() - (*rows.first)->row();
  169. cluster->setCols(nCols);
  170. cluster->setRows(nRows);
  171. // Add the pixel hits to the cluster, sum up the charge and
  172. // calculate the centre of gravity.
  173. unsigned int tot = 0;
  174. double charge = 0.;
  175. double xLocal = 0.;
  176. double yLocal = 0.;
  177. for (auto it = pixels.cbegin(), end = pixels.cend(); it != end; ++it) {
  178. cluster->addToHits(*it);
  179. tot += (*it)->ToT();
  180. const double q = (*it)->charge();
  181. charge += q;
  182. double x = 0.;
  183. double y = 0.;
  184. geomSvc()->pixelToPoint((*it)->scol(), (*it)->row(), plane, x, y);
  185. xLocal += x * q;
  186. yLocal += y * q;
  187. }
  188. cluster->setToT(tot);
  189. cluster->setCharge(charge);
  190. // Assume that the cluster time is the time of the earliest hit.
  191. cluster->setTime(pixels.front()->time());
  192. cluster->setHtime(pixels.front()->htime());
  193. // Calculate the local and global coordinates.
  194. xLocal /= charge;
  195. yLocal /= charge;
  196. // Apply eta-corrections.
  197. etaCorrection(xLocal, yLocal, nCols, nRows, plane);
  198. cluster->setXloc(xLocal);
  199. cluster->setYloc(yLocal);
  200. const Gaudi::XYZPoint pLocal(xLocal, yLocal, 0.);
  201. const auto pGlobal = geomSvc()->localToGlobal(pLocal, plane);
  202. cluster->setX(pGlobal.x());
  203. cluster->setY(pGlobal.y());
  204. cluster->setZ(pGlobal.z());
  205. // Assign error estimates.
  206. setClusterError(cluster);
  207. // Add the cluster to the container.
  208. clusters->insert(cluster);
  209. }
  210.  
  211. //=============================================================================
  212. // Calculate the cluster error estimate.
  213. //=============================================================================
  214. void TbClustering::setClusterError(LHCb::TbCluster* cluster) const {
  215.  
  216. // Initialise the errors to some default value (for large cluster widths).
  217. double dx = 0.1;
  218. double dy = 0.1;
  219. switch (m_clusterErrorMethod) {
  220. case 0: {
  221. constexpr std::array<double, 5> sigmas = {0.0027, 0.0035, 0.016, 0.045, 0.065};
  222. if (cluster->cols() < 6) dx = sigmas[cluster->cols() - 1];
  223. if (cluster->rows() < 6) dy = sigmas[cluster->rows() - 1];
  224. break;
  225. }
  226. case 2:
  227. if (cluster->size() <= 6) {
  228. dx = dy = 0.004;
  229. } else {
  230. // HS: not sure this makes sense
  231. const double sigmaHit = (1. / sqrt(12.)) * Tb::PixelPitch;
  232. dx = sigmaHit * cluster->cols();
  233. dy = sigmaHit * cluster->rows();
  234. }
  235. break;
  236. default:
  237. dx = dy = 0.004;
  238. break;
  239. }
  240. cluster->setXErr(dx);
  241. cluster->setYErr(dy);
  242. cluster->setWx(1. / (dx * dx));
  243. cluster->setWy(1. / (dy * dy));
  244. }
  245.  
  246.  
  247. //=============================================================================
  248. // Read the eta correction parameters from file.
  249. //=============================================================================
  250. void TbClustering::readEta(const std::string& filename) {
  251.  
  252. TbCondFile f(filename);
  253. if (!f.is_open()) {
  254. warning() << "Cannot open " << filename << endmsg;
  255. return;
  256. }
  257. unsigned int indexPlane = 999;
  258. unsigned int indexXy = 0;
  259. unsigned int indexCol = 0;
  260. unsigned int nLines = 0;
  261. while (!f.eof()) {
  262. std::string line = "";
  263. ++nLines;
  264. if (!f.getLine(line)) continue;
  265. if (line.find("Plane") != std::string::npos) {
  266. indexPlane = 999;
  267. std::string key = "";
  268. std::string id = "";
  269. std::string xy = "";
  270. int cols = 0;
  271. if (!f.split(line, ' ', key, id, xy, cols)) {
  272. warning() << "Error reading line " << nLines << endmsg;
  273. continue;
  274. }
  275. // Find the index corresponding to the given plane name.
  276. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  277. if (id == geomSvc()->modules()[i]->id()) {
  278. indexPlane = i;
  279. break;
  280. }
  281. }
  282. if (indexPlane == 999) {
  283. warning() << "Module " << id
  284. << " is not in the alignment file" << endmsg;
  285. continue;
  286. }
  287. // Check whether this set of parameters is for the x or the y direction.
  288. if (xy == "x" || xy == "X") {
  289. indexXy = 0;
  290. } else if (xy == "y" || xy == "Y") {
  291. indexXy = 1;
  292. } else {
  293. warning() << "Unexpected direction specifier " << xy
  294. << " (expected x or y). Skip this parameter set." << endmsg;
  295. indexPlane = 999;
  296. continue;
  297. }
  298. // Check the cluster width.
  299. if (cols < 2) {
  300. warning() << "Unexpected cluster width (" << cols
  301. << "). Skip this parameter set." << endmsg;
  302. indexPlane = 999;
  303. continue;
  304. }
  305. indexCol = cols - 2;
  306. if (m_eta[indexPlane][indexXy].size() < indexCol + 1) {
  307. m_eta[indexPlane][indexXy].resize(indexCol + 1);
  308. }
  309. const std::string rowcol = indexXy == 0 ? " columns." : " rows.";
  310. info() << "Reading eta-correction parameters for plane " << indexPlane
  311. << " for clusters covering " << cols << rowcol << endmsg;
  312. m_eta[indexPlane][indexXy][indexCol].clear();
  313. } else {
  314. if (indexPlane == 999) continue;
  315. // Read the intra-pixel interval and the coefficients of the polynomial.
  316. double xmin = 0., xmax = 0.;
  317. double a = 0., b = 0., c = 0.;
  318. if (!f.split(line, ' ', xmin, xmax, a, b, c)) {
  319. warning() << "Error reading line " << nLines << endmsg;
  320. continue;
  321. }
  322. EtaCorrection corr;
  323. corr.xmin = xmin;
  324. corr.xmax = xmax;
  325. corr.coefficients = {a, b, c};
  326. m_eta[indexPlane][indexXy][indexCol].push_back(corr);
  327. }
  328. }
  329.  
  330. }
  331.  
  332. //=============================================================================
  333. // Calculate the eta correction terms.
  334. //=============================================================================
  335. void TbClustering::etaCorrection(double& xLocal, double& yLocal,
  336. const unsigned int nCols,
  337. const unsigned int nRows,
  338. const unsigned int plane) const {
  339.  
  340. // No point to continue in case of single-pixel clusters.
  341. if (nCols * nRows == 1) return;
  342. // Calculate the intra-pixel coordinates of the cluster.
  343. unsigned int scol = 0, row = 0;
  344. geomSvc()->pointToPixel(xLocal, yLocal, plane, scol, row);
  345. double x0 = 0., y0 = 0.;
  346. geomSvc()->pixelToPoint(scol, row, plane, x0, y0);
  347. // TODO: this only works for single-chip assemblies.
  348. const double x = xLocal - x0 + 0.5 * Tb::PixelPitch;
  349. const double y = yLocal - y0 + 0.5 * Tb::PixelPitch;
  350. // Calculate the correction terms.
  351. if (nCols > 1 && m_eta[plane][0].size() > nCols - 2) {
  352. const auto& corrections = m_eta[plane][0][nCols - 2];
  353. for (const auto& corr : corrections) {
  354. if (x >= corr.xmin && x < corr.xmax) {
  355. xLocal += corr.coefficients[0] + corr.coefficients[1] * x +
  356. corr.coefficients[2] * x * x;
  357. break;
  358. }
  359. }
  360. }
  361. if (nRows > 1 && m_eta[plane][1].size() > nRows - 2) {
  362. const auto& corrections = m_eta[plane][1][nRows - 2];
  363. for (const auto& corr : corrections) {
  364. if (y >= corr.xmin && y < corr.xmax) {
  365. yLocal += corr.coefficients[0] + corr.coefficients[1] * y +
  366. corr.coefficients[2] * y * y;
  367. break;
  368. }
  369. }
  370. }
  371. }