[u/mrichter/AliRoot.git] / PWGPP / TRD / info / AliTRDtrackInfo.cxx

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

/* $Id: AliTRDtrackInfo.cxx 27496 2008-07-22 08:35:45Z cblume $ */

////////////////////////////////////////////////////////////////////////////
//                                                                        //
//  Reconstruction QA                                                     //
//                                                                        //
//  Authors:                                                              //
//    Alex Bercuci <A.Bercuci@gsi.de>                                     //
//    Markus Fasel <M.Fasel@gsi.de>                                       //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

#include "TDatabasePDG.h"
#include "TPDGCode.h"
#include "TVectorT.h"

#include "AliTrackReference.h"
#include "AliTrackPointArray.h"
#include "AliExternalTrackParam.h"

#include "AliTRDseedV1.h"
#include "AliTRDtrackV1.h"
#include "AliTRDgeometry.h"
#include "AliTRDtrackerV1.h"

#include "AliTRDtrackInfo.h"

ClassImp(AliTRDtrackInfo)
ClassImp(AliTRDtrackInfo::AliMCinfo)
ClassImp(AliTRDtrackInfo::AliESDinfo)
Double_t AliTRDtrackInfo::AliMCinfo::fgKalmanStep = 2.;
Bool_t AliTRDtrackInfo::AliMCinfo::fgKalmanUpdate = kTRUE;

//___________________________________________________
AliTRDtrackInfo::AliTRDtrackInfo():
  TObject()
  ,fNClusters(0)
  ,fTRDtrack(NULL)
  ,fMC(NULL)
  ,fESD()
{
  //
  // Default constructor
  //
}


//___________________________________________________
AliTRDtrackInfo::AliTRDtrackInfo(const AliTRDtrackInfo &trdInfo):
  TObject((const TObject&)trdInfo)  
  ,fNClusters(trdInfo.fNClusters)
  ,fTRDtrack(NULL)
  ,fMC(NULL)
  ,fESD(trdInfo.fESD)
{
  //
  // copy Entries
  //

  if(trdInfo.fMC) fMC = new AliMCinfo(*trdInfo.fMC);
  SetTrack(trdInfo.fTRDtrack);
}

//___________________________________________________
AliTRDtrackInfo::AliMCinfo::AliMCinfo()
  :fLabel(0)
  ,fTRDlabel(0)
  ,fPDG(0)
  ,fNTrackRefs(0)
  ,fEta(-999.)
  ,fPhi(-999.)
  ,fPt(-1.)
{
  // Set 0-Pointers
  memset(fTrackRefs, 0, sizeof(AliTrackReference *) * 12);
}

//___________________________________________________
AliTRDtrackInfo::AliMCinfo::AliMCinfo(const AliMCinfo &mc)
  :fLabel(mc.fLabel)
  ,fTRDlabel(mc.fTRDlabel)
  ,fPDG(mc.fPDG)
  ,fNTrackRefs(mc.fNTrackRefs)
  ,fEta(mc.fEta)
  ,fPhi(mc.fPhi)
  ,fPt(mc.fPt)
{
  //
  // Constructor
  //

  memset(fTrackRefs, 0, sizeof(AliTrackReference *) * 12);
  for(Int_t ien = 0; ien < 12; ien++){
    if(mc.fTrackRefs[ien])
      fTrackRefs[ien] = new AliTrackReference(*(mc.fTrackRefs[ien]));
  }
}

//________________________________________________________
Int_t AliTRDtrackInfo::AliMCinfo::GetPID() const
{
// Translate pdg code to PID index

  switch(fPDG){
  case kElectron: 
  case kPositron: return AliPID::kElectron;  
  case kMuonPlus:
  case kMuonMinus: return AliPID::kMuon;  
  case kPiPlus: 
  case kPiMinus: return AliPID::kPion;  
  case kKPlus: 
  case kKMinus: return AliPID::kKaon;
  case kProton: 
  case kProtonBar: return AliPID::kProton;
  } 
  return -1;
}

//___________________________________________________
AliTRDtrackInfo::AliESDinfo::AliESDinfo()
  :fSteer(0)
  ,fId(-1)
  ,fStatus(0)
  ,fKinkIndex(0)
  ,fTPCncls(0)
  ,fTPCdedx(0.)
  ,fTOFbeta(0.)
  ,fTOFbc(0)
  ,fTRDpidQuality(0)
  ,fTRDnSlices(0)
  ,fPt(0.)
  ,fPhi(-999.)
  ,fEta(-999.)
  ,fTRDslices(NULL)
  ,fOP(NULL)
  ,fTPCout(NULL)
  ,fITSout(NULL)
  ,fTPArray(NULL)
{
  //
  // Constructor
  //

  memset(fTRDr, 0, AliPID::kSPECIES*sizeof(Double32_t));
  memset(fTRDv0pid, 0, AliPID::kSPECIES*sizeof(Int_t));
}

//___________________________________________________
AliTRDtrackInfo::AliESDinfo::AliESDinfo(const AliESDinfo &esd)
  :fSteer(esd.fSteer)
  ,fId(esd.fId)
  ,fStatus(esd.fStatus)
  ,fKinkIndex(esd.fKinkIndex)
  ,fTPCncls(esd.fTPCncls)
  ,fTPCdedx(esd.fTPCdedx)
  ,fTOFbeta(esd.fTOFbeta)
  ,fTOFbc(esd.fTOFbc)
  ,fTRDpidQuality(esd.fTRDpidQuality)
  ,fTRDnSlices(esd.fTRDnSlices)
  ,fPt(esd.fPt)
  ,fPhi(esd.fPhi)
  ,fEta(esd.fEta)
  ,fTRDslices(NULL)
  ,fOP(NULL)
  ,fTPCout(NULL)
  ,fITSout(NULL)
  ,fTPArray(NULL)
{
  //
  // Constructor
  //

  memcpy(fTRDr, esd.fTRDr, AliPID::kSPECIES*sizeof(Double32_t));
  memcpy(fTRDv0pid, esd.fTRDv0pid, AliPID::kSPECIES*sizeof(Int_t));

  if(fTRDnSlices){
    fTRDslices = new Double32_t[fTRDnSlices];
    memcpy(fTRDslices, esd.fTRDslices, fTRDnSlices*sizeof(Double32_t));
  }
  SetOuterParam(esd.fOP);
  SetTPCoutParam(esd.fTPCout);
  SetITSoutParam(esd.fITSout);
  SetTrackPointArray(esd.fTPArray);
}


//___________________________________________________
AliTRDtrackInfo::~AliTRDtrackInfo()
{
  //
  // Destructor
  //

  if(fMC) delete fMC;
  if(fTRDtrack) delete fTRDtrack;
}

//___________________________________________________
AliTRDtrackInfo::AliMCinfo::~AliMCinfo()
{
  //
  // Destructor
  //

  fNTrackRefs = 0;
  for(Int_t ien = 0; ien < 12; ien++){
    if(fTrackRefs[ien]) delete fTrackRefs[ien];
    fTrackRefs[ien] = NULL;
  }
}

