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
  41 #include "AliMpArea.h"
  42
  43 #include "AliLog.h"
  44
  45 #include <Riostream.h>
  46 #include <TMath.h>
  47 #include <TMatrixD.h>
  48
  49 /// \cond CLASSIMP
  50 ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context
  51 /// \endcond
  52
  53   //__________________________________________________________________________
  54 AliMUONTrackReconstructorK::AliMUONTrackReconstructorK(AliMUONVClusterServer& clusterServer)
  55   : AliMUONVTrackReconstructor(clusterServer)
  56 {
  57   /// Constructor
  58 }
  59
  60   //__________________________________________________________________________
  61 AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK()
  62 {
  63 /// Destructor
  64 }
  65
  66   //__________________________________________________________________________
  67 void AliMUONTrackReconstructorK::MakeTrackCandidates(AliMUONVClusterStore& clusterStore)
  68 {
  69   /// To make track candidates (assuming linear propagation if the flag fgkMakeTrackCandidatesFast is set to kTRUE):
  70   /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4.
  71   /// Good candidates are made of at least three clusters.
  72   /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
  73
  74   TClonesArray *segments;
  75   AliMUONTrack *track;
  76   Int_t iCandidate = 0;
  77   Bool_t clusterFound;
  78
  79   AliDebug(1,"Enter MakeTrackCandidates");
  80
  81   // Unless we're doing combined tracking, we'll clusterize all stations at once
  82   Int_t firstChamber(0);
  83   Int_t lastChamber(9);
  84
  85   if (AliMUONReconstructor::GetRecoParam()->CombineClusterTrackReco()) {
  86     // ... Here's the exception : ask the clustering to reconstruct
  87     // clusters *only* in station 4 and 5 for combined tracking
  88     firstChamber = 6;
  89   }
  90
  91   for (Int_t i = firstChamber; i <= lastChamber; ++i )
  92   {
  93     fClusterServer.Clusterize(i, clusterStore, AliMpArea());
  94   }
  95
  96   // Loop over stations(1..) 5 and 4 and make track candidates
  97   for (Int_t istat=4; istat>=3; istat--) {
  98
  99     // Make segments in the station
 100     segments = MakeSegmentsInStation(clusterStore, istat);
 101
 102     // Loop over segments
 103     for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++) {
 104       AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
 105
 106       // Transform segments to tracks and put them at the end of fRecTracksPtr
 107       track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack((AliMUONObjectPair*)((*segments)[iseg]));
 108       fNRecTracks++;
 109
 110       // Printout for debuging
 111       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 112         cout<<endl<<"Track parameter covariances at first cluster:"<<endl;
 113         ((AliMUONTrackParam*) (track->GetTrackParamAtCluster()->First()))->GetCovariances().Print();
 114       }
 115
 116       // Look for compatible cluster(s) in the other station
 117       if (AliMUONReconstructor::GetRecoParam()->MakeTrackCandidatesFast())
 118         clusterFound = FollowLinearTrackInStation(*track, clusterStore, 7-istat);
 119       else {
 120         // First recompute track parameters and covariances on station(1..) 5 using Kalman filter
 121         // (to make sure all tracks are treated in the same way)
 122         if (istat == 4) RetraceTrack(*track, kFALSE);
 123         clusterFound = FollowTrackInStation(*track, clusterStore, 7-istat);
 124       }
 125
 126       // Remove track if no cluster found
 127       if (!clusterFound) {
 128         fRecTracksPtr->Remove(track);
 129         fNRecTracks--;
 130       }
 131
 132     }
 133     // delete the array of segments
 134     delete segments;
 135   }
 136
 137
 138   // Retrace tracks using Kalman filter and remove bad ones
 139   if (AliMUONReconstructor::GetRecoParam()->MakeTrackCandidatesFast()) {
 140     fRecTracksPtr->Compress(); // this is essential before checking tracks
 141
 142     Int_t currentNRecTracks = fNRecTracks;
 143     for (Int_t iRecTrack = 0; iRecTrack < currentNRecTracks; iRecTrack++) {
 144       track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
 145
 146       // Recompute track parameters and covariances using Kalman filter
 147       RetraceTrack(*track,kTRUE);
 148
 149       // Remove the track if the normalized chi2 is too high
 150       if (track->GetNormalizedChi2() > AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
 151                                        AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking()) {
 152         fRecTracksPtr->Remove(track);
 153         fNRecTracks--;
 154       }
 155
 156     }
 157
 158   }
 159
 160   fRecTracksPtr->Compress(); // this is essential before checking tracks
 161
 162   // Keep all different tracks or only the best ones as required
 163   if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) RemoveIdenticalTracks();
 164   else RemoveDoubleTracks();
 165
 166   AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
 167
 168 }
 169
 170   //__________________________________________________________________________
 171 void AliMUONTrackReconstructorK::RetraceTrack(AliMUONTrack &trackCandidate, Bool_t resetSeed)
 172 {
 173   /// Re-run the kalman filter from the most downstream cluster to the most uptream one
 174   AliDebug(1,"Enter RetraceTrack");
 175
 176   AliMUONTrackParam* lastTrackParam = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtCluster()->Last();
 177
 178   // Reset the "seed" (= track parameters and their covariances at last cluster) if required
 179   if (resetSeed) {
 180
 181     // parameters at last cluster
 182     AliMUONVCluster* cluster2 = lastTrackParam->GetClusterPtr();
 183     Double_t x2 = cluster2->GetX();
 184     Double_t y2 = cluster2->GetY();
 185     Double_t z2 = cluster2->GetZ();
 186
 187     // parameters at last but one cluster
 188     AliMUONTrackParam* previousTrackParam = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtCluster()->Before(lastTrackParam);
 189     AliMUONVCluster* cluster1 = previousTrackParam->GetClusterPtr();
 190     Double_t x1 = cluster1->GetX();
 191     Double_t y1 = cluster1->GetY();
 192     Double_t z1 = cluster1->GetZ();
 193
 194     // reset track parameters
 195     Double_t dZ = z1 - z2;
 196     lastTrackParam->SetNonBendingCoor(x2);
 197     lastTrackParam->SetBendingCoor(y2);
 198     lastTrackParam->SetZ(z2);
 199     lastTrackParam->SetNonBendingSlope((x1 - x2) / dZ);
 200     lastTrackParam->SetBendingSlope((y1 - y2) / dZ);
 201     Double_t bendingImpact = y2 - z2 * lastTrackParam->GetBendingSlope();
 202     lastTrackParam->SetInverseBendingMomentum(1. / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact));
 203
 204     // => Reset track parameter covariances at last cluster (as if the other clusters did not exist)
 205     TMatrixD lastParamCov(5,5);
 206     lastParamCov.Zero();
 207     // Non bending plane
 208     lastParamCov(0,0) = cluster2->GetErrX2();
 209     lastParamCov(1,1) = 100. * ( cluster1->GetErrX2() + cluster2->GetErrX2() ) / dZ / dZ;
 210     // Bending plane
 211     lastParamCov(2,2) = cluster2->GetErrY2();
 212     lastParamCov(3,3) = 100. * ( cluster1->GetErrY2() + cluster2->GetErrY2() ) / dZ / dZ;
 213     // Inverse bending momentum (50% error)
 214     lastParamCov(4,4) = 0.5*lastTrackParam->GetInverseBendingMomentum() * 0.5*lastTrackParam->GetInverseBendingMomentum();
 215     lastTrackParam->SetCovariances(lastParamCov);
 216
 217     // Reset the track chi2
 218     lastTrackParam->SetTrackChi2(0.);
 219
 220   }
 221
 222   // Redo the tracking
 223   RetracePartialTrack(trackCandidate, lastTrackParam);
 224
 225 }
 226
 227   //__________________________________________________________________________
 228 void AliMUONTrackReconstructorK::RetracePartialTrack(AliMUONTrack &trackCandidate, const AliMUONTrackParam* startingTrackParam)
 229 {
 230   /// Re-run the kalman filter from the cluster attached to startingTrackParam to the most uptream cluster
 231   AliDebug(1,"Enter RetracePartialTrack");
 232
 233   // Printout for debuging
 234   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 235     cout << "RetracePartialTrack: track chi2 before re-tracking: " << trackCandidate.GetGlobalChi2() << endl;
 236   }
 237
 238   // Reset the track chi2
 239   trackCandidate.SetGlobalChi2(startingTrackParam->GetTrackChi2());
 240
 241   // loop over attached clusters until the first one and recompute track parameters and covariances using kalman filter
 242   Int_t expectedChamber = startingTrackParam->GetClusterPtr()->GetChamberId() - 1;
 243   Int_t currentChamber;
 244   Double_t addChi2TrackAtCluster;
 245   AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtCluster()->Before(startingTrackParam);
 246   while (trackParamAtCluster) {
 247
 248     // reset track parameters and their covariances
 249     trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
 250     trackParamAtCluster->SetZ(startingTrackParam->GetZ());
 251     trackParamAtCluster->SetCovariances(startingTrackParam->GetCovariances());
 252
 253     // add MCS effect
 254     AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
 255
 256     // reset propagator for smoother
 257     if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) trackParamAtCluster->ResetPropagator();
 258
 259     // add MCS in missing chambers if any (at most 2 chambers can be missing according to tracking criteria)
 260     currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
 261     while (currentChamber < expectedChamber) {
 262       // extrapolation to the missing chamber (update the propagator)
 263       AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, AliMUONConstants::DefaultChamberZ(expectedChamber),
 264                                        AliMUONReconstructor::GetRecoParam()->UseSmoother());
 265       // add MCS effect
 266       AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
 267       expectedChamber--;
 268     }
 269
 270     // extrapolation to the plane of the cluster attached to the current trackParamAtCluster (update the propagator)
 271     AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ(),
 272                                      AliMUONReconstructor::GetRecoParam()->UseSmoother());
 273
 274     if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) {
 275       // save extrapolated parameters for smoother
 276       trackParamAtCluster->SetExtrapParameters(trackParamAtCluster->GetParameters());
 277
 278       // save extrapolated covariance matrix for smoother
 279       trackParamAtCluster->SetExtrapCovariances(trackParamAtCluster->GetCovariances());
 280     }
 281
 282     // Compute new track parameters including "clusterCh2" using kalman filter
 283     addChi2TrackAtCluster = RunKalmanFilter(*trackParamAtCluster);
 284
 285     // Update the track chi2
 286     trackCandidate.SetGlobalChi2(trackCandidate.GetGlobalChi2() + addChi2TrackAtCluster);
 287     trackParamAtCluster->SetTrackChi2(trackCandidate.GetGlobalChi2());
 288
 289     // prepare next step
 290     expectedChamber = currentChamber - 1;
 291     startingTrackParam = trackParamAtCluster;
 292     trackParamAtCluster = (AliMUONTrackParam*) (trackCandidate.GetTrackParamAtCluster()->Before(startingTrackParam));
 293   }
 294
 295   // Printout for debuging
 296   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 297     cout << "RetracePartialTrack: track chi2 after re-tracking: " << trackCandidate.GetGlobalChi2() << endl;
 298   }
 299
 300 }
 301
 302   //__________________________________________________________________________
 303 void AliMUONTrackReconstructorK::FollowTracks(AliMUONVClusterStore& clusterStore)
 304 {
 305   /// Follow tracks in stations(1..) 3, 2 and 1
 306   AliDebug(1,"Enter FollowTracks");
 307
 308   AliMUONTrack *track;
 309   Int_t currentNRecTracks;
 310   Bool_t clusterFound;
 311
 312   for (Int_t station = 2; station >= 0; station--) {
 313
 314     // Save the actual number of reconstructed track in case of
 315     // tracks are added or suppressed during the tracking procedure
 316     // !! Do not compress fRecTracksPtr until the end of the loop over tracks !!
 317     currentNRecTracks = fNRecTracks;
 318
 319     for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) {
 320       AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1));
 321
 322       track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
 323
 324       // Printout for debuging
 325       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 326         cout<<endl<<"Track parameter covariances at first cluster:"<<endl;
 327         ((AliMUONTrackParam*) (track->GetTrackParamAtCluster()->First()))->GetCovariances().Print();
 328       }
 329
 330       // Look for compatible cluster(s) in station(0..) "station"
 331       clusterFound = FollowTrackInStation(*track, clusterStore, station);
 332
 333       // Try to recover track if required
 334       if (!clusterFound && AliMUONReconstructor::GetRecoParam()->RecoverTracks())
 335         clusterFound = RecoverTrack(*track, clusterStore, station);
 336
 337       // remove track if no cluster found
 338       if (!clusterFound) {
 339         fRecTracksPtr->Remove(track);
 340         fNRecTracks--;
 341       }
 342
 343     }
 344
 345     fRecTracksPtr->Compress(); // this is essential before checking tracks
 346
 347     // Keep only the best tracks if required
 348     if (!AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
 349
 350   }
 351
 352 }
 353
 354   //__________________________________________________________________________
 355 Bool_t AliMUONTrackReconstructorK::FollowTrackInStation(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation)
 356 {
 357   /// Follow trackCandidate in station(0..) nextStation and search for compatible cluster(s)
 358   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
 359   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
 360   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
 361   ///         Remove the obsolete "trackCandidate" at the end.
 362   /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
 363   /// return kTRUE if new cluster(s) have been found (otherwise return kFALSE)
 364   AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1));
 365
 366   // Order the chamber according to the propagation direction (tracking starts with chamber 2):
 367   // - nextStation == station(1...) 5 => forward propagation
 368   // - nextStation < station(1...) 5 => backward propagation
 369   Int_t ch1, ch2;
 370   if (nextStation==4) {
 371     ch1 = 2*nextStation+1;
 372     ch2 = 2*nextStation;
 373   } else {
 374     ch1 = 2*nextStation;
 375     ch2 = 2*nextStation+1;
 376   }
 377
 378   Double_t chi2OfCluster;
 379   Double_t maxChi2OfCluster = 2. * AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
 380                                    AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
 381   Double_t addChi2TrackAtCluster1;
 382   Double_t addChi2TrackAtCluster2;
 383   Double_t bestAddChi2TrackAtCluster1 = 1.e10;
 384   Double_t bestAddChi2TrackAtCluster2 = 1.e10;
 385   Bool_t foundOneCluster = kFALSE;
 386   Bool_t foundTwoClusters = kFALSE;
 387   AliMUONTrack *newTrack = 0x0;
 388   AliMUONVCluster *clusterCh1, *clusterCh2;
 389   AliMUONTrackParam extrapTrackParam;
 390   AliMUONTrackParam extrapTrackParamAtCh;
 391   AliMUONTrackParam extrapTrackParamAtCluster1;
 392   AliMUONTrackParam extrapTrackParamAtCluster2;
 393   AliMUONTrackParam bestTrackParamAtCluster1;
 394   AliMUONTrackParam bestTrackParamAtCluster2;
 395
 396   Int_t nClusters = clusterStore.GetSize();
 397   Bool_t *clusterCh1Used = new Bool_t[nClusters];
 398   for (Int_t i = 0; i < nClusters; i++) clusterCh1Used[i] = kFALSE;
 399   Int_t iCluster1;
 400
 401   // Get track parameters according to the propagation direction
 402
 403   if (nextStation==4) extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last();
 404   else extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First();
 405
 406   // Add MCS effect
 407   AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
 408
 409   // reset propagator for smoother
 410   if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) extrapTrackParamAtCh.ResetPropagator();
 411
 412   // Add MCS in the missing chamber if any (only 1 chamber can be missing according to tracking criteria)
 413   if (ch1 < ch2 && extrapTrackParamAtCh.GetClusterPtr()->GetChamberId() > ch2 + 1) {
 414     // extrapolation to the missing chamber
 415     AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2 + 1),
 416                                      AliMUONReconstructor::GetRecoParam()->UseSmoother());
 417     // add MCS effect
 418     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
 419   }
 420
 421   //Extrapolate trackCandidate to chamber "ch2"
 422   AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2),
 423                                    AliMUONReconstructor::GetRecoParam()->UseSmoother());
 424
 425   // Printout for debuging
 426   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 427     cout<<endl<<"Track parameter covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch2)<<":"<<endl;
 428     extrapTrackParamAtCh.GetCovariances().Print();
 429   }
 430
 431   // Printout for debuging
 432   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 433     cout << "FollowTrackInStation: look for clusters in chamber(1..): " << ch2+1 << endl;
 434   }
 435
 436   // Ask the clustering to reconstruct new clusters around the track position in the current station
 437   // except for station 4 and 5 that are already entirely clusterized
 438   if (AliMUONReconstructor::GetRecoParam()->CombineClusterTrackReco()) {
 439     if (nextStation < 3) AskForNewClustersInStation(extrapTrackParamAtCh, clusterStore, nextStation);
 440   }
 441
 442   // Create iterators to loop over clusters in both chambers
 443   TIter nextInCh1(clusterStore.CreateChamberIterator(ch1,ch1));
 444   TIter nextInCh2(clusterStore.CreateChamberIterator(ch2,ch2));
 445
 446   // look for candidates in chamber 2
 447   while ( ( clusterCh2 = static_cast<AliMUONVCluster*>(nextInCh2()) ) ) {
 448
 449     // try to add the current cluster fast
 450     if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh2)) continue;
 451
 452     // try to add the current cluster accuratly
 453     chi2OfCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh2, extrapTrackParamAtCluster2,
 454                                   AliMUONReconstructor::GetRecoParam()->UseSmoother());
 455
 456     // if good chi2 then try to attach a cluster in the other chamber too
 457     if (chi2OfCluster < maxChi2OfCluster) {
 458       Bool_t foundSecondCluster = kFALSE;
 459
 460       // Printout for debuging
 461       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 462         cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch2+1
 463              << " (Chi2 = " << chi2OfCluster << ")" << endl;
 464         cout << "                      look for second clusters in chamber(1..): " << ch1+1 << " ..." << endl;
 465       }
 466
 467       if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) {
 468         // save extrapolated parameters for smoother
 469         extrapTrackParamAtCluster2.SetExtrapParameters(extrapTrackParamAtCluster2.GetParameters());
 470
 471         // save extrapolated covariance matrix for smoother
 472         extrapTrackParamAtCluster2.SetExtrapCovariances(extrapTrackParamAtCluster2.GetCovariances());
 473       }
 474
 475       // Compute new track parameters including "clusterCh2" using kalman filter
 476       addChi2TrackAtCluster2 = RunKalmanFilter(extrapTrackParamAtCluster2);
 477
 478       // copy new track parameters for next step
 479       extrapTrackParam = extrapTrackParamAtCluster2;
 480
 481       // add MCS effect
 482       AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
 483
 484       // reset propagator for smoother
 485       if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) extrapTrackParam.ResetPropagator();
 486
 487       //Extrapolate track parameters to chamber "ch1"
 488       AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(ch1),
 489                                        AliMUONReconstructor::GetRecoParam()->UseSmoother());
 490
 491       // reset cluster iterator of chamber 1
 492       nextInCh1.Reset();
 493       iCluster1 = -1;
 494
 495       // look for second candidates in chamber 1
 496       while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
 497         iCluster1++;
 498
 499         // try to add the current cluster fast
 500         if (!TryOneClusterFast(extrapTrackParam, clusterCh1)) continue;
 501
 502         // try to add the current cluster accuratly
 503         chi2OfCluster = TryOneCluster(extrapTrackParam, clusterCh1, extrapTrackParamAtCluster1,
 504                                           AliMUONReconstructor::GetRecoParam()->UseSmoother());
 505
 506         // if good chi2 then consider to add the 2 clusters to the "trackCandidate"
 507         if (chi2OfCluster < maxChi2OfCluster) {
 508           foundSecondCluster = kTRUE;
 509           foundTwoClusters = kTRUE;
 510
 511           // Printout for debuging
 512           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 513             cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch1+1
 514                  << " (Chi2 = " << chi2OfCluster << ")" << endl;
 515           }
 516
 517           if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) {
 518             // save extrapolated parameters for smoother
 519             extrapTrackParamAtCluster1.SetExtrapParameters(extrapTrackParamAtCluster1.GetParameters());
 520
 521             // save extrapolated covariance matrix for smoother
 522             extrapTrackParamAtCluster1.SetExtrapCovariances(extrapTrackParamAtCluster1.GetCovariances());
 523           }
 524
 525           // Compute new track parameters including "clusterCh1" using kalman filter
 526           addChi2TrackAtCluster1 = RunKalmanFilter(extrapTrackParamAtCluster1);
 527
 528           if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
 529             // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new clusters
 530             newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 531             UpdateTrack(*newTrack,extrapTrackParamAtCluster1,extrapTrackParamAtCluster2,addChi2TrackAtCluster1,addChi2TrackAtCluster2);
 532             fNRecTracks++;
 533
 534             // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
 535             // (going in the right direction)
 536             if (nextStation == 4) RetraceTrack(*newTrack,kTRUE);
 537
 538             // Tag clusterCh1 as used
 539             clusterCh1Used[iCluster1] = kTRUE;
 540
 541             // Printout for debuging
 542             if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 543               cout << "FollowTrackInStation: added two clusters in station(1..): " << nextStation+1 << endl;
 544               if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 545             }
 546
 547           } else if (addChi2TrackAtCluster1+addChi2TrackAtCluster2 < bestAddChi2TrackAtCluster1+bestAddChi2TrackAtCluster2) {
 548             // keep track of the best couple of clusters
 549             bestAddChi2TrackAtCluster1 = addChi2TrackAtCluster1;
 550             bestAddChi2TrackAtCluster2 = addChi2TrackAtCluster2;
 551             bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
 552             bestTrackParamAtCluster2 = extrapTrackParamAtCluster2;
 553           }
 554
 555         }
 556
 557       }
 558
 559       // if no clusterCh1 found then consider to add clusterCh2 only
 560       if (!foundSecondCluster) {
 561         foundOneCluster = kTRUE;
 562
 563         if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
 564           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
 565           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 566           UpdateTrack(*newTrack,extrapTrackParamAtCluster2,addChi2TrackAtCluster2);
 567           fNRecTracks++;
 568
 569           // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
 570           // (going in the right direction)
 571           if (nextStation == 4) RetraceTrack(*newTrack,kTRUE);
 572
 573           // Printout for debuging
 574           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 575             cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch2+1 << endl;
 576             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 577           }
 578
 579         } else if (!foundTwoClusters && addChi2TrackAtCluster2 < bestAddChi2TrackAtCluster1) {
 580           // keep track of the best single cluster except if a couple of clusters has already been found
 581           bestAddChi2TrackAtCluster1 = addChi2TrackAtCluster2;
 582           bestTrackParamAtCluster1 = extrapTrackParamAtCluster2;
 583         }
 584
 585       }
 586
 587     }
 588
 589   }
 590
 591   // look for candidates in chamber 1 not already attached to a track
 592   // if we want to keep all possible tracks or if no good couple of clusters has been found
 593   if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks() || !foundTwoClusters) {
 594
 595     // Printout for debuging
 596     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 597       cout << "FollowTrackInStation: look for single clusters in chamber(1..): " << ch1+1 << endl;
 598     }
 599
 600     // add MCS effect for next step
 601     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
 602
 603     //Extrapolate trackCandidate to chamber "ch1"
 604     AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch1),
 605                                      AliMUONReconstructor::GetRecoParam()->UseSmoother());
 606
 607     // reset cluster iterator of chamber 1
 608     nextInCh1.Reset();
 609     iCluster1 = -1;
 610
 611     // look for second candidates in chamber 1
 612     while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
 613       iCluster1++;
 614
 615       if (clusterCh1Used[iCluster1]) continue; // Skip clusters already used
 616
 617       // try to add the current cluster fast
 618       if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh1)) continue;
 619
 620       // try to add the current cluster accuratly
 621       chi2OfCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh1, extrapTrackParamAtCluster1,
 622                                     AliMUONReconstructor::GetRecoParam()->UseSmoother());
 623
 624       // if good chi2 then consider to add clusterCh1
 625       // We do not try to attach a cluster in the other chamber too since it has already been done above
 626       if (chi2OfCluster < maxChi2OfCluster) {
 627         foundOneCluster = kTRUE;
 628
 629         // Printout for debuging
 630         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 631           cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch1+1
 632                << " (Chi2 = " << chi2OfCluster << ")" << endl;
 633         }
 634
 635         if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) {
 636           // save extrapolated parameters for smoother
 637           extrapTrackParamAtCluster1.SetExtrapParameters(extrapTrackParamAtCluster1.GetParameters());
 638
 639           // save extrapolated covariance matrix for smoother
 640           extrapTrackParamAtCluster1.SetExtrapCovariances(extrapTrackParamAtCluster1.GetCovariances());
 641         }
 642
 643         // Compute new track parameters including "clusterCh1" using kalman filter
 644         addChi2TrackAtCluster1 = RunKalmanFilter(extrapTrackParamAtCluster1);
 645
 646         if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
 647           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
 648           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 649           UpdateTrack(*newTrack,extrapTrackParamAtCluster1,addChi2TrackAtCluster1);
 650           fNRecTracks++;
 651
 652           // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
 653           // (going in the right direction)
 654           if (nextStation == 4) RetraceTrack(*newTrack,kTRUE);
 655
 656           // Printout for debuging
 657           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 658             cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch1+1 << endl;
 659             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 660           }
 661
 662         } else if (addChi2TrackAtCluster1 < bestAddChi2TrackAtCluster1) {
 663           // keep track of the best single cluster except if a couple of clusters has already been found
 664           bestAddChi2TrackAtCluster1 = addChi2TrackAtCluster1;
 665           bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
 666         }
 667
 668       }
 669
 670     }
 671
 672   }
 673
 674   // fill out the best track if required else clean up the fRecTracksPtr array
 675   if (!AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
 676     if (foundTwoClusters) {
 677       UpdateTrack(trackCandidate,bestTrackParamAtCluster1,bestTrackParamAtCluster2,bestAddChi2TrackAtCluster1,bestAddChi2TrackAtCluster2);
 678
 679       // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
 680       // (going in the right direction)
 681       if (nextStation == 4) RetraceTrack(trackCandidate,kTRUE);
 682
 683       // Printout for debuging
 684       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 685         cout << "FollowTrackInStation: added the two best clusters in station(1..): " << nextStation+1 << endl;
 686         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 687       }
 688
 689     } else if (foundOneCluster) {
 690       UpdateTrack(trackCandidate,bestTrackParamAtCluster1,bestAddChi2TrackAtCluster1);
 691
 692       // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
 693       // (going in the right direction)
 694       if (nextStation == 4) RetraceTrack(trackCandidate,kTRUE);
 695
 696       // Printout for debuging
 697       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 698         cout << "FollowTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster1.GetClusterPtr()->GetChamberId()+1 << endl;
 699         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 700       }
 701
 702     } else {
 703       delete [] clusterCh1Used;
 704       return kFALSE;
 705     }
 706   } else if (foundOneCluster || foundTwoClusters) {
 707
 708     // remove obsolete track
 709     fRecTracksPtr->Remove(&trackCandidate);
 710     fNRecTracks--;
 711
 712   } else {
 713     delete [] clusterCh1Used;
 714     return kFALSE;
 715   }
 716   delete [] clusterCh1Used;
 717   return kTRUE;
 718
 719 }
 720
 721   //__________________________________________________________________________
 722 Double_t AliMUONTrackReconstructorK::RunKalmanFilter(AliMUONTrackParam &trackParamAtCluster)
 723 {
 724   /// Compute new track parameters and their covariances including new cluster using kalman filter
 725   /// return the additional track chi2
 726   AliDebug(1,"Enter RunKalmanFilter");
 727
 728   // Get actual track parameters (p)
 729   TMatrixD param(trackParamAtCluster.GetParameters());
 730
 731   // Get new cluster parameters (m)
 732   AliMUONVCluster *cluster = trackParamAtCluster.GetClusterPtr();
 733   TMatrixD clusterParam(5,1);
 734   clusterParam.Zero();
 735   clusterParam(0,0) = cluster->GetX();
 736   clusterParam(2,0) = cluster->GetY();
 737
 738   // Compute the actual parameter weight (W)
 739   TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
 740   if (paramWeight.Determinant() != 0) {
 741     paramWeight.Invert();
 742   } else {
 743     AliWarning(" Determinant = 0");
 744     return 1.e10;
 745   }
 746
 747   // Compute the new cluster weight (U)
 748   TMatrixD clusterWeight(5,5);
 749   clusterWeight.Zero();
 750   clusterWeight(0,0) = 1. / cluster->GetErrX2();
 751   clusterWeight(2,2) = 1. / cluster->GetErrY2();
 752
 753   // Compute the new parameters covariance matrix ( (W+U)^-1 )
 754   TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
 755   if (newParamCov.Determinant() != 0) {
 756     newParamCov.Invert();
 757   } else {
 758     AliWarning(" Determinant = 0");
 759     return 1.e10;
 760   }
 761
 762   // Save the new parameters covariance matrix
 763   trackParamAtCluster.SetCovariances(newParamCov);
 764
 765   // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
 766   TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
 767   TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp); // U(m-p)
 768   TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
 769   newParam += param; // ((W+U)^-1)U(m-p) + p
 770
 771   // Save the new parameters
 772   trackParamAtCluster.SetParameters(newParam);
 773
 774   // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
 775   tmp = newParam; // p'
 776   tmp -= param; // (p'-p)
 777   TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp); // W(p'-p)
 778   TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
 779   tmp = newParam; // p'
 780   tmp -= clusterParam; // (p'-m)
 781   TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp); // U(p'-m)
 782   addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
 783
 784   return addChi2Track(0,0);
 785
 786 }
 787
 788   //__________________________________________________________________________
 789 void AliMUONTrackReconstructorK::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster, Double_t addChi2)
 790 {
 791   /// Add 1 cluster to the track candidate
 792   /// Update chi2 of the track
 793
 794   // Flag cluster as being not removable
 795   trackParamAtCluster.SetRemovable(kFALSE);
 796   trackParamAtCluster.SetLocalChi2(0.); // --> Local chi2 not used
 797
 798   // Update the track chi2 into trackParamAtCluster
 799   trackParamAtCluster.SetTrackChi2(track.GetGlobalChi2() + addChi2);
 800
 801   // Update the chi2 of the new track
 802   track.SetGlobalChi2(trackParamAtCluster.GetTrackChi2());
 803
 804   // Update array of TrackParamAtCluster
 805   track.AddTrackParamAtCluster(trackParamAtCluster,*(trackParamAtCluster.GetClusterPtr()));
 806
 807 }
 808
 809   //__________________________________________________________________________
 810 void AliMUONTrackReconstructorK::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster1, AliMUONTrackParam &trackParamAtCluster2,
 811                                              Double_t addChi2AtCluster1, Double_t addChi2AtCluster2)
 812 {
 813   /// Add 2 clusters to the track candidate (order is important)
 814   /// Update track and local chi2
 815
 816   // Update local chi2 at first cluster
 817   AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr();
 818   Double_t deltaX = trackParamAtCluster1.GetNonBendingCoor() - cluster1->GetX();
 819   Double_t deltaY = trackParamAtCluster1.GetBendingCoor() - cluster1->GetY();
 820   Double_t localChi2AtCluster1 = deltaX*deltaX / cluster1->GetErrX2() +
 821                                  deltaY*deltaY / cluster1->GetErrY2();
 822   trackParamAtCluster1.SetLocalChi2(localChi2AtCluster1);
 823
 824   // Flag first cluster as being removable
 825   trackParamAtCluster1.SetRemovable(kTRUE);
 826
 827   // Update local chi2 at second cluster
 828   AliMUONVCluster* cluster2 = trackParamAtCluster2.GetClusterPtr();
 829   AliMUONTrackParam extrapTrackParamAtCluster2(trackParamAtCluster1);
 830   AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParamAtCluster2, trackParamAtCluster2.GetZ());
 831   deltaX = extrapTrackParamAtCluster2.GetNonBendingCoor() - cluster2->GetX();
 832   deltaY = extrapTrackParamAtCluster2.GetBendingCoor() - cluster2->GetY();
 833   Double_t localChi2AtCluster2 = deltaX*deltaX / cluster2->GetErrX2() +
 834                                  deltaY*deltaY / cluster2->GetErrY2();
 835   trackParamAtCluster2.SetLocalChi2(localChi2AtCluster2);
 836
 837   // Flag second cluster as being removable
 838   trackParamAtCluster2.SetRemovable(kTRUE);
 839
 840   // Update the track chi2 into trackParamAtCluster2
 841   trackParamAtCluster2.SetTrackChi2(track.GetGlobalChi2() + addChi2AtCluster2);
 842
 843   // Update the track chi2 into trackParamAtCluster1
 844   trackParamAtCluster1.SetTrackChi2(trackParamAtCluster2.GetTrackChi2() + addChi2AtCluster1);
 845
 846   // Update the chi2 of the new track
 847   track.SetGlobalChi2(trackParamAtCluster1.GetTrackChi2());
 848
 849   // Update array of trackParamAtCluster
 850   track.AddTrackParamAtCluster(trackParamAtCluster1,*cluster1);
 851   track.AddTrackParamAtCluster(trackParamAtCluster2,*cluster2);
 852
 853 }
 854
 855   //__________________________________________________________________________
 856 Bool_t AliMUONTrackReconstructorK::RecoverTrack(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation)
 857 {
 858   /// Try to recover the track candidate in the next station
 859   /// by removing the worst of the two clusters attached in the current station
 860   /// Return kTRUE if recovering succeeds
 861   AliDebug(1,"Enter RecoverTrack");
 862
 863   // Do not try to recover track until we have attached cluster(s) on station(1..) 3
 864   if (nextStation > 1) return kFALSE;
 865
 866   Int_t worstClusterNumber = -1;
 867   Double_t localChi2, worstLocalChi2 = 0.;
 868
 869   // Look for the cluster to remove
 870   for (Int_t clusterNumber = 0; clusterNumber < 2; clusterNumber++) {
 871     AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(clusterNumber);
 872
 873     // check if current cluster is removable
 874     if (!trackParamAtCluster->IsRemovable()) return kFALSE;
 875
 876     // Pick up cluster with the worst chi2
 877     localChi2 = trackParamAtCluster->GetLocalChi2();
 878     if (localChi2 > worstLocalChi2) {
 879       worstLocalChi2 = localChi2;
 880       worstClusterNumber = clusterNumber;
 881     }
 882   }
 883
 884   // Reset best cluster as being NOT removable
 885   ((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt((worstClusterNumber+1)%2))->SetRemovable(kFALSE);
 886
 887   // Remove the worst cluster
 888   trackCandidate.RemoveTrackParamAtCluster((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(worstClusterNumber));
 889
 890   // Re-calculate track parameters at the (new) first cluster
 891   RetracePartialTrack(trackCandidate,(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(1));
 892
 893   // Look for new cluster(s) in next station
 894   return FollowTrackInStation(trackCandidate, clusterStore, nextStation);
 895
 896 }
 897
 898   //__________________________________________________________________________
 899 Bool_t AliMUONTrackReconstructorK::RunSmoother(AliMUONTrack &track)
 900 {
 901   /// Compute new track parameters and their covariances using smoother
 902   AliDebug(1,"Enter UseSmoother");
 903
 904   AliMUONTrackParam *previousTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->First();
 905
 906   // Smoothed parameters and covariances at first cluster = filtered parameters and covariances
 907   previousTrackParam->SetSmoothParameters(previousTrackParam->GetParameters());
 908   previousTrackParam->SetSmoothCovariances(previousTrackParam->GetCovariances());
 909
 910   // Compute local chi2 at first cluster
 911   AliMUONVCluster *cluster = previousTrackParam->GetClusterPtr();
 912   Double_t dX = cluster->GetX() - previousTrackParam->GetNonBendingCoor();
 913   Double_t dY = cluster->GetY() - previousTrackParam->GetBendingCoor();
 914   Double_t localChi2 = dX * dX / cluster->GetErrX2() + dY * dY / cluster->GetErrY2();
 915
 916   // Save local chi2 at first cluster
 917   previousTrackParam->SetLocalChi2(localChi2);
 918
 919   AliMUONTrackParam *currentTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(previousTrackParam);
 920   while (currentTrackParam) {
 921
 922     // Get variables
 923     const TMatrixD &extrapParameters          = previousTrackParam->GetExtrapParameters();  // X(k+1 k)
 924     const TMatrixD &filteredParameters        = currentTrackParam->GetParameters();         // X(k k)
 925     const TMatrixD &previousSmoothParameters  = previousTrackParam->GetSmoothParameters();  // X(k+1 n)
 926     const TMatrixD &propagator                = previousTrackParam->GetPropagator();        // F(k)
 927     const TMatrixD &extrapCovariances         = previousTrackParam->GetExtrapCovariances(); // C(k+1 k)
 928     const TMatrixD &filteredCovariances       = currentTrackParam->GetCovariances();        // C(k k)
 929     const TMatrixD &previousSmoothCovariances = previousTrackParam->GetSmoothCovariances(); // C(k+1 n)
 930
 931     // Compute smoother gain: A(k) = C(kk) * F(f)^t * (C(k+1 k))^-1
 932     TMatrixD extrapWeight(extrapCovariances);
 933     if (extrapWeight.Determinant() != 0) {
 934       extrapWeight.Invert(); // (C(k+1 k))^-1
 935     } else {
 936       AliWarning(" Determinant = 0");
 937       return kFALSE;
 938     }
 939     TMatrixD smootherGain(filteredCovariances,TMatrixD::kMultTranspose,propagator); // C(kk) * F(f)^t
 940     smootherGain *= extrapWeight; // C(kk) * F(f)^t * (C(k+1 k))^-1
 941
 942     // Compute smoothed parameters: X(k n) = X(k k) + A(k) * (X(k+1 n) - X(k+1 k))
 943     TMatrixD tmpParam(previousSmoothParameters,TMatrixD::kMinus,extrapParameters); // X(k+1 n) - X(k+1 k)
 944     TMatrixD smoothParameters(smootherGain,TMatrixD::kMult,tmpParam); // A(k) * (X(k+1 n) - X(k+1 k))
 945     smoothParameters += filteredParameters; // X(k k) + A(k) * (X(k+1 n) - X(k+1 k))
 946
 947     // Save smoothed parameters
 948     currentTrackParam->SetSmoothParameters(smoothParameters);
 949
 950     // Compute smoothed covariances: C(k n) = C(k k) + A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
 951     TMatrixD tmpCov(previousSmoothCovariances,TMatrixD::kMinus,extrapCovariances); // C(k+1 n) - C(k+1 k)
 952     TMatrixD tmpCov2(tmpCov,TMatrixD::kMultTranspose,smootherGain); // (C(k+1 n) - C(k+1 k)) * (A(k))^t
 953     TMatrixD smoothCovariances(smootherGain,TMatrixD::kMult,tmpCov2); // A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
 954     smoothCovariances += filteredCovariances; // C(k k) + A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
 955
 956     // Save smoothed covariances
 957     currentTrackParam->SetSmoothCovariances(smoothCovariances);
 958
 959     // Compute smoothed residual: r(k n) = cluster - X(k n)
 960     cluster = currentTrackParam->GetClusterPtr();
 961     TMatrixD smoothResidual(2,1);
 962     smoothResidual.Zero();
 963     smoothResidual(0,0) = cluster->GetX() - smoothParameters(0,0);
 964     smoothResidual(1,0) = cluster->GetY() - smoothParameters(2,0);
 965
 966     // Compute weight of smoothed residual: W(k n) = (clusterCov - C(k n))^-1
 967     TMatrixD smoothResidualWeight(2,2);
 968     smoothResidualWeight(0,0) = cluster->GetErrX2() - smoothCovariances(0,0);
 969     smoothResidualWeight(0,1) = - smoothCovariances(0,2);
 970     smoothResidualWeight(1,0) = - smoothCovariances(2,0);
 971     smoothResidualWeight(1,1) = cluster->GetErrY2() - smoothCovariances(2,2);
 972     if (smoothResidualWeight.Determinant() != 0) {
 973       smoothResidualWeight.Invert();
 974     } else {
 975       AliWarning(" Determinant = 0");
 976       return kFALSE;
 977     }
 978
 979     // Compute local chi2 = (r(k n))^t * W(k n) * r(k n)
 980     TMatrixD tmpChi2(smoothResidual,TMatrixD::kTransposeMult,smoothResidualWeight); // (r(k n))^t * W(k n)
 981     TMatrixD localChi2(tmpChi2,TMatrixD::kMult,smoothResidual); // (r(k n))^t * W(k n) * r(k n)
 982
 983     // Save local chi2
 984     currentTrackParam->SetLocalChi2(localChi2(0,0));
 985
 986     previousTrackParam = currentTrackParam;
 987     currentTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(previousTrackParam);
 988   }
 989
 990   return kTRUE;
 991
 992 }
 993
 994   //__________________________________________________________________________
 995 void AliMUONTrackReconstructorK::ComplementTracks(const AliMUONVClusterStore& clusterStore)
 996 {
 997   /// Complete tracks by adding missing clusters (if there is an overlap between
 998   /// two detection elements, the track may have two clusters in the same chamber)
 999   /// Recompute track parameters and covariances at each clusters
1000   AliDebug(1,"Enter ComplementTracks");
1001
1002   Int_t chamberId, detElemId;
1003   Double_t chi2OfCluster, addChi2TrackAtCluster, bestAddChi2TrackAtCluster;
1004   Double_t maxChi2OfCluster = 2. * AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
1005                                    AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
1006   Bool_t foundOneCluster, trackModified;
1007   AliMUONVCluster *cluster;
1008   AliMUONTrackParam *trackParam, *previousTrackParam, *nextTrackParam, trackParamAtCluster, bestTrackParamAtCluster;
1009
1010   // Remove double track to complete only "good" tracks
1011   RemoveDoubleTracks();
1012
1013   AliMUONTrack *track = (AliMUONTrack*) fRecTracksPtr->First();
1014   while (track) {
1015     trackModified = kFALSE;
1016
1017     trackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->First();
1018     previousTrackParam = trackParam;
1019     while (trackParam) {
1020       foundOneCluster = kFALSE;
1021       bestAddChi2TrackAtCluster = 1.e10;
1022       chamberId = trackParam->GetClusterPtr()->GetChamberId();
1023       detElemId = trackParam->GetClusterPtr()->GetDetElemId();
1024
1025       // prepare nextTrackParam before adding new cluster because of the sorting
1026       nextTrackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->After(trackParam);
1027
1028       // Create iterators to loop over clusters in current chamber
1029       TIter nextInCh(clusterStore.CreateChamberIterator(chamberId,chamberId));
1030
1031       // look for one second candidate in the same chamber
1032       while ( ( cluster = static_cast<AliMUONVCluster*>(nextInCh()) ) ) {
1033
1034         // look for a cluster in another detection element
1035         if (cluster->GetDetElemId() == detElemId) continue;
1036
1037         // try to add the current cluster fast
1038         if (!TryOneClusterFast(*trackParam, cluster)) continue;
1039
1040         // try to add the current cluster accurately
1041         // never use track parameters at last cluster because the covariance matrix is meaningless
1042         if (nextTrackParam) chi2OfCluster = TryOneCluster(*trackParam, cluster, trackParamAtCluster);
1043         else chi2OfCluster = TryOneCluster(*previousTrackParam, cluster, trackParamAtCluster);
1044
1045         // if good chi2 then consider to add this cluster to the track
1046         if (chi2OfCluster < maxChi2OfCluster) {
1047
1048           // Compute local track parameters including current cluster using kalman filter
1049           addChi2TrackAtCluster = RunKalmanFilter(trackParamAtCluster);
1050
1051           // keep track of the best cluster
1052           if (addChi2TrackAtCluster < bestAddChi2TrackAtCluster) {
1053             bestAddChi2TrackAtCluster = addChi2TrackAtCluster;
1054             bestTrackParamAtCluster = trackParamAtCluster;
1055             foundOneCluster = kTRUE;
1056           }
1057
1058         }
1059
1060       }
1061
1062       // add new cluster if any
1063       if (foundOneCluster) {
1064         trackParam->SetRemovable(kTRUE);
1065         bestTrackParamAtCluster.SetRemovable(kTRUE);
1066         track->AddTrackParamAtCluster(bestTrackParamAtCluster,*(bestTrackParamAtCluster.GetClusterPtr()));
1067         trackModified = kTRUE;
1068       }
1069
1070       previousTrackParam = trackParam;
1071       trackParam = nextTrackParam;
1072     }
1073
1074     // re-compute track parameters using kalman filter if needed
1075     if (trackModified) RetraceTrack(*track,kTRUE);
1076
1077     track = (AliMUONTrack*) fRecTracksPtr->After(track);
1078   }
1079
1080 }
1081
1082 //__________________________________________________________________________
1083 void AliMUONTrackReconstructorK::ImproveTrack(AliMUONTrack &track)
1084 {
1085   /// Improve the given track by removing removable clusters with local chi2 highter than the defined cut
1086   /// Removable clusters are identified by the method AliMUONTrack::TagRemovableClusters()
1087   /// Recompute track parameters and covariances at the remaining clusters
1088   AliDebug(1,"Enter ImproveTrack");
1089
1090   Double_t localChi2, worstLocalChi2;
1091   AliMUONTrackParam *trackParamAtCluster, *worstTrackParamAtCluster, *nextTrackParam;
1092   Bool_t smoothed;
1093   Double_t sigmaCut2 = AliMUONReconstructor::GetRecoParam()->GetSigmaCutForImprovement() *
1094                        AliMUONReconstructor::GetRecoParam()->GetSigmaCutForImprovement();
1095
1096   while (!track.IsImproved()) {
1097
1098     // identify removable clusters
1099     track.TagRemovableClusters();
1100
1101     // Run smoother if required
1102     smoothed = kFALSE;
1103     if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) smoothed = RunSmoother(track);
1104
1105     // Use standard procedure to compute local chi2 if smoother not required or not working
1106     if (!smoothed) {
1107
1108       // Update track parameters and covariances
1109       track.UpdateCovTrackParamAtCluster();
1110
1111       // Compute local chi2 of each clusters
1112       track.ComputeLocalChi2(kTRUE);
1113     }
1114
1115     // Look for the cluster to remove
1116     worstTrackParamAtCluster = 0x0;
1117     worstLocalChi2 = 0.;
1118     trackParamAtCluster = (AliMUONTrackParam*)track.GetTrackParamAtCluster()->First();
1119     while (trackParamAtCluster) {
1120
1121       // save parameters into smooth parameters in case of smoother did not work properly
1122       if (AliMUONReconstructor::GetRecoParam()->UseSmoother() && !smoothed) {
1123         trackParamAtCluster->SetSmoothParameters(trackParamAtCluster->GetParameters());
1124         trackParamAtCluster->SetSmoothCovariances(trackParamAtCluster->GetCovariances());
1125       }
1126
1127       // Pick up cluster with the worst chi2
1128       localChi2 = trackParamAtCluster->GetLocalChi2();
1129       if (localChi2 > worstLocalChi2) {
1130         worstLocalChi2 = localChi2;
1131         worstTrackParamAtCluster = trackParamAtCluster;
1132       }
1133
1134       trackParamAtCluster = (AliMUONTrackParam*)track.GetTrackParamAtCluster()->After(trackParamAtCluster);
1135     }
1136
1137     // Check if worst cluster found
1138     if (!worstTrackParamAtCluster) {
1139       AliWarning("Bad local chi2 values?");
1140       break;
1141     }
1142
1143     // Check whether the worst chi2 is under requirement or not
1144     if (worstLocalChi2 < 2. * sigmaCut2) { // 2 because 2 quantities in chi2
1145       track.SetImproved(kTRUE);
1146       break;
1147     }
1148
1149     // if the worst cluster is not removable then stop improvement
1150     if (!worstTrackParamAtCluster->IsRemovable()) break;
1151
1152     // get track parameters at cluster next to the one to be removed
1153     nextTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(worstTrackParamAtCluster);
1154
1155     // Remove the worst cluster
1156     track.RemoveTrackParamAtCluster(worstTrackParamAtCluster);
1157
1158     // Re-calculate track parameters
1159     // - from the cluster immediately downstream the one suppressed
1160     // - or from the begining - if parameters have been re-computed using the standard method (kalman parameters have been lost)
1161     //                  - or if the removed cluster was the last one
1162     if (smoothed && nextTrackParam) RetracePartialTrack(track,nextTrackParam);
1163     else RetraceTrack(track,kTRUE);
1164
1165     // Printout for debuging
1166     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1167       cout << "ImproveTracks: track " << fRecTracksPtr->IndexOf(&track)+1 << " improved " << endl;
1168     }
1169
1170   }
1171
1172 }
1173
1174 //__________________________________________________________________________
1175 void AliMUONTrackReconstructorK::FinalizeTrack(AliMUONTrack &track)
1176 {
1177   /// Update track parameters and covariances at each attached cluster
1178   /// using smoother if required, if not already done
1179
1180   AliMUONTrackParam *trackParamAtCluster;
1181   Bool_t smoothed = kFALSE;
1182
1183   // update track parameters (using smoother if required) if not already done
1184   if (!track.IsImproved()) {
1185     smoothed = kFALSE;
1186     if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) smoothed = RunSmoother(track);
1187     if (!smoothed) track.UpdateCovTrackParamAtCluster();
1188   } else smoothed = AliMUONReconstructor::GetRecoParam()->UseSmoother();
1189
1190   // copy smoothed parameters and covariances if any
1191   if (smoothed) {
1192
1193     trackParamAtCluster = (AliMUONTrackParam*) (track.GetTrackParamAtCluster()->First());
1194     while (trackParamAtCluster) {
1195
1196       trackParamAtCluster->SetParameters(trackParamAtCluster->GetSmoothParameters());
1197       trackParamAtCluster->SetCovariances(trackParamAtCluster->GetSmoothCovariances());
1198
1199       trackParamAtCluster = (AliMUONTrackParam*) (track.GetTrackParamAtCluster()->After(trackParamAtCluster));
1200     }
1201
1202   }
1203
1204 }
1205
1206   //__________________________________________________________________________
1207 Bool_t AliMUONTrackReconstructorK::RefitTrack(AliMUONTrack &track)
1208 {
1209   /// re-fit the given track
1210
1211   // check validity of the track
1212   if (!track.IsValid()) {
1213     AliWarning("the track does not contain enough clusters --> unable to refit");
1214     return kFALSE;
1215   }
1216
1217   // re-compute track parameters and covariances using Kalman filter
1218   RetraceTrack(track,kTRUE);
1219
1220   // Improve the reconstructed tracks if required
1221   if (AliMUONReconstructor::GetRecoParam()->ImproveTracks()) ImproveTrack(track);
1222
1223   // Fill AliMUONTrack data members
1224   FinalizeTrack(track);
1225
1226   return kTRUE;
1227
1228 }
1229