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