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