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