Kinematics TLorentsVector added to 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 #include "TLorentzVector.h"
36 #include "TDatabasePDG.h"
37
38 ClassImp(AliMuonForwardTrack)
39
40 //====================================================================================================================================================
41
42 AliMuonForwardTrack::AliMuonForwardTrack():
43   AliMUONTrack(),
44   fMUONTrack(0),
45   fMCTrackRef(0),
46   fMFTClusters(0),
47   fNWrongClustersMC(-1),
48   fTrackMCId(-1),
49   fKinem(0,0,0,0),
50   fParamCovMatrix(5,5)
51 {
52
53   // default constructor
54   
55   for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = kFALSE;
56   for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
57     fParentMCLabel[iParent] = -1;
58     fParentPDGCode[iParent] =  0;
59   }
60   fMFTClusters = new TClonesArray("AliMFTCluster");
61   fMFTClusters -> SetOwner(kTRUE);
62
63 }
64
65 //====================================================================================================================================================
66
67 AliMuonForwardTrack::AliMuonForwardTrack(AliMUONTrack *MUONTrack):
68   AliMUONTrack(),
69   fMUONTrack(0),
70   fMCTrackRef(0),
71   fMFTClusters(0),
72   fNWrongClustersMC(-1),
73   fTrackMCId(-1),
74   fKinem(0,0,0,0),
75   fParamCovMatrix(5,5)
76 {
77
78   SetMUONTrack(MUONTrack);
79   for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = kFALSE;
80   for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
81     fParentMCLabel[iParent] = -1;
82     fParentPDGCode[iParent] =  0;
83   }
84   fMFTClusters = new TClonesArray("AliMFTCluster");
85   fMFTClusters -> SetOwner(kTRUE);
86
87 }
88
89 //====================================================================================================================================================
90
91 AliMuonForwardTrack::AliMuonForwardTrack(const AliMuonForwardTrack& track): 
92   AliMUONTrack(track),
93   fMUONTrack(0x0),
94   fMCTrackRef(0x0),
95   fMFTClusters(0x0),
96   fNWrongClustersMC(track.fNWrongClustersMC),
97   fTrackMCId(track.fTrackMCId),
98   fKinem(track.fKinem),
99   fParamCovMatrix(track.fParamCovMatrix)
100 {
101
102   // copy constructor
103   fMUONTrack        = new AliMUONTrack(*(track.fMUONTrack));
104   if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
105   fMFTClusters      = new TClonesArray(*(track.fMFTClusters));
106   fMFTClusters->SetOwner(kTRUE);
107   for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
108   for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
109     fParentMCLabel[iParent] = (track.fParentMCLabel)[iParent];
110     fParentPDGCode[iParent] = (track.fParentPDGCode)[iParent];
111   }
112   
113 }
114
115 //====================================================================================================================================================
116
117 AliMuonForwardTrack& AliMuonForwardTrack::operator=(const AliMuonForwardTrack& track) {
118
119   // Asignment operator
120
121   // check assignement to self
122   if (this == &track) return *this;
123
124   // base class assignement
125   AliMUONTrack::operator=(track);
126   
127   // clear memory
128   Clear("");
129   
130   fMUONTrack        = new AliMUONTrack(*(track.fMUONTrack));
131   if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
132   fMFTClusters      = new TClonesArray(*(track.fMFTClusters));
133   fMFTClusters->SetOwner(kTRUE);
134   fNWrongClustersMC = track.fNWrongClustersMC;
135   fTrackMCId        = track.fTrackMCId;
136   fKinem            = track.fKinem;
137   fParamCovMatrix   = track.fParamCovMatrix;
138
139   for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
140   for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
141     fParentMCLabel[iParent] = (track.fParentMCLabel)[iParent];
142     fParentPDGCode[iParent] = (track.fParentPDGCode)[iParent];
143   }
144   
145   return *this;
146
147 }
148
149 //====================================================================================================================================================
150
151 void AliMuonForwardTrack::Clear(const Option_t* /*opt*/) {
152
153   // Clear arrays
154   fMFTClusters -> Delete(); 
155   delete fMFTClusters; fMFTClusters = 0x0;
156   delete fMUONTrack;   fMUONTrack   = 0x0;
157   delete fMCTrackRef;  fMCTrackRef  = 0x0;
158   
159 }
160
161 //====================================================================================================================================================
162
163 AliMuonForwardTrack::~AliMuonForwardTrack() {
164
165   delete fMUONTrack;
166   delete fMCTrackRef;
167   fMFTClusters -> Delete();
168   delete fMFTClusters;
169   
170 }
171
172 //====================================================================================================================================================
173
174 void AliMuonForwardTrack::SetMUONTrack(AliMUONTrack *MUONTrack) {
175
176   if (fMUONTrack) {
177     AliInfo("fMUONTrack already exists, nothing will be done");
178     return;
179   }
180
181   fMUONTrack = MUONTrack;
182   SetGlobalChi2(fMUONTrack->GetGlobalChi2());
183   
184 }
185
186 //====================================================================================================================================================
187
188 void AliMuonForwardTrack::SetMCTrackRef(TParticle *MCTrackRef) {
189
190   if (fMCTrackRef) {
191     AliInfo("fMCTrackRef already exists, nothing will be done");
192     return;
193   }
194
195   fMCTrackRef = MCTrackRef;
196   
197 }
198
199 //====================================================================================================================================================
200
201 void AliMuonForwardTrack::AddTrackParamAtMFTCluster(AliMUONTrackParam &trackParam, AliMFTCluster &mftCluster) {
202
203   AliDebug(1, Form("Before adding: this->fMFTClusters=%p has %d entries", this->fMFTClusters, this->fMFTClusters->GetEntries()));
204   Int_t iMFTCluster = this->fMFTClusters->GetEntries();
205   AliDebug(1, Form("mftCluster->GetX() = %f  mftCluster->GetY() = %f  mftCluster->GetErrX() = %f  cmftCluster->GetErrY() = %f", 
206            mftCluster.GetX(), mftCluster.GetY(), mftCluster.GetErrX(), mftCluster.GetErrY()));
207   AliMUONVCluster *muonCluster = (AliMUONVCluster*) mftCluster.CreateMUONCluster();
208   AliDebug(1, Form("Created MUON cluster %p", muonCluster));
209   trackParam.SetUniqueID(iMFTCluster);    // we profit of this slot to store the reference to the corresponding MFTCluster
210   AliDebug(1, Form("Now adding muonCluster %p and trackParam %p",muonCluster, &trackParam));
211   AddTrackParamAtCluster(trackParam, *muonCluster, kTRUE);
212   // we pass the parameters this->GetTrackParamAtCluster()->First() to the Kalman Filter algorithm: they will be updated!!
213   Double_t chi2Kalman = RunKalmanFilter(*(AliMUONTrackParam*)(GetTrackParamAtCluster()->First()));
214   AliDebug(1, Form("Adding Kalman chi2 = %f to global chi2 = %f", chi2Kalman, GetGlobalChi2()));
215   Double_t newGlobalChi2 = GetGlobalChi2() + chi2Kalman;
216   mftCluster.SetLocalChi2(chi2Kalman);
217   mftCluster.SetTrackChi2(newGlobalChi2);
218   new ((*(this->fMFTClusters))[iMFTCluster]) AliMFTCluster(mftCluster);
219   AliDebug(1, Form("GetTrackParamAtCluster() = %p  has %d entries while this->fMFTClusters=%p has %d entries",
220                    GetTrackParamAtCluster(), GetTrackParamAtCluster()->GetEntries(), this->fMFTClusters, this->fMFTClusters->GetEntries()));
221   AliDebug(1, Form("muonCluster->GetZ() = %f, trackParam->GetZ() = %f",muonCluster->GetZ(), trackParam.GetZ()));
222   SetGlobalChi2(newGlobalChi2);
223   AliDebug(1, Form("New global chi2 = %f", GetGlobalChi2()));
224   ((AliMUONTrackParam*) GetTrackParamAtCluster()->First())->SetTrackChi2(newGlobalChi2);
225   for (Int_t iPar=0; iPar<GetTrackParamAtCluster()->GetEntries(); iPar++) {
226     AliDebug(1, Form("GetTrackParamAtCluster()->At(%d)->GetClusterPtr() = %p",
227                      iPar, ((AliMUONTrackParam*)(GetTrackParamAtCluster()->At(iPar)))->GetClusterPtr()));
228   }
229
230 }
231
232 //====================================================================================================================================================
233
234 AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMUONCluster(Int_t iMUONCluster) {
235
236   if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
237     AliError("Invalid MUON cluster index. NULL pointer will be returned");
238     return NULL;
239   }
240
241   AliMUONTrackParam *trackParam = (AliMUONTrackParam*) fMUONTrack->GetTrackParamAtCluster()->At(iMUONCluster); 
242
243   return trackParam;
244
245 }
246
247 //====================================================================================================================================================
248
249 AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMFTCluster(Int_t iMFTCluster) {
250
251   if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
252     AliError("Invalid MFT cluster index. NULL pointer will be returned");
253     return NULL;
254   }
255
256   AliMUONTrackParam *trackParam = (AliMUONTrackParam*) GetTrackParamAtCluster()->At(iMFTCluster); 
257
258   return trackParam;
259
260 }
261
262 //====================================================================================================================================================
263
264 AliMUONVCluster* AliMuonForwardTrack::GetMUONCluster(Int_t iMUONCluster) {
265
266   if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
267     AliError("Invalid MUON cluster index. NULL pointer will be returned");
268     return NULL;
269   }
270
271   AliMUONTrackParam *trackParam = GetTrackParamAtMUONCluster(iMUONCluster);
272   AliMUONVCluster *muonCluster = trackParam->GetClusterPtr();
273
274   return muonCluster;
275
276 }
277
278 //====================================================================================================================================================
279
280 AliMFTCluster* AliMuonForwardTrack::GetMFTCluster(Int_t iMFTCluster) {
281
282   if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
283     AliError(Form("Invalid MFT cluster index (%d). GetNMFTClusters()=%d. NULL pointer will be returned", iMFTCluster, GetNMFTClusters()));
284     return NULL;
285   }
286
287   AliMUONTrackParam *trackParam = GetTrackParamAtMFTCluster(iMFTCluster);
288   AliMFTCluster *mftCluster = (AliMFTCluster*) this->fMFTClusters->At(trackParam->GetUniqueID());
289
290   return mftCluster;
291
292 }
293
294 //====================================================================================================================================================
295
296 Double_t AliMuonForwardTrack::RunKalmanFilter(AliMUONTrackParam &trackParamAtCluster) {
297
298   AliDebug(1, Form("Running Kalman filter for parameters %p (z = %f) and cluster %p (z = %f)", 
299                    &trackParamAtCluster, trackParamAtCluster.GetZ(), trackParamAtCluster.GetClusterPtr(), trackParamAtCluster.GetClusterPtr()->GetZ()));
300   
301   // Compute new track parameters and their covariances including new cluster using kalman filter
302   // return the *additional* track chi2
303   
304   // Get actual track parameters (p)
305   TMatrixD param(trackParamAtCluster.GetParameters());
306   
307   // Get new cluster parameters (m)
308   AliMUONVCluster *cluster = trackParamAtCluster.GetClusterPtr();
309   AliDebug(1, Form("cluster->GetX() = %f  cluster->GetY() = %f  cluster->GetErrX() = %f  cluster->GetErrY() = %f", 
310                    cluster->GetX(), cluster->GetY(), cluster->GetErrX(), cluster->GetErrY()));
311   TMatrixD clusterParam(5,1);
312   clusterParam.Zero();
313   clusterParam(0,0) = cluster->GetX();
314   clusterParam(2,0) = cluster->GetY();
315
316   // Compute the current parameter weight (W)
317   TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
318   if (paramWeight.Determinant() != 0) {
319     paramWeight.Invert();
320   } else {
321     Warning("RunKalmanFilter"," Determinant = 0");
322     return 1.e10;
323   }
324   
325   // Compute the new cluster weight (U)
326   TMatrixD clusterWeight(5,5);
327   clusterWeight.Zero();
328   clusterWeight(0,0) = 1. / cluster->GetErrX2();
329   clusterWeight(2,2) = 1. / cluster->GetErrY2();
330
331   // Compute the new parameters covariance matrix ( (W+U)^-1 )
332   TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
333   if (newParamCov.Determinant() != 0) {
334     newParamCov.Invert();
335   } 
336   else {
337     Warning("RunKalmanFilter"," Determinant = 0");
338     return 1.e10;
339   }
340   
341   // Save the new parameters covariance matrix
342   trackParamAtCluster.SetCovariances(newParamCov);
343   
344   // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
345   TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
346   TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp);    // U(m-p)
347   TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
348   newParam += param;                                   // ((W+U)^-1)U(m-p) + p
349
350   // Save the new parameters
351   trackParamAtCluster.SetParameters(newParam);
352   
353   // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
354   tmp = newParam; // p'
355   tmp -= param;   // (p'-p)
356   TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp);           // W(p'-p)
357   TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
358   tmp = newParam;      // p'
359   tmp -= clusterParam; // (p'-m)
360   TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp);            // U(p'-m)
361   addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
362   
363   AliDebug(1,Form("Adding Kalman chi2 = %f",addChi2Track(0,0)));
364
365   return addChi2Track(0,0);
366   
367 }
368
369 //====================================================================================================================================================
370
371 Double_t AliMuonForwardTrack::GetWeightedOffset(Double_t x, Double_t y, Double_t z) {
372
373   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
374   AliMUONTrackExtrap::ExtrapToZCov(param, z);
375
376   TMatrixD cov(5,5);
377   cov = param->GetCovariances();
378
379   TMatrixD covCoordinates(2,2);
380   covCoordinates(0,0) = cov(0,0);
381   covCoordinates(0,1) = cov(0,2);
382   covCoordinates(1,0) = cov(2,0);
383   covCoordinates(1,1) = cov(2,2);
384   
385   TMatrixD covCoordinatesInverse = covCoordinates.Invert();
386
387   Double_t dX = param->GetNonBendingCoor() - x;
388   Double_t dY = param->GetBendingCoor()    - y;
389   
390   Double_t weightedOffset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) + 
391                                              dY*dY*covCoordinatesInverse(1,1) + 
392                                              2.*dX*dY*covCoordinatesInverse(0,1)));
393
394   return weightedOffset;
395
396 }
397
398 //====================================================================================================================================================
399
400 Double_t AliMuonForwardTrack::GetOffsetX(Double_t x, Double_t z) {
401
402   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
403   AliMUONTrackExtrap::ExtrapToZCov(param, z);
404   Double_t dX = param->GetNonBendingCoor() - x;
405   return dX;
406
407 }
408
409 //====================================================================================================================================================
410
411 Double_t AliMuonForwardTrack::GetOffsetY(Double_t y, Double_t z) {
412
413   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
414   AliMUONTrackExtrap::ExtrapToZCov(param, z);
415   Double_t dY = param->GetBendingCoor() - y;
416   return dY;
417
418 }
419
420 //====================================================================================================================================================
421
422 Double_t AliMuonForwardTrack::GetOffset(Double_t x, Double_t y, Double_t z) {
423
424   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
425   AliMUONTrackExtrap::ExtrapToZCov(param, z);
426   Double_t dX = param->GetNonBendingCoor() - x;
427   Double_t dY = param->GetBendingCoor()    - y;
428   return TMath::Sqrt(dX*dX + dY*dY);
429
430 }
431
432 //====================================================================================================================================================
433
434 Bool_t AliMuonForwardTrack::IsFromResonance() {
435
436   Bool_t result = kFALSE;
437
438   if ( GetParentPDGCode(0) ==    113 ||
439        GetParentPDGCode(0) ==    221 ||
440        GetParentPDGCode(0) ==    223 ||
441        GetParentPDGCode(0) ==    331 ||
442        GetParentPDGCode(0) ==    333 ||
443        GetParentPDGCode(0) ==    443 ||
444        GetParentPDGCode(0) == 100443 ||
445        GetParentPDGCode(0) ==    553 ||
446        GetParentPDGCode(0) == 100553 ) result = kTRUE;
447   
448   if (result) AliDebug(1, Form("Muon comes from a resonance %d", GetParentPDGCode(0)));
449   
450   return result; 
451   
452 }
453
454 //====================================================================================================================================================
455
456 Bool_t AliMuonForwardTrack::IsFromCharm() {
457
458   Bool_t result = kFALSE;
459
460   if ( GetParentPDGCode(0) ==   411 ||
461        GetParentPDGCode(0) ==   421 ||
462        GetParentPDGCode(0) == 10411 ||
463        GetParentPDGCode(0) == 10421 ||
464        GetParentPDGCode(0) ==   413 ||
465        GetParentPDGCode(0) ==   423 ||
466        GetParentPDGCode(0) == 10413 ||
467        GetParentPDGCode(0) == 10423 ||
468        GetParentPDGCode(0) == 20413 ||
469        GetParentPDGCode(0) == 20423 ||
470        GetParentPDGCode(0) ==   415 ||
471        GetParentPDGCode(0) ==   425 ||
472        GetParentPDGCode(0) ==   431 ||
473        GetParentPDGCode(0) == 10431 ||
474        GetParentPDGCode(0) ==   433 ||
475        GetParentPDGCode(0) == 10433 ||
476        GetParentPDGCode(0) == 20433 ||
477        GetParentPDGCode(0) ==   435 ) result = kTRUE;
478   
479   if (result) AliDebug(1, Form("Muon comes from a charmed hadron %d", GetParentPDGCode(0)));
480   
481   return result; 
482   
483 }
484
485 //====================================================================================================================================================
486
487 Bool_t AliMuonForwardTrack::IsFromBeauty() {
488
489   Bool_t result = kFALSE;
490
491   if ( GetParentPDGCode(0) ==   511 ||
492        GetParentPDGCode(0) ==   521 ||
493        GetParentPDGCode(0) == 10511 ||
494        GetParentPDGCode(0) == 10521 ||
495        GetParentPDGCode(0) ==   513 ||
496        GetParentPDGCode(0) ==   523 ||
497        GetParentPDGCode(0) == 10513 ||
498        GetParentPDGCode(0) == 10523 ||
499        GetParentPDGCode(0) == 20513 ||
500        GetParentPDGCode(0) == 20523 ||
501        GetParentPDGCode(0) ==   515 ||
502        GetParentPDGCode(0) ==   525 ||
503        GetParentPDGCode(0) ==   531 ||
504        GetParentPDGCode(0) == 10531 ||
505        GetParentPDGCode(0) ==   533 ||
506        GetParentPDGCode(0) == 10533 ||
507        GetParentPDGCode(0) == 20533 ||
508        GetParentPDGCode(0) ==   535 ||
509        GetParentPDGCode(0) ==   541 ||
510        GetParentPDGCode(0) == 10541 ||
511        GetParentPDGCode(0) ==   543 ||
512        GetParentPDGCode(0) == 10543 ||
513        GetParentPDGCode(0) == 20543 ||
514        GetParentPDGCode(0) ==   545 ) result = kTRUE;
515   
516   if (result) AliDebug(1, Form("Muon comes from a beauty hadron %d", GetParentPDGCode(0)));
517   
518   return result; 
519   
520 }
521
522 //====================================================================================================================================================
523
524 Bool_t AliMuonForwardTrack::IsFromBackground() {
525
526   Bool_t result = kFALSE;
527
528   if (!IsFromResonance() && !IsFromCharm() && !IsFromBeauty()) result = kTRUE;
529
530   if (result) AliDebug(1, Form("Muon comes from a background source %d", GetParentPDGCode(0)));
531
532   return result;
533
534 }
535
536 //====================================================================================================================================================
537
538 void AliMuonForwardTrack::EvalKinem(Double_t z) {
539
540   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
541   AliMUONTrackExtrap::ExtrapToZCov(param, z);
542
543   Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
544   Double_t energy = TMath::Sqrt(param->P()*param->P() + mMu*mMu);
545   fKinem.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
546
547   TMatrixD cov(5,5);
548   fParamCovMatrix = param->GetCovariances();
549
550 }
551
552 //====================================================================================================================================================
553