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