[u/mrichter/AliRoot.git] / STEER / AliAODTrack.cxx

/**************************************************************************
 * Copyright(c) 1998-2007, 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.                  *
 **************************************************************************/

/* $Id$ */

//-------------------------------------------------------------------------
//     AOD track implementation of AliVParticle
//     Author: Markus Oldenburg, CERN
//     Markus.Oldenburg@cern.ch
//-------------------------------------------------------------------------

#include "AliLog.h"
#include "AliAODTrack.h"

ClassImp(AliAODTrack)

//______________________________________________________________________________
AliAODTrack::AliAODTrack() : 
  AliVParticle(),
  fChi2perNDF(-999.),
  fChi2MatchTrigger(0.),
  fFlags(0),
  fLabel(-999),
  fITSMuonClusterMap(0),
  fFilterMap(0),
  fID(-999),
  fCharge(-99),
  fType(kUndef),
  fCovMatrix(NULL),
  fDetPid(NULL),
  fProdVertex(NULL)
{
  // default constructor

  SetP();
  SetPosition((Float_t*)NULL);
  SetPID((Float_t*)NULL);
}

//______________________________________________________________________________
AliAODTrack::AliAODTrack(Short_t id,
			 Int_t label, 
			 Double_t p[3],
			 Bool_t cartesian,
			 Double_t x[3],
			 Bool_t isDCA,
			 Double_t covMatrix[21],
			 Short_t charge,
			 UChar_t itsClusMap,
			 Double_t pid[10],
			 AliAODVertex *prodVertex,
			 Bool_t usedForVtxFit,
			 Bool_t usedForPrimVtxFit,
			 AODTrk_t ttype,
			 UInt_t selectInfo) :
  AliVParticle(),
  fChi2perNDF(-999.),
  fChi2MatchTrigger(0.),
  fFlags(0),
  fLabel(label),
  fITSMuonClusterMap(itsClusMap),
  fFilterMap(selectInfo),
  fID(id),
  fCharge(charge),
  fType(ttype),
  fCovMatrix(NULL),
  fDetPid(NULL),
  fProdVertex(prodVertex)
{
  // constructor
 
  SetP(p, cartesian);
  SetPosition(x, isDCA);
  SetUsedForVtxFit(usedForVtxFit);
  SetUsedForPrimVtxFit(usedForPrimVtxFit);
  if(covMatrix) SetCovMatrix(covMatrix);
  SetPID(pid);

}

//______________________________________________________________________________
AliAODTrack::AliAODTrack(Short_t id,
			 Int_t label, 
			 Float_t p[3],
			 Bool_t cartesian,
			 Float_t x[3],
			 Bool_t isDCA,
			 Float_t covMatrix[21],
			 Short_t charge,
			 UChar_t itsClusMap,
			 Float_t pid[10],
			 AliAODVertex *prodVertex,
			 Bool_t usedForVtxFit,
			 Bool_t usedForPrimVtxFit,
			 AODTrk_t ttype,
			 UInt_t selectInfo) :
  AliVParticle(),
  fChi2perNDF(-999.),
  fChi2MatchTrigger(0.),
  fFlags(0),
  fLabel(label),
  fITSMuonClusterMap(itsClusMap),
  fFilterMap(selectInfo),
  fID(id),
  fCharge(charge),
  fType(ttype),
  fCovMatrix(NULL),
  fDetPid(NULL),
  fProdVertex(prodVertex)
{
  // constructor
 
  SetP(p, cartesian);
  SetPosition(x, isDCA);
  SetUsedForVtxFit(usedForVtxFit);
  SetUsedForPrimVtxFit(usedForPrimVtxFit);
  if(covMatrix) SetCovMatrix(covMatrix);
  SetPID(pid);
}

//______________________________________________________________________________
AliAODTrack::~AliAODTrack() 
{
  // destructor
  delete fCovMatrix;
}


//______________________________________________________________________________
AliAODTrack::AliAODTrack(const AliAODTrack& trk) :
  AliVParticle(trk),
  fChi2perNDF(trk.fChi2perNDF),
  fChi2MatchTrigger(trk.fChi2MatchTrigger),
  fFlags(trk.fFlags),
  fLabel(trk.fLabel),
  fITSMuonClusterMap(trk.fITSMuonClusterMap),
  fFilterMap(trk.fFilterMap),
  fID(trk.fID),
  fCharge(trk.fCharge),
  fType(trk.fType),
  fCovMatrix(NULL),
  fDetPid(NULL),
  fProdVertex(trk.fProdVertex)
{
  // Copy constructor

  trk.GetP(fMomentum);
  trk.GetPosition(fPosition);
  SetUsedForVtxFit(trk.GetUsedForVtxFit());
  SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
  if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
  if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
  SetPID(trk.fPID);
}

//______________________________________________________________________________
AliAODTrack& AliAODTrack::operator=(const AliAODTrack& trk)
{
  // Assignment operator
  if(this!=&trk) {

    AliVParticle::operator=(trk);

    trk.GetP(fMomentum);
    trk.GetPosition(fPosition);
    trk.GetPID(fPID);

    fChi2perNDF = trk.fChi2perNDF;
    fChi2MatchTrigger = trk.fChi2MatchTrigger;

    fFlags = trk.fFlags;
    fLabel = trk.fLabel;    
    
    fITSMuonClusterMap = trk.fITSMuonClusterMap;
    fFilterMap = trk.fFilterMap;

    fID = trk.fID;

    fCharge = trk.fCharge;
    fType = trk.fType;

    delete fCovMatrix;
    if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
    else fCovMatrix=NULL;
    fProdVertex = trk.fProdVertex;

    SetUsedForVtxFit(trk.GetUsedForVtxFit());
    SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());

    delete fDetPid;
    if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
    else fDetPid=NULL;
  }

  return *this;
}

