]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONVTrackReconstructor.cxx
Fix histo name
[u/mrichter/AliRoot.git] / MUON / AliMUONVTrackReconstructor.cxx
index 902bb863c0ccc7d07856dffb7485f93d21a11164..fc4247b7220c54eb1eb5103878a9145d2b3e7721 100644 (file)
@@ -100,7 +100,8 @@ AliMUONVTrackReconstructor::AliMUONVTrackReconstructor(const AliMUONRecoParam* r
 fRecTracksPtr(0x0),
 fNRecTracks(0),
 fClusterServer(clusterServer),
-fkRecoParam(recoParam)
+fkRecoParam(recoParam),
+fMaxMCSAngle2(0.)
 {
   /// Constructor for class AliMUONVTrackReconstructor
   /// WARNING: if clusterServer=0x0, no clusterization will be possible at this level
@@ -110,6 +111,13 @@ fkRecoParam(recoParam)
   
   // set the magnetic field for track extrapolations
   AliMUONTrackExtrap::SetField();
+  
+  // set the maximum MCS angle in chamber from the minimum acceptable momentum
+  AliMUONTrackParam param;
+  Double_t inverseBendingP = (GetRecoParam()->GetMinBendingMomentum() > 0.) ? 1./GetRecoParam()->GetMinBendingMomentum() : 1.;
+  param.SetInverseBendingMomentum(inverseBendingP);
+  fMaxMCSAngle2 = AliMUONTrackExtrap::GetMCSAngle2(param, AliMUONConstants::ChamberThicknessInX0(), 1.);
+  
 }
 
   //__________________________________________________________________________
@@ -138,20 +146,24 @@ void AliMUONVTrackReconstructor::EventReconstruct(AliMUONVClusterStore& clusterS
   // Reset array of tracks
   ResetTracks();
   
-  // Look for candidates from clusters in stations(1..) 4 and 5
-  MakeTrackCandidates(clusterStore);
+  // Look for candidates from clusters in stations(1..) 4 and 5 (abort in case of failure)
+  if (!MakeTrackCandidates(clusterStore)) return;
   
-  // Look for extra candidates from clusters in stations(1..) 4 and 5
-  if (GetRecoParam()->MakeMoreTrackCandidates()) MakeMoreTrackCandidates(clusterStore);
+  // Look for extra candidates from clusters in stations(1..) 4 and 5 (abort in case of failure)
+  if (GetRecoParam()->MakeMoreTrackCandidates()) {
+    if (!MakeMoreTrackCandidates(clusterStore)) return;
+  }
   
   // Stop tracking if no candidate found
   if (fRecTracksPtr->GetEntriesFast() == 0) return;
   
-  // Follow tracks in stations(1..) 3, 2 and 1
-  FollowTracks(clusterStore);
+  // Follow tracks in stations(1..) 3, 2 and 1 (abort in case of failure)
+  if (!FollowTracks(clusterStore)) return;
   
   // Complement the reconstructed tracks
-  if (GetRecoParam()->ComplementTracks()) ComplementTracks(clusterStore);
+  if (GetRecoParam()->ComplementTracks()) {
+    if (ComplementTracks(clusterStore)) RemoveIdenticalTracks();
+  }
   
   // Improve the reconstructed tracks
   if (GetRecoParam()->ImproveTracks()) ImproveTracks();
@@ -167,22 +179,100 @@ void AliMUONVTrackReconstructor::EventReconstruct(AliMUONVClusterStore& clusterS
   for (Int_t i=0; i<fNRecTracks; ++i) 
   {
     AliMUONTrack * track = (AliMUONTrack*) fRecTracksPtr->At(i);
-    track->SetUniqueID(i+1);
-    trackStore.Add(*track);
+    if (track->GetGlobalChi2() < AliMUONTrack::MaxChi2()) {
+      track->SetUniqueID(i+1);
+      trackStore.Add(*track);
+    } else AliWarning("problem occur somewhere during track refitting --> discard track");
   }
 }
 
-  //__________________________________________________________________________
+//__________________________________________________________________________
+Bool_t AliMUONVTrackReconstructor::IsAcceptable(AliMUONTrackParam &trackParam)
+{
+  /// Return kTRUE if the track is within given limits on momentum/angle/origin
+  
+  const TMatrixD& kParamCov = trackParam.GetCovariances();
+  Int_t chamber = trackParam.GetClusterPtr()->GetChamberId();
+  Double_t z = trackParam.GetZ();
+  Double_t sigmaCut = GetRecoParam()->GetSigmaCutForTracking();
+  
+  // MCS dipersion
+  Double_t angleMCS2 = 0.;
+  if (AliMUONTrackExtrap::IsFieldON() && chamber < 6)
+    angleMCS2 = AliMUONTrackExtrap::GetMCSAngle2(trackParam, AliMUONConstants::ChamberThicknessInX0(), 1.);
+  else angleMCS2 = fMaxMCSAngle2;
+  Double_t impactMCS2 = 0.;
+  if (!GetRecoParam()->SelectOnTrackSlope()) for (Int_t iCh=0; iCh<=chamber; iCh++)
+    impactMCS2 += AliMUONConstants::DefaultChamberZ(chamber) * AliMUONConstants::DefaultChamberZ(chamber) * angleMCS2;
+  
+  // ------ track selection in non bending direction ------
+  if (GetRecoParam()->SelectOnTrackSlope()) {
+    
+    // check if non bending slope is within tolerances
+    Double_t nonBendingSlopeErr = TMath::Sqrt(kParamCov(1,1) + (chamber + 1.) * angleMCS2);
+    if ((TMath::Abs(trackParam.GetNonBendingSlope()) - sigmaCut * nonBendingSlopeErr) > GetRecoParam()->GetMaxNonBendingSlope()) return kFALSE;
+    
+  } else {
+    
+    // or check if non bending impact parameter is within tolerances
+    Double_t nonBendingImpactParam = TMath::Abs(trackParam.GetNonBendingCoor() - z * trackParam.GetNonBendingSlope());
+    Double_t nonBendingImpactParamErr = TMath::Sqrt(kParamCov(0,0) + z * z * kParamCov(1,1) - 2. * z * kParamCov(0,1) + impactMCS2);
+    if ((nonBendingImpactParam - sigmaCut * nonBendingImpactParamErr) > (3. * GetRecoParam()->GetNonBendingVertexDispersion())) return kFALSE;
+    
+  }
+  
+  // ------ track selection in bending direction ------
+  if (AliMUONTrackExtrap::IsFieldON()) { // depending whether the field is ON or OFF
+    
+    // check if bending momentum is within tolerances
+    Double_t bendingMomentum = TMath::Abs(1. / trackParam.GetInverseBendingMomentum());
+    Double_t bendingMomentumErr = TMath::Sqrt(kParamCov(4,4)) * bendingMomentum * bendingMomentum;
+    if (chamber < 6 && (bendingMomentum + sigmaCut * bendingMomentumErr) < GetRecoParam()->GetMinBendingMomentum()) return kFALSE;
+    else if ((bendingMomentum + 3. * bendingMomentumErr) < GetRecoParam()->GetMinBendingMomentum()) return kFALSE;
+    
+  } else {
+    
+    if (GetRecoParam()->SelectOnTrackSlope()) {
+      
+      // check if bending slope is within tolerances
+      Double_t bendingSlopeErr = TMath::Sqrt(kParamCov(3,3) + (chamber + 1.) * angleMCS2);
+      if ((TMath::Abs(trackParam.GetBendingSlope()) - sigmaCut * bendingSlopeErr) > GetRecoParam()->GetMaxBendingSlope()) return kFALSE;
+      
+    } else {
+      
+      // or check if bending impact parameter is within tolerances
+      Double_t bendingImpactParam = TMath::Abs(trackParam.GetBendingCoor() - z * trackParam.GetBendingSlope());
+      Double_t bendingImpactParamErr = TMath::Sqrt(kParamCov(2,2) + z * z * kParamCov(3,3) - 2. * z * kParamCov(2,3) + impactMCS2);
+      if ((bendingImpactParam - sigmaCut * bendingImpactParamErr) > (3. * GetRecoParam()->GetBendingVertexDispersion())) return kFALSE;
+      
+    }
+    
+  }
+  
+  return kTRUE;
+  
+}
+
+//__________________________________________________________________________
 TClonesArray* AliMUONVTrackReconstructor::MakeSegmentsBetweenChambers(const AliMUONVClusterStore& clusterStore, Int_t ch1, Int_t ch2)
 {
   /// To make the list of segments from the list of clusters in the 2 given chambers.
   /// Return a new TClonesArray of segments.
   /// It is the responsibility of the user to delete it afterward.
   AliDebug(1,Form("Enter MakeSegmentsBetweenChambers (1..) %d-%d", ch1+1, ch2+1));
+  AliCodeTimerAuto("");
   
   AliMUONVCluster *cluster1, *cluster2;
   AliMUONObjectPair *segment;
-  Double_t nonBendingSlope = 0, bendingSlope = 0, impactParam = 0., bendingMomentum = 0.; // to avoid compilation warning
+  Double_t z1 = 0., z2 = 0., dZ = 0.;
+  Double_t nonBendingSlope = 0., nonBendingSlopeErr = 0., nonBendingImpactParam = 0., nonBendingImpactParamErr = 0.;
+  Double_t bendingSlope = 0., bendingSlopeErr = 0., bendingImpactParam = 0., bendingImpactParamErr = 0., bendingImpactParamErr2 = 0.;
+  Double_t bendingMomentum = 0., bendingMomentumErr = 0.;
+  Double_t bendingVertexDispersion2 = GetRecoParam()->GetBendingVertexDispersion() * GetRecoParam()->GetBendingVertexDispersion();
+  Double_t impactMCS2 = 0; // maximum impact parameter dispersion**2 due to MCS in chamber
+  if (!GetRecoParam()->SelectOnTrackSlope() || AliMUONTrackExtrap::IsFieldON()) for (Int_t iCh=0; iCh<=ch1; iCh++)
+    impactMCS2 += AliMUONConstants::DefaultChamberZ(iCh) * AliMUONConstants::DefaultChamberZ(iCh) * fMaxMCSAngle2;
+  Double_t sigmaCut = GetRecoParam()->GetSigmaCutForTracking();
   
   // Create iterators to loop over clusters in both chambers
   TIter nextInCh1(clusterStore.CreateChamberIterator(ch1,ch1));
@@ -193,41 +283,64 @@ TClonesArray* AliMUONVTrackReconstructor::MakeSegmentsBetweenChambers(const AliM
   
   // Loop over clusters in the first chamber of the station
   while ( ( cluster1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
+    z1 = cluster1->GetZ();
     
     // reset cluster iterator of chamber 2
     nextInCh2.Reset();
     
     // Loop over clusters in the second chamber of the station
     while ( ( cluster2 = static_cast<AliMUONVCluster*>(nextInCh2()) ) ) {
+      z2 = cluster2->GetZ();
+      dZ = z1 - z2;
       
-      // non bending slope
-      nonBendingSlope = (cluster1->GetX() - cluster2->GetX()) / (cluster1->GetZ() - cluster2->GetZ());
-      
-      // check if non bending slope is within tolerances
-      if (TMath::Abs(nonBendingSlope) > GetRecoParam()->GetMaxNonBendingSlope()) continue;
-      
-      // bending slope
-      bendingSlope = (cluster1->GetY() - cluster2->GetY()) / (cluster1->GetZ() - cluster2->GetZ());
-      
-      // check the bending momentum of the bending slope depending if the field is ON or OFF
-      if (AliMUONTrackExtrap::IsFieldON()) {
+      // ------ track selection in non bending direction ------
+      nonBendingSlope = (cluster1->GetX() - cluster2->GetX()) / dZ;
+      if (GetRecoParam()->SelectOnTrackSlope()) {
        
-       // impact parameter
-       impactParam = cluster1->GetY() - cluster1->GetZ() * bendingSlope;
+       // check if non bending slope is within tolerances
+       nonBendingSlopeErr = TMath::Sqrt((cluster1->GetErrX2() + cluster2->GetErrX2()) / dZ / dZ + (ch1 + 1.) * fMaxMCSAngle2);
+       if ((TMath::Abs(nonBendingSlope) - sigmaCut * nonBendingSlopeErr) > GetRecoParam()->GetMaxNonBendingSlope()) continue;
        
-       // absolute value of bending momentum
-       bendingMomentum = TMath::Abs(AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(impactParam));
+      } else {
+       
+       // or check if non bending impact parameter is within tolerances
+       nonBendingImpactParam = TMath::Abs(cluster1->GetX() - cluster1->GetZ() * nonBendingSlope);
+       nonBendingImpactParamErr = TMath::Sqrt((z1 * z1 * cluster2->GetErrX2() + z2 * z2 * cluster1->GetErrX2()) / dZ / dZ + impactMCS2);
+       if ((nonBendingImpactParam - sigmaCut * nonBendingImpactParamErr) > (3. * GetRecoParam()->GetNonBendingVertexDispersion())) continue;
+       
+      }
+      
+      // ------ track selection in bending direction ------
+      bendingSlope = (cluster1->GetY() - cluster2->GetY()) / dZ;
+      if (AliMUONTrackExtrap::IsFieldON()) { // depending whether the field is ON or OFF
        
        // check if bending momentum is within tolerances
-       if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() ||
-           bendingMomentum > GetRecoParam()->GetMaxBendingMomentum()) continue;
+       bendingImpactParam = cluster1->GetY() - cluster1->GetZ() * bendingSlope;
+       bendingImpactParamErr2 = (z1 * z1 * cluster2->GetErrY2() + z2 * z2 * cluster1->GetErrY2()) / dZ / dZ + impactMCS2;
+       bendingMomentum = TMath::Abs(AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpactParam));
+       bendingMomentumErr = TMath::Sqrt((bendingVertexDispersion2 + bendingImpactParamErr2) /
+                                        bendingImpactParam / bendingImpactParam + 0.01) * bendingMomentum;
+       if ((bendingMomentum + 3. * bendingMomentumErr) < GetRecoParam()->GetMinBendingMomentum()) continue;
        
       } else {
        
-       // check if non bending slope is within tolerances
-       if (TMath::Abs(bendingSlope) > GetRecoParam()->GetMaxBendingSlope()) continue;
-      
+       if (GetRecoParam()->SelectOnTrackSlope()) {
+         
+         // check if bending slope is within tolerances
+         bendingSlopeErr = TMath::Sqrt((cluster1->GetErrY2() + cluster2->GetErrY2()) / dZ / dZ + (ch1 + 1.) * fMaxMCSAngle2);
+         if ((TMath::Abs(bendingSlope) - sigmaCut * bendingSlopeErr) > GetRecoParam()->GetMaxBendingSlope()) continue;
+         
+       } else {
+         
+         // or check if bending impact parameter is within tolerances
+         bendingImpactParam = TMath::Abs(cluster1->GetY() - cluster1->GetZ() * bendingSlope);
+         bendingImpactParamErr = TMath::Sqrt((z1 * z1 * cluster2->GetErrY2() + z2 * z2 * cluster1->GetErrY2()) / dZ / dZ + impactMCS2);
+         if ((bendingImpactParam - sigmaCut * bendingImpactParamErr) > (3. * GetRecoParam()->GetBendingVertexDispersion())) continue;
+         
+       }
+       
       }
+      
       // make new segment
       segment = new ((*segments)[segments->GetLast()+1]) AliMUONObjectPair(cluster1, cluster2, kFALSE, kFALSE);
       
@@ -251,23 +364,86 @@ TClonesArray* AliMUONVTrackReconstructor::MakeSegmentsBetweenChambers(const AliM
   return segments;
 }
 
+  //__________________________________________________________________________
+void AliMUONVTrackReconstructor::RemoveUsedSegments(TClonesArray& segments)
+{
+  /// To remove pairs of clusters already attached to a track
+  AliDebug(1,"Enter RemoveUsedSegments");
+  Int_t nSegments = segments.GetEntriesFast();
+  Int_t nTracks = fRecTracksPtr->GetEntriesFast();
+  AliMUONObjectPair *segment;
+  AliMUONTrack *track;
+  AliMUONVCluster *cluster, *cluster1, *cluster2;
+  Bool_t foundCluster1, foundCluster2, removeSegment;
+  
+  // Loop over segments
+  for (Int_t iSegment=0; iSegment<nSegments; iSegment++) {
+    segment = (AliMUONObjectPair*) segments.UncheckedAt(iSegment);
+    
+    cluster1 = (AliMUONVCluster*) segment->First();
+    cluster2 = (AliMUONVCluster*) segment->Second();
+    removeSegment = kFALSE;
+    
+    // Loop over tracks
+    for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
+      track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iTrack);
+      
+      // skip empty slot
+      if (!track) continue;
+      
+      foundCluster1 = kFALSE;
+      foundCluster2 = kFALSE;
+      
+      // Loop over clusters
+      Int_t nClusters = track->GetNClusters();
+      for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
+        cluster = ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iCluster))->GetClusterPtr();
+       
+       // check if both clusters are in that track
+       if (cluster == cluster1) foundCluster1 = kTRUE;
+       else if (cluster == cluster2) foundCluster2 = kTRUE;
+       
+       if (foundCluster1 && foundCluster2) {
+         removeSegment = kTRUE;
+         break;
+       }
+       
+      }
+      
+      if (removeSegment) break;
+      
+    }
+    
+    if (removeSegment) segments.RemoveAt(iSegment);
+      
+  }
+  
+  segments.Compress();
+  
+  // Printout for debug
+  AliDebug(1,Form("NSegments =  %d ", segments.GetEntriesFast()));
+}
+
   //__________________________________________________________________________
 void AliMUONVTrackReconstructor::RemoveIdenticalTracks()
 {
   /// To remove identical tracks:
   /// Tracks are considered identical if they have all their clusters in common.
   /// One keeps the track with the larger number of clusters if need be
-  AliMUONTrack *track1, *track2, *trackToRemove;
+  AliMUONTrack *track1, *track2;
+  Int_t nTracks = fRecTracksPtr->GetEntriesFast();
   Int_t clustersInCommon, nClusters1, nClusters2;
-  Bool_t removedTrack1;
   // Loop over first track of the pair
-  track1 = (AliMUONTrack*) fRecTracksPtr->First();
-  while (track1) {
-    removedTrack1 = kFALSE;
+  for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
+    track1 = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iTrack1);
+    // skip empty slot
+    if (!track1) continue;
     nClusters1 = track1->GetNClusters();
     // Loop over second track of the pair
-    track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
-    while (track2) {
+    for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
+      track2 = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iTrack2);
+      // skip empty slot
+      if (!track2) continue;
       nClusters2 = track2->GetNClusters();
       // number of clusters in common between two tracks
       clustersInCommon = track1->ClustersInCommon(track2);
@@ -276,27 +452,18 @@ void AliMUONVTrackReconstructor::RemoveIdenticalTracks()
         // decide which track to remove
         if (nClusters2 > nClusters1) {
          // remove track1 and continue the first loop with the track next to track1
-         trackToRemove = track1;
-         track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
-          fRecTracksPtr->Remove(trackToRemove);
-         fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
+          fRecTracksPtr->RemoveAt(iTrack1);
          fNRecTracks--;
-         removedTrack1 = kTRUE;
          break;
        } else {
          // remove track2 and continue the second loop with the track next to track2
-         trackToRemove = track2;
-         track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
-         fRecTracksPtr->Remove(trackToRemove);
-         fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
+         fRecTracksPtr->RemoveAt(iTrack2);
          fNRecTracks--;
         }
-      } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
+      }
     } // track2
