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