9 #include <boost/io/ios_state.hpp>
62 if (m_nclusters <= m_tools->clustersmin() ||
70 if (tsos) dtsos.push_back(tsos);
77 if (tsos) dtsos.push_back(tsos);
83 if (tsos) dtsos.push_back(tsos);
116 if (tsos) dtsos.push_back(tsos);
123 if (tsos) dtsos.push_back(tsos);
129 if (tsos) dtsos.push_back(tsos);
161 if (tsos) dtsos.push_back(tsos);
165 if (tsos) dtsos.push_back(tsos);
167 int lastClusterElement = 0;
171 lastClusterElement = j;
175 if( lastClusterElement==0 || lastClusterElement==
i )
return dtsos;
182 if (tsos) dtsos.push_back(tsos);
188 if (tsos) dtsos.push_back(tsos);
207 if (tsos) dtsos.push_back(tsos);
214 if (tsos) dtsos.push_back(tsos);
220 if (tsos) dtsos.push_back(tsos);
254 if (tsos) dtsos.push_back(tsos);
258 if (tsos) dtsos.push_back(tsos);
260 int lastClusterElement = 0;
264 lastClusterElement = j;
268 if( lastClusterElement==0 || lastClusterElement==
i )
return dtsos;
275 if (tsos) dtsos.push_back(tsos);
281 if (tsos) dtsos.push_back(tsos);
292 return std::make_unique<Trk::FitQuality>(xi2, (
m_ndf - 5));
300 (std::multimap<const Trk::PrepRawData*,const Trk::Track*>& map)
const
304 std::multimap<const Trk::PrepRawData*,const Trk::Track*>::const_iterator
305 ti,
t[100],te = map.end();
313 t[
n] = map.find(prd[
n]);
314 if (
t[
n]==te)
return true;
318 t[
n] = map.find(prd[
n]);
319 if (
t[
n]==te)
return true;
326 for (
int i=0;
i!=
n; ++
i) {
328 for (ti=
t[
i]; ti!=te; ++ti) {
329 if ( (*ti).first != prd[
i] )
break;
330 int ncl = (*ti).second->measurementsOnTrack()->size();
331 if (ncl > nclmax) nclmax = ncl;
333 if (nclt > nclmax)
return true;
342 std::ostream& InDet::operator <<
354 boost::io::ios_all_saver ias(
out);
357 out<<
"Trajectory does not exist"<<std::endl; ias.restore();
361 out<<
"Trajectory is wrong"<<std::endl; ias.restore();
365 out<<
"|--------------------------------------------------------------------------------------------------------|"
378 <<std::setw(2)<<
m_ndf<<
" weighted clusters and quality = "<<std::setw(12)<<std::setprecision(5)<<
quality()
381 out<<
"| Has number of holes before, inside, after and gap= "
390 out<<
"|---|--|---|-----|-|-|---------|---------|---------|---------|----------|---------|-|--|--|--|-|-|-|-|-|-|---------|"
392 out<<
"| # |SS| D| Ncl |C|O| Xi2F | Xi2B | Az.Angle| Radius | pT(GeV) | dZ/dR |N|In|Lf|Lb|S|D|H|G|H|G| Step |"
394 out<<
"|---|--|---|-----|-|-|---------|---------|---------|---------|----------|---------|-|--|--|--|-|-|-|-|-|-|---------|"
401 std::string DET =
"D ";
403 std::string DE =
" ";
407 if (D->
isBarrel()) DET =
"Pb";
else DET =
"Pe";
409 else if (D->
isSCT()) {
410 if (D->
isBarrel()) DET =
"Sb";
else DET =
"Se";
430 <<std::setw(2)<<DET <<
"|"
431 <<std::setw(5)<<
c <<
"|"
450 ra = sqrt(gp.x()*gp.x()+gp.y()*gp.y());
451 fa = atan2(gp.y(),gp.x());
458 ra = sqrt(gp.x()*gp.x()+gp.y()*gp.y());
459 fa = atan2(gp.y(),gp.x());
469 ra = sqrt(gp.x()*gp.x()+gp.y()*gp.y());
470 fa = atan2(gp.y(),gp.x());
477 ra = sqrt(gp.x()*gp.x()+gp.y()*gp.y());
478 fa = atan2(gp.y(),gp.x());
495 ra = sqrt(gp.x()*gp.x()+gp.y()*gp.y());
496 fa = atan2(gp.y(),gp.x());
497 pt =
SM.momentum().perp();
501 out<<std::setw( 9)<<std::setprecision(4)<<fa <<
"|";
502 out<<std::setw( 9)<<std::setprecision(4)<<ra <<
"|";
503 out<<std::setw(10)<<std::setprecision(4)<<
pt*.001<<
"|";
504 out<<std::setw( 9)<<std::setprecision(4)<<tz <<
"|";
539 out<<
"|---|--|---|-----|-|-|---------|---------|---------|---------|----------|---------|-|--|--|--|-|-|-|-|-|-|---------|"
551 std::vector<const InDet::SiCluster*> & Cl,
552 std::vector<const InDet::SiDetElementBoundaryLink_xk*>& DE,
553 const EventContext & ctx)
557 InDet::SiClusterCollection::const_iterator sib,sie;
559 std::vector<const InDet::SiCluster*>
::iterator s=Cl.begin();
563 if(!
m_elements[
n].firstTrajectorElement(Tp,ctx))
return 0.;
565 for(++
r;
r!=
re; ++
r) {
585 std::vector<const InDet::SiCluster*> & lSiCluster,
586 std::vector<const InDet::SiDetElementBoundaryLink_xk*>& DE ,
588 const EventContext & ctx )
608 double Rdead = 142.5;
617 if (lSiCluster.size() < 2)
return false;
626 for (iter_boundaryLink=DE.begin(); iter_boundaryLink!=endBoundaryLinks; ++iter_boundaryLink) {
636 if (detectorElement->
isPixel()) {
642 if(not initDeadMaterial) {
644 double R = pla->
center().perp();
646 initDeadMaterial =
true;
657 InDet::PixelClusterCollection::const_iterator iter_PixelClusterColl, iter_PixelClusterCollEnd;
661 if (clustersOnElement!=
nullptr && clustersOnElement->begin()!=clustersOnElement->end()) {
664 iter_PixelClusterColl = clustersOnElement->begin();
665 iter_PixelClusterCollEnd = clustersOnElement->end();
669 for (iter_cluster=lSiCluster.begin(); iter_cluster!=lSiCluster.end(); ++iter_cluster) {
671 if ((*iter_cluster)->detectorElement()==detectorElement) {
685 theCluster=(*iter_cluster);
687 iter_cluster=lSiCluster.erase(iter_cluster);
696 bool valid_set =
m_elements[
m_nElements].
set(1,(*iter_boundaryLink),iter_PixelClusterColl,iter_PixelClusterCollEnd,theCluster,ctx);
704 bool valid_set =
m_elements[
m_nElements].
set(0,(*iter_boundaryLink),iter_PixelClusterColl,iter_PixelClusterCollEnd,theCluster,ctx);
722 InDet::SCT_ClusterCollection::const_iterator iter_stripClusterColl, iter_StripClusterCollEnd;
727 if (clustersOnElement!=
nullptr && clustersOnElement->begin()!=clustersOnElement->end()) {
729 iter_stripClusterColl = clustersOnElement->begin();
730 iter_StripClusterCollEnd = clustersOnElement->end();
733 for (iter_cluster=lSiCluster.begin(); iter_cluster!=lSiCluster.end(); ++iter_cluster) {
734 if ((*iter_cluster)->detectorElement()==detectorElement) {
749 theCluster=(*iter_cluster);
750 iter_cluster=lSiCluster.erase(iter_cluster);
757 bool valid_set =
m_elements[
m_nElements].
set(1,(*iter_boundaryLink),iter_stripClusterColl,iter_StripClusterCollEnd,theCluster,ctx);
765 bool valid_set =
m_elements[
m_nElements].
set(0,(*iter_boundaryLink),iter_stripClusterColl,iter_StripClusterCollEnd,theCluster,ctx);
783 if (!
m_elements[
up].firstTrajectorElement(Tp, ctx))
return false;
789 else if (theCluster) {
808 if (ndfwrong > 3)
return false;
844 if (!lSiCluster.empty()) {
918 std::vector<const InDet::SiDetElementBoundaryLink_xk*> & DE ,
919 std::multimap<const Trk::PrepRawData*,const Trk::Track*>&
PT ,
920 std::vector<const InDet::SiCluster*> & lSiCluster,
921 const EventContext& ctx)
926 std::multimap<double,const InDet::SiCluster*> xi2cluster;
929 std::multimap<const Trk::PrepRawData*,const Trk::Track*>::const_iterator
t, te =
PT.end();
934 for (iter_boundaryLink=DE.begin(); iter_boundaryLink!=endBoundaryLinks; ++iter_boundaryLink) {
939 bool sct = detectorElement->
isSCT();
942 InDet::PixelClusterCollection::const_iterator sib, sie;
945 if (
w!=
nullptr &&
w->begin()!=
w->end()) {
951 if (!
m_elements[0].ForwardPropagationForClusterSeach(
m_nElements,Tp,(*iter_boundaryLink),sib,sie,ctx))
return false;
953 InDet::SCT_ClusterCollection::const_iterator sib, sie;
956 if (
w!=
nullptr &&
w->begin()!=
w->end()) {
962 if (!
m_elements[0].ForwardPropagationForClusterSeach(
m_nElements,Tp,(*iter_boundaryLink),sib,sie,ctx))
return false;
971 if (
t!=te && (*t).second->measurementsOnTrack()->size() >= 10)
continue;
976 if (
x <= xi2Cut) xi2cluster.insert(std::make_pair(
x,
m_elements[0].linkF(
i).cluster()));
982 if (xi2cluster.size() < 3)
return false;
986 for (; xc!=xce; ++xc) {
987 lSiCluster.push_back((*xc).second);
988 (*xc).second->detectorElement()->isSCT() ?
m_ndf+=1 :
m_ndf+=2;
989 if (
m_ndf >= ndfCut )
break;
1002 const std::vector<Amg::Vector3D> & Gp ,
1003 std::vector<const InDet::SiDetElementBoundaryLink_xk*> & DE ,
1004 std::multimap<const Trk::PrepRawData*,const Trk::Track*>&
PT ,
1005 std::vector<const InDet::SiCluster*> & lSiCluster)
1008 std::vector<Amg::Vector3D>::const_iterator
g,gb = Gp.begin(), ge = Gp.end();
1009 InDet::PixelClusterCollection::const_iterator pib, pie;
1010 InDet::SCT_ClusterCollection::const_iterator sib, sie;
1011 std::multimap<const Trk::PrepRawData*,const Trk::Track*>::const_iterator
t, te =
PT.end();
1015 double pv[ 5]={0.,0.,0.,0.,0.};
1016 double cv[15]={ .1 ,
1020 0. , 0., 0., 0.,.00001};
1022 double xi2Cut = 10.;
1025 for (; iter_boundaryLink!=endBoundaryLinks; ++iter_boundaryLink) {
1034 double Ax[3] = {tr(0,0),tr(1,0),tr(2,0)};
1035 double Ay[3] = {tr(0,1),tr(1,1),tr(2,1)};
1036 double Az[3] = {tr(0,2),tr(1,2),tr(2,2)};
1037 double x0 = tr(0,3);
1038 double y0 = tr(1,3);
1039 double z0 = tr(2,3);
1040 double zcut = .001 ;
1042 bool sct =
d->isSCT();
1045 if (
w!=
nullptr &&
w->begin()!=
w->end()) {
1054 if (
w!=
nullptr &&
w->begin()!=
w->end()) {
1062 for (
g=gb;
g!=ge; ++
g) {
1064 double dx = (*g).x()-x0;
1065 double dy = (*g).y()-y0;
1066 double dz = (*g).z()-
z0;
1067 double z =
dx*Az[0]+
dy*Az[1]+dz*Az[2];
1068 if (std::abs(
z) > zcut)
continue;
1070 pv[0] =
dx*Ax[0]+
dy*Ax[1]+dz*Ax[2];
1071 pv[1] =
dx*Ay[0]+
dy*Ay[1]+dz*Ay[2];
1078 if (!
c ||
m_elements[0].xi2F() > xi2Cut)
continue;
1081 if (
t!=te && (*t).second->measurementsOnTrack()->size() >= 10)
continue;
1084 lSiCluster.push_back(
c);
1190 if (L==0)
return true;
1201 const int itm = itmax-1 ;
1213 double Xi2best = 0. ;
1217 for (;
it!=itmax; ++
it) {
1220 int lastElementWithExpHit =
F;
1222 for (--
F;
F>=0; --
F) {
1230 lastElementWithExpHit =
F;
1233 else if (Ef.
inside() < 0) {
1234 lastElementWithExpHit =
F;
1239 if (Ef.
ndist() > ndcut ||
1241 (
nm == nclbest && Ef.
xi2totalB() > Xi2best) )
break;
1258 if ( (
q > qbest) || (
q==qbest &&
X < Xi2best ) ) {
1269 if (fl==0 && nd < ndcut) ndcut = nd;
1271 if (fl!=0 || nd > 0 ||
np < 3) {
1274 for (
int i=
l;
i!=L; ++
i) {
1279 MPbest[++lbest] =
i;
1282 XI2B[nbest] = Ei.
xi2B();
1284 TE[nbest++] = lbest;
1285 ndfbest += Ei.
ndf();
1294 for (
int i=fl;
i!=L; ++
i) {
1298 if (Ei.
inside() <= 0 && ++lbest >=0 ) {
1299 MPbest[lbest] = lbest;
1302 XI2B[nbest] = Ei.
xi2B();
1304 TE[nbest++] = lbest;
1305 ndfbest += Ei.
ndf();
1334 for (;
l < L; ++
l) {
1340 if (Ei.
dist() < -2. && Ei.
ndist() > ndcut-1 )
continue;
1343 if (
nm < nclbest || (
nm == nclbest && Xn > Xi2best))
continue;
1351 if (
it == itmax) --
it;
1352 if (!nbest)
return true;
1373 if (itbest==
it)
return true;
1375 for (
int n = L-1-
dn;
n>=0; --
n) {
1379 for (;
m>=0; --
m)
if (TE[
m]==
n)
break;
1387 if (--nbest==0)
break;
1428 for (++extensionStartIndex; extensionStartIndex<=
m_lastElement; ++extensionStartIndex) {
1440 for (++extensionStartIndex; extensionStartIndex<=
m_lastElement; ++extensionStartIndex) {
1462 --extensionStartIndex;
1466 if ( extensionStartIndex== lastElementOnTraj)
return true;
1482 const int itm = itmax-1 ;
1503 int lbest = extensionStartIndex ;
1505 int index_currentElement = extensionStartIndex ;
1507 int M = extensionStartIndex ;
1509 MP [M] = extensionStartIndex ;
1510 double Xi2best = 0. ;
1512 const double dfmax = 2.2 ;
1521 for (; iteration!=itmax; ++iteration) {
1523 int lastElementWithSeenHit = index_currentElement;
1524 int lastElementWithExpHit = index_currentElement;
1526 int index_previousElement = index_currentElement;
1528 int mLastCluster = M;
1530 int Cm = nclbest-lastElementOnTraj;
1532 bool haveHole =
false;
1536 for (++index_currentElement; index_currentElement!=
m_nElements; ++index_currentElement) {
1547 if (!currentElement.
isBarrel() || index_previousElement!=index_currentElement-1)
break;
1552 int index_auxElement = index_currentElement;
1553 for (; index_auxElement!=
m_nElements; ++index_auxElement) {
1562 index_currentElement = index_auxElement-1;
1571 index_previousElement = index_currentElement;
1576 MP[++M] = index_currentElement;
1579 if (currentElement.
cluster()) {
1582 double df = std::abs(currentElement.
parametersUF().parameters()[2]-f0);
1586 if (
df > dfmax)
break;
1589 lastElementWithExpHit = index_currentElement;
1590 lastElementWithSeenHit = index_currentElement;
1596 else if (currentElement.
inside() < 0 ) {
1597 lastElementWithExpHit=index_currentElement;
1599 if (currentElement.
nholesF() > maxholes || currentElement.
dholesF() > maxdholes)
break;
1604 double df = std::abs(currentElement.
parametersPF().parameters()[2]-f0);
1606 if (
df > dfmax )
break;
1613 double df = std::abs(currentElement.
parametersPF().parameters()[2]-f0);
1615 if (
df > dfmax)
break;
1619 int nm = currentElement.
nclustersF()-index_currentElement;
1621 if ( currentElement.
ndist() > ndcut
1623 || (
nm == Cm && currentElement.
xi2totalF() > Xi2best)
1643 int lastElementProcessed = index_currentElement;
1645 if (lastElementProcessed==
m_nElements) --lastElementProcessed;
1655 if ( (
q > qbest) || (
q==qbest &&
X < Xi2best ) ) {
1670 lbest = extensionStartIndex ;
1677 if (lastElementProcessed==lastElementOnTraj && nd < ndcut) ndcut = nd;
1680 for (
int j=extensionStartIndex+1; j<=mLastCluster; ++j) {
1687 MPbest[++lbest] =
i;
1692 TE[nbest++] = lbest;
1694 ndfbest += Ei.
ndf();
1702 if ( (nclbest >= 14 && !haveHole) || (lastElementProcessed==lastElementOnTraj && ndbest == 0))
break;
1706 index_currentElement = -1;
1710 int nb = lastElementOnTraj-nclbest-1;
1721 for (
int j=mLastCluster; j!=extensionStartIndex; --j) {
1738 if (Ei.
dist() < -2. && Ei.
ndist() > ndcut-1)
continue;
1745 if (
nm < 0 || (
nm == 0 && Xn > Xi2best))
continue;
1747 index_currentElement =
i;
1754 if (index_currentElement < 0 )
break;
1761 if (iteration == itmax) --iteration;
1764 if (!nbest)
return true;
1782 if (itbest==iteration)
return true;
1787 index_currentElement = -1;
1796 for (;
m!=nbest; ++
m)
if (TE[
m]==
n)
break;
1802 if (index_currentElement<0) {
1804 index_currentElement=
n;
1815 if (++
mb == nbest)
break;
1823 if (index_currentElement<0) {
1824 index_currentElement =
n;
1840 if (index_currentElement < 0 ||
m_lastElement == index_currentElement) {
1846 for (++index_currentElement; index_currentElement<=
m_lastElement; ++index_currentElement) {
1860 (std::vector<const InDet::SiCluster*>& Cl)
const
1963 else if ((
step*stp) < 0.) {
order =
false;
break;}
1966 if (
order)
return true;
1971 double rad = gp.x()*gp.x()+gp.y()*gp.y();
1980 double R = gp.x()*gp.x()+gp.y()*gp.y();
1981 if (R <
rad)
return false;
2020 for (;
n<=LA; ++
n) {
2029 nc =
true; so =
false;
2041 for (;
n<=LA; ++
n) {
2050 nc =
true; so =
false;
2062 for (;
n<= LA; ++
n) {
2074 for (;
m>=
n ; --
m) {
2191 if (fE==lE ||
m_nclusters+m_nclustersNoAdd < m_tools->clustersmin())
return -100.;
2218 if (tsos) dtsos.push_back(tsos);
2232 if (
s<=1)
return 0.;
2243 bool prevIsSctHole =
false;
2252 bool isPix = theElement.
ndf() == 2;
2253 bool isSCTHole=
false;
2258 prevIsSctHole =
false;
2266 switch (boundaryStatus){
2299 if (isSCTHole && prevIsSctHole){
2300 prevIsSctHole =
false;
2304 prevIsSctHole = isSCTHole;