Implement the ability to tilt the TPC envelope in sPHENIX

offline/QA/Tracking/TpcSeedsQA.cc

@@ -73,7 +73,7 @@ int TpcSeedsQA::InitRun(PHCompositeNode *topNode)
   // global position wrapper
   m_globalPositionWrapper.loadNodes(topNode);
 
-  m_clusterMover.initialize_geometry(g4geom);
+  m_clusterMover.initialize_geometry(g4geom, actsgeom);
   m_clusterMover.set_verbosity(0);
 
   auto *hm = QAHistManagerDef::getHistoManager();

offline/packages/TrackingDiagnostics/TrackResiduals.cc

@@ -118,7 +118,9 @@ int TrackResiduals::InitRun(PHCompositeNode* topNode)
   m_globalPositionWrapper.set_suppressCrossing(m_convertSeeds);
   // clusterMover needs the correct radii of the TPC layers
   auto *tpccellgeo = findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
-  m_clusterMover.initialize_geometry(tpccellgeo);
+
+  auto *geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
+  m_clusterMover.initialize_geometry(tpccellgeo, geometry);
   m_clusterMover.set_verbosity(0);
 
   auto *se = Fun4AllServer::instance();

offline/packages/tpc/TpcClusterMover.cc

@@ -45,13 +45,21 @@ TpcClusterMover::TpcClusterMover()
   }
 }
 
