[u/mrichter/AliRoot.git] / MFT / AliMFTAnalysisTools.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//====================================================================================================================================================
//
//      Support class for various common operation on MFT objects
//
//      Contact author: antonio.uras@cern.ch
//
//====================================================================================================================================================

#include "AliMUONTrackParam.h"
#include "AliMUONTrackExtrap.h"
#include "AliAODTrack.h"
#include "AliAODDimuon.h"
#include "TLorentzVector.h"
#include "AliMFTConstants.h"
#include "TDatabasePDG.h"
#include "TMath.h"
#include "AliLog.h"
#include "TObjArray.h"
#include "TDecompLU.h"

#include "AliMFTAnalysisTools.h"

ClassImp(AliMFTAnalysisTools)

//====================================================================================================================================================

Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2]) {

  if (!(muon->Pz()!=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))) );

  AliMUONTrackExtrap::ExtrapToZ(param, z);
  xy[0] = param->GetNonBendingCoor();
  xy[1] = param->GetBendingCoor();

  delete param;

  return kTRUE;

}

//====================================================================================================================================================

Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2], TLorentzVector &kinem) {

  if (!(muon->Pz()!=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))) );

  AliMUONTrackExtrap::ExtrapToZ(param, z);
  xy[0] = param->GetNonBendingCoor();
  xy[1] = param->GetBendingCoor();

  Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
  Double_t energy = TMath::Sqrt(massMu*massMu + param->Px()*param->Px() + param->Py()*param->Py() + param->Pz()*param->Pz());
  
  kinem.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);

  delete param;

  return kTRUE;

}

//====================================================================================================================================================

Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2], TLorentzVector &kinem, TMatrixD &cov) {

  if (!(muon->Pz()!=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 -> SetCovariances(ConvertCovMatrixAOD2MUON(muon));

  AliMUONTrackExtrap::ExtrapToZCov(param, z);
  xy[0] = param->GetNonBendingCoor();
  xy[1] = param->GetBendingCoor();

  Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
  Double_t energy = TMath::Sqrt(massMu*massMu + param->Px()*param->Px() + param->Py()*param->Py() + param->Pz()*param->Pz());
  
  kinem.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);

  cov = param->GetCovariances();

  delete param;

  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);
  
  offset = TMath::Sqrt((xv-xy[0])*(xv-xy[0]) + (yv-xy[1])*(yv-xy[1]));

  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);
  TMatrixD cov(5,5);

  ExtrapAODMuonToZ(muon, zv, xy, kinem, cov);

  TMatrixD covCoordinates(2,2);
  covCoordinates(0,0) = cov(0,0);
  covCoordinates(0,1) = cov(0,2);
  covCoordinates(1,0) = cov(2,0);
  covCoordinates(1,1) = cov(2,2);
  
  if (TDecompLU::InvertLU(covCoordinates,covCoordinates.GetTol(),0)) {

    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;

}

//====================================================================================================================================================

Double_t AliMFTAnalysisTools::GetPseudoProperDecayTimeXY(Double_t xVtx,  Double_t yVtx, 
							 Double_t xDimu, Double_t yDimu, 
							 Double_t mDimu, Double_t ptDimu) {
  
  // pseudo-proper decay time of a particle produced in the primary vertex and decaying into a dimuon (+ X)
  // evaluated using the transverse degree of freedom of the decay topology 

  if (ptDimu != 0) {
    Double_t decayLengthXY = TMath::Sqrt((xVtx-xDimu)*(xVtx-xDimu)+(yVtx-yDimu)*(yVtx-yDimu));
    return (decayLengthXY * mDimu/ptDimu)/TMath::Ccgs()*1E12;   // in ps
  }
  
  return -99999999;
  
}

//====================================================================================================================================================

Double_t AliMFTAnalysisTools::GetPseudoProperDecayTimeZ(Double_t zVtx, 
							Double_t zDimu, 
							Double_t mDimu, Double_t pzDimu) {

  // pseudo-proper decay time of a particle produced in the primary vertex and decaying into a dimuon (+ X)
  // evaluated using the longitudinal degree of freedom of the decay topology 

  if (pzDimu != 0) {
    Double_t decayLengthZ = zDimu - zVtx;
    return (decayLengthZ * mDimu/pzDimu)/TMath::Ccgs()*1E12;  // in ps
  }

  return -99999999;

}

//====================================================================================================================================================

