analysis code for harmonic decomposition paper
authoraadare <aadare@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 2 Feb 2012 21:07:24 +0000 (21:07 +0000)
committeraadare <aadare@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 2 Feb 2012 21:07:24 +0000 (21:07 +0000)
PWGCF/Correlations/DPhi/FourierDecomposition/lrc/CreatePaperFigs.C [new file with mode: 0644]
PWGCF/Correlations/DPhi/FourierDecomposition/lrc/FourierPlus.C [new file with mode: 0644]
PWGCF/Correlations/DPhi/FourierDecomposition/lrc/ProtonProtonAna.C [new file with mode: 0644]
PWGCF/Correlations/DPhi/FourierDecomposition/lrc/RunFourierPlus.C [new file with mode: 0644]

diff --git a/PWGCF/Correlations/DPhi/FourierDecomposition/lrc/CreatePaperFigs.C b/PWGCF/Correlations/DPhi/FourierDecomposition/lrc/CreatePaperFigs.C
new file mode 100644 (file)
index 0000000..891a0c7
--- /dev/null
@@ -0,0 +1,467 @@
+// --- Switches and settings ---
+bool saveOutput = 1;
+bool resave = 0;  // Sometimes pdf gets corrupted...retry then quit.
+const char* tag = "-Dec2"; // identifier appended to filenames
+const char* outputDir = "../analysis-output";
+const char* inputHistFile = "$MYDATA/gsi07_etamin08-withmeanpt.root";
+TString grFile   = Form("%s/root-objects/graphs%s.root",
+                       outputDir, tag);
+TString canvFile = Form("%s/root-objects/canvases%s.root",
+                       outputDir, tag);
+TString pdfFile  = Form("%s/plots/all-figs%s",
+                       outputDir, tag);
+
+// --- Globals ---
+TFile* pwg2File   = new TFile("comparisons/PWG2points.root");
+TFile* alexv2File = new TFile("comparisons/v2_Eta04_pass2.root");
+TFile* alexv3File = new TFile("comparisons/v3_EP_eta04.root");
+TLatex ltx;
+TObjArray* cList = new TObjArray();
+
+// --- Functions ---
+void PlotCollection(int k, int i, int j, TString opt);
+void SaveCanvases(TObjArray* canvases, const char* fileName);
+TObjArray* GetObjectsFromFile(TFile& file, TString clname, TString dir = "");
+void SaveCanvasesFromFile(const char* rootFile, const char* targetDir, const char* fileType);
+TGraphErrors* PWG2v(int n, double cen1, double cen2, double etaGap=1.);
+
+void CreatePaperFigs()
+{
+  ltx.SetNDC();
+  int NCENT = 8;
+  bool isBatch = gROOT->IsBatch();
+  gROOT->LoadMacro("./common/Utils.C+");
+  PlotUtils utils;
+  gROOT->LoadMacro("FourierPlus.C+g");
+  initialize(inputHistFile);
+
+  if (resave) {
+    SaveCanvasesFromFile(canvFile.Data(), Form("%s/plots", outputDir), "pdf");
+    return;
+  }
+
+  int ncb_set = 6;
+  int centbins_set[] = { CentBin(0,2), 
+                        CentBin(2,10), 
+                        CentBin(10,20), 
+                        CentBin(20,30), 
+                        CentBin(30,40), 
+                        CentBin(40,50) };
+  
+  // SingleDraw2D(CentBin(0,10), PtBin(2, 2.5), PtBin(1.5, 2));
+  // SingleDraw2D(CentBin(0,10), PtBin(3, 4), PtBin(2, 2.5));
+  // //  cList->Add(SingleDraw2D(CentBin(0,2), PtBin(2, 2.5), PtBin(1.5, 2)));
+  // SingleDraw(CentBin(0,2), PtBin(6, 8), PtBin(1., 1.5), "color");
+  // SingleDraw(CentBin(0,10), PtBin(6, 8), PtBin(1., 1.5), "color");
+
+    // cList->Add(Drawv2to5(ncb_set, centbins_set, 0, 999, "") );
+  //  DrawAgreement(CentBin(0, 10), "contour");
+  //DrawChi2(CentBin(0, 10), "");
+
+  int t1[] = {PtBin(2., 2.5)};
+  int a1[] = {PtBin(1.5, 2.)};
+  int t2[] = {PtBin(8,15)};
+  int a2[] = {PtBin(6, 8)};
+  int c2[] = {CentBin(40,50), CentBin(0,20)};
+  int c3[] = {CentBin(0,2)};
+  
+  int ncb = 5;
+  int centbins[] = { CentBin(40,50), 
+                    CentBin(20,30), 
+                    CentBin(10,20), 
+                    CentBin(2,10), 
+                    CentBin(0,2) };
+
+  // int cbin05[] = { CentBin(0,5) };
+  // int cbin010[] = { CentBin(0,10) };
+  
+  // cList->Add(DrawVnDeltaVsN(1,  t1, 1, a1, 1,cbin05, "05_nobars_ext"));
+  // cList->Add(DrawVnDeltaVsN(1,  t1, 1, a1, 1,cbin010, "010_nobars_ext"));
+  
+  // cList->Add(DrawVnDeltaVsN(1,  t1, 1, a1, ncb, centbins, "nobars_ext")); // all centralities
+  // cList->Add(DrawVnDeltaVsN(1,  t1, 1, a1, 1,   c3, "cent02_nobars_ext")); // just 0-2%
+
+  // cList->Add(DrawVnDeltaVsN(1,  t1, 1, a1, ncb, centbins,
+  //                       "sine_ext")); // for sys errors.
+  // cList->Add(DrawVnDeltaVsN(1,  t2, 1, a2, 2, c2, "sine_ext"));
+  // cList->Add(DrawVnDeltaVsN(1, t2, 1, a2, 2, c2, "nobars_ext"));
+  
+  // DrawQ(CentBin(0,10), 2);
+  // DrawQ(CentBin(0,10), 3);
+  //  return;
+
+  // For CERN courier article
+  cList->Add(SingleDrawPlain(CentBin(0,2), PtBin(2, 2.5), PtBin(1.5, 2), 
+                            "harmonics_sum") );
+
+  int aa1[1] = {PtBin(0.5, 0.75)};
+  int aa2[3] = {PtBin(0.5, 0.75), PtBin(1.5, 2), PtBin(2.5, 3)};
+  int cc1[1] = {CentBin(0,20)};
+  int cc2[5] = {CentBin(0,2), CentBin(2,10), CentBin(10,20),
+  CentBin(20,30), CentBin(30,40) };
+
+  //  cList->Add( DrawChi2vsPtTrig(1, aa2, 1, cc1, "") );
+
+  if (1) { // Draw Correlation functions
+
+    if (0) {// Draw 3 CFs overlaid together
+      int ptt = PtBin(2., 2.5);
+      int pta = PtBin(1.5, 2.);
+      TLegend* ml = new TLegend(0.5, 0.58, 0.98, 0.8, "Pb-Pb 2.76 TeV");
+      ml->SetFillColor(kNone);
+      ml->SetBorderSize(0);
+      int col[] = {18, kAzure-9, kBlack, kYellow, kOrange, kRed };
+      cList->Add(SingleDrawPlain(CentBin(3,5), ptt, pta, "ALL"));
+      TH1* h0 = (TH1*)(gPad->GetListOfPrimitives())->FindObject("h");
+      utils.set_hist_props(h0, kNone, kNone, kNone, kOpenCircle, 1.4);
+      TH1* hUltCnt[3];
+      int cntBins[] = {CentBin(3,5), CentBin(2,3), CentBin(0,1)};
+      for (int n=0; n<3; n++) {
+       hUltCnt[n] = Hist("ALL", "cA", ptt, pta, cntBins[n], Form("_%d", n));
+       utils.set_hist_props(hUltCnt[n], kBlack, col[n], col[n], kFullCircle, 1.4);
+       hUltCnt[n]->SetLineWidth(2); 
+       hUltCnt[n]->Draw("same");
+       TH1* hc = (TH1*)hUltCnt[n]->Clone();
+       utils.set_hist_props(hc, kBlack, kBlack, kBlack, kOpenCircle, 1.4);
+       hc->Draw("same");
+      }
+      for (int n=2; n>=0; n--) {
+       ml->AddEntry(hUltCnt[n], centLabel(cntBins[n]), "epl");
+      }
+      ml->Draw();
+    }
+    if (0) { // "plain" CFs (no lower ratio panel)
+      int centEvol[] = { CentBin(0,5), 
+                        CentBin(0,1), 
+                        CentBin(1,2), 
+                        CentBin(2,3), 
+                        CentBin(3,4), 
+                        CentBin(4,5)};
+      for (int k=0;k<6;k++) {
+       int i = PtBin(2., 2.5), j = PtBin(1.5, 2);
+       int kk = centEvol[k]; 
+       
+       cList->Add(SingleDrawPlain(kk,i,j, "ALL"));
+       TH1* h0 = (TH1*)(gPad->GetListOfPrimitives())->FindObject("h");
+       h0->GetYaxis()->SetRangeUser(0.98, 1.03);
+       
+       cList->Add(SingleDrawPlain(kk,i,j, "harmonics_ALL"));
+       TH1* h0 = (TH1*)(gPad->GetListOfPrimitives())->FindObject("h");
+       h0->GetYaxis()->SetRangeUser(0.98, 1.03);
+      }
+    }
+
+    // New: High ptt, low pta...
+    PlotCollection(CentBin(0,2), PtBin(6, 8), PtBin(1., 1.5), "2D");
+    PlotCollection(CentBin(0,10), PtBin(6, 8), PtBin(1., 1.5), "2D");
+  
+    PlotCollection(CentBin(0,1),  PtBin(2., 2.5), PtBin(1.5,2.), "");
+    PlotCollection(CentBin(0,2),  PtBin(2., 2.5), PtBin(1.5,2.), "");
+    PlotCollection(CentBin(0,10), PtBin(2, 2.5),  PtBin(1.5, 2), "2D");
+    PlotCollection(CentBin(0,10), PtBin(3.0, 4.), PtBin(2.,2.5), "2D");
+    PlotCollection(CentBin(0,20), PtBin(8., 15.), PtBin(6.,8), "2D");
+    PlotCollection(CentBin(30,40),PtBin(6., 8.),  PtBin(1,1.5), "");
+    cList->Add(SingleDraw2D(CentBin(0,20),PtBin(8, 15), PtBin(6, 8), "zoom"));
+    ltx.DrawLatex(0.4, 0.7, "#splitline{zoomed to }{0 < C(#Delta#phi) < 5}");
+  }
+  //  return;
+
+  if (1) { // v_n{GF} for v1 to v5
+    cList->Add(Drawv1to5(ncb_set, centbins_set, 0, PtBin(2.0, 2.5), "fitbelow2.5") );
+    cList->Add(Drawv1to5(ncb_set, centbins_set, 0, 999, "") );
+    cList->Add(Drawv2to5(ncb_set, centbins_set, 0, 999, "") );
+    cout << "Drawv1to5() ok" << endl;
+  }
+
+  if (1) { // Global fit plots
+    DrawQ(CentBin(0,10), 6); // bug somewhere...w/o this, the next canvas is messed up
+    for (int n=1; n<=5; n++) {
+      cList->Add( DrawQ(CentBin(0,10), n) );
+    }
+    for (int cb=0; cb<ncb_set; cb++) {
+      for (int n=1; n<=5; n++) {
+       int k = centbins_set[cb];
+       cout << Form("DrawQ() %s, n = %d / 5", centLabel(k), n) << endl;
+       cList->Add( DrawQ(k, n) );
+      }
+    }
+  }
+  
+  if (1) { // "Power spectrum" plots
+    int t1[] = {PtBin(2., 2.5)};
+    int t2[] = {PtBin(8,15)};
+    int t3[] = {PtBin(6,8)};
+    int a1[] = {PtBin(1.5, 2.)};
+    int a2[] = {PtBin(6, 8)};
+    int a3[] = {PtBin(1., 1.5)};
+
+    int c2[] = {CentBin(40,50), CentBin(0,20)};
+    int c3[] = {CentBin(0,2)};
+
+    int ncb = 5;
+    int centbins[] = { CentBin(40,50), 
+                      CentBin(20,30), 
+                      CentBin(10,20), 
+                      CentBin(2,10), 
+                      CentBin(0,2) };
+
+    int ncb2 = 4;
+    int centbins2[] = { CentBin(40,50), 
+                      CentBin(20,30), 
+                      CentBin(10,20), 
+                      CentBin(0,10)};
+
+    
+    // VnDelta vs n...
+    if (0) {
+    // Fig 2,3 with color bars
+    cList->Add(DrawVnDeltaVsN(1,t1, 1, a1, ncb, centbins, "")); // all
+    cList->Add(DrawVnDeltaVsN(1,t1, 1, a1, 1,   c3, "cent02")); // just 0-2%
+    cList->Add(DrawVnDeltaVsN(1,t2, 1, a2, 2, c2, "ext"));
+    }
+    // And without
+    cList->Add(DrawVnDeltaVsN(1,t1, 1, a1, ncb, centbins, "nobars"));
+    cList->Add(DrawVnDeltaVsN(1,t1, 1, a1, 1,   c3, "cent02_nobars"));
+    cList->Add(DrawVnDeltaVsN(1,t2, 1, a2, 2, c2, "ext_nobars"));
+
+    // up to n=12
+    cList->Add(DrawVnDeltaVsN(1,t1, 1, a1, ncb, centbins,"sine_ext")); // sys err.
+    cList->Add(DrawVnDeltaVsN(1,t2, 1, a2, 2, c2, "sine_ext"));
+
+
+    // High ptt x low pta
+    cList->Add(DrawVnDeltaVsN(1,a2, 1, a1, ncb2, centbins2, ""));
+    cList->Add(DrawVnDeltaVsN(1,t3, 1, a3, ncb2, centbins2, "nobars"));
+    cList->Add(DrawVnDeltaVsN(1,t2, 1, a3, ncb2, centbins2, "nobars"));
+
+    // Additional
+    cList->Add(DrawVnDeltaVsN(1,t1, 1, a1, ncb2, centbins2, "sine"));
+
+
+    // vn{GF} vs n...
+    cList->Add(DrawGlobalvnVsN(1, t1, ncb, centbins, ""));
+    cList->Add(DrawGlobalvnVsN(1, a1, ncb, centbins, ""));
+    //    cList->Add(DrawVnDeltaVsN(1, t1, 1, a1, ncb_set, centbins_set, "zoom"));
+
+  }
+
+  if (1) { // VnDelta vs. trigger pt -- 1 canvas/centrality bin
+    int cbins[] = { CentBin(0,10) , CentBin(10,20), CentBin(20,30), 
+                   CentBin(30,40), CentBin(40,50) };
+    int ptabins02[] = {PtBin(0.25, 0.5), PtBin(1.0, 1.5)};
+    int ptabins10[] = {PtBin(0.25, 0.5), PtBin(1.0, 1.5), PtBin(2.0, 2.5)};
+
+    // cList->Add(DrawVnVsPtTrig(CentBin(0,2), PtBin(0.25, 0.5),
+    // PtBin(1.0, 1.5)));
+    cList->Add(DrawVnVsPtTrig(CentBin(0,2), 2, ptabins02, ""));
+
+    for (int i=0; i<5; i++) {
+
+      // cList->Add(DrawVnVsPtTrig(cbins[i], PtBin(0.25, 0.5), PtBin(2.5,
+      // 3.0)));
+      cList->Add(DrawVnVsPtTrig(cbins[i], 3, ptabins10, ""));
+    }
+  }
+
+ if (1) { // individual v_n{GF} canvases
+    // v1 comparison with Luzum
+    int ncbv1 = 2;
+    int centbinsv1[] = {CentBin(0,10), CentBin(40,50)};
+
+    //    Just our points
+    for (int n=1; n<=5; n++) {
+      cList->Add(DrawVnFromGlobalFit(n,0, 999, ncb_set,centbins_set,""));
+      cList->Add(DrawVnFromGlobalFit(n,0, PtBin(2.0, 2.5), ncb_set,centbins_set,""));
+    }
+    int ncb = 2;
+    int centbins[] = {CentBin(0,2), CentBin(30, 40)};
+    for (int n=2; n<=5; n++) {
+      TCanvas* c = DrawVnFromGlobalFit(n,0,999,ncb,centbins,"_pwg2");
+      TGraphErrors* g0002 = PWG2v(n, 0,  2 );
+      TGraphErrors* g3040 = PWG2v(n, 30, 40);
+      utils.set_tgraph_props(g0002, kBlack, kBlack, kOpenSquare, 1.5);
+      utils.set_tgraph_props(g3040, kRed,   kRed,   kFullSquare, 0.8);
+      g0002->Draw("epsame");
+      g3040->Draw("epsame");
+      if (0)
+       if(n==2 || n==3)
+         Alexv(n, 30, 40)->Draw("epsame"); // Great agmt, but I need <pt>
+      // ltx.SetTextSize(0.04);
+      // ltx.DrawLatex(0.17, 0.95, Form("Global fit"));
+      // ltx.DrawLatex(0.45, 0.95, Form("S.P. (CERN-PH-EP-2011-073)"));
+      // ltx.SetTextSize(0.05);
+      // ltx.DrawLatex(0.45, 0.88, Form("1.0 < |#Delta#eta| < 1.6", n));
+      TObject* obj = c->GetListOfPrimitives()->FindObject("lv");
+      TLegend* leg = obj;
+      leg->AddEntry(g0002, Form("SP 0-2%%", n), "epl");
+      leg->AddEntry(g3040, Form("SP 30-40%%", n), "epl");
+      cList->Add(c);
+    }
+  }
+
+  if (saveOutput) { 
+    SaveGraphs(grFile.Data());
+    SaveCanvases(cList, canvFile.Data()); 
+
+    if (1) // individual pdfs + multipage pdf
+      SaveCanvasesFromFile(canvFile.Data(), Form("%s/plots", outputDir), "pdf");
+    if (0) // individual plot macros
+      SaveCanvasesFromFile(canvFile.Data(), Form("%s/plots", outputDir), "C");
+
+  }
+  
+  return; // end CreatePaperFigs()
+}
+
+TGraphErrors* Alexv(int n, double cen1, double cen2)
+{
+  // 0-5% -> 0, 5-10% ->1, 10-20%->2
+  int bin = -1;
+  if (cen1==0  && cen2==5)   bin = 0;
+  else if (cen1==5  && cen2==10)  bin = 1;
+  else if (cen1==10  && cen2==20) bin = 2;
+  else if (cen1==20  && cen2==30) bin = 3;
+  else if (cen1==30  && cen2==40) bin = 4;
+  else if (cen1==40  && cen2==50) bin = 5;
+  else if (cen1==50  && cen2==60) bin = 6;
+  else if (cen1==60  && cen2==70) bin = 7;
+  else if (cen1==70  && cen2==80) bin = 8;
+
+  if (bin<0)
+    return 0;
+
+  TFile* f = 0;
+  if (n==2)
+    f = alexv2File;
+  else if (n==3)
+    f = alexv3File;
+  else
+    return 0;
+
+  const char* xtra = cen1==30 ? "ALICE" : "";
+  TGraphErrors* g = 0;
+  const char* name = "";
+
+  if (n==2)
+    name = Form("v2_all_Eta04_hist_%d", bin);
+  if (n==3)
+    name = Form("v3_eta04_%d", bin);
+
+   if (0) cout << name << endl;
+
+  g = (TGraphErrors*)f->Get(name);
+  if (!g)
+    Error("Alexv()", "%s not found");
+
+  return g;
+}
+
+TGraphErrors* PWG2v(int n, double cen1, double cen2, double etaGap)
+{
+  const char* xtra = cen1==30 ? "ALICE" : "";
+  TGraphErrors* g = 0;
+   const char* name = Form("GrSP_%02d%02d%s_v%d_etaGap%.0f", 
+                          (int)cen1, (int)cen2, xtra, n, 10*etaGap);
+
+   if (0) cout << name << endl;
+
+  g = (TGraphErrors*)pwg2File->Get(name);
+  if (!g)
+    Error("PWG2v()", "%s not found");
+  
+  return g;
+}
+
+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 PlotCollection(int k, int i, int j, TString opt)
+{
+  // cList->Add(SingleDrawPlain(k, i, j, "pl"));
+  cList->Add(SingleDrawPlain(k, i, j, "harmonics"));
+  cList->Add(SingleDraw(k, i, j, "color"));
+  //  cList->Add(SingleDraw(k, i, j, "gray"));
+  //cList->Add(SingleDraw(k, i, j, "gray_upto10"));
+  cList->Add(SingleDraw(k, i, j, "global"));
+  if (opt.Contains("2D")) {
+    cList->Add(SingleDraw2D(k, i, j, "pl"));
+  }
+  return;
+}
+
+void SaveCanvasesFromFile(const char* rootFile, const char* targetDir, 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", base.Data(), tag), "pdf");
+  }
+  
+  return;
+}
+
+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;
+}
diff --git a/PWGCF/Correlations/DPhi/FourierDecomposition/lrc/FourierPlus.C b/PWGCF/Correlations/DPhi/FourierDecomposition/lrc/FourierPlus.C
new file mode 100644 (file)
index 0000000..df7c2f9
--- /dev/null
@@ -0,0 +1,4168 @@
+
+#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;
+}
diff --git a/PWGCF/Correlations/DPhi/FourierDecomposition/lrc/ProtonProtonAna.C b/PWGCF/Correlations/DPhi/FourierDecomposition/lrc/ProtonProtonAna.C
new file mode 100644 (file)
index 0000000..8399951
--- /dev/null
@@ -0,0 +1,92 @@
+// --- Global Switches and settings ---
+bool saveOutput = 1;
+bool resave = 0;  // Sometimes pdf gets corrupted...retry then quit.
+const char* tag = "-pp-nov28-etamin12"; // identifier appended to filenames
+const char* outputDir = "../test";
+const char* inputHistFile = "$MYDATA/lhc11a_etamin12.root";
+TString grFile   = Form("%s/root-objects/graphs%s.root",
+                       outputDir, tag);
+TString canvFile = Form("%s/root-objects/canvases%s.root",
+                       outputDir, tag);
+TString pdfFile  = Form("%s/plots/all-figs%s",
+                       outputDir, tag);
+int NCENT = 1;
+TLatex ltx;
+TObjArray* cList = new TObjArray();
+bool isBatch = gROOT->IsBatch();
+
+void PlotCollection(int k, int i, int j, TString opt);
+
+int Setup()
+{
+  ltx.SetNDC();
+  gROOT->LoadMacro("~/Dropbox/ALICE/common/Utils.C+");
+  gROOT->LoadMacro("~/Dropbox/ALICE/common/IOUtilFns.C");
+  gROOT->LoadMacro("FourierPlus.C+g");
+  initialize(inputHistFile);
+  SetProtonProtonFlag(1);
+  SetMinRidgeDEta(1.2);
+  return 0;
+}
+
+void ProtonProtonAna()
+{
+  Setup();
+  cout << IsProtonProton() << endl;
+  if (resave) {
+    SaveCanvasesFromFile(canvFile.Data(), Form("%s/plots", outputDir), tag, "pdf");
+    return;
+  }
+
+  int nt = 5, na = 6;
+  int t1[] = {PtBin(1., 1.5), PtBin(2., 2.5), PtBin(4,5), PtBin(6,8), PtBin(8,15)};
+  int a1[] = {PtBin(0.25., 0.5), PtBin(0.5., 0.75), PtBin(1., 1.5), PtBin(2., 2.5), 
+             PtBin(3, 4), PtBin(6, 8)};
+
+  //    cList->Add( DrawQ(0, 1, Form("ptabin%dto%d_%s",PtBin(0.25,0.5),PtBin(4,5), "ptcons")));
+
+    if (1) {
+      for (int n=1; n<=5; n++) {
+       cList->Add( DrawQ(0, n, Form("ptabin%dto%d_%s",PtBin(0.25,0.5),PtBin(8,15), "RIDGE")));
+      }
+      for (int n=1; n<=5; n++) {
+       cList->Add( DrawQ(0, n, Form("ptabin%dto%d_%s",PtBin(0.25, 0.5),PtBin(1,1.5), "RIDGE")));
+      }
+    }
+    
+  int ppcb[] = {CentBin(0,0)};
+  cList->Add(Drawv1to5(1, ppcb, PtBin(0.25, 0.5), PtBin(1., 1.5), "ptcons" ));
+
+  if (1) {
+    for (int i=0; i<nt; i++) {
+      for (int j=0; j<na; j++) {
+       int ti = t1[i], aj = a1[j];
+       if (ti >= aj)
+         cList->Add(PlotCollection(ti, aj, "2D"));
+      }
+    }
+  }
+  
+
+  if (saveOutput) {
+    SaveGraphs(grFile.Data());
+    SaveCanvases(cList, canvFile.Data()); 
+
+    if (1) // individual pdfs + multipage pdf
+      SaveCanvasesFromFile(canvFile.Data(), Form("%s/plots", outputDir), tag, "pdf");
+    if (0) // individual plot macros
+      SaveCanvasesFromFile(canvFile.Data(), Form("%s/plots", outputDir), tag, "C");
+  }
+
+  return;
+}
+
+void PlotCollection(int i, int j, TString opt)
+{
+  //  cList->Add(SingleDraw(0, i, j, ""));
+  if (opt.Contains("2D")) {
+    cList->Add(SingleDraw2D(0, i, j, "pl"));
+  }
+  cList->Add(SingleDrawPlain(0, i, j, "harmonics_sum"));
+  return;
+}
diff --git a/PWGCF/Correlations/DPhi/FourierDecomposition/lrc/RunFourierPlus.C b/PWGCF/Correlations/DPhi/FourierDecomposition/lrc/RunFourierPlus.C
new file mode 100644 (file)
index 0000000..1969e25
--- /dev/null
@@ -0,0 +1,188 @@
+// --- Global Switches and settings ---
+bool saveOutput = 1;
+bool resave = 1;  // Sometimes pdf gets corrupted...retry then quit.
+const char* tag = "-v3.4"; // identifier appended to filenames
+const char* outputDir = "../test";
+const char* inputHistFile = "$MYDATA/gsi07_etamin08-withmeanpt.root";
+TString grFile   = Form("%s/root-objects/graphs%s.root",
+                       outputDir, tag);
+TString canvFile = Form("%s/root-objects/canvases%s.root",
+                       outputDir, tag);
+TString pdfFile  = Form("%s/plots/all-figs%s",
+                       outputDir, tag);
+int NCENT = 8;
+TLatex ltx;
+TObjArray* cList = new TObjArray();
+bool isBatch = gROOT->IsBatch();
+
+void PlotCollection(int k, int i, int j, TString opt);
+
+int Setup()
+{
+  ltx.SetNDC();
+  gROOT->LoadMacro("~/Dropbox/ALICE/common/Utils.C+");
+  gROOT->LoadMacro("~/Dropbox/ALICE/common/IOUtilFns.C");
+  gROOT->LoadMacro("FourierPlus.C+g");
+  initialize(inputHistFile);
+  return 0;
+}
+
+void RunFourierPlus()
+{
+  Setup();
+
+  if (resave) {
+    SaveCanvasesFromFile(canvFile.Data(), Form("%s/plots", outputDir), tag, "pdf");
+    return;
+  }
+
+  int t1[] = {PtBin(2., 2.5)};
+  int t2[] = {PtBin(8,15)};
+  int t3[] = {PtBin(6,8)};
+  int t4[] = {PtBin(4,5)};
+  int a1[] = {PtBin(1.5, 2.)};
+  int a2[] = {PtBin(6, 8)};
+  int a3[] = {PtBin(1., 1.5)};
+  int a4[] = {PtBin(3, 4)};
+  int a5[] = {PtBin(0.5., 0.75)};
+  int c2[] = {CentBin(40,50),CentBin(0,20)};
+  int c3[] = {CentBin(0,2)};
+  int c4[] = {CentBin(0,10)};
+  int c5[] = {CentBin(0,5)};
+
+  int ncb_set = 5;
+  int centbins_set[] = { CentBin(0,2), 
+                        CentBin(2,10), 
+                        CentBin(10,20), 
+                        CentBin(20,30), 
+                        CentBin(40,50) };
+  int ncb_1 = 6;
+  int centbins_1[] = { CentBin(0,2), 
+                      CentBin(2,10), 
+                      CentBin(10,20), 
+                      CentBin(20,30), 
+                      CentBin(30,40), 
+                      CentBin(40,50) };
+  int ncb_2 = 2;
+  int centbins_2[] = { CentBin(0,20), 
+                      CentBin(40,50) };
+  int ncb_3 = 2;
+  int centbins_3[] = { CentBin(0,10), 
+                      CentBin(40,50) };
+  int ncb_4 = 5;
+  int centbins_4[] = { CentBin(40,50), 
+                    CentBin(20,30), 
+                    CentBin(10,20), 
+                    CentBin(2,10), 
+                    CentBin(0,2) };
+  
+  // For CERN courier article
+  cList->Add(SingleDrawPlain(CentBin(0,2), PtBin(2, 2.5), PtBin(1.5, 2), 
+                            "harmonics_sum") );
+  
+  // Correlation functions
+  PlotCollection(CentBin(0,1),  PtBin(2., 2.5), PtBin(1.5,2.), "");
+  PlotCollection(CentBin(0,2),  PtBin(2., 2.5), PtBin(1.5,2.), "");
+  PlotCollection(CentBin(0,10), PtBin(2, 2.5),  PtBin(1.5, 2), "2D");
+  PlotCollection(CentBin(0,10), PtBin(3.0, 4.), PtBin(2.,2.5), "2D");
+  PlotCollection(CentBin(0,20), PtBin(8., 15.), PtBin(6.,8), "2D");
+  PlotCollection(CentBin(30,40),PtBin(6., 8.),  PtBin(1,1.5), "");
+  cList->Add(SingleDraw2D(CentBin(0,20),PtBin(8, 15), PtBin(6, 8), "zoom"));
+  ltx.DrawLatex(0.4, 0.7, "#splitline{zoomed to }{0 < C(#Delta#phi) < 5}");
+  
+  if (1) { // v_n{GF} for v1 to v5
+    cList->Add(Drawv1to5(ncb_1, centbins_1, 0, PtBin(2.0, 2.5), "fitbelow2.5") );
+    cList->Add(Drawv1to5(ncb_1, centbins_1, 0, 999, "") );
+    cList->Add(Drawv2to5(ncb_1, centbins_1, 0, 999, "") );
+    cout << "Drawv1to5() ok" << endl;
+  }
+  //  cList->Add(Drawv2to5(ncb_set, centbins_set, 0, 999, "") );  
+
+  if (1) { // v_n{GF} for v1 to v5 - superimposed with two fit ranges
+    cList->Add(Drawv1to5(ncb_2, centbins_2, PtBin(5., 6.), 999));
+    cList->Add(Drawv1to5(ncb_2, centbins_2, 0, PtBin(3., 4.), "split") );
+    cList->Add(Drawv1to5(ncb_2, centbins_2, 0, PtBin(4., 5.), "split") );
+    cout << "Drawv1to5() ok" << endl;
+  }
+
+  if (1) { // Global fit plots
+    DrawQ(CentBin(0,10), 6); // bug somewhere...w/o this, the next canvas is messed up
+    for (int n=1; n<=5; n++) {
+      cList->Add( DrawQ(CentBin(0,10), n) );
+    }
+    for (int cb=0; cb<ncb_set; cb++) {
+      for (int n=1; n<=5; n++) {
+       int k = centbins_1[cb];
+       cout << Form("DrawQ() %s, n = %d / 5", centLabel(k), n) << endl;
+       cList->Add( DrawQ(k, n) );
+      }
+    }
+  }
+  if (1) { // Global fit plots at high pt only
+    for (int n=1; n<=5; n++) {
+      cout << Form("DrawQ() %s, n = %d / 5", "highptfit", n) << endl;
+      cList->Add( DrawQ(CentBin(0,10), n, "highptfit") );   // pta > 5
+      cList->Add( DrawQ(CentBin(0,20), n, "highptfit") );   // pta > 5
+    }
+    cout << "DrawQ() ok" << endl;
+  }
+
+  cList->Add( DrawQ(CentBin(0,20), 1, 
+                   Form("ptabin%dto%d_%s",PtBin(0.25,0.5),PtBin(0.75,1.0), "RIDGE")));
+  cList->Add( DrawQ(CentBin(0,20), 1, 
+                   Form("ptabin%dto%d_%s",PtBin(2.,2.5),PtBin(3,4), "RIDGE")));
+  
+  // individual v_n{GF} canvases
+  cList->Add(DrawVnFromGlobalFit(1,0,999,ncb_3,centbins_3,"multi"));
+  cList->Add(DrawVnFromGlobalFit(2,0,999,ncb_3,centbins_3,"multi"));
+
+  // VnDelta vs n...
+  cList->Add(DrawVnDeltaVsN(1,t1, 1, a1, ncb_4, centbins_4, "nobars"));
+  cList->Add(DrawVnDeltaVsN(1,t1, 1, a1, 1,   c3, "cent02_nobars"));
+  cList->Add(DrawVnDeltaVsN(1,t2, 1, a2, 2, c2, "ext_nobars"));
+  cout << "DrawVnDeltaVsN() ok" << endl;
+  
+  // // vn{GF} vs n...
+  // cList->Add(DrawGlobalvnVsN(1, t1, ncb_4, centbins_4, ""));
+  // cList->Add(DrawGlobalvnVsN(1, a1, ncb_4, centbins_4, ""));
+
+  if (1) { // VnDelta vs. trigger pt -- 1 canvas/centrality bin
+    int cbins[] = { CentBin(0,10) , CentBin(10,20), CentBin(20,30), 
+                   CentBin(30,40), CentBin(40,50) };
+    int ptabins02[] = {PtBin(0.25, 0.5), PtBin(1.0, 1.5)};
+    int ptabins10[] = {PtBin(0.25, 0.5), PtBin(1.0, 1.5), PtBin(2.0, 2.5)};
+
+    cList->Add(DrawVnVsPtTrig(CentBin(0,2), 2, ptabins02, ""));
+
+    for (int i=0; i<5; i++) {
+      cList->Add(DrawVnVsPtTrig(cbins[i], 3, ptabins10, ""));
+    }
+      cout << Form("DrawVnVsPtTrig() ok") << endl;
+  }
+
+  if (saveOutput) {
+    SaveGraphs(grFile.Data());
+    cout << Form("SaveGraphs() ok") << endl;
+    SaveCanvases(cList, canvFile.Data()); 
+    cout << Form("SaveCanvases() ok") << endl;
+
+    if (1) // individual pdfs + multipage pdf
+      SaveCanvasesFromFile(canvFile.Data(), Form("%s/plots", outputDir), tag, "pdf");
+    if (0) // individual plot macros
+      SaveCanvasesFromFile(canvFile.Data(), Form("%s/plots", outputDir), tag, "C");
+  }
+
+  return;
+}
+
+
+void PlotCollection(int k, int i, int j, TString opt)
+{
+  cList->Add(SingleDrawPlain(k, i, j, "harmonics"));
+  cList->Add(SingleDraw(k, i, j, "color"));
+  cList->Add(SingleDraw(k, i, j, "global"));
+  if (opt.Contains("2D")) {
+    cList->Add(SingleDraw2D(k, i, j, "pl"));
+  }
+  return;
+}