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