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