PWGCF/Correlations/macros/dphicorrelations/fit.C

   1 #include "TF2.h"
   2 #include "TH2F.h"
   3 #include "TPad.h"
   4 #include "TCanvas.h"
   5 #include "TMath.h"
   6 #include "TGraphErrors.h"
   7 #include "TFile.h"
   8 #include "TLatex.h"
   9 #include "TROOT.h"
  10 #include "TLegend.h"
  11 #include "TLegendEntry.h"
  12 #include "TMultiGraph.h"
  13 #include "TPaveText.h"
  14
  15 void AddPoint(TGraphErrors* graph, Float_t x, Float_t y, Float_t xe, Float_t ye)
  16 {
  17         graph->SetPoint(graph->GetN(), x, y);
  18         graph->SetPointError(graph->GetN() - 1, xe, ye);
  19 }
  20
  21 void DrawLatex(Float_t x, Float_t y, Int_t color, const char* text, Float_t textSize = 0.06)
  22 {
  23   TLatex* latex = new TLatex(x, y, text);
  24   latex->SetNDC();
  25   latex->SetTextSize(textSize);
  26   latex->SetTextColor(color);
  27   latex->Draw();
  28 }
  29
  30 // 0    1    2    3     4     5     6      7      8          9          10         11         12         13         14          15
  31 // norm,dphi,deta,norm2,dphi2,deta2,chi2_1,chi2_2,moment2phi,moment2eta,moment3phi,moment3eta,moment4phi,moment4eta,kurtosisphi,kurtosiseta,
  32 //            16   17   18   19    20    21    22     23     24         25         26         27         28         29         30          31
  33 // second fit:norm,dphi,deta,norm2,dphi2,deta2,chi2_1,chi2_2,moment2phi,moment2eta,moment3phi,moment3eta,moment4phi,moment4eta,kurtosisphi,kurtosiseta,
  34 const Int_t NGraphs = 32;
  35 const Int_t NHists = 6;
  36 TGraphErrors*** graphs = 0;
  37 TList* labelList = 0;
  38
  39 void CreateGraphStructure()
  40 {
  41   graphs = new TGraphErrors**[NGraphs];
  42   for (Int_t i=0; i<NGraphs; i++)
  43   {
  44     graphs[i] = new TGraphErrors*[NHists];
  45     for (Int_t j=0; j<NHists; j++)
  46       graphs[i][j] = new TGraphErrors;
  47   }
  48
  49   labelList = new TList;
  50 }
  51
  52 void WriteGraphs()
  53 {
  54   TFile::Open("graphs.root", "RECREATE");
  55   for (Int_t i=0; i<NGraphs; i++)
  56     for (Int_t j=0; j<NHists; j++)
  57       graphs[i][j]->Write(Form("graph_%d_%d", i, j));
  58
  59   labelList->Write("labelList", TObject::kSingleKey);
  60
  61   gFile->Close();
  62 }
  63
  64 void ReadGraphs(const char* fileName = "graphs.root")
  65 {
  66   CreateGraphStructure();
  67   TFile::Open(fileName);
  68   for (Int_t i=0; i<NGraphs; i++)
  69     for (Int_t j=0; j<NHists; j++)
  70       graphs[i][j] = (TGraphErrors*) gFile->Get(Form("graph_%d_%d", i, j));
  71
  72   labelList = (TList*) gFile->Get("labelList");
  73 }
  74
  75 /*
  76 // old one with 1 Gaussian
  77 Double_t DeltaPhiWidth2DFitFunction(Double_t *x, Double_t *par)
  78 {
  79   // params: 0: gaussian amplitude, 1: phi width, 2: eta width
  80   //         3..bins+2 constants as fct of phi
  81
  82   Float_t phi = x[0];
  83   if (phi < 0)
  84     phi = -phi;
  85   if (phi > TMath::Pi())
  86     phi = TMath::TwoPi() - phi;
  87   Int_t phiBin = (Int_t) (phi / TMath::Pi() * 36);
  88 //   phiBin = 0;
  89
  90   return par[3+phiBin]+par[0]*TMath::Exp(-0.5*((x[0]/par[1])*(x[0]/par[1])+(x[1]/par[2])*(x[1]/par[2])));
  91 }*/
  92
  93 const Double_t k1OverSqrtTwoPi = 1.0 / TMath::Sqrt(TMath::TwoPi());
  94
  95 Double_t DeltaPhiWidth2DFitFunction(Double_t *x, Double_t *par)
  96 {
  97   // params: 0: gaussian amplitude, 1: phi width I, 2: eta width I
  98   //         3: fraction for first Gaussian, 4: phi width II, 5: eta width II
  99   //         6..bins+5 constants as fct of phi
 100
 101   Float_t phi = x[0];
 102   if (phi < 0)
 103     phi = -phi;
 104   if (phi > TMath::Pi())
 105     phi = TMath::TwoPi() - phi;
 106   Int_t phiBin = (Int_t) (phi / TMath::Pi() * 36);
 107 //   phiBin = 0;
 108
 109   return par[6+phiBin] + par[0]*(
 110       par[3]/TMath::TwoPi()/par[1]/par[2] *
 111         TMath::Exp(-0.5*((x[0]/par[1])*(x[0]/par[1])+(x[1]/par[2])*(x[1]/par[2])))
 112       + (1-par[3])/TMath::TwoPi()/par[4]/par[5]
 113         * TMath::Exp(-0.5*((x[0]/par[4])*(x[0]/par[4])+(x[1]/par[5])*(x[1]/par[5])))
 114     );
 115 }
 116
 117 void SubtractEtaGap(TH2* hist, Float_t etaLimit, Float_t outerLimit, Bool_t scale, Bool_t drawEtaGapDist = kFALSE)
 118 {
 119   TString histName(hist->GetName());
 120
 121   TH1D* etaGap = hist->ProjectionX(histName + "_1", TMath::Max(1, hist->GetYaxis()->FindBin(-outerLimit + 0.01)), hist->GetYaxis()->FindBin(-etaLimit - 0.01));
 122   Int_t etaBins = hist->GetYaxis()->FindBin(-etaLimit - 0.01) - TMath::Max(1, hist->GetYaxis()->FindBin(-outerLimit + 0.01)) + 1;
 123
 124   TH1D* tracksTmp = hist->ProjectionX(histName + "_2", hist->GetYaxis()->FindBin(etaLimit + 0.01), TMath::Min(hist->GetYaxis()->GetNbins(), hist->GetYaxis()->FindBin(outerLimit - 0.01)));
 125   etaBins += TMath::Min(hist->GetYaxis()->GetNbins(), hist->GetYaxis()->FindBin(outerLimit - 0.01)) - hist->GetYaxis()->FindBin(etaLimit + 0.01) + 1;
 126
 127   etaGap->Add(tracksTmp);
 128
 129   // get per bin result
 130   etaGap->Scale(1.0 / etaBins);
 131
 132   if (drawEtaGapDist)
 133   {
 134     TH1D* centralRegion = hist->ProjectionX(histName + "_3", hist->GetYaxis()->FindBin(-etaLimit + 0.01), hist->GetYaxis()->FindBin(etaLimit - 0.01));
 135
 136     centralRegion->Scale(1.0 / (hist->GetYaxis()->FindBin(etaLimit - 0.01) - hist->GetYaxis()->FindBin(-etaLimit + 0.01) + 1));
 137
 138     TCanvas* c = new TCanvas("SubtractEtaGap", "SubtractEtaGap", 800, 800);
 139     centralRegion->SetStats(0);
 140     centralRegion->Draw();
 141     etaGap->DrawCopy("SAME")->SetLineColor(2);
 142     c->SaveAs("etagap_proj.eps");
 143   }
 144
 145 //   new TCanvas; etaGap->DrawCopy();
 146
 147   TH2* histTmp2D = (TH2*) hist->Clone("histTmp2D");
 148   histTmp2D->Reset();
 149
 150   for (Int_t xbin=1; xbin<=histTmp2D->GetNbinsX(); xbin++)
 151     for (Int_t y=1; y<=histTmp2D->GetNbinsY(); y++)
 152       histTmp2D->SetBinContent(xbin, y, etaGap->GetBinContent(xbin));
 153
 154   if (scale)
 155   {
 156     // mixed event does not reproduce away-side perfectly
 157     // --> extract scaling factor on the away-side from ratios of eta gap and central region
 158     TH1D* centralRegion = hist->ProjectionX(histName + "_3", hist->GetYaxis()->FindBin(-etaLimit + 0.01), hist->GetYaxis()->FindBin(etaLimit - 0.01));
 159     etaBins = hist->GetYaxis()->FindBin(etaLimit - 0.01) - hist->GetYaxis()->FindBin(-etaLimit + 0.01) + 1;
 160     centralRegion->Scale(1.0 / etaBins);
 161
 162 //     new TCanvas; centralRegion->DrawCopy(); etaGap->SetLineColor(2); etaGap->DrawCopy("SAME");
 163     centralRegion->Divide(etaGap);
 164 //     new TCanvas; centralRegion->Draw();
 165     centralRegion->Fit("pol0", "0", "", TMath::Pi() - 1, TMath::Pi() + 1);
 166     Float_t scalingFactor = centralRegion->GetFunction("pol0")->GetParameter(0);
 167     Printf("  scalingFactor = %f", scalingFactor);
 168     histTmp2D->Scale(scalingFactor);
 169   }
 170
 171 //   new TCanvas; hist->DrawCopy("SURF1");
 172
 173   hist->Add(histTmp2D, -1);
 174 }
 175
 176 /*void FitDeltaPhiEtaGap2D(TH2* hist, Bool_t scale,  TCanvas* canvas, Int_t canvasPos, TGraphErrors** width, Float_t x, Float_t yPosChi2, TGraphErrors* chi2_1, TGraphErrors* chi2_2)
 177 {
 178   Float_t etaLimit = 1.0;
 179   Float_t outerLimit = 1.6;
 180
 181   SubtractEtaGap(hist, etaLimit, outerLimit, scale);
 182
 183 //   new TCanvas; hist->DrawCopy("SURF1");
 184
 185   hist->GetYaxis()->SetRangeUser(-outerLimit+0.01, outerLimit-0.01);
 186   hist->GetXaxis()->SetRangeUser(-TMath::Pi() / 2 + 0.01, TMath::Pi() * 0.5 - 0.01);
 187
 188   canvas->cd(canvasPos);
 189   hist->SetStats(0);
 190   hist->DrawCopy("SURF1");
 191
 192   Float_t min = hist->GetMinimum();
 193   Float_t max = hist->GetMaximum();
 194
 195   // ranges are to exclude eta gap region from fit
 196   TF2* func = new TF2("func", "[0]+[1]*exp(-0.5*((x/[2])**2+(y/[3])**2))", -0.5 * TMath::Pi(), 0.5 * TMath::Pi(), -outerLimit, outerLimit);
 197   func->SetParameters(0, 1, 0.3, 0.3);
 198   func->SetParLimits(1, 0, 10);
 199   func->SetParLimits(2, 0.05, 1);
 200   func->SetParLimits(3, 0.05, 1);
 201   func->FixParameter(0, 0);
 202
 203 //   TF2* func = new TF2("func", DeltaPhiWidth2DFitFunction, -TMath::Pi() / 2, TMath::Pi() * 1.5, -1, 1, 4);
 204 //   func->SetParameters(1, 0.3, 0.3, 0);
 205 //   func->SetParLimits(0, 0, 10);
 206 //   func->SetParLimits(1, 0.1, 10);
 207 //   func->SetParLimits(2, 0.1, 10);
 208
 209   hist->Fit(func, "0R", "");
 210   hist->Fit(func, "I0R", "");
 211
 212 //   func->SetRange(-0.5 * TMath::Pi(), 0.5 * TMath::Pi(), -outerLimit, outerLimit);
 213
 214   canvas->cd(canvasPos + 1);
 215   TH2* funcHist = (TH2*) hist->Clone("funcHist");
 216   funcHist->Reset();
 217   funcHist->Add(func);
 218   funcHist->SetMinimum(min);
 219   funcHist->SetMaximum(max);
 220   funcHist->Draw("SURF1");
 221
 222   canvas->cd(canvasPos + 2);
 223   hist->Add(func, -1);
 224   hist->SetMinimum(min);
 225   hist->SetMaximum(max);
 226   hist->DrawCopy("SURF1");
 227
 228   width[0]->SetPoint(width[0]->GetN(), x, TMath::Abs(func->GetParameter(2)));
 229   width[0]->SetPointError(width[0]->GetN()-1, 0, func->GetParError(2));
 230
 231   width[1]->SetPoint(width[1]->GetN(), x, TMath::Abs(func->GetParameter(3)));
 232   width[1]->SetPointError(width[1]->GetN()-1, 0, func->GetParError(3));
 233
 234   Float_t chi2 = 0;
 235   Int_t ndf = 0;
 236   for (Int_t i=hist->GetXaxis()->FindBin(-0.8); i<=hist->GetXaxis()->FindBin(0.8); i++)
 237     for (Int_t j=hist->GetYaxis()->FindBin(-0.8); j<=hist->GetYaxis()->FindBin(0.8); j++)
 238     {
 239       if (hist->GetBinError(i, j) > 0)
 240       {
 241         chi2 += TMath::Power(hist->GetBinContent(i, j) / hist->GetBinError(i, j), 2);
 242         ndf++;
 243       }
 244     }
 245   ndf -= func->GetNumberFreeParameters();
 246
 247   printf("#chi^{2}/ndf = %.1f/%d = %.1f  ", func->GetChisquare(), func->GetNDF(), func->GetChisquare() / func->GetNDF());
 248   Printf("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf);
 249
 250   DrawLatex(0.5, yPosChi2, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", func->GetChisquare(), func->GetNDF(), func->GetChisquare() / func->GetNDF()));
 251   DrawLatex(0.5, yPosChi2 - 0.05, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf));
 252
 253   chi2_1->SetPoint(chi2_1->GetN(), x, func->GetChisquare() / func->GetNDF());
 254   chi2_2->SetPoint(chi2_2->GetN(), x, chi2 / ndf);
 255 }
 256 */
 257
 258 void CalculateMomentsKurtosis(Float_t momentCalcLimit, TH1* proj, Int_t graphIDStart, Int_t histId, Float_t x)
 259 {
 260   //   momentCalcLimit = 0.6;
 261   for (Int_t n=2; n <= 4; n++)
 262   {
 263     Float_t moment = 0;
 264     Float_t sum = 0;
 265     for (Int_t bin = proj->GetXaxis()->FindBin(-momentCalcLimit+0.01); bin <= proj->GetXaxis()->FindBin(momentCalcLimit-0.01); bin++)
 266     {
 267       moment += proj->GetBinContent(bin) * TMath::Power(proj->GetMean() - proj->GetXaxis()->GetBinCenter(bin), n);
 268       sum += proj->GetBinContent(bin);
 269     }
 270     moment /= sum;
 271
 272     Float_t error = 0;
 273     for (Int_t bin = proj->GetXaxis()->FindBin(-momentCalcLimit+0.01); bin <= proj->GetXaxis()->FindBin(momentCalcLimit-0.01); bin++)
 274     {
 275       error += proj->GetBinError(bin) * proj->GetBinError(bin) *
 276         TMath::Power(TMath::Power(proj->GetMean() - proj->GetXaxis()->GetBinCenter(bin), n) / sum
 277           - moment / sum, 2);
 278     }
 279
 280     AddPoint(graphs[graphIDStart+(n-2)*2][histId], x, moment, 0, TMath::Sqrt(error));
 281   //     Printf("%d %f +- %f <-> %f +- %f", n, moment, error, projx2->GetRMS() * projx2->GetRMS(), 2 * projx2->GetRMSError() * projx2->GetRMSError());
 282     }
 283
 284   proj->GetXaxis()->SetRangeUser(-momentCalcLimit, momentCalcLimit);
 285   AddPoint(graphs[graphIDStart+6][histId], x, proj->GetKurtosis(1), 0, proj->GetKurtosis(11));
 286 }
 287
 288 void FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Float_t x, Float_t yPosChi2, Bool_t quick, Int_t histId, Int_t limits, Bool_t twoTrack)
 289 {
 290   Float_t etaLimit = 1.0;
 291   Float_t outerLimit = 1.59;
 292   Float_t sigmaFitLimit = 0.1 - limits * 0.05;
 293   Float_t etaFitUpperLimit = 0.8;
 294   Float_t initSigma = 0.6;
 295   if (histId == 2) // pp
 296   {
 297     etaFitUpperLimit = 0.6;
 298     initSigma = 0.4;
 299   }
 300
 301   hist->GetYaxis()->SetRangeUser(-outerLimit+0.01, outerLimit-0.01);
 302   hist->GetXaxis()->SetRangeUser(-TMath::Pi() / 2 + 0.01, TMath::Pi() * 0.5 - 0.01);
 303
 304   Float_t mean = hist->Integral(hist->GetYaxis()->FindBin(-TMath::Pi() / 2), hist->GetYaxis()->FindBin(TMath::Pi() / 2), hist->GetYaxis()->FindBin(1.0), hist->GetYaxis()->FindBin(outerLimit)) / (hist->GetYaxis()->FindBin(TMath::Pi() / 2) - hist->GetYaxis()->FindBin(-TMath::Pi() / 2)) / (hist->GetYaxis()->FindBin(outerLimit) - hist->GetYaxis()->FindBin(1.0) + 1);
 305 //   Printf("%f", mean);
 306
 307   // sums
 308   TGraphErrors* sumSummary = new TGraphErrors;
 309   TGraphErrors* phiWidthSummary = new TGraphErrors;
 310   TGraphErrors* etaWidthSummary = new TGraphErrors;
 311
 312   canvas->cd(canvasPos++);
 313   hist->SetStats(0);
 314   hist->DrawCopy("SURF1");
 315   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), hist->Integral());
 316
 317   Float_t min = hist->GetMinimum();
 318   Float_t max = hist->GetMaximum();
 319
 320   Int_t bins = hist->GetNbinsX() / 2 / 2;
 321
 322   TF2* func = new TF2("func", DeltaPhiWidth2DFitFunction, -TMath::Pi() / 2, TMath::Pi() * 0.5, -outerLimit, outerLimit, bins+6);
 323   func->SetParameters(1, 0.3, 0.3, 0.25, initSigma, initSigma);
 324   for (Int_t i=6; i<bins+6; i++)
 325     func->SetParameter(i, mean);
 326
 327   func->SetParLimits(0, 0, 10);
 328   func->SetParLimits(1, sigmaFitLimit, 1);
 329   func->SetParLimits(2, sigmaFitLimit, etaFitUpperLimit);
 330   func->SetParLimits(3, 0.1, 0.9);
 331   func->SetParLimits(4, sigmaFitLimit, 1);
 332   func->SetParLimits(5, sigmaFitLimit, etaFitUpperLimit);
 333
 334   Int_t fitResult = hist->Fit(func, "0R", "");
 335   Printf("Fit result: %d", fitResult);
 336
 337   // if both parameters are within 1%, refit with 1 Gaussian only
 338   if (TMath::Abs(1.0 - func->GetParameter(1) / func->GetParameter(4)) < 0.01 && TMath::Abs(1.0 - func->GetParameter(2) / func->GetParameter(5)) < 0.01)
 339   {
 340     Printf("Parameters within 1%%. Refitting with 1 Gaussian...");
 341
 342     func->SetParLimits(3, 1, 1);
 343     func->FixParameter(3, 1);
 344     func->FixParameter(4, sigmaFitLimit);
 345     func->FixParameter(5, sigmaFitLimit);
 346
 347     fitResult = hist->Fit(func, "0R", "");
 348     Printf("Fit result: %d", fitResult);
 349   }
 350
 351   Int_t first = 1;
 352   Int_t second = 4;
 353   if (func->GetParameter(1) < func->GetParameter(4))
 354   {
 355     first = 4;
 356     second = 1;
 357   }
 358   //dphi
 359   AddPoint(graphs[1+16][histId], x, TMath::Abs(func->GetParameter(first)), 0, func->GetParError(first));
 360   AddPoint(graphs[4+16][histId], x, TMath::Abs(func->GetParameter(second)), 0, func->GetParError(second));
 361
 362   //deta
 363   AddPoint(graphs[2+16][histId], x, TMath::Abs(func->GetParameter(first+1)), 0, func->GetParError(first+1));
 364   AddPoint(graphs[5+16][histId], x, TMath::Abs(func->GetParameter(second+1)), 0, func->GetParError(second+1));
 365
 366   // norm
 367   AddPoint(graphs[0+16][histId], x, func->GetParameter(0), 0, func->GetParError(0));
 368   if (first < second)
 369     AddPoint(graphs[3+16][histId], x, func->GetParameter(3), 0, func->GetParError(3));
 370   else
 371     AddPoint(graphs[3+16][histId], x, 1.0 - func->GetParameter(3), 0, func->GetParError(3));
 372
 373   //   hist->Fit(func, "0RI", "");
 374
 375   //   width[0]->SetPoint(width[0]->GetN(), x, TMath::Abs(func->GetParameter(1)));
 376   //   width[0]->SetPointError(width[0]->GetN()-1, 0, func->GetParError(1));
 377   //
 378   //   width[1]->SetPoint(width[1]->GetN(), x, TMath::Abs(func->GetParameter(2)));
 379   //   width[1]->SetPointError(width[1]->GetN()-1, 0, func->GetParError(2));
 380
 381   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), func->GetParameter(1), 0, func->GetParError(1));
 382   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), func->GetParameter(4), 0, func->GetParError(4));
 383
 384   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), func->GetParameter(2), 0, func->GetParError(2));
 385   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), func->GetParameter(5), 0, func->GetParError(5));
 386
 387   canvas->cd(canvasPos++);
 388   TH2* funcHist = (TH2*) hist->Clone("funcHist");
 389   funcHist->Reset();
 390   funcHist->Add(func);
 391   funcHist->SetMinimum(min);
 392   funcHist->SetMaximum(max);
 393   funcHist->Draw("SURF1");
 394   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), funcHist->Integral());
 395
 396   canvas->cd(canvasPos++);
 397   TH2* residuals = (TH2*) hist->Clone("residuals");
 398   residuals->Add(func, -1);
 399   residuals->SetMinimum(-(max - min) / 2);
 400   residuals->SetMaximum((max - min) / 2);
 401   residuals->Draw("SURF1");
 402
 403   Float_t chi2 = 0;
 404   Int_t ndf = 0;
 405   for (Int_t i=hist->GetXaxis()->FindBin(-0.8); i<=hist->GetXaxis()->FindBin(0.8); i++)
 406     for (Int_t j=hist->GetYaxis()->FindBin(-etaLimit); j<=hist->GetYaxis()->FindBin(etaLimit); j++)
 407     {
 408       if (residuals->GetBinError(i, j) > 0)
 409       {
 410         chi2 += TMath::Power(residuals->GetBinContent(i, j) / residuals->GetBinError(i, j), 2);
 411         ndf++;
 412       }
 413     }
 414   ndf -= func->GetNumberFreeParameters();
 415
 416   if (func->GetNDF() > 0)
 417   {
 418     printf("#chi^{2}/ndf = %.1f/%d = %.1f  ", func->GetChisquare(), func->GetNDF(), func->GetChisquare() / func->GetNDF());
 419     DrawLatex(0.2, yPosChi2, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", func->GetChisquare(), func->GetNDF(), func->GetChisquare() / func->GetNDF()));
 420   }
 421   if (ndf)
 422   {
 423     Printf("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf);
 424     DrawLatex(0.2, yPosChi2 - 0.05, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf));
 425   }
 426
 427   // draw gaussian only
 428   TF2* funcClone = new TF2("funcClone", DeltaPhiWidth2DFitFunction, -TMath::Pi() / 2, TMath::Pi() * 0.5, -outerLimit, outerLimit, bins+6);
 429   for (Int_t i=0; i<6; i++)
 430     funcClone->SetParameter(i, func->GetParameter(i));
 431   for (Int_t i=6; i<bins+6; i++)
 432     funcClone->SetParameter(i, 0);
 433 //   funcClone->Print();
 434   canvas->cd(canvasPos++);
 435   funcHist = (TH2*) hist->Clone("funcHistb");
 436   funcHist->Reset();
 437   funcHist->Add(funcClone);
 438   funcHist->SetMinimum(-(max - min) / 2);
 439   funcHist->SetMaximum((max - min) / 2);
 440   funcHist->Draw("SURF1");
 441
 442   // eta gap subtraction
 443   canvas->cd(canvasPos++);
 444   func->SetParameter(0, 0);
 445   TH2* subtractFlow = (TH2*) hist->Clone("subtractFlow");
 446   subtractFlow->Add(func, -1);
 447   subtractFlow->SetMinimum(-(max - min) / 2);
 448   subtractFlow->SetMaximum((max - min) / 2);
 449   subtractFlow->DrawCopy("SURF1");
 450   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), subtractFlow->Integral());
 451
 452   canvas->cd(canvasPos++);
 453   SubtractEtaGap(hist, etaLimit, outerLimit, kFALSE);
 454   hist->SetMinimum(-(max - min) / 2);
 455   hist->SetMaximum((max - min) / 2);
 456   hist->DrawCopy("SURF1");
 457   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), hist->Integral());
 458
 459   if (!quick)
 460   {
 461     canvas->cd(canvasPos++);
 462     TH2* difference = (TH2*) hist->Clone("difference");
 463     difference->Add(subtractFlow, -1);
 464     difference->SetMinimum(-(max - min) / 2);
 465     difference->SetMaximum((max - min) / 2);
 466     difference->DrawCopy("SURF1");
 467
 468     canvas->cd(canvasPos++);
 469     TF2* func2 = new TF2("func2", "[0]+[1]*exp(-0.5*((x/[2])**2+(y/[3])**2))", -0.5 * TMath::Pi(), 0.5 * TMath::Pi(), -outerLimit, outerLimit);
 470     func2->SetParameters(0, 1, 0.3, 0.3);
 471     func2->SetParLimits(1, 0, 10);
 472     func2->SetParLimits(2, sigmaFitLimit, 1);
 473     func2->SetParLimits(3, sigmaFitLimit, 1);
 474     func2->FixParameter(0, 0);
 475
 476     hist->Fit(func2, "0R", "");
 477     //   hist->Fit(func2, "I0R", "");
 478
 479     AddPoint(phiWidthSummary, phiWidthSummary->GetN(), func2->GetParameter(2), 0, func2->GetParError(2));
 480     AddPoint(etaWidthSummary, etaWidthSummary->GetN(), func2->GetParameter(3), 0, func2->GetParError(3));
 481
 482     TH2* funcHist2 = (TH2*) hist->Clone("funcHist2");
 483     funcHist2->Reset();
 484     funcHist2->Add(func2);
 485     funcHist2->SetMinimum(-(max - min) / 2);
 486     funcHist2->SetMaximum((max - min) / 2);
 487     funcHist2->Draw("SURF1");
 488     sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), funcHist2->Integral());
 489
 490     if (func2->GetNDF() > 0)
 491     {
 492       Printf("#chi^{2}/ndf = %.1f/%d = %.1f", func2->GetChisquare(), func2->GetNDF(), func2->GetChisquare() / func2->GetNDF());
 493       DrawLatex(0.2, yPosChi2, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", func2->GetChisquare(), func2->GetNDF(), func2->GetChisquare() / func2->GetNDF()));
 494     }
 495
 496     canvas->cd(canvasPos++);
 497     TH2* residuals2 = (TH2*) hist->Clone("residuals");
 498     residuals2->Add(funcHist2, -1);
 499     residuals2->SetMinimum(-(max - min) / 2);
 500     residuals2->SetMaximum((max - min) / 2);
 501     residuals2->Draw("SURF1");
 502   }
 503
 504   Float_t momentFitLimit = 0.8;
 505   TH1* projx2 = hist->ProjectionX(Form("%s_projx2", hist->GetName()), hist->GetYaxis()->FindBin(-momentFitLimit+0.01), hist->GetYaxis()->FindBin(momentFitLimit-0.01));
 506   projx2->GetXaxis()->SetRangeUser(-1, 1);
 507
 508   TH1* projy2 = hist->ProjectionY(Form("%s_projy2", hist->GetName()), hist->GetXaxis()->FindBin(-momentFitLimit+0.01), hist->GetXaxis()->FindBin(momentFitLimit-0.01));
 509   projy2->GetXaxis()->SetRangeUser(-1, 1);
 510
 511   CalculateMomentsKurtosis(momentFitLimit, projx2, 8, histId, x);
 512   CalculateMomentsKurtosis(momentFitLimit, projy2, 9, histId, x);
 513
 514   //   return;
 515
 516   TH1* projx1 = subtractFlow->ProjectionX(Form("%s_projx1", hist->GetName()), hist->GetYaxis()->FindBin(-0.8), hist->GetYaxis()->FindBin(0.8));
 517   TH1* projy1 = subtractFlow->ProjectionY(Form("%s_projy1", hist->GetName()), hist->GetXaxis()->FindBin(-0.8), hist->GetXaxis()->FindBin(0.8));
 518
 519   CalculateMomentsKurtosis(momentFitLimit, projx1, 8+16, histId, x);
 520   CalculateMomentsKurtosis(momentFitLimit, projy1, 9+16, histId, x);
 521
 522 //   TF1* twoGauss = new TF1("twoGauss", "gaus(0)+gaus(3)", -2, 2);
 523 //   twoGauss->SetParameters(1, 0, 0.3, 1, 0, 0.6);
 524 //   twoGauss->FixParameter(1, 0);
 525 //   twoGauss->FixParameter(4, 0);
 526 //   twoGauss->SetLineColor(4);
 527 //   projx1->Fit("twoGauss", "I+", "SAME");
 528
 529   canvas->cd(canvasPos++);
 530   projx1->Draw();
 531   projx1->Fit("gaus", "I");
 532
 533   projx2->SetLineColor(2);
 534   projx2->Draw("SAME");
 535   projx2->Fit("gaus", "I+", "SAME");
 536   projx2->GetFunction("gaus")->SetLineColor(2);
 537
 538   canvas->cd(canvasPos++);
 539   projy1->Draw();
 540   projy1->Fit("gaus", "I");
 541
 542   projy2->SetLineColor(2);
 543   projy2->Draw("SAME");
 544   projy2->Fit("gaus", "I+", "SAME");
 545   projy2->GetFunction("gaus")->SetLineColor(2);
 546
 547   // 1d fit (lots of params)
 548   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), projx2->GetFunction("gaus")->GetParameter(2), 0, projx2->GetFunction("gaus")->GetParError(2));
 549
 550   // 1d fit (eta gap subtraction)
 551   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), projx1->GetFunction("gaus")->GetParameter(2), 0, projx1->GetFunction("gaus")->GetParError(2));
 552
 553   // 1d fit (lots of params)
 554   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), projy2->GetFunction("gaus")->GetParameter(2), 0, projy2->GetFunction("gaus")->GetParError(2));
 555
 556   // 1d fit (eta gap subtraction)
 557   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), projy1->GetFunction("gaus")->GetParameter(2), 0, projy1->GetFunction("gaus")->GetParError(2));
 558
 559   // set errors large for bins potentially affected by two-track effects
 560   if (0 && twoTrack)
 561   {
 562     Printf("NOTE : Skipping bins at (0, 0)");
 563 //     for (Int_t binx = hist->GetXaxis()->FindBin(-0.25); binx <= hist->GetXaxis()->FindBin(0.25); binx++)
 564     for (Int_t binx = hist->GetXaxis()->FindBin(-0.01); binx <= hist->GetXaxis()->FindBin(0.01); binx++)
 565       for (Int_t biny = hist->GetYaxis()->FindBin(-0.01); biny <= hist->GetYaxis()->FindBin(0.01); biny++)
 566       {
 567 //      hist->SetBinContent(binx, biny,  0);
 568         hist->SetBinError(binx, biny,  1e5);
 569       }
 570   }
 571
 572   // 2d fit with two gaussians
 573   canvas->cd(canvasPos++);
 574 //   TF2* func3 = new TF2("func3", "[0]+[1]*exp(-0.5*((x/[2])**2+(y/[3])**2))+[4]*exp(-0.5*((x/[5])**2+(y/[6])**2))", -0.5 * TMath::Pi(), 0.5 * TMath::Pi(), -outerLimit, outerLimit);
 575 //   func3->SetParameters(0, 1, 0.3, 0.3, 1, 0.6, 0.6);
 576 //   func3->SetParLimits(4, 0, 10);
 577   Float_t etaFitLimit = outerLimit;
 578 //   Float_t etaFitLimit = 0.5;
 579   TF2* func3 = new TF2("func3", "[0]+[1]*([4]/TMath::TwoPi()/[2]/[3]*exp(-0.5*((x/[2])**2+(y/[3])**2))+(1-[4])/TMath::TwoPi()/[5]/[6]*exp(-0.5*((x/[5])**2+(y/[6])**2)))", -0.5 * TMath::Pi(), 0.5 * TMath::Pi(), -etaFitLimit, etaFitLimit);
 580   func3->SetParameters(0, 1, 0.3, 0.3, 0.25, initSigma, initSigma);
 581   func3->SetParLimits(4, 0.1, 0.9);
 582
 583   func3->SetParLimits(1, 0, 10);
 584   func3->SetParLimits(2, sigmaFitLimit, 1);
 585   func3->SetParLimits(3, sigmaFitLimit, etaFitUpperLimit);
 586   func3->SetParLimits(5, sigmaFitLimit, 1);
 587   func3->SetParLimits(6, sigmaFitLimit, etaFitUpperLimit);
 588   func3->FixParameter(0, 0);
 589
 590   if (0 && histId == 0)
 591   {
 592     // central --> go to 1 Gaussian only
 593     func3->SetParLimits(4, 1, 1);
 594     func3->FixParameter(4, 1);
 595     func3->FixParameter(5, sigmaFitLimit);
 596     func3->FixParameter(6, sigmaFitLimit);
 597   }
 598
 599   // set errors 20% of the value for bins potentially affected by two-track effects
 600 //   hist->SetBinError(hist->GetXaxis()->FindBin(0.0001),  hist->GetYaxis()->FindBin(0.0001),  0.1 * hist->GetBinContent(hist->GetXaxis()->FindBin(0.0001),  hist->GetYaxis()->FindBin(0.0001)));
 601 //   hist->SetBinError(hist->GetXaxis()->FindBin(0.0001),  hist->GetYaxis()->FindBin(-0.0001), 0.1 * hist->GetBinContent(hist->GetXaxis()->FindBin(0.0001),  hist->GetYaxis()->FindBin(-0.0001)));
 602 //   hist->SetBinError(hist->GetXaxis()->FindBin(-0.0001), hist->GetYaxis()->FindBin(0.0001),  0.1 * hist->GetBinContent(hist->GetXaxis()->FindBin(-0.0001), hist->GetYaxis()->FindBin(0.0001)));
 603 //   hist->SetBinError(hist->GetXaxis()->FindBin(-0.0001), hist->GetYaxis()->FindBin(-0.0001), 0.1 * hist->GetBinContent(hist->GetXaxis()->FindBin(-0.0001), hist->GetYaxis()->FindBin(-0.0001)));
 604
 605   fitResult = hist->Fit(func3, "0R", "");
 606   Printf("Fit result: %d", fitResult);
 607
 608   // if both parameters are within 1%, refit with 1 Gaussian only
 609   if (TMath::Abs(1.0 - func3->GetParameter(2) / func3->GetParameter(5)) < 0.01 && TMath::Abs(1.0 - func3->GetParameter(3) / func3->GetParameter(6)) < 0.01)
 610   {
 611     Printf("Parameters within 1%%. Refitting with 1 Gaussian...");
 612
 613     func3->SetParLimits(4, 1, 1);
 614     func3->FixParameter(4, 1);
 615     func3->FixParameter(5, sigmaFitLimit);
 616     func3->FixParameter(6, sigmaFitLimit);
 617
 618     fitResult = hist->Fit(func3, "0R", "");
 619     Printf("Fit result: %d", fitResult);
 620   }
 621
 622 //   hist->Fit(func3, "I0R", "");
 623
 624   first = 2;
 625   second = 5;
 626   if (func3->GetParameter(2) < func3->GetParameter(5))
 627   {
 628     first = 5;
 629     second = 2;
 630   }
 631   //dphi
 632   AddPoint(graphs[1][histId], x, TMath::Abs(func3->GetParameter(first)), 0, func3->GetParError(first));
 633   AddPoint(graphs[4][histId], x, TMath::Abs(func3->GetParameter(second)), 0, func3->GetParError(second));
 634
 635   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), TMath::Abs(func3->GetParameter(first)), 0, func3->GetParError(first));
 636   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), TMath::Abs(func3->GetParameter(second)), 0, func3->GetParError(second));
 637
 638   //deta
 639   AddPoint(graphs[2][histId], x, TMath::Abs(func3->GetParameter(first+1)), 0, func3->GetParError(first+1));
 640   AddPoint(graphs[5][histId], x, TMath::Abs(func3->GetParameter(second+1)), 0, func3->GetParError(second+1));
 641
 642   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), TMath::Abs(func3->GetParameter(first+1)), 0, func3->GetParError(first+1));
 643   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), TMath::Abs(func3->GetParameter(second+1)), 0, func3->GetParError(second+1));
 644
 645   // norm
 646   AddPoint(graphs[0][histId], x, func3->GetParameter(1), 0, func3->GetParError(1));
 647   if (first < second)
 648     AddPoint(graphs[3][histId], x, func3->GetParameter(4), 0, func3->GetParError(4));
 649   else
 650     AddPoint(graphs[3][histId], x, 1.0 - func3->GetParameter(4), 0, func3->GetParError(4));
 651
 652   TH2* funcHist3 = (TH2*) hist->Clone("funcHist3");
 653   funcHist3->Reset();
 654   funcHist3->Add(func3);
 655   funcHist3->SetMinimum(-(max - min) / 2);
 656   funcHist3->SetMaximum((max - min) / 2);
 657   funcHist3->Draw("SURF1");
 658   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), funcHist3->Integral());
 659
 660   canvas->cd(canvasPos++);
 661   TH2* residuals3 = (TH2*) hist->Clone("residuals");
 662   residuals3->Add(funcHist3, -1);
 663   residuals3->SetMinimum(-(max - min) / 2);
 664   residuals3->SetMaximum((max - min) / 2);
 665   residuals3->Draw("SURF1");
 666
 667   chi2 = 0;
 668   ndf = 0;
 669   for (Int_t i=hist->GetXaxis()->FindBin(-0.8); i<=hist->GetXaxis()->FindBin(0.8); i++)
 670     for (Int_t j=hist->GetYaxis()->FindBin(-etaLimit); j<=hist->GetYaxis()->FindBin(etaLimit); j++)
 671     {
 672       if (residuals3->GetBinError(i, j) > 0)
 673       {
 674         chi2 += TMath::Power(residuals3->GetBinContent(i, j) / residuals3->GetBinError(i, j), 2);
 675         ndf++;
 676       }
 677     }
 678   ndf -= func3->GetNumberFreeParameters();
 679
 680   if (func3->GetNDF() > 0)
 681   {
 682     printf("#chi^{2}/ndf = %.1f/%d = %.1f  ", func3->GetChisquare(), func3->GetNDF(), func3->GetChisquare() / func3->GetNDF());
 683     DrawLatex(0.2, yPosChi2, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", func3->GetChisquare(), func3->GetNDF(), func3->GetChisquare() / func3->GetNDF()));
 684     AddPoint(graphs[6][histId], x, func3->GetChisquare() / func3->GetNDF(), 0, 0);
 685   }
 686   if (ndf)
 687   {
 688     Printf("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf);
 689     DrawLatex(0.2, yPosChi2 - 0.05, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf));
 690     AddPoint(graphs[7][histId], x, chi2 / ndf, 0, 0);
 691   }
 692
 693   canvas->cd(canvasPos++);
 694   phiWidthSummary->SetMarkerStyle(20);
 695   phiWidthSummary->Draw("AP");
 696   gPad->SetGridy();
 697
 698   etaWidthSummary->SetMarkerStyle(21);
 699   etaWidthSummary->SetLineColor(2);
 700   etaWidthSummary->SetMarkerColor(2);
 701   etaWidthSummary->Draw("PSAME");
 702
 703   phiWidthSummary->GetYaxis()->SetRangeUser(0, 0.9);
 704
 705   canvas->cd(canvasPos++);
 706   sumSummary->Draw("*A");
 707   gPad->SetGridy();
 708 }
 709
 710 void AnalyzeDeltaPhiEtaGap2D(const char* fileName, Int_t method = 1)
 711 {
 712   gROOT->SetBatch(kTRUE);
 713   if (!gROOT->IsBatch())
 714   {
 715     Printf("Not in batch mode. Exiting!");
 716     return;
 717   }
 718
 719   CreateGraphStructure();
 720
 721   TFile::Open(fileName);
 722
 723   Int_t maxLeadingPt = 3;
 724   Int_t maxAssocPt = 6;
 725
 726   for (Int_t histId = 0; histId < NHists; histId++)
 727   {
 728     for (Int_t i=0; i<maxLeadingPt; i++)
 729     {
 730       for (Int_t j=1; j<maxAssocPt; j++)
 731       {
 732 //      if (i != 1 || j != 2)
 733 //        continue;
 734
 735         TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
 736         if (!hist1)
 737           continue;
 738
 739         if (hist1->GetEntries() < 1e4)
 740         {
 741           Printf("Only %f entries. Skipping...", hist1->GetEntries());
 742           continue;
 743         }
 744
 745         TCanvas* canvas = new TCanvas(Form("DeltaPhi_%d_%d_%d_%d", histId, i, j, method), Form("DeltaPhi_%d_%d_%d_%d", histId, i, j, method), 1400, 1100);
 746         canvas->Divide(5, 3);
 747
 748         Printf("\n\n>>> %d %d %d", i, j, histId);
 749
 750         for (Int_t k=1; k<=3; k++)
 751         {
 752           canvas->cd(3 * j + k);
 753           gPad->SetLeftMargin(0.15);
 754           gPad->SetBottomMargin(0.2);
 755           gPad->SetTopMargin(0.01);
 756           gPad->SetRightMargin(0.01);
 757         }
 758
 759         if (histId == 0)
 760           labelList->Add(new TObjString(hist1->GetTitle()));
 761
 762 //      hist1->Rebin2D(2, 1);
 763
 764         Float_t xPos = j*8+i;
 765
 766 /*      if (method == 0)
 767           FitDeltaPhiEtaGap2D((TH2*) hist1, kFALSE, canvas, 1, width, xPos, 0.9, graphs[6][histId], graphs[7][histId]);
 768         else*/
 769         FitDeltaPhi2DOneFunction((TH2*) hist1, canvas, 1, xPos, 0.9, kTRUE, histId, (i > 0) ? 1 : 0, (j <= 2 && histId != 2));
 770
 771 //      canvas->SaveAs(Form("%s.png", canvas->GetName()));
 772 //      delete canvas;
 773 //      break;
 774 // return;
 775       }
 776
 777 //       break;
 778     }
 779
 780 //     break;
 781   }
 782
 783   WriteGraphs();
 784 }
 785
 786 void AnalyzeDeltaPhiEtaGap2DExample(const char* fileName, Int_t i, Int_t j, Int_t histId, Bool_t drawDetails = kFALSE)
 787 {
 788   CreateGraphStructure();
 789
 790   TFile::Open(fileName);
 791
 792   TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
 793   if (!hist1)
 794     return;
 795
 796   Printf("Entries: %f %s", hist1->GetEntries(), hist1->GetTitle());
 797
 798   hist1->SetTitle("");
 799   hist1->GetYaxis()->SetRangeUser(-1.79, 1.79);
 800   hist1->GetXaxis()->SetTitleOffset(1.5);
 801   hist1->GetYaxis()->SetTitleOffset(2);
 802   hist1->SetStats(kFALSE);
 803
 804   if (hist1->GetEntries() < 1e4)
 805     return;
 806
 807   TCanvas* canvas = new TCanvas(Form("DeltaPhi_%d_%d_%d_%d", histId, i, j, 1), Form("DeltaPhi_%d_%d_%d_%d", histId, i, j, 1), 1400, 1100);
 808   canvas->Divide(5, 3);
 809
 810   Float_t xPos = j*8+i;
 811
 812   FitDeltaPhi2DOneFunction((TH2*) hist1, canvas, 1, xPos, 0.9, kFALSE, histId, (i > 0) ? 1 : 0, (j <= 2 && histId != 2));
 813
 814   if (!drawDetails)
 815     return;
 816
 817   TVirtualPad* pad = canvas->cd(6);
 818   TCanvas* c = new TCanvas("c", "c", 800, 800);
 819   pad->SetPad(0.05, 0, 1, 1);
 820   pad->Draw();
 821   c->SaveAs("fit_subtracted.eps");
 822
 823   pad = canvas->cd(12);
 824   c = new TCanvas("c2", "c2", 800, 800);
 825   pad->SetPad(0.05, 0, 1, 1);
 826   pad->Draw();
 827   c->SaveAs("fit_fit1.eps");
 828
 829   pad = canvas->cd(13);
 830   c = new TCanvas("c3", "c3", 800, 800);
 831   gPad->SetLeftMargin(0.15);
 832   pad->SetPad(0.05, 0, 1, 1);
 833   pad->Draw();
 834   c->SaveAs("fit_residual1.eps");
 835
 836   pad = canvas->cd(4);
 837   c = new TCanvas("c4", "c4", 800, 800);
 838   pad->SetPad(0.05, 0, 1, 1);
 839   pad->Draw();
 840   c->SaveAs("fit_fit2.eps");
 841
 842   pad = canvas->cd(3);
 843   c = new TCanvas("c5", "c5", 800, 800);
 844   pad->SetPad(0.05, 0, 1, 1);
 845   pad->Draw();
 846   c->SaveAs("fit_residual2.eps");
 847
 848   pad = canvas->cd(7);
 849   c = new TCanvas("c6", "c6", 800, 800);
 850   pad->SetPad(0.05, 0, 1, 1);
 851   pad->Draw();
 852   c->SaveAs("fit_differences.eps");
 853 }
 854
 855 void SqrtAll(Int_t nHists, TGraphErrors** graph)
 856 {
 857   for (Int_t histId = 0; histId < nHists; histId++)
 858   {
 859     for (Int_t i=0; i<graph[histId]->GetN(); i++)
 860     {
 861       if (graph[histId]->GetY()[i] > 0)
 862       {
 863         graph[histId]->GetEY()[i] *= 0.5 / TMath::Sqrt(graph[histId]->GetY()[i]);
 864         graph[histId]->GetY()[i] = TMath::Sqrt(graph[histId]->GetY()[i]);
 865       }
 866     }
 867   }
 868 }
 869
 870 Int_t marker[6] = { 24, 25, 26, 30, 27, 28 };
 871 Int_t marker2[6] = { 20, 21, 22, 29, 33, 34 };
 872 const char* labels[6] = { "0-5%", "60-90%", "pp", "20-30%", "10-20%", "40-60%" };
 873
 874 void Draw(const char* canvasName, Int_t nHists, TGraphErrors** graph, TGraphErrors** graph2 = 0, Bool_t normalize = kFALSE, Float_t min = 0, Float_t max = 0)
 875 {
 876   Int_t colors[6] = { 1, 2, 4, 6, 7, 8 };
 877   Bool_t found = kTRUE;
 878   TCanvas* c1 = (TCanvas*) gROOT->GetListOfCanvases()->FindObject(canvasName);
 879   if (!c1)
 880   {
 881     c1 = new TCanvas(canvasName, canvasName, 800, 600);
 882 //     colors[0] = 1;
 883     found = kFALSE;
 884   }
 885   c1->cd();
 886
 887   if (normalize && graph2)
 888   {
 889     // relative norms
 890     for (Int_t histId = 0; histId < nHists; histId++)
 891     {
 892       for (Int_t i=0; i<graph[histId]->GetN(); i++)
 893       {
 894         Float_t sum = graph[histId]->GetY()[i] + graph2[histId]->GetY()[i];
 895         graph[histId]->GetY()[i] /= sum;
 896         graph2[histId]->GetY()[i] /= sum;
 897         graph[histId]->GetEY()[i] /= sum;
 898         graph2[histId]->GetEY()[i] /= sum;
 899       }
 900     }
 901   }
 902
 903   for (Int_t histId = 0; histId < nHists; histId++)
 904   {
 905     graph[histId]->SetMarkerStyle((found) ? marker2[histId] : marker[histId]);
 906     graph[histId]->SetMarkerColor(colors[histId]);
 907     graph[histId]->SetLineColor(colors[histId]);
 908     if (normalize)
 909       graph[histId]->GetYaxis()->SetRangeUser(0, 1);
 910     if (max > min)
 911       graph[histId]->GetYaxis()->SetRangeUser(min, max);
 912     graph[histId]->GetXaxis()->SetRangeUser(7, 38);
 913     graph[histId]->DrawClone((histId == 0 && !found) ? "AP" : "PSAME");
 914     DrawLatex(0.7, 0.8 - 0.05 * histId, colors[histId], labels[histId]);
 915   }
 916
 917   if (graph2)
 918   {
 919     for (Int_t histId = 0; histId < nHists; histId++)
 920     {
 921       graph2[histId]->SetMarkerStyle(marker2[histId]);
 922       graph2[histId]->SetMarkerColor(colors[histId]);
 923       graph2[histId]->SetLineColor(colors[histId]);
 924       graph2[histId]->DrawClone("PSAME");
 925     }
 926   }
 927
 928   gPad->SetGridx();
 929   gPad->SetGridy();
 930   c1->SaveAs(Form("%s.png", canvasName));
 931 }
 932
 933 void CalculateRMS(Int_t nHists, TGraphErrors** graph, TGraphErrors** graph2, TGraphErrors** weight)
 934 {
 935   // rms
 936   for (Int_t histId = 0; histId < nHists; histId++)
 937   {
 938     for (Int_t i=0; i<graph[histId]->GetN(); i++)
 939     {
 940       Float_t rms = graph[histId]->GetY()[i] * weight[histId]->GetY()[i] + graph2[histId]->GetY()[i] * (1.0 - weight[histId]->GetY()[i]);
 941       graph[histId]->GetY()[i] = rms;
 942       // error**2 = weight**2 * error_sigma**2 + (weight-1)**2 * error_sigma2**2 + (sigma1 + sigma2)**2 * error_weight**2
 943       // TODO this neglects some correlations
 944       graph[histId]->GetEY()[i] = TMath::Sqrt(TMath::Power(weight[histId]->GetY()[i] * graph[histId]->GetEY()[i], 2) +
 945         TMath::Power((weight[histId]->GetY()[i] - 1) * graph2[histId]->GetEY()[i], 2) + TMath::Power((graph[histId]->GetY()[i] + graph2[histId]->GetY()[i]) * weight[histId]->GetEY()[i], 2));
 946     }
 947   }
 948 }
 949
 950 void DrawRMS(const char* canvasName, Int_t nHists, TGraphErrors** graph, TGraphErrors** graph2, TGraphErrors** weight, Float_t min = 0, Float_t max = 0)
 951 {
 952 Int_t colors[6] = { 2, 2, 4, 6, 7, 8 };
 953   Bool_t found = kTRUE;
 954   TCanvas* c1 = (TCanvas*) gROOT->GetListOfCanvases()->FindObject(canvasName);
 955   if (!c1)
 956   {
 957     c1 = new TCanvas(canvasName, canvasName, 800, 600);
 958     colors[0] = 1;
 959     found = kFALSE;
 960   }
 961   c1->cd();
 962
 963   for (Int_t histId = 0; histId < nHists; histId++)
 964   {
 965     graph[histId]->SetMarkerStyle(marker[histId]);
 966     graph[histId]->SetMarkerColor(colors[histId]);
 967     graph[histId]->SetLineColor(colors[histId]);
 968     if (max > min)
 969       graph[histId]->GetYaxis()->SetRangeUser(min, max);
 970     graph[histId]->DrawClone((histId == 0 && !found) ? "AP" : "PSAME");
 971     DrawLatex(0.7, 0.8 - 0.05 * histId, colors[histId], labels[histId]);
 972
 973     graph2[histId]->SetMarkerStyle(marker2[histId]);
 974     graph2[histId]->SetMarkerColor(colors[histId]);
 975     graph2[histId]->SetLineColor(colors[histId]);
 976     graph2[histId]->DrawClone("PSAME");
 977   }
 978
 979   CalculateRMS(nHists, graph, graph2, weight);
 980
 981   // rms
 982   for (Int_t histId = 0; histId < nHists; histId++)
 983     graph[histId]->DrawClone("*SAME");
 984
 985   gPad->SetGridx();
 986   gPad->SetGridy();
 987   c1->SaveAs(Form("%s.png", canvasName));
 988 }
 989
 990 TGraphErrors* GetCentralityGraph(Int_t nHists, TGraphErrors** graph, Int_t point)
 991 {
 992   TGraphErrors* graphcentrality = new TGraphErrors;
 993
 994   Float_t centralityAxisMapping[] = { 5, 75, 100, 30, 15, 50 };
 995   Float_t centralityAxisMappingE[] = { 5, 15, 0, 10, 5, 10 };
 996
 997   for (Int_t histId = 0; histId < nHists; histId++)
 998   {
 999     Int_t j = 0;
1000     while (j < graph[histId]->GetN() && TMath::Abs(graph[0]->GetX()[point] - graph[histId]->GetX()[j]) > 0.1)
1001       j++;
1002     if (j == graph[histId]->GetN())
1003       continue;
1004
1005 //       Printf("%d %d %f %f %f %f", point, histId, graph[0]->GetX()[point], graph[histId]->GetX()[j], graph[0]->GetY()[point], graph[histId]->GetY()[j]);
1006     if (TMath::Abs(graph[0]->GetX()[point] - graph[histId]->GetX()[j]) > 0.1)
1007       continue;
1008     AddPoint(graphcentrality, centralityAxisMapping[histId], graph[histId]->GetY()[j], centralityAxisMappingE[histId], graph[histId]->GetEY()[j]);
1009 //       Printf("%f %f", graph[histId]->GetY()[i], graph[histId]->GetEY()[i]);
1010 //       graphcentrality->Print();
1011   }
1012
1013   graphcentrality->Sort();
1014   return graphcentrality;
1015 }
1016
1017 void PrepareGraphs(Int_t nHists, TGraphErrors** graph, TGraphErrors** systematic, TMultiGraph** multiGraph, TMultiGraph** multiGraphSyst)
1018 {
1019   Int_t colors[] =  { 1, 2, 3, 4, 6, 7, 8, 9 };
1020   Int_t markers[] = { 20, 21, 22, 23, 24, 26, 27, 28 };
1021   Int_t fillStyle[] = { 3001, 3002, 3003, 3004, 3005, 3006, 3007, 3008, 3009 };
1022   Int_t count = 0;
1023
1024   Int_t beginWith = 0;
1025
1026   if (*multiGraph == 0)
1027     *multiGraph = new TMultiGraph;
1028   if (*multiGraphSyst == 0)
1029     *multiGraphSyst = new TMultiGraph;
1030
1031   for (Int_t i=beginWith; i<graph[0]->GetN(); i++)
1032   {
1033     Int_t xPos = (Int_t) graph[0]->GetX()[i];
1034     if (xPos == 24 || xPos == 33 || xPos == 34 || xPos == 41 || xPos == 26)
1035       continue;
1036     Printf("%d %i", i, xPos);
1037
1038     TGraphErrors* graphcentrality = GetCentralityGraph(nHists, graph, i);
1039
1040     if (systematic)
1041     {
1042       TGraphErrors* graphsystematics = GetCentralityGraph(nHists, systematic, i);
1043
1044       if (graphcentrality->GetN() != graphsystematics->GetN())
1045       {
1046         Printf("Different number of points %d %d", graphcentrality->GetN(), graphsystematics->GetN());
1047         return;
1048       }
1049
1050       for (Int_t j=0; j<=graphsystematics->GetN(); j++)
1051       {
1052         graphsystematics->GetEY()[j] = TMath::Abs(graphcentrality->GetY()[j] - graphsystematics->GetY()[j]) / 2;
1053         graphsystematics->GetY()[j]  = (graphcentrality->GetY()[j] + graphsystematics->GetY()[j]) / 2;
1054         graphsystematics->GetEX()[j] = 1;
1055       }
1056
1057 //       graphsystematics->SetFillColor(kGray);
1058       graphsystematics->SetFillColor(colors[count]);
1059 //       graphsystematics->SetFillStyle(1001);
1060       graphsystematics->SetFillStyle(fillStyle[count]);
1061       graphsystematics->SetMarkerStyle(0);
1062       graphsystematics->SetLineColor(0);
1063       (*multiGraphSyst)->Add(graphsystematics, "2");
1064     }
1065
1066     graphcentrality->SetMarkerStyle(markers[count]);
1067     graphcentrality->SetMarkerColor(colors[count]);
1068     graphcentrality->SetLineColor(colors[count]);
1069
1070     (*multiGraph)->Add(graphcentrality);
1071
1072     TString label = (labelList && i < labelList->GetEntries()) ? labelList->At(i)->GetName() : "";
1073     if (label.Length() > 0)
1074     {
1075       TObjArray* tokens = label.Tokenize("-");
1076       label.Form("%s-%s", tokens->At(0)->GetName(),tokens->At(1)->GetName());
1077       label.ReplaceAll(".00", "");
1078       label.ReplaceAll(".0", "");
1079     }
1080     graphcentrality->SetTitle(label);
1081
1082     count++;
1083   }
1084 }
1085
1086 void DrawCentrality(const char* canvasName, Int_t nHists, TGraphErrors** graph, Float_t min = 0, Float_t max = 0, const char* yLabel = "", TGraphErrors** systematic = 0, TGraphErrors** graph2 = 0, TGraphErrors** systematic2 = 0)
1087 {
1088   Bool_t found = kTRUE;
1089   TCanvas* c1 = (TCanvas*) gROOT->GetListOfCanvases()->FindObject(canvasName);
1090   if (!c1)
1091   {
1092     c1 = new TCanvas(canvasName, canvasName, 800, 600);
1093     found = kFALSE;
1094   }
1095   c1->cd();
1096
1097   TMultiGraph* multiGraph = 0;
1098   TMultiGraph* multiGraphSyst = 0;
1099
1100   PrepareGraphs(nHists, graph, systematic, &multiGraph, &multiGraphSyst);
1101   TLegend* legend = new TLegend(0.65, 0.65, 0.95, 0.95);
1102   legend->SetFillColor(0);
1103   for (Int_t i=0; i<multiGraph->GetListOfGraphs()->GetEntries(); i++)
1104     legend->AddEntry(multiGraph->GetListOfGraphs()->At(i), 0, "PL");
1105
1106   if (graph2)
1107   {
1108     TMultiGraph* multiGraph2 = 0;
1109     PrepareGraphs(nHists, graph2, systematic2, &multiGraph2, &multiGraphSyst);
1110     for (Int_t i=0; i<multiGraph2->GetListOfGraphs()->GetEntries(); i++)
1111       ((TGraph*)multiGraph2->GetListOfGraphs()->At(i))->SetLineWidth(2);
1112
1113     while (multiGraph2->GetListOfGraphs()->GetEntries() > multiGraph->GetListOfGraphs()->GetEntries())
1114       multiGraph2->GetListOfGraphs()->RemoveAt(multiGraph->GetListOfGraphs()->GetEntries());
1115
1116
1117     multiGraphSyst->Add(multiGraph2, "LX");
1118   }
1119
1120   multiGraphSyst->Add(multiGraph, (found) ? "LX" : "P");
1121
1122   TString drawString("A");
1123   multiGraphSyst->Draw(drawString);
1124   multiGraphSyst->GetXaxis()->SetTitle("Centrality");
1125   multiGraphSyst->GetYaxis()->SetTitle(yLabel);
1126   if (max > min)
1127     multiGraphSyst->GetYaxis()->SetRangeUser(min, max);
1128
1129   legend->Draw();
1130
1131   gPad->SetGridx();
1132   gPad->SetGridy();
1133   c1->SaveAs(Form("results/%s.png", canvasName));
1134 //   c1->SaveAs(Form("%s.eps", canvasName));
1135 }
1136
1137 void DrawResults(const char* fileName = "graphs.root")
1138 {
1139   ReadGraphs(fileName);
1140
1141 //   return;
1142
1143   Int_t nHists = 6;
1144
1145   if (1)
1146   {
1147     Draw("width_phi", nHists, graphs[1], graphs[4], kFALSE, 0, 0.8);
1148     Draw("width_eta", nHists, graphs[2], graphs[5], kFALSE, 0, 0.8);
1149     Draw("norm", nHists, graphs[0], graphs[3], kFALSE, 0, 1);
1150     Draw("chi2", nHists, graphs[6], graphs[7], kFALSE, 0.5, 2.5);
1151     DrawRMS("width_phi_rms", nHists, graphs[1], graphs[4], graphs[3], 0, 0.8);
1152     DrawRMS("width_eta_rms", nHists, graphs[2], graphs[5], graphs[3], 0, 0.8);
1153   }
1154
1155   Draw("moment2", nHists, graphs[8], graphs[9], kFALSE, 0, 0.3);
1156   Draw("moment3", nHists, graphs[10], graphs[11], kFALSE, -0.05, 0.05);
1157   Draw("moment4", nHists, graphs[12], graphs[13], kFALSE, -0.1, 0.2);
1158   Draw("kurtosis", nHists, graphs[14], graphs[15], kFALSE, -2, 2);
1159   Draw("kurtosisphi", nHists, graphs[14], 0, kFALSE, -2, 2);
1160   Draw("kurtosiseta", nHists, graphs[15], 0, kFALSE, -2, 2);
1161 }
1162
1163 void DrawResultsCentrality(const char* fileName = "graphs.root")
1164 {
1165   ReadGraphs(fileName);
1166
1167   Int_t nHists = 6;
1168
1169   if (1)
1170   {
1171     DrawCentrality("norm", nHists, graphs[0], 0, 0.2, "N (a.u.)", graphs[0+16]);
1172
1173     DrawCentrality("width_phi1_centrality", nHists, graphs[1], 0, 0.8, "#sigma_{#phi, 1} (rad.)", graphs[1+16]);
1174 //     return;
1175     DrawCentrality("width_phi2_centrality", nHists, graphs[4], 0, 0.8, "#sigma_{#phi, 2} (rad.)", graphs[4+16]);
1176     DrawCentrality("width_eta1_centrality", nHists, graphs[2], 0, 0.8, "#sigma_{#eta, 1} (rad.)", graphs[2+16]);
1177     DrawCentrality("width_eta2_centrality", nHists, graphs[5], 0, 0.8, "#sigma_{#eta, 2} (rad.)", graphs[5+16]);
1178
1179     CalculateRMS(nHists, graphs[1], graphs[4], graphs[3]);
1180     CalculateRMS(nHists, graphs[2], graphs[5], graphs[3]);
1181
1182     CalculateRMS(nHists, graphs[1+16], graphs[4+16], graphs[3+16]);
1183     CalculateRMS(nHists, graphs[2+16], graphs[5+16], graphs[3+16]);
1184
1185     DrawCentrality("phi_rms_centrality", nHists, graphs[1], 0, 0.8, "#sigma_{#phi} (fit) (rad.)", graphs[1+16]);
1186     DrawCentrality("eta_rms_centrality", nHists, graphs[2], 0, 0.8, "#sigma_{#eta} (fit)", graphs[2+16]);
1187
1188     DrawCentrality("chi2_1", nHists, graphs[6], 0.5, 2.5, "#chi^2/ndf (full region)");
1189     DrawCentrality("chi2_2", nHists, graphs[7], 0.5, 2.5, "#chi^2/ndf (peak region)");
1190
1191     // sqrt(moment2) = sigma
1192     SqrtAll(nHists, graphs[8]);
1193     SqrtAll(nHists, graphs[9]);
1194     SqrtAll(nHists, graphs[8+16]);
1195     SqrtAll(nHists, graphs[9+16]);
1196     DrawCentrality("sigma_phi", nHists, graphs[8], 0, 0.8, "#sigma_{#phi} (rad.)", graphs[8+16]);
1197     DrawCentrality("sigma_eta", nHists, graphs[9], 0, 0.8, "#sigma_{#eta}", graphs[9+16]);
1198   }
1199
1200   DrawCentrality("kurtosisphi_centrality", nHists, graphs[14], -2, 4, "Kurtosis #phi", graphs[14+16]);
1201   DrawCentrality("kurtosiseta_centrality", nHists, graphs[15], -2, 4, "Kurtosis #eta", graphs[15+16]);
1202 }
1203
1204 void DrawRMSCentrality(const char* fileNameData, const char* fileNameHijing)
1205 {
1206   Int_t nHists = 6;
1207
1208   ReadGraphs(fileNameHijing);
1209   if (1)
1210   {
1211     Int_t n=0;
1212     while (n++ < graphs[0][0]->GetN())
1213     {
1214   //     Printf("%d %f", n, graphs[0][0]->GetX()[n]); continue;
1215       if (graphs[0][0]->GetX()[n] == 42)
1216       {
1217         for (Int_t i=0; i<NGraphs; i++)
1218           graphs[i][0]->RemovePoint(n);
1219         break;
1220       }
1221     }
1222   }
1223 //   return;
1224
1225   CalculateRMS(nHists, graphs[1], graphs[4], graphs[3]);
1226   CalculateRMS(nHists, graphs[2], graphs[5], graphs[3]);
1227
1228   CalculateRMS(nHists, graphs[1+16], graphs[4+16], graphs[3+16]);
1229   CalculateRMS(nHists, graphs[2+16], graphs[5+16], graphs[3+16]);
1230
1231 //   DrawCentrality("phi_rms_centrality", nHists, graphs[1], 0, 0.8, "#sigma_{#phi} (fit) (rad.)");
1232 //   return;
1233
1234   TGraphErrors*** graphs1 = graphs;
1235
1236   ReadGraphs(fileNameData);
1237
1238   //Remove high trigger pT:18.26,34,42
1239   if (1)
1240   {
1241     Int_t n=0;
1242     while (n++ < graphs[0][0]->GetN())
1243     {
1244   //     Printf("%d %f", n, graphs[0][0]->GetX()[n]); continue;
1245       if (graphs[0][0]->GetX()[n] == 18 || graphs[0][0]->GetX()[n] == 10) //  || graphs[0][0]->GetX()[n] == 34 || graphs[0][0]->GetX()[n] == 42
1246       {
1247         for (Int_t i=0; i<NGraphs; i++)
1248         {
1249           graphs[i][0]->RemovePoint(n);
1250         }
1251         n--;
1252       }
1253     }
1254   }
1255
1256   CalculateRMS(nHists, graphs[1], graphs[4], graphs[3]);
1257   CalculateRMS(nHists, graphs[2], graphs[5], graphs[3]);
1258
1259   CalculateRMS(nHists, graphs[1+16], graphs[4+16], graphs[3+16]);
1260   CalculateRMS(nHists, graphs[2+16], graphs[5+16], graphs[3+16]);
1261
1262   DrawCentrality("phi_rms_centrality_mc", nHists, graphs[1], 0, 0.8, "#sigma_{#phi} (fit) (rad.)", graphs[1+16], graphs1[1]);
1263   DrawCentrality("eta_rms_centrality_mc", nHists, graphs[2], 0, 0.8, "#sigma_{#eta} (fit)", graphs[2+16], graphs1[2]);
1264
1265   DrawCentrality("kurtosisphi_centrality_mc", nHists, graphs[14], -2, 4, "Kurtosis #phi", graphs[14+16], graphs1[14]);
1266   DrawCentrality("kurtosiseta_centrality_mc", nHists, graphs[15], -2, 4, "Kurtosis #eta", graphs[15+16], graphs1[15]);
1267 }
1268
1269 void CompareRMS(const char* graphFileName1, const char* graphFileName2)
1270 {
1271   ReadGraphs(graphFileName1);
1272
1273   Int_t nHists = 3;
1274
1275   DrawRMS("width_phi_rms1", nHists, graphs[1], graphs[4], graphs[3], 0, 0.8);
1276   DrawRMS("width_eta_rms1", nHists, graphs[2], graphs[5], graphs[3], 0, 0.8);
1277
1278   Draw("width_phi_rms_compare", nHists, graphs[1], 0, kFALSE, 0, 0.8);
1279   Draw("width_eta_rms_compare", nHists, graphs[2], 0, kFALSE, 0, 0.8);
1280
1281   ReadGraphs(graphFileName2);
1282   DrawRMS("width_phi_rms2", nHists, graphs[1], graphs[4], graphs[3], 0, 0.8);
1283   DrawRMS("width_eta_rms2", nHists, graphs[2], graphs[5], graphs[3], 0, 0.8);
1284
1285   Draw("width_phi_rms_compare", nHists, graphs[1], 0, kFALSE, 0, 0.8);
1286   Draw("width_eta_rms_compare", nHists, graphs[2], 0, kFALSE, 0, 0.8);
1287 }
1288
1289 void CompareEtaPhi(const char* graphFileName1)
1290 {
1291   ReadGraphs(graphFileName1);
1292
1293   Int_t nHists = 6;
1294
1295   CalculateRMS(nHists, graphs[1], graphs[4], graphs[3]);
1296   CalculateRMS(nHists, graphs[2], graphs[5], graphs[3]);
1297
1298   CalculateRMS(nHists, graphs[1+16], graphs[4+16], graphs[3+16]);
1299   CalculateRMS(nHists, graphs[2+16], graphs[5+16], graphs[3+16]);
1300
1301   DrawCentrality("rms_centrality", nHists, graphs[1], 0, 0.8, "#sigma_{#phi} (fit) (rad.) / #sigma_{#eta}", graphs[1+16], graphs[2], graphs[2+16]);
1302   DrawCentrality("rms_centrality_nosyst", nHists, graphs[1], 0, 0.8, "#sigma_{#phi} (fit) (rad.) / #sigma_{#eta}", 0, graphs[2], 0);
1303 }
1304
1305 void DrawExamples(const char* histFileName, const char* graphFileName, Bool_t drawFunc = kTRUE)
1306 {
1307   Float_t etaLimit = 1.0;
1308   Float_t outerLimit = 1.79;
1309   Float_t projectLimit = 0.8;
1310
1311   ReadGraphs(graphFileName);
1312
1313   TFile::Open(histFileName);
1314
1315   Int_t j=1;
1316
1317   Int_t exColors[] = { 1, 2, 4, 6 };
1318
1319   for (Int_t i=0; i<6; i++)
1320   {
1321     TCanvas* c = new TCanvas(Form("c_%d", i), Form("c_%d", i), 1000, 600);
1322     c->Divide(2, 1);
1323     Int_t nHists = 3;
1324     for (Int_t histId = 0; histId < nHists; histId++)
1325     {
1326       TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
1327       if (!hist)
1328         continue;
1329
1330       Float_t xPos = j*8+i;
1331       Int_t n = 0;
1332       while ( n < graphs[0][histId]->GetN() && graphs[0][histId]->GetX()[n] != xPos)
1333         n++;
1334       if (n == graphs[0][histId]->GetN())
1335       {
1336         Printf("ERROR: Pos in graph not found: %d %d", i, histId);
1337         continue;
1338       }
1339
1340       SubtractEtaGap(hist, etaLimit, outerLimit, kFALSE, kFALSE);
1341
1342       c->cd(1);
1343       TH1* proj = hist->ProjectionX(Form("%s_proj1", hist->GetName()), hist->GetYaxis()->FindBin(-projectLimit+0.01), hist->GetYaxis()->FindBin(projectLimit-0.01));
1344       proj->SetLineColor(exColors[histId]);
1345       proj->GetXaxis()->SetRangeUser(-1.49, 1.49);
1346       proj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 1.2, proj->GetMaximum() * 1.5);
1347       proj->Draw((histId == 0) ? "" : "SAME");
1348
1349       if (drawFunc)
1350       {
1351         // integral over y
1352         TF1* func = new TF1("func", "[0]+[1]*([4]/[2]/TMath::Sqrt(TMath::TwoPi())*exp(-0.5*(x/[2])**2)+(1-[4])/[5]/TMath::Sqrt(TMath::TwoPi())*exp(-0.5*(x/[5])**2))", -0.5 * TMath::Pi(), 0.5 * TMath::Pi());
1353         func->SetParameter(0, 0);
1354         for (Int_t k=0; k<6; k++)
1355           func->SetParameter(k+1, graphs[k][histId]->GetY()[n]);
1356         // scale by bin width (to compare with projection)
1357         func->SetParameter(1, func->GetParameter(1) / hist->GetYaxis()->GetBinWidth(1));
1358         func->SetLineColor(exColors[histId]);
1359         func->SetLineWidth(1);
1360         func->DrawCopy("SAME");
1361
1362         // draw contributions
1363         Float_t scale = func->GetParameter(1);
1364         Float_t weighting = func->GetParameter(4);
1365         func->SetParameter(1, scale * weighting);
1366         func->SetParameter(4, 1);
1367         func->SetLineStyle(2);
1368         func->DrawCopy("SAME");
1369
1370         func->SetParameter(1, scale * (1.0-weighting));
1371         func->SetParameter(4, 0);
1372         func->SetLineStyle(3);
1373         func->DrawCopy("SAME");
1374       }
1375
1376       DrawLatex(0.15, 0.8 - 0.05 * histId, exColors[histId], labels[histId]);
1377
1378       c->cd(2);
1379       proj = hist->ProjectionY(Form("%s_proj2b", hist->GetName()), hist->GetXaxis()->FindBin(-projectLimit+0.01), hist->GetXaxis()->FindBin(projectLimit-0.01));
1380       proj->SetLineColor(exColors[histId]);
1381       proj->GetXaxis()->SetRangeUser(-1.49, 1.49);
1382       proj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 1.2, proj->GetMaximum() * 2);
1383       proj->Draw((histId == 0) ? "" : "SAME");
1384
1385       if (0)
1386       {
1387         proj = hist->ProjectionY(Form("%s_proj2", hist->GetName()), hist->GetXaxis()->FindBin(-0.5 * TMath::Pi()), hist->GetXaxis()->FindBin(0.5 * TMath::Pi()));
1388         proj->SetLineColor(exColors[histId]);
1389         proj->GetXaxis()->SetRangeUser(-1.2, 1.2);
1390         proj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 1.2, proj->GetMaximum() * 2);
1391         proj->SetLineStyle(2);
1392         proj->Draw("SAME");
1393       }
1394
1395       if (drawFunc)
1396       {
1397         // integral over x
1398         TF1* func = new TF1("func", "[0]+[1]*([4]/[3]/TMath::Sqrt(TMath::TwoPi())*exp(-0.5*(x/[3])**2)+(1-[4])/[6]/TMath::Sqrt(TMath::TwoPi())*exp(-0.5*(x/[6])**2))", -outerLimit, outerLimit);
1399         func->SetParameter(0, 0);
1400         for (Int_t k=0; k<6; k++)
1401           func->SetParameter(k+1, graphs[k][histId]->GetY()[n]);
1402         // scale by bin width (to compare with projection)
1403         func->SetParameter(1, func->GetParameter(1) / hist->GetXaxis()->GetBinWidth(1));
1404         func->SetLineColor(exColors[histId]);
1405         func->SetLineWidth(1);
1406         func->DrawCopy("SAME");
1407
1408         // draw contributions
1409         Float_t scale = func->GetParameter(1);
1410         Float_t weighting = func->GetParameter(4);
1411         func->SetParameter(1, scale * weighting);
1412         func->SetParameter(4, 1);
1413         func->SetLineStyle(2);
1414         func->DrawCopy("SAME");
1415
1416         func->SetParameter(1, scale * (1.0-weighting));
1417         func->SetParameter(4, 0);
1418         func->SetLineStyle(3);
1419         func->DrawCopy("SAME");
1420       }
1421
1422       DrawLatex(0.15, 0.8 - 0.05 * histId, exColors[histId], labels[histId]);
1423     }
1424   }
1425 }
1426
1427 void DrawExample(const char* histFileName, Int_t i, Int_t j, Int_t histId, TH1** projPhi, TH1** projEta)
1428 {
1429   Float_t etaLimit = 1.0;
1430   Float_t outerLimit = 1.6;
1431   Float_t projectLimit = 0.8;
1432
1433   TFile::Open(histFileName);
1434
1435   TCanvas* c = new TCanvas(Form("ex_%d_%d_%d", i, j, histId), Form("ex_%d_%d_%d", i, j, histId), 1200, 400);
1436   c->Divide(3, 1);
1437
1438   TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
1439   if (!hist)
1440     return;
1441
1442   TString label(hist->GetTitle());
1443   label.ReplaceAll(".00", ".0");
1444   TObjArray* objArray = label.Tokenize("-");
1445   TPaveText* paveText = new TPaveText(0.52, 0.72, 0.95, 0.9, "BRNDC");
1446   paveText->SetTextSize(0.04);
1447   paveText->SetFillColor(0);
1448   paveText->SetShadowColor(0);
1449   paveText->AddText(objArray->At(0)->GetName());
1450   paveText->AddText(objArray->At(1)->GetName());
1451   if (objArray->GetEntries() == 4)
1452     paveText->AddText(Form("Pb-Pb %s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()));
1453   else
1454     paveText->AddText(objArray->At(2)->GetName());
1455
1456   c->cd(1);
1457   hist->GetXaxis()->SetRangeUser(-TMath::Pi() / 2, TMath::Pi() / 2);
1458   hist->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
1459   hist->GetXaxis()->SetTitleOffset(1.5);
1460   hist->GetYaxis()->SetTitleOffset(1.5);
1461   hist->SetStats(0);
1462   hist->SetTitle("a) Correlation");
1463   hist->DrawCopy("SURF1");
1464   paveText->Draw();
1465
1466   SubtractEtaGap(hist, etaLimit, outerLimit, kFALSE, kFALSE);
1467   c->cd(2);
1468   hist->SetTitle("b) #eta-gap subtracted");
1469   hist->DrawCopy("SURF1");
1470
1471   c->cd(3);
1472   TH1* proj = hist->ProjectionX(Form("%s_proj1", hist->GetName()), hist->GetYaxis()->FindBin(-projectLimit+0.01), hist->GetYaxis()->FindBin(projectLimit-0.01));
1473   TH1* proj2 = hist->ProjectionY(Form("%s_proj2b", hist->GetName()), hist->GetXaxis()->FindBin(-projectLimit+0.01), hist->GetXaxis()->FindBin(projectLimit-0.01));
1474
1475   proj->SetStats(0);
1476   proj->SetTitle("c) Projections");
1477   proj->GetXaxis()->SetRangeUser(-1.49, 1.49);
1478   proj->GetXaxis()->SetTitleOffset(1);
1479   proj->GetYaxis()->SetRangeUser(proj->GetMinimum(), proj->GetMaximum() * 1.4);
1480   TH1* copy = proj->DrawCopy();
1481   copy->GetXaxis()->SetTitle(Form("%s / %s", proj->GetXaxis()->GetTitle(), proj2->GetXaxis()->GetTitle()));
1482   DrawLatex(0.3, 0.85, 1, Form("#Delta#phi projection in |#Delta#eta| < %.1f", projectLimit), 0.04);
1483
1484   proj2->SetLineColor(2);
1485   proj2->SetStats(0);
1486   proj2->GetXaxis()->SetTitleOffset(1);
1487   proj2->GetYaxis()->SetRangeUser(proj2->GetMinimum(), proj2->GetMaximum() * 1.6);
1488   proj2->DrawCopy("SAME");
1489   DrawLatex(0.3, 0.80, 2, Form("#Delta#eta projection in |#Delta#phi| < %.1f", projectLimit), 0.04);
1490
1491   proj->SetTitle(label);
1492   proj2->SetTitle(label);
1493
1494   *projPhi = proj;
1495   *projEta = proj2;
1496
1497   c->SaveAs(Form("ex/%s.png", c->GetName()));
1498   c->SaveAs(Form("ex/%s.eps", c->GetName()));
1499
1500   c = new TCanvas(Form("ex_%d_%d_%d_a", i, j, histId), Form("ex_%d_%d_%d_a", i, j, histId), 400, 400);
1501   hist->SetTitle("");
1502   hist->DrawCopy("SURF1");
1503   paveText->Draw();
1504   c->SaveAs(Form("ex/%s.png", c->GetName()));
1505   c->SaveAs(Form("ex/%s.eps", c->GetName()));
1506 }
1507
1508 void DrawExampleAll(const char* histFileName)
1509 {
1510   Int_t colors[] = { 1, 2, 4 };
1511
1512   Float_t exampleI[] = { 0, 1, 2 };
1513   Float_t exampleJ[] = { 1, 2, 3 };
1514
1515   TH1* projectionsPhi[9];
1516   TH1* projectionsEta[9];
1517
1518   Int_t count = 0;
1519   for (Int_t i=0; i<3; i++)
1520   {
1521     for (Int_t histId = 0; histId<3; histId++)
1522     {
1523       DrawExample(histFileName, exampleI[i], exampleJ[i], histId, &projectionsPhi[count], &projectionsEta[count]);
1524       count++;
1525     }
1526
1527     TCanvas* c = new TCanvas(Form("centralities_%d", i), Form("centralities_%d", i), 800, 400);
1528     c->Divide(2, 1);
1529
1530     TLegend* legend = new TLegend(0.62, 0.7, 0.88, 0.88);
1531     legend->SetFillColor(0);
1532
1533     for (Int_t histId = 0; histId<3; histId++)
1534     {
1535       c->cd(1);
1536       TH1* clone = projectionsPhi[count-3+histId]->DrawCopy((histId > 0) ? "SAME" : "");
1537       clone->SetLineColor(colors[histId]);
1538
1539       TString label(clone->GetTitle());
1540       TObjArray* objArray = label.Tokenize("-");
1541       clone->SetTitle(Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()));
1542
1543       legend->AddEntry(clone, (objArray->GetEntries() == 4) ? Form("%s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()) : objArray->At(2)->GetName(), "L");
1544
1545       c->cd(2);
1546       clone = projectionsEta[count-3+histId]->DrawCopy((histId > 0) ? "SAME" : "");
1547       clone->SetTitle(Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()));
1548       clone->SetLineColor(colors[histId]);
1549     }
1550
1551     c->cd(1);
1552     legend->Draw();
1553     c->SaveAs(Form("ex/%s.png", c->GetName()));
1554     c->SaveAs(Form("ex/%s.eps", c->GetName()));
1555   }
1556
1557   for (Int_t histId = 0; histId<3; histId++)
1558   {
1559     TCanvas* c = new TCanvas(Form("pt_%d", histId), Form("pt_%d", histId), 800, 400);
1560     c->Divide(2, 1);
1561
1562     TLegend* legend = new TLegend(0.15, 0.7, 0.88, 0.88);
1563     legend->SetFillColor(0);
1564
1565     for (Int_t i=2; i>=0; i--)
1566     {
1567       c->cd(1);
1568       TH1* clone = projectionsPhi[i*3+histId]->DrawCopy((i < 2) ? "SAME" : "");
1569       clone->SetLineColor(colors[i]);
1570       clone->GetYaxis()->SetRangeUser(clone->GetMinimum(), clone->GetMaximum() * 1.2);
1571
1572       TString label(clone->GetTitle());
1573       TObjArray* objArray = label.Tokenize("-");
1574       clone->SetTitle((objArray->GetEntries() == 4) ? Form("%s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()) : objArray->At(2)->GetName());
1575
1576       legend->GetListOfPrimitives()->AddFirst(new TLegendEntry(clone, Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()), "L"));
1577
1578       c->cd(2);
1579       clone = projectionsEta[i*3+histId]->DrawCopy((i < 2) ? "SAME" : "");
1580       clone->SetTitle((objArray->GetEntries() == 4) ? Form("%s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()) : objArray->At(2)->GetName());
1581       clone->SetLineColor(colors[i]);
1582 //       clone->GetYaxis()->SetRangeUser(clone->GetMinimum(), clone->GetMaximum() * 1.2);
1583     }
1584
1585     c->cd(1);
1586     legend->Draw();
1587
1588     c->SaveAs(Form("ex/%s.png", c->GetName()));
1589     c->SaveAs(Form("ex/%s.eps", c->GetName()));
1590   }
1591 }
1592
1593 void DrawFullCentralityDependence(const char* histFileName)
1594 {
1595   Int_t colors[] = { 1, 2, 4, 3, 5, 6 };
1596
1597   TH1* projectionsPhi[9];
1598   TH1* projectionsEta[9];
1599
1600   Int_t count = 0;
1601   for (Int_t histId = 0; histId<6; histId++)
1602   {
1603     DrawExample(histFileName, 0, 1, histId, &projectionsPhi[count], &projectionsEta[count]);
1604     count++;
1605   }
1606
1607   TCanvas* c = new TCanvas(Form("centralities_%d", 0), Form("centralities_%d", 0), 800, 400);
1608   c->Divide(2, 1);
1609
1610   TLegend* legend = new TLegend(0.62, 0.7, 0.88, 0.88);
1611   legend->SetFillColor(0);
1612
1613   for (Int_t histId = 0; histId<6; histId++)
1614   {
1615     c->cd(1);
1616     TH1* clone = projectionsPhi[count-6+histId]->DrawCopy((histId > 0) ? "SAME" : "");
1617     clone->SetLineColor(colors[histId]);
1618
1619     TString label(clone->GetTitle());
1620     TObjArray* objArray = label.Tokenize("-");
1621     clone->SetTitle(Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()));
1622
1623     legend->AddEntry(clone, (objArray->GetEntries() == 4) ? Form("%s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()) : objArray->At(2)->GetName(), "L");
1624
1625     c->cd(2);
1626     clone = projectionsEta[count-6+histId]->DrawCopy((histId > 0) ? "SAME" : "");
1627     clone->SetTitle(Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()));
1628     clone->SetLineColor(colors[histId]);
1629   }
1630
1631   c->cd(1);
1632   legend->Draw();
1633   c->SaveAs(Form("ex/%s.png", c->GetName()));
1634   c->SaveAs(Form("ex/%s.eps", c->GetName()));
1635 }
1636
1637 void CompareExamples(const char* histFileName1, const char* histFileName2, Int_t i, Int_t j, Int_t histId, Bool_t twoTrack = kFALSE)
1638 {
1639   Float_t etaLimit = 1.0;
1640   Float_t outerLimit = 1.79;
1641   Float_t projectLimit = 0.8;
1642
1643   Int_t exColors[] = { 1, 2, 4, 6 };
1644
1645   TCanvas* c = new TCanvas(Form("c_%d", i), Form("c_%d", i), 1000, 600);
1646   c->Divide(2, 1);
1647
1648   for (Int_t k=0; k<2; k++)
1649   {
1650     if (k == 0)
1651       TFile::Open(histFileName1);
1652     else
1653       TFile::Open(histFileName2);
1654
1655     TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
1656     if (!hist)
1657       continue;
1658 //     if (k == 1)
1659 //       hist->Scale(2);
1660 //     hist->Rebin2D(2, 2);
1661
1662 //     new TCanvas; hist->DrawCopy("COLZ"); return;
1663
1664     SubtractEtaGap(hist, etaLimit, outerLimit, kFALSE);
1665
1666   //     new TCanvas; hist->DrawCopy("COLZ");
1667
1668     // remove bins potentially affected by two-track effects
1669     if (twoTrack)
1670     {
1671       Printf("NOTE : Skipping bins at (0, 0)");
1672       for (Int_t binx = hist->GetXaxis()->FindBin(-0.25); binx <= hist->GetXaxis()->FindBin(0.25); binx++)
1673         for (Int_t biny = hist->GetYaxis()->FindBin(-0.01); biny <= hist->GetYaxis()->FindBin(0.01); biny++)
1674         {
1675           hist->SetBinContent(binx, biny,  0);
1676           hist->SetBinError(binx, biny, 0);
1677         }
1678     }
1679
1680     c->cd(1);
1681     TH1* proj = hist->ProjectionX(Form("%s_proj1%d", hist->GetName(), k), hist->GetYaxis()->FindBin(-projectLimit+0.01), hist->GetYaxis()->FindBin(projectLimit-0.01));
1682     proj->SetLineColor(exColors[k]);
1683     proj->GetXaxis()->SetRangeUser(-1.49, 1.49);
1684     proj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 1.2, proj->GetMaximum() * 1.5);
1685     proj->Draw((k == 0) ? "" : "SAME");
1686
1687     c->cd(2);
1688     proj = hist->ProjectionY(Form("%s_proj2_%d", hist->GetName(), k), hist->GetXaxis()->FindBin(-projectLimit+0.01), hist->GetXaxis()->FindBin(projectLimit-0.01));
1689     proj->SetLineColor(exColors[k]);
1690     proj->GetXaxis()->SetRangeUser(-1.49, 1.49);
1691     proj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 1.2, proj->GetMaximum() * 2);
1692     proj->Draw((k == 0) ? "" : "SAME");
1693   }
1694 }
1695
1696 void TestMomentCode()
1697 {
1698   Float_t momentLimit = 1;
1699   TH1* proj = new TH1F("hist", "", 100, -momentLimit, momentLimit);
1700   TF1* func = new TF1("func", "gaus(0)", -momentLimit, momentLimit);
1701   func->SetParameters(1, 0, 0.2);
1702   proj->FillRandom("func", 100);
1703   proj->Sumw2();
1704   proj->Draw();
1705
1706   for (Int_t n=2; n <= 4; n++)
1707   {
1708     Float_t moment = 0;
1709     Float_t sum = 0;
1710     for (Int_t bin = proj->GetXaxis()->FindBin(-momentLimit); bin <= proj->GetXaxis()->FindBin(momentLimit); bin++)
1711     {
1712       moment += proj->GetBinContent(bin) * TMath::Power(proj->GetMean() - proj->GetXaxis()->GetBinCenter(bin), n);
1713       sum += proj->GetBinContent(bin);
1714     }
1715     moment /= sum;
1716
1717     Float_t error = 0;
1718     for (Int_t bin = proj->GetXaxis()->FindBin(-momentLimit); bin <= proj->GetXaxis()->FindBin(momentLimit); bin++)
1719     {
1720       error += proj->GetBinError(bin) * proj->GetBinError(bin) *
1721         TMath::Power(TMath::Power(proj->GetMean() - proj->GetXaxis()->GetBinCenter(bin), n) / sum
1722           - moment / sum, 2);
1723     }
1724
1725     Printf("%d %f %f", n, moment, TMath::Sqrt(error));
1726   }
1727 }
1728
1729 void DivideGraphs(TGraphErrors* graph1, TGraphErrors* graph2)
1730 {
1731   graph1->Print();
1732   graph2->Print();
1733
1734   for (Int_t bin1 = 0; bin1 < graph1->GetN(); bin1++)
1735   {
1736     Float_t x = graph1->GetX()[bin1];
1737
1738     Int_t bin2 = 0;
1739     for (bin2 = 0; bin2<graph2->GetN(); bin2++)
1740       if (graph2->GetX()[bin2] >= x)
1741         break;
1742
1743     if (bin2 == graph2->GetN())
1744             bin2--;
1745
1746     if (bin2 > 0)
1747       if (TMath::Abs(graph2->GetX()[bin2-1] - x) < TMath::Abs(graph2->GetX()[bin2] - x))
1748         bin2--;
1749
1750     if (graph1->GetY()[bin1] == 0 || graph2->GetY()[bin2] == 0 || bin2 == graph2->GetN())
1751     {
1752       Printf("%d %d removed", bin1, bin2);
1753       graph1->RemovePoint(bin1--);
1754       continue;
1755     }
1756
1757     Float_t graph2Extrapolated = graph2->GetY()[bin2];
1758     if (TMath::Abs(x - graph2->GetX()[bin2]) > 0.0001)
1759     {
1760       Printf("%f %f %d %d not found", x, graph2->GetX()[bin2], bin1, bin2);
1761       graph1->RemovePoint(bin1--);
1762       continue;
1763     }
1764
1765     Float_t value = graph1->GetY()[bin1] / graph2Extrapolated;
1766     Float_t error = value * TMath::Sqrt(TMath::Power(graph1->GetEY()[bin1] / graph1->GetY()[bin1], 2) + TMath::Power(graph2->GetEY()[bin2] / graph2->GetY()[bin2], 2));
1767
1768     graph1->GetY()[bin1] = value;
1769     graph1->GetEY()[bin1] = error;
1770
1771     Printf("%d %d %f %f %f %f", bin1, bin2, x, graph2Extrapolated, value, error);
1772   }
1773 }
1774
1775 void CompareGraph(const char* fileName1, const char* fileName2, Int_t graph, Int_t centrality)
1776 {
1777   ReadGraphs(fileName1);
1778   TGraphErrors* graph1 = (TGraphErrors*) graphs[graph][centrality]->Clone("graph1");
1779
1780   ReadGraphs(fileName2);
1781   TGraphErrors* graph2 = (TGraphErrors*) graphs[graph][centrality]->Clone("graph2");
1782
1783   graph1->Sort();
1784   graph2->Sort();
1785
1786   TCanvas* c = new TCanvas(Form("%s_%s", fileName1, fileName2), Form("%s_%s", fileName1, fileName2), 800, 800);
1787   c->Divide(1, 2);
1788
1789   c->cd(1);
1790   gPad->SetGridx();
1791   gPad->SetGridy();
1792   graph1->SetMarkerStyle(24);
1793   graph1->GetXaxis()->SetRangeUser(7, 38);
1794   graph1->DrawClone("AP");
1795   graph2->SetLineColor(2);
1796   graph2->SetMarkerColor(2);
1797   graph2->SetMarkerStyle(25);
1798   graph2->DrawClone("PSAME");
1799
1800   c->cd(2);
1801   gPad->SetGridx();
1802   gPad->SetGridy();
1803   DivideGraphs(graph1, graph2);
1804   graph1->DrawClone("AP");
1805 }
1806
1807 void Compare2D(const char* fileName1, const char* fileName2, const char* histName)
1808 {
1809   TFile::Open(fileName1);
1810   TH1* hist1 = (TH1*) gFile->Get(histName);
1811
1812   TFile::Open(fileName2);
1813   TH1* hist2 = (TH1*) gFile->Get(histName);
1814
1815   hist1->Divide(hist2);
1816   hist1->GetZaxis()->SetRangeUser(0.5, 1.5);
1817   hist1->Draw("COLZ");
1818 }
1819
1820 void TestTwoGaussian()
1821 {
1822 //   TF1* func = new TF1("func", "[0]/TMath::Sqrt(TMath::TwoPi())/[2]*exp(-0.5*((x-[1])/[2])**2) + [3]/TMath::Sqrt(TMath::TwoPi())/[5]*exp(-0.5*((x-[4])/[5])**2)", -2, 2);
1823   TF1* func = new TF1("func", "[0]*exp(-0.5*((x-[1])/[2])**2) + [3]*exp(-0.5*((x-[4])/[5])**2)", -2, 2);
1824   func->SetParameters(0.25, 0, 0.5, 0.75, 0, 0.2);
1825
1826   TH1* hist = new TH1F("hist", "", 100, -2, 2);
1827   hist->FillRandom("func", 10000000);
1828
1829   hist->Draw();
1830
1831   TH1* proj = hist;
1832
1833   Float_t momentLimit = 1.99;
1834   for (Int_t n=2; n <= 4; n++)
1835   {
1836     Float_t moment = 0;
1837     Float_t sum = 0;
1838     for (Int_t bin = proj->GetXaxis()->FindBin(-momentLimit); bin <= proj->GetXaxis()->FindBin(momentLimit); bin++)
1839     {
1840       moment += proj->GetBinContent(bin) * TMath::Power(proj->GetMean() - proj->GetXaxis()->GetBinCenter(bin), n);
1841       sum += proj->GetBinContent(bin);
1842     }
1843     moment /= sum;
1844
1845     Float_t error = 0;
1846     for (Int_t bin = proj->GetXaxis()->FindBin(-momentLimit); bin <= proj->GetXaxis()->FindBin(momentLimit); bin++)
1847     {
1848       error += proj->GetBinError(bin) * proj->GetBinError(bin) *
1849         TMath::Power(TMath::Power(proj->GetMean() - proj->GetXaxis()->GetBinCenter(bin), n) / sum
1850           - moment / sum, 2);
1851     }
1852
1853     Printf("%d %f %f", n, moment, TMath::Sqrt(error));
1854   }
1855 }
1856
1857 void DrawEtaGapExample(const char* histFileName, Int_t i = 1, Int_t j = 2)
1858 {
1859   Float_t etaLimit = 1.0;
1860   Float_t outerLimit = 1.79;
1861   Float_t projectLimit = 0.8;
1862
1863   TFile::Open(histFileName);
1864
1865   TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, 0));
1866   if (!hist)
1867     return;
1868
1869   TCanvas* c = new TCanvas("c", "c", 800, 800);
1870   gPad->SetLeftMargin(0.15);
1871   TH2* clone = (TH2*) hist->DrawCopy("SURF1");
1872   clone->SetTitle("");
1873   clone->GetYaxis()->SetRangeUser(-1.79, 1.79);
1874   clone->GetXaxis()->SetTitleOffset(1.5);
1875   clone->GetYaxis()->SetTitleOffset(2);
1876   clone->SetStats(kFALSE);
1877   c->SaveAs("etagap_raw.eps");
1878
1879   SubtractEtaGap(hist, etaLimit, outerLimit, kFALSE, kTRUE);
1880
1881   c = new TCanvas("c3", "c3", 800, 800);
1882   TH1* proj = hist->ProjectionX(Form("%s_proj1%d", hist->GetName(), 0), hist->GetYaxis()->FindBin(-projectLimit+0.01), hist->GetYaxis()->FindBin(projectLimit-0.01));
1883   proj->GetXaxis()->SetRangeUser(-1.49, 1.49);
1884   proj->SetStats(0);
1885   proj->GetXaxis()->SetTitle(Form("%s / %s", proj->GetXaxis()->GetTitle(), "#Delta#eta"));
1886   proj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 1.2, proj->GetMaximum() * 1.5);
1887   proj->Draw("");
1888
1889   proj = hist->ProjectionY(Form("%s_proj2_%d", hist->GetName(), 0), hist->GetXaxis()->FindBin(-projectLimit+0.01), hist->GetXaxis()->FindBin(projectLimit-0.01));
1890   proj->GetXaxis()->SetRangeUser(-1.49, 1.49);
1891   proj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 1.2, proj->GetMaximum() * 2);
1892   proj->SetLineColor(2);
1893   proj->Draw("SAME");
1894   c->SaveAs("etagap_subtracted_proj.eps");
1895
1896   c = new TCanvas("c2", "c2", 800, 800);
1897   gPad->SetLeftMargin(0.15);
1898   hist->SetTitle("");
1899   hist->GetYaxis()->SetRangeUser(-1.79, 1.79);
1900   hist->GetXaxis()->SetTitleOffset(1.5);
1901   hist->GetYaxis()->SetTitleOffset(2);
1902   hist->SetStats(kFALSE);
1903   hist->Draw("SURF1");
1904   c->SaveAs("etagap_subtracted.eps");
1905 }