[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);
  SetXYAtDCA(-999., -999.);
  SetPxPyPzAtDCA(-999., -999., -999.);
  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(0),
  fFilterMap(selectInfo),
  fID(id),
  fCharge(charge),
  fType(ttype),
  fCovMatrix(NULL),
  fDetPid(NULL),
  fProdVertex(prodVertex)
{
  // constructor
 
  SetP(p, cartesian);
  SetPosition(x, isDCA);
  SetXYAtDCA(-999., -999.);
  SetPxPyPzAtDCA(-999., -999., -999.);
  SetUsedForVtxFit(usedForVtxFit);
  SetUsedForPrimVtxFit(usedForPrimVtxFit);
  if(covMatrix) SetCovMatrix(covMatrix);
  SetPID(pid);
  SetITSClusterMap(itsClusMap);
}

//______________________________________________________________________________
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(0),
  fFilterMap(selectInfo),
  fID(id),
  fCharge(charge),
  fType(ttype),
  fCovMatrix(NULL),
  fDetPid(NULL),
  fProdVertex(prodVertex)
{
  // constructor
 
  SetP(p, cartesian);
  SetPosition(x, isDCA);
  SetXYAtDCA(-999., -999.);
  SetPxPyPzAtDCA(-999., -999., -999.);
  SetUsedForVtxFit(usedForVtxFit);
  SetUsedForPrimVtxFit(usedForPrimVtxFit);
  if(covMatrix) SetCovMatrix(covMatrix);
  SetPID(pid);
  SetITSClusterMap(itsClusMap);
}

//______________________________________________________________________________
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);
  SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
  SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
  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);

    SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
    SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
    
    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));
  }
}

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