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