#include <TMath.h>
#include <TF1.h>
#include <TGraph.h>
+#include <THnBase.h>
#include <AliVTrack.h>
#include <AliVCluster.h>
#include <AliAODPid.h>
#include "AliDielectronVarManager.h"
+#include "AliDielectronVarCuts.h"
#include "AliDielectronPID.h"
ClassImp(AliDielectronPID)
-TGraph *AliDielectronPID::fgFitCorr=0x0;
+TGraph *AliDielectronPID::fgFitCorr=0x0;
Double_t AliDielectronPID::fgCorr=0.0;
Double_t AliDielectronPID::fgCorrdEdx=1.0;
-TF1 *AliDielectronPID::fgFunEtaCorr=0x0;
-TGraph *AliDielectronPID::fgdEdxRunCorr=0x0;
+TF1 *AliDielectronPID::fgFunEtaCorr=0x0;
+TH1 *AliDielectronPID::fgFunCntrdCorr=0x0;
+TH1 *AliDielectronPID::fgFunWdthCorr=0x0;
+TGraph *AliDielectronPID::fgdEdxRunCorr=0x0;
AliDielectronPID::AliDielectronPID() :
AliAnalysisCuts(),
fSigmaFunLow[icut]=0;
fSigmaFunUp[icut]=0;
fFunSigma[icut]=0x0;
+ fVarCuts[icut]=0x0;
+ fMapElectronCutLow[icut]=0x0;
}
}
fSigmaFunLow[icut]=0;
fSigmaFunUp[icut]=0;
fFunSigma[icut]=0x0;
+ fVarCuts[icut]=0x0;
+ fMapElectronCutLow[icut]=0x0;
}
}
AddCut(det,type,nSigmaLow,nSigmaUp,0,0,exclude,pidBitType,-1);
}
+//______________________________________________
+void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, THnBase * const histLow, Double_t nSigmaUp,
+ Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
+ UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
+{
+ //
+ // cut using a THnBase as a lower cut
+ //
+ if(histLow==0x0){
+ AliError("A valid histogram is required for the lower cut. Not adding the cut!");
+ return;
+ }
+ fMapElectronCutLow[fNcuts]=histLow;
+ AddCut(det,type,0.,nSigmaUp,min,max,exclude,pidBitType,var);
+}
+//______________________________________________
+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)
{
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;
+ // test var range. In case min==max do not cut
+ 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]){
selected = IsSelectedITS(part,icut);
break;
case kTPC:
- selected = IsSelectedTPC(part,icut);
+ selected = IsSelectedTPC(part,icut,values);
break;
case kTRD:
selected = IsSelectedTRD(part,icut);
}
//______________________________________________
-Bool_t AliDielectronPID::IsSelectedTPC(AliVTrack * const part, Int_t icut)
+Bool_t AliDielectronPID::IsSelectedTPC(AliVTrack * const part, Int_t icut, Double_t *values)
{
//
// TPC part of the PID check
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();
Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTPC(part, fPartType[icut]);
+ // eta corrections
if (fPartType[icut]==AliPID::kElectron){
+ // old way 1D
numberOfSigmas-=fgCorr;
+ // via functions (1-3D)
+ numberOfSigmas-=GetCntrdCorr(part);
+ numberOfSigmas/=GetWdthCorr(part);
}
-
- // test if we are supposed to use a function for the cut
- if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom);
- if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom);
-
- Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
+
+ // matching of MC and data //
+
+ // test if we are supposed to use a function for the cut (matching via fits)
+ if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(values[AliDielectronVarManager::kPIn]);
+ if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(values[AliDielectronVarManager::kPIn]);
+
+ // test if we are using cut maps (in calibrated units)
+ Double_t lowElectronCut = fNsigmaLow[icut];
+ Double_t upElectronCut = fNsigmaUp[icut];
+
+ if((fPartType[icut]==AliPID::kElectron) && (fMapElectronCutLow[icut] )) {
+ Double_t *vals = new Double_t[fMapElectronCutLow[icut]->GetNdimensions()];//={-1};
+ // get array of values for the corresponding dimensions using axis names
+ for(Int_t idim=0; idim<fMapElectronCutLow[icut]->GetNdimensions(); idim++) {
+ vals[idim] = values[AliDielectronVarManager::GetValueType(fMapElectronCutLow[icut]->GetAxis(idim)->GetName())];
+ }
+ // find bin for values (w/o creating it in case it is not filled)
+ Long_t bin = fMapElectronCutLow[icut]->GetBin(vals,kFALSE);
+ if(bin>0) lowElectronCut=fMapElectronCutLow[icut]->GetBinContent(bin);
+ else lowElectronCut=100;
+ //printf("low cut %.3f \t for %d dimensional cut map \n",lowElectronCut,fMapElectronCutLow[icut]->GetNdimensions());
+ delete [] vals;
+ }
+
+
+ Bool_t selected=((numberOfSigmas>=lowElectronCut)&&(numberOfSigmas<=upElectronCut))^fExclude[icut];
return selected;
}
else if (def==13) {
// TPC electron inclusion
// TOF electron inclusion if available
- AddCut(kTOF,AliPID::kElectron,-4.,4.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
- AddCut(kTPC,AliPID::kElectron,-3.5,3.5);
+ 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);
}
}
if (!fgFunEtaCorr) return 1;
return fgFunEtaCorr->Eval(track->Eta());
}
+
+//______________________________________________
+Double_t AliDielectronPID::GetPIDCorr(const AliVTrack *track, TH1 *hist)
+{
+ //
+ // return correction value
+ //
+ //TODO: think about taking an values array as argument to reduce # var fills
+
+ //Fill only event and vparticle values (otherwise we end up in a circle)
+ Double_t values[AliDielectronVarManager::kNMaxValues];
+ AliDielectronVarManager::FillVarVParticle(track,values);
+
+ TF1 *fun = (TF1*)hist->GetListOfFunctions()->At(0);
+ Int_t dim=fun->GetNdim();
+
+ Double_t var[3] = {0.,0.,0.};
+ if(dim>0) var[0] = values[hist->GetXaxis()->GetUniqueID()];
+ if(dim>1) var[1] = values[hist->GetYaxis()->GetUniqueID()];
+ if(dim>2) var[2] = values[hist->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