CA tracker - updates
[u/mrichter/AliRoot.git] / MFT / AliMFTAnalysisTools.cxx
CommitLineData
bffc7f8c 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// Support class for various common operation on MFT objects
19//
20// Contact author: antonio.uras@cern.ch
21//
22//====================================================================================================================================================
23
24#include "AliMUONTrackParam.h"
25#include "AliMUONTrackExtrap.h"
26#include "AliAODTrack.h"
27#include "AliAODDimuon.h"
28#include "TLorentzVector.h"
29#include "AliMFTConstants.h"
30#include "TDatabasePDG.h"
31#include "TMath.h"
32#include "AliLog.h"
19393921 33#include "TObjArray.h"
bffc7f8c 34
35#include "AliMFTAnalysisTools.h"
36
37ClassImp(AliMFTAnalysisTools)
38
39//====================================================================================================================================================
40
41Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2]) {
42
43 if (!(muon->Pz()!=0)) return kFALSE;
44
45 AliMUONTrackParam *param = new AliMUONTrackParam();
46
47 param -> SetNonBendingCoor(muon->XAtDCA());
48 param -> SetBendingCoor(muon->YAtDCA());
49 param -> SetZ(AliMFTConstants::fZEvalKinem);
50 param -> SetNonBendingSlope(muon->Px()/muon->Pz());
51 param -> SetBendingSlope(muon->Py()/muon->Pz());
52 param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) );
53
54 AliMUONTrackExtrap::ExtrapToZ(param, z);
55 xy[0] = param->GetNonBendingCoor();
56 xy[1] = param->GetBendingCoor();
57
58 delete param;
59
60 return kTRUE;
61
62}
63
64//====================================================================================================================================================
65
66Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2], TLorentzVector &kinem) {
67
68 if (!(muon->Pz()!=0)) return kFALSE;
69
70 AliMUONTrackParam *param = new AliMUONTrackParam();
71
72 param -> SetNonBendingCoor(muon->XAtDCA());
73 param -> SetBendingCoor(muon->YAtDCA());
74 param -> SetZ(AliMFTConstants::fZEvalKinem);
75 param -> SetNonBendingSlope(muon->Px()/muon->Pz());
76 param -> SetBendingSlope(muon->Py()/muon->Pz());
77 param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) );
78
79 AliMUONTrackExtrap::ExtrapToZ(param, z);
80 xy[0] = param->GetNonBendingCoor();
81 xy[1] = param->GetBendingCoor();
82
83 Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
84 Double_t energy = TMath::Sqrt(massMu*massMu + param->Px()*param->Px() + param->Py()*param->Py() + param->Pz()*param->Pz());
85
86 kinem.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
87
88 delete param;
89
90 return kTRUE;
91
92}
93
94//====================================================================================================================================================
95
96Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2], TLorentzVector &kinem, TMatrixD &cov) {
97
98 if (!(muon->Pz()!=0)) return kFALSE;
99
100 AliMUONTrackParam *param = new AliMUONTrackParam();
101
102 param -> SetNonBendingCoor(muon->XAtDCA());
103 param -> SetBendingCoor(muon->YAtDCA());
104 param -> SetZ(AliMFTConstants::fZEvalKinem);
105 param -> SetNonBendingSlope(muon->Px()/muon->Pz());
106 param -> SetBendingSlope(muon->Py()/muon->Pz());
107 param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) );
108
109 param -> SetCovariances(ConvertCovMatrixAOD2MUON(muon));
110
111 AliMUONTrackExtrap::ExtrapToZCov(param, z);
112 xy[0] = param->GetNonBendingCoor();
113 xy[1] = param->GetBendingCoor();
114
115 Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
116 Double_t energy = TMath::Sqrt(massMu*massMu + param->Px()*param->Px() + param->Py()*param->Py() + param->Pz()*param->Pz());
117
118 kinem.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
119
120 cov = param->GetCovariances();
121
122 delete param;
123
124 return kTRUE;
125
126}
127
128//====================================================================================================================================================
129
130Double_t AliMFTAnalysisTools::GetAODMuonOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv) {
131
132 Double_t xy[2] = {0};
133 ExtrapAODMuonToZ(muon, zv, xy);
134
135 return TMath::Sqrt((xv-xy[0])*(xv-xy[0]) + (yv-xy[1])*(yv-xy[1]));
136
137}
138
139//====================================================================================================================================================
140
141Double_t AliMFTAnalysisTools::GetAODMuonWeightedOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv) {
142
143 Double_t xy[2] = {0};
144 TLorentzVector kinem(0,0,0,0);
145 TMatrixD cov(5,5);
146
147 ExtrapAODMuonToZ(muon, zv, xy, kinem, cov);
148
149 TMatrixD covCoordinates(2,2);
150 covCoordinates(0,0) = cov(0,0);
151 covCoordinates(0,1) = cov(0,2);
152 covCoordinates(1,0) = cov(2,0);
153 covCoordinates(1,1) = cov(2,2);
154
155 TMatrixD covCoordinatesInverse = covCoordinates.Invert();
156
157 Double_t dX = xy[0] - xv;
158 Double_t dY = xy[1] - yv;
159
160 Double_t weightedOffset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) +
161 dY*dY*covCoordinatesInverse(1,1) +
162 2.*dX*dY*covCoordinatesInverse(0,1)));
163
164 return weightedOffset;
165
166}
167
168//====================================================================================================================================================
169
170Double_t AliMFTAnalysisTools::GetPseudoProperDecayTimeXY(Double_t xVtx, Double_t yVtx,
171 Double_t xDimu, Double_t yDimu,
172 Double_t mDimu, Double_t ptDimu) {
173
174 // pseudo-proper decay time of a particle produced in the primary vertex and decaying into a dimuon (+ X)
175 // evaluated using the transverse degree of freedom of the decay topology
176
177 if (ptDimu != 0) {
178 Double_t decayLengthXY = TMath::Sqrt((xVtx-xDimu)*(xVtx-xDimu)+(yVtx-yDimu)*(yVtx-yDimu));
179 return (decayLengthXY * mDimu/ptDimu)/TMath::Ccgs()*1E12; // in ps
180 }
181
182 return -99999999;
183
184}
185
186//====================================================================================================================================================
187
188Double_t AliMFTAnalysisTools::GetPseudoProperDecayTimeZ(Double_t zVtx,
189 Double_t zDimu,
190 Double_t mDimu, Double_t pzDimu) {
191
192 // pseudo-proper decay time of a particle produced in the primary vertex and decaying into a dimuon (+ X)
193 // evaluated using the longitudinal degree of freedom of the decay topology
194
195 if (pzDimu != 0) {
196 Double_t decayLengthZ = zDimu - zVtx;
197 return (decayLengthZ * mDimu/pzDimu)/TMath::Ccgs()*1E12; // in ps
198 }
199
200 return -99999999;
201
202}
203
204//====================================================================================================================================================
205
19393921 206Bool_t AliMFTAnalysisTools::CalculatePCA(AliAODDimuon *dimuon, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
207
208 TObjArray *muons = new TObjArray();
209 muons -> Add(dimuon->GetMu(0));
210 muons -> Add(dimuon->GetMu(1));
07d2587b 211
212 Bool_t result = CalculatePCA(muons, pca, pcaQuality, kinem);
213 delete muons;
214 return result;
19393921 215
216}
217
218//====================================================================================================================================================
219
220Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
bffc7f8c 221
222 const Int_t nMuons = muons->GetEntriesFast();
223 if (nMuons<2 || nMuons>AliMFTConstants::fNMaxMuonsForPCA) {
224 printf("W-AliMFTAnalysisTools::CalculatePCA: number of muons not valid\n");
225 return kFALSE;
226 }
227
228 Double_t fXPointOfClosestApproach=0, fYPointOfClosestApproach=0, fZPointOfClosestApproach=0;
07d2587b 229
6ffd24bb 230 AliAODTrack *muon[AliMFTConstants::fNMaxMuonsForPCA] = {0};
231 AliMUONTrackParam *param[AliMFTConstants::fNMaxMuonsForPCA] = {0};
bffc7f8c 232
233 // Finding AliMUONTrackParam objects for each muon
234
235 for (Int_t iMu=0; iMu<nMuons; iMu++) {
236 muon[iMu] = (AliAODTrack*) muons->At(iMu);
07d2587b 237 if (TMath::Abs(muon[iMu]->Pz())<1.e-6) {
238 for(Int_t i=0;i<iMu;i++) delete param[i];
239 return kFALSE;
240 }
bffc7f8c 241 param[iMu] = new AliMUONTrackParam();
242 param[iMu] -> SetNonBendingCoor(muon[iMu]->XAtDCA());
243 param[iMu] -> SetBendingCoor(muon[iMu]->YAtDCA());
244 param[iMu] -> SetZ(0.);
245 param[iMu] -> SetNonBendingSlope(muon[iMu]->Px()/muon[iMu]->Pz());
246 param[iMu] -> SetBendingSlope(muon[iMu]->Py()/muon[iMu]->Pz());
247 param[iMu] -> SetInverseBendingMomentum( muon[iMu]->Charge() * (1./muon[iMu]->Pz()) / (TMath::Sqrt(1+TMath::Power(muon[iMu]->Py()/muon[iMu]->Pz(),2))) );
248 }
07d2587b 249
bffc7f8c 250 // here we want to understand in which direction we have to search the minimum...
251
252 Double_t step = 1.; // initial step, in cm
253 Double_t startPoint = 0.;
07d2587b 254
bffc7f8c 255 Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step};
256
6ffd24bb 257 TVector3 **points = new TVector3*[AliMFTConstants::fNMaxMuonsForPCA];
bffc7f8c 258
259 for (Int_t i=0; i<3; i++) {
260 for (Int_t iMu=0; iMu<nMuons; iMu++) {
261 AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
262 points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
07d2587b 263 }
bffc7f8c 264 r[i] = GetDistanceBetweenPoints(points,nMuons);
07d2587b 265 for (Int_t iMu=0; iMu<nMuons; iMu++) delete points[iMu];
bffc7f8c 266 }
07d2587b 267
bffc7f8c 268 Int_t researchDirection = 0;
269
270 if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
271 else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative
272 else if (r[0]<r[1] && r[1]>r[2]) {
273 printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima)\n");
07d2587b 274 for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
275 delete points;
bffc7f8c 276 return kFALSE;
07d2587b 277 }
278
bffc7f8c 279 while (TMath::Abs(researchDirection)>0.5) {
07d2587b 280
bffc7f8c 281 if (researchDirection>0) {
282 z[0] = z[1];
283 z[1] = z[2];
284 z[2] = z[1]+researchDirection*step;
285 }
286 else {
287 z[2] = z[1];
288 z[1] = z[0];
289 z[0] = z[1]+researchDirection*step;
290 }
291 if (TMath::Abs(z[0])>900.) {
292 printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima in the fiducial region)\n");
07d2587b 293 for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
294 delete points;
bffc7f8c 295 return kFALSE;
296 }
297
bffc7f8c 298 for (Int_t i=0; i<3; i++) {
299 for (Int_t iMu=0; iMu<nMuons; iMu++) {
07d2587b 300 AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
301 points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
bffc7f8c 302 }
303 r[i] = GetDistanceBetweenPoints(points,nMuons);
07d2587b 304 for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu];
bffc7f8c 305 }
306 researchDirection=0;
307 if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
308 else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative
07d2587b 309
bffc7f8c 310 }
07d2587b 311
312 // now we now that the minimum is between z[0] and z[2] and we search for it
bffc7f8c 313
314 step *= 0.5;
315 while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) {
316 z[0] = z[1]-step;
317 z[2] = z[1]+step;
318 for (Int_t i=0; i<3; i++) {
319 for (Int_t iMu=0; iMu<nMuons; iMu++) {
07d2587b 320 AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
321 points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
bffc7f8c 322 }
323 r[i] = GetDistanceBetweenPoints(points,nMuons);
07d2587b 324 for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu];
bffc7f8c 325 }
326 if (r[0]<r[1]) z[1] = z[0];
327 else if (r[2]<r[1]) z[1] = z[2];
328 else step *= 0.5;
329 }
330
bffc7f8c 331 // Once z of minimum is found, we evaluate the x and y coordinates by averaging over the contributing tracks
332
333 fZPointOfClosestApproach = z[1];
334 fXPointOfClosestApproach = 0.;
335 fYPointOfClosestApproach = 0.;
336 for (Int_t iMu=0; iMu<nMuons; iMu++) {
337 AliMUONTrackExtrap::ExtrapToZ(param[iMu], fZPointOfClosestApproach);
338 fXPointOfClosestApproach += param[iMu]->GetNonBendingCoor();
339 fYPointOfClosestApproach += param[iMu]->GetBendingCoor();
340 }
341 fXPointOfClosestApproach /= Double_t(nMuons);
342 fYPointOfClosestApproach /= Double_t(nMuons);
343
344 pca[0] = fXPointOfClosestApproach;
345 pca[1] = fYPointOfClosestApproach;
346 pca[2] = fZPointOfClosestApproach;
347
348 // Evaluating the kinematics of the N-muon
07d2587b 349
bffc7f8c 350 Double_t pTot[3] = {0};
351 Double_t ene = 0.;
352 Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
353 for (Int_t iMu=0; iMu<nMuons; iMu++) {
354 pTot[0] += param[iMu]->Px();
355 pTot[1] += param[iMu]->Py();
356 pTot[2] += param[iMu]->Pz();
357 ene += TMath::Sqrt(massMu*massMu + param[iMu]->Px()*param[iMu]->Px() + param[iMu]->Py()*param[iMu]->Py() + param[iMu]->Pz()*param[iMu]->Pz());
358 }
359
360 kinem.SetPxPyPzE(pTot[0], pTot[1], pTot[2], ene);
07d2587b 361
bffc7f8c 362 // Evaluating the PCA quality of the N-muon
07d2587b 363
bffc7f8c 364 Double_t sum=0.,squareSum=0.;
365 for (Int_t iMu=0; iMu<nMuons; iMu++) {
366 Double_t wOffset = AliMFTAnalysisTools::GetAODMuonWeightedOffset(muon[iMu],fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach);
367 Double_t f = TMath::Exp(-0.5 * wOffset);
368 sum += f;
369 squareSum += f*f;
370 }
371 if (sum > 0.) pcaQuality = (sum-squareSum/sum) / (nMuons-1);
372 else pcaQuality = 0.;
07d2587b 373
374 for(Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
375 delete points;
bffc7f8c 376 return kTRUE;
377
378}
379
380//=========================================================================================================================
381
382Double_t AliMFTAnalysisTools::GetDistanceBetweenPoints(TVector3 **points, Int_t nPoints) {
383
384 if (nPoints>AliMFTConstants::fNMaxMuonsForPCA) {
385 printf("W-AliMFTAnalysisTools::GetDistanceBetweenPoints: number of points not valid\n");
386 return 1.e9;
387 }
388
389 if (nPoints<2) return 0.;
390 if (nPoints<3) return TMath::Sqrt( (points[0]->X()-points[1]->X()) * (points[0]->X()-points[1]->X()) +
391 (points[0]->Y()-points[1]->Y()) * (points[0]->Y()-points[1]->Y()) +
392 (points[0]->Z()-points[1]->Z()) * (points[0]->Z()-points[1]->Z()) );
393
394 const Int_t nEdgesMax = ((AliMFTConstants::fNMaxMuonsForPCA) * (AliMFTConstants::fNMaxMuonsForPCA - 1)) / 2;
395
396 Int_t startID[nEdgesMax] = {0};
397 Int_t stopID[nEdgesMax] = {0};
398 Double_t edgeLength[nEdgesMax] = {0};
399
400 Bool_t pointStatus[AliMFTConstants::fNMaxMuonsForPCA] = {0};
401
402 Int_t nEdges=0;
403 for (Int_t i=0; i<nPoints-1; i++) {
404 for (Int_t j=i+1; j<nPoints; j++) {
405 edgeLength[nEdges] = TMath::Sqrt( (points[i]->X()-points[j]->X()) * (points[i]->X()-points[j]->X()) +
406 (points[i]->Y()-points[j]->Y()) * (points[i]->Y()-points[j]->Y()) +
407 (points[i]->Z()-points[j]->Z()) * (points[i]->Z()-points[j]->Z()) );
408 stopID[nEdges] = i;
409 startID[nEdges] = j;
410 nEdges++;
411 }
412 }
413
414 // Order Edges
415
416 Double_t min = 0;
417 Int_t iMin = 0;
418
419 for (Int_t iEdge=0; iEdge<nEdges-1; iEdge++) {
420 min = edgeLength[iEdge];
421 iMin = iEdge;
422 for (Int_t j=iEdge+1; j<nEdges; j++) {
423 if (edgeLength[j]<min) {
424 min = edgeLength[j];
425 iMin = j;
426 }
427 }
428
429 if (iMin != iEdge) {
430
431 Double_t edgeLengthMin = edgeLength[iMin];
432 Int_t startIDmin = startID[iMin];
433 Int_t stopIDmin = stopID[iMin];
434
435 edgeLength[iMin] = edgeLength[iEdge];
436 startID[iMin] = startID[iEdge];
437 stopID[iMin] = stopID[iEdge];
438
439 edgeLength[iEdge] = edgeLengthMin;
440 startID[iEdge] = startIDmin;
441 stopID[iEdge] = stopIDmin;
442
443 }
444
445 }
446
447 // Connect
448
449 Double_t length = 0.;
450 for (Int_t i=0; i<nEdges; i++) {
451 if (!(pointStatus[startID[i]] && pointStatus[stopID[i]])) {
452 pointStatus[startID[i]] = kTRUE;
453 pointStatus[stopID[i]] = kTRUE;
454 length += edgeLength[i];
455 }
456 }
457
458 return length;
459
460}
461
462//====================================================================================================================================================
463
464void AliMFTAnalysisTools::ConvertCovMatrixMUON2AOD(const TMatrixD& covMUON, Double_t covAOD[21]) {
465
466 // Converts the cov matrix from the MUON format (TMatrixD) to the AOD one (Double_t[21])
467 //
468 // Cov(x,x) ... : cv[0]
469 // Cov(x,slopeX) ... : cv[1] cv[2]
470 // Cov(x,y) ... : cv[3] cv[4] cv[5]
471 // Cov(x,slopeY) ... : cv[6] cv[7] cv[8] cv[9]
472 // Cov(x,invP_yz) ... : cv[10] cv[11] cv[12] cv[13] cv[14]
473 // not-used ... : cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
474
475 covAOD[0] = covMUON(0,0);
476
477 covAOD[1] = covMUON(1,0);
478 covAOD[2] = covMUON(1,1);
479
480 covAOD[3] = covMUON(2,0);
481 covAOD[4] = covMUON(2,1);
482 covAOD[5] = covMUON(2,2);
483
484 covAOD[6] = covMUON(3,0);
485 covAOD[7] = covMUON(3,1);
486 covAOD[8] = covMUON(3,2);
487 covAOD[9] = covMUON(3,3);
488
489 covAOD[10] = covMUON(4,0);
490 covAOD[11] = covMUON(4,1);
491 covAOD[12] = covMUON(4,2);
492 covAOD[13] = covMUON(4,3);
493 covAOD[14] = covMUON(4,4);
494
495 covAOD[15] = 0;
496 covAOD[16] = 0;
497 covAOD[17] = 0;
498 covAOD[18] = 0;
499 covAOD[19] = 0;
500 covAOD[20] = 0;
501
502}
503
504//====================================================================================================================================================
505
506const TMatrixD AliMFTAnalysisTools::ConvertCovMatrixAOD2MUON(AliAODTrack *muon) {
507
508 Double_t covAOD[21] = {0};
509 muon -> GetCovarianceXYZPxPyPz(covAOD);
510
511 TMatrixD covMUON(5,5);
512
513 covMUON(0,0) = covAOD[0];
514
515 covMUON(1,0) = covAOD[1];
516 covMUON(1,1) = covAOD[2];
517
518 covMUON(2,0) = covAOD[3];
519 covMUON(2,1) = covAOD[4];
520 covMUON(2,2) = covAOD[5];
521
522 covMUON(3,0) = covAOD[6];
523 covMUON(3,1) = covAOD[7];
524 covMUON(3,2) = covAOD[8];
525 covMUON(3,3) = covAOD[9];
526
527 covMUON(4,0) = covAOD[10];
528 covMUON(4,1) = covAOD[11];
529 covMUON(4,2) = covAOD[12];
530 covMUON(4,3) = covAOD[13];
531 covMUON(4,4) = covAOD[14];
532
533 return covMUON;
534
535}
536
537//====================================================================================================================================================
538