//______________________________________________________________________________
Double_t AliAODTrack::M(AODTrkPID_t pid) const
{
  // Returns the mass.
  // Masses for nuclei don't exist in the PDG tables, therefore they were put by hand.

  switch (pid) {

  case kElectron :
    return 0.000510999; //TDatabasePDG::Instance()->GetParticle(11/*::kElectron*/)->Mass();
    break;

  case kMuon :
    return 0.1056584; //TDatabasePDG::Instance()->GetParticle(13/*::kMuonMinus*/)->Mass();
    break;

  case kPion :
    return 0.13957; //TDatabasePDG::Instance()->GetParticle(211/*::kPiPlus*/)->Mass();
    break;

  case kKaon :
    return 0.4937; //TDatabasePDG::Instance()->GetParticle(321/*::kKPlus*/)->Mass();
    break;

  case kProton :
    return 0.9382720; //TDatabasePDG::Instance()->GetParticle(2212/*::kProton*/)->Mass();
    break;

  case kDeuteron :
    return 1.8756; //TDatabasePDG::Instance()->GetParticle(1000010020)->Mass();
    break;

  case kTriton :
    return 2.8089; //TDatabasePDG::Instance()->GetParticle(1000010030)->Mass();
    break;

  case kHelium3 :
    return 2.8084; //TDatabasePDG::Instance()->GetParticle(1000020030)->Mass();
    break;

  case kAlpha :
    return 3.7274; //TDatabasePDG::Instance()->GetParticle(1000020040)->Mass();
    break;

  case kUnknown :
    return -999.;
    break;

  default :
    return -999.;
  }
}

//______________________________________________________________________________
Double_t AliAODTrack::E(AODTrkPID_t pid) const
{
  // Returns the energy of the particle of a given pid.
  
  if (pid != kUnknown) { // particle was identified
    Double_t m = M(pid);
    return TMath::Sqrt(P()*P() + m*m);
  } else { // pid unknown
    return -999.;
  }
}

//______________________________________________________________________________
Double_t AliAODTrack::Y(AODTrkPID_t pid) const
{
  // Returns the rapidity of a particle of a given pid.
  
  if (pid != kUnknown) { // particle was identified
    Double_t e = E(pid);
    Double_t pz = Pz();
    if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
      return 0.5*TMath::Log((e+pz)/(e-pz));
    } else { // energy not known or equal to pz
      return -999.;
    }
  } else { // pid unknown
    return -999.;
  }
}

//______________________________________________________________________________
Double_t AliAODTrack::Y(Double_t m) const
{
  // Returns the rapidity of a particle of a given mass.
  
  if (m >= 0.) { // mass makes sense
    Double_t e = E(m);
    Double_t pz = Pz();
    if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
      return 0.5*TMath::Log((e+pz)/(e-pz));
    } else { // energy not known or equal to pz
      return -999.;
    }
  } else { // pid unknown
    return -999.;
  }
}

//______________________________________________________________________________
AliAODTrack::AODTrkPID_t AliAODTrack::GetMostProbablePID() const 
{
  // Returns the most probable PID array element.
  
  Int_t nPID = 10;
  if (fPID) {
    AODTrkPID_t loc = kUnknown;
    Double_t max = 0.;
    Bool_t allTheSame = kTRUE;
    
    for (Int_t iPID = 0; iPID < nPID; iPID++) {
      if (fPID[iPID] >= max) {
	if (fPID[iPID] > max) {
	  allTheSame = kFALSE;
	  max = fPID[iPID];
	  loc = (AODTrkPID_t)iPID;
	} else {
	  allTheSame = kTRUE;
	}
      }
    }
    
    return allTheSame ? kUnknown : loc;
  } else {
    return kUnknown;
  }
}

//______________________________________________________________________________
void AliAODTrack::ConvertAliPIDtoAODPID()
{
  // Converts AliPID array.
  // The numbering scheme is the same for electrons, muons, pions, kaons, and protons.
  // Everything else has to be set to zero.

  fPID[kDeuteron] = 0.;
  fPID[kTriton]   = 0.;
  fPID[kHelium3]  = 0.;
  fPID[kAlpha]    = 0.;
  fPID[kUnknown]  = 0.;
  
  return;
}


//______________________________________________________________________________
template <class T> void AliAODTrack::SetP(const T *p, const Bool_t cartesian) 
{
  // Set the momentum

  if (p) {
    if (cartesian) {
      Double_t pt2 = p[0]*p[0] + p[1]*p[1];
      Double_t pp  = TMath::Sqrt(pt2 + p[2]*p[2]);
      
      fMomentum[0] = TMath::Sqrt(pt2); // pt
      fMomentum[1] = (pt2 != 0.) ? TMath::Pi()+TMath::ATan2(-p[1], -p[0]) : -999; // phi
      fMomentum[2] = (pp != 0.) ? TMath::ACos(p[2] / pp) : -999.; // theta
    } else {
      fMomentum[0] = p[0];  // pt
      fMomentum[1] = p[1];  // phi
      fMomentum[2] = p[2];  // theta
    }
  } else {
    fMomentum[0] = -999.;
    fMomentum[1] = -999.;
    fMomentum[2] = -999.;
  }
}

