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