[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrack.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//====================================================================================================================================================
//
//      Description of an ALICE muon forward track, combining the information of the Muon Spectrometer and the Muon Forward Tracker
//
//      Contact author: antonio.uras@cern.ch
//
//====================================================================================================================================================

#include "AliLog.h"
#include "AliMUONTrack.h"
#include "AliMFTCluster.h"
#include "AliMUONVCluster.h"
#include "AliMUONTrackParam.h"
#include "AliMUONTrackExtrap.h"
#include "TClonesArray.h"
#include "TMatrixD.h"
#include "TParticle.h"
#include "AliMuonForwardTrack.h"
#include "AliMFTConstants.h"

ClassImp(AliMuonForwardTrack)

//====================================================================================================================================================

AliMuonForwardTrack::AliMuonForwardTrack():
  AliMUONTrack(),
  fMUONTrack(0),
  fMCTrackRef(0),
  fMFTClusters(0),
  fNWrongClustersMC(-1),
  fTrackMCId(-1)
{

  // default constructor
  
  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = kFALSE;
  for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
    fParentMCLabel[iParent] = -1;
    fParentPDGCode[iParent] =  0;
  }
  fMFTClusters = new TClonesArray("AliMFTCluster");
  fMFTClusters -> SetOwner(kTRUE);

}

//====================================================================================================================================================

AliMuonForwardTrack::AliMuonForwardTrack(AliMUONTrack *MUONTrack):
  AliMUONTrack(),
  fMUONTrack(0),
  fMCTrackRef(0),
  fMFTClusters(0),
  fNWrongClustersMC(-1),
  fTrackMCId(-1)
{

  SetMUONTrack(MUONTrack);
  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = kFALSE;
  for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
    fParentMCLabel[iParent] = -1;
    fParentPDGCode[iParent] =  0;
  }
  fMFTClusters = new TClonesArray("AliMFTCluster");
  fMFTClusters -> SetOwner(kTRUE);

}

//====================================================================================================================================================

AliMuonForwardTrack::AliMuonForwardTrack(const AliMuonForwardTrack& track): 
  AliMUONTrack(track),
  fMUONTrack(0x0),
  fMCTrackRef(0x0),
  fMFTClusters(0x0),
  fNWrongClustersMC(track.fNWrongClustersMC),
  fTrackMCId(track.fTrackMCId)
{

  // copy constructor
  fMUONTrack        = new AliMUONTrack(*(track.fMUONTrack));
  if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
  fMFTClusters      = new TClonesArray(*(track.fMFTClusters));
  fMFTClusters->SetOwner(kTRUE);
  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
  for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
    fParentMCLabel[iParent] = (track.fParentMCLabel)[iParent];
    fParentPDGCode[iParent] = (track.fParentPDGCode)[iParent];
  }
  
}

//====================================================================================================================================================

AliMuonForwardTrack& AliMuonForwardTrack::operator=(const AliMuonForwardTrack& track) {

  // Asignment operator

  // check assignement to self
  if (this == &track) return *this;

  // base class assignement
  AliMUONTrack::operator=(track);
  
  // clear memory
  Clear("");
  
  fMUONTrack        = new AliMUONTrack(*(track.fMUONTrack));
  if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
  fMFTClusters      = new TClonesArray(*(track.fMFTClusters));
  fMFTClusters->SetOwner(kTRUE);
  fNWrongClustersMC = track.fNWrongClustersMC;
  fTrackMCId        = track.fTrackMCId;

  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
  for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
    fParentMCLabel[iParent] = (track.fParentMCLabel)[iParent];
    fParentPDGCode[iParent] = (track.fParentPDGCode)[iParent];
  }
  
  return *this;

}

//====================================================================================================================================================

void AliMuonForwardTrack::Clear(const Option_t* /*opt*/) {

  // Clear arrays
  delete fMUONTrack;   fMUONTrack   = 0x0;
  delete fMCTrackRef;  fMCTrackRef  = 0x0;
  delete fMFTClusters; fMFTClusters = 0x0;
  
}

//====================================================================================================================================================

AliMuonForwardTrack::~AliMuonForwardTrack() {

  delete fMUONTrack;
  delete fMCTrackRef;
  fMFTClusters -> Delete();
  delete fMFTClusters;
  
}

//====================================================================================================================================================

void AliMuonForwardTrack::SetMUONTrack(AliMUONTrack *MUONTrack) {

  if (fMUONTrack) {
    AliInfo("fMUONTrack already exists, nothing will be done");
    return;
  }

  fMUONTrack = MUONTrack;
  SetGlobalChi2(fMUONTrack->GetGlobalChi2());
  
}

//====================================================================================================================================================

void AliMuonForwardTrack::SetMCTrackRef(TParticle *MCTrackRef) {

  if (fMCTrackRef) {
    AliInfo("fMCTrackRef already exists, nothing will be done");
    return;
  }

  fMCTrackRef = MCTrackRef;
  
}

