Merge branch 'feature-movesplit'

[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / FourierDecomposition / lrc / FourierPlus.C

#include "TROOT.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TNtuple.h"
#include "TFile.h"
#include "TKey.h"
#include "TRandom.h"
#include "TMath.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TBox.h"
#include "TGraphErrors.h"
#include "TStyle.h"
#include "TPolyLine.h"
#include "TLegend.h"
#include "TLine.h"
#include "TSystem.h"
#include "TObjArray.h"
#include "TGraph2D.h"
#include "TGaxis.h"
#include "TArrow.h"
#include "TASImage.h"
#include "TFitResult.h"
#include "TStopwatch.h"
#include "TProfile2D.h"
#include "TRegexp.h"
#include <iostream>
#include <fstream>
#include <vector>

// Plot utilities class - loaded by rootlogon.C + .rootrc
#include "Utils.h"
PlotUtils utils;
const char* inputFileName = ""; // Assign in macro
const char* cfDef = "cA"; // s, m, cA or cB. Only for phi projections.
double minRidgeDEta = 0.8;
double maxRidgeDEta = 1.799;
const int gNMax = 12; // Global const -- largest n in VnDelta

// Don't do global fit to points outside these pt bin
// indices. Globally scoped for availability in vntvna(). Reassigned
// within individual functions.
double minFitPtBin = 0, maxFitPtBin = 999;

// Switches
//============================================================
bool usePtBinCenter = 1;
bool ispp = 0;

// Functions
//============================================================
void SetProtonProtonFlag(int val) { ispp = val; }
void SetMinRidgeDEta(int val) { minRidgeDEta = val; }
bool IsProtonProton() { return ispp; }
void SetCFDef(char* def) { cfDef = def; }
void SaveGraphs(const char* outputFileName, TString opt = "");
TF1* DoGlobalFit(int n, int k, int includedPtBin1 = 0, int includedPtBin2 = 999,
		 const char* region = "RIDGE", const char* cfDef = "cA", TString opt = "");
void VnDelta(int n, const TH1 &h, double &vnd, double &vnd_err, TString opt="");
void VnDelta(int n, const TF1 &f, double &vnd, double &vnd_err, TString opt="");
void MakeVnDVsQGraphs(int n1, int n2, int k, const char* region, const char* corrtype, TString opt="");
TF1* Harmonic(TH1* h, int n, TString opt);
TF1* HarmonicSum(TH1* h, int n1, int n2, TString opt="");
TH1* Hist(const char* region, const char* type, int i, int j, int k, TString opt="");
TH2* EtaPhiHist(int cent, int ptt, int pta);
TH2* PttPtaHist();
TH2F* Agreement2DHist(int k, int n);
void ijkFromHistName(TH1* h, int& i, int& j, int& k);

TGraphErrors* VnDeltaVsN(int i, int j, int k, int nMax, TString opt);
TGraphErrors* VnDeltaNFVsN(int i, int j, int k, int nMax, TString opt);
TGraphErrors* GlobalvnVsN(int i, int k, int nMax, TString opt);
TGraphErrors* VnDeltaVsPtAssoc(int i, int k, int n, TString opt="");
TGraphErrors* VnDeltaVsPtTrig(int j, int k, int n, TString opt="");
TGraphErrors* AgreementVsPtSum(int k, int n);
TGraphErrors* VnDVsQ(int n, int k, const char* region = "RIDGE", const char* cfDef = "cA", TString opt="");
TGraphErrors* GlobalFitVsPtAssoc(int i, int k, int n, TString opt=""); // v_n(ptt)v_n(pta) vs. assoc pt
TGraphErrors* vnGFvsPt(int n, int k, int ipt1, int ipt2,
		       const char* region = "RIDGE", 
		       const char* corrtype = "cA", TString opt = "");
TF1* GlobalFitVsQ(int n, int k, int ipt1=0, int ipt2=999, const char* region="RIDGE",
		  const char* corrtype="cA", TString opt=""); // The real fit function
//TF1* GlobalFitVsQSubset(int n, int k, int i, TString opt = "");
TGraph* Luzumv1(int cent, const char* pid);
double MomCorrection(int i, int j, int k);

// Global fit function - give to TF1 constructor
double vntvna(double *x, double *par);
double Chi2(TH1* h, TF1* f, TString opt = ""); // "", "ndf" or "prob"
double ReducedChi2(int i, int k, int n);
double SineError(int i, int j, int k); // RMS of <sin n D phi> vs. n about zero
double vnGF(int n, int k, int ptBin, int includedPtBin1 = 0, int includedPtBin2 = 999,
	    const char* region = "RIDGE", const char* cfDef = "cA", TString opt = "");

// Show +/- solutions in chi2 space by varying one fixed parameter (par)
TCanvas* FitStudy(int icent, int n, int par);
TCanvas* DrawQ(int k, int n, TString opt = "");
TCanvas* SingleDraw(int cent, int ptt, int pta, TString opt="");
TCanvas* SingleDrawEta(int cent, int ptt, int pta, TString opt);
TCanvas* SingleDrawPlain(int cent, int ptt, int pta, TString opt);
TCanvas* DrawVnVsPtTrig(int k, int npta, int ptabins[], TString opt = "");
TCanvas* SingleDraw2D(int cent, int ptt, int pta, TString opt="");
TCanvas* MultiDraw(int cent, TString opt);
TCanvas* DrawVnFromGlobalFit(int n, int ipt1=0, int ipt2=999, int ncb = 999, 
			     int centbins[] = 0, TString opt=""); // little v_n
TCanvas* TwoPanelVert(double ydiv, const char* canvName, const char* canvTitle);
TCanvas* DrawChi2vsPtTrig(int npta, int ptabins[], int ncb, int centbins[], TString opt);
TCanvas* DrawAgreement(int icent, TString opt = "");
TCanvas* Drawv1to5(int ncb, int centbins[], int ipt1=0, int ipt2=999, TString opt="");
TCanvas* Drawv2to5(int ncb, int centbins[], int ipt1=0, int ipt2=999, TString opt="");
TCanvas* DrawGlobalvnVsN(int npt, int ptbins[], int ncb, int centbins[], TString opt="");

// functions needed for above routine
void initialize(const char* ifName = inputFileName);
void AddPrelimStampToCurrentPad(float x1, float y1, float x2, float y2, TString opt="");

void SaveCanvases(TObjArray* canvases, const char* fileName);
void SaveCanvasesFromFile(const char* rootFile, const char* targetDir, const char* tag, const char* fileType);
TObjArray* GetObjectsFromFile(TFile& file, TString clname, TString dir="");


// Global variables
//============================================================
TFile *fin=0;
int maxcent;
int centlow[100];  // low and high percentile for each selection
int centhigh[100];
int maxpt;
double ptmean[100];
double ptlow[100];
double pthigh[100];

int colorsPale[] = {kOrange-3, kGreen-6, kAzure-9, kGray+1, kRed-9, kRed-7, kBlue-9, kViolet-9};
int colorsDark[] = {kOrange-8, kGreen+1, kAzure+2, kBlack,  kCyan+1, kRed+0, kBlue+1, kViolet+1};
int centColors[] = {kBlack, kRed, kOrange-3, kGreen+1, kCyan+2, kBlue-1, kViolet, kGreen+1, kOrange-3, kRed+0,  kOrange-8, kRed+1, kBlue+1};
int centColorsPale[] = {kGray, kRed-9, kOrange-4, kGreen-7, kCyan-7, kBlue-9, kViolet-9, kGreen, kOrange, kRed-4,  kOrange-9, kRed-1, kBlue-1};

TLatex ltx;  // NDC
TLatex ltx2; // absolute

// Global container for graphs, TF1s, etc.
TObjArray* gList = new TObjArray();
TObjArray* gFitList = new TObjArray();
TF1* axisFn = 0;

TGraphErrors *gi, *gj, *gk;
char* trigLabel(int i);
char* asscLabel(int j);
char* centLabel(int k);
double MeanPt(int i, int j, int k, TString t_or_a = "t" /* "t" or "a" */, 
	      TString opt="" /* "" or "squared" */);
int CentBin(double cntLo, double cntHi); 
int PtBin(double ptLo, double ptHi);
int CentColor(double cen1, double cen2);
int CentColorPale(double cen1, double cen2);
int PtColor(double p1, double p2);
int PtColorPale(double p1, double p2);
TStopwatch timer;

// The business
//============================================================

void SaveGraphs(const char* outputFileName, TString opt)
{
  const char* mode = "recreate";
  if (opt.Contains("update"))
    mode = "update";

  TFile* outputFile = new TFile(outputFileName, mode);
  gList->Write();
  gFitList->Write();
  outputFile->Close();
  return;
}

void initialize(const char* ifName) {

  // TH1F files with histograms for:
  // NSJET_y_PTTRIG_PTASSOC_CENT (tracks with |delta-eta| < 0.8)
  // RIDGE_y_PTTRIG_PTASSOC_CENT (tracks with |delta-eta| = 0.8 - 1.6)

  if (fin) // don't keep re-initializing
    return;

  fin = new TFile(ifName);
  if (!fin) {
    Error("FourierPlus - initialize()", "%s not found", ifName);
    gSystem->Exit(1);
  }
  
  TList* info = (TList*)fin->Get("FileInfo");
  if (info)
    info->Print();

  gi = (TGraphErrors*) fin->Get("TrigPtBins");
  gj = (TGraphErrors*) fin->Get("AsscPtBins");
  gk = (TGraphErrors*) fin->Get("EvCentBins");

  maxpt = gi->GetN();
  maxcent = gk->GetN();
  gList->Add(gi);   gList->Add(gj);   gList->Add(gk);

  if (gj->GetN() != maxpt)
    Warning("initialize()", 
	    "Trigger and assc. binning are not symmetric: %d vs. %d", 
	    maxpt, gj->GetN() );
  
  for (int k=0; k<maxcent; k++) {
    centlow[k] = gk->GetX()[k] - gk->GetEX()[k];
    centhigh[k] = gk->GetX()[k] + gk->GetEX()[k];
    //    cout << centlow[k] << " " << centhigh[k] << endl;
  }

  for (int i=0; i<maxpt; i++) {
    ptlow[i] = gi->GetX()[i] - gi->GetEX()[i];
    pthigh[i] = gi->GetX()[i] + gi->GetEX()[i];
    //    cout << ptlow[i] << " " << pthigh[i] << endl;
  }

  // OK, leave hard-coded for now....
  ptmean[0] = 0.66; // 0.5-1 GeV
  ptmean[1] = 1.30; // 1-2 GeV
  ptmean[2] = 2.30; // 2-3 GeV
  ptmean[3] = 3.30; // 3-4 GeV
  ptmean[4] = 4.60; // 4-6 GeV
  ptmean[5] = 6.60; // 6-8 GeV
  ptmean[6] = 11.7; // 8-15 GeV

  if (usePtBinCenter)
    for (int i=0; i<maxpt; i++) {
      ptmean[i] = gi->GetX()[i];
    }

  gROOT->SetStyle("Plain");
  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  ltx.SetNDC();
}

TH1* Hist(const char* region, const char* type, int i, int j, int k, TString opt)
{
  // Options: 
  // "sys": return measured pts. w/systematic errors. 
  // "hi_sys": meas. pts + sys, unmodified stat. errors.
  // "lo_sys": meas. pts - sys, unmodified stat. errors.

  TH1* h_in=0, *h_out=0;
  const char* name = Form("PbPb/%s_%s_%d_%d_%d",region,type,i,j,k);
  h_in = (TH1*)fin->Get(name);
  if (!h_in) {
    Error("Hist()","%s not found, exiting.", name);
    gSystem->Exit(-1);
  }

  h_out = (TH1*)h_in->Clone(Form("%s_%s", name, opt.Data()));
  
  // Stat --> sys error bars. Error estimated as difference in Fourier
  // series (0<n<6) between same-event dist. and CF.
  if (opt.Contains("sys")) {
    const char* s = Form("PbPb/%s_s_%d_%d_%d",region,i,j,k);
    TH1* htmp = (TH1*)fin->Get(s);
    TH1* hs = 0;
    if (htmp)
      hs = (TH1*)htmp->Clone();
    else
      Error("Hist()", "No hist %s", s);

    double c = double(hs->GetNbinsX())/hs->Integral();
    hs->Scale(c);

    TF1* fs = HarmonicSum(hs,  1,5);
    TF1* fc = HarmonicSum(h_in,1,5);

    for (int n=1; n<=hs->GetNbinsX(); n++) {
      double x = hs->GetBinCenter(n);
      double err = TMath::Abs(fs->Eval(x) - fc->Eval(x));
      double bc = h_in->GetBinContent(n);
      
      if (opt.Contains("hi_sys"))
	h_out->SetBinContent(n, bc+err);
      else if (opt.Contains("lo_sys"))
	h_out->SetBinContent(n, bc-err);
      else {
	h_out->SetBinError(n, err);
	h_out->SetFillColor(kGray);
      }
    }
    
    delete htmp;
    delete hs;
    delete fs;
    delete fc;
  }

  return h_out;
}

void ijkFromHistName(TH1* h, int& i, int& j, int& k)
{
  TString tok, name(h->GetName());
  Ssiz_t from = 0;
  int counter = 0;
  while (name.Tokenize(tok, from, "_")) {
    if (counter==2) i = tok.Atoi();
    if (counter==3) j = tok.Atoi();
    if (counter==4) k = tok.Atoi();
    counter++;
  }
  if (0)
    cout << i<<j<<k << endl;

  return;
}

double PtBinCenterSum(int ptt, int pta)
{
  if ((ptt<0 || ptt>=maxpt) || (pta>ptt))
    Error("PtBinCenterSum()","Bad arg(s): ptt %d pta %d", ptt, pta);
 
  double pti = gi->GetX()[ptt];
  double ptj = gj->GetX()[pta];
  return pti+ptj;
}

int GlobalIndex(int ptt, int pta)
{
  // Find the global index from the trig and assc. pt bins.
  // It is useful that TMath::Binomial(n,k) returns 0 when n < k.
  if (pta > ptt) return -1;
  return TMath::Binomial(ptt+1, 2) + pta;
}

void ijFromGlobalIndex(int q, int& ptt, int& pta)
{
  // Assign ptt and pta from the global index q. This function is
  // effectively the inverse of GlobalIndex(i,j).
  for (int i=0; i<maxpt; i++) {
    for (int j=0; j<=i; j++) {
      if (GlobalIndex(i,j)==q) {
	ptt = i; pta = j;
	return;
      }
    }
  }
  return;
}

// Needs to be updated! Has hard-coded numbers like 27.9. -AMA 11/15/2011
TCanvas* FitStudy(int icent, int n, int par)
{
  // First get the VnDelta vs q graph...
  TGraphErrors* gVn = VnDVsQ(n, icent);
  
  // Then setup fit function
  TF1 *fn = new TF1(Form("fn_cent%dto%d_n%d", centlow[icent],centhigh[icent], n), 
		    vntvna, 0.0, 27.9, maxpt);
  fn->SetNpx(1000);
  fn->SetLineColor(kRed);

  TCanvas* c = new TCanvas(Form("chi2_cent%dto%d_n%d", centlow[icent],centhigh[icent], n), 
			   Form("chi2 test for %s, varying v_{%d}(%.1f)",
				centLabel(icent), n, ptmean[par]), 
			   700, 1200);
  c->Divide(1, 3, 0.001, 0.001);
  
  TGraphErrors *g = new TGraphErrors();
  TGraphErrors *chi2 = new TGraphErrors();
  chi2->SetTitle(Form("Total #chi^{2} for n=%d;v_{%d}(p_{T}=%.2f);total #chi^{2}", n,n, ptmean[par]));
  utils.set_tgraph_props(g, kGray, kGray, kFullCircle, 1.3);
  utils.set_tgraph_props(chi2, kBlue, kBlue, kFullCircle, 0.8);
  g->SetLineWidth(3);

  double bestvn     = vnGF(n,icent,par,0,999,"RIDGE","cA","");
  double bestvn_err = vnGF(n,icent,par,0,999,"RIDGE","cA","err");

  int nsteps = 500;

  double parmax = TMath::Abs(2*(bestvn + bestvn_err));
  double parmin = -parmax;

  // First panel: step thru one fixed vn
  c->cd(1);
  TH1F* hf1 = gPad->DrawFrame(0.0, 1.2*parmin, 12.5, 1.2*parmax);
  hf1->SetTitle(Form("title;p_{T};v_{%d}(p_{T})", n));

  for (int step=0; step<=nsteps; step++) {
    double stepsize =  (parmax-parmin)/nsteps;
    double val = parmin + step*stepsize;
    
    for (int ip=0;ip<maxpt;ip++) fn->SetParameter(ip,0.0);
    fn->FixParameter(par, val);

    gVn->Fit(fn, "QR");
    for (int i=0;i<maxpt;i++) {
      g->SetPoint(i, ptmean[i], fn->GetParameter(i));
      g->SetPointError(i, 0.0, fn->GetParError(i));
    }
    if (step%50==0)
      g->DrawClone("elpsame");

    // store chi2 at this val
    chi2->SetPoint(step, val, fn->GetChisquare());
  }
  // Finally, draw once with no fixed parameters...
  utils.set_tgraph_props(g, kBlack, kBlack, kFullCircle, 1.3);
  for (int ip=0;ip<maxpt;ip++) fn->ReleaseParameter(ip);
  gVn->Fit(fn,"QR");
  for (int i=0;i<maxpt;i++) {
    g->SetPoint(i, ptmean[i], fn->GetParameter(i));
    g->SetPointError(i, 0.0, fn->GetParError(i));
  }
  g->DrawClone("elpsame");
  ltx.SetTextSize(0.08);
  ltx.DrawLatex(0.6, 0.8, Form("v_{%d}, %s", n, centLabel(icent)));

  // Plot chi2 over full range
  c->cd(2);
  c->Update();
  chi2->DrawClone("alp");

  // and zoomed to see 1 sigma interval
  c->cd(3);
  c->Update();
  double ymin = 1e9;//chi2->GetHistogram()->GetMinimum();
  for (int j=0; j<chi2->GetN(); j++)
    if (chi2->GetY()[j] < ymin)
      ymin = chi2->GetY()[j];
  double nsigma = n==2? 3 : 3;
  TH1F* hf3 = gPad->DrawFrame(parmin, ymin - 1, parmax, ymin + nsigma);
  hf3->SetTitle(Form("Total #chi^{2} for n=%d;v_{%d}(p_{T}=%.2f);total #chi^{2}", n,n, ptmean[par]));
  chi2->DrawClone("lpsame");
  TLine best, onesig;
  best.SetLineWidth(3);
  onesig.SetLineWidth(3);
  onesig.SetLineColor(kGray);
  best.DrawLine(parmin, ymin, parmax, ymin);
  onesig.DrawLine(parmin, ymin+1, parmax, ymin+1);
  TLatex lbl;
  lbl.DrawLatex(parmin + 0.05*(parmax-parmin), ymin+0.1, Form("Best v_{%d}(%.2f)",n,ptmean[par]));
  lbl.DrawLatex(parmin + 0.05*(parmax-parmin), ymin+1.1, "+ 1#sigma");

  return c;
}

TCanvas* MultiDraw(int cent, TString opt)
{
  // TODO: make this function much more flexible. Pass in an array of
  //  bins or histos and m,n ints to plot in an m x n grid.
  TCanvas* c = new TCanvas(Form("cfmulti_%s_%d",opt.Data(),cent),
			   Form("Corr. Functions - %s", centLabel(cent)),
			   1200, 900);
  c->Divide(maxpt, maxpt-2, 0, 0);

  int ipad=1;

  // Combine first 3 rows (i=0,1,2) to save space
  // Then draw the rest as usual
  for (int i=0; i<maxpt; i++) {
    for (int j=0; j<maxpt; j++) {
      if(j<=i) {
	TCanvas* cs = 0;
	if (opt.Contains("2D")) 
	  cs = SingleDraw2D(cent, i, j, "");
	else if (opt.Contains("deta"))
 	  cs = SingleDrawEta(cent, i, j, "");
	else
	  cs = SingleDraw(cent, i, j, "");
	c->cd(ipad);
	cs->DrawClonePad();
	ipad++;
	if (ipad==maxpt) ipad++;
      }
      else if (i > 2) ipad++;
    }
  }

  // Add some labeling in empty pads
  ltx.SetTextSize(0.15);
  c->cd(1*maxpt);
  ltx.DrawLatex(0.1, 0.8, "p_{T}^{trig}: 0.5-1,");
  ltx.DrawLatex(0.1, 0.6, "1-2, 2-3 GeV/c");  
  c->cd(2*maxpt);
  ltx.DrawLatex(0.1, 0.8, "p_{T}^{trig}:");
  ltx.DrawLatex(0.1, 0.6, "3-4 GeV/c");  
  c->cd(3*maxpt);
  ltx.DrawLatex(0.1, 0.8, "p_{T}^{trig}:");
  ltx.DrawLatex(0.1, 0.6, "4-6 GeV/c");  
  c->cd(4*maxpt);
  ltx.DrawLatex(0.1, 0.8, "p_{T}^{trig}:");
  ltx.DrawLatex(0.1, 0.6, "6-8 GeV/c");  

  c->cd(2*maxpt-2);
  ltx.SetTextSize(0.2);
  ltx.DrawLatex(0.1, 0.7, "Pb+Pb");
  ltx.DrawLatex(0.1, 0.4, "LHC10h");
  ltx.DrawLatex(0.1, 0.1, centLabel(cent));  

  return c;
}

TH2* EtaPhiHist(int cent, int ptt, int pta)
{
 initialize();
  TH2* htmp = (TH2*)fin->Get(Form("PbPb/ETAPHI_c_%1d_%1d_%1d",ptt,pta,cent));
  if (!htmp) {
    Error("EtaPhiHist()",
	  "Problem getting pointer to histogram. cent %d ptt %d pta %d", 
	  cent, ptt, pta);
    return 0;
  }
  TH2* h = (TH2*)htmp->Clone(Form("EtaPhi_c_%1d_%1d_%1d",ptt,pta,cent));
  TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis(), *az = h->GetZaxis();
  ax->SetTitle("#Delta#phi");
  ay->SetTitle("#Delta#eta");
  az->SetTitle("C(#Delta#phi, #Delta#eta)");
  ay->SetRangeUser(-maxRidgeDEta, maxRidgeDEta);
  ax->SetTitleOffset(1.3);
  ay->SetTitleOffset(1.3);
  az->SetTitleOffset(1.4);

  return h;
}

TCanvas* SingleDraw2D(int cent, int ptt, int pta, TString opt)
{
  initialize();
  TCanvas* c = 0; // The canvas to be returned
  double labelScale = 1.3;
  double nDiv = 5;
  
  // See if we have already made this once
  const char* name;
  if (opt=="") 
    name = Form("etaphi_cent%dto%d_%d_%d",
		centlow[cent],centhigh[cent],ptt,pta );
  else
    name = Form("etaphi_cent%dto%d_%d_%d_%s",
		centlow[cent],centhigh[cent],ptt,pta, opt.Data() );

  const char* title = Form("dphideta_cent%dto%d_pTtrig%.2gto%.2g_pTassoc%.2gto%.2g_%s",
			   centlow[cent],centhigh[cent],
			   ptlow[ptt], pthigh[ptt],
			   ptlow[pta], pthigh[pta],
			   opt.Data());

  if (gDirectory->FindObject(name)) {
    c = (TCanvas*) gDirectory->FindObject(name);
    return c;
  }
  else
    c = new TCanvas(name, title, 1);

  TH2* h = EtaPhiHist(cent, ptt,pta);
  if (!h) {
    Error("SingleDraw2D()",
	  "Problem getting pointer to histogram. cent %d ptt %d pta %d", 
	  cent, ptt, pta);
    return 0;
  }
  utils.make_nice_axes(c, h, labelScale, nDiv, nDiv);
  TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis(), *az = h->GetZaxis();
  ax->SetTitle("#Delta#phi [rad]");
  ay->SetTitle("#Delta#eta");
  az->SetTitle("C(#Delta#phi, #Delta#eta)");
  ay->SetRangeUser(-maxRidgeDEta, maxRidgeDEta);
  ax->SetTitleOffset(1.3);
  ay->SetTitleOffset(1.3);
  az->SetTitleOffset(1.4);

  double czlo=0, czhi=5.0;
  if (opt.Contains("zoom")) 
    az->SetRangeUser(czlo, czhi);

  if (opt.Contains("colz"))
    h->Draw("colz");
  else
    h->Draw("surf1");

  ltx.SetTextSize(0.05);
  if (opt.Contains("colz")) {
    if (opt.Contains("zoom")) 
      ltx.DrawLatex(0.6, 0.9, Form("C(#Delta#phi) (C(#Delta#phi) < %.3g)", czhi));
    else
      ltx.DrawLatex(0.6, 0.9, Form("C(#Delta#phi) (full range)"));
  }
  
  // Label pt, centrality, energy
  TString ptStr(Form("#splitline"
		     "{%.3g < p_{T}^{t} < %.3g GeV/c}"
		     "{%.3g < p_{T}^{a} < %.3g GeV/c}",
		     ptlow[ptt], pthigh[ptt],
		     ptlow[pta], pthigh[pta]));
  ltx.SetTextSize(0.07);
  ltx.DrawLatex(0.01, 0.86, ptStr.Data());
  if (ispp)
    ltx.DrawLatex(0.65, 0.9, Form("#splitline{pp 2.76 TeV}{}"));
  else
    ltx.DrawLatex(0.65, 0.9, Form("#splitline{Pb-Pb 2.76 TeV}{%s}", centLabel(cent)));

  return c;
}

TCanvas* SingleDrawEta(int cent, int ptt, int pta, TString opt)
{
  TCanvas* c = 0;
  int lineWidth  = opt.Contains("small") ? 1 : 4;
  double mrkSize = opt.Contains("small") ? 0.5 : 1.5;

  const char* name = Form("deta_cent%dto%d_%d_%d", 
			  centlow[cent],centhigh[cent],ptt,pta);

  // See if we have already made this once
  if (gDirectory->FindObject(name)) {
    c = (TCanvas*) gDirectory->FindObject(name);
    Info("SingleDrawEta()", "%s already created", c->GetName());
    return c;
  }
  else
    c = new TCanvas(name, name, 600, 600);

  TH1 *hNS = 0, *hAS = 0;    // divide-then-project
  TH1 *hNS2 = 0, *hAS2 = 0;  // project-then-divide

  bool includeProjThenDiv = 0;

  if (1) { // Recreate CFs using divide-then-project method

    TH2* hc = (TH2*)fin->Get(Form("PbPb/ETAPHI_c_%1d_%1d_%1d",ptt,pta,cent));
    if (!hc) {
      Error("SingleDrawEta()", "Problem getting hc\n");
      gSystem->Exit(-1);
    }

    int bin1 = hc->GetXaxis()->FindBin(-TMath::PiOver2()+0.001);
    int bin2 = hc->GetXaxis()->FindBin(TMath::PiOver2()-0.001);
    int bin3 = hc->GetXaxis()->FindBin(TMath::PiOver2()+0.001);
    int bin4 = hc->GetXaxis()->FindBin(3*TMath::PiOver2()-0.001);
    
    bool smallEtaRange = false;
    if (smallEtaRange) {
      bin1 = hc->GetXaxis()->FindBin(-TMath::PiOver2()+0.001);
      bin2 = hc->GetXaxis()->FindBin(TMath::PiOver2()-0.001);
      bin3 = hc->GetXaxis()->FindBin(TMath::Pi()-0.5);
      bin4 = hc->GetXaxis()->FindBin(TMath::Pi()+0.5);
    }

    const char* ns_name = Form("ns_deta_cent%dto%d_%d_%d", 
			       centlow[cent],centhigh[cent],ptt,pta);
    const char* as_name = Form("as_deta_cent%dto%d_%d_%d", 
			       centlow[cent],centhigh[cent],ptt,pta);
    
    hNS = hc->ProjectionY(ns_name, bin1, bin2, "e");
    hAS = hc->ProjectionY(as_name, bin3, bin4, "e");

    hNS->Scale(1./(bin2-bin1+1));
    hAS->Scale(1./(bin4-bin3+1));

  }

  if (includeProjThenDiv) { // As created by Project.C (project-then-divide)
    const char* ns = Form("PbPb/ETA_NS_c_%1d_%1d_%1d",ptt,pta,cent);
    const char* as = Form("PbPb/ETA_AS_c_%1d_%1d_%1d",ptt,pta,cent);
    hNS2 = (TH1*)fin->Get(ns);
    hAS2 = (TH1*)fin->Get(as);
    if (!hNS2) {
      Error("SingleDrawEta()", "Problem getting %s\n", ns);
      gSystem->Exit(-1);
    }
    if (!hAS2) {
      Error("SingleDrawEta()", "Problem getting %s\n", as);
      gSystem->Exit(-1);
    }
  }

  utils.set_hist_props(hNS, kBlack, kNone, kBlack, kFullDotLarge, mrkSize);
  utils.set_hist_props(hAS, kRed, kNone, kRed, kFullDotLarge, mrkSize);
  utils.set_ylimits(hNS, hAS);
  hNS->SetLineWidth(lineWidth);
  hAS->SetLineWidth(lineWidth);

  // hNS->Draw("hist");
  // hAS->Draw("histsame");
  hNS->DrawClone("ep");
  hAS->DrawClone("epsame");

  if (includeProjThenDiv) {
    utils.set_hist_props(hNS2, kGray+2, kNone, kGray+2, kFullDotLarge, mrkSize);
    utils.set_hist_props(hAS2, kCyan, kNone, kCyan, kFullDotLarge, mrkSize);
    hNS2->SetLineWidth(lineWidth);
    hAS2->SetLineWidth(lineWidth);

    hNS2->DrawClone("epsame");
    hAS2->DrawClone("epsame");
  }

  return c;
}

TCanvas* SingleDrawPlain(int cent, int ptt, int pta, TString opt)
{
  initialize();
  const char* region = "RIDGE";
  const char* yTitle = Form("%.2g < |#Delta#eta| < %.2g", 
			    minRidgeDEta, maxRidgeDEta);
  if (opt.Contains("ALL")) {
    region = "ALL";
    yTitle = Form("%.2g < |#Delta#eta| < %.2g", 0.0, maxRidgeDEta);
  }
  if (opt.Contains("NSJET")) {
    region = "NSJET";
    yTitle = Form("%.2g < |#Delta#eta| < %.2g", 0.0, minRidgeDEta);
  }
  
  const char* name = Form("%s_cent%dto%d_%d_%d_%s", region,
			  centlow[cent],centhigh[cent],ptt,pta, opt.Data());
  const char* title = Form("dphiPlain_%s_cent%dto%d_pTtrig%.2gto%.2g_pTassoc%.2gto%.2g_%s",
			   region, 
			   centlow[cent],centhigh[cent],
			   ptlow[ptt], pthigh[ptt],
			   ptlow[pta], pthigh[pta],
			   opt.Data());
  TCanvas* c = 0;
  if (gDirectory->FindObject(name)) {

    c = (TCanvas*) gDirectory->FindObject(name);
    Info("SingleDrawPlain()", "%s already created", c->GetName());
    return c;
  }
  else
    c = new TCanvas(name, title, 700, 700);
  c->cd();

  TH1* h_sys = Hist(region, cfDef, ptt,pta,cent, "sys");
  TH1* h     = Hist(region, cfDef, ptt,pta,cent);
  h->SetName("h");

  TAxis* ax = h->GetXaxis();
  TAxis* ay = h->GetYaxis();

  h->SetTitle(Form(";#Delta#phi [rad];C(#Delta#phi),  %s", yTitle));
  h->SetLineWidth(2);
  ax->SetNdivisions(208);
  ay->SetNdivisions(208);
  ax->SetLabelSize(0.04);
  ay->SetLabelSize(0.04);
  ax->SetTitleSize(0.05);
  ay->SetTitleSize(0.05);
  ax->SetTitleOffset(1.5);
  ay->SetTitleOffset(-1.9);
  c->SetFrameLineWidth(2);
  c->SetLeftMargin(0.2);
  c->SetBottomMargin(0.15);
  c->SetTopMargin(0.01);
  c->SetRightMargin(0.01);
  ax->CenterTitle();
  ay->CenterTitle();
  ax->SetTicks("+-");
  ay->SetTicks("+-");

  utils.set_ylimits(h, h, 0.05);
  utils.set_hist_props(h, kBlack, kNone, kBlack, kFullDotLarge, 1.4);
  h->Draw();
  h_sys->Draw("e2psame");
  h->Draw("same");
  if (opt.Contains("harmonics"))
    for (int n=1; n<6; n++)
      Harmonic(h, n, "draw");
  if (opt.Contains("sum"))
    HarmonicSum(h,1,5,"draw");
  
  ltx.SetTextSize(0.045);
  double top = 0.92, gap = 0.08;
  double x1 = (TString(region)=="RIDGE" && (ptlow[ptt] > 4 || ispp))? 0.25 : 0.5;
  ltx.DrawLatex(x1, top-0*gap, Form("p_{T}^{t} %s GeV/c", trigLabel(ptt)));  
  ltx.DrawLatex(x1, top-1*gap, Form("p_{T}^{a} %s GeV/c", asscLabel(pta)));  
  if (ispp)
  ltx.DrawLatex(x1, top-2*gap, Form("pp 2.76 TeV"));  
  else
  ltx.DrawLatex(x1, top-2*gap, Form("Pb-Pb 2.76 TeV, %s", centLabel(cent)));  
  if (0)
    ltx.DrawLatex(0.22, 0.92, Form("|#Delta#eta| > %.1f", minRidgeDEta));  

  // ltx.SetTextColor(kBlack);
  // ltx.SetTextSize(0.03);
  // ltx.DrawLatex(0.7, 0.18, "statistical error only");  
 
  return c;
}

TF1* Harmonic(TH1* h, int n, TString opt)
{
  TF1* f = 0;
  double vndelta=0, vndelta_err=0;
  double phiMin = -0.5*TMath::Pi();
  double phiMax = +1.5*TMath::Pi();
  double b0 = h->Integral() / h->GetNbinsX();
  int color = opt.Contains("gray")? kGray : colorsPale[n-1];

  VnDelta(n, *h, vndelta, vndelta_err);

  f = new TF1(Form("%scos%ddphi", h->GetName(), n), 
	      "[0]*(1. + 2.*[1]*cos([2]*x))",
	      phiMin, phiMax);
    f->SetParameters(b0, vndelta, n);
    f->SetLineColor(color);
    f->SetLineWidth(3);
    if (opt.Contains("draw")) 
      f->Draw("same");

    return f;
}

TF1* HarmonicSum(TH1* h, int n1, int n2, TString opt)
{
  if (n1<1)
    Warning("HarmonicSum()", "n1=%d but should be > 0", n1);

  double phiMin = -0.5*TMath::Pi();
  double phiMax = +1.5*TMath::Pi();

  // Function definition string
  TString fdef("[0]*(1.0 + 2.0*(");
  for (int n=n1; n<=n2; n++) {
    fdef.Append(Form(" [%d]*cos(%d*x) ", n, n));
    fdef.Append( n<n2 ? "+" : "))");
  }

  // Normalization
  TF1 *ff = new TF1("ff", fdef.Data(), phiMin,phiMax);
  ff->SetParameter(0, h->Integral()/h->GetNbinsX());

  // Set parameters 
  for (int n=n1; n<=n2; n++) {

    // Sum of vnt x vna from global fit
    if (opt.Contains("global") ) {
      int i=-1, j=-1, k=-1;
      ijkFromHistName(h,i,j,k);
      double vnt = vnGF(n,k,i,0,999,"RIDGE","cA","");
      double vna = vnGF(n,k,j,0,999,"RIDGE","cA","");
      ff->SetParameter(n, vnt*vna); // TODO: consider including errors
      ff->SetLineColor(kRed+1);
    }
    else { // sum of directly-extracted pair moments
      double Vn=0, VnErr=0;
      VnDelta(n, *h, Vn, VnErr);
      ff->SetParameter(n, Vn);
      ff->SetLineColor(kBlack); // ff->SetLineColor(kAzure);
      //      ff->SetLineWidth(4);
      ff->SetLineStyle(kDashed);
      if (n2==10) {
	ff->SetLineColor(kGreen);
      }
    } 
  }

  if (opt.Contains("draw")) 
    ff->Draw("same");
  
  return ff;
}

TCanvas* SingleDraw(int cent, int ptt, int pta, TString opt)
{
  int lineWidth  = opt.Contains("small") ? 1 : 4;
  double mrkSize = opt.Contains("small") ? 0.5 : 1.5;
  TString dEtaLabel = "";

  const char* region = "RIDGE";
  if (opt.Contains("ALL")) {
    region = "ALL";
    dEtaLabel = Form("%.2g < |#Delta#eta| < %.2g", 0.0, maxRidgeDEta);
  }
  if (opt.Contains("NSJET")) {
    region = "NSJET";
    dEtaLabel = Form("%.2g < |#Delta#eta| < %.2g", 0.0, minRidgeDEta);
  }
  if (TString(region)=="RIDGE")
    dEtaLabel = Form("%.2g < |#Delta#eta| < %.2g", minRidgeDEta, maxRidgeDEta);
  
  initialize();

  TCanvas* c1 = 0; // The canvas to be returned

  // See if we have already made this once
  TString name = Form("dphi_etamin%02d_cent%dto%d_%d_%d", (int)(10*minRidgeDEta), 
			  centlow[cent],centhigh[cent],ptt,pta);
  TString title = Form("dphiFourier_%s_cent%dto%d_pTtrig%.2gto%.2g_pTassoc%.2gto%.2g",
			   region, 
			   centlow[cent],centhigh[cent],
			   ptlow[ptt], pthigh[ptt],
			   ptlow[pta], pthigh[pta]);

  if (!opt.IsNull()) {
    name.Append("_"); name.Append(opt);
    title.Append("_"); title.Append(opt);
  }

  if (gDirectory->FindObject(name.Data())) {

    c1 = (TCanvas*) gDirectory->FindObject(name.Data());
    Info("SingleDraw()", "%s already created", c1->GetName());
    return c1;
  }
  else
    c1 = TwoPanelVert(0.3, name.Data(), title.Data());

  c1->cd(1);
  gPad->SetTopMargin(0.1);
  c1->cd();

  TH1* h = Hist(region, cfDef, ptt,pta,cent);
  TH1* h_sys = Hist(region, cfDef, ptt,pta,cent, "sys");
  utils.set_hist_props(h, kBlack, kNone, kBlack, kFullDotLarge, mrkSize);
  h->SetLineWidth(opt.Contains("small") ? 1 : 2);
  h->GetXaxis()->SetTitle("#Delta#phi [rad]");
  h->GetYaxis()->SetTitle("C(#Delta#phi)");

  TF1* ff5 = HarmonicSum(h, 1, 5, "");
  ff5->SetLineWidth(lineWidth);
  TF1* fflowg = HarmonicSum(h, 1, 5, "global");
  fflowg->SetLineWidth(lineWidth);
  TF1* ff10 = HarmonicSum(h, 1, 10, "");
  ff10->SetLineWidth(lineWidth);
  TF1* ff12 = HarmonicSum(h, 1, 12, "");
  ff12->SetLineWidth(lineWidth);

  c1->cd(1);

  int nXticks = 404, nYticks = 210;
  double pad1Scale = 1.6;
  double pad2Scale = 2.3;
  utils.make_nice_axes(c1, h, pad1Scale, nXticks, nYticks, 0.05, 0.01);
  h->GetXaxis()->SetTitleOffset(1.2);
  h->GetYaxis()->SetTitleOffset(1.5);
  utils.set_ylimits(h,h);

  h->Draw("ep");
  h_sys->Draw("e2psame");
  h->Draw("epsame");

  if (opt.Contains("gray"))
    for (int n=1; n<6; n++)
      Harmonic(h, n, "gray, draw");
  else if (opt.Contains("color"))
    for (int n=1; n<6; n++)
      Harmonic(h, n, "draw");
  
  if (opt.Contains("upto10")) {
    ff10->SetLineWidth(lineWidth+2);
    ff10->SetLineColor(kBlack);
    ff10->DrawClone("same");
    ff10->SetLineWidth(lineWidth);
    ff10->SetLineColor(kGreen);
    ff10->SetLineStyle(kDashed);
    ff10->Draw("same");
  }

  //  ff10->Draw("same");
  ff5->Draw("same");
  h->Draw("p,e,same");
  
  if (opt.Contains("global")) {
    fflowg->Draw("same");
  }

  ltx.SetTextSize(0.06);
  double top = 0.82, gap = 0.08;
  double x1 = (TString(region)=="RIDGE" && ptlow[ptt] > 6)? 0.25 : 0.6;
  if (ispp)
    ltx.DrawLatex(0.2, 0.92, Form("pp 2.76 TeV"));  
  else
    ltx.DrawLatex(0.2, 0.92, Form("Pb-Pb 2.76 TeV, %s central", centLabel(cent)));  
  ltx.DrawLatex(x1, top-0*gap, Form("%.3g < p_{T}^{t} < %.3g GeV/c",
				    ptlow[ptt], pthigh[ptt]));  
  ltx.DrawLatex(x1, top-1*gap, Form("%.3g < p_{T}^{a} < %.3g GeV/c", 
				    ptlow[pta], pthigh[pta]));  
  // ltx.DrawLatex(x1, top-0*gap, Form("p_{T}^{trig} %s GeV/c", trigLabel(ptt)));  
  // ltx.DrawLatex(x1, top-1*gap, Form("p_{T}^{assoc} %s GeV/c", asscLabel(pta)));  
  ltx.DrawLatex(x1, top-2*gap, dEtaLabel.Data());  


  c1->cd(2);

  TH1F *hsub = static_cast <TH1F *> (h->Clone());
  TH1F *hsub10 = static_cast <TH1F *> (h->Clone());
  TH1F *hsubg = (TH1F*) (h->Clone());

  utils.set_hist_props(hsub10,  kBlack, kNone, kBlack, kOpenCircle, mrkSize);
  utils.set_hist_props(hsub,  kAzure, kNone, kAzure,   kFullSquare, mrkSize);
  utils.set_hist_props(hsubg, kRed+1, kNone, kRed+1, kFullSquare, mrkSize);

  for (int i=1;i<=hsub->GetNbinsX();i++) {
    double val, err;
    val = h->GetBinContent(i) * h->GetBinWidth(i);
    err = h->GetBinError(i)   * h->GetBinWidth(i);
    double dphi1 = h->GetXaxis()->GetBinLowEdge(i);
    double dphi2 = h->GetXaxis()->GetBinUpEdge(i);
    double fSum5, fSum10, fSumGlobal;
    fSum5      = ff5->Integral(dphi1, dphi2);
    fSum10     = ff10->Integral(dphi1, dphi2);
    fSumGlobal = fflowg->Integral(dphi1, dphi2);
    hsub->SetBinContent(i,   val/fSum5);
    hsub->SetBinError(i,     err/fSum5);
    hsub10->SetBinContent(i, val/fSum10);
    hsub10->SetBinError(i,   err/fSum10);
    hsubg->SetBinContent(i,  val/fSumGlobal);
    hsubg->SetBinError(i,    err/fSumGlobal);
  }
  utils.make_nice_axes(c1, hsub, pad2Scale, nXticks, 5, 0.05, 0.01);
  if (opt.Contains("global"))
    utils.set_ylimits(hsub, hsubg, 0.2);
  else 
    utils.set_ylimits(hsub, hsub10, 0.2);
  
  hsub->SetTitle(Form(";#Delta#phi [rad];ratio"));
  //  hsub->SetTitle(Form(";#Delta#phi [rad];data/#Sigma_{n=1}^{%s}", opt.Contains("global")?"5":"N"));
  hsub->GetXaxis()->SetTitleOffset(1.25);
  hsub->GetYaxis()->SetTitleOffset(0.7);
  hsub->GetXaxis()->SetTicks("+-");   hsub->GetXaxis()->SetTickLength(0.08);
  hsub10->SetTitle(";#Delta#phi [rad];data/#Sigma_{n=1}^{N}");
  hsub10->GetXaxis()->SetTitleOffset(1.25);
  hsub10->GetYaxis()->SetTitleOffset(0.7);
  hsubg->SetTitle(";#Delta#phi [rad];data/calc");
  hsubg->GetXaxis()->SetTitleOffset(1.25);
  hsubg->GetYaxis()->SetTitleOffset(0.7);

  TLegend* leg = new TLegend(0.2, 0.01, 0.38, 0.38);
  leg->SetFillColor(kNone);
  leg->SetBorderSize(1);
  if (opt.Contains("global")) {
    leg->AddEntry(hsub, "#LTcos(n#Delta#phi)#GT", "epl");
    leg->SetBorderSize(0);
  }
  else
    leg->AddEntry(hsub, "1#leqn#leq5", "epl");

  TLine unity;
  hsub->Draw("ep");

  if (opt.Contains("global")) {
    hsubg->Draw("epsame");
    leg->AddEntry(hsubg, "v(p_{Tt})v(p_{Ta})", "epl");
  }
  else if (opt.Contains("upto10")) {
    TH1* hsub10green = (TH1*)hsub10->Clone();
    utils.set_hist_props(hsub10green,  kBlack, kNone, kGreen, kFullCircle, mrkSize);
    hsub10green->Draw("epsame");
    hsub10->Draw("epsame");
    leg->AddEntry(hsub10green, "1#leqn#leq10", "epl");
  }

  unity.DrawLine(-TMath::PiOver2(), 1.0, 3*TMath::PiOver2(), 1.0);

  c1->cd(2);
  if (opt.Contains("global"))
    leg->Draw();

  if (1) {
    // calculate chi^2
    TF1* flat = new TF1("flat", "1.0", -TMath::PiOver2(), 3*TMath::PiOver2());
    double chi2 = Chi2(hsub, flat);
    double chi2_10 = Chi2(hsub10, flat);
    int ndof = hsub->GetNbinsX();
    delete flat;
    c1->cd(1);
    ltx.DrawLatex(0.65, 0.05, Form("#chi^{2}/ndf = %.3g / %d", chi2, ndof-1));
    if (0)
      cout << "Chi2[5] = " << chi2 << " Chi2[10] = " << chi2_10 
	   << " ndof = " << ndof-1 << " (Prob = " << TMath::Prob(chi2,ndof-6) 
	   << " , " << TMath::Prob(chi2_10,ndof-11) << " )" << endl;
  }
  return c1;
}

double Chi2(TH1* h, TF1* f, TString opt) {
  double chi2 = 0;
  double ndf = 0;
  for (int i=1;i<=h->GetNbinsX();i++) {
    double fy = f->Eval(h->GetBinCenter(i));
    chi2 += pow( (h->GetBinContent(i) - fy)/h->GetBinError(i) , 2.0);
    ndf++;
  }
  ndf -= 1;

  if (opt.Contains("ndf") && ndf > 0)
    chi2 /= ndf;

  if (opt.Contains("prob") )
    return TMath::Prob(chi2, ndf);
 
  return chi2; 
}

double vntvna(double *x, double *par) 
{
  // This function is called many times during a fit. For a 1-d fit,
  // x[0] is the x-position of the nth point,
  // i.e. graph->GetX()[n]. For a single-valued TGraph, it can be used
  // to indicate which point the fitter is currently working on.  The
  // par array contains maxpt free parameters that are continually
  // updated during the fit process.

  // q is the global index
  int q = (int)x[0]; // Assume points are at q + 0.5.
  int i = 999; 
  int j = 999;
  int nPoints = (maxpt*(maxpt+1))/2;
  TF1::fgRejectPoint = kFALSE;

  // Assign i and j from q
  ijFromGlobalIndex(q, i, j);

  if (q < 0 || q >= nPoints) {
    Error("vntvna()", 
	  "Problem finding global index. x[0]=%f, q=%d, i=%d, j=%d",
	  x[0], q, i, j);
    
    return 999;
  }
  
  // Exclude points from fit. Warning: TF1::fgRejectPoint is a global
  // flag checked by all TF1s, beware of side-effects.
  if (j < minFitPtBin || j > maxFitPtBin) {
    TF1::RejectPoint();
    if (0)
      cout << Form("Rejecting x[0],q,i,j = %f %d %d %d", 
		   x[0], q, i, j)
	   << endl;
    return 0;
  }
  
  double v2t = par[i];
  double v2a = par[j];
  int n = par[maxpt];                       // Fourier coeff. index
  int k = par[maxpt+1];                     // cent. bin index
  double fitval = v2t * v2a;

  // Add extra fit term for momentum conservation
  // It is 1/<\sum pt2> * <ptt><pta>
  // From ML & JYO, PRL 106, 102301 (2011) eqs 3-7
  double mptt=0, mpta=0, v1fac=0;
  if (n==1) {
    mptt = MeanPt(i, j, k, "t");
    mpta = MeanPt(i, j, k, "a");
    v1fac = par[maxpt+2];

    fitval += v1fac * mptt * mpta;  
  }
  
  // Odd vnd's < 0 at high pt. // can I replace the below by saying if
  // (vndelta < 0) then make fitval < 0?
  if (j >= PtBin(5., 6.) && (n%2) ) 
    fitval = -fitval;

  return fitval;
}

void VnDelta(int n, const TH1 &h, double &vnd, double &vnd_err, TString opt)
{
  // Extract fourier coefficient n directly from h. If opt is "hi" or
  // "lo", the coefficient is calculated for the points + or - their
  // errors. "hi" and "lo" are meant for an h with systematic (not
  // statistical) errors, so vnd_err is not calculated.
  double VN = 0.0;
  double norm = 0.0;

  // Calculate and set coeff. vnd
  for (int ib=1; ib<=h.GetNbinsX(); ib++) {
    double deltaphi = h.GetBinCenter(ib);
    double weight   = h.GetBinContent(ib);

    if (opt.Contains("hi"))
      weight += h.GetBinError(ib);
    else if (opt.Contains("lo"))
      weight -= h.GetBinError(ib);

    if (opt.Contains("sin") )
      VN += sin(n * deltaphi ) * weight;
    else
      VN += cos(n * deltaphi ) * weight;
    
    norm += weight;
  }
  
  if (norm==0) 
    ::Error("VnDelta()", "Div/0 error");

  VN = VN / norm;

  // // Mom. cons correction?????????????????????????????????
  // if (n==1) VN += 0.002*ptmean[i]*ptmean[j];

  vnd = VN;
  if (opt.Contains("hi") || opt.Contains("lo")) {
    vnd_err = 0;
    return;
  }

  // Statistical uncertainty vnd_err
  double quad_sum_uncertainty = 0.0;
  for (int ib=1; ib<=h.GetNbinsX(); ib++) {
    double deltaphi = h.GetBinCenter(ib);
    double dfdwi = cos(n * deltaphi)/norm - VN/norm;

    if (opt.Contains("sin") )
      dfdwi = sin(n * deltaphi)/norm - VN/norm;

    double sigma = h.GetBinError(ib);
    quad_sum_uncertainty += dfdwi*dfdwi * sigma*sigma;
  }

  vnd_err = TMath::Sqrt(quad_sum_uncertainty);
  return;
}

void VnDelta(int n, const TF1 &f, double &vnd, double &vnd_err, TString opt)
{
  // Extract fourier coefficient n from f.
  double VN = 0.0;
  double norm = 0.0;
  int nSteps = 100;
  double x1=0, x2=0;
  f.GetRange(x1, x2);
  double dx = (x2-x1)/nSteps;

  if (!opt.IsNull() ) cout<<opt.Data()<<endl;

  // Calculate and set coeff. vnd
  for (int ib=0; ib<nSteps; ib++) {

    double deltaphi = x1 + ib*dx; 
    double weight   = f.Eval(deltaphi);

    VN += cos(n * deltaphi ) * weight;
    norm += weight;
  }
  
  if (norm==0) 
    Error("VnDelta()", "Div/0 error");

  VN = VN / norm;
  vnd = VN;
  vnd_err = 0; // for now
  return;
}



// The argument ydiv is the division between pads as a fraction of the
// height.
TCanvas* TwoPanelVert(double ydiv, const char* canvName, const char* canvTitle)
{
  TCanvas* c = new TCanvas(canvName, canvTitle, 800, 800);
  // Extents of pads inside canvas
  double x1=0.01, x2=0.99, y1=0.01, y2=0.99;
  double lm=0.2, rm=0.02;

  // Division between top and bottom pads
  double ysplit = y1 + ydiv*(y2-y1);

  // Divide first, adjust afterward.
  c->Divide(1,2);
  TPad* ptop = (TPad*)(c->GetPad(1));
  TPad* pbot = (TPad*)(c->GetPad(2));

  ptop->SetPad(x1, ysplit, x2, y2);
  ptop->SetMargin(lm, rm, 0.01, 0.01); // left, right, bottom, top.
  ptop->SetFrameLineWidth(2);

  pbot->SetPad(x1, y1, x2, ysplit);
  pbot->SetMargin(lm, rm, 0.4, 0.01); // left, right, bottom, top.
  pbot->SetFrameLineWidth(2);
  pbot->SetTicks();

  c->cd();
  return c;
}

char* trigLabel(int i)
{
  if (i<0 || i>maxpt)
    return Form("Error: no trig pt bin %d", i);
  
  double x1 = gi->GetX()[i] - gi->GetEX()[i];
  double x2 = gi->GetX()[i] + gi->GetEX()[i];

  // int prec1 = 0, prec2 = 0;
  // if (x1 - (int)x1 > 0)
  //   prec1 = 1;
  // if (x2 - (int)x2 > 0)
  //   prec2 = 1;

  return Form("%.3g-%.3g", x1, x2);
}

char* asscLabel(int j)
{
  if (j<0 || j>maxpt)
    return Form("Error: no assc pt bin %d", j);
  
  double x1 = gj->GetX()[j] - gj->GetEX()[j];
  double x2 = gj->GetX()[j] + gj->GetEX()[j];

  /* Stupid coding...
  int prec1 = 0, prec2 = 0;

  if (x1 - (int)x1 > 0)
    prec1 = 1;
  if (x2 - (int)x2 > 0)
    prec2 = 1;
  */

  return Form("%.3g-%.3g", x1, x2);
}

char* centLabel(int k)
{
  if (k<0 || k>maxcent)
    return Form("Error: no cent bin %d", k);
  
  double c1 = gk->GetX()[k] - gk->GetEX()[k];
  double c2 = gk->GetX()[k] + gk->GetEX()[k];

  return Form("%.0f-%.0f%%", c1, c2);
}


void VnLimits(double &minVn_k, double &maxVn_k, int k, int imin, int imax)
{
  // Find limits for each centrality bin k between imin (inclusive)
  // and imax (also inclusive) to plot a common y axis.
  //  for (int k=0; k < maxcent; k++) {
    for (int n=1; n<=5; n++) {
      for (int i=imin; i<=imax; i++) {
	TGraphErrors* gc = VnDeltaVsPtAssoc(i, k, n);
	if (!gc) {
	  Error("VnLimits()", "No graph i %d k %d n %d", i,k,n );
	  continue;
	}
	for (int j=0; j<gc->GetN(); j++) {
	  double vn = gc->GetY()[j];
	  double evn = gc->GetEY()[j];
	  if (vn+evn > maxVn_k) maxVn_k = vn+evn;
	  if (vn-evn < minVn_k) minVn_k = vn-evn;
	}
      }
    }

  return;
}

TCanvas* DrawQ(int k, int n, TString opt)
{
  TStopwatch ts1, ts2;
  int fitCurveColor = kRed-4;
  int divColor = kBlue;
  const char* region = "RIDGE";
  if (opt.Contains("ALL"))
    region = "ALL";

  TString title = Form("globalfit_%d_cent%dto%d",
			   n, centlow[k],centhigh[k]);
  TString name = Form("V%dvsQ_etamin%02d_cent%dto%d",
		      n, (int)(10*minRidgeDEta), 
		      centlow[k],centhigh[k]);
  
  if (!opt.IsNull()) {
    name.Append("_");
    name.Append(opt.Data());
    title.Append("_");
    title.Append(opt.Data());
  }

  TCanvas* c = TwoPanelVert(0.4, name.Data(), title.Data());
  int cwidth = opt.Contains("highptfit") ? 500 : 1200; 
  c->SetCanvasSize(cwidth, 450);
  c->SetWindowSize(cwidth+50, 500);

  if (opt.Contains("highptfit")) {
    c->GetPad(1)->SetLeftMargin(0.2);
    c->GetPad(2)->SetLeftMargin(0.2);
    c->GetPad(1)->SetRightMargin(0.05);
    c->GetPad(2)->SetRightMargin(0.05);

  }

  TLegend* lq = new TLegend(0.14, 0.56, 0.25, 0.83);
  lq->SetFillColor(kNone);
  lq->SetBorderSize(0);

  // g:  global VnDelta vs global index (stat errors).
  // gs: global VnDelta vs global index (sys errors).
  // fn: global fit TF1.
  // r:  ratio of g/fit.
  TString opt1 = "";
  if (opt.Contains("ALL"))
    opt1.Append("_ALL");
  if (opt.Contains("ptcons"))
    opt1.Append("_ptcons");

  TGraphErrors* g  = VnDVsQ(n, k, region, cfDef, opt1);
  opt1.Append("_sys");
  TGraphErrors* gs = VnDVsQ(n, k, region, cfDef, opt1);

  int ptaLo = 0, ptaHi = 999;
  if (opt.Contains("lowptfit")) {
    ptaHi = PtBin(3., 4.);
  }
  if (opt.Contains("highptfit")) {
    ptaLo = PtBin(5., 6.);
  }
  // If custom pta range was passed in, e.g.
  // Form("ptabin%dto%d", PtBin(1.5, 2.0), PtBin(8, 15))
  if (opt.Contains("ptabin")) {
    TRegexp re("ptabin[0-9]+to[0-9]+");
    TRegexp re1("ptabin[0-9]+");
    TRegexp re2("to[0-9]+");
    TString rx = opt(re), pta1s = rx(re1), pta2s = rx(re2);
    pta1s.ReplaceAll("ptabin","");
    pta2s.ReplaceAll("to","");
    ptaLo = pta1s.Atoi();
    ptaHi = pta2s.Atoi();
  }

  TF1* fn    = GlobalFitVsQ(n, k, ptaLo, ptaHi, region, cfDef, "");
  TF1* fn_hi = GlobalFitVsQ(n, k, ptaLo, ptaHi, region, cfDef, "hi_sys");
  TF1* fn_lo = GlobalFitVsQ(n, k, ptaLo, ptaHi, region, cfDef, "lo_sys");

  // Systematic band on ratio line at unity
  TH1F *r_sys = new TH1F(Form("r_sys_%s", title.Data()), "r_sys", 
			 g->GetN()-0.5, 0, g->GetN()-0.5);
  utils.set_hist_props(r_sys, kBlack, kGray, kBlack, kDot, 1.0);

  TGraphErrors* r = new TGraphErrors();
  TString rname(g->GetName());
  rname.Append("_ratio");
  r->SetName(rname.Data());

  if (!fn) {
    Error("DrawQ()","Problem getting/creating global fit function");
    gSystem->Exit(-1);
  }

  // Set visual properties
  utils.set_tgraph_props(g, kBlack, kBlack, kFullSquare, 1.0);
  utils.set_tgraph_props(gs, kBlack, kBlack, kDot, 0.5);
  gs->SetFillColor(kGray);
  utils.set_tgraph_props(r, kBlack, kBlack, kFullCircle, 1.0);
  g->SetLineWidth(1);
  r->SetLineWidth(1);
  fn->SetLineWidth(1);
  fn->SetLineColor(fitCurveColor);
  fn_hi->SetLineStyle(kDashed);
  fn_lo->SetLineStyle(kDashed);
  fn_hi->SetLineWidth(1);
  fn_lo->SetLineWidth(1);

  // Compute ratio for r graph
  for (int q=0; q<g->GetN(); q++) {
    int ii=0, jj=0;
    ijFromGlobalIndex(q, ii, jj);

    double fy = fn->Eval(q+0.5);
    double y = fy==0? 1e9 : g->GetY()[q]/fy;
    double rerr = fn_hi->Eval(q+0.5) / fy - 1;

    if (opt.Contains("highptfit") && n%2 && g->GetY()[q]!=0)
      y = fy/g->GetY()[q];

    // No error band when no fit point!
    if (jj < ptaLo || jj > ptaHi)
      rerr = 0;

    r->SetPoint(q, g->GetX()[q], y);
    r->SetPointError(q, 0.0, TMath::Abs(g->GetEY()[q]/fy));

    r_sys->SetBinContent(q+1, 1.0);
    r_sys->SetBinError(q+1, rerr);
  }

  // Set up plot ranges for top panel (g and f)
  int nx = g->GetN();
  double x1,y1,x2,y2;
  double sf = 1; // scale factor
  fn->GetRange(x1,x2);
  if (opt.Contains("highptfit")) {
    x1 = 45;
  }
  y1 = sf*fn->GetMinimum(x1,x2);
  y2 = sf*fn->GetMaximum(x1,x2);
  double marg = 0.4*(y2-y1);
  //  if (opt.Contains("lowptfit") || opt.Contains("highptfit") ) marg = 2*(y2-y1);
  x1 -= 0.5;
  x2 += 0.5;
  y1 -= marg;
  y2 += marg;

  // Set up plot ranges for lower panel (r)
  double ry1=0.21, ry2=1.79;
  if (opt.Contains("highptfit")) {
    ry1 = 0.01;
    ry2 = 1.99;
  }

  // Get / make frame histos for top and bottom
  TH2F* hf1 = 0; // (TH2F*)gDirectory->Get("qhf1");
  TH2F* hf2 = 0; // (TH2F*)gDirectory->Get("qhf2");
  TString name1 = Form("qhf1_%s", title.Data());
  TString name2 = Form("qhf2_%s", title.Data());
  // TString sv = n==1 ? "#LTcos(#Delta#phi)#GT" :
  // Form("#LTcos(%d#Delta#phi)#GT", n);
  TString sv = Form("V_{%d#Delta}", n);
  TString sgf = Form("(v_{%d}^{t}v_{%d}^{a})_{fit}", n, n);

  if (!hf1)
    hf1 = new TH2F(name1.Data(), "hf1", nx+1, x1, x2, 1000, y1, y2);
  if (!hf2)
    hf2 = new TH2F(name2.Data(), "hf2", nx+1, x1, x2, 1000, ry1, ry2);

  // Axes
  TAxis* ax1 = hf1->GetXaxis();
  TAxis* ay1 = hf1->GetYaxis();
  TAxis* ax2 = hf2->GetXaxis();
  TAxis* ay2 = hf2->GetYaxis();
  hf1->SetBins(nx+1, x1, x2, 100, y1, y2);
  hf2->SetBins(nx+1, x1, x2, 100, ry1, ry2);
  hf1->SetAxisRange(x1, x2, "X");  
  hf1->SetAxisRange(y1, y2, "Y");  
  hf2->SetAxisRange(x1, x2, "X");  
  hf2->SetAxisRange(ry1, ry2, "Y");  
  ax1->SetLabelOffset(0.02);
  ax1->SetLabelSize(0.04);  ay1->SetLabelSize(0.08);
  ax2->SetLabelSize(0.18);   ay2->SetLabelSize(0.12);
  ax1->SetNdivisions(1);   ax2->SetNdivisions(1);
  ay1->SetNdivisions(208);  ay2->SetNdivisions(opt.Contains("highptfit")?104:210);
  ax1->SetTickLength(0.02);

  //ay1->SetTitle(Form("%s and #color[%d]{%s}",sv.Data(), fitCurveColor, sgf.Data()));
  ay1->SetTitle(Form("%s", sv.Data()));
  ay2->SetTitle(Form("#frac{%s}{fit}", sv.Data()));
  //  ay2->SetTitle(Form("#frac{%s}{%s}", sv.Data(), sgf.Data()));
  // ax2->SetTitle(Form("#splitline{associated p_{T}^{a}}{#color[%d]{   trigger p_{T}^{t}}} [GeV/c]", divColor));
  ax2->SetTitle(Form("#color[%d]{p_{T}^{t}}, p_{T}^{a} [GeV/c]", divColor));
  ay1->SetTitleOffset(opt.Contains("highptfit") ? 1. : 0.45);
  ay1->SetTitleSize(0.09);
  ay1->CenterTitle();
  ay2->SetTitleOffset(opt.Contains("highptfit") ? 0.6 : 0.26);
  ay2->SetTitleSize(0.14);
  ay2->CenterTitle();
  ax2->SetLabelOffset(2.1);
  ax2->SetTitleOffset(1.3);
  ax2->SetTitleSize(0.14);
  ax2->CenterTitle();

  double leftMarg = opt.Contains("highptfit") ? 0.2 : 0.08;
  c->cd(1);
  gPad->SetFrameLineWidth(1);
  gPad->SetLeftMargin(leftMarg);
  hf1->Draw();
  TLine zero;
  zero.DrawLine(x1, 0., x2, 0.);

  
  // NEW========================================
  for (int i=0; i<gi->GetN(); i++) {
    int q1 = GlobalIndex(i, 1), q2 = GlobalIndex(i, i);
    TF1* ifn = 0, *ifn_hi = 0, *ifn_lo = 0;
    //    TF1* _ifn = 0, *_ifn_hi = 0, *_ifn_lo = 0; // temp, unscaled.

    // This works, but has lines between \ptt groups :(
    // TH1F *hfn=0, *hfn_hi=0, *hfn_lo=0;
    // int nq = g->GetN();
    // hfn = new TH1F(Form("%s_%d", fn->GetName(), i), Form("%s_%d", fn->GetName(), i),
    // 		   nq, 0, nq);
    // hfn->SetLineColor(kMagenta);
    // for (int q=0; q<g->GetN(); q++) {
    //   hfn->SetBinContent(q+1, sf*fn->Eval(q+0.5));
    // }

    ifn = (TF1*)       fn->Clone(Form("%s_%d", fn->GetName(), i));
    ifn_hi = (TF1*) fn_hi->Clone(Form("%s_%d", fn_hi->GetName(), i));
    ifn_lo = (TF1*) fn_lo->Clone(Form("%s_%d", fn_lo->GetName(), i));

    ifn->SetName(Form("%s_%d", fn->GetName(), i));
    ifn_hi->SetName(Form("%s_%d", fn_hi->GetName(), i));
    ifn_lo->SetName(Form("%s_%d", fn_lo->GetName(), i));

    double r1 = q1-0.9, r2 = q2+0.92;
    ifn->SetRange(r1, r2);
    ifn_hi->SetRange(r1, r2);
    ifn_lo->SetRange(r1, r2);


    // for (int ip=0; ip<ifn->GetNpar(); ip++) {
    //   double par = ifn->GetParameter(ip);
    //   ifn->SetParameter(ip, 1000*par);
    //   cout<<par<< " " << ifn->GetParameter(ip) << endl;
    // }


    // ifn    = new TF1(Form("%s_%d",    fn->GetName(), i), Form("100*%s",   _ifn->GetName()));
    // ifn_hi = new TF1(Form("%s_%d", fn_hi->GetName(), i), Form("100*%s", _ifn_hi->GetName()));
    // ifn_lo = new TF1(Form("%s_%d", fn_lo->GetName(), i), Form("100*%s", _ifn_lo->GetName()));

    ifn_hi->Draw("same");
    ifn_lo->Draw("same");
    ifn->Draw("same");
    //    hfn->Draw("same");
  }
  // NEW========================================

  for (int i=0; i<gi->GetN(); i++) {
    char* newName = Form("%s_ptt%d", g->GetName(), i);
    TGraphErrors* gcol = (TGraphErrors*)g->Clone(newName);
    int col = PtColor(ptlow[i], pthigh[i]);
    utils.set_tgraph_props(gcol, kBlack, col, kFullSquare, 1.0);
    gcol->SetLineWidth(1);

    newName = Form("%s_ptt%d", gs->GetName(), i);
    TGraphErrors* g_sy = (TGraphErrors*)gs->Clone(newName);
    int col_sy = PtColorPale(ptlow[i], pthigh[i]);
    utils.set_tgraph_props(g_sy, kBlack, col_sy, kDot, 0.5);
    g_sy->SetFillColor(col_sy);
    
    for (int q=0; q<gcol->GetN(); q++) {
      int i2=0, j=0;
      ijFromGlobalIndex(q,i2,j);
      if (i2!=i) {
	gcol->SetPoint(q, gcol->GetX()[q], -99.);
	g_sy->SetPoint(q, gcol->GetX()[q], -99.);
      }

      // Scale points by sf
      gcol->SetPoint(q, gcol->GetX()[q], sf*gcol->GetY()[q]);
      g_sy->SetPoint(q, g_sy->GetX()[q], sf*g_sy->GetY()[q]);
      gcol->SetPointError(q, gcol->GetEX()[q], sf*gcol->GetEY()[q]);
      g_sy->SetPointError(q, g_sy->GetEX()[q], sf*g_sy->GetEY()[q]);

    }
    g_sy->Draw("e2psame");
    gcol->Draw("ezpsame");
  }
  
  c->cd(2);
  gPad->SetFrameLineWidth(1);
  gPad->SetLeftMargin(leftMarg);
  hf2->Draw();
  r_sys->Draw("e2psame");
  
  // Overlay lines, arrows, text, etc.
  TLine div, div_a; // division ticks for trigger and assoc
  TArrow a;
  double tipSize = 0.008; // size of arrowhead
  a.SetLineWidth(1);
  div.SetLineWidth(2);
  div.SetLineColor(fitCurveColor);
  div.DrawLine(x1, 1.0, x2, 1.0);
  r->Draw("ezpsame");
  gList->Add(r);

  for (int i=0; i<gi->GetN(); i++) {
    char* newName = Form("%s_ptt%d", r->GetName(), i);
    TGraphErrors* rcol = (TGraphErrors*)r->Clone(newName);
    int col = PtColor(ptlow[i], pthigh[i]);
    utils.set_tgraph_props(rcol, kBlack, col, kFullCircle, 1.0);
    for (int q=0; q<rcol->GetN(); q++) {
      int i2=0, j=0;
      ijFromGlobalIndex(q,i2,j);
      if (i2!=i) rcol->SetPoint(q, rcol->GetX()[q], -99.);
    }
    rcol->SetLineWidth(1);
    rcol->Draw("ezpsame");
  }
  
  div.SetLineColor(divColor);
  div.SetLineStyle(kSolid);

  // Label centrality and fourier index
  c->cd(1);
  ltx.SetTextSize(0.14);
  ltx.DrawLatex(opt.Contains("highptfit")?0.25:0.14, 0.85, Form("n = %d", n));

  if (!opt.Contains("highptfit")) {
  if (n<=2) {
    ltx.SetTextSize(0.1);
    if (ispp)
    ltx.DrawLatex(0.14, 0.70, Form("pp 2.76 TeV"));
    else {
      ltx.DrawLatex(0.14, 0.70, Form("Pb-Pb 2.76 TeV"));
      ltx.DrawLatex(0.14, 0.61, Form("%d-%d%%", centlow[k],centhigh[k]));
    }
    if (opt.Contains("ptcons")) {
    ltx.DrawLatex(0.14, 0.58, Form("#LT#sum p_{T}^{2}#GT = 4.76 GeV^{2}"));
    }

  }
  if (n==3) {
    lq->AddEntry(g, Form("V_{n#Delta}"), "epl");
    lq->AddEntry(fn, Form("fit"), "l");
    lq->Draw();
  }
  // Always show
    ltx.SetTextSize(0.07);
  double ptaFitLow = ptlow[ptaLo];
  double ptaFitHigh = ptaHi<999? pthigh[ptaHi] : pthigh[maxpt-1];
  //  if (!opt.Contains("highptfit"))
    ltx.DrawLatex(0.14, 0.09, 
		  Form("Global fit range: %3.2g < p_{T}^{a} < min(p_{T}^{t}, %3.2g GeV/c)",
		       ptaFitLow, ptaFitHigh));
  }

  for (int q=0; q<g->GetN(); q++) {
    int i=0, j=0;
    ijFromGlobalIndex(q,i,j);
    double drop = 0.35;
    double ptt = gi->GetX()[i] + gi->GetEX()[i];
    double pta = gj->GetX()[j] + gj->GetEX()[j];
    double tickLength = 0.03*(y2-y1);

    // Draw divisions between trigger pt intervals...
    bool skipTicks = opt.Contains("highptfit") && (pta < 5.);
    if (i==j && !skipTicks) {
      //      if (ptt > 5) {
	c->cd(1);
	div.DrawLine(q+1, y1, q+1, y1+tickLength);
	c->cd(2);
	tickLength = drop-0.1;
	div.DrawLine(q+1, ry2-tickLength, q+1, ry2+tickLength);
	div.DrawLine(q+1, ry1-tickLength, q+1, ry1+tickLength);
	//      }

      // Alignment: align = 10*HorizontalAlign + VerticalAlign
      // For horizontal alignment the following convention applies:
      // 1=left adjusted, 2=centered, 3=right adjusted
      // For vertical alignment the following convention applies:
      // 1=bottom adjusted, 2=centered, 3=top adjusted
      ltx2.SetTextAlign(23);
      ltx2.SetTextSize(0.12);
      ltx2.SetTextColor(divColor);

      if (ptt != 0.75) // too cluttered otherwise 
	ltx2.DrawLatex(q+1, ry1-drop, Form("%.2g", ptt) ); // pt_trig labels
      }
    
    if (pta < ptt && (pta == (int)pta || pta == 0.5) && !skipTicks) {
	div_a.DrawLine(q+1, ry1-drop/8, q+1, ry1+drop/4);
	ltx2.SetTextColor(kBlack);
	ltx2.SetTextFont(52); // italic
	ltx2.DrawLatex(q+1, ry1-drop/4, Form("%.2g", pta) ); // pt_assoc labels
	ltx2.SetTextColor(divColor);
	ltx2.SetTextFont(62); // 62 is default
    }

    // Draw arrows to off-scale points
    c->cd(1);
    double y = g->GetY()[q];
    double ymax = hf1->GetYaxis()->GetXmax();
    double ymin = hf1->GetYaxis()->GetXmin();
    double len = 0.06*(ymax-ymin);
    double gap = 0.01*(ymax-ymin);

    if (y > ymax)
      a.DrawArrow(q+0.5, ymax-len-gap, q+0.5, ymax-gap, tipSize, ">");
    if (y < ymin)
      a.DrawArrow(q+0.5, ymin+len+gap, q+0.5, ymin+gap, tipSize, ">");

    // Do the same for the ratio.
    // If fit has not been applied to some points, don't draw the arrow.
    c->cd(2);
    //    gPad->SetTickx();
    y = r->GetY()[q];
    ymax = hf2->GetYaxis()->GetXmax();
    ymin = hf2->GetYaxis()->GetXmin();
    len = 0.06/0.4*(ymax-ymin); // 0.4 from TwoPanelVert ctor
    gap = 0.02*(ymax-ymin);

    if (0) { // Draw at upper or lower edge of panel
      if (y > ymax)
	a.DrawArrow(q+0.5, ymax-len-gap, q+0.5, ymax-gap, tipSize, ">");
      if (y < ymin)
	a.DrawArrow(q+0.5, ymin+len+gap, q+0.5, ymin+gap, tipSize, ">");
    }
    // Draw with arrow based from ref. line.
    // No arrow when no fit point!
    if (j < ptaLo || j > ptaHi)
      continue;
    if (y > ymax)
      a.DrawArrow(q+0.5, 1.0, q+0.5, 1.0+len, tipSize, ">");
    if (y < ymin)
      a.DrawArrow(q+0.5, 1.0, q+0.5, 1.0-len, tipSize, ">");
  } // end q loop
  
  //  ax2->LabelsOption("h");
  c->cd();
  c->Update();
  return c;
}

TCanvas* DrawVn(int k, int ptt1, int ptt2)
{
  // Draw Fourier coefficients on one canvas for each centrality bin k
  // for trigger pt bins from ptt1 to (and including) ptt2.
  // The x-axis is associated pt.
  const char* name = Form("vndelta_etamin%02d_cent%dto%d_trig%dto%d", 				 
			  (int)(10*minRidgeDEta), centlow[k],centhigh[k], ptt1, ptt2);
  TCanvas* cv = new TCanvas(name, name, 1200, 500);
  TLegend* lv = new TLegend(0.09, 0.72, 0.27, 0.88); // w.r.t. cv
  lv->SetFillColor(kNone);
  lv->SetBorderSize(0);

  double lv2bottom = ptt2<5? 0.65 : 0.8;
  TLegend* lv2 = new TLegend(0.82, lv2bottom, 0.94, 0.95, "p_{T}^{t} (GeV/c)");
  lv2->SetFillColor(kNone);
  lv2->SetBorderSize(0);
  utils.padsetup(cv, 5, 1, "", 0.12, 0.01, 0.03, 0.2);

  // Find limits
  double minVn_k=1000, maxVn_k=-1000;	  
  VnLimits(minVn_k, maxVn_k, k, ptt1, ptt2);

    for (int n=1; n<=5; n++) {
      cv->cd(n);
      double vertMargin = 0.1*(maxVn_k - minVn_k);
      double xmax = gj->GetX()[ptt2] + gj->GetEX()[ptt2];
      double xmin =  -0.43;
      TH1F* hf = gPad->DrawFrame(xmin, minVn_k - vertMargin,
				 xmax, maxVn_k + vertMargin);
      int nXticks = 404;
      int nYticks = 210;
      if (n==1) {
	utils.make_nice_axes(cv, hf, 2., nXticks, nYticks, -0.02, 0.02);
	hf->SetLabelOffset(-0.01, "X");
      }
      else {
	utils.make_nice_axes(cv, hf, 3.4, nXticks, nYticks, -0.075, 0.0);
	hf->SetLabelOffset(-0.057, "X");
	hf->SetTickLength(0.02, "X");
	hf->SetTickLength(0.05, "Y");
      }
      hf->GetXaxis()->SetTicks("+-");
      hf->GetYaxis()->SetTicks("+-");
      
      // Loop backwards to improve visibility of low-pt points
      for (int i=ptt2; i>=ptt1; i--) {
	
	TGraphErrors *gc=0;
	gc = VnDeltaVsPtAssoc(i, k, n);
	if (!gc) {
	  Error("DrawVn()", "Problem finding graph n%d, i%d, k%d" ,n,i,k);
	  continue;
	}
	
	utils.set_tgraph_props(gc, colorsDark[i], colorsDark[i], kFullSquare, 1.6);
	
	TGraphErrors* gfit = GlobalFitVsPtAssoc(i, k, n);
	
	gfit->Draw("plsame");
	gc->Draw("epsame");

	if (n==1)
	  lv2->AddEntry(gc,trigLabel(i), "p");

      }
    }

    cv->cd();
    TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1);
    overlay->SetFillStyle(4000); // transparent
    overlay->Draw();
    overlay->cd();

    ltx.SetTextSize(0.07);
    lv2->Draw();
    ltx.DrawLatex(0.14, 0.76, Form("%s", centLabel(k)));

    for (int n=1; n<=4; n++) {
      cv->cd(n);
      ltx.SetTextSize(n==1? 0.1 : 0.15);
      ltx.DrawLatex((n==1)? 0.45 : 0.1, 0.86, Form("V_{%d#Delta}", n));
    }
    cv->cd(5);
    ltx.DrawLatex(0.72, 0.86, Form("V_{%d#Delta}", 5));
    overlay->cd();
    ltx.SetTextSize(0.075);
    ltx.DrawLatex(0.4, 0.05, "associated p_{T} [GeV/c]");
    ltx.SetTextAngle(90);
    ltx.DrawLatex(0.03, 0.4, "V_{n#Delta}(p_{T}^{assoc.})");
    ltx.SetTextAngle(0);

  return cv;
}