-    if (removedTrack1) continue;
-    track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
   } // track1
-  return;
+  fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
 }
 
   //__________________________________________________________________________
@@ -307,46 +474,40 @@ void AliMUONVTrackReconstructor::RemoveDoubleTracks()
   /// which has the smaller number of clusters are in common with the other track.
   /// Among two identical tracks, one keeps the track with the larger number of clusters
   /// or, if these numbers are equal, the track with the minimum chi2.
-  AliMUONTrack *track1, *track2, *trackToRemove;
-  Int_t clustersInCommon, nClusters1, nClusters2;
-  Bool_t removedTrack1;
+  AliMUONTrack *track1, *track2;
+  Int_t nTracks = fRecTracksPtr->GetEntriesFast();
+  Int_t clustersInCommon2, nClusters1, nClusters2;
   // Loop over first track of the pair
-  track1 = (AliMUONTrack*) fRecTracksPtr->First();
-  while (track1) {
-    removedTrack1 = kFALSE;
+  for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
+    track1 = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iTrack1);
+    // skip empty slot
+    if (!track1) continue;
     nClusters1 = track1->GetNClusters();
     // Loop over second track of the pair
-    track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
-    while (track2) {
+    for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
+      track2 = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iTrack2);
+      // skip empty slot
+      if (!track2) continue;
       nClusters2 = track2->GetNClusters();
       // number of clusters in common between two tracks