//___________________________________________________
AliTRDtrackInfo::AliESDinfo::~AliESDinfo()
{
  //
  // Destructor
  //

  if(fTRDnSlices){
    delete [] fTRDslices;
    fTRDslices = NULL;
    fTRDnSlices = 0;
  }
  if(fOP) delete fOP; fOP = NULL;
  if(fTPCout) delete fTPCout; fTPCout = NULL;
  if(fITSout) delete fITSout; fITSout = NULL;
  if(fTPArray) delete fTPArray;
}

//___________________________________________________
void AliTRDtrackInfo::AliESDinfo::Delete(const Option_t *){
  //
  // Delete Pointer members 
  // 
  if(fTRDnSlices){
    delete [] fTRDslices;
    fTRDslices = NULL;
    fTRDnSlices = 0;
  }
  if(fOP) delete fOP; fOP = NULL;
  if(fTPCout) delete fTPCout; fTPCout = NULL;
  if(fITSout) delete fITSout; fITSout = NULL;
  if(fTPArray) delete fTPArray; fTPArray = NULL;
}


//___________________________________________________
AliTRDtrackInfo& AliTRDtrackInfo::operator=(const AliTRDtrackInfo &trdInfo)
{
  //
  // = Operator
  //
  if(this == &trdInfo) return *this;

  fNClusters  = trdInfo.fNClusters;
  fESD        = trdInfo.fESD;

  if(trdInfo.fMC){
    if(!fMC) fMC = new AliMCinfo(*trdInfo.fMC);
    else{
      fMC->~AliMCinfo();
      new(fMC) AliMCinfo(*trdInfo.fMC);
    }
  } else {if(fMC) delete fMC; fMC = NULL;}

  SetTrack(trdInfo.fTRDtrack);

  return *this;
}

//___________________________________________________
AliTRDtrackInfo::AliMCinfo& AliTRDtrackInfo::AliMCinfo::operator=(const AliMCinfo &mc)
{
  //
  // Assignment operator
  //

  if(this == &mc) return *this;
  fLabel      = mc.fLabel;
  fTRDlabel   = mc.fTRDlabel;
  fPDG        = mc.fPDG;
  fNTrackRefs = mc.fNTrackRefs;

  AliTrackReference **itr = &fTrackRefs[0];
  AliTrackReference* const *jtr = &mc.fTrackRefs[0];
  for(Int_t ien = 0; ien < 2*AliTRDgeometry::kNlayer; ien++, itr++, jtr++){
    if((*jtr)){
      if(!(*itr)) (*itr) = new AliTrackReference(*(*jtr));
      else{
        (*itr)->~AliTrackReference();
        new(&(*itr)) AliTrackReference(*(*jtr));
      }
    } else {if((*itr)) delete (*itr); (*itr) = NULL;}
  }
  return *this;
}

//___________________________________________________
AliTRDtrackInfo::AliESDinfo& AliTRDtrackInfo::AliESDinfo::operator=(const AliESDinfo &esd)
{
  //
  // Assignment operator
  //

  if(this == &esd) return *this;
  fSteer       = esd.fSteer;
  fId          = esd.fId;
  fStatus      = esd.fStatus;
  fKinkIndex   = esd.fKinkIndex;
  fTPCncls     = esd.fTPCncls;
  fTPCdedx     = esd.fTPCdedx;
  fTOFbeta     = esd.fTOFbeta;
  fTOFbc       = esd.fTOFbc;
  fTRDpidQuality= esd.fTRDpidQuality;
  fTRDnSlices  = esd.fTRDnSlices;
  fPt          = esd.fPt;
  fPhi         = esd.fPhi;
  fEta         = esd.fEta;
  
  memcpy(fTRDr, esd.fTRDr, AliPID::kSPECIES*sizeof(Double32_t));
  memcpy(fTRDv0pid, esd.fTRDv0pid, AliPID::kSPECIES*sizeof(Int_t));

  if(fTRDnSlices){
    if(!fTRDslices) fTRDslices = new Double32_t[fTRDnSlices];
    memcpy(fTRDslices, esd.fTRDslices, fTRDnSlices*sizeof(Double32_t));
  }
  SetOuterParam(esd.fOP);
  SetTPCoutParam(esd.fTPCout);
  SetITSoutParam(esd.fITSout);
  SetTrackPointArray(esd.fTPArray);

  return *this;
}

//___________________________________________________
void AliTRDtrackInfo::Delete(const Option_t *)
{
  //
  // Delete
  //

  AliDebug(2, Form("track[%p] mc[%p]", (void*)fTRDtrack, (void*)fMC));
  fNClusters  = 0;
  if(fMC) delete fMC; fMC = NULL;
  fESD.Delete(NULL);
  if(fTRDtrack) delete fTRDtrack; fTRDtrack = NULL;
}

//___________________________________________________
void AliTRDtrackInfo::SetTrack(const AliTRDtrackV1 *track)
{
  //
  // Set the TRD track
  //

  if(track){
    if(!fTRDtrack) fTRDtrack = new AliTRDtrackV1(*track);
    else{
      fTRDtrack->~AliTRDtrackV1();
      new(fTRDtrack) AliTRDtrackV1(*track);
    }
    if(track->IsOwner()) fTRDtrack->SetOwner();
  } else {
    if(fTRDtrack) delete fTRDtrack; fTRDtrack = NULL;
  }
}

//___________________________________________________
void AliTRDtrackInfo::AliESDinfo::SetTrackPointArray(const AliTrackPointArray *tps)
{
  //
  // Set the track point array for alignment task
  //


  if(tps){
    if(!fTPArray) fTPArray = new AliTrackPointArray(*tps);
    else{
      fTPArray->~AliTrackPointArray();
      new(fTPArray) AliTrackPointArray(*tps);
    }
  } else {
    if(fTPArray) delete fTPArray; fTPArray = NULL;
  }
}

//___________________________________________________
void AliTRDtrackInfo::AddTrackRef(const AliTrackReference *tref)
{
  //
  // Add track reference
  //

  if(fMC->fNTrackRefs >= 2*AliTRDgeometry::kNlayer){ 
    SetCurved();
    return;
  }
  // Make a copy for the object in order to avoid ownership problems
  fMC->fTrackRefs[fMC->fNTrackRefs++] = new AliTrackReference(*tref);
}

//___________________________________________________
AliTrackReference* AliTRDtrackInfo::GetTrackRef(Int_t idx) const
{
//
// Returns a track reference
//
  if(!fMC) return NULL;
  return (idx>=0 && idx < 12) ? fMC->fTrackRefs[idx] : NULL;
}