TCanvas* DrawVnVsPtTrig(int k, int npta, int ptabins[], TString opt)
{
  TString title = Form("VnDeltaVsPtTrig_cent%dto%d_pTassoc%.2gto%.2g",
		       centlow[k],centhigh[k],
		       ptlow[ptabins[0]], pthigh[ptabins[npta-1]]);
  if (!opt.IsNull()) {
    title.Append("_");
    title.Append(opt);
  }
  TCanvas* cv = new TCanvas(title.Data(), title.Data(), 1400, 500);
  double lv2bot = 0.9 - 0.07*npta;
  TLegend* lv2 = new TLegend(0.82, lv2bot, 0.98, 0.95, "p_{T}^{a} (GeV/c)");
  lv2->SetFillColor(kNone);
  lv2->SetBorderSize(0);
  utils.padsetup(cv, 5, 1, "", 0.12, 0.01, 0.03, 0.2);

  // Start without knowing y-limits, adjust below.
  double minVn_k=-1., maxVn_k=1.;	  
  TObjArray* graphs = new TObjArray();
  
  for (int n=1; n<=5; n++) {
    cv->cd(n);
    double vertMargin = 0.1*(maxVn_k - minVn_k);
    double xmax = 11.1;
    double xmin = -0.43;

    TH1F* hf = gPad->DrawFrame(xmin, minVn_k - 4*vertMargin,
			       xmax, maxVn_k + vertMargin);
    TLine l;
    l.DrawLine(xmin, 0, xmax, 0);
    
    int nXticks = 404;
    int nYticks = 210;
    if (n==1) {
      utils.make_nice_axes(cv, hf, 2., nXticks, nYticks, -0.02, 0.02);
      hf->SetLabelOffset(-0.01, "X");
    }
    else {
      utils.make_nice_axes(cv, hf, 3.4, nXticks, nYticks, -0.075, 0.0);
      hf->SetLabelOffset(-0.057, "X");
      hf->SetTickLength(0.02, "X");
      hf->SetTickLength(0.05, "Y");
    }
    hf->GetXaxis()->SetTicks("+-");
    hf->GetYaxis()->SetTicks("+-");

    for (int jpta=0; jpta<npta; jpta++) {
      int j = ptabins[jpta];
      TGraphErrors *gc = VnDeltaVsPtTrig(j, k, n);
      if (!gc) {
	Error("DrawVnVsPtTrig()", "Problem finding gc n%d, i%d, k%d" ,n,j,k);
	continue;
      }
      TGraphErrors *gs = VnDeltaVsPtTrig(j, k, n, "sys");
      if (!gs) {
	Error("DrawVnVsPtTrig()", "Problem finding gs n%d, i%d, k%d" ,n,j,k);
	continue;
      }

      int col = PtColor(ptlow[j], pthigh[j]);
      utils.set_tgraph_props(gc, col, col, kFullSquare, 1.2);
      gs->SetFillColor(PtColorPale(ptlow[j], pthigh[j]));
      TGraphErrors* open = (TGraphErrors*)gc->Clone();
      utils.set_tgraph_props(open, kBlack, kBlack, kOpenSquare, 1.2);
	
      gs->Draw("e2psame");
      gc->Draw("epsame");
      open->Draw("epsame");

      if (n==1)
	lv2->AddEntry(gc,asscLabel(j), "elp");

      graphs->Add(gc); // store graph array to find y limits for plot
    }
  }

  // Adjust y limits
  double newYmin = 0, newYmax = 0, newSpace = 0;
  for (int m=0; m<graphs->GetEntries(); m++) {
    TGraphErrors* tg = (TGraphErrors*)graphs->At(m);
    double yhi = TMath::MaxElement(tg->GetN()-1, tg->GetY());
    double ylo = TMath::MinElement(tg->GetN()-1, tg->GetY());
    double ey = TMath::MaxElement(tg->GetN()-1, tg->GetEY());
    ylo -= ey;
    if (yhi > newYmax) newYmax = yhi;
    if (ylo < newYmin) newYmin = ylo;
  }
  newSpace = 0.1*(newYmax - newYmin);
  for (int n=1; n<=5; n++) {
    cv->cd(n);
    TH1F* hf = (TH1F*)gPad->FindObject("hframe");
    hf->GetYaxis()->SetRangeUser(newYmin-3*newSpace, newYmax+newSpace);
    gPad->Modified(); // Signal to redraw
  }
 
  cv->cd();
  TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1);
  overlay->SetFillStyle(4000); // transparent
  overlay->Draw();
  overlay->cd();
  lv2->Draw();
  ltx.SetTextSize(0.05);
  ltx.DrawLatex(0.14, 0.88, Form("%s 2.76 TeV", ispp?"pp":"Pb-Pb"));
  ltx.SetTextSize(0.07);
  if (!ispp)
    ltx.DrawLatex(0.14, 0.8, Form("%s", centLabel(k)));

  for (int n=1; n<=5; n++) {
    cv->cd(n);
    ltx.SetTextSize(n==1? 0.1 : 0.15);
    if (n==1)
      ltx.DrawLatex(0.44, 0.24, Form("n=%d", n));
    else
      ltx.DrawLatex(0.04, 0.24, Form("n=%d", n));

  }
  overlay->cd();
  ltx.SetTextSize(0.075);
  ltx.DrawLatex(0.48, 0.05, "trigger p_{T}^{t} [GeV/c]");
  ltx.SetTextAngle(90);
  ltx.DrawLatex(0.06, 0.48, "V_{n#Delta}  [10^{-2}]");
  ltx.SetTextAngle(0);
  return cv;
}

