MFT/MFTbase/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 #include "TRandom.h"
  36
  37 #include "AliMFTAnalysisTools.h"
  38
  39 ClassImp(AliMFTAnalysisTools)
  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
  56   AliMUONTrackExtrap::ExtrapToZ(param, z);
  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
  81   AliMUONTrackExtrap::ExtrapToZ(param, z);
  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
 115   AliMUONTrackExtrap::ExtrapToZCov(param, z);
 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;
 213   while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) {
 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
 342 Double_t AliMFTAnalysisTools::GetPseudoProperDecayTimeZ(Double_t zVtx,
 343                                                         Double_t zDimu,
 344                                                         Double_t mDimu, Double_t pzDimu) {
 345
 346   // pseudo-proper decay time of a particle produced in the primary vertex and decaying into a dimuon (+ X)
 347   // evaluated using the longitudinal degree of freedom of the decay topology
 348
 349   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};
 385   AliMUONTrackParam *param[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;
 475   while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) {
 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
 680 const TMatrixD AliMFTAnalysisTools::ConvertCovMatrixAOD2MUON(AliAODTrack *muon) {
 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
 713 Bool_t AliMFTAnalysisTools::IsCorrectMatch(AliAODTrack *muon) {
 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 //====================================================================================================================================================