//====================================================================================================================================================

void AliMuonForwardTrack::AddTrackParamAtMFTCluster(AliMUONTrackParam &trackParam, AliMFTCluster &mftCluster) {

  AliDebug(1, Form("Before adding: this->fMFTClusters=%p has %d entries", this->fMFTClusters, this->fMFTClusters->GetEntries()));
  Int_t iMFTCluster = this->fMFTClusters->GetEntries();
  AliDebug(1, Form("mftCluster->GetX() = %f  mftCluster->GetY() = %f  mftCluster->GetErrX() = %f  cmftCluster->GetErrY() = %f", 
	   mftCluster.GetX(), mftCluster.GetY(), mftCluster.GetErrX(), mftCluster.GetErrY()));
  AliMUONVCluster *muonCluster = (AliMUONVCluster*) mftCluster.CreateMUONCluster();
  AliDebug(1, Form("Created MUON cluster %p", muonCluster));
  trackParam.SetUniqueID(iMFTCluster);    // we profit of this slot to store the reference to the corresponding MFTCluster
  AliDebug(1, Form("Now adding muonCluster %p and trackParam %p",muonCluster, &trackParam));
  AddTrackParamAtCluster(trackParam, *muonCluster, kTRUE);
  // we pass the parameters this->GetTrackParamAtCluster()->First() to the Kalman Filter algorithm: they will be updated!!
  Double_t chi2Kalman = RunKalmanFilter(*(AliMUONTrackParam*)(GetTrackParamAtCluster()->First()));
  AliDebug(1, Form("Adding Kalman chi2 = %f to global chi2 = %f", chi2Kalman, GetGlobalChi2()));
  Double_t newGlobalChi2 = GetGlobalChi2() + chi2Kalman;
  mftCluster.SetLocalChi2(chi2Kalman);
  mftCluster.SetTrackChi2(newGlobalChi2);
  new ((*(this->fMFTClusters))[iMFTCluster]) AliMFTCluster(mftCluster);
  AliDebug(1, Form("GetTrackParamAtCluster() = %p  has %d entries while this->fMFTClusters=%p has %d entries",
		   GetTrackParamAtCluster(), GetTrackParamAtCluster()->GetEntries(), this->fMFTClusters, this->fMFTClusters->GetEntries()));
  AliDebug(1, Form("muonCluster->GetZ() = %f, trackParam->GetZ() = %f",muonCluster->GetZ(), trackParam.GetZ()));
  SetGlobalChi2(newGlobalChi2);
  AliDebug(1, Form("New global chi2 = %f", GetGlobalChi2()));
  ((AliMUONTrackParam*) GetTrackParamAtCluster()->First())->SetTrackChi2(newGlobalChi2);
  for (Int_t iPar=0; iPar<GetTrackParamAtCluster()->GetEntries(); iPar++) {
    AliDebug(1, Form("GetTrackParamAtCluster()->At(%d)->GetClusterPtr() = %p",
		     iPar, ((AliMUONTrackParam*)(GetTrackParamAtCluster()->At(iPar)))->GetClusterPtr()));
  }

}

//====================================================================================================================================================

AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMUONCluster(Int_t iMUONCluster) {

  if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
    AliError("Invalid MUON cluster index. NULL pointer will be returned");
    return NULL;
  }

  AliMUONTrackParam *trackParam = (AliMUONTrackParam*) fMUONTrack->GetTrackParamAtCluster()->At(iMUONCluster); 

  return trackParam;

}

//====================================================================================================================================================

AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMFTCluster(Int_t iMFTCluster) {

  if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
    AliError("Invalid MFT cluster index. NULL pointer will be returned");
    return NULL;
  }

  AliMUONTrackParam *trackParam = (AliMUONTrackParam*) GetTrackParamAtCluster()->At(iMFTCluster); 

  return trackParam;

}

//====================================================================================================================================================

AliMUONVCluster* AliMuonForwardTrack::GetMUONCluster(Int_t iMUONCluster) {

  if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
    AliError("Invalid MUON cluster index. NULL pointer will be returned");
    return NULL;
  }

  AliMUONTrackParam *trackParam = GetTrackParamAtMUONCluster(iMUONCluster);
  AliMUONVCluster *muonCluster = trackParam->GetClusterPtr();

  return muonCluster;

}

//====================================================================================================================================================

AliMFTCluster* AliMuonForwardTrack::GetMFTCluster(Int_t iMFTCluster) {

  if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
    AliError(Form("Invalid MFT cluster index (%d). GetNMFTClusters()=%d. NULL pointer will be returned", iMFTCluster, GetNMFTClusters()));
    return NULL;
  }

  AliMUONTrackParam *trackParam = GetTrackParamAtMFTCluster(iMFTCluster);
  AliMFTCluster *mftCluster = (AliMFTCluster*) this->fMFTClusters->At(trackParam->GetUniqueID());

  return mftCluster;

}

