]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MFT/AliMFTAnalysisTools.cxx
added methods to define the set of cuts used in the V0 and cascade re-vertexers ...
[u/mrichter/AliRoot.git] / MFT / AliMFTAnalysisTools.cxx
index 7782fde3d3019dd5fe5094ffdec5a097b3b92b6d..da4132c56c74855e126be48174a4d0fd45cc534e 100644 (file)
@@ -31,6 +31,7 @@
 #include "TMath.h"
 #include "AliLog.h"
 #include "TObjArray.h"
+#include "TDecompLU.h"
 
 #include "AliMFTAnalysisTools.h"
 
@@ -40,16 +41,16 @@ ClassImp(AliMFTAnalysisTools)
 
 Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2]) {
 
-  if (!(muon->Pz()!=0)) return kFALSE;
+  if (!(muon->PzAtDCA()!=0)) return kFALSE;
 
   AliMUONTrackParam *param = new AliMUONTrackParam();
 
   param -> SetNonBendingCoor(muon->XAtDCA());
   param -> SetBendingCoor(muon->YAtDCA());
   param -> SetZ(AliMFTConstants::fZEvalKinem);
-  param -> SetNonBendingSlope(muon->Px()/muon->Pz());
-  param -> SetBendingSlope(muon->Py()/muon->Pz());
-  param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) );
+  param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
+  param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
+  param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon->PyAtDCA()/muon->PzAtDCA(),2))) );
 
   AliMUONTrackExtrap::ExtrapToZ(param, z);
   xy[0] = param->GetNonBendingCoor();
@@ -65,16 +66,16 @@ Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Doub
 
 Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2], TLorentzVector &kinem) {
 
-  if (!(muon->Pz()!=0)) return kFALSE;
+  if (!(muon->PzAtDCA()!=0)) return kFALSE;
 
   AliMUONTrackParam *param = new AliMUONTrackParam();
 
   param -> SetNonBendingCoor(muon->XAtDCA());
   param -> SetBendingCoor(muon->YAtDCA());
   param -> SetZ(AliMFTConstants::fZEvalKinem);
-  param -> SetNonBendingSlope(muon->Px()/muon->Pz());
-  param -> SetBendingSlope(muon->Py()/muon->Pz());
-  param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) );
+  param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
+  param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
+  param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon->PyAtDCA()/muon->PzAtDCA(),2))) );
 
   AliMUONTrackExtrap::ExtrapToZ(param, z);
   xy[0] = param->GetNonBendingCoor();
@@ -95,16 +96,18 @@ Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Doub
 
 Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2], TLorentzVector &kinem, TMatrixD &cov) {
 
-  if (!(muon->Pz()!=0)) return kFALSE;
+  // Extrapolate muon to a given z providing the corresponding (x,y) position and updating kinematics and covariance matrix
+
+  if (!(muon->PzAtDCA()!=0)) return kFALSE;
 
   AliMUONTrackParam *param = new AliMUONTrackParam();
 
   param -> SetNonBendingCoor(muon->XAtDCA());
   param -> SetBendingCoor(muon->YAtDCA());
   param -> SetZ(AliMFTConstants::fZEvalKinem);
-  param -> SetNonBendingSlope(muon->Px()/muon->Pz());
-  param -> SetBendingSlope(muon->Py()/muon->Pz());
-  param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) );
+  param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
+  param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
+  param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon->PyAtDCA()/muon->PzAtDCA(),2))) );
 
   param -> SetCovariances(ConvertCovMatrixAOD2MUON(muon));
 
@@ -127,18 +130,142 @@ Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Doub
 
 //====================================================================================================================================================
 
