2D D meson efficiency correction + updates for running on MC samples

diff --git a/PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.cxx b/PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.cxx
index 54a1304972a2c69c21ad1ecc7692fea944f71743..dfacaabeb1aba96f747cb875a7a91ba3dd01186c 100644 (file)
--- a/PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.cxx
+++ b/PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.cxx
@@ -52,6 +52,7 @@
 #include "AliAnalysisTaskSE.h"
 #include "AliAnalysisTaskSED0Correlations.h"
 #include "AliNormalizationCounter.h"
+#include "AliVertexingHFUtils.h"
 
 using std::cout;
 using std::endl;
@@ -66,7 +67,6 @@ AliAnalysisTaskSE(),
   fBinLimsCorr(),
   fPtThreshLow(),
   fPtThreshUp(), 
-  fD0Eff(),
   fEvents(0),
   fAlreadyFilled(kFALSE),
   fOutputMass(0),
@@ -90,7 +90,10 @@ AliAnalysisTaskSE(),
   fSys(0),
   fEtaForCorrel(0),
   fIsRejectSDDClusters(0),
-  fFillGlobal(kTRUE)
+  fFillGlobal(kTRUE),
+  fMultEv(0.),
+  fSoftPiCut(kTRUE),
+  fMEAxisThresh(kFALSE)
 {
   // Default constructor
 
@@ -103,7 +106,6 @@ AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const char *nam
   fBinLimsCorr(),
   fPtThreshLow(),
   fPtThreshUp(),
-  fD0Eff(),
   fEvents(0),
   fAlreadyFilled(kFALSE),
   fOutputMass(0),
@@ -127,7 +129,10 @@ AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const char *nam
   fSys(0),
   fEtaForCorrel(0),
   fIsRejectSDDClusters(0),
-  fFillGlobal(kTRUE)
+  fFillGlobal(kTRUE),
+  fMultEv(0.),
+  fSoftPiCut(kTRUE),
+  fMEAxisThresh(kFALSE)
 {
   // Default constructor
 
@@ -158,7 +163,6 @@ AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const AliAnalys
   fBinLimsCorr(source.fBinLimsCorr),
   fPtThreshLow(source.fPtThreshLow),
   fPtThreshUp(source.fPtThreshUp),
-  fD0Eff(source.fD0Eff),
   fEvents(source.fEvents),
   fAlreadyFilled(source.fAlreadyFilled),
   fOutputMass(source.fOutputMass),
@@ -182,7 +186,10 @@ AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const AliAnalys
   fSys(source.fSys),
   fEtaForCorrel(source.fEtaForCorrel),
   fIsRejectSDDClusters(source.fIsRejectSDDClusters),
-  fFillGlobal(source.fFillGlobal)
+  fFillGlobal(source.fFillGlobal),
+  fMultEv(source.fMultEv),
+  fSoftPiCut(source.fSoftPiCut),
+  fMEAxisThresh(source.fMEAxisThresh)
 {
   // Copy constructor
 }
@@ -239,7 +246,6 @@ AliAnalysisTaskSED0Correlations& AliAnalysisTaskSED0Correlations::operator=(cons
   fBinLimsCorr = orig.fBinLimsCorr;
   fPtThreshLow = orig.fPtThreshLow;
   fPtThreshUp = orig.fPtThreshUp;
-  fD0Eff = orig.fD0Eff;
   fEvents = orig.fEvents;
   fAlreadyFilled = orig.fAlreadyFilled;
   fOutputMass = orig.fOutputMass;
@@ -264,6 +270,9 @@ AliAnalysisTaskSED0Correlations& AliAnalysisTaskSED0Correlations::operator=(cons
   fEtaForCorrel = orig.fEtaForCorrel;
   fIsRejectSDDClusters = orig.fIsRejectSDDClusters;
   fFillGlobal = orig.fFillGlobal;
+  fMultEv = orig.fMultEv;
+  fSoftPiCut = orig.fSoftPiCut;
+  fMEAxisThresh = orig.fMEAxisThresh;
 
   return *this; //returns pointer of the class
 }
@@ -300,9 +309,9 @@ void AliAnalysisTaskSED0Correlations::UserCreateOutputObjects()
   fCorrelatorTr = new AliHFCorrelator("CorrelatorTr",fCutsTracks,fSys);
   fCorrelatorKc = new AliHFCorrelator("CorrelatorKc",fCutsTracks,fSys);
   fCorrelatorK0 = new AliHFCorrelator("CorrelatorK0",fCutsTracks,fSys);
-  fCorrelatorTr->SetDeltaPhiInterval(-TMath::Pi()/2+TMath::Pi()/32.,3*TMath::Pi()/2+TMath::Pi()/32.);// set the Delta Phi Interval you want 
-  fCorrelatorKc->SetDeltaPhiInterval(-TMath::Pi()/2+TMath::Pi()/32.,3*TMath::Pi()/2+TMath::Pi()/32.);
-  fCorrelatorK0->SetDeltaPhiInterval(-TMath::Pi()/2+TMath::Pi()/32.,3*TMath::Pi()/2+TMath::Pi()/32.);
+  fCorrelatorTr->SetDeltaPhiInterval(-TMath::Pi()/2,3*TMath::Pi()/2);// set the Delta Phi Interval you want (in this case -0.5Pi to 1.5 Pi)
+  fCorrelatorKc->SetDeltaPhiInterval(-TMath::Pi()/2,3*TMath::Pi()/2);
+  fCorrelatorK0->SetDeltaPhiInterval(-TMath::Pi()/2,3*TMath::Pi()/2);
   fCorrelatorTr->SetEventMixing(fMixing);// sets the analysis on a single event (kFALSE) or mixed events (kTRUE)
   fCorrelatorKc->SetEventMixing(fMixing);
   fCorrelatorK0->SetEventMixing(fMixing);
@@ -341,39 +350,40 @@ void AliAnalysisTaskSED0Correlations::UserCreateOutputObjects()
 
   TString nameMass=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassWg=" ",nameSgnWg=" ", nameBkgWg=" ", nameRflWg=" ";
 
+//for origin c case (or for data)
   for(Int_t i=0;i<fCutsD0->GetNPtBins();i++){
 
-    nameMass="histMass_";
+    nameMass="histMass_"; if(fReadMC) nameMass+="c_";
     nameMass+=i;
-    nameMassWg="histMass_WeigD0Eff_";
+    nameMassWg="histMass_WeigD0Eff_"; if(fReadMC) nameMassWg+="c_";
     nameMassWg+=i;
-    nameSgn="histSgn_";
+    nameSgn="histSgn_"; if(fReadMC) nameSgn+="c_";
     nameSgn+=i;
-    nameSgnWg="histSgn_WeigD0Eff_";
+    nameSgnWg="histSgn_WeigD0Eff_"; if(fReadMC) nameSgnWg+="c_";
     nameSgnWg+=i;
-    nameBkg="histBkg_";
+    nameBkg="histBkg_"; if(fReadMC) nameBkg+="c_";
     nameBkg+=i;
-    nameBkgWg="histBkg_WeigD0Eff_";
+    nameBkgWg="histBkg_WeigD0Eff_"; if(fReadMC) nameBkgWg+="c_";
     nameBkgWg+=i;
-    nameRfl="histRfl_";
+    nameRfl="histRfl_"; if(fReadMC) nameRfl+="c_";
     nameRfl+=i;
-    nameRflWg="histRfl_WeigD0Eff_";
+    nameRflWg="histRfl_WeigD0Eff_"; if(fReadMC) nameRflWg+="c_";
     nameRflWg+=i;
 
     //histograms of invariant mass distributions
 
     //MC signal
     if(fReadMC){
-      TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
-      TH1F* tmpStWg = new TH1F(nameSgnWg.Data(), "D^{0} invariant mass - MC; M [GeV] - weight 1/D0eff; Entries",120,1.5648,2.1648);
+      TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass c - MC; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpStWg = new TH1F(nameSgnWg.Data(), "D^{0} invariant mass c - MC; M [GeV] - weight 1/D0eff; Entries",120,1.5648,2.1648);
       tmpSt->Sumw2();
       tmpStWg->Sumw2();
 
       //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
-      TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
-      TH1F* tmpRtWg = new TH1F(nameRflWg.Data(), "Reflected signal invariant mass - MC - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
-      TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
-      TH1F* tmpBtWg = new TH1F(nameBkgWg.Data(), "Background invariant mass - MC - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass c - MC; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpRtWg = new TH1F(nameRflWg.Data(), "Reflected signal invariant mass c - MC - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass c - MC; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpBtWg = new TH1F(nameBkgWg.Data(), "Background invariant mass c - MC - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
       tmpBt->Sumw2();
       tmpBtWg->Sumw2();
       tmpRt->Sumw2();
@@ -387,15 +397,48 @@ void AliAnalysisTaskSED0Correlations::UserCreateOutputObjects()
     }
 
     //mass
-    TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",120,1.5648,2.1648);
+    TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass c; M [GeV]; Entries",120,1.5648,2.1648);
     tmpMt->Sumw2();
     fOutputMass->Add(tmpMt);
     //mass weighted by 1/D0eff
-    TH1F* tmpMtwg = new TH1F(nameMassWg.Data(),"D^{0} invariant mass - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
+    TH1F* tmpMtwg = new TH1F(nameMassWg.Data(),"D^{0} invariant mass c - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
     tmpMtwg->Sumw2();
     fOutputMass->Add(tmpMtwg);
   }
 
+//for origin b case (no Bkg and Mass histos, here for weights you should use c+b efficiencies, while on data (on MC they're useless))
+  for(Int_t i=0;i<fCutsD0->GetNPtBins();i++){
+
+    nameSgn="histSgn_b_";
+    nameSgn+=i;
+    nameSgnWg="histSgn_WeigD0Eff_b_";
+    nameSgnWg+=i;
+    nameRfl="histRfl_b_";
+    nameRfl+=i;
+    nameRflWg="histRfl_WeigD0Eff_b_";
+    nameRflWg+=i;
+
+    //histograms of invariant mass distributions
+
+    //MC signal
+    if(fReadMC){
+      TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass b - MC; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpStWg = new TH1F(nameSgnWg.Data(), "D^{0} invariant mass b - MC; M [GeV] - weight 1/D0eff; Entries",120,1.5648,2.1648);
+      tmpSt->Sumw2();
+      tmpStWg->Sumw2();
+
+      //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
+      TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass b - MC; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpRtWg = new TH1F(nameRflWg.Data(), "Reflected signal invariant mass b - MC - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
+      tmpRt->Sumw2();
+      tmpRtWg->Sumw2();
+      fOutputMass->Add(tmpSt);
+      fOutputMass->Add(tmpStWg);
+      fOutputMass->Add(tmpRt);
+      fOutputMass->Add(tmpRtWg);
+    }
+  }
+
   const char* nameoutput=GetOutputSlot(2)->GetContainer()->GetName();
 
   fNentries=new TH1F(nameoutput, "Integral(1,2) = number of AODs *** Integral(2,3) = number of candidates selected with cuts *** Integral(3,4) = number of D0 selected with cuts *** Integral(4,5) = events with good vertex ***  Integral(5,6) = pt out of bounds", 20,-0.5,19.5);
@@ -461,6 +504,8 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
 
   TClonesArray *inputArray=0;
 
+  fMultEv = 0.; //reset event multiplicity
+
   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.    
@@ -546,7 +591,7 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
     }
                 
     if(NMCevents && !isMCeventgood){
-      std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;
+      if(fDebug>2)std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;
       return; 
     }
     fNentries->Fill(19); //event with particular production type                
@@ -634,12 +679,15 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
 
   Int_t nSelectedloose=0,nSelectedtight=0;  
   
+  //Fill Event Multiplicity (needed only in Reco)
+  fMultEv = (Double_t)(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.));
+
   //RecoD0 case ************************************************
   if(fRecoD0) {
 
     for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
       AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
-
+ 
       if(d->Pt()<2.) continue; //to save time and merging memory...
 
       if(d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
@@ -671,12 +719,18 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
         fCorrelatorK0->SetD0Properties(d,fIsSelectedCandidate);
 
         if(!fReadMC) {
-          if (TMath::Abs(d->Eta())<fEtaForCorrel) CalculateCorrelations(d); //correlations on real data
+          if (TMath::Abs(d->Eta())<fEtaForCorrel) {
+           CalculateCorrelations(d); //correlations on real data
+           if(!fMixing) ((TH1F*)fOutputStudy->FindObject(Form("hMultiplEvt_Bin%d",ptbin)))->Fill(fMultEv);
+         }
         } else { //correlations on MC -> association of selected D0 to MCinfo with MCtruth
           if (TMath::Abs(d->Eta())<fEtaForCorrel) {
             Int_t pdgDgD0toKpi[2]={321,211};
            Int_t labD0 = d->MatchToMC(421,mcArray,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not
-            if (labD0>-1) CalculateCorrelations(d,labD0,mcArray);
+            if (labD0>-1) {
+             CalculateCorrelations(d,labD0,mcArray);
+             if(!fMixing) ((TH1F*)fOutputStudy->FindObject(Form("hMultiplEvt_Bin%d",ptbin)))->Fill(fMultEv); //Fill multiplicity histo
+           }
           }
         }
 
@@ -721,15 +775,20 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
          /*   Int_t mother = mcPart->GetMother();
            AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
             if(!mcMoth) continue;
-           if(TMath::Abs(mcMoth->GetPdgCode())==413) continue;
-       */
+           if(TMath::Abs(mcMoth->GetPdgCode())==413) continue;*/
+
             if (mcPart->GetPdgCode()==421) fIsSelectedCandidate = 1;
            else fIsSelectedCandidate = 2;
 
-           TString fillthis="histSgn_"; fillthis+=ptbin;
+           TString fillthis="histSgn_"; 
+           if(CheckD0Origin(mcArray,mcPart)==4) fillthis+="c_";
+           else if(CheckD0Origin(mcArray,mcPart)==5) fillthis+="b_";
+            else continue;
+            fillthis+=ptbin;
            ((TH1F*)(fOutputMass->FindObject(fillthis)))->Fill(1.864);
          
             CalculateCorrelationsMCKine(mcPart,mcArray);
+            if(!fMixing) ((TH1F*)fOutputStudy->FindObject(Form("hMultiplEvt_Bin%d",ptbin)))->Fill(fMultEv); //Fill multiplicity histo
           }
         }
       }
@@ -782,33 +841,53 @@ void AliAnalysisTaskSED0Correlations::FillMassHists(AliAODRecoDecayHF2Prong *par
   if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
 
     if(fReadMC){
-      if(labD0>=0) {
-       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
+      if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==4) {
+       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
        Int_t pdgD0 = partD0->GetPdgCode();
        if (pdgD0==421){ //D0
-         fillthis="histSgn_";
+         fillthis="histSgn_c_";
          fillthis+=ptbin;
          ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
-          fillthis="histSgn_WeigD0Eff_";
+          fillthis="histSgn_WeigD0Eff_c_";
           fillthis+=ptbin;
-          Double_t effD0 = fD0Eff.at(ptbin);
+          Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
-       } else{ //it was a D0bar
-         fillthis="histRfl_";
+       } else{ //it was a D0bar
+         fillthis="histRfl_c_";
          fillthis+=ptbin;
-         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
-          fillthis="histRfl_WeigD0Eff_";
+         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
+          fillthis="histRfl_WeigD0Eff_c_";
           fillthis+=ptbin;
-          Double_t effD0 = fD0Eff.at(ptbin);
+          Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
-       }
+       }
+      } else if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==5) {
+         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
+         Int_t pdgD0 = partD0->GetPdgCode();
+         if (pdgD0==421){ //D0
+           fillthis="histSgn_b_";
+           fillthis+=ptbin;
+           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
+            fillthis="histSgn_WeigD0Eff_b_";
+            fillthis+=ptbin;
+            Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
+            ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
+         } else{ //it was a D0bar
+           fillthis="histRfl_b_";
+           fillthis+=ptbin;
+           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
+            fillthis="histRfl_WeigD0Eff_b_";
+            fillthis+=ptbin;
+            Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
+            ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
+         }
       } else {//background
-       fillthis="histBkg_";
+       fillthis="histBkg_c_";
        fillthis+=ptbin;
        ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
-        fillthis="histBkg_WeigD0Eff_";
+        fillthis="histBkg_WeigD0Eff_c_";
         fillthis+=ptbin;
-        Double_t effD0 = fD0Eff.at(ptbin);
+        Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
       }
     }else{
@@ -817,7 +896,7 @@ void AliAnalysisTaskSED0Correlations::FillMassHists(AliAODRecoDecayHF2Prong *par
       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
       fillthis="histMass_WeigD0Eff_";
       fillthis+=ptbin;
-      Double_t effD0 = fD0Eff.at(ptbin);
+      Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
     }
      
@@ -825,45 +904,65 @@ void AliAnalysisTaskSED0Correlations::FillMassHists(AliAODRecoDecayHF2Prong *par
   if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
 
     if(fReadMC){
-      if(labD0>=0) {
-       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
+      if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==4) {
+       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
        Int_t pdgD0 = partD0->GetPdgCode();
-
-       if (pdgD0==-421){ //D0bar
-         fillthis="histSgn_";
+       if (pdgD0==-421){ //D0
+         fillthis="histSgn_c_";
          fillthis+=ptbin;
          ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
-          fillthis="histSgn_WeigD0Eff_";
+          fillthis="histSgn_WeigD0Eff_c_";
           fillthis+=ptbin;
-          Double_t effD0 = fD0Eff.at(ptbin);
+          Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
-       } else{
-         fillthis="histRfl_";
+       } else{ //it was a D0bar
+         fillthis="histRfl_c_";
          fillthis+=ptbin;
          ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
-          fillthis="histRfl_WeigD0Eff_";
+          fillthis="histRfl_WeigD0Eff_c_";
           fillthis+=ptbin;
-          Double_t effD0 = fD0Eff.at(ptbin);
+          Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
-       }
-      } else {//background or LS
-       fillthis="histBkg_";
+       }
+      } else if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==5) {
+         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
+         Int_t pdgD0 = partD0->GetPdgCode();
+         if (pdgD0==-421){ //D0
+           fillthis="histSgn_b_";
+           fillthis+=ptbin;
+           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
+            fillthis="histSgn_WeigD0Eff_b_";
+            fillthis+=ptbin;
+            Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
+            ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
+         } else{ //it was a D0bar
+           fillthis="histRfl_b_";
+           fillthis+=ptbin;
+           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
+            fillthis="histRfl_WeigD0Eff_b_";
+            fillthis+=ptbin;
+            Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
+            ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
+         }
+      } else {//background
+       fillthis="histBkg_c_";
        fillthis+=ptbin;
        ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
-        fillthis="histBkg_WeigD0Eff_";
+        fillthis="histBkg_WeigD0Eff_c_";
         fillthis+=ptbin;
-        Double_t effD0 = fD0Eff.at(ptbin);
+        Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
       }
     }else{
       fillthis="histMass_";
       fillthis+=ptbin;
-      ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
+      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
       fillthis="histMass_WeigD0Eff_";
       fillthis+=ptbin;
-      Double_t effD0 = fD0Eff.at(ptbin);
+      Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
     }
+
   }
 
   return;
@@ -965,13 +1064,13 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
   //These for limits in THnSparse (one per bin, same limits). 
   //Vars: DeltaPhi, InvMass, PtTrack, Displacement, DeltaEta
   Int_t nBinsPhi[5] = {32,150,6,3,16};
-  Double_t binMinPhi[5] = {-TMath::Pi()/2.+TMath::Pi()/32.,1.6,0.,0.,-1.6};  //is the minimum for all the bins
-  Double_t binMaxPhi[5] = {3*TMath::Pi()/2.+TMath::Pi()/32.,2.2,3.0,3.,1.6};  //is the maximum for all the bins
+  Double_t binMinPhi[5] = {-TMath::Pi()/2.,1.6,0.,0.,-1.6};  //is the minimum for all the bins
+  Double_t binMaxPhi[5] = {3.*TMath::Pi()/2.,2.2,3.0,3.,1.6};  //is the maximum for all the bins
 
   //Vars: DeltaPhi, InvMass, DeltaEta
-  Int_t nBinsMix[3] = {32,150,16};
-  Double_t binMinMix[3] = {-TMath::Pi()/2.+TMath::Pi()/32.,1.6,-1.6};  //is the minimum for all the bins
-  Double_t binMaxMix[3] = {3*TMath::Pi()/2.+TMath::Pi()/32.,2.2,1.6};  //is the maximum for all the bins
+  Int_t nBinsMix[4] = {32,150,16,6};
+  Double_t binMinMix[4] = {-TMath::Pi()/2.,1.6,-1.6,0.};  //is the minimum for all the bins
+  Double_t binMaxMix[4] = {3.*TMath::Pi()/2.,2.2,1.6,3.};  //is the maximum for all the bins
 
   for(Int_t i=0;i<fNPtBinsCorr;i++){
 
@@ -1113,7 +1212,7 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
       namePlot="hPhi_Lead_Bin";
       namePlot+=i;
 
-      THnSparseF *hCorrLead = new THnSparseF(namePlot.Data(), "Leading particle correlations; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+      THnSparseF *hCorrLead = new THnSparseF(namePlot.Data(), "Leading particle correlations; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
       hCorrLead->Sumw2();
       fOutputCorr->Add(hCorrLead);
 
@@ -1121,35 +1220,35 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
         namePlot="hPhi_Lead_From_c_Bin";
         namePlot+=i;
 
-        THnSparseF *hCorrLead_c = new THnSparseF(namePlot.Data(), "Leading particle correlations - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+        THnSparseF *hCorrLead_c = new THnSparseF(namePlot.Data(), "Leading particle correlations - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
         hCorrLead_c->Sumw2();
         fOutputCorr->Add(hCorrLead_c);
   
         namePlot="hPhi_Lead_From_b_Bin";
         namePlot+=i;
   
-        THnSparseF *hCorrLead_b = new THnSparseF(namePlot.Data(), "Leading particle correlations - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+        THnSparseF *hCorrLead_b = new THnSparseF(namePlot.Data(), "Leading particle correlations - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
         hCorrLead_b->Sumw2();
         fOutputCorr->Add(hCorrLead_b);
   
         namePlot="hPhi_Lead_HF_From_c_Bin";
         namePlot+=i;
   
-        THnSparseF *hCorrLead_HF_c = new THnSparseF(namePlot.Data(), "Leading particle correlations HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+        THnSparseF *hCorrLead_HF_c = new THnSparseF(namePlot.Data(), "Leading particle correlations HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
         hCorrLead_HF_c->Sumw2();
         fOutputCorr->Add(hCorrLead_HF_c);
   
         namePlot="hPhi_Lead_HF_From_b_Bin";
         namePlot+=i;
   
-        THnSparseF *hCorrLead_HF_b = new THnSparseF(namePlot.Data(), "Leading particle correlations HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+        THnSparseF *hCorrLead_HF_b = new THnSparseF(namePlot.Data(), "Leading particle correlations HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
         hCorrLead_HF_b->Sumw2();
         fOutputCorr->Add(hCorrLead_HF_b);
 
         namePlot="hPhi_Lead_NonHF_Bin";
         namePlot+=i;
   
-        THnSparseF *hCorrLead_Non = new THnSparseF(namePlot.Data(), "Leading particle correlations - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+        THnSparseF *hCorrLead_Non = new THnSparseF(namePlot.Data(), "Leading particle correlations - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
         hCorrLead_Non->Sumw2();
         fOutputCorr->Add(hCorrLead_Non);
       }
@@ -1158,38 +1257,38 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
       namePlot="hPhi_Weig_Bin";
       namePlot+=i;
   
-      THnSparseF *hCorrWeig = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted); #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+      THnSparseF *hCorrWeig = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted); #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
       fOutputCorr->Add(hCorrWeig);
   
       if (fReadMC) {
         namePlot="hPhi_Weig_From_c_Bin";
         namePlot+=i;
   
-        THnSparseF *hCorrWeig_c = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+        THnSparseF *hCorrWeig_c = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
         fOutputCorr->Add(hCorrWeig_c);
   
         namePlot="hPhi_Weig_From_b_Bin";
         namePlot+=i;
   
-        THnSparseF *hCorrWeig_b = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+        THnSparseF *hCorrWeig_b = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
         fOutputCorr->Add(hCorrWeig_b);
   
         namePlot="hPhi_Weig_HF_From_c_Bin";
         namePlot+=i;
   
-        THnSparseF *hCorrWeig_HF_c = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+        THnSparseF *hCorrWeig_HF_c = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
         fOutputCorr->Add(hCorrWeig_HF_c);
   
         namePlot="hPhi_Weig_HF_From_b_Bin";
         namePlot+=i;
   
-        THnSparseF *hCorrWeig_HF_b = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+        THnSparseF *hCorrWeig_HF_b = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
         fOutputCorr->Add(hCorrWeig_HF_b);
 
         namePlot="hPhi_Weig_NonHF_Bin";
         namePlot+=i;
   
-        THnSparseF *hCorrWeig_Non = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+        THnSparseF *hCorrWeig_Non = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
         fOutputCorr->Add(hCorrWeig_Non);
       }
 
@@ -1216,6 +1315,13 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
     hDstarPions->GetXaxis()->SetBinLabel(2,"Rejected");
     hDstarPions->SetMinimum(0);
     fOutputStudy->Add(hDstarPions); 
+
+    //Events multiplicity
+    Double_t yAxisMult[13] = {0, 4, 8, 12, 16, 20, 28, 36, 44, 100}; 
+    namePlot = "hMultiplEvt_Bin"; namePlot+=i;
+    TH1F *hMultEv = new TH1F(namePlot.Data(), "Event multiplicity",9,yAxisMult);
+    hMultEv->SetMinimum(0);
+    fOutputStudy->Add(hMultEv);
  
     }
 
@@ -1224,21 +1330,21 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
       namePlot="hPhi_K0_Bin";
       namePlot+=i;namePlot+="_EvMix";
 
-      THnSparseF *hPhiK_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+      THnSparseF *hPhiK_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
       hPhiK_EvMix->Sumw2();
       fOutputCorr->Add(hPhiK_EvMix);
 
       namePlot="hPhi_Kcharg_Bin";
       namePlot+=i;namePlot+="_EvMix";
   
-      THnSparseF *hPhiH_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+      THnSparseF *hPhiH_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
       hPhiH_EvMix->Sumw2();
       fOutputCorr->Add(hPhiH_EvMix);
 
       namePlot="hPhi_Charg_Bin";
       namePlot+=i;namePlot+="_EvMix";
 
-      THnSparseF *hPhiC_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+      THnSparseF *hPhiC_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
       hPhiC_EvMix->Sumw2();
       fOutputCorr->Add(hPhiC_EvMix);
     }
@@ -1452,6 +1558,14 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
       hOrigin_D0->GetXaxis()->SetBinLabel(1,"From c");
       hOrigin_D0->GetXaxis()->SetBinLabel(2,"From b");
       fOutputStudy->Add(hOrigin_D0);
+
+      //primary tracks (Kine & Reco)
+      namePlot="hPhysPrim_Bin";  namePlot+=i;
+      TH1F *hPhysPrim = new TH1F(namePlot.Data(), "Origin of hadrons",2,0.,2.);
+      hPhysPrim->SetMinimum(0);
+      hPhysPrim->GetXaxis()->SetBinLabel(1,"OK");
+      hPhysPrim->GetXaxis()->SetBinLabel(2,"NO");
+      fOutputStudy->Add(hPhysPrim);
     }
   }
 }
@@ -1524,7 +1638,7 @@ void AliAnalysisTaskSED0Correlations::CalculateCorrelations(AliAODRecoDecayHF2Pr
       AliReducedParticle* track = fCorrelatorTr->GetAssociatedParticle();
 
       if(!fMixing) {      
-       if(!track->CheckSoftPi()) { //removal of soft pions
+       if(fSoftPiCut && !track->CheckSoftPi()) { //removal of soft pions
           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
           continue;
@@ -1537,11 +1651,24 @@ void AliAnalysisTaskSED0Correlations::CalculateCorrelations(AliAODRecoDecayHF2Pr
       }
       if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
 
+      if(fReadMC) {
+        AliAODMCParticle* trkKine = (AliAODMCParticle*)mcArray->At(track->GetLabel());
+        if (!trkKine) continue;
+        if (!trkKine->IsPhysicalPrimary()) {
+         ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(1.);  
+         continue; //reject the Reco track if correspondent Kine track is not primary
+        } else ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(0.);
+      }
+
       Double_t effTr = track->GetWeight(); //extract track efficiency
-      Double_t effD0 = fD0Eff.at(ptbin); 
+      Double_t effD0 = 1.;
+      if(fReadMC) {
+        if(origD0==4) effD0 = fCutsTracks->GetTrigWeight(d->Pt(),fMultEv);
+        if(origD0==5) effD0 = fCutsTracks->GetTrigWeightB(d->Pt(),fMultEv);
+      } else effD0 = fCutsTracks->GetTrigWeight(d->Pt(),fMultEv);
       Double_t eff = effTr*effD0;
- 
-      FillSparsePlots(mcArray,mInv,origD0,PDGD0,track,ptbin,kTrack,1./eff); //fills for charged tracks, weight = 1./eff
+
+      FillSparsePlots(mcArray,mInv,origD0,PDGD0,track,ptbin,kTrack,1./eff); //fills for charged tracks
 
       if(!fMixing) N_Charg++;
 
@@ -1552,7 +1679,10 @@ void AliAnalysisTaskSED0Correlations::CalculateCorrelations(AliAODRecoDecayHF2Pr
         lead[0] = fCorrelatorTr->GetDeltaPhi();
         lead[1] = fCorrelatorTr->GetDeltaEta();
         lead[2] = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
-       lead[3] = 1./track->GetWeight(); //weight is 1./efficiency
+        if(fReadMC) {
+         if(origD0==4) lead[3] = 1./(track->GetWeight()*fCutsTracks->GetTrigWeight(d->Pt(),fMultEv)); //weight is 1./efficiency
+         if(origD0==5) lead[3] = 1./(track->GetWeight()*fCutsTracks->GetTrigWeightB(d->Pt(),fMultEv)); //weight is 1./efficiency
+       } else lead[3] = 1./(track->GetWeight()*fCutsTracks->GetTrigWeight(d->Pt(),fMultEv));
         highPt = track->Pt();
       }
 
@@ -1583,7 +1713,7 @@ void AliAnalysisTaskSED0Correlations::CalculateCorrelations(AliAODRecoDecayHF2Pr
       AliReducedParticle* kCharg = fCorrelatorKc->GetAssociatedParticle();
 
       if(!fMixing) {      
-       if(!kCharg->CheckSoftPi()) { //removal of soft pions
+       if(fSoftPiCut && !kCharg->CheckSoftPi()) { //removal of soft pions
           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
           continue;
@@ -1635,20 +1765,20 @@ void AliAnalysisTaskSED0Correlations::CalculateCorrelations(AliAODRecoDecayHF2Pr
     } // end of charged kaons loop
   } //end of event loop for fCorrelatorK0
 
-  Double_t fillSpLeadD0[3] = {lead[0],mD0,lead[1]};
-  Double_t fillSpLeadD0bar[3] = {lead[0],mD0bar,lead[1]};
+  Double_t fillSpLeadD0[4] = {lead[0],mD0,lead[1],0.4}; //dummy value for threshold of leading!
+  Double_t fillSpLeadD0bar[4] = {lead[0],mD0bar,lead[1],0.4};
 
   //leading track correlations fill
   if(!fMixing) {
     if(fReadMC) {
-      if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0
+      if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3)) { //D0
         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]); //c and b D0
         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]); //c or b D0
         if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]);  
         if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]);  
         if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]);  //non HF  
       }
-      if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
+      if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421 && fIsSelectedCandidate>1) { //D0bar
         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]); //c or b D0
         if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]);  
