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