0a5897f015d7dc26301754ecba2221a6f70030e8
[u/mrichter/AliRoot.git] / 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
35 ClassImp(AliMuonForwardTrack)
36
37 //====================================================================================================================================================
38
39 AliMuonForwardTrack::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
54 AliMuonForwardTrack::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
69 AliMuonForwardTrack::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
83 AliMuonForwardTrack& 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
108 void 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
121 void 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
134 void 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
164 AliMUONTrackParam* 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
179 AliMUONTrackParam* 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
194 AliMUONVCluster* 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
210 AliMFTCluster* 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
226 Double_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
301 Double_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
330 Double_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
341 Double_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
352 Double_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