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