MUON/AliMUONTrack.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 AliMUONTrack
  20 //-------------------
  21 // Reconstructed track in ALICE dimuon spectrometer
  22 //-----------------------------------------------------------------------------
  23
  24 #include "AliMUONTrack.h"
  25
  26 #include "AliMUONReconstructor.h"
  27 #include "AliMUONVCluster.h"
  28 #include "AliMUONVClusterStore.h"
  29 #include "AliMUONObjectPair.h"
  30 #include "AliMUONTrackExtrap.h"
  31 #include "AliMUONConstants.h"
  32 #include "AliMUONTrackParam.h"
  33
  34 #include "AliLog.h"
  35
  36 #include <TMath.h>
  37
  38 #include <Riostream.h>
  39
  40 /// \cond CLASSIMP
  41 ClassImp(AliMUONTrack) // Class implementation in ROOT context
  42 /// \endcond
  43
  44 //__________________________________________________________________________
  45 AliMUONTrack::AliMUONTrack()
  46   : TObject(),
  47     fTrackParamAtCluster(0x0),
  48     fFitWithVertex(kFALSE),
  49     fVertexErrXY2(),
  50     fFitWithMCS(kFALSE),
  51     fClusterWeightsNonBending(0x0),
  52     fClusterWeightsBending(0x0),
  53     fGlobalChi2(-1.),
  54     fImproved(kFALSE),
  55     fMatchTrigger(-1),
  56     floTrgNum(-1),
  57     fChi2MatchTrigger(0.),
  58     fTrackID(-1),
  59     fTrackParamAtVertex(0x0),
  60     fHitsPatternInTrigCh(0),
  61     fLocalTrigger(0)
  62 {
  63   /// Default constructor
  64   fVertexErrXY2[0] = 0.;
  65   fVertexErrXY2[1] = 0.;
  66 }
  67
  68   //__________________________________________________________________________
  69 AliMUONTrack::AliMUONTrack(AliMUONObjectPair *segment, Double_t bendingVertexDispersion)
  70   : TObject(),
  71     fTrackParamAtCluster(new TClonesArray("AliMUONTrackParam",10)),
  72     fFitWithVertex(kFALSE),
  73     fVertexErrXY2(),
  74     fFitWithMCS(kFALSE),
  75     fClusterWeightsNonBending(0x0),
  76     fClusterWeightsBending(0x0),
  77     fGlobalChi2(0.),
  78     fImproved(kFALSE),
  79     fMatchTrigger(-1),
  80     floTrgNum(-1),
  81     fChi2MatchTrigger(0.),
  82     fTrackID(-1),
  83     fTrackParamAtVertex(0x0),
  84     fHitsPatternInTrigCh(0),
  85     fLocalTrigger(0)
  86 {
  87   /// Constructor from two clusters
  88
  89   fVertexErrXY2[0] = 0.;
  90   fVertexErrXY2[1] = 0.;
  91
  92   // Pointers to clusters from the segment
  93   AliMUONVCluster* firstCluster = (AliMUONVCluster*) segment->First();
  94   AliMUONVCluster* lastCluster = (AliMUONVCluster*) segment->Second();
  95
  96   // Compute track parameters
  97   Double_t z1 = firstCluster->GetZ();
  98   Double_t z2 = lastCluster->GetZ();
  99   Double_t dZ = z1 - z2;
 100   // Non bending plane
 101   Double_t nonBendingCoor1 = firstCluster->GetX();
 102   Double_t nonBendingCoor2 = lastCluster->GetX();
 103   Double_t nonBendingSlope = (nonBendingCoor1 - nonBendingCoor2) / dZ;
 104   // Bending plane
 105   Double_t bendingCoor1 = firstCluster->GetY();
 106   Double_t bendingCoor2 = lastCluster->GetY();
 107   Double_t bendingSlope = (bendingCoor1 - bendingCoor2) / dZ;
 108   // Inverse bending momentum
 109   Double_t bendingImpact = bendingCoor1 - z1 * bendingSlope;
 110   Double_t inverseBendingMomentum = 1. / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact);
 111
 112   // Set track parameters at first cluster
 113   AliMUONTrackParam trackParamAtFirstCluster;
 114   trackParamAtFirstCluster.SetZ(z1);
 115   trackParamAtFirstCluster.SetNonBendingCoor(nonBendingCoor1);
 116   trackParamAtFirstCluster.SetNonBendingSlope(nonBendingSlope);
 117   trackParamAtFirstCluster.SetBendingCoor(bendingCoor1);
 118   trackParamAtFirstCluster.SetBendingSlope(bendingSlope);
 119   trackParamAtFirstCluster.SetInverseBendingMomentum(inverseBendingMomentum);
 120
 121   // Set track parameters at last cluster
 122   AliMUONTrackParam trackParamAtLastCluster;
 123   trackParamAtLastCluster.SetZ(z2);
 124   trackParamAtLastCluster.SetNonBendingCoor(nonBendingCoor2);
 125   trackParamAtLastCluster.SetNonBendingSlope(nonBendingSlope);
 126   trackParamAtLastCluster.SetBendingCoor(bendingCoor2);
 127   trackParamAtLastCluster.SetBendingSlope(bendingSlope);
 128   trackParamAtLastCluster.SetInverseBendingMomentum(inverseBendingMomentum);
 129
 130   // Compute and set track parameters covariances at first cluster
 131   TMatrixD paramCov(5,5);
 132   paramCov.Zero();
 133   // Non bending plane
 134   paramCov(0,0) = firstCluster->GetErrX2();
 135   paramCov(0,1) = firstCluster->GetErrX2() / dZ;
 136   paramCov(1,0) = paramCov(0,1);
 137   paramCov(1,1) = ( firstCluster->GetErrX2() + lastCluster->GetErrX2() ) / dZ / dZ;
 138   // Bending plane
 139   paramCov(2,2) = firstCluster->GetErrY2();
 140   paramCov(2,3) = firstCluster->GetErrY2() / dZ;
 141   paramCov(3,2) = paramCov(2,3);
 142   paramCov(3,3) = ( firstCluster->GetErrY2() + lastCluster->GetErrY2() ) / dZ / dZ;
 143   // Inverse bending momentum (vertex resolution + bending slope resolution + 10% error on dipole parameters+field)
 144   if (AliMUONTrackExtrap::IsFieldON()) {
 145     paramCov(4,4) = ( ( bendingVertexDispersion*bendingVertexDispersion +
 146                        (z1 * z1 * lastCluster->GetErrY2() + z2 * z2 * firstCluster->GetErrY2()) / dZ / dZ) /
 147                      bendingImpact / bendingImpact + 0.1 * 0.1) * inverseBendingMomentum * inverseBendingMomentum ;
 148     paramCov(2,4) = - z2 * firstCluster->GetErrY2() * inverseBendingMomentum / bendingImpact / dZ;
 149     paramCov(4,2) = paramCov(2,4);
 150     paramCov(3,4) = - (z1 * lastCluster->GetErrY2() + z2 * firstCluster->GetErrY2()) * inverseBendingMomentum / bendingImpact / dZ / dZ;
 151     paramCov(4,3) = paramCov(3,4);
 152   } else paramCov(4,4) = inverseBendingMomentum*inverseBendingMomentum;
 153   trackParamAtFirstCluster.SetCovariances(paramCov);
 154
 155   // Compute and set track parameters covariances at last cluster
 156   // Non bending plane
 157   paramCov(0,0) = lastCluster->GetErrX2();
 158   paramCov(0,1) = - lastCluster->GetErrX2() / dZ;
 159   paramCov(1,0) = paramCov(0,1);
 160   // Bending plane
 161   paramCov(2,2) = lastCluster->GetErrY2();
 162   paramCov(2,3) = - lastCluster->GetErrY2() / dZ;
 163   paramCov(3,2) = paramCov(2,3);
 164   // Inverse bending momentum (vertex resolution + bending slope resolution + 10% error on dipole parameters+field)
 165   if (AliMUONTrackExtrap::IsFieldON()) {
 166     paramCov(2,4) = z1 * lastCluster->GetErrY2() * inverseBendingMomentum / bendingImpact / dZ;
 167     paramCov(4,2) = paramCov(2,4);
 168   }
 169   trackParamAtLastCluster.SetCovariances(paramCov);
 170
 171   // Add track parameters at clusters
 172   AddTrackParamAtCluster(trackParamAtFirstCluster,*firstCluster);
 173   AddTrackParamAtCluster(trackParamAtLastCluster,*lastCluster);
 174
 175 }
 176
 177 //__________________________________________________________________________
 178 AliMUONTrack::AliMUONTrack(const AliMUONTrack& track)
 179   : TObject(track),
 180     fTrackParamAtCluster(0x0),
 181     fFitWithVertex(track.fFitWithVertex),
 182     fVertexErrXY2(),
 183     fFitWithMCS(track.fFitWithMCS),
 184     fClusterWeightsNonBending(0x0),
 185     fClusterWeightsBending(0x0),
 186     fGlobalChi2(track.fGlobalChi2),
 187     fImproved(track.fImproved),
 188     fMatchTrigger(track.fMatchTrigger),
 189     floTrgNum(track.floTrgNum),
 190     fChi2MatchTrigger(track.fChi2MatchTrigger),
 191     fTrackID(track.fTrackID),
 192     fTrackParamAtVertex(0x0),
 193     fHitsPatternInTrigCh(track.fHitsPatternInTrigCh),
 194     fLocalTrigger(track.fLocalTrigger)
 195 {
 196   ///copy constructor
 197
 198   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
 199   if (track.fTrackParamAtCluster) {
 200     fTrackParamAtCluster = new TClonesArray("AliMUONTrackParam",10);
 201     AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->First();
 202     while (trackParamAtCluster) {
 203       new ((*fTrackParamAtCluster)[GetNClusters()]) AliMUONTrackParam(*trackParamAtCluster);
 204       trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->After(trackParamAtCluster);
 205     }
 206   }
 207
 208   // copy vertex resolution square used during the tracking procedure
 209   fVertexErrXY2[0] = track.fVertexErrXY2[0];
 210   fVertexErrXY2[1] = track.fVertexErrXY2[1];
 211
 212   // copy cluster weights matrices if any
 213   if (track.fClusterWeightsNonBending) fClusterWeightsNonBending = new TMatrixD(*(track.fClusterWeightsNonBending));
 214   if (track.fClusterWeightsBending) fClusterWeightsBending = new TMatrixD(*(track.fClusterWeightsBending));
 215
 216   // copy track parameters at vertex if any
 217   if (track.fTrackParamAtVertex) fTrackParamAtVertex = new AliMUONTrackParam(*(track.fTrackParamAtVertex));
 218
 219 }
 220
 221   //__________________________________________________________________________
 222 AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& track)
 223 {
 224   /// Asignment operator
 225   // check assignement to self
 226   if (this == &track)
 227     return *this;
 228
 229   // base class assignement
 230   TObject::operator=(track);
 231
 232   // clear memory
 233   Clear();
 234
 235   // necessary to make a copy of the objects and not only the pointers in TClonesArray
 236   if (track.fTrackParamAtCluster) {
 237     fTrackParamAtCluster = new TClonesArray("AliMUONTrackParam",10);
 238     AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->First();
 239     while (trackParamAtCluster) {
 240       new ((*fTrackParamAtCluster)[GetNClusters()]) AliMUONTrackParam(*trackParamAtCluster);
 241       trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->After(trackParamAtCluster);
 242     }
 243   }
 244
 245   // copy cluster weights matrix if any
 246   if (track.fClusterWeightsNonBending) {
 247     if (fClusterWeightsNonBending) {
 248       fClusterWeightsNonBending->ResizeTo(*(track.fClusterWeightsNonBending));
 249       *fClusterWeightsNonBending = *(track.fClusterWeightsNonBending);
 250     } else fClusterWeightsNonBending = new TMatrixD(*(track.fClusterWeightsNonBending));
 251   }
 252
 253   // copy cluster weights matrix if any
 254   if (track.fClusterWeightsBending) {
 255     if (fClusterWeightsBending) {
 256       fClusterWeightsBending->ResizeTo(*(track.fClusterWeightsBending));
 257       *fClusterWeightsBending = *(track.fClusterWeightsBending);
 258     } else fClusterWeightsBending = new TMatrixD(*(track.fClusterWeightsBending));
 259   }
 260
 261   // copy track parameters at vertex if any
 262   if (track.fTrackParamAtVertex) {
 263     if (fTrackParamAtVertex) *fTrackParamAtVertex = *(track.fTrackParamAtVertex);
 264     else fTrackParamAtVertex = new AliMUONTrackParam(*(track.fTrackParamAtVertex));
 265   }
 266
 267   fFitWithVertex      =  track.fFitWithVertex;
 268   fVertexErrXY2[0]    =  track.fVertexErrXY2[0];
 269   fVertexErrXY2[1]    =  track.fVertexErrXY2[1];
 270   fFitWithMCS         =  track.fFitWithMCS;
 271   fGlobalChi2         =  track.fGlobalChi2;
 272   fImproved           =  track.fImproved;
 273   fMatchTrigger       =  track.fMatchTrigger;
 274   floTrgNum           =  track.floTrgNum;
 275   fChi2MatchTrigger   =  track.fChi2MatchTrigger;
 276   fTrackID            =  track.fTrackID;
 277   fHitsPatternInTrigCh = track.fHitsPatternInTrigCh;
 278   fLocalTrigger        = track.fLocalTrigger;
 279
 280   return *this;
 281 }
 282
 283   //__________________________________________________________________________
 284 AliMUONTrack::~AliMUONTrack()
 285 {
 286   /// Destructor
 287   delete fTrackParamAtCluster;
 288   delete fClusterWeightsNonBending;
 289   delete fClusterWeightsBending;
 290   delete fTrackParamAtVertex;
 291 }
 292
 293   //__________________________________________________________________________
 294 void AliMUONTrack::Clear(Option_t* opt)
 295 {
 296   /// Clear arrays
 297   if (opt && opt[0] == 'C' && fTrackParamAtCluster) fTrackParamAtCluster->Clear("C");
 298   else {
 299     delete fTrackParamAtCluster;
 300     fTrackParamAtCluster = 0x0;
 301   }
 302   delete fClusterWeightsNonBending; fClusterWeightsNonBending = 0x0;
 303   delete fClusterWeightsBending; fClusterWeightsBending = 0x0;
 304   delete fTrackParamAtVertex; fTrackParamAtVertex = 0x0;
 305 }
 306
 307   //__________________________________________________________________________
 308 void AliMUONTrack::Reset()
 309 {
 310   /// Reset to default values
 311   SetUniqueID(0);
 312   fFitWithVertex = kFALSE;
 313   fVertexErrXY2[0] = 0.;
 314   fVertexErrXY2[1] = 0.;
 315   fFitWithMCS = kFALSE;
 316   fGlobalChi2 = -1.;
 317   fImproved = kFALSE;
 318   fMatchTrigger = -1;
 319   floTrgNum = -1;
 320   fChi2MatchTrigger = 0.;
 321   fTrackID = -1;
 322   fHitsPatternInTrigCh = 0;
 323   fLocalTrigger = 0;
 324   delete fTrackParamAtCluster; fTrackParamAtCluster = 0x0;
 325   delete fClusterWeightsNonBending; fClusterWeightsNonBending = 0x0;
 326   delete fClusterWeightsBending; fClusterWeightsBending = 0x0;
 327   delete fTrackParamAtVertex; fTrackParamAtVertex = 0x0;
 328 }
 329
 330   //__________________________________________________________________________
 331 TClonesArray* AliMUONTrack::GetTrackParamAtCluster() const
 332 {
 333   /// return array of track parameters at cluster (create it if needed)
 334   if (!fTrackParamAtCluster) fTrackParamAtCluster = new TClonesArray("AliMUONTrackParam",10);
 335   return fTrackParamAtCluster;
 336 }
 337
 338   //__________________________________________________________________________
 339 void AliMUONTrack::AddTrackParamAtCluster(const AliMUONTrackParam &trackParam, AliMUONVCluster &cluster, Bool_t copy)
 340 {
 341   /// Copy given track parameters into a new TrackParamAtCluster
 342   /// Link parameters with the associated cluster
 343   /// If copy=kTRUE: the cluster is copied then passed the trackParam which become its owner
 344   ///     otherwise: make sure to do not delete the cluster until it is used by the track
 345
 346   // check chamber ID of the associated cluster
 347   if (cluster.GetChamberId() < 0 || cluster.GetChamberId() > AliMUONConstants::NTrackingCh()) {
 348     AliError(Form("Chamber ID of the associated cluster is not valid (ChamberId=%d)",cluster.GetChamberId()));
 349     return;
 350   }
 351
 352   // check whether track parameters are given at the correct cluster z position
 353   if (cluster.GetZ() != trackParam.GetZ()) {
 354     AliError("track parameters are given at a different z position than the one of the associated cluster");
 355     return;
 356   }
 357
 358   // add parameters to the array of track parameters
 359   if (!fTrackParamAtCluster) fTrackParamAtCluster = new TClonesArray("AliMUONTrackParam",10);
 360   AliMUONTrackParam* trackParamAtCluster = new ((*fTrackParamAtCluster)[GetNClusters()]) AliMUONTrackParam(trackParam);
 361
 362   // link parameters with the associated cluster or its copy
 363   if (copy) {
 364     AliMUONVCluster *clusterCopy = static_cast<AliMUONVCluster*>(cluster.Clone());
 365     trackParamAtCluster->SetClusterPtr(clusterCopy, kTRUE);
 366   } else trackParamAtCluster->SetClusterPtr(&cluster);
 367
 368   // sort the array of track parameters
 369   fTrackParamAtCluster->Sort();
 370 }
 371
 372   //__________________________________________________________________________
 373 void AliMUONTrack::RemoveTrackParamAtCluster(AliMUONTrackParam *trackParam)
 374 {
 375   /// Remove trackParam from the array of TrackParamAtCluster
 376   if (!fTrackParamAtCluster || !fTrackParamAtCluster->Remove(trackParam)) {
 377     AliWarning("object to remove does not exist in array fTrackParamAtCluster");
 378     return;
 379   }
 380
 381   fTrackParamAtCluster->Compress();
 382 }
 383
 384   //__________________________________________________________________________
 385 void AliMUONTrack::UpdateTrackParamAtCluster()
 386 {
 387   /// Update track parameters at each attached cluster
 388
 389   if (GetNClusters() == 0) {
 390     AliWarning("no cluster attached to the track");
 391     return;
 392   }
 393
 394   AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->First();
 395   AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(startingTrackParam);
 396   while (trackParamAtCluster) {
 397
 398     // reset track parameters and their covariances
 399     trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
 400     trackParamAtCluster->SetZ(startingTrackParam->GetZ());
 401
 402     // extrapolation to the given z
 403     AliMUONTrackExtrap::ExtrapToZ(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ());
 404
 405     // prepare next step
 406     startingTrackParam = trackParamAtCluster;
 407     trackParamAtCluster = (AliMUONTrackParam*) (fTrackParamAtCluster->After(trackParamAtCluster));
 408   }
 409
 410 }
 411
 412   //__________________________________________________________________________
 413 void AliMUONTrack::UpdateCovTrackParamAtCluster()
 414 {
 415   /// Update track parameters and their covariances at each attached cluster
 416   /// Include effects of multiple scattering in chambers
 417
 418   if (GetNClusters() == 0) {
 419     AliWarning("no cluster attached to the track");
 420     return;
 421   }
 422
 423   AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->First();
 424   AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(startingTrackParam);
 425   Int_t expectedChamber = startingTrackParam->GetClusterPtr()->GetChamberId() + 1;
 426   Int_t currentChamber;
 427   while (trackParamAtCluster) {
 428
 429     // reset track parameters and their covariances
 430     trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
 431     trackParamAtCluster->SetZ(startingTrackParam->GetZ());
 432     trackParamAtCluster->SetCovariances(startingTrackParam->GetCovariances());
 433
 434     // add MCS effect
 435     AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
 436
 437     // add MCS in missing chambers if any
 438     currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
 439     while (currentChamber > expectedChamber) {
 440       // extrapolation to the missing chamber
 441       AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, AliMUONConstants::DefaultChamberZ(expectedChamber));
 442       // add MCS effect
 443       AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
 444       expectedChamber++;
 445     }
 446
 447     // extrapolation to the z of the current cluster
 448     AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ());
 449
 450     // prepare next step
 451     expectedChamber = currentChamber + 1;
 452     startingTrackParam = trackParamAtCluster;
 453     trackParamAtCluster = (AliMUONTrackParam*) (fTrackParamAtCluster->After(trackParamAtCluster));
 454   }
 455
 456 }
 457
 458   //__________________________________________________________________________
 459 Bool_t AliMUONTrack::IsValid(UInt_t requestedStationMask)
 460 {
 461   /// check the validity of the current track:
 462   /// at least one cluster per requested station
 463   /// and at least 2 chambers in stations 4 & 5 that contain cluster(s)
 464
 465   Int_t nClusters = GetNClusters();
 466   AliMUONTrackParam *trackParam;
 467   Int_t currentCh, currentSt, previousCh = -1, nChHitInSt45 = 0;
 468   UInt_t presentStationMask(0);
 469
 470   // first loop over clusters
 471   for (Int_t i = 0; i < nClusters; i++) {
 472     trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i);
 473
 474     currentCh = trackParam->GetClusterPtr()->GetChamberId();
 475     currentSt = currentCh/2;
 476
 477     // build present station mask
 478     presentStationMask |= ( 1 << currentSt );
 479
 480     // count the number of chambers in station 4 & 5 that contain cluster(s)
 481     if (currentCh > 5 && currentCh != previousCh) {
 482       nChHitInSt45++;
 483       previousCh = currentCh;
 484     }
 485
 486   }
 487
 488   return (((requestedStationMask & presentStationMask) == requestedStationMask) && (nChHitInSt45 >= 2));
 489 }
 490
 491   //__________________________________________________________________________
 492 void AliMUONTrack::TagRemovableClusters(UInt_t requestedStationMask) {
 493   /// Identify clusters that can be removed from the track,
 494   /// with the only requirements to have at least 1 cluster per requested station
 495   /// and at least 2 chambers over 4 in stations 4 & 5 that contain cluster(s)
 496
 497   Int_t nClusters = GetNClusters();
 498   AliMUONTrackParam *trackParam, *nextTrackParam;
 499   Int_t currentCh, nextCh, currentSt, nextSt, previousCh = -1, nChHitInSt45 = 0;
 500
 501   // first loop over clusters
 502   for (Int_t i = 0; i < nClusters; i++) {
 503     trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i);
 504
 505     currentCh = trackParam->GetClusterPtr()->GetChamberId();
 506     currentSt = currentCh/2;
 507
 508     // reset flags to kFALSE for all clusters in required station
 509     if ((1 << currentSt) & requestedStationMask) trackParam->SetRemovable(kFALSE);
 510     else trackParam->SetRemovable(kTRUE);
 511
 512     // count the number of chambers in station 4 & 5 that contain cluster(s)
 513     if (currentCh > 5 && currentCh != previousCh) {
 514       nChHitInSt45++;
 515       previousCh = currentCh;
 516     }
 517
 518   }
 519
 520   // second loop over clusters
 521   for (Int_t i = 0; i < nClusters; i++) {
 522     trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i);
 523
 524     currentCh = trackParam->GetClusterPtr()->GetChamberId();
 525     currentSt = currentCh/2;
 526
 527     // make sure they are more than 2 clusters in 2 different chambers of stations 4 & 5
 528     // but 2 clusters in he same chamber will still be flagged as removable
 529     if (nChHitInSt45 < 3 && currentSt > 2) {
 530
 531       if (i == nClusters-1) {
 532
 533         trackParam->SetRemovable(kFALSE);
 534
 535       } else {
 536
 537         nextTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i+1);
 538         nextCh = nextTrackParam->GetClusterPtr()->GetChamberId();
 539
 540         // set clusters in the same chamber as being removable
 541         if (nextCh == currentCh) {
 542           trackParam->SetRemovable(kTRUE);
 543           nextTrackParam->SetRemovable(kTRUE);
 544           i++; // skip cluster already checked
 545         } else {
 546           trackParam->SetRemovable(kFALSE);
 547         }
 548
 549       }
 550
 551     } else {
 552
 553       // skip clusters already flag as removable
 554       if (trackParam->IsRemovable()) continue;
 555
 556       // loop over next track parameters
 557       for (Int_t j = i+1; j < nClusters; j++) {
 558         nextTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(j);
 559
 560         nextCh = nextTrackParam->GetClusterPtr()->GetChamberId();
 561         nextSt = nextCh/2;
 562
 563         // set clusters in the same station as being removable
 564         if (nextSt == currentSt) {
 565           trackParam->SetRemovable(kTRUE);
 566           nextTrackParam->SetRemovable(kTRUE);
 567           i++; // skip cluster already checked
 568         }
 569
 570       }
 571
 572     }
 573
 574   }
 575
 576 }
 577
 578   //__________________________________________________________________________
 579 Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS)
 580 {
 581   /// Compute each cluster contribution to the chi2 of the track
 582   /// accounting for multiple scattering or not according to the flag
 583   /// - Also recompute the weight matrices of the attached clusters if accountForMCS=kTRUE
 584   /// - Assume that track parameters at each cluster are corrects
 585   /// - Return kFALSE if computation failed
 586   AliDebug(1,"Enter ComputeLocalChi2");
 587
 588   if (!fTrackParamAtCluster) {
 589     AliWarning("no cluster attached to this track");
 590     return kFALSE;
 591   }
 592
 593   if (accountForMCS) { // Compute local chi2 taking into account multiple scattering effects
 594
 595     // Compute MCS covariance matrix only once
 596     Int_t nClusters = GetNClusters();
 597     TMatrixD mcsCovariances(nClusters,nClusters);
 598     ComputeMCSCovariances(mcsCovariances);
 599
 600     // Make sure cluster weights are consistent with following calculations
 601     if (!ComputeClusterWeights(&mcsCovariances)) {
 602       AliWarning("cannot take into account the multiple scattering effects");
 603       return ComputeLocalChi2(kFALSE);
 604     }
 605
 606     // Compute chi2 of the track
 607     Double_t globalChi2 = ComputeGlobalChi2(kTRUE);
 608     if (globalChi2 < 0.) return kFALSE;
 609
 610     // Loop over removable clusters and compute their local chi2
 611     AliMUONTrackParam* trackParamAtCluster1;
 612     AliMUONVCluster *cluster, *discardedCluster;
 613     Int_t iCluster1, iCluster2, iCurrentCluster1, iCurrentCluster2;
 614     TMatrixD clusterWeightsNB(nClusters-1,nClusters-1);
 615     TMatrixD clusterWeightsB(nClusters-1,nClusters-1);
 616     Double_t *dX = new Double_t[nClusters-1];
 617     Double_t *dY = new Double_t[nClusters-1];
 618     Double_t globalChi2b;
 619     AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->First();
 620     while (trackParamAtCluster) {
 621
 622       discardedCluster = trackParamAtCluster->GetClusterPtr();
 623
 624       // Recompute cluster weights without the current cluster
 625       if (!ComputeClusterWeights(clusterWeightsNB, clusterWeightsB, &mcsCovariances, discardedCluster)) {
 626         AliWarning("cannot take into account the multiple scattering effects");
 627         delete [] dX;
 628         delete [] dY;
 629         return ComputeLocalChi2(kFALSE);
 630       }
 631
 632       // Compute track chi2 without the current cluster
 633       globalChi2b = 0.;
 634       iCurrentCluster1 = 0;
 635       for (iCluster1 = 0; iCluster1 < nClusters ; iCluster1++) {
 636         trackParamAtCluster1 = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1);
 637         cluster = trackParamAtCluster1->GetClusterPtr();
 638
 639         if (cluster == discardedCluster) continue;
 640
 641         // Compute and save residuals
 642         dX[iCurrentCluster1] = cluster->GetX() - trackParamAtCluster1->GetNonBendingCoor();
 643         dY[iCurrentCluster1] = cluster->GetY() - trackParamAtCluster1->GetBendingCoor();
 644
 645         iCurrentCluster2 = 0;
 646         for (iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) {
 647           cluster = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
 648
 649           if (cluster == discardedCluster) continue;
 650
 651           // Add contribution from covariances
 652           globalChi2b += (clusterWeightsNB(iCurrentCluster1, iCurrentCluster2) +
 653                           clusterWeightsNB(iCurrentCluster2, iCurrentCluster1)) * dX[iCurrentCluster1] * dX[iCurrentCluster2] +
 654                          (clusterWeightsB(iCurrentCluster1, iCurrentCluster2) +
 655                           clusterWeightsB(iCurrentCluster2, iCurrentCluster1)) * dY[iCurrentCluster1] * dY[iCurrentCluster2];
 656
 657           iCurrentCluster2++;
 658         }
 659
 660         // Add contribution from variances
 661         globalChi2b += clusterWeightsNB(iCurrentCluster1, iCurrentCluster1) * dX[iCurrentCluster1] * dX[iCurrentCluster1] +
 662                        clusterWeightsB(iCurrentCluster1, iCurrentCluster1) * dY[iCurrentCluster1] * dY[iCurrentCluster1];
 663
 664         iCurrentCluster1++;
 665       }
 666
 667       // Set local chi2
 668       trackParamAtCluster->SetLocalChi2(globalChi2 - globalChi2b);
 669
 670       trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(trackParamAtCluster);
 671     }
 672
 673     delete [] dX;
 674     delete [] dY;
 675
 676   } else { // without multiple scattering effects
 677
 678     AliMUONVCluster *discardedCluster;
 679     Double_t dX, dY;
 680     AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->First();
 681     while (trackParamAtCluster) {
 682
 683       discardedCluster = trackParamAtCluster->GetClusterPtr();
 684
 685       // Compute residuals
 686       dX = discardedCluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
 687       dY = discardedCluster->GetY() - trackParamAtCluster->GetBendingCoor();
 688
 689       // Set local chi2
 690       trackParamAtCluster->SetLocalChi2(dX * dX / discardedCluster->GetErrX2() + dY * dY / discardedCluster->GetErrY2());
 691
 692       trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(trackParamAtCluster);
 693     }
 694
 695   }
 696
 697   return kTRUE;
 698
 699 }
 700
 701   //__________________________________________________________________________
 702 Double_t AliMUONTrack::ComputeGlobalChi2(Bool_t accountForMCS)
 703 {
 704   /// Compute the chi2 of the track accounting for multiple scattering or not according to the flag
 705   /// - Assume that track parameters at each cluster are corrects
 706   /// - Assume the cluster weights matrices are corrects
 707   /// - Return negative value if chi2 computation failed
 708   AliDebug(1,"Enter ComputeGlobalChi2");
 709
 710   if (!fTrackParamAtCluster) {
 711     AliWarning("no cluster attached to this track");
 712     return 1.e10;
 713   }
 714
 715   Double_t chi2 = 0.;
 716
 717   if (accountForMCS) {
 718
 719     // Check the weight matrices. If weight matrices are not available compute chi2 without MCS
 720     if (!fClusterWeightsNonBending || !fClusterWeightsBending) {
 721       AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
 722       return ComputeGlobalChi2(kFALSE);
 723     }
 724     Int_t nClusters = GetNClusters();
 725     if (fClusterWeightsNonBending->GetNrows() != nClusters || fClusterWeightsBending->GetNcols() != nClusters) {
 726       AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
 727       return ComputeGlobalChi2(kFALSE);
 728     }
 729
 730     // Compute chi2
 731     AliMUONVCluster *cluster;
 732     Double_t *dX = new Double_t[nClusters];
 733     Double_t *dY = new Double_t[nClusters];
 734     AliMUONTrackParam* trackParamAtCluster;
 735     for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) {
 736       trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1);
 737       cluster = trackParamAtCluster->GetClusterPtr();
 738       dX[iCluster1] = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
 739       dY[iCluster1] = cluster->GetY() - trackParamAtCluster->GetBendingCoor();
 740       for (Int_t iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) {
 741         chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster2) + (*fClusterWeightsNonBending)(iCluster2, iCluster1)) * dX[iCluster1] * dX[iCluster2] +
 742                 ((*fClusterWeightsBending)(iCluster1, iCluster2) + (*fClusterWeightsBending)(iCluster2, iCluster1)) * dY[iCluster1] * dY[iCluster2];
 743       }
 744       chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster1) * dX[iCluster1] * dX[iCluster1]) +
 745               ((*fClusterWeightsBending)(iCluster1, iCluster1) * dY[iCluster1] * dY[iCluster1]);
 746     }
 747     delete [] dX;
 748     delete [] dY;
 749
 750   } else {
 751
 752     AliMUONVCluster *cluster;
 753     Double_t dX, dY;
 754     AliMUONTrackParam* trackParamAtCluster;
 755     Int_t nClusters = GetNClusters();
 756     for (Int_t iCluster = 0; iCluster < nClusters ; iCluster++) {
 757       trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster);
 758       cluster = trackParamAtCluster->GetClusterPtr();
 759       dX = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
 760       dY = cluster->GetY() - trackParamAtCluster->GetBendingCoor();
 761       chi2 += dX * dX / cluster->GetErrX2() + dY * dY / cluster->GetErrY2();
 762     }
 763
 764   }
 765
 766   return chi2;
 767
 768 }
 769
 770   //__________________________________________________________________________
 771 Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD* mcsCovariances)
 772 {
 773   /// Compute the weight matrices of the attached clusters, in non bending and bending direction,
 774   /// accounting for multiple scattering correlations and cluster resolution
 775   /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
 776   /// - Assume that track parameters at each cluster are corrects
 777   /// - Return kFALSE if computation failed
 778   AliDebug(1,"Enter ComputeClusterWeights1");
 779
 780   if (!fTrackParamAtCluster) {
 781     AliWarning("no cluster attached to this track");
 782     return kFALSE;
 783   }
 784
 785   // Alocate memory
 786   Int_t nClusters = GetNClusters();
 787   if (!fClusterWeightsNonBending) fClusterWeightsNonBending = new TMatrixD(nClusters,nClusters);
 788   if (!fClusterWeightsBending) fClusterWeightsBending = new TMatrixD(nClusters,nClusters);
 789
 790   // Compute weights matrices
 791   if (!ComputeClusterWeights(*fClusterWeightsNonBending, *fClusterWeightsBending, mcsCovariances)) return kFALSE;
 792
 793   return kTRUE;
 794
 795 }
 796
 797   //__________________________________________________________________________
 798 Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD& clusterWeightsNB, TMatrixD& clusterWeightsB,
 799                                            TMatrixD* mcsCovariances, AliMUONVCluster* discardedCluster) const
 800 {
 801   /// Compute the weight matrices, in non bending and bending direction,
 802   /// of the other attached clusters assuming the discarded one does not exist
 803   /// accounting for multiple scattering correlations and cluster resolution
 804   /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
 805   /// - Return kFALSE if computation failed
 806   AliDebug(1,"Enter ComputeClusterWeights2");
 807
 808   // Check MCS covariance matrix and recompute it if need
 809   Int_t nClusters = GetNClusters();
 810   Bool_t deleteMCSCov = kFALSE;
 811   if (!mcsCovariances) {
 812     mcsCovariances = new TMatrixD(nClusters,nClusters);
 813     deleteMCSCov = kTRUE;
 814     ComputeMCSCovariances(*mcsCovariances);
 815   }
 816
 817   // Resize the weights matrices; alocate memory
 818   if (discardedCluster) {
 819     clusterWeightsNB.ResizeTo(nClusters-1,nClusters-1);
 820     clusterWeightsB.ResizeTo(nClusters-1,nClusters-1);
 821   } else {
 822     clusterWeightsNB.ResizeTo(nClusters,nClusters);
 823     clusterWeightsB.ResizeTo(nClusters,nClusters);
 824   }
 825
 826   // Define variables
 827   AliMUONVCluster *cluster1, *cluster2;
 828   Int_t iCurrentCluster1, iCurrentCluster2;
 829
 830   // Compute the covariance matrices
 831   iCurrentCluster1 = 0;
 832   for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) {
 833     cluster1 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1))->GetClusterPtr();
 834
 835     if (cluster1 == discardedCluster) continue;
 836
 837     // Loop over next clusters
 838     iCurrentCluster2 = iCurrentCluster1;
 839     for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) {
 840       cluster2 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
 841
 842       if (cluster2 == discardedCluster) continue;
 843
 844       // Fill with MCS covariances
 845       clusterWeightsNB(iCurrentCluster1, iCurrentCluster2) = (*mcsCovariances)(iCluster1,iCluster2);
 846
 847       // Equal contribution from multiple scattering in non bending and bending directions
 848       clusterWeightsB(iCurrentCluster1, iCurrentCluster2) = clusterWeightsNB(iCurrentCluster1, iCurrentCluster2);
 849
 850       // Add contribution from cluster resolution to diagonal element and symmetrize the matrix
 851       if (iCurrentCluster1 == iCurrentCluster2) {
 852
 853         // In non bending plane
 854         clusterWeightsNB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrX2();
 855         // In bending plane
 856         clusterWeightsB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrY2();
 857
 858       } else {
 859
 860         // In non bending plane
 861         clusterWeightsNB(iCurrentCluster2, iCurrentCluster1) = clusterWeightsNB(iCurrentCluster1, iCurrentCluster2);
 862         // In bending plane
 863         clusterWeightsB(iCurrentCluster2, iCurrentCluster1) = clusterWeightsB(iCurrentCluster1, iCurrentCluster2);
 864
 865       }
 866
 867       iCurrentCluster2++;
 868     }
 869
 870     iCurrentCluster1++;
 871   }
 872
 873   // Inversion of covariance matrices to get the weights
 874   if (clusterWeightsNB.Determinant() != 0 && clusterWeightsB.Determinant() != 0) {
 875     clusterWeightsNB.Invert();
 876     clusterWeightsB.Invert();
 877   } else {
 878     AliWarning(" Determinant = 0");
 879     clusterWeightsNB.ResizeTo(0,0);
 880     clusterWeightsB.ResizeTo(0,0);
 881     if(deleteMCSCov) delete mcsCovariances;
 882     return kFALSE;
 883   }
 884
 885   if(deleteMCSCov) delete mcsCovariances;
 886
 887   return kTRUE;
 888
 889 }
 890
 891   //__________________________________________________________________________
 892 void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const
 893 {
 894   /// Compute the multiple scattering covariance matrix
 895   /// (assume that track parameters at each cluster are corrects)
 896   AliDebug(1,"Enter ComputeMCSCovariances");
 897
 898   // Reset the size of the covariance matrix if needed
 899   Int_t nClusters = GetNClusters();
 900   if (mcsCovariances.GetNrows() != nClusters) mcsCovariances.ResizeTo(nClusters,nClusters);
 901
 902   // Define variables
 903   Int_t nChambers = AliMUONConstants::NTrackingCh();
 904   AliMUONTrackParam* trackParamAtCluster;
 905   AliMUONTrackParam extrapTrackParam;
 906   Int_t currentChamber = 0, expectedChamber = 0, size = 0;
 907   Double_t *mcsAngle2 = new Double_t[2*nChambers];
 908   Double_t *zMCS = new Double_t[2*nChambers];
 909   Int_t *indices = new Int_t[2*nClusters];
 910
 911   // Compute multiple scattering dispersion angle at each chamber
 912   // and save the z position where it is calculated
 913   for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
 914     trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster);
 915
 916     // look for missing chambers if any
 917     currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
 918     while (currentChamber > expectedChamber) {
 919
 920       // Save the z position where MCS dispersion is calculated
 921       zMCS[size] = AliMUONConstants::DefaultChamberZ(expectedChamber);
 922
 923       // Do not take into account MCS in chambers prior the first cluster
 924       if (iCluster > 0) {
 925
 926         // Get track parameters at missing chamber z
 927         extrapTrackParam = *trackParamAtCluster;
 928         AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, zMCS[size]);
 929
 930         // Save multiple scattering dispersion angle in missing chamber
 931         mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
 932
 933       } else mcsAngle2[size] = 0.;
 934
 935       expectedChamber++;
 936       size++;
 937     }
 938
 939     // Save z position where MCS dispersion is calculated
 940     zMCS[size] = trackParamAtCluster->GetZ();
 941
 942     // Save multiple scattering dispersion angle in current chamber
 943     mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(*trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
 944
 945     // Save indice in zMCS array corresponding to the current cluster
 946     indices[iCluster] = size;
 947
 948     expectedChamber = currentChamber + 1;
 949     size++;
 950   }
 951
 952   // complete array of z if last cluster is on the last but one chamber
 953   if (currentChamber != nChambers-1) zMCS[size++] = AliMUONConstants::DefaultChamberZ(nChambers-1);
 954
 955   // Compute the covariance matrix
 956   for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) {
 957
 958     for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) {
 959
 960       // Initialization to 0 (diagonal plus upper triangular part)
 961       mcsCovariances(iCluster1,iCluster2) = 0.;
 962
 963       // Compute contribution from multiple scattering in upstream chambers
 964       for (Int_t k = 0; k < indices[iCluster1]; k++) {
 965         mcsCovariances(iCluster1,iCluster2) += (zMCS[indices[iCluster1]] - zMCS[k]) * (zMCS[indices[iCluster2]] - zMCS[k]) * mcsAngle2[k];
 966       }
 967
 968       // Symetrize the matrix
 969       mcsCovariances(iCluster2,iCluster1) = mcsCovariances(iCluster1,iCluster2);
 970     }
 971
 972   }
 973
 974   delete [] mcsAngle2;
 975   delete [] zMCS;
 976   delete [] indices;
 977
 978 }
 979
 980   //__________________________________________________________________________
 981 Int_t AliMUONTrack::ClustersInCommon(AliMUONTrack* track) const
 982 {
 983   /// Returns the number of clusters in common between the current track ("this")
 984   /// and the track pointed to by "track".
 985   if (!fTrackParamAtCluster || !this->fTrackParamAtCluster) return 0;
 986   Int_t nCluster1 = this->GetNClusters();
 987   Int_t nCluster2 = track->GetNClusters();
 988   Int_t clustersInCommon = 0;
 989   AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2;
 990   // Loop over clusters of first track
 991   for(Int_t iCluster1 = 0; iCluster1 < nCluster1; iCluster1++) {
 992     trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->UncheckedAt(iCluster1);
 993     // Loop over clusters of second track
 994     for(Int_t iCluster2 = 0; iCluster2 < nCluster2; iCluster2++) {
 995       trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->UncheckedAt(iCluster2);
 996       // Increment "clustersInCommon" if both trackParamAtCluster1 & 2 point to the same cluster
 997       if ((trackParamAtCluster1->GetClusterPtr()) == (trackParamAtCluster2->GetClusterPtr())) {
 998         clustersInCommon++;
 999         break;
1000       }
1001     }
1002   }
1003   return clustersInCommon;
1004 }
1005
1006   //__________________________________________________________________________
1007 Int_t AliMUONTrack::ClustersInCommonInSt345(AliMUONTrack* track) const
1008 {
1009   /// Returns the number of clusters in common on stations 3, 4 and 5
1010   /// between the current track ("this") and the track pointed to by "track".
1011   if (!fTrackParamAtCluster || !this->fTrackParamAtCluster) return 0;
1012   Int_t nCluster1 = this->GetNClusters();
1013   Int_t nCluster2 = track->GetNClusters();
1014   Int_t clustersInCommon = 0;
1015   AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2;
1016   // Loop over clusters of first track
1017   for(Int_t iCluster1 = 0; iCluster1 < nCluster1; iCluster1++) {
1018     trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->UncheckedAt(iCluster1);
1019     if (trackParamAtCluster1->GetClusterPtr()->GetChamberId() < 4) continue;
1020     // Loop over clusters of second track
1021     for(Int_t iCluster2 = 0; iCluster2 < nCluster2; iCluster2++) {
1022       trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->UncheckedAt(iCluster2);
1023       if (trackParamAtCluster2->GetClusterPtr()->GetChamberId() < 4) continue;
1024       // Increment "clustersInCommon" if both trackParamAtCluster1 & 2 point to the same cluster
1025       if ((trackParamAtCluster1->GetClusterPtr()) == (trackParamAtCluster2->GetClusterPtr())) {
1026         clustersInCommon++;
1027         break;
1028       }
1029     }
1030   }
1031   return clustersInCommon;
1032 }
1033
1034   //__________________________________________________________________________
1035 Int_t AliMUONTrack::GetNDF() const
1036 {
1037   /// return the number of degrees of freedom
1038
1039   Int_t ndf = 2 * GetNClusters() - 5;
1040   return (ndf > 0) ? ndf : 0;
1041 }
1042
1043   //__________________________________________________________________________
1044 Double_t AliMUONTrack::GetNormalizedChi2() const
1045 {
1046   /// return the chi2 value divided by the number of degrees of freedom (or FLT_MAX if ndf <= 0)
1047
1048   Double_t ndf = (Double_t) GetNDF();
1049   return (ndf > 0.) ? fGlobalChi2 / ndf : FLT_MAX;
1050 }
1051
1052   //__________________________________________________________________________
1053 Int_t AliMUONTrack::CompatibleTrack(AliMUONTrack* track, Double_t sigmaCut, Bool_t compatibleCluster[10]) const
1054 {
1055   /// for each chamber: return kTRUE (kFALSE) if clusters are compatible (not compatible).
1056   /// nMatchClusters = number of clusters of "this" track matched with one cluster of track "track"
1057   AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2;
1058   AliMUONVCluster *cluster1, *cluster2;
1059   Double_t chi2, dX, dY;
1060   Double_t chi2Max = sigmaCut * sigmaCut;
1061
1062   // initialization
1063   Int_t nMatchClusters = 0;
1064   for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) compatibleCluster[ch] = kFALSE;
1065
1066   if (!fTrackParamAtCluster || !this->fTrackParamAtCluster) return nMatchClusters;
1067
1068   // Loop over clusters of first track
1069   trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->First();
1070   while (trackParamAtCluster1) {
1071
1072     cluster1 = trackParamAtCluster1->GetClusterPtr();
1073
1074     // Loop over clusters of second track
1075     trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->First();
1076     while (trackParamAtCluster2) {
1077
1078       cluster2 = trackParamAtCluster2->GetClusterPtr();
1079
1080       //prepare next step
1081       trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->After(trackParamAtCluster2);
1082
1083       // check DE Id
1084       if (cluster1->GetDetElemId() != cluster2->GetDetElemId()) continue;
1085
1086       // check local chi2
1087       dX = cluster1->GetX() - cluster2->GetX();
1088       dY = cluster1->GetY() - cluster2->GetY();
1089       chi2 = dX * dX / (cluster1->GetErrX2() + cluster2->GetErrX2()) + dY * dY / (cluster1->GetErrY2() + cluster2->GetErrY2());
1090       if (chi2 > 2. * chi2Max) continue; // 2 because 2 quantities in chi2
1091
1092       compatibleCluster[cluster1->GetChamberId()] = kTRUE;
1093       nMatchClusters++;
1094       break;
1095     }
1096
1097     trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->After(trackParamAtCluster1);
1098   }
1099
1100   return nMatchClusters;
1101 }
1102
1103 //__________________________________________________________________________
1104 void AliMUONTrack::SetTrackParamAtVertex(const AliMUONTrackParam* trackParam)
1105 {
1106   /// set track parameters at vertex
1107   if (trackParam == 0x0) return;
1108   if (fTrackParamAtVertex) *fTrackParamAtVertex = *trackParam;
1109   else fTrackParamAtVertex = new AliMUONTrackParam(*trackParam);
1110 }
1111
1112 //__________________________________________________________________________
1113 void AliMUONTrack::RecursiveDump() const
1114 {
1115   /// Recursive dump of AliMUONTrack, i.e. with dump of trackParamAtCluster and attached clusters
1116   AliMUONTrackParam *trackParamAtCluster;
1117   AliMUONVCluster *cluster;
1118   cout << "Recursive dump of Track: " << this << endl;
1119   // Track
1120   this->Dump();
1121   for (Int_t iCluster = 0; iCluster < GetNClusters(); iCluster++) {
1122     trackParamAtCluster = (AliMUONTrackParam*) ((*fTrackParamAtCluster)[iCluster]);
1123     // trackParamAtCluster
1124     cout << "trackParamAtCluster: " << trackParamAtCluster << " (index: " << iCluster << ")" << endl;
1125     trackParamAtCluster->Dump();
1126     cluster = trackParamAtCluster->GetClusterPtr();
1127     // cluster
1128     cout << "cluster: " << cluster << endl;
1129     cluster->Print();
1130   }
1131   return;
1132 }
1133
1134 //_____________________________________________-
1135 void AliMUONTrack::Print(Option_t*) const
1136 {
1137   /// Printing Track information
1138
1139   cout << "<AliMUONTrack> No.Clusters=" << setw(2)   << GetNClusters() <<
1140       ", Match2Trig=" << setw(1) << GetMatchTrigger()  <<
1141       ", LoTrgNum=" << setw(3) << GetLoTrgNum()  <<
1142     ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) <<  GetChi2MatchTrigger();
1143   cout << Form(" HitTriggerPattern %x",fHitsPatternInTrigCh);
1144   cout << Form(" MClabel=%d",fTrackID) << endl;
1145   if (fTrackParamAtCluster) fTrackParamAtCluster->First()->Print("FULL");
1146 }
1147
1148 //__________________________________________________________________________
1149 void AliMUONTrack::SetLocalTrigger(Int_t loCirc, Int_t loStripX, Int_t loStripY, Int_t loDev, Int_t loLpt, Int_t loHpt)
1150 {
1151   /// pack the local trigger information and store
1152
1153   if (loCirc < 0) return;
1154
1155   fLocalTrigger = 0;
1156   fLocalTrigger += loCirc;
1157   fLocalTrigger += loStripX << 8;
1158   fLocalTrigger += loStripY << 13;
1159   fLocalTrigger += loDev    << 17;
1160   fLocalTrigger += loLpt    << 22;
1161   fLocalTrigger += loHpt    << 24;
1162
1163 }
1164
1165 //__________________________________________________________________________
1166 void AliMUONTrack::FindMCLabel()
1167 {
1168   /// Determine the MC label from the label of the attached clusters and fill fMCLabel data member:
1169   /// More than 50% of clusters, including 1 before and 1 after the dipole, must share the same label
1170
1171   Int_t nClusters = GetNClusters();
1172   Int_t halfCluster = nClusters/2;
1173
1174   // reset MC label
1175   fTrackID = -1;
1176
1177   // loop over first clusters (if nClusters left < (nClusters-halfCluster) the conditions cannot be fulfilled)
1178   for (Int_t iCluster1 = 0; iCluster1 < nClusters-halfCluster; iCluster1++) {
1179     AliMUONVCluster* cluster1 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1))->GetClusterPtr();
1180
1181     // if the first cluster is not on station 1 or 2 the conditions cannot be fulfilled
1182     if (cluster1->GetChamberId() > 3) return;
1183
1184     Int_t label1 = cluster1->GetMCLabel();
1185     if (label1 < 0) continue;
1186
1187     Int_t nIdenticalLabel = 1;
1188
1189     // Loop over next clusters
1190     for (Int_t iCluster2 = iCluster1+1; iCluster2 < nClusters; iCluster2++) {
1191       AliMUONVCluster* cluster2 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
1192
1193       if (cluster2->GetMCLabel() != label1) continue;
1194
1195       nIdenticalLabel++;
1196
1197       // stop as soon as conditions are fulfilled
1198       if (nIdenticalLabel > halfCluster && cluster2->GetChamberId() > 5) {
1199         fTrackID = label1;
1200         return;
1201       }
1202
1203     }
1204
1205   }
1206
1207 }
1208