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