Newer
Older
TB_Chris / TbAlgorithms / src / .svn / text-base / TbVertexTracking.cpp.svn-base
  1. #include <algorithm>
  2. #include <fstream>
  3.  
  4. // Gaudi
  5. #include "GaudiKernel/PhysicalConstants.h"
  6.  
  7. // Tb/TbKernel
  8. #include "TbKernel/TbFunctors.h"
  9. #include "TbKernel/TbConstants.h"
  10.  
  11. // Local
  12. #include "TbVertexTracking.h"
  13.  
  14. DECLARE_ALGORITHM_FACTORY(TbVertexTracking)
  15.  
  16. //=============================================================================
  17. // Standard constructor
  18. //=============================================================================
  19. TbVertexTracking::TbVertexTracking(const std::string& name,
  20. ISvcLocator* pSvcLocator)
  21. : TbAlgorithm(name, pSvcLocator),
  22. m_tracks(nullptr),
  23. m_trackFit(nullptr),
  24. m_clusterFinder(nullptr),
  25. m_event(0) {
  26.  
  27. // Declare algorithm properties - see header for description.
  28. declareProperty("ClusterLocation",
  29. m_clusterLocation = LHCb::TbClusterLocation::Default);
  30. declareProperty("TrackLocation",
  31. m_trackLocation = LHCb::TbTrackLocation::Default);
  32. declareProperty("TrackFitTool", m_trackFitTool = "TbTrackFit");
  33.  
  34. declareProperty("TimeWindow", m_twindow = 150. * Gaudi::Units::ns);
  35. declareProperty("MinNClusters", m_MinNClusters = 7);
  36. declareProperty("SearchVolumeFillAlgorithm",
  37. m_ClusterFinderSearchAlgorithm = "adap_seq");
  38. declareProperty("SearchPlanes", m_PlaneSearchOrder = {4, 3, 5});
  39. declareProperty("ClusterSizeCut", m_clusterSizeCut = 1000);
  40. declareProperty("MinNClustersRepeat", m_MinNClustersRepeat = 3);
  41. declareProperty("RadialCut", m_radialCut = 0.2);
  42.  
  43. declareProperty("ViewerOutput", m_viewerOutput = false);
  44. declareProperty("ViewerEventNum", m_viewerEvent = 100);
  45. declareProperty("CombatRun", m_combatRun = false);
  46. declareProperty("MaxChi2", m_ChiSqRedCut = 200.);
  47. declareProperty("DoVertexting", m_doVertexting = false);
  48. declareProperty("DoRepeat", m_doRepeat = true);
  49. declareProperty("AngleCut", m_angleCut = 0.2);
  50.  
  51. declareProperty("VertexDelR", m_vertexDelR = 0.1);
  52. declareProperty("VertexDelT", m_vertexDelT = 10);
  53. m_currentAngleCut = m_angleCut;
  54. }
  55.  
  56. //=============================================================================
  57. // Destructor
  58. //=============================================================================
  59. TbVertexTracking::~TbVertexTracking() {}
  60.  
  61. //=============================================================================
  62. // Initialization
  63. //=============================================================================
  64. StatusCode TbVertexTracking::initialize() {
  65.  
  66. // Initialise the base class.
  67. StatusCode sc = TbAlgorithm::initialize();
  68. if (sc.isFailure()) return sc;
  69. // Setup the track fit tool.
  70. m_trackFit = tool<ITbTrackFit>(m_trackFitTool, "Fitter", this);
  71. if (!m_trackFit) return Error("Cannot retrieve track fit tool.");
  72. // Set up the cluster finder.
  73. m_clusterFinder =
  74. tool<ITbClusterFinder>("TbClusterFinder", "ClusterFinder", this);
  75. if (!m_clusterFinder) return Error("Cannot retrieve cluster finder tool.");
  76. m_endCluster.resize(m_nPlanes);
  77. m_vertexedCluster.resize(m_nPlanes);
  78. m_volumed.resize(m_nPlanes);
  79. initialStateVsFitStateTx =
  80. book2D("initialStateVsFitStateTx", "initialStateVsFitStateTx", -0.001,
  81. 0.001, 200, -0.001, 0.001, 200);
  82. initialStateVsFitStateTy =
  83. book2D("initialStateVsFitStateTy", "initialStateVsFitStateTy", -0.001,
  84. 0.001, 200, -0.001, 0.001, 200);
  85. return StatusCode::SUCCESS;
  86. }
  87.  
  88. //=============================================================================
  89. // Main execution
  90. //=============================================================================
  91. StatusCode TbVertexTracking::execute() {
  92. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  93. if (m_event == m_viewerEvent && m_viewerOutput) {
  94. std::ofstream myfile;
  95. myfile.open("/afs/cern.ch/user/d/dsaunder/KeplerViewerData.dat");
  96. myfile << "# Output\n";
  97. myfile.close();
  98. }
  99.  
  100. const std::string clusterLocation = m_clusterLocation + std::to_string(i);
  101. LHCb::TbClusters* clusters = getIfExists<LHCb::TbClusters>(clusterLocation);
  102. if (!clusters) {
  103. error() << "No clusters in " << clusterLocation << endmsg;
  104. return StatusCode::FAILURE;
  105. }
  106. m_endCluster[i].clear();
  107. m_endCluster[i].resize(clusters->size(), false);
  108. m_vertexedCluster[i].clear();
  109. m_vertexedCluster[i].resize(clusters->size(), false);
  110. m_volumed[i].clear();
  111. m_volumed[i].resize(clusters->size(), false);
  112. // Store the cluster iterators in the cluster finder.
  113. m_clusterFinder->setClusters(clusters, i);
  114. }
  115.  
  116. // Check if there is already a track container.
  117. m_tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
  118. if (!m_tracks) {
  119. // Create a track container and transfer its ownership to the TES.
  120. m_tracks = new LHCb::TbTracks();
  121. put(m_tracks, m_trackLocation);
  122. }
  123. m_vertices = new LHCb::TbVertices();
  124. put(m_vertices, LHCb::TbVertexLocation::Default);
  125.  
  126. // Advantagous to filter out the thin tracks.
  127. m_currentAngleCut = 0.01;
  128. performVertexTracking();
  129. m_currentAngleCut = m_angleCut;
  130. performVertexTracking();
  131.  
  132. if (m_doRepeat) {
  133. // Allow lower demands on track size.
  134. unsigned int tempSizeHolder = m_MinNClusters;
  135. m_MinNClusters = m_MinNClustersRepeat;
  136. performVertexTracking();
  137. m_MinNClusters = tempSizeHolder;
  138. }
  139.  
  140. timeOrderTracks();
  141.  
  142. if (m_doVertexting) collectIntoVertices();
  143. /*
  144. LHCb::TbTracks::iterator itrack;
  145. for (itrack = m_tracks->begin(); itrack < m_tracks->end() - 1; itrack++) {
  146. plot((*(itrack + 1))->htime() - (*itrack)->htime(), "TimeBetweenTracks",
  147. "TimeBetweenTracks", 0.0, 5000, 500);
  148. }
  149. */
  150. counter("NumberOfTracks") += m_tracks->size();
  151. counter("NumberOfVertices") += m_vertices->size();
  152. if (m_event == m_viewerEvent && m_viewerOutput) outputViewerData();
  153. m_event++;
  154. return StatusCode::SUCCESS;
  155. }
  156.  
  157. //=============================================================================
  158. // Initialization
  159. //=============================================================================
  160. void TbVertexTracking::collectIntoVertices() {
  161. // Iterate over tracks.
  162. LHCb::TbVertex* tempVertex = nullptr;
  163. LHCb::TbTracks::iterator itrack;
  164. for (itrack = m_tracks->begin(); itrack < m_tracks->end(); itrack++) {
  165. tempVertex = nullptr;
  166. if ((*itrack)->vertexed()) continue;
  167. LHCb::TbTracks::iterator jtrack;
  168. for (jtrack = itrack; jtrack < m_tracks->end(); jtrack++) {
  169. if ((*jtrack)->vertexed() || itrack == jtrack) continue;
  170.  
  171. // Check time difference.
  172. if (fabs((*itrack)->htime() - (*jtrack)->htime()) < m_vertexDelT) {
  173. // Assume decay happened near a detector.
  174. for (unsigned int iplane = 0; iplane < m_nPlanes; iplane++) {
  175. // Work out the spatial separation.
  176. Gaudi::XYZPoint intercept1 = geomSvc()->intercept((*itrack), iplane);
  177. Gaudi::XYZPoint intercept2 = geomSvc()->intercept((*jtrack), iplane);
  178. double delr2 = pow(intercept1.x() - intercept2.x(), 2);
  179. delr2 += pow(intercept1.y() - intercept2.y(), 2);
  180. if (pow(delr2, 0.5) < m_vertexDelR) {
  181.  
  182. // We have a vertex, although should also add some cluster counting.
  183. if (tempVertex == nullptr) {
  184. tempVertex = new LHCb::TbVertex();
  185. tempVertex->addToTracks(*itrack);
  186. (*itrack)->setVertexed(true);
  187. tempVertex->setX(intercept1.x());
  188. tempVertex->setY(intercept1.y());
  189. tempVertex->setZ(intercept1.z());
  190. tempVertex->setHtime((*itrack)->htime());
  191. tempVertex->setInteractionPlane(iplane);
  192.  
  193. SmartRefVector<LHCb::TbCluster>& clusters =
  194. const_cast<SmartRefVector<LHCb::TbCluster>&>(
  195. (*itrack)->clusters());
  196. for (unsigned int i = 0; i < clusters.size(); i++) {
  197. m_vertexedCluster[clusters[i]->plane()][clusters[i]->key()] =
  198. true;
  199. unsigned int plane = clusters[i]->plane();
  200. if (plane == 0 || plane == 1 || plane == 2)
  201. (*itrack)->setParentVertex(true);
  202. }
  203. }
  204.  
  205. tempVertex->addToTracks(*jtrack);
  206. (*jtrack)->setVertexed(true);
  207. SmartRefVector<LHCb::TbCluster>& clusters =
  208. const_cast<SmartRefVector<LHCb::TbCluster>&>(
  209. (*jtrack)->clusters());
  210. for (unsigned int i = 0; i < clusters.size(); i++) {
  211. m_vertexedCluster[clusters[i]->plane()][clusters[i]->key()] =
  212. true;
  213. unsigned int plane = clusters[i]->plane();
  214. if (plane == 0 || plane == 1 || plane == 2)
  215. (*jtrack)->setParentVertex(true);
  216. }
  217. break;
  218. }
  219. }
  220. }
  221. if ((*jtrack)->htime() - (*itrack)->htime() > m_vertexDelT) break;
  222. }
  223. if (tempVertex != nullptr) m_vertices->insert(tempVertex);
  224. }
  225.  
  226. // Remove tracks forming vertices from m_tracks.
  227. LHCb::TbVertices::iterator ivertex;
  228. for (ivertex = m_vertices->begin(); ivertex != m_vertices->end(); ivertex++) {
  229. SmartRefVector<LHCb::TbTrack>& tracks =
  230. const_cast<SmartRefVector<LHCb::TbTrack>&>((*ivertex)->tracks());
  231. for (unsigned int i = 0; i < tracks.size(); i++)
  232. m_tracks->remove(tracks[i]);
  233. }
  234. }
  235.  
  236. //=============================================================================
  237. // Initialization
  238. //=============================================================================
  239. void TbVertexTracking::performVertexTracking() {
  240. // Iterate over search planes.
  241. for (const auto& plane : m_PlaneSearchOrder) {
  242.  
  243. if (masked(plane) || m_clusterFinder->empty(plane)) continue;
  244. auto ic = m_clusterFinder->first(plane);
  245. const auto end = m_clusterFinder->end(plane);
  246. for (; ic != end; ++ic) {
  247. // Check if too big, and if already tracked.
  248. if ((*ic)->size() > m_clusterSizeCut) continue;
  249. if ((*ic)->associated() && !m_endCluster[plane][(*ic)->key()]) continue;
  250. formTrack(*ic);
  251. }
  252. }
  253. }
  254.  
  255. //=============================================================================
  256. // Initialization
  257. //=============================================================================
  258. void TbVertexTracking::evalHoughState(LHCb::TbCluster* seed,
  259. LHCb::TbCluster* cluster2,
  260. LHCb::TbState* tempInitialState) {
  261. if (m_event == m_viewerEvent && m_viewerOutput)
  262. outputHoughState(seed, cluster2);
  263. LHCb::TbTrack* track = new LHCb::TbTrack();
  264. track->addToClusters(seed);
  265. track->addToClusters(cluster2);
  266. track->setFirstState(*tempInitialState);
  267. bool sizeSuitable = fillTrack(track, seed, cluster2);
  268. if (!sizeSuitable) {
  269. counter("NumberOfSizeRejectedTracks") += 1;
  270. delete track;
  271. } else {
  272. m_trackFit->fit(track);
  273. initialStateVsFitStateTx->fill(track->firstState().tx(),
  274. tempInitialState->tx());
  275. initialStateVsFitStateTy->fill(track->firstState().ty(),
  276. tempInitialState->ty());
  277. plot(track->firstState().tx() - tempInitialState->tx(),
  278. "initialStateMinusFittedStateX", "initialStateMinusFittedStateX",
  279. -0.005, 0.005, 200);
  280.  
  281. // Consider popping one clusters if its going to fail.
  282. int popID = -1;
  283. double chi = track->chi2PerNdof();
  284. if (track->clusters().size() > m_MinNClusters + 1 &&
  285. track->chi2PerNdof() > m_ChiSqRedCut) {
  286. for (unsigned int i = 0; i < m_nPlanes; i++) {
  287. m_trackFit->maskPlane(i);
  288. m_trackFit->fit(track);
  289. if (track->chi2PerNdof() < chi) {
  290. chi = track->chi2PerNdof();
  291. popID = i;
  292. }
  293. m_trackFit->unmaskPlane(i);
  294. }
  295. if (popID != -1) {
  296. for (auto cluster : track->clusters()) {
  297. if (int(cluster->plane()) == popID) {
  298. track->removeFromClusters(cluster);
  299. }
  300. }
  301. }
  302. m_trackFit->fit(track);
  303. if (track->chi2PerNdof() < m_ChiSqRedCut)
  304. counter("NumberOfPopRecoveredTracks") += 1;
  305. }
  306.  
  307. if (track->chi2PerNdof() < m_ChiSqRedCut) {
  308. SmartRefVector<LHCb::TbCluster>& clusters =
  309. const_cast<SmartRefVector<LHCb::TbCluster>&>(track->clusters());
  310. auto earliest_hit =
  311. std::min_element(clusters.begin(), clusters.end(),
  312. TbFunctors::LessByTime<const LHCb::TbCluster*>());
  313.  
  314. if (timingSvc()->beforeOverlap((*earliest_hit)->time()) || m_combatRun) {
  315. auto farthest_hit =
  316. std::max_element(clusters.begin(), clusters.end(),
  317. TbFunctors::LessByZ<const LHCb::TbCluster*>());
  318. if ((*farthest_hit)->plane() != m_nPlanes - 1) {
  319. m_endCluster[(*farthest_hit)->plane()][(*farthest_hit)->key()] = true;
  320. }
  321. auto nearest_hit =
  322. std::min_element(clusters.begin(), clusters.end(),
  323. TbFunctors::LessByZ<const LHCb::TbCluster*>());
  324. if ((*nearest_hit)->plane() != 0) {
  325. m_endCluster[(*nearest_hit)->plane()][(*nearest_hit)->key()] = true;
  326. }
  327.  
  328. m_tracks->insert(track);
  329. track->setTime(timingSvc()->localToGlobal(track->htime()));
  330. track->setAssociated(true);
  331. }
  332. } else {
  333. delete track;
  334. counter("NumberOfChiRejectedTracks") += 1;
  335. }
  336. }
  337. }
  338.  
  339. //=============================================================================
  340. //
  341. //=============================================================================
  342. bool TbVertexTracking::fillTrack(LHCb::TbTrack* track,
  343. LHCb::TbCluster* seedCluster,
  344. LHCb::TbCluster* cluster2) {
  345. unsigned int nAllowedGaps = m_nPlanes - m_MinNClusters;
  346. unsigned int nGaps = 0;
  347. for (unsigned int iplane = 0; iplane < m_nPlanes; iplane++) {
  348. bool foundCluster = false;
  349. if (m_clusterFinder->empty(iplane) || masked(iplane) ||
  350. iplane == seedCluster->plane() || iplane == cluster2->plane()) {
  351. // nGaps++;
  352. continue;
  353. }
  354. auto ic =
  355. m_clusterFinder->getIterator(seedCluster->htime() - m_twindow, iplane);
  356. const auto end = m_clusterFinder->end(iplane);
  357. for (; ic != end; ++ic) {
  358. if ((*ic)->size() > m_clusterSizeCut) continue;
  359. Gaudi::XYZPoint intercept = geomSvc()->intercept(track, iplane);
  360. double delr2 = pow(intercept.y() - (*ic)->y(), 2);
  361. delr2 += pow(intercept.x() - (*ic)->x(), 2);
  362.  
  363. plot(intercept.x() - (*ic)->x(), "initialResidualX", "initialResidualX",
  364. -0.05, 0.05, 200);
  365. plot(intercept.y() - (*ic)->y(), "initialResidualY", "initialResidualY",
  366. -0.05, 0.05, 200);
  367.  
  368. if (m_event == m_viewerEvent && m_viewerOutput)
  369. outputPatternRecog(intercept.x(), intercept.y(), intercept.z(),
  370. seedCluster->htime());
  371.  
  372. if (pow(delr2, 0.5) < m_radialCut) {
  373. if (!(*ic)->associated() ||
  374. m_endCluster[(*ic)->plane()][(*ic)->key()]) {
  375. m_volumed[(*ic)->plane()][(*ic)->key()] = true;
  376. track->addToClusters(*ic);
  377. m_trackFit->fit(track);
  378. foundCluster = true;
  379. break;
  380. }
  381. }
  382. if ((*ic)->htime() > seedCluster->htime() + m_twindow) break;
  383. }
  384. if (!foundCluster) nGaps++;
  385. if (nGaps == nAllowedGaps) return false;
  386. }
  387. if (track->clusters().size() < m_MinNClusters) return false;
  388. return true;
  389. }
  390.  
  391. //=============================================================================
  392. //
  393. //=============================================================================
  394. void TbVertexTracking::formTrack(LHCb::TbCluster* seedCluster) {
  395. bool madeTrack = false;
  396. // Form pairs of clusters from either side of seed plane.
  397. // Eval each pair as its formed.
  398. for (unsigned int iplane = seedCluster->plane() - 1;
  399. iplane != seedCluster->plane() + 2; iplane++) {
  400. if (m_clusterFinder->empty(iplane) || iplane == seedCluster->plane() ||
  401. masked(iplane))
  402. continue;
  403. auto ic =
  404. m_clusterFinder->getIterator(seedCluster->htime() - m_twindow, iplane);
  405. const auto end = m_clusterFinder->end(iplane);
  406. for (; ic != end; ++ic) {
  407.  
  408. // Find the gradients between the pair.
  409. if ((*ic)->size() > m_clusterSizeCut || (*ic)->associated()) continue;
  410. double tx =
  411. (seedCluster->x() - (*ic)->x()) / (seedCluster->z() - (*ic)->z());
  412. double ty =
  413. (seedCluster->y() - (*ic)->y()) / (seedCluster->z() - (*ic)->z());
  414. double x0 = seedCluster->x() - tx * seedCluster->z();
  415. double y0 = seedCluster->y() - ty * seedCluster->z();
  416. Gaudi::SymMatrix4x4 cov;
  417. LHCb::TbState fstate(Gaudi::Vector4(x0, y0, tx, ty), cov, 0., 0);
  418. plot(tx, "initialTx", "initialTx", -0.1, 0.1, 200);
  419. plot(ty, "initialTy", "initialTy", -0.1, 0.1, 200);
  420.  
  421. counter("NumberOfFormedHoughStates") += 1;
  422. if (fabs(tx) < 0.01 && fabs(ty) < 0.01) counter("nThinStates") += 1;
  423. if (fabs(tx) < m_currentAngleCut && fabs(ty) < m_currentAngleCut) {
  424. evalHoughState(seedCluster, (*ic), &fstate);
  425. if (seedCluster->associated() &&
  426. !m_endCluster[seedCluster->plane()][seedCluster->key()]) {
  427. madeTrack = true;
  428. break;
  429. }
  430. } else
  431. counter("nAngleCutStates") += 1;
  432. if ((*ic)->htime() > seedCluster->htime() + m_twindow) break;
  433. }
  434. if (madeTrack) break;
  435. }
  436. }
  437.  
  438. //=============================================================================
  439. // Viewer outputs
  440. //=============================================================================
  441. void TbVertexTracking::outputVertices() {
  442. std::ofstream myfile;
  443. myfile.open("/afs/cern.ch/user/d/dsaunder/KeplerViewerData.dat",
  444. std::ofstream::app);
  445. for (LHCb::TbVertex* vertex : *m_vertices) {
  446. myfile << "Vertex " << vertex->x() << " " << vertex->y() << " "
  447. << vertex->z() << " " << vertex->htime() << " "
  448. << vertex->tracks().size() << " ";
  449. for (unsigned int i = 0; i < vertex->tracks().size(); i++) {
  450. myfile << vertex->tracks()[i]->firstState().tx() << " "
  451. << vertex->tracks()[i]->firstState().x() << " "
  452. << vertex->tracks()[i]->firstState().ty() << " "
  453. << vertex->tracks()[i]->firstState().y() << " ";
  454. if (vertex->tracks()[i]->parentVertex())
  455. myfile << "1 ";
  456. else
  457. myfile << "0 ";
  458. }
  459. myfile << "\n";
  460. }
  461. myfile.close();
  462. }
  463.  
  464. //=============================================================================
  465. // Viewer output
  466. //=============================================================================
  467. void TbVertexTracking::outputPatternRecog(double x, double y, double z,
  468. double ht) {
  469. std::ofstream myfile;
  470. myfile.open("/afs/cern.ch/user/d/dsaunder/KeplerViewerData.dat",
  471. std::ofstream::app);
  472. myfile << "PatternRecogCircle " << x << " " << y << " " << z << " " << ht
  473. << " " << m_radialCut << "\n";
  474. myfile.close();
  475. }
  476.  
  477. //=============================================================================
  478. // Viewer output
  479. //=============================================================================
  480. void TbVertexTracking::outputHoughState(LHCb::TbCluster* c1,
  481. LHCb::TbCluster* c2) {
  482. std::ofstream myfile;
  483. myfile.open("/afs/cern.ch/user/d/dsaunder/KeplerViewerData.dat",
  484. std::ofstream::app);
  485. myfile << "HoughState " << c1->x() << " " << c1->y() << " " << c1->z() << " "
  486. << c2->x() << " " << c2->y() << " " << c2->z() << " " << c1->htime()
  487. << " " << c2->htime() << " \n";
  488. myfile.close();
  489. }
  490.  
  491. //=============================================================================
  492. // Viewer output
  493. //=============================================================================
  494. void TbVertexTracking::outputViewerData() {
  495. std::ofstream myfile;
  496. myfile.open("/afs/cern.ch/user/d/dsaunder/KeplerViewerData.dat",
  497. std::ofstream::app);
  498. // First output the chips.
  499. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  500. myfile << "Chip ";
  501. Gaudi::XYZPoint posn1(0., 14.08, 0.);
  502. Gaudi::XYZPoint posn = geomSvc()->localToGlobal(posn1, i);
  503. myfile << posn.x() << " " << posn.y() << " " << posn.z() << " ";
  504.  
  505. Gaudi::XYZPoint posn2(14.08, 14.08, 0.);
  506. posn = geomSvc()->localToGlobal(posn2, i);
  507. myfile << posn.x() << " " << posn.y() << " " << posn.z() << " ";
  508.  
  509. Gaudi::XYZPoint posn3(14.08, 0., 0.);
  510. posn = geomSvc()->localToGlobal(posn3, i);
  511. myfile << posn.x() << " " << posn.y() << " " << posn.z() << " ";
  512.  
  513. Gaudi::XYZPoint posn4(0., 0., 0.);
  514. posn = geomSvc()->localToGlobal(posn4, i);
  515. myfile << posn.x() << " " << posn.y() << " " << posn.z() << " ";
  516. myfile << "\n";
  517. }
  518.  
  519. // Tracks.
  520. for (LHCb::TbTrack* track : *m_tracks) {
  521. myfile << "Track " << track->firstState().tx() << " "
  522. << track->firstState().x() << " " << track->firstState().ty() << " "
  523. << track->firstState().y() << " " << track->htime() << "\n";
  524. }
  525.  
  526. // Clusters.
  527. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  528. const auto end = m_clusterFinder->end(i);
  529. for (auto it = m_clusterFinder->first(i); it != end; ++it) {
  530. myfile << "Cluster " << (*it)->x() << " " << (*it)->y() << " "
  531. << (*it)->z() << " " << (*it)->htime() << " ";
  532. const bool endCluster = m_endCluster[(*it)->plane()][(*it)->key()];
  533. const bool vertexed = m_vertexedCluster[(*it)->plane()][(*it)->key()];
  534. if (endCluster && vertexed) {
  535. myfile << "5 \n";
  536. } else if (vertexed) {
  537. myfile << "4 \n";
  538. } else if (endCluster) {
  539. myfile << "3 \n";
  540. } else if ((*it)->associated()) {
  541. myfile << "2 \n";
  542. } else if (m_volumed[(*it)->plane()][(*it)->key()]) {
  543. myfile << "1 \n";
  544. } else {
  545. myfile << "0 \n";
  546. }
  547.  
  548. // Its hits.
  549. for (auto hit : (*it)->hits()) {
  550. myfile << "Pixel ";
  551. double xLocal = 0.;
  552. double yLocal = 0.;
  553. geomSvc()->pixelToPoint(hit->scol(), hit->row(), i, xLocal, yLocal);
  554. Gaudi::XYZPoint pLocal(xLocal - 0.5 * Tb::PixelPitch,
  555. yLocal - 0.5 * Tb::PixelPitch, 0.);
  556. Gaudi::XYZPoint posn = geomSvc()->localToGlobal(pLocal, i);
  557. myfile << posn.x() << " " << posn.y() << " " << posn.z() << " ";
  558.  
  559. Gaudi::XYZPoint posn2(pLocal.x() + 0.055, pLocal.y(), 0.);
  560. posn = geomSvc()->localToGlobal(posn2, i);
  561. myfile << posn.x() << " " << posn.y() << " " << posn.z() << " ";
  562.  
  563. Gaudi::XYZPoint posn3(pLocal.x() + 0.055, pLocal.y() + 0.055, 0.);
  564. posn = geomSvc()->localToGlobal(posn3, i);
  565. myfile << posn.x() << " " << posn.y() << " " << posn.z() << " ";
  566.  
  567. Gaudi::XYZPoint posn4(pLocal.x(), pLocal.y() + 0.055, 0.);
  568. posn = geomSvc()->localToGlobal(posn4, i);
  569. myfile << posn.x() << " " << posn.y() << " " << posn.z() << " ";
  570.  
  571. myfile << hit->htime() << " " << hit->ToT() << "\n";
  572. }
  573. }
  574. }
  575. myfile.close();
  576. outputVertices();
  577. }
  578.  
  579. //=============================================================================
  580. // Track time ordering - honeycomb
  581. //=============================================================================
  582. void TbVertexTracking::timeOrderTracks() {
  583.  
  584. const double s_factor = 1.3;
  585. LHCb::TbTrack* track1;
  586. LHCb::TbTrack* track2;
  587. unsigned int gap = m_tracks->size() / s_factor;
  588. bool swapped = false;
  589.  
  590. // Start the swap loops.
  591. while (gap > 1 || swapped) {
  592. if (gap > 1) gap /= s_factor;
  593. swapped = false; // Reset per swap loop.
  594.  
  595. // Do the swaps.
  596. LHCb::TbTracks::iterator itt;
  597. for (itt = m_tracks->begin(); itt < m_tracks->end() - gap; ++itt) {
  598. track1 = *itt;
  599. track2 = *(itt + gap);
  600. if (track1->time() > track2->time()) {
  601. // Do the swap.
  602. std::iter_swap(itt, itt + gap);
  603. swapped = true;
  604. }
  605. }
  606. }
  607. }