-      clustersInCommon = track1->ClustersInCommon(track2);
+      clustersInCommon2 = 2 * track1->ClustersInCommon(track2);
       // check for identical tracks
-      if (((nClusters1 < nClusters2) && (2 * clustersInCommon > nClusters1)) || (2 * clustersInCommon > nClusters2)) {
+      if (clustersInCommon2 > nClusters1 || clustersInCommon2 > nClusters2) {
         // decide which track to remove
         if ((nClusters1 > nClusters2) || ((nClusters1 == nClusters2) && (track1->GetGlobalChi2() <= track2->GetGlobalChi2()))) {
          // remove track2 and continue the second loop with the track next to track2
-         trackToRemove = track2;
-         track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
-         fRecTracksPtr->Remove(trackToRemove);
-         fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
+         fRecTracksPtr->RemoveAt(iTrack2);
          fNRecTracks--;
         } else {
          // else remove track1 and continue the first loop with the track next to track1
-         trackToRemove = track1;
-         track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
-          fRecTracksPtr->Remove(trackToRemove);
-         fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
+          fRecTracksPtr->RemoveAt(iTrack1);
          fNRecTracks--;
-         removedTrack1 = kTRUE;
          break;
         }
-      } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
+      }
     } // track2
-    if (removedTrack1) continue;
-    track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
   } // track1
