const Int_t numberOfCentralityBins = 10;
TString centralityArray[numberOfCentralityBins] = {"0-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-1","1-2"};

const Int_t gRebin = 1;
void drawBalanceFunction2DPsi(const char* filename = "AnalysisResultsPsi.root", 
			      Int_t gCentrality = 1,
			      Int_t gBit = -1,
			      const char* gCentralityEstimator = 0x0,
			      Bool_t kShowShuffled = kFALSE, 
			      Bool_t kShowMixed = kTRUE, 
			      Double_t psiMin = -0.5, Double_t psiMax = 0.5,
			      Double_t ptTriggerMin = -1.,
			      Double_t ptTriggerMax = -1.,
			      Double_t ptAssociatedMin = -1.,
			      Double_t ptAssociatedMax = -1.,
			      Bool_t k2pMethod = kFALSE) {
  //Macro that draws the BF distributions for each centrality bin
  //for reaction plane dependent analysis
  //Author: Panos.Christakoglou@nikhef.nl
  //Load the PWG2ebye library
  gSystem->Load("libANALYSIS.so");
  gSystem->Load("libANALYSISalice.so");
  gSystem->Load("libEventMixing.so");
  gSystem->Load("libCORRFW.so");
  gSystem->Load("libPWGTools.so");
  gSystem->Load("libPWGCFebye.so");

  gROOT->LoadMacro("~/SetPlotStyle.C");
  SetPlotStyle();
  gStyle->SetPalette(1,0);

  //Prepare the objects and return them
  TList *listBF = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0);
  TList *listBFShuffled = NULL;
  if(kShowShuffled) listBFShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1);
  TList *listBFMixed = NULL;
  if(kShowMixed) listBFMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2);
  if(!listBF) {
    Printf("The TList object was not created");
    return;
  }
  else 
    draw(listBF,listBFShuffled,listBFMixed,gCentrality,psiMin,psiMax,
	 ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
	 k2pMethod);  
}

//______________________________________________________//
TList *GetListOfObjects(const char* filename,
			Int_t gCentrality,
			Int_t gBit,
			const char *gCentralityEstimator,
			Int_t kData = 1) {
  //Get the TList objects (QA, bf, bf shuffled)
  TList *listBF = 0x0;
  
  //Open the file
  TFile *f = TFile::Open(filename,"UPDATE");
  if((!f)||(!f->IsOpen())) {
    Printf("The file %s is not found. Aborting...",filename);
    return listBF;
  }
  //f->ls();
  
  TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
  if(!dir) {   
    Printf("The TDirectoryFile is not found. Aborting...",filename);
    return listBF;
  }
  //dir->ls();
  
  TString listBFName;
  if(kData == 0) {
    //cout<<"no shuffling - no mixing"<<endl;
    listBFName = "listBFPsi_";
  }
  else if(kData == 1) {
    //cout<<"shuffling - no mixing"<<endl;
    listBFName = "listBFPsiShuffled_";
  }
  else if(kData == 2) {
    //cout<<"no shuffling - mixing"<<endl;
    listBFName = "listBFPsiMixed_";
  }
  listBFName += centralityArray[gCentrality-1];
  if(gBit > -1) {
    listBFName += "_Bit"; listBFName += gBit; }
  if(gCentralityEstimator) {
    listBFName += "_"; listBFName += gCentralityEstimator;}

  // histograms were already retrieved (in first iteration)
  if(dir->Get(Form("%s_histograms",listBFName.Data()))){
    //cout<<"second iteration"<<endl;
    listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
  }

  // histograms were not yet retrieved (this is the first iteration)
  else{

    listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
    cout<<"======================================================="<<endl;
    cout<<"List name (check): "<<listBFName.Data()<<endl;
    cout<<"List name: "<<listBF->GetName()<<endl;
    //listBF->ls();
    
    //Get the histograms
    TString histoName;
    if(kData == 0)
      histoName = "fHistPV0M";
    else if(kData == 1)
      histoName = "fHistP_shuffleV0M";
    else if(kData == 2)
      histoName = "fHistPV0M";
    AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));  
    if(!fHistP) {
      Printf("fHistP %s not found!!!",histoName.Data());
      break;
    }
    fHistP->FillParent(); fHistP->DeleteContainers();
    
    if(kData == 0)
      histoName = "fHistNV0M";
    if(kData == 1)
      histoName = "fHistN_shuffleV0M";
    if(kData == 2)
      histoName = "fHistNV0M";
    AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
    if(!fHistN) {
      Printf("fHistN %s not found!!!",histoName.Data());
      break;
    }
    fHistN->FillParent(); fHistN->DeleteContainers();
    
    if(kData == 0)
      histoName = "fHistPNV0M";
    if(kData == 1)
      histoName = "fHistPN_shuffleV0M";
    if(kData == 2)
      histoName = "fHistPNV0M";
    AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
    if(!fHistPN) {
      Printf("fHistPN %s not found!!!",histoName.Data());
      break;
    }
    fHistPN->FillParent(); fHistPN->DeleteContainers();
    
    if(kData == 0)
      histoName = "fHistNPV0M";
    if(kData == 1)
      histoName = "fHistNP_shuffleV0M";
    if(kData == 2)
      histoName = "fHistNPV0M";
    AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
    if(!fHistNP) {
      Printf("fHistNP %s not found!!!",histoName.Data());
      break;
    }
    fHistNP->FillParent(); fHistNP->DeleteContainers();
    
    if(kData == 0)
      histoName = "fHistPPV0M";
    if(kData == 1)
      histoName = "fHistPP_shuffleV0M";
    if(kData == 2)
      histoName = "fHistPPV0M";
    AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
    if(!fHistPP) {
      Printf("fHistPP %s not found!!!",histoName.Data());
      break;
    }
    fHistPP->FillParent(); fHistPP->DeleteContainers();
    
    if(kData == 0)
      histoName = "fHistNNV0M";
    if(kData == 1)
      histoName = "fHistNN_shuffleV0M";
    if(kData == 2)
      histoName = "fHistNNV0M";
    AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
    if(!fHistNN) {
      Printf("fHistNN %s not found!!!",histoName.Data());
      break;
    }
    fHistNN->FillParent(); fHistNN->DeleteContainers();
    
    dir->cd();
    listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
    
  }// first iteration
  
  f->Close();
  
  return listBF;
}

