MFT module
[u/mrichter/AliRoot.git] / MFT / MFTbase / 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"
7bafb008 34#include "TDecompLU.h"
bffc7f8c 35
36#include "AliMFTAnalysisTools.h"
37
38ClassImp(AliMFTAnalysisTools)
39
40//====================================================================================================================================================
41
42Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2]) {
43
2666aca2 44 if (!(muon->PzAtDCA()!=0)) return kFALSE;
bffc7f8c 45
46 AliMUONTrackParam *param = new AliMUONTrackParam();
47
48 param -> SetNonBendingCoor(muon->XAtDCA());
49 param -> SetBendingCoor(muon->YAtDCA());
50 param -> SetZ(AliMFTConstants::fZEvalKinem);
2666aca2 51 param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
52 param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
53 param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon->PyAtDCA()/muon->PzAtDCA(),2))) );
bffc7f8c 54
55 AliMUONTrackExtrap::ExtrapToZ(param, z);
56 xy[0] = param->GetNonBendingCoor();
57 xy[1] = param->GetBendingCoor();
58
59 delete param;
60
61 return kTRUE;
62
63}
64
65//====================================================================================================================================================
66
67Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2], TLorentzVector &kinem) {
68
2666aca2 69 if (!(muon->PzAtDCA()!=0)) return kFALSE;
bffc7f8c 70
71 AliMUONTrackParam *param = new AliMUONTrackParam();
72
73 param -> SetNonBendingCoor(muon->XAtDCA());
74 param -> SetBendingCoor(muon->YAtDCA());
75 param -> SetZ(AliMFTConstants::fZEvalKinem);
2666aca2 76 param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
77 param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
78 param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon->PyAtDCA()/muon->PzAtDCA(),2))) );
bffc7f8c 79
80 AliMUONTrackExtrap::ExtrapToZ(param, z);
81 xy[0] = param->GetNonBendingCoor();
82 xy[1] = param->GetBendingCoor();
83
84 Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
85 Double_t energy = TMath::Sqrt(massMu*massMu + param->Px()*param->Px() + param->Py()*param->Py() + param->Pz()*param->Pz());
86
87 kinem.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
88
89 delete param;
90
91 return kTRUE;
92
93}
94
95//====================================================================================================================================================
96
97Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2], TLorentzVector &kinem, TMatrixD &cov) {
98
9664859a 99 // Extrapolate muon to a given z providing the corresponding (x,y) position and updating kinematics and covariance matrix
100
2666aca2 101 if (!(muon->PzAtDCA()!=0)) return kFALSE;
bffc7f8c 102
103 AliMUONTrackParam *param = new AliMUONTrackParam();
104
105 param -> SetNonBendingCoor(muon->XAtDCA());
106 param -> SetBendingCoor(muon->YAtDCA());
107 param -> SetZ(AliMFTConstants::fZEvalKinem);
2666aca2 108 param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
109 param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
110 param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon->PyAtDCA()/muon->PzAtDCA(),2))) );
bffc7f8c 111
112 param -> SetCovariances(ConvertCovMatrixAOD2MUON(muon));
113
114 AliMUONTrackExtrap::ExtrapToZCov(param, z);
115 xy[0] = param->GetNonBendingCoor();
116 xy[1] = param->GetBendingCoor();
117
118 Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
119 Double_t energy = TMath::Sqrt(massMu*massMu + param->Px()*param->Px() + param->Py()*param->Py() + param->Pz()*param->Pz());
120
121 kinem.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
122
123 cov = param->GetCovariances();
124
125 delete param;
126
127 return kTRUE;
128
129}
130
131//====================================================================================================================================================
132
9664859a 133Bool_t AliMFTAnalysisTools::ExtrapAODMuonToXY(AliAODTrack *muon, Double_t xy[2], Double_t &zFinal, TLorentzVector &kinem, TMatrixD &cov) {
134
135 // Find the point of closest approach between the muon and the direction parallel to the z-axis defined by the given (x,y)
136 // Provide the z of the above point as weel as the updated kinematics and covariance matrix
137
138 // We look for the above-defined PCA
139
9664859a 140 AliMUONTrackParam *param = new AliMUONTrackParam();
141 param -> SetNonBendingCoor(muon->XAtDCA());
142 param -> SetBendingCoor(muon->YAtDCA());
143 param -> SetZ(AliMFTConstants::fZEvalKinem);
2666aca2 144 param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
145 param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
146 param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon->PyAtDCA()/muon->PzAtDCA(),2))) );
9664859a 147
148 // here we want to understand in which direction we have to search the minimum...
149
150 Double_t step = 1.; // initial step, in cm
151 Double_t startPoint = 0.;
152
2666aca2 153 Double_t r[3]={0}, z[3]={startPoint-step, startPoint, startPoint+step};
9664859a 154
155 TVector3 **points = new TVector3*[2]; // points[0] for the muon, points[1] for the direction parallel to the z-axis defined by the given (x,y)
156
157 for (Int_t i=0; i<3; i++) {
158 AliMUONTrackExtrap::ExtrapToZ(param, z[i]);
159 points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]);
160 points[1] = new TVector3(xy[0],xy[1],z[i]);
161 r[i] = GetDistanceBetweenPoints(points,2);
162 for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu];
163 }
164
165 Int_t researchDirection = 0;
166
167 if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
168 else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative
169 else if (r[0]<r[1] && r[1]>r[2]) {
170 printf("E-AliMFTAnalysisTools::ExtrapAODMuonToXY: Point of closest approach cannot be found (no minima)\n");
171 delete param;
172 delete points;
173 return kFALSE;
174 }
175
176 while (TMath::Abs(researchDirection)>0.5) {
177
178 if (researchDirection>0) {
179 z[0] = z[1];
180 z[1] = z[2];
181 z[2] = z[1]+researchDirection*step;
182 }
183 else {
184 z[2] = z[1];
185 z[1] = z[0];
186 z[0] = z[1]+researchDirection*step;
187 }
188 if (TMath::Abs(z[0])>900.) {
189 printf("E-AliMFTAnalysisTools::ExtrapAODMuonToXY: Point of closest approach cannot be found (no minima in the fiducial region)\n");
190 delete param;
191 delete points;
192 return kFALSE;
193 }
194
195 for (Int_t i=0; i<3; i++) {
196 AliMUONTrackExtrap::ExtrapToZ(param, z[i]);
197 points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]);
198 points[1] = new TVector3(xy[0],xy[1],z[i]);
199 r[i] = GetDistanceBetweenPoints(points,2);
200 for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu];
201 }
202
203 researchDirection=0;
204 if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
205 else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative
206
207 }
208
209 // now we now that the minimum is between z[0] and z[2] and we search for it
210
211 step *= 0.5;
212 while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) {
213 z[0] = z[1]-step;
214 z[2] = z[1]+step;
215 for (Int_t i=0; i<3; i++) {
216 AliMUONTrackExtrap::ExtrapToZ(param, z[i]);
217 points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]);
218 points[1] = new TVector3(xy[0],xy[1],z[i]);
219 r[i] = GetDistanceBetweenPoints(points,2);
220 for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu];
221 }
222 if (r[0]<r[1]) z[1] = z[0];
223 else if (r[2]<r[1]) z[1] = z[2];
224 else step *= 0.5;
225 }
226
227 zFinal = z[1];
228
229 Double_t xyMuon[2] = {0};
230 ExtrapAODMuonToZ(muon, zFinal, xyMuon, kinem, cov);
231
232 return kTRUE;
233
234}
235
236//====================================================================================================================================================
237
a8b6b8e9 238Bool_t AliMFTAnalysisTools::GetAODMuonOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
bffc7f8c 239
240 Double_t xy[2] = {0};
241 ExtrapAODMuonToZ(muon, zv, xy);
242
a8b6b8e9 243 offset = TMath::Sqrt((xv-xy[0])*(xv-xy[0]) + (yv-xy[1])*(yv-xy[1]));
244
245 return kTRUE;
bffc7f8c 246
247}
248
249//====================================================================================================================================================
250
9664859a 251Bool_t AliMFTAnalysisTools::GetAODMuonOffsetZ(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
252
253 Double_t xy[2] = {xv, yv};
254 Double_t zFinal = 0;
255 TLorentzVector kinem(0,0,0,0);
256 TMatrixD cov(5,5);
257
258 ExtrapAODMuonToXY(muon, xy, zFinal, kinem, cov);
259
260 offset = TMath::Abs(zFinal - zv);
261
262 return kTRUE;
263
264}
265
266//====================================================================================================================================================
267
a8b6b8e9 268Bool_t AliMFTAnalysisTools::GetAODMuonWeightedOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
bffc7f8c 269
270 Double_t xy[2] = {0};
271 TLorentzVector kinem(0,0,0,0);
272 TMatrixD cov(5,5);
273
274 ExtrapAODMuonToZ(muon, zv, xy, kinem, cov);
275
276 TMatrixD covCoordinates(2,2);
277 covCoordinates(0,0) = cov(0,0);
278 covCoordinates(0,1) = cov(0,2);
279 covCoordinates(1,0) = cov(2,0);
280 covCoordinates(1,1) = cov(2,2);
281
2666aca2 282 if (covCoordinates.Determinant() < covCoordinates.GetTol()) return kFALSE;
283
7bafb008 284 if (TDecompLU::InvertLU(covCoordinates,covCoordinates.GetTol(),0)) {
285
286 TMatrixD covCoordinatesInverse = covCoordinates;
287 Double_t dX = xy[0] - xv;
288 Double_t dY = xy[1] - yv;
289
a8b6b8e9 290 offset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) +
291 dY*dY*covCoordinatesInverse(1,1) +
292 2.*dX*dY*covCoordinatesInverse(0,1)));
7bafb008 293
a8b6b8e9 294 return kTRUE;
7bafb008 295
296 }
bffc7f8c 297
a8b6b8e9 298 return kFALSE;
bffc7f8c 299
300}
301
302//====================================================================================================================================================
303
304Double_t AliMFTAnalysisTools::GetPseudoProperDecayTimeXY(Double_t xVtx, Double_t yVtx,
305 Double_t xDimu, Double_t yDimu,
306 Double_t mDimu, Double_t ptDimu) {
307
308 // pseudo-proper decay time of a particle produced in the primary vertex and decaying into a dimuon (+ X)
309 // evaluated using the transverse degree of freedom of the decay topology
310
311 if (ptDimu != 0) {
312 Double_t decayLengthXY = TMath::Sqrt((xVtx-xDimu)*(xVtx-xDimu)+(yVtx-yDimu)*(yVtx-yDimu));
313 return (decayLengthXY * mDimu/ptDimu)/TMath::Ccgs()*1E12; // in ps
314 }
315
316 return -99999999;
317
318}
319
320//====================================================================================================================================================
321
322Double_t AliMFTAnalysisTools::GetPseudoProperDecayTimeZ(Double_t zVtx,
323 Double_t zDimu,
324 Double_t mDimu, Double_t pzDimu) {
325
326 // pseudo-proper decay time of a particle produced in the primary vertex and decaying into a dimuon (+ X)
327 // evaluated using the longitudinal degree of freedom of the decay topology
328
329 if (pzDimu != 0) {
330 Double_t decayLengthZ = zDimu - zVtx;
331 return (decayLengthZ * mDimu/pzDimu)/TMath::Ccgs()*1E12; // in ps
332 }
333
334 return -99999999;
335
336}
337
338//====================================================================================================================================================
339
19393921 340Bool_t AliMFTAnalysisTools::CalculatePCA(AliAODDimuon *dimuon, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
341
342 TObjArray *muons = new TObjArray();
343 muons -> Add(dimuon->GetMu(0));
344 muons -> Add(dimuon->GetMu(1));
07d2587b 345
346 Bool_t result = CalculatePCA(muons, pca, pcaQuality, kinem);
347 delete muons;
348 return result;
19393921 349
350}
351
352//====================================================================================================================================================
353
354Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
bffc7f8c 355
356 const Int_t nMuons = muons->GetEntriesFast();
357 if (nMuons<2 || nMuons>AliMFTConstants::fNMaxMuonsForPCA) {
358 printf("W-AliMFTAnalysisTools::CalculatePCA: number of muons not valid\n");
359 return kFALSE;
360 }
2666aca2 361
bffc7f8c 362 Double_t fXPointOfClosestApproach=0, fYPointOfClosestApproach=0, fZPointOfClosestApproach=0;
07d2587b 363
6ffd24bb 364 AliAODTrack *muon[AliMFTConstants::fNMaxMuonsForPCA] = {0};
365 AliMUONTrackParam *param[AliMFTConstants::fNMaxMuonsForPCA] = {0};
bffc7f8c 366
367 // Finding AliMUONTrackParam objects for each muon
368
369 for (Int_t iMu=0; iMu<nMuons; iMu++) {
370 muon[iMu] = (AliAODTrack*) muons->At(iMu);
2666aca2 371 if (TMath::Abs(muon[iMu]->PzAtDCA())<1.e-6) {
07d2587b 372 for(Int_t i=0;i<iMu;i++) delete param[i];
373 return kFALSE;
374 }
bffc7f8c 375 param[iMu] = new AliMUONTrackParam();
376 param[iMu] -> SetNonBendingCoor(muon[iMu]->XAtDCA());
377 param[iMu] -> SetBendingCoor(muon[iMu]->YAtDCA());
9664859a 378 param[iMu] -> SetZ(AliMFTConstants::fZEvalKinem);
2666aca2 379 param[iMu] -> SetNonBendingSlope(muon[iMu]->PxAtDCA()/muon[iMu]->PzAtDCA());
380 param[iMu] -> SetBendingSlope(muon[iMu]->PyAtDCA()/muon[iMu]->PzAtDCA());
381 param[iMu] -> SetInverseBendingMomentum( muon[iMu]->Charge() * (1./muon[iMu]->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon[iMu]->PyAtDCA()/muon[iMu]->PzAtDCA(),2))) );
bffc7f8c 382 }
07d2587b 383
bffc7f8c 384 // here we want to understand in which direction we have to search the minimum...
385
386 Double_t step = 1.; // initial step, in cm
387 Double_t startPoint = 0.;
07d2587b 388
2666aca2 389 Double_t r[3]={0}, z[3]={startPoint-step, startPoint, startPoint+step};
bffc7f8c 390
6ffd24bb 391 TVector3 **points = new TVector3*[AliMFTConstants::fNMaxMuonsForPCA];
bffc7f8c 392
393 for (Int_t i=0; i<3; i++) {
394 for (Int_t iMu=0; iMu<nMuons; iMu++) {
2666aca2 395 // if (TMath::Abs(param[iMu]->GetInverseBendingMomentum())<1.) {
396 // printf("W-AliMFTAnalysisTools::CalculatePCA: Evoiding floating point exception in PCA finding\n");
397 // return kFALSE;
398 // }
bffc7f8c 399 AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
400 points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
2666aca2 401 }
bffc7f8c 402 r[i] = GetDistanceBetweenPoints(points,nMuons);
07d2587b 403 for (Int_t iMu=0; iMu<nMuons; iMu++) delete points[iMu];
bffc7f8c 404 }
07d2587b 405
bffc7f8c 406 Int_t researchDirection = 0;
407
408 if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
409 else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative
410 else if (r[0]<r[1] && r[1]>r[2]) {
411 printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima)\n");
07d2587b 412 for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
413 delete points;
bffc7f8c 414 return kFALSE;
07d2587b 415 }
416
bffc7f8c 417 while (TMath::Abs(researchDirection)>0.5) {
07d2587b 418
bffc7f8c 419 if (researchDirection>0) {
420 z[0] = z[1];
421 z[1] = z[2];
422 z[2] = z[1]+researchDirection*step;
423 }
424 else {
425 z[2] = z[1];
426 z[1] = z[0];
427 z[0] = z[1]+researchDirection*step;
428 }
429 if (TMath::Abs(z[0])>900.) {
430 printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima in the fiducial region)\n");
07d2587b 431 for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
432 delete points;
bffc7f8c 433 return kFALSE;
434 }
435
bffc7f8c 436 for (Int_t i=0; i<3; i++) {
437 for (Int_t iMu=0; iMu<nMuons; iMu++) {
07d2587b 438 AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
439 points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
bffc7f8c 440 }
441 r[i] = GetDistanceBetweenPoints(points,nMuons);
07d2587b 442 for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu];
bffc7f8c 443 }
444 researchDirection=0;
445 if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
446 else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative
07d2587b 447
bffc7f8c 448 }
07d2587b 449
450 // now we now that the minimum is between z[0] and z[2] and we search for it
bffc7f8c 451
2666aca2 452 Int_t nSteps = 0;
453
bffc7f8c 454 step *= 0.5;
455 while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) {
456 z[0] = z[1]-step;
457 z[2] = z[1]+step;
458 for (Int_t i=0; i<3; i++) {
459 for (Int_t iMu=0; iMu<nMuons; iMu++) {
07d2587b 460 AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
461 points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
bffc7f8c 462 }
463 r[i] = GetDistanceBetweenPoints(points,nMuons);
2666aca2 464 for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu];
bffc7f8c 465 }
2666aca2 466 // printf("Step #%d : %f %f %f\n",nSteps,r[0],r[1],r[2]);
bffc7f8c 467 if (r[0]<r[1]) z[1] = z[0];
468 else if (r[2]<r[1]) z[1] = z[2];
469 else step *= 0.5;
2666aca2 470 nSteps++;
bffc7f8c 471 }
2666aca2 472
473 // if (TMath::Abs(z[1]-1.)<0.1) {
474 // printf("Minimum found in %f in %d steps. Step = %f. p1->X() = %f, p2->X() = %f, p1->Y() = %f, p2->Y() = %f\n",
475 // z[1], nSteps, step, points[0]->X(), points[1]->X(), points[0]->Y(), points[1]->Y());
476 // printf("m1->X() = %f, m2->X() = %f, m1->Y() = %f, m2->Y() = %f\n",
477 // muon[0]->XAtDCA(),muon[1]->XAtDCA(), muon[0]->YAtDCA(),muon[1]->YAtDCA());
478 // }
479
bffc7f8c 480 // Once z of minimum is found, we evaluate the x and y coordinates by averaging over the contributing tracks
481
482 fZPointOfClosestApproach = z[1];
483 fXPointOfClosestApproach = 0.;
484 fYPointOfClosestApproach = 0.;
485 for (Int_t iMu=0; iMu<nMuons; iMu++) {
486 AliMUONTrackExtrap::ExtrapToZ(param[iMu], fZPointOfClosestApproach);
487 fXPointOfClosestApproach += param[iMu]->GetNonBendingCoor();
488 fYPointOfClosestApproach += param[iMu]->GetBendingCoor();
489 }
490 fXPointOfClosestApproach /= Double_t(nMuons);
491 fYPointOfClosestApproach /= Double_t(nMuons);
492
493 pca[0] = fXPointOfClosestApproach;
494 pca[1] = fYPointOfClosestApproach;
495 pca[2] = fZPointOfClosestApproach;
496
497 // Evaluating the kinematics of the N-muon
07d2587b 498
bffc7f8c 499 Double_t pTot[3] = {0};
500 Double_t ene = 0.;
501 Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
502 for (Int_t iMu=0; iMu<nMuons; iMu++) {
503 pTot[0] += param[iMu]->Px();
504 pTot[1] += param[iMu]->Py();
505 pTot[2] += param[iMu]->Pz();
506 ene += TMath::Sqrt(massMu*massMu + param[iMu]->Px()*param[iMu]->Px() + param[iMu]->Py()*param[iMu]->Py() + param[iMu]->Pz()*param[iMu]->Pz());
507 }
508
509 kinem.SetPxPyPzE(pTot[0], pTot[1], pTot[2], ene);
07d2587b 510
bffc7f8c 511 // Evaluating the PCA quality of the N-muon
07d2587b 512
bffc7f8c 513 Double_t sum=0.,squareSum=0.;
514 for (Int_t iMu=0; iMu<nMuons; iMu++) {
a8b6b8e9 515 Double_t wOffset = 0;
516 if (!GetAODMuonWeightedOffset(muon[iMu],fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach, wOffset)) {
517 for(Int_t jMu=0;jMu<nMuons;jMu++) delete param[jMu];
518 delete points;
519 return kFALSE;
520 }
bffc7f8c 521 Double_t f = TMath::Exp(-0.5 * wOffset);
522 sum += f;
523 squareSum += f*f;
524 }
525 if (sum > 0.) pcaQuality = (sum-squareSum/sum) / (nMuons-1);
526 else pcaQuality = 0.;
07d2587b 527
528 for(Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
529 delete points;
bffc7f8c 530 return kTRUE;
531
532}
533
534//=========================================================================================================================
535
536Double_t AliMFTAnalysisTools::GetDistanceBetweenPoints(TVector3 **points, Int_t nPoints) {
537
538 if (nPoints>AliMFTConstants::fNMaxMuonsForPCA) {
539 printf("W-AliMFTAnalysisTools::GetDistanceBetweenPoints: number of points not valid\n");
540 return 1.e9;
541 }
542
543 if (nPoints<2) return 0.;
544 if (nPoints<3) return TMath::Sqrt( (points[0]->X()-points[1]->X()) * (points[0]->X()-points[1]->X()) +
2666aca2 545 (points[0]->Y()-points[1]->Y()) * (points[0]->Y()-points[1]->Y()) );
546 // (points[0]->Z()-points[1]->Z()) * (points[0]->Z()-points[1]->Z()) );
bffc7f8c 547
548 const Int_t nEdgesMax = ((AliMFTConstants::fNMaxMuonsForPCA) * (AliMFTConstants::fNMaxMuonsForPCA - 1)) / 2;
549
550 Int_t startID[nEdgesMax] = {0};
551 Int_t stopID[nEdgesMax] = {0};
552 Double_t edgeLength[nEdgesMax] = {0};
553
554 Bool_t pointStatus[AliMFTConstants::fNMaxMuonsForPCA] = {0};
555
556 Int_t nEdges=0;
557 for (Int_t i=0; i<nPoints-1; i++) {
558 for (Int_t j=i+1; j<nPoints; j++) {
559 edgeLength[nEdges] = TMath::Sqrt( (points[i]->X()-points[j]->X()) * (points[i]->X()-points[j]->X()) +
560 (points[i]->Y()-points[j]->Y()) * (points[i]->Y()-points[j]->Y()) +
561 (points[i]->Z()-points[j]->Z()) * (points[i]->Z()-points[j]->Z()) );
562 stopID[nEdges] = i;
563 startID[nEdges] = j;
564 nEdges++;
565 }
566 }
567
568 // Order Edges
569
570 Double_t min = 0;
571 Int_t iMin = 0;
572
573 for (Int_t iEdge=0; iEdge<nEdges-1; iEdge++) {
574 min = edgeLength[iEdge];
575 iMin = iEdge;
576 for (Int_t j=iEdge+1; j<nEdges; j++) {
577 if (edgeLength[j]<min) {
578 min = edgeLength[j];
579 iMin = j;
580 }
581 }
582
583 if (iMin != iEdge) {
584
585 Double_t edgeLengthMin = edgeLength[iMin];
586 Int_t startIDmin = startID[iMin];
587 Int_t stopIDmin = stopID[iMin];
588
589 edgeLength[iMin] = edgeLength[iEdge];
590 startID[iMin] = startID[iEdge];
591 stopID[iMin] = stopID[iEdge];
592
593 edgeLength[iEdge] = edgeLengthMin;
594 startID[iEdge] = startIDmin;
595 stopID[iEdge] = stopIDmin;
596
597 }
598
599 }
600
601 // Connect
602
603 Double_t length = 0.;
604 for (Int_t i=0; i<nEdges; i++) {
605 if (!(pointStatus[startID[i]] && pointStatus[stopID[i]])) {
606 pointStatus[startID[i]] = kTRUE;
607 pointStatus[stopID[i]] = kTRUE;
608 length += edgeLength[i];
609 }
610 }
611
612 return length;
613
614}
615
616//====================================================================================================================================================
617
618void AliMFTAnalysisTools::ConvertCovMatrixMUON2AOD(const TMatrixD& covMUON, Double_t covAOD[21]) {
619
620 // Converts the cov matrix from the MUON format (TMatrixD) to the AOD one (Double_t[21])
621 //
622 // Cov(x,x) ... : cv[0]
623 // Cov(x,slopeX) ... : cv[1] cv[2]
624 // Cov(x,y) ... : cv[3] cv[4] cv[5]
625 // Cov(x,slopeY) ... : cv[6] cv[7] cv[8] cv[9]
626 // Cov(x,invP_yz) ... : cv[10] cv[11] cv[12] cv[13] cv[14]
627 // not-used ... : cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
628
629 covAOD[0] = covMUON(0,0);
630
631 covAOD[1] = covMUON(1,0);
632 covAOD[2] = covMUON(1,1);
633
634 covAOD[3] = covMUON(2,0);
635 covAOD[4] = covMUON(2,1);
636 covAOD[5] = covMUON(2,2);
637
638 covAOD[6] = covMUON(3,0);
639 covAOD[7] = covMUON(3,1);
640 covAOD[8] = covMUON(3,2);
641 covAOD[9] = covMUON(3,3);
642
643 covAOD[10] = covMUON(4,0);
644 covAOD[11] = covMUON(4,1);
645 covAOD[12] = covMUON(4,2);
646 covAOD[13] = covMUON(4,3);
647 covAOD[14] = covMUON(4,4);
648
649 covAOD[15] = 0;
650 covAOD[16] = 0;
651 covAOD[17] = 0;
652 covAOD[18] = 0;
653 covAOD[19] = 0;
654 covAOD[20] = 0;
655
656}
657
658//====================================================================================================================================================
659
660const TMatrixD AliMFTAnalysisTools::ConvertCovMatrixAOD2MUON(AliAODTrack *muon) {
661
662 Double_t covAOD[21] = {0};
663 muon -> GetCovarianceXYZPxPyPz(covAOD);
664
665 TMatrixD covMUON(5,5);
666
667 covMUON(0,0) = covAOD[0];
668
669 covMUON(1,0) = covAOD[1];
670 covMUON(1,1) = covAOD[2];
671
672 covMUON(2,0) = covAOD[3];
673 covMUON(2,1) = covAOD[4];
674 covMUON(2,2) = covAOD[5];
675
676 covMUON(3,0) = covAOD[6];
677 covMUON(3,1) = covAOD[7];
678 covMUON(3,2) = covAOD[8];
679 covMUON(3,3) = covAOD[9];
680
681 covMUON(4,0) = covAOD[10];
682 covMUON(4,1) = covAOD[11];
683 covMUON(4,2) = covAOD[12];
684 covMUON(4,3) = covAOD[13];
685 covMUON(4,4) = covAOD[14];
686
687 return covMUON;
688
689}
690
691//====================================================================================================================================================
692
0ebb458a 693Bool_t AliMFTAnalysisTools::IsCorrectMatch(AliAODTrack *muon) {
694
695 for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) if (IsWrongCluster(muon, iPlane)) return kFALSE;
696 return kTRUE;
697
698}
699
700//====================================================================================================================================================
701
702TString AliMFTAnalysisTools::GetGenerator(Int_t label, AliAODMCHeader* header) {
703
704 // get the name of the generator that produced a given particle
705
706 Int_t partCounter = 0;
707 TList *genHeaders = header->GetCocktailHeaders();
708 Int_t nGenHeaders = genHeaders->GetEntries();
709
710 for (Int_t i=0; i<nGenHeaders; i++){
711 AliGenEventHeader *gh = (AliGenEventHeader*) genHeaders->At(i);
712 TString genName = gh->GetName();
713 Int_t nPart = gh->NProduced();
714 if (label>=partCounter && label<(partCounter+nPart)) return genName;
715 partCounter += nPart;
716 }
717
718 TString empty="";
719 return empty;
720
721}
722
723//====================================================================================================================================================
724
725void AliMFTAnalysisTools::GetTrackPrimaryGenerator(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC, TString &nameGen) {
726
727 // method to check if a track comes from a given generator
728
729 Int_t label = TMath::Abs(track->GetLabel());
730 nameGen = GetGenerator(label,header);
731
732 // In case the particle is not primary nameGen will contain blank spaces. In this case, we search backward for the primary which originated the chain
733
734 while (nameGen.IsWhitespace()) {
735 AliAODMCParticle *mcPart = (AliAODMCParticle*) arrayMC->At(label);
736 if (!mcPart) {
737 printf("AliMFTAnalysisTools::GetTrackPrimaryGenerator - BREAK: No valid AliAODMCParticle at label %i\n",label);
738 break;
739 }
740 Int_t motherLabel = mcPart->GetMother();
741 if (motherLabel < 0) {
742 printf("AliMFTAnalysisTools::GetTrackPrimaryGenerator - BREAK: Reached primary particle without valid mother\n");
743 break;
744 }
745 label = motherLabel;
746 nameGen = GetGenerator(label,header);
747 }
748
749 return;
750
751}
752
753//====================================================================================================================================================
754
2666aca2 755Bool_t AliMFTAnalysisTools::IsTrackInjected(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC) {
0ebb458a 756
757 // method to check if a track comes from the signal event or from the underlying Hijing event
758
759 TString nameGen;
760
761 GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
762
763 if (nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
764
765 return kTRUE;
2666aca2 766
0ebb458a 767}
768
769//====================================================================================================================================================
770
2666aca2 771Bool_t AliMFTAnalysisTools::TranslateMuon(AliAODTrack *muon, Double_t vtxInitial[3], Double_t vtxFinal[3]) {
772
773 if (!(muon->PzAtDCA()!=0)) return kFALSE;
774
775 AliMUONTrackParam *param = new AliMUONTrackParam();
776
777 Double_t deltaVtx[3] = {0};
778 for (Int_t i=0; i<3; i++) deltaVtx[i] = vtxInitial[i] - vtxFinal[i];
779
780 param -> SetNonBendingCoor(muon->XAtDCA());
781 param -> SetBendingCoor(muon->YAtDCA());
782 param -> SetZ(AliMFTConstants::fZEvalKinem);
783 param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
784 param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
785 param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon->PyAtDCA()/muon->PzAtDCA(),2))) );
786
787 // This will be interpreted as if the track (produced in an event having the primary vertex in vtxInitial)
788 // were produced in an event having the primary vertex in vtxFinal
789
790 AliMUONTrackExtrap::ExtrapToZ(param, deltaVtx[2]);
791 muon->SetXYAtDCA(param->GetNonBendingCoor() - deltaVtx[0], param->GetBendingCoor() - deltaVtx[1]);
792 muon->SetPxPyPzAtDCA(param->Px(), param->Py(), param->Pz());
793
794 delete param;
795
796 return kTRUE;
797
798}
799
800//====================================================================================================================================================
801
802Bool_t AliMFTAnalysisTools::TranslateMuonToOrigin(AliAODTrack *muon, Double_t vtx[3]) {
803
804 Double_t origin[3] = {0,0,0};
805
806 return TranslateMuon(muon, vtx, origin);
807
808}
809
810//====================================================================================================================================================
f67f98b6 811
812Bool_t AliMFTAnalysisTools::IsPDGCharm(Int_t pdgCode) {
813
814 pdgCode = TMath::Abs(pdgCode/100);
815 if (pdgCode>9) pdgCode /= 10;
816 if (pdgCode == 4 ) return kTRUE;
817 else return kFALSE;
818
819}
820
821//====================================================================================================================================================
822
823Bool_t AliMFTAnalysisTools::IsPDGBeauty(Int_t pdgCode) {
824
825 pdgCode = TMath::Abs(pdgCode/100);
826 if (pdgCode>9) pdgCode /= 10;
827 if (pdgCode == 5) return kTRUE;
828 else return kFALSE;
829
830}
831
832//====================================================================================================================================================
833
834Bool_t AliMFTAnalysisTools::IsPDGResonance(Int_t pdgCode) {
835
836 Int_t id = pdgCode%100000;
837 return (!((id-id%10)%110));
838
839}
840
841//====================================================================================================================================================