-  return;
+  fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
 }
 
   //__________________________________________________________________________
@@ -370,7 +531,8 @@ void AliMUONVTrackReconstructor::RemoveConnectedTracks(Bool_t inSt345)
     while (track2) {
       nClusters2 = track2->GetNClusters();
       // number of clusters in common between two tracks
-      clustersInCommon = track1->ClustersInCommon(track2, inSt345);
+      if (inSt345) clustersInCommon = track1->ClustersInCommonInSt345(track2);
+      else clustersInCommon = track1->ClustersInCommon(track2);
       // check for identical tracks
       if (clustersInCommon > 0) {
         // decide which track to remove
@@ -414,18 +576,20 @@ void AliMUONVTrackReconstructor::AskForNewClustersInChamber(const AliMUONTrackPa
   
   // extrapolate track parameters to the chamber
   AliMUONTrackParam extrapTrackParam(trackParam);
-  AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(chamber));
+  if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(chamber))) return;
   
-  // build the searching area using the track resolution and the maximum-distance-to-track value
+  // build the searching area using the track and chamber resolutions and the maximum-distance-to-track value
   const TMatrixD& kParamCov = extrapTrackParam.GetCovariances();
-  Double_t errX2 = kParamCov(0,0) + kMaxDZ * kMaxDZ * kParamCov(1,1) + 2. * kMaxDZ * TMath::Abs(kParamCov(0,1));
-  Double_t errY2 = kParamCov(2,2) + kMaxDZ * kMaxDZ * kParamCov(3,3) + 2. * kMaxDZ * TMath::Abs(kParamCov(2,3));
+  Double_t errX2 = kParamCov(0,0) + kMaxDZ * kMaxDZ * kParamCov(1,1) + 2. * kMaxDZ * TMath::Abs(kParamCov(0,1)) +
+                   GetRecoParam()->GetDefaultNonBendingReso(chamber) * GetRecoParam()->GetDefaultNonBendingReso(chamber);
+  Double_t errY2 = kParamCov(2,2) + kMaxDZ * kMaxDZ * kParamCov(3,3) + 2. * kMaxDZ * TMath::Abs(kParamCov(2,3)) +
+                  GetRecoParam()->GetDefaultBendingReso(chamber) * GetRecoParam()->GetDefaultBendingReso(chamber);
   Double_t dX = TMath::Abs(trackParam.GetNonBendingSlope()) * kMaxDZ +
                GetRecoParam()->GetMaxNonBendingDistanceToTrack() +
