packages/trackreco/PHTpcClusterMover.cc

0001 #include "PHTpcClusterMover.h"
0002
0003 /// Tracking includes
0004
0005 #include <trackbase/TrackFitUtils.h>
0006 #include <trackbase/TrkrClusterContainerv4.h>
0007 #include <trackbase/TrkrClusterv3.h>  // for TrkrCluster
0008 #include <trackbase/TrkrDefs.h>       // for cluskey, getLayer, TrkrId
0009 #include <trackbase_historic/ActsTransformations.h>
0010 #include <trackbase_historic/SvtxTrack.h>  // for SvtxTrack, SvtxTrack::C...
0011 #include <trackbase_historic/SvtxTrackMap.h>
0012
0013 #include <g4detectors/PHG4TpcCylinderGeom.h>
0014 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0015
0016 #include <fun4all/Fun4AllReturnCodes.h>
0017
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/getClass.h>
0020 #include <phool/phool.h>
0021
0022 #include <TF1.h>
0023
0024 #include <cmath>     // for sqrt, fabs, atan2, cos
0025 #include <iostream>  // for operator<<, basic_ostream
0026 #include <map>       // for map
0027 #include <set>       // for _Rb_tree_const_iterator
0028 #include <utility>   // for pair, make_pair
0029
0030 //____________________________________________________________________________..
0031 PHTpcClusterMover::PHTpcClusterMover(const std::string &name)
0032   : SubsysReco(name)
0033 {
0034 }
0035
0036 //____________________________________________________________________________..
0037 int PHTpcClusterMover::InitRun(PHCompositeNode *topNode)
0038 {
0039   int ret = GetNodes(topNode);
0040   if (ret != Fun4AllReturnCodes::EVENT_OK)
0041   {
0042     return ret;
0043   }
0044   for (int layer = 7; layer < 7 + 48; layer++)
0045   {
0046     PHG4TpcCylinderGeom *GeoLayer = _tpc_geom_container->GetLayerCellGeom(layer);
0047     std::cout << "PHTpcClusterMover:: layer = " << layer << " layer_radius " << GeoLayer->get_radius() << std::endl;
0048     layer_radius[layer - 7] = GeoLayer->get_radius();
0049   }
0050
0051   return ret;
0052 }
0053
0054 //____________________________________________________________________________..
0055 int PHTpcClusterMover::process_event(PHCompositeNode * /*topNode*/)
0056 {
0057   if (Verbosity() > 0)
0058   {
0059     std::cout << PHWHERE << " track map size " << _track_map->size() << std::endl;
0060   }
0061
0062   // loop over the tracks
0063   for (auto &phtrk_iter : *_track_map)
0064   {
0065     _track = phtrk_iter.second;
0066
0067     const auto crossing = _track->get_crossing();
0068     assert(crossing != SHRT_MAX);
0069
0070     if (Verbosity() >= 1)
0071     {
0072       std::cout << std::endl
0073                 << __LINE__
0074                 << ": Processing track itrack: " << phtrk_iter.first
0075                 << ": nhits: " << _track->size_cluster_keys()
0076                 << ": Total tracks: " << _track_map->size()
0077                 << ": phi: " << _track->get_phi()
0078                 << std::endl;
0079     }
0080
0081     // Get the TPC clusters for this track and correct them for distortions
0082     std::vector<Acts::Vector3> globalClusterPositions;
0083     std::map<TrkrDefs::cluskey, Acts::Vector3> tpc_clusters;
0084
0085     for (auto key_iter = _track->begin_cluster_keys(); key_iter != _track->end_cluster_keys(); ++key_iter)
0086     {
0087       const TrkrDefs::cluskey &cluster_key = *key_iter;
0088       const unsigned int trkrId = TrkrDefs::getTrkrId(cluster_key);
0089
0090       // non Tpc clusters are copied unchanged to the new map if they are present
0091       // This is needed for truth seeding case, where the tracks already have silicon clusters
0092       if (trkrId != TrkrDefs::tpcId)
0093       {
0094         // check if clusters has not been inserted already
0095         if (_corrected_cluster_map->findCluster(cluster_key))
0096         {
0097           continue;
0098         }
0099
0100         // get cluster from original map
0101         auto cluster = _cluster_map->findCluster(cluster_key);
0102         if (!cluster)
0103         {
0104           continue;
0105         }
0106
0107         // create a copy
0108         auto newclus = new TrkrClusterv3;
0109         newclus->CopyFrom(cluster);
0110
0111         // insert in corrected map
0112         _corrected_cluster_map->addClusterSpecifyKey(cluster_key, newclus);
0113         continue;
0114       }
0115
0116       // get the cluster in 3D coordinates
0117       auto tpc_clus = _cluster_map->findCluster(cluster_key);
0118       const auto global = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(cluster_key, tpc_clus, crossing);
0119
0120       // Store the corrected 3D cluster positions
0121       globalClusterPositions.push_back(global);
0122       tpc_clusters.insert(std::make_pair(cluster_key, global));
0123     }
0124
0125     // need at least 3 clusters to fit a circle
0126     if (globalClusterPositions.size() < 3)
0127     {
0128       if (Verbosity() > 3)
0129       {
0130         std::cout << PHWHERE << "  -- skip this tpc track, not enough clusters: " << globalClusterPositions.size() << std::endl;
0131       }
0132       continue;  // skip to the next TPC track
0133     }
0134
0135     // fit a circle to the clusters
0136     const auto [R, X0, Y0] = TrackFitUtils::circle_fit_by_taubin(globalClusterPositions);
0137     if (Verbosity() > 10)
0138     {
0139       std::cout << " Fitted circle has R " << R << " X0 " << X0 << " Y0 " << Y0 << std::endl;
0140     }
0141
0142     // toss tracks for which the fitted circle could not have come from the vertex
0143     // if(R < 30.0) continue;
0144
0145     // get the straight line representing the z trajectory in the form of z vs radius
0146     const auto [A, B] = TrackFitUtils::line_fit(globalClusterPositions);
0147     if (Verbosity() > 10)
0148     {
0149       std::cout << " Fitted line has A " << A << " B " << B << std::endl;
0150     }
0151
0152     // Now we need to move each cluster associated with this track to the readout layer radius
0153     for (const auto &[cluskey, global] : tpc_clusters)
0154     {
0155       const unsigned int layer = TrkrDefs::getLayer(cluskey);
0156
0157       // get circle position at target surface radius
0158       double target_radius = layer_radius[layer - 7];
0159       int ret = get_circle_circle_intersection(target_radius, R, X0, Y0, global[0], global[1], _x_proj, _y_proj);
0160       if (ret == Fun4AllReturnCodes::ABORTEVENT)
0161       {
0162         continue;  // skip to next cluster
0163       }
0164       // z projection is unique
0165       _z_proj = B + A * target_radius;
0166
0167       // get circle position at cluster radius
0168       double cluster_radius = sqrt(global[0] * global[0] + global[1] * global[1]);
0169       ret = get_circle_circle_intersection(cluster_radius, R, X0, Y0, global[0], global[1], _x_start, _y_start);
0170       if (ret == Fun4AllReturnCodes::ABORTEVENT)
0171       {
0172         continue;  // skip to next cluster
0173       }
0174       // z projection is unique
0175       _z_start = B + A * cluster_radius;
0176
0177       // calculate dx, dy, dz along circle trajectory from cluster radius to surface radius
0178       double xnew = global[0] - (_x_start - _x_proj);
0179       double ynew = global[1] - (_y_start - _y_proj);
0180       double znew = global[2] - (_z_start - _z_proj);
0181
0182       // now move the cluster to the surface radius
0183       // we keep the cluster key fixed, change the surface if necessary
0184       // write the new cluster position local coordinates on the surface
0185
0186       Acts::Vector3 global_new(xnew, ynew, znew);
0187
0188       TrkrDefs::subsurfkey subsurfkey;
0189       TrkrDefs::hitsetkey tpcHitSetKey = TrkrDefs::getHitSetKeyFromClusKey(cluskey);
0190       Surface surface = _tGeometry->get_tpc_surface_from_coords(
0191           tpcHitSetKey,
0192           global_new,
0193           subsurfkey);
0194
0195       if (!surface)
0196       {
0197         /// If the surface can't be found, we can't track with it. So
0198         /// just continue and don't modify the cluster to the container
0199         std::cout << PHWHERE << "Failed to find surface for cluster " << cluskey << std::endl;
0200         continue;
0201       }
0202
0203       // get the original cluster
0204       TrkrCluster *cluster = _cluster_map->findCluster(cluskey);
0205
0206       // put the corrected cluster in the new cluster map
0207
0208       // ghost tracks can have repeat clusters, so check if cluster already moved
0209       if (_corrected_cluster_map->findCluster(cluskey))
0210       {
0211         continue;
0212       }
0213
0214       // create new cluster
0215       auto newclus = new TrkrClusterv3;
0216
0217       // copy from source
0218       newclus->CopyFrom(cluster);
0219
0220       // assign subsurface key
0221       newclus->setSubSurfKey(subsurfkey);
0222
0223       // get local coordinates
0224       Acts::Vector3 normal = surface->normal(_tGeometry->geometry().getGeoContext(), Acts::Vector3(1,1,1),Acts::Vector3(1,1,1));
0225       auto local = surface->globalToLocal(_tGeometry->geometry().getGeoContext(),
0226                                           global * Acts::UnitConstants::cm,
0227                                           normal);
0228
0229       Acts::Vector2 localPos;
0230       if (local.ok())
0231       {
0232         localPos = local.value() / Acts::UnitConstants::cm;
0233       }
0234       else
0235       {
0236         /// otherwise take the manual calculation
0237         Acts::Vector3 center = surface->center(_tGeometry->geometry().getGeoContext()) / Acts::UnitConstants::cm;
0238         double clusRadius = sqrt(xnew * xnew + ynew * ynew);
0239         double clusphi = atan2(ynew, xnew);
0240         double rClusPhi = clusRadius * clusphi;
0241         double surfRadius = sqrt(center(0) * center(0) + center(1) * center(1));
0242         double surfPhiCenter = atan2(center[1], center[0]);
0243         double surfRphiCenter = surfPhiCenter * surfRadius;
0244         double surfZCenter = center[2];
0245
0246         localPos(0) = rClusPhi - surfRphiCenter;
0247         localPos(1) = znew - surfZCenter;
0248       }
0249
0250       if (Verbosity() > 4)
0251       {
0252         std::cout << "*** cluster_radius " << cluster_radius << " cluster x,y,z: " << global[0] << "  " << global[1] << "  " << global[2] << std::endl;
0253         std::cout << "    projection_radius " << target_radius << " proj x,y,z: " << _x_proj << "  " << _y_proj << "  " << _z_proj << std::endl;
0254         std::cout << "    traj_start_radius " << cluster_radius << " start x,y,z: " << _x_start << "  " << _y_start << "  " << _z_start << std::endl;
0255         std::cout << "    moved_clus_radius " << target_radius << " final x,y,z: " << xnew << "  " << ynew << "  " << znew << std::endl;
0256       }
0257
0258       // assign to new cluster
0259       newclus->setLocalX(localPos(0));
0260       newclus->setLocalY(localPos(1));
0261
0262       // insert in map
0263       _corrected_cluster_map->addClusterSpecifyKey(cluskey, newclus);
0264     }
0265
0266     // For normal reconstruction, the silicon clusters  for this track will be copied over after the matching is done
0267   }
0268
0269   return Fun4AllReturnCodes::EVENT_OK;
0270 }
0271
0272 int PHTpcClusterMover::End(PHCompositeNode * /*topNode*/)
0273 {
0274   return Fun4AllReturnCodes::EVENT_OK;
0275 }
0276
0277 int PHTpcClusterMover::GetNodes(PHCompositeNode *topNode)
0278 {
0279   // tpc geometry
0280   _tpc_geom_container = findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0281   if (!_tpc_geom_container)
0282   {
0283     std::cout << PHWHERE << " ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0284     return Fun4AllReturnCodes::ABORTEVENT;
0285   }
0286
0287   // clusters
0288   _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0289   if (!_cluster_map)
0290   {
0291     std::cout << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
0292     return Fun4AllReturnCodes::ABORTEVENT;
0293   }
0294
0295   // tracks
0296   _track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0297   if (!_track_map)
0298   {
0299     std::cout << PHWHERE << " ERROR: Can't find SvtxTrackMap: " << std::endl;
0300     return Fun4AllReturnCodes::ABORTEVENT;
0301   }
0302
0303   // geometry
0304   _tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0305   if (!_tGeometry)
0306   {
0307     std::cout << PHWHERE << "Error, can't find acts tracking geometry" << std::endl;
0308     return Fun4AllReturnCodes::ABORTEVENT;
0309   }
0310
0311   // tpc global position wrapper
0312   m_globalPositionWrapper.loadNodes(topNode);
0313
0314   // create the node for distortion corrected clusters, if it does not already exist
0315   _corrected_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
0316   if (!_corrected_cluster_map)
0317   {
0318     std::cout << "Creating node CORRECTED_TRKR_CLUSTER" << std::endl;
0319
0320     PHNodeIterator iter(topNode);
0321
0322     // Looking for the DST node
0323     PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0324     if (!dstNode)
0325     {
0326       std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0327       return Fun4AllReturnCodes::ABORTRUN;
0328     }
0329     PHNodeIterator dstiter(dstNode);
0330     PHCompositeNode *DetNode =
0331         dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0332     if (!DetNode)
0333     {
0334       DetNode = new PHCompositeNode("TRKR");
0335       dstNode->addNode(DetNode);
0336     }
0337
0338     _corrected_cluster_map = new TrkrClusterContainerv4;
0339     PHIODataNode<PHObject> *TrkrClusterContainerNode =
0340         new PHIODataNode<PHObject>(_corrected_cluster_map, "CORRECTED_TRKR_CLUSTER", "PHObject");
0341     DetNode->addNode(TrkrClusterContainerNode);
0342   }
0343   else
0344   {
0345     _corrected_cluster_map->Reset();
0346   }
0347
0348   return Fun4AllReturnCodes::EVENT_OK;
0349 }
0350
0351 int PHTpcClusterMover::get_circle_circle_intersection(double target_radius, double R, double X0, double Y0, double xclus, double yclus, double &x, double &y)
0352 {
0353   // finds the intersection of the fitted circle with the cylinder having radius = target_radius
0354   const auto [xplus, yplus, xminus, yminus] = TrackFitUtils::circle_circle_intersection(target_radius, R, X0, Y0);
0355
0356   // We only need to check xplus for failure, skip this TPC cluster in that case
0357   if (std::isnan(xplus))
0358   {
0359     if (Verbosity() > 0)
0360     {
0361       std::cout << " circle/circle intersection calculation failed, skip this cluster" << std::endl;
0362       std::cout << " target_radius " << target_radius << " fitted R " << R << " fitted X0 " << X0 << " fitted Y0 " << Y0 << std::endl;
0363     }
0364     return Fun4AllReturnCodes::ABORTEVENT;  // skip to next cluster
0365   }
0366
0367   // we can figure out which solution is correct based on the cluster position in the TPC
0368   if (fabs(xclus - xplus) < 5.0 && fabs(yclus - yplus) < 5.0)  // 5 cm, large and arbitrary
0369   {
0370     x = xplus;
0371     y = yplus;
0372   }
0373   else
0374   {
0375     x = xminus;
0376     y = yminus;
0377   }
0378   return Fun4AllReturnCodes::EVENT_OK;
0379 }