da4132c56c74855e126be48174a4d0fd45cc534e
[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->PzAtDCA()!=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->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))) );
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->PzAtDCA()!=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->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))) );
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->PzAtDCA()!=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->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))) );
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->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))) );
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-step, startPoint, startPoint+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 (covCoordinates.Determinant() < covCoordinates.GetTol()) return kFALSE;
283
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     
290     offset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) +
291                               dY*dY*covCoordinatesInverse(1,1) +
292                               2.*dX*dY*covCoordinatesInverse(0,1)));
293     
294     return kTRUE;
295
296   }
297   
298   return kFALSE;
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
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));
345   
346   Bool_t result = CalculatePCA(muons, pca, pcaQuality, kinem);
347   delete muons;
348   return result;
349
350 }
351
352 //====================================================================================================================================================
353
354 Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
355   
356   const Int_t nMuons = muons->GetEntriesFast();
357   if (nMuons<2 || nMuons>AliMFTConstants::fNMaxMuonsForPCA) {
358     printf("W-AliMFTAnalysisTools::CalculatePCA: number of muons not valid\n");
359     return kFALSE;
360   }
361
362   Double_t fXPointOfClosestApproach=0, fYPointOfClosestApproach=0, fZPointOfClosestApproach=0;
363
364   AliAODTrack *muon[AliMFTConstants::fNMaxMuonsForPCA]        = {0};
365   AliMUONTrackParam *param[AliMFTConstants::fNMaxMuonsForPCA] = {0};
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);
371     if (TMath::Abs(muon[iMu]->PzAtDCA())<1.e-6) {
372       for(Int_t i=0;i<iMu;i++) delete param[i];
373       return kFALSE;
374     }
375     param[iMu] = new AliMUONTrackParam();
376     param[iMu] -> SetNonBendingCoor(muon[iMu]->XAtDCA());
377     param[iMu] -> SetBendingCoor(muon[iMu]->YAtDCA());
378     param[iMu] -> SetZ(AliMFTConstants::fZEvalKinem);
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))) );
382   }
383   
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.;
388   
389   Double_t r[3]={0}, z[3]={startPoint-step, startPoint, startPoint+step};
390   
391   TVector3 **points = new TVector3*[AliMFTConstants::fNMaxMuonsForPCA];
392   
393   for (Int_t i=0; i<3; i++) {
394     for (Int_t iMu=0; iMu<nMuons; iMu++) {
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       // }
399       AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
400       points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
401     }
402     r[i] = GetDistanceBetweenPoints(points,nMuons);
403     for (Int_t iMu=0; iMu<nMuons; iMu++) delete points[iMu]; 
404   }
405   
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");
412     for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
413     delete points;
414     return kFALSE;
415   }     
416   
417   while (TMath::Abs(researchDirection)>0.5) {
418       
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");
431       for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
432       delete points;
433       return kFALSE;
434     }
435     
436     for (Int_t i=0; i<3; i++) {
437       for (Int_t iMu=0; iMu<nMuons; iMu++) {
438         AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
439         points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
440       }      
441       r[i] = GetDistanceBetweenPoints(points,nMuons);
442       for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu];
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
447     
448   }
449   
450   // now we now that the minimum is between z[0] and z[2] and we search for it
451   
452   Int_t nSteps = 0;
453
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++) {
460         AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
461         points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
462       }      
463       r[i] = GetDistanceBetweenPoints(points,nMuons);
464       for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu];
465     }
466     //    printf("Step #%d : %f  %f  %f\n",nSteps,r[0],r[1],r[2]);
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;
470     nSteps++;
471   }
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
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
498   
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);
510   
511   // Evaluating the PCA quality of the N-muon
512   
513   Double_t sum=0.,squareSum=0.;
514   for (Int_t iMu=0; iMu<nMuons; iMu++) {
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     }
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.;
527   
528   for(Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
529   delete points;
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()) +
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()) );
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
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
755 Bool_t AliMFTAnalysisTools::IsTrackInjected(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC) {
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;
766
767 }
768
769 //====================================================================================================================================================
770
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 //====================================================================================================================================================
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 //====================================================================================================================================================