//______________________________________________________________________________
template <class T> void AliAODTrack::SetPosition(const T *x, const Bool_t dca) 
{
  // set the position

  if (x) {
    if (!dca) {
      ResetBit(kIsDCA);

      fPosition[0] = x[0];
      fPosition[1] = x[1];
      fPosition[2] = x[2];
    } else {
      SetBit(kIsDCA);
      // don't know any better yet
      fPosition[0] = -999.;
      fPosition[1] = -999.;
      fPosition[2] = -999.;
    }
  } else {
    ResetBit(kIsDCA);

    fPosition[0] = -999.;
    fPosition[1] = -999.;
    fPosition[2] = -999.;
  }
}

//______________________________________________________________________________
void AliAODTrack::SetDCA(Double_t d, Double_t z) 
{
  // set the dca
  fPosition[0] = d;
  fPosition[1] = z;
  fPosition[2] = 0.;
  SetBit(kIsDCA);
}

//______________________________________________________________________________
void AliAODTrack::Print(Option_t* /* option */) const
{
  // prints information about AliAODTrack

  printf("Object name: %s   Track type: %s\n", GetName(), GetTitle()); 
  printf("        px = %f\n", Px());
  printf("        py = %f\n", Py());
  printf("        pz = %f\n", Pz());
  printf("        pt = %f\n", Pt());
  printf("      1/pt = %f\n", OneOverPt());
  printf("     theta = %f\n", Theta());
  printf("       phi = %f\n", Phi());
  printf("  chi2/NDF = %f\n", Chi2perNDF());
  printf("    charge = %d\n", Charge());
}

void AliAODTrack::SetMatchTrigger(Int_t matchTrig){
//
// Set the MUON trigger information
  switch(matchTrig){
    case 0: // 0 track does not match trigger
      fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
      break;
    case 1: // 1 track match but does not pass pt cut
      fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x40000000;
      break;
    case 2: // 2 track match Low pt cut
      fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x80000000;
      break;
    case 3: // 3 track match High pt cut
      fITSMuonClusterMap=fITSMuonClusterMap|0xc0000000;
      break;
    default:
      fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
      AliWarning(Form("unknown case for matchTrig: %d\n",matchTrig));
  }
}

void AliAODTrack::SetHitsPatternInTrigCh(UShort_t hitsPatternInTrigCh){
//
// Set the MUON hit pattern (1 bit per chamber) 
  fITSMuonClusterMap=(fITSMuonClusterMap&0xffff00ff)|(hitsPatternInTrigCh<<8);
}

Int_t AliAODTrack::HitsMT(Int_t istation, Int_t iplane, Char_t *cathode){
//
// Retrieve hit information for MUON identified by  (station, plane, cathode)
  if(cathode){
    if(cathode[0]=='x'||cathode[0]=='X'){
      if(istation==1){
        if(iplane==1)
	  return (fITSMuonClusterMap&0x8000)?1:0;
	else if(iplane==2)
	  return (fITSMuonClusterMap&0x4000)?1:0;
	else
	  return 0;
      }else if(istation==2){
        if(iplane==1)
	  return (fITSMuonClusterMap&0x2000)?1:0;
	else if(iplane==2)
	  return (fITSMuonClusterMap&0x1000)?1:0;
	else
	  return 0;
      }else{
        return 0;
      }
    }else if(cathode[0]=='y'||cathode[0]=='Y'){
      if(istation==1){
        if(iplane==1)
	  return (fITSMuonClusterMap&0x0800)?1:0;
	else if(iplane==2)
	  return (fITSMuonClusterMap&0x0400)?1:0;
	else
	  return 0;
      }else if(istation==2){
        if(iplane==1)
	  return (fITSMuonClusterMap&0x0200)?1:0;
	else if(iplane==2)
	  return (fITSMuonClusterMap&0x0100)?1:0;
	else
	  return 0;
      }else{
        return 0;
      }
    }else{
      return 0;
    }
  }else{
    if(istation==1){
      if(iplane==1)
	return (HitsMT(1,1,"X")||HitsMT(1,1,"Y"))?1:0;
      else if(iplane==2)
	return (HitsMT(1,2,"X")||HitsMT(1,2,"Y"))?1:0;
      else
	return 0;
    }else if(istation==2){
      if(iplane==1)
	return (HitsMT(2,1,"X")||HitsMT(2,1,"Y"))?1:0;
      else if(iplane==2)
	return (HitsMT(2,2,"X")||HitsMT(2,2,"Y"))?1:0;
      else
	return 0;
    }else{
      return 0;
    }
  }
}

Int_t AliAODTrack::HitsMuonChamber(Int_t MuonChamber){
// Retrieve hit information for MUON Chamber
  switch(MuonChamber){
    case 11:
      return HitsMT(1,1);
    case 12:
      return HitsMT(1,2);
    case 13:
      return HitsMT(2,1);
    case 14:
      return HitsMT(2,2);
    default:
      printf("Unknown MUON chamber: %d\n",MuonChamber);
      return 0;
  }
}


