- Added track parameters at DCA in AliESDMuonTrack
[u/mrichter/AliRoot.git] / MUON / AliMUONTrack.cxx
index adfd54c..718e2ea 100644 (file)
 #include "AliMUONTrack.h"
 
 #include "AliMUONTrackParam.h" 
-#include "AliMUONHitForRec.h" 
+#include "AliMUONRawClusterV2.h" 
 #include "AliMUONObjectPair.h" 
 #include "AliMUONConstants.h"
 #include "AliMUONTrackExtrap.h" 
 
 #include "AliLog.h"
+#include "AliESDMuonTrack.h"
+#include "AliESDMuonCluster.h"
 
 #include <TMath.h>
 #include <TMatrixD.h>
@@ -43,134 +45,131 @@ ClassImp(AliMUONTrack) // Class implementation in ROOT context
 //__________________________________________________________________________
 AliMUONTrack::AliMUONTrack()
   : TObject(),
-    fTrackParamAtVertex(),
-    fTrackParamAtHit(0x0),
-    fHitForRecAtHit(0x0),
-    fNTrackHits(0),
+    fTrackParamAtCluster(new TClonesArray("AliMUONTrackParam",10)),
     fFitWithVertex(kFALSE),
-    fVertex(0x0),
+    fVertexErrXY2(),
     fFitWithMCS(kFALSE),
-    fHitWeightsNonBending(0x0),
-    fHitWeightsBending(0x0),
+    fClusterWeightsNonBending(0x0),
+    fClusterWeightsBending(0x0),
     fGlobalChi2(-1.),
     fImproved(kFALSE),
     fMatchTrigger(-1),
     floTrgNum(-1),
     fChi2MatchTrigger(0.),
     fTrackID(0),
+    fTrackParamAtVertex(0x0),
     fHitsPatternInTrigCh(0),
     fLocalTrigger(0)
 {
   /// Default constructor
+  fTrackParamAtCluster->SetOwner(kTRUE);
+  fVertexErrXY2[0] = 0.;
+  fVertexErrXY2[1] = 0.;
 }
 
   //__________________________________________________________________________
 AliMUONTrack::AliMUONTrack(AliMUONObjectPair *segment)
   : TObject(),
-    fTrackParamAtVertex(),
-    fTrackParamAtHit(0x0),
-    fHitForRecAtHit(0x0),
-    fNTrackHits(0),
+    fTrackParamAtCluster(new TClonesArray("AliMUONTrackParam",10)),
     fFitWithVertex(kFALSE),
-    fVertex(0x0),
+    fVertexErrXY2(),
     fFitWithMCS(kFALSE),
-    fHitWeightsNonBending(0x0),
-    fHitWeightsBending(0x0),
+    fClusterWeightsNonBending(0x0),
+    fClusterWeightsBending(0x0),
     fGlobalChi2(0.),
     fImproved(kFALSE),
     fMatchTrigger(-1),
     floTrgNum(-1),    
     fChi2MatchTrigger(0.),
     fTrackID(0),
+    fTrackParamAtVertex(0x0),
     fHitsPatternInTrigCh(0),
     fLocalTrigger(0)
 {
-  /// Constructor from thw hitForRec's
-
-  fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
-  fTrackParamAtHit->SetOwner(kTRUE);
-  fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
-  fHitForRecAtHit->SetOwner(kTRUE);
+  /// Constructor from two clusters
+  fTrackParamAtCluster->SetOwner(kTRUE);
   
-  if (!segment) return; //AZ
+  fVertexErrXY2[0] = 0.;
+  fVertexErrXY2[1] = 0.;
   
-  // Pointers to hits from the segment
-  AliMUONHitForRec* hit1 = (AliMUONHitForRec*) segment->First();
-  AliMUONHitForRec* hit2 = (AliMUONHitForRec*) segment->Second();
+  // Pointers to clusters from the segment
+  AliMUONVCluster* cluster1 = (AliMUONVCluster*) segment->First();
+  AliMUONVCluster* cluster2 = (AliMUONVCluster*) segment->Second();
   
   // check sorting in -Z (spectro z<0)
-  if (hit1->GetZ() < hit2->GetZ()) {
-    hit1 = hit2;
-    hit2 = (AliMUONHitForRec*) segment->First();
+  if (cluster1->GetZ() < cluster2->GetZ()) {
+    cluster1 = cluster2;
+    cluster2 = (AliMUONVCluster*) segment->First();
   }
   
-  // order the hits into the track according to the station the segment belong to
-  //(the hit first attached is the one from which we will start the tracking procedure)
-  if (hit1->GetChamberNumber() == 8) {
-    AddTrackParamAtHit(0,hit1);
-    AddTrackParamAtHit(0,hit2);
+  // order the clusters into the track according to the station the segment belong to
+  //(the cluster first attached is the one from which we will start the tracking procedure)
+  AliMUONVCluster *firstCluster, *lastCluster;
+  if (cluster1->GetChamberId() == 8) {
+    firstCluster = cluster1;
+    lastCluster = cluster2;
   } else {
-    AddTrackParamAtHit(0,hit2);
-    AddTrackParamAtHit(0,hit1);
+    firstCluster = cluster2;
+    lastCluster = cluster1;
   }
   
-  AliMUONTrackParam* trackParamAtFirstHit = (AliMUONTrackParam*) fTrackParamAtHit->First();
-  AliMUONHitForRec* firstHit = trackParamAtFirstHit->GetHitForRecPtr();
-  AliMUONTrackParam* trackParamAtLastHit = (AliMUONTrackParam*) fTrackParamAtHit->Last();
-  AliMUONHitForRec* lastHit = trackParamAtLastHit->GetHitForRecPtr();
-  
-  
   // Compute track parameters
-  Double_t dZ = firstHit->GetZ() - lastHit->GetZ();
+  Double_t z1 = firstCluster->GetZ();
+  Double_t z2 = lastCluster->GetZ();
+  Double_t dZ = z1 - z2;
   // Non bending plane
-  Double_t nonBendingCoor1 = firstHit->GetNonBendingCoor();
-  Double_t nonBendingCoor2 = lastHit->GetNonBendingCoor();
+  Double_t nonBendingCoor1 = firstCluster->GetX();
+  Double_t nonBendingCoor2 = lastCluster->GetX();
   Double_t nonBendingSlope = (nonBendingCoor1 - nonBendingCoor2) / dZ;
   // Bending plane
-  Double_t bendingCoor1 = firstHit->GetBendingCoor();
-  Double_t bendingCoor2 = lastHit->GetBendingCoor();
+  Double_t bendingCoor1 = firstCluster->GetY();
+  Double_t bendingCoor2 = lastCluster->GetY();
   Double_t bendingSlope = (bendingCoor1 - bendingCoor2) / dZ;
   // Inverse bending momentum
-  Double_t bendingImpact = bendingCoor1 - firstHit->GetZ() * bendingSlope;
+  Double_t bendingImpact = bendingCoor1 - z1 * bendingSlope;
   Double_t inverseBendingMomentum = 1. / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact);
   
   
-  // Set track parameters at first hit
-  trackParamAtFirstHit->SetNonBendingCoor(nonBendingCoor1);
-  trackParamAtFirstHit->SetNonBendingSlope(nonBendingSlope);
-  trackParamAtFirstHit->SetBendingCoor(bendingCoor1);
-  trackParamAtFirstHit->SetBendingSlope(bendingSlope);
-  trackParamAtFirstHit->SetInverseBendingMomentum(inverseBendingMomentum);
+  // Set track parameters at first cluster
+  AliMUONTrackParam trackParamAtFirstCluster;
+  trackParamAtFirstCluster.SetZ(z1);
+  trackParamAtFirstCluster.SetNonBendingCoor(nonBendingCoor1);
+  trackParamAtFirstCluster.SetNonBendingSlope(nonBendingSlope);
+  trackParamAtFirstCluster.SetBendingCoor(bendingCoor1);
+  trackParamAtFirstCluster.SetBendingSlope(bendingSlope);
+  trackParamAtFirstCluster.SetInverseBendingMomentum(inverseBendingMomentum);
   
   
-  // Set track parameters at last hit
-  trackParamAtLastHit->SetNonBendingCoor(nonBendingCoor2);
-  trackParamAtLastHit->SetNonBendingSlope(nonBendingSlope);
-  trackParamAtLastHit->SetBendingCoor(bendingCoor2);
-  trackParamAtLastHit->SetBendingSlope(bendingSlope);
-  trackParamAtLastHit->SetInverseBendingMomentum(inverseBendingMomentum);
+  // Set track parameters at last cluster
+  AliMUONTrackParam trackParamAtLastCluster;
+  trackParamAtLastCluster.SetZ(z2);
+  trackParamAtLastCluster.SetNonBendingCoor(nonBendingCoor2);
+  trackParamAtLastCluster.SetNonBendingSlope(nonBendingSlope);
+  trackParamAtLastCluster.SetBendingCoor(bendingCoor2);
+  trackParamAtLastCluster.SetBendingSlope(bendingSlope);
+  trackParamAtLastCluster.SetInverseBendingMomentum(inverseBendingMomentum);
   
   
-  // Compute and set track parameters covariances at first hit
+  // Compute and set track parameters covariances at first cluster
   TMatrixD paramCov1(5,5);
   paramCov1.Zero();
   // Non bending plane
-  paramCov1(0,0) = firstHit->GetNonBendingReso2();
-  paramCov1(0,1) = firstHit->GetNonBendingReso2() / dZ;
+  paramCov1(0,0) = firstCluster->GetErrX2();
+  paramCov1(0,1) = firstCluster->GetErrX2() / dZ;
   paramCov1(1,0) = paramCov1(0,1);
-  paramCov1(1,1) = ( firstHit->GetNonBendingReso2() + lastHit->GetNonBendingReso2() ) / dZ / dZ;
+  paramCov1(1,1) = ( firstCluster->GetErrX2() + lastCluster->GetErrX2() ) / dZ / dZ;
   // Bending plane
-  paramCov1(2,2) = firstHit->GetBendingReso2();
-  paramCov1(2,3) = firstHit->GetBendingReso2() / dZ;
+  paramCov1(2,2) = firstCluster->GetErrY2();
+  paramCov1(2,3) = firstCluster->GetErrY2() / dZ;
   paramCov1(3,2) = paramCov1(2,3);
-  paramCov1(3,3) = ( firstHit->GetBendingReso2() + lastHit->GetBendingReso2() ) / dZ / dZ;
+  paramCov1(3,3) = ( firstCluster->GetErrY2() + lastCluster->GetErrY2() ) / dZ / dZ;
   // Inverse bending momentum (50% error)
   paramCov1(4,4) = 0.5*inverseBendingMomentum * 0.5*inverseBendingMomentum;
   // Set covariances
-  trackParamAtFirstHit->SetCovariances(paramCov1);
+  trackParamAtFirstCluster.SetCovariances(paramCov1);
   
   
-  // Compute and set track parameters covariances at last hit (as if the first hit did not exist)
+  // Compute and set track parameters covariances at last cluster (as if the first cluster did not exist)
   TMatrixD paramCov2(5,5);
   paramCov2.Zero();
   // Non bending plane
@@ -182,65 +181,160 @@ AliMUONTrack::AliMUONTrack(AliMUONObjectPair *segment)
   // Inverse bending momentum
   paramCov2(4,4) = paramCov1(4,4);
   // Set covariances
-  trackParamAtLastHit->SetCovariances(paramCov2);
+  trackParamAtLastCluster.SetCovariances(paramCov2);
+  
+  // Flag clusters as being removable
+  trackParamAtFirstCluster.SetRemovable(kTRUE);
+  trackParamAtLastCluster.SetRemovable(kTRUE);
+  
+  // Add track parameters at clusters
+  AddTrackParamAtCluster(trackParamAtFirstCluster,*firstCluster);
+  AddTrackParamAtCluster(trackParamAtLastCluster,*lastCluster);
+  
+}
+
+//__________________________________________________________________________
+AliMUONTrack::AliMUONTrack(AliESDMuonTrack &esdTrack)
+  : TObject(),
+    fTrackParamAtCluster(new TClonesArray("AliMUONTrackParam",10)),
+    fFitWithVertex(kFALSE),
+    fVertexErrXY2(),
+    fFitWithMCS(kFALSE),
+    fClusterWeightsNonBending(0x0),
+    fClusterWeightsBending(0x0),
+    fGlobalChi2(esdTrack.GetChi2()),
+    fImproved(kFALSE),
+    fMatchTrigger(esdTrack.GetMatchTrigger()),
+    floTrgNum(-1),
+    fChi2MatchTrigger(esdTrack.GetChi2MatchTrigger()),
+    fTrackID(0),
+    fTrackParamAtVertex(new AliMUONTrackParam()),
+    fHitsPatternInTrigCh(esdTrack.GetHitsPatternInTrigCh()),
+    fLocalTrigger(0)
+{
+  /// Constructor from ESD muon track
+  /// Compute track parameters and covariances at each cluster if available
+  /// or store parameters and covariances at first (fake) cluster only if not
+  
+  fTrackParamAtCluster->SetOwner(kTRUE);
   
+  fVertexErrXY2[0] = 0.;
+  fVertexErrXY2[1] = 0.;
   
-  // Flag first hit as being removable
-  trackParamAtFirstHit->SetRemovable(kTRUE);
+  // global info
+  SetLocalTrigger(esdTrack.LoCircuit(), esdTrack.LoStripX(), esdTrack.LoStripY(),
+                 esdTrack.LoDev(), esdTrack.LoLpt(), esdTrack.LoHpt());
   
-  // Flag last hit as being removable
-  trackParamAtLastHit->SetRemovable(kTRUE);
+  // track parameters at vertex
+  fTrackParamAtVertex->GetParamFrom(esdTrack);
+  
+  // track parameters at first cluster
+  AliMUONTrackParam param;
+  param.GetParamFromUncorrected(esdTrack);
+  param.GetCovFrom(esdTrack);
+  
+  // fill fTrackParamAtCluster with track parameters at each cluster if available
+  if(esdTrack.ClustersStored()) {
+    
+    AliMUONRawClusterV2 cluster;
+    
+    // loop over ESD clusters
+    AliESDMuonCluster *esdCluster = (AliESDMuonCluster*) esdTrack.GetClusters().First();
+    while (esdCluster) {
+      
+      // copy cluster information
+      cluster.SetUniqueID(esdCluster->GetUniqueID());
+      cluster.SetXYZ(esdCluster->GetX(), esdCluster->GetY(), esdCluster->GetZ());
+      cluster.SetErrXY(esdCluster->GetErrX(), esdCluster->GetErrY());
+      
+      // only set the Z parameter to avoid error in the AddTrackParamAtCluster(...) method
+      param.SetZ(cluster.GetZ());
+      
+      // add common track parameters at current cluster
+      AddTrackParamAtCluster(param, cluster, kTRUE);
+      
+      esdCluster = (AliESDMuonCluster*) esdTrack.GetClusters().After(esdCluster);
+    }
+    
+    // sort array of track parameters at clusters
+    fTrackParamAtCluster->Sort();
+    
+    // check that parameters stored in ESD are parameters at the most upstream cluster
+    // (convert Z position values to Float_t because of Double32_t in ESD track)
+    AliMUONTrackParam *firstTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->First();
+    if (((Float_t)firstTrackParam->GetZ()) != ((Float_t)esdTrack.GetZUncorrected())) {
+      
+      AliError("track parameters are not given at the most upstream stored cluster");
+      fTrackParamAtCluster->Clear("C");
+      
+    } else {
+      
+      // Compute track parameters and covariances at each cluster from info at the first one
+      UpdateCovTrackParamAtCluster();
+      
+    }
+    
+  }
+  
+  // fill fTrackParamAtCluster with track parameters at first (fake) cluster
+  // if first cluster not found or clusters not available
+  if (GetNClusters() == 0) {
+    
+    // get number of the first hit chamber (according to the MUONClusterMap if not empty)
+    Int_t firstCh = 0;
+    if (esdTrack.GetMuonClusterMap() != 0) while (!esdTrack.IsInMuonClusterMap(firstCh)) firstCh++;
+    else firstCh = AliMUONConstants::ChamberNumber(param.GetZ());
+    
+    // produce fake cluster at this chamber
+    AliMUONRawClusterV2 fakeCluster(firstCh, 0, 0);
+    fakeCluster.SetXYZ(param.GetNonBendingCoor(), param.GetBendingCoor(), param.GetZ());
+    fakeCluster.SetErrXY(0., 0.);
+    
+    // add track parameters at first (fake) cluster
+    AddTrackParamAtCluster(param, fakeCluster, kTRUE);
+    
+  }
   
 }
 
-  //__________________________________________________________________________
-AliMUONTrack::AliMUONTrack (const AliMUONTrack& track)
+//__________________________________________________________________________
+AliMUONTrack::AliMUONTrack(const AliMUONTrack& track)
   : TObject(track),
-    fTrackParamAtVertex(track.fTrackParamAtVertex),
-    fTrackParamAtHit(0x0),
-    fHitForRecAtHit(0x0),
-    fNTrackHits(track.fNTrackHits),
+    fTrackParamAtCluster(new TClonesArray("AliMUONTrackParam",10)),
     fFitWithVertex(track.fFitWithVertex),
-    fVertex(0x0),
+    fVertexErrXY2(),
     fFitWithMCS(track.fFitWithMCS),
-    fHitWeightsNonBending(0x0),
-    fHitWeightsBending(0x0),
+    fClusterWeightsNonBending(0x0),
+    fClusterWeightsBending(0x0),
     fGlobalChi2(track.fGlobalChi2),
     fImproved(track.fImproved),
     fMatchTrigger(track.fMatchTrigger),
     floTrgNum(track.floTrgNum),    
     fChi2MatchTrigger(track.fChi2MatchTrigger),
     fTrackID(track.fTrackID),
+    fTrackParamAtVertex(0x0),
     fHitsPatternInTrigCh(track.fHitsPatternInTrigCh),
     fLocalTrigger(track.fLocalTrigger)
 {
   ///copy constructor
-  Int_t maxIndex = 0;
   
   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
-  if (track.fTrackParamAtHit) {
-    maxIndex = (track.fTrackParamAtHit)->GetEntriesFast();
-    fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",maxIndex);
-    for (Int_t index = 0; index < maxIndex; index++) {
-      new ((*fTrackParamAtHit)[index]) AliMUONTrackParam(*(AliMUONTrackParam*)track.fTrackParamAtHit->At(index));
-    }
+  AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->First();
+  while (trackParamAtCluster) {
+    new ((*fTrackParamAtCluster)[GetNClusters()]) AliMUONTrackParam(*trackParamAtCluster);
+    trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->After(trackParamAtCluster);
   }
   
-  // necessary to make a copy of the objects and not only the pointers in TClonesArray.
-  if (track.fHitForRecAtHit) {
-    maxIndex = (track.fHitForRecAtHit)->GetEntriesFast();
-    fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",maxIndex);
-    for (Int_t index = 0; index < maxIndex; index++) {
-      new ((*fHitForRecAtHit)[index]) AliMUONHitForRec(*(AliMUONHitForRec*)track.fHitForRecAtHit->At(index));
-    }
-  }
+  // copy vertex resolution square used during the tracking procedure
+  fVertexErrXY2[0] = track.fVertexErrXY2[0];
+  fVertexErrXY2[1] = track.fVertexErrXY2[1];
   
-  // copy vertex used during the tracking procedure if any
-  if (track.fVertex) fVertex = new AliMUONHitForRec(*(track.fVertex));
+  // copy cluster weights matrices if any
+  if (track.fClusterWeightsNonBending) fClusterWeightsNonBending = new TMatrixD(*(track.fClusterWeightsNonBending));
+  if (track.fClusterWeightsBending) fClusterWeightsBending = new TMatrixD(*(track.fClusterWeightsBending));
   
-  // copy hit weights matrices if any
-  if (track.fHitWeightsNonBending) fHitWeightsNonBending = new TMatrixD(*(track.fHitWeightsNonBending));
-  if (track.fHitWeightsBending) fHitWeightsBending = new TMatrixD(*(track.fHitWeightsBending));
+  // copy track parameters at vertex if any
+  if (track.fTrackParamAtVertex) fTrackParamAtVertex = new AliMUONTrackParam(*(track.fTrackParamAtVertex));
   
 }
 
@@ -255,71 +349,48 @@ AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& track)
   // base class assignement
   TObject::operator=(track);
 
-  fTrackParamAtVertex = track.fTrackParamAtVertex;
-
-  Int_t maxIndex = 0;
-  
-  // necessary to make a copy of the objects and not only the pointers in TClonesArray.
-  if (track.fTrackParamAtHit) {
-    if (fTrackParamAtHit) fTrackParamAtHit->Clear();
-    else fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
-    maxIndex = (track.fTrackParamAtHit)->GetEntriesFast();
-    for (Int_t index = 0; index < maxIndex; index++) {
-      new ((*fTrackParamAtHit)[fTrackParamAtHit->GetEntriesFast()])
-       AliMUONTrackParam(*(AliMUONTrackParam*)(track.fTrackParamAtHit)->At(index));
-    }
-  } else if (fTrackParamAtHit) {
-    delete fTrackParamAtHit;
-    fTrackParamAtHit = 0x0;
-  }
-
   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
-  if (track.fHitForRecAtHit) {
-    if (fHitForRecAtHit) fHitForRecAtHit->Clear();
-    else fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
-    maxIndex = (track.fHitForRecAtHit)->GetEntriesFast();
-    for (Int_t index = 0; index < maxIndex; index++) {
-      new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()])
-       AliMUONHitForRec(*(AliMUONHitForRec*)(track.fHitForRecAtHit)->At(index));
-    }
-  } else if (fHitForRecAtHit) {
-    delete fHitForRecAtHit;
-    fHitForRecAtHit = 0x0;
+  fTrackParamAtCluster = new TClonesArray("AliMUONTrackParam",10);
+  AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->First();
+  while (trackParamAtCluster) {
+    new ((*fTrackParamAtCluster)[GetNClusters()]) AliMUONTrackParam(*trackParamAtCluster);
+    trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->After(trackParamAtCluster);
   }
   
-  // copy vertex used during the tracking procedure if any.
-  if (track.fVertex) {
-    if (fVertex) *fVertex = *(track.fVertex);
-    else fVertex = new AliMUONHitForRec(*(track.fVertex));
-  } else if (fVertex) {
-    delete fVertex;
-    fVertex = 0x0;
+  // copy cluster weights matrix if any
+  if (track.fClusterWeightsNonBending) {
+    if (fClusterWeightsNonBending) {
+      fClusterWeightsNonBending->ResizeTo(*(track.fClusterWeightsNonBending));
+      *fClusterWeightsNonBending = *(track.fClusterWeightsNonBending);
+    } else fClusterWeightsNonBending = new TMatrixD(*(track.fClusterWeightsNonBending));
+  } else if (fClusterWeightsNonBending) {
+    delete fClusterWeightsNonBending;
+    fClusterWeightsNonBending = 0x0;
   }
   
-  // copy hit weights matrix if any
-  if (track.fHitWeightsNonBending) {
-    if (fHitWeightsNonBending) {
-      fHitWeightsNonBending->ResizeTo(*(track.fHitWeightsNonBending));
-      *fHitWeightsNonBending = *(track.fHitWeightsNonBending);
-    } else fHitWeightsNonBending = new TMatrixD(*(track.fHitWeightsNonBending));
-  } else if (fHitWeightsNonBending) {
-    delete fHitWeightsNonBending;
-    fHitWeightsNonBending = 0x0;
+  // copy cluster weights matrix if any
+  if (track.fClusterWeightsBending) {
+    if (fClusterWeightsBending) {
+      fClusterWeightsBending->ResizeTo(*(track.fClusterWeightsBending));
+      *fClusterWeightsBending = *(track.fClusterWeightsBending);
+    } else fClusterWeightsBending = new TMatrixD(*(track.fClusterWeightsBending));
+  } else if (fClusterWeightsBending) {
+    delete fClusterWeightsBending;
+    fClusterWeightsBending = 0x0;
   }
   
-  // copy hit weights matrix if any
-  if (track.fHitWeightsBending) {
-    if (fHitWeightsBending) {
-      fHitWeightsBending->ResizeTo(*(track.fHitWeightsBending));
-      *fHitWeightsBending = *(track.fHitWeightsBending);
-    } else fHitWeightsBending = new TMatrixD(*(track.fHitWeightsBending));
-  } else if (fHitWeightsBending) {
-    delete fHitWeightsBending;
-    fHitWeightsBending = 0x0;
+  // copy track parameters at vertex if any
+  if (track.fTrackParamAtVertex) {
+    if (fTrackParamAtVertex) *fTrackParamAtVertex = *(track.fTrackParamAtVertex);
+    else fTrackParamAtVertex = new AliMUONTrackParam(*(track.fTrackParamAtVertex));
+  } else if (fTrackParamAtVertex) {
+    delete fTrackParamAtVertex;
+    fTrackParamAtVertex = 0x0;
   }
   
-  fNTrackHits         =  track.fNTrackHits;
   fFitWithVertex      =  track.fFitWithVertex;
+  fVertexErrXY2[0]    =  track.fVertexErrXY2[0];
+  fVertexErrXY2[1]    =  track.fVertexErrXY2[1];
   fFitWithMCS         =  track.fFitWithMCS;
   fGlobalChi2         =  track.fGlobalChi2;
   fImproved           =  track.fImproved;
@@ -337,178 +408,156 @@ AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& track)
 AliMUONTrack::~AliMUONTrack()
 {
   /// Destructor
-  delete fTrackParamAtHit;
-  delete fHitForRecAtHit;
-  delete fVertex;
-  delete fHitWeightsNonBending;
-  delete fHitWeightsBending;
+  delete fTrackParamAtCluster;
+  delete fClusterWeightsNonBending;
+  delete fClusterWeightsBending;
+  delete fTrackParamAtVertex;
 }
 
   //__________________________________________________________________________
 void AliMUONTrack::Clear(Option_t* opt)
 {
   /// Clear arrays
-  if ( fTrackParamAtHit ) fTrackParamAtHit->Clear(opt);
-  if ( fHitForRecAtHit ) fHitForRecAtHit->Clear(opt);
-  delete fVertex; fVertex = 0x0;
-  delete fHitWeightsNonBending; fHitWeightsNonBending = 0x0;
-  delete fHitWeightsBending; fHitWeightsBending = 0x0;
+  fTrackParamAtCluster->Clear(opt);
+  delete fClusterWeightsNonBending; fClusterWeightsNonBending = 0x0;
+  delete fClusterWeightsBending; fClusterWeightsBending = 0x0;
+  delete fTrackParamAtVertex; fTrackParamAtVertex = 0x0;
 }
 
   //__________________________________________________________________________
-void AliMUONTrack::AddTrackParamAtHit(const AliMUONTrackParam *trackParam, AliMUONHitForRec *hitForRec)
+void AliMUONTrack::AddTrackParamAtCluster(const AliMUONTrackParam &trackParam, AliMUONVCluster &cluster, Bool_t copy)
 {
-  /// Add TrackParamAtHit if "trackParam" != NULL
-  /// else create empty TrackParamAtHit and set the z position to the one of "hitForRec" if any
-  /// Update link to HitForRec if "hitForRec" != NULL
-  if (!fTrackParamAtHit) {
-    fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);  
-    fNTrackHits = 0;
-  }
-  AliMUONTrackParam* trackParamAtHit;
-  if (trackParam) {
-    trackParamAtHit = new ((*fTrackParamAtHit)[fNTrackHits]) AliMUONTrackParam(*trackParam);
-    if (hitForRec) {
-      if (hitForRec->GetZ() != trackParam->GetZ())
-        AliWarning("Added track parameters at a different z position than the one of the attached hit");
-    }
-  } else {
-    trackParamAtHit = new ((*fTrackParamAtHit)[fNTrackHits]) AliMUONTrackParam();
-    if (hitForRec) trackParamAtHit->SetZ(hitForRec->GetZ());
-  }
-  if (hitForRec) trackParamAtHit->SetHitForRecPtr(hitForRec);
-  fNTrackHits++;
-}
-
-  //__________________________________________________________________________
-void AliMUONTrack::RemoveTrackParamAtHit(AliMUONTrackParam *trackParam)
-{
-  /// Remove trackParam from the array of TrackParamAtHit
-  if (!fTrackParamAtHit) {
-    AliWarning("array fTrackParamAtHit does not exist");
+  /// Copy given track parameters into a new TrackParamAtCluster
+  /// Link parameters with the associated cluster
+  /// If copy=kTRUE: the cluster is copied then passed the trackParam which become its owner 
+  ///     otherwise: make sure to do not delete the cluster until it is used by the track
+  
+  // check chamber ID of the associated cluster
+  if (cluster.GetChamberId() < 0 || cluster.GetChamberId() > AliMUONConstants::NTrackingCh()) {
+    AliError(Form("Chamber ID of the associated cluster is not valid (ChamberId=%d)",cluster.GetChamberId()));
     return;
   }
   
-  if (!fTrackParamAtHit->Remove(trackParam)) {
-    AliWarning("object to remove does not exist in array fTrackParamAtHit");
+  // check whether track parameters are given at the correct cluster z position
+  if (cluster.GetZ() != trackParam.GetZ()) {
+    AliError("track parameters are given at a different z position than the one of the associated cluster");
     return;
   }
   
-  fTrackParamAtHit->Compress();
-  fNTrackHits--;
+  // add parameters to the array of track parameters
+  AliMUONTrackParam* trackParamAtCluster = new ((*fTrackParamAtCluster)[GetNClusters()]) AliMUONTrackParam(trackParam);
+  
+  // link parameters with the associated cluster or its copy
+  if (copy) {
+    AliMUONVCluster *clusterCopy = static_cast<AliMUONVCluster*>(cluster.Clone());
+    trackParamAtCluster->SetClusterPtr(clusterCopy, kTRUE);
+  } else trackParamAtCluster->SetClusterPtr(&cluster);
 }
 
   //__________________________________________________________________________
-void AliMUONTrack::AddHitForRecAtHit(const AliMUONHitForRec *hitForRec) 
+void AliMUONTrack::RemoveTrackParamAtCluster(AliMUONTrackParam *trackParam)
 {
-  /// Add hitForRec to the array of hitForRec at hit
-  if (!fHitForRecAtHit)
-    fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10); 
-  
-  if (!hitForRec)
-    AliFatal("AliMUONTrack::AddHitForRecAtHit: hitForRec == NULL");
+  /// Remove trackParam from the array of TrackParamAtCluster
+  if (!fTrackParamAtCluster->Remove(trackParam)) {
+    AliWarning("object to remove does not exist in array fTrackParamAtCluster");
+    return;
+  }
   
-  new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()]) AliMUONHitForRec(*hitForRec);
+  fTrackParamAtCluster->Compress();
 }
 
   //__________________________________________________________________________
-void AliMUONTrack::UpdateTrackParamAtHit()
+void AliMUONTrack::UpdateTrackParamAtCluster()
 {
-  /// Update track parameters at each attached hit
+  /// Update track parameters at each attached cluster
   
-  if (fNTrackHits == 0) {
-    AliWarning("no hit attached to the track");
+  if (GetNClusters() == 0) {
+    AliWarning("no cluster attached to the track");
     return;
   }
   
-  Double_t z;
-  AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtHit->First();
-  AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->After(startingTrackParam);
-  while (trackParamAtHit) {
-    
-    // save current z
-    z = trackParamAtHit->GetZ();
+  AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->First();
+  AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(startingTrackParam);
+  while (trackParamAtCluster) {
     
     // reset track parameters and their covariances
-    trackParamAtHit->SetParameters(startingTrackParam->GetParameters());
-    trackParamAtHit->SetZ(startingTrackParam->GetZ());
+    trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
+    trackParamAtCluster->SetZ(startingTrackParam->GetZ());
     
     // extrapolation to the given z
-    AliMUONTrackExtrap::ExtrapToZ(trackParamAtHit, z);
+    AliMUONTrackExtrap::ExtrapToZ(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ());
     
     // prepare next step
-    startingTrackParam = trackParamAtHit;
-    trackParamAtHit = (AliMUONTrackParam*) (fTrackParamAtHit->After(trackParamAtHit));
+    startingTrackParam = trackParamAtCluster;
+    trackParamAtCluster = (AliMUONTrackParam*) (fTrackParamAtCluster->After(trackParamAtCluster));
   }
 
 }
 
   //__________________________________________________________________________
-void AliMUONTrack::UpdateCovTrackParamAtHit()
+void AliMUONTrack::UpdateCovTrackParamAtCluster()
 {
-  /// Update track parameters and their covariances at each attached hit
+  /// Update track parameters and their covariances at each attached cluster
+  /// Include effects of multiple scattering in chambers
   
-  if (fNTrackHits == 0) {
-    AliWarning("no hit attached to the track");
+  if (GetNClusters() == 0) {
+    AliWarning("no cluster attached to the track");
     return;
   }
   
-  Double_t z;
-  AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtHit->First();
-  AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->After(startingTrackParam);
-  while (trackParamAtHit) {
-    
-    // save current z
-    z = trackParamAtHit->GetZ();
+  AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->First();
+  AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(startingTrackParam);
+  Int_t expectedChamber = startingTrackParam->GetClusterPtr()->GetChamberId() + 1;
+  Int_t currentChamber;
+  while (trackParamAtCluster) {
     
     // reset track parameters and their covariances
-    trackParamAtHit->SetParameters(startingTrackParam->GetParameters());
-    trackParamAtHit->SetZ(startingTrackParam->GetZ());
-    trackParamAtHit->SetCovariances(startingTrackParam->GetCovariances());
+    trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
+    trackParamAtCluster->SetZ(startingTrackParam->GetZ());
+    trackParamAtCluster->SetCovariances(startingTrackParam->GetCovariances());
     
-    // extrapolation to the given z
-    AliMUONTrackExtrap::ExtrapToZCov(trackParamAtHit, z);
+    // add MCS effect
+    AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
+    
+    // add MCS in missing chambers if any
+    currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
+    while (currentChamber > expectedChamber) {
+      // extrapolation to the missing chamber
+      AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, AliMUONConstants::DefaultChamberZ(expectedChamber));
+      // add MCS effect
+      AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
+      expectedChamber++;
+    }
+    
+    // extrapolation to the z of the current cluster
+    AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ());
     
     // prepare next step
-    startingTrackParam = trackParamAtHit;
-    trackParamAtHit = (AliMUONTrackParam*) (fTrackParamAtHit->After(trackParamAtHit));
+    expectedChamber = currentChamber + 1;
+    startingTrackParam = trackParamAtCluster;
+    trackParamAtCluster = (AliMUONTrackParam*) (fTrackParamAtCluster->After(trackParamAtCluster));
   }
-
-}
-
-  //__________________________________________________________________________
-void AliMUONTrack::SetVertex(const AliMUONHitForRec* vertex)
-{
-  /// Set the vertex used during the tracking procedure
-  if (!fVertex) fVertex = new AliMUONHitForRec(*vertex);
-  else *fVertex = *vertex;
+  
 }
 
-
   //__________________________________________________________________________
 Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS)
 {
-  /// Compute the removable hit contribution to the chi2 of the track
+  /// Compute each cluster contribution to the chi2 of the track
   /// accounting for multiple scattering or not according to the flag
-  /// - Also recompute the weight matrices of the attached hits if accountForMCS=kTRUE
-  /// - Assume that track parameters at each hit are corrects
+  /// - Also recompute the weight matrices of the attached clusters if accountForMCS=kTRUE
+  /// - Assume that track parameters at each cluster are corrects
   /// - Return kFALSE if computation failed
   
-  // Check hits (if the first one exist, assume that the other ones exit too!)
-  AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->First();
-  if (!trackParamAtHit->GetHitForRecPtr()) {
-    AliWarning("hit is missing");
-    return kFALSE;
-  }
-  
   if (accountForMCS) { // Compute local chi2 taking into account multiple scattering effects
       
     // Compute MCS covariance matrix only once
-    TMatrixD mcsCovariances(fNTrackHits,fNTrackHits);
+    Int_t nClusters = GetNClusters();
+    TMatrixD mcsCovariances(nClusters,nClusters);
     ComputeMCSCovariances(mcsCovariances);
     
-    // Make sure hit weights are consistent with following calculations
-    if (!ComputeHitWeights(&mcsCovariances)) {
+    // Make sure cluster weights are consistent with following calculations
+    if (!ComputeClusterWeights(&mcsCovariances)) {
       AliWarning("cannot take into account the multiple scattering effects");
       return ComputeLocalChi2(kFALSE);
     }
@@ -517,64 +566,67 @@ Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS)
     Double_t globalChi2 = ComputeGlobalChi2(kTRUE);
     if (globalChi2 < 0.) return kFALSE;
     
-    // Loop over removable hits and compute their local chi2
-    AliMUONTrackParam* trackParamAtHit1;
-    AliMUONHitForRec *hitForRec, *discardedHit;
-    Int_t hitNumber1, hitNumber2, currentHitNumber1, currentHitNumber2;
-    TMatrixD hitWeightsNB(fNTrackHits-1,fNTrackHits-1);
-    TMatrixD hitWeightsB(fNTrackHits-1,fNTrackHits-1);
-    Double_t *dX = new Double_t[fNTrackHits-1];
-    Double_t *dY = new Double_t[fNTrackHits-1];
+    // Loop over removable clusters and compute their local chi2
+    AliMUONTrackParam* trackParamAtCluster1;
+    AliMUONVCluster *cluster, *discardedCluster;
+    Int_t iCluster1, iCluster2, iCurrentCluster1, iCurrentCluster2;
+    TMatrixD clusterWeightsNB(nClusters-1,nClusters-1);
+    TMatrixD clusterWeightsB(nClusters-1,nClusters-1);
+    Double_t *dX = new Double_t[nClusters-1];
+    Double_t *dY = new Double_t[nClusters-1];
     Double_t globalChi2b;
-    while (trackParamAtHit) {
+    AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->First();
+    while (trackParamAtCluster) {
       
-      discardedHit = trackParamAtHit->GetHitForRecPtr();
+      discardedCluster = trackParamAtCluster->GetClusterPtr();
       
-      // Recompute hit weights without the current hit
-      if (!ComputeHitWeights(hitWeightsNB, hitWeightsB, &mcsCovariances, discardedHit)) {
+      // Recompute cluster weights without the current cluster
+      if (!ComputeClusterWeights(clusterWeightsNB, clusterWeightsB, &mcsCovariances, discardedCluster)) {
        AliWarning("cannot take into account the multiple scattering effects");
-       ComputeLocalChi2(kFALSE);
+       delete [] dX;
+       delete [] dY;
+       return ComputeLocalChi2(kFALSE);
       }
       
-      // Compute track chi2 without the current hit
+      // Compute track chi2 without the current cluster
       globalChi2b = 0.;
-      currentHitNumber1 = 0;
-      for (hitNumber1 = 0; hitNumber1 < fNTrackHits ; hitNumber1++) { 
-       trackParamAtHit1 = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber1);
-       hitForRec = trackParamAtHit1->GetHitForRecPtr();
+      iCurrentCluster1 = 0;
+      for (iCluster1 = 0; iCluster1 < nClusters ; iCluster1++) { 
+       trackParamAtCluster1 = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1);
+       cluster = trackParamAtCluster1->GetClusterPtr();
         
-        if (hitForRec == discardedHit) continue;
+        if (cluster == discardedCluster) continue;
         
         // Compute and save residuals
-       dX[currentHitNumber1] = hitForRec->GetNonBendingCoor() - trackParamAtHit1->GetNonBendingCoor();
-       dY[currentHitNumber1] = hitForRec->GetBendingCoor() - trackParamAtHit1->GetBendingCoor();
+       dX[iCurrentCluster1] = cluster->GetX() - trackParamAtCluster1->GetNonBendingCoor();
+       dY[iCurrentCluster1] = cluster->GetY() - trackParamAtCluster1->GetBendingCoor();
         
-        currentHitNumber2 = 0;
-       for (hitNumber2 = 0; hitNumber2 < hitNumber1; hitNumber2++) {
-         hitForRec = ((AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber2))->GetHitForRecPtr();
+        iCurrentCluster2 = 0;
+       for (iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) {
+         cluster = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
           
-          if (hitForRec == discardedHit) continue;
+          if (cluster == discardedCluster) continue;
           
           // Add contribution from covariances
-          globalChi2b += (hitWeightsNB(currentHitNumber1, currentHitNumber2) +
-                         hitWeightsNB(currentHitNumber2, currentHitNumber1)) * dX[currentHitNumber1] * dX[currentHitNumber2] +
-                        (hitWeightsB(currentHitNumber1, currentHitNumber2) +
-                         hitWeightsB(currentHitNumber2, currentHitNumber1)) * dY[currentHitNumber1] * dY[currentHitNumber2];
+          globalChi2b += (clusterWeightsNB(iCurrentCluster1, iCurrentCluster2) +
+                         clusterWeightsNB(iCurrentCluster2, iCurrentCluster1)) * dX[iCurrentCluster1] * dX[iCurrentCluster2] +
+                        (clusterWeightsB(iCurrentCluster1, iCurrentCluster2) +
+                         clusterWeightsB(iCurrentCluster2, iCurrentCluster1)) * dY[iCurrentCluster1] * dY[iCurrentCluster2];
           
-          currentHitNumber2++;
+          iCurrentCluster2++;
        }
         
         // Add contribution from variances
-       globalChi2b += hitWeightsNB(currentHitNumber1, currentHitNumber1) * dX[currentHitNumber1] * dX[currentHitNumber1] +
-                      hitWeightsB(currentHitNumber1, currentHitNumber1) * dY[currentHitNumber1] * dY[currentHitNumber1];
+       globalChi2b += clusterWeightsNB(iCurrentCluster1, iCurrentCluster1) * dX[iCurrentCluster1] * dX[iCurrentCluster1] +
+                      clusterWeightsB(iCurrentCluster1, iCurrentCluster1) * dY[iCurrentCluster1] * dY[iCurrentCluster1];
        
-        currentHitNumber1++;
+        iCurrentCluster1++;
       }
 
       // Set local chi2
-      trackParamAtHit->SetLocalChi2(globalChi2 - globalChi2b);
+      trackParamAtCluster->SetLocalChi2(globalChi2 - globalChi2b);
       
-      trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->After(trackParamAtHit);
+      trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(trackParamAtCluster);
     }
     
     delete [] dX;
@@ -582,20 +634,21 @@ Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS)
     
   } else { // without multiple scattering effects
     
-    AliMUONHitForRec *discardedHit;
+    AliMUONVCluster *discardedCluster;
     Double_t dX, dY;
-    while (trackParamAtHit) {
+    AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->First();
+    while (trackParamAtCluster) {
       
-      discardedHit = trackParamAtHit->GetHitForRecPtr();
+      discardedCluster = trackParamAtCluster->GetClusterPtr();
       
       // Compute residuals
-      dX = discardedHit->GetNonBendingCoor() - trackParamAtHit->GetNonBendingCoor();
-      dY = discardedHit->GetBendingCoor() - trackParamAtHit->GetBendingCoor();
+      dX = discardedCluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
+      dY = discardedCluster->GetY() - trackParamAtCluster->GetBendingCoor();
       
       // Set local chi2
-      trackParamAtHit->SetLocalChi2(dX * dX / discardedHit->GetNonBendingReso2() + dY * dY / discardedHit->GetBendingReso2());
+      trackParamAtCluster->SetLocalChi2(dX * dX / discardedCluster->GetErrX2() + dY * dY / discardedCluster->GetErrY2());
     
-    trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->After(trackParamAtHit);
+      trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(trackParamAtCluster);
     }
   
   }
@@ -608,61 +661,57 @@ Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS)
 Double_t AliMUONTrack::ComputeGlobalChi2(Bool_t accountForMCS)
 {
   /// Compute the chi2 of the track accounting for multiple scattering or not according to the flag
-  /// - Assume that track parameters at each hit are corrects
-  /// - Assume the hits weights matrices are corrects
+  /// - Assume that track parameters at each cluster are corrects
+  /// - Assume the cluster weights matrices are corrects
   /// - Return negative value if chi2 computation failed
   
-  // Check hits (if the first one exist, assume that the other ones exit too!)
-  AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->First();
-  if (!trackParamAtHit->GetHitForRecPtr()) {
-    AliWarning("hit is missing");
-    return -1.;
-  }
-  
   Double_t chi2 = 0.;
   
   if (accountForMCS) {
     
     // Check the weight matrices. If weight matrices are not available compute chi2 without MCS
-    if (!fHitWeightsNonBending || !fHitWeightsBending) {
-      AliWarning("hit weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
+    if (!fClusterWeightsNonBending || !fClusterWeightsBending) {
+      AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
       return ComputeGlobalChi2(kFALSE);
     }
-    if (fHitWeightsNonBending->GetNrows() != fNTrackHits || fHitWeightsBending->GetNcols() != fNTrackHits) {
-      AliWarning("hit weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
+    Int_t nClusters = GetNClusters();
+    if (fClusterWeightsNonBending->GetNrows() != nClusters || fClusterWeightsBending->GetNcols() != nClusters) {
+      AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
       return ComputeGlobalChi2(kFALSE);
     }
     
     // Compute chi2
-    AliMUONHitForRec *hitForRec;
-    Double_t *dX = new Double_t[fNTrackHits];
-    Double_t *dY = new Double_t[fNTrackHits];
-    Int_t hitNumber1, hitNumber2;
-    for (hitNumber1 = 0; hitNumber1 < fNTrackHits ; hitNumber1++) { 
-      trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber1);
-      hitForRec = trackParamAtHit->GetHitForRecPtr();
-      dX[hitNumber1] = hitForRec->GetNonBendingCoor() - trackParamAtHit->GetNonBendingCoor();
-      dY[hitNumber1] = hitForRec->GetBendingCoor() - trackParamAtHit->GetBendingCoor();
-      for (hitNumber2 = 0; hitNumber2 < hitNumber1; hitNumber2++) {
-        chi2 += ((*fHitWeightsNonBending)(hitNumber1, hitNumber2) + (*fHitWeightsNonBending)(hitNumber2, hitNumber1)) * dX[hitNumber1] * dX[hitNumber2] +
-               ((*fHitWeightsBending)(hitNumber1, hitNumber2) + (*fHitWeightsBending)(hitNumber2, hitNumber1)) * dY[hitNumber1] * dY[hitNumber2];
+    AliMUONVCluster *cluster;
+    Double_t *dX = new Double_t[nClusters];
+    Double_t *dY = new Double_t[nClusters];
+    AliMUONTrackParam* trackParamAtCluster;
+    for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) { 
+      trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1);
+      cluster = trackParamAtCluster->GetClusterPtr();
+      dX[iCluster1] = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
+      dY[iCluster1] = cluster->GetY() - trackParamAtCluster->GetBendingCoor();
+      for (Int_t iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) {
+        chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster2) + (*fClusterWeightsNonBending)(iCluster2, iCluster1)) * dX[iCluster1] * dX[iCluster2] +
+               ((*fClusterWeightsBending)(iCluster1, iCluster2) + (*fClusterWeightsBending)(iCluster2, iCluster1)) * dY[iCluster1] * dY[iCluster2];
       }
-      chi2 += ((*fHitWeightsNonBending)(hitNumber1, hitNumber1) * dX[hitNumber1] * dX[hitNumber1]) +
-             ((*fHitWeightsBending)(hitNumber1, hitNumber1) * dY[hitNumber1] * dY[hitNumber1]);
+      chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster1) * dX[iCluster1] * dX[iCluster1]) +
+             ((*fClusterWeightsBending)(iCluster1, iCluster1) * dY[iCluster1] * dY[iCluster1]);
     }
     delete [] dX;
     delete [] dY;
     
   } else {
     
-    AliMUONHitForRec *hitForRec;
+    AliMUONVCluster *cluster;
     Double_t dX, dY;
-    for (Int_t hitNumber = 0; hitNumber < fNTrackHits ; hitNumber++) { 
-      trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber);
-      hitForRec = trackParamAtHit->GetHitForRecPtr();
-      dX = hitForRec->GetNonBendingCoor() - trackParamAtHit->GetNonBendingCoor();
-      dY = hitForRec->GetBendingCoor() - trackParamAtHit->GetBendingCoor();
-      chi2 += dX * dX / hitForRec->GetNonBendingReso2() + dY * dY / hitForRec->GetBendingReso2();
+    AliMUONTrackParam* trackParamAtCluster;
+    Int_t nClusters = GetNClusters();
+    for (Int_t iCluster = 0; iCluster < nClusters ; iCluster++) { 
+      trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster);
+      cluster = trackParamAtCluster->GetClusterPtr();
+      dX = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
+      dY = cluster->GetY() - trackParamAtCluster->GetBendingCoor();
+      chi2 += dX * dX / cluster->GetErrX2() + dY * dY / cluster->GetErrY2();
     }
     
   }
@@ -672,114 +721,109 @@ Double_t AliMUONTrack::ComputeGlobalChi2(Bool_t accountForMCS)
 }
 
   //__________________________________________________________________________
-Bool_t AliMUONTrack::ComputeHitWeights(TMatrixD* mcsCovariances)
+Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD* mcsCovariances)
 {
-  /// Compute the weight matrices of the attached hits, in non bending and bending direction,
-  /// accounting for multiple scattering correlations and hits resolution
+  /// Compute the weight matrices of the attached clusters, in non bending and bending direction,
+  /// accounting for multiple scattering correlations and cluster resolution
   /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
-  /// - Assume that track parameters at each hit are corrects
+  /// - Assume that track parameters at each cluster are corrects
   /// - Return kFALSE if computation failed
   
   // Alocate memory
-  if (!fHitWeightsNonBending) fHitWeightsNonBending = new TMatrixD(fNTrackHits,fNTrackHits);
-  if (!fHitWeightsBending) fHitWeightsBending = new TMatrixD(fNTrackHits,fNTrackHits);
-  
-  // Check hits (if the first one exist, assume that the other ones exit too!)
-  if (!((AliMUONTrackParam*) fTrackParamAtHit->First())->GetHitForRecPtr()) {
-    AliWarning("hit is missing");
-    fHitWeightsNonBending->ResizeTo(0,0);
-    fHitWeightsBending->ResizeTo(0,0);
-    return kFALSE;
-  }
+  Int_t nClusters = GetNClusters();
+  if (!fClusterWeightsNonBending) fClusterWeightsNonBending = new TMatrixD(nClusters,nClusters);
+  if (!fClusterWeightsBending) fClusterWeightsBending = new TMatrixD(nClusters,nClusters);
   
   // Compute weights matrices
-  if (!ComputeHitWeights(*fHitWeightsNonBending, *fHitWeightsBending, mcsCovariances)) return kFALSE;
+  if (!ComputeClusterWeights(*fClusterWeightsNonBending, *fClusterWeightsBending, mcsCovariances)) return kFALSE;
   
   return kTRUE;
   
 }
 
   //__________________________________________________________________________
-Bool_t AliMUONTrack::ComputeHitWeights(TMatrixD& hitWeightsNB, TMatrixD& hitWeightsB, TMatrixD* mcsCovariances, AliMUONHitForRec* discardedHit) const
+Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD& clusterWeightsNB, TMatrixD& clusterWeightsB,
+                                          TMatrixD* mcsCovariances, AliMUONVCluster* discardedCluster) const
 {
   /// Compute the weight matrices, in non bending and bending direction,
-  /// of the other attached hits assuming the discarded one does not exist
-  /// accounting for multiple scattering correlations and hits resolution
+  /// of the other attached clusters assuming the discarded one does not exist
+  /// accounting for multiple scattering correlations and cluster resolution
   /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
   /// - Return kFALSE if computation failed
   
   // Check MCS covariance matrix and recompute it if need
+  Int_t nClusters = GetNClusters();
   Bool_t deleteMCSCov = kFALSE;
   if (!mcsCovariances) {
-    mcsCovariances = new TMatrixD(fNTrackHits,fNTrackHits);
+    mcsCovariances = new TMatrixD(nClusters,nClusters);
     deleteMCSCov = kTRUE;
     ComputeMCSCovariances(*mcsCovariances);
   }
   
   // Resize the weights matrices; alocate memory
-  if (discardedHit) {
-    hitWeightsNB.ResizeTo(fNTrackHits-1,fNTrackHits-1);
-    hitWeightsB.ResizeTo(fNTrackHits-1,fNTrackHits-1);
+  if (discardedCluster) {
+    clusterWeightsNB.ResizeTo(nClusters-1,nClusters-1);
+    clusterWeightsB.ResizeTo(nClusters-1,nClusters-1);
   } else {
-    hitWeightsNB.ResizeTo(fNTrackHits,fNTrackHits);
-    hitWeightsB.ResizeTo(fNTrackHits,fNTrackHits);
+    clusterWeightsNB.ResizeTo(nClusters,nClusters);
+    clusterWeightsB.ResizeTo(nClusters,nClusters);
   }
   
   // Define variables
-  AliMUONHitForRec *hitForRec1, *hitForRec2;
-  Int_t currentHitNumber1, currentHitNumber2;
+  AliMUONVCluster *cluster1, *cluster2;
+  Int_t iCurrentCluster1, iCurrentCluster2;
   
   // Compute the covariance matrices
-  currentHitNumber1 = 0;
-  for (Int_t hitNumber1 = 0; hitNumber1 < fNTrackHits; hitNumber1++) { 
-    hitForRec1 = ((AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber1))->GetHitForRecPtr();
+  iCurrentCluster1 = 0;
+  for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) { 
+    cluster1 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1))->GetClusterPtr();
     
-    if (hitForRec1 == discardedHit) continue;
+    if (cluster1 == discardedCluster) continue;
     
-    // Loop over next hits
-    currentHitNumber2 = currentHitNumber1;
-    for (Int_t hitNumber2 = hitNumber1; hitNumber2 < fNTrackHits; hitNumber2++) {
-      hitForRec2 = ((AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber2))->GetHitForRecPtr();
+    // Loop over next clusters
+    iCurrentCluster2 = iCurrentCluster1;
+    for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) {
+      cluster2 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
       
-      if (hitForRec2 == discardedHit) continue;
+      if (cluster2 == discardedCluster) continue;
       
       // Fill with MCS covariances
-      hitWeightsNB(currentHitNumber1, currentHitNumber2) = (*mcsCovariances)(hitNumber1,hitNumber2);
+      clusterWeightsNB(iCurrentCluster1, iCurrentCluster2) = (*mcsCovariances)(iCluster1,iCluster2);
       
       // Equal contribution from multiple scattering in non bending and bending directions
-      hitWeightsB(currentHitNumber1, currentHitNumber2) = hitWeightsNB(currentHitNumber1, currentHitNumber2);
+      clusterWeightsB(iCurrentCluster1, iCurrentCluster2) = clusterWeightsNB(iCurrentCluster1, iCurrentCluster2);
       
-      // Add contribution from hit resolution to diagonal element and symmetrize the matrix
-      if (currentHitNumber1 == currentHitNumber2) {
+      // Add contribution from cluster resolution to diagonal element and symmetrize the matrix
+      if (iCurrentCluster1 == iCurrentCluster2) {
        
        // In non bending plane
-        hitWeightsNB(currentHitNumber1, currentHitNumber1) += hitForRec1->GetNonBendingReso2();
+        clusterWeightsNB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrX2();
        // In bending plane
-       hitWeightsB(currentHitNumber1, currentHitNumber1) += hitForRec1->GetBendingReso2();
+       clusterWeightsB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrY2();
        
       } else {
        
        // In non bending plane
-       hitWeightsNB(currentHitNumber2, currentHitNumber1) = hitWeightsNB(currentHitNumber1, currentHitNumber2);
+       clusterWeightsNB(iCurrentCluster2, iCurrentCluster1) = clusterWeightsNB(iCurrentCluster1, iCurrentCluster2);
        // In bending plane
-       hitWeightsB(currentHitNumber2, currentHitNumber1) = hitWeightsB(currentHitNumber1, currentHitNumber2);
+       clusterWeightsB(iCurrentCluster2, iCurrentCluster1) = clusterWeightsB(iCurrentCluster1, iCurrentCluster2);
        
       }
       
-      currentHitNumber2++;
+      iCurrentCluster2++;
     }
     
-    currentHitNumber1++;
+    iCurrentCluster1++;
   }
     
   // Inversion of covariance matrices to get the weights
-  if (hitWeightsNB.Determinant() != 0 && hitWeightsB.Determinant() != 0) {
-    hitWeightsNB.Invert();
-    hitWeightsB.Invert();
+  if (clusterWeightsNB.Determinant() != 0 && clusterWeightsB.Determinant() != 0) {
+    clusterWeightsNB.Invert();
+    clusterWeightsB.Invert();
   } else {
     AliWarning(" Determinant = 0");
-    hitWeightsNB.ResizeTo(0,0);
-    hitWeightsB.ResizeTo(0,0);
+    clusterWeightsNB.ResizeTo(0,0);
+    clusterWeightsB.ResizeTo(0,0);
     if(deleteMCSCov) delete mcsCovariances;
     return kFALSE;
   }
@@ -794,40 +838,39 @@ Bool_t AliMUONTrack::ComputeHitWeights(TMatrixD& hitWeightsNB, TMatrixD& hitWeig
 void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const
 {
   /// Compute the multiple scattering covariance matrix
-  /// (assume that track parameters at each hit are corrects)
+  /// (assume that track parameters at each cluster are corrects)
   
   // Reset the size of the covariance matrix if needed
-  if (mcsCovariances.GetNrows() != fNTrackHits) mcsCovariances.ResizeTo(fNTrackHits,fNTrackHits);
+  Int_t nClusters = GetNClusters();
+  if (mcsCovariances.GetNrows() != nClusters) mcsCovariances.ResizeTo(nClusters,nClusters);
   
   // Define variables
   Int_t nChambers = AliMUONConstants::NTrackingCh();
-  AliMUONTrackParam* trackParamAtHit;
-  AliMUONHitForRec *hitForRec;
+  AliMUONTrackParam* trackParamAtCluster;
   AliMUONTrackParam extrapTrackParam;
   Int_t currentChamber = 0, expectedChamber = 0, size = 0;
   Double_t *mcsAngle2 = new Double_t[2*nChambers];
   Double_t *zMCS = new Double_t[2*nChambers];
-  Int_t *indices = new Int_t[2*fNTrackHits];
+  Int_t *indices = new Int_t[2*nClusters];
   
   // Compute multiple scattering dispersion angle at each chamber
   // and save the z position where it is calculated
-  for (Int_t hitNumber = 0; hitNumber < fNTrackHits; hitNumber++) {
-    trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber);
-    hitForRec = trackParamAtHit->GetHitForRecPtr();
+  for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
+    trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster);
     
     // look for missing chambers if any
-    currentChamber = hitForRec->GetChamberNumber();
+    currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
     while (currentChamber > expectedChamber) {
       
       // Save the z position where MCS dispersion is calculated
       zMCS[size] = AliMUONConstants::DefaultChamberZ(expectedChamber);
       
-      // Do not take into account MCS in chambers prior the first hit
-      if (hitNumber > 0) {
+      // Do not take into account MCS in chambers prior the first cluster
+      if (iCluster > 0) {
         
         // Get track parameters at missing chamber z
-        extrapTrackParam = *trackParamAtHit;
-        AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, zMCS[expectedChamber]);
+        extrapTrackParam = *trackParamAtCluster;
+        AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, zMCS[size]);
         
         // Save multiple scattering dispersion angle in missing chamber
         mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
@@ -839,36 +882,36 @@ void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const
     }
     
     // Save z position where MCS dispersion is calculated
-    zMCS[size] = trackParamAtHit->GetZ();
+    zMCS[size] = trackParamAtCluster->GetZ();
     
     // Save multiple scattering dispersion angle in current chamber
-    mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(*trackParamAtHit,AliMUONConstants::ChamberThicknessInX0(),1.);
+    mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(*trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
     
     // Save indice in zMCS array corresponding to the current cluster
-    indices[hitNumber] = size;
+    indices[iCluster] = size;
     
     expectedChamber = currentChamber + 1;
     size++;
   }
   
-  // complete array of z if last hit is on the last but one chamber
+  // complete array of z if last cluster is on the last but one chamber
   if (currentChamber != nChambers-1) zMCS[size++] = AliMUONConstants::DefaultChamberZ(nChambers-1);
   
   // Compute the covariance matrix
-  for (Int_t hitNumber1 = 0; hitNumber1 < fNTrackHits; hitNumber1++) { 
+  for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) { 
     
-    for (Int_t hitNumber2 = hitNumber1; hitNumber2 < fNTrackHits; hitNumber2++) {
+    for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) {
       
       // Initialization to 0 (diagonal plus upper triangular part)
-      mcsCovariances(hitNumber1,hitNumber2) = 0.;
+      mcsCovariances(iCluster1,iCluster2) = 0.;
       
       // Compute contribution from multiple scattering in upstream chambers
-      for (Int_t k = 0; k < indices[hitNumber1]; k++) {        
-       mcsCovariances(hitNumber1,hitNumber2) += (zMCS[indices[hitNumber1]] - zMCS[k]) * (zMCS[indices[hitNumber2]] - zMCS[k]) * mcsAngle2[k];
+      for (Int_t k = 0; k < indices[iCluster1]; k++) {         
+       mcsCovariances(iCluster1,iCluster2) += (zMCS[indices[iCluster1]] - zMCS[k]) * (zMCS[indices[iCluster2]] - zMCS[k]) * mcsAngle2[k];
       }
       
       // Symetrize the matrix
-      mcsCovariances(hitNumber2,hitNumber1) = mcsCovariances(hitNumber1,hitNumber2);
+      mcsCovariances(iCluster2,iCluster1) = mcsCovariances(iCluster1,iCluster2);
     }
     
   }
@@ -880,28 +923,28 @@ void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const
 }
 
   //__________________________________________________________________________
-Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* track) const
+Int_t AliMUONTrack::ClustersInCommon(AliMUONTrack* track) const
 {
-  /// Returns the number of hits in common between the current track ("this")
+  /// Returns the number of clusters in common between the current track ("this")
   /// and the track pointed to by "track".
-  Int_t hitsInCommon = 0;
-  AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2;
-  // Loop over hits of first track
-  trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->First();
-  while (trackParamAtHit1) {
-    // Loop over hits of second track
-    trackParamAtHit2 = (AliMUONTrackParam*) track->fTrackParamAtHit->First();
-    while (trackParamAtHit2) {
-      // Increment "hitsInCommon" if both TrackParamAtHits point to the same HitForRec
-      if ((trackParamAtHit1->GetHitForRecPtr()) == (trackParamAtHit2->GetHitForRecPtr())) {
-        hitsInCommon++;
+  Int_t clustersInCommon = 0;
+  AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2;
+  // Loop over clusters of first track
+  trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->First();
+  while (trackParamAtCluster1) {
+    // Loop over clusters of second track
+    trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->First();
+    while (trackParamAtCluster2) {
+      // Increment "clustersInCommon" if both trackParamAtCluster1 & 2 point to the same cluster
+      if ((trackParamAtCluster1->GetClusterPtr()) == (trackParamAtCluster2->GetClusterPtr())) {
+        clustersInCommon++;
        break;
       }
-      trackParamAtHit2 = (AliMUONTrackParam*) track->fTrackParamAtHit->After(trackParamAtHit2);
-    } // trackParamAtHit2
-    trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->After(trackParamAtHit1);
-  } // trackParamAtHit1
-  return hitsInCommon;
+      trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->After(trackParamAtCluster2);
+    } // trackParamAtCluster2
+    trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->After(trackParamAtCluster1);
+  } // trackParamAtCluster1
+  return clustersInCommon;
 }
 
   //__________________________________________________________________________
@@ -909,68 +952,90 @@ Double_t AliMUONTrack::GetNormalizedChi2() const
 {
   /// return the chi2 value divided by the number of degrees of freedom (or 1.e10 if ndf < 0)
   
-  Double_t numberOfDegFree = (2. * fNTrackHits - 5.);
+  Double_t numberOfDegFree = (2. * GetNClusters() - 5.);
   if (numberOfDegFree > 0.) return fGlobalChi2 / numberOfDegFree;
   else return 1.e10;
 }
 
   //__________________________________________________________________________
-Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack * track, Double_t sigma2Cut) const
+Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack *track, Double_t sigmaCut) const
 {
-  /// Return kTRUE/kFALSE for each chamber if hit is compatible or not 
-  TClonesArray *hitArray, *thisHitArray;
-  AliMUONHitForRec *hit, *thisHit;
-  Int_t chamberNumber;
-  Float_t deltaZ;
-  Float_t deltaZMax = 1.; // 1 cm
-  Float_t chi2 = 0;
-  Bool_t *nCompHit = new Bool_t[AliMUONConstants::NTrackingCh()]; 
-
-  for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
-    nCompHit[ch] = kFALSE;
-  }
+  /// for each chamber: return kTRUE (kFALSE) if clusters are compatible (not compatible)
+  AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2;
+  AliMUONVCluster *cluster1, *cluster2;
+  Double_t chi2, dX, dY, dZ;
+  Double_t chi2Max = sigmaCut * sigmaCut;
+  Double_t dZMax = 1.; // 1 cm
+  
+  Bool_t *compatibleCluster = new Bool_t[AliMUONConstants::NTrackingCh()]; 
+  for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) compatibleCluster[ch] = kFALSE;
 