-Double_t AliMFTAnalysisTools::GetAODMuonOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv) {
+Bool_t AliMFTAnalysisTools::ExtrapAODMuonToXY(AliAODTrack *muon, Double_t xy[2], Double_t &zFinal, TLorentzVector &kinem, TMatrixD &cov) {
+
+  // Find the point of closest approach between the muon and the direction parallel to the z-axis defined by the given (x,y)
+  // Provide the z of the above point as weel as the updated kinematics and covariance matrix
+
+  // We look for the above-defined PCA
+
+  AliMUONTrackParam *param = new AliMUONTrackParam();
+  param -> SetNonBendingCoor(muon->XAtDCA());
+  param -> SetBendingCoor(muon->YAtDCA());
+  param -> SetZ(AliMFTConstants::fZEvalKinem);
+  param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
+  param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
+  param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon->PyAtDCA()/muon->PzAtDCA(),2))) );
+  
+  // here we want to understand in which direction we have to search the minimum...
+  
+  Double_t step = 1.;  // initial step, in cm
+  Double_t startPoint = 0.;
+  
+  Double_t r[3]={0}, z[3]={startPoint-step, startPoint, startPoint+step};
+  
+  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)
+  
+  for (Int_t i=0; i<3; i++) {
+    AliMUONTrackExtrap::ExtrapToZ(param, z[i]);
+    points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]);
+    points[1] = new TVector3(xy[0],xy[1],z[i]);
+    r[i] = GetDistanceBetweenPoints(points,2);
+    for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu];
+  }
+  
+  Int_t researchDirection = 0;
+  
+  if      (r[0]>r[1] && r[1]>r[2]) researchDirection = +1;   // towards z positive
+  else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1;   // towards z negative
+  else if (r[0]<r[1] && r[1]>r[2]) {
+    printf("E-AliMFTAnalysisTools::ExtrapAODMuonToXY: Point of closest approach cannot be found (no minima)\n");
+    delete param;
+    delete points;
+    return kFALSE;
+  }
+  
+  while (TMath::Abs(researchDirection)>0.5) {
+      
+    if (researchDirection>0) {
+      z[0] = z[1];
+      z[1] = z[2];
+      z[2] = z[1]+researchDirection*step;
+    }
+    else {
+      z[2] = z[1];
+      z[1] = z[0];
+      z[0] = z[1]+researchDirection*step;
+    }
+    if (TMath::Abs(z[0])>900.) {
+      printf("E-AliMFTAnalysisTools::ExtrapAODMuonToXY: Point of closest approach cannot be found (no minima in the fiducial region)\n");
+      delete param;
+      delete points;
+      return kFALSE;
+    }
+    
+    for (Int_t i=0; i<3; i++) {
+      AliMUONTrackExtrap::ExtrapToZ(param, z[i]);
+      points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]);
+      points[1] = new TVector3(xy[0],xy[1],z[i]);
+      r[i] = GetDistanceBetweenPoints(points,2);
+      for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu];
+    }
+
+    researchDirection=0;
+    if      (r[0]>r[1] && r[1]>r[2]) researchDirection = +1;   // towards z positive
+    else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1;   // towards z negative
+    
+  }
+  
+  // now we now that the minimum is between z[0] and z[2] and we search for it
+  
+  step *= 0.5;
+  while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) {
+    z[0] = z[1]-step;
+    z[2] = z[1]+step;
+    for (Int_t i=0; i<3; i++) {
+      AliMUONTrackExtrap::ExtrapToZ(param, z[i]);
+      points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]);
+      points[1] = new TVector3(xy[0],xy[1],z[i]);
+      r[i] = GetDistanceBetweenPoints(points,2);
+      for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu];
+    }
+    if      (r[0]<r[1]) z[1] = z[0];
+    else if (r[2]<r[1]) z[1] = z[2];
+    else step *= 0.5;
+  }
+  
+  zFinal = z[1];
+
+  Double_t xyMuon[2] = {0};
+  ExtrapAODMuonToZ(muon, zFinal, xyMuon, kinem, cov);
+
+  return kTRUE;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTAnalysisTools::GetAODMuonOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
 
   Double_t xy[2] = {0};
   ExtrapAODMuonToZ(muon, zv, xy);
   
-  return TMath::Sqrt((xv-xy[0])*(xv-xy[0]) + (yv-xy[1])*(yv-xy[1]));
+  offset = TMath::Sqrt((xv-xy[0])*(xv-xy[0]) + (yv-xy[1])*(yv-xy[1]));
+
+  return kTRUE;
 
 }
 
 //====================================================================================================================================================
 