Bool_t AliAODTrack::PropagateTo(Double_t xk, Double_t b) {
  //----------------------------------------------------------------
  // Propagate this track to the plane X=xk (cm) in the field "b" (kG)
  // This is in local coordinates!!!
  //----------------------------------------------------------------

  Double_t alpha = 0.;
  Double_t localP[3] = {Px(), Py(), Pz()}; // set global (sic!) p
  Global2LocalMomentum(localP, Charge(), alpha); // convert global to local momentum

  AliAODVertex *origin = (AliAODVertex*)fProdVertex.GetObject();
  Double_t localX[3] = {origin->GetX(), origin->GetY(), origin->GetZ()}; // set global (sic!) location of first track point
  Global2LocalPosition(localX, alpha); // convert global to local position

  Double_t &fX = localX[0];

  Double_t dx=xk-fX;
  if (TMath::Abs(dx)<=kAlmost0)  return kTRUE;

  Double_t crv=localP[0]*b*kB2C;
  if (TMath::Abs(b) < kAlmost0Field) crv=0.;

  Double_t f1=localP[1], f2=f1 + crv*dx;
  if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
  if (TMath::Abs(f2) >= kAlmost1) return kFALSE;

  Double_t &fP0=localX[1], &fP1=localX[2], &fP2=localP[0], &fP3=localP[1], &fP4=localP[2];
  /* covariance matrix to be fixed! 
  Double_t 
  &fC00=fC[0],
  &fC10=fC[1],   &fC11=fC[2],  
  &fC20=fC[3],   &fC21=fC[4],   &fC22=fC[5],
  &fC30=fC[6],   &fC31=fC[7],   &fC32=fC[8],   &fC33=fC[9],  
  &fC40=fC[10],  &fC41=fC[11],  &fC42=fC[12],  &fC43=fC[13], &fC44=fC[14];
  */
  Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);

  fX=xk;
  fP0 += dx*(f1+f2)/(r1+r2);
  fP1 += dx*(r2 + f2*(f1+f2)/(r1+r2))*fP3;
  fP2 += dx*crv;

  //f = F - 1
   
  //Double_t f02=    dx/(r1*r1*r1);            
  Double_t cc=crv/fP4;
  Double_t f04=0.5*dx*dx/(r1*r1*r1);         f04*=cc;
  //Double_t f12=    dx*fP3*f1/(r1*r1*r1);
  Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1);  f14*=cc;
  //Double_t f13=    dx/r1;
  Double_t f24=    dx;                       f24*=cc;
  
  /* covariance matrix to be fixed!
  //b = C*ft
  Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
  Double_t b02=f24*fC40;
  Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
  Double_t b12=f24*fC41;
  Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
  Double_t b22=f24*fC42;
  Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
  Double_t b42=f24*fC44;
  Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
  Double_t b32=f24*fC43;
  
  //a = f*b = f*C*ft
  Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
  Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
  Double_t a22=f24*b42;

  //F*C*Ft = C + (b + bt + a)
  fC00 += b00 + b00 + a00;
  fC10 += b10 + b01 + a01; 
  fC20 += b20 + b02 + a02;
  fC30 += b30;
  fC40 += b40;
  fC11 += b11 + b11 + a11;
  fC21 += b21 + b12 + a12;
  fC31 += b31; 
  fC41 += b41;
  fC22 += b22 + b22 + a22;
  fC32 += b32;
  fC42 += b42;
  */
  
  Local2GlobalMomentum(localP, alpha); // convert local to global momentum
  SetP(localP);

  return kTRUE;
}
Commit	Line	Data
df9db588	1	/**************************************************************************
	2	* Copyright(c) 1998-2007, 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	/* $Id$ */
	17
	18	//-------------------------------------------------------------------------
0c5f89fb	19	// AOD track implementation of AliVParticle
df9db588	20	// Author: Markus Oldenburg, CERN
0c5f89fb	21	// Markus.Oldenburg@cern.ch
df9db588	22	//-------------------------------------------------------------------------
df9db588	23
8ac4fa64	24	#include "AliLog.h"
df9db588	25	#include "AliAODTrack.h"
	26
	27	ClassImp(AliAODTrack)
	28
	29	//______________________________________________________________________________
	30	AliAODTrack::AliAODTrack() :
9861edc0	31	AliVParticle(),
1912763f	32	fChi2perNDF(-999.),
9333290e	33	fChi2MatchTrigger(0.),
6efb741f	34	fFlags(0),
df9db588	35	fLabel(-999),
1912763f	36	fITSMuonClusterMap(0),
9333290e	37	fFilterMap(0),
02153d58	38	fID(-999),
9333290e	39	fCharge(-99),
e1c744ca	40	fType(kUndef),
9333290e	41	fCovMatrix(NULL),
7be1db84	42	fDetPid(NULL),
7be1db84	43	fProdVertex(NULL)
df9db588	44	{
	45	// default constructor
	46
	47	SetP();
	48	SetPosition((Float_t*)NULL);
	49	SetPID((Float_t*)NULL);
	50	}
	51
	52	//______________________________________________________________________________
02153d58	53	AliAODTrack::AliAODTrack(Short_t id,
df9db588	54	Int_t label,
	55	Double_t p[3],
	56	Bool_t cartesian,
	57	Double_t x[3],
	58	Bool_t isDCA,
	59	Double_t covMatrix[21],
	60	Short_t charge,
	61	UChar_t itsClusMap,
	62	Double_t pid[10],
	63	AliAODVertex *prodVertex,
1912763f	64	Bool_t usedForVtxFit,
dc825b15	65	Bool_t usedForPrimVtxFit,
ec40c484	66	AODTrk_t ttype,
ec40c484	67	UInt_t selectInfo) :
9861edc0	68	AliVParticle(),
1912763f	69	fChi2perNDF(-999.),
9333290e	70	fChi2MatchTrigger(0.),
6efb741f	71	fFlags(0),
df9db588	72	fLabel(label),
1912763f	73	fITSMuonClusterMap(itsClusMap),
9333290e	74	fFilterMap(selectInfo),
02153d58	75	fID(id),
9333290e	76	fCharge(charge),
e1c744ca	77	fType(ttype),
9333290e	78	fCovMatrix(NULL),
7be1db84	79	fDetPid(NULL),
9333290e	80	fProdVertex(prodVertex)
df9db588	81	{
	82	// constructor
	83
	84	SetP(p, cartesian);
	85	SetPosition(x, isDCA);
1912763f	86	SetUsedForVtxFit(usedForVtxFit);
dc825b15	87	SetUsedForPrimVtxFit(usedForPrimVtxFit);
df9db588	88	if(covMatrix) SetCovMatrix(covMatrix);
	89	SetPID(pid);
	90
	91	}
	92
	93	//______________________________________________________________________________