-  thisHitArray = this->GetHitForRecAtHit();
-
-  hitArray =  track->GetHitForRecAtHit();
-
-  for (Int_t iHthis = 0; iHthis < thisHitArray->GetEntriesFast(); iHthis++) {
-    thisHit = (AliMUONHitForRec*) thisHitArray->At(iHthis);
-    chamberNumber = thisHit->GetChamberNumber();
-    if (chamberNumber < 0 || chamberNumber > AliMUONConstants::NTrackingCh()) continue; 
-    nCompHit[chamberNumber] = kFALSE;
-    for (Int_t iH = 0; iH < hitArray->GetEntriesFast(); iH++) {
-      hit = (AliMUONHitForRec*) hitArray->At(iH);
-      deltaZ = TMath::Abs(thisHit->GetZ() - hit->GetZ());
-      chi2 = thisHit->NormalizedChi2WithHitForRec(hit,sigma2Cut); // set cut to 4 sigmas
-      if (chi2 < 3. && deltaZ < deltaZMax) {
-       nCompHit[chamberNumber] = kTRUE;
-       break;
-      }
-    }  
+  // Loop over clusters of first track
+  trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->First();
+  while (trackParamAtCluster1) {
+    
+    cluster1 = trackParamAtCluster1->GetClusterPtr();
+    
+    // Loop over clusters of second track
+    trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->First();
+    while (trackParamAtCluster2) {
+      
+      cluster2 = trackParamAtCluster2->GetClusterPtr();
+      
+      //prepare next step
+      trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->After(trackParamAtCluster2);
+      
+      // z direction
+      dZ = cluster1->GetZ() - cluster2->GetZ();
+      if (dZ > dZMax) continue;
+      
+      // non bending direction
+      dX = cluster1->GetX() - cluster2->GetX();
+      chi2 = dX * dX / (cluster1->GetErrX2() + cluster2->GetErrX2());
+      if (chi2 > chi2Max) continue;
+      
+      // bending direction
+      dY = cluster1->GetY() - cluster2->GetY();
+      chi2 = dY * dY / (cluster1->GetErrY2() + cluster2->GetErrY2());
+      if (chi2 > chi2Max) continue;
+      
+      compatibleCluster[cluster1->GetChamberId()] = kTRUE;
+      break;
+    }
+    
+    trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->After(trackParamAtCluster1);
   }
   
-  return nCompHit;
+  return compatibleCluster;
 }
 
