]>
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 | ||
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 | ||
67 | Bool_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 | ||
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 | ||
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 | 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); | |
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 | 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 | ||
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 | ||
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 | } | |
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 | ||
536 | Double_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 | ||
618 | void 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 | ||
660 | const 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 | 693 | Bool_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 | ||
702 | TString 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 | ||
725 | void 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 | 755 | Bool_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 | 771 | Bool_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 | ||
802 | Bool_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 | |
812 | Bool_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 | ||
823 | Bool_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 | ||
834 | Bool_t AliMFTAnalysisTools::IsPDGResonance(Int_t pdgCode) { | |
835 | ||
836 | Int_t id = pdgCode%100000; | |
837 | return (!((id-id%10)%110)); | |
838 | ||
839 | } | |
840 | ||
841 | //==================================================================================================================================================== |