Bool_t AliMFTAnalysisTools::CalculatePCA(AliAODDimuon *dimuon, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {

  TObjArray *muons = new TObjArray();
  muons -> Add(dimuon->GetMu(0));
  muons -> Add(dimuon->GetMu(1));
  
  Bool_t result = CalculatePCA(muons, pca, pcaQuality, kinem);
  delete muons;
  return result;

}

//====================================================================================================================================================

Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
  
  const Int_t nMuons = muons->GetEntriesFast();
  if (nMuons<2 || nMuons>AliMFTConstants::fNMaxMuonsForPCA) {
    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};

  // Finding AliMUONTrackParam objects for each muon
  
  for (Int_t iMu=0; iMu<nMuons; iMu++) {
    muon[iMu] = (AliAODTrack*) muons->At(iMu);
    if (TMath::Abs(muon[iMu]->Pz())<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))) );
  }
  
  // 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};
  
  TVector3 **points = new TVector3*[AliMFTConstants::fNMaxMuonsForPCA];
  
  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]);
  	  	}
    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];
      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::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]);
      }      
      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 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++) {
      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]);
      }      
      r[i] = GetDistanceBetweenPoints(points,nMuons);
      for (Int_t iMu=0;iMu<nMuons;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;
  }
  
  // Once z of minimum is found, we evaluate the x and y coordinates by averaging over the contributing tracks
  
  fZPointOfClosestApproach = z[1];
  fXPointOfClosestApproach = 0.;
  fYPointOfClosestApproach = 0.;
  for (Int_t iMu=0; iMu<nMuons; iMu++) {
    AliMUONTrackExtrap::ExtrapToZ(param[iMu], fZPointOfClosestApproach);
    fXPointOfClosestApproach += param[iMu]->GetNonBendingCoor();
    fYPointOfClosestApproach += param[iMu]->GetBendingCoor();
  }
  fXPointOfClosestApproach /= Double_t(nMuons);
  fYPointOfClosestApproach /= Double_t(nMuons);
  
  pca[0] = fXPointOfClosestApproach;
  pca[1] = fYPointOfClosestApproach;
  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();
  for (Int_t iMu=0; iMu<nMuons; iMu++) {
    pTot[0] += param[iMu]->Px();
    pTot[1] += param[iMu]->Py();
    pTot[2] += param[iMu]->Pz();
    ene += TMath::Sqrt(massMu*massMu + param[iMu]->Px()*param[iMu]->Px() + param[iMu]->Py()*param[iMu]->Py() + param[iMu]->Pz()*param[iMu]->Pz());
  }
  
  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 = 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<nMuons;iMu++) delete param[iMu];
  delete points;
  return kTRUE;
  
}

//=========================================================================================================================

Double_t AliMFTAnalysisTools::GetDistanceBetweenPoints(TVector3 **points, Int_t nPoints) {
  
  if (nPoints>AliMFTConstants::fNMaxMuonsForPCA) {
    printf("W-AliMFTAnalysisTools::GetDistanceBetweenPoints: number of points not valid\n");
    return 1.e9;
  }
  
  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()) );

  const Int_t nEdgesMax = ((AliMFTConstants::fNMaxMuonsForPCA) * (AliMFTConstants::fNMaxMuonsForPCA - 1)) / 2;
  
  Int_t startID[nEdgesMax]       = {0};
  Int_t stopID[nEdgesMax]        = {0};
  Double_t edgeLength[nEdgesMax] = {0};

  Bool_t pointStatus[AliMFTConstants::fNMaxMuonsForPCA] = {0};
  
  Int_t nEdges=0;
  for (Int_t i=0; i<nPoints-1; i++) {
    for (Int_t j=i+1; j<nPoints; j++) {
      edgeLength[nEdges] = TMath::Sqrt( (points[i]->X()-points[j]->X()) * (points[i]->X()-points[j]->X()) +
                                        (points[i]->Y()-points[j]->Y()) * (points[i]->Y()-points[j]->Y()) +
                                        (points[i]->Z()-points[j]->Z()) * (points[i]->Z()-points[j]->Z()) );
      stopID[nEdges]  = i;
      startID[nEdges] = j;
      nEdges++;
    }
  }
 
  // Order Edges

  Double_t min = 0;
  Int_t   iMin = 0;

  for (Int_t iEdge=0; iEdge<nEdges-1; iEdge++) {
    min  = edgeLength[iEdge];
    iMin = iEdge;
    for (Int_t j=iEdge+1; j<nEdges; j++) {
      if (edgeLength[j]<min) {
        min  = edgeLength[j];
        iMin = j;
      }
    }
    
    if (iMin != iEdge) {

      Double_t edgeLengthMin = edgeLength[iMin];
      Int_t startIDmin = startID[iMin];
      Int_t stopIDmin  = stopID[iMin];
      
      edgeLength[iMin] = edgeLength[iEdge];
      startID[iMin]    = startID[iEdge];
      stopID[iMin]     = stopID[iEdge];

      edgeLength[iEdge] = edgeLengthMin;
      startID[iEdge]    = startIDmin;
      stopID[iEdge]     = stopIDmin;

    }
    
  }
  
  // Connect

  Double_t length = 0.;
  for (Int_t i=0; i<nEdges; i++) {
    if (!(pointStatus[startID[i]] && pointStatus[stopID[i]])) {
      pointStatus[startID[i]] = kTRUE;
      pointStatus[stopID[i]]  = kTRUE;
      length += edgeLength[i];
    }
  }
  
  return length;
  
}