@@ -1744,8 +1874,14 @@ void AliAnalysisTaskSED0Correlations::CalculateCorrelationsMCKine(AliAODMCPartic
       AliAODMCParticle *trkMC = (AliAODMCParticle*)mcArray->At(track->GetLabel());
       if(!trkMC) continue;
 
-      if (!trkMC->IsPhysicalPrimary()) continue; //reject material budget, or other fake tracks
+      if (!trkMC->IsPhysicalPrimary()) {  //reject material budget, or other fake tracks
+       ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(1.);  
+       continue;
+      } else ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(0.);
+
       if (IsDDaughter(d,trkMC,mcArray)) continue;
+      if (fSoftPiCut && IsSoftPion_MCKine(d,trkMC,mcArray)) continue; //remove soft pions (if requestes, e.g. for templates)
+
       FillSparsePlots(mcArray,mInv,origD0,PDGD0,track,ptbin,kTrack); //fills for charged tracks
 
       //retrieving leading info...
@@ -1832,18 +1968,18 @@ void AliAnalysisTaskSED0Correlations::CalculateCorrelationsMCKine(AliAODMCPartic
     } // end of charged kaons loop
   } //end of event loop for fCorrelatorK0
 
-  Double_t fillSpLeadMC[3] = {lead[0],mD0,lead[1]}; //mD0 = mD0bar = 1.864
+  Double_t fillSpLeadMC[4] = {lead[0],mD0,lead[1],0.4}; //mD0 = mD0bar = 1.864
 
   //leading track correlations fill
   if(!fMixing) {
-    if(d->GetPdgCode()==421) { //D0
+    if(d->GetPdgCode()==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3)) { //D0
       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadMC); //c and b D0
       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC); //c or b D0
       if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);  
       if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);  
       if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadMC);  //non HF
     }