//====================================================================================================================================================

Double_t AliMuonForwardTrack::RunKalmanFilter(AliMUONTrackParam &trackParamAtCluster) {

  AliDebug(1, Form("Running Kalman filter for parameters %p (z = %f) and cluster %p (z = %f)", 
		   &trackParamAtCluster, trackParamAtCluster.GetZ(), trackParamAtCluster.GetClusterPtr(), trackParamAtCluster.GetClusterPtr()->GetZ()));
  
  // Compute new track parameters and their covariances including new cluster using kalman filter
  // return the *additional* track chi2
  
  // Get actual track parameters (p)
  TMatrixD param(trackParamAtCluster.GetParameters());
  
  // Get new cluster parameters (m)
  AliMUONVCluster *cluster = trackParamAtCluster.GetClusterPtr();
  AliDebug(1, Form("cluster->GetX() = %f  cluster->GetY() = %f  cluster->GetErrX() = %f  cluster->GetErrY() = %f", 
		   cluster->GetX(), cluster->GetY(), cluster->GetErrX(), cluster->GetErrY()));
  TMatrixD clusterParam(5,1);
  clusterParam.Zero();
  clusterParam(0,0) = cluster->GetX();
  clusterParam(2,0) = cluster->GetY();

  // Compute the current parameter weight (W)
  TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
  if (paramWeight.Determinant() != 0) {
    paramWeight.Invert();
  } else {
    Warning("RunKalmanFilter"," Determinant = 0");
    return 1.e10;
  }
  
  // Compute the new cluster weight (U)
  TMatrixD clusterWeight(5,5);
  clusterWeight.Zero();
  clusterWeight(0,0) = 1. / cluster->GetErrX2();
  clusterWeight(2,2) = 1. / cluster->GetErrY2();

  // Compute the new parameters covariance matrix ( (W+U)^-1 )
  TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
  if (newParamCov.Determinant() != 0) {
    newParamCov.Invert();
  } 
  else {
    Warning("RunKalmanFilter"," Determinant = 0");
    return 1.e10;
  }
  
  // Save the new parameters covariance matrix
  trackParamAtCluster.SetCovariances(newParamCov);
  
  // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
  TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
  TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp);    // U(m-p)
  TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
  newParam += param;                                   // ((W+U)^-1)U(m-p) + p

  // Save the new parameters
  trackParamAtCluster.SetParameters(newParam);
  
  // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
  tmp = newParam; // p'
  tmp -= param;   // (p'-p)
  TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp);           // W(p'-p)
  TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
  tmp = newParam;      // p'
  tmp -= clusterParam; // (p'-m)
  TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp);            // U(p'-m)
  addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
  
  AliDebug(1,Form("Adding Kalman chi2 = %f",addChi2Track(0,0)));

  return addChi2Track(0,0);
  
}

//====================================================================================================================================================

Double_t AliMuonForwardTrack::GetWeightedOffset(Double_t x, Double_t y, Double_t z) {

  AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
  AliMUONTrackExtrap::ExtrapToZCov(param, z);

  TMatrixD cov(5,5);
  cov = param->GetCovariances();

  TMatrixD covCoordinates(2,2);
  covCoordinates(0,0) = cov(0,0);
  covCoordinates(0,1) = cov(0,2);
  covCoordinates(1,0) = cov(2,0);
  covCoordinates(1,1) = cov(2,2);
  
  TMatrixD covCoordinatesInverse = covCoordinates.Invert();

  Double_t dX = param->GetNonBendingCoor() - x;
  Double_t dY = param->GetBendingCoor()    - y;
  
  Double_t weightedOffset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) + 
					     dY*dY*covCoordinatesInverse(1,1) + 
					     2.*dX*dY*covCoordinatesInverse(0,1)));

  return weightedOffset;

}

//====================================================================================================================================================

Double_t AliMuonForwardTrack::GetOffsetX(Double_t x, Double_t z) {

  AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
  AliMUONTrackExtrap::ExtrapToZCov(param, z);
  Double_t dX = param->GetNonBendingCoor() - x;
  return dX;

}

//====================================================================================================================================================

Double_t AliMuonForwardTrack::GetOffsetY(Double_t y, Double_t z) {

  AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
  AliMUONTrackExtrap::ExtrapToZCov(param, z);
  Double_t dY = param->GetBendingCoor() - y;
  return dY;

}

//====================================================================================================================================================

Double_t AliMuonForwardTrack::GetOffset(Double_t x, Double_t y, Double_t z) {

  AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
  AliMUONTrackExtrap::ExtrapToZCov(param, z);
  Double_t dX = param->GetNonBendingCoor() - x;
  Double_t dY = param->GetBendingCoor()    - y;
  return TMath::Sqrt(dX*dX + dY*dY);

}

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