From: fprino Date: Wed, 30 Jul 2014 12:56:02 +0000 (+0200) Subject: Updated PID for D0 at low pt (Christian Mohler, Andrea Rossi) X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=f2221737b0c5a45a1f027b370d26dce420d48102;p=u%2Fmrichter%2FAliRoot.git Updated PID for D0 at low pt (Christian Mohler, Andrea Rossi) --- diff --git a/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.cxx b/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.cxx index 24cd4c92702..e556b880a22 100644 --- a/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.cxx +++ b/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.cxx @@ -46,7 +46,7 @@ ClassImp(AliAnalysisTaskCombinHF) //________________________________________________________________________ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(): AliAnalysisTaskSE(), - fOutput(0x0), + fOutput(0x0), fHistNEvents(0x0), fHistTrackStatus(0x0), fHistCheckOrigin(0x0), @@ -74,8 +74,9 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(): fTrackCutsKaon(0x0), fPidHF(new AliAODPidHF()), fAnalysisCuts(0x0), - fMinMass(1.750), + fMinMass(1.720), fMaxMass(2.150), + fMaxPt(10.), fEtaAccCut(0.9), fPtAccCut(0.1), fNRotations(9), @@ -89,9 +90,13 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(): fPromptFeeddown(kPrompt), fGoUpToQuark(kTRUE), fFullAnalysis(0), - fmaxPforIDPion(0.6), - fmaxPforIDKaon(2.), - fKeepNegID(kFALSE) + fPIDstrategy(knSigma), + fmaxPforIDPion(0.8), + fmaxPforIDKaon(2.), + fKeepNegID(kFALSE), + fPIDselCaseZero(0), + fBayesThresKaon(0.4), + fBayesThresPion(0.4) { // default constructor } @@ -99,7 +104,7 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(): //________________________________________________________________________ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analysiscuts): AliAnalysisTaskSE("DmesonCombin"), - fOutput(0x0), + fOutput(0x0), fHistNEvents(0x0), fHistTrackStatus(0x0), fHistCheckOrigin(0x0), @@ -127,8 +132,9 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analy fTrackCutsKaon(0x0), fPidHF(new AliAODPidHF()), fAnalysisCuts(analysiscuts), - fMinMass(1.750), + fMinMass(1.720), fMaxMass(2.150), + fMaxPt(10.), fEtaAccCut(0.9), fPtAccCut(0.1), fNRotations(9), @@ -142,16 +148,19 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analy fPromptFeeddown(1), fGoUpToQuark(kTRUE), fFullAnalysis(0), - fmaxPforIDPion(0.6), - fmaxPforIDKaon(2.), - fKeepNegID(kFALSE) - + fPIDstrategy(knSigma), + fmaxPforIDPion(0.8), + fmaxPforIDKaon(2.), + fKeepNegID(kFALSE), + fPIDselCaseZero(0), + fBayesThresKaon(0.4), + fBayesThresPion(0.4) { // standard constructor - + DefineOutput(1,TList::Class()); //My private output DefineOutput(2,AliNormalizationCounter::Class()); - } +} //________________________________________________________________________ AliAnalysisTaskCombinHF::~AliAnalysisTaskCombinHF() @@ -170,7 +179,7 @@ AliAnalysisTaskCombinHF::~AliAnalysisTaskCombinHF() delete fPtVsYGenLimAcc; delete fPtVsYGenAcc; delete fPtVsYReco; - delete fMassVsPtVsY; + delete fMassVsPtVsY; delete fMassVsPtVsYLSpp; delete fMassVsPtVsYLSmm; delete fMassVsPtVsYRot; @@ -194,11 +203,11 @@ void AliAnalysisTaskCombinHF::UserCreateOutputObjects() // Create the output container // if(fDebug > 1) printf("AnalysisTaskCombinHF::UserCreateOutputObjects() \n"); - + fOutput = new TList(); fOutput->SetOwner(); fOutput->SetName("OutputHistos"); - + fHistNEvents = new TH1F("hNEvents", "number of events ",8,-0.5,7.5); fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal"); fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected"); @@ -208,12 +217,12 @@ void AliAnalysisTaskCombinHF::UserCreateOutputObjects() fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for vertex out of accept"); fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for pileup events"); fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events"); - + fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE); fHistNEvents->Sumw2(); fHistNEvents->SetMinimum(0); fOutput->Add(fHistNEvents); - + fHistTrackStatus = new TH1F("hTrackStatus", "",8,-0.5,7.5); fHistTrackStatus->GetXaxis()->SetBinLabel(1,"Not OK"); fHistTrackStatus->GetXaxis()->SetBinLabel(2,"Track OK"); @@ -227,134 +236,136 @@ void AliAnalysisTaskCombinHF::UserCreateOutputObjects() fHistTrackStatus->Sumw2(); fHistTrackStatus->SetMinimum(0); fOutput->Add(fHistTrackStatus); - + + Int_t nPtBins = fMaxPt*10; + if(fReadMC){ - + fHistCheckOrigin=new TH1F("hCheckOrigin","",7,-1.5,5.5); fHistCheckOrigin->Sumw2(); fHistCheckOrigin->SetMinimum(0); fOutput->Add(fHistCheckOrigin); - + fHistCheckOriginSel=new TH1F("hCheckOriginSel","",7,-1.5,5.5); fHistCheckOriginSel->Sumw2(); fHistCheckOriginSel->SetMinimum(0); fOutput->Add(fHistCheckOriginSel); - + fHistCheckDecChan=new TH1F("hCheckDecChan","",7,-2.5,4.5); fHistCheckDecChan->Sumw2(); fHistCheckDecChan->SetMinimum(0); fOutput->Add(fHistCheckDecChan); - + fHistCheckDecChanAcc=new TH1F("hCheckDecChanAcc","",7,-2.5,4.5); fHistCheckDecChanAcc->Sumw2(); fHistCheckDecChanAcc->SetMinimum(0); fOutput->Add(fHistCheckDecChanAcc); - - fPtVsYGen= new TH2F("hPtVsYGen","",20,0.,10.,20,-1.,1.); + + fPtVsYGen= new TH2F("hPtVsYGen","",nPtBins,0.,fMaxPt,20,-1.,1.); fPtVsYGen->Sumw2(); fPtVsYGen->SetMinimum(0); fOutput->Add(fPtVsYGen); - - fPtVsYGenLimAcc= new TH2F("hPtVsYGenLimAcc","",20,0.,10.,20,-1.,1.); + + fPtVsYGenLimAcc= new TH2F("hPtVsYGenLimAcc","",nPtBins,0.,fMaxPt,20,-1.,1.); fPtVsYGenLimAcc->Sumw2(); fPtVsYGenLimAcc->SetMinimum(0); fOutput->Add(fPtVsYGenLimAcc); - - fPtVsYGenAcc= new TH2F("hPtVsYGenAcc","",20,0.,10.,20,-1.,1.); + + fPtVsYGenAcc= new TH2F("hPtVsYGenAcc","",nPtBins,0.,fMaxPt,20,-1.,1.); fPtVsYGenAcc->Sumw2(); fPtVsYGenAcc->SetMinimum(0); fOutput->Add(fPtVsYGenAcc); - - fPtVsYReco= new TH2F("hPtVsYReco","",20,0.,10.,20,-1.,1.); + + fPtVsYReco= new TH2F("hPtVsYReco","",nPtBins,0.,fMaxPt,20,-1.,1.); fPtVsYReco->Sumw2(); fPtVsYReco->SetMinimum(0); fOutput->Add(fPtVsYReco); } - - - Int_t nMassBins=fMaxMass*1000.-fMinMass*1000.; + + + Int_t nMassBins=fMaxMass*1000.-fMinMass*1000.; Double_t maxm=fMinMass+nMassBins*0.001; - fMassVsPtVsY=new TH3F("hMassVsPtVsY","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.); + fMassVsPtVsY=new TH3F("hMassVsPtVsY","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.); fMassVsPtVsY->Sumw2(); fMassVsPtVsY->SetMinimum(0); fOutput->Add(fMassVsPtVsY); - - fMassVsPtVsYRot=new TH3F("hMassVsPtVsYRot","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.); + + fMassVsPtVsYRot=new TH3F("hMassVsPtVsYRot","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.); fMassVsPtVsYRot->Sumw2(); fMassVsPtVsYRot->SetMinimum(0); fOutput->Add(fMassVsPtVsYRot); - + if(fMeson==kDzero){ - fMassVsPtVsYLSpp=new TH3F("hMassVsPtVsYLSpp","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.); + fMassVsPtVsYLSpp=new TH3F("hMassVsPtVsYLSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.); fMassVsPtVsYLSpp->Sumw2(); fMassVsPtVsYLSpp->SetMinimum(0); fOutput->Add(fMassVsPtVsYLSpp); - fMassVsPtVsYLSmm=new TH3F("hMassVsPtVsYLSmm","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.); + fMassVsPtVsYLSmm=new TH3F("hMassVsPtVsYLSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.); fMassVsPtVsYLSmm->Sumw2(); fMassVsPtVsYLSmm->SetMinimum(0); fOutput->Add(fMassVsPtVsYLSmm); } - - fMassVsPtVsYSig=new TH3F("hMassVsPtVsYSig","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.); + + fMassVsPtVsYSig=new TH3F("hMassVsPtVsYSig","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.); fMassVsPtVsYSig->Sumw2(); fMassVsPtVsYSig->SetMinimum(0); fOutput->Add(fMassVsPtVsYSig); - - fMassVsPtVsYRefl=new TH3F("hMassVsPtVsYRefl","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.); + + fMassVsPtVsYRefl=new TH3F("hMassVsPtVsYRefl","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.); fMassVsPtVsYRefl->Sumw2(); fMassVsPtVsYRefl->SetMinimum(0); fOutput->Add(fMassVsPtVsYRefl); - - fMassVsPtVsYBkg=new TH3F("hMassVsPtVsYBkg","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.); + + fMassVsPtVsYBkg=new TH3F("hMassVsPtVsYBkg","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.); fMassVsPtVsYBkg->Sumw2(); fMassVsPtVsYBkg->SetMinimum(0); fOutput->Add(fMassVsPtVsYBkg); - + fNSelected=new TH1F("hNSelected","",100,-0.5,99.5); fNSelected->Sumw2(); fNSelected->SetMinimum(0); fOutput->Add(fNSelected); - + fNormRotated=new TH1F("hNormRotated","",11,-0.5,10.5); fNormRotated->Sumw2(); fNormRotated->SetMinimum(0); fOutput->Add(fNormRotated); - + fDeltaMass=new TH1F("hDeltaMass","",100,-0.4,0.4); fDeltaMass->Sumw2(); fDeltaMass->SetMinimum(0); fOutput->Add(fDeltaMass); - + Int_t binSparseDMassRot[5]={nMassBins,100,24,40,20}; Double_t edgeLowSparseDMassRot[5]={fMinMass,-0.4,0.,-4.,0}; Double_t edgeHighSparseDMassRot[5]={maxm,0.4,12.,4.,3.14}; fDeltaMassFullAnalysis=new THnSparseF("fDeltaMassFullAnalysis","fDeltaMassFullAnalysis;inv mass (GeV/c);#Delta inv mass (GeV/c) ; p_{T}^{D} (GeV/c); #Delta p_{T} (GeV/c); daughter angle (2prongs) (rad);",5,binSparseDMassRot,edgeLowSparseDMassRot,edgeHighSparseDMassRot); fOutput->Add(fDeltaMassFullAnalysis); - + //Counter for Normalization fCounter = new AliNormalizationCounter("NormalizationCounter"); fCounter->Init(); - - PostData(1,fOutput); - PostData(2,fCounter); + + PostData(1,fOutput); + PostData(2,fCounter); } //________________________________________________________________________ void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){ //Build the 3-track combinatorics (+-+ and -+-) for D+->Kpipi decays - + AliAODEvent *aod = dynamic_cast (InputEvent()); if(!aod && AODEvent() && IsStandardAOD()) { - // In case there is an AOD handler writing a standard AOD, use the AOD - // event in memory rather than the input (ESD) event. + // In case there is an AOD handler writing a standard AOD, use the AOD + // event in memory rather than the input (ESD) event. aod = dynamic_cast (AODEvent()); } if(!aod){ printf("AliAnalysisTaskCombinHF::UserExec: AOD not found!\n"); return; } - - // fix for temporary bug in ESDfilter + + // fix for temporary bug in ESDfilter // the AODs with null vertex pointer didn't pass the PhysSel if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return; AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); @@ -362,23 +373,23 @@ void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){ AliPIDResponse *pidResp=inputHandler->GetPIDResponse(); fPidHF->SetPidResponse(pidResp); - + fHistNEvents->Fill(0); // count event // Post the data already here PostData(1,fOutput); - + fCounter->StoreEvent(aod,fAnalysisCuts,fReadMC); - + Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod); if(fAnalysisCuts->IsEventRejectedDueToTrigger())fHistNEvents->Fill(2); if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex())fHistNEvents->Fill(3); if(fAnalysisCuts->IsEventRejectedDueToVertexContributors())fHistNEvents->Fill(4); if(fAnalysisCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion())fHistNEvents->Fill(5); // if(fAnalysisCuts->IsEventRejectedDueToPileup())fHistNEvents->Fill(6); - if(fAnalysisCuts->IsEventRejectedDueToCentrality()) fHistNEvents->Fill(7); - + if(fAnalysisCuts->IsEventRejectedDueToCentrality()) fHistNEvents->Fill(7); + if(!isEvSel)return; - + fHistNEvents->Fill(1); TClonesArray *arrayMC=0; @@ -398,21 +409,41 @@ void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){ } FillGenHistos(arrayMC); } - - + + Int_t ntracks=aod->GetNTracks(); - + // select and flag tracks UChar_t* status = new UChar_t[ntracks]; for(Int_t iTr=0; iTrGetTrack(iTr); if(IsTrackSelected(track)) status[iTr]+=1; - if(IsKaon(track)) status[iTr]+=2; - if(IsPion(track)) status[iTr]+=4; + + // PID + if (fPIDstrategy == knSigma) { + // nsigma PID + if(IsKaon(track)) status[iTr]+=2; + if(IsPion(track)) status[iTr]+=4; + } + else if (fPIDstrategy == kBayesianMaxProb || fPIDstrategy == kBayesianThres) { + // Bayesian PID + Double_t *weights = new Double_t[AliPID::kSPECIES]; + fPidHF->GetPidCombined()->ComputeProbabilities(track, fPidHF->GetPidResponse(), weights); + if (fPIDstrategy == kBayesianMaxProb) { + if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kKaon]) status[iTr] += 2; + if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kPion]) status[iTr] += 4; + } + if (fPIDstrategy == kBayesianThres) { + if (weights[AliPID::kKaon] > fBayesThresKaon) status[iTr] += 2; + if (weights[AliPID::kPion] > fBayesThresPion) status[iTr] += 4; + } + delete[] weights; + } + fHistTrackStatus->Fill(status[iTr]); } - + // build the combinatorics Int_t nSelected=0; Int_t nFiltered=0; @@ -430,16 +461,16 @@ void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){ Double_t tmpp[3]; Double_t px[3],py[3],pz[3]; Int_t dgLabels[3]; - + for(Int_t iTr1=0; iTr1GetTrack(iTr1); if((status[iTr1] & 1)==0) continue; if((status[iTr1] & 2)==0) continue; Int_t chargeK=trK->Charge(); trK->GetPxPyPz(tmpp); - px[0] = tmpp[0]; - py[0] = tmpp[1]; - pz[0] = tmpp[2]; + px[0] = tmpp[0]; + py[0] = tmpp[1]; + pz[0] = tmpp[2]; dgLabels[0]=trK->GetLabel(); for(Int_t iTr2=0; iTr2GetTrack(iTr2); Int_t chargePi1=trPi1->Charge(); trPi1->GetPxPyPz(tmpp); - px[1] = tmpp[0]; - py[1] = tmpp[1]; - pz[1] = tmpp[2]; + px[1] = tmpp[0]; + py[1] = tmpp[1]; + pz[1] = tmpp[2]; dgLabels[1]=trPi1->GetLabel(); if(chargePi1==chargeK){ - if(fMeson==kDzero) FillLSHistos(421,2,tmpRD2,px,py,pz,pdg0,chargePi1); - continue; + if(fMeson==kDzero) FillLSHistos(421,2,tmpRD2,px,py,pz,pdg0,chargePi1); + continue; } if(fMeson==kDzero){ - nFiltered++; - v2->AddDaughter(trK); - v2->AddDaughter(trPi1); - tmpRD2->SetSecondaryVtx(v2); - Bool_t ok=FillHistos(421,2,tmpRD2,px,py,pz,pdg0,arrayMC,dgLabels); - v2->RemoveDaughters(); - if(ok) nSelected++; + nFiltered++; + v2->AddDaughter(trK); + v2->AddDaughter(trPi1); + tmpRD2->SetSecondaryVtx(v2); + Bool_t ok=FillHistos(421,2,tmpRD2,px,py,pz,pdg0,arrayMC,dgLabels); + v2->RemoveDaughters(); + if(ok) nSelected++; }else{ - for(Int_t iTr3=iTr2+1; iTr3GetTrack(iTr3); - Int_t chargePi2=trPi2->Charge(); - if(chargePi2==chargeK) continue; - nFiltered++; - trPi2->GetPxPyPz(tmpp); - px[2] = tmpp[0]; - py[2] = tmpp[1]; - pz[2] = tmpp[2]; - dgLabels[2]=trPi2->GetLabel(); - v3->AddDaughter(trK); - v3->AddDaughter(trPi1); - v3->AddDaughter(trPi2); - tmpRD3->SetSecondaryVtx(v3); - Bool_t ok=FillHistos(411,3,tmpRD3,px,py,pz,pdgp,arrayMC,dgLabels); - v3->RemoveDaughters(); - if(ok) nSelected++; - } + for(Int_t iTr3=iTr2+1; iTr3GetTrack(iTr3); + Int_t chargePi2=trPi2->Charge(); + if(chargePi2==chargeK) continue; + nFiltered++; + trPi2->GetPxPyPz(tmpp); + px[2] = tmpp[0]; + py[2] = tmpp[1]; + pz[2] = tmpp[2]; + dgLabels[2]=trPi2->GetLabel(); + v3->AddDaughter(trK); + v3->AddDaughter(trPi1); + v3->AddDaughter(trPi2); + tmpRD3->SetSecondaryVtx(v3); + Bool_t ok=FillHistos(411,3,tmpRD3,px,py,pz,pdgp,arrayMC,dgLabels); + v3->RemoveDaughters(); + if(ok) nSelected++; + } } } } - + delete [] status; delete v2; delete v3; delete tmpRD2; delete tmpRD3; - + fNSelected->Fill(nSelected); - + fCounter->StoreCandidates(aod,nFiltered,kTRUE); fCounter->StoreCandidates(aod,nSelected,kFALSE); - - PostData(1,fOutput); - PostData(2,fCounter); - + + PostData(1,fOutput); + PostData(2,fCounter); + return; } //________________________________________________________________________ void AliAnalysisTaskCombinHF::FillLSHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, Int_t charge){ // Fill histos for LS candidates - + tmpRD->SetPxPyPzProngs(nProngs,px,py,pz); Double_t pt = tmpRD->Pt(); Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau); @@ -547,23 +578,23 @@ void AliAnalysisTaskCombinHF::FillGenHistos(TClonesArray* arrayMC){ Bool_t isGoodDecay=kFALSE; Int_t labDau[4]={-1,-1,-1,-1}; if(fMeson==kDzero){ - deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,part,labDau); - if(part->GetNDaughters()!=2) continue; - if(deca==1) isGoodDecay=kTRUE; - }else if(fMeson==kDplus){ - deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,part,labDau); - if(deca>0) isGoodDecay=kTRUE; + deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,part,labDau); + if(part->GetNDaughters()!=2) continue; + if(deca==1) isGoodDecay=kTRUE; + }else if(fMeson==kDplus){ + deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,part,labDau); + if(deca>0) isGoodDecay=kTRUE; } fHistCheckDecChan->Fill(deca); if(labDau[0]==-1){ - // printf(Form("Meson %d Label of daughters not filled correctly -- %d\n",fMeson,isGoodDecay)); - continue; //protection against unfilled array of labels + // printf(Form("Meson %d Label of daughters not filled correctly -- %d\n",fMeson,isGoodDecay)); + continue; //protection against unfilled array of labels } Bool_t isInAcc=CheckAcceptance(arrayMC,nProng,labDau); if(isInAcc) fHistCheckDecChanAcc->Fill(deca); if(isGoodDecay){ - Double_t ptgen=part->Pt(); - Double_t ygen=part->Y(); + Double_t ptgen=part->Pt(); + Double_t ygen=part->Y(); if(fAnalysisCuts->IsInFiducialAcceptance(ptgen,ygen)){ fPtVsYGen->Fill(ptgen,ygen); if(TMath::Abs(ygen)<0.5) fPtVsYGenLimAcc->Fill(ptgen,ygen); @@ -577,47 +608,46 @@ void AliAnalysisTaskCombinHF::FillGenHistos(TClonesArray* arrayMC){ //________________________________________________________________________ Bool_t AliAnalysisTaskCombinHF::FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, TClonesArray *arrayMC, Int_t* dgLabels){ // Fill histos for candidates with proper charge sign - + Bool_t accept=kFALSE; - + tmpRD->SetPxPyPzProngs(nProngs,px,py,pz); Double_t pt = tmpRD->Pt(); Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau); Double_t mass=TMath::Sqrt(minv2); - + if(minv2>fMinMass*fMinMass && minv2Y(pdgD); if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){ fMassVsPtVsY->Fill(mass,pt,rapid); accept=kTRUE; if(fReadMC){ - Int_t signPdg[3]={0,0,0}; - for(Int_t iii=0; iiiMatchToMC(pdgD,arrayMC,nProngs,signPdg); - if(labD>=0){ - AliAODMCParticle* part = dynamic_cast(arrayMC->At(TMath::Abs(dgLabels[0]))); - if(part){ - Int_t pdgCode = TMath::Abs( part->GetPdgCode() ); - if(pdgCode==321){ - AliAODMCParticle* dmes = dynamic_cast(arrayMC->At(labD)); - if(dmes){ - Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,dmes,fGoUpToQuark); - if((fPromptFeeddown==kFeeddown && orig==5)|| (fPromptFeeddown==kPrompt && orig==4) || (fPromptFeeddown==kBoth && orig>=4)) { - fPtVsYReco->Fill(dmes->Pt(),dmes->Y()); - } - } - fMassVsPtVsYSig->Fill(mass,pt,rapid); - }else{ - fMassVsPtVsYRefl->Fill(mass,pt,rapid); - } - } - }else{ - fMassVsPtVsYBkg->Fill(mass,pt,rapid); - } + Int_t signPdg[3]={0,0,0}; + for(Int_t iii=0; iiiMatchToMC(pdgD,arrayMC,nProngs,signPdg); + if(labD>=0){ + AliAODMCParticle* part = dynamic_cast(arrayMC->At(TMath::Abs(dgLabels[0]))); + if(part){ + Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,part,fGoUpToQuark); + if((fPromptFeeddown==kFeeddown && orig==5)|| (fPromptFeeddown==kPrompt && orig==4) || (fPromptFeeddown==kBoth && orig>=4)) { + + Int_t pdgCode = TMath::Abs( part->GetPdgCode() ); + if(pdgCode==321){ + fMassVsPtVsYSig->Fill(mass,pt,rapid); + AliAODMCParticle* dmes = dynamic_cast(arrayMC->At(labD)); + fPtVsYReco->Fill(dmes->Pt(),dmes->Y()); + }else{ + fMassVsPtVsYRefl->Fill(mass,pt,rapid); + } + } + } + }else{ + fMassVsPtVsYBkg->Fill(mass,pt,rapid); + } } } } - + Int_t nRotated=0; Double_t massRot=0;// calculated later only if candidate is acceptable Double_t angleProngXY; @@ -626,11 +656,11 @@ Bool_t AliAnalysisTaskCombinHF::FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoD angleProngXY=TMath::ACos(((px[0]+px[1])*px[2]+(py[0]+py[1])*py[2])/TMath::Sqrt(((px[0]+px[1])*(px[0]+px[1])+(py[0]+py[1])*(py[0]+py[1]))*(px[2]*px[2]+py[2]*py[2]))); } Double_t ptOrig=pt; - - + + Double_t rotStep=(fMaxAngleForRot-fMinAngleForRot)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot Double_t rotStep3=(fMaxAngleForRot3-fMinAngleForRot3)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot - + for(Int_t irot=0; irotfMinMass*fMinMass && minv2Y(pdgD); if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){ - massRot=TMath::Sqrt(minv2); - fMassVsPtVsYRot->Fill(massRot,pt,rapid); - nRotated++; - fDeltaMass->Fill(massRot-mass); - if(fFullAnalysis){ - Double_t pointRot[5]={mass,massRot-mass,ptOrig,pt-ptOrig,angleProngXY}; - fDeltaMassFullAnalysis->Fill(pointRot); - } + massRot=TMath::Sqrt(minv2); + fMassVsPtVsYRot->Fill(massRot,pt,rapid); + nRotated++; + fDeltaMass->Fill(massRot-mass); + if(fFullAnalysis){ + Double_t pointRot[5]={mass,massRot-mass,ptOrig,pt-ptOrig,angleProngXY}; + fDeltaMassFullAnalysis->Fill(pointRot); + } } } px[0]=tmpx; @@ -668,59 +698,119 @@ Bool_t AliAnalysisTaskCombinHF::FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoD } } fNormRotated->Fill(nRotated); - + return accept; - + } //________________________________________________________________________ Bool_t AliAnalysisTaskCombinHF::IsTrackSelected(AliAODTrack* track){ // track selection cuts - + if(track->Charge()==0) return kFALSE; if(track->GetID()<0&&!fKeepNegID)return kFALSE; if(!(track->TestFilterMask(fFilterMask))) return kFALSE; - if(!SelectAODTrack(track,fTrackCutsAll)) return kFALSE; + if(!SelectAODTrack(track,fTrackCutsAll)) return kFALSE; return kTRUE; } + //________________________________________________________________________ Bool_t AliAnalysisTaskCombinHF::IsKaon(AliAODTrack* track){ // kaon selection cuts - + if(!fPidHF) return kTRUE; - Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon); - + Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon); + Double_t mom=track->P(); if(SelectAODTrack(track,fTrackCutsKaon)) { if(isKaon>=1) return kTRUE; - else if(isKaon>=0 && track->P()>fmaxPforIDKaon)return kTRUE; + if(isKaon<=-1) return kFALSE; + switch(fPIDselCaseZero){// isKaon=0 + case 0: + { + return kTRUE;// always accept + } + break; + case 1: + { + if(isKaon>=0 && track->P()>fmaxPforIDKaon)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDKaon + } + break; + case 2: + { + if(track->P()>fmaxPforIDKaon)return kTRUE; + AliPIDResponse *pidResp=fPidHF->GetPidResponse();// more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment + Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kKaon); + if(nsigma>-2.&& nsigma<3. && mom<0.6)isKaon=1; + else if(nsigma>-1.&& nsigma<3.&& mom<0.8)isKaon=1; + if(isKaon==1)return kTRUE; + } + break; + default: + { + AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero)); + return kFALSE;// actually case 0 could be set as the default and return kTRUE + } + } } + return kFALSE; } -//________________________________________________________________________ +//_______________________________________________________________________ Bool_t AliAnalysisTaskCombinHF::IsPion(AliAODTrack* track){ // pion selection cuts - + if(!fPidHF) return kTRUE; Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion); + Double_t mom=track->P(); if(SelectAODTrack(track,fTrackCutsPion)) { if(isPion>=1) return kTRUE; - else if(isPion>=0 && track->P()>fmaxPforIDPion)return kTRUE; + if(isPion<=-1) return kFALSE; + switch(fPIDselCaseZero){// isPion=0 + case 0: + { + return kTRUE;// always accept + } + break; + case 1: + { + if(track->P()>fmaxPforIDPion)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDPion + } + break; + case 2: + { + // more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment + if(track->P()>fmaxPforIDPion)return kTRUE; + AliPIDResponse *pidResp=fPidHF->GetPidResponse(); + Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kPion); + if(nsigma<2.&& nsigma>-3. && mom<0.6)isPion=1; + else if(nsigma<1. && nsigma> -3. && mom<0.8)isPion=1; + if(isPion==1)return kTRUE; + } + break; + default: + { + AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero)); + return kFALSE;// actually case 0 could be set as the default and return kTRUE + } + } } + return kFALSE; } + //________________________________________________________________________ Bool_t AliAnalysisTaskCombinHF::SelectAODTrack(AliAODTrack *track, AliESDtrackCuts *cuts){ // AOD track selection - + if(!cuts) return kTRUE; - + AliESDtrack esdTrack(track); // set the TPC cluster info esdTrack.SetTPCClusterMap(track->GetTPCClusterMap()); esdTrack.SetTPCSharedMap(track->GetTPCSharedMap()); esdTrack.SetTPCPointsF(track->GetTPCNclsF()); - if(!cuts->IsSelected(&esdTrack)) return kFALSE; - - return kTRUE; + if(!cuts->IsSelected(&esdTrack)) return kFALSE; + + return kTRUE; } //_________________________________________________________________ @@ -743,7 +833,7 @@ void AliAnalysisTaskCombinHF::Terminate(Option_t */*option*/) // if(fDebug > 1) printf("AliAnalysisTaskCombinHF: Terminate() \n"); fOutput = dynamic_cast (GetOutputData(1)); - if (!fOutput) { + if (!fOutput) { printf("ERROR: fOutput not available\n"); return; } diff --git a/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.h b/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.h index 2aa1b3f29d0..121b25b4538 100644 --- a/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.h +++ b/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.h @@ -4,7 +4,7 @@ /* Copyright(c) 1998-2018, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ -/* $Id: $ */ +/* $Id: $ */ //************************************************************************* // Class AliAnalysisTaskCombinHF @@ -23,18 +23,18 @@ class AliAnalysisTaskCombinHF : public AliAnalysisTaskSE { - public: - +public: + AliAnalysisTaskCombinHF(); AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analysiscuts); virtual ~AliAnalysisTaskCombinHF(); - + virtual void UserCreateOutputObjects(); virtual void Init(){}; virtual void LocalInit() {Init();} virtual void UserExec(Option_t *option); virtual void Terminate(Option_t *option); - + void SetReadMC(Bool_t read){fReadMC=read;} void SelectPromptD(){fPromptFeeddown=kPrompt;} void SelectFeeddownD(){fPromptFeeddown=kFeeddown;} @@ -60,31 +60,46 @@ class AliAnalysisTaskCombinHF : public AliAnalysisTaskSE void SetRDHFCuts(AliRDHFCuts* cuts){ fAnalysisCuts=cuts; } - void SetFilterMask(UInt_t mask=16){fFilterMask=mask;} + void SetFilterMask(UInt_t mask=16){fFilterMask=mask;} void SetAnalysisLevel(Int_t level){fFullAnalysis=level;} void ConfigureRotation(Int_t n, Double_t phimin, Double_t phimax){ fNRotations=n; fMinAngleForRot=phimin; fMaxAngleForRot=phimax; } + void SetMassWindow(Double_t minMass, Double_t maxMass){fMinMass=minMass; fMaxMass=maxMass;} + void SetMaxPt(Double_t maxPt){fMaxPt=maxPt;} + void SetEtaAccCut(Double_t etacut){fEtaAccCut=etacut;} + void SetPtAccCut(Double_t ptcut){fPtAccCut=ptcut;} + + void SetPIDstrategy(Int_t strat){fPIDstrategy=strat;} + void SetMaxPforIDPion(Double_t maxpIdPion){fmaxPforIDPion=maxpIdPion;} + void SetMaxPforIDKaon(Double_t maxpIdKaon){fmaxPforIDKaon=maxpIdKaon;} + void SetPIDselCaseZero(Int_t strat){fPIDselCaseZero=strat;} + void SetBayesThres(Double_t thresKaon, Double_t thresPion){ + fBayesThresKaon=thresKaon; + fBayesThresPion=thresPion; + } + Bool_t IsTrackSelected(AliAODTrack* track); Bool_t IsKaon(AliAODTrack* track); Bool_t IsPion(AliAODTrack* track); Bool_t SelectAODTrack(AliAODTrack *track, AliESDtrackCuts *cuts); + Bool_t FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, TClonesArray *arrayMC, Int_t* dgLabels); void FillLSHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, Int_t charge); void FillGenHistos(TClonesArray* arrayMC); Bool_t CheckAcceptance(TClonesArray* arrayMC, Int_t nProng, Int_t *labDau); + enum EMesonSpecies {kDzero, kDplus, kDstar, kDs}; enum EPrompFd {kNone,kPrompt,kFeeddown,kBoth}; - void SetMaxPforIDPion(Double_t maxpIdPion){fmaxPforIDPion=maxpIdPion;} - void SetMaxPforIDKaon(Double_t maxpIdKaon){fmaxPforIDKaon=maxpIdKaon;} - - private: - + enum EPIDstrategy {knSigma, kBayesianMaxProb, kBayesianThres}; + +private: + AliAnalysisTaskCombinHF(const AliAnalysisTaskCombinHF &source); - AliAnalysisTaskCombinHF& operator=(const AliAnalysisTaskCombinHF& source); - + AliAnalysisTaskCombinHF& operator=(const AliAnalysisTaskCombinHF& source); + TList *fOutput; //! list send on output slot 0 TH1F *fHistNEvents; //!hist. for No. of events TH1F *fHistTrackStatus; //!hist. of status of tracks @@ -107,37 +122,42 @@ class AliAnalysisTaskCombinHF : public AliAnalysisTaskSE TH1F *fNormRotated; //! hist. rotated/selected D+ TH1F *fDeltaMass; //! hist. mass difference after rotations THnSparse *fDeltaMassFullAnalysis; //! hist. mass difference after rotations with more details - + UInt_t fFilterMask; // FilterMask AliESDtrackCuts* fTrackCutsAll; // track selection AliESDtrackCuts* fTrackCutsPion; // pion track selection AliESDtrackCuts* fTrackCutsKaon; // kaon track selection AliAODPidHF* fPidHF; // PID configuration AliRDHFCuts *fAnalysisCuts; // Cuts for candidates + Double_t fMinMass; // minimum value of invariant mass Double_t fMaxMass; // maximum value of invariant mass - + Double_t fMaxPt; // maximum pT value for inv. mass histograms Double_t fEtaAccCut; // eta limits for acceptance step Double_t fPtAccCut; // pt limits for acceptance step - - + Int_t fNRotations; // number of rotations Double_t fMinAngleForRot; // minimum angle for track rotation Double_t fMaxAngleForRot; // maximum angle for track rotation Double_t fMinAngleForRot3; // minimum angle for track rotation (3rd prong) Double_t fMaxAngleForRot3; // maximum angle for track rotation (3rd prong) - + AliNormalizationCounter *fCounter;//!Counter for normalization - + Int_t fMeson; // mesonSpecies (see enum) Bool_t fReadMC; // flag for access to MC - Int_t fPromptFeeddown; // flag to select prompt (1), feeddown (2) or all (3) + Int_t fPromptFeeddown; // flag to select prompt (1), feeddown (2) or all (3) Bool_t fGoUpToQuark; // flag for definition of c,b origin Int_t fFullAnalysis; // flag to set analysis level (0 is the fastest) - + + Int_t fPIDstrategy; // knSigma, kBayesianMaxProb, kBayesianThres Double_t fmaxPforIDPion; // flag for upper p limit for id band for pion Double_t fmaxPforIDKaon; // flag for upper p limit for id band for kaon Bool_t fKeepNegID; // flag to keep also track with negative ID (default kFALSE, change it only if you know what you are doing) + Int_t fPIDselCaseZero; // flag to change PID strategy + Double_t fBayesThresKaon; // threshold for kaon identification via Bayesian PID + Double_t fBayesThresPion; // threshold for pion identification via Bayesian PID + ClassDef(AliAnalysisTaskCombinHF,3); // D+ task from AOD tracks }; diff --git a/PWGHF/vertexingHF/macros/AddTaskCombinHF.C b/PWGHF/vertexingHF/macros/AddTaskCombinHF.C index 38ec6707349..9cec59e23ac 100644 --- a/PWGHF/vertexingHF/macros/AddTaskCombinHF.C +++ b/PWGHF/vertexingHF/macros/AddTaskCombinHF.C @@ -1,4 +1,18 @@ -AliAnalysisTaskCombinHF *AddTaskCombinHF(Int_t meson=0, TString containerStr="",Bool_t readMC=kTRUE, TString cutObjFile="",TString cutObjNam="") +AliAnalysisTaskCombinHF *AddTaskCombinHF(Int_t meson = 0, + Bool_t readMC = kTRUE, + TString containerStr = "", + TString cutObjFile = "", + TString cutObjNam = "", + Int_t filterMask = 1, + Double_t ptcut = 0.3, + Double_t etacut = 0.9, + Int_t pidStrategy = 0, + Int_t casePID = 2, + Double_t bayesThresKaon = 0.4, + Double_t bayesThresPion = 0.4, + Double_t minMass = 1.6, + Double_t maxMass = 2.15, + Double_t maxPt = 20.) { AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); @@ -46,16 +60,30 @@ AliAnalysisTaskCombinHF *AddTaskCombinHF(Int_t meson=0, TString containerStr="", } dTask->SetReadMC(readMC); dTask->SetDebugLevel(0); + + dTask->SetFilterMask(filterMask); + dTask->SetPtAccCut(ptcut); + dTask->SetEtaAccCut(etacut); + + // mass and pt range for histograms + dTask->SetMassWindow(minMass, maxMass); + dTask->SetMaxPt(maxPt); + + // PID settings dTask->SetPIDHF(pid); - + dTask->SetPIDstrategy(pidStrategy); + dTask->SetPIDselCaseZero(casePID); + dTask->SetBayesThres(bayesThresKaon, bayesThresPion); + mgr->AddTask(dTask); + // Create containers for input/output TString mesname="Dzero"; if(meson==1) mesname="Dplus"; - TString inname = Form("cinput%s",mesname.Data()); + TString inname = Form("cinput%s%s",mesname.Data(),containerStr.Data()); TString outname = Form("coutput%s%s",mesname.Data(),containerStr.Data()); TString normname = Form("coutput%sNorm%s",mesname.Data(),containerStr.Data());