-Double_t AliMFTAnalysisTools::GetAODMuonWeightedOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv) {
+Bool_t AliMFTAnalysisTools::GetAODMuonOffsetZ(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
+
+  Double_t xy[2] = {xv, yv};
+  Double_t zFinal = 0;
+  TLorentzVector kinem(0,0,0,0);
+  TMatrixD cov(5,5);
+
+  ExtrapAODMuonToXY(muon, xy, zFinal, kinem, cov);
+
+  offset = TMath::Abs(zFinal - zv);
+
+  return kTRUE;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTAnalysisTools::GetAODMuonWeightedOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
 
   Double_t xy[2] = {0};
   TLorentzVector kinem(0,0,0,0);
@@ -152,16 +279,23 @@ Double_t AliMFTAnalysisTools::GetAODMuonWeightedOffset(AliAODTrack *muon, Double
   covCoordinates(1,0) = cov(2,0);
   covCoordinates(1,1) = cov(2,2);
   
-  TMatrixD covCoordinatesInverse = covCoordinates.Invert();
+  if (covCoordinates.Determinant() < covCoordinates.GetTol()) return kFALSE;
 
-  Double_t dX = xy[0] - xv;
-  Double_t dY = xy[1] - yv;
-  
-  Double_t weightedOffset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) + 
-                                            dY*dY*covCoordinatesInverse(1,1) + 
-                                            2.*dX*dY*covCoordinatesInverse(0,1)));
+  if (TDecompLU::InvertLU(covCoordinates,covCoordinates.GetTol(),0)) {
 
-  return weightedOffset;
+    TMatrixD covCoordinatesInverse = covCoordinates;
+    Double_t dX = xy[0] - xv;
+    Double_t dY = xy[1] - yv;
+    
+    offset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) +
+                             dY*dY*covCoordinatesInverse(1,1) +
+                             2.*dX*dY*covCoordinatesInverse(0,1)));
+    
+    return kTRUE;
+
+  }
+  
+  return kFALSE;
 
 }
 
@@ -208,8 +342,10 @@ Bool_t AliMFTAnalysisTools::CalculatePCA(AliAODDimuon *dimuon, Double_t *pca, Do
   TObjArray *muons = new TObjArray();
   muons -> Add(dimuon->GetMu(0));
   muons -> Add(dimuon->GetMu(1));
-
-  return CalculatePCA(muons, pca, pcaQuality, kinem);
+  
+  Bool_t result = CalculatePCA(muons, pca, pcaQuality, kinem);
+  delete muons;
+  return result;
 
 }
 
@@ -222,9 +358,9 @@ Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double
     printf("W-AliMFTAnalysisTools::CalculatePCA: number of muons not valid\n");
     return kFALSE;
   }
-  
+
   Double_t fXPointOfClosestApproach=0, fYPointOfClosestApproach=0, fZPointOfClosestApproach=0;
-  
+
   AliAODTrack *muon[AliMFTConstants::fNMaxMuonsForPCA]        = {0};
   AliMUONTrackParam *param[AliMFTConstants::fNMaxMuonsForPCA] = {0};
 
@@ -232,44 +368,54 @@ Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double
   
   for (Int_t iMu=0; iMu<nMuons; iMu++) {
     muon[iMu] = (AliAODTrack*) muons->At(iMu);
-    if (TMath::Abs(muon[iMu]->Pz())<1.e-6) return kFALSE;
+    if (TMath::Abs(muon[iMu]->PzAtDCA())<1.e-6) {
+      for(Int_t i=0;i<iMu;i++) delete param[i];
+      return kFALSE;
+    }
     param[iMu] = new AliMUONTrackParam();
     param[iMu] -> SetNonBendingCoor(muon[iMu]->XAtDCA());
     param[iMu] -> SetBendingCoor(muon[iMu]->YAtDCA());
-    param[iMu] -> SetZ(0.);
-    param[iMu] -> SetNonBendingSlope(muon[iMu]->Px()/muon[iMu]->Pz());
-    param[iMu] -> SetBendingSlope(muon[iMu]->Py()/muon[iMu]->Pz());
-    param[iMu] -> SetInverseBendingMomentum( muon[iMu]->Charge() * (1./muon[iMu]->Pz()) / (TMath::Sqrt(1+TMath::Power(muon[iMu]->Py()/muon[iMu]->Pz(),2))) );
+    param[iMu] -> SetZ(AliMFTConstants::fZEvalKinem);
+    param[iMu] -> SetNonBendingSlope(muon[iMu]->PxAtDCA()/muon[iMu]->PzAtDCA());
+    param[iMu] -> SetBendingSlope(muon[iMu]->PyAtDCA()/muon[iMu]->PzAtDCA());
+    param[iMu] -> SetInverseBendingMomentum( muon[iMu]->Charge() * (1./muon[iMu]->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon[iMu]->PyAtDCA()/muon[iMu]->PzAtDCA(),2))) );
   }