-               GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errX2);
+               GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errX2);
   Double_t dY = TMath::Abs(trackParam.GetBendingSlope()) * kMaxDZ +
                GetRecoParam()->GetMaxBendingDistanceToTrack() +
-               GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errY2);
+               GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errY2);
   AliMpArea area(extrapTrackParam.GetNonBendingCoor(), 
                  extrapTrackParam.GetBendingCoor(),
                  dX, dY);
@@ -454,11 +618,11 @@ Double_t AliMUONVTrackReconstructor::TryOneCluster(const AliMUONTrackParam &trac
 /// return trackParamAtCluster
   
   // extrapolate track parameters and covariances at the z position of the tested cluster
+  // and set pointer to cluster into trackParamAtCluster
   trackParamAtCluster = trackParam;
-  AliMUONTrackExtrap::ExtrapToZCov(&trackParamAtCluster, cluster->GetZ(), updatePropagator);
-  
-  // set pointer to cluster into trackParamAtCluster
   trackParamAtCluster.SetClusterPtr(cluster);
+  if (!AliMUONTrackExtrap::ExtrapToZCov(&trackParamAtCluster, cluster->GetZ(), updatePropagator))
+    return 2.*AliMUONTrack::MaxChi2();
   
   // Set differences between trackParam and cluster in the bending and non bending directions
   Double_t dX = cluster->GetX() - trackParamAtCluster.GetNonBendingCoor();
@@ -468,9 +632,12 @@ Double_t AliMUONVTrackReconstructor::TryOneCluster(const AliMUONTrackParam &trac
   const TMatrixD& kParamCov = trackParamAtCluster.GetCovariances();
   Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2();
   Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2();
+  Double_t covXY   = kParamCov(0,2);
+  Double_t det     = sigmaX2 * sigmaY2 - covXY * covXY;
   
   // Compute chi2
-  return dX * dX / sigmaX2 + dY * dY / sigmaY2;
+  if (det == 0.) return 2.*AliMUONTrack::MaxChi2();
+  return (dX * dX * sigmaY2 + dY * dY * sigmaX2 - 2. * dX * dY * covXY) / det;
   
 }
 
@@ -478,7 +645,7 @@ Double_t AliMUONVTrackReconstructor::TryOneCluster(const AliMUONTrackParam &trac
 Bool_t AliMUONVTrackReconstructor::TryOneClusterFast(const AliMUONTrackParam &trackParam, const AliMUONVCluster* cluster)
 {
 /// Test the compatibility between the track and the cluster
-/// given the track resolution + the maximum-distance-to-track value
+/// given the track and cluster resolutions + the maximum-distance-to-track value
 /// and assuming linear propagation of the track:
 /// return kTRUE if they are compatibles
   
@@ -486,12 +653,12 @@ Bool_t AliMUONVTrackReconstructor::TryOneClusterFast(const AliMUONTrackParam &tr
   Double_t dX = cluster->GetX() - (trackParam.GetNonBendingCoor() + trackParam.GetNonBendingSlope() * dZ);
   Double_t dY = cluster->GetY() - (trackParam.GetBendingCoor() + trackParam.GetBendingSlope() * dZ);
   const TMatrixD& kParamCov = trackParam.GetCovariances();
-  Double_t errX2 = kParamCov(0,0) + dZ * dZ * kParamCov(1,1) + 2. * dZ * kParamCov(0,1);
-  Double_t errY2 = kParamCov(2,2) + dZ * dZ * kParamCov(3,3) + 2. * dZ * kParamCov(2,3);
+  Double_t errX2 = kParamCov(0,0) + dZ * dZ * kParamCov(1,1) + 2. * dZ * kParamCov(0,1) + cluster->GetErrX2();
+  Double_t errY2 = kParamCov(2,2) + dZ * dZ * kParamCov(3,3) + 2. * dZ * kParamCov(2,3) + cluster->GetErrY2();
 
-  Double_t dXmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errX2) +
+  Double_t dXmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errX2) +
                    GetRecoParam()->GetMaxNonBendingDistanceToTrack();
-  Double_t dYmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errY2) +
+  Double_t dYmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errY2) +
                   GetRecoParam()->GetMaxBendingDistanceToTrack();
   
   if (TMath::Abs(dX) > dXmax || TMath::Abs(dY) > dYmax) return kFALSE;
