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