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