MUON/AliMUONTrackReconstructor.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 AliMUONTrackReconstructor
  20 /// MUON track reconstructor using the original method
  21 ///
  22 /// This class contains as data:
  23 /// - the parameters for the track reconstruction
  24 ///
  25 /// It contains as methods, among others:
  26 /// - MakeTracks to build the tracks
  27 //-----------------------------------------------------------------------------
  28
  29 #include "AliMUONTrackReconstructor.h"
  30
  31 #include "AliMUONConstants.h"
  32 #include "AliMUONVCluster.h"
  33 #include "AliMUONVClusterServer.h"
  34 #include "AliMUONVClusterStore.h"
  35 #include "AliMUONTrack.h"
  36 #include "AliMUONTrackParam.h"
  37 #include "AliMUONTrackExtrap.h"
  38 #include "AliMUONRecoParam.h"
  39
  40 #include "AliMpArea.h"
  41
  42 #include "AliLog.h"
  43
  44 #include <TMinuit.h>
  45 #include <Riostream.h>
  46 #include <TMath.h>
  47 #include <TMatrixD.h>
  48 #include <TClonesArray.h>
  49
  50 // Functions to be minimized with Minuit
  51 void TrackChi2(Int_t &nParam, Double_t *gradient, Double_t &chi2, Double_t *param, Int_t flag);
  52
  53 /// \cond CLASSIMP
  54 ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context
  55 /// \endcond
  56
  57   //__________________________________________________________________________
  58 AliMUONTrackReconstructor::AliMUONTrackReconstructor(const AliMUONRecoParam* recoParam, AliMUONVClusterServer* clusterServer)
  59   : AliMUONVTrackReconstructor(recoParam,clusterServer)
  60 {
  61   /// Constructor
  62 }
  63
  64   //__________________________________________________________________________
  65 AliMUONTrackReconstructor::~AliMUONTrackReconstructor()
  66 {
  67 /// Destructor
  68 }
  69
  70   //__________________________________________________________________________
  71 Bool_t AliMUONTrackReconstructor::MakeTrackCandidates(AliMUONVClusterStore& clusterStore)
  72 {
  73   /// To make track candidates (assuming linear propagation if the flag fgkMakeTrackCandidatesFast is set to kTRUE):
  74   /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4.
  75   /// Good candidates are made of at least three clusters.
  76   /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
  77
  78   TClonesArray *segments;
  79   AliMUONTrack *track;
  80   Int_t iCandidate = 0;
  81   Bool_t clusterFound;
  82
  83   AliDebug(1,"Enter MakeTrackCandidates");
  84
  85   // Unless we're doing combined tracking, we'll clusterize all stations at once
  86   Int_t firstChamber(0);
  87   Int_t lastChamber(9);
  88
  89   if (GetRecoParam()->CombineClusterTrackReco()) {
  90     // ... Here's the exception : ask the clustering to reconstruct
  91     // clusters *only* in station 4 and 5 for combined tracking
  92     firstChamber = 6;
  93   }
  94
  95   for (Int_t i = firstChamber; i <= lastChamber; ++i )
  96   {
  97     if (fClusterServer && GetRecoParam()->UseChamber(i)) fClusterServer->Clusterize(i, clusterStore, AliMpArea(), GetRecoParam());
  98   }
  99
 100   // Loop over stations(1..) 5 and 4 and make track candidates
 101   for (Int_t istat=4; istat>=3; istat--) {
 102
 103     // Make segments in the station
 104     segments = MakeSegmentsBetweenChambers(clusterStore, 2*istat, 2*istat+1);
 105
 106     // Loop over segments
 107     for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++)
 108     {
 109       AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
 110
 111       // Transform segments to tracks and put them at the end of fRecTracksPtr
 112       track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack((AliMUONObjectPair*)((*segments)[iseg]),GetRecoParam()->GetBendingVertexDispersion());
 113       fNRecTracks++;
 114
 115       // Look for compatible cluster(s) in the other station
 116       if (GetRecoParam()->MakeTrackCandidatesFast())
 117         clusterFound = FollowLinearTrackInStation(*track, clusterStore, 7-istat);
 118       else clusterFound = FollowTrackInStation(*track, clusterStore, 7-istat);
 119
 120       // Remove track if no cluster found on a requested station
 121       // or abort tracking if there are too many candidates
 122       if (GetRecoParam()->RequestStation(7-istat)) {
 123         if (!clusterFound) {
 124           fRecTracksPtr->Remove(track);
 125           fNRecTracks--;
 126         } else if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
 127           AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks));
 128           delete segments;
 129           return kFALSE;
 130         }
 131       } else {
 132         if ((fNRecTracks + segments->GetEntriesFast() - iseg - 1) > GetRecoParam()->GetMaxTrackCandidates()) {
 133           AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks + segments->GetEntriesFast() - iseg - 1));
 134           delete segments;
 135           return kFALSE;
 136         }
 137       }
 138
 139     }
 140
 141     // delete the array of segments
 142     delete segments;
 143   }
 144
 145   // Keep all different tracks or only the best ones as required
 146   if (GetRecoParam()->TrackAllTracks()) RemoveIdenticalTracks();
 147   else RemoveDoubleTracks();
 148
 149   AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
 150
 151   return kTRUE;
 152
 153 }
 154
 155   //__________________________________________________________________________
 156 Bool_t AliMUONTrackReconstructor::MakeMoreTrackCandidates(AliMUONVClusterStore& clusterStore)
 157 {
 158   /// To make extra track candidates (assuming linear propagation if the flag fgkMakeTrackCandidatesFast is set to kTRUE):
 159   /// clustering is supposed to be already done
 160   /// Start with segments made of 1 cluster in each of the stations 4 and 5 then follow track in remaining chambers.
 161   /// Good candidates are made of at least three clusters if both station are requested (two otherwise).
 162   /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
 163
 164   TClonesArray *segments;
 165   AliMUONObjectPair *segment;
 166   AliMUONTrack *track;
 167   Int_t iCandidate = 0, iCurrentTrack, nCurrentTracks;
 168   Bool_t clusterFound;
 169
 170   AliDebug(1,"Enter MakeMoreTrackCandidates");
 171
 172   // Double loop over chambers in stations(1..) 4 and 5 to make track candidates
 173   for (Int_t ich1 = 6; ich1 <= 7; ich1++) {
 174     for (Int_t ich2 = 8; ich2 <= 9; ich2++) {
 175
 176       // Make segments in the station
 177       segments = MakeSegmentsBetweenChambers(clusterStore, ich1, ich2);
 178
 179       /// Remove segments already attached to a track
 180       RemoveUsedSegments(*segments);
 181
 182       // Loop over segments
 183       for (Int_t iSegment=0; iSegment<segments->GetEntriesFast(); iSegment++)
 184       {
 185         AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
 186         segment = (AliMUONObjectPair*) segments->UncheckedAt(iSegment);
 187
 188         // Transform segments to tracks and put them at the end of fRecTracksPtr
 189         iCurrentTrack = fRecTracksPtr->GetLast()+1;
 190         track = new ((*fRecTracksPtr)[iCurrentTrack]) AliMUONTrack(segment,GetRecoParam()->GetBendingVertexDispersion());
 191         fNRecTracks++;
 192
 193         // Look for compatible cluster(s) in the second chamber of station 5
 194         clusterFound = FollowLinearTrackInChamber(*track, clusterStore, 17-ich2);
 195
 196         // skip the original track in case it has been removed
 197         if (GetRecoParam()->TrackAllTracks() && clusterFound) iCurrentTrack++;
 198
 199         // loop over every new tracks
 200         nCurrentTracks = fRecTracksPtr->GetLast()+1;
 201         while (iCurrentTrack < nCurrentTracks) {
 202           track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iCurrentTrack);
 203
 204           // Look for compatible cluster(s) in the second chamber of station 4
 205           FollowLinearTrackInChamber(*track, clusterStore, 13-ich1);
 206
 207           iCurrentTrack++;
 208         }
 209
 210         // abort tracking if there are too many candidates
 211         if ((fNRecTracks + segments->GetEntriesFast() - iSegment - 1) > GetRecoParam()->GetMaxTrackCandidates()) {
 212           AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks + segments->GetEntriesFast() - iSegment - 1));
 213           delete segments;
 214           return kFALSE;
 215         }
 216
 217       }
 218
 219       // delete the array of segments
 220       delete segments;
 221     }
 222   }
 223
 224   // Keep only the best tracks if required
 225   if (!GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
 226   else fRecTracksPtr->Compress();
 227
 228   AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
 229
 230   return kTRUE;
 231
 232 }
 233
 234   //__________________________________________________________________________
 235 Bool_t AliMUONTrackReconstructor::FollowTracks(AliMUONVClusterStore& clusterStore)
 236 {
 237   /// Follow tracks in stations(1..) 3, 2 and 1
 238   AliDebug(1,"Enter FollowTracks");
 239
 240   AliMUONTrack *track, *nextTrack;
 241   AliMUONTrackParam *trackParam, *nextTrackParam;
 242   Int_t currentNRecTracks;
 243
 244   Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForTracking() *
 245                        GetRecoParam()->GetSigmaCutForTracking();
 246
 247   for (Int_t station = 2; station >= 0; station--) {
 248
 249     // Save the actual number of reconstructed track in case of
 250     // tracks are added or suppressed during the tracking procedure
 251     // !! Do not compress fRecTracksPtr until the end of the loop over tracks !!
 252     currentNRecTracks = fNRecTracks;
 253
 254     for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) {
 255       AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1));
 256
 257       track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
 258
 259       // Fit the track:
 260       // Do not take into account the multiple scattering to speed up the fit
 261       // Calculate the track parameter covariance matrix
 262       // If there is no cluster out of station 4 or 5 then use the vertex to better constrain the fit
 263       if (((AliMUONTrackParam*) track->GetTrackParamAtCluster()->First())->GetClusterPtr()->GetChamberId() > 5)
 264         Fit(*track, kFALSE, kTRUE, kTRUE);
 265       else Fit(*track, kFALSE, kFALSE, kTRUE);
 266
 267       // remove tracks out of limits
 268       if (!IsAcceptable(*((AliMUONTrackParam*)track->GetTrackParamAtCluster()->First()))) {
 269         fRecTracksPtr->Remove(track);
 270         fNRecTracks--;
 271         continue;
 272       }
 273
 274       // remove track if the normalized chi2 is too high
 275       if (track->GetNormalizedChi2() > sigmaCut2) {
 276         fRecTracksPtr->Remove(track);
 277         fNRecTracks--;
 278         continue;
 279       }
 280
 281       // save parameters from fit into smoothed parameters to complete track afterward
 282       if (GetRecoParam()->ComplementTracks()) {
 283
 284         if (station==2) { // save track parameters on stations 4 and 5
 285
 286           // extrapolate track parameters and covariances at each cluster
 287           // remove the track in case of failure
 288           if (!track->UpdateCovTrackParamAtCluster()) {
 289             fRecTracksPtr->Remove(track);
 290             fNRecTracks--;
 291             continue;
 292           }
 293
 294           // save them
 295           trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
 296           while (trackParam) {
 297             trackParam->SetSmoothParameters(trackParam->GetParameters());
 298             trackParam->SetSmoothCovariances(trackParam->GetCovariances());
 299             trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(trackParam);
 300           }
 301
 302         } else { // or save track parameters on last station only
 303
 304           trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
 305           if (trackParam->GetClusterPtr()->GetChamberId() < 2*(station+2)) {
 306
 307             // save parameters from fit
 308             trackParam->SetSmoothParameters(trackParam->GetParameters());
 309             trackParam->SetSmoothCovariances(trackParam->GetCovariances());
 310
 311             // save parameters extrapolated to the second chamber of the same station if it has been hit
 312             nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(trackParam);
 313             if (nextTrackParam->GetClusterPtr()->GetChamberId() < 2*(station+2)) {
 314
 315               // reset parameters and covariances
 316               nextTrackParam->SetParameters(trackParam->GetParameters());
 317               nextTrackParam->SetZ(trackParam->GetZ());
 318               nextTrackParam->SetCovariances(trackParam->GetCovariances());
 319
 320               // extrapolate them to the z of the corresponding cluster
 321               // remove the track in case of failure
 322               if (!AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetClusterPtr()->GetZ())) {
 323                 fRecTracksPtr->Remove(track);
 324                 fNRecTracks--;
 325                 continue;
 326               }
 327
 328               // save them
 329               nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters());
 330               nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances());
 331
 332             }
 333
 334           }
 335
 336         }
 337
 338       }
 339
 340       // Look for compatible cluster(s) in station(0..) "station"
 341       if (!FollowTrackInStation(*track, clusterStore, station)) {
 342
 343         // Try to recover track if required
 344         if (GetRecoParam()->RecoverTracks()) {
 345
 346           // work on a copy of the track if this station is not required
 347           // to keep the case where no cluster is reconstructed as a possible candidate
 348           if (!GetRecoParam()->RequestStation(station)) {
 349             track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*track);
 350             fNRecTracks++;
 351           }
 352
 353           // try to recover
 354           if (!RecoverTrack(*track, clusterStore, station)) {
 355             // remove track if no cluster found
 356             fRecTracksPtr->Remove(track);
 357             fNRecTracks--;
 358           }
 359
 360         } else if (GetRecoParam()->RequestStation(station)) {
 361           // remove track if no cluster found
 362           fRecTracksPtr->Remove(track);
 363           fNRecTracks--;
 364         }
 365
 366       }
 367
 368       // abort tracking if there are too many candidates
 369       if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
 370         AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks));
 371         return kFALSE;
 372       }
 373
 374     }
 375
 376     // Compress fRecTracksPtr for the next step
 377     fRecTracksPtr->Compress();
 378
 379     // Keep only the best tracks if required
 380     if (!GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
 381
 382   }
 383
 384   // Last fit of track candidates with all stations
 385   // Take into account the multiple scattering and remove bad tracks
 386   Int_t trackIndex = -1;
 387   track = (AliMUONTrack*) fRecTracksPtr->First();
 388   while (track) {
 389
 390     trackIndex++;
 391     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
 392
 393     Fit(*track, kTRUE, kFALSE, kTRUE);
 394
 395     // Printout for debuging
 396     if (AliLog::GetGlobalDebugLevel() >= 3) {
 397       cout << "FollowTracks: track candidate(0..) " << trackIndex << " after final fit" << endl;
 398       track->RecursiveDump();
 399     }
 400
 401     // Remove the track if the normalized chi2 is too high
 402     if (track->GetNormalizedChi2() > sigmaCut2) {
 403       fRecTracksPtr->Remove(track);
 404       fNRecTracks--;
 405       track = nextTrack;
 406       continue;
 407     }
 408
 409     // save parameters from fit into smoothed parameters to complete track afterward
 410     if (GetRecoParam()->ComplementTracks()) {
 411
 412       trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
 413       if (trackParam->GetClusterPtr()->GetChamberId() < 2) {
 414
 415         // save parameters from fit
 416         trackParam->SetSmoothParameters(trackParam->GetParameters());
 417         trackParam->SetSmoothCovariances(trackParam->GetCovariances());
 418
 419         // save parameters extrapolated to the second chamber of the same station if it has been hit
 420         nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(trackParam);
 421         if (nextTrackParam->GetClusterPtr()->GetChamberId() < 2) {
 422
 423           // reset parameters and covariances
 424           nextTrackParam->SetParameters(trackParam->GetParameters());
 425           nextTrackParam->SetZ(trackParam->GetZ());
 426           nextTrackParam->SetCovariances(trackParam->GetCovariances());
 427
 428           // extrapolate them to the z of the corresponding cluster
 429           // remove the track in case of failure
 430           if (!AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetClusterPtr()->GetZ())) {
 431             fRecTracksPtr->Remove(track);
 432             fNRecTracks--;
 433             track = nextTrack;
 434             continue;
 435           }
 436
 437           // save them
 438           nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters());
 439           nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances());
 440
 441         }
 442
 443       }
 444
 445     }
 446
 447     track = nextTrack;
 448
 449   }
 450
 451   fRecTracksPtr->Compress();
 452
 453   return kTRUE;
 454
 455 }
 456
 457   //__________________________________________________________________________
 458 Bool_t AliMUONTrackReconstructor::FollowTrackInChamber(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextChamber)
 459 {
 460   /// Follow trackCandidate in chamber(0..) nextChamber and search for compatible cluster(s)
 461   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
 462   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
 463   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
 464   ///         Remove the obsolete "trackCandidate" at the end.
 465   /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
 466   AliDebug(1,Form("Enter FollowTrackInChamber(1..) %d", nextChamber+1));
 467
 468   Double_t chi2WithOneCluster = AliMUONTrack::MaxChi2();
 469   Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
 470                                         GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
 471   Double_t bestChi2WithOneCluster = maxChi2WithOneCluster;
 472   Bool_t foundOneCluster = kFALSE;
 473   AliMUONTrack *newTrack = 0x0;
 474   AliMUONVCluster *cluster;
 475   AliMUONTrackParam extrapTrackParam;
 476   AliMUONTrackParam extrapTrackParamAtCh;
 477   AliMUONTrackParam extrapTrackParamAtCluster;
 478   AliMUONTrackParam bestTrackParamAtCluster;
 479
 480   // Get track parameters according to the propagation direction
 481   if (nextChamber > 7) extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last();
 482   else extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First();
 483
 484   // Printout for debuging
 485   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 486     cout<<endl<<"Track parameters and covariances at first cluster:"<<endl;
 487     extrapTrackParamAtCh.GetParameters().Print();
 488     extrapTrackParamAtCh.GetCovariances().Print();
 489   }
 490
 491   // Add MCS effect
 492   Int_t currentChamber = extrapTrackParamAtCh.GetClusterPtr()->GetChamberId();
 493   AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.);
 494
 495   // Add MCS in the missing chamber(s) if any
 496   while (currentChamber > nextChamber + 1) {
 497     // extrapolation to the missing chamber
 498     currentChamber--;
 499     if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(currentChamber))) return kFALSE;
 500     // add MCS effect
 501     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.);
 502   }
 503
 504   //Extrapolate trackCandidate to chamber
 505   if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(nextChamber))) return kFALSE;
 506
 507   // Printout for debuging
 508   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 509     cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(nextChamber)<<":"<<endl;
 510     extrapTrackParamAtCh.GetParameters().Print();
 511     extrapTrackParamAtCh.GetCovariances().Print();
 512   }
 513
 514   // Printout for debuging
 515   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 516     cout << "FollowTrackInStation: look for clusters in chamber(1..): " << nextChamber+1 << endl;
 517   }
 518
 519   // Ask the clustering to reconstruct new clusters around the track position in the current chamber
 520   // except for station 4 and 5 that are already entirely clusterized
 521   if (GetRecoParam()->CombineClusterTrackReco()) {
 522     if (nextChamber < 6) AskForNewClustersInChamber(extrapTrackParamAtCh, clusterStore, nextChamber);
 523   }
 524
 525   // Create iterators to loop over clusters in both chambers
 526   TIter next(clusterStore.CreateChamberIterator(nextChamber,nextChamber));
 527
 528   // look for cluster in chamber
 529   while ( ( cluster = static_cast<AliMUONVCluster*>(next()) ) ) {
 530
 531     // try to add the current cluster fast
 532     if (!TryOneClusterFast(extrapTrackParamAtCh, cluster)) continue;
 533
 534     // try to add the current cluster accuratly
 535     chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCh, cluster, extrapTrackParamAtCluster);
 536
 537     // if good chi2 then consider to add cluster
 538     if (chi2WithOneCluster < maxChi2WithOneCluster) {
 539       foundOneCluster = kTRUE;
 540
 541       // Printout for debuging
 542       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 543         cout << "FollowTrackInStation: found one cluster in chamber(1..): " << nextChamber+1
 544         << " (Chi2 = " << chi2WithOneCluster << ")" << endl;
 545         cluster->Print();
 546       }
 547
 548       if (GetRecoParam()->TrackAllTracks()) {
 549         // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
 550         newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 551         UpdateTrack(*newTrack,extrapTrackParamAtCluster);
 552         fNRecTracks++;
 553
 554         // Printout for debuging
 555         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 556           cout << "FollowTrackInStation: added one cluster in chamber(1..): " << nextChamber+1 << endl;
 557           if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 558         }
 559
 560       } else if (chi2WithOneCluster < bestChi2WithOneCluster) {
 561         // keep track of the best single cluster except if a couple of clusters has already been found
 562         bestChi2WithOneCluster = chi2WithOneCluster;
 563         bestTrackParamAtCluster = extrapTrackParamAtCluster;
 564       }
 565
 566     }
 567
 568   }
 569
 570   // fill out the best track if required else clean up the fRecTracksPtr array
 571   if (!GetRecoParam()->TrackAllTracks()) {
 572     if (foundOneCluster) {
 573       UpdateTrack(trackCandidate,bestTrackParamAtCluster);
 574
 575       // Printout for debuging
 576       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 577         cout << "FollowTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster.GetClusterPtr()->GetChamberId()+1 << endl;
 578         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 579       }
 580
 581     } else return kFALSE;
 582
 583   } else if (foundOneCluster) {
 584
 585     // remove obsolete track
 586     fRecTracksPtr->Remove(&trackCandidate);
 587     fNRecTracks--;
 588
 589   } else return kFALSE;
 590
 591   return kTRUE;
 592
 593 }
 594
 595   //__________________________________________________________________________
 596 Bool_t AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation)
 597 {
 598   /// Follow trackCandidate in station(0..) nextStation and search for compatible cluster(s)
 599   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
 600   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
 601   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
 602   ///         Remove the obsolete "trackCandidate" at the end.
 603   /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
 604   AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1));
 605
 606   // Order the chamber according to the propagation direction (tracking starts with chamber 2):
 607   // - nextStation == station(1...) 5 => forward propagation
 608   // - nextStation < station(1...) 5 => backward propagation
 609   Int_t ch1, ch2;
 610   if (nextStation==4) {
 611     ch1 = 2*nextStation+1;
 612     ch2 = 2*nextStation;
 613   } else {
 614     ch1 = 2*nextStation;
 615     ch2 = 2*nextStation+1;
 616   }
 617
 618   Double_t chi2WithOneCluster = AliMUONTrack::MaxChi2();
 619   Double_t chi2WithTwoClusters = AliMUONTrack::MaxChi2();
 620   Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
 621                                         GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
 622   Double_t maxChi2WithTwoClusters = 4. * GetRecoParam()->GetSigmaCutForTracking() *
 623                                          GetRecoParam()->GetSigmaCutForTracking(); // 4 because 4 quantities in chi2
 624   Double_t bestChi2WithOneCluster = maxChi2WithOneCluster;
 625   Double_t bestChi2WithTwoClusters = maxChi2WithTwoClusters;
 626   Bool_t foundOneCluster = kFALSE;
 627   Bool_t foundTwoClusters = kFALSE;
 628   AliMUONTrack *newTrack = 0x0;
 629   AliMUONVCluster *clusterCh1, *clusterCh2;
 630   AliMUONTrackParam extrapTrackParam;
 631   AliMUONTrackParam extrapTrackParamAtCh;
 632   AliMUONTrackParam extrapTrackParamAtCluster1;
 633   AliMUONTrackParam extrapTrackParamAtCluster2;
 634   AliMUONTrackParam bestTrackParamAtCluster1;
 635   AliMUONTrackParam bestTrackParamAtCluster2;
 636
 637   // Get track parameters according to the propagation direction
 638   if (nextStation==4) extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last();
 639   else extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First();
 640
 641   // Printout for debuging
 642   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 643     cout<<endl<<"Track parameters and covariances at first cluster:"<<endl;
 644     extrapTrackParamAtCh.GetParameters().Print();
 645     extrapTrackParamAtCh.GetCovariances().Print();
 646   }
 647
 648   // Add MCS effect
 649   Int_t currentChamber = extrapTrackParamAtCh.GetClusterPtr()->GetChamberId();
 650   AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.);
 651
 652   // Add MCS in the missing chamber(s) if any
 653   while (ch1 < ch2 && currentChamber > ch2 + 1) {
 654     // extrapolation to the missing chamber
 655     currentChamber--;
 656     if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(currentChamber))) return kFALSE;
 657     // add MCS effect
 658     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.);
 659   }
 660
 661   //Extrapolate trackCandidate to chamber "ch2"
 662   if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2))) return kFALSE;
 663
 664   // Printout for debuging
 665   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 666     cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch2)<<":"<<endl;
 667     extrapTrackParamAtCh.GetParameters().Print();
 668     extrapTrackParamAtCh.GetCovariances().Print();
 669   }
 670
 671   // Printout for debuging
 672   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 673     cout << "FollowTrackInStation: look for clusters in chamber(1..): " << ch2+1 << endl;
 674   }
 675
 676   // Ask the clustering to reconstruct new clusters around the track position in the current station
 677   // except for station 4 and 5 that are already entirely clusterized
 678   if (GetRecoParam()->CombineClusterTrackReco()) {
 679     if (nextStation < 3) AskForNewClustersInStation(extrapTrackParamAtCh, clusterStore, nextStation);
 680   }
 681
 682   Int_t nClusters = clusterStore.GetSize();
 683   Bool_t *clusterCh1Used = new Bool_t[nClusters];
 684   for (Int_t i = 0; i < nClusters; i++) clusterCh1Used[i] = kFALSE;
 685   Int_t iCluster1;
 686
 687   // Create iterators to loop over clusters in both chambers
 688   TIter nextInCh1(clusterStore.CreateChamberIterator(ch1,ch1));
 689   TIter nextInCh2(clusterStore.CreateChamberIterator(ch2,ch2));
 690
 691   // look for candidates in chamber 2
 692   while ( ( clusterCh2 = static_cast<AliMUONVCluster*>(nextInCh2()) ) ) {
 693
 694     // try to add the current cluster fast
 695     if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh2)) continue;
 696
 697     // try to add the current cluster accuratly
 698     chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh2, extrapTrackParamAtCluster2);
 699
 700     // if good chi2 then try to attach a cluster in the other chamber too
 701     if (chi2WithOneCluster < maxChi2WithOneCluster) {
 702       Bool_t foundSecondCluster = kFALSE;
 703
 704       // Printout for debuging
 705       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 706         cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch2+1
 707              << " (Chi2 = " << chi2WithOneCluster << ")" << endl;
 708         clusterCh2->Print();
 709         cout << "                      look for second clusters in chamber(1..): " << ch1+1 << " ..." << endl;
 710       }
 711
 712       // add MCS effect for next step
 713       AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCluster2,AliMUONConstants::ChamberThicknessInX0(ch2),-1.);
 714
 715       // copy new track parameters for next step
 716       extrapTrackParam = extrapTrackParamAtCluster2;
 717
 718       //Extrapolate track parameters to chamber "ch1"
 719       Bool_t normalExtrap = AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(ch1));
 720
 721       // reset cluster iterator of chamber 1
 722       nextInCh1.Reset();
 723       iCluster1 = -1;
 724
 725       // look for second candidates in chamber 1
 726       if (normalExtrap) while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
 727         iCluster1++;
 728
 729         // try to add the current cluster fast
 730         if (!TryOneClusterFast(extrapTrackParam, clusterCh1)) continue;
 731
 732         // try to add the current cluster accurately
 733         chi2WithTwoClusters = TryTwoClusters(extrapTrackParamAtCluster2, clusterCh1, extrapTrackParamAtCluster1);
 734
 735         // if good chi2 then create a new track by adding the 2 clusters to the "trackCandidate"
 736         if (chi2WithTwoClusters < maxChi2WithTwoClusters) {
 737           foundSecondCluster = kTRUE;
 738           foundTwoClusters = kTRUE;
 739
 740           // Printout for debuging
 741           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 742             cout << "FollowTrackInStation: found second cluster in chamber(1..): " << ch1+1
 743                  << " (Global Chi2 = " << chi2WithTwoClusters << ")" << endl;
 744             clusterCh1->Print();
 745           }
 746
 747           if (GetRecoParam()->TrackAllTracks()) {
 748             // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new clusters
 749             newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 750             UpdateTrack(*newTrack,extrapTrackParamAtCluster1,extrapTrackParamAtCluster2);
 751             fNRecTracks++;
 752
 753             // Tag clusterCh1 as used
 754             clusterCh1Used[iCluster1] = kTRUE;
 755
 756             // Printout for debuging
 757             if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 758               cout << "FollowTrackInStation: added two clusters in station(1..): " << nextStation+1 << endl;
 759               if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 760             }
 761
 762           } else if (chi2WithTwoClusters < bestChi2WithTwoClusters) {
 763             // keep track of the best couple of clusters
 764             bestChi2WithTwoClusters = chi2WithTwoClusters;
 765             bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
 766             bestTrackParamAtCluster2 = extrapTrackParamAtCluster2;
 767           }
 768
 769         }
 770
 771       }
 772
 773       // if no clusterCh1 found then consider to add clusterCh2 only
 774       if (!foundSecondCluster) {
 775         foundOneCluster = kTRUE;
 776
 777         if (GetRecoParam()->TrackAllTracks()) {
 778           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
 779           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 780           UpdateTrack(*newTrack,extrapTrackParamAtCluster2);
 781           fNRecTracks++;
 782
 783           // Printout for debuging
 784           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 785             cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch2+1 << endl;
 786             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 787           }
 788
 789         } else if (!foundTwoClusters && chi2WithOneCluster < bestChi2WithOneCluster) {
 790           // keep track of the best single cluster except if a couple of clusters has already been found
 791           bestChi2WithOneCluster = chi2WithOneCluster;
 792           bestTrackParamAtCluster1 = extrapTrackParamAtCluster2;
 793         }
 794
 795       }
 796
 797     }
 798
 799   }
 800
 801   // look for candidates in chamber 1 not already attached to a track
 802   // if we want to keep all possible tracks or if no good couple of clusters has been found
 803   if (GetRecoParam()->TrackAllTracks() || !foundTwoClusters) {
 804
 805     // add MCS effect for next step
 806     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(ch2),-1.);
 807
 808     //Extrapolate trackCandidate to chamber "ch1"
 809     Bool_t normalExtrap = AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch1));
 810
 811     // Printout for debuging
 812     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 813       cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch1)<<":"<<endl;
 814       extrapTrackParamAtCh.GetParameters().Print();
 815       extrapTrackParamAtCh.GetCovariances().Print();
 816     }
 817
 818     // Printout for debuging
 819     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 820       cout << "FollowTrackInStation: look for single clusters in chamber(1..): " << ch1+1 << endl;
 821     }
 822
 823     // reset cluster iterator of chamber 1
 824     nextInCh1.Reset();
 825     iCluster1 = -1;
 826
 827     // look for second candidates in chamber 1
 828     if (normalExtrap) while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
 829       iCluster1++;
 830
 831       if (clusterCh1Used[iCluster1]) continue; // Skip cluster already used
 832
 833       // try to add the current cluster fast
 834       if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh1)) continue;
 835
 836       // try to add the current cluster accuratly
 837       chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh1, extrapTrackParamAtCluster1);
 838
 839       // if good chi2 then consider to add clusterCh1
 840       // We do not try to attach a cluster in the other chamber too since it has already been done above
 841       if (chi2WithOneCluster < maxChi2WithOneCluster) {
 842         foundOneCluster = kTRUE;
 843
 844         // Printout for debuging
 845         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 846           cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch1+1
 847                << " (Chi2 = " << chi2WithOneCluster << ")" << endl;
 848           clusterCh1->Print();
 849         }
 850
 851         if (GetRecoParam()->TrackAllTracks()) {
 852           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
 853           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 854           UpdateTrack(*newTrack,extrapTrackParamAtCluster1);
 855           fNRecTracks++;
 856
 857           // Printout for debuging
 858           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 859             cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch1+1 << endl;
 860             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 861           }
 862
 863         } else if (chi2WithOneCluster < bestChi2WithOneCluster) {
 864           // keep track of the best single cluster except if a couple of clusters has already been found
 865           bestChi2WithOneCluster = chi2WithOneCluster;
 866           bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
 867         }
 868
 869       }
 870
 871     }
 872
 873   }
 874
 875   // fill out the best track if required else clean up the fRecTracksPtr array
 876   if (!GetRecoParam()->TrackAllTracks()) {
 877     if (foundTwoClusters) {
 878       UpdateTrack(trackCandidate,bestTrackParamAtCluster1,bestTrackParamAtCluster2);
 879
 880       // Printout for debuging
 881       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 882         cout << "FollowTrackInStation: added the two best clusters in station(1..): " << nextStation+1 << endl;
 883         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 884       }
 885
 886     } else if (foundOneCluster) {
 887       UpdateTrack(trackCandidate,bestTrackParamAtCluster1);
 888
 889       // Printout for debuging
 890       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 891         cout << "FollowTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster1.GetClusterPtr()->GetChamberId()+1 << endl;
 892         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 893       }
 894
 895     } else {
 896       delete [] clusterCh1Used;
 897       return kFALSE;
 898     }
 899
 900   } else if (foundOneCluster || foundTwoClusters) {
 901
 902     // remove obsolete track
 903     fRecTracksPtr->Remove(&trackCandidate);
 904     fNRecTracks--;
 905
 906   } else {
 907     delete [] clusterCh1Used;
 908     return kFALSE;
 909   }
 910
 911   delete [] clusterCh1Used;
 912   return kTRUE;
 913
 914 }
 915
 916   //__________________________________________________________________________
 917 Double_t AliMUONTrackReconstructor::TryTwoClusters(const AliMUONTrackParam &trackParamAtCluster1, AliMUONVCluster* cluster2,
 918                                                    AliMUONTrackParam &trackParamAtCluster2)
 919 {
 920 /// Test the compatibility between the track and the 2 clusters together (using trackParam's covariance matrix):
 921 /// return the corresponding Chi2 accounting for covariances between the 2 clusters
 922 /// return trackParamAtCluster1 & 2
 923
 924   // extrapolate track parameters at the z position of the second cluster (no need to extrapolate the covariances)
 925   // and set pointer to cluster into trackParamAtCluster
 926   trackParamAtCluster2.SetParameters(trackParamAtCluster1.GetParameters());
 927   trackParamAtCluster2.SetZ(trackParamAtCluster1.GetZ());
 928   trackParamAtCluster2.SetClusterPtr(cluster2);
 929   if (!AliMUONTrackExtrap::ExtrapToZ(&trackParamAtCluster2, cluster2->GetZ())) return 2.*AliMUONTrack::MaxChi2();
 930
 931   // Set differences between track and the 2 clusters in the bending and non bending directions
 932   AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr();
 933   TMatrixD dPos(4,1);
 934   dPos(0,0) = cluster1->GetX() - trackParamAtCluster1.GetNonBendingCoor();
 935   dPos(1,0) = cluster1->GetY() - trackParamAtCluster1.GetBendingCoor();
 936   dPos(2,0) = cluster2->GetX() - trackParamAtCluster2.GetNonBendingCoor();
 937   dPos(3,0) = cluster2->GetY() - trackParamAtCluster2.GetBendingCoor();
 938
 939   // Calculate the error matrix from the track parameter covariances at first cluster
 940   TMatrixD error(4,4);
 941   error.Zero();
 942   if (trackParamAtCluster1.CovariancesExist()) {
 943     // Save track parameters at first cluster
 944     TMatrixD paramAtCluster1Save(trackParamAtCluster1.GetParameters());
 945
 946     // Save track coordinates at second cluster
 947     Double_t nonBendingCoor2 = trackParamAtCluster2.GetNonBendingCoor();
 948     Double_t bendingCoor2    = trackParamAtCluster2.GetBendingCoor();
 949
 950     // copy track parameters at first cluster for jacobian calculation
 951     AliMUONTrackParam trackParam(trackParamAtCluster1);
 952
 953     // Get the pointer to the parameter covariance matrix at first cluster
 954     const TMatrixD& kParamCov = trackParam.GetCovariances();
 955
 956     // Calculate the jacobian related to the transformation between track parameters
 957     // at first cluster and track coordinates at the 2 cluster z-positions
 958     TMatrixD jacob(4,5);
 959     jacob.Zero();
 960     // first derivative at the first cluster:
 961     jacob(0,0) = 1.; // dx1/dx
 962     jacob(1,2) = 1.; // dy1/dy
 963     // first derivative at the second cluster:
 964     TMatrixD dParam(5,1);
 965     Double_t direction[5] = {-1.,-1.,1.,1.,-1.};
 966     for (Int_t i=0; i<5; i++) {
 967       // Skip jacobian calculation for parameters with no associated error
 968       if (kParamCov(i,i) == 0.) continue;
 969       // Small variation of parameter i only
 970       for (Int_t j=0; j<5; j++) {
 971         if (j==i) {
 972           dParam(j,0) = TMath::Sqrt(kParamCov(i,i));
 973           dParam(j,0) *= TMath::Sign(1.,direction[j]*paramAtCluster1Save(j,0)); // variation always in the same direction
 974         } else dParam(j,0) = 0.;
 975       }
 976
 977       // Set new track parameters at first cluster
 978       trackParam.SetParameters(paramAtCluster1Save);
 979       trackParam.AddParameters(dParam);
 980       trackParam.SetZ(cluster1->GetZ());
 981
 982       // Extrapolate new track parameters to the z position of the second cluster
 983       if (!AliMUONTrackExtrap::ExtrapToZ(&trackParam,cluster2->GetZ())) return 2.*AliMUONTrack::MaxChi2();
 984
 985       // Calculate the jacobian
 986       jacob(2,i) = (trackParam.GetNonBendingCoor() - nonBendingCoor2) / dParam(i,0); // dx2/dParami
 987       jacob(3,i) = (trackParam.GetBendingCoor()    - bendingCoor2   ) / dParam(i,0); // dy2/dParami
 988     }
 989
 990     // Calculate the error matrix
 991     TMatrixD tmp(jacob,TMatrixD::kMult,kParamCov);
 992     error = TMatrixD(tmp,TMatrixD::kMultTranspose,jacob);
 993   }
 994
 995   // Add cluster resolution to the error matrix
 996   error(0,0) += cluster1->GetErrX2();
 997   error(1,1) += cluster1->GetErrY2();
 998   error(2,2) += cluster2->GetErrX2();
 999   error(3,3) += cluster2->GetErrY2();
1000
1001   // invert the error matrix for Chi2 calculation
1002   if (error.Determinant() != 0) {
1003     error.Invert();
1004   } else {
1005     AliWarning(" Determinant error=0");
1006     return 2.*AliMUONTrack::MaxChi2();
1007   }
1008
1009   // Compute the Chi2 value
1010   TMatrixD tmp2(dPos,TMatrixD::kTransposeMult,error);
1011   TMatrixD result(tmp2,TMatrixD::kMult,dPos);
1012
1013   return result(0,0);
1014
1015 }
1016
1017   //__________________________________________________________________________
1018 void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster)
1019 {
1020   /// Add 1 cluster to the track candidate
1021   /// Update chi2 of the track
1022
1023   // Compute local chi2
1024   AliMUONVCluster* cluster = trackParamAtCluster.GetClusterPtr();
1025   Double_t deltaX = trackParamAtCluster.GetNonBendingCoor() - cluster->GetX();
1026   Double_t deltaY = trackParamAtCluster.GetBendingCoor() - cluster->GetY();
1027   Double_t localChi2 = deltaX*deltaX / cluster->GetErrX2() +
1028                        deltaY*deltaY / cluster->GetErrY2();
1029
1030   // Flag cluster as being not removable
1031   if (GetRecoParam()->RequestStation(cluster->GetChamberId()/2))
1032     trackParamAtCluster.SetRemovable(kFALSE);
1033   else trackParamAtCluster.SetRemovable(kTRUE);
1034   trackParamAtCluster.SetLocalChi2(0.); // --> Local chi2 not used
1035
1036   // Update the chi2 of the new track
1037   track.SetGlobalChi2(track.GetGlobalChi2() + localChi2);
1038
1039   // Update TrackParamAtCluster
1040   track.AddTrackParamAtCluster(trackParamAtCluster,*cluster);
1041
1042 }
1043
1044   //__________________________________________________________________________
1045 void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster1, AliMUONTrackParam &trackParamAtCluster2)
1046 {
1047   /// Add 2 clusters to the track candidate
1048   /// Update track and local chi2
1049
1050   // Update local chi2 at first cluster
1051   AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr();
1052   Double_t deltaX = trackParamAtCluster1.GetNonBendingCoor() - cluster1->GetX();
1053   Double_t deltaY = trackParamAtCluster1.GetBendingCoor() - cluster1->GetY();
1054   Double_t localChi2AtCluster1 = deltaX*deltaX / cluster1->GetErrX2() +
1055                                  deltaY*deltaY / cluster1->GetErrY2();
1056   trackParamAtCluster1.SetLocalChi2(localChi2AtCluster1);
1057
1058   // Flag first cluster as being removable
1059   trackParamAtCluster1.SetRemovable(kTRUE);
1060
1061   // Update local chi2 at second cluster
1062   AliMUONVCluster* cluster2 = trackParamAtCluster2.GetClusterPtr();
1063   deltaX = trackParamAtCluster2.GetNonBendingCoor() - cluster2->GetX();
1064   deltaY = trackParamAtCluster2.GetBendingCoor() - cluster2->GetY();
1065   Double_t localChi2AtCluster2 = deltaX*deltaX / cluster2->GetErrX2() +
1066                                  deltaY*deltaY / cluster2->GetErrY2();
1067   trackParamAtCluster2.SetLocalChi2(localChi2AtCluster2);
1068
1069   // Flag first cluster as being removable
1070   trackParamAtCluster2.SetRemovable(kTRUE);
1071
1072   // Update the chi2 of the new track
1073   track.SetGlobalChi2(track.GetGlobalChi2() + localChi2AtCluster1 + localChi2AtCluster2);
1074
1075   // Update TrackParamAtCluster
1076   track.AddTrackParamAtCluster(trackParamAtCluster1,*cluster1);
1077   track.AddTrackParamAtCluster(trackParamAtCluster2,*cluster2);
1078
1079 }
1080
1081   //__________________________________________________________________________
1082 Bool_t AliMUONTrackReconstructor::RecoverTrack(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation)
1083 {
1084   /// Try to recover the track candidate in the next station
1085   /// by removing the worst of the two clusters attached in the current station
1086   /// Return kTRUE if recovering succeeds
1087   AliDebug(1,"Enter RecoverTrack");
1088
1089   // Do not try to recover track until we have attached cluster(s) on station(1..) 3
1090   if (nextStation > 1) return kFALSE;
1091
1092   Int_t worstClusterNumber = -1;
1093   Double_t localChi2, worstLocalChi2 = -1.;
1094
1095   // Look for the cluster to remove
1096   for (Int_t clusterNumber = 0; clusterNumber < 2; clusterNumber++) {
1097     AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(clusterNumber);
1098
1099     // check if current cluster is in the previous station
1100     if (trackParamAtCluster->GetClusterPtr()->GetChamberId()/2 != nextStation+1) break;
1101
1102     // check if current cluster is removable
1103     if (!trackParamAtCluster->IsRemovable()) return kFALSE;
1104
1105     // reset the current cluster as beig not removable if it is on a required station
1106     if (GetRecoParam()->RequestStation(nextStation+1)) trackParamAtCluster->SetRemovable(kFALSE);
1107
1108     // Pick up cluster with the worst chi2
1109     localChi2 = trackParamAtCluster->GetLocalChi2();
1110     if (localChi2 > worstLocalChi2) {
1111       worstLocalChi2 = localChi2;
1112       worstClusterNumber = clusterNumber;
1113     }
1114
1115   }
1116
1117   // check if worst cluster found
1118   if (worstClusterNumber < 0) return kFALSE;
1119
1120   // Remove the worst cluster
1121   trackCandidate.RemoveTrackParamAtCluster((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(worstClusterNumber));
1122
1123   // Re-fit the track:
1124   // Do not take into account the multiple scattering to speed up the fit
1125   // Calculate the track parameter covariance matrix
1126   Fit(trackCandidate, kFALSE, kFALSE, kTRUE);
1127
1128   // skip track if the normalized chi2 is too high
1129   if (trackCandidate.GetNormalizedChi2() > GetRecoParam()->GetSigmaCutForTracking() * GetRecoParam()->GetSigmaCutForTracking()) return kFALSE;
1130
1131   // skip track out of limits
1132   if (!IsAcceptable(*((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First()))) return kFALSE;
1133
1134   // Look for new cluster(s) in next station
1135   return FollowTrackInStation(trackCandidate,clusterStore,nextStation);
1136
1137 }
1138
1139   //__________________________________________________________________________
1140 void AliMUONTrackReconstructor::SetVertexErrXY2ForFit(AliMUONTrack &trackCandidate)
1141 {
1142   /// Compute the vertex resolution square from natural vertex dispersion and
1143   /// multiple scattering effets according to trackCandidate path in absorber
1144   /// It is necessary to account for multiple scattering effects here instead of during the fit of
1145   /// the "trackCandidate" to do not influence the result by changing track resolution at vertex
1146   AliDebug(1,"Enter SetVertexForFit");
1147
1148   Double_t nonBendingReso2 = GetRecoParam()->GetNonBendingVertexDispersion() *
1149                              GetRecoParam()->GetNonBendingVertexDispersion();
1150   Double_t bendingReso2 = GetRecoParam()->GetBendingVertexDispersion() *
1151                           GetRecoParam()->GetBendingVertexDispersion();
1152
1153   // add multiple scattering effets
1154   AliMUONTrackParam paramAtVertex(*((AliMUONTrackParam*)(trackCandidate.GetTrackParamAtCluster()->First())));
1155   paramAtVertex.DeleteCovariances(); // to be sure to account only for multiple scattering
1156   if (!AliMUONTrackExtrap::ExtrapToZ(&paramAtVertex,AliMUONConstants::AbsZEnd())) {
1157     nonBendingReso2 = 0.;
1158     bendingReso2 = 0.;
1159   } else {
1160     AliMUONTrackExtrap::ExtrapToVertexUncorrected(&paramAtVertex,0.);
1161     const TMatrixD& kParamCov = paramAtVertex.GetCovariances();
1162     nonBendingReso2 += kParamCov(0,0);
1163     bendingReso2 += kParamCov(2,2);
1164   }
1165
1166   // Set the vertex resolution square
1167   trackCandidate.SetVertexErrXY2(nonBendingReso2,bendingReso2);
1168 }
1169
1170   //__________________________________________________________________________
1171 void AliMUONTrackReconstructor::Fit(AliMUONTrack &track, Bool_t includeMCS, Bool_t fitWithVertex, Bool_t calcCov)
1172 {
1173   /// Fit the track
1174   /// w/wo multiple Coulomb scattering according to "includeMCS".
1175   /// w/wo constraining the vertex according to "fitWithVertex".
1176   /// calculating or not the covariance matrix according to "calcCov".
1177   AliDebug(1,"Enter Fit");
1178
1179   Double_t benC, errorParam, invBenP, nonBenC, x, y;
1180   AliMUONTrackParam *trackParam;
1181   Double_t arg[1], fedm, errdef, globalChi2;
1182   Int_t npari, nparx;
1183   Int_t status, covStatus;
1184
1185   // Instantiate gMinuit if not already done
1186   if (!gMinuit) gMinuit = new TMinuit(6);
1187   // Clear MINUIT parameters
1188   gMinuit->mncler();
1189   // Give the fitted track to MINUIT
1190   gMinuit->SetObjectFit(&track);
1191   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
1192     // Define print level
1193     arg[0] = 1;
1194     gMinuit->mnexcm("SET PRI", arg, 1, status);
1195     // Print covariance matrix
1196     gMinuit->mnexcm("SHO COV", arg, 0, status);
1197   } else {
1198     arg[0] = -1;
1199     gMinuit->mnexcm("SET PRI", arg, 1, status);
1200   }
1201   // No warnings
1202   gMinuit->mnexcm("SET NOW", arg, 0, status);
1203   // Define strategy
1204   //arg[0] = 2;
1205   //gMinuit->mnexcm("SET STR", arg, 1, status);
1206
1207   // set flag w/wo multiple scattering according to "includeMCS"
1208   track.FitWithMCS(includeMCS);
1209   if (includeMCS) {
1210     // compute cluster weights only once
1211     if (!track.UpdateTrackParamAtCluster() || !track.ComputeClusterWeights()) {
1212       AliWarning("cannot take into account the multiple scattering effects");
1213       track.FitWithMCS(kFALSE);
1214     }
1215   }
1216
1217   track.FitWithVertex(fitWithVertex);
1218   if (fitWithVertex) SetVertexErrXY2ForFit(track);
1219
1220   // Set fitting function
1221   gMinuit->SetFCN(TrackChi2);
1222
1223   // Set fitted parameters (!! The order is very important for the covariance matrix !!)
1224   // Mandatory limits to avoid NaN values of parameters
1225   trackParam = (AliMUONTrackParam*) (track.GetTrackParamAtCluster()->First());
1226   Double_t maxIBM = 1. / GetRecoParam()->GetMinBendingMomentum();
1227   gMinuit->mnparm(0, "X", trackParam->GetNonBendingCoor(), 0.03, -500.0, 500.0, status);
1228   gMinuit->mnparm(1, "NonBenS", trackParam->GetNonBendingSlope(), 0.001, -1., 1., status);
1229   gMinuit->mnparm(2, "Y", trackParam->GetBendingCoor(), 0.10, -500.0, 500.0, status);
1230   gMinuit->mnparm(3, "BenS", trackParam->GetBendingSlope(), 0.001, -1.5, 1.5, status);
1231   gMinuit->mnparm(4, "InvBenP", trackParam->GetInverseBendingMomentum(), 0.003, -maxIBM, maxIBM, status);
1232
1233   // minimization
1234   gMinuit->mnexcm("MIGRAD", arg, 0, status);
1235
1236   // Calculate the covariance matrix more accurately if required
1237   if (calcCov) gMinuit->mnexcm("HESSE", arg, 0, status);
1238
1239   // get results into "invBenP", "benC", "nonBenC" ("x", "y")
1240   gMinuit->GetParameter(0, x, errorParam);
1241   trackParam->SetNonBendingCoor(x);
1242   gMinuit->GetParameter(1, nonBenC, errorParam);
1243   trackParam->SetNonBendingSlope(nonBenC);
1244   gMinuit->GetParameter(2, y, errorParam);
1245   trackParam->SetBendingCoor(y);
1246   gMinuit->GetParameter(3, benC, errorParam);
1247   trackParam->SetBendingSlope(benC);
1248   gMinuit->GetParameter(4, invBenP, errorParam);
1249   trackParam->SetInverseBendingMomentum(invBenP);
1250
1251   // global result of the fit
1252   gMinuit->mnstat(globalChi2, fedm, errdef, npari, nparx, covStatus);
1253   track.SetGlobalChi2(globalChi2);
1254
1255   // Get the covariance matrix if required
1256   if (calcCov) {
1257     // Covariance matrix according to HESSE status
1258     // If problem then keep only the diagonal terms (variances)
1259     Double_t matrix[5][5];
1260     for (Int_t i=0; i<5; i++) for (Int_t j=0; j<5; j++) matrix[i][j] = 0.;
1261     gMinuit->mnemat(&matrix[0][0],5);
1262     if (covStatus == 3) trackParam->SetCovariances(matrix);
1263     else trackParam->SetVariances(matrix);
1264   } else trackParam->DeleteCovariances();
1265
1266 }
1267
1268   //__________________________________________________________________________
1269 void TrackChi2(Int_t & /*nParam*/, Double_t * /*gradient*/, Double_t &chi2, Double_t *param, Int_t /*flag*/)
1270 {
1271   /// Return the "Chi2" to be minimized with Minuit for track fitting.
1272   /// Assumes that the trackParamAtCluster are sorted according to increasing Z.
1273   /// Track parameters at each cluster are updated accordingly.
1274   /// Vertex is used according to the flag "trackBeingFitted->GetFitWithVertex()".
1275   /// Multiple Coulomb scattering is taken into account according to the flag "trackBeingFitted->GetFitWithMCS()".
1276
1277   AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit();
1278   AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) trackBeingFitted->GetTrackParamAtCluster()->First();
1279   Double_t dX, dY;
1280   chi2 = 0.; // initialize chi2
1281
1282   // update track parameters
1283   trackParamAtCluster->SetNonBendingCoor(param[0]);
1284   trackParamAtCluster->SetNonBendingSlope(param[1]);
1285   trackParamAtCluster->SetBendingCoor(param[2]);
1286   trackParamAtCluster->SetBendingSlope(param[3]);
1287   trackParamAtCluster->SetInverseBendingMomentum(param[4]);
1288   if (!trackBeingFitted->UpdateTrackParamAtCluster()) {
1289     chi2 = 2.*AliMUONTrack::MaxChi2();
1290     return;
1291   }
1292
1293   // Take the vertex into account in the fit if required
1294   if (trackBeingFitted->FitWithVertex()) {
1295     Double_t nonBendingReso2,bendingReso2;
1296     trackBeingFitted->GetVertexErrXY2(nonBendingReso2,bendingReso2);
1297     AliMUONTrackParam paramAtVertex(*trackParamAtCluster);
1298     if (nonBendingReso2 != 0. && bendingReso2 != 0. && AliMUONTrackExtrap::ExtrapToZ(&paramAtVertex, 0.)) { // vextex position = (0,0,0)
1299       dX = paramAtVertex.GetNonBendingCoor();
1300       dY = paramAtVertex.GetBendingCoor();
1301       chi2 += dX * dX / nonBendingReso2 + dY * dY / bendingReso2;
1302     } else {
1303       chi2 = 2.*AliMUONTrack::MaxChi2();
1304       return;
1305     }
1306   }
1307
1308   // compute chi2 w/wo multiple scattering
1309   chi2 += trackBeingFitted->ComputeGlobalChi2(trackBeingFitted->FitWithMCS());
1310
1311 }
1312
1313   //__________________________________________________________________________
1314 Bool_t AliMUONTrackReconstructor::ComplementTracks(const AliMUONVClusterStore& clusterStore)
1315 {
1316   /// Complete tracks by adding missing clusters (if there is an overlap between
1317   /// two detection elements, the track may have two clusters in the same chamber).
1318   /// Re-fit track parameters and covariances.
1319   /// Return kTRUE if one or more tracks have been complemented.
1320   AliDebug(1,"Enter ComplementTracks");
1321
1322   Int_t chamberId, detElemId;
1323   Double_t chi2OfCluster, bestChi2OfCluster;
1324   Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForTracking() *
1325                        GetRecoParam()->GetSigmaCutForTracking();
1326   Bool_t foundOneCluster, trackModified, hasChanged = kFALSE;
1327   AliMUONVCluster* cluster;
1328   AliMUONTrackParam *trackParam, *nextTrackParam, copyOfTrackParam, trackParamAtCluster, bestTrackParamAtCluster;
1329
1330   AliMUONTrack *track = (AliMUONTrack*) fRecTracksPtr->First();
1331   while (track) {
1332     trackModified = kFALSE;
1333
1334     trackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->First();
1335     while (trackParam) {
1336       foundOneCluster = kFALSE;
1337       bestChi2OfCluster = 2. * sigmaCut2; // 2 because 2 quantities in chi2
1338       chamberId = trackParam->GetClusterPtr()->GetChamberId();
1339       detElemId = trackParam->GetClusterPtr()->GetDetElemId();
1340
1341       // prepare nextTrackParam before adding new cluster because of the sorting
1342       nextTrackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->After(trackParam);
1343
1344       // recover track parameters from local fit and put them into a copy of trackParam
1345       copyOfTrackParam.SetZ(trackParam->GetZ());
1346       copyOfTrackParam.SetParameters(trackParam->GetSmoothParameters());
1347       copyOfTrackParam.SetCovariances(trackParam->GetSmoothCovariances());
1348
1349       // Create iterators to loop over clusters in current chamber
1350       TIter nextInCh(clusterStore.CreateChamberIterator(chamberId,chamberId));
1351
1352       // look for one second candidate in the same chamber
1353       while ( ( cluster = static_cast<AliMUONVCluster*>(nextInCh()) ) ) {
1354
1355         // look for a cluster in another detection element
1356         if (cluster->GetDetElemId() == detElemId) continue;
1357
1358         // try to add the current cluster fast
1359         if (!TryOneClusterFast(copyOfTrackParam, cluster)) continue;
1360
1361         // try to add the current cluster accurately
1362         chi2OfCluster = TryOneCluster(copyOfTrackParam, cluster, trackParamAtCluster);
1363
1364         // if better chi2 then prepare to add this cluster to the track
1365         if (chi2OfCluster < bestChi2OfCluster) {
1366           bestChi2OfCluster = chi2OfCluster;
1367           bestTrackParamAtCluster = trackParamAtCluster;
1368           foundOneCluster = kTRUE;
1369         }
1370
1371       }
1372
1373       // add new cluster if any
1374       if (foundOneCluster) {
1375
1376         // Printout for debuging
1377         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1378           cout << "ComplementTracks: found one cluster in chamber(1..): " << chamberId+1 << endl;
1379           bestTrackParamAtCluster.GetClusterPtr()->Print();
1380           cout<<endl<<"Track parameters and covariances at cluster:"<<endl;
1381           bestTrackParamAtCluster.GetParameters().Print();
1382           bestTrackParamAtCluster.GetCovariances().Print();
1383         }
1384
1385         trackParam->SetRemovable(kTRUE);
1386         bestTrackParamAtCluster.SetRemovable(kTRUE);
1387         track->AddTrackParamAtCluster(bestTrackParamAtCluster,*(bestTrackParamAtCluster.GetClusterPtr()));
1388         trackModified = kTRUE;
1389         hasChanged = kTRUE;
1390       }
1391
1392       trackParam = nextTrackParam;
1393     }
1394
1395     // re-fit track parameters if needed
1396     if (trackModified) Fit(*track, kTRUE, kFALSE, kTRUE);
1397
1398     track = (AliMUONTrack*) fRecTracksPtr->After(track);
1399   }
1400
1401   return hasChanged;
1402
1403 }
1404
1405   //__________________________________________________________________________
1406 void AliMUONTrackReconstructor::ImproveTrack(AliMUONTrack &track)
1407 {
1408   /// Improve the given track by removing clusters with local chi2 highter than the defined cut
1409   /// Recompute track parameters and covariances at the remaining clusters
1410   AliDebug(1,"Enter ImproveTrack");
1411
1412   Double_t localChi2, worstLocalChi2;
1413   AliMUONTrackParam *trackParamAtCluster, *worstTrackParamAtCluster;
1414   Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForImprovement() *
1415                        GetRecoParam()->GetSigmaCutForImprovement();
1416
1417   while (!track.IsImproved()) {
1418
1419     // identify removable clusters
1420     track.TagRemovableClusters(GetRecoParam()->RequestedStationMask());
1421
1422     // Update track parameters and covariances
1423     if (!track.UpdateCovTrackParamAtCluster()) {
1424       AliWarning("unable to update track parameters and covariances --> stop improvement");
1425       break;
1426     }
1427
1428     // Compute local chi2 of each clusters
1429     track.ComputeLocalChi2(kTRUE);
1430
1431     // Look for the cluster to remove
1432     worstTrackParamAtCluster = NULL;
1433     worstLocalChi2 = 0.;
1434     trackParamAtCluster = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->First();
1435     while (trackParamAtCluster) {
1436
1437       // Pick up cluster with the worst chi2
1438       localChi2 = trackParamAtCluster->GetLocalChi2();
1439       if (localChi2 > worstLocalChi2) {
1440         worstLocalChi2 = localChi2;
1441         worstTrackParamAtCluster = trackParamAtCluster;
1442       }
1443
1444     trackParamAtCluster = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(trackParamAtCluster);
1445     }
1446
1447     // Check if worst cluster found
1448     if (!worstTrackParamAtCluster) {
1449       AliWarning("Bad local chi2 values?");
1450       break;
1451     }
1452
1453     // Check whether the worst chi2 is under requirement or not
1454     if (worstLocalChi2 < 2. * sigmaCut2) { // 2 because 2 quantities in chi2
1455       track.SetImproved(kTRUE);
1456       break;
1457     }
1458
1459     // if the worst cluster is not removable then stop improvement
1460     if (!worstTrackParamAtCluster->IsRemovable()) break;
1461
1462     // Remove the worst cluster
1463     track.RemoveTrackParamAtCluster(worstTrackParamAtCluster);
1464
1465     // Re-fit the track:
1466     // Take into account the multiple scattering
1467     // Calculate the track parameter covariance matrix
1468     Fit(track, kTRUE, kFALSE, kTRUE);
1469
1470     // Printout for debuging
1471     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1472       cout << "ImproveTracks: track " << fRecTracksPtr->IndexOf(&track)+1 << " improved " << endl;
1473     }
1474
1475   }
1476
1477 }
1478
1479   //__________________________________________________________________________
1480 Bool_t AliMUONTrackReconstructor::FinalizeTrack(AliMUONTrack &track)
1481 {
1482   /// Recompute track parameters and covariances at each attached cluster
1483   /// from those at the first one, if not already done
1484   AliDebug(1,"Enter FinalizeTrack");
1485   if (!track.IsImproved() && !track.UpdateCovTrackParamAtCluster()) {
1486     AliWarning("finalization failed due to extrapolation problem");
1487     return kFALSE;
1488   }
1489   return kTRUE;
1490 }
1491
1492   //__________________________________________________________________________
1493 Bool_t AliMUONTrackReconstructor::RefitTrack(AliMUONTrack &track, Bool_t enableImprovement)
1494 {
1495   /// re-fit the given track
1496
1497   // check validity of the track
1498   if (track.GetNClusters() < 3) {
1499     AliWarning("the track does not contain enough clusters --> unable to refit");
1500     return kFALSE;
1501   }
1502
1503   // reset the seed (i.e. parameters at first cluster) before fitting
1504   AliMUONTrackParam* firstTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->First();
1505   if (firstTrackParam->GetInverseBendingMomentum() == 0.) {
1506     AliWarning("track parameters at first chamber are not initialized --> unable to refit");
1507     return kFALSE;
1508   }
1509
1510   // compute track parameters at each cluster from parameters at the first one
1511   // necessary to compute multiple scattering effect during refitting
1512   if (!track.UpdateTrackParamAtCluster()) {
1513     AliWarning("bad track refitting due to extrapolation failure");
1514     return kFALSE;
1515   }
1516
1517   // Re-fit the track:
1518   // Take into account the multiple scattering
1519   // Calculate the track parameter covariance matrix
1520   Fit(track, kTRUE, kFALSE, kTRUE);
1521
1522   // Improve the reconstructed tracks if required
1523   track.SetImproved(kFALSE);
1524   if (enableImprovement && GetRecoParam()->ImproveTracks()) ImproveTrack(track);
1525
1526   // Fill AliMUONTrack data members
1527   if (track.GetGlobalChi2() < AliMUONTrack::MaxChi2()) return FinalizeTrack(track);
1528   else {
1529     AliWarning("track not finalized due to extrapolation failure");
1530     return kFALSE;
1531   }
1532
1533 }
1534