packages/MvtxPrototype2/MvtxStandaloneTracking.cc

0001 #include "MvtxStandaloneTracking.h"
0002
0003 #include <trackbase/TrkrClusterContainer.h>
0004 #include <trackbase/TrkrClusterv1.h>
0005 #include <trackbase/TrkrDefs.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, "TRKR_CLUSTER");
0049   if (!clusters_)
0050   {
0051     std::cout << PHWHERE << "ERROR: Can't find node TRKR_CLUSTER" << 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::mvtxId, 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::mvtxId, 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(), mzy);
0121
0122       // -- loop over all clusters in the third desired layer
0123       TrkrClusterContainer::ConstRange clusrange2 =
0124         clusters_->getClusters(TrkrDefs::TrkrId::mvtxId, 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             bool found_clu3 = false;
0145             TrkrClusterContainer::ConstRange clusrange3 =
0146               clusters_->getClusters(TrkrDefs::TrkrId::mvtxId, lyrs.at(3));
0147             for ( TrkrClusterContainer::ConstIterator iter3 = clusrange3.first;
0148                   iter3 != clusrange3.second;
0149                   ++iter3)
0150             {
0151
0152               // check that the projection is within our window
0153               if ( fabs((iter3->second)->getX() - CalcProjection((iter3->second)->getY(), mxy, bxy)) < window_x_ &&
0154                    fabs((iter3->second)->getZ() - CalcProjection((iter3->second)->getY(), mzy, bzy)) < window_z_ )
0155               {
0156                 found_clu3 = true;
0157                 MvtxTrack tmp_trk = trk;
0158                 (tmp_trk.ClusterList).push_back(iter3->second);
0159                 trklst.push_back(tmp_trk);
0160               }
0161             }
0162             if ( !found_clu3 ) trklst.push_back(trk);
0163           }
0164           else
0165           {
0166             // else we're done
0167             trklst.push_back(trk);
0168           }
0169         }
0170       } // clusrange 2
0171     } // clusrange 1
0172   } // clusrange 0
0173
0174 }
0175
0176
0177 void
0178 MvtxStandaloneTracking::RunGhostRejection(MvtxTrackList &trklst)
0179 {
0180   if ( verbosity_ > 0 )
0181     std::cout << PHWHERE << " Running Ghost Rejection on "
0182               << trklst.size() << " tracks" << std::endl;
0183
0184   // --- First, make a map of all cluster keys & track index
0185   std::multimap<TrkrDefs::cluskey, unsigned int> key_trk_map;
0186
0187   for (unsigned int itrk = 0; itrk < trklst.size(); itrk++)
0188   {
0189     for (unsigned int iclus = 0; iclus < trklst.at(itrk).ClusterList.size(); iclus++)
0190     {
0191       TrkrDefs::cluskey ckey = trklst.at(itrk).ClusterList.at(iclus)->getClusKey();
0192       key_trk_map.insert(std::make_pair(ckey, itrk));
0193     } // iclus
0194   } // itrk
0195
0196   // --- find clusters associated with more than one track and pick the best track
0197   std::set<unsigned int> remove_set;
0198   for ( auto iter = key_trk_map.begin(); iter != key_trk_map.end(); ++iter)
0199   {
0200     int ntrk = key_trk_map.count(iter->first);
0201     if ( ntrk > 1 ) //cluster belong to more than one track
0202     {
0203       if ( verbosity_ > 1 )
0204       {
0205         std::cout << PHWHERE << " Tracks sharing cluster:" << ntrk << std::endl;
0206       }
0207
0208       // get the upper bound for this key
0209       auto upiter = key_trk_map.upper_bound(iter->first);
0210
0211       // iterate over common clusters and get the best track
0212       double chi2_best = 9999.;
0213       unsigned int idx_best = 0;
0214       int ntrk = 0;
0215       for ( auto jter = iter; jter != upiter; ++jter)
0216       {
0217         ntrk++;
0218         double chi2 = trklst.at(jter->second).chi2_xy + trklst.at(jter->second).chi2_zy;
0219         if ( chi2 < chi2_best && chi2 > 0 )
0220         {
0221           chi2_best = chi2;
0222           idx_best = jter->second;
0223         }
0224       }
0225
0226       for ( auto jter = iter; jter != upiter; ++jter)
0227       {
0228         if ( verbosity_ > 1 )
0229         {
0230           double chi2 = trklst.at(jter->second).chi2_xy + trklst.at(jter->second).chi2_zy;
0231           std::cout << "     "
0232                     << " trk idx: " << jter->second
0233                     << " chi2:" << chi2
0234                     << " m_xy:" << trklst.at(jter->second).m_xy;
0235
0236           if ( jter->second == idx_best)
0237           {
0238             std::cout << "  <--- BEST" << std::endl;
0239           }
0240           else
0241           {
0242             std::cout << std::endl;
0243           }
0244         }
0245         // remove all pairs that aren't the best
0246         if ( jter->second != idx_best )
0247         {
0248           remove_set.insert(jter->second);
0249         }
0250       }
0251     }
0252
0253   } // iter
0254
0255
0256   // --- Now remove tracks from the track list
0257   // reverse iterate so we don't have to keep track of indeces changed by removal
0258   if ( verbosity_ > 0 )
0259   {
0260     std::cout << PHWHERE << " List of idx to remove:";
0261     for ( auto it = remove_set.begin(); it != remove_set.end(); ++it)
0262     {
0263       std::cout << " " << *it;
0264     }
0265     std::cout << std::endl;
0266   }
0267
0268
0269   for ( auto rit = remove_set.rbegin(); rit != remove_set.rend(); ++rit)
0270     trklst.erase(trklst.begin() + *rit);
0271
0272   // --- done
0273
0274   if ( verbosity_ > 1 )
0275   {
0276     std::cout << PHWHERE << " Remaining tracks:" << std::endl;
0277     for ( unsigned int i = 0; i < trklst.size(); i++)
0278     {
0279       double chi2 = trklst.at(i).chi2_xy + trklst.at(i).chi2_zy;
0280       std::cout << "    chi2:" << chi2 << " m_xy:" << trklst.at(i).m_xy << std::endl;
0281     }
0282   }
0283
0284   return;
0285 }
0286
0287
0288 TVectorD
0289 MvtxStandaloneTracking::SolveGLS(TMatrixD &X, TVectorD &y, TMatrixD &L)
0290 {
0291   // Simple generalized linear least-squares fitter.
0292   // Solve y(X) = beta'*X + error(L) for beta (vector of regression coefs.)
0293   // L is the inverse covariance matrix for y.
0294   // Least-squares solution involves solving X' L X beta = X' L y
0295
0296   TMatrixD XT(X); XT.T();
0297   TMatrixD A = XT * L * X;
0298   TVectorD b = XT * L * y;
0299
0300   // Now solve A*beta = b using SVD. Decompose A = U S V'.
0301   TDecompSVD svd(A);
0302   TMatrixD UT = svd.GetU(); UT.T();
0303
0304   // Construct Moore-Penrose pseudoinverse of S
0305   TVectorD s = svd.GetSig();
0306   TMatrixD Sd(s.GetNrows(), s.GetNrows());
0307   for (int i = 0; i < s.GetNrows(); i++)
0308     Sd(i, i) = s(i) > 0 ? 1. / s(i) : 0.;
0309
0310   TVectorD beta = svd.GetV() * Sd * UT * b;
0311
0312   return beta;
0313 }
0314
0315 void
0316 MvtxStandaloneTracking::TrackFitXY(MvtxTrack &trk)
0317 {
0318   // Longitudinal/polar component
0319   // Fit track using a straight line x(y') = x0 + c*y'.
0320   // Assigns residuals, z-intercept, and polar angle.
0321
0322   // m = # measurements; n = # parameters.
0323   int m = (trk.ClusterList).size(), n = 2;
0324
0325   TMatrixD X(m, n);
0326   TMatrixD Cinv(m, m);
0327   TVectorD y(m);
0328   TVectorD x(m);
0329
0330   for (int iclus = 0; iclus < m; iclus++)
0331   {
0332     TrkrCluster *clus = (trk.ClusterList).at(iclus);
0333
0334     x(iclus) = clus->getX();
0335     y(iclus) = clus->getY();
0336
0337     X(iclus, 0) = 1;
0338     X(iclus, 1) = y(iclus);
0339
0340     Cinv(iclus, iclus) = 2.0 * sqrt(clus->getSize(0, 0));
0341   }
0342
0343   TVectorD beta = SolveGLS(X, x, Cinv);
0344   double x0 = beta(0), c = beta(1);
0345
0346   trk.m_xy = c;
0347   trk.b_xy = x0;
0348
0349   double chi2  = 0;
0350   for (int iclus = 0; iclus < m; iclus++)
0351   {
0352     TrkrCluster *clus = (trk.ClusterList).at(iclus);
0353
0354     double cx = clus->getX();
0355     double cy = clus->getY();
0356
0357     double px = cy * c + x0;
0358
0359     double dx = cx - px;
0360     trk.dx.push_back(dx);
0361
0362     chi2 += pow(dx, 2) / clus->getError(0, 0);
0363   }
0364   chi2 /= double(m - 2);
0365   trk.chi2_xy = chi2;
0366
0367   return;
0368 }
0369
0370 void
0371 MvtxStandaloneTracking::TrackFitZY(MvtxTrack &trk)
0372 {
0373   // Longitudinal/polar component
0374   // Fit track using a straight line z(y') = z0 + c*y'.
0375   // Assigns residuals, z-intercept, and polar angle.
0376
0377   // m = # measurements; n = # parameters.
0378   int m = (trk.ClusterList).size(), n = 2;
0379
0380   TMatrixD X(m, n);
0381   TMatrixD Cinv(m, m);
0382   TVectorD y(m);
0383   TVectorD z(m);
0384
0385   for (int iclus = 0; iclus < m; iclus++)
0386   {
0387     TrkrCluster *clus = (trk.ClusterList).at(iclus);
0388
0389     z(iclus) = clus->getZ();
0390     y(iclus) = clus->getY();
0391
0392     X(iclus, 0) = 1;
0393     X(iclus, 1) = y(iclus);
0394
0395     Cinv(iclus, iclus) = 2.0 * sqrt(clus->getSize(2, 2));
0396   }
0397
0398   TVectorD beta = SolveGLS(X, z, Cinv);
0399   double z0 = beta(0), c = beta(1);
0400
0401   trk.m_zy = c;
0402   trk.b_zy = z0;
0403
0404   double chi2 = 0;
0405   for (int iclus = 0; iclus < m; iclus++)
0406   {
0407     TrkrCluster *clus = (trk.ClusterList).at(iclus);
0408
0409     double cz = clus->getZ();
0410     double cy = clus->getY();
0411
0412     double pz = cy * c + z0;
0413
0414     double dz = cz - pz;
0415     trk.dz.push_back(dz);
0416
0417     chi2 += pow(dz, 2) / clus->getError(2, 2);
0418   }
0419   chi2 /= double(m - 2);
0420   trk.chi2_zy = chi2;
0421
0422   return;
0423 }
0424
0425 void
0426 MvtxStandaloneTracking::PrintTrackCandidates(MvtxTrackList &trklst)
0427 {
0428   std::cout << "===================================================" << std::endl;
0429   std::cout << "== " << PHWHERE << " Found " << trklst.size() << " Track Candidates" << std::endl;
0430   std::cout << "===================================================" << std::endl;
0431   for ( unsigned int i = 0; i < trklst.size(); i++)
0432   {
0433     std::cout << "== " << i << std::endl;
0434     for ( unsigned int j = 0; j < trklst.at(i).ClusterList.size(); j++)
0435     {
0436       std::cout << "    clus " << j
0437                 << " key:0x" << std::hex << trklst.at(i).ClusterList.at(j)->getClusKey() << std::dec
0438                 << " (" << trklst.at(i).ClusterList.at(j)->getX()
0439                 << ", " << trklst.at(i).ClusterList.at(j)->getY()
0440                 << ", " << trklst.at(i).ClusterList.at(j)->getZ()
0441                 << ")"
0442                 << " dx:" << trklst.at(i).dx.at(j)
0443                 << " dz:" << trklst.at(i).dz.at(j)
0444                 << std::endl;
0445     }
0446     std::cout << "    xy"
0447               << " m:" << trklst.at(i).m_xy
0448               << " b:" << trklst.at(i).b_xy
0449               << " chi2:" << trklst.at(i).chi2_xy
0450               << std::endl;
0451     std::cout << "    zy"
0452               << " m:" << trklst.at(i).m_zy
0453               << " b:" << trklst.at(i).b_zy
0454               << " chi2:" << trklst.at(i).chi2_zy
0455               << std::endl;
0456   } // i
0457   std::cout << "===================================================" << std::endl;
0458 }
0459
0460 double
0461 MvtxStandaloneTracking::CalcSlope(double x0, double y0, double x1, double y1)
0462 {
0463   return (y1 - y0) / (x1 - x0);
0464 }
0465
0466 double
0467 MvtxStandaloneTracking::CalcIntecept(double x0, double y0, double m)
0468 {
0469   return y0 - x0 * m;
0470 }
0471
0472 double
0473 MvtxStandaloneTracking::CalcProjection(double x, double m, double b)
0474 {
0475   return m * x + b;
0476 }
0477