#include <AliLog.h>
#include <AliExternalTrackParam.h>
#include <AliPIDResponse.h>
+#include <AliTRDPIDResponse.h>
#include <AliESDtrack.h> //!!!!! Remove once Eta correction is treated in the tender
+#include <AliAODTrack.h>
+#include <AliAODPid.h>
#include "AliDielectronVarManager.h"
+#include "AliDielectronVarCuts.h"
#include "AliDielectronPID.h"
Double_t AliDielectronPID::fgCorr=0.0;
Double_t AliDielectronPID::fgCorrdEdx=1.0;
TF1 *AliDielectronPID::fgFunEtaCorr=0x0;
+TF1 *AliDielectronPID::fgFunCntrdCorr=0x0;
+TF1 *AliDielectronPID::fgFunWdthCorr=0x0;
TGraph *AliDielectronPID::fgdEdxRunCorr=0x0;
AliDielectronPID::AliDielectronPID() :
fFunLowerCut[icut]=0x0;
fRequirePIDbit[icut]=0;
fActiveCuts[icut]=-1;
+ fSigmaFunLow[icut]=0;
+ fSigmaFunUp[icut]=0;
+ fFunSigma[icut]=0x0;
+ fVarCuts[icut]=0x0;
}
}
fFunLowerCut[icut]=0x0;
fRequirePIDbit[icut]=0;
fActiveCuts[icut]=0;
+ fSigmaFunLow[icut]=0;
+ fSigmaFunUp[icut]=0;
+ fFunSigma[icut]=0x0;
+ fVarCuts[icut]=0x0;
}
}
fRequirePIDbit[fNcuts]=pidBitType;
fActiveCuts[fNcuts]=(var==-1 ? AliDielectronVarManager::kP : var);
- // AliInfo(Form("Add PID cut %d: sigma [% .1f,% .1f] \t cut [% .1f,% .f] \t var %d->%s \n",
- // fNcuts,nSigmaLow,nSigmaUp,min,max,fActiveCuts[fNcuts],AliDielectronVarManager::GetValueName(fActiveCuts[fNcuts])));
+ AliDebug(1,Form("Add PID cut %d: sigma [% .1f,% .1f] \t cut [% .1f,% .f] \t var %d->%s \n",
+ fNcuts,nSigmaLow,nSigmaUp,min,max,fActiveCuts[fNcuts],AliDielectronVarManager::GetValueName(fActiveCuts[fNcuts])));
++fNcuts;
AddCut(det,type,0.,0.,min,max,exclude,pidBitType,var);
}
+//______________________________________________
+void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp,
+ Double_t min, Double_t max, Bool_t exclude,
+ UInt_t pidBitType, TF1 * const funSigma)
+{
+ //
+ // cut using a TF1 as lower cut
+ //
+ if (funSigma==0x0){
+ AliError("A valid function is required for the sigma cut. Not adding the cut!");
+ return;
+ }
+ fFunSigma[fNcuts]=funSigma;
+ fSigmaFunLow[fNcuts]=min;
+ fSigmaFunUp[fNcuts]=max;
+
+ AddCut(det,type,nSigmaLow,nSigmaUp,0,0,exclude,pidBitType,-1);
+}
+
+//______________________________________________
+void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp,
+ AliDielectronVarCuts *var, Bool_t exclude/*=kFALSE*/,
+ UInt_t pidBitType/*=AliDielectronPID::kRequire*/)
+{
+ //
+ // Add a pid nsigma cut
+ // use response of detector 'det' in the ranges for variables defined in var
+ // use a sigma band between 'nSigmaLow' and 'nSigmaUp'
+ // if nSigmaUp==-99999. then nSigmaLow will be uesd as a symmetric band +- nSigmaLow
+ // specify whether to 'exclude' the given band
+ //
+ if(!var) return;
+ if (fNcuts>=kNmaxPID){
+ AliError(Form("only %d pid cut ranges allowed",kNmaxPID));
+ return;
+ }
+ if (TMath::Abs(nSigmaUp+99999.)<1e-20){
+ nSigmaUp=TMath::Abs(nSigmaLow);
+ nSigmaLow=-1.*nSigmaUp;
+ }
+ fDetType[fNcuts]=det;
+ fPartType[fNcuts]=type;
+ fNsigmaLow[fNcuts]=nSigmaLow;
+ fNsigmaUp[fNcuts]=nSigmaUp;
+ fExclude[fNcuts]=exclude;
+ fRequirePIDbit[fNcuts]=pidBitType;
+ fVarCuts[fNcuts]=var;
+
+ AliDebug(1,Form("Add PID cut %d: sigma [% .1f,% .1f] \n",
+ fNcuts,nSigmaLow,nSigmaUp));
+
+ ++fNcuts;
+
+}
+
//______________________________________________
Bool_t AliDielectronPID::IsSelected(TObject* track)
{
//loop over all cuts
AliVTrack *part=static_cast<AliVTrack*>(track);
AliESDtrack *esdTrack=0x0;
+ AliAODTrack *aodTrack=0x0;
Double_t origdEdx=-1;
// apply ETa correction, remove once this is in the tender
- if( (part->IsA() == AliESDtrack::Class()) && fgFunEtaCorr ){
+ if( (part->IsA() == AliESDtrack::Class()) ){
esdTrack=static_cast<AliESDtrack*>(part);
origdEdx=esdTrack->GetTPCsignal();
esdTrack->SetTPCsignal(origdEdx/GetEtaCorr(esdTrack)/fgCorrdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
+ } else if ( (part->IsA() == AliAODTrack::Class()) ){
+ aodTrack=static_cast<AliAODTrack*>(track);
+ AliAODPid *pid=const_cast<AliAODPid*>(aodTrack->GetDetPid());
+ if (pid){
+ origdEdx=pid->GetTPCsignal();
+ pid->SetTPCsignal(origdEdx/GetEtaCorr(aodTrack)/fgCorrdEdx);
+ }
}
-
//Fill values
Double_t values[AliDielectronVarManager::kNMaxValues];
Double_t min=fmin[icut];
Double_t max=fmax[icut];
Double_t val=values[fActiveCuts[icut]];
-
+
// test var range. In case min==max do not cut
- if ( (TMath::Abs(min-max)>1e-20) && (val<min || val>=max) ) continue;
-
+ if ( fVarCuts[icut] ) {
+ if ( !fVarCuts[icut]->IsSelected(part) ) continue;
+ }
+ else if ( ( (TMath::Abs(min-max)>1e-20) && (val<min || val>=max) ) ) {
+ continue;
+ }
+
+ // check if fFunSigma is set, then check if 'part' is in sigma range of the function
+ if(fFunSigma[icut]){
+ val= fPIDResponse->NumberOfSigmasTPC(part, fPartType[icut]);
+ if (fPartType[icut]==AliPID::kElectron){
+ val-=fgCorr;
+ }
+ min= fFunSigma[icut]->Eval(part->GetTPCmomentum())+fSigmaFunLow[icut];
+ max= fFunSigma[icut]->Eval(part->GetTPCmomentum())+fSigmaFunUp[icut];
+ if(val<min || val>=max) continue;
+ }
+
switch (fDetType[icut]){
case kITS:
selected = IsSelectedITS(part,icut);
case kTRDeleEff:
selected = IsSelectedTRDeleEff(part,icut);
break;
+ case kTRDeleEff2D:
+ selected = IsSelectedTRDeleEff(part,icut,AliTRDPIDResponse::kLQ2D);
+ break;
case kTOF:
selected = IsSelectedTOF(part,icut);
break;
}
if (!selected) {
if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
+ else if (aodTrack){
+ AliAODPid *pid=const_cast<AliAODPid*>(aodTrack->GetDetPid());
+ if (pid) pid->SetTPCsignal(origdEdx);
+ }
return kFALSE;
}
}
if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
+ else if (aodTrack){
+ AliAODPid *pid=const_cast<AliAODPid*>(aodTrack->GetDetPid());
+ if (pid) pid->SetTPCsignal(origdEdx);
+ }
return selected;
}
// ITS part of the PID check
// Don't accept the track if there was no pid bit set
//
-
- if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kITSpid)) return kFALSE;
- if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kITSpid)) return kTRUE;
+ AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kITS,part);
+ if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
+ if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
Double_t mom=part->P();
// TPC part of the PID check
// Don't accept the track if there was no pid bit set
//
- if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTPCpid)) return kFALSE;
- if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTPCpid)) return kTRUE;
+ AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC,part);
+ if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
+ if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
Double_t mom=part->GetTPCmomentum();
if (fPartType[icut]==AliPID::kElectron){
numberOfSigmas-=fgCorr;
+ numberOfSigmas-=GetCntrdCorr(part);
+ numberOfSigmas/=GetWdthCorr(part);
}
// test if we are supposed to use a function for the cut
// the TRD checks on the probabilities.
//
- if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kFALSE;
- if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kTRUE;
+ AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTRD,part);
+ if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
+ if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
+
if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable && (part->GetTRDntrackletsPID()<4)) return kTRUE;
Double_t p[AliPID::kSPECIES]={0.};
}
//______________________________________________
-Bool_t AliDielectronPID::IsSelectedTRDeleEff(AliVTrack * const part, Int_t icut)
+Bool_t AliDielectronPID::IsSelectedTRDeleEff(AliVTrack * const part, Int_t icut, AliTRDPIDResponse::ETRDPIDMethod PIDmethod)
{
//
// TRD part of the pid check using electron efficiency requirement
// and the lower limit regarded as the requested electron efficiency
//
- if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kFALSE;
- if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kTRUE;
+ AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTRD,part);
+ if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
+ if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
+
+ Double_t centrality = -1.;
+ if(part->IsA() == AliESDtrack::Class())
+ centrality=(const_cast<AliESDEvent*>( (static_cast<const AliESDtrack*>(part))->GetESDEvent()) )->GetCentrality()->GetCentralityPercentile("V0M");
+ if(part->IsA() == AliAODTrack::Class())
+ centrality=(const_cast<AliAODEvent*>( (static_cast<const AliAODTrack*>(part))->GetAODEvent()) )->GetCentrality()->GetCentralityPercentile("V0M");
- Bool_t selected=fPIDResponse->IdentifiedAsElectronTRD(part,fNsigmaLow[icut])^fExclude[icut];
+ Bool_t selected=fPIDResponse->IdentifiedAsElectronTRD(part,fNsigmaLow[icut], centrality, PIDmethod)^fExclude[icut];
return selected;
}
// TOF part of the PID check
// Don't accept the track if there was no pid bit set
//
-
- if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTOFpid)) return kFALSE;
- if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTOFpid)) return kTRUE;
+ AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,part);
+ if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
+ if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTOF(part, fPartType[icut]);
AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
}
+ else if (def==13) {
+ // TPC electron inclusion
+ // TOF electron inclusion if available
+ AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
+ AddCut(kTPC,AliPID::kElectron,10.);
+ }
+ else if (def==14) {
+ // TRD 1D 90% elec eff, 4-6 tracklets
+ // TPC electron inclusion
+ // TOF electron inclusion if available
+ AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,3.5,6.,kFALSE,
+ AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDpidQuality);
+ AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
+ AddCut(kTPC,AliPID::kElectron,10.);
+ }
+ else if (def==15) {
+ // TRD 1D 90% elec eff, 4-6 tracklets, chi2 < 2
+ // TPC electron inclusion
+ // TOF electron inclusion if available
+ AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,3.5,6.,kFALSE,
+ AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDpidQuality);
+ AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,0.,2.,kFALSE,
+ AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDchi2);
+ AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
+ AddCut(kTPC,AliPID::kElectron,10.);
+ }
+ else if (def==16) {
+ // TPC electron inclusion
+ // TOF electron inclusion if available
+ AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
+ AddCut(kTPC,AliPID::kElectron,-3.5,+3.5);
+ }
}
}
//______________________________________________
-Double_t AliDielectronPID::GetEtaCorr(AliVTrack *track)
+Double_t AliDielectronPID::GetEtaCorr(const AliVTrack *track)
{
//
// return eta correction
if (!fgFunEtaCorr) return 1;
return fgFunEtaCorr->Eval(track->Eta());
}
+
+//______________________________________________
+void AliDielectronPID::SetCentroidCorrFunction(TF1 *fun, UInt_t varx, UInt_t vary, UInt_t varz)
+{
+ fun->GetHistogram()->GetXaxis()->SetUniqueID(varx);
+ fun->GetHistogram()->GetYaxis()->SetUniqueID(vary);
+ fun->GetHistogram()->GetZaxis()->SetUniqueID(varz);
+ fgFunCntrdCorr=fun;
+}
+//______________________________________________
+void AliDielectronPID::SetWidthCorrFunction(TF1 *fun, UInt_t varx, UInt_t vary, UInt_t varz)
+{
+ fun->GetHistogram()->GetXaxis()->SetUniqueID(varx);
+ fun->GetHistogram()->GetYaxis()->SetUniqueID(vary);
+ fun->GetHistogram()->GetZaxis()->SetUniqueID(varz);
+ fgFunWdthCorr=fun;
+}
+
+//______________________________________________
+Double_t AliDielectronPID::GetPIDCorr(const AliVTrack *track, TF1 *fun)
+{
+ //
+ // return correction value
+ //
+
+ //Fill only event and vparticle values (otherwise we end up in a circle)
+ Double_t values[AliDielectronVarManager::kNMaxValues];
+ AliDielectronVarManager::FillVarVParticle(track,values);
+
+ Int_t dim=fun->GetNdim();
+ Double_t var[3] = {0.,0.,0.};
+ if(dim>0) var[0] = values[fun->GetHistogram()->GetXaxis()->GetUniqueID()];
+ if(dim>1) var[1] = values[fun->GetHistogram()->GetYaxis()->GetUniqueID()];
+ if(dim>2) var[2] = values[fun->GetHistogram()->GetZaxis()->GetUniqueID()];
+ Double_t corr = fun->Eval(var[0],var[1],var[2]);
+ // printf(" %d-dim CORR value: %f (track %p) \n",dim,corr,track);
+ return corr;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff