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