@@ -511,7 +678,7 @@ Double_t AliMUONVTrackReconstructor::TryTwoClustersFast(const AliMUONTrackParam
   
   // extrapolate linearly track parameters and covariances at the z position of the second cluster
   trackParamAtCluster2 = trackParamAtCluster1;
-  AliMUONTrackExtrap::LinearExtrapToZ(&trackParamAtCluster2, cluster2->GetZ());
+  AliMUONTrackExtrap::LinearExtrapToZCov(&trackParamAtCluster2, cluster2->GetZ());
   
   // set pointer to cluster2 into trackParamAtCluster2
   trackParamAtCluster2.SetClusterPtr(cluster2);
@@ -537,7 +704,7 @@ Double_t AliMUONVTrackReconstructor::TryTwoClustersFast(const AliMUONTrackParam
   // Compute chi2
   Double_t detX = sigma2X1 * sigma2X2 - covX1X2 * covX1X2;
   Double_t detY = sigma2Y1 * sigma2Y2 - covY1Y2 * covY1Y2;
-  if (detX == 0. || detY == 0.) return 1.e10;
+  if (detX == 0. || detY == 0.) return 2.*AliMUONTrack::MaxChi2();
   return   (dX1 * dX1 * sigma2X2 + dX2 * dX2 * sigma2X1 - 2. * dX1 * dX2 * covX1X2) / detX
         + (dY1 * dY1 * sigma2Y2 + dY2 * dY2 * sigma2Y1 - 2. * dY1 * dY2 * covY1Y2) / detY;
   
@@ -556,7 +723,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInChamber(AliMUONTrack &trac
   /// return kTRUE if new cluster(s) have been found (otherwise return kFALSE)
   AliDebug(1,Form("Enter FollowLinearTrackInChamber(1..) %d", nextChamber+1));
   
-  Double_t chi2WithOneCluster = 1.e10;
+  Double_t chi2WithOneCluster = AliMUONTrack::MaxChi2();
   Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
                                        GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
   Double_t bestChi2WithOneCluster = maxChi2WithOneCluster;
@@ -597,7 +764,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInChamber(AliMUONTrack &trac
     
     // try to add the current cluster accuratly
     extrapTrackParamAtCluster = trackParam;
-    AliMUONTrackExtrap::LinearExtrapToZ(&extrapTrackParamAtCluster, cluster->GetZ());
+    AliMUONTrackExtrap::LinearExtrapToZCov(&extrapTrackParamAtCluster, cluster->GetZ());
     chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCluster, cluster, extrapTrackParamAtCluster);
     
     // if good chi2 then consider to add cluster
@@ -618,6 +785,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInChamber(AliMUONTrack &trac
          extrapTrackParamAtCluster.SetRemovable(kFALSE);
        else extrapTrackParamAtCluster.SetRemovable(kTRUE);
        newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster,*cluster);
+       newTrack->SetGlobalChi2(trackCandidate.GetGlobalChi2()+chi2WithOneCluster);
        fNRecTracks++;
        
        // Printout for debuging
@@ -643,6 +811,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInChamber(AliMUONTrack &trac
        bestTrackParamAtCluster.SetRemovable(kFALSE);
       else bestTrackParamAtCluster.SetRemovable(kTRUE);
       trackCandidate.AddTrackParamAtCluster(bestTrackParamAtCluster,*(bestTrackParamAtCluster.GetClusterPtr()));
+      trackCandidate.SetGlobalChi2(trackCandidate.GetGlobalChi2()+bestChi2WithOneCluster);
       
       // Printout for debuging
       if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
@@ -689,8 +858,8 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac
     ch2 = 2*nextStation+1;
   }
   
-  Double_t chi2WithOneCluster = 1.e10;
-  Double_t chi2WithTwoClusters = 1.e10;
+  Double_t chi2WithOneCluster = AliMUONTrack::MaxChi2();
+  Double_t chi2WithTwoClusters = AliMUONTrack::MaxChi2();
   Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
                                        GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
   Double_t maxChi2WithTwoClusters = 4. * GetRecoParam()->GetSigmaCutForTracking() *
@@ -743,7 +912,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac
     
     // try to add the current cluster accuratly
     extrapTrackParamAtCluster2 = trackParam;
-    AliMUONTrackExtrap::LinearExtrapToZ(&extrapTrackParamAtCluster2, clusterCh2->GetZ());
+    AliMUONTrackExtrap::LinearExtrapToZCov(&extrapTrackParamAtCluster2, clusterCh2->GetZ());
     chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCluster2, clusterCh2, extrapTrackParamAtCluster2);
     
     // if good chi2 then try to attach a cluster in the other chamber too
@@ -794,6 +963,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac
            newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster1,*clusterCh1);
            extrapTrackParamAtCluster2.SetRemovable(kTRUE);
            newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster2,*clusterCh2);
