#include <fstream> // ROOT #include "Math/Translation3D.h" #include "Math/RotationZYX.h" // Gaudi #include "GaudiKernel/MsgStream.h" #include "GaudiKernel/Service.h" // Tb/TbEvent #include "Event/TbTrack.h" // Local #include "TbKernel/TbCondFile.h" #include "TbKernel/ITbDataSvc.h" #include "TbKernel/TbFunctors.h" #include "TbKernel/TbConstants.h" #include "TbGeometrySvc.h" DECLARE_SERVICE_FACTORY(TbGeometrySvc) //============================================================================ // Constructor //============================================================================ TbGeometrySvc::TbGeometrySvc(const std::string& name, ISvcLocator* svc) : base_class(name, svc), m_modules(), m_moduleIndex(), m_deviceIndex(), m_planes(), m_nDevices(0), m_msg(nullptr) {} //============================================================================ // Destructor //============================================================================ TbGeometrySvc::~TbGeometrySvc() { // Delete the modules. for (auto it = m_modules.begin(), end = m_modules.end(); it != end; ++it) { if (*it) delete (*it); } // Delete the message service. if (m_msg) delete m_msg; } //============================================================================ // Initialisation //============================================================================ StatusCode TbGeometrySvc::initialize() { StatusCode sc = Service::initialize(); if (!sc.isSuccess()) return sc; const std::string& filename = Gaudi::svcLocator()->service<ITbDataSvc>("TbDataSvc")->getAlignmentFile(); msg() << MSG::INFO << "Importing alignment conditions from " << filename << endmsg; // Import geometry conditions from alignment file. if (!readConditions(filename, m_modules)) { msg() << MSG::ERROR << "Cannot import alignment conditions" << endmsg; return StatusCode::FAILURE; } // Sort the modules by z-position. std::sort(m_modules.begin(), m_modules.end(), TbFunctors::LessByZ<const TbModule*>()); // Loop over the modules. for (auto it = m_modules.begin(), end = m_modules.end(); it != end; ++it) { // Map the name of the module to its index in the list. const unsigned index = it - m_modules.begin(); const std::string name = (*it)->id(); m_moduleIndex[name] = index; for (auto ic = (*it)->chips().begin(); ic != (*it)->chips().end(); ++ic) m_planes[ic->id] = index; } // Cache the x-coordinates of the pixel centres for sensor tiles. m_xTriple.resize(3 * Tb::NCols); for (unsigned int chip = 0; chip < 3; ++chip) { const double x0 = chip * (Tb::NCols + 2) * Tb::PixelPitch; const unsigned int offset = chip * Tb::NCols; for (unsigned int col = 0; col < Tb::NCols; ++col) { const unsigned int scol = offset + col; double x = x0 + (0.5 + col) * Tb::PixelPitch; if (scol == 256 || col == 512) { x -= 0.5 * Tb::PixelPitch; } else if (scol == 255 || col == 511) { x += 0.5 * Tb::PixelPitch; } m_xTriple[scol] = x; } } printAlignment(m_modules); return StatusCode::SUCCESS; } //============================================================================ // Finalisation //============================================================================ StatusCode TbGeometrySvc::finalize() { return Service::finalize(); } //============================================================================ // Calculate the local coordinates (pixel centre) of a given pixel //============================================================================ bool TbGeometrySvc::pixelToPoint(const unsigned int scol, const unsigned int row, const unsigned int plane, double& x, double& y) { if (m_modules[plane]->type() == TbModule::Tpx3) { x = (scol + 0.5) * Tb::PixelPitch; y = (row + 0.5) * Tb::PixelPitch; return true; } else if (m_modules[plane]->type() == TbModule::Tpx3Triple) { x = m_xTriple[scol]; y = (row + 0.5) * Tb::PixelPitch; return true; } return false; } //============================================================================ // Calculate the pixel and inter-pixel position of a given local point //============================================================================ bool TbGeometrySvc::pointToPixel(const double x, const double y, const unsigned int plane, unsigned int& scol, unsigned int& row) { if (x < 0. || y < 0.) return false; if (m_modules[plane]->type() == TbModule::Tpx3) { const double fcol = x / Tb::PixelPitch; const double frow = y / Tb::PixelPitch; scol = int(fcol); row = int(frow); if (scol > 255 || row > 255) return false; return true; } else if (m_modules[plane]->type() == TbModule::Tpx3Triple) { const double chipSize = Tb::NCols * Tb::PixelPitch; const double interChipDistance = 2 * Tb::PixelPitch; double x0 = 0.; for (unsigned int i = 0; i < 3; ++i) { const double xl = x - x0; if (xl < chipSize + 0.5 * interChipDistance) { unsigned int col = 0; if (xl > 0.) { col = int(xl / Tb::PixelPitch); if (col >= Tb::NCols) col = Tb::NCols - 1; } scol = col + i * Tb::NCols; row = int(y / Tb::PixelPitch); if (row >= Tb::NRows) row = Tb::NRows - 1; return true; } x0 += chipSize + interChipDistance; } } return false; } //============================================================================ // Calculate intercept of track with detector plane //============================================================================ Gaudi::XYZPoint TbGeometrySvc::intercept(const LHCb::TbTrack* track, const unsigned int i) { const unsigned int nStates = track->states().size(); if (nStates > 1) { auto states = track->states(); const double zPlane = m_modules[i]->centre().z(); if (zPlane < states.front().z()) { return intercept(states.front().position(), states.front().slopes(), i); } else if (zPlane > states.back().z()) { return intercept(states.back().position(), states.back().slopes(), i); } for (auto it = states.cbegin(), end = states.cend(); it != end; ++it) { if (zPlane < (*it).z()) { const auto p = std::prev(it); return intercept((*p).position(), (*p).slopes(), i); } } } return intercept(track->firstState().position(), track->firstState().slopes(), i); } //============================================================================ // Import geometry conditions from alignment file //============================================================================ bool TbGeometrySvc::readConditions(const std::string& filename, std::vector<TbModule*>& modules) { if (filename == "") return false; TbCondFile f(filename); if (!f.is_open()) return false; while (!f.eof()) { std::string line = ""; if (!f.getLine(line)) continue; std::string id = ""; double dx(0.), dy(0.), dz(0.); double rx(0.), ry(0.), rz(0.); f.split(line, ' ', id, dx, dy, dz, rx, ry, rz); TbModule* m = new TbModule(); std::string c1, c2, c3; f.split(id, ',', c1, c2, c3); m->setId(id); m->setAlignment(dx, dy, dz, rx, ry, rz, 0., 0., 0., 0., 0., 0.); m_deviceIndex[c1] = m_nDevices++; if (c1 != "" && c2 != "" && c3 != "") { m->setType(TbModule::Tpx3Triple); m->addChip(c1); m->addChip(c2); m->addChip(c3); m_deviceIndex[c2] = m_nDevices++; m_deviceIndex[c3] = m_nDevices++; } else { m->setType(TbModule::Tpx3); m->addChip(id); } modules.push_back(m); } f.close(); return true; } //============================================================================ // Return the module for a given chip name //============================================================================ TbModule* TbGeometrySvc::module(const std::string& id) { if (m_moduleIndex.count(id) < 1) { msg() << MSG::ERROR << "Module " << id << " not found" << endmsg; return nullptr; } return m_modules[m_moduleIndex[id]]; } //============================================================================ // Print the geometry conditions for each module //============================================================================ void TbGeometrySvc::printAlignment(const std::vector<TbModule*>& modules) { for (auto it = modules.begin(), end = modules.end(); it != end; ++it) { const std::string name = (*it)->id(); const unsigned index = m_moduleIndex[name]; // Print out the geometry conditions. msg() << MSG::INFO << format("%2i", index) << format(" %-15s", name.c_str()) << format(" %10.3f", (*it)->x()) << format(" %10.3f", (*it)->y()) << format(" %10.3f", (*it)->z()) << format(" %10.3f", (*it)->rotX()) << format(" %10.3f", (*it)->rotY()) << format(" %10.3f", (*it)->rotZ()) << endmsg; } }