MUON/AliMUONTrackReconstructorK.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //-----------------------------------------------------------------------------
  19 /// \class AliMUONTrackReconstructorK
  20 ///
  21 /// MUON track reconstructor using the kalman method
  22 ///
  23 /// This class contains as data:
  24 /// - the parameters for the track reconstruction
  25 ///
  26 /// It contains as methods, among others:
  27 /// - MakeTracks to build the tracks
  28 ///
  29 //-----------------------------------------------------------------------------
  30
  31 #include "AliMUONTrackReconstructorK.h"
  32
  33 #include "AliMUONConstants.h"
  34 #include "AliMUONVCluster.h"
  35 #include "AliMUONVClusterServer.h"
  36 #include "AliMUONVClusterStore.h"
  37 #include "AliMUONTrack.h"
  38 #include "AliMUONTrackParam.h"
  39 #include "AliMUONTrackExtrap.h"
  40 #include "AliMUONRecoParam.h"
  41
  42 #include "AliMpArea.h"
  43
  44 #include "AliLog.h"
  45
  46 #include <Riostream.h>
  47 #include <TMath.h>
  48 #include <TMatrixD.h>
  49 #include <TClonesArray.h>
  50
  51 /// \cond CLASSIMP
  52 ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context
  53 /// \endcond
  54
  55   //__________________________________________________________________________
  56 AliMUONTrackReconstructorK::AliMUONTrackReconstructorK(const AliMUONRecoParam* recoParam, AliMUONVClusterServer* clusterServer)
  57   : AliMUONVTrackReconstructor(recoParam, clusterServer)
  58 {
  59   /// Constructor
  60 }
  61
  62   //__________________________________________________________________________
  63 AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK()
  64 {
  65 /// Destructor
  66 }
  67
  68   //__________________________________________________________________________
  69 Bool_t AliMUONTrackReconstructorK::MakeTrackCandidates(AliMUONVClusterStore& clusterStore)
  70 {
  71   /// To make track candidates (assuming linear propagation if AliMUONRecoParam::MakeTrackCandidatesFast() return kTRUE):
  72   /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4.
  73   /// Good candidates are made of at least three clusters if both stations are requested (two otherwise).
  74   /// Keep only best candidates or all of them according to the flag AliMUONRecoParam::TrackAllTracks().
  75
  76   TClonesArray *segments;
  77   AliMUONObjectPair *segment;
  78   AliMUONTrack *track;
  79   Int_t iCandidate = 0;
  80   Bool_t clusterFound;
  81
  82   AliDebug(1,"Enter MakeTrackCandidates");
  83
  84   // Unless we're doing combined tracking, we'll clusterize all stations at once
  85   Int_t firstChamber(0);
  86   Int_t lastChamber(9);
  87
  88   if (GetRecoParam()->CombineClusterTrackReco()) {
  89     // ... Here's the exception : ask the clustering to reconstruct
  90     // clusters *only* in station 4 and 5 for combined tracking
  91     firstChamber = 6;
  92   }
  93
  94   for (Int_t i = firstChamber; i <= lastChamber; ++i )
  95   {
  96     if (fClusterServer && GetRecoParam()->UseChamber(i)) fClusterServer->Clusterize(i, clusterStore, AliMpArea(), GetRecoParam());
  97   }
  98
  99   // Loop over stations(1..) 5 and 4 and make track candidates
 100   for (Int_t istat=4; istat>=3; istat--) {
 101
 102     // Make segments in the station
 103     segments = MakeSegmentsBetweenChambers(clusterStore, 2*istat, 2*istat+1);
 104
 105     // Loop over segments
 106     for (Int_t iSegment=0; iSegment<segments->GetEntriesFast(); iSegment++) {
 107       AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
 108       segment = (AliMUONObjectPair*) segments->UncheckedAt(iSegment);
 109
 110       // Transform segments to tracks and put them at the end of fRecTracksPtr
 111       track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(segment,GetRecoParam()->GetBendingVertexDispersion());
 112       fNRecTracks++;
 113
 114       // Look for compatible cluster(s) in the other station
 115       if (GetRecoParam()->MakeTrackCandidatesFast()) clusterFound = FollowLinearTrackInStation(*track, clusterStore, 7-istat);
 116       else clusterFound = FollowTrackInStation(*track, clusterStore, 7-istat);
 117
 118       // Remove track if no cluster found on a requested station
 119       // or abort tracking if there are too many candidates
 120       if (GetRecoParam()->RequestStation(7-istat)) {
 121         if (!clusterFound) {
 122           fRecTracksPtr->Remove(track);
 123           fNRecTracks--;
 124         } else if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
 125           AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks));
 126           delete segments;
 127           return kFALSE;
 128         }
 129       } else {
 130         if ((fNRecTracks + segments->GetEntriesFast() - iSegment - 1) > GetRecoParam()->GetMaxTrackCandidates()) {
 131           AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks + segments->GetEntriesFast() - iSegment - 1));
 132           delete segments;
 133           return kFALSE;
 134         }
 135       }
 136
 137     }
 138
 139     // delete the array of segments
 140     delete segments;
 141   }
 142
 143   // Keep all different tracks if required
 144   if (GetRecoParam()->TrackAllTracks()) RemoveIdenticalTracks();
 145
 146   // Retrace tracks using Kalman filter and select them if needed
 147   Int_t nCurrentTracks = fRecTracksPtr->GetLast()+1;
 148   for (Int_t iRecTrack = 0; iRecTrack < nCurrentTracks; iRecTrack++) {
 149     track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
 150
 151     // skip empty slots
 152     if(!track) continue;
 153
 154     // retrace tracks using Kalman filter and remove the ones for which extrap failed or that are out of limits
 155     if (!RetraceTrack(*track,kTRUE) || !IsAcceptable(*((AliMUONTrackParam*)track->GetTrackParamAtCluster()->First()))) {
 156       fRecTracksPtr->Remove(track);
 157       fNRecTracks--;
 158     }
 159
 160   }
 161
 162   // Keep only the best tracks if required
 163   if (!GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
 164   else fRecTracksPtr->Compress();
 165
 166   AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
 167
 168   return kTRUE;
 169
 170 }
 171
 172   //__________________________________________________________________________
 173 Bool_t AliMUONTrackReconstructorK::MakeMoreTrackCandidates(AliMUONVClusterStore& clusterStore)
 174 {
 175   /// To make extra track candidates assuming linear propagation:
 176   /// clustering is supposed to be already done
 177   /// Start with segments made of 1 cluster in each of the stations 4 and 5 then follow track in remaining chambers.
 178   /// Good candidates are made of at least three clusters if both stations are requested (two otherwise).
 179   /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
 180
 181   TClonesArray *segments;
 182   AliMUONObjectPair *segment;
 183   AliMUONTrack *track;
 184   Int_t iCandidate = 0, iCurrentTrack, nCurrentTracks;
 185   Int_t initialNRecTracks = fNRecTracks;
 186   Bool_t clusterFound;
 187
 188   AliDebug(1,"Enter MakeMoreTrackCandidates");
 189
 190   // Double loop over chambers in stations(1..) 4 and 5 to make track candidates
 191   for (Int_t ich1 = 6; ich1 <= 7; ich1++) {
 192     for (Int_t ich2 = 8; ich2 <= 9; ich2++) {
 193
 194       // Make segments between ch1 and ch2
 195       segments = MakeSegmentsBetweenChambers(clusterStore, ich1, ich2);
 196
 197       /// Remove segments already attached to a track
 198       RemoveUsedSegments(*segments);
 199
 200       // Loop over segments
 201       for (Int_t iSegment=0; iSegment<segments->GetEntriesFast(); iSegment++) {
 202         AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
 203         segment = (AliMUONObjectPair*) segments->UncheckedAt(iSegment);
 204
 205         // Transform segments to tracks and put them at the end of fRecTracksPtr
 206         iCurrentTrack = fRecTracksPtr->GetLast()+1;
 207         track = new ((*fRecTracksPtr)[iCurrentTrack]) AliMUONTrack(segment,GetRecoParam()->GetBendingVertexDispersion());
 208         fNRecTracks++;
 209
 210         // Look for compatible cluster(s) in the second chamber of station 5
 211         clusterFound = FollowLinearTrackInChamber(*track, clusterStore, 17-ich2);
 212
 213         // skip the original track in case it has been removed
 214         if (GetRecoParam()->TrackAllTracks() && clusterFound) iCurrentTrack++;
 215
 216         // loop over every new tracks
 217         nCurrentTracks = fRecTracksPtr->GetLast()+1;
 218         while (iCurrentTrack < nCurrentTracks) {
 219           track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iCurrentTrack);
 220
 221           // Look for compatible cluster(s) in the second chamber of station 4
 222           FollowLinearTrackInChamber(*track, clusterStore, 13-ich1);
 223
 224           iCurrentTrack++;
 225         }
 226
 227         // abort tracking if there are too many candidates
 228         if ((fNRecTracks + segments->GetEntriesFast() - iSegment - 1) > GetRecoParam()->GetMaxTrackCandidates()) {
 229           AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks + segments->GetEntriesFast() - iSegment - 1));
 230           delete segments;
 231           return kFALSE;
 232         }
 233
 234       }
 235
 236       // delete the array of segments
 237       delete segments;
 238     }
 239   }
 240
 241   // Retrace tracks using Kalman filter (also compute track chi2) and select them
 242   nCurrentTracks = fRecTracksPtr->GetLast()+1;
 243   for (Int_t iRecTrack = initialNRecTracks; iRecTrack < nCurrentTracks; iRecTrack++) {
 244     track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
 245
 246     // skip empty slots
 247     if(!track) continue;
 248
 249     // retrace tracks using Kalman filter and remove the ones for which extrap failed or that are out of limits
 250     if (!RetraceTrack(*track,kTRUE) || !IsAcceptable(*((AliMUONTrackParam*)track->GetTrackParamAtCluster()->First()))) {
 251       fRecTracksPtr->Remove(track);
 252       fNRecTracks--;
 253     }
 254
 255   }
 256
 257   // Keep only the best tracks if required
 258   if (!GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
 259   else fRecTracksPtr->Compress();
 260
 261   AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
 262
 263   return kTRUE;
 264
 265 }
 266
 267   //__________________________________________________________________________
 268 Bool_t AliMUONTrackReconstructorK::RetraceTrack(AliMUONTrack &trackCandidate, Bool_t resetSeed)
 269 {
 270   /// Re-run the kalman filter from the most downstream cluster to the most uptream one
 271   /// Return kFALSE in case of failure (i.e. extrapolation problem)
 272   AliDebug(1,"Enter RetraceTrack");
 273
 274   AliMUONTrackParam* lastTrackParam = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtCluster()->Last();
 275
 276   // Reset the "seed" (= track parameters and their covariances at last cluster) if required
 277   if (resetSeed) {
 278
 279     // parameters at last cluster
 280     AliMUONVCluster* cluster2 = lastTrackParam->GetClusterPtr();
 281     Double_t x2 = cluster2->GetX();
 282     Double_t y2 = cluster2->GetY();
 283     Double_t z2 = cluster2->GetZ();
 284
 285     // parameters at last but one cluster
 286     AliMUONTrackParam* previousTrackParam = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtCluster()->Before(lastTrackParam);
 287     AliMUONVCluster* cluster1 = previousTrackParam->GetClusterPtr();
 288     // make sure it is on the previous chamber (can have 2 clusters in the same chamber after "ComplementTrack")
 289     if (cluster2->GetChamberId() == cluster1->GetChamberId()) {
 290       previousTrackParam = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtCluster()->Before(previousTrackParam);
 291       cluster1 = previousTrackParam->GetClusterPtr();
 292     }
 293     Double_t x1 = cluster1->GetX();
 294     Double_t y1 = cluster1->GetY();
 295     Double_t z1 = cluster1->GetZ();
 296
 297     // reset track parameters
 298     Double_t dZ = z1 - z2;
 299     lastTrackParam->SetNonBendingCoor(x2);
 300     lastTrackParam->SetBendingCoor(y2);
 301     lastTrackParam->SetZ(z2);
 302     lastTrackParam->SetNonBendingSlope((x1 - x2) / dZ);
 303     lastTrackParam->SetBendingSlope((y1 - y2) / dZ);
 304     Double_t bendingImpact = y2 - z2 * lastTrackParam->GetBendingSlope();
 305     Double_t inverseBendingMomentum = 1. / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact);
 306     lastTrackParam->SetInverseBendingMomentum(inverseBendingMomentum);
 307
 308     // => Reset track parameter covariances at last cluster (as if the other clusters did not exist)
 309     TMatrixD lastParamCov(5,5);
 310     lastParamCov.Zero();
 311     // Non bending plane
 312     lastParamCov(0,0) = cluster2->GetErrX2();
 313     lastParamCov(0,1) = - cluster2->GetErrX2() / dZ;
 314     lastParamCov(1,0) = lastParamCov(0,1);
 315     lastParamCov(1,1) = ( 1000. * cluster1->GetErrX2() + cluster2->GetErrX2() ) / dZ / dZ;
 316     // Bending plane
 317     lastParamCov(2,2) = cluster2->GetErrY2();
 318     lastParamCov(2,3) = - cluster2->GetErrY2() / dZ;
 319     lastParamCov(3,2) = lastParamCov(2,3);
 320     lastParamCov(3,3) = ( 1000. * cluster1->GetErrY2() + cluster2->GetErrY2() ) / dZ / dZ;
 321     // Inverse bending momentum (vertex resolution + bending slope resolution + 10% error on dipole parameters+field)
 322     if (AliMUONTrackExtrap::IsFieldON()) {
 323       lastParamCov(4,4) = ((GetRecoParam()->GetBendingVertexDispersion() *
 324                             GetRecoParam()->GetBendingVertexDispersion() +
 325                             (z1 * z1 * cluster2->GetErrY2() + z2 * z2 * 1000. * cluster1->GetErrY2()) / dZ / dZ) /
 326                            bendingImpact / bendingImpact + 0.1 * 0.1) * inverseBendingMomentum * inverseBendingMomentum;
 327       lastParamCov(2,4) = z1 * cluster2->GetErrY2() * inverseBendingMomentum / bendingImpact / dZ;
 328       lastParamCov(4,2) = lastParamCov(2,4);
 329       lastParamCov(3,4) = - (z1 * cluster2->GetErrY2() + z2 * 1000. * cluster1->GetErrY2()) *
 330       inverseBendingMomentum / bendingImpact / dZ / dZ;
 331       lastParamCov(4,3) = lastParamCov(3,4);
 332     } else lastParamCov(4,4) = inverseBendingMomentum*inverseBendingMomentum;
 333     lastTrackParam->SetCovariances(lastParamCov);
 334
 335     // Reset the track chi2
 336     lastTrackParam->SetTrackChi2(0.);
 337
 338   }
 339
 340   // Redo the tracking
 341   return RetracePartialTrack(trackCandidate, lastTrackParam);
 342
 343 }
 344
 345   //__________________________________________________________________________
 346 Bool_t AliMUONTrackReconstructorK::RetracePartialTrack(AliMUONTrack &trackCandidate, const AliMUONTrackParam* startingTrackParam)
 347 {
 348   /// Re-run the kalman filter from the cluster attached to startingTrackParam to the most uptream cluster
 349   /// Return kFALSE in case of failure (i.e. extrapolation problem)
 350   AliDebug(1,"Enter RetracePartialTrack");
 351
 352   // Printout for debuging
 353   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 354     cout << "RetracePartialTrack: track chi2 before re-tracking: " << trackCandidate.GetGlobalChi2() << endl;
 355   }
 356
 357   // Reset the track chi2
 358   trackCandidate.SetGlobalChi2(startingTrackParam->GetTrackChi2());
 359
 360   // loop over attached clusters until the first one and recompute track parameters and covariances using kalman filter
 361   Bool_t extrapStatus = kTRUE;
 362   Int_t expectedChamber = startingTrackParam->GetClusterPtr()->GetChamberId() - 1;
 363   Int_t currentChamber;
 364   Double_t addChi2TrackAtCluster;
 365   AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtCluster()->Before(startingTrackParam);
 366   while (trackParamAtCluster) {
 367
 368     // reset track parameters and their covariances
 369     trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
 370     trackParamAtCluster->SetZ(startingTrackParam->GetZ());
 371     trackParamAtCluster->SetCovariances(startingTrackParam->GetCovariances());
 372
 373     // add MCS effect
 374     AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(expectedChamber+1),-1.);
 375
 376     // reset propagator for smoother
 377     if (GetRecoParam()->UseSmoother()) trackParamAtCluster->ResetPropagator();
 378
 379     // add MCS in missing chambers if any
 380     currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
 381     while (currentChamber < expectedChamber) {
 382       // extrapolation to the missing chamber (update the propagator)
 383       if (!AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, AliMUONConstants::DefaultChamberZ(expectedChamber),
 384                                             GetRecoParam()->UseSmoother())) extrapStatus = kFALSE;
 385       // add MCS effect
 386       AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(expectedChamber),-1.);
 387       expectedChamber--;
 388     }
 389
 390     // extrapolation to the plane of the cluster attached to the current trackParamAtCluster (update the propagator)
 391     if (!AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ(),
 392                                           GetRecoParam()->UseSmoother())) extrapStatus = kFALSE;
 393
 394     if (GetRecoParam()->UseSmoother()) {
 395       // save extrapolated parameters for smoother
 396       trackParamAtCluster->SetExtrapParameters(trackParamAtCluster->GetParameters());
 397
 398       // save extrapolated covariance matrix for smoother
 399       trackParamAtCluster->SetExtrapCovariances(trackParamAtCluster->GetCovariances());
 400     }
 401
 402     // Compute new track parameters using kalman filter
 403     addChi2TrackAtCluster = RunKalmanFilter(*trackParamAtCluster);
 404
 405     // Update the track chi2
 406     trackCandidate.SetGlobalChi2(trackCandidate.GetGlobalChi2() + addChi2TrackAtCluster);
 407     trackParamAtCluster->SetTrackChi2(trackCandidate.GetGlobalChi2());
 408
 409     // prepare next step
 410     expectedChamber = currentChamber - 1;
 411     startingTrackParam = trackParamAtCluster;
 412     trackParamAtCluster = (AliMUONTrackParam*) (trackCandidate.GetTrackParamAtCluster()->Before(startingTrackParam));
 413   }
 414
 415   // Printout for debuging
 416   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 417     cout << "RetracePartialTrack: track chi2 after re-tracking: " << trackCandidate.GetGlobalChi2() << endl;
 418   }
 419
 420   // set global chi2 to max value in case of problem during track extrapolation
 421   if (!extrapStatus) trackCandidate.SetGlobalChi2(2.*AliMUONTrack::MaxChi2());
 422   return extrapStatus;
 423
 424 }
 425
 426   //__________________________________________________________________________
 427 Bool_t AliMUONTrackReconstructorK::FollowTracks(AliMUONVClusterStore& clusterStore)
 428 {
 429   /// Follow tracks in stations(1..) 3, 2 and 1
 430   AliDebug(1,"Enter FollowTracks");
 431
 432   AliMUONTrack *track;
 433   Int_t currentNRecTracks;
 434
 435   for (Int_t station = 2; station >= 0; station--) {
 436
 437     // Save the actual number of reconstructed track in case of
 438     // tracks are added or suppressed during the tracking procedure
 439     // !! Do not compress fRecTracksPtr until the end of the loop over tracks !!
 440     currentNRecTracks = fNRecTracks;
 441
 442     for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) {
 443       AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1));
 444
 445       track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
 446
 447       // Look for compatible cluster(s) in station(0..) "station"
 448       if (!FollowTrackInStation(*track, clusterStore, station)) {
 449
 450         // Try to recover track if required
 451         if (GetRecoParam()->RecoverTracks()) {
 452
 453           // work on a copy of the track if this station is not required
 454           // to keep the case where no cluster is reconstructed as a possible candidate
 455           if (!GetRecoParam()->RequestStation(station)) {
 456             track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*track);
 457             fNRecTracks++;
 458           }
 459
 460           // try to recover
 461           if (!RecoverTrack(*track, clusterStore, station)) {
 462             // remove track if no cluster found
 463             fRecTracksPtr->Remove(track);
 464             fNRecTracks--;
 465           }
 466
 467         } else if (GetRecoParam()->RequestStation(station)) {
 468           // remove track if no cluster found
 469           fRecTracksPtr->Remove(track);
 470           fNRecTracks--;
 471         }
 472
 473       }
 474
 475       // abort tracking if there are too many candidates
 476       if (GetRecoParam()->RequestStation(station)) {
 477         if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
 478           AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks));
 479           return kFALSE;
 480         }
 481       } else {
 482         if ((fNRecTracks + currentNRecTracks - iRecTrack - 1) > GetRecoParam()->GetMaxTrackCandidates()) {
 483           AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks + currentNRecTracks - iRecTrack - 1));
 484           return kFALSE;
 485         }
 486       }
 487
 488     }
 489
 490     fRecTracksPtr->Compress(); // this is essential before checking tracks
 491
 492     // Keep only the best tracks if required
 493     if (!GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
 494
 495   }
 496
 497   return kTRUE;
 498
 499 }
 500
 501   //__________________________________________________________________________
 502 Bool_t AliMUONTrackReconstructorK::FollowTrackInChamber(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextChamber)
 503 {
 504   /// Follow trackCandidate in chamber(0..) nextChamber and search for compatible cluster(s)
 505   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
 506   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
 507   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
 508   ///         Remove the obsolete "trackCandidate" at the end.
 509   /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
 510   /// return kTRUE if new cluster(s) have been found (otherwise return kFALSE)
 511   AliDebug(1,Form("Enter FollowTrackInChamber(1..) %d", nextChamber+1));
 512
 513   Double_t chi2OfCluster;
 514   Double_t maxChi2OfCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
 515                                    GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
 516   Double_t addChi2TrackAtCluster;
 517   Double_t bestAddChi2TrackAtCluster = AliMUONTrack::MaxChi2();
 518   Bool_t foundOneCluster = kFALSE;
 519   AliMUONTrack *newTrack = 0x0;
 520   AliMUONVCluster *cluster;
 521   AliMUONTrackParam extrapTrackParamAtCh;
 522   AliMUONTrackParam extrapTrackParamAtCluster;
 523   AliMUONTrackParam bestTrackParamAtCluster;
 524
 525   // Get track parameters according to the propagation direction
 526   if (nextChamber > 7) extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last();
 527   else extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First();
 528
 529   // Printout for debuging
 530   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 531     cout<<endl<<"Track parameters and covariances at first cluster:"<<endl;
 532     extrapTrackParamAtCh.GetParameters().Print();
 533     extrapTrackParamAtCh.GetCovariances().Print();
 534   }
 535
 536   // Add MCS effect
 537   Int_t currentChamber = extrapTrackParamAtCh.GetClusterPtr()->GetChamberId();
 538   AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.);
 539
 540   // reset propagator for smoother
 541   if (GetRecoParam()->UseSmoother()) extrapTrackParamAtCh.ResetPropagator();
 542
 543   // Add MCS in the missing chamber(s) if any
 544   while (currentChamber > nextChamber + 1) {
 545     // extrapolation to the missing chamber
 546     currentChamber--;
 547     if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(currentChamber),
 548                                           GetRecoParam()->UseSmoother())) return kFALSE;
 549     // add MCS effect
 550     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.);
 551   }
 552
 553   //Extrapolate trackCandidate to chamber
 554   if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(nextChamber),
 555                                         GetRecoParam()->UseSmoother())) return kFALSE;
 556
 557   // Printout for debuging
 558   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 559     cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(nextChamber)<<":"<<endl;
 560     extrapTrackParamAtCh.GetParameters().Print();
 561     extrapTrackParamAtCh.GetCovariances().Print();
 562   }
 563
 564   // Printout for debuging
 565   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 566     cout << "FollowTrackInChamber: look for clusters in chamber(1..): " << nextChamber+1 << endl;
 567   }
 568
 569   // Ask the clustering to reconstruct new clusters around the track position in the current chamber
 570   // except for station 4 and 5 that are already entirely clusterized
 571   if (GetRecoParam()->CombineClusterTrackReco()) {
 572     if (nextChamber < 6) AskForNewClustersInChamber(extrapTrackParamAtCh, clusterStore, nextChamber);
 573   }
 574
 575   // Create iterators to loop over clusters in both chambers
 576   TIter next(clusterStore.CreateChamberIterator(nextChamber,nextChamber));
 577
 578   // look for candidates in chamber
 579   while ( ( cluster = static_cast<AliMUONVCluster*>(next()) ) ) {
 580
 581     // try to add the current cluster fast
 582     if (!TryOneClusterFast(extrapTrackParamAtCh, cluster)) continue;
 583
 584     // try to add the current cluster accuratly
 585     chi2OfCluster = TryOneCluster(extrapTrackParamAtCh, cluster, extrapTrackParamAtCluster,
 586                                   GetRecoParam()->UseSmoother());
 587
 588     // if good chi2 then consider to add cluster
 589     if (chi2OfCluster < maxChi2OfCluster) {
 590
 591       if (GetRecoParam()->UseSmoother()) {
 592         // save extrapolated parameters for smoother
 593         extrapTrackParamAtCluster.SetExtrapParameters(extrapTrackParamAtCluster.GetParameters());
 594
 595         // save extrapolated covariance matrix for smoother
 596         extrapTrackParamAtCluster.SetExtrapCovariances(extrapTrackParamAtCluster.GetCovariances());
 597       }
 598
 599       // Compute new track parameters including new cluster using kalman filter
 600       addChi2TrackAtCluster = RunKalmanFilter(extrapTrackParamAtCluster);
 601
 602       // skip track out of limits
 603       if (!IsAcceptable(extrapTrackParamAtCluster)) continue;
 604
 605       // remember a cluster was found
 606       foundOneCluster = kTRUE;
 607
 608       // Printout for debuging
 609       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 610         cout << "FollowTrackInChamber: found one cluster in chamber(1..): " << nextChamber+1
 611         << " (Chi2 = " << chi2OfCluster << ")" << endl;
 612         cluster->Print();
 613       }
 614
 615       if (GetRecoParam()->TrackAllTracks()) {
 616         // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
 617         newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 618         UpdateTrack(*newTrack,extrapTrackParamAtCluster,addChi2TrackAtCluster);
 619         fNRecTracks++;
 620
 621         // Printout for debuging
 622         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 623           cout << "FollowTrackInChamber: added one cluster in chamber(1..): " << nextChamber+1 << endl;
 624           if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 625         }
 626
 627       } else if (addChi2TrackAtCluster < bestAddChi2TrackAtCluster) {
 628         // keep track of the best cluster
 629         bestAddChi2TrackAtCluster = addChi2TrackAtCluster;
 630         bestTrackParamAtCluster = extrapTrackParamAtCluster;
 631       }
 632
 633     }
 634
 635   }
 636
 637   // fill out the best track if required else clean up the fRecTracksPtr array
 638   if (!GetRecoParam()->TrackAllTracks()) {
 639     if (foundOneCluster) {
 640       UpdateTrack(trackCandidate,bestTrackParamAtCluster,bestAddChi2TrackAtCluster);
 641
 642       // Printout for debuging
 643       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 644         cout << "FollowTrackInChamber: added the best cluster in chamber(1..): " << bestTrackParamAtCluster.GetClusterPtr()->GetChamberId()+1 << endl;
 645         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 646       }
 647
 648     } else return kFALSE;
 649
 650   } else if (foundOneCluster) {
 651
 652     // remove obsolete track
 653     fRecTracksPtr->Remove(&trackCandidate);
 654     fNRecTracks--;
 655
 656   } else return kFALSE;
 657
 658   return kTRUE;
 659
 660 }
 661
 662   //__________________________________________________________________________
 663 Bool_t AliMUONTrackReconstructorK::FollowTrackInStation(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation)
 664 {
 665   /// Follow trackCandidate in station(0..) nextStation and search for compatible cluster(s)
 666   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
 667   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
 668   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
 669   ///         Remove the obsolete "trackCandidate" at the end.
 670   /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
 671   /// return kTRUE if new cluster(s) have been found (otherwise return kFALSE)
 672   AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1));
 673
 674   // Order the chamber according to the propagation direction (tracking starts with chamber 2):
 675   // - nextStation == station(1...) 5 => forward propagation
 676   // - nextStation < station(1...) 5 => backward propagation
 677   Int_t ch1, ch2;
 678   if (nextStation==4) {
 679     ch1 = 2*nextStation+1;
 680     ch2 = 2*nextStation;
 681   } else {
 682     ch1 = 2*nextStation;
 683     ch2 = 2*nextStation+1;
 684   }
 685
 686   Double_t chi2OfCluster;
 687   Double_t maxChi2OfCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
 688                                    GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
 689   Double_t addChi2TrackAtCluster1;
 690   Double_t addChi2TrackAtCluster2;
 691   Double_t bestAddChi2TrackAtCluster1 = AliMUONTrack::MaxChi2();
 692   Double_t bestAddChi2TrackAtCluster2 = AliMUONTrack::MaxChi2();
 693   Bool_t foundOneCluster = kFALSE;
 694   Bool_t foundTwoClusters = kFALSE;
 695   AliMUONTrack *newTrack = 0x0;
 696   AliMUONVCluster *clusterCh1, *clusterCh2;
 697   AliMUONTrackParam extrapTrackParam;
 698   AliMUONTrackParam extrapTrackParamAtCh;
 699   AliMUONTrackParam extrapTrackParamAtCluster1;
 700   AliMUONTrackParam extrapTrackParamAtCluster2;
 701   AliMUONTrackParam bestTrackParamAtCluster1;
 702   AliMUONTrackParam bestTrackParamAtCluster2;
 703
 704   // Get track parameters according to the propagation direction
 705   if (nextStation==4) extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last();
 706   else extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First();
 707
 708   // Printout for debuging
 709   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 710     cout<<endl<<"Track parameters and covariances at first cluster:"<<endl;
 711     extrapTrackParamAtCh.GetParameters().Print();
 712     extrapTrackParamAtCh.GetCovariances().Print();
 713   }
 714
 715   // Add MCS effect
 716   Int_t currentChamber = extrapTrackParamAtCh.GetClusterPtr()->GetChamberId();
 717   AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.);
 718
 719   // reset propagator for smoother
 720   if (GetRecoParam()->UseSmoother()) extrapTrackParamAtCh.ResetPropagator();
 721
 722   // Add MCS in the missing chamber(s) if any
 723   while (ch1 < ch2 && currentChamber > ch2 + 1) {
 724     // extrapolation to the missing chamber
 725     currentChamber--;
 726     if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(currentChamber),
 727                                           GetRecoParam()->UseSmoother())) return kFALSE;
 728     // add MCS effect
 729     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.);
 730   }
 731
 732   //Extrapolate trackCandidate to chamber "ch2"
 733   if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2),
 734                                         GetRecoParam()->UseSmoother())) return kFALSE;
 735
 736   // Printout for debuging
 737   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 738     cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch2)<<":"<<endl;
 739     extrapTrackParamAtCh.GetParameters().Print();
 740     extrapTrackParamAtCh.GetCovariances().Print();
 741   }
 742
 743   // Printout for debuging
 744   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 745     cout << "FollowTrackInStation: look for clusters in chamber(1..): " << ch2+1 << endl;
 746   }
 747
 748   // Ask the clustering to reconstruct new clusters around the track position in the current station
 749   // except for station 4 and 5 that are already entirely clusterized
 750   if (GetRecoParam()->CombineClusterTrackReco()) {
 751     if (nextStation < 3) AskForNewClustersInStation(extrapTrackParamAtCh, clusterStore, nextStation);
 752   }
 753
 754   Int_t nClusters = clusterStore.GetSize();
 755   Bool_t *clusterCh1Used = new Bool_t[nClusters];
 756   for (Int_t i = 0; i < nClusters; i++) clusterCh1Used[i] = kFALSE;
 757   Int_t iCluster1;
 758
 759   // Create iterators to loop over clusters in both chambers
 760   TIter nextInCh1(clusterStore.CreateChamberIterator(ch1,ch1));
 761   TIter nextInCh2(clusterStore.CreateChamberIterator(ch2,ch2));
 762
 763   // look for candidates in chamber 2
 764   while ( ( clusterCh2 = static_cast<AliMUONVCluster*>(nextInCh2()) ) ) {
 765
 766     // try to add the current cluster fast
 767     if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh2)) continue;
 768
 769     // try to add the current cluster accuratly
 770     chi2OfCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh2, extrapTrackParamAtCluster2,
 771                                   GetRecoParam()->UseSmoother());
 772
 773     // if good chi2 then try to attach a cluster in the other chamber too
 774     if (chi2OfCluster < maxChi2OfCluster) {
 775
 776       if (GetRecoParam()->UseSmoother()) {
 777         // save extrapolated parameters for smoother
 778         extrapTrackParamAtCluster2.SetExtrapParameters(extrapTrackParamAtCluster2.GetParameters());
 779
 780         // save extrapolated covariance matrix for smoother
 781         extrapTrackParamAtCluster2.SetExtrapCovariances(extrapTrackParamAtCluster2.GetCovariances());
 782       }
 783
 784       // Compute new track parameters including "clusterCh2" using kalman filter
 785       addChi2TrackAtCluster2 = RunKalmanFilter(extrapTrackParamAtCluster2);
 786
 787       // skip track out of limits
 788       if (!IsAcceptable(extrapTrackParamAtCluster2)) continue;
 789
 790       // Printout for debuging
 791       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 792         cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch2+1
 793         << " (Chi2 = " << chi2OfCluster << ")" << endl;
 794         clusterCh2->Print();
 795         cout << "                      look for second clusters in chamber(1..): " << ch1+1 << " ..." << endl;
 796       }
 797
 798       // copy new track parameters for next step
 799       extrapTrackParam = extrapTrackParamAtCluster2;
 800
 801       // add MCS effect
 802       AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(ch2),-1.);
 803
 804       // reset propagator for smoother
 805       if (GetRecoParam()->UseSmoother()) extrapTrackParam.ResetPropagator();
 806
 807       //Extrapolate track parameters to chamber "ch1"
 808       Bool_t normalExtrap = AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(ch1),
 809                                                              GetRecoParam()->UseSmoother());
 810
 811       // reset cluster iterator of chamber 1
 812       nextInCh1.Reset();
 813       iCluster1 = -1;
 814
 815       // look for second candidates in chamber 1
 816       Bool_t foundSecondCluster = kFALSE;
 817       if (normalExtrap) while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
 818         iCluster1++;
 819
 820         // try to add the current cluster fast
 821         if (!TryOneClusterFast(extrapTrackParam, clusterCh1)) continue;
 822
 823         // try to add the current cluster accuratly
 824         chi2OfCluster = TryOneCluster(extrapTrackParam, clusterCh1, extrapTrackParamAtCluster1,
 825                                       GetRecoParam()->UseSmoother());
 826
 827         // if good chi2 then consider to add the 2 clusters to the "trackCandidate"
 828         if (chi2OfCluster < maxChi2OfCluster) {
 829
 830           if (GetRecoParam()->UseSmoother()) {
 831             // save extrapolated parameters for smoother
 832             extrapTrackParamAtCluster1.SetExtrapParameters(extrapTrackParamAtCluster1.GetParameters());
 833
 834             // save extrapolated covariance matrix for smoother
 835             extrapTrackParamAtCluster1.SetExtrapCovariances(extrapTrackParamAtCluster1.GetCovariances());
 836           }
 837
 838           // Compute new track parameters including "clusterCh1" using kalman filter
 839           addChi2TrackAtCluster1 = RunKalmanFilter(extrapTrackParamAtCluster1);
 840
 841           // skip track out of limits
 842           if (!IsAcceptable(extrapTrackParamAtCluster1)) continue;
 843
 844           // remember a second cluster was found
 845           foundSecondCluster = kTRUE;
 846           foundTwoClusters = kTRUE;
 847
 848           // Printout for debuging
 849           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 850             cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch1+1
 851             << " (Chi2 = " << chi2OfCluster << ")" << endl;
 852             clusterCh1->Print();
 853           }
 854
 855           if (GetRecoParam()->TrackAllTracks()) {
 856             // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new clusters
 857             newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 858             UpdateTrack(*newTrack,extrapTrackParamAtCluster1,extrapTrackParamAtCluster2,addChi2TrackAtCluster1,addChi2TrackAtCluster2);
 859             fNRecTracks++;
 860
 861             // Tag clusterCh1 as used
 862             clusterCh1Used[iCluster1] = kTRUE;
 863
 864             // Printout for debuging
 865             if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 866               cout << "FollowTrackInStation: added two clusters in station(1..): " << nextStation+1 << endl;
 867               if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 868             }
 869
 870           } else if (addChi2TrackAtCluster1+addChi2TrackAtCluster2 < bestAddChi2TrackAtCluster1+bestAddChi2TrackAtCluster2) {
 871             // keep track of the best couple of clusters
 872             bestAddChi2TrackAtCluster1 = addChi2TrackAtCluster1;
 873             bestAddChi2TrackAtCluster2 = addChi2TrackAtCluster2;
 874             bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
 875             bestTrackParamAtCluster2 = extrapTrackParamAtCluster2;
 876           }
 877
 878         }
 879
 880       }
 881
 882       // if no clusterCh1 found then consider to add clusterCh2 only
 883       if (!foundSecondCluster) {
 884
 885         // remember a cluster was found
 886         foundOneCluster = kTRUE;
 887
 888         if (GetRecoParam()->TrackAllTracks()) {
 889           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
 890           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 891           UpdateTrack(*newTrack,extrapTrackParamAtCluster2,addChi2TrackAtCluster2);
 892           fNRecTracks++;
 893
 894           // Printout for debuging
 895           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 896             cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch2+1 << endl;
 897             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 898           }
 899
 900         } else if (!foundTwoClusters && addChi2TrackAtCluster2 < bestAddChi2TrackAtCluster1) {
 901           // keep track of the best single cluster except if a couple of clusters has already been found
 902           bestAddChi2TrackAtCluster1 = addChi2TrackAtCluster2;
 903           bestTrackParamAtCluster1 = extrapTrackParamAtCluster2;
 904         }
 905
 906       }
 907
 908     }
 909
 910   }
 911
 912   // look for candidates in chamber 1 not already attached to a track
 913   // if we want to keep all possible tracks or if no good couple of clusters has been found
 914   if (GetRecoParam()->TrackAllTracks() || !foundTwoClusters) {
 915
 916     // add MCS effect for next step
 917     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(ch2),-1.);
 918
 919     //Extrapolate trackCandidate to chamber "ch1"
 920     Bool_t normalExtrap = AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch1),
 921                                                            GetRecoParam()->UseSmoother());
 922
 923     // Printout for debuging
 924     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 925       cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch1)<<":"<<endl;
 926       extrapTrackParamAtCh.GetParameters().Print();
 927       extrapTrackParamAtCh.GetCovariances().Print();
 928     }
 929
 930     // Printout for debuging
 931     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 932       cout << "FollowTrackInStation: look for single clusters in chamber(1..): " << ch1+1 << endl;
 933     }
 934
 935     // reset cluster iterator of chamber 1
 936     nextInCh1.Reset();
 937     iCluster1 = -1;
 938
 939     // look for second candidates in chamber 1
 940     if (normalExtrap) while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
 941       iCluster1++;
 942
 943       if (clusterCh1Used[iCluster1]) continue; // Skip clusters already used
 944
 945       // try to add the current cluster fast
 946       if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh1)) continue;
 947
 948       // try to add the current cluster accuratly
 949       chi2OfCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh1, extrapTrackParamAtCluster1,
 950                                     GetRecoParam()->UseSmoother());
 951
 952       // if good chi2 then consider to add clusterCh1
 953       // We do not try to attach a cluster in the other chamber too since it has already been done above
 954       if (chi2OfCluster < maxChi2OfCluster) {
 955
 956         if (GetRecoParam()->UseSmoother()) {
 957           // save extrapolated parameters for smoother
 958           extrapTrackParamAtCluster1.SetExtrapParameters(extrapTrackParamAtCluster1.GetParameters());
 959
 960           // save extrapolated covariance matrix for smoother
 961           extrapTrackParamAtCluster1.SetExtrapCovariances(extrapTrackParamAtCluster1.GetCovariances());
 962         }
 963
 964         // Compute new track parameters including "clusterCh1" using kalman filter
 965         addChi2TrackAtCluster1 = RunKalmanFilter(extrapTrackParamAtCluster1);
 966
 967         // skip track out of limits
 968         if (!IsAcceptable(extrapTrackParamAtCluster1)) continue;
 969
 970         // remember a cluster was found
 971         foundOneCluster = kTRUE;
 972
 973         // Printout for debuging
 974         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 975           cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch1+1
 976           << " (Chi2 = " << chi2OfCluster << ")" << endl;
 977           clusterCh1->Print();
 978         }
 979
 980         if (GetRecoParam()->TrackAllTracks()) {
 981           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
 982           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 983           UpdateTrack(*newTrack,extrapTrackParamAtCluster1,addChi2TrackAtCluster1);
 984           fNRecTracks++;
 985
 986           // Printout for debuging
 987           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 988             cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch1+1 << endl;
 989             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 990           }
 991
 992         } else if (addChi2TrackAtCluster1 < bestAddChi2TrackAtCluster1) {
 993           // keep track of the best single cluster except if a couple of clusters has already been found
 994           bestAddChi2TrackAtCluster1 = addChi2TrackAtCluster1;
 995           bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
 996         }
 997
 998       }
 999
