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