1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //====================================================================================================================================================
18 // Description of an ALICE muon forward track, combining the information of the Muon Spectrometer and the Muon Forward Tracker
20 // Contact author: antonio.uras@cern.ch
22 //====================================================================================================================================================
25 #include "AliMUONTrack.h"
26 #include "AliMFTCluster.h"
27 #include "AliMUONVCluster.h"
28 #include "AliMUONTrackParam.h"
29 #include "AliMUONTrackExtrap.h"
30 #include "TClonesArray.h"
32 #include "TParticle.h"
33 #include "AliMuonForwardTrack.h"
34 #include "AliMFTConstants.h"
36 ClassImp(AliMuonForwardTrack)
38 //====================================================================================================================================================
40 AliMuonForwardTrack::AliMuonForwardTrack():
45 fNWrongClustersMC(-1),
49 // default constructor
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;
56 fMFTClusters = new TClonesArray("AliMFTCluster");
57 fMFTClusters -> SetOwner(kTRUE);
61 //====================================================================================================================================================
63 AliMuonForwardTrack::AliMuonForwardTrack(AliMUONTrack *MUONTrack):
68 fNWrongClustersMC(-1),
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;
78 fMFTClusters = new TClonesArray("AliMFTCluster");
79 fMFTClusters -> SetOwner(kTRUE);
83 //====================================================================================================================================================
85 AliMuonForwardTrack::AliMuonForwardTrack(const AliMuonForwardTrack& track):
90 fNWrongClustersMC(track.fNWrongClustersMC),
91 fTrackMCId(track.fTrackMCId)
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];
107 //====================================================================================================================================================
109 AliMuonForwardTrack& AliMuonForwardTrack::operator=(const AliMuonForwardTrack& track) {
111 // Asignment operator
113 // check assignement to self
114 if (this == &track) return *this;
116 // base class assignement
117 AliMUONTrack::operator=(track);
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;
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];
139 //====================================================================================================================================================
141 void AliMuonForwardTrack::Clear(const Option_t* /*opt*/) {
144 delete fMUONTrack; fMUONTrack = 0x0;
145 delete fMCTrackRef; fMCTrackRef = 0x0;
146 delete fMFTClusters; fMFTClusters = 0x0;
150 //====================================================================================================================================================
152 AliMuonForwardTrack::~AliMuonForwardTrack() {
156 fMFTClusters -> Delete();
161 //====================================================================================================================================================
163 void AliMuonForwardTrack::SetMUONTrack(AliMUONTrack *MUONTrack) {
166 AliInfo("fMUONTrack already exists, nothing will be done");
170 fMUONTrack = MUONTrack;
171 SetGlobalChi2(fMUONTrack->GetGlobalChi2());
175 //====================================================================================================================================================
177 void AliMuonForwardTrack::SetMCTrackRef(TParticle *MCTrackRef) {
180 AliInfo("fMCTrackRef already exists, nothing will be done");
184 fMCTrackRef = MCTrackRef;
188 //====================================================================================================================================================
190 void AliMuonForwardTrack::AddTrackParamAtMFTCluster(AliMUONTrackParam &trackParam, AliMFTCluster &mftCluster) {
192 AliDebug(1, Form("Before adding: this->fMFTClusters=%p has %d entries", this->fMFTClusters, this->fMFTClusters->GetEntries()));
193 Int_t iMFTCluster = this->fMFTClusters->GetEntries();
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);
201 // we pass the parameters this->GetTrackParamAtCluster()->First() to the Kalman Filter algorithm: they will be updated!!
202 Double_t chi2Kalman = RunKalmanFilter(*(AliMUONTrackParam*)(GetTrackParamAtCluster()->First()));
203 AliDebug(1, Form("Adding Kalman chi2 = %f to global chi2 = %f", chi2Kalman, GetGlobalChi2()));
204 Double_t newGlobalChi2 = GetGlobalChi2() + chi2Kalman;
205 mftCluster.SetLocalChi2(chi2Kalman);
206 mftCluster.SetTrackChi2(newGlobalChi2);
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()));
210 AliDebug(1, Form("muonCluster->GetZ() = %f, trackParam->GetZ() = %f",muonCluster->GetZ(), trackParam.GetZ()));
211 SetGlobalChi2(newGlobalChi2);
212 AliDebug(1, Form("New global chi2 = %f", GetGlobalChi2()));
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()));
221 //====================================================================================================================================================
223 AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMUONCluster(Int_t iMUONCluster) {
225 if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
226 AliError("Invalid MUON cluster index. NULL pointer will be returned");
230 AliMUONTrackParam *trackParam = (AliMUONTrackParam*) fMUONTrack->GetTrackParamAtCluster()->At(iMUONCluster);
236 //====================================================================================================================================================
238 AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMFTCluster(Int_t iMFTCluster) {
240 if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
241 AliError("Invalid MFT cluster index. NULL pointer will be returned");
245 AliMUONTrackParam *trackParam = (AliMUONTrackParam*) GetTrackParamAtCluster()->At(iMFTCluster);
251 //====================================================================================================================================================
253 AliMUONVCluster* AliMuonForwardTrack::GetMUONCluster(Int_t iMUONCluster) {
255 if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
256 AliError("Invalid MUON cluster index. NULL pointer will be returned");
260 AliMUONTrackParam *trackParam = GetTrackParamAtMUONCluster(iMUONCluster);
261 AliMUONVCluster *muonCluster = trackParam->GetClusterPtr();
267 //====================================================================================================================================================
269 AliMFTCluster* AliMuonForwardTrack::GetMFTCluster(Int_t iMFTCluster) {
271 if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
272 AliError(Form("Invalid MFT cluster index (%d). GetNMFTClusters()=%d. NULL pointer will be returned", iMFTCluster, GetNMFTClusters()));
276 AliMUONTrackParam *trackParam = GetTrackParamAtMFTCluster(iMFTCluster);
277 AliMFTCluster *mftCluster = (AliMFTCluster*) this->fMFTClusters->At(trackParam->GetUniqueID());
283 //====================================================================================================================================================
285 Double_t AliMuonForwardTrack::RunKalmanFilter(AliMUONTrackParam &trackParamAtCluster) {
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()));
290 // Compute new track parameters and their covariances including new cluster using kalman filter
291 // return the *additional* track chi2
293 // Get actual track parameters (p)
294 TMatrixD param(trackParamAtCluster.GetParameters());
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);
302 clusterParam(0,0) = cluster->GetX();
303 clusterParam(2,0) = cluster->GetY();
305 // Compute the current parameter weight (W)
306 TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
307 if (paramWeight.Determinant() != 0) {
308 paramWeight.Invert();
310 Warning("RunKalmanFilter"," Determinant = 0");
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();
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();
326 Warning("RunKalmanFilter"," Determinant = 0");
330 // Save the new parameters covariance matrix
331 trackParamAtCluster.SetCovariances(newParamCov);
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
339 // Save the new parameters
340 trackParamAtCluster.SetParameters(newParam);
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)
352 AliDebug(1,Form("Adding Kalman chi2 = %f",addChi2Track(0,0)));
354 return addChi2Track(0,0);
358 //====================================================================================================================================================
360 Double_t AliMuonForwardTrack::GetWeightedOffset(Double_t x, Double_t y, Double_t z) {
362 AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
363 AliMUONTrackExtrap::ExtrapToZCov(param, z);
366 cov = param->GetCovariances();
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);
374 TMatrixD covCoordinatesInverse = covCoordinates.Invert();
376 Double_t dX = param->GetNonBendingCoor() - x;
377 Double_t dY = param->GetBendingCoor() - y;
379 Double_t weightedOffset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) +
380 dY*dY*covCoordinatesInverse(1,1) +
381 2.*dX*dY*covCoordinatesInverse(0,1)));
383 return weightedOffset;
387 //====================================================================================================================================================
389 Double_t AliMuonForwardTrack::GetOffsetX(Double_t x, Double_t z) {
391 AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
392 AliMUONTrackExtrap::ExtrapToZCov(param, z);
393 Double_t dX = param->GetNonBendingCoor() - x;
398 //====================================================================================================================================================
400 Double_t AliMuonForwardTrack::GetOffsetY(Double_t y, Double_t z) {
402 AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
403 AliMUONTrackExtrap::ExtrapToZCov(param, z);
404 Double_t dY = param->GetBendingCoor() - y;
409 //====================================================================================================================================================
411 Double_t AliMuonForwardTrack::GetOffset(Double_t x, Double_t y, Double_t z) {
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(dX*dX + dY*dY);
421 //====================================================================================================================================================