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