//______________________________________________________//
void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
	  Int_t gCentrality, Double_t psiMin, Double_t psiMax,
	  Double_t ptTriggerMin, Double_t ptTriggerMax,
	  Double_t ptAssociatedMin, Double_t ptAssociatedMax,
	  Bool_t k2pMethod = kFALSE) {  
  //balance function
  AliTHn *hP = NULL;
  AliTHn *hN = NULL;
  AliTHn *hPN = NULL;
  AliTHn *hNP = NULL;
  AliTHn *hPP = NULL;
  AliTHn *hNN = NULL;
  //listBF->ls();
  //Printf("=================");
  hP = (AliTHn*) listBF->FindObject("fHistPV0M");
  hP->SetName("gHistP");
  hN = (AliTHn*) listBF->FindObject("fHistNV0M");
  hN->SetName("gHistN");
  hPN = (AliTHn*) listBF->FindObject("fHistPNV0M");
  hPN->SetName("gHistPN");
  hNP = (AliTHn*) listBF->FindObject("fHistNPV0M");
  hNP->SetName("gHistNP");
  hPP = (AliTHn*) listBF->FindObject("fHistPPV0M");
  hPP->SetName("gHistPP");
  hNN = (AliTHn*) listBF->FindObject("fHistNNV0M");
  hNN->SetName("gHistNN");

  AliBalancePsi *b = new AliBalancePsi();
  b->SetHistNp(hP);
  b->SetHistNn(hN);
  b->SetHistNpn(hPN);
  b->SetHistNnp(hNP);
  b->SetHistNpp(hPP);
  b->SetHistNnn(hNN);

  //balance function shuffling
  AliTHn *hPShuffled = NULL;
  AliTHn *hNShuffled = NULL;
  AliTHn *hPNShuffled = NULL;
  AliTHn *hNPShuffled = NULL;
  AliTHn *hPPShuffled = NULL;
  AliTHn *hNNShuffled = NULL;
  if(listBFShuffled) {
    //listBFShuffled->ls();
    
    hPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistP_shuffleV0M");
    hPShuffled->SetName("gHistPShuffled");
    hNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistN_shuffleV0M");
    hNShuffled->SetName("gHistNShuffled");
    hPNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPN_shuffleV0M");
    hPNShuffled->SetName("gHistPNShuffled");
    hNPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNP_shuffleV0M");
    hNPShuffled->SetName("gHistNPShuffled");
    hPPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPP_shuffleV0M");
    hPPShuffled->SetName("gHistPPShuffled");
    hNNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNN_shuffleV0M");
    hNNShuffled->SetName("gHistNNShuffled");
    
    AliBalancePsi *bShuffled = new AliBalancePsi();
    bShuffled->SetHistNp(hPShuffled);
    bShuffled->SetHistNn(hNShuffled);
    bShuffled->SetHistNpn(hPNShuffled);
    bShuffled->SetHistNnp(hNPShuffled);
    bShuffled->SetHistNpp(hPPShuffled);
    bShuffled->SetHistNnn(hNNShuffled);
  }

  //balance function mixing
  AliTHn *hPMixed = NULL;
  AliTHn *hNMixed = NULL;
  AliTHn *hPNMixed = NULL;
  AliTHn *hNPMixed = NULL;
  AliTHn *hPPMixed = NULL;
  AliTHn *hNNMixed = NULL;

  if(listBFMixed) {
    //listBFMixed->ls();

    hPMixed = (AliTHn*) listBFMixed->FindObject("fHistPV0M");
    hPMixed->SetName("gHistPMixed");
    hNMixed = (AliTHn*) listBFMixed->FindObject("fHistNV0M");
    hNMixed->SetName("gHistNMixed");
    hPNMixed = (AliTHn*) listBFMixed->FindObject("fHistPNV0M");
    hPNMixed->SetName("gHistPNMixed");
    hNPMixed = (AliTHn*) listBFMixed->FindObject("fHistNPV0M");
    hNPMixed->SetName("gHistNPMixed");
    hPPMixed = (AliTHn*) listBFMixed->FindObject("fHistPPV0M");
    hPPMixed->SetName("gHistPPMixed");
    hNNMixed = (AliTHn*) listBFMixed->FindObject("fHistNNV0M");
    hNNMixed->SetName("gHistNNMixed");
    
    AliBalancePsi *bMixed = new AliBalancePsi();
    bMixed->SetHistNp(hPMixed);
    bMixed->SetHistNn(hNMixed);
    bMixed->SetHistNpn(hPNMixed);
    bMixed->SetHistNnp(hNPMixed);
    bMixed->SetHistNpp(hPPMixed);
    bMixed->SetHistNnn(hNNMixed);
  }

  TH2D *gHistBalanceFunction;
  TH2D *gHistBalanceFunctionSubtracted;
  TH2D *gHistBalanceFunctionShuffled;
  TH2D *gHistBalanceFunctionMixed;
  TString histoTitle, pngName;
  
  histoTitle = "Centrality: "; 
  histoTitle += centralityArray[gCentrality-1]; 
  histoTitle += "%";
  if((psiMin == -0.5)&&(psiMax == 0.5))
    histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
  else if((psiMin == 0.5)&&(psiMax == 1.5))
    histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
  else if((psiMin == 1.5)&&(psiMax == 2.5))
    histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
  else 
    histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
  
  if(k2pMethod) 
    if(bMixed)
      gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
    else{
      cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
      return;
    }
  else
    gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
  gHistBalanceFunction->SetTitle(histoTitle.Data());
  gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
  gHistBalanceFunction->SetName("gHistBalanceFunction");

  if(listBFShuffled) {
    
    if(k2pMethod) 
      if(bMixed)
	gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
      else{
	cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
	return;
      }
    else
      gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
    gHistBalanceFunctionShuffled->SetTitle(histoTitle.Data());
    gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
    gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
  }

  if(listBFMixed) {
    if(k2pMethod) 
      if(bMixed)
	gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
      else{
	cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
	return;
      }
    else
      gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
    gHistBalanceFunctionMixed->SetTitle(histoTitle.Data());
    gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
    gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
  
    gHistBalanceFunctionSubtracted = dynamic_cast<TH2D *>(gHistBalanceFunction->Clone());
    gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
    gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
    gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
    gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
  }

  //Draw the results
  TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
  c1->SetFillColor(10); 
  c1->SetHighLightColor(10);
  c1->SetLeftMargin(0.15);
  gHistBalanceFunction->DrawCopy("lego2");
  gPad->SetTheta(30); // default is 30
  gPad->SetPhi(-60); // default is 30
  gPad->Update();  
  TCanvas *c1a = new TCanvas("c1a","",600,0,600,500);
  c1a->SetFillColor(10); 
  c1a->SetHighLightColor(10);
  c1a->SetLeftMargin(0.15);
  gHistBalanceFunction->DrawCopy("colz");

  if(listBFShuffled) {
    TCanvas *c2 = new TCanvas("c2","",100,100,600,500);
    c2->SetFillColor(10); 
    c2->SetHighLightColor(10);
    c2->SetLeftMargin(0.15);
    gHistBalanceFunctionShuffled->DrawCopy("lego2");
    gPad->SetTheta(30); // default is 30
    gPad->SetPhi(-60); // default is 30
    gPad->Update();  
    TCanvas *c2a = new TCanvas("c2a","",700,100,600,500);
    c2a->SetFillColor(10); 
    c2a->SetHighLightColor(10);
    c2a->SetLeftMargin(0.15);
    gHistBalanceFunctionShuffled->DrawCopy("colz");
  }

  if(listBFMixed) {
    TCanvas *c3 = new TCanvas("c3","",200,200,600,500);
    c3->SetFillColor(10); 
    c3->SetHighLightColor(10);
    c3->SetLeftMargin(0.15);
    gHistBalanceFunctionMixed->DrawCopy("lego2");
    gPad->SetTheta(30); // default is 30
    gPad->SetPhi(-60); // default is 30
    gPad->Update();  
    TCanvas *c3a = new TCanvas("c3a","",800,200,600,500);
    c3a->SetFillColor(10); 
    c3a->SetHighLightColor(10);
    c3a->SetLeftMargin(0.15);
    gHistBalanceFunctionMixed->DrawCopy("colz");

    TCanvas *c4 = new TCanvas("c4","",300,300,600,500);
    c4->SetFillColor(10); 
    c4->SetHighLightColor(10);
    c4->SetLeftMargin(0.15);
    gHistBalanceFunctionSubtracted->DrawCopy("lego2");
    gPad->SetTheta(30); // default is 30
    gPad->SetPhi(-60); // default is 30
    gPad->Update();  
    TCanvas *c4a = new TCanvas("c4a","",900,300,600,500);
    c4a->SetFillColor(10); 
    c4a->SetHighLightColor(10);
    c4a->SetLeftMargin(0.15);
    gHistBalanceFunctionSubtracted->DrawCopy("colz");

    fitbalanceFunction(gCentrality, psiMin = -0.5, psiMax,
		       ptTriggerMin, ptTriggerMax,
		       ptAssociatedMin, ptAssociatedMax,
		       gHistBalanceFunctionSubtracted);
  }

  TString newFileName = "balanceFunction2D.Centrality"; 
  newFileName += gCentrality; newFileName += ".Psi";
  if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
  else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
  else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
  else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
  else newFileName += "All.PttFrom";
  newFileName += Form("%.1f",ptTriggerMin); newFileName += "To"; 
  newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
  newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To"; 
  newFileName += Form("%.1f",ptAssociatedMax); 
  if(k2pMethod) newFileName += "_2pMethod";
  newFileName += ".root";

  TFile *fOutput = new TFile(newFileName.Data(),"recreate");
  fOutput->cd();
  /*hP->Write(); hN->Write();
  hPN->Write(); hNP->Write();
  hPP->Write(); hNN->Write();
  hPShuffled->Write(); hNShuffled->Write();
  hPNShuffled->Write(); hNPShuffled->Write();
  hPPShuffled->Write(); hNNShuffled->Write();
  hPMixed->Write(); hNMixed->Write();
  hPNMixed->Write(); hNPMixed->Write();
  hPPMixed->Write(); hNNMixed->Write();*/
  gHistBalanceFunction->Write();
  if(listBFShuffled) gHistBalanceFunctionShuffled->Write();
  if(listBFMixed) {
    gHistBalanceFunctionMixed->Write();
    gHistBalanceFunctionSubtracted->Write();
  }
  fOutput->Close();
}