-
+  
   // here we want to understand in which direction we have to search the minimum...
   
   Double_t step = 1.;  // initial step, in cm
   Double_t startPoint = 0.;
-
-  Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step};
+  
+  Double_t r[3]={0}, z[3]={startPoint-step, startPoint, startPoint+step};
   
   TVector3 **points = new TVector3*[AliMFTConstants::fNMaxMuonsForPCA];
   
   for (Int_t i=0; i<3; i++) {
     for (Int_t iMu=0; iMu<nMuons; iMu++) {
+      // if (TMath::Abs(param[iMu]->GetInverseBendingMomentum())<1.) {
+      //       printf("W-AliMFTAnalysisTools::CalculatePCA: Evoiding floating point exception in PCA finding\n");
+      //       return kFALSE;
+      // }
       AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
       points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
     }
     r[i] = GetDistanceBetweenPoints(points,nMuons);
+    for (Int_t iMu=0; iMu<nMuons; iMu++) delete points[iMu]; 
   }
-
+  
   Int_t researchDirection = 0;
   
   if      (r[0]>r[1] && r[1]>r[2]) researchDirection = +1;   // towards z positive
   else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1;   // towards z negative
   else if (r[0]<r[1] && r[1]>r[2]) {
     printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima)\n");
+    for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
+    delete points;
     return kFALSE;
-  }
-
+  }    
+  
   while (TMath::Abs(researchDirection)>0.5) {
-
+      
     if (researchDirection>0) {
       z[0] = z[1];
       z[1] = z[2];
@@ -282,42 +428,54 @@ Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double
     }
     if (TMath::Abs(z[0])>900.) {
       printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima in the fiducial region)\n");
+      for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
+      delete points;
       return kFALSE;
     }
     
     for (Int_t i=0; i<3; i++) {
       for (Int_t iMu=0; iMu<nMuons; iMu++) {
-        AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
-        points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
+       AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
+       points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
       }      
       r[i] = GetDistanceBetweenPoints(points,nMuons);
+      for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu];
     }
     researchDirection=0;
     if      (r[0]>r[1] && r[1]>r[2]) researchDirection = +1;   // towards z positive
     else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1;   // towards z negative
-
+    
   }
-
-  // now we know that the minimum is between z[0] and z[2] and we search for it
   
+  // now we now that the minimum is between z[0] and z[2] and we search for it
+  
+  Int_t nSteps = 0;
+
   step *= 0.5;
   while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) {
     z[0] = z[1]-step;
     z[2] = z[1]+step;
     for (Int_t i=0; i<3; i++) {
       for (Int_t iMu=0; iMu<nMuons; iMu++) {
-        AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
-        points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
+       AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
+       points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
       }      
       r[i] = GetDistanceBetweenPoints(points,nMuons);
+      for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu];
     }
+    //    printf("Step #%d : %f  %f  %f\n",nSteps,r[0],r[1],r[2]);
     if      (r[0]<r[1]) z[1] = z[0];
     else if (r[2]<r[1]) z[1] = z[2];
     else step *= 0.5;
+    nSteps++;
   }
