]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/vertexingHF/charmFlow/Extractv2from2Dhistos.C
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / charmFlow / Extractv2from2Dhistos.C
index 8680819be2dd1069143fe7f274890ebb914af9af..650762c19aee724a1dae300fdb298e5011a72d00 100644 (file)
 #include <TGraphErrors.h>
 #include <TGraphAsymmErrors.h>
 
-
 #include "AliAODRecoDecayHF.h"
 #include "AliRDHFCuts.h"
 #include "AliRDHFCutsDplustoKpipi.h"
+#include "AliRDHFCutsDStartoKpipi.h"
 #include "AliRDHFCutsD0toKpi.h"
 #include "AliHFMassFitter.h"
-#include "AliEventPlaneResolution.h"
+#include "AliEventPlaneResolutionHandler.h"
+#include "AliVertexingHFUtils.h"
 #endif
 
 
 // Common variables: to be configured by the user
-TString filname="AnalysisResultsgoodruns.root";
-TString listname="coutputv2D0Std";//"coutputv2D0RAACuts"; //"coutputv2D0Std";
-TString filcutsname="AnalysisResultsgoodruns.root";
+TString filename="AnalysisResults_train63.root";
+TString suffix="v2Dplus3050Cut4upcutPIDTPC";
+TString partname="Dplus";
 Int_t minCent=30;
 Int_t maxCent=50;
-Int_t mesonPDG=421;
-
-const Int_t nFinalPtBins=3;
-Double_t ptlims[nFinalPtBins+1]={2.,5.,8.,12.};
+//evPlane flag from AliEventPlaneResolutionHandler: 
+//kTPCFullEta, kTPCPosEta,kVZERO,kVZEROA,kVZEROC
+Int_t evPlane=AliEventPlaneResolutionHandler::kTPCPosEta;
+//resolution flag fromAliEventPlaneResolutionHandler:
+//kTwoRandSub,kTwoChargeSub,kTwoEtaSub,kThreeSub,kThreeSubTPCGap
+Int_t evPlaneRes=AliEventPlaneResolutionHandler::kTwoEtaSub;
+Bool_t useNcollWeight=kFALSE;
+
+const Int_t nFinalPtBins=4;
+Double_t ptlims[nFinalPtBins+1]={3.,4.,6.,8.,12.};
 Double_t sigmaRangeForSig=2.5;
 Double_t sigmaRangeForBkg=4.5;
-Int_t rebinHistoSideBands[nFinalPtBins]={2,2,2};
+Int_t rebinHistoSideBands[nFinalPtBins]={2,2,2,2};
 Bool_t useConstantvsBkgVsMass=kFALSE;
-Int_t rebinHistov2Mass[nFinalPtBins]={2,2,5};
+Int_t rebinHistov2Mass[nFinalPtBins]={2,2,2,2};
 Int_t factor4refl=0;
-Int_t minPtBin[nFinalPtBins]={-1,-1,-1};
-Int_t maxPtBin[nFinalPtBins]={-1,-1,-1};
+Int_t minPtBin[nFinalPtBins]={-1,-1,-1,-1};
+Int_t maxPtBin[nFinalPtBins]={-1,-1,-1,-1};
+Bool_t saveAllCanvas=kFALSE;
+
+// systematic errors for D+, 30-50 TPC eventplane
+Double_t systErrMeth1[nFinalPtBins]={
+  (0.247359-0.152503)/2.,
+  (0.337073-0.290978)/2.,
+  (0.286888-0.222045)/2.,
+  (-0.078888+0.178259)/2.
+};
+
+Double_t systErrMeth2[nFinalPtBins]={
+  (0.151019-0.066624)/2.,
+  (0.272252-0.232044)/2.,
+  (0.246019-0.143002)/2.,
+  (-0.026077+0.292956)/2.
+};
 
 // systematic errors for 2-5, 5-8 and 8-12 (no run-by-run weights)
 /*
@@ -63,6 +86,7 @@ Double_t systErrMeth2[nFinalPtBins]={
 */
 
 // systematic errors for 2-5, 5-8 and 8-12 (93 runs with run-by-run weights)