//___________________________________________________
AliTrackReference* AliTRDtrackInfo::GetTrackRef(const AliTRDseedV1* const tracklet) const
{
//
// Returns a track reference
//
  if(!fMC) return NULL;
  Double_t cw = AliTRDgeometry::CamHght() + AliTRDgeometry::CdrHght();
  AliTrackReference * const* jtr = &(fMC->fTrackRefs[0]);
  for(Int_t itr = 0; itr < fMC->fNTrackRefs; itr++, ++jtr){
    if(!(*jtr)) break;   
    if(TMath::Abs(tracklet->GetX0() - (*jtr)->LocalX()) < cw) return (*jtr);
  }
  return NULL;
}

//___________________________________________________
Int_t AliTRDtrackInfo::GetNumberOfClusters() const
{
  //
  // Returns the number of clusters
  //

  Int_t n = 0;
  if(!fTRDtrack) return 0;
  if(fTRDtrack->GetNumberOfTracklets() == 0) return n;
  AliTRDseedV1 *tracklet = NULL;
  for(Int_t ip=0; ip<AliTRDgeometry::kNlayer; ip++){
    if(!(tracklet = const_cast<AliTRDseedV1 *>(fTRDtrack->GetTracklet(ip)))) continue;
    n+=tracklet->GetN();
  }
  return n;
}


//___________________________________________________
void  AliTRDtrackInfo::AliESDinfo::SetOuterParam(const AliExternalTrackParam *op)
{
  //
  // Set outer track parameters
  //

  if(op){
    if(fOP){
      fOP->~AliExternalTrackParam();
      // RS: Constructor from VTrack was used instead of Constructor from AliExternalTrackParam
      new(fOP) AliExternalTrackParam(*op);
    } else fOP = new AliExternalTrackParam(*op);
  } else {
    if(fOP) delete fOP; fOP = NULL;
  }
}

//___________________________________________________
void  AliTRDtrackInfo::AliESDinfo::SetITSoutParam(const AliExternalTrackParam *op)
{
  //
  // Set TPCout track parameters
  //

  if(op){
    if(fITSout){
      fITSout->~AliExternalTrackParam();
      // RS: Constructor from VTrack was used instead of Constructor from AliExternalTrackParam
      new(fITSout) AliExternalTrackParam(*op);
    } else fITSout = new AliExternalTrackParam(*op);
  } else {
    if(fITSout) delete fITSout; fITSout = NULL;
  }
}

//___________________________________________________
void  AliTRDtrackInfo::AliESDinfo::SetTPCoutParam(const AliExternalTrackParam *op)
{
  //
  // Set TPCout track parameters
  //

  if(op){
    if(fTPCout){
      fTPCout->~AliExternalTrackParam();
      // RS: Constructor from VTrack was used instead of Constructor from AliExternalTrackParam
      new(fTPCout) AliExternalTrackParam(*op);
    } else fTPCout = new AliExternalTrackParam(*op);
  } else {
    if(fTPCout) delete fTPCout; fTPCout = NULL;
  }
}

//___________________________________________________
Int_t AliTRDtrackInfo::GetNTracklets() const
{
  //
  // Return the number of tracklets
  //

  if(!fTRDtrack) return 0;
  return fTRDtrack->GetNumberOfTracklets();
}

//___________________________________________________
void AliTRDtrackInfo::SetSlices(Int_t n, Double32_t *s)
{
  //
  // Set the slices
  //
  if(fESD.fTRDnSlices != n){
    fESD.fTRDnSlices = 0;
    delete [] fESD.fTRDslices;
    fESD.fTRDslices = NULL;
  }

  if(!fESD.fTRDnSlices){
    fESD.fTRDnSlices = n;
    fESD.fTRDslices = new Double32_t[fESD.fTRDnSlices];
  }

  memcpy(fESD.fTRDslices, s, n*sizeof(Double32_t));
}
 
//___________________________________________________
Bool_t AliTRDtrackInfo::AliMCinfo::GetDirections(Float_t &x0, Float_t &y0, Float_t &z0, Float_t &dydx, Float_t &dzdx, Float_t &pt, Float_t &p, Float_t &eta, Float_t &phi, UChar_t &status) const
{
// Check for 2 track ref for the tracklet defined bythe radial position x0
// The "status" is a bit map and gives a more informative output in case of failure:
//   - 0 : everything is OK
//   - BIT(0) : 0 track refs found
//   - BIT(1) : 1 track reference found
//   - BIT(2) : dx <= 0 between track references
//   - BIT(3) : dx > 0 && dx < 3.7 - tangent tracks 

  status = 0;
  Double_t cw = AliTRDgeometry::CamHght() + AliTRDgeometry::CdrHght();
  Int_t nFound = 0;
  AliTrackReference *tr[2] = {NULL, NULL};
  AliTrackReference * const* jtr = &fTrackRefs[0];
  for(Int_t itr = 0; itr < fNTrackRefs; itr++, ++jtr){
    if(!(*jtr)) break;
/*
    if(fDebugLevel>=5) printf("\t\tref[%2d] x[%6.3f]\n", itr, (*jtr)->LocalX());*/
    if(TMath::Abs(x0 - (*jtr)->LocalX()) > cw) continue;
    tr[nFound++] = (*jtr);
    if(nFound == 2) break;
  } 
  if(nFound < 2){ 
    //AliDebug(1, Form("Missing track ref x0[%6.3f] nref[%d]", x0, nFound));
    if(!nFound) SETBIT(status, 0);
    else SETBIT(status, 1);
    return kFALSE;
  }
  pt=tr[1]->Pt();
  p =tr[1]->P();
  if(pt < 1.e-3) return kFALSE;

  Double_t dx = tr[1]->LocalX() - tr[0]->LocalX();
  if(dx <= 0.){
    AliWarningGeneral("AliTRDtrackInfo::AliMCinfo::GetDirections()", Form("Track ref with wrong radial distances refX0[%6.3f] refX1[%6.3f]", tr[0]->LocalX(), tr[1]->LocalX()));
    SETBIT(status, 2);
    return kFALSE;
  }
  if(TMath::Abs(dx-AliTRDgeometry::CamHght()-AliTRDgeometry::CdrHght())>1.E-3) SETBIT(status, 3); 

  dydx = (tr[1]->LocalY() - tr[0]->LocalY()) / dx;
  dzdx = (tr[1]->Z() - tr[0]->Z()) / dx;
  //Float_t dx0 = tr[1]->LocalX() - x0;
  y0   =  tr[1]->LocalY()/* - dydx*dx0*/;
  z0   =  tr[1]->Z()/* - dzdx*dx0*/;
  x0   =  tr[1]->LocalX();
  eta  =  -TMath::Log(TMath::Tan(0.5 * tr[1]->Theta()));
  phi  =  TMath::ATan2(tr[1]->Y(), tr[1]->X());
  return kTRUE;
}

