]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Updated PID for D0 at low pt (Christian Mohler, Andrea Rossi)
authorfprino <prino@to.infn.it>
Wed, 30 Jul 2014 12:56:02 +0000 (14:56 +0200)
committerfprino <prino@to.infn.it>
Wed, 30 Jul 2014 12:56:35 +0000 (14:56 +0200)
PWGHF/vertexingHF/AliAnalysisTaskCombinHF.cxx
PWGHF/vertexingHF/AliAnalysisTaskCombinHF.h
PWGHF/vertexingHF/macros/AddTaskCombinHF.C

index 24cd4c92702490c40ef6dddc930700521461ba1c..e556b880a222ca5e7439a5a12d92ad1c913830c3 100644 (file)
@@ -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<AliAODEvent*> (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<AliAODEvent*> (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; iTr<ntracks; iTr++){
     status[iTr]=0;
     AliAODTrack* track=aod->GetTrack(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; iTr1<ntracks; iTr1++){
     AliAODTrack* trK=aod->GetTrack(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; iTr2<ntracks; iTr2++){
       if((status[iTr2] & 1)==0) continue;
@@ -448,69 +479,69 @@ void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){
       AliAODTrack* trPi1=aod->GetTrack(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; iTr3<ntracks; iTr3++){
-         if((status[iTr3] & 1)==0) continue;
-         if((status[iTr3] & 4)==0) continue;
-         if(iTr1==iTr3) continue;
-         AliAODTrack* trPi2=aod->GetTrack(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; iTr3<ntracks; iTr3++){
+          if((status[iTr3] & 1)==0) continue;
+          if((status[iTr3] & 4)==0) continue;
+          if(iTr1==iTr3) continue;
+          AliAODTrack* trPi2=aod->GetTrack(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 && minv2<fMaxMass*fMaxMass){
     Double_t rapid = tmpRD->Y(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; iii<nProngs; iii++) signPdg[iii]=pdgdau[iii];
-       Int_t labD = tmpRD->MatchToMC(pdgD,arrayMC,nProngs,signPdg);
-       if(labD>=0){
-         AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(TMath::Abs(dgLabels[0])));
-         if(part){
-           Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
-           if(pdgCode==321){
-             AliAODMCParticle* dmes =  dynamic_cast<AliAODMCParticle*>(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; iii<nProngs; iii++) signPdg[iii]=pdgdau[iii];
+        Int_t labD = tmpRD->MatchToMC(pdgD,arrayMC,nProngs,signPdg);
+        if(labD>=0){
+          AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(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<AliAODMCParticle*>(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; irot<fNRotations; irot++){
     Double_t phirot=fMinAngleForRot+rotStep*irot;
     Double_t tmpx=px[0];
@@ -650,14 +680,14 @@ Bool_t AliAnalysisTaskCombinHF::FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoD
     if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
       Double_t rapid = tmpRD->Y(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<TList*> (GetOutputData(1));
-  if (!fOutput) {     
+  if (!fOutput) {
     printf("ERROR: fOutput not available\n");
     return;
   }
index 2aa1b3f29d0e402a92531302091779798b2d4f72..121b25b4538e1548fff8bd71b5c2117074c950fa 100644 (file)
@@ -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
 
 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
 };
 
index 38ec6707349bf66a51f20c0ed54ddf47634ad3e4..9cec59e23ac6f39bed7f68b5c81b0403e8994b00 100644 (file)
@@ -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());