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