From e2aa82b63a0b3cb33e5f82eac76874943fd84554 Mon Sep 17 00:00:00 2001 From: fprino Date: Sun, 7 Apr 2013 20:58:52 +0000 Subject: [PATCH] Option to use Bayesian PID added for D0 (Jeremy) --- PWGHF/vertexingHF/AliAnalysisTaskSED0Mass.cxx | 200 +++++++-- PWGHF/vertexingHF/AliAnalysisTaskSED0Mass.h | 13 +- PWGHF/vertexingHF/AliRDHFCutsD0toKpi.cxx | 388 +++++++++++++++++- PWGHF/vertexingHF/AliRDHFCutsD0toKpi.h | 41 +- PWGHF/vertexingHF/macros/AddTaskD0Mass.C | 24 +- 5 files changed, 613 insertions(+), 53 deletions(-) diff --git a/PWGHF/vertexingHF/AliAnalysisTaskSED0Mass.cxx b/PWGHF/vertexingHF/AliAnalysisTaskSED0Mass.cxx index 41d0ad705bc..e64ac6a5c37 100644 --- a/PWGHF/vertexingHF/AliAnalysisTaskSED0Mass.cxx +++ b/PWGHF/vertexingHF/AliAnalysisTaskSED0Mass.cxx @@ -23,6 +23,7 @@ // Authors: A.Dainese, andrea.dainese@lnl.infn.it // Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass) // Carmelo Di Giglio, carmelo.digiglio@ba.infn.it (like sign) +// Jeremy Wilkinson, jwilkinson@physi.uni-heidelberg.de (weighted Bayesian) ///////////////////////////////////////////////////////////// #include @@ -85,7 +86,10 @@ AliAnalysisTaskSE(), fUseSelectionBit(kTRUE), fWriteVariableTree(kFALSE), fVariablesTree(0), - fCandidateVariables() + fCandidateVariables(), + fPIDCheck(kFALSE), + fDrawDetSignal(kFALSE), + fDetSignal(0) { // Default constructor for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0; @@ -118,7 +122,10 @@ AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0t fUseSelectionBit(kTRUE), fWriteVariableTree(kFALSE), fVariablesTree(0), - fCandidateVariables() + fCandidateVariables(), + fPIDCheck(kFALSE), + fDrawDetSignal(kFALSE), + fDetSignal(0) { // Default constructor @@ -141,6 +148,8 @@ AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0t DefineOutput(6,TList::Class()); //My private output // Output slot #7 keeps a tree of the candidate variables after track selection DefineOutput(7,TTree::Class()); //My private outpu + // Output slot #8 writes into a TList container (Detector signals) + DefineOutput(8, TList::Class()); //My private output } //________________________________________________________________________ @@ -177,6 +186,10 @@ AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass() delete fVariablesTree; fVariablesTree = 0; } + if (fDetSignal) { + delete fDetSignal; + fDetSignal = 0; + } } @@ -219,6 +232,10 @@ void AliAnalysisTaskSED0Mass::UserCreateOutputObjects() fDistr->SetOwner(); fDistr->SetName("distributionslist"); + fDetSignal = new TList(); + fDetSignal->SetOwner(); + fDetSignal->SetName("detectorsignals"); + TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" "; TString nameMassPt="", nameSgnPt="", nameBkgPt="", nameRflPt=""; Int_t nbins2dPt=60; Float_t binInPt=0., binFinPt=30.; @@ -708,7 +725,7 @@ void AliAnalysisTaskSED0Mass::UserCreateOutputObjects() fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1"); fNentries->GetXaxis()->SetBinLabel(6,"no daughter"); if(fSys==0) fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)"); - if(fFillVarHists){ + if(fFillVarHists || fPIDCheck){ fNentries->GetXaxis()->SetBinLabel(8,"PID=0"); fNentries->GetXaxis()->SetBinLabel(9,"PID=1"); fNentries->GetXaxis()->SetBinLabel(10,"PID=2"); @@ -736,16 +753,42 @@ void AliAnalysisTaskSED0Mass::UserCreateOutputObjects() Int_t nVar = 15; fCandidateVariables = new Double_t [nVar]; TString * fCandidateVariableNames = new TString[nVar]; - fCandidateVariableNames[0]="massD0"; fCandidateVariableNames[1]="massD0bar"; fCandidateVariableNames[2]="pt"; - fCandidateVariableNames[3]="dca"; fCandidateVariableNames[4]="costhsD0"; fCandidateVariableNames[5]="costhsD0bar"; - fCandidateVariableNames[6]="ptk"; fCandidateVariableNames[7]="ptpi"; - fCandidateVariableNames[8]="d0k"; fCandidateVariableNames[9]="d0pi"; fCandidateVariableNames[10]="d0xd0"; - fCandidateVariableNames[11]="costhp"; fCandidateVariableNames[12]="costhpxy"; - fCandidateVariableNames[13]="lxy"; fCandidateVariableNames[14]="specialcuts"; + fCandidateVariableNames[0] = "massD0"; + fCandidateVariableNames[1] = "massD0bar"; + fCandidateVariableNames[2] = "pt"; + fCandidateVariableNames[3] = "dca"; + fCandidateVariableNames[4] = "costhsD0"; + fCandidateVariableNames[5] = "costhsD0bar"; + fCandidateVariableNames[6] = "ptk"; + fCandidateVariableNames[7] = "ptpi"; + fCandidateVariableNames[8] = "d0k"; + fCandidateVariableNames[9] = "d0pi"; + fCandidateVariableNames[10] = "d0xd0"; + fCandidateVariableNames[11] = "costhp"; + fCandidateVariableNames[12] = "costhpxy"; + fCandidateVariableNames[13] = "lxy"; + fCandidateVariableNames[14] = "specialcuts"; for(Int_t ivar=0; ivarBranch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/d",fCandidateVariableNames[ivar].Data())); } + // + // Output slot 8 : List for detector response histograms + // + fDetSignal = new TList(); + if (fDrawDetSignal) { + TH2F *TOFSigBefPID = new TH2F("TOFSigBefPID", "TOF signal of daughters before PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3); + TH2F *TOFSigAftPID = new TH2F("TOFSigAftPID", "TOF signal after PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3); + + TH2F *TPCSigBefPID = new TH2F("TPCSigBefPID", "TPC dE/dx before PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500); + TH2F *TPCSigAftPID = new TH2F("TPCSigAftPID", "TPC dE/dx after PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500); + + fDetSignal->Add(TOFSigBefPID); + fDetSignal->Add(TOFSigAftPID); + fDetSignal->Add(TPCSigBefPID); + fDetSignal->Add(TPCSigAftPID); + } + // Post the data PostData(1,fOutputMass); PostData(2,fDistr); @@ -753,6 +796,7 @@ void AliAnalysisTaskSED0Mass::UserCreateOutputObjects() PostData(5,fCounter); PostData(6,fOutputMassPt); PostData(7,fVariablesTree); + PostData(8, fDetSignal); return; } @@ -934,7 +978,21 @@ void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/) FillVarHists(aod,d,mcArray,fCuts,fDistr); } + if (fDrawDetSignal) { + DrawDetSignal(d, fDetSignal); + } + FillMassHists(d,mcArray,mcHeader,fCuts,fOutputMass); + if (fPIDCheck) { + Int_t isSelectedPIDfill = 3; + if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(d); //0 rejected,1 D0,2 Dobar, 3 both + + if (isSelectedPIDfill == 0)fNentries->Fill(7); + if (isSelectedPIDfill == 1)fNentries->Fill(8); + if (isSelectedPIDfill == 2)fNentries->Fill(9); + if (isSelectedPIDfill == 3)fNentries->Fill(10); + } + } fDaughterTracks.Clear(); @@ -949,6 +1007,45 @@ void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/) PostData(5,fCounter); PostData(6,fOutputMassPt); PostData(7,fVariablesTree); + PostData(8, fDetSignal); + + return; +} + +//____________________________________________________________________________ +void AliAnalysisTaskSED0Mass::DrawDetSignal(AliAODRecoDecayHF2Prong *part, TList *ListDetSignal) +{ + // + // Function called in UserExec for drawing detector signal histograms: + // + fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(0), 0); + fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(1), 1); + + AliESDtrack *esdtrack1 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0)); + AliESDtrack *esdtrack2 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1)); + + + //For filling detector histograms + Int_t isSelectedPIDfill = 3; + if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both + + + //fill "before PID" histos with every daughter + ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal()); + ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal()); + ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal()); + ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal()); + + if (isSelectedPIDfill != 0) { //fill "After PID" with everything that's not rejected + ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal()); + ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal()); + ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal()); + ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal()); + + } + + delete esdtrack1; + delete esdtrack2; return; } @@ -966,11 +1063,13 @@ void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Pr //Double_t pt = d->Pt(); //mother pt Int_t isSelectedPID=3; if(!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPID=cuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both - if (isSelectedPID==0)fNentries->Fill(7); - if (isSelectedPID==1)fNentries->Fill(8); - if (isSelectedPID==2)fNentries->Fill(9); - if (isSelectedPID==3)fNentries->Fill(10); - //fNentries->Fill(8+isSelectedPID); + if (!fPIDCheck) { //if fPIDCheck, this is already filled elsewhere + if (isSelectedPID==0)fNentries->Fill(7); + if (isSelectedPID==1)fNentries->Fill(8); + if (isSelectedPID==2)fNentries->Fill(9); + if (isSelectedPID==3)fNentries->Fill(10); + //fNentries->Fill(8+isSelectedPID); + } if(fCutOnDistr && !fIsSelectedCandidate) return; //printf("\nif no cuts or cuts passed\n"); @@ -1575,6 +1674,15 @@ void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClon Bool_t isPrimary=kTRUE; if (fReadMC) labD0 = part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx) + //Define weights for Bayesian (if appropriate) + + Double_t weigD0=1.; + Double_t weigD0bar=1.; + if (fCuts->GetCombPID() && (fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeight || fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeightNoFilter)) { + weigD0=fCuts->GetWeightsNegative()[AliPID::kKaon] * fCuts->GetWeightsPositive()[AliPID::kPion]; + weigD0bar=fCuts->GetWeightsPositive()[AliPID::kKaon] * fCuts->GetWeightsNegative()[AliPID::kPion]; + } + //count candidates selected by cuts fNentries->Fill(1); //count true D0 selected by cuts @@ -1601,17 +1709,17 @@ void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClon // cout<<"Fill S with D0"<FindObject(fillthis)))->Fill(invmassD0); + ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0); if(fFillPtHist){ fillthismasspt="histSgnPt"; - ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt); + ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0); } if(fFillImpParHist){ - if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse); + if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0); else { - fHistMassPtImpParTC[2]->Fill(arrayForSparse); - fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue); + fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0); + fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0); } } @@ -1619,32 +1727,32 @@ void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClon if(TMath::Abs(invmassD0 - mPDG) < 0.027 && fFillVarHists){ fillthis="histSgn27_"; fillthis+=ptbin; - ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0); + ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0); } } } else{ //it was a D0bar fillthis="histRfl_"; fillthis+=ptbin; - ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0); + ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0); if(fFillPtHist){ fillthismasspt="histRflPt"; // cout << " Filling "<FindObject(fillthismasspt)))->Fill(invmassD0,pt); + ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0); } } } else {//background fillthis="histBkg_"; fillthis+=ptbin; - ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0); + ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0); if(fFillPtHist){ fillthismasspt="histBkgPt"; // cout << " Filling "<FindObject(fillthismasspt)))->Fill(invmassD0,pt); + ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0); } - if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse); + if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0); } @@ -1654,16 +1762,17 @@ void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClon // cout<<"Filling "<FindObject(fillthis)))->Fill(invmassD0); + ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0); + if(fFillPtHist){ fillthismasspt="histMassPt"; // cout<<"Filling "<FindObject(fillthismasspt)))->Fill(invmassD0,pt); + ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0); } if(fFillImpParHist) { // cout << "Filling fHistMassPtImpParTC[0]"<Fill(arrayForSparse); + fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0); } } @@ -1688,7 +1797,7 @@ void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClon if (pdgD0==-421){ //D0bar fillthis="histSgn_"; fillthis+=ptbin; - ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar); + ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar); // if (TMath::Abs(invmassD0bar - mPDG) < 0.027){ // fillthis="histSgn27_"; // fillthis+=ptbin; @@ -1698,52 +1807,54 @@ void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClon if(fFillPtHist){ fillthismasspt="histSgnPt"; // cout<<" Filling "<< fillthismasspt << endl; - ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt); + ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar); } if(fFillImpParHist){ // cout << " Filling impact parameter thnsparse"<Fill(arrayForSparse); + if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0bar); else { - fHistMassPtImpParTC[2]->Fill(arrayForSparse); - fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue); + fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0bar); + fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0bar); } } } else{ fillthis="histRfl_"; fillthis+=ptbin; - ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar); + ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar); if(fFillPtHist){ fillthismasspt="histRflPt"; // cout << " Filling "<FindObject(fillthismasspt)))->Fill(invmassD0bar,pt); + ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar); } } } else {//background or LS fillthis="histBkg_"; fillthis+=ptbin; - ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar); + ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar); if(fFillPtHist){ fillthismasspt="histBkgPt"; // cout<<" Filling "<< fillthismasspt << endl; - ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt); + ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar); } - if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse); + if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0bar); } }else{ fillthis="histMass_"; fillthis+=ptbin; // printf("Fill mass with D0bar"); - ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar); + + ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar,weigD0bar); + if(fFillPtHist){ fillthismasspt="histMassPt"; // cout<<" Filling "<< fillthismasspt << endl; - ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt); + ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar); } - if(fFillImpParHist) fHistMassPtImpParTC[0]->Fill(arrayForSparse); + if(fFillImpParHist) fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0bar); } } @@ -1844,6 +1955,13 @@ void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/) printf("ERROR: fCounter not available\n"); return; } + if (fDrawDetSignal) { + fDetSignal = dynamic_cast(GetOutputData(8)); + if (!fDetSignal) { + printf("ERROR: fDetSignal not available\n"); + return; + } + } Int_t nptbins=fCuts->GetNPtBins(); for(Int_t ipt=0;iptKpi + ClassDef(AliAnalysisTaskSED0Mass,18); // AliAnalysisTaskSE for D0->Kpi }; #endif diff --git a/PWGHF/vertexingHF/AliRDHFCutsD0toKpi.cxx b/PWGHF/vertexingHF/AliRDHFCutsD0toKpi.cxx index b04d4500291..340bdebe8cd 100644 --- a/PWGHF/vertexingHF/AliRDHFCutsD0toKpi.cxx +++ b/PWGHF/vertexingHF/AliRDHFCutsD0toKpi.cxx @@ -49,7 +49,11 @@ fDefaultPID(kFALSE), fUseKF(kFALSE), fPtLowPID(2.), fPtMaxSpecialCuts(9999.), -fmaxPtrackForPID(9999) +fmaxPtrackForPID(9999), +fCombPID(kFALSE), +fWeightsPositive(new Double_t[AliPID::kSPECIES]), +fWeightsNegative(new Double_t[AliPID::kSPECIES]), +fBayesianStrategy(0) { // // Default Constructor @@ -104,7 +108,11 @@ AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) : fUseKF(source.fUseKF), fPtLowPID(source.fPtLowPID), fPtMaxSpecialCuts(source.fPtMaxSpecialCuts), - fmaxPtrackForPID(source.fmaxPtrackForPID) + fmaxPtrackForPID(source.fmaxPtrackForPID), + fCombPID(source.fCombPID), + fWeightsPositive(source.fWeightsPositive), + fWeightsNegative(source.fWeightsNegative), + fBayesianStrategy(source.fBayesianStrategy) { // // Copy constructor @@ -127,10 +135,31 @@ AliRDHFCutsD0toKpi &AliRDHFCutsD0toKpi::operator=(const AliRDHFCutsD0toKpi &sour 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) { // @@ -371,11 +400,38 @@ Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEve 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; @@ -682,6 +738,7 @@ Int_t AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* d) 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; @@ -1147,6 +1204,7 @@ void AliRDHFCutsD0toKpi::SetUseKF(Bool_t useKF) } +//--------------------------------------------------------------------------- void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() { // //STANDARD CUTS USED FOR 2010 pp analysis @@ -1272,6 +1330,7 @@ void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() { } +//--------------------------------------------------------------------------- void AliRDHFCutsD0toKpi::SetStandardCutsPP2011_276TeV() { // // STANDARD CUTS USED FOR 2011 pp analysis at 2.76TeV @@ -1397,6 +1456,7 @@ void AliRDHFCutsD0toKpi::SetStandardCutsPP2011_276TeV() { } +//--------------------------------------------------------------------------- void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() { // // Final CUTS USED FOR 2010 PbPb analysis of central events @@ -1537,6 +1597,7 @@ void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() { } +//--------------------------------------------------------------------------- void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010Peripherals() { // CUTS USED FOR D0 RAA for the CENTRALITY RANGE 20-80% @@ -1673,6 +1734,7 @@ void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010Peripherals() { } +//--------------------------------------------------------------------------- void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2011() { // Default 2010 PbPb cut object @@ -1690,3 +1752,319 @@ void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2011() { 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); + } + + } +} + + + diff --git a/PWGHF/vertexingHF/AliRDHFCutsD0toKpi.h b/PWGHF/vertexingHF/AliRDHFCutsD0toKpi.h index 57d49614d63..149e2142e89 100644 --- a/PWGHF/vertexingHF/AliRDHFCutsD0toKpi.h +++ b/PWGHF/vertexingHF/AliRDHFCutsD0toKpi.h @@ -23,7 +23,7 @@ class AliRDHFCutsD0toKpi : public AliRDHFCuts AliRDHFCutsD0toKpi(const char* name="CutsD0toKpi"); - virtual ~AliRDHFCutsD0toKpi(){} + virtual ~AliRDHFCutsD0toKpi(); AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi& source); AliRDHFCutsD0toKpi& operator=(const AliRDHFCutsD0toKpi& source); @@ -39,6 +39,10 @@ class AliRDHFCutsD0toKpi : public AliRDHFCuts {return IsSelected(obj,selectionLevel,0);} virtual Int_t IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod); + virtual Int_t IsSelectedCombPID(AliAODRecoDecayHF* d); + virtual Int_t IsSelectedSimpleBayesianPID(AliAODRecoDecayHF* d); + virtual void CalculateBayesianWeights(AliAODRecoDecayHF* d); + Float_t GetMassCut(Int_t iPtBin=0) const { return (GetCuts() ? fCutsRD[GetGlobalIndex(0,iPtBin)] : 1.e6);} Float_t GetDCACut(Int_t iPtBin=0) const { return (GetCuts() ? fCutsRD[GetGlobalIndex(1,iPtBin)] : 1.e6);} Int_t CombineSelectionLevels(Int_t selectionvalTrack,Int_t selectionvalCand,Int_t selectionvalPID)const; @@ -63,6 +67,31 @@ class AliRDHFCutsD0toKpi : public AliRDHFCuts Double_t GetPtForPIDtight()const {return fPtLowPID;} void SetUseKF(Bool_t useKF); Bool_t GetIsUsedKF() const {return fUseKF;} + void SetWeightsPositive(Double_t* weights){ + for (Int_t i = 0; iKpi + Bool_t fCombPID; //switch for Bayesian + + Double_t* fWeightsPositive; //Bayesian weights for positive track + Double_t* fWeightsNegative; //Bayesian weights for negative track + Int_t fBayesianStrategy; + + ClassDef(AliRDHFCutsD0toKpi,9); // class for cuts on AOD reconstructed D0->Kpi }; #endif diff --git a/PWGHF/vertexingHF/macros/AddTaskD0Mass.C b/PWGHF/vertexingHF/macros/AddTaskD0Mass.C index a4d969f13c0..1d9d45d88bf 100644 --- a/PWGHF/vertexingHF/macros/AddTaskD0Mass.C +++ b/PWGHF/vertexingHF/macros/AddTaskD0Mass.C @@ -4,7 +4,8 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read Float_t minC=0, Float_t maxC=0, TString finDirname="Loose", TString finname="",TString finObjname="D0toKpiCuts", Bool_t flagAOD049=kFALSE, - Bool_t FillMassPt=false, Bool_t FillImpPar=false) + Bool_t FillMassPt=false, Bool_t FillImpPar=false, + Bool_t DrawDetSignal=false, Bool_t PIDCheck=false) { // // AddTask for the AliAnalysisTaskSE for D0 candidates @@ -22,7 +23,7 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read return NULL; } - TString filename="",out1name="",out2name="",out3name="",out4name="",out5name="",out6name="",out7name="",inname=""; + TString filename="",out1name="",out2name="",out3name="",out4name="",out5name="",out6name="",out7name="",out8name="", inname=""; filename = AliAnalysisManager::GetCommonFileName(); filename += ":PWG3_D2H_"; if(flag==0){ @@ -65,6 +66,16 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read out7name ="coutputVarTree"; + //detector signal hists + + out8name="detectorSignals"; + if(cutOnDistr) out8name+="C"; + if(flagD0D0bar==1)out8name+="D0"; + if(flagD0D0bar==2)out8name+="D0bar"; + + + + inname="cinputmassD0_0"; if(cutOnDistr) inname+="C"; if(flagD0D0bar==1)inname+="D0"; @@ -109,6 +120,7 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read out5name += finDirname.Data(); out6name += finDirname.Data(); out7name += finDirname.Data(); + out8name += finDirname.Data(); inname += finDirname.Data(); //setting my cut values @@ -172,6 +184,7 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read out5name+=centr; out6name+=centr; out7name+=centr; + out8name+=centr; inname+=centr; // Aanalysis task @@ -190,7 +203,8 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read massD0Task->SetFillPtHistos(FillMassPt); massD0Task->SetFillImpactParameterHistos(FillImpPar); - + massD0Task->SetDrawDetSignal(DrawDetSignal); + massD0Task->SetPIDCheck(PIDCheck); // massD0Task->SetRejectSDDClusters(kTRUE); // massD0Task->SetWriteVariableTree(kTRUE); @@ -209,6 +223,9 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read AliAnalysisDataContainer *coutputmassD05 = mgr->CreateContainer(out5name,AliNormalizationCounter::Class(),AliAnalysisManager::kOutputContainer, filename.Data()); //counter AliAnalysisDataContainer *coutputmassD06 = mgr->CreateContainer(out6name,TList::Class(),AliAnalysisManager::kOutputContainer, filename.Data()); //mass vs pt vs impt par AliAnalysisDataContainer *coutputmassD07 = mgr->CreateContainer(out7name,TTree::Class(),AliAnalysisManager::kOutputContainer, filename.Data()); //mass vs pt vs impt par + AliAnalysisDataContainer *coutputmassD08 = mgr->CreateContainer(out8name,TList::Class(),AliAnalysisManager::kOutputContainer, filename.Data()); //dedx + + mgr->ConnectInput(massD0Task,0,mgr->GetCommonInputContainer()); @@ -219,6 +236,7 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read mgr->ConnectOutput(massD0Task,5,coutputmassD05); mgr->ConnectOutput(massD0Task,6,coutputmassD06); mgr->ConnectOutput(massD0Task,7,coutputmassD07); + mgr->ConnectOutput(massD0Task,8,coutputmassD08); return massD0Task; -- 2.39.3