TGraphErrors* AgreementVsPtSum(int k, int n)
{
  TGraphErrors* gVn = VnDVsQ(n, k);
  TF1* fVn = GlobalFitVsQ(n,k);
  TGraphErrors* g = new TGraphErrors();
  int N = gVn->GetN(); // And f better have the same N.
  int i = 0, j = 0;
  if(N<=0) {
    Error("AgreementVsPtSum()", "Invalid number of points: %d", N);
    gSystem->Exit(1);
  }

  for (int q=0; q<N; q++) {
    ijFromGlobalIndex(q,i,j);
    double gy = gVn->GetY()[q];
    double fy = fVn->Eval(q+0.5);
    double r = TMath::Abs((gy-fy)/gy); 
    g->SetPoint(q, PtBinCenterSum(i, j), r);
  }
  g->SetMarkerStyle(kFullCircle);
  return g;
}

TGraph2D* Agreement2D(int k, int n)
{
  TGraphErrors* gVn = VnDVsQ(n, k);
  TF1* fVn = GlobalFitVsQ(n,k);
  TGraph2D* g = new TGraph2D();
  int N = gVn->GetN(); // And f better have the same N.
  int i = 0, j = 0;
  if(N<=0) {
    Error("AgreementVsPtSum()", "Invalid number of points: %d", N);
    gSystem->Exit(1);
  }

  for (int q=0; q<N; q++) {
    ijFromGlobalIndex(q,i,j);
    double gy = gVn->GetY()[q];
    double fy = fVn->Eval(q+0.5);
    double r = TMath::Abs((gy-fy)/fy); 
    g->SetPoint(q, ptmean[j], ptmean[i], r);
  }
  return g;
}