02153d58	94	AliAODTrack::AliAODTrack(Short_t id,
df9db588	95	Int_t label,
	96	Float_t p[3],
	97	Bool_t cartesian,
	98	Float_t x[3],
	99	Bool_t isDCA,
	100	Float_t covMatrix[21],
	101	Short_t charge,
	102	UChar_t itsClusMap,
	103	Float_t pid[10],
	104	AliAODVertex *prodVertex,
1912763f	105	Bool_t usedForVtxFit,
dc825b15	106	Bool_t usedForPrimVtxFit,
ec40c484	107	AODTrk_t ttype,
ec40c484	108	UInt_t selectInfo) :
9861edc0	109	AliVParticle(),
1912763f	110	fChi2perNDF(-999.),
9333290e	111	fChi2MatchTrigger(0.),
6efb741f	112	fFlags(0),
df9db588	113	fLabel(label),
1912763f	114	fITSMuonClusterMap(itsClusMap),
9333290e	115	fFilterMap(selectInfo),
02153d58	116	fID(id),
9333290e	117	fCharge(charge),
e1c744ca	118	fType(ttype),
9333290e	119	fCovMatrix(NULL),
7be1db84	120	fDetPid(NULL),
9333290e	121	fProdVertex(prodVertex)
df9db588	122	{
	123	// constructor
	124
	125	SetP(p, cartesian);
	126	SetPosition(x, isDCA);
1912763f	127	SetUsedForVtxFit(usedForVtxFit);
dc825b15	128	SetUsedForPrimVtxFit(usedForPrimVtxFit);
df9db588	129	if(covMatrix) SetCovMatrix(covMatrix);
df9db588	130	SetPID(pid);
df9db588	131	}
df9db588	132
df9db588	133	//______________________________________________________________________________
	134	AliAODTrack::~AliAODTrack()
	135	{
	136	// destructor
	137	delete fCovMatrix;
	138	}
	139
	140
	141	//______________________________________________________________________________
	142	AliAODTrack::AliAODTrack(const AliAODTrack& trk) :
9861edc0	143	AliVParticle(trk),
1912763f	144	fChi2perNDF(trk.fChi2perNDF),
9333290e	145	fChi2MatchTrigger(trk.fChi2MatchTrigger),
6efb741f	146	fFlags(trk.fFlags),
df9db588	147	fLabel(trk.fLabel),
1912763f	148	fITSMuonClusterMap(trk.fITSMuonClusterMap),
9333290e	149	fFilterMap(trk.fFilterMap),
02153d58	150	fID(trk.fID),
9333290e	151	fCharge(trk.fCharge),
e1c744ca	152	fType(trk.fType),
9333290e	153	fCovMatrix(NULL),
7be1db84	154	fDetPid(NULL),
9333290e	155	fProdVertex(trk.fProdVertex)
df9db588	156	{
	157	// Copy constructor
	158
	159	trk.GetP(fMomentum);
	160	trk.GetPosition(fPosition);
1912763f	161	SetUsedForVtxFit(trk.GetUsedForVtxFit());
dc825b15	162	SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
5d62ce04	163	if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
7be1db84	164	if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
df9db588	165	SetPID(trk.fPID);
df9db588	166	}
	167
	168	//______________________________________________________________________________
	169	AliAODTrack& AliAODTrack::operator=(const AliAODTrack& trk)
	170	{
	171	// Assignment operator
	172	if(this!=&trk) {
	173
9861edc0	174	AliVParticle::operator=(trk);
df9db588	175
	176	trk.GetP(fMomentum);
	177	trk.GetPosition(fPosition);
	178	trk.GetPID(fPID);
	179
1912763f	180	fChi2perNDF = trk.fChi2perNDF;
9333290e	181	fChi2MatchTrigger = trk.fChi2MatchTrigger;
df9db588	182
6efb741f	183	fFlags = trk.fFlags;
df9db588	184	fLabel = trk.fLabel;
df9db588	185
9333290e	186	fITSMuonClusterMap = trk.fITSMuonClusterMap;
	187	fFilterMap = trk.fFilterMap;
	188
6efb741f	189	fID = trk.fID;
6efb741f	190
9333290e	191	fCharge = trk.fCharge;
	192	fType = trk.fType;
	193
df9db588	194	delete fCovMatrix;
5d62ce04	195	if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
df9db588	196	else fCovMatrix=NULL;
	197	fProdVertex = trk.fProdVertex;
	198
1912763f	199	SetUsedForVtxFit(trk.GetUsedForVtxFit());
dc825b15	200	SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
7be1db84	201
	202	delete fDetPid;
	203	if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
	204	else fDetPid=NULL;
df9db588	205	}
	206
	207	return *this;
	208	}
	209
