12 #include <TDecompSVD.h>
23 , ghostrejection_(
false)
37 if ( lyrs.size() < 3 || lyrs.size() > 4 )
39 std::cout <<
PHWHERE <<
"ERROR: Inappropriate number of input layers - " << lyrs.size() << std::endl;
48 clusters_ = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
51 std::cout <<
PHWHERE <<
"ERROR: Can't find node TRKR_CLUSTER" << std::endl;
62 std::cout <<
PHWHERE <<
" Finished associating clusters. N candidates:" << trklst.size() << std::endl;
66 for (
unsigned int itrk = 0; itrk < trklst.size(); itrk++)
98 iter0 != clusrange0.second;
106 iter1 != clusrange1.second;
126 iter2 != clusrange2.second;
137 (trk.ClusterList).push_back(iter0->second);
138 (trk.ClusterList).push_back(iter1->second);
139 (trk.ClusterList).push_back(iter2->second);
142 if ( lyrs.size() > 3 )
144 bool found_clu3 =
false;
148 iter3 != clusrange3.second;
157 MvtxTrack tmp_trk = trk;
158 (tmp_trk.ClusterList).push_back(iter3->second);
159 trklst.push_back(tmp_trk);
162 if ( !found_clu3 ) trklst.push_back(trk);
167 trklst.push_back(trk);
181 std::cout <<
PHWHERE <<
" Running Ghost Rejection on "
182 << trklst.size() <<
" tracks" << std::endl;
185 std::multimap<TrkrDefs::cluskey, unsigned int> key_trk_map;
187 for (
unsigned int itrk = 0; itrk < trklst.size(); itrk++)
189 for (
unsigned int iclus = 0; iclus < trklst.at(itrk).ClusterList.size(); iclus++)
192 key_trk_map.insert(std::make_pair(ckey, itrk));
197 std::set<unsigned int> remove_set;
198 for (
auto iter = key_trk_map.begin(); iter != key_trk_map.end(); ++iter)
200 int ntrk = key_trk_map.count(iter->first);
205 std::cout <<
PHWHERE <<
" Tracks sharing cluster:" << ntrk << std::endl;
209 auto upiter = key_trk_map.upper_bound(iter->first);
212 double chi2_best = 9999.;
213 unsigned int idx_best = 0;
215 for (
auto jter = iter; jter != upiter; ++jter)
218 double chi2 = trklst.at(jter->second).chi2_xy + trklst.at(jter->second).chi2_zy;
219 if ( chi2 < chi2_best && chi2 > 0 )
222 idx_best = jter->second;
226 for (
auto jter = iter; jter != upiter; ++jter)
230 double chi2 = trklst.at(jter->second).chi2_xy + trklst.at(jter->second).chi2_zy;
232 <<
" trk idx: " << jter->second
234 <<
" m_xy:" << trklst.at(jter->second).m_xy;
236 if ( jter->second == idx_best)
238 std::cout <<
" <--- BEST" << std::endl;
242 std::cout << std::endl;
246 if ( jter->second != idx_best )
248 remove_set.insert(jter->second);
260 std::cout <<
PHWHERE <<
" List of idx to remove:";
261 for (
auto it = remove_set.begin();
it != remove_set.end(); ++
it)
263 std::cout <<
" " << *
it;
265 std::cout << std::endl;
269 for (
auto rit = remove_set.rbegin(); rit != remove_set.rend(); ++rit)
270 trklst.erase(trklst.begin() + *rit);
276 std::cout << PHWHERE <<
" Remaining tracks:" << std::endl;
277 for (
unsigned int i = 0;
i < trklst.size();
i++)
279 double chi2 = trklst.at(
i).chi2_xy + trklst.at(
i).chi2_zy;
280 std::cout <<
" chi2:" << chi2 <<
" m_xy:" << trklst.at(
i).m_xy << std::endl;
296 TMatrixD XT(X); XT.T();
297 TMatrixD
A = XT * L *
X;
298 TVectorD
b = XT * L *
y;
302 TMatrixD UT = svd.GetU(); UT.T();
305 TVectorD
s = svd.GetSig();
306 TMatrixD Sd(s.GetNrows(), s.GetNrows());
307 for (
int i = 0;
i < s.GetNrows();
i++)
308 Sd(
i,
i) =
s(
i) > 0 ? 1. /
s(
i) : 0.;
310 TVectorD beta = svd.GetV() * Sd * UT *
b;
323 int m = (trk.ClusterList).
size(),
n = 2;
330 for (
int iclus = 0; iclus <
m; iclus++)
334 x(iclus) = clus->
getX();
335 y(iclus) = clus->
getY();
338 X(iclus, 1) =
y(iclus);
340 Cinv(iclus, iclus) = 2.0 * sqrt(clus->
getSize(0, 0));
344 double x0 = beta(0),
c = beta(1);
350 for (
int iclus = 0; iclus <
m; iclus++)
354 double cx = clus->
getX();
355 double cy = clus->
getY();
357 double px = cy *
c +
x0;
360 trk.dx.push_back(dx);
362 chi2 += pow(dx, 2) / clus->
getError(0, 0);
378 int m = (trk.ClusterList).
size(),
n = 2;
385 for (
int iclus = 0; iclus <
m; iclus++)
389 z(iclus) = clus->
getZ();
390 y(iclus) = clus->
getY();
393 X(iclus, 1) =
y(iclus);
395 Cinv(iclus, iclus) = 2.0 * sqrt(clus->
getSize(2, 2));
399 double z0 = beta(0),
c = beta(1);
405 for (
int iclus = 0; iclus <
m; iclus++)
409 double cz = clus->
getZ();
410 double cy = clus->
getY();
412 double pz = cy *
c + z0;
415 trk.dz.push_back(dz);
417 chi2 += pow(dz, 2) / clus->
getError(2, 2);
428 std::cout <<
"===================================================" << std::endl;
429 std::cout <<
"== " << PHWHERE <<
" Found " << trklst.size() <<
" Track Candidates" << std::endl;
430 std::cout <<
"===================================================" << std::endl;
431 for (
unsigned int i = 0;
i < trklst.size();
i++)
433 std::cout <<
"== " <<
i << std::endl;
434 for (
unsigned int j = 0;
j < trklst.at(
i).ClusterList.size();
j++)
436 std::cout <<
" clus " <<
j
437 <<
" key:0x" << std::hex << trklst.at(
i).ClusterList.at(
j)->getClusKey() << std::dec
438 <<
" (" << trklst.at(
i).ClusterList.at(
j)->getX()
439 <<
", " << trklst.at(
i).ClusterList.at(
j)->getY()
440 <<
", " << trklst.at(
i).ClusterList.at(
j)->getZ()
442 <<
" dx:" << trklst.at(
i).dx.at(
j)
443 <<
" dz:" << trklst.at(
i).dz.at(
j)
447 <<
" m:" << trklst.at(
i).m_xy
448 <<
" b:" << trklst.at(
i).b_xy
449 <<
" chi2:" << trklst.at(
i).chi2_xy
452 <<
" m:" << trklst.at(
i).m_zy
453 <<
" b:" << trklst.at(
i).b_zy
454 <<
" chi2:" << trklst.at(
i).chi2_zy
457 std::cout <<
"===================================================" << std::endl;
463 return (y1 - y0) / (x1 -
x0);