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