4697e4fb	210	//______________________________________________________________________________
	211	Double_t AliAODTrack::M(AODTrkPID_t pid) const
	212	{
	213	// Returns the mass.
9861edc0	214	// Masses for nuclei don't exist in the PDG tables, therefore they were put by hand.
4697e4fb	215
	216	switch (pid) {
	217
	218	case kElectron :
9861edc0	219	return 0.000510999; //TDatabasePDG::Instance()->GetParticle(11/::kElectron/)->Mass();
4697e4fb	220	break;
	221
	222	case kMuon :
9861edc0	223	return 0.1056584; //TDatabasePDG::Instance()->GetParticle(13/::kMuonMinus/)->Mass();
4697e4fb	224	break;
	225
	226	case kPion :
9861edc0	227	return 0.13957; //TDatabasePDG::Instance()->GetParticle(211/::kPiPlus/)->Mass();
4697e4fb	228	break;
	229
	230	case kKaon :
9861edc0	231	return 0.4937; //TDatabasePDG::Instance()->GetParticle(321/::kKPlus/)->Mass();
4697e4fb	232	break;
	233
	234	case kProton :
9861edc0	235	return 0.9382720; //TDatabasePDG::Instance()->GetParticle(2212/::kProton/)->Mass();
4697e4fb	236	break;
	237
	238	case kDeuteron :
9861edc0	239	return 1.8756; //TDatabasePDG::Instance()->GetParticle(1000010020)->Mass();
4697e4fb	240	break;
	241
	242	case kTriton :
9861edc0	243	return 2.8089; //TDatabasePDG::Instance()->GetParticle(1000010030)->Mass();
4697e4fb	244	break;
	245
	246	case kHelium3 :
9861edc0	247	return 2.8084; //TDatabasePDG::Instance()->GetParticle(1000020030)->Mass();
4697e4fb	248	break;
	249
	250	case kAlpha :
9861edc0	251	return 3.7274; //TDatabasePDG::Instance()->GetParticle(1000020040)->Mass();
4697e4fb	252	break;
	253
	254	case kUnknown :
	255	return -999.;
	256	break;
	257
	258	default :
	259	return -999.;
	260	}
	261	}
	262
	263	//______________________________________________________________________________
	264	Double_t AliAODTrack::E(AODTrkPID_t pid) const
	265	{
	266	// Returns the energy of the particle of a given pid.
	267
	268	if (pid != kUnknown) { // particle was identified
	269	Double_t m = M(pid);
	270	return TMath::Sqrt(P()P() + mm);
	271	} else { // pid unknown
	272	return -999.;
	273	}
	274	}
	275
	276	//______________________________________________________________________________
	277	Double_t AliAODTrack::Y(AODTrkPID_t pid) const
	278	{
9861edc0	279	// Returns the rapidity of a particle of a given pid.
4697e4fb	280
	281	if (pid != kUnknown) { // particle was identified
	282	Double_t e = E(pid);
	283	Double_t pz = Pz();
	284	if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
	285	return 0.5*TMath::Log((e+pz)/(e-pz));
	286	} else { // energy not known or equal to pz
	287	return -999.;
	288	}
	289	} else { // pid unknown
	290	return -999.;
	291	}
	292	}
	293
	294	//______________________________________________________________________________
	295	Double_t AliAODTrack::Y(Double_t m) const
	296	{
9861edc0	297	// Returns the rapidity of a particle of a given mass.
4697e4fb	298
	299	if (m >= 0.) { // mass makes sense
	300	Double_t e = E(m);
	301	Double_t pz = Pz();
	302	if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
	303	return 0.5*TMath::Log((e+pz)/(e-pz));
	304	} else { // energy not known or equal to pz
	305	return -999.;
	306	}
	307	} else { // pid unknown
	308	return -999.;
	309	}
	310	}
	311
	312	//______________________________________________________________________________
	313	AliAODTrack::AODTrkPID_t AliAODTrack::GetMostProbablePID() const
	314	{
	315	// Returns the most probable PID array element.
	316
	317	Int_t nPID = 10;
	318	if (fPID) {
	319	AODTrkPID_t loc = kUnknown;
	320	Double_t max = 0.;
	321	Bool_t allTheSame = kTRUE;
	322
	323	for (Int_t iPID = 0; iPID < nPID; iPID++) {
	324	if (fPID[iPID] >= max) {
	325	if (fPID[iPID] > max) {
	326	allTheSame = kFALSE;
	327	max = fPID[iPID];
	328	loc = (AODTrkPID_t)iPID;
	329	} else {
	330	allTheSame = kTRUE;
	331	}
	332	}
	333	}
	334
	335	return allTheSame ? kUnknown : loc;
	336	} else {
	337	return kUnknown;
	338	}
	339	}
	340
	341	//______________________________________________________________________________
	342	void AliAODTrack::ConvertAliPIDtoAODPID()
	343	{
	344	// Converts AliPID array.
	345	// The numbering scheme is the same for electrons, muons, pions, kaons, and protons.
	346	// Everything else has to be set to zero.
	347
	348	fPID[kDeuteron] = 0.;
8a1418dc	349	fPID[kTriton] = 0.;
	350	fPID[kHelium3] = 0.;
	351	fPID[kAlpha] = 0.;
	352	fPID[kUnknown] = 0.;
4697e4fb	353
	354	return;
	355	}
	356
	357
