diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.h b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.h
index 4b06836f77ae05afcb9bafea75352fa8ffcf1591..843316c22409df7a9bd12f75e276d89292890691 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.h
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.h
@@ -47,7 +47,8 @@ namespace InDet{
const int& detstatus () const {return m_detstatus;}
const int& inside () const {return m_inside; }
- const int& ndist () const {return m_ndist ; }
+ /// number of crossed without hit - dead + holes
+ const int& ndist () const {return m_nMissing ; }
const int& nlinksF () const {return m_nlinksForward; }
const int& nlinksB () const {return m_nlinksBackward; }
const int& nholesF () const {return m_nholesForward; }
@@ -418,7 +419,7 @@ namespace InDet{
int m_status ;
int m_detstatus ; //!< 0 (no clusters)
int m_inside ;
- int m_ndist ;
+ int m_nMissing ;
int m_nlinksForward ;
int m_nlinksBackward ;
int m_nholesForward ;
diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.icc b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.icc
index bd1cee22895658dae4f1634dffdcdee4a80adc98..d1c79e2be1cb004542d07b6b0564f92143ce4c47 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.icc
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.icc
@@ -11,24 +11,33 @@
template<typename T>
int InDet::SiTrajectoryElement_xk::searchClustersSub
(Trk::PatternTrackParameters& Tp, SiClusterLink_xk* L) {
+ /// Case a): No PRD association in use
if (not m_useassoTool) {
+ /// a1): No stereo
if (not m_stereo) {
if (m_ndf == 2) {
return searchClustersWithoutStereoPIX<T>(Tp, L);
} else {
return searchClustersWithoutStereoSCT<T>(Tp, L);
}
- } else {
+ }
+ /// a2): With stereo
+ else {
return searchClustersWithStereo<T>(Tp, L);
}
- } else {
+ }
+ /// Case b) Using PRD association tool
+ else {
+ /// b1): no stereo
if (not m_stereo) {
if (m_ndf == 2) {
return searchClustersWithoutStereoAssPIX<T>(Tp, L,*m_prdToTrackMap);
} else {
return searchClustersWithoutStereoAssSCT<T>(Tp, L,*m_prdToTrackMap);
}
- } else {
+ }
+ /// b2): with stereo
+ else {
return searchClustersWithStereoAss<T>(Tp, L, *m_prdToTrackMap);
}
}
@@ -44,22 +53,27 @@ int InDet::SiTrajectoryElement_xk::searchClustersWithStereo
{
if (m_detstatus<=0) return 0;
- int nl = 0;
+ int nLinksFound = 0;
+ /// predicted local position
double P0 = Tp.par()[0];
double P1 = Tp.par()[1];
- double PV0 = Tp.cov()[0];
- double PV1 = Tp.cov()[1];
- double PV2 = Tp.cov()[2];
- double Xc = m_xi2maxlink;
- double Xl = m_xi2maxlink;
+ /// covariance for prediction
+ double PV0 = Tp.cov()[0]; /// X,X
+ double PV1 = Tp.cov()[1]; /// X,Y
+ double PV2 = Tp.cov()[2]; /// Y,Y
+ /// chi² cut values
+ double Xc = m_xi2maxlink;
+ double bestX2 = m_xi2maxlink;
double Xm = m_xi2max ;
+ /// holder for best cluster seen not pasing cuts
const InDet::SiCluster* cl = nullptr;
T* sibegin = std::any_cast<T>(&m_sibegin);
T* siend = std::any_cast<T>(&m_siend);
if (sibegin==nullptr or siend==nullptr) return 0;
+ /// loop over all clusters on this element
for (T p=*sibegin; p!=*siend; ++p) {
const InDet::SiCluster* c = static_cast<const InDet::SiCluster*>(*p);
const Amg::Vector2D& M = c->localPosition();
@@ -68,31 +82,73 @@ int InDet::SiTrajectoryElement_xk::searchClustersWithStereo
double MV0 = V(0,0);
double MV1 = V(1,0);
double MV2 = V(1,1);
+ /// add the two covariances
double v0 = MV0+PV0;
double v1 = MV1+PV1;
double v2 = MV2+PV2;
+ /// get the distance from the cluster to the predicted location
double r0 = M[0]-P0;
double r1 = M[1]-P1;
- double d = v0*v2-v1*v1; if(d<=0.) continue; d=1./d;
+ /// determinant of the cov matrix
+ double d = v0*v2-v1*v1;
+ if(d<=0.) continue;
+ d=1./d;
+
+ /// chi² calculation
double x = (r0*(r0*v2-r1*v1)+r1*(r1*v0-r0*v1))*d;
- if(x > Xc) continue;
+ if(x > Xc) continue; /// if we don't satisfy the minimum, don't bother further
+ /// refine y estimate
r1 = fabs(r1+d*((PV1*v2-PV2*v1)*r0+(PV2*v0-PV1*v1)*r1));
x -= (r1*r1)/MV2 ;
r1 -= m_halflength ;
- if(r1 > 0. && (x+=((r1*r1)/PV2)) > Xc) continue;
-
+ /// update the chi² estimate, and check again the cut
+ if(r1 > 0.){
+ x+=(r1*r1)/PV2;
+ if (x > Xc) continue;
+ }
+ /// if we satisfy the chi2 criterion:
if(x < Xm) {
+ /// create a cluster link to this cluster
InDet::SiClusterLink_xk l(c,x);
- for(int i=0; i!=nl; ++i) L[i].Comparison(l);
- if(nl<10) L[nl++]=l; else Xm=L[9].xi2();
+
+ /// This misleadingly named method actually is a sorting
+ /// Mechanism: Our current link will take the place of the one
+ /// it replaces in ascending chi², and we will be
+ /// left with a link to the worst chi² one
+ for(int i=0; i!=nLinksFound; ++i) L[i].Comparison(l);
+ /// Now our link points not the the cluster we just processed but the worst
+ /// so far
+ /// If we have less than 10 saved, save it in the next slot
+ if(nLinksFound<10){
+ L[nLinksFound++]=l;
+ } else {
+ /// otherwise, don't add it, and update the chi2 cut to
+ /// the worst one in our list, to get rid of further
+ /// candidates early on.
+ Xm=L[9].xi2();
+ }
+ /// finally, update the chi2 cut - don't bother with anything worse
+ /// than the current candidate by 6 units in chi2
Xc = Xm+6.;
+ } /// done with branch if we satisfy the chi2 criterion
+ /// if we have no links yet and this is the best chi2 so far
+ else if(!nLinksFound && x < bestX2) {
+ /// update best chi2
+ bestX2 = x;
+ /// update max cut
+ Xc = x+6.;
+ /// update best cluster pointer
+ cl = c;
}
- else if(!nl && x < Xl) {Xl = x; Xc = x+6.; cl = c;}
- }
- if(cl && !nl) {L[nl++].Set(cl,Xl);}
- return nl;
+ } /// end cluster loop
+ /// If we didn't find any cluster satisfying the chi2 cut,
+ /// but have at least one cluster,
+ /// we write one link to the best (of the bad) candidate
+ if(cl && !nLinksFound) {L[nLinksFound++].Set(cl,bestX2);}
+ /// return number of found links
+ return nLinksFound;
}
///////////////////////////////////////////////////////////////////
@@ -105,49 +161,87 @@ int InDet::SiTrajectoryElement_xk::searchClustersWithoutStereoPIX
{
if (m_detstatus<=0) return 0;
- int nl = 0;
+ int nLinksFound = 0;
+ /// predicted local position
double P0 = Tp.par()[0];
double P1 = Tp.par()[1];
+ /// covariance of prediction
double PV0 = Tp.cov()[0];
double PV1 = Tp.cov()[1];
double PV2 = Tp.cov()[2];
+ /// chi2 cuts
double Xc = m_xi2maxlink;
double Xm = m_xi2max ;
+ /// best cluster seen
const InDet::SiCluster* cl = nullptr;
T* sibegin = std::any_cast<T>(&m_sibegin);
T* siend = std::any_cast<T>(&m_siend);
if (sibegin==nullptr or siend==nullptr) return 0;
+ /// loop over all hits
for (T p=*sibegin; p!=*siend; ++p) {
const InDet::SiCluster* c = static_cast<const InDet::SiCluster*>(*p);
const Amg::Vector2D& M = c->localPosition();
+ /// estimate the covariance by the pitch
+ /// Factors 1/sqrt(12) will be added later
double MV0 = c->width().phiR();
double MV2 = c->width().z ();
- double r0 = M[0]-P0, r02 = r0*r0;
- double r1 = M[1]-P1, r12 = r1*r1;
-
- double v0 = .08333*(MV0*MV0)+PV0; if(r02 >(Xc*v0)) continue;
- double v2 = .08333*(MV2*MV2)+PV2; if(r12 >(Xc*v2)) continue;
+ /// squared distance between prediction and hit location,
+ double r0 = M[0]-P0,
+ r02 = r0*r0;
+ double r1 = M[1]-P1,
+ r12 = r1*r1;
+
+ /// sum cov on X - using 1/sqrt(12) update of hit cov
+ double v0 = .08333*(MV0*MV0)+PV0;
+ /// bail out if X alone sufficient to fail chi2
+ if(r02 >(Xc*v0)) continue;
+
+ /// sum cov on Y - using 1/sqrt(12) update of hit cov
+ double v2 = .08333*(MV2*MV2)+PV2;
+ /// bail out if Y alone sufficient to fail chi2
+ if(r12 >(Xc*v2)) continue;
double v1 = PV1;
- double d = v0*v2-v1*v1; if( d<=0. ) continue;
+ /// determinant of cov matrix
+ double d = v0*v2-v1*v1;
+ if( d<=0. ) continue;
+ /// chi2
double x = (r02*v2+r12*v0-(r0*r1)*(2.*v1))/d;
+ /// skip clusters failing the cut
if(x>Xc) continue;
+ /// if we satisfy the tighter requirement:
if(x < Xm) {
+ /// set up a cluster link
InDet::SiClusterLink_xk l(c,x);
- for(int i=0; i!=nl; ++i) L[i].Comparison(l);
- if(nl<10) L[nl++]=l; else Xm=L[9].xi2();
+ /// This misleadingly named method actually is a sorting
+ /// Mechanism: Our current link will take the place of the one
+ /// it replaces in ascending chi², and we will be
+ /// left with a link to the worst chi² one
+ for(int i=0; i!=nLinksFound; ++i) L[i].Comparison(l);
+ /// and place the worst link (remember that l is now a different link, see above)!
+ /// If within capacity, increment counter and add
+ if(nLinksFound<10) L[nLinksFound++]=l;
+ /// otherwise increment cut to avoid bothering with lower qualities in the future
+ else Xm=L[9].xi2();
+ /// update looser chi2 cut to match the tighter one - we found at least one good cluster
Xc = Xm;
}
- else if(!nl) {Xc = x; cl = c;}
+ /// if we found no clusters yet satisfying the chi2 cut, use this one
+ else if(!nLinksFound) {
+ Xc = x; cl = c;
+ }
+ } /// end of cluster loop
+ /// if we found no cluster satisfying the chi2, we return one link to the best (least bad) we saw
+ if(cl && !nLinksFound) {
+ L[nLinksFound++].Set(cl,Xc);
}
- if(cl && !nl) {L[nl++].Set(cl,Xc);}
- return nl;
+ return nLinksFound;
}
///////////////////////////////////////////////////////////////////
@@ -160,57 +254,84 @@ int InDet::SiTrajectoryElement_xk::searchClustersWithoutStereoSCT
{
if (m_detstatus<=0) return 0;
- int nl = 0;
+ int nLinksFound = 0;
+ /// predicted local position
double P0 = Tp.par()[0];
double P1 = Tp.par()[1];
+ /// and cov of prediction
double PV0 = Tp.cov()[0];
double PV1 = Tp.cov()[1];
double PV2 = Tp.cov()[2];
+ /// chi2 cuts
double Xc = m_xi2maxlink;
double Xm = m_xi2max ;
-
+ /// best cluster seen
const InDet::SiCluster* cl = nullptr;
T* sibegin = std::any_cast<T>(&m_sibegin);
T* siend = std::any_cast<T>(&m_siend);
if (sibegin==nullptr or siend==nullptr) return 0;
+ /// loop over all hits on this element
for (T p=*sibegin; p!=*siend; ++p) {
const InDet::SiCluster* c = static_cast<const InDet::SiCluster*>(*p);
const Amg::Vector2D& M = c->localPosition();
+ /// estimate covariance by pitch
double MV0 = c->width().phiR() ;
+ /// include 1/sqrt(12) factor and add to contribution from pred
double v0 = .08333*(MV0*MV0)+PV0;
+ /// get residual in X
double r0 = M[0]-P0;
+ /// one over cov
double d = 1./v0;
+ /// chi2
double x = (r0*r0)*d;
+ /// stop if chi2 too large
if(x>Xc) continue;
-
+ /// estimate of Y residual
double dP1 = (P1-M[1])+PV1*d*r0;
+ /// if residual beyond half-length: assume hit is on very end of strip and
+ /// re-evaluate chi2 with a penalty term on the Y coordinate
if(fabs(dP1) > m_halflength) {
-
- double r1;
- dP1 > m_halflength ? r1 = m_halflength-P1 : r1 = -(m_halflength+P1);
-
+ double r1 = (dP1 > m_halflength ?m_halflength-P1 : -(m_halflength+P1));
double v1 = PV1;
double v2 = PV2;
- d = v0*v2-v1*v1 ; if(d<=0.) continue;
+ /// determinant of 2D covariance (only contribution from hit on X,X)
+ d = v0*v2-v1*v1 ;
+ if(d<=0.) continue;
+ /// updated chi2 estimate, with penalty for being at end of strip
x = (r0*(r0*v2-r1*v1)+r1*(r1*v0-r0*v1))/d;
+ /// and re-apply cut
if(x>Xc) continue;
- }
+ } /// end of branch for Y residual beyond half-length
+ /// if we satisfy the minimum criterion:
if(x < Xm) {
+ /// create a cluster link for our firend
InDet::SiClusterLink_xk l(c,x);
- for(int i=0; i!=nl; ++i) L[i].Comparison(l);
- if(nl<10) L[nl++]=l; else Xm=L[9].xi2();
+ /// This misleadingly named method actually is a sorting
+ /// Mechanism: Our current link will take the place of the one
+ /// it replaces in ascending chi², and we will be
+ /// left with a link to the worst chi² one
+ for(int i=0; i!=nLinksFound; ++i) L[i].Comparison(l);
+ /// and place the worst link so far at the end if space.
+ if(nLinksFound<10) L[nLinksFound++]=l;
+ /// otherwise, update chi2 cut to not bother with worse clusters anymore
+ else Xm=L[9].xi2();
+ /// update looser chi2 to match the tighter one - we have at least one good link now
Xc = Xm;
}
- else if(!nl) {Xc = x; cl = c;}
+ /// if we didn't find any good links yet, keep this one as best so far
+ /// update looser chi2 cut to only keep better candidates in the future
+ else if(!nLinksFound) {Xc = x; cl = c;}
}
- if(cl && !nl) {L[nl++].Set(cl,Xc);}
- return nl;
+ /// if we found no clusters satisfying the chi2 cut, we return
+ /// a single link to the best (least bad) candidate
+ if(cl && !nLinksFound) {L[nLinksFound++].Set(cl,Xc);}
+ return nLinksFound;
}
///////////////////////////////////////////////////////////////////
@@ -221,9 +342,12 @@ template <typename T>
int InDet::SiTrajectoryElement_xk::searchClustersWithStereoAss
(Trk::PatternTrackParameters& Tp,InDet::SiClusterLink_xk* L, const Trk::PRDtoTrackMap &prd_to_track_map)
{
+
+ /// for detailed documentation, please see searchClustersWithStereo
+
if (m_detstatus<=0) return 0;
- int nl = 0;
+ int nLinksFound = 0;
double P0 = Tp.par()[0];
double P1 = Tp.par()[1];
double PV0 = Tp.cov()[0];
@@ -241,7 +365,9 @@ int InDet::SiTrajectoryElement_xk::searchClustersWithStereoAss
for (T p=*sibegin; p!=*siend; ++p) {
const InDet::SiCluster* c = static_cast<const InDet::SiCluster*>(*p);
+ /// check the PRD association map, skip clusters already in use
if (prd_to_track_map.isUsed(*c)) continue;
+
const Amg::Vector2D& M = c->localPosition();
const Amg::MatrixX& V = c->localCovariance();
@@ -267,15 +393,19 @@ int InDet::SiTrajectoryElement_xk::searchClustersWithStereoAss
if(x < Xm) {
InDet::SiClusterLink_xk l(c,x);
- for(int i=0; i!=nl; ++i) L[i].Comparison(l);
- if(nl<10) L[nl++]=l;
+ /// This misleadingly named method actually is a sorting
+ /// Mechanism: Our current link will take the place of the one
+ /// it replaces in ascending chi², and we will be
+ /// left with a link to the worst chi² one
+ for(int i=0; i!=nLinksFound; ++i) L[i].Comparison(l);
+ if(nLinksFound<10) L[nLinksFound++]=l;
else Xm=L[9].xi2();
Xc = Xm+6.;
}
- else if(!nl && x < Xl) {Xl = x; Xc = x+6.; cl = c;}
+ else if(!nLinksFound && x < Xl) {Xl = x; Xc = x+6.; cl = c;}
}
- if(cl && !nl) {L[nl++].Set(cl,Xl);}
- return nl;
+ if(cl && !nLinksFound) {L[nLinksFound++].Set(cl,Xl);}
+ return nLinksFound;
}
///////////////////////////////////////////////////////////////////
@@ -286,9 +416,10 @@ template <typename T>
int InDet::SiTrajectoryElement_xk::searchClustersWithoutStereoAssPIX
(Trk::PatternTrackParameters& Tp,InDet::SiClusterLink_xk* L, const Trk::PRDtoTrackMap &prd_to_track_map)
{
+ /// for detailed documentation, please see searchClustersWithoutStereoPIX
if (m_detstatus<=0) return 0;
- int nl = 0;
+ int nLinksFound = 0;
double P0 = Tp.par()[0];
double P1 = Tp.par()[1];
double PV0 = Tp.cov()[0];
@@ -324,14 +455,18 @@ int InDet::SiTrajectoryElement_xk::searchClustersWithoutStereoAssPIX
if(x < Xm) {
InDet::SiClusterLink_xk l(c,x);
- for(int i=0; i!=nl; ++i) L[i].Comparison(l);
- if(nl<10) L[nl++]=l; else Xm=L[9].xi2();
+ /// This misleadingly named method actually is a sorting
+ /// Mechanism: Our current link will take the place of the one
+ /// it replaces in ascending chi², and we will be
+ /// left with a link to the worst chi² one
+ for(int i=0; i!=nLinksFound; ++i) L[i].Comparison(l);
+ if(nLinksFound<10) L[nLinksFound++]=l; else Xm=L[9].xi2();
Xc = Xm;
}
- else if(!nl) {Xc = x; cl = c;}
+ else if(!nLinksFound) {Xc = x; cl = c;}
}
- if(cl && !nl) {L[nl++].Set(cl,Xc);}
- return nl;
+ if(cl && !nLinksFound) {L[nLinksFound++].Set(cl,Xc);}
+ return nLinksFound;
}
///////////////////////////////////////////////////////////////////
@@ -342,9 +477,10 @@ template <typename T>
int InDet::SiTrajectoryElement_xk::searchClustersWithoutStereoAssSCT
(Trk::PatternTrackParameters& Tp,InDet::SiClusterLink_xk* L, const Trk::PRDtoTrackMap &prd_to_track_map)
{
+ /// for detailed documentation, please see searchClustersWithoutStereoSCT
if (m_detstatus<=0) return 0;
- int nl = 0;
+ int nLinksFound = 0;
double P0 = Tp.par()[0];
double P1 = Tp.par()[1];
double PV0 = Tp.cov()[0];
@@ -388,14 +524,18 @@ int InDet::SiTrajectoryElement_xk::searchClustersWithoutStereoAssSCT
if(x < Xm) {
InDet::SiClusterLink_xk l(c,x);
- for(int i=0; i!=nl; ++i) L[i].Comparison(l);
- if(nl<10) L[nl++]=l; else Xm=L[9].xi2();
+ /// This misleadingly named method actually is a sorting
+ /// Mechanism: Our current link will take the place of the one
+ /// it replaces in ascending chi², and we will be
+ /// left with a link to the worst chi² one
+ for(int i=0; i!=nLinksFound; ++i) L[i].Comparison(l);
+ if(nLinksFound<10) L[nLinksFound++]=l; else Xm=L[9].xi2();
Xc = Xm;
}
- else if(!nl) {Xc = x; cl = c;}
+ else if(!nLinksFound) {Xc = x; cl = c;}
}
- if(cl && !nl) {L[nl++].Set(cl,Xc);}
- return nl;
+ if(cl && !nLinksFound) {L[nLinksFound++].Set(cl,Xc);}
+ return nLinksFound;
}
///////////////////////////////////////////////////////////////////
@@ -422,7 +562,7 @@ void InDet::SiTrajectoryElement_xk::set
if(m_tools->fieldTool().magneticFieldMode()!=0) m_fieldMode = true;
m_status = 0 ;
m_detstatus = st ;
- m_ndist = 0 ;
+ m_nMissing = 0 ;
m_nlinksForward = 0 ;
m_nlinksBackward = 0 ;
m_nholesForward = 0 ;
diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectoryElement_xk.cxx b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectoryElement_xk.cxx
index 2c95135dd7ea5f165a74057e08b0e3386de42e43..5b20bf15667ce0c843cd98c2dca8bb011a93d99a 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectoryElement_xk.cxx
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectoryElement_xk.cxx
@@ -103,7 +103,7 @@ bool InDet::SiTrajectoryElement_xk::firstTrajectorElement
m_xi2totalForward = 0. ;
m_status = 1 ;
m_inside = -1 ;
- m_ndist = 0 ;
+ m_nMissing = 0 ;
m_nlinksForward = 0 ;
m_nholesForward = 0 ;
m_dholesForward = 0 ;
@@ -133,7 +133,7 @@ bool InDet::SiTrajectoryElement_xk::firstTrajectorElement()
m_xi2totalForward = 0. ;
m_status = 1 ;
m_inside = -1 ;
- m_ndist = 0 ;
+ m_nMissing = 0 ;
m_nlinksForward = 0 ;
m_nholesForward = 0 ;
m_dholesForward = 0 ;
@@ -159,7 +159,7 @@ bool InDet::SiTrajectoryElement_xk::lastTrajectorElement()
m_status = 3 ;
m_inside = -1 ;
- m_ndist = 0 ;
+ m_nMissing = 0 ;
m_nlinksBackward = 0 ;
m_nholesBackward = 0 ;
m_dholesBackward = 0 ;
@@ -206,7 +206,7 @@ bool InDet::SiTrajectoryElement_xk::ForwardPropagationWithoutSearch
/// copy information from starting TE
m_nclustersForward = TE.m_nclustersForward;
m_nholesForward = TE.m_nholesForward ;
- m_ndist = TE.m_ndist ;
+ m_nMissing = TE.m_nMissing ;
m_ndfForward = TE.m_ndfForward ;
m_xi2totalForward = TE.m_xi2totalForward ;
m_step += TE.m_step ;
@@ -237,7 +237,7 @@ bool InDet::SiTrajectoryElement_xk::ForwardPropagationWithoutSearch
++m_dholesForward;
}
/// for bond gaps, increment ndist
- if(m_dist < -2.) ++m_ndist;
+ if(m_dist < -2.) ++m_nMissing;
}
}
@@ -267,13 +267,15 @@ bool InDet::SiTrajectoryElement_xk::ForwardPropagationWithoutSearch
bool InDet::SiTrajectoryElement_xk::ForwardPropagationWithSearch
(InDet::SiTrajectoryElement_xk& TE)
{
- // Track propagation
- //
+ /// Track propagation
+ /// as usual, if the previous element has a cluster, propagate the updated parameters.
+ /// otherwise the predicted ones.
if(TE.m_cluster) {
TE.m_parametersUpdatedForward.addNoise (TE.m_noise,Trk::alongMomentum);
if(!propagate(TE.m_parametersUpdatedForward,m_parametersPredForward,m_step)) return false;
TE.m_parametersUpdatedForward.removeNoise(TE.m_noise,Trk::alongMomentum);
+ /// double hole running counter is reset if the prev element has a cluster
m_dholesForward = 0;
}
else {
@@ -284,47 +286,73 @@ bool InDet::SiTrajectoryElement_xk::ForwardPropagationWithSearch
m_dholesForward = TE.m_dholesForward;
}
- m_status = 1 ;
+ /// reset old fwd propagation info
+ m_status = 1 ; /// note the socially distant semicolons!
m_nlinksForward = 0 ;
- m_clusterOld = m_cluster ;
- m_cluster = 0 ;
- m_clusterNoAdd = 0 ;
+ /// current cluster --> old cluster, reset current
+ m_clusterOld = m_cluster ;
+ m_cluster = nullptr ;
+ m_clusterNoAdd = nullptr ;
m_xi2Forward = 10000. ;
+ /// re-evaluate if we are expected to intersect this element based on the updated prediction
m_inside = m_detlink->intersect(m_parametersPredForward,m_dist);
+ /// copy running counts from previous element
m_nholesForward = TE.m_nholesForward ;
- m_ndist = TE.m_ndist ;
+ m_nMissing = TE.m_nMissing ;
m_nclustersForward = TE.m_nclustersForward ;
m_ndfForward = TE.m_ndfForward ;
m_xi2totalForward = TE.m_xi2totalForward ;
- m_step += TE.m_step ;
+ m_step += TE.m_step ;
- if(m_inside > 0) {noiseInitiate(); return true;}
+ /// if we are not passing through this element, reset noise production and return
+ if(m_inside > 0) {
+ noiseInitiate(); return true;
+ }
- if((m_nlinksForward=searchClusters(m_parametersPredForward,m_linkForward))) {
-
+ /// now perform cluster search.
+ m_nlinksForward=searchClusters(m_parametersPredForward,m_linkForward);
+ /// if we found any:
+ if(m_nlinksForward!=0) {
+ /// set chi2 to the one for the first candidate
m_xi2Forward = m_linkForward[0].xi2();
+ /// if the chi2 is acceptable:
if (m_xi2Forward <= m_xi2max ) {
-
+ /// use first cluster as cluster on this element
m_cluster = m_linkForward[0].cluster();
+ /// try to add it to the trajectory, update our parameters
if(!addCluster(m_parametersPredForward,m_parametersUpdatedForward)) return false;
+ /// update noise based on this cluster
noiseProduction(1,m_parametersUpdatedForward);
- ++m_nclustersForward; m_xi2totalForward+=m_xi2Forward; m_ndfForward+=m_ndf;
+ /// increment running cluster count for forward trajectory
+ ++m_nclustersForward;
+ /// increment total chi2 for forward trajectory
+ m_xi2totalForward+=m_xi2Forward;
+ /// and the degrees of freedom as well
+ m_ndfForward+=m_ndf;
}
+ /// if we have an outlier:
else if(m_xi2Forward <= m_xi2maxNoAdd) {
-
+ /// update outlier cluster member
m_clusterNoAdd = m_linkForward[0].cluster();
+ /// and reset noise
noiseProduction(1,m_parametersPredForward);
}
- }
+ } /// end of branch for found clusters
else {
+ /// reset noise if nothing was found
noiseProduction(1,m_parametersPredForward);
}
-
+ /// if we are in an active module (with nonzero clusters) and didn't find anything:
if(m_detstatus >=0 && !m_cluster) {
- if(m_inside < 0 && !m_clusterNoAdd) {++m_nholesForward; ++m_dholesForward;}
- if(m_dist < -2. ) ++m_ndist;
- }
+ /// if we expect a hit, and didn't see an outlier: increment hole counts
+ if(m_inside < 0 && !m_clusterNoAdd) {
+ ++m_nholesForward;
+ ++m_dholesForward;
+ }
+ /// if we are well within bounds, increment nDist (may be bond-gap case)
+ if(m_dist < -2. ) ++m_nMissing;
+ } /// end of case covering no found cluster
return true;
}
@@ -375,7 +403,7 @@ bool InDet::SiTrajectoryElement_xk::BackwardPropagationFilter
m_xi2Backward = 10000.;
m_inside = m_detlink->intersect(m_parametersPredBackward,m_dist);
m_nholesBackward = TE.m_nholesBackward ;
- m_ndist = TE.m_ndist ;
+ m_nMissing = TE.m_nMissing ;
m_nclustersBackward = TE.m_nclustersBackward ;
m_npixelsBackward = TE.m_npixelsBackward ;
m_ndfBackward = TE.m_ndfBackward ;
@@ -407,7 +435,7 @@ bool InDet::SiTrajectoryElement_xk::BackwardPropagationFilter
if(m_detstatus >=0 && !m_cluster){
if(m_inside < 0 && !m_clusterNoAdd) {++m_nholesBackward; ++m_dholesBackward;}
- if(m_dist < -2.) ++m_ndist;
+ if(m_dist < -2.) ++m_nMissing;
}
return true;
}
@@ -450,7 +478,7 @@ bool InDet::SiTrajectoryElement_xk::BackwardPropagationSmoother
m_clusterNoAdd = 0 ;
m_xi2Backward = 10000. ;
m_nholesBackward = TE.m_nholesBackward ;
- m_ndist = TE.m_ndist ;
+ m_nMissing = TE.m_nMissing ;
m_nclustersBackward = TE.m_nclustersBackward ;
m_npixelsBackward = TE.m_npixelsBackward ;
m_ndfBackward = TE.m_ndfBackward ;
@@ -482,7 +510,7 @@ bool InDet::SiTrajectoryElement_xk::BackwardPropagationSmoother
}
if(m_detstatus >=0 && !m_cluster) {
if(m_inside < 0 && !m_clusterNoAdd) {++m_nholesBackward; ++m_dholesBackward;}
- if(m_dist < -2.) ++m_ndist;
+ if(m_dist < -2.) ++m_nMissing;
}
return true;
}
@@ -498,7 +526,7 @@ bool InDet::SiTrajectoryElement_xk::BackwardPropagationSmoother
if(!m_cluster) {
if(m_inside < 0 && !m_clusterNoAdd) {++m_nholesBackward; ++m_dholesBackward;}
- if(m_dist < -2.) ++m_ndist;
+ if(m_dist < -2.) ++m_nMissing;
}
return true;
@@ -529,7 +557,7 @@ bool InDet::SiTrajectoryElement_xk::addNextClusterB()
m_nlinksBackward = 0;
m_cluster = 0; --m_nclustersBackward; m_ndfBackward-=m_ndf; if(m_ndf==2) --m_npixelsBackward;
if(m_inside < 0 ) {++m_nholesBackward; ++m_dholesBackward;} m_xi2Backward = 0.;
- if(m_dist < -2.) ++m_ndist;
+ if(m_dist < -2.) ++m_nMissing;
}
return true;
}
@@ -559,7 +587,7 @@ bool InDet::SiTrajectoryElement_xk::addNextClusterF()
m_nlinksForward = 0;
m_cluster = 0; --m_nclustersForward; m_ndfForward-=m_ndf;
if(m_inside < 0) {++m_nholesForward; ++m_dholesForward;} m_xi2Forward = 0.;
- if(m_dist < -2.) ++m_ndist;
+ if(m_dist < -2.) ++m_nMissing;
}
return true;
}
@@ -578,7 +606,7 @@ bool InDet::SiTrajectoryElement_xk::addNextClusterB
m_npixelsBackward = TE.m_npixelsBackward ;
m_ndfBackward = TE.m_ndfBackward ;
m_nholesBackward = TE.m_nholesBackward ;
- m_ndist = TE.m_ndist ;
+ m_nMissing = TE.m_nMissing ;
m_xi2totalBackward = TE.m_xi2totalBackward ;
TE.m_cluster ? m_dholesBackward = 0 : m_dholesBackward = TE.m_dholesBackward;
@@ -592,7 +620,7 @@ bool InDet::SiTrajectoryElement_xk::addNextClusterB
else {
if(m_inside < 0) {++m_nholesBackward; ++m_dholesBackward;} m_xi2Backward = 0.;
- if(m_dist < -2.) ++m_ndist;
+ if(m_dist < -2.) ++m_nMissing;
}
return true;
}
@@ -612,7 +640,7 @@ bool InDet::SiTrajectoryElement_xk::addNextClusterF
m_nclustersForward = TE.m_nclustersForward;
m_ndfForward = TE.m_ndfForward ;
m_nholesForward = TE.m_nholesForward ;
- m_ndist = TE.m_ndist ;
+ m_nMissing = TE.m_nMissing ;
m_xi2totalForward = TE.m_xi2totalForward ;
/// if the previous element has a hit, reset incremental
/// hole counter
@@ -645,8 +673,8 @@ bool InDet::SiTrajectoryElement_xk::addNextClusterF
}
/// no chi2 contribution here
m_xi2Forward = 0.;
- /// if bond gap, increment ndist
- if(m_dist < -2.) ++m_ndist;
+ /// if crossed
+ if(m_dist < -2.) ++m_nMissing;
}
return true;
}
@@ -1574,7 +1602,7 @@ InDet::SiTrajectoryElement_xk::SiTrajectoryElement_xk()
m_status = 0 ;
m_nlinksForward = 0 ;
m_nlinksBackward = 0 ;
- m_ndist = 0 ;
+ m_nMissing = 0 ;
m_radlength = .03;
m_radlengthN = .03;
m_energylose = .4 ;
@@ -1636,7 +1664,7 @@ InDet::SiTrajectoryElement_xk& InDet::SiTrajectoryElement_xk::operator =
m_status = E.m_status ;
m_detstatus = E.m_detstatus ;
m_inside = E.m_inside ;
- m_ndist = E.m_ndist ;
+ m_nMissing = E.m_nMissing ;
m_stereo = E.m_stereo ;
m_detelement = E.m_detelement ;
m_detlink = E.m_detlink ;
@@ -1748,21 +1776,35 @@ bool InDet::SiTrajectoryElement_xk::isNextClusterHoleB(bool& cl,double& X)
return true;
}
+/// checks if removing this cluster from the forward propagation would
+/// result in a critical number of holes or double holes.
+/// Fills "cl" with true if we would still find another cluster after removing
+/// the preferred one. Fills X with the chi square we would have on removal.
+/// Return true if we would *not* get a problematic number of holes
bool InDet::SiTrajectoryElement_xk::isNextClusterHoleF(bool& cl,double& X)
{
cl = false ;
+ /// chi2 if we drop this cluster
X = m_xi2totalForward-m_xi2Forward;
-
+ /// not used in SiTrajectory?
if(m_detstatus == 2) return false;
+ /// if we have another good cluster on this element which could replace
+ /// the current one were it dropped:
if(m_nlinksForward > 1 && m_linkForward[1].xi2() <= m_xi2max) {
+ /// add the chi2 of the second forward link to the second arg
X+=m_linkForward[1].xi2();
- cl = true; return true;
+ /// set the return-arg signifying we would still see a cluster and return
+ cl = true;
+ return true;
}
-
+ /// if we expect a hit:
if(m_inside < 0) {
+ /// return true if we are below the max allowed holes (so we can add one)
if(m_nholesForward < m_maxholes && m_dholesForward < m_maxdholes) return true;
+ /// otherwise return false
return false;
}
+ /// if we do not expect a hit return true
return true;
}
@@ -1784,7 +1826,7 @@ void InDet::SiTrajectoryElement_xk::setParametersF(Trk::PatternTrackParameters&
void InDet::SiTrajectoryElement_xk::setNdist(int n)
{
- m_ndist = n;
+ m_nMissing = n;
}
/////////////////////////////////////////////////////////////////////////////////
@@ -1926,6 +1968,7 @@ void InDet::SiTrajectoryElement_xk::bremNoiseModel()
///////////////////////////////////////////////////////////////////
int InDet::SiTrajectoryElement_xk::searchClusters(Trk::PatternTrackParameters& Tp, SiClusterLink_xk* L) {
+ /// delegate work to subsystem-specific template implementation
if (m_itType==PixelClusterColl) return searchClustersSub<InDet::PixelClusterCollection::const_iterator>(Tp, L);
if (m_itType==SCT_ClusterColl) return searchClustersSub<InDet::SCT_ClusterCollection::const_iterator>(Tp, L);
return searchClustersSub<InDet::SiClusterCollection::const_iterator>(Tp, L);
diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectory_xk.cxx b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectory_xk.cxx
index 12f281ce1f09a021b5b5219a85b476bc1254025b..d4fb95a6b74ce4673213694a132a952f7a4d2e1c 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectory_xk.cxx
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectory_xk.cxx
@@ -1074,7 +1074,7 @@ bool InDet::SiTrajectory_xk::backwardExtension(int itmax)
for (; it!=itmax; ++it) {
int l = F;
- int p = F;
+ int lastElementWithExpHit = F;
for (--F; F>=0; --F) {
@@ -1084,10 +1084,10 @@ bool InDet::SiTrajectory_xk::backwardExtension(int itmax)
if (!Ef.BackwardPropagationFilter(El)) break;
if (Ef.cluster()) {
- p=F; l=F;
+ lastElementWithExpHit=F; l=F;
}
else if (Ef.inside() < 0 ) {
- p=F; if (Ef.nholesB() > maxholes || Ef.dholesB() > maxdholes) break;
+ lastElementWithExpHit=F; if (Ef.nholesB() > maxholes || Ef.dholesB() > maxdholes) break;
}
int nm = Ef.nclustersB()+F;
if (Ef.ndist() > ndcut || nm < nclbest || (nm == nclbest && Ef.xi2totalB() > Xi2best) ) break;
@@ -1112,7 +1112,7 @@ bool InDet::SiTrajectory_xk::backwardExtension(int itmax)
nbest = 0 ;
ndfbest = 0 ;
hbest = nh ;
- hbestb = m_elements[m_elementsMap[p]].nholesB()-nh ;
+ hbestb = m_elements[m_elementsMap[lastElementWithExpHit]].nholesB()-nh ;
itbest = it ;
Xi2best = X ;
PA = m_elements[m_elementsMap[l]].parametersUB();
@@ -1231,28 +1231,29 @@ bool InDet::SiTrajectory_xk::forwardExtension(bool smoother,int itmax)
if (m_firstElement >= m_lastElement) return false;
- /// last element with a known hit
- int L = m_lastElement;
+ /// index at which we start the extension
+ int extensionStartIndex = m_lastElement;
/// last element on the trajectory
int lastElementOnTraj = m_nElements-1;
- /// smoothing step, if requested
+ /// smoothing step, if requested.
+ /// Will re-evaluate the extension start and parameters there
if (smoother) {
- L = m_firstElement;
+ extensionStartIndex = m_firstElement;
/// This checks if we already ran a forward and backward propagation within the current
/// hits on the trajectory. Will enter this branch if we have not.
- if (m_elements[m_elementsMap[L]].difference()) {
+ if (m_elements[m_elementsMap[extensionStartIndex]].difference()) {
/// in this case, we will simply set up for forward propagation, projecting the predicted
/// state forward without searching for new hits yet
/// prepare first element for forward propagation
- if (!m_elements[m_elementsMap[L]].firstTrajectorElement()) return false;
+ if (!m_elements[m_elementsMap[extensionStartIndex]].firstTrajectorElement()) return false;
/// for all following elements with until the last with a known hit:
- for (++L; L<=m_lastElement; ++L) {
- InDet::SiTrajectoryElement_xk& previousElement = m_elements[m_elementsMap[L-1]];
- InDet::SiTrajectoryElement_xk& thisElement = m_elements[m_elementsMap[L ]];
+ for (++extensionStartIndex; extensionStartIndex<=m_lastElement; ++extensionStartIndex) {
+ InDet::SiTrajectoryElement_xk& previousElement = m_elements[m_elementsMap[extensionStartIndex-1]];
+ InDet::SiTrajectoryElement_xk& thisElement = m_elements[m_elementsMap[extensionStartIndex ]];
/// propagate forward the state from the previous element
if (!thisElement.ForwardPropagationWithoutSearch(previousElement)) return false;
}
@@ -1262,10 +1263,10 @@ bool InDet::SiTrajectory_xk::forwardExtension(bool smoother,int itmax)
else{
bool diff = false;
/// for all elements starting from the one after our first hit, up to the last one with a hit:
- for (++L; L<=m_lastElement; ++L) {
+ for (++extensionStartIndex; extensionStartIndex<=m_lastElement; ++extensionStartIndex) {
- InDet::SiTrajectoryElement_xk& previousElement = m_elements[m_elementsMap[L-1]];
- InDet::SiTrajectoryElement_xk& thisElement = m_elements[m_elementsMap[L ]];
+ InDet::SiTrajectoryElement_xk& previousElement = m_elements[m_elementsMap[extensionStartIndex-1]];
+ InDet::SiTrajectoryElement_xk& thisElement = m_elements[m_elementsMap[extensionStartIndex ]];
/// this branch is entered while the range of hits we are looking at was still covered
/// by the previously made backward smoothing
if (!diff) {
@@ -1282,153 +1283,306 @@ bool InDet::SiTrajectory_xk::forwardExtension(bool smoother,int itmax)
/// propagate forward
if (!thisElement.ForwardPropagationWithoutSearch(previousElement)) return false;
}
- }
- }
- --L;
- }
- if ( L== lastElementOnTraj) return true;
+ } /// end loop from second to last element with hit
+ } /// end case of existing smoothing
+ --extensionStartIndex; /// decrement extensionStartIndex to get back to the last element with a hit on it
+ } /// end smoother
+
+ /// if we are already at the last element on the trajectory, we have no-where to extend further
+ if ( extensionStartIndex== lastElementOnTraj) return true;
// Search best forward trajectory prolongation
- //
+ /// These represent the indices used for the extension. Only contains
+ /// the elements beyond the start (so not the piece already existing).
int MP [300] ;
int MPbest[300] ;
+ /// Trajectory elements of clusters on the best candidate
int TE [100] ;
+ /// Clusters on the best candidate
const InDet::SiCluster* CL [100] ;
+ /// hole and double-hole cuts
int maxholes = m_tools->maxholes () ;
int maxdholes = m_tools->maxdholes() ;
+ /// max iteration we expect to reach
const int itm = itmax-1 ;
- int it = 0 ;
+ /// current iteration
+ int iteration = 0 ;
+ /// index of best iteration so far
int itbest = 0 ;
+ /// best quality seen so far
int qbest =-100 ;
+ /// number of hits on the best extension
int nbest = 0 ;
int ndfbest = 0 ;
+ /// holes for best extension (only up to last hit)
int hbest = 0 ;
+ /// holes beyond last hit for best extension
int hbeste = 0 ;
+ /// number of missing expected hits on best extension (holes + deads)
int ndbest = 0 ;
+ /// max holes plus deads
int ndcut = 3 ;
+ /// best number of clusters
int nclbest = 0 ;
- int lbest = L ;
- int F = L ;
- int M = L ;
- MP [M] = L ;
+ /// start index for best extension
+ int lbest = extensionStartIndex ;
+ /// index for looping over detector elements
+ int index_currentElement = extensionStartIndex ;
+ /// current index in the MP array - runs along all elements on the candidate track
+ int M = extensionStartIndex ;
+ /// first point of current trajectory: start of extension
+ MP [M] = extensionStartIndex ;
double Xi2best = 0. ;
+ /// maximum dPhi
const double dfmax = 2.2 ;
+
+ /// this is the at the location of the first hit on track
double f0 = m_elements[m_elementsMap[m_firstElement]].parametersUF().par()[2];
- m_elements[m_elementsMap[F]].setNdist(0);
- for (; it!=itmax; ++it) {
+ m_elements[m_elementsMap[index_currentElement]].setNdist(0);
- int l = F;
- int p = F;
- int Fs = F;
- int Ml = M;
+ /// start to iterate
+ for (; iteration!=itmax; ++iteration) {
+
+ int lastElementWithSeenHit = index_currentElement;
+ int lastElementWithExpHit = index_currentElement;
+ /// index of the previous detector element in the forward propagation
+ int index_previousElement = index_currentElement;
+ /// Last index in the M array where we found a cluster
+ int mLastCluster = M;
+ /// missing clusters on the best candidate compared to the maximum possible
int Cm = nclbest-lastElementOnTraj;
- bool h = false;
+ /// hole flag
+ bool haveHole = false;
- for (++F; F!=m_nElements; ++F) {
+ /// loop over detector elements - the first time this is called, we start at the first
+ /// one beyond the one that saw the last hit
+ for (++index_currentElement; index_currentElement!=m_nElements; ++index_currentElement) {
- InDet::SiTrajectoryElement_xk& El = m_elements[m_elementsMap[Fs]];
- InDet::SiTrajectoryElement_xk& Ef = m_elements[m_elementsMap[F ]];
-
- if (!Ef.ForwardPropagationWithSearch(El)) {
-
- if (!Ef.isBarrel() || Fs!=F-1) break;
-
- int f = F;
- for (; f!=m_nElements; ++f) {
- if (!m_elements[m_elementsMap[f]].isBarrel() ) break;
+ ///
+ InDet::SiTrajectoryElement_xk& prevElement = m_elements[m_elementsMap[index_previousElement]];
+ InDet::SiTrajectoryElement_xk& currentElement = m_elements[m_elementsMap[index_currentElement ]];
+
+ /// propagate forward to the current element, and search for matching clusters
+ if (!currentElement.ForwardPropagationWithSearch(prevElement)) {
+ /// In case of a forward propagation error, try to recover cases due to barrel-endcap transitions
+
+ /// we are already in the endcaps, or if the incompatible elements the failure are not subsequent: abooooort!
+ if (!currentElement.isBarrel() || index_previousElement!=index_currentElement-1) break;
+
+ /// in the barrel or if the previous element is the one just before this, loop over the remaining elements
+ /// and look for the first one not in the barrel.
+ /// This detects barrel-endcap transitions
+ int index_auxElement = index_currentElement;
+ for (; index_auxElement!=m_nElements; ++index_auxElement) {
+ if (!m_elements[m_elementsMap[index_auxElement]].isBarrel() ) break;
}
- if (f==m_nElements) break;
-
- F = f-1; continue;
- }
- else {Fs = F;} MP[++M] = F;
+ /// we didn't find any endcap elements? Likely a bad attempt, so abort
+ if (index_auxElement==m_nElements) break;
- if (Ef.cluster() ) {
- double df = fabs(Ef.parametersUF().par()[2]-f0); if (df > pi) df = pi2-df;
+ /// If we found a barrel-endcap transition,
+ /// we continue our loop at the first endcap element
+ /// along the trajectory.
+ index_currentElement = index_auxElement-1;
+ continue;
+ } /// end of propagation failure handling
+
+ /// NOTE: The 'else' here is mainly for clarity - due to the 'continue'/'break' statements
+ /// in the 'if' above, we only reach here if the condition above (propagation failure) is not met.
+ else {
+ /// if the forward propagation was succesful:
+ /// update the 'previous element' to the current one
+ index_previousElement = index_currentElement;
+ }
+ /// Reminder: we only reach this point in case of a succesful forward propagation
+
+ /// add the current element to the MP list, increment M index
+ MP[++M] = index_currentElement;
+
+ /// if we found a matching cluster at this element:
+ if (currentElement.cluster() ) {
+ /// get the dPhi w.r.t the first hit
+ double df = fabs(currentElement.parametersUF().par()[2]-f0);
+ /// wrap into 0...pi space
+ if (df > pi) df = pi2-df;
+ /// and bail out if the track is curving too extremely in phi
if (df > dfmax) break;
- p=F; l=F; Ml = M; h=false;
+ /// update indices for last hit
+ lastElementWithExpHit=index_currentElement;
+ lastElementWithSeenHit=index_currentElement;
+ mLastCluster = M;
+ haveHole=false;
}
-
- else if (Ef.inside() < 0 ) {
- p=F; if (Ef.nholesF() > maxholes || Ef.dholesF() > maxdholes) break; h=true;
-
- double df = fabs(Ef.parametersPF().par()[2]-f0); if (df > pi) df = pi2-df;
+ /// we did not find a compatible hit, but we would expect one --> hole!
+ else if (currentElement.inside() < 0 ) {
+ lastElementWithExpHit=index_currentElement;
+ /// check if we exceeded the maximum number of holes or double holes. If yes, abort!
+ if (currentElement.nholesF() > maxholes || currentElement.dholesF() > maxdholes) break;
+ haveHole=true;
+ /// and do the dPhi check again, but using the predicted parameters (as no hit)
+ double df = fabs(currentElement.parametersPF().par()[2]-f0);
+ if (df > pi) df = pi2-df;
if (df > dfmax ) break;
}
+ /// we did not find a hit, but also do not cross the active surface (no hit expected)
else {
- double df = fabs(Ef.parametersPF().par()[2]-f0); if (df > pi) df = pi2-df;
+ /// apply only the dPhi check
+ double df = fabs(currentElement.parametersPF().par()[2]-f0);
+ if (df > pi) df = pi2-df;
if (df > dfmax) break;
}
- int nm = Ef.nclustersF()-F;
- if (Ef.ndist() > ndcut || nm < Cm || (nm == Cm && Ef.xi2totalF() > Xi2best)) break;
- }
-
- int m = m_elementsMap[l];
+ /// number of missing clusters so far (for whatever reason)
+ /// note: negative number (0 == all possible collected)
+ int nm = currentElement.nclustersF()-index_currentElement;
+ /// Now, check a number of reasons to exit early:
+ if ( currentElement.ndist() > ndcut /// if we have too many deads
+ || nm < Cm /// or we have already missed more hits than in the best so far
+ || (nm == Cm && currentElement.xi2totalF() > Xi2best) /// or we have missed the same #hits but have a worse chi2
+ ) break;
+ } /// end of loop over elements
+
+ /// now we have assembled a full candidate for the extended trajectory.
+
+ /// get the index of the element where we last saw a hit
+ int m = m_elementsMap[lastElementWithSeenHit];
+ /// get the number of clusters we saw at that point
int nc = m_elements[m].nclustersF();
+ /// and the number of holes (this means we exclude holes after the last hit!)
int nh = m_elements[m].nholesF();
- m_nHolesAfter = m_elements[m_elementsMap[p]].nholesF()-nh;
-
- if (it==0 && nc==m_nclusters) return true;
-
- int fl = F; if (fl==m_nElements) --fl;
- int nd = m_elements[m_elementsMap[fl]].ndist();
+ /// set the number of holes after the last hit - total holes minus holes until last hit
+ m_nHolesAfter = m_elements[m_elementsMap[lastElementWithExpHit]].nholesF()-nh;
+
+ /// if we are in the first iteration, and have not found any new clusters
+ /// beyond the ones already collected in initialize(), return.
+ if (iteration==0 && nc==m_nclusters) return true;
+
+ /// get the last element we reached before exiting the loop
+ int lastElementProcessed = index_currentElement;
+ /// if we looped over all, make it the last element there is
+ if (lastElementProcessed==m_nElements) --lastElementProcessed;
+
+ /// number of missed expected hits seen
+ int nd = m_elements[m_elementsMap[lastElementProcessed]].ndist();
+ /// clusters minus holes before last hit
+ /// This is used as a quality estimator
int q = nc-nh;
+ /// total chi2
double X = m_elements[m].xi2totalF();
+ /// if the quality is better than the current best (chi2 as tie-breaker):
if ( (q > qbest) || (q==qbest && X < Xi2best ) ) {
+ /// update the quality flags
qbest = q;
+ ///reset the running index of clusters on best track (will now be populated)
nbest = 0;
+ /// and the degree-of-freedom count
ndfbest = 0;
+ /// update best candidate hole counts
hbest = nh;
- hbeste = m_elements[m_elementsMap[p]].nholesF()-nh;
- itbest = it;
- lbest = L ;
+ /// as well as holes after last hit
+ hbeste = m_elements[m_elementsMap[lastElementWithExpHit]].nholesF()-nh;
+ /// store the iteration number...
+ itbest = iteration;
+ /// extensionStartIndex-value used for the best candidate
+ /// (starting point)
+ lbest = extensionStartIndex ;
+ /// the chi2...
Xi2best = X ;
- ndbest = nd; if (fl==lastElementOnTraj && nd < ndcut) ndcut = nd;
-
- for (int j=L+1; j<=Ml; ++j) {
-
+ /// the number of missed expected hits
+ ndbest = nd;
+ /// if our track went through the full trajectory and we collect a small number
+ /// of missed expected hits: update the cut for the future (need to be better)
+ if (lastElementProcessed==lastElementOnTraj && nd < ndcut) ndcut = nd;
+
+ /// now, loop over the elements on this possible extension again, up to the last cluster
+ for (int j=extensionStartIndex+1; j<=mLastCluster; ++j) {
+ /// get the corresponding element
int i = MP[j];
InDet::SiTrajectoryElement_xk& Ei = m_elements[m_elementsMap[i]];
-
+ /// if we are inside or within tolerance
if (Ei.inside() <= 0) {
+ /// then store this index in MPbest (meaning we skip any element that isn't crossed)
MPbest[++lbest] = i;
- if (Ei.cluster()) {CL[nbest]=Ei.cluster(); TE[nbest++]=lbest; ndfbest+=Ei.ndf();}
+ /// if we have a hit here,
+ if (Ei.cluster()) {
+ /// store the cluster and l-index
+ CL[nbest]=Ei.cluster();
+ TE[nbest++]=lbest;
+ /// and increment NDF count
+ ndfbest+=Ei.ndf();
+ }
}
}
+ /// best total cluster count: clusters we had so far + newly found ones
nclbest = m_nclusters+nbest;
- if ( (nclbest >= 14 && !h) || (fl==lastElementOnTraj && ndbest == 0)) break;
- }
-
- F = -1; bool cl = false; int nb = lastElementOnTraj-nclbest-1; double Xn;
-
- for (int j=Ml; j!=L; --j) {
-
+ /// if we have an amazing track with at least 14(!) clusters and no holes, or if we passed along the full trajectory
+ /// without any missed expected hits, we stop iterating, as we can not improve
+ if ( (nclbest >= 14 && !haveHole) || (lastElementProcessed==lastElementOnTraj && ndbest == 0)) break;
+ } /// end of branch for quality better than current best
+
+ /// set index_currentElement negative (use as iteration exit condition below)
+ index_currentElement = -1;
+ /// boolean to check if we have two clusters on a given element
+ bool cl = false;
+ /// helper number - used below
+ int nb = lastElementOnTraj-nclbest-1;
+ /// chi square holder
+ double Xn;
+
+ /// Now, we prepare for the next iteration.
+ /// This is done by walking backward along the found trajectory and seeing if, by not using one cluster,
+ /// we could theoretically come to a better solution if doing so would give us hits everywhere
+ /// else along the remaining
+ /// trajectory. The first such case is used to start the next iteration.
+
+ /// Start reverse loop
+ for (int j=mLastCluster; j!=extensionStartIndex; --j) {
+ /// corresponding element
int i = MP[j];
InDet::SiTrajectoryElement_xk& Ei = m_elements[m_elementsMap[i]];
+ /// if we are not on the last element, have a hit, and removing the hit would not make us fail the cuts:
if (i!=lastElementOnTraj && Ei.cluster() && Ei.isNextClusterHoleF(cl,Xn)) {
+ /// this is Ei.nclustersF() + [steps from lastElementOnTraj to i] - nclbest - 1
+ /// ==> change in number of hits that would be possible w.r.t the best track
+ /// if ALL following elements had hits, if we removed this hit (to get a better extension)
int nm = nb-i+Ei.nclustersF();
- if (!cl) {if (Ei.dist() < -2. && Ei.ndist() > ndcut-1) continue;}
+ /// if removing this hit would result in us losing a hit on track (no backup on hand):
+ if (!cl) {
+ /// if we would gain a missed expected hit by removing this one and fail the missing hit cut due to this,
+ /// keep iterating and don't do anything here
+ if (Ei.dist() < -2. && Ei.ndist() > ndcut-1) continue;
+ }
+ /// if we would keep a hit even if we remove this cluster, increment the counter to account for this
else ++nm;
-
+
+ /// if we can not possibly improve (nm < 0) or if we can not improve and the chi2 can not
+ /// go below the best one, keep looping, look for another cluster to pull out
if (nm < 0 || (nm == 0 && Xn > Xi2best)) continue;
- F=i; M=j; break;
+ /// otherwise, start the next iteration here!
+ index_currentElement=i;
+ M=j;
+ break;
}
- }
+ } /// end of reverse loop
- if (F < 0 ) break;
- if (it!=itm) if (!m_elements[m_elementsMap[F]].addNextClusterF()) break;
+ /// if we did not find any candidate where we could improve by removing a hit, stop iteration
+ if (index_currentElement < 0 ) break;
+ /// if we are not in the last iteration, final preparation for iterating by making our new start
+ /// point skip the cluster used there so far.
+ if (iteration!=itm && !m_elements[m_elementsMap[index_currentElement]].addNextClusterF()) break;
}
- if (it == itmax) --it;
+ /// if reaching max iterations, reset iteration counter to the last iteration that was run
+ if (iteration == itmax) --iteration;
+ /// exit if no good track found
if (!nbest) return true;
+ /// set members to counters of best track
m_nholes = hbest ;
m_nHolesAfter = hbeste ;
m_nclusters += nbest ;
@@ -1436,44 +1590,84 @@ bool InDet::SiTrajectory_xk::forwardExtension(bool smoother,int itmax)
m_lastElement = TE[nbest-1];
m_nElements = m_lastElement+1;
+ /// if the indices do not match, there is an element along the way which we do not traverse on the sensitive area
if (m_lastElement != MPbest[m_lastElement]) {
- for (int n = L+1; n<=m_lastElement; ++n) m_elementsMap[n]=m_elementsMap[MPbest[n]];
+ /// update the element map, removing unused elements
+ for (int n = extensionStartIndex+1; n<=m_lastElement; ++n){
+ m_elementsMap[n]=m_elementsMap[MPbest[n]];
+ }
}
- if (itbest==it) return true;
+ /// if the last iteration was the best, all is good and we return
+ if (itbest==iteration) return true;
+ ///
int mb = 0;
- F = -1;
- for (int n = L+1; n!=m_nElements; ++n) {
+ /// reset index_currentElement again
+ index_currentElement = -1;
+ /// loop over all detector elements beyond the start point
+ for (int n = extensionStartIndex+1; n!=m_nElements; ++n) {
+ /// get the corresponding element
InDet::SiTrajectoryElement_xk& En = m_elements[m_elementsMap[n]];
int m = mb;
+ /// find the index in the hits-on-track array corresponding to this DE
for (; m!=nbest; ++m) if (TE[m]==n) break;
-
+ /// if we found this element in the array:
if (m!=nbest) {
+ /// if we have since switched to a different cluster on this element compared to the one in the best track:
if (CL[m]!=En.cluster()) {
- if (F<0) {F=n; En.addNextClusterF(m_elements[m_elementsMap[n-1]],CL[m]);}
- else { En.setCluster(CL[m]); }
- }
+ /// if this is the time we need to update a cluster:
+ if (index_currentElement<0) {
+ /// update the index
+ index_currentElement=n;
+ /// set the cluster to be used at this element to the one used for the best extension
+ /// and propagate info from previous element
+ En.addNextClusterF(m_elements[m_elementsMap[n-1]],CL[m]);
+ }
+ /// if we already updated one element: Can't use addNextClusterF directly,
+ /// need to clean up later. For now, just update the cluster info itself
+ else {
+ En.setCluster(CL[m]);
+ }
+ } /// end branch for different active cluster on DE
+ /// exit if we have dealt with all of the clusters now
if (++mb == nbest) break;
- }
+ } /// end of case if we found this element in the cluster-on-track array
+ /// case if a detector element is not on best trajectory
else {
+ /// if this element has an active cluster:
if (En.cluster()) {
- if (F<0) {F=n; En.addNextClusterF(m_elements[m_elementsMap[n-1]],0);}
- else { En.setCluster(0); }
+ /// if this is the first we need to update a cluster, can do the full work
+ if (index_currentElement<0) {
+ index_currentElement=n;
+ /// tell this DE to not use a cluster, and propagate info from prev. element
+ En.addNextClusterF(m_elements[m_elementsMap[n-1]],0);
+ }
+ /// if we already updated one element: Can't use addNextClusterF directly,
+ /// need to clean up later. For now, just update the cluster info itself
+ else {
+ /// tell this DE to not use a cluster
+ En.setCluster(0);
+ }
}
}
- }
- if (F < 0 || m_lastElement == F) {
+ } /// end of loop over all DE beyond start point
+
+ /// if we did not have to update any cluster, or if we only had to
+ /// change the last hit: all good, can exit
+ if (index_currentElement < 0 || m_lastElement == index_currentElement) {
return true;
}
- for (++F; F<=m_lastElement; ++F) {
-
- InDet::SiTrajectoryElement_xk& El = m_elements[m_elementsMap[F-1]];
- InDet::SiTrajectoryElement_xk& Ef = m_elements[m_elementsMap[F ]];
- if (!Ef.ForwardPropagationWithoutSearch(El)) return false;
+ /// otherwise, if we had to update a cluster along the way, need one more forward propagation run
+ /// to make sure all the counters and chi2 etc reflect the best track's cluster configuration
+ for (++index_currentElement; index_currentElement<=m_lastElement; ++index_currentElement) {
+ InDet::SiTrajectoryElement_xk& prevElement = m_elements[m_elementsMap[index_currentElement-1]];
+ InDet::SiTrajectoryElement_xk& currentElement = m_elements[m_elementsMap[index_currentElement ]];
+ if (!currentElement.ForwardPropagationWithoutSearch(prevElement)) return false;
}
+ /// now we can finally exit
return true;
}