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