//====================================================================================================================================================

void AliMFTAnalysisTools::ConvertCovMatrixMUON2AOD(const TMatrixD& covMUON, Double_t covAOD[21]) {

  // Converts the cov matrix from the MUON format (TMatrixD) to the AOD one (Double_t[21])
  // 
  // Cov(x,x)       ... :   cv[0]
  // Cov(x,slopeX)  ... :   cv[1]  cv[2]
  // Cov(x,y)       ... :   cv[3]  cv[4]  cv[5]
  // Cov(x,slopeY)  ... :   cv[6]  cv[7]  cv[8]  cv[9]
  // Cov(x,invP_yz) ... :   cv[10] cv[11] cv[12] cv[13] cv[14]
  // not-used       ... :   cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]

  covAOD[0]  = covMUON(0,0);

  covAOD[1]  = covMUON(1,0);
  covAOD[2]  = covMUON(1,1);

  covAOD[3]  = covMUON(2,0);  
  covAOD[4]  = covMUON(2,1);  
  covAOD[5]  = covMUON(2,2);  

  covAOD[6]  = covMUON(3,0);  
  covAOD[7]  = covMUON(3,1);  
  covAOD[8]  = covMUON(3,2);  
  covAOD[9]  = covMUON(3,3);  

  covAOD[10] = covMUON(4,0);  
  covAOD[11] = covMUON(4,1);  
  covAOD[12] = covMUON(4,2);  
  covAOD[13] = covMUON(4,3);  
  covAOD[14] = covMUON(4,4);  

  covAOD[15] = 0;  
  covAOD[16] = 0;  
  covAOD[17] = 0;  
  covAOD[18] = 0;  
  covAOD[19] = 0;  
  covAOD[20] = 0;  

}

//====================================================================================================================================================

const TMatrixD AliMFTAnalysisTools::ConvertCovMatrixAOD2MUON(AliAODTrack *muon) {

  Double_t covAOD[21] = {0};
  muon -> GetCovarianceXYZPxPyPz(covAOD);

  TMatrixD covMUON(5,5);

  covMUON(0,0) = covAOD[0];
  	         		
  covMUON(1,0) = covAOD[1];
  covMUON(1,1) = covAOD[2];
  	         		
  covMUON(2,0) = covAOD[3];  
  covMUON(2,1) = covAOD[4];  
  covMUON(2,2) = covAOD[5];  
  	         		
  covMUON(3,0) = covAOD[6];  
  covMUON(3,1) = covAOD[7];  
  covMUON(3,2) = covAOD[8];  
  covMUON(3,3) = covAOD[9];  
  	         		
  covMUON(4,0) = covAOD[10];  
  covMUON(4,1) = covAOD[11];  
  covMUON(4,2) = covAOD[12];  
  covMUON(4,3) = covAOD[13];  
  covMUON(4,4) = covAOD[14]; 

  return covMUON;

}

//====================================================================================================================================================
Commit	Line	Data
	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(massMumassMu + 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(massMumassMu + 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(dXdX*covCoordinatesInverse(0,0) +
	165	dYdYcovCoordinatesInverse(1,1) +
	166	2.dXdY*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(massMumassMu + 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