Newer
Older
TB_Chris / TbKernel / src / .svn / text-base / TbGeometrySvc.cpp.svn-base
  1. #include <fstream>
  2.  
  3. // ROOT
  4. #include "Math/Translation3D.h"
  5. #include "Math/RotationZYX.h"
  6.  
  7. // Gaudi
  8. #include "GaudiKernel/MsgStream.h"
  9. #include "GaudiKernel/Service.h"
  10.  
  11. // Tb/TbEvent
  12. #include "Event/TbTrack.h"
  13.  
  14. // Local
  15. #include "TbKernel/TbCondFile.h"
  16. #include "TbKernel/ITbDataSvc.h"
  17. #include "TbKernel/TbFunctors.h"
  18. #include "TbKernel/TbConstants.h"
  19. #include "TbGeometrySvc.h"
  20.  
  21. DECLARE_SERVICE_FACTORY(TbGeometrySvc)
  22.  
  23. //============================================================================
  24. // Constructor
  25. //============================================================================
  26. TbGeometrySvc::TbGeometrySvc(const std::string& name, ISvcLocator* svc)
  27. : base_class(name, svc),
  28. m_modules(),
  29. m_moduleIndex(),
  30. m_deviceIndex(),
  31. m_planes(),
  32. m_nDevices(0),
  33. m_msg(nullptr) {}
  34.  
  35. //============================================================================
  36. // Destructor
  37. //============================================================================
  38. TbGeometrySvc::~TbGeometrySvc() {
  39.  
  40. // Delete the modules.
  41. for (auto it = m_modules.begin(), end = m_modules.end(); it != end; ++it) {
  42. if (*it) delete (*it);
  43. }
  44. // Delete the message service.
  45. if (m_msg) delete m_msg;
  46. }
  47.  
  48. //============================================================================
  49. // Initialisation
  50. //============================================================================
  51. StatusCode TbGeometrySvc::initialize() {
  52.  
  53. StatusCode sc = Service::initialize();
  54. if (!sc.isSuccess()) return sc;
  55.  
  56. const std::string& filename =
  57. Gaudi::svcLocator()->service<ITbDataSvc>("TbDataSvc")->getAlignmentFile();
  58. msg() << MSG::INFO << "Importing alignment conditions from " << filename
  59. << endmsg;
  60. // Import geometry conditions from alignment file.
  61. if (!readConditions(filename, m_modules)) {
  62. msg() << MSG::ERROR << "Cannot import alignment conditions" << endmsg;
  63. return StatusCode::FAILURE;
  64. }
  65. // Sort the modules by z-position.
  66. std::sort(m_modules.begin(), m_modules.end(),
  67. TbFunctors::LessByZ<const TbModule*>());
  68. // Loop over the modules.
  69. for (auto it = m_modules.begin(), end = m_modules.end(); it != end; ++it) {
  70. // Map the name of the module to its index in the list.
  71. const unsigned index = it - m_modules.begin();
  72. const std::string name = (*it)->id();
  73. m_moduleIndex[name] = index;
  74. for (auto ic = (*it)->chips().begin(); ic != (*it)->chips().end(); ++ic)
  75. m_planes[ic->id] = index;
  76. }
  77. // Cache the x-coordinates of the pixel centres for sensor tiles.
  78. m_xTriple.resize(3 * Tb::NCols);
  79. for (unsigned int chip = 0; chip < 3; ++chip) {
  80. const double x0 = chip * (Tb::NCols + 2) * Tb::PixelPitch;
  81. const unsigned int offset = chip * Tb::NCols;
  82. for (unsigned int col = 0; col < Tb::NCols; ++col) {
  83. const unsigned int scol = offset + col;
  84. double x = x0 + (0.5 + col) * Tb::PixelPitch;
  85. if (scol == 256 || col == 512) {
  86. x -= 0.5 * Tb::PixelPitch;
  87. } else if (scol == 255 || col == 511) {
  88. x += 0.5 * Tb::PixelPitch;
  89. }
  90. m_xTriple[scol] = x;
  91. }
  92. }
  93. printAlignment(m_modules);
  94. return StatusCode::SUCCESS;
  95. }
  96.  
  97. //============================================================================
  98. // Finalisation
  99. //============================================================================
  100. StatusCode TbGeometrySvc::finalize() { return Service::finalize(); }
  101.  
  102. //============================================================================
  103. // Calculate the local coordinates (pixel centre) of a given pixel
  104. //============================================================================
  105. bool TbGeometrySvc::pixelToPoint(const unsigned int scol,
  106. const unsigned int row,
  107. const unsigned int plane, double& x,
  108. double& y) {
  109.  
  110. if (m_modules[plane]->type() == TbModule::Tpx3) {
  111. x = (scol + 0.5) * Tb::PixelPitch;
  112. y = (row + 0.5) * Tb::PixelPitch;
  113. return true;
  114. } else if (m_modules[plane]->type() == TbModule::Tpx3Triple) {
  115. x = m_xTriple[scol];
  116. y = (row + 0.5) * Tb::PixelPitch;
  117. return true;
  118. }
  119. return false;
  120. }
  121.  
  122. //============================================================================
  123. // Calculate the pixel and inter-pixel position of a given local point
  124. //============================================================================
  125. bool TbGeometrySvc::pointToPixel(const double x, const double y,
  126. const unsigned int plane, unsigned int& scol,
  127. unsigned int& row) {
  128.  
  129. if (x < 0. || y < 0.) return false;
  130. if (m_modules[plane]->type() == TbModule::Tpx3) {
  131. const double fcol = x / Tb::PixelPitch;
  132. const double frow = y / Tb::PixelPitch;
  133. scol = int(fcol);
  134. row = int(frow);
  135. if (scol > 255 || row > 255) return false;
  136. return true;
  137. } else if (m_modules[plane]->type() == TbModule::Tpx3Triple) {
  138. const double chipSize = Tb::NCols * Tb::PixelPitch;
  139. const double interChipDistance = 2 * Tb::PixelPitch;
  140. double x0 = 0.;
  141. for (unsigned int i = 0; i < 3; ++i) {
  142. const double xl = x - x0;
  143. if (xl < chipSize + 0.5 * interChipDistance) {
  144. unsigned int col = 0;
  145. if (xl > 0.) {
  146. col = int(xl / Tb::PixelPitch);
  147. if (col >= Tb::NCols) col = Tb::NCols - 1;
  148. }
  149. scol = col + i * Tb::NCols;
  150. row = int(y / Tb::PixelPitch);
  151. if (row >= Tb::NRows) row = Tb::NRows - 1;
  152. return true;
  153. }
  154. x0 += chipSize + interChipDistance;
  155. }
  156. }
  157. return false;
  158. }
  159.  
  160. //============================================================================
  161. // Calculate intercept of track with detector plane
  162. //============================================================================
  163. Gaudi::XYZPoint TbGeometrySvc::intercept(const LHCb::TbTrack* track,
  164. const unsigned int i) {
  165.  
  166. const unsigned int nStates = track->states().size();
  167. if (nStates > 1) {
  168. auto states = track->states();
  169. const double zPlane = m_modules[i]->centre().z();
  170. if (zPlane < states.front().z()) {
  171. return intercept(states.front().position(), states.front().slopes(), i);
  172. } else if (zPlane > states.back().z()) {
  173. return intercept(states.back().position(), states.back().slopes(), i);
  174. }
  175. for (auto it = states.cbegin(), end = states.cend(); it != end; ++it) {
  176. if (zPlane < (*it).z()) {
  177. const auto p = std::prev(it);
  178. return intercept((*p).position(), (*p).slopes(), i);
  179. }
  180. }
  181. }
  182. return intercept(track->firstState().position(), track->firstState().slopes(),
  183. i);
  184. }
  185.  
  186. //============================================================================
  187. // Import geometry conditions from alignment file
  188. //============================================================================
  189. bool TbGeometrySvc::readConditions(const std::string& filename,
  190. std::vector<TbModule*>& modules) {
  191. if (filename == "") return false;
  192. TbCondFile f(filename);
  193. if (!f.is_open()) return false;
  194. while (!f.eof()) {
  195. std::string line = "";
  196. if (!f.getLine(line)) continue;
  197. std::string id = "";
  198. double dx(0.), dy(0.), dz(0.);
  199. double rx(0.), ry(0.), rz(0.);
  200. f.split(line, ' ', id, dx, dy, dz, rx, ry, rz);
  201.  
  202. TbModule* m = new TbModule();
  203. std::string c1, c2, c3;
  204. f.split(id, ',', c1, c2, c3);
  205. m->setId(id);
  206. m->setAlignment(dx, dy, dz, rx, ry, rz, 0., 0., 0., 0., 0., 0.);
  207. m_deviceIndex[c1] = m_nDevices++;
  208. if (c1 != "" && c2 != "" && c3 != "") {
  209. m->setType(TbModule::Tpx3Triple);
  210. m->addChip(c1);
  211. m->addChip(c2);
  212. m->addChip(c3);
  213. m_deviceIndex[c2] = m_nDevices++;
  214. m_deviceIndex[c3] = m_nDevices++;
  215. } else {
  216. m->setType(TbModule::Tpx3);
  217. m->addChip(id);
  218. }
  219. modules.push_back(m);
  220. }
  221. f.close();
  222. return true;
  223. }
  224.  
  225. //============================================================================
  226. // Return the module for a given chip name
  227. //============================================================================
  228. TbModule* TbGeometrySvc::module(const std::string& id) {
  229.  
  230. if (m_moduleIndex.count(id) < 1) {
  231. msg() << MSG::ERROR << "Module " << id << " not found" << endmsg;
  232. return nullptr;
  233. }
  234.  
  235. return m_modules[m_moduleIndex[id]];
  236. }
  237.  
  238. //============================================================================
  239. // Print the geometry conditions for each module
  240. //============================================================================
  241. void TbGeometrySvc::printAlignment(const std::vector<TbModule*>& modules) {
  242.  
  243. for (auto it = modules.begin(), end = modules.end(); it != end; ++it) {
  244. const std::string name = (*it)->id();
  245. const unsigned index = m_moduleIndex[name];
  246. // Print out the geometry conditions.
  247. msg() << MSG::INFO << format("%2i", index) << format(" %-15s", name.c_str())
  248. << format(" %10.3f", (*it)->x()) << format(" %10.3f", (*it)->y())
  249. << format(" %10.3f", (*it)->z()) << format(" %10.3f", (*it)->rotX())
  250. << format(" %10.3f", (*it)->rotY())
  251. << format(" %10.3f", (*it)->rotZ()) << endmsg;
  252. }
  253. }