]>
Commit | Line | Data |
---|---|---|
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 | ||
38 | ClassImp(AliMFTAnalysisTools) | |
39 | ||
40 | //==================================================================================================================================================== | |
41 | ||
42 | Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2]) { | |
43 | ||
44 | if (!(muon->Pz()!=0)) return kFALSE; | |
45 | ||
46 | AliMUONTrackParam *param = new AliMUONTrackParam(); | |
47 | ||
48 | param -> SetNonBendingCoor(muon->XAtDCA()); | |
49 | param -> SetBendingCoor(muon->YAtDCA()); | |
50 | param -> SetZ(AliMFTConstants::fZEvalKinem); | |
51 | param -> SetNonBendingSlope(muon->Px()/muon->Pz()); | |
52 | param -> SetBendingSlope(muon->Py()/muon->Pz()); | |
53 | param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) ); | |
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 | ||
67 | Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2], TLorentzVector &kinem) { | |
68 | ||
69 | if (!(muon->Pz()!=0)) return kFALSE; | |
70 | ||
71 | AliMUONTrackParam *param = new AliMUONTrackParam(); | |
72 | ||
73 | param -> SetNonBendingCoor(muon->XAtDCA()); | |
74 | param -> SetBendingCoor(muon->YAtDCA()); | |
75 | param -> SetZ(AliMFTConstants::fZEvalKinem); | |
76 | param -> SetNonBendingSlope(muon->Px()/muon->Pz()); | |
77 | param -> SetBendingSlope(muon->Py()/muon->Pz()); | |
78 | param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) ); | |
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 | ||
97 | Bool_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 | ||
bffc7f8c | 101 | if (!(muon->Pz()!=0)) return kFALSE; |
102 | ||
103 | AliMUONTrackParam *param = new AliMUONTrackParam(); | |
104 | ||
105 | param -> SetNonBendingCoor(muon->XAtDCA()); | |
106 | param -> SetBendingCoor(muon->YAtDCA()); | |
107 | param -> SetZ(AliMFTConstants::fZEvalKinem); | |
108 | param -> SetNonBendingSlope(muon->Px()/muon->Pz()); | |
109 | param -> SetBendingSlope(muon->Py()/muon->Pz()); | |
110 | param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) ); | |
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 | 133 | Bool_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 | ||
140 | Double_t xPCA=0, yPCA=0, zPCA=0; | |
141 | ||
142 | AliMUONTrackParam *param = new AliMUONTrackParam(); | |
143 | param -> SetNonBendingCoor(muon->XAtDCA()); | |
144 | param -> SetBendingCoor(muon->YAtDCA()); | |
145 | param -> SetZ(AliMFTConstants::fZEvalKinem); | |
146 | param -> SetNonBendingSlope(muon->Px()/muon->Pz()); | |
147 | param -> SetBendingSlope(muon->Py()/muon->Pz()); | |
148 | param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) ); | |
149 | ||
150 | // here we want to understand in which direction we have to search the minimum... | |
151 | ||
152 | Double_t step = 1.; // initial step, in cm | |
153 | Double_t startPoint = 0.; | |
154 | ||
155 | Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step}; | |
156 | ||
157 | 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) | |
158 | ||
159 | for (Int_t i=0; i<3; i++) { | |
160 | AliMUONTrackExtrap::ExtrapToZ(param, z[i]); | |
161 | points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]); | |
162 | points[1] = new TVector3(xy[0],xy[1],z[i]); | |
163 | r[i] = GetDistanceBetweenPoints(points,2); | |
164 | for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu]; | |
165 | } | |
166 | ||
167 | Int_t researchDirection = 0; | |
168 | ||
169 | if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive | |
170 | else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative | |
171 | else if (r[0]<r[1] && r[1]>r[2]) { | |
172 | printf("E-AliMFTAnalysisTools::ExtrapAODMuonToXY: Point of closest approach cannot be found (no minima)\n"); | |
173 | delete param; | |
174 | delete points; | |
175 | return kFALSE; | |
176 | } | |
177 | ||
178 | while (TMath::Abs(researchDirection)>0.5) { | |
179 | ||
180 | if (researchDirection>0) { | |
181 | z[0] = z[1]; | |
182 | z[1] = z[2]; | |
183 | z[2] = z[1]+researchDirection*step; | |
184 | } | |
185 | else { | |
186 | z[2] = z[1]; | |
187 | z[1] = z[0]; | |
188 | z[0] = z[1]+researchDirection*step; | |
189 | } | |
190 | if (TMath::Abs(z[0])>900.) { | |
191 | printf("E-AliMFTAnalysisTools::ExtrapAODMuonToXY: Point of closest approach cannot be found (no minima in the fiducial region)\n"); | |
192 | delete param; | |
193 | delete points; | |
194 | return kFALSE; | |
195 | } | |
196 | ||
197 | for (Int_t i=0; i<3; i++) { | |
198 | AliMUONTrackExtrap::ExtrapToZ(param, z[i]); | |
199 | points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]); | |
200 | points[1] = new TVector3(xy[0],xy[1],z[i]); | |
201 | r[i] = GetDistanceBetweenPoints(points,2); | |
202 | for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu]; | |
203 | } | |
204 | ||
205 | researchDirection=0; | |
206 | if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive | |
207 | else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative | |
208 | ||
209 | } | |
210 | ||
211 | // now we now that the minimum is between z[0] and z[2] and we search for it | |
212 | ||
213 | step *= 0.5; | |
214 | while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) { | |
215 | z[0] = z[1]-step; | |
216 | z[2] = z[1]+step; | |
217 | for (Int_t i=0; i<3; i++) { | |
218 | AliMUONTrackExtrap::ExtrapToZ(param, z[i]); | |
219 | points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]); | |
220 | points[1] = new TVector3(xy[0],xy[1],z[i]); | |
221 | r[i] = GetDistanceBetweenPoints(points,2); | |
222 | for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu]; | |
223 | } | |
224 | if (r[0]<r[1]) z[1] = z[0]; | |
225 | else if (r[2]<r[1]) z[1] = z[2]; | |
226 | else step *= 0.5; | |
227 | } | |
228 | ||
229 | zFinal = z[1]; | |
230 | ||
231 | Double_t xyMuon[2] = {0}; | |
232 | ExtrapAODMuonToZ(muon, zFinal, xyMuon, kinem, cov); | |
233 | ||
234 | return kTRUE; | |
235 | ||
236 | } | |
237 | ||
238 | //==================================================================================================================================================== | |
239 | ||
a8b6b8e9 | 240 | Bool_t AliMFTAnalysisTools::GetAODMuonOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) { |
bffc7f8c | 241 | |
242 | Double_t xy[2] = {0}; | |
243 | ExtrapAODMuonToZ(muon, zv, xy); | |
244 | ||
a8b6b8e9 | 245 | offset = TMath::Sqrt((xv-xy[0])*(xv-xy[0]) + (yv-xy[1])*(yv-xy[1])); |
246 | ||
247 | return kTRUE; | |
bffc7f8c | 248 | |
249 | } | |
250 | ||
251 | //==================================================================================================================================================== | |
252 | ||
9664859a | 253 | Bool_t AliMFTAnalysisTools::GetAODMuonOffsetZ(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) { |
254 | ||
255 | Double_t xy[2] = {xv, yv}; | |
256 | Double_t zFinal = 0; | |
257 | TLorentzVector kinem(0,0,0,0); | |
258 | TMatrixD cov(5,5); | |
259 | ||
260 | ExtrapAODMuonToXY(muon, xy, zFinal, kinem, cov); | |
261 | ||
262 | offset = TMath::Abs(zFinal - zv); | |
263 | ||
264 | return kTRUE; | |
265 | ||
266 | } | |
267 | ||
268 | //==================================================================================================================================================== | |
269 | ||
a8b6b8e9 | 270 | Bool_t AliMFTAnalysisTools::GetAODMuonWeightedOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) { |
bffc7f8c | 271 | |
272 | Double_t xy[2] = {0}; | |
273 | TLorentzVector kinem(0,0,0,0); | |
274 | TMatrixD cov(5,5); | |
275 | ||
276 | ExtrapAODMuonToZ(muon, zv, xy, kinem, cov); | |
277 | ||
278 | TMatrixD covCoordinates(2,2); | |
279 | covCoordinates(0,0) = cov(0,0); | |
280 | covCoordinates(0,1) = cov(0,2); | |
281 | covCoordinates(1,0) = cov(2,0); | |
282 | covCoordinates(1,1) = cov(2,2); | |
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 | ||
304 | Double_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 | ||
322 | Double_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 | 340 | Bool_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 | ||
354 | Bool_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 | } | |
361 | ||
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); | |
07d2587b | 371 | if (TMath::Abs(muon[iMu]->Pz())<1.e-6) { |
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); |
bffc7f8c | 379 | param[iMu] -> SetNonBendingSlope(muon[iMu]->Px()/muon[iMu]->Pz()); |
380 | param[iMu] -> SetBendingSlope(muon[iMu]->Py()/muon[iMu]->Pz()); | |
381 | param[iMu] -> SetInverseBendingMomentum( muon[iMu]->Charge() * (1./muon[iMu]->Pz()) / (TMath::Sqrt(1+TMath::Power(muon[iMu]->Py()/muon[iMu]->Pz(),2))) ); | |
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 | |
bffc7f8c | 389 | Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step}; |
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++) { | |
395 | AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]); | |
396 | points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]); | |
07d2587b | 397 | } |
bffc7f8c | 398 | r[i] = GetDistanceBetweenPoints(points,nMuons); |
07d2587b | 399 | for (Int_t iMu=0; iMu<nMuons; iMu++) delete points[iMu]; |
bffc7f8c | 400 | } |
07d2587b | 401 | |
bffc7f8c | 402 | Int_t researchDirection = 0; |
403 | ||
404 | if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive | |
405 | else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative | |
406 | else if (r[0]<r[1] && r[1]>r[2]) { | |
407 | printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima)\n"); | |
07d2587b | 408 | for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu]; |
409 | delete points; | |
bffc7f8c | 410 | return kFALSE; |
07d2587b | 411 | } |
412 | ||
bffc7f8c | 413 | while (TMath::Abs(researchDirection)>0.5) { |
07d2587b | 414 | |
bffc7f8c | 415 | if (researchDirection>0) { |
416 | z[0] = z[1]; | |
417 | z[1] = z[2]; | |
418 | z[2] = z[1]+researchDirection*step; | |
419 | } | |
420 | else { | |
421 | z[2] = z[1]; | |
422 | z[1] = z[0]; | |
423 | z[0] = z[1]+researchDirection*step; | |
424 | } | |
425 | if (TMath::Abs(z[0])>900.) { | |
426 | printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima in the fiducial region)\n"); | |
07d2587b | 427 | for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu]; |
428 | delete points; | |
bffc7f8c | 429 | return kFALSE; |
430 | } | |
431 | ||
bffc7f8c | 432 | for (Int_t i=0; i<3; i++) { |
433 | for (Int_t iMu=0; iMu<nMuons; iMu++) { | |
07d2587b | 434 | AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]); |
435 | points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]); | |
bffc7f8c | 436 | } |
437 | r[i] = GetDistanceBetweenPoints(points,nMuons); | |
07d2587b | 438 | for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu]; |
bffc7f8c | 439 | } |
440 | researchDirection=0; | |
441 | if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive | |
442 | else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative | |
07d2587b | 443 | |
bffc7f8c | 444 | } |
07d2587b | 445 | |
446 | // now we now that the minimum is between z[0] and z[2] and we search for it | |
bffc7f8c | 447 | |
448 | step *= 0.5; | |
449 | while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) { | |
450 | z[0] = z[1]-step; | |
451 | z[2] = z[1]+step; | |
452 | for (Int_t i=0; i<3; i++) { | |
453 | for (Int_t iMu=0; iMu<nMuons; iMu++) { | |
07d2587b | 454 | AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]); |
455 | points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]); | |
bffc7f8c | 456 | } |
457 | r[i] = GetDistanceBetweenPoints(points,nMuons); | |
07d2587b | 458 | for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu]; |
bffc7f8c | 459 | } |
460 | if (r[0]<r[1]) z[1] = z[0]; | |
461 | else if (r[2]<r[1]) z[1] = z[2]; | |
462 | else step *= 0.5; | |
463 | } | |
464 | ||
bffc7f8c | 465 | // Once z of minimum is found, we evaluate the x and y coordinates by averaging over the contributing tracks |
466 | ||
467 | fZPointOfClosestApproach = z[1]; | |
468 | fXPointOfClosestApproach = 0.; | |
469 | fYPointOfClosestApproach = 0.; | |
470 | for (Int_t iMu=0; iMu<nMuons; iMu++) { | |
471 | AliMUONTrackExtrap::ExtrapToZ(param[iMu], fZPointOfClosestApproach); | |
472 | fXPointOfClosestApproach += param[iMu]->GetNonBendingCoor(); | |
473 | fYPointOfClosestApproach += param[iMu]->GetBendingCoor(); | |
474 | } | |
475 | fXPointOfClosestApproach /= Double_t(nMuons); | |
476 | fYPointOfClosestApproach /= Double_t(nMuons); | |
477 | ||
478 | pca[0] = fXPointOfClosestApproach; | |
479 | pca[1] = fYPointOfClosestApproach; | |
480 | pca[2] = fZPointOfClosestApproach; | |
481 | ||
482 | // Evaluating the kinematics of the N-muon | |
07d2587b | 483 | |
bffc7f8c | 484 | Double_t pTot[3] = {0}; |
485 | Double_t ene = 0.; | |
486 | Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass(); | |
487 | for (Int_t iMu=0; iMu<nMuons; iMu++) { | |
488 | pTot[0] += param[iMu]->Px(); | |
489 | pTot[1] += param[iMu]->Py(); | |
490 | pTot[2] += param[iMu]->Pz(); | |
491 | ene += TMath::Sqrt(massMu*massMu + param[iMu]->Px()*param[iMu]->Px() + param[iMu]->Py()*param[iMu]->Py() + param[iMu]->Pz()*param[iMu]->Pz()); | |
492 | } | |
493 | ||
494 | kinem.SetPxPyPzE(pTot[0], pTot[1], pTot[2], ene); | |
07d2587b | 495 | |
bffc7f8c | 496 | // Evaluating the PCA quality of the N-muon |
07d2587b | 497 | |
bffc7f8c | 498 | Double_t sum=0.,squareSum=0.; |
499 | for (Int_t iMu=0; iMu<nMuons; iMu++) { | |
a8b6b8e9 | 500 | Double_t wOffset = 0; |
501 | if (!GetAODMuonWeightedOffset(muon[iMu],fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach, wOffset)) { | |
502 | for(Int_t jMu=0;jMu<nMuons;jMu++) delete param[jMu]; | |
503 | delete points; | |
504 | return kFALSE; | |
505 | } | |
bffc7f8c | 506 | Double_t f = TMath::Exp(-0.5 * wOffset); |
507 | sum += f; | |
508 | squareSum += f*f; | |
509 | } | |
510 | if (sum > 0.) pcaQuality = (sum-squareSum/sum) / (nMuons-1); | |
511 | else pcaQuality = 0.; | |
07d2587b | 512 | |
513 | for(Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu]; | |
514 | delete points; | |
bffc7f8c | 515 | return kTRUE; |
516 | ||
517 | } | |
518 | ||
519 | //========================================================================================================================= | |
520 | ||
521 | Double_t AliMFTAnalysisTools::GetDistanceBetweenPoints(TVector3 **points, Int_t nPoints) { | |
522 | ||
523 | if (nPoints>AliMFTConstants::fNMaxMuonsForPCA) { | |
524 | printf("W-AliMFTAnalysisTools::GetDistanceBetweenPoints: number of points not valid\n"); | |
525 | return 1.e9; | |
526 | } | |
527 | ||
528 | if (nPoints<2) return 0.; | |
529 | if (nPoints<3) return TMath::Sqrt( (points[0]->X()-points[1]->X()) * (points[0]->X()-points[1]->X()) + | |
530 | (points[0]->Y()-points[1]->Y()) * (points[0]->Y()-points[1]->Y()) + | |
531 | (points[0]->Z()-points[1]->Z()) * (points[0]->Z()-points[1]->Z()) ); | |
532 | ||
533 | const Int_t nEdgesMax = ((AliMFTConstants::fNMaxMuonsForPCA) * (AliMFTConstants::fNMaxMuonsForPCA - 1)) / 2; | |
534 | ||
535 | Int_t startID[nEdgesMax] = {0}; | |
536 | Int_t stopID[nEdgesMax] = {0}; | |
537 | Double_t edgeLength[nEdgesMax] = {0}; | |
538 | ||
539 | Bool_t pointStatus[AliMFTConstants::fNMaxMuonsForPCA] = {0}; | |
540 | ||
541 | Int_t nEdges=0; | |
542 | for (Int_t i=0; i<nPoints-1; i++) { | |
543 | for (Int_t j=i+1; j<nPoints; j++) { | |
544 | edgeLength[nEdges] = TMath::Sqrt( (points[i]->X()-points[j]->X()) * (points[i]->X()-points[j]->X()) + | |
545 | (points[i]->Y()-points[j]->Y()) * (points[i]->Y()-points[j]->Y()) + | |
546 | (points[i]->Z()-points[j]->Z()) * (points[i]->Z()-points[j]->Z()) ); | |
547 | stopID[nEdges] = i; | |
548 | startID[nEdges] = j; | |
549 | nEdges++; | |
550 | } | |
551 | } | |
552 | ||
553 | // Order Edges | |
554 | ||
555 | Double_t min = 0; | |
556 | Int_t iMin = 0; | |
557 | ||
558 | for (Int_t iEdge=0; iEdge<nEdges-1; iEdge++) { | |
559 | min = edgeLength[iEdge]; | |
560 | iMin = iEdge; | |
561 | for (Int_t j=iEdge+1; j<nEdges; j++) { | |
562 | if (edgeLength[j]<min) { | |
563 | min = edgeLength[j]; | |
564 | iMin = j; | |
565 | } | |
566 | } | |
567 | ||
568 | if (iMin != iEdge) { | |
569 | ||
570 | Double_t edgeLengthMin = edgeLength[iMin]; | |
571 | Int_t startIDmin = startID[iMin]; | |
572 | Int_t stopIDmin = stopID[iMin]; | |
573 | ||
574 | edgeLength[iMin] = edgeLength[iEdge]; | |
575 | startID[iMin] = startID[iEdge]; | |
576 | stopID[iMin] = stopID[iEdge]; | |
577 | ||
578 | edgeLength[iEdge] = edgeLengthMin; | |
579 | startID[iEdge] = startIDmin; | |
580 | stopID[iEdge] = stopIDmin; | |
581 | ||
582 | } | |
583 | ||
584 | } | |
585 | ||
586 | // Connect | |
587 | ||
588 | Double_t length = 0.; | |
589 | for (Int_t i=0; i<nEdges; i++) { | |
590 | if (!(pointStatus[startID[i]] && pointStatus[stopID[i]])) { | |
591 | pointStatus[startID[i]] = kTRUE; | |
592 | pointStatus[stopID[i]] = kTRUE; | |
593 | length += edgeLength[i]; | |
594 | } | |
595 | } | |
596 | ||
597 | return length; | |
598 | ||
599 | } | |
600 | ||
601 | //==================================================================================================================================================== | |
602 | ||
603 | void AliMFTAnalysisTools::ConvertCovMatrixMUON2AOD(const TMatrixD& covMUON, Double_t covAOD[21]) { | |
604 | ||
605 | // Converts the cov matrix from the MUON format (TMatrixD) to the AOD one (Double_t[21]) | |
606 | // | |
607 | // Cov(x,x) ... : cv[0] | |
608 | // Cov(x,slopeX) ... : cv[1] cv[2] | |
609 | // Cov(x,y) ... : cv[3] cv[4] cv[5] | |
610 | // Cov(x,slopeY) ... : cv[6] cv[7] cv[8] cv[9] | |
611 | // Cov(x,invP_yz) ... : cv[10] cv[11] cv[12] cv[13] cv[14] | |
612 | // not-used ... : cv[15] cv[16] cv[17] cv[18] cv[19] cv[20] | |
613 | ||
614 | covAOD[0] = covMUON(0,0); | |
615 | ||
616 | covAOD[1] = covMUON(1,0); | |
617 | covAOD[2] = covMUON(1,1); | |
618 | ||
619 | covAOD[3] = covMUON(2,0); | |
620 | covAOD[4] = covMUON(2,1); | |
621 | covAOD[5] = covMUON(2,2); | |
622 | ||
623 | covAOD[6] = covMUON(3,0); | |
624 | covAOD[7] = covMUON(3,1); | |
625 | covAOD[8] = covMUON(3,2); | |
626 | covAOD[9] = covMUON(3,3); | |
627 | ||
628 | covAOD[10] = covMUON(4,0); | |
629 | covAOD[11] = covMUON(4,1); | |
630 | covAOD[12] = covMUON(4,2); | |
631 | covAOD[13] = covMUON(4,3); | |
632 | covAOD[14] = covMUON(4,4); | |
633 | ||
634 | covAOD[15] = 0; | |
635 | covAOD[16] = 0; | |
636 | covAOD[17] = 0; | |
637 | covAOD[18] = 0; | |
638 | covAOD[19] = 0; | |
639 | covAOD[20] = 0; | |
640 | ||
641 | } | |
642 | ||
643 | //==================================================================================================================================================== | |
644 | ||
645 | const TMatrixD AliMFTAnalysisTools::ConvertCovMatrixAOD2MUON(AliAODTrack *muon) { | |
646 | ||
647 | Double_t covAOD[21] = {0}; | |
648 | muon -> GetCovarianceXYZPxPyPz(covAOD); | |
649 | ||
650 | TMatrixD covMUON(5,5); | |
651 | ||
652 | covMUON(0,0) = covAOD[0]; | |
653 | ||
654 | covMUON(1,0) = covAOD[1]; | |
655 | covMUON(1,1) = covAOD[2]; | |
656 | ||
657 | covMUON(2,0) = covAOD[3]; | |
658 | covMUON(2,1) = covAOD[4]; | |
659 | covMUON(2,2) = covAOD[5]; | |
660 | ||
661 | covMUON(3,0) = covAOD[6]; | |
662 | covMUON(3,1) = covAOD[7]; | |
663 | covMUON(3,2) = covAOD[8]; | |
664 | covMUON(3,3) = covAOD[9]; | |
665 | ||
666 | covMUON(4,0) = covAOD[10]; | |
667 | covMUON(4,1) = covAOD[11]; | |
668 | covMUON(4,2) = covAOD[12]; | |
669 | covMUON(4,3) = covAOD[13]; | |
670 | covMUON(4,4) = covAOD[14]; | |
671 | ||
672 | return covMUON; | |
673 | ||
674 | } | |
675 | ||
676 | //==================================================================================================================================================== | |
677 |