AliRoot Core  a565103 (a565103)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliMFTAnalysisTools.cxx
Go to the documentation of this file.
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 
40 
41 //====================================================================================================================================================
42 
43 Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2]) {
44 
45  if (!(muon->PzAtDCA()!=0)) return kFALSE;
46 
47  AliMUONTrackParam *param = new AliMUONTrackParam();
48 
49  param -> SetNonBendingCoor(muon->XAtDCA());
50  param -> SetBendingCoor(muon->YAtDCA());
51  param -> SetZ(AliMFTConstants::fZEvalKinem);
52  param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
53  param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
54  param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon->PyAtDCA()/muon->PzAtDCA(),2))) );
55 
57  xy[0] = param->GetNonBendingCoor();
58  xy[1] = param->GetBendingCoor();
59 
60  delete param;
61 
62  return kTRUE;
63 
64 }
65 
66 //====================================================================================================================================================
67 
68 Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2], TLorentzVector &kinem) {
69 
70  if (!(muon->PzAtDCA()!=0)) return kFALSE;
71 
72  AliMUONTrackParam *param = new AliMUONTrackParam();
73 
74  param -> SetNonBendingCoor(muon->XAtDCA());
75  param -> SetBendingCoor(muon->YAtDCA());
76  param -> SetZ(AliMFTConstants::fZEvalKinem);
77  param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
78  param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
79  param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon->PyAtDCA()/muon->PzAtDCA(),2))) );
80 
82  xy[0] = param->GetNonBendingCoor();
83  xy[1] = param->GetBendingCoor();
84 
85  Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
86  Double_t energy = TMath::Sqrt(massMu*massMu + param->Px()*param->Px() + param->Py()*param->Py() + param->Pz()*param->Pz());
87 
88  kinem.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
89 
90  delete param;
91 
92  return kTRUE;
93 
94 }
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  if (!(muon->PzAtDCA()!=0)) return kFALSE;
103 
104  AliMUONTrackParam *param = new AliMUONTrackParam();
105 
106  param -> SetNonBendingCoor(muon->XAtDCA());
107  param -> SetBendingCoor(muon->YAtDCA());
108  param -> SetZ(AliMFTConstants::fZEvalKinem);
109  param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
110  param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
111  param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon->PyAtDCA()/muon->PzAtDCA(),2))) );
112 
113  param -> SetCovariances(ConvertCovMatrixAOD2MUON(muon));
114 
116  xy[0] = param->GetNonBendingCoor();
117  xy[1] = param->GetBendingCoor();
118 
119  Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
120  Double_t energy = TMath::Sqrt(massMu*massMu + param->Px()*param->Px() + param->Py()*param->Py() + param->Pz()*param->Pz());
121 
122  kinem.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
123 
124  cov = param->GetCovariances();
125 
126  delete param;
127 
128  return kTRUE;
129 
130 }
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  AliMUONTrackParam *param = new AliMUONTrackParam();
142  param -> SetNonBendingCoor(muon->XAtDCA());
143  param -> SetBendingCoor(muon->YAtDCA());
144  param -> SetZ(AliMFTConstants::fZEvalKinem);
145  param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
146  param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
147  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  Double_t r[3]={0}, z[3]={startPoint-step, startPoint, startPoint+step};
155 
156  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  for (Int_t i=0; i<3; i++) {
159  AliMUONTrackExtrap::ExtrapToZ(param, z[i]);
160  points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]);
161  points[1] = new TVector3(xy[0],xy[1],z[i]);
162  r[i] = GetDistanceBetweenPoints(points,2);
163  for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu];
164  }
165 
166  Int_t researchDirection = 0;
167 
168  if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
169  else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative
170  else if (r[0]<r[1] && r[1]>r[2]) {
171  printf("E-AliMFTAnalysisTools::ExtrapAODMuonToXY: Point of closest approach cannot be found (no minima)\n");
172  delete param;
173  delete points;
174  return kFALSE;
175  }
176 
177  while (TMath::Abs(researchDirection)>0.5) {
178 
179  if (researchDirection>0) {
180  z[0] = z[1];
181  z[1] = z[2];
182  z[2] = z[1]+researchDirection*step;
183  }
184  else {
185  z[2] = z[1];
186  z[1] = z[0];
187  z[0] = z[1]+researchDirection*step;
188  }
189  if (TMath::Abs(z[0])>900.) {
190  printf("E-AliMFTAnalysisTools::ExtrapAODMuonToXY: Point of closest approach cannot be found (no minima in the fiducial region)\n");
191  delete param;
192  delete points;
193  return kFALSE;
194  }
195 
196  for (Int_t i=0; i<3; i++) {
197  AliMUONTrackExtrap::ExtrapToZ(param, z[i]);
198  points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]);
199  points[1] = new TVector3(xy[0],xy[1],z[i]);
200  r[i] = GetDistanceBetweenPoints(points,2);
201  for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu];
202  }
203 
204  researchDirection=0;
205  if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
206  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;
214  z[0] = z[1]-step;
215  z[2] = z[1]+step;
216  for (Int_t i=0; i<3; i++) {
217  AliMUONTrackExtrap::ExtrapToZ(param, z[i]);
218  points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]);
219  points[1] = new TVector3(xy[0],xy[1],z[i]);
220  r[i] = GetDistanceBetweenPoints(points,2);
221  for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu];
222  }
223  if (r[0]<r[1]) z[1] = z[0];
224  else if (r[2]<r[1]) z[1] = z[2];
225  else step *= 0.5;
226  }
227 
228  zFinal = z[1];
229 
230  Double_t xyMuon[2] = {0};
231  ExtrapAODMuonToZ(muon, zFinal, xyMuon, kinem, cov);
232 
233  return kTRUE;
234 
235 }
236 
237 //====================================================================================================================================================
238 
239 Bool_t AliMFTAnalysisTools::GetAODMuonOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
240 
241  Double_t xy[2] = {0};
242  ExtrapAODMuonToZ(muon, zv, xy);
243 
244  offset = TMath::Sqrt((xv-xy[0])*(xv-xy[0]) + (yv-xy[1])*(yv-xy[1]));
245 
246  return kTRUE;
247 
248 }
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  Double_t xy[2] = {0};
258  ExtrapAODMuonToZ(muon, zv, xy);
259 
260  xy[0] = gRandom->Gaus(xy[0], smearOffsetX);
261  xy[1] = gRandom->Gaus(xy[1], smearOffsetY);
262 
263  offset = TMath::Sqrt((xv-xy[0])*(xv-xy[0]) + (yv-xy[1])*(yv-xy[1]));
264 
265  return kTRUE;
266 
267 }
268 
269 //====================================================================================================================================================
270 
271 Bool_t AliMFTAnalysisTools::GetAODMuonOffsetZ(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
272 
273  Double_t xy[2] = {xv, yv};
274  Double_t zFinal = 0;
275  TLorentzVector kinem(0,0,0,0);
276  TMatrixD cov(5,5);
277 
278  ExtrapAODMuonToXY(muon, xy, zFinal, kinem, cov);
279 
280  offset = TMath::Abs(zFinal - zv);
281 
282  return kTRUE;
283 
284 }
285 
286 //====================================================================================================================================================
287 
288 Bool_t AliMFTAnalysisTools::GetAODMuonWeightedOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
289 
290  Double_t xy[2] = {0};
291  TLorentzVector kinem(0,0,0,0);
292  TMatrixD cov(5,5);
293 
294  ExtrapAODMuonToZ(muon, zv, xy, kinem, cov);
295 
296  TMatrixD covCoordinates(2,2);
297  covCoordinates(0,0) = cov(0,0);
298  covCoordinates(0,1) = cov(0,2);
299  covCoordinates(1,0) = cov(2,0);
300  covCoordinates(1,1) = cov(2,2);
301 
302  if (covCoordinates.Determinant() < covCoordinates.GetTol()) return kFALSE;
303 
304  if (TDecompLU::InvertLU(covCoordinates,covCoordinates.GetTol(),0)) {
305 
306  TMatrixD covCoordinatesInverse = covCoordinates;
307  Double_t dX = xy[0] - xv;
308  Double_t dY = xy[1] - yv;
309 
310  offset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) +
311  dY*dY*covCoordinatesInverse(1,1) +
312  2.*dX*dY*covCoordinatesInverse(0,1)));
313 
314  return kTRUE;
315 
316  }
317 
318  return kFALSE;
319 
320 }
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  if (ptDimu != 0) {
332  Double_t decayLengthXY = TMath::Sqrt((xVtx-xDimu)*(xVtx-xDimu)+(yVtx-yDimu)*(yVtx-yDimu));
333  return (decayLengthXY * mDimu/ptDimu)/TMath::Ccgs()*1E12; // in ps
334  }
335 
336  return -99999999;
337 
338 }
339 
340 //====================================================================================================================================================
341 
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  if (pzDimu != 0) {
350  Double_t decayLengthZ = zDimu - zVtx;
351  return (decayLengthZ * mDimu/pzDimu)/TMath::Ccgs()*1E12; // in ps
352  }
353 
354  return -99999999;
355 
356 }
357 
358 //====================================================================================================================================================
359 
360 Bool_t AliMFTAnalysisTools::CalculatePCA(AliAODDimuon *dimuon, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
361 
362  TObjArray *muons = new TObjArray();
363  muons -> Add(dimuon->GetMu(0));
364  muons -> Add(dimuon->GetMu(1));
365 
366  Bool_t result = CalculatePCA(muons, pca, pcaQuality, kinem);
367  delete muons;
368  return result;
369 
370 }
371 
372 //====================================================================================================================================================
373 
374 Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
375 
376  const Int_t nMuons = muons->GetEntriesFast();
377  if (nMuons<2 || nMuons>AliMFTConstants::fNMaxMuonsForPCA) {
378  printf("W-AliMFTAnalysisTools::CalculatePCA: number of muons not valid\n");
379  return kFALSE;
380  }
381 
382  Double_t fXPointOfClosestApproach=0, fYPointOfClosestApproach=0, fZPointOfClosestApproach=0;
383 
384  AliAODTrack *muon[AliMFTConstants::fNMaxMuonsForPCA] = {0};
386 
387  // Finding AliMUONTrackParam objects for each muon
388 
389  for (Int_t iMu=0; iMu<nMuons; iMu++) {
390  muon[iMu] = (AliAODTrack*) muons->At(iMu);
391  if (TMath::Abs(muon[iMu]->PzAtDCA())<1.e-6) {
392  for(Int_t i=0;i<iMu;i++) delete param[i];
393  return kFALSE;
394  }
395  param[iMu] = new AliMUONTrackParam();
396  param[iMu] -> SetNonBendingCoor(muon[iMu]->XAtDCA());
397  param[iMu] -> SetBendingCoor(muon[iMu]->YAtDCA());
398  param[iMu] -> SetZ(AliMFTConstants::fZEvalKinem);
399  param[iMu] -> SetNonBendingSlope(muon[iMu]->PxAtDCA()/muon[iMu]->PzAtDCA());
400  param[iMu] -> SetBendingSlope(muon[iMu]->PyAtDCA()/muon[iMu]->PzAtDCA());
401  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  Double_t r[3]={0}, z[3]={startPoint-step, startPoint, startPoint+step};
410 
411  TVector3 **points = new TVector3*[AliMFTConstants::fNMaxMuonsForPCA];
412 
413  for (Int_t i=0; i<3; i++) {
414  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  AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
420  points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
421  }
422  r[i] = GetDistanceBetweenPoints(points,nMuons);
423  for (Int_t iMu=0; iMu<nMuons; iMu++) delete points[iMu];
424  }
425 
426  Int_t researchDirection = 0;
427 
428  if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
429  else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative
430  else if (r[0]<r[1] && r[1]>r[2]) {
431  printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima)\n");
432  for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
433  delete points;
434  return kFALSE;
435  }
436 
437  while (TMath::Abs(researchDirection)>0.5) {
438 
439  if (researchDirection>0) {
440  z[0] = z[1];
441  z[1] = z[2];
442  z[2] = z[1]+researchDirection*step;
443  }
444  else {
445  z[2] = z[1];
446  z[1] = z[0];
447  z[0] = z[1]+researchDirection*step;
448  }
449  if (TMath::Abs(z[0])>900.) {
450  printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima in the fiducial region)\n");
451  for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
452  delete points;
453  return kFALSE;
454  }
455 
456  for (Int_t i=0; i<3; i++) {
457  for (Int_t iMu=0; iMu<nMuons; iMu++) {
458  AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
459  points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
460  }
461  r[i] = GetDistanceBetweenPoints(points,nMuons);
462  for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu];
463  }
464  researchDirection=0;
465  if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
466  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;
476  z[0] = z[1]-step;
477  z[2] = z[1]+step;
478  for (Int_t i=0; i<3; i++) {
479  for (Int_t iMu=0; iMu<nMuons; iMu++) {
480  AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
481  points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
482  }
483  r[i] = GetDistanceBetweenPoints(points,nMuons);
484  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  if (r[0]<r[1]) z[1] = z[0];
488  else if (r[2]<r[1]) z[1] = z[2];
489  else step *= 0.5;
490  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  for (Int_t iMu=0; iMu<nMuons; iMu++) {
506  AliMUONTrackExtrap::ExtrapToZ(param[iMu], fZPointOfClosestApproach);
507  fXPointOfClosestApproach += param[iMu]->GetNonBendingCoor();
508  fYPointOfClosestApproach += param[iMu]->GetBendingCoor();
509  }
510  fXPointOfClosestApproach /= Double_t(nMuons);
511  fYPointOfClosestApproach /= Double_t(nMuons);
512 
513  pca[0] = fXPointOfClosestApproach;
514  pca[1] = fYPointOfClosestApproach;
515  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  Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
522  for (Int_t iMu=0; iMu<nMuons; iMu++) {
523  pTot[0] += param[iMu]->Px();
524  pTot[1] += param[iMu]->Py();
525  pTot[2] += param[iMu]->Pz();
526  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  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  for (Int_t iMu=0; iMu<nMuons; iMu++) {
535  Double_t wOffset = 0;
536  if (!GetAODMuonWeightedOffset(muon[iMu],fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach, wOffset)) {
537  for(Int_t jMu=0;jMu<nMuons;jMu++) delete param[jMu];
538  delete points;
539  return kFALSE;
540  }
541  Double_t f = TMath::Exp(-0.5 * wOffset);
542  sum += f;
543  squareSum += f*f;
544  }
545  if (sum > 0.) pcaQuality = (sum-squareSum/sum) / (nMuons-1);
546  else pcaQuality = 0.;
547 
548  for(Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
549  delete points;
550  return kTRUE;
551 
552 }
553 
554 //=========================================================================================================================
555 
556 Double_t AliMFTAnalysisTools::GetDistanceBetweenPoints(TVector3 **points, Int_t nPoints) {
557 
558  if (nPoints>AliMFTConstants::fNMaxMuonsForPCA) {
559  printf("W-AliMFTAnalysisTools::GetDistanceBetweenPoints: number of points not valid\n");
560  return 1.e9;
561  }
562 
563  if (nPoints<2) return 0.;
564  if (nPoints<3) return TMath::Sqrt( (points[0]->X()-points[1]->X()) * (points[0]->X()-points[1]->X()) +
565  (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  Int_t startID[nEdgesMax] = {0};
571  Int_t stopID[nEdgesMax] = {0};
572  Double_t edgeLength[nEdgesMax] = {0};
573 
574  Bool_t pointStatus[AliMFTConstants::fNMaxMuonsForPCA] = {0};
575 
576  Int_t nEdges=0;
577  for (Int_t i=0; i<nPoints-1; i++) {
578  for (Int_t j=i+1; j<nPoints; j++) {
579  edgeLength[nEdges] = TMath::Sqrt( (points[i]->X()-points[j]->X()) * (points[i]->X()-points[j]->X()) +
580  (points[i]->Y()-points[j]->Y()) * (points[i]->Y()-points[j]->Y()) +
581  (points[i]->Z()-points[j]->Z()) * (points[i]->Z()-points[j]->Z()) );
582  stopID[nEdges] = i;
583  startID[nEdges] = j;
584  nEdges++;
585  }
586  }
587 
588  // Order Edges
589 
590  Double_t min = 0;
591  Int_t iMin = 0;
592 
593  for (Int_t iEdge=0; iEdge<nEdges-1; iEdge++) {
594  min = edgeLength[iEdge];
595  iMin = iEdge;
596  for (Int_t j=iEdge+1; j<nEdges; j++) {
597  if (edgeLength[j]<min) {
598  min = edgeLength[j];
599  iMin = j;
600  }
601  }
602 
603  if (iMin != iEdge) {
604 
605  Double_t edgeLengthMin = edgeLength[iMin];
606  Int_t startIDmin = startID[iMin];
607  Int_t stopIDmin = stopID[iMin];
608 
609  edgeLength[iMin] = edgeLength[iEdge];
610  startID[iMin] = startID[iEdge];
611  stopID[iMin] = stopID[iEdge];
612 
613  edgeLength[iEdge] = edgeLengthMin;
614  startID[iEdge] = startIDmin;
615  stopID[iEdge] = stopIDmin;
616 
617  }
618 
619  }
620 
621  // Connect
622 
623  Double_t length = 0.;
624  for (Int_t i=0; i<nEdges; i++) {
625  if (!(pointStatus[startID[i]] && pointStatus[stopID[i]])) {
626  pointStatus[startID[i]] = kTRUE;
627  pointStatus[stopID[i]] = kTRUE;
628  length += edgeLength[i];
629  }
630  }
631 
632  return length;
633 
634 }
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  covAOD[0] = covMUON(0,0);
650 
651  covAOD[1] = covMUON(1,0);
652  covAOD[2] = covMUON(1,1);
653 
654  covAOD[3] = covMUON(2,0);
655  covAOD[4] = covMUON(2,1);
656  covAOD[5] = covMUON(2,2);
657 
658  covAOD[6] = covMUON(3,0);
659  covAOD[7] = covMUON(3,1);
660  covAOD[8] = covMUON(3,2);
661  covAOD[9] = covMUON(3,3);
662 
663  covAOD[10] = covMUON(4,0);
664  covAOD[11] = covMUON(4,1);
665  covAOD[12] = covMUON(4,2);
666  covAOD[13] = covMUON(4,3);
667  covAOD[14] = covMUON(4,4);
668 
669  covAOD[15] = 0;
670  covAOD[16] = 0;
671  covAOD[17] = 0;
672  covAOD[18] = 0;
673  covAOD[19] = 0;
674  covAOD[20] = 0;
675 
676 }
677 
678 //====================================================================================================================================================
679 
681 
682  Double_t covAOD[21] = {0};
683  muon -> GetCovarianceXYZPxPyPz(covAOD);
684 
685  TMatrixD covMUON(5,5);
686 
687  covMUON(0,0) = covAOD[0];
688 
689  covMUON(1,0) = covAOD[1];
690  covMUON(1,1) = covAOD[2];
691 
692  covMUON(2,0) = covAOD[3];
693  covMUON(2,1) = covAOD[4];
694  covMUON(2,2) = covAOD[5];
695 
696  covMUON(3,0) = covAOD[6];
697  covMUON(3,1) = covAOD[7];
698  covMUON(3,2) = covAOD[8];
699  covMUON(3,3) = covAOD[9];
700 
701  covMUON(4,0) = covAOD[10];
702  covMUON(4,1) = covAOD[11];
703  covMUON(4,2) = covAOD[12];
704  covMUON(4,3) = covAOD[13];
705  covMUON(4,4) = covAOD[14];
706 
707  return covMUON;
708 
709 }
710 
711 //====================================================================================================================================================
712 
714 
715  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) if (IsWrongCluster(muon, iPlane)) return kFALSE;
716  return kTRUE;
717 
718 }
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  TList *genHeaders = header->GetCocktailHeaders();
728  Int_t nGenHeaders = genHeaders->GetEntries();
729 
730  for (Int_t i=0; i<nGenHeaders; i++){
731  AliGenEventHeader *gh = (AliGenEventHeader*) genHeaders->At(i);
732  TString genName = gh->GetName();
733  Int_t nPart = gh->NProduced();
734  if (label>=partCounter && label<(partCounter+nPart)) return genName;
735  partCounter += nPart;
736  }
737 
738  TString empty="";
739  return empty;
740 
741 }
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  Int_t label = TMath::Abs(track->GetLabel());
750  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  while (nameGen.IsWhitespace()) {
755  AliAODMCParticle *mcPart = (AliAODMCParticle*) arrayMC->At(label);
756  if (!mcPart) {
757  printf("AliMFTAnalysisTools::GetTrackPrimaryGenerator - BREAK: No valid AliAODMCParticle at label %i\n",label);
758  break;
759  }
760  Int_t motherLabel = mcPart->GetMother();
761  if (motherLabel < 0) {
762  printf("AliMFTAnalysisTools::GetTrackPrimaryGenerator - BREAK: Reached primary particle without valid mother\n");
763  break;
764  }
765  label = motherLabel;
766  nameGen = GetGenerator(label,header);
767  }
768 
769  return;
770 
771 }
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  TString nameGen;
780 
781  GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
782 
783  if (nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
784 
785  return kTRUE;
786 
787 }
788 
789 //====================================================================================================================================================
790 
791 Bool_t AliMFTAnalysisTools::TranslateMuon(AliAODTrack *muon, Double_t vtxInitial[3], Double_t vtxFinal[3]) {
792 
793  if (!(muon->PzAtDCA()!=0)) return kFALSE;
794 
795  AliMUONTrackParam *param = new AliMUONTrackParam();
796 
797  Double_t deltaVtx[3] = {0};
798  for (Int_t i=0; i<3; i++) deltaVtx[i] = vtxInitial[i] - vtxFinal[i];
799 
800  param -> SetNonBendingCoor(muon->XAtDCA());
801  param -> SetBendingCoor(muon->YAtDCA());
802  param -> SetZ(AliMFTConstants::fZEvalKinem);
803  param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
804  param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
805  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  AliMUONTrackExtrap::ExtrapToZ(param, deltaVtx[2]);
811  muon->SetXYAtDCA(param->GetNonBendingCoor() - deltaVtx[0], param->GetBendingCoor() - deltaVtx[1]);
812  muon->SetPxPyPzAtDCA(param->Px(), param->Py(), param->Pz());
813 
814  delete param;
815 
816  return kTRUE;
817 
818 }
819 
820 //====================================================================================================================================================
821 
822 Bool_t AliMFTAnalysisTools::TranslateMuonToOrigin(AliAODTrack *muon, Double_t vtx[3]) {
823 
824  Double_t origin[3] = {0,0,0};
825 
826  return TranslateMuon(muon, vtx, origin);
827 
828 }
829 
830 //====================================================================================================================================================
831 
832 Bool_t AliMFTAnalysisTools::IsPDGCharm(Int_t pdgCode) {
833 
834  pdgCode = TMath::Abs(pdgCode/100);
835  if (pdgCode>9) pdgCode /= 10;
836  if (pdgCode == 4 ) return kTRUE;
837  else return kFALSE;
838 
839 }
840 
841 //====================================================================================================================================================
842 
843 Bool_t AliMFTAnalysisTools::IsPDGBeauty(Int_t pdgCode) {
844 
845  pdgCode = TMath::Abs(pdgCode/100);
846  if (pdgCode>9) pdgCode /= 10;
847  if (pdgCode == 5) return kTRUE;
848  else return kFALSE;
849 
850 }
851 
852 //====================================================================================================================================================
853 
854 Bool_t AliMFTAnalysisTools::IsPDGResonance(Int_t pdgCode) {
855 
856  Int_t id = pdgCode%100000;
857  return (!((id-id%10)%110));
858 
859 }
860 
861 //====================================================================================================================================================
static Bool_t CalculatePCA(AliAODDimuon *dimuon, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem)
static const Double_t fPrecisionPointOfClosestApproach
precision (along z) for the research of the point of closest approach for a dimuon ...
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
static Bool_t TranslateMuonToOrigin(AliAODTrack *muon, Double_t vtx[3])
static Double_t GetPseudoProperDecayTimeXY(Double_t xVtx, Double_t yVtx, Double_t xDimu, Double_t yDimu, Double_t mDimu, Double_t ptDimu)
static Bool_t GetAODMuonOffsetSmeared(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t smearOffsetX, Double_t smearOffsetY, Double_t &offset)
#define TObjArray
static Bool_t GetAODMuonOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset)
Double_t GetBendingCoor() const
return bending coordinate (cm)
void Add(AliMUONVStore &destStore, const AliMUONVStore &srcStore)
Definition: MCHOCCda.cxx:68
TFile f("CalibObjects.root")
static Bool_t GetAODMuonOffsetZ(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset)
static Bool_t ExtrapAODMuonToXY(AliAODTrack *muon, Double_t xy[2], Double_t &zFinal, TLorentzVector &kinem, TMatrixD &cov)
static const TMatrixD ConvertCovMatrixAOD2MUON(AliAODTrack *muon)
static const Int_t fNMaxPlanes
static void GetTrackPrimaryGenerator(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC, TString &nameGen)
Track parameters in ALICE dimuon spectrometer.
AliTPCfastTrack * track
static Bool_t ExtrapToZ(AliMUONTrackParam *trackParam, Double_t zEnd)
static Bool_t IsCorrectMatch(AliAODTrack *muon)
void sum()
void dY()
Definition: CalibAlign.C:547
static Double_t GetPseudoProperDecayTimeZ(Double_t zVtx, Double_t zDimu, Double_t mDimu, Double_t pzDimu)
static Bool_t TranslateMuon(AliAODTrack *muon, Double_t vtxInitial[3], Double_t vtxFinal[3])
static Bool_t ExtrapToZCov(AliMUONTrackParam *trackParam, Double_t zEnd, Bool_t updatePropagator=kFALSE)
Double_t Py() const
static Bool_t IsPDGBeauty(Int_t pdgCode)
static TString GetGenerator(Int_t label, AliAODMCHeader *header)
static void ConvertCovMatrixMUON2AOD(const TMatrixD &covMUON, Double_t covAOD[21])
static const Double_t fZEvalKinem
ClassImp(AliMFTAnalysisTools) Bool_t AliMFTAnalysisTools
static Bool_t IsPDGResonance(Int_t pdgCode)
static Double_t GetDistanceBetweenPoints(TVector3 **points, Int_t nPoints)
static Bool_t IsWrongCluster(AliAODTrack *muon, Int_t iPlane)
static Bool_t IsTrackInjected(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC)
static const Int_t fNMaxMuonsForPCA
Double_t Px() const
AliMUON * muon()
Double_t GetNonBendingCoor() const
return non bending coordinate (cm)
Double_t Pz() const
static Bool_t GetAODMuonWeightedOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset)
static Bool_t ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2])
static Bool_t IsPDGCharm(Int_t pdgCode)
const TMatrixD & GetCovariances() const