Memory leak fixed in AliMFTAnalysisTools
[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
35 #include "AliMFTAnalysisTools.h"
36
37 ClassImp(AliMFTAnalysisTools)
38
39 //====================================================================================================================================================
40
41 Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2]) {
42
43   if (!(muon->Pz()!=0)) return kFALSE;
44
45   AliMUONTrackParam *param = new AliMUONTrackParam();
46
47   param -> SetNonBendingCoor(muon->XAtDCA());
48   param -> SetBendingCoor(muon->YAtDCA());
49   param -> SetZ(AliMFTConstants::fZEvalKinem);
50   param -> SetNonBendingSlope(muon->Px()/muon->Pz());
51   param -> SetBendingSlope(muon->Py()/muon->Pz());
52   param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) );
53
54   AliMUONTrackExtrap::ExtrapToZ(param, z);
55   xy[0] = param->GetNonBendingCoor();
56   xy[1] = param->GetBendingCoor();
57
58   delete param;
59
60   return kTRUE;
61
62 }
63
64 //====================================================================================================================================================
65
66 Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2], TLorentzVector &kinem) {
67
68   if (!(muon->Pz()!=0)) return kFALSE;
69
70   AliMUONTrackParam *param = new AliMUONTrackParam();
71
72   param -> SetNonBendingCoor(muon->XAtDCA());
73   param -> SetBendingCoor(muon->YAtDCA());
74   param -> SetZ(AliMFTConstants::fZEvalKinem);
75   param -> SetNonBendingSlope(muon->Px()/muon->Pz());
76   param -> SetBendingSlope(muon->Py()/muon->Pz());
77   param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) );
78
79   AliMUONTrackExtrap::ExtrapToZ(param, z);
80   xy[0] = param->GetNonBendingCoor();
81   xy[1] = param->GetBendingCoor();
82
83   Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
84   Double_t energy = TMath::Sqrt(massMu*massMu + param->Px()*param->Px() + param->Py()*param->Py() + param->Pz()*param->Pz());
85   
86   kinem.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
87
88   delete param;
89
90   return kTRUE;
91
92 }
93
94 //====================================================================================================================================================
95
96 Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2], TLorentzVector &kinem, TMatrixD &cov) {
97
98   if (!(muon->Pz()!=0)) return kFALSE;
99
100   AliMUONTrackParam *param = new AliMUONTrackParam();
101
102   param -> SetNonBendingCoor(muon->XAtDCA());
103   param -> SetBendingCoor(muon->YAtDCA());
104   param -> SetZ(AliMFTConstants::fZEvalKinem);
105   param -> SetNonBendingSlope(muon->Px()/muon->Pz());
106   param -> SetBendingSlope(muon->Py()/muon->Pz());
107   param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) );
108
109   param -> SetCovariances(ConvertCovMatrixAOD2MUON(muon));
110
111   AliMUONTrackExtrap::ExtrapToZCov(param, z);
112   xy[0] = param->GetNonBendingCoor();
113   xy[1] = param->GetBendingCoor();
114
115   Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
116   Double_t energy = TMath::Sqrt(massMu*massMu + param->Px()*param->Px() + param->Py()*param->Py() + param->Pz()*param->Pz());
117   
118   kinem.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
119
120   cov = param->GetCovariances();
121
122   delete param;
123
124   return kTRUE;
125
126 }
127
128 //====================================================================================================================================================
129
130 Double_t AliMFTAnalysisTools::GetAODMuonOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv) {
131
132   Double_t xy[2] = {0};
133   ExtrapAODMuonToZ(muon, zv, xy);
134   
135   return TMath::Sqrt((xv-xy[0])*(xv-xy[0]) + (yv-xy[1])*(yv-xy[1]));
136
137 }
138
139 //====================================================================================================================================================
140
141 Double_t AliMFTAnalysisTools::GetAODMuonWeightedOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv) {
142
143   Double_t xy[2] = {0};
144   TLorentzVector kinem(0,0,0,0);
145   TMatrixD cov(5,5);
146
147   ExtrapAODMuonToZ(muon, zv, xy, kinem, cov);
148
149   TMatrixD covCoordinates(2,2);
150   covCoordinates(0,0) = cov(0,0);
151   covCoordinates(0,1) = cov(0,2);
152   covCoordinates(1,0) = cov(2,0);
153   covCoordinates(1,1) = cov(2,2);
154   
155   TMatrixD covCoordinatesInverse = covCoordinates.Invert();
156
157   Double_t dX = xy[0] - xv;
158   Double_t dY = xy[1] - yv;
159   
160   Double_t weightedOffset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) + 
161                                              dY*dY*covCoordinatesInverse(1,1) + 
162                                              2.*dX*dY*covCoordinatesInverse(0,1)));
163
164   return weightedOffset;
165
166 }
167
168 //====================================================================================================================================================
169
170 Double_t AliMFTAnalysisTools::GetPseudoProperDecayTimeXY(Double_t xVtx,  Double_t yVtx, 
171                                                          Double_t xDimu, Double_t yDimu, 
172                                                          Double_t mDimu, Double_t ptDimu) {
173   
174   // pseudo-proper decay time of a particle produced in the primary vertex and decaying into a dimuon (+ X)
175   // evaluated using the transverse degree of freedom of the decay topology 
176
177   if (ptDimu != 0) {
178     Double_t decayLengthXY = TMath::Sqrt((xVtx-xDimu)*(xVtx-xDimu)+(yVtx-yDimu)*(yVtx-yDimu));
179     return (decayLengthXY * mDimu/ptDimu)/TMath::Ccgs()*1E12;   // in ps
180   }
181   
182   return -99999999;
183   
184 }
185
186 //====================================================================================================================================================
187
188 Double_t AliMFTAnalysisTools::GetPseudoProperDecayTimeZ(Double_t zVtx, 
189                                                         Double_t zDimu, 
190                                                         Double_t mDimu, Double_t pzDimu) {
191
192   // pseudo-proper decay time of a particle produced in the primary vertex and decaying into a dimuon (+ X)
193   // evaluated using the longitudinal degree of freedom of the decay topology 
194
195   if (pzDimu != 0) {
196     Double_t decayLengthZ = zDimu - zVtx;
197     return (decayLengthZ * mDimu/pzDimu)/TMath::Ccgs()*1E12;  // in ps
198   }
199
200   return -99999999;
201
202 }
203
204 //====================================================================================================================================================
205
206 Bool_t AliMFTAnalysisTools::CalculatePCA(AliAODDimuon *dimuon, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
207
208   TObjArray *muons = new TObjArray();
209   muons -> Add(dimuon->GetMu(0));
210   muons -> Add(dimuon->GetMu(1));
211   
212   Bool_t result = CalculatePCA(muons, pca, pcaQuality, kinem);
213   delete muons;
214   return result;
215
216 }
217
218 //====================================================================================================================================================
219
220 Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
221   
222   const Int_t nMuons = muons->GetEntriesFast();
223   if (nMuons<2 || nMuons>AliMFTConstants::fNMaxMuonsForPCA) {
224     printf("W-AliMFTAnalysisTools::CalculatePCA: number of muons not valid\n");
225     return kFALSE;
226   }
227   
228   Double_t fXPointOfClosestApproach=0, fYPointOfClosestApproach=0, fZPointOfClosestApproach=0;
229
230   AliAODTrack *muon[AliMFTConstants::fNMaxMuonsForPCA]        = {0};
231   AliMUONTrackParam *param[AliMFTConstants::fNMaxMuonsForPCA] = {0};
232
233   // Finding AliMUONTrackParam objects for each muon
234   
235   for (Int_t iMu=0; iMu<nMuons; iMu++) {
236     muon[iMu] = (AliAODTrack*) muons->At(iMu);
237     if (TMath::Abs(muon[iMu]->Pz())<1.e-6) {
238       for(Int_t i=0;i<iMu;i++) delete param[i];
239       return kFALSE;
240     }
241     param[iMu] = new AliMUONTrackParam();
242     param[iMu] -> SetNonBendingCoor(muon[iMu]->XAtDCA());
243     param[iMu] -> SetBendingCoor(muon[iMu]->YAtDCA());
244     param[iMu] -> SetZ(0.);
245     param[iMu] -> SetNonBendingSlope(muon[iMu]->Px()/muon[iMu]->Pz());
246     param[iMu] -> SetBendingSlope(muon[iMu]->Py()/muon[iMu]->Pz());
247     param[iMu] -> SetInverseBendingMomentum( muon[iMu]->Charge() * (1./muon[iMu]->Pz()) / (TMath::Sqrt(1+TMath::Power(muon[iMu]->Py()/muon[iMu]->Pz(),2))) );
248   }
249   
250   // here we want to understand in which direction we have to search the minimum...
251   
252   Double_t step = 1.;  // initial step, in cm
253   Double_t startPoint = 0.;
254   
255   Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step};
256   
257   TVector3 **points = new TVector3*[AliMFTConstants::fNMaxMuonsForPCA];
258   
259   for (Int_t i=0; i<3; i++) {
260     for (Int_t iMu=0; iMu<nMuons; iMu++) {
261       AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
262       points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
263                 }
264     r[i] = GetDistanceBetweenPoints(points,nMuons);
265     for (Int_t iMu=0; iMu<nMuons; iMu++) delete points[iMu]; 
266   }
267   
268   Int_t researchDirection = 0;
269   
270   if      (r[0]>r[1] && r[1]>r[2]) researchDirection = +1;   // towards z positive
271   else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1;   // towards z negative
272   else if (r[0]<r[1] && r[1]>r[2]) {
273     printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima)\n");
274     for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
275     delete points;
276     return kFALSE;
277   }     
278   
279   while (TMath::Abs(researchDirection)>0.5) {
280       
281     if (researchDirection>0) {
282       z[0] = z[1];
283       z[1] = z[2];
284       z[2] = z[1]+researchDirection*step;
285     }
286     else {
287       z[2] = z[1];
288       z[1] = z[0];
289       z[0] = z[1]+researchDirection*step;
290     }
291     if (TMath::Abs(z[0])>900.) {
292       printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima in the fiducial region)\n");
293       for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
294       delete points;
295       return kFALSE;
296     }
297     
298     for (Int_t i=0; i<3; i++) {
299       for (Int_t iMu=0; iMu<nMuons; iMu++) {
300         AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
301         points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
302       }      
303       r[i] = GetDistanceBetweenPoints(points,nMuons);
304       for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu];
305     }
306     researchDirection=0;
307     if      (r[0]>r[1] && r[1]>r[2]) researchDirection = +1;   // towards z positive
308     else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1;   // towards z negative
309     
310   }
311   
312   // now we now that the minimum is between z[0] and z[2] and we search for it
313   
314   step *= 0.5;
315   while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) {
316     z[0] = z[1]-step;
317     z[2] = z[1]+step;
318     for (Int_t i=0; i<3; i++) {
319       for (Int_t iMu=0; iMu<nMuons; iMu++) {
320         AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
321         points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
322       }      
323       r[i] = GetDistanceBetweenPoints(points,nMuons);
324       for (Int_t iMu=0;iMu<nMuons;iMu++)        delete points[iMu];
325     }
326     if      (r[0]<r[1]) z[1] = z[0];
327     else if (r[2]<r[1]) z[1] = z[2];
328     else step *= 0.5;
329   }
330   
331   // Once z of minimum is found, we evaluate the x and y coordinates by averaging over the contributing tracks
332   
333   fZPointOfClosestApproach = z[1];
334   fXPointOfClosestApproach = 0.;
335   fYPointOfClosestApproach = 0.;
336   for (Int_t iMu=0; iMu<nMuons; iMu++) {
337     AliMUONTrackExtrap::ExtrapToZ(param[iMu], fZPointOfClosestApproach);
338     fXPointOfClosestApproach += param[iMu]->GetNonBendingCoor();
339     fYPointOfClosestApproach += param[iMu]->GetBendingCoor();
340   }
341   fXPointOfClosestApproach /= Double_t(nMuons);
342   fYPointOfClosestApproach /= Double_t(nMuons);
343   
344   pca[0] = fXPointOfClosestApproach;
345   pca[1] = fYPointOfClosestApproach;
346   pca[2] = fZPointOfClosestApproach;
347   
348   // Evaluating the kinematics of the N-muon
349   
350   Double_t pTot[3] = {0};
351   Double_t ene = 0.;
352   Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
353   for (Int_t iMu=0; iMu<nMuons; iMu++) {
354     pTot[0] += param[iMu]->Px();
355     pTot[1] += param[iMu]->Py();
356     pTot[2] += param[iMu]->Pz();
357     ene += TMath::Sqrt(massMu*massMu + param[iMu]->Px()*param[iMu]->Px() + param[iMu]->Py()*param[iMu]->Py() + param[iMu]->Pz()*param[iMu]->Pz());
358   }
359   
360   kinem.SetPxPyPzE(pTot[0], pTot[1], pTot[2], ene);
361   
362   // Evaluating the PCA quality of the N-muon
363   
364   Double_t sum=0.,squareSum=0.;
365   for (Int_t iMu=0; iMu<nMuons; iMu++) {
366     Double_t wOffset = AliMFTAnalysisTools::GetAODMuonWeightedOffset(muon[iMu],fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach);
367     Double_t f = TMath::Exp(-0.5 * wOffset);
368     sum += f;
369     squareSum += f*f;
370   }
371   if (sum > 0.) pcaQuality =  (sum-squareSum/sum) / (nMuons-1);
372   else pcaQuality = 0.;
373   
374   for(Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
375   delete points;
376   return kTRUE;
377   
378 }
379
380 //=========================================================================================================================
381
382 Double_t AliMFTAnalysisTools::GetDistanceBetweenPoints(TVector3 **points, Int_t nPoints) {
383   
384   if (nPoints>AliMFTConstants::fNMaxMuonsForPCA) {
385     printf("W-AliMFTAnalysisTools::GetDistanceBetweenPoints: number of points not valid\n");
386     return 1.e9;
387   }
388   
389   if (nPoints<2) return 0.;
390   if (nPoints<3) return TMath::Sqrt( (points[0]->X()-points[1]->X()) * (points[0]->X()-points[1]->X()) +
391                                      (points[0]->Y()-points[1]->Y()) * (points[0]->Y()-points[1]->Y()) +
392                                      (points[0]->Z()-points[1]->Z()) * (points[0]->Z()-points[1]->Z()) );
393
394   const Int_t nEdgesMax = ((AliMFTConstants::fNMaxMuonsForPCA) * (AliMFTConstants::fNMaxMuonsForPCA - 1)) / 2;
395   
396   Int_t startID[nEdgesMax]       = {0};
397   Int_t stopID[nEdgesMax]        = {0};
398   Double_t edgeLength[nEdgesMax] = {0};
399
400   Bool_t pointStatus[AliMFTConstants::fNMaxMuonsForPCA] = {0};
401   
402   Int_t nEdges=0;
403   for (Int_t i=0; i<nPoints-1; i++) {
404     for (Int_t j=i+1; j<nPoints; j++) {
405       edgeLength[nEdges] = TMath::Sqrt( (points[i]->X()-points[j]->X()) * (points[i]->X()-points[j]->X()) +
406                                         (points[i]->Y()-points[j]->Y()) * (points[i]->Y()-points[j]->Y()) +
407                                         (points[i]->Z()-points[j]->Z()) * (points[i]->Z()-points[j]->Z()) );
408       stopID[nEdges]  = i;
409       startID[nEdges] = j;
410       nEdges++;
411     }
412   }
413  
414   // Order Edges
415
416   Double_t min = 0;
417   Int_t   iMin = 0;
418
419   for (Int_t iEdge=0; iEdge<nEdges-1; iEdge++) {
420     min  = edgeLength[iEdge];
421     iMin = iEdge;
422     for (Int_t j=iEdge+1; j<nEdges; j++) {
423       if (edgeLength[j]<min) {
424         min  = edgeLength[j];
425         iMin = j;
426       }
427     }
428     
429     if (iMin != iEdge) {
430
431       Double_t edgeLengthMin = edgeLength[iMin];
432       Int_t startIDmin = startID[iMin];
433       Int_t stopIDmin  = stopID[iMin];
434       
435       edgeLength[iMin] = edgeLength[iEdge];
436       startID[iMin]    = startID[iEdge];
437       stopID[iMin]     = stopID[iEdge];
438
439       edgeLength[iEdge] = edgeLengthMin;
440       startID[iEdge]    = startIDmin;
441       stopID[iEdge]     = stopIDmin;
442
443     }
444     
445   }
446   
447   // Connect
448
449   Double_t length = 0.;
450   for (Int_t i=0; i<nEdges; i++) {
451     if (!(pointStatus[startID[i]] && pointStatus[stopID[i]])) {
452       pointStatus[startID[i]] = kTRUE;
453       pointStatus[stopID[i]]  = kTRUE;
454       length += edgeLength[i];
455     }
456   }
457   
458   return length;
459   
460 }
461
462 //====================================================================================================================================================
463
464 void AliMFTAnalysisTools::ConvertCovMatrixMUON2AOD(const TMatrixD& covMUON, Double_t covAOD[21]) {
465
466   // Converts the cov matrix from the MUON format (TMatrixD) to the AOD one (Double_t[21])
467   // 
468   // Cov(x,x)       ... :   cv[0]
469   // Cov(x,slopeX)  ... :   cv[1]  cv[2]
470   // Cov(x,y)       ... :   cv[3]  cv[4]  cv[5]
471   // Cov(x,slopeY)  ... :   cv[6]  cv[7]  cv[8]  cv[9]
472   // Cov(x,invP_yz) ... :   cv[10] cv[11] cv[12] cv[13] cv[14]
473   // not-used       ... :   cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
474
475   covAOD[0]  = covMUON(0,0);
476
477   covAOD[1]  = covMUON(1,0);
478   covAOD[2]  = covMUON(1,1);
479
480   covAOD[3]  = covMUON(2,0);  
481   covAOD[4]  = covMUON(2,1);  
482   covAOD[5]  = covMUON(2,2);  
483
484   covAOD[6]  = covMUON(3,0);  
485   covAOD[7]  = covMUON(3,1);  
486   covAOD[8]  = covMUON(3,2);  
487   covAOD[9]  = covMUON(3,3);  
488
489   covAOD[10] = covMUON(4,0);  
490   covAOD[11] = covMUON(4,1);  
491   covAOD[12] = covMUON(4,2);  
492   covAOD[13] = covMUON(4,3);  
493   covAOD[14] = covMUON(4,4);  
494
495   covAOD[15] = 0;  
496   covAOD[16] = 0;  
497   covAOD[17] = 0;  
498   covAOD[18] = 0;  
499   covAOD[19] = 0;  
500   covAOD[20] = 0;  
501
502 }
503
504 //====================================================================================================================================================
505
506 const TMatrixD AliMFTAnalysisTools::ConvertCovMatrixAOD2MUON(AliAODTrack *muon) {
507
508   Double_t covAOD[21] = {0};
509   muon -> GetCovarianceXYZPxPyPz(covAOD);
510
511   TMatrixD covMUON(5,5);
512
513   covMUON(0,0) = covAOD[0];
514                                 
515   covMUON(1,0) = covAOD[1];
516   covMUON(1,1) = covAOD[2];
517                                 
518   covMUON(2,0) = covAOD[3];  
519   covMUON(2,1) = covAOD[4];  
520   covMUON(2,2) = covAOD[5];  
521                                 
522   covMUON(3,0) = covAOD[6];  
523   covMUON(3,1) = covAOD[7];  
524   covMUON(3,2) = covAOD[8];  
525   covMUON(3,3) = covAOD[9];  
526                                 
527   covMUON(4,0) = covAOD[10];  
528   covMUON(4,1) = covAOD[11];  
529   covMUON(4,2) = covAOD[12];  
530   covMUON(4,3) = covAOD[13];  
531   covMUON(4,4) = covAOD[14]; 
532
533   return covMUON;
534
535 }
536
537 //====================================================================================================================================================
538