+/*
 Double_t systErrMeth1[nFinalPtBins]={
   (0.23-0.10)/2.,
   (0.152-0.078)/2.,
@@ -72,7 +96,9 @@ Double_t systErrMeth2[nFinalPtBins]={
   (0.164-0.097)/2.,
   (0.110-0.012)/2.,
   (0.131-0.036)/2.
-};
+;
+*/
+
 
 
 // systematic errors for 2-5, 5-8 and 8-12 (93 runs with run-by-run weights, RAA cuts)
@@ -102,34 +128,56 @@ void LoadMassHistos(TList* lst, TH2F** hMassDphi);
 Bool_t DefinePtBins(TDirectoryFile* df);
 Double_t v2vsMass(Double_t *x, Double_t *par);
 Double_t GetEPResolution(TList* lst, Double_t &rcflow, Double_t &rcfhigh);
+Double_t ComputeEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3);
 Bool_t FitMassSpectrum(TH1F* hMass, TPad* pad);
 TH1F* DoSideBands(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin, TCanvas* c1, Int_t optBkg=0);
 TF1* DoFitv2vsMass(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin, TCanvas* c2, Int_t optErrors=0, Bool_t useConst=kTRUE);
 
 
+//______________________________________________________________________________
 void Extractv2from2Dhistos(){
+  // main function: computes v2 with side band and v2(M) methods
 
   gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
 
-  TFile* filcuts=new TFile(filcutsname.Data());
-  TDirectoryFile* dfcuts=(TDirectoryFile*)filcuts->Get("PWG3_D2H_HFv2");
-  Bool_t binOK=DefinePtBins(dfcuts);
+  TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
+  TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
+
+  TFile *f = TFile::Open(filename.Data());
+  if(!f){
+    printf("file %s not found, please check file name\n",filename.Data());
+    return;
+  }
+  TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
+  if(!dir){
+    printf("Directory %s not found, please check dir name\n",dirname.Data());
+    return;
+  }
+  Bool_t binOK=DefinePtBins(dir);
   if(!binOK){
     printf("ERROR: mismatch in pt binning\n");
     return;
   }
-
-  TFile* fil=new TFile(filname.Data());
-  TDirectoryFile* df=(TDirectoryFile*)fil->Get("PWG3_D2H_HFv2");
-  TList* lst=(TList*)df->Get(listname.Data());
+  TList *lst =(TList*)dir->Get(listname.Data());
   if(!lst){
-    printf("ERROR: list %s not found in file\n",listname.Data());
+    printf("list %s not found in file, please check list name\n",listname.Data());
     return;
   }
-  Double_t rcfmin,rcfmax;
-  Double_t resolFull=GetEPResolution(lst,rcfmin,rcfmax);
-  Double_t resolSyst=(rcfmax-rcfmin)/2./resolFull;
-  printf("Relative Systematic Error on RCF=%f\n",resolSyst);
+  // Double_t rcfmin,rcfmax;
+  // Double_t resolFull=GetEPResolution(lst,rcfmin,rcfmax);
+  // Double_t resolSyst=(rcfmax-rcfmin)/2./resolFull;
+  // printf("Relative Systematic Error on RCF=%f\n",resolSyst);
+  AliEventPlaneResolutionHandler* epres=new AliEventPlaneResolutionHandler();
+  epres->SetEventPlane(evPlane);
+  epres->SetResolutionOption(evPlaneRes);
+  epres->SetUseNcollWeights(useNcollWeight);
+  Double_t resolFull=epres->GetEventPlaneResolution(minCent,maxCent);
+  Double_t rcfmin=epres->GetEventPlaneResolution(minCent,minCent+2.5);
+  Double_t rcfmax=epres->GetEventPlaneResolution(maxCent-2.5,maxCent);
+  Double_t resolSyst=TMath::Abs(rcfmax-rcfmin)/2./resolFull;
+  delete epres;
+  printf("Event plane resolution %f\n",resolFull);
 
   TH2F** hMassDphi=new TH2F*[nFinalPtBins];
   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
@@ -177,8 +225,12 @@ void Extractv2from2Dhistos(){
     v2M2[iFinalPtBin]=v2obsM2/resolFull;
     errv2M2[iFinalPtBin]=errv2obsM2/resolFull;
     printf("v2 = %f +- %f\n",v2M2[iFinalPtBin],errv2M2[iFinalPtBin]);
-    c1[iFinalPtBin]->SaveAs(Form("cMeth1PtBin%d.root",iFinalPtBin));
-    c2[iFinalPtBin]->SaveAs(Form("cMeth2PtBin%d.root",iFinalPtBin));
+    if(saveAllCanvas){
+      c1[iFinalPtBin]->SaveAs(Form("cMeth1PtBin%d.root",iFinalPtBin));
+      c2[iFinalPtBin]->SaveAs(Form("cMeth2PtBin%d.root",iFinalPtBin));
+      c1[iFinalPtBin]->SaveAs(Form("cMeth1PtBin%d.eps",iFinalPtBin));
+      c2[iFinalPtBin]->SaveAs(Form("cMeth2PtBin%d.eps",iFinalPtBin));
+   }
   }
 
   printf("\n--------- Summary ------------\n");
@@ -239,9 +291,10 @@ void Extractv2from2Dhistos(){
    ent->SetTextColor(hv2m2->GetMarkerColor());
    leg2->Draw();
    cv2->Update();
-   cv2->SaveAs("Dzero-v2-2Dmethods.gif");
+   cv2->SaveAs(Form("v2%s-2Dmethods-%s.png",partname.Data(),suffix.Data()));
+   cv2->SaveAs(Form("v2%s-2Dmethods-%s.eps",partname.Data(),suffix.Data()));
 
-   TFile* outfil=new TFile("Dzero-v2-2Dmethods.root","recreate");
+   TFile* outfil=new TFile(Form("v2Output2DMeth_%s_%d_%d_%s.root",partname.Data(),minCent,maxCent,suffix.Data()),"recreate");
    outfil->cd();
    hv2m1->Write();
    hv2m2->Write();
@@ -251,6 +304,7 @@ void Extractv2from2Dhistos(){
 }
 
 
+//______________________________________________________________________________
 Double_t v2vsMass(Double_t *x, Double_t *par){
   // Fit function for signal+background
   // par[0] = S/B at the mass peak
@@ -264,6 +318,7 @@ Double_t v2vsMass(Double_t *x, Double_t *par){
   return par[3]*fracsig+par[4]*fracbkg;
 }
 
+//______________________________________________________________________________
 Double_t v2vsMassLin(Double_t *x, Double_t *par){
   // Fit function for signal+background
   // par[0] = S/B at the mass peak
@@ -278,6 +333,7 @@ Double_t v2vsMassLin(Double_t *x, Double_t *par){
   return par[3]*fracsig+(par[4]+par[5]*x[0])*fracbkg;
 }
 
+//______________________________________________________________________________
 Bool_t FitMassSpectrum(TH1F* hMass, TPad* pad){
 
   Int_t nMassBins=hMass->GetNbinsX();
@@ -285,7 +341,17 @@ Bool_t FitMassSpectrum(TH1F* hMass, TPad* pad){
   hMaxBin=nMassBins-2;
   Double_t hmin=hMass->GetBinLowEdge(hMinBin);
   Double_t hmax=hMass->GetBinLowEdge(hMaxBin)+hMass->GetBinWidth(hMaxBin);
-  Float_t massD=TDatabasePDG::Instance()->GetParticle(mesonPDG)->Mass();
+  Double_t massD;
+  if(partname.Contains("Dzero")) {
+    massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
+  }else if(partname.Contains("Dplus")){
+    massD=TDatabasePDG::Instance()->GetParticle(411)->Mass();
+  }else if(partname.Contains("Dstar")) {
+    massD=(TDatabasePDG::Instance()->GetParticle(413)->Mass() - TDatabasePDG::Instance()->GetParticle(421)->Mass()); 
+  }else{
+    printf("Wrong particle name\n");
+    massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();    
+  }
     
   AliHFMassFitter* fitter=new AliHFMassFitter(hMass,hmin,hmax,2,0,0);
   fitter->SetReflectionSigmaFactor(factor4refl);
@@ -307,6 +373,7 @@ Bool_t FitMassSpectrum(TH1F* hMass, TPad* pad){
   return kTRUE;
 }
 
+//______________________________________________________________________________
 TH1F* DoSideBands(Int_t iFinalPtBin,
                  TH2F* hMassDphi, 
                  TH1F* hMass, 
@@ -470,6 +537,7 @@ TH1F* DoSideBands(Int_t iFinalPtBin,
 }
 
 
+//______________________________________________________________________________
 TF1* DoFitv2vsMass(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin, TCanvas* c2, Int_t optErrors, Bool_t useConst){
 
   Int_t npars=fSB->GetNpar();
@@ -496,7 +564,6 @@ TF1* DoFitv2vsMass(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin,
   TH1F* hAverCos2Phi=new TH1F(Form("hAverCos2PhiBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
   TH1F* hFractionSig=new TH1F(Form("hFractionSigBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
   TH1F* hFractionBkg=new TH1F(Form("hFractionBkgBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
-
   for(Int_t iBin=1; iBin<=hMassDphi->GetNbinsY(); iBin++){
     TH1F* htemp=(TH1F*)hMassDphi->ProjectionX("htemp",iBin,iBin);
     hAverCos2Phi->SetBinContent(iBin,htemp->GetMean());
@@ -538,6 +605,7 @@ TF1* DoFitv2vsMass(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin,
   fv2->FixParameter(1,massSB);
   fv2->FixParameter(2,sigmaSB);
 
+  printf("LOOP %d\n",rebin);
   if((hAverCos2Phi->GetNbinsX()%rebin)==0){
     hAverCos2Phi->Rebin(rebin);
     hAverCos2Phi->Scale(1./(Double_t)rebin);
@@ -591,64 +659,90 @@ TF1* DoFitv2vsMass(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin,
 }
 
 
+//______________________________________________________________________________
 Double_t GetEPResolution(TList* lst, Double_t &rcflow, Double_t &rcfhigh){
   // Event plane resolution syst err (from wide centrality bin
   TH1F* hResolSubAB=0x0;
+  TH1F* hResolSubBC=0x0;
+  TH1F* hResolSubAC=0x0;
   Double_t xmin=1.;
   Double_t xmax=-1.;
   TGraphAsymmErrors* grSingle=new TGraphAsymmErrors(0);
   TGraphAsymmErrors* grInteg=new TGraphAsymmErrors(0);
   Int_t iPt=0;
-  for(Int_t iHisC=minCent; iHisC<=maxCent-5; iHisC+=5){    
-    TString hisnameEP=Form("hEvPlaneResocentr%d_%d",iHisC,iHisC+5);
-    TH1F* hResolSubABsing=(TH1F*)lst->FindObject(hisnameEP.Data());
-    Double_t resolFull=AliEventPlaneResolution::GetFullEvResol(hResolSubABsing,1);
-    Double_t resolFullmin=AliEventPlaneResolution::GetFullEvResolLowLim(hResolSubABsing,1);
-    Double_t resolFullmax=AliEventPlaneResolution::GetFullEvResolHighLim(hResolSubABsing,1);
-    grSingle->SetPoint(iPt,iHisC+2.5,resolFull);
-    grSingle->SetPointEXlow(iPt,2.5);
-    grSingle->SetPointEXhigh(iPt,2.5);
+  Int_t minCentTimesTen=minCent*10;
+  Int_t maxCentTimesTen=maxCent*10;  
+  Double_t binHalfWid=25./20.;
+  for(Int_t iHisC=minCentTimesTen; iHisC<=maxCentTimesTen-25; iHisC+=25){    
+    TString hisnameAB=Form("hEvPlaneResocentr%d_%d",iHisC,iHisC+25);
+    TString hisnameBC=Form("hEvPlaneReso2centr%d_%d",iHisC,iHisC+25);
+    TString hisnameAC=Form("hEvPlaneReso3centr%d_%d",iHisC,iHisC+25);
+    TH1F* hResolSubABsing=(TH1F*)lst->FindObject(hisnameAB.Data());
+    cout<<hResolSubABsing->GetName()<<endl;
+    TH1F* hResolSubBCsing=0x0;
+    TH1F* hResolSubACsing=0x0;
+    if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
+       evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){
+      hResolSubBCsing=(TH1F*)lst->FindObject(hisnameBC.Data());
+      cout<<hResolSubBCsing->GetName()<<endl;
+      hResolSubACsing=(TH1F*)lst->FindObject(hisnameAC.Data());
+      cout<<hResolSubACsing->GetName()<<endl;
+    }
+    Double_t error; 
+    Double_t resolFull=ComputeEventPlaneResolution(error,hResolSubABsing,hResolSubBCsing,hResolSubACsing);
+    Double_t resolFullmin=resolFull-error;
+    Double_t resolFullmax=resolFull+error;
+    
+    Double_t binCentr=(Double_t)iHisC/10.+binHalfWid;
+    grSingle->SetPoint(iPt,binCentr,resolFull);
+    grSingle->SetPointEXlow(iPt,binHalfWid);
+    grSingle->SetPointEXhigh(iPt,binHalfWid);
     grSingle->SetPointEYlow(iPt,resolFullmax-resolFull);
     grSingle->SetPointEYhigh(iPt,resolFull-resolFullmin);
     ++iPt;
     if(resolFullmin<xmin) xmin=resolFullmin;
     if(resolFullmax>xmax) xmax=resolFullmax;
-    if(iHisC==minCent){
+    if(iHisC==minCentTimesTen){
       hResolSubAB=(TH1F*)hResolSubABsing->Clone("hResolSubAB");
+      if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
+        evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){
+       hResolSubBC=(TH1F*)hResolSubBCsing->Clone("hResolSubAB");
+       hResolSubAC=(TH1F*)hResolSubACsing->Clone("hResolSubAB");
+      }
     }else{
       hResolSubAB->Add(hResolSubABsing);
+      if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
+        evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){
+       hResolSubBC->Add(hResolSubBCsing);
+       hResolSubAC->Add(hResolSubACsing);
+      }
     }
-    printf("Adding histogram %s\n",hisnameEP.Data());
+    printf("Adding histogram %s\n",hResolSubAB->GetName());
   }
   rcflow=xmin;
   rcfhigh=xmax;
 
-  Double_t resolSub=AliEventPlaneResolution::GetSubEvResol(hResolSubAB);
-  printf("Resolution on sub event  = %.4f\n",resolSub);
-  Double_t chisub=AliEventPlaneResolution::FindChi(resolSub,1);
-  printf("Chi Subevent             = %.4f\n",chisub);
-  Double_t chifull=chisub*TMath::Sqrt(2);
-  printf("Chi Full Event           = %.4f\n",chifull);
-  Double_t resolFull=AliEventPlaneResolution::GetFullEvResol(hResolSubAB,1);
-  Double_t resolFullmin=AliEventPlaneResolution::GetFullEvResolLowLim(hResolSubAB,1);
-  Double_t resolFullmax=AliEventPlaneResolution::GetFullEvResolHighLim(hResolSubAB,1);
-
-  AliEventPlaneResolution *resol=new AliEventPlaneResolution(1);
-  resol->SetSubEventHisto(hResolSubAB);  
-  Double_t resolFull2=resol->GetFullEvResol();
-  printf("Resolution on full event = %.4f %.4f\n",resolFull,resolFull2);
+
+  Double_t error; 
+  Double_t resolFull=ComputeEventPlaneResolution(error,hResolSubAB,hResolSubBC,hResolSubAC);
 
   grInteg->SetPoint(0,40.,resolFull);
   grInteg->SetPointEXlow(0,10);
   grInteg->SetPointEXhigh(0,10.);
-  grInteg->SetPointEYlow(0,resolFullmax-resolFull);
-  grInteg->SetPointEYhigh(0,resolFull-resolFullmin);
+  grInteg->SetPointEYlow(0,error);
+  grInteg->SetPointEYhigh(0,error);
   
   TCanvas* cEP=new TCanvas("cEP","EvPlaneResol");
   cEP->Divide(1,2);
   cEP->cd(1);
   hResolSubAB->Draw();
+  if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
+     evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){
+    hResolSubBC->SetLineColor(2);
+    hResolSubBC->Draw("same");
+    hResolSubAC->SetLineColor(4);
+    hResolSubAC->Draw("same");
+  }
   TLatex* tres=new TLatex(0.15,0.7,Form("Resolution on full event = %.4f\n",resolFull));
   tres->SetNDC();
   tres->Draw();
@@ -666,13 +760,60 @@ Double_t GetEPResolution(TList* lst, Double_t &rcflow, Double_t &rcfhigh){
   
 }
   
+//______________________________________________________________________________
+Double_t ComputeEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3){
+  Double_t resolFull;
+  if(evPlaneRes==AliEventPlaneResolutionHandler::kTwoRandSub ||
+     evPlaneRes==AliEventPlaneResolutionHandler::kTwoChargeSub){
+    resolFull=AliVertexingHFUtils::GetFullEvResol(hsubev1);
+    error = TMath::Abs(resolFull-AliVertexingHFUtils::GetFullEvResolLowLim(hsubev1));
+  }else if(evPlaneRes==AliEventPlaneResolutionHandler::kTwoEtaSub){
+    if(evPlane==AliEventPlaneResolutionHandler::kTPCFullEta){
+      resolFull=AliVertexingHFUtils::GetFullEvResol(hsubev1);
+      error = TMath::Abs(resolFull-AliVertexingHFUtils::GetFullEvResolLowLim(hsubev1));
+    }else if(evPlane==AliEventPlaneResolutionHandler::kTPCPosEta){
+      resolFull=AliVertexingHFUtils::GetSubEvResol(hsubev1);
+      error = TMath::Abs(resolFull-AliVertexingHFUtils::GetSubEvResolLowLim(hsubev1));      
+    }
+  }else{
+    Double_t resolSub[3];
+    Double_t errors[3];
+    TH1F* hevplresos[3];
+    hevplresos[0]=hsubev1;
+    hevplresos[1]=hsubev2;
+    hevplresos[2]=hsubev3;
+    for(Int_t ires=0;ires<3;ires++){
+      resolSub[ires]=hevplresos[ires]->GetMean();
+      errors[ires]=hevplresos[ires]->GetMeanError();
+    }
+    Double_t lowlim[3];for(Int_t ie=0;ie<3;ie++)lowlim[ie]=TMath::Abs(resolSub[ie]-errors[ie]);
+    if(evPlane==AliEventPlaneResolutionHandler::kVZEROC ||
+       evPlane==AliEventPlaneResolutionHandler::kVZERO){
+      resolFull=TMath::Sqrt(resolSub[1]*resolSub[2]/resolSub[0]);
+      error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[1]/lowlim[0]);
+    }
+    else if(evPlane==AliEventPlaneResolutionHandler::kVZEROA){
+      resolFull=TMath::Sqrt(resolSub[0]*resolSub[2]/resolSub[1]);
+      error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[0]/lowlim[1]);
+    }
+    else if(evPlane==AliEventPlaneResolutionHandler::kTPCFullEta ||
+           evPlane==AliEventPlaneResolutionHandler::kTPCPosEta){
+      resolFull=TMath::Sqrt(resolSub[0]*resolSub[1]/resolSub[2]);
+      error=resolFull-TMath::Sqrt(lowlim[0]*lowlim[1]/lowlim[2]);
+    }
+  }
+  return resolFull;
+}
 
+//______________________________________________________________________________
 void LoadMassHistos(TList* lst, TH2F** hMassDphi){
 
-  for(Int_t iHisC=minCent; iHisC<=maxCent-5; iHisC+=5){    
+  Int_t minCentTimesTen=minCent*10;
+  Int_t maxCentTimesTen=maxCent*10;  
+  for(Int_t iHisC=minCentTimesTen; iHisC<=maxCentTimesTen-25; iHisC+=25){    
     for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
       for(Int_t iPtBin=minPtBin[iFinalPtBin]; iPtBin<=maxPtBin[iFinalPtBin]; iPtBin++){
-       TString hisname=Form("hMc2phi_pt%dcentr%d_%d",iPtBin,iHisC,iHisC+5);
+       TString hisname=Form("hMc2deltaphi_pt%dcentr%d_%d",iPtBin,iHisC,iHisC+25);
        TH2F* htmp=(TH2F*)lst->FindObject(hisname.Data());
        if(hMassDphi[iFinalPtBin]==0x0){
          hMassDphi[iFinalPtBin]=(TH2F*)htmp->Clone(Form("hMassCos2DphiBin%d",iFinalPtBin));
@@ -685,10 +826,14 @@ void LoadMassHistos(TList* lst, TH2F** hMassDphi){
   }
 }
 
+//______________________________________________________________________________
 Bool_t DefinePtBins(TDirectoryFile* df){
-  AliRDHFCuts *d0cuts;//=(AliRDHFCutsD0toKpi*)df->Get("D0toKpiCuts");
-  if(mesonPDG==421)d0cuts=(AliRDHFCutsD0toKpi*)df->Get("D0toKpiCuts");
-  if(mesonPDG==411)d0cuts=(AliRDHFCutsDplustoKpipi*)df->Get("AnalysisCuts");
+
+  AliRDHFCuts *d0cuts;
+  if(partname.Contains("Dzero")) d0cuts=(AliRDHFCutsD0toKpi*)df->Get(df->GetListOfKeys()->At(2)->GetName());
+  if(partname.Contains("Dplus")) d0cuts=(AliRDHFCutsDplustoKpipi*)df->Get(df->GetListOfKeys()->At(2)->GetName());
+  if(partname.Contains("Dstar")) d0cuts=((AliRDHFCutsDStartoKpipi*)df->Get(df->GetListOfKeys()->At(2)->GetName()));
+
   Int_t nPtBinsCuts=d0cuts->GetNPtBins();
   printf("Number of pt bins for cut object = %d\n",nPtBinsCuts);
   Float_t *ptlimsCuts=d0cuts->GetPtBinLimits();
@@ -712,24 +857,42 @@ Bool_t DefinePtBins(TDirectoryFile* df){
 }
 
 
+//______________________________________________________________________________
 void SystForSideBands(){
+  // Compute systematic ucnertainty for the side band subtraction method
   gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
 
-  TFile* filcuts=new TFile(filcutsname.Data());
-  TDirectoryFile* dfcuts=(TDirectoryFile*)filcuts->Get("PWG3_D2H_HFv2");
-  Bool_t binOK=DefinePtBins(dfcuts);
-  if(!binOK) return;
+  TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
+  TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
 
-  TFile* fil=new TFile(filname.Data());
-  TDirectoryFile* df=(TDirectoryFile*)fil->Get("PWG3_D2H_HFv2");
-  TList* lst=(TList*)df->Get(listname.Data());
+  TFile *f = TFile::Open(filename.Data());
+  if(!f){
+    printf("file %s not found, please check file name\n",filename.Data());
+    return;
+  }
+  TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
+  if(!dir){
+    printf("Directory %s not found, please check dir name\n",dirname.Data());
+    return;
+  }
+  Bool_t binOK=DefinePtBins(dir);
+  if(!binOK){
+    printf("ERROR: mismatch in pt binning\n");
+    return;
+  }
+  TList *lst =(TList*)dir->Get(listname.Data());
   if(!lst){
-    printf("ERROR: list %s not found in file\n",listname.Data());
+    printf("list %s not found in file, please check list name\n",listname.Data());
     return;
   }
 
-  Double_t rcfmin,rcfmax;
-  Double_t resolFull=GetEPResolution(lst,rcfmin,rcfmax);
+  AliEventPlaneResolutionHandler* epres=new AliEventPlaneResolutionHandler();
+  epres->SetEventPlane(evPlane);
+  epres->SetResolutionOption(evPlaneRes);
+  epres->SetUseNcollWeights(useNcollWeight);
+  Double_t resolFull=epres->GetEventPlaneResolution(minCent,maxCent);
+  delete epres;
+  printf("Event plane resolution %f\n",resolFull);
   
   TH2F** hMassDphi=new TH2F*[nFinalPtBins];
   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
@@ -915,25 +1078,43 @@ void SystForSideBands(){
   }
 }
 
+//______________________________________________________________________________
 void SystForFitv2Mass(){
+  // Compute systematic ucnertainty for the v2(M) method
   gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
 
-  TFile* filcuts=new TFile(filcutsname.Data());
-  TDirectoryFile* dfcuts=(TDirectoryFile*)filcuts->Get("PWG3_D2H_HFv2");
-  Bool_t binOK=DefinePtBins(dfcuts);
-  if(!binOK) return;
+  TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
+  TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
 
-  TFile* fil=new TFile(filname.Data());
-  TDirectoryFile* df=(TDirectoryFile*)fil->Get("PWG3_D2H_HFv2");
-  TList* lst=(TList*)df->Get(listname.Data());
+  TFile *f = TFile::Open(filename.Data());
+  if(!f){
+    printf("file %s not found, please check file name\n",filename.Data());
+    return;
+  }
+  TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
+  if(!dir){
+    printf("Directory %s not found, please check dir name\n",dirname.Data());
+    return;
+  }
+  Bool_t binOK=DefinePtBins(dir);
+  if(!binOK){
+    printf("ERROR: mismatch in pt binning\n");
+    return;
+  }
+  TList *lst =(TList*)dir->Get(listname.Data());
   if(!lst){
-    printf("ERROR: list %s not found in file\n",listname.Data());
+    printf("list %s not found in file, please check list name\n",listname.Data());
     return;
   }
 
-  Double_t rcfmin,rcfmax;
-  Double_t resolFull=GetEPResolution(lst,rcfmin,rcfmax);
-  
+  AliEventPlaneResolutionHandler* epres=new AliEventPlaneResolutionHandler();
+  epres->SetEventPlane(evPlane);
+  epres->SetResolutionOption(evPlaneRes);
+  epres->SetUseNcollWeights(useNcollWeight);
+  Double_t resolFull=epres->GetEventPlaneResolution(minCent,maxCent);
+  delete epres;
+  printf("Event plane resolution %f\n",resolFull);
+
   TH2F** hMassDphi=new TH2F*[nFinalPtBins];
   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
     hMassDphi[iFinalPtBin]=0x0;