marcro updates
[u/mrichter/AliRoot.git] / PWGCF / Correlations / macros / dphicorrelations / fit.C
1 #include "TF2.h"
2 #include "TH2F.h"
3 #include "TH3.h"
4 #include "TPad.h"
5 #include "TCanvas.h"
6 #include "TMath.h"
7 #include "TGraphErrors.h"
8 #include "TGraphAsymmErrors.h"
9 #include "TFile.h"
10 #include "TLatex.h"
11 #include "TROOT.h"
12 #include "TLegend.h"
13 #include "TLegendEntry.h"
14 #include "TMultiGraph.h"
15 #include "TPaveText.h"
16 #include "TRandom.h"
17 #include "TSystem.h"
18 #include "TASImage.h"
19 #include "TObjString.h"
20 #include "TFitResult.h"
21 #include "fstream"
22 #include "TMatrixTBase.h"
23
24 void AddPoint(TGraphErrors* graph, Float_t x, Float_t y, Float_t xe, Float_t ye)
25 {
26         graph->SetPoint(graph->GetN(), x, y);
27         graph->SetPointError(graph->GetN() - 1, xe, ye);
28 }
29
30 void DrawLatex(Float_t x, Float_t y, Int_t color, const char* text, Float_t textSize = 0.06)
31 {
32   TLatex* latex = new TLatex(x, y, text);
33   latex->SetNDC();
34   latex->SetTextSize(textSize);
35   latex->SetTextColor(color);
36   latex->Draw();
37 }
38
39 // 0    1    2    3     4     5     6      7      8          9          10     11     12         13         14          15
40 // norm,dphi,deta,norm2,dphi2,deta2,chi2_1,chi2_2,moment2phi,moment2eta,phirms,etarms,moment4phi,moment4eta,kurtosisphi,kurtosiseta,
41 //            16   17   18   19    20    21    22     23     24         25         26     27     28         29         30          31         
42 // second fit:norm,dphi,deta,norm2,dphi2,deta2,chi2_1,chi2_2,moment2phi,moment2eta,phirms,etarms,moment4phi,moment4eta,kurtosisphi,kurtosiseta,
43 // 32     33              34      
44 // IAAFit,Yield(integral),IAAHist,
45 //        35      36       37    38    39      40       41    42    43       44
46 // 1D fit:normphi,norm2phi,dphi1,dphi2,normeta,norm2eta,deta1,deta2,chi2phi,chi2eta
47 // if fitting two 1D Gaussians phi rms is stored in 37 and eta rms in 41 without calling the CalculateRMS function
48 const Int_t NGraphs = 45;
49 const Int_t NHists = 6*4; // pt index
50 TGraphErrors*** graphs = 0;
51 const char* kCorrFuncTitle = "1/N_{trig} dN_{assoc}/d#Delta#etad#Delta#varphi (1/rad.)";
52 const char* kProjYieldTitlePhi = "1/N_{trig} dN_{assoc}/d#Delta#varphi (1/rad.)";
53 const char* kProjYieldTitleEta = "1/N_{trig} dN_{assoc}/d#Delta#eta";
54 const char* kProjYieldTitlePhiOrEta = "1/N_{trig} dN_{assoc}/d#Delta#varphi (1/rad.) , dN_{assoc}/d#Delta#eta";
55 TString fgFolder = "tmpresults";
56 const char* fitLabel = "fit";
57
58 void CreateGraphStructure()
59 {
60   graphs = new TGraphErrors**[NGraphs];
61   for (Int_t i=0; i<NGraphs; i++)
62   {
63     graphs[i] = new TGraphErrors*[NHists];
64     for (Int_t j=0; j<NHists; j++)
65       graphs[i][j] = new TGraphErrors;
66   }
67 }
68
69 void WriteGraphs(const char* outputFileName = "graphs.root")
70 {
71   TFile::Open(outputFileName, "RECREATE");
72   for (Int_t i=0; i<NGraphs; i++)
73     for (Int_t j=0; j<NHists; j++)
74       graphs[i][j]->Write(Form("graph_%d_%d", i, j));
75
76   gFile->Close();
77 }
78
79 void ReadGraphs(const char* fileName = "graphs.root")
80 {
81   CreateGraphStructure();
82   TFile* file = TFile::Open(fileName);
83   for (Int_t i=0; i<NGraphs; i++)
84     for (Int_t j=0; j<NHists; j++)
85       graphs[i][j] = (TGraphErrors*) file->Get(Form("graph_%d_%d", i, j));
86 }
87
88 void DivideGraphs(TGraphErrors* graph1, TGraphErrors* graph2)
89 {
90 //   graph1->Print();
91 //   graph2->Print();
92
93   for (Int_t bin1 = 0; bin1 < graph1->GetN(); bin1++)
94   {
95     Float_t x = graph1->GetX()[bin1];
96
97     Int_t bin2 = 0;
98     for (bin2 = 0; bin2<graph2->GetN(); bin2++)
99       if (graph2->GetX()[bin2] >= x)
100         break;
101
102     if (bin2 == graph2->GetN())
103             bin2--;
104
105     if (bin2 > 0)
106       if (TMath::Abs(graph2->GetX()[bin2-1] - x) < TMath::Abs(graph2->GetX()[bin2] - x))
107         bin2--;
108
109     if (graph1->GetY()[bin1] == 0 || graph2->GetY()[bin2] == 0 || bin2 == graph2->GetN())
110     {
111       Printf("%d %d removed", bin1, bin2);
112       graph1->RemovePoint(bin1--);
113       continue;
114     }
115
116     Float_t graph2Extrapolated = graph2->GetY()[bin2];
117     if (TMath::Abs(x - graph2->GetX()[bin2]) > 0.0001)
118     {
119       Printf("%f %f %d %d not found", x, graph2->GetX()[bin2], bin1, bin2);
120       graph1->RemovePoint(bin1--);
121       continue;
122     }
123
124     Float_t value = graph1->GetY()[bin1] / graph2Extrapolated;
125     Float_t error = value * TMath::Sqrt(TMath::Power(graph1->GetEY()[bin1] / graph1->GetY()[bin1], 2) + TMath::Power(graph2->GetEY()[bin2] / graph2->GetY()[bin2], 2));
126
127     graph1->GetY()[bin1] = value;
128     graph1->GetEY()[bin1] = error;
129
130 //     Printf("%d %d %f %f %f %f", bin1, bin2, x, graph2Extrapolated, value, error);
131   }
132 }
133
134 void SubtractGraphs(TGraphErrors* graph1, TGraphErrors* graph2)
135 {
136 //   graph1->Print();
137 //   graph2->Print();
138
139   for (Int_t bin1 = 0; bin1 < graph1->GetN(); bin1++)
140   {
141     Float_t x = graph1->GetX()[bin1];
142
143     Int_t bin2 = 0;
144     for (bin2 = 0; bin2<graph2->GetN(); bin2++)
145       if (graph2->GetX()[bin2] >= x)
146         break;
147
148     if (bin2 == graph2->GetN())
149             bin2--;
150
151     if (bin2 > 0)
152       if (TMath::Abs(graph2->GetX()[bin2-1] - x) < TMath::Abs(graph2->GetX()[bin2] - x))
153         bin2--;
154
155     if (graph1->GetY()[bin1] == 0 || graph2->GetY()[bin2] == 0 || bin2 == graph2->GetN())
156     {
157       Printf("%d %d removed", bin1, bin2);
158       graph1->RemovePoint(bin1--);
159       continue;
160     }
161
162     Float_t graph2Extrapolated = graph2->GetY()[bin2];
163     if (TMath::Abs(x - graph2->GetX()[bin2]) > 0.0001)
164     {
165       Printf("%f %f %d %d not found", x, graph2->GetX()[bin2], bin1, bin2);
166       graph1->RemovePoint(bin1--);
167       continue;
168     }
169
170     Float_t value = graph1->GetY()[bin1] - graph2Extrapolated;
171     Float_t error = TMath::Sqrt(TMath::Power(graph1->GetEY()[bin1], 2) + TMath::Power(graph2->GetEY()[bin2], 2));
172
173     graph1->GetY()[bin1] = value;
174     graph1->GetEY()[bin1] = error;
175
176 //     Printf("%d %d %f %f %f %f", bin1, bin2, x, graph2Extrapolated, value, error);
177   }
178 }
179
180 /* 
181 // old one with 1 Gaussian
182 Double_t DeltaPhiWidth2DFitFunction(Double_t *x, Double_t *par)
183 {
184   // params: 0: gaussian amplitude, 1: phi width, 2: eta width
185   //         3..bins+2 constants as fct of phi
186   
187   Float_t phi = x[0];
188   if (phi < 0)
189     phi = -phi;
190   if (phi > TMath::Pi())
191     phi = TMath::TwoPi() - phi;
192   Int_t phiBin = (Int_t) (phi / TMath::Pi() * 36);
193 //   phiBin = 0;
194   
195   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])));
196 }*/
197
198 void DrawALICELogo(Bool_t prel, Float_t x1, Float_t y1, Float_t x2, Float_t y2, Bool_t debug = kFALSE)
199 {
200   // correct for aspect ratio of figure plus aspect ratio of pad (coordinates are NDC!)
201 //   Printf("%d %f %d %f", gPad->GetCanvas()->GetWindowHeight(), gPad->GetHNDC(), gPad->GetCanvas()->GetWindowWidth(), gPad->GetWNDC());
202 //   x2 = x1 + (y2 - y1) * (620. / 671) * gPad->GetCanvas()->GetWindowHeight() * gPad->GetHNDC() / (gPad->GetWNDC() * gPad->GetCanvas()->GetWindowWidth());
203   x2 = x1 + (y2 - y1) * (466. / 523) * gPad->GetWh() * gPad->GetHNDC() / (gPad->GetWNDC() * gPad->GetWw());
204 //   x2 = x1 + (y2 - y1) * (620. / 671) * gPad->GetWh() / gPad->GetWw();
205
206 //   Printf("%f %f %f %f", x1, x2, y1, y2);
207   
208   TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo", x1, y1, x2, y2);
209   if (debug)
210     myPadLogo->SetFillColor(2); // color to first figure out where is the pad then comment !
211   myPadLogo->SetLeftMargin(0);
212   myPadLogo->SetTopMargin(0);
213   myPadLogo->SetRightMargin(0);
214   myPadLogo->SetBottomMargin(0);
215   myPadLogo->Draw();
216   myPadLogo->cd();
217 //   TASImage *myAliceLogo = new TASImage("~/alice_logo_transparent.png");
218   TASImage *myAliceLogo = new TASImage((prel) ? "~/alice_logo_preliminary.eps" : "~/alice_logo_performance.eps");
219   myAliceLogo->Draw();
220 }
221
222 void logotest()
223 {
224   TCanvas* c = new TCanvas("c", "c", 800, 200);
225   c->Divide(3, 1);
226   c->cd(1);
227   DrawALICELogo(1, 0.1, 0.1, 0.9, 0.9);
228   c->SaveAs("test.eps");  
229 }
230
231 const Double_t k1OverSqrtTwoPi = 1.0 / TMath::Sqrt(TMath::TwoPi());
232
233 Double_t DeltaPhiWidth2DFitFunction(Double_t *x, Double_t *par)
234 {
235   // params: 0: gaussian amplitude, 1: phi width I, 2: eta width I
236   //         3: fraction for first Gaussian, 4: phi width II, 5: eta width II
237   //         6..bins+5 constants as fct of phi
238   
239   Float_t phi = x[0];
240   if (phi < 0)
241     phi = -phi;
242   if (phi > TMath::Pi())
243     phi = TMath::TwoPi() - phi;
244   Int_t phiBin = (Int_t) (phi / TMath::Pi() * 36);
245 //   phiBin = 0;
246   
247   return par[6+phiBin] + par[0]*( 
248       par[3]/TMath::TwoPi()/par[1]/par[2] * 
249         TMath::Exp(-0.5*((x[0]/par[1])*(x[0]/par[1])+(x[1]/par[2])*(x[1]/par[2]))) 
250       + (1-par[3])/TMath::TwoPi()/par[4]/par[5] 
251         * TMath::Exp(-0.5*((x[0]/par[4])*(x[0]/par[4])+(x[1]/par[5])*(x[1]/par[5]))) 
252     );
253 }
254
255 TH2* SubtractEtaGap(TH2* hist, Float_t etaLimit, Float_t outerLimit, Bool_t scale, Bool_t drawEtaGapDist = kFALSE)
256 {
257   TString histName(hist->GetName());
258   Int_t etaBins = 0;
259
260   TH1D* etaGap = hist->ProjectionX(histName + "_1", TMath::Max(1, hist->GetYaxis()->FindBin(-outerLimit + 0.01)), hist->GetYaxis()->FindBin(-etaLimit - 0.01));
261 //   Printf("%f", etaGap->GetEntries());
262   if (etaGap->GetEntries() > 0)
263     etaBins += hist->GetYaxis()->FindBin(-etaLimit - 0.01) - TMath::Max(1, hist->GetYaxis()->FindBin(-outerLimit + 0.01)) + 1;
264
265   TH1D* tracksTmp = hist->ProjectionX(histName + "_2", hist->GetYaxis()->FindBin(etaLimit + 0.01), TMath::Min(hist->GetYaxis()->GetNbins(), hist->GetYaxis()->FindBin(outerLimit - 0.01)));
266 //   Printf("%f", tracksTmp->GetEntries());
267   if (tracksTmp->GetEntries() > 0)
268     etaBins += TMath::Min(hist->GetYaxis()->GetNbins(), hist->GetYaxis()->FindBin(outerLimit - 0.01)) - hist->GetYaxis()->FindBin(etaLimit + 0.01) + 1;
269   
270   etaGap->Add(tracksTmp);
271
272   // get per bin result
273   if (etaBins > 0)
274     etaGap->Scale(1.0 / etaBins);
275  
276   if (drawEtaGapDist)
277   {
278     TH1D* centralRegion = hist->ProjectionX(histName + "_3", hist->GetYaxis()->FindBin(-etaLimit + 0.01), hist->GetYaxis()->FindBin(etaLimit - 0.01));
279     
280 //    centralRegion->Scale(1.0 / (hist->GetYaxis()->FindBin(etaLimit - 0.01) - hist->GetYaxis()->FindBin(-etaLimit + 0.01) + 1));
281     centralRegion->Scale(hist->GetXaxis()->GetBinWidth(1));
282
283     TCanvas* c = new TCanvas("SubtractEtaGap", "SubtractEtaGap", 800, 800);
284     gPad->SetLeftMargin(0.13);
285     centralRegion->SetStats(0);
286     TString label(centralRegion->GetTitle());
287     label.ReplaceAll(".00", " GeV/c");
288     label.ReplaceAll(".0", " GeV/c");
289     centralRegion->SetTitle(label);
290     centralRegion->SetLineColor(3);
291     centralRegion->Draw();
292     centralRegion->GetYaxis()->SetTitle(kProjYieldTitlePhi);
293     centralRegion->GetYaxis()->SetTitleOffset(1.6);
294     TH1* copy = etaGap->DrawCopy("SAME");
295     copy->Scale(hist->GetXaxis()->GetBinWidth(1));
296     copy->Scale((hist->GetYaxis()->FindBin(etaLimit - 0.01) - hist->GetYaxis()->FindBin(-etaLimit + 0.01) + 1));
297     copy->SetLineColor(2);
298     TLegend* legend = new TLegend(0.41, 0.73, 0.69, 0.85);
299     legend->SetFillColor(0);
300     legend->AddEntry(centralRegion, Form("|#Delta#eta| < %.1f", etaLimit), "L");
301     legend->AddEntry(copy, Form("%.1f < |#Delta#eta| < %.1f (scaled)", etaLimit, outerLimit), "L");
302     legend->Draw();
303     
304     DrawLatex(0.705, 0.62, 1, "Pb-Pb 2.76 TeV", 0.025);
305     DrawLatex(0.705, 0.58, 1, "Stat. unc. only", 0.025);
306     
307     DrawALICELogo(kTRUE, 0.7, 0.65, 0.9, 0.85);
308     
309     c->SaveAs("note/etagap_proj.eps");
310     c->SaveAs("note/etagap_proj.png");
311   }
312   
313 //   new TCanvas; etaGap->DrawCopy();
314   
315   TH2* histTmp2D = (TH2*) hist->Clone("histTmp2D");
316   histTmp2D->Reset();
317   
318   for (Int_t xbin=1; xbin<=histTmp2D->GetNbinsX(); xbin++)
319     for (Int_t y=1; y<=histTmp2D->GetNbinsY(); y++)
320       histTmp2D->SetBinContent(xbin, y, etaGap->GetBinContent(xbin));
321     
322   if (scale)
323   {
324     // mixed event does not reproduce away-side perfectly
325     // --> extract scaling factor on the away-side from ratios of eta gap and central region
326     TH1D* centralRegion = hist->ProjectionX(histName + "_3", hist->GetYaxis()->FindBin(-etaLimit + 0.01), hist->GetYaxis()->FindBin(etaLimit - 0.01));
327     etaBins = hist->GetYaxis()->FindBin(etaLimit - 0.01) - hist->GetYaxis()->FindBin(-etaLimit + 0.01) + 1;
328     centralRegion->Scale(1.0 / etaBins);
329     
330 //     new TCanvas; centralRegion->DrawCopy(); etaGap->SetLineColor(2); etaGap->DrawCopy("SAME");
331     centralRegion->Divide(etaGap);
332 //     new TCanvas; centralRegion->Draw();
333     centralRegion->Fit("pol0", "0", "", TMath::Pi() - 1, TMath::Pi() + 1);
334     Float_t scalingFactor = centralRegion->GetFunction("pol0")->GetParameter(0);
335     Printf("  scalingFactor = %f", scalingFactor);
336     histTmp2D->Scale(scalingFactor);
337   }
338     
339 //   new TCanvas; hist->DrawCopy("SURF1");
340
341   hist->Add(histTmp2D, -1);  
342   return histTmp2D;
343 }
344
345 void SubtractEtaGap1D(TH1* projPhi, TH1* projPhiSubtractPositive, TH1* projPhiSubtractNegative, TH1* projEta, Float_t etaLimit, Float_t outerLimit, Bool_t draw = kFALSE )
346 {
347   if (draw)
348   {
349     TCanvas* c1 = new TCanvas("c1", "eta", 800, 600);
350     c1->cd();
351     projEta->GetYaxis()->SetRangeUser(0,1.1);
352     projEta->DrawCopy(); 
353   }
354   Float_t background = 0;
355   Float_t backgroundError = 0;
356   Int_t nBins = 0;
357
358   for (Int_t i=projEta->FindBin(-outerLimit + 0.01); i<=projEta->FindBin(-etaLimit - 0.01); i++)
359   {
360     nBins++;
361     background += projEta->GetBinContent(i);
362     backgroundError += TMath::Power(projEta->GetBinError(i),2);
363   }
364   for (Int_t i=projEta->FindBin(etaLimit + 0.01); i<=projEta->FindBin(outerLimit - 0.01); i++)
365   { 
366     nBins++;
367     background += projEta->GetBinContent(i);
368     backgroundError += TMath::Power(projEta->GetBinError(i),2);
369   }
370   background = background/nBins;
371   backgroundError = TMath::Sqrt(backgroundError)/nBins;
372   for (Int_t i=1; i<=projEta->GetNbinsX(); i++)
373   {
374     Float_t content = projEta->GetBinContent(i);
375     Float_t contentError = projEta->GetBinError(i);
376     projEta->SetBinContent(i,content-background);
377     projEta->SetBinError(i,TMath::Sqrt(TMath::Power(contentError,2)+TMath::Power(backgroundError,2)));
378   }
379   if (draw) projEta->DrawCopy("SAME");
380   background = 0;
381   backgroundError = 0;
382   if (draw)
383   {
384     TCanvas* c2 = new TCanvas("c2", "phi", 800, 600);
385     c2->cd();
386     projPhi->GetYaxis()->SetRangeUser(0,1.5);
387     projPhi->DrawCopy();
388     projPhiSubtractPositive->DrawCopy("SAME");
389     projPhiSubtractNegative->DrawCopy("SAME");
390   }
391   projPhiSubtractNegative->Add(projPhiSubtractPositive);
392   projPhiSubtractNegative->Scale(1.0/2.0);
393   projPhi->Add(projPhiSubtractNegative,-1);
394   if (draw) projPhi->DrawCopy("SAME");
395 }
396
397 /*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)
398 {
399   Float_t etaLimit = 1.0;
400   Float_t outerLimit = 1.6;
401   
402   SubtractEtaGap(hist, etaLimit, outerLimit, scale);
403
404 //   new TCanvas; hist->DrawCopy("SURF1");
405
406   hist->GetYaxis()->SetRangeUser(-outerLimit+0.01, outerLimit-0.01);
407   hist->GetXaxis()->SetRangeUser(-TMath::Pi() / 2 + 0.01, TMath::Pi() * 0.5 - 0.01);
408   
409   canvas->cd(canvasPos);
410   hist->SetStats(0);
411   hist->DrawCopy("SURF1");
412   
413   Float_t min = hist->GetMinimum();
414   Float_t max = hist->GetMaximum();
415   
416   // ranges are to exclude eta gap region from fit
417   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);
418   func->SetParameters(0, 1, 0.3, 0.3);
419   func->SetParLimits(1, 0, 10);
420   func->SetParLimits(2, 0.05, 1);
421   func->SetParLimits(3, 0.05, 1);
422   func->FixParameter(0, 0);
423
424 //   TF2* func = new TF2("func", DeltaPhiWidth2DFitFunction, -TMath::Pi() / 2, TMath::Pi() * 1.5, -1, 1, 4);
425 //   func->SetParameters(1, 0.3, 0.3, 0);
426 //   func->SetParLimits(0, 0, 10);
427 //   func->SetParLimits(1, 0.1, 10);
428 //   func->SetParLimits(2, 0.1, 10);
429   
430   hist->Fit(func, "0R", "");
431   hist->Fit(func, "I0R", "");
432
433 //   func->SetRange(-0.5 * TMath::Pi(), 0.5 * TMath::Pi(), -outerLimit, outerLimit);
434   
435   canvas->cd(canvasPos + 1);
436   TH2* funcHist = (TH2*) hist->Clone("funcHist");
437   funcHist->Reset();
438   funcHist->Add(func);
439   funcHist->SetMinimum(min);
440   funcHist->SetMaximum(max);
441   funcHist->Draw("SURF1");
442   
443   canvas->cd(canvasPos + 2);
444   hist->Add(func, -1);
445   hist->SetMinimum(min);
446   hist->SetMaximum(max);
447   hist->DrawCopy("SURF1");
448   
449   width[0]->SetPoint(width[0]->GetN(), x, TMath::Abs(func->GetParameter(2)));
450   width[0]->SetPointError(width[0]->GetN()-1, 0, func->GetParError(2));
451     
452   width[1]->SetPoint(width[1]->GetN(), x, TMath::Abs(func->GetParameter(3)));
453   width[1]->SetPointError(width[1]->GetN()-1, 0, func->GetParError(3));
454
455   Float_t chi2 = 0;
456   Int_t ndf = 0;
457   for (Int_t i=hist->GetXaxis()->FindBin(-0.8); i<=hist->GetXaxis()->FindBin(0.8); i++)
458     for (Int_t j=hist->GetYaxis()->FindBin(-0.8); j<=hist->GetYaxis()->FindBin(0.8); j++)
459     {
460       if (hist->GetBinError(i, j) > 0)
461       {
462         chi2 += TMath::Power(hist->GetBinContent(i, j) / hist->GetBinError(i, j), 2);
463         ndf++;
464       }
465     }
466   ndf -= func->GetNumberFreeParameters();
467   
468   printf("#chi^{2}/ndf = %.1f/%d = %.1f  ", func->GetChisquare(), func->GetNDF(), func->GetChisquare() / func->GetNDF());
469   Printf("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf);
470
471   DrawLatex(0.5, yPosChi2, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", func->GetChisquare(), func->GetNDF(), func->GetChisquare() / func->GetNDF()));
472   DrawLatex(0.5, yPosChi2 - 0.05, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf));
473
474   chi2_1->SetPoint(chi2_1->GetN(), x, func->GetChisquare() / func->GetNDF());
475   chi2_2->SetPoint(chi2_2->GetN(), x, chi2 / ndf);
476 }
477 */
478
479 void CalculateMomentsKurtosis(Float_t momentCalcLimit, TH1* proj, Int_t graphIDStart, Int_t histId, Float_t x, Float_t xE)
480 {
481   //   momentCalcLimit = 0.6;
482   for (Int_t n=2; n <= 4; n+=2) // skipping moment 3
483   {
484     Float_t moment = 0;
485     Float_t sum = 0;
486     for (Int_t bin = proj->GetXaxis()->FindBin(-momentCalcLimit); bin <= proj->GetXaxis()->FindBin(momentCalcLimit); bin++)
487     {
488       moment += proj->GetBinContent(bin) * TMath::Power(proj->GetMean() - proj->GetXaxis()->GetBinCenter(bin), n);
489       sum += proj->GetBinContent(bin);
490     }
491     Printf("%f %f", moment, sum);
492     moment /= sum;
493     
494     Float_t error = 0;
495     for (Int_t bin = proj->GetXaxis()->FindBin(-momentCalcLimit); bin <= proj->GetXaxis()->FindBin(momentCalcLimit); bin++)
496     {
497       error += proj->GetBinError(bin) * proj->GetBinError(bin) * 
498         TMath::Power(TMath::Power(proj->GetMean() - proj->GetXaxis()->GetBinCenter(bin), n) / sum 
499           - moment / sum, 2);
500     }
501     AddPoint(graphs[graphIDStart+(n-2)*2][histId], x, moment, xE, TMath::Sqrt(error));
502     Printf("%d %d %f +- %f <-> %f +- %f", n, graphIDStart+(n-2)*2, moment, TMath::Sqrt(error), proj->GetRMS() * proj->GetRMS(), 2 * proj->GetRMSError() / proj->GetRMS() * proj->GetRMS() * proj->GetRMS());
503   }
504   
505   proj->GetXaxis()->SetRangeUser(-momentCalcLimit, momentCalcLimit);
506   AddPoint(graphs[graphIDStart+6][histId], x, proj->GetKurtosis(1), xE, proj->GetKurtosis(11));
507 }
508
509 void AddRMS(TGraphErrors* graph, Float_t x, Float_t xE, TF1* func, TMatrixDSym& cov, Int_t sigma1, Int_t sigma2, Int_t weight)
510 {
511   Float_t a = TMath::Abs(func->GetParameter(sigma1));
512   Float_t b = TMath::Abs(func->GetParameter(sigma2));
513   Float_t w = func->GetParameter(weight);
514
515   Float_t rms = a * w + b * (1.0 - w);
516
517   /*
518   func->Print(); Printf("%f %f %f", a, b, w); cov.Print();
519   Printf("%f %f %f %f %f %f", TMath::Power(w * func->GetParError(sigma1), 2),
520     TMath::Power((1.0 - w) * func->GetParError(sigma2), 2),
521     TMath::Power((a - b) * func->GetParError(weight), 2),
522     2 * (a - b) * w * cov(weight, sigma1),
523     2 * (a - b) * (1.0 - w) * cov(weight, sigma2),
524     2 * w * (1.0 - w) * cov(sigma1, sigma2));
525   */
526   
527   Float_t sigma_rms = 
528     TMath::Power(w * func->GetParError(sigma1), 2) + 
529     TMath::Power((1.0 - w) * func->GetParError(sigma2), 2) + 
530     TMath::Power((a - b) * func->GetParError(weight), 2) + 
531     2 * (a - b) * w * cov(weight, sigma1) +
532     2 * (a - b) * (1.0 - w) * cov(weight, sigma2) +
533     2 * w * (1.0 - w) * cov(sigma1, sigma2);
534     
535    if (sigma_rms < 0)
536    {
537      Printf("WARNING: error is negative!");
538      
539      sigma_rms = TMath::Power(w * func->GetParError(sigma1), 2) + 
540      TMath::Power((1.0 - w) * func->GetParError(sigma2), 2) + 
541      TMath::Power((a - b) * func->GetParError(weight), 2);
542    }
543     
544   sigma_rms = TMath::Sqrt(sigma_rms);
545   
546   Printf("RMS: %f +- %f", rms, sigma_rms);
547   AddPoint(graph, x, rms, xE, sigma_rms);
548 }
549
550 void AddKurtosis(TGraphErrors* graph, Float_t x, Float_t xE, TF1* func, TMatrixDSym& cov, Int_t sigma1, Int_t sigma2, Int_t weight)
551 {
552   Float_t a = TMath::Abs(func->GetParameter(sigma1));
553   Float_t b = TMath::Abs(func->GetParameter(sigma2));
554   Float_t w = func->GetParameter(weight);
555
556   Float_t rms = a * w + b * (1.0 - w);
557   Float_t kurtosis = 3 * (w * TMath::Power(a, 4) + (1.0 - w) * TMath::Power(b, 4)) / rms / rms - 3;
558   
559   Float_t der_a = 12*a*a*a*w / TMath::Power(b*b+(a*a-b*b)*w, 2) - (12*a*w*(b*b*b*b+(a*a*a*a-b*b*b*b)*w))/TMath::Power(b*b+(a*a-b*b)*w, 3);
560   Float_t der_b = (3*(4*b*b*b-4*b*b*b*w))/TMath::Power((b*b+(a*a-b*b)*w), 2) - (6*(2*b-2*b*w)*(b*b*b*b+(a*a*a*a-b*b*b*b)* w))/TMath::Power(b*b+(a*a-b*b)*w, 3);
561   Float_t der_w = (3*(a*a*a*a-b*b*b*b))/TMath::Power(b*b+(a*a-b*b)*w, 2) - (6*(a*a-b*b)*(b*b*b*b+(a*a*a*a-b*b*b*b)*w))/TMath::Power(b*b+(a*a-b*b)*w,3);
562   
563   Float_t sigma_kurtosis = 
564     TMath::Power(der_a * func->GetParError(sigma1), 2) + 
565     TMath::Power(der_b * func->GetParError(sigma2), 2) + 
566     TMath::Power(der_w * func->GetParError(weight), 2) + 
567     2 * der_w * der_a * cov(weight, sigma1) +
568     2 * der_w * der_b * cov(weight, sigma2) +
569     2 * der_a * der_b * cov(sigma1, sigma2);
570     
571 //   if (sigma_rms < 0)
572 //   {
573 //     Printf("WARNING: error is negative!");
574 //     
575 //     sigma_rms = TMath::Power(w * func->GetParError(sigma1), 2) + 
576 //     TMath::Power((1.0 - w) * func->GetParError(sigma2), 2) + 
577 //     TMath::Power((a - b) * func->GetParError(weight), 2);
578 //   }
579     
580   sigma_kurtosis = TMath::Sqrt(sigma_kurtosis);
581   
582   Printf("Kurtosis: %f +- %f", kurtosis, sigma_kurtosis);
583   AddPoint(graph, x, kurtosis, xE, sigma_kurtosis);
584 }
585
586 Float_t kEtaLimit = 1.0;
587 Float_t kOuterLimit = 1.59;
588
589
590 Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int_t graphID, Float_t x, Float_t xE, Float_t yPosChi2, Bool_t quick, Int_t histId, Int_t limits, Bool_t /*twoTrack*/)
591 {
592   Bool_t success = kTRUE;
593   Float_t etaLimit = kEtaLimit;
594   Float_t outerLimit = kOuterLimit;
595
596   // for pt,T < 3
597   if (graphID <= 10)
598     etaLimit = 1.4;
599   if (graphID <= 15)
600     etaLimit = 1.2;
601   
602   // for pT,a > 4
603   if (graphID == 18 || graphID == 23 || graphID == 24)
604   {
605     etaLimit = 0.5;
606     outerLimit = 0.99;
607   }
608   Float_t sigmaFitLimit = 0.1 - limits * 0.05;
609   Float_t etaFitUpperLimit = 0.8;
610   Float_t initSigma = 0.6;
611   if (histId == 2) // pp
612   {
613     etaFitUpperLimit = 0.6;
614     initSigma = 0.4;
615   }
616   if ((graphID == 10 && histId == 4) || (graphID == 5 && histId == 0))
617     etaFitUpperLimit += 0.1;
618   
619   Printf("Limits: etaLimit = %f outerLimit = %f sigmaFitLimit = %f etaFitUpperLimit = %f initSigma = %f", etaLimit, outerLimit, sigmaFitLimit, etaFitUpperLimit, initSigma);
620
621   // set errors large for bins around (0,0)
622   if (1 && graphID <= 12)
623   {
624     Float_t exclusionRegion = 0.19;
625     if (graphID == 0)
626       exclusionRegion = 0.29;
627     if (graphID == 10 || graphID == 11 || graphID == 12)
628       exclusionRegion = 0.075;
629     
630     Printf("NOTE : Skipping bins at (0, 0)");
631     for (Int_t binx = hist->GetXaxis()->FindBin(-exclusionRegion); binx <= hist->GetXaxis()->FindBin(exclusionRegion); binx++)
632       for (Int_t biny = hist->GetYaxis()->FindBin(-exclusionRegion); biny <= hist->GetYaxis()->FindBin(exclusionRegion); biny++)
633       {
634 //      hist->SetBinContent(binx, biny,  0);
635         hist->SetBinError(binx, biny,  1e5);
636       }
637   }
638
639   hist->GetYaxis()->SetRangeUser(-outerLimit+0.01, outerLimit-0.01);
640   hist->GetXaxis()->SetRangeUser(-TMath::Pi() / 2 + 0.01, TMath::Pi() * 0.5 - 0.01);
641   
642   Float_t mean = hist->Integral(hist->GetYaxis()->FindBin(-TMath::Pi() / 2), hist->GetYaxis()->FindBin(TMath::Pi() / 2), hist->GetYaxis()->FindBin(etaLimit), hist->GetYaxis()->FindBin(outerLimit)) / (hist->GetYaxis()->FindBin(TMath::Pi() / 2) - hist->GetYaxis()->FindBin(-TMath::Pi() / 2)) / (hist->GetYaxis()->FindBin(outerLimit) - hist->GetYaxis()->FindBin(etaLimit) + 1);
643 //   Printf("%f", mean);
644
645   // sums
646   TGraphErrors* sumSummary = new TGraphErrors;
647   TGraphErrors* phiWidthSummary = new TGraphErrors;
648   TGraphErrors* etaWidthSummary = new TGraphErrors;
649   
650   canvas->cd(canvasPos++);
651   hist->SetStats(0);
652   hist->DrawCopy("SURF1");
653   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), hist->Integral());
654   
655   Float_t min = hist->GetMinimum();
656   Float_t max = hist->GetMaximum();
657
658   Int_t bins = hist->GetNbinsX() / 2 / 2;
659     
660   TF2* func = new TF2("func", DeltaPhiWidth2DFitFunction, -TMath::Pi() / 2, TMath::Pi() * 0.5, -outerLimit, outerLimit, bins+6);
661   func->SetParameters(hist->GetBinContent(hist->GetXaxis()->FindBin(0.0), hist->GetYaxis()->FindBin(0.0)) - mean, 0.3, 0.3, 0.25, initSigma, initSigma);
662   for (Int_t i=6; i<bins+6; i++)
663     func->SetParameter(i, mean);
664  
665  
666   if (1)
667   {
668     // STEP 1: fit only flow using one delta eta side
669     for (Int_t i=0; i<6; i++)
670       func->FixParameter(i, func->GetParameter(i));
671     hist->GetYaxis()->SetRangeUser(etaLimit+0.01, outerLimit-0.01);
672     Int_t fitResult = hist->Fit(func, "0", "");
673     Printf("Fit result: %d; Chi2/ndf: %f/%d", fitResult, func->GetChisquare(), func->GetNDF());
674     if (fitResult != 0)
675       success = kFALSE;
676
677     // STEP2 : fit only Gaussian in central region
678     for (Int_t i=0; i<6; i++)
679       func->ReleaseParameter(i);
680     //   func->SetParameters(1, 0.3, 0.3, 0.25, initSigma, initSigma);
681     func->SetParLimits(0, 0, 10);
682     func->SetParLimits(1, sigmaFitLimit, 0.6);
683     func->SetParLimits(2, sigmaFitLimit, etaFitUpperLimit);
684     func->SetParLimits(3, 0.05, 0.95);
685     func->SetParLimits(4, sigmaFitLimit, 0.601);
686     func->SetParLimits(5, sigmaFitLimit, etaFitUpperLimit);
687     for (Int_t i=6; i<bins+6; i++)
688       func->FixParameter(i, func->GetParameter(i));
689     hist->GetYaxis()->SetRangeUser(-etaLimit+0.01, etaLimit-0.01);
690     fitResult = hist->Fit(func, "0", "");
691     Printf("Fit result: %d; Chi2/ndf: %f/%d", fitResult, func->GetChisquare(), func->GetNDF());
692     if (fitResult != 0)
693       success = kFALSE;
694
695     // STEP3: fit everything, with limits
696     for (Int_t i=6; i<bins+6; i++)
697     {
698       func->ReleaseParameter(i);
699       func->SetParLimits(i, func->GetParameter(i) * 0.8, func->GetParameter(i) * 1.2);
700     }
701     hist->GetYaxis()->SetRangeUser(-outerLimit+0.01, outerLimit-0.01);
702
703     // STEP4: fit everything, without limits
704     for (Int_t i=6; i<bins+6; i++)
705       func->SetParLimits(i, 0, 0);
706     fitResult = hist->Fit(func, "0", "");
707     Printf("Fit result: %d; Chi2/ndf: %f/%d", fitResult, func->GetChisquare(), func->GetNDF());
708     if (fitResult != 0)
709       success = kFALSE;
710   }
711   TFitResultPtr fitResultPtr = hist->Fit(func, "S0", "");
712   TMatrixDSym cov = fitResultPtr->GetCovarianceMatrix();
713   Printf("Fit result: %d; Chi2/ndf: %f/%d", (Int_t) fitResultPtr, func->GetChisquare(), func->GetNDF());
714   if (fitResultPtr != 0)
715     success = kFALSE;
716
717   Printf("Trying 1 Gaussian...");
718   TF2* func_clone = new TF2("func_clone", DeltaPhiWidth2DFitFunction, -TMath::Pi() / 2, TMath::Pi() * 0.5, -outerLimit, outerLimit, bins+6);
719   for (Int_t i=0; i<bins+6; i++)
720   {
721     func_clone->SetParameter(i, func->GetParameter(i));
722     Double_t parmin, parmax;
723     func->GetParLimits(i, parmin, parmax);
724     func_clone->SetParLimits(i, parmin, parmax);
725   }
726   func_clone->SetParLimits(3, 1, 1);
727   func_clone->FixParameter(3, 1);
728   func_clone->FixParameter(4, sigmaFitLimit);
729   func_clone->FixParameter(5, sigmaFitLimit);
730   TFitResultPtr fitResultPtr2 = hist->Fit(func_clone, "S0R", "");
731   Printf("Fit result: %d", (Int_t) fitResultPtr2);
732   
733   // if both parameters are within 1%, refit with 1 Gaussian only
734   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)
735   {
736     Printf("Parameters within 1%%. Using result from 1 Gaussian...");
737     
738     func = func_clone;
739     cov = fitResultPtr2->GetCovarianceMatrix();
740     
741     if (fitResultPtr2 != 0)
742       success = kFALSE;
743   }
744   else if (func_clone->GetChisquare() * 0.99 < func->GetChisquare())
745   {
746     Printf("1 Gaussian fit has a at maximum 1%% (%f %f) larger chi2. Using results from 1 Gaussian...", func_clone->GetChisquare(), func->GetChisquare());
747     
748     func = func_clone;
749     cov = fitResultPtr2->GetCovarianceMatrix();
750     
751     if (fitResultPtr2 != 0)
752       success = kFALSE;
753   } 
754   
755   Int_t first = 1;
756   Int_t second = 4;
757   if (func->GetParameter(1) < func->GetParameter(4))
758   {
759     first = 4;
760     second = 1;
761   }
762   //dphi
763   AddPoint(graphs[1+16][graphID], x, TMath::Abs(func->GetParameter(first)), xE, func->GetParError(first));
764   AddPoint(graphs[4+16][graphID], x, TMath::Abs(func->GetParameter(second)), xE, func->GetParError(second));
765
766   //deta
767   AddPoint(graphs[2+16][graphID], x, TMath::Abs(func->GetParameter(first+1)), xE, func->GetParError(first+1));
768   AddPoint(graphs[5+16][graphID], x, TMath::Abs(func->GetParameter(second+1)), xE, func->GetParError(second+1));
769   
770   // norm
771   AddPoint(graphs[0+16][graphID], x, func->GetParameter(0), xE, func->GetParError(0));
772   if (first < second)
773     AddPoint(graphs[3+16][graphID], x, func->GetParameter(3), xE, func->GetParError(3));
774   else
775     AddPoint(graphs[3+16][graphID], x, 1.0 - func->GetParameter(3), xE, func->GetParError(3));    
776   
777   //rms
778   AddRMS(graphs[10+16][graphID], x, xE, func, cov, 1, 4, 3);
779   AddRMS(graphs[11+16][graphID], x, xE, func, cov, 1+1, 4+1, 3);
780   AddKurtosis(graphs[14+16][graphID], x, xE, func, cov, 1, 4, 3);
781   AddKurtosis(graphs[15+16][graphID], x, xE, func, cov, 1+1, 4+1, 3);
782   
783   //   hist->Fit(func, "0RI", "");
784
785   //   width[0]->SetPoint(width[0]->GetN(), x, TMath::Abs(func->GetParameter(1)));
786   //   width[0]->SetPointError(width[0]->GetN()-1, 0, func->GetParError(1));
787   //     
788   //   width[1]->SetPoint(width[1]->GetN(), x, TMath::Abs(func->GetParameter(2)));
789   //   width[1]->SetPointError(width[1]->GetN()-1, 0, func->GetParError(2));
790
791   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), func->GetParameter(1), 0, func->GetParError(1));
792   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), func->GetParameter(4), 0, func->GetParError(4));
793   
794   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), func->GetParameter(2), 0, func->GetParError(2));
795   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), func->GetParameter(5), 0, func->GetParError(5));
796
797   canvas->cd(canvasPos++);
798   TH2* funcHist = (TH2*) hist->Clone("funcHist");
799   funcHist->Reset();
800   funcHist->Add(func);
801   funcHist->SetMinimum(min);
802   funcHist->SetMaximum(max);
803   funcHist->Draw("SURF1");
804   
805   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), funcHist->Integral());
806   
807   canvas->cd(canvasPos++);
808   TH2* residuals = (TH2*) hist->Clone("residuals");
809   residuals->Add(func, -1);
810   residuals->SetMinimum(-(max - min) / 2);
811   residuals->SetMaximum((max - min) / 2);
812   residuals->Draw("SURF1");
813
814   Double_t chi2 = 0;
815   Int_t ndf = 0;
816   for (Int_t i=hist->GetXaxis()->FindBin(-0.8); i<=hist->GetXaxis()->FindBin(0.8); i++)
817     for (Int_t j=hist->GetYaxis()->FindBin(-etaLimit); j<=hist->GetYaxis()->FindBin(etaLimit); j++)
818     {
819       if (residuals->GetBinError(i, j) > 0)
820       {
821         chi2 += TMath::Power(residuals->GetBinContent(i, j) / residuals->GetBinError(i, j), 2);
822         ndf++;
823       }
824     }
825   ndf -= func->GetNumberFreeParameters();
826   
827   if (func->GetNDF() > 0)
828   {
829     printf("#chi^{2}/ndf = %.1f/%d = %.1f  ", func->GetChisquare(), func->GetNDF(), func->GetChisquare() / func->GetNDF());
830     DrawLatex(0.2, yPosChi2, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", func->GetChisquare(), func->GetNDF(), func->GetChisquare() / func->GetNDF()));
831     AddPoint(graphs[6+16][graphID], x, func->GetChisquare() / func->GetNDF(), xE, 0);
832   }
833   if (ndf)
834   {
835     Printf("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf);
836     DrawLatex(0.2, yPosChi2 - 0.05, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf));
837     AddPoint(graphs[7+16][graphID], x, chi2 / ndf, xE, 0);
838   }
839   
840   // draw gaussian only
841   TF2* funcClone = new TF2("funcClone", DeltaPhiWidth2DFitFunction, -TMath::Pi() / 2, TMath::Pi() * 0.5, -outerLimit, outerLimit, bins+6);
842   for (Int_t i=0; i<6; i++)
843     funcClone->SetParameter(i, func->GetParameter(i));
844   for (Int_t i=6; i<bins+6; i++)
845     funcClone->SetParameter(i, 0);
846 //   funcClone->Print();
847   canvas->cd(canvasPos++);
848   funcHist = (TH2*) hist->Clone("funcHistb");
849   funcHist->Reset();
850   funcHist->Add(funcClone);
851   funcHist->SetMinimum(0);
852   funcHist->SetMaximum(max - min);
853   funcHist->Draw("SURF1");
854   
855   // eta gap subtraction
856   canvas->cd(canvasPos++);
857   func->SetParameter(0, 0);
858   TH2* subtractFlow = (TH2*) hist->Clone("subtractFlow");
859   subtractFlow->Add(func, -1);
860   subtractFlow->SetMinimum(0);
861   subtractFlow->SetMaximum(max - min);
862   subtractFlow->DrawCopy("SURF1");
863   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), subtractFlow->Integral());
864 /*
865   Float_t etaLimitSubtract = 0;
866   Float_t outerLimitSubtract = 0;
867   if (graphID >= 15)
868   {
869     etaLimitSubtract = 0.5;
870     outerLimitSubtract = 0.99;
871   }
872   else 
873   {
874     etaLimitSubtract = etaLimit;
875     outerLimitSubtract = outerLimit;
876   }
877 */
878
879   Float_t momentFitLimit = 0.8 - 1e-4;
880   if (graphID >= 20)
881     momentFitLimit = 0.4 - 1e-4;
882
883   TH1* projx3 = hist->ProjectionX(Form("%s_projx3", hist->GetName()), hist->GetYaxis()->FindBin(-etaLimit+0.01), hist->GetYaxis()->FindBin(etaLimit-0.01));
884   Float_t nBins = hist->GetYaxis()->FindBin(etaLimit-0.01) - hist->GetYaxis()->FindBin(-etaLimit+0.01)+1;
885   projx3->Scale(1.0/nBins);
886   
887   TH1* projx3SubtractNegative = hist->ProjectionX(Form("%s_projx3SubtractPositive", hist->GetName()), hist->GetYaxis()->FindBin(-outerLimit+0.01), hist->GetYaxis()->FindBin(-etaLimit-0.01));
888   nBins = hist->GetYaxis()->FindBin(-etaLimit-0.01) - hist->GetYaxis()->FindBin(-outerLimit+0.01)+1;
889   projx3SubtractNegative->Scale(1.0/nBins);
890
891   TH1* projx3SubtractPositive = hist->ProjectionX(Form("%s_projx3SubtractNegative ", hist->GetName()), hist->GetYaxis()->FindBin(etaLimit+0.01), hist->GetYaxis()->FindBin(outerLimit-0.01));
892   nBins = hist->GetYaxis()->FindBin(outerLimit-0.01) - hist->GetYaxis()->FindBin(etaLimit+0.01)+1;
893   projx3SubtractPositive->Scale(1.0/nBins);
894
895   TH1* projy3 = hist->ProjectionY(Form("%s_projy3", hist->GetName()), hist->GetXaxis()->FindBin(-momentFitLimit+0.01), hist->GetXaxis()->FindBin(momentFitLimit-0.01));
896   nBins = hist->GetXaxis()->FindBin(momentFitLimit-0.01) - hist->GetXaxis()->FindBin(-momentFitLimit+0.01)+1;
897   projy3->Scale(1.0/nBins);
898
899   SubtractEtaGap1D(projx3, projx3SubtractPositive, projx3SubtractNegative, projy3, etaLimit, outerLimit);
900   projy3->GetXaxis()->SetRangeUser(-momentFitLimit+0.01,momentFitLimit-0.01);
901   projx3->GetXaxis()->SetRangeUser(-momentFitLimit+0.01, momentFitLimit-0.01);
902   SubtractEtaGap(hist, etaLimit, outerLimit, kFALSE);
903 //  return 0;
904   canvas->cd(canvasPos++);
905   hist->SetMinimum(-(max - min) / 2);
906   hist->SetMaximum((max - min) / 2);
907   hist->DrawCopy("SURF1");
908   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), hist->Integral());
909
910   if (!quick)
911   {
912     canvas->cd(canvasPos++);
913     TH2* difference = (TH2*) hist->Clone("difference");
914     difference->Add(subtractFlow, -1);
915     difference->SetMinimum(-(max - min) / 2);
916     difference->SetMaximum((max - min) / 2);
917     difference->DrawCopy("SURF1");
918
919     canvas->cd(canvasPos++);
920     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);
921     func2->SetParameters(0, 1, 0.3, 0.3);
922     func2->SetParLimits(1, 0, 10);
923     func2->SetParLimits(2, sigmaFitLimit, 1);
924     func2->SetParLimits(3, sigmaFitLimit, 1);
925     func2->FixParameter(0, 0);
926     
927     hist->Fit(func2, "0R", "");
928     //   hist->Fit(func2, "I0R", "");
929     
930     AddPoint(phiWidthSummary, phiWidthSummary->GetN(), func2->GetParameter(2), 0, func2->GetParError(2));
931     AddPoint(etaWidthSummary, etaWidthSummary->GetN(), func2->GetParameter(3), 0, func2->GetParError(3));
932   
933     TH2* funcHist2 = (TH2*) hist->Clone("funcHist2");
934     funcHist2->Reset();
935     funcHist2->Add(func2);
936     funcHist2->SetMinimum(-(max - min) / 2);
937     funcHist2->SetMaximum((max - min) / 2);
938     funcHist2->Draw("SURF1");
939     sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), funcHist2->Integral());
940     
941     if (func2->GetNDF() > 0)
942     {
943       Printf("#chi^{2}/ndf = %.1f/%d = %.1f", func2->GetChisquare(), func2->GetNDF(), func2->GetChisquare() / func2->GetNDF());
944       DrawLatex(0.2, yPosChi2, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", func2->GetChisquare(), func2->GetNDF(), func2->GetChisquare() / func2->GetNDF()));
945     }
946
947     canvas->cd(canvasPos++);
948     TH2* residuals2 = (TH2*) hist->Clone("residuals");
949     residuals2->Add(funcHist2, -1);
950     residuals2->SetMinimum(-(max - min) / 2);
951     residuals2->SetMaximum((max - min) / 2);
952     residuals2->Draw("SURF1");
953   }
954
955
956   TH1* projx2 = hist->ProjectionX(Form("%s_projx2", hist->GetName()), hist->GetYaxis()->FindBin(-etaLimit+0.01), hist->GetYaxis()->FindBin(etaLimit-0.01));
957   projx2->GetXaxis()->SetRangeUser(-momentFitLimit, momentFitLimit);
958   nBins = hist->GetYaxis()->FindBin(etaLimit-0.01) - hist->GetYaxis()->FindBin(-etaLimit+0.01)+1;
959   projx2->Scale(1.0/nBins);
960   
961   TH1* projy2 = hist->ProjectionY(Form("%s_projy2", hist->GetName()), hist->GetXaxis()->FindBin(-momentFitLimit+0.01), hist->GetXaxis()->FindBin(momentFitLimit-0.01));
962   projy2->GetXaxis()->SetRangeUser(-momentFitLimit, momentFitLimit);
963   nBins = hist->GetXaxis()->FindBin(momentFitLimit-0.01) - hist->GetXaxis()->FindBin(-momentFitLimit+0.01)+1;
964   projy2->Scale(1.0/nBins);
965   
966 //   CalculateMomentsKurtosis(momentFitLimit, projx2, 8, graphID, x, xE);
967 //   CalculateMomentsKurtosis(momentFitLimit, projy2, 9, graphID, x, xE);
968
969 //     return success;
970   
971   TH1* projx1 = subtractFlow->ProjectionX(Form("%s_projx1", hist->GetName()), hist->GetYaxis()->FindBin(-0.79), hist->GetYaxis()->FindBin(0.79));
972   projx1->GetXaxis()->SetRangeUser(-momentFitLimit, momentFitLimit);
973
974   TH1* projy1 = subtractFlow->ProjectionY(Form("%s_projy1", hist->GetName()), hist->GetXaxis()->FindBin(-0.79), hist->GetXaxis()->FindBin(0.79));
975   projy1->GetXaxis()->SetRangeUser(-momentFitLimit, momentFitLimit);
976
977 //   CalculateMomentsKurtosis(momentFitLimit, projx1, 8+16, graphID, x, xE);
978 //   CalculateMomentsKurtosis(momentFitLimit, projy1, 9+16, graphID, x, xE);
979
980 //   TF1* twoGauss = new TF1("twoGauss", "gaus(0)+gaus(3)", -2, 2);
981 //   twoGauss->SetParameters(1, 0, 0.3, 1, 0, 0.6);
982 //   twoGauss->FixParameter(1, 0);
983 //   twoGauss->FixParameter(4, 0);
984 //   twoGauss->SetLineColor(4);
985 //   projx1->Fit("twoGauss", "I+", "SAME");
986   
987   canvas->cd(canvasPos++);
988   projx1->Draw();
989   projx1->Fit("gaus", "I");
990   projx1->GetFunction("gaus")->SetLineColor(1);
991   projx1->GetYaxis()->SetRangeUser(0, projx1->GetMaximum() * 1.05);
992
993   projx2->SetLineColor(2);
994   projx2->Draw("SAME");
995   projx2->Fit("gaus", "I+0", "SAME");
996   projx2->GetFunction("gaus")->SetLineColor(2);
997
998   canvas->cd(canvasPos--);
999   projy1->Draw();
1000   projy1->Fit("gaus", "I");
1001   projy1->GetFunction("gaus")->SetLineColor(1);
1002   projy1->GetYaxis()->SetRangeUser(0, projy1->GetMaximum() * 1.05);
1003   
1004   projy2->SetLineColor(2);
1005   projy2->Draw("SAME");
1006   projy2->Fit("gaus", "I+0", "SAME");
1007   projy2->GetFunction("gaus")->SetLineColor(2);
1008   
1009   // 1d fit (eta gap subtraction)
1010   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), projx2->GetFunction("gaus")->GetParameter(2), 0, projx2->GetFunction("gaus")->GetParError(2));
1011   
1012   // 1d fit (lots of params)
1013   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), projx1->GetFunction("gaus")->GetParameter(2), 0, projx1->GetFunction("gaus")->GetParError(2));
1014
1015   // 1d fit (eta gap subtraction)
1016   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), projy2->GetFunction("gaus")->GetParameter(2), 0, projy2->GetFunction("gaus")->GetParError(2));
1017   
1018   // 1d fit (lots of params)
1019   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), projy1->GetFunction("gaus")->GetParameter(2), 0, projy1->GetFunction("gaus")->GetParError(2));
1020
1021   Float_t etaFitLimit = outerLimit;
1022
1023   Bool_t oneGaussian = kTRUE;
1024
1025   if (!oneGaussian) cout << endl << "Fitting two 1D Gaussians" << endl;
1026   else cout << endl << "Fitting one 1D Gaussians" << endl;
1027
1028   TF1* func4phi = new TF1("func4phi", "[0]+[1]*([4]/TMath::Sqrt(TMath::TwoPi())/[2]*exp(-0.5*((x/[2])**2))+(1-[4])/TMath::Sqrt(TMath::TwoPi())/[3]*exp(-0.5*((x/[3])**2)))", -0.5 * TMath::Pi()+0.01, 0.5 * TMath::Pi()-0.01);
1029   TF1* func4eta = new TF1("func4eta", "[0]+[1]*([4]/TMath::Sqrt(TMath::TwoPi())/[2]*exp(-0.5*((x/[2])**2))+(1-[4])/TMath::Sqrt(TMath::TwoPi())/[3]*exp(-0.5*((x/[3])**2)))", -etaFitLimit+0.01, etaFitLimit-0.01);
1030
1031   func4phi->SetParLimits(1, 0, 10);
1032   func4phi->SetParLimits(3, sigmaFitLimit, 0.7);
1033   func4phi->FixParameter(0, 0);
1034   if (!oneGaussian)
1035   {
1036     func4phi->SetParLimits(4, 0, 1);
1037     func4phi->SetParLimits(2, sigmaFitLimit, 0.7);
1038   }
1039   else
1040   {
1041     func4phi->FixParameter(4, 0);
1042     func4phi->FixParameter(2, 0);
1043   }
1044
1045   func4eta->SetParLimits(1, 0, 10);
1046   func4eta->SetParLimits(3, sigmaFitLimit, etaFitUpperLimit);
1047   func4eta->FixParameter(0, 0);
1048   if (!oneGaussian)
1049   {
1050     func4eta->SetParLimits(4, 0, 1);
1051     func4eta->SetParLimits(2, sigmaFitLimit, etaFitUpperLimit);
1052   }
1053   else
1054   {
1055     func4eta->FixParameter(4, 0);
1056     func4eta->FixParameter(2, 0);
1057   }
1058   canvas->cd(canvasPos++);
1059   func4phi->SetLineColor(4);
1060   fitResultPtr = projx3->Fit(func4phi, "SIR", "SAME");
1061   TMatrixDSym cov2 = fitResultPtr->GetCovarianceMatrix();
1062
1063   projx3->Draw("SAME");
1064   projx3->SetLineColor(4);
1065
1066   if (func4phi->GetNDF() > 0)
1067   {
1068     printf("#chi^{2}/ndf = %.1f/%d = %.1f  ", func4phi->GetChisquare(), func4phi->GetNDF(), func4phi->GetChisquare() / func4phi->GetNDF());
1069     AddPoint(graphs[43][graphID], x, func4phi->GetChisquare() / func4phi->GetNDF(), xE, 0);
1070   }
1071   chi2 = 0;
1072   ndf = 0;
1073
1074   for (Int_t i=projx3->FindBin(-momentFitLimit); i<=projx3->FindBin(momentFitLimit);i++)
1075   {
1076 //    Float_t temp = chi2;
1077     if (projx3->GetBinError(i) > 0)
1078     {
1079       chi2 += TMath::Power((projx3->GetBinContent(i)-func4phi->Eval(projx3->GetBinCenter(i)))/projx3->GetBinError(i),2);
1080       //    cerr << i << "\t" << projx3->GetBinCenter(i) << "\t" << chi2-temp << endl;
1081       ndf++;
1082     }
1083   }
1084   if (!oneGaussian) ndf = ndf - 4;
1085   else ndf = ndf - 2;
1086   printf("Calculated #chi^{2}/ndf = %.1f/%d = %.1f  ", chi2, ndf, chi2/ndf);
1087
1088   first = 2;
1089   second = 3;
1090   if (func4phi->GetParameter(2) < func4phi->GetParameter(3))
1091   {
1092     first = 3;
1093     second = 2;
1094   }
1095   
1096   if (!oneGaussian)
1097   {
1098     Float_t vector[3];
1099     for (Int_t i=2; i<5;i++)
1100       vector[i-2] = cov2[i][2]*func4phi->GetParameter(4) + cov2[i][3]*(1-func4phi->GetParameter(4))+cov2[i][4]*(TMath::Abs(func4phi->GetParameter(2))-TMath::Abs(func4phi->GetParameter(3)));
1101
1102     Float_t sigma = TMath::Sqrt(vector[0]*func4phi->GetParameter(4) + vector[1]*(1-func4phi->GetParameter(4))+ vector[2]*(TMath::Abs(func4phi->GetParameter(2))-TMath::Abs(func4phi->GetParameter(3))));
1103     Float_t rms = TMath::Abs(func4phi->GetParameter(2))*func4phi->GetParameter(4)+TMath::Abs(func4phi->GetParameter(3))*(1-func4phi->GetParameter(4));
1104     AddPoint(graphs[37][graphID], x, rms, xE, sigma);
1105   }
1106   else AddPoint(graphs[37][graphID], x, TMath::Abs(func4phi->GetParameter(first)), xE, func4phi->GetParError(first));
1107   AddPoint(graphs[38][graphID], x, 0, xE, 0);
1108   AddPoint(graphs[36][graphID], x, 1, xE, 0);
1109
1110   AddPoint(graphs[35][graphID], x, func4phi->GetParameter(1), xE, func4phi->GetParError(1));
1111
1112   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), TMath::Abs(func4phi->GetParameter(first)), 0, func4phi->GetParError(first));
1113   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), TMath::Abs(func4phi->GetParameter(second)), 0, func4phi->GetParError(second));
1114
1115   Float_t scale = func4phi->GetParameter(1);
1116   Float_t weighting = func4phi->GetParameter(4);
1117   func4phi->SetParameter(1, scale * weighting);
1118   func4phi->SetParameter(4, 1);
1119   func4phi->SetLineStyle(2);
1120   func4phi->DrawCopy("SAME");
1121
1122   func4phi->SetParameter(1, scale * (1.0-weighting));
1123   func4phi->SetParameter(4, 0);
1124   func4phi->SetLineStyle(3);
1125   func4phi->DrawCopy("SAME");
1126
1127   fitResultPtr->Print("V");
1128   if (fitResultPtr != 0)
1129     success = kFALSE;
1130
1131   canvas->cd(canvasPos++);
1132   func4eta->SetLineColor(4);
1133   fitResultPtr = projy3->Fit(func4eta, "SIR", "SAME");
1134   cov2 = fitResultPtr->GetCovarianceMatrix();
1135
1136   projy3->Draw("SAME");
1137   projy3->SetLineColor(4);
1138
1139   chi2 = 0;
1140   ndf = 0;
1141   if (func4eta->GetNDF() > 0)
1142   {
1143     printf("#chi^{2}/ndf = %.1f/%d = %.1f  ", func4eta->GetChisquare(), func4eta->GetNDF(), func4eta->GetChisquare() / func4eta->GetNDF());
1144     AddPoint(graphs[44][graphID], x, func4eta->GetChisquare() / func4eta->GetNDF(), xE, 0);
1145   }
1146   for (Int_t i=projy3->FindBin(-momentFitLimit); i<=projy3->FindBin(momentFitLimit);i++)
1147   {
1148 //    Float_t temp = chi2;
1149     if (projy3->GetBinError(i) <= 0)
1150       continue;
1151     
1152     chi2 += TMath::Power((projy3->GetBinContent(i)-func4eta->Eval(projy3->GetBinCenter(i)))/projy3->GetBinError(i),2);
1153 //    cerr << i << "\t" << projy3->GetBinCenter(i) << "\t" << chi2-temp << endl;
1154     ndf++;
1155   }
1156   if (!oneGaussian) ndf = ndf - 4;
1157   else ndf = ndf - 2;
1158   printf("Calculated #chi^{2}/ndf = %.1f/%d = %.1f  ", chi2, ndf, chi2/ndf);
1159
1160   first = 2;
1161   second = 3;
1162   if (func4eta->GetParameter(2) < func4eta->GetParameter(3))
1163   {
1164     first = 3;
1165     second = 2;
1166   }
1167
1168   AddPoint(graphs[39][graphID], x, func4eta->GetParameter(1), xE, func4eta->GetParError(1));
1169   if (!oneGaussian)
1170   {
1171     Float_t vector[3];
1172     for (Int_t i=2; i<5;i++)
1173       vector[i-2] = cov2[i][2]*func4eta->GetParameter(4) + cov2[i][3]*(1-func4eta->GetParameter(4))+cov2[i][4]*(TMath::Abs(func4eta->GetParameter(2))-TMath::Abs(func4eta->GetParameter(3)));
1174
1175     Float_t sigma = TMath::Sqrt(vector[0]*func4eta->GetParameter(4) + vector[1]*(1-func4eta->GetParameter(4))+ vector[2]*(TMath::Abs(func4eta->GetParameter(2))-TMath::Abs(func4eta->GetParameter(3))));
1176     Float_t rms = TMath::Abs(func4eta->GetParameter(2))*func4eta->GetParameter(4)+TMath::Abs(func4eta->GetParameter(3))*(1-func4eta->GetParameter(4));
1177     AddPoint(graphs[41][graphID], x, rms, xE, sigma);
1178   }
1179   else AddPoint(graphs[41][graphID], x, TMath::Abs(func4eta->GetParameter(first)), xE, func4eta->GetParError(first));
1180   AddPoint(graphs[42][graphID], x, 0, xE, 0);
1181   AddPoint(graphs[40][graphID], x, 1, xE, 0);
1182
1183   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), TMath::Abs(func4eta->GetParameter(first)), 0, func4eta->GetParError(first));
1184   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), TMath::Abs(func4eta->GetParameter(second)), 0, func4eta->GetParError(second));
1185
1186   scale = func4eta->GetParameter(1);
1187   weighting = func4eta->GetParameter(4);
1188   func4eta->SetParameter(1, scale * weighting);
1189   func4eta->SetParameter(4, 1);
1190   func4eta->SetLineStyle(2);
1191   func4eta->DrawCopy("SAME");
1192
1193   func4eta->SetParameter(1, scale * (1.0-weighting));
1194   func4eta->SetParameter(4, 0);
1195   func4eta->SetLineStyle(3);
1196   func4eta->DrawCopy("SAME");
1197   
1198   fitResultPtr->Print("V");
1199   if (fitResultPtr != 0)
1200     success = kFALSE;
1201
1202 //  return 0;
1203
1204   // 2d fit with two gaussians
1205   canvas->cd(canvasPos++);
1206 //   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);
1207 //   func3->SetParameters(0, 1, 0.3, 0.3, 1, 0.6, 0.6);
1208 //   func3->SetParLimits(4, 0, 10);
1209 //   Float_t etaFitLimit = 0.5;
1210   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);
1211
1212   func3->SetParameters(0, hist->GetBinContent(hist->GetXaxis()->FindBin(0.0), hist->GetYaxis()->FindBin(0.0)), 0.3, 0.3, 0.25, initSigma, initSigma);
1213
1214   func3->FixParameter(0, 0);
1215
1216   func3->SetParLimits(4, 0.05, 0.95);
1217   func3->SetParLimits(1, 0, 10);
1218   func3->SetParLimits(2, sigmaFitLimit, 0.6);
1219   func3->SetParLimits(3, sigmaFitLimit, etaFitUpperLimit);
1220   func3->SetParLimits(5, sigmaFitLimit, 0.601);
1221   func3->SetParLimits(6, sigmaFitLimit, etaFitUpperLimit);
1222
1223 /*
1224   sigmaFitLimit = 0.09;
1225   etaFitUpperLimit = 0.12;
1226   func3->SetParLimits(4, 0.1, 0.9);
1227   func3->SetParLimits(1, 0, 10);
1228   func3->SetParLimits(2, sigmaFitLimit, 0.12);
1229   func3->SetParLimits(3, sigmaFitLimit, etaFitUpperLimit);
1230   func3->SetParLimits(5, sigmaFitLimit, 0.12);
1231   func3->SetParLimits(6, sigmaFitLimit, etaFitUpperLimit);
1232   func3->FixParameter(0, 0);
1233 */
1234 //   for (Int_t i=0; i<6; i++)
1235 //     func3->SetParameter(i+1, func->GetParameter(i));
1236
1237   fitResultPtr = hist->Fit(func3, "0RS", "");
1238   Printf("Fit result: %d", (Int_t) fitResultPtr);
1239 //   return 0;
1240   if (fitResultPtr != 0)
1241     success = kFALSE;
1242   TMatrixDSym cov3 = fitResultPtr->GetCovarianceMatrix();
1243   cov3.Print();
1244   
1245   Printf("Testing 1 Gaussian...");
1246
1247   TF2* func3_clone = (TF2*) func3->Clone("func3_clone");
1248   
1249   func3_clone->SetParLimits(4, 1, 1);
1250   func3_clone->FixParameter(4, 1);
1251   func3_clone->FixParameter(5, sigmaFitLimit);
1252   func3_clone->FixParameter(6, sigmaFitLimit);
1253   
1254   fitResultPtr2 = hist->Fit(func3_clone, "0RS", "");
1255   Printf("Fit result: %d", (Int_t) fitResultPtr2);
1256   
1257   // if both parameters were within 1%, use 1 Gaussian parameters
1258   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)
1259   {
1260     Printf("Parameters within 1%%. Using results from 1 Gaussian...");
1261     
1262     if (fitResultPtr2 != 0)
1263       success = kFALSE;
1264     func3 = func3_clone;
1265     cov3 = fitResultPtr2->GetCovarianceMatrix();
1266 //     cov3.Print();
1267   }
1268   else if (func3_clone->GetChisquare() * 0.99 < func3->GetChisquare())
1269   {
1270     Printf("1 Gaussian fit has a at maximum 1%% (%f/%d %f/%d) larger chi2. Using results from 1 Gaussian...", func3_clone->GetChisquare(), func3_clone->GetNDF(), func3->GetChisquare(), func3_clone->GetNDF());
1271     
1272     if (fitResultPtr2 != 0)
1273       success = kFALSE;
1274     func3 = func3_clone;
1275     cov3 = fitResultPtr2->GetCovarianceMatrix();
1276   }
1277   
1278 //   hist->Fit(func3, "I0R", "");
1279
1280   first = 2;
1281   second = 5;
1282   if (func3->GetParameter(2) < func3->GetParameter(5))
1283   {
1284     first = 5;
1285     second = 2;
1286   }
1287   //dphi
1288   AddPoint(graphs[1][graphID], x, TMath::Abs(func3->GetParameter(first)), xE, func3->GetParError(first));
1289   AddPoint(graphs[4][graphID], x, TMath::Abs(func3->GetParameter(second)), xE, func3->GetParError(second));
1290
1291   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), TMath::Abs(func3->GetParameter(first)), 0, func3->GetParError(first));
1292   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), TMath::Abs(func3->GetParameter(second)), 0, func3->GetParError(second));
1293     
1294   //deta
1295   AddPoint(graphs[2][graphID], x, TMath::Abs(func3->GetParameter(first+1)), xE, func3->GetParError(first+1));
1296   AddPoint(graphs[5][graphID], x, TMath::Abs(func3->GetParameter(second+1)), xE, func3->GetParError(second+1));
1297   
1298   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), TMath::Abs(func3->GetParameter(first+1)), 0, func3->GetParError(first+1));
1299   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), TMath::Abs(func3->GetParameter(second+1)), 0, func3->GetParError(second+1));
1300   
1301   // norm
1302   AddPoint(graphs[0][graphID], x, func3->GetParameter(1), xE, func3->GetParError(1));
1303   if (first < second)
1304     AddPoint(graphs[3][graphID], x, func3->GetParameter(4), xE, func3->GetParError(4));
1305   else
1306     AddPoint(graphs[3][graphID], x, 1.0 - func3->GetParameter(4), xE, func3->GetParError(4));
1307
1308   AddRMS(graphs[10][graphID], x, xE, func3, cov3, 2, 5, 4);
1309   AddRMS(graphs[11][graphID], x, xE, func3, cov3, 2+1, 5+1, 4);
1310   
1311   TH2* funcHist3 = (TH2*) hist->Clone("funcHist3");
1312   funcHist3->Reset();
1313   funcHist3->Add(func3);
1314   funcHist3->SetMinimum(-(max - min) / 2);
1315   funcHist3->SetMaximum((max - min) / 2);
1316   funcHist3->Draw("SURF1");
1317   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), funcHist3->Integral());
1318   
1319   canvas->cd(canvasPos++);
1320   TH2* residuals3 = (TH2*) hist->Clone("residuals");
1321   residuals3->Add(funcHist3, -1);
1322   residuals3->SetMinimum(-(max - min) / 2);
1323   residuals3->SetMaximum((max - min) / 2);
1324   residuals3->Draw("SURF1");
1325
1326   chi2 = 0;
1327   ndf = 0;
1328   for (Int_t i=hist->GetXaxis()->FindBin(-0.8); i<=hist->GetXaxis()->FindBin(0.8); i++)
1329     for (Int_t j=hist->GetYaxis()->FindBin(-etaLimit); j<=hist->GetYaxis()->FindBin(etaLimit); j++)
1330     {
1331       if (residuals3->GetBinError(i, j) > 0)
1332       {
1333         chi2 += TMath::Power(residuals3->GetBinContent(i, j) / residuals3->GetBinError(i, j), 2);
1334         ndf++;
1335       }
1336     }
1337   ndf -= func3->GetNumberFreeParameters();
1338   
1339   if (func3->GetNDF() > 0)
1340   {
1341     printf("#chi^{2}/ndf = %.1f/%d = %.1f  ", func3->GetChisquare(), func3->GetNDF(), func3->GetChisquare() / func3->GetNDF());
1342     DrawLatex(0.2, yPosChi2, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", func3->GetChisquare(), func3->GetNDF(), func3->GetChisquare() / func3->GetNDF()));
1343     AddPoint(graphs[6][graphID], x, func3->GetChisquare() / func3->GetNDF(), xE, 0);
1344   }
1345   if (ndf)
1346   {
1347     Printf("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf);
1348     DrawLatex(0.2, yPosChi2 - 0.05, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf));
1349     AddPoint(graphs[7][graphID], x, chi2 / ndf, xE, 0);
1350   }
1351
1352   canvas->cd(canvasPos++);
1353   phiWidthSummary->SetMarkerStyle(20);
1354   phiWidthSummary->Draw("AP");
1355   gPad->SetGridy();
1356     
1357   etaWidthSummary->SetMarkerStyle(21);
1358   etaWidthSummary->SetLineColor(2);
1359   etaWidthSummary->SetMarkerColor(2);
1360   etaWidthSummary->Draw("PSAME");  
1361   
1362   phiWidthSummary->GetYaxis()->SetRangeUser(0, 0.9);
1363   
1364   canvas->cd(canvasPos++);
1365   sumSummary->Draw("*A");
1366   gPad->SetGridy();
1367   
1368   Printf("Finished with %d", success);
1369
1370   return success;
1371 }
1372
1373 void AnalyzeDeltaPhiEtaGap2D(const char* fileName, const char* outputFileName = "graphs.root")
1374 {
1375   gROOT->SetBatch(kTRUE);
1376   if (!gROOT->IsBatch())
1377   {
1378     Printf("Not in batch mode. Exiting!");
1379     return;
1380   }
1381   
1382   CreateGraphStructure();
1383
1384   TFile::Open(fileName);
1385   TList checkList;
1386   
1387   Int_t maxLeadingPt = 4;
1388   Int_t maxAssocPt = 6;
1389
1390
1391   for (Int_t i=0; i<maxLeadingPt; i++)
1392   {
1393     for (Int_t j=1; j<maxAssocPt; j++)
1394     {
1395       // only process when first is filled
1396       TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, 0));
1397       if (!hist1)
1398       {
1399         cout << "Hist1 does not exist." << endl;
1400         continue;
1401       }
1402 //      if (hist1->GetEntries() < 1e4)
1403       if (hist1->GetEntries() < 10)
1404       {
1405         Printf("%d %d Only %f entries. Skipping...", i, j, hist1->GetEntries());
1406         continue;
1407       }
1408       
1409       for (Int_t histId = 0; histId < 8; histId++)
1410       {
1411 //      if (i != 1 || j != 2)
1412 //        continue;
1413
1414         // skip 4.0 < p_{T,trig} < 8.0 - 6.00 < p_{T,assoc} < 8.00 above 40% and pp
1415         if (i == 2 && j == 5 && (histId == 1 || histId == 2 || histId == 5))
1416           continue;
1417         
1418         hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
1419         if (!hist1)
1420         {
1421           cout << "Hist1 does not exist." << endl;
1422           continue;
1423         }
1424         
1425 //      if (hist1->GetEntries() < 1e4)
1426         if (hist1->GetEntries() < 10)
1427         {
1428           Printf("%d %d %d Only %f entries. Skipping...", i, j, histId, hist1->GetEntries());
1429           continue;
1430         }
1431         
1432         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);
1433         canvas->Divide(5, 3);
1434         
1435         for (Int_t k=1; k<=3; k++)
1436         {
1437           canvas->cd(3 * j + k);
1438           gPad->SetLeftMargin(0.15);
1439           gPad->SetBottomMargin(0.2);
1440           gPad->SetTopMargin(0.01);
1441           gPad->SetRightMargin(0.01);
1442         }
1443         
1444         Int_t graphID = i * (maxAssocPt - 1) + j - 1;
1445
1446         Printf("\n\n>>> %d %d %d %d", i, j, histId, graphID);
1447         
1448         if (histId == 0)
1449           for (Int_t k=0; k<NGraphs; k++)
1450             graphs[k][graphID]->SetTitle(hist1->GetTitle());
1451         
1452         Float_t centralityAxisMapping[] = { 5, 70, 100, 30, 15, 50 };
1453         Float_t centralityAxisMappingE[] = { 5, 10, 0, 10, 5, 10 };
1454
1455         Bool_t success = FitDeltaPhi2DOneFunction((TH2*) hist1, canvas, 1, graphID, centralityAxisMapping[histId], centralityAxisMappingE[histId], 0.9, kTRUE, histId, (i > 0) ? 1 : 0, (j <= 2 && histId != 2));
1456         if (!success)
1457           checkList.Add(new TObjString(Form("AnalyzeDeltaPhiEtaGap2DExample(\"%s\", %d, %d, %d)", fileName, i, j, histId)));
1458         
1459 //      canvas->SaveAs(Form("%s.png", canvas->GetName()));
1460 //      delete canvas;
1461 //      break;
1462 // return;
1463       }
1464       
1465 //       break;
1466     }
1467     
1468 //     break;
1469   }
1470   
1471   WriteGraphs(outputFileName);
1472
1473   for (Int_t i=0; i<checkList.GetEntries(); i++)
1474     Printf("%s", checkList.At(i)->GetName());
1475 }
1476
1477 void AnalyzeDeltaPhiEtaGap2DExample(const char* fileName, Int_t i, Int_t j, Int_t histId, Bool_t drawDetails = kFALSE)
1478 {
1479   CreateGraphStructure();
1480
1481   TFile::Open(fileName);
1482   
1483   TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
1484   if (!hist1)
1485     return;
1486   
1487   Printf("Entries: %f %s", hist1->GetEntries(), hist1->GetTitle());
1488
1489   TString label(hist1->GetTitle());
1490   if (drawDetails)
1491     hist1->SetTitle("");
1492   hist1->GetYaxis()->SetRangeUser(-1.79, 1.79);
1493   hist1->GetXaxis()->SetTitleOffset(1.5);
1494   hist1->GetYaxis()->SetTitleOffset(2);
1495   hist1->GetZaxis()->SetTitle(kCorrFuncTitle);
1496   hist1->GetZaxis()->SetTitleOffset(1.8);
1497   hist1->SetStats(kFALSE);
1498
1499 //  if (hist1->GetEntries() < 1e4)
1500   if (hist1->GetEntries() < 10)
1501     return;
1502   
1503   // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
1504   if (drawDetails)
1505     hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
1506
1507   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);
1508   canvas->Divide(5, 3);
1509   
1510   Int_t graphID = i * (6 - 1) + j - 1;
1511   FitDeltaPhi2DOneFunction((TH2*) hist1, canvas, 1, graphID, 0, 0, 0.9, kFALSE, histId, (i > 0) ? 1 : 0, (j <= 2 && histId != 2));
1512   
1513   if (!drawDetails)
1514     return;
1515   
1516   TVirtualPad* pad = canvas->cd(6);
1517   TCanvas* c = new TCanvas("c", "c", 1500, 1000);
1518   c->Divide(3, 2);
1519   c->cd(1);
1520   pad->SetPad(0, 0, 1, 1);
1521   pad->SetLeftMargin(0.15);
1522   pad->Draw();
1523   pad->cd();
1524
1525   label.ReplaceAll(".00", " GeV/c");
1526   label.ReplaceAll(".0", " GeV/c");
1527   TObjArray* objArray = label.Tokenize("-");
1528   TPaveText* paveText = new TPaveText(0.52, 0.72, 0.95, 0.95, "BRNDC");
1529   paveText->SetTextSize(0.035);
1530   paveText->SetFillColor(0);
1531   paveText->SetShadowColor(0);
1532   paveText->AddText(objArray->At(0)->GetName());
1533   paveText->AddText(objArray->At(1)->GetName());
1534   if (objArray->GetEntries() == 4)
1535     paveText->AddText(Form("Pb-Pb 2.76 TeV %s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()));
1536   else
1537     paveText->AddText(Form("%s 2.76 TeV", objArray->At(2)->GetName()));
1538   paveText->AddText("|#eta| < 0.9");
1539   paveText->Draw();
1540   
1541   c->cd(1);
1542 //   DrawALICELogo(kTRUE, 0.2, 0.7, 0.4, 0.9);
1543   
1544 //   return;
1545     
1546   //   c->SaveAs("fit_subtracted.eps");
1547
1548   pad = canvas->cd(12);
1549   c->cd(2);
1550   pad->SetPad(0, 0, 1, 1);
1551   pad->SetLeftMargin(0.15);
1552   pad->Draw();
1553 //   c->SaveAs("fit_fit1.eps");
1554   
1555   pad = canvas->cd(13);
1556   c->cd(3);
1557   gPad->SetLeftMargin(0.15);
1558   pad->SetPad(0, 0, 1, 1);
1559   pad->SetLeftMargin(0.15);
1560   pad->Draw();
1561   c->SaveAs("fit_residual1.eps");
1562   
1563   TCanvas* c2 = new TCanvas("c3b", "c3b", 800, 800);
1564   gPad->SetLeftMargin(0.15);
1565   gPad->SetGridy();
1566   hist1 = (TH2*) pad->GetListOfPrimitives()->First();
1567   for (i=0; i<10; i++)
1568   {
1569     TH1* proj = hist1->ProjectionX(Form("p_%d", i), 11+i*2, 12+i*2);
1570     proj->Scale(0.5);
1571     
1572     proj->Add(new TF1("func", "1", -10, 10), -0.5 + 0.1 * i);
1573     proj->SetStats(0);
1574     proj->GetXaxis()->SetTitleOffset(1.0);
1575     proj->GetYaxis()->SetTitleOffset(1.4);
1576     proj->GetYaxis()->SetNdivisions(512);
1577     proj->SetYTitle("Residuals (projected / shifted)");
1578     proj->GetYaxis()->SetRangeUser(-0.59, 0.49);
1579     proj->Draw((i == 0) ? "" : "SAME");    
1580   }
1581   c2->SaveAs("fit_residual1_proj.eps");
1582
1583   pad = canvas->cd(4);
1584   c->cd(4);
1585   pad->SetPad(0, 0, 1, 1);
1586   pad->SetLeftMargin(0.15);
1587   pad->Draw();
1588 //   c->SaveAs("fit_fit2.eps");
1589
1590   pad = canvas->cd(3);
1591   c->cd(5);
1592   pad->SetPad(0, 0, 1, 1);
1593   pad->SetLeftMargin(0.15);
1594   pad->Draw();
1595 //   c->SaveAs("fit_residual2.eps");
1596
1597   c2 = new TCanvas("c5b", "c5b", 800, 800);
1598   gPad->SetLeftMargin(0.15);
1599   gPad->SetGridy();
1600   hist1 = (TH2*) pad->GetListOfPrimitives()->First();
1601   for (i=0; i<10; i++)
1602   {
1603     TH1* proj = hist1->ProjectionX(Form("p2_%d", i), 11+i*2, 12+i*2);
1604     proj->Scale(0.5);
1605     
1606     proj->Add(new TF1("func", "1", -10, 10), -0.5 + 0.1 * i);
1607     proj->SetStats(0);
1608     proj->GetXaxis()->SetTitleOffset(1.0);
1609     proj->GetYaxis()->SetTitleOffset(1.4);
1610     proj->GetYaxis()->SetNdivisions(512);
1611     proj->SetYTitle("Residuals (projected / shifted)");
1612     proj->GetYaxis()->SetRangeUser(-0.59, 0.49);
1613     proj->Draw((i == 0) ? "" : "SAME");    
1614   }
1615   c2->SaveAs("fit_residual2_proj.eps");
1616
1617   pad = canvas->cd(7);
1618   c->cd(6);
1619   pad->SetPad(0, 0, 1, 1);
1620   pad->SetLeftMargin(0.15);
1621   pad->Draw();
1622   
1623   c->SaveAs("fit_example.eps");
1624   c->SaveAs("fit_example.png");
1625 }
1626
1627 void SqrtAll(Int_t nHists, TGraphErrors** graph)
1628 {
1629   for (Int_t histId = 0; histId < nHists; histId++)
1630   {
1631     for (Int_t i=0; i<graph[histId]->GetN(); i++)
1632     {
1633       if (graph[histId]->GetY()[i] > 0)
1634       {
1635         graph[histId]->GetEY()[i] *= 0.5 / TMath::Sqrt(graph[histId]->GetY()[i]);
1636         graph[histId]->GetY()[i] = TMath::Sqrt(graph[histId]->GetY()[i]);
1637       }
1638       else
1639       {
1640         Printf("WARNING negative value %d %d, removing from histogram", histId, i);
1641         graph[histId]->RemovePoint(i);
1642         i--;
1643       }
1644     }
1645   }
1646 }
1647
1648 void SqrtAll2(Int_t nHists, TGraphErrors** graph, TGraphErrors** graph2)
1649 {
1650   for (Int_t histId = 0; histId < nHists; histId++)
1651   {
1652     if (!graph[histId])
1653       continue;
1654  
1655     for (Int_t i=0; i<graph[histId]->GetN(); i++)
1656     {
1657       if (graph[histId]->GetY()[i] > 0 && graph2[histId]->GetY()[i] > 0)
1658       {
1659         graph[histId]->GetEY()[i] *= 0.5 / TMath::Sqrt(graph[histId]->GetY()[i]);
1660         graph[histId]->GetY()[i] = TMath::Sqrt(graph[histId]->GetY()[i]);
1661
1662         graph2[histId]->GetEY()[i] *= 0.5 / TMath::Sqrt(graph2[histId]->GetY()[i]);
1663         graph2[histId]->GetY()[i] = TMath::Sqrt(graph2[histId]->GetY()[i]);
1664       }
1665       else
1666       {
1667         Printf("WARNING negative value %d %d, removing from histogram", histId, i);
1668         graph[histId]->RemovePoint(i);
1669         graph2[histId]->RemovePoint(i);
1670         i--;
1671       }
1672     }
1673   }
1674 }
1675
1676 Int_t marker[6] = { 24, 25, 26, 30, 27, 28 };
1677 Int_t marker2[6] = { 20, 21, 22, 29, 33, 34 };
1678 const char* labels[6] = { "0-10%", "60-70%", "pp", "20-40%", "10-20%", "40-60%" };
1679
1680 void CalculateRMS(Int_t nHists, TGraphErrors** graph, TGraphErrors** graph2, TGraphErrors** weight)
1681 {
1682   // rms
1683   for (Int_t histId = 0; histId < nHists; histId++)
1684   {
1685     if (!graph[histId])
1686       continue;
1687     
1688     if (graph[histId]->GetN() != graph2[histId]->GetN() || graph[histId]->GetN() != weight[histId]->GetN())
1689     {
1690       Printf("E-CalculateRMS: Different number of points for hist %d: %d %d %d", histId, graph[histId]->GetN(), graph2[histId]->GetN(), weight[histId]->GetN());
1691       continue;
1692     }
1693     
1694     for (Int_t i=0; i<graph[histId]->GetN(); i++)
1695     {
1696       Float_t rms = graph[histId]->GetY()[i] * weight[histId]->GetY()[i] + graph2[histId]->GetY()[i] * (1.0 - weight[histId]->GetY()[i]);
1697       graph[histId]->GetY()[i] = rms;
1698       // error**2 = weight**2 * error_sigma**2 + (weight-1)**2 * error_sigma2**2 + (sigma1 + sigma2)**2 * error_weight**2
1699       // TODO this neglects some correlations
1700       graph[histId]->GetEY()[i] = TMath::Sqrt(
1701         TMath::Power(weight[histId]->GetY()[i] * graph[histId]->GetEY()[i], 2) + 
1702         TMath::Power((weight[histId]->GetY()[i] - 1) * graph2[histId]->GetEY()[i], 2) + 
1703         TMath::Power((graph[histId]->GetY()[i] - graph2[histId]->GetY()[i]) * weight[histId]->GetEY()[i], 2)
1704       );
1705     }
1706   }
1707 }
1708
1709 Bool_t SkipGraph(Int_t i)
1710 {
1711   return (i == 14);
1712 //  return (i == 1 || i == 7 || i == 9 || i == 13 || i == 14);
1713 //   return (i == 7 || i == 9 || i == 13 || i == 14); //For the STAR comparison
1714 }
1715
1716 Int_t skipGraphList[] =  { -1, -1, -1 };
1717 Bool_t SkipGraphForThisPlot(Int_t graphId)
1718 {
1719   for (Int_t i=0; i<3; i++)
1720     if (skipGraphList[i] == graphId)
1721       return kTRUE;
1722   return kFALSE;
1723 }
1724
1725 TGraphAsymmErrors* FixGraph(TGraphErrors* graph, Float_t shift)
1726 {
1727   graph->Sort();
1728 //   graph->Print();
1729   
1730   if (graph->GetN() > 4 && graph->GetX()[4] == 75)
1731   {
1732     graph->GetX()[4] = 65;
1733     graph->GetEX()[4] = 5;
1734   }
1735   
1736   TGraphAsymmErrors* res = new TGraphAsymmErrors;
1737   res->SetTitle(graph->GetTitle());
1738   res->GetXaxis()->SetTitle(graph->GetXaxis()->GetTitle());
1739   res->GetYaxis()->SetTitle(graph->GetYaxis()->GetTitle());
1740
1741   for (Int_t i=0; i<graph->GetN(); i++)
1742   {
1743     res->SetPoint(res->GetN(), graph->GetX()[i] + shift, graph->GetY()[i]);
1744     res->SetPointError(res->GetN()-1, graph->GetEX()[i], graph->GetEX()[i], graph->GetEY()[i], graph->GetEY()[i]);
1745     
1746     // asymmetric if space
1747     if (graph->GetEX()[i] > TMath::Abs(shift))
1748     {
1749       res->GetEXlow()[i] += shift;
1750       res->GetEXhigh()[i] -= shift;
1751     }
1752   }
1753   
1754   return res;
1755 }
1756
1757 void GraphRemovePoint(TGraphAsymmErrors* graph1, TGraphAsymmErrors* graph2)
1758 {
1759   if (graph1->GetN() != graph2->GetN())
1760   {
1761     graph1->Sort();
1762     graph2->Sort();
1763     Int_t Nbin = 0;
1764     if (graph1->GetN() < graph2->GetN()) Nbin = graph1->GetN();
1765     else Nbin = graph2->GetN();
1766     for (Int_t bin1 = 0; bin1 < Nbin; bin1++)
1767     {
1768       Float_t x1 = graph1->GetX()[bin1]; 
1769       Float_t x2 = graph2->GetX()[bin1];
1770       if (x1 != x2)
1771       {
1772         if (x1<x2) {graph1->RemovePoint(bin1--);cout << x1 << "\t" << x2 << "\t" << bin1 << "\t" << graph1->GetN() << "\t" << graph2->GetN()<< endl;}
1773         else {graph2->RemovePoint(bin1--);cout << x1 << "\t" << x2 << "\t" << bin1 << "\t" << graph1->GetN() << "\t" << graph2->GetN()<< endl;}
1774       }
1775     }
1776   }
1777 }
1778
1779 void PrepareGraphs(Int_t nHists, TGraphErrors** graph, TGraphErrors** systematicA, TGraphErrors** systematicB, TMultiGraph** multiGraph, TMultiGraph** multiGraphSyst, Int_t uncertaintyID, Float_t offset = 0)
1780 {
1781   Int_t colors[16] =  { 1, 3, 2, 6, 4, 7, 8, 9, 11, 12, 28, 30, 36, 40, 46 };
1782   Int_t markers[16] = { 20, 21, 22, 23, 24, 25, 26, 27, 28, 30, 31, 32, 33, 34, 2, 5};
1783 //  Int_t fillStyle[11] = { 3003, 3003, 3003, 3003, 3003, 3003, 3003, 3003, 3003, 3003, 3003 };
1784   Int_t fillStyle[11] = { 3001, 3002, 3003, 3004, 3005, 3006, 3007, 3008, 3009, 3010, 3011 };
1785   Int_t count = 0;
1786   
1787   if (*multiGraph == 0)
1788     *multiGraph = new TMultiGraph;
1789   if (*multiGraphSyst == 0)
1790     *multiGraphSyst = new TMultiGraph;
1791   
1792   Float_t shift = -2 + offset;
1793   for (Int_t i=0; i<nHists; i++)
1794   {
1795     if (SkipGraph(i))
1796       continue;
1797     
1798     if (SkipGraphForThisPlot(i))
1799     {
1800       count++;
1801       continue;
1802     }
1803     
1804     TGraphAsymmErrors* graphcentrality = FixGraph((TGraphErrors*) graph[i]->Clone(), shift);
1805     if (graphcentrality->GetN() <= 0)
1806       continue;
1807     
1808     if (systematicA)
1809     {
1810       TGraphAsymmErrors* graphsystematicsA = FixGraph((TGraphErrors*) systematicA[i]->Clone(), shift);
1811       
1812       if (graphcentrality->GetN() != graphsystematicsA->GetN())
1813       {
1814         Printf("Different number of points %d %d: %s", graphcentrality->GetN(), graphsystematicsA->GetN(), graphcentrality->GetTitle());
1815         if (graphcentrality->GetN() == 0 || graphsystematicsA->GetN() == 0) continue;
1816         GraphRemovePoint(graphcentrality,graphsystematicsA);
1817       }
1818       
1819       TGraphErrors* graphsystematicsB = 0;
1820       if (systematicB)
1821       {
1822         graphsystematicsB = (TGraphErrors*) systematicB[i]->Clone();
1823         graphsystematicsB->Sort();
1824       
1825         if (graphcentrality->GetN() != graphsystematicsB->GetN())
1826         {
1827           Printf("Different number of points %d %d", graphcentrality->GetN(), graphsystematicsB->GetN());
1828           continue;
1829         }
1830       }
1831       
1832       for (Int_t j=0; j<graphsystematicsA->GetN(); j++)
1833       {
1834         // uncertaintyID
1835         Double_t yMin = graphcentrality->GetY()[j];
1836         Double_t yMax = graphcentrality->GetY()[j];
1837         
1838         if (uncertaintyID == 1)
1839         {
1840           // 5%
1841           yMin *= 0.95;
1842           yMax *= 1.05;
1843         }
1844         else if (uncertaintyID == 2)
1845         {
1846           yMin -= 0.20;
1847           yMax += 0.20;
1848         }
1849         
1850         if (graphsystematicsA->GetY()[j] < graphcentrality->GetY()[j])
1851           yMin = graphcentrality->GetY()[j] - TMath::Sqrt(TMath::Power(yMin - graphcentrality->GetY()[j], 2) + TMath::Power(graphsystematicsA->GetY()[j] - graphcentrality->GetY()[j], 2));
1852         if (graphsystematicsA->GetY()[j] > graphcentrality->GetY()[j])
1853           yMax = graphcentrality->GetY()[j] + TMath::Sqrt(TMath::Power(yMax - graphcentrality->GetY()[j], 2) + TMath::Power(graphsystematicsA->GetY()[j] - graphcentrality->GetY()[j], 2));
1854         
1855         if (graphsystematicsB)
1856         {
1857           if (graphsystematicsB->GetY()[j] < graphcentrality->GetY()[j])
1858             yMin = graphcentrality->GetY()[j] - TMath::Sqrt(TMath::Power(yMin - graphcentrality->GetY()[j], 2) + TMath::Power(graphsystematicsB->GetY()[j] - graphcentrality->GetY()[j], 2));
1859           if (graphsystematicsB->GetY()[j] > graphcentrality->GetY()[j])
1860             yMax = graphcentrality->GetY()[j] + TMath::Sqrt(TMath::Power(yMax - graphcentrality->GetY()[j], 2) + TMath::Power(graphsystematicsB->GetY()[j] - graphcentrality->GetY()[j], 2));
1861         }
1862
1863         graphsystematicsA->GetEYlow()[j] = graphcentrality->GetY()[j] - yMin;
1864         graphsystematicsA->GetEYhigh()[j] = yMax - graphcentrality->GetY()[j];
1865         graphsystematicsA->GetY()[j]  = graphcentrality->GetY()[j];
1866         
1867         graphsystematicsA->GetEXlow()[j] = 1;
1868         graphsystematicsA->GetEXhigh()[j] = graphsystematicsA->GetEXlow()[j];
1869       }
1870       
1871 //       graphsystematicsA->SetFillColor(kGray);
1872       graphsystematicsA->SetFillColor(colors[count]);
1873 //       graphsystematicsA->SetFillStyle(1001);
1874       graphsystematicsA->SetFillStyle(fillStyle[count]);
1875       graphsystematicsA->SetMarkerStyle(0);
1876       graphsystematicsA->SetLineColor(0);
1877       (*multiGraphSyst)->Add(graphsystematicsA, "2");
1878     }
1879     
1880     TString label = graphcentrality->GetTitle();
1881     if (label.Length() > 0)
1882     {
1883       TObjArray* tokens = label.Tokenize("-");
1884       label.Form("%s-%s", tokens->At(0)->GetName(),tokens->At(1)->GetName());
1885       label.ReplaceAll(".00", "");
1886       label.ReplaceAll(".0", "");
1887       label.ReplaceAll("p_{T,trig}", "p_{T,t}");
1888       label.ReplaceAll("p_{T,assoc}", "p_{T,a}");
1889     }
1890     label += " GeV/c";
1891     graphcentrality->SetTitle(label);
1892     Printf("%d %s %d", i, label.Data(), colors[count]);
1893     
1894     graphcentrality->SetMarkerStyle(markers[count]);
1895     graphcentrality->SetMarkerColor(colors[count]);
1896     graphcentrality->SetLineColor(colors[count]);
1897     
1898     (*multiGraph)->Add(graphcentrality);
1899     shift += 1;
1900    
1901     count++;
1902   }
1903 }
1904
1905
1906 Bool_t drawLogo = kFALSE;
1907 const char* MCLabel = 0;
1908
1909 void DrawCentrality(const char* canvasName, Int_t nHists, TGraphErrors** graph, Float_t min = 0, Float_t max = 0, const char* yLabel = "", TGraphErrors** systematicA = 0, TGraphErrors** systematicB = 0, TGraphErrors** graph2 = 0, TGraphErrors** systematic2 = 0, TGraphErrors** graph3 = 0, Int_t uncertaintyID = -1)
1910 {
1911 //   Bool_t found = kTRUE;
1912   TCanvas* c1 = 0;//(TCanvas*) gROOT->GetListOfCanvases()->FindObject(canvasName);
1913   if (!c1)
1914   {
1915     c1 = new TCanvas(canvasName, canvasName, 800, 600);
1916 //     found = kFALSE;
1917   }
1918   c1->cd();
1919   
1920   TMultiGraph* multiGraph = 0;
1921   TMultiGraph* multiGraph2 = 0;
1922   TMultiGraph* multiGraphSyst = 0;
1923   
1924   PrepareGraphs(nHists, graph, systematicA, systematicB, &multiGraph, &multiGraphSyst, uncertaintyID);
1925   if (!multiGraph->GetListOfGraphs())
1926     return;
1927
1928   if (graph2)
1929   {
1930     PrepareGraphs(nHists, graph2, systematic2, 0, &multiGraph2, &multiGraphSyst, uncertaintyID, 0.5);
1931     for (Int_t i=0; i<multiGraph2->GetListOfGraphs()->GetEntries(); i++)
1932       ((TGraph*)multiGraph2->GetListOfGraphs()->At(i))->SetLineWidth(2);
1933     
1934     while (multiGraph2->GetListOfGraphs()->GetEntries() > multiGraph->GetListOfGraphs()->GetEntries())
1935       multiGraph2->GetListOfGraphs()->RemoveAt(multiGraph->GetListOfGraphs()->GetEntries());
1936     
1937     TMultiGraph* multiGraph3 = 0;
1938     if (graph3)
1939     {
1940       PrepareGraphs(nHists, graph3, 0, 0, &multiGraph3, &multiGraphSyst, uncertaintyID);
1941       if (multiGraph3->GetListOfGraphs())
1942       {
1943         for (Int_t i=0; i<multiGraph3->GetListOfGraphs()->GetEntries(); i++)
1944         {
1945           ((TGraph*)multiGraph3->GetListOfGraphs()->At(i))->SetLineWidth(2);
1946           ((TGraph*)multiGraph3->GetListOfGraphs()->At(i))->SetLineStyle(2);
1947         }
1948         
1949         while (multiGraph3->GetListOfGraphs()->GetEntries() > multiGraph->GetListOfGraphs()->GetEntries())
1950           multiGraph3->GetListOfGraphs()->RemoveAt(multiGraph->GetListOfGraphs()->GetEntries());
1951       }
1952     }
1953       
1954 //     multiGraphSyst->Add(multiGraph2, "LX");
1955     
1956 //     multiGraph->Add(multiGraph2);
1957     
1958     if (graph3)
1959       multiGraph2->Add(multiGraph3, "LX");
1960   }
1961
1962 //   multiGraphSyst->Add(multiGraph, (found) ? "LX" : "P");
1963   
1964   // draw by hand, multigraph draw spoils the pngs for some reason, same thing happens with this long code, have changed to gif below...
1965   
1966   TGraphErrors* first = (TGraphErrors*) multiGraph->GetListOfGraphs()->First();
1967   TH2* dummy = new TH2F("dummy", ";Centrality", 100, -2, 115, 100, min, max);
1968   dummy->SetStats(0);
1969   dummy->GetYaxis()->SetTitle(yLabel);
1970   dummy->GetYaxis()->SetTitleOffset(1.2);
1971   dummy->DrawCopy();
1972   
1973   if (multiGraphSyst->GetListOfGraphs())
1974   {
1975     multiGraphSyst->GetListOfGraphs()->First()->Draw((!first) ? "A2" : "2SAME");
1976     if (!first)
1977       first = (TGraphErrors*) multiGraphSyst->GetListOfGraphs()->First();
1978     
1979     for (Int_t i=1; i<multiGraphSyst->GetListOfGraphs()->GetEntries(); i++)
1980       multiGraphSyst->GetListOfGraphs()->At(i)->Draw("2SAME");
1981   }
1982   
1983   if (multiGraph2 && multiGraph2->GetListOfGraphs())
1984   {
1985     for (Int_t i=0; i<multiGraph2->GetListOfGraphs()->GetEntries(); i++)
1986     {
1987       TGraphAsymmErrors* graphTmp = (TGraphAsymmErrors*) multiGraph2->GetListOfGraphs()->At(i);
1988       for (Int_t j=0; j<graphTmp->GetN(); j++)
1989       {
1990         graphTmp->GetEXlow()[j] = 0;
1991         graphTmp->GetEXhigh()[j] = 0;
1992       }
1993       graphTmp->SetMarkerSize(0);
1994     }
1995
1996     multiGraph2->GetListOfGraphs()->First()->Draw((!first) ? "ALZ" : "LXSAME");
1997     if (!first)
1998       first = (TGraphErrors*) multiGraph2->GetListOfGraphs()->First();
1999     
2000     for (Int_t i=1; i<multiGraph2->GetListOfGraphs()->GetEntries(); i++)
2001       multiGraph2->GetListOfGraphs()->At(i)->Draw("LZSAME");
2002   }
2003
2004   if (multiGraph && multiGraph->GetListOfGraphs())
2005   {
2006     multiGraph->GetListOfGraphs()->First()->Draw((!first) ? "AP" : "PSAME");
2007     if (!first)
2008       first = (TGraphErrors*) multiGraph->GetListOfGraphs()->First();
2009     
2010     for (Int_t i=1; i<multiGraph->GetListOfGraphs()->GetEntries(); i++)
2011       multiGraph->GetListOfGraphs()->At(i)->Draw("PSAME");
2012   }
2013
2014 //  TPaveText* paveText = new TPaveText(0.826, 0.059, 0.895, 0.094, "BRNDC");
2015   TPaveText* paveText = new TPaveText(0.77, 0.059, 0.84, 0.094, "BRNDC");
2016   paveText->SetTextSize(0.04);
2017   paveText->SetFillColor(0);
2018   paveText->SetShadowColor(0);
2019   paveText->SetLineColor(0);
2020   paveText->AddText("pp");
2021   paveText->Draw();
2022     
2023   TLegend* legend = new TLegend(0.55, 0.65, 0.95, 0.95);
2024   legend->SetFillColor(0);
2025   legend->SetTextSize(0.03);
2026   
2027   if (graph2)
2028   {
2029     legend->SetNColumns(2);
2030     legend->SetHeader(MCLabel);
2031   }
2032   
2033   for (Int_t i=0; i<multiGraph->GetListOfGraphs()->GetEntries(); i++)
2034   {
2035     if (graph2)
2036     {
2037       legend->AddEntry(multiGraph->GetListOfGraphs()->At(i), " ", "P");
2038       legend->AddEntry(multiGraph2->GetListOfGraphs()->At(i), 0, "L");
2039     }
2040     else
2041       legend->AddEntry(multiGraph->GetListOfGraphs()->At(i), 0, "P");
2042   }
2043   legend->Draw();
2044   
2045   gPad->SetGridx();
2046   gPad->SetGridy();
2047   
2048   DrawLatex(0.13, 0.85,  1, "Pb-Pb #sqrt{s_{NN}} = 2.76 TeV", 0.03);
2049   DrawLatex(0.13, 0.81,  1, "pp #sqrt{s} = 2.76 TeV", 0.03);
2050   DrawLatex(0.13, 0.77, 1, "|#eta| < 0.9", 0.03);
2051   TString text;
2052   TString text2;
2053   if (TString(canvasName).BeginsWith("sigma_phi") || TString(canvasName).BeginsWith("kurtosisphi"))
2054   {
2055     text = "Projected within |#Delta#eta| < 0.80";
2056     text2 = "Calculated within |#Delta#varphi| < 0.87";
2057   }
2058   if (TString(canvasName).BeginsWith("sigma_eta") || TString(canvasName).BeginsWith("kurtosiseta"))
2059   {
2060     text = "Projected within |#Delta#varphi| < 0.87";
2061     text2 = "Calculated within |#Delta#eta| < 0.80";
2062   }
2063   if (text.Length() > 0)
2064   {
2065     DrawLatex(0.13, 0.73, 1, text, 0.03);
2066     DrawLatex(0.13, 0.69, 1, text2, 0.03);
2067   }
2068   
2069   if (0 && MCLabel)
2070   {
2071     DrawLatex(0.13, 0.18, 1, "Points: Data", 0.03);
2072     DrawLatex(0.13, 0.14, 1, Form("Lines: %s", MCLabel), 0.03);
2073   }
2074   
2075   if (drawLogo)
2076     DrawALICELogo(kTRUE, 0.41, 0.75, 0.54, 0.95);
2077   
2078   gSystem->mkdir(fgFolder, kTRUE);
2079   c1->SaveAs(Form("%s/%s.png", fgFolder.Data(), canvasName));
2080 //   c1->SaveAs(Form("%s/%s.eps", fgFolder.Data(), canvasName));
2081 }
2082
2083 void CalculateRMSSigma(TGraphErrors*** graphsTmp = 0)
2084 {
2085   if (!graphsTmp)
2086     graphsTmp = graphs;
2087   
2088   CalculateRMS(NHists, graphsTmp[1], graphsTmp[4], graphsTmp[3]);
2089   CalculateRMS(NHists, graphsTmp[2], graphsTmp[5], graphsTmp[3]);
2090
2091   CalculateRMS(NHists, graphsTmp[1+16], graphsTmp[4+16], graphsTmp[3+16]);
2092   CalculateRMS(NHists, graphsTmp[2+16], graphsTmp[5+16], graphsTmp[3+16]);
2093
2094   // sqrt(moment2) = sigma
2095   SqrtAll2(NHists, graphsTmp[8], graphsTmp[8+16]);
2096   SqrtAll2(NHists, graphsTmp[9], graphsTmp[9+16]);
2097 }
2098
2099 // void CombineSyst(const char* systFile)
2100 // {
2101 //   TGraphErrors*** current = graphs;
2102 //   
2103 //   ReadGraphs(systFile);
2104 //   CalculateRMSSigma();
2105 //   TGraphErrors*** graphsSyst = graphs;
2106 //   graphs = current;
2107 //   
2108 //   const Int_t NGraphList = 6;
2109 //   Int_t graphList[] = { 1, 2, 8, 9, 14, 15 };
2110 //   for (Int_t i=0; i<NGraphList; i++)
2111 //   {
2112 //     for (Int_t j=0; j<NHists; j++)
2113 //     {
2114 //       
2115 //     graphs[i+16][j]
2116 //     graphsSyst[i][j]
2117 //     }
2118 //   }
2119 //   
2120 // 
2121 // }
2122
2123 void DrawResultsCentrality(const char* fileName = "graphs.root", const char* fileNameWingRemoved = 0, Int_t offset = 0)
2124 {
2125   TGraphErrors*** graphsWingRemoved = 0;
2126   if (fileNameWingRemoved)
2127   {
2128     ReadGraphs(fileNameWingRemoved);
2129     graphsWingRemoved = graphs;
2130   }
2131   
2132   ReadGraphs(fileName);
2133   
2134   Int_t nHists = 24; //NHists;
2135   
2136   if (1)
2137   {
2138 /*    DrawCentrality("norm", nHists, graphs[0+offset], 0, 0.2, "N (a.u.)", (graphsWingRemoved) ? graphsWingRemoved[0+offset] : 0);
2139     
2140 //     return;
2141     
2142     DrawCentrality("width_phi1_centrality", nHists, graphs[1+offset], 0, 0.8, "#sigma_{#Delta#varphi, 1} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[1+offset] : 0);
2143 //     return;
2144     DrawCentrality("width_phi2_centrality", nHists, graphs[4+offset], 0, 0.8, "#sigma_{#Delta#varphi, 2} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[4+offset] : 0);
2145     DrawCentrality("width_eta1_centrality", nHists, graphs[2+offset], 0, 0.8, "#sigma_{#Delta#eta, 1} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[2+offset] : 0);
2146     DrawCentrality("width_eta2_centrality", nHists, graphs[5+offset], 0, 0.8, "#sigma_{#Delta#eta, 2} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[5+offset] : 0);
2147 */    
2148 /*    DrawCentrality("norm_phi(1D)", nHists, graphs[35], 0, 1, "N (a.u.)", (graphsWingRemoved) ? graphsWingRemoved[0+offset] : 0);
2149     DrawCentrality("norm_eta(1D)", nHists, graphs[39], 0, 1, "N (a.u.)", (graphsWingRemoved) ? graphsWingRemoved[0+offset] : 0);
2150 */    
2151 //     DrawCentrality("width_phi1_centrality(1D)", nHists, graphs[37], 0, 0.8, "#sigma_{#Delta#varphi, 1} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[1+offset] : 0);
2152 //    DrawCentrality("width_phi2_centrality(1D)", nHists, graphs[38], 0, 0.8, "#sigma_{#Delta#varphi, 2} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[4+offset] : 0);
2153 //     DrawCentrality("width_eta1_centrality(1D)", nHists, graphs[41], 0, 0.8, "#sigma_{#Delta#eta, 1} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[2+offset] : 0);
2154 //    DrawCentrality("width_eta2_centrality(1D)", nHists, graphs[42], 0, 0.8, "#sigma_{#Delta#eta, 2} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[5+offset] : 0);
2155 //    DrawCentrality("norm2phi(1D)", nHists, graphs[36], 0, 1.5, "N (a.u.)", (graphsWingRemoved) ? graphsWingRemoved[0+offset] : 0);
2156 //    DrawCentrality("norm2eta(1D)", nHists, graphs[40], 0, 1.5, "N (a.u.)", (graphsWingRemoved) ? graphsWingRemoved[0+offset] : 0);
2157
2158     
2159     Bool_t uncertainties = kFALSE;
2160
2161     CalculateRMSSigma();
2162     if (graphsWingRemoved)
2163       CalculateRMSSigma(graphsWingRemoved);
2164 //     DrawCentrality("phi_rms_centrality_old", nHists, graphs[1+offset], 0, 0.8, Form("#sigma_{#Delta#varphi} (%s) (rad.)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[1+offset] : 0, 0, 0, 0, 0, (uncertainties) ? 1 : -1);
2165 //     DrawCentrality("eta_rms_centrality_old", nHists, graphs[2+offset], 0, 0.8 + 0.4 / 16 * offset, Form("#sigma_{#Delta#eta} (%s)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[2+offset] : 0, 0, 0, 0, 0, (uncertainties) ? 1 : -1);
2166
2167     DrawCentrality("phi_rms_centrality", nHists, graphs[10+offset], 0, 0.8, Form("#sigma_{#Delta#varphi} (%s) (rad.)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[10+offset] : 0, 0, 0, 0, 0, (uncertainties) ? 1 : -1);
2168     DrawCentrality("eta_rms_centrality", nHists, graphs[11+offset], 0, 0.8 + 0.4 / 16 * offset, Form("#sigma_{#Delta#eta} (%s)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[11+offset] : 0, 0, 0, 0, 0, (uncertainties) ? 1 : -1);
2169
2170 //     return;
2171 //   DrawCentrality("phi_rms_centrality(1D)", nHists, graphs[37], 0, 0.8, Form("#sigma_{#Delta#varphi} (%s) (rad.)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[1+offset] : 0, 0, 0, 0, 0, 1);
2172
2173 //   DrawCentrality("eta_rms_centrality(1D)", nHists, graphs[41], 0, 0.8 + 0.4 / 16 * offset, Form("#sigma_{#Delta#eta} (%s)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[2+offset] : 0, 0, 0, 0, 0, 1);
2174 //   DrawCentrality("chi2_phi(1D)", nHists, graphs[43], 0, 150, "#chi^{2}/ndf");
2175 //   DrawCentrality("chi2_eta(1D)", nHists, graphs[44], 0, 150, "#chi^{2}/ndf");
2176     DrawCentrality("chi2_1", nHists, graphs[6+offset], 0.5, 5, "#chi^{2}/ndf (full region)", (graphsWingRemoved) ? graphsWingRemoved[6+offset] : 0);
2177     DrawCentrality("chi2_2", nHists, graphs[7+offset], 0.5, 5, "#chi^{2}/ndf (peak region)", (graphsWingRemoved) ? graphsWingRemoved[7+offset] : 0);
2178     
2179     DrawCentrality("sigma_phi", nHists, graphs[8+offset], 0, 1, "#sigma_{#Delta#varphi} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[8+offset] : 0, 0, 0, 0, 0, (uncertainties) ? 1 : -1);
2180     DrawCentrality("sigma_eta", nHists, graphs[9+offset], 0, 1, "#sigma_{#Delta#eta}", (graphsWingRemoved) ? graphsWingRemoved[9+offset] : 0, 0, 0, 0, 0, (uncertainties) ? 1 : -1);
2181
2182 /*    
2183 //     DrawCentrality("moment3_phi", nHists, graphs[10+offset], -0.5, 0.5, "moment3 #varphi (rad.)", (graphsWingRemoved) ? graphsWingRemoved[10+offset] : 0, 0, 0, 0, 0, 0);
2184 //     DrawCentrality("moment3_eta", nHists, graphs[11+offset], -0.5, 0.5, "moment3 #eta", (graphsWingRemoved) ? graphsWingRemoved[11+offset] : 0, 0, 0, 0, 0, 0);
2185 */
2186 /*    DrawCentrality("IAAFit", nHists, graphs[32+offset], 0, 20, "I_AA", (graphsWingRemoved) ? graphsWingRemoved[32+offset] : 0);
2187     DrawCentrality("Yield", nHists, graphs[33+offset], 0, 20, "Yield", (graphsWingRemoved) ? graphsWingRemoved[33+offset] : 0);
2188     DrawCentrality("IAAHist", nHists, graphs[34+offset], 0, 20, "I_AA", (graphsWingRemoved) ? graphsWingRemoved[34+offset] : 0);
2189     DrawCentrality("IAA", nHists, graphs[34+offset], 0, 20, "I_AA", graphs[32+offset]);
2190 */
2191   }
2192   
2193  DrawCentrality("kurtosisphi_centrality", 12, graphs[14+offset], -1.5, 2.5, "Kurtosis #Delta#varphi", (graphsWingRemoved) ? graphsWingRemoved[14+offset] : 0, 0, 0, 0, 0, 2);
2194  DrawCentrality("kurtosiseta_centrality", 12, graphs[15+offset], -1.5, 2.5, "Kurtosis #Delta#eta", (graphsWingRemoved) ? graphsWingRemoved[15+offset] : 0, 0, 0, 0, 0, 2);
2195 }
2196
2197 /*
2198 void DrawResultsCentrality(const char* fileName = "graphs.root", const char* fileNameWingRemoved = 0, Int_t offset = 0)
2199 {
2200   TGraphErrors*** graphsWingRemoved = 0;
2201   if (fileNameWingRemoved)
2202   {
2203     ReadGraphs(fileNameWingRemoved);
2204     graphsWingRemoved = graphs;
2205   }
2206   
2207   ReadGraphs(fileName);
2208   
2209   Int_t nHists = 13; //NHists;
2210   
2211   if (1)
2212   {
2213     DrawCentrality("norm", nHists, graphs[0], 0, 0.2, "N (a.u.)", graphs[0+16], (graphsWingRemoved) ? graphsWingRemoved[0] : 0);
2214     
2215 //     return;
2216     
2217     DrawCentrality("width_phi1_centrality", nHists, graphs[1], 0, 0.8, "#sigma_{#varphi, 1} (rad.)", graphs[1+16], (graphsWingRemoved) ? graphsWingRemoved[1] : 0);
2218 //     return;
2219     DrawCentrality("width_phi2_centrality", nHists, graphs[4], 0, 0.8, "#sigma_{#varphi, 2} (rad.)", graphs[4+16], (graphsWingRemoved) ? graphsWingRemoved[4] : 0);
2220     DrawCentrality("width_eta1_centrality", nHists, graphs[2], 0, 0.8, "#sigma_{#eta, 1} (rad.)", graphs[2+16], (graphsWingRemoved) ? graphsWingRemoved[2] : 0);
2221     DrawCentrality("width_eta2_centrality", nHists, graphs[5], 0, 0.8, "#sigma_{#eta, 2} (rad.)", graphs[5+16], (graphsWingRemoved) ? graphsWingRemoved[5] : 0);
2222     
2223     CalculateRMSSigma();
2224     if (graphsWingRemoved)
2225       CalculateRMSSigma(graphsWingRemoved);
2226
2227     DrawCentrality("phi_rms_centrality", nHists, graphs[1], 0, 0.8, "#sigma_{#varphi} (fit) (rad.)", graphs[1+16], (graphsWingRemoved) ? graphsWingRemoved[1] : 0, 0, 0, 0, 1);
2228
2229     DrawCentrality("eta_rms_centrality", nHists, graphs[2], 0, 0.8, "#sigma_{#eta} (fit)", graphs[2+16], (graphsWingRemoved) ? graphsWingRemoved[2] : 0, 0, 0, 0, 1);
2230
2231     DrawCentrality("chi2_1", nHists, graphs[6], 0.5, 2.5, "#chi^{2}/ndf (full region)");
2232     DrawCentrality("chi2_2", nHists, graphs[7], 0.5, 2.5, "#chi^{2}/ndf (peak region)");
2233     
2234     DrawCentrality("sigma_phi", nHists, graphs[8], 0.1, 0.5, "#sigma_{#varphi} (rad.)", graphs[8+16], (graphsWingRemoved) ? graphsWingRemoved[8] : 0, 0, 0, 0, 1);
2235     DrawCentrality("sigma_eta", nHists, graphs[9], 0.1, 0.5, "#sigma_{#eta}", graphs[9+16], (graphsWingRemoved) ? graphsWingRemoved[9] : 0, 0, 0, 0, 1);
2236   }
2237   
2238   DrawCentrality("kurtosisphi_centrality", 12, graphs[14], -2, 2, "Kurtosis #varphi", graphs[14+16], (graphsWingRemoved) ? graphsWingRemoved[14] : 0, 0, 0, 0, 2);
2239   DrawCentrality("kurtosiseta_centrality", 12, graphs[15], -2, 2, "Kurtosis #eta", graphs[15+16], (graphsWingRemoved) ? graphsWingRemoved[15] : 0, 0, 0, 0, 2);
2240 }
2241 */
2242
2243 void MCComparison(const char* fileNameData, const char* fileNameWingRemoved, const char* fileNameHijing, const char* fileNameAMPT, Int_t offset = 0)
2244 {
2245   Int_t nHists = 14; //NHists;
2246
2247   ReadGraphs(fileNameWingRemoved);
2248   CalculateRMSSigma();
2249   TGraphErrors*** graphsWingRemoved = graphs;
2250
2251   ReadGraphs(fileNameHijing);
2252   CalculateRMSSigma();
2253   TGraphErrors*** graphs1 = graphs;
2254
2255   CreateGraphStructure();
2256   TGraphErrors*** graphs2 = graphs;
2257   if (fileNameAMPT)
2258   {
2259     ReadGraphs(fileNameAMPT);
2260     CalculateRMSSigma();
2261     graphs2 = graphs;
2262   }
2263
2264   ReadGraphs(fileNameData);
2265   CalculateRMSSigma();
2266
2267   if (0)
2268   {
2269     Int_t graphList[] = { 10, 11, 8, 9, 14, 15 };
2270     for (Int_t i=0; i<6; i++)
2271     {
2272       graphs[graphList[i]+offset][5] = new TGraphErrors;
2273       graphs[graphList[i]+offset][10] = new TGraphErrors;
2274       graphs1[graphList[i]+offset][5] = new TGraphErrors;
2275       graphs1[graphList[i]+offset][10] = new TGraphErrors;
2276     }
2277   }
2278   
2279   DrawCentrality("phi_rms_centrality_mc", nHists, graphs[10+offset], 0, 0.8, Form("#sigma_{#Delta#varphi} (%s) (rad.)", fitLabel), graphsWingRemoved[1+offset], 0, graphs1[10+offset], 0, graphs2[10+offset], 1);
2280   
2281 //   return;
2282   
2283   DrawCentrality("eta_rms_centrality_mc", nHists, graphs[11+offset], 0, 0.8 + 0.4 / 16 * offset, Form("#sigma_{#Delta#eta} (%s)", fitLabel), graphsWingRemoved[11+offset], 0, graphs1[11+offset], 0, graphs2[2+offset], 1);
2284   
2285   DrawCentrality("sigma_phi_mc", nHists, graphs[8+offset], 0, 0.6, "#sigma_{#Delta#varphi} (rad.)", graphsWingRemoved[8+offset], 0, graphs1[8+offset], 0, graphs2[8+offset], 1);
2286   DrawCentrality("sigma_eta_mc", nHists, graphs[9+offset], 0, 0.6, "#sigma_{#Delta#eta}", graphsWingRemoved[9+offset], 0, graphs1[9+offset], 0, graphs2[9+offset], 1);
2287
2288   DrawCentrality("kurtosisphi_centrality_mc", nHists, graphs[14+offset], -1.5, 3.5, "Kurtosis #Delta#varphi", graphsWingRemoved[14+offset], 0, graphs1[14+offset], 0, graphs2[14+offset], 2);
2289   DrawCentrality("kurtosiseta_centrality_mc", nHists, graphs[15+offset], -1.5, 3.5, "Kurtosis #Delta#eta", graphsWingRemoved[15+offset], 0, graphs1[15+offset], 0, graphs2[15+offset], 2);
2290
2291 /*  DrawCentrality("phi_rms_centrality_mc", nHists, graphs[1], 0, 0.8, "#sigma_{#varphi} (fit) (rad.)", graphs[1+16], graphsWingRemoved[1],  graphs1[1], 0, graphs2[1], 1);
2292   DrawCentrality("eta_rms_centrality_mc", nHists, graphs[2], 0, 0.8, "#sigma_{#eta} (fit)", graphs[2+16], graphsWingRemoved[2], graphs1[2], 0, graphs2[2], 1);
2293
2294   DrawCentrality("sigma_phi", nHists, graphs[8], 0.1, 0.5, "#sigma_{#varphi} (rad.)", graphs[8+16], graphsWingRemoved[8], graphs1[8], 0, graphs2[8], 1);
2295   DrawCentrality("sigma_eta", nHists, graphs[9], 0.1, 0.5, "#sigma_{#eta}", graphs[9+16], graphsWingRemoved[9], graphs1[9], 0, graphs2[9], 1);
2296
2297   DrawCentrality("kurtosisphi_centrality_mc", nHists, graphs[14], -2, 2, "Kurtosis #varphi", graphs[14+16], graphsWingRemoved[14], graphs1[14], 0, graphs2[14], 2);
2298   DrawCentrality("kurtosiseta_centrality_mc", nHists, graphs[15], -2, 2, "Kurtosis #eta", graphs[15+16], graphsWingRemoved[15], graphs1[15], 0, graphs2[15], 2);*/
2299 }
2300
2301 void ShowWingEffect(const char* fileNameData, const char* fileNameWingRemoved, Int_t offset = 0)
2302 {
2303   Int_t nHists = 24; //NHists;
2304
2305   ReadGraphs(fileNameWingRemoved);
2306   CalculateRMSSigma();
2307   TGraphErrors*** graphs1 = graphs;
2308
2309   ReadGraphs(fileNameData);
2310   CalculateRMSSigma();
2311   
2312   DrawCentrality("phi_rms_centrality_mc", nHists, graphs[1+offset], 0, 0.8, "#sigma_{#varphi} (fit) (rad.)", graphs1[1+offset]);
2313   DrawCentrality("eta_rms_centrality_mc", nHists, graphs[2+offset], 0, 0.8, "#sigma_{#eta} (fit)", graphs1[2+offset]);
2314
2315   DrawCentrality("sigma_phi", nHists, graphs[8+offset], 0, 0.8, "#sigma_{#varphi} (rad.)", graphs1[8+offset]);
2316   DrawCentrality("sigma_eta", nHists, graphs[9+offset], 0, 0.8, "#sigma_{#eta}", graphs1[9+offset]);
2317
2318   DrawCentrality("kurtosisphi_centrality_mc", nHists, graphs[14+offset], -2, 4, "Kurtosis #varphi", graphs1[14+offset]);
2319   DrawCentrality("kurtosiseta_centrality_mc", nHists, graphs[15+offset], -2, 4, "Kurtosis #eta", graphs1[15+offset]);
2320 }
2321
2322 void ShowFitEffect(const char* fileNameData)
2323 {
2324   Int_t nHists = NHists;
2325
2326   ReadGraphs(fileNameData);
2327   CalculateRMSSigma();
2328   
2329   DrawCentrality("phi_rms_centrality_mc", nHists, graphs[1], 0, 0.8, "#sigma_{#varphi} (fit) (rad.)", graphs[1+16]);
2330   DrawCentrality("eta_rms_centrality_mc", nHists, graphs[2], 0, 0.8, "#sigma_{#eta} (fit)", graphs[2+16]);
2331
2332   DrawCentrality("sigma_phi", nHists, graphs[8], 0, 0.8, "#sigma_{#varphi} (rad.)", graphs[8+16]);
2333   DrawCentrality("sigma_eta", nHists, graphs[9], 0, 0.8, "#sigma_{#eta}", graphs[9+16]);
2334
2335   DrawCentrality("kurtosisphi_centrality_mc", nHists, graphs[14], -2, 4, "Kurtosis #varphi", graphs[14+16]);
2336   DrawCentrality("kurtosiseta_centrality_mc", nHists, graphs[15], -2, 4, "Kurtosis #eta", graphs[15+16]);
2337
2338   DrawCentrality("chi2_1", nHists, graphs[6], 0.5, 5, "chi2_1", graphs[6+16]);
2339   DrawCentrality("chi2_2", nHists, graphs[7], 0.5, 5, "chi2_2", graphs[7+16]);
2340 }
2341
2342 void CompareSigmas(const char* fileNameData)
2343 {
2344   // compare sigma from fit with sigma from direct calculation but taking for the first the limited range of 0.8 into account
2345   
2346   Int_t nHists = 12; //NHists;
2347
2348   ReadGraphs(fileNameData);
2349   CalculateRMSSigma();
2350   TGraphErrors*** graphsOriginal = graphs;
2351   
2352   ReadGraphs(fileNameData);
2353   
2354   // rms
2355   for (Int_t histId = 0; histId < nHists; histId++)
2356   {
2357     if (!graphs[0][histId])
2358       continue;
2359     
2360     for (Int_t i=0; i<graphs[0][histId]->GetN(); i++)
2361     {
2362       // phi
2363       TF1* func = new TF1("func", "[0]+[1]*([4]/[2]/TMath::Sqrt(TMath::TwoPi())*exp(-0.5*(x/[2])**2)*TMath::Erf([7]/TMath::Sqrt(2)/[3])+(1-[4])/[5]/TMath::Sqrt(TMath::TwoPi())*exp(-0.5*(x/[5])**2)*TMath::Erf([7]/TMath::Sqrt(2)/[6]))", -0.5 * TMath::Pi(), 0.5 * TMath::Pi());
2364       func->SetParameter(0, 0);
2365       for (Int_t k=0; k<6; k++)
2366         func->SetParameter(k+1, graphs[k][histId]->GetY()[i]);
2367       func->SetParameter(7, 0.8); // project Limit
2368       
2369       graphs[8][histId]->GetY()[i] = TMath::Sqrt(func->CentralMoment(2, -0.87, 0.87));
2370       
2371       // eta
2372       func = new TF1("func", "[0]+[1]*([4]/[3]/TMath::Sqrt(TMath::TwoPi())*exp(-0.5*(x/[3])**2)*TMath::Erf([7]/TMath::Sqrt(2)/[2])+(1-[4])/[6]/TMath::Sqrt(TMath::TwoPi())*exp(-0.5*(x/[6])**2)*TMath::Erf([7]/TMath::Sqrt(2)/[5]))", -kOuterLimit, kOuterLimit);
2373       func->SetParameter(0, 0);
2374       for (Int_t k=0; k<6; k++)
2375         func->SetParameter(k+1, graphs[k][histId]->GetY()[i]);
2376       func->SetParameter(7, 0.87); // project Limit
2377       
2378       graphs[9][histId]->GetY()[i] = TMath::Sqrt(func->CentralMoment(2, -0.8, 0.8));
2379     }
2380   }
2381   
2382   DrawCentrality("phi_comparison1", nHists, graphsOriginal[1], 0.2, 0.8, "#sigma_{#Delta#varphi} (fit)", 0, 0, graphs[8]);
2383   DrawCentrality("phi_comparison2", nHists, graphsOriginal[8], 0.2, 0.5, "#sigma_{#Delta#varphi} (fit)", 0, 0, graphs[8]);
2384   DrawCentrality("eta_comparison1", nHists, graphsOriginal[2], 0.2, 0.8, "#sigma_{#Delta#eta} (fit)", 0, 0, graphs[9]);
2385   DrawCentrality("eta_comparison2", nHists, graphsOriginal[9], 0.2, 0.5, "#sigma_{#Delta#eta} (fit)", 0, 0, graphs[9]);
2386 }  
2387
2388 Float_t** ExtractSystematics(const char* baseFile, const char* systFile, Int_t offset)
2389 {
2390   ReadGraphs(baseFile);
2391   CalculateRMSSigma();
2392  
2393   TGraphErrors*** graphsBase = graphs;
2394   if (systFile)
2395   {
2396     ReadGraphs(systFile);
2397     CalculateRMSSigma();
2398   }
2399   
2400   // calculate syst unc for these graphs
2401   const Int_t NGraphList = 6;
2402   Int_t graphList[] = { 1, 2, 8, 9, 14, 15 };
2403
2404   Float_t** results = new Float_t*[NGraphList];
2405
2406   for (Int_t i=0; i<NGraphList; i++)
2407   {
2408     results[i] = new Float_t[2];
2409
2410     Bool_t kurtosis = kFALSE;
2411     if (graphList[i] == 14 || graphList[i] == 15)
2412       kurtosis = kTRUE;
2413     
2414     TH1* hist = 0;
2415     if (!kurtosis)
2416       hist = new TH1F(Form("hist_%d", i), "", 50, 0.5, 1.5);  
2417     else
2418       hist = new TH1F(Form("hist_%d", i), "", 50, -1, 1);  
2419   
2420     TCanvas* c = new TCanvas(Form("%s_%d_%d", systFile, i, 0), Form("%s_%d_%d", systFile, i, 0), 1000, 1000);
2421     c->Divide(4, 4);
2422
2423     Int_t count = 1;
2424     for (Int_t j=0; j<NHists; j++)
2425     {
2426       if (SkipGraph(j))
2427         continue;
2428       
2429       // for kurtosis we show only up to 12
2430       if (kurtosis)
2431         if (j >= 12)
2432           continue;
2433       
2434       TGraphErrors* graph1 = graphsBase[graphList[i]+offset][j];
2435       TGraphErrors* graph2 = (systFile) ? graphs[graphList[i]+offset][j] : graphsBase[graphList[i]+16][j];
2436       
2437       if (graph1->GetN() == 0)
2438         continue;
2439       
2440       graph1->Sort();
2441       graph2->Sort();
2442       
2443       c->cd(count++);
2444       gPad->SetGridx();
2445       gPad->SetGridy();
2446       graph1->SetMarkerStyle(24);
2447       graph1->DrawClone("AP");
2448       graph2->SetLineColor(2);
2449       graph2->SetMarkerColor(2);
2450       graph2->SetMarkerStyle(25);
2451       graph2->DrawClone("LSAME");
2452       
2453       c->cd(count++);
2454       gPad->SetGridx();
2455       gPad->SetGridy();
2456       
2457       // reset errors on graph2 
2458 /*      for (Int_t k=0; k<graph2->GetN(); k++)
2459         graph2->GetEY()[k] = 0;*/
2460       
2461       if (!kurtosis)
2462       {
2463         DivideGraphs(graph1, graph2);
2464         ((TGraphErrors*) graph1->DrawClone("AP"))->GetYaxis()->SetRangeUser(0.8, 1.2);
2465       }
2466       else
2467       {
2468         SubtractGraphs(graph1, graph2);
2469         ((TGraphErrors*) graph1->DrawClone("AP"));
2470       } 
2471       
2472       for (Int_t k=0; k<graph1->GetN(); k++)
2473       {
2474 //      hist->Fill(graph1->GetY()[k], 1.0 / (graph1->GetEY()[k] / graph1->GetY()[k]));
2475         if (kurtosis)
2476           hist->Fill(TMath::Abs(graph1->GetY()[k]));
2477         else
2478           hist->Fill(graph1->GetY()[k]);
2479       }
2480       
2481       if (count == 37)
2482         break;
2483     }
2484     
2485     c->SaveAs(Form("syst/%s_%d.eps", systFile, i));
2486     
2487     new TCanvas;
2488     hist->Draw();
2489     hist->Sumw2();
2490     
2491     Float_t mean = -1;
2492     Float_t sigma = -1;
2493     if (0)
2494     {
2495       hist->Fit("gaus", "Q");
2496       mean = hist->GetFunction("gaus")->GetParameter(1);
2497       sigma = hist->GetFunction("gaus")->GetParameter(2);
2498       
2499       Printf("%d: %.3f +- %.3f %.3f | %.3f +- %.3f", i, mean, sigma, sigma / TMath::Sqrt(hist->GetEntries()), hist->GetMean(), hist->GetMean(11));
2500     }
2501     else
2502     {
2503       mean = hist->GetMean(1);
2504       sigma = hist->GetMean(11);
2505
2506       Printf("%d: %.3f +- %.3f", i, hist->GetMean(), hist->GetMean(11));
2507     }
2508     
2509     if (!kurtosis)
2510     {
2511       mean -= 1;
2512       mean *= 100;
2513       sigma *= 100;
2514     }
2515     results[i][0] = mean;
2516     results[i][1] = sigma / TMath::Sqrt(hist->GetEntries());
2517   }
2518   
2519   return results;
2520 }
2521
2522 void GetProjections(TH2* hist, TH1** projPhi, TH1** projEta, Int_t k)
2523 {
2524   Float_t projectLimit = 0.8;
2525
2526   *projPhi = hist->ProjectionX(Form("%s_proj1%d", hist->GetName(), k), hist->GetYaxis()->FindBin(-projectLimit+0.01), hist->GetYaxis()->FindBin(projectLimit-0.01));
2527     
2528   *projEta = hist->ProjectionY(Form("%s_proj2_%d", hist->GetName(), k), hist->GetXaxis()->FindBin(-projectLimit+0.01), hist->GetXaxis()->FindBin(projectLimit-0.01));
2529 }
2530
2531 // uncorrected ones
2532 const char* systBaseFile = "dphi_corr_2d_120429.root";
2533 const char* systTrackCuts1 = "dphi_corr_120425_hybrid.root";
2534 const char* systTrackCuts2 = "dphi_corr_120430_raa.root";
2535 const char* systVertex = "dphi_corr_2d_120425_vertex.root";
2536 const char* systResonances = "dphi_corr_120502_resonances.root";
2537 const char* systTTR = "dphi_corr_2d_120430_widettr.root";
2538
2539 // corrected ones
2540 const char* systBaseFileCorrected = "dphi_corr_2d_120508.root";
2541 const char* systTrackCuts1Corrected = "dphi_corr_120507_hybrid.root";
2542 const char* systTrackCuts2Corrected = "dphi_corr_120507_raa.root";
2543 const char* systWingRemovedCorrected = "dphi_corr_2d_120508_wingremoved.root";
2544
2545 void ExtractSystematicsProjections(Int_t mode, Float_t rangeBegin = -0.4, Float_t rangeEnd = 0.4, Bool_t draw = kFALSE)
2546 {
2547   if (!draw)
2548     gROOT->SetBatch(kTRUE);
2549   
2550   Int_t maxLeadingPt = 4;
2551   Int_t maxAssocPt = 6;
2552
2553   Int_t NEffects = 3;
2554   
2555   TFile* file1 = 0;
2556   const char** systFiles = 0;
2557   if (mode == 2)
2558   {
2559     NEffects = 5;
2560     file1 = TFile::Open(systBaseFile);
2561     static const char* systFilesTmp[5] = { systVertex, systResonances, systTTR, systTrackCuts1Corrected, systTrackCuts2Corrected };
2562     systFiles = systFilesTmp;
2563   }
2564   else if (mode == 0)
2565   {
2566     file1 = TFile::Open(systBaseFile);
2567     static const char* systFilesTmp[3] = { systVertex, systResonances, systTTR };
2568     systFiles = (const char**) systFilesTmp;
2569   }
2570   else if (mode == 1)
2571   {
2572     NEffects = 3;
2573     file1 = TFile::Open(systBaseFileCorrected);
2574     static const char* systFilesTmp[3] = { systTrackCuts1Corrected, systTrackCuts2Corrected, systWingRemovedCorrected };
2575     systFiles = systFilesTmp;
2576   }
2577   else if (mode == 3)
2578   {
2579     NEffects = 3;
2580     file1 = TFile::Open(systBaseFileCorrected);
2581     static const char* systFilesTmp[3] = { systBaseFileCorrected, systBaseFileCorrected, systBaseFileCorrected };
2582     systFiles = systFilesTmp;
2583   }
2584   
2585   TH2* uncertainty[4];
2586   for (Int_t i=0; i<4; i++)
2587     uncertainty[i] = new TH2F(Form("uncertainty_%d", i), Form("%s %s;ptbin;centrality", (i % 2 == 0) ? "#varphi" : "#eta", (i < 2) ?"average" : "deviation"), 30, 0, 30, 6, 0, 6);
2588   
2589   Int_t count = 0;
2590   Int_t ptBin = 0;
2591   for (Int_t i=0; i<maxLeadingPt-1; i++)
2592   {
2593     for (Int_t j=1; j<maxAssocPt-1; j++)
2594     {
2595       if (i == 0 && j == 2)
2596         continue;
2597       if (i == 1 && j == 3)
2598         continue;
2599       if (i == 2 && j == 4)
2600         continue;
2601       
2602       ptBin++;
2603       for (Int_t histId = 0; histId < NHists; histId++)
2604       {
2605         TH2* hist1 = (TH2*) file1->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
2606         if (!hist1)
2607           continue;
2608         
2609         if (hist1->GetEntries() < 1e4)
2610         {
2611           Printf("Only %f entries. Skipping...", hist1->GetEntries());
2612           continue;
2613         }
2614         
2615         Printf("%d %d %d %s", i, j, histId, hist1->GetTitle());
2616         
2617         TH1* proj1[2];
2618         TH1* proj2[10][2];
2619         
2620         SubtractEtaGap(hist1, kEtaLimit, kOuterLimit, kFALSE);
2621         GetProjections(hist1, &proj1[0], &proj1[1], count++);
2622
2623         TCanvas* c = new TCanvas(Form("c_%d_%d_%d", mode, ptBin, histId), Form("c_%d_%d_%d", mode, ptBin, histId), 1200, 1000);
2624         c->Divide(3, 2);
2625           
2626         for (Int_t k=0; k<2; k++)
2627         {
2628           c->cd(1+k*3);
2629           proj1[k]->SetStats(0);
2630           proj1[k]->GetXaxis()->SetRangeUser(-1.5, 1.5);
2631           proj1[k]->GetYaxis()->SetRangeUser(proj1[k]->GetMinimum() * 1.1, proj1[k]->GetMaximum() * 1.4);
2632           proj1[k]->DrawCopy();
2633         }
2634         
2635         Double_t maxAverage[2] = { 0, 0 };
2636         Double_t maxDev[2] = { 0, 0 };
2637   
2638         for (Int_t n=0; n<NEffects; n++)
2639         {
2640           if (histId == 2 && mode == 0 && n > 0)
2641             continue;
2642           if (histId == 2 && mode == 1 && n > 0)
2643             continue;
2644 //        Printf("%d %s", n, systFiles[n]);
2645           TFile* file2 = TFile::Open(systFiles[n]);
2646           hist1 = (TH2*) file2->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
2647           if (!hist1)
2648             continue;
2649           
2650           if (mode == 3 && n == 0) // change eta limits
2651             SubtractEtaGap(hist1, kEtaLimit-0.2, kOuterLimit, kFALSE);
2652           if (mode == 3 && n == 1) // change eta limits
2653             SubtractEtaGap(hist1, kEtaLimit+0.2, kOuterLimit, kFALSE);
2654           if (mode == 3 && n == 2) // change eta limits
2655             SubtractEtaGap(hist1, kEtaLimit, kOuterLimit-0.2, kFALSE);
2656           else
2657             SubtractEtaGap(hist1, kEtaLimit, kOuterLimit, kFALSE);
2658           GetProjections(hist1, &proj2[n][0], &proj2[n][1], count++);
2659           
2660           for (Int_t k=0; k<2; k++)
2661           {
2662             c->cd(1+k*3);
2663             proj2[n][k]->SetLineColor(n+2);
2664             /* TH1* copy = */ proj2[n][k]->DrawCopy("SAME");
2665             
2666             c->cd(2+k*3);
2667             gPad->SetGridx();
2668             gPad->SetGridy();
2669             TH1* ratio = (TH1*) proj2[n][k]->Clone(Form("%s_ratio", proj2[n][k]->GetName()));
2670             ratio->SetStats(0);
2671             ratio->Divide(proj1[k]);
2672             ratio->GetXaxis()->SetRangeUser(-1.5, 1.5);
2673             ratio->Fit("pol0", "0Q", "", -1, 1);
2674             Double_t average = ratio->GetFunction("pol0")->GetParameter(0);
2675 //          Printf("Average is %f", average);
2676             maxAverage[k] = TMath::Max(maxAverage[k], TMath::Abs(average - 1));
2677 //          ratio->Scale(1.0 / average);
2678 //          copy->Scale(1.0 / average);
2679             ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
2680             ratio->DrawCopy((n == 0) ? "" : "SAME");
2681             
2682             if (0)
2683             {
2684               Double_t sum = 0;
2685               Double_t sumCount = 0;
2686               for (Int_t bin = ratio->FindBin(rangeBegin); bin <= ratio->FindBin(rangeEnd); bin++)
2687               {
2688                 sum += TMath::Abs(ratio->GetBinContent(bin) - 1);
2689                 sumCount++;
2690               }
2691               maxDev[k] = TMath::Max(maxDev[k], sum / sumCount);
2692             }
2693
2694             c->cd(3+k*3);
2695             gPad->SetGridx();
2696             gPad->SetGridy();
2697             TH1* diff = (TH1*) proj2[n][k]->Clone(Form("%s_diff", proj2[n][k]->GetName()));
2698             diff->Scale(1.0 / average);
2699             diff->Add(proj1[k], -1);
2700             diff->SetStats(0);
2701             diff->GetXaxis()->SetRangeUser(-1.5, 1.5);
2702             diff->GetYaxis()->SetRangeUser(-0.05, 0.05);
2703             diff->DrawCopy((n == 0) ? "" : "SAME");
2704             
2705             if (0)
2706             {
2707               for (Int_t bin = diff->FindBin(rangeBegin); bin <= diff->FindBin(rangeEnd); bin++)
2708                 maxDev[k] = TMath::Max(maxDev[k], TMath::Abs(diff->GetBinContent(bin)));
2709             }
2710             else
2711             {
2712               Double_t sum = 0;
2713               Double_t sumCount = 0;
2714               for (Int_t bin = diff->FindBin(rangeBegin); bin <= diff->FindBin(rangeEnd); bin++)
2715               {
2716                 sum += TMath::Abs(diff->GetBinContent(bin));
2717                 sumCount++;
2718               }
2719               maxDev[k] = TMath::Max(maxDev[k], sum / sumCount);
2720             }
2721           }
2722           
2723           delete file2;
2724         }
2725         
2726         // normalize to peak strength
2727 /*      for (Int_t k=0; k<2; k++)
2728           maxDev[k] /= proj1[k]->Integral(proj1[k]->GetXaxis()->FindBin(-0.2), proj1[k]->GetXaxis()->FindBin(0.2)) / (proj1[k]->GetXaxis()->FindBin(0.2) - proj1[k]->GetXaxis()->FindBin(-0.2) + 1);*/
2729         
2730         Int_t centralityAxisMapping[] = { 0, 4, 3, 5, 1, 2 };
2731         uncertainty[0]->SetBinContent(ptBin, centralityAxisMapping[histId]+1, maxAverage[0]);
2732         uncertainty[1]->SetBinContent(ptBin, centralityAxisMapping[histId]+1, maxAverage[1]);
2733         uncertainty[2]->SetBinContent(ptBin, centralityAxisMapping[histId]+1, maxDev[0]);
2734         uncertainty[3]->SetBinContent(ptBin, centralityAxisMapping[histId]+1, maxDev[1]);
2735         
2736         c->SaveAs(Form("syst_corr/%s.eps", c->GetName()));
2737       }
2738       
2739       if (draw)
2740         break;
2741     }
2742     if (draw)
2743       break;
2744   }
2745   
2746   TFile::Open("uncertainty.root", "RECREATE");
2747   for (Int_t i=0; i<4; i++)
2748     uncertainty[i]->Write();
2749   gFile->Close();
2750 }
2751
2752 void DrawSystematicsProjections()
2753 {
2754   TH2* uncertainty[4];
2755   
2756   TFile::Open("uncertainty.root");
2757   for (Int_t i=0; i<4; i++)
2758     uncertainty[i] = (TH2*) gFile->Get(Form("uncertainty_%d", i));
2759   
2760   TCanvas* cResult = new TCanvas("cResult", "cResult", 1000, 1000);
2761   cResult->Divide(2, 2);
2762   for (Int_t i=0; i<4; i++)
2763   {
2764     cResult->cd(i+1);
2765     gPad->SetRightMargin(0.15);
2766     uncertainty[i]->SetStats(0);
2767     uncertainty[i]->Draw("colz");
2768   }
2769
2770   TCanvas* cResult2 = new TCanvas("cResult2", "cResult2", 1000, 1000);
2771   cResult2->Divide(2, 2);
2772   for (Int_t i=0; i<4; i++)
2773   {
2774     cResult2->cd(i+1);
2775     gPad->SetRightMargin(0.15);
2776     Int_t color = 1;
2777     for (Int_t j=1; j<=uncertainty[i]->GetNbinsX(); j++)
2778     {
2779       TH1* proj = uncertainty[i]->ProjectionY(Form("p_%d_%d", i, j), j, j);
2780       if (proj->Integral() > 0)
2781       {
2782         proj->SetLineColor(color++);
2783         proj->SetStats(0);
2784         proj->GetYaxis()->SetRangeUser(0, 0.5);
2785         proj->Draw((color == 2) ? "" : "SAME");
2786         
2787         Printf("%d %d: %.3f %.3f", i, color, proj->GetBinContent(1), proj->Integral(2, 6) / 5);
2788       }
2789     }
2790   }
2791 }
2792
2793 void BuildSystematicFiles()
2794 {
2795   Printf("Hope you know what you are doing... 5 seconds to abort...");
2796   gSystem->Sleep(5000);
2797
2798   AnalyzeDeltaPhiEtaGap2D(systBaseFile, "syst_base.root");
2799
2800   kEtaLimit = 0.8;
2801   AnalyzeDeltaPhiEtaGap2D(systBaseFile, "syst_eta08.root");
2802   
2803   kEtaLimit = 1.2;
2804   AnalyzeDeltaPhiEtaGap2D(systBaseFile, "syst_eta12.root");
2805   kEtaLimit = 1.0;
2806   
2807   kOuterLimit = 1.39;
2808   AnalyzeDeltaPhiEtaGap2D(systBaseFile, "syst_outer14.root");
2809   kOuterLimit = 1.59;
2810   
2811   AnalyzeDeltaPhiEtaGap2D(systTrackCuts1, "syst_hybrid.root");
2812   AnalyzeDeltaPhiEtaGap2D(systTrackCuts2, "syst_raa.root");
2813   AnalyzeDeltaPhiEtaGap2D(systVertex, "syst_vertex.root");
2814   AnalyzeDeltaPhiEtaGap2D(systResonances, "syst_resonances.root");
2815   AnalyzeDeltaPhiEtaGap2D(systTTR, "syst_widettr.root");
2816 }
2817
2818 void ExtractSystematicsAll(Int_t offset = 0)
2819 {
2820   gROOT->SetBatch(kTRUE);
2821   
2822   const Int_t NEffects = 9;
2823   
2824 //   const char* defaultFile = "graphs_120425.root";
2825 //   const char* systFiles[] = { "graphs_120425_eta08.root", "graphs_120425_eta12.root", "graphs_120425_outer14.root", "graphs_hybrid_120427.root", "graphs_120430_raa.root", "graphs_120425_vertex.root", "graphs_120502_resonances.root", "graphs_120502_widettr.root", "graphs_120425_wingremoved.root" };
2826   const char* defaultFile = "syst_base.root";
2827   const char* systFiles[] = { "syst_eta08.root", "syst_eta12.root", "syst_outer14.root", "syst_hybrid.root", "syst_raa.root", "syst_vertex.root", "syst_resonances.root", "syst_widettr.root", "graphs_120429_wingremoved.root" };
2828   const char* effectIdString[] = { "1a", "1b", "1c", "2a", "2b", "3", "4", "5", "6" };
2829
2830   Float_t** results[NEffects+3];
2831   
2832   for (Int_t i=0; i<NEffects; i++)
2833   {
2834     results[i] = ExtractSystematics(defaultFile, systFiles[i], offset);
2835   }
2836   
2837   const Int_t NParameters = 6;
2838   const char* names[] = { "$\\sigma_{\\Dphi}$ (fit)", "$\\sigma_{\\Deta}$ (fit)", "$\\sigma_{\\Dphi}$", "$\\sigma_{\\Deta}$", "Kurtosis $\\Dphi$", "Kurtosis $\\Deta$" };
2839   
2840   for (Int_t j=0; j<NParameters; j++)
2841   {
2842     printf("%s \t & ", names[j]);
2843     for (Int_t i=0; i<NEffects; i++)
2844     {
2845       if (j < 4)
2846         printf("$%.1f\\%% \\pm %.1f\\%%$ \t & ", results[i][j][0], results[i][j][1]);
2847       else
2848         printf("$%.2f \\pm %.2f$ \t & ", results[i][j][0], results[i][j][1]);
2849     }
2850     Printf("\\\\");
2851   }
2852
2853   // put together 0-2 (into NEffects)
2854   results[NEffects] = new Float_t*[NParameters];
2855   for (Int_t j=0; j<NParameters; j++)
2856   {
2857     results[NEffects][j] = new Float_t[2];
2858     Float_t mean = 0;
2859     Float_t sigma = 0;
2860     
2861     printf("%s \t:", names[j]);
2862     for (Int_t i=0; i<=2; i++)
2863     {
2864       mean = TMath::Max(mean, TMath::Abs(results[i][j][0]));
2865       printf("%.1f%% ", results[i][j][0]);
2866     }
2867     printf("--> %.1f%% \t", mean);
2868     results[NEffects][j][0] = mean;
2869       
2870     for (Int_t i=0; i<=2; i++)
2871     {
2872       sigma = TMath::Max(sigma, TMath::Abs(results[i][j][1]));
2873       printf("%.1f%% ", results[i][j][1]);
2874     }
2875     Printf("--> %.1f%% \t", sigma);
2876     results[NEffects][j][1] = sigma;
2877   }
2878
2879   // put together 3-4 (into NEffects+1)
2880   results[NEffects+1] = new Float_t*[NParameters];
2881   for (Int_t j=0; j<NParameters; j++)
2882   {
2883     results[NEffects+1][j] = new Float_t[2];
2884     
2885     Float_t mean = 0;
2886     Float_t sigma = 0;
2887     
2888     printf("%s \t:", names[j]);
2889     for (Int_t i=3; i<=4; i++)
2890     {
2891       mean = TMath::Max(mean, TMath::Abs(results[i][j][0]));
2892       printf("%.1f%% ", results[i][j][0]);
2893     }
2894     printf("--> %.1f%% \t", mean);
2895     results[NEffects+1][j][0] = mean;
2896       
2897     for (Int_t i=3; i<=4; i++)
2898     {
2899       sigma = TMath::Max(sigma, TMath::Abs(results[i][j][1]));
2900       printf("%.1f%% ", results[i][j][1]);
2901     }
2902     Printf("--> %.1f%% \t", sigma);
2903     results[NEffects+1][j][1] = sigma;
2904   }
2905   
2906   // combine in quadrature (only mean)
2907   results[NEffects+2] = new Float_t*[NParameters];
2908   for (Int_t j=0; j<NParameters; j++)
2909   {
2910     results[NEffects+2][j] = new Float_t[2];
2911     Float_t mean = 0;
2912   
2913     printf("%s \t:", names[j]);
2914     for (Int_t i=5; i<NEffects+2; i++)
2915     {
2916       mean += results[i][j][0] * results[i][j][0];
2917       printf("%.1f%% ", results[i][j][0]);
2918     }
2919     mean = TMath::Sqrt(mean);
2920     Printf("--> %.1f%% \t", mean);
2921     results[NEffects+2][j][0] = mean;
2922   }
2923     
2924   for (Int_t j=0; j<NParameters; j++)
2925   {
2926     printf("%s \t & ", names[j]);
2927     for (Int_t i=0; i<NEffects+3; i++)
2928     {
2929       if (i == NEffects || i == NEffects+1)
2930         continue;
2931       if (j < 4)
2932         printf("$%.1f\\%%$ \t ", results[i][j][0]);
2933       else
2934         printf("$%.2f$ \t ", results[i][j][0]);
2935       if (i < NEffects+2)
2936         printf("& ");
2937     }
2938     Printf("\\\\");
2939   }
2940   
2941   // generate LaTeX code
2942   for (Int_t i=0; i<NEffects; i++)
2943   {
2944     for (Int_t j=0; j<NParameters; j++)
2945     {
2946       Printf("\\bfig[ht]");
2947       Printf("  \\includegraphics[width=\\linewidth]{syst/%s_%d.eps}", systFiles[i], j);
2948       Printf("  \\caption{Systematic effect (%s) - Parameter %s}", effectIdString[i], names[j]);
2949       Printf("\\efig");
2950     }
2951     Printf("\\clearpage");
2952   }
2953 }
2954
2955 void CompareEtaPhi(const char* fileName, const char* fileNameWingRemoved)
2956 {
2957   Int_t offset = 0;
2958   
2959   TGraphErrors*** graphsWingRemoved = 0;
2960   if (fileNameWingRemoved)
2961   {
2962     ReadGraphs(fileNameWingRemoved);
2963     CalculateRMSSigma();
2964     graphsWingRemoved = graphs;
2965   }
2966   
2967   ReadGraphs(fileName);
2968   
2969   Int_t nHists = 12; //NHists;
2970   skipGraphList[0] = 5;
2971   skipGraphList[1] = 6;
2972   skipGraphList[2] = 10;
2973   
2974   CalculateRMSSigma();
2975   
2976   drawLogo = 1;
2977   
2978   MCLabel = "#sigma_{#Delta#varphi}     #sigma_{#Delta#eta}";
2979
2980   DrawCentrality("rms_centrality", nHists, graphs[1+offset], 0.2, 0.8, "#sigma_{#Delta#varphi} (fit) (rad.) , #sigma_{#Delta#eta} (fit)", graphsWingRemoved[1+offset], 0, graphs[2+offset], graphsWingRemoved[2+offset], 0, 1);
2981 //   DrawCentrality("rms_centrality_nosyst", nHists, graphs[1], 0, 0.8, "#sigma_{#varphi} (fit) (rad.) / #sigma_{#eta}", 0, 0,  graphs[2], 0);
2982 }
2983
2984 void DrawExamples(const char* histFileName, const char* graphFileName, Int_t i = 0, Int_t j = 1, Bool_t drawFunc = kTRUE)
2985 {
2986   Float_t etaLimit = kEtaLimit;
2987   Float_t outerLimit = kOuterLimit;
2988   Float_t projectLimit = 0.8;
2989
2990   ReadGraphs(graphFileName);
2991
2992   TFile::Open(histFileName);
2993   
2994   Int_t exColors[] = { 1, 2, 4, 3, 5, 6 };
2995   
2996   Int_t graphID = i * (6 - 1) + j - 1;
2997
2998   TCanvas* c = new TCanvas(Form("c_%d", i), Form("c_%d", i), 1200, 600);
2999   c->Divide(2, 1);
3000   Int_t nHists = 3;
3001   for (Int_t histId = 0; histId < nHists; histId++)
3002   {
3003     TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
3004     if (!hist)
3005       continue;
3006
3007     if (graphs[0][graphID]->GetN() < histId)
3008     {
3009       Printf("ERROR: Pos in graph not found: %d %d", i, histId);
3010       continue;
3011     }
3012     
3013     SubtractEtaGap(hist, etaLimit, outerLimit, kFALSE, kFALSE);
3014   
3015     c->cd(1);
3016     TH1* proj = hist->ProjectionX(Form("%s_proj1", hist->GetName()), hist->GetYaxis()->FindBin(-projectLimit+0.01), hist->GetYaxis()->FindBin(projectLimit-0.01));
3017     proj->SetLineColor(exColors[histId]);
3018     proj->SetStats(0);
3019     proj->GetXaxis()->SetRangeUser(-1.49, 1.49);
3020     proj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 1.2, proj->GetMaximum() * 1.5);
3021     proj->GetYaxis()->SetTitle(kProjYieldTitlePhi);
3022     TString label(proj->GetTitle());
3023     TObjArray* objArray = label.Tokenize("-");
3024     proj->SetTitle(Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()));
3025     proj->Draw((histId == 0) ? "" : "SAME");
3026   
3027     if (drawFunc)
3028     {
3029       // integral over y
3030       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());
3031       func->SetParameter(0, 0);
3032       for (Int_t k=0; k<6; k++)
3033         func->SetParameter(k+1, graphs[k][graphID]->GetY()[histId]);
3034       // scale by bin width (to compare with projection)
3035       func->SetParameter(1, func->GetParameter(1) / hist->GetYaxis()->GetBinWidth(1));
3036       func->SetLineColor(exColors[histId]);
3037       func->SetLineWidth(1);
3038       func->DrawCopy("SAME");
3039       
3040       // draw contributions
3041       Float_t scale = func->GetParameter(1);
3042       Float_t weighting = func->GetParameter(4);
3043       func->SetParameter(1, scale * weighting);
3044       func->SetParameter(4, 1);
3045       func->SetLineStyle(2);
3046       func->DrawCopy("SAME");
3047
3048       func->SetParameter(1, scale * (1.0-weighting));
3049       func->SetParameter(4, 0);
3050       func->SetLineStyle(3);
3051       func->DrawCopy("SAME");
3052     }
3053
3054     DrawLatex(0.15, 0.8 - 0.05 * histId, exColors[histId], labels[histId]);
3055
3056     c->cd(2);
3057     proj = hist->ProjectionY(Form("%s_proj2b", hist->GetName()), hist->GetXaxis()->FindBin(-projectLimit+0.01), hist->GetXaxis()->FindBin(projectLimit-0.01));
3058     proj->SetLineColor(exColors[histId]);
3059     proj->GetXaxis()->SetRangeUser(-1.49, 1.49);
3060     proj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 1.2, proj->GetMaximum() * 2);
3061     proj->GetYaxis()->SetTitle(kProjYieldTitleEta);
3062     proj->SetStats(0);
3063     proj->SetTitle(Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()));
3064     proj->Draw((histId == 0) ? "" : "SAME");
3065     
3066     if (0)
3067     {
3068       proj = hist->ProjectionY(Form("%s_proj2", hist->GetName()), hist->GetXaxis()->FindBin(-0.5 * TMath::Pi()), hist->GetXaxis()->FindBin(0.5 * TMath::Pi()));
3069       proj->SetLineColor(exColors[histId]);
3070       proj->GetXaxis()->SetRangeUser(-1.2, 1.2);
3071       proj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 1.2, proj->GetMaximum() * 2);
3072       proj->SetLineStyle(2);
3073       proj->Draw("SAME");
3074     }
3075
3076     if (drawFunc)
3077     {
3078       // integral over x
3079       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);
3080       func->SetParameter(0, 0);
3081       for (Int_t k=0; k<6; k++)
3082         func->SetParameter(k+1, graphs[k][graphID]->GetY()[histId]);
3083       // scale by bin width (to compare with projection)
3084       func->SetParameter(1, func->GetParameter(1) / hist->GetXaxis()->GetBinWidth(1));
3085       func->SetLineColor(exColors[histId]);
3086       func->SetLineWidth(1);
3087       func->DrawCopy("SAME");
3088
3089       // draw contributions
3090       Float_t scale = func->GetParameter(1);
3091       Float_t weighting = func->GetParameter(4);
3092       func->SetParameter(1, scale * weighting);
3093       func->SetParameter(4, 1);
3094       func->SetLineStyle(2);
3095       func->DrawCopy("SAME");
3096
3097       func->SetParameter(1, scale * (1.0-weighting));
3098       func->SetParameter(4, 0);
3099       func->SetLineStyle(3);
3100       func->DrawCopy("SAME");   
3101     }
3102
3103     DrawLatex(0.15, 0.8 - 0.05 * histId, exColors[histId], labels[histId]);
3104   }
3105 }
3106
3107 TH1* GetSystUnc(TH1* hist, Int_t i, Int_t j, Int_t histId, Int_t etaPhi)
3108 {
3109 //   Float_t uncertainty = (histId == 0) ? 0.1 : 0.08;
3110
3111   // track cuts
3112   //   2 2: 0.048 0.020
3113   //   2 3: 0.078 0.031
3114   //   2 4: 0.023 0.008
3115   //   2 5: 0.159 0.056
3116   //   2 6: 0.043 0.015
3117   //   2 7: 0.011 0.006
3118   //   3 2: 0.050 0.022
3119   //   3 3: 0.089 0.033
3120   //   3 4: 0.020 0.008
3121   //   3 5: 0.183 0.056
3122   //   3 6: 0.042 0.013
3123   //   3 7: 0.010 0.005
3124   
3125   // others
3126   //   2 2: 0.023 0.007
3127   //   2 3: 0.024 0.015
3128   //   2 4: 0.008 0.005
3129   //   2 5: 0.075 0.024
3130   //   2 6: 0.019 0.010
3131   //   2 7: 0.007 0.004
3132   //   3 2: 0.029 0.009
3133   //   3 3: 0.032 0.012
3134   //   3 4: 0.005 0.003
3135   //   3 5: 0.086 0.013
3136   //   3 6: 0.015 0.006
3137   //   3 7: 0.008 0.003
3138
3139   Float_t uncertainty = 1;
3140   if (i == 0 && j == 1)
3141     uncertainty = 0.021;
3142   if (i == 1 && j == 1)
3143     uncertainty = 0.032;
3144   if (i == 1 && j == 2)
3145     uncertainty = 0.008;
3146   if (i == 2 && j == 1)
3147     uncertainty = 0.056;
3148   if (i == 2 && j == 2)
3149     uncertainty = 0.014;
3150   if (i == 2 && j == 3)
3151     uncertainty = 0.006;
3152   
3153   if (histId == 0)
3154     uncertainty *= 2;
3155   
3156   TH1* systUnc = (TH1*) hist->Clone(Form("%s_syst", hist->GetName()));
3157   for (Int_t n=1; n<=systUnc->GetNbinsX(); n++)
3158     systUnc->SetBinError(n, uncertainty);
3159   
3160   systUnc->SetFillColor(hist->GetLineColor());
3161   systUnc->SetFillStyle((etaPhi == 0) ? 3004 : 3005);
3162   systUnc->SetMarkerStyle(0);
3163   systUnc->SetLineColor(0);
3164   
3165   return systUnc;
3166 }
3167
3168 void Draw2DExamples(const char* histFileName)
3169 {
3170   Int_t is[] = { 0, 0, 1, 1, 2, 2, 3 };
3171   Int_t js[] = { 1, 2, 2, 3, 4, 5, 6 };
3172   Int_t cs[] = { 0, 3, 5, 1 };
3173   Float_t outerLimit = kOuterLimit;
3174
3175   TFile::Open(histFileName);
3176
3177   TCanvas* canvas = new TCanvas("2d", "2d", 1800, 600);
3178   canvas->Divide(7, 4);
3179   Int_t cID = 1;
3180
3181   for (Int_t c=0; c<4; c++)
3182   {
3183     for (Int_t i=0; i<7; i++)
3184     {
3185       canvas->cd(cID++);
3186
3187       TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", is[i], js[i], cs[c]));
3188       if (!hist)
3189         continue;
3190
3191       // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
3192       hist->Scale(1.0 / hist->GetYaxis()->GetBinWidth(1));
3193       hist->GetZaxis()->SetTitle(kCorrFuncTitle);
3194       hist->GetZaxis()->SetTitleOffset(1.9);
3195       
3196       hist->Rebin2D(2, 2);
3197       hist->Scale(0.25);
3198   
3199       TString label(hist->GetTitle());
3200       label.ReplaceAll(".00", " GeV/c");
3201       label.ReplaceAll(".0", " GeV/c");
3202       TObjArray* objArray = label.Tokenize("-");
3203       TPaveText* paveText = new TPaveText(0.52, 0.72, 0.95, 0.95, "BRNDC");
3204       paveText->SetTextSize(0.035);
3205       paveText->SetFillColor(0);
3206       paveText->SetShadowColor(0);
3207       paveText->AddText(objArray->At(0)->GetName());
3208       paveText->AddText(objArray->At(1)->GetName());
3209       if (objArray->GetEntries() == 4)
3210         paveText->AddText(Form("Pb-Pb 2.76 TeV %s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()));
3211       else
3212         paveText->AddText(Form("%s 2.76 TeV", objArray->At(2)->GetName()));
3213       paveText->AddText("|#eta| < 0.9");
3214       
3215       gPad->SetLeftMargin(0.15);
3216       hist->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
3217       hist->GetXaxis()->SetTitleOffset(1.5);
3218       hist->GetYaxis()->SetTitleOffset(1.7);
3219       hist->SetStats(0);
3220       hist->SetTitle("a) Correlation");
3221       TH2* clone = (TH2*) hist->Clone(Form("%s_clone", hist->GetName()));
3222     //   clone->GetXaxis()->SetRangeUser(-TMath::Pi() / 2, TMath::Pi() / 2);
3223       clone->Draw("SURF1");
3224       paveText->Draw();
3225     }
3226   }
3227 }
3228
3229 Bool_t disableUncertainties = kTRUE;
3230
3231 void DrawExample(const char* histFileName, Int_t i, Int_t j, Int_t histId, TH1** projPhi = 0, TH1** projEta = 0, TH1** projPhiSyst = 0 , TH1** projEtaSyst = 0, Bool_t Ratio = kFALSE)
3232 {
3233   Float_t etaLimit = kEtaLimit;
3234   Float_t outerLimit = kOuterLimit;
3235   Float_t projectLimit = 0.8;
3236   Int_t maxAssocPt = 6;
3237   Int_t graphID = i * (maxAssocPt - 1) + j - 1;
3238
3239   TFile::Open(histFileName);
3240   
3241   TCanvas* c = new TCanvas(Form("ex_%d_%d_%d", i, j, histId), Form("ex_%d_%d_%d", i, j, histId), 1800, 600);
3242   c->Divide(3, 1);
3243
3244   TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
3245   if (!hist)
3246     return;
3247   // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
3248   hist->Scale(1.0 / hist->GetYaxis()->GetBinWidth(1));