//___________________________________________________
Bool_t AliTRDtrackInfo::AliMCinfo::PropagateKalman(
      TVectorD *x, TVectorD *y, TVectorD *z,
      TVectorD *dx, TVectorD *dy, TVectorD *dz,
      TVectorD *pt, TVectorD *dpt, TVectorD *budget, TVectorD *c, Double_t mass) const
{
// Propagate Kalman from the first TRD track reference to 
// last one and save residuals in the y, z and pt.
// 
// This is to calibrate the dEdx and MS corrections

  if(!fNTrackRefs) return kFALSE;
  for(Int_t itr=kNTrackRefs; itr--;){
    (*x)[itr] = 0.;(*y)[itr] = 0.;(*z)[itr] = 0.;
    (*dx)[itr] = -1.; (*dy)[itr] = 100.; (*dz)[itr] = 100.; (*dpt)[itr] = 100.;
  }

  // Initialize TRD track to the first track reference
  AliTrackReference *tr(NULL);
  Int_t itr(0); while(!(tr = fTrackRefs[itr])) itr++;
  if(tr->Pt()<1.e-3) return kFALSE;

  AliTRDtrackV1 tt;
  Double_t xyz[3]={tr->X(),tr->Y(),tr->Z()};
  Double_t pxyz[3]={tr->Px(),tr->Py(),tr->Pz()};
  Double_t var[6] = {1.e-4, 1.e-4, 1.e-4, 1.e-4, 1.e-4, 1.e-4};
  Double_t cov[21]={
    var[0],  0.,  0.,  0.,  0.,  0.,
         var[1],  0.,  0.,  0.,  0.,
              var[2],  0.,  0.,  0.,
                   var[3],  0.,  0.,
                        var[4],  0.,
                             var[5]
  };
  TDatabasePDG db;
  const TParticlePDG *pdg=db.GetParticle(fPDG);
  if(!pdg){
    AliWarningGeneral("AliTRDtrackInfo::AliMCinfo::PropagateKalman()", Form("PDG entry missing for code %d. References for track %d", fPDG, fNTrackRefs));
    return kFALSE;
  }
  tt.Set(xyz, pxyz, cov, Short_t(pdg->Charge()));
  if(mass<0){ // mass 0 use PDG
    tt.SetMass(pdg->Mass());
  } else { // use rec value
    tt.SetMass(mass);
  }

//  Double_t bg(tr->P()/pdg->Mass());
//  printf("\n\nNEW track PDG[%d] bg[%f] x[%f]\n", fPDG, bg, TMath::Log(bg));
  Double_t x0(tr->LocalX());
  const Double_t *cc(NULL);
  for(Int_t ip=0; itr<fNTrackRefs; itr++){
    if(!(tr = fTrackRefs[itr])) continue;
//    printf("ip[%d] det[%d]\n", ip, tr->DetectorId());
    if(!AliTRDtrackerV1::PropagateToX(tt, tr->LocalX(), fgKalmanStep)) continue;

    //if(update) ...
    tt.GetXYZ(xyz);
    (*x)[ip]   = xyz[0];
    (*y)[ip]   = xyz[1];
    (*z)[ip]   = xyz[2];
    (*dx)[ip]  = tt.GetX() - x0;
    (*dy)[ip]  = tt.GetY() - tr->LocalY();
    (*dz)[ip]  = tt.GetZ() - tr->Z();
    (*pt)[ip]  = tr->Pt();
    (*dpt)[ip] = tt.Pt()- tr->Pt();
//    printf("pt : kalman[%e] MC[%e]\n", tt.Pt(), tr->Pt());
    (*budget)[ip] = tt.GetBudget(0);
    cc = tt.GetCovariance();
//    printf("dx[%5.2f] sy[%f]\n", (*dx)[ip], TMath::Sqrt(cc[0]));
    for(Int_t ic(0), jp(ip*15); ic<15; ic++, jp++) (*c)[jp]=cc[ic];
    ip++;
  }
  return kTRUE;
}
Commit	Line	Data
1ee39b3a	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id: AliTRDtrackInfo.cxx 27496 2008-07-22 08:35:45Z cblume $ */
	17
	18	////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Reconstruction QA //
	21	// //
	22	// Authors: //
	23	// Alex Bercuci <A.Bercuci@gsi.de> //
	24	// Markus Fasel <M.Fasel@gsi.de> //
	25	// //
	26	////////////////////////////////////////////////////////////////////////////
	27
	28	#include "TDatabasePDG.h"
b9ddd472	29	#include "TPDGCode.h"
4226db3e	30	#include "TVectorT.h"
1ee39b3a	31
1ee39b3a	32	#include "AliTrackReference.h"
eb05d549	33	#include "AliTrackPointArray.h"
1ee39b3a	34	#include "AliExternalTrackParam.h"
	35
	36	#include "AliTRDseedV1.h"
	37	#include "AliTRDtrackV1.h"
	38	#include "AliTRDgeometry.h"
	39	#include "AliTRDtrackerV1.h"
	40
	41	#include "AliTRDtrackInfo.h"
	42
	43	ClassImp(AliTRDtrackInfo)
	44	ClassImp(AliTRDtrackInfo::AliMCinfo)
	45	ClassImp(AliTRDtrackInfo::AliESDinfo)
5066aa9a	46	Double_t AliTRDtrackInfo::AliMCinfo::fgKalmanStep = 2.;
3ceb45ae	47	Bool_t AliTRDtrackInfo::AliMCinfo::fgKalmanUpdate = kTRUE;
1ee39b3a	48
	49	//___________________________________________________
	50	AliTRDtrackInfo::AliTRDtrackInfo():
	51	TObject()
	52	,fNClusters(0)
f232df0d	53	,fTRDtrack(NULL)
f232df0d	54	,fMC(NULL)
1ee39b3a	55	,fESD()
	56	{
	57	//
	58	// Default constructor
	59	//
	60	}
	61
	62
	63	//___________________________________________________
	64	AliTRDtrackInfo::AliTRDtrackInfo(const AliTRDtrackInfo &trdInfo):
	65	TObject((const TObject&)trdInfo)
	66	,fNClusters(trdInfo.fNClusters)
f232df0d	67	,fTRDtrack(NULL)
a340fe39	68	,fMC(NULL)
1ee39b3a	69	,fESD(trdInfo.fESD)
	70	{
	71	//
	72	// copy Entries
	73	//
	74
	75	if(trdInfo.fMC) fMC = new AliMCinfo(*trdInfo.fMC);
eb05d549	76	SetTrack(trdInfo.fTRDtrack);
1ee39b3a	77	}
	78
	79	//___________________________________________________
	80	AliTRDtrackInfo::AliMCinfo::AliMCinfo()
	81	:fLabel(0)
eb05d549	82	,fTRDlabel(0)
1ee39b3a	83	,fPDG(0)
1ee39b3a	84	,fNTrackRefs(0)
95d47440	85	,fEta(-999.)
	86	,fPhi(-999.)
	87	,fPt(-1.)
1ee39b3a	88	{
	89	// Set 0-Pointers
	90	memset(fTrackRefs, 0, sizeof(AliTrackReference ) 12);
	91	}
	92
	93	//___________________________________________________
	94	AliTRDtrackInfo::AliMCinfo::AliMCinfo(const AliMCinfo &mc)
	95	:fLabel(mc.fLabel)
