#include <TSystem.h>
#include "AliLog.h"
+
+#include "AliVTrack.h"
#include "AliTRDPIDResponseObject.h"
#include "AliTRDPIDResponse.h"
-#include "AliTRDTKDInterpolator.h"
+//#include "AliTRDTKDInterpolator.h"
+#include "AliTRDNDFast.h"
+#include "AliTRDdEdxParams.h"
ClassImp(AliTRDPIDResponse)
AliTRDPIDResponse::AliTRDPIDResponse():
TObject()
,fkPIDResponseObject(NULL)
+ ,fkTRDdEdxParams(NULL)
,fGainNormalisationFactor(1.)
{
//
AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
TObject(ref)
,fkPIDResponseObject(NULL)
+ ,fkTRDdEdxParams(NULL)
,fGainNormalisationFactor(ref.fGainNormalisationFactor)
{
//
TObject::operator=(ref);
fGainNormalisationFactor = ref.fGainNormalisationFactor;
fkPIDResponseObject = ref.fkPIDResponseObject;
-
+ fkTRDdEdxParams = ref.fkTRDdEdxParams;
+
return *this;
}
//
// Destructor
//
- if(IsOwner()) delete fkPIDResponseObject;
+ if(IsOwner()) {
+ delete fkPIDResponseObject;
+ delete fkTRDdEdxParams;
+ }
}
//____________________________________________________________
return kTRUE;
}
+
+
+//____________________________________________________________
+Double_t AliTRDPIDResponse::GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType type) const
+{
+ //
+ //calculate the TRD nSigma
+ //
+
+ const Double_t badval = -9999;
+ Double_t info[5]; for(int i=0; i<5; i++){info[i]=badval;}
+
+ const Double_t delta = GetSignalDelta(track, type, kFALSE, info);
+
+ const Double_t mean = info[0];
+ const Double_t res = info[1];
+ if(res<0){
+ return badval;
+ }
+
+ const Double_t sigma = mean*res;
+ const Double_t eps = 1e-12;
+ return delta/(sigma + eps);
+}
+
+//____________________________________________________________
+Double_t AliTRDPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/, Double_t *info/*=0x0*/) const
+{
+ //
+ //calculate the TRD signal difference w.r.t. the expected
+ //output other information in info[]
+ //
+
+ const Double_t badval = -9999;
+
+ if(!track){
+ return badval;
+ }
+
+ Double_t pTRD = -999;
+ for(Int_t ich=0; ich<6; ich++){
+ pTRD = track->GetTRDmomentum(ich);
+ if(pTRD>0)
+ break;
+ }
+ if(pTRD<0){
+ return badval;
+ }
+
+ if(!fkTRDdEdxParams){
+ AliError("fkTRDdEdxParams null");
+ return -99999;
+ }
+
+ const Double_t nch = track->GetTRDNchamberdEdx();
+ const Double_t ncls = track->GetTRDNclusterdEdx();
+
+ const TVectorF meanvec = fkTRDdEdxParams->GetMeanParameter(type, nch, ncls);
+ if(meanvec.GetNrows()==0){
+ return badval;
+ }
+
+ const TVectorF resvec = fkTRDdEdxParams->GetSigmaParameter(type, nch, ncls);
+ if(resvec.GetNrows()==0){
+ return badval;
+ }
+
+ const Float_t *meanpar = meanvec.GetMatrixArray();
+ const Float_t *respar = resvec.GetMatrixArray();
+
+ //============================================================================================<<<<<<<<<<<<<
+
+ const Double_t bg = pTRD/AliPID::ParticleMass(type);
+ const Double_t expsig = MeandEdxTR(&bg, meanpar);
+
+ if(info){
+ info[0]= expsig;
+ info[1]= ResolutiondEdxTR(&ncls, respar);
+ }
+
+ const Double_t eps = 1e-10;
+
+ if(ratio){
+ return track->GetTRDsignal()/(expsig + eps);
+ }
+ else{
+ return track->GetTRDsignal() - expsig;
+ }
+}
+
+
+Double_t AliTRDPIDResponse::ResolutiondEdxTR(const Double_t * xx, const Float_t * par)
+{
+ //
+ //resolution
+ //npar=3
+ //
+
+ const Double_t ncls = xx[0];
+
+ return par[0]+par[1]*TMath::Power(ncls, par[2]);
+}
+
+Double_t AliTRDPIDResponse::MeandEdxTR(const Double_t * xx, const Float_t * pin)
+{
+ //
+ //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
+ //npar = 8 = 3+5
+ //
+
+ Float_t ptr[4]={0};
+ for(int ii=0; ii<3; ii++){
+ ptr[ii+1]=pin[ii];
+ }
+ return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
+}
+
+Double_t AliTRDPIDResponse::MeanTR(const Double_t * xx, const Float_t * par)
+{
+ //
+ //LOGISTIC parametrization for TR, in unit of MIP
+ //npar = 4
+ //
+
+ const Double_t bg = xx[0];
+ const Double_t gamma = sqrt(1+bg*bg);
+
+ const Double_t p0 = TMath::Abs(par[1]);
+ const Double_t p1 = TMath::Abs(par[2]);
+ const Double_t p2 = TMath::Abs(par[3]);
+
+ const Double_t zz = TMath::Log(gamma);
+ const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
+
+ return par[0]+tryield;
+}
+
+Double_t AliTRDPIDResponse::MeandEdx(const Double_t * xx, const Float_t * par)
+{
+ //
+ //ALEPH parametrization for dEdx
+ //npar = 5
+ //
+
+ const Double_t bg = xx[0];
+ const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
+
+ const Double_t p0 = TMath::Abs(par[0]);
+ const Double_t p1 = TMath::Abs(par[1]);
+
+ const Double_t p2 = TMath::Abs(par[2]);
+
+ const Double_t p3 = TMath::Abs(par[3]);
+ const Double_t p4 = TMath::Abs(par[4]);
+
+ const Double_t aa = TMath::Power(beta, p3);
+
+ const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
+
+ return (p1-aa-bb)*p0/aa;
+
+}
+
+
//____________________________________________________________
-Bool_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES],ETRDPIDMethod PIDmethod,Bool_t kNorm) const
+Int_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES],ETRDPIDMethod PIDmethod,Bool_t kNorm) const
{
//
// Calculate TRD likelihood values for the track based on dedx and
// kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
//
// Return value
- // true if calculation success
+ // number of tracklets used for PID, 0 if no PID
//
AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
if(!fkPIDResponseObject){
- AliWarning("Missing reference data. PID calculation not possible.");
- return kFALSE;
+ AliDebug(3,"Missing reference data. PID calculation not possible.");
+ return 0;
}
for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
Double_t prLayer[AliPID::kSPECIES];
Double_t dE[10], s(0.);
+ Int_t ntrackletsPID=0;
for(Int_t il(kNlayer); il--;){
memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;
-
- s=0.;
+ s=0.;
Bool_t filled=kTRUE;
for(Int_t is(AliPID::kSPECIES); is--;){
- if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
+ //if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod);
AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
if(prLayer[is]<1.e-30){
if(kNorm) prLayer[is] /= s;
prob[is] *= prLayer[is];
}
+ ntrackletsPID++;
}
- if(!kNorm) return kTRUE;
+ if(!kNorm) return ntrackletsPID;
s=0.;
for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
if(s<1.e-30){
AliDebug(2, "Null total prob.");
- return kFALSE;
+ return 0;
}
for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
- return kTRUE;
+ return ntrackletsPID;
}
//____________________________________________________________
AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
Double_t probLayer = 0.;
+
Float_t pLower, pUpper;
- switch(PIDmethod){
- case kNN: // NN
+ AliTRDNDFast *refUpper = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,PIDmethod)),
+ *refLower = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,PIDmethod));
+
+ // Do Interpolation exept for underflow and overflow
+ if(refLower && refUpper){
+ Double_t probLower = refLower->Eval(dEdx);
+ Double_t probUpper = refUpper->Eval(dEdx);
+
+ probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
+ } else if(refLower){
+ // underflow
+ probLayer = refLower->Eval(dEdx);
+ } else if(refUpper){
+ // overflow
+ probLayer = refUpper->Eval(dEdx);
+ } else {
+ AliError("No references available");
+ }
+ AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
+
+ return probLayer;
+
+/* old implementation
+
+switch(PIDmethod){
+case kNN: // NN
break;
case kLQ2D: // 2D LQ
{
if(species==0||species==2){ // references only for electrons and pions
- Double_t error;
+ Double_t error = 0.;
Double_t point[kNslicesLQ2D];
for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
- AliTRDTKDInterpolator *refLower = dynamic_cast<AliTRDTKDInterpolator*>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
-
- if(refLower){
+ AliTRDTKDInterpolator *refUpper = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ2D)),
+ *refLower = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
+ // Do Interpolation exept for underflow and overflow
+ if(refLower && refUpper){
+ Double_t probLower=0,probUpper=0;
+ refLower->Eval(point,probLower,error);
+ refUpper->Eval(point,probUpper,error);
+ probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
+ } else if(refLower){
+ // underflow
refLower->Eval(point,probLayer,error);
- }
- else {
+ } else if(refUpper){
+ // overflow
+ refUpper->Eval(point,probLayer,error);
+ } else {
AliError("No references available");
}
AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
break;
default:
break;
- }
- return probLayer;
+ }
+ return probLayer;
+ */
+
}
//____________________________________________________________
out[0]=0;
out[1]=0;
for(Int_t islice = 0; islice < nSlice; islice++){
+ if(in[islice]<=0){out[0]=0;out[1]=0;return kFALSE;} // Require that all slices are filled
if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
else out[1]+= in[islice];
}
+ // normalize signal to number of slices
+ out[0]*=1./Double_t(fkPIDResponseObject->GetNSlicesQ0());
+ out[1]*=1./Double_t(nSlice-fkPIDResponseObject->GetNSlicesQ0());
if(out[0] < 1e-6) return kFALSE;
AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
break;
out[0]= 0.;
for(Int_t islice = 0; islice < nSlice; islice++)
if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor; // Protect against negative values for slices having no dE/dx information
+ out[0]*=1./Double_t(nSlice);
if(out[0] < 1e-6) return kFALSE;
AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
break;
// If the function fails when the params are not accessible, the function returns true
//
if(!fkPIDResponseObject){
- AliError("No PID Param object available");
+ AliDebug(3,"No PID Param object available");
return kTRUE;
}
Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
return kTRUE;
}
+
+
+//____________________________________________________________
+Bool_t AliTRDPIDResponse::SetdEdxParams(const AliTRDdEdxParams * par)
+{
+ fkTRDdEdxParams = par;
+ return kTRUE;
+}