PWGHF/vertexingHF/charmFlow/DmesonsFlowAnalysis.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <Riostream.h>
   3 #include <TFile.h>
   4 #include <TString.h>
   5 #include <TH2F.h>
   6 #include <TH1F.h>
   7 #include <TF1.h>
   8 #include <TGraph.h>
   9 #include <TGraphErrors.h>
  10 #include <TMultiGraph.h>
  11 #include <TDirectoryFile.h>
  12 #include <TList.h>
  13 #include <TCanvas.h>
  14 #include <TLegend.h>
  15 #include <TLegendEntry.h>
  16 #include <TPaveText.h>
  17 #include <TStyle.h>
  18 #include <TASImage.h>
  19 #include <TPad.h>
  20 #include <TROOT.h>
  21 #include <TDatabasePDG.h>
  22 #include <TParameter.h>
  23 #include <TLatex.h>
  24 #include <TGraphAsymmErrors.h>
  25
  26 #include <AliHFMassFitter.h>
  27 #include "AliRDHFCutsD0toKpi.h"
  28 #include "AliVertexingHFUtils.h"
  29 #include "AliRDHFCutsDplustoKpipi.h"
  30 #include "AliRDHFCutsDStartoKpipi.h"
  31 #include "AliEventPlaneResolutionHandler.h"
  32 #include <AliVertexingHFUtils.h>
  33
  34 #endif
  35
  36 /* $Id$ */
  37
  38 //methods for the analysis of AliAnalysisTaskSEHFv2 output
  39 //Authors: Chiara Bianchin cbianchi@pd.infn.it
  40 //         Giacomo Ortona  ortona@to.infn.it
  41 //         Francesco Prino prino@to.infn.it
  42
  43 //global variables to be set
  44 TString filename="AnalysisResults_train63.root";
  45 TString suffix="v2Dplus3050Cut4upcutPIDTPC";
  46 TString partname="Dplus";
  47 Int_t minCent=30;
  48 Int_t maxCent=50;
  49 //evPlane flag from AliEventPlaneResolutionHandler:
  50 //kTPCFullEta, kTPCPosEta,kVZERO,kVZEROA,kVZEROC
  51 Int_t evPlane=AliEventPlaneResolutionHandler::kTPCPosEta;
  52 //resolution flag fromAliEventPlaneResolutionHandler:
  53 //kTwoRandSub,kTwoChargeSub,kTwoEtaSub,kThreeSub,kThreeSubTPCGap
  54 Int_t evPlaneRes=AliEventPlaneResolutionHandler::kTwoEtaSub;
  55 Bool_t useNcollWeight=kFALSE;
  56 // pt and phi binning
  57 const Int_t nptbinsnew=4;
  58 Float_t ptbinsnew[nptbinsnew+1]={3.,4.,6.,8.,12.};
  59 const Int_t nphibins=4;
  60 Float_t phibinslim[nphibins+1]={0,TMath::Pi()/4,TMath::Pi()/2,3*TMath::Pi()/4,TMath::Pi()};
  61 // mass fit configuration
  62 Int_t rebin[nptbinsnew]={4,4,4,4};
  63 Int_t typeb=0;  // Background: 0=expo, 1=linear, 2=pol2
  64 Bool_t fixAlsoMass=kFALSE;
  65 Double_t minMassForFit=1.6;
  66 Double_t maxMassForFit=2.2;
  67 Float_t massRangeForCounting=0.1; // GeV
  68
  69 Bool_t gnopng=kFALSE; //don't save in png format (only root and eps)
  70 Int_t minPtBin[nptbinsnew]={-1,-1,-1,-1};
  71 Int_t maxPtBin[nptbinsnew]={-1,-1,-1,-1};
  72 Double_t effInOverEffOut=1.03;
  73 Double_t massD;
  74
  75 //methods
  76 TList *LoadMassHistos(TList *inputlist,Bool_t inoutanis);
  77 TList* LoadResolutionHistos(TList *inputlist);
  78 Int_t FindPtBin(Int_t nbins, Float_t* array,Float_t value);
  79 void FillSignalGraph(TList *histlist,TGraphAsymmErrors **gSignal,TGraphAsymmErrors **gSignalfs, TGraphAsymmErrors **gSignalBC1, TGraphAsymmErrors **gSignalBC2, Bool_t inoutanis);
  80 void DmesonsFlowAnalysis(Bool_t inoutanis=kTRUE);
  81 void DrawEventPlane();
  82 Double_t GetEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3);
  83
  84
  85
  86 //method implementation
  87 //_________________________________________________________________
  88 Double_t GetEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3){
  89   Double_t resolFull=1.;
  90   if(evPlaneRes==AliEventPlaneResolutionHandler::kTwoRandSub ||
  91      evPlaneRes==AliEventPlaneResolutionHandler::kTwoChargeSub){
  92     resolFull=AliVertexingHFUtils::GetFullEvResol(hsubev1);
  93     error = TMath::Abs(resolFull-AliVertexingHFUtils::GetFullEvResolLowLim(hsubev1));
  94   }else if(evPlaneRes==AliEventPlaneResolutionHandler::kTwoEtaSub){
  95     if(evPlane==AliEventPlaneResolutionHandler::kTPCFullEta){
  96       resolFull=AliVertexingHFUtils::GetFullEvResol(hsubev1);
  97       error = TMath::Abs(resolFull-AliVertexingHFUtils::GetFullEvResolLowLim(hsubev1));
  98     }else if(evPlane==AliEventPlaneResolutionHandler::kTPCPosEta){
  99       resolFull=AliVertexingHFUtils::GetSubEvResol(hsubev1);
 100       error = TMath::Abs(resolFull-AliVertexingHFUtils::GetSubEvResolLowLim(hsubev1));
 101     }
 102   }else{
 103     Double_t resolSub[3];
 104     Double_t errors[3];
 105     TH1F* hevplresos[3];
 106     hevplresos[0]=hsubev1;
 107     hevplresos[1]=hsubev2;
 108     hevplresos[2]=hsubev3;
 109     for(Int_t ires=0;ires<3;ires++){
 110       resolSub[ires]=hevplresos[ires]->GetMean();
 111       errors[ires]=hevplresos[ires]->GetMeanError();
 112     }
 113     Double_t lowlim[3];for(Int_t ie=0;ie<3;ie++)lowlim[ie]=TMath::Abs(resolSub[ie]-errors[ie]);
 114     if(evPlane==AliEventPlaneResolutionHandler::kVZEROC ||
 115        evPlane==AliEventPlaneResolutionHandler::kVZERO){
 116       resolFull=TMath::Sqrt(resolSub[1]*resolSub[2]/resolSub[0]);
 117       error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[1]/lowlim[0]);
 118     }
 119     else if(evPlane==AliEventPlaneResolutionHandler::kVZEROA){
 120       resolFull=TMath::Sqrt(resolSub[0]*resolSub[2]/resolSub[1]);
 121       error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[0]/lowlim[1]);
 122     }
 123     else if(evPlane==AliEventPlaneResolutionHandler::kTPCFullEta ||
 124             evPlane==AliEventPlaneResolutionHandler::kTPCPosEta){
 125       resolFull=TMath::Sqrt(resolSub[0]*resolSub[1]/resolSub[2]);
 126       error=resolFull-TMath::Sqrt(lowlim[0]*lowlim[1]/lowlim[2]);
 127     }
 128   }
 129   return resolFull;
 130 }
 131
 132 //____________________________________________________________________
 133 TList* LoadResolutionHistos(TList *inputlist){
 134
 135   TList *outlist = new TList();
 136   outlist->SetName("eventplanehistlist");
 137
 138   const Int_t nBins=20;
 139   Double_t ncoll[nBins]={1790.77,1578.44,1394.82,1236.17,1095.08,
 140                          969.86,859.571,759.959,669.648,589.588,
 141                          516.039,451.409,392.853,340.493,294.426,
 142                          252.385,215.484,183.284,155.101,130.963};
 143
 144   Int_t minCentTimesTen=minCent*10;
 145   Int_t maxCentTimesTen=maxCent*10;
 146   TGraphErrors* gResolVsCent=new TGraphErrors(0);
 147   Int_t iPt=0;
 148   Int_t nSubRes=1;
 149   if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
 150      evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap)nSubRes=3;
 151   TString namereso[3]={"Reso","Reso2","Reso3"};
 152   TString suffixcentr=Form("centr%d_%d",minCent,maxCent);
 153   TH2F* hevpls=(TH2F*)inputlist->FindObject(Form("hEvPlanecentr%d_%d",minCentTimesTen,minCentTimesTen+25));
 154   hevpls->SetName(Form("hEvPlane%s",suffixcentr.Data()));
 155   hevpls->SetTitle(Form("Event Plane angle %s",suffixcentr.Data()));
 156     //    TH1F* hevplresos=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",mincentr,mincentr+5));QUI!!!
 157   TH1F* hevplresos[3];
 158   Int_t ncBin=minCentTimesTen/25;
 159
 160   for(Int_t ires=0;ires<nSubRes;ires++){
 161     hevplresos[ires]=(TH1F*)inputlist->FindObject(Form("hEvPlane%scentr%d_%d",namereso[ires].Data(),minCentTimesTen,minCentTimesTen+25));
 162     if(hevplresos[ires]){
 163       hevplresos[ires]->SetName(Form("hEvPlane%s%s",namereso[ires].Data(),suffixcentr.Data()));
 164       hevplresos[ires]->SetTitle(Form("Event Plane Resolution %s%s",namereso[ires].Data(),suffixcentr.Data()));
 165       if(useNcollWeight){
 166         printf("Centr %d Bin %d  Ncoll %f\n",minCentTimesTen,ncBin,ncoll[ncBin]);
 167         hevplresos[ires]->Scale(ncoll[ncBin]);
 168       }
 169     }
 170   }
 171   Double_t error;
 172   Double_t lowestRes=1;
 173   Double_t highestRes=0;
 174   Double_t resolBin=GetEventPlaneResolution(error,hevplresos[0],hevplresos[1],hevplresos[2]);
 175   if(resolBin<lowestRes) lowestRes=resolBin;
 176   if(resolBin>highestRes) highestRes=resolBin;
 177
 178   Double_t binHalfWid=25./20.;
 179   Double_t binCentr=(Double_t)minCentTimesTen/10.+binHalfWid;
 180   gResolVsCent->SetPoint(iPt,binCentr,resolBin);
 181   gResolVsCent->SetPointError(iPt,binHalfWid,error);
 182   ++iPt;
 183
 184   for(Int_t icentr=minCentTimesTen+25;icentr<maxCentTimesTen;icentr=icentr+25){
 185     TH2F* h=(TH2F*)inputlist->FindObject(Form("hEvPlanecentr%d_%d",icentr,icentr+25));
 186     if(h)hevpls->Add(h);
 187     else cout<<"skipping ev plane "<<icentr<<"_"<<icentr+5<<endl;
 188     TH1F* htmpresos[3];
 189     for(Int_t ires=0;ires<nSubRes;ires++){
 190       htmpresos[ires]=(TH1F*)inputlist->FindObject(Form("hEvPlane%scentr%d_%d",namereso[ires].Data(),icentr,icentr+25));
 191       if(!htmpresos[ires])cout<<"skipping ev pl reso "<<icentr<<"_"<<icentr+25<<endl;
 192     }
 193     resolBin=GetEventPlaneResolution(error,htmpresos[0],htmpresos[1],htmpresos[2]);
 194     if(resolBin<lowestRes) lowestRes=resolBin;
 195     if(resolBin>highestRes) highestRes=resolBin;
 196     binCentr=(Double_t)icentr/10.+binHalfWid;
 197     gResolVsCent->SetPoint(iPt,binCentr,resolBin);
 198     gResolVsCent->SetPointError(iPt,binHalfWid,error);
 199     ++iPt;
 200     ncBin=icentr/25;
 201     for(Int_t ires=0;ires<nSubRes;ires++){
 202       if(htmpresos[ires]){
 203         if(useNcollWeight){
 204           printf("Centr %d Bin %d  Ncoll %f\n",icentr,ncBin,ncoll[ncBin]);
 205           htmpresos[ires]->Scale(ncoll[ncBin]);
 206         }
 207         hevplresos[ires]->Add(htmpresos[ires]);
 208       }
 209     }
 210   }
 211   outlist->Add(hevpls->Clone());
 212   for(Int_t ires=0;ires<nSubRes;ires++){
 213     if(hevplresos[ires]) outlist->Add(hevplresos[ires]->Clone());
 214   }
 215   gResolVsCent->SetName("gResolVsCent");
 216   gResolVsCent->SetTitle("Resolution vs. Centrality");
 217   outlist->Add(gResolVsCent->Clone());
 218   return outlist;
 219 }
 220
 221 //__________________________________________________________
 222 TList *LoadMassHistos(TList *inputlist,Bool_t inoutanis){
 223   // printf("Start load histos\n");
 224   //  const Int_t nptbins=cutobj->GetNPtBins();
 225   TList *outlist = new TList();
 226   outlist->SetName("azimuthalhistoslist");
 227
 228   Int_t nphi=nphibins;
 229   if(inoutanis)nphi=2;
 230   Int_t minCentTimesTen=minCent*10;
 231   Int_t maxCentTimesTen=maxCent*10;
 232
 233   //Create 2D histogram in final pt bins
 234   for(Int_t iFinalPtBin=0; iFinalPtBin<nptbinsnew; iFinalPtBin++){
 235     for(Int_t iphi=0;iphi<nphi;iphi++){
 236       TH1F *hMass=0x0;//=new TH1F();
 237       for(Int_t iPtBin=minPtBin[iFinalPtBin]; iPtBin<=maxPtBin[iFinalPtBin];iPtBin++){
 238         for(Int_t iHisC=minCentTimesTen; iHisC<=maxCentTimesTen-25; iHisC+=25){
 239           TString hisname=Form("hMdeltaphi_pt%dcentr%d_%d",iPtBin,iHisC,iHisC+25);
 240           TH2F* htmp=(TH2F*)inputlist->FindObject(hisname.Data());
 241           printf("---> Histogram: %s\n",htmp->GetName());
 242           Int_t startX=htmp->FindBin(phibinslim[iphi]);
 243           Int_t endX=htmp->FindBin(phibinslim[iphi+1]-0.0001); // -0.0001 to be sure that the upper limit of the bin is properly set
 244           TH1F *h1tmp;
 245           if(inoutanis){
 246             if(iphi==0){
 247               Int_t firstBin=htmp->FindBin(0);
 248               Int_t lastBin=htmp->FindBin(TMath::Pi()/4.-0.0001); // -0.0001 to be sure that the upper limit of the bin is pi/4
 249               h1tmp=(TH1F*)htmp->ProjectionY(Form("hMass%d_phi0",iPtBin),firstBin,lastBin);
 250               printf("In-plane, Range: bins %d-%d -> phi %f - %f\n",firstBin,lastBin,htmp->GetXaxis()->GetBinLowEdge(firstBin),htmp->GetXaxis()->GetBinUpEdge(lastBin));
 251               firstBin=htmp->FindBin(3.*TMath::Pi()/4.);
 252               lastBin=htmp->FindBin(TMath::Pi()-0.0001); // -0.0001 to be sure that the upper limit of the bin is pi
 253               h1tmp->Add((TH1F*)htmp->ProjectionY(Form("hMass%d",iPtBin),firstBin,lastBin));
 254               printf("In-plane, Range: bins %d-%d -> phi %f - %f\n",firstBin,lastBin,htmp->GetXaxis()->GetBinLowEdge(firstBin),htmp->GetXaxis()->GetBinUpEdge(lastBin));
 255             }else{
 256               Int_t firstBin=htmp->FindBin(TMath::Pi()/4.);
 257               Int_t lastBin=htmp->FindBin(3.*TMath::Pi()/4.-0.0001); // -0.0001 to be sure that the upper limit of the bin is pi/4
 258               h1tmp=(TH1F*)htmp->ProjectionY(Form("hMass%d_phi1",iPtBin),firstBin,lastBin);
 259               printf("Out-of-plane, Range: bins %d-%d -> phi %f - %f\n",firstBin,lastBin,htmp->GetXaxis()->GetBinLowEdge(firstBin),htmp->GetXaxis()->GetBinUpEdge(lastBin));
 260             }
 261           }else{
 262             h1tmp=(TH1F*)htmp->ProjectionY(Form("hMass%d_phi%d",iPtBin,iphi),startX,endX);
 263           }
 264           if(hMass==0)hMass=(TH1F*)h1tmp->Clone();
 265           else hMass->Add((TH1F*)h1tmp->Clone());
 266         }
 267       }
 268       hMass->SetTitle(Form("hMass_pt%d_phi%d",iFinalPtBin,iphi));
 269       hMass->SetName(Form("hMass_pt%d_phi%d",iFinalPtBin,iphi));
 270       outlist->Add(hMass->Clone());
 271       //      hMass->DrawClone();
 272       delete hMass;
 273       hMass=0x0;
 274     }
 275   }
 276   return outlist;
 277 }
 278 //______________________________________________________________
 279 Bool_t DefinePtBins(AliRDHFCuts *cutobj){
 280   Int_t nPtBinsCuts=cutobj->GetNPtBins();
 281   Float_t *ptlimsCuts=cutobj->GetPtBinLimits();
 282   //  for(Int_t iPt=0; iPt<nPtBinsCuts; iPt++) printf(" %d %f-%f\n",iPt,ptlimsCuts[iPt],ptlimsCuts[iPt+1]);
 283   for(Int_t iPtCuts=0; iPtCuts<nPtBinsCuts; iPtCuts++){
 284     for(Int_t iFinalPtBin=0; iFinalPtBin<nptbinsnew; iFinalPtBin++){
 285       if(TMath::Abs(ptlimsCuts[iPtCuts]-ptbinsnew[iFinalPtBin])<0.0001){
 286         minPtBin[iFinalPtBin]=iPtCuts;
 287         if(iFinalPtBin>0) maxPtBin[iFinalPtBin-1]=iPtCuts-1;
 288       }
 289     }
 290     if(TMath::Abs(ptlimsCuts[iPtCuts]-ptbinsnew[nptbinsnew])<0.0001) maxPtBin[nptbinsnew-1]=iPtCuts-1;
 291   }
 292   if(TMath::Abs(ptbinsnew[nptbinsnew]-ptlimsCuts[nPtBinsCuts])<0.0001) maxPtBin[nptbinsnew-1]=nPtBinsCuts-1;
 293   for(Int_t iFinalPtBin=0; iFinalPtBin<nptbinsnew; iFinalPtBin++){
 294     printf("Pt bins to be merged: %d %d\n",minPtBin[iFinalPtBin],maxPtBin[iFinalPtBin]);
 295     if(minPtBin[iFinalPtBin]<0 || maxPtBin[iFinalPtBin]<0) return kFALSE;
 296   }
 297
 298   return kTRUE;
 299 }
 300 //______________________________________________________________
 301 Int_t GetPadNumber(Int_t ix,Int_t iy){
 302   return (iy)*nptbinsnew+ix+1;
 303 }
 304 //______________________________________________________________
 305 void FillSignalGraph(TList *histlist,TGraphAsymmErrors **gSignal,TGraphAsymmErrors **gSignalfs, TGraphAsymmErrors **gSignalBC1, TGraphAsymmErrors **gSignalBC2, Bool_t inoutanis){
 306
 307   Int_t nphi=nphibins;
 308   if(inoutanis)nphi=2;
 309
 310   //Canvases for drawing histograms
 311   TCanvas *cDeltaPhi = new TCanvas("cinvmassdeltaphi","Invariant mass distributions",1200,700);
 312   TCanvas *cDeltaPhifs = new TCanvas("cinvmassdeltaphifs","Invariant mass distributions - fit with fixed sigma",1200,700);
 313   TCanvas *cPhiInteg = new TCanvas("cinvmass","Invariant mass distributions - #phi integrated",1200,350);
 314   cDeltaPhi->Divide(nptbinsnew,nphi);
 315   cDeltaPhifs->Divide(nptbinsnew,nphi);
 316   cPhiInteg->Divide(nptbinsnew,1);
 317   Int_t nMassBins;
 318   Double_t hmin,hmax;
 319   for(Int_t ipt=0;ipt<nptbinsnew;ipt++){
 320     TH1F *histtofitfullsigma=(TH1F*)histlist->FindObject(Form("hMass_pt%d_phi0",ipt))->Clone();
 321     for(Int_t iphi=0;iphi<nphi;iphi++){
 322       Int_t ipad=GetPadNumber(ipt,iphi);
 323       Double_t signal=0,esignal=0;
 324       TH1F *histtofit=(TH1F*)histlist->FindObject(Form("hMass_pt%d_phi%d",ipt,iphi))->Clone();
 325       if(iphi>0)histtofitfullsigma->Add((TH1F*)histtofit->Clone());
 326       if(!histtofit){
 327         gSignal[ipt]->SetPoint(iphi,iphi,signal);
 328         gSignal[ipt]->SetPointError(iphi,0,0,esignal,esignal);
 329         return;
 330       }
 331       histtofit->SetTitle(Form("%.1f<p_{t}<%.1f, #phi%d",ptbinsnew[ipt],ptbinsnew[ipt+1],iphi));
 332       nMassBins=histtofit->GetNbinsX();
 333       hmin=TMath::Max(minMassForFit,histtofit->GetBinLowEdge(2));
 334       hmax=TMath::Min(maxMassForFit,histtofit->GetBinLowEdge(nMassBins-2));
 335       histtofit->Rebin(rebin[ipt]);
 336       AliHFMassFitter fitter(histtofit,hmin,hmax,1,typeb);
 337       fitter.SetInitialGaussianMean(massD);
 338       //      fitter.SetInitialGaussianSigma(0.012);
 339       Bool_t ok=fitter.MassFitter(kFALSE);
 340       if(ok){
 341         fitter.DrawHere(cDeltaPhi->cd(ipad),3,1);
 342         fitter.Signal(3,signal,esignal);
 343       }
 344       gSignal[ipt]->SetPoint(iphi,iphi,signal);
 345       gSignal[ipt]->SetPointError(iphi,0,0,esignal,esignal);
 346       TF1* fB1=fitter.GetBackgroundFullRangeFunc();
 347       TF1* fB2=fitter.GetBackgroundRecalcFunc();
 348       Float_t minBinSum=histtofit->FindBin(massD-massRangeForCounting);
 349       Float_t maxBinSum=histtofit->FindBin(massD+massRangeForCounting);
 350       Float_t cntSig1=0.;
 351       Float_t cntSig2=0.;
 352       Float_t cntErr=0.;
 353       for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
 354         Float_t bkg1=fB1 ? fB1->Eval(histtofit->GetBinCenter(iMB)) : 0;
 355         Float_t bkg2=fB2 ? fB2->Eval(histtofit->GetBinCenter(iMB)) : 0;
 356         //      printf("PtBin %d MassBin%d BC %f %f %f\n",ipt,iMB,histtofit->GetBinContent(iMB),bkg1,bkg2);
 357         cntSig1+=(histtofit->GetBinContent(iMB)-bkg1);
 358         cntSig2+=(histtofit->GetBinContent(iMB)-bkg2);
 359         cntErr+=(histtofit->GetBinContent(iMB));
 360       }
 361       cntErr=TMath::Sqrt(cntErr);
 362       gSignalBC1[ipt]->SetPoint(iphi,iphi,cntSig1);
 363       gSignalBC1[ipt]->SetPointError(iphi,0,0,cntErr,cntErr);
 364       gSignalBC2[ipt]->SetPoint(iphi,iphi,cntSig2);
 365       gSignalBC2[ipt]->SetPointError(iphi,0,0,cntErr,cntErr);
 366     }
 367     //fit for fixed sigma
 368     histtofitfullsigma->SetTitle(Form("%.1f<p_{t}<%.1f",ptbinsnew[ipt],ptbinsnew[ipt+1]));
 369     nMassBins=histtofitfullsigma->GetNbinsX();
 370     hmin=TMath::Max(minMassForFit,histtofitfullsigma->GetBinLowEdge(2));
 371     hmax=TMath::Min(maxMassForFit,histtofitfullsigma->GetBinLowEdge(nMassBins-2));
 372     histtofitfullsigma->Rebin(rebin[ipt]);
 373     AliHFMassFitter fitter(histtofitfullsigma,hmin,hmax,1,typeb);
 374     fitter.SetInitialGaussianMean(massD);
 375     //      fitter.SetInitialGaussianSigma(0.012);
 376     //fitter.SetFixGaussianSigma(0.012);
 377     Bool_t ok=fitter.MassFitter(kFALSE);
 378     if(ok){
 379       fitter.DrawHere(cPhiInteg->cd(ipt+1),3,1);
 380     }
 381     Double_t sigma=fitter.GetSigma();
 382     Double_t massFromFit=fitter.GetMean();
 383     for(Int_t iphi=0;iphi<nphi;iphi++){
 384       Int_t ipad=GetPadNumber(ipt,iphi);
 385       TH1F *histtofit=(TH1F*)histlist->FindObject(Form("hMass_pt%d_phi%d",ipt,iphi))->Clone();
 386       histtofit->SetTitle(Form("%.1f<p_{t}<%.1f, #phi%d",ptbinsnew[ipt],ptbinsnew[ipt+1],iphi));
 387       nMassBins=histtofit->GetNbinsX();
 388       hmin=TMath::Max(minMassForFit,histtofit->GetBinLowEdge(2));
 389       hmax=TMath::Min(maxMassForFit,histtofit->GetBinLowEdge(nMassBins-2));
 390       histtofit->Rebin(rebin[ipt]);
 391       AliHFMassFitter fitter2(histtofit,hmin,hmax,1,typeb);
 392       fitter2.SetInitialGaussianMean(massD);
 393       fitter2.SetFixGaussianSigma(sigma);
 394       if(fixAlsoMass) fitter2.SetFixGaussianMean(massFromFit);
 395       Bool_t ok2=fitter2.MassFitter(kFALSE);
 396       Double_t signal=0,esignal=0;
 397       if(ok2){
 398         fitter2.DrawHere(cDeltaPhifs->cd(ipad),3,1);
 399         fitter2.Signal(3,signal,esignal);
 400       }
 401       gSignalfs[ipt]->SetPoint(iphi,iphi,signal);
 402       gSignalfs[ipt]->SetPointError(iphi,0,0,esignal,esignal);
 403     }
 404   }//end loop on pt bin
 405
 406
 407   cDeltaPhi->SaveAs(Form("InvMassDeltaPhi_%s.eps",suffix.Data()));
 408   cDeltaPhifs->SaveAs(Form("InvMassDeltaPhi_fs_%s.eps",suffix.Data()));
 409   cDeltaPhifs->SaveAs(Form("InvMassDeltaPhi_fs_%s.root",suffix.Data()));
 410   cPhiInteg->SaveAs(Form("InvMassfullphi_%s.eps",suffix.Data()));
 411
 412   if(!gnopng){
 413     cDeltaPhi->SaveAs(Form("InvMassDeltaPhi_%s.png",suffix.Data()));
 414     cDeltaPhifs->SaveAs(Form("InvMassDeltaPhi_fs_%s.png",suffix.Data()));
 415     cPhiInteg->SaveAs(Form("InvMassfullphi_%s.png",suffix.Data()));
 416   }
 417   //cDeltaPhifs->DrawClone();
 418   //cDeltaPhifs->Close();
 419
 420 }
 421 //______________________________________________________________
 422 void DmesonsFlowAnalysis(Bool_t inoutanis){
 423   TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
 424   TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
 425
 426   AliRDHFCuts *cutsobj=0x0;
 427   //Load input data from AliAnalysisTaskSEHFv2
 428   TFile *f = TFile::Open(filename.Data());
 429   if(!f){
 430     printf("file %s not found, please check file name\n",filename.Data());return;
 431   }
 432   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
 433   if(!dir){
 434     printf("Directory %s not found, please check dir name\n",dirname.Data());return;
 435   }
 436   if(partname.Contains("Dzero")) {
 437     cutsobj=((AliRDHFCutsD0toKpi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
 438     massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 439   }
 440   if(partname.Contains("Dplus")){
 441     cutsobj=((AliRDHFCutsDplustoKpipi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
 442     massD=TDatabasePDG::Instance()->GetParticle(411)->Mass();
 443   }
 444   if(partname.Contains("Dstar")) {
 445     cutsobj=((AliRDHFCutsDStartoKpipi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
 446     massD=(TDatabasePDG::Instance()->GetParticle(413)->Mass() - TDatabasePDG::Instance()->GetParticle(421)->Mass());
 447   }
 448
 449   TList *list =(TList*)dir->Get(listname.Data());
 450   if(!list){
 451     printf("list %s not found in file, please check list name\n",listname.Data());return;
 452   }
 453   if(!cutsobj){
 454     printf("cut object not found in file, please check keylist number\n");return;
 455   }
 456   //Define new pt bins
 457   if(!DefinePtBins(cutsobj)){
 458     printf("cut not define pt bins\n");return;
 459   }
 460
 461   //Load mass histograms corresponding to the required centrality, pt range and phi binning
 462   printf("Load mass histos \n");
 463   TList *histlist=LoadMassHistos(list,inoutanis);
 464   TString aniss="";
 465   if(inoutanis)aniss+="anis";
 466   histlist->SaveAs(Form("v2Histograms_%d_%d_%s_%s.root",minCent,maxCent,aniss.Data(),suffix.Data()),"RECREATE");
 467
 468   Int_t nphi=nphibins;
 469   if(inoutanis)nphi=2;
 470
 471
 472   printf("average pt for pt bin \n");
 473   //average pt for pt bin
 474   AliVertexingHFUtils *utils=new AliVertexingHFUtils();
 475   Int_t minCentTimesTen=minCent*10;
 476   Int_t maxCentTimesTen=maxCent*10;
 477   TH2F* hmasspt=(TH2F*)list->FindObject(Form("hMPtCandcentr%d_%d",minCentTimesTen,minCentTimesTen+25));
 478   for(Int_t icent=minCentTimesTen+25;icent<maxCentTimesTen;icent=icent+25)hmasspt->Add((TH2F*)list->FindObject(Form("hMPtCandcentr%d_%d",icent,icent+25)));
 479   Float_t averagePt[nptbinsnew];
 480   Float_t errorPt[nptbinsnew];
 481   for(Int_t ipt=0;ipt<nptbinsnew;ipt++){
 482     Int_t binMin=hmasspt->FindBin(ptbinsnew[ipt]);
 483     Int_t binMax=hmasspt->FindBin(ptbinsnew[ipt+1]-0.001);
 484     if(TMath::Abs(hmasspt->GetXaxis()->GetBinLowEdge(binMin)-ptbinsnew[ipt])>0.001 ||
 485        TMath::Abs(hmasspt->GetXaxis()->GetBinUpEdge(binMax)-ptbinsnew[ipt+1])>0.001){
 486       printf("Error in pt bin limits for projection!\n");
 487       return;
 488     }
 489     TH1F *histtofit = (TH1F*)hmasspt->ProjectionY("_py",binMin,binMax);
 490     Int_t nMassBins=histtofit->GetNbinsX();
 491     Double_t hmin=histtofit->GetBinLowEdge(2); // need wide range for <pt>
 492     Double_t hmax=histtofit->GetBinLowEdge(nMassBins-2); // need wide range for <pt>
 493     AliHFMassFitter fitter(histtofit,hmin,hmax,1);
 494     fitter.MassFitter(kFALSE);
 495     Float_t massFromFit=fitter.GetMean();
 496     Float_t sigmaFromFit=fitter.GetSigma();
 497     TF1* funcB2=fitter.GetBackgroundRecalcFunc();
 498     utils->AveragePt(averagePt[ipt],errorPt[ipt],ptbinsnew[ipt],ptbinsnew[ipt+1],hmasspt,massFromFit,sigmaFromFit,funcB2,2.5,4.5,0.,3.,1);
 499   }
 500   printf("Average pt\n");
 501   for(Int_t ipt=0;ipt<nptbinsnew;ipt++) printf("%f +- %f\n",averagePt[ipt],errorPt[ipt]);
 502
 503   printf("Fill TGraphs for signal \n");
 504   //Fill TGraphs for signal
 505   TGraphAsymmErrors *gSignal[nptbinsnew];
 506   TGraphAsymmErrors *gSignalfs[nptbinsnew];
 507   TGraphAsymmErrors *gSignalBC1[nptbinsnew];
 508   TGraphAsymmErrors *gSignalBC2[nptbinsnew];
 509   for(Int_t i=0;i<nptbinsnew;i++){
 510     gSignal[i]=new TGraphAsymmErrors(nphi);
 511     gSignal[i]->SetName(Form("gasigpt%d",i));
 512     gSignal[i]->SetTitle(Form("Signal %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
 513     gSignal[i]->SetMarkerStyle(25);
 514     gSignalfs[i]=new TGraphAsymmErrors(nphi);
 515     gSignalfs[i]->SetName(Form("gasigfspt%d",i));
 516     gSignalfs[i]->SetTitle(Form("Signal (fixed sigma) %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
 517     gSignalfs[i]->SetMarkerStyle(21);
 518     gSignalBC1[i]=new TGraphAsymmErrors(nphi);
 519     gSignalBC1[i]->SetName(Form("gasigBC1pt%d",i));
 520     gSignalBC1[i]->SetTitle(Form("Signal (BC1) %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
 521     gSignalBC2[i]=new TGraphAsymmErrors(nphi);
 522     gSignalBC2[i]->SetName(Form("gasigBC2pt%d",i));
 523     gSignalBC2[i]->SetTitle(Form("Signal (BC2) %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
 524   }
 525   FillSignalGraph(histlist,gSignal,gSignalfs,gSignalBC1,gSignalBC2,inoutanis);
 526
 527   //EP resolution
 528   // TList* resolhist=LoadResolutionHistos(list);
 529   // TString suffixcentr=Form("centr%d_%d",minCent,maxCent);
 530   // TH1F* hevplresos[3];
 531   // TString namereso[3]={"Reso","Reso2","Reso3"};
 532   // Int_t nSubRes=1;
 533   // if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
 534   //    evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap)nSubRes=3;
 535   // for(Int_t ires=0;ires<nSubRes;ires++){
 536   //   hevplresos[ires]=(TH1F*)resolhist->FindObject(Form("hEvPlane%s%s",namereso[ires].Data(),suffixcentr.Data()));
 537   // }
 538   // Double_t errorres;
 539   // Double_t resol=GetEventPlaneResolution(errorres,hevplresos[0],hevplresos[1],hevplresos[2]);
 540   AliEventPlaneResolutionHandler* epres=new AliEventPlaneResolutionHandler();
 541   epres->SetEventPlane(evPlane);
 542   epres->SetResolutionOption(evPlaneRes);
 543   epres->SetUseNcollWeights(useNcollWeight);
 544   Double_t resol=epres->GetEventPlaneResolution(minCent,maxCent);
 545   delete epres;
 546   printf("Event plane resolution %f\n",resol);
 547
 548   printf("Compute v2 \n");
 549   //compute v2
 550   TCanvas *cv2fs =new TCanvas("v2_fs","v2 Fit with fixed sigma");
 551   TCanvas *cv2 =new TCanvas("v2","v2 - systematic on yield extraction");
 552   TGraphAsymmErrors *gv2=new TGraphAsymmErrors(nptbinsnew);
 553   TGraphAsymmErrors *gv2fs=new TGraphAsymmErrors(nptbinsnew);
 554   TGraphAsymmErrors *gv2BC1=new TGraphAsymmErrors(nptbinsnew);
 555   TGraphAsymmErrors *gv2BC2=new TGraphAsymmErrors(nptbinsnew);
 556   TGraphAsymmErrors *grelSystEff=new TGraphAsymmErrors(nptbinsnew);
 557
 558   gv2->SetName("gav2");
 559   gv2->SetTitle("v_{2}(p_{t})");
 560   gv2fs->SetName("gav2fs");
 561   gv2fs->SetTitle("v_{2}(p_{t}) (fit with fixed sigma)");
 562   gv2BC1->SetName("gav2BC1");
 563   gv2BC1->SetTitle("v_{2}(p_{t}) (bin counting 1)");
 564   gv2BC2->SetName("gav2BC2");
 565   gv2BC2->SetTitle("v_{2}(p_{t}) (bin counting 2)");
 566   grelSystEff->SetName("grelSystEff");
 567   grelSystEff->SetTitle("SystErrEff");
 568
 569   //Prepare output file
 570   TFile *fout=new TFile(Form("v2Output_%d_%d_%s_%s.root",minCent,maxCent,aniss.Data(),suffix.Data()),"RECREATE");
 571
 572   for(Int_t ipt=0;ipt<nptbinsnew;ipt++){
 573     fout->cd();
 574     gSignal[ipt]->Write();
 575     gSignalfs[ipt]->Write();
 576     gSignalBC1[ipt]->Write();
 577     gSignalBC2[ipt]->Write();
 578
 579     if(inoutanis){
 580       //v2 from in-out anisotropy
 581       Double_t *y,*yfs,*ybc1,*ybc2;
 582       y=gSignal[ipt]->GetY();
 583       yfs=gSignalfs[ipt]->GetY();
 584       ybc1=gSignalBC1[ipt]->GetY();
 585       ybc2=gSignalBC2[ipt]->GetY();
 586
 587       Double_t nIn=y[0];
 588       Double_t nOut=y[1];
 589       Double_t enIn=gSignal[ipt]->GetErrorY(0);
 590       Double_t enOut=gSignal[ipt]->GetErrorY(1);
 591       Double_t anis=(nIn-nOut)/(nIn+nOut);
 592       //      Double_t eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
 593       Double_t eAnis=2./((nIn+nOut)*(nIn+nOut))*TMath::Sqrt(nIn*nIn*enOut*enOut+nOut*nOut*enIn*enIn);
 594       Double_t v2=anis*TMath::Pi()/4./resol;
 595       Double_t ev2=eAnis*TMath::Pi()/4./resol;
 596       gv2->SetPoint(ipt,averagePt[ipt],v2);
 597       gv2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
 598
 599       nIn=yfs[0];
 600       nOut=yfs[1];
 601       enIn=gSignalfs[ipt]->GetErrorY(0);
 602       enOut=gSignalfs[ipt]->GetErrorY(1);
 603       anis=(nIn-nOut)/(nIn+nOut);
 604       //     eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
 605       eAnis=2./((nIn+nOut)*(nIn+nOut))*TMath::Sqrt(nIn*nIn*enOut*enOut+nOut*nOut*enIn*enIn);
 606       v2=anis*TMath::Pi()/4./resol;
 607       ev2=eAnis*TMath::Pi()/4./resol;
 608       gv2fs->SetPoint(ipt,averagePt[ipt],v2);
 609       gv2fs->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
 610       Double_t anis1=(nIn-nOut*effInOverEffOut)/(nIn+nOut*effInOverEffOut);
 611       Double_t anis2=(nIn*effInOverEffOut-nOut)/(nIn*effInOverEffOut+nOut);
 612       Double_t systEffUp=0.,systEffDown=0.;
 613       if(anis1>anis && anis1>anis2) systEffUp=anis1/anis;
 614       if(anis2>anis && anis2>anis1) systEffUp=anis2/anis;
 615       if(anis1<anis && anis1<anis2) systEffDown=anis1/anis;
 616       if(anis2<anis && anis2<anis1) systEffDown=anis2/anis;
 617       printf(" Bin %d <pt>=%.3f  v2=%f+-%f systEff=%f %f\n",ipt,averagePt[ipt],v2,ev2,systEffUp*v2,systEffDown*v2);
 618       grelSystEff->SetPoint(ipt,averagePt[ipt],v2);
 619       grelSystEff->SetPointError(ipt,0.4,0.4,v2*(1-systEffDown),v2*(systEffUp-1));
 620
 621       nIn=ybc1[0];
 622       nOut=ybc1[1];
 623       enIn=gSignalBC1[ipt]->GetErrorY(0);
 624       enOut=gSignalBC1[ipt]->GetErrorY(1);
 625       anis=(nIn-nOut)/(nIn+nOut);
 626       //     eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
 627       eAnis=2./((nIn+nOut)*(nIn+nOut))*TMath::Sqrt(nIn*nIn*enOut*enOut+nOut*nOut*enIn*enIn);
 628       v2=anis*TMath::Pi()/4./resol;
 629       ev2=eAnis*TMath::Pi()/4./resol;
 630       gv2BC1->SetPoint(ipt,averagePt[ipt],v2);
 631       gv2BC1->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
 632
 633       nIn=ybc2[0];
 634       nOut=ybc2[1];
 635       enIn=gSignalBC2[ipt]->GetErrorY(0);
 636       enOut=gSignalBC2[ipt]->GetErrorY(1);
 637       anis=(nIn-nOut)/(nIn+nOut);
 638       //     eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
 639       eAnis=2./((nIn+nOut)*(nIn+nOut))*TMath::Sqrt(nIn*nIn*enOut*enOut+nOut*nOut*enIn*enIn);
 640       v2=anis*TMath::Pi()/4./resol;
 641       ev2=eAnis*TMath::Pi()/4./resol;
 642       gv2BC2->SetPoint(ipt,averagePt[ipt],v2);
 643       gv2BC2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
 644
 645    }else{
 646       TF1 *flowFunc = new TF1("flow","[0]*(1.+2.*[1]*TMath::Cos(2.*x))");
 647       //v2 from fit to Deltaphi distribution
 648       gSignal[ipt]->Fit(flowFunc);
 649       Double_t v2 = flowFunc->GetParameter(1)/resol;
 650       Double_t ev2=flowFunc->GetParError(1)/resol;
 651       gv2->SetPoint(ipt,averagePt[ipt],v2);
 652       gv2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
 653
 654       gSignalfs[ipt]->Fit(flowFunc);
 655       v2 = flowFunc->GetParameter(1)/resol;
 656       ev2=flowFunc->GetParError(1)/resol;
 657       gv2fs->SetPoint(ipt,averagePt[ipt],v2);
 658       gv2fs->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
 659
 660       gSignalBC1[ipt]->Fit(flowFunc);
 661       v2 = flowFunc->GetParameter(1)/resol;
 662       ev2=flowFunc->GetParError(1)/resol;
 663       gv2BC1->SetPoint(ipt,averagePt[ipt],v2);
 664       gv2BC1->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
 665
 666       gSignalBC2[ipt]->Fit(flowFunc);
 667       v2 = flowFunc->GetParameter(1)/resol;
 668       ev2=flowFunc->GetParError(1)/resol;
 669       gv2BC2->SetPoint(ipt,averagePt[ipt],v2);
 670       gv2BC2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
 671     }
 672   }
 673
 674   gv2->SetMarkerStyle(22);
 675   gv2->SetMarkerColor(2);
 676   gv2->SetLineColor(2);
 677   gv2fs->SetMarkerStyle(20);
 678   gv2fs->SetMarkerColor(1);
 679   gv2fs->SetLineColor(1);
 680   gv2BC1->SetMarkerStyle(23);
 681   gv2BC1->SetMarkerColor(kGreen+1);
 682   gv2BC1->SetLineColor(kGreen+1);
 683   gv2BC2->SetMarkerStyle(26);
 684   gv2BC2->SetMarkerColor(kGreen+1);
 685   gv2BC2->SetLineColor(kGreen+1);
 686
 687   cv2->cd();
 688   gv2fs->Draw("AP");
 689   gv2fs->SetMinimum(-0.2);
 690   gv2fs->SetMaximum(0.6);
 691   gv2fs->GetYaxis()->SetTitle("v_{2}");
 692   gv2fs->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 693   gv2->Draw("PSAME");
 694   gv2BC1->Draw("PSAME");
 695   gv2BC2->Draw("PSAME");
 696   TLegend* leg=new TLegend(0.65,0.7,0.89,0.89);
 697   leg->SetFillStyle(0);
 698   leg->SetBorderSize(0);
 699   leg->SetTextFont(42);
 700   leg->AddEntry(gv2fs,"Fixed sigma","P")->SetTextColor(gv2fs->GetMarkerColor());
 701   leg->AddEntry(gv2,"Free sigma","P")->SetTextColor(gv2->GetMarkerColor());
 702   leg->AddEntry(gv2BC1,"Bin counting 1","P")->SetTextColor(gv2BC1->GetMarkerColor());
 703   leg->AddEntry(gv2BC2,"Bin counting 2","P")->SetTextColor(gv2BC2->GetMarkerColor());
 704   leg->Draw();
 705   cv2fs->cd();
 706   gv2fs->Draw("AP");
 707   gv2fs->GetYaxis()->SetTitle("v_{2}");
 708   gv2fs->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 709
 710   fout->cd();
 711   gv2->Write();
 712   gv2fs->Write();
 713   gv2BC1->Write();
 714   gv2BC2->Write();
 715   grelSystEff->Write();
 716 }
 717 //___________________________________________________________
 718 Int_t FindPtBin(Int_t nbins, Float_t* array,Float_t value){
 719   for (Int_t i=0;i<nbins;i++){
 720     if(value>=array[i] && value<array[i+1]){
 721       return i;
 722     }
 723   }
 724   cout<<value<< " out of range "<<array[0]<<", "<<array[nbins]<<endl;
 725   return -1;
 726 }
 727
 728 //__________________________________________________________
 729 void DrawEventPlane(){
 730   TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
 731   TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
 732   TFile *f = TFile::Open(filename.Data());
 733   if(!f){
 734     printf("file %s not found, please check file name\n",filename.Data());return;
 735   }
 736   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
 737   if(!dir){
 738     printf("Directory %s not found, please check dir name\n",dirname.Data());return;
 739   }
 740   TList *list =(TList*)dir->Get(listname.Data());
 741   if(!list){
 742     printf("list %s not found in file, please check list name\n",listname.Data());return;
 743   }
 744   if(evPlane==AliEventPlaneResolutionHandler::kTPCPosEta &&
 745      evPlaneRes==AliEventPlaneResolutionHandler::kTwoEtaSub){
 746     printf("Wrong setting of event plane resolution forced it to kTwoSub\n");
 747     evPlaneRes=AliEventPlaneResolutionHandler::kTwoRandSub;
 748   }
 749   TList* resolhist=LoadResolutionHistos(list);
 750   TGraphErrors* gResolVsCent=(TGraphErrors*)resolhist->FindObject("gResolVsCent");
 751   Double_t lowestRes=1;
 752   Double_t highestRes=0;
 753   for(Int_t ipt=0; ipt<gResolVsCent->GetN(); ipt++){
 754     Double_t x,resolBin;
 755     gResolVsCent->GetPoint(ipt,x,resolBin);
 756     if(resolBin<lowestRes) lowestRes=resolBin;
 757     if(resolBin>highestRes) highestRes=resolBin;
 758   }
 759   TString suffixcentr=Form("centr%d_%d",minCent,maxCent);
 760   TH2F* hevpls=(TH2F*)resolhist->FindObject(Form("hEvPlane%s",suffixcentr.Data()));
 761   TH1F* hevplresos[3];
 762   TString namereso[3]={"Reso","Reso2","Reso3"};
 763   Int_t nSubRes=1;
 764   if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
 765      evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap)nSubRes=3;
 766   for(Int_t ires=0;ires<nSubRes;ires++){
 767     hevplresos[ires]=(TH1F*)resolhist->FindObject(Form("hEvPlane%s%s",namereso[ires].Data(),suffixcentr.Data()));
 768   }
 769
 770   gStyle->SetCanvasColor(0);
 771   gStyle->SetTitleFillColor(0);
 772   gStyle->SetStatColor(0);
 773   gStyle->SetOptStat(0);
 774   //gStyle->SetPalette(1);
 775   gStyle->SetOptFit(1);
 776   gStyle->SetLabelFont(42,"XYZ");
 777   gStyle->SetTextFont(42);
 778   gStyle->SetStatFont(42);
 779   //  gStyle->SetTitleFont(42);
 780   gStyle->SetStatY(0.89);
 781   gStyle->SetStatX(0.89);
 782   //  gStyle->SetTitleAlign(22);
 783   gStyle->SetTitleFont(42,"xyzg");
 784
 785   TH1F* htpc = (TH1F*)hevpls->ProjectionX();
 786   TH1F* hv0 = (TH1F*)hevpls->ProjectionY();
 787   TH1F* hUsedPlane;
 788   if(evPlane==AliEventPlaneResolutionHandler::kVZERO ||
 789      evPlane==AliEventPlaneResolutionHandler::kVZEROA ||
 790      evPlane==AliEventPlaneResolutionHandler::kVZEROC) hUsedPlane=(TH1F*)hv0->Clone("hUsedPlane");
 791   else hUsedPlane=(TH1F*)htpc->Clone("hUsedPlane");
 792
 793   TCanvas* cvtotevpl=new TCanvas(Form("cvtotevpl%s",suffixcentr.Data()),Form("Ev plane %s",suffixcentr.Data()),1200,400);
 794   cvtotevpl->Divide(3,1);
 795   cvtotevpl->cd(1);
 796   hevpls->SetTitleFont(42);
 797   hevpls->Draw("COLZ");
 798   cvtotevpl->cd(2);
 799   htpc->Draw();
 800   htpc->Fit("pol0");
 801   cvtotevpl->cd(3);
 802   hv0->Draw();
 803   hv0->Fit("pol0");
 804
 805   TCanvas* cvfitevpl=new TCanvas(Form("cvditevpl%s",suffixcentr.Data()),Form("Fit to Ev plane %s",suffixcentr.Data()));
 806   cvfitevpl->cd();
 807   hUsedPlane->Draw();
 808   TF1* four2=new TF1("four2","[0]*(1+2*[1]*cos(2*x))",hUsedPlane->GetXaxis()->GetXmin(),hUsedPlane->GetXaxis()->GetXmax());
 809   TF1* four2s=new TF1("four2s","[0]*(1+2*[1]*sin(2*x))",hUsedPlane->GetXaxis()->GetXmin(),hUsedPlane->GetXaxis()->GetXmax());
 810   hUsedPlane->Fit(four2s);
 811   hUsedPlane->Fit(four2);
 812   four2->SetLineColor(1);
 813   four2s->SetLineColor(4);
 814   four2->Draw("same");
 815   four2s->Draw("same");
 816   hUsedPlane->SetMaximum(hUsedPlane->GetMaximum()*1.07);
 817   TLatex* tsin=new TLatex(0.15,0.84,Form("1+2*(%.4f)*sin(2*#Phi_{EP})",four2s->GetParameter(1)));
 818   tsin->SetNDC();
 819   tsin->SetTextColor(4);
 820   tsin->Draw();
 821   TLatex* tcos=new TLatex(0.15,0.77,Form("1+2*(%.4f)*cos(2*#Phi_{EP})",four2->GetParameter(1)));
 822   tcos->SetNDC();
 823   tcos->SetTextColor(1);
 824   tcos->Draw();
 825
 826   // Compute the second Fourier component for since and cosine
 827   Double_t aveCos2Phi=0.;
 828   Double_t aveSin2Phi=0.;
 829   Double_t counts=0.;
 830   for(Int_t i=1; i<=hUsedPlane->GetNbinsX(); i++){
 831     Double_t co=hUsedPlane->GetBinContent(i);
 832     Double_t phi=hUsedPlane->GetBinCenter(i);
 833     aveCos2Phi+=co*TMath::Cos(2*phi);
 834     aveSin2Phi+=co*TMath::Sin(2*phi);
 835     counts+=co;
 836   }
 837   if(counts>0){
 838     aveCos2Phi/=counts;
 839     aveSin2Phi/=counts;
 840   }
 841   printf("\n------ Fourier 2nd components of EP distribution ------\n");
 842   printf("<cos(2*Psi_EP)>=%f\n",aveCos2Phi);
 843   printf("<sin(2*Psi_EP)>=%f\n",aveSin2Phi);
 844   printf("\n");
 845
 846
 847   TCanvas* cvtotevplreso=new TCanvas(Form("cvtotevplreso%s",suffixcentr.Data()),Form("Ev plane Resolution %s",suffixcentr.Data()));
 848   cvtotevplreso->cd();
 849   hevplresos[0]->SetTitle("TPC cos2(#Psi_{A}-#Psi_{B})");
 850   hevplresos[0]->Draw();
 851   gStyle->SetOptStat(0);
 852   if(nSubRes>1){
 853     hevplresos[1]->SetLineColor(2);
 854     hevplresos[2]->SetLineColor(3);
 855     hevplresos[0]->SetTitle("cos2(#Psi_{TPC}-#Psi_{V0A})");
 856     hevplresos[1]->SetTitle("cos2(#Psi_{TPC}-#Psi_{V0C})");
 857     hevplresos[2]->SetTitle("cos2(#Psi_{V0A}-#Psi_{V0C})");
 858     hevplresos[1]->Draw("SAME");
 859     hevplresos[2]->Draw("SAME");
 860   }
 861   TLegend *leg = cvtotevplreso->BuildLegend();
 862   leg->SetLineColor(0);
 863   leg->SetFillColor(0);
 864
 865   Double_t error;
 866   Double_t resolFull=GetEventPlaneResolution(error,hevplresos[0],hevplresos[1],hevplresos[2]);
 867
 868   TPaveText* pvreso=new TPaveText(0.1,0.35,0.6,0.48,"NDC");
 869   pvreso->SetBorderSize(0);
 870   pvreso->SetFillStyle(0);
 871   pvreso->AddText(Form("Number of events = %.0f\n",hevplresos[0]->GetEntries()));
 872   pvreso->AddText(Form("Resolution on full event = %.4f#pm%.4f\n",resolFull,error));
 873   pvreso->Draw();
 874
 875   gResolVsCent->SetMarkerStyle(20);
 876   gResolVsCent->GetXaxis()->SetTitle("Centrality");
 877   gResolVsCent->GetYaxis()->SetTitle("EP Resolution");
 878
 879   Int_t colorave=kBlue-7;
 880   if(useNcollWeight) colorave=kOrange+1;
 881
 882   TCanvas* cresvscent=new TCanvas("cresvscent","ResolVsCent");
 883   cresvscent->cd();
 884   gResolVsCent->Draw("AP");
 885   TLine* lin=new TLine(gResolVsCent->GetXaxis()->GetXmin(),resolFull,gResolVsCent->GetXaxis()->GetXmax(),resolFull);
 886   lin->SetLineColor(colorave);
 887   lin->SetLineWidth(2);
 888   lin->Draw();
 889   TLatex* taveres=new TLatex(gResolVsCent->GetXaxis()->GetXmin()+0.5,resolFull*1.01,Form("<R_{2}>=%.4f",resolFull));
 890   taveres->SetTextColor(colorave);
 891   taveres->Draw();
 892   printf("\n\n---- Syst from centrality dependence of resolution ----\n");
 893   printf("MinRes=%f  MaxRes=%f AveRes=%f ---> Syst: -%f%% +%f%%\n",lowestRes,highestRes,resolFull,(resolFull-lowestRes)/resolFull,(highestRes-resolFull)/resolFull);
 894
 895   TFile* fout=new TFile(Form("EvPlanecentr%d-%d.root",minCent,maxCent),"recreate");
 896   fout->cd();
 897   hevpls->Write();
 898   for(Int_t ires=0;ires<nSubRes;ires++)hevplresos[ires]->Write();
 899   gResolVsCent->Write();
 900
 901
 902 }
 903