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 //====================================================================================================================================================