#include "AliKFVertex.h"
#include "AliKFParticle.h"
+using std::cout;
+using std::endl;
+
ClassImp(AliRDHFCutsD0toKpi)
//--------------------------------------------------------------------------
fLowPt(kTRUE),
fDefaultPID(kFALSE),
fUseKF(kFALSE),
-fPtLowPID(2.)
+fPtLowPID(2.),
+fPtMaxSpecialCuts(9999.),
+fmaxPtrackForPID(9999),
+fCombPID(kFALSE),
+fWeightsPositive(new Double_t[AliPID::kSPECIES]),
+fWeightsNegative(new Double_t[AliPID::kSPECIES]),
+fBayesianStrategy(0)
{
//
// Default Constructor
Float_t limits[2]={0,999999999.};
SetPtBins(2,limits);
+ for (Int_t i = 0; i<AliPID::kSPECIES; i++) {
+ fWeightsPositive[i] = 0.;
+ fWeightsNegative[i] = 0.;
+ }
+
}
//--------------------------------------------------------------------------
AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) :
fLowPt(source.fLowPt),
fDefaultPID(source.fDefaultPID),
fUseKF(source.fUseKF),
- fPtLowPID(source.fPtLowPID)
+ fPtLowPID(source.fPtLowPID),
+ fPtMaxSpecialCuts(source.fPtMaxSpecialCuts),
+ fmaxPtrackForPID(source.fmaxPtrackForPID),
+ fCombPID(source.fCombPID),
+ fWeightsPositive(source.fWeightsPositive),
+ fWeightsNegative(source.fWeightsNegative),
+ fBayesianStrategy(source.fBayesianStrategy)
{
//
// Copy constructor
fDefaultPID=source.fDefaultPID;
fUseKF=source.fUseKF;
fPtLowPID=source.fPtLowPID;
-
+ fPtMaxSpecialCuts=source.fPtMaxSpecialCuts;
+ fmaxPtrackForPID=source.fmaxPtrackForPID;
+ fCombPID = source.fCombPID;
+ fWeightsPositive = source.fWeightsPositive;
+ fWeightsNegative = source.fWeightsNegative;
+ fBayesianStrategy = source.fBayesianStrategy;
return *this;
}
+
+//---------------------------------------------------------------------------
+AliRDHFCutsD0toKpi::~AliRDHFCutsD0toKpi()
+{
+//
+// Destructor
+//
+ if (fWeightsNegative) {
+ delete[] fWeightsNegative;
+ fWeightsNegative = 0;
+ }
+ if (fWeightsPositive) {
+ delete[] fWeightsNegative;
+ fWeightsNegative = 0;
+ }
+}
+
//---------------------------------------------------------------------------
void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters,AliAODEvent *aod) {
//
}
//PrintAll();
AliAODRecoDecayHF2Prong* d=(AliAODRecoDecayHF2Prong*)obj;
-
if(!d){
cout<<"AliAODRecoDecayHF2Prong null"<<endl;
return 0;
if(ptD<fMinPtCand) return 0;
if(ptD>fMaxPtCand) return 0;
- if(d->HasBadDaughters()) return 0;
+ if(fUseTrackSelectionWithFilterBits && d->HasBadDaughters()) return 0;
// returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
Int_t returnvaluePID=3;
// call special cuts
Int_t special=1;
- if(fUseSpecialCuts) special=IsSelectedSpecialCuts(d);
+ if(fUseSpecialCuts && (pt<fPtMaxSpecialCuts)) special=IsSelectedSpecialCuts(d);
if(!special) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
// unset recalculated primary vertex when not needed any more
if(selectionLevel==AliRDHFCuts::kAll ||
selectionLevel==AliRDHFCuts::kCandidate ||
selectionLevel==AliRDHFCuts::kPID) {
- returnvaluePID = IsSelectedPID(d);
- fIsSelectedPID=returnvaluePID;
- if(!returnvaluePID) return 0;
+ if (!fCombPID) {
+ returnvaluePID = IsSelectedPID(d);
+ fIsSelectedPID=returnvaluePID;
+ if(!returnvaluePID) return 0;
+ } else {
+ switch (fBayesianStrategy) {
+ case kBayesSimple: {
+ returnvaluePID = IsSelectedSimpleBayesianPID(d);
+ break;
+ }
+ case kBayesWeightNoFilter: {
+ CalculateBayesianWeights(d);
+ returnvaluePID = 3;
+ break;
+ }
+ case(kBayesWeight): {
+ returnvaluePID = IsSelectedCombPID(d);
+ break;
+ }
+ case(kBayesMomentum): {
+ returnvaluePID = IsSelectedCombPID(d);
+ break;
+ }
+ default: {
+ returnvaluePID = IsSelectedCombPID(d);
+ break;
+ }
+ }
+ }
}
+
Int_t returnvalueComb=CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
if(!returnvalueComb) return 0;
// Checking if D0 is in fiducial acceptance region
//
+ if(fMaxRapidityCand>-998.){
+ if(TMath::Abs(y) > fMaxRapidityCand) return kFALSE;
+ else return kTRUE;
+ }
+
if(pt > 5.) {
// applying cut for pt > 5 GeV
AliDebug(2,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt));
if(!fUsePID) return 3;
if(fDefaultPID) return IsSelectedPIDdefault(d);
+ if (fCombPID) return IsSelectedCombPID(d); //to use Bayesian method if applicable
fWhyRejection=0;
Int_t isD0D0barPID[2]={1,2};
Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
checkPIDInfo[daught]=kFALSE;
continue;
}
+ Double_t pProng=aodtrack->P();
// identify kaon
- combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
-
+ if(pProng<fmaxPtrackForPID){
+ combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
+ }
// identify pion
-
- if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
- combinedPID[daught][1]=0;
- }else{
- fPidHF->SetTOF(kFALSE);
- combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
- fPidHF->SetTOF(kTRUE);
- fPidHF->SetCompat(kTRUE);
- }
-
+ if(pProng<fmaxPtrackForPID){
+ if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
+ combinedPID[daught][1]=0;
+ }else{
+ fPidHF->SetTOF(kFALSE);
+ combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
+ fPidHF->SetTOF(kTRUE);
+ fPidHF->SetCompat(kTRUE);
+ }
+ }
if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){ // if not a K- and not a pi- both D0 and D0bar excluded
isD0D0barPID[0]=0;
// identify kaon
combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
- Double_t ptProng=aodtrack->P();
-
- if(ptProng<0.6){
+ if(pProng<0.6){
fPidHF->SetCompat(kFALSE);
combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
fPidHF->SetCompat(kTRUE);
fPidHF->SetSigmaForTPC(sigmaTPCpi);
combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
fPidHF->SetTOF(kTRUE);
- if(ptProng<0.8){
+ if(pProng<0.8){
Bool_t isTPCpion=fPidHF->IsPionRaw(aodtrack,"TPC");
if(isTPCpion){
combinedPID[daught][1]=1;
}
+//---------------------------------------------------------------------------
void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() {
//
//STANDARD CUTS USED FOR 2010 pp analysis
// EVENT CUTS
SetMinVtxContr(1);
-
+ // MAX Z-VERTEX CUT
+ SetMaxVtxZ(10.);
// TRACKS ON SINGLE TRACKS
AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
esdTrackCuts->SetPtRange(0.3,1.e10);
AddTrackCuts(esdTrackCuts);
+
const Int_t nptbins =14;
const Double_t ptmax = 9999.;
pidObj->SetPCompatTOF(1.5);
pidObj->SetSigmaForTPCCompat(3.);
pidObj->SetSigmaForTOFCompat(3.);
+ pidObj->SetOldPid(kTRUE);
SetPidHF(pidObj);
SetUsePID(kTRUE);
SetUseDefaultPID(kFALSE);
+ SetLowPt(kFALSE);
PrintAll();
}
+//---------------------------------------------------------------------------
+void AliRDHFCutsD0toKpi::SetStandardCutsPP2011_276TeV() {
+ //
+ // STANDARD CUTS USED FOR 2011 pp analysis at 2.76TeV
+ //
+
+ SetName("D0toKpiCutsStandard");
+ SetTitle("Standard Cuts for D0 analysis in pp2011 at 2.76TeV run");
+
+ //
+ // Track cuts
+ //
+ AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
+ esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
+ //default
+ esdTrackCuts->SetRequireTPCRefit(kTRUE);
+ esdTrackCuts->SetRequireITSRefit(kTRUE);
+ esdTrackCuts->SetEtaRange(-0.8,0.8);
+ esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
+ AliESDtrackCuts::kAny);
+ // default is kBoth, otherwise kAny
+ esdTrackCuts->SetMinDCAToVertexXY(0.);
+ esdTrackCuts->SetPtRange(0.3,1.e10);
+
+ esdTrackCuts->SetMaxDCAToVertexXY(1.);
+ esdTrackCuts->SetMaxDCAToVertexZ(1.);
+ esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
+
+ AddTrackCuts(esdTrackCuts);
+
+
+ const Int_t nvars=11;
+ const Int_t nptbins=13;
+ Float_t ptbins[nptbins+1];
+ ptbins[0]=0.;
+ ptbins[1]=0.5;
+ ptbins[2]=1.;
+ ptbins[3]=2.;
+ ptbins[4]=3.;
+ ptbins[5]=4.;
+ ptbins[6]=5.;
+ ptbins[7]=6.;
+ ptbins[8]=8.;
+ ptbins[9]=12.;
+ ptbins[10]=16.;
+ ptbins[11]=20.;
+ ptbins[12]=24.;
+ ptbins[13]=9999.;
+
+ SetPtBins(nptbins+1,ptbins);
+
+
+ Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,0.04,0.75,0.3,0.3,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* pt<0.5*/
+ {0.400,0.04,0.75,0.3,0.3,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* 0.5<pt<1*/
+ {0.400,0.03,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.8,0.,0.},/* 1<pt<2 */
+ {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.0003,0.9,0.,0.},/* 2<pt<3 */
+ {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.0002,0.9,0.,0.},/* 3<pt<4 */
+ {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.00015,0.9,0.,0.},/* 4<pt<5 */
+ {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.0001,0.9,0.,0.},/* 5<pt<6 */
+ {0.400,0.09,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* 6<pt<8 */
+ {0.400,0.06,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.00001,0.85,0.,0.},/* 8<pt<12 */
+ {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 12<pt<16 */
+ {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 16<pt<20 */
+ {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 20<pt<24 */
+ {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.}};/* pt>24 */
+
+ //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
+ Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
+ for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
+ for (Int_t ibin=0;ibin<nptbins;ibin++){
+ for (Int_t ivar = 0; ivar<nvars; ivar++){
+ cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
+ }
+ }
+ SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
+ for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
+ delete [] cutsMatrixTransposeStand;
+
+
+ //pid settings
+ AliAODPidHF* pidObj=new AliAODPidHF();
+ Int_t mode=1;
+ const Int_t nlims=2;
+ Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
+ Bool_t compat=kTRUE; //effective only for this mode
+ Bool_t asym=kTRUE;
+ Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
+ pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
+ pidObj->SetMatch(mode);
+ pidObj->SetPLimit(plims,nlims);
+ pidObj->SetSigma(sigmas);
+ pidObj->SetCompat(compat);
+ pidObj->SetTPC(kTRUE);
+ pidObj->SetTOF(kTRUE);
+ pidObj->SetOldPid(kTRUE);
+ SetPidHF(pidObj);
+
+ SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
+
+ SetUsePID(kTRUE);
+
+ //activate pileup rejection (for pp)
+ SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
+
+ //Do recalculate the vertex
+ SetRemoveDaughtersFromPrim(kTRUE);
+
+ // Cut on the zvtx
+ SetMaxVtxZ(10.);
+
+ // Use the kFirst selection for tracks with candidate D with pt<2
+ SetSelectCandTrackSPDFirst(kTRUE,2.);
+
+ // Use special cuts for D candidates with pt<2
+ SetUseSpecialCuts(kTRUE);
+ SetMaximumPtSpecialCuts(2.);
+
+ PrintAll();
+
+ delete pidObj;
+ pidObj=NULL;
+
+ return;
+}
+
+
+//---------------------------------------------------------------------------
void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
//
- //PRELIMINARY CUTS USED FOR 2010 PbPb analysis
- //... EVOLVING SOON
+ // Final CUTS USED FOR 2010 PbPb analysis of central events
+ // REMEMBER TO EVENTUALLY SWITCH ON
+ // SetUseAOD049(kTRUE);
//
SetName("D0toKpiCutsStandard");
SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");
- // PILE UP REJECTION
- //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
+
// CENTRALITY SELECTION
SetMinCentrality(0.);
SetMaxCentrality(80.);
SetUseCentrality(AliRDHFCuts::kCentV0M);
+
// EVENT CUTS
SetMinVtxContr(1);
// MAX Z-VERTEX CUT
SetMaxVtxZ(10.);
-
+ // PILE UP
+ SetOptPileup(AliRDHFCuts::kNoPileupSelection);
+ // PILE UP REJECTION
+ //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
+
// TRACKS ON SINGLE TRACKS
AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
AddTrackCuts(esdTrackCuts);
+ // SPD k FIRST for Pt<3 GeV/c
+ SetSelectCandTrackSPDFirst(kTRUE, 3);
// CANDIDATE CUTS
const Int_t nptbins =13;
SetPtBins(nptbins+1,ptbins);
SetMinPtCandidate(2.);
- Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.85,0.,5.},/* pt<0.5*/
- {0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-35000.*1E-8,0.9,0.
,5.},/* 0.5<pt<1*/
- {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-43000.*1E-8,0.85,0.
,5.},/* 1<pt<2 */
- {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-4
0000.*1E-8,0.95,0.998,5.},/* 2<pt<3 */
+ Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.85,0.99,2.},/* pt<0.5*/
+ {0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-35000.*1E-8,0.9,0.
99,2.},/* 0.5<pt<1*/
+ {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-43000.*1E-8,0.85,0.
995,8.},/* 1<pt<2 */
+ {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-4
5000.*1E-8,0.95,0.998,7.},/* 2<pt<3 */
{0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-36000.*1E-8,0.95,0.998,5.},/* 3<pt<4 */
{0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-27000.*1E-8,0.95,0.998,5.},/* 4<pt<5 */
{0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-21000.*1E-8,0.92,0.998,5.},/* 5<pt<6 */
{0.400,270.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-14000.*1E-8,0.88,0.998,5.},/* 6<pt<8 */
{0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.998,5.},/* 8<pt<12 */
- {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.83,0.998,
5.},/* 12<pt<16 */
- {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.82,0.99
8,5.},/* 16<pt<20 */
- {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.81,0.99
8,5.},/* 20<pt<24 */
- {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.8,0.998,5.}};/* pt>24 */
+ {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.83,0.998,
8.},/* 12<pt<16 */
+ {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.82,0.99
5,6.},/* 16<pt<20 */
+ {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.81,0.99
5,6.},/* 20<pt<24 */
+ {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.8,0.99,2.}};/* pt>24 */
//CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
pidObj->SetPCompatTOF(2.);
pidObj->SetSigmaForTPCCompat(3.);
pidObj->SetSigmaForTOFCompat(3.);
+ pidObj->SetOldPid(kTRUE);
+
+
+ SetPidHF(pidObj);
+ SetUsePID(kTRUE);
+ SetUseDefaultPID(kFALSE);
+
+
+ PrintAll();
+
+
+ delete pidObj;
+ pidObj=NULL;
+
+ return;
+
+}
+
+//---------------------------------------------------------------------------
+void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010Peripherals() {
+ // CUTS USED FOR D0 RAA for the CENTRALITY RANGE 20-80%
+
+
+ SetName("D0toKpiCutsStandard");
+ SetTitle("Standard Cuts for D0 analysis in PbPb2010 run, for peripheral events");
+
+
+ // CENTRALITY SELECTION
+ SetMinCentrality(40.);
+ SetMaxCentrality(80.);
+ SetUseCentrality(AliRDHFCuts::kCentV0M);
+
+ // EVENT CUTS
+ SetMinVtxContr(1);
+ // MAX Z-VERTEX CUT
+ SetMaxVtxZ(10.);
+ // PILE UP
+ SetOptPileup(AliRDHFCuts::kNoPileupSelection);
+ // PILE UP REJECTION
+ //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
+
+ // TRACKS ON SINGLE TRACKS
+ AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
+ esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
+ esdTrackCuts->SetRequireTPCRefit(kTRUE);
+ esdTrackCuts->SetRequireITSRefit(kTRUE);
+ // esdTrackCuts->SetMinNClustersITS(4);
+ esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
+ esdTrackCuts->SetMinDCAToVertexXY(0.);
+ esdTrackCuts->SetEtaRange(-0.8,0.8);
+ esdTrackCuts->SetPtRange(0.5,1.e10);
+
+ esdTrackCuts->SetMaxDCAToVertexXY(1.);
+ esdTrackCuts->SetMaxDCAToVertexZ(1.);
+ esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0025*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
+
+
+ AddTrackCuts(esdTrackCuts);
+ // SPD k FIRST for Pt<3 GeV/c
+ SetSelectCandTrackSPDFirst(kTRUE, 3);
+
+ // CANDIDATE CUTS
+ const Int_t nptbins =13;
+ const Double_t ptmax = 9999.;
+ const Int_t nvars=11;
+ Float_t ptbins[nptbins+1];
+ ptbins[0]=0.;
+ ptbins[1]=0.5;
+ ptbins[2]=1.;
+ ptbins[3]=2.;
+ ptbins[4]=3.;
+ ptbins[5]=4.;
+ ptbins[6]=5.;
+ ptbins[7]=6.;
+ ptbins[8]=8.;
+ ptbins[9]=12.;
+ ptbins[10]=16.;
+ ptbins[11]=20.;
+ ptbins[12]=24.;
+ ptbins[13]=ptmax;
+ SetGlobalIndex(nvars,nptbins);
+ SetPtBins(nptbins+1,ptbins);
+ SetMinPtCandidate(0.);
+
+ Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-20000.*1E-8,0.7,0.993,2.},/* pt<0.5*/
+ {0.400,400.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.85,0.993,2.},/* 0.5<pt<1*/
+ {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.85,0.995,6.},/* 1<pt<2 */
+ {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.95,0.991,5.},/* 2<pt<3 */
+ {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-36000.*1E-8,0.95,0.993,5.},/* 3<pt<4 */
+ {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-27000.*1E-8,0.95,0.995,5.},/* 4<pt<5 */
+ {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-21000.*1E-8,0.92,0.998,5.},/* 5<pt<6 */
+ {0.400,270.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-14000.*1E-8,0.88,0.998,5.},/* 6<pt<8 */
+ {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.995,5.},/* 8<pt<12 */
+ {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.83,0.995,4.},/* 12<pt<16 */
+ {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.82,0.995,4.},/* 16<pt<20 */
+ {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.81,0.995,4.},/* 20<pt<24 */
+ {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.8,0.995,4.}};/* pt>24 */
+
+
+ //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
+ Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
+ for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
+
+ for (Int_t ibin=0;ibin<nptbins;ibin++){
+ for (Int_t ivar = 0; ivar<nvars; ivar++){
+ cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
+ }
+ }
+
+ SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
+ SetUseSpecialCuts(kTRUE);
+ SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
+ for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
+ delete [] cutsMatrixTransposeStand;
+ cutsMatrixTransposeStand=NULL;
+
+ // PID SETTINGS
+ AliAODPidHF* pidObj=new AliAODPidHF();
+ //pidObj->SetName("pid4D0");
+ Int_t mode=1;
+ const Int_t nlims=2;
+ Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
+ Bool_t compat=kTRUE; //effective only for this mode
+ Bool_t asym=kTRUE;
+ Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
+ pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
+ pidObj->SetMatch(mode);
+ pidObj->SetPLimit(plims,nlims);
+ pidObj->SetSigma(sigmas);
+ pidObj->SetCompat(compat);
+ pidObj->SetTPC(kTRUE);
+ pidObj->SetTOF(kTRUE);
+ pidObj->SetPCompatTOF(2.);
+ pidObj->SetSigmaForTPCCompat(3.);
+ pidObj->SetSigmaForTOFCompat(3.);
+ pidObj->SetOldPid(kTRUE);
SetPidHF(pidObj);
SetUsePID(kTRUE);
SetUseDefaultPID(kFALSE);
+ SetLowPt(kTRUE,2.);
PrintAll();
pidObj=NULL;
return;
+
+}
+
+
+//---------------------------------------------------------------------------
+void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2011() {
+
+ // Default 2010 PbPb cut object
+ SetStandardCutsPbPb2010();
+ AliAODPidHF *pidobj=GetPidHF();
+
+ pidobj->SetOldPid(kFALSE);
+
+ //
+ // Enable all 2011 PbPb run triggers
+ //
+ SetTriggerClass("");
+ ResetMaskAndEnableMBTrigger();
+ EnableCentralTrigger();
+ EnableSemiCentralTrigger();
}
+
+//---------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::IsSelectedCombPID(AliAODRecoDecayHF* d)
+{
+ // ############################################################
+ //
+ // Apply Bayesian PID selection
+ // Implementation of Bayesian PID: Jeremy Wilkinson
+ //
+ // ############################################################
+
+ if (!fUsePID || !d) return 3;
+
+ switch (fBayesianStrategy) {
+
+ case kBayesWeightNoFilter: {
+ CalculateBayesianWeights(d); //If weighted, no filtering: Calculate weights, return as unidentified
+ return 3;
+ break;
+ }
+ case kBayesSimple: {
+ return IsSelectedSimpleBayesianPID(d); // If simple, go to simple method
+ break;
+ }
+ case(kBayesWeight): {
+ break;
+ }
+ case(kBayesMomentum): {
+ break; //If weighted or momentum method, stay in this function
+ }
+
+ }
+
+ // Int_t isD0D0barPID[2]={1,2};
+
+
+
+
+ Int_t isD0 = 0;
+ Int_t isD0bar = 0;
+ Int_t returnvalue = 0;
+
+ Int_t isPosKaon = 0, isNegKaon = 0;
+ // Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
+ // same for prong 2
+ // values convention 0 = not identified (but compatible) || No PID (->hasPID flag)
+ // 1 = identified
+ // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
+ // Initial hypothesis: unknown (but compatible)
+
+
+ // combinedPID[0][0] = 0; // First daughter, kaon
+ // combinedPID[0][1] = 0; // First daughter, pion
+ // combinedPID[1][0] = 0; // Second daughter, kaon
+ // combinedPID[1][1] = 0; // Second daughter, pion
+
+ Bool_t checkPIDInfo[2] = {kTRUE, kTRUE};
+ AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
+
+ if ((aodtrack[0]->Charge()*aodtrack[1]->Charge()) != -1) {
+ return 0; //reject if charges not opposite
+ }
+ for (Int_t daught = 0; daught < 2; daught++) {
+ //Loop over prongs
+
+ if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
+ checkPIDInfo[daught] = kFALSE;
+ continue;
+ }
+ CalculateBayesianWeights(d);
+ // Double_t prob0[AliPID::kSPECIES];
+ // fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), prob0);
+ ///Three possible Bayesian methods: Pick one (only first implemented so far. TODO: implement other two methods)
+ ///A: Standard max. probability method
+ if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kKaon]) { //If highest probability lies with kaon
+ isPosKaon = 1; //flag [daught] as a kaon
+ }
+
+ if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kKaon]) { //If highest probability lies with kaon
+ isNegKaon = 1; //flag [daught] as a kaon
+ }
+ ///B: Method based on probability greater than prior
+ /* Double_t kaonpriorprob;
+ TH1D *kaonprior = fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon);
+ kaonpriorprob = kaonprior->GetBinContent(kaonprior->FindBin(aodtrack[daught]->P()));
+
+ if (prob0[AliPID::kKaon] > kaonpriorprob){combinedPID[daught][0] = 1;}
+
+ */
+ ///C: Probability threshold (TODO: Add necessary threshold variables; currently unused)
+ /*
+ if(prob0[AliPID::kKaon] >= fMinProbKaon) {combinedPID[daught][0] = 1;}
+
+ */
+
+ }
+
+
+ //Loop over daughters ends here
+
+ if (!checkPIDInfo[0] && !checkPIDInfo[1]) {
+ return 0; //Reject if both daughters lack both TPC and TOF info
+ }
+ //Momentum-based selection
+
+
+ if (isNegKaon && isPosKaon) { // If both are kaons, reject
+ isD0 = 0;
+ isD0bar = 0;
+ } else if (isNegKaon) { //If negative kaon present, D0
+ isD0 = 1;
+ } else if (isPosKaon) { //If positive kaon present, D0bar
+ isD0bar = 1;
+ } else { //If neither ID'd as kaon, subject to extra tests
+ isD0 = 1;
+ isD0bar = 1;
+ if (aodtrack[0]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[0], "TOF"))) { //If it's low momentum and not ID'd as a kaon, we assume it must be a pion
+ if (aodtrack[0]->Charge() == -1) {
+ isD0 = 0; //Check charge and reject based on pion hypothesis
+ }
+ if (aodtrack[0]->Charge() == 1) {
+ isD0bar = 0;
+ }
+ }
+ if (aodtrack[1]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[1], "TOF"))) {
+ if (aodtrack[1]->Charge() == -1) {
+ isD0 = 0;
+ }
+ if (aodtrack[1]->Charge() == 1) {
+ isD0bar = 0;
+ }
+ }
+ }
+
+
+
+ if (isD0 && isD0bar) {
+ returnvalue = 3;
+ }
+ if (isD0 && !isD0bar) {
+ returnvalue = 1;
+ }
+ if (!isD0 && isD0bar) {
+ returnvalue = 2;
+ }
+ if (!isD0 && !isD0bar) {
+ returnvalue = 0;
+ }
+
+ return returnvalue;
+
+
+
+ /*
+
+ //OLD CODE
+ if (combinedPID[0][0] == 1 && combinedPID[1][1] == 1) { //First track is kaon, second is pion
+ if (aodtrack[0]->Charge() == -1) isD0 = 1; //Kaon negative: particle is D0
+ else isD0bar = 1; //Kaon positive: particle is anti
+ } else if (combinedPID[1][0] == 1 && combinedPID[0][1] == 1) { //First track is pion, second is kaon
+ if (aodtrack[1]->Charge() == -1) isD0 = 1; //Kaon negative: particle is D0
+ else isD0bar = 1; //Kaon positive: particle is anti
+ }
+
+ else if (combinedPID[0][0] == 1 && combinedPID[1][0] == 1) {
+ isD0 = 1; //If both are kaons, leave to top. cuts
+ isD0bar = 1;
+ }
+
+
+ else {
+ isD0 = 1;
+ isD0bar = 1;
+ if (combinedPID[0][0] == 0 && aodtrack[0]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[0], "TOF"))) { //If it's low momentum and not ID'd as a kaon, we assume it must be a pion
+ if (aodtrack[0]->Charge() == -1) {
+ isD0 = 0; //Check charge and reject depending on pion hypothesis
+ }
+ if (aodtrack[0]->Charge() == 1) {
+ isD0bar = 0;
+ }
+ }
+ if (combinedPID[1][0] == 0 && aodtrack[1]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[1], "TOF"))) {
+ if (aodtrack[1]->Charge() == -1) {
+ isD0 = 0;
+ }
+ if (aodtrack[1]->Charge() == 1) {
+ isD0bar = 0;
+ }
+ }
+ }*/
+ /* if(combinedPID[0][1] == 1 && combinedPID[1][1] == 1) { //If both are ID'd as pions:
+ if ((aodtrack[0]->P()) <1.0 && (aodtrack[1]->P() < 1.0)) { //Check both momenta. If both low, we require one kaon.
+ isD0 == 0; isD0bar == 0;} //Both pions, both low momenta: reject
+ else {isD0 == 1; isD0bar == 1;} //Both pions, both high momenta: Allow as both (filter with topological cuts)
+ }
+ */
+
+
+}
+
+
+//---------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::IsSelectedSimpleBayesianPID(AliAODRecoDecayHF* d)
+
+{
+ //Apply Bayesian PID selection; "simple" method.
+ //Looks for a kaon via max. probability. If neither daughter is a kaon, candidate is rejected.
+
+ Int_t isPosKaon = 0, isNegKaon = 0;
+ Int_t returnvalue; //0 = rejected, 1 = D0, 2 = D0bar. This version will not output 3 (both).
+ Bool_t checkPIDInfo[2] = {kTRUE, kTRUE};
+
+ AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
+ if ((aodtrack[0]->Charge() * aodtrack[1]->Charge()) != -1) {
+ return 0; //Reject if charges do not oppose
+ }
+ for (Int_t daught = 0; daught < 2; daught++) {
+ //Loop over prongs to check PID information
+
+ if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
+ checkPIDInfo[daught] = kFALSE;
+ }
+
+ }
+ //Loop over daughters ends here.
+
+ if (!checkPIDInfo[0] && !checkPIDInfo[1]) {
+ return 0;
+ }
+
+ CalculateBayesianWeights(d);
+
+
+ // identify kaon
+ ///Three possible Bayesian methods: Pick one (only first implemented so far. TODO: implement other two methods)
+ ///A: Standard max. probability method
+
+
+ if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kKaon]) { //If highest probability lies with kaon
+ isPosKaon = 1; //flag [daught] as a kaon
+ }
+
+ if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kKaon]) { //If highest probability lies with kaon
+ isNegKaon = 1; //flag [daught] as a kaon
+ }
+
+
+ ///B: Method based on probability greater than prior
+ /* Double_t kaonpriorprob;
+ TH1D *kaonprior = fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon);
+ kaonpriorprob = kaonprior->GetBinContent(kaonprior->FindBin(aodtrack[daught]->P()));
+
+ if (prob0[AliPID::kKaon] > kaonpriorprob){combinedPID[daught][0] = 1;}
+
+ */
+ ///C: Probability threshold (TODO: Add necessary threshold variables; currently unused)
+ /*
+ if(prob0[AliPID::kKaon] >= fMinProbKaon) {combinedPID[daught][0] = 1;}
+
+ */
+
+
+
+
+ if (isPosKaon && isNegKaon) { //If both are ID'd as kaons, reject
+ returnvalue = 0;
+ } else if (isNegKaon) { //If negative kaon, D0
+ returnvalue = 1;
+ } else if (isPosKaon) { //If positive kaon, D0-bar
+ returnvalue = 2;
+ } else {
+ returnvalue = 0; //If neither kaon, reject
+ }
+
+ return returnvalue;
+}
+
+
+
+
+//---------------------------------------------------------------------------
+void AliRDHFCutsD0toKpi::CalculateBayesianWeights(AliAODRecoDecayHF* d)
+
+{
+ //Function to compute weights for Bayesian method
+ Double_t *prob0;
+ prob0 = new Double_t[AliPID::kSPECIES];
+
+
+ AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
+ if ((aodtrack[0]->Charge() * aodtrack[1]->Charge()) != -1) {
+ return; //Reject if charges do not oppose
+ }
+ for (Int_t daught = 0; daught < 2; daught++) {
+ //Loop over prongs
+
+ if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
+ continue;
+ }
+
+
+ // identify kaon
+ fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), prob0);
+
+ if (aodtrack[daught]->Charge() == +1) {
+ SetWeightsPositive(prob0);
+ }
+ if (aodtrack[daught]->Charge() == -1) {
+ SetWeightsNegative(prob0);
+ }
+
+ }
+}
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff