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   fMFTClusters -> SetOwner(kTRUE);
  58
  59 }
  60
  61 //====================================================================================================================================================
  62
  63 AliMuonForwardTrack::AliMuonForwardTrack(AliMUONTrack *MUONTrack):
  64   AliMUONTrack(),
  65   fMUONTrack(0),
  66   fMCTrackRef(0),
  67   fMFTClusters(0),
  68   fNWrongClustersMC(-1),
  69   fTrackMCId(-1)
  70 {
  71
  72   SetMUONTrack(MUONTrack);
  73   for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = kFALSE;
  74   for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
  75     fParentMCLabel[iParent] = -1;
  76     fParentPDGCode[iParent] =  0;
  77   }
  78   fMFTClusters = new TClonesArray("AliMFTCluster");
  79   fMFTClusters -> SetOwner(kTRUE);
  80
  81 }
  82
  83 //====================================================================================================================================================
  84
  85 AliMuonForwardTrack::AliMuonForwardTrack(const AliMuonForwardTrack& track):
  86   AliMUONTrack(track),
  87   fMUONTrack(0x0),
  88   fMCTrackRef(0x0),
  89   fMFTClusters(0x0),
  90   fNWrongClustersMC(track.fNWrongClustersMC),
  91   fTrackMCId(track.fTrackMCId)
  92 {
  93
  94   // copy constructor
  95   fMUONTrack        = new AliMUONTrack(*(track.fMUONTrack));
  96   if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
  97   fMFTClusters      = new TClonesArray(*(track.fMFTClusters));
  98   fMFTClusters->SetOwner(kTRUE);
  99   for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
 100   for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
 101     fParentMCLabel[iParent] = (track.fParentMCLabel)[iParent];
 102     fParentPDGCode[iParent] = (track.fParentPDGCode)[iParent];
 103   }
 104
 105 }
 106
 107 //====================================================================================================================================================
 108
 109 AliMuonForwardTrack& AliMuonForwardTrack::operator=(const AliMuonForwardTrack& track) {
 110
 111   // Asignment operator
 112
 113   // check assignement to self
 114   if (this == &track) return *this;
 115
 116   // base class assignement
 117   AliMUONTrack::operator=(track);
 118
 119   // clear memory
 120   Clear("");
 121
 122   fMUONTrack        = new AliMUONTrack(*(track.fMUONTrack));
 123   if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
 124   fMFTClusters      = new TClonesArray(*(track.fMFTClusters));
 125   fMFTClusters->SetOwner(kTRUE);
 126   fNWrongClustersMC = track.fNWrongClustersMC;
 127   fTrackMCId        = track.fTrackMCId;
 128
 129   for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
 130   for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
 131     fParentMCLabel[iParent] = (track.fParentMCLabel)[iParent];
 132     fParentPDGCode[iParent] = (track.fParentPDGCode)[iParent];
 133   }
 134
 135   return *this;
 136
 137 }
 138
 139 //====================================================================================================================================================
 140
 141 void AliMuonForwardTrack::Clear(const Option_t* /*opt*/) {
 142
 143   // Clear arrays
 144   fMFTClusters -> Delete();
 145   delete fMFTClusters; fMFTClusters = 0x0;
 146   delete fMUONTrack;   fMUONTrack   = 0x0;
 147   delete fMCTrackRef;  fMCTrackRef  = 0x0;
 148
 149 }
 150
 151 //====================================================================================================================================================
 152
 153 AliMuonForwardTrack::~AliMuonForwardTrack() {
 154
 155   delete fMUONTrack;
 156   delete fMCTrackRef;
 157   fMFTClusters -> Delete();
 158   delete fMFTClusters;
 159
 160 }
 161
 162 //====================================================================================================================================================
 163
 164 void AliMuonForwardTrack::SetMUONTrack(AliMUONTrack *MUONTrack) {
 165
 166   if (fMUONTrack) {
 167     AliInfo("fMUONTrack already exists, nothing will be done");
 168     return;
 169   }
 170
 171   fMUONTrack = MUONTrack;
 172   SetGlobalChi2(fMUONTrack->GetGlobalChi2());
 173
 174 }
 175
 176 //====================================================================================================================================================
 177
 178 void AliMuonForwardTrack::SetMCTrackRef(TParticle *MCTrackRef) {
 179
 180   if (fMCTrackRef) {
 181     AliInfo("fMCTrackRef already exists, nothing will be done");
 182     return;
 183   }
 184
 185   fMCTrackRef = MCTrackRef;
 186
 187 }
 188
 189 //====================================================================================================================================================
 190
 191 void AliMuonForwardTrack::AddTrackParamAtMFTCluster(AliMUONTrackParam &trackParam, AliMFTCluster &mftCluster) {
 192
 193   AliDebug(1, Form("Before adding: this->fMFTClusters=%p has %d entries", this->fMFTClusters, this->fMFTClusters->GetEntries()));
 194   Int_t iMFTCluster = this->fMFTClusters->GetEntries();
 195   AliDebug(1, Form("mftCluster->GetX() = %f  mftCluster->GetY() = %f  mftCluster->GetErrX() = %f  cmftCluster->GetErrY() = %f",
 196            mftCluster.GetX(), mftCluster.GetY(), mftCluster.GetErrX(), mftCluster.GetErrY()));
 197   AliMUONVCluster *muonCluster = (AliMUONVCluster*) mftCluster.CreateMUONCluster();
 198   AliDebug(1, Form("Created MUON cluster %p", muonCluster));
 199   trackParam.SetUniqueID(iMFTCluster);    // we profit of this slot to store the reference to the corresponding MFTCluster
 200   AliDebug(1, Form("Now adding muonCluster %p and trackParam %p",muonCluster, &trackParam));
 201   AddTrackParamAtCluster(trackParam, *muonCluster, kTRUE);
 202   // we pass the parameters this->GetTrackParamAtCluster()->First() to the Kalman Filter algorithm: they will be updated!!
 203   Double_t chi2Kalman = RunKalmanFilter(*(AliMUONTrackParam*)(GetTrackParamAtCluster()->First()));
 204   AliDebug(1, Form("Adding Kalman chi2 = %f to global chi2 = %f", chi2Kalman, GetGlobalChi2()));
 205   Double_t newGlobalChi2 = GetGlobalChi2() + chi2Kalman;
 206   mftCluster.SetLocalChi2(chi2Kalman);
 207   mftCluster.SetTrackChi2(newGlobalChi2);
 208   new ((*(this->fMFTClusters))[iMFTCluster]) AliMFTCluster(mftCluster);
 209   AliDebug(1, Form("GetTrackParamAtCluster() = %p  has %d entries while this->fMFTClusters=%p has %d entries",
 210                    GetTrackParamAtCluster(), GetTrackParamAtCluster()->GetEntries(), this->fMFTClusters, this->fMFTClusters->GetEntries()));
 211   AliDebug(1, Form("muonCluster->GetZ() = %f, trackParam->GetZ() = %f",muonCluster->GetZ(), trackParam.GetZ()));
 212   SetGlobalChi2(newGlobalChi2);
 213   AliDebug(1, Form("New global chi2 = %f", GetGlobalChi2()));
 214   ((AliMUONTrackParam*) GetTrackParamAtCluster()->First())->SetTrackChi2(newGlobalChi2);
 215   for (Int_t iPar=0; iPar<GetTrackParamAtCluster()->GetEntries(); iPar++) {
 216     AliDebug(1, Form("GetTrackParamAtCluster()->At(%d)->GetClusterPtr() = %p",
 217                      iPar, ((AliMUONTrackParam*)(GetTrackParamAtCluster()->At(iPar)))->GetClusterPtr()));
 218   }
 219
 220 }
 221
 222 //====================================================================================================================================================
 223
 224 AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMUONCluster(Int_t iMUONCluster) {
 225
 226   if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
 227     AliError("Invalid MUON cluster index. NULL pointer will be returned");
 228     return NULL;
 229   }
 230
 231   AliMUONTrackParam *trackParam = (AliMUONTrackParam*) fMUONTrack->GetTrackParamAtCluster()->At(iMUONCluster);
 232
 233   return trackParam;
 234
 235 }
 236
 237 //====================================================================================================================================================
 238
 239 AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMFTCluster(Int_t iMFTCluster) {
 240
 241   if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
 242     AliError("Invalid MFT cluster index. NULL pointer will be returned");
 243     return NULL;
 244   }
 245
 246   AliMUONTrackParam *trackParam = (AliMUONTrackParam*) GetTrackParamAtCluster()->At(iMFTCluster);
 247
 248   return trackParam;
 249
 250 }
 251
 252 //====================================================================================================================================================
 253
 254 AliMUONVCluster* AliMuonForwardTrack::GetMUONCluster(Int_t iMUONCluster) {
 255
 256   if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
 257     AliError("Invalid MUON cluster index. NULL pointer will be returned");
 258     return NULL;
 259   }
 260
 261   AliMUONTrackParam *trackParam = GetTrackParamAtMUONCluster(iMUONCluster);
 262   AliMUONVCluster *muonCluster = trackParam->GetClusterPtr();
 263
 264   return muonCluster;
 265
 266 }
 267
 268 //====================================================================================================================================================
 269
 270 AliMFTCluster* AliMuonForwardTrack::GetMFTCluster(Int_t iMFTCluster) {
 271
 272   if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
 273     AliError(Form("Invalid MFT cluster index (%d). GetNMFTClusters()=%d. NULL pointer will be returned", iMFTCluster, GetNMFTClusters()));
 274     return NULL;
 275   }
 276
 277   AliMUONTrackParam *trackParam = GetTrackParamAtMFTCluster(iMFTCluster);
 278   AliMFTCluster *mftCluster = (AliMFTCluster*) this->fMFTClusters->At(trackParam->GetUniqueID());
 279
 280   return mftCluster;
 281
 282 }
 283
 284 //====================================================================================================================================================
 285
 286 Double_t AliMuonForwardTrack::RunKalmanFilter(AliMUONTrackParam &trackParamAtCluster) {
 287
 288   AliDebug(1, Form("Running Kalman filter for parameters %p (z = %f) and cluster %p (z = %f)",
 289                    &trackParamAtCluster, trackParamAtCluster.GetZ(), trackParamAtCluster.GetClusterPtr(), trackParamAtCluster.GetClusterPtr()->GetZ()));
 290
 291   // Compute new track parameters and their covariances including new cluster using kalman filter
 292   // return the *additional* track chi2
 293
 294   // Get actual track parameters (p)
 295   TMatrixD param(trackParamAtCluster.GetParameters());
 296
 297   // Get new cluster parameters (m)
 298   AliMUONVCluster *cluster = trackParamAtCluster.GetClusterPtr();
 299   AliDebug(1, Form("cluster->GetX() = %f  cluster->GetY() = %f  cluster->GetErrX() = %f  cluster->GetErrY() = %f",
 300                    cluster->GetX(), cluster->GetY(), cluster->GetErrX(), cluster->GetErrY()));
 301   TMatrixD clusterParam(5,1);
 302   clusterParam.Zero();
 303   clusterParam(0,0) = cluster->GetX();
 304   clusterParam(2,0) = cluster->GetY();
 305
 306   // Compute the current parameter weight (W)
 307   TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
 308   if (paramWeight.Determinant() != 0) {
 309     paramWeight.Invert();
 310   } else {
 311     Warning("RunKalmanFilter"," Determinant = 0");
 312     return 1.e10;
 313   }
 314
 315   // Compute the new cluster weight (U)
 316   TMatrixD clusterWeight(5,5);
 317   clusterWeight.Zero();
 318   clusterWeight(0,0) = 1. / cluster->GetErrX2();
 319   clusterWeight(2,2) = 1. / cluster->GetErrY2();
 320
 321   // Compute the new parameters covariance matrix ( (W+U)^-1 )
 322   TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
 323   if (newParamCov.Determinant() != 0) {
 324     newParamCov.Invert();
 325   }
 326   else {
 327     Warning("RunKalmanFilter"," Determinant = 0");
 328     return 1.e10;
 329   }
 330
 331   // Save the new parameters covariance matrix
 332   trackParamAtCluster.SetCovariances(newParamCov);
 333
 334   // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
 335   TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
 336   TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp);    // U(m-p)
 337   TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
 338   newParam += param;                                   // ((W+U)^-1)U(m-p) + p
 339
 340   // Save the new parameters
 341   trackParamAtCluster.SetParameters(newParam);
 342
 343   // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
 344   tmp = newParam; // p'
 345   tmp -= param;   // (p'-p)
 346   TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp);           // W(p'-p)
 347   TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
 348   tmp = newParam;      // p'
 349   tmp -= clusterParam; // (p'-m)
 350   TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp);            // U(p'-m)
 351   addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
 352
 353   AliDebug(1,Form("Adding Kalman chi2 = %f",addChi2Track(0,0)));
 354
 355   return addChi2Track(0,0);
 356
 357 }
 358
 359 //====================================================================================================================================================
 360
 361 Double_t AliMuonForwardTrack::GetWeightedOffset(Double_t x, Double_t y, Double_t z) {
 362
 363   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
 364   AliMUONTrackExtrap::ExtrapToZCov(param, z);
 365
 366   TMatrixD cov(5,5);
 367   cov = param->GetCovariances();
 368
 369   TMatrixD covCoordinates(2,2);
 370   covCoordinates(0,0) = cov(0,0);
 371   covCoordinates(0,1) = cov(0,2);
 372   covCoordinates(1,0) = cov(2,0);
 373   covCoordinates(1,1) = cov(2,2);
 374
 375   TMatrixD covCoordinatesInverse = covCoordinates.Invert();
 376
 377   Double_t dX = param->GetNonBendingCoor() - x;
 378   Double_t dY = param->GetBendingCoor()    - y;
 379
 380   Double_t weightedOffset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) +
 381                                              dY*dY*covCoordinatesInverse(1,1) +
 382                                              2.*dX*dY*covCoordinatesInverse(0,1)));
 383
 384   return weightedOffset;
 385
 386 }
 387
 388 //====================================================================================================================================================
 389
 390 Double_t AliMuonForwardTrack::GetOffsetX(Double_t x, Double_t z) {
 391
 392   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
 393   AliMUONTrackExtrap::ExtrapToZCov(param, z);
 394   Double_t dX = param->GetNonBendingCoor() - x;
 395   return dX;
 396
 397 }
 398
 399 //====================================================================================================================================================
 400
 401 Double_t AliMuonForwardTrack::GetOffsetY(Double_t y, Double_t z) {
 402
 403   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
 404   AliMUONTrackExtrap::ExtrapToZCov(param, z);
 405   Double_t dY = param->GetBendingCoor() - y;
 406   return dY;
 407
 408 }
 409
 410 //====================================================================================================================================================
 411
 412 Double_t AliMuonForwardTrack::GetOffset(Double_t x, Double_t y, Double_t z) {
 413
 414   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
 415   AliMUONTrackExtrap::ExtrapToZCov(param, z);
 416   Double_t dX = param->GetNonBendingCoor() - x;
 417   Double_t dY = param->GetBendingCoor()    - y;
 418   return TMath::Sqrt(dX*dX + dY*dY);
 419
 420 }
 421
 422 //====================================================================================================================================================
 423
 424 Bool_t AliMuonForwardTrack::IsFromResonance() {
 425
 426   Bool_t result = kFALSE;
 427
 428   if ( GetParentPDGCode(0) ==    113 ||
 429        GetParentPDGCode(0) ==    221 ||
 430        GetParentPDGCode(0) ==    223 ||
 431        GetParentPDGCode(0) ==    331 ||
 432        GetParentPDGCode(0) ==    333 ||
 433        GetParentPDGCode(0) ==    443 ||
 434        GetParentPDGCode(0) == 100443 ||
 435        GetParentPDGCode(0) ==    553 ||
 436        GetParentPDGCode(0) == 100553 ) result = kTRUE;
 437
 438   if (result) AliDebug(1, Form("Muon comes from a resonance %d", GetParentPDGCode(0)));
 439
 440   return result;
 441
 442 }
 443
 444 //====================================================================================================================================================
 445
 446 Bool_t AliMuonForwardTrack::IsFromCharm() {
 447
 448   Bool_t result = kFALSE;
 449
 450   if ( GetParentPDGCode(0) ==   411 ||
 451        GetParentPDGCode(0) ==   421 ||
 452        GetParentPDGCode(0) == 10411 ||
 453        GetParentPDGCode(0) == 10421 ||
 454        GetParentPDGCode(0) ==   413 ||
 455        GetParentPDGCode(0) ==   423 ||
 456        GetParentPDGCode(0) == 10413 ||
 457        GetParentPDGCode(0) == 10423 ||
 458        GetParentPDGCode(0) == 20413 ||
 459        GetParentPDGCode(0) == 20423 ||
 460        GetParentPDGCode(0) ==   415 ||
 461        GetParentPDGCode(0) ==   425 ||
 462        GetParentPDGCode(0) ==   431 ||
 463        GetParentPDGCode(0) == 10431 ||
 464        GetParentPDGCode(0) ==   433 ||
 465        GetParentPDGCode(0) == 10433 ||
 466        GetParentPDGCode(0) == 20433 ||
 467        GetParentPDGCode(0) ==   435 ) result = kTRUE;
 468
 469   if (result) AliDebug(1, Form("Muon comes from a charmed hadron %d", GetParentPDGCode(0)));
 470
 471   return result;
 472
 473 }
 474
 475 //====================================================================================================================================================
 476
 477 Bool_t AliMuonForwardTrack::IsFromBeauty() {
 478
 479   Bool_t result = kFALSE;
 480
 481   if ( GetParentPDGCode(0) ==   511 ||
 482        GetParentPDGCode(0) ==   521 ||
 483        GetParentPDGCode(0) == 10511 ||
 484        GetParentPDGCode(0) == 10521 ||
 485        GetParentPDGCode(0) ==   513 ||
 486        GetParentPDGCode(0) ==   523 ||
 487        GetParentPDGCode(0) == 10513 ||
 488        GetParentPDGCode(0) == 10523 ||
 489        GetParentPDGCode(0) == 20513 ||
 490        GetParentPDGCode(0) == 20523 ||
 491        GetParentPDGCode(0) ==   515 ||
 492        GetParentPDGCode(0) ==   525 ||
 493        GetParentPDGCode(0) ==   531 ||
 494        GetParentPDGCode(0) == 10531 ||
 495        GetParentPDGCode(0) ==   533 ||
 496        GetParentPDGCode(0) == 10533 ||
 497        GetParentPDGCode(0) == 20533 ||
 498        GetParentPDGCode(0) ==   535 ||
 499        GetParentPDGCode(0) ==   541 ||
 500        GetParentPDGCode(0) == 10541 ||
 501        GetParentPDGCode(0) ==   543 ||
 502        GetParentPDGCode(0) == 10543 ||
 503        GetParentPDGCode(0) == 20543 ||
 504        GetParentPDGCode(0) ==   545 ) result = kTRUE;
 505
 506   if (result) AliDebug(1, Form("Muon comes from a beauty hadron %d", GetParentPDGCode(0)));
 507
 508   return result;
 509
 510 }
 511
 512 //====================================================================================================================================================
 513
 514 Bool_t AliMuonForwardTrack::IsFromBackground() {
 515
 516   Bool_t result = kFALSE;
 517
 518   if (!IsFromResonance() && !IsFromCharm() && !IsFromBeauty()) result = kTRUE;
 519
 520   if (result) AliDebug(1, Form("Muon comes from a background source %d", GetParentPDGCode(0)));
 521
 522   return result;
 523
 524 }
 525
 526 //====================================================================================================================================================