TH2F* Agreement2DHist(int k, int n)
{
  TGraphErrors* gVn = VnDVsQ(n, k);
  TF1* fVn = GlobalFitVsQ(n,k);

  // Pt binning array
  double ptarr[100];
  for (int i=0; i<maxpt; i++) ptarr[i] = gi->GetX()[i] - gi->GetEX()[i];
  ptarr[maxpt] = gi->GetX()[maxpt-1] + gi->GetEX()[maxpt-1];

  int N = gVn->GetN(); // And f better have the same N.
  int i = 0, j = 0;
  if(N<=0) {
    Error("AgreementVsPtSum()", "Invalid number of points: %d", N);
    gSystem->Exit(1);
  }
  const char* name = Form("agmt_cent%dto%d_n%d", 
			  centlow[k],centhigh[k], n);
  TH2F* h = new TH2F(name, name, maxpt, ptarr, maxpt, ptarr);
  h->SetTitle("title;p_{T}^{a};p_{T}^{t};");

  for (int q=0; q<N; q++) {
    ijFromGlobalIndex(q,i,j);
    double gy = gVn->GetY()[q];
    double fy = fVn->Eval(q+0.5);
    
    // double r = fy==0? 0 : TMath::Abs(gy-fy) / fy; 
    double r = gy==0? 0 : TMath::Abs((gy-fy) / gy); 
    //    double r = fy==0? 0 : TMath::Abs(gy / fy); 
    h->SetBinContent(j+1, i+1, r);
  }
  return h;
}