1000     }
1001
1002   }
1003
1004   // fill out the best track if required else clean up the fRecTracksPtr array
1005   if (!GetRecoParam()->TrackAllTracks()) {
1006     if (foundTwoClusters) {
1007       UpdateTrack(trackCandidate,bestTrackParamAtCluster1,bestTrackParamAtCluster2,bestAddChi2TrackAtCluster1,bestAddChi2TrackAtCluster2);
1008
1009       // Printout for debuging
1010       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1011         cout << "FollowTrackInStation: added the two best clusters in station(1..): " << nextStation+1 << endl;
1012         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
1013       }
1014
1015     } else if (foundOneCluster) {
1016       UpdateTrack(trackCandidate,bestTrackParamAtCluster1,bestAddChi2TrackAtCluster1);
1017
1018       // Printout for debuging
1019       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1020         cout << "FollowTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster1.GetClusterPtr()->GetChamberId()+1 << endl;
1021         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
1022       }
1023
1024     } else {
1025       delete [] clusterCh1Used;
1026       return kFALSE;
1027     }
1028
1029   } else if (foundOneCluster || foundTwoClusters) {
1030
1031     // remove obsolete track
1032     fRecTracksPtr->Remove(&trackCandidate);
1033     fNRecTracks--;
1034
1035   } else {
1036     delete [] clusterCh1Used;
1037     return kFALSE;
1038   }
1039
1040   delete [] clusterCh1Used;
1041   return kTRUE;
1042
1043 }
1044
1045   //__________________________________________________________________________
1046 Double_t AliMUONTrackReconstructorK::RunKalmanFilter(AliMUONTrackParam &trackParamAtCluster)
1047 {
1048   /// Compute new track parameters and their covariances including new cluster using kalman filter
1049   /// return the additional track chi2
1050   AliDebug(1,"Enter RunKalmanFilter");
1051
1052   // Get actual track parameters (p)
1053   TMatrixD param(trackParamAtCluster.GetParameters());
1054
1055   // Get new cluster parameters (m)
1056   AliMUONVCluster *cluster = trackParamAtCluster.GetClusterPtr();
1057   TMatrixD clusterParam(5,1);
1058   clusterParam.Zero();
1059   clusterParam(0,0) = cluster->GetX();
1060   clusterParam(2,0) = cluster->GetY();
1061
1062   // Compute the actual parameter weight (W)
1063   TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
1064   if (paramWeight.Determinant() != 0) {
1065     paramWeight.Invert();
1066   } else {
1067     AliWarning(" Determinant = 0");
1068     return 2.*AliMUONTrack::MaxChi2();
1069   }
1070
1071   // Compute the new cluster weight (U)
1072   TMatrixD clusterWeight(5,5);
1073   clusterWeight.Zero();
1074   clusterWeight(0,0) = 1. / cluster->GetErrX2();
1075   clusterWeight(2,2) = 1. / cluster->GetErrY2();
1076
1077   // Compute the new parameters covariance matrix ( (W+U)^-1 )
1078   TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
1079   if (newParamCov.Determinant() != 0) {
1080     newParamCov.Invert();
1081   } else {
1082     AliWarning(" Determinant = 0");
1083     return 2.*AliMUONTrack::MaxChi2();
1084   }
1085
1086   // Save the new parameters covariance matrix
1087   trackParamAtCluster.SetCovariances(newParamCov);
1088
1089   // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
1090   TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
1091   TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp); // U(m-p)
1092   TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
1093   newParam += param; // ((W+U)^-1)U(m-p) + p
1094
1095   // Save the new parameters
1096   trackParamAtCluster.SetParameters(newParam);
1097
1098   // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
1099   tmp = newParam; // p'
1100   tmp -= param; // (p'-p)
1101   TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp); // W(p'-p)
1102   TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
1103   tmp = newParam; // p'
1104   tmp -= clusterParam; // (p'-m)
1105   TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp); // U(p'-m)
1106   addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
1107
1108   return addChi2Track(0,0);
1109
1110 }
1111
1112   //__________________________________________________________________________
1113 void AliMUONTrackReconstructorK::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster, Double_t addChi2)
1114 {
1115   /// Add 1 cluster to the track candidate
1116   /// Update chi2 of the track
1117
1118   // Flag cluster as being (not) removable
1119   if (GetRecoParam()->RequestStation(trackParamAtCluster.GetClusterPtr()->GetChamberId()/2))
1120     trackParamAtCluster.SetRemovable(kFALSE);
1121   else trackParamAtCluster.SetRemovable(kTRUE);
1122   trackParamAtCluster.SetLocalChi2(0.); // --> Local chi2 not used
1123
1124   // Update the track chi2 into trackParamAtCluster
1125   trackParamAtCluster.SetTrackChi2(track.GetGlobalChi2() + addChi2);
1126
1127   // Update the chi2 of the new track
1128   track.SetGlobalChi2(trackParamAtCluster.GetTrackChi2());
1129
1130   // Update array of TrackParamAtCluster
1131   track.AddTrackParamAtCluster(trackParamAtCluster,*(trackParamAtCluster.GetClusterPtr()));
1132
1133 }
1134
1135   //__________________________________________________________________________
1136 void AliMUONTrackReconstructorK::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster1, AliMUONTrackParam &trackParamAtCluster2,
1137                                              Double_t addChi2AtCluster1, Double_t addChi2AtCluster2)
1138 {
1139   /// Add 2 clusters to the track candidate (order is important)
1140   /// Update track and local chi2
1141
1142   // Update local chi2 at first cluster
1143   AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr();
1144   Double_t deltaX = trackParamAtCluster1.GetNonBendingCoor() - cluster1->GetX();
1145   Double_t deltaY = trackParamAtCluster1.GetBendingCoor() - cluster1->GetY();
1146   Double_t localChi2AtCluster1 = deltaX*deltaX / cluster1->GetErrX2() +
1147                                  deltaY*deltaY / cluster1->GetErrY2();
1148   trackParamAtCluster1.SetLocalChi2(localChi2AtCluster1);
1149
1150   // Flag first cluster as being removable
1151   trackParamAtCluster1.SetRemovable(kTRUE);
1152
1153   // Update local chi2 at second cluster
1154   AliMUONVCluster* cluster2 = trackParamAtCluster2.GetClusterPtr();
1155   AliMUONTrackParam extrapTrackParamAtCluster2(trackParamAtCluster1);
1156   AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParamAtCluster2, trackParamAtCluster2.GetZ());
1157   deltaX = extrapTrackParamAtCluster2.GetNonBendingCoor() - cluster2->GetX();
1158   deltaY = extrapTrackParamAtCluster2.GetBendingCoor() - cluster2->GetY();
1159   Double_t localChi2AtCluster2 = deltaX*deltaX / cluster2->GetErrX2() +
1160                                  deltaY*deltaY / cluster2->GetErrY2();
1161   trackParamAtCluster2.SetLocalChi2(localChi2AtCluster2);
1162
1163   // Flag second cluster as being removable
1164   trackParamAtCluster2.SetRemovable(kTRUE);
1165
1166   // Update the track chi2 into trackParamAtCluster2
1167   trackParamAtCluster2.SetTrackChi2(track.GetGlobalChi2() + addChi2AtCluster2);
1168
1169   // Update the track chi2 into trackParamAtCluster1
1170   trackParamAtCluster1.SetTrackChi2(trackParamAtCluster2.GetTrackChi2() + addChi2AtCluster1);
1171
1172   // Update the chi2 of the new track
1173   track.SetGlobalChi2(trackParamAtCluster1.GetTrackChi2());
1174
1175   // Update array of trackParamAtCluster
1176   track.AddTrackParamAtCluster(trackParamAtCluster1,*cluster1);
1177   track.AddTrackParamAtCluster(trackParamAtCluster2,*cluster2);
1178
1179 }
1180
1181   //__________________________________________________________________________
1182 Bool_t AliMUONTrackReconstructorK::RecoverTrack(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation)
1183 {
1184   /// Try to recover the track candidate in the next station
1185   /// by removing the worst of the two clusters attached in the current station
1186   /// Return kTRUE if recovering succeeds
1187   AliDebug(1,"Enter RecoverTrack");
1188
1189   // Do not try to recover track until we have attached cluster(s) on station(1..) 3
1190   if (nextStation > 1) return kFALSE;
1191
1192   Int_t worstClusterNumber = -1;
1193   Double_t localChi2, worstLocalChi2 = -1.;
1194
1195   // Look for the cluster to remove
1196   for (Int_t clusterNumber = 0; clusterNumber < 2; clusterNumber++) {
1197     AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(clusterNumber);
1198
1199     // check if current cluster is in the previous station
1200     if (trackParamAtCluster->GetClusterPtr()->GetChamberId()/2 != nextStation+1) break;
1201
1202     // check if current cluster is removable
1203     if (!trackParamAtCluster->IsRemovable()) return kFALSE;
1204
1205     // reset the current cluster as being not removable if it is on a required station
1206     if (GetRecoParam()->RequestStation(nextStation+1)) trackParamAtCluster->SetRemovable(kFALSE);
1207
1208     // Pick up cluster with the worst chi2
1209     localChi2 = trackParamAtCluster->GetLocalChi2();
1210     if (localChi2 > worstLocalChi2) {
1211       worstLocalChi2 = localChi2;
1212       worstClusterNumber = clusterNumber;
1213     }
1214   }
1215
1216   // check if worst cluster found
1217   if (worstClusterNumber < 0) return kFALSE;
1218
1219   // Remove the worst cluster
1220   trackCandidate.RemoveTrackParamAtCluster((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(worstClusterNumber));
1221
1222   // Re-calculate track parameters at the (new) first cluster
1223   if (!RetracePartialTrack(trackCandidate,(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(1))) return kFALSE;
1224
1225   // skip track out of limits
1226   if (!IsAcceptable(*((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First()))) return kFALSE;
1227
1228   // Look for new cluster(s) in next station
1229   return FollowTrackInStation(trackCandidate, clusterStore, nextStation);
1230
1231 }
1232
1233   //__________________________________________________________________________
1234 Bool_t AliMUONTrackReconstructorK::RunSmoother(AliMUONTrack &track)
1235 {
1236   /// Compute new track parameters and their covariances using smoother
1237   AliDebug(1,"Enter UseSmoother");
1238
1239   AliMUONTrackParam *previousTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->First();
1240
1241   // Smoothed parameters and covariances at first cluster = filtered parameters and covariances
1242   previousTrackParam->SetSmoothParameters(previousTrackParam->GetParameters());
1243   previousTrackParam->SetSmoothCovariances(previousTrackParam->GetCovariances());
1244
1245   AliMUONTrackParam *currentTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(previousTrackParam);
1246
1247   // Save local chi2 at first cluster = last additional chi2 provided by Kalman
1248   previousTrackParam->SetLocalChi2(previousTrackParam->GetTrackChi2() - currentTrackParam->GetTrackChi2());
1249
1250   // if the track contains only 2 clusters simply copy the filtered parameters
1251   if (track.GetNClusters() == 2) {
1252     currentTrackParam->SetSmoothParameters(currentTrackParam->GetParameters());
1253     currentTrackParam->SetSmoothCovariances(currentTrackParam->GetCovariances());
1254     currentTrackParam->SetLocalChi2(currentTrackParam->GetTrackChi2());
1255     return kTRUE;
1256   }
1257
1258   while (currentTrackParam) {
1259
1260     // Get variables
1261     const TMatrixD &extrapParameters          = previousTrackParam->GetExtrapParameters();  // X(k+1 k)
1262     const TMatrixD &filteredParameters        = currentTrackParam->GetParameters();         // X(k k)
1263     const TMatrixD &previousSmoothParameters  = previousTrackParam->GetSmoothParameters();  // X(k+1 n)
1264     const TMatrixD &propagator                = previousTrackParam->GetPropagator();        // F(k)
1265     const TMatrixD &extrapCovariances         = previousTrackParam->GetExtrapCovariances(); // C(k+1 k)
1266     const TMatrixD &filteredCovariances       = currentTrackParam->GetCovariances();        // C(k k)
1267     const TMatrixD &previousSmoothCovariances = previousTrackParam->GetSmoothCovariances(); // C(k+1 n)
1268
1269     // Compute smoother gain: A(k) = C(kk) * F(k)^t * (C(k+1 k))^-1
1270     TMatrixD extrapWeight(extrapCovariances);
1271     if (extrapWeight.Determinant() != 0) {
1272       extrapWeight.Invert(); // (C(k+1 k))^-1
1273     } else {
1274       AliWarning(" Determinant = 0");
1275       return kFALSE;
1276     }
1277     TMatrixD smootherGain(filteredCovariances,TMatrixD::kMultTranspose,propagator); // C(kk) * F(k)^t
1278     smootherGain *= extrapWeight; // C(kk) * F(k)^t * (C(k+1 k))^-1
1279
1280     // Compute smoothed parameters: X(k n) = X(k k) + A(k) * (X(k+1 n) - X(k+1 k))
1281     TMatrixD tmpParam(previousSmoothParameters,TMatrixD::kMinus,extrapParameters); // X(k+1 n) - X(k+1 k)
1282     TMatrixD smoothParameters(smootherGain,TMatrixD::kMult,tmpParam); // A(k) * (X(k+1 n) - X(k+1 k))
1283     smoothParameters += filteredParameters; // X(k k) + A(k) * (X(k+1 n) - X(k+1 k))
1284
1285     // Save smoothed parameters
1286     currentTrackParam->SetSmoothParameters(smoothParameters);
1287
1288     // Compute smoothed covariances: C(k n) = C(k k) + A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
1289     TMatrixD tmpCov(previousSmoothCovariances,TMatrixD::kMinus,extrapCovariances); // C(k+1 n) - C(k+1 k)
1290     TMatrixD tmpCov2(tmpCov,TMatrixD::kMultTranspose,smootherGain); // (C(k+1 n) - C(k+1 k)) * (A(k))^t
1291     TMatrixD smoothCovariances(smootherGain,TMatrixD::kMult,tmpCov2); // A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
1292     smoothCovariances += filteredCovariances; // C(k k) + A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
1293
1294     // Save smoothed covariances
1295     currentTrackParam->SetSmoothCovariances(smoothCovariances);
1296
1297     // Compute smoothed residual: r(k n) = cluster - X(k n)
1298     AliMUONVCluster* cluster = currentTrackParam->GetClusterPtr();
1299     TMatrixD smoothResidual(2,1);
1300     smoothResidual.Zero();
1301     smoothResidual(0,0) = cluster->GetX() - smoothParameters(0,0);
1302     smoothResidual(1,0) = cluster->GetY() - smoothParameters(2,0);
1303
1304     // Compute weight of smoothed residual: W(k n) = (clusterCov - C(k n))^-1
1305     TMatrixD smoothResidualWeight(2,2);
1306     smoothResidualWeight(0,0) = cluster->GetErrX2() - smoothCovariances(0,0);
1307     smoothResidualWeight(0,1) = - smoothCovariances(0,2);
1308     smoothResidualWeight(1,0) = - smoothCovariances(2,0);
1309     smoothResidualWeight(1,1) = cluster->GetErrY2() - smoothCovariances(2,2);
1310     if (smoothResidualWeight.Determinant() != 0) {
1311       smoothResidualWeight.Invert();
1312     } else {
1313       AliWarning(" Determinant = 0");
1314       return kFALSE;
1315     }
1316
1317     // Compute local chi2 = (r(k n))^t * W(k n) * r(k n)
1318     TMatrixD tmpChi2(smoothResidual,TMatrixD::kTransposeMult,smoothResidualWeight); // (r(k n))^t * W(k n)
1319     TMatrixD localChi2(tmpChi2,TMatrixD::kMult,smoothResidual); // (r(k n))^t * W(k n) * r(k n)
1320
1321     // Save local chi2
1322     currentTrackParam->SetLocalChi2(localChi2(0,0));
1323
1324     previousTrackParam = currentTrackParam;
1325     currentTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(previousTrackParam);
1326   }
1327
1328   return kTRUE;
1329
1330 }
1331
1332   //__________________________________________________________________________
1333 Bool_t AliMUONTrackReconstructorK::ComplementTracks(const AliMUONVClusterStore& clusterStore)
1334 {
1335   /// Complete tracks by adding missing clusters (if there is an overlap between
1336   /// two detection elements, the track may have two clusters in the same chamber).
1337   /// Recompute track parameters and covariances at each clusters.
1338   /// Remove tracks getting abnormal (i.e. extrapolation failed...) after being complemented.
1339   /// Return kTRUE if one or more tracks have been complemented or removed.
1340   AliDebug(1,"Enter ComplementTracks");
1341
1342   Int_t chamberId, detElemId;
1343   Double_t chi2OfCluster, addChi2TrackAtCluster, bestAddChi2TrackAtCluster;
1344   Double_t maxChi2OfCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
1345                                    GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
1346   Bool_t foundOneCluster, trackModified, hasChanged = kFALSE;
1347   AliMUONVCluster *cluster;
1348   AliMUONTrackParam *trackParam, *previousTrackParam, *nextTrackParam, trackParamAtCluster, bestTrackParamAtCluster;
1349   AliMUONTrack *nextTrack;
1350
1351   AliMUONTrack *track = (AliMUONTrack*) fRecTracksPtr->First();
1352   while (track) {
1353     trackModified = kFALSE;
1354
1355     trackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->First();
1356     previousTrackParam = trackParam;
1357     while (trackParam) {
1358       foundOneCluster = kFALSE;
1359       bestAddChi2TrackAtCluster = AliMUONTrack::MaxChi2();
1360       chamberId = trackParam->GetClusterPtr()->GetChamberId();
1361       detElemId = trackParam->GetClusterPtr()->GetDetElemId();
1362
1363       // prepare nextTrackParam before adding new cluster because of the sorting
1364       nextTrackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->After(trackParam);
1365
1366       // Create iterators to loop over clusters in current chamber
1367       TIter nextInCh(clusterStore.CreateChamberIterator(chamberId,chamberId));
1368
1369       // look for one second candidate in the same chamber
1370       while ( ( cluster = static_cast<AliMUONVCluster*>(nextInCh()) ) ) {
1371
1372         // look for a cluster in another detection element
1373         if (cluster->GetDetElemId() == detElemId) continue;
1374
1375         // try to add the current cluster fast
1376         if (!TryOneClusterFast(*trackParam, cluster)) continue;
1377
1378         // try to add the current cluster accurately
1379         // never use track parameters at last cluster because the covariance matrix is meaningless
1380         if (nextTrackParam) chi2OfCluster = TryOneCluster(*trackParam, cluster, trackParamAtCluster);
1381         else chi2OfCluster = TryOneCluster(*previousTrackParam, cluster, trackParamAtCluster);
1382
1383         // if good chi2 then consider to add this cluster to the track
1384         if (chi2OfCluster < maxChi2OfCluster) {
1385
1386           // Compute local track parameters including current cluster using kalman filter
1387           addChi2TrackAtCluster = RunKalmanFilter(trackParamAtCluster);
1388
1389           // keep track of the best cluster
1390           if (addChi2TrackAtCluster < bestAddChi2TrackAtCluster) {
1391             bestAddChi2TrackAtCluster = addChi2TrackAtCluster;
1392             bestTrackParamAtCluster = trackParamAtCluster;
1393             foundOneCluster = kTRUE;
1394           }
1395
1396         }
1397
1398       }
1399
1400       // add new cluster if any
1401       if (foundOneCluster) {
1402
1403         // Printout for debuging
1404         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1405           cout << "ComplementTracks: found one cluster in chamber(1..): " << chamberId+1 << endl;
1406           bestTrackParamAtCluster.GetClusterPtr()->Print();
1407           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
1408             cout<<endl<<"Track parameters and covariances at cluster:"<<endl;
1409             bestTrackParamAtCluster.GetParameters().Print();
1410             bestTrackParamAtCluster.GetCovariances().Print();
1411           }
1412         }
1413
1414         trackParam->SetRemovable(kTRUE);
1415         bestTrackParamAtCluster.SetRemovable(kTRUE);
1416         track->AddTrackParamAtCluster(bestTrackParamAtCluster,*(bestTrackParamAtCluster.GetClusterPtr()));
1417         trackModified = kTRUE;
1418         hasChanged = kTRUE;
1419       }
1420
1421       previousTrackParam = trackParam;
1422       trackParam = nextTrackParam;
1423     }
1424
1425     // prepare next track
1426     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track);
1427
1428     // re-compute track parameters using kalman filter if needed
1429     if (trackModified && !RetraceTrack(*track,kTRUE)) {
1430       AliWarning("track modified but problem occur during refitting --> remove track");
1431       fRecTracksPtr->Remove(track);
1432       fNRecTracks--;
1433     }
1434
1435     track = nextTrack;
1436   }
1437
1438   return hasChanged;
1439
1440 }
1441
1442 //__________________________________________________________________________
1443 void AliMUONTrackReconstructorK::ImproveTrack(AliMUONTrack &track)
1444 {
1445   /// Improve the given track by removing removable clusters with local chi2 highter than the defined cut
1446   /// Removable clusters are identified by the method AliMUONTrack::TagRemovableClusters()
1447   /// Recompute track parameters and covariances at the remaining clusters
1448   /// and if something goes wrong (i.e. extrapolation failed...) set track chi2 to max value
1449   AliDebug(1,"Enter ImproveTrack");
1450
1451   Double_t localChi2, worstLocalChi2;
1452   AliMUONTrackParam *trackParamAtCluster, *worstTrackParamAtCluster, *nextTrackParam, *next2nextTrackParam;
1453   Int_t nextChamber, next2nextChamber;
1454   Bool_t smoothed;
1455   Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForImprovement() *
1456                        GetRecoParam()->GetSigmaCutForImprovement();
1457
1458   while (!track.IsImproved()) {
1459
1460     // identify removable clusters
1461     track.TagRemovableClusters(GetRecoParam()->RequestedStationMask());
1462
1463     // Run smoother if required
1464     smoothed = kFALSE;
1465     if (GetRecoParam()->UseSmoother()) smoothed = RunSmoother(track);
1466
1467     // Use standard procedure to compute local chi2 if smoother not required or not working
1468     if (!smoothed) {
1469
1470       // Update track parameters and covariances
1471       if (!track.UpdateCovTrackParamAtCluster()) {
1472         AliWarning("unable to update track parameters and covariances --> stop improvement");
1473         // restore the kalman parameters
1474         RetraceTrack(track,kTRUE);
1475         break;
1476       }
1477
1478       // Compute local chi2 of each clusters
1479       track.ComputeLocalChi2(kTRUE);
1480     }
1481
1482     // Look for the cluster to remove
1483     worstTrackParamAtCluster = 0x0;
1484     worstLocalChi2 = -1.;
1485     trackParamAtCluster = (AliMUONTrackParam*)track.GetTrackParamAtCluster()->First();
1486     while (trackParamAtCluster) {
1487
1488       // save parameters into smooth parameters in case of smoother did not work properly
1489       if (GetRecoParam()->UseSmoother() && !smoothed) {
1490         trackParamAtCluster->SetSmoothParameters(trackParamAtCluster->GetParameters());
1491         trackParamAtCluster->SetSmoothCovariances(trackParamAtCluster->GetCovariances());
1492       }
1493
1494       // Pick up cluster with the worst chi2
1495       localChi2 = trackParamAtCluster->GetLocalChi2();
1496       if (localChi2 > worstLocalChi2) {
1497         worstLocalChi2 = localChi2;
1498         worstTrackParamAtCluster = trackParamAtCluster;
1499       }
1500
1501       trackParamAtCluster = (AliMUONTrackParam*)track.GetTrackParamAtCluster()->After(trackParamAtCluster);
1502     }
1503
1504     // Check whether the worst chi2 is under requirement or not
1505     if (worstLocalChi2 < 2. * sigmaCut2) { // 2 because 2 quantities in chi2
1506       track.SetImproved(kTRUE);
1507       break;
1508     }
1509
1510     // if the worst cluster is not removable then stop improvement
1511     if (!worstTrackParamAtCluster->IsRemovable()) {
1512       // restore the kalman parameters in case they have been lost
1513       if (!smoothed) RetraceTrack(track,kTRUE);
1514       break;
1515     }
1516
1517     // get track parameters at cluster next to the one to be removed
1518     nextTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(worstTrackParamAtCluster);
1519
1520     // Remove the worst cluster
1521     track.RemoveTrackParamAtCluster(worstTrackParamAtCluster);
1522
1523     // Re-calculate track parameters
1524     // - from the cluster immediately downstream the one suppressed
1525     // - or from the begining - if parameters have been re-computed using the standard method (kalman parameters have been lost)
1526     //                        - or if the removed cluster was used to compute the tracking seed
1527     Bool_t normalExtrap;
1528     if (smoothed && nextTrackParam) {
1529
1530       nextChamber = nextTrackParam->GetClusterPtr()->GetChamberId();
1531       next2nextTrackParam = nextTrackParam;
1532       do {
1533
1534         next2nextChamber = next2nextTrackParam->GetClusterPtr()->GetChamberId();
1535         next2nextTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(next2nextTrackParam);
1536
1537       } while (next2nextTrackParam && (next2nextChamber == nextChamber));
1538
1539       if (next2nextChamber == nextChamber) normalExtrap = RetraceTrack(track,kTRUE);
1540       else normalExtrap = RetracePartialTrack(track,nextTrackParam);
1541
1542     } else normalExtrap = RetraceTrack(track,kTRUE);
1543
1544     // stop in case of extrapolation problem
1545     if (!normalExtrap) {
1546       AliWarning("track partially improved but problem occur during refitting --> stop improvement");
1547       break;
1548     }
1549
1550     // Printout for debuging
1551     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1552       cout << "ImproveTracks: track " << fRecTracksPtr->IndexOf(&track)+1 << " improved " << endl;
1553     }
1554
1555   }
1556
1557 }
1558
1559 //__________________________________________________________________________
1560 Bool_t AliMUONTrackReconstructorK::FinalizeTrack(AliMUONTrack &track)
1561 {
1562   /// Update track parameters and covariances at each attached cluster
1563   /// using smoother if required, if not already done
1564   /// return kFALSE if the track cannot be extrapolated uo to the last chamber
1565
1566   AliMUONTrackParam *trackParamAtCluster;
1567   Bool_t smoothed = kFALSE;
1568
1569   // update track parameters (using smoother if required) if not already done
1570   if (track.IsImproved()) smoothed = GetRecoParam()->UseSmoother();
1571   else {
1572     if (GetRecoParam()->UseSmoother()) smoothed = RunSmoother(track);
1573     if (!smoothed) {
1574       if (track.UpdateCovTrackParamAtCluster()) track.ComputeLocalChi2(kTRUE);
1575       else {
1576         AliWarning("finalization failed due to extrapolation problem");
1577         return kFALSE;
1578       }
1579     }
1580   }
1581
1582   // copy smoothed parameters and covariances if any
1583   if (smoothed) {
1584
1585     trackParamAtCluster = (AliMUONTrackParam*) (track.GetTrackParamAtCluster()->First());
1586     while (trackParamAtCluster) {
1587
1588       trackParamAtCluster->SetParameters(trackParamAtCluster->GetSmoothParameters());
1589       trackParamAtCluster->SetCovariances(trackParamAtCluster->GetSmoothCovariances());
1590
1591       trackParamAtCluster = (AliMUONTrackParam*) (track.GetTrackParamAtCluster()->After(trackParamAtCluster));
1592     }
1593
1594   }
1595
1596   return kTRUE;
1597
1598 }
1599
1600   //__________________________________________________________________________
1601 Bool_t AliMUONTrackReconstructorK::RefitTrack(AliMUONTrack &track, Bool_t enableImprovement)
1602 {
1603   /// re-fit the given track
1604   AliDebug(1,"Enter RefitTrack");
1605
1606   // check validity of the track (i.e. at least 2 chambers hit on stations 4 and 5)
1607   if (!track.IsValid(0)) {
1608     AliWarning("the track is not valid --> unable to refit");
1609     return kFALSE;
1610   }
1611
1612   // re-compute track parameters and covariances using Kalman filter
1613   if (!RetraceTrack(track,kTRUE)) {
1614     AliWarning("bad track refitting due to extrapolation failure");
1615     return kFALSE;
1616   }
1617
1618   // Improve the reconstructed tracks if required
1619   track.SetImproved(kFALSE);
1620   if (enableImprovement && GetRecoParam()->ImproveTracks()) ImproveTrack(track);
1621
1622   // Fill AliMUONTrack data members
1623   if (track.GetGlobalChi2() < AliMUONTrack::MaxChi2()) return FinalizeTrack(track);
1624   else {
1625     AliWarning("track not finalized due to extrapolation failure");
1626     return kFALSE;
1627   }
1628
1629 }
1630