ROOT logo
#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TFile.h>
#include <TString.h>
#include <TH2F.h>
#include <TH1F.h>
#include <TF1.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TMultiGraph.h>
#include <TDirectoryFile.h>
#include <TList.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TLegendEntry.h>
#include <TPaveText.h>
#include <TStyle.h>
#include <TASImage.h>
#include <TPad.h>
#include <TROOT.h>
#include <TDatabasePDG.h>
#include <TParameter.h>
#include <TLatex.h>
#include <TGraphAsymmErrors.h>

#include <AliHFMassFitter.h>
#include "AliRDHFCutsD0toKpi.h"
#include "AliVertexingHFUtils.h"
#include "AliRDHFCutsDplustoKpipi.h"
#include "AliRDHFCutsDStartoKpipi.h"
#include "AliEventPlaneResolutionHandler.h"
#include <AliVertexingHFUtils.h>

#endif

/* $Id$ */

//methods for the analysis of AliAnalysisTaskSEHFv2 output
//Authors: Chiara Bianchin cbianchi@pd.infn.it
//         Giacomo Ortona  ortona@to.infn.it
//         Francesco Prino prino@to.infn.it

//global variables to be set
TString filename="AnalysisResults_train63.root";
TString suffix="v2Dplus3050Cut4upcutPIDTPC";
TString partname="Dplus";
Int_t minCent=30;
Int_t maxCent=50;
//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;
// pt and phi binning
const Int_t nptbinsnew=4;
Float_t ptbinsnew[nptbinsnew+1]={3.,4.,6.,8.,12.};
const Int_t nphibins=4;
Float_t phibinslim[nphibins+1]={0,TMath::Pi()/4,TMath::Pi()/2,3*TMath::Pi()/4,TMath::Pi()};
// mass fit configuration
Int_t rebin[nptbinsnew]={4,4,4,4};
Int_t typeb=0;  // Background: 0=expo, 1=linear, 2=pol2
Bool_t fixAlsoMass=kFALSE;
Double_t minMassForFit=1.6;
Double_t maxMassForFit=2.2;
Float_t massRangeForCounting=0.1; // GeV

Bool_t gnopng=kFALSE; //don't save in png format (only root and eps)
Int_t minPtBin[nptbinsnew]={-1,-1,-1,-1};
Int_t maxPtBin[nptbinsnew]={-1,-1,-1,-1};
Double_t effInOverEffOut=1.03;
Double_t massD;

//methods
TList *LoadMassHistos(TList *inputlist,Bool_t inoutanis);
TList* LoadResolutionHistos(TList *inputlist);
Int_t FindPtBin(Int_t nbins, Float_t* array,Float_t value);
void FillSignalGraph(TList *histlist,TGraphAsymmErrors **gSignal,TGraphAsymmErrors **gSignalfs, TGraphAsymmErrors **gSignalBC1, TGraphAsymmErrors **gSignalBC2, Bool_t inoutanis);
void DmesonsFlowAnalysis(Bool_t inoutanis=kTRUE);
void DrawEventPlane();
Double_t GetEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3);



//method implementation
//_________________________________________________________________
Double_t GetEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3){
  Double_t resolFull=1.;
  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;
}

//____________________________________________________________________
TList* LoadResolutionHistos(TList *inputlist){

  TList *outlist = new TList();
  outlist->SetName("eventplanehistlist");

  const Int_t nBins=20;
  Double_t ncoll[nBins]={1790.77,1578.44,1394.82,1236.17,1095.08,
			 969.86,859.571,759.959,669.648,589.588,
			 516.039,451.409,392.853,340.493,294.426,
			 252.385,215.484,183.284,155.101,130.963};

  Int_t minCentTimesTen=minCent*10;
  Int_t maxCentTimesTen=maxCent*10;
  TGraphErrors* gResolVsCent=new TGraphErrors(0);
  Int_t iPt=0;
  Int_t nSubRes=1;
  if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
     evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap)nSubRes=3;
  TString namereso[3]={"Reso","Reso2","Reso3"};
  TString suffixcentr=Form("centr%d_%d",minCent,maxCent);
  TH2F* hevpls=(TH2F*)inputlist->FindObject(Form("hEvPlanecentr%d_%d",minCentTimesTen,minCentTimesTen+25));
  hevpls->SetName(Form("hEvPlane%s",suffixcentr.Data()));
  hevpls->SetTitle(Form("Event Plane angle %s",suffixcentr.Data()));
    //    TH1F* hevplresos=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",mincentr,mincentr+5));QUI!!!
  TH1F* hevplresos[3];
  Int_t ncBin=minCentTimesTen/25;
  
  for(Int_t ires=0;ires<nSubRes;ires++){
    hevplresos[ires]=(TH1F*)inputlist->FindObject(Form("hEvPlane%scentr%d_%d",namereso[ires].Data(),minCentTimesTen,minCentTimesTen+25));
    if(hevplresos[ires]){
      hevplresos[ires]->SetName(Form("hEvPlane%s%s",namereso[ires].Data(),suffixcentr.Data()));
      hevplresos[ires]->SetTitle(Form("Event Plane Resolution %s%s",namereso[ires].Data(),suffixcentr.Data()));
      if(useNcollWeight){
	printf("Centr %d Bin %d  Ncoll %f\n",minCentTimesTen,ncBin,ncoll[ncBin]);
	hevplresos[ires]->Scale(ncoll[ncBin]);
      }
    }
  }
  Double_t error;
  Double_t lowestRes=1;
  Double_t highestRes=0;
  Double_t resolBin=GetEventPlaneResolution(error,hevplresos[0],hevplresos[1],hevplresos[2]);
  if(resolBin<lowestRes) lowestRes=resolBin;
  if(resolBin>highestRes) highestRes=resolBin;

  Double_t binHalfWid=25./20.;
  Double_t binCentr=(Double_t)minCentTimesTen/10.+binHalfWid;
  gResolVsCent->SetPoint(iPt,binCentr,resolBin);
  gResolVsCent->SetPointError(iPt,binHalfWid,error);
  ++iPt;
  
  for(Int_t icentr=minCentTimesTen+25;icentr<maxCentTimesTen;icentr=icentr+25){
    TH2F* h=(TH2F*)inputlist->FindObject(Form("hEvPlanecentr%d_%d",icentr,icentr+25));
    if(h)hevpls->Add(h);
    else cout<<"skipping ev plane "<<icentr<<"_"<<icentr+5<<endl;
    TH1F* htmpresos[3];
    for(Int_t ires=0;ires<nSubRes;ires++){
      htmpresos[ires]=(TH1F*)inputlist->FindObject(Form("hEvPlane%scentr%d_%d",namereso[ires].Data(),icentr,icentr+25));
      if(!htmpresos[ires])cout<<"skipping ev pl reso "<<icentr<<"_"<<icentr+25<<endl;
    }
    resolBin=GetEventPlaneResolution(error,htmpresos[0],htmpresos[1],htmpresos[2]);
    if(resolBin<lowestRes) lowestRes=resolBin;
    if(resolBin>highestRes) highestRes=resolBin;
    binCentr=(Double_t)icentr/10.+binHalfWid;
    gResolVsCent->SetPoint(iPt,binCentr,resolBin);
    gResolVsCent->SetPointError(iPt,binHalfWid,error);
    ++iPt;
    ncBin=icentr/25;
    for(Int_t ires=0;ires<nSubRes;ires++){
      if(htmpresos[ires]){
	if(useNcollWeight){
	  printf("Centr %d Bin %d  Ncoll %f\n",icentr,ncBin,ncoll[ncBin]);
	  htmpresos[ires]->Scale(ncoll[ncBin]);
	}
	hevplresos[ires]->Add(htmpresos[ires]);
      }
    }
  }
  outlist->Add(hevpls->Clone());
  for(Int_t ires=0;ires<nSubRes;ires++){
    if(hevplresos[ires]) outlist->Add(hevplresos[ires]->Clone());
  }
  gResolVsCent->SetName("gResolVsCent");
  gResolVsCent->SetTitle("Resolution vs. Centrality");
  outlist->Add(gResolVsCent->Clone());
  return outlist;
}

