packages/trackbase/ActsGeometry.cc

0001 #include "ActsGeometry.h"
0002 #include <Acts/Definitions/Algebra.hpp>
0003 #include "TpcDefs.h"
0004 #include "TrkrCluster.h"
0005 #include "alignmentTransformationContainer.h"
0006 #include <phool/sphenix_constants.h>
0007
0008 namespace
0009 {
0010   /// square
0011   template <class T>
0012   inline constexpr T square(const T& x)
0013   {
0014     return x * x;
0015   }
0016
0017   /// get radius from coordinates
0018   template <class T>
0019   T radius(const T& x, const T& y)
0020   {
0021     return std::sqrt(square(x) + square(y));
0022   }
0023 }  // namespace
0024
0025 //________________________________________________________________________________________________
0026 Acts::Vector3 ActsGeometry::getGlobalPosition(TrkrDefs::cluskey key, TrkrCluster* cluster) const
0027 {
0028   Acts::Vector3 glob;
0029
0030   const auto trkrid = TrkrDefs::getTrkrId(key);
0031   if (trkrid == TrkrDefs::tpcId)
0032   {
0033     return getGlobalPositionTpc(key, cluster);
0034   }
0035
0036   /// If silicon/TPOT, the transform is one-to-one since the surface is planar
0037
0038   auto surface = m_surfMaps.getSurface(key, cluster);
0039
0040   if (!surface)
0041   {
0042     std::cerr << "Couldn't identify cluster surface. Returning NAN"
0043               << std::endl;
0044     glob(0) = NAN;
0045     glob(1) = NAN;
0046     glob(2) = NAN;
0047     return glob;
0048   }
0049
0050   Acts::Vector2 local(cluster->getLocalX(), cluster->getLocalY());
0051   Acts::Vector3 global;
0052   global = surface->localToGlobal(m_tGeometry.getGeoContext(),
0053                                   local * Acts::UnitConstants::cm,
0054                                   Acts::Vector3(1, 1, 1));
0055   global /= Acts::UnitConstants::cm;
0056
0057   return global;
0058 }
0059
0060 //________________________________________________________________________________________________
0061 Acts::Vector3 ActsGeometry::getGlobalPositionTpc(const TrkrDefs::hitsetkey& hitsetkey,
0062 const TrkrDefs::hitkey& hitkey, const float& phi, const float& rad,
0063 const float& clockPeriod) const
0064 {
0065   Acts::Vector3 glob;
0066   const auto trkrid = TrkrDefs::getTrkrId(hitsetkey);
0067   if (trkrid != TrkrDefs::tpcId)
0068   {
0069     std::cout << "ActsGeometry::getGlobalPositionTpc -  this is the wrong global transform for silicon or MM's clusters! Returning zero" << std::endl;
0070     return glob;
0071   }
0072
0073   auto tbin = TpcDefs::getTBin(hitkey);
0074   double zdriftlength = tbin * clockPeriod * _drift_velocity;  // cm
0075   double zloc =  _max_driftlength / 2.0 - zdriftlength;                   // local z relative to surface center (for north side):
0076   unsigned int side = TpcDefs::getSide(hitsetkey);
0077   if (side == 0)
0078   {
0079     zloc = -zloc;
0080   }
0081
0082   float surfaceZCenter = _max_driftlength/2.0 + _CM_halfwidth;
0083
0084   auto x = rad * std::cos(phi);
0085   auto y = rad * std::sin(phi);
0086   auto z = surfaceZCenter + zloc;
0087   glob.x() = x;
0088   glob.y() = y;
0089   glob.z() = z;
0090   return glob;
0091 }
0092
0093 //________________________________________________________________________________________________
0094 Acts::Vector3 ActsGeometry::getGlobalPositionTpc(TrkrDefs::cluskey key, TrkrCluster* cluster) const
0095 {
0096   Acts::Vector3 glob;
0097
0098   // This method is for the TPC only
0099   const auto trkrid = TrkrDefs::getTrkrId(key);
0100   if (trkrid != TrkrDefs::tpcId)
0101   {
0102     std::cout << "ActsGeometry::getGlobalPositionTpc -  this is the wrong global transform for silicon or MM's clusters! Returning zero" << std::endl;
0103     return glob;
0104   }
0105
0106   auto surface = m_surfMaps.getSurface(key, cluster);
0107
0108   /*
0109   std::cout << " getGlobalPositionTpc transform is: " << std::endl
0110         <<  surface->transform(m_tGeometry.getGeoContext()).matrix()
0111         << std::endl;
0112   alignmentTransformationContainer::use_alignment = false;
0113   Acts::Vector3 ideal_center = surface->center(m_tGeometry.getGeoContext()) * 0.1;
0114   alignmentTransformationContainer::use_alignment = true;
0115   Acts::Vector3 sensorCenter = surface->center(m_tGeometry.getGeoContext()) * 0.1;  // cm
0116   std::cout << "  ideal surface center: " << ideal_center << std::endl;
0117   std::cout << "  aligned surface center: " << sensorCenter << std::endl;
0118   */
0119
0120   if (!surface)
0121   {
0122     std::cerr << "Couldn't identify cluster surface. Returning NAN"
0123               << std::endl;
0124     glob(0) = NAN;
0125     glob(1) = NAN;
0126     glob(2) = NAN;
0127     return glob;
0128   }
0129
0130   Acts::Vector2 local = getLocalCoords(key, cluster);  // no crossing correction here
0131
0132   glob = surface->localToGlobal(m_tGeometry.getGeoContext(),
0133                                 local * Acts::UnitConstants::cm,
0134                                 Acts::Vector3(1, 1, 1));
0135   glob /= Acts::UnitConstants::cm;
0136
0137
0138   // std::cout << "  local " << local << std::endl;
0139   // std::cout << "  glob " << glob << std::endl;
0140
0141
0142   return glob;
0143 }
0144
0145 Surface ActsGeometry::get_tpc_surface_from_coords(
0146     TrkrDefs::hitsetkey hitsetkey,
0147     Acts::Vector3 world,
0148     TrkrDefs::subsurfkey& subsurfkey) const
0149 {
0150   unsigned int layer = TrkrDefs::getLayer(hitsetkey);
0151   unsigned int side = TpcDefs::getSide(hitsetkey);
0152
0153   auto mapIter = m_surfMaps.m_tpcSurfaceMap.find(layer);
0154
0155   if (mapIter == m_surfMaps.m_tpcSurfaceMap.end())
0156   {
0157     std::cout << "Error: hitsetkey not found in ActsGeometry::get_tpc_surface_from_coords, hitsetkey = "
0158               << hitsetkey << std::endl;
0159     return nullptr;
0160   }
0161   double world_phi = atan2(world[1], world[0]);
0162
0163   const auto& surf_vec = mapIter->second;
0164   unsigned int surf_index = 999;
0165
0166   // Predict which surface index this phi and side will correspond to
0167   // assumes that the vector elements are ordered positive z, -pi to pi, then negative z, -pi to pi
0168   // we use TPC side from the hitsetkey, since z can be either sign in northa nd south, depending on crossing
0169   double fraction = (world_phi + M_PI) / (2.0 * M_PI);
0170
0171   double rounded_nsurf = std::round((double) (surf_vec.size() / 2) * fraction - 0.5);  // NOLINT
0172   unsigned int nsurfm = (unsigned int) rounded_nsurf;
0173
0174   if (side == 0)
0175   {
0176     nsurfm += surf_vec.size() / 2;
0177   }
0178
0179   unsigned int nsurf = nsurfm % surf_vec.size();
0180   Surface this_surf = surf_vec[nsurf];
0181
0182   auto vec3d = this_surf->center(m_tGeometry.getGeoContext());
0183   std::vector<double> surf_center = {vec3d(0) / 10.0, vec3d(1) / 10.0, vec3d(2) / 10.0};  // convert from mm to cm
0184   double surf_phi = atan2(surf_center[1], surf_center[0]);
0185   double surfStepPhi = m_tGeometry.tpcSurfStepPhi;
0186
0187   if ((world_phi > surf_phi - surfStepPhi / 2.0 && world_phi < surf_phi + surfStepPhi / 2.0))
0188   {
0189     surf_index = nsurf;
0190     subsurfkey = nsurf;
0191   }
0192   else
0193   {
0194     // check for the periodic boundary condition
0195     auto firstsurf = *surf_vec.begin();
0196     auto firstsurfcenter = firstsurf->center(geometry().getGeoContext());
0197     float firstsurf_phi = atan2(firstsurfcenter[1], firstsurfcenter[0]);
0198     if (world_phi < firstsurf_phi - surfStepPhi / 2.0)
0199     {
0200       world_phi += 2.0 * M_PI;
0201     }
0202     //check a few surfaces around this one
0203     for( int i = -1; i <= 1; i++)
0204     {
0205       if(i==0) // already tried this one
0206       {
0207         continue;
0208       }
0209       unsigned int new_nsurf = (nsurf+i) % surf_vec.size();
0210       this_surf = surf_vec[new_nsurf];
0211       vec3d = this_surf->center(geometry().getGeoContext());
0212       surf_center = {vec3d(0) / 10.0, vec3d(1) / 10.0, vec3d(2) / 10.0};  // convert from mm to cm
0213       surf_phi = atan2(surf_center[1], surf_center[0]);
0214
0215       if ((world_phi > surf_phi - surfStepPhi / 2.0 && world_phi < surf_phi + surfStepPhi / 2.0))
0216       {
0217         surf_index = new_nsurf;
0218         subsurfkey = new_nsurf;
0219         return surf_vec[surf_index];
0220       }
0221     }
0222     return nullptr;
0223   }
0224
0225   return surf_vec[surf_index];
0226 }
0227
0228 //________________________________________________________________________________________________
0229 Acts::Transform3 ActsGeometry::makeAffineTransform(Acts::Vector3 rot, Acts::Vector3 trans) const
0230 {
0231   Acts::Transform3 actsAffine;
0232
0233   Eigen::AngleAxisd alpha(rot(0), Eigen::Vector3d::UnitX());
0234   Eigen::AngleAxisd beta(rot(1), Eigen::Vector3d::UnitY());
0235   Eigen::AngleAxisd gamma(rot(2), Eigen::Vector3d::UnitZ());
0236   Eigen::Quaternion<double> q = gamma * beta * alpha;
0237   actsAffine.linear() = q.matrix();
0238
0239   Eigen::Vector3d translation(trans(0), trans(1), trans(2));
0240   actsAffine.translation() = translation;
0241
0242   return actsAffine;
0243 }
0244
0245 //________________________________________________________________________________________________
0246 Acts::Vector2 ActsGeometry::getLocalCoords(TrkrDefs::cluskey key, TrkrCluster* cluster) const
0247 {
0248   short int crossing = 0;
0249   Acts::Vector2 local = getLocalCoords(key, cluster, crossing);
0250   return local;
0251 }
0252
0253
0254 //________________________________________________________________________________________________
0255 Acts::Vector2 ActsGeometry::getLocalCoords(TrkrDefs::cluskey key, TrkrCluster* cluster, short int crossing) const
0256 {
0257   Acts::Vector2 local;
0258
0259   const auto trkrid = TrkrDefs::getTrkrId(key);
0260   if (trkrid == TrkrDefs::tpcId)
0261   {
0262     double crossing_tzero_correction = crossing * sphenix_constants::time_between_crossings;
0263     double tcorrected = cluster->getLocalY() +  _tpc_tzero + _sampa_tzero_bias - crossing_tzero_correction;
0264     double zdriftlength = tcorrected * _drift_velocity;
0265     double zloc = _max_driftlength/2.0 - zdriftlength;         // local z relative to surface center (for north side):
0266     unsigned int side = TpcDefs::getSide(key);
0267     if (side == 0)
0268     {
0269       zloc = -zloc;
0270     }
0271     local(0) = cluster->getLocalX();
0272     local(1) = zloc;
0273
0274     /*
0275     std::cout << " clust " << cluster->getLocalY() << " tpc tzero " << _tpc_tzero << " sampa tbias " << _sampa_tzero_bias
0276           << " crossing tzero correction " << crossing_tzero_correction << " corrected clust " << tcorrected
0277           << " drift vel " << _drift_velocity
0278           << " crossing " << crossing << " crossing period " << sphenix_constants::time_between_crossings
0279           << " maxdriftlength " << _max_driftlength << " zdriftlength " << zdriftlength << " zloc " << zloc << std::endl;
0280     */
0281
0282   }
0283   else
0284   {
0285     local(0) = cluster->getLocalX();
0286     local(1) = cluster->getLocalY();
0287   }
0288
0289   return local;
0290 }