-  
-  //  for (Int_t iMuon=0; iMuon<AliMFTConstants::fNMaxMuonsForPCA; iMuon++) if (points[iMuon]) delete points[iMuon];
-  delete [] points;
+
+  // if (TMath::Abs(z[1]-1.)<0.1) {
+  //   printf("Minimum found in %f in %d steps. Step = %f. p1->X() = %f, p2->X() = %f, p1->Y() = %f, p2->Y() = %f\n",
+  //      z[1], nSteps, step, points[0]->X(), points[1]->X(), points[0]->Y(), points[1]->Y());
+  //   printf("m1->X() = %f, m2->X() = %f, m1->Y() = %f, m2->Y() = %f\n",
+  //      muon[0]->XAtDCA(),muon[1]->XAtDCA(), muon[0]->YAtDCA(),muon[1]->YAtDCA());
+  // }
 
   // Once z of minimum is found, we evaluate the x and y coordinates by averaging over the contributing tracks
   
@@ -337,7 +495,7 @@ Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double
   pca[2] = fZPointOfClosestApproach;
   
   // Evaluating the kinematics of the N-muon
-
+  
   Double_t pTot[3] = {0};
   Double_t ene = 0.;
   Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
@@ -349,22 +507,26 @@ Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double
   }
   
   kinem.SetPxPyPzE(pTot[0], pTot[1], pTot[2], ene);
-
+  
   // Evaluating the PCA quality of the N-muon
-
+  
   Double_t sum=0.,squareSum=0.;
   for (Int_t iMu=0; iMu<nMuons; iMu++) {
-    Double_t wOffset = AliMFTAnalysisTools::GetAODMuonWeightedOffset(muon[iMu],fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach);
+    Double_t wOffset = 0;
+    if (!GetAODMuonWeightedOffset(muon[iMu],fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach, wOffset)) {
+      for(Int_t jMu=0;jMu<nMuons;jMu++) delete param[jMu];
+      delete points;
+      return kFALSE;
+    }
     Double_t f = TMath::Exp(-0.5 * wOffset);
     sum += f;
     squareSum += f*f;
   }
   if (sum > 0.) pcaQuality =  (sum-squareSum/sum) / (nMuons-1);
   else pcaQuality = 0.;
-
-  for (Int_t iMu=0; iMu<AliMFTConstants::fNMaxMuonsForPCA; iMu++) if (param[iMu]) delete param[iMu];
-  //  delete param;
-
+  
+  for(Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
+  delete points;
   return kTRUE;
   
 }
@@ -380,8 +542,8 @@ Double_t AliMFTAnalysisTools::GetDistanceBetweenPoints(TVector3 **points, Int_t
   
   if (nPoints<2) return 0.;
   if (nPoints<3) return TMath::Sqrt( (points[0]->X()-points[1]->X()) * (points[0]->X()-points[1]->X()) +
-                                    (points[0]->Y()-points[1]->Y()) * (points[0]->Y()-points[1]->Y()) +
-                                    (points[0]->Z()-points[1]->Z()) * (points[0]->Z()-points[1]->Z()) );
+                                    (points[0]->Y()-points[1]->Y()) * (points[0]->Y()-points[1]->Y()) );
+  //                                (points[0]->Z()-points[1]->Z()) * (points[0]->Z()-points[1]->Z()) );
 
   const Int_t nEdgesMax = ((AliMFTConstants::fNMaxMuonsForPCA) * (AliMFTConstants::fNMaxMuonsForPCA - 1)) / 2;
   
@@ -528,3 +690,152 @@ const TMatrixD AliMFTAnalysisTools::ConvertCovMatrixAOD2MUON(AliAODTrack *muon)
 
 //====================================================================================================================================================
 