-    if(d->GetPdgCode()==-421) { //D0bar
+    if(d->GetPdgCode()==-421 && fIsSelectedCandidate>1) { //D0bar
       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadMC);
       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC); //c or b D0
       if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);  
@@ -1920,8 +2056,8 @@ void AliAnalysisTaskSED0Correlations::FillSparsePlots(TClonesArray* mcArray, Dou
     //variables for filling histos
     Double_t fillSpPhiD0[5] = {deltaphi,mD0,ptTrack,d0Track,deltaeta};
     Double_t fillSpPhiD0bar[5] = {deltaphi,mD0bar,ptTrack,d0Track,deltaeta};
-    Double_t fillSpWeigD0[3] = {deltaphi,mD0,deltaeta};
-    Double_t fillSpWeigD0bar[3] = {deltaphi,mD0bar,deltaeta};
+    Double_t fillSpWeigD0[4] = {deltaphi,mD0,deltaeta,ptTrack};
+    Double_t fillSpWeigD0bar[4] = {deltaphi,mD0bar,deltaeta,ptTrack};
 
     if(fReadMC == 0) {
       //sparse fill for data (tracks, K+-, K0) + weighted
@@ -1944,7 +2080,7 @@ void AliAnalysisTaskSED0Correlations::FillSparsePlots(TClonesArray* mcArray, Dou
       if(origD0==4) {orig = "_From_c";} else {orig = "_From_b";}
 
       //sparse fill for data (tracks, K+-, K0) + weighted
-      if(PdgD0==421) { //D0 (from MCTruth)
+      if(PdgD0==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3)) { //D0 (from MCTruth)
          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
          if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
@@ -1956,7 +2092,7 @@ void AliAnalysisTaskSED0Correlations::FillSparsePlots(TClonesArray* mcArray, Dou
          if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
          if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_NonHF_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
       }
-      if(PdgD0==-421) { //D0bar
+      if(PdgD0==-421 && fIsSelectedCandidate>1) { //D0bar
          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
          if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
@@ -1987,10 +2123,16 @@ void AliAnalysisTaskSED0Correlations::FillSparsePlots(TClonesArray* mcArray, Dou
     Double_t EtaLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ptbin)))->GetAxis(2)->GetXmax();
     if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
     if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
+    Double_t ptLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ptbin)))->GetAxis(3)->GetXmax(); //all plots have same axes...
+    if(ptTrack > ptLim_Sparse) ptTrack = ptLim_Sparse-0.01;
 
     //variables for filling histos
-    Double_t fillSpPhiD0[3] = {deltaphi,mD0,deltaeta};
-    Double_t fillSpPhiD0bar[3] = {deltaphi,mD0bar,deltaeta};
+    Double_t fillSpPhiD0[4] = {deltaphi,mD0,deltaeta,0.4}; //dummy for ME threshold! unless explicitly set by flag...
+    Double_t fillSpPhiD0bar[4] = {deltaphi,mD0bar,deltaeta,0.4};
+    if(fMEAxisThresh) {
+      fillSpPhiD0[3] = ptTrack;
+      fillSpPhiD0bar[3] = ptTrack;
+    }
 
     if(fReadMC == 0) {
       //sparse fill for data (tracks, K+-, K0)
@@ -1999,8 +2141,8 @@ void AliAnalysisTaskSED0Correlations::FillSparsePlots(TClonesArray* mcArray, Dou
     }
     if(fReadMC == 1) {
       //sparse fill for data (tracks, K+-, K0)
-      if(PdgD0==421)  ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
-      if(PdgD0==-421) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
+      if(PdgD0==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3))  ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
+      if(PdgD0==-421 && fIsSelectedCandidate>1) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
     }//end MC case
 
   } //end of ME fill