eb05d549	96	,fTRDlabel(mc.fTRDlabel)
1ee39b3a	97	,fPDG(mc.fPDG)
1ee39b3a	98	,fNTrackRefs(mc.fNTrackRefs)
95d47440	99	,fEta(mc.fEta)
	100	,fPhi(mc.fPhi)
	101	,fPt(mc.fPt)
1ee39b3a	102	{
	103	//
	104	// Constructor
	105	//
	106
	107	memset(fTrackRefs, 0, sizeof(AliTrackReference ) 12);
	108	for(Int_t ien = 0; ien < 12; ien++){
	109	if(mc.fTrackRefs[ien])
	110	fTrackRefs[ien] = new AliTrackReference(*(mc.fTrackRefs[ien]));
	111	}
	112	}
	113
b9ddd472	114	//________________________________________________________
	115	Int_t AliTRDtrackInfo::AliMCinfo::GetPID() const
	116	{
61f6b45e	117	// Translate pdg code to PID index
61f6b45e	118
b9ddd472	119	switch(fPDG){
	120	case kElectron:
	121	case kPositron: return AliPID::kElectron;
	122	case kMuonPlus:
	123	case kMuonMinus: return AliPID::kMuon;
	124	case kPiPlus:
	125	case kPiMinus: return AliPID::kPion;
	126	case kKPlus:
	127	case kKMinus: return AliPID::kKaon;
	128	case kProton:
	129	case kProtonBar: return AliPID::kProton;
	130	}
	131	return -1;
	132	}
	133
1ee39b3a	134	//___________________________________________________
1ee39b3a	135	AliTRDtrackInfo::AliESDinfo::AliESDinfo()
5047978d	136	:fSteer(0)
b9ddd472	137	,fId(-1)
1ee39b3a	138	,fStatus(0)
	139	,fKinkIndex(0)
	140	,fTPCncls(0)
2b33fd9c	141	,fTPCdedx(0.)
f036f7ec	142	,fTOFbeta(0.)
3ceb45ae	143	,fTOFbc(0)
1ee39b3a	144	,fTRDpidQuality(0)
1ee39b3a	145	,fTRDnSlices(0)
eb05d549	146	,fPt(0.)
	147	,fPhi(-999.)
	148	,fEta(-999.)
f232df0d	149	,fTRDslices(NULL)
f232df0d	150	,fOP(NULL)
35983729	151	,fTPCout(NULL)
b9058d72	152	,fITSout(NULL)
eb05d549	153	,fTPArray(NULL)
1ee39b3a	154	{
	155	//
	156	// Constructor
	157	//
	158
	159	memset(fTRDr, 0, AliPID::kSPECIES*sizeof(Double32_t));
d80a6a00	160	memset(fTRDv0pid, 0, AliPID::kSPECIES*sizeof(Int_t));
1ee39b3a	161	}
	162
	163	//___________________________________________________
	164	AliTRDtrackInfo::AliESDinfo::AliESDinfo(const AliESDinfo &esd)
5047978d	165	:fSteer(esd.fSteer)
b9ddd472	166	,fId(esd.fId)
1ee39b3a	167	,fStatus(esd.fStatus)
	168	,fKinkIndex(esd.fKinkIndex)
	169	,fTPCncls(esd.fTPCncls)
2b33fd9c	170	,fTPCdedx(esd.fTPCdedx)
f036f7ec	171	,fTOFbeta(esd.fTOFbeta)
3ceb45ae	172	,fTOFbc(esd.fTOFbc)
1ee39b3a	173	,fTRDpidQuality(esd.fTRDpidQuality)
1ee39b3a	174	,fTRDnSlices(esd.fTRDnSlices)
eb05d549	175	,fPt(esd.fPt)
	176	,fPhi(esd.fPhi)
	177	,fEta(esd.fEta)
