Methods added to distinguish muons from charm and beauty
[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   fMFTClusters -> Delete(); 
145   delete fMFTClusters; fMFTClusters = 0x0;
146   delete fMUONTrack;   fMUONTrack   = 0x0;
147   delete fMCTrackRef;  fMCTrackRef  = 0x0;
148   
149 }
150
151 //====================================================================================================================================================
152
153 AliMuonForwardTrack::~AliMuonForwardTrack() {
154
155   delete fMUONTrack;
156   delete fMCTrackRef;
157   fMFTClusters -> Delete();
158   delete fMFTClusters;
159   
160 }
161
162 //====================================================================================================================================================
163
164 void AliMuonForwardTrack::SetMUONTrack(AliMUONTrack *MUONTrack) {
165
166   if (fMUONTrack) {
167     AliInfo("fMUONTrack already exists, nothing will be done");
168     return;
169   }
170
171   fMUONTrack = MUONTrack;
172   SetGlobalChi2(fMUONTrack->GetGlobalChi2());
173   
174 }
175
176 //====================================================================================================================================================
177
178 void AliMuonForwardTrack::SetMCTrackRef(TParticle *MCTrackRef) {
179
180   if (fMCTrackRef) {
181     AliInfo("fMCTrackRef already exists, nothing will be done");
182     return;
183   }
184
185   fMCTrackRef = MCTrackRef;
186   
187 }
188
189 //====================================================================================================================================================
190
191 void AliMuonForwardTrack::AddTrackParamAtMFTCluster(AliMUONTrackParam &trackParam, AliMFTCluster &mftCluster) {
192
193   AliDebug(1, Form("Before adding: this->fMFTClusters=%p has %d entries", this->fMFTClusters, this->fMFTClusters->GetEntries()));
194   Int_t iMFTCluster = this->fMFTClusters->GetEntries();
195   AliDebug(1, Form("mftCluster->GetX() = %f  mftCluster->GetY() = %f  mftCluster->GetErrX() = %f  cmftCluster->GetErrY() = %f", 
196            mftCluster.GetX(), mftCluster.GetY(), mftCluster.GetErrX(), mftCluster.GetErrY()));
197   AliMUONVCluster *muonCluster = (AliMUONVCluster*) mftCluster.CreateMUONCluster();
198   AliDebug(1, Form("Created MUON cluster %p", muonCluster));
199   trackParam.SetUniqueID(iMFTCluster);    // we profit of this slot to store the reference to the corresponding MFTCluster
200   AliDebug(1, Form("Now adding muonCluster %p and trackParam %p",muonCluster, &trackParam));
201   AddTrackParamAtCluster(trackParam, *muonCluster, kTRUE);
202   // we pass the parameters this->GetTrackParamAtCluster()->First() to the Kalman Filter algorithm: they will be updated!!
203   Double_t chi2Kalman = RunKalmanFilter(*(AliMUONTrackParam*)(GetTrackParamAtCluster()->First()));
204   AliDebug(1, Form("Adding Kalman chi2 = %f to global chi2 = %f", chi2Kalman, GetGlobalChi2()));
205   Double_t newGlobalChi2 = GetGlobalChi2() + chi2Kalman;
206   mftCluster.SetLocalChi2(chi2Kalman);
207   mftCluster.SetTrackChi2(newGlobalChi2);
208   new ((*(this->fMFTClusters))[iMFTCluster]) AliMFTCluster(mftCluster);
209   AliDebug(1, Form("GetTrackParamAtCluster() = %p  has %d entries while this->fMFTClusters=%p has %d entries",
210                    GetTrackParamAtCluster(), GetTrackParamAtCluster()->GetEntries(), this->fMFTClusters, this->fMFTClusters->GetEntries()));
211   AliDebug(1, Form("muonCluster->GetZ() = %f, trackParam->GetZ() = %f",muonCluster->GetZ(), trackParam.GetZ()));
212   SetGlobalChi2(newGlobalChi2);
213   AliDebug(1, Form("New global chi2 = %f", GetGlobalChi2()));
214   ((AliMUONTrackParam*) GetTrackParamAtCluster()->First())->SetTrackChi2(newGlobalChi2);
215   for (Int_t iPar=0; iPar<GetTrackParamAtCluster()->GetEntries(); iPar++) {
216     AliDebug(1, Form("GetTrackParamAtCluster()->At(%d)->GetClusterPtr() = %p",
217                      iPar, ((AliMUONTrackParam*)(GetTrackParamAtCluster()->At(iPar)))->GetClusterPtr()));
218   }
219
220 }
221
222 //====================================================================================================================================================
223
224 AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMUONCluster(Int_t iMUONCluster) {
225
226   if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
227     AliError("Invalid MUON cluster index. NULL pointer will be returned");
228     return NULL;
229   }
230
231   AliMUONTrackParam *trackParam = (AliMUONTrackParam*) fMUONTrack->GetTrackParamAtCluster()->At(iMUONCluster); 
232
233   return trackParam;
234
235 }
236
237 //====================================================================================================================================================
238
239 AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMFTCluster(Int_t iMFTCluster) {
240
241   if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
242     AliError("Invalid MFT cluster index. NULL pointer will be returned");
243     return NULL;
244   }
245
246   AliMUONTrackParam *trackParam = (AliMUONTrackParam*) GetTrackParamAtCluster()->At(iMFTCluster); 
247
248   return trackParam;
249
250 }
251
252 //====================================================================================================================================================
253
254 AliMUONVCluster* AliMuonForwardTrack::GetMUONCluster(Int_t iMUONCluster) {
255
256   if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
257     AliError("Invalid MUON cluster index. NULL pointer will be returned");
258     return NULL;
259   }
260
261   AliMUONTrackParam *trackParam = GetTrackParamAtMUONCluster(iMUONCluster);
262   AliMUONVCluster *muonCluster = trackParam->GetClusterPtr();
263
264   return muonCluster;
265
266 }
267
268 //====================================================================================================================================================
269
270 AliMFTCluster* AliMuonForwardTrack::GetMFTCluster(Int_t iMFTCluster) {
271
272   if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
273     AliError(Form("Invalid MFT cluster index (%d). GetNMFTClusters()=%d. NULL pointer will be returned", iMFTCluster, GetNMFTClusters()));
274     return NULL;
275   }
276
277   AliMUONTrackParam *trackParam = GetTrackParamAtMFTCluster(iMFTCluster);
278   AliMFTCluster *mftCluster = (AliMFTCluster*) this->fMFTClusters->At(trackParam->GetUniqueID());
279
280   return mftCluster;
281
282 }
283
284 //====================================================================================================================================================
285
286 Double_t AliMuonForwardTrack::RunKalmanFilter(AliMUONTrackParam &trackParamAtCluster) {
287
288   AliDebug(1, Form("Running Kalman filter for parameters %p (z = %f) and cluster %p (z = %f)", 
289                    &trackParamAtCluster, trackParamAtCluster.GetZ(), trackParamAtCluster.GetClusterPtr(), trackParamAtCluster.GetClusterPtr()->GetZ()));
290   
291   // Compute new track parameters and their covariances including new cluster using kalman filter
292   // return the *additional* track chi2
293   
294   // Get actual track parameters (p)
295   TMatrixD param(trackParamAtCluster.GetParameters());
296   
297   // Get new cluster parameters (m)
298   AliMUONVCluster *cluster = trackParamAtCluster.GetClusterPtr();
299   AliDebug(1, Form("cluster->GetX() = %f  cluster->GetY() = %f  cluster->GetErrX() = %f  cluster->GetErrY() = %f", 
300                    cluster->GetX(), cluster->GetY(), cluster->GetErrX(), cluster->GetErrY()));
301   TMatrixD clusterParam(5,1);
302   clusterParam.Zero();
303   clusterParam(0,0) = cluster->GetX();
304   clusterParam(2,0) = cluster->GetY();
305
306   // Compute the current parameter weight (W)
307   TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
308   if (paramWeight.Determinant() != 0) {
309     paramWeight.Invert();
310   } else {
311     Warning("RunKalmanFilter"," Determinant = 0");
312     return 1.e10;
313   }
314   
315   // Compute the new cluster weight (U)
316   TMatrixD clusterWeight(5,5);
317   clusterWeight.Zero();
318   clusterWeight(0,0) = 1. / cluster->GetErrX2();
319   clusterWeight(2,2) = 1. / cluster->GetErrY2();
320
321   // Compute the new parameters covariance matrix ( (W+U)^-1 )
322   TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
323   if (newParamCov.Determinant() != 0) {
324     newParamCov.Invert();
325   } 
326   else {
327     Warning("RunKalmanFilter"," Determinant = 0");
328     return 1.e10;
329   }
330   
331   // Save the new parameters covariance matrix
332   trackParamAtCluster.SetCovariances(newParamCov);
333   
334   // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
335   TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
336   TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp);    // U(m-p)
337   TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
338   newParam += param;                                   // ((W+U)^-1)U(m-p) + p
339
340   // Save the new parameters
341   trackParamAtCluster.SetParameters(newParam);
342   
343   // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
344   tmp = newParam; // p'
345   tmp -= param;   // (p'-p)
346   TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp);           // W(p'-p)
347   TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
348   tmp = newParam;      // p'
349   tmp -= clusterParam; // (p'-m)
350   TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp);            // U(p'-m)
351   addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
352   
353   AliDebug(1,Form("Adding Kalman chi2 = %f",addChi2Track(0,0)));
354
355   return addChi2Track(0,0);
356   
357 }
358
359 //====================================================================================================================================================
360
361 Double_t AliMuonForwardTrack::GetWeightedOffset(Double_t x, Double_t y, Double_t z) {
362
363   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
364   AliMUONTrackExtrap::ExtrapToZCov(param, z);
365
366   TMatrixD cov(5,5);
367   cov = param->GetCovariances();
368
369   TMatrixD covCoordinates(2,2);
370   covCoordinates(0,0) = cov(0,0);
371   covCoordinates(0,1) = cov(0,2);
372   covCoordinates(1,0) = cov(2,0);
373   covCoordinates(1,1) = cov(2,2);
374   
375   TMatrixD covCoordinatesInverse = covCoordinates.Invert();
376
377   Double_t dX = param->GetNonBendingCoor() - x;
378   Double_t dY = param->GetBendingCoor()    - y;
379   
380   Double_t weightedOffset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) + 
381                                              dY*dY*covCoordinatesInverse(1,1) + 
382                                              2.*dX*dY*covCoordinatesInverse(0,1)));
383
384   return weightedOffset;
385
386 }
387
388 //====================================================================================================================================================
389
390 Double_t AliMuonForwardTrack::GetOffsetX(Double_t x, Double_t z) {
391
392   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
393   AliMUONTrackExtrap::ExtrapToZCov(param, z);
394   Double_t dX = param->GetNonBendingCoor() - x;
395   return dX;
396
397 }
398
399 //====================================================================================================================================================
400
401 Double_t AliMuonForwardTrack::GetOffsetY(Double_t y, Double_t z) {
402
403   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
404   AliMUONTrackExtrap::ExtrapToZCov(param, z);
405   Double_t dY = param->GetBendingCoor() - y;
406   return dY;
407
408 }
409
410 //====================================================================================================================================================
411
412 Double_t AliMuonForwardTrack::GetOffset(Double_t x, Double_t y, Double_t z) {
413
414   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
415   AliMUONTrackExtrap::ExtrapToZCov(param, z);
416   Double_t dX = param->GetNonBendingCoor() - x;
417   Double_t dY = param->GetBendingCoor()    - y;
418   return TMath::Sqrt(dX*dX + dY*dY);
419
420 }
421
422 //====================================================================================================================================================
423
424 Bool_t AliMuonForwardTrack::IsFromResonance() {
425
426   Bool_t result = kFALSE;
427
428   if ( GetParentPDGCode(0) ==    113 ||
429        GetParentPDGCode(0) ==    221 ||
430        GetParentPDGCode(0) ==    223 ||
431        GetParentPDGCode(0) ==    331 ||
432        GetParentPDGCode(0) ==    333 ||
433        GetParentPDGCode(0) ==    443 ||
434        GetParentPDGCode(0) == 100443 ||
435        GetParentPDGCode(0) ==    553 ||
436        GetParentPDGCode(0) == 100553 ) result = kTRUE;
437   
438   if (result) AliDebug(1, Form("Muon comes from a resonance %d", GetParentPDGCode(0)));
439   
440   return result; 
441   
442 }
443
444 //====================================================================================================================================================
445
446 Bool_t AliMuonForwardTrack::IsFromCharm() {
447
448   Bool_t result = kFALSE;
449
450   if ( GetParentPDGCode(0) ==   411 ||
451        GetParentPDGCode(0) ==   421 ||
452        GetParentPDGCode(0) == 10411 ||
453        GetParentPDGCode(0) == 10421 ||
454        GetParentPDGCode(0) ==   413 ||
455        GetParentPDGCode(0) ==   423 ||
456        GetParentPDGCode(0) == 10413 ||
457        GetParentPDGCode(0) == 10423 ||
458        GetParentPDGCode(0) == 20413 ||
459        GetParentPDGCode(0) == 20423 ||
460        GetParentPDGCode(0) ==   415 ||
461        GetParentPDGCode(0) ==   425 ||
462        GetParentPDGCode(0) ==   431 ||
463        GetParentPDGCode(0) == 10431 ||
464        GetParentPDGCode(0) ==   433 ||
465        GetParentPDGCode(0) == 10433 ||
466        GetParentPDGCode(0) == 20433 ||
467        GetParentPDGCode(0) ==   435 ) result = kTRUE;
468   
469   if (result) AliDebug(1, Form("Muon comes from a charmed hadron %d", GetParentPDGCode(0)));
470   
471   return result; 
472   
473 }
474
475 //====================================================================================================================================================
476
477 Bool_t AliMuonForwardTrack::IsFromBeauty() {
478
479   Bool_t result = kFALSE;
480
481   if ( GetParentPDGCode(0) ==   511 ||
482        GetParentPDGCode(0) ==   521 ||
483        GetParentPDGCode(0) == 10511 ||
484        GetParentPDGCode(0) == 10521 ||
485        GetParentPDGCode(0) ==   513 ||
486        GetParentPDGCode(0) ==   523 ||
487        GetParentPDGCode(0) == 10513 ||
488        GetParentPDGCode(0) == 10523 ||
489        GetParentPDGCode(0) == 20513 ||
490        GetParentPDGCode(0) == 20523 ||
491        GetParentPDGCode(0) ==   515 ||
492        GetParentPDGCode(0) ==   525 ||
493        GetParentPDGCode(0) ==   531 ||
494        GetParentPDGCode(0) == 10531 ||
495        GetParentPDGCode(0) ==   533 ||
496        GetParentPDGCode(0) == 10533 ||
497        GetParentPDGCode(0) == 20533 ||
498        GetParentPDGCode(0) ==   535 ||
499        GetParentPDGCode(0) ==   541 ||
500        GetParentPDGCode(0) == 10541 ||
501        GetParentPDGCode(0) ==   543 ||
502        GetParentPDGCode(0) == 10543 ||
503        GetParentPDGCode(0) == 20543 ||
504        GetParentPDGCode(0) ==   545 ) result = kTRUE;
505   
506   if (result) AliDebug(1, Form("Muon comes from a beauty hadron %d", GetParentPDGCode(0)));
507   
508   return result; 
509   
510 }
511
512 //====================================================================================================================================================
513
514 Bool_t AliMuonForwardTrack::IsFromBackground() {
515
516   Bool_t result = kFALSE;
517
518   if (!IsFromResonance() && !IsFromCharm() && !IsFromBeauty()) result = kTRUE;
519
520   if (result) AliDebug(1, Form("Muon comes from a background source %d", GetParentPDGCode(0)));
521
522   return result;
523
524 }
525
526 //====================================================================================================================================================