df9db588	358	//______________________________________________________________________________
	359	template <class T> void AliAODTrack::SetP(const T *p, const Bool_t cartesian)
	360	{
8a1418dc	361	// Set the momentum
df9db588	362
	363	if (p) {
	364	if (cartesian) {
16b65f2a	365	Double_t pt2 = p[0]p[0] + p[1]p[1];
0c5f89fb	366	Double_t pp = TMath::Sqrt(pt2 + p[2]*p[2]);
df9db588	367
16b65f2a	368	fMomentum[0] = TMath::Sqrt(pt2); // pt
b1a9edc8	369	fMomentum[1] = (pt2 != 0.) ? TMath::Pi()+TMath::ATan2(-p[1], -p[0]) : -999; // phi
0c5f89fb	370	fMomentum[2] = (pp != 0.) ? TMath::ACos(p[2] / pp) : -999.; // theta
df9db588	371	} else {
16b65f2a	372	fMomentum[0] = p[0]; // pt
df9db588	373	fMomentum[1] = p[1]; // phi
	374	fMomentum[2] = p[2]; // theta
	375	}
	376	} else {
	377	fMomentum[0] = -999.;
	378	fMomentum[1] = -999.;
	379	fMomentum[2] = -999.;
	380	}
	381	}
	382
	383	//______________________________________________________________________________
	384	template <class T> void AliAODTrack::SetPosition(const T *x, const Bool_t dca)
	385	{
	386	// set the position
	387
	388	if (x) {
	389	if (!dca) {
	390	ResetBit(kIsDCA);
	391
	392	fPosition[0] = x[0];
	393	fPosition[1] = x[1];
	394	fPosition[2] = x[2];
	395	} else {
	396	SetBit(kIsDCA);
	397	// don't know any better yet
	398	fPosition[0] = -999.;
	399	fPosition[1] = -999.;
	400	fPosition[2] = -999.;
	401	}
	402	} else {
	403	ResetBit(kIsDCA);
	404
	405	fPosition[0] = -999.;
	406	fPosition[1] = -999.;
	407	fPosition[2] = -999.;
	408	}
	409	}
	410
	411	//______________________________________________________________________________
	412	void AliAODTrack::SetDCA(Double_t d, Double_t z)
	413	{
	414	// set the dca
	415	fPosition[0] = d;
	416	fPosition[1] = z;
	417	fPosition[2] = 0.;
	418	SetBit(kIsDCA);
	419	}
	420
	421	//______________________________________________________________________________
	422	void AliAODTrack::Print(Option_t* /* option */) const
	423	{
	424	// prints information about AliAODTrack
	425
	426	printf("Object name: %s Track type: %s\n", GetName(), GetTitle());
	427	printf(" px = %f\n", Px());
	428	printf(" py = %f\n", Py());
	429	printf(" pz = %f\n", Pz());
	430	printf(" pt = %f\n", Pt());
	431	printf(" 1/pt = %f\n", OneOverPt());
	432	printf(" theta = %f\n", Theta());
	433	printf(" phi = %f\n", Phi());
1912763f	434	printf(" chi2/NDF = %f\n", Chi2perNDF());
df9db588	435	printf(" charge = %d\n", Charge());
df9db588	436	}
df9db588	437
8ac4fa64	438	void AliAODTrack::SetMatchTrigger(Int_t matchTrig){
8a1418dc	439	//
8a1418dc	440	// Set the MUON trigger information
8ac4fa64	441	switch(matchTrig){
e1c744ca	442	case 0: // 0 track does not match trigger
	443	fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
	444	break;
	445	case 1: // 1 track match but does not pass pt cut
	446	fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)\|0x40000000;
	447	break;
	448	case 2: // 2 track match Low pt cut
	449	fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)\|0x80000000;
	450	break;
	451	case 3: // 3 track match High pt cut
	452	fITSMuonClusterMap=fITSMuonClusterMap\|0xc0000000;
	453	break;
	454	default:
	455	fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
8ac4fa64	456	AliWarning(Form("unknown case for matchTrig: %d\n",matchTrig));
e1c744ca	457	}
	458	}
	459
	460	void AliAODTrack::SetHitsPatternInTrigCh(UShort_t hitsPatternInTrigCh){
8a1418dc	461	//
8a1418dc	462	// Set the MUON hit pattern (1 bit per chamber)
e1c744ca	463	fITSMuonClusterMap=(fITSMuonClusterMap&0xffff00ff)\|(hitsPatternInTrigCh<<8);
	464	}
	465
	466	Int_t AliAODTrack::HitsMT(Int_t istation, Int_t iplane, Char_t *cathode){
8a1418dc	467	//
8a1418dc	468	// Retrieve hit information for MUON identified by (station, plane, cathode)
e1c744ca	469	if(cathode){
	470	if(cathode[0]=='x'\|\|cathode[0]=='X'){
	471	if(istation==1){
	472	if(iplane==1)
	473	return (fITSMuonClusterMap&0x8000)?1:0;
	474	else if(iplane==2)
	475	return (fITSMuonClusterMap&0x4000)?1:0;
	476	else
	477	return 0;
	478	}else if(istation==2){
	479	if(iplane==1)
	480	return (fITSMuonClusterMap&0x2000)?1:0;
	481	else if(iplane==2)
	482	return (fITSMuonClusterMap&0x1000)?1:0;
	483	else
	484	return 0;
	485	}else{
	486	return 0;
	487	}
	488	}else if(cathode[0]=='y'\|\|cathode[0]=='Y'){
	489	if(istation==1){
	490	if(iplane==1)
	491	return (fITSMuonClusterMap&0x0800)?1:0;
	492	else if(iplane==2)
	493	return (fITSMuonClusterMap&0x0400)?1:0;
	494	else
	495	return 0;
	496	}else if(istation==2){
	497	if(iplane==1)
	498	return (fITSMuonClusterMap&0x0200)?1:0;
	499	else if(iplane==2)
	500	return (fITSMuonClusterMap&0x0100)?1:0;
	501	else
	502	return 0;
	503	}else{
	504	return 0;
	505	}
	506	}else{
	507	return 0;
	508	}
	509	}else{
	510	if(istation==1){
	511	if(iplane==1)
	512	return (HitsMT(1,1,"X")\|\|HitsMT(1,1,"Y"))?1:0;
	513	else if(iplane==2)
	514	return (HitsMT(1,2,"X")\|\|HitsMT(1,2,"Y"))?1:0;
	515	else
	516	return 0;
	517	}else if(istation==2){
	518	if(iplane==1)
	519	return (HitsMT(2,1,"X")\|\|HitsMT(2,1,"Y"))?1:0;
	520	else if(iplane==2)
	521	return (HitsMT(2,2,"X")\|\|HitsMT(2,2,"Y"))?1:0;
	522	else
	523	return 0;
	524	}else{
	525	return 0;
	526	}
	527	}
	528	}
	529
	530	Int_t AliAODTrack::HitsMuonChamber(Int_t MuonChamber){
8a1418dc	531	// Retrieve hit information for MUON Chamber
e1c744ca	532	switch(MuonChamber){
	533	case 11:
	534	return HitsMT(1,1);
	535	case 12:
	536	return HitsMT(1,2);
	537	case 13:
	538	return HitsMT(2,1);
	539	case 14:
	540	return HitsMT(2,2);
	541	default:
	542	printf("Unknown MUON chamber: %d\n",MuonChamber);
	543	return 0;
	544	}
	545	}