+Bool_t AliMFTAnalysisTools::IsCorrectMatch(AliAODTrack *muon) {
+
+  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) if (IsWrongCluster(muon, iPlane)) return kFALSE;
+  return kTRUE;
+
+}
+
+//====================================================================================================================================================
+
+TString AliMFTAnalysisTools::GetGenerator(Int_t label, AliAODMCHeader* header) {
+
+  // get the name of the generator that produced a given particle
+
+  Int_t partCounter = 0;
+  TList *genHeaders = header->GetCocktailHeaders();
+  Int_t nGenHeaders = genHeaders->GetEntries();
+
+  for (Int_t i=0; i<nGenHeaders; i++){
+    AliGenEventHeader *gh = (AliGenEventHeader*) genHeaders->At(i);
+    TString genName = gh->GetName();
+    Int_t nPart = gh->NProduced();
+    if (label>=partCounter && label<(partCounter+nPart)) return genName;
+    partCounter += nPart;
+  }
+
+  TString empty="";
+  return empty;
+
+}
+
+//====================================================================================================================================================
+
+void AliMFTAnalysisTools::GetTrackPrimaryGenerator(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC, TString &nameGen) {
+
+  // method to check if a track comes from a given generator
+
+  Int_t label = TMath::Abs(track->GetLabel());
+  nameGen = GetGenerator(label,header);
+  
+  // 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
+  
+  while (nameGen.IsWhitespace()) {
+    AliAODMCParticle *mcPart = (AliAODMCParticle*) arrayMC->At(label);
+    if (!mcPart) {
+      printf("AliMFTAnalysisTools::GetTrackPrimaryGenerator - BREAK: No valid AliAODMCParticle at label %i\n",label);
+      break;
+    }
+    Int_t motherLabel = mcPart->GetMother();
+    if (motherLabel < 0) {
+      printf("AliMFTAnalysisTools::GetTrackPrimaryGenerator - BREAK: Reached primary particle without valid mother\n");
+      break;
+    }
+    label = motherLabel;
+    nameGen = GetGenerator(label,header);
+  }
+  
+  return;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTAnalysisTools::IsTrackInjected(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC) {
+
+  // method to check if a track comes from the signal event or from the underlying Hijing event
+
+  TString nameGen;
+
+  GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
+  
+  if (nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
+  
+  return kTRUE;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTAnalysisTools::TranslateMuon(AliAODTrack *muon, Double_t vtxInitial[3], Double_t vtxFinal[3]) {
+
+  if (!(muon->PzAtDCA()!=0)) return kFALSE;
+
+  AliMUONTrackParam *param = new AliMUONTrackParam();
+
+  Double_t deltaVtx[3] = {0};
+  for (Int_t i=0; i<3; i++) deltaVtx[i] = vtxInitial[i] - vtxFinal[i];
+
+  param -> SetNonBendingCoor(muon->XAtDCA());
+  param -> SetBendingCoor(muon->YAtDCA());
+  param -> SetZ(AliMFTConstants::fZEvalKinem);
+  param -> SetNonBendingSlope(muon->PxAtDCA()/muon->PzAtDCA());
+  param -> SetBendingSlope(muon->PyAtDCA()/muon->PzAtDCA());
+  param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->PzAtDCA()) / (TMath::Sqrt(1+TMath::Power(muon->PyAtDCA()/muon->PzAtDCA(),2))) );
+
+  // This will be interpreted as if the track (produced in an event having the primary vertex in vtxInitial) 
+  // were produced in an event having the primary vertex in vtxFinal
+
+  AliMUONTrackExtrap::ExtrapToZ(param, deltaVtx[2]);
+  muon->SetXYAtDCA(param->GetNonBendingCoor() - deltaVtx[0], param->GetBendingCoor() - deltaVtx[1]);
+  muon->SetPxPyPzAtDCA(param->Px(), param->Py(), param->Pz());
+
+  delete param;
+
+  return kTRUE;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTAnalysisTools::TranslateMuonToOrigin(AliAODTrack *muon, Double_t vtx[3]) {
+
+  Double_t origin[3] = {0,0,0};
+
+  return TranslateMuon(muon, vtx, origin);
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTAnalysisTools::IsPDGCharm(Int_t pdgCode) {
+
+  pdgCode = TMath::Abs(pdgCode/100);
+  if (pdgCode>9) pdgCode /= 10;
+  if (pdgCode == 4 ) return kTRUE;
+  else return kFALSE;
+  
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTAnalysisTools::IsPDGBeauty(Int_t pdgCode) {
+
+  pdgCode = TMath::Abs(pdgCode/100);
+  if (pdgCode>9) pdgCode /= 10;
+  if (pdgCode == 5) return kTRUE;
+  else return kFALSE;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTAnalysisTools::IsPDGResonance(Int_t pdgCode) {
+
+  Int_t id = pdgCode%100000;
+  return (!((id-id%10)%110));
+
+} 
+
+//====================================================================================================================================================