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