c683ddc2	546
	547
	548
	549	Bool_t AliAODTrack::PropagateTo(Double_t xk, Double_t b) {
	550	//----------------------------------------------------------------
	551	// Propagate this track to the plane X=xk (cm) in the field "b" (kG)
	552	// This is in local coordinates!!!
	553	//----------------------------------------------------------------
	554
	555	Double_t alpha = 0.;
	556	Double_t localP[3] = {Px(), Py(), Pz()}; // set global (sic!) p
	557	Global2LocalMomentum(localP, Charge(), alpha); // convert global to local momentum
	558
	559	AliAODVertex origin = (AliAODVertex)fProdVertex.GetObject();
	560	Double_t localX[3] = {origin->GetX(), origin->GetY(), origin->GetZ()}; // set global (sic!) location of first track point
	561	Global2LocalPosition(localX, alpha); // convert global to local position
	562
	563	Double_t &fX = localX[0];
	564
	565	Double_t dx=xk-fX;
	566	if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
	567
	568	Double_t crv=localP[0]bkB2C;
	569	if (TMath::Abs(b) < kAlmost0Field) crv=0.;
	570
	571	Double_t f1=localP[1], f2=f1 + crv*dx;
	572	if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
	573	if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
	574
	575	Double_t &fP0=localX[1], &fP1=localX[2], &fP2=localP[0], &fP3=localP[1], &fP4=localP[2];
	576	/* covariance matrix to be fixed!
	577	Double_t
	578	&fC00=fC[0],
	579	&fC10=fC[1], &fC11=fC[2],
	580	&fC20=fC[3], &fC21=fC[4], &fC22=fC[5],
	581	&fC30=fC[6], &fC31=fC[7], &fC32=fC[8], &fC33=fC[9],
	582	&fC40=fC[10], &fC41=fC[11], &fC42=fC[12], &fC43=fC[13], &fC44=fC[14];
	583	*/
	584	Double_t r1=TMath::Sqrt(1.- f1f1), r2=TMath::Sqrt(1.- f2f2);
	585
	586	fX=xk;
	587	fP0 += dx*(f1+f2)/(r1+r2);
	588	fP1 += dx(r2 + f2(f1+f2)/(r1+r2))*fP3;
	589	fP2 += dx*crv;
	590
	591	//f = F - 1
	592
d067f357	593	//Double_t f02= dx/(r1r1r1);
d067f357	594	Double_t cc=crv/fP4;
c683ddc2	595	Double_t f04=0.5dxdx/(r1r1r1); f04*=cc;
d067f357	596	//Double_t f12= dxfP3f1/(r1r1r1);
c683ddc2	597	Double_t f14=0.5dxdxfP3f1/(r1r1r1); f14*=cc;
d067f357	598	//Double_t f13= dx/r1;
c683ddc2	599	Double_t f24= dx; f24*=cc;
	600
	601	/* covariance matrix to be fixed!
	602	//b = C*ft
	603	Double_t b00=f02fC20 + f04fC40, b01=f12fC20 + f14fC40 + f13*fC30;
	604	Double_t b02=f24*fC40;
	605	Double_t b10=f02fC21 + f04fC41, b11=f12fC21 + f14fC41 + f13*fC31;
	606	Double_t b12=f24*fC41;
	607	Double_t b20=f02fC22 + f04fC42, b21=f12fC22 + f14fC42 + f13*fC32;
	608	Double_t b22=f24*fC42;
	609	Double_t b40=f02fC42 + f04fC44, b41=f12fC42 + f14fC44 + f13*fC43;
	610	Double_t b42=f24*fC44;
	611	Double_t b30=f02fC32 + f04fC43, b31=f12fC32 + f14fC43 + f13*fC33;
	612	Double_t b32=f24*fC43;
	613
	614	//a = fb = fC*ft
	615	Double_t a00=f02b20+f04b40,a01=f02b21+f04b41,a02=f02b22+f04b42;
	616	Double_t a11=f12b21+f14b41+f13b31,a12=f12b22+f14b42+f13b32;
	617	Double_t a22=f24*b42;
	618
	619	//FCFt = C + (b + bt + a)
	620	fC00 += b00 + b00 + a00;
	621	fC10 += b10 + b01 + a01;
	622	fC20 += b20 + b02 + a02;
	623	fC30 += b30;
	624	fC40 += b40;
	625	fC11 += b11 + b11 + a11;
	626	fC21 += b21 + b12 + a12;
	627	fC31 += b31;
	628	fC41 += b41;
	629	fC22 += b22 + b22 + a22;
	630	fC32 += b32;
	631	fC42 += b42;
	632	*/
	633
	634	Local2GlobalMomentum(localP, alpha); // convert local to global momentum
	635	SetP(localP);
	636
	637	return kTRUE;
	638	}