#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliTracker.h"
#include "AliTRDpidESD.h"
#include "AliTRDgeometry.h"
#include "AliTRDcalibDB.h"
#include "Cal/AliTRDCalPID.h"
ClassImp(AliTRDpidESD)
Bool_t AliTRDpidESD::fgCheckTrackStatus = kTRUE;
Bool_t AliTRDpidESD::fgCheckKinkStatus = kFALSE;
Int_t AliTRDpidESD::fgMinPlane = 0;
AliTRDpidESD::AliTRDpidESD()
:TObject()
,fTrack(NULL)
{
}
AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p)
:TObject(p)
,fTrack(NULL)
{
((AliTRDpidESD &) p).Copy(*this);
}
AliTRDpidESD::~AliTRDpidESD()
{
if(fTrack) delete fTrack;
}
AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p)
{
if (this != &p) ((AliTRDpidESD &) p).Copy(*this);
return *this;
}
void AliTRDpidESD::Copy(TObject &p) const
{
((AliTRDpidESD &) p).fgCheckTrackStatus = fgCheckTrackStatus;
((AliTRDpidESD &) p).fgCheckKinkStatus = fgCheckKinkStatus;
((AliTRDpidESD &) p).fgMinPlane = fgMinPlane;
((AliTRDpidESD &) p).fTrack = NULL;
}
Int_t AliTRDpidESD::MakePID(AliESDEvent * const event)
{
AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
if (!calibration) {
AliErrorGeneral("AliTRDpidESD::MakePID()"
,"No access to calibration data");
return -1;
}
const AliTRDCalPID *pd = calibration->GetPIDObject(AliTRDpidUtil::kLQ);
if (!pd) {
AliErrorGeneral("AliTRDpidESD::MakePID()"
,"No access to AliTRDCalPID");
return -1;
}
Double_t p[10];
AliESDtrack *t = NULL;
Float_t dedx[AliTRDCalPID::kNSlicesLQ], dEdx;
Int_t timebin;
Float_t mom, length;
Int_t nPlanePID;
for (Int_t i=0; i<event->GetNumberOfTracks(); i++) {
t = event->GetTrack(i);
if(!CheckTrack(t)) continue;
if (t->GetTRDsignal()<1.e-5) continue;
mom = 0.;
length = 0.;
nPlanePID = 0;
for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
p[iSpecies] = 1./AliPID::kSPECIES;
for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNlayer; iLayer++) {
for(int iSlice=0; iSlice<AliTRDCalPID::kNSlicesLQ; iSlice++)
dedx[iSlice] = t->GetTRDslice(iLayer, iSlice);
dEdx = t->GetTRDslice(iLayer, -1);
timebin = t->GetTRDTimBin(iLayer);
if ((dEdx <= 0.) || (timebin <= -1.)) continue;
if(!RecalculateTrackSegmentKine(t, iLayer, mom, length)){
}
nPlanePID++;
for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length, iLayer);
}
}
if(nPlanePID == 0) continue;
Double_t probTotal = 0.;
for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += p[iSpecies];
if(probTotal <= 0.){
AliWarning(Form("The total probability (%e) over all species <= 0 in ESD track %d."
, probTotal, i));
AliWarning("This may be caused by some error in reference data.");
AliWarning("Calculation continues but results might be corrupted.");
continue;
}
for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal;
t->SetTRDpid(p);
t->SetTRDntracklets(nPlanePID<<3);
}
return 0;
}
Bool_t AliTRDpidESD::CheckTrack(AliESDtrack * const t)
{
if (fgCheckTrackStatus) {
if (((t->GetStatus() & AliESDtrack::kTRDout ) == 0) &&
((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) return kFALSE;
}
if (fgCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE;
return kTRUE;
}
Bool_t AliTRDpidESD::RecalculateTrackSegmentKine(AliESDtrack * const esd
, Int_t layer
, Float_t &mom
, Float_t &length)
{
const Float_t kAmHalfWidth = AliTRDgeometry::AmThick() / 2.;
const Float_t kDrWidth = AliTRDgeometry::DrThick();
const Float_t kTime0 = AliTRDgeometry::GetTime0(layer);
length = 2 * kAmHalfWidth + kDrWidth;
const AliExternalTrackParam *op = esd->GetOuterParam();
if(!op){
mom = esd->GetP();
return kFALSE;
}
AliExternalTrackParam *param = NULL;
if(!fTrack){
fTrack = new AliExternalTrackParam(*op);
param = fTrack;
} else param = new(fTrack) AliExternalTrackParam(*op);
Double_t xyz0[3];
op->GetXYZ(xyz0);
Float_t field = AliTracker::GetBz(xyz0);
Double_t s, t;
if(!param->PropagateTo(kTime0-kAmHalfWidth-kDrWidth, field)){
mom = op->GetP();
s = op->GetSnp();
t = op->GetTgl();
if (s < 1.) length /= TMath::Sqrt((1.-s)*(1.+s) / (1. + t*t));
return kFALSE;
}
mom = param->GetP();
s = param->GetSnp();
t = param->GetTgl();
if (s < 1.) length /= TMath::Sqrt((1.-s)*(1.+s) / (1. + t*t));
Double_t alpha = param->GetAlpha();
if(!param->PropagateTo(kTime0+kAmHalfWidth, field)) return kFALSE;
if(TMath::Abs(alpha-param->GetAlpha())>.01) return kFALSE;
return kTRUE;
}