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