MUON/AliMUONVTrackReconstructor.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 AliMUONVTrackReconstructor
  20 /// Virtual MUON track reconstructor in ALICE (class renamed from AliMUONEventReconstructor)
  21 ///
  22 /// This class contains as data a pointer to the array of reconstructed tracks
  23 ///
  24 /// It contains as methods, among others:
  25 /// * EventReconstruct to build the muon tracks
  26 /// * EventReconstructTrigger to build the trigger tracks
  27 /// * ValidateTracksWithTrigger to match tracker/trigger tracks
  28 ///
  29 /// Several options and adjustable parameters are available for both KALMAN and ORIGINAL
  30 /// tracking algorithms. They can be changed through the AliMUONRecoParam object
  31 /// set in the reconstruction macro or read from the CDB
  32 /// (see methods in AliMUONRecoParam.h file for details)
  33 ///
  34 /// Main parameters and options are:
  35 /// - *fgkSigmaToCutForTracking* : quality cut used to select new clusters to be
  36 ///   attached to the track candidate and to select good tracks.
  37 /// - *fgkMakeTrackCandidatesFast* : if this flag is set to 'true', the track candidates
  38 ///   are made assuming linear propagation between stations 4 and 5.
  39 /// - *fgkTrackAllTracks* : according to the value of this flag, in case that several
  40 ///   new clusters pass the quality cut, either we consider all the possibilities
  41 ///   (duplicating tracks) or we attach only the best cluster.
  42 /// - *fgkRecoverTracks* : if this flag is set to 'true', we try to recover the tracks
  43 ///   lost during the tracking by removing the worst of the 2 clusters attached in the
  44 ///   previous station.
  45 /// - *fgkImproveTracks* : if this flag is set to 'true', we try to improve the quality
  46 ///   of the tracks at the end of the tracking by removing clusters that do not pass
  47 ///   new quality cut (the track is removed is it does not contain enough cluster anymore).
  48 /// - *fgkComplementTracks* : if this flag is set to 'true', we try to improve the quality
  49 ///   of the tracks at the end of the tracking by adding potentially missing clusters
  50 ///   (we may have 2 clusters in the same chamber because of the overlapping of detection
  51 ///   elements, which is not handle by the tracking algorithm).
  52 /// - *fgkSigmaToCutForImprovement* : quality cut used when we try to improve the
  53 ///   quality of the tracks.
  54 ///
  55 ///  \author Philippe Pillot
  56 //-----------------------------------------------------------------------------
  57
  58 #include "AliMUONVTrackReconstructor.h"
  59
  60 #include "AliMUONConstants.h"
  61 #include "AliMUONObjectPair.h"
  62 #include "AliMUONTriggerTrack.h"
  63 #include "AliMUONTriggerCircuit.h"
  64 #include "AliMUONLocalTrigger.h"
  65 #include "AliMUONGlobalTrigger.h"
  66 #include "AliMUONTrack.h"
  67 #include "AliMUONTrackParam.h"
  68 #include "AliMUONTrackExtrap.h"
  69 #include "AliMUONTrackHitPattern.h"
  70 #include "AliMUONVTrackStore.h"
  71 #include "AliMUONVClusterStore.h"
  72 #include "AliMUONVCluster.h"
  73 #include "AliMUONVClusterServer.h"
  74 #include "AliMUONVTriggerStore.h"
  75 #include "AliMUONVTriggerTrackStore.h"
  76 #include "AliMUONRecoParam.h"
  77
  78 #include "AliMpDEManager.h"
  79 #include "AliMpArea.h"
  80
  81 #include "AliMpDDLStore.h"
  82 #include "AliMpVSegmentation.h"
  83 #include "AliMpSegmentation.h"
  84 #include "AliMpPad.h"
  85
  86 #include "AliLog.h"
  87 #include "AliCodeTimer.h"
  88 #include "AliTracker.h"
  89
  90 #include <TClonesArray.h>
  91 #include <TMath.h>
  92 #include <TMatrixD.h>
  93 #include <TVector2.h>
  94
  95 #include <Riostream.h>
  96
  97 /// \cond CLASSIMP
  98 ClassImp(AliMUONVTrackReconstructor) // Class implementation in ROOT context
  99 /// \endcond
 100
 101   //__________________________________________________________________________
 102 AliMUONVTrackReconstructor::AliMUONVTrackReconstructor(const AliMUONRecoParam* recoParam,
 103                                                        AliMUONVClusterServer* clusterServer)
 104 : TObject(),
 105 fRecTracksPtr(0x0),
 106 fNRecTracks(0),
 107 fClusterServer(clusterServer),
 108 fkRecoParam(recoParam),
 109 fMaxMCSAngle2(0x0)
 110 {
 111   /// Constructor for class AliMUONVTrackReconstructor
 112   /// WARNING: if clusterServer=0x0, no clusterization will be possible at this level
 113
 114   // Memory allocation for the TClonesArray of reconstructed tracks
 115   fRecTracksPtr = new TClonesArray("AliMUONTrack", 100);
 116
 117   // set the magnetic field for track extrapolations
 118   AliMUONTrackExtrap::SetField();
 119
 120   // set the maximum MCS angle in chamber from the minimum acceptable momentum
 121   AliMUONTrackParam param;
 122   Double_t inverseBendingP = (GetRecoParam()->GetMinBendingMomentum() > 0.) ? 1./GetRecoParam()->GetMinBendingMomentum() : 1.;
 123   param.SetInverseBendingMomentum(inverseBendingP);
 124   fMaxMCSAngle2 = new Double_t [AliMUONConstants::NTrackingCh()];
 125   for (Int_t iCh=0; iCh<AliMUONConstants::NTrackingCh(); iCh++)
 126     fMaxMCSAngle2[iCh] = AliMUONTrackExtrap::GetMCSAngle2(param, AliMUONConstants::ChamberThicknessInX0(iCh), 1.);
 127
 128 }
 129
 130   //__________________________________________________________________________
 131 AliMUONVTrackReconstructor::~AliMUONVTrackReconstructor()
 132 {
 133   /// Destructor for class AliMUONVTrackReconstructor
 134   delete fRecTracksPtr;
 135   delete[] fMaxMCSAngle2;
 136 }
 137
 138   //__________________________________________________________________________
 139 void AliMUONVTrackReconstructor::ResetTracks()
 140 {
 141   /// To reset the TClonesArray of reconstructed tracks
 142   if (fRecTracksPtr) fRecTracksPtr->Clear("C");
 143   fNRecTracks = 0;
 144   return;
 145 }
 146
 147   //__________________________________________________________________________
 148 void AliMUONVTrackReconstructor::EventReconstruct(AliMUONVClusterStore& clusterStore, AliMUONVTrackStore& trackStore)
 149 {
 150   /// To reconstruct one event
 151   AliDebug(1,"");
 152   AliCodeTimerAuto("",0);
 153
 154   // Reset array of tracks
 155   ResetTracks();
 156
 157   // Look for candidates from clusters in stations(1..) 4 and 5 (abort in case of failure)
 158   if (!MakeTrackCandidates(clusterStore)) return;
 159
 160   // Look for extra candidates from clusters in stations(1..) 4 and 5 (abort in case of failure)
 161   if (GetRecoParam()->MakeMoreTrackCandidates()) {
 162     if (!MakeMoreTrackCandidates(clusterStore)) return;
 163   }
 164
 165   // Stop tracking if no candidate found
 166   if (fRecTracksPtr->GetEntriesFast() == 0) return;
 167
 168   // Follow tracks in stations(1..) 3, 2 and 1 (abort in case of failure)
 169   if (!FollowTracks(clusterStore)) return;
 170
 171   // Complement the reconstructed tracks
 172   if (GetRecoParam()->ComplementTracks()) {
 173     if (ComplementTracks(clusterStore)) RemoveIdenticalTracks();
 174   }
 175
 176   // Improve the reconstructed tracks
 177   if (GetRecoParam()->ImproveTracks()) ImproveTracks();
 178
 179   // Remove connected tracks
 180   if (GetRecoParam()->RemoveConnectedTracksInSt12()) RemoveConnectedTracks(kFALSE);
 181   else RemoveConnectedTracks(kTRUE);
 182
 183   // Fill AliMUONTrack data members
 184   Finalize();
 185
 186   // Add tracks to MUON data container
 187   for (Int_t i=0; i<fNRecTracks; ++i)
 188   {
 189     AliMUONTrack * track = (AliMUONTrack*) fRecTracksPtr->At(i);
 190     if (track->GetGlobalChi2() < AliMUONTrack::MaxChi2()) {
 191       track->SetUniqueID(i+1);
 192       trackStore.Add(*track);
 193     } else AliWarning("problem occur somewhere during track refitting --> discard track");
 194   }
 195 }
 196
 197 //__________________________________________________________________________
 198 Bool_t AliMUONVTrackReconstructor::IsAcceptable(AliMUONTrackParam &trackParam)
 199 {
 200   /// Return kTRUE if the track is within given limits on momentum/angle/origin
 201
 202   const TMatrixD& kParamCov = trackParam.GetCovariances();
 203   Int_t chamber = trackParam.GetClusterPtr()->GetChamberId();
 204   Double_t z = trackParam.GetZ();
 205   Double_t sigmaCut = GetRecoParam()->GetSigmaCutForTracking();
 206
 207   // MCS dipersion
 208   Double_t angleMCS2 = 0.;
 209   Double_t impactMCS2 = 0.;
 210   if (AliMUONTrackExtrap::IsFieldON() && chamber < 6) {
 211
 212     // track momentum is known
 213     for (Int_t iCh=0; iCh<=chamber; iCh++) {
 214       Double_t localMCS2 = AliMUONTrackExtrap::GetMCSAngle2(trackParam, AliMUONConstants::ChamberThicknessInX0(iCh), 1.);
 215       angleMCS2 += localMCS2;
 216       impactMCS2 += AliMUONConstants::DefaultChamberZ(chamber) * AliMUONConstants::DefaultChamberZ(chamber) * localMCS2;
 217     }
 218
 219   } else {
 220
 221     // track momentum is unknown
 222     for (Int_t iCh=0; iCh<=chamber; iCh++) {
 223       angleMCS2 += fMaxMCSAngle2[iCh];
 224       impactMCS2 += AliMUONConstants::DefaultChamberZ(chamber) * AliMUONConstants::DefaultChamberZ(chamber) * fMaxMCSAngle2[iCh];
 225     }
 226
 227   }
 228
 229   // ------ track selection in non bending direction ------
 230   if (GetRecoParam()->SelectOnTrackSlope()) {
 231
 232     // check if non bending slope is within tolerances
 233     Double_t nonBendingSlopeErr = TMath::Sqrt(kParamCov(1,1) + angleMCS2);
 234     if ((TMath::Abs(trackParam.GetNonBendingSlope()) - sigmaCut * nonBendingSlopeErr) > GetRecoParam()->GetMaxNonBendingSlope()) return kFALSE;
 235
 236   } else {
 237
 238     // or check if non bending impact parameter is within tolerances
 239     Double_t nonBendingImpactParam = TMath::Abs(trackParam.GetNonBendingCoor() - z * trackParam.GetNonBendingSlope());
 240     Double_t nonBendingImpactParamErr = TMath::Sqrt(kParamCov(0,0) + z * z * kParamCov(1,1) - 2. * z * kParamCov(0,1) + impactMCS2);
 241     if ((nonBendingImpactParam - sigmaCut * nonBendingImpactParamErr) > (3. * GetRecoParam()->GetNonBendingVertexDispersion())) return kFALSE;
 242
 243   }
 244
 245   // ------ track selection in bending direction ------
 246   if (AliMUONTrackExtrap::IsFieldON()) { // depending whether the field is ON or OFF
 247
 248     // check if bending momentum is within tolerances
 249     Double_t bendingMomentum = TMath::Abs(1. / trackParam.GetInverseBendingMomentum());
 250     Double_t bendingMomentumErr = TMath::Sqrt(kParamCov(4,4)) * bendingMomentum * bendingMomentum;
 251     if (chamber < 6 && (bendingMomentum + sigmaCut * bendingMomentumErr) < GetRecoParam()->GetMinBendingMomentum()) return kFALSE;
 252     else if ((bendingMomentum + 3. * bendingMomentumErr) < GetRecoParam()->GetMinBendingMomentum()) return kFALSE;
 253
 254   } else {
 255
 256     if (GetRecoParam()->SelectOnTrackSlope()) {
 257
 258       // check if bending slope is within tolerances
 259       Double_t bendingSlopeErr = TMath::Sqrt(kParamCov(3,3) + angleMCS2);
 260       if ((TMath::Abs(trackParam.GetBendingSlope()) - sigmaCut * bendingSlopeErr) > GetRecoParam()->GetMaxBendingSlope()) return kFALSE;
 261
 262     } else {
 263
 264       // or check if bending impact parameter is within tolerances
 265       Double_t bendingImpactParam = TMath::Abs(trackParam.GetBendingCoor() - z * trackParam.GetBendingSlope());
 266       Double_t bendingImpactParamErr = TMath::Sqrt(kParamCov(2,2) + z * z * kParamCov(3,3) - 2. * z * kParamCov(2,3) + impactMCS2);
 267       if ((bendingImpactParam - sigmaCut * bendingImpactParamErr) > (3. * GetRecoParam()->GetBendingVertexDispersion())) return kFALSE;
 268
 269     }
 270
 271   }
 272
 273   return kTRUE;
 274
 275 }
 276
 277 //__________________________________________________________________________
 278 TClonesArray* AliMUONVTrackReconstructor::MakeSegmentsBetweenChambers(const AliMUONVClusterStore& clusterStore, Int_t ch1, Int_t ch2)
 279 {
 280   /// To make the list of segments from the list of clusters in the 2 given chambers.
 281   /// Return a new TClonesArray of segments.
 282   /// It is the responsibility of the user to delete it afterward.
 283   AliDebug(1,Form("Enter MakeSegmentsBetweenChambers (1..) %d-%d", ch1+1, ch2+1));
 284   AliCodeTimerAuto("",0);
 285
 286   AliMUONVCluster *cluster1, *cluster2;
 287   AliMUONObjectPair *segment;
 288   Double_t z1 = 0., z2 = 0., dZ = 0.;
 289   Double_t nonBendingSlope = 0., nonBendingSlopeErr = 0., nonBendingImpactParam = 0., nonBendingImpactParamErr = 0.;
 290   Double_t bendingSlope = 0., bendingSlopeErr = 0., bendingImpactParam = 0., bendingImpactParamErr = 0., bendingImpactParamErr2 = 0.;
 291   Double_t bendingMomentum = 0., bendingMomentumErr = 0.;
 292   Double_t bendingVertexDispersion2 = GetRecoParam()->GetBendingVertexDispersion() * GetRecoParam()->GetBendingVertexDispersion();
 293   Double_t angleMCS2 = 0.; // maximum angular dispersion**2 due to MCS in chamber
 294   Double_t impactMCS2 = 0.; // maximum impact parameter dispersion**2 due to MCS in chamber
 295   for (Int_t iCh=0; iCh<=ch1; iCh++) {
 296     angleMCS2 += fMaxMCSAngle2[iCh];
 297     impactMCS2 += AliMUONConstants::DefaultChamberZ(iCh) * AliMUONConstants::DefaultChamberZ(iCh) * fMaxMCSAngle2[iCh];
 298   }
 299   Double_t sigmaCut = GetRecoParam()->GetSigmaCutForTracking();
 300
 301   // Create iterators to loop over clusters in both chambers
 302   TIter nextInCh1(clusterStore.CreateChamberIterator(ch1,ch1));
 303   TIter nextInCh2(clusterStore.CreateChamberIterator(ch2,ch2));
 304
 305   // list of segments
 306   TClonesArray *segments = new TClonesArray("AliMUONObjectPair", 100);
 307
 308   // Loop over clusters in the first chamber of the station
 309   while ( ( cluster1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
 310     z1 = cluster1->GetZ();
 311
 312     // reset cluster iterator of chamber 2
 313     nextInCh2.Reset();
 314
 315     // Loop over clusters in the second chamber of the station
 316     while ( ( cluster2 = static_cast<AliMUONVCluster*>(nextInCh2()) ) ) {
 317       z2 = cluster2->GetZ();
 318       dZ = z1 - z2;
 319
 320       // ------ track selection in non bending direction ------
 321       nonBendingSlope = (cluster1->GetX() - cluster2->GetX()) / dZ;
 322       if (GetRecoParam()->SelectOnTrackSlope()) {
 323
 324         // check if non bending slope is within tolerances
 325         nonBendingSlopeErr = TMath::Sqrt((cluster1->GetErrX2() + cluster2->GetErrX2()) / dZ / dZ + angleMCS2);
 326         if ((TMath::Abs(nonBendingSlope) - sigmaCut * nonBendingSlopeErr) > GetRecoParam()->GetMaxNonBendingSlope()) continue;
 327
 328       } else {
 329
 330         // or check if non bending impact parameter is within tolerances
 331         nonBendingImpactParam = TMath::Abs(cluster1->GetX() - cluster1->GetZ() * nonBendingSlope);
 332         nonBendingImpactParamErr = TMath::Sqrt((z1 * z1 * cluster2->GetErrX2() + z2 * z2 * cluster1->GetErrX2()) / dZ / dZ + impactMCS2);
 333         if ((nonBendingImpactParam - sigmaCut * nonBendingImpactParamErr) > (3. * GetRecoParam()->GetNonBendingVertexDispersion())) continue;
 334
 335       }
 336
 337       // ------ track selection in bending direction ------
 338       bendingSlope = (cluster1->GetY() - cluster2->GetY()) / dZ;
 339       if (AliMUONTrackExtrap::IsFieldON()) { // depending whether the field is ON or OFF
 340
 341         // check if bending momentum is within tolerances
 342         bendingImpactParam = cluster1->GetY() - cluster1->GetZ() * bendingSlope;
 343         bendingImpactParamErr2 = (z1 * z1 * cluster2->GetErrY2() + z2 * z2 * cluster1->GetErrY2()) / dZ / dZ + impactMCS2;
 344         bendingMomentum = TMath::Abs(AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpactParam));
 345         bendingMomentumErr = TMath::Sqrt((bendingVertexDispersion2 + bendingImpactParamErr2) /
 346                                          bendingImpactParam / bendingImpactParam + 0.01) * bendingMomentum;
 347         if ((bendingMomentum + 3. * bendingMomentumErr) < GetRecoParam()->GetMinBendingMomentum()) continue;
 348
 349       } else {
 350
 351         if (GetRecoParam()->SelectOnTrackSlope()) {
 352
 353           // check if bending slope is within tolerances
 354           bendingSlopeErr = TMath::Sqrt((cluster1->GetErrY2() + cluster2->GetErrY2()) / dZ / dZ + angleMCS2);
 355           if ((TMath::Abs(bendingSlope) - sigmaCut * bendingSlopeErr) > GetRecoParam()->GetMaxBendingSlope()) continue;
 356
 357         } else {
 358
 359           // or check if bending impact parameter is within tolerances
 360           bendingImpactParam = TMath::Abs(cluster1->GetY() - cluster1->GetZ() * bendingSlope);
 361           bendingImpactParamErr = TMath::Sqrt((z1 * z1 * cluster2->GetErrY2() + z2 * z2 * cluster1->GetErrY2()) / dZ / dZ + impactMCS2);
 362           if ((bendingImpactParam - sigmaCut * bendingImpactParamErr) > (3. * GetRecoParam()->GetBendingVertexDispersion())) continue;
 363
 364         }
 365
 366       }
 367
 368       // make new segment
 369       segment = new ((*segments)[segments->GetLast()+1]) AliMUONObjectPair(cluster1, cluster2, kFALSE, kFALSE);
 370
 371       // Printout for debug
 372       if (AliLog::GetGlobalDebugLevel() > 1) {
 373         cout << "segmentIndex(0...): " << segments->GetLast() << endl;
 374         segment->Dump();
 375         cout << "Cluster in first chamber" << endl;
 376         cluster1->Print();
 377         cout << "Cluster in second chamber" << endl;
 378         cluster2->Print();
 379       }
 380
 381     }
 382
 383   }
 384
 385   // Printout for debug
 386   AliDebug(1,Form("chambers%d-%d: NSegments =  %d ", ch1+1, ch2+1, segments->GetEntriesFast()));
 387
 388   return segments;
 389 }
 390
 391   //__________________________________________________________________________
 392 void AliMUONVTrackReconstructor::RemoveUsedSegments(TClonesArray& segments)
 393 {
 394   /// To remove pairs of clusters already attached to a track
 395   AliDebug(1,"Enter RemoveUsedSegments");
 396   Int_t nSegments = segments.GetEntriesFast();
 397   Int_t nTracks = fRecTracksPtr->GetEntriesFast();
 398   AliMUONObjectPair *segment;
 399   AliMUONTrack *track;
 400   AliMUONVCluster *cluster, *cluster1, *cluster2;
 401   Bool_t foundCluster1, foundCluster2, removeSegment;
 402
 403   // Loop over segments
 404   for (Int_t iSegment=0; iSegment<nSegments; iSegment++) {
 405     segment = (AliMUONObjectPair*) segments.UncheckedAt(iSegment);
 406
 407     cluster1 = (AliMUONVCluster*) segment->First();
 408     cluster2 = (AliMUONVCluster*) segment->Second();
 409     removeSegment = kFALSE;
 410
 411     // Loop over tracks
 412     for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 413       track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iTrack);
 414
 415       // skip empty slot
 416       if (!track) continue;
 417
 418       foundCluster1 = kFALSE;
 419       foundCluster2 = kFALSE;
 420
 421       // Loop over clusters
 422       Int_t nClusters = track->GetNClusters();
 423       for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
 424         cluster = ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iCluster))->GetClusterPtr();
 425
 426         // check if both clusters are in that track
 427         if (cluster == cluster1) foundCluster1 = kTRUE;
 428         else if (cluster == cluster2) foundCluster2 = kTRUE;
 429
 430         if (foundCluster1 && foundCluster2) {
 431           removeSegment = kTRUE;
 432           break;
 433         }
 434
 435       }
 436
 437       if (removeSegment) break;
 438
 439     }
 440
 441     if (removeSegment) segments.RemoveAt(iSegment);
 442
 443   }
 444
 445   segments.Compress();
 446
 447   // Printout for debug
 448   AliDebug(1,Form("NSegments =  %d ", segments.GetEntriesFast()));
 449 }
 450
 451   //__________________________________________________________________________
 452 void AliMUONVTrackReconstructor::RemoveIdenticalTracks()
 453 {
 454   /// To remove identical tracks:
 455   /// Tracks are considered identical if they have all their clusters in common.
 456   /// One keeps the track with the larger number of clusters if need be
 457   AliMUONTrack *track1, *track2;
 458   Int_t nTracks = fRecTracksPtr->GetEntriesFast();
 459   Int_t clustersInCommon, nClusters1, nClusters2;
 460   // Loop over first track of the pair
 461   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
 462     track1 = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iTrack1);
 463     // skip empty slot
 464     if (!track1) continue;
 465     nClusters1 = track1->GetNClusters();
 466     // Loop over second track of the pair
 467     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
 468       track2 = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iTrack2);
 469       // skip empty slot
 470       if (!track2) continue;
 471       nClusters2 = track2->GetNClusters();
 472       // number of clusters in common between two tracks
 473       clustersInCommon = track1->ClustersInCommon(track2);
 474       // check for identical tracks
 475       if ((clustersInCommon == nClusters1) || (clustersInCommon == nClusters2)) {
 476         // decide which track to remove
 477         if (nClusters2 > nClusters1) {
 478           // remove track1 and continue the first loop with the track next to track1
 479           fRecTracksPtr->RemoveAt(iTrack1);
 480           fNRecTracks--;
 481           break;
 482         } else {
 483           // remove track2 and continue the second loop with the track next to track2
 484           fRecTracksPtr->RemoveAt(iTrack2);
 485           fNRecTracks--;
 486         }
 487       }
 488     } // track2
 489   } // track1
 490   fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
 491 }
 492
 493   //__________________________________________________________________________
 494 void AliMUONVTrackReconstructor::RemoveDoubleTracks()
 495 {
 496   /// To remove double tracks:
 497   /// Tracks are considered identical if more than half of the clusters of the track
 498   /// which has the smaller number of clusters are in common with the other track.
 499   /// Among two identical tracks, one keeps the track with the larger number of clusters
 500   /// or, if these numbers are equal, the track with the minimum chi2.
 501   AliMUONTrack *track1, *track2;
 502   Int_t nTracks = fRecTracksPtr->GetEntriesFast();
 503   Int_t clustersInCommon2, nClusters1, nClusters2;
 504   // Loop over first track of the pair
 505   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
 506     track1 = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iTrack1);
 507     // skip empty slot
 508     if (!track1) continue;
 509     nClusters1 = track1->GetNClusters();
 510     // Loop over second track of the pair
 511     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
 512       track2 = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iTrack2);
 513       // skip empty slot
 514       if (!track2) continue;
 515       nClusters2 = track2->GetNClusters();
 516       // number of clusters in common between two tracks
 517       clustersInCommon2 = 2 * track1->ClustersInCommon(track2);
 518       // check for identical tracks
 519       if (clustersInCommon2 > nClusters1 || clustersInCommon2 > nClusters2) {
 520         // decide which track to remove
 521         if ((nClusters1 > nClusters2) || ((nClusters1 == nClusters2) && (track1->GetGlobalChi2() <= track2->GetGlobalChi2()))) {
 522           // remove track2 and continue the second loop with the track next to track2
 523           fRecTracksPtr->RemoveAt(iTrack2);
 524           fNRecTracks--;
 525         } else {
 526           // else remove track1 and continue the first loop with the track next to track1
 527           fRecTracksPtr->RemoveAt(iTrack1);
 528           fNRecTracks--;
 529           break;
 530         }
 531       }
 532     } // track2
 533   } // track1
 534   fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
 535 }
 536
 537   //__________________________________________________________________________
 538 void AliMUONVTrackReconstructor::RemoveConnectedTracks(Bool_t inSt345)
 539 {
 540   /// To remove double tracks:
 541   /// Tracks are considered identical if they share 1 cluster or more.
 542   /// If inSt345=kTRUE only stations 3, 4 and 5 are considered.
 543   /// Among two identical tracks, one keeps the track with the larger number of clusters
 544   /// or, if these numbers are equal, the track with the minimum chi2.
 545   AliMUONTrack *track1, *track2, *trackToRemove;
 546   Int_t clustersInCommon, nClusters1, nClusters2;
 547   Bool_t removedTrack1;
 548   // Loop over first track of the pair
 549   track1 = (AliMUONTrack*) fRecTracksPtr->First();
 550   while (track1) {
 551     removedTrack1 = kFALSE;
 552     nClusters1 = track1->GetNClusters();
 553     // Loop over second track of the pair
 554     track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
 555     while (track2) {
 556       nClusters2 = track2->GetNClusters();
 557       // number of clusters in common between two tracks
 558       if (inSt345) clustersInCommon = track1->ClustersInCommonInSt345(track2);
 559       else clustersInCommon = track1->ClustersInCommon(track2);
 560       // check for identical tracks
 561       if (clustersInCommon > 0) {
 562         // decide which track to remove
 563         if ((nClusters1 > nClusters2) || ((nClusters1 == nClusters2) && (track1->GetGlobalChi2() <= track2->GetGlobalChi2()))) {
 564           // remove track2 and continue the second loop with the track next to track2
 565           trackToRemove = track2;
 566           track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
 567           fRecTracksPtr->Remove(trackToRemove);
 568           fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
 569           fNRecTracks--;
 570         } else {
 571           // else remove track1 and continue the first loop with the track next to track1
 572           trackToRemove = track1;
 573           track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
 574           fRecTracksPtr->Remove(trackToRemove);
 575           fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
 576           fNRecTracks--;
 577           removedTrack1 = kTRUE;
 578           break;
 579         }
 580       } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
 581     } // track2
 582     if (removedTrack1) continue;
 583     track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
 584   } // track1
 585   return;
 586 }
 587
 588   //__________________________________________________________________________
 589 void AliMUONVTrackReconstructor::AskForNewClustersInChamber(const AliMUONTrackParam &trackParam,
 590                                                             AliMUONVClusterStore& clusterStore, Int_t chamber)
 591 {
 592   /// Ask the clustering to reconstruct new clusters around the track candidate position
 593
 594   // check if the current chamber is useable
 595   if (!fClusterServer || !GetRecoParam()->UseChamber(chamber)) return;
 596
 597   // maximum distance between the center of the chamber and a detection element
 598   // (accounting for the inclination of the chamber)
 599   static const Double_t kMaxDZ = 15.; // 15 cm
 600
 601   // extrapolate track parameters to the chamber
 602   AliMUONTrackParam extrapTrackParam(trackParam);
 603   if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(chamber))) return;
 604
 605   // build the searching area using the track and chamber resolutions and the maximum-distance-to-track value
 606   const TMatrixD& kParamCov = extrapTrackParam.GetCovariances();
 607   Double_t errX2 = kParamCov(0,0) + kMaxDZ * kMaxDZ * kParamCov(1,1) + 2. * kMaxDZ * TMath::Abs(kParamCov(0,1)) +
 608                    GetRecoParam()->GetDefaultNonBendingReso(chamber) * GetRecoParam()->GetDefaultNonBendingReso(chamber);
 609   Double_t errY2 = kParamCov(2,2) + kMaxDZ * kMaxDZ * kParamCov(3,3) + 2. * kMaxDZ * TMath::Abs(kParamCov(2,3)) +
 610                    GetRecoParam()->GetDefaultBendingReso(chamber) * GetRecoParam()->GetDefaultBendingReso(chamber);
 611   Double_t dX = TMath::Abs(trackParam.GetNonBendingSlope()) * kMaxDZ +
 612                 GetRecoParam()->GetMaxNonBendingDistanceToTrack() +
 613                 GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errX2);
 614   Double_t dY = TMath::Abs(trackParam.GetBendingSlope()) * kMaxDZ +
 615                 GetRecoParam()->GetMaxBendingDistanceToTrack() +
 616                 GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errY2);
 617   AliMpArea area(extrapTrackParam.GetNonBendingCoor(),
 618                  extrapTrackParam.GetBendingCoor(),
 619                  dX, dY);
 620
 621   // ask to cluterize in the given area of the given chamber
 622   fClusterServer->Clusterize(chamber, clusterStore, area, GetRecoParam());
 623
 624 }
 625
 626   //__________________________________________________________________________
 627 void AliMUONVTrackReconstructor::AskForNewClustersInStation(const AliMUONTrackParam &trackParam,
 628                                                             AliMUONVClusterStore& clusterStore, Int_t station)
 629 {
 630   /// Ask the clustering to reconstruct new clusters around the track candidate position
 631   /// in the 2 chambers of the given station
 632   AskForNewClustersInChamber(trackParam, clusterStore, 2*station+1);
 633   AskForNewClustersInChamber(trackParam, clusterStore, 2*station);
 634 }
 635
 636   //__________________________________________________________________________
 637 Double_t AliMUONVTrackReconstructor::TryOneCluster(const AliMUONTrackParam &trackParam, AliMUONVCluster* cluster,
 638                                                    AliMUONTrackParam &trackParamAtCluster, Bool_t updatePropagator)
 639 {
 640 /// Test the compatibility between the track and the cluster (using trackParam's covariance matrix):
 641 /// return the corresponding Chi2
 642 /// return trackParamAtCluster
 643
 644   // extrapolate track parameters and covariances at the z position of the tested cluster
 645   // and set pointer to cluster into trackParamAtCluster
 646   trackParamAtCluster = trackParam;
 647   trackParamAtCluster.SetClusterPtr(cluster);
 648   if (!AliMUONTrackExtrap::ExtrapToZCov(&trackParamAtCluster, cluster->GetZ(), updatePropagator))
 649     return 2.*AliMUONTrack::MaxChi2();
 650
 651   // Set differences between trackParam and cluster in the bending and non bending directions
 652   Double_t dX = cluster->GetX() - trackParamAtCluster.GetNonBendingCoor();
 653   Double_t dY = cluster->GetY() - trackParamAtCluster.GetBendingCoor();
 654
 655   // Calculate errors and covariances
 656   const TMatrixD& kParamCov = trackParamAtCluster.GetCovariances();
 657   Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2();
 658   Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2();
 659   Double_t covXY   = kParamCov(0,2);
 660   Double_t det     = sigmaX2 * sigmaY2 - covXY * covXY;
 661
 662   // Compute chi2
 663   if (det == 0.) return 2.*AliMUONTrack::MaxChi2();
 664   return (dX * dX * sigmaY2 + dY * dY * sigmaX2 - 2. * dX * dY * covXY) / det;
 665
 666 }
 667
 668   //__________________________________________________________________________
 669 Bool_t AliMUONVTrackReconstructor::TryOneClusterFast(const AliMUONTrackParam &trackParam, const AliMUONVCluster* cluster)
 670 {
 671 /// Test the compatibility between the track and the cluster
 672 /// given the track and cluster resolutions + the maximum-distance-to-track value
 673 /// and assuming linear propagation of the track:
 674 /// return kTRUE if they are compatibles
 675
 676   Double_t dZ = cluster->GetZ() - trackParam.GetZ();
 677   Double_t dX = cluster->GetX() - (trackParam.GetNonBendingCoor() + trackParam.GetNonBendingSlope() * dZ);
 678   Double_t dY = cluster->GetY() - (trackParam.GetBendingCoor() + trackParam.GetBendingSlope() * dZ);
 679   const TMatrixD& kParamCov = trackParam.GetCovariances();
 680   Double_t errX2 = kParamCov(0,0) + dZ * dZ * kParamCov(1,1) + 2. * dZ * kParamCov(0,1) + cluster->GetErrX2();
 681   Double_t errY2 = kParamCov(2,2) + dZ * dZ * kParamCov(3,3) + 2. * dZ * kParamCov(2,3) + cluster->GetErrY2();
 682
 683   Double_t dXmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errX2) +
 684                    GetRecoParam()->GetMaxNonBendingDistanceToTrack();
 685   Double_t dYmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errY2) +
 686                    GetRecoParam()->GetMaxBendingDistanceToTrack();
 687
 688   if (TMath::Abs(dX) > dXmax || TMath::Abs(dY) > dYmax) return kFALSE;
 689
 690   return kTRUE;
 691
 692 }
 693
 694   //__________________________________________________________________________
 695 Double_t AliMUONVTrackReconstructor::TryTwoClustersFast(const AliMUONTrackParam &trackParamAtCluster1, AliMUONVCluster* cluster2,
 696                                                         AliMUONTrackParam &trackParamAtCluster2)
 697 {
 698 /// Test the compatibility between the track and the 2 clusters together (using trackParam's covariance matrix)
 699 /// assuming linear propagation between the two clusters:
 700 /// return the corresponding Chi2 accounting for covariances between the 2 clusters
 701 /// return trackParamAtCluster2
 702
 703   // extrapolate linearly track parameters and covariances at the z position of the second cluster
 704   trackParamAtCluster2 = trackParamAtCluster1;
 705   AliMUONTrackExtrap::LinearExtrapToZCov(&trackParamAtCluster2, cluster2->GetZ());
 706
 707   // set pointer to cluster2 into trackParamAtCluster2
 708   trackParamAtCluster2.SetClusterPtr(cluster2);
 709
 710   // Set differences between track and clusters in the bending and non bending directions
 711   AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr();
 712   Double_t dX1 = cluster1->GetX() - trackParamAtCluster1.GetNonBendingCoor();
 713   Double_t dX2 = cluster2->GetX() - trackParamAtCluster2.GetNonBendingCoor();
 714   Double_t dY1 = cluster1->GetY() - trackParamAtCluster1.GetBendingCoor();
 715   Double_t dY2 = cluster2->GetY() - trackParamAtCluster2.GetBendingCoor();
 716
 717   // Calculate errors and covariances
 718   const TMatrixD& kParamCov1 = trackParamAtCluster1.GetCovariances();
 719   const TMatrixD& kParamCov2 = trackParamAtCluster2.GetCovariances();
 720   Double_t dZ = trackParamAtCluster2.GetZ() - trackParamAtCluster1.GetZ();
 721   Double_t sigma2X1 = kParamCov1(0,0) + cluster1->GetErrX2();
 722   Double_t sigma2X2 = kParamCov2(0,0) + cluster2->GetErrX2();
 723   Double_t covX1X2  = kParamCov1(0,0) + dZ * kParamCov1(0,1);
 724   Double_t sigma2Y1 = kParamCov1(2,2) + cluster1->GetErrY2();
 725   Double_t sigma2Y2 = kParamCov2(2,2) + cluster2->GetErrY2();
 726   Double_t covY1Y2  = kParamCov1(2,2) + dZ * kParamCov1(2,3);
 727
 728   // Compute chi2
 729   Double_t detX = sigma2X1 * sigma2X2 - covX1X2 * covX1X2;
 730   Double_t detY = sigma2Y1 * sigma2Y2 - covY1Y2 * covY1Y2;
 731   if (detX == 0. || detY == 0.) return 2.*AliMUONTrack::MaxChi2();
 732   return   (dX1 * dX1 * sigma2X2 + dX2 * dX2 * sigma2X1 - 2. * dX1 * dX2 * covX1X2) / detX
 733          + (dY1 * dY1 * sigma2Y2 + dY2 * dY2 * sigma2Y1 - 2. * dY1 * dY2 * covY1Y2) / detY;
 734
 735 }
 736
 737   //__________________________________________________________________________
 738 Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInChamber(AliMUONTrack &trackCandidate, const AliMUONVClusterStore& clusterStore,
 739                                                               Int_t nextChamber)
 740 {
 741   /// Follow trackCandidate in chamber(0..) nextChamber assuming linear propagation, and search for compatible cluster(s)
 742   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
 743   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
 744   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
 745   ///         Remove the obsolete "trackCandidate" at the end.
 746   /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
 747   /// return kTRUE if new cluster(s) have been found (otherwise return kFALSE)
 748   AliDebug(1,Form("Enter FollowLinearTrackInChamber(1..) %d", nextChamber+1));
 749
 750   Double_t chi2WithOneCluster = AliMUONTrack::MaxChi2();
 751   Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
 752                                         GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
 753   Double_t bestChi2WithOneCluster = maxChi2WithOneCluster;
 754   Bool_t foundOneCluster = kFALSE;
 755   AliMUONTrack *newTrack = 0x0;
 756   AliMUONVCluster *cluster;
 757   AliMUONTrackParam trackParam;
 758   AliMUONTrackParam extrapTrackParamAtCluster;
 759   AliMUONTrackParam bestTrackParamAtCluster;
 760
 761   // Get track parameters according to the propagation direction
 762   if (nextChamber > 7) trackParam = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last();
 763   else trackParam = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First();
 764
 765   // Printout for debuging
 766   if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 767     cout<<endl<<"Track parameters and covariances at first cluster:"<<endl;
 768     trackParam.GetParameters().Print();
 769     trackParam.GetCovariances().Print();
 770   }
 771
 772   // Add MCS effect
 773   AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(trackParam.GetClusterPtr()->GetChamberId()),-1.);
 774
 775   // Printout for debuging
 776   if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 777     cout << "FollowLinearTrackInChamber: look for cluster in chamber(1..): " << nextChamber+1 << endl;
 778   }
 779
 780   // Create iterators to loop over clusters in chamber
 781   TIter next(clusterStore.CreateChamberIterator(nextChamber,nextChamber));
 782
 783   // look for candidates in chamber
 784   while ( ( cluster = static_cast<AliMUONVCluster*>(next()) ) ) {
 785
 786     // try to add the current cluster fast
 787     if (!TryOneClusterFast(trackParam, cluster)) continue;
 788
 789     // try to add the current cluster accuratly
 790     extrapTrackParamAtCluster = trackParam;
 791     AliMUONTrackExtrap::LinearExtrapToZCov(&extrapTrackParamAtCluster, cluster->GetZ());
 792     chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCluster, cluster, extrapTrackParamAtCluster);
 793
 794     // if good chi2 then consider to add cluster
 795     if (chi2WithOneCluster < maxChi2WithOneCluster) {
 796       foundOneCluster = kTRUE;
 797
 798       // Printout for debuging
 799       if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 800         cout << "FollowLinearTrackInChamber: found one cluster in chamber(1..): " << nextChamber+1
 801         << " (Chi2 = " << chi2WithOneCluster << ")" << endl;
 802         cluster->Print();
 803       }
 804
 805       if (GetRecoParam()->TrackAllTracks()) {
 806         // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
 807         newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 808         if (GetRecoParam()->RequestStation(nextChamber/2))
 809           extrapTrackParamAtCluster.SetRemovable(kFALSE);
 810         else extrapTrackParamAtCluster.SetRemovable(kTRUE);
 811         newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster,*cluster);
 812         newTrack->SetGlobalChi2(trackCandidate.GetGlobalChi2()+chi2WithOneCluster);
 813         fNRecTracks++;
 814
 815         // Printout for debuging
 816         if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 817           cout << "FollowLinearTrackInChamber: added one cluster in chamber(1..): " << nextChamber+1 << endl;
 818           if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 819         }
 820
 821       } else if (chi2WithOneCluster < bestChi2WithOneCluster) {
 822         // keep track of the best cluster
 823         bestChi2WithOneCluster = chi2WithOneCluster;
 824         bestTrackParamAtCluster = extrapTrackParamAtCluster;
 825       }
 826
 827     }
 828
 829   }
 830
 831   // fill out the best track if required else clean up the fRecTracksPtr array
 832   if (!GetRecoParam()->TrackAllTracks()) {
 833     if (foundOneCluster) {
 834       if (GetRecoParam()->RequestStation(nextChamber/2))
 835         bestTrackParamAtCluster.SetRemovable(kFALSE);
 836       else bestTrackParamAtCluster.SetRemovable(kTRUE);
 837       trackCandidate.AddTrackParamAtCluster(bestTrackParamAtCluster,*(bestTrackParamAtCluster.GetClusterPtr()));
 838       trackCandidate.SetGlobalChi2(trackCandidate.GetGlobalChi2()+bestChi2WithOneCluster);
 839
 840       // Printout for debuging
 841       if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 842         cout << "FollowLinearTrackInChamber: added the best cluster in chamber(1..): " << bestTrackParamAtCluster.GetClusterPtr()->GetChamberId()+1 << endl;
 843         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 844       }
 845
 846     } else return kFALSE;
 847
 848   } else if (foundOneCluster) {
 849
 850     // remove obsolete track
 851     fRecTracksPtr->Remove(&trackCandidate);
 852     fNRecTracks--;
 853
 854   } else return kFALSE;
 855
 856   return kTRUE;
 857
 858 }
 859
 860 //__________________________________________________________________________
 861 Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trackCandidate, const AliMUONVClusterStore& clusterStore,
 862                                                               Int_t nextStation)
 863 {
 864   /// Follow trackCandidate in station(0..) nextStation assuming linear propagation, and search for compatible cluster(s)
 865   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
 866   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
 867   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
 868   ///         Remove the obsolete "trackCandidate" at the end.
 869   /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
 870   /// return kTRUE if new cluster(s) have been found (otherwise return kFALSE)
 871   AliDebug(1,Form("Enter FollowLinearTrackInStation(1..) %d", nextStation+1));
 872
 873   // Order the chamber according to the propagation direction (tracking starts with chamber 2):
 874   // - nextStation == station(1...) 5 => forward propagation
 875   // - nextStation < station(1...) 5 => backward propagation
 876   Int_t ch1, ch2;
 877   if (nextStation==4) {
 878     ch1 = 2*nextStation+1;
 879     ch2 = 2*nextStation;
 880   } else {
 881     ch1 = 2*nextStation;
 882     ch2 = 2*nextStation+1;
 883   }
 884
 885   Double_t chi2WithOneCluster = AliMUONTrack::MaxChi2();
 886   Double_t chi2WithTwoClusters = AliMUONTrack::MaxChi2();
 887   Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
 888                                         GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
 889   Double_t maxChi2WithTwoClusters = 4. * GetRecoParam()->GetSigmaCutForTracking() *
 890                                          GetRecoParam()->GetSigmaCutForTracking(); // 4 because 4 quantities in chi2
 891   Double_t bestChi2WithOneCluster = maxChi2WithOneCluster;
 892   Double_t bestChi2WithTwoClusters = maxChi2WithTwoClusters;
 893   Bool_t foundOneCluster = kFALSE;
 894   Bool_t foundTwoClusters = kFALSE;
 895   AliMUONTrack *newTrack = 0x0;
 896   AliMUONVCluster *clusterCh1, *clusterCh2;
 897   AliMUONTrackParam trackParam;
 898   AliMUONTrackParam extrapTrackParamAtCluster1;
 899   AliMUONTrackParam extrapTrackParamAtCluster2;
 900   AliMUONTrackParam bestTrackParamAtCluster1;
 901   AliMUONTrackParam bestTrackParamAtCluster2;
 902
 903   Int_t nClusters = clusterStore.GetSize();
 904   Bool_t *clusterCh1Used = new Bool_t[nClusters];
 905   for (Int_t i = 0; i < nClusters; i++) clusterCh1Used[i] = kFALSE;
 906   Int_t iCluster1;
 907
 908   // Get track parameters according to the propagation direction
 909   if (nextStation==4) trackParam = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last();
 910   else trackParam = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First();
 911
 912   // Printout for debuging
 913   if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 914     cout<<endl<<"Track parameters and covariances at first cluster:"<<endl;
 915     trackParam.GetParameters().Print();
 916     trackParam.GetCovariances().Print();
 917   }
 918
 919   // Add MCS effect
 920   AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(trackParam.GetClusterPtr()->GetChamberId()),-1.);
 921
 922   // Printout for debuging
 923   if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 924     cout << "FollowLinearTrackInStation: look for clusters in chamber(1..): " << ch2+1 << endl;
 925   }
 926
 927   // Create iterators to loop over clusters in both chambers
 928   TIter nextInCh1(clusterStore.CreateChamberIterator(ch1,ch1));
 929   TIter nextInCh2(clusterStore.CreateChamberIterator(ch2,ch2));
 930
 931   // look for candidates in chamber 2
 932   while ( ( clusterCh2 = static_cast<AliMUONVCluster*>(nextInCh2()) ) ) {
 933
 934     // try to add the current cluster fast
 935     if (!TryOneClusterFast(trackParam, clusterCh2)) continue;
 936
 937     // try to add the current cluster accuratly
 938     extrapTrackParamAtCluster2 = trackParam;
 939     AliMUONTrackExtrap::LinearExtrapToZCov(&extrapTrackParamAtCluster2, clusterCh2->GetZ());
 940     chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCluster2, clusterCh2, extrapTrackParamAtCluster2);
 941
 942     // if good chi2 then try to attach a cluster in the other chamber too
 943     if (chi2WithOneCluster < maxChi2WithOneCluster) {
 944       Bool_t foundSecondCluster = kFALSE;
 945
 946       // Printout for debuging
 947       if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 948         cout << "FollowLinearTrackInStation: found one cluster in chamber(1..): " << ch2+1
 949              << " (Chi2 = " << chi2WithOneCluster << ")" << endl;
 950         clusterCh2->Print();
 951         cout << "                      look for second clusters in chamber(1..): " << ch1+1 << " ..." << endl;
 952       }
 953
 954       // add MCS effect
 955       AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCluster2,AliMUONConstants::ChamberThicknessInX0(ch2),-1.);
 956
 957       // reset cluster iterator of chamber 1
 958       nextInCh1.Reset();
 959       iCluster1 = -1;
 960
 961       // look for second candidates in chamber 1
 962       while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
 963         iCluster1++;
 964
 965         // try to add the current cluster fast
 966         if (!TryOneClusterFast(extrapTrackParamAtCluster2, clusterCh1)) continue;
 967
 968         // try to add the current cluster in addition to the one found in the previous chamber
 969         chi2WithTwoClusters = TryTwoClustersFast(extrapTrackParamAtCluster2, clusterCh1, extrapTrackParamAtCluster1);
 970
 971         // if good chi2 then consider to add the 2 clusters to the "trackCandidate"
 972         if (chi2WithTwoClusters < maxChi2WithTwoClusters) {
 973           foundSecondCluster = kTRUE;
 974           foundTwoClusters = kTRUE;
 975
 976           // Printout for debuging
 977           if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 978             cout << "FollowLinearTrackInStation: found one cluster in chamber(1..): " << ch1+1
 979                  << " (Chi2 = " << chi2WithTwoClusters << ")" << endl;
 980             clusterCh1->Print();
 981           }
 982
 983           if (GetRecoParam()->TrackAllTracks()) {
 984             // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new clusters
 985             newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 986             extrapTrackParamAtCluster1.SetRemovable(kTRUE);
 987             newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster1,*clusterCh1);
 988             extrapTrackParamAtCluster2.SetRemovable(kTRUE);
 989             newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster2,*clusterCh2);
 990             newTrack->SetGlobalChi2(newTrack->GetGlobalChi2()+chi2WithTwoClusters);
 991             fNRecTracks++;
 992
 993             // Tag clusterCh1 as used
 994             clusterCh1Used[iCluster1] = kTRUE;
 995
 996             // Printout for debuging
 997             if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 998               cout << "FollowLinearTrackInStation: added two clusters in station(1..): " << nextStation+1 << endl;
 999               if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