f232df0d	178	,fTRDslices(NULL)
f232df0d	179	,fOP(NULL)
35983729	180	,fTPCout(NULL)
b9058d72	181	,fITSout(NULL)
eb05d549	182	,fTPArray(NULL)
1ee39b3a	183	{
	184	//
	185	// Constructor
	186	//
	187
	188	memcpy(fTRDr, esd.fTRDr, AliPID::kSPECIES*sizeof(Double32_t));
d80a6a00	189	memcpy(fTRDv0pid, esd.fTRDv0pid, AliPID::kSPECIES*sizeof(Int_t));
1ee39b3a	190
	191	if(fTRDnSlices){
	192	fTRDslices = new Double32_t[fTRDnSlices];
	193	memcpy(fTRDslices, esd.fTRDslices, fTRDnSlices*sizeof(Double32_t));
	194	}
eb05d549	195	SetOuterParam(esd.fOP);
	196	SetTPCoutParam(esd.fTPCout);
	197	SetITSoutParam(esd.fITSout);
	198	SetTrackPointArray(esd.fTPArray);
1ee39b3a	199	}
	200
	201
	202	//___________________________________________________
	203	AliTRDtrackInfo::~AliTRDtrackInfo()
	204	{
	205	//
	206	// Destructor
	207	//
	208
	209	if(fMC) delete fMC;
	210	if(fTRDtrack) delete fTRDtrack;
	211	}
	212
	213	//___________________________________________________
	214	AliTRDtrackInfo::AliMCinfo::~AliMCinfo()
	215	{
	216	//
	217	// Destructor
	218	//
	219
	220	fNTrackRefs = 0;
	221	for(Int_t ien = 0; ien < 12; ien++){
	222	if(fTrackRefs[ien]) delete fTrackRefs[ien];
f232df0d	223	fTrackRefs[ien] = NULL;
1ee39b3a	224	}
	225	}
	226
	227	//___________________________________________________
	228	AliTRDtrackInfo::AliESDinfo::~AliESDinfo()
	229	{
	230	//
	231	// Destructor
	232	//
	233
	234	if(fTRDnSlices){
	235	delete [] fTRDslices;
f232df0d	236	fTRDslices = NULL;
1ee39b3a	237	fTRDnSlices = 0;
1ee39b3a	238	}
f232df0d	239	if(fOP) delete fOP; fOP = NULL;
35983729	240	if(fTPCout) delete fTPCout; fTPCout = NULL;
b9058d72	241	if(fITSout) delete fITSout; fITSout = NULL;
eb05d549	242	if(fTPArray) delete fTPArray;
1ee39b3a	243	}
1ee39b3a	244
8778aca3	245	//___________________________________________________
	246	void AliTRDtrackInfo::AliESDinfo::Delete(const Option_t *){
	247	//
	248	// Delete Pointer members
	249	//
	250	if(fTRDnSlices){
	251	delete [] fTRDslices;
	252	fTRDslices = NULL;
	253	fTRDnSlices = 0;
	254	}
	255	if(fOP) delete fOP; fOP = NULL;
35983729	256	if(fTPCout) delete fTPCout; fTPCout = NULL;
b9058d72	257	if(fITSout) delete fITSout; fITSout = NULL;
eb05d549	258	if(fTPArray) delete fTPArray; fTPArray = NULL;
8778aca3	259	}
8778aca3	260
1ee39b3a	261
	262	//___________________________________________________
	263	AliTRDtrackInfo& AliTRDtrackInfo::operator=(const AliTRDtrackInfo &trdInfo)
	264	{
	265	//
	266	// = Operator
	267	//
01ccc21a	268	if(this == &trdInfo) return *this;
1ee39b3a	269
1ee39b3a	270	fNClusters = trdInfo.fNClusters;
1881bb2b	271	fESD = trdInfo.fESD;
1ee39b3a	272
1ee39b3a	273	if(trdInfo.fMC){
2c33fb46	274	if(!fMC) fMC = new AliMCinfo(*trdInfo.fMC);
	275	else{
	276	fMC->~AliMCinfo();
1ee39b3a	277	new(fMC) AliMCinfo(*trdInfo.fMC);
2c33fb46	278	}
1881bb2b	279	} else {if(fMC) delete fMC; fMC = NULL;}
1ee39b3a	280
eb05d549	281	SetTrack(trdInfo.fTRDtrack);
1ee39b3a	282
	283	return *this;
	284	}
	285
	286	//___________________________________________________
	287	AliTRDtrackInfo::AliMCinfo& AliTRDtrackInfo::AliMCinfo::operator=(const AliMCinfo &mc)
	288	{
	289	//
	290	// Assignment operator
	291	//
	292
01ccc21a	293	if(this == &mc) return *this;
1ee39b3a	294	fLabel = mc.fLabel;
eb05d549	295	fTRDlabel = mc.fTRDlabel;
1ee39b3a	296	fPDG = mc.fPDG;
	297	fNTrackRefs = mc.fNTrackRefs;
	298
	299	AliTrackReference **itr = &fTrackRefs[0];
	300	AliTrackReference* const *jtr = &mc.fTrackRefs[0];
eb05d549	301	for(Int_t ien = 0; ien < 2*AliTRDgeometry::kNlayer; ien++, itr++, jtr++){
1ee39b3a	302	if((*jtr)){
1ee39b3a	303	if(!(itr)) (itr) = new AliTrackReference((jtr));
2c33fb46	304	else{
	305	(*itr)->~AliTrackReference();
	306	new(&(itr)) AliTrackReference((*jtr));
	307	}
1881bb2b	308	} else {if((itr)) delete (itr); (*itr) = NULL;}
1ee39b3a	309	}
	310	return *this;
	311	}
	312
	313	//___________________________________________________
	314	AliTRDtrackInfo::AliESDinfo& AliTRDtrackInfo::AliESDinfo::operator=(const AliESDinfo &esd)
	315	{
	316	//
	317	// Assignment operator
	318	//
	319
01ccc21a	320	if(this == &esd) return *this;
5047978d	321	fSteer = esd.fSteer;
1ee39b3a	322	fId = esd.fId;
01ccc21a	323	fStatus = esd.fStatus;
	324	fKinkIndex = esd.fKinkIndex;
	325	fTPCncls = esd.fTPCncls;
2b33fd9c	326	fTPCdedx = esd.fTPCdedx;
f036f7ec	327	fTOFbeta = esd.fTOFbeta;
01ccc21a	328	fTOFbc = esd.fTOFbc;
01ccc21a	329	fTRDpidQuality= esd.fTRDpidQuality;
1ee39b3a	330	fTRDnSlices = esd.fTRDnSlices;
eb05d549	331	fPt = esd.fPt;
	332	fPhi = esd.fPhi;
	333	fEta = esd.fEta;
178d284a	334
1ee39b3a	335	memcpy(fTRDr, esd.fTRDr, AliPID::kSPECIES*sizeof(Double32_t));
01ccc21a	336	memcpy(fTRDv0pid, esd.fTRDv0pid, AliPID::kSPECIES*sizeof(Int_t));
1ee39b3a	337
1ee39b3a	338	if(fTRDnSlices){
56838f9d	339	if(!fTRDslices) fTRDslices = new Double32_t[fTRDnSlices];
1ee39b3a	340	memcpy(fTRDslices, esd.fTRDslices, fTRDnSlices*sizeof(Double32_t));
1ee39b3a	341	}
eb05d549	342	SetOuterParam(esd.fOP);
	343	SetTPCoutParam(esd.fTPCout);
	344	SetITSoutParam(esd.fITSout);
	345	SetTrackPointArray(esd.fTPArray);
1ee39b3a	346
	347	return *this;
	348	}
	349
	350	//___________________________________________________
	351	void AliTRDtrackInfo::Delete(const Option_t *)
	352	{
	353	//
	354	// Delete
	355	//
	356
2c33fb46	357	AliDebug(2, Form("track[%p] mc[%p]", (void)fTRDtrack, (void)fMC));
1ee39b3a	358	fNClusters = 0;
f232df0d	359	if(fMC) delete fMC; fMC = NULL;
8778aca3	360	fESD.Delete(NULL);
f232df0d	361	if(fTRDtrack) delete fTRDtrack; fTRDtrack = NULL;
1ee39b3a	362	}
	363
	364	//___________________________________________________
	365	void AliTRDtrackInfo::SetTrack(const AliTRDtrackV1 *track)
	366	{
	367	//
	368	// Set the TRD track
	369	//
	370
eb05d549	371	if(track){
	372	if(!fTRDtrack) fTRDtrack = new AliTRDtrackV1(*track);
	373	else{
	374	fTRDtrack->~AliTRDtrackV1();
	375	new(fTRDtrack) AliTRDtrackV1(*track);
	376	}
	377	if(track->IsOwner()) fTRDtrack->SetOwner();
	378	} else {
	379	if(fTRDtrack) delete fTRDtrack; fTRDtrack = NULL;
	380	}
	381	}
	382
	383	//___________________________________________________
	384	void AliTRDtrackInfo::AliESDinfo::SetTrackPointArray(const AliTrackPointArray *tps)
	385	{
	386	//
	387	// Set the track point array for alignment task
	388	//
	389
	390
	391	if(tps){
	392	if(!fTPArray) fTPArray = new AliTrackPointArray(*tps);
	393	else{
	394	fTPArray->~AliTrackPointArray();
	395	new(fTPArray) AliTrackPointArray(*tps);
	396	}
	397	} else {
	398	if(fTPArray) delete fTPArray; fTPArray = NULL;
2c33fb46	399	}
1ee39b3a	400	}
	401
	402	//___________________________________________________
	403	void AliTRDtrackInfo::AddTrackRef(const AliTrackReference *tref)
	404	{
	405	//
	406	// Add track reference
	407	//
	408
eb05d549	409	if(fMC->fNTrackRefs >= 2*AliTRDgeometry::kNlayer){
1ee39b3a	410	SetCurved();
	411	return;
	412	}
	413	// Make a copy for the object in order to avoid ownership problems
	414	fMC->fTrackRefs[fMC->fNTrackRefs++] = new AliTrackReference(*tref);
	415	}
	416
	417	//___________________________________________________
	418	AliTrackReference* AliTRDtrackInfo::GetTrackRef(Int_t idx) const
	419	{
	420	//
	421	// Returns a track reference
	422	//
f232df0d	423	if(!fMC) return NULL;
f232df0d	424	return (idx>=0 && idx < 12) ? fMC->fTrackRefs[idx] : NULL;
1ee39b3a	425	}
	426
	427	//___________________________________________________
