MUON/AliMUONTrack.cxx

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