PWGCF/Correlations/DPhi/FourierDecomposition/lrc/FourierPlus.C

   1
   2 #include "TROOT.h"
   3 #include "TF1.h"
   4 #include "TH1.h"
   5 #include "TH2.h"
   6 #include "TNtuple.h"
   7 #include "TFile.h"
   8 #include "TKey.h"
   9 #include "TRandom.h"
  10 #include "TMath.h"
  11 #include "TGraph.h"
  12 #include "TCanvas.h"
  13 #include "TLatex.h"
  14 #include "TBox.h"
  15 #include "TGraphErrors.h"
  16 #include "TStyle.h"
  17 #include "TPolyLine.h"
  18 #include "TLegend.h"
  19 #include "TLine.h"
  20 #include "TSystem.h"
  21 #include "TObjArray.h"
  22 #include "TGraph2D.h"
  23 #include "TGaxis.h"
  24 #include "TArrow.h"
  25 #include "TASImage.h"
  26 #include "TFitResult.h"
  27 #include "TStopwatch.h"
  28 #include "TProfile2D.h"
  29 #include "TRegexp.h"
  30 #include <iostream>
  31 #include <fstream>
  32 #include <vector>
  33
  34 // Plot utilities class - loaded by rootlogon.C + .rootrc
  35 #include "Utils.h"
  36 PlotUtils utils;
  37 const char* inputFileName = ""; // Assign in macro
  38 const char* cfDef = "cA"; // s, m, cA or cB. Only for phi projections.
  39 double minRidgeDEta = 0.8;
  40 double maxRidgeDEta = 1.799;
  41 const int gNMax = 12; // Global const -- largest n in VnDelta
  42
  43 // Don't do global fit to points outside these pt bin
  44 // indices. Globally scoped for availability in vntvna(). Reassigned
  45 // within individual functions.
  46 double minFitPtBin = 0, maxFitPtBin = 999;
  47
  48 // Switches
  49 //============================================================
  50 bool usePtBinCenter = 1;
  51 bool ispp = 0;
  52
  53 // Functions
  54 //============================================================
  55 void SetProtonProtonFlag(int val) { ispp = val; }
  56 void SetMinRidgeDEta(int val) { minRidgeDEta = val; }
  57 bool IsProtonProton() { return ispp; }
  58 void SetCFDef(char* def) { cfDef = def; }
  59 void SaveGraphs(const char* outputFileName, TString opt = "");
  60 TF1* DoGlobalFit(int n, int k, int includedPtBin1 = 0, int includedPtBin2 = 999,
  61                  const char* region = "RIDGE", const char* cfDef = "cA", TString opt = "");
  62 void VnDelta(int n, const TH1 &h, double &vnd, double &vnd_err, TString opt="");
  63 void VnDelta(int n, const TF1 &f, double &vnd, double &vnd_err, TString opt="");
  64 void MakeVnDVsQGraphs(int n1, int n2, int k, const char* region, const char* corrtype, TString opt="");
  65 TF1* Harmonic(TH1* h, int n, TString opt);
  66 TF1* HarmonicSum(TH1* h, int n1, int n2, TString opt="");
  67 TH1* Hist(const char* region, const char* type, int i, int j, int k, TString opt="");
  68 TH2* EtaPhiHist(int cent, int ptt, int pta);
  69 TH2* PttPtaHist();
  70 TH2F* Agreement2DHist(int k, int n);
  71 void ijkFromHistName(TH1* h, int& i, int& j, int& k);
  72
  73 TGraphErrors* VnDeltaVsN(int i, int j, int k, int nMax, TString opt);
  74 TGraphErrors* VnDeltaNFVsN(int i, int j, int k, int nMax, TString opt);
  75 TGraphErrors* GlobalvnVsN(int i, int k, int nMax, TString opt);
  76 TGraphErrors* VnDeltaVsPtAssoc(int i, int k, int n, TString opt="");
  77 TGraphErrors* VnDeltaVsPtTrig(int j, int k, int n, TString opt="");
  78 TGraphErrors* AgreementVsPtSum(int k, int n);
  79 TGraphErrors* VnDVsQ(int n, int k, const char* region = "RIDGE", const char* cfDef = "cA", TString opt="");
  80 TGraphErrors* GlobalFitVsPtAssoc(int i, int k, int n, TString opt=""); // v_n(ptt)v_n(pta) vs. assoc pt
  81 TGraphErrors* vnGFvsPt(int n, int k, int ipt1, int ipt2,
  82                        const char* region = "RIDGE",
  83                        const char* corrtype = "cA", TString opt = "");
  84 TF1* GlobalFitVsQ(int n, int k, int ipt1=0, int ipt2=999, const char* region="RIDGE",
  85                   const char* corrtype="cA", TString opt=""); // The real fit function
  86 //TF1* GlobalFitVsQSubset(int n, int k, int i, TString opt = "");
  87 TGraph* Luzumv1(int cent, const char* pid);
  88 double MomCorrection(int i, int j, int k);
  89
  90 // Global fit function - give to TF1 constructor
  91 double vntvna(double *x, double *par);
  92 double Chi2(TH1* h, TF1* f, TString opt = ""); // "", "ndf" or "prob"
  93 double ReducedChi2(int i, int k, int n);
  94 double SineError(int i, int j, int k); // RMS of <sin n D phi> vs. n about zero
  95 double vnGF(int n, int k, int ptBin, int includedPtBin1 = 0, int includedPtBin2 = 999,
  96             const char* region = "RIDGE", const char* cfDef = "cA", TString opt = "");
  97
  98 // Show +/- solutions in chi2 space by varying one fixed parameter (par)
  99 TCanvas* FitStudy(int icent, int n, int par);
 100 TCanvas* DrawQ(int k, int n, TString opt = "");
 101 TCanvas* SingleDraw(int cent, int ptt, int pta, TString opt="");
 102 TCanvas* SingleDrawEta(int cent, int ptt, int pta, TString opt);
 103 TCanvas* SingleDrawPlain(int cent, int ptt, int pta, TString opt);
 104 TCanvas* DrawVnVsPtTrig(int k, int npta, int ptabins[], TString opt = "");
 105 TCanvas* SingleDraw2D(int cent, int ptt, int pta, TString opt="");
 106 TCanvas* MultiDraw(int cent, TString opt);
 107 TCanvas* DrawVnFromGlobalFit(int n, int ipt1=0, int ipt2=999, int ncb = 999,
 108                              int centbins[] = 0, TString opt=""); // little v_n
 109 TCanvas* TwoPanelVert(double ydiv, const char* canvName, const char* canvTitle);
 110 TCanvas* DrawChi2vsPtTrig(int npta, int ptabins[], int ncb, int centbins[], TString opt);
 111 TCanvas* DrawAgreement(int icent, TString opt = "");
 112 TCanvas* Drawv1to5(int ncb, int centbins[], int ipt1=0, int ipt2=999, TString opt="");
 113 TCanvas* Drawv2to5(int ncb, int centbins[], int ipt1=0, int ipt2=999, TString opt="");
 114 TCanvas* DrawGlobalvnVsN(int npt, int ptbins[], int ncb, int centbins[], TString opt="");
 115
 116 // functions needed for above routine
 117 void initialize(const char* ifName = inputFileName);
 118 void AddPrelimStampToCurrentPad(float x1, float y1, float x2, float y2, TString opt="");
 119
 120 void SaveCanvases(TObjArray* canvases, const char* fileName);
 121 void SaveCanvasesFromFile(const char* rootFile, const char* targetDir, const char* tag, const char* fileType);
 122 TObjArray* GetObjectsFromFile(TFile& file, TString clname, TString dir="");
 123
 124
 125 // Global variables
 126 //============================================================
 127 TFile *fin=0;
 128 int maxcent;
 129 int centlow[100];  // low and high percentile for each selection
 130 int centhigh[100];
 131 int maxpt;
 132 double ptmean[100];
 133 double ptlow[100];
 134 double pthigh[100];
 135
 136 int colorsPale[] = {kOrange-3, kGreen-6, kAzure-9, kGray+1, kRed-9, kRed-7, kBlue-9, kViolet-9};
 137 int colorsDark[] = {kOrange-8, kGreen+1, kAzure+2, kBlack,  kCyan+1, kRed+0, kBlue+1, kViolet+1};
 138 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};
 139 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};
 140
 141 TLatex ltx;  // NDC
 142 TLatex ltx2; // absolute
 143
 144 // Global container for graphs, TF1s, etc.
 145 TObjArray* gList = new TObjArray();
 146 TObjArray* gFitList = new TObjArray();
 147 TF1* axisFn = 0;
 148
 149 TGraphErrors *gi, *gj, *gk;
 150 char* trigLabel(int i);
 151 char* asscLabel(int j);
 152 char* centLabel(int k);
 153 double MeanPt(int i, int j, int k, TString t_or_a = "t" /* "t" or "a" */,
 154               TString opt="" /* "" or "squared" */);
 155 int CentBin(double cntLo, double cntHi);
 156 int PtBin(double ptLo, double ptHi);
 157 int CentColor(double cen1, double cen2);
 158 int CentColorPale(double cen1, double cen2);
 159 int PtColor(double p1, double p2);
 160 int PtColorPale(double p1, double p2);
 161 TStopwatch timer;
 162
 163 // The business
 164 //============================================================
 165
 166 void SaveGraphs(const char* outputFileName, TString opt)
 167 {
 168   const char* mode = "recreate";
 169   if (opt.Contains("update"))
 170     mode = "update";
 171
 172   TFile* outputFile = new TFile(outputFileName, mode);
 173   gList->Write();
 174   gFitList->Write();
 175   outputFile->Close();
 176   return;
 177 }
 178
 179 void initialize(const char* ifName) {
 180
 181   // TH1F files with histograms for:
 182   // NSJET_y_PTTRIG_PTASSOC_CENT (tracks with |delta-eta| < 0.8)
 183   // RIDGE_y_PTTRIG_PTASSOC_CENT (tracks with |delta-eta| = 0.8 - 1.6)
 184
 185   if (fin) // don't keep re-initializing
 186     return;
 187
 188   fin = new TFile(ifName);
 189   if (!fin) {
 190     Error("FourierPlus - initialize()", "%s not found", ifName);
 191     gSystem->Exit(1);
 192   }
 193
 194   TList* info = (TList*)fin->Get("FileInfo");
 195   if (info)
 196     info->Print();
 197
 198   gi = (TGraphErrors*) fin->Get("TrigPtBins");
 199   gj = (TGraphErrors*) fin->Get("AsscPtBins");
 200   gk = (TGraphErrors*) fin->Get("EvCentBins");
 201
 202   maxpt = gi->GetN();
 203   maxcent = gk->GetN();
 204   gList->Add(gi);   gList->Add(gj);   gList->Add(gk);
 205
 206   if (gj->GetN() != maxpt)
 207     Warning("initialize()",
 208             "Trigger and assc. binning are not symmetric: %d vs. %d",
 209             maxpt, gj->GetN() );
 210
 211   for (int k=0; k<maxcent; k++) {
 212     centlow[k] = gk->GetX()[k] - gk->GetEX()[k];
 213     centhigh[k] = gk->GetX()[k] + gk->GetEX()[k];
 214     //    cout << centlow[k] << " " << centhigh[k] << endl;
 215   }
 216
 217   for (int i=0; i<maxpt; i++) {
 218     ptlow[i] = gi->GetX()[i] - gi->GetEX()[i];
 219     pthigh[i] = gi->GetX()[i] + gi->GetEX()[i];
 220     //    cout << ptlow[i] << " " << pthigh[i] << endl;
 221   }
 222
 223   // OK, leave hard-coded for now....
 224   ptmean[0] = 0.66; // 0.5-1 GeV
 225   ptmean[1] = 1.30; // 1-2 GeV
 226   ptmean[2] = 2.30; // 2-3 GeV
 227   ptmean[3] = 3.30; // 3-4 GeV
 228   ptmean[4] = 4.60; // 4-6 GeV
 229   ptmean[5] = 6.60; // 6-8 GeV
 230   ptmean[6] = 11.7; // 8-15 GeV
 231
 232   if (usePtBinCenter)
 233     for (int i=0; i<maxpt; i++) {
 234       ptmean[i] = gi->GetX()[i];
 235     }
 236
 237   gROOT->SetStyle("Plain");
 238   gStyle->SetOptTitle(0);
 239   gStyle->SetOptStat(0);
 240   ltx.SetNDC();
 241 }
 242
 243 TH1* Hist(const char* region, const char* type, int i, int j, int k, TString opt)
 244 {
 245   // Options:
 246   // "sys": return measured pts. w/systematic errors.
 247   // "hi_sys": meas. pts + sys, unmodified stat. errors.
 248   // "lo_sys": meas. pts - sys, unmodified stat. errors.
 249
 250   TH1* h_in=0, *h_out=0;
 251   const char* name = Form("PbPb/%s_%s_%d_%d_%d",region,type,i,j,k);
 252   h_in = (TH1*)fin->Get(name);
 253   if (!h_in) {
 254     Error("Hist()","%s not found, exiting.", name);
 255     gSystem->Exit(-1);
 256   }
 257
 258   h_out = (TH1*)h_in->Clone(Form("%s_%s", name, opt.Data()));
 259
 260   // Stat --> sys error bars. Error estimated as difference in Fourier
 261   // series (0<n<6) between same-event dist. and CF.
 262   if (opt.Contains("sys")) {
 263     const char* s = Form("PbPb/%s_s_%d_%d_%d",region,i,j,k);
 264     TH1* htmp = (TH1*)fin->Get(s);
 265     TH1* hs = 0;
 266     if (htmp)
 267       hs = (TH1*)htmp->Clone();
 268     else
 269       Error("Hist()", "No hist %s", s);
 270
 271     double c = double(hs->GetNbinsX())/hs->Integral();
 272     hs->Scale(c);
 273
 274     TF1* fs = HarmonicSum(hs,  1,5);
 275     TF1* fc = HarmonicSum(h_in,1,5);
 276
 277     for (int n=1; n<=hs->GetNbinsX(); n++) {
 278       double x = hs->GetBinCenter(n);
 279       double err = TMath::Abs(fs->Eval(x) - fc->Eval(x));
 280       double bc = h_in->GetBinContent(n);
 281
 282       if (opt.Contains("hi_sys"))
 283         h_out->SetBinContent(n, bc+err);
 284       else if (opt.Contains("lo_sys"))
 285         h_out->SetBinContent(n, bc-err);
 286       else {
 287         h_out->SetBinError(n, err);
 288         h_out->SetFillColor(kGray);
 289       }
 290     }
 291
 292     delete htmp;
 293     delete hs;
 294     delete fs;
 295     delete fc;
 296   }
 297
 298   return h_out;
 299 }
 300
 301 void ijkFromHistName(TH1* h, int& i, int& j, int& k)
 302 {
 303   TString tok, name(h->GetName());
 304   Ssiz_t from = 0;
 305   int counter = 0;
 306   while (name.Tokenize(tok, from, "_")) {
 307     if (counter==2) i = tok.Atoi();
 308     if (counter==3) j = tok.Atoi();
 309     if (counter==4) k = tok.Atoi();
 310     counter++;
 311   }
 312   if (0)
 313     cout << i<<j<<k << endl;
 314
 315   return;
 316 }
 317
 318 double PtBinCenterSum(int ptt, int pta)
 319 {
 320   if ((ptt<0 || ptt>=maxpt) || (pta>ptt))
 321     Error("PtBinCenterSum()","Bad arg(s): ptt %d pta %d", ptt, pta);
 322
 323   double pti = gi->GetX()[ptt];
 324   double ptj = gj->GetX()[pta];
 325   return pti+ptj;
 326 }
 327
 328 int GlobalIndex(int ptt, int pta)
 329 {
 330   // Find the global index from the trig and assc. pt bins.
 331   // It is useful that TMath::Binomial(n,k) returns 0 when n < k.
 332   if (pta > ptt) return -1;
 333   return TMath::Binomial(ptt+1, 2) + pta;
 334 }
 335
 336 void ijFromGlobalIndex(int q, int& ptt, int& pta)
 337 {
 338   // Assign ptt and pta from the global index q. This function is
 339   // effectively the inverse of GlobalIndex(i,j).
 340   for (int i=0; i<maxpt; i++) {
 341     for (int j=0; j<=i; j++) {
 342       if (GlobalIndex(i,j)==q) {
 343         ptt = i; pta = j;
 344         return;
 345       }
 346     }
 347   }
 348   return;
 349 }
 350
 351 // Needs to be updated! Has hard-coded numbers like 27.9. -AMA 11/15/2011
 352 TCanvas* FitStudy(int icent, int n, int par)
 353 {
 354   // First get the VnDelta vs q graph...
 355   TGraphErrors* gVn = VnDVsQ(n, icent);
 356
 357   // Then setup fit function
 358   TF1 *fn = new TF1(Form("fn_cent%dto%d_n%d", centlow[icent],centhigh[icent], n),
 359                     vntvna, 0.0, 27.9, maxpt);
 360   fn->SetNpx(1000);
 361   fn->SetLineColor(kRed);
 362
 363   TCanvas* c = new TCanvas(Form("chi2_cent%dto%d_n%d", centlow[icent],centhigh[icent], n),
 364                            Form("chi2 test for %s, varying v_{%d}(%.1f)",
 365                                 centLabel(icent), n, ptmean[par]),
 366                            700, 1200);
 367   c->Divide(1, 3, 0.001, 0.001);
 368
 369   TGraphErrors *g = new TGraphErrors();
 370   TGraphErrors *chi2 = new TGraphErrors();
 371   chi2->SetTitle(Form("Total #chi^{2} for n=%d;v_{%d}(p_{T}=%.2f);total #chi^{2}", n,n, ptmean[par]));
 372   utils.set_tgraph_props(g, kGray, kGray, kFullCircle, 1.3);
 373   utils.set_tgraph_props(chi2, kBlue, kBlue, kFullCircle, 0.8);
 374   g->SetLineWidth(3);
 375
 376   double bestvn     = vnGF(n,icent,par,0,999,"RIDGE","cA","");
 377   double bestvn_err = vnGF(n,icent,par,0,999,"RIDGE","cA","err");
 378
 379   int nsteps = 500;
 380
 381   double parmax = TMath::Abs(2*(bestvn + bestvn_err));
 382   double parmin = -parmax;
 383
 384   // First panel: step thru one fixed vn
 385   c->cd(1);
 386   TH1F* hf1 = gPad->DrawFrame(0.0, 1.2*parmin, 12.5, 1.2*parmax);
 387   hf1->SetTitle(Form("title;p_{T};v_{%d}(p_{T})", n));
 388
 389   for (int step=0; step<=nsteps; step++) {
 390     double stepsize =  (parmax-parmin)/nsteps;
 391     double val = parmin + step*stepsize;
 392
 393     for (int ip=0;ip<maxpt;ip++) fn->SetParameter(ip,0.0);
 394     fn->FixParameter(par, val);
 395
 396     gVn->Fit(fn, "QR");
 397     for (int i=0;i<maxpt;i++) {
 398       g->SetPoint(i, ptmean[i], fn->GetParameter(i));
 399       g->SetPointError(i, 0.0, fn->GetParError(i));
 400     }
 401     if (step%50==0)
 402       g->DrawClone("elpsame");
 403
 404     // store chi2 at this val
 405     chi2->SetPoint(step, val, fn->GetChisquare());
 406   }
 407   // Finally, draw once with no fixed parameters...
 408   utils.set_tgraph_props(g, kBlack, kBlack, kFullCircle, 1.3);
 409   for (int ip=0;ip<maxpt;ip++) fn->ReleaseParameter(ip);
 410   gVn->Fit(fn,"QR");
 411   for (int i=0;i<maxpt;i++) {
 412     g->SetPoint(i, ptmean[i], fn->GetParameter(i));
 413     g->SetPointError(i, 0.0, fn->GetParError(i));
 414   }
 415   g->DrawClone("elpsame");
 416   ltx.SetTextSize(0.08);
 417   ltx.DrawLatex(0.6, 0.8, Form("v_{%d}, %s", n, centLabel(icent)));
 418
 419   // Plot chi2 over full range
 420   c->cd(2);
 421   c->Update();
 422   chi2->DrawClone("alp");
 423
 424   // and zoomed to see 1 sigma interval
 425   c->cd(3);
 426   c->Update();
 427   double ymin = 1e9;//chi2->GetHistogram()->GetMinimum();
 428   for (int j=0; j<chi2->GetN(); j++)
 429     if (chi2->GetY()[j] < ymin)
 430       ymin = chi2->GetY()[j];
 431   double nsigma = n==2? 3 : 3;
 432   TH1F* hf3 = gPad->DrawFrame(parmin, ymin - 1, parmax, ymin + nsigma);
 433   hf3->SetTitle(Form("Total #chi^{2} for n=%d;v_{%d}(p_{T}=%.2f);total #chi^{2}", n,n, ptmean[par]));
 434   chi2->DrawClone("lpsame");
 435   TLine best, onesig;
 436   best.SetLineWidth(3);
 437   onesig.SetLineWidth(3);
 438   onesig.SetLineColor(kGray);
 439   best.DrawLine(parmin, ymin, parmax, ymin);
 440   onesig.DrawLine(parmin, ymin+1, parmax, ymin+1);
 441   TLatex lbl;
 442   lbl.DrawLatex(parmin + 0.05*(parmax-parmin), ymin+0.1, Form("Best v_{%d}(%.2f)",n,ptmean[par]));
 443   lbl.DrawLatex(parmin + 0.05*(parmax-parmin), ymin+1.1, "+ 1#sigma");
 444
 445   return c;
 446 }
 447
 448 TCanvas* MultiDraw(int cent, TString opt)
 449 {
 450   // TODO: make this function much more flexible. Pass in an array of
 451   //  bins or histos and m,n ints to plot in an m x n grid.
 452   TCanvas* c = new TCanvas(Form("cfmulti_%s_%d",opt.Data(),cent),
 453                            Form("Corr. Functions - %s", centLabel(cent)),
 454                            1200, 900);
 455   c->Divide(maxpt, maxpt-2, 0, 0);
 456
 457   int ipad=1;
 458
 459   // Combine first 3 rows (i=0,1,2) to save space
 460   // Then draw the rest as usual
 461   for (int i=0; i<maxpt; i++) {
 462     for (int j=0; j<maxpt; j++) {
 463       if(j<=i) {
 464         TCanvas* cs = 0;
 465         if (opt.Contains("2D"))
 466           cs = SingleDraw2D(cent, i, j, "");
 467         else if (opt.Contains("deta"))
 468           cs = SingleDrawEta(cent, i, j, "");
 469         else
 470           cs = SingleDraw(cent, i, j, "");
 471         c->cd(ipad);
 472         cs->DrawClonePad();
 473         ipad++;
 474         if (ipad==maxpt) ipad++;
 475       }
 476       else if (i > 2) ipad++;
 477     }
 478   }
 479
 480   // Add some labeling in empty pads
 481   ltx.SetTextSize(0.15);
 482   c->cd(1*maxpt);
 483   ltx.DrawLatex(0.1, 0.8, "p_{T}^{trig}: 0.5-1,");
 484   ltx.DrawLatex(0.1, 0.6, "1-2, 2-3 GeV/c");
 485   c->cd(2*maxpt);
 486   ltx.DrawLatex(0.1, 0.8, "p_{T}^{trig}:");
 487   ltx.DrawLatex(0.1, 0.6, "3-4 GeV/c");
 488   c->cd(3*maxpt);
 489   ltx.DrawLatex(0.1, 0.8, "p_{T}^{trig}:");
 490   ltx.DrawLatex(0.1, 0.6, "4-6 GeV/c");
 491   c->cd(4*maxpt);
 492   ltx.DrawLatex(0.1, 0.8, "p_{T}^{trig}:");
 493   ltx.DrawLatex(0.1, 0.6, "6-8 GeV/c");
 494
 495   c->cd(2*maxpt-2);
 496   ltx.SetTextSize(0.2);
 497   ltx.DrawLatex(0.1, 0.7, "Pb+Pb");
 498   ltx.DrawLatex(0.1, 0.4, "LHC10h");
 499   ltx.DrawLatex(0.1, 0.1, centLabel(cent));
 500
 501   return c;
 502 }
 503
 504 TH2* EtaPhiHist(int cent, int ptt, int pta)
 505 {
 506  initialize();
 507   TH2* htmp = (TH2*)fin->Get(Form("PbPb/ETAPHI_c_%1d_%1d_%1d",ptt,pta,cent));
 508   if (!htmp) {
 509     Error("EtaPhiHist()",
 510           "Problem getting pointer to histogram. cent %d ptt %d pta %d",
 511           cent, ptt, pta);
 512     return 0;
 513   }
 514   TH2* h = (TH2*)htmp->Clone(Form("EtaPhi_c_%1d_%1d_%1d",ptt,pta,cent));
 515   TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis(), *az = h->GetZaxis();
 516   ax->SetTitle("#Delta#phi");
 517   ay->SetTitle("#Delta#eta");
 518   az->SetTitle("C(#Delta#phi, #Delta#eta)");
 519   ay->SetRangeUser(-maxRidgeDEta, maxRidgeDEta);
 520   ax->SetTitleOffset(1.3);
 521   ay->SetTitleOffset(1.3);
 522   az->SetTitleOffset(1.4);
 523
 524   return h;
 525 }
 526
 527 TCanvas* SingleDraw2D(int cent, int ptt, int pta, TString opt)
 528 {
 529   initialize();
 530   TCanvas* c = 0; // The canvas to be returned
 531   double labelScale = 1.3;
 532   double nDiv = 5;
 533
 534   // See if we have already made this once
 535   const char* name;
 536   if (opt=="")
 537     name = Form("etaphi_cent%dto%d_%d_%d",
 538                 centlow[cent],centhigh[cent],ptt,pta );
 539   else
 540     name = Form("etaphi_cent%dto%d_%d_%d_%s",
 541                 centlow[cent],centhigh[cent],ptt,pta, opt.Data() );
 542
 543   const char* title = Form("dphideta_cent%dto%d_pTtrig%.2gto%.2g_pTassoc%.2gto%.2g_%s",
 544                            centlow[cent],centhigh[cent],
 545                            ptlow[ptt], pthigh[ptt],
 546                            ptlow[pta], pthigh[pta],
 547                            opt.Data());
 548
 549   if (gDirectory->FindObject(name)) {
 550     c = (TCanvas*) gDirectory->FindObject(name);
 551     return c;
 552   }
 553   else
 554     c = new TCanvas(name, title, 1);
 555
 556   TH2* h = EtaPhiHist(cent, ptt,pta);
 557   if (!h) {
 558     Error("SingleDraw2D()",
 559           "Problem getting pointer to histogram. cent %d ptt %d pta %d",
 560           cent, ptt, pta);
 561     return 0;
 562   }
 563   utils.make_nice_axes(c, h, labelScale, nDiv, nDiv);
 564   TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis(), *az = h->GetZaxis();
 565   ax->SetTitle("#Delta#phi [rad]");
 566   ay->SetTitle("#Delta#eta");
 567   az->SetTitle("C(#Delta#phi, #Delta#eta)");
 568   ay->SetRangeUser(-maxRidgeDEta, maxRidgeDEta);
 569   ax->SetTitleOffset(1.3);
 570   ay->SetTitleOffset(1.3);
 571   az->SetTitleOffset(1.4);
 572
 573   double czlo=0, czhi=5.0;
 574   if (opt.Contains("zoom"))
 575     az->SetRangeUser(czlo, czhi);
 576
 577   if (opt.Contains("colz"))
 578     h->Draw("colz");
 579   else
 580     h->Draw("surf1");
 581
 582   ltx.SetTextSize(0.05);
 583   if (opt.Contains("colz")) {
 584     if (opt.Contains("zoom"))
 585       ltx.DrawLatex(0.6, 0.9, Form("C(#Delta#phi) (C(#Delta#phi) < %.3g)", czhi));
 586     else
 587       ltx.DrawLatex(0.6, 0.9, Form("C(#Delta#phi) (full range)"));
 588   }
 589
 590   // Label pt, centrality, energy
 591   TString ptStr(Form("#splitline"
 592                      "{%.3g < p_{T}^{t} < %.3g GeV/c}"
 593                      "{%.3g < p_{T}^{a} < %.3g GeV/c}",
 594                      ptlow[ptt], pthigh[ptt],
 595                      ptlow[pta], pthigh[pta]));
 596   ltx.SetTextSize(0.07);
 597   ltx.DrawLatex(0.01, 0.86, ptStr.Data());
 598   if (ispp)
 599     ltx.DrawLatex(0.65, 0.9, Form("#splitline{pp 2.76 TeV}{}"));
 600   else
 601     ltx.DrawLatex(0.65, 0.9, Form("#splitline{Pb-Pb 2.76 TeV}{%s}", centLabel(cent)));
 602
 603   return c;
 604 }
 605
 606 TCanvas* SingleDrawEta(int cent, int ptt, int pta, TString opt)
 607 {
 608   TCanvas* c = 0;
 609   int lineWidth  = opt.Contains("small") ? 1 : 4;
 610   double mrkSize = opt.Contains("small") ? 0.5 : 1.5;
 611
 612   const char* name = Form("deta_cent%dto%d_%d_%d",
 613                           centlow[cent],centhigh[cent],ptt,pta);
 614
 615   // See if we have already made this once
 616   if (gDirectory->FindObject(name)) {
 617     c = (TCanvas*) gDirectory->FindObject(name);
 618     Info("SingleDrawEta()", "%s already created", c->GetName());
 619     return c;
 620   }
 621   else
 622     c = new TCanvas(name, name, 600, 600);
 623
 624   TH1 *hNS = 0, *hAS = 0;    // divide-then-project
 625   TH1 *hNS2 = 0, *hAS2 = 0;  // project-then-divide
 626
 627   bool includeProjThenDiv = 0;
 628
 629   if (1) { // Recreate CFs using divide-then-project method
 630
 631     TH2* hc = (TH2*)fin->Get(Form("PbPb/ETAPHI_c_%1d_%1d_%1d",ptt,pta,cent));
 632     if (!hc) {
 633       Error("SingleDrawEta()", "Problem getting hc\n");
 634       gSystem->Exit(-1);
 635     }
 636
 637     int bin1 = hc->GetXaxis()->FindBin(-TMath::PiOver2()+0.001);
 638     int bin2 = hc->GetXaxis()->FindBin(TMath::PiOver2()-0.001);
 639     int bin3 = hc->GetXaxis()->FindBin(TMath::PiOver2()+0.001);
 640     int bin4 = hc->GetXaxis()->FindBin(3*TMath::PiOver2()-0.001);
 641
 642     bool smallEtaRange = false;
 643     if (smallEtaRange) {
 644       bin1 = hc->GetXaxis()->FindBin(-TMath::PiOver2()+0.001);
 645       bin2 = hc->GetXaxis()->FindBin(TMath::PiOver2()-0.001);
 646       bin3 = hc->GetXaxis()->FindBin(TMath::Pi()-0.5);
 647       bin4 = hc->GetXaxis()->FindBin(TMath::Pi()+0.5);
 648     }
 649
 650     const char* ns_name = Form("ns_deta_cent%dto%d_%d_%d",
 651                                centlow[cent],centhigh[cent],ptt,pta);
 652     const char* as_name = Form("as_deta_cent%dto%d_%d_%d",
 653                                centlow[cent],centhigh[cent],ptt,pta);
 654
 655     hNS = hc->ProjectionY(ns_name, bin1, bin2, "e");
 656     hAS = hc->ProjectionY(as_name, bin3, bin4, "e");
 657
 658     hNS->Scale(1./(bin2-bin1+1));
 659     hAS->Scale(1./(bin4-bin3+1));
 660
 661   }
 662
 663   if (includeProjThenDiv) { // As created by Project.C (project-then-divide)
 664     const char* ns = Form("PbPb/ETA_NS_c_%1d_%1d_%1d",ptt,pta,cent);
 665     const char* as = Form("PbPb/ETA_AS_c_%1d_%1d_%1d",ptt,pta,cent);
 666     hNS2 = (TH1*)fin->Get(ns);
 667     hAS2 = (TH1*)fin->Get(as);
 668     if (!hNS2) {
 669       Error("SingleDrawEta()", "Problem getting %s\n", ns);
 670       gSystem->Exit(-1);
 671     }
 672     if (!hAS2) {
 673       Error("SingleDrawEta()", "Problem getting %s\n", as);
 674       gSystem->Exit(-1);
 675     }
 676   }
 677
 678   utils.set_hist_props(hNS, kBlack, kNone, kBlack, kFullDotLarge, mrkSize);
 679   utils.set_hist_props(hAS, kRed, kNone, kRed, kFullDotLarge, mrkSize);
 680   utils.set_ylimits(hNS, hAS);
 681   hNS->SetLineWidth(lineWidth);
 682   hAS->SetLineWidth(lineWidth);
 683
 684   // hNS->Draw("hist");
 685   // hAS->Draw("histsame");
 686   hNS->DrawClone("ep");
 687   hAS->DrawClone("epsame");
 688
 689   if (includeProjThenDiv) {
 690     utils.set_hist_props(hNS2, kGray+2, kNone, kGray+2, kFullDotLarge, mrkSize);
 691     utils.set_hist_props(hAS2, kCyan, kNone, kCyan, kFullDotLarge, mrkSize);
 692     hNS2->SetLineWidth(lineWidth);
 693     hAS2->SetLineWidth(lineWidth);
 694
 695     hNS2->DrawClone("epsame");
 696     hAS2->DrawClone("epsame");
 697   }
 698
 699   return c;
 700 }
 701
 702 TCanvas* SingleDrawPlain(int cent, int ptt, int pta, TString opt)
 703 {
 704   initialize();
 705   const char* region = "RIDGE";
 706   const char* yTitle = Form("%.2g < |#Delta#eta| < %.2g",
 707                             minRidgeDEta, maxRidgeDEta);
 708   if (opt.Contains("ALL")) {
 709     region = "ALL";
 710     yTitle = Form("%.2g < |#Delta#eta| < %.2g", 0.0, maxRidgeDEta);
 711   }
 712   if (opt.Contains("NSJET")) {
 713     region = "NSJET";
 714     yTitle = Form("%.2g < |#Delta#eta| < %.2g", 0.0, minRidgeDEta);
 715   }
 716
 717   const char* name = Form("%s_cent%dto%d_%d_%d_%s", region,
 718                           centlow[cent],centhigh[cent],ptt,pta, opt.Data());
 719   const char* title = Form("dphiPlain_%s_cent%dto%d_pTtrig%.2gto%.2g_pTassoc%.2gto%.2g_%s",
 720                            region,
 721                            centlow[cent],centhigh[cent],
 722                            ptlow[ptt], pthigh[ptt],
 723                            ptlow[pta], pthigh[pta],
 724                            opt.Data());
 725   TCanvas* c = 0;
 726   if (gDirectory->FindObject(name)) {
 727
 728     c = (TCanvas*) gDirectory->FindObject(name);
 729     Info("SingleDrawPlain()", "%s already created", c->GetName());
 730     return c;
 731   }
 732   else
 733     c = new TCanvas(name, title, 700, 700);
 734   c->cd();
 735
 736   TH1* h_sys = Hist(region, cfDef, ptt,pta,cent, "sys");
 737   TH1* h     = Hist(region, cfDef, ptt,pta,cent);
 738   h->SetName("h");
 739
 740   TAxis* ax = h->GetXaxis();
 741   TAxis* ay = h->GetYaxis();
 742
 743   h->SetTitle(Form(";#Delta#phi [rad];C(#Delta#phi),  %s", yTitle));
 744   h->SetLineWidth(2);
 745   ax->SetNdivisions(208);
 746   ay->SetNdivisions(208);
 747   ax->SetLabelSize(0.04);
 748   ay->SetLabelSize(0.04);
 749   ax->SetTitleSize(0.05);
 750   ay->SetTitleSize(0.05);
 751   ax->SetTitleOffset(1.5);
 752   ay->SetTitleOffset(-1.9);
 753   c->SetFrameLineWidth(2);
 754   c->SetLeftMargin(0.2);
 755   c->SetBottomMargin(0.15);
 756   c->SetTopMargin(0.01);
 757   c->SetRightMargin(0.01);
 758   ax->CenterTitle();
 759   ay->CenterTitle();
 760   ax->SetTicks("+-");
 761   ay->SetTicks("+-");
 762
 763   utils.set_ylimits(h, h, 0.05);
 764   utils.set_hist_props(h, kBlack, kNone, kBlack, kFullDotLarge, 1.4);
 765   h->Draw();
 766   h_sys->Draw("e2psame");
 767   h->Draw("same");
 768   if (opt.Contains("harmonics"))
 769     for (int n=1; n<6; n++)
 770       Harmonic(h, n, "draw");
 771   if (opt.Contains("sum"))
 772     HarmonicSum(h,1,5,"draw");
 773
 774   ltx.SetTextSize(0.045);
 775   double top = 0.92, gap = 0.08;
 776   double x1 = (TString(region)=="RIDGE" && (ptlow[ptt] > 4 || ispp))? 0.25 : 0.5;
 777   ltx.DrawLatex(x1, top-0*gap, Form("p_{T}^{t} %s GeV/c", trigLabel(ptt)));
 778   ltx.DrawLatex(x1, top-1*gap, Form("p_{T}^{a} %s GeV/c", asscLabel(pta)));
 779   if (ispp)
 780   ltx.DrawLatex(x1, top-2*gap, Form("pp 2.76 TeV"));
 781   else
 782   ltx.DrawLatex(x1, top-2*gap, Form("Pb-Pb 2.76 TeV, %s", centLabel(cent)));
 783   if (0)
 784     ltx.DrawLatex(0.22, 0.92, Form("|#Delta#eta| > %.1f", minRidgeDEta));
 785
 786   // ltx.SetTextColor(kBlack);
 787   // ltx.SetTextSize(0.03);
 788   // ltx.DrawLatex(0.7, 0.18, "statistical error only");
 789
 790   return c;
 791 }
 792
 793 TF1* Harmonic(TH1* h, int n, TString opt)
 794 {
 795   TF1* f = 0;
 796   double vndelta=0, vndelta_err=0;
 797   double phiMin = -0.5*TMath::Pi();
 798   double phiMax = +1.5*TMath::Pi();
 799   double b0 = h->Integral() / h->GetNbinsX();
 800   int color = opt.Contains("gray")? kGray : colorsPale[n-1];
 801
 802   VnDelta(n, *h, vndelta, vndelta_err);
 803
 804   f = new TF1(Form("%scos%ddphi", h->GetName(), n),
 805               "[0]*(1. + 2.*[1]*cos([2]*x))",
 806               phiMin, phiMax);
 807     f->SetParameters(b0, vndelta, n);
 808     f->SetLineColor(color);
 809     f->SetLineWidth(3);
 810     if (opt.Contains("draw"))
 811       f->Draw("same");
 812
 813     return f;
 814 }
 815
 816 TF1* HarmonicSum(TH1* h, int n1, int n2, TString opt)
 817 {
 818   if (n1<1)
 819     Warning("HarmonicSum()", "n1=%d but should be > 0", n1);
 820
 821   double phiMin = -0.5*TMath::Pi();
 822   double phiMax = +1.5*TMath::Pi();
 823
 824   // Function definition string
 825   TString fdef("[0]*(1.0 + 2.0*(");
 826   for (int n=n1; n<=n2; n++) {
 827     fdef.Append(Form(" [%d]*cos(%d*x) ", n, n));
 828     fdef.Append( n<n2 ? "+" : "))");
 829   }
 830
 831   // Normalization
 832   TF1 *ff = new TF1("ff", fdef.Data(), phiMin,phiMax);
 833   ff->SetParameter(0, h->Integral()/h->GetNbinsX());
 834
 835   // Set parameters
 836   for (int n=n1; n<=n2; n++) {
 837
 838     // Sum of vnt x vna from global fit
 839     if (opt.Contains("global") ) {
 840       int i=-1, j=-1, k=-1;
 841       ijkFromHistName(h,i,j,k);
 842       double vnt = vnGF(n,k,i,0,999,"RIDGE","cA","");
 843       double vna = vnGF(n,k,j,0,999,"RIDGE","cA","");
 844       ff->SetParameter(n, vnt*vna); // TODO: consider including errors
 845       ff->SetLineColor(kRed+1);
 846     }
 847     else { // sum of directly-extracted pair moments
 848       double Vn=0, VnErr=0;
 849       VnDelta(n, *h, Vn, VnErr);
 850       ff->SetParameter(n, Vn);
 851       ff->SetLineColor(kBlack); // ff->SetLineColor(kAzure);
 852       //      ff->SetLineWidth(4);
 853       ff->SetLineStyle(kDashed);
 854       if (n2==10) {
 855         ff->SetLineColor(kGreen);
 856       }
 857     }
 858   }
 859
 860   if (opt.Contains("draw"))
 861     ff->Draw("same");
 862
 863   return ff;
 864 }
 865
 866 TCanvas* SingleDraw(int cent, int ptt, int pta, TString opt)
 867 {
 868   int lineWidth  = opt.Contains("small") ? 1 : 4;
 869   double mrkSize = opt.Contains("small") ? 0.5 : 1.5;
 870   TString dEtaLabel = "";
 871
 872   const char* region = "RIDGE";
 873   if (opt.Contains("ALL")) {
 874     region = "ALL";
 875     dEtaLabel = Form("%.2g < |#Delta#eta| < %.2g", 0.0, maxRidgeDEta);
 876   }
 877   if (opt.Contains("NSJET")) {
 878     region = "NSJET";
 879     dEtaLabel = Form("%.2g < |#Delta#eta| < %.2g", 0.0, minRidgeDEta);
 880   }
 881   if (TString(region)=="RIDGE")
 882     dEtaLabel = Form("%.2g < |#Delta#eta| < %.2g", minRidgeDEta, maxRidgeDEta);
 883
 884   initialize();
 885
 886   TCanvas* c1 = 0; // The canvas to be returned
 887
 888   // See if we have already made this once
 889   TString name = Form("dphi_etamin%02d_cent%dto%d_%d_%d", (int)(10*minRidgeDEta),
 890                           centlow[cent],centhigh[cent],ptt,pta);
 891   TString title = Form("dphiFourier_%s_cent%dto%d_pTtrig%.2gto%.2g_pTassoc%.2gto%.2g",
 892                            region,
 893                            centlow[cent],centhigh[cent],
 894                            ptlow[ptt], pthigh[ptt],
 895                            ptlow[pta], pthigh[pta]);
 896
 897   if (!opt.IsNull()) {
 898     name.Append("_"); name.Append(opt);
 899     title.Append("_"); title.Append(opt);
 900   }
 901
 902   if (gDirectory->FindObject(name.Data())) {
 903
 904     c1 = (TCanvas*) gDirectory->FindObject(name.Data());
 905     Info("SingleDraw()", "%s already created", c1->GetName());
 906     return c1;
 907   }
 908   else
 909     c1 = TwoPanelVert(0.3, name.Data(), title.Data());
 910
 911   c1->cd(1);
 912   gPad->SetTopMargin(0.1);
 913   c1->cd();
 914
 915   TH1* h = Hist(region, cfDef, ptt,pta,cent);
 916   TH1* h_sys = Hist(region, cfDef, ptt,pta,cent, "sys");
 917   utils.set_hist_props(h, kBlack, kNone, kBlack, kFullDotLarge, mrkSize);
 918   h->SetLineWidth(opt.Contains("small") ? 1 : 2);
 919   h->GetXaxis()->SetTitle("#Delta#phi [rad]");
 920   h->GetYaxis()->SetTitle("C(#Delta#phi)");
 921
 922   TF1* ff5 = HarmonicSum(h, 1, 5, "");
 923   ff5->SetLineWidth(lineWidth);
 924   TF1* fflowg = HarmonicSum(h, 1, 5, "global");
 925   fflowg->SetLineWidth(lineWidth);
 926   TF1* ff10 = HarmonicSum(h, 1, 10, "");
 927   ff10->SetLineWidth(lineWidth);
 928   TF1* ff12 = HarmonicSum(h, 1, 12, "");
 929   ff12->SetLineWidth(lineWidth);
 930
 931   c1->cd(1);
 932
 933   int nXticks = 404, nYticks = 210;
 934   double pad1Scale = 1.6;
 935   double pad2Scale = 2.3;
 936   utils.make_nice_axes(c1, h, pad1Scale, nXticks, nYticks, 0.05, 0.01);
 937   h->GetXaxis()->SetTitleOffset(1.2);
 938   h->GetYaxis()->SetTitleOffset(1.5);
 939   utils.set_ylimits(h,h);
 940
 941   h->Draw("ep");
 942   h_sys->Draw("e2psame");
 943   h->Draw("epsame");
 944
 945   if (opt.Contains("gray"))
 946     for (int n=1; n<6; n++)
 947       Harmonic(h, n, "gray, draw");
 948   else if (opt.Contains("color"))
 949     for (int n=1; n<6; n++)
 950       Harmonic(h, n, "draw");
 951
 952   if (opt.Contains("upto10")) {
 953     ff10->SetLineWidth(lineWidth+2);
 954     ff10->SetLineColor(kBlack);
 955     ff10->DrawClone("same");
 956     ff10->SetLineWidth(lineWidth);
 957     ff10->SetLineColor(kGreen);
 958     ff10->SetLineStyle(kDashed);
 959     ff10->Draw("same");
 960   }
 961
 962   //  ff10->Draw("same");
 963   ff5->Draw("same");
 964   h->Draw("p,e,same");
 965
 966   if (opt.Contains("global")) {
 967     fflowg->Draw("same");
 968   }
 969
 970   ltx.SetTextSize(0.06);
 971   double top = 0.82, gap = 0.08;
 972   double x1 = (TString(region)=="RIDGE" && ptlow[ptt] > 6)? 0.25 : 0.6;
 973   if (ispp)
 974     ltx.DrawLatex(0.2, 0.92, Form("pp 2.76 TeV"));
 975   else
 976     ltx.DrawLatex(0.2, 0.92, Form("Pb-Pb 2.76 TeV, %s central", centLabel(cent)));
 977   ltx.DrawLatex(x1, top-0*gap, Form("%.3g < p_{T}^{t} < %.3g GeV/c",
 978                                     ptlow[ptt], pthigh[ptt]));
 979   ltx.DrawLatex(x1, top-1*gap, Form("%.3g < p_{T}^{a} < %.3g GeV/c",
 980                                     ptlow[pta], pthigh[pta]));
 981   // ltx.DrawLatex(x1, top-0*gap, Form("p_{T}^{trig} %s GeV/c", trigLabel(ptt)));
 982   // ltx.DrawLatex(x1, top-1*gap, Form("p_{T}^{assoc} %s GeV/c", asscLabel(pta)));
 983   ltx.DrawLatex(x1, top-2*gap, dEtaLabel.Data());
 984
 985
 986   c1->cd(2);
 987
 988   TH1F *hsub = static_cast <TH1F *> (h->Clone());
 989   TH1F *hsub10 = static_cast <TH1F *> (h->Clone());
 990   TH1F *hsubg = (TH1F*) (h->Clone());
 991
 992   utils.set_hist_props(hsub10,  kBlack, kNone, kBlack, kOpenCircle, mrkSize);
 993   utils.set_hist_props(hsub,  kAzure, kNone, kAzure,   kFullSquare, mrkSize);
 994   utils.set_hist_props(hsubg, kRed+1, kNone, kRed+1, kFullSquare, mrkSize);
 995
 996   for (int i=1;i<=hsub->GetNbinsX();i++) {
 997     double val, err;
 998     val = h->GetBinContent(i) * h->GetBinWidth(i);
 999     err = h->GetBinError(i)   * h->GetBinWidth(i);
1000     double dphi1 = h->GetXaxis()->GetBinLowEdge(i);
1001     double dphi2 = h->GetXaxis()->GetBinUpEdge(i);
1002     double fSum5, fSum10, fSumGlobal;
1003     fSum5      = ff5->Integral(dphi1, dphi2);
1004     fSum10     = ff10->Integral(dphi1, dphi2);
1005     fSumGlobal = fflowg->Integral(dphi1, dphi2);
1006     hsub->SetBinContent(i,   val/fSum5);
1007     hsub->SetBinError(i,     err/fSum5);
1008     hsub10->SetBinContent(i, val/fSum10);
1009     hsub10->SetBinError(i,   err/fSum10);
1010     hsubg->SetBinContent(i,  val/fSumGlobal);
1011     hsubg->SetBinError(i,    err/fSumGlobal);
1012   }
1013   utils.make_nice_axes(c1, hsub, pad2Scale, nXticks, 5, 0.05, 0.01);
1014   if (opt.Contains("global"))
1015     utils.set_ylimits(hsub, hsubg, 0.2);
1016   else
1017     utils.set_ylimits(hsub, hsub10, 0.2);
1018
1019   hsub->SetTitle(Form(";#Delta#phi [rad];ratio"));
1020   //  hsub->SetTitle(Form(";#Delta#phi [rad];data/#Sigma_{n=1}^{%s}", opt.Contains("global")?"5":"N"));
1021   hsub->GetXaxis()->SetTitleOffset(1.25);
1022   hsub->GetYaxis()->SetTitleOffset(0.7);
1023   hsub->GetXaxis()->SetTicks("+-");   hsub->GetXaxis()->SetTickLength(0.08);
1024   hsub10->SetTitle(";#Delta#phi [rad];data/#Sigma_{n=1}^{N}");
1025   hsub10->GetXaxis()->SetTitleOffset(1.25);
1026   hsub10->GetYaxis()->SetTitleOffset(0.7);
1027   hsubg->SetTitle(";#Delta#phi [rad];data/calc");
1028   hsubg->GetXaxis()->SetTitleOffset(1.25);
1029   hsubg->GetYaxis()->SetTitleOffset(0.7);
1030
1031   TLegend* leg = new TLegend(0.2, 0.01, 0.38, 0.38);
1032   leg->SetFillColor(kNone);
1033   leg->SetBorderSize(1);
1034   if (opt.Contains("global")) {
1035     leg->AddEntry(hsub, "#LTcos(n#Delta#phi)#GT", "epl");
1036     leg->SetBorderSize(0);
1037   }
1038   else
1039     leg->AddEntry(hsub, "1#leqn#leq5", "epl");
1040
1041   TLine unity;
1042   hsub->Draw("ep");
1043
1044   if (opt.Contains("global")) {
1045     hsubg->Draw("epsame");
1046     leg->AddEntry(hsubg, "v(p_{Tt})v(p_{Ta})", "epl");
1047   }
1048   else if (opt.Contains("upto10")) {
1049     TH1* hsub10green = (TH1*)hsub10->Clone();
1050     utils.set_hist_props(hsub10green,  kBlack, kNone, kGreen, kFullCircle, mrkSize);
1051     hsub10green->Draw("epsame");
1052     hsub10->Draw("epsame");
1053     leg->AddEntry(hsub10green, "1#leqn#leq10", "epl");
1054   }
1055
1056   unity.DrawLine(-TMath::PiOver2(), 1.0, 3*TMath::PiOver2(), 1.0);
1057
1058   c1->cd(2);
1059   if (opt.Contains("global"))
1060     leg->Draw();
1061
1062   if (1) {
1063     // calculate chi^2
1064     TF1* flat = new TF1("flat", "1.0", -TMath::PiOver2(), 3*TMath::PiOver2());
1065     double chi2 = Chi2(hsub, flat);
1066     double chi2_10 = Chi2(hsub10, flat);
1067     int ndof = hsub->GetNbinsX();
1068     delete flat;
1069     c1->cd(1);
1070     ltx.DrawLatex(0.65, 0.05, Form("#chi^{2}/ndf = %.3g / %d", chi2, ndof-1));
1071     if (0)
1072       cout << "Chi2[5] = " << chi2 << " Chi2[10] = " << chi2_10
1073            << " ndof = " << ndof-1 << " (Prob = " << TMath::Prob(chi2,ndof-6)
1074            << " , " << TMath::Prob(chi2_10,ndof-11) << " )" << endl;
1075   }
1076   return c1;
1077 }
1078
1079 double Chi2(TH1* h, TF1* f, TString opt) {
1080   double chi2 = 0;
1081   double ndf = 0;
1082   for (int i=1;i<=h->GetNbinsX();i++) {
1083     double fy = f->Eval(h->GetBinCenter(i));
1084     chi2 += pow( (h->GetBinContent(i) - fy)/h->GetBinError(i) , 2.0);
1085     ndf++;
1086   }
1087   ndf -= 1;
1088
1089   if (opt.Contains("ndf") && ndf > 0)
1090     chi2 /= ndf;
1091
1092   if (opt.Contains("prob") )
1093     return TMath::Prob(chi2, ndf);
1094
1095   return chi2;
1096 }
1097
1098 double vntvna(double *x, double *par)
1099 {
1100   // This function is called many times during a fit. For a 1-d fit,
1101   // x[0] is the x-position of the nth point,
1102   // i.e. graph->GetX()[n]. For a single-valued TGraph, it can be used
1103   // to indicate which point the fitter is currently working on.  The
1104   // par array contains maxpt free parameters that are continually
1105   // updated during the fit process.
1106
1107   // q is the global index
1108   int q = (int)x[0]; // Assume points are at q + 0.5.
1109   int i = 999;
1110   int j = 999;
1111   int nPoints = (maxpt*(maxpt+1))/2;
1112   TF1::fgRejectPoint = kFALSE;
1113
1114   // Assign i and j from q
1115   ijFromGlobalIndex(q, i, j);
1116
1117   if (q < 0 || q >= nPoints) {
1118     Error("vntvna()",
1119           "Problem finding global index. x[0]=%f, q=%d, i=%d, j=%d",
1120           x[0], q, i, j);
1121
1122     return 999;
1123   }
1124
1125   // Exclude points from fit. Warning: TF1::fgRejectPoint is a global
1126   // flag checked by all TF1s, beware of side-effects.
1127   if (j < minFitPtBin || j > maxFitPtBin) {
1128     TF1::RejectPoint();
1129     if (0)
1130       cout << Form("Rejecting x[0],q,i,j = %f %d %d %d",
1131                    x[0], q, i, j)
1132            << endl;
1133     return 0;
1134   }
1135
1136   double v2t = par[i];
1137   double v2a = par[j];
1138   int n = par[maxpt];                       // Fourier coeff. index
1139   int k = par[maxpt+1];                     // cent. bin index
1140   double fitval = v2t * v2a;
1141
1142   // Add extra fit term for momentum conservation
1143   // It is 1/<\sum pt2> * <ptt><pta>
1144   // From ML & JYO, PRL 106, 102301 (2011) eqs 3-7
1145   double mptt=0, mpta=0, v1fac=0;
1146   if (n==1) {
1147     mptt = MeanPt(i, j, k, "t");
1148     mpta = MeanPt(i, j, k, "a");
1149     v1fac = par[maxpt+2];
1150
1151     fitval += v1fac * mptt * mpta;
1152   }
1153
1154   // Odd vnd's < 0 at high pt. // can I replace the below by saying if
1155   // (vndelta < 0) then make fitval < 0?
1156   if (j >= PtBin(5., 6.) && (n%2) )
1157     fitval = -fitval;
1158
1159   return fitval;
1160 }
1161
1162 void VnDelta(int n, const TH1 &h, double &vnd, double &vnd_err, TString opt)
1163 {
1164   // Extract fourier coefficient n directly from h. If opt is "hi" or
1165   // "lo", the coefficient is calculated for the points + or - their
1166   // errors. "hi" and "lo" are meant for an h with systematic (not
1167   // statistical) errors, so vnd_err is not calculated.
1168   double VN = 0.0;
1169   double norm = 0.0;
1170
1171   // Calculate and set coeff. vnd
1172   for (int ib=1; ib<=h.GetNbinsX(); ib++) {
1173     double deltaphi = h.GetBinCenter(ib);
1174     double weight   = h.GetBinContent(ib);
1175
1176     if (opt.Contains("hi"))
1177       weight += h.GetBinError(ib);
1178     else if (opt.Contains("lo"))
1179       weight -= h.GetBinError(ib);
1180
1181     if (opt.Contains("sin") )
1182       VN += sin(n * deltaphi ) * weight;
1183     else
1184       VN += cos(n * deltaphi ) * weight;
1185
1186     norm += weight;
1187   }
1188
1189   if (norm==0)
1190     ::Error("VnDelta()", "Div/0 error");
1191
1192   VN = VN / norm;
1193
1194   // // Mom. cons correction?????????????????????????????????
1195   // if (n==1) VN += 0.002*ptmean[i]*ptmean[j];
1196
1197   vnd = VN;
1198   if (opt.Contains("hi") || opt.Contains("lo")) {
1199     vnd_err = 0;
1200     return;
1201   }
1202
1203   // Statistical uncertainty vnd_err
1204   double quad_sum_uncertainty = 0.0;
1205   for (int ib=1; ib<=h.GetNbinsX(); ib++) {
1206     double deltaphi = h.GetBinCenter(ib);
1207     double dfdwi = cos(n * deltaphi)/norm - VN/norm;
1208
1209     if (opt.Contains("sin") )
1210       dfdwi = sin(n * deltaphi)/norm - VN/norm;
1211
1212     double sigma = h.GetBinError(ib);
1213     quad_sum_uncertainty += dfdwi*dfdwi * sigma*sigma;
1214   }
1215
1216   vnd_err = TMath::Sqrt(quad_sum_uncertainty);
1217   return;
1218 }
1219
1220 void VnDelta(int n, const TF1 &f, double &vnd, double &vnd_err, TString opt)
1221 {
1222   // Extract fourier coefficient n from f.
1223   double VN = 0.0;
1224   double norm = 0.0;
1225   int nSteps = 100;
1226   double x1=0, x2=0;
1227   f.GetRange(x1, x2);
1228   double dx = (x2-x1)/nSteps;
1229
1230   if (!opt.IsNull() ) cout<<opt.Data()<<endl;
1231
1232   // Calculate and set coeff. vnd
1233   for (int ib=0; ib<nSteps; ib++) {
1234
1235     double deltaphi = x1 + ib*dx;
1236     double weight   = f.Eval(deltaphi);
1237
1238     VN += cos(n * deltaphi ) * weight;
1239     norm += weight;
1240   }
1241
1242   if (norm==0)
1243     Error("VnDelta()", "Div/0 error");
1244
1245   VN = VN / norm;
1246   vnd = VN;
1247   vnd_err = 0; // for now
1248   return;
1249 }
1250
1251
1252
1253 // The argument ydiv is the division between pads as a fraction of the
1254 // height.
1255 TCanvas* TwoPanelVert(double ydiv, const char* canvName, const char* canvTitle)
1256 {
1257   TCanvas* c = new TCanvas(canvName, canvTitle, 800, 800);
1258   // Extents of pads inside canvas
1259   double x1=0.01, x2=0.99, y1=0.01, y2=0.99;
1260   double lm=0.2, rm=0.02;
1261
1262   // Division between top and bottom pads
1263   double ysplit = y1 + ydiv*(y2-y1);
1264
1265   // Divide first, adjust afterward.
1266   c->Divide(1,2);
1267   TPad* ptop = (TPad*)(c->GetPad(1));
1268   TPad* pbot = (TPad*)(c->GetPad(2));
1269
1270   ptop->SetPad(x1, ysplit, x2, y2);
1271   ptop->SetMargin(lm, rm, 0.01, 0.01); // left, right, bottom, top.
1272   ptop->SetFrameLineWidth(2);
1273
1274   pbot->SetPad(x1, y1, x2, ysplit);
1275   pbot->SetMargin(lm, rm, 0.4, 0.01); // left, right, bottom, top.
1276   pbot->SetFrameLineWidth(2);
1277   pbot->SetTicks();
1278
1279   c->cd();
1280   return c;
1281 }
1282
1283 char* trigLabel(int i)
1284 {
1285   if (i<0 || i>maxpt)
1286     return Form("Error: no trig pt bin %d", i);
1287
1288   double x1 = gi->GetX()[i] - gi->GetEX()[i];
1289   double x2 = gi->GetX()[i] + gi->GetEX()[i];
1290
1291   // int prec1 = 0, prec2 = 0;
1292   // if (x1 - (int)x1 > 0)
1293   //   prec1 = 1;
1294   // if (x2 - (int)x2 > 0)
1295   //   prec2 = 1;
1296
1297   return Form("%.3g-%.3g", x1, x2);
1298 }
1299
1300 char* asscLabel(int j)
1301 {
1302   if (j<0 || j>maxpt)
1303     return Form("Error: no assc pt bin %d", j);
1304
1305   double x1 = gj->GetX()[j] - gj->GetEX()[j];
1306   double x2 = gj->GetX()[j] + gj->GetEX()[j];
1307
1308   /* Stupid coding...
1309   int prec1 = 0, prec2 = 0;
1310
1311   if (x1 - (int)x1 > 0)
1312     prec1 = 1;
1313   if (x2 - (int)x2 > 0)
1314     prec2 = 1;
1315   */
1316
1317   return Form("%.3g-%.3g", x1, x2);
1318 }
1319
1320 char* centLabel(int k)
1321 {
1322   if (k<0 || k>maxcent)
1323     return Form("Error: no cent bin %d", k);
1324
1325   double c1 = gk->GetX()[k] - gk->GetEX()[k];
1326   double c2 = gk->GetX()[k] + gk->GetEX()[k];
1327
1328   return Form("%.0f-%.0f%%", c1, c2);
1329 }
1330
1331
1332 void VnLimits(double &minVn_k, double &maxVn_k, int k, int imin, int imax)
1333 {
1334   // Find limits for each centrality bin k between imin (inclusive)
1335   // and imax (also inclusive) to plot a common y axis.
1336   //  for (int k=0; k < maxcent; k++) {
1337     for (int n=1; n<=5; n++) {
1338       for (int i=imin; i<=imax; i++) {
1339         TGraphErrors* gc = VnDeltaVsPtAssoc(i, k, n);
1340         if (!gc) {
1341           Error("VnLimits()", "No graph i %d k %d n %d", i,k,n );
1342           continue;
1343         }
1344         for (int j=0; j<gc->GetN(); j++) {
1345           double vn = gc->GetY()[j];
1346           double evn = gc->GetEY()[j];
1347           if (vn+evn > maxVn_k) maxVn_k = vn+evn;
1348           if (vn-evn < minVn_k) minVn_k = vn-evn;
1349         }
1350       }
1351     }
1352
1353   return;
1354 }
1355
1356 TCanvas* DrawQ(int k, int n, TString opt)
1357 {
1358   TStopwatch ts1, ts2;
1359   int fitCurveColor = kRed-4;
1360   int divColor = kBlue;
1361   const char* region = "RIDGE";
1362   if (opt.Contains("ALL"))
1363     region = "ALL";
1364
1365   TString title = Form("globalfit_%d_cent%dto%d",
1366                            n, centlow[k],centhigh[k]);
1367   TString name = Form("V%dvsQ_etamin%02d_cent%dto%d",
1368                       n, (int)(10*minRidgeDEta),
1369                       centlow[k],centhigh[k]);
1370
1371   if (!opt.IsNull()) {
1372     name.Append("_");
1373     name.Append(opt.Data());
1374     title.Append("_");
1375     title.Append(opt.Data());
1376   }
1377
1378   TCanvas* c = TwoPanelVert(0.4, name.Data(), title.Data());
1379   int cwidth = opt.Contains("highptfit") ? 500 : 1200;
1380   c->SetCanvasSize(cwidth, 450);
1381   c->SetWindowSize(cwidth+50, 500);
1382
1383   if (opt.Contains("highptfit")) {
1384     c->GetPad(1)->SetLeftMargin(0.2);
1385     c->GetPad(2)->SetLeftMargin(0.2);
1386     c->GetPad(1)->SetRightMargin(0.05);
1387     c->GetPad(2)->SetRightMargin(0.05);
1388
1389   }
1390
1391   TLegend* lq = new TLegend(0.14, 0.56, 0.25, 0.83);
1392   lq->SetFillColor(kNone);
1393   lq->SetBorderSize(0);
1394
1395   // g:  global VnDelta vs global index (stat errors).
1396   // gs: global VnDelta vs global index (sys errors).
1397   // fn: global fit TF1.
1398   // r:  ratio of g/fit.
1399   TString opt1 = "";
1400   if (opt.Contains("ALL"))
1401     opt1.Append("_ALL");
1402   if (opt.Contains("ptcons"))
1403     opt1.Append("_ptcons");
1404
1405   TGraphErrors* g  = VnDVsQ(n, k, region, cfDef, opt1);
1406   opt1.Append("_sys");
1407   TGraphErrors* gs = VnDVsQ(n, k, region, cfDef, opt1);
1408
1409   int ptaLo = 0, ptaHi = 999;
1410   if (opt.Contains("lowptfit")) {
1411     ptaHi = PtBin(3., 4.);
1412   }
1413   if (opt.Contains("highptfit")) {
1414     ptaLo = PtBin(5., 6.);
1415   }
1416   // If custom pta range was passed in, e.g.
1417   // Form("ptabin%dto%d", PtBin(1.5, 2.0), PtBin(8, 15))
1418   if (opt.Contains("ptabin")) {
1419     TRegexp re("ptabin[0-9]+to[0-9]+");
1420     TRegexp re1("ptabin[0-9]+");
1421     TRegexp re2("to[0-9]+");
1422     TString rx = opt(re), pta1s = rx(re1), pta2s = rx(re2);
1423     pta1s.ReplaceAll("ptabin","");
1424     pta2s.ReplaceAll("to","");
1425     ptaLo = pta1s.Atoi();
1426     ptaHi = pta2s.Atoi();
1427   }
1428
1429   TF1* fn    = GlobalFitVsQ(n, k, ptaLo, ptaHi, region, cfDef, "");
1430   TF1* fn_hi = GlobalFitVsQ(n, k, ptaLo, ptaHi, region, cfDef, "hi_sys");
1431   TF1* fn_lo = GlobalFitVsQ(n, k, ptaLo, ptaHi, region, cfDef, "lo_sys");
1432
1433   // Systematic band on ratio line at unity
1434   TH1F *r_sys = new TH1F(Form("r_sys_%s", title.Data()), "r_sys",
1435                          g->GetN()-0.5, 0, g->GetN()-0.5);
1436   utils.set_hist_props(r_sys, kBlack, kGray, kBlack, kDot, 1.0);
1437
1438   TGraphErrors* r = new TGraphErrors();
1439   TString rname(g->GetName());
1440   rname.Append("_ratio");
1441   r->SetName(rname.Data());
1442
1443   if (!fn) {
1444     Error("DrawQ()","Problem getting/creating global fit function");
1445     gSystem->Exit(-1);
1446   }
1447
1448   // Set visual properties
1449   utils.set_tgraph_props(g, kBlack, kBlack, kFullSquare, 1.0);
1450   utils.set_tgraph_props(gs, kBlack, kBlack, kDot, 0.5);
1451   gs->SetFillColor(kGray);
1452   utils.set_tgraph_props(r, kBlack, kBlack, kFullCircle, 1.0);
1453   g->SetLineWidth(1);
1454   r->SetLineWidth(1);
1455   fn->SetLineWidth(1);
1456   fn->SetLineColor(fitCurveColor);
1457   fn_hi->SetLineStyle(kDashed);
1458   fn_lo->SetLineStyle(kDashed);
1459   fn_hi->SetLineWidth(1);
1460   fn_lo->SetLineWidth(1);
1461
1462   // Compute ratio for r graph
1463   for (int q=0; q<g->GetN(); q++) {
1464     int ii=0, jj=0;
1465     ijFromGlobalIndex(q, ii, jj);
1466
1467     double fy = fn->Eval(q+0.5);
1468     double y = fy==0? 1e9 : g->GetY()[q]/fy;
1469     double rerr = fn_hi->Eval(q+0.5) / fy - 1;
1470
1471     if (opt.Contains("highptfit") && n%2 && g->GetY()[q]!=0)
1472       y = fy/g->GetY()[q];
1473
1474     // No error band when no fit point!
1475     if (jj < ptaLo || jj > ptaHi)
1476       rerr = 0;
1477
1478     r->SetPoint(q, g->GetX()[q], y);
1479     r->SetPointError(q, 0.0, TMath::Abs(g->GetEY()[q]/fy));
1480
1481     r_sys->SetBinContent(q+1, 1.0);
1482     r_sys->SetBinError(q+1, rerr);
1483   }
1484
1485   // Set up plot ranges for top panel (g and f)
1486   int nx = g->GetN();
1487   double x1,y1,x2,y2;
1488   double sf = 1; // scale factor
1489   fn->GetRange(x1,x2);
1490   if (opt.Contains("highptfit")) {
1491     x1 = 45;
1492   }
1493   y1 = sf*fn->GetMinimum(x1,x2);
1494   y2 = sf*fn->GetMaximum(x1,x2);
1495   double marg = 0.4*(y2-y1);
1496   //  if (opt.Contains("lowptfit") || opt.Contains("highptfit") ) marg = 2*(y2-y1);
1497   x1 -= 0.5;
1498   x2 += 0.5;
1499   y1 -= marg;
1500   y2 += marg;
1501
1502   // Set up plot ranges for lower panel (r)
1503   double ry1=0.21, ry2=1.79;
1504   if (opt.Contains("highptfit")) {
1505     ry1 = 0.01;
1506     ry2 = 1.99;
1507   }
1508
1509   // Get / make frame histos for top and bottom
1510   TH2F* hf1 = 0; // (TH2F*)gDirectory->Get("qhf1");
1511   TH2F* hf2 = 0; // (TH2F*)gDirectory->Get("qhf2");
1512   TString name1 = Form("qhf1_%s", title.Data());
1513   TString name2 = Form("qhf2_%s", title.Data());
1514   // TString sv = n==1 ? "#LTcos(#Delta#phi)#GT" :
1515   // Form("#LTcos(%d#Delta#phi)#GT", n);
1516   TString sv = Form("V_{%d#Delta}", n);
1517   TString sgf = Form("(v_{%d}^{t}v_{%d}^{a})_{fit}", n, n);
1518
1519   if (!hf1)
1520     hf1 = new TH2F(name1.Data(), "hf1", nx+1, x1, x2, 1000, y1, y2);
1521   if (!hf2)
1522     hf2 = new TH2F(name2.Data(), "hf2", nx+1, x1, x2, 1000, ry1, ry2);
1523
1524   // Axes
1525   TAxis* ax1 = hf1->GetXaxis();
1526   TAxis* ay1 = hf1->GetYaxis();
1527   TAxis* ax2 = hf2->GetXaxis();
1528   TAxis* ay2 = hf2->GetYaxis();
1529   hf1->SetBins(nx+1, x1, x2, 100, y1, y2);
1530   hf2->SetBins(nx+1, x1, x2, 100, ry1, ry2);
1531   hf1->SetAxisRange(x1, x2, "X");
1532   hf1->SetAxisRange(y1, y2, "Y");
1533   hf2->SetAxisRange(x1, x2, "X");
1534   hf2->SetAxisRange(ry1, ry2, "Y");
1535   ax1->SetLabelOffset(0.02);
1536   ax1->SetLabelSize(0.04);  ay1->SetLabelSize(0.08);
1537   ax2->SetLabelSize(0.18);   ay2->SetLabelSize(0.12);
1538   ax1->SetNdivisions(1);   ax2->SetNdivisions(1);
1539   ay1->SetNdivisions(208);  ay2->SetNdivisions(opt.Contains("highptfit")?104:210);
1540   ax1->SetTickLength(0.02);
1541
1542   //ay1->SetTitle(Form("%s and #color[%d]{%s}",sv.Data(), fitCurveColor, sgf.Data()));
1543   ay1->SetTitle(Form("%s", sv.Data()));
1544   ay2->SetTitle(Form("#frac{%s}{fit}", sv.Data()));
1545   //  ay2->SetTitle(Form("#frac{%s}{%s}", sv.Data(), sgf.Data()));
1546   // ax2->SetTitle(Form("#splitline{associated p_{T}^{a}}{#color[%d]{   trigger p_{T}^{t}}} [GeV/c]", divColor));
1547   ax2->SetTitle(Form("#color[%d]{p_{T}^{t}}, p_{T}^{a} [GeV/c]", divColor));
1548   ay1->SetTitleOffset(opt.Contains("highptfit") ? 1. : 0.45);
1549   ay1->SetTitleSize(0.09);
1550   ay1->CenterTitle();
1551   ay2->SetTitleOffset(opt.Contains("highptfit") ? 0.6 : 0.26);
1552   ay2->SetTitleSize(0.14);
1553   ay2->CenterTitle();
1554   ax2->SetLabelOffset(2.1);
1555   ax2->SetTitleOffset(1.3);
1556   ax2->SetTitleSize(0.14);
1557   ax2->CenterTitle();
1558
1559   double leftMarg = opt.Contains("highptfit") ? 0.2 : 0.08;
1560   c->cd(1);
1561   gPad->SetFrameLineWidth(1);
1562   gPad->SetLeftMargin(leftMarg);
1563   hf1->Draw();
1564   TLine zero;
1565   zero.DrawLine(x1, 0., x2, 0.);
1566
1567
1568   // NEW========================================
1569   for (int i=0; i<gi->GetN(); i++) {
1570     int q1 = GlobalIndex(i, 1), q2 = GlobalIndex(i, i);
1571     TF1* ifn = 0, *ifn_hi = 0, *ifn_lo = 0;
1572     //    TF1* _ifn = 0, *_ifn_hi = 0, *_ifn_lo = 0; // temp, unscaled.
1573
1574     // This works, but has lines between \ptt groups :(
1575     // TH1F *hfn=0, *hfn_hi=0, *hfn_lo=0;
1576     // int nq = g->GetN();
1577     // hfn = new TH1F(Form("%s_%d", fn->GetName(), i), Form("%s_%d", fn->GetName(), i),
1578     //             nq, 0, nq);
1579     // hfn->SetLineColor(kMagenta);
1580     // for (int q=0; q<g->GetN(); q++) {
1581     //   hfn->SetBinContent(q+1, sf*fn->Eval(q+0.5));
1582     // }
1583
1584     ifn = (TF1*)       fn->Clone(Form("%s_%d", fn->GetName(), i));
1585     ifn_hi = (TF1*) fn_hi->Clone(Form("%s_%d", fn_hi->GetName(), i));
1586     ifn_lo = (TF1*) fn_lo->Clone(Form("%s_%d", fn_lo->GetName(), i));
1587
1588     ifn->SetName(Form("%s_%d", fn->GetName(), i));
1589     ifn_hi->SetName(Form("%s_%d", fn_hi->GetName(), i));
1590     ifn_lo->SetName(Form("%s_%d", fn_lo->GetName(), i));
1591
1592     double r1 = q1-0.9, r2 = q2+0.92;
1593     ifn->SetRange(r1, r2);
1594     ifn_hi->SetRange(r1, r2);
1595     ifn_lo->SetRange(r1, r2);
1596
1597
1598     // for (int ip=0; ip<ifn->GetNpar(); ip++) {
1599     //   double par = ifn->GetParameter(ip);
1600     //   ifn->SetParameter(ip, 1000*par);
1601     //   cout<<par<< " " << ifn->GetParameter(ip) << endl;
1602     // }
1603
1604
1605     // ifn    = new TF1(Form("%s_%d",    fn->GetName(), i), Form("100*%s",   _ifn->GetName()));
1606     // ifn_hi = new TF1(Form("%s_%d", fn_hi->GetName(), i), Form("100*%s", _ifn_hi->GetName()));
1607     // ifn_lo = new TF1(Form("%s_%d", fn_lo->GetName(), i), Form("100*%s", _ifn_lo->GetName()));
1608
1609     ifn_hi->Draw("same");
1610     ifn_lo->Draw("same");
1611     ifn->Draw("same");
1612     //    hfn->Draw("same");
1613   }
1614   // NEW========================================
1615
1616   for (int i=0; i<gi->GetN(); i++) {
1617     char* newName = Form("%s_ptt%d", g->GetName(), i);
1618     TGraphErrors* gcol = (TGraphErrors*)g->Clone(newName);
1619     int col = PtColor(ptlow[i], pthigh[i]);
1620     utils.set_tgraph_props(gcol, kBlack, col, kFullSquare, 1.0);
1621     gcol->SetLineWidth(1);
1622
1623     newName = Form("%s_ptt%d", gs->GetName(), i);
1624     TGraphErrors* g_sy = (TGraphErrors*)gs->Clone(newName);
1625     int col_sy = PtColorPale(ptlow[i], pthigh[i]);
1626     utils.set_tgraph_props(g_sy, kBlack, col_sy, kDot, 0.5);
1627     g_sy->SetFillColor(col_sy);
1628
1629     for (int q=0; q<gcol->GetN(); q++) {
1630       int i2=0, j=0;
1631       ijFromGlobalIndex(q,i2,j);
1632       if (i2!=i) {
1633         gcol->SetPoint(q, gcol->GetX()[q], -99.);
1634         g_sy->SetPoint(q, gcol->GetX()[q], -99.);
1635       }
1636
1637       // Scale points by sf
1638       gcol->SetPoint(q, gcol->GetX()[q], sf*gcol->GetY()[q]);
1639       g_sy->SetPoint(q, g_sy->GetX()[q], sf*g_sy->GetY()[q]);
1640       gcol->SetPointError(q, gcol->GetEX()[q], sf*gcol->GetEY()[q]);
1641       g_sy->SetPointError(q, g_sy->GetEX()[q], sf*g_sy->GetEY()[q]);
1642
1643     }
1644     g_sy->Draw("e2psame");
1645     gcol->Draw("ezpsame");
1646   }
1647
1648   c->cd(2);
1649   gPad->SetFrameLineWidth(1);
1650   gPad->SetLeftMargin(leftMarg);
1651   hf2->Draw();
1652   r_sys->Draw("e2psame");
1653
1654   // Overlay lines, arrows, text, etc.
1655   TLine div, div_a; // division ticks for trigger and assoc
1656   TArrow a;
1657   double tipSize = 0.008; // size of arrowhead
1658   a.SetLineWidth(1);
1659   div.SetLineWidth(2);
1660   div.SetLineColor(fitCurveColor);
1661   div.DrawLine(x1, 1.0, x2, 1.0);
1662   r->Draw("ezpsame");
1663   gList->Add(r);
1664
1665   for (int i=0; i<gi->GetN(); i++) {
1666     char* newName = Form("%s_ptt%d", r->GetName(), i);
1667     TGraphErrors* rcol = (TGraphErrors*)r->Clone(newName);
1668     int col = PtColor(ptlow[i], pthigh[i]);
1669     utils.set_tgraph_props(rcol, kBlack, col, kFullCircle, 1.0);
1670     for (int q=0; q<rcol->GetN(); q++) {
1671       int i2=0, j=0;
1672       ijFromGlobalIndex(q,i2,j);
1673       if (i2!=i) rcol->SetPoint(q, rcol->GetX()[q], -99.);
1674     }
1675     rcol->SetLineWidth(1);
1676     rcol->Draw("ezpsame");
1677   }
1678
1679   div.SetLineColor(divColor);
1680   div.SetLineStyle(kSolid);
1681
1682   // Label centrality and fourier index
1683   c->cd(1);
1684   ltx.SetTextSize(0.14);
1685   ltx.DrawLatex(opt.Contains("highptfit")?0.25:0.14, 0.85, Form("n = %d", n));
1686
1687   if (!opt.Contains("highptfit")) {
1688   if (n<=2) {
1689     ltx.SetTextSize(0.1);
1690     if (ispp)
1691     ltx.DrawLatex(0.14, 0.70, Form("pp 2.76 TeV"));
1692     else {
1693       ltx.DrawLatex(0.14, 0.70, Form("Pb-Pb 2.76 TeV"));
1694       ltx.DrawLatex(0.14, 0.61, Form("%d-%d%%", centlow[k],centhigh[k]));
1695     }
1696     if (opt.Contains("ptcons")) {
1697     ltx.DrawLatex(0.14, 0.58, Form("#LT#sum p_{T}^{2}#GT = 4.76 GeV^{2}"));
1698     }
1699
1700   }
1701   if (n==3) {
1702     lq->AddEntry(g, Form("V_{n#Delta}"), "epl");
1703     lq->AddEntry(fn, Form("fit"), "l");
1704     lq->Draw();
1705   }
1706   // Always show
1707     ltx.SetTextSize(0.07);
1708   double ptaFitLow = ptlow[ptaLo];
1709   double ptaFitHigh = ptaHi<999? pthigh[ptaHi] : pthigh[maxpt-1];
1710   //  if (!opt.Contains("highptfit"))
1711     ltx.DrawLatex(0.14, 0.09,
1712                   Form("Global fit range: %3.2g < p_{T}^{a} < min(p_{T}^{t}, %3.2g GeV/c)",
1713                        ptaFitLow, ptaFitHigh));
1714   }
1715
1716   for (int q=0; q<g->GetN(); q++) {
1717     int i=0, j=0;
1718     ijFromGlobalIndex(q,i,j);
1719     double drop = 0.35;
1720     double ptt = gi->GetX()[i] + gi->GetEX()[i];
1721     double pta = gj->GetX()[j] + gj->GetEX()[j];
1722     double tickLength = 0.03*(y2-y1);
1723
1724     // Draw divisions between trigger pt intervals...
1725     bool skipTicks = opt.Contains("highptfit") && (pta < 5.);
1726     if (i==j && !skipTicks) {
1727       //      if (ptt > 5) {
1728         c->cd(1);
1729         div.DrawLine(q+1, y1, q+1, y1+tickLength);
1730         c->cd(2);
1731         tickLength = drop-0.1;
1732         div.DrawLine(q+1, ry2-tickLength, q+1, ry2+tickLength);
1733         div.DrawLine(q+1, ry1-tickLength, q+1, ry1+tickLength);
1734         //      }
1735
1736       // Alignment: align = 10*HorizontalAlign + VerticalAlign
1737       // For horizontal alignment the following convention applies:
1738       // 1=left adjusted, 2=centered, 3=right adjusted
1739       // For vertical alignment the following convention applies:
1740       // 1=bottom adjusted, 2=centered, 3=top adjusted
1741       ltx2.SetTextAlign(23);
1742       ltx2.SetTextSize(0.12);
1743       ltx2.SetTextColor(divColor);
1744
1745       if (ptt != 0.75) // too cluttered otherwise
1746         ltx2.DrawLatex(q+1, ry1-drop, Form("%.2g", ptt) ); // pt_trig labels
1747       }
1748
1749     if (pta < ptt && (pta == (int)pta || pta == 0.5) && !skipTicks) {
1750         div_a.DrawLine(q+1, ry1-drop/8, q+1, ry1+drop/4);
1751         ltx2.SetTextColor(kBlack);
1752         ltx2.SetTextFont(52); // italic
1753         ltx2.DrawLatex(q+1, ry1-drop/4, Form("%.2g", pta) ); // pt_assoc labels
1754         ltx2.SetTextColor(divColor);
1755         ltx2.SetTextFont(62); // 62 is default
1756     }
1757
1758     // Draw arrows to off-scale points
1759     c->cd(1);
1760     double y = g->GetY()[q];
1761     double ymax = hf1->GetYaxis()->GetXmax();
1762     double ymin = hf1->GetYaxis()->GetXmin();
1763     double len = 0.06*(ymax-ymin);
1764     double gap = 0.01*(ymax-ymin);
1765
1766     if (y > ymax)
1767       a.DrawArrow(q+0.5, ymax-len-gap, q+0.5, ymax-gap, tipSize, ">");
1768     if (y < ymin)
1769       a.DrawArrow(q+0.5, ymin+len+gap, q+0.5, ymin+gap, tipSize, ">");
1770
1771     // Do the same for the ratio.
1772     // If fit has not been applied to some points, don't draw the arrow.
1773     c->cd(2);
1774     //    gPad->SetTickx();
1775     y = r->GetY()[q];
1776     ymax = hf2->GetYaxis()->GetXmax();
1777     ymin = hf2->GetYaxis()->GetXmin();
1778     len = 0.06/0.4*(ymax-ymin); // 0.4 from TwoPanelVert ctor
1779     gap = 0.02*(ymax-ymin);
1780
1781     if (0) { // Draw at upper or lower edge of panel
1782       if (y > ymax)
1783         a.DrawArrow(q+0.5, ymax-len-gap, q+0.5, ymax-gap, tipSize, ">");
1784       if (y < ymin)
1785         a.DrawArrow(q+0.5, ymin+len+gap, q+0.5, ymin+gap, tipSize, ">");
1786     }
1787     // Draw with arrow based from ref. line.
1788     // No arrow when no fit point!
1789     if (j < ptaLo || j > ptaHi)
1790       continue;
1791     if (y > ymax)
1792       a.DrawArrow(q+0.5, 1.0, q+0.5, 1.0+len, tipSize, ">");
1793     if (y < ymin)
1794       a.DrawArrow(q+0.5, 1.0, q+0.5, 1.0-len, tipSize, ">");
1795   } // end q loop
1796
1797   //  ax2->LabelsOption("h");
1798   c->cd();
1799   c->Update();
1800   return c;
1801 }
1802
1803 TCanvas* DrawVn(int k, int ptt1, int ptt2)
1804 {
1805   // Draw Fourier coefficients on one canvas for each centrality bin k
1806   // for trigger pt bins from ptt1 to (and including) ptt2.
1807   // The x-axis is associated pt.
1808   const char* name = Form("vndelta_etamin%02d_cent%dto%d_trig%dto%d",
1809                           (int)(10*minRidgeDEta), centlow[k],centhigh[k], ptt1, ptt2);
1810   TCanvas* cv = new TCanvas(name, name, 1200, 500);
1811   TLegend* lv = new TLegend(0.09, 0.72, 0.27, 0.88); // w.r.t. cv
1812   lv->SetFillColor(kNone);
1813   lv->SetBorderSize(0);
1814
1815   double lv2bottom = ptt2<5? 0.65 : 0.8;
1816   TLegend* lv2 = new TLegend(0.82, lv2bottom, 0.94, 0.95, "p_{T}^{t} (GeV/c)");
1817   lv2->SetFillColor(kNone);
1818   lv2->SetBorderSize(0);
1819   utils.padsetup(cv, 5, 1, "", 0.12, 0.01, 0.03, 0.2);
1820
1821   // Find limits
1822   double minVn_k=1000, maxVn_k=-1000;
1823   VnLimits(minVn_k, maxVn_k, k, ptt1, ptt2);
1824
1825     for (int n=1; n<=5; n++) {
1826       cv->cd(n);
1827       double vertMargin = 0.1*(maxVn_k - minVn_k);
1828       double xmax = gj->GetX()[ptt2] + gj->GetEX()[ptt2];
1829       double xmin =  -0.43;
1830       TH1F* hf = gPad->DrawFrame(xmin, minVn_k - vertMargin,
1831                                  xmax, maxVn_k + vertMargin);
1832       int nXticks = 404;
1833       int nYticks = 210;
1834       if (n==1) {
1835         utils.make_nice_axes(cv, hf, 2., nXticks, nYticks, -0.02, 0.02);
1836         hf->SetLabelOffset(-0.01, "X");
1837       }
1838       else {
1839         utils.make_nice_axes(cv, hf, 3.4, nXticks, nYticks, -0.075, 0.0);
1840         hf->SetLabelOffset(-0.057, "X");
1841         hf->SetTickLength(0.02, "X");
1842         hf->SetTickLength(0.05, "Y");
1843       }
1844       hf->GetXaxis()->SetTicks("+-");
1845       hf->GetYaxis()->SetTicks("+-");
1846
1847       // Loop backwards to improve visibility of low-pt points
1848       for (int i=ptt2; i>=ptt1; i--) {
1849
1850         TGraphErrors *gc=0;
1851         gc = VnDeltaVsPtAssoc(i, k, n);
1852         if (!gc) {
1853           Error("DrawVn()", "Problem finding graph n%d, i%d, k%d" ,n,i,k);
1854           continue;
1855         }
1856
1857         utils.set_tgraph_props(gc, colorsDark[i], colorsDark[i], kFullSquare, 1.6);
1858
1859         TGraphErrors* gfit = GlobalFitVsPtAssoc(i, k, n);
1860
1861         gfit->Draw("plsame");
1862         gc->Draw("epsame");
1863
1864         if (n==1)
1865           lv2->AddEntry(gc,trigLabel(i), "p");
1866
1867       }
1868     }
1869
1870     cv->cd();
1871     TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1);
1872     overlay->SetFillStyle(4000); // transparent
1873     overlay->Draw();
1874     overlay->cd();
1875
1876     ltx.SetTextSize(0.07);
1877     lv2->Draw();
1878     ltx.DrawLatex(0.14, 0.76, Form("%s", centLabel(k)));
1879
1880     for (int n=1; n<=4; n++) {
1881       cv->cd(n);
1882       ltx.SetTextSize(n==1? 0.1 : 0.15);
1883       ltx.DrawLatex((n==1)? 0.45 : 0.1, 0.86, Form("V_{%d#Delta}", n));
1884     }
1885     cv->cd(5);
1886     ltx.DrawLatex(0.72, 0.86, Form("V_{%d#Delta}", 5));
1887     overlay->cd();
1888     ltx.SetTextSize(0.075);
1889     ltx.DrawLatex(0.4, 0.05, "associated p_{T} [GeV/c]");
1890     ltx.SetTextAngle(90);
1891     ltx.DrawLatex(0.03, 0.4, "V_{n#Delta}(p_{T}^{assoc.})");
1892     ltx.SetTextAngle(0);
1893
1894   return cv;
1895 }
1896
1897 TCanvas* DrawVnVsPtTrig(int k, int npta, int ptabins[], TString opt)
1898 {
1899   TString title = Form("VnDeltaVsPtTrig_cent%dto%d_pTassoc%.2gto%.2g",
1900                        centlow[k],centhigh[k],
1901                        ptlow[ptabins[0]], pthigh[ptabins[npta-1]]);
1902   if (!opt.IsNull()) {
1903     title.Append("_");
1904     title.Append(opt);
1905   }
1906   TCanvas* cv = new TCanvas(title.Data(), title.Data(), 1400, 500);
1907   double lv2bot = 0.9 - 0.07*npta;
1908   TLegend* lv2 = new TLegend(0.82, lv2bot, 0.98, 0.95, "p_{T}^{a} (GeV/c)");
1909   lv2->SetFillColor(kNone);
1910   lv2->SetBorderSize(0);
1911   utils.padsetup(cv, 5, 1, "", 0.12, 0.01, 0.03, 0.2);
1912
1913   // Start without knowing y-limits, adjust below.
1914   double minVn_k=-1., maxVn_k=1.;
1915   TObjArray* graphs = new TObjArray();
1916
1917   for (int n=1; n<=5; n++) {
1918     cv->cd(n);
1919     double vertMargin = 0.1*(maxVn_k - minVn_k);
1920     double xmax = 11.1;
1921     double xmin = -0.43;
1922
1923     TH1F* hf = gPad->DrawFrame(xmin, minVn_k - 4*vertMargin,
1924                                xmax, maxVn_k + vertMargin);
1925     TLine l;
1926     l.DrawLine(xmin, 0, xmax, 0);
1927
1928     int nXticks = 404;
1929     int nYticks = 210;
1930     if (n==1) {
1931       utils.make_nice_axes(cv, hf, 2., nXticks, nYticks, -0.02, 0.02);
1932       hf->SetLabelOffset(-0.01, "X");
1933     }
1934     else {
1935       utils.make_nice_axes(cv, hf, 3.4, nXticks, nYticks, -0.075, 0.0);
1936       hf->SetLabelOffset(-0.057, "X");
1937       hf->SetTickLength(0.02, "X");
1938       hf->SetTickLength(0.05, "Y");
1939     }
1940     hf->GetXaxis()->SetTicks("+-");
1941     hf->GetYaxis()->SetTicks("+-");
1942
1943     for (int jpta=0; jpta<npta; jpta++) {
1944       int j = ptabins[jpta];
1945       TGraphErrors *gc = VnDeltaVsPtTrig(j, k, n);
1946       if (!gc) {
1947         Error("DrawVnVsPtTrig()", "Problem finding gc n%d, i%d, k%d" ,n,j,k);
1948         continue;
1949       }
1950       TGraphErrors *gs = VnDeltaVsPtTrig(j, k, n, "sys");
1951       if (!gs) {
1952         Error("DrawVnVsPtTrig()", "Problem finding gs n%d, i%d, k%d" ,n,j,k);
1953         continue;
1954       }
1955
1956       int col = PtColor(ptlow[j], pthigh[j]);
1957       utils.set_tgraph_props(gc, col, col, kFullSquare, 1.2);
1958       gs->SetFillColor(PtColorPale(ptlow[j], pthigh[j]));
1959       TGraphErrors* open = (TGraphErrors*)gc->Clone();
1960       utils.set_tgraph_props(open, kBlack, kBlack, kOpenSquare, 1.2);
1961
1962       gs->Draw("e2psame");
1963       gc->Draw("epsame");
1964       open->Draw("epsame");
1965
1966       if (n==1)
1967         lv2->AddEntry(gc,asscLabel(j), "elp");
1968
1969       graphs->Add(gc); // store graph array to find y limits for plot
1970     }
1971   }
1972
1973   // Adjust y limits
1974   double newYmin = 0, newYmax = 0, newSpace = 0;
1975   for (int m=0; m<graphs->GetEntries(); m++) {
1976     TGraphErrors* tg = (TGraphErrors*)graphs->At(m);
1977     double yhi = TMath::MaxElement(tg->GetN()-1, tg->GetY());
1978     double ylo = TMath::MinElement(tg->GetN()-1, tg->GetY());
1979     double ey = TMath::MaxElement(tg->GetN()-1, tg->GetEY());
1980     ylo -= ey;
1981     if (yhi > newYmax) newYmax = yhi;
1982     if (ylo < newYmin) newYmin = ylo;
1983   }
1984   newSpace = 0.1*(newYmax - newYmin);
1985   for (int n=1; n<=5; n++) {
1986     cv->cd(n);
1987     TH1F* hf = (TH1F*)gPad->FindObject("hframe");
1988     hf->GetYaxis()->SetRangeUser(newYmin-3*newSpace, newYmax+newSpace);
1989     gPad->Modified(); // Signal to redraw
1990   }
1991
1992   cv->cd();
1993   TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1);
1994   overlay->SetFillStyle(4000); // transparent
1995   overlay->Draw();
1996   overlay->cd();
1997   lv2->Draw();
1998   ltx.SetTextSize(0.05);
1999   ltx.DrawLatex(0.14, 0.88, Form("%s 2.76 TeV", ispp?"pp":"Pb-Pb"));
2000   ltx.SetTextSize(0.07);
2001   if (!ispp)
2002     ltx.DrawLatex(0.14, 0.8, Form("%s", centLabel(k)));
2003
2004   for (int n=1; n<=5; n++) {
2005     cv->cd(n);
2006     ltx.SetTextSize(n==1? 0.1 : 0.15);
2007     if (n==1)
2008       ltx.DrawLatex(0.44, 0.24, Form("n=%d", n));
2009     else
2010       ltx.DrawLatex(0.04, 0.24, Form("n=%d", n));
2011
2012   }
2013   overlay->cd();
2014   ltx.SetTextSize(0.075);
2015   ltx.DrawLatex(0.48, 0.05, "trigger p_{T}^{t} [GeV/c]");
2016   ltx.SetTextAngle(90);
2017   ltx.DrawLatex(0.06, 0.48, "V_{n#Delta}  [10^{-2}]");
2018   ltx.SetTextAngle(0);
2019   return cv;
2020 }
2021
2022 TGraphErrors* AgreementVsPtSum(int k, int n)
2023 {
2024   TGraphErrors* gVn = VnDVsQ(n, k);
2025   TF1* fVn = GlobalFitVsQ(n,k);
2026   TGraphErrors* g = new TGraphErrors();
2027   int N = gVn->GetN(); // And f better have the same N.
2028   int i = 0, j = 0;
2029   if(N<=0) {
2030     Error("AgreementVsPtSum()", "Invalid number of points: %d", N);
2031     gSystem->Exit(1);
2032   }
2033
2034   for (int q=0; q<N; q++) {
2035     ijFromGlobalIndex(q,i,j);
2036     double gy = gVn->GetY()[q];
2037     double fy = fVn->Eval(q+0.5);
2038     double r = TMath::Abs((gy-fy)/gy);
2039     g->SetPoint(q, PtBinCenterSum(i, j), r);
2040   }
2041   g->SetMarkerStyle(kFullCircle);
2042   return g;
2043 }
2044
2045 TGraph2D* Agreement2D(int k, int n)
2046 {
2047   TGraphErrors* gVn = VnDVsQ(n, k);
2048   TF1* fVn = GlobalFitVsQ(n,k);
2049   TGraph2D* g = new TGraph2D();
2050   int N = gVn->GetN(); // And f better have the same N.
2051   int i = 0, j = 0;
2052   if(N<=0) {
2053     Error("AgreementVsPtSum()", "Invalid number of points: %d", N);
2054     gSystem->Exit(1);
2055   }
2056
2057   for (int q=0; q<N; q++) {
2058     ijFromGlobalIndex(q,i,j);
2059     double gy = gVn->GetY()[q];
2060     double fy = fVn->Eval(q+0.5);
2061     double r = TMath::Abs((gy-fy)/fy);
2062     g->SetPoint(q, ptmean[j], ptmean[i], r);
2063   }
2064   return g;
2065 }
2066
2067 TH2F* Agreement2DHist(int k, int n)
2068 {
2069   TGraphErrors* gVn = VnDVsQ(n, k);
2070   TF1* fVn = GlobalFitVsQ(n,k);
2071
2072   // Pt binning array
2073   double ptarr[100];
2074   for (int i=0; i<maxpt; i++) ptarr[i] = gi->GetX()[i] - gi->GetEX()[i];
2075   ptarr[maxpt] = gi->GetX()[maxpt-1] + gi->GetEX()[maxpt-1];
2076
2077   int N = gVn->GetN(); // And f better have the same N.
2078   int i = 0, j = 0;
2079   if(N<=0) {
2080     Error("AgreementVsPtSum()", "Invalid number of points: %d", N);
2081     gSystem->Exit(1);
2082   }
2083   const char* name = Form("agmt_cent%dto%d_n%d",
2084                           centlow[k],centhigh[k], n);
2085   TH2F* h = new TH2F(name, name, maxpt, ptarr, maxpt, ptarr);
2086   h->SetTitle("title;p_{T}^{a};p_{T}^{t};");
2087
2088   for (int q=0; q<N; q++) {
2089     ijFromGlobalIndex(q,i,j);
2090     double gy = gVn->GetY()[q];
2091     double fy = fVn->Eval(q+0.5);
2092
2093     // double r = fy==0? 0 : TMath::Abs(gy-fy) / fy;
2094     double r = gy==0? 0 : TMath::Abs((gy-fy) / gy);
2095     //    double r = fy==0? 0 : TMath::Abs(gy / fy);
2096     h->SetBinContent(j+1, i+1, r);
2097   }
2098   return h;
2099 }
2100
2101 TH2* PttPtaHist()
2102 {
2103   // Pt binning array
2104   double ptarr[100];
2105   for (int i=0; i<maxpt; i++) ptarr[i] = gi->GetX()[i] - gi->GetEX()[i];
2106   ptarr[maxpt] = gi->GetX()[maxpt-1] + gi->GetEX()[maxpt-1];
2107
2108   int i = 0, j = 0;
2109   const char* name = Form("PttPtaHist");
2110   TH2F* h = new TH2F(name, name, maxpt, ptarr, maxpt, ptarr);
2111   h->SetTitle("title;p_{T}^{a} (GeV/c);p_{T}^{t} (GeV/c);");
2112   h->GetXaxis()->CenterTitle();
2113   h->GetYaxis()->CenterTitle();
2114   int N = maxpt * (maxpt + 1) / 2;
2115   for (int q=0; q<N; q++) {
2116     ijFromGlobalIndex(q,i,j);
2117     h->SetBinContent(j+1, i+1, 1.);
2118   }
2119   return h;
2120 }
2121
2122 TGraphErrors* Chi2VsTrigPt(int k, int n)
2123 {
2124   TGraphErrors* g = new TGraphErrors();
2125   for (int i=0; i<maxpt; i++) {
2126     g->SetPoint(i, ptmean[i], ReducedChi2(i, k, n));
2127   }
2128   g->SetMarkerStyle(kFullCircle);
2129   return g;
2130 }
2131
2132 TGraphErrors*
2133 CorrFnChi2VsTrigPt(int j, int k, int n1, int n2, TString opt)
2134 {
2135   TGraphErrors* g = new TGraphErrors();
2136   int ip=0;
2137   for (int i=j; i<maxpt; i++) {
2138     TH1* h = Hist("RIDGE", "cA", i, j, k);
2139     double x = MeanPt(i, j, k, "t");
2140     double y = Chi2(h, HarmonicSum(h, n1, n2, "global"), opt);
2141     g->SetPoint(ip, x, y);
2142     ip++;
2143   }
2144   g->SetMarkerStyle(kFullCircle);
2145   return g;
2146 }
2147
2148 double ReducedChi2(int i, int k, int n)
2149 {
2150   TGraphErrors* g = VnDeltaVsPtAssoc(i, k, n);
2151   TGraphErrors* f = GlobalFitVsPtAssoc(i, k, n);
2152   int N = g->GetN(); // And f better have the same N.
2153   double chi2 = 0;
2154
2155   // Value (or one-sigma error if opt=="err") of v_n
2156   //  double vnerr = Bestvn(k,n,i,"err");
2157   double vnerr = vnGF(n,k,i,0,999,"RIDGE","cA","err");
2158
2159   if(N<=0 || N!= f->GetN()) {
2160     Error("ReducedChi2()", "Invalid number of points: %d", N);
2161     gSystem->Exit(1);
2162   }
2163
2164   for (int j=0; j<N; j++) {
2165     double errj = g->GetEY()[j];
2166
2167     chi2 +=
2168       (g->GetY()[j] - f->GetY()[j]) * (g->GetY()[j] - f->GetY()[j]) /
2169       pow(errj + vnerr, 2);
2170   }
2171   chi2 /= N;
2172   return chi2;
2173 }
2174
2175 TGraphErrors* VnDeltaVsPtAssoc(int i, int k, int n, TString opt)
2176 {
2177   // Create a TGraph from VnDelta points
2178   TGraphErrors* gALL = VnDVsQ(n, k, "RIDGE", "cA", opt);
2179   TGraphErrors* g = new TGraphErrors();
2180   TString gname(Form("VnDeltaVsPtAssoc_ptt%d_cent%d_n%d",i,k,n));
2181   gname.Append(opt);
2182   g->SetName(gname.Data());
2183
2184   for (int j=0; j<=i; j++) {
2185     int q = GlobalIndex(i, j);
2186     g->SetPoint(j, gj->GetX()[j], gALL->GetY()[q] );
2187     g->SetPointError(j, 0,  gALL->GetEY()[q] );
2188     utils.set_tgraph_props(g, colorsDark[i], colorsDark[i], kFullSquare, 1.0);
2189     g->SetLineWidth(2);
2190   }
2191
2192   return g;
2193 }
2194
2195 TGraphErrors* VnDeltaVsPtTrig(int j, int k, int n, TString opt)
2196 {
2197   // Create a TGraph from VnDelta points
2198   TGraphErrors* gALL = VnDVsQ(n, k, "RIDGE", "cA", opt);
2199   TGraphErrors* g = new TGraphErrors();
2200   TString gname(Form("VnDeltaVsPtTrig_pta%d_cent%dto%d_n%d",
2201                      j,centlow[k],centhigh[k],n));
2202   gname.Append(opt);
2203   g->SetName(gname.Data());
2204
2205   int ip=0;
2206   for (int i=j; i<maxpt; i++) {
2207     int q = GlobalIndex(i, j);
2208
2209     double xerr = 0.;
2210     if (opt.Contains("sys")) {
2211       if (gi->GetEX()[i] < 0.126)
2212         xerr = gi->GetEX()[i];
2213       else
2214         xerr = 0.250;
2215     }
2216     double mpt = MeanPt(i, j, k, "t");
2217     double sf = 100; // Scale y-value; make sure to label accordingly!!
2218     g->SetPoint(ip, mpt, sf*gALL->GetY()[q] );
2219     g->SetPointError(ip, xerr,  sf*gALL->GetEY()[q] );
2220
2221     utils.set_tgraph_props(g, colorsDark[j], colorsDark[j], kFullSquare, 1.0);
2222     g->SetLineWidth(2);
2223     ip++;
2224   }
2225
2226   return g;
2227 }
2228
2229 TGraphErrors* VnDeltaVsN(int i, int j, int k, int nMax, TString opt)
2230 {
2231   // Create a TGraph from VnDelta points
2232   // Stat errors by default. Call with opt "sys" to get sys errors.
2233   // Call with opt "sine" to get residual <sin n Delta phi> for sys. est.
2234   int q = GlobalIndex(i, j);
2235   TString gname(Form("VnDeltaVsN_ptt%.2gto%.2g_pta%.2gto%.2g_cent%dto%d",
2236                      ptlow[i],pthigh[i],ptlow[j],pthigh[j],
2237                      centlow[k],centhigh[k]));
2238   if (!opt.IsNull())  {
2239     gname.Append("_");
2240     gname.Append(opt);
2241   }
2242   TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data());
2243   if (g) {
2244     return g;
2245   }
2246   else {
2247     g = new TGraphErrors();
2248     g->SetName(gname.Data());
2249   }
2250
2251   int ip=0;
2252   gStyle->SetEndErrorSize(5);
2253   for (int n=1; n<=nMax; n++) {
2254
2255     TGraphErrors* gALL = VnDVsQ(n, k, "RIDGE", "cA", opt);
2256     double xerr = opt.Contains("sys") ? 0. : 0.0; // 0.2
2257     g->SetPoint(ip, n, 100*gALL->GetY()[q] );
2258     g->SetPointError(ip, xerr,  100*gALL->GetEY()[q] );
2259     if (0)
2260       cout<<gALL->GetName()<< " " << g->GetY()[ip] << endl;
2261
2262     ip++;
2263   }
2264
2265   utils.set_tgraph_props(g, kBlack, kBlack, kFullCircle, 1.2);
2266   g->SetLineWidth(2);
2267   gList->Add(g);
2268   return g;
2269 }
2270
2271 TGraphErrors* GlobalvnVsN(int i, int k, int nMax, TString opt)
2272 {
2273   // Create a TGraph from vn{GF} points
2274   // Stat errors by default. Call with opt "sys" to get sys errors.
2275
2276   TString gname(Form("GlobalvnVsN_pt%.2gto%.2g_cent%dto%d",
2277                      ptlow[i],pthigh[i], centlow[k],centhigh[k]));
2278   if (!opt.IsNull())  {
2279     gname.Append("_");
2280     gname.Append(opt);
2281   }
2282   TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data());
2283   if (g) {
2284     return g;
2285   }
2286   else {
2287     g = new TGraphErrors();
2288     g->SetName(gname.Data());
2289   }
2290   g->SetName(gname.Data());
2291
2292   int ip=0;
2293   for (int n=1; n<=nMax; n++) {
2294
2295     // point, stat err, sys err.
2296     double vnp = vnGF(n, k, i, 0, 999, "RIDGE", "cA", "");
2297     double vne = vnGF(n, k, i, 0, 999, "RIDGE", "cA", "err");
2298     double vns = vnGF(n, k, i, 0, 999, "RIDGE", "cA", "sys");
2299     double xerr = opt.Contains("sys") ? 0. : 0.0; // 0.2
2300     double yerr = opt.Contains("sys") ? vns : vne;
2301
2302     g->SetPoint(ip, n, vnp );
2303     g->SetPointError(ip, xerr, yerr);
2304     ip++;
2305   }
2306
2307   utils.set_tgraph_props(g, kBlack, kBlack, kFullCircle, 1.2);
2308   g->SetLineWidth(2);
2309
2310   return g;
2311 }
2312
2313 TGraphErrors* VnDeltaNFVsN(int i, int j, int k, int nMax, TString opt)
2314 {
2315
2316   TGraphErrors* VnD = VnDeltaVsN(i,j,k,nMax,opt);
2317   TGraphErrors* vnt = GlobalvnVsN(i,k,nMax,opt);
2318   TGraphErrors* vna = GlobalvnVsN(j,k,nMax,opt);
2319   TString gname(Form("VnDeltaNFVsN_ptt%.2gto%.2g_pta%.2gto%.2g_cent%dto%d",
2320                      ptlow[i],pthigh[i],ptlow[j],pthigh[j],
2321                      centlow[k],centhigh[k]));
2322   if (!opt.IsNull())  {
2323     gname.Append("_");
2324     gname.Append(opt);
2325   }
2326   TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data());
2327   if (g) {
2328     return g;
2329   }
2330   else {
2331     g = new TGraphErrors();
2332     g->SetName(gname.Data());
2333   }
2334   g->SetName(gname.Data());
2335
2336   int ip=0;
2337   for (int n=1; n<=nMax; n++) {
2338
2339     // Uncertainty vne is stat or sys depending on opt.
2340     double vnp = VnD->GetY()[ip] - 100*vnt->GetY()[ip]*vna->GetY()[ip];
2341     double vne = VnD->GetEY()[ip] - 100*vnt->GetEY()[ip]*vna->GetEY()[ip];
2342     g->SetPoint(ip, n, vnp );
2343     g->SetPointError(ip, 0.0, vne);
2344     ip++;
2345   }
2346
2347   utils.set_tgraph_props(g, kBlack, kBlack, kOpenSquare, 1.2);
2348   g->SetLineWidth(2);
2349   return g;
2350 }
2351
2352
2353 TF1* GlobalFitVsQ(int n, int k, int ipt1, int ipt2, const char* region,
2354                   const char* corrtype, TString opt)
2355 {
2356
2357   // vn{GF}(trig) x vn{GF}(assc) vs. q, where the points between (and
2358   // including) ipt1 and ipt2 were used in the global fit.
2359   // Note: gfname must be identical format to fnName in GlobalFitVsQ().
2360   TF1* fitfn = 0;
2361   TString errType("meas");
2362   if (opt.Contains("hi_sys"))
2363     errType = "hi";
2364   if (opt.Contains("lo_sys"))
2365     errType = "lo";
2366   TString gfname =
2367     Form("globalFitFn_cent%dto%d_n%d_ptbins%dto%d_%s_corrtype_%s_%s",
2368          centlow[k],centhigh[k], n, ipt1, ipt2,
2369          region, corrtype, errType.Data());
2370   fitfn = (TF1*)gFitList->FindObject(gfname.Data());
2371   if (fitfn) {
2372     if (0)
2373       Info("GlobalFitVsQ()", "Returning cached TF1 %s", gfname.Data());
2374     return fitfn;
2375   }
2376
2377   // If nonexistent, create (& store) a new one
2378   fitfn = DoGlobalFit(n,k,ipt1,ipt2,region,corrtype,opt);
2379   //  fitfn = (TF1*)gFitList->FindObject(gfname.Data());
2380
2381   if (!fitfn)
2382     Warning("GlobalFitVsQ()",
2383             "no fitfn %s,\neven after calling DoGlobalFit()",
2384             gfname.Data());
2385
2386   return fitfn;
2387 }
2388
2389 TF1* DoGlobalFit(int n, int k, int ipt1, int ipt2,
2390                  const char* region, const char* corrtype, TString opt)
2391 {
2392   // fnName must be identical format to gfname in GlobalFitVsQ().
2393   TF1 *globalFitFn = 0;
2394   TGraphErrors* g = 0;
2395   TString errType("meas");
2396   if (opt.Contains("hi_sys"))
2397     errType = "hi";
2398   if (opt.Contains("lo_sys"))
2399     errType = "lo";
2400   TString fnName =
2401     Form("globalFitFn_cent%dto%d_n%d_ptbins%dto%d_%s_corrtype_%s_%s",
2402          centlow[k],centhigh[k], n, ipt1, ipt2,
2403          region, corrtype, errType.Data());
2404
2405   g = VnDVsQ(n, k, region, corrtype, opt);
2406
2407   if (!g) {
2408     Error("DoGlobalFit()", "No graph found");
2409     return 0;
2410   }
2411
2412   // Assign global variables that will get checked in the vntvna
2413   // function during the fit.
2414   minFitPtBin = ipt1;
2415   maxFitPtBin = ipt2;
2416
2417   // Set up the global fit function...
2418   double qmax = g->GetN() - 0.4999;
2419   int nPars = maxpt + 3;
2420   globalFitFn = new TF1(fnName.Data(), vntvna, 0.0, qmax, nPars);
2421
2422   /* Not needed now, I just flip the v1 points later...AMA 11/15/2011 */
2423   // Each parameter has an equally good minimum at + or - v_nbest. The
2424   // starting value is chosen to steer the fit to the positive "root".
2425   // v_1 is special because it is the only coefficient that crosses
2426   // zero at low pt, so it gets a negative starting point. Pars are
2427   // vn(ipt), then last one is n.
2428   //  double startVal = (n==1)? -0.02 : 0.01;
2429
2430   double startVal = 0.01;
2431   for (int ip=0; ip<maxpt; ip++) {
2432     globalFitFn->SetParameter(ip, startVal);
2433     globalFitFn->SetParName(ip, Form("v_%d(%.3g)",n,ptmean[ip]));
2434   }
2435
2436   // Pass n into vntvna(). Access as par[maxpt].
2437   globalFitFn->FixParameter(maxpt, double(n));
2438
2439   // Pass k into vntvna(). Access as par[maxpt+1].
2440   globalFitFn->FixParameter(maxpt+1, double(k));
2441
2442   // Mom. conservation term prefactor. Access as par[maxpt+2].
2443   if (n==1) {
2444     globalFitFn->FixParameter(maxpt+2, 0.);
2445     //    globalFitFn->SetParLimits(maxpt+2, 0., 1e-3);
2446   }
2447   else
2448     globalFitFn->FixParameter(maxpt+2, 0.);
2449
2450   globalFitFn->SetNpx(1000);
2451
2452   TFitResultPtr result = g->Fit(globalFitFn,"QRNS"); // fast, takes ~15 ms
2453
2454   if (0) {
2455     cout<<globalFitFn->GetName();
2456     result->Print();
2457   }
2458   // TMinuit status:
2459   // migradResult + 10*minosResult + 100*hesseResult + 1000*improveResult
2460   // 0:ok, <0:user error, >0: see above.
2461   int fitStatus = result;
2462   if (fitStatus)
2463     Warning("DoGlobalFit()", "TMinuit status %d", fitStatus);
2464
2465   gFitList->Add(globalFitFn);
2466
2467   return globalFitFn;
2468 }
2469
2470 double vnGF(int n, int k, int ptBin, int ipt1, int ipt2,
2471             const char* region, const char* corrtype, TString opt)
2472 {
2473
2474   double vn, err;
2475   TF1* gfn = GlobalFitVsQ(n,k,ipt1,ipt2,region,corrtype,opt);
2476   if (!gfn)
2477     Error("vnGF()", "!gfn");
2478   vn = gfn->GetParameter(ptBin); // ptBin is the q index
2479   err = gfn->GetParError(ptBin);
2480
2481   if (opt.Contains("sys")) {
2482     TF1* hi = GlobalFitVsQ(n,k,ipt1,ipt2,region,corrtype,"hi_sys");
2483     TF1* lo = GlobalFitVsQ(n,k,ipt1,ipt2,region,corrtype,"lo_sys");
2484     double vhi = hi->GetParameter(ptBin);
2485     double vlo = lo->GetParameter(ptBin);
2486     return 0.5*TMath::Abs(vhi-vlo);
2487   }
2488
2489   // v1 can be negative, otherwise ensure the positive solution is chosen
2490   // if (n>1 && vn < 0)
2491   //   vn *= -1;
2492
2493   // if (ptlow[ipt1] > 4.) // !!!!!!!!!!!!!!!!!!!!!!!!!
2494   //   vn = TMath::Abs(vn);
2495
2496   // if (n==1)
2497   //   vn += MomCorrection(i, j, k);
2498
2499   return opt.Contains("err") ? err : vn;
2500 }
2501
2502 TGraphErrors* vnGFvsPt(int n, int k, int ipt1, int ipt2,
2503                        const char* region, const char* corrtype, TString opt)
2504 {
2505
2506   // vn{GF} vs. pt, where the points between (and including) ipt1 and
2507   // ipt2 were used in the global fit.
2508   bool isHighPtFit = ipt1 > 0 && ipt2==999;
2509   bool flipSign = n==1;  // Flip v1?
2510   if (isHighPtFit)
2511     flipSign = 0;
2512   if (opt.Contains("keepsign"))
2513     flipSign = 0;
2514   TString name(Form("v%dGFvsPt_fitpt%dto%d_cent%dto%d",
2515                   n, ipt1, ipt2, centlow[k], centhigh[k]));
2516   if (opt.Contains("sys"))
2517     name.Append("_sys");
2518   if (opt.Contains("ptcons"))
2519     name.Append("_ptcons");
2520
2521   TGraphErrors* g = (TGraphErrors*)gList->FindObject(name.Data());
2522   if (g)
2523     return g;
2524
2525   g = new TGraphErrors();
2526   g->SetName(name);
2527
2528   for (int i=0; i<maxpt; i++) {
2529     TString opt1 = "";
2530     if (opt.Contains("ptcons"))
2531       opt1.Append("_ptcons");
2532     double vn = vnGF(n,k,i,ipt1,ipt2,region,corrtype,opt1);
2533
2534     if (flipSign)
2535       vn = -vn;
2536
2537     //    if (isHighPtFit && i<ipt1)
2538     if (i<ipt1)
2539       vn = 0.;
2540
2541     double mpt = MeanPt(i, 0, k, "t");
2542     if (i<ipt1)
2543       mpt = -10.;
2544     g->SetPoint(i, mpt, vn);
2545
2546     if (opt.Contains("sys")) {
2547       opt1.Append("_sys");
2548       double ex  = (gi->GetEX()[i] < 0.126) ? gi->GetEX()[i] : 0.250;
2549       double vns = vnGF(n,k,i,ipt1,ipt2,region,corrtype, opt1);
2550       g->SetPointError(i, ex, vns);
2551       //      g->SetPointError(i, gi->GetEX()[i], vns); // too wide
2552     }
2553     else {
2554       opt1.Append("_err");
2555       double vne = vnGF(n,k,i,ipt1,ipt2,region,corrtype,opt1);
2556       g->SetPointError(i, 0.0, vne);
2557     }
2558   }
2559
2560   int col = CentColor(centlow[k], centhigh[k]);
2561   int col_pale = CentColorPale(centlow[k], centhigh[k]);
2562   g->SetLineWidth(2);
2563   g->SetFillColor(col_pale); // for sys
2564   utils.set_tgraph_props(g, col, col, kFullCircle, 1.2);
2565
2566   if (0) cout<<opt.Data()<<endl;
2567   gList->Add(g);
2568   return g;
2569 }
2570
2571 TGraph* Luzumv1(int cent, const char* pid)
2572 {
2573   // v1 for pions vs pt from viscous hydro calc by Matt Luzum
2574   // valid cent args 0-4 - 10% bins.
2575   // pid = "pi", "K", "p".
2576   const char* dat = Form("v1LHCGlauber%ssv08.dat", pid);
2577   TString fmt("%lg "); // pt
2578   for (int ky=0; ky<5; ky++)
2579     fmt.Append( ky==cent? "%lg " : "%*lg "); // centrality - skip if *lg
2580   return new TGraph(dat, fmt.Data());
2581 }
2582
2583 TGraphErrors* GlobalFitVsPtAssoc(int i, int k, int n, TString opt)
2584 {
2585   // Create a TGraph of vnt{gf} x vna{gf} vs. pt_assoc from gfit
2586   TF1* gfit = GlobalFitVsQ(n, k, 0, 999, "RIDGE", cfDef, opt);
2587   TGraphErrors* g = new TGraphErrors();
2588   TString gname(Form("GlobalFitVsPtAssoc_ptt%d_cent%d_n%d",i,k,n));
2589   if (!opt.IsNull())
2590     gname.Append(opt);
2591   g->SetName(gname);
2592
2593   for (int j=0; j<=i; j++) {
2594     int q = GlobalIndex(i, j);
2595     g->SetPoint(j, gj->GetX()[j], gfit->Eval(q + 0.5) );
2596     g->SetPointError(j, 0, 0 );
2597     utils.set_tgraph_props(g, colorsPale[i], colorsPale[i], kFullCircle, 1.0);
2598     g->SetLineWidth(4);
2599   }
2600   return g;
2601 }
2602
2603 double SineError(int i, int j, int k)
2604 {
2605   double rms = 0;
2606   double resid_sin=0, rs_err=0;
2607   TH1* hst  = Hist("RIDGE", "cA", i, j, k, "");
2608   double arr[100] = {0};
2609   //  arr[0] = 0;
2610   for (int n=1; n<=gNMax; n++) {
2611     VnDelta(n, *hst,  resid_sin, rs_err, "sin");
2612     arr[n-1] = resid_sin;
2613   }
2614
2615   rms = TMath::RMS(gNMax, arr) / 2;
2616   if (0)
2617     cout << " SineError: " << rms << endl;
2618
2619   return rms;
2620 }
2621
2622 TCanvas* DrawVnDeltaVsN(int nptt, int pttbins[], int npta,
2623                         int ptabins[], int ncb, int centbins[], TString opt)
2624 {
2625
2626   double t1=ptlow[pttbins[0]], t2=pthigh[pttbins[nptt-1]];
2627   double a1=ptlow[ptabins[0]], a2=pthigh[ptabins[npta-1]];
2628   TString cname = Form("VnDeltaVsN_ptt%.2gto%.2g_pta%.2gto%.2g",t1,t2,a1,a2);
2629   if (!opt.IsNull()) {
2630     cname.Append("_");
2631     cname.Append(opt);
2632   }
2633   TCanvas* c = new TCanvas(cname.Data(), cname.Data(), 650,700);
2634   double legbottom = 0.9 - 0.07*ncb;
2635   TLegend* leg = new TLegend(0.62, legbottom, 0.98, 0.98, "Centrality");
2636   leg->SetFillColor(kNone);
2637   leg->SetBorderSize(0);
2638   leg->SetName("leg");
2639
2640   TObjArray* graphs = new TObjArray();
2641   double xmin = 0.4;
2642   double xmax = opt.Contains("ext") ? gNMax + 0.5 : 8.5;
2643   int nXticks =  10, nYticks = 210;
2644   TLine l;
2645   TH1F* hf = gPad->DrawFrame(xmin, 0.0001, xmax, 0.26);
2646   //  bool scale1000 = 0;
2647
2648   hf->SetTitle(Form(";n;V_{n#Delta} [10^{-2}]"));
2649   if (opt.Contains("sine"))
2650     hf->SetTitle(Form(";n;#LTsin(n#Delta#phi)#GT [10^{-2}]"));
2651
2652   // if (opt.Contains("ext") )
2653   //   if (!opt.Contains("sine")) {
2654   //     scale1000 = 0;
2655   //     hf->SetTitle(Form(";n;V_{n#Delta}"));
2656   // }
2657
2658   l.SetLineStyle(kDashed);
2659
2660   utils.make_nice_axes(c, hf, 1.3, nXticks, nYticks, 1.4,-1.4);
2661   hf->SetLabelOffset(0.0, "X");
2662   hf->SetTitleOffset(-1.96, "Y");
2663
2664   hf->GetXaxis()->SetTicks("+-");
2665   hf->GetYaxis()->SetTicks("+-");
2666   l.DrawLine(hf->GetXaxis()->GetXmin(), 0, hf->GetXaxis()->GetXmax(),0);
2667   gPad->SetLeftMargin(0.22);
2668   gPad->SetTopMargin(0.01);
2669   gPad->SetRightMargin(0.01);
2670
2671   for (int i=0; i<nptt; i++) {
2672     for (int j=0; j<npta; j++) {
2673       for(int k=0; k<ncb; k++) {
2674
2675         // Stat errors on g, sys errors on gs.
2676         TString opt1 = opt.Contains("sine") ? "sine" : "";
2677         TGraphErrors* g =
2678           VnDeltaVsN(pttbins[i], ptabins[j], centbins[k], gNMax, opt1);
2679         TGraphErrors* gs =
2680           VnDeltaVsN(pttbins[i], ptabins[j], centbins[k], gNMax, "sys");
2681         TGraphErrors* gc = (TGraphErrors*)g->Clone();
2682         int dark = CentColor(centlow[centbins[k]], centhigh[centbins[k]]);
2683         int pale = CentColorPale(centlow[centbins[k]], centhigh[centbins[k]]);
2684         utils.set_tgraph_props(g, dark, dark, kFullCircle, 1.5);
2685         utils.set_tgraph_props(gs, dark, dark, kDot, 2.);
2686         utils.set_tgraph_props(gc, kBlack, kBlack, kOpenCircle, 1.5);
2687         g->SetFillColor(pale);
2688         gs->SetFillColor(dark);
2689
2690         TGraphErrors* gnf = 0;
2691         TGraphErrors* gnf_sys = 0;
2692         TGraphErrors* gnf_solid = 0;
2693         if (opt.Contains("vnf")) {  // Non-factorizing part
2694           gnf =
2695             VnDeltaNFVsN(pttbins[i], ptabins[j], centbins[k], gNMax, "");
2696           gnf_sys =
2697             VnDeltaNFVsN(pttbins[i], ptabins[j], centbins[k], gNMax, "sys");
2698
2699           // Just for appearance...
2700           gnf_solid = (TGraphErrors*) gnf->Clone(Form("%s_solid", gnf->GetName()));
2701           utils.set_tgraph_props(gnf_solid, pale, pale, kFullSquare, gnf->GetMarkerSize());
2702         }
2703         // if (scale1000) {
2704         //   for (int ip=0; ip<g->GetN(); ip++) {
2705         //     g->SetPoint(ip,   g->GetX()[ip], 1000* g->GetY()[ip]);
2706         //     gs->SetPoint(ip, gs->GetX()[ip], 1000*gs->GetY()[ip]);
2707         //     gc->SetPoint(ip, gc->GetX()[ip], 1000*gc->GetY()[ip]);
2708
2709         //     g->SetPointError(ip,   g->GetEX()[ip], 1000* g->GetEY()[ip]);
2710         //     gs->SetPointError(ip, gs->GetEX()[ip], 1000*gs->GetEY()[ip]);
2711         //     gc->SetPointError(ip, gc->GetEX()[ip], 1000*gc->GetEY()[ip]);
2712         //   }
2713         // }
2714
2715         if (opt.Contains("sine")) {
2716           g->Fit("pol0", "Q");
2717           TF1* gfit = g->GetFunction("pol0");
2718           gfit->SetLineColor(dark);
2719         }
2720
2721         else if (!opt.Contains("nobars"))
2722           g->Draw("Bsame");  // bars
2723
2724         if (opt.Contains("vnf")) {  // Non-factorizing part
2725           utils.draw_errorbox(gnf_sys, pale, 0.3);
2726           gnf_solid->Draw("epsame");
2727           gnf->Draw("epsame");
2728         }
2729
2730         g->Draw("psame");  // colored markers
2731         gc->Draw("epsame");  // open markers w/black stat error bars
2732
2733         if (!opt.Contains("sine"))
2734           utils.draw_errorbox(gs, dark, 0.3);
2735
2736         graphs->Add(g);
2737         leg->AddEntry(g,Form("%d-%d%%",centlow[centbins[k]],centhigh[centbins[k]]),"l,p");
2738
2739         if (opt.Contains("vnf")) {  // Non-factorizing part
2740           leg->AddEntry(gnf_solid,"NF V_{n#Delta}","l,p");
2741           //      leg->AddEntry(gnf_solid,"[v_{n}(p_{T}^{t}) #times v_{n}(p_{T}^{a})]_{GF}","l,p");
2742           leg->SetHeader("");
2743         }
2744
2745       }
2746     }
2747   }
2748
2749   // Adjust y limits
2750   double newYmin = 0, newYmax = 0, newSpace = 0;
2751   for (int m=0; m<graphs->GetEntries(); m++) {
2752     TGraphErrors* tg = (TGraphErrors*)graphs->At(m);
2753     double yhi = TMath::MaxElement(tg->GetN()-1, tg->GetY());
2754     double ylo = TMath::MinElement(tg->GetN()-1, tg->GetY());
2755     if (yhi > newYmax) newYmax = yhi;
2756     if (ylo < newYmin) newYmin = ylo;
2757   }
2758
2759   newSpace = 0.1*(newYmax - newYmin);
2760
2761   if (opt.Contains("ext"))
2762     hf->GetYaxis()->SetRangeUser(newYmin-newSpace, newYmax+4*newSpace);
2763   else
2764     hf->GetYaxis()->SetRangeUser(newYmin-newSpace, newYmax+newSpace);
2765
2766   // if (opt.Contains("sine")) {
2767   //   hf->GetYaxis()->SetRangeUser(0.0001, newYmax+4*newSpace);
2768   //   gPad->SetLogy();
2769   // }
2770
2771   if (!opt.Contains("sine"))
2772     leg->Draw();
2773
2774   // Draw text
2775   ltx.SetTextSize(0.04); // 0.05
2776
2777   if ( opt.Contains("sine") )
2778     legbottom = 0.99;
2779   else if (opt.Contains("ext") && !opt.Contains("vnf") )
2780     legbottom -= 0.4;
2781
2782   double xstart = 0.6;
2783   ltx.DrawLatex(xstart, legbottom-1*0.07, Form("%.2g < p_{T}^{t} < %.2g GeV/c", t1, t2));
2784   ltx.DrawLatex(xstart, legbottom-2*0.07, Form("%.2g < p_{T}^{a} < %.2g GeV/c", a1, a2));
2785
2786   gPad->Update();
2787   gPad->Modified();
2788   return c;
2789 }
2790
2791 TCanvas* DrawGlobalvnVsN(int npt, int ptbins[], int ncb, int centbins[], TString opt)
2792 {
2793   double t1=ptlow[ptbins[0]], t2=pthigh[ptbins[npt-1]];
2794   TString cname = Form("GlobalvnVsN_pt%.2gto%.2g",t1,t2);
2795   if (!opt.IsNull()) {
2796     cname.Append("_");
2797     cname.Append(opt);
2798   }
2799   TCanvas* c = new TCanvas(cname.Data(), cname.Data(), 600,700);
2800   double legbottom = 0.9 - 0.07*ncb;
2801   TLegend* leg = new TLegend(0.62, legbottom, 0.98, 0.98, "Centrality");
2802   leg->SetFillColor(kNone);
2803   leg->SetBorderSize(0);
2804   leg->SetName("leg");
2805
2806   TObjArray* graphs = new TObjArray();
2807   double xmin = opt.Contains("zoom") ? 2.6 : 0.4;
2808   double xmax = 8.5;
2809   int nXticks =  10, nYticks = 210;
2810   TLine l;
2811   TH1F* hf = gPad->DrawFrame(xmin, -0.06, xmax, 0.26);
2812   hf->SetTitle(Form(";n;v_{n}{GF}"));
2813   l.SetLineStyle(kDashed);
2814
2815   utils.make_nice_axes(c, hf, 1.3, nXticks, nYticks, 1.4,-1.4);
2816   hf->SetLabelOffset(0.0, "X");
2817   if (opt.Contains("zoom"))
2818     hf->SetTitleOffset(-2.3, "Y");
2819   else
2820     hf->SetTitleOffset(-1.9, "Y");
2821
2822   hf->GetXaxis()->SetTicks("+-");
2823   hf->GetYaxis()->SetTicks("+-");
2824   l.DrawLine(hf->GetXaxis()->GetXmin(), 0, hf->GetXaxis()->GetXmax(),0);
2825   if (opt.Contains("zoom"))
2826     gPad->SetLeftMargin(0.25);
2827   else
2828     gPad->SetLeftMargin(0.2);
2829   gPad->SetTopMargin(0.01);
2830   gPad->SetRightMargin(0.01);
2831
2832   for (int i=0; i<npt; i++) {
2833     for(int k=0; k<ncb; k++) {
2834
2835       // Stat errors on g, sys errors on gs.
2836       TGraphErrors* g  = GlobalvnVsN(ptbins[i], centbins[k], gNMax, "");
2837       TGraphErrors* gs = GlobalvnVsN(ptbins[i], centbins[k], gNMax, "sys");
2838       TGraphErrors* gc = (TGraphErrors*)g->Clone();
2839       int dark = CentColor(centlow[centbins[k]], centhigh[centbins[k]]);
2840       int pale = CentColorPale(centlow[centbins[k]], centhigh[centbins[k]]);
2841       utils.set_tgraph_props(g, dark, dark, kFullCircle, 1.5);
2842       //      utils.set_tgraph_props(gs, dark, dark, kDot, 2.);
2843       utils.set_tgraph_props(gc, kBlack, kBlack, kOpenCircle, 1.5);
2844       g->SetFillColor(pale);
2845       //      gs->SetFillColor(dark);
2846       g->Draw("Bsame");  // colored bars
2847       utils.draw_errorbox(gs, dark, 0.3);
2848       g->Draw("psame");  // colored markers
2849       //      gs->Draw("e[]psame"); // sys error filled rectangles
2850       gc->Draw("epsame");  // open markers w/black stat error bars
2851       graphs->Add(g);
2852       leg->AddEntry(g,Form("%d-%d%%",centlow[centbins[k]],centhigh[centbins[k]]),"l,p");
2853     }
2854   }
2855
2856   // Adjust y limits
2857   double newYmin = 0, newYmax = 0, newSpace = 0;
2858   for (int m=0; m<graphs->GetEntries(); m++) {
2859     TGraphErrors* tg = (TGraphErrors*)graphs->At(m);
2860     double yhi = TMath::MaxElement(tg->GetN()-1, tg->GetY());
2861     double ylo = TMath::MinElement(tg->GetN()-1, tg->GetY());
2862     if (yhi > newYmax) newYmax = yhi;
2863     if (ylo < newYmin) newYmin = ylo;
2864   }
2865
2866   newSpace = 0.1*(newYmax - newYmin);
2867   hf->GetYaxis()->SetRangeUser(newYmin-newSpace, newYmax+newSpace);
2868   if (opt.Contains("zoom"))
2869     hf->GetYaxis()->SetRangeUser(-0.0008, 0.0085);
2870
2871   leg->Draw();
2872   c->Update();
2873
2874   // Draw text
2875   // if (opt.Contains("zoom"))
2876   //   ltx.SetTextSize(0.03); // 0.04
2877   // else
2878
2879   // ltx.SetTextSize(0.04); // 0.05
2880
2881   double xstart = (opt.Contains("zoom")) ? 0.7 : 0.6;
2882   ltx.DrawLatex(xstart, legbottom-1*0.07, Form("%.2g < p_{T} < %.2g GeV/c", t1, t2));
2883   // ltx.DrawLatex(xstart, legbottom-2*0.07, Form("%.2g < p_{T}^{assoc} < %.2g GeV/c", a1, a2));
2884
2885   return c;
2886 }
2887
2888
2889 TCanvas* DrawVnFromGlobalFit(int n, int ipt1, int ipt2, int ncb,
2890                              int centbins[], TString opt)
2891 {
2892   // Draw vn for one n on one canvas - many centralities
2893   initialize();
2894   double maxPtFit = ipt2==999 ? pthigh[maxpt-1] : pthigh[ipt2];
2895   const char* cname = Form("v%d_global%s_ptfit%.2gto%.2g",
2896                            n, opt.Data(), ptlow[ipt1], maxPtFit);
2897   const char* title = Form("v%d_global%s_ptfit%.2gto%.2g",
2898                            n, opt.Data(), ptlow[ipt1], maxPtFit);
2899   double lvleft = opt.Contains("multi")? 0.74 : 0.82;
2900   double lvbot = opt.Contains("multi")? 0.25 : 0.56;
2901   int cvHeight = opt.Contains("multi")? 500 : 500;
2902   TCanvas* cv = new TCanvas(cname, title, 700, cvHeight);
2903   TLegend* lv = new TLegend(lvleft, lvbot, 0.99, 0.97);
2904   lv->SetFillColor(kNone); lv->SetName("lv");
2905   lv->SetMargin(0.3);
2906   lv->SetBorderSize(1);
2907   lv->SetHeader("#splitline{Fit range}{(GeV/c):}");
2908   utils.padsetup(cv, 1, 1, "", 0.12, 0.01, 0.03, 0.2); // ######
2909   //  gPad->SetRightMargin(0.15);
2910
2911
2912   TPad* inset = 0;
2913
2914   TObjArray* graphs = new TObjArray();
2915   double xmin = 0., xmax = 14.8; // 8.3;
2916   double ymin = -0.06, ymax = 0.26;
2917   if (opt.Contains("multi")) {
2918     ymin = -1.19;
2919     ymax = 0.52;
2920   }
2921   TH1F* hf = gPad->DrawFrame(xmin, ymin, xmax, ymax);
2922   hf->SetTitle(Form(";p_{T}^{t} (GeV/c);v_{%d}{GF}",n));
2923   int nXticks =  210;
2924   int nYticks = 210;
2925   TLine l;
2926   l.SetLineStyle(kDashed);
2927   if (n==1) {
2928     cv->cd();
2929     inset = new TPad("inset", "inset", 0.18, 0.18, 0.45, 0.59, kNone, 1, 0);
2930     inset->Draw();
2931     inset->cd();
2932     TH1F* hfi = gPad->DrawFrame(-0.05, -0.04, 4.2, 0.16);
2933     TAxis* hx = hfi->GetXaxis();
2934     TAxis* hy = hfi->GetYaxis();
2935     hx->SetLabelSize(0.1);
2936     hy->SetLabelSize(0.1);
2937     hx->SetNdivisions(104);
2938     hy->SetNdivisions(104);
2939     gPad->SetBottomMargin(0.2);
2940     gPad->SetLeftMargin(0.2);
2941     gPad->SetRightMargin(0.01);
2942     gPad->SetTopMargin(0.01);
2943     //    utils.make_nice_axes(inset, hfi, 3, nXticks, nYticks, 1.4,-1.4);
2944     l.DrawLine(-0.05, 0, 4.2, 0);
2945     cv->cd();
2946   }
2947
2948   utils.make_nice_axes(cv, hf, 1.3, nXticks, nYticks, 1.4,-1.4);
2949   hf->SetLabelOffset(0.0, "X");
2950   hf->SetTitleOffset(1.5, "Y");
2951   hf->GetXaxis()->SetTicks("+");
2952   //  hf->GetYaxis()->SetTicks("+-");
2953
2954   TGraph* v1h[5][3];
2955   if (opt.Contains("luzum")) {
2956     int lSty[] = {7, kSolid, kSolid};
2957     int lCol[] = {kBlack, kCyan, kBlue, kGreen+1, kRed};
2958     TString pidString[] = {"pi", "K", "p"};
2959     for (int m=0; m<5; m++) {
2960       for (int pid=0; pid<3; pid++) {
2961         v1h[m][pid] = Luzumv1(m, pidString[pid].Data());
2962         TGraph* g1 = v1h[m][pid];
2963         g1->SetLineStyle(lSty[pid]);
2964         g1->SetLineColor(lCol[m]);
2965         g1->SetLineWidth(2);
2966       }
2967     }
2968     hf->GetXaxis()->SetRangeUser(0, 4.8);
2969     l.DrawLine(0, 0, 4.8, 0);
2970     v1h[0][0]->Draw("plsame");
2971     v1h[0][2]->Draw("plsame");
2972     v1h[4][0]->Draw("plsame");
2973     v1h[4][2]->Draw("plsame");
2974   }
2975   else
2976     l.DrawLine(hf->GetXaxis()->GetXmin(), 0, hf->GetXaxis()->GetXmax(),0);
2977
2978   for (int cb=0; cb<ncb; cb++) {
2979     int k = centbins[cb];
2980
2981     if (!opt.Contains("multi")) {
2982       // last 2 args are included pt bins
2983       TGraphErrors* gtmp = vnGFvsPt(n, k, ipt1, ipt2);
2984
2985       // clone to modify without ripples
2986       TGraphErrors* gc = (TGraphErrors*) gtmp->Clone();
2987       if (!gc) {
2988         Error("DrawVnFromGlobalFit()", "Problem making graph" );
2989         continue;
2990       }
2991       TGraphErrors* gtmp_sys = vnGFvsPt(n, k, ipt1, ipt2, "RIDGE", cfDef, "sys");
2992       TGraphErrors* gs = (TGraphErrors*) gtmp_sys->Clone();
2993
2994       if (opt.Contains("luzum")) {
2995         if (centlow[k]==20)
2996           utils.set_tgraph_props(gc, kBlue, kBlue, kFullCircle, 1.2);
2997         if (centlow[k]==40)
2998           utils.set_tgraph_props(gc, kRed, kRed, kFullCircle, 1.2);
2999       }
3000
3001       utils.draw_errorbox(gs, gs->GetFillColor());
3002       gc->Draw("epsame");
3003       graphs->Add(gc);
3004
3005       TGraphErrors* open = (TGraphErrors*) gc->Clone(Form("%s_open", gc->GetName()));
3006       utils.set_tgraph_props(open, kBlack, kBlack, kOpenCircle, 1.2);
3007       open->Draw("epsame");
3008
3009       lv->AddEntry(gc,Form("%d-%d%%",centlow[k],centhigh[k]),"l,p");
3010     }
3011     // NEW
3012     if (opt.Contains("multi")) {
3013       int b1,b2,b3,b4;
3014       b1=PtBin(0.25,0.5);
3015       b2=PtBin(0.75,1.0);
3016       b3=PtBin(2.0, 2.5);
3017       b4=PtBin(3.0, 4.0);
3018       TGraphErrors* gg1 = vnGFvsPt(n, k, b1, b2, "RIDGE", cfDef, "");
3019       TGraphErrors* gg2 = vnGFvsPt(n, k, b3, b4, "RIDGE", cfDef, "keepsign");
3020       TGraphErrors* gs1 = vnGFvsPt(n, k, b1, b2, "RIDGE", cfDef, "sys");
3021       TGraphErrors* gs2 = vnGFvsPt(n, k, b3, b4, "RIDGE", cfDef, "sys_keepsign");
3022       utils.set_tgraph_props(gg1, gg1->GetLineColor(), gg1->GetMarkerColor(), kOpenCircle, 1.0);
3023       utils.set_tgraph_props(gg2, gg2->GetLineColor(), gg2->GetMarkerColor(), kFullSquare, 1.0);
3024
3025       utils.draw_errorbox(gs1, gs1->GetFillColor());
3026       utils.draw_errorbox(gs2, gs2->GetFillColor());
3027       gg1->Draw("epsame");
3028       gg2->Draw("epsame");
3029       graphs->Add(gs1);
3030       graphs->Add(gs2);
3031       const char* s1 = Form("#splitline{%.2g < p_{T}^{a} < %.2g  }{%d-%d%%}",
3032                             ptlow[b1], pthigh[b2], centlow[k],centhigh[k]);
3033       const char* s2 = Form("#splitline{%.2g < p_{T}^{a} < %.2g  }{%d-%d%%}",
3034                             ptlow[b3], pthigh[b4], centlow[k],centhigh[k]);
3035       lv->AddEntry(gg1, s1, "l,p");
3036       lv->AddEntry(gg2, s2, "l,p");
3037
3038
3039       if (n==1) {
3040         inset->cd();
3041         gg1->Draw("epsame");
3042         gg2->Draw("epsame");
3043         cv->cd();
3044       }
3045
3046     }
3047
3048   } // centrality bin loop
3049
3050   if (opt.Contains("luzum")) {
3051     lv->AddEntry(v1h[0][0], Form("v.h. (#pi)"),"l,p");
3052     lv->AddEntry(v1h[0][2], Form("v.h. (p)"), "l,p");
3053   }
3054   lv->Draw();
3055   ltx.SetTextSize(0.05);
3056
3057   if (opt.Contains("multi")) {
3058     gPad->SetTopMargin(0.02);
3059     if (ispp)
3060     ltx.DrawLatex(0.19, 0.92, Form("%s", "#splitline{pp}{2.76 TeV}"));
3061     else
3062     ltx.DrawLatex(0.19, 0.92, Form("%s", "#splitline{Pb-Pb}{2.76 TeV}"));
3063     ltx.DrawLatex(0.49, 0.94, Form("%.2g < |#Delta#eta| < 1.8", minRidgeDEta) );
3064
3065   }
3066   else {
3067     ltx.DrawLatex(0.17, 0.88, Form("%.2g < |#Delta#eta| < 1.8", minRidgeDEta) );
3068     ltx.DrawLatex(0.17, 0.94, Form("%.2g < fit p_{T}^{a} < %.2g GeV/c",
3069                                    ptlow[ipt1], maxPtFit));
3070     ltx.DrawLatex(0.19, 0.78, Form("%s", "#splitline{Pb-Pb}{2.76 TeV}"));
3071   }
3072   double newYmin = 0, newYmax = 0, newSpace = 0;
3073   for (int m=0; m<graphs->GetEntries(); m++) {
3074     TGraphErrors* tg = (TGraphErrors*)graphs->At(m);
3075     double yhi = TMath::MaxElement(tg->GetN()-1, tg->GetY());
3076     double ylo = TMath::MinElement(tg->GetN()-1, tg->GetY());
3077     if (yhi > newYmax) newYmax = yhi;
3078     if (ylo < newYmin) newYmin = ylo;
3079   }
3080
3081   if (opt.Contains("multi")) {
3082     if(n==1)
3083       newYmin *= 2.2;
3084     if(n==2)
3085       newYmin = 0.02;
3086   }
3087
3088   newSpace = 0.2*(newYmax - newYmin);
3089   //  if (opt.Contains("multi")) newSpace /= 2;
3090   if (n==5) newSpace *= 2;
3091   hf->GetYaxis()->SetRangeUser(newYmin-newSpace, newYmax+newSpace);
3092   cv->Modified();
3093   cv->Update();
3094   cv->Draw();
3095   return cv;
3096 }
3097
3098 TCanvas* Drawv1to5(int ncb, int centbins[], int ipt1, int ipt2, TString opt)
3099 {
3100   if (!opt.IsNull())
3101     cout<<opt.Data()<<endl;
3102
3103   initialize();
3104   int cent1 = centlow[0];
3105   int cent2 = centhigh[centbins[ncb-1]];
3106   double pmax = 15.;
3107   if (ipt2 < 999.)
3108     pmax = pthigh[ipt2];
3109   const char* cname = Form("vn_etamin%02d_cent%dto%d_fitptbin%dto%d",
3110                            (int)(10*minRidgeDEta), cent1, cent2, ipt1, ipt2);
3111   const char* title = Form("v1to5gf_etamin%02d_cent%dto%d_fitpt%.2gto%.2g",
3112                            (int)(10*minRidgeDEta), cent1, cent2, ptlow[ipt1], pmax);
3113   TCanvas* cv = new TCanvas(cname, title, 1400, 500);
3114   TLegend* lv = new TLegend(0.09, 0.72, 0.27, 0.88);
3115   lv->SetFillColor(kNone);
3116   lv->SetBorderSize(0);
3117   TLine zero;
3118   double lv2bottom = 0.62;
3119   TLegend* lv2 = new TLegend(0.82, lv2bottom, 0.93, 0.95, "Centrality");
3120   lv2->SetFillColor(kNone);
3121   lv2->SetBorderSize(0);
3122   utils.padsetup(cv, 5, 1, "", 0.12, 0.01, 0.03, 0.2);
3123
3124   for (int n=1; n<=5; n++) {
3125     Printf("Drawv1to5() n = %d / 5", n);
3126     cv->cd(n);
3127     double xmin = -0.43, xmax =  8.3;
3128     double ymin = -0.09, ymax =  0.26;
3129     if (ptlow[ipt1] > 4.) {
3130       xmin = 4.9;
3131       ymin = -0.1;
3132     }
3133     if (pthigh[ipt2] > 6.) {
3134       xmax = 11.8;
3135       ymax = 1;
3136     }
3137     if (opt.Contains("split")) {
3138       xmin = -0.43;
3139       xmax = 11.8;
3140       ymin = -0.2;
3141       ymax = 1;
3142     }
3143
3144     TH1F* hf = gPad->DrawFrame(xmin, ymin, xmax, ymax);
3145     int nXticks = 206;
3146     int nYticks = 208;
3147     if (n==1) {
3148       utils.make_nice_axes(cv, hf, 2., nXticks, nYticks, -0.02, 0.02);
3149       hf->SetLabelOffset(-0.01, "X");
3150     }
3151     else {
3152       utils.make_nice_axes(cv, hf, 3.4, nXticks, nYticks, -0.075, 0.0);
3153       hf->SetLabelOffset(-0.057, "X");
3154       hf->SetTickLength(0.02, "X");
3155       hf->SetTickLength(0.05, "Y");
3156     }
3157     hf->GetXaxis()->SetTicks("+-");
3158     hf->GetYaxis()->SetTicks("+-");
3159     zero.SetLineStyle(kDashed);
3160     zero.DrawLine(xmin, 0., xmax, 0.);
3161
3162     // Loop thru every cent bin. Skip non-requested ones.
3163     // for (int k=0; k<maxcent; k++) {
3164     //   bool centOK = 1;
3165     //   if (ncb != 999) {
3166     //  centOK = 0;
3167     //  for (int cb=0; cb<ncb; cb++)
3168     //    if(k==centbins[cb])
3169     //      centOK = 1;
3170     //   }
3171     //   if (!centOK)
3172     //  continue;
3173
3174     for (int cb=0; cb<ncb; cb++) {
3175       int k = centbins[cb];
3176       TGraphErrors* gc = vnGFvsPt(n, k, ipt1, ipt2);
3177       TGraphErrors* gs = vnGFvsPt(n, k, ipt1, ipt2, "RIDGE", cfDef, "sys");
3178
3179       if (!gc) {
3180         Error("Drawv1to5()", "Problem finding graph");
3181         continue;
3182       }
3183
3184       TGraphErrors* gc_hipt=0;
3185       TGraphErrors* gs_hipt=0;
3186       if (opt.Contains("highptfit") || opt.Contains("split")) {
3187         gc_hipt = vnGFvsPt(n, k, ipt2+1, 999, "highptfit");
3188         gs_hipt = vnGFvsPt(n, k, ipt2+1, 999, "RIDGE", cfDef, "sys");
3189         gc_hipt->SetMarkerStyle(kOpenCircle);
3190       }
3191
3192       int pale = gs->GetFillColor(); // it is CentColor(centlow[k], centhigh[k]).
3193       utils.draw_errorbox(gs, pale);
3194       gc->Draw("epsame");
3195
3196       if (opt.Contains("split")) {
3197         utils.draw_errorbox(gs_hipt, pale);
3198         gc_hipt->Draw("epsame");
3199       }
3200
3201       TGraphErrors* open = (TGraphErrors*) gc->Clone(Form("%s_open", gc->GetName()));
3202       utils.set_tgraph_props(open, kBlack, kBlack, kOpenCircle, 1.2);
3203       open->Draw("epsame");
3204
3205       if (n==1)
3206         lv2->AddEntry(gc,Form("%d-%d%%",centlow[k],centhigh[k]),"l,p");
3207
3208     } // cent bin k
3209   } // Fourier moment n
3210
3211   cv->cd();
3212   TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1);
3213   overlay->SetFillStyle(4000); // transparent
3214   overlay->Draw();
3215   overlay->cd();
3216
3217   ltx.SetTextSize(0.07);
3218   lv2->Draw();
3219
3220   for (int n=1; n<=5; n++) {
3221     cv->cd(n);
3222     ltx.SetTextSize(n==1? 0.1 : 0.15);
3223     ltx.DrawLatex((n==1)? 0.85 : 0.7, 0.86, Form("v_{%d}", n));
3224
3225     if (0 && n==1)
3226       AddPrelimStampToCurrentPad(0.45,0.7,0.75,0.95, "");
3227   }
3228   overlay->cd();
3229   ltx.SetTextSize(0.075);
3230   ltx.DrawLatex(0.5, 0.05, "p_{T} [GeV/c]");
3231   ltx.SetTextAngle(90);
3232   ltx.DrawLatex(0.05, 0.5, "v_{n }(p_{T})");
3233   ltx.SetTextAngle(0);
3234   ltx.SetTextSize(0.05);
3235   ltx.DrawLatex(0.14, 0.88, Form("%s", "#splitline{Pb-Pb}{2.76 TeV}"));
3236   ltx.SetTextSize(0.05);
3237   ltx.DrawLatex(0.13, 0.24, Form("%.2g < |#Delta#eta| < 1.8", minRidgeDEta) );
3238
3239   const char* frange = Form("%.2g < fit p_{T}^{a} < %.2g GeV",
3240                             ptlow[ipt1], pmax);
3241   const char* frange1 = Form("#splitline{solid:}{%.2g < fit p_{T}^{a} < %.2g GeV}",
3242                             ptlow[ipt1], pthigh[ipt2]);
3243   const char* frange2 = Form("#splitline{open:}{%.2g < fit p_{T}^{a} < %.2g GeV}",
3244                             ptlow[ipt2+1], pthigh[maxpt-1]);
3245
3246   if (opt.Contains("split")) {
3247     ltx.SetTextSize(0.04);
3248     ltx.DrawLatex(0.31, 0.8, frange1);
3249     ltx.DrawLatex(0.31, 0.7, frange2);
3250   }
3251   else
3252     ltx.DrawLatex(0.3, 0.27, frange);
3253
3254   return cv;
3255 }
3256
3257 TCanvas* Drawv2to5(int ncb, int centbins[], int ipt1, int ipt2, TString opt)
3258 {
3259   if (!opt.IsNull())
3260     cout<<opt.Data()<<endl;
3261
3262   initialize();
3263   int cent1 = centlow[0];
3264   int cent2 = centhigh[centbins[ncb-1]];
3265   double pmax = 15.;
3266   if (ipt2 < 999.)
3267     pmax = pthigh[ipt2];
3268   const char* cname = Form("v2to5_etamin%02d_cent%dto%d_fitptbin%dto%d",
3269                            (int)(10*minRidgeDEta), cent1, cent2, ipt1, ipt2);
3270   const char* title = Form("v2to5gf_etamin%02d_cent%dto%d_fitpt%.2gto%.2g",
3271                            (int)(10*minRidgeDEta), cent1, cent2, ptlow[ipt1], pmax);
3272   TCanvas* cv = new TCanvas(cname, title, 1400, 500);
3273   TLine zero;
3274   double lv2bottom = 0.58; // 0.5;
3275   TLegend* lv2 = new TLegend(0.77, lv2bottom, 0.9, 0.95, "Centrality");
3276   lv2->SetFillColor(kNone);
3277   lv2->SetBorderSize(0);
3278   utils.padsetup(cv, 4, 1, "", 0.08, 0.01, 0.03, 0.2);
3279
3280   for (int n=2; n<=5; n++) {
3281     Printf("Drawv2to5() n = %d", n);
3282     cv->cd(n-1);
3283     double xmin = -0.2;
3284     double xmax =  7.6;
3285     TH1F* hf = gPad->DrawFrame(xmin, -0.12, xmax, 0.26); //was y = -0.03, 0.26
3286     int nXticks = 206;
3287     int nYticks = 208;
3288     if (n==2) {
3289       utils.make_nice_axes(cv, hf, 2., nXticks, nYticks, -0.02, 0.02);
3290       hf->SetLabelOffset(0.0, "X"); //       hf->SetLabelOffset(-0.01, "X");
3291     }
3292     else {
3293       utils.make_nice_axes(cv, hf, 2.6, nXticks, nYticks, -0.075, 0.0);
3294       hf->SetLabelOffset(-0.02, "X");
3295       hf->SetTickLength(0.02, "X");
3296       hf->SetTickLength(0.05, "Y");
3297     }
3298     hf->GetXaxis()->SetTicks("+-");
3299     hf->GetYaxis()->SetTicks("+-");
3300     zero.SetLineStyle(kDashed);
3301     zero.DrawLine(xmin, 0., xmax, 0.);
3302
3303     for (int cb=0; cb<ncb; cb++) {
3304       int k = centbins[cb];
3305
3306       TGraphErrors* gc = vnGFvsPt(n, k, ipt1, ipt2);
3307       TGraphErrors* gs = vnGFvsPt(n, k, ipt1, ipt2, "RIDGE", cfDef, "sys");
3308
3309       if (!gc) {
3310         Error("Drawv2to5()", "Problem finding graph");
3311         continue;
3312       }
3313
3314       int pale = gs->GetFillColor(); // it is CentColor(centlow[k], centhigh[k]).
3315       utils.draw_errorbox(gs, pale);
3316       //      gs->Draw("e2psame");
3317       gc->Draw("epsame");
3318       TGraphErrors* open = (TGraphErrors*) gc->Clone(Form("%s_open", gc->GetName()));
3319       utils.set_tgraph_props(open, kBlack, kBlack, kOpenCircle, 1.2);
3320       open->Draw("epsame");
3321
3322       if (n==2)
3323         lv2->AddEntry(gc,Form("%d-%d%%",centlow[k],centhigh[k]),"l,p");
3324
3325     } // cent bin k
3326   } // Fourier moment n
3327
3328   cv->cd();
3329   TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1);
3330   overlay->SetFillStyle(4000); // transparent
3331   overlay->Draw();
3332   overlay->cd();
3333
3334   ltx.SetTextSize(0.07);
3335   lv2->Draw();
3336
3337   for (int n=2; n<=5; n++) {
3338     cv->cd(n-1);
3339     ltx.SetTextSize(n==2? 0.1 : 0.15);
3340     ltx.DrawLatex((n==2)? 0.85 : 0.7, 0.86, Form("v_{%d}", n));
3341   }
3342   overlay->cd();
3343   ltx.SetTextSize(0.075);
3344   ltx.DrawLatex(0.48, 0.04, "p_{T} [GeV/c]");
3345   ltx.SetTextAngle(90);
3346   ltx.DrawLatex(0.02, 0.5, "v_{n }{GF}");
3347   ltx.SetTextAngle(0);
3348   ltx.SetTextSize(0.05);
3349   ltx.DrawLatex(0.33, 0.88, Form("%s", "#splitline{Pb-Pb}{2.76 TeV}"));
3350   ltx.SetTextSize(0.05);
3351   ltx.DrawLatex(0.33, 0.75, Form("%.2g < |#Delta#eta| < 1.8", minRidgeDEta) );
3352
3353   const char* frange = Form("#splitline{Global fit range}{%.2g < p_{T} < %.2g GeV}",
3354                             ptlow[ipt1], pmax);
3355   ltx.DrawLatex(0.56, 0.75, frange);
3356
3357   return cv;
3358 }
3359
3360 TCanvas* DrawChi2vsPtTrig(int npta, int ptabins[], int ncb, int centbins[], TString opt)
3361 {
3362
3363   // TODO: options-- global or 1to5, chi2/N or prob
3364   double a1=ptlow[ptabins[0]], a2=pthigh[ptabins[npta-1]];
3365   int c1=centlow[centbins[0]], c2=centhigh[centbins[ncb-1]];
3366   TString cname = Form("chi2_etamin%02d_pta%.2gto%.2g_cent%dto%d",
3367                        (int)(10*minRidgeDEta), a1, a2, c1, c2 );
3368
3369   if (!opt.IsNull() ) {
3370     cname.Append("_");
3371     cname.Append(opt.Data());
3372   }
3373
3374   TCanvas* cv = new TCanvas(cname, cname, 1);
3375   TLegend* lv = new TLegend(0.09, 0.72, 0.27, 0.88); // w.r.t. cv
3376   lv->SetFillColor(kNone);
3377   lv->SetBorderSize(0);
3378
3379   double lv2bottom = 0.62;
3380   TLegend* lv2 = new TLegend(0.82, lv2bottom, 0.93, 0.95, "Centrality"); // w.r.t. cv
3381   lv2->SetFillColor(kNone);
3382   lv2->SetBorderSize(0);
3383   utils.padsetup(cv, 1, 1, "", 0.12, 0.01, 0.03, 0.2);
3384
3385   double xmin = -0.43;
3386   double xmax =  12.5;
3387   TH1F* hf = gPad->DrawFrame(xmin, -0.09, xmax, 26);
3388   // int nXticks = 210;
3389   // int nYticks = 210;
3390   hf->GetXaxis()->SetTicks("+-");
3391   hf->GetYaxis()->SetTicks("+-");
3392
3393 // IN PROGRESS---------------
3394     for (int j=0; j<npta; j++) {
3395       for(int k=0; k<ncb; k++) {
3396
3397         TGraphErrors *gc = CorrFnChi2VsTrigPt(ptabins[j], centbins[k], 1, 5, "ndf");
3398         if (!gc) {
3399           Error("DrawChi2()", "!gc");
3400           continue;
3401         }
3402         int col = CentColor(centlow[centbins[k]], centhigh[centbins[k]]);
3403         utils.set_tgraph_props(gc, col, col, kFullCircle, 1.6);
3404         gc->Draw("epsame");
3405         lv2->AddEntry(gc, centLabel(centbins[k]), "l,p");
3406       }
3407     }
3408
3409   // cv->cd();
3410   // TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1);
3411   // overlay->SetFillStyle(4000); // transparent
3412   // overlay->Draw();
3413   // overlay->cd();
3414
3415   // ltx.SetTextSize(0.07);
3416   // lv2->Draw();
3417   // overlay->cd();
3418   // ltx.SetTextSize(0.075);
3419   // ltx.DrawLatex(0.5, 0.05, "trigger p_{T} [GeV/c]");
3420   // ltx.SetTextAngle(90);
3421   // ltx.DrawLatex(0.05, 0.5, "#chi^{2} / NDF");
3422   // ltx.SetTextAngle(0);
3423
3424   return cv;
3425 }
3426
3427
3428 TCanvas* DrawAgreement(int icent, TString opt)
3429 {
3430   // Draw Chi2 on one canvas for each
3431   // centrality bin k.
3432
3433   initialize();
3434   int cent1 = centlow[icent];
3435   int cent2 = centhigh[icent];
3436
3437   if (icent==999) {
3438     cent1 = 0;
3439     cent2 = 90;
3440   }
3441   int k = icent;
3442   TString cname = Form("agmt_etamin%02d_cent%dto%d_%s",
3443                        (int)(10*minRidgeDEta), cent1, cent2, opt.Data());
3444   if (!opt.IsNull()) {
3445     cname.Append("_");
3446     cname.Append(opt.Data());
3447   }
3448
3449   TCanvas* cv = new TCanvas(cname.Data(), cname.Data(), 1400, 500);
3450   TLegend* lv = new TLegend(0.09, 0.72, 0.27, 0.88);
3451   lv->SetFillColor(kNone);
3452   lv->SetBorderSize(0);
3453
3454   double lv2bottom = 0.62;
3455   TLegend* lv2 = new TLegend(0.82, lv2bottom, 0.93, 0.95, "Centrality");
3456   lv2->SetFillColor(kNone);
3457   lv2->SetBorderSize(0);
3458   cv->Divide(5,1);
3459
3460   for (int n=1; n<=5; n++) {
3461     cv->cd(n);
3462
3463       TH2F *gc = Agreement2DHist(k,n);
3464       if (!gc) {
3465         Error("DrawAgreement()", "Problem finding graph" );
3466         continue;
3467       }
3468
3469       //      gc->GetZaxis()->SetRangeUser(1.e-1, 1000);
3470
3471       gPad->SetLogx(); gPad->SetLogy(); gPad->SetLogz();
3472
3473       if (opt.Contains("lego")) {
3474         gc->SetTitleOffset(0.8, "X");
3475         gc->GetXaxis()->SetLabelSize(0.07);
3476         gc->GetYaxis()->SetLabelSize(0.07);
3477         gc->GetZaxis()->SetLabelSize(0.07);
3478         gc->GetXaxis()->SetTitleSize(0.1);
3479         gc->GetYaxis()->SetTitleSize(0.1);
3480         gc->GetZaxis()->SetTitleSize(0.1);
3481         gc->GetXaxis()->SetNdivisions(8);
3482         gc->GetYaxis()->SetNdivisions(8);
3483
3484         gPad->SetTopMargin(0.2);
3485         gPad->SetLeftMargin(0.12);
3486         gPad->SetRightMargin(0.01);
3487         gc->Draw("lego2");
3488       }
3489       else if (opt.Contains("contour"))
3490         gc->Draw("colz");
3491
3492       if (n==1) {
3493         lv2->AddEntry(gc, Form("%d-%d%%",centlow[k],centhigh[k]),"l,p");
3494       }
3495
3496   } // Fourier moment n
3497
3498   cv->cd();
3499   TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1);
3500   overlay->SetFillStyle(4000); // transparent
3501   overlay->Draw();
3502   overlay->cd();
3503
3504   ltx.SetTextSize(0.07);
3505   const char* str = Form("#left|V_{n#Delta} / (v_{n}^{t} v_{n}^{a})#right|");
3506   ltx.DrawLatex(0.42, 0.92, str);
3507   ltx.DrawLatex(0.91, 0.92, Form("%s", centLabel(k)));
3508
3509   for (int n=1; n<=5; n++) {
3510     cv->cd(n);
3511     ltx.SetTextSize(0.10);
3512     ltx.DrawLatex(0.05, 0.76, Form("n = %d", n));
3513   }
3514   return cv;
3515 }
3516
3517 void AddPrelimStampToCurrentPad(float x1, float y1, float x2, float y2, TString opt)
3518 {
3519   TASImage* logo = new TASImage("../commonfigs/alice-prelim.png"); // 287x323
3520   TPad* logopad = new TPad("logopad", "logo", x1,y1,x2,y2);
3521   if (opt.Contains("border")) {
3522     logopad->SetFillColor(kRed);
3523     logopad->SetBorderSize(2);
3524   }
3525   logopad->Draw();
3526   logopad->cd();
3527   logo->Draw("");
3528   return;
3529 }
3530
3531 int CentBin(double cntLo, double cntHi)
3532 {
3533   int bin = -123;
3534   for (int k=0; k<maxcent; k++) {
3535     if (centlow[k]==cntLo && centhigh[k]==cntHi)
3536       bin = k;
3537   }
3538   if (bin < 0 && 0)
3539     Warning("CentBin","No bin for %.2g - %.2g", cntLo, cntHi);
3540
3541   return bin;
3542 }
3543
3544 int PtBin(double ptLo, double ptHi)
3545 {
3546   int bin = -123;
3547   for (int i=0; i<maxpt; i++) {
3548     if (ptlow[i]==ptLo && pthigh[i]==ptHi)
3549       bin = i;
3550   }
3551   if (bin < 0)
3552     Warning("PtBin","No bin for %.2g - %.2g", ptLo, ptHi);
3553
3554   return bin;
3555 }
3556
3557 double MomCorrection(int i, int j, int k)
3558 {
3559   if(0)cout<<i<<j<<k<<endl;
3560   // For momentum (non)conservation in v1.
3561   // Multiplicitave weight: w = w_t*w_a.
3562   // w_t = ptt - <pt^2>/<pt>
3563   // w_a = pta - <pt^2>/<pt>
3564
3565   //  double c = 1./3 * -0.00185;
3566
3567   //  (ptmean[i]-)* ptmean[j]
3568
3569   double correction =  1./3 * -0.00185 * ptmean[i]* ptmean[j];
3570   return correction;
3571 }
3572
3573 void MakeVnDVsQGraphs(int n1, int n2, int k, const char* region, const char* corrtype, TString opt)
3574 {
3575   // The first two errType constants have points set at the central measured value of
3576   // VnDelta, with statistical and systematic error bars
3577   // respectively. The last two have points set at the upper and lower
3578   // systematic value, with statistical error bars set to be the same
3579   // as for meas_stat.
3580   const int nErrTypes = 4;
3581   enum errType {MEAS_STAT, MEAS_SYS, HI, LO};
3582   //  const char* errTypeStr[] = {"meas_stat","meas_sys","hi","lo"};
3583
3584   // Fill a huge array with all the coefficients, then use it to set the tgraphs.
3585   const int nN = gNMax + 1; // 1-10, don't use 0.
3586   const int nI = maxpt+1;
3587   const int nJ = maxpt+1;
3588   const int nT = nErrTypes;
3589   double vVal[nN][nI][nJ][nT];
3590   double vErr[nN][nI][nJ][nT];
3591   double vMix[nN][nI][nJ][nT]; // Just the acceptance correction unc.
3592
3593   double vSin[nN][nI][nJ][nT]; // <sin n delta phi> for sys. estimation
3594   double vSinErr[nN][nI][nJ][nT];
3595
3596   for (int i=0; i<maxpt; i++) {
3597     for (int j=0; j<=i; j++) {
3598
3599       double esin  = 0;  // residual sine error
3600       double ebw   = 0;  // bin width error \propto n / (# dphi bins)
3601       double ecnt  = 0;  // centrality determination
3602       double emix  = 0;  // acc. correction error
3603       double etrk  = 0;  // tracking resolution
3604       double ev1   = 0;  // v1
3605       double esys  = 0;  // quadrature sum
3606
3607       TH1* hst=0, *hsys=0; // stat and sys error, but same points
3608
3609       hst  = Hist(region, corrtype, i, j, k, "");
3610       hsys = Hist(region, corrtype, i, j, k, "sys");
3611
3612       if (!hst) Error("MakeVnDVsQGraphs()","!hst");
3613       if (!hsys) Error("MakeVnDVsQGraphs()","!hsys");
3614
3615       esin = SineError(i,j,k); // resid_sin / 2;
3616
3617       for (int n=n1; n<=n2; n++) { // Order of coeff.
3618
3619         // VnDelta on measured points, with statistical errors.
3620         double val=0, err=0, resid_sin=0, rs_err=0;
3621         VnDelta(n, *hst,  val, err, "");
3622         VnDelta(n, *hst,  resid_sin, rs_err, "sin");
3623
3624         // VnDelta on high and low sys values. Stat err. not used.
3625         double val_hi=0, xxx=0;
3626         VnDelta(n, *hsys, val_hi, xxx, "hi");
3627         double val_lo=0, yyy=0;
3628         VnDelta(n, *hsys, val_lo, yyy, "lo");
3629
3630         // Momentum conservation correction
3631         if (n==1 && opt.Contains("ptcons")) { // eqs. 3-7, PRL 106, 102301 (2011)
3632           double factor =  ispp ? 0.212 : 0; // 7e-5; <------ trial & error: 7e-5 must be close for 0-20%
3633           val    += factor*ptmean[i]*ptmean[j];
3634           val_hi += factor*ptmean[i]*ptmean[j];
3635           val_lo += factor*ptmean[i]*ptmean[j];
3636         }
3637
3638         int nbins = hst->GetNbinsX();
3639         double mean = (val_hi + val_lo + val) / 3.;
3640         double ms =
3641           (val_hi-mean)*(val_hi-mean) +
3642           (val_lo-mean)*(val_lo-mean) +
3643           (val   -mean)*(val   -mean);
3644         double rms = TMath::Sqrt(ms/2); // sqrt(1/(N-1) sum (x_i - <x>)^2 )
3645
3646
3647         ebw = n/nbins/TMath::Sqrt(12)*val; // This is about 0.008 * n * val
3648         etrk = 0.01*ptmean[j]*val; // see comments in comparisons/VnDcomparison.C
3649         ecnt = 0.01*val;
3650         emix = rms/2; // take half not to double-count stat err. // 0.25*TMath::Abs(val_hi - val_lo);
3651         ev1 = (n==1) ? 0.001*ptmean[i]*ptmean[j]*val : 0.;
3652
3653         esys = TMath::Sqrt(ecnt*ecnt + emix*emix + ebw*ebw + etrk*etrk + ev1*ev1 + esin*esin);
3654
3655         double pointVal = val;
3656         double pointErr = err;
3657
3658         vVal[n][i][j][MEAS_STAT] = pointVal;
3659         vErr[n][i][j][MEAS_STAT] = pointErr;
3660         vVal[n][i][j][MEAS_SYS] = pointVal;
3661         vErr[n][i][j][MEAS_SYS] = esys;
3662         vVal[n][i][j][HI] = pointVal + esys;
3663         vErr[n][i][j][HI] = pointErr;
3664         vVal[n][i][j][LO] = pointVal - esys;
3665         vErr[n][i][j][LO] = pointErr;
3666         vSin[n][i][j][MEAS_STAT] = resid_sin;
3667         vSinErr[n][i][j][MEAS_STAT] = rs_err;
3668         vMix[n][i][j][MEAS_SYS] = emix;
3669
3670       } // n
3671     } // j
3672   } // i
3673
3674   // Now set the graphs from the arrays.
3675   for (int n=n1; n<=n2; n++) {
3676     for (int t=0; t<nErrTypes; t++) {
3677
3678       TString gname = Form("V%dDeltaVsGlobalIndex_cent%dto%d",
3679                            n,centlow[k],centhigh[k]);
3680       if (t==HI)
3681         gname.Append("hi");
3682       else if (t==LO)
3683         gname.Append("lo");
3684       else if (t==MEAS_SYS)
3685         gname.Append("meas_sys");
3686       else
3687         gname.Append("meas_stat");
3688
3689       if (opt.Contains("sine"))
3690         gname.Append("_sine");
3691
3692       if (opt.Contains("ALL"))
3693         gname.Append("_ALL");
3694
3695       TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data());
3696       if (!g) {
3697         g = new TGraphErrors();
3698         g->SetName(gname.Data());
3699       }
3700
3701       TGraphErrors* gmix = 0;
3702       if (t==MEAS_SYS) {
3703         gmix = new TGraphErrors();
3704         gname.Append("_mix");
3705         gmix->SetName(gname.Data());
3706       }
3707
3708       double xer = 0.4; // was 0.5
3709       for (int i=0; i<maxpt; i++) {
3710         for (int j=0; j<=i; j++) {
3711           int q = GlobalIndex(i, j);
3712           g->SetPoint(q, q+0.5, vVal[n][i][j][t]);
3713           g->SetPointError(q, xer, vErr[n][i][j][t]);
3714
3715           if (opt.Contains("sine")) {
3716             g->SetPoint(q, q+0.5, vSin[n][i][j][MEAS_STAT]);
3717             g->SetPointError(q, xer, vSinErr[n][i][j][MEAS_STAT]);
3718           }
3719
3720           if (t==MEAS_SYS) {
3721             gmix->SetPoint(q, q+0.5, vVal[n][i][j][t]);
3722             gmix->SetPointError(q, xer, vMix[n][i][j][t]);
3723           }
3724
3725         }
3726       }
3727
3728       gList->Add(g);
3729       if (gmix) gList->Add(gmix);
3730     }
3731   }
3732
3733   return;
3734 }
3735
3736 TGraphErrors* VnDVsQ(int n, int k, const char* region, const char* corrtype, TString opt)
3737 {
3738
3739   // Graph of pair fourier coefficients VnDelta vs. the linear
3740   // global index q. Arguments:
3741   //
3742   // n, k:   fourier index, centrality bin
3743   // region: "RIDGE", "NSJET" or "ALL"
3744   // corrtype:  "s", "m", "cA", "cB"
3745   //
3746   // Select opt "hi_sys" or "lo_sys" for hi and low pts
3747   // with central stat errors and "sys" for the central points
3748   // with sys. errors. Just like Hist().
3749
3750   TString gname = Form("V%dDeltaVsGlobalIndex_cent%dto%d",
3751                        n,centlow[k],centhigh[k]);
3752
3753   if (opt.Contains("hi"))
3754     gname.Append("hi");
3755   else if (opt.Contains("lo"))
3756     gname.Append("lo");
3757   else if (opt.Contains("sys"))
3758     gname.Append("meas_sys");
3759   else
3760     gname.Append("meas_stat");
3761
3762   if (opt.Contains("sine"))
3763     gname.Append("_sine");
3764
3765   if (opt.Contains("ALL"))
3766     gname.Append("_ALL");
3767
3768   TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data());
3769   if (g) {
3770     if (0)
3771       Info("VnDVsQ()", "Returning preexisting graph");
3772     return g;
3773   }
3774
3775   TString opt1 = (opt.Contains("sine")) ? "sine" : "";
3776   if (opt.Contains("ALL"))
3777     opt1.Append("_ALL");
3778   if (opt.Contains("ptcons"))
3779     opt1.Append("_ptcons");
3780
3781   MakeVnDVsQGraphs(1, gNMax, k, region, corrtype, opt1);
3782   g = (TGraphErrors*) gList->FindObject(gname.Data());
3783   if (!g) {
3784     Error("VnDVsQ()", "No graph %s", gname.Data());
3785     gSystem->Exit(-1);
3786   }
3787   return g;
3788 }
3789
3790 double MeanPt(int i, int j, int k, TString t_or_a, TString opt)
3791 {
3792   // For combined centralities, map k's to the closest available bin
3793   // double cb1[] =  {0, 0, 0,2,2, 1,3,0,    0, 1, 2, 3, 4,  5, 10, 20, 30, 40, 60 };
3794   // double cb2[] =  {10,20,2,5,10,3,5,5,    1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 90 };
3795 // OBJ: TProfile2D      hMeanPtTrgCen1  cen=1 (0.5) : 0 at: 0x102af06e0
3796 //  OBJ: TProfile2D     hMeanPtAssCen1  cen=1 (0.5) : 0 at: 0x108308e60
3797 //  OBJ: TProfile2D     hMean2PtTrgCen1 cen=1 (0.5) : 0 at: 0x108309480
3798 //  OBJ: TProfile2D     hMean2PtAssCen1 cen=1 (0.5) : 0 at: 0x108309ad0
3799 //  OBJ: TProfile2D     hMeanPtTrgCen2  cen=2 (1.5) : 0 at: 0x10830a120
3800 //  OBJ: TProfile2D     hMeanPtAssCen2  cen=2 (1.5) : 0 at: 0x10830a770
3801 //  OBJ: TProfile2D     hMean2PtTrgCen2 cen=2 (1.5) : 0 at: 0x10830adc0
3802 //  OBJ: TProfile2D     hMean2PtAssCen2 cen=2 (1.5) : 0 at: 0x10830b410
3803 //  OBJ: TProfile2D     hMeanPtTrgCen3  cen=3 (2.5) : 0 at: 0x10830ba60
3804 //  OBJ: TProfile2D     hMeanPtAssCen3  cen=3 (2.5) : 0 at: 0x10830c0b0
3805 //  OBJ: TProfile2D     hMean2PtTrgCen3 cen=3 (2.5) : 0 at: 0x10830c700
3806 //  OBJ: TProfile2D     hMean2PtAssCen3 cen=3 (2.5) : 0 at: 0x10830cd50
3807 //  OBJ: TProfile2D     hMeanPtTrgCen4  cen=4 (3.5) : 0 at: 0x10830d3a0
3808 //  OBJ: TProfile2D     hMeanPtAssCen4  cen=4 (3.5) : 0 at: 0x10830d9f0
3809 //  OBJ: TProfile2D     hMean2PtTrgCen4 cen=4 (3.5) : 0 at: 0x10830e040
3810 //  OBJ: TProfile2D     hMean2PtAssCen4 cen=4 (3.5) : 0 at: 0x10830e690
3811 //  OBJ: TProfile2D     hMeanPtTrgCen5  cen=5 (4.5) : 0 at: 0x10830ece0
3812 //  OBJ: TProfile2D     hMeanPtAssCen5  cen=5 (4.5) : 0 at: 0x10830f330
3813 //  OBJ: TProfile2D     hMean2PtTrgCen5 cen=5 (4.5) : 0 at: 0x10830f980
3814 //  OBJ: TProfile2D     hMean2PtAssCen5 cen=5 (4.5) : 0 at: 0x10830ffd0
3815 //  OBJ: TProfile2D     hMeanPtTrgCen6  cen=6 (7.5) : 0 at: 0x108310620
3816 //  OBJ: TProfile2D     hMeanPtAssCen6  cen=6 (7.5) : 0 at: 0x108310c70
3817 //  OBJ: TProfile2D     hMean2PtTrgCen6 cen=6 (7.5) : 0 at: 0x1083112c0
3818 //  OBJ: TProfile2D     hMean2PtAssCen6 cen=6 (7.5) : 0 at: 0x108311910
3819 //  OBJ: TProfile2D     hMeanPtTrgCen7  cen=7 (15.0) : 0 at: 0x108311f60
3820 //  OBJ: TProfile2D     hMeanPtAssCen7  cen=7 (15.0) : 0 at: 0x1083125b0
3821 //  OBJ: TProfile2D     hMean2PtTrgCen7 cen=7 (15.0) : 0 at: 0x108312c00
3822 //  OBJ: TProfile2D     hMean2PtAssCen7 cen=7 (15.0) : 0 at: 0x108313250
3823 //  OBJ: TProfile2D     hMeanPtTrgCen8  cen=8 (25.0) : 0 at: 0x1083138a0
3824 //  OBJ: TProfile2D     hMeanPtAssCen8  cen=8 (25.0) : 0 at: 0x108313ef0
3825 //  OBJ: TProfile2D     hMean2PtTrgCen8 cen=8 (25.0) : 0 at: 0x108314540
3826 //  OBJ: TProfile2D     hMean2PtAssCen8 cen=8 (25.0) : 0 at: 0x108314b90
3827 //  OBJ: TProfile2D     hMeanPtTrgCen9  cen=9 (35.0) : 0 at: 0x1083151e0
3828 //  OBJ: TProfile2D     hMeanPtAssCen9  cen=9 (35.0) : 0 at: 0x108315830
3829 //  OBJ: TProfile2D     hMean2PtTrgCen9 cen=9 (35.0) : 0 at: 0x108315e80
3830 //  OBJ: TProfile2D     hMean2PtAssCen9 cen=9 (35.0) : 0 at: 0x1083164d0
3831 //  OBJ: TProfile2D     hMeanPtTrgCen10 cen=10 (45.0) : 0 at: 0x108316b20
3832 //  OBJ: TProfile2D     hMeanPtAssCen10 cen=10 (45.0) : 0 at: 0x108317170
3833 //  OBJ: TProfile2D     hMean2PtTrgCen10        cen=10 (45.0) : 0 at: 0x1083177c0
3834 //  OBJ: TProfile2D     hMean2PtAssCen10        cen=10 (45.0) : 0 at: 0x108317e10
3835 //  OBJ: TProfile2D     hMeanPtTrgCen11 cen=11 (55.0) : 0 at: 0x108318460
3836 //  OBJ: TProfile2D     hMeanPtAssCen11 cen=11 (55.0) : 0 at: 0x108318ab0
3837 //  OBJ: TProfile2D     hMean2PtTrgCen11        cen=11 (55.0) : 0 at: 0x108319100
3838 //  OBJ: TProfile2D     hMean2PtAssCen11        cen=11 (55.0) : 0 at: 0x108319750
3839 //  OBJ: TProfile2D     hMeanPtTrgCen12 cen=12 (75.0) : 0 at: 0x108319da0
3840 //  OBJ: TProfile2D     hMeanPtAssCen12 cen=12 (75.0) : 0 at: 0x10831a3f0
3841 //  OBJ: TProfile2D     hMean2PtTrgCen12        cen=12 (75.0) : 0 at: 0x10831aa40
3842 //  OBJ: TProfile2D     hMean2PtAssCen12        cen=12 (75.0) : 0 at: 0x10831b090
3843
3844   int k2 = -1;
3845   if (k == CentBin(0,  1))   k2 = 1;
3846   else if (k == CentBin(1,  2))   k2 = 2;
3847   else if (k == CentBin(2,  3))   k2 = 3;
3848   else if (k == CentBin(3,  4))   k2 = 4;
3849   else if (k == CentBin(4,  5))   k2 = 5;
3850   else if (k == CentBin(5, 10))   k2 = 6;
3851   else if (k == CentBin(10,20))   k2 = 7;
3852   else if (k == CentBin(20,30))   k2 = 8;
3853   else if (k == CentBin(30,40))   k2 = 9;
3854   else if (k == CentBin(40,50))   k2 = 10;
3855   else if (k == CentBin(50,60))   k2 = 11;
3856   else if (k == CentBin(60,90))   k2 = 12;
3857   else if (k == CentBin(0, 10))   k2 = 5;    //CentBin(4, 5);
3858   else if (k == CentBin(0, 20))   k2 = 6;    //CentBin(5, 10);
3859   else if (k == CentBin(0,  2))   k2 = 1;    //CentBin(0, 1);
3860   else if (k == CentBin(2,  5))   k2 = 4;    //CentBin(3, 4);
3861   else if (k == CentBin(2, 10))   k2 = 6;    //CentBin(5, 10);
3862   else if (k == CentBin(1,  3))   k2 = 2;    //CentBin(1, 2);
3863   else if (k == CentBin(3,  5))   k2 = 4;    //CentBin(3, 4);
3864   else if (k == CentBin(0,  5))   k2 = 3;    //CentBin(2, 3);
3865
3866   else if (k == CentBin(0, 0))   k2 = 1;    //pp
3867   else {
3868     Error("MeanPt()", "cent bin %d (%d-%d%%) not assigned", k, centlow[k], centhigh[k]);
3869   }
3870
3871   const char* m = opt.Contains("squared") ? "Mean2" : "Mean";
3872
3873   TString prName_t = Form("h%sPtTrgCen%d", m, k2);
3874   TString prName_a = Form("h%sPtAssCen%d", m, k2);
3875   TProfile2D* hmpt_t = (TProfile2D*)fin->Get(prName_t.Data());
3876   TProfile2D* hmpt_a = (TProfile2D*)fin->Get(prName_a.Data());
3877
3878   if (!hmpt_t) {
3879     Error("MeanPt", "Problem finding TProfile %s. Input %d", prName_t.Data(), k);
3880     gSystem->Exit(-1);
3881   }
3882   //  cout << hmpt_t->GetTitle() << " " << centlow[k] << " " << centhigh[k] << endl;
3883
3884   double mptt = hmpt_t->GetBinContent(i+1, j+1);
3885   double mpta = hmpt_a->GetBinContent(i+1, j+1);
3886
3887   if (0) {
3888     // x axis is trig pt,  y axis assc pt
3889     TAxis* ax = hmpt_t->GetXaxis();
3890     TAxis* ay = hmpt_t->GetYaxis();
3891     Printf("trig %.3g - %.3g, assc: %.3g - %.3g, <pt_t,a> %.3g, %.3g",
3892            ax->GetBinLowEdge(i+1),
3893            ax->GetBinUpEdge(i+1),
3894            ay->GetBinLowEdge(j+1),
3895            ay->GetBinUpEdge(j+1),
3896            mptt, mpta );
3897   }
3898
3899   if (t_or_a.Contains("t"))
3900     return mptt;
3901   if (t_or_a.Contains("a"))
3902     return mpta;
3903   else
3904     return -1;
3905 }
3906
3907
3908 int CentColor(double cen1, double cen2)
3909 {
3910   if (cen1==0 && cen2==10)
3911     return kGray+3;
3912   if (cen1==0 && cen2==20)
3913     return kMagenta+3;
3914   if (cen1==0 && cen2==2)
3915     return kBlack;
3916   if (cen1==2 && cen2==10)
3917     return kRed;
3918   if (cen1==10 && cen2==20)
3919     return kOrange-3;
3920   if (cen1==20 && cen2==30)
3921     return kGreen+1;
3922   if (cen1==30 && cen2==40)
3923     return kYellow+2;
3924   if (cen1==40 && cen2==50)
3925     return kAzure;
3926   if (cen1==50 && cen2==60)
3927     return kBlue-1;
3928   if (cen1==60 && cen2==90)
3929     return kViolet;
3930
3931   return kBlack;
3932 }
3933
3934 int CentColorPale(double cen1, double cen2)
3935 {
3936   if (cen1==0 && cen2==10)
3937     return kGray+1;
3938   if (cen1==0 && cen2==20)
3939     return kPink-1; // kPink-6;
3940   if (cen1==0 && cen2==2)
3941     return kGray+1;
3942   if (cen1==2 && cen2==10)
3943     return kRed-9;
3944   if (cen1==10 && cen2==20)
3945     return kOrange-4;
3946   if (cen1==20 && cen2==30)
3947     return kGreen-7;
3948   if (cen1==30 && cen2==40)
3949     return kYellow-8;
3950   if (cen1==40 && cen2==50)
3951     return kAzure-4;
3952   if (cen1==50 && cen2==60)
3953     return kBlue-9;
3954   if (cen1==60 && cen2==90)
3955     return kViolet-9;
3956
3957   return kBlack;
3958 }
3959
3960 int PtColor(double p1, double p2)
3961 {
3962  if (p1==0.25 && p2==0.5)
3963     return kGray+1;
3964  if (p1==0.5 && p2==0.75)
3965     return kYellow+3;
3966  if (p1==0.75 && p2==1.0)
3967    return kPink-6;
3968  if (p1==1.0 && p2==1.5)
3969    return kBlack;
3970  if (p1==1.5 && p2==2.0)
3971    return kRed;
3972  if (p1==2.0 && p2==2.5)
3973    return kOrange-3;
3974  if (p1==2.5 && p2==3)
3975    return kGreen+2;
3976  if (p1==3 && p2==4)
3977    return kBlue;
3978  if (p1==4 && p2==5)
3979    return kViolet + 2;
3980  if (p1==5 && p2==6)
3981    return kViolet+10;
3982  if (p1==4 && p2==6)
3983    return kViolet + 2;
3984  if (p1==6 && p2==8)
3985    return kOrange-7;
3986  if (p1==8 && p2==15)
3987    return kOrange+9;
3988
3989  return kBlack;
3990 }
3991
3992 int PtColorPale(double p1, double p2)
3993 {
3994  if (p1==0.25 && p2==0.5)
3995     return kGray;
3996  if (p1==0.5 && p2==0.75)
3997     return kYellow+1;
3998  if (p1==0.75 && p2==1.0)
3999    return kPink-4;
4000  if (p1==1.0 && p2==1.5)
4001    return kGray+1;
4002  if (p1==1.5 && p2==2.0)
4003    return kRed-7;
4004  if (p1==2.0 && p2==2.5)
4005    return kOrange-4;
4006  if (p1==2.5 && p2==3)
4007    return kSpring-4;
4008  if (p1==3 && p2==4)
4009    return kAzure+6;
4010  if (p1==4 && p2==5)
4011    return kViolet+1;
4012  if (p1==5 && p2==6)
4013    return kViolet+8;
4014  if (p1==4 && p2==6)
4015    return kViolet+1;
4016  if (p1==6 && p2==8)
4017    return kOrange-9;
4018  if (p1==8 && p2==15)
4019    return kOrange+6;
4020
4021  return kGray;
4022 }
4023
4024 int PtaColor(double p1, double p2)
4025 {
4026  if (p1==0.25 && p2==0.5)
4027     return kGray+1;
4028  if (p1==0.5 && p2==0.75)
4029    return kBlack;
4030  if (p1==0.75 && p2==1.0)
4031    return kRed+1;
4032  if (p1==1.0 && p2==1.5)
4033    return kOrange-3;//kBlack;
4034  if (p1==1.5 && p2==2.0)
4035    return kGreen+2;//kRed;
4036  if (p1==2.0 && p2==2.5)
4037    return kAzure-3;//kOrange-3;
4038  if (p1==2.5 && p2==3)
4039    return kViolet + 2;//kBlue;
4040  if (p1==3 && p2==4)
4041    return kViolet + 10;//kBlue;
4042  if (p1==4 && p2==5)
4043    return kMagenta+3; //kViolet + 2;
4044  if (p1==5 && p2==6)
4045    return kViolet;
4046  if (p1==4 && p2==6)
4047    return kMagenta-5;
4048  if (p1==6 && p2==8)
4049    return kGray+2;
4050  if (p1==8 && p2==15)
4051    return kGray+3;
4052
4053  return kBlack;
4054 }
4055
4056 int PtaColorPale(double p1, double p2)
4057 {
4058  if (p1==0.25 && p2==0.5)
4059     return kGray;
4060  if (p1==0.5 && p2==0.75)
4061    return kGray+2;
4062  if (p1==0.75 && p2==1.0)
4063    return kPink-2;
4064  if (p1==1.0 && p2==1.5)
4065    return kOrange-4;
4066  if (p1==1.5 && p2==2.0)
4067    return kGreen-3;
4068  if (p1==2.0 && p2==2.5)
4069    return kAzure-4;
4070  if (p1==2.5 && p2==3)
4071    return kViolet;
4072  if (p1==3 && p2==4)
4073    return kViolet + 2;
4074  if (p1==4 && p2==5)
4075    return kMagenta+3;
4076  if (p1==5 && p2==6)
4077    return kViolet+10;
4078  if (p1==4 && p2==6)
4079    return kMagenta-5;
4080  if (p1==6 && p2==8)
4081    return kGray+2;
4082  if (p1==8 && p2==15)
4083    return kGray+3;
4084
4085  return kGray;
4086 }
4087
4088 void SaveCanvases(TObjArray* canvases, const char* fileName)
4089 {
4090   TFile* f = new TFile(fileName, "recreate");
4091
4092   for (int n=0; n<canvases->GetEntries(); n++) {
4093     TCanvas* c = (TCanvas*)canvases->At(n);
4094     c->Write(c->GetTitle());
4095   }
4096   f->Close();
4097   return;
4098 }
4099
4100 void SaveCanvasesFromFile(const char* rootFile, const char* targetDir, const char* tag, const char* fileType)
4101 {
4102   // Get a list of canvases from rootFile into array, then save each
4103   // to its own file in targetDir/. fileType = "eps", "pdf", "C",
4104   // "png", etc. Not all formats have been tested.
4105   gStyle->SetOptTitle(0);
4106   gStyle->SetOptStat(0);
4107   TString name = "";
4108   TString base(targetDir);
4109   TFile *cFile = new TFile(rootFile, "read");
4110   cout << cFile->GetName() << endl;
4111   TObjArray* cList = GetObjectsFromFile(*cFile, "TCanvas");
4112
4113   for (int n=0; n<cList->GetEntries(); n++) {
4114     TCanvas* c = (TCanvas*)cList->At(n);
4115     if (c) {
4116       name = "";
4117       name = base;
4118       name += TString("/");
4119       name += TString(fileType);
4120       name += TString("/");
4121       name += TString(c->GetTitle());
4122       name += TString(".");
4123       name += TString(fileType);
4124       cout<<name.Data()<<endl;
4125
4126       c->Draw();
4127       c->Modified();
4128       c->Update();
4129       c->SaveAs(name.Data());
4130     }
4131     else
4132       Error("SaveCanvasesFromFile()", "!c");
4133   }
4134
4135   if (1) {
4136     utils.print_pdf(cList, Form("%s/all-figs%s", targetDir, tag), "pdf");
4137   }
4138
4139   return;
4140 }
4141
4142 //void SaveCanvases(TObjArray* canvases, const char* fileName)
4143 //void SaveCanvasesFromFile(const char* rootFile, const char* targetDir, const char* fileType)
4144 TObjArray* GetObjectsFromFile(TFile& file, TString clname, TString dir)
4145 {
4146   file.cd(dir.Data());
4147
4148   TObjArray* objList = new TObjArray();
4149   TIter next(gDirectory->GetListOfKeys());
4150   TKey *key;
4151
4152   while ((key=(TKey*)next())) {
4153     TString className(key->GetClassName());
4154     TString keyName(key->GetName());
4155     if (1)
4156       printf("%10s %20s\n", className.Data(), keyName.Data());
4157
4158     if (className.Contains(clname)) {
4159       objList->Add(gDirectory->Get(keyName.Data()));
4160     }
4161   }
4162
4163   cout << objList->GetEntries() << " objects retrieved from "
4164        << file.GetName()  << "/" << gDirectory->GetName()
4165        << endl;
4166
4167   return objList;
4168 }