61f6b45e	428	AliTrackReference* AliTRDtrackInfo::GetTrackRef(const AliTRDseedV1* const tracklet) const
1ee39b3a	429	{
	430	//
	431	// Returns a track reference
	432	//
f232df0d	433	if(!fMC) return NULL;
1ee39b3a	434	Double_t cw = AliTRDgeometry::CamHght() + AliTRDgeometry::CdrHght();
	435	AliTrackReference * const* jtr = &(fMC->fTrackRefs[0]);
	436	for(Int_t itr = 0; itr < fMC->fNTrackRefs; itr++, ++jtr){
	437	if(!(*jtr)) break;
	438	if(TMath::Abs(tracklet->GetX0() - (jtr)->LocalX()) < cw) return (jtr);
	439	}
f232df0d	440	return NULL;
1ee39b3a	441	}
	442
	443	//___________________________________________________
	444	Int_t AliTRDtrackInfo::GetNumberOfClusters() const
	445	{
	446	//
	447	// Returns the number of clusters
	448	//
	449
	450	Int_t n = 0;
	451	if(!fTRDtrack) return 0;
	452	if(fTRDtrack->GetNumberOfTracklets() == 0) return n;
f232df0d	453	AliTRDseedV1 *tracklet = NULL;
eb05d549	454	for(Int_t ip=0; ip<AliTRDgeometry::kNlayer; ip++){
1ee39b3a	455	if(!(tracklet = const_cast<AliTRDseedV1 *>(fTRDtrack->GetTracklet(ip)))) continue;
	456	n+=tracklet->GetN();
	457	}
	458	return n;
	459	}
	460
	461
	462	//___________________________________________________
eb05d549	463	void AliTRDtrackInfo::AliESDinfo::SetOuterParam(const AliExternalTrackParam *op)
1ee39b3a	464	{
	465	//
	466	// Set outer track parameters
	467	//
	468
eb05d549	469	if(op){
	470	if(fOP){
	471	fOP->~AliExternalTrackParam();
	472	// RS: Constructor from VTrack was used instead of Constructor from AliExternalTrackParam
	473	new(fOP) AliExternalTrackParam(*op);
	474	} else fOP = new AliExternalTrackParam(*op);
	475	} else {
	476	if(fOP) delete fOP; fOP = NULL;
	477	}
1ee39b3a	478	}
1ee39b3a	479
b9058d72	480	//___________________________________________________
eb05d549	481	void AliTRDtrackInfo::AliESDinfo::SetITSoutParam(const AliExternalTrackParam *op)
b9058d72	482	{
	483	//
	484	// Set TPCout track parameters
	485	//
	486
eb05d549	487	if(op){
	488	if(fITSout){
	489	fITSout->~AliExternalTrackParam();
	490	// RS: Constructor from VTrack was used instead of Constructor from AliExternalTrackParam
	491	new(fITSout) AliExternalTrackParam(*op);
	492	} else fITSout = new AliExternalTrackParam(*op);
	493	} else {
	494	if(fITSout) delete fITSout; fITSout = NULL;
	495	}
b9058d72	496	}
b9058d72	497
35983729	498	//___________________________________________________
eb05d549	499	void AliTRDtrackInfo::AliESDinfo::SetTPCoutParam(const AliExternalTrackParam *op)
35983729	500	{
	501	//
	502	// Set TPCout track parameters
	503	//
	504
eb05d549	505	if(op){
	506	if(fTPCout){
	507	fTPCout->~AliExternalTrackParam();
	508	// RS: Constructor from VTrack was used instead of Constructor from AliExternalTrackParam
	509	new(fTPCout) AliExternalTrackParam(*op);
	510	} else fTPCout = new AliExternalTrackParam(*op);
	511	} else {
	512	if(fTPCout) delete fTPCout; fTPCout = NULL;
	513	}
35983729	514	}
35983729	515
1ee39b3a	516	//___________________________________________________
	517	Int_t AliTRDtrackInfo::GetNTracklets() const
	518	{
	519	//
	520	// Return the number of tracklets
	521	//
	522
e2927481	523	if(!fTRDtrack) return 0;
1ee39b3a	524	return fTRDtrack->GetNumberOfTracklets();
	525	}
	526
	527	//___________________________________________________
	528	void AliTRDtrackInfo::SetSlices(Int_t n, Double32_t *s)
	529	{
	530	//
	531	// Set the slices
	532	//
1ee39b3a	533	if(fESD.fTRDnSlices != n){
	534	fESD.fTRDnSlices = 0;
	535	delete [] fESD.fTRDslices;
f232df0d	536	fESD.fTRDslices = NULL;
1ee39b3a	537	}
	538
	539	if(!fESD.fTRDnSlices){
	540	fESD.fTRDnSlices = n;
f232df0d	541	fESD.fTRDslices = new Double32_t[fESD.fTRDnSlices];
1ee39b3a	542	}
	543
	544	memcpy(fESD.fTRDslices, s, n*sizeof(Double32_t));
	545	}
	546
	547	//___________________________________________________
