packages/MvtxPrototype1/MvtxStandaloneTracking.C

0001 #include "MvtxStandaloneTracking.h"
0002
0003 #include <trackbase/TrkrClusterContainer.h>
0004 #include <trackbase/TrkrClusterv1.h>
0005 #include <trackbase/TrkrDefUtil.h>
0006
0007 #include <phool/PHCompositeNode.h>
0008 #include <phool/PHIODataNode.h>
0009 #include <phool/PHNodeIterator.h>
0010 #include <phool/getClass.h>
0011
0012 #include <TDecompSVD.h>
0013 #include <TMath.h>
0014
0015 #include <iostream>
0016 #include <utility>
0017 #include <stdio.h>
0018 #include <map>
0019
0020 MvtxStandaloneTracking::MvtxStandaloneTracking()
0021   : window_x_(10)
0022   , window_z_(10)
0023   , ghostrejection_(false)
0024   , verbosity_(0)
0025 {
0026
0027 }
0028
0029 MvtxStandaloneTracking::~MvtxStandaloneTracking()
0030 {
0031 }
0032
0033 void
0034 MvtxStandaloneTracking::RunTracking(PHCompositeNode* topNode, MvtxTrackList &trklst, std::vector<int> &lyrs)
0035 {
0036   // check for appropriate layers
0037   if ( lyrs.size() < 3 || lyrs.size() > 4 )
0038   {
0039     std::cout << PHWHERE << "ERROR: Inappropriate number of input layers - " << lyrs.size() << std::endl;
0040     return;
0041   }
0042
0043   trklst.clear();
0044
0045   //------
0046   //--- get cluster container
0047   //------
0048   clusters_ = findNode::getClass<TrkrClusterContainer>(topNode, "TrkrClusterContainer");
0049   if (!clusters_)
0050   {
0051     std::cout << PHWHERE << "ERROR: Can't find node TrkrClusterContainer" << std::endl;
0052     return;
0053   }
0054
0055
0056   //------
0057   //--- associate clusters
0058   //------
0059   AssociateClusters(trklst, lyrs);
0060
0061   if ( verbosity_ > 0 )
0062     std::cout << PHWHERE << " Finished associating clusters. N candidates:" << trklst.size() << std::endl;
0063   //------
0064   //--- fit tracks
0065   //------
0066   for ( unsigned int itrk = 0; itrk < trklst.size(); itrk++)
0067   {
0068     TrackFitXY(trklst.at(itrk));
0069     TrackFitZY(trklst.at(itrk));
0070   } // itrk
0071
0072   if ( verbosity_ > 1 )
0073     PrintTrackCandidates(trklst);
0074
0075
0076   //------
0077   // --- choose best tracks
0078   //------
0079   if ( ghostrejection_ && trklst.size() > 1 )
0080     RunGhostRejection(trklst);
0081
0082   // --- done
0083   return;
0084 }
0085
0086 void
0087 MvtxStandaloneTracking::AssociateClusters(MvtxTrackList &trklst, std::vector<int> &lyrs)
0088 {
0089
0090   // --- utility class
0091   TrkrDefUtil util;
0092
0093
0094   // --- loop over all clusters in the first desired layer
0095   TrkrClusterContainer::ConstRange clusrange0 =
0096     clusters_->GetClusters(TrkrDefs::TRKRID::mvtx_id, lyrs.at(0));
0097   for ( TrkrClusterContainer::ConstIterator iter0 = clusrange0.first;
0098         iter0 != clusrange0.second;
0099         ++iter0)
0100   {
0101
0102     // -- loop over all clusters in the second desired layer
0103     TrkrClusterContainer::ConstRange clusrange1 =
0104       clusters_->GetClusters(TrkrDefs::TRKRID::mvtx_id, lyrs.at(1));
0105     for ( TrkrClusterContainer::ConstIterator iter1 = clusrange1.first;
0106           iter1 != clusrange1.second;
0107           ++iter1)
0108     {
0109
0110       // get clusters
0111       TrkrCluster* clus0 = iter0->second;
0112       TrkrCluster* clus1 = iter1->second;
0113
0114       // calculate slope & interecept in xy plane
0115       double mxy = CalcSlope(clus0->GetY(), clus0->GetX(), clus1->GetY(), clus1->GetX());
0116       double bxy = CalcIntecept(clus0->GetY(), clus0->GetX(), mxy);
0117
0118       // calculate slope & interecept in zy plane
0119       double mzy = CalcSlope(clus0->GetY(), clus0->GetZ(), clus1->GetY(), clus1->GetZ());
0120       double bzy = CalcIntecept(clus0->GetY(), clus0->GetZ(), mxy);
0121
0122       // -- loop over all clusters in the third desired layer
0123       TrkrClusterContainer::ConstRange clusrange2 =
0124         clusters_->GetClusters(TrkrDefs::TRKRID::mvtx_id, lyrs.at(2));
0125       for ( TrkrClusterContainer::ConstIterator iter2 = clusrange2.first;
0126             iter2 != clusrange2.second;
0127             ++iter2)
0128       {
0129
0130         // check that the projection is within our window
0131         if ( fabs((iter2->second)->GetX() - CalcProjection((iter2->second)->GetY(), mxy, bxy)) < window_x_ &&
0132              fabs((iter2->second)->GetZ() - CalcProjection((iter2->second)->GetY(), mzy, bzy)) < window_z_ )
0133         {
0134
0135           // make the candidate
0136           MvtxTrack trk;
0137           (trk.ClusterList).push_back(iter0->second);
0138           (trk.ClusterList).push_back(iter1->second);
0139           (trk.ClusterList).push_back(iter2->second);
0140
0141           // if there's another layer, require it
0142           if ( lyrs.size() > 3 )
0143           {
0144             TrkrClusterContainer::ConstRange clusrange3 =
0145               clusters_->GetClusters(TrkrDefs::TRKRID::mvtx_id, lyrs.at(3));
0146             for ( TrkrClusterContainer::ConstIterator iter3 = clusrange3.first;
0147                   iter3 != clusrange3.second;
0148                   ++iter3)
0149             {
0150
0151               // check that the projection is within our window
0152               if ( fabs((iter3->second)->GetX() - CalcProjection((iter3->second)->GetY(), mxy, bxy)) < window_x_ &&
0153                    fabs((iter3->second)->GetZ() - CalcProjection((iter3->second)->GetY(), mzy, bzy)) < window_z_ )
0154               {
0155
0156                 (trk.ClusterList).push_back(iter3->second);
0157                 trklst.push_back(trk);
0158               }
0159             }
0160           }
0161           // else we're done
0162           else
0163           {
0164             trklst.push_back(trk);
0165           }
0166
0167         }
0168       } // clusrange 2
0169     } // clusrange 1
0170   } // clusrange 0
0171
0172 }
0173
0174
0175 void
0176 MvtxStandaloneTracking::RunGhostRejection(MvtxTrackList &trklst)
0177 {
0178   if ( verbosity_ > 0 )
0179     std::cout << PHWHERE << " Running Ghost Rejection on "
0180               << trklst.size() << " tracks" << std::endl;
0181
0182   // --- First, make a map of all cluster keys & track index
0183   std::multimap<TrkrDefs::cluskey, unsigned int> key_trk_map;
0184
0185   for (unsigned int itrk = 0; itrk < trklst.size(); itrk++)
0186   {
0187     for (unsigned int iclus = 0; iclus < trklst.at(itrk).ClusterList.size(); iclus++)
0188     {
0189       TrkrDefs::cluskey ckey = trklst.at(itrk).ClusterList.at(iclus)->GetClusKey();
0190       key_trk_map.insert(std::make_pair(ckey, itrk));
0191     } // iclus
0192   } // itrk
0193
0194   // --- find clusters associated with more than one track and pick the best track
0195   std::set<unsigned int> remove_set;
0196   for ( auto iter = key_trk_map.begin(); iter != key_trk_map.end(); ++iter)
0197   {
0198     // get the upper bound for this key
0199     auto upiter = key_trk_map.upper_bound(iter->first);
0200
0201     // iterate over common clusters and get the best track
0202     double chi2_best = 9999.;
0203     unsigned int idx_best = 0;
0204     int ntrk = 0;
0205     for ( auto jter = iter; jter != upiter; ++jter)
0206     {
0207       ntrk++;
0208       double chi2 = trklst.at(jter->second).chi2_xy + trklst.at(jter->second).chi2_zy;
0209       if ( chi2 < chi2_best && chi2 > 0 )
0210       {
0211         chi2_best = chi2;
0212         idx_best = jter->second;
0213       }
0214     }
0215
0216     // FOR TESTING
0217     if ( ntrk > 1 && verbosity_ > 1 )
0218     {
0219       std::cout << PHWHERE << " Tracks sharing cluster:" << ntrk << std::endl;
0220       for ( auto jter = iter; jter != upiter; ++jter)
0221       {
0222         double chi2 = trklst.at(jter->second).chi2_xy + trklst.at(jter->second).chi2_zy;
0223         std::cout << "     "
0224                   << " trk idx: " << jter->second
0225                   << " chi2:" << chi2
0226                   << " m_xy:" << trklst.at(jter->second).m_xy;
0227
0228         if ( jter->second == idx_best)
0229         {
0230           std::cout << "  <--- BEST" << std::endl;
0231         }
0232         else
0233         {
0234           std::cout << std::endl;
0235         }
0236       }
0237     }
0238
0239     // remove all pairs that aren't the best
0240     for ( auto jter = iter; jter != upiter; ++jter)
0241     {
0242       if ( jter->second != idx_best )
0243       {
0244         remove_set.insert(jter->second);
0245       }
0246     }
0247   } // iter
0248
0249
0250   // --- Now remove tracks from the track list
0251   // reverse iterate so we don't have to keep track of indeces changed by removal
0252   if ( verbosity_ > 0 )
0253   {
0254     std::cout << PHWHERE << " List of idx to remove:";
0255     for ( auto it = remove_set.begin(); it != remove_set.end(); ++it)
0256     {
0257       std::cout << " " << *it;
0258     }
0259     std::cout << std::endl;
0260   }
0261
0262
0263   for ( auto rit = remove_set.rbegin(); rit != remove_set.rend(); ++rit)
0264     trklst.erase(trklst.begin() + *rit);
0265
0266   // --- done
0267
0268   if ( verbosity_ > 1 )
0269   {
0270     std::cout << PHWHERE << " Remaining tracks:" << std::endl;
0271     for ( unsigned int i = 0; i < trklst.size(); i++)
0272     {
0273       double chi2 = trklst.at(i).chi2_xy + trklst.at(i).chi2_zy;
0274       std::cout << "    chi2:" << chi2 << " m_xy:" << trklst.at(i).m_xy << std::endl;
0275     }
0276   }
0277
0278   return;
0279 }
0280
0281
0282 TVectorD
0283 MvtxStandaloneTracking::SolveGLS(TMatrixD &X, TVectorD &y, TMatrixD &L)
0284 {
0285   // Simple generalized linear least-squares fitter.
0286   // Solve y(X) = beta'*X + error(L) for beta (vector of regression coefs.)
0287   // L is the inverse covariance matrix for y.
0288   // Least-squares solution involves solving X' L X beta = X' L y
0289
0290   TMatrixD XT(X); XT.T();
0291   TMatrixD A = XT * L * X;
0292   TVectorD b = XT * L * y;
0293
0294   // Now solve A*beta = b using SVD. Decompose A = U S V'.
0295   TDecompSVD svd(A);
0296   TMatrixD UT = svd.GetU(); UT.T();
0297
0298   // Construct Moore-Penrose pseudoinverse of S
0299   TVectorD s = svd.GetSig();
0300   TMatrixD Sd(s.GetNrows(), s.GetNrows());
0301   for (int i = 0; i < s.GetNrows(); i++)
0302     Sd(i, i) = s(i) > 0 ? 1. / s(i) : 0.;
0303
0304   TVectorD beta = svd.GetV() * Sd * UT * b;
0305
0306   return beta;
0307 }
0308
0309 void
0310 MvtxStandaloneTracking::TrackFitXY(MvtxTrack &trk)
0311 {
0312   // Longitudinal/polar component
0313   // Fit track using a straight line x(y') = x0 + c*y'.
0314   // Assigns residuals, z-intercept, and polar angle.
0315
0316   // m = # measurements; n = # parameters.
0317   int m = (trk.ClusterList).size(), n = 2;
0318
0319   TMatrixD X(m, n);
0320   TMatrixD Cinv(m, m);
0321   TVectorD y(m);
0322   TVectorD x(m);
0323
0324   for (int iclus = 0; iclus < m; iclus++)
0325   {
0326     TrkrCluster *clus = (trk.ClusterList).at(iclus);
0327
0328     x(iclus) = clus->GetX();
0329     y(iclus) = clus->GetY();
0330
0331     X(iclus, 0) = 1;
0332     X(iclus, 1) = y(iclus);
0333
0334     Cinv(iclus, iclus) = 2.0 * sqrt(clus->GetSize(0, 0));
0335   }
0336
0337   TVectorD beta = SolveGLS(X, x, Cinv);
0338   double x0 = beta(0), c = beta(1);
0339
0340   trk.m_xy = c;
0341   trk.b_xy = x0;
0342
0343   double chi2  = 0;
0344   for (int iclus = 0; iclus < m; iclus++)
0345   {
0346     TrkrCluster *clus = (trk.ClusterList).at(iclus);
0347
0348     double cx = clus->GetX();
0349     double cy = clus->GetY();
0350
0351     double px = cy * c + x0;
0352
0353     double dx = cx - px;
0354     trk.dx.push_back(dx);
0355
0356     chi2 += pow(dx, 2) / clus->GetError(0, 0);
0357   }
0358   chi2 /= double(m - 2);
0359   trk.chi2_xy = chi2;
0360
0361   return;
0362 }
0363
0364 void
0365 MvtxStandaloneTracking::TrackFitZY(MvtxTrack &trk)
0366 {
0367   // Longitudinal/polar component
0368   // Fit track using a straight line z(y') = z0 + c*y'.
0369   // Assigns residuals, z-intercept, and polar angle.
0370
0371   // m = # measurements; n = # parameters.
0372   int m = (trk.ClusterList).size(), n = 2;
0373
0374   TMatrixD X(m, n);
0375   TMatrixD Cinv(m, m);
0376   TVectorD y(m);
0377   TVectorD z(m);
0378
0379   for (int iclus = 0; iclus < m; iclus++)
0380   {
0381     TrkrCluster *clus = (trk.ClusterList).at(iclus);
0382
0383     z(iclus) = clus->GetZ();
0384     y(iclus) = clus->GetY();
0385
0386     X(iclus, 0) = 1;
0387     X(iclus, 1) = y(iclus);
0388
0389     Cinv(iclus, iclus) = 2.0 * sqrt(clus->GetSize(2, 2));
0390   }
0391
0392   TVectorD beta = SolveGLS(X, z, Cinv);
0393   double z0 = beta(0), c = beta(1);
0394
0395   trk.m_zy = c;
0396   trk.b_zy = z0;
0397
0398   double chi2 = 0;
0399   for (int iclus = 0; iclus < m; iclus++)
0400   {
0401     TrkrCluster *clus = (trk.ClusterList).at(iclus);
0402
0403     double cz = clus->GetZ();
0404     double cy = clus->GetY();
0405
0406     double pz = cy * c + z0;
0407
0408     double dz = cz - pz;
0409     trk.dz.push_back(dz);
0410
0411     chi2 += pow(dz, 2) / clus->GetError(2, 2);
0412   }
0413   chi2 /= double(m - 2);
0414   trk.chi2_zy = chi2;
0415
0416   return;
0417 }
0418
0419 void
0420 MvtxStandaloneTracking::PrintTrackCandidates(MvtxTrackList &trklst)
0421 {
0422   std::cout << "===================================================" << std::endl;
0423   std::cout << "== " << PHWHERE << " Found " << trklst.size() << " Track Candidates" << std::endl;
0424   std::cout << "===================================================" << std::endl;
0425   for ( unsigned int i = 0; i < trklst.size(); i++)
0426   {
0427     std::cout << "== " << i << std::endl;
0428     for ( unsigned int j = 0; j < trklst.at(i).ClusterList.size(); j++)
0429     {
0430       std::cout << "    clus " << j
0431                 << " key:0x" << std::hex << trklst.at(i).ClusterList.at(j)->GetClusKey() << std::dec
0432                 << " (" << trklst.at(i).ClusterList.at(j)->GetX()
0433                 << ", " << trklst.at(i).ClusterList.at(j)->GetY()
0434                 << ", " << trklst.at(i).ClusterList.at(j)->GetZ()
0435                 << ")"
0436                 << " dx:" << trklst.at(i).dx.at(j)
0437                 << " dz:" << trklst.at(i).dz.at(j)
0438                 << std::endl;
0439     }
0440     std::cout << "    xy"
0441               << " m:" << trklst.at(i).m_xy
0442               << " b:" << trklst.at(i).b_xy
0443               << " chi2:" << trklst.at(i).chi2_xy
0444               << std::endl;
0445     std::cout << "    zy"
0446               << " m:" << trklst.at(i).m_zy
0447               << " b:" << trklst.at(i).b_zy
0448               << " chi2:" << trklst.at(i).chi2_zy
0449               << std::endl;
0450   } // i
0451   std::cout << "===================================================" << std::endl;
0452 }
0453
0454 double
0455 MvtxStandaloneTracking::CalcSlope(double x0, double y0, double x1, double y1)
0456 {
0457   return (y1 - y0) / (x1 - x0);
0458 }
0459
0460 double
0461 MvtxStandaloneTracking::CalcIntecept(double x0, double y0, double m)
0462 {
0463   return y0 - x0 * m;
0464 }
0465
0466 double
0467 MvtxStandaloneTracking::CalcProjection(double x, double m, double b)
0468 {
0469   return m * x + b;
0470 }
0471