//__________________________________________________________
TList *LoadMassHistos(TList *inputlist,Bool_t inoutanis){
  // printf("Start load histos\n");
  //  const Int_t nptbins=cutobj->GetNPtBins();
  TList *outlist = new TList();
  outlist->SetName("azimuthalhistoslist");
  
  Int_t nphi=nphibins;
  if(inoutanis)nphi=2;
  Int_t minCentTimesTen=minCent*10;
  Int_t maxCentTimesTen=maxCent*10;
  
  //Create 2D histogram in final pt bins
  for(Int_t iFinalPtBin=0; iFinalPtBin<nptbinsnew; iFinalPtBin++){
    for(Int_t iphi=0;iphi<nphi;iphi++){
      TH1F *hMass=0x0;//=new TH1F();
      for(Int_t iPtBin=minPtBin[iFinalPtBin]; iPtBin<=maxPtBin[iFinalPtBin];iPtBin++){
	for(Int_t iHisC=minCentTimesTen; iHisC<=maxCentTimesTen-25; iHisC+=25){    
	  TString hisname=Form("hMdeltaphi_pt%dcentr%d_%d",iPtBin,iHisC,iHisC+25);
	  TH2F* htmp=(TH2F*)inputlist->FindObject(hisname.Data());
	  printf("---> Histogram: %s\n",htmp->GetName());
	  Int_t startX=htmp->FindBin(phibinslim[iphi]);
	  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
	  TH1F *h1tmp;
	  if(inoutanis){
	    if(iphi==0){
	      Int_t firstBin=htmp->FindBin(0);
	      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
	      h1tmp=(TH1F*)htmp->ProjectionY(Form("hMass%d_phi0",iPtBin),firstBin,lastBin);
	      printf("In-plane, Range: bins %d-%d -> phi %f - %f\n",firstBin,lastBin,htmp->GetXaxis()->GetBinLowEdge(firstBin),htmp->GetXaxis()->GetBinUpEdge(lastBin));
	      firstBin=htmp->FindBin(3.*TMath::Pi()/4.);
	      lastBin=htmp->FindBin(TMath::Pi()-0.0001); // -0.0001 to be sure that the upper limit of the bin is pi
	      h1tmp->Add((TH1F*)htmp->ProjectionY(Form("hMass%d",iPtBin),firstBin,lastBin));
	      printf("In-plane, Range: bins %d-%d -> phi %f - %f\n",firstBin,lastBin,htmp->GetXaxis()->GetBinLowEdge(firstBin),htmp->GetXaxis()->GetBinUpEdge(lastBin));
	    }else{
	      Int_t firstBin=htmp->FindBin(TMath::Pi()/4.);
	      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
	      h1tmp=(TH1F*)htmp->ProjectionY(Form("hMass%d_phi1",iPtBin),firstBin,lastBin);
	      printf("Out-of-plane, Range: bins %d-%d -> phi %f - %f\n",firstBin,lastBin,htmp->GetXaxis()->GetBinLowEdge(firstBin),htmp->GetXaxis()->GetBinUpEdge(lastBin));
	    }
	  }else{
	    h1tmp=(TH1F*)htmp->ProjectionY(Form("hMass%d_phi%d",iPtBin,iphi),startX,endX);
	  }
	  if(hMass==0)hMass=(TH1F*)h1tmp->Clone();
	  else hMass->Add((TH1F*)h1tmp->Clone());
	}
      }
      hMass->SetTitle(Form("hMass_pt%d_phi%d",iFinalPtBin,iphi));
      hMass->SetName(Form("hMass_pt%d_phi%d",iFinalPtBin,iphi));
      outlist->Add(hMass->Clone());
      //      hMass->DrawClone();
      delete hMass;
      hMass=0x0;
    }
  }
  return outlist;
}
//______________________________________________________________
Bool_t DefinePtBins(AliRDHFCuts *cutobj){
  Int_t nPtBinsCuts=cutobj->GetNPtBins();
  Float_t *ptlimsCuts=cutobj->GetPtBinLimits();
  //  for(Int_t iPt=0; iPt<nPtBinsCuts; iPt++) printf(" %d %f-%f\n",iPt,ptlimsCuts[iPt],ptlimsCuts[iPt+1]);
  for(Int_t iPtCuts=0; iPtCuts<nPtBinsCuts; iPtCuts++){
    for(Int_t iFinalPtBin=0; iFinalPtBin<nptbinsnew; iFinalPtBin++){  
      if(TMath::Abs(ptlimsCuts[iPtCuts]-ptbinsnew[iFinalPtBin])<0.0001){ 
	minPtBin[iFinalPtBin]=iPtCuts;
	if(iFinalPtBin>0) maxPtBin[iFinalPtBin-1]=iPtCuts-1;
      }
    }
    if(TMath::Abs(ptlimsCuts[iPtCuts]-ptbinsnew[nptbinsnew])<0.0001) maxPtBin[nptbinsnew-1]=iPtCuts-1;
  }
  if(TMath::Abs(ptbinsnew[nptbinsnew]-ptlimsCuts[nPtBinsCuts])<0.0001) maxPtBin[nptbinsnew-1]=nPtBinsCuts-1;
  for(Int_t iFinalPtBin=0; iFinalPtBin<nptbinsnew; iFinalPtBin++){
    printf("Pt bins to be merged: %d %d\n",minPtBin[iFinalPtBin],maxPtBin[iFinalPtBin]);
    if(minPtBin[iFinalPtBin]<0 || maxPtBin[iFinalPtBin]<0) return kFALSE;
  }

  return kTRUE;
}
//______________________________________________________________
Int_t GetPadNumber(Int_t ix,Int_t iy){
  return (iy)*nptbinsnew+ix+1;
}
//______________________________________________________________
void FillSignalGraph(TList *histlist,TGraphAsymmErrors **gSignal,TGraphAsymmErrors **gSignalfs, TGraphAsymmErrors **gSignalBC1, TGraphAsymmErrors **gSignalBC2, Bool_t inoutanis){

  Int_t nphi=nphibins;
  if(inoutanis)nphi=2;
  
  //Canvases for drawing histograms
  TCanvas *cDeltaPhi = new TCanvas("cinvmassdeltaphi","Invariant mass distributions",1200,700);
  TCanvas *cDeltaPhifs = new TCanvas("cinvmassdeltaphifs","Invariant mass distributions - fit with fixed sigma",1200,700);
  TCanvas *cPhiInteg = new TCanvas("cinvmass","Invariant mass distributions - #phi integrated",1200,350);
  cDeltaPhi->Divide(nptbinsnew,nphi);
  cDeltaPhifs->Divide(nptbinsnew,nphi);
  cPhiInteg->Divide(nptbinsnew,1);
  Int_t nMassBins;
  Double_t hmin,hmax;
  for(Int_t ipt=0;ipt<nptbinsnew;ipt++){    
    TH1F *histtofitfullsigma=(TH1F*)histlist->FindObject(Form("hMass_pt%d_phi0",ipt))->Clone();
    for(Int_t iphi=0;iphi<nphi;iphi++){
      Int_t ipad=GetPadNumber(ipt,iphi);
      Double_t signal=0,esignal=0;
      TH1F *histtofit=(TH1F*)histlist->FindObject(Form("hMass_pt%d_phi%d",ipt,iphi))->Clone();
      if(iphi>0)histtofitfullsigma->Add((TH1F*)histtofit->Clone());
      if(!histtofit){
	gSignal[ipt]->SetPoint(iphi,iphi,signal);
	gSignal[ipt]->SetPointError(iphi,0,0,esignal,esignal);
	return;
      }
      histtofit->SetTitle(Form("%.1f<p_{t}<%.1f, #phi%d",ptbinsnew[ipt],ptbinsnew[ipt+1],iphi));
      nMassBins=histtofit->GetNbinsX();
      hmin=TMath::Max(minMassForFit,histtofit->GetBinLowEdge(2));
      hmax=TMath::Min(maxMassForFit,histtofit->GetBinLowEdge(nMassBins-2));
      histtofit->Rebin(rebin[ipt]);
      AliHFMassFitter fitter(histtofit,hmin,hmax,1,typeb);
      fitter.SetInitialGaussianMean(massD);
      //      fitter.SetInitialGaussianSigma(0.012);
      Bool_t ok=fitter.MassFitter(kFALSE);
      if(ok){
	fitter.DrawHere(cDeltaPhi->cd(ipad),3,1);
	fitter.Signal(3,signal,esignal);
      }
      gSignal[ipt]->SetPoint(iphi,iphi,signal);
      gSignal[ipt]->SetPointError(iphi,0,0,esignal,esignal);
      TF1* fB1=fitter.GetBackgroundFullRangeFunc();
      TF1* fB2=fitter.GetBackgroundRecalcFunc();
      Float_t minBinSum=histtofit->FindBin(massD-massRangeForCounting);
      Float_t maxBinSum=histtofit->FindBin(massD+massRangeForCounting);
      Float_t cntSig1=0.;
      Float_t cntSig2=0.;
      Float_t cntErr=0.;
      for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
	Float_t bkg1=fB1 ? fB1->Eval(histtofit->GetBinCenter(iMB)) : 0;
	Float_t bkg2=fB2 ? fB2->Eval(histtofit->GetBinCenter(iMB)) : 0;
	//	printf("PtBin %d MassBin%d BC %f %f %f\n",ipt,iMB,histtofit->GetBinContent(iMB),bkg1,bkg2);
	cntSig1+=(histtofit->GetBinContent(iMB)-bkg1);
	cntSig2+=(histtofit->GetBinContent(iMB)-bkg2);
	cntErr+=(histtofit->GetBinContent(iMB));
      }
      cntErr=TMath::Sqrt(cntErr);
      gSignalBC1[ipt]->SetPoint(iphi,iphi,cntSig1);
      gSignalBC1[ipt]->SetPointError(iphi,0,0,cntErr,cntErr);
      gSignalBC2[ipt]->SetPoint(iphi,iphi,cntSig2);
      gSignalBC2[ipt]->SetPointError(iphi,0,0,cntErr,cntErr);
    }
    //fit for fixed sigma
    histtofitfullsigma->SetTitle(Form("%.1f<p_{t}<%.1f",ptbinsnew[ipt],ptbinsnew[ipt+1]));
    nMassBins=histtofitfullsigma->GetNbinsX();
    hmin=TMath::Max(minMassForFit,histtofitfullsigma->GetBinLowEdge(2));
    hmax=TMath::Min(maxMassForFit,histtofitfullsigma->GetBinLowEdge(nMassBins-2));
    histtofitfullsigma->Rebin(rebin[ipt]);
    AliHFMassFitter fitter(histtofitfullsigma,hmin,hmax,1,typeb);
    fitter.SetInitialGaussianMean(massD);
    //      fitter.SetInitialGaussianSigma(0.012);
    //fitter.SetFixGaussianSigma(0.012);
    Bool_t ok=fitter.MassFitter(kFALSE);
    if(ok){
      fitter.DrawHere(cPhiInteg->cd(ipt+1),3,1);
    }
    Double_t sigma=fitter.GetSigma();
    Double_t massFromFit=fitter.GetMean();
    for(Int_t iphi=0;iphi<nphi;iphi++){
      Int_t ipad=GetPadNumber(ipt,iphi);
      TH1F *histtofit=(TH1F*)histlist->FindObject(Form("hMass_pt%d_phi%d",ipt,iphi))->Clone();
      histtofit->SetTitle(Form("%.1f<p_{t}<%.1f, #phi%d",ptbinsnew[ipt],ptbinsnew[ipt+1],iphi));
      nMassBins=histtofit->GetNbinsX();
      hmin=TMath::Max(minMassForFit,histtofit->GetBinLowEdge(2));
      hmax=TMath::Min(maxMassForFit,histtofit->GetBinLowEdge(nMassBins-2));
      histtofit->Rebin(rebin[ipt]);
      AliHFMassFitter fitter2(histtofit,hmin,hmax,1,typeb);
      fitter2.SetInitialGaussianMean(massD);
      fitter2.SetFixGaussianSigma(sigma);
      if(fixAlsoMass) fitter2.SetFixGaussianMean(massFromFit);
      Bool_t ok2=fitter2.MassFitter(kFALSE);
      Double_t signal=0,esignal=0;
      if(ok2){
	fitter2.DrawHere(cDeltaPhifs->cd(ipad),3,1);
	fitter2.Signal(3,signal,esignal);
      }
      gSignalfs[ipt]->SetPoint(iphi,iphi,signal);
      gSignalfs[ipt]->SetPointError(iphi,0,0,esignal,esignal);
    }
  }//end loop on pt bin


  cDeltaPhi->SaveAs(Form("InvMassDeltaPhi_%s.eps",suffix.Data()));
  cDeltaPhifs->SaveAs(Form("InvMassDeltaPhi_fs_%s.eps",suffix.Data()));
  cDeltaPhifs->SaveAs(Form("InvMassDeltaPhi_fs_%s.root",suffix.Data()));
  cPhiInteg->SaveAs(Form("InvMassfullphi_%s.eps",suffix.Data()));

  if(!gnopng){
    cDeltaPhi->SaveAs(Form("InvMassDeltaPhi_%s.png",suffix.Data()));
    cDeltaPhifs->SaveAs(Form("InvMassDeltaPhi_fs_%s.png",suffix.Data()));
    cPhiInteg->SaveAs(Form("InvMassfullphi_%s.png",suffix.Data()));
  }
  //cDeltaPhifs->DrawClone();
  //cDeltaPhifs->Close();
 
}
//______________________________________________________________
void DmesonsFlowAnalysis(Bool_t inoutanis){
  TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
  TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());

  AliRDHFCuts *cutsobj=0x0;
  //Load input data from AliAnalysisTaskSEHFv2
  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;
  }
  if(partname.Contains("Dzero")) {
    cutsobj=((AliRDHFCutsD0toKpi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
    massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
  }
  if(partname.Contains("Dplus")){
    cutsobj=((AliRDHFCutsDplustoKpipi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
    massD=TDatabasePDG::Instance()->GetParticle(411)->Mass();
  }
  if(partname.Contains("Dstar")) {
    cutsobj=((AliRDHFCutsDStartoKpipi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
    massD=(TDatabasePDG::Instance()->GetParticle(413)->Mass() - TDatabasePDG::Instance()->GetParticle(421)->Mass()); 
  }

  TList *list =(TList*)dir->Get(listname.Data());
  if(!list){
    printf("list %s not found in file, please check list name\n",listname.Data());return;
  }
  if(!cutsobj){
    printf("cut object not found in file, please check keylist number\n");return;
  }
  //Define new pt bins
  if(!DefinePtBins(cutsobj)){
    printf("cut not define pt bins\n");return;
  }
  
  //Load mass histograms corresponding to the required centrality, pt range and phi binning
  printf("Load mass histos \n");
  TList *histlist=LoadMassHistos(list,inoutanis);
  TString aniss="";
  if(inoutanis)aniss+="anis";
  histlist->SaveAs(Form("v2Histograms_%d_%d_%s_%s.root",minCent,maxCent,aniss.Data(),suffix.Data()),"RECREATE");

  Int_t nphi=nphibins;
  if(inoutanis)nphi=2;


  printf("average pt for pt bin \n");
  //average pt for pt bin
  AliVertexingHFUtils *utils=new AliVertexingHFUtils();
  Int_t minCentTimesTen=minCent*10;
  Int_t maxCentTimesTen=maxCent*10;
  TH2F* hmasspt=(TH2F*)list->FindObject(Form("hMPtCandcentr%d_%d",minCentTimesTen,minCentTimesTen+25));
  for(Int_t icent=minCentTimesTen+25;icent<maxCentTimesTen;icent=icent+25)hmasspt->Add((TH2F*)list->FindObject(Form("hMPtCandcentr%d_%d",icent,icent+25)));
  Float_t averagePt[nptbinsnew];
  Float_t errorPt[nptbinsnew];
  for(Int_t ipt=0;ipt<nptbinsnew;ipt++){
    Int_t binMin=hmasspt->FindBin(ptbinsnew[ipt]);
    Int_t binMax=hmasspt->FindBin(ptbinsnew[ipt+1]-0.001);
    if(TMath::Abs(hmasspt->GetXaxis()->GetBinLowEdge(binMin)-ptbinsnew[ipt])>0.001 || 
       TMath::Abs(hmasspt->GetXaxis()->GetBinUpEdge(binMax)-ptbinsnew[ipt+1])>0.001){
      printf("Error in pt bin limits for projection!\n");
      return;
    }
    TH1F *histtofit = (TH1F*)hmasspt->ProjectionY("_py",binMin,binMax);
    Int_t nMassBins=histtofit->GetNbinsX();
    Double_t hmin=histtofit->GetBinLowEdge(2); // need wide range for <pt>
    Double_t hmax=histtofit->GetBinLowEdge(nMassBins-2); // need wide range for <pt>
    AliHFMassFitter fitter(histtofit,hmin,hmax,1);
    fitter.MassFitter(kFALSE);
    Float_t massFromFit=fitter.GetMean();
    Float_t sigmaFromFit=fitter.GetSigma();
    TF1* funcB2=fitter.GetBackgroundRecalcFunc();
    utils->AveragePt(averagePt[ipt],errorPt[ipt],ptbinsnew[ipt],ptbinsnew[ipt+1],hmasspt,massFromFit,sigmaFromFit,funcB2,2.5,4.5,0.,3.,1);
  }
  printf("Average pt\n");
  for(Int_t ipt=0;ipt<nptbinsnew;ipt++) printf("%f +- %f\n",averagePt[ipt],errorPt[ipt]);

  printf("Fill TGraphs for signal \n");
  //Fill TGraphs for signal
  TGraphAsymmErrors *gSignal[nptbinsnew];
  TGraphAsymmErrors *gSignalfs[nptbinsnew];
  TGraphAsymmErrors *gSignalBC1[nptbinsnew];
  TGraphAsymmErrors *gSignalBC2[nptbinsnew];
  for(Int_t i=0;i<nptbinsnew;i++){
    gSignal[i]=new TGraphAsymmErrors(nphi);
    gSignal[i]->SetName(Form("gasigpt%d",i));
    gSignal[i]->SetTitle(Form("Signal %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
    gSignal[i]->SetMarkerStyle(25);
    gSignalfs[i]=new TGraphAsymmErrors(nphi);
    gSignalfs[i]->SetName(Form("gasigfspt%d",i));
    gSignalfs[i]->SetTitle(Form("Signal (fixed sigma) %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
    gSignalfs[i]->SetMarkerStyle(21);
    gSignalBC1[i]=new TGraphAsymmErrors(nphi);
    gSignalBC1[i]->SetName(Form("gasigBC1pt%d",i));
    gSignalBC1[i]->SetTitle(Form("Signal (BC1) %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
    gSignalBC2[i]=new TGraphAsymmErrors(nphi);
    gSignalBC2[i]->SetName(Form("gasigBC2pt%d",i));
    gSignalBC2[i]->SetTitle(Form("Signal (BC2) %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
  }
  FillSignalGraph(histlist,gSignal,gSignalfs,gSignalBC1,gSignalBC2,inoutanis);

  //EP resolution
  // TList* resolhist=LoadResolutionHistos(list);
  // TString suffixcentr=Form("centr%d_%d",minCent,maxCent);
  // TH1F* hevplresos[3];
  // TString namereso[3]={"Reso","Reso2","Reso3"};
  // Int_t nSubRes=1;
  // if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
  //    evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap)nSubRes=3;
  // for(Int_t ires=0;ires<nSubRes;ires++){
  //   hevplresos[ires]=(TH1F*)resolhist->FindObject(Form("hEvPlane%s%s",namereso[ires].Data(),suffixcentr.Data()));
  // }
  // Double_t errorres;
  // Double_t resol=GetEventPlaneResolution(errorres,hevplresos[0],hevplresos[1],hevplresos[2]);
  AliEventPlaneResolutionHandler* epres=new AliEventPlaneResolutionHandler();
  epres->SetEventPlane(evPlane);
  epres->SetResolutionOption(evPlaneRes);
  epres->SetUseNcollWeights(useNcollWeight);
  Double_t resol=epres->GetEventPlaneResolution(minCent,maxCent);
  delete epres;
  printf("Event plane resolution %f\n",resol);

  printf("Compute v2 \n");
  //compute v2
  TCanvas *cv2fs =new TCanvas("v2_fs","v2 Fit with fixed sigma");
  TCanvas *cv2 =new TCanvas("v2","v2 - systematic on yield extraction");
  TGraphAsymmErrors *gv2=new TGraphAsymmErrors(nptbinsnew);
  TGraphAsymmErrors *gv2fs=new TGraphAsymmErrors(nptbinsnew);
  TGraphAsymmErrors *gv2BC1=new TGraphAsymmErrors(nptbinsnew);
  TGraphAsymmErrors *gv2BC2=new TGraphAsymmErrors(nptbinsnew);
  TGraphAsymmErrors *grelSystEff=new TGraphAsymmErrors(nptbinsnew);

  gv2->SetName("gav2");
  gv2->SetTitle("v_{2}(p_{t})");
  gv2fs->SetName("gav2fs");
  gv2fs->SetTitle("v_{2}(p_{t}) (fit with fixed sigma)");
  gv2BC1->SetName("gav2BC1");
  gv2BC1->SetTitle("v_{2}(p_{t}) (bin counting 1)");
  gv2BC2->SetName("gav2BC2");
  gv2BC2->SetTitle("v_{2}(p_{t}) (bin counting 2)");
  grelSystEff->SetName("grelSystEff");
  grelSystEff->SetTitle("SystErrEff");

  //Prepare output file
  TFile *fout=new TFile(Form("v2Output_%d_%d_%s_%s.root",minCent,maxCent,aniss.Data(),suffix.Data()),"RECREATE");

  for(Int_t ipt=0;ipt<nptbinsnew;ipt++){
    fout->cd();
    gSignal[ipt]->Write();
    gSignalfs[ipt]->Write();
    gSignalBC1[ipt]->Write();
    gSignalBC2[ipt]->Write();

    if(inoutanis){
      //v2 from in-out anisotropy
      Double_t *y,*yfs,*ybc1,*ybc2;
      y=gSignal[ipt]->GetY();
      yfs=gSignalfs[ipt]->GetY();
      ybc1=gSignalBC1[ipt]->GetY();
      ybc2=gSignalBC2[ipt]->GetY();

      Double_t nIn=y[0];
      Double_t nOut=y[1];
      Double_t enIn=gSignal[ipt]->GetErrorY(0);
      Double_t enOut=gSignal[ipt]->GetErrorY(1);
      Double_t anis=(nIn-nOut)/(nIn+nOut);
      //      Double_t eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
      Double_t eAnis=2./((nIn+nOut)*(nIn+nOut))*TMath::Sqrt(nIn*nIn*enOut*enOut+nOut*nOut*enIn*enIn);
      Double_t v2=anis*TMath::Pi()/4./resol;
      Double_t ev2=eAnis*TMath::Pi()/4./resol;
      gv2->SetPoint(ipt,averagePt[ipt],v2);
      gv2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);

      nIn=yfs[0];
      nOut=yfs[1];
      enIn=gSignalfs[ipt]->GetErrorY(0);
      enOut=gSignalfs[ipt]->GetErrorY(1);
      anis=(nIn-nOut)/(nIn+nOut);
      //     eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
      eAnis=2./((nIn+nOut)*(nIn+nOut))*TMath::Sqrt(nIn*nIn*enOut*enOut+nOut*nOut*enIn*enIn);
      v2=anis*TMath::Pi()/4./resol;
      ev2=eAnis*TMath::Pi()/4./resol;
      gv2fs->SetPoint(ipt,averagePt[ipt],v2);
      gv2fs->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
      Double_t anis1=(nIn-nOut*effInOverEffOut)/(nIn+nOut*effInOverEffOut);
      Double_t anis2=(nIn*effInOverEffOut-nOut)/(nIn*effInOverEffOut+nOut);
      Double_t systEffUp=0.,systEffDown=0.;
      if(anis1>anis && anis1>anis2) systEffUp=anis1/anis;
      if(anis2>anis && anis2>anis1) systEffUp=anis2/anis;
      if(anis1<anis && anis1<anis2) systEffDown=anis1/anis;
      if(anis2<anis && anis2<anis1) systEffDown=anis2/anis;
      printf(" Bin %d <pt>=%.3f  v2=%f+-%f systEff=%f %f\n",ipt,averagePt[ipt],v2,ev2,systEffUp*v2,systEffDown*v2);
      grelSystEff->SetPoint(ipt,averagePt[ipt],v2);
      grelSystEff->SetPointError(ipt,0.4,0.4,v2*(1-systEffDown),v2*(systEffUp-1));

      nIn=ybc1[0];
      nOut=ybc1[1];
      enIn=gSignalBC1[ipt]->GetErrorY(0);
      enOut=gSignalBC1[ipt]->GetErrorY(1);
      anis=(nIn-nOut)/(nIn+nOut);
      //     eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
      eAnis=2./((nIn+nOut)*(nIn+nOut))*TMath::Sqrt(nIn*nIn*enOut*enOut+nOut*nOut*enIn*enIn);
      v2=anis*TMath::Pi()/4./resol;
      ev2=eAnis*TMath::Pi()/4./resol;
      gv2BC1->SetPoint(ipt,averagePt[ipt],v2);
      gv2BC1->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);

      nIn=ybc2[0];
      nOut=ybc2[1];
      enIn=gSignalBC2[ipt]->GetErrorY(0);
      enOut=gSignalBC2[ipt]->GetErrorY(1);
      anis=(nIn-nOut)/(nIn+nOut);
      //     eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
      eAnis=2./((nIn+nOut)*(nIn+nOut))*TMath::Sqrt(nIn*nIn*enOut*enOut+nOut*nOut*enIn*enIn);
      v2=anis*TMath::Pi()/4./resol;
      ev2=eAnis*TMath::Pi()/4./resol;
      gv2BC2->SetPoint(ipt,averagePt[ipt],v2);
      gv2BC2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);

   }else{
      TF1 *flowFunc = new TF1("flow","[0]*(1.+2.*[1]*TMath::Cos(2.*x))");
      //v2 from fit to Deltaphi distribution
      gSignal[ipt]->Fit(flowFunc);
      Double_t v2 = flowFunc->GetParameter(1)/resol;
      Double_t ev2=flowFunc->GetParError(1)/resol;
      gv2->SetPoint(ipt,averagePt[ipt],v2);
      gv2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);

      gSignalfs[ipt]->Fit(flowFunc);
      v2 = flowFunc->GetParameter(1)/resol;
      ev2=flowFunc->GetParError(1)/resol;
      gv2fs->SetPoint(ipt,averagePt[ipt],v2);
      gv2fs->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);

      gSignalBC1[ipt]->Fit(flowFunc);
      v2 = flowFunc->GetParameter(1)/resol;
      ev2=flowFunc->GetParError(1)/resol;
      gv2BC1->SetPoint(ipt,averagePt[ipt],v2);
      gv2BC1->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);

      gSignalBC2[ipt]->Fit(flowFunc);
      v2 = flowFunc->GetParameter(1)/resol;
      ev2=flowFunc->GetParError(1)/resol;
      gv2BC2->SetPoint(ipt,averagePt[ipt],v2);
      gv2BC2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
    }
  }

  gv2->SetMarkerStyle(22);
  gv2->SetMarkerColor(2);
  gv2->SetLineColor(2);
  gv2fs->SetMarkerStyle(20);
  gv2fs->SetMarkerColor(1);
  gv2fs->SetLineColor(1);
  gv2BC1->SetMarkerStyle(23);
  gv2BC1->SetMarkerColor(kGreen+1);
  gv2BC1->SetLineColor(kGreen+1);
  gv2BC2->SetMarkerStyle(26);
  gv2BC2->SetMarkerColor(kGreen+1);
  gv2BC2->SetLineColor(kGreen+1);
  
  cv2->cd();
  gv2fs->Draw("AP");
  gv2fs->SetMinimum(-0.2);
  gv2fs->SetMaximum(0.6);
  gv2fs->GetYaxis()->SetTitle("v_{2}");
  gv2fs->GetXaxis()->SetTitle("p_{t} [GeV/c]");
  gv2->Draw("PSAME");
  gv2BC1->Draw("PSAME");
  gv2BC2->Draw("PSAME");
  TLegend* leg=new TLegend(0.65,0.7,0.89,0.89);
  leg->SetFillStyle(0);
  leg->SetBorderSize(0);
  leg->SetTextFont(42);
  leg->AddEntry(gv2fs,"Fixed sigma","P")->SetTextColor(gv2fs->GetMarkerColor());
  leg->AddEntry(gv2,"Free sigma","P")->SetTextColor(gv2->GetMarkerColor());
  leg->AddEntry(gv2BC1,"Bin counting 1","P")->SetTextColor(gv2BC1->GetMarkerColor());
  leg->AddEntry(gv2BC2,"Bin counting 2","P")->SetTextColor(gv2BC2->GetMarkerColor());
  leg->Draw();
  cv2fs->cd();
  gv2fs->Draw("AP");
  gv2fs->GetYaxis()->SetTitle("v_{2}");
  gv2fs->GetXaxis()->SetTitle("p_{t} [GeV/c]");

  fout->cd();
  gv2->Write();
  gv2fs->Write();
  gv2BC1->Write();
  gv2BC2->Write();
  grelSystEff->Write();
}
//___________________________________________________________
Int_t FindPtBin(Int_t nbins, Float_t* array,Float_t value){
  for (Int_t i=0;i<nbins;i++){
    if(value>=array[i] && value<array[i+1]){
      return i;
    }
  }
  cout<<value<< " out of range "<<array[0]<<", "<<array[nbins]<<endl;
  return -1;
}

//__________________________________________________________
void DrawEventPlane(){
  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;
  }
  TList *list =(TList*)dir->Get(listname.Data());
  if(!list){
    printf("list %s not found in file, please check list name\n",listname.Data());return;
  } 
  if(evPlane==AliEventPlaneResolutionHandler::kTPCPosEta &&
     evPlaneRes==AliEventPlaneResolutionHandler::kTwoEtaSub){
    printf("Wrong setting of event plane resolution forced it to kTwoSub\n");
    evPlaneRes=AliEventPlaneResolutionHandler::kTwoRandSub;
  }
  TList* resolhist=LoadResolutionHistos(list);
  TGraphErrors* gResolVsCent=(TGraphErrors*)resolhist->FindObject("gResolVsCent");
  Double_t lowestRes=1;
  Double_t highestRes=0;
  for(Int_t ipt=0; ipt<gResolVsCent->GetN(); ipt++){    
    Double_t x,resolBin;
    gResolVsCent->GetPoint(ipt,x,resolBin);
    if(resolBin<lowestRes) lowestRes=resolBin;
    if(resolBin>highestRes) highestRes=resolBin;
  }
  TString suffixcentr=Form("centr%d_%d",minCent,maxCent);
  TH2F* hevpls=(TH2F*)resolhist->FindObject(Form("hEvPlane%s",suffixcentr.Data()));
  TH1F* hevplresos[3];
  TString namereso[3]={"Reso","Reso2","Reso3"};
  Int_t nSubRes=1;
  if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
     evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap)nSubRes=3;
  for(Int_t ires=0;ires<nSubRes;ires++){
    hevplresos[ires]=(TH1F*)resolhist->FindObject(Form("hEvPlane%s%s",namereso[ires].Data(),suffixcentr.Data()));
  }

  gStyle->SetCanvasColor(0);
  gStyle->SetTitleFillColor(0);
  gStyle->SetStatColor(0);
  gStyle->SetOptStat(0);
  //gStyle->SetPalette(1);
  gStyle->SetOptFit(1);
  gStyle->SetLabelFont(42,"XYZ");
  gStyle->SetTextFont(42);
  gStyle->SetStatFont(42);
  //  gStyle->SetTitleFont(42);
  gStyle->SetStatY(0.89);
  gStyle->SetStatX(0.89);
  //  gStyle->SetTitleAlign(22);
  gStyle->SetTitleFont(42,"xyzg");

  TH1F* htpc = (TH1F*)hevpls->ProjectionX();
  TH1F* hv0 = (TH1F*)hevpls->ProjectionY();
  TH1F* hUsedPlane;
  if(evPlane==AliEventPlaneResolutionHandler::kVZERO ||
     evPlane==AliEventPlaneResolutionHandler::kVZEROA ||
     evPlane==AliEventPlaneResolutionHandler::kVZEROC) hUsedPlane=(TH1F*)hv0->Clone("hUsedPlane");
  else hUsedPlane=(TH1F*)htpc->Clone("hUsedPlane");

  TCanvas* cvtotevpl=new TCanvas(Form("cvtotevpl%s",suffixcentr.Data()),Form("Ev plane %s",suffixcentr.Data()),1200,400);
  cvtotevpl->Divide(3,1);
  cvtotevpl->cd(1);
  hevpls->SetTitleFont(42);
  hevpls->Draw("COLZ");
  cvtotevpl->cd(2);
  htpc->Draw();
  htpc->Fit("pol0");
  cvtotevpl->cd(3);
  hv0->Draw();
  hv0->Fit("pol0");
  
  TCanvas* cvfitevpl=new TCanvas(Form("cvditevpl%s",suffixcentr.Data()),Form("Fit to Ev plane %s",suffixcentr.Data()));
  cvfitevpl->cd();
  hUsedPlane->Draw();
  TF1* four2=new TF1("four2","[0]*(1+2*[1]*cos(2*x))",hUsedPlane->GetXaxis()->GetXmin(),hUsedPlane->GetXaxis()->GetXmax());
  TF1* four2s=new TF1("four2s","[0]*(1+2*[1]*sin(2*x))",hUsedPlane->GetXaxis()->GetXmin(),hUsedPlane->GetXaxis()->GetXmax());
  hUsedPlane->Fit(four2s);
  hUsedPlane->Fit(four2);
  four2->SetLineColor(1);
  four2s->SetLineColor(4);
  four2->Draw("same");
  four2s->Draw("same");
  hUsedPlane->SetMaximum(hUsedPlane->GetMaximum()*1.07);
  TLatex* tsin=new TLatex(0.15,0.84,Form("1+2*(%.4f)*sin(2*#Phi_{EP})",four2s->GetParameter(1)));
  tsin->SetNDC();
  tsin->SetTextColor(4);
  tsin->Draw();
  TLatex* tcos=new TLatex(0.15,0.77,Form("1+2*(%.4f)*cos(2*#Phi_{EP})",four2->GetParameter(1)));
  tcos->SetNDC();
  tcos->SetTextColor(1);
  tcos->Draw();
  
  // Compute the second Fourier component for since and cosine
  Double_t aveCos2Phi=0.;
  Double_t aveSin2Phi=0.;
  Double_t counts=0.;
  for(Int_t i=1; i<=hUsedPlane->GetNbinsX(); i++){
    Double_t co=hUsedPlane->GetBinContent(i);
    Double_t phi=hUsedPlane->GetBinCenter(i);
    aveCos2Phi+=co*TMath::Cos(2*phi);
    aveSin2Phi+=co*TMath::Sin(2*phi);
    counts+=co;
  }
  if(counts>0){ 
    aveCos2Phi/=counts;
    aveSin2Phi/=counts;
  }
  printf("\n------ Fourier 2nd components of EP distribution ------\n");
  printf("<cos(2*Psi_EP)>=%f\n",aveCos2Phi);
  printf("<sin(2*Psi_EP)>=%f\n",aveSin2Phi);
  printf("\n");


  TCanvas* cvtotevplreso=new TCanvas(Form("cvtotevplreso%s",suffixcentr.Data()),Form("Ev plane Resolution %s",suffixcentr.Data()));
  cvtotevplreso->cd();
  hevplresos[0]->SetTitle("TPC cos2(#Psi_{A}-#Psi_{B})");
  hevplresos[0]->Draw();
  gStyle->SetOptStat(0);
  if(nSubRes>1){
    hevplresos[1]->SetLineColor(2);
    hevplresos[2]->SetLineColor(3);
    hevplresos[0]->SetTitle("cos2(#Psi_{TPC}-#Psi_{V0A})");
    hevplresos[1]->SetTitle("cos2(#Psi_{TPC}-#Psi_{V0C})");
    hevplresos[2]->SetTitle("cos2(#Psi_{V0A}-#Psi_{V0C})");
    hevplresos[1]->Draw("SAME");
    hevplresos[2]->Draw("SAME");
  }
  TLegend *leg = cvtotevplreso->BuildLegend();
  leg->SetLineColor(0);
  leg->SetFillColor(0);

  Double_t error;
  Double_t resolFull=GetEventPlaneResolution(error,hevplresos[0],hevplresos[1],hevplresos[2]);
  
  TPaveText* pvreso=new TPaveText(0.1,0.35,0.6,0.48,"NDC");
  pvreso->SetBorderSize(0);
  pvreso->SetFillStyle(0);
  pvreso->AddText(Form("Number of events = %.0f\n",hevplresos[0]->GetEntries()));
  pvreso->AddText(Form("Resolution on full event = %.4f#pm%.4f\n",resolFull,error));
  pvreso->Draw();

  gResolVsCent->SetMarkerStyle(20);
  gResolVsCent->GetXaxis()->SetTitle("Centrality");
  gResolVsCent->GetYaxis()->SetTitle("EP Resolution");
 
  Int_t colorave=kBlue-7;
  if(useNcollWeight) colorave=kOrange+1;

  TCanvas* cresvscent=new TCanvas("cresvscent","ResolVsCent");
  cresvscent->cd();
  gResolVsCent->Draw("AP");
  TLine* lin=new TLine(gResolVsCent->GetXaxis()->GetXmin(),resolFull,gResolVsCent->GetXaxis()->GetXmax(),resolFull);
  lin->SetLineColor(colorave);
  lin->SetLineWidth(2);
  lin->Draw();
  TLatex* taveres=new TLatex(gResolVsCent->GetXaxis()->GetXmin()+0.5,resolFull*1.01,Form("<R_{2}>=%.4f",resolFull));
  taveres->SetTextColor(colorave);
  taveres->Draw();
  printf("\n\n---- Syst from centrality dependence of resolution ----\n");
  printf("MinRes=%f  MaxRes=%f AveRes=%f ---> Syst: -%f%% +%f%%\n",lowestRes,highestRes,resolFull,(resolFull-lowestRes)/resolFull,(highestRes-resolFull)/resolFull);

  TFile* fout=new TFile(Form("EvPlanecentr%d-%d.root",minCent,maxCent),"recreate");
  fout->cd();
  hevpls->Write();
  for(Int_t ires=0;ires<nSubRes;ires++)hevplresos[ires]->Write();
  gResolVsCent->Write();


}

 DmesonsFlowAnalysis.C:1
 DmesonsFlowAnalysis.C:2
 DmesonsFlowAnalysis.C:3
 DmesonsFlowAnalysis.C:4
 DmesonsFlowAnalysis.C:5
 DmesonsFlowAnalysis.C:6
 DmesonsFlowAnalysis.C:7
 DmesonsFlowAnalysis.C:8
 DmesonsFlowAnalysis.C:9
 DmesonsFlowAnalysis.C:10
 DmesonsFlowAnalysis.C:11
 DmesonsFlowAnalysis.C:12
 DmesonsFlowAnalysis.C:13
 DmesonsFlowAnalysis.C:14
 DmesonsFlowAnalysis.C:15
 DmesonsFlowAnalysis.C:16
 DmesonsFlowAnalysis.C:17
 DmesonsFlowAnalysis.C:18
 DmesonsFlowAnalysis.C:19
 DmesonsFlowAnalysis.C:20
 DmesonsFlowAnalysis.C:21
 DmesonsFlowAnalysis.C:22
 DmesonsFlowAnalysis.C:23
 DmesonsFlowAnalysis.C:24
 DmesonsFlowAnalysis.C:25
 DmesonsFlowAnalysis.C:26
 DmesonsFlowAnalysis.C:27
 DmesonsFlowAnalysis.C:28
 DmesonsFlowAnalysis.C:29
 DmesonsFlowAnalysis.C:30
 DmesonsFlowAnalysis.C:31
 DmesonsFlowAnalysis.C:32
 DmesonsFlowAnalysis.C:33
 DmesonsFlowAnalysis.C:34
 DmesonsFlowAnalysis.C:35
 DmesonsFlowAnalysis.C:36
 DmesonsFlowAnalysis.C:37
 DmesonsFlowAnalysis.C:38
 DmesonsFlowAnalysis.C:39
 DmesonsFlowAnalysis.C:40
 DmesonsFlowAnalysis.C:41
 DmesonsFlowAnalysis.C:42
 DmesonsFlowAnalysis.C:43
 DmesonsFlowAnalysis.C:44
 DmesonsFlowAnalysis.C:45
 DmesonsFlowAnalysis.C:46
 DmesonsFlowAnalysis.C:47
 DmesonsFlowAnalysis.C:48
 DmesonsFlowAnalysis.C:49
 DmesonsFlowAnalysis.C:50
 DmesonsFlowAnalysis.C:51
 DmesonsFlowAnalysis.C:52
 DmesonsFlowAnalysis.C:53
 DmesonsFlowAnalysis.C:54
 DmesonsFlowAnalysis.C:55
 DmesonsFlowAnalysis.C:56
 DmesonsFlowAnalysis.C:57
 DmesonsFlowAnalysis.C:58
 DmesonsFlowAnalysis.C:59
 DmesonsFlowAnalysis.C:60
 DmesonsFlowAnalysis.C:61
 DmesonsFlowAnalysis.C:62
 DmesonsFlowAnalysis.C:63
 DmesonsFlowAnalysis.C:64
 DmesonsFlowAnalysis.C:65
 DmesonsFlowAnalysis.C:66
 DmesonsFlowAnalysis.C:67
 DmesonsFlowAnalysis.C:68
 DmesonsFlowAnalysis.C:69
 DmesonsFlowAnalysis.C:70
 DmesonsFlowAnalysis.C:71
 DmesonsFlowAnalysis.C:72
 DmesonsFlowAnalysis.C:73
 DmesonsFlowAnalysis.C:74
 DmesonsFlowAnalysis.C:75
 DmesonsFlowAnalysis.C:76
 DmesonsFlowAnalysis.C:77
 DmesonsFlowAnalysis.C:78
 DmesonsFlowAnalysis.C:79
 DmesonsFlowAnalysis.C:80
 DmesonsFlowAnalysis.C:81
 DmesonsFlowAnalysis.C:82
 DmesonsFlowAnalysis.C:83
 DmesonsFlowAnalysis.C:84
 DmesonsFlowAnalysis.C:85
 DmesonsFlowAnalysis.C:86
 DmesonsFlowAnalysis.C:87
 DmesonsFlowAnalysis.C:88
 DmesonsFlowAnalysis.C:89
 DmesonsFlowAnalysis.C:90
 DmesonsFlowAnalysis.C:91
 DmesonsFlowAnalysis.C:92
 DmesonsFlowAnalysis.C:93
 DmesonsFlowAnalysis.C:94
 DmesonsFlowAnalysis.C:95
 DmesonsFlowAnalysis.C:96
 DmesonsFlowAnalysis.C:97
 DmesonsFlowAnalysis.C:98
 DmesonsFlowAnalysis.C:99
 DmesonsFlowAnalysis.C:100
 DmesonsFlowAnalysis.C:101
 DmesonsFlowAnalysis.C:102
 DmesonsFlowAnalysis.C:103
 DmesonsFlowAnalysis.C:104
 DmesonsFlowAnalysis.C:105
 DmesonsFlowAnalysis.C:106
 DmesonsFlowAnalysis.C:107
 DmesonsFlowAnalysis.C:108
 DmesonsFlowAnalysis.C:109
 DmesonsFlowAnalysis.C:110
 DmesonsFlowAnalysis.C:111
 DmesonsFlowAnalysis.C:112
 DmesonsFlowAnalysis.C:113
 DmesonsFlowAnalysis.C:114
 DmesonsFlowAnalysis.C:115
 DmesonsFlowAnalysis.C:116
 DmesonsFlowAnalysis.C:117
 DmesonsFlowAnalysis.C:118
 DmesonsFlowAnalysis.C:119
 DmesonsFlowAnalysis.C:120
 DmesonsFlowAnalysis.C:121
 DmesonsFlowAnalysis.C:122
 DmesonsFlowAnalysis.C:123
 DmesonsFlowAnalysis.C:124
 DmesonsFlowAnalysis.C:125
 DmesonsFlowAnalysis.C:126
 DmesonsFlowAnalysis.C:127
 DmesonsFlowAnalysis.C:128
 DmesonsFlowAnalysis.C:129
 DmesonsFlowAnalysis.C:130
 DmesonsFlowAnalysis.C:131
 DmesonsFlowAnalysis.C:132
 DmesonsFlowAnalysis.C:133
 DmesonsFlowAnalysis.C:134
 DmesonsFlowAnalysis.C:135
 DmesonsFlowAnalysis.C:136
 DmesonsFlowAnalysis.C:137
 DmesonsFlowAnalysis.C:138
 DmesonsFlowAnalysis.C:139
 DmesonsFlowAnalysis.C:140
 DmesonsFlowAnalysis.C:141
 DmesonsFlowAnalysis.C:142
 DmesonsFlowAnalysis.C:143
 DmesonsFlowAnalysis.C:144
 DmesonsFlowAnalysis.C:145
 DmesonsFlowAnalysis.C:146
 DmesonsFlowAnalysis.C:147
 DmesonsFlowAnalysis.C:148
 DmesonsFlowAnalysis.C:149
 DmesonsFlowAnalysis.C:150
 DmesonsFlowAnalysis.C:151
 DmesonsFlowAnalysis.C:152
 DmesonsFlowAnalysis.C:153
 DmesonsFlowAnalysis.C:154
 DmesonsFlowAnalysis.C:155
 DmesonsFlowAnalysis.C:156
 DmesonsFlowAnalysis.C:157
 DmesonsFlowAnalysis.C:158
 DmesonsFlowAnalysis.C:159
 DmesonsFlowAnalysis.C:160
 DmesonsFlowAnalysis.C:161
 DmesonsFlowAnalysis.C:162
 DmesonsFlowAnalysis.C:163
 DmesonsFlowAnalysis.C:164
 DmesonsFlowAnalysis.C:165
 DmesonsFlowAnalysis.C:166
 DmesonsFlowAnalysis.C:167
 DmesonsFlowAnalysis.C:168
 DmesonsFlowAnalysis.C:169
 DmesonsFlowAnalysis.C:170
 DmesonsFlowAnalysis.C:171
 DmesonsFlowAnalysis.C:172
 DmesonsFlowAnalysis.C:173
 DmesonsFlowAnalysis.C:174
 DmesonsFlowAnalysis.C:175
 DmesonsFlowAnalysis.C:176
 DmesonsFlowAnalysis.C:177
 DmesonsFlowAnalysis.C:178
 DmesonsFlowAnalysis.C:179
 DmesonsFlowAnalysis.C:180
 DmesonsFlowAnalysis.C:181
 DmesonsFlowAnalysis.C:182
 DmesonsFlowAnalysis.C:183
 DmesonsFlowAnalysis.C:184
 DmesonsFlowAnalysis.C:185
 DmesonsFlowAnalysis.C:186
 DmesonsFlowAnalysis.C:187
 DmesonsFlowAnalysis.C:188
 DmesonsFlowAnalysis.C:189
 DmesonsFlowAnalysis.C:190
 DmesonsFlowAnalysis.C:191
 DmesonsFlowAnalysis.C:192
 DmesonsFlowAnalysis.C:193
 DmesonsFlowAnalysis.C:194
 DmesonsFlowAnalysis.C:195
 DmesonsFlowAnalysis.C:196
 DmesonsFlowAnalysis.C:197
 DmesonsFlowAnalysis.C:198
 DmesonsFlowAnalysis.C:199
 DmesonsFlowAnalysis.C:200
 DmesonsFlowAnalysis.C:201
 DmesonsFlowAnalysis.C:202
 DmesonsFlowAnalysis.C:203
 DmesonsFlowAnalysis.C:204
 DmesonsFlowAnalysis.C:205
 DmesonsFlowAnalysis.C:206
 DmesonsFlowAnalysis.C:207
 DmesonsFlowAnalysis.C:208
 DmesonsFlowAnalysis.C:209
 DmesonsFlowAnalysis.C:210
 DmesonsFlowAnalysis.C:211
 DmesonsFlowAnalysis.C:212
 DmesonsFlowAnalysis.C:213
 DmesonsFlowAnalysis.C:214
 DmesonsFlowAnalysis.C:215
 DmesonsFlowAnalysis.C:216
 DmesonsFlowAnalysis.C:217
 DmesonsFlowAnalysis.C:218
 DmesonsFlowAnalysis.C:219
 DmesonsFlowAnalysis.C:220
 DmesonsFlowAnalysis.C:221
 DmesonsFlowAnalysis.C:222
 DmesonsFlowAnalysis.C:223
 DmesonsFlowAnalysis.C:224
 DmesonsFlowAnalysis.C:225
 DmesonsFlowAnalysis.C:226
 DmesonsFlowAnalysis.C:227
 DmesonsFlowAnalysis.C:228
 DmesonsFlowAnalysis.C:229
 DmesonsFlowAnalysis.C:230
 DmesonsFlowAnalysis.C:231
 DmesonsFlowAnalysis.C:232
 DmesonsFlowAnalysis.C:233
 DmesonsFlowAnalysis.C:234
 DmesonsFlowAnalysis.C:235
 DmesonsFlowAnalysis.C:236
 DmesonsFlowAnalysis.C:237
 DmesonsFlowAnalysis.C:238
 DmesonsFlowAnalysis.C:239
 DmesonsFlowAnalysis.C:240
 DmesonsFlowAnalysis.C:241
 DmesonsFlowAnalysis.C:242
 DmesonsFlowAnalysis.C:243
 DmesonsFlowAnalysis.C:244
 DmesonsFlowAnalysis.C:245
 DmesonsFlowAnalysis.C:246
 DmesonsFlowAnalysis.C:247
 DmesonsFlowAnalysis.C:248
 DmesonsFlowAnalysis.C:249
 DmesonsFlowAnalysis.C:250
 DmesonsFlowAnalysis.C:251
 DmesonsFlowAnalysis.C:252
 DmesonsFlowAnalysis.C:253
 DmesonsFlowAnalysis.C:254
 DmesonsFlowAnalysis.C:255
 DmesonsFlowAnalysis.C:256
 DmesonsFlowAnalysis.C:257
 DmesonsFlowAnalysis.C:258
 DmesonsFlowAnalysis.C:259
 DmesonsFlowAnalysis.C:260
 DmesonsFlowAnalysis.C:261
 DmesonsFlowAnalysis.C:262
 DmesonsFlowAnalysis.C:263
 DmesonsFlowAnalysis.C:264
 DmesonsFlowAnalysis.C:265
 DmesonsFlowAnalysis.C:266
 DmesonsFlowAnalysis.C:267
 DmesonsFlowAnalysis.C:268
 DmesonsFlowAnalysis.C:269
 DmesonsFlowAnalysis.C:270
 DmesonsFlowAnalysis.C:271
 DmesonsFlowAnalysis.C:272
 DmesonsFlowAnalysis.C:273
 DmesonsFlowAnalysis.C:274
 DmesonsFlowAnalysis.C:275
 DmesonsFlowAnalysis.C:276
 DmesonsFlowAnalysis.C:277
 DmesonsFlowAnalysis.C:278
 DmesonsFlowAnalysis.C:279
 DmesonsFlowAnalysis.C:280
 DmesonsFlowAnalysis.C:281
 DmesonsFlowAnalysis.C:282
 DmesonsFlowAnalysis.C:283
 DmesonsFlowAnalysis.C:284
 DmesonsFlowAnalysis.C:285
 DmesonsFlowAnalysis.C:286
 DmesonsFlowAnalysis.C:287
 DmesonsFlowAnalysis.C:288
 DmesonsFlowAnalysis.C:289
 DmesonsFlowAnalysis.C:290
 DmesonsFlowAnalysis.C:291
 DmesonsFlowAnalysis.C:292
 DmesonsFlowAnalysis.C:293
 DmesonsFlowAnalysis.C:294
 DmesonsFlowAnalysis.C:295
 DmesonsFlowAnalysis.C:296
 DmesonsFlowAnalysis.C:297
 DmesonsFlowAnalysis.C:298
 DmesonsFlowAnalysis.C:299
 DmesonsFlowAnalysis.C:300
 DmesonsFlowAnalysis.C:301
 DmesonsFlowAnalysis.C:302
 DmesonsFlowAnalysis.C:303
 DmesonsFlowAnalysis.C:304
 DmesonsFlowAnalysis.C:305
 DmesonsFlowAnalysis.C:306
 DmesonsFlowAnalysis.C:307
 DmesonsFlowAnalysis.C:308
 DmesonsFlowAnalysis.C:309
 DmesonsFlowAnalysis.C:310
 DmesonsFlowAnalysis.C:311
 DmesonsFlowAnalysis.C:312
 DmesonsFlowAnalysis.C:313
 DmesonsFlowAnalysis.C:314
 DmesonsFlowAnalysis.C:315
 DmesonsFlowAnalysis.C:316
 DmesonsFlowAnalysis.C:317
 DmesonsFlowAnalysis.C:318
 DmesonsFlowAnalysis.C:319
 DmesonsFlowAnalysis.C:320
 DmesonsFlowAnalysis.C:321
 DmesonsFlowAnalysis.C:322
 DmesonsFlowAnalysis.C:323
 DmesonsFlowAnalysis.C:324
 DmesonsFlowAnalysis.C:325
 DmesonsFlowAnalysis.C:326
 DmesonsFlowAnalysis.C:327
 DmesonsFlowAnalysis.C:328
 DmesonsFlowAnalysis.C:329
 DmesonsFlowAnalysis.C:330
 DmesonsFlowAnalysis.C:331
 DmesonsFlowAnalysis.C:332
 DmesonsFlowAnalysis.C:333
 DmesonsFlowAnalysis.C:334
 DmesonsFlowAnalysis.C:335
 DmesonsFlowAnalysis.C:336
 DmesonsFlowAnalysis.C:337
 DmesonsFlowAnalysis.C:338
 DmesonsFlowAnalysis.C:339
 DmesonsFlowAnalysis.C:340
 DmesonsFlowAnalysis.C:341
 DmesonsFlowAnalysis.C:342
 DmesonsFlowAnalysis.C:343
 DmesonsFlowAnalysis.C:344
 DmesonsFlowAnalysis.C:345
 DmesonsFlowAnalysis.C:346
 DmesonsFlowAnalysis.C:347
 DmesonsFlowAnalysis.C:348
 DmesonsFlowAnalysis.C:349
 DmesonsFlowAnalysis.C:350
 DmesonsFlowAnalysis.C:351
 DmesonsFlowAnalysis.C:352
 DmesonsFlowAnalysis.C:353
 DmesonsFlowAnalysis.C:354
 DmesonsFlowAnalysis.C:355
 DmesonsFlowAnalysis.C:356
 DmesonsFlowAnalysis.C:357
 DmesonsFlowAnalysis.C:358
 DmesonsFlowAnalysis.C:359
 DmesonsFlowAnalysis.C:360
 DmesonsFlowAnalysis.C:361
 DmesonsFlowAnalysis.C:362
 DmesonsFlowAnalysis.C:363
 DmesonsFlowAnalysis.C:364
 DmesonsFlowAnalysis.C:365
 DmesonsFlowAnalysis.C:366
 DmesonsFlowAnalysis.C:367
 DmesonsFlowAnalysis.C:368
 DmesonsFlowAnalysis.C:369
 DmesonsFlowAnalysis.C:370
 DmesonsFlowAnalysis.C:371
 DmesonsFlowAnalysis.C:372
 DmesonsFlowAnalysis.C:373
 DmesonsFlowAnalysis.C:374
 DmesonsFlowAnalysis.C:375
 DmesonsFlowAnalysis.C:376
 DmesonsFlowAnalysis.C:377
 DmesonsFlowAnalysis.C:378
 DmesonsFlowAnalysis.C:379
 DmesonsFlowAnalysis.C:380
 DmesonsFlowAnalysis.C:381
 DmesonsFlowAnalysis.C:382
 DmesonsFlowAnalysis.C:383
 DmesonsFlowAnalysis.C:384
 DmesonsFlowAnalysis.C:385
 DmesonsFlowAnalysis.C:386
 DmesonsFlowAnalysis.C:387
 DmesonsFlowAnalysis.C:388
 DmesonsFlowAnalysis.C:389
 DmesonsFlowAnalysis.C:390
 DmesonsFlowAnalysis.C:391
 DmesonsFlowAnalysis.C:392
 DmesonsFlowAnalysis.C:393
 DmesonsFlowAnalysis.C:394
 DmesonsFlowAnalysis.C:395
 DmesonsFlowAnalysis.C:396
 DmesonsFlowAnalysis.C:397
 DmesonsFlowAnalysis.C:398
 DmesonsFlowAnalysis.C:399
 DmesonsFlowAnalysis.C:400
 DmesonsFlowAnalysis.C:401
 DmesonsFlowAnalysis.C:402
 DmesonsFlowAnalysis.C:403
 DmesonsFlowAnalysis.C:404
 DmesonsFlowAnalysis.C:405
 DmesonsFlowAnalysis.C:406
 DmesonsFlowAnalysis.C:407
 DmesonsFlowAnalysis.C:408
 DmesonsFlowAnalysis.C:409
 DmesonsFlowAnalysis.C:410
 DmesonsFlowAnalysis.C:411
 DmesonsFlowAnalysis.C:412
 DmesonsFlowAnalysis.C:413
 DmesonsFlowAnalysis.C:414
 DmesonsFlowAnalysis.C:415
 DmesonsFlowAnalysis.C:416
 DmesonsFlowAnalysis.C:417
 DmesonsFlowAnalysis.C:418
 DmesonsFlowAnalysis.C:419
 DmesonsFlowAnalysis.C:420
 DmesonsFlowAnalysis.C:421
 DmesonsFlowAnalysis.C:422
 DmesonsFlowAnalysis.C:423
 DmesonsFlowAnalysis.C:424
 DmesonsFlowAnalysis.C:425
 DmesonsFlowAnalysis.C:426
 DmesonsFlowAnalysis.C:427
 DmesonsFlowAnalysis.C:428
 DmesonsFlowAnalysis.C:429
 DmesonsFlowAnalysis.C:430
 DmesonsFlowAnalysis.C:431
 DmesonsFlowAnalysis.C:432
 DmesonsFlowAnalysis.C:433
 DmesonsFlowAnalysis.C:434
 DmesonsFlowAnalysis.C:435
 DmesonsFlowAnalysis.C:436
 DmesonsFlowAnalysis.C:437
 DmesonsFlowAnalysis.C:438
 DmesonsFlowAnalysis.C:439
 DmesonsFlowAnalysis.C:440
 DmesonsFlowAnalysis.C:441
 DmesonsFlowAnalysis.C:442
 DmesonsFlowAnalysis.C:443
 DmesonsFlowAnalysis.C:444
 DmesonsFlowAnalysis.C:445
 DmesonsFlowAnalysis.C:446
 DmesonsFlowAnalysis.C:447
 DmesonsFlowAnalysis.C:448
 DmesonsFlowAnalysis.C:449
 DmesonsFlowAnalysis.C:450
 DmesonsFlowAnalysis.C:451
 DmesonsFlowAnalysis.C:452
 DmesonsFlowAnalysis.C:453
 DmesonsFlowAnalysis.C:454
 DmesonsFlowAnalysis.C:455
 DmesonsFlowAnalysis.C:456
 DmesonsFlowAnalysis.C:457
 DmesonsFlowAnalysis.C:458
 DmesonsFlowAnalysis.C:459
 DmesonsFlowAnalysis.C:460
 DmesonsFlowAnalysis.C:461
 DmesonsFlowAnalysis.C:462
 DmesonsFlowAnalysis.C:463
 DmesonsFlowAnalysis.C:464
 DmesonsFlowAnalysis.C:465
 DmesonsFlowAnalysis.C:466
 DmesonsFlowAnalysis.C:467
 DmesonsFlowAnalysis.C:468
 DmesonsFlowAnalysis.C:469
 DmesonsFlowAnalysis.C:470
 DmesonsFlowAnalysis.C:471
 DmesonsFlowAnalysis.C:472
 DmesonsFlowAnalysis.C:473
 DmesonsFlowAnalysis.C:474
 DmesonsFlowAnalysis.C:475
 DmesonsFlowAnalysis.C:476
 DmesonsFlowAnalysis.C:477
 DmesonsFlowAnalysis.C:478
 DmesonsFlowAnalysis.C:479
 DmesonsFlowAnalysis.C:480
 DmesonsFlowAnalysis.C:481
 DmesonsFlowAnalysis.C:482
 DmesonsFlowAnalysis.C:483
 DmesonsFlowAnalysis.C:484
 DmesonsFlowAnalysis.C:485
 DmesonsFlowAnalysis.C:486
 DmesonsFlowAnalysis.C:487
 DmesonsFlowAnalysis.C:488
 DmesonsFlowAnalysis.C:489
 DmesonsFlowAnalysis.C:490
 DmesonsFlowAnalysis.C:491
 DmesonsFlowAnalysis.C:492
 DmesonsFlowAnalysis.C:493
 DmesonsFlowAnalysis.C:494
 DmesonsFlowAnalysis.C:495
 DmesonsFlowAnalysis.C:496
 DmesonsFlowAnalysis.C:497
 DmesonsFlowAnalysis.C:498
 DmesonsFlowAnalysis.C:499
 DmesonsFlowAnalysis.C:500
 DmesonsFlowAnalysis.C:501
 DmesonsFlowAnalysis.C:502
 DmesonsFlowAnalysis.C:503
 DmesonsFlowAnalysis.C:504
 DmesonsFlowAnalysis.C:505
 DmesonsFlowAnalysis.C:506
 DmesonsFlowAnalysis.C:507
 DmesonsFlowAnalysis.C:508
 DmesonsFlowAnalysis.C:509
 DmesonsFlowAnalysis.C:510
 DmesonsFlowAnalysis.C:511
 DmesonsFlowAnalysis.C:512
 DmesonsFlowAnalysis.C:513
 DmesonsFlowAnalysis.C:514
 DmesonsFlowAnalysis.C:515
 DmesonsFlowAnalysis.C:516
 DmesonsFlowAnalysis.C:517
 DmesonsFlowAnalysis.C:518
 DmesonsFlowAnalysis.C:519
 DmesonsFlowAnalysis.C:520
 DmesonsFlowAnalysis.C:521
 DmesonsFlowAnalysis.C:522
 DmesonsFlowAnalysis.C:523
 DmesonsFlowAnalysis.C:524
 DmesonsFlowAnalysis.C:525
 DmesonsFlowAnalysis.C:526
 DmesonsFlowAnalysis.C:527
 DmesonsFlowAnalysis.C:528
 DmesonsFlowAnalysis.C:529
 DmesonsFlowAnalysis.C:530
 DmesonsFlowAnalysis.C:531
 DmesonsFlowAnalysis.C:532
 DmesonsFlowAnalysis.C:533
 DmesonsFlowAnalysis.C:534
 DmesonsFlowAnalysis.C:535
 DmesonsFlowAnalysis.C:536
 DmesonsFlowAnalysis.C:537
 DmesonsFlowAnalysis.C:538
 DmesonsFlowAnalysis.C:539
 DmesonsFlowAnalysis.C:540
 DmesonsFlowAnalysis.C:541
 DmesonsFlowAnalysis.C:542
 DmesonsFlowAnalysis.C:543
 DmesonsFlowAnalysis.C:544
 DmesonsFlowAnalysis.C:545
 DmesonsFlowAnalysis.C:546
 DmesonsFlowAnalysis.C:547
 DmesonsFlowAnalysis.C:548
 DmesonsFlowAnalysis.C:549
 DmesonsFlowAnalysis.C:550
 DmesonsFlowAnalysis.C:551
 DmesonsFlowAnalysis.C:552
 DmesonsFlowAnalysis.C:553
 DmesonsFlowAnalysis.C:554
 DmesonsFlowAnalysis.C:555
 DmesonsFlowAnalysis.C:556
 DmesonsFlowAnalysis.C:557
 DmesonsFlowAnalysis.C:558
 DmesonsFlowAnalysis.C:559
 DmesonsFlowAnalysis.C:560
 DmesonsFlowAnalysis.C:561
 DmesonsFlowAnalysis.C:562
 DmesonsFlowAnalysis.C:563
 DmesonsFlowAnalysis.C:564
 DmesonsFlowAnalysis.C:565
 DmesonsFlowAnalysis.C:566
 DmesonsFlowAnalysis.C:567
 DmesonsFlowAnalysis.C:568
 DmesonsFlowAnalysis.C:569
 DmesonsFlowAnalysis.C:570
 DmesonsFlowAnalysis.C:571
 DmesonsFlowAnalysis.C:572
 DmesonsFlowAnalysis.C:573
 DmesonsFlowAnalysis.C:574
 DmesonsFlowAnalysis.C:575
 DmesonsFlowAnalysis.C:576
 DmesonsFlowAnalysis.C:577
 DmesonsFlowAnalysis.C:578
 DmesonsFlowAnalysis.C:579
 DmesonsFlowAnalysis.C:580
 DmesonsFlowAnalysis.C:581
 DmesonsFlowAnalysis.C:582
 DmesonsFlowAnalysis.C:583
 DmesonsFlowAnalysis.C:584
 DmesonsFlowAnalysis.C:585
 DmesonsFlowAnalysis.C:586
 DmesonsFlowAnalysis.C:587
 DmesonsFlowAnalysis.C:588
 DmesonsFlowAnalysis.C:589
 DmesonsFlowAnalysis.C:590
 DmesonsFlowAnalysis.C:591
 DmesonsFlowAnalysis.C:592
 DmesonsFlowAnalysis.C:593
 DmesonsFlowAnalysis.C:594
 DmesonsFlowAnalysis.C:595
 DmesonsFlowAnalysis.C:596
 DmesonsFlowAnalysis.C:597
 DmesonsFlowAnalysis.C:598
 DmesonsFlowAnalysis.C:599
 DmesonsFlowAnalysis.C:600
 DmesonsFlowAnalysis.C:601
 DmesonsFlowAnalysis.C:602
 DmesonsFlowAnalysis.C:603
 DmesonsFlowAnalysis.C:604
 DmesonsFlowAnalysis.C:605
 DmesonsFlowAnalysis.C:606
 DmesonsFlowAnalysis.C:607
 DmesonsFlowAnalysis.C:608
 DmesonsFlowAnalysis.C:609
 DmesonsFlowAnalysis.C:610
 DmesonsFlowAnalysis.C:611
 DmesonsFlowAnalysis.C:612
 DmesonsFlowAnalysis.C:613
 DmesonsFlowAnalysis.C:614
 DmesonsFlowAnalysis.C:615
 DmesonsFlowAnalysis.C:616
 DmesonsFlowAnalysis.C:617
 DmesonsFlowAnalysis.C:618
 DmesonsFlowAnalysis.C:619
 DmesonsFlowAnalysis.C:620
 DmesonsFlowAnalysis.C:621
 DmesonsFlowAnalysis.C:622
 DmesonsFlowAnalysis.C:623
 DmesonsFlowAnalysis.C:624
 DmesonsFlowAnalysis.C:625
 DmesonsFlowAnalysis.C:626
 DmesonsFlowAnalysis.C:627
 DmesonsFlowAnalysis.C:628
 DmesonsFlowAnalysis.C:629
 DmesonsFlowAnalysis.C:630
 DmesonsFlowAnalysis.C:631
 DmesonsFlowAnalysis.C:632
 DmesonsFlowAnalysis.C:633
 DmesonsFlowAnalysis.C:634
 DmesonsFlowAnalysis.C:635
 DmesonsFlowAnalysis.C:636
 DmesonsFlowAnalysis.C:637
 DmesonsFlowAnalysis.C:638
 DmesonsFlowAnalysis.C:639
 DmesonsFlowAnalysis.C:640
 DmesonsFlowAnalysis.C:641
 DmesonsFlowAnalysis.C:642
 DmesonsFlowAnalysis.C:643
 DmesonsFlowAnalysis.C:644
 DmesonsFlowAnalysis.C:645
 DmesonsFlowAnalysis.C:646
 DmesonsFlowAnalysis.C:647
 DmesonsFlowAnalysis.C:648
 DmesonsFlowAnalysis.C:649
 DmesonsFlowAnalysis.C:650
 DmesonsFlowAnalysis.C:651
 DmesonsFlowAnalysis.C:652
 DmesonsFlowAnalysis.C:653
 DmesonsFlowAnalysis.C:654
 DmesonsFlowAnalysis.C:655
 DmesonsFlowAnalysis.C:656
 DmesonsFlowAnalysis.C:657
 DmesonsFlowAnalysis.C:658
 DmesonsFlowAnalysis.C:659
 DmesonsFlowAnalysis.C:660
 DmesonsFlowAnalysis.C:661
 DmesonsFlowAnalysis.C:662
 DmesonsFlowAnalysis.C:663
 DmesonsFlowAnalysis.C:664
 DmesonsFlowAnalysis.C:665
 DmesonsFlowAnalysis.C:666
 DmesonsFlowAnalysis.C:667
 DmesonsFlowAnalysis.C:668
 DmesonsFlowAnalysis.C:669
 DmesonsFlowAnalysis.C:670
 DmesonsFlowAnalysis.C:671
 DmesonsFlowAnalysis.C:672
 DmesonsFlowAnalysis.C:673
 DmesonsFlowAnalysis.C:674
 DmesonsFlowAnalysis.C:675
 DmesonsFlowAnalysis.C:676
 DmesonsFlowAnalysis.C:677
 DmesonsFlowAnalysis.C:678
 DmesonsFlowAnalysis.C:679
 DmesonsFlowAnalysis.C:680
 DmesonsFlowAnalysis.C:681
 DmesonsFlowAnalysis.C:682
 DmesonsFlowAnalysis.C:683
 DmesonsFlowAnalysis.C:684
 DmesonsFlowAnalysis.C:685
 DmesonsFlowAnalysis.C:686
 DmesonsFlowAnalysis.C:687
 DmesonsFlowAnalysis.C:688
 DmesonsFlowAnalysis.C:689
 DmesonsFlowAnalysis.C:690
 DmesonsFlowAnalysis.C:691
 DmesonsFlowAnalysis.C:692
 DmesonsFlowAnalysis.C:693
 DmesonsFlowAnalysis.C:694
 DmesonsFlowAnalysis.C:695
 DmesonsFlowAnalysis.C:696
 DmesonsFlowAnalysis.C:697
 DmesonsFlowAnalysis.C:698
 DmesonsFlowAnalysis.C:699
 DmesonsFlowAnalysis.C:700
 DmesonsFlowAnalysis.C:701
 DmesonsFlowAnalysis.C:702
 DmesonsFlowAnalysis.C:703
 DmesonsFlowAnalysis.C:704
 DmesonsFlowAnalysis.C:705
 DmesonsFlowAnalysis.C:706
 DmesonsFlowAnalysis.C:707
 DmesonsFlowAnalysis.C:708
 DmesonsFlowAnalysis.C:709
 DmesonsFlowAnalysis.C:710
 DmesonsFlowAnalysis.C:711
 DmesonsFlowAnalysis.C:712
 DmesonsFlowAnalysis.C:713
 DmesonsFlowAnalysis.C:714
 DmesonsFlowAnalysis.C:715
 DmesonsFlowAnalysis.C:716
 DmesonsFlowAnalysis.C:717
 DmesonsFlowAnalysis.C:718
 DmesonsFlowAnalysis.C:719
 DmesonsFlowAnalysis.C:720
 DmesonsFlowAnalysis.C:721
 DmesonsFlowAnalysis.C:722
 DmesonsFlowAnalysis.C:723
 DmesonsFlowAnalysis.C:724
 DmesonsFlowAnalysis.C:725
 DmesonsFlowAnalysis.C:726
 DmesonsFlowAnalysis.C:727
 DmesonsFlowAnalysis.C:728
 DmesonsFlowAnalysis.C:729
 DmesonsFlowAnalysis.C:730
 DmesonsFlowAnalysis.C:731
 DmesonsFlowAnalysis.C:732
 DmesonsFlowAnalysis.C:733
 DmesonsFlowAnalysis.C:734
 DmesonsFlowAnalysis.C:735
 DmesonsFlowAnalysis.C:736
 DmesonsFlowAnalysis.C:737
 DmesonsFlowAnalysis.C:738
 DmesonsFlowAnalysis.C:739
 DmesonsFlowAnalysis.C:740
 DmesonsFlowAnalysis.C:741
 DmesonsFlowAnalysis.C:742
 DmesonsFlowAnalysis.C:743
 DmesonsFlowAnalysis.C:744
 DmesonsFlowAnalysis.C:745
 DmesonsFlowAnalysis.C:746
 DmesonsFlowAnalysis.C:747
 DmesonsFlowAnalysis.C:748
 DmesonsFlowAnalysis.C:749
 DmesonsFlowAnalysis.C:750
 DmesonsFlowAnalysis.C:751
 DmesonsFlowAnalysis.C:752
 DmesonsFlowAnalysis.C:753
 DmesonsFlowAnalysis.C:754
 DmesonsFlowAnalysis.C:755
 DmesonsFlowAnalysis.C:756
 DmesonsFlowAnalysis.C:757
 DmesonsFlowAnalysis.C:758
 DmesonsFlowAnalysis.C:759
 DmesonsFlowAnalysis.C:760
 DmesonsFlowAnalysis.C:761
 DmesonsFlowAnalysis.C:762
 DmesonsFlowAnalysis.C:763
 DmesonsFlowAnalysis.C:764
 DmesonsFlowAnalysis.C:765
 DmesonsFlowAnalysis.C:766
 DmesonsFlowAnalysis.C:767
 DmesonsFlowAnalysis.C:768
 DmesonsFlowAnalysis.C:769
 DmesonsFlowAnalysis.C:770
 DmesonsFlowAnalysis.C:771
 DmesonsFlowAnalysis.C:772
 DmesonsFlowAnalysis.C:773
 DmesonsFlowAnalysis.C:774
 DmesonsFlowAnalysis.C:775
 DmesonsFlowAnalysis.C:776
 DmesonsFlowAnalysis.C:777
 DmesonsFlowAnalysis.C:778
 DmesonsFlowAnalysis.C:779
 DmesonsFlowAnalysis.C:780
 DmesonsFlowAnalysis.C:781
 DmesonsFlowAnalysis.C:782
 DmesonsFlowAnalysis.C:783
 DmesonsFlowAnalysis.C:784
 DmesonsFlowAnalysis.C:785
 DmesonsFlowAnalysis.C:786
 DmesonsFlowAnalysis.C:787
 DmesonsFlowAnalysis.C:788
 DmesonsFlowAnalysis.C:789
 DmesonsFlowAnalysis.C:790
 DmesonsFlowAnalysis.C:791
 DmesonsFlowAnalysis.C:792
 DmesonsFlowAnalysis.C:793
 DmesonsFlowAnalysis.C:794
 DmesonsFlowAnalysis.C:795
 DmesonsFlowAnalysis.C:796
 DmesonsFlowAnalysis.C:797
 DmesonsFlowAnalysis.C:798
 DmesonsFlowAnalysis.C:799
 DmesonsFlowAnalysis.C:800
 DmesonsFlowAnalysis.C:801
 DmesonsFlowAnalysis.C:802
 DmesonsFlowAnalysis.C:803
 DmesonsFlowAnalysis.C:804
 DmesonsFlowAnalysis.C:805
 DmesonsFlowAnalysis.C:806
 DmesonsFlowAnalysis.C:807
 DmesonsFlowAnalysis.C:808
 DmesonsFlowAnalysis.C:809
 DmesonsFlowAnalysis.C:810
 DmesonsFlowAnalysis.C:811
 DmesonsFlowAnalysis.C:812
 DmesonsFlowAnalysis.C:813
 DmesonsFlowAnalysis.C:814
 DmesonsFlowAnalysis.C:815
 DmesonsFlowAnalysis.C:816
 DmesonsFlowAnalysis.C:817
 DmesonsFlowAnalysis.C:818
 DmesonsFlowAnalysis.C:819
 DmesonsFlowAnalysis.C:820
 DmesonsFlowAnalysis.C:821
 DmesonsFlowAnalysis.C:822
 DmesonsFlowAnalysis.C:823
 DmesonsFlowAnalysis.C:824
 DmesonsFlowAnalysis.C:825
 DmesonsFlowAnalysis.C:826
 DmesonsFlowAnalysis.C:827
 DmesonsFlowAnalysis.C:828
 DmesonsFlowAnalysis.C:829
 DmesonsFlowAnalysis.C:830
 DmesonsFlowAnalysis.C:831
 DmesonsFlowAnalysis.C:832
 DmesonsFlowAnalysis.C:833
 DmesonsFlowAnalysis.C:834
 DmesonsFlowAnalysis.C:835
 DmesonsFlowAnalysis.C:836
 DmesonsFlowAnalysis.C:837
 DmesonsFlowAnalysis.C:838
 DmesonsFlowAnalysis.C:839
 DmesonsFlowAnalysis.C:840
 DmesonsFlowAnalysis.C:841
 DmesonsFlowAnalysis.C:842
 DmesonsFlowAnalysis.C:843
 DmesonsFlowAnalysis.C:844
 DmesonsFlowAnalysis.C:845
 DmesonsFlowAnalysis.C:846
 DmesonsFlowAnalysis.C:847
 DmesonsFlowAnalysis.C:848
 DmesonsFlowAnalysis.C:849
 DmesonsFlowAnalysis.C:850
 DmesonsFlowAnalysis.C:851
 DmesonsFlowAnalysis.C:852
 DmesonsFlowAnalysis.C:853
 DmesonsFlowAnalysis.C:854
 DmesonsFlowAnalysis.C:855
 DmesonsFlowAnalysis.C:856
 DmesonsFlowAnalysis.C:857
 DmesonsFlowAnalysis.C:858
 DmesonsFlowAnalysis.C:859
 DmesonsFlowAnalysis.C:860
 DmesonsFlowAnalysis.C:861
 DmesonsFlowAnalysis.C:862
 DmesonsFlowAnalysis.C:863
 DmesonsFlowAnalysis.C:864
 DmesonsFlowAnalysis.C:865
 DmesonsFlowAnalysis.C:866
 DmesonsFlowAnalysis.C:867
 DmesonsFlowAnalysis.C:868
 DmesonsFlowAnalysis.C:869
 DmesonsFlowAnalysis.C:870
 DmesonsFlowAnalysis.C:871
 DmesonsFlowAnalysis.C:872
 DmesonsFlowAnalysis.C:873
 DmesonsFlowAnalysis.C:874
 DmesonsFlowAnalysis.C:875
 DmesonsFlowAnalysis.C:876
 DmesonsFlowAnalysis.C:877
 DmesonsFlowAnalysis.C:878
 DmesonsFlowAnalysis.C:879
 DmesonsFlowAnalysis.C:880
 DmesonsFlowAnalysis.C:881
 DmesonsFlowAnalysis.C:882
 DmesonsFlowAnalysis.C:883
 DmesonsFlowAnalysis.C:884
 DmesonsFlowAnalysis.C:885
 DmesonsFlowAnalysis.C:886
 DmesonsFlowAnalysis.C:887
 DmesonsFlowAnalysis.C:888
 DmesonsFlowAnalysis.C:889
 DmesonsFlowAnalysis.C:890
 DmesonsFlowAnalysis.C:891
 DmesonsFlowAnalysis.C:892
 DmesonsFlowAnalysis.C:893
 DmesonsFlowAnalysis.C:894
 DmesonsFlowAnalysis.C:895
 DmesonsFlowAnalysis.C:896
 DmesonsFlowAnalysis.C:897
 DmesonsFlowAnalysis.C:898
 DmesonsFlowAnalysis.C:899
 DmesonsFlowAnalysis.C:900
 DmesonsFlowAnalysis.C:901
 DmesonsFlowAnalysis.C:902
 DmesonsFlowAnalysis.C:903
 DmesonsFlowAnalysis.C:904