packages/tpc/TpcClusterMover.cc

0001 /*!
0002  * \file TpcClusterMover.cc
0003  * \brief moves distortion corrected clusters back to their TPC surface
0004  * \author Tony Frawley, April 2022
0005  */
0006
0007 #include "TpcClusterMover.h"
0008
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010 #include <trackbase/TrackFitUtils.h>
0011
0012 #include <g4detectors/PHG4TpcCylinderGeom.h>
0013 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0014 #include <climits>
0015 #include <cmath>
0016 #include <iostream>
0017
0018 namespace
0019 {
0020   [[maybe_unused]] std::ostream& operator << (std::ostream& out, const Acts::Vector3& v )
0021   {
0022     out << "(" << v.x() << ", " << v.y() << ", " << v.z() << ")";
0023     return out;
0024   }
0025 }
0026
0027 TpcClusterMover::TpcClusterMover()
0028 {
0029   // initialize layer radii
0030   inner_tpc_spacing = (mid_tpc_min_radius - inner_tpc_min_radius) / 16.0;
0031   mid_tpc_spacing = (outer_tpc_min_radius - mid_tpc_min_radius) / 16.0;
0032   outer_tpc_spacing = (outer_tpc_max_radius - outer_tpc_min_radius) / 16.0;
0033   for (int i = 0; i < 16; ++i)
0034   {
0035     layer_radius[i] = inner_tpc_min_radius + (double) i * inner_tpc_spacing + 0.5 * inner_tpc_spacing;
0036   }
0037   for (int i = 0; i < 16; ++i)
0038   {
0039     layer_radius[i + 16] = mid_tpc_min_radius + (double) i * mid_tpc_spacing + 0.5 * mid_tpc_spacing;
0040   }
0041   for (int i = 0; i < 16; ++i)
0042   {
0043     layer_radius[i + 32] = outer_tpc_min_radius + (double) i * outer_tpc_spacing + 0.5 * outer_tpc_spacing;
0044   }
0045 }
0046
0047 void TpcClusterMover::initialize_geometry(PHG4TpcCylinderGeomContainer *cellgeo)
0048 {
0049   if (_verbosity > 0)
0050   {
0051     std::cout << "TpcClusterMover: Initializing layer radii for Tpc from cell geometry object" << std::endl;
0052   }
0053
0054   int layer = 0;
0055   PHG4TpcCylinderGeomContainer::ConstRange layerrange = cellgeo->get_begin_end();
0056   for (PHG4TpcCylinderGeomContainer::ConstIterator layeriter = layerrange.first;
0057        layeriter != layerrange.second;
0058        ++layeriter)
0059   {
0060     layer_radius[layer] = layeriter->second->get_radius();
0061     layer++;
0062   }
0063 }
0064
0065 //____________________________________________________________________________..
0066 std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>> TpcClusterMover::processTrack(const std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>>& global_in)
0067 {
0068
0069   // Get the global positions of the TPC clusters for this track, already corrected for distortions, and move them to the surfaces
0070   // The input object contains all clusters for the track
0071
0072   std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>> global_moved;
0073
0074   std::vector<Acts::Vector3> tpc_global_vec;
0075   std::vector<TrkrDefs::cluskey> tpc_cluskey_vec;
0076
0077   for (const auto& [ckey,global]:global_in)
0078   {
0079     const auto trkrid = TrkrDefs::getTrkrId(ckey);
0080     if (trkrid == TrkrDefs::tpcId)
0081     {
0082       tpc_cluskey_vec.push_back(ckey);
0083       tpc_global_vec.push_back(global);
0084     }
0085     else
0086     {
0087       // si clusters stay where they are
0088       global_moved.emplace_back(ckey,global);
0089     }
0090   }
0091
0092   // need at least 3 clusters to fit a circle
0093   if (tpc_global_vec.size() < 3)
0094   {
0095     if (_verbosity > 0)
0096     {
0097       std::cout << "  -- skip this tpc track, not enough clusters: " << tpc_global_vec.size() << std::endl;
0098     }
0099     return global_in;
0100   }
0101
0102   // fit a circle to the TPC clusters
0103   const auto [R, X0, Y0] = TrackFitUtils::circle_fit_by_taubin(tpc_global_vec);
0104
0105   // get the straight line representing the z trajectory in the form of z vs radius
0106   const auto [A, B] = TrackFitUtils::line_fit(tpc_global_vec);
0107
0108   // Now we need to move each TPC cluster associated with this track to the readout layer radius
0109   for (unsigned int i = 0; i < tpc_global_vec.size(); ++i)
0110   {
0111     TrkrDefs::cluskey cluskey = tpc_cluskey_vec[i];
0112     unsigned int layer = TrkrDefs::getLayer(cluskey);
0113     Acts::Vector3 global = tpc_global_vec[i];
0114
0115     // get circle position at target surface radius
0116     double target_radius = layer_radius[layer - 7];
0117     int ret = get_circle_circle_intersection(target_radius, R, X0, Y0, global[0], global[1], _x_proj, _y_proj);
0118     if (ret == Fun4AllReturnCodes::ABORTEVENT)
0119     {
0120       continue;  // skip to next cluster
0121     }
0122     // z projection is unique
0123     _z_proj = B + A * target_radius;
0124
0125     // get circle position at cluster radius
0126     double cluster_radius = sqrt(global[0] * global[0] + global[1] * global[1]);
0127     ret = get_circle_circle_intersection(cluster_radius, R, X0, Y0, global[0], global[1], _x_start, _y_start);
0128     if (ret == Fun4AllReturnCodes::ABORTEVENT)
0129     {
0130       continue;  // skip to next cluster
0131     }
0132     // z projection is unique
0133     _z_start = B + A * cluster_radius;
0134
0135     // calculate dx, dy, dz along circle trajectory from cluster radius to surface radius
0136     double xnew = global[0] - (_x_start - _x_proj);
0137     double ynew = global[1] - (_y_start - _y_proj);
0138     double znew = global[2] - (_z_start - _z_proj);
0139
0140     // now move the cluster to the surface radius
0141     // we keep the cluster key fixed, change the surface if necessary
0142
0143     Acts::Vector3 global_new(xnew, ynew, znew);
0144
0145     // add the new position and surface to the return object
0146     global_moved.emplace_back(cluskey, global_new);
0147
0148     if (_verbosity > 2)
0149     {
0150       std::cout << "Cluster " << cluskey << " xstart " << _x_start << " xproj " << _x_proj << " ystart " << _y_start << " yproj " << _y_proj
0151                 << " zstart " << _z_start << " zproj " << _z_proj << std::endl;
0152       std::cout << " layer " << layer << " layer radius " << target_radius << " cluster radius " << cluster_radius << std::endl;
0153       std::cout << "  global in " << global[0] << "  " << global[1] << "  " << global[2] << std::endl;
0154       std::cout << "  global new " << global_new[0] << "  " << global_new[1] << "  " << global_new[2] << std::endl;
0155     }
0156   }
0157
0158   return global_moved;
0159 }
0160
0161 int TpcClusterMover::get_circle_circle_intersection(double target_radius, double R, double X0, double Y0, double xclus, double yclus, double &x, double &y)
0162 {
0163   // finds the intersection of the fitted circle with the cylinder having radius = target_radius
0164   const auto [xplus, yplus, xminus, yminus] = TrackFitUtils::circle_circle_intersection(target_radius, R, X0, Y0);
0165
0166   // We only need to check xplus for failure, skip this TPC cluster in that case
0167   if (std::isnan(xplus))
0168   {
0169     {
0170       if (_verbosity > 1)
0171       {
0172         std::cout << " circle/circle intersection calculation failed, skip this cluster" << std::endl;
0173         std::cout << " target_radius " << target_radius << " fitted R " << R << " fitted X0 " << X0 << " fitted Y0 " << Y0 << std::endl;
0174       }
0175     }
0176     return Fun4AllReturnCodes::ABORTEVENT;  // skip to next cluster
0177   }
0178
0179   // we can figure out which solution is correct based on the cluster position in the TPC
0180   if (fabs(xclus - xplus) < 5.0 && fabs(yclus - yplus) < 5.0)  // 5 cm, large and arbitrary
0181   {
0182     x = xplus;
0183     y = yplus;
0184   }
0185   else
0186   {
0187     x = xminus;
0188     y = yminus;
0189   }
0190   return Fun4AllReturnCodes::EVENT_OK;
0191 }