packages/trackreco/PHMicromegasTpcTrackMatching.cc

0001 #include "PHMicromegasTpcTrackMatching.h"
0002
0003 //#include "PHTrackPropagating.h"     // for PHTrackPropagating
0004
0005 #include <g4detectors/PHG4CylinderGeomContainer.h>
0006
0007 #include <micromegas/CylinderGeomMicromegas.h>
0008 #include <micromegas/MicromegasDefs.h>
0009
0010 /// Tracking includes
0011 #include <trackbase/ActsGeometry.h>
0012 #include <trackbase/TpcDefs.h>
0013 #include <trackbase/TrackFitUtils.h>
0014 #include <trackbase/TrkrCluster.h>  // for TrkrCluster
0015 #include <trackbase/TrkrClusterContainer.h>
0016 #include <trackbase/TrkrClusterIterationMapv1.h>
0017 #include <trackbase/TrkrClusterv3.h>  // for TrkrCluster
0018 #include <trackbase/TrkrDefs.h>       // for cluskey, getLayer, TrkrId
0019 #include <trackbase_historic/TrackSeedHelper.h>
0020 #include <trackbase_historic/SvtxTrack.h>
0021 #include <trackbase_historic/TrackSeed.h>
0022 #include <trackbase_historic/TrackSeedContainer.h>
0023
0024 #include <tpc/TpcDistortionCorrectionContainer.h>
0025
0026 #include <fun4all/Fun4AllReturnCodes.h>
0027 #include <fun4all/SubsysReco.h>  // for SubsysReco
0028
0029 #include <phool/getClass.h>
0030 #include <phool/phool.h>
0031
0032 #include <TF1.h>
0033 #include <TVector3.h>
0034
0035 #include <array>
0036 #include <cmath>     // for sqrt, std::abs, atan2, cos
0037 #include <iostream>  // for operator<<, basic_ostream
0038 #include <map>       // for map
0039 #include <set>       // for _Rb_tree_const_iterator
0040 #include <utility>   // for pair, make_pair
0041 #include <optional>  // for line-circle function
0042
0043 namespace
0044 {
0045
0046   //! convenience square method
0047   template <class T>
0048   inline constexpr T square(const T& x)
0049   {
0050     return x * x;
0051   }
0052
0053   //! get radius from x and y
0054   template <class T>
0055   inline constexpr T get_r(const T& x, const T& y)
0056   {
0057     return std::sqrt(square(x) + square(y));
0058   }
0059
0060   //! assemble list of cluster keys associated to seeds
0061   std::vector<TrkrDefs::cluskey> get_cluster_keys( const std::vector<TrackSeed*>& seeds )
0062   {
0063     std::vector<TrkrDefs::cluskey> out;
0064     for( const auto& seed: seeds )
0065     {
0066       if( seed )
0067       { std::copy( seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter( out ) ); }
0068     }
0069     return out;
0070   }
0071
0072   // To avoid confusion it is good to wrapp the angle around 2π
0073   double normalize_angle(double phi)
0074   {
0075     while (phi < 0) phi += 2 * M_PI;
0076     while (phi >= 2 * M_PI) phi -= 2 * M_PI;
0077     return phi;
0078   }
0079   bool phi_in_range(double phi, double min, double max)
0080   {
0081     phi = normalize_angle(phi);
0082     min = normalize_angle(min);
0083     max = normalize_angle(max);
0084     if (min < max)
0085       return (phi >= min && phi <= max);
0086     else
0087     return (phi >= min || phi <= max);  // wrapped around 2π
0088   }
0089
0090   /// calculate intersection from a line to the tile plane in 3d. return true on success
0091   /**
0092    * Plane is defined as (p - ptile).ntile = 0
0093    * Line is defined as p = p0 + v*t
0094    * We solve for t and then substitute in p to find the line plane intersection
0095   */
0096   bool line_plane_intersection(
0097     const TVector3& p0,
0098     const TVector3& v,
0099     const TVector3& ptile,
0100     const TVector3& ntile,
0101     TVector3& intersect)
0102   {
0103     double denom = ntile.Dot(v);
0104     if (std::abs(denom) < 1e-6) return false;  // line and plane are parallel
0105
0106     double t = ntile.Dot(ptile - p0) / denom;
0107     intersect = p0 + t * v;
0108     return true;
0109   }
0110
0111   /// calculate intersection from a helix to the tile plane in 3d. return true on success
0112   /**
0113    * Plane is defined as (p-ptile).ntile = 0
0114    * Helix is parameterize with p = (x(t), y(t), z(t)) = (X0 + R*cos(t), Y0 + R*sin(t), slope_rz*R(t) + intersect_rz)
0115    * We substitue p and find the root of t using the Newton Raphson method for the equation:
0116    * nx*R*cos(t) + ny*R*sin(t) + nz*slope_rz*R(t) + C = 0
0117    * Where C = nx*(X0-x0) + ny*(Y0-y0) + nz*(intersect_rz-z0)
0118    * and R(t) = sqrt((X0 + R*cos(t))^2 + (Y0 + R*sin(t))^2)
0119    * Finally, the Newton-Raphson method is used to calculate the t parameter
0120   */
0121
0122   bool helix_plane_intersection(
0123     double t_min,
0124     double t_max,
0125     double zmin,
0126     double zmax,
0127     double R,
0128     double X0,
0129     double Y0,
0130     double intersect_rz,
0131     double slope_rz,
0132     const TVector3& ptile,
0133     const TVector3& ntile,
0134     TVector3& intersect)
0135   {
0136
0137     // Number of iterations and tolerance for Newton Raphson method
0138     const int max_iter = 10;
0139     const double tol = 1e-6; // microns level precision
0140
0141     // Defines C
0142     double C = ntile.X() * (X0 - ptile.X()) + ntile.Y() * (Y0 - ptile.Y()) + ntile.Z() * (intersect_rz - ptile.Z());
0143
0144
0145     // Defines the function and the corresponding derivative to be used in the Newton Raphson method
0146     auto f = [&](double t) {
0147
0148     double xt = X0 + R * cos(t);
0149     double yt = Y0 + R * sin(t);
0150     double Rt = sqrt(xt*xt + yt*yt);
0151
0152     return ntile.X() * R * std::cos(t) +
0153            ntile.Y() * R * std::sin(t) +
0154            ntile.Z() * slope_rz * Rt  +
0155            C;
0156     };
0157
0158     auto df = [&](double t) {
0159
0160     double xt = X0 + R * cos(t);
0161     double yt = Y0 + R * sin(t);
0162     double Rt = sqrt(xt*xt + yt*yt);
0163
0164     return -ntile.X() * R * sin(t) +
0165             ntile.Y() * R * cos(t) +
0166             ntile.Z() * R * slope_rz * (Y0 * cos(t) - X0 * sin(t))/Rt ;
0167     };
0168
0169     // Start of Newton-Raphson iterations
0170     auto solve_from = [&](double t_seed, TVector3& result) -> bool
0171     {
0172
0173       double t = t_seed;
0174
0175       for (int i = 0; i < max_iter; ++i)
0176       {
0177         double ft = f(t);
0178         double dft = df(t);
0179         if (std::abs(dft) < 1e-8){
0180           return false; // avoid division by near-zero
0181         }
0182         double t_new = t - ft / dft;
0183
0184         double x = X0 + R * std::cos(t_new);
0185         double y = Y0 + R * std::sin(t_new);
0186         double Rt_n = std::sqrt(x * x + y * y);
0187         double z = slope_rz * Rt_n + intersect_rz;
0188         double phi = std::atan2(y, x);
0189
0190
0191         TVector3 candidate_intersect(x, y, z);
0192         //Tolerance in phi
0193     bool phi_ok = phi_in_range(phi, t_min - 1e-4, t_max + 1e-4);
0194         //Tolerance in z
0195         bool z_ok = (z >= zmin - 1e-4 && z <= zmax + 1e-4);
0196         bool proj_ok = (std::fabs(ntile.Dot(candidate_intersect - ptile)) <= 0.05);
0197
0198     // Passes the checks for the projection inside the tile acceptance
0199         if (std::abs(t_new - t) < tol && proj_ok && phi_ok && z_ok)
0200         {
0201       result = candidate_intersect;
0202       return true;
0203         }
0204         t = t_new;
0205       }
0206
0207       return false;
0208
0209     };
0210
0211
0212     auto wrap = [&](double t) {
0213       while (t > t_max) t -= 2 * M_PI;
0214       while (t < t_min) t += 2 * M_PI;
0215       return t;
0216     };
0217
0218     std::vector<double> t_seeds;
0219     double t_center = 0.5 * (t_min + t_max);
0220     double delta = 2.0 * M_PI / 3.0;
0221
0222     // Wrap the angle
0223     for (int i = 0; i < 3; ++i)
0224     {
0225       double t = wrap(t_center + i * delta);
0226       t_seeds.push_back(t);
0227     }
0228
0229     // Looks for the solution within the tile acceptance in three different phi seeds in the Newton-Raphson (helix_plane could have more than one solution)
0230     for (double t_seed : t_seeds)
0231     {
0232       if (solve_from(t_seed, intersect)) return true;
0233     }
0234
0235     return false;
0236
0237   }
0238
0239   bool line_line_intersection(
0240       double m, double b,
0241       double x0, double y0, double nx, double ny,
0242       double& xplus, double& yplus, double& xminus, double& yminus)
0243   {
0244     if (ny == 0)
0245     {
0246       // vertical lines are defined by ny=0 and x = x0
0247       xplus = xminus = x0;
0248
0249       // calculate y accordingly
0250       yplus = yminus = m * x0 + b;
0251     }
0252     else
0253     {
0254
0255       double denom = nx + ny*m;
0256       if(denom == 0) {
0257         return false; // lines are parallel and there is no intersection
0258       }
0259
0260       double x = (nx*x0 + ny*y0 - ny*b)/denom;
0261       double y = m*x + b;
0262       // a straight line has a unique intersection point
0263       xplus = xminus = x;
0264       yplus = yminus = y;
0265
0266     }
0267
0268     return true;
0269   }
0270
0271
0272   // streamer of TVector3
0273   [[maybe_unused]] inline std::ostream& operator<<(std::ostream& out, const TVector3& vector)
0274   {
0275     out << "( " << vector.x() << "," << vector.y() << "," << vector.z() << ")";
0276     return out;
0277   }
0278
0279 }  // namespace
0280
0281 //____________________________________________________________________________..
0282 PHMicromegasTpcTrackMatching::PHMicromegasTpcTrackMatching(const std::string& name)
0283   : SubsysReco(name)
0284 {
0285 }
0286 //____________________________________________________________________________..
0287 int PHMicromegasTpcTrackMatching::Init(PHCompositeNode* /* topNode */)
0288 {
0289   std::cout
0290             << "PHMicromegasTpcTrackMatching::Init - "
0291             << " rphi_search_win inner layer " << _rphi_search_win[0]
0292             << " z_search_win inner layer " << _z_search_win[0]
0293             << " rphi_search_win outer layer " << _rphi_search_win[1]
0294             << " z_search_win outer layer " << _z_search_win[1]
0295             << std::endl;
0296
0297   std::cout << "PHMicromegasTpcTrackMatching::Init - _use_silicon: " << _use_silicon << std::endl;
0298   std::cout << "PHMicromegasTpcTrackMatching::Init - _zero_field: " << _zero_field << std::endl;
0299   std::cout << "PHMicromegasTpcTrackMatching::Init - _min_tpc_layer: " << _min_tpc_layer << std::endl;
0300   std::cout << "PHMicromegasTpcTrackMatching::Init - _max_tpc_layer: " << _max_tpc_layer << std::endl;
0301
0302   return Fun4AllReturnCodes::EVENT_OK;
0303 }
0304
0305 //____________________________________________________________________________..
0306 int PHMicromegasTpcTrackMatching::InitRun(PHCompositeNode* topNode)
0307 {
0308   // load micromegas geometry
0309   _geomContainerMicromegas = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
0310   if (!_geomContainerMicromegas)
0311   {
0312     std::cout << PHWHERE << "Could not find CYLINDERGEOM_MICROMEGAS_FULL." << std::endl;
0313     return Fun4AllReturnCodes::ABORTRUN;
0314   }
0315
0316   // ensures there are at least two micromegas layers
0317   if (_geomContainerMicromegas->get_NLayers() != _n_mm_layers)
0318   {
0319     std::cout << PHWHERE << "Inconsistent number of micromegas layers." << std::endl;
0320     return Fun4AllReturnCodes::ABORTRUN;
0321   }
0322
0323   // get first micromegas layer
0324   _min_mm_layer = static_cast<CylinderGeomMicromegas*>(_geomContainerMicromegas->GetFirstLayerGeom())->get_layer();
0325
0326   int ret = GetNodes(topNode);
0327   if (ret != Fun4AllReturnCodes::EVENT_OK)
0328   {
0329     return ret;
0330   }
0331
0332   return ret;
0333 }
0334
0335 //____________________________________________________________________________..
0336 int PHMicromegasTpcTrackMatching::process_event(PHCompositeNode* topNode)
0337 {
0338   if (_n_iteration > 0)
0339   {
0340     _iteration_map = findNode::getClass<TrkrClusterIterationMapv1>(topNode, "CLUSTER_ITERATION_MAP");
0341     if (!_iteration_map)
0342     {
0343       std::cerr << PHWHERE << "Cluster Iteration Map missing, aborting." << std::endl;
0344       return Fun4AllReturnCodes::ABORTEVENT;
0345     }
0346   }
0347
0348   // We will add the micromegas cluster to the TPC tracks already on the node tree
0349
0350   _event++;
0351
0352   if (Verbosity() > 0)
0353   {
0354     std::cout << PHWHERE << " Event " << _event << " Seed track map size " << _svtx_seed_map->size() << std::endl;
0355   }
0356
0357   // loop over the seed tracks - these are the seeds formed from matched tpc and silicon track seeds
0358   for (unsigned int seedID = 0;
0359        seedID != _svtx_seed_map->size(); ++seedID)
0360   {
0361     auto seed = _svtx_seed_map->get(seedID);
0362     auto siID = seed->get_silicon_seed_index();
0363     auto tracklet_si = _si_track_map->get(siID);
0364
0365     short int crossing = 0;
0366     if (_pp_mode)
0367     {
0368       if (!tracklet_si)
0369       {
0370         continue;  // cannot use tracks not matched to silicon because crossing is unknown
0371       }
0372
0373       crossing = tracklet_si->get_crossing();
0374       if (crossing == SHRT_MAX)
0375       {
0376         if (Verbosity() > 0)
0377         {
0378           std::cout << " svtx seed " << seedID << " with si seed " << siID
0379                     << " crossing not defined: crossing = " << crossing << " skip this track" << std::endl;
0380         }
0381         continue;
0382       }
0383     }
0384
0385     auto tpcID = seed->get_tpc_seed_index();
0386
0387     auto tracklet_tpc = _tpc_track_map->get(tpcID);
0388     if (!tracklet_tpc)
0389     {
0390       continue;
0391     }
0392
0393     if (Verbosity() >= 1)
0394     {
0395       std::cout << std::endl
0396                 << __LINE__
0397                 << ": Processing TPC seed track: " << tpcID
0398                 << ": crossing: " << crossing
0399                 << ": nhits: " << tracklet_tpc->size_cluster_keys()
0400                 << ": Total TPC tracks: " << _tpc_track_map->size()
0401                 << ": phi: " << tracklet_tpc->get_phi()
0402                 << std::endl;
0403     }
0404
0405     // Get the outermost TPC clusters for this tracklet
0406     std::vector<Acts::Vector3> clusGlobPos;
0407     std::vector<Acts::Vector3> clusGlobPos_silicon;
0408     std::vector<Acts::Vector3> clusGlobPos_mvtx;
0409     std::vector<Acts::Vector3> clusGlobPos_intt;
0410
0411     bool has_micromegas = false;
0412
0413     // try extrapolate track to Micromegas and find corresponding tile
0414     /* this is all copied from PHMicromegasTpcTrackMatching */
0415     const auto cluster_keys = get_cluster_keys({tracklet_tpc,tracklet_si});
0416     for( const auto& cluster_key:cluster_keys )
0417     {
0418
0419       // detector id and layer
0420       const auto detid = TrkrDefs::getTrkrId(cluster_key);
0421       switch( detid )
0422       {
0423
0424         case TrkrDefs::tpcId:
0425         {
0426           // layer
0427           const unsigned int layer = TrkrDefs::getLayer(cluster_key);
0428
0429           // check layer range
0430           if( layer < _min_tpc_layer || layer >= _max_tpc_layer )
0431           { continue; }
0432
0433           // get matching
0434           const auto cluster = _cluster_map->findCluster(cluster_key);
0435           clusGlobPos.push_back( m_globalPositionWrapper.getGlobalPositionDistortionCorrected(cluster_key, cluster, crossing) );
0436           break;
0437         }
0438
0439         case TrkrDefs::mvtxId:
0440         {
0441           // get matching
0442           const auto cluster = _cluster_map->findCluster(cluster_key);
0443           const auto global_position = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(cluster_key, cluster, crossing);
0444           clusGlobPos_silicon.push_back( global_position );
0445           clusGlobPos_mvtx.push_back( global_position );
0446           break;
0447         }
0448
0449         case TrkrDefs::inttId:
0450         {
0451           // get matching
0452           const auto cluster = _cluster_map->findCluster(cluster_key);
0453           const auto global_position = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(cluster_key, cluster, crossing);
0454           clusGlobPos_silicon.push_back( global_position );
0455       clusGlobPos_intt.push_back( global_position );
0456           break;
0457         }
0458
0459         case TrkrDefs::micromegasId:
0460         {
0461           /*
0462            * micromegas clusters already associated to seed
0463            * skip
0464            */
0465           has_micromegas = true;
0466           break;
0467         }
0468
0469         default:
0470         break;
0471       }
0472
0473     }
0474
0475     if( has_micromegas )
0476     {
0477       if( Verbosity() )
0478       {
0479         std::cout
0480           << "PHMicromegasTpcTrackMatching::process_event -"
0481           << " Micromegas hits already associated to TPC seed."
0482           << " Skipping this track"
0483           << std::endl;
0484       }
0485       continue;
0486     }
0487
0488     // check number of clusters
0489     if( _use_silicon )
0490     {
0491
0492       if( clusGlobPos_mvtx.size()<3 ) { continue; }
0493  //     if( clusGlobPos_intt.size()<2 ) { continue; }
0494
0495     } else {
0496
0497       if( clusGlobPos.size()<3 ) { continue; }
0498
0499     }
0500
0501     // r,z linear fit
0502     /* applies to both field ON and field OFF configurations */
0503     const auto [slope_rz, intersect_rz] = _use_silicon ?
0504       TrackFitUtils::line_fit(clusGlobPos_mvtx):
0505       TrackFitUtils::line_fit(clusGlobPos);
0506
0507     if (Verbosity() > 10)
0508     {
0509       std::cout << " r,z fitted line has slope_rz " << slope_rz << " intersect_rz " << intersect_rz << std::endl;
0510     }
0511
0512     // x,y straight fit paramers
0513     double slope_xy = 0, intersect_xy = 0;
0514
0515     // circle fit parameters
0516     double R = 0, X0 = 0, Y0 = 0;
0517
0518     if(_zero_field) {
0519
0520       if (Verbosity() > 10)
0521       {
0522         std::cout << "zero field is ON, starting TPC clusters linear fit" << std::endl;
0523       }
0524
0525       // x,y straight fit
0526       std::tie( slope_xy, intersect_xy ) = _use_silicon ?
0527         TrackFitUtils::line_fit_xy(clusGlobPos_silicon):
0528         TrackFitUtils::line_fit_xy(clusGlobPos);
0529
0530       if (Verbosity() > 10)
0531       {
0532         std::cout << " zero field x,y fit has slope_xy " << slope_xy << " intersect_xy " << intersect_xy << std::endl;
0533       }
0534
0535     } else {
0536
0537       if(Verbosity() > 10)
0538       {
0539         std::cout << "zero field is OFF, starting TPC clusters circle fit" << std::endl;
0540       }
0541
0542       // x,y circle
0543       std::tie( R, X0, Y0 ) = _use_silicon ?
0544         TrackFitUtils::circle_fit_by_taubin(clusGlobPos_silicon):
0545         TrackFitUtils::circle_fit_by_taubin(clusGlobPos);
0546
0547       // toss tracks for which the fitted circle could not have come from the vertex
0548       if (R < 40.0)
0549       {
0550         continue;
0551       }
0552
0553     }
0554
0555     // loop over micromegas layer
0556     for (unsigned int imm = 0; imm < _n_mm_layers; ++imm)
0557     {
0558       // get micromegas geometry object
0559       const unsigned int layer = _min_mm_layer + imm;
0560       const auto layergeom = static_cast<CylinderGeomMicromegas*>(_geomContainerMicromegas->GetLayerGeom(layer));
0561       const auto layer_radius = layergeom->get_radius();
0562
0563       // get intersection to track
0564       auto [xplus, yplus, xminus, yminus] =
0565         _zero_field ?
0566         TrackFitUtils::line_circle_intersection(layer_radius, slope_xy, intersect_xy):
0567         TrackFitUtils::circle_circle_intersection(layer_radius, R, X0, Y0);
0568
0569       if (Verbosity() > 10)
0570       {
0571         std::cout << "xplus: " << xplus << " yplus " << yplus << " xminus " << xminus << " yminus " << std::endl;
0572       }
0573
0574       if (!std::isfinite(xplus))
0575        {
0576          if (Verbosity() > 10)
0577          {
0578            std::cout << PHWHERE << " circle/circle intersection calculation failed, skip this case" << std::endl;
0579            std::cout << PHWHERE << " mm_radius " << layer_radius << " fitted R " << R << " fitted X0 " << X0 << " fitted Y0 " << Y0 << std::endl;
0580          }
0581
0582          continue;
0583       }
0584       // we can figure out which solution is correct based on the last cluster position in the TPC
0585       const double last_clus_phi = _use_silicon ?
0586         std::atan2(clusGlobPos_silicon.back()(1), clusGlobPos_silicon.back()(0)):
0587         std::atan2(clusGlobPos.back()(1), clusGlobPos.back()(0));
0588       double phi_plus = std::atan2(yplus, xplus);
0589       double phi_minus = std::atan2(yminus, xminus);
0590
0591       // calculate z
0592       double r = layer_radius;
0593       double z = intersect_rz + slope_rz * r;
0594
0595       // select the angle that is the closest to last cluster
0596       // store phi, apply coarse space charge corrections in calibration mode
0597       double phi = std::abs(last_clus_phi - phi_plus) < std::abs(last_clus_phi - phi_minus) ? phi_plus : phi_minus;
0598
0599       // create cylinder intersection point in world coordinates
0600       const TVector3 world_intersection_cylindrical(r * std::cos(phi), r * std::sin(phi), z);
0601
0602       // find matching tile
0603       int tileid = layergeom->find_tile_cylindrical(world_intersection_cylindrical);
0604       if (tileid < 0)
0605       {
0606         continue;
0607       }
0608
0609       // get tile center and norm vector
0610       const auto tile_center = layergeom->get_world_from_local_coords(tileid, _tGeometry, {0, 0});
0611       const double x0 = tile_center.x();
0612       const double y0 = tile_center.y();
0613       const double z0 = tile_center.z();
0614
0615       TVector3 ptile(x0, y0, z0);
0616
0617       const auto tile_norm = layergeom->get_world_from_local_vect(tileid, _tGeometry, {0, 0, 1});
0618       const double nx = tile_norm.x();
0619       const double ny = tile_norm.y();
0620       const double nz = tile_norm.z();
0621
0622       TVector3 ntile(nx, ny, nz);
0623
0624       TVector3 intersection;
0625       double x;
0626       double y;
0627
0628       if( Verbosity() > 0 )
0629       {
0630         std::cout << "tile " << tileid << " layer " << layer <<  " nx " << nx << " ny " << ny << " nz " << nz << " x0 " << x0 << " y0 " << y0 << " z0 " << z0 << std::endl;
0631       }
0632
0633       if(_zero_field) {
0634
0635     // finds the x,y coordinates in the line fit
0636
0637     if (!line_line_intersection(slope_xy, intersect_xy, x0, y0, nx, ny, xplus, yplus, xminus, yminus))
0638         {
0639           if (Verbosity() > 10)
0640           {
0641             std::cout << PHWHERE << "line_line_intersection - failed" << std::endl;
0642           }
0643           continue;
0644         }
0645
0646     // calculates z0_line with these coordinates (this is not the real intersection in z) and asigns a p0 vector in the line (consider that xplus = xminus and yplus = yminus)
0647     double x0_line = xplus;
0648     double y0_line = yplus;
0649         double r0 = get_r(x0_line, y0_line);
0650         double z0_line = slope_rz * r0 + intersect_rz;
0651
0652         TVector3 p0(x0_line, y0_line, z0_line);
0653
0654     // calculates a unit vector in the direction of the line with the slope_xy and intersect_xy considering that dx/dx=1
0655         double dy_dx = slope_xy;
0656         double y_line = slope_xy * x0_line + intersect_xy;
0657         double r_line = std::sqrt(x0_line * x0_line + y_line * y_line);
0658         double dr_dx = (x0_line + y_line * dy_dx) / r_line;
0659         double dz_dx = slope_rz * dr_dx;
0660
0661         TVector3 v(1.0, dy_dx, dz_dx);
0662         v = v.Unit();
0663
0664         // calculates the real intersection to the tile
0665     if (!line_plane_intersection(p0, v, ptile, ntile, intersection))
0666         {
0667           if (Verbosity() > 10)
0668           {
0669             std::cout << PHWHERE << "line_plane_intersection - failed" << std::endl;
0670           }
0671           continue;
0672         }
0673
0674     x = intersection.X();
0675         y = intersection.Y();
0676         z = intersection.Z();
0677
0678     // looking for projections outside the tile
0679     const double zmin = layergeom->get_zmin();
0680         const double zmax = layergeom->get_zmax();
0681     if (z < zmin || z > zmax)
0682         {
0683           if (Verbosity() > 10)
0684           {
0685             std::cout << PHWHERE << "Intersection outside tile Z bounds: z = " << z
0686                   << ", zmin = " << zmin << ", zmax = " << zmax << std::endl;
0687           }
0688           continue; // reject this projection
0689         }
0690
0691       } else {
0692
0693     auto phi_range = layergeom->get_phi_range(tileid, _tGeometry);
0694         double t_min = phi_range.first;
0695         double t_max = phi_range.second;
0696     const double zmin = layergeom->get_zmin();
0697         const double zmax = layergeom->get_zmax();
0698
0699         // calculates the real intersection to tile
0700     if (!helix_plane_intersection(t_min, t_max, zmin, zmax, R, X0, Y0, intersect_rz, slope_rz, ptile, ntile, intersection))
0701     {
0702           if (Verbosity() > 0)
0703       {
0704         std::cout << PHWHERE << "helix_plane_intersection - failed" << std::endl;
0705       }
0706       continue;
0707     }
0708
0709     x = intersection.X();
0710         y = intersection.Y();
0711         z = intersection.Z();
0712
0713       }
0714
0715       /*
0716        * create planar intersection point in world coordinates
0717        * this is the position to be compared to the clusters
0718        */
0719       const TVector3 world_intersection_planar(x, y, z);
0720
0721       // convert to tile local reference frame, apply SC correction
0722       const auto local_intersection_planar = layergeom->get_local_from_world_coords(tileid, _tGeometry, world_intersection_planar);
0723
0724       // store segmentation type
0725       const auto segmentation_type = layergeom->get_segmentation_type();
0726
0727       // generate tilesetid and get corresponding clusters
0728       const auto tilesetid = MicromegasDefs::genHitSetKey(layer, segmentation_type, tileid);
0729       const auto mm_clusrange = _cluster_map->getClusters(tilesetid);
0730
0731       // do nothing if cluster range is empty
0732       if( mm_clusrange.first == mm_clusrange.second )
0733       { continue; }
0734
0735       // keep track of cluster with smallest distance to local intersection
0736       double drphi_min = 0;
0737       double dz_min = 0;
0738       TrkrDefs::cluskey ckey_min = 0;
0739       bool first = true;
0740       for (auto clusiter = mm_clusrange.first; clusiter != mm_clusrange.second; ++clusiter)
0741       {
0742         const auto& [ckey, cluster] = *clusiter;
0743         if (_iteration_map)
0744         {
0745           if (_iteration_map->getIteration(ckey) > 0)
0746           {
0747             continue;
0748           }
0749         }
0750
0751         // compute residuals and store
0752     /* in local tile coordinate, x is along rphi, and z is along y) */
0753         const double drphi = local_intersection_planar.x() - cluster->getLocalX();
0754         const double dz = local_intersection_planar.y() - cluster->getLocalY();
0755         switch( segmentation_type )
0756         {
0757           case MicromegasDefs::SegmentationType::SEGMENTATION_PHI:
0758           {
0759             // reject if outside of strip boundary
0760             if( std::abs(dz)>_z_search_win[imm] )
0761             { continue; }
0762
0763             // keep as best if closer to projection
0764             if( first || std::abs(drphi) < std::abs(drphi_min) )
0765             {
0766               first = false;
0767               drphi_min = drphi;
0768               dz_min = dz;
0769               ckey_min = ckey;
0770             }
0771             break;
0772           }
0773
0774           case MicromegasDefs::SegmentationType::SEGMENTATION_Z:
0775           {
0776             // reject if outside of strip boundary
0777             if( std::abs(drphi)>_rphi_search_win[imm] )
0778             { continue; }
0779
0780             // keep as best if closer to projection
0781             if( first || std::abs(dz) < std::abs(dz_min) )
0782             {
0783               first = false;
0784               drphi_min = drphi;
0785               dz_min = dz;
0786               ckey_min = ckey;
0787             }
0788             break;
0789           }
0790         }
0791
0792         // prints out a line that can be grep-ed from the output file to feed to a display macro
0793         // compare to cuts and add to track if matching
0794         if( _test_windows && std::abs(drphi) < _rphi_search_win[imm] && std::abs(dz) < _z_search_win[imm])
0795         {
0796
0797           // cluster rphi and z
0798       const auto glob = _tGeometry->getGlobalPosition(ckey, cluster);
0799           const double mm_clus_rphi = get_r(glob.x(), glob.y()) * std::atan2(glob.y(), glob.x());
0800           const double mm_clus_z = glob.z();
0801
0802           // projection phi and z, without correction
0803           const double rphi_proj = get_r(world_intersection_planar.x(), world_intersection_planar.y()) * std::atan2(world_intersection_planar.y(), world_intersection_planar.x());
0804           const double z_proj = world_intersection_planar.z();
0805
0806           /*
0807            * Note: drphi and dz might not match the difference of the rphi and z quoted values. This is because
0808            * 1/ drphi and dz are actually calculated in Tile's local reference frame, not in world coordinates
0809            * 2/ drphi also includes SC distortion correction, which the world coordinates don't
0810           */
0811       std::cout
0812             << "  Try_mms: " << (int) layer
0813             << " drphi " << drphi
0814             << " dz " << dz
0815             << " mm_clus_rphi " << mm_clus_rphi << " mm_clus_z " << mm_clus_z
0816             << " rphi_proj " << rphi_proj << " z_proj " << z_proj
0817             << " pt " << tracklet_tpc->get_pt()
0818             << " charge " << tracklet_tpc->get_charge()
0819             << std::endl;
0820         }
0821       }  // end loop over clusters
0822
0823       // compare to cuts and add to track if matching
0824       if( (!first) && ckey_min > 0 && std::abs(drphi_min) < _rphi_search_win[imm] && std::abs(dz_min) < _z_search_win[imm])
0825       {
0826         tracklet_tpc->insert_cluster_key(ckey_min);
0827         if (Verbosity() > 0)
0828         {
0829           std::cout << " Match to MM's found for seedID " << seedID << " tpcID " << tpcID << " siID " << siID << std::endl;
0830         }
0831       }
0832
0833     }  // end loop over Micromegas layers
0834
0835     if (Verbosity() > 3)
0836     {
0837       tracklet_tpc->identify();
0838     }
0839   }
0840
0841   if (Verbosity() > 0)
0842   {
0843     std::cout << " Final seed map size " << _svtx_seed_map->size() << std::endl;
0844   }
0845
0846   return Fun4AllReturnCodes::EVENT_OK;
0847 }
0848
0849 //_________________________________________________________________________________________________
0850 int PHMicromegasTpcTrackMatching::End(PHCompositeNode* /*unused*/)
0851 {
0852   return Fun4AllReturnCodes::EVENT_OK;
0853 }
0854
0855 //_________________________________________________________________________________________________
0856 int PHMicromegasTpcTrackMatching::GetNodes(PHCompositeNode* topNode)
0857 {
0858   // all clusters
0859   if (_use_truth_clusters)
0860   {
0861     _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER_TRUTH");
0862   }
0863   else
0864   {
0865     _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0866   }
0867
0868   if (!_cluster_map)
0869   {
0870     std::cerr << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
0871     return Fun4AllReturnCodes::ABORTEVENT;
0872   }
0873
0874   _tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0875   if (!_tGeometry)
0876   {
0877     std::cerr << PHWHERE << "No acts tracking geometry, can't continue."
0878               << std::endl;
0879     return Fun4AllReturnCodes::ABORTEVENT;
0880   }
0881
0882   _svtx_seed_map = findNode::getClass<TrackSeedContainer>(topNode, "SvtxTrackSeedContainer");
0883   if (!_svtx_seed_map)
0884   {
0885     std::cerr << PHWHERE << " ERROR: Can't find "
0886               << "SvtxTrackSeedContainer" << std::endl;
0887     return Fun4AllReturnCodes::ABORTEVENT;
0888   }
0889
0890   _tpc_track_map = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0891   if (!_tpc_track_map)
0892   {
0893     std::cerr << PHWHERE << " ERROR: Can't find "
0894               << "TpcTrackSeedContainer" << std::endl;
0895     return Fun4AllReturnCodes::ABORTEVENT;
0896   }
0897
0898   _si_track_map = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0899   if (!_si_track_map)
0900   {
0901     std::cerr << PHWHERE << " ERROR: Can't find "
0902               << "SiliconTrackSeedContainer" << std::endl;
0903     return Fun4AllReturnCodes::ABORTEVENT;
0904   }
0905
0906   // micromegas geometry
0907   _geomContainerMicromegas = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
0908   if (!_geomContainerMicromegas)
0909   {
0910     std::cout << PHWHERE << "Could not find CYLINDERGEOM_MICROMEGAS_FULL." << std::endl;
0911     return Fun4AllReturnCodes::ABORTEVENT;
0912   }
0913
0914   // global position wrapper
0915   m_globalPositionWrapper.loadNodes(topNode);
0916
0917   return Fun4AllReturnCodes::EVENT_OK;
0918 }
0919
0920 std::vector<TrkrDefs::cluskey> PHMicromegasTpcTrackMatching::getTrackletClusterList(TrackSeed* tracklet)
0921 {
0922   std::vector<TrkrDefs::cluskey> cluskey_vec;
0923   for (auto clusIter = tracklet->begin_cluster_keys();
0924        clusIter != tracklet->end_cluster_keys();
0925        ++clusIter)
0926   {
0927     auto key = *clusIter;
0928     auto cluster = _cluster_map->findCluster(key);
0929     if (!cluster)
0930     {
0931       if(Verbosity() > 1)
0932       {
0933         std::cout << PHWHERE << "Failed to get cluster with key " << key << std::endl;
0934       }
0935       continue;
0936     }
0937
0938     /// Make a safety check for clusters that couldn't be attached to a surface
0939     auto surf = _tGeometry->maps().getSurface(key, cluster);
0940     if (!surf)
0941     {
0942       continue;
0943     }
0944
0945     // drop some bad layers in the TPC completely
0946     unsigned int layer = TrkrDefs::getLayer(key);
0947     if (layer == 7 || layer == 22 || layer == 23 || layer == 38 || layer == 39)
0948     {
0949       continue;
0950     }
0951
0952     /* if (layer > 2 && layer < 7) */
0953     /* { */
0954       /* continue; */
0955     /* } */
0956
0957
0958
0959     cluskey_vec.push_back(key);
0960   }  // end loop over clusters for this track
0961   return cluskey_vec;
0962 }