d65930879fbd6c4891ae717ab87adfde716e9cb2
[u/mrichter/AliRoot.git] / MFT / AliMFTAnalysisTools.cxx
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
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
99   // Extrapolate muon to a given z providing the corresponding (x,y) position and updating kinematics and covariance matrix
100
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
133 Bool_t AliMFTAnalysisTools::ExtrapAODMuonToXY(AliAODTrack *muon, Double_t xy[2], Double_t &zFinal, TLorentzVector &kinem, TMatrixD &cov) {
134
135   // Find the point of closest approach between the muon and the direction parallel to the z-axis defined by the given (x,y)
136   // Provide the z of the above point as weel as the updated kinematics and covariance matrix
137
138   // We look for the above-defined PCA
139
140   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
238 Bool_t AliMFTAnalysisTools::GetAODMuonOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
239
240   Double_t xy[2] = {0};
241   ExtrapAODMuonToZ(muon, zv, xy);
242   
243   offset = TMath::Sqrt((xv-xy[0])*(xv-xy[0]) + (yv-xy[1])*(yv-xy[1]));
244
245   return kTRUE;
246
247 }
248
249 //====================================================================================================================================================
250
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
268 Bool_t AliMFTAnalysisTools::GetAODMuonWeightedOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
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   
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     
288     offset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) +
289                               dY*dY*covCoordinatesInverse(1,1) +
290                               2.*dX*dY*covCoordinatesInverse(0,1)));
291     
292     return kTRUE;
293
294   }
295   
296   return kFALSE;
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
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));
343   
344   Bool_t result = CalculatePCA(muons, pca, pcaQuality, kinem);
345   delete muons;
346   return result;
347
348 }
349
350 //====================================================================================================================================================
351
352 Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
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;
361
362   AliAODTrack *muon[AliMFTConstants::fNMaxMuonsForPCA]        = {0};
363   AliMUONTrackParam *param[AliMFTConstants::fNMaxMuonsForPCA] = {0};
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);
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     }
373     param[iMu] = new AliMUONTrackParam();
374     param[iMu] -> SetNonBendingCoor(muon[iMu]->XAtDCA());
375     param[iMu] -> SetBendingCoor(muon[iMu]->YAtDCA());
376     param[iMu] -> SetZ(AliMFTConstants::fZEvalKinem);
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   }
381   
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.;
386   
387   Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step};
388   
389   TVector3 **points = new TVector3*[AliMFTConstants::fNMaxMuonsForPCA];
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]);
395                 }
396     r[i] = GetDistanceBetweenPoints(points,nMuons);
397     for (Int_t iMu=0; iMu<nMuons; iMu++) delete points[iMu]; 
398   }
399   
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");
406     for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
407     delete points;
408     return kFALSE;
409   }     
410   
411   while (TMath::Abs(researchDirection)>0.5) {
412       
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");
425       for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
426       delete points;
427       return kFALSE;
428     }
429     
430     for (Int_t i=0; i<3; i++) {
431       for (Int_t iMu=0; iMu<nMuons; iMu++) {
432         AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
433         points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
434       }      
435       r[i] = GetDistanceBetweenPoints(points,nMuons);
436       for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu];
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
441     
442   }
443   
444   // now we now that the minimum is between z[0] and z[2] and we search for it
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++) {
452         AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
453         points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
454       }      
455       r[i] = GetDistanceBetweenPoints(points,nMuons);
456       for (Int_t iMu=0;iMu<nMuons;iMu++)        delete points[iMu];
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   
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
481   
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);
493   
494   // Evaluating the PCA quality of the N-muon
495   
496   Double_t sum=0.,squareSum=0.;
497   for (Int_t iMu=0; iMu<nMuons; iMu++) {
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     }
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.;
510   
511   for(Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
512   delete points;
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
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