0f2c4c4f	548	Bool_t AliTRDtrackInfo::AliMCinfo::GetDirections(Float_t &x0, Float_t &y0, Float_t &z0, Float_t &dydx, Float_t &dzdx, Float_t &pt, Float_t &p, Float_t &eta, Float_t &phi, UChar_t &status) const
1ee39b3a	549	{
	550	// Check for 2 track ref for the tracklet defined bythe radial position x0
	551	// The "status" is a bit map and gives a more informative output in case of failure:
	552	// - 0 : everything is OK
	553	// - BIT(0) : 0 track refs found
	554	// - BIT(1) : 1 track reference found
	555	// - BIT(2) : dx <= 0 between track references
	556	// - BIT(3) : dx > 0 && dx < 3.7 - tangent tracks
	557
	558	status = 0;
	559	Double_t cw = AliTRDgeometry::CamHght() + AliTRDgeometry::CdrHght();
	560	Int_t nFound = 0;
f232df0d	561	AliTrackReference *tr[2] = {NULL, NULL};
1ee39b3a	562	AliTrackReference * const* jtr = &fTrackRefs[0];
	563	for(Int_t itr = 0; itr < fNTrackRefs; itr++, ++jtr){
	564	if(!(*jtr)) break;
	565	/*
	566	if(fDebugLevel>=5) printf("\t\tref[%2d] x[%6.3f]\n", itr, (jtr)->LocalX());/
	567	if(TMath::Abs(x0 - (*jtr)->LocalX()) > cw) continue;
	568	tr[nFound++] = (*jtr);
	569	if(nFound == 2) break;
	570	}
	571	if(nFound < 2){
e02d9eb3	572	//AliDebug(1, Form("Missing track ref x0[%6.3f] nref[%d]", x0, nFound));
1ee39b3a	573	if(!nFound) SETBIT(status, 0);
	574	else SETBIT(status, 1);
	575	return kFALSE;
	576	}
61f6b45e	577	pt=tr[1]->Pt();
0f2c4c4f	578	p =tr[1]->P();
61f6b45e	579	if(pt < 1.e-3) return kFALSE;
1ee39b3a	580
	581	Double_t dx = tr[1]->LocalX() - tr[0]->LocalX();
	582	if(dx <= 0.){
	583	AliWarningGeneral("AliTRDtrackInfo::AliMCinfo::GetDirections()", Form("Track ref with wrong radial distances refX0[%6.3f] refX1[%6.3f]", tr[0]->LocalX(), tr[1]->LocalX()));
	584	SETBIT(status, 2);
	585	return kFALSE;
	586	}
	587	if(TMath::Abs(dx-AliTRDgeometry::CamHght()-AliTRDgeometry::CdrHght())>1.E-3) SETBIT(status, 3);
	588
	589	dydx = (tr[1]->LocalY() - tr[0]->LocalY()) / dx;
	590	dzdx = (tr[1]->Z() - tr[0]->Z()) / dx;
	591	//Float_t dx0 = tr[1]->LocalX() - x0;
	592	y0 = tr[1]->LocalY()/* - dydxdx0/;
	593	z0 = tr[1]->Z()/* - dzdxdx0/;
	594	x0 = tr[1]->LocalX();
3ceb45ae	595	eta = -TMath::Log(TMath::Tan(0.5 * tr[1]->Theta()));
01ccc21a	596	phi = TMath::ATan2(tr[1]->Y(), tr[1]->X());
1ee39b3a	597	return kTRUE;
	598	}
	599
	600	//___________________________________________________
3ceb45ae	601	Bool_t AliTRDtrackInfo::AliMCinfo::PropagateKalman(
	602	TVectorD x, TVectorD y, TVectorD *z,
	603	TVectorD dx, TVectorD dy, TVectorD *dz,
	604	TVectorD pt, TVectorD dpt, TVectorD budget, TVectorD c, Double_t mass) const
1ee39b3a	605	{
	606	// Propagate Kalman from the first TRD track reference to
	607	// last one and save residuals in the y, z and pt.
	608	//
	609	// This is to calibrate the dEdx and MS corrections
	610
3ceb45ae	611	if(!fNTrackRefs) return kFALSE;
1ee39b3a	612	for(Int_t itr=kNTrackRefs; itr--;){
3ceb45ae	613	(x)[itr] = 0.;(y)[itr] = 0.;(*z)[itr] = 0.;
4226db3e	614	(dx)[itr] = -1.; (dy)[itr] = 100.; (dz)[itr] = 100.; (dpt)[itr] = 100.;
1ee39b3a	615	}
1ee39b3a	616
1ee39b3a	617	// Initialize TRD track to the first track reference
4226db3e	618	AliTrackReference *tr(NULL);
4226db3e	619	Int_t itr(0); while(!(tr = fTrackRefs[itr])) itr++;
3ceb45ae	620	if(tr->Pt()<1.e-3) return kFALSE;
1ee39b3a	621
	622	AliTRDtrackV1 tt;
	623	Double_t xyz[3]={tr->X(),tr->Y(),tr->Z()};
	624	Double_t pxyz[3]={tr->Px(),tr->Py(),tr->Pz()};
	625	Double_t var[6] = {1.e-4, 1.e-4, 1.e-4, 1.e-4, 1.e-4, 1.e-4};
	626	Double_t cov[21]={
	627	var[0], 0., 0., 0., 0., 0.,
	628	var[1], 0., 0., 0., 0.,
	629	var[2], 0., 0., 0.,
	630	var[3], 0., 0.,
	631	var[4], 0.,
	632	var[5]
	633	};
	634	TDatabasePDG db;
	635	const TParticlePDG *pdg=db.GetParticle(fPDG);
	636	if(!pdg){
	637	AliWarningGeneral("AliTRDtrackInfo::AliMCinfo::PropagateKalman()", Form("PDG entry missing for code %d. References for track %d", fPDG, fNTrackRefs));
3ceb45ae	638	return kFALSE;
1ee39b3a	639	}
1ee39b3a	640	tt.Set(xyz, pxyz, cov, Short_t(pdg->Charge()));
3ceb45ae	641	if(mass<0){ // mass 0 use PDG
	642	tt.SetMass(pdg->Mass());
	643	} else { // use rec value
	644	tt.SetMass(mass);
	645	}
	646
	647	// Double_t bg(tr->P()/pdg->Mass());
	648	// printf("\n\nNEW track PDG[%d] bg[%f] x[%f]\n", fPDG, bg, TMath::Log(bg));
4226db3e	649	Double_t x0(tr->LocalX());
	650	const Double_t *cc(NULL);
	651	for(Int_t ip=0; itr<fNTrackRefs; itr++){
	652	if(!(tr = fTrackRefs[itr])) continue;
3ceb45ae	653	// printf("ip[%d] det[%d]\n", ip, tr->DetectorId());
5066aa9a	654	if(!AliTRDtrackerV1::PropagateToX(tt, tr->LocalX(), fgKalmanStep)) continue;
1ee39b3a	655
1ee39b3a	656	//if(update) ...
3ceb45ae	657	tt.GetXYZ(xyz);
	658	(*x)[ip] = xyz[0];
	659	(*y)[ip] = xyz[1];
	660	(*z)[ip] = xyz[2];
4226db3e	661	(*dx)[ip] = tt.GetX() - x0;
	662	(*dy)[ip] = tt.GetY() - tr->LocalY();
	663	(*dz)[ip] = tt.GetZ() - tr->Z();
5066aa9a	664	(*pt)[ip] = tr->Pt();
4226db3e	665	(*dpt)[ip] = tt.Pt()- tr->Pt();
3ceb45ae	666	// printf("pt : kalman[%e] MC[%e]\n", tt.Pt(), tr->Pt());
3ceb45ae	667	(*budget)[ip] = tt.GetBudget(0);
1ee39b3a	668	cc = tt.GetCovariance();
3ceb45ae	669	// printf("dx[%5.2f] sy[%f]\n", (*dx)[ip], TMath::Sqrt(cc[0]));
3ceb45ae	670	for(Int_t ic(0), jp(ip15); ic<15; ic++, jp++) (c)[jp]=cc[ic];
1ee39b3a	671	ip++;
1ee39b3a	672	}
3ceb45ae	673	return kTRUE;
1ee39b3a	674	}