TH2* PttPtaHist()
{
  // Pt binning array
  double ptarr[100];
  for (int i=0; i<maxpt; i++) ptarr[i] = gi->GetX()[i] - gi->GetEX()[i];
  ptarr[maxpt] = gi->GetX()[maxpt-1] + gi->GetEX()[maxpt-1];

  int i = 0, j = 0;
  const char* name = Form("PttPtaHist");
  TH2F* h = new TH2F(name, name, maxpt, ptarr, maxpt, ptarr);
  h->SetTitle("title;p_{T}^{a} (GeV/c);p_{T}^{t} (GeV/c);");
  h->GetXaxis()->CenterTitle();
  h->GetYaxis()->CenterTitle();
  int N = maxpt * (maxpt + 1) / 2;
  for (int q=0; q<N; q++) {
    ijFromGlobalIndex(q,i,j);
    h->SetBinContent(j+1, i+1, 1.);
  }
  return h;
}

TGraphErrors* Chi2VsTrigPt(int k, int n)
{
  TGraphErrors* g = new TGraphErrors();
  for (int i=0; i<maxpt; i++) {
    g->SetPoint(i, ptmean[i], ReducedChi2(i, k, n));
  }
  g->SetMarkerStyle(kFullCircle);
  return g;
}

TGraphErrors* 
CorrFnChi2VsTrigPt(int j, int k, int n1, int n2, TString opt)
{
  TGraphErrors* g = new TGraphErrors();
  int ip=0;
  for (int i=j; i<maxpt; i++) {
    TH1* h = Hist("RIDGE", "cA", i, j, k);
    double x = MeanPt(i, j, k, "t");
    double y = Chi2(h, HarmonicSum(h, n1, n2, "global"), opt);
    g->SetPoint(ip, x, y);
    ip++;
  }
  g->SetMarkerStyle(kFullCircle);
  return g;
}

double ReducedChi2(int i, int k, int n)
{
  TGraphErrors* g = VnDeltaVsPtAssoc(i, k, n);
  TGraphErrors* f = GlobalFitVsPtAssoc(i, k, n);
  int N = g->GetN(); // And f better have the same N.
  double chi2 = 0;

  // Value (or one-sigma error if opt=="err") of v_n
  //  double vnerr = Bestvn(k,n,i,"err");
  double vnerr = vnGF(n,k,i,0,999,"RIDGE","cA","err");

  if(N<=0 || N!= f->GetN()) {
    Error("ReducedChi2()", "Invalid number of points: %d", N);
    gSystem->Exit(1);
  }

  for (int j=0; j<N; j++) {
    double errj = g->GetEY()[j];

    chi2 += 
      (g->GetY()[j] - f->GetY()[j]) * (g->GetY()[j] - f->GetY()[j]) / 
      pow(errj + vnerr, 2);
  }
  chi2 /= N;
  return chi2;
}

TGraphErrors* VnDeltaVsPtAssoc(int i, int k, int n, TString opt)
{
  // Create a TGraph from VnDelta points
  TGraphErrors* gALL = VnDVsQ(n, k, "RIDGE", "cA", opt);
  TGraphErrors* g = new TGraphErrors();
  TString gname(Form("VnDeltaVsPtAssoc_ptt%d_cent%d_n%d",i,k,n));
  gname.Append(opt);
  g->SetName(gname.Data());

  for (int j=0; j<=i; j++) {
    int q = GlobalIndex(i, j);
    g->SetPoint(j, gj->GetX()[j], gALL->GetY()[q] );
    g->SetPointError(j, 0,  gALL->GetEY()[q] );
    utils.set_tgraph_props(g, colorsDark[i], colorsDark[i], kFullSquare, 1.0); 
    g->SetLineWidth(2);
  }
  
  return g;
}

TGraphErrors* VnDeltaVsPtTrig(int j, int k, int n, TString opt)
{
  // Create a TGraph from VnDelta points
  TGraphErrors* gALL = VnDVsQ(n, k, "RIDGE", "cA", opt);
  TGraphErrors* g = new TGraphErrors();
  TString gname(Form("VnDeltaVsPtTrig_pta%d_cent%dto%d_n%d",
		     j,centlow[k],centhigh[k],n));
  gname.Append(opt);
  g->SetName(gname.Data());

  int ip=0;
  for (int i=j; i<maxpt; i++) {
    int q = GlobalIndex(i, j);
    
    double xerr = 0.;
    if (opt.Contains("sys")) {
      if (gi->GetEX()[i] < 0.126)
	xerr = gi->GetEX()[i];
      else
	xerr = 0.250;
    }
    double mpt = MeanPt(i, j, k, "t"); 
    double sf = 100; // Scale y-value; make sure to label accordingly!!
    g->SetPoint(ip, mpt, sf*gALL->GetY()[q] ); 
    g->SetPointError(ip, xerr,  sf*gALL->GetEY()[q] );

    utils.set_tgraph_props(g, colorsDark[j], colorsDark[j], kFullSquare, 1.0); 
    g->SetLineWidth(2);
    ip++;
  }
  
  return g;
}

TGraphErrors* VnDeltaVsN(int i, int j, int k, int nMax, TString opt)
{
  // Create a TGraph from VnDelta points
  // Stat errors by default. Call with opt "sys" to get sys errors.
  // Call with opt "sine" to get residual <sin n Delta phi> for sys. est.
  int q = GlobalIndex(i, j);
  TString gname(Form("VnDeltaVsN_ptt%.2gto%.2g_pta%.2gto%.2g_cent%dto%d",
		     ptlow[i],pthigh[i],ptlow[j],pthigh[j],
		     centlow[k],centhigh[k]));
  if (!opt.IsNull())  {
    gname.Append("_");
    gname.Append(opt);
  }
  TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data());
  if (g) {
    return g;
  } 
  else {
    g = new TGraphErrors();
    g->SetName(gname.Data());
  }

  int ip=0;
  gStyle->SetEndErrorSize(5);
  for (int n=1; n<=nMax; n++) {

    TGraphErrors* gALL = VnDVsQ(n, k, "RIDGE", "cA", opt);
    double xerr = opt.Contains("sys") ? 0. : 0.0; // 0.2
    g->SetPoint(ip, n, 100*gALL->GetY()[q] );
    g->SetPointError(ip, xerr,  100*gALL->GetEY()[q] );
    if (0)
      cout<<gALL->GetName()<< " " << g->GetY()[ip] << endl;

    ip++;
  }

  utils.set_tgraph_props(g, kBlack, kBlack, kFullCircle, 1.2); 
  g->SetLineWidth(2);
  gList->Add(g);
  return g;
}

TGraphErrors* GlobalvnVsN(int i, int k, int nMax, TString opt)
{
  // Create a TGraph from vn{GF} points
  // Stat errors by default. Call with opt "sys" to get sys errors.

  TString gname(Form("GlobalvnVsN_pt%.2gto%.2g_cent%dto%d",
		     ptlow[i],pthigh[i], centlow[k],centhigh[k]));
  if (!opt.IsNull())  {
    gname.Append("_");
    gname.Append(opt);
  }
  TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data());
  if (g) {
    return g;
  } 
  else {
    g = new TGraphErrors();
    g->SetName(gname.Data());
  }
  g->SetName(gname.Data());

  int ip=0;
  for (int n=1; n<=nMax; n++) {

    // point, stat err, sys err.
    double vnp = vnGF(n, k, i, 0, 999, "RIDGE", "cA", "");
    double vne = vnGF(n, k, i, 0, 999, "RIDGE", "cA", "err");
    double vns = vnGF(n, k, i, 0, 999, "RIDGE", "cA", "sys");
    double xerr = opt.Contains("sys") ? 0. : 0.0; // 0.2
    double yerr = opt.Contains("sys") ? vns : vne;

    g->SetPoint(ip, n, vnp );
    g->SetPointError(ip, xerr, yerr);
    ip++;
  }

  utils.set_tgraph_props(g, kBlack, kBlack, kFullCircle, 1.2); 
  g->SetLineWidth(2);

  return g;
}

TGraphErrors* VnDeltaNFVsN(int i, int j, int k, int nMax, TString opt)
{

  TGraphErrors* VnD = VnDeltaVsN(i,j,k,nMax,opt);
  TGraphErrors* vnt = GlobalvnVsN(i,k,nMax,opt);
  TGraphErrors* vna = GlobalvnVsN(j,k,nMax,opt);
  TString gname(Form("VnDeltaNFVsN_ptt%.2gto%.2g_pta%.2gto%.2g_cent%dto%d",
		     ptlow[i],pthigh[i],ptlow[j],pthigh[j],
		     centlow[k],centhigh[k]));
  if (!opt.IsNull())  {
    gname.Append("_");
    gname.Append(opt);
  }
  TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data());
  if (g) {
    return g;
  } 
  else {
    g = new TGraphErrors();
    g->SetName(gname.Data());
  }
  g->SetName(gname.Data());

  int ip=0;
  for (int n=1; n<=nMax; n++) {

    // Uncertainty vne is stat or sys depending on opt.
    double vnp = VnD->GetY()[ip] - 100*vnt->GetY()[ip]*vna->GetY()[ip];
    double vne = VnD->GetEY()[ip] - 100*vnt->GetEY()[ip]*vna->GetEY()[ip];
    g->SetPoint(ip, n, vnp );
    g->SetPointError(ip, 0.0, vne);
    ip++;
  }

  utils.set_tgraph_props(g, kBlack, kBlack, kOpenSquare, 1.2); 
  g->SetLineWidth(2);
  return g;
}


TF1* GlobalFitVsQ(int n, int k, int ipt1, int ipt2, const char* region, 
		  const char* corrtype, TString opt)
{

  // vn{GF}(trig) x vn{GF}(assc) vs. q, where the points between (and
  // including) ipt1 and ipt2 were used in the global fit.
  // Note: gfname must be identical format to fnName in GlobalFitVsQ().
  TF1* fitfn = 0;
  TString errType("meas");
  if (opt.Contains("hi_sys"))
    errType = "hi";
  if (opt.Contains("lo_sys"))
    errType = "lo";
  TString gfname = 
    Form("globalFitFn_cent%dto%d_n%d_ptbins%dto%d_%s_corrtype_%s_%s", 
	 centlow[k],centhigh[k], n, ipt1, ipt2, 
	 region, corrtype, errType.Data());
  fitfn = (TF1*)gFitList->FindObject(gfname.Data());
  if (fitfn) {
    if (0)
      Info("GlobalFitVsQ()", "Returning cached TF1 %s", gfname.Data());
    return fitfn;
  }

  // If nonexistent, create (& store) a new one
  fitfn = DoGlobalFit(n,k,ipt1,ipt2,region,corrtype,opt);
  //  fitfn = (TF1*)gFitList->FindObject(gfname.Data());
  
  if (!fitfn)
    Warning("GlobalFitVsQ()", 
	    "no fitfn %s,\neven after calling DoGlobalFit()", 
	    gfname.Data());

  return fitfn;
}

TF1* DoGlobalFit(int n, int k, int ipt1, int ipt2, 
		 const char* region, const char* corrtype, TString opt)
{
  // fnName must be identical format to gfname in GlobalFitVsQ().
  TF1 *globalFitFn = 0;
  TGraphErrors* g = 0;
  TString errType("meas");
  if (opt.Contains("hi_sys"))
    errType = "hi";
  if (opt.Contains("lo_sys"))
    errType = "lo";
  TString fnName = 
    Form("globalFitFn_cent%dto%d_n%d_ptbins%dto%d_%s_corrtype_%s_%s", 
	 centlow[k],centhigh[k], n, ipt1, ipt2, 
	 region, corrtype, errType.Data());

  g = VnDVsQ(n, k, region, corrtype, opt);

  if (!g) {
    Error("DoGlobalFit()", "No graph found");
    return 0;
  }

  // Assign global variables that will get checked in the vntvna
  // function during the fit.
  minFitPtBin = ipt1;
  maxFitPtBin = ipt2;

  // Set up the global fit function...  
  double qmax = g->GetN() - 0.4999;
  int nPars = maxpt + 3;
  globalFitFn = new TF1(fnName.Data(), vntvna, 0.0, qmax, nPars);
  
  /* Not needed now, I just flip the v1 points later...AMA 11/15/2011 */
  // Each parameter has an equally good minimum at + or - v_nbest. The
  // starting value is chosen to steer the fit to the positive "root".
  // v_1 is special because it is the only coefficient that crosses
  // zero at low pt, so it gets a negative starting point. Pars are
  // vn(ipt), then last one is n.
  //  double startVal = (n==1)? -0.02 : 0.01;

  double startVal = 0.01;
  for (int ip=0; ip<maxpt; ip++) { 
    globalFitFn->SetParameter(ip, startVal);
    globalFitFn->SetParName(ip, Form("v_%d(%.3g)",n,ptmean[ip]));
  }

  // Pass n into vntvna(). Access as par[maxpt].
  globalFitFn->FixParameter(maxpt, double(n));

  // Pass k into vntvna(). Access as par[maxpt+1].
  globalFitFn->FixParameter(maxpt+1, double(k));

  // Mom. conservation term prefactor. Access as par[maxpt+2].
  if (n==1) {
    globalFitFn->FixParameter(maxpt+2, 0.);
    //    globalFitFn->SetParLimits(maxpt+2, 0., 1e-3);
  }
  else
    globalFitFn->FixParameter(maxpt+2, 0.);
  
  globalFitFn->SetNpx(1000);

  TFitResultPtr result = g->Fit(globalFitFn,"QRNS"); // fast, takes ~15 ms

  if (0) {
    cout<<globalFitFn->GetName();
    result->Print();
  }
  // TMinuit status: 
  // migradResult + 10*minosResult + 100*hesseResult + 1000*improveResult
  // 0:ok, <0:user error, >0: see above.
  int fitStatus = result;
  if (fitStatus)
    Warning("DoGlobalFit()", "TMinuit status %d", fitStatus);

  gFitList->Add(globalFitFn);

  return globalFitFn;
}

double vnGF(int n, int k, int ptBin, int ipt1, int ipt2,
	    const char* region, const char* corrtype, TString opt)
{

  double vn, err;
  TF1* gfn = GlobalFitVsQ(n,k,ipt1,ipt2,region,corrtype,opt);
  if (!gfn) 
    Error("vnGF()", "!gfn");  
  vn = gfn->GetParameter(ptBin); // ptBin is the q index
  err = gfn->GetParError(ptBin);
  
  if (opt.Contains("sys")) {
    TF1* hi = GlobalFitVsQ(n,k,ipt1,ipt2,region,corrtype,"hi_sys");
    TF1* lo = GlobalFitVsQ(n,k,ipt1,ipt2,region,corrtype,"lo_sys");
    double vhi = hi->GetParameter(ptBin);
    double vlo = lo->GetParameter(ptBin);
    return 0.5*TMath::Abs(vhi-vlo);
  }

  // v1 can be negative, otherwise ensure the positive solution is chosen
  // if (n>1 && vn < 0)
  //   vn *= -1;

  // if (ptlow[ipt1] > 4.) // !!!!!!!!!!!!!!!!!!!!!!!!!
  //   vn = TMath::Abs(vn);

  // if (n==1) 
  //   vn += MomCorrection(i, j, k);

  return opt.Contains("err") ? err : vn;
}

TGraphErrors* vnGFvsPt(int n, int k, int ipt1, int ipt2, 
		       const char* region, const char* corrtype, TString opt)
{

  // vn{GF} vs. pt, where the points between (and including) ipt1 and
  // ipt2 were used in the global fit.
  bool isHighPtFit = ipt1 > 0 && ipt2==999;
  bool flipSign = n==1;  // Flip v1?
  if (isHighPtFit) 
    flipSign = 0;
  if (opt.Contains("keepsign"))
    flipSign = 0;
  TString name(Form("v%dGFvsPt_fitpt%dto%d_cent%dto%d",
		  n, ipt1, ipt2, centlow[k], centhigh[k]));
  if (opt.Contains("sys"))
    name.Append("_sys");
  if (opt.Contains("ptcons"))
    name.Append("_ptcons");

  TGraphErrors* g = (TGraphErrors*)gList->FindObject(name.Data());
  if (g)
    return g;
  
  g = new TGraphErrors();
  g->SetName(name);
  
  for (int i=0; i<maxpt; i++) {
    TString opt1 = "";
    if (opt.Contains("ptcons"))
      opt1.Append("_ptcons");
    double vn = vnGF(n,k,i,ipt1,ipt2,region,corrtype,opt1);
    
    if (flipSign)
      vn = -vn;
    
    //    if (isHighPtFit && i<ipt1)
    if (i<ipt1)
      vn = 0.;

    double mpt = MeanPt(i, 0, k, "t"); 
    if (i<ipt1) 
      mpt = -10.;
    g->SetPoint(i, mpt, vn);
    
    if (opt.Contains("sys")) {
      opt1.Append("_sys");
      double ex  = (gi->GetEX()[i] < 0.126) ? gi->GetEX()[i] : 0.250;
      double vns = vnGF(n,k,i,ipt1,ipt2,region,corrtype, opt1);
      g->SetPointError(i, ex, vns);
      //      g->SetPointError(i, gi->GetEX()[i], vns); // too wide 
    }
    else {
      opt1.Append("_err");
      double vne = vnGF(n,k,i,ipt1,ipt2,region,corrtype,opt1);
      g->SetPointError(i, 0.0, vne);
    } 
  }

  int col = CentColor(centlow[k], centhigh[k]);
  int col_pale = CentColorPale(centlow[k], centhigh[k]);
  g->SetLineWidth(2);
  g->SetFillColor(col_pale); // for sys
  utils.set_tgraph_props(g, col, col, kFullCircle, 1.2);

  if (0) cout<<opt.Data()<<endl;
  gList->Add(g);
  return g;
}

TGraph* Luzumv1(int cent, const char* pid) 
{
  // v1 for pions vs pt from viscous hydro calc by Matt Luzum
  // valid cent args 0-4 - 10% bins.
  // pid = "pi", "K", "p".
  const char* dat = Form("v1LHCGlauber%ssv08.dat", pid);
  TString fmt("%lg "); // pt
  for (int ky=0; ky<5; ky++) 
    fmt.Append( ky==cent? "%lg " : "%*lg "); // centrality - skip if *lg
  return new TGraph(dat, fmt.Data());
}

TGraphErrors* GlobalFitVsPtAssoc(int i, int k, int n, TString opt)
{
  // Create a TGraph of vnt{gf} x vna{gf} vs. pt_assoc from gfit
  TF1* gfit = GlobalFitVsQ(n, k, 0, 999, "RIDGE", cfDef, opt);
  TGraphErrors* g = new TGraphErrors();
  TString gname(Form("GlobalFitVsPtAssoc_ptt%d_cent%d_n%d",i,k,n));
  if (!opt.IsNull())
    gname.Append(opt);
  g->SetName(gname);

  for (int j=0; j<=i; j++) {
    int q = GlobalIndex(i, j);
    g->SetPoint(j, gj->GetX()[j], gfit->Eval(q + 0.5) );
    g->SetPointError(j, 0, 0 );
    utils.set_tgraph_props(g, colorsPale[i], colorsPale[i], kFullCircle, 1.0); 
    g->SetLineWidth(4);
  }
  return g;
}

double SineError(int i, int j, int k)
{
  double rms = 0;
  double resid_sin=0, rs_err=0;
  TH1* hst  = Hist("RIDGE", "cA", i, j, k, ""); 	
  double arr[100] = {0};
  //  arr[0] = 0;
  for (int n=1; n<=gNMax; n++) {
    VnDelta(n, *hst,  resid_sin, rs_err, "sin");
    arr[n-1] = resid_sin;
  }

  rms = TMath::RMS(gNMax, arr) / 2;
  if (0)
    cout << " SineError: " << rms << endl;

  return rms;
}