1000             }
1001
1002           } else if (chi2WithTwoClusters < bestChi2WithTwoClusters) {
1003             // keep track of the best couple of clusters
1004             bestChi2WithTwoClusters = chi2WithTwoClusters;
1005             bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
1006             bestTrackParamAtCluster2 = extrapTrackParamAtCluster2;
1007           }
1008
1009         }
1010
1011       }
1012
1013       // if no cluster found in chamber1 then consider to add clusterCh2 only
1014       if (!foundSecondCluster) {
1015         foundOneCluster = kTRUE;
1016
1017         if (GetRecoParam()->TrackAllTracks()) {
1018           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
1019           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
1020           if (GetRecoParam()->RequestStation(nextStation))
1021             extrapTrackParamAtCluster2.SetRemovable(kFALSE);
1022           else extrapTrackParamAtCluster2.SetRemovable(kTRUE);
1023           newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster2,*clusterCh2);
1024           newTrack->SetGlobalChi2(newTrack->GetGlobalChi2()+chi2WithOneCluster);
1025           fNRecTracks++;
1026
1027           // Printout for debuging
1028           if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1029             cout << "FollowLinearTrackInStation: added one cluster in chamber(1..): " << ch2+1 << endl;
1030             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
1031           }
1032
1033         } else if (!foundTwoClusters && chi2WithOneCluster < bestChi2WithOneCluster) {
1034           // keep track of the best cluster except if a couple of clusters has already been found
1035           bestChi2WithOneCluster = chi2WithOneCluster;
1036           bestTrackParamAtCluster1 = extrapTrackParamAtCluster2;
1037         }
1038
1039       }
1040
1041     }
1042
1043   }
1044
1045   // look for candidates in chamber 1 not already attached to a track
1046   // if we want to keep all possible tracks or if no good couple of clusters has been found
1047   if (GetRecoParam()->TrackAllTracks() || !foundTwoClusters) {
1048
1049     // Printout for debuging
1050     if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1051       cout << "FollowLinearTrackInStation: look for single cluster in chamber(1..): " << ch1+1 << endl;
1052     }
1053
1054     //Extrapolate trackCandidate to chamber "ch2"
1055     AliMUONTrackExtrap::LinearExtrapToZCov(&trackParam, AliMUONConstants::DefaultChamberZ(ch2));
1056
1057     // add MCS effect for next step
1058     AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(ch2),-1.);
1059
1060     // reset cluster iterator of chamber 1
1061     nextInCh1.Reset();
1062     iCluster1 = -1;
1063
1064     // look for second candidates in chamber 1
1065     while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
1066       iCluster1++;
1067
1068       if (clusterCh1Used[iCluster1]) continue; // Skip clusters already used
1069
1070       // try to add the current cluster fast
1071       if (!TryOneClusterFast(trackParam, clusterCh1)) continue;
1072
1073       // try to add the current cluster accuratly
1074       extrapTrackParamAtCluster1 = trackParam;
1075       AliMUONTrackExtrap::LinearExtrapToZCov(&extrapTrackParamAtCluster1, clusterCh1->GetZ());
1076       chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCluster1, clusterCh1, extrapTrackParamAtCluster1);
1077
1078       // if good chi2 then consider to add clusterCh1
1079       // We do not try to attach a cluster in the other chamber too since it has already been done above
1080       if (chi2WithOneCluster < maxChi2WithOneCluster) {
1081         foundOneCluster = kTRUE;
1082
1083         // Printout for debuging
1084         if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1085           cout << "FollowLinearTrackInStation: found one cluster in chamber(1..): " << ch1+1
1086                << " (Chi2 = " << chi2WithOneCluster << ")" << endl;
1087           clusterCh1->Print();
1088         }
1089
1090         if (GetRecoParam()->TrackAllTracks()) {
1091           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
1092           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
1093           if (GetRecoParam()->RequestStation(nextStation))
1094             extrapTrackParamAtCluster1.SetRemovable(kFALSE);
1095           else extrapTrackParamAtCluster1.SetRemovable(kTRUE);
1096           newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster1,*clusterCh1);
1097           newTrack->SetGlobalChi2(newTrack->GetGlobalChi2()+chi2WithOneCluster);
1098           fNRecTracks++;
1099
1100           // Printout for debuging
1101           if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1102             cout << "FollowLinearTrackInStation: added one cluster in chamber(1..): " << ch1+1 << endl;
1103             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
1104           }
1105
1106         } else if (chi2WithOneCluster < bestChi2WithOneCluster) {
1107           // keep track of the best cluster except if a couple of clusters has already been found
1108           bestChi2WithOneCluster = chi2WithOneCluster;
1109           bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
1110         }
1111
1112       }
1113
1114     }
1115
1116   }
1117
1118   // fill out the best track if required else clean up the fRecTracksPtr array
1119   if (!GetRecoParam()->TrackAllTracks()) {
1120     if (foundTwoClusters) {
1121       bestTrackParamAtCluster1.SetRemovable(kTRUE);
1122       trackCandidate.AddTrackParamAtCluster(bestTrackParamAtCluster1,*(bestTrackParamAtCluster1.GetClusterPtr()));
1123       bestTrackParamAtCluster2.SetRemovable(kTRUE);
1124       trackCandidate.AddTrackParamAtCluster(bestTrackParamAtCluster2,*(bestTrackParamAtCluster2.GetClusterPtr()));
1125       trackCandidate.SetGlobalChi2(trackCandidate.GetGlobalChi2()+bestChi2WithTwoClusters);
1126
1127       // Printout for debuging
1128       if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1129         cout << "FollowLinearTrackInStation: added the two best clusters in station(1..): " << nextStation+1 << endl;
1130         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
1131       }
1132
1133     } else if (foundOneCluster) {
1134       if (GetRecoParam()->RequestStation(nextStation))
1135         bestTrackParamAtCluster1.SetRemovable(kFALSE);
1136       else bestTrackParamAtCluster1.SetRemovable(kTRUE);
1137       trackCandidate.AddTrackParamAtCluster(bestTrackParamAtCluster1,*(bestTrackParamAtCluster1.GetClusterPtr()));
1138       trackCandidate.SetGlobalChi2(trackCandidate.GetGlobalChi2()+bestChi2WithOneCluster);
1139
1140       // Printout for debuging
1141       if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1142         cout << "FollowLinearTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster1.GetClusterPtr()->GetChamberId()+1 << endl;
1143         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
1144       }
1145
1146     } else {
1147       delete [] clusterCh1Used;
1148       return kFALSE;
1149     }
1150
1151   } else if (foundOneCluster || foundTwoClusters) {
1152
1153     // remove obsolete track
1154     fRecTracksPtr->Remove(&trackCandidate);
1155     fNRecTracks--;
1156
1157   } else {
1158     delete [] clusterCh1Used;
1159     return kFALSE;
1160   }
1161
1162   delete [] clusterCh1Used;
1163   return kTRUE;
1164
1165 }
1166
1167 //__________________________________________________________________________
1168 void AliMUONVTrackReconstructor::ImproveTracks()
1169 {
1170   /// Improve tracks by removing clusters with local chi2 highter than the defined cut
1171   /// Recompute track parameters and covariances at the remaining clusters
1172   AliDebug(1,"Enter ImproveTracks");
1173
1174   AliMUONTrack *track, *nextTrack;
1175
1176   track = (AliMUONTrack*) fRecTracksPtr->First();
1177   while (track) {
1178
1179     // prepare next track in case the actual track is suppressed
1180     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track);
1181
1182     ImproveTrack(*track);
1183
1184     // remove track if improvement failed
1185     if (!track->IsImproved()) {
1186       fRecTracksPtr->Remove(track);
1187       fNRecTracks--;
1188     }
1189
1190     track = nextTrack;
1191   }
1192
1193   // compress the array in case of some tracks have been removed
1194   fRecTracksPtr->Compress();
1195
1196 }
1197
1198 //__________________________________________________________________________
1199 void AliMUONVTrackReconstructor::Finalize()
1200 {
1201   /// Recompute track parameters and covariances at each attached cluster from those at the first one
1202   /// Set the label pointing to the corresponding MC track
1203   /// Remove the track if finalization failed
1204
1205   AliMUONTrack *track, *nextTrack;
1206   Bool_t trackRemoved = kFALSE;
1207
1208   track = (AliMUONTrack*) fRecTracksPtr->First();
1209   while (track) {
1210
1211     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track);
1212
1213     if (FinalizeTrack(*track)) track->FindMCLabel();
1214     else {
1215       fRecTracksPtr->Remove(track);
1216       fNRecTracks--;
1217       trackRemoved = kTRUE;
1218     }
1219
1220     track = nextTrack;
1221
1222   }
1223
1224   // compress array of tracks if needed
1225   if (trackRemoved) fRecTracksPtr->Compress();
1226
1227 }
1228
1229 //__________________________________________________________________________
1230 void AliMUONVTrackReconstructor::ValidateTracksWithTrigger(AliMUONVTrackStore& trackStore,
1231                                                            const AliMUONVTriggerTrackStore& triggerTrackStore,
1232                                                            const AliMUONVTriggerStore& triggerStore,
1233                                                            const AliMUONTrackHitPattern& trackHitPattern)
1234 {
1235   /// Try to match track from tracking system with trigger track
1236   AliCodeTimerAuto("",0);
1237
1238   trackHitPattern.ExecuteValidation(trackStore, triggerTrackStore, triggerStore);
1239 }
1240
1241   //__________________________________________________________________________
1242 void AliMUONVTrackReconstructor::EventReconstructTrigger(const AliMUONTriggerCircuit& circuit,
1243                                                          const AliMUONVTriggerStore& triggerStore,
1244                                                          AliMUONVTriggerTrackStore& triggerTrackStore)
1245 {
1246   /// To make the trigger tracks from Local Trigger
1247   AliDebug(1, "");
1248   AliCodeTimerAuto("",0);
1249
1250   AliMUONGlobalTrigger* globalTrigger = triggerStore.Global();
1251
1252   UChar_t gloTrigPat = 0;
1253
1254   if (globalTrigger)
1255   {
1256     gloTrigPat = globalTrigger->GetGlobalResponse();
1257   }
1258
1259   TIter next(triggerStore.CreateIterator());
1260   AliMUONLocalTrigger* locTrg(0x0);
1261
1262   const Double_t kTrigNonBendReso = AliMUONConstants::TriggerNonBendingReso();
1263   const Double_t kTrigBendReso = AliMUONConstants::TriggerBendingReso();
1264   const Double_t kSqrt12 = TMath::Sqrt(12.);
1265
1266   AliMUONTriggerTrack triggerTrack;
1267   TMatrixD trigCov(3,3);
1268
1269   while ( ( locTrg = static_cast<AliMUONLocalTrigger*>(next()) ) )
1270   {
1271     if ( locTrg->IsTrigX() && locTrg->IsTrigY() )
1272     { // make Trigger Track if trigger in X and Y
1273
1274       Int_t localBoardId = locTrg->LoCircuit();
1275
1276       Float_t y11 = circuit.GetY11Pos(localBoardId, locTrg->LoStripX());
1277       Float_t z11 = circuit.GetZ11Pos(localBoardId, locTrg->LoStripX());
1278       // need first to convert deviation to [0-30]
1279       // (see AliMUONLocalTriggerBoard::LocalTrigger)
1280       Int_t deviation = locTrg->GetDeviation();
1281       Int_t stripX21 = locTrg->LoStripX()+deviation+1;
1282       Float_t y21 = circuit.GetY21Pos(localBoardId, stripX21);
1283       Float_t z21 = circuit.GetZ21Pos(localBoardId, stripX21);
1284       Float_t x11 = circuit.GetX11Pos(localBoardId, locTrg->LoStripY());
1285
1286       AliDebug(1, Form(" MakeTriggerTrack %3d %2d %2d %2d (%f %f %f) (%f %f)\n",locTrg->LoCircuit(),
1287                        locTrg->LoStripX(),locTrg->LoStripX()+deviation+1,locTrg->LoStripY(),x11, y11, z11, y21, z21));
1288
1289
1290       Double_t deltaZ = z11 - z21;
1291
1292       Float_t slopeX = x11/z11;
1293       Float_t slopeY = (y11-y21) / deltaZ;
1294
1295       Float_t sigmaX = circuit.GetX11Width(localBoardId, locTrg->LoStripY()) / kSqrt12;
1296       Float_t sigmaY = circuit.GetY11Width(localBoardId, locTrg->LoStripX()) / kSqrt12;
1297       Float_t sigmaY21 = circuit.GetY21Width(localBoardId, locTrg->LoStripX()) / kSqrt12;
1298
1299       trigCov.Zero();
1300       trigCov(0,0) = kTrigNonBendReso * kTrigNonBendReso + sigmaX * sigmaX;
1301       trigCov(1,1) = kTrigBendReso * kTrigBendReso + sigmaY * sigmaY;
1302       trigCov(2,2) =
1303         (2. * kTrigBendReso * kTrigBendReso + sigmaY * sigmaY + sigmaY21 * sigmaY21 ) / deltaZ / deltaZ;
1304       trigCov(1,2) = trigCov(2,1) = trigCov(1,1) / deltaZ;
1305
1306       triggerTrack.SetX11(x11);
1307       triggerTrack.SetY11(y11);
1308       triggerTrack.SetZ11(z11);
1309       triggerTrack.SetZ21(z21);
1310       triggerTrack.SetSlopeX(slopeX);
1311       triggerTrack.SetSlopeY(slopeY);
1312       triggerTrack.SetGTPattern(gloTrigPat);
1313       triggerTrack.SetLoTrgNum(localBoardId);
1314       triggerTrack.SetCovariances(trigCov);
1315
1316       triggerTrackStore.Add(triggerTrack);
1317     } // board is fired
1318   } // end of loop on Local Trigger
1319 }