-  //__________________________________________________________________________
-void AliMUONTrack::RecursiveDump(void) const
+//__________________________________________________________________________
+void AliMUONTrack::SetTrackParamAtVertex(const AliMUONTrackParam* trackParam)
+{
+  /// set track parameters at vertex
+  if (trackParam == 0x0) return;
+  if (fTrackParamAtVertex) *fTrackParamAtVertex = *trackParam;
+  else fTrackParamAtVertex = new AliMUONTrackParam(*trackParam);
+}
+
+//__________________________________________________________________________
+void AliMUONTrack::RecursiveDump() const
 {
-  /// Recursive dump of AliMUONTrack, i.e. with dump of TrackParamAtHit's and attached HitForRec's
-  AliMUONTrackParam *trackParamAtHit;
-  AliMUONHitForRec *hitForRec;
+  /// Recursive dump of AliMUONTrack, i.e. with dump of trackParamAtCluster and attached clusters
+  AliMUONTrackParam *trackParamAtCluster;
+  AliMUONVCluster *cluster;
   cout << "Recursive dump of Track: " << this << endl;
   // Track
   this->Dump();
-  for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) {
-    trackParamAtHit = (AliMUONTrackParam*) ((*fTrackParamAtHit)[trackHitIndex]);
-    // TrackHit
-    cout << "TrackParamAtHit: " << trackParamAtHit << " (index: " << trackHitIndex << ")" << endl;
-    trackParamAtHit->Dump();
-    hitForRec = trackParamAtHit->GetHitForRecPtr();
-    // HitForRec
-    cout << "HitForRec: " << hitForRec << endl;
-    hitForRec->Dump();
+  for (Int_t iCluster = 0; iCluster < GetNClusters(); iCluster++) {
+    trackParamAtCluster = (AliMUONTrackParam*) ((*fTrackParamAtCluster)[iCluster]);
+    // trackParamAtCluster
+    cout << "trackParamAtCluster: " << trackParamAtCluster << " (index: " << iCluster << ")" << endl;
+    trackParamAtCluster->Dump();
+    cluster = trackParamAtCluster->GetClusterPtr();
+    // cluster
+    cout << "cluster: " << cluster << endl;
+    cluster->Print();
   }
   return;
 }
@@ -980,12 +1045,12 @@ void AliMUONTrack::Print(Option_t*) const
 {
   /// Printing Track information 
 
-  cout << "<AliMUONTrack> No.Clusters=" << setw(2)   << GetNTrackHits() << 
+  cout << "<AliMUONTrack> No.Clusters=" << setw(2)   << GetNClusters() << 
       ", Match2Trig=" << setw(1) << GetMatchTrigger()  << 
       ", LoTrgNum=" << setw(3) << GetLoTrgNum()  << 
     ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) <<  GetChi2MatchTrigger();
   cout << Form(" HitTriggerPattern %x",fHitsPatternInTrigCh) << endl;
-  GetTrackParamAtHit()->First()->Print("FULL");
+  fTrackParamAtCluster->First()->Print("FULL");
 }
 
 //__________________________________________________________________________