@@ -2152,6 +2294,37 @@ Bool_t AliAnalysisTaskSED0Correlations::SelectV0(AliAODv0* v0, AliAODVertex *vtx
   return kTRUE;
 }
 
+//---------------------------------------------------------------------------
+Bool_t AliAnalysisTaskSED0Correlations::IsSoftPion_MCKine(AliAODMCParticle* d, AliAODMCParticle* track, TClonesArray* arrayMC) const
+{
+  //
+  // Removes soft pions in Kine
+
+  //Daughter removal in MCKine case
+  Bool_t isSoftPi = kFALSE;
+  Int_t labelD0 = d->GetLabel();
+
+  Int_t mother = track->GetMother();
+  if(!mother) return isSoftPi; //safety check
+
+  AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother)); //it's the mother of the track!
+  if(TMath::Abs(mcMoth->GetPdgCode())==413 && mcMoth->GetNDaughters()==2) { //mother is D* with 2 daughs
+    Int_t labdau1 = mcMoth->GetDaughter(0);
+    Int_t labdau2 = mcMoth->GetDaughter(1);
+    AliAODMCParticle* dau1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(labdau1));
+    AliAODMCParticle* dau2 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(labdau2));
+    if(!dau1 || !dau2) return isSoftPi; //safety check
+    if(dau1->GetLabel()==labelD0 || dau2->GetLabel()==labelD0) { //one of the daughs is the D0 trigger
+      if((TMath::Abs(dau1->GetPdgCode())==421 && TMath::Abs(dau2->GetPdgCode())==211)||(TMath::Abs(dau1->GetPdgCode())==211 && TMath::Abs(dau2->GetPdgCode())==421)) {
+       isSoftPi = kTRUE; //ok, soft pion was found
+       return isSoftPi;
+      }
+    }
+  } 
+
+  return isSoftPi;
+}
+
 //________________________________________________________________________
 void AliAnalysisTaskSED0Correlations::PrintBinsAndLimits() {
 
@@ -2170,10 +2343,6 @@ void AliAnalysisTaskSED0Correlations::PrintBinsAndLimits() {
   for (int i=0; i<fNPtBinsCorr; i++) {
     cout << fPtThreshUp.at(i) << "; ";
   }
-  cout << "D0 efficiencies:----\n";
-  for (int i=0; i<fNPtBinsCorr; i++) {
-    cout << fD0Eff.at(i) << "; ";
-  }
   cout << "\n--------------------------\n";
   cout << "D0 Eta cut for Correl = "<<fEtaForCorrel<<"\n";
   cout << "--------------------------\n";
@@ -2183,5 +2352,9 @@ void AliAnalysisTaskSED0Correlations::PrintBinsAndLimits() {
   cout << "--------------------------\n";
   cout << "Ev Mixing = "<<fMixing<<"\n";
   cout << "--------------------------\n";
+  cout << "ME thresh axis = "<<fMEAxisThresh<<"\n";
+  cout << "--------------------------\n";
+  cout << "Soft Pi Cut = "<<fSoftPiCut<<"\n";
+  cout << "--------------------------\n";
 }
 

diff --git a/PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.h b/PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.h
index e50f5abb1db2291fc4bc242f9e3f836ec44e4a46..dc5a4ebce75eb6b11355210f5c79f246875aef8e 100644 (file)
--- a/PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.h
+++ b/PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.h
@@ -4,7 +4,7 @@
 /* Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
-/* $Id$ */ 
+/* $Id$ */
 
 //*************************************************************************
 // Class AliAnalysisTaskSED0Correlations
@@ -54,6 +54,8 @@ class AliAnalysisTaskSED0Correlations : public AliAnalysisTaskSE
   void SetSystem(Int_t sys){fSys=sys;}
   void SetRejectSDDClusters(Bool_t flag) {fIsRejectSDDClusters=flag; }
   void SetFillGlobalPlots(Bool_t fill=kTRUE){fFillGlobal=fill;}
+  void SetSoftPiFlag(Bool_t piflag) {fSoftPiCut=piflag;}
+  void SetMEAxisThresh(Bool_t methresh) {fMEAxisThresh=methresh;}
 
   Int_t  GetReadMC() const {return fReadMC;}
   Int_t  GetMCReconstructedTracks() const {return fRecoTr;}
@@ -64,6 +66,9 @@ class AliAnalysisTaskSED0Correlations : public AliAnalysisTaskSE
   Bool_t GetRejectSDDClusters() const {return fIsRejectSDDClusters;}
   Bool_t GetFillGlobalPlots() const {return fFillGlobal;}
   Double_t GetEtaForCorrel() {return fEtaForCorrel;}
+  Double_t GetMultEv() {return fMultEv;}
+  Bool_t GetSoftPiFlag() const {return fSoftPiCut;}
+  Bool_t GetMEAxisThresh() const {return fMEAxisThresh;}
 
   //correlations setters/printers
   void SetNPtBinsCorr(Int_t nbins) {fNPtBinsCorr = nbins;}
@@ -71,7 +76,6 @@ class AliAnalysisTaskSED0Correlations : public AliAnalysisTaskSE
   void SetPtBinsLimsCorr(Float_t* ptlims) {for(int i=0;i<=fNPtBinsCorr;i++) {fBinLimsCorr.push_back((Double_t)ptlims[i]);}}
   void SetPtTreshLow(Double_t* pttreshlow) {for(int i=0;i<fNPtBinsCorr;i++) {fPtThreshLow.push_back(pttreshlow[i]);}}
   void SetPtTreshUp(Double_t* pttreshup) {for(int i=0;i<fNPtBinsCorr;i++) {fPtThreshUp.push_back(pttreshup[i]);}}
-  void SetD0Efficiency(Double_t* effD0) {for(int i=0;i<fNPtBinsCorr;i++) {fD0Eff.push_back(effD0[i]);}}
 
   void PrintBinsAndLimits();
   Int_t PtBinCorr(Double_t pt) const;
@@ -95,12 +99,12 @@ class AliAnalysisTaskSED0Correlations : public AliAnalysisTaskSE
   Int_t CheckTrackOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const;
   Bool_t IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track, TClonesArray* mcArray) const;
   Bool_t SelectV0(AliAODv0* v0, AliAODVertex *vtx, Int_t option, Int_t idArrayV0[][2]) const;
+  Bool_t IsSoftPion_MCKine(AliAODMCParticle* d, AliAODMCParticle* track, TClonesArray* arrayMC) const;
 
   Int_t             fNPtBinsCorr;        // number of pt bins per correlations
   std::vector<Double_t>  fBinLimsCorr;        // limits of pt bins per correlations
   std::vector<Double_t>  fPtThreshLow;        // pT treshold of hadrons - low
   std::vector<Double_t>  fPtThreshUp;         // pT treshold of hadrons - up
-  std::vector<Double_t>  fD0Eff;              // D0 efficiencies
 
   Int_t     fEvents;                   // EventCounter
   Bool_t    fAlreadyFilled;            // D0 in an event already analyzed (for track distribution plots)
@@ -126,8 +130,11 @@ class AliAnalysisTaskSED0Correlations : public AliAnalysisTaskSE
   Double_t  fEtaForCorrel;             // cut for D0 eta to enable correlation with associated particles
   Bool_t    fIsRejectSDDClusters;      // flag to reject events with SDD clusters
   Bool_t    fFillGlobal;               // flag to fill global plots (in loops on tracks and V0 for each event)
+  Double_t  fMultEv;                   // event multiplicity (for trigger eff)
+  Bool_t    fSoftPiCut;                        // flag to activate soft pion cut on Data
+  Bool_t    fMEAxisThresh;             // flag to fill threshold axis in ME plots
 
-  ClassDef(AliAnalysisTaskSED0Correlations,3); // AliAnalysisTaskSE for D0->Kpi
+  ClassDef(AliAnalysisTaskSED0Correlations,4); // AliAnalysisTaskSE for D0->Kpi
 };
 
 #endif

