Bugs fixed
[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 #include "AliMFTConstants.h"
35
36 ClassImp(AliMuonForwardTrack)
37
38 //====================================================================================================================================================
39
40 AliMuonForwardTrack::AliMuonForwardTrack():
41   AliMUONTrack(),
42   fMUONTrack(0),
43   fMCTrackRef(0),
44   fMFTClusters(0),
45   fNWrongClustersMC(-1),
46   fTrackMCId(-1)
47 {
48
49   // default constructor
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   }
56   fMFTClusters = new TClonesArray("AliMFTCluster");
57
58 }
59
60 //====================================================================================================================================================
61
62 AliMuonForwardTrack::AliMuonForwardTrack(AliMUONTrack *MUONTrack):
63   AliMUONTrack(),
64   fMUONTrack(0),
65   fMCTrackRef(0),
66   fMFTClusters(0),
67   fNWrongClustersMC(-1),
68   fTrackMCId(-1)
69 {
70
71   SetMUONTrack(MUONTrack);
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   }
77   fMFTClusters = new TClonesArray("AliMFTCluster");
78
79 }
80
81 //====================================================================================================================================================
82
83 AliMuonForwardTrack::AliMuonForwardTrack(const AliMuonForwardTrack& track): 
84   AliMUONTrack(track),
85   fMUONTrack(0x0),
86   fMCTrackRef(0x0),
87   fMFTClusters(0x0),
88   fNWrongClustersMC(track.fNWrongClustersMC),
89   fTrackMCId(track.fTrackMCId)
90 {
91
92   // copy constructor
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   }
101   
102 }
103
104 //====================================================================================================================================================
105
106 AliMuonForwardTrack& 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   
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;
123   fTrackMCId        = track.fTrackMCId;
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   }
130   
131   return *this;
132
133 }
134
135 //====================================================================================================================================================
136
137 AliMuonForwardTrack::~AliMuonForwardTrack() {
138
139   delete fMUONTrack;
140   delete fMCTrackRef;
141   delete fMFTClusters;
142   
143 }
144
145 //====================================================================================================================================================
146
147 void AliMuonForwardTrack::SetMUONTrack(AliMUONTrack *MUONTrack) {
148
149   if (fMUONTrack) {
150     AliInfo("fMUONTrack already exists, nothing will be done");
151     return;
152   }
153
154   fMUONTrack = MUONTrack;
155   SetGlobalChi2(fMUONTrack->GetGlobalChi2());
156   
157 }
158
159 //====================================================================================================================================================
160
161 void 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
174 void AliMuonForwardTrack::AddTrackParamAtMFTCluster(AliMUONTrackParam &trackParam, AliMFTCluster &mftCluster) {
175
176   AliDebug(1, Form("Before adding: this->fMFTClusters=%p has %d entries", this->fMFTClusters, this->fMFTClusters->GetEntries()));
177   Int_t iMFTCluster = this->fMFTClusters->GetEntries();
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);
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()));
187   AliDebug(1, Form("Adding Kalman chi2 = %f to global chi2 = %f", chi2Kalman, GetGlobalChi2()));
188   Double_t newGlobalChi2 = GetGlobalChi2() + chi2Kalman;
189   mftCluster.SetLocalChi2(chi2Kalman);
190   mftCluster.SetTrackChi2(newGlobalChi2);
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()));
194   AliDebug(1, Form("muonCluster->GetZ() = %f, trackParam->GetZ() = %f",muonCluster->GetZ(), trackParam.GetZ()));
195   SetGlobalChi2(newGlobalChi2);
196   AliDebug(1, Form("New global chi2 = %f", GetGlobalChi2()));
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
207 AliMUONTrackParam* 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
222 AliMUONTrackParam* 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
237 AliMUONVCluster* 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
253 AliMFTCluster* AliMuonForwardTrack::GetMFTCluster(Int_t iMFTCluster) {
254
255   if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
256     AliError(Form("Invalid MFT cluster index (%d). GetNMFTClusters()=%d. NULL pointer will be returned", iMFTCluster, GetNMFTClusters()));
257     return NULL;
258   }
259
260   AliMUONTrackParam *trackParam = GetTrackParamAtMFTCluster(iMFTCluster);
261   AliMFTCluster *mftCluster = (AliMFTCluster*) this->fMFTClusters->At(trackParam->GetUniqueID());
262
263   return mftCluster;
264
265 }
266
267 //====================================================================================================================================================
268
269 Double_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
344 Double_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
373 Double_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
384 Double_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
395 Double_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