+           newTrack->SetGlobalChi2(newTrack->GetGlobalChi2()+chi2WithTwoClusters);
            fNRecTracks++;
            
            // Tag clusterCh1 as used
@@ -827,6 +997,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac
            extrapTrackParamAtCluster2.SetRemovable(kFALSE);
          else extrapTrackParamAtCluster2.SetRemovable(kTRUE);
          newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster2,*clusterCh2);
+         newTrack->SetGlobalChi2(newTrack->GetGlobalChi2()+chi2WithOneCluster);
          fNRecTracks++;
          
          // Printout for debuging
@@ -857,7 +1028,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac
     }
     
     //Extrapolate trackCandidate to chamber "ch2"
-    AliMUONTrackExtrap::LinearExtrapToZ(&trackParam, AliMUONConstants::DefaultChamberZ(ch2));
+    AliMUONTrackExtrap::LinearExtrapToZCov(&trackParam, AliMUONConstants::DefaultChamberZ(ch2));
     
     // add MCS effect for next step
     AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
@@ -877,7 +1048,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac
       
       // try to add the current cluster accuratly
       extrapTrackParamAtCluster1 = trackParam;
-      AliMUONTrackExtrap::LinearExtrapToZ(&extrapTrackParamAtCluster1, clusterCh1->GetZ());
+      AliMUONTrackExtrap::LinearExtrapToZCov(&extrapTrackParamAtCluster1, clusterCh1->GetZ());
       chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCluster1, clusterCh1, extrapTrackParamAtCluster1);
       
       // if good chi2 then consider to add clusterCh1
