Newer
Older
TB_Chris / TbAlgorithms / src / .svn / text-base / TbSimpleTracking.cpp.svn-base
  1. // Gaudi
  2. #include "GaudiKernel/PhysicalConstants.h"
  3.  
  4. // Tb/TbKernel
  5. #include "TbKernel/TbFunctors.h"
  6. #include "TbKernel/TbModule.h"
  7.  
  8. // Local
  9. #include "TbSimpleTracking.h"
  10.  
  11. DECLARE_ALGORITHM_FACTORY(TbSimpleTracking)
  12.  
  13. //=============================================================================
  14. // Standard constructor
  15. //=============================================================================
  16. TbSimpleTracking::TbSimpleTracking(const std::string& name,
  17. ISvcLocator* pSvcLocator)
  18. : TbAlgorithm(name, pSvcLocator),
  19. m_trackFit(nullptr) {
  20.  
  21. declareProperty("ClusterLocation",
  22. m_clusterLocation = LHCb::TbClusterLocation::Default);
  23. declareProperty("TrackLocation",
  24. m_trackLocation = LHCb::TbTrackLocation::Default);
  25. declareProperty("TrackFitTool", m_trackFitTool = "TbTrackFit");
  26. declareProperty("TimeWindow", m_timeWindow = 10. * Gaudi::Units::ns);
  27. declareProperty("MinPlanes", m_nMinPlanes = 8);
  28. declareProperty("MaxDistance", m_maxDist = 0. * Gaudi::Units::mm);
  29. declareProperty("MaxOpeningAngle", m_maxAngle = 0.01);
  30. declareProperty("RecheckTrack", m_recheckTrack = true);
  31. declareProperty("RemoveOutliers", m_removeOutliers = true);
  32. declareProperty("ChargeCutLow", m_chargeCutLow = 0);
  33. declareProperty("MaxClusterSize", m_maxClusterSize = 10);
  34. declareProperty("MaxClusterWidth", m_maxClusterWidth = 3);
  35. declareProperty("MaxOccupancy", m_maxOccupancy = 8);
  36. declareProperty("DoOccupancyCut", m_doOccupancyCut = false);
  37. declareProperty("Monitoring", m_monitoring = false);
  38. }
  39.  
  40. //=============================================================================
  41. // Initialization
  42. //=============================================================================
  43. StatusCode TbSimpleTracking::initialize() {
  44. // Initialise the base class.
  45. StatusCode sc = TbAlgorithm::initialize();
  46. if (sc.isFailure()) return sc;
  47.  
  48. m_clusters.resize(m_nPlanes, nullptr);
  49. // Setup the track fit tool.
  50. m_trackFit = tool<ITbTrackFit>(m_trackFitTool, "Fitter", this);
  51. if (!m_trackFit) return Error("Cannot retrieve track fit tool.");
  52.  
  53. m_nMaxPlanes = 0;
  54. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  55. if (!masked(i)) ++m_nMaxPlanes;
  56. }
  57. if (m_nMaxPlanes == 0) return Error("All planes are masked.");
  58. return StatusCode::SUCCESS;
  59. }
  60.  
  61. //=============================================================================
  62. // Main execution
  63. //=============================================================================
  64. StatusCode TbSimpleTracking::execute() {
  65.  
  66. // Get the clusters on each plane.
  67. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  68. if (masked(i)) continue;
  69. const std::string clusterLocation = m_clusterLocation + std::to_string(i);
  70. LHCb::TbClusters* clusters = getIfExists<LHCb::TbClusters>(clusterLocation);
  71. if (!clusters) {
  72. return Error("No clusters in " + LHCb::TbClusterLocation::Default);
  73. }
  74. m_clusters[i] = clusters;
  75. }
  76.  
  77. if (m_doOccupancyCut) findHighOccupancies();
  78. unsigned int nKilledOccupancyCut = 0;
  79. // Check if there is already a track container.
  80. LHCb::TbTracks* tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
  81. if (!tracks) {
  82. // Create a track container and transfer its ownership to the TES.
  83. tracks = new LHCb::TbTracks();
  84. put(tracks, m_trackLocation);
  85. }
  86. LHCb::TbTrack* track = nullptr;
  87. const unsigned int lastPlane = m_nPlanes - 2;
  88. for (unsigned int i = 0; i < lastPlane; ++i) {
  89. if (masked(i)) continue;
  90. // Get the next valid plane.
  91. unsigned int j = i + 1;
  92. while (masked(j) && j <= lastPlane) ++j;
  93. if (j > lastPlane) break;
  94. LHCb::TbClusters::iterator first0 = m_clusters[i]->begin();
  95. LHCb::TbClusters::iterator first1 = m_clusters[j]->begin();
  96. LHCb::TbClusters::iterator end0 = m_clusters[i]->end();
  97. LHCb::TbClusters::iterator end1 = m_clusters[j]->end();
  98. for (LHCb::TbClusters::iterator it0 = first0; it0 != end0; ++it0) {
  99. // Skip clusters already assigned to tracks.
  100. if ((*it0)->associated()) continue;
  101. // Skip clusters with too low ToT.
  102. if ((*it0)->ToT() < m_chargeCutLow) continue;
  103. // Skip large clusters.
  104. if ((*it0)->size() > m_maxClusterSize) continue;
  105. if (((*it0)->cols() > m_maxClusterWidth) ||
  106. ((*it0)->rows() > m_maxClusterWidth)) continue;
  107. const double t0 = (*it0)->htime();
  108. const double tMin = t0 - m_timeWindow;
  109. const double tMax = t0 + m_timeWindow;
  110. // Loop over the hits in the second plane.
  111. for (LHCb::TbClusters::iterator it1 = first1; it1 != end1; ++it1) {
  112. const double t1 = (*it1)->htime();
  113. // Skip hits below the time limit.
  114. if (t1 < tMin) {
  115. // Update the search range.
  116. first1 = it1 + 1;
  117. continue;
  118. }
  119. // Stop search when above the time limit.
  120. if (t1 > tMax) break;
  121. // Skip clusters already assigned to tracks.
  122. if ((*it1)->associated()) continue;
  123. // Skip clusters with too low charge.
  124. if ((*it1)->ToT() < m_chargeCutLow) continue;
  125. // Skip large clusters.
  126. if ((*it1)->size() > m_maxClusterSize) continue;
  127. if (((*it1)->cols() > m_maxClusterWidth) ||
  128. ((*it1)->rows() > m_maxClusterWidth)) continue;
  129. if (!track)
  130. track = new LHCb::TbTrack();
  131. else
  132. track->clearClusters();
  133.  
  134. track->addToClusters(*it0);
  135. track->addToClusters(*it1);
  136. if (!extendTrack(track, true)) continue;
  137. // Fit the track.
  138. m_trackFit->fit(track);
  139. if (m_recheckTrack) recheckTrack(track);
  140. // Check if there are enough measurements on the track.
  141. if (track->clusters().size() < m_nMinPlanes) continue;
  142. // Check if there are enough unused clusters on the track.
  143. /*
  144. unsigned int nUnused = 0;
  145. for (const LHCb::TbCluster* cluster : track->clusters()) {
  146. if (!cluster->associated()) ++nUnused;
  147. }
  148. if (nUnused < 5) continue;
  149. */
  150. if (!lowClusterOccupancy(track->htime())) {
  151. ++nKilledOccupancyCut;
  152. continue;
  153. }
  154. auto clusters = const_cast<SmartRefVector<LHCb::TbCluster>&>(track->clusters());
  155. // Sort the clusters by plane.
  156. std::sort(clusters.begin(), clusters.end(),
  157. TbFunctors::LessByZ<const LHCb::TbCluster*>());
  158. // Tag the clusters as used.
  159. for (auto itc = clusters.begin(), end = clusters.end(); itc != end;
  160. ++itc) {
  161. (*itc)->setAssociated(true);
  162. }
  163. // Calculate the global timestamp and add the track to the container.
  164. track->setTime(timingSvc()->localToGlobal(track->htime()));
  165. tracks->insert(track);
  166. track = nullptr;
  167. }
  168. }
  169. }
  170. if (track) delete track;
  171.  
  172. // Sort the tracks by time.
  173. std::sort(tracks->begin(), tracks->end(),
  174. TbFunctors::LessByTime<const LHCb::TbTrack*>());
  175. // Fill the counters.
  176. counter("Number of tracks") += tracks->size();
  177. counter("Number of tracks removed by occupancy cut") += nKilledOccupancyCut;
  178. appendTrackingEfficiencies();
  179.  
  180. return StatusCode::SUCCESS;
  181. }
  182.  
  183. //=============================================================================
  184. // Extrapolate and try to add clusters to a given seed track
  185. //=============================================================================
  186. bool TbSimpleTracking::extendTrack(LHCb::TbTrack* track, const bool fwd) {
  187.  
  188. unsigned int nAdded = 0;
  189. // Count planes without matching clusters.
  190. unsigned int nMissed = 0;
  191. const LHCb::TbCluster* c1 = track->clusters().front();
  192. const LHCb::TbCluster* c2 = track->clusters().back();
  193. if (!fwd) {
  194. c2 = c1;
  195. c1 = track->clusters().at(1);
  196. }
  197. double t = c2->htime();
  198. // Calculate the track slopes based on the first two clusters on the track.
  199. double td = 1. / (c2->z() - c1->z());
  200. double tx = (c2->x() - c1->x()) * td;
  201. double ty = (c2->y() - c1->y()) * td;
  202. unsigned int plane = c2->plane();
  203. if (fwd) {
  204. plane += 1;
  205. } else {
  206. plane -= 1;
  207. }
  208. while (plane < m_nPlanes) {
  209. // Calculate the extrapolated position on this plane.
  210. const double dz = geomSvc()->modules().at(plane)->z() - c1->z();
  211. const double xPred = c1->x() + tx * dz;
  212. const double yPred = c1->y() + ty * dz;
  213. // Calculate the tolerance window.
  214. const double tol = fabs(dz * m_maxAngle) + m_maxDist;
  215. // Try adding clusters on this plane.
  216. const LHCb::TbCluster* c3 = bestCluster(plane, xPred, yPred, t, tol);
  217. if (c3) {
  218. // Add the cluster.
  219. track->addToClusters(c3);
  220. ++nAdded;
  221. // Reset the missed plane counter.
  222. nMissed = 0;
  223. // Update the pair of clusters to be used for extrapolating.
  224. if ((c3->cols() <= m_maxClusterWidth) &&
  225. (c3->rows() <= m_maxClusterWidth)) {
  226. c1 = c2;
  227. c2 = c3;
  228. // Update the track slopes.
  229. td = 1. / (c2->z() - c1->z());
  230. tx = (c2->x() - c1->x()) * td;
  231. ty = (c2->y() - c1->y()) * td;
  232. }
  233. // Update the predicted timestamp.
  234. t = c3->htime();
  235. } else {
  236. // No matching cluster.
  237. ++nMissed;
  238. if (nMissed > 1) break;
  239. }
  240. if (fwd) {
  241. ++plane;
  242. } else {
  243. if (plane == 0) break;
  244. --plane;
  245. }
  246. }
  247. return nAdded > 0;
  248. }
  249.  
  250. //=============================================================================
  251. // Get the closest cluster to a given point on the plane
  252. //=============================================================================
  253. const LHCb::TbCluster* TbSimpleTracking::bestCluster(const unsigned int plane,
  254. const double xPred, const double yPred,
  255. const double tPred, const double tol) {
  256. if (masked(plane) || m_clusters[plane]->empty()) return nullptr;
  257. // Get the first cluster within the time window.
  258. const double tMin = tPred - m_timeWindow;
  259. LHCb::TbClusters::iterator end = m_clusters[plane]->end();
  260. LHCb::TbClusters::iterator it =
  261. std::lower_bound(m_clusters[plane]->begin(), end, tMin, lowerBound());
  262. if (it == end) return nullptr;
  263. const double tMax = tPred + m_timeWindow;
  264. LHCb::TbCluster* best = nullptr;
  265. double dtBest = m_timeWindow;
  266. for (; it != end; ++it) {
  267. const double t = (*it)->htime();
  268. if ((*it)->ToT() < m_chargeCutLow) continue;
  269. if ((*it)->size() > m_maxClusterSize) continue;
  270. // Stop searching when outside the time window.
  271. if (t > tMax) break;
  272. const double dt = fabs(t - tPred);
  273. if (m_monitoring) {
  274. const std::string title = "Plane" + std::to_string(plane);
  275. plot(t - tPred, "delT" + title, title, -100, 100, 200);
  276. }
  277. if (dt < dtBest) {
  278. // Check if the cluster is within the spatial tolerance window.
  279. const double dx = xPred - (*it)->x();
  280. const double dy = yPred - (*it)->y();
  281. if (m_monitoring) {
  282. const std::string title = "Plane " + std::to_string(plane);
  283. plot(dx, "delXY/" + title, title, -3, 3, 200);
  284. plot(dy, "delXY/" + title, title, -3, 3, 200);
  285. }
  286. if (fabs(dx) > tol || fabs(dy) > tol) continue;
  287. // Update the best cluster.
  288. dtBest = dt;
  289. best = *it;
  290. }
  291. }
  292. return best;
  293. }
  294.  
  295. //=============================================================================
  296. // Remove outliers and try to find additional clusters on a track candidate.
  297. //=============================================================================
  298. void TbSimpleTracking::recheckTrack(LHCb::TbTrack* track) {
  299.  
  300. const unsigned int nClusters = track->clusters().size();
  301. if (nClusters < 3 || nClusters >= m_nMaxPlanes) return;
  302. // Spatial window for adding clusters to the track.
  303. const double tol = 0.1;
  304. // Keep track of the planes which already have a cluster on the track.
  305. std::vector<bool> planeAssociated(m_nPlanes, false);
  306. // Get the straight-line fit parameters.
  307. const double x0 = track->firstState().x();
  308. const double y0 = track->firstState().y();
  309. const double tx = track->firstState().tx();
  310. const double ty = track->firstState().ty();
  311. auto clusters = track->clusters();
  312. const unsigned int nSizeBefore = track->size();
  313. for (auto it = clusters.begin(); it != clusters.end(); ) {
  314. auto cluster = *it;
  315. const unsigned int plane = cluster->plane();
  316. ++it;
  317. // Remove outliers if requested.
  318. const double dx = x0 + tx * cluster->z() - cluster->x();
  319. const double dy = y0 + ty * cluster->z() - cluster->y();
  320. if (m_removeOutliers && (fabs(dx) > tol || fabs(dy) > tol)) {
  321. track->removeFromClusters(cluster);
  322. } else {
  323. planeAssociated[plane] = true;
  324. }
  325. }
  326. const unsigned int nSizeAfter = track->size();
  327. plot(nSizeBefore - nSizeAfter, "NOutliersRemoved", "NOutliersRemoved", -0.5, 10.5, 11);
  328. for (unsigned int i = m_nPlanes; i-- > 0;) {
  329. if (masked(i) || planeAssociated[i]) continue;
  330. // Calculate the track intercept on the plane.
  331. const Gaudi::XYZPoint interceptG = geomSvc()->intercept(track, i);
  332. const LHCb::TbCluster* cluster = bestCluster(
  333. i, interceptG.x(), interceptG.y(), track->htime(), tol);
  334. if (cluster) {
  335. track->addToClusters(cluster);
  336. m_trackFit->fit(track);
  337. }
  338. }
  339. }
  340.  
  341. //=============================================================================
  342.  
  343. void TbSimpleTracking::findHighOccupancies() {
  344.  
  345. m_htimesHighOccupancy.clear();
  346. // Find the first non-masked plane.
  347. unsigned int seedPlane = 0;
  348. while (masked(seedPlane)) ++seedPlane;
  349. for (LHCb::TbClusters::const_iterator iSeed = m_clusters[seedPlane]->begin();
  350. iSeed != m_clusters[seedPlane]->end(); ++iSeed) {
  351. uint nClustersPerSeed = 0;
  352. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  353. if (masked(i)) continue;
  354. for (LHCb::TbClusters::const_iterator it = m_clusters[i]->begin();
  355. it != m_clusters[i]->end(); ++it) {
  356. if (fabs((*iSeed)->htime() - (*it)->htime()) < m_timeWindow)
  357. nClustersPerSeed++;
  358. }
  359. }
  360. if (nClustersPerSeed > m_maxOccupancy)
  361. m_htimesHighOccupancy.push_back((*iSeed)->htime());
  362. }
  363. }
  364.  
  365. //=============================================================================
  366.  
  367. bool TbSimpleTracking::lowClusterOccupancy(const double t) const {
  368. if (!m_doOccupancyCut) return true;
  369. for (unsigned int i = 0; i < m_htimesHighOccupancy.size(); ++i) {
  370. if (fabs(t - m_htimesHighOccupancy[i]) < m_timeWindow) return false;
  371. }
  372. return true;
  373. }
  374.  
  375. //=============================================================================
  376.  
  377. void TbSimpleTracking::appendTrackingEfficiencies() {
  378. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  379. if (masked(i)) continue;
  380. for (LHCb::TbClusters::const_iterator it = m_clusters[i]->begin();
  381. it != m_clusters[i]->end(); ++it) {
  382. counter("nClusters")++;
  383. if ((*it)->size() <= m_maxClusterSize &&
  384. (*it)->ToT() >= m_chargeCutLow &&
  385. lowClusterOccupancy((*it)->htime())) {
  386. counter("nClusters selected for tracking")++;
  387. if ((*it)->associated()) counter("nClusters associated")++;
  388. }
  389. }
  390. }
  391. }
  392.  
  393. //=============================================================================