MUON/AliMUONTrack.cxx

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