diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.h b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.h
index a4d8f2b378b295c9347693076456d1c38b1415fa..4b06836f77ae05afcb9bafea75352fa8ffcf1591 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.h
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.h
@@ -48,21 +48,21 @@ namespace InDet{
const int& detstatus () const {return m_detstatus;}
const int& inside () const {return m_inside; }
const int& ndist () const {return m_ndist ; }
- const int& nlinksF () const {return m_nlinksF; }
- const int& nlinksB () const {return m_nlinksB; }
- const int& nholesF () const {return m_nholesF; }
- const int& nholesB () const {return m_nholesB; }
- const int& dholesF () const {return m_dholesF; }
- const int& dholesB () const {return m_dholesB; }
- const int& nclustersF () const {return m_nclustersF;}
- const int& nclustersB () const {return m_nclustersB;}
- const int& npixelsB () const {return m_npixelsB; }
+ const int& nlinksF () const {return m_nlinksForward; }
+ const int& nlinksB () const {return m_nlinksBackward; }
+ const int& nholesF () const {return m_nholesForward; }
+ const int& nholesB () const {return m_nholesBackward; }
+ const int& dholesF () const {return m_dholesForward; }
+ const int& dholesB () const {return m_dholesBackward; }
+ const int& nclustersF () const {return m_nclustersForward;}
+ const int& nclustersB () const {return m_nclustersBackward;}
+ const int& npixelsB () const {return m_npixelsBackward; }
const bool& stereo () const {return m_stereo; }
const int& status () const {return m_status; }
const int& noiseModel () const {return m_noisemodel;}
const int& ndf () const {return m_ndf; }
- const int& ndfF () const {return m_ndfF; }
- const int& ndfB () const {return m_ndfB; }
+ const int& ndfF () const {return m_ndfForward; }
+ const int& ndfB () const {return m_ndfBackward; }
const int& ntsos () const {return m_ntsos; }
const double& step () const {return m_step; }
const double& dist () const {return m_dist; }
@@ -111,22 +111,26 @@ namespace InDet{
void setNdist(int);
int numberClusters() const;
- const double& xi2F () const {return m_xi2F ;}
- const double& xi2B () const {return m_xi2B ;}
- const double& xi2totalF () const {return m_xi2totalF ;}
- const double& xi2totalB () const {return m_xi2totalB ;}
+ const double& xi2F () const {return m_xi2Forward ;}
+ const double& xi2B () const {return m_xi2Backward ;}
+ const double& xi2totalF () const {return m_xi2totalForward ;}
+ const double& xi2totalB () const {return m_xi2totalBackward ;}
const double& radlength () const {return m_radlength ;}
const double& energylose() const {return m_energylose;}
-
- const Trk::PatternTrackParameters& parametersPF() const {return m_parametersPF;}
- const Trk::PatternTrackParameters& parametersUF() const {return m_parametersUF;}
- const Trk::PatternTrackParameters& parametersPB() const {return m_parametersPB;}
- const Trk::PatternTrackParameters& parametersUB() const {return m_parametersUB;}
+ /// track parameters for forward filter / smoother
+ /// predicted
+ const Trk::PatternTrackParameters& parametersPF() const {return m_parametersPredForward;}
+ /// updated
+ const Trk::PatternTrackParameters& parametersUF() const {return m_parametersUpdatedForward;}
+ /// predicted
+ const Trk::PatternTrackParameters& parametersPB() const {return m_parametersPredBackward;}
+ /// observed
+ const Trk::PatternTrackParameters& parametersUB() const {return m_parametersUpdatedBackward;}
const Trk::PatternTrackParameters& parametersSM() const {return m_parametersSM;}
const Trk::Surface* surface() const {return m_surface;}
- const InDet::SiClusterLink_xk& linkF (int i) const {return m_linkF[i];}
- const InDet::SiClusterLink_xk& linkB (int i) const {return m_linkB[i];}
+ const InDet::SiClusterLink_xk& linkF (int i) const {return m_linkForward[i];}
+ const InDet::SiClusterLink_xk& linkB (int i) const {return m_linkBackward[i];}
//@}
///////////////////////////////////////////////////////////////////
@@ -238,12 +242,8 @@ namespace InDet{
int searchClustersWithStereoAss (Trk::PatternTrackParameters&,SiClusterLink_xk*, const Trk::PRDtoTrackMap &);
//@}
- ///////////////////////////////////////////////////////////////////
- /// @name Is difference between forward and backward propagation
- ///////////////////////////////////////////////////////////////////
- //@{
+ /// check for a difference between forward and back propagation
bool difference() const;
- //@}
///////////////////////////////////////////////////////////////////
// @name TrackStateOnSurface production
@@ -265,9 +265,13 @@ namespace InDet{
///////////////////////////////////////////////////////////////////
// @name Initiate state
+ /// inputPars: input parameters.
+ /// outputPars: output parameters.
+ /// Gives outputPars the x and (for pix) the y of the associated cluster, and the corresponding cov elements,
+ /// and copies all other elements else from inputPars
///////////////////////////////////////////////////////////////////
//@{
- bool initiateState(Trk::PatternTrackParameters&,Trk::PatternTrackParameters&);
+ bool initiateState(Trk::PatternTrackParameters& inputPars,Trk::PatternTrackParameters& outputPars);
//@}
///////////////////////////////////////////////////////////////////
@@ -303,34 +307,91 @@ namespace InDet{
/// @name Propagate parameters with covariance
///////////////////////////////////////////////////////////////////
//@{
- bool propagate
- (Trk::PatternTrackParameters &,
- Trk::PatternTrackParameters &,
- double &);
+ /// Will propagate the startingParameters from their reference
+ /// to the surface associated with this element and write the result into
+ /// the second argument. The third argument records the step length traversed.
+ bool propagate(Trk::PatternTrackParameters & startingParameters,
+ Trk::PatternTrackParameters & outParameters,
+ double & step);
//@}
////////////////////////////////////////////////////////////////////
- /// @name Propagate parameters without covariance
+ /// @name Propagate parameters without covariance update.
+ /// Will propagate the startingParameters from their reference
+ /// to the surface associated with this element and write the result into
+ /// the second argument. The third argument records the step length traversed.
///////////////////////////////////////////////////////////////////
//@{
- bool propagateParameters
- (Trk::PatternTrackParameters&,
- Trk::PatternTrackParameters&,
- double &);
+ bool propagateParameters(Trk::PatternTrackParameters& startingParameters,
+ Trk::PatternTrackParameters& outParameters,
+ double & step);
//@}
///////////////////////////////////////////////////////////////////
/// @name Work methods for propagation
///////////////////////////////////////////////////////////////////
//@{
- void transformPlaneToGlobal
- (bool,Trk::PatternTrackParameters&,double*);
- bool transformGlobalToPlane
- (bool,double*,Trk::PatternTrackParameters&,Trk::PatternTrackParameters&);
- bool rungeKuttaToPlane
- (bool,double*);
- bool straightLineStepToPlane
- (bool,double*);
+
+
+ /////////////////////////////////////////////////////////////////////////////////
+ /// Tramsform from plane to global
+ /// Will take the surface and parameters from localParameters
+ /// and populate globalPars with the result.
+ /// globalPars needs to be at least 46 elements long.
+ /// Convention:
+ /// /dL0 /dL1 /dPhi /dThe /dCM
+ /// X ->globalPars[0] dX / globalPars[ 7] globalPars[14] globalPars[21] globalPars[28] globalPars[35]
+ /// Y ->globalPars[1] dY / globalPars[ 8] globalPars[15] globalPars[22] globalPars[29] globalPars[36]
+ /// Z ->globalPars[2] dZ / globalPars[ 9] globalPars[16] globalPars[23] globalPars[30] globalPars[37]
+ /// Ax ->globalPars[3] dAx/ globalPars[10] globalPars[17] globalPars[24] globalPars[31] globalPars[38]
+ /// Ay ->globalPars[4] dAy/ globalPars[11] globalPars[18] globalPars[25] globalPars[32] globalPars[39]
+ /// Az ->globalPars[5] dAz/ globalPars[12] globalPars[19] globalPars[26] globalPars[33] globalPars[40]
+ /// CM ->globalPars[6] dCM/ globalPars[13] globalPars[20] globalPars[27] globalPars[34] globalPars[41]
+ /// + dA/dS (42-44) and step (45)
+ // /////////////////////////////////////////////////////////////////////////////////
+ void transformPlaneToGlobal(bool,
+ Trk::PatternTrackParameters& localParameters,
+ double* globalPars);
+
+
+ /////////////////////////////////////////////////////////////////////////////////
+ /// Tramsform from global to plane surface
+ /// Will take the global parameters in globalPars, the surface from this trajectory
+ /// element and populate outputParameters with the result.
+ /// If covariance update ("updateJacobian") is requested, also transport the covariance
+ /// from "startingParameters" using the jacobian stored in the globalPars and write this
+ /// into outputParameters.
+ /// Convention for globalPars:
+ /// /dL0 /dL1 /dPhi /dThe /dCM
+ /// X ->globalPars[0] dX / globalPars[ 7] globalPars[14] globalPars[21] globalPars[28] globalPars[35]
+ /// Y ->globalPars[1] dY / globalPars[ 8] globalPars[15] globalPars[22] globalPars[29] globalPars[36]
+ /// Z ->globalPars[2] dZ / globalPars[ 9] globalPars[16] globalPars[23] globalPars[30] globalPars[37]
+ /// Ax ->globalPars[3] dAx/ globalPars[10] globalPars[17] globalPars[24] globalPars[31] globalPars[38]
+ /// Ay ->globalPars[4] dAy/ globalPars[11] globalPars[18] globalPars[25] globalPars[32] globalPars[39]
+ /// Az ->globalPars[5] dAz/ globalPars[12] globalPars[19] globalPars[26] globalPars[33] globalPars[40]
+ /// CM ->globalPars[6] dCM/ globalPars[13] globalPars[20] globalPars[27] globalPars[34] globalPars[41]
+ /// + dA/dS (42-44) and step (45)
+ // /////////////////////////////////////////////////////////////////////////////////
+ bool transformGlobalToPlane(bool updateJacobian,
+ double* globalPars ,
+ Trk::PatternTrackParameters& startingParameters,
+ Trk::PatternTrackParameters& outputParameters);
+
+
+ /////////////////////////////////////////////////////////////////////////////////
+ /// Runge Kutta step to plane
+ /// Updates the "globalPars" array, which is also used to
+ /// pass the input.
+ /// This array follows the convention outlined above (46 components)
+ /////////////////////////////////////////////////////////////////////////////////
+ bool rungeKuttaToPlane(bool updateJacobian, double* globalPars);
+ /////////////////////////////////////////////////////////////////////////////////
+ /// Straight line step to plane
+ /// Updates the "globalPars" array, which is also used to
+ /// pass the input.
+ /// This array follows the convention outlined above (46 components)
+ /////////////////////////////////////////////////////////////////////////////////
+ bool straightLineStepToPlane(bool updateJacobian, double* globalPars);
//@}
private:
@@ -350,31 +411,35 @@ namespace InDet{
bool m_utsos[3] ;
bool m_fieldMode ;
bool m_useassoTool = false ;
+ /// status flag. Start as 0.
+ /// set to 1 if set up for forward
+ /// set to 2 if set up for backward
+ /// set to 3 by lastTrajectoryElement or BackwardPropagationSmoother if m_cluster
int m_status ;
int m_detstatus ; //!< 0 (no clusters)
int m_inside ;
int m_ndist ;
- int m_nlinksF ;
- int m_nlinksB ;
- int m_nholesF ;
- int m_nholesB ;
- int m_dholesF ;
- int m_dholesB ;
- int m_nclustersF ;
- int m_nclustersB ;
- int m_npixelsB ;
+ int m_nlinksForward ;
+ int m_nlinksBackward ;
+ int m_nholesForward ;
+ int m_nholesBackward ;
+ int m_dholesForward ;
+ int m_dholesBackward ;
+ int m_nclustersForward ;
+ int m_nclustersBackward ;
+ int m_npixelsBackward ;
int m_noisemodel ;
int m_ndf ;
- int m_ndfF ;
- int m_ndfB ;
+ int m_ndfForward ;
+ int m_ndfBackward ;
int m_ntsos ;
int m_maxholes ;
int m_maxdholes ;
double m_dist ;
- double m_xi2F ;
- double m_xi2B ;
- double m_xi2totalF ;
- double m_xi2totalB ;
+ double m_xi2Forward ;
+ double m_xi2Backward ;
+ double m_xi2totalForward ;
+ double m_xi2totalBackward ;
double m_radlength ;
double m_radlengthN ;
double m_energylose ;
@@ -384,8 +449,8 @@ namespace InDet{
double m_xi2maxNoAdd ;
double m_xi2maxlink ;
double m_xi2multi ;
- double m_Tr[13] ;
- double m_A [ 3] ;
+ double m_localTransform[13] ; /// the transform for this element
+ double m_localDir [ 3] ;
const InDetDD::SiDetectorElement* m_detelement ;
const InDet::SiDetElementBoundaryLink_xk* m_detlink ;
@@ -405,13 +470,22 @@ namespace InDet{
const InDet::SiCluster* m_cluster ;
const InDet::SiCluster* m_clusterOld ;
const InDet::SiCluster* m_clusterNoAdd;
- Trk::PatternTrackParameters m_parametersPF;
- Trk::PatternTrackParameters m_parametersUF;
- Trk::PatternTrackParameters m_parametersPB;
- Trk::PatternTrackParameters m_parametersUB;
+
+ /// Pattern track parameters
+
+ /// For forward filtering / smoothing
+ /// Predicted state, forward
+ Trk::PatternTrackParameters m_parametersPredForward;
+ /// Updated state, forward
+ Trk::PatternTrackParameters m_parametersUpdatedForward;
+ /// For backward filtering / smoothing
+ /// Predicted state, backward
+ Trk::PatternTrackParameters m_parametersPredBackward;
+ /// Updated state, backward
+ Trk::PatternTrackParameters m_parametersUpdatedBackward;
Trk::PatternTrackParameters m_parametersSM;
- InDet::SiClusterLink_xk m_linkF[10] ;
- InDet::SiClusterLink_xk m_linkB[10] ;
+ InDet::SiClusterLink_xk m_linkForward[10] ;
+ InDet::SiClusterLink_xk m_linkBackward[10] ;
Trk::NoiseOnSurface m_noise ;
const InDet::SiTools_xk* m_tools ;
MagField::AtlasFieldCache m_fieldCache;
diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.icc b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.icc
index bf339d646c7cdaa20695dd1feab5bb0385b812fd..bd1cee22895658dae4f1634dffdcdee4a80adc98 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.icc
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectoryElement_xk.icc
@@ -423,31 +423,31 @@ void InDet::SiTrajectoryElement_xk::set
m_status = 0 ;
m_detstatus = st ;
m_ndist = 0 ;
- m_nlinksF = 0 ;
- m_nlinksB = 0 ;
- m_nholesF = 0 ;
- m_nholesB = 0 ;
- m_dholesF = 0 ;
- m_dholesB = 0 ;
- m_nclustersF = 0 ;
- m_nclustersB = 0 ;
- m_npixelsB = 0 ;
- m_ndfF = 0 ;
- m_ndfB = 0 ;
- m_ntsos = 0 ;
- m_detelement = dl->detElement() ;
- m_detlink = dl ;
- m_surface = &m_detelement->surface();
- m_sibegin = sb ;
- m_siend = se ;
- m_cluster = si ;
- m_clusterOld = si ;
- m_clusterNoAdd = 0 ;
- m_stereo = false ;
- m_xi2F = 10000. ;
- m_xi2B = 10000. ;
- m_xi2totalF = 0. ;
- m_xi2totalB = 0. ;
+ m_nlinksForward = 0 ;
+ m_nlinksBackward = 0 ;
+ m_nholesForward = 0 ;
+ m_nholesBackward = 0 ;
+ m_dholesForward = 0 ;
+ m_dholesBackward = 0 ;
+ m_nclustersForward = 0 ;
+ m_nclustersBackward = 0 ;
+ m_npixelsBackward = 0 ;
+ m_ndfForward = 0 ;
+ m_ndfBackward = 0 ;
+ m_ntsos = 0 ;
+ m_detelement = dl->detElement() ;
+ m_detlink = dl ;
+ m_surface = &m_detelement->surface();
+ m_sibegin = sb ;
+ m_siend = se ;
+ m_cluster = si ;
+ m_clusterOld = si ;
+ m_clusterNoAdd = 0 ;
+ m_stereo = false ;
+ m_xi2Forward = 10000. ;
+ m_xi2Backward = 10000. ;
+ m_xi2totalForward = 0. ;
+ m_xi2totalBackward = 0. ;
m_tools->electron() ? m_xi2max = m_tools->xi2maxBrem() : m_xi2max = m_tools->xi2max();
m_halflength = 0. ;
m_detelement->isSCT() ? m_ndf=1 : m_ndf=2;
@@ -472,12 +472,14 @@ void InDet::SiTrajectoryElement_xk::set
const Amg::Transform3D& tr = m_surface->transform();
- m_Tr[ 0] = tr(0,0); m_Tr[ 1]=tr(1,0); m_Tr[ 2]=tr(2,0);
- m_Tr[ 3] = tr(0,1); m_Tr[ 4]=tr(1,1); m_Tr[ 5]=tr(2,1);
- m_Tr[ 6] = tr(0,2); m_Tr[ 7]=tr(1,2); m_Tr[ 8]=tr(2,2);
- m_Tr[ 9] = tr(0,3); m_Tr[10]=tr(1,3); m_Tr[11]=tr(2,3);
- m_Tr[12] = m_Tr[ 9]*m_Tr[ 6]+m_Tr[10]*m_Tr[ 7]+m_Tr[11]*m_Tr[ 8];
- m_A[0] = 1.; m_A[1] = 0.; m_A[2] = 0.;
+ m_localTransform[ 0] = tr(0,0); m_localTransform[ 1]=tr(1,0); m_localTransform[ 2]=tr(2,0);
+ m_localTransform[ 3] = tr(0,1); m_localTransform[ 4]=tr(1,1); m_localTransform[ 5]=tr(2,1);
+ m_localTransform[ 6] = tr(0,2); m_localTransform[ 7]=tr(1,2); m_localTransform[ 8]=tr(2,2);
+ m_localTransform[ 9] = tr(0,3); m_localTransform[10]=tr(1,3); m_localTransform[11]=tr(2,3);
+ m_localTransform[12] = m_localTransform[ 9]*m_localTransform[ 6]+m_localTransform[10]*m_localTransform[ 7]+m_localTransform[11]*m_localTransform[ 8];
+ m_localDir[0] = 1.;
+ m_localDir[1] = 0.;
+ m_localDir[2] = 0.;
return;
}
@@ -517,9 +519,9 @@ bool InDet::SiTrajectoryElement_xk::ForwardPropagationForClusterSeach
if(!n) {
Trk::PatternTrackParameters Tp; if(!Tp.production(&Tpa)) return false;
- if(!propagateParameters(Tp,m_parametersPF,m_step)) return false;
+ if(!propagateParameters(Tp,m_parametersPredForward,m_step)) return false;
- if(!m_parametersPF.iscovariance()) {
+ if(!m_parametersPredForward.iscovariance()) {
double cv[15]={ .02 ,
0., .02,
@@ -527,13 +529,13 @@ bool InDet::SiTrajectoryElement_xk::ForwardPropagationForClusterSeach
0., 0., 0.,.000001,
0., 0., 0., 0.,.000000001};
- m_parametersPF.setCovariance(cv);
+ m_parametersPredForward.setCovariance(cv);
}
}
else {
- if(!propagate(m_parametersPF,m_parametersPF,m_step)) return false;
+ if(!propagate(m_parametersPredForward,m_parametersPredForward,m_step)) return false;
}
- m_nlinksF=searchClusters(m_parametersPF,m_linkF);
+ m_nlinksForward=searchClusters(m_parametersPredForward,m_linkForward);
return true;
}
@@ -570,7 +572,7 @@ void InDet::SiTrajectoryElement_xk::CloseClusterSeach
m_stereo = true : m_stereo = false;
if(m_detstatus && m_ndf == 1) m_halflength = (*sb)->width().z()*.5;
- if(m_detlink->intersect(Tpa,m_dist) > 0 || !searchClusters(Tpa,m_linkF)) return;
- m_cluster = m_linkF[0].cluster();
- m_xi2F = m_linkF[0].xi2 ();
+ if(m_detlink->intersect(Tpa,m_dist) > 0 || !searchClusters(Tpa,m_linkForward)) return;
+ m_cluster = m_linkForward[0].cluster();
+ m_xi2Forward = m_linkForward[0].xi2 ();
}
diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectory_xk.h b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectory_xk.h
index c7dc76f89882e81e03592ac597f74f9f276948e6..3d6d500894dc185fb5428a31b4495f599e90677a 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectory_xk.h
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectory_xk.h
@@ -48,13 +48,13 @@ namespace InDet{
const int& nholes () const {return m_nholes ;}
const int& dholes () const {return m_dholes ;}
- const int& nholesb () const {return m_nholesb ;}
- const int& nholese () const {return m_nholese ;}
+ const int& nHolesBefore () const {return m_nHolesBefore ;}
+ const int& nHolesAfter () const {return m_nHolesAfter ;}
const int& nclusters () const {return m_nclusters ;}
const int& ndf () const {return m_ndf ;}
const int& nclustersNoAdd() const {return m_nclustersNoAdd;}
const int& nElements () const {return m_nElements ;}
- const int& naElements () const {return m_naElements ;}
+ const int& naElements () const {return m_nActiveElements ;}
const int& difference () const {return m_difference ;}
const int& elementsMap(int& i) const {return m_elementsMap[i];}
@@ -120,24 +120,28 @@ namespace InDet{
// Protected Data
///////////////////////////////////////////////////////////////////
- int m_firstElement ;
- int m_lastElement ;
- int m_nclusters ; // (NCL)
+ int m_firstElement ; /// index of the first element where we have
+ /// a cluster
+ int m_lastElement ; /// index of the last element where we have
+ /// a cluster
+ int m_nclusters ; /// Number of clusters on trajectory
int m_nclustersNoAdd ; // (NCL)
int m_difference ; // forward-bacward diff
- int m_nholesb ; // holes before
- int m_nholese ; // holes after
+ int m_nHolesBefore ; // holes before
+ int m_nHolesAfter ; // holes after
int m_nholes ; // holes
int m_dholes ; // dholes
- int m_naElements ; //
- int m_nElements ; // nindex
+ int m_nActiveElements ; /// count active elements
+ int m_nElements ; // index
int m_elementsMap[300]; // index
int m_ndfcut ; //
int m_ndf ; //
int m_ntos ; //
int m_atos[100] ; //
int m_itos[100] ; //
- SiTrajectoryElement_xk m_elements [300]; // mC
+ SiTrajectoryElement_xk m_elements [300]; /// Trajectory elements on this trajectory.
+ /// Each one corresponds to one detector element on
+ /// the search road
const InDet::SiTools_xk* m_tools ; //
///////////////////////////////////////////////////////////////////
diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectory_xk.icc b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectory_xk.icc
index b70bd1017f528fb4c6987d4286620847d2193728..d22a87102c40524864b79bbe0e8319b3d2975c8d 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectory_xk.icc
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiTrajectory_xk.icc
@@ -17,11 +17,11 @@ inline InDet::SiTrajectory_xk::SiTrajectory_xk()
m_nclusters = 0 ;
m_nclustersNoAdd = 0 ;
m_difference = 0 ;
- m_nholesb = 0 ;
- m_nholese = 0 ;
+ m_nHolesBefore = 0 ;
+ m_nHolesAfter = 0 ;
m_nholes = 0 ;
m_dholes = 0 ;
- m_naElements = 0 ;
+ m_nActiveElements = 0 ;
m_ndfcut = 0 ;
m_ndf = 0 ;
m_ntos = 0 ;
@@ -42,11 +42,11 @@ inline InDet::SiTrajectory_xk& InDet::SiTrajectory_xk::operator =
m_ndf = T.m_ndf ;
m_ntos = T.m_ntos ;
m_nclustersNoAdd = T.m_nclustersNoAdd ;
- m_nholesb = T.m_nholesb ;
- m_nholese = T.m_nholese ;
+ m_nHolesBefore = T.m_nHolesBefore ;
+ m_nHolesAfter = T.m_nHolesAfter ;
m_nholes = T.m_nholes ;
m_dholes = T.m_dholes ;
- m_naElements = T.m_naElements ;
+ m_nActiveElements = T.m_nActiveElements ;
m_nElements = T.m_nElements ;
m_tools = T.m_tools ;
diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectoryElement_xk.cxx b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectoryElement_xk.cxx
index 1235156547776a795a607d45ba72839970fafe38..2c95135dd7ea5f165a74057e08b0e3386de42e43 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectoryElement_xk.cxx
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectoryElement_xk.cxx
@@ -51,52 +51,65 @@ void InDet::SiTrajectoryElement_xk::setParameters()
///////////////////////////////////////////////////////////////////
bool InDet::SiTrajectoryElement_xk::firstTrajectorElement
-(const Trk::TrackParameters& Tpa)
+(const Trk::TrackParameters& startingParameters)
{
+ /// if we don't have a cluster, something went wrong!
if(!m_cluster) return false;
- Trk::PatternTrackParameters Tp; if(!Tp.production(&Tpa)) return false;
+ /// generate pattern track parameters from the starting parameters
+ Trk::PatternTrackParameters startingPatternPars;
+ if(!startingPatternPars.production(&startingParameters)) return false;
- const Trk::Surface* pl = Tp.associatedSurface();
+ /// get the surface belonging to our staring pars
+ const Trk::Surface* pl = startingPatternPars.associatedSurface();
- // Track propagation if need
- //
- if(m_surface==pl) m_parametersPF = Tp;
- else if(!propagate(Tp,m_parametersPF,m_step)) return false;
+ /// Track propagation if needed
+ /// if we are already on the correct surface, we can assign the params as they ar
+ if(m_surface==pl) m_parametersPredForward = startingPatternPars;
+ /// otherwise, we need to propagate to "our" surface
+ else if(!propagate(startingPatternPars,m_parametersPredForward,m_step)) return false;
// Initiate track parameters without initial covariance
//
- double cv[15]={ 1. ,
+ double cv[15]={
+ 1. ,
0. , 1.,
0. , 0.,.001,
0. , 0., 0.,.001,
0. , 0., 0., 0.,.00001};
+ /// update starting cov for DBM case
if(m_detelement->isDBM()) {
- double tn = tan(Tp.par()[3]);
+ double tn = tan(startingPatternPars.par()[3]);
cv[ 5] = .0001 ;
cv[14] = (tn*tn*1.e-6) ;
}
-
- m_parametersPF.setCovariance(cv);
- initiateState(m_parametersPF,m_parametersUF);
-
- noiseProduction(1,m_parametersUF);
-
- m_dist = -10. ;
- m_step = 0. ;
- m_xi2F = 0. ;
- m_xi2totalF = 0. ;
- m_status = 1 ;
- m_inside = -1 ;
- m_ndist = 0 ;
- m_nlinksF = 0 ;
- m_nholesF = 0 ;
- m_dholesF = 0 ;
- m_clusterNoAdd = 0 ;
- m_nclustersF = 1 ;
- m_ndfF = m_ndf;
+ /// write our covariance into the parameters
+ m_parametersPredForward.setCovariance(cv);
+ /// and initiate our state. parametersUF will be constrained
+ /// to the local cluster coordinates (with its covariance),
+ /// while taking the momentum and direction from the input
+ /// parameters
+ initiateState(m_parametersPredForward,m_parametersUpdatedForward);
+
+ /// add noise to the parameter estimate (scattering, eloss)
+ noiseProduction(1,m_parametersUpdatedForward);
+
+ /// and init other values
+ m_dist = -10. ;
+ m_step = 0. ;
+ m_xi2Forward = 0. ;
+ m_xi2totalForward = 0. ;
+ m_status = 1 ;
+ m_inside = -1 ;
+ m_ndist = 0 ;
+ m_nlinksForward = 0 ;
+ m_nholesForward = 0 ;
+ m_dholesForward = 0 ;
+ m_clusterNoAdd = 0 ;
+ m_nclustersForward = 1 ;
+ m_ndfForward = m_ndf;
return true;
}
@@ -108,24 +121,25 @@ bool InDet::SiTrajectoryElement_xk::firstTrajectorElement()
{
if(!m_cluster || !m_status) return false;
- if(m_status > 1 ) m_parametersPF = m_parametersUB;
-
- m_parametersPF.diagonalization(100.);
- if(!initiateState(m_parametersPF,m_parametersUF)) return false;
- noiseProduction(1,m_parametersUF);
-
- m_dist = -10. ;
- m_xi2F = 0. ;
- m_xi2totalF = 0. ;
- m_status = 1 ;
- m_inside = -1 ;
- m_ndist = 0 ;
- m_nlinksF = 0 ;
- m_nholesF = 0 ;
- m_dholesF = 0 ;
- m_clusterNoAdd = 0 ;
- m_nclustersF = 1 ;
- m_ndfF = m_ndf;
+ if(m_status > 1 ) m_parametersPredForward = m_parametersUpdatedBackward;
+
+ /// reset the cov (off-diagonal to zero, diagonal multiplied by 100)
+ m_parametersPredForward.diagonalization(100.);
+ if(!initiateState(m_parametersPredForward,m_parametersUpdatedForward)) return false;
+ noiseProduction(1,m_parametersUpdatedForward);
+
+ m_dist = -10. ;
+ m_xi2Forward = 0. ;
+ m_xi2totalForward = 0. ;
+ m_status = 1 ;
+ m_inside = -1 ;
+ m_ndist = 0 ;
+ m_nlinksForward = 0 ;
+ m_nholesForward = 0 ;
+ m_dholesForward = 0 ;
+ m_clusterNoAdd = 0 ;
+ m_nclustersForward = 1 ;
+ m_ndfForward = m_ndf;
return true;
}
@@ -136,25 +150,25 @@ bool InDet::SiTrajectoryElement_xk::firstTrajectorElement()
bool InDet::SiTrajectoryElement_xk::lastTrajectorElement()
{
if(m_status==0 || !m_cluster) return false;
- noiseProduction(1,m_parametersUF);
-
- m_parametersSM = m_parametersPF;
- m_parametersPB = m_parametersUF;
- m_parametersPB.diagonalization(100.);
- if(!initiateState(m_parametersPB,m_parametersUB)) return false;
-
- m_status = 3 ;
- m_inside = -1 ;
- m_ndist = 0 ;
- m_nlinksB = 0 ;
- m_nholesB = 0 ;
- m_dholesB = 0 ;
- m_clusterNoAdd = 0 ;
- m_nclustersB = 1 ; m_ndf == 2 ? m_npixelsB = 1 : m_npixelsB = 0;
- m_ndfB = m_ndf ;
- m_xi2B = m_xi2F;
- m_xi2totalB = m_xi2F;
- m_dist = -10.;
+ noiseProduction(1,m_parametersUpdatedForward);
+
+ m_parametersSM = m_parametersPredForward;
+ m_parametersPredBackward = m_parametersUpdatedForward;
+ m_parametersPredBackward.diagonalization(100.);
+ if(!initiateState(m_parametersPredBackward,m_parametersUpdatedBackward)) return false;
+
+ m_status = 3 ;
+ m_inside = -1 ;
+ m_ndist = 0 ;
+ m_nlinksBackward = 0 ;
+ m_nholesBackward = 0 ;
+ m_dholesBackward = 0 ;
+ m_clusterNoAdd = 0 ;
+ m_nclustersBackward = 1 ; m_ndf == 2 ? m_npixelsBackward = 1 : m_npixelsBackward = 0;
+ m_ndfBackward = m_ndf ;
+ m_xi2Backward = m_xi2Forward;
+ m_xi2totalBackward = m_xi2Forward;
+ m_dist = -10.;
return true;
}
@@ -166,59 +180,81 @@ bool InDet::SiTrajectoryElement_xk::lastTrajectorElement()
bool InDet::SiTrajectoryElement_xk::ForwardPropagationWithoutSearch
(InDet::SiTrajectoryElement_xk& TE)
{
- // Track propagation
- //
+ /// Track propagation
+ /// If the starting trajectory element has a cluster:
if(TE.m_cluster) {
-
- TE.m_parametersUF.addNoise (TE.m_noise,Trk::alongMomentum);
- if(!propagate(TE.m_parametersUF,m_parametersPF,m_step)) return false;
- TE.m_parametersUF.removeNoise(TE.m_noise,Trk::alongMomentum);
- m_dholesF = 0;
+ /// add noise to the track parameters of the starting TE
+ TE.m_parametersUpdatedForward.addNoise (TE.m_noise,Trk::alongMomentum);
+ /// propagate to the current element, plug into m_parametersPredForward
+ if(!propagate(TE.m_parametersUpdatedForward,m_parametersPredForward,m_step)) return false;
+ /// and restore the starting TE again
+ TE.m_parametersUpdatedForward.removeNoise(TE.m_noise,Trk::alongMomentum);
+ /// reset the double hole count to zero - we had a cluster
+ m_dholesForward = 0;
}
else {
-
- TE.m_parametersPF.addNoise (TE.m_noise,Trk::alongMomentum);
- if(!propagate(TE.m_parametersPF,m_parametersPF,m_step)) return false;
- TE.m_parametersPF.removeNoise(TE.m_noise,Trk::alongMomentum);
- m_dholesF = TE.m_dholesF;
- }
-
- m_nclustersF = TE.m_nclustersF;
- m_nholesF = TE.m_nholesF ;
+ /// add noise to the starting pars
+ TE.m_parametersPredForward.addNoise (TE.m_noise,Trk::alongMomentum);
+ /// propagate, plug into m_parametersPredForward
+ if(!propagate(TE.m_parametersPredForward,m_parametersPredForward,m_step)) return false;
+ /// restore starting TE
+ TE.m_parametersPredForward.removeNoise(TE.m_noise,Trk::alongMomentum);
+ /// double hole count is copied from the previous one
+ m_dholesForward = TE.m_dholesForward;
+ }
+
+ /// copy information from starting TE
+ m_nclustersForward = TE.m_nclustersForward;
+ m_nholesForward = TE.m_nholesForward ;
m_ndist = TE.m_ndist ;
- m_ndfF = TE.m_ndfF ;
- m_xi2totalF = TE.m_xi2totalF ;
+ m_ndfForward = TE.m_ndfForward ;
+ m_xi2totalForward = TE.m_xi2totalForward ;
m_step += TE.m_step ;
- // Track update
- //
+ /// Track update
+ /// If we have a cluster on this element
if( m_cluster) {
-
- if(!addCluster(m_parametersPF,m_parametersUF,m_xi2F)) return false;
- m_inside = -1; ++m_nclustersF; m_xi2totalF+=m_xi2F; m_ndfF+=m_ndf;
- }
+ /// add it to refine the parameter estimate
+ if(!addCluster(m_parametersPredForward,m_parametersUpdatedForward,m_xi2Forward)) return false;
+ /// set m_inside to signify this is a hit on track
+ m_inside = -1;
+ /// increment cluster count
+ ++m_nclustersForward;
+ /// increment chi² and degrees of freedom
+ m_xi2totalForward+=m_xi2Forward;
+ m_ndfForward+=m_ndf;
+ }
+ /// if we have no cluster on this element
else {
-
- m_inside = m_detlink->intersect(m_parametersPF,m_dist);
-
+ /// check the intersection of the propagated track with this surface
+ m_inside = m_detlink->intersect(m_parametersPredForward,m_dist);
+ /// if we have hits on this element
if( m_detstatus >=0) {
- if(m_inside < 0 ) {++m_nholesF; ++m_dholesF;}
+ /// and expect to be inside the active area
+ if(m_inside < 0 ) {
+ /// then increment the hole counters
+ ++m_nholesForward;
+ ++m_dholesForward;
+ }
+ /// for bond gaps, increment ndist
if(m_dist < -2.) ++m_ndist;
}
}
- // Noise production
- //
+ /// Noise production
+ /// If the track passes through our element
if(m_inside<=0) {
-
- if(m_cluster) noiseProduction(1,m_parametersUF);
- else noiseProduction(1,m_parametersPF);
+ /// if we have a cluster, update noise based on the updated parameters
+ if(m_cluster) noiseProduction(1,m_parametersUpdatedForward);
+ /// otherwise, update noise based on the predicted parameters
+ else noiseProduction(1,m_parametersPredForward);
}
+ /// otherwise, reset noise
else {
noiseInitiate();
}
m_status = 1;
- m_nlinksF = 0;
+ m_nlinksForward = 0;
m_clusterNoAdd = 0;
return true;
}
@@ -235,58 +271,58 @@ bool InDet::SiTrajectoryElement_xk::ForwardPropagationWithSearch
//
if(TE.m_cluster) {
- TE.m_parametersUF.addNoise (TE.m_noise,Trk::alongMomentum);
- if(!propagate(TE.m_parametersUF,m_parametersPF,m_step)) return false;
- TE.m_parametersUF.removeNoise(TE.m_noise,Trk::alongMomentum);
- m_dholesF = 0;
+ 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);
+ m_dholesForward = 0;
}
else {
- TE.m_parametersPF.addNoise (TE.m_noise,Trk::alongMomentum);
- if(!propagate(TE.m_parametersPF,m_parametersPF,m_step)) return false;
- TE.m_parametersPF.removeNoise(TE.m_noise,Trk::alongMomentum);
- m_dholesF = TE.m_dholesF;
+ TE.m_parametersPredForward.addNoise (TE.m_noise,Trk::alongMomentum);
+ if(!propagate(TE.m_parametersPredForward,m_parametersPredForward,m_step)) return false;
+ TE.m_parametersPredForward.removeNoise(TE.m_noise,Trk::alongMomentum);
+ m_dholesForward = TE.m_dholesForward;
}
m_status = 1 ;
- m_nlinksF = 0 ;
+ m_nlinksForward = 0 ;
m_clusterOld = m_cluster ;
m_cluster = 0 ;
m_clusterNoAdd = 0 ;
- m_xi2F = 10000. ;
- m_inside = m_detlink->intersect(m_parametersPF,m_dist);
- m_nholesF = TE.m_nholesF ;
+ m_xi2Forward = 10000. ;
+ m_inside = m_detlink->intersect(m_parametersPredForward,m_dist);
+ m_nholesForward = TE.m_nholesForward ;
m_ndist = TE.m_ndist ;
- m_nclustersF = TE.m_nclustersF ;
- m_ndfF = TE.m_ndfF ;
- m_xi2totalF = TE.m_xi2totalF ;
+ m_nclustersForward = TE.m_nclustersForward ;
+ m_ndfForward = TE.m_ndfForward ;
+ m_xi2totalForward = TE.m_xi2totalForward ;
m_step += TE.m_step ;
if(m_inside > 0) {noiseInitiate(); return true;}
- if((m_nlinksF=searchClusters(m_parametersPF,m_linkF))) {
+ if((m_nlinksForward=searchClusters(m_parametersPredForward,m_linkForward))) {
- m_xi2F = m_linkF[0].xi2();
+ m_xi2Forward = m_linkForward[0].xi2();
- if (m_xi2F <= m_xi2max ) {
+ if (m_xi2Forward <= m_xi2max ) {
- m_cluster = m_linkF[0].cluster();
- if(!addCluster(m_parametersPF,m_parametersUF)) return false;
- noiseProduction(1,m_parametersUF);
- ++m_nclustersF; m_xi2totalF+=m_xi2F; m_ndfF+=m_ndf;
+ m_cluster = m_linkForward[0].cluster();
+ if(!addCluster(m_parametersPredForward,m_parametersUpdatedForward)) return false;
+ noiseProduction(1,m_parametersUpdatedForward);
+ ++m_nclustersForward; m_xi2totalForward+=m_xi2Forward; m_ndfForward+=m_ndf;
}
- else if(m_xi2F <= m_xi2maxNoAdd) {
+ else if(m_xi2Forward <= m_xi2maxNoAdd) {
- m_clusterNoAdd = m_linkF[0].cluster();
- noiseProduction(1,m_parametersPF);
+ m_clusterNoAdd = m_linkForward[0].cluster();
+ noiseProduction(1,m_parametersPredForward);
}
}
else {
- noiseProduction(1,m_parametersPF);
+ noiseProduction(1,m_parametersPredForward);
}
if(m_detstatus >=0 && !m_cluster) {
- if(m_inside < 0 && !m_clusterNoAdd) {++m_nholesF; ++m_dholesF;}
+ if(m_inside < 0 && !m_clusterNoAdd) {++m_nholesForward; ++m_dholesForward;}
if(m_dist < -2. ) ++m_ndist;
}
return true;
@@ -305,72 +341,72 @@ bool InDet::SiTrajectoryElement_xk::BackwardPropagationFilter
if(TE.m_cluster) {
- TE.m_parametersUB.addNoise (TE.m_noise,Trk::oppositeMomentum);
- if(!propagate(TE.m_parametersUB,m_parametersPB,m_step)) return false;
- TE.m_parametersUB.removeNoise(TE.m_noise,Trk::oppositeMomentum);
- m_dholesB = 0;
+ TE.m_parametersUpdatedBackward.addNoise (TE.m_noise,Trk::oppositeMomentum);
+ if(!propagate(TE.m_parametersUpdatedBackward,m_parametersPredBackward,m_step)) return false;
+ TE.m_parametersUpdatedBackward.removeNoise(TE.m_noise,Trk::oppositeMomentum);
+ m_dholesBackward = 0;
}
else {
- TE.m_parametersPB.addNoise (TE.m_noise,Trk::oppositeMomentum);
- if(!propagate(TE.m_parametersPB,m_parametersPB,m_step)) return false;
- TE.m_parametersPB.removeNoise(TE.m_noise,Trk::oppositeMomentum);
- m_dholesB = TE.m_dholesB;
+ TE.m_parametersPredBackward.addNoise (TE.m_noise,Trk::oppositeMomentum);
+ if(!propagate(TE.m_parametersPredBackward,m_parametersPredBackward,m_step)) return false;
+ TE.m_parametersPredBackward.removeNoise(TE.m_noise,Trk::oppositeMomentum);
+ m_dholesBackward = TE.m_dholesBackward;
}
}
else {
if(TE.m_cluster) {
- if(!propagate(TE.m_parametersUB,m_parametersPB,m_step)) return false;
- m_dholesB = 0;
+ if(!propagate(TE.m_parametersUpdatedBackward,m_parametersPredBackward,m_step)) return false;
+ m_dholesBackward = 0;
}
else {
- if(!propagate(TE.m_parametersPB,m_parametersPB,m_step)) return false;
- m_dholesB = TE.m_dholesB;
+ if(!propagate(TE.m_parametersPredBackward,m_parametersPredBackward,m_step)) return false;
+ m_dholesBackward = TE.m_dholesBackward;
}
}
m_status = 2;
- m_nlinksB = 0;
+ m_nlinksBackward = 0;
m_clusterOld = m_cluster;
m_cluster = 0;
m_clusterNoAdd = 0;
- m_xi2B = 10000.;
- m_inside = m_detlink->intersect(m_parametersPB,m_dist);
- m_nholesB = TE.m_nholesB ;
+ m_xi2Backward = 10000.;
+ m_inside = m_detlink->intersect(m_parametersPredBackward,m_dist);
+ m_nholesBackward = TE.m_nholesBackward ;
m_ndist = TE.m_ndist ;
- m_nclustersB = TE.m_nclustersB ;
- m_npixelsB = TE.m_npixelsB ;
- m_ndfB = TE.m_ndfB ;
- m_xi2totalB = TE.m_xi2totalB ;
+ m_nclustersBackward = TE.m_nclustersBackward ;
+ m_npixelsBackward = TE.m_npixelsBackward ;
+ m_ndfBackward = TE.m_ndfBackward ;
+ m_xi2totalBackward = TE.m_xi2totalBackward ;
m_step += TE.m_step ;
if(m_inside >0 ) {noiseInitiate(); return true;}
- if((m_nlinksB = searchClusters(m_parametersPB,m_linkB))) {
+ if((m_nlinksBackward = searchClusters(m_parametersPredBackward,m_linkBackward))) {
- m_xi2B = m_linkB[0].xi2();
+ m_xi2Backward = m_linkBackward[0].xi2();
- if (m_xi2B <= m_xi2max ) {
+ if (m_xi2Backward <= m_xi2max ) {
- m_cluster = m_linkB[0].cluster();
- if(!addCluster(m_parametersPB,m_parametersUB)) return false;
- noiseProduction(-1,m_parametersUB);
- ++m_nclustersB; m_xi2totalB+=m_xi2B; m_ndfB+=m_ndf; if(m_ndf==2) ++m_npixelsB;
+ m_cluster = m_linkBackward[0].cluster();
+ if(!addCluster(m_parametersPredBackward,m_parametersUpdatedBackward)) return false;
+ noiseProduction(-1,m_parametersUpdatedBackward);
+ ++m_nclustersBackward; m_xi2totalBackward+=m_xi2Backward; m_ndfBackward+=m_ndf; if(m_ndf==2) ++m_npixelsBackward;
}
- else if(m_xi2B <= m_xi2maxNoAdd) {
+ else if(m_xi2Backward <= m_xi2maxNoAdd) {
- m_clusterNoAdd = m_linkB[0].cluster();
- noiseProduction(-1,m_parametersPB);
+ m_clusterNoAdd = m_linkBackward[0].cluster();
+ noiseProduction(-1,m_parametersPredBackward);
}
}
else {
- noiseProduction(-1,m_parametersPB);
+ noiseProduction(-1,m_parametersPredBackward);
}
if(m_detstatus >=0 && !m_cluster){
- if(m_inside < 0 && !m_clusterNoAdd) {++m_nholesB; ++m_dholesB;}
+ if(m_inside < 0 && !m_clusterNoAdd) {++m_nholesBackward; ++m_dholesBackward;}
if(m_dist < -2.) ++m_ndist;
}
return true;
@@ -389,63 +425,63 @@ bool InDet::SiTrajectoryElement_xk::BackwardPropagationSmoother
double step;
if(TE.m_cluster) {
- if(!propagate(TE.m_parametersUB,m_parametersPB,step)) return false;
- m_dholesB = 0;
+ if(!propagate(TE.m_parametersUpdatedBackward,m_parametersPredBackward,step)) return false;
+ m_dholesBackward = 0;
}
else {
- if(!propagate(TE.m_parametersPB,m_parametersPB,step)) return false;
- m_dholesB = TE.m_dholesB;
+ if(!propagate(TE.m_parametersPredBackward,m_parametersPredBackward,step)) return false;
+ m_dholesBackward = TE.m_dholesBackward;
}
- m_parametersPB.addNoise(m_noise,Trk::oppositeMomentum);
+ m_parametersPredBackward.addNoise(m_noise,Trk::oppositeMomentum);
// Forward-backward predict parameters
//
- if(!combineStates(m_parametersPB,m_parametersPF,m_parametersSM)) return false;
+ if(!combineStates(m_parametersPredBackward,m_parametersPredForward,m_parametersSM)) return false;
m_cluster ? m_status = 3 : m_status = 2;
double Xi2max = m_xi2max; if( isTwoSpacePointsSeed) Xi2max*=2.;
m_inside = m_detlink->intersect(m_parametersSM,m_dist);
- m_nlinksB = 0 ;
- m_clusterOld = m_cluster ;
- m_cluster = 0 ;
- m_clusterNoAdd = 0 ;
- m_xi2B = 10000. ;
- m_nholesB = TE.m_nholesB ;
- m_ndist = TE.m_ndist ;
- m_nclustersB = TE.m_nclustersB ;
- m_npixelsB = TE.m_npixelsB ;
- m_ndfB = TE.m_ndfB ;
- m_xi2totalB = TE.m_xi2totalB ;
+ m_nlinksBackward = 0 ;
+ m_clusterOld = m_cluster ;
+ m_cluster = 0 ;
+ m_clusterNoAdd = 0 ;
+ m_xi2Backward = 10000. ;
+ m_nholesBackward = TE.m_nholesBackward ;
+ m_ndist = TE.m_ndist ;
+ m_nclustersBackward = TE.m_nclustersBackward ;
+ m_npixelsBackward = TE.m_npixelsBackward ;
+ m_ndfBackward = TE.m_ndfBackward ;
+ m_xi2totalBackward = TE.m_xi2totalBackward ;
//m_step += TE.m_step ;
if(m_inside> 0 ) return true;
// For not first cluster on trajectory
//
- if(!m_clusterOld || m_ndfF != m_ndf || m_ndfF+m_ndfB >6) {
+ if(!m_clusterOld || m_ndfForward != m_ndf || m_ndfForward+m_ndfBackward >6) {
- if((m_nlinksB = searchClusters(m_parametersSM,m_linkB))) {
+ if((m_nlinksBackward = searchClusters(m_parametersSM,m_linkBackward))) {
- m_xi2B = m_linkB[0].xi2();
+ m_xi2Backward = m_linkBackward[0].xi2();
- if (m_xi2B <= Xi2max) {
+ if (m_xi2Backward <= Xi2max) {
- m_cluster = m_linkB[0].cluster();
+ m_cluster = m_linkBackward[0].cluster();
- if(!addCluster(m_parametersPB,m_parametersUB)) return false;
- ++m_nclustersB; m_xi2totalB+=m_xi2B; m_ndfB+=m_ndf; if(m_ndf==2) ++m_npixelsB;
+ if(!addCluster(m_parametersPredBackward,m_parametersUpdatedBackward)) return false;
+ ++m_nclustersBackward; m_xi2totalBackward+=m_xi2Backward; m_ndfBackward+=m_ndf; if(m_ndf==2) ++m_npixelsBackward;
}
- else if(m_xi2B <= m_xi2maxNoAdd) {
+ else if(m_xi2Backward <= m_xi2maxNoAdd) {
- m_clusterNoAdd = m_linkB[0].cluster();
+ m_clusterNoAdd = m_linkBackward[0].cluster();
}
}
if(m_detstatus >=0 && !m_cluster) {
- if(m_inside < 0 && !m_clusterNoAdd) {++m_nholesB; ++m_dholesB;}
+ if(m_inside < 0 && !m_clusterNoAdd) {++m_nholesBackward; ++m_dholesBackward;}
if(m_dist < -2.) ++m_ndist;
}
return true;
@@ -454,14 +490,14 @@ bool InDet::SiTrajectoryElement_xk::BackwardPropagationSmoother
// For first cluster of short trajectory
//
m_cluster = m_clusterOld;
- if(!addCluster(m_parametersPB,m_parametersUB,m_xi2B)) return false;
+ if(!addCluster(m_parametersPredBackward,m_parametersUpdatedBackward,m_xi2Backward)) return false;
- if (m_xi2B <= Xi2max ) {++m_nclustersB; m_xi2totalB+=m_xi2B; m_ndfB+=m_ndf; if(m_ndf==2) ++m_npixelsB;}
- else if(m_xi2B <= m_xi2maxNoAdd) {m_cluster = 0; m_clusterNoAdd = m_clusterOld; }
+ if (m_xi2Backward <= Xi2max ) {++m_nclustersBackward; m_xi2totalBackward+=m_xi2Backward; m_ndfBackward+=m_ndf; if(m_ndf==2) ++m_npixelsBackward;}
+ else if(m_xi2Backward <= m_xi2maxNoAdd) {m_cluster = 0; m_clusterNoAdd = m_clusterOld; }
else {m_cluster = 0; }
if(!m_cluster) {
- if(m_inside < 0 && !m_clusterNoAdd) {++m_nholesB; ++m_dholesB;}
+ if(m_inside < 0 && !m_clusterNoAdd) {++m_nholesBackward; ++m_dholesBackward;}
if(m_dist < -2.) ++m_ndist;
}
@@ -474,25 +510,25 @@ bool InDet::SiTrajectoryElement_xk::BackwardPropagationSmoother
bool InDet::SiTrajectoryElement_xk::addNextClusterB()
{
- if(m_nlinksB <= 0) return false;
+ if(m_nlinksBackward <= 0) return false;
- if(m_cluster) m_xi2totalB-=m_xi2B;
+ if(m_cluster) m_xi2totalBackward-=m_xi2Backward;
- if(m_nlinksB > 1 && m_linkB[1].xi2() <= m_xi2max) {
+ if(m_nlinksBackward > 1 && m_linkBackward[1].xi2() <= m_xi2max) {
int n = 0;
- for(; n!=m_nlinksB-1; ++n) m_linkB[n]=m_linkB[n+1];
- m_nlinksB=n;
+ for(; n!=m_nlinksBackward-1; ++n) m_linkBackward[n]=m_linkBackward[n+1];
+ m_nlinksBackward=n;
- m_cluster = m_linkB[0].cluster();
- m_xi2B = m_linkB[0].xi2() ;
- m_xi2totalB+= m_xi2B ;
- if(!addCluster(m_parametersPB,m_parametersUB)) return false;
+ m_cluster = m_linkBackward[0].cluster();
+ m_xi2Backward = m_linkBackward[0].xi2() ;
+ m_xi2totalBackward+= m_xi2Backward ;
+ if(!addCluster(m_parametersPredBackward,m_parametersUpdatedBackward)) return false;
}
else {
- m_nlinksB = 0;
- m_cluster = 0; --m_nclustersB; m_ndfB-=m_ndf; if(m_ndf==2) --m_npixelsB;
- if(m_inside < 0 ) {++m_nholesB; ++m_dholesB;} m_xi2B = 0.;
+ 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;
}
return true;
@@ -504,25 +540,25 @@ bool InDet::SiTrajectoryElement_xk::addNextClusterB()
bool InDet::SiTrajectoryElement_xk::addNextClusterF()
{
- if(m_nlinksF <= 0) return false;
+ if(m_nlinksForward <= 0) return false;
- if(m_cluster) m_xi2totalF-=m_xi2F;
+ if(m_cluster) m_xi2totalForward-=m_xi2Forward;
- if(m_nlinksF > 1 && m_linkF[1].xi2() <= m_xi2max) {
+ if(m_nlinksForward > 1 && m_linkForward[1].xi2() <= m_xi2max) {
int n = 0;
- for(; n!=m_nlinksF-1; ++n) m_linkF[n]=m_linkF[n+1];
- m_nlinksF=n;
+ for(; n!=m_nlinksForward-1; ++n) m_linkForward[n]=m_linkForward[n+1];
+ m_nlinksForward=n;
- m_cluster = m_linkF[0].cluster();
- m_xi2F = m_linkF[0].xi2() ;
- m_xi2totalF+= m_xi2F ;
- if(!addCluster(m_parametersPF,m_parametersUF)) return false;
+ m_cluster = m_linkForward[0].cluster();
+ m_xi2Forward = m_linkForward[0].xi2() ;
+ m_xi2totalForward+= m_xi2Forward ;
+ if(!addCluster(m_parametersPredForward,m_parametersUpdatedForward)) return false;
}
else {
- m_nlinksF = 0;
- m_cluster = 0; --m_nclustersF; m_ndfF-=m_ndf;
- if(m_inside < 0) {++m_nholesF; ++m_dholesF;} m_xi2F = 0.;
+ 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;
}
return true;
@@ -538,24 +574,24 @@ bool InDet::SiTrajectoryElement_xk::addNextClusterB
m_clusterOld = m_cluster ;
m_cluster = Cl ;
- m_nclustersB = TE.m_nclustersB;
- m_npixelsB = TE.m_npixelsB ;
- m_ndfB = TE.m_ndfB ;
- m_nholesB = TE.m_nholesB ;
+ m_nclustersBackward = TE.m_nclustersBackward;
+ m_npixelsBackward = TE.m_npixelsBackward ;
+ m_ndfBackward = TE.m_ndfBackward ;
+ m_nholesBackward = TE.m_nholesBackward ;
m_ndist = TE.m_ndist ;
- m_xi2totalB = TE.m_xi2totalB ;
+ m_xi2totalBackward = TE.m_xi2totalBackward ;
- TE.m_cluster ? m_dholesB = 0 : m_dholesB = TE.m_dholesB;
+ TE.m_cluster ? m_dholesBackward = 0 : m_dholesBackward = TE.m_dholesBackward;
if(Cl) {
- if(!addCluster(m_parametersPB,m_parametersUB,m_xi2B)) return false;
- m_inside = -1; noiseProduction(-1,m_parametersUB);
- ++m_nclustersB; m_xi2totalB+=m_xi2B; m_ndfB+=m_ndf; if(m_ndf==2) ++m_npixelsB;
+ if(!addCluster(m_parametersPredBackward,m_parametersUpdatedBackward,m_xi2Backward)) return false;
+ m_inside = -1; noiseProduction(-1,m_parametersUpdatedBackward);
+ ++m_nclustersBackward; m_xi2totalBackward+=m_xi2Backward; m_ndfBackward+=m_ndf; if(m_ndf==2) ++m_npixelsBackward;
}
else {
- if(m_inside < 0) {++m_nholesB; ++m_dholesB;} m_xi2B = 0.;
+ if(m_inside < 0) {++m_nholesBackward; ++m_dholesBackward;} m_xi2Backward = 0.;
if(m_dist < -2.) ++m_ndist;
}
return true;
@@ -568,26 +604,48 @@ bool InDet::SiTrajectoryElement_xk::addNextClusterB
bool InDet::SiTrajectoryElement_xk::addNextClusterF
(InDet::SiTrajectoryElement_xk& TE,const InDet::SiCluster* Cl)
{
-
+ /// write the current cluster (if any) into the old cluster slot
m_clusterOld = m_cluster ;
- m_cluster = Cl ;
- m_nclustersF = TE.m_nclustersF;
- m_ndfF = TE.m_ndfF ;
- m_nholesF = TE.m_nholesF ;
+ /// and set the current cluster to the one assigned
+ m_cluster = Cl ; /// can be null!
+ /// copy running counts from the previous element
+ m_nclustersForward = TE.m_nclustersForward;
+ m_ndfForward = TE.m_ndfForward ;
+ m_nholesForward = TE.m_nholesForward ;
m_ndist = TE.m_ndist ;
- m_xi2totalF = TE.m_xi2totalF ;
-
- TE.m_cluster ? m_dholesF = 0 : m_dholesF = TE.m_dholesF;
-
- if(Cl) {
-
- if(!addCluster(m_parametersPF,m_parametersUF,m_xi2F)) return false;
- m_inside = -1; noiseProduction(1,m_parametersUF);
- ++m_nclustersF; m_xi2totalF+=m_xi2F; m_ndfF+=m_ndf;
+ m_xi2totalForward = TE.m_xi2totalForward ;
+ /// if the previous element has a hit, reset incremental
+ /// hole counter
+ if (TE.m_cluster){
+ m_dholesForward = 0;
+ }
+ /// otherwise start from the running count
+ else{
+ m_dholesForward = TE.m_dholesForward;
}
+ /// if we have a non-null cluster to add
+ if(Cl) {
+ /// add it to the track, and update our parameters using it
+ if(!addCluster(m_parametersPredForward,m_parametersUpdatedForward,m_xi2Forward)) return false;
+ /// update m_inside to signify we are on the module
+ m_inside = -1;
+ /// update noise to current parameters
+ noiseProduction(1,m_parametersUpdatedForward);
+ /// update hit counts, chi2, ndf
+ ++m_nclustersForward;
+ m_xi2totalForward+=m_xi2Forward;
+ m_ndfForward+=m_ndf;
+ }
+ /// no hit (nullptr cluster)
else {
-
- if(m_inside < 0) {++m_nholesF; ++m_dholesF;} m_xi2F = 0.;
+ /// if we are within boundaries, add to hole counts
+ if(m_inside < 0) {
+ ++m_nholesForward;
+ ++m_dholesForward;
+ }
+ /// no chi2 contribution here
+ m_xi2Forward = 0.;
+ /// if bond gap, increment ndist
if(m_dist < -2.) ++m_ndist;
}
return true;
@@ -609,8 +667,8 @@ InDet::SiTrajectoryElement_xk::trackStateOnSurface (bool change,bool cov,bool mu
const Trk::FitQualityOnSurface* fq = 0;
const Trk::MeasurementBase * ro = 0;
- if (m_status == 1) fq = new Trk::FitQualityOnSurface(m_xi2F,m_ndf);
- else fq = new Trk::FitQualityOnSurface(m_xi2B,m_ndf);
+ if (m_status == 1) fq = new Trk::FitQualityOnSurface(m_xi2Forward,m_ndf);
+ else fq = new Trk::FitQualityOnSurface(m_xi2Backward,m_ndf);
std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> pat(0);
@@ -633,19 +691,19 @@ InDet::SiTrajectoryElement_xk::trackStateOnSurface (bool change,bool cov,bool mu
m_tsos[0] = sos; m_utsos[0] = true; m_ntsos = 1;
- if(multi && m_cluster && m_ndf==2 && m_nlinksB > 1) {
+ if(multi && m_cluster && m_ndf==2 && m_nlinksBackward > 1) {
- for(int i=1; i!= m_nlinksB; ++i) {
+ for(int i=1; i!= m_nlinksBackward; ++i) {
- if(m_linkB[i].xi2() > m_xi2multi) break;
+ if(m_linkBackward[i].xi2() > m_xi2multi) break;
const Trk::TrackParameters * tpn = 0;
if(!change) tpn = trackParameters (cov,Q);
else tpn = trackParametersWithNewDirection(cov,Q);
if(!tpn) break;
- const Trk::FitQualityOnSurface * fqn = new Trk::FitQualityOnSurface(m_linkB[i].xi2(),m_ndf);
- const Trk::MeasurementBase * ron = m_riotool->correct(*m_linkB[i].cluster(),*tp);
+ const Trk::FitQualityOnSurface * fqn = new Trk::FitQualityOnSurface(m_linkBackward[i].xi2(),m_ndf);
+ const Trk::MeasurementBase * ron = m_riotool->correct(*m_linkBackward[i].cluster(),*tp);
const Trk::MaterialEffectsOnTrack* men = new Trk::MaterialEffectsOnTrack(*me);
m_tsos [m_ntsos] = new Trk::TrackStateOnSurface(ron,tpn,fqn,men,pat);
m_utsos[m_ntsos] = false;
@@ -688,8 +746,8 @@ InDet::SiTrajectoryElement_xk::trackSimpleStateOnSurface
Amg::MatrixX cv = cl->localCovariance();
const Trk::FitQualityOnSurface* fq = 0;
- if (m_status == 1) fq = new Trk::FitQualityOnSurface(m_xi2F,m_ndf);
- else fq = new Trk::FitQualityOnSurface(m_xi2B,m_ndf);
+ if (m_status == 1) fq = new Trk::FitQualityOnSurface(m_xi2Forward,m_ndf);
+ else fq = new Trk::FitQualityOnSurface(m_xi2Backward,m_ndf);
if(m_ndf == 1) {
@@ -711,7 +769,7 @@ InDet::SiTrajectoryElement_xk::trackSimpleStateOnSurface
Trk::TrackStateOnSurface*
InDet::SiTrajectoryElement_xk::trackPerigeeStateOnSurface ()
{
- if(m_parametersUB.associatedSurface()!=m_surface) return 0;
+ if(m_parametersUpdatedBackward.associatedSurface()!=m_surface) return 0;
double step ;
Trk::PatternTrackParameters Tp;
@@ -719,7 +777,7 @@ InDet::SiTrajectoryElement_xk::trackPerigeeStateOnSurface ()
Trk::PerigeeSurface per;
bool Q = m_proptool->propagate
- (m_parametersUB,per,Tp,Trk::anyDirection,m_tools->fieldTool(),step,Trk::pion);
+ (m_parametersUpdatedBackward,per,Tp,Trk::anyDirection,m_tools->fieldTool(),step,Trk::pion);
if(Q) {
std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
@@ -740,26 +798,32 @@ const Trk::TrackParameters*
InDet::SiTrajectoryElement_xk::trackParameters(bool cov,int Q)
{
if (m_status == 1) {
- if(m_cluster) return m_parametersUF.convert(cov);
- else return m_parametersPF.convert(cov);
+ if(m_cluster) return m_parametersUpdatedForward.convert(cov);
+ else return m_parametersPredForward.convert(cov);
}
else if(m_status == 2) {
- if(m_cluster) return m_parametersUB.convert(cov);
- else return m_parametersPB.convert(cov);
+ if(m_cluster) return m_parametersUpdatedBackward.convert(cov);
+ else return m_parametersPredBackward.convert(cov);
}
else if(m_status == 3) {
if(Q==0) {
if(m_cluster) {
- if(addCluster(m_parametersSM,m_parametersSM)) return m_parametersSM.convert(cov);
- else if(m_parametersUB.cov()[14] < m_parametersPF.cov()[14]) return m_parametersUB.convert(cov);
- else return m_parametersPF.convert(cov);
+ if(addCluster(m_parametersSM,m_parametersSM)){
+ return m_parametersSM.convert(cov);
+ }
+ else if(m_parametersUpdatedBackward.cov()[14] < m_parametersPredForward.cov()[14]){
+ return m_parametersUpdatedBackward.convert(cov);
+ }
+ else{
+ return m_parametersPredForward.convert(cov);
+ }
}
- else return m_parametersSM.convert(cov);
+ else return m_parametersSM.convert(cov);
}
- if(Q==1) {if(m_cluster) return m_parametersUB.convert(cov);}
- if(Q==2) {if(m_cluster) return m_parametersUF.convert(cov);}
+ if(Q==1) {if(m_cluster) return m_parametersUpdatedBackward.convert(cov);}
+ if(Q==2) {if(m_cluster) return m_parametersUpdatedForward.convert(cov);}
}
return 0;
}
@@ -777,26 +841,49 @@ void InDet::SiTrajectoryElement_xk::noiseProduction
int Model = m_noisemodel;
if(Model < 1 || Model > 2) return;
- double q = fabs(Tp.par()[4]);
- double s = fabs(m_A[0]*m_Tr[6]+m_A[1]*m_Tr[7]+m_A[2]*m_Tr[8]);
- s < .05 ? s = 20. : s = 1./s;
+ double q = fabs(Tp.par()[4]); /// qoverp
+ /// projection of direction normal to surface
+ double s = fabs(m_localDir[0]*m_localTransform[6]+
+ m_localDir[1]*m_localTransform[7]+
+ m_localDir[2]*m_localTransform[8]);
+ if (s < .05) s = 20.;
+ else s = 1./s;
+
+ /// number of radiation lengths when crossing this surface
m_radlengthN = s*m_radlength;
+ /// 134 is a very magic number related to multiple scattering,
+ /// used frequently throught the track finder chain.
+ /// Here we use it to estimate the scattering impact on the errors
double covariancePola = (134.*m_radlengthN)*(q*q);
- double d = (1.-m_A[2])*(1.+m_A[2]); if(d < 1.e-5) d = 1.e-5;
+
+ /// 1 - Az² --> squared transverse direction component
+ double d = (1.-m_localDir[2])*(1.+m_localDir[2]);
+ /// if too small, set to a reasonable minimum
+ if(d < 1.e-5) d = 1.e-5;
+ /// azimuthal component: scale with 1/d
double covarianceAzim = covariancePola/d;
- double covarianceIMom,correctionIMom;
+ double covarianceIMom;
+ double correctionIMom;
+
+ /// muon model
if(Model==1) {
+ /// estimate e-loss
double dp = m_energylose*q*s;
+ /// impact on covariance
covarianceIMom = (.2*dp*dp)*(q*q);
+ /// impact on momentum
correctionIMom = 1.-dp;
}
+ /// electron model - much more generous
else {
correctionIMom = .70;
covarianceIMom = (correctionIMom-1.)*(correctionIMom-1.)*(q*q);
}
+ /// for forward direction, invert correction
if(Dir>0) correctionIMom = 1./correctionIMom;
+ /// store in noise object
m_noise.set(covarianceAzim,covariancePola,covarianceIMom,correctionIMom);
}
@@ -811,26 +898,26 @@ const Trk::TrackParameters*
InDet::SiTrajectoryElement_xk::trackParametersWithNewDirection(bool cov,int Q)
{
if (m_status == 1) {
- if(m_cluster) return trackParameters(m_parametersUF,cov);
- else return trackParameters(m_parametersPF,cov);
+ if(m_cluster) return trackParameters(m_parametersUpdatedForward,cov);
+ else return trackParameters(m_parametersPredForward,cov);
}
else if(m_status == 2) {
- if(m_cluster) return trackParameters(m_parametersUB,cov);
- else return trackParameters(m_parametersPB,cov);
+ if(m_cluster) return trackParameters(m_parametersUpdatedBackward,cov);
+ else return trackParameters(m_parametersPredBackward,cov);
}
else if(m_status == 3) {
if(Q==0) {
if(m_cluster) {
if(addCluster(m_parametersSM,m_parametersSM)) return trackParameters(m_parametersSM,cov);
- else if(m_parametersUB.cov()[14] < m_parametersPF.cov()[14]) return trackParameters(m_parametersUB,cov);
- else return trackParameters(m_parametersPF,cov);
+ else if(m_parametersUpdatedBackward.cov()[14] < m_parametersPredForward.cov()[14]) return trackParameters(m_parametersUpdatedBackward,cov);
+ else return trackParameters(m_parametersPredForward,cov);
}
else return trackParameters(m_parametersSM,cov);
}
- if(Q==1) {if(m_cluster) return trackParameters(m_parametersUB,cov);}
- if(Q==2) {if(m_cluster) return trackParameters(m_parametersUF,cov);}
+ if(Q==1) {if(m_cluster) return trackParameters(m_parametersUpdatedBackward,cov);}
+ if(Q==2) {if(m_cluster) return trackParameters(m_parametersUpdatedForward,cov);}
}
return 0;
}
@@ -858,12 +945,12 @@ double InDet::SiTrajectoryElement_xk::step
Trk::PatternTrackParameters Ta,Tb;
if (TE.m_status == 1) {
- if(TE.m_cluster) Ta = TE.m_parametersUF;
- else Ta = TE.m_parametersPF;
+ if(TE.m_cluster) Ta = TE.m_parametersUpdatedForward;
+ else Ta = TE.m_parametersPredForward;
}
else if(TE.m_status == 2) {
- if(TE.m_cluster) Ta = TE.m_parametersUB;
- else Ta = TE.m_parametersPB;
+ if(TE.m_cluster) Ta = TE.m_parametersUpdatedBackward;
+ else Ta = TE.m_parametersPredBackward;
}
else if(TE.m_status == 3) {
Ta = TE.m_parametersSM;
@@ -881,12 +968,12 @@ double InDet::SiTrajectoryElement_xk::step
Amg::Vector3D InDet::SiTrajectoryElement_xk::globalPosition()
{
if (m_status == 1) {
- if(m_cluster) return m_parametersUF.position();
- else return m_parametersPF.position();
+ if(m_cluster) return m_parametersUpdatedForward.position();
+ else return m_parametersPredForward.position();
}
else if(m_status == 2) {
- if(m_cluster) return m_parametersUB.position();
- else return m_parametersPB.position();
+ if(m_cluster) return m_parametersUpdatedBackward.position();
+ else return m_parametersPredBackward.position();
}
else if(m_status == 3) {
@@ -917,12 +1004,12 @@ double InDet::SiTrajectoryElement_xk::stepToPerigee()
Amg::Vector3D P;
if(m_cluster) {
- M = m_parametersUB.momentum();
- P = m_parametersUB.position();
+ M = m_parametersUpdatedBackward.momentum();
+ P = m_parametersUpdatedBackward.position();
}
else {
- M = m_parametersPB.momentum();
- P = m_parametersPB.position();
+ M = m_parametersPredBackward.momentum();
+ P = m_parametersPredBackward.position();
}
return -(P[0]*M[0]+P[1]*M[1]);
@@ -934,7 +1021,10 @@ double InDet::SiTrajectoryElement_xk::stepToPerigee()
void InDet::SiTrajectoryElement_xk::eraseClusterForwardPropagation()
{
- m_cluster=0; --m_nclustersF; m_xi2totalF-=m_xi2F; m_ndfF-=m_ndf;
+ m_cluster=0;
+ --m_nclustersForward;
+ m_xi2totalForward-=m_xi2Forward;
+ m_ndfForward-=m_ndf;
}
///////////////////////////////////////////////////////////////////
@@ -956,7 +1046,7 @@ double InDet::SiTrajectoryElement_xk::quality(int& holes) const
}
double w,X,Xc = m_xi2max+2.;
- m_status == 1 ? X = m_xi2F : X = m_xi2B;
+ m_status == 1 ? X = m_xi2Forward : X = m_xi2Backward;
m_ndf == 2 ? w = 1.2*(Xc-X*.5) : w = Xc-X ; if(w < -1.) w = -1.;
holes = 0;
@@ -965,211 +1055,351 @@ double InDet::SiTrajectoryElement_xk::quality(int& holes) const
/////////////////////////////////////////////////////////////////////////////////
// Main function for pattern track parameters and covariance matrix propagation
-// to PlaneSurface Ta->pSu = Tb with step to surface calculation
+// to PlaneSurface.
+/// Start from 'startingParameters', propagate to current surface.
+/// Write into outputParameters
+/// and wrote step to surface into StepLength
/////////////////////////////////////////////////////////////////////////////////
-
bool
-InDet::SiTrajectoryElement_xk::propagate(Trk::PatternTrackParameters & Ta,
- Trk::PatternTrackParameters & Tb,
- double & S ) {
- bool useJac = true;
- if(!Ta.iscovariance()) useJac = false;
- double P[64];
- transformPlaneToGlobal(useJac,Ta,P);
- if( m_fieldMode) {if(!rungeKuttaToPlane (useJac,P)) return false;}
- else {if(!straightLineStepToPlane(useJac,P)) return false;}
- S = P[45]; return transformGlobalToPlane(useJac,P,Ta,Tb);
+InDet::SiTrajectoryElement_xk::propagate(Trk::PatternTrackParameters & startingParameters,
+ Trk::PatternTrackParameters & outputParameters,
+ double & StepLength ) {
+ bool useJac = (startingParameters.iscovariance());
+ double globalParameters[64]; /// helper field to store global frame parameters
+ /// Convention:
+ /// 0-2: Position
+ /// 3-5: Direction
+ /// 6: qoverp
+ /// 7 - 41: jacobian
+ /// 42-44: extension of jacobian for dA/ds
+ /// 45: step length
+
+ /// transform starting parameters to global, write into globalParameters
+ transformPlaneToGlobal(useJac,startingParameters,globalParameters);
+ /// if running with field, use Runge Kutta. Will update globalParameters to the result in global coordinates
+ if( m_fieldMode) {
+ if(!rungeKuttaToPlane (useJac,globalParameters)) return false;
+ }
+ /// otherwise, we can use the straight line
+ else{
+ if(!straightLineStepToPlane(useJac,globalParameters)) return false;
+ }
+ /// Update step length
+ StepLength = globalParameters[45];
+ /// and finally transform from global back to local, and write into the output parameters.
+ return transformGlobalToPlane(useJac,globalParameters,startingParameters,outputParameters);
}
/////////////////////////////////////////////////////////////////////////////////
-// Main function for pattern track parameters and covariance matrix propagation
-// to PlaneSurface Ta->pSu = Tb with step to surface calculation
+// Main function for pattern track parameters propagation without covariance.
+/// Start from 'startingParameters', propagate to current surface.
+/// Write into outputParameters
+/// and wrote step to surface into StepLength
/////////////////////////////////////////////////////////////////////////////////
bool
-InDet::SiTrajectoryElement_xk::propagateParameters(Trk::PatternTrackParameters & Ta,
- Trk::PatternTrackParameters & Tb, double & S ) {
+InDet::SiTrajectoryElement_xk::propagateParameters(Trk::PatternTrackParameters & startingParameters,
+ Trk::PatternTrackParameters & outputParameters,
+ double & StepLength ) {
bool useJac = false;
- double P[64];
- transformPlaneToGlobal(useJac,Ta,P);
+ double globalParameters[64]; /// helper field to store global frame parameters
+ /// Convention:
+ /// 0-2: Position
+ /// 3-5: Direction
+ /// 6: qoverp
+ /// 7 - 41: jacobian
+ /// 42-44: extension of jacobian for dA/ds
+ /// 45: step length
+ transformPlaneToGlobal(useJac,startingParameters,globalParameters);
+ /// if running with field, use Runge Kutta.
+ /// Will update globalParameters to the result in global coordinates
if( m_fieldMode) {
- if(!rungeKuttaToPlane (useJac,P)) return false;
- } else {
- if(!straightLineStepToPlane(useJac,P)) return false;
+ if(!rungeKuttaToPlane (useJac,globalParameters)) return false;
+ }
+ /// otherwise, we can use the straight line
+ else {
+ if(!straightLineStepToPlane(useJac,globalParameters)) return false;
}
- S = P[45];
- return transformGlobalToPlane(useJac,P,Ta,Tb);
+ /// Update step length
+ StepLength = globalParameters[45];
+ /// and finally transform from global back to local, and write into the output parameters.
+ return transformGlobalToPlane(useJac,globalParameters,startingParameters,outputParameters);
}
/////////////////////////////////////////////////////////////////////////////////
-// Tramsform from plane to global
-/////////////////////////////////////////////////////////////////////////////////
+/// Tramsform from plane to global
+/// Will take the surface and parameters from localParameters
+/// and populate globalPars with the result.
+/// globalPars needs to be at least 46 elements long.
+/// Convention:
+/// /dL0 /dL1 /dPhi /dThe /dCM
+/// X ->globalPars[0] dX / globalPars[ 7] globalPars[14] globalPars[21] globalPars[28] globalPars[35]
+/// Y ->globalPars[1] dY / globalPars[ 8] globalPars[15] globalPars[22] globalPars[29] globalPars[36]
+/// Z ->globalPars[2] dZ / globalPars[ 9] globalPars[16] globalPars[23] globalPars[30] globalPars[37]
+/// Ax ->globalPars[3] dAx/ globalPars[10] globalPars[17] globalPars[24] globalPars[31] globalPars[38]
+/// Ay ->globalPars[4] dAy/ globalPars[11] globalPars[18] globalPars[25] globalPars[32] globalPars[39]
+/// Az ->globalPars[5] dAz/ globalPars[12] globalPars[19] globalPars[26] globalPars[33] globalPars[40]
+/// CM ->globalPars[6] dCM/ globalPars[13] globalPars[20] globalPars[27] globalPars[34] globalPars[41]
+// /////////////////////////////////////////////////////////////////////////////////
void InDet::SiTrajectoryElement_xk::transformPlaneToGlobal(bool useJac,
- Trk::PatternTrackParameters& Tp,double* P) {
- double Sf,Cf,Ce,Se;
- sincos(Tp.par()[2],&Sf,&Cf);
- sincos(Tp.par()[3],&Se,&Ce);
- const auto pSurface=Tp.associatedSurface();
- if (not pSurface){
+ Trk::PatternTrackParameters& localParameters,
+ double* globalPars) {
+ /// obtain trigonometric functions required for transform
+ double sinPhi,cosPhi,cosTheta,sintheta;
+ sincos(localParameters.par()[2],&sinPhi,&cosPhi);
+ sincos(localParameters.par()[3],&sintheta,&cosTheta);
+ /// get the surface corresponding to the local parameters
+ const Trk::Surface* pSurface=localParameters.associatedSurface();
+ if (!pSurface){
throw(std::runtime_error("TrackParameters associated surface is null pointer in InDet::SiTrajectoryElement_xk::transformPlaneToGlobal"));
}
+ /// get the associated transform from the surface
const Amg::Transform3D& T = pSurface->transform();
+ /// Get the local x and y axis in the global frame
double Ax[3] = {T(0,0),T(1,0),T(2,0)};
double Ay[3] = {T(0,1),T(1,1),T(2,1)};
- P[ 0] = Tp.par()[0]*Ax[0]+Tp.par()[1]*Ay[0]+T(0,3); // X
- P[ 1] = Tp.par()[0]*Ax[1]+Tp.par()[1]*Ay[1]+T(1,3); // Y
- P[ 2] = Tp.par()[0]*Ax[2]+Tp.par()[1]*Ay[2]+T(2,3); // Z
- P[ 3] = Cf*Se; // Ax
- P[ 4] = Sf*Se; // Ay
- P[ 5] = Ce; // Az
- P[ 6] = Tp.par()[4]; // CM
- if(fabs(P[6])<1.e-20) {P[6] < 0. ? P[6]=-1.e-20 : P[6]= 1.e-20;}
+ /// position
+ globalPars[ 0] = localParameters.par()[0]*Ax[0]+localParameters.par()[1]*Ay[0]+T(0,3); // X
+ globalPars[ 1] = localParameters.par()[0]*Ax[1]+localParameters.par()[1]*Ay[1]+T(1,3); // Y
+ globalPars[ 2] = localParameters.par()[0]*Ax[2]+localParameters.par()[1]*Ay[2]+T(2,3); // Z
+ /// direction vectors
+ globalPars[ 3] = cosPhi*sintheta; // Ax
+ globalPars[ 4] = sinPhi*sintheta; // Ay
+ globalPars[ 5] = cosTheta;
+ /// qoeverp // Az
+ globalPars[ 6] = localParameters.par()[4]; // CM
+ /// for very high momenta, truncate to avoid zero
+ if(std::abs(globalPars[6])<1.e-20) {
+ if (globalPars[6] < 0){
+ globalPars[6]=-1.e-20;
+ }
+ else globalPars[6]= 1.e-20;
+ }
+ /// Update jacobian if requested
if(useJac) {
// /dL1 | /dL2 | /dPhi | /dThe | /dCM |
- P[ 7] = Ax[0]; P[14] = Ay[0]; P[21] = 0.; P[28] = 0.; P[35] = 0.; // dX /
- P[ 8] = Ax[1]; P[15] = Ay[1]; P[22] = 0.; P[29] = 0.; P[36] = 0.; // dY /
- P[ 9] = Ax[2]; P[16] = Ay[2]; P[23] = 0.; P[30] = 0.; P[37] = 0.; // dZ /
- P[10] = 0.; P[17] = 0.; P[24] =-P[4]; P[31] = Cf*Ce; P[38] = 0.; // dAx/
- P[11] = 0.; P[18] = 0.; P[25] = P[3]; P[32] = Sf*Ce; P[39] = 0.; // dAy/
- P[12] = 0.; P[19] = 0.; P[26] = 0.; P[33] = -Se; P[40] = 0.; // dAz/
- P[42] = 0.; P[43] = 0.; P[44] = 0.; // d(Ax,Ay,Az)/ds
- }
- P[45] = 0.;
+ globalPars[ 7] = Ax[0]; globalPars[14] = Ay[0]; globalPars[21] = 0.; globalPars[28] = 0.; globalPars[35] = 0.; // dX /
+ globalPars[ 8] = Ax[1]; globalPars[15] = Ay[1]; globalPars[22] = 0.; globalPars[29] = 0.; globalPars[36] = 0.; // dY /
+ globalPars[ 9] = Ax[2]; globalPars[16] = Ay[2]; globalPars[23] = 0.; globalPars[30] = 0.; globalPars[37] = 0.; // dZ /
+ globalPars[10] = 0.; globalPars[17] = 0.; globalPars[24] =-globalPars[4]; globalPars[31] = cosPhi*cosTheta; globalPars[38] = 0.; // dAx/
+ globalPars[11] = 0.; globalPars[18] = 0.; globalPars[25] = globalPars[3]; globalPars[32] = sinPhi*cosTheta; globalPars[39] = 0.; // dAy/
+ globalPars[12] = 0.; globalPars[19] = 0.; globalPars[26] = 0.; globalPars[33] = -sintheta; globalPars[40] = 0.; // dAz/
+
+ /// d(Ax,Ay,Az)/ds
+ globalPars[42] = 0.;
+ globalPars[43] = 0.;
+ globalPars[44] = 0.;
+ }
+ /// Step length is initialised to zero - no propagation done yet
+ globalPars[45] = 0.;
}
/////////////////////////////////////////////////////////////////////////////////
-// Tramsform from global to plane
-/////////////////////////////////////////////////////////////////////////////////
-
+/// Tramsform from global to plane surface
+/// Will take the global parameters in globalPars, the surface from this trajectory
+/// element and populate outputParameters with the result.
+/// Convention for globalPars:
+/// /dL0 /dL1 /dPhi /dThe /dCM
+/// X ->globalPars[0] dX / globalPars[ 7] globalPars[14] globalPars[21] globalPars[28] globalPars[35]
+/// Y ->globalPars[1] dY / globalPars[ 8] globalPars[15] globalPars[22] globalPars[29] globalPars[36]
+/// Z ->globalPars[2] dZ / globalPars[ 9] globalPars[16] globalPars[23] globalPars[30] globalPars[37]
+/// Ax ->globalPars[3] dAx/ globalPars[10] globalPars[17] globalPars[24] globalPars[31] globalPars[38]
+/// Ay ->globalPars[4] dAy/ globalPars[11] globalPars[18] globalPars[25] globalPars[32] globalPars[39]
+/// Az ->globalPars[5] dAz/ globalPars[12] globalPars[19] globalPars[26] globalPars[33] globalPars[40]
+/// CM ->globalPars[6] dCM/ globalPars[13] globalPars[20] globalPars[27] globalPars[34] globalPars[41]
+// /////////////////////////////////////////////////////////////////////////////////
bool InDet::SiTrajectoryElement_xk::transformGlobalToPlane
-(bool useJac,double* P,Trk::PatternTrackParameters& Ta,Trk::PatternTrackParameters& Tb)
+(bool useJac,double* globalPars,Trk::PatternTrackParameters& startingParameters,Trk::PatternTrackParameters& outputParameters)
{
-
- double Ax[3] = {m_Tr[0],m_Tr[1],m_Tr[2]};
- double Ay[3] = {m_Tr[3],m_Tr[4],m_Tr[5]};
- double Az[3] = {m_Tr[6],m_Tr[7],m_Tr[8]};
- double d [3] = {P[0]-m_Tr[ 9],P[1]-m_Tr[10],P[2]-m_Tr[11]};
-
- double p[5] = {d[0]*Ax[0]+d[1]*Ax[1]+d[2]*Ax[2],
- d[0]*Ay[0]+d[1]*Ay[1]+d[2]*Ay[2],
- atan2(P[4],P[3]) ,
- acos(P[5]) ,
- P[6] };
-
- Tb.setParameters(m_surface,p); m_A[0] = P[3]; m_A[1] = P[4]; m_A[2] = P[5];
+ /// local x,y,z axes in global frame
+ double Ax[3] = {m_localTransform[0],m_localTransform[1],m_localTransform[2]};
+ double Ay[3] = {m_localTransform[3],m_localTransform[4],m_localTransform[5]};
+ double Az[3] = {m_localTransform[6],m_localTransform[7],m_localTransform[8]};
+ /// distance of global position to origin of local frame
+ double d [3] = {globalPars[0]-m_localTransform[ 9],
+ globalPars[1]-m_localTransform[10],
+ globalPars[2]-m_localTransform[11]};
+
+ /// local parameters obtained by transforming from global
+ double p[5] = {
+ d[0]*Ax[0]+d[1]*Ax[1]+d[2]*Ax[2], /// local X
+ d[0]*Ay[0]+d[1]*Ay[1]+d[2]*Ay[2], /// local Y
+ atan2(globalPars[4],globalPars[3]), /// phi (from global direction x,y)
+ acos(globalPars[5]), /// theta (from global cosTheta)
+ globalPars[6] /// qoverp
+ };
+ /// write into output parameters. Assign our surface to them
+ outputParameters.setParameters(m_surface,p);
+ /// update local direction vector
+ m_localDir[0] = globalPars[3];
+ m_localDir[1] = globalPars[4];
+ m_localDir[2] = globalPars[5];
+
+ /// if we don't need the cov, we are done
if(!useJac) return true;
- // Condition trajectory on surface
- //
- double A = Az[0]*P[3]+Az[1]*P[4]+Az[2]*P[5]; if(A!=0.) A=1./A;
- double s0 = Az[0]*P[ 7]+Az[1]*P[ 8]+Az[2]*P[ 9];
- double s1 = Az[0]*P[14]+Az[1]*P[15]+Az[2]*P[16];
- double s2 = Az[0]*P[21]+Az[1]*P[22]+Az[2]*P[23];
- double s3 = Az[0]*P[28]+Az[1]*P[29]+Az[2]*P[30];
- double s4 = Az[0]*P[35]+Az[1]*P[36]+Az[2]*P[37];
- double T0 =(Ax[0]*P[ 3]+Ax[1]*P[ 4]+Ax[2]*P[ 5])*A;
- double T1 =(Ay[0]*P[ 3]+Ay[1]*P[ 4]+Ay[2]*P[ 5])*A;
- double n = 1./P[6];
+ /// Condition trajectory on surface
+ double A = Az[0]*globalPars[3]+Az[1]*globalPars[4]+Az[2]*globalPars[5];
+ if(A!=0.) A=1./A;
+ double s0 = Az[0]*globalPars[ 7]+Az[1]*globalPars[ 8]+Az[2]*globalPars[ 9];
+ double s1 = Az[0]*globalPars[14]+Az[1]*globalPars[15]+Az[2]*globalPars[16];
+ double s2 = Az[0]*globalPars[21]+Az[1]*globalPars[22]+Az[2]*globalPars[23];
+ double s3 = Az[0]*globalPars[28]+Az[1]*globalPars[29]+Az[2]*globalPars[30];
+ double s4 = Az[0]*globalPars[35]+Az[1]*globalPars[36]+Az[2]*globalPars[37];
+ double T0 =(Ax[0]*globalPars[ 3]+Ax[1]*globalPars[ 4]+Ax[2]*globalPars[ 5])*A;
+ double T1 =(Ay[0]*globalPars[ 3]+Ay[1]*globalPars[ 4]+Ay[2]*globalPars[ 5])*A;
+ double n = 1./globalPars[6];
double Jac[21];
// Jacobian production
//
- Jac[ 0] = (Ax[0]*P[ 7]+Ax[1]*P[ 8])+(Ax[2]*P[ 9]-s0*T0); // dL0/dL0
- Jac[ 1] = (Ax[0]*P[14]+Ax[1]*P[15])+(Ax[2]*P[16]-s1*T0); // dL0/dL1
- Jac[ 2] = (Ax[0]*P[21]+Ax[1]*P[22])+(Ax[2]*P[23]-s2*T0); // dL0/dPhi
- Jac[ 3] = (Ax[0]*P[28]+Ax[1]*P[29])+(Ax[2]*P[30]-s3*T0); // dL0/dThe
- Jac[ 4] =((Ax[0]*P[35]+Ax[1]*P[36])+(Ax[2]*P[37]-s4*T0))*n; // dL0/dCM
-
- Jac[ 5] = (Ay[0]*P[ 7]+Ay[1]*P[ 8])+(Ay[2]*P[ 9]-s0*T1); // dL1/dL0
- Jac[ 6] = (Ay[0]*P[14]+Ay[1]*P[15])+(Ay[2]*P[16]-s1*T1); // dL1/dL1
- Jac[ 7] = (Ay[0]*P[21]+Ay[1]*P[22])+(Ay[2]*P[23]-s2*T1); // dL1/dPhi
- Jac[ 8] = (Ay[0]*P[28]+Ay[1]*P[29])+(Ay[2]*P[30]-s3*T1); // dL1/dThe
- Jac[ 9] =((Ay[0]*P[35]+Ay[1]*P[36])+(Ay[2]*P[37]-s4*T1))*n; // dL1/dCM
-
- double P3,P4, C = P[3]*P[3]+P[4]*P[4];
- if(C > 1.e-20) {C= 1./C ; P3 = P[3]*C; P4 =P[4]*C; C =-sqrt(C);}
- else {C=-1.e10; P3 = 1. ; P4 =0. ; }
-
- double T2 =(P3*P[43]-P4*P[42])*A;
- double C44 = C*P[44] *A;
-
- Jac[10] = P3*P[11]-P4*P[10]-s0*T2; // dPhi/dL0
- Jac[11] = P3*P[18]-P4*P[17]-s1*T2; // dPhi/dL1
- Jac[12] = P3*P[25]-P4*P[24]-s2*T2; // dPhi/dPhi
- Jac[13] = P3*P[32]-P4*P[31]-s3*T2; // dPhi/dThe
- Jac[14] =(P3*P[39]-P4*P[38]-s4*T2)*n; // dPhi/dCM
-
- Jac[15] = C*P[12]-s0*C44; // dThe/dL0
- Jac[16] = C*P[19]-s1*C44; // dThe/dL1
- Jac[17] = C*P[26]-s2*C44; // dThe/dPhi
- Jac[18] = C*P[33]-s3*C44; // dThe/dThe
- Jac[19] =(C*P[40]-s4*C44)*n; // dThe/dCM
+ Jac[ 0] = (Ax[0]*globalPars[ 7]+Ax[1]*globalPars[ 8])+(Ax[2]*globalPars[ 9]-s0*T0); // dL0/dL0
+ Jac[ 1] = (Ax[0]*globalPars[14]+Ax[1]*globalPars[15])+(Ax[2]*globalPars[16]-s1*T0); // dL0/dL1
+ Jac[ 2] = (Ax[0]*globalPars[21]+Ax[1]*globalPars[22])+(Ax[2]*globalPars[23]-s2*T0); // dL0/dPhi
+ Jac[ 3] = (Ax[0]*globalPars[28]+Ax[1]*globalPars[29])+(Ax[2]*globalPars[30]-s3*T0); // dL0/dThe
+ Jac[ 4] =((Ax[0]*globalPars[35]+Ax[1]*globalPars[36])+(Ax[2]*globalPars[37]-s4*T0))*n; // dL0/dCM
+
+ Jac[ 5] = (Ay[0]*globalPars[ 7]+Ay[1]*globalPars[ 8])+(Ay[2]*globalPars[ 9]-s0*T1); // dL1/dL0
+ Jac[ 6] = (Ay[0]*globalPars[14]+Ay[1]*globalPars[15])+(Ay[2]*globalPars[16]-s1*T1); // dL1/dL1
+ Jac[ 7] = (Ay[0]*globalPars[21]+Ay[1]*globalPars[22])+(Ay[2]*globalPars[23]-s2*T1); // dL1/dPhi
+ Jac[ 8] = (Ay[0]*globalPars[28]+Ay[1]*globalPars[29])+(Ay[2]*globalPars[30]-s3*T1); // dL1/dThe
+ Jac[ 9] =((Ay[0]*globalPars[35]+Ay[1]*globalPars[36])+(Ay[2]*globalPars[37]-s4*T1))*n; // dL1/dCM
+
+ double P3=0;
+ double P4=0;
+ /// transverse direction component
+ double C = globalPars[3]*globalPars[3]+globalPars[4]*globalPars[4];
+ if(C > 1.e-20) {
+ C= 1./C ;
+ /// unit direction vector in the X,Y plane
+ P3 = globalPars[3]*C;
+ P4 =globalPars[4]*C;
+ C =-sqrt(C);
+ }
+ else{
+ C=-1.e10;
+ P3 = 1.;
+ P4 =0.;
+ }
+
+ double T2 =(P3*globalPars[43]-P4*globalPars[42])*A;
+ double C44 = C*globalPars[44] *A;
+
+ Jac[10] = P3*globalPars[11]-P4*globalPars[10]-s0*T2; // dPhi/dL0
+ Jac[11] = P3*globalPars[18]-P4*globalPars[17]-s1*T2; // dPhi/dL1
+ Jac[12] = P3*globalPars[25]-P4*globalPars[24]-s2*T2; // dPhi/dPhi
+ Jac[13] = P3*globalPars[32]-P4*globalPars[31]-s3*T2; // dPhi/dThe
+ Jac[14] =(P3*globalPars[39]-P4*globalPars[38]-s4*T2)*n; // dPhi/dCM
+
+ Jac[15] = C*globalPars[12]-s0*C44; // dThe/dL0
+ Jac[16] = C*globalPars[19]-s1*C44; // dThe/dL1
+ Jac[17] = C*globalPars[26]-s2*C44; // dThe/dPhi
+ Jac[18] = C*globalPars[33]-s3*C44; // dThe/dThe
+ Jac[19] =(C*globalPars[40]-s4*C44)*n; // dThe/dCM
Jac[20] = 1.; // dCM /dCM
- Tb.newCovarianceMatrix(Ta,Jac);
- const double* t = &Tb.cov()[0];
+ /// covariance matrix production using jacobian - CovNEW = J*CovOLD*Jt
+ outputParameters.newCovarianceMatrix(startingParameters,Jac);
+ /// check for negative diagonals in the cov
+ const double* t = &outputParameters.cov()[0];
if(t[0]<=0. || t[2]<=0. || t[5]<=0. || t[9]<=0. || t[14]<=0.) return false;
return true;
}
/////////////////////////////////////////////////////////////////////////////////
-// Runge Kutta step to plane
+/// Runge Kutta step to plane
+/// Updates the "globalPars" array, which is also used to
+/// pass the input.
/////////////////////////////////////////////////////////////////////////////////
bool InDet::SiTrajectoryElement_xk::rungeKuttaToPlane
-(bool Jac,double* P)
+(bool Jac,double* globalPars)
{
+ /// minimum step length
const double Smin = .1 ;
+ /// step length below which we are allowed to use
+ /// a helical approximation
const double Shel = 5. ;
+ /// tolerance before reducing step length
const double dlt = .001 ;
- if(fabs(P[6]) > .05) return false;
+ /// Minimum momentum check
+ if(fabs(globalPars[6]) > .05) return false;
+
int it = 0;
- double* R = &P[ 0]; // Coordinates
- double* A = &P[ 3]; // Directions
- double* sA = &P[42];
- double Pi = 149.89626*P[6]; // Invert mometum/2.
- double Pa = fabs (P[6]);
-
- double a = A[0]*m_Tr[6]+A[1]*m_Tr[7]+A[2]*m_Tr[8] ; if(a==0.) return false;
- double S = ((m_Tr[12]-R[0]*m_Tr[6])-(R[1]*m_Tr[7]+R[2]*m_Tr[8]))/a;
+ double* R = &globalPars[ 0]; // Coordinates
+
+ double* A = &globalPars[ 3]; // Directions
+
+ double* sA = &globalPars[42]; /// dA/ds
+ double Pi = 149.89626*globalPars[6]; /// This is a conversion from p to half-curvature (1/2R), when multiplied with the orthogonal field component
+ /// Curvature = 2 * Pi x B x sin(angle)
+ double Pa = fabs (globalPars[6]); /// abs(qoverp)
+
+ /// projection of global direction onto local Z-axis
+ double a = A[0]*m_localTransform[6]+A[1]*m_localTransform[7]+A[2]*m_localTransform[8];
+ if(a==0.) return false;
+ /// distance orthogonal to surface in units of direction z-component
+ double S = ((m_localTransform[12]-R[0]*m_localTransform[6])-(R[1]*m_localTransform[7]+R[2]*m_localTransform[8]))/a;
double S0 = fabs(S) ;
+ /// if we are sufficiently close in Z already, add a
+ /// final straight line step corresponding to the remaining difference
+ /// to get the estimate on the surface, update step length, and return the result
if(S0 <= Smin) {
- R[0]+=(A[0]*S); R[1]+=(A[1]*S); R[2]+=(A[2]*S); P[45]+=S; return true;
+ R[0]+=(A[0]*S);
+ R[1]+=(A[1]*S);
+ R[2]+=(A[2]*S);
+ globalPars[45]+=S;
+ return true;
}
+ /// if we are too far away,
+ /// Limit the step size to 0.3 * p
else if( (Pa*S0) > .3) {
- S > 0. ? S = .3/Pa : S=-.3/Pa;
+ if (S >0) S = 0.3 / Pa;
+ else S = -0.3/Pa;
}
bool ste = false;
- double f0[3],f[3];
+ /// holders for the B-field
+ double f0[3];
+ double f[3];
+ /// get the field at our reference point
m_fieldCache.getFieldZR(R,f0);
while(true) {
- bool Helix = false; if(fabs(S) < Shel) Helix = true;
- double S3=(1./3.)*S, S4=.25*S, PS2=Pi*S;
+ bool Helix = false;
+ /// if the step is sifficiently small, we are allowed to use a helical model
+ if(fabs(S) < Shel) Helix = true;
+ double S3=(1./3.)*S; /// S/3
+ double S4=.25*S; /// S/4
+ double PS2=Pi*S; /// 2 S/curvature radius, when multiplied with the appropriate field
- // First point
- //
- double H0[3] = {f0[0]*PS2, f0[1]*PS2, f0[2]*PS2};
- double A0 = A[1]*H0[2]-A[2]*H0[1] ;
+ /// First point
+ ///
+ double H0[3] = {f0[0]*PS2, /// equip field with conversion factor to 2S/curvature
+ f0[1]*PS2,
+ f0[2]*PS2};
+
+ double A0 = A[1]*H0[2]-A[2]*H0[1] ; /// v cross B
double B0 = A[2]*H0[0]-A[0]*H0[2] ;
double C0 = A[0]*H0[1]-A[1]*H0[0] ;
+
double A2 = A0+A[0] ;
double B2 = B0+A[1] ;
double C2 = C0+A[2] ;
+
double A1 = A2+A[0] ;
double B1 = B2+A[1] ;
double C1 = C2+A[2] ;
@@ -1177,20 +1407,31 @@ bool InDet::SiTrajectoryElement_xk::rungeKuttaToPlane
// Second point
//
if(!Helix) {
- double gP[3]={R[0]+A1*S4, R[1]+B1*S4, R[2]+C1*S4};
+ double gP[3]={R[0]+A1*S4, /// updated position for field eval
+ R[1]+B1*S4,
+ R[2]+C1*S4};
- m_fieldCache.getFieldZR(gP,f);
+ m_fieldCache.getFieldZR(gP,f); /// field at updated position
}
- else {f[0]=f0[0]; f[1]=f0[1]; f[2]=f0[2];}
+ /// if assuming helix, homogenous field
+ else {
+ f[0]=f0[0];
+ f[1]=f0[1];
+ f[2]=f0[2];
+ }
- double H1[3] = {f[0]*PS2,f[1]*PS2,f[2]*PS2};
+ double H1[3] = {f[0]*PS2,
+ f[1]*PS2,
+ f[2]*PS2};
double A3 = (A[0]+B2*H1[2])-C2*H1[1] ;
double B3 = (A[1]+C2*H1[0])-A2*H1[2] ;
double C3 = (A[2]+A2*H1[1])-B2*H1[0] ;
+
double A4 = (A[0]+B3*H1[2])-C3*H1[1] ;
double B4 = (A[1]+C3*H1[0])-A3*H1[2] ;
double C4 = (A[2]+A3*H1[1])-B3*H1[0] ;
+
double A5 = 2.*A4-A[0] ;
double B5 = 2.*B4-A[1] ;
double C5 = 2.*C4-A[2] ;
@@ -1198,14 +1439,23 @@ bool InDet::SiTrajectoryElement_xk::rungeKuttaToPlane
// Last point
//
if(!Helix) {
- double gP[3]={R[0]+S*A4, R[1]+S*B4, R[2]+S*C4};
+ double gP[3]={R[0]+S*A4,
+ R[1]+S*B4,
+ R[2]+S*C4};
m_fieldCache.getFieldZR(gP,f);
}
- else {f[0]=f0[0]; f[1]=f0[1]; f[2]=f0[2];}
+ else{
+ f[0]=f0[0];
+ f[1]=f0[1];
+ f[2]=f0[2];
+ }
- double H2[3] = {f[0]*PS2,f[1]*PS2,f[2]*PS2};
+ double H2[3] = {f[0]*PS2,
+ f[1]*PS2,
+ f[2]*PS2};
+
double A6 = B5*H2[2]-C5*H2[1] ;
double B6 = C5*H2[0]-A5*H2[2] ;
double C6 = A5*H2[1]-B5*H2[0] ;
@@ -1214,12 +1464,15 @@ bool InDet::SiTrajectoryElement_xk::rungeKuttaToPlane
//
if(!ste) {
double EST = fabs((A1+A6)-(A3+A4))+fabs((B1+B6)-(B3+B4))+fabs((C1+C6)-(C3+C4));
- if(EST>dlt) {S*=.6; continue;}
+ if(EST>dlt) {
+ S*=.6;
+ continue;
+ }
}
- // Parameters calculation
- //
- if((!ste && S0 > fabs(S)*100.) || fabs(P[45]+=S) > 2000.) return false;
+ /// Parameters calculation
+ /// if the step grew too small, or we traveled more than 2 meters, abort, this went badly wrong...
+ if((!ste && S0 > fabs(S)*100.) || fabs(globalPars[45]+=S) > 2000.) return false;
ste = true;
double A0arr[3]{A0,B0,C0};
@@ -1228,7 +1481,7 @@ bool InDet::SiTrajectoryElement_xk::rungeKuttaToPlane
double A6arr[3]{A6,B6,C6};
if(Jac) {
- Trk::propJacobian(P,H0,H1,H2,A,A0arr,A3arr,A4arr,A6arr,S3);
+ Trk::propJacobian(globalPars,H0,H1,H2,A,A0arr,A3arr,A4arr,A6arr,S3);
}
R[0]+=(A2+A3+A4)*S3;
@@ -1241,26 +1494,46 @@ bool InDet::SiTrajectoryElement_xk::rungeKuttaToPlane
double D = 1./sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]);
A[0]*=D; A[1]*=D; A[2]*=D;
- // New step estimation
- //
- double a = A[0]*m_Tr[6]+A[1]*m_Tr[7]+A[2]*m_Tr[8]; if(a==0.) return false;
- double Sn = ((m_Tr[12]-R[0]*m_Tr[6])-(R[1]*m_Tr[7]+R[2]*m_Tr[8]))/a;
+ /// New step estimation
+ ///
+ double a = A[0]*m_localTransform[6]+A[1]*m_localTransform[7]+A[2]*m_localTransform[8];
+ if(a==0.) return false;
+ double Sn = ((m_localTransform[12]-R[0]*m_localTransform[6])-(R[1]*m_localTransform[7]+R[2]*m_localTransform[8]))/a;
double aSn = fabs(Sn);
+ /// if the remaining step is below threshold,
+ /// we can finish with a small linear step
if(aSn <= Smin) {
- double Sl = 2./S; sA[0] = A6*Sl; sA[1] = B6*Sl; sA[2] = C6*Sl;
- R[0]+=(A[0]*Sn); R[1]+=(A[1]*Sn); R[2]+=(A[2]*Sn); P[45]+=Sn; return true;
+ double Sl = 2./S;
+ sA[0] = A6*Sl;
+ sA[1] = B6*Sl;
+ sA[2] = C6*Sl;
+
+ R[0]+=(A[0]*Sn);
+ R[1]+=(A[1]*Sn);
+ R[2]+=(A[2]*Sn);
+ globalPars[45]+=Sn;
+
+ return true;
}
double aS = fabs(S);
-
+
+ /// if we ran past the surface and changed signs more than twice, abort
if ( S*Sn < 0. ) {
if(++it > 2) return false;
- if (aSn < aS) S = Sn; else S =-S;
+ /// if absolute step size is decreasing, we keep going
+ if (aSn < aS) S = Sn;
+ /// if step increased, try the other direction
+ else S =-S;
}
+ /// otherwise refine step length if shorter than before
else if( aSn < aS ) S = Sn;
-
- f0[0]=f[0]; f0[1]=f[1]; f0[2]=f[2];
+
+ /// update field for next iteration
+ f0[0]=f[0];
+ f0[1]=f[1];
+ f0[2]=f[2];
}
return false;
}
@@ -1270,14 +1543,15 @@ bool InDet::SiTrajectoryElement_xk::rungeKuttaToPlane
/////////////////////////////////////////////////////////////////////////////////
bool InDet::SiTrajectoryElement_xk::straightLineStepToPlane
-(bool Jac,double* P)
+(bool Jac,double* globalPars)
{
- double* R = &P[ 0]; // Start coordinates
- double* A = &P[ 3]; // Start directions
-
- double a = A[0]*m_Tr[6]+A[1]*m_Tr[7]+A[2]*m_Tr[8]; if(a==0.) return false;
- double S = ((m_Tr[12]-R[0]*m_Tr[6])-(R[1]*m_Tr[7]+R[2]*m_Tr[8]))/a;
- P[45] = S;
+ double* R = &globalPars[ 0]; // Start coordinates
+ double* A = &globalPars[ 3]; // Start directions
+ ///
+ double a = A[0]*m_localTransform[6]+A[1]*m_localTransform[7]+A[2]*m_localTransform[8];
+ if(a==0.) return false;
+ double S = ((m_localTransform[12]-R[0]*m_localTransform[6])-(R[1]*m_localTransform[7]+R[2]*m_localTransform[8]))/a;
+ globalPars[45] = S;
// Track parameters in last point
//
@@ -1287,8 +1561,8 @@ bool InDet::SiTrajectoryElement_xk::straightLineStepToPlane
//
for(int i=7; i<42; i+=7) {
- double* dR = &P[i ];
- double* dA = &P[i+3];
+ double* dR = &globalPars[i ];
+ double* dA = &globalPars[i+3];
dR[0]+=(dA[0]*S); dR[1]+=(dA[1]*S); dR[2]+=(dA[2]*S);
}
return true;
@@ -1296,54 +1570,54 @@ bool InDet::SiTrajectoryElement_xk::straightLineStepToPlane
InDet::SiTrajectoryElement_xk::SiTrajectoryElement_xk()
{
- m_detstatus =-1 ;
- m_status = 0 ;
- m_nlinksF = 0 ;
- m_nlinksB = 0 ;
- m_ndist = 0 ;
- m_radlength = .03;
- m_radlengthN = .03;
- m_energylose = .4 ;
- m_tools = 0 ;
- m_noisemodel = 0 ;
+ m_detstatus =-1 ;
+ m_status = 0 ;
+ m_nlinksForward = 0 ;
+ m_nlinksBackward = 0 ;
+ m_ndist = 0 ;
+ m_radlength = .03;
+ m_radlengthN = .03;
+ m_energylose = .4 ;
+ m_tools = 0 ;
+ m_noisemodel = 0 ;
m_covariance.resize(2,2);
m_covariance<<0.,0.,0.,0.;
- m_ndf = 0 ;
- m_ndfF = 0 ;
- m_ndfB = 0 ;
- m_ntsos = 0 ;
- m_maxholes = 0 ;
- m_maxdholes = 0 ;
- m_xi2F = 0.;
- m_xi2B = 0.;
- m_xi2totalF = 0.;
- m_xi2totalB = 0.;
- m_halflength = 0.;
- m_step = 0.;
- m_xi2max = 0.;
- m_dist = 0.;
- m_xi2maxNoAdd = 0.;
- m_xi2maxlink = 0.;
- m_xi2multi = 0.;
- m_detelement = 0 ;
- m_detlink = 0 ;
- m_surface = 0 ;
- m_cluster = 0 ;
- m_clusterOld = 0 ;
- m_clusterNoAdd= 0 ;
- m_updatorTool = 0 ;
- m_proptool = 0 ;
- m_riotool = 0 ;
- m_inside = 0 ;
- m_nholesF = 0 ;
- m_nholesB = 0 ;
- m_dholesF = 0 ;
- m_dholesB = 0 ;
- m_nclustersF = 0 ;
- m_nclustersB = 0 ;
- m_npixelsB = 0 ;
- m_stereo = false;
- m_fieldMode = false;
+ m_ndf = 0 ;
+ m_ndfBackward = 0 ;
+ m_ndfForward = 0 ;
+ m_ntsos = 0 ;
+ m_maxholes = 0 ;
+ m_maxdholes = 0 ;
+ m_xi2Forward = 0.;
+ m_xi2Backward = 0.;
+ m_xi2totalForward = 0.;
+ m_xi2totalBackward = 0.;
+ m_halflength = 0.;
+ m_step = 0.;
+ m_xi2max = 0.;
+ m_dist = 0.;
+ m_xi2maxNoAdd = 0.;
+ m_xi2maxlink = 0.;
+ m_xi2multi = 0.;
+ m_detelement = 0 ;
+ m_detlink = 0 ;
+ m_surface = 0 ;
+ m_cluster = 0 ;
+ m_clusterOld = 0 ;
+ m_clusterNoAdd = 0 ;
+ m_updatorTool = 0 ;
+ m_proptool = 0 ;
+ m_riotool = 0 ;
+ m_inside = 0 ;
+ m_nholesForward = 0 ;
+ m_nholesBackward = 0 ;
+ m_dholesForward = 0 ;
+ m_dholesBackward = 0 ;
+ m_nclustersForward = 0 ;
+ m_nclustersBackward = 0 ;
+ m_npixelsBackward = 0 ;
+ m_stereo = false ;
+ m_fieldMode = false ;
m_tsos[0]=m_tsos[1]=m_tsos[2]=0;
}
@@ -1358,54 +1632,54 @@ InDet::SiTrajectoryElement_xk& InDet::SiTrajectoryElement_xk::operator =
{
if(&E==this) return(*this);
- m_fieldMode = E.m_fieldMode ;
- m_status = E.m_status ;
- m_detstatus = E.m_detstatus ;
- m_inside = E.m_inside ;
- m_ndist = E.m_ndist ;
- m_stereo = E.m_stereo ;
- m_detelement = E.m_detelement ;
- m_detlink = E.m_detlink ;
- m_surface = E.m_surface ;
- m_sibegin = E.m_sibegin ;
- m_siend = E.m_siend ;
- m_cluster = E.m_cluster ;
- m_clusterOld = E.m_clusterOld ;
- m_clusterNoAdd = E.m_clusterNoAdd;
- m_parametersPF = E.m_parametersPF;
- m_parametersUF = E.m_parametersUF;
- m_parametersPB = E.m_parametersPB;
- m_parametersUB = E.m_parametersUB;
- m_parametersSM = E.m_parametersSM;
- m_dist = E.m_dist ;
- m_xi2F = E.m_xi2F ;
- m_xi2B = E.m_xi2B ;
- m_xi2totalF = E.m_xi2totalF ;
- m_xi2totalB = E.m_xi2totalB ;
- m_radlength = E.m_radlength ;
- m_radlengthN = E.m_radlengthN ;
- m_energylose = E.m_energylose ;
- m_halflength = E.m_halflength ;
- m_step = E.m_step ;
- m_nlinksF = E.m_nlinksF ;
- m_nlinksB = E.m_nlinksB ;
- m_nholesF = E.m_nholesF ;
- m_nholesB = E.m_nholesB ;
- m_dholesF = E.m_dholesF ;
- m_dholesB = E.m_dholesB ;
- m_noisemodel = E.m_noisemodel ;
- m_ndf = E.m_ndf ;
- m_ndfF = E.m_ndfF ;
- m_ndfB = E.m_ndfB ;
- m_ntsos = E.m_ntsos ;
- m_nclustersF = E.m_nclustersF ;
- m_nclustersB = E.m_nclustersB ;
- m_npixelsB = E.m_npixelsB ;
- m_noise = E.m_noise ;
- m_tools = E.m_tools ;
- m_covariance = E.m_covariance ;
- for(int i=0; i!=m_nlinksF; ++i) {m_linkF[i]=E.m_linkF[i];}
- for(int i=0; i!=m_nlinksB; ++i) {m_linkB[i]=E.m_linkB[i];}
+ m_fieldMode = E.m_fieldMode ;
+ m_status = E.m_status ;
+ m_detstatus = E.m_detstatus ;
+ m_inside = E.m_inside ;
+ m_ndist = E.m_ndist ;
+ m_stereo = E.m_stereo ;
+ m_detelement = E.m_detelement ;
+ m_detlink = E.m_detlink ;
+ m_surface = E.m_surface ;
+ m_sibegin = E.m_sibegin ;
+ m_siend = E.m_siend ;
+ m_cluster = E.m_cluster ;
+ m_clusterOld = E.m_clusterOld ;
+ m_clusterNoAdd = E.m_clusterNoAdd;
+ m_parametersPredForward = E.m_parametersPredForward;
+ m_parametersUpdatedForward = E.m_parametersUpdatedForward;
+ m_parametersPredBackward = E.m_parametersPredBackward;
+ m_parametersUpdatedBackward = E.m_parametersUpdatedBackward;
+ m_parametersSM = E.m_parametersSM;
+ m_dist = E.m_dist ;
+ m_xi2Forward = E.m_xi2Forward ;
+ m_xi2Backward = E.m_xi2Backward ;
+ m_xi2totalForward = E.m_xi2totalForward ;
+ m_xi2totalBackward = E.m_xi2totalBackward;
+ m_radlength = E.m_radlength ;
+ m_radlengthN = E.m_radlengthN ;
+ m_energylose = E.m_energylose ;
+ m_halflength = E.m_halflength ;
+ m_step = E.m_step ;
+ m_nlinksForward = E.m_nlinksForward ;
+ m_nlinksBackward = E.m_nlinksBackward ;
+ m_nholesForward = E.m_nholesForward ;
+ m_nholesBackward = E.m_nholesBackward ;
+ m_dholesForward = E.m_dholesForward ;
+ m_dholesBackward = E.m_dholesBackward ;
+ m_noisemodel = E.m_noisemodel ;
+ m_ndf = E.m_ndf ;
+ m_ndfForward = E.m_ndfForward ;
+ m_ndfBackward = E.m_ndfBackward ;
+ m_ntsos = E.m_ntsos ;
+ m_nclustersForward = E.m_nclustersForward ;
+ m_nclustersBackward = E.m_nclustersBackward ;
+ m_npixelsBackward = E.m_npixelsBackward ;
+ m_noise = E.m_noise ;
+ m_tools = E.m_tools ;
+ m_covariance = E.m_covariance ;
+ for(int i=0; i!=m_nlinksForward; ++i) {m_linkForward[i]=E.m_linkForward[i];}
+ for(int i=0; i!=m_nlinksBackward; ++i) {m_linkBackward[i]=E.m_linkBackward[i];}
for(int i=0; i!=m_ntsos ; ++i) {m_tsos [i]=E.m_tsos [i];}
for(int i=0; i!=m_ntsos ; ++i) {m_utsos[i]=E.m_utsos [i];}
return(*this);
@@ -1460,15 +1734,15 @@ bool InDet::SiTrajectoryElement_xk::difference() const
bool InDet::SiTrajectoryElement_xk::isNextClusterHoleB(bool& cl,double& X)
{
cl = false ;
- X = m_xi2totalB-m_xi2B;
+ X = m_xi2totalBackward-m_xi2Backward;
- if(m_nlinksB > 1 && m_linkB[1].xi2() <= m_xi2max) {
- X+=m_linkB[1].xi2();
+ if(m_nlinksBackward > 1 && m_linkBackward[1].xi2() <= m_xi2max) {
+ X+=m_linkBackward[1].xi2();
cl = true; return true;
}
if(m_inside < 0) {
- if(m_nholesB < m_maxholes && m_dholesB < m_maxdholes) return true;
+ if(m_nholesBackward < m_maxholes && m_dholesBackward < m_maxdholes) return true;
return false;
}
return true;
@@ -1477,16 +1751,16 @@ bool InDet::SiTrajectoryElement_xk::isNextClusterHoleB(bool& cl,double& X)
bool InDet::SiTrajectoryElement_xk::isNextClusterHoleF(bool& cl,double& X)
{
cl = false ;
- X = m_xi2totalF-m_xi2F;
+ X = m_xi2totalForward-m_xi2Forward;
if(m_detstatus == 2) return false;
- if(m_nlinksF > 1 && m_linkF[1].xi2() <= m_xi2max) {
- X+=m_linkF[1].xi2();
+ if(m_nlinksForward > 1 && m_linkForward[1].xi2() <= m_xi2max) {
+ X+=m_linkForward[1].xi2();
cl = true; return true;
}
if(m_inside < 0) {
- if(m_nholesF < m_maxholes && m_dholesF < m_maxdholes) return true;
+ if(m_nholesForward < m_maxholes && m_dholesForward < m_maxdholes) return true;
return false;
}
return true;
@@ -1500,12 +1774,12 @@ void InDet::SiTrajectoryElement_xk::setCluster(const InDet::SiCluster* Cl)
void InDet::SiTrajectoryElement_xk::setParametersB(Trk::PatternTrackParameters& P)
{
- m_parametersUB = P;
+ m_parametersUpdatedBackward = P;
}
void InDet::SiTrajectoryElement_xk::setParametersF(Trk::PatternTrackParameters& P)
{
- m_parametersUF = P;
+ m_parametersUpdatedForward = P;
}
void InDet::SiTrajectoryElement_xk::setNdist(int n)
@@ -1521,11 +1795,12 @@ bool InDet::SiTrajectoryElement_xk::addCluster
(Trk::PatternTrackParameters& Ta,Trk::PatternTrackParameters& Tb,double& Xi2)
{
int N;
+ /// non-stereo measurements
if(!m_stereo) {
-
+ /// set m_covariance to the cluster errors
patternCovariances
(m_cluster,m_covariance(0,0),m_covariance(1,0),m_covariance(1,1));
-
+ /// update state
if(m_detelement->isSCT()) {
return m_updatorTool->addToStateOneDimension
(Ta,m_cluster->localPosition(),m_covariance,Tb,Xi2,N);
@@ -1533,6 +1808,7 @@ bool InDet::SiTrajectoryElement_xk::addCluster
return m_updatorTool->addToState
(Ta,m_cluster->localPosition(),m_covariance,Tb,Xi2,N);
}
+ /// for stereo, use the cluster local covariance directly
return m_updatorTool->addToStateOneDimension
(Ta,m_cluster->localPosition(),m_cluster->localCovariance(),Tb,Xi2,N);
}
@@ -1585,9 +1861,11 @@ void InDet::SiTrajectoryElement_xk::noiseInitiate()
///////////////////////////////////////////////////////////////////
bool InDet::SiTrajectoryElement_xk::initiateState
-(Trk::PatternTrackParameters& Ta,Trk::PatternTrackParameters& Tb)
+(Trk::PatternTrackParameters& inputPars,Trk::PatternTrackParameters& outputPars)
{
- return Tb.initiate(Ta,m_cluster->localPosition(),m_cluster->localCovariance());
+ /// Gives Tb the x and (for pix) the y of the cluster, and the corresponding cov elements,
+ /// and copies everything else from inputPars
+ return outputPars.initiate(inputPars,m_cluster->localPosition(),m_cluster->localCovariance());
}
///////////////////////////////////////////////////////////////////
@@ -1598,11 +1876,18 @@ void InDet::SiTrajectoryElement_xk::patternCovariances
(const InDet::SiCluster* c,double& covX,double& covXY,double& covY)
{
const Amg::MatrixX& v = c->localCovariance();
- covX = c->width().phiR(); covX*=(covX*.08333); if(covX < v(0,0)) covX=v(0,0);
+ covX = c->width().phiR();
+ covX*=(covX*.08333); /// sigma ~pitch / sqrt(12)
+ if(covX < v(0,0)) covX=v(0,0); /// if larger error in cluster covariance, replace covx by it
covXY = 0.;
- if(m_ndf==1) {covY=v(1,1);}
+ if(m_ndf==1) { /// strips: take y error from cov
+ covY=v(1,1);
+ }
else {
- covY=c->width().z(); covY*=(covY*.08333); if(covY < v(1,1)) covY=v(1,1);
+ /// pixels: again take either pitch/sqrt(12) or the cluster cov, whichever is larger
+ covY=c->width().z();
+ covY*=(covY*.08333);
+ if(covY < v(1,1)) covY=v(1,1);
}
}
diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectory_xk.cxx b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectory_xk.cxx
index 6d30898c9aed94aedd2763b068f84a3809e17832..12f281ce1f09a021b5b5219a85b476bc1254025b 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectory_xk.cxx
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectory_xk.cxx
@@ -303,7 +303,7 @@ std::ostream& InDet::SiTrajectory_xk::dump( std::ostream& out ) const
out<<"| Has"<<std::setw(3)<<m_nElements
<<" ("
- <<std::setw(3)<<m_naElements
+ <<std::setw(3)<<m_nActiveElements
<<")"
<<" elements and "
<<std::setw(2)<<m_nclusters+m_nclustersNoAdd<<" ("
@@ -312,9 +312,9 @@ std::ostream& InDet::SiTrajectory_xk::dump( std::ostream& out ) const
<<" |"
<<std::endl;
out<<"| Has number of holes before, inside, after and gap= "
- <<std::setw(2)<<m_nholesb
+ <<std::setw(2)<<m_nHolesBefore
<<std::setw(2)<<m_nholes
- <<std::setw(2)<<m_nholese
+ <<std::setw(2)<<m_nHolesAfter
<<std::setw(3)<<m_dholes<<" |"
<<std::endl;
out<<"| F B |"
@@ -476,7 +476,7 @@ std::ostream& InDet::SiTrajectory_xk::dump( std::ostream& out ) const
}
///////////////////////////////////////////////////////////////////
-// Build initiate trajectory
+// Initiate trajectory
///////////////////////////////////////////////////////////////////
bool InDet::SiTrajectory_xk::initialize
@@ -490,14 +490,15 @@ bool InDet::SiTrajectory_xk::initialize
bool & rquality ,
const EventContext & ctx )
{
+ /// reset state
m_nholes = 0;
- m_nholesb = 0;
- m_nholese = 0;
+ m_nHolesBefore = 0;
+ m_nHolesAfter = 0;
m_dholes = 0;
m_nclusters = 0;
m_nclustersNoAdd = 0;
m_nElements = 0;
- m_naElements = 0;
+ m_nActiveElements = 0;
m_firstElement = -100;
m_lastElement = 0;
m_ndfcut = 0;
@@ -506,133 +507,235 @@ bool InDet::SiTrajectory_xk::initialize
int ndfwrong = 0;
double Xi2cut = 2.*m_tools->xi2max();
- std::list<const InDet::SiCluster*>::iterator c;
+ std::list<const InDet::SiCluster*>::iterator iter_cluster;
if (lSiCluster.size() < 2) return false;
- std::vector<const InDet::SiDetElementBoundaryLink_xk*>::iterator r,re=DE.end();
+ std::vector<const InDet::SiDetElementBoundaryLink_xk*>::iterator iter_boundaryLink,endBoundaryLinks=DE.end();
int up = 0;
int last = 0;
- for (r=DE.begin(); r!=re; ++r) {
-
- const InDetDD::SiDetectorElement* de = (*r)->detElement();
- IdentifierHash id = de->identifyHash();
- const InDet::SiCluster* sic = nullptr;
-
- if (de->isPixel()) {
- if (PIX) {
- InDet::PixelClusterCollection::const_iterator sib, sie;
- const InDet::PixelClusterCollection *w = (*PIXc).indexFindPtr(id);
-
- if (w!=nullptr && w->begin()!=w->end()) {
-
- sib = w->begin(); sie = w->end();
-
- for (c=lSiCluster.begin(); c!=lSiCluster.end(); ++c) {
- if ((*c)->detectorElement()==de) {
- if (!m_nclusters) m_firstElement = m_nElements;
- else m_lastElement = m_nElements;
+ /// loop over all detector elements and assign the starting set of silicon clusters
+ /// at the right indices
+ for (iter_boundaryLink=DE.begin(); iter_boundaryLink!=endBoundaryLinks; ++iter_boundaryLink) {
+
+ /// get the associated detector element from the boundary link
+ const InDetDD::SiDetectorElement* detectorElement = (*iter_boundaryLink)->detElement();
+ IdentifierHash id = detectorElement->identifyHash();
+
+ /// book a placeholder for a potential cluster on this element
+ const InDet::SiCluster* theCluster = nullptr;
+
+ /// First case: Current element is a pixel module
+ if (detectorElement->isPixel()) {
+ if (PIX) { /// check if we are configured to use pixels!
+
+ InDet::PixelClusterCollection::const_iterator iter_PixelClusterColl, iter_PixelClusterCollEnd;
+ /// check for the pixel clusters on the given detector element, using the ID hash
+ const InDet::PixelClusterCollection *clustersOnElement = (*PIXc).indexFindPtr(id);
+ /// if we have any hits:
+ if (clustersOnElement!=nullptr && clustersOnElement->begin()!=clustersOnElement->end()) {
+
+ /// set iterators to the local cluster collection
+ iter_PixelClusterColl = clustersOnElement->begin();
+ iter_PixelClusterCollEnd = clustersOnElement->end();
+
+ /// Loop over the passed list with Si clusters to initiate the track with - these
+ /// are for example the clusters on our seed in inside-out- or backtracking.
+ for (iter_cluster=lSiCluster.begin(); iter_cluster!=lSiCluster.end(); ++iter_cluster) {
+ ///if this cluster is on the current element...
+ if ((*iter_cluster)->detectorElement()==detectorElement) {
+ /// if it is the first cluster we see, set the first element to the current index
+ if (m_nclusters==0){
+ m_firstElement = m_nElements;
+ }
+ /// otherwise, set the last element to the current index (will eventually point to the final cluster)
+ else{
+ m_lastElement = m_nElements;
+ }
+ /// increment cluster counter
++m_nclusters;
+ /// add 2 to the number of degrees of freedom counter (Pix is 2D)
m_ndfcut+=2;
- sic=(*c);
- lSiCluster.erase(c);
+ /// set our cluster pointer to point to this cluster
+ theCluster=(*iter_cluster);
+ /// remove the cluster from the list
+ lSiCluster.erase(iter_cluster);
+ /// and exit the loop over Si clusters. We can do this
+ /// because no two clusters are allowed to be on the same
+ /// detector element
break;
}
}
- m_elements[m_nElements].set(1,(*r),sib,sie,sic,ctx);
- ++m_naElements;
- } else if (m_naElements) {
- m_elements[m_nElements].set(0,(*r),sib,sie,sic,ctx);
- } else {
+ /// done, now we know if one of the existing clusters is on this element
+ /// set status = 1 (there are hits on this module), give it the boundary link and the space points on this element. Finally, also give it the cluster, if we found one.
+ m_elements[m_nElements].set(1,(*iter_boundaryLink),iter_PixelClusterColl,iter_PixelClusterCollEnd,theCluster,ctx);
+ /// and increment the counter of active (nonzero hits) elements
+ ++m_nActiveElements;
+ }
+ /// this branch is the case of a pixel module with no hits on it, if we have previously had an active element
+ else if (m_nActiveElements) {
+ /// here we set a status of 0.
+ m_elements[m_nElements].set(0,(*iter_boundaryLink),iter_PixelClusterColl,iter_PixelClusterCollEnd,theCluster,ctx);
+ }
+ /// this branch is taken if we have not yet found an active element and there are no hits on this module. No need to already start the trajectory!
+ else {
+
continue;
}
-
- m_elementsMap[m_nElements] = m_nElements; if (++m_nElements==300) break;
- }
- } else if (SCT) {
- InDet::SCT_ClusterCollection::const_iterator sib, sie;
- const InDet::SCT_ClusterCollection *w = (*SCTc).indexFindPtr(id);
-
- if (w!=nullptr && w->begin()!=w->end()) {
-
- sib = w->begin(); sie = w->end();
-
- for (c=lSiCluster.begin(); c!=lSiCluster.end(); ++c) {
- if ((*c)->detectorElement()==de) {
- if (!m_nclusters) m_firstElement = m_nElements;
- else m_lastElement = m_nElements;
+ /// map the index to itself. Always useful. Don't worry, it will get
+ /// more interesting later...
+ m_elementsMap[m_nElements] = m_nElements;
+ /// if we exceed the bounds of our array, bail out.
+ /// Also increment the detector element index
+ if (++m_nElements==300) break;
+ } /// end of check on pixel flag
+ } /// end of pixel case
+ /// case 2: Strip module
+ else if (SCT) {
+ InDet::SCT_ClusterCollection::const_iterator iter_stripClusterColl, iter_StripClusterCollEnd;
+ /// again, we fetch the clusters for this detector element
+ const InDet::SCT_ClusterCollection *clustersOnElement = (*SCTc).indexFindPtr(id);
+
+ /// do we have any?
+ if (clustersOnElement!=nullptr && clustersOnElement->begin()!=clustersOnElement->end()) {
+
+ iter_stripClusterColl = clustersOnElement->begin();
+ iter_StripClusterCollEnd = clustersOnElement->end();
+
+ ///
+ for (iter_cluster=lSiCluster.begin(); iter_cluster!=lSiCluster.end(); ++iter_cluster) {
+ if ((*iter_cluster)->detectorElement()==detectorElement) {
+ /// is this the first cluster we found?
+ if (m_nclusters==0){
+ /// then mark it as the first element
+ m_firstElement = m_nElements;
+ }
+ /// otherwise, mark it as the last element
+ else{
+ m_lastElement = m_nElements;
+ }
+ /// increment cluster counter
++m_nclusters;
+ /// for SCT, add one degree of freedom (1D measurement)
m_ndfcut+=1;
- sic=(*c);
- lSiCluster.erase(c);
+ /// and update the cluster pointer, before cleaning up
+ theCluster=(*iter_cluster);
+ lSiCluster.erase(iter_cluster);
+ /// remember - only one cluster per detector element is possible due to
+ /// upstream filtering. So we can exit when we found one.
break;
}
}
- m_elements[m_nElements].set(1,(*r),sib,sie,sic,ctx);
- ++m_naElements;
- } else if (m_naElements) {
- m_elements[m_nElements].set(0,(*r),sib,sie,sic,ctx);
- } else {
+ /// Now, set up the trajectory element (det status 1) as in the pixel case
+ m_elements[m_nElements].set(1,(*iter_boundaryLink),iter_stripClusterColl,iter_StripClusterCollEnd,theCluster,ctx);
+ /// and increment the active element count
+ ++m_nActiveElements;
+ }
+ /// branch if no clusters on module and previously seen active element
+ else if (m_nActiveElements) {
+ /// set an corresponding element to detstatus = 0
+ m_elements[m_nElements].set(0,(*iter_boundaryLink),iter_stripClusterColl,iter_StripClusterCollEnd,theCluster,ctx);
+ }
+ /// branch for no clusters and no active elements seen so far
+ else {
+ /// skip this one
continue;
}
-
+ /// update elements map
m_elementsMap[m_nElements] = m_nElements;
+ /// and array boundary checking (yuck!!)
+ /// Also increment the detector element index
if (++m_nElements==300) break;
- }
+ } /// end of SCT if-statement
+ /// if the element we are currently processing is the first one where we saw a cluster:
if (m_firstElement == m_nElements-1) {
- up = m_nElements-1;
- if (!m_elements[up].firstTrajectorElement(Tp)) return false;
- } else if (sic) {
+ up = m_nElements-1; /// set the upper populated point to the current index
+ if (!m_elements[up].firstTrajectorElement(Tp)) return false; /// and update the current trajectory
+ /// element with the track parameters
+ /// we obtained upstream for
+ /// the starting surface
+ }
+ /// if the element we are currently processing saw a cluster but is not the first one
+ else if (theCluster) {
+ /// run forward propagation from the last element with a cluster to this one
if (!m_elements[m_nElements-1].ForwardPropagationWithoutSearch(m_elements[up])) {
return false;
}
-
+ /// update index of last element that had a cluster to point to this one
up = m_nElements-1;
+ /// if the chi2 looks good for the forward extension step, update index of last good cluster
if (m_elements[m_nElements-1].xi2F() <= Xi2cut) {
last=up;
- } else {
+ }
+ /// if it does not look good
+ else {
+ /// add the ndf to "ndfwrong"
ndfwrong+=m_elements[m_nElements-1].ndf();
+ /// if we have collected more than 3 badly fitting DoF (2 pix or 1 Pix + 1 SCT or 3 SCT), bail out
if (ndfwrong > 3) return false;
+ /// reduce ndf for cut by the bad degrees of freedom
m_ndfcut-=m_elements[m_nElements-1].ndf();
+ /// reduce cluster count, don't include this track
--m_nclusters;
+ /// and erase the cluster again
m_elements[m_nElements-1].eraseClusterForwardPropagation();
}
}
+ /// Case if we already have clusters from the seed on track,
+ /// and some clusters left to put on it, but no seed cluster on this element
+ /// Corresponds to have a hole in the seed.
else if (m_nclusters && lSiCluster.size()) {
-
+ /// propagate to the current DE from the last one where we had a cluster from the seed
+ /// if the propagation fails
if (!m_elements[m_nElements-1].ForwardPropagationWithoutSearch(m_elements[up])) {
+ /// if we have a cluster here, something went wrong in the upstream logic
if (m_elements[m_nElements-1].cluster()) return false;
+ /// otherwise, remove this guy from consideration.
+ /// It will be overwritten by the next one we see
--m_nElements;
- if (m_elements[m_nElements-1].detstatus()) --m_naElements;
- } else {
+ /// also adapt the number of active elements if needed
+ if (m_elements[m_nElements-1].detstatus()) --m_nActiveElements;
+ }
+ /// if the propagation succeeded
+ else {
+ /// increment hole count if we expect a hit
if (m_elements[m_nElements-1].inside()<0) ++m_nholes;
+ /// update the index of the last element
up = m_nElements-1;
}
- }
- }
+ } /// end of case of no hit and expecting a hit
+ } /// end of loop over boundary links
+ /// did we manage to assign all our seed hits to elements on the road?
if (lSiCluster.size()) {
+ /// no? Then our search road was badly chosen. Return this info for logging and bail out
rquality = false;
return false;
}
+ /// if some seed hits are badly fitting and we have less than 6 remaining good DoF,
+ /// we give up
if (ndfwrong && m_ndfcut < 6) return false;
+ /// update degrees of freedom to the current count of good ones
m_ndf = m_ndfcut;
+ /// truncate the ndfcut variable to 6
if (m_ndfcut > 6) m_ndfcut = 6;
- // Kill empty trajectory elements in end of trajectory
- //
+ /// Kill empty trajectory elements at the end of our trajectory
int n = m_nElements-1;
+ /// count downwards
for (; n>0; --n) {
+ /// and if the element is active, stop looping
if (m_elements[n].detstatus()>=0) break;
}
+ /// the index where we aborted is the last one where we
+ /// found an active element.
m_nElements = n+1;
- // Find last detector element with clusters
- //
+ /// Find last detector element with clusters
for (; n>0; --n) {
if (m_elements[n].detstatus() == 1) {
m_elements[n].lastActive();
@@ -640,24 +743,33 @@ bool InDet::SiTrajectory_xk::initialize
}
}
- // Kill uncrossed detector elements
- //
+ /// Kill uncrossed detector elements
+ /// this repopulates the elementsMap we started to fill earlier
int m = m_firstElement+1;
m_lastElement = last ;
+ /// loop from the element after the one with the first hit to the one
+ /// with the last hit from the seed
for (n = m; n!=m_lastElement; ++n) {
InDet::SiTrajectoryElement_xk& En = m_elements[m_elementsMap[n]];
+ /// if we either have a cluster on track or at least cross the element,
+ /// plug it into the m_elementsMap at the next position
if (En.cluster() || En.inside() <= 0) m_elementsMap[m++] = m_elementsMap[n];
}
+ /// now m_elementsMap contains the interesting elements
+ /// if we kicked out some elements:
if (m!=n) {
+ /// update the index of the last element
m_lastElement = m;
- for (; n!=m_nElements; ++n) m_elementsMap[m++] = m_elementsMap[n];
+ /// then add the remaining ones beyond the last cluster back to the map
+ for (; n!=m_nElements; ++n){
+ m_elementsMap[m++] = m_elementsMap[n];
+ }
m_nElements = m;
}
if (!m_tools->bremNoise()) return true;
- // Change noise model for last trajectory elements
- //
+ /// If we run with brem, update noise model for last trajectory elements
for (n=m_lastElement; n!=m_nElements; ++n) {
m_elements[m_elementsMap[n]].bremNoiseModel();
}
@@ -681,18 +793,18 @@ bool InDet::SiTrajectory_xk::trackParametersToClusters
std::multimap<double,const InDet::SiCluster*> xi2cluster;
- std::vector<const InDet::SiDetElementBoundaryLink_xk*>::iterator r,re=DE.end();
+ std::vector<const InDet::SiDetElementBoundaryLink_xk*>::iterator iter_boundaryLink,endBoundaryLinks=DE.end();
std::multimap<const Trk::PrepRawData*,const Trk::Track*>::const_iterator t, te =PT.end();
double xi2Cut = .5;
int ndfCut = 6;
- for (r=DE.begin(); r!=re; ++r) {
+ for (iter_boundaryLink=DE.begin(); iter_boundaryLink!=endBoundaryLinks; ++iter_boundaryLink) {
- const InDetDD::SiDetectorElement* de = (*r)->detElement();
- IdentifierHash id = de->identifyHash();
+ const InDetDD::SiDetectorElement* detectorElement = (*iter_boundaryLink)->detElement();
+ IdentifierHash id = detectorElement->identifyHash();
- bool sct = de->isSCT();
+ bool sct = detectorElement->isSCT();
if (!sct) {
InDet::PixelClusterCollection::const_iterator sib, sie;
@@ -704,7 +816,7 @@ bool InDet::SiTrajectory_xk::trackParametersToClusters
} else {
continue;
}
- if (!m_elements[0].ForwardPropagationForClusterSeach(m_nElements,Tp,(*r),sib,sie)) return false;
+ if (!m_elements[0].ForwardPropagationForClusterSeach(m_nElements,Tp,(*iter_boundaryLink),sib,sie)) return false;
} else {
InDet::SCT_ClusterCollection::const_iterator sib, sie;
const InDet::SCT_ClusterCollection *w = (*SCTc).indexFindPtr(id);
@@ -715,7 +827,7 @@ bool InDet::SiTrajectory_xk::trackParametersToClusters
} else {
continue;
}
- if (!m_elements[0].ForwardPropagationForClusterSeach(m_nElements,Tp,(*r),sib,sie)) return false;
+ if (!m_elements[0].ForwardPropagationForClusterSeach(m_nElements,Tp,(*iter_boundaryLink),sib,sie)) return false;
}
for (int i=0; i!=m_elements[0].nlinksF(); ++i) {
@@ -762,7 +874,7 @@ bool InDet::SiTrajectory_xk::globalPositionsToClusters
std::multimap<const Trk::PrepRawData*,const Trk::Track*>& PT ,
std::list<const InDet::SiCluster*> & lSiCluster)
{
- std::vector<const InDet::SiDetElementBoundaryLink_xk*>::iterator r = DE.begin(), re = DE.end();
+ std::vector<const InDet::SiDetElementBoundaryLink_xk*>::iterator iter_boundaryLink = DE.begin(), endBoundaryLinks = DE.end();
std::list<Amg::Vector3D>::const_iterator g,gb = Gp.begin(), ge = Gp.end();
InDet::PixelClusterCollection::const_iterator pib, pie;
InDet::SCT_ClusterCollection::const_iterator sib, sie;
@@ -780,9 +892,9 @@ bool InDet::SiTrajectory_xk::globalPositionsToClusters
double xi2Cut = 10.;
m_ndf = 0 ;
- for (; r!=re; ++r) {
+ for (; iter_boundaryLink!=endBoundaryLinks; ++iter_boundaryLink) {
- const InDetDD::SiDetectorElement* d = (*r)->detElement();
+ const InDetDD::SiDetectorElement* d = (*iter_boundaryLink)->detElement();
IdentifierHash id = d->identifyHash ();
const Trk::Surface* su = &d->surface();
const Trk::PlaneSurface* pla = static_cast<const Trk::PlaneSurface*>(su);
@@ -828,8 +940,8 @@ bool InDet::SiTrajectory_xk::globalPositionsToClusters
Tp.setParametersWithCovariance(su,pv,cv);
- if (!sct) m_elements[0].CloseClusterSeach(Tp, (*r), pib, pie);
- else m_elements[0].CloseClusterSeach(Tp, (*r), sib, sie);
+ if (!sct) m_elements[0].CloseClusterSeach(Tp, (*iter_boundaryLink), pib, pie);
+ else m_elements[0].CloseClusterSeach(Tp, (*iter_boundaryLink), sib, sie);
const InDet::SiCluster* c = m_elements[0].cluster();
if (!c || m_elements[0].xi2F() > xi2Cut) continue;
if (sct) {
@@ -901,7 +1013,7 @@ bool InDet::SiTrajectory_xk::backwardSmoother(bool TWO)
n = m_elementsMap[ n ] ;
m_nclusters = m_elements[m].nclustersB() ;
m_nholes = m_elements[m].nholesB () ;
- m_nholesb = m_elements[n].nholesB ()-m_nholes;
+ m_nHolesBefore = m_elements[n].nholesB ()-m_nholes;
m_ndf = m_elements[m].ndfB() ;
if (m_ndf < m_ndfcut) return false;
@@ -912,14 +1024,14 @@ bool InDet::SiTrajectory_xk::backwardSmoother(bool TWO)
for (; m!=m_lastElement; ++m) {
InDet::SiTrajectoryElement_xk& Em = m_elements[m_elementsMap[m]];
- if (Em.inside() > 0) {if (Em.detstatus() > 0) --m_naElements; }
+ if (Em.inside() > 0) {if (Em.detstatus() > 0) --m_nActiveElements; }
else {m_elementsMap[n++] = m_elementsMap[m];}
}
m_lastElement = n;
// Test number trajector elemenst with clusters
//
- if (m_naElements < m_tools->clustersmin() && m_nholes+m_nholese) return false;
+ if (m_nActiveElements < m_tools->clustersmin() && m_nholes+m_nHolesAfter) return false;
if (n!=m) {
for (; m!=m_nElements; ++m) m_elementsMap[n++] = m_elementsMap[m];
m_nElements = n;
@@ -1072,7 +1184,7 @@ bool InDet::SiTrajectory_xk::backwardExtension(int itmax)
if (!nbest) return true;
m_nholes = hbest ;
- m_nholesb = hbestb;
+ m_nHolesBefore = hbestb;
m_nclusters += nbest ;
m_ndf += ndfbest;
m_firstElement = TE[0] ;
@@ -1114,48 +1226,67 @@ bool InDet::SiTrajectory_xk::backwardExtension(int itmax)
bool InDet::SiTrajectory_xk::forwardExtension(bool smoother,int itmax)
{
- const double pi2 = 2.*M_PI, pi = M_PI;
-
+ const double pi2 = 2.*M_PI;
+ const double pi = M_PI;
+
if (m_firstElement >= m_lastElement) return false;
- int L = m_lastElement, lElement = m_nElements-1;
+ /// last element with a known hit
+ int L = m_lastElement;
+ /// last element on the trajectory
+ int lastElementOnTraj = m_nElements-1;
+ /// smoothing step, if requested
if (smoother) {
L = 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[L]].firstTrajectorElement()) return false;
+ /// 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;
+ /// for all following elements with until the last with a known hit:
for (++L; L<=m_lastElement; ++L) {
-
- InDet::SiTrajectoryElement_xk& El = m_elements[m_elementsMap[L-1]];
- InDet::SiTrajectoryElement_xk& Ef = m_elements[m_elementsMap[L ]];
- if (!Ef.ForwardPropagationWithoutSearch(El)) return false;
+ InDet::SiTrajectoryElement_xk& previousElement = m_elements[m_elementsMap[L-1]];
+ InDet::SiTrajectoryElement_xk& thisElement = m_elements[m_elementsMap[L ]];
+ /// propagate forward the state from the previous element
+ if (!thisElement.ForwardPropagationWithoutSearch(previousElement)) return false;
}
}
- else {
-
+ /// this branch is entered if the first element is already in agreement between forward and back
+ /// so we already ran a smoothing step for the current hits on the trajectory
+ 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) {
- InDet::SiTrajectoryElement_xk& El = m_elements[m_elementsMap[L-1]];
- InDet::SiTrajectoryElement_xk& Ef = m_elements[m_elementsMap[L ]];
-
+ InDet::SiTrajectoryElement_xk& previousElement = m_elements[m_elementsMap[L-1]];
+ InDet::SiTrajectoryElement_xk& thisElement = m_elements[m_elementsMap[L ]];
+ /// this branch is entered while the range of hits we are looking at was still covered
+ /// by the previously made backward smoothing
if (!diff) {
- if ((diff = Ef.difference())) {
- if (!Ef.addNextClusterF(El,Ef.cluster())) return false;
+ /// update flag if we are still covered by the smoothing
+ diff = thisElement.difference();
+ /// if not (first element not covered):
+ if (diff) {
+ /// re-add cluster - this will update the counters and chi2 based on the refined upstream elements
+ if (!thisElement.addNextClusterF(previousElement,thisElement.cluster())) return false;
}
}
+ /// case if we have entered a piece of the tracks where there was no back smoothing yet
else {
- if (!Ef.ForwardPropagationWithoutSearch(El)) return false;
+ /// propagate forward
+ if (!thisElement.ForwardPropagationWithoutSearch(previousElement)) return false;
}
}
}
--L;
}
- if ( L== lElement) return true;
+ if ( L== lastElementOnTraj) return true;
// Search best forward trajectory prolongation
//
@@ -1193,7 +1324,7 @@ bool InDet::SiTrajectory_xk::forwardExtension(bool smoother,int itmax)
int p = F;
int Fs = F;
int Ml = M;
- int Cm = nclbest-lElement;
+ int Cm = nclbest-lastElementOnTraj;
bool h = false;
for (++F; F!=m_nElements; ++F) {
@@ -1238,7 +1369,7 @@ bool InDet::SiTrajectory_xk::forwardExtension(bool smoother,int itmax)
int m = m_elementsMap[l];
int nc = m_elements[m].nclustersF();
int nh = m_elements[m].nholesF();
- m_nholese = m_elements[m_elementsMap[p]].nholesF()-nh;
+ m_nHolesAfter = m_elements[m_elementsMap[p]].nholesF()-nh;
if (it==0 && nc==m_nclusters) return true;
@@ -1255,7 +1386,7 @@ bool InDet::SiTrajectory_xk::forwardExtension(bool smoother,int itmax)
itbest = it;
lbest = L ;
Xi2best = X ;
- ndbest = nd; if (fl==lElement && nd < ndcut) ndcut = nd;
+ ndbest = nd; if (fl==lastElementOnTraj && nd < ndcut) ndcut = nd;
for (int j=L+1; j<=Ml; ++j) {
@@ -1268,17 +1399,17 @@ bool InDet::SiTrajectory_xk::forwardExtension(bool smoother,int itmax)
}
}
nclbest = m_nclusters+nbest;
- if ( (nclbest >= 14 && !h) || (fl==lElement && ndbest == 0)) break;
+ if ( (nclbest >= 14 && !h) || (fl==lastElementOnTraj && ndbest == 0)) break;
}
- F = -1; bool cl = false; int nb = lElement-nclbest-1; double Xn;
+ F = -1; bool cl = false; int nb = lastElementOnTraj-nclbest-1; double Xn;
for (int j=Ml; j!=L; --j) {
int i = MP[j];
InDet::SiTrajectoryElement_xk& Ei = m_elements[m_elementsMap[i]];
- if (i!=lElement && Ei.cluster() && Ei.isNextClusterHoleF(cl,Xn)) {
+ if (i!=lastElementOnTraj && Ei.cluster() && Ei.isNextClusterHoleF(cl,Xn)) {
int nm = nb-i+Ei.nclustersF();
@@ -1299,7 +1430,7 @@ bool InDet::SiTrajectory_xk::forwardExtension(bool smoother,int itmax)
if (!nbest) return true;
m_nholes = hbest ;
- m_nholese = hbeste ;
+ m_nHolesAfter = hbeste ;
m_nclusters += nbest ;
m_ndf += ndfbest ;
m_lastElement = TE[nbest-1];
@@ -1484,7 +1615,7 @@ void InDet::SiTrajectory_xk::sortStep()
int e = m_elementsMap[n];
if (m_elements[e].cluster()) break;
if (m_elements[e].clusterNoAdd()) --m_nclustersNoAdd;
- else if (m_elements[e].inside() < 0 && m_elements[e].detstatus()>=0) {--m_nholes; ++m_nholesb;}
+ else if (m_elements[e].inside() < 0 && m_elements[e].detstatus()>=0) {--m_nholes; ++m_nHolesBefore;}
}
// Search last detector elements with cluster
@@ -1495,7 +1626,7 @@ void InDet::SiTrajectory_xk::sortStep()
int e = m_elementsMap[m];
if (m_elements[e].cluster()) break;
if (m_elements[e].clusterNoAdd()) --m_nclustersNoAdd;
- else if (m_elements[e].inside() < 0 && m_elements[e].detstatus()>=0) {--m_nholes; ++m_nholese;}
+ else if (m_elements[e].inside() < 0 && m_elements[e].detstatus()>=0) {--m_nholes; ++m_nHolesAfter;}
}
m_firstElement = n;
m_lastElement = m;
diff --git a/InnerDetector/InDetRecTools/SiCombinatorialTrackFinderTool_xk/src/SiCombinatorialTrackFinder_xk.cxx b/InnerDetector/InDetRecTools/SiCombinatorialTrackFinderTool_xk/src/SiCombinatorialTrackFinder_xk.cxx
index d7a3ca2f1a38843b8b8da18f21a11743042377dc..3b19b87a565247cd985b2db71c1c4c1968656936 100644
--- a/InnerDetector/InDetRecTools/SiCombinatorialTrackFinderTool_xk/src/SiCombinatorialTrackFinder_xk.cxx
+++ b/InnerDetector/InDetRecTools/SiCombinatorialTrackFinderTool_xk/src/SiCombinatorialTrackFinder_xk.cxx
@@ -296,13 +296,13 @@ void InDet::SiCombinatorialTrackFinder_xk::newEvent
newEvent(ctx, data);
data.trackinfo() = info;
- // Get track qulaity cuts information
- //
+ /// Get track quality cuts information from argument and write it into
+ /// the event data
getTrackQualityCuts(data, Cuts);
data.heavyIon() = false;
data.cosmicTrack() = 0;
-
+ /// update pattern recognition flags in the event data based on track info arg
if (info.patternRecoInfo(Trk::TrackInfo::SiSpacePointsSeedMaker_Cosmic)) {
data.cosmicTrack() = 1;
} else if (info.patternRecoInfo(Trk::TrackInfo::SiSpacePointsSeedMaker_HeavyIon)) {
@@ -345,7 +345,9 @@ const std::list<Trk::Track*>& InDet::SiCombinatorialTrackFinder_xk::getTracks
data.statistic().fill(false);
+ /// turn off brem fit & electron flags
data.tools().setBremNoise(false, false);
+ /// remove existing tracks
data.tracks().erase(data.tracks().begin(), data.tracks().end());
++data.inputseeds();
@@ -353,21 +355,23 @@ const std::list<Trk::Track*>& InDet::SiCombinatorialTrackFinder_xk::getTracks
return data.tracks();
}
- // Get track qulaity cuts information
- //
+ // Get track quality cuts information and write them into the event data... again?
getTrackQualityCuts(data, Cuts);
std::multimap<const Trk::PrepRawData*, const Trk::Track*> PT;
+ ///try to find the tracks
EStat_t FT = findTrack(data, Tp, Sp, Gp, DE, PT, ctx);
+ /// if we didn't find anything, bail out
if(FT!=Success) {
data.statistic()[FT] = true;
return data.tracks();
}
+
+ /// sort in step order
data.trajectory().sortStep();
- // Trk::Track production
- //
+
Trk::Track* t = convertToTrack(data);
++data.findtracks();
data.tracks().push_back(t);
@@ -550,15 +554,14 @@ InDet::SiCombinatorialTrackFinder_xk::EStat_t InDet::SiCombinatorialTrackFinder_
std::multimap<const Trk::PrepRawData*,const Trk::Track*>& PT,
const EventContext& ctx) const
{
+ /// init event data
if (not data.isInitialized()) initializeCombinatorialData(ctx, data);
- // List detector element links preparation
- //
+ /// populate a list of boundary links for the detector elements on our search road
std::vector<const InDet::SiDetElementBoundaryLink_xk*> DEL;
detectorElementLinks(DE, DEL,ctx);
- // Retrieve cached pointers to SG collections, or create the cache
- //
+ /// Retrieve cached pointers to SG collections, or create the cache
const InDet::PixelClusterContainer* p_pixcontainer = data.pixContainer();
if (m_usePIX && !p_pixcontainer) {
SG::ReadHandle<InDet::PixelClusterContainer> pixcontainer(m_pixcontainerkey,ctx);
@@ -572,77 +575,106 @@ InDet::SiCombinatorialTrackFinder_xk::EStat_t InDet::SiCombinatorialTrackFinder_
data.setSctContainer(p_sctcontainer);
}
- // List cluster preparation
- //
+ /// Cluster list preparationn
std::list<const InDet::SiCluster*> Cl;
- bool TWO = false;
+ bool isTwoPointSeed = false;
+ /// in inside-out track finding, Sp.size() is typically 3
+ /// for TRT-seeded backtracking, it is 2
+ /// both applications go into this branch
if (Sp.size() > 1) {
+
+ /// returns false if two clusters are on the same detector element
if (!spacePointsToClusters(Sp,Cl)) {
return TwoCluster;
}
- if (Sp.size()<=2) TWO = true;
- } else if (Gp.size() > 2) {
+ if (Sp.size()==2) isTwoPointSeed = true;
+ }
+ /// use case if we have a set of global positions rather than space points to start from
+ else if (Gp.size() > 2) {
if (!data.trajectory().globalPositionsToClusters(p_pixcontainer, p_sctcontainer, Gp, DEL, PT, Cl)) return TwoCluster;
} else {
+ /// use case if we have neither space-points nor global posittions, but track parameters to start from
if (!data.trajectory().trackParametersToClusters(p_pixcontainer, p_sctcontainer, Tp, DEL, PT, Cl)) return TwoCluster;
}
++data.goodseeds();
- // Build initial trajectory
- //
+ /// Build initial trajectory
bool Qr;
+ /// This will initialize the trajectory using the clusters we have and the parameter estimate
bool Q = data.trajectory().initialize(m_usePIX, m_useSCT, p_pixcontainer, p_sctcontainer, Tp, Cl, DEL, Qr,ctx);
+ /// if the initialisation fails (indicating this is probably a bad attempt) and we are running with
+ /// global positions instead of seeding
if (!Q && Sp.size() < 2 && Gp.size() > 3) {
-
+ /// reset our cluster list
Cl.clear();
+ /// try again using the clusters from the track parameters only
if (!data.trajectory().trackParametersToClusters(p_pixcontainer, p_sctcontainer, Tp, DEL, PT, Cl)) return TwoCluster;
+
if (!data.trajectory().initialize(m_usePIX, m_useSCT, p_pixcontainer, p_sctcontainer, Tp, Cl, DEL, Qr,ctx)) return TwoCluster;
+ /// if it worked now, set the quality flag to true
+
Q = Qr = true;
}
- if (!Qr){++data.roadbug(); return WrongRoad;}
+ /// this can never happen?!
+ if (!Qr){
+ ++data.roadbug();
+ return WrongRoad;
+ }
+ /// this can never happen either?!
if (!Q) return WrongInit;
+
++data.inittracks();
+ /// if the last cluster on track is in the pixels, this is assumed to come from a pixel seed
bool pixseed = data.trajectory().isLastPixel();
+ /// max #iterations
int itmax = 30;
if (data.simpleTrack()) itmax = 10;
if (data.heavyIon()) itmax = 50;
- // Track finding
- //
- if (pixseed) { // Strategy for pixel seeds
+ /// Track finding
+ if (pixseed) { /// Strategy for pixel seeds
if (!data.trajectory().forwardExtension (false,itmax)) return CantFindTrk;
if (!data.trajectory().backwardSmoother (false) ) return CantFindTrk;
if (!data.trajectory().backwardExtension(itmax) ) return CantFindTrk;
-
+ /// refine if needed
if (data.trajectory().difference() > 0) {
if (!data.trajectory().forwardFilter() ) return CantFindTrk;
if (!data.trajectory().backwardSmoother (false) ) return CantFindTrk;
}
int na = data.trajectory().nclustersNoAdd();
+ /// check if we found enough clusters
if (data.trajectory().nclusters()+na < data.nclusmin() || data.trajectory().ndf() < data.nwclusmin()) return CantFindTrk;
- } else { // Strategy for mixed seeds
- if (!data.trajectory().backwardSmoother(TWO) ) return CantFindTrk;
+ }
+ /// case of a strip seed or mixed PPS
+ else { // Strategy for mixed seeds
+ if (!data.trajectory().backwardSmoother(isTwoPointSeed) ) return CantFindTrk;
if (!data.trajectory().backwardExtension(itmax) ) return CantFindTrk;
if (!data.trajectory().forwardExtension(true,itmax)) return CantFindTrk;
+ /// first application of hit cut
int na = data.trajectory().nclustersNoAdd();
if (data.trajectory().nclusters()+na < data.nclusmin() || data.trajectory().ndf() < data.nwclusmin()) return CantFindTrk;
+ /// backward smooting
if (!data.trajectory().backwardSmoother(false) ) return CantFindTrk;
+ /// apply hit cut again following smoothing step
na = data.trajectory().nclustersNoAdd();
if (data.trajectory().nclusters()+na < data.nclusmin() || data.trajectory().ndf() < data.nwclusmin()) return CantFindTrk;
+ /// refine if needed
if (data.trajectory().difference() > 0) {
if (!data.trajectory().forwardFilter() ) return CantFindTrk;
if (!data.trajectory().backwardSmoother (false)) return CantFindTrk;
}
}
-
+ /// quality cut
if (data.trajectory().qualityOptimization() < (m_qualityCut*data.nclusmin()) ) return CantFindTrk;
+
if (data.trajectory().pTfirst () < data.pTmin() && data.trajectory().nclusters() < data.nclusmin() ) return CantFindTrk;
+
if (data.trajectory().nclusters() < data.nclusminb() || data.trajectory().ndf () < data.nwclusmin()) return CantFindTrk;
return Success;
@@ -700,12 +732,16 @@ void InDet::SiCombinatorialTrackFinder_xk::magneticFieldInit()
bool InDet::SiCombinatorialTrackFinder_xk::spacePointsToClusters
(const std::vector<const Trk::SpacePoint*>& Sp, std::list<const InDet::SiCluster*>& Sc) const
{
+ /// loop over all SP
for (const Trk::SpacePoint* s: Sp) {
+ /// get the first cluster on an SP
const Trk::PrepRawData* p = s->clusterList().first;
if (p) {
+ /// add to list
const InDet::SiCluster* c = static_cast<const InDet::SiCluster*>(p);
if (c) Sc.push_back(c);
}
+ /// for strips, also make sure to pick up the second one!
p = s->clusterList().second;
if (p) {
const InDet::SiCluster* c = static_cast<const InDet::SiCluster*>(p);
@@ -713,18 +749,18 @@ bool InDet::SiCombinatorialTrackFinder_xk::spacePointsToClusters
}
}
- // Detector elments test
- //
- std::list<const InDet::SiCluster*>::iterator c = Sc.begin(), cn, ce = Sc.end();
+ /// Detector elments test
+ std::list<const InDet::SiCluster*>::iterator cluster = Sc.begin(), nextCluster, endClusters = Sc.end();
- for (; c!=ce; ++c) {
+ /// here we reject cases where two subsequent clusters are on the same detector element
+ for (; cluster!=endClusters; ++cluster) {
- const InDetDD::SiDetectorElement* de = (*c)->detectorElement();
+ const InDetDD::SiDetectorElement* de = (*cluster)->detectorElement();
- cn = c;
- ++cn;
- for (; cn!=ce; ++cn) {
- if (de == (*cn)->detectorElement()) return false;
+ nextCluster = cluster;
+ ++nextCluster;
+ for (; nextCluster!=endClusters; ++nextCluster) {
+ if (de == (*nextCluster)->detectorElement()) return false;
}
}
@@ -829,8 +865,8 @@ void InDet::SiCombinatorialTrackFinder_xk::getTrackQualityCuts
void InDet::SiCombinatorialTrackFinder_xk::initializeCombinatorialData(const EventContext& ctx, SiCombinatorialTrackFinderData_xk& data) const {
- // Add conditions object to SiCombinatorialTrackFinderData to be able to access the field cache for each new event
- // Get conditions object for field cache
+ /// Add conditions object to SiCombinatorialTrackFinderData to be able to access the field cache for each new event
+ /// Get conditions object for field cache
SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCondObjInputKey, ctx};
const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
if (fieldCondObj == nullptr) {
@@ -839,7 +875,7 @@ void InDet::SiCombinatorialTrackFinder_xk::initializeCombinatorialData(const Eve
}
data.setFieldCondObj(fieldCondObj);
- // Must have set fieldCondObj BEFORE calling setTools because fieldCondObj is used there
+ /// Must have set fieldCondObj BEFORE calling setTools because fieldCondObj is used there
data.setTools(&*m_proptool,
&*m_updatortool,
&*m_riocreator,
diff --git a/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/SiTrackMakerTool_xk/SiTrackMaker_xk.h b/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/SiTrackMakerTool_xk/SiTrackMaker_xk.h
index d33d9e3866807fa4ba85d650e64d2add5484e9f1..2b4a4a100ee747966cf4b5941914370942e450dd 100644
--- a/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/SiTrackMakerTool_xk/SiTrackMaker_xk.h
+++ b/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/SiTrackMakerTool_xk/SiTrackMaker_xk.h
@@ -176,21 +176,21 @@ namespace InDet{
/////////////////////////////////////////////////////////////////////
mutable std::mutex m_counterMutex;
- mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_totalInputSeeds ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<double>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_totalUsedSeeds ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_totalNoTrackPar ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_totalBremSeeds ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_twoClusters ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_wrongRoad ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_wrongInit ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_noTrack ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_notNewTrack ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_bremAttempt ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_outputTracks ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_extraTracks ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_bremTracks ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_seedsWithTrack ATLAS_THREAD_SAFE;
- mutable std::array<std::atomic<double>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_deSize ATLAS_THREAD_SAFE;
+ mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_totalInputSeeds ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<double>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_totalUsedSeeds ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_totalNoTrackPar ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_totalBremSeeds ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_twoClusters ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_wrongRoad ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_wrongInit ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_noTrack ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_notNewTrack ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_bremAttempt ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_outputTracks ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_extraTracks ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_bremTracks ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<int>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_seedsWithTrack ATLAS_THREAD_SAFE {};
+ mutable std::array<std::atomic<double>,SiCombinatorialTrackFinderData_xk::kNSeedTypes> m_deSize ATLAS_THREAD_SAFE {};
mutable std::vector<std::vector<double>> m_usedSeedsEta ATLAS_THREAD_SAFE;
mutable std::vector<std::vector<double>> m_seedsWithTracksEta ATLAS_THREAD_SAFE;
@@ -255,6 +255,14 @@ namespace InDet{
MsgStream& dumpStatistics(MsgStream &out) const;
MsgStream& dumpconditions(MsgStream& out) const;
MsgStream& dumpevent(SiTrackMakerEventData_xk& data, MsgStream& out) const;
+
+ /// helper for working with the stat arrays
+ template <typename T, size_t N,size_t M> void resetCounter(std::array<std::array<T,M>,N> & a) const{
+ for (auto & subarr : a) resetCounter(subarr);
+ }
+ template <typename T, size_t N> void resetCounter(std::array<T,N> & a) const{
+ std::fill(a.begin(),a.end(),0);
+ }
};
} // end of name space
diff --git a/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/src/SiTrackMaker_xk.cxx b/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/src/SiTrackMaker_xk.cxx
index efc016c513edf12dbb006642654ead7be817d083..ac26da25831ace94f4454b8eb03941a7b3a60bc9 100644
--- a/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/src/SiTrackMaker_xk.cxx
+++ b/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/src/SiTrackMaker_xk.cxx
@@ -27,24 +27,7 @@
InDet::SiTrackMaker_xk::SiTrackMaker_xk
(const std::string& t,const std::string& n,const IInterface* p)
- : base_class(t, n, p),
- m_totalInputSeeds{ },
- m_totalUsedSeeds{ },
- m_totalNoTrackPar{ },
- m_totalBremSeeds{ },
- m_twoClusters{ },
- m_wrongRoad{ },
- m_wrongInit{ },
- m_noTrack{ },
- m_notNewTrack{ },
- m_bremAttempt{ },
- m_outputTracks{ },
- m_extraTracks{ },
- m_bremTracks{ },
- m_seedsWithTrack{ },
- m_deSize{ },
- m_usedSeedsEta( SiCombinatorialTrackFinderData_xk::kNSeedTypes, std::vector<double>(SiCombinatorialTrackFinderData_xk::kNRapidityRanges, 0.) ),
- m_seedsWithTracksEta( SiCombinatorialTrackFinderData_xk::kNSeedTypes, std::vector<double>(SiCombinatorialTrackFinderData_xk::kNRapidityRanges, 0.) )
+ : base_class(t, n, p)
{
}
@@ -55,18 +38,19 @@ InDet::SiTrackMaker_xk::SiTrackMaker_xk
StatusCode InDet::SiTrackMaker_xk::initialize()
{
- // Get beam geometry
- //
+ /// Get beam geometry
+ ///
if (not m_beamSpotKey.empty()) {
ATH_CHECK( m_beamSpotKey.initialize() );
}
+ /// read the config string for the field mode
if (m_fieldmode == "NoField") m_fieldModeEnum = Trk::NoField;
- else if (m_fieldmode == "MapSolenoid") m_fieldModeEnum = Trk::FastField;
+ else if (m_fieldmode == "MapSolenoid") m_fieldModeEnum = Trk::FastField; /// this one is the default
else m_fieldModeEnum = Trk::FullField;
- // Get detector elements road maker tool
- //
+ /// Get detector elements road maker tool
+ ///
if ( m_roadmaker.retrieve().isFailure() ) {
ATH_MSG_FATAL( "Failed to retrieve tool " << m_roadmaker );
return StatusCode::FAILURE;
@@ -74,8 +58,8 @@ StatusCode InDet::SiTrackMaker_xk::initialize()
ATH_MSG_INFO( "Retrieved tool " << m_roadmaker );
}
- // Get combinatorial track finder tool
- //
+ /// Get combinatorial track finder tool
+ ///
if ( m_tracksfinder.retrieve().isFailure() ) {
ATH_MSG_FATAL( "Failed to retrieve tool " << m_tracksfinder );
return StatusCode::FAILURE;
@@ -83,8 +67,9 @@ StatusCode InDet::SiTrackMaker_xk::initialize()
ATH_MSG_INFO( "Retrieved tool " << m_tracksfinder );
}
- // Get seed to track conversion tool
- //
+ /// Get seed to track conversion tool
+ /// This is used if we want to write out the seeds for
+ /// performance studies
if (m_seedsegmentsWrite) {
if (m_seedtrack.retrieve().isFailure()) {
ATH_MSG_FATAL( "Failed to retrieve tool " << m_seedtrack );
@@ -97,6 +82,7 @@ StatusCode InDet::SiTrackMaker_xk::initialize()
m_seedtrack.disable();
}
+ /// flag for HI running
m_heavyion = false;
@@ -121,31 +107,37 @@ StatusCode InDet::SiTrackMaker_xk::initialize()
m_trackinfo.setPatternRecognitionInfo(Trk::TrackInfo::SiSPSeededFinder );
}
+ /// this overrides the seed filter property for cosmic reco, in case it is not already set to 3 in the config.
+ /// For p-p, we usually operate either 1 (LRT, some special configs) or 2 (standard inside-out)
if (m_cosmicTrack) m_seedsfilter = 3;
+ /// this is on by default in offline tracking
ATH_CHECK(m_caloCluster.initialize(m_useBremModel && m_useCaloSeeds));
+ /// useSSSfilter is on by default in offline, m_useHClusSeed is off by default.
ATH_CHECK(m_caloHad.initialize( !m_useSSSfilter && m_useHClusSeed) );
+ /// pt cut can never be below 20 MeV
if (m_pTmin < 20.) m_pTmin = 20.;
-
- for (int i=0; i!=SiCombinatorialTrackFinderData_xk::kNSeedTypes; ++i) {
- m_totalInputSeeds[i] = 0;
- m_totalNoTrackPar[i] = 0;
- m_totalUsedSeeds[i] = 0.;
- m_outputTracks[i] = 0;
- m_totalBremSeeds[i] = 0;
- m_bremTracks[i] = 0;
- m_twoClusters[i] = 0;
- m_wrongRoad[i] = 0;
- m_wrongInit[i] = 0;
- m_noTrack[i] = 0;
- m_notNewTrack[i] = 0;
- m_bremAttempt[i] = 0;
- m_seedsWithTrack[i] = 0;
- m_deSize[i] = 0.;
- for(int j=0; j!=SiCombinatorialTrackFinderData_xk::kNRapidityRanges ;++j) { m_usedSeedsEta[i][j]=0.; m_seedsWithTracksEta[i][j]=0.;}
- }
+ /// initialize counters
+ resetCounter(m_totalInputSeeds);
+ resetCounter(m_totalNoTrackPar);
+ resetCounter(m_totalUsedSeeds);
+ resetCounter(m_outputTracks);
+ resetCounter(m_totalBremSeeds);
+ resetCounter(m_bremTracks);
+ resetCounter(m_twoClusters);
+ resetCounter(m_wrongRoad);
+ resetCounter(m_wrongInit);
+ resetCounter(m_noTrack);
+ resetCounter(m_notNewTrack);
+ resetCounter(m_bremAttempt);
+ resetCounter(m_seedsWithTrack);
+ resetCounter(m_deSize);
+ m_usedSeedsEta.resize(SiCombinatorialTrackFinderData_xk::kNSeedTypes);
+ m_seedsWithTracksEta.resize(SiCombinatorialTrackFinderData_xk::kNSeedTypes);
+ std::fill(m_usedSeedsEta.begin(), m_usedSeedsEta.end(), std::vector<double>(SiCombinatorialTrackFinderData_xk::kNRapidityRanges, 0.));
+ std::fill(m_seedsWithTracksEta.begin(), m_seedsWithTracksEta.end(), std::vector<double>(SiCombinatorialTrackFinderData_xk::kNRapidityRanges, 0.));
////////////////////////////////////////////////////////////////////////////////
ATH_CHECK( m_fieldCondObjInputKey.initialize());
@@ -470,6 +462,8 @@ MsgStream& InDet::SiTrackMaker_xk::dumpevent(SiTrackMakerEventData_xk& data, Msg
void InDet::SiTrackMaker_xk::newEvent(const EventContext& ctx, SiTrackMakerEventData_xk& data, bool PIX, bool SCT) const
{
+
+ /// initialize beam position
data.xybeam()[0] = 0.;
data.xybeam()[1] = 0.;
if (not m_beamSpotKey.empty()) {
@@ -479,33 +473,33 @@ void InDet::SiTrackMaker_xk::newEvent(const EventContext& ctx, SiTrackMakerEvent
data.xybeam()[1] = beamSpotHandle->beamPos()[1];
}
}
-
+ /// propagate pixel / strip usage to the event data object
data.pix() = PIX and m_usePix;
data.sct() = SCT and m_useSct;
- bool simpleTrack = false;
- InDet::TrackQualityCuts trackquality = setTrackQualityCuts(simpleTrack);
+ /// build a holder for the configured track quality cuts
+ InDet::TrackQualityCuts trackquality = setTrackQualityCuts(false);
- // New event for track finder tool
- //
+ /// Setup New event for track finder tool
m_tracksfinder->newEvent(ctx, data.combinatorialData(), m_trackinfo, trackquality);
- // Erase cluster to track association
- //
+ /// Erase cluster to track association
+ /// m_seedsfilter is non-zero in the vast majority of all applications,
+ /// so this is usually done
if (m_seedsfilter) data.clusterTrack().clear();
- // Erase statistic information
- //
+ /// Erase statistic information
+ ///
data.inputseeds() = 0;
data.goodseeds() = 0;
data.findtracks() = 0;
- for(int i=0; i!=SiCombinatorialTrackFinderData_xk::kNStatAllTypes; ++i) { for(int k = 0; k!=SiCombinatorialTrackFinderData_xk::kNSeedTypes; ++k) data.summaryStatAll()[i][k]; }
- for(int i=0; i!=SiCombinatorialTrackFinderData_xk::kNStatEtaTypes; ++i) { for(int k = 0; k!=SiCombinatorialTrackFinderData_xk::kNSeedTypes; ++k) { for(int r=0; r!=SiCombinatorialTrackFinderData_xk::kNRapidityRanges; ++r) data.summaryStatUsedInTrack()[i][k][r]; } }
+ resetCounter(data.summaryStatAll());
+ resetCounter(data.summaryStatUsedInTrack());
- // Retrieve
- //
+ /// retrieve calo seeds for brem fit
if (m_useBremModel && m_useCaloSeeds) {
+ /// reset old event data entries if any
data.caloF().clear();
data.caloR().clear();
data.caloZ().clear();
@@ -513,6 +507,7 @@ void InDet::SiTrackMaker_xk::newEvent(const EventContext& ctx, SiTrackMakerEvent
if (!m_caloCluster.key().empty()) {
SG::ReadHandle<CaloClusterROI_Collection> calo_cluster(m_caloCluster, ctx);
if (calo_cluster.isValid()) {
+ /// loop over the cluster ROIs and write them into the event data
for (const Trk::CaloClusterROI *c : * calo_cluster) {
data.caloF().push_back( c->globalPosition().phi ());
data.caloR().push_back( c->globalPosition().perp());
@@ -522,7 +517,9 @@ void InDet::SiTrackMaker_xk::newEvent(const EventContext& ctx, SiTrackMakerEvent
}
}
+ /// retrieve hadronic seeds for SSS seed filter
if (!m_useSSSfilter && m_useHClusSeed) {
+ /// reset old event data entries if any
data.hadF().clear();
data.hadR().clear();
data.hadZ().clear();
@@ -531,6 +528,7 @@ void InDet::SiTrackMaker_xk::newEvent(const EventContext& ctx, SiTrackMakerEvent
SG::ReadHandle<CaloClusterROI_Collection> calo_had(m_caloHad, ctx);
if (calo_had.isValid()) {
for (const Trk::CaloClusterROI *c : * calo_had) {
+ /// loop over the cluster ROIs and write them into the event data
data.hadF().push_back( c->globalPosition().phi ());
data.hadR().push_back( c->globalPosition().perp());
data.hadZ().push_back( c->globalPosition().z ());
@@ -538,6 +536,7 @@ void InDet::SiTrackMaker_xk::newEvent(const EventContext& ctx, SiTrackMakerEvent
}
}
}
+ /// if we want to write out the seeds, also call newEvent for the seed-to-track converter
if (m_seedsegmentsWrite) m_seedtrack->newEvent(data.conversionData(), m_trackinfo, m_patternName);
}
@@ -577,11 +576,10 @@ void InDet::SiTrackMaker_xk::newTrigEvent(const EventContext& ctx, SiTrackMakerE
void InDet::SiTrackMaker_xk::endEvent(SiTrackMakerEventData_xk& data) const
{
- // End event for track finder tool
- //
+ /// End event for track finder tool
m_tracksfinder->endEvent(data.combinatorialData());
- //correction to exclude memory fragmentation
+ /// correction to exclude memory fragmentation
data.clusterTrack().clear();
// end event for seed to track tool
@@ -598,22 +596,22 @@ void InDet::SiTrackMaker_xk::endEvent(SiTrackMakerEventData_xk& data) const
}
}
}
- for (int K = 0; K != SiCombinatorialTrackFinderData_xk::kNSeedTypes; ++K) {
- m_totalInputSeeds[K] += data.summaryStatAll()[kTotalInputSeeds][K];
- m_totalUsedSeeds[K] = m_totalUsedSeeds[K] + data.summaryStatAll()[kTotalUsedSeeds][K];
- m_totalNoTrackPar[K] += data.summaryStatAll()[kTotalNoTrackPar][K];
- m_totalBremSeeds[K] += data.summaryStatAll()[kTotalBremSeeds][K];
- m_twoClusters[K] += data.summaryStatAll()[kTwoClusters][K];
- m_wrongRoad[K] += data.summaryStatAll()[kWrongRoad][K];
- m_wrongInit[K] += data.summaryStatAll()[kWrongInit][K];
- m_noTrack[K] += data.summaryStatAll()[kNoTrack][K];
- m_notNewTrack[K] += data.summaryStatAll()[kNotNewTrack][K];
- m_bremAttempt[K] += data.summaryStatAll()[kBremAttempt][K];
- m_outputTracks[K] += data.summaryStatAll()[kOutputTracks][K];
- m_extraTracks[K] += data.summaryStatAll()[kExtraTracks][K];
- m_bremTracks[K] += data.summaryStatAll()[kBremTracks][K];
- m_seedsWithTrack[K] += data.summaryStatAll()[kSeedsWithTracks][K];
- m_deSize[K] = m_deSize[K] + data.summaryStatAll()[kDESize][K];
+ for (int K = 0; K != SiCombinatorialTrackFinderData_xk::kNSeedTypes; ++K) {
+ m_totalInputSeeds[K] += data.summaryStatAll()[kTotalInputSeeds][K];
+ m_totalUsedSeeds[K] = m_totalUsedSeeds[K] + data.summaryStatAll()[kTotalUsedSeeds][K];
+ m_totalNoTrackPar[K] += data.summaryStatAll()[kTotalNoTrackPar][K];
+ m_totalBremSeeds[K] += data.summaryStatAll()[kTotalBremSeeds][K];
+ m_twoClusters[K] += data.summaryStatAll()[kTwoClusters][K];
+ m_wrongRoad[K] += data.summaryStatAll()[kWrongRoad][K];
+ m_wrongInit[K] += data.summaryStatAll()[kWrongInit][K];
+ m_noTrack[K] += data.summaryStatAll()[kNoTrack][K];
+ m_notNewTrack[K] += data.summaryStatAll()[kNotNewTrack][K];
+ m_bremAttempt[K] += data.summaryStatAll()[kBremAttempt][K];
+ m_outputTracks[K] += data.summaryStatAll()[kOutputTracks][K];
+ m_extraTracks[K] += data.summaryStatAll()[kExtraTracks][K];
+ m_bremTracks[K] += data.summaryStatAll()[kBremTracks][K];
+ m_seedsWithTrack[K] += data.summaryStatAll()[kSeedsWithTracks][K];
+ m_deSize[K] = m_deSize[K] + data.summaryStatAll()[kDESize][K];
}
// Print event information
@@ -631,26 +629,35 @@ void InDet::SiTrackMaker_xk::endEvent(SiTrackMakerEventData_xk& data) const
std::list<Trk::Track*> InDet::SiTrackMaker_xk::getTracks
(const EventContext& ctx, SiTrackMakerEventData_xk& data, const std::vector<const Trk::SpacePoint*>& Sp) const
{
+ /// incremenet seed counter
++data.inputseeds();
+ /// 0 or 3 typically (PPP or SSS), other numbers indicate PPS/PSS (1/2)
int K = kindSeed(Sp);
+ /// eta of the seed, rounded down to leading digit via int-conversion
int r = rapidity(Sp);
+ /// more counter incrementation
++data.summaryStatAll()[kTotalInputSeeds][K];
+ /// prepare output list
std::list<Trk::Track*> tracks;
+ /// if we run the SI track maker without using the Si, this becomes a trivial task...
if (!data.pix() && !data.sct()) return tracks;
- bool good;
- !m_seedsfilter ? good=true : good=newSeed(data, Sp);
-
- if (!good) return tracks;
+ /// check seed quality.
+ bool isGoodSeed{true};
+ /// this checks if all of the clusters on the seed are already on one single existing track.
+ /// If not, we consider the seed to be "good"
+ if (m_seedsfilter) isGoodSeed=newSeed(data, Sp);
+ if (!isGoodSeed) return tracks;
data.dbm() = isDBMSeeds(*Sp.begin());
// Get AtlasFieldCache
MagField::AtlasFieldCache fieldCache;
+ /// read the B-field cache
SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCondObjInputKey, ctx};
const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
if (fieldCondObj == nullptr) {
@@ -659,61 +666,78 @@ std::list<Trk::Track*> InDet::SiTrackMaker_xk::getTracks
}
fieldCondObj->getInitializedCache (fieldCache);
- // Get initial parameters estimation
- //
- bool sss = false;
+ /// Get initial parameters estimation from our seed
const Trk::TrackParameters* Tp = nullptr;
if (data.dbm()) Tp = getAtaPlaneDBM(fieldCache, data, Sp);
- else Tp = getAtaPlane(fieldCache, data, sss && m_useHClusSeed, Sp);
- if (!Tp)
- {
+ else Tp = getAtaPlane(fieldCache, data, false, Sp);
+ /// if we failed to get the initial parameters, we bail out.
+ /// Can happen in certain pathological cases (e.g. malformed strip hits),
+ /// or if we would be running with calo-ROI strip seeds (we aren't)
+ if (!Tp) {
++data.summaryStatAll()[kTotalNoTrackPar][K];
return tracks;
}
+ /// otherwise, increment the 'good seeds' counter
++data.goodseeds();
- // Get detector elements road
- //
+ /// Now, obtain a search road of detector elements.
+ /// This is done by extrapolating our estimated starting parameters through the detector
+ /// and collecting all detector elements reasonably close to the projected trajectory.
+ /// This will populate the 'DE" list.
std::list<const InDetDD::SiDetectorElement*> DE;
if (!m_cosmicTrack) m_roadmaker->detElementsRoad(ctx, fieldCache, *Tp,Trk::alongMomentum, DE, data.roadMakerData());
else m_roadmaker->detElementsRoad(ctx, fieldCache, *Tp,Trk::oppositeMomentum,DE, data.roadMakerData());
+ /// if we don't use all of pix and SCT, filter our list, erasing any that don't fit our requirements
if (!data.pix() || !data.sct() || data.dbm()) detectorElementsSelection(data, DE);
+ /// if we did not find sufficient detector elements to fulfill the minimum cluster requirement,
+ /// bail out. We will not be able to build a track satisfying the cuts.
if ( static_cast<int>(DE.size()) < m_nclusmin) {
delete Tp;
return tracks;
}
+ /// update statistics tables - we have sufficient detector elements to have a chance of finding a track!
data.summaryStatAll()[kDESize][K] += double(DE.size());
++data.summaryStatAll()[kTotalUsedSeeds][K];
+ /// prepare a list of global positions
std::list<Amg::Vector3D> Gp;
+ /// update another counter
++data.summaryStatUsedInTrack()[kUsedSeedsEta][K][r];
- // Find possible list of tracks using space points space points information
- //
+ /// Find possible list of tracks using space points space points information
+ ///
if (!m_useBremModel) {
tracks = m_tracksfinder->getTracks (data.combinatorialData(), *Tp, Sp, Gp, DE, data.clusterTrack(),ctx);
} else if (!m_useCaloSeeds) {
++data.summaryStatAll()[kTotalBremSeeds][K];
tracks = m_tracksfinder->getTracksWithBrem(data.combinatorialData(), *Tp, Sp, Gp, DE, data.clusterTrack(), false,ctx);
- } else if (isCaloCompatible(data)) {
+ }
+ /// Note: The branch below is the one taken in ATLAS default inside-out tracking for run-3
+ else if (isCaloCompatible(data)) {
+
++data.summaryStatAll()[kTotalBremSeeds][K];
tracks = m_tracksfinder->getTracksWithBrem(data.combinatorialData(), *Tp, Sp, Gp, DE, data.clusterTrack(), true,ctx);
} else {
tracks = m_tracksfinder->getTracks (data.combinatorialData(), *Tp, Sp, Gp, DE, data.clusterTrack(),ctx);
}
- std::array<bool,SiCombinatorialTrackFinderData_xk::kNCombStats> inf{0,0,0,0,0,0}; m_tracksfinder->fillStatistic(data.combinatorialData(),inf);
+ /// update stat tables
+ std::array<bool,SiCombinatorialTrackFinderData_xk::kNCombStats> inf{0,0,0,0,0,0};
+ m_tracksfinder->fillStatistic(data.combinatorialData(),inf);
for (size_t p =0; p<inf.size(); ++p){
if(inf[p]) ++data.summaryStatAll()[m_indexToEnum[p]][K];
}
+ /// update the cluster-track-map to allow to filter any
+ /// upcoming seeds with hits that are already taken
if (m_seedsfilter) {
std::list<Trk::Track*>::iterator t = tracks.begin();
while (t!=tracks.end()) {
+ /// require sufficient free clusters on track
if (!isNewTrack(data, *t)) {
delete (*t);
tracks.erase(t++);
@@ -814,17 +838,23 @@ const Trk::TrackParameters* InDet::SiTrackMaker_xk::getAtaPlane
bool sss,
const std::vector<const Trk::SpacePoint*>& SP) const
{
+ /// we need at least three space points on the seed.
if (SP.size() < 3) return nullptr;
+ /// get the first cluster on the first hit
const Trk::PrepRawData* cl = SP[0]->clusterList().first;
if (!cl) return nullptr;
+ /// and use the surface from this cluster as our reference plane
const Trk::PlaneSurface* pla =
static_cast<const Trk::PlaneSurface*>(&cl->detectorElement()->surface());
if (!pla) return nullptr;
+ /// write the global positions into arrays. This includes an improved position estimate
+ /// for strip spacepoints. If this improvement fails, the method can return false -> then we abort
double p0[3],p1[3],p2[3];
if (!globalPositions(*(SP[0]),*(SP[1]),*(SP[2]),p0,p1,p2)) return nullptr;
+ /// translate second and third SP w.r.t first one
double x0 = p0[0] ;
double y0 = p0[1] ;
double z0 = p0[2] ;
@@ -834,62 +864,92 @@ const Trk::TrackParameters* InDet::SiTrackMaker_xk::getAtaPlane
double y2 = p2[1]-y0;
double z2 = p2[2]-z0;
+ /// distance of second SP to first in transverse plane
+ /// Also happens to be u-coordinate of second SP in conformal mapping
double u1 = 1./sqrt(x1*x1+y1*y1) ;
+ /// denominator for conformal mapping
double rn = x2*x2+y2*y2 ;
double r2 = 1./rn ;
+ /// coordinate system for conformal mapping - this is local x
double a = x1*u1 ;
double b = y1*u1 ;
+ /// u-coordinate of third SP in conformal mapping
double u2 = (a*x2+b*y2)*r2 ;
+ /// v-coordinate of third SP in conformal mapping
double v2 = (a*y2-b*x2)*r2 ;
- double A = v2/(u2-u1) ;
- double B = 2.*(v2-A*u2) ;
- double C = B/sqrt(1.+A*A) ;
- double T = z2*sqrt(r2)/(1.+.04*C*C*rn);
+ /// A,B are slope and intercept of the straight line in the u,v plane
+ /// connecting the three points.
+ double A = v2/(u2-u1) ; /// Keep in mind that v1 == 0
+ double B = 2.*(v2-A*u2) ; /// From inserting A into linear equation. Note that Igor sneaks in a factor two here
+ double C = B/sqrt(1.+A*A) ; /// Curvature estimate. (2R)²=(1+A²)/b² => 1/2R = b/sqrt(1+A²) = B / sqrt(1+A²).
+ double T = z2*sqrt(r2)/(1.+.04*C*C*rn); /// estimate of the track dz/dr (1/tanTheta), corrected for curvature effects
const Amg::Transform3D& Tp = pla->transform();
- double Ax[3] = {Tp(0,0),Tp(1,0),Tp(2,0)};
+ /// local x of the surface in the global frame
+ double Ax[3] = {Tp(0,0),Tp(1,0),Tp(2,0)};
+ /// local y of the surface in the global frame
double Ay[3] = {Tp(0,1),Tp(1,1),Tp(2,1)};
+ /// centre of the surface in the global frame
double D [3] = {Tp(0,3),Tp(1,3),Tp(2,3)};
-
+ /// location of the first SP w.r.t centre of the surface
double d[3] = {x0-D[0],y0-D[1],z0-D[2]};
-
+ /// local x and y - project onto local axes
data.par()[0] = d[0]*Ax[0]+d[1]*Ax[1]+d[2]*Ax[2];
data.par()[1] = d[0]*Ay[0]+d[1]*Ay[1]+d[2]*Ay[2];
+ /// silently switch off field if solenoid is off
Trk::MagneticFieldMode fieldModeEnum(m_fieldModeEnum);
if (!fieldCache.solenoidOn()) fieldModeEnum = Trk::NoField;
+
Trk::MagneticFieldProperties fieldprop(fieldModeEnum);
+ /// if we are not running with "no field":
if (fieldprop.magneticFieldMode() > 0) {
+ /// get the field at our first SP
double H[3],gP[3] ={x0,y0,z0};
fieldCache.getFieldZR(gP, H);
+ /// field is in kiloTesla. So here we check for more
+ /// than 0.1 Tesla
if (fabs(H[2])>.0001) {
+ /// phi estimate
data.par()[2] = atan2(b+a*A,a-b*A);
+ /// theta estimate
data.par()[3] = atan2(1.,T) ;
+ /// inverse transverse momentum estimate
data.par()[5] = -C/(300.*H[2]) ;
} else {
- T = z2*sqrt(r2) ;
- data.par()[2] = atan2(y2,x2);
+ /// if we have low field, use a
+ /// straight-line estimate
+ T = z2*sqrt(r2) ; /// note: Now no curvature correction
+ data.par()[2] = atan2(y2,x2);
data.par()[3] = atan2(1.,T) ;
- data.par()[5] = 1./m_pTmin ;
+ data.par()[5] = 1./m_pTmin ; /// no pt estimate, assume min pt
}
- } else {
+ }
+ /// treat absence of solenoid like the low-field case
+ else {
T = z2*sqrt(r2) ;
data.par()[2] = atan2(y2,x2);
data.par()[3] = atan2(1.,T) ;
data.par()[5] = 1./m_pTmin ;
}
+ /// apply the pt on the initial parameter estimate, with some margin
if (fabs(data.par()[5])*m_pTmin > 1.1) return nullptr;
- data.par()[4] = data.par()[5]/sqrt(1.+T*T);
- data.par()[6] = x0 ;
+ /// update qoverp
+ data.par()[4] = data.par()[5]/sqrt(1.+T*T); /// qoverp from qoverpt and theta
+ data.par()[6] = x0 ; /// ref point = first SP
data.par()[7] = y0 ;
data.par()[8] = z0 ;
+ /// never done in main ATLAS tracking. Would otherwise check if the seed is compatible with a
+ /// hadronic ROI
if (sss && !isHadCaloCompatible(data)) return nullptr;
+ /// now we can return the initial track parameters we built, parameterised using the ref surface.
+ /// Pass a nullptr for the covariance
return pla->createTrackParameters(data.par()[0],data.par()[1],data.par()[2],data.par()[3],data.par()[4],0);
}
@@ -1061,48 +1121,65 @@ void InDet::SiTrackMaker_xk::detectorElementsSelection(SiTrackMakerEventData_xk&
bool InDet::SiTrackMaker_xk::newSeed(SiTrackMakerEventData_xk& data, const std::vector<const Trk::SpacePoint*>& Sp) const
{
std::multiset<const Trk::Track*> trackseed;
- std::multimap<const Trk::PrepRawData*,const Trk::Track*>::const_iterator pt,pte = data.clusterTrack().end();
- std::vector<const Trk::SpacePoint*>::const_iterator s=Sp.begin(),se=Sp.end();
+ std::multimap<const Trk::PrepRawData*,const Trk::Track*>::const_iterator iter_clusterOnTrack,iter_clusterOnTrackEnd = data.clusterTrack().end();
+ /// counter for clusters on track
size_t n = 0;
- for (;s!=se; ++s) {
- const Trk::PrepRawData* p = (*s)->clusterList().first;
-
- for (pt = data.clusterTrack().find(p); pt!=pte; ++pt) {
- if ((*pt).first!=p) break;
- trackseed.insert((*pt).second);
+ for (const Trk::SpacePoint* spacePoint : Sp) {
+
+ /// start by checking the first cluster - always needed
+ const Trk::PrepRawData* prd = spacePoint->clusterList().first;
+
+ /// lookup if the cluster is already used by another track
+ for (iter_clusterOnTrack = data.clusterTrack().find(prd); iter_clusterOnTrack!=iter_clusterOnTrackEnd; ++iter_clusterOnTrack) {
+ if ((*iter_clusterOnTrack).first!=prd) break;
+ /// add tracks consuming our seed space-points to the trackseed list.
+ trackseed.insert((*iter_clusterOnTrack).second);
}
+ /// increment cluster counter
++n;
- p = (*s)->clusterList().second;
- if (!p) continue;
-
- for (pt = data.clusterTrack().find(p); pt!=pte; ++pt) {
- if ((*pt).first!=p) break;
- trackseed.insert((*pt).second);
+ /// now, prepare to check also the second cluster on any strip seed
+ prd = spacePoint->clusterList().second;
+ /// if we don't have one, nothing to do
+ if (!prd) continue;
+ /// otherwise, same game as before. Note that a track consuming both clusters on a strip hit is counted twice into the map
+ for (iter_clusterOnTrack = data.clusterTrack().find(prd); iter_clusterOnTrack!=iter_clusterOnTrackEnd; ++iter_clusterOnTrack) {
+ if ((*iter_clusterOnTrack).first!=prd) break;
+ trackseed.insert((*iter_clusterOnTrack).second);
}
+ /// incremenent counter again
++n;
-
}
- if(trackseed.size() < n) return true;
- if( m_heavyion && n==3 ) return true;
+ /// check if at least on cluster is not already used by any track.
+ /// This works since the multiset allows adding the same track multiple times
+ /// If this is the case, we accept the seed.
+ if(trackseed.size() < n) return true;
+ /// in the case of HI reco, we accept any 3-cluster (PPP) seed.
+ if( m_heavyion && n==3 ) return true;
- std::multiset<const Trk::Track*>::iterator t = trackseed.begin(), te = trackseed.end();
-
- const Trk::Track* tr = (*t) ;
-
- size_t nt = 1 ;
-
- for(++t; t!=te; ++t) {
-
- if((*t) != tr) {
- tr = (*t);
- nt = 1;
+ /// Now we look for the track consuming the largest number of clusters
+
+ /// This is done by looping over all tracks using any of our clusters,
+ /// and counting the appearance of each track in the multiset.
+ /// If one single track contains all of the clusters (--> is included n times),
+ /// we reject this seed.
+ const Trk::Track* currentTrack {nullptr};
+ size_t clustersOnCurrent = 1;
+ /// loop over the list of tracks
+ for(const Trk::Track* track : trackseed) {
+ /// if this is a new track, reset the counter
+ if(track != currentTrack) {
+ currentTrack = track;
+ clustersOnCurrent = 1;
continue ;
}
- if(++nt == n) return false;
+ /// otherwise increment the counter.
+ /// If this track has all clusters from the seed on it, reject the event
+ if(++clustersOnCurrent == n) return false;
}
- return nt!=n;
+ /// If we have no single track 'eating' all of our clusters, we accept the seed
+ return clustersOnCurrent!=n;
}
///////////////////////////////////////////////////////////////////
@@ -1206,6 +1283,7 @@ bool InDet::SiTrackMaker_xk::globalPositions
double* p0, double* p1, double* p2) const
{
+ /// first, fill the arrays with the global positions of the 3 points
p0[0] = s0.globalPosition().x();
p0[1] = s0.globalPosition().y();
p0[2] = s0.globalPosition().z();
@@ -1218,12 +1296,16 @@ bool InDet::SiTrackMaker_xk::globalPositions
p2[1] = s2.globalPosition().y();
p2[2] = s2.globalPosition().z();
+ /// for PPP seeds, we are done
if (!s0.clusterList().second && !s1.clusterList().second && !s2.clusterList().second) return true;
+ /// for SSS, we need some extra work
double dir0[3],dir1[3],dir2[3];
globalDirections(p0,p1,p2,dir0,dir1,dir2);
+ /// try to refine the position estimate for strip SP using the cluster information
+ /// in combination with the direction estimate
if (s0.clusterList().second && !globalPosition(s0,dir0,p0)) return false;
if (s1.clusterList().second && !globalPosition(s1,dir1,p1)) return false;
if (s2.clusterList().second && !globalPosition(s2,dir2,p2)) return false;
@@ -1235,15 +1317,20 @@ bool InDet::SiTrackMaker_xk::globalPositions
// Calculation global position for space points
///////////////////////////////////////////////////////////////////
+/// This is a refinement of the global position for strip space-points.
+/// It uses the direction estimate to fix the hit position along the
+/// strip axis (non-sensitive direction).
bool InDet::SiTrackMaker_xk::globalPosition
(const Trk::SpacePoint& sp, double* dir,double* p) const
{
+ /// pick up the two components of the space point
const Trk::PrepRawData* c0 = sp.clusterList().first;
const Trk::PrepRawData* c1 = sp.clusterList().second;
const InDetDD::SiDetectorElement* de0 = static_cast<const InDet::SiCluster*>(c0)->detectorElement();
const InDetDD::SiDetectorElement* de1 = static_cast<const InDet::SiCluster*>(c1)->detectorElement();
+ /// get the two ends of the strip in the global frame for each of the two clusters
Amg::Vector2D localPos = c0->localPosition();
std::pair<Amg::Vector3D,Amg::Vector3D> e0
(de0->endsOfStrip(InDetDD::SiLocalPosition(localPos.y(),localPos.x(),0.)));
@@ -1252,28 +1339,56 @@ bool InDet::SiTrackMaker_xk::globalPosition
std::pair<Amg::Vector3D,Amg::Vector3D> e1
(de1->endsOfStrip(InDetDD::SiLocalPosition(localPos.y(),localPos.x(),0.)));
+ /// get the "strip axis" in the global frame for each of the two clusters
+ /// these define two planes in which the strips are sensitive
double a0[3] = {e0.second.x()-e0.first.x(), e0.second.y()-e0.first.y(), e0.second.z()-e0.first.z()};
double a1[3] = {e1.second.x()-e1.first.x(), e1.second.y()-e1.first.y(), e1.second.z()-e1.first.z()};
+ /// get the connection vector between the bottom ends of the two strips
double dr[3] = {e1.first .x()-e0.first.x(), e1.first .y()-e0.first.y(), e1.first .z()-e0.first.z()};
+ /// divide max distance by the first strip length
double d0 = m_distmax/sqrt(a0[0]*a0[0]+a0[1]*a0[1]+a0[2]*a0[2]);
- // u = a1 x dir and v = a0 x dir
- //
+ /// Get the cross products of the direction vector and the strip axes.
+ /// This gives us a normal representation of the planes containing both.
double u[3] = {a1[1]*dir[2]-a1[2]*dir[1],a1[2]*dir[0]-a1[0]*dir[2],a1[0]*dir[1]-a1[1]*dir[0]};
double v[3] = {a0[1]*dir[2]-a0[2]*dir[1],a0[2]*dir[0]-a0[0]*dir[2],a0[0]*dir[1]-a0[1]*dir[0]};
+ /// The two planes are still only defined up to a shift along the normal.
+ /// Now, we fix this degree of freedom by requiring that the u-plane (perpendicular to momentum and second strip direction)
+ /// should actually contain the second strip itself.
+ /// To do so, we virtually 'hold' the u-plane at its intersection plane with the first strip,
+ /// starting with it touching the very bottom of the first strip, and then move the intersection
+ /// point along the strip. We can thus parameterise the shift by the distance we need to walk this way
+
+ /// Find the maximum distance we can go until we have walked past the full length of the first strip.
+ /// Equivalent to the projection of the full strip length on the u-plane normal.
double du = a0[0]*u[0]+a0[1]*u[1]+a0[2]*u[2];
+ /// if no such component, bail out - no solution exists
if (du==0. ) return false;
- //these need checking; original code calculated dv and used
- //s1 = -(dr[0]*v[0]+dr[1]*v[1]+dr[2]*v[2])/dv;
+ /// The distance we need to walk to contain the second strip in the u-plane is
+ /// is obtained by projecting the separation vector between the strip starting points
+ /// onto the normal.
+ /// (remember, we virtually start our shift with the u-plane intersecting the first strip
+ /// at its starting point).
+ /// We normalise the shifts to the maximum allowed shift.
double s0 = (dr[0]*u[0]+dr[1]*u[1]+dr[2]*u[2])/du;
+ /// this is playing the same game, but for shifting the v-plane w.r.t the
+ /// version that intersects with the bottom of the second strip until we contain
+ /// the first strip in it.
+ /// du has the same length as dv (symmetry of the problem)
double s1 = (dr[0]*v[0]+dr[1]*v[1]+dr[2]*v[2])/du;
+ /// check that the result is valid. We want the final fixed planes to still
+ /// intersect the respective other strip within some tolerance.
+ /// Keep in mind that s=0 --> intersect start of strip, s == 1 --> intersect end of strip.
+ /// We apply a tolerance beyond each end
if (s0 < -d0 || s0 > 1.+d0 || s1 < -d0 || s1 > 1.+d0) return false;
+ /// now, update the position estimate as the intersection of the
+ /// updated u-plane with the first strip.
p[0] = e0.first.x()+s0*a0[0];
p[1] = e0.first.y()+s0*a0[1];
p[2] = e0.first.z()+s0*a0[2];
@@ -1361,39 +1476,76 @@ bool InDet::SiTrackMaker_xk::isDBMSeeds(const Trk::SpacePoint* s) const
// Calculation global direction for positions of space points
///////////////////////////////////////////////////////////////////
+/// Here, we derive the local track direction in the space-points as the
+/// tangents to the estimated trajectory (assuming a circle in x-y and straight
+/// line in r-z)
void InDet::SiTrackMaker_xk::globalDirections
(double* p0, double* p1, double* p2, double* d0, double* d1, double* d2) const
{
+ /// transform transverse coordinates relative to the first SP
double x01 = p1[0]-p0[0] ;
double y01 = p1[1]-p0[1] ;
double x02 = p2[0]-p0[0] ;
double y02 = p2[1]-p0[1] ;
+ /// Now apply the same transform used in the seed maker:
+ /// x,y --> u:=x/(x²+y²); v:=y/(x²+y²);
+ /// in a frame centered around first SP
+
+ /// set x axis as direction unit vector in transverse plane, from first to second SP
double d01 = x01*x01+y01*y01 ;
double x1 = sqrt(d01) ;
- double u01 = 1./x1 ;
+ double u01 = 1./x1 ; /// reasoning: u = sqrt(d01)/d01, and v01 = 0
double a = x01*u01 ;
double b = y01*u01 ;
+ /// decompose the SP3-SP1 separation into components parallel...
double x2 = a*x02+b*y02 ;
+ /// and orthogonal to our new local x axis
double y2 = a*y02-b*x02 ;
+
+ /// squared distance third to first SP the transverse plane
double d02 = x2*x2+y2*y2 ;
+ /// u,v coordinates of third SP in the new frame
double u02 = x2/d02 ;
double v02 = y2/d02 ;
-
- double A0 = v02 /(u02-u01) ;
- double B0 = 2.*(v02-A0*u02) ;
+ /// Now obtain A,B parameters from circle parameterisation in the new frame:
+ /// V = (-x0/y0) x U + 1/(2y0) =: A x U + B.
+ /// A ~= delta V / delta U,
+ /// B = V - A x U, add a factor 2 for utility
+ double A0 = v02 /(u02-u01) ; /// v01 is 0, as the direction from SP 1 to 2 is the u-axis
+ double B0 = 2.*(v02-A0*u02) ;
+
+ /// Now we can resolve the related geometry.
+ /// Note that A0, in the original frame,
+ /// is related to the angle Phi1 between the tangent the
+ /// central SP and the tendon between the first and second
+ /// SP, alpha, as tan(alpha) = A0.
double C2 = (1.-B0*y2) ;
double S2 = (A0+B0*x2) ;
- double T = (p2[2]-p0[2])/sqrt(d02);
- double Se = 1./sqrt(1.+T*T) ;
- double Ce = Se*T ;
- Se /= sqrt(1.+A0*A0) ;
- double Sa = Se*a ;
- double Sb = Se*b ;
-
- d0[0] = Sa -Sb*A0; d0[1]= Sa*A0+Sb ; d0[2]=Ce;
- d1[0] = Sa +Sb*A0; d1[1]= Sb -Sa*A0; d1[2]=Ce;
- d2[0] = Sa*C2-Sb*S2; d2[1]= Sa*S2+Sb*C2; d2[2]=Ce;
+ double T = (p2[2]-p0[2])/sqrt(d02); /// dz/dr along the seed = 1 / tan(theta)
+ double sinTheta = 1./sqrt(1.+T*T) ; /// 1/sqrt(1+1/tan²(theta)) -> sin(theta)
+ double cosTheta = sinTheta*T ; /// cos theta (= sin(theta)/tan(theta))
+ double sinThetaCosAlpha = sinTheta / sqrt(1.+A0*A0) ; /// multiply with cos(alpha) via A0
+ double Sa = sinThetaCosAlpha*a ;
+ double Sb = sinThetaCosAlpha*b ;
+
+ /// (a,b) parameterises the direction corresponding to the tendon between the first two SP.
+ /// Now, go from there to the local tangents by removing or adding the angle 'alpha' from before, derived from A0.
+ /// direction at first SP - rotated by - 1 x alpha w.r.t the tendon defining a and b
+ /// The structure below comes from the formulae for sin / cos of a sum of two angles with (a,b) = r(cos phi / - sin phi)
+ d0[0] = Sa -Sb*A0; /// px0: a sin theta cos alpha - b sin theta sin alpha
+ d0[1]= Sb +Sa*A0; /// py0: b sin theta cos alpha + a sin theta sin alpha
+ d0[2]=cosTheta; /// pz0: cos theta
+
+ /// direction at second SP - rotated by + 1 x alpha w.r.t the tendon
+ d1[0] = Sa +Sb*A0; /// px1: a sin theta cos alpha + b sin theta sin alpha
+ d1[1]= Sb -Sa*A0; /// py1: b sin theta cos alpha - a sin theta sin alpha
+ d1[2]=cosTheta; /// pz1: cos theta
+
+ /// direction at third SP
+ d2[0] = Sa*C2-Sb*S2; /// px2: a sin theta cos alpha * C2 - b sin theta cos alpha S2
+ d2[1]= Sb*C2+Sa*S2; /// py2: b sin theta cos alpha * C2 + a sin theta cos alpha S2
+ d2[2]=cosTheta;
}