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