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   //__________________________________________________________________________
  52 AliMUONTrackReconstructor::AliMUONTrackReconstructor()
  53   : AliMUONVTrackReconstructor()
  54 {
  55   /// Constructor for class AliMUONTrackReconstructor
  56   AliInfo("*** Original tracking ***");
  57 }
  58
  59   //__________________________________________________________________________
  60 AliMUONTrackReconstructor::~AliMUONTrackReconstructor()
  61 {
  62 /// Destructor
  63 }
  64
  65   //__________________________________________________________________________
  66 void AliMUONTrackReconstructor::MakeTrackCandidates()
  67 {
  68   /// To make track candidates (assuming linear propagation if the flag fgkMakeTrackCandidatesFast is set to kTRUE):
  69   /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4.
  70   /// Good candidates are made of at least three hitForRec's.
  71   /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
  72
  73   TClonesArray *segments;
  74   AliMUONTrack *track;
  75   Int_t iCandidate = 0;
  76   Bool_t hitFound;
  77
  78   AliDebug(1,"Enter MakeTrackCandidates");
  79
  80   // Loop over stations(1..) 5 and 4 and make track candidates
  81   for (Int_t istat=4; istat>=3; istat--) {
  82
  83     // Make segments in the station
  84     segments = MakeSegmentsInStation(istat);
  85
  86     // Loop over segments
  87     for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++)
  88     {
  89       AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
  90
  91       // Transform segments to tracks and put them at the end of fRecTracksPtr
  92       track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack((AliMUONObjectPair*)((*segments)[iseg]));
  93       fNRecTracks++;
  94
  95       // Printout for debuging
  96       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2))
  97       {
  98         cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl;
  99         ((AliMUONTrackParam*) track->GetTrackParamAtHit()->First())->GetCovariances().Print();
 100       }
 101
 102       // Look for compatible hitForRec(s) in the other station
 103       if (AliMUONReconstructor::GetRecoParam()->MakeTrackCandidatesFast()) hitFound = FollowLinearTrackInStation(*track,7-istat);
 104       else hitFound = FollowTrackInStation(*track,7-istat);
 105
 106       // Remove track if no hit found
 107       if (!hitFound) {
 108         fRecTracksPtr->Remove(track);
 109         fNRecTracks--;
 110       }
 111
 112     }
 113
 114     // delete the array of segments
 115     delete segments;
 116   }
 117
 118   fRecTracksPtr->Compress(); // this is essential before checking tracks
 119
 120   // Keep all different tracks or only the best ones as required
 121   if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) RemoveIdenticalTracks();
 122   else RemoveDoubleTracks();
 123
 124   AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
 125
 126 }
 127
 128   //__________________________________________________________________________
 129 void AliMUONTrackReconstructor::FollowTracks()
 130 {
 131   /// Follow tracks in stations(1..) 3, 2 and 1
 132   AliDebug(1,"Enter FollowTracks");
 133
 134   AliMUONTrack *track, *nextTrack;
 135   AliMUONTrackParam *trackParam, *nextTrackParam;
 136   Int_t currentNRecTracks;
 137   Bool_t hitFound;
 138
 139   Double_t sigmaCut2 = AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
 140                        AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking();
 141
 142   for (Int_t station = 2; station >= 0; station--) {
 143
 144     // Save the actual number of reconstructed track in case of
 145     // tracks are added or suppressed during the tracking procedure
 146     // !! Do not compress fRecTracksPtr until the end of the loop over tracks !!
 147     currentNRecTracks = fNRecTracks;
 148
 149     for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) {
 150       AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1));
 151
 152       track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
 153
 154       // Fit the track:
 155       // Do not take into account the multiple scattering to speed up the fit
 156       // Calculate the track parameter covariance matrix
 157       // If "station" is station(1..) 3 then use the vertex to better constrain the fit
 158       if (station==2) {
 159         SetVertexForFit(*track);
 160         track->SetFitWithVertex(kTRUE);
 161       } else track->SetFitWithVertex(kFALSE);
 162       Fit(*track, kFALSE, kTRUE);
 163
 164       // Remove the track if the normalized chi2 is too high
 165       if (track->GetNormalizedChi2() > sigmaCut2) {
 166         fRecTracksPtr->Remove(track);
 167         fNRecTracks--;
 168         continue;
 169       }
 170
 171       // save parameters from fit into smoothed parameters to complete track afterward
 172       if (AliMUONReconstructor::GetRecoParam()->ComplementTracks()) {
 173
 174         if (station==2) { // save track parameters on stations 4 and 5
 175
 176           // extrapolate track parameters and covariances at each cluster
 177           track->UpdateCovTrackParamAtHit();
 178
 179           // save them
 180           trackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First();
 181           while (trackParam) {
 182             trackParam->SetSmoothParameters(trackParam->GetParameters());
 183             trackParam->SetSmoothCovariances(trackParam->GetCovariances());
 184             trackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->After(trackParam);
 185           }
 186
 187         } else { // or save track parameters on last station only
 188
 189           // save parameters from fit
 190           trackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First();
 191           trackParam->SetSmoothParameters(trackParam->GetParameters());
 192           trackParam->SetSmoothCovariances(trackParam->GetCovariances());
 193
 194           // save parameters extrapolated to the second chamber of the same station if it has been hit
 195           nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->After(trackParam);
 196           if (nextTrackParam->GetHitForRecPtr()->GetChamberNumber() < 2*(station+2)) {
 197
 198             // reset parameters and covariances
 199             nextTrackParam->SetParameters(trackParam->GetParameters());
 200             nextTrackParam->SetZ(trackParam->GetZ());
 201             nextTrackParam->SetCovariances(trackParam->GetCovariances());
 202
 203             // extrapolate them to the z of the corresponding cluster
 204             AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetHitForRecPtr()->GetZ());
 205
 206             // save them
 207             nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters());
 208             nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances());
 209
 210           }
 211
 212         }
 213
 214       }
 215
 216       // Printout for debuging
 217       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 218         cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl;
 219         ((AliMUONTrackParam*) track->GetTrackParamAtHit()->First())->GetCovariances().Print();
 220       }
 221
 222       // Look for compatible hitForRec in station(0..) "station"
 223       hitFound = FollowTrackInStation(*track,station);
 224
 225       // Try to recover track if required
 226       if (!hitFound && AliMUONReconstructor::GetRecoParam()->RecoverTracks()) hitFound = RecoverTrack(*track,station);
 227
 228       // remove track if no hit found
 229       if (!hitFound) {
 230         fRecTracksPtr->Remove(track);
 231         fNRecTracks--;
 232       }
 233
 234     }
 235
 236     // Compress fRecTracksPtr for the next step
 237     fRecTracksPtr->Compress();
 238
 239     // Keep only the best tracks if required
 240     if (!AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
 241
 242   }
 243
 244   // Last fit of track candidates with all station
 245   // Take into account the multiple scattering and remove bad tracks
 246   Int_t trackIndex = -1;
 247   track = (AliMUONTrack*) fRecTracksPtr->First();
 248   while (track) {
 249
 250     trackIndex++;
 251     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
 252
 253     track->SetFitWithVertex(kFALSE); // just to be sure
 254     Fit(*track, kTRUE, kTRUE);
 255
 256     // Printout for debuging
 257     if (AliLog::GetGlobalDebugLevel() >= 3) {
 258       cout << "FollowTracks: track candidate(0..) " << trackIndex << " after final fit" << endl;
 259       track->RecursiveDump();
 260     }
 261
 262     // Remove the track if the normalized chi2 is too high
 263     if (track->GetNormalizedChi2() > sigmaCut2) {
 264       fRecTracksPtr->Remove(track);
 265       fNRecTracks--;
 266     }
 267
 268     // save parameters from fit into smoothed parameters to complete track afterward
 269     if (AliMUONReconstructor::GetRecoParam()->ComplementTracks()) {
 270
 271       // save parameters from fit
 272       trackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First();
 273       trackParam->SetSmoothParameters(trackParam->GetParameters());
 274       trackParam->SetSmoothCovariances(trackParam->GetCovariances());
 275
 276       // save parameters extrapolated to the second chamber of the same station if it has been hit
 277       nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->After(trackParam);
 278       if (nextTrackParam->GetHitForRecPtr()->GetChamberNumber() < 2) {
 279
 280         // reset parameters and covariances
 281         nextTrackParam->SetParameters(trackParam->GetParameters());
 282         nextTrackParam->SetZ(trackParam->GetZ());
 283         nextTrackParam->SetCovariances(trackParam->GetCovariances());
 284
 285         // extrapolate them to the z of the corresponding cluster
 286         AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetHitForRecPtr()->GetZ());
 287
 288         // save them
 289         nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters());
 290         nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances());
 291
 292       }
 293
 294     }
 295
 296     track = nextTrack;
 297
 298   }
 299
 300   fRecTracksPtr->Compress();
 301
 302 }
 303
 304   //__________________________________________________________________________
 305 Bool_t AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack &trackCandidate, Int_t nextStation)
 306 {
 307   /// Follow trackCandidate in station(0..) nextStation and search for compatible HitForRec(s)
 308   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
 309   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
 310   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
 311   ///         Remove the obsolete "trackCandidate" at the end.
 312   /// kFALSE: add only the best hit(s) to the "trackCandidate". Try to add a couple of hits in priority.
 313   AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1));
 314
 315   // Order the chamber according to the propagation direction (tracking starts with chamber 2):
 316   // - nextStation == station(1...) 5 => forward propagation
 317   // - nextStation < station(1...) 5 => backward propagation
 318   Int_t ch1, ch2;
 319   if (nextStation==4) {
 320     ch1 = 2*nextStation+1;
 321     ch2 = 2*nextStation;
 322   } else {
 323     ch1 = 2*nextStation;
 324     ch2 = 2*nextStation+1;
 325   }
 326
 327   Double_t chi2WithOneHitForRec = 1.e10;
 328   Double_t chi2WithTwoHitForRec = 1.e10;
 329   Double_t maxChi2WithOneHitForRec = 2. * AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
 330                                           AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
 331   Double_t maxChi2WithTwoHitForRec = 4. * AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
 332                                           AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking(); // 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 (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
 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 (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
 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 (AliMUONReconstructor::GetRecoParam()->TrackAllTracks() || !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 (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
 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 (!AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
 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 = AliMUONReconstructor::GetRecoParam()->GetNonBendingVertexDispersion() *
 808                              AliMUONReconstructor::GetRecoParam()->GetNonBendingVertexDispersion();
 809   Double_t bendingReso2 = AliMUONReconstructor::GetRecoParam()->GetBendingVertexDispersion() *
 810                           AliMUONReconstructor::GetRecoParam()->GetBendingVertexDispersion();
 811   // add multiple scattering effets
 812   AliMUONTrackParam paramAtVertex(*((AliMUONTrackParam*)(trackCandidate.GetTrackParamAtHit()->First())));
 813   paramAtVertex.DeleteCovariances(); // to be sure to account only for multiple scattering
 814   AliMUONTrackExtrap::ExtrapToVertexUncorrected(&paramAtVertex,0.);
 815   const TMatrixD& kParamCov = paramAtVertex.GetCovariances();
 816   nonBendingReso2 += kParamCov(0,0);
 817   bendingReso2 += kParamCov(2,2);
 818   // Set the vertex
 819   AliMUONHitForRec vertex; // Coordinates set to (0.,0.,0.) by default
 820   vertex.SetNonBendingReso2(nonBendingReso2);
 821   vertex.SetBendingReso2(bendingReso2);
 822   trackCandidate.SetVertex(&vertex);
 823 }
 824
 825   //__________________________________________________________________________
 826 void AliMUONTrackReconstructor::Fit(AliMUONTrack &track, Bool_t includeMCS, Bool_t calcCov)
 827 {
 828   /// Fit the track "track" w/wo multiple Coulomb scattering according to "includeMCS".
 829
 830   Double_t benC, errorParam, invBenP, nonBenC, x, y;
 831   AliMUONTrackParam *trackParam;
 832   Double_t arg[1], fedm, errdef, fitFMin;
 833   Int_t npari, nparx;
 834   Int_t status, covStatus;
 835
 836   // Instantiate gMinuit if not already done
 837   if (!gMinuit) gMinuit = new TMinuit(6);
 838   // Clear MINUIT parameters
 839   gMinuit->mncler();
 840   // Give the fitted track to MINUIT
 841   gMinuit->SetObjectFit(&track);
 842   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 843     // Define print level
 844     arg[0] = 1;
 845     gMinuit->mnexcm("SET PRI", arg, 1, status);
 846     // Print covariance matrix
 847     gMinuit->mnexcm("SHO COV", arg, 0, status);
 848   } else {
 849     arg[0] = -1;
 850     gMinuit->mnexcm("SET PRI", arg, 1, status);
 851   }
 852   // No warnings
 853   gMinuit->mnexcm("SET NOW", arg, 0, status);
 854   // Define strategy
 855   //arg[0] = 2;
 856   //gMinuit->mnexcm("SET STR", arg, 1, status);
 857
 858   // set flag w/wo multiple scattering according to "includeMCS"
 859   track.SetFitWithMCS(includeMCS);
 860   if (includeMCS) {
 861     // compute hit weights only once
 862     if (!track.ComputeHitWeights()) {
 863       AliWarning("cannot take into account the multiple scattering effects");
 864       track.SetFitWithMCS(kFALSE);
 865     }
 866   }
 867
 868   // Set fitting function
 869   gMinuit->SetFCN(TrackChi2);
 870
 871   // Set fitted parameters (!! The order is very important for the covariance matrix !!)
 872   trackParam = (AliMUONTrackParam*) (track.GetTrackParamAtHit()->First());
 873   // could be tried with no limits for the search (min=max=0) ????
 874   // mandatory limits in non Bending to avoid NaN values of parameters
 875   gMinuit->mnparm(0, "X", trackParam->GetNonBendingCoor(), 0.03, -500.0, 500.0, status);
 876   gMinuit->mnparm(1, "NonBenS", trackParam->GetNonBendingSlope(), 0.001, -0.5, 0.5, status);
 877   // mandatory limits in Bending to avoid NaN values of parameters
 878   gMinuit->mnparm(2, "Y", trackParam->GetBendingCoor(), 0.10, -500.0, 500.0, status);
 879   gMinuit->mnparm(3, "BenS", trackParam->GetBendingSlope(), 0.001, -0.5, 0.5, status);
 880   gMinuit->mnparm(4, "InvBenP", trackParam->GetInverseBendingMomentum(), 0.003, -0.4, 0.4, status);
 881
 882   // minimization
 883   gMinuit->mnexcm("MIGRAD", arg, 0, status);
 884
 885   // Calculate the covariance matrix more accurately if required
 886   if (calcCov) gMinuit->mnexcm("HESSE", arg, 0, status);
 887
 888   // get results into "invBenP", "benC", "nonBenC" ("x", "y")
 889   gMinuit->GetParameter(0, x, errorParam);
 890   trackParam->SetNonBendingCoor(x);
 891   gMinuit->GetParameter(1, nonBenC, errorParam);
 892   trackParam->SetNonBendingSlope(nonBenC);
 893   gMinuit->GetParameter(2, y, errorParam);
 894   trackParam->SetBendingCoor(y);
 895   gMinuit->GetParameter(3, benC, errorParam);
 896   trackParam->SetBendingSlope(benC);
 897   gMinuit->GetParameter(4, invBenP, errorParam);
 898   trackParam->SetInverseBendingMomentum(invBenP);
 899
 900   // global result of the fit
 901   gMinuit->mnstat(fitFMin, fedm, errdef, npari, nparx, covStatus);
 902   track.SetFitFMin(fitFMin);
 903
 904   // Get the covariance matrix if required
 905   if (calcCov) {
 906     // Covariance matrix according to HESSE status
 907     // If problem then keep only the diagonal terms (variances)
 908     Double_t matrix[5][5];
 909     gMinuit->mnemat(&matrix[0][0],5);
 910     if (covStatus == 3) trackParam->SetCovariances(matrix);
 911     else trackParam->SetVariances(matrix);
 912   } else trackParam->DeleteCovariances();
 913
 914 }
 915
 916   //__________________________________________________________________________
 917 void TrackChi2(Int_t & /*nParam*/, Double_t * /*gradient*/, Double_t &chi2, Double_t *param, Int_t /*flag*/)
 918 {
 919   /// Return the "Chi2" to be minimized with Minuit for track fitting.
 920   /// Assumes that the track hits are sorted according to increasing Z.
 921   /// Track parameters at each TrackHit are updated accordingly.
 922   /// Vertex is used according to the flag "trackBeingFitted->GetFitWithVertex()".
 923   /// Multiple Coulomb scattering is taken into account according to the flag "trackBeingFitted->GetFitWithMCS()".
 924
 925   AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit();
 926   AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) trackBeingFitted->GetTrackParamAtHit()->First();
 927   Double_t dX, dY;
 928   chi2 = 0.; // initialize chi2
 929
 930   // update track parameters
 931   trackParamAtHit->SetNonBendingCoor(param[0]);
 932   trackParamAtHit->SetNonBendingSlope(param[1]);
 933   trackParamAtHit->SetBendingCoor(param[2]);
 934   trackParamAtHit->SetBendingSlope(param[3]);
 935   trackParamAtHit->SetInverseBendingMomentum(param[4]);
 936   trackBeingFitted->UpdateTrackParamAtHit();
 937
 938   // Take the vertex into account in the fit if required
 939   if (trackBeingFitted->GetFitWithVertex()) {
 940     AliMUONTrackParam paramAtVertex(*trackParamAtHit);
 941     AliMUONTrackExtrap::ExtrapToZ(&paramAtVertex, 0.);
 942     AliMUONHitForRec *vertex = trackBeingFitted->GetVertex();
 943     if (!vertex) {
 944       cout<<"Error in TrackChi2: Want to use the vertex in tracking but it has not been created!!"<<endl;
 945       exit(-1);
 946     }
 947     dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor();
 948     dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor();
 949     chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2();
 950   }
 951
 952   // compute chi2 w/wo multiple scattering
 953   chi2 += trackBeingFitted->ComputeGlobalChi2(trackBeingFitted->GetFitWithMCS());
 954
 955 }
 956
 957   //__________________________________________________________________________
 958 void AliMUONTrackReconstructor::ComplementTracks()
 959 {
 960   /// Complete tracks by adding missing clusters (if there is an overlap between
 961   /// two detection elements, the track may have two clusters in the same chamber)
 962   /// Re-fit track parameters and covariances
 963   AliDebug(1,"Enter ComplementTracks");
 964
 965   Int_t chamberId, detElemId;
 966   Double_t chi2OfHitForRec, bestChi2OfHitForRec;
 967   Bool_t foundOneHit, trackModified;
 968   AliMUONHitForRec* hitForRec;
 969   AliMUONTrackParam* trackParam, *nextTrackParam;
 970   AliMUONTrackParam copyOfTrackParam;
 971   AliMUONTrackParam trackParamAtHit;
 972   AliMUONTrackParam bestTrackParamAtHit;
 973   Double_t sigmaCut2 = AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
 974                        AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking();
 975
 976   // Remove double track to complete only "good" tracks
 977   RemoveDoubleTracks();
 978
 979   AliMUONTrack *track = (AliMUONTrack*) fRecTracksPtr->First();
 980   while (track) {
 981     trackModified = kFALSE;
 982
 983     trackParam = (AliMUONTrackParam*)track->GetTrackParamAtHit()->First();
 984     while (trackParam) {
 985       foundOneHit = kFALSE;
 986       bestChi2OfHitForRec = 2. * sigmaCut2; // 2 because 2 quantities in chi2
 987
 988       // prepare nextTrackParam before adding new cluster because of the sorting
 989       nextTrackParam = (AliMUONTrackParam*)track->GetTrackParamAtHit()->After(trackParam);
 990
 991       chamberId = trackParam->GetHitForRecPtr()->GetChamberNumber();
 992       detElemId = trackParam->GetHitForRecPtr()->GetDetElemId();
 993
 994       // recover track parameters from local fit and put them into a copy of trackParam
 995       copyOfTrackParam.SetZ(trackParam->GetZ());
 996       copyOfTrackParam.SetParameters(trackParam->GetSmoothParameters());
 997       copyOfTrackParam.SetCovariances(trackParam->GetSmoothCovariances());
 998
 999       // look for one second candidate in the same chamber
1000       for (Int_t hit = 0; hit < fNHitsForRecPerChamber[chamberId]; hit++) {
1001
1002         hitForRec = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[chamberId]+hit);
1003
1004         // look for a cluster in another detection element
1005         if (hitForRec->GetDetElemId() == detElemId) continue;
1006
1007         // try to add the current hit fast
1008         if (!TryOneHitForRecFast(copyOfTrackParam, hitForRec)) continue;
1009
1010         // try to add the current hit accurately
1011         chi2OfHitForRec = TryOneHitForRec(copyOfTrackParam, hitForRec, trackParamAtHit);
1012
1013         // if better chi2 then prepare to add this cluster to the track
1014         if (chi2OfHitForRec < bestChi2OfHitForRec) {
1015           bestChi2OfHitForRec = chi2OfHitForRec;
1016           bestTrackParamAtHit = trackParamAtHit;
1017           foundOneHit = kTRUE;
1018         }
1019
1020       }
1021
1022       // add new cluster if any
1023       if (foundOneHit) {
1024         UpdateTrack(*track,bestTrackParamAtHit);
1025         bestTrackParamAtHit.SetAloneInChamber(kFALSE);
1026         trackParam->SetAloneInChamber(kFALSE);
1027         trackModified = kTRUE;
1028       }
1029
1030       trackParam = nextTrackParam;
1031     }
1032
1033     // re-fit track parameters if needed
1034     if (trackModified) Fit(*track, kTRUE, kTRUE);
1035
1036     track = (AliMUONTrack*) fRecTracksPtr->After(track);
1037   }
1038
1039 }
1040
1041   //__________________________________________________________________________
1042 void AliMUONTrackReconstructor::ImproveTracks()
1043 {
1044   /// Improve tracks by removing clusters with local chi2 highter than the defined cut
1045   /// Recompute track parameters and covariances at the remaining clusters
1046   AliDebug(1,"Enter ImproveTracks");
1047
1048   Double_t localChi2, worstLocalChi2;
1049   Int_t worstChamber, previousChamber;
1050   AliMUONTrack *track, *nextTrack;
1051   AliMUONTrackParam *trackParamAtHit, *worstTrackParamAtHit, *previousTrackParam, *nextTrackParam;
1052   Double_t sigmaCut2 = AliMUONReconstructor::GetRecoParam()->GetSigmaCutForImprovement() *
1053                        AliMUONReconstructor::GetRecoParam()->GetSigmaCutForImprovement();
1054
1055   // Remove double track to improve only "good" tracks
1056   RemoveDoubleTracks();
1057
1058   track = (AliMUONTrack*) fRecTracksPtr->First();
1059   while (track) {
1060
1061     // prepare next track in case the actual track is suppressed
1062     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track);
1063
1064     while (!track->IsImproved()) {
1065
1066       // Update track parameters and covariances
1067       track->UpdateCovTrackParamAtHit();
1068
1069       // Compute local chi2 of each hits
1070       track->ComputeLocalChi2(kTRUE);
1071
1072       // Look for the hit to remove
1073       worstTrackParamAtHit = NULL;
1074       worstLocalChi2 = 0.;
1075       trackParamAtHit = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First();
1076       while (trackParamAtHit) {
1077
1078         // Pick up hit with the worst chi2
1079         localChi2 = trackParamAtHit->GetLocalChi2();
1080         if (localChi2 > worstLocalChi2) {
1081           worstLocalChi2 = localChi2;
1082           worstTrackParamAtHit = trackParamAtHit;
1083         }
1084
1085       trackParamAtHit = (AliMUONTrackParam*) track->GetTrackParamAtHit()->After(trackParamAtHit);
1086       }
1087
1088       // Check if bad hit found
1089       if (!worstTrackParamAtHit) {
1090         track->SetImproved(kTRUE);
1091         break;
1092       }
1093
1094       // Check whether the worst chi2 is under requirement or not
1095       if (worstLocalChi2 < 2. * sigmaCut2) { // 2 because 2 quantities in chi2
1096         track->SetImproved(kTRUE);
1097         break;
1098       }
1099
1100       // if the worst hit is not removable then remove the entire track
1101       if (!worstTrackParamAtHit->IsRemovable() && worstTrackParamAtHit->IsAloneInChamber()) {
1102         fRecTracksPtr->Remove(track);
1103         fNRecTracks--;
1104         break;
1105       }
1106
1107       // Reset the second hit in the same station as being not removable
1108       // or reset the second hit in the same chamber as being alone
1109       worstChamber = worstTrackParamAtHit->GetHitForRecPtr()->GetChamberNumber();
1110       previousTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->Before(worstTrackParamAtHit);
1111       nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->After(worstTrackParamAtHit);
1112       if (worstTrackParamAtHit->IsAloneInChamber()) { // Worst hit removable and alone in chamber
1113
1114         if (worstChamber%2 == 0) { // Modify flags in next chamber
1115
1116           nextTrackParam->SetRemovable(kFALSE);
1117           if (!nextTrackParam->IsAloneInChamber()) // Make sure both hits in second chamber are not removable anymore
1118             ((AliMUONTrackParam*) track->GetTrackParamAtHit()->After(nextTrackParam))->SetRemovable(kFALSE);
1119
1120         } else { // Modify flags in previous chamber
1121
1122           previousTrackParam->SetRemovable(kFALSE);
1123           if (!previousTrackParam->IsAloneInChamber()) // Make sure both hits in second chamber are not removable anymore
1124             ((AliMUONTrackParam*) track->GetTrackParamAtHit()->Before(previousTrackParam))->SetRemovable(kFALSE);
1125
1126         }
1127
1128       } else { // Worst hit not alone in its chamber
1129
1130         if (previousTrackParam) previousChamber = previousTrackParam->GetHitForRecPtr()->GetChamberNumber();
1131         else previousChamber = -1;
1132
1133         if (previousChamber == worstChamber) { // the second hit on the same chamber is the previous one
1134
1135           previousTrackParam->SetAloneInChamber(kTRUE);
1136           // transfert the removability to the second hit
1137           if (worstTrackParamAtHit->IsRemovable()) previousTrackParam->SetRemovable(kTRUE);
1138
1139         } else { // the second hit on the same chamber is the next one
1140
1141           nextTrackParam->SetAloneInChamber(kTRUE);
1142           // transfert the removability to the second hit
1143           if (worstTrackParamAtHit->IsRemovable()) nextTrackParam->SetRemovable(kTRUE);
1144
1145         }
1146
1147       }
1148
1149       // Remove the worst hit
1150       track->RemoveTrackParamAtHit(worstTrackParamAtHit);
1151
1152       // Re-fit the track:
1153       // Take into account the multiple scattering
1154       // Calculate the track parameter covariance matrix
1155       Fit(*track, kTRUE, kTRUE);
1156
1157       // Printout for debuging
1158       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1159         cout << "ImproveTracks: track " << fRecTracksPtr->IndexOf(track)+1 << " improved " << endl;
1160       }
1161
1162     }
1163
1164     track = nextTrack;
1165   }
1166
1167   // compress the array in case of some tracks have been removed
1168   fRecTracksPtr->Compress();
1169
1170 }
1171
1172   //__________________________________________________________________________
1173 void AliMUONTrackReconstructor::Finalize()
1174 {
1175   /// Fill AliMUONTrack's fHitForRecAtHit array
1176   /// Recompute track parameters and covariances at each attached cluster from those at the first one
1177
1178   AliMUONTrack *track;
1179   AliMUONTrackParam *trackParamAtHit;
1180
1181   track = (AliMUONTrack*) fRecTracksPtr->First();
1182   while (track) {
1183
1184     // update track parameters if not already done
1185     if (!track->IsImproved()) track->UpdateCovTrackParamAtHit();
1186
1187     trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
1188     while (trackParamAtHit) {
1189
1190       // update array of track hit
1191       track->AddHitForRecAtHit(trackParamAtHit->GetHitForRecPtr());
1192
1193       trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->After(trackParamAtHit));
1194     }
1195
1196     track = (AliMUONTrack*) fRecTracksPtr->After(track);
1197
1198   }
1199
1200 }
1201
1202