+/// Test the compatibility between the track and the 2 clusters together (using trackParam's covariance matrix):
+/// return the corresponding Chi2 accounting for covariances between the 2 clusters
+/// return trackParamAtCluster1 & 2
+
+ // extrapolate track parameters at the z position of the second cluster (no need to extrapolate the covariances)
+ trackParamAtCluster2.SetParameters(trackParamAtCluster1.GetParameters());
+ trackParamAtCluster2.SetZ(trackParamAtCluster1.GetZ());
+ AliMUONTrackExtrap::ExtrapToZ(&trackParamAtCluster2, cluster2->GetZ());
+
+ // set pointer to cluster2 into trackParamAtCluster2
+ trackParamAtCluster2.SetClusterPtr(cluster2);
+
+ // Set differences between track and the 2 clusters in the bending and non bending directions
+ AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr();
+ TMatrixD dPos(4,1);
+ dPos(0,0) = cluster1->GetX() - trackParamAtCluster1.GetNonBendingCoor();
+ dPos(1,0) = cluster1->GetY() - trackParamAtCluster1.GetBendingCoor();
+ dPos(2,0) = cluster2->GetX() - trackParamAtCluster2.GetNonBendingCoor();
+ dPos(3,0) = cluster2->GetY() - trackParamAtCluster2.GetBendingCoor();
+
+ // Calculate the error matrix from the track parameter covariances at first cluster
+ TMatrixD error(4,4);
+ error.Zero();
+ if (trackParamAtCluster1.CovariancesExist()) {
+ // Save track parameters at first cluster
+ TMatrixD paramAtCluster1Save(trackParamAtCluster1.GetParameters());
+
+ // Save track coordinates at second cluster
+ Double_t nonBendingCoor2 = trackParamAtCluster2.GetNonBendingCoor();
+ Double_t bendingCoor2 = trackParamAtCluster2.GetBendingCoor();
+
+ // copy track parameters at first cluster for jacobian calculation
+ AliMUONTrackParam trackParam(trackParamAtCluster1);
+
+ // add MCS effect to the covariance matrix at first cluster
+ AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
+
+ // Get the pointer to the parameter covariance matrix at first cluster
+ const TMatrixD& kParamCov = trackParam.GetCovariances();
+
+ // Calculate the jacobian related to the transformation between track parameters
+ // at first cluster and track coordinates at the 2 cluster z-positions
+ TMatrixD jacob(4,5);
+ jacob.Zero();
+ // first derivative at the first cluster:
+ jacob(0,0) = 1.; // dx1/dx
+ jacob(1,2) = 1.; // dy1/dy
+ // first derivative at the second cluster:
+ TMatrixD dParam(5,1);
+ Double_t direction[5] = {-1.,-1.,1.,1.,-1.};
+ for (Int_t i=0; i<5; i++) {
+ // Skip jacobian calculation for parameters with no associated error
+ if (kParamCov(i,i) == 0.) continue;
+ // Small variation of parameter i only
+ for (Int_t j=0; j<5; j++) {
+ if (j==i) {
+ dParam(j,0) = TMath::Sqrt(kParamCov(i,i));
+ dParam(j,0) *= TMath::Sign(1.,direction[j]*paramAtCluster1Save(j,0)); // variation always in the same direction
+ } else dParam(j,0) = 0.;
+ }
+
+ // Set new track parameters at first cluster
+ trackParam.SetParameters(paramAtCluster1Save);
+ trackParam.AddParameters(dParam);
+ trackParam.SetZ(cluster1->GetZ());
+
+ // Extrapolate new track parameters to the z position of the second cluster
+ AliMUONTrackExtrap::ExtrapToZ(&trackParam,cluster2->GetZ());
+
+ // Calculate the jacobian
+ jacob(2,i) = (trackParam.GetNonBendingCoor() - nonBendingCoor2) / dParam(i,0); // dx2/dParami
+ jacob(3,i) = (trackParam.GetBendingCoor() - bendingCoor2 ) / dParam(i,0); // dy2/dParami
+ }
+
+ // Calculate the error matrix
+ TMatrixD tmp(jacob,TMatrixD::kMult,kParamCov);
+ error = TMatrixD(tmp,TMatrixD::kMultTranspose,jacob);
+ }
+
+ // Add cluster resolution to the error matrix
+ error(0,0) += cluster1->GetErrX2();
+ error(1,1) += cluster1->GetErrY2();
+ error(2,2) += cluster2->GetErrX2();
+ error(3,3) += cluster2->GetErrY2();
+
+ // invert the error matrix for Chi2 calculation
+ if (error.Determinant() != 0) {
+ error.Invert();
+ } else {
+ AliWarning(" Determinant error=0");
+ return 1.e10;
+ }
+
+ // Compute the Chi2 value
+ TMatrixD tmp2(dPos,TMatrixD::kTransposeMult,error);
+ TMatrixD result(tmp2,TMatrixD::kMult,dPos);
+
+ return result(0,0);
+
+}
+
+ //__________________________________________________________________________
+void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster)
+{
+ /// Add 1 cluster to the track candidate
+ /// Update chi2 of the track
+
+ // Compute local chi2
+ AliMUONVCluster* cluster = trackParamAtCluster.GetClusterPtr();
+ Double_t deltaX = trackParamAtCluster.GetNonBendingCoor() - cluster->GetX();
+ Double_t deltaY = trackParamAtCluster.GetBendingCoor() - cluster->GetY();
+ Double_t localChi2 = deltaX*deltaX / cluster->GetErrX2() +
+ deltaY*deltaY / cluster->GetErrY2();
+
+ // Flag cluster as being not removable
+ if (GetRecoParam()->RequestStation(cluster->GetChamberId()/2))
+ trackParamAtCluster.SetRemovable(kFALSE);
+ else trackParamAtCluster.SetRemovable(kTRUE);
+ trackParamAtCluster.SetLocalChi2(0.); // --> Local chi2 not used
+
+ // Update the chi2 of the new track
+ track.SetGlobalChi2(track.GetGlobalChi2() + localChi2);
+
+ // Update TrackParamAtCluster
+ track.AddTrackParamAtCluster(trackParamAtCluster,*cluster);
+
+}
+
+ //__________________________________________________________________________
+void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster1, AliMUONTrackParam &trackParamAtCluster2)
+{
+ /// Add 2 clusters to the track candidate
+ /// Update track and local chi2
+
+ // Update local chi2 at first cluster
+ AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr();
+ Double_t deltaX = trackParamAtCluster1.GetNonBendingCoor() - cluster1->GetX();
+ Double_t deltaY = trackParamAtCluster1.GetBendingCoor() - cluster1->GetY();
+ Double_t localChi2AtCluster1 = deltaX*deltaX / cluster1->GetErrX2() +
+ deltaY*deltaY / cluster1->GetErrY2();
+ trackParamAtCluster1.SetLocalChi2(localChi2AtCluster1);
+
+ // Flag first cluster as being removable
+ trackParamAtCluster1.SetRemovable(kTRUE);
+
+ // Update local chi2 at second cluster
+ AliMUONVCluster* cluster2 = trackParamAtCluster2.GetClusterPtr();
+ deltaX = trackParamAtCluster2.GetNonBendingCoor() - cluster2->GetX();
+ deltaY = trackParamAtCluster2.GetBendingCoor() - cluster2->GetY();
+ Double_t localChi2AtCluster2 = deltaX*deltaX / cluster2->GetErrX2() +
+ deltaY*deltaY / cluster2->GetErrY2();
+ trackParamAtCluster2.SetLocalChi2(localChi2AtCluster2);
+
+ // Flag first cluster as being removable
+ trackParamAtCluster2.SetRemovable(kTRUE);
+
+ // Update the chi2 of the new track
+ track.SetGlobalChi2(track.GetGlobalChi2() + localChi2AtCluster1 + localChi2AtCluster2);
+
+ // Update TrackParamAtCluster
+ track.AddTrackParamAtCluster(trackParamAtCluster1,*cluster1);
+ track.AddTrackParamAtCluster(trackParamAtCluster2,*cluster2);
+
+}
+
+ //__________________________________________________________________________
+Bool_t AliMUONTrackReconstructor::RecoverTrack(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation)
+{
+ /// Try to recover the track candidate in the next station
+ /// by removing the worst of the two clusters attached in the current station
+ /// Return kTRUE if recovering succeeds
+ AliDebug(1,"Enter RecoverTrack");
+
+ // Do not try to recover track until we have attached cluster(s) on station(1..) 3
+ if (nextStation > 1) return kFALSE;
+
+ Int_t worstClusterNumber = -1;
+ Double_t localChi2, worstLocalChi2 = -1.;
+
+ // Look for the cluster to remove
+ for (Int_t clusterNumber = 0; clusterNumber < 2; clusterNumber++) {
+ AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(clusterNumber);
+
+ // check if current cluster is in the previous station
+ if (trackParamAtCluster->GetClusterPtr()->GetChamberId()/2 != nextStation+1) break;
+
+ // check if current cluster is removable
+ if (!trackParamAtCluster->IsRemovable()) return kFALSE;
+
+ // reset the current cluster as beig not removable if it is on a required station
+ if (GetRecoParam()->RequestStation(nextStation+1)) trackParamAtCluster->SetRemovable(kFALSE);
+
+ // Pick up cluster with the worst chi2
+ localChi2 = trackParamAtCluster->GetLocalChi2();
+ if (localChi2 > worstLocalChi2) {
+ worstLocalChi2 = localChi2;
+ worstClusterNumber = clusterNumber;
+ }
+
+ }
+
+ // check if worst cluster found
+ if (worstClusterNumber < 0) return kFALSE;
+
+ // Remove the worst cluster
+ trackCandidate.RemoveTrackParamAtCluster((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(worstClusterNumber));
+
+ // Re-fit the track:
+ // Do not take into account the multiple scattering to speed up the fit
+ // Calculate the track parameter covariance matrix
+ Fit(trackCandidate, kFALSE, kFALSE, kTRUE);
+
+ // skip track out of limits
+ if (!IsAcceptable(*((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First()))) return kFALSE;
+
+ // Look for new cluster(s) in next station
+ return FollowTrackInStation(trackCandidate,clusterStore,nextStation);
+
+}
+
+ //__________________________________________________________________________
+void AliMUONTrackReconstructor::SetVertexErrXY2ForFit(AliMUONTrack &trackCandidate)
+{
+ /// Compute the vertex resolution square from natural vertex dispersion and