MFT/AliMuonForwardTrack.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 //====================================================================================================================================================
  17 //
  18 //      Description of an ALICE muon forward track, combining the information of the Muon Spectrometer and the Muon Forward Tracker
  19 //
  20 //      Contact author: antonio.uras@cern.ch
  21 //
  22 //====================================================================================================================================================
  23
  24 #include "AliLog.h"
  25 #include "AliMUONTrack.h"
  26 #include "AliMFTCluster.h"
  27 #include "AliMUONVCluster.h"
  28 #include "AliMUONTrackParam.h"
  29 #include "AliMUONTrackExtrap.h"
  30 #include "TClonesArray.h"
  31 #include "TMatrixD.h"
  32 #include "TParticle.h"
  33 #include "AliMuonForwardTrack.h"
  34
  35 ClassImp(AliMuonForwardTrack)
  36
  37 //====================================================================================================================================================
  38
  39 AliMuonForwardTrack::AliMuonForwardTrack():
  40   AliMUONTrack(),
  41   fMUONTrack(0),
  42   fMCTrackRef(0),
  43   fMFTClusters(0)
  44 {
  45
  46   // default constructor
  47   for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) fPlaneExists[iPlane] = kFALSE;
  48   fMFTClusters = new TClonesArray("AliMFTCluster");
  49
  50 }
  51
  52 //====================================================================================================================================================
  53
  54 AliMuonForwardTrack::AliMuonForwardTrack(AliMUONTrack *MUONTrack):
  55   AliMUONTrack(),
  56   fMUONTrack(0),
  57   fMCTrackRef(0),
  58   fMFTClusters(0)
  59 {
  60
  61   SetMUONTrack(MUONTrack);
  62   for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) fPlaneExists[iPlane] = kFALSE;
  63   fMFTClusters = new TClonesArray("AliMFTCluster");
  64
  65 }
  66
  67 //====================================================================================================================================================
  68
  69 AliMuonForwardTrack::AliMuonForwardTrack(const AliMuonForwardTrack& track):
  70   AliMUONTrack(track),
  71   fMUONTrack(track.fMUONTrack),
  72   fMCTrackRef(track.fMCTrackRef),
  73   fMFTClusters(track.fMFTClusters)
  74 {
  75
  76   // copy constructor
  77   for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
  78
  79 }
  80
  81 //====================================================================================================================================================
  82
  83 AliMuonForwardTrack& AliMuonForwardTrack::operator=(const AliMuonForwardTrack& track) {
  84
  85   // Asignment operator
  86
  87   // check assignement to self
  88   if (this == &track) return *this;
  89
  90   // base class assignement
  91   AliMUONTrack::operator=(track);
  92
  93   // clear memory
  94   Clear();
  95
  96   fMUONTrack    = track.fMUONTrack;
  97   fMCTrackRef   = track.fMCTrackRef;
  98   fMFTClusters  = track.fMFTClusters;
  99
 100   for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
 101
 102   return *this;
 103
 104 }
 105
 106 //====================================================================================================================================================
 107
 108 void AliMuonForwardTrack::SetMUONTrack(AliMUONTrack *MUONTrack) {
 109
 110   if (fMUONTrack) {
 111     AliInfo("fMUONTrack already exists, nothing will be done");
 112     return;
 113   }
 114
 115   fMUONTrack = MUONTrack;
 116
 117 }
 118
 119 //====================================================================================================================================================
 120
 121 void AliMuonForwardTrack::SetMCTrackRef(TParticle *MCTrackRef) {
 122
 123   if (fMCTrackRef) {
 124     AliInfo("fMCTrackRef already exists, nothing will be done");
 125     return;
 126   }
 127
 128   fMCTrackRef = MCTrackRef;
 129
 130 }
 131
 132 //====================================================================================================================================================
 133
 134 void AliMuonForwardTrack::AddTrackParamAtMFTCluster(AliMUONTrackParam &trackParam, AliMFTCluster &mftCluster) {
 135
 136   AliDebug(1, Form("fMFTClusters=%p has %d entries", fMFTClusters, fMFTClusters->GetEntries()));
 137   Int_t iMFTCluster = fMFTClusters->GetEntries();
 138   AliDebug(1, Form("mftCluster->GetX() = %f  mftCluster->GetY() = %f  mftCluster->GetErrX() = %f  cmftCluster->GetErrY() = %f",
 139            mftCluster.GetX(), mftCluster.GetY(), mftCluster.GetErrX(), mftCluster.GetErrY()));
 140   AliMUONVCluster *muonCluster = (AliMUONVCluster*) mftCluster.CreateMUONCluster();
 141   AliDebug(1, Form("Created MUON cluster %p", muonCluster));
 142   trackParam.SetUniqueID(iMFTCluster);    // we profit of this slot to store the reference to the corresponding MFTCluster
 143   AliDebug(1, Form("Now adding muonCluster %p and trackParam %p",muonCluster, &trackParam));
 144   AddTrackParamAtCluster(trackParam, *muonCluster, kTRUE);
 145   AliDebug(1, Form("GetTrackParamAtCluster() = %p  has %d entries",GetTrackParamAtCluster(), GetTrackParamAtCluster()->GetEntries()));
 146   // we pass the parameters this->GetTrackParamAtCluster()->First() to the Kalman Filter algorithm: they will be updated!!
 147   Double_t chi2Kalman = RunKalmanFilter(*(AliMUONTrackParam*)(GetTrackParamAtCluster()->First()));
 148   Double_t newGlobalChi2 = GetGlobalChi2() + chi2Kalman;
 149   mftCluster.SetLocalChi2(chi2Kalman);
 150   mftCluster.SetTrackChi2(newGlobalChi2);
 151   new ((*fMFTClusters)[iMFTCluster]) AliMFTCluster(mftCluster);
 152   AliDebug(1, Form("muonCluster->GetZ() = %f, trackParam->GetZ() = %f",muonCluster->GetZ(), trackParam.GetZ()));
 153   SetGlobalChi2(newGlobalChi2);
 154   ((AliMUONTrackParam*) GetTrackParamAtCluster()->First())->SetTrackChi2(newGlobalChi2);
 155   for (Int_t iPar=0; iPar<GetTrackParamAtCluster()->GetEntries(); iPar++) {
 156     AliDebug(1, Form("GetTrackParamAtCluster()->At(%d)->GetClusterPtr() = %p",
 157                      iPar, ((AliMUONTrackParam*)(GetTrackParamAtCluster()->At(iPar)))->GetClusterPtr()));
 158   }
 159
 160 }
 161
 162 //====================================================================================================================================================
 163
 164 AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMUONCluster(Int_t iMUONCluster) {
 165
 166   if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
 167     AliError("Invalid MUON cluster index. NULL pointer will be returned");
 168     return NULL;
 169   }
 170
 171   AliMUONTrackParam *trackParam = (AliMUONTrackParam*) fMUONTrack->GetTrackParamAtCluster()->At(iMUONCluster);
 172
 173   return trackParam;
 174
 175 }
 176
 177 //====================================================================================================================================================
 178
 179 AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMFTCluster(Int_t iMFTCluster) {
 180
 181   if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
 182     AliError("Invalid MFT cluster index. NULL pointer will be returned");
 183     return NULL;
 184   }
 185
 186   AliMUONTrackParam *trackParam = (AliMUONTrackParam*) GetTrackParamAtCluster()->At(iMFTCluster);
 187
 188   return trackParam;
 189
 190 }
 191
 192 //====================================================================================================================================================
 193
 194 AliMUONVCluster* AliMuonForwardTrack::GetMUONCluster(Int_t iMUONCluster) {
 195
 196   if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
 197     AliError("Invalid MUON cluster index. NULL pointer will be returned");
 198     return NULL;
 199   }
 200
 201   AliMUONTrackParam *trackParam = GetTrackParamAtMUONCluster(iMUONCluster);
 202   AliMUONVCluster *muonCluster = trackParam->GetClusterPtr();
 203
 204   return muonCluster;
 205
 206 }
 207
 208 //====================================================================================================================================================
 209
 210 AliMFTCluster* AliMuonForwardTrack::GetMFTCluster(Int_t iMFTCluster) {
 211
 212   if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
 213     AliError("Invalid MFT cluster index. NULL pointer will be returned");
 214     return NULL;
 215   }
 216
 217   AliMUONTrackParam *trackParam = GetTrackParamAtMFTCluster(iMFTCluster);
 218   AliMFTCluster *mftCluster = (AliMFTCluster*) fMFTClusters->At(trackParam->GetUniqueID());
 219
 220   return mftCluster;
 221
 222 }
 223
 224 //====================================================================================================================================================
 225
 226 Double_t AliMuonForwardTrack::RunKalmanFilter(AliMUONTrackParam &trackParamAtCluster) {
 227
 228   AliDebug(1, Form("Running Kalman filter for parameters %p (z = %f) and cluster %p (z = %f)",
 229                    &trackParamAtCluster, trackParamAtCluster.GetZ(), trackParamAtCluster.GetClusterPtr(), trackParamAtCluster.GetClusterPtr()->GetZ()));
 230
 231   // Compute new track parameters and their covariances including new cluster using kalman filter
 232   // return the *additional* track chi2
 233
 234   // Get actual track parameters (p)
 235   TMatrixD param(trackParamAtCluster.GetParameters());
 236
 237   // Get new cluster parameters (m)
 238   AliMUONVCluster *cluster = trackParamAtCluster.GetClusterPtr();
 239   AliDebug(1, Form("cluster->GetX() = %f  cluster->GetY() = %f  cluster->GetErrX() = %f  cluster->GetErrY() = %f",
 240                    cluster->GetX(), cluster->GetY(), cluster->GetErrX(), cluster->GetErrY()));
 241   TMatrixD clusterParam(5,1);
 242   clusterParam.Zero();
 243   clusterParam(0,0) = cluster->GetX();
 244   clusterParam(2,0) = cluster->GetY();
 245
 246   // Compute the current parameter weight (W)
 247   TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
 248   if (paramWeight.Determinant() != 0) {
 249     paramWeight.Invert();
 250   } else {
 251     Warning("RunKalmanFilter"," Determinant = 0");
 252     return 1.e10;
 253   }
 254
 255   // Compute the new cluster weight (U)
 256   TMatrixD clusterWeight(5,5);
 257   clusterWeight.Zero();
 258   clusterWeight(0,0) = 1. / cluster->GetErrX2();
 259   clusterWeight(2,2) = 1. / cluster->GetErrY2();
 260
 261   // Compute the new parameters covariance matrix ( (W+U)^-1 )
 262   TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
 263   if (newParamCov.Determinant() != 0) {
 264     newParamCov.Invert();
 265   }
 266   else {
 267     Warning("RunKalmanFilter"," Determinant = 0");
 268     return 1.e10;
 269   }
 270
 271   // Save the new parameters covariance matrix
 272   trackParamAtCluster.SetCovariances(newParamCov);
 273
 274   // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
 275   TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
 276   TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp);    // U(m-p)
 277   TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
 278   newParam += param;                                   // ((W+U)^-1)U(m-p) + p
 279
 280   // Save the new parameters
 281   trackParamAtCluster.SetParameters(newParam);
 282
 283   // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
 284   tmp = newParam; // p'
 285   tmp -= param;   // (p'-p)
 286   TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp);           // W(p'-p)
 287   TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
 288   tmp = newParam;      // p'
 289   tmp -= clusterParam; // (p'-m)
 290   TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp);            // U(p'-m)
 291   addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
 292
 293   AliDebug(1,Form("Adding Kalman chi2 = %f",addChi2Track(0,0)));
 294
 295   return addChi2Track(0,0);
 296
 297 }
 298
 299 //====================================================================================================================================================
 300
 301 Double_t AliMuonForwardTrack::GetWeightedOffset(Double_t x, Double_t y, Double_t z) {
 302
 303   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
 304   AliMUONTrackExtrap::ExtrapToZCov(param, z);
 305
 306   TMatrixD cov(5,5);
 307   cov = param->GetCovariances();
 308
 309   TMatrixD covCoordinates(2,2);
 310   covCoordinates(0,0) = cov(0,0);
 311   covCoordinates(0,1) = cov(0,2);
 312   covCoordinates(1,0) = cov(2,0);
 313   covCoordinates(1,1) = cov(2,2);
 314
 315   TMatrixD covCoordinatesInverse = covCoordinates.Invert();
 316
 317   Double_t dX = param->GetNonBendingCoor() - x;
 318   Double_t dY = param->GetBendingCoor()    - y;
 319
 320   Double_t weightedOffset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) +
 321                                              dY*dY*covCoordinatesInverse(1,1) +
 322                                              2.*dX*dY*covCoordinatesInverse(0,1)));
 323
 324   return weightedOffset;
 325
 326 }
 327
 328 //====================================================================================================================================================
 329
 330 Double_t AliMuonForwardTrack::GetOffsetX(Double_t x, Double_t z) {
 331
 332   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
 333   AliMUONTrackExtrap::ExtrapToZCov(param, z);
 334   Double_t dX = param->GetNonBendingCoor() - x;
 335   return dX;
 336
 337 }
 338
 339 //====================================================================================================================================================
 340
 341 Double_t AliMuonForwardTrack::GetOffsetY(Double_t y, Double_t z) {
 342
 343   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
 344   AliMUONTrackExtrap::ExtrapToZCov(param, z);
 345   Double_t dY = param->GetBendingCoor() - y;
 346   return dY;
 347
 348 }
 349
 350 //====================================================================================================================================================
 351
 352 Double_t AliMuonForwardTrack::GetOffset(Double_t x, Double_t y, Double_t z) {
 353
 354   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
 355   AliMUONTrackExtrap::ExtrapToZCov(param, z);
 356   Double_t dX = param->GetNonBendingCoor() - x;
 357   Double_t dY = param->GetBendingCoor()    - y;
 358   return TMath::Sqrt(dX*dX + dY*dY);
 359
 360 }
 361
 362 //====================================================================================================================================================
 363