[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;
  delete fDetPid;
}


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