//____________________________________________________________//
void fitbalanceFunction(Int_t gCentrality = 1,
			Double_t psiMin = -0.5, Double_t psiMax = 3.5,
			Double_t ptTriggerMin = -1.,
			Double_t ptTriggerMax = -1.,
			Double_t ptAssociatedMin = -1.,
			Double_t ptAssociatedMax = -1.,
			TH2D *gHist) {

  cout<<"FITTING FUNCTION"<<endl;

  //near side peak: [1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))
  //away side ridge: [5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))
  //longitudinal ridge: [8]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[9]),2)),[10]))
  //wing structures: [11]*TMath::Power(x,2)
  //flow contribution (v1 up to v4): 2.*([12]*TMath::Cos(y) + [13]*TMath::Cos(2.*y) + [14]*TMath::Cos(3.*y) + [15]*TMath::Cos(4.*y))
  TF2 *gFitFunction = new TF2("gFitFunction","[0]+[1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))+[5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))+[8]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[17])/[9]),2)),[10]))+[8]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[17])/[9]),2)),[10]))+[11]*TMath::Power(x,2)+2.*[12]*([13]*TMath::Cos(y) + [14]*TMath::Cos(2.*y) + [15]*TMath::Cos(3.*y) + [16]*TMath::Cos(4.*y))",-2.0,2.0,-TMath::Pi()/2.,3.*TMath::Pi()/2.); 
  gFitFunction->SetName("gFitFunction");


  //Normalization
  gFitFunction->SetParName(0,"N1"); 
  //near side peak
  gFitFunction->SetParName(1,"N_{near side}");gFitFunction->FixParameter(1,0.0);
  gFitFunction->SetParName(2,"Sigma_{near side}(delta eta)"); gFitFunction->FixParameter(2,0.0);
  gFitFunction->SetParName(3,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(3,0.0);
  gFitFunction->SetParName(4,"Exponent_{near side}"); gFitFunction->FixParameter(4,1.0);
  //away side ridge
  gFitFunction->SetParName(5,"N_{away side}"); gFitFunction->FixParameter(5,0.0);
  gFitFunction->SetParName(6,"Sigma_{away side}(delta phi)"); gFitFunction->FixParameter(6,0.0);
  gFitFunction->SetParName(7,"Exponent_{away side}"); gFitFunction->FixParameter(7,1.0);
  //longitudinal ridge
  gFitFunction->SetParName(8,"N_{long. ridge}"); gFitFunction->SetParameter(8,0.05);//
  gFitFunction->FixParameter(8,0.0);
  gFitFunction->SetParName(9,"Sigma_{long. ridge}(delta eta)"); gFitFunction->FixParameter(9,0.0);
  gFitFunction->SetParName(10,"Exponent_{long. ridge}"); gFitFunction->FixParameter(10,1.0);
  //wing structures
  gFitFunction->SetParName(11,"N_{wing}"); gFitFunction->FixParameter(11,0.0);
  //flow contribution
  gFitFunction->SetParName(12,"N_{flow}"); gFitFunction->FixParameter(12,0.0);
  gFitFunction->SetParName(13,"V1"); gFitFunction->FixParameter(13,0.0);
  gFitFunction->SetParName(14,"V2"); gFitFunction->FixParameter(14,0.0);
  gFitFunction->SetParName(15,"V3"); gFitFunction->FixParameter(15,0.0);
  gFitFunction->SetParName(16,"V4"); gFitFunction->FixParameter(16,0.0);
  gFitFunction->SetParName(17,"Offset"); gFitFunction->FixParameter(17,0.0);

  // just release the near side peak
  gFitFunction->ReleaseParameter(0);gFitFunction->SetParameter(0,1.0);
  gFitFunction->ReleaseParameter(1);gFitFunction->SetParameter(1,0.3);gFitFunction->SetParLimits(1,0.0,100);
  gFitFunction->ReleaseParameter(2);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
  gFitFunction->ReleaseParameter(3);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);

  gHist->Fit("gFitFunction","nm");


  //Cloning the histogram
  TString histoName = gHist->GetName(); histoName += "Fit"; 
  TH2D *gHistFit = new TH2D(histoName.Data(),";#Delta#eta;#Delta#varphi (rad);C(#Delta#eta,#Delta#varphi)",gHist->GetNbinsX(),gHist->GetXaxis()->GetXmin(),gHist->GetXaxis()->GetXmax(),gHist->GetNbinsY(),gHist->GetYaxis()->GetXmin(),gHist->GetYaxis()->GetXmax());
  TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
  gHistResidual->SetName("gHistResidual");
  gHistResidual->Sumw2();

  for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
    for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
      gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,gFitFunction->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
    }
  }
  gHistResidual->Add(gHistFit,-1);

  //Write to output file
  TString newFileName = "balanceFunctionFit";
  newFileName += ".Centrality";  
  newFileName += gCentrality; newFileName += ".Psi";
  if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
  else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
  else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
  else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
  else newFileName += "All.PttFrom";
  newFileName += Form("%.1f",ptTriggerMin); newFileName += "To"; 
  newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
  newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To"; 
  newFileName += Form("%.1f",ptAssociatedMax); 
  newFileName += ".root";
  TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
  gHist->Write();
  gHistFit->Write();
  gHistResidual->Write();
  gFitFunction->Write();
  newFile->Close();
  

}



//____________________________________________________________//
void drawBFPsi2D(const char* lhcPeriod = "LHC11h",
		 Int_t gCentrality = 1,
		 Bool_t kShowShuffled = kFALSE, 
		 Bool_t kShowMixed = kFALSE, 
		 Double_t psiMin = -0.5, Double_t psiMax = 0.5,
		 Double_t ptTriggerMin = -1.,
		 Double_t ptTriggerMax = -1.,
		 Double_t ptAssociatedMin = -1.,
		 Double_t ptAssociatedMax = -1.) {
  //Macro that draws the BF distributions for each centrality bin
  //for reaction plane dependent analysis
  //Author: Panos.Christakoglou@nikhef.nl
  TGaxis::SetMaxDigits(3);

  //Get the input file
  TString filename = lhcPeriod; filename +="/PttFrom";
  filename += Form("%.1f",ptTriggerMin); filename += "To"; 
  filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
  filename += Form("%.1f",ptAssociatedMin); filename += "To"; 
  filename += Form("%.1f",ptAssociatedMax); 
  filename += "/balanceFunction2D.Centrality"; 
  filename += gCentrality; filename += ".Psi";
  if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
  else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
  else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
  else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
  else filename += "All.PttFrom";
  filename += Form("%.1f",ptTriggerMin); filename += "To"; 
  filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
  filename += Form("%.1f",ptAssociatedMin); filename += "To"; 
  filename += Form("%.1f",ptAssociatedMax);  filename += ".root";  

  //Open the file
  TFile *f = TFile::Open(filename.Data());
  if((!f)||(!f->IsOpen())) {
    Printf("The file %s is not found. Aborting...",filename);
    return listBF;
  }
  //f->ls();
  
  //Raw balance function
  TH1D *gHistBalanceFunction = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunction"));
  gHistBalanceFunction->SetStats(kFALSE);
  gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
  gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
  gHistBalanceFunction->GetZaxis()->SetNdivisions(10);
  gHistBalanceFunction->GetXaxis()->SetTitleOffset(1.3);
  gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
  gHistBalanceFunction->GetZaxis()->SetTitleOffset(1.3);
  gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #eta");
  gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
  gHistBalanceFunction->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");

  //Shuffled balance function
  if(kShowShuffled) {
    TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
    gHistBalanceFunctionShuffled->SetStats(kFALSE);
    gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
    gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
    gHistBalanceFunctionShuffled->GetZaxis()->SetNdivisions(10);
    gHistBalanceFunctionShuffled->GetXaxis()->SetTitleOffset(1.3);
    gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
    gHistBalanceFunctionShuffled->GetZaxis()->SetTitleOffset(1.3);
    gHistBalanceFunctionShuffled->GetXaxis()->SetTitle("#Delta #eta");
    gHistBalanceFunctionShuffled->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
    gHistBalanceFunctionShuffled->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
  }

  //Mixed balance function
  if(kShowMixed) {
    TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
    gHistBalanceFunctionMixed->SetStats(kFALSE);
    gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
    gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
    gHistBalanceFunctionMixed->GetZaxis()->SetNdivisions(10);
    gHistBalanceFunctionMixed->GetXaxis()->SetTitleOffset(1.3);
    gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
    gHistBalanceFunctionMixed->GetZaxis()->SetTitleOffset(1.3);
    gHistBalanceFunctionMixed->GetXaxis()->SetTitle("#Delta #eta");
    gHistBalanceFunctionMixed->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
    gHistBalanceFunctionMixed->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
  }

  //Subtracted balance function
  if(kShowMixed) {
    TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
    gHistBalanceFunctionSubtracted->SetStats(kFALSE);
    gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
    gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
    gHistBalanceFunctionSubtracted->GetZaxis()->SetNdivisions(10);
    gHistBalanceFunctionSubtracted->GetXaxis()->SetTitleOffset(1.3);
    gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
    gHistBalanceFunctionSubtracted->GetZaxis()->SetTitleOffset(1.3);
    gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #eta");
    gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
    gHistBalanceFunctionSubtracted->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
  }

  TString pngName;
  
  TString centralityLatex = "Centrality: ";
  centralityLatex += centralityArray[gCentrality-1]; 
  centralityLatex += "%";

  TString psiLatex;
  if((psiMin == -0.5)&&(psiMax == 0.5))
    psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}"; 
  else if((psiMin == 0.5)&&(psiMax == 1.5))
    psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}"; 
  else if((psiMin == 1.5)&&(psiMax == 2.5))
    psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}"; 
  else 
    psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}"; 
 
  TString pttLatex = Form("%.1f",ptTriggerMin);
  pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
  pttLatex += " GeV/c";

  TString ptaLatex = Form("%.1f",ptAssociatedMin);
  ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
  ptaLatex += " GeV/c";

  TLatex *latexInfo1 = new TLatex();
  latexInfo1->SetNDC();
  latexInfo1->SetTextSize(0.045);
  latexInfo1->SetTextColor(1);

  //Draw the results
  TCanvas *c1 = new TCanvas("c1","Raw balance function 2D",0,0,600,500);
  c1->SetFillColor(10); c1->SetHighLightColor(10);
  c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
  gHistBalanceFunction->SetTitle("Raw balance function");
  gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
  gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
  gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
  gHistBalanceFunction->DrawCopy("lego2");
  gPad->SetTheta(30); // default is 30
  gPad->SetPhi(-60); // default is 30
  gPad->Update();  

  latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
  latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
  latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
  latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());

  if(kShowShuffled) {
    TCanvas *c2 = new TCanvas("c2","Shuffled balance function 2D",100,100,600,500);
    c2->SetFillColor(10); c2->SetHighLightColor(10);
    c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
    gHistBalanceFunctionShuffled->SetTitle("Shuffled events");
    gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.4);
    gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
    gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
    gHistBalanceFunctionShuffled->DrawCopy("lego2");
    gPad->SetTheta(30); // default is 30
    gPad->SetPhi(-60); // default is 30
    gPad->Update();  

    latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
    latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
    latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
    latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
  }

  if(kShowMixed) {
    TCanvas *c3 = new TCanvas("c3","Mixed balance function 2D",200,200,600,500);
    c3->SetFillColor(10); c3->SetHighLightColor(10);
    c3->SetLeftMargin(0.17); c3->SetTopMargin(0.05);
    gHistBalanceFunctionMixed->SetTitle("Mixed events");
    gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.4);
    gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
    gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
    gHistBalanceFunctionMixed->DrawCopy("lego2");
    gPad->SetTheta(30); // default is 30
    gPad->SetPhi(-60); // default is 30
    gPad->Update();  

    latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
    latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
    latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
    latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
  }

  if(kShowMixed) {
    TCanvas *c4 = new TCanvas("c4","Subtracted balance function 2D",300,300,600,500);
    c4->SetFillColor(10); c4->SetHighLightColor(10);
    c4->SetLeftMargin(0.17); c4->SetTopMargin(0.05);
    gHistBalanceFunctionSubtracted->SetTitle("Subtracted balance function");
    gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
    gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
    gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
    gHistBalanceFunctionSubtracted->DrawCopy("lego2");
    gPad->SetTheta(30); // default is 30
    gPad->SetPhi(-60); // default is 30
    gPad->Update();  

    latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
    latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
    latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
    latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
  }
}