-void TpcClusterMover::initialize_geometry(PHG4TpcGeomContainer* cellgeo)
+void TpcClusterMover::initialize_geometry(PHG4TpcGeomContainer* cellgeo, ActsGeometry *tGeometry)
 {
   if (_verbosity > 0)
   {
-    std::cout << "TpcClusterMover: Initializing layer radii for Tpc from cell geometry object" << std::endl;
+    std::cout << "TpcClusterMover: Getting ActsGeometry, and getting layer radii for Tpc from cell geometry object" << std::endl;
   }
 
+  if(!tGeometry || !cellgeo)
+    {
+      std::cout << PHWHERE << " Failed to get ActsGeometry or TPC cell geometry, cannot continue - quit!" << std::endl;
+	exit(1);
+    }
+
+  _tGeometry = tGeometry;
+
   int layer = 0;
   PHG4TpcGeomContainer::ConstRange layerrange = cellgeo->get_begin_end();
   for (PHG4TpcGeomContainer::ConstIterator layeriter = layerrange.first;
@@ -67,8 +75,9 @@ void TpcClusterMover::initialize_geometry(PHG4TpcGeomContainer* cellgeo)
 std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>> TpcClusterMover::processTrack(const std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>>& global_in)
 {
   // Get the global positions of the TPC clusters for this track, already corrected for distortions, and move them to the surfaces
-  // The input object contains all clusters for the track
-
+  // The input object contains all clusters for the track in world coordinates
+  // The surface radii are in envelope coordinates, we transform the positions to envelope coordinates
+
   std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>> global_moved;
 
   std::vector<Acts::Vector3> tpc_global_vec;
@@ -80,7 +89,8 @@ std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>> TpcClusterMover::proces
     if (trkrid == TrkrDefs::tpcId)
     {
       tpc_cluskey_vec.push_back(ckey);
-      tpc_global_vec.push_back(global);
+      Acts::Vector3 env_global = _tGeometry->transformTpcWorldToEnvelope(global);
+      tpc_global_vec.push_back(env_global);
     }
     else
     {
@@ -140,9 +150,9 @@ std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>> TpcClusterMover::proces
     // now move the cluster to the surface radius
     // we keep the cluster key fixed, change the surface if necessary
 
-    Acts::Vector3 global_new(xnew, ynew, znew);
-
-    // add the new position and surface to the return object
+    Acts::Vector3 env_global_new(xnew, ynew, znew);
+    // now we transform back to global coordinates and add the new position and surface to the return object
+    Acts::Vector3 global_new = _tGeometry->transformTpcEnvelopeToWorld(env_global_new);
     global_moved.emplace_back(cluskey, global_new);
 
     if (_verbosity > 2)

offline/packages/tpc/TpcClusterMover.h

@@ -24,7 +24,7 @@ class TpcClusterMover
 
   //! Updates the assumed default geometry below to that contained in the
   //! cell geo
-  void initialize_geometry(PHG4TpcGeomContainer *cellgeo);
+  void initialize_geometry(PHG4TpcGeomContainer *cellgeo, ActsGeometry *tGeometry);
 
  private:
   int get_circle_circle_intersection(double target_radius, double R, double X0, double Y0, double xclus, double yclus, double &x, double &y) const;
@@ -48,6 +48,8 @@ class TpcClusterMover
   double outer_tpc_spacing = 0.0;
 
   int _verbosity = 0;
+
+ ActsGeometry *_tGeometry  = nullptr;
 };
 
 #endif

offline/packages/tpc/TpcClusterizer.cc

@@ -30,7 +30,8 @@
 #include <fun4all/Fun4AllReturnCodes.h>
 #include <fun4all/SubsysReco.h>  // for SubsysReco
 
-#include <g4detectors/PHG4TpcGeomv1.h>
+//#include <g4detectors/PHG4TpcGeomv1.h>
+#include <g4detectors/PHG4TpcGeom.h>
 #include <g4detectors/PHG4TpcGeomContainer.h>
 
 #include <Acts/Definitions/Units.hpp>
@@ -806,9 +807,9 @@ namespace
     // This is local position
     double clusiphi = iphi_sum / adc_sum;
     double clusit = it_sum / adc_sum;
-
-    // This is the global position
-    double clusphi = my_data.layergeom->get_phi(clusiphi, my_data.side);
+    
+    // this is the phi position in the TPC envelope
+    double env_clusphi = my_data.layergeom->get_phi(clusiphi, my_data.side);    
     double clust = t_sum / adc_sum;
 
     // ADC of centroid bin
@@ -835,34 +836,33 @@ namespace
 
     if (my_data.layergeom->get_phistep() > 0)
     {
-      padphase = (clusphi - maxphi) / my_data.layergeom->get_phistep();
+      padphase = (env_clusphi - maxphi) / my_data.layergeom->get_phistep();
     }
 
     if (my_data.layergeom->get_zstep() > 0)
     {
       tbinphase = (clust - maxt) / my_data.layergeom->get_zstep();
     }
 
-    double clusx = radius * cos(clusphi);
-    double clusy = radius * sin(clusphi);
+    // these positions are in the tpc_envelope
+    double env_clusx = radius * cos(env_clusphi);
+    double env_clusy = radius * sin(env_clusphi);
 
     // needed for surface identification
     double zdriftlength = clust * my_data.tGeometry->get_drift_velocity();
     // convert z drift length to z position in the TPC
-    double clusz = my_data.m_tdriftmax * my_data.tGeometry->get_drift_velocity() - zdriftlength;
+    double env_clusz = my_data.m_tdriftmax * my_data.tGeometry->get_drift_velocity() - zdriftlength;
     if (my_data.side == 0)
     {
-      clusz = -clusz;
+      env_clusz = -env_clusz;
     }
-    //  std::cout << " side " << my_data.side << " clusz " << clusz << " clust " << clust << " driftmax " << my_data.m_tdriftmax << std::endl;
     const double phi_cov = (iphi2_sum / adc_sum - square(clusiphi)) * pow(my_data.layergeom->get_phistep(), 2);
     const double t_cov = t2_sum / adc_sum - square(clust);
 
-    // Get the surface key to find the surface from the
-    // TrkrDefs::hitsetkey tpcHitSetKey = TpcDefs::genHitSetKey(my_data.layer, my_data.sector, my_data.side);
-    Acts::Vector3 global(clusx, clusy, clusz);
+    // get_tpc_surface_from_coords expects a world global position
+    Acts::Vector3 env_global(env_clusx, env_clusy, env_clusz);
+    Acts::Vector3 global = my_data.tGeometry->transformTpcEnvelopeToWorld(env_global);
     TrkrDefs::subsurfkey subsurfkey = 0;
-
     Surface surface = my_data.tGeometry->get_tpc_surface_from_coords(
         tpcHitSetKey,
         global,
@@ -875,6 +875,8 @@ namespace
       hitkeyvec.clear();
       return;
     }
+    // Acts::Vector3 surfcent = surface->center(my_data.tGeometry->geometry().getGeoContext()) / Acts::UnitConstants::cm;
+    //    std::cout << "     surf center = " << surfcent.x() << "  " << surfcent.y() << "  " << surfcent.z() << std::endl;
 
     // Estimate the errors
     // Blow up error on single pixel clusters by a factor 3 to compensate for threshold effects
@@ -961,6 +963,7 @@ namespace
       b_made_cluster = true;
     }
 
+    // This code needs to be reviewed in case of a non-zero TPC tilt - ADF 6/16/26
     if (use_nn && clus_base && training_hits)
     {
       try
@@ -978,9 +981,11 @@ namespace
         double nn_z = training_hits->z + std::clamp(ten_pos[0][1][0].item<double>(), -(double) nd, (double) nd) * training_hits->zstep;
         double nn_x = radius * std::cos(nn_phi);
         double nn_y = radius * std::sin(nn_phi);
-        Acts::Vector3 nn_global(nn_x, nn_y, nn_z);
+
+        Acts::Vector3 nn_env_global(nn_x, nn_y, nn_z);
+        Acts::Vector3 nn_global = my_data.tGeometry->transformTpcEnvelopeToWorld(nn_env_global);
         nn_global *= Acts::UnitConstants::cm;
-        Acts::Vector3 nn_local = surface->localToGlobalTransform(my_data.tGeometry->geometry().geoContext).inverse() * nn_global;
+         Acts::Vector3 nn_local = surface->localToGlobalTransform(my_data.tGeometry->geometry().geoContext).inverse() * nn_global;
         nn_local /= Acts::UnitConstants::cm;
         double nn_t = my_data.m_tdriftmax - std::fabs(nn_z) / my_data.tGeometry->get_drift_velocity();
         clus_base->setLocalX(nn_local(0));
@@ -1524,16 +1529,11 @@ int TpcClusterizer::InitRun(PHCompositeNode *topNode)
 
   AdcClockPeriod = geom->GetFirstLayerCellGeom()->get_zstep();
 
-  std::cout << "FirstLayerCellGeomv1 streamer: " << std::endl;  
-  auto *g1 = static_cast<PHG4TpcGeomv1*> (geom->GetFirstLayerCellGeom()); // cast because << not in the base class
-  std::cout << *g1 << std::endl;
-  std::cout << "LayerCellGeomv1 streamer for layer 24: " << std::endl;
-  auto *g2 = static_cast<PHG4TpcGeomv1*> (geom->GetLayerCellGeom(24)); // cast because << not in the base class
-  std::cout << *g2 << std::endl;
-  std::cout << "LayerCellGeomv1 streamer for layer 40: " << std::endl;  
-  auto *g3 = static_cast<PHG4TpcGeomv1*> (geom->GetLayerCellGeom(40)); // cast because << not in the base class
-  std::cout << *g3 << std::endl;
-
+  // the identify now contains all information from the streamer for v2
+  geom->GetFirstLayerCellGeom()->identify();
+  geom->GetLayerCellGeom(24)->identify();
+  geom->GetLayerCellGeom(40)->identify();
+
   if (m_maskDeadChannels)
   {
     m_deadChannelMap.clear();

offline/packages/trackbase/ActsGeometry.cc

@@ -7,6 +7,10 @@
 
 #include <Acts/Definitions/Algebra.hpp>
 
+#include <Eigen/Dense>
+#include <Eigen/Geometry>
+#include <Eigen/LU>
+
 namespace
 {
   /// square
@@ -149,9 +153,16 @@ Surface ActsGeometry::get_tpc_surface_from_coords(
     Acts::Vector3 world,
     TrkrDefs::subsurfkey& subsurfkey) const
 {
+  // Assume that the world coordinates are in the sPHENIX frame, where the TPC is tilted
+  // We convert the position to tpc envelope coordinates, where we know where everything is
+  Acts::Vector3 world_envelope = transformTpcWorldToEnvelope(world);
+  double world_phi = atan2(world_envelope[1], world_envelope[0]);
+
   unsigned int layer = TrkrDefs::getLayer(hitsetkey);
   unsigned int side = TpcDefs::getSide(hitsetkey);
-
+  unsigned int sector = TpcDefs::getSectorId(hitsetkey);
+
+  // returns an iterator to all of the surfaces for this layer
   auto mapIter = m_surfMaps.m_tpcSurfaceMap.find(layer);
 
   if (mapIter == m_surfMaps.m_tpcSurfaceMap.end())
@@ -160,45 +171,82 @@ Surface ActsGeometry::get_tpc_surface_from_coords(
               << hitsetkey << std::endl;
     return nullptr;
   }
-  double world_phi = atan2(world[1], world[0]);
 
   const auto& surf_vec = mapIter->second;
   unsigned int surf_index = 999;
 
+  // Apparently, tilting the TPC leads to the surfaces not being sorted in phi in the outer layers
+  // just test all surfaces in each layer for now
+  for(unsigned int isurf = 0; isurf < surf_vec.size(); ++isurf)
+    {
+      Surface this_surf = surf_vec[isurf];
+      auto surf_center = this_surf->center(m_tGeometry.getGeoContext());
+      surf_center /= 10.0;   // convert from mm to cm
+      //this surface center includes the TPC tilt used in PHG4TpcDetector construction, transform it to tpc envelope coordinates
+      Acts::Vector3 surf_center_envelope = transformTpcWorldToEnvelope(surf_center);
+      double surf_phi = atan2(surf_center_envelope[1], surf_center_envelope[0]);  
+      double surfStepPhi = m_tGeometry.tpcSurfStepPhi;
+
+      const double dphi = std::atan2(std::sin(world_phi - surf_phi), std::cos(world_phi - surf_phi));
+      if (std::abs(dphi) < surfStepPhi / 2.0)
+	{
+	  if(surf_center_envelope.z() < 0 && side != 0) { continue; }
+	  if(surf_center_envelope.z() > 0 && side != 1) { continue; }
+	  surf_index = isurf;
+	  subsurfkey = isurf;
+	  break;
+	}            
+    }  
+
+  if(surf_index == 999)
+    {
+    std::cout << "Error: surface not found in ActsGeometry::get_tpc_surface_from_coords "
+              << " layer " << layer << " side " << side << " sector " << sector << " world_phi " << world_phi << "  world[0]  " << world[0] << " world[1] " << world[1] << " hitsetkey " << hitsetkey << std::endl;
+    return nullptr;
+    }
+
+  return surf_vec[surf_index];
+
+  /*
   // Predict which surface index this phi and side will correspond to
   // assumes that the vector elements are ordered positive z, -pi to pi, then negative z, -pi to pi
   // we use TPC side from the hitsetkey, since z can be either sign in north and south, depending on crossing
   double fraction = (world_phi + M_PI) / (2.0 * M_PI);
 
   double rounded_nsurf = std::round((double) (surf_vec.size() / 2) * fraction - 0.5);  // NOLINT
   unsigned int nsurfm = (unsigned int) rounded_nsurf;
-
+ std::cout << " surf_vec.size " << surf_vec.size() << " rounded_nsurf " << rounded_nsurf << " initial nsurfm " << nsurfm << std::endl;
+
   if (side == 0)
   {
     nsurfm += surf_vec.size() / 2;
   }
   unsigned int nsurf = nsurfm % surf_vec.size();
   Surface this_surf = surf_vec[nsurf];
-  //std::cout << "    world_phi " << world_phi << " fraction " << fraction << " rounded_nsurf " << rounded_nsurf << " nsurfm " << nsurfm << " nsurf " << nsurf << std::endl;
-
-  auto vec3d = this_surf->center(m_tGeometry.getGeoContext());
-  std::vector<double> surf_center = {vec3d(0) / 10.0, vec3d(1) / 10.0, vec3d(2) / 10.0};  // convert from mm to cm
-  double surf_phi = atan2(surf_center[1], surf_center[0]);
+  auto surf_center = this_surf->center(m_tGeometry.getGeoContext());
+  surf_center /= 10.0;   // convert from mm to cm
+  //this surface center is from the default geometry, which includes the TPC tilt used in PHG4TpcDetector construction
+  // transform it to tpc envelope coordinates
+  Acts::Vector3 surf_center_envelope = m_tpc_world_envelope_transform * surf_center;
+
+  double surf_phi = atan2(surf_center_envelope[1], surf_center_envelope[0]);  
   double surfStepPhi = m_tGeometry.tpcSurfStepPhi;
-  //  std::cout << "    surf_phi " << surf_phi << " surfStepPhi " << surfStepPhi  << " nsurf " << nsurf << std::endl;
 
-  if ((world_phi > surf_phi - surfStepPhi / 2.0 && world_phi < surf_phi + surfStepPhi / 2.0))
+  if ((world_phi > surf_phi - surfStepPhi / 2.0) && (world_phi < surf_phi + surfStepPhi / 2.0))
   {
     surf_index = nsurf;
     subsurfkey = nsurf;
+    std::cout << "success, found nsurf = " << nsurf << std::endl;
   }
   else
   {
     // check for the periodic boundary condition
     auto firstsurf = *surf_vec.begin();
-    auto firstsurfcenter = firstsurf->center(geometry().getGeoContext());
-    float firstsurf_phi = atan2(firstsurfcenter[1], firstsurfcenter[0]);
-    if (world_phi < firstsurf_phi - surfStepPhi / 2.0)
+    auto firstsurfcenter = firstsurf->center(m_tGeometry.getGeoContext());
+    firstsurfcenter /= 10.0;
+    auto firstsurfcenter_envelope = m_tpc_world_envelope_transform * firstsurfcenter;
+    double firstsurf_phi = atan2(firstsurfcenter_envelope[1], firstsurfcenter_envelope[0]);
+    if (world_phi < -M_PI)
     {
       world_phi += 2.0 * M_PI;
     }
@@ -211,21 +259,26 @@ Surface ActsGeometry::get_tpc_surface_from_coords(
       }
       unsigned int new_nsurf = (nsurf+i) % surf_vec.size();
       this_surf = surf_vec[new_nsurf];
-      vec3d = this_surf->center(geometry().getGeoContext());
-      surf_center = {vec3d(0) / 10.0, vec3d(1) / 10.0, vec3d(2) / 10.0};  // convert from mm to cm
-      surf_phi = atan2(surf_center[1], surf_center[0]);
-      //std::cout << "    new world_phi " << world_phi << " new surf_phi " << surf_phi  << " new_nsurf " << new_nsurf << std::endl;
-      if ((world_phi > surf_phi - surfStepPhi / 2.0 && world_phi < surf_phi + surfStepPhi / 2.0))
+      surf_center = this_surf->center(m_tGeometry.getGeoContext());
+      surf_center /= 10.0;
+      surf_center_envelope = m_tpc_world_envelope_transform * surf_center;
+      surf_phi = atan2(surf_center_envelope[1], surf_center_envelope[0]);
+      double this_philow =  surf_phi - surfStepPhi / 2.0;
+      double this_phihigh =  surf_phi + surfStepPhi / 2.0;
+
+      if ((world_phi > this_philow) && (world_phi < this_phihigh))
       {
+	std::cout << "success, found nsurf = " << new_nsurf << std::endl;	
         surf_index = new_nsurf;
         subsurfkey = new_nsurf;
         return surf_vec[surf_index];
       }
     }
     return nullptr;
   }
+  */
+
 
-  return surf_vec[surf_index];
 }
 
 //________________________________________________________________________________________________
@@ -291,3 +344,17 @@ Acts::Vector2 ActsGeometry::getLocalCoords(TrkrDefs::cluskey key, TrkrCluster* c
 
   return local;
 }
+
+  Acts::Vector3  ActsGeometry::transformTpcWorldToEnvelope(const Acts::Vector3& world) const
+  {
+    Acts::Vector3 envelope = m_tpc_world_envelope_transform * world;
+
+    return envelope;
+  }
+
+  Acts::Vector3  ActsGeometry::transformTpcEnvelopeToWorld(const Acts::Vector3& envelope) const
+  {
+    Acts::Vector3 world = m_tpc_world_envelope_transform.inverse() * envelope;
+
+    return world;
+  }