diff --git a/PWGHF/correlationHF/macros/AddTaskD0Correlations.C b/PWGHF/correlationHF/macros/AddTaskD0Correlations.C
index fd225111791c6007073633dc1644607ddba15f87..f739ce306a395a89956ef92a0d8b5f6f045bb46b 100644 (file)
--- a/PWGHF/correlationHF/macros/AddTaskD0Correlations.C
+++ b/PWGHF/correlationHF/macros/AddTaskD0Correlations.C
@@ -1,6 +1,4 @@
-AliAnalysisTaskSED0Correlations *AddTaskD0Correlations(Bool_t readMC=kFALSE, Bool_t mixing=kFALSE, Bool_t recoTrMC=kFALSE, Bool_t recoD0MC = kFALSE, Double_t etacorr=1.5, Int_t 
-system=0/*0=pp,1=PbPb*/, Int_t flagD0D0bar=0, Float_t minC=0, Float_t maxC=0, TString finDirname="Output", TString cutsfilename="D0toKpiCuts.root", TString cutsfilename2="AssocPartCuts_fBit0_woITS.root", TString 
-cutsD0name="D0toKpiCuts", TString cutsTrkname="AssociatedTrkCuts", TString effName = "3D_eff_wo_ITScls2_f0_p8eta.root", Bool_t flagAOD049=kFALSE, Int_t standardbins=1, Bool_t stdcuts=kFALSE, TString effD0name="D0Eff_From_c_wLimAcc.root")
+AliAnalysisTaskSED0Correlations *AddTaskD0Correlations(Bool_t readMC=kFALSE, Bool_t mixing=kFALSE, Bool_t recoTrMC=kFALSE, Bool_t recoD0MC = kFALSE,  Bool_t flagsoftpicut = kTRUE, Bool_t MEthresh = kFALSE, TString cutsfilename="D0toKpiCuts.root", TString cutsfilename2="AssocPartCuts_fBit0_woITS_EffMap2D.root", TString effD0namec="D0Eff_From_c_2D.root", TString effD0nameb="D0Eff_From_b_2D.root", TString effName = "3D_eff_wo_ITScls2_f0_p8eta_FineBins.root", TString cutsD0name="D0toKpiCuts", TString cutsTrkname="AssociatedTrkCuts", Double_t etacorr=1.5, Int_t system=0/*0=pp,1=PbPb*/, Int_t flagD0D0bar=0, Float_t minC=0, Float_t maxC=0, TString finDirname="Output", Bool_t flagAOD049=kFALSE, Int_t standardbins=1, Bool_t stdcuts=kFALSE)
 {
   //
   // AddTask for the AliAnalysisTaskSE for D0 candidates
@@ -79,14 +77,12 @@ cutsD0name="D0toKpiCuts", TString cutsTrkname="AssociatedTrkCuts", TString effNa
     //     printf("    d0d0  [cm^2] < %f\n",fD0CorrCuts[7]);
     //     printf("    cosThetaPoint    > %f\n",fD0CorrCuts[8]);
 
-  TFile* filecuts;
-  filecuts=TFile::Open(cutsfilename.Data());
+  TFile* filecuts=new TFile(cutsfilename.Data());
   if(!filecuts->IsOpen()){
     cout<<"Input file not found for D0 cuts: using std cut object"<<endl;
     stdcuts=kTRUE;
   }
-  TFile* filecuts2;
-  filecuts2=TFile::Open(cutsfilename2.Data());
+  TFile* filecuts2=new TFile(cutsfilename2.Data());
   if(!filecuts2->IsOpen()){
     cout<<"Input file not found for tracks cuts!"<<endl;
     return;
@@ -122,8 +118,7 @@ cutsD0name="D0toKpiCuts", TString cutsTrkname="AssociatedTrkCuts", TString effNa
   }
 
   //Load efficiency map
-  TFile* fileeff;
-  fileeff=TFile::Open(effname.Data());
+  TFile* fileeff=new TFile(effName.Data());
   if(!fileeff->IsOpen()){
     cout<<"Input file not found for efficiency! Exiting..."<<endl;
     return;
@@ -131,6 +126,26 @@ cutsD0name="D0toKpiCuts", TString cutsTrkname="AssociatedTrkCuts", TString effNa
   TCanvas *c = (TCanvas*)fileeff->Get("c");
   TH3D *h3D = (TH3D*)c->FindObject("heff_rebin");
 
+  //Load D0 efficiency map
+  TFile* fileeffD0c=new TFile(effD0namec.Data());
+  if(!fileeffD0c->IsOpen()){
+    cout<<"Input file not found for efficiency! Exiting..."<<endl;
+    return;
+  }
+  TCanvas *cc = (TCanvas*)fileeffD0c->Get("c1");
+  TH2D *hEffD0c = (TH2D*)cc->FindObject("h_Eff");
+
+  //Load D0 efficiency map from b
+  if(readMC) {
+    TFile* fileeffD0b=new TFile(effD0nameb.Data());
+    if(!fileeffD0b->IsOpen()){
+      cout<<"Input file not found for efficiency! Exiting..."<<endl;
+      return;
+    }
+    TCanvas *cb = (TCanvas*)fileeffD0b->Get("c1");
+    TH2D *hEffD0b = (TH2D*)cb->FindObject("h_Eff");
+  }
+
   //Cuts for correlated tracks/K0
   AliHFAssociatedTrackCuts* corrCuts=new AliHFAssociatedTrackCuts();
   corrCuts = (AliHFAssociatedTrackCuts*)filecuts2->Get(cutsTrkname.Data());
@@ -140,7 +155,9 @@ cutsD0name="D0toKpiCuts", TString cutsTrkname="AssociatedTrkCuts", TString effNa
   }
   corrCuts->SetTrackCutsNames();
   corrCuts->SetvZeroCutsNames();
-  corrCuts->SetEfficiencyWeightMap(h3D);
+  if(recoD0MC) corrCuts->SetEfficiencyWeightMap(h3D); //data and MC Reco
+  if(recoD0MC) corrCuts->SetTriggerEffWeightMap(hEffD0c); //data and MC Reco
+  if(recoD0MC && readMC) corrCuts->SetTriggerEffWeightMapB(hEffD0b); //MC Reco
   corrCuts->PrintAll();
 
   TString centr="";
@@ -166,6 +183,8 @@ cutsD0name="D0toKpiCuts", TString cutsTrkname="AssociatedTrkCuts", TString effNa
   massD0Task->SetFillOnlyD0D0bar(flagD0D0bar);
   massD0Task->SetSystem(system); //0=pp, 1=PbPb
   massD0Task->SetEtaForCorrel(etacorr);
+  massD0Task->SetSoftPiFlag(flagsoftpicut);
+  massD0Task->SetMEAxisThresh(MEthresh);
 
 //*********************
 //correlation settings
@@ -185,22 +204,6 @@ cutsD0name="D0toKpiCuts", TString cutsTrkname="AssociatedTrkCuts", TString effNa
   massD0Task->SetPtTreshLow(pttreshlow);
   massD0Task->SetPtTreshUp(pttreshup);
 
-//********************
-//D0 efficiency
-//********************
-
-  TFile* fileeffD0;
-  fileeffD0=TFile::Open(effD0name.Data());
-  if(!fileeffD0->IsOpen()){
-    cout<<"Input file not found for efficiency! Exiting..."<<endl;
-    return;
-  }
-  TCanvas *c1 = (TCanvas*)fileeffD0->Get("c1");
-  TH1D *hEff = (TH1D*)c1->FindObject("h_Eff");
-  Double_t effD0values[15] = {1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.}; //only 4 to 11 bins are filled with efficiencies for now!
-  for(int i=3; i<11; i++) {effD0values[i]=hEff->GetBinContent(i+1);}
-  massD0Task->SetD0Efficiency(effD0values);
-
   //  massD0Task->SetRejectSDDClusters(kTRUE);
   //  massD0Task->SetWriteVariableTree(kTRUE);