]>
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 | ||
9664859a | 140 | AliMUONTrackParam *param = new AliMUONTrackParam(); |
141 | param -> SetNonBendingCoor(muon->XAtDCA()); | |
142 | param -> SetBendingCoor(muon->YAtDCA()); | |
143 | param -> SetZ(AliMFTConstants::fZEvalKinem); | |
144 | param -> SetNonBendingSlope(muon->Px()/muon->Pz()); | |
145 | param -> SetBendingSlope(muon->Py()/muon->Pz()); | |
146 | param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) ); | |
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 | ||
153 | Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step}; | |
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 | 238 | Bool_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 | 251 | Bool_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 | 268 | Bool_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 | ||
7bafb008 | 282 | if (TDecompLU::InvertLU(covCoordinates,covCoordinates.GetTol(),0)) { |
283 | ||
284 | TMatrixD covCoordinatesInverse = covCoordinates; | |
285 | Double_t dX = xy[0] - xv; | |
286 | Double_t dY = xy[1] - yv; | |
287 | ||
a8b6b8e9 | 288 | offset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) + |
289 | dY*dY*covCoordinatesInverse(1,1) + | |
290 | 2.*dX*dY*covCoordinatesInverse(0,1))); | |
7bafb008 | 291 | |
a8b6b8e9 | 292 | return kTRUE; |
7bafb008 | 293 | |
294 | } | |
bffc7f8c | 295 | |
a8b6b8e9 | 296 | return kFALSE; |
bffc7f8c | 297 | |
298 | } | |
299 | ||
300 | //==================================================================================================================================================== | |
301 | ||
302 | Double_t AliMFTAnalysisTools::GetPseudoProperDecayTimeXY(Double_t xVtx, Double_t yVtx, | |
303 | Double_t xDimu, Double_t yDimu, | |
304 | Double_t mDimu, Double_t ptDimu) { | |
305 | ||
306 | // pseudo-proper decay time of a particle produced in the primary vertex and decaying into a dimuon (+ X) | |
307 | // evaluated using the transverse degree of freedom of the decay topology | |
308 | ||
309 | if (ptDimu != 0) { | |
310 | Double_t decayLengthXY = TMath::Sqrt((xVtx-xDimu)*(xVtx-xDimu)+(yVtx-yDimu)*(yVtx-yDimu)); | |
311 | return (decayLengthXY * mDimu/ptDimu)/TMath::Ccgs()*1E12; // in ps | |
312 | } | |
313 | ||
314 | return -99999999; | |
315 | ||
316 | } | |
317 | ||
318 | //==================================================================================================================================================== | |
319 | ||
320 | Double_t AliMFTAnalysisTools::GetPseudoProperDecayTimeZ(Double_t zVtx, | |
321 | Double_t zDimu, | |
322 | Double_t mDimu, Double_t pzDimu) { | |
323 | ||
324 | // pseudo-proper decay time of a particle produced in the primary vertex and decaying into a dimuon (+ X) | |
325 | // evaluated using the longitudinal degree of freedom of the decay topology | |
326 | ||
327 | if (pzDimu != 0) { | |
328 | Double_t decayLengthZ = zDimu - zVtx; | |
329 | return (decayLengthZ * mDimu/pzDimu)/TMath::Ccgs()*1E12; // in ps | |
330 | } | |
331 | ||
332 | return -99999999; | |
333 | ||
334 | } | |
335 | ||
336 | //==================================================================================================================================================== | |
337 | ||
19393921 | 338 | Bool_t AliMFTAnalysisTools::CalculatePCA(AliAODDimuon *dimuon, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) { |
339 | ||
340 | TObjArray *muons = new TObjArray(); | |
341 | muons -> Add(dimuon->GetMu(0)); | |
342 | muons -> Add(dimuon->GetMu(1)); | |
07d2587b | 343 | |
344 | Bool_t result = CalculatePCA(muons, pca, pcaQuality, kinem); | |
345 | delete muons; | |
346 | return result; | |
19393921 | 347 | |
348 | } | |
349 | ||
350 | //==================================================================================================================================================== | |
351 | ||
352 | Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) { | |
bffc7f8c | 353 | |
354 | const Int_t nMuons = muons->GetEntriesFast(); | |
355 | if (nMuons<2 || nMuons>AliMFTConstants::fNMaxMuonsForPCA) { | |
356 | printf("W-AliMFTAnalysisTools::CalculatePCA: number of muons not valid\n"); | |
357 | return kFALSE; | |
358 | } | |
359 | ||
360 | Double_t fXPointOfClosestApproach=0, fYPointOfClosestApproach=0, fZPointOfClosestApproach=0; | |
07d2587b | 361 | |
6ffd24bb | 362 | AliAODTrack *muon[AliMFTConstants::fNMaxMuonsForPCA] = {0}; |
363 | AliMUONTrackParam *param[AliMFTConstants::fNMaxMuonsForPCA] = {0}; | |
bffc7f8c | 364 | |
365 | // Finding AliMUONTrackParam objects for each muon | |
366 | ||
367 | for (Int_t iMu=0; iMu<nMuons; iMu++) { | |
368 | muon[iMu] = (AliAODTrack*) muons->At(iMu); | |
07d2587b | 369 | if (TMath::Abs(muon[iMu]->Pz())<1.e-6) { |
370 | for(Int_t i=0;i<iMu;i++) delete param[i]; | |
371 | return kFALSE; | |
372 | } | |
bffc7f8c | 373 | param[iMu] = new AliMUONTrackParam(); |
374 | param[iMu] -> SetNonBendingCoor(muon[iMu]->XAtDCA()); | |
375 | param[iMu] -> SetBendingCoor(muon[iMu]->YAtDCA()); | |
9664859a | 376 | param[iMu] -> SetZ(AliMFTConstants::fZEvalKinem); |
bffc7f8c | 377 | param[iMu] -> SetNonBendingSlope(muon[iMu]->Px()/muon[iMu]->Pz()); |
378 | param[iMu] -> SetBendingSlope(muon[iMu]->Py()/muon[iMu]->Pz()); | |
379 | param[iMu] -> SetInverseBendingMomentum( muon[iMu]->Charge() * (1./muon[iMu]->Pz()) / (TMath::Sqrt(1+TMath::Power(muon[iMu]->Py()/muon[iMu]->Pz(),2))) ); | |
380 | } | |
07d2587b | 381 | |
bffc7f8c | 382 | // here we want to understand in which direction we have to search the minimum... |
383 | ||
384 | Double_t step = 1.; // initial step, in cm | |
385 | Double_t startPoint = 0.; | |
07d2587b | 386 | |
bffc7f8c | 387 | Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step}; |
388 | ||
6ffd24bb | 389 | TVector3 **points = new TVector3*[AliMFTConstants::fNMaxMuonsForPCA]; |
bffc7f8c | 390 | |
391 | for (Int_t i=0; i<3; i++) { | |
392 | for (Int_t iMu=0; iMu<nMuons; iMu++) { | |
393 | AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]); | |
394 | points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]); | |
07d2587b | 395 | } |
bffc7f8c | 396 | r[i] = GetDistanceBetweenPoints(points,nMuons); |
07d2587b | 397 | for (Int_t iMu=0; iMu<nMuons; iMu++) delete points[iMu]; |
bffc7f8c | 398 | } |
07d2587b | 399 | |
bffc7f8c | 400 | Int_t researchDirection = 0; |
401 | ||
402 | if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive | |
403 | else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative | |
404 | else if (r[0]<r[1] && r[1]>r[2]) { | |
405 | printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima)\n"); | |
07d2587b | 406 | for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu]; |
407 | delete points; | |
bffc7f8c | 408 | return kFALSE; |
07d2587b | 409 | } |
410 | ||
bffc7f8c | 411 | while (TMath::Abs(researchDirection)>0.5) { |
07d2587b | 412 | |
bffc7f8c | 413 | if (researchDirection>0) { |
414 | z[0] = z[1]; | |
415 | z[1] = z[2]; | |
416 | z[2] = z[1]+researchDirection*step; | |
417 | } | |
418 | else { | |
419 | z[2] = z[1]; | |
420 | z[1] = z[0]; | |
421 | z[0] = z[1]+researchDirection*step; | |
422 | } | |
423 | if (TMath::Abs(z[0])>900.) { | |
424 | printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima in the fiducial region)\n"); | |
07d2587b | 425 | for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu]; |
426 | delete points; | |
bffc7f8c | 427 | return kFALSE; |
428 | } | |
429 | ||
bffc7f8c | 430 | for (Int_t i=0; i<3; i++) { |
431 | for (Int_t iMu=0; iMu<nMuons; iMu++) { | |
07d2587b | 432 | AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]); |
433 | points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]); | |
bffc7f8c | 434 | } |
435 | r[i] = GetDistanceBetweenPoints(points,nMuons); | |
07d2587b | 436 | for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu]; |
bffc7f8c | 437 | } |
438 | researchDirection=0; | |
439 | if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive | |
440 | else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative | |
07d2587b | 441 | |
bffc7f8c | 442 | } |
07d2587b | 443 | |
444 | // now we now that the minimum is between z[0] and z[2] and we search for it | |
bffc7f8c | 445 | |
446 | step *= 0.5; | |
447 | while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) { | |
448 | z[0] = z[1]-step; | |
449 | z[2] = z[1]+step; | |
450 | for (Int_t i=0; i<3; i++) { | |
451 | for (Int_t iMu=0; iMu<nMuons; iMu++) { | |
07d2587b | 452 | AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]); |
453 | points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]); | |
bffc7f8c | 454 | } |
455 | r[i] = GetDistanceBetweenPoints(points,nMuons); | |
07d2587b | 456 | for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu]; |
bffc7f8c | 457 | } |
458 | if (r[0]<r[1]) z[1] = z[0]; | |
459 | else if (r[2]<r[1]) z[1] = z[2]; | |
460 | else step *= 0.5; | |
461 | } | |
462 | ||
bffc7f8c | 463 | // Once z of minimum is found, we evaluate the x and y coordinates by averaging over the contributing tracks |
464 | ||
465 | fZPointOfClosestApproach = z[1]; | |
466 | fXPointOfClosestApproach = 0.; | |
467 | fYPointOfClosestApproach = 0.; | |
468 | for (Int_t iMu=0; iMu<nMuons; iMu++) { | |
469 | AliMUONTrackExtrap::ExtrapToZ(param[iMu], fZPointOfClosestApproach); | |
470 | fXPointOfClosestApproach += param[iMu]->GetNonBendingCoor(); | |
471 | fYPointOfClosestApproach += param[iMu]->GetBendingCoor(); | |
472 | } | |
473 | fXPointOfClosestApproach /= Double_t(nMuons); | |
474 | fYPointOfClosestApproach /= Double_t(nMuons); | |
475 | ||
476 | pca[0] = fXPointOfClosestApproach; | |
477 | pca[1] = fYPointOfClosestApproach; | |
478 | pca[2] = fZPointOfClosestApproach; | |
479 | ||
480 | // Evaluating the kinematics of the N-muon | |
07d2587b | 481 | |
bffc7f8c | 482 | Double_t pTot[3] = {0}; |
483 | Double_t ene = 0.; | |
484 | Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass(); | |
485 | for (Int_t iMu=0; iMu<nMuons; iMu++) { | |
486 | pTot[0] += param[iMu]->Px(); | |
487 | pTot[1] += param[iMu]->Py(); | |
488 | pTot[2] += param[iMu]->Pz(); | |
489 | ene += TMath::Sqrt(massMu*massMu + param[iMu]->Px()*param[iMu]->Px() + param[iMu]->Py()*param[iMu]->Py() + param[iMu]->Pz()*param[iMu]->Pz()); | |
490 | } | |
491 | ||
492 | kinem.SetPxPyPzE(pTot[0], pTot[1], pTot[2], ene); | |
07d2587b | 493 | |
bffc7f8c | 494 | // Evaluating the PCA quality of the N-muon |
07d2587b | 495 | |
bffc7f8c | 496 | Double_t sum=0.,squareSum=0.; |
497 | for (Int_t iMu=0; iMu<nMuons; iMu++) { | |
a8b6b8e9 | 498 | Double_t wOffset = 0; |
499 | if (!GetAODMuonWeightedOffset(muon[iMu],fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach, wOffset)) { | |
500 | for(Int_t jMu=0;jMu<nMuons;jMu++) delete param[jMu]; | |
501 | delete points; | |
502 | return kFALSE; | |
503 | } | |
bffc7f8c | 504 | Double_t f = TMath::Exp(-0.5 * wOffset); |
505 | sum += f; | |
506 | squareSum += f*f; | |
507 | } | |
508 | if (sum > 0.) pcaQuality = (sum-squareSum/sum) / (nMuons-1); | |
509 | else pcaQuality = 0.; | |
07d2587b | 510 | |
511 | for(Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu]; | |
512 | delete points; | |
bffc7f8c | 513 | return kTRUE; |
514 | ||
515 | } | |
516 | ||
517 | //========================================================================================================================= | |
518 | ||
519 | Double_t AliMFTAnalysisTools::GetDistanceBetweenPoints(TVector3 **points, Int_t nPoints) { | |
520 | ||
521 | if (nPoints>AliMFTConstants::fNMaxMuonsForPCA) { | |
522 | printf("W-AliMFTAnalysisTools::GetDistanceBetweenPoints: number of points not valid\n"); | |
523 | return 1.e9; | |
524 | } | |
525 | ||
526 | if (nPoints<2) return 0.; | |
527 | if (nPoints<3) return TMath::Sqrt( (points[0]->X()-points[1]->X()) * (points[0]->X()-points[1]->X()) + | |
528 | (points[0]->Y()-points[1]->Y()) * (points[0]->Y()-points[1]->Y()) + | |
529 | (points[0]->Z()-points[1]->Z()) * (points[0]->Z()-points[1]->Z()) ); | |
530 | ||
531 | const Int_t nEdgesMax = ((AliMFTConstants::fNMaxMuonsForPCA) * (AliMFTConstants::fNMaxMuonsForPCA - 1)) / 2; | |
532 | ||
533 | Int_t startID[nEdgesMax] = {0}; | |
534 | Int_t stopID[nEdgesMax] = {0}; | |
535 | Double_t edgeLength[nEdgesMax] = {0}; | |
536 | ||
537 | Bool_t pointStatus[AliMFTConstants::fNMaxMuonsForPCA] = {0}; | |
538 | ||
539 | Int_t nEdges=0; | |
540 | for (Int_t i=0; i<nPoints-1; i++) { | |
541 | for (Int_t j=i+1; j<nPoints; j++) { | |
542 | edgeLength[nEdges] = TMath::Sqrt( (points[i]->X()-points[j]->X()) * (points[i]->X()-points[j]->X()) + | |
543 | (points[i]->Y()-points[j]->Y()) * (points[i]->Y()-points[j]->Y()) + | |
544 | (points[i]->Z()-points[j]->Z()) * (points[i]->Z()-points[j]->Z()) ); | |
545 | stopID[nEdges] = i; | |
546 | startID[nEdges] = j; | |
547 | nEdges++; | |
548 | } | |
549 | } | |
550 | ||
551 | // Order Edges | |
552 | ||
553 | Double_t min = 0; | |
554 | Int_t iMin = 0; | |
555 | ||
556 | for (Int_t iEdge=0; iEdge<nEdges-1; iEdge++) { | |
557 | min = edgeLength[iEdge]; | |
558 | iMin = iEdge; | |
559 | for (Int_t j=iEdge+1; j<nEdges; j++) { | |
560 | if (edgeLength[j]<min) { | |
561 | min = edgeLength[j]; | |
562 | iMin = j; | |
563 | } | |
564 | } | |
565 | ||
566 | if (iMin != iEdge) { | |
567 | ||
568 | Double_t edgeLengthMin = edgeLength[iMin]; | |
569 | Int_t startIDmin = startID[iMin]; | |
570 | Int_t stopIDmin = stopID[iMin]; | |
571 | ||
572 | edgeLength[iMin] = edgeLength[iEdge]; | |
573 | startID[iMin] = startID[iEdge]; | |
574 | stopID[iMin] = stopID[iEdge]; | |
575 | ||
576 | edgeLength[iEdge] = edgeLengthMin; | |
577 | startID[iEdge] = startIDmin; | |
578 | stopID[iEdge] = stopIDmin; | |
579 | ||
580 | } | |
581 | ||
582 | } | |
583 | ||
584 | // Connect | |
585 | ||
586 | Double_t length = 0.; | |
587 | for (Int_t i=0; i<nEdges; i++) { | |
588 | if (!(pointStatus[startID[i]] && pointStatus[stopID[i]])) { | |
589 | pointStatus[startID[i]] = kTRUE; | |
590 | pointStatus[stopID[i]] = kTRUE; | |
591 | length += edgeLength[i]; | |
592 | } | |
593 | } | |
594 | ||
595 | return length; | |
596 | ||
597 | } | |
598 | ||
599 | //==================================================================================================================================================== | |
600 | ||
601 | void AliMFTAnalysisTools::ConvertCovMatrixMUON2AOD(const TMatrixD& covMUON, Double_t covAOD[21]) { | |
602 | ||
603 | // Converts the cov matrix from the MUON format (TMatrixD) to the AOD one (Double_t[21]) | |
604 | // | |
605 | // Cov(x,x) ... : cv[0] | |
606 | // Cov(x,slopeX) ... : cv[1] cv[2] | |
607 | // Cov(x,y) ... : cv[3] cv[4] cv[5] | |
608 | // Cov(x,slopeY) ... : cv[6] cv[7] cv[8] cv[9] | |
609 | // Cov(x,invP_yz) ... : cv[10] cv[11] cv[12] cv[13] cv[14] | |
610 | // not-used ... : cv[15] cv[16] cv[17] cv[18] cv[19] cv[20] | |
611 | ||
612 | covAOD[0] = covMUON(0,0); | |
613 | ||
614 | covAOD[1] = covMUON(1,0); | |
615 | covAOD[2] = covMUON(1,1); | |
616 | ||
617 | covAOD[3] = covMUON(2,0); | |
618 | covAOD[4] = covMUON(2,1); | |
619 | covAOD[5] = covMUON(2,2); | |
620 | ||
621 | covAOD[6] = covMUON(3,0); | |
622 | covAOD[7] = covMUON(3,1); | |
623 | covAOD[8] = covMUON(3,2); | |
624 | covAOD[9] = covMUON(3,3); | |
625 | ||
626 | covAOD[10] = covMUON(4,0); | |
627 | covAOD[11] = covMUON(4,1); | |
628 | covAOD[12] = covMUON(4,2); | |
629 | covAOD[13] = covMUON(4,3); | |
630 | covAOD[14] = covMUON(4,4); | |
631 | ||
632 | covAOD[15] = 0; | |
633 | covAOD[16] = 0; | |
634 | covAOD[17] = 0; | |
635 | covAOD[18] = 0; | |
636 | covAOD[19] = 0; | |
637 | covAOD[20] = 0; | |
638 | ||
639 | } | |
640 | ||
641 | //==================================================================================================================================================== | |
642 | ||
643 | const TMatrixD AliMFTAnalysisTools::ConvertCovMatrixAOD2MUON(AliAODTrack *muon) { | |
644 | ||
645 | Double_t covAOD[21] = {0}; | |
646 | muon -> GetCovarianceXYZPxPyPz(covAOD); | |
647 | ||
648 | TMatrixD covMUON(5,5); | |
649 | ||
650 | covMUON(0,0) = covAOD[0]; | |
651 | ||
652 | covMUON(1,0) = covAOD[1]; | |
653 | covMUON(1,1) = covAOD[2]; | |
654 | ||
655 | covMUON(2,0) = covAOD[3]; | |
656 | covMUON(2,1) = covAOD[4]; | |
657 | covMUON(2,2) = covAOD[5]; | |
658 | ||
659 | covMUON(3,0) = covAOD[6]; | |
660 | covMUON(3,1) = covAOD[7]; | |
661 | covMUON(3,2) = covAOD[8]; | |
662 | covMUON(3,3) = covAOD[9]; | |
663 | ||
664 | covMUON(4,0) = covAOD[10]; | |
665 | covMUON(4,1) = covAOD[11]; | |
666 | covMUON(4,2) = covAOD[12]; | |
667 | covMUON(4,3) = covAOD[13]; | |
668 | covMUON(4,4) = covAOD[14]; | |
669 | ||
670 | return covMUON; | |
671 | ||
672 | } | |
673 | ||
674 | //==================================================================================================================================================== | |
675 | ||
0ebb458a | 676 | Bool_t AliMFTAnalysisTools::IsCorrectMatch(AliAODTrack *muon) { |
677 | ||
678 | for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) if (IsWrongCluster(muon, iPlane)) return kFALSE; | |
679 | return kTRUE; | |
680 | ||
681 | } | |
682 | ||
683 | //==================================================================================================================================================== | |
684 | ||
685 | TString AliMFTAnalysisTools::GetGenerator(Int_t label, AliAODMCHeader* header) { | |
686 | ||
687 | // get the name of the generator that produced a given particle | |
688 | ||
689 | Int_t partCounter = 0; | |
690 | TList *genHeaders = header->GetCocktailHeaders(); | |
691 | Int_t nGenHeaders = genHeaders->GetEntries(); | |
692 | ||
693 | for (Int_t i=0; i<nGenHeaders; i++){ | |
694 | AliGenEventHeader *gh = (AliGenEventHeader*) genHeaders->At(i); | |
695 | TString genName = gh->GetName(); | |
696 | Int_t nPart = gh->NProduced(); | |
697 | if (label>=partCounter && label<(partCounter+nPart)) return genName; | |
698 | partCounter += nPart; | |
699 | } | |
700 | ||
701 | TString empty=""; | |
702 | return empty; | |
703 | ||
704 | } | |
705 | ||
706 | //==================================================================================================================================================== | |
707 | ||
708 | void AliMFTAnalysisTools::GetTrackPrimaryGenerator(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC, TString &nameGen) { | |
709 | ||
710 | // method to check if a track comes from a given generator | |
711 | ||
712 | Int_t label = TMath::Abs(track->GetLabel()); | |
713 | nameGen = GetGenerator(label,header); | |
714 | ||
715 | // 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 | |
716 | ||
717 | while (nameGen.IsWhitespace()) { | |
718 | AliAODMCParticle *mcPart = (AliAODMCParticle*) arrayMC->At(label); | |
719 | if (!mcPart) { | |
720 | printf("AliMFTAnalysisTools::GetTrackPrimaryGenerator - BREAK: No valid AliAODMCParticle at label %i\n",label); | |
721 | break; | |
722 | } | |
723 | Int_t motherLabel = mcPart->GetMother(); | |
724 | if (motherLabel < 0) { | |
725 | printf("AliMFTAnalysisTools::GetTrackPrimaryGenerator - BREAK: Reached primary particle without valid mother\n"); | |
726 | break; | |
727 | } | |
728 | label = motherLabel; | |
729 | nameGen = GetGenerator(label,header); | |
730 | } | |
731 | ||
732 | return; | |
733 | ||
734 | } | |
735 | ||
736 | //==================================================================================================================================================== | |
737 | ||
738 | Bool_t AliMFTAnalysisTools::IsTrackInjected(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC){ | |
739 | ||
740 | // method to check if a track comes from the signal event or from the underlying Hijing event | |
741 | ||
742 | TString nameGen; | |
743 | ||
744 | GetTrackPrimaryGenerator(track,header,arrayMC,nameGen); | |
745 | ||
746 | if (nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE; | |
747 | ||
748 | return kTRUE; | |
749 | } | |
750 | ||
751 | //==================================================================================================================================================== | |
752 |