]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MFT/AliMuonForwardTrack.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Methods added to distinguish muons from charm and beauty
[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrack.cxx
CommitLineData
820b4d9e 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"
d4643a10 34#include "AliMFTConstants.h"
820b4d9e 35
36ClassImp(AliMuonForwardTrack)
37
38//====================================================================================================================================================
39
40AliMuonForwardTrack::AliMuonForwardTrack():
41 AliMUONTrack(),
42 fMUONTrack(0),
43 fMCTrackRef(0),
d4643a10 44 fMFTClusters(0),
7e3dd1af 45 fNWrongClustersMC(-1),
46 fTrackMCId(-1)
820b4d9e 47{
48
49 // default constructor
d4643a10 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 }
820b4d9e 56 fMFTClusters = new TClonesArray("AliMFTCluster");
274c2dce 57 fMFTClusters -> SetOwner(kTRUE);
820b4d9e 58
59}
60
61//====================================================================================================================================================
62
63AliMuonForwardTrack::AliMuonForwardTrack(AliMUONTrack *MUONTrack):
64 AliMUONTrack(),
65 fMUONTrack(0),
66 fMCTrackRef(0),
d4643a10 67 fMFTClusters(0),
7e3dd1af 68 fNWrongClustersMC(-1),
69 fTrackMCId(-1)
820b4d9e 70{
71
72 SetMUONTrack(MUONTrack);
d4643a10 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 }
820b4d9e 78 fMFTClusters = new TClonesArray("AliMFTCluster");
274c2dce 79 fMFTClusters -> SetOwner(kTRUE);
820b4d9e 80
81}
82
83//====================================================================================================================================================
84
85AliMuonForwardTrack::AliMuonForwardTrack(const AliMuonForwardTrack& track):
86 AliMUONTrack(track),
d4643a10 87 fMUONTrack(0x0),
88 fMCTrackRef(0x0),
89 fMFTClusters(0x0),
7e3dd1af 90 fNWrongClustersMC(track.fNWrongClustersMC),
91 fTrackMCId(track.fTrackMCId)
820b4d9e 92{
93
94 // copy constructor
d4643a10 95 fMUONTrack = new AliMUONTrack(*(track.fMUONTrack));
96 if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
97 fMFTClusters = new TClonesArray(*(track.fMFTClusters));
274c2dce 98 fMFTClusters->SetOwner(kTRUE);
d4643a10 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 }
820b4d9e 104
105}
106
107//====================================================================================================================================================
108
109AliMuonForwardTrack& 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
274c2dce 120 Clear("");
820b4d9e 121
d4643a10 122 fMUONTrack = new AliMUONTrack(*(track.fMUONTrack));
123 if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
124 fMFTClusters = new TClonesArray(*(track.fMFTClusters));
274c2dce 125 fMFTClusters->SetOwner(kTRUE);
d4643a10 126 fNWrongClustersMC = track.fNWrongClustersMC;
7e3dd1af 127 fTrackMCId = track.fTrackMCId;
d4643a10 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 }
820b4d9e 134
135 return *this;
136
137}
138
139//====================================================================================================================================================
140
274c2dce 141void AliMuonForwardTrack::Clear(const Option_t* /*opt*/) {
142
143 // Clear arrays
3b0abc02 144 fMFTClusters -> Delete();
145 delete fMFTClusters; fMFTClusters = 0x0;
274c2dce 146 delete fMUONTrack; fMUONTrack = 0x0;
147 delete fMCTrackRef; fMCTrackRef = 0x0;
274c2dce 148
149}
150
151//====================================================================================================================================================
152
53b30119 153AliMuonForwardTrack::~AliMuonForwardTrack() {
154
155 delete fMUONTrack;
156 delete fMCTrackRef;
b2741f03 157 fMFTClusters -> Delete();
53b30119 158 delete fMFTClusters;
159
160}
161
162//====================================================================================================================================================
163
820b4d9e 164void AliMuonForwardTrack::SetMUONTrack(AliMUONTrack *MUONTrack) {
165
166 if (fMUONTrack) {
167 AliInfo("fMUONTrack already exists, nothing will be done");
168 return;
169 }
170
171 fMUONTrack = MUONTrack;
d4643a10 172 SetGlobalChi2(fMUONTrack->GetGlobalChi2());
820b4d9e 173
174}
175
176//====================================================================================================================================================
177
178void 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
191void AliMuonForwardTrack::AddTrackParamAtMFTCluster(AliMUONTrackParam &trackParam, AliMFTCluster &mftCluster) {
192
d4643a10 193 AliDebug(1, Form("Before adding: this->fMFTClusters=%p has %d entries", this->fMFTClusters, this->fMFTClusters->GetEntries()));
194 Int_t iMFTCluster = this->fMFTClusters->GetEntries();
820b4d9e 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);
820b4d9e 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()));
d4643a10 204 AliDebug(1, Form("Adding Kalman chi2 = %f to global chi2 = %f", chi2Kalman, GetGlobalChi2()));
820b4d9e 205 Double_t newGlobalChi2 = GetGlobalChi2() + chi2Kalman;
206 mftCluster.SetLocalChi2(chi2Kalman);
207 mftCluster.SetTrackChi2(newGlobalChi2);
d4643a10 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()));
820b4d9e 211 AliDebug(1, Form("muonCluster->GetZ() = %f, trackParam->GetZ() = %f",muonCluster->GetZ(), trackParam.GetZ()));
212 SetGlobalChi2(newGlobalChi2);
d4643a10 213 AliDebug(1, Form("New global chi2 = %f", GetGlobalChi2()));
820b4d9e 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
224AliMUONTrackParam* 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
239AliMUONTrackParam* 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
254AliMUONVCluster* 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
270AliMFTCluster* AliMuonForwardTrack::GetMFTCluster(Int_t iMFTCluster) {
271
272 if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
d4643a10 273 AliError(Form("Invalid MFT cluster index (%d). GetNMFTClusters()=%d. NULL pointer will be returned", iMFTCluster, GetNMFTClusters()));
820b4d9e 274 return NULL;
275 }
276
277 AliMUONTrackParam *trackParam = GetTrackParamAtMFTCluster(iMFTCluster);
d4643a10 278 AliMFTCluster *mftCluster = (AliMFTCluster*) this->fMFTClusters->At(trackParam->GetUniqueID());
820b4d9e 279
280 return mftCluster;
281
282}
283
284//====================================================================================================================================================
285
286Double_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
361Double_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
390Double_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
401Double_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
412Double_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
3a775f35 424Bool_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
446Bool_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
477Bool_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
514Bool_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//====================================================================================================================================================