TCanvas* DrawVnDeltaVsN(int nptt, int pttbins[], int npta, 
			int ptabins[], int ncb, int centbins[], TString opt)
{

  double t1=ptlow[pttbins[0]], t2=pthigh[pttbins[nptt-1]];
  double a1=ptlow[ptabins[0]], a2=pthigh[ptabins[npta-1]];
  TString cname = Form("VnDeltaVsN_ptt%.2gto%.2g_pta%.2gto%.2g",t1,t2,a1,a2);
  if (!opt.IsNull()) {
    cname.Append("_");
    cname.Append(opt);
  }
  TCanvas* c = new TCanvas(cname.Data(), cname.Data(), 650,700);
  double legbottom = 0.9 - 0.07*ncb;
  TLegend* leg = new TLegend(0.62, legbottom, 0.98, 0.98, "Centrality");
  leg->SetFillColor(kNone);
  leg->SetBorderSize(0);
  leg->SetName("leg");

  TObjArray* graphs = new TObjArray();
  double xmin = 0.4;
  double xmax = opt.Contains("ext") ? gNMax + 0.5 : 8.5;
  int nXticks =  10, nYticks = 210;
  TLine l;
  TH1F* hf = gPad->DrawFrame(xmin, 0.0001, xmax, 0.26);
  //  bool scale1000 = 0;

  hf->SetTitle(Form(";n;V_{n#Delta} [10^{-2}]"));
  if (opt.Contains("sine"))
    hf->SetTitle(Form(";n;#LTsin(n#Delta#phi)#GT [10^{-2}]"));

  // if (opt.Contains("ext") )  
  //   if (!opt.Contains("sine")) {
  //     scale1000 = 0;
  //     hf->SetTitle(Form(";n;V_{n#Delta}"));
  // }

  l.SetLineStyle(kDashed);

  utils.make_nice_axes(c, hf, 1.3, nXticks, nYticks, 1.4,-1.4);
  hf->SetLabelOffset(0.0, "X");
  hf->SetTitleOffset(-1.96, "Y");

  hf->GetXaxis()->SetTicks("+-");
  hf->GetYaxis()->SetTicks("+-");
  l.DrawLine(hf->GetXaxis()->GetXmin(), 0, hf->GetXaxis()->GetXmax(),0);
  gPad->SetLeftMargin(0.22);
  gPad->SetTopMargin(0.01);
  gPad->SetRightMargin(0.01);

  for (int i=0; i<nptt; i++) {
    for (int j=0; j<npta; j++) {
      for(int k=0; k<ncb; k++) {

	// Stat errors on g, sys errors on gs.
	TString opt1 = opt.Contains("sine") ? "sine" : "";
	TGraphErrors* g = 
	  VnDeltaVsN(pttbins[i], ptabins[j], centbins[k], gNMax, opt1);
	TGraphErrors* gs = 
	  VnDeltaVsN(pttbins[i], ptabins[j], centbins[k], gNMax, "sys");
	TGraphErrors* gc = (TGraphErrors*)g->Clone();
	int dark = CentColor(centlow[centbins[k]], centhigh[centbins[k]]);
	int pale = CentColorPale(centlow[centbins[k]], centhigh[centbins[k]]);
	utils.set_tgraph_props(g, dark, dark, kFullCircle, 1.5); 
	utils.set_tgraph_props(gs, dark, dark, kDot, 2.); 
	utils.set_tgraph_props(gc, kBlack, kBlack, kOpenCircle, 1.5);
	g->SetFillColor(pale);
	gs->SetFillColor(dark);

	TGraphErrors* gnf = 0;
	TGraphErrors* gnf_sys = 0;
	TGraphErrors* gnf_solid = 0;
	if (opt.Contains("vnf")) {  // Non-factorizing part
	  gnf = 
	    VnDeltaNFVsN(pttbins[i], ptabins[j], centbins[k], gNMax, "");
	  gnf_sys = 
	    VnDeltaNFVsN(pttbins[i], ptabins[j], centbins[k], gNMax, "sys");

	  // Just for appearance...
	  gnf_solid = (TGraphErrors*) gnf->Clone(Form("%s_solid", gnf->GetName()));
	  utils.set_tgraph_props(gnf_solid, pale, pale, kFullSquare, gnf->GetMarkerSize());
	}
	// if (scale1000) {
	//   for (int ip=0; ip<g->GetN(); ip++) {
	//     g->SetPoint(ip,   g->GetX()[ip], 1000* g->GetY()[ip]);
	//     gs->SetPoint(ip, gs->GetX()[ip], 1000*gs->GetY()[ip]);
	//     gc->SetPoint(ip, gc->GetX()[ip], 1000*gc->GetY()[ip]);

	//     g->SetPointError(ip,   g->GetEX()[ip], 1000* g->GetEY()[ip]);
	//     gs->SetPointError(ip, gs->GetEX()[ip], 1000*gs->GetEY()[ip]);
	//     gc->SetPointError(ip, gc->GetEX()[ip], 1000*gc->GetEY()[ip]);
	//   }
	// }

	if (opt.Contains("sine")) {
	  g->Fit("pol0", "Q");
	  TF1* gfit = g->GetFunction("pol0");
	  gfit->SetLineColor(dark);
	}

	else if (!opt.Contains("nobars"))
	  g->Draw("Bsame");  // bars

	if (opt.Contains("vnf")) {  // Non-factorizing part
	  utils.draw_errorbox(gnf_sys, pale, 0.3);
	  gnf_solid->Draw("epsame");
	  gnf->Draw("epsame");
	}

	g->Draw("psame");  // colored markers
	gc->Draw("epsame");  // open markers w/black stat error bars

	if (!opt.Contains("sine"))
	  utils.draw_errorbox(gs, dark, 0.3);

	graphs->Add(g);
	leg->AddEntry(g,Form("%d-%d%%",centlow[centbins[k]],centhigh[centbins[k]]),"l,p");

	if (opt.Contains("vnf")) {  // Non-factorizing part
	  leg->AddEntry(gnf_solid,"NF V_{n#Delta}","l,p");
	  //	  leg->AddEntry(gnf_solid,"[v_{n}(p_{T}^{t}) #times v_{n}(p_{T}^{a})]_{GF}","l,p");
	  leg->SetHeader(""); 
	}

      }
    }
  }

  // Adjust y limits
  double newYmin = 0, newYmax = 0, newSpace = 0;
  for (int m=0; m<graphs->GetEntries(); m++) {
    TGraphErrors* tg = (TGraphErrors*)graphs->At(m);
    double yhi = TMath::MaxElement(tg->GetN()-1, tg->GetY());
    double ylo = TMath::MinElement(tg->GetN()-1, tg->GetY());
    if (yhi > newYmax) newYmax = yhi;
    if (ylo < newYmin) newYmin = ylo;
  }

  newSpace = 0.1*(newYmax - newYmin);

  if (opt.Contains("ext"))
    hf->GetYaxis()->SetRangeUser(newYmin-newSpace, newYmax+4*newSpace);
  else
    hf->GetYaxis()->SetRangeUser(newYmin-newSpace, newYmax+newSpace);

  // if (opt.Contains("sine")) {
  //   hf->GetYaxis()->SetRangeUser(0.0001, newYmax+4*newSpace);
  //   gPad->SetLogy();
  // }

  if (!opt.Contains("sine")) 
    leg->Draw();
 
  // Draw text
  ltx.SetTextSize(0.04); // 0.05

  if ( opt.Contains("sine") )  
    legbottom = 0.99;
  else if (opt.Contains("ext") && !opt.Contains("vnf") )  
    legbottom -= 0.4;
  
  double xstart = 0.6;
  ltx.DrawLatex(xstart, legbottom-1*0.07, Form("%.2g < p_{T}^{t} < %.2g GeV/c", t1, t2));
  ltx.DrawLatex(xstart, legbottom-2*0.07, Form("%.2g < p_{T}^{a} < %.2g GeV/c", a1, a2));
  
  gPad->Update();
  gPad->Modified();
  return c; 
}

TCanvas* DrawGlobalvnVsN(int npt, int ptbins[], int ncb, int centbins[], TString opt)
{
  double t1=ptlow[ptbins[0]], t2=pthigh[ptbins[npt-1]];
  TString cname = Form("GlobalvnVsN_pt%.2gto%.2g",t1,t2);
  if (!opt.IsNull()) {
    cname.Append("_");
    cname.Append(opt);
  }
  TCanvas* c = new TCanvas(cname.Data(), cname.Data(), 600,700);
  double legbottom = 0.9 - 0.07*ncb;
  TLegend* leg = new TLegend(0.62, legbottom, 0.98, 0.98, "Centrality");
  leg->SetFillColor(kNone);
  leg->SetBorderSize(0);
  leg->SetName("leg");

  TObjArray* graphs = new TObjArray();
  double xmin = opt.Contains("zoom") ? 2.6 : 0.4;
  double xmax = 8.5;
  int nXticks =  10, nYticks = 210;
  TLine l;
  TH1F* hf = gPad->DrawFrame(xmin, -0.06, xmax, 0.26);
  hf->SetTitle(Form(";n;v_{n}{GF}"));
  l.SetLineStyle(kDashed);

  utils.make_nice_axes(c, hf, 1.3, nXticks, nYticks, 1.4,-1.4);
  hf->SetLabelOffset(0.0, "X");
  if (opt.Contains("zoom"))
    hf->SetTitleOffset(-2.3, "Y");
  else
    hf->SetTitleOffset(-1.9, "Y");

  hf->GetXaxis()->SetTicks("+-");
  hf->GetYaxis()->SetTicks("+-");
  l.DrawLine(hf->GetXaxis()->GetXmin(), 0, hf->GetXaxis()->GetXmax(),0);
  if (opt.Contains("zoom"))
    gPad->SetLeftMargin(0.25);
  else
    gPad->SetLeftMargin(0.2);
  gPad->SetTopMargin(0.01);
  gPad->SetRightMargin(0.01);

  for (int i=0; i<npt; i++) {
    for(int k=0; k<ncb; k++) {

      // Stat errors on g, sys errors on gs.
      TGraphErrors* g  = GlobalvnVsN(ptbins[i], centbins[k], gNMax, "");
      TGraphErrors* gs = GlobalvnVsN(ptbins[i], centbins[k], gNMax, "sys");
      TGraphErrors* gc = (TGraphErrors*)g->Clone();
      int dark = CentColor(centlow[centbins[k]], centhigh[centbins[k]]);
      int pale = CentColorPale(centlow[centbins[k]], centhigh[centbins[k]]);
      utils.set_tgraph_props(g, dark, dark, kFullCircle, 1.5); 
      //      utils.set_tgraph_props(gs, dark, dark, kDot, 2.); 
      utils.set_tgraph_props(gc, kBlack, kBlack, kOpenCircle, 1.5);
      g->SetFillColor(pale);
      //      gs->SetFillColor(dark);
      g->Draw("Bsame");  // colored bars
      utils.draw_errorbox(gs, dark, 0.3);
      g->Draw("psame");  // colored markers
      //      gs->Draw("e[]psame"); // sys error filled rectangles
      gc->Draw("epsame");  // open markers w/black stat error bars
      graphs->Add(g);
      leg->AddEntry(g,Form("%d-%d%%",centlow[centbins[k]],centhigh[centbins[k]]),"l,p");
    }
  }

  // Adjust y limits
  double newYmin = 0, newYmax = 0, newSpace = 0;
  for (int m=0; m<graphs->GetEntries(); m++) {
    TGraphErrors* tg = (TGraphErrors*)graphs->At(m);
    double yhi = TMath::MaxElement(tg->GetN()-1, tg->GetY());
    double ylo = TMath::MinElement(tg->GetN()-1, tg->GetY());
    if (yhi > newYmax) newYmax = yhi;
    if (ylo < newYmin) newYmin = ylo;
  }

  newSpace = 0.1*(newYmax - newYmin);
  hf->GetYaxis()->SetRangeUser(newYmin-newSpace, newYmax+newSpace);
  if (opt.Contains("zoom"))
    hf->GetYaxis()->SetRangeUser(-0.0008, 0.0085);

  leg->Draw();
  c->Update();
 
  // Draw text
  // if (opt.Contains("zoom"))
  //   ltx.SetTextSize(0.03); // 0.04
  // else

  // ltx.SetTextSize(0.04); // 0.05

  double xstart = (opt.Contains("zoom")) ? 0.7 : 0.6;
  ltx.DrawLatex(xstart, legbottom-1*0.07, Form("%.2g < p_{T} < %.2g GeV/c", t1, t2));
  // ltx.DrawLatex(xstart, legbottom-2*0.07, Form("%.2g < p_{T}^{assoc} < %.2g GeV/c", a1, a2));

  return c; 
}


TCanvas* DrawVnFromGlobalFit(int n, int ipt1, int ipt2, int ncb, 
			     int centbins[], TString opt)
{
  // Draw vn for one n on one canvas - many centralities
  initialize();
  double maxPtFit = ipt2==999 ? pthigh[maxpt-1] : pthigh[ipt2];
  const char* cname = Form("v%d_global%s_ptfit%.2gto%.2g",
			   n, opt.Data(), ptlow[ipt1], maxPtFit);
  const char* title = Form("v%d_global%s_ptfit%.2gto%.2g",
			   n, opt.Data(), ptlow[ipt1], maxPtFit);
  double lvleft = opt.Contains("multi")? 0.74 : 0.82;
  double lvbot = opt.Contains("multi")? 0.25 : 0.56;
  int cvHeight = opt.Contains("multi")? 500 : 500;
  TCanvas* cv = new TCanvas(cname, title, 700, cvHeight);
  TLegend* lv = new TLegend(lvleft, lvbot, 0.99, 0.97);
  lv->SetFillColor(kNone); lv->SetName("lv");
  lv->SetMargin(0.3);
  lv->SetBorderSize(1);
  lv->SetHeader("#splitline{Fit range}{(GeV/c):}");
  utils.padsetup(cv, 1, 1, "", 0.12, 0.01, 0.03, 0.2); // ######
  //  gPad->SetRightMargin(0.15);


  TPad* inset = 0;

  TObjArray* graphs = new TObjArray();
  double xmin = 0., xmax = 14.8; // 8.3;
  double ymin = -0.06, ymax = 0.26;
  if (opt.Contains("multi")) {
    ymin = -1.19;
    ymax = 0.52;
  }
  TH1F* hf = gPad->DrawFrame(xmin, ymin, xmax, ymax);
  hf->SetTitle(Form(";p_{T}^{t} (GeV/c);v_{%d}{GF}",n));
  int nXticks =  210;
  int nYticks = 210;
  TLine l;
  l.SetLineStyle(kDashed);
  if (n==1) {
    cv->cd();
    inset = new TPad("inset", "inset", 0.18, 0.18, 0.45, 0.59, kNone, 1, 0);
    inset->Draw();
    inset->cd();
    TH1F* hfi = gPad->DrawFrame(-0.05, -0.04, 4.2, 0.16);
    TAxis* hx = hfi->GetXaxis();     
    TAxis* hy = hfi->GetYaxis();
    hx->SetLabelSize(0.1);
    hy->SetLabelSize(0.1);
    hx->SetNdivisions(104);
    hy->SetNdivisions(104);
    gPad->SetBottomMargin(0.2);
    gPad->SetLeftMargin(0.2);
    gPad->SetRightMargin(0.01);
    gPad->SetTopMargin(0.01);
    //    utils.make_nice_axes(inset, hfi, 3, nXticks, nYticks, 1.4,-1.4);
    l.DrawLine(-0.05, 0, 4.2, 0);
    cv->cd();
  }

  utils.make_nice_axes(cv, hf, 1.3, nXticks, nYticks, 1.4,-1.4);
  hf->SetLabelOffset(0.0, "X");
  hf->SetTitleOffset(1.5, "Y");
  hf->GetXaxis()->SetTicks("+");
  //  hf->GetYaxis()->SetTicks("+-");

  TGraph* v1h[5][3];  
  if (opt.Contains("luzum")) {
    int lSty[] = {7, kSolid, kSolid}; 
    int lCol[] = {kBlack, kCyan, kBlue, kGreen+1, kRed}; 
    TString pidString[] = {"pi", "K", "p"};
    for (int m=0; m<5; m++) {
      for (int pid=0; pid<3; pid++) {
	v1h[m][pid] = Luzumv1(m, pidString[pid].Data()); 	
	TGraph* g1 = v1h[m][pid];
	g1->SetLineStyle(lSty[pid]);
	g1->SetLineColor(lCol[m]);
	g1->SetLineWidth(2);
      }
    }
    hf->GetXaxis()->SetRangeUser(0, 4.8);
    l.DrawLine(0, 0, 4.8, 0);
    v1h[0][0]->Draw("plsame");
    v1h[0][2]->Draw("plsame");
    v1h[4][0]->Draw("plsame");
    v1h[4][2]->Draw("plsame");
  }
  else
    l.DrawLine(hf->GetXaxis()->GetXmin(), 0, hf->GetXaxis()->GetXmax(),0);

  for (int cb=0; cb<ncb; cb++) {
    int k = centbins[cb];

    if (!opt.Contains("multi")) {    
      // last 2 args are included pt bins
      TGraphErrors* gtmp = vnGFvsPt(n, k, ipt1, ipt2);

      // clone to modify without ripples
      TGraphErrors* gc = (TGraphErrors*) gtmp->Clone();
      if (!gc) {
	Error("DrawVnFromGlobalFit()", "Problem making graph" );
	continue;
      }
      TGraphErrors* gtmp_sys = vnGFvsPt(n, k, ipt1, ipt2, "RIDGE", cfDef, "sys");
      TGraphErrors* gs = (TGraphErrors*) gtmp_sys->Clone();

      if (opt.Contains("luzum")) {
	if (centlow[k]==20)
	  utils.set_tgraph_props(gc, kBlue, kBlue, kFullCircle, 1.2);
	if (centlow[k]==40)
	  utils.set_tgraph_props(gc, kRed, kRed, kFullCircle, 1.2);
      }

      utils.draw_errorbox(gs, gs->GetFillColor());
      gc->Draw("epsame");
      graphs->Add(gc);

      TGraphErrors* open = (TGraphErrors*) gc->Clone(Form("%s_open", gc->GetName()));
      utils.set_tgraph_props(open, kBlack, kBlack, kOpenCircle, 1.2);
      open->Draw("epsame");

      lv->AddEntry(gc,Form("%d-%d%%",centlow[k],centhigh[k]),"l,p");
    }
    // NEW
    if (opt.Contains("multi")) {
      int b1,b2,b3,b4;
      b1=PtBin(0.25,0.5);
      b2=PtBin(0.75,1.0);
      b3=PtBin(2.0, 2.5);
      b4=PtBin(3.0, 4.0);
      TGraphErrors* gg1 = vnGFvsPt(n, k, b1, b2, "RIDGE", cfDef, "");
      TGraphErrors* gg2 = vnGFvsPt(n, k, b3, b4, "RIDGE", cfDef, "keepsign");
      TGraphErrors* gs1 = vnGFvsPt(n, k, b1, b2, "RIDGE", cfDef, "sys");
      TGraphErrors* gs2 = vnGFvsPt(n, k, b3, b4, "RIDGE", cfDef, "sys_keepsign");
      utils.set_tgraph_props(gg1, gg1->GetLineColor(), gg1->GetMarkerColor(), kOpenCircle, 1.0);
      utils.set_tgraph_props(gg2, gg2->GetLineColor(), gg2->GetMarkerColor(), kFullSquare, 1.0);
      
      utils.draw_errorbox(gs1, gs1->GetFillColor());
      utils.draw_errorbox(gs2, gs2->GetFillColor());
      gg1->Draw("epsame");
      gg2->Draw("epsame");
      graphs->Add(gs1);
      graphs->Add(gs2);
      const char* s1 = Form("#splitline{%.2g < p_{T}^{a} < %.2g  }{%d-%d%%}",
			    ptlow[b1], pthigh[b2], centlow[k],centhigh[k]);
      const char* s2 = Form("#splitline{%.2g < p_{T}^{a} < %.2g  }{%d-%d%%}",
			    ptlow[b3], pthigh[b4], centlow[k],centhigh[k]);
      lv->AddEntry(gg1, s1, "l,p");
      lv->AddEntry(gg2, s2, "l,p");

      
      if (n==1) {
	inset->cd();
	gg1->Draw("epsame");
	gg2->Draw("epsame");
	cv->cd();
      }
      
    }

  } // centrality bin loop

  if (opt.Contains("luzum")) {
    lv->AddEntry(v1h[0][0], Form("v.h. (#pi)"),"l,p");
    lv->AddEntry(v1h[0][2], Form("v.h. (p)"), "l,p");
  }
  lv->Draw();
  ltx.SetTextSize(0.05);

  if (opt.Contains("multi")) {
    gPad->SetTopMargin(0.02);
    if (ispp)
    ltx.DrawLatex(0.19, 0.92, Form("%s", "#splitline{pp}{2.76 TeV}"));
    else
    ltx.DrawLatex(0.19, 0.92, Form("%s", "#splitline{Pb-Pb}{2.76 TeV}"));
    ltx.DrawLatex(0.49, 0.94, Form("%.2g < |#Delta#eta| < 1.8", minRidgeDEta) ); 

  }
  else {
    ltx.DrawLatex(0.17, 0.88, Form("%.2g < |#Delta#eta| < 1.8", minRidgeDEta) ); 
    ltx.DrawLatex(0.17, 0.94, Form("%.2g < fit p_{T}^{a} < %.2g GeV/c",
				   ptlow[ipt1], maxPtFit));
    ltx.DrawLatex(0.19, 0.78, Form("%s", "#splitline{Pb-Pb}{2.76 TeV}"));
  }
  double newYmin = 0, newYmax = 0, newSpace = 0;
  for (int m=0; m<graphs->GetEntries(); m++) {
    TGraphErrors* tg = (TGraphErrors*)graphs->At(m);
    double yhi = TMath::MaxElement(tg->GetN()-1, tg->GetY());
    double ylo = TMath::MinElement(tg->GetN()-1, tg->GetY());
    if (yhi > newYmax) newYmax = yhi;
    if (ylo < newYmin) newYmin = ylo;
  }

  if (opt.Contains("multi")) {
    if(n==1) 
      newYmin *= 2.2;
    if(n==2) 
      newYmin = 0.02;
  }

  newSpace = 0.2*(newYmax - newYmin);
  //  if (opt.Contains("multi")) newSpace /= 2;
  if (n==5) newSpace *= 2;
  hf->GetYaxis()->SetRangeUser(newYmin-newSpace, newYmax+newSpace);
  cv->Modified();
  cv->Update();
  cv->Draw();
  return cv;
}