@@ -899,6 +1070,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac
            extrapTrackParamAtCluster1.SetRemovable(kFALSE);
          else extrapTrackParamAtCluster1.SetRemovable(kTRUE);
          newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster1,*clusterCh1);
+         newTrack->SetGlobalChi2(newTrack->GetGlobalChi2()+chi2WithOneCluster);
          fNRecTracks++;
          
          // Printout for debuging
@@ -926,6 +1098,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac
       trackCandidate.AddTrackParamAtCluster(bestTrackParamAtCluster1,*(bestTrackParamAtCluster1.GetClusterPtr()));
       bestTrackParamAtCluster2.SetRemovable(kTRUE);
       trackCandidate.AddTrackParamAtCluster(bestTrackParamAtCluster2,*(bestTrackParamAtCluster2.GetClusterPtr()));
+      trackCandidate.SetGlobalChi2(trackCandidate.GetGlobalChi2()+bestChi2WithTwoClusters);
       
       // Printout for debuging
       if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
@@ -938,6 +1111,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac
        bestTrackParamAtCluster1.SetRemovable(kFALSE);
       else bestTrackParamAtCluster1.SetRemovable(kTRUE);
       trackCandidate.AddTrackParamAtCluster(bestTrackParamAtCluster1,*(bestTrackParamAtCluster1.GetClusterPtr()));
+      trackCandidate.SetGlobalChi2(trackCandidate.GetGlobalChi2()+bestChi2WithOneCluster);
       
       // Printout for debuging
       if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
@@ -1002,20 +1176,30 @@ void AliMUONVTrackReconstructor::Finalize()
 {
   /// Recompute track parameters and covariances at each attached cluster from those at the first one
   /// Set the label pointing to the corresponding MC track
+  /// Remove the track if finalization failed
   
-  AliMUONTrack *track;
+  AliMUONTrack *track, *nextTrack;
+  Bool_t trackRemoved = kFALSE;
   
   track = (AliMUONTrack*) fRecTracksPtr->First();
   while (track) {
     
-    FinalizeTrack(*track);
+    nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track);
     
-    track->FindMCLabel();
+    if (FinalizeTrack(*track)) track->FindMCLabel();
+    else {
+      fRecTracksPtr->Remove(track);
+      fNRecTracks--;
+      trackRemoved = kTRUE;
+    }
     
-    track = (AliMUONTrack*) fRecTracksPtr->After(track);
+    track = nextTrack;
     
   }
   
+  // compress array of tracks if needed
+  if (trackRemoved) fRecTracksPtr->Compress();
+  
 }
 
 //__________________________________________________________________________
@@ -1079,6 +1263,9 @@ void AliMUONVTrackReconstructor::EventReconstructTrigger(const AliMUONTriggerCir
       
       Float_t thetax = TMath::ATan2( x11 , z11 );
       Float_t thetay = TMath::ATan2( (y21-y11) , (z21-z11) );
+
+      CorrectThetaRange(thetax);
+      CorrectThetaRange(thetay);
       
       triggerTrack.SetX11(x11);
       triggerTrack.SetY11(y11);
@@ -1092,3 +1279,13 @@ void AliMUONVTrackReconstructor::EventReconstructTrigger(const AliMUONTriggerCir
   } // end of loop on Local Trigger
 }
 
+//__________________________________________________________________________
+void AliMUONVTrackReconstructor::CorrectThetaRange(Float_t& theta)
+{
+  /// The angles of the trigger tracks, obtained with ATan2, 
+  /// have values around +pi and -pi. On the contrary, the angles 
+  /// used in the tracker tracks have values around 0.
+  /// This function sets the same range for the trigger tracks angles.
+  if (theta < -TMath::PiOver2()) theta += TMath::Pi();
+  else if(theta > TMath::PiOver2()) theta -= TMath::Pi();
+}