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 surfaceZCenter = 52.89;                 // this is where G4 thinks the surface center is in cm
0075   double zdriftlength = tbin * clockPeriod * _drift_velocity;  // cm
0076   double zloc = surfaceZCenter - zdriftlength;                   // local z relative to surface center (for north side):
0077   unsigned int side = TpcDefs::getSide(hitsetkey);
0078   if (side == 0)
0079   {
0080     zloc = -zloc;
0081   }
0082
0083   auto x = rad * std::cos(phi);
0084   auto y = rad * std::sin(phi);
0085   auto z = surfaceZCenter + zloc;
0086   glob.x() = x;
0087   glob.y() = y;
0088   glob.z() = z;
0089   return glob;
0090 }
0091
0092 //________________________________________________________________________________________________
0093 Acts::Vector3 ActsGeometry::getGlobalPositionTpc(TrkrDefs::cluskey key, TrkrCluster* cluster) const
0094 {
0095   Acts::Vector3 glob;
0096
0097   // This method is for the TPC only
0098   const auto trkrid = TrkrDefs::getTrkrId(key);
0099   if (trkrid != TrkrDefs::tpcId)
0100   {
0101     std::cout << "ActsGeometry::getGlobalPositionTpc -  this is the wrong global transform for silicon or MM's clusters! Returning zero" << std::endl;
0102     return glob;
0103   }
0104
0105   auto surface = m_surfMaps.getSurface(key, cluster);
0106
0107   if (!surface)
0108   {
0109     std::cerr << "Couldn't identify cluster surface. Returning NAN"
0110               << std::endl;
0111     glob(0) = NAN;
0112     glob(1) = NAN;
0113     glob(2) = NAN;
0114     return glob;
0115   }
0116
0117   double surfaceZCenter = 52.89;                                 // this is where G4 thinks the surface center is in cm
0118   double zdriftlength = cluster->getLocalY() * _drift_velocity;  // cm
0119   double zloc = surfaceZCenter - zdriftlength;                   // local z relative to surface center (for north side):
0120   unsigned int side = TpcDefs::getSide(key);
0121   if (side == 0)
0122   {
0123     zloc = -zloc;
0124   }
0125   Acts::Vector2 local(cluster->getLocalX(), zloc);
0126   glob = surface->localToGlobal(m_tGeometry.getGeoContext(),
0127                                 local * Acts::UnitConstants::cm,
0128                                 Acts::Vector3(1, 1, 1));
0129   glob /= Acts::UnitConstants::cm;
0130
0131   return glob;
0132 }
0133
0134 Surface ActsGeometry::get_tpc_surface_from_coords(
0135     TrkrDefs::hitsetkey hitsetkey,
0136     Acts::Vector3 world,
0137     TrkrDefs::subsurfkey& subsurfkey) const
0138 {
0139   unsigned int layer = TrkrDefs::getLayer(hitsetkey);
0140   unsigned int side = TpcDefs::getSide(hitsetkey);
0141
0142   auto mapIter = m_surfMaps.m_tpcSurfaceMap.find(layer);
0143
0144   if (mapIter == m_surfMaps.m_tpcSurfaceMap.end())
0145   {
0146     std::cout << "Error: hitsetkey not found in ActsGeometry::get_tpc_surface_from_coords, hitsetkey = "
0147               << hitsetkey << std::endl;
0148     return nullptr;
0149   }
0150   double world_phi = atan2(world[1], world[0]);
0151
0152   const auto& surf_vec = mapIter->second;
0153   unsigned int surf_index = 999;
0154
0155   // Predict which surface index this phi and side will correspond to
0156   // assumes that the vector elements are ordered positive z, -pi to pi, then negative z, -pi to pi
0157   // we use TPC side from the hitsetkey, since z can be either sign in northa nd south, depending on crossing
0158   double fraction = (world_phi + M_PI) / (2.0 * M_PI);
0159
0160   double rounded_nsurf = std::round((double) (surf_vec.size() / 2) * fraction - 0.5);  // NOLINT
0161   unsigned int nsurfm = (unsigned int) rounded_nsurf;
0162
0163   if (side == 0)
0164   {
0165     nsurfm += surf_vec.size() / 2;
0166   }
0167
0168   unsigned int nsurf = nsurfm % surf_vec.size();
0169   Surface this_surf = surf_vec[nsurf];
0170
0171   auto vec3d = this_surf->center(m_tGeometry.getGeoContext());
0172   std::vector<double> surf_center = {vec3d(0) / 10.0, vec3d(1) / 10.0, vec3d(2) / 10.0};  // convert from mm to cm
0173   double surf_phi = atan2(surf_center[1], surf_center[0]);
0174   double surfStepPhi = m_tGeometry.tpcSurfStepPhi;
0175
0176   if ((world_phi > surf_phi - surfStepPhi / 2.0 && world_phi < surf_phi + surfStepPhi / 2.0))
0177   {
0178     surf_index = nsurf;
0179     subsurfkey = nsurf;
0180   }
0181   else
0182   {
0183     // check for the periodic boundary condition
0184     auto firstsurf = *surf_vec.begin();
0185     auto firstsurfcenter = firstsurf->center(geometry().getGeoContext());
0186     float firstsurf_phi = atan2(firstsurfcenter[1], firstsurfcenter[0]);
0187     if (world_phi < firstsurf_phi - surfStepPhi / 2.0)
0188     {
0189       world_phi += 2.0 * M_PI;
0190     }
0191     //check a few surfaces around this one
0192     for( int i = -1; i <= 1; i++)
0193     {
0194       if(i==0) // already tried this one
0195       {
0196         continue;
0197       }
0198       unsigned int new_nsurf = (nsurf+i) % surf_vec.size();
0199       this_surf = surf_vec[new_nsurf];
0200       vec3d = this_surf->center(geometry().getGeoContext());
0201       surf_center = {vec3d(0) / 10.0, vec3d(1) / 10.0, vec3d(2) / 10.0};  // convert from mm to cm
0202       surf_phi = atan2(surf_center[1], surf_center[0]);
0203
0204       if ((world_phi > surf_phi - surfStepPhi / 2.0 && world_phi < surf_phi + surfStepPhi / 2.0))
0205       {
0206         surf_index = new_nsurf;
0207         subsurfkey = new_nsurf;
0208         return surf_vec[surf_index];
0209       }
0210     }
0211     return nullptr;
0212   }
0213
0214   return surf_vec[surf_index];
0215 }
0216
0217 //________________________________________________________________________________________________
0218 Acts::Transform3 ActsGeometry::makeAffineTransform(Acts::Vector3 rot, Acts::Vector3 trans) const
0219 {
0220   Acts::Transform3 actsAffine;
0221
0222   Eigen::AngleAxisd alpha(rot(0), Eigen::Vector3d::UnitX());
0223   Eigen::AngleAxisd beta(rot(1), Eigen::Vector3d::UnitY());
0224   Eigen::AngleAxisd gamma(rot(2), Eigen::Vector3d::UnitZ());
0225   Eigen::Quaternion<double> q = gamma * beta * alpha;
0226   actsAffine.linear() = q.matrix();
0227
0228   Eigen::Vector3d translation(trans(0), trans(1), trans(2));
0229   actsAffine.translation() = translation;
0230
0231   return actsAffine;
0232 }
0233
0234 //________________________________________________________________________________________________
0235 Acts::Vector2 ActsGeometry::getLocalCoords(TrkrDefs::cluskey key, TrkrCluster* cluster) const
0236 {
0237   short int crossing = 0;
0238   Acts::Vector2 local = getLocalCoords(key, cluster, crossing);
0239   return local;
0240 }
0241
0242
0243 //________________________________________________________________________________________________
0244 Acts::Vector2 ActsGeometry::getLocalCoords(TrkrDefs::cluskey key, TrkrCluster* cluster, short int crossing) const
0245 {
0246   Acts::Vector2 local;
0247
0248   const auto trkrid = TrkrDefs::getTrkrId(key);
0249   if (trkrid == TrkrDefs::tpcId)
0250   {
0251     double crossing_tzero_correction = crossing * sphenix_constants::time_between_crossings;
0252     double surfaceZCenter = 52.89;                                 // this is where G4 thinks the surface center is in cm
0253     double zdriftlength = (cluster->getLocalY() - crossing_tzero_correction) * _drift_velocity;  // cm
0254     double zloc = surfaceZCenter - zdriftlength;         // local z relative to surface center (for north side):
0255     unsigned int side = TpcDefs::getSide(key);
0256     if (side == 0)
0257     {
0258       zloc = -zloc;
0259     }
0260     local(0) = cluster->getLocalX();
0261     local(1) = zloc;
0262   }
0263   else
0264   {
0265     local(0) = cluster->getLocalX();
0266     local(1) = cluster->getLocalY();
0267   }
0268
0269   return local;
0270 }