TCanvas* Drawv1to5(int ncb, int centbins[], int ipt1, int ipt2, TString opt)
{
  if (!opt.IsNull())
    cout<<opt.Data()<<endl;

  initialize();
  int cent1 = centlow[0];
  int cent2 = centhigh[centbins[ncb-1]];
  double pmax = 15.;
  if (ipt2 < 999.)
    pmax = pthigh[ipt2];
  const char* cname = Form("vn_etamin%02d_cent%dto%d_fitptbin%dto%d", 
			   (int)(10*minRidgeDEta), cent1, cent2, ipt1, ipt2);
  const char* title = Form("v1to5gf_etamin%02d_cent%dto%d_fitpt%.2gto%.2g", 
			   (int)(10*minRidgeDEta), cent1, cent2, ptlow[ipt1], pmax);
  TCanvas* cv = new TCanvas(cname, title, 1400, 500);
  TLegend* lv = new TLegend(0.09, 0.72, 0.27, 0.88);
  lv->SetFillColor(kNone);
  lv->SetBorderSize(0);
  TLine zero;
  double lv2bottom = 0.62;
  TLegend* lv2 = new TLegend(0.82, lv2bottom, 0.93, 0.95, "Centrality");
  lv2->SetFillColor(kNone);
  lv2->SetBorderSize(0);
  utils.padsetup(cv, 5, 1, "", 0.12, 0.01, 0.03, 0.2);

  for (int n=1; n<=5; n++) {
    Printf("Drawv1to5() n = %d / 5", n);
    cv->cd(n);
    double xmin = -0.43, xmax =  8.3;
    double ymin = -0.09, ymax =  0.26;
    if (ptlow[ipt1] > 4.) {
      xmin = 4.9;
      ymin = -0.1;
    }
    if (pthigh[ipt2] > 6.) {
      xmax = 11.8;
      ymax = 1;
    }
    if (opt.Contains("split")) {
      xmin = -0.43;
      xmax = 11.8;
      ymin = -0.2;
      ymax = 1;
    }

    TH1F* hf = gPad->DrawFrame(xmin, ymin, xmax, ymax);
    int nXticks = 206;
    int nYticks = 208;
    if (n==1) {
      utils.make_nice_axes(cv, hf, 2., nXticks, nYticks, -0.02, 0.02);
      hf->SetLabelOffset(-0.01, "X");
    }
    else {
      utils.make_nice_axes(cv, hf, 3.4, nXticks, nYticks, -0.075, 0.0);
      hf->SetLabelOffset(-0.057, "X");
      hf->SetTickLength(0.02, "X");
      hf->SetTickLength(0.05, "Y");
    }
    hf->GetXaxis()->SetTicks("+-");
    hf->GetYaxis()->SetTicks("+-");
    zero.SetLineStyle(kDashed);
    zero.DrawLine(xmin, 0., xmax, 0.);

    // Loop thru every cent bin. Skip non-requested ones.
    // for (int k=0; k<maxcent; k++) {
    //   bool centOK = 1;
    //   if (ncb != 999) {
    // 	centOK = 0;
    // 	for (int cb=0; cb<ncb; cb++) 
    // 	  if(k==centbins[cb])
    // 	    centOK = 1;
    //   }
    //   if (!centOK)
    // 	continue;

    for (int cb=0; cb<ncb; cb++) {
      int k = centbins[cb];
      TGraphErrors* gc = vnGFvsPt(n, k, ipt1, ipt2);
      TGraphErrors* gs = vnGFvsPt(n, k, ipt1, ipt2, "RIDGE", cfDef, "sys");

      if (!gc) {
	Error("Drawv1to5()", "Problem finding graph");
	continue;
      }

      TGraphErrors* gc_hipt=0;
      TGraphErrors* gs_hipt=0;
      if (opt.Contains("highptfit") || opt.Contains("split")) {
	gc_hipt = vnGFvsPt(n, k, ipt2+1, 999, "highptfit");
	gs_hipt = vnGFvsPt(n, k, ipt2+1, 999, "RIDGE", cfDef, "sys");
	gc_hipt->SetMarkerStyle(kOpenCircle);
      }

      int pale = gs->GetFillColor(); // it is CentColor(centlow[k], centhigh[k]).
      utils.draw_errorbox(gs, pale);
      gc->Draw("epsame");

      if (opt.Contains("split")) {
	utils.draw_errorbox(gs_hipt, pale);
	gc_hipt->Draw("epsame");
      }
      
      TGraphErrors* open = (TGraphErrors*) gc->Clone(Form("%s_open", gc->GetName()));
      utils.set_tgraph_props(open, kBlack, kBlack, kOpenCircle, 1.2);
      open->Draw("epsame");

      if (n==1)
	lv2->AddEntry(gc,Form("%d-%d%%",centlow[k],centhigh[k]),"l,p");

    } // cent bin k
  } // Fourier moment n
  
  cv->cd();
  TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1);
  overlay->SetFillStyle(4000); // transparent
  overlay->Draw();
  overlay->cd();
  
  ltx.SetTextSize(0.07);
  lv2->Draw();
  
  for (int n=1; n<=5; n++) {
    cv->cd(n);
    ltx.SetTextSize(n==1? 0.1 : 0.15);
    ltx.DrawLatex((n==1)? 0.85 : 0.7, 0.86, Form("v_{%d}", n));

    if (0 && n==1)
      AddPrelimStampToCurrentPad(0.45,0.7,0.75,0.95, "");
  }
  overlay->cd();
  ltx.SetTextSize(0.075);
  ltx.DrawLatex(0.5, 0.05, "p_{T} [GeV/c]");
  ltx.SetTextAngle(90);
  ltx.DrawLatex(0.05, 0.5, "v_{n }(p_{T})");
  ltx.SetTextAngle(0);
  ltx.SetTextSize(0.05);
  ltx.DrawLatex(0.14, 0.88, Form("%s", "#splitline{Pb-Pb}{2.76 TeV}"));
  ltx.SetTextSize(0.05);
  ltx.DrawLatex(0.13, 0.24, Form("%.2g < |#Delta#eta| < 1.8", minRidgeDEta) ); 
  
  const char* frange = Form("%.2g < fit p_{T}^{a} < %.2g GeV",
  			    ptlow[ipt1], pmax);  
  const char* frange1 = Form("#splitline{solid:}{%.2g < fit p_{T}^{a} < %.2g GeV}",
			    ptlow[ipt1], pthigh[ipt2]);  
  const char* frange2 = Form("#splitline{open:}{%.2g < fit p_{T}^{a} < %.2g GeV}",
			    ptlow[ipt2+1], pthigh[maxpt-1]);  

  if (opt.Contains("split")) {
    ltx.SetTextSize(0.04);
    ltx.DrawLatex(0.31, 0.8, frange1);
    ltx.DrawLatex(0.31, 0.7, frange2);
  }
  else
    ltx.DrawLatex(0.3, 0.27, frange);

  return cv;
}

TCanvas* Drawv2to5(int ncb, int centbins[], int ipt1, int ipt2, TString opt)
{
  if (!opt.IsNull())
    cout<<opt.Data()<<endl;

  initialize();
  int cent1 = centlow[0];
  int cent2 = centhigh[centbins[ncb-1]];
  double pmax = 15.;
  if (ipt2 < 999.)
    pmax = pthigh[ipt2];
  const char* cname = Form("v2to5_etamin%02d_cent%dto%d_fitptbin%dto%d", 
			   (int)(10*minRidgeDEta), cent1, cent2, ipt1, ipt2);
  const char* title = Form("v2to5gf_etamin%02d_cent%dto%d_fitpt%.2gto%.2g", 
			   (int)(10*minRidgeDEta), cent1, cent2, ptlow[ipt1], pmax);
  TCanvas* cv = new TCanvas(cname, title, 1400, 500);
  TLine zero;
  double lv2bottom = 0.58; // 0.5;
  TLegend* lv2 = new TLegend(0.77, lv2bottom, 0.9, 0.95, "Centrality");
  lv2->SetFillColor(kNone);
  lv2->SetBorderSize(0);
  utils.padsetup(cv, 4, 1, "", 0.08, 0.01, 0.03, 0.2);

  for (int n=2; n<=5; n++) {
    Printf("Drawv2to5() n = %d", n);
    cv->cd(n-1);
    double xmin = -0.2;
    double xmax =  7.6;
    TH1F* hf = gPad->DrawFrame(xmin, -0.12, xmax, 0.26); //was y = -0.03, 0.26
    int nXticks = 206;
    int nYticks = 208;
    if (n==2) {
      utils.make_nice_axes(cv, hf, 2., nXticks, nYticks, -0.02, 0.02);
      hf->SetLabelOffset(0.0, "X"); //       hf->SetLabelOffset(-0.01, "X");
    }
    else {
      utils.make_nice_axes(cv, hf, 2.6, nXticks, nYticks, -0.075, 0.0);
      hf->SetLabelOffset(-0.02, "X");
      hf->SetTickLength(0.02, "X");
      hf->SetTickLength(0.05, "Y");
    }
    hf->GetXaxis()->SetTicks("+-");
    hf->GetYaxis()->SetTicks("+-");
    zero.SetLineStyle(kDashed);
    zero.DrawLine(xmin, 0., xmax, 0.);

    for (int cb=0; cb<ncb; cb++) {
      int k = centbins[cb];
      
      TGraphErrors* gc = vnGFvsPt(n, k, ipt1, ipt2);
      TGraphErrors* gs = vnGFvsPt(n, k, ipt1, ipt2, "RIDGE", cfDef, "sys");

      if (!gc) {
	Error("Drawv2to5()", "Problem finding graph");
	continue;
      }

      int pale = gs->GetFillColor(); // it is CentColor(centlow[k], centhigh[k]).
      utils.draw_errorbox(gs, pale);
      //      gs->Draw("e2psame");
      gc->Draw("epsame");
      TGraphErrors* open = (TGraphErrors*) gc->Clone(Form("%s_open", gc->GetName()));
      utils.set_tgraph_props(open, kBlack, kBlack, kOpenCircle, 1.2);
      open->Draw("epsame");

      if (n==2)
	lv2->AddEntry(gc,Form("%d-%d%%",centlow[k],centhigh[k]),"l,p");

    } // cent bin k
  } // Fourier moment n
  
  cv->cd();
  TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1);
  overlay->SetFillStyle(4000); // transparent
  overlay->Draw();
  overlay->cd();
  
  ltx.SetTextSize(0.07);
  lv2->Draw();
  
  for (int n=2; n<=5; n++) {
    cv->cd(n-1);
    ltx.SetTextSize(n==2? 0.1 : 0.15);
    ltx.DrawLatex((n==2)? 0.85 : 0.7, 0.86, Form("v_{%d}", n));
  }
  overlay->cd();
  ltx.SetTextSize(0.075);
  ltx.DrawLatex(0.48, 0.04, "p_{T} [GeV/c]");
  ltx.SetTextAngle(90);
  ltx.DrawLatex(0.02, 0.5, "v_{n }{GF}");
  ltx.SetTextAngle(0);
  ltx.SetTextSize(0.05);
  ltx.DrawLatex(0.33, 0.88, Form("%s", "#splitline{Pb-Pb}{2.76 TeV}"));
  ltx.SetTextSize(0.05);
  ltx.DrawLatex(0.33, 0.75, Form("%.2g < |#Delta#eta| < 1.8", minRidgeDEta) ); 
  
  const char* frange = Form("#splitline{Global fit range}{%.2g < p_{T} < %.2g GeV}", 
			    ptlow[ipt1], pmax);  
  ltx.DrawLatex(0.56, 0.75, frange);
  
  return cv;
}

TCanvas* DrawChi2vsPtTrig(int npta, int ptabins[], int ncb, int centbins[], TString opt)
{
    
  // TODO: options-- global or 1to5, chi2/N or prob
  double a1=ptlow[ptabins[0]], a2=pthigh[ptabins[npta-1]];
  int c1=centlow[centbins[0]], c2=centhigh[centbins[ncb-1]];
  TString cname = Form("chi2_etamin%02d_pta%.2gto%.2g_cent%dto%d",
		       (int)(10*minRidgeDEta), a1, a2, c1, c2 );
  
  if (!opt.IsNull() ) {
    cname.Append("_");
    cname.Append(opt.Data());
  }

  TCanvas* cv = new TCanvas(cname, cname, 1);
  TLegend* lv = new TLegend(0.09, 0.72, 0.27, 0.88); // w.r.t. cv
  lv->SetFillColor(kNone);
  lv->SetBorderSize(0);

  double lv2bottom = 0.62;
  TLegend* lv2 = new TLegend(0.82, lv2bottom, 0.93, 0.95, "Centrality"); // w.r.t. cv
  lv2->SetFillColor(kNone);
  lv2->SetBorderSize(0);
  utils.padsetup(cv, 1, 1, "", 0.12, 0.01, 0.03, 0.2);
  
  double xmin = -0.43;
  double xmax =  12.5;
  TH1F* hf = gPad->DrawFrame(xmin, -0.09, xmax, 26);
  // int nXticks = 210;
  // int nYticks = 210;
  hf->GetXaxis()->SetTicks("+-");
  hf->GetYaxis()->SetTicks("+-");
    
// IN PROGRESS---------------
    for (int j=0; j<npta; j++) {
      for(int k=0; k<ncb; k++) {

	TGraphErrors *gc = CorrFnChi2VsTrigPt(ptabins[j], centbins[k], 1, 5, "ndf");
	if (!gc) {
	  Error("DrawChi2()", "!gc");
	  continue;
	}    
	int col = CentColor(centlow[centbins[k]], centhigh[centbins[k]]);
	utils.set_tgraph_props(gc, col, col, kFullCircle, 1.6);
	gc->Draw("epsame");
	lv2->AddEntry(gc, centLabel(centbins[k]), "l,p");
      }
    }
  
  // cv->cd();
  // TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1);
  // overlay->SetFillStyle(4000); // transparent
  // overlay->Draw();
  // overlay->cd();
  
  // ltx.SetTextSize(0.07);
  // lv2->Draw();
  // overlay->cd();
  // ltx.SetTextSize(0.075);
  // ltx.DrawLatex(0.5, 0.05, "trigger p_{T} [GeV/c]");
  // ltx.SetTextAngle(90);
  // ltx.DrawLatex(0.05, 0.5, "#chi^{2} / NDF");
  // ltx.SetTextAngle(0);

  return cv;
}


TCanvas* DrawAgreement(int icent, TString opt)
{
  // Draw Chi2 on one canvas for each
  // centrality bin k.

  initialize();
  int cent1 = centlow[icent];
  int cent2 = centhigh[icent];

  if (icent==999) {
    cent1 = 0;
    cent2 = 90;
  }
  int k = icent;
  TString cname = Form("agmt_etamin%02d_cent%dto%d_%s", 
		       (int)(10*minRidgeDEta), cent1, cent2, opt.Data());
  if (!opt.IsNull()) {
    cname.Append("_");
    cname.Append(opt.Data());
  }

  TCanvas* cv = new TCanvas(cname.Data(), cname.Data(), 1400, 500);
  TLegend* lv = new TLegend(0.09, 0.72, 0.27, 0.88);
  lv->SetFillColor(kNone);
  lv->SetBorderSize(0);

  double lv2bottom = 0.62;
  TLegend* lv2 = new TLegend(0.82, lv2bottom, 0.93, 0.95, "Centrality");
  lv2->SetFillColor(kNone);
  lv2->SetBorderSize(0);
  cv->Divide(5,1);

  for (int n=1; n<=5; n++) {
    cv->cd(n);

      TH2F *gc = Agreement2DHist(k,n);
      if (!gc) {
	Error("DrawAgreement()", "Problem finding graph" );
	continue;
      }

      //      gc->GetZaxis()->SetRangeUser(1.e-1, 1000);

      gPad->SetLogx(); gPad->SetLogy(); gPad->SetLogz();

      if (opt.Contains("lego")) {
	gc->SetTitleOffset(0.8, "X");
	gc->GetXaxis()->SetLabelSize(0.07);
	gc->GetYaxis()->SetLabelSize(0.07);
	gc->GetZaxis()->SetLabelSize(0.07);
	gc->GetXaxis()->SetTitleSize(0.1);
	gc->GetYaxis()->SetTitleSize(0.1);
	gc->GetZaxis()->SetTitleSize(0.1);
	gc->GetXaxis()->SetNdivisions(8);
	gc->GetYaxis()->SetNdivisions(8);
	
	gPad->SetTopMargin(0.2);
	gPad->SetLeftMargin(0.12);
	gPad->SetRightMargin(0.01);
	gc->Draw("lego2");
      }
      else if (opt.Contains("contour"))
	gc->Draw("colz");

      if (n==1) {
	lv2->AddEntry(gc, Form("%d-%d%%",centlow[k],centhigh[k]),"l,p");
      }
      
  } // Fourier moment n
  
  cv->cd();
  TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1);
  overlay->SetFillStyle(4000); // transparent
  overlay->Draw();
  overlay->cd();
  
  ltx.SetTextSize(0.07);
  const char* str = Form("#left|V_{n#Delta} / (v_{n}^{t} v_{n}^{a})#right|");
  ltx.DrawLatex(0.42, 0.92, str);
  ltx.DrawLatex(0.91, 0.92, Form("%s", centLabel(k)));
  
  for (int n=1; n<=5; n++) {
    cv->cd(n);
    ltx.SetTextSize(0.10);
    ltx.DrawLatex(0.05, 0.76, Form("n = %d", n));
  }
  return cv;
}

void AddPrelimStampToCurrentPad(float x1, float y1, float x2, float y2, TString opt)
{
  TASImage* logo = new TASImage("../commonfigs/alice-prelim.png"); // 287x323
  TPad* logopad = new TPad("logopad", "logo", x1,y1,x2,y2);
  if (opt.Contains("border")) {
    logopad->SetFillColor(kRed);
    logopad->SetBorderSize(2);
  }
  logopad->Draw();  
  logopad->cd(); 
  logo->Draw("");
  return;
}

int CentBin(double cntLo, double cntHi) 
{
  int bin = -123;
  for (int k=0; k<maxcent; k++) {
    if (centlow[k]==cntLo && centhigh[k]==cntHi)
      bin = k;
  }
  if (bin < 0 && 0)
    Warning("CentBin","No bin for %.2g - %.2g", cntLo, cntHi);

  return bin;
}

int PtBin(double ptLo, double ptHi) 
{
  int bin = -123;
  for (int i=0; i<maxpt; i++) {
    if (ptlow[i]==ptLo && pthigh[i]==ptHi)
      bin = i;
  }
  if (bin < 0)
    Warning("PtBin","No bin for %.2g - %.2g", ptLo, ptHi);

  return bin;
}

double MomCorrection(int i, int j, int k)
{
  if(0)cout<<i<<j<<k<<endl;
  // For momentum (non)conservation in v1.
  // Multiplicitave weight: w = w_t*w_a.
  // w_t = ptt - <pt^2>/<pt>
  // w_a = pta - <pt^2>/<pt>

  //  double c = 1./3 * -0.00185;

  //  (ptmean[i]-)* ptmean[j]

  double correction =  1./3 * -0.00185 * ptmean[i]* ptmean[j];
  return correction;
}

void MakeVnDVsQGraphs(int n1, int n2, int k, const char* region, const char* corrtype, TString opt)
{
  // The first two errType constants have points set at the central measured value of
  // VnDelta, with statistical and systematic error bars
  // respectively. The last two have points set at the upper and lower
  // systematic value, with statistical error bars set to be the same
  // as for meas_stat.
  const int nErrTypes = 4;
  enum errType {MEAS_STAT, MEAS_SYS, HI, LO};
  //  const char* errTypeStr[] = {"meas_stat","meas_sys","hi","lo"};
  
  // Fill a huge array with all the coefficients, then use it to set the tgraphs.
  const int nN = gNMax + 1; // 1-10, don't use 0.
  const int nI = maxpt+1;
  const int nJ = maxpt+1;
  const int nT = nErrTypes;
  double vVal[nN][nI][nJ][nT];
  double vErr[nN][nI][nJ][nT];
  double vMix[nN][nI][nJ][nT]; // Just the acceptance correction unc.

  double vSin[nN][nI][nJ][nT]; // <sin n delta phi> for sys. estimation
  double vSinErr[nN][nI][nJ][nT];

  for (int i=0; i<maxpt; i++) {
    for (int j=0; j<=i; j++) {

      double esin  = 0;  // residual sine error
      double ebw   = 0;  // bin width error \propto n / (# dphi bins)
      double ecnt  = 0;  // centrality determination
      double emix  = 0;  // acc. correction error
      double etrk  = 0;  // tracking resolution
      double ev1   = 0;  // v1
      double esys  = 0;  // quadrature sum

      TH1* hst=0, *hsys=0; // stat and sys error, but same points

      hst  = Hist(region, corrtype, i, j, k, "");
      hsys = Hist(region, corrtype, i, j, k, "sys"); 

      if (!hst) Error("MakeVnDVsQGraphs()","!hst");
      if (!hsys) Error("MakeVnDVsQGraphs()","!hsys");

      esin = SineError(i,j,k); // resid_sin / 2;

      for (int n=n1; n<=n2; n++) { // Order of coeff.
	
	// VnDelta on measured points, with statistical errors.
	double val=0, err=0, resid_sin=0, rs_err=0;
	VnDelta(n, *hst,  val, err, "");
	VnDelta(n, *hst,  resid_sin, rs_err, "sin");
	
	// VnDelta on high and low sys values. Stat err. not used.
	double val_hi=0, xxx=0;
	VnDelta(n, *hsys, val_hi, xxx, "hi");
	double val_lo=0, yyy=0;
	VnDelta(n, *hsys, val_lo, yyy, "lo");

	// Momentum conservation correction
	if (n==1 && opt.Contains("ptcons")) { // eqs. 3-7, PRL 106, 102301 (2011)
	  double factor =  ispp ? 0.212 : 0; // 7e-5; <------ trial & error: 7e-5 must be close for 0-20%
	  val    += factor*ptmean[i]*ptmean[j];
	  val_hi += factor*ptmean[i]*ptmean[j];
	  val_lo += factor*ptmean[i]*ptmean[j];
	}

	int nbins = hst->GetNbinsX();
	double mean = (val_hi + val_lo + val) / 3.;
	double ms = 
	  (val_hi-mean)*(val_hi-mean) + 
	  (val_lo-mean)*(val_lo-mean) + 
	  (val   -mean)*(val   -mean);
	double rms = TMath::Sqrt(ms/2); // sqrt(1/(N-1) sum (x_i - <x>)^2 )


	ebw = n/nbins/TMath::Sqrt(12)*val; // This is about 0.008 * n * val
	etrk = 0.01*ptmean[j]*val; // see comments in comparisons/VnDcomparison.C 
	ecnt = 0.01*val;
	emix = rms/2; // take half not to double-count stat err. // 0.25*TMath::Abs(val_hi - val_lo);
	ev1 = (n==1) ? 0.001*ptmean[i]*ptmean[j]*val : 0.;

	esys = TMath::Sqrt(ecnt*ecnt + emix*emix + ebw*ebw + etrk*etrk + ev1*ev1 + esin*esin);
	
	double pointVal = val;
	double pointErr = err;

	vVal[n][i][j][MEAS_STAT] = pointVal;
	vErr[n][i][j][MEAS_STAT] = pointErr;
	vVal[n][i][j][MEAS_SYS] = pointVal;
	vErr[n][i][j][MEAS_SYS] = esys;
	vVal[n][i][j][HI] = pointVal + esys;
	vErr[n][i][j][HI] = pointErr;
	vVal[n][i][j][LO] = pointVal - esys;
	vErr[n][i][j][LO] = pointErr;
	vSin[n][i][j][MEAS_STAT] = resid_sin;
	vSinErr[n][i][j][MEAS_STAT] = rs_err;
	vMix[n][i][j][MEAS_SYS] = emix;

      } // n
    } // j
  } // i

  // Now set the graphs from the arrays.
  for (int n=n1; n<=n2; n++) {
    for (int t=0; t<nErrTypes; t++) {
  
      TString gname = Form("V%dDeltaVsGlobalIndex_cent%dto%d",
			   n,centlow[k],centhigh[k]);
      if (t==HI)
	gname.Append("hi");
      else if (t==LO)
	gname.Append("lo");
      else if (t==MEAS_SYS)
	gname.Append("meas_sys");
      else
	gname.Append("meas_stat");

      if (opt.Contains("sine"))
	gname.Append("_sine");

      if (opt.Contains("ALL"))
	gname.Append("_ALL");
      
      TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data());
      if (!g) {
	g = new TGraphErrors();
	g->SetName(gname.Data());
      }

      TGraphErrors* gmix = 0;
      if (t==MEAS_SYS) {
	gmix = new TGraphErrors();
	gname.Append("_mix");
	gmix->SetName(gname.Data());
      }
      
      double xer = 0.4; // was 0.5
      for (int i=0; i<maxpt; i++) {
	for (int j=0; j<=i; j++) {
	  int q = GlobalIndex(i, j);
	  g->SetPoint(q, q+0.5, vVal[n][i][j][t]);
	  g->SetPointError(q, xer, vErr[n][i][j][t]);

	  if (opt.Contains("sine")) {
	    g->SetPoint(q, q+0.5, vSin[n][i][j][MEAS_STAT]);
	    g->SetPointError(q, xer, vSinErr[n][i][j][MEAS_STAT]);
	  }

	  if (t==MEAS_SYS) {
	    gmix->SetPoint(q, q+0.5, vVal[n][i][j][t]);
	    gmix->SetPointError(q, xer, vMix[n][i][j][t]);
	  }

	}
      }
      
      gList->Add(g);
      if (gmix) gList->Add(gmix);
    }
  }
  
  return;
}

TGraphErrors* VnDVsQ(int n, int k, const char* region, const char* corrtype, TString opt)
{

  // Graph of pair fourier coefficients VnDelta vs. the linear
  // global index q. Arguments:
  // 
  // n, k:   fourier index, centrality bin
  // region: "RIDGE", "NSJET" or "ALL"
  // corrtype:  "s", "m", "cA", "cB"
  //
  // Select opt "hi_sys" or "lo_sys" for hi and low pts
  // with central stat errors and "sys" for the central points
  // with sys. errors. Just like Hist().

  TString gname = Form("V%dDeltaVsGlobalIndex_cent%dto%d",
		       n,centlow[k],centhigh[k]);

  if (opt.Contains("hi"))
    gname.Append("hi");
  else if (opt.Contains("lo"))
    gname.Append("lo");
  else if (opt.Contains("sys"))
    gname.Append("meas_sys");
  else
    gname.Append("meas_stat");

  if (opt.Contains("sine"))
    gname.Append("_sine");

  if (opt.Contains("ALL"))
    gname.Append("_ALL");
  
  TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data()); 
  if (g) {
    if (0)
      Info("VnDVsQ()", "Returning preexisting graph");
    return g;
  }

  TString opt1 = (opt.Contains("sine")) ? "sine" : "";
  if (opt.Contains("ALL"))
    opt1.Append("_ALL");
  if (opt.Contains("ptcons"))
    opt1.Append("_ptcons");

  MakeVnDVsQGraphs(1, gNMax, k, region, corrtype, opt1);
  g = (TGraphErrors*) gList->FindObject(gname.Data()); 
  if (!g) {
    Error("VnDVsQ()", "No graph %s", gname.Data());
    gSystem->Exit(-1);
  }
  return g;
}

double MeanPt(int i, int j, int k, TString t_or_a, TString opt)
{
  // For combined centralities, map k's to the closest available bin
  // double cb1[] =  {0, 0, 0,2,2, 1,3,0,    0, 1, 2, 3, 4,  5, 10, 20, 30, 40, 60 };
  // double cb2[] =  {10,20,2,5,10,3,5,5,    1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 90 };
// OBJ: TProfile2D	hMeanPtTrgCen1	cen=1 (0.5) : 0 at: 0x102af06e0
//  OBJ: TProfile2D	hMeanPtAssCen1	cen=1 (0.5) : 0 at: 0x108308e60
//  OBJ: TProfile2D	hMean2PtTrgCen1	cen=1 (0.5) : 0 at: 0x108309480
//  OBJ: TProfile2D	hMean2PtAssCen1	cen=1 (0.5) : 0 at: 0x108309ad0
//  OBJ: TProfile2D	hMeanPtTrgCen2	cen=2 (1.5) : 0 at: 0x10830a120
//  OBJ: TProfile2D	hMeanPtAssCen2	cen=2 (1.5) : 0 at: 0x10830a770
//  OBJ: TProfile2D	hMean2PtTrgCen2	cen=2 (1.5) : 0 at: 0x10830adc0
//  OBJ: TProfile2D	hMean2PtAssCen2	cen=2 (1.5) : 0 at: 0x10830b410
//  OBJ: TProfile2D	hMeanPtTrgCen3	cen=3 (2.5) : 0 at: 0x10830ba60
//  OBJ: TProfile2D	hMeanPtAssCen3	cen=3 (2.5) : 0 at: 0x10830c0b0
//  OBJ: TProfile2D	hMean2PtTrgCen3	cen=3 (2.5) : 0 at: 0x10830c700
//  OBJ: TProfile2D	hMean2PtAssCen3	cen=3 (2.5) : 0 at: 0x10830cd50
//  OBJ: TProfile2D	hMeanPtTrgCen4	cen=4 (3.5) : 0 at: 0x10830d3a0
//  OBJ: TProfile2D	hMeanPtAssCen4	cen=4 (3.5) : 0 at: 0x10830d9f0
//  OBJ: TProfile2D	hMean2PtTrgCen4	cen=4 (3.5) : 0 at: 0x10830e040
//  OBJ: TProfile2D	hMean2PtAssCen4	cen=4 (3.5) : 0 at: 0x10830e690
//  OBJ: TProfile2D	hMeanPtTrgCen5	cen=5 (4.5) : 0 at: 0x10830ece0
//  OBJ: TProfile2D	hMeanPtAssCen5	cen=5 (4.5) : 0 at: 0x10830f330
//  OBJ: TProfile2D	hMean2PtTrgCen5	cen=5 (4.5) : 0 at: 0x10830f980
//  OBJ: TProfile2D	hMean2PtAssCen5	cen=5 (4.5) : 0 at: 0x10830ffd0
//  OBJ: TProfile2D	hMeanPtTrgCen6	cen=6 (7.5) : 0 at: 0x108310620
//  OBJ: TProfile2D	hMeanPtAssCen6	cen=6 (7.5) : 0 at: 0x108310c70
//  OBJ: TProfile2D	hMean2PtTrgCen6	cen=6 (7.5) : 0 at: 0x1083112c0
//  OBJ: TProfile2D	hMean2PtAssCen6	cen=6 (7.5) : 0 at: 0x108311910
//  OBJ: TProfile2D	hMeanPtTrgCen7	cen=7 (15.0) : 0 at: 0x108311f60
//  OBJ: TProfile2D	hMeanPtAssCen7	cen=7 (15.0) : 0 at: 0x1083125b0
//  OBJ: TProfile2D	hMean2PtTrgCen7	cen=7 (15.0) : 0 at: 0x108312c00
//  OBJ: TProfile2D	hMean2PtAssCen7	cen=7 (15.0) : 0 at: 0x108313250
//  OBJ: TProfile2D	hMeanPtTrgCen8	cen=8 (25.0) : 0 at: 0x1083138a0
//  OBJ: TProfile2D	hMeanPtAssCen8	cen=8 (25.0) : 0 at: 0x108313ef0
//  OBJ: TProfile2D	hMean2PtTrgCen8	cen=8 (25.0) : 0 at: 0x108314540
//  OBJ: TProfile2D	hMean2PtAssCen8	cen=8 (25.0) : 0 at: 0x108314b90
//  OBJ: TProfile2D	hMeanPtTrgCen9	cen=9 (35.0) : 0 at: 0x1083151e0
//  OBJ: TProfile2D	hMeanPtAssCen9	cen=9 (35.0) : 0 at: 0x108315830
//  OBJ: TProfile2D	hMean2PtTrgCen9	cen=9 (35.0) : 0 at: 0x108315e80
//  OBJ: TProfile2D	hMean2PtAssCen9	cen=9 (35.0) : 0 at: 0x1083164d0
//  OBJ: TProfile2D	hMeanPtTrgCen10	cen=10 (45.0) : 0 at: 0x108316b20
//  OBJ: TProfile2D	hMeanPtAssCen10	cen=10 (45.0) : 0 at: 0x108317170
//  OBJ: TProfile2D	hMean2PtTrgCen10	cen=10 (45.0) : 0 at: 0x1083177c0
//  OBJ: TProfile2D	hMean2PtAssCen10	cen=10 (45.0) : 0 at: 0x108317e10
//  OBJ: TProfile2D	hMeanPtTrgCen11	cen=11 (55.0) : 0 at: 0x108318460
//  OBJ: TProfile2D	hMeanPtAssCen11	cen=11 (55.0) : 0 at: 0x108318ab0
//  OBJ: TProfile2D	hMean2PtTrgCen11	cen=11 (55.0) : 0 at: 0x108319100
//  OBJ: TProfile2D	hMean2PtAssCen11	cen=11 (55.0) : 0 at: 0x108319750
//  OBJ: TProfile2D	hMeanPtTrgCen12	cen=12 (75.0) : 0 at: 0x108319da0
//  OBJ: TProfile2D	hMeanPtAssCen12	cen=12 (75.0) : 0 at: 0x10831a3f0
//  OBJ: TProfile2D	hMean2PtTrgCen12	cen=12 (75.0) : 0 at: 0x10831aa40
//  OBJ: TProfile2D	hMean2PtAssCen12	cen=12 (75.0) : 0 at: 0x10831b090

  int k2 = -1;
  if (k == CentBin(0,  1))   k2 = 1;  
  else if (k == CentBin(1,  2))   k2 = 2;  
  else if (k == CentBin(2,  3))   k2 = 3;  
  else if (k == CentBin(3,  4))   k2 = 4;  
  else if (k == CentBin(4,  5))   k2 = 5;  
  else if (k == CentBin(5, 10))   k2 = 6;  
  else if (k == CentBin(10,20))   k2 = 7;  
  else if (k == CentBin(20,30))   k2 = 8;  
  else if (k == CentBin(30,40))   k2 = 9;  
  else if (k == CentBin(40,50))   k2 = 10; 
  else if (k == CentBin(50,60))   k2 = 11; 
  else if (k == CentBin(60,90))   k2 = 12; 
  else if (k == CentBin(0, 10))   k2 = 5;    //CentBin(4, 5);   
  else if (k == CentBin(0, 20))   k2 = 6;    //CentBin(5, 10);  
  else if (k == CentBin(0,  2))   k2 = 1;    //CentBin(0, 1);   
  else if (k == CentBin(2,  5))   k2 = 4;    //CentBin(3, 4);   
  else if (k == CentBin(2, 10))   k2 = 6;    //CentBin(5, 10);  
  else if (k == CentBin(1,  3))   k2 = 2;    //CentBin(1, 2);   
  else if (k == CentBin(3,  5))   k2 = 4;    //CentBin(3, 4);   
  else if (k == CentBin(0,  5))   k2 = 3;    //CentBin(2, 3);   

  else if (k == CentBin(0, 0))   k2 = 1;    //pp
  else {
    Error("MeanPt()", "cent bin %d (%d-%d%%) not assigned", k, centlow[k], centhigh[k]);
  }

  const char* m = opt.Contains("squared") ? "Mean2" : "Mean";

  TString prName_t = Form("h%sPtTrgCen%d", m, k2);
  TString prName_a = Form("h%sPtAssCen%d", m, k2);
  TProfile2D* hmpt_t = (TProfile2D*)fin->Get(prName_t.Data());
  TProfile2D* hmpt_a = (TProfile2D*)fin->Get(prName_a.Data());

  if (!hmpt_t) {
    Error("MeanPt", "Problem finding TProfile %s. Input %d", prName_t.Data(), k);
    gSystem->Exit(-1);
  }
  //  cout << hmpt_t->GetTitle() << " " << centlow[k] << " " << centhigh[k] << endl;

  double mptt = hmpt_t->GetBinContent(i+1, j+1);
  double mpta = hmpt_a->GetBinContent(i+1, j+1);

  if (0) {
    // x axis is trig pt,  y axis assc pt
    TAxis* ax = hmpt_t->GetXaxis();
    TAxis* ay = hmpt_t->GetYaxis();
    Printf("trig %.3g - %.3g, assc: %.3g - %.3g, <pt_t,a> %.3g, %.3g",
	   ax->GetBinLowEdge(i+1), 
	   ax->GetBinUpEdge(i+1), 
	   ay->GetBinLowEdge(j+1), 
	   ay->GetBinUpEdge(j+1),
	   mptt, mpta );
  }

  if (t_or_a.Contains("t"))
    return mptt;
  if (t_or_a.Contains("a"))
    return mpta;
  else
    return -1;
}


int CentColor(double cen1, double cen2) 
{
  if (cen1==0 && cen2==10)
    return kGray+3;
  if (cen1==0 && cen2==20)
    return kMagenta+3;
  if (cen1==0 && cen2==2)
    return kBlack;
  if (cen1==2 && cen2==10)
    return kRed;
  if (cen1==10 && cen2==20)
    return kOrange-3;
  if (cen1==20 && cen2==30)
    return kGreen+1;
  if (cen1==30 && cen2==40)
    return kYellow+2;
  if (cen1==40 && cen2==50)
    return kAzure;
  if (cen1==50 && cen2==60)
    return kBlue-1;
  if (cen1==60 && cen2==90)
    return kViolet;

  return kBlack;
}

int CentColorPale(double cen1, double cen2) 
{
  if (cen1==0 && cen2==10)
    return kGray+1;
  if (cen1==0 && cen2==20)
    return kPink-1; // kPink-6;
  if (cen1==0 && cen2==2)
    return kGray+1;
  if (cen1==2 && cen2==10)
    return kRed-9;
  if (cen1==10 && cen2==20)
    return kOrange-4;
  if (cen1==20 && cen2==30)
    return kGreen-7;
  if (cen1==30 && cen2==40)
    return kYellow-8;
  if (cen1==40 && cen2==50)
    return kAzure-4;
  if (cen1==50 && cen2==60)
    return kBlue-9;
  if (cen1==60 && cen2==90)
    return kViolet-9;

  return kBlack;
}

int PtColor(double p1, double p2)
{
 if (p1==0.25 && p2==0.5)
    return kGray+1;
 if (p1==0.5 && p2==0.75)
    return kYellow+3;
 if (p1==0.75 && p2==1.0)
   return kPink-6;
 if (p1==1.0 && p2==1.5)
   return kBlack;
 if (p1==1.5 && p2==2.0)
   return kRed;
 if (p1==2.0 && p2==2.5)
   return kOrange-3;
 if (p1==2.5 && p2==3)
   return kGreen+2;
 if (p1==3 && p2==4)
   return kBlue;
 if (p1==4 && p2==5)
   return kViolet + 2;
 if (p1==5 && p2==6)
   return kViolet+10;
 if (p1==4 && p2==6)
   return kViolet + 2;
 if (p1==6 && p2==8)
   return kOrange-7;
 if (p1==8 && p2==15)
   return kOrange+9;

 return kBlack; 
}

int PtColorPale(double p1, double p2)
{
 if (p1==0.25 && p2==0.5)
    return kGray;
 if (p1==0.5 && p2==0.75)
    return kYellow+1;
 if (p1==0.75 && p2==1.0)
   return kPink-4;
 if (p1==1.0 && p2==1.5)
   return kGray+1;
 if (p1==1.5 && p2==2.0)
   return kRed-7;
 if (p1==2.0 && p2==2.5)
   return kOrange-4;
 if (p1==2.5 && p2==3)
   return kSpring-4;
 if (p1==3 && p2==4)
   return kAzure+6;
 if (p1==4 && p2==5)
   return kViolet+1;
 if (p1==5 && p2==6)
   return kViolet+8;
 if (p1==4 && p2==6)
   return kViolet+1;
 if (p1==6 && p2==8)
   return kOrange-9;
 if (p1==8 && p2==15)
   return kOrange+6;

 return kGray; 
}

int PtaColor(double p1, double p2)
{
 if (p1==0.25 && p2==0.5)
    return kGray+1;
 if (p1==0.5 && p2==0.75)
   return kBlack;
 if (p1==0.75 && p2==1.0)
   return kRed+1;
 if (p1==1.0 && p2==1.5)
   return kOrange-3;//kBlack;
 if (p1==1.5 && p2==2.0)
   return kGreen+2;//kRed;
 if (p1==2.0 && p2==2.5)
   return kAzure-3;//kOrange-3;
 if (p1==2.5 && p2==3)
   return kViolet + 2;//kBlue;
 if (p1==3 && p2==4)
   return kViolet + 10;//kBlue;
 if (p1==4 && p2==5)
   return kMagenta+3; //kViolet + 2;
 if (p1==5 && p2==6)
   return kViolet;
 if (p1==4 && p2==6)
   return kMagenta-5;
 if (p1==6 && p2==8)
   return kGray+2;
 if (p1==8 && p2==15)
   return kGray+3;

 return kBlack; 
}

int PtaColorPale(double p1, double p2)
{
 if (p1==0.25 && p2==0.5)
    return kGray;
 if (p1==0.5 && p2==0.75)
   return kGray+2;
 if (p1==0.75 && p2==1.0)
   return kPink-2;
 if (p1==1.0 && p2==1.5)
   return kOrange-4;
 if (p1==1.5 && p2==2.0)
   return kGreen-3;
 if (p1==2.0 && p2==2.5)
   return kAzure-4;
 if (p1==2.5 && p2==3)
   return kViolet;
 if (p1==3 && p2==4)
   return kViolet + 2;
 if (p1==4 && p2==5)
   return kMagenta+3;
 if (p1==5 && p2==6)
   return kViolet+10;
 if (p1==4 && p2==6)
   return kMagenta-5;
 if (p1==6 && p2==8)
   return kGray+2;
 if (p1==8 && p2==15)
   return kGray+3;

 return kGray; 
}

void SaveCanvases(TObjArray* canvases, const char* fileName)
{
  TFile* f = new TFile(fileName, "recreate");

  for (int n=0; n<canvases->GetEntries(); n++) {
    TCanvas* c = (TCanvas*)canvases->At(n);
    c->Write(c->GetTitle());
  }
  f->Close();
  return;
}

void SaveCanvasesFromFile(const char* rootFile, const char* targetDir, const char* tag, const char* fileType)
{
  // Get a list of canvases from rootFile into array, then save each
  // to its own file in targetDir/. fileType = "eps", "pdf", "C",
  // "png", etc. Not all formats have been tested.
  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  TString name = "";
  TString base(targetDir);
  TFile *cFile = new TFile(rootFile, "read");
  cout << cFile->GetName() << endl;
  TObjArray* cList = GetObjectsFromFile(*cFile, "TCanvas");

  for (int n=0; n<cList->GetEntries(); n++) {
    TCanvas* c = (TCanvas*)cList->At(n);
    if (c) {
      name = "";
      name = base;
      name += TString("/");
      name += TString(fileType);
      name += TString("/");
      name += TString(c->GetTitle());
      name += TString(".");
      name += TString(fileType);
      cout<<name.Data()<<endl;

      c->Draw();
      c->Modified();
      c->Update();
      c->SaveAs(name.Data());
    }
    else
      Error("SaveCanvasesFromFile()", "!c");
  }

  if (1) {
    utils.print_pdf(cList, Form("%s/all-figs%s", targetDir, tag), "pdf");
  }
  
  return;
}

//void SaveCanvases(TObjArray* canvases, const char* fileName)
//void SaveCanvasesFromFile(const char* rootFile, const char* targetDir, const char* fileType)
TObjArray* GetObjectsFromFile(TFile& file, TString clname, TString dir)
{
  file.cd(dir.Data());

  TObjArray* objList = new TObjArray();
  TIter next(gDirectory->GetListOfKeys());
  TKey *key;
  
  while ((key=(TKey*)next())) {
    TString className(key->GetClassName());
    TString keyName(key->GetName());
    if (1) 
      printf("%10s %20s\n", className.Data(), keyName.Data());
    
    if (className.Contains(clname)) {
      objList->Add(gDirectory->Get(keyName.Data()));
    }
  }

  cout << objList->GetEntries() << " objects retrieved from "
       << file.GetName()  << "/" << gDirectory->GetName() 
       << endl;

  return objList;
}