#include "AliKFVertex.h"
#include "AliKFParticle.h"
+using std::cout;
+using std::endl;
+
ClassImp(AliRDHFCutsD0toKpi)
//--------------------------------------------------------------------------
-AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const char* name) :
-AliRDHFCuts(name),
-fUseSpecialCuts(kFALSE),
-fLowPt(kTRUE),
-fDefaultPID(kFALSE),
-fUseKF(kFALSE),
-fPtLowPID(2.),
-fPtMaxSpecialCuts(9999.)
+AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const char* name) :
+ AliRDHFCuts(name),
+ fUseSpecialCuts(kFALSE),
+ fLowPt(kTRUE),
+ fDefaultPID(kFALSE),
+ fUseKF(kFALSE),
+ fPtLowPID(2.),
+ fPtMaxSpecialCuts(9999.),
+ fmaxPtrackForPID(9999),
+ fCombPID(kFALSE),
+ fnSpecies(AliPID::kSPECIES),
+ fWeightsPositive(0),
+ fWeightsNegative(0),
+ fProbThreshold(0.5),
+ fBayesianStrategy(0),
+ fBayesianCondition(0)
{
//
// Default Constructor
Float_t limits[2]={0,999999999.};
SetPtBins(2,limits);
+ fWeightsNegative = new Double_t[AliPID::kSPECIES];
+ fWeightsPositive = new Double_t[AliPID::kSPECIES];
+
+ for (Int_t i = 0; i<AliPID::kSPECIES; i++) {
+ fWeightsPositive[i] = 0.;
+ fWeightsNegative[i] = 0.;
+ }
+
}
//--------------------------------------------------------------------------
AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) :
fDefaultPID(source.fDefaultPID),
fUseKF(source.fUseKF),
fPtLowPID(source.fPtLowPID),
- fPtMaxSpecialCuts(source.fPtMaxSpecialCuts)
+ fPtMaxSpecialCuts(source.fPtMaxSpecialCuts),
+ fmaxPtrackForPID(source.fmaxPtrackForPID),
+ fCombPID(source.fCombPID),
+ fnSpecies(source.fnSpecies),
+ fWeightsPositive(source.fWeightsPositive),
+ fWeightsNegative(source.fWeightsNegative),
+ fProbThreshold(source.fProbThreshold),
+ fBayesianStrategy(source.fBayesianStrategy),
+ fBayesianCondition(source.fBayesianCondition)
{
//
// Copy constructor
fUseKF=source.fUseKF;
fPtLowPID=source.fPtLowPID;
fPtMaxSpecialCuts=source.fPtMaxSpecialCuts;
-
+ fmaxPtrackForPID=source.fmaxPtrackForPID;
+ fCombPID = source.fCombPID;
+ fnSpecies = source.fnSpecies;
+ fWeightsPositive = source.fWeightsPositive;
+ fWeightsNegative = source.fWeightsNegative;
+ fProbThreshold = source.fProbThreshold;
+ fBayesianStrategy = source.fBayesianStrategy;
+ fBayesianCondition = source.fBayesianCondition;
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;
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 {
+ returnvaluePID = IsSelectedCombPID(d);
+ if(!returnvaluePID) return 0;
+ }
+ }
+
+
Int_t returnvalueComb=CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
// 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;
Bool_t checkPIDInfo[2]={kTRUE,kTRUE};
Double_t sigma_tmp[3]={fPidHF->GetSigma(0),fPidHF->GetSigma(1),fPidHF->GetSigma(2)};
+ Bool_t isTOFused=fPidHF->GetTOF(),isCompat=fPidHF->GetCompat();
for(Int_t daught=0;daught<2;daught++){
//Loop con prongs
AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
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);
+ if(isTOFused)fPidHF->SetTOF(kTRUE);
+ if(isCompat)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;
isD0D0barPID[1]=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
esdTrackCuts->SetPtRange(0.3,1.e10);
AddTrackCuts(esdTrackCuts);
-
+ delete esdTrackCuts;
+ esdTrackCuts=NULL;
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);
}
+//___________________________________________________________________________
+void AliRDHFCutsD0toKpi::SetStandardCutsPP2010vsMult() {
+ //
+ // Cuts for 2010 pp 7 TeV data analysis in Ntracklets bins (vs multiplicity)
+ //
+ SetName("D0toKpiCuts");
+ SetTitle("Cuts for D0 analysis in 2010-data pp 7 TeV vs multiplicity");
+ //
+ // 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);
+
+ AddTrackCuts(esdTrackCuts);
+ delete esdTrackCuts;
+ esdTrackCuts=NULL;
+
+
+ //
+ // Cut values per pt bin
+ //
+ const Int_t nvars=11;
+ const Int_t nptbins=14;
+ Float_t* ptbins;
+ ptbins=new Float_t[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]=7.;
+ ptbins[9]=8.;
+ ptbins[10]=12.;
+ ptbins[11]=16.;
+ ptbins[12]=20.;
+ ptbins[13]=24.;
+ ptbins[14]=9999.;
+
+ SetPtBins(nptbins+1,ptbins);
+
+ //setting cut values
+ Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,350.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,3.},/* pt<0.5*/
+ {0.400,350.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,3.},/* 0.5<pt<1*/
+ {0.400,300.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-30000.*1E-8,0.80,0.,4.},/* 1<pt<2 */
+ {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-20000.*1E-8,0.85,0.,4.},/* 2<pt<3 */
+ {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-15000.*1E-8,0.85,0.,4.},/* 3<pt<4 */
+ {0.400,300.*1E-4,0.75,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.875,0.,0.},/* 4<pt<5 */
+ {0.400,300.*1E-4,0.75,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.875,0.,0.},/* 5<pt<6 */
+ {0.400,400.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 6<pt<7 */
+ {0.400,400.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 7<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;
+
+ SetUseSpecialCuts(kTRUE);
+ SetMaximumPtSpecialCuts(8.);
+
+ //Do recalculate the vertex
+ SetRemoveDaughtersFromPrim(kTRUE);
+
+ //
+ // Pid settings
+ //
+ Bool_t pidflag=kTRUE;
+ SetUsePID(pidflag);
+ if(pidflag) cout<<"PID is used"<<endl;
+ else cout<<"PID is not used"<<endl;
+ //
+ 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->SetPCompatTOF(1.5);
+ pidObj->SetSigmaForTPCCompat(3.);
+ pidObj->SetSigmaForTOFCompat(3.);
+ pidObj->SetOldPid(kTRUE);
+ // pidObj->SetOldPid(kFALSE);
+ SetPidHF(pidObj);
+
+ SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
+
+ SetLowPt(kFALSE);
+
+ //activate pileup rejection (for pp)
+ SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
+
+
+ PrintAll();
+
+ delete pidObj;
+ pidObj=NULL;
+ delete [] ptbins;
+ ptbins=NULL;
+
+ return;
+}
+
+//---------------------------------------------------------------------------
void AliRDHFCutsD0toKpi::SetStandardCutsPP2011_276TeV() {
//
// STANDARD CUTS USED FOR 2011 pp analysis at 2.76TeV
esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
AddTrackCuts(esdTrackCuts);
+ delete esdTrackCuts;
+ esdTrackCuts=NULL;
const Int_t nvars=11;
pidObj->SetCompat(compat);
pidObj->SetTPC(kTRUE);
pidObj->SetTOF(kTRUE);
+ pidObj->SetOldPid(kTRUE);
SetPidHF(pidObj);
SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
}
+//---------------------------------------------------------------------------
void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
//
// Final CUTS USED FOR 2010 PbPb analysis of central events
AddTrackCuts(esdTrackCuts);
+ delete esdTrackCuts;
+ esdTrackCuts=NULL;
// SPD k FIRST for Pt<3 GeV/c
SetSelectCandTrackSPDFirst(kTRUE, 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,
2.},/* 1<pt<2 */
+ {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,-45000.*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 */
pidObj->SetPCompatTOF(2.);
pidObj->SetSigmaForTPCCompat(3.);
pidObj->SetSigmaForTOFCompat(3.);
+ pidObj->SetOldPid(kTRUE);
SetPidHF(pidObj);
}
+//---------------------------------------------------------------------------
void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010Peripherals() {
// CUTS USED FOR D0 RAA for the CENTRALITY RANGE 20-80%
AddTrackCuts(esdTrackCuts);
+ delete esdTrackCuts;
+ esdTrackCuts=NULL;
// SPD k FIRST for Pt<3 GeV/c
SetSelectCandTrackSPDFirst(kTRUE, 3);
pidObj->SetPCompatTOF(2.);
pidObj->SetSigmaForTPCCompat(3.);
pidObj->SetSigmaForTOFCompat(3.);
+ pidObj->SetOldPid(kTRUE);
SetPidHF(pidObj);
SetUsePID(kTRUE);
}
+//---------------------------------------------------------------------------
void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2011() {
// Default 2010 PbPb cut object
SetStandardCutsPbPb2010();
+ AliAODPidHF *pidobj=GetPidHF();
+
+ pidobj->SetOldPid(kFALSE);
//
// Enable all 2011 PbPb run triggers
EnableSemiCentralTrigger();
}
+
+//---------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::IsSelectedCombPID(AliAODRecoDecayHF* d)
+{
+ // ############################################################
+ //
+ // Apply Bayesian PID selection
+ // Implementation of Bayesian PID: Jeremy Wilkinson
+ //
+ // ############################################################
+
+ if (!fUsePID || !d) return 3;
+
+
+ if (fBayesianStrategy == kBayesWeightNoFilter) {
+ //WeightNoFilter: Accept all particles (no PID cut) but fill mass histos with weights in task
+ CalculateBayesianWeights(d);
+ return 3;
+ }
+
+
+
+ Int_t isD0 = 0;
+ Int_t isD0bar = 0;
+ Int_t returnvalue = 0;
+
+ Int_t isPosKaon = 0, isNegKaon = 0, isPosPion = 0, isNegPion = 0;
+
+ //Bayesian methods used here check for ID of kaon, and whether it is positive or negative.
+
+
+ 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
+ }
+ Double_t momentumpositive=0., momentumnegative=0.; //Used in "prob > prior" method
+ for (Int_t daught = 0; daught < 2; daught++) {
+ //Loop over prongs
+
+ if (aodtrack[daught]->Charge() == -1) {
+ momentumnegative = aodtrack[daught]->P();
+ }
+ if (aodtrack[daught]->Charge() == +1) {
+ momentumpositive = aodtrack[daught]->P();
+ }
+ if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
+ checkPIDInfo[daught] = kFALSE;
+ continue;
+ }
+
+ }
+
+ //Loop over daughters ends here
+
+ if (!checkPIDInfo[0] && !checkPIDInfo[1]) {
+ return 0; //Reject if both daughters lack both TPC and TOF info
+ }
+
+
+ CalculateBayesianWeights(d); //Calculates all Bayesian probabilities for both positive and negative tracks
+ // Double_t prob0[AliPID::kSPECIES];
+ // fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), prob0);
+ ///Three possible Bayesian probability cuts: Picked using SetBayesianCondition(int).
+ switch (fBayesianCondition) {
+ ///A: Standard max. probability method (accept most likely species)
+ case kMaxProb:
+ 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, fWeightsPositive) == fWeightsPositive[AliPID::kPion]) { //If highest probability lies with pion
+ isPosPion = 1; //flag [daught] as a pion
+ }
+
+ if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kKaon]) { //If highest probability lies with kaon
+ isNegKaon = 1; //flag [daught] as a kaon
+ }
+
+ if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kPion]) { //If highest probability lies with kaon
+ isNegPion = 1; //flag [daught] as a pion
+ }
+
+
+ break;
+ ///B: Accept if probability greater than prior
+ case kAbovePrior:
+
+ if (fWeightsNegative[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
+ GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumnegative)))) { //Retrieves relevant prior, gets value at momentum
+ isNegKaon = 1;
+ }
+ if (fWeightsPositive[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
+ GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumpositive)))) { //Retrieves relevant prior, gets value at momentum
+ isPosKaon = 1;
+ }
+
+ break;
+
+ ///C: Accept if probability greater than user-defined threshold
+ case kThreshold:
+ if (fWeightsNegative[AliPID::kKaon] > fProbThreshold) {
+ isNegKaon = 1;
+ }
+ if (fWeightsNegative[AliPID::kPion] > fProbThreshold) {
+ isNegPion = 1;
+ }
+
+ if (fWeightsPositive[AliPID::kKaon] > fProbThreshold) {
+ isPosKaon = 1;
+ }
+
+ if (fWeightsPositive[AliPID::kPion] > fProbThreshold) {
+ isPosPion = 1;
+ }
+
+ break;
+ }
+
+
+ //Momentum-based selection (also applied to filtered weighted method)
+
+ if (fBayesianStrategy == kBayesMomentum || fBayesianCondition == kBayesWeight) {
+ if (isNegKaon && isPosKaon) { // If both are kaons, reject
+ isD0 = 1;
+ isD0bar = 1;
+ } else if (isNegKaon && isPosPion) { //If negative kaon present, D0
+ isD0 = 1;
+ } else if (isPosKaon && isNegPion) { //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;
+ }
+ }
+
+ //Simple Bayesian
+ if (fBayesianStrategy == kBayesSimple) {
+
+ if (isPosKaon && isNegKaon) { //If both are ID'd as kaons, accept as possible
+ returnvalue = 3;
+ } else if (isNegKaon && isPosPion) { //If negative kaon, D0
+ returnvalue = 1;
+ } else if (isPosKaon && isNegPion) { //If positive kaon, D0-bar
+ returnvalue = 2;
+ } else if (isPosPion && isNegPion) {
+ returnvalue = 0; //If neither kaon, reject
+ } else {returnvalue = 0;} //default
+
+ }
+
+ return returnvalue;
+
+
+
+}
+
+
+
+//---------------------------------------------------------------------------
+void AliRDHFCutsD0toKpi::CalculateBayesianWeights(AliAODRecoDecayHF* d)
+
+{
+ //Function to compute weights for Bayesian method
+
+ 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, define weights
+ if (aodtrack[daught]->Charge() == +1) {
+ fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsPositive);
+ }
+
+ if (aodtrack[daught]->Charge() == -1) {
+ fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsNegative);
+ }
+ }
+}
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff