ea99de74c9513f53b32d183ea9ec511453abf713
[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,moment3phi,moment3eta,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,moment3phi,moment3eta,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++)
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 Float_t kEtaLimit = 1.0;
510 Float_t kOuterLimit = 1.59;
511
512
513 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)
514 {
515   Bool_t success = kTRUE;
516   Float_t etaLimit = kEtaLimit;
517   Float_t outerLimit = kOuterLimit;
518
519   if (graphID >= 15)
520   {
521     etaLimit = 0.5;
522     outerLimit = 0.99;
523   }
524   Float_t sigmaFitLimit = 0.1 - limits * 0.05;
525   Float_t etaFitUpperLimit = 0.8;
526   Float_t initSigma = 0.6;
527   if (histId == 2) // pp
528   {
529     etaFitUpperLimit = 0.6;
530     initSigma = 0.4;
531   }
532
533   hist->GetYaxis()->SetRangeUser(-outerLimit+0.01, outerLimit-0.01);
534   hist->GetXaxis()->SetRangeUser(-TMath::Pi() / 2 + 0.01, TMath::Pi() * 0.5 - 0.01);
535   
536   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);
537 //   Printf("%f", mean);
538
539   // sums
540   TGraphErrors* sumSummary = new TGraphErrors;
541   TGraphErrors* phiWidthSummary = new TGraphErrors;
542   TGraphErrors* etaWidthSummary = new TGraphErrors;
543   
544   canvas->cd(canvasPos++);
545   hist->SetStats(0);
546   hist->DrawCopy("SURF1");
547   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), hist->Integral());
548   
549   Float_t min = hist->GetMinimum();
550   Float_t max = hist->GetMaximum();
551
552   Int_t bins = hist->GetNbinsX() / 2 / 2;
553     
554   TF2* func = new TF2("func", DeltaPhiWidth2DFitFunction, -TMath::Pi() / 2, TMath::Pi() * 0.5, -outerLimit, outerLimit, bins+6);
555   func->SetParameters(hist->GetBinContent(hist->GetXaxis()->FindBin(0.0), hist->GetYaxis()->FindBin(0.0)) - mean, 0.3, 0.3, 0.25, initSigma, initSigma);
556   for (Int_t i=6; i<bins+6; i++)
557     func->SetParameter(i, mean);
558  
559  
560   if (1)
561   {
562     // STEP 1: fit only flow using one delta eta side
563     for (Int_t i=0; i<6; i++)
564       func->FixParameter(i, func->GetParameter(i));
565     hist->GetYaxis()->SetRangeUser(etaLimit+0.01, outerLimit-0.01);
566     Int_t fitResult = hist->Fit(func, "0", "");
567     Printf("Fit result: %d; Chi2/ndf: %f/%d", fitResult, func->GetChisquare(), func->GetNDF());
568     if (fitResult != 0)
569       success = kFALSE;
570
571     // STEP2 : fit only Gaussian in central region
572     for (Int_t i=0; i<6; i++)
573       func->ReleaseParameter(i);
574     //   func->SetParameters(1, 0.3, 0.3, 0.25, initSigma, initSigma);
575     func->SetParLimits(0, 0, 10);
576     func->SetParLimits(1, sigmaFitLimit, 0.6);
577     func->SetParLimits(2, sigmaFitLimit, etaFitUpperLimit + (((graphID == 5 && histId == 4) || (graphID == 0 && histId == 0)) ? 0.1 : 0));
578     func->SetParLimits(3, 0.1, 0.9);
579     func->SetParLimits(4, sigmaFitLimit, 0.6);
580     func->SetParLimits(5, sigmaFitLimit, etaFitUpperLimit + (((graphID == 5 && histId == 4) || (graphID == 0 && histId == 0)) ? 0.1 : 0));
581     for (Int_t i=6; i<bins+6; i++)
582       func->FixParameter(i, func->GetParameter(i));
583     hist->GetYaxis()->SetRangeUser(-etaLimit+0.01, etaLimit-0.01);
584     fitResult = hist->Fit(func, "0", "");
585     Printf("Fit result: %d; Chi2/ndf: %f/%d", fitResult, func->GetChisquare(), func->GetNDF());
586     if (fitResult != 0)
587       success = kFALSE;
588
589     // STEP3: fit everything, with limits
590     for (Int_t i=6; i<bins+6; i++)
591     {
592       func->ReleaseParameter(i);
593       func->SetParLimits(i, func->GetParameter(i) * 0.8, func->GetParameter(i) * 1.2);
594     }
595     hist->GetYaxis()->SetRangeUser(-outerLimit+0.01, outerLimit-0.01);
596
597     // STEP4: fit everything, without limits
598     for (Int_t i=6; i<bins+6; i++)
599       func->SetParLimits(i, 0, 0);
600     fitResult = hist->Fit(func, "0", "");
601     Printf("Fit result: %d; Chi2/ndf: %f/%d", fitResult, func->GetChisquare(), func->GetNDF());
602     if (fitResult != 0)
603       success = kFALSE;
604   }
605   Int_t fitResult = hist->Fit(func, "0", "");
606   Printf("Fit result: %d; Chi2/ndf: %f/%d", fitResult, func->GetChisquare(), func->GetNDF());
607   if (fitResult != 0)
608     success = kFALSE;
609
610   Printf("Trying 1 Gaussian...");
611   TF2* func_clone = new TF2("func_clone", DeltaPhiWidth2DFitFunction, -TMath::Pi() / 2, TMath::Pi() * 0.5, -outerLimit, outerLimit, bins+6);
612   for (Int_t i=0; i<bins+6; i++)
613   {
614     func_clone->SetParameter(i, func->GetParameter(i));
615     Double_t parmin, parmax;
616     func->GetParLimits(i, parmin, parmax);
617     func_clone->SetParLimits(i, parmin, parmax);
618   }
619   func_clone->SetParLimits(3, 1, 1);
620   func_clone->FixParameter(3, 1);
621   func_clone->FixParameter(4, sigmaFitLimit);
622   func_clone->FixParameter(5, sigmaFitLimit);
623   fitResult = hist->Fit(func_clone, "0R", "");
624   Printf("Fit result: %d", fitResult);
625   
626   // if both parameters are within 1%, refit with 1 Gaussian only
627   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)
628   {
629     Printf("Parameters within 1%%. Using result from 1 Gaussian...");
630     
631     func = func_clone;
632     
633     if (fitResult != 0)
634       success = kFALSE;
635   }
636   else if (func_clone->GetChisquare() * 0.99 < func->GetChisquare())
637   {
638     Printf("1 Gaussian fit has a at maximum 1%% (%f %f) larger chi2. Using results from 1 Gaussian...", func_clone->GetChisquare(), func->GetChisquare());
639     
640     func = func_clone;
641     
642     if (fitResult != 0)
643       success = kFALSE;
644   } 
645   
646   Int_t first = 1;
647   Int_t second = 4;
648   if (func->GetParameter(1) < func->GetParameter(4))
649   {
650     first = 4;
651     second = 1;
652   }
653   //dphi
654   AddPoint(graphs[1+16][graphID], x, TMath::Abs(func->GetParameter(first)), xE, func->GetParError(first));
655   AddPoint(graphs[4+16][graphID], x, TMath::Abs(func->GetParameter(second)), xE, func->GetParError(second));
656
657   //deta
658   AddPoint(graphs[2+16][graphID], x, TMath::Abs(func->GetParameter(first+1)), xE, func->GetParError(first+1));
659   AddPoint(graphs[5+16][graphID], x, TMath::Abs(func->GetParameter(second+1)), xE, func->GetParError(second+1));
660   
661   // norm
662   AddPoint(graphs[0+16][graphID], x, func->GetParameter(0), xE, func->GetParError(0));
663   if (first < second)
664     AddPoint(graphs[3+16][graphID], x, func->GetParameter(3), xE, func->GetParError(3));
665   else
666     AddPoint(graphs[3+16][graphID], x, 1.0 - func->GetParameter(3), xE, func->GetParError(3));    
667   
668   //   hist->Fit(func, "0RI", "");
669
670   //   width[0]->SetPoint(width[0]->GetN(), x, TMath::Abs(func->GetParameter(1)));
671   //   width[0]->SetPointError(width[0]->GetN()-1, 0, func->GetParError(1));
672   //     
673   //   width[1]->SetPoint(width[1]->GetN(), x, TMath::Abs(func->GetParameter(2)));
674   //   width[1]->SetPointError(width[1]->GetN()-1, 0, func->GetParError(2));
675
676   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), func->GetParameter(1), 0, func->GetParError(1));
677   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), func->GetParameter(4), 0, func->GetParError(4));
678   
679   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), func->GetParameter(2), 0, func->GetParError(2));
680   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), func->GetParameter(5), 0, func->GetParError(5));
681
682   canvas->cd(canvasPos++);
683   TH2* funcHist = (TH2*) hist->Clone("funcHist");
684   funcHist->Reset();
685   funcHist->Add(func);
686   funcHist->SetMinimum(min);
687   funcHist->SetMaximum(max);
688   funcHist->Draw("SURF1");
689   
690   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), funcHist->Integral());
691   
692   canvas->cd(canvasPos++);
693   TH2* residuals = (TH2*) hist->Clone("residuals");
694   residuals->Add(func, -1);
695   residuals->SetMinimum(-(max - min) / 2);
696   residuals->SetMaximum((max - min) / 2);
697   residuals->Draw("SURF1");
698
699   Float_t chi2 = 0;
700   Int_t ndf = 0;
701   for (Int_t i=hist->GetXaxis()->FindBin(-0.8); i<=hist->GetXaxis()->FindBin(0.8); i++)
702     for (Int_t j=hist->GetYaxis()->FindBin(-etaLimit); j<=hist->GetYaxis()->FindBin(etaLimit); j++)
703     {
704       if (residuals->GetBinError(i, j) > 0)
705       {
706         chi2 += TMath::Power(residuals->GetBinContent(i, j) / residuals->GetBinError(i, j), 2);
707         ndf++;
708       }
709     }
710   ndf -= func->GetNumberFreeParameters();
711   
712   if (func->GetNDF() > 0)
713   {
714     printf("#chi^{2}/ndf = %.1f/%d = %.1f  ", func->GetChisquare(), func->GetNDF(), func->GetChisquare() / func->GetNDF());
715     DrawLatex(0.2, yPosChi2, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", func->GetChisquare(), func->GetNDF(), func->GetChisquare() / func->GetNDF()));
716     AddPoint(graphs[6+16][graphID], x, func->GetChisquare() / func->GetNDF(), xE, 0);
717   }
718   if (ndf)
719   {
720     Printf("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf);
721     DrawLatex(0.2, yPosChi2 - 0.05, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf));
722     AddPoint(graphs[7+16][graphID], x, chi2 / ndf, xE, 0);
723   }
724   
725   // draw gaussian only
726   TF2* funcClone = new TF2("funcClone", DeltaPhiWidth2DFitFunction, -TMath::Pi() / 2, TMath::Pi() * 0.5, -outerLimit, outerLimit, bins+6);
727   for (Int_t i=0; i<6; i++)
728     funcClone->SetParameter(i, func->GetParameter(i));
729   for (Int_t i=6; i<bins+6; i++)
730     funcClone->SetParameter(i, 0);
731 //   funcClone->Print();
732   canvas->cd(canvasPos++);
733   funcHist = (TH2*) hist->Clone("funcHistb");
734   funcHist->Reset();
735   funcHist->Add(funcClone);
736   funcHist->SetMinimum(-(max - min) / 2);
737   funcHist->SetMaximum((max - min) / 2);
738   funcHist->Draw("SURF1");
739   
740   // eta gap subtraction
741   canvas->cd(canvasPos++);
742   func->SetParameter(0, 0);
743   TH2* subtractFlow = (TH2*) hist->Clone("subtractFlow");
744   subtractFlow->Add(func, -1);
745   subtractFlow->SetMinimum(-(max - min) / 2);
746   subtractFlow->SetMaximum((max - min) / 2);
747   subtractFlow->DrawCopy("SURF1");
748   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), subtractFlow->Integral());
749 /*
750   Float_t etaLimitSubtract = 0;
751   Float_t outerLimitSubtract = 0;
752   if (graphID >= 15)
753   {
754     etaLimitSubtract = 0.5;
755     outerLimitSubtract = 0.99;
756   }
757   else 
758   {
759     etaLimitSubtract = etaLimit;
760     outerLimitSubtract = outerLimit;
761   }
762 */
763
764   Float_t momentFitLimit = 0.8 - 1e-4;
765   if (graphID >= 15)
766     momentFitLimit = 0.4 - 1e-4;
767
768   TH1* projx3 = hist->ProjectionX(Form("%s_projx3", hist->GetName()), hist->GetYaxis()->FindBin(-etaLimit+0.01), hist->GetYaxis()->FindBin(etaLimit-0.01));
769   Float_t nBins = hist->GetYaxis()->FindBin(etaLimit-0.01) - hist->GetYaxis()->FindBin(-etaLimit+0.01)+1;
770   projx3->Scale(1.0/nBins);
771   
772   TH1* projx3SubtractNegative = hist->ProjectionX(Form("%s_projx3SubtractPositive", hist->GetName()), hist->GetYaxis()->FindBin(-outerLimit+0.01), hist->GetYaxis()->FindBin(-etaLimit-0.01));
773   nBins = hist->GetYaxis()->FindBin(-etaLimit-0.01) - hist->GetYaxis()->FindBin(-outerLimit+0.01)+1;
774   projx3SubtractNegative->Scale(1.0/nBins);
775
776   TH1* projx3SubtractPositive = hist->ProjectionX(Form("%s_projx3SubtractNegative ", hist->GetName()), hist->GetYaxis()->FindBin(etaLimit+0.01), hist->GetYaxis()->FindBin(outerLimit-0.01));
777   nBins = hist->GetYaxis()->FindBin(outerLimit-0.01) - hist->GetYaxis()->FindBin(etaLimit+0.01)+1;
778   projx3SubtractPositive->Scale(1.0/nBins);
779
780   TH1* projy3 = hist->ProjectionY(Form("%s_projy3", hist->GetName()), hist->GetXaxis()->FindBin(-momentFitLimit+0.01), hist->GetXaxis()->FindBin(momentFitLimit-0.01));
781   nBins = hist->GetXaxis()->FindBin(momentFitLimit-0.01) - hist->GetXaxis()->FindBin(-momentFitLimit+0.01)+1;
782   projy3->Scale(1.0/nBins);
783
784   SubtractEtaGap1D(projx3, projx3SubtractPositive, projx3SubtractNegative, projy3, etaLimit, outerLimit);
785   projy3->GetXaxis()->SetRangeUser(-momentFitLimit+0.01,momentFitLimit-0.01);
786   projx3->GetXaxis()->SetRangeUser(-momentFitLimit+0.01, momentFitLimit-0.01);
787   SubtractEtaGap(hist, etaLimit, outerLimit, kFALSE);
788 //  return 0;
789   canvas->cd(canvasPos++);
790   hist->SetMinimum(-(max - min) / 2);
791   hist->SetMaximum((max - min) / 2);
792   hist->DrawCopy("SURF1");
793   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), hist->Integral());
794
795   if (!quick)
796   {
797     canvas->cd(canvasPos++);
798     TH2* difference = (TH2*) hist->Clone("difference");
799     difference->Add(subtractFlow, -1);
800     difference->SetMinimum(-(max - min) / 2);
801     difference->SetMaximum((max - min) / 2);
802     difference->DrawCopy("SURF1");
803
804     canvas->cd(canvasPos++);
805     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);
806     func2->SetParameters(0, 1, 0.3, 0.3);
807     func2->SetParLimits(1, 0, 10);
808     func2->SetParLimits(2, sigmaFitLimit, 1);
809     func2->SetParLimits(3, sigmaFitLimit, 1);
810     func2->FixParameter(0, 0);
811     
812     hist->Fit(func2, "0R", "");
813     //   hist->Fit(func2, "I0R", "");
814     
815     AddPoint(phiWidthSummary, phiWidthSummary->GetN(), func2->GetParameter(2), 0, func2->GetParError(2));
816     AddPoint(etaWidthSummary, etaWidthSummary->GetN(), func2->GetParameter(3), 0, func2->GetParError(3));
817   
818     TH2* funcHist2 = (TH2*) hist->Clone("funcHist2");
819     funcHist2->Reset();
820     funcHist2->Add(func2);
821     funcHist2->SetMinimum(-(max - min) / 2);
822     funcHist2->SetMaximum((max - min) / 2);
823     funcHist2->Draw("SURF1");
824     sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), funcHist2->Integral());
825     
826     if (func2->GetNDF() > 0)
827     {
828       Printf("#chi^{2}/ndf = %.1f/%d = %.1f", func2->GetChisquare(), func2->GetNDF(), func2->GetChisquare() / func2->GetNDF());
829       DrawLatex(0.2, yPosChi2, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", func2->GetChisquare(), func2->GetNDF(), func2->GetChisquare() / func2->GetNDF()));
830     }
831
832     canvas->cd(canvasPos++);
833     TH2* residuals2 = (TH2*) hist->Clone("residuals");
834     residuals2->Add(funcHist2, -1);
835     residuals2->SetMinimum(-(max - min) / 2);
836     residuals2->SetMaximum((max - min) / 2);
837     residuals2->Draw("SURF1");
838   }
839
840
841   TH1* projx2 = hist->ProjectionX(Form("%s_projx2", hist->GetName()), hist->GetYaxis()->FindBin(-etaLimit+0.01), hist->GetYaxis()->FindBin(etaLimit-0.01));
842   projx2->GetXaxis()->SetRangeUser(-momentFitLimit, momentFitLimit);
843   nBins = hist->GetYaxis()->FindBin(etaLimit-0.01) - hist->GetYaxis()->FindBin(-etaLimit+0.01)+1;
844   projx2->Scale(1.0/nBins);
845   
846   TH1* projy2 = hist->ProjectionY(Form("%s_projy2", hist->GetName()), hist->GetXaxis()->FindBin(-momentFitLimit+0.01), hist->GetXaxis()->FindBin(momentFitLimit-0.01));
847   projy2->GetXaxis()->SetRangeUser(-momentFitLimit, momentFitLimit);
848   nBins = hist->GetXaxis()->FindBin(momentFitLimit-0.01) - hist->GetXaxis()->FindBin(-momentFitLimit+0.01)+1;
849   projy2->Scale(1.0/nBins);
850   
851   CalculateMomentsKurtosis(momentFitLimit, projx2, 8, graphID, x, xE);
852   CalculateMomentsKurtosis(momentFitLimit, projy2, 9, graphID, x, xE);
853
854 //     return success;
855   
856   TH1* projx1 = subtractFlow->ProjectionX(Form("%s_projx1", hist->GetName()), hist->GetYaxis()->FindBin(-0.79), hist->GetYaxis()->FindBin(0.79));
857   projx1->GetXaxis()->SetRangeUser(-momentFitLimit, momentFitLimit);
858
859   TH1* projy1 = subtractFlow->ProjectionY(Form("%s_projy1", hist->GetName()), hist->GetXaxis()->FindBin(-0.79), hist->GetXaxis()->FindBin(0.79));
860   projy1->GetXaxis()->SetRangeUser(-momentFitLimit, momentFitLimit);
861
862   CalculateMomentsKurtosis(momentFitLimit, projx1, 8+16, graphID, x, xE);
863   CalculateMomentsKurtosis(momentFitLimit, projy1, 9+16, graphID, x, xE);
864
865 //   TF1* twoGauss = new TF1("twoGauss", "gaus(0)+gaus(3)", -2, 2);
866 //   twoGauss->SetParameters(1, 0, 0.3, 1, 0, 0.6);
867 //   twoGauss->FixParameter(1, 0);
868 //   twoGauss->FixParameter(4, 0);
869 //   twoGauss->SetLineColor(4);
870 //   projx1->Fit("twoGauss", "I+", "SAME");
871   
872   canvas->cd(canvasPos++);
873   projx1->Draw();
874   projx1->Fit("gaus", "I");
875   projx1->GetFunction("gaus")->SetLineColor(1);
876   projx1->GetYaxis()->SetRangeUser(0, projx1->GetMaximum() * 1.05);
877
878   projx2->SetLineColor(2);
879   projx2->Draw("SAME");
880   projx2->Fit("gaus", "I+0", "SAME");
881   projx2->GetFunction("gaus")->SetLineColor(2);
882
883   canvas->cd(canvasPos--);
884   projy1->Draw();
885   projy1->Fit("gaus", "I");
886   projy1->GetFunction("gaus")->SetLineColor(1);
887   projy1->GetYaxis()->SetRangeUser(0, projy1->GetMaximum() * 1.05);
888   
889   projy2->SetLineColor(2);
890   projy2->Draw("SAME");
891   projy2->Fit("gaus", "I+0", "SAME");
892   projy2->GetFunction("gaus")->SetLineColor(2);
893   
894   // 1d fit (eta gap subtraction)
895   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), projx2->GetFunction("gaus")->GetParameter(2), 0, projx2->GetFunction("gaus")->GetParError(2));
896   
897   // 1d fit (lots of params)
898   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), projx1->GetFunction("gaus")->GetParameter(2), 0, projx1->GetFunction("gaus")->GetParError(2));
899
900   // 1d fit (eta gap subtraction)
901   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), projy2->GetFunction("gaus")->GetParameter(2), 0, projy2->GetFunction("gaus")->GetParError(2));
902   
903   // 1d fit (lots of params)
904   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), projy1->GetFunction("gaus")->GetParameter(2), 0, projy1->GetFunction("gaus")->GetParError(2));
905
906   // set errors large for bins potentially affected by two-track effects
907   if (0 && twoTrack)
908   {
909     Printf("NOTE : Skipping bins at (0, 0)");
910 //     for (Int_t binx = hist->GetXaxis()->FindBin(-0.25); binx <= hist->GetXaxis()->FindBin(0.25); binx++)
911     for (Int_t binx = hist->GetXaxis()->FindBin(-0.01); binx <= hist->GetXaxis()->FindBin(0.01); binx++)
912       for (Int_t biny = hist->GetYaxis()->FindBin(-0.01); biny <= hist->GetYaxis()->FindBin(0.01); biny++)
913       {
914 //      hist->SetBinContent(binx, biny,  0);
915         hist->SetBinError(binx, biny,  1e5);
916       }
917   }
918
919   Float_t etaFitLimit = outerLimit;
920
921   Bool_t oneGaussian = kTRUE;
922
923   if (!oneGaussian) cout << endl << "Fitting two 1D Gaussians" << endl;
924   else cout << endl << "Fitting one 1D Gaussians" << endl;
925
926   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);
927   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);
928
929   func4phi->SetParLimits(1, 0, 10);
930   func4phi->SetParLimits(3, sigmaFitLimit, 0.7);
931   func4phi->FixParameter(0, 0);
932   if (!oneGaussian)
933   {
934     func4phi->SetParLimits(4, 0, 1);
935     func4phi->SetParLimits(2, sigmaFitLimit, 0.7);
936   }
937   else
938   {
939     func4phi->FixParameter(4, 0);
940     func4phi->FixParameter(2, 0);
941   }
942
943   func4eta->SetParLimits(1, 0, 10);
944   func4eta->SetParLimits(3, sigmaFitLimit, etaFitUpperLimit);
945   func4eta->FixParameter(0, 0);
946   if (!oneGaussian)
947   {
948     func4eta->SetParLimits(4, 0, 1);
949     func4eta->SetParLimits(2, sigmaFitLimit, etaFitUpperLimit);
950   }
951   else
952   {
953     func4eta->FixParameter(4, 0);
954     func4eta->FixParameter(2, 0);
955   }
956   canvas->cd(canvasPos++);
957   func4phi->SetLineColor(4);
958   TFitResultPtr fitResultPtr = projx3->Fit(func4phi, "SIR", "SAME");
959   TMatrixDSym cov = fitResultPtr->GetCovarianceMatrix();
960
961   projx3->Draw("SAME");
962   projx3->SetLineColor(4);
963
964   if (func4phi->GetNDF() > 0)
965   {
966     printf("#chi^{2}/ndf = %.1f/%d = %.1f  ", func4phi->GetChisquare(), func4phi->GetNDF(), func4phi->GetChisquare() / func4phi->GetNDF());
967     AddPoint(graphs[43][graphID], x, func4phi->GetChisquare() / func4phi->GetNDF(), xE, 0);
968   }
969   chi2 = 0;
970   ndf = 0;
971
972   for (Int_t i=projx3->FindBin(-momentFitLimit); i<=projx3->FindBin(momentFitLimit);i++)
973   {
974 //    Float_t temp = chi2;
975     if (projx3->GetBinError(i) > 0)
976     {
977       chi2 += TMath::Power((projx3->GetBinContent(i)-func4phi->Eval(projx3->GetBinCenter(i)))/projx3->GetBinError(i),2);
978       //    cerr << i << "\t" << projx3->GetBinCenter(i) << "\t" << chi2-temp << endl;
979       ndf++;
980     }
981   }
982   if (!oneGaussian) ndf = ndf - 4;
983   else ndf = ndf - 2;
984   printf("Calculated #chi^{2}/ndf = %.1f/%d = %.1f  ", chi2, ndf, chi2/ndf);
985
986   first = 2;
987   second = 3;
988   if (func4phi->GetParameter(2) < func4phi->GetParameter(3))
989   {
990     first = 3;
991     second = 2;
992   }
993   
994   if (!oneGaussian)
995   {
996     Float_t vector[3];
997     for (Int_t i=2; i<5;i++)
998       vector[i-2] = cov[i][2]*func4phi->GetParameter(4) + cov[i][3]*(1-func4phi->GetParameter(4))+cov[i][4]*(TMath::Abs(func4phi->GetParameter(2))-TMath::Abs(func4phi->GetParameter(3)));
999
1000     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))));
1001     Float_t rms = TMath::Abs(func4phi->GetParameter(2))*func4phi->GetParameter(4)+TMath::Abs(func4phi->GetParameter(3))*(1-func4phi->GetParameter(4));
1002     AddPoint(graphs[37][graphID], x, rms, xE, sigma);
1003   }
1004   else AddPoint(graphs[37][graphID], x, TMath::Abs(func4phi->GetParameter(first)), xE, func4phi->GetParError(first));
1005   AddPoint(graphs[38][graphID], x, 0, xE, 0);
1006   AddPoint(graphs[36][graphID], x, 1, xE, 0);
1007
1008   AddPoint(graphs[35][graphID], x, func4phi->GetParameter(1), xE, func4phi->GetParError(1));
1009
1010   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), TMath::Abs(func4phi->GetParameter(first)), 0, func4phi->GetParError(first));
1011   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), TMath::Abs(func4phi->GetParameter(second)), 0, func4phi->GetParError(second));
1012
1013   Float_t scale = func4phi->GetParameter(1);
1014   Float_t weighting = func4phi->GetParameter(4);
1015   func4phi->SetParameter(1, scale * weighting);
1016   func4phi->SetParameter(4, 1);
1017   func4phi->SetLineStyle(2);
1018   func4phi->DrawCopy("SAME");
1019
1020   func4phi->SetParameter(1, scale * (1.0-weighting));
1021   func4phi->SetParameter(4, 0);
1022   func4phi->SetLineStyle(3);
1023   func4phi->DrawCopy("SAME");
1024
1025   fitResultPtr->Print("V");
1026   if (fitResultPtr != 0)
1027     success = kFALSE;
1028
1029   canvas->cd(canvasPos++);
1030   func4eta->SetLineColor(4);
1031   fitResultPtr = projy3->Fit(func4eta, "SIR", "SAME");
1032   cov = fitResultPtr->GetCovarianceMatrix();
1033
1034   projy3->Draw("SAME");
1035   projy3->SetLineColor(4);
1036
1037   chi2 = 0;
1038   ndf = 0;
1039   if (func4eta->GetNDF() > 0)
1040   {
1041     printf("#chi^{2}/ndf = %.1f/%d = %.1f  ", func4eta->GetChisquare(), func4eta->GetNDF(), func4eta->GetChisquare() / func4eta->GetNDF());
1042     AddPoint(graphs[44][graphID], x, func4eta->GetChisquare() / func4eta->GetNDF(), xE, 0);
1043   }
1044   for (Int_t i=projy3->FindBin(-momentFitLimit); i<=projy3->FindBin(momentFitLimit);i++)
1045   {
1046 //    Float_t temp = chi2;
1047     chi2 += TMath::Power((projy3->GetBinContent(i)-func4eta->Eval(projy3->GetBinCenter(i)))/projy3->GetBinError(i),2);
1048 //    cerr << i << "\t" << projy3->GetBinCenter(i) << "\t" << chi2-temp << endl;
1049     ndf++;
1050   }
1051   if (!oneGaussian) ndf = ndf - 4;
1052   else ndf = ndf - 2;
1053   printf("Calculated #chi^{2}/ndf = %.1f/%d = %.1f  ", chi2, ndf, chi2/ndf);
1054
1055   first = 2;
1056   second = 3;
1057   if (func4eta->GetParameter(2) < func4eta->GetParameter(3))
1058   {
1059     first = 3;
1060     second = 2;
1061   }
1062
1063   AddPoint(graphs[39][graphID], x, func4eta->GetParameter(1), xE, func4eta->GetParError(1));
1064   if (!oneGaussian)
1065   {
1066     Float_t vector[3];
1067     for (Int_t i=2; i<5;i++)
1068       vector[i-2] = cov[i][2]*func4eta->GetParameter(4) + cov[i][3]*(1-func4eta->GetParameter(4))+cov[i][4]*(TMath::Abs(func4eta->GetParameter(2))-TMath::Abs(func4eta->GetParameter(3)));
1069
1070     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))));
1071     Float_t rms = TMath::Abs(func4eta->GetParameter(2))*func4eta->GetParameter(4)+TMath::Abs(func4eta->GetParameter(3))*(1-func4eta->GetParameter(4));
1072     AddPoint(graphs[41][graphID], x, rms, xE, sigma);
1073   }
1074   else AddPoint(graphs[41][graphID], x, TMath::Abs(func4eta->GetParameter(first)), xE, func4eta->GetParError(first));
1075   AddPoint(graphs[42][graphID], x, 0, xE, 0);
1076   AddPoint(graphs[40][graphID], x, 1, xE, 0);
1077
1078   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), TMath::Abs(func4eta->GetParameter(first)), 0, func4eta->GetParError(first));
1079   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), TMath::Abs(func4eta->GetParameter(second)), 0, func4eta->GetParError(second));
1080
1081   scale = func4eta->GetParameter(1);
1082   weighting = func4eta->GetParameter(4);
1083   func4eta->SetParameter(1, scale * weighting);
1084   func4eta->SetParameter(4, 1);
1085   func4eta->SetLineStyle(2);
1086   func4eta->DrawCopy("SAME");
1087
1088   func4eta->SetParameter(1, scale * (1.0-weighting));
1089   func4eta->SetParameter(4, 0);
1090   func4eta->SetLineStyle(3);
1091   func4eta->DrawCopy("SAME");
1092   
1093   fitResultPtr->Print("V");
1094   if (fitResultPtr != 0)
1095     success = kFALSE;
1096
1097 //  return 0;
1098
1099   // 2d fit with two gaussians
1100   canvas->cd(canvasPos++);
1101 //   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);
1102 //   func3->SetParameters(0, 1, 0.3, 0.3, 1, 0.6, 0.6);
1103 //   func3->SetParLimits(4, 0, 10);
1104 //   Float_t etaFitLimit = 0.5;
1105   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);
1106
1107   func3->SetParLimits(4, 0.1, 0.9);
1108   func3->SetParLimits(1, 0, 10);
1109   func3->SetParLimits(2, sigmaFitLimit, 0.7);
1110   func3->SetParLimits(3, sigmaFitLimit, etaFitUpperLimit);
1111   func3->SetParLimits(5, sigmaFitLimit, 0.7);
1112   func3->SetParLimits(6, sigmaFitLimit, etaFitUpperLimit);
1113   
1114   func3->FixParameter(0, 0);
1115 /*
1116   sigmaFitLimit = 0.09;
1117   etaFitUpperLimit = 0.12;
1118   func3->SetParLimits(4, 0.1, 0.9);
1119   func3->SetParLimits(1, 0, 10);
1120   func3->SetParLimits(2, sigmaFitLimit, 0.12);
1121   func3->SetParLimits(3, sigmaFitLimit, etaFitUpperLimit);
1122   func3->SetParLimits(5, sigmaFitLimit, 0.12);
1123   func3->SetParLimits(6, sigmaFitLimit, etaFitUpperLimit);
1124   func3->FixParameter(0, 0);
1125 */
1126   func3->SetParameters(0, hist->GetBinContent(hist->GetXaxis()->FindBin(0.0), hist->GetYaxis()->FindBin(0.0)) - mean, 0.3, 0.3, 0.25, initSigma, initSigma);
1127 //   for (Int_t i=0; i<6; i++)
1128 //     func3->SetParameter(i+1, func->GetParameter(i));
1129
1130   if (0 && histId == 0)
1131   {
1132     // central --> go to 1 Gaussian only
1133     func3->SetParLimits(4, 1, 1);
1134     func3->FixParameter(4, 1);
1135     func3->FixParameter(5, sigmaFitLimit);
1136     func3->FixParameter(6, sigmaFitLimit);
1137   }
1138   
1139   // set errors 20% of the value for bins potentially affected by two-track effects
1140 //   hist->SetBinError(hist->GetXaxis()->FindBin(0.0001),  hist->GetYaxis()->FindBin(0.0001),  0.1 * hist->GetBinContent(hist->GetXaxis()->FindBin(0.0001),  hist->GetYaxis()->FindBin(0.0001)));
1141 //   hist->SetBinError(hist->GetXaxis()->FindBin(0.0001),  hist->GetYaxis()->FindBin(-0.0001), 0.1 * hist->GetBinContent(hist->GetXaxis()->FindBin(0.0001),  hist->GetYaxis()->FindBin(-0.0001)));
1142 //   hist->SetBinError(hist->GetXaxis()->FindBin(-0.0001), hist->GetYaxis()->FindBin(0.0001),  0.1 * hist->GetBinContent(hist->GetXaxis()->FindBin(-0.0001), hist->GetYaxis()->FindBin(0.0001)));
1143 //   hist->SetBinError(hist->GetXaxis()->FindBin(-0.0001), hist->GetYaxis()->FindBin(-0.0001), 0.1 * hist->GetBinContent(hist->GetXaxis()->FindBin(-0.0001), hist->GetYaxis()->FindBin(-0.0001)));
1144
1145   fitResult = hist->Fit(func3, "0R", "");
1146   Printf("Fit result: %d", fitResult);
1147   if (fitResult != 0)
1148     success = kFALSE;
1149   
1150   Printf("Testing 1 Gaussian...");
1151
1152   TF2* func3_clone = (TF2*) func3->Clone("func3_clone");
1153   
1154   func3_clone->SetParLimits(4, 1, 1);
1155   func3_clone->FixParameter(4, 1);
1156   func3_clone->FixParameter(5, sigmaFitLimit);
1157   func3_clone->FixParameter(6, sigmaFitLimit);
1158   
1159   fitResult = hist->Fit(func3_clone, "0R", "");
1160   Printf("Fit result: %d", fitResult);
1161   
1162   // if both parameters were within 1%, use 1 Gaussian parameters
1163   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)
1164   {
1165     Printf("Parameters within 1%%. Using results from 1 Gaussian...");
1166     
1167     if (fitResult != 0)
1168       success = kFALSE;
1169     func3 = func3_clone;
1170   }
1171   else if (func3_clone->GetChisquare() * 0.99 < func3->GetChisquare())
1172   {
1173     Printf("1 Gaussian fit has a at maximum 1%% (%f %f) larger chi2. Using results from 1 Gaussian...", func3_clone->GetChisquare(), func3->GetChisquare());
1174     
1175     if (fitResult != 0)
1176       success = kFALSE;
1177     func3 = func3_clone;
1178   }
1179   
1180 //   hist->Fit(func3, "I0R", "");
1181
1182   first = 2;
1183   second = 5;
1184   if (func3->GetParameter(2) < func3->GetParameter(5))
1185   {
1186     first = 5;
1187     second = 2;
1188   }
1189   //dphi
1190   AddPoint(graphs[1][graphID], x, TMath::Abs(func3->GetParameter(first)), xE, func3->GetParError(first));
1191   AddPoint(graphs[4][graphID], x, TMath::Abs(func3->GetParameter(second)), xE, func3->GetParError(second));
1192
1193   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), TMath::Abs(func3->GetParameter(first)), 0, func3->GetParError(first));
1194   AddPoint(phiWidthSummary, phiWidthSummary->GetN(), TMath::Abs(func3->GetParameter(second)), 0, func3->GetParError(second));
1195     
1196   //deta
1197   AddPoint(graphs[2][graphID], x, TMath::Abs(func3->GetParameter(first+1)), xE, func3->GetParError(first+1));
1198   AddPoint(graphs[5][graphID], x, TMath::Abs(func3->GetParameter(second+1)), xE, func3->GetParError(second+1));
1199   
1200   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), TMath::Abs(func3->GetParameter(first+1)), 0, func3->GetParError(first+1));
1201   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), TMath::Abs(func3->GetParameter(second+1)), 0, func3->GetParError(second+1));
1202   
1203   // norm
1204   AddPoint(graphs[0][graphID], x, func3->GetParameter(1), xE, func3->GetParError(1));
1205   if (first < second)
1206     AddPoint(graphs[3][graphID], x, func3->GetParameter(4), xE, func3->GetParError(4));
1207   else
1208     AddPoint(graphs[3][graphID], x, 1.0 - func3->GetParameter(4), xE, func3->GetParError(4));
1209
1210   TH2* funcHist3 = (TH2*) hist->Clone("funcHist3");
1211   funcHist3->Reset();
1212   funcHist3->Add(func3);
1213   funcHist3->SetMinimum(-(max - min) / 2);
1214   funcHist3->SetMaximum((max - min) / 2);
1215   funcHist3->Draw("SURF1");
1216   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), funcHist3->Integral());
1217   
1218   canvas->cd(canvasPos++);
1219   TH2* residuals3 = (TH2*) hist->Clone("residuals");
1220   residuals3->Add(funcHist3, -1);
1221   residuals3->SetMinimum(-(max - min) / 2);
1222   residuals3->SetMaximum((max - min) / 2);
1223   residuals3->Draw("SURF1");
1224
1225   chi2 = 0;
1226   ndf = 0;
1227   for (Int_t i=hist->GetXaxis()->FindBin(-0.8); i<=hist->GetXaxis()->FindBin(0.8); i++)
1228     for (Int_t j=hist->GetYaxis()->FindBin(-etaLimit); j<=hist->GetYaxis()->FindBin(etaLimit); j++)
1229     {
1230       if (residuals3->GetBinError(i, j) > 0)
1231       {
1232         chi2 += TMath::Power(residuals3->GetBinContent(i, j) / residuals3->GetBinError(i, j), 2);
1233         ndf++;
1234       }
1235     }
1236   ndf -= func3->GetNumberFreeParameters();
1237   
1238   if (func3->GetNDF() > 0)
1239   {
1240     printf("#chi^{2}/ndf = %.1f/%d = %.1f  ", func3->GetChisquare(), func3->GetNDF(), func3->GetChisquare() / func3->GetNDF());
1241     DrawLatex(0.2, yPosChi2, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", func3->GetChisquare(), func3->GetNDF(), func3->GetChisquare() / func3->GetNDF()));
1242     AddPoint(graphs[6][graphID], x, func3->GetChisquare() / func3->GetNDF(), xE, 0);
1243   }
1244   if (ndf)
1245   {
1246     Printf("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf);
1247     DrawLatex(0.2, yPosChi2 - 0.05, 1, Form("#chi^{2}/ndf = %.1f/%d = %.1f", chi2, ndf, chi2 / ndf));
1248     AddPoint(graphs[7][graphID], x, chi2 / ndf, xE, 0);
1249   }
1250
1251   canvas->cd(canvasPos++);
1252   phiWidthSummary->SetMarkerStyle(20);
1253   phiWidthSummary->Draw("AP");
1254   gPad->SetGridy();
1255     
1256   etaWidthSummary->SetMarkerStyle(21);
1257   etaWidthSummary->SetLineColor(2);
1258   etaWidthSummary->SetMarkerColor(2);
1259   etaWidthSummary->Draw("PSAME");  
1260   
1261   phiWidthSummary->GetYaxis()->SetRangeUser(0, 0.9);
1262   
1263   canvas->cd(canvasPos++);
1264   sumSummary->Draw("*A");
1265   gPad->SetGridy();
1266   
1267   Printf("Finished with %d", success);
1268
1269   return success;
1270 }
1271
1272 void AnalyzeDeltaPhiEtaGap2D(const char* fileName, const char* outputFileName = "graphs.root")
1273 {
1274   gROOT->SetBatch(kTRUE);
1275   if (!gROOT->IsBatch())
1276   {
1277     Printf("Not in batch mode. Exiting!");
1278     return;
1279   }
1280   
1281   CreateGraphStructure();
1282
1283   TFile::Open(fileName);
1284   TList checkList;
1285   
1286   Int_t maxLeadingPt = 4;
1287   Int_t maxAssocPt = 6;
1288
1289
1290   for (Int_t i=0; i<maxLeadingPt; i++)
1291   {
1292     for (Int_t j=1; j<maxAssocPt; j++)
1293     {
1294       // only process when first is filled
1295       TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, 0));
1296       if (!hist1)
1297       {
1298         cout << "Hist1 does not exist." << endl;
1299         continue;
1300       }
1301 //      if (hist1->GetEntries() < 1e4)
1302       if (hist1->GetEntries() < 10)
1303       {
1304         Printf("%d %d Only %f entries. Skipping...", i, j, hist1->GetEntries());
1305         continue;
1306       }
1307       
1308       for (Int_t histId = 0; histId < NHists; histId++)
1309       {
1310 //      if (i != 1 || j != 2)
1311 //        continue;
1312     
1313         hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
1314         if (!hist1)
1315         {
1316           cout << "Hist1 does not exist." << endl;
1317           continue;
1318         }
1319         
1320 //      if (hist1->GetEntries() < 1e4)
1321         if (hist1->GetEntries() < 10)
1322         {
1323           Printf("%d %d %d Only %f entries. Skipping...", i, j, histId, hist1->GetEntries());
1324           continue;
1325         }
1326         
1327         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);
1328         canvas->Divide(5, 3);
1329         
1330         for (Int_t k=1; k<=3; k++)
1331         {
1332           canvas->cd(3 * j + k);
1333           gPad->SetLeftMargin(0.15);
1334           gPad->SetBottomMargin(0.2);
1335           gPad->SetTopMargin(0.01);
1336           gPad->SetRightMargin(0.01);
1337         }
1338         
1339         Int_t graphID = i * (maxAssocPt - 1) + j - 1;
1340
1341         Printf("\n\n>>> %d %d %d %d", i, j, histId, graphID);
1342         
1343         if (histId == 0)
1344           for (Int_t k=0; k<NGraphs; k++)
1345             graphs[k][graphID]->SetTitle(hist1->GetTitle());
1346         
1347         Float_t centralityAxisMapping[] = { 5, 65, 100, 30, 15, 50 };
1348         Float_t centralityAxisMappingE[] = { 5, 5, 0, 10, 5, 10 };
1349
1350         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));
1351         if (!success)
1352           checkList.Add(new TObjString(Form("AnalyzeDeltaPhiEtaGap2DExample(\"%s\", %d, %d, %d)", fileName, i, j, histId)));
1353         
1354 //      canvas->SaveAs(Form("%s.png", canvas->GetName()));
1355 //      delete canvas;
1356 //      break;
1357 // return;
1358       }
1359       
1360 //       break;
1361     }
1362     
1363 //     break;
1364   }
1365   
1366   WriteGraphs(outputFileName);
1367
1368   for (Int_t i=0; i<checkList.GetEntries(); i++)
1369     Printf("%s", checkList.At(i)->GetName());
1370 }
1371
1372 void AnalyzeDeltaPhiEtaGap2DExample(const char* fileName, Int_t i, Int_t j, Int_t histId, Bool_t drawDetails = kFALSE)
1373 {
1374   CreateGraphStructure();
1375
1376   TFile::Open(fileName);
1377   
1378   TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
1379   if (!hist1)
1380     return;
1381   
1382   Printf("Entries: %f %s", hist1->GetEntries(), hist1->GetTitle());
1383
1384   TString label(hist1->GetTitle());
1385   if (drawDetails)
1386     hist1->SetTitle("");
1387   hist1->GetYaxis()->SetRangeUser(-1.79, 1.79);
1388   hist1->GetXaxis()->SetTitleOffset(1.5);
1389   hist1->GetYaxis()->SetTitleOffset(2);
1390   hist1->GetZaxis()->SetTitle(kCorrFuncTitle);
1391   hist1->GetZaxis()->SetTitleOffset(1.8);
1392   hist1->SetStats(kFALSE);
1393
1394 //  if (hist1->GetEntries() < 1e4)
1395   if (hist1->GetEntries() < 10)
1396     return;
1397   
1398   // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
1399   if (drawDetails)
1400     hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
1401
1402   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);
1403   canvas->Divide(5, 3);
1404   
1405   Int_t graphID = i * (6 - 1) + j - 1;
1406   FitDeltaPhi2DOneFunction((TH2*) hist1, canvas, 1, graphID, 0, 0, 0.9, kFALSE, histId, (i > 0) ? 1 : 0, (j <= 2 && histId != 2));
1407   
1408   if (!drawDetails)
1409     return;
1410   
1411   TVirtualPad* pad = canvas->cd(6);
1412   TCanvas* c = new TCanvas("c", "c", 1500, 1000);
1413   c->Divide(3, 2);
1414   c->cd(1);
1415   pad->SetPad(0, 0, 1, 1);
1416   pad->SetLeftMargin(0.15);
1417   pad->Draw();
1418   pad->cd();
1419
1420   label.ReplaceAll(".00", " GeV/c");
1421   label.ReplaceAll(".0", " GeV/c");
1422   TObjArray* objArray = label.Tokenize("-");
1423   TPaveText* paveText = new TPaveText(0.52, 0.72, 0.95, 0.95, "BRNDC");
1424   paveText->SetTextSize(0.035);
1425   paveText->SetFillColor(0);
1426   paveText->SetShadowColor(0);
1427   paveText->AddText(objArray->At(0)->GetName());
1428   paveText->AddText(objArray->At(1)->GetName());
1429   if (objArray->GetEntries() == 4)
1430     paveText->AddText(Form("Pb-Pb 2.76 TeV %s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()));
1431   else
1432     paveText->AddText(Form("%s 2.76 TeV", objArray->At(2)->GetName()));
1433   paveText->AddText("|#eta| < 0.9");
1434   paveText->Draw();
1435   
1436   c->cd(1);
1437   DrawALICELogo(kTRUE, 0.2, 0.7, 0.4, 0.9);
1438   
1439 //   return;
1440     
1441   //   c->SaveAs("fit_subtracted.eps");
1442
1443   pad = canvas->cd(12);
1444   c->cd(2);
1445   pad->SetPad(0, 0, 1, 1);
1446   pad->SetLeftMargin(0.15);
1447   pad->Draw();
1448 //   c->SaveAs("fit_fit1.eps");
1449   
1450   pad = canvas->cd(13);
1451   c->cd(3);
1452   gPad->SetLeftMargin(0.15);
1453   pad->SetPad(0, 0, 1, 1);
1454   pad->SetLeftMargin(0.15);
1455   pad->Draw();
1456   c->SaveAs("fit_residual1.eps");
1457   
1458   TCanvas* c2 = new TCanvas("c3b", "c3b", 800, 800);
1459   gPad->SetLeftMargin(0.15);
1460   gPad->SetGridy();
1461   hist1 = (TH2*) pad->GetListOfPrimitives()->First();
1462   for (i=0; i<10; i++)
1463   {
1464     TH1* proj = hist1->ProjectionX(Form("p_%d", i), 11+i*2, 12+i*2);
1465     proj->Scale(0.5);
1466     
1467     proj->Add(new TF1("func", "1", -10, 10), -0.5 + 0.1 * i);
1468     proj->SetStats(0);
1469     proj->GetXaxis()->SetTitleOffset(1.0);
1470     proj->GetYaxis()->SetTitleOffset(1.4);
1471     proj->GetYaxis()->SetNdivisions(512);
1472     proj->SetYTitle("Residuals (projected / shifted)");
1473     proj->GetYaxis()->SetRangeUser(-0.59, 0.49);
1474     proj->Draw((i == 0) ? "" : "SAME");    
1475   }
1476   c2->SaveAs("fit_residual1_proj.eps");
1477
1478   pad = canvas->cd(4);
1479   c->cd(4);
1480   pad->SetPad(0, 0, 1, 1);
1481   pad->SetLeftMargin(0.15);
1482   pad->Draw();
1483 //   c->SaveAs("fit_fit2.eps");
1484
1485   pad = canvas->cd(3);
1486   c->cd(5);
1487   pad->SetPad(0, 0, 1, 1);
1488   pad->SetLeftMargin(0.15);
1489   pad->Draw();
1490 //   c->SaveAs("fit_residual2.eps");
1491
1492   c2 = new TCanvas("c5b", "c5b", 800, 800);
1493   gPad->SetLeftMargin(0.15);
1494   gPad->SetGridy();
1495   hist1 = (TH2*) pad->GetListOfPrimitives()->First();
1496   for (i=0; i<10; i++)
1497   {
1498     TH1* proj = hist1->ProjectionX(Form("p2_%d", i), 11+i*2, 12+i*2);
1499     proj->Scale(0.5);
1500     
1501     proj->Add(new TF1("func", "1", -10, 10), -0.5 + 0.1 * i);
1502     proj->SetStats(0);
1503     proj->GetXaxis()->SetTitleOffset(1.0);
1504     proj->GetYaxis()->SetTitleOffset(1.4);
1505     proj->GetYaxis()->SetNdivisions(512);
1506     proj->SetYTitle("Residuals (projected / shifted)");
1507     proj->GetYaxis()->SetRangeUser(-0.59, 0.49);
1508     proj->Draw((i == 0) ? "" : "SAME");    
1509   }
1510   c2->SaveAs("fit_residual2_proj.eps");
1511
1512   pad = canvas->cd(7);
1513   c->cd(6);
1514   pad->SetPad(0, 0, 1, 1);
1515   pad->SetLeftMargin(0.15);
1516   pad->Draw();
1517   
1518   c->SaveAs("fit_example.eps");
1519   c->SaveAs("fit_example.png");
1520 }
1521
1522 void SqrtAll(Int_t nHists, TGraphErrors** graph)
1523 {
1524   for (Int_t histId = 0; histId < nHists; histId++)
1525   {
1526     for (Int_t i=0; i<graph[histId]->GetN(); i++)
1527     {
1528       if (graph[histId]->GetY()[i] > 0)
1529       {
1530         graph[histId]->GetEY()[i] *= 0.5 / TMath::Sqrt(graph[histId]->GetY()[i]);
1531         graph[histId]->GetY()[i] = TMath::Sqrt(graph[histId]->GetY()[i]);
1532       }
1533       else
1534       {
1535         Printf("WARNING negative value %d %d, removing from histogram", histId, i);
1536         graph[histId]->RemovePoint(i);
1537         i--;
1538       }
1539     }
1540   }
1541 }
1542
1543 void SqrtAll2(Int_t nHists, TGraphErrors** graph, TGraphErrors** graph2)
1544 {
1545   for (Int_t histId = 0; histId < nHists; histId++)
1546   {
1547     if (!graph[histId])
1548       continue;
1549  
1550     for (Int_t i=0; i<graph[histId]->GetN(); i++)
1551     {
1552       if (graph[histId]->GetY()[i] > 0 && graph2[histId]->GetY()[i] > 0)
1553       {
1554         graph[histId]->GetEY()[i] *= 0.5 / TMath::Sqrt(graph[histId]->GetY()[i]);
1555         graph[histId]->GetY()[i] = TMath::Sqrt(graph[histId]->GetY()[i]);
1556
1557         graph2[histId]->GetEY()[i] *= 0.5 / TMath::Sqrt(graph2[histId]->GetY()[i]);
1558         graph2[histId]->GetY()[i] = TMath::Sqrt(graph2[histId]->GetY()[i]);
1559       }
1560       else
1561       {
1562         Printf("WARNING negative value %d %d, removing from histogram", histId, i);
1563         graph[histId]->RemovePoint(i);
1564         graph2[histId]->RemovePoint(i);
1565         i--;
1566       }
1567     }
1568   }
1569 }
1570
1571 Int_t marker[6] = { 24, 25, 26, 30, 27, 28 };
1572 Int_t marker2[6] = { 20, 21, 22, 29, 33, 34 };
1573 const char* labels[6] = { "0-10%", "60-70%", "pp", "20-40%", "10-20%", "40-60%" };
1574
1575 void CalculateRMS(Int_t nHists, TGraphErrors** graph, TGraphErrors** graph2, TGraphErrors** weight)
1576 {
1577   // rms
1578   for (Int_t histId = 0; histId < nHists; histId++)
1579   {
1580     if (!graph[histId])
1581       continue;
1582     
1583     if (graph[histId]->GetN() != graph2[histId]->GetN() || graph[histId]->GetN() != weight[histId]->GetN())
1584     {
1585       Printf("E-CalculateRMS: Different number of points for hist %d: %d %d %d", histId, graph[histId]->GetN(), graph2[histId]->GetN(), weight[histId]->GetN());
1586       continue;
1587     }
1588     
1589     for (Int_t i=0; i<graph[histId]->GetN(); i++)
1590     {
1591       Float_t rms = graph[histId]->GetY()[i] * weight[histId]->GetY()[i] + graph2[histId]->GetY()[i] * (1.0 - weight[histId]->GetY()[i]);
1592       graph[histId]->GetY()[i] = rms;
1593       // error**2 = weight**2 * error_sigma**2 + (weight-1)**2 * error_sigma2**2 + (sigma1 + sigma2)**2 * error_weight**2
1594       // TODO this neglects some correlations
1595       graph[histId]->GetEY()[i] = TMath::Sqrt(TMath::Power(weight[histId]->GetY()[i] * graph[histId]->GetEY()[i], 2) + 
1596         TMath::Power((weight[histId]->GetY()[i] - 1) * graph2[histId]->GetEY()[i], 2) + TMath::Power((graph[histId]->GetY()[i] + graph2[histId]->GetY()[i]) * weight[histId]->GetEY()[i], 2));
1597     }
1598   }
1599 }
1600
1601 Bool_t SkipGraph(Int_t i)
1602 {
1603 //  return (i == 1 || i == 7 || i == 9 || i == 13 || i == 14);
1604   return (i == 7 || i == 9 || i == 13 || i == 14); //For the STAR comparison
1605 }
1606
1607 Int_t skipGraphList[] =  { -1, -1, -1 };
1608 Bool_t SkipGraphForThisPlot(Int_t graphId)
1609 {
1610   for (Int_t i=0; i<3; i++)
1611     if (skipGraphList[i] == graphId)
1612       return kTRUE;
1613   return kFALSE;
1614 }
1615
1616 TGraphAsymmErrors* FixGraph(TGraphErrors* graph, Float_t shift)
1617 {
1618   graph->Sort();
1619 //   graph->Print();
1620   
1621   if (graph->GetN() > 4 && graph->GetX()[4] == 75)
1622   {
1623     graph->GetX()[4] = 65;
1624     graph->GetEX()[4] = 5;
1625   }
1626   
1627   TGraphAsymmErrors* res = new TGraphAsymmErrors;
1628   res->SetTitle(graph->GetTitle());
1629   res->GetXaxis()->SetTitle(graph->GetXaxis()->GetTitle());
1630   res->GetYaxis()->SetTitle(graph->GetYaxis()->GetTitle());
1631
1632   for (Int_t i=0; i<graph->GetN(); i++)
1633   {
1634     res->SetPoint(res->GetN(), graph->GetX()[i] + shift, graph->GetY()[i]);
1635     res->SetPointError(res->GetN()-1, graph->GetEX()[i], graph->GetEX()[i], graph->GetEY()[i], graph->GetEY()[i]);
1636     
1637     // asymmetric if space
1638     if (graph->GetEX()[i] > TMath::Abs(shift))
1639     {
1640       res->GetEXlow()[i] += shift;
1641       res->GetEXhigh()[i] -= shift;
1642     }
1643   }
1644   
1645   return res;
1646 }
1647
1648 void GraphRemovePoint(TGraphAsymmErrors* graph1, TGraphAsymmErrors* graph2)
1649 {
1650   if (graph1->GetN() != graph2->GetN())
1651   {
1652     graph1->Sort();
1653     graph2->Sort();
1654     Int_t Nbin = 0;
1655     if (graph1->GetN() < graph2->GetN()) Nbin = graph1->GetN();
1656     else Nbin = graph2->GetN();
1657     for (Int_t bin1 = 0; bin1 < Nbin; bin1++)
1658     {
1659       Float_t x1 = graph1->GetX()[bin1]; 
1660       Float_t x2 = graph2->GetX()[bin1];
1661       if (x1 != x2)
1662       {
1663         if (x1<x2) {graph1->RemovePoint(bin1--);cout << x1 << "\t" << x2 << "\t" << bin1 << "\t" << graph1->GetN() << "\t" << graph2->GetN()<< endl;}
1664         else {graph2->RemovePoint(bin1--);cout << x1 << "\t" << x2 << "\t" << bin1 << "\t" << graph1->GetN() << "\t" << graph2->GetN()<< endl;}
1665       }
1666     }
1667   }
1668 }
1669
1670 void PrepareGraphs(Int_t nHists, TGraphErrors** graph, TGraphErrors** systematicA, TGraphErrors** systematicB, TMultiGraph** multiGraph, TMultiGraph** multiGraphSyst, Int_t uncertaintyID, Float_t offset = 0)
1671 {
1672   Int_t colors[16] =  { 1, 3, 2, 6, 4, 7, 8, 9, 11, 12, 28, 30, 36, 40, 46 };
1673   Int_t markers[16] = { 20, 21, 22, 23, 24, 25, 26, 27, 28, 30, 31, 32, 33, 34, 2, 5};
1674 //  Int_t fillStyle[11] = { 3003, 3003, 3003, 3003, 3003, 3003, 3003, 3003, 3003, 3003, 3003 };
1675   Int_t fillStyle[11] = { 3001, 3002, 3003, 3004, 3005, 3006, 3007, 3008, 3009, 3010, 3011 };
1676   Int_t count = 0;
1677   
1678   if (*multiGraph == 0)
1679     *multiGraph = new TMultiGraph;
1680   if (*multiGraphSyst == 0)
1681     *multiGraphSyst = new TMultiGraph;
1682   
1683   Float_t shift = -2 + offset;
1684   for (Int_t i=0; i<nHists; i++)
1685   {
1686     if (SkipGraph(i))
1687       continue;
1688     
1689     if (SkipGraphForThisPlot(i))
1690     {
1691       count++;
1692       continue;
1693     }
1694     
1695     TGraphAsymmErrors* graphcentrality = FixGraph((TGraphErrors*) graph[i]->Clone(), shift);
1696     if (graphcentrality->GetN() <= 0)
1697       continue;
1698     
1699     if (systematicA)
1700     {
1701       TGraphAsymmErrors* graphsystematicsA = FixGraph((TGraphErrors*) systematicA[i]->Clone(), shift);
1702       
1703       if (graphcentrality->GetN() != graphsystematicsA->GetN())
1704       {
1705         Printf("Different number of points %d %d: %s", graphcentrality->GetN(), graphsystematicsA->GetN(), graphcentrality->GetTitle());
1706         if (graphcentrality->GetN() == 0 || graphsystematicsA->GetN() == 0) continue;
1707         GraphRemovePoint(graphcentrality,graphsystematicsA);
1708       }
1709       
1710       TGraphErrors* graphsystematicsB = 0;
1711       if (systematicB)
1712       {
1713         graphsystematicsB = (TGraphErrors*) systematicB[i]->Clone();
1714         graphsystematicsB->Sort();
1715       
1716         if (graphcentrality->GetN() != graphsystematicsB->GetN())
1717         {
1718           Printf("Different number of points %d %d", graphcentrality->GetN(), graphsystematicsB->GetN());
1719           continue;
1720         }
1721       }
1722       
1723       for (Int_t j=0; j<graphsystematicsA->GetN(); j++)
1724       {
1725         // uncertaintyID
1726         Double_t yMin = graphcentrality->GetY()[j];
1727         Double_t yMax = graphcentrality->GetY()[j];
1728         
1729         if (uncertaintyID == 1)
1730         {
1731           // 5%
1732           yMin *= 0.95;
1733           yMax *= 1.05;
1734         }
1735         else if (uncertaintyID == 2)
1736         {
1737           yMin -= 0.20;
1738           yMax += 0.20;
1739         }
1740         
1741         if (graphsystematicsA->GetY()[j] < graphcentrality->GetY()[j])
1742           yMin = graphcentrality->GetY()[j] - TMath::Sqrt(TMath::Power(yMin - graphcentrality->GetY()[j], 2) + TMath::Power(graphsystematicsA->GetY()[j] - graphcentrality->GetY()[j], 2));
1743         if (graphsystematicsA->GetY()[j] > graphcentrality->GetY()[j])
1744           yMax = graphcentrality->GetY()[j] + TMath::Sqrt(TMath::Power(yMax - graphcentrality->GetY()[j], 2) + TMath::Power(graphsystematicsA->GetY()[j] - graphcentrality->GetY()[j], 2));
1745         
1746         if (graphsystematicsB)
1747         {
1748           if (graphsystematicsB->GetY()[j] < graphcentrality->GetY()[j])
1749             yMin = graphcentrality->GetY()[j] - TMath::Sqrt(TMath::Power(yMin - graphcentrality->GetY()[j], 2) + TMath::Power(graphsystematicsB->GetY()[j] - graphcentrality->GetY()[j], 2));
1750           if (graphsystematicsB->GetY()[j] > graphcentrality->GetY()[j])
1751             yMax = graphcentrality->GetY()[j] + TMath::Sqrt(TMath::Power(yMax - graphcentrality->GetY()[j], 2) + TMath::Power(graphsystematicsB->GetY()[j] - graphcentrality->GetY()[j], 2));
1752         }
1753
1754         graphsystematicsA->GetEYlow()[j] = graphcentrality->GetY()[j] - yMin;
1755         graphsystematicsA->GetEYhigh()[j] = yMax - graphcentrality->GetY()[j];
1756         graphsystematicsA->GetY()[j]  = graphcentrality->GetY()[j];
1757         
1758         graphsystematicsA->GetEXlow()[j] = 1;
1759         graphsystematicsA->GetEXhigh()[j] = graphsystematicsA->GetEXlow()[j];
1760       }
1761       
1762 //       graphsystematicsA->SetFillColor(kGray);
1763       graphsystematicsA->SetFillColor(colors[count]);
1764 //       graphsystematicsA->SetFillStyle(1001);
1765       graphsystematicsA->SetFillStyle(fillStyle[count]);
1766       graphsystematicsA->SetMarkerStyle(0);
1767       graphsystematicsA->SetLineColor(0);
1768       (*multiGraphSyst)->Add(graphsystematicsA, "2");
1769     }
1770     
1771     TString label = graphcentrality->GetTitle();
1772     if (label.Length() > 0)
1773     {
1774       TObjArray* tokens = label.Tokenize("-");
1775       label.Form("%s-%s", tokens->At(0)->GetName(),tokens->At(1)->GetName());
1776       label.ReplaceAll(".00", "");
1777       label.ReplaceAll(".0", "");
1778       label.ReplaceAll("p_{T,trig}", "p_{T,t}");
1779       label.ReplaceAll("p_{T,assoc}", "p_{T,a}");
1780     }
1781     label += " GeV/c";
1782     graphcentrality->SetTitle(label);
1783     Printf("%d %s %d", i, label.Data(), colors[count]);
1784     
1785     graphcentrality->SetMarkerStyle(markers[count]);
1786     graphcentrality->SetMarkerColor(colors[count]);
1787     graphcentrality->SetLineColor(colors[count]);
1788     
1789     (*multiGraph)->Add(graphcentrality);
1790     shift += 1;
1791    
1792     count++;
1793   }
1794 }
1795
1796
1797 Bool_t drawLogo = kFALSE;
1798 const char* MCLabel = 0;
1799
1800 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)
1801 {
1802   Bool_t found = kTRUE;
1803   TCanvas* c1 = 0;//(TCanvas*) gROOT->GetListOfCanvases()->FindObject(canvasName);
1804   if (!c1)
1805   {
1806     c1 = new TCanvas(canvasName, canvasName, 800, 600);
1807     found = kFALSE;
1808   }
1809   c1->cd();
1810   
1811   TMultiGraph* multiGraph = 0;
1812   TMultiGraph* multiGraph2 = 0;
1813   TMultiGraph* multiGraphSyst = 0;
1814   
1815   PrepareGraphs(nHists, graph, systematicA, systematicB, &multiGraph, &multiGraphSyst, uncertaintyID);
1816   if (!multiGraph->GetListOfGraphs())
1817     return;
1818
1819   if (graph2)
1820   {
1821     PrepareGraphs(nHists, graph2, systematic2, 0, &multiGraph2, &multiGraphSyst, uncertaintyID, 0.5);
1822     for (Int_t i=0; i<multiGraph2->GetListOfGraphs()->GetEntries(); i++)
1823       ((TGraph*)multiGraph2->GetListOfGraphs()->At(i))->SetLineWidth(2);
1824     
1825     while (multiGraph2->GetListOfGraphs()->GetEntries() > multiGraph->GetListOfGraphs()->GetEntries())
1826       multiGraph2->GetListOfGraphs()->RemoveAt(multiGraph->GetListOfGraphs()->GetEntries());
1827     
1828     TMultiGraph* multiGraph3 = 0;
1829     if (graph3)
1830     {
1831       PrepareGraphs(nHists, graph3, 0, 0, &multiGraph3, &multiGraphSyst, uncertaintyID);
1832       if (multiGraph3->GetListOfGraphs())
1833       {
1834         for (Int_t i=0; i<multiGraph3->GetListOfGraphs()->GetEntries(); i++)
1835         {
1836           ((TGraph*)multiGraph3->GetListOfGraphs()->At(i))->SetLineWidth(2);
1837           ((TGraph*)multiGraph3->GetListOfGraphs()->At(i))->SetLineStyle(2);
1838         }
1839         
1840         while (multiGraph3->GetListOfGraphs()->GetEntries() > multiGraph->GetListOfGraphs()->GetEntries())
1841           multiGraph3->GetListOfGraphs()->RemoveAt(multiGraph->GetListOfGraphs()->GetEntries());
1842       }
1843     }
1844       
1845 //     multiGraphSyst->Add(multiGraph2, "LX");
1846     
1847 //     multiGraph->Add(multiGraph2);
1848     
1849     if (graph3)
1850       multiGraph2->Add(multiGraph3, "LX");
1851   }
1852
1853 //   multiGraphSyst->Add(multiGraph, (found) ? "LX" : "P");
1854   
1855   // draw by hand, multigraph draw spoils the pngs for some reason, same thing happens with this long code, have changed to gif below...
1856   
1857   TGraphErrors* first = (TGraphErrors*) multiGraph->GetListOfGraphs()->First();
1858   TH2* dummy = new TH2F("dummy", ";Centrality", 100, -2, 115, 100, min, max);
1859   dummy->SetStats(0);
1860   dummy->GetYaxis()->SetTitle(yLabel);
1861   dummy->GetYaxis()->SetTitleOffset(1.2);
1862   dummy->DrawCopy();
1863   
1864   if (multiGraphSyst->GetListOfGraphs())
1865   {
1866     multiGraphSyst->GetListOfGraphs()->First()->Draw((!first) ? "A2" : "2SAME");
1867     if (!first)
1868       first = (TGraphErrors*) multiGraphSyst->GetListOfGraphs()->First();
1869     
1870     for (Int_t i=1; i<multiGraphSyst->GetListOfGraphs()->GetEntries(); i++)
1871       multiGraphSyst->GetListOfGraphs()->At(i)->Draw("2SAME");
1872   }
1873   
1874   if (multiGraph2 && multiGraph2->GetListOfGraphs())
1875   {
1876     for (Int_t i=0; i<multiGraph2->GetListOfGraphs()->GetEntries(); i++)
1877     {
1878       TGraphAsymmErrors* graphTmp = (TGraphAsymmErrors*) multiGraph2->GetListOfGraphs()->At(i);
1879       for (Int_t j=0; j<graphTmp->GetN(); j++)
1880       {
1881         graphTmp->GetEXlow()[j] = 0;
1882         graphTmp->GetEXhigh()[j] = 0;
1883       }
1884       graphTmp->SetMarkerSize(0);
1885     }
1886
1887     multiGraph2->GetListOfGraphs()->First()->Draw((!first) ? "ALZ" : "LXSAME");
1888     if (!first)
1889       first = (TGraphErrors*) multiGraph2->GetListOfGraphs()->First();
1890     
1891     for (Int_t i=1; i<multiGraph2->GetListOfGraphs()->GetEntries(); i++)
1892       multiGraph2->GetListOfGraphs()->At(i)->Draw("LZSAME");
1893   }
1894
1895   if (multiGraph && multiGraph->GetListOfGraphs())
1896   {
1897     multiGraph->GetListOfGraphs()->First()->Draw((!first) ? "AP" : "PSAME");
1898     if (!first)
1899       first = (TGraphErrors*) multiGraph->GetListOfGraphs()->First();
1900     
1901     for (Int_t i=1; i<multiGraph->GetListOfGraphs()->GetEntries(); i++)
1902       multiGraph->GetListOfGraphs()->At(i)->Draw("PSAME");
1903   }
1904
1905 //  TPaveText* paveText = new TPaveText(0.826, 0.059, 0.895, 0.094, "BRNDC");
1906   TPaveText* paveText = new TPaveText(0.77, 0.059, 0.84, 0.094, "BRNDC");
1907   paveText->SetTextSize(0.04);
1908   paveText->SetFillColor(0);
1909   paveText->SetShadowColor(0);
1910   paveText->SetLineColor(0);
1911   paveText->AddText("pp");
1912   paveText->Draw();
1913     
1914   TLegend* legend = new TLegend(0.55, 0.65, 0.95, 0.95);
1915   legend->SetFillColor(0);
1916   legend->SetTextSize(0.03);
1917   
1918   if (graph2)
1919   {
1920     legend->SetNColumns(2);
1921     legend->SetHeader(MCLabel);
1922   }
1923   
1924   for (Int_t i=0; i<multiGraph->GetListOfGraphs()->GetEntries(); i++)
1925   {
1926     if (graph2)
1927     {
1928       legend->AddEntry(multiGraph->GetListOfGraphs()->At(i), " ", "P");
1929       legend->AddEntry(multiGraph2->GetListOfGraphs()->At(i), 0, "L");
1930     }
1931     else
1932       legend->AddEntry(multiGraph->GetListOfGraphs()->At(i), 0, "P");
1933   }
1934   legend->Draw();
1935   
1936   gPad->SetGridx();
1937   gPad->SetGridy();
1938   
1939   DrawLatex(0.13, 0.85,  1, "Pb-Pb #sqrt{s_{NN}} = 2.76 TeV", 0.03);
1940   DrawLatex(0.13, 0.81,  1, "pp #sqrt{s} = 2.76 TeV", 0.03);
1941   DrawLatex(0.13, 0.77, 1, "|#eta| < 0.9", 0.03);
1942   TString text;
1943   TString text2;
1944   if (TString(canvasName).BeginsWith("sigma_phi") || TString(canvasName).BeginsWith("kurtosisphi"))
1945   {
1946     text = "Projected within |#Delta#eta| < 0.80";
1947     text2 = "Calculated within |#Delta#varphi| < 0.87";
1948   }
1949   if (TString(canvasName).BeginsWith("sigma_eta") || TString(canvasName).BeginsWith("kurtosiseta"))
1950   {
1951     text = "Projected within |#Delta#varphi| < 0.87";
1952     text2 = "Calculated within |#Delta#eta| < 0.80";
1953   }
1954   if (text.Length() > 0)
1955   {
1956     DrawLatex(0.13, 0.73, 1, text, 0.03);
1957     DrawLatex(0.13, 0.69, 1, text2, 0.03);
1958   }
1959   
1960   if (0 && MCLabel)
1961   {
1962     DrawLatex(0.13, 0.18, 1, "Points: Data", 0.03);
1963     DrawLatex(0.13, 0.14, 1, Form("Lines: %s", MCLabel), 0.03);
1964   }
1965   
1966   if (drawLogo)
1967     DrawALICELogo(kTRUE, 0.41, 0.75, 0.54, 0.95);
1968   
1969   gSystem->mkdir(fgFolder, kTRUE);
1970   c1->SaveAs(Form("%s/%s.gif", fgFolder.Data(), canvasName));
1971   c1->SaveAs(Form("%s/%s.eps", fgFolder.Data(), canvasName));
1972 }
1973
1974 void CalculateRMSSigma(TGraphErrors*** graphsTmp = 0)
1975 {
1976   if (!graphsTmp)
1977     graphsTmp = graphs;
1978   
1979   CalculateRMS(NHists, graphsTmp[1], graphsTmp[4], graphsTmp[3]);
1980   CalculateRMS(NHists, graphsTmp[2], graphsTmp[5], graphsTmp[3]);
1981
1982   CalculateRMS(NHists, graphsTmp[1+16], graphsTmp[4+16], graphsTmp[3+16]);
1983   CalculateRMS(NHists, graphsTmp[2+16], graphsTmp[5+16], graphsTmp[3+16]);
1984
1985   // sqrt(moment2) = sigma
1986   SqrtAll2(NHists, graphsTmp[8], graphsTmp[8+16]);
1987   SqrtAll2(NHists, graphsTmp[9], graphsTmp[9+16]);
1988 }
1989
1990 // void CombineSyst(const char* systFile)
1991 // {
1992 //   TGraphErrors*** current = graphs;
1993 //   
1994 //   ReadGraphs(systFile);
1995 //   CalculateRMSSigma();
1996 //   TGraphErrors*** graphsSyst = graphs;
1997 //   graphs = current;
1998 //   
1999 //   const Int_t NGraphList = 6;
2000 //   Int_t graphList[] = { 1, 2, 8, 9, 14, 15 };
2001 //   for (Int_t i=0; i<NGraphList; i++)
2002 //   {
2003 //     for (Int_t j=0; j<NHists; j++)
2004 //     {
2005 //       
2006 //     graphs[i+16][j]
2007 //     graphsSyst[i][j]
2008 //     }
2009 //   }
2010 //   
2011 // 
2012 // }
2013
2014 void DrawResultsCentrality(const char* fileName = "graphs.root", const char* fileNameWingRemoved = 0, Int_t offset = 0)
2015 {
2016   TGraphErrors*** graphsWingRemoved = 0;
2017   if (fileNameWingRemoved)
2018   {
2019     ReadGraphs(fileNameWingRemoved);
2020     graphsWingRemoved = graphs;
2021   }
2022   
2023   ReadGraphs(fileName);
2024   
2025   Int_t nHists = 24; //NHists;
2026   
2027   if (1)
2028   {
2029 /*    DrawCentrality("norm", nHists, graphs[0+offset], 0, 0.2, "N (a.u.)", (graphsWingRemoved) ? graphsWingRemoved[0+offset] : 0);
2030     
2031 //     return;
2032     
2033     DrawCentrality("width_phi1_centrality", nHists, graphs[1+offset], 0, 0.8, "#sigma_{#Delta#varphi, 1} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[1+offset] : 0);
2034 //     return;
2035     DrawCentrality("width_phi2_centrality", nHists, graphs[4+offset], 0, 0.8, "#sigma_{#Delta#varphi, 2} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[4+offset] : 0);
2036     DrawCentrality("width_eta1_centrality", nHists, graphs[2+offset], 0, 0.8, "#sigma_{#Delta#eta, 1} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[2+offset] : 0);
2037     DrawCentrality("width_eta2_centrality", nHists, graphs[5+offset], 0, 0.8, "#sigma_{#Delta#eta, 2} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[5+offset] : 0);
2038 */    
2039 /*    DrawCentrality("norm_phi(1D)", nHists, graphs[35], 0, 1, "N (a.u.)", (graphsWingRemoved) ? graphsWingRemoved[0+offset] : 0);
2040     DrawCentrality("norm_eta(1D)", nHists, graphs[39], 0, 1, "N (a.u.)", (graphsWingRemoved) ? graphsWingRemoved[0+offset] : 0);
2041 */    
2042 //     DrawCentrality("width_phi1_centrality(1D)", nHists, graphs[37], 0, 0.8, "#sigma_{#Delta#varphi, 1} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[1+offset] : 0);
2043 //    DrawCentrality("width_phi2_centrality(1D)", nHists, graphs[38], 0, 0.8, "#sigma_{#Delta#varphi, 2} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[4+offset] : 0);
2044 //     DrawCentrality("width_eta1_centrality(1D)", nHists, graphs[41], 0, 0.8, "#sigma_{#Delta#eta, 1} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[2+offset] : 0);
2045 //    DrawCentrality("width_eta2_centrality(1D)", nHists, graphs[42], 0, 0.8, "#sigma_{#Delta#eta, 2} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[5+offset] : 0);
2046 //    DrawCentrality("norm2phi(1D)", nHists, graphs[36], 0, 1.5, "N (a.u.)", (graphsWingRemoved) ? graphsWingRemoved[0+offset] : 0);
2047 //    DrawCentrality("norm2eta(1D)", nHists, graphs[40], 0, 1.5, "N (a.u.)", (graphsWingRemoved) ? graphsWingRemoved[0+offset] : 0);
2048
2049     CalculateRMSSigma();
2050     if (graphsWingRemoved)
2051       CalculateRMSSigma(graphsWingRemoved);
2052
2053     DrawCentrality("phi_rms_centrality", nHists, graphs[1+offset], 0, 0.8, Form("#sigma_{#Delta#varphi} (%s) (rad.)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[1+offset] : 0, 0, 0, 0, 0, 1);
2054
2055     DrawCentrality("eta_rms_centrality", 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, 1);
2056
2057 //     return;
2058 //   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);
2059
2060 //   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);
2061 //   DrawCentrality("chi2_phi(1D)", nHists, graphs[43], 0, 150, "#chi^{2}/ndf");
2062 //   DrawCentrality("chi2_eta(1D)", nHists, graphs[44], 0, 150, "#chi^{2}/ndf");
2063 /*    
2064     DrawCentrality("chi2_1", nHists, graphs[6+offset], 0.5, 5, "#chi^{2}/ndf (full region)");
2065     DrawCentrality("chi2_2", nHists, graphs[7+offset], 0.5, 5, "#chi^{2}/ndf (peak region)");
2066     
2067     DrawCentrality("sigma_phi", nHists, graphs[8+offset], 0, 1, "#sigma_{#Delta#varphi} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[8+offset] : 0, 0, 0, 0, 0, 1);
2068     DrawCentrality("sigma_eta", nHists, graphs[9+offset], 0, 1, "#sigma_{#Delta#eta}", (graphsWingRemoved) ? graphsWingRemoved[9+offset] : 0, 0, 0, 0, 0, 1);
2069
2070 //     DrawCentrality("moment3_phi", nHists, graphs[10+offset], -0.5, 0.5, "moment3 #varphi (rad.)", (graphsWingRemoved) ? graphsWingRemoved[10+offset] : 0, 0, 0, 0, 0, 0);
2071 //     DrawCentrality("moment3_eta", nHists, graphs[11+offset], -0.5, 0.5, "moment3 #eta", (graphsWingRemoved) ? graphsWingRemoved[11+offset] : 0, 0, 0, 0, 0, 0);
2072 */
2073 /*    DrawCentrality("IAAFit", nHists, graphs[32+offset], 0, 20, "I_AA", (graphsWingRemoved) ? graphsWingRemoved[32+offset] : 0);
2074     DrawCentrality("Yield", nHists, graphs[33+offset], 0, 20, "Yield", (graphsWingRemoved) ? graphsWingRemoved[33+offset] : 0);
2075     DrawCentrality("IAAHist", nHists, graphs[34+offset], 0, 20, "I_AA", (graphsWingRemoved) ? graphsWingRemoved[34+offset] : 0);
2076     DrawCentrality("IAA", nHists, graphs[34+offset], 0, 20, "I_AA", graphs[32+offset]);
2077 */
2078   }
2079   
2080 //  DrawCentrality("kurtosisphi_centrality", 12, graphs[14+offset], -1.5, 2.5, "Kurtosis #Delta#varphi", (graphsWingRemoved) ? graphsWingRemoved[14+offset] : 0, 0, 0, 0, 0, 2);
2081 //  DrawCentrality("kurtosiseta_centrality", 12, graphs[15+offset], -1.5, 2.5, "Kurtosis #Delta#eta", (graphsWingRemoved) ? graphsWingRemoved[15+offset] : 0, 0, 0, 0, 0, 2);
2082 }
2083
2084 /*
2085 void DrawResultsCentrality(const char* fileName = "graphs.root", const char* fileNameWingRemoved = 0, Int_t offset = 0)
2086 {
2087   TGraphErrors*** graphsWingRemoved = 0;
2088   if (fileNameWingRemoved)
2089   {
2090     ReadGraphs(fileNameWingRemoved);
2091     graphsWingRemoved = graphs;
2092   }
2093   
2094   ReadGraphs(fileName);
2095   
2096   Int_t nHists = 13; //NHists;
2097   
2098   if (1)
2099   {
2100     DrawCentrality("norm", nHists, graphs[0], 0, 0.2, "N (a.u.)", graphs[0+16], (graphsWingRemoved) ? graphsWingRemoved[0] : 0);
2101     
2102 //     return;
2103     
2104     DrawCentrality("width_phi1_centrality", nHists, graphs[1], 0, 0.8, "#sigma_{#varphi, 1} (rad.)", graphs[1+16], (graphsWingRemoved) ? graphsWingRemoved[1] : 0);
2105 //     return;
2106     DrawCentrality("width_phi2_centrality", nHists, graphs[4], 0, 0.8, "#sigma_{#varphi, 2} (rad.)", graphs[4+16], (graphsWingRemoved) ? graphsWingRemoved[4] : 0);
2107     DrawCentrality("width_eta1_centrality", nHists, graphs[2], 0, 0.8, "#sigma_{#eta, 1} (rad.)", graphs[2+16], (graphsWingRemoved) ? graphsWingRemoved[2] : 0);
2108     DrawCentrality("width_eta2_centrality", nHists, graphs[5], 0, 0.8, "#sigma_{#eta, 2} (rad.)", graphs[5+16], (graphsWingRemoved) ? graphsWingRemoved[5] : 0);
2109     
2110     CalculateRMSSigma();
2111     if (graphsWingRemoved)
2112       CalculateRMSSigma(graphsWingRemoved);
2113
2114     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);
2115
2116     DrawCentrality("eta_rms_centrality", nHists, graphs[2], 0, 0.8, "#sigma_{#eta} (fit)", graphs[2+16], (graphsWingRemoved) ? graphsWingRemoved[2] : 0, 0, 0, 0, 1);
2117
2118     DrawCentrality("chi2_1", nHists, graphs[6], 0.5, 2.5, "#chi^{2}/ndf (full region)");
2119     DrawCentrality("chi2_2", nHists, graphs[7], 0.5, 2.5, "#chi^{2}/ndf (peak region)");
2120     
2121     DrawCentrality("sigma_phi", nHists, graphs[8], 0.1, 0.5, "#sigma_{#varphi} (rad.)", graphs[8+16], (graphsWingRemoved) ? graphsWingRemoved[8] : 0, 0, 0, 0, 1);
2122     DrawCentrality("sigma_eta", nHists, graphs[9], 0.1, 0.5, "#sigma_{#eta}", graphs[9+16], (graphsWingRemoved) ? graphsWingRemoved[9] : 0, 0, 0, 0, 1);
2123   }
2124   
2125   DrawCentrality("kurtosisphi_centrality", 12, graphs[14], -2, 2, "Kurtosis #varphi", graphs[14+16], (graphsWingRemoved) ? graphsWingRemoved[14] : 0, 0, 0, 0, 2);
2126   DrawCentrality("kurtosiseta_centrality", 12, graphs[15], -2, 2, "Kurtosis #eta", graphs[15+16], (graphsWingRemoved) ? graphsWingRemoved[15] : 0, 0, 0, 0, 2);
2127 }
2128 */
2129
2130 void MCComparison(const char* fileNameData, const char* fileNameWingRemoved, const char* fileNameHijing, const char* fileNameAMPT, Int_t offset = 0)
2131 {
2132   Int_t nHists = 12; //NHists;
2133
2134   ReadGraphs(fileNameWingRemoved);
2135   CalculateRMSSigma();
2136   TGraphErrors*** graphsWingRemoved = graphs;
2137
2138   ReadGraphs(fileNameHijing);
2139   CalculateRMSSigma();
2140   TGraphErrors*** graphs1 = graphs;
2141
2142   CreateGraphStructure();
2143   TGraphErrors*** graphs2 = graphs;
2144   if (fileNameAMPT)
2145   {
2146     ReadGraphs(fileNameAMPT);
2147     CalculateRMSSigma();
2148     graphs2 = graphs;
2149   }
2150
2151   ReadGraphs(fileNameData);
2152   CalculateRMSSigma();
2153
2154   if (0)
2155   {
2156     Int_t graphList[] = { 1, 2, 8, 9, 14, 15 };
2157     for (Int_t i=0; i<6; i++)
2158     {
2159       graphs[graphList[i]+offset][5] = new TGraphErrors;
2160       graphs[graphList[i]+offset][10] = new TGraphErrors;
2161       graphs1[graphList[i]+offset][5] = new TGraphErrors;
2162       graphs1[graphList[i]+offset][10] = new TGraphErrors;
2163     }
2164   }
2165   
2166   DrawCentrality("phi_rms_centrality_mc", nHists, graphs[1+offset], 0.2, 0.8, Form("#sigma_{#Delta#varphi} (%s) (rad.)", fitLabel), graphsWingRemoved[1+offset], 0, graphs1[1+offset], 0, graphs2[1+offset], 1);
2167   
2168 //   return;
2169   
2170   DrawCentrality("eta_rms_centrality_mc", nHists, graphs[2+offset], 0.2, 0.8 + 0.4 / 16 * offset, Form("#sigma_{#Delta#eta} (%s)", fitLabel), graphsWingRemoved[2+offset], 0, graphs1[2+offset], 0, graphs2[2+offset], 1);
2171   
2172   DrawCentrality("sigma_phi_mc", nHists, graphs[8+offset], 0.2, 0.6, "#sigma_{#Delta#varphi} (rad.)", graphsWingRemoved[8+offset], 0, graphs1[8+offset], 0, graphs2[8+offset], 1);
2173   DrawCentrality("sigma_eta_mc", nHists, graphs[9+offset], 0.2, 0.6, "#sigma_{#Delta#eta}", graphsWingRemoved[9+offset], 0, graphs1[9+offset], 0, graphs2[9+offset], 1);
2174
2175   DrawCentrality("kurtosisphi_centrality_mc", nHists, graphs[14+offset], -1.5, 2.5, "Kurtosis #Delta#varphi", graphsWingRemoved[14+offset], 0, graphs1[14+offset], 0, graphs2[14+offset], 2);
2176   DrawCentrality("kurtosiseta_centrality_mc", nHists, graphs[15+offset], -1.5, 2.5, "Kurtosis #Delta#eta", graphsWingRemoved[15+offset], 0, graphs1[15+offset], 0, graphs2[15+offset], 2);
2177
2178 /*  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);
2179   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);
2180
2181   DrawCentrality("sigma_phi", nHists, graphs[8], 0.1, 0.5, "#sigma_{#varphi} (rad.)", graphs[8+16], graphsWingRemoved[8], graphs1[8], 0, graphs2[8], 1);
2182   DrawCentrality("sigma_eta", nHists, graphs[9], 0.1, 0.5, "#sigma_{#eta}", graphs[9+16], graphsWingRemoved[9], graphs1[9], 0, graphs2[9], 1);
2183
2184   DrawCentrality("kurtosisphi_centrality_mc", nHists, graphs[14], -2, 2, "Kurtosis #varphi", graphs[14+16], graphsWingRemoved[14], graphs1[14], 0, graphs2[14], 2);
2185   DrawCentrality("kurtosiseta_centrality_mc", nHists, graphs[15], -2, 2, "Kurtosis #eta", graphs[15+16], graphsWingRemoved[15], graphs1[15], 0, graphs2[15], 2);*/
2186 }
2187
2188 void ShowWingEffect(const char* fileNameData, const char* fileNameWingRemoved, Int_t offset = 0)
2189 {
2190   Int_t nHists = 24; //NHists;
2191
2192   ReadGraphs(fileNameWingRemoved);
2193   CalculateRMSSigma();
2194   TGraphErrors*** graphs1 = graphs;
2195
2196   ReadGraphs(fileNameData);
2197   CalculateRMSSigma();
2198   
2199   DrawCentrality("phi_rms_centrality_mc", nHists, graphs[1+offset], 0, 0.8, "#sigma_{#varphi} (fit) (rad.)", graphs1[1+offset]);
2200   DrawCentrality("eta_rms_centrality_mc", nHists, graphs[2+offset], 0, 0.8, "#sigma_{#eta} (fit)", graphs1[2+offset]);
2201
2202   DrawCentrality("sigma_phi", nHists, graphs[8+offset], 0, 0.8, "#sigma_{#varphi} (rad.)", graphs1[8+offset]);
2203   DrawCentrality("sigma_eta", nHists, graphs[9+offset], 0, 0.8, "#sigma_{#eta}", graphs1[9+offset]);
2204
2205   DrawCentrality("kurtosisphi_centrality_mc", nHists, graphs[14+offset], -2, 4, "Kurtosis #varphi", graphs1[14+offset]);
2206   DrawCentrality("kurtosiseta_centrality_mc", nHists, graphs[15+offset], -2, 4, "Kurtosis #eta", graphs1[15+offset]);
2207 }
2208
2209 void ShowFitEffect(const char* fileNameData)
2210 {
2211   Int_t nHists = 12; //NHists;
2212
2213   ReadGraphs(fileNameData);
2214   CalculateRMSSigma();
2215   
2216   DrawCentrality("phi_rms_centrality_mc", nHists, graphs[1], 0, 0.8, "#sigma_{#varphi} (fit) (rad.)", graphs[1+16]);
2217   DrawCentrality("eta_rms_centrality_mc", nHists, graphs[2], 0, 0.8, "#sigma_{#eta} (fit)", graphs[2+16]);
2218
2219   DrawCentrality("sigma_phi", nHists, graphs[8], 0, 0.8, "#sigma_{#varphi} (rad.)", graphs[8+16]);
2220   DrawCentrality("sigma_eta", nHists, graphs[9], 0, 0.8, "#sigma_{#eta}", graphs[9+16]);
2221
2222   DrawCentrality("kurtosisphi_centrality_mc", nHists, graphs[14], -2, 4, "Kurtosis #varphi", graphs[14+16]);
2223   DrawCentrality("kurtosiseta_centrality_mc", nHists, graphs[15], -2, 4, "Kurtosis #eta", graphs[15+16]);
2224 }
2225
2226 void CompareSigmas(const char* fileNameData)
2227 {
2228   // compare sigma from fit with sigma from direct calculation but taking for the first the limited range of 0.8 into account
2229   
2230   Int_t nHists = 12; //NHists;
2231
2232   ReadGraphs(fileNameData);
2233   CalculateRMSSigma();
2234   TGraphErrors*** graphsOriginal = graphs;
2235   
2236   ReadGraphs(fileNameData);
2237   
2238   // rms
2239   for (Int_t histId = 0; histId < nHists; histId++)
2240   {
2241     if (!graphs[0][histId])
2242       continue;
2243     
2244     for (Int_t i=0; i<graphs[0][histId]->GetN(); i++)
2245     {
2246       // phi
2247       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());
2248       func->SetParameter(0, 0);
2249       for (Int_t k=0; k<6; k++)
2250         func->SetParameter(k+1, graphs[k][histId]->GetY()[i]);
2251       func->SetParameter(7, 0.8); // project Limit
2252       
2253       graphs[8][histId]->GetY()[i] = TMath::Sqrt(func->CentralMoment(2, -0.87, 0.87));
2254       
2255       // eta
2256       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);
2257       func->SetParameter(0, 0);
2258       for (Int_t k=0; k<6; k++)
2259         func->SetParameter(k+1, graphs[k][histId]->GetY()[i]);
2260       func->SetParameter(7, 0.87); // project Limit
2261       
2262       graphs[9][histId]->GetY()[i] = TMath::Sqrt(func->CentralMoment(2, -0.8, 0.8));
2263     }
2264   }
2265   
2266   DrawCentrality("phi_comparison1", nHists, graphsOriginal[1], 0.2, 0.8, "#sigma_{#Delta#varphi} (fit)", 0, 0, graphs[8]);
2267   DrawCentrality("phi_comparison2", nHists, graphsOriginal[8], 0.2, 0.5, "#sigma_{#Delta#varphi} (fit)", 0, 0, graphs[8]);
2268   DrawCentrality("eta_comparison1", nHists, graphsOriginal[2], 0.2, 0.8, "#sigma_{#Delta#eta} (fit)", 0, 0, graphs[9]);
2269   DrawCentrality("eta_comparison2", nHists, graphsOriginal[9], 0.2, 0.5, "#sigma_{#Delta#eta} (fit)", 0, 0, graphs[9]);
2270 }  
2271
2272 Float_t** ExtractSystematics(const char* baseFile, const char* systFile, Int_t offset)
2273 {
2274   ReadGraphs(baseFile);
2275   CalculateRMSSigma();
2276  
2277   TGraphErrors*** graphsBase = graphs;
2278   if (systFile)
2279   {
2280     ReadGraphs(systFile);
2281     CalculateRMSSigma();
2282   }
2283   
2284   // calculate syst unc for these graphs
2285   const Int_t NGraphList = 6;
2286   Int_t graphList[] = { 1, 2, 8, 9, 14, 15 };
2287
2288   Float_t** results = new Float_t*[NGraphList];
2289
2290   for (Int_t i=0; i<NGraphList; i++)
2291   {
2292     results[i] = new Float_t[2];
2293
2294     Bool_t kurtosis = kFALSE;
2295     if (graphList[i] == 14 || graphList[i] == 15)
2296       kurtosis = kTRUE;
2297     
2298     TH1* hist = 0;
2299     if (!kurtosis)
2300       hist = new TH1F(Form("hist_%d", i), "", 50, 0.5, 1.5);  
2301     else
2302       hist = new TH1F(Form("hist_%d", i), "", 50, -1, 1);  
2303   
2304     TCanvas* c = new TCanvas(Form("%s_%d_%d", systFile, i, 0), Form("%s_%d_%d", systFile, i, 0), 1000, 1000);
2305     c->Divide(4, 4);
2306
2307     Int_t count = 1;
2308     for (Int_t j=0; j<NHists; j++)
2309     {
2310       if (SkipGraph(j))
2311         continue;
2312       
2313       // for kurtosis we show only up to 12
2314       if (kurtosis)
2315         if (j >= 12)
2316           continue;
2317       
2318       TGraphErrors* graph1 = graphsBase[graphList[i]+offset][j];
2319       TGraphErrors* graph2 = (systFile) ? graphs[graphList[i]+offset][j] : graphsBase[graphList[i]+16][j];
2320       
2321       if (graph1->GetN() == 0)
2322         continue;
2323       
2324       graph1->Sort();
2325       graph2->Sort();
2326       
2327       c->cd(count++);
2328       gPad->SetGridx();
2329       gPad->SetGridy();
2330       graph1->SetMarkerStyle(24);
2331       graph1->DrawClone("AP");
2332       graph2->SetLineColor(2);
2333       graph2->SetMarkerColor(2);
2334       graph2->SetMarkerStyle(25);
2335       graph2->DrawClone("LSAME");
2336       
2337       c->cd(count++);
2338       gPad->SetGridx();
2339       gPad->SetGridy();
2340       
2341       // reset errors on graph2 
2342 /*      for (Int_t k=0; k<graph2->GetN(); k++)
2343         graph2->GetEY()[k] = 0;*/
2344       
2345       if (!kurtosis)
2346       {
2347         DivideGraphs(graph1, graph2);
2348         ((TGraphErrors*) graph1->DrawClone("AP"))->GetYaxis()->SetRangeUser(0.8, 1.2);
2349       }
2350       else
2351       {
2352         SubtractGraphs(graph1, graph2);
2353         ((TGraphErrors*) graph1->DrawClone("AP"));
2354       } 
2355       
2356       for (Int_t k=0; k<graph1->GetN(); k++)
2357       {
2358 //      hist->Fill(graph1->GetY()[k], 1.0 / (graph1->GetEY()[k] / graph1->GetY()[k]));
2359         if (kurtosis)
2360           hist->Fill(TMath::Abs(graph1->GetY()[k]));
2361         else
2362           hist->Fill(graph1->GetY()[k]);
2363       }
2364       
2365       if (count == 37)
2366         break;
2367     }
2368     
2369     c->SaveAs(Form("syst/%s_%d.eps", systFile, i));
2370     
2371     new TCanvas;
2372     hist->Draw();
2373     hist->Sumw2();
2374     
2375     Float_t mean = -1;
2376     Float_t sigma = -1;
2377     if (0)
2378     {
2379       hist->Fit("gaus", "Q");
2380       mean = hist->GetFunction("gaus")->GetParameter(1);
2381       sigma = hist->GetFunction("gaus")->GetParameter(2);
2382       
2383       Printf("%d: %.3f +- %.3f %.3f | %.3f +- %.3f", i, mean, sigma, sigma / TMath::Sqrt(hist->GetEntries()), hist->GetMean(), hist->GetMean(11));
2384     }
2385     else
2386     {
2387       mean = hist->GetMean(1);
2388       sigma = hist->GetMean(11);
2389
2390       Printf("%d: %.3f +- %.3f", i, hist->GetMean(), hist->GetMean(11));
2391     }
2392     
2393     if (!kurtosis)
2394     {
2395       mean -= 1;
2396       mean *= 100;
2397       sigma *= 100;
2398     }
2399     results[i][0] = mean;
2400     results[i][1] = sigma / TMath::Sqrt(hist->GetEntries());
2401   }
2402   
2403   return results;
2404 }
2405
2406 void GetProjections(TH2* hist, TH1** projPhi, TH1** projEta, Int_t k)
2407 {
2408   Float_t projectLimit = 0.8;
2409
2410   *projPhi = hist->ProjectionX(Form("%s_proj1%d", hist->GetName(), k), hist->GetYaxis()->FindBin(-projectLimit+0.01), hist->GetYaxis()->FindBin(projectLimit-0.01));
2411     
2412   *projEta = hist->ProjectionY(Form("%s_proj2_%d", hist->GetName(), k), hist->GetXaxis()->FindBin(-projectLimit+0.01), hist->GetXaxis()->FindBin(projectLimit-0.01));
2413 }
2414
2415 // uncorrected ones
2416 const char* systBaseFile = "dphi_corr_2d_120429.root";
2417 const char* systTrackCuts1 = "dphi_corr_120425_hybrid.root";
2418 const char* systTrackCuts2 = "dphi_corr_120430_raa.root";
2419 const char* systVertex = "dphi_corr_2d_120425_vertex.root";
2420 const char* systResonances = "dphi_corr_120502_resonances.root";
2421 const char* systTTR = "dphi_corr_2d_120430_widettr.root";
2422
2423 // corrected ones
2424 const char* systBaseFileCorrected = "dphi_corr_2d_120508.root";
2425 const char* systTrackCuts1Corrected = "dphi_corr_120507_hybrid.root";
2426 const char* systTrackCuts2Corrected = "dphi_corr_120507_raa.root";
2427 const char* systWingRemovedCorrected = "dphi_corr_2d_120508_wingremoved.root";
2428
2429 void ExtractSystematicsProjections(Int_t mode, Float_t rangeBegin = -0.4, Float_t rangeEnd = 0.4, Bool_t draw = kFALSE)
2430 {
2431   if (!draw)
2432     gROOT->SetBatch(kTRUE);
2433   
2434   Int_t maxLeadingPt = 4;
2435   Int_t maxAssocPt = 6;
2436
2437   Int_t NEffects = 3;
2438   
2439   TFile* file1 = 0;
2440   const char** systFiles = 0;
2441   if (mode == 2)
2442   {
2443     NEffects = 5;
2444     file1 = TFile::Open(systBaseFile);
2445     static const char* systFilesTmp[5] = { systVertex, systResonances, systTTR, systTrackCuts1Corrected, systTrackCuts2Corrected };
2446     systFiles = systFilesTmp;
2447   }
2448   else if (mode == 0)
2449   {
2450     file1 = TFile::Open(systBaseFile);
2451     static const char* systFilesTmp[3] = { systVertex, systResonances, systTTR };
2452     systFiles = (const char**) systFilesTmp;
2453   }
2454   else if (mode == 1)
2455   {
2456     NEffects = 3;
2457     file1 = TFile::Open(systBaseFileCorrected);
2458     static const char* systFilesTmp[3] = { systTrackCuts1Corrected, systTrackCuts2Corrected, systWingRemovedCorrected };
2459     systFiles = systFilesTmp;
2460   }
2461   else if (mode == 3)
2462   {
2463     NEffects = 3;
2464     file1 = TFile::Open(systBaseFileCorrected);
2465     static const char* systFilesTmp[3] = { systBaseFileCorrected, systBaseFileCorrected, systBaseFileCorrected };
2466     systFiles = systFilesTmp;
2467   }
2468   
2469   TH2* uncertainty[4];
2470   for (Int_t i=0; i<4; i++)
2471     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);
2472   
2473   Int_t count = 0;
2474   Int_t ptBin = 0;
2475   for (Int_t i=0; i<maxLeadingPt-1; i++)
2476   {
2477     for (Int_t j=1; j<maxAssocPt-1; j++)
2478     {
2479       if (i == 0 && j == 2)
2480         continue;
2481       if (i == 1 && j == 3)
2482         continue;
2483       if (i == 2 && j == 4)
2484         continue;
2485       
2486       ptBin++;
2487       for (Int_t histId = 0; histId < NHists; histId++)
2488       {
2489         TH2* hist1 = (TH2*) file1->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
2490         if (!hist1)
2491           continue;
2492         
2493         if (hist1->GetEntries() < 1e4)
2494         {
2495           Printf("Only %f entries. Skipping...", hist1->GetEntries());
2496           continue;
2497         }
2498         
2499         Printf("%d %d %d %s", i, j, histId, hist1->GetTitle());
2500         
2501         TH1* proj1[2];
2502         TH1* proj2[10][2];
2503         
2504         SubtractEtaGap(hist1, kEtaLimit, kOuterLimit, kFALSE);
2505         GetProjections(hist1, &proj1[0], &proj1[1], count++);
2506
2507         TCanvas* c = new TCanvas(Form("c_%d_%d_%d", mode, ptBin, histId), Form("c_%d_%d_%d", mode, ptBin, histId), 1200, 1000);
2508         c->Divide(3, 2);
2509           
2510         for (Int_t k=0; k<2; k++)
2511         {
2512           c->cd(1+k*3);
2513           proj1[k]->SetStats(0);
2514           proj1[k]->GetXaxis()->SetRangeUser(-1.5, 1.5);
2515           proj1[k]->GetYaxis()->SetRangeUser(proj1[k]->GetMinimum() * 1.1, proj1[k]->GetMaximum() * 1.4);
2516           proj1[k]->DrawCopy();
2517         }
2518         
2519         Double_t maxAverage[2] = { 0, 0 };
2520         Double_t maxDev[2] = { 0, 0 };
2521   
2522         for (Int_t n=0; n<NEffects; n++)
2523         {
2524           if (histId == 2 && mode == 0 && n > 0)
2525             continue;
2526           if (histId == 2 && mode == 1 && n > 0)
2527             continue;
2528 //        Printf("%d %s", n, systFiles[n]);
2529           TFile* file2 = TFile::Open(systFiles[n]);
2530           hist1 = (TH2*) file2->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
2531           if (!hist1)
2532             continue;
2533           
2534           if (mode == 3 && n == 0) // change eta limits
2535             SubtractEtaGap(hist1, kEtaLimit-0.2, kOuterLimit, kFALSE);
2536           if (mode == 3 && n == 1) // change eta limits
2537             SubtractEtaGap(hist1, kEtaLimit+0.2, kOuterLimit, kFALSE);
2538           if (mode == 3 && n == 2) // change eta limits
2539             SubtractEtaGap(hist1, kEtaLimit, kOuterLimit-0.2, kFALSE);
2540           else
2541             SubtractEtaGap(hist1, kEtaLimit, kOuterLimit, kFALSE);
2542           GetProjections(hist1, &proj2[n][0], &proj2[n][1], count++);
2543           
2544           for (Int_t k=0; k<2; k++)
2545           {
2546             c->cd(1+k*3);
2547             proj2[n][k]->SetLineColor(n+2);
2548             /* TH1* copy = */ proj2[n][k]->DrawCopy("SAME");
2549             
2550             c->cd(2+k*3);
2551             gPad->SetGridx();
2552             gPad->SetGridy();
2553             TH1* ratio = (TH1*) proj2[n][k]->Clone(Form("%s_ratio", proj2[n][k]->GetName()));
2554             ratio->SetStats(0);
2555             ratio->Divide(proj1[k]);
2556             ratio->GetXaxis()->SetRangeUser(-1.5, 1.5);
2557             ratio->Fit("pol0", "0Q", "", -1, 1);
2558             Double_t average = ratio->GetFunction("pol0")->GetParameter(0);
2559 //          Printf("Average is %f", average);
2560             maxAverage[k] = TMath::Max(maxAverage[k], TMath::Abs(average - 1));
2561 //          ratio->Scale(1.0 / average);
2562 //          copy->Scale(1.0 / average);
2563             ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
2564             ratio->DrawCopy((n == 0) ? "" : "SAME");
2565             
2566             if (0)
2567             {
2568               Double_t sum = 0;
2569               Double_t sumCount = 0;
2570               for (Int_t bin = ratio->FindBin(rangeBegin); bin <= ratio->FindBin(rangeEnd); bin++)
2571               {
2572                 sum += TMath::Abs(ratio->GetBinContent(bin) - 1);
2573                 sumCount++;
2574               }
2575               maxDev[k] = TMath::Max(maxDev[k], sum / sumCount);
2576             }
2577
2578             c->cd(3+k*3);
2579             gPad->SetGridx();
2580             gPad->SetGridy();
2581             TH1* diff = (TH1*) proj2[n][k]->Clone(Form("%s_diff", proj2[n][k]->GetName()));
2582             diff->Scale(1.0 / average);
2583             diff->Add(proj1[k], -1);
2584             diff->SetStats(0);
2585             diff->GetXaxis()->SetRangeUser(-1.5, 1.5);
2586             diff->GetYaxis()->SetRangeUser(-0.05, 0.05);
2587             diff->DrawCopy((n == 0) ? "" : "SAME");
2588             
2589             if (0)
2590             {
2591               for (Int_t bin = diff->FindBin(rangeBegin); bin <= diff->FindBin(rangeEnd); bin++)
2592                 maxDev[k] = TMath::Max(maxDev[k], TMath::Abs(diff->GetBinContent(bin)));
2593             }
2594             else
2595             {
2596               Double_t sum = 0;
2597               Double_t sumCount = 0;
2598               for (Int_t bin = diff->FindBin(rangeBegin); bin <= diff->FindBin(rangeEnd); bin++)
2599               {
2600                 sum += TMath::Abs(diff->GetBinContent(bin));
2601                 sumCount++;
2602               }
2603               maxDev[k] = TMath::Max(maxDev[k], sum / sumCount);
2604             }
2605           }
2606           
2607           delete file2;
2608         }
2609         
2610         // normalize to peak strength
2611 /*      for (Int_t k=0; k<2; k++)
2612           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);*/
2613         
2614         Int_t centralityAxisMapping[] = { 0, 4, 3, 5, 1, 2 };
2615         uncertainty[0]->SetBinContent(ptBin, centralityAxisMapping[histId]+1, maxAverage[0]);
2616         uncertainty[1]->SetBinContent(ptBin, centralityAxisMapping[histId]+1, maxAverage[1]);
2617         uncertainty[2]->SetBinContent(ptBin, centralityAxisMapping[histId]+1, maxDev[0]);
2618         uncertainty[3]->SetBinContent(ptBin, centralityAxisMapping[histId]+1, maxDev[1]);
2619         
2620         c->SaveAs(Form("syst_corr/%s.eps", c->GetName()));
2621       }
2622       
2623       if (draw)
2624         break;
2625     }
2626     if (draw)
2627       break;
2628   }
2629   
2630   TFile::Open("uncertainty.root", "RECREATE");
2631   for (Int_t i=0; i<4; i++)
2632     uncertainty[i]->Write();
2633   gFile->Close();
2634 }
2635
2636 void DrawSystematicsProjections()
2637 {
2638   TH2* uncertainty[4];
2639   
2640   TFile::Open("uncertainty.root");
2641   for (Int_t i=0; i<4; i++)
2642     uncertainty[i] = (TH2*) gFile->Get(Form("uncertainty_%d", i));
2643   
2644   TCanvas* cResult = new TCanvas("cResult", "cResult", 1000, 1000);
2645   cResult->Divide(2, 2);
2646   for (Int_t i=0; i<4; i++)
2647   {
2648     cResult->cd(i+1);
2649     gPad->SetRightMargin(0.15);
2650     uncertainty[i]->SetStats(0);
2651     uncertainty[i]->Draw("colz");
2652   }
2653
2654   TCanvas* cResult2 = new TCanvas("cResult2", "cResult2", 1000, 1000);
2655   cResult2->Divide(2, 2);
2656   for (Int_t i=0; i<4; i++)
2657   {
2658     cResult2->cd(i+1);
2659     gPad->SetRightMargin(0.15);
2660     Int_t color = 1;
2661     for (Int_t j=1; j<=uncertainty[i]->GetNbinsX(); j++)
2662     {
2663       TH1* proj = uncertainty[i]->ProjectionY(Form("p_%d_%d", i, j), j, j);
2664       if (proj->Integral() > 0)
2665       {
2666         proj->SetLineColor(color++);
2667         proj->SetStats(0);
2668         proj->GetYaxis()->SetRangeUser(0, 0.5);
2669         proj->Draw((color == 2) ? "" : "SAME");
2670         
2671         Printf("%d %d: %.3f %.3f", i, color, proj->GetBinContent(1), proj->Integral(2, 6) / 5);
2672       }
2673     }
2674   }
2675 }
2676
2677 void BuildSystematicFiles()
2678 {
2679   Printf("Hope you know what you are doing... 5 seconds to abort...");
2680   gSystem->Sleep(5000);
2681
2682   AnalyzeDeltaPhiEtaGap2D(systBaseFile, "syst_base.root");
2683
2684   kEtaLimit = 0.8;
2685   AnalyzeDeltaPhiEtaGap2D(systBaseFile, "syst_eta08.root");
2686   
2687   kEtaLimit = 1.2;
2688   AnalyzeDeltaPhiEtaGap2D(systBaseFile, "syst_eta12.root");
2689   kEtaLimit = 1.0;
2690   
2691   kOuterLimit = 1.39;
2692   AnalyzeDeltaPhiEtaGap2D(systBaseFile, "syst_outer14.root");
2693   kOuterLimit = 1.59;
2694   
2695   AnalyzeDeltaPhiEtaGap2D(systTrackCuts1, "syst_hybrid.root");
2696   AnalyzeDeltaPhiEtaGap2D(systTrackCuts2, "syst_raa.root");
2697   AnalyzeDeltaPhiEtaGap2D(systVertex, "syst_vertex.root");
2698   AnalyzeDeltaPhiEtaGap2D(systResonances, "syst_resonances.root");
2699   AnalyzeDeltaPhiEtaGap2D(systTTR, "syst_widettr.root");
2700 }
2701
2702 void ExtractSystematicsAll(Int_t offset = 0)
2703 {
2704   gROOT->SetBatch(kTRUE);
2705   
2706   const Int_t NEffects = 9;
2707   
2708 //   const char* defaultFile = "graphs_120425.root";
2709 //   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" };
2710   const char* defaultFile = "syst_base.root";
2711   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" };
2712   const char* effectIdString[] = { "1a", "1b", "1c", "2a", "2b", "3", "4", "5", "6" };
2713
2714   Float_t** results[NEffects+3];
2715   
2716   for (Int_t i=0; i<NEffects; i++)
2717   {
2718     results[i] = ExtractSystematics(defaultFile, systFiles[i], offset);
2719   }
2720   
2721   const Int_t NParameters = 6;
2722   const char* names[] = { "$\\sigma_{\\Dphi}$ (fit)", "$\\sigma_{\\Deta}$ (fit)", "$\\sigma_{\\Dphi}$", "$\\sigma_{\\Deta}$", "Kurtosis $\\Dphi$", "Kurtosis $\\Deta$" };
2723   
2724   for (Int_t j=0; j<NParameters; j++)
2725   {
2726     printf("%s \t & ", names[j]);
2727     for (Int_t i=0; i<NEffects; i++)
2728     {
2729       if (j < 4)
2730         printf("$%.1f\\%% \\pm %.1f\\%%$ \t & ", results[i][j][0], results[i][j][1]);
2731       else
2732         printf("$%.2f \\pm %.2f$ \t & ", results[i][j][0], results[i][j][1]);
2733     }
2734     Printf("\\\\");
2735   }
2736
2737   // put together 0-2 (into NEffects)
2738   results[NEffects] = new Float_t*[NParameters];
2739   for (Int_t j=0; j<NParameters; j++)
2740   {
2741     results[NEffects][j] = new Float_t[2];
2742     Float_t mean = 0;
2743     Float_t sigma = 0;
2744     
2745     printf("%s \t:", names[j]);
2746     for (Int_t i=0; i<=2; i++)
2747     {
2748       mean = TMath::Max(mean, TMath::Abs(results[i][j][0]));
2749       printf("%.1f%% ", results[i][j][0]);
2750     }
2751     printf("--> %.1f%% \t", mean);
2752     results[NEffects][j][0] = mean;
2753       
2754     for (Int_t i=0; i<=2; i++)
2755     {
2756       sigma = TMath::Max(sigma, TMath::Abs(results[i][j][1]));
2757       printf("%.1f%% ", results[i][j][1]);
2758     }
2759     Printf("--> %.1f%% \t", sigma);
2760     results[NEffects][j][1] = sigma;
2761   }
2762
2763   // put together 3-4 (into NEffects+1)
2764   results[NEffects+1] = new Float_t*[NParameters];
2765   for (Int_t j=0; j<NParameters; j++)
2766   {
2767     results[NEffects+1][j] = new Float_t[2];
2768     
2769     Float_t mean = 0;
2770     Float_t sigma = 0;
2771     
2772     printf("%s \t:", names[j]);
2773     for (Int_t i=3; i<=4; i++)
2774     {
2775       mean = TMath::Max(mean, TMath::Abs(results[i][j][0]));
2776       printf("%.1f%% ", results[i][j][0]);
2777     }
2778     printf("--> %.1f%% \t", mean);
2779     results[NEffects+1][j][0] = mean;
2780       
2781     for (Int_t i=3; i<=4; i++)
2782     {
2783       sigma = TMath::Max(sigma, TMath::Abs(results[i][j][1]));
2784       printf("%.1f%% ", results[i][j][1]);
2785     }
2786     Printf("--> %.1f%% \t", sigma);
2787     results[NEffects+1][j][1] = sigma;
2788   }
2789   
2790   // combine in quadrature (only mean)
2791   results[NEffects+2] = new Float_t*[NParameters];
2792   for (Int_t j=0; j<NParameters; j++)
2793   {
2794     results[NEffects+2][j] = new Float_t[2];
2795     Float_t mean = 0;
2796   
2797     printf("%s \t:", names[j]);
2798     for (Int_t i=5; i<NEffects+2; i++)
2799     {
2800       mean += results[i][j][0] * results[i][j][0];
2801       printf("%.1f%% ", results[i][j][0]);
2802     }
2803     mean = TMath::Sqrt(mean);
2804     Printf("--> %.1f%% \t", mean);
2805     results[NEffects+2][j][0] = mean;
2806   }
2807     
2808   for (Int_t j=0; j<NParameters; j++)
2809   {
2810     printf("%s \t & ", names[j]);
2811     for (Int_t i=0; i<NEffects+3; i++)
2812     {
2813       if (i == NEffects || i == NEffects+1)
2814         continue;
2815       if (j < 4)
2816         printf("$%.1f\\%%$ \t ", results[i][j][0]);
2817       else
2818         printf("$%.2f$ \t ", results[i][j][0]);
2819       if (i < NEffects+2)
2820         printf("& ");
2821     }
2822     Printf("\\\\");
2823   }
2824   
2825   // generate LaTeX code
2826   for (Int_t i=0; i<NEffects; i++)
2827   {
2828     for (Int_t j=0; j<NParameters; j++)
2829     {
2830       Printf("\\bfig[ht]");
2831       Printf("  \\includegraphics[width=\\linewidth]{syst/%s_%d.eps}", systFiles[i], j);
2832       Printf("  \\caption{Systematic effect (%s) - Parameter %s}", effectIdString[i], names[j]);
2833       Printf("\\efig");
2834     }
2835     Printf("\\clearpage");
2836   }
2837 }
2838
2839 void CompareEtaPhi(const char* fileName, const char* fileNameWingRemoved)
2840 {
2841   Int_t offset = 0;
2842   
2843   TGraphErrors*** graphsWingRemoved = 0;
2844   if (fileNameWingRemoved)
2845   {
2846     ReadGraphs(fileNameWingRemoved);
2847     CalculateRMSSigma();
2848     graphsWingRemoved = graphs;
2849   }
2850   
2851   ReadGraphs(fileName);
2852   
2853   Int_t nHists = 12; //NHists;
2854   skipGraphList[0] = 5;
2855   skipGraphList[1] = 6;
2856   skipGraphList[2] = 10;
2857   
2858   CalculateRMSSigma();
2859   
2860   drawLogo = 1;
2861   
2862   MCLabel = "#sigma_{#Delta#varphi}     #sigma_{#Delta#eta}";
2863
2864   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);
2865 //   DrawCentrality("rms_centrality_nosyst", nHists, graphs[1], 0, 0.8, "#sigma_{#varphi} (fit) (rad.) / #sigma_{#eta}", 0, 0,  graphs[2], 0);
2866 }
2867
2868 void DrawExamples(const char* histFileName, const char* graphFileName, Int_t i = 0, Int_t j = 1, Bool_t drawFunc = kTRUE)
2869 {
2870   Float_t etaLimit = kEtaLimit;
2871   Float_t outerLimit = kOuterLimit;
2872   Float_t projectLimit = 0.8;
2873
2874   ReadGraphs(graphFileName);
2875
2876   TFile::Open(histFileName);
2877   
2878   Int_t exColors[] = { 1, 2, 4, 3, 5, 6 };
2879   
2880   Int_t graphID = i * (6 - 1) + j - 1;
2881
2882   TCanvas* c = new TCanvas(Form("c_%d", i), Form("c_%d", i), 1200, 600);
2883   c->Divide(2, 1);
2884   Int_t nHists = 3;
2885   for (Int_t histId = 0; histId < nHists; histId++)
2886   {
2887     TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
2888     if (!hist)
2889       continue;
2890
2891     if (graphs[0][graphID]->GetN() < histId)
2892     {
2893       Printf("ERROR: Pos in graph not found: %d %d", i, histId);
2894       continue;
2895     }
2896     
2897     SubtractEtaGap(hist, etaLimit, outerLimit, kFALSE, kFALSE);
2898   
2899     c->cd(1);
2900     TH1* proj = hist->ProjectionX(Form("%s_proj1", hist->GetName()), hist->GetYaxis()->FindBin(-projectLimit+0.01), hist->GetYaxis()->FindBin(projectLimit-0.01));
2901     proj->SetLineColor(exColors[histId]);
2902     proj->SetStats(0);
2903     proj->GetXaxis()->SetRangeUser(-1.49, 1.49);
2904     proj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 1.2, proj->GetMaximum() * 1.5);
2905     proj->GetYaxis()->SetTitle(kProjYieldTitlePhi);
2906     TString label(proj->GetTitle());
2907     TObjArray* objArray = label.Tokenize("-");
2908     proj->SetTitle(Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()));
2909     proj->Draw((histId == 0) ? "" : "SAME");
2910   
2911     if (drawFunc)
2912     {
2913       // integral over y
2914       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());
2915       func->SetParameter(0, 0);
2916       for (Int_t k=0; k<6; k++)
2917         func->SetParameter(k+1, graphs[k][graphID]->GetY()[histId]);
2918       // scale by bin width (to compare with projection)
2919       func->SetParameter(1, func->GetParameter(1) / hist->GetYaxis()->GetBinWidth(1));
2920       func->SetLineColor(exColors[histId]);
2921       func->SetLineWidth(1);
2922       func->DrawCopy("SAME");
2923       
2924       // draw contributions
2925       Float_t scale = func->GetParameter(1);
2926       Float_t weighting = func->GetParameter(4);
2927       func->SetParameter(1, scale * weighting);
2928       func->SetParameter(4, 1);
2929       func->SetLineStyle(2);
2930       func->DrawCopy("SAME");
2931
2932       func->SetParameter(1, scale * (1.0-weighting));
2933       func->SetParameter(4, 0);
2934       func->SetLineStyle(3);
2935       func->DrawCopy("SAME");
2936     }
2937
2938     DrawLatex(0.15, 0.8 - 0.05 * histId, exColors[histId], labels[histId]);
2939
2940     c->cd(2);
2941     proj = hist->ProjectionY(Form("%s_proj2b", hist->GetName()), hist->GetXaxis()->FindBin(-projectLimit+0.01), hist->GetXaxis()->FindBin(projectLimit-0.01));
2942     proj->SetLineColor(exColors[histId]);
2943     proj->GetXaxis()->SetRangeUser(-1.49, 1.49);
2944     proj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 1.2, proj->GetMaximum() * 2);
2945     proj->GetYaxis()->SetTitle(kProjYieldTitleEta);
2946     proj->SetStats(0);
2947     proj->SetTitle(Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()));
2948     proj->Draw((histId == 0) ? "" : "SAME");
2949     
2950     if (0)
2951     {
2952       proj = hist->ProjectionY(Form("%s_proj2", hist->GetName()), hist->GetXaxis()->FindBin(-0.5 * TMath::Pi()), hist->GetXaxis()->FindBin(0.5 * TMath::Pi()));
2953       proj->SetLineColor(exColors[histId]);
2954       proj->GetXaxis()->SetRangeUser(-1.2, 1.2);
2955       proj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 1.2, proj->GetMaximum() * 2);
2956       proj->SetLineStyle(2);
2957       proj->Draw("SAME");
2958     }
2959
2960     if (drawFunc)
2961     {
2962       // integral over x
2963       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);
2964       func->SetParameter(0, 0);
2965       for (Int_t k=0; k<6; k++)
2966         func->SetParameter(k+1, graphs[k][graphID]->GetY()[histId]);
2967       // scale by bin width (to compare with projection)
2968       func->SetParameter(1, func->GetParameter(1) / hist->GetXaxis()->GetBinWidth(1));
2969       func->SetLineColor(exColors[histId]);
2970       func->SetLineWidth(1);
2971       func->DrawCopy("SAME");
2972
2973       // draw contributions
2974       Float_t scale = func->GetParameter(1);
2975       Float_t weighting = func->GetParameter(4);
2976       func->SetParameter(1, scale * weighting);
2977       func->SetParameter(4, 1);
2978       func->SetLineStyle(2);
2979       func->DrawCopy("SAME");
2980
2981       func->SetParameter(1, scale * (1.0-weighting));
2982       func->SetParameter(4, 0);
2983       func->SetLineStyle(3);
2984       func->DrawCopy("SAME");   
2985     }
2986
2987     DrawLatex(0.15, 0.8 - 0.05 * histId, exColors[histId], labels[histId]);
2988   }
2989 }
2990
2991 TH1* GetSystUnc(TH1* hist, Int_t i, Int_t j, Int_t histId, Int_t etaPhi)
2992 {
2993 //   Float_t uncertainty = (histId == 0) ? 0.1 : 0.08;
2994
2995   // track cuts
2996   //   2 2: 0.048 0.020
2997   //   2 3: 0.078 0.031
2998   //   2 4: 0.023 0.008
2999   //   2 5: 0.159 0.056
3000   //   2 6: 0.043 0.015
3001   //   2 7: 0.011 0.006
3002   //   3 2: 0.050 0.022
3003   //   3 3: 0.089 0.033
3004   //   3 4: 0.020 0.008
3005   //   3 5: 0.183 0.056
3006   //   3 6: 0.042 0.013
3007   //   3 7: 0.010 0.005
3008   
3009   // others
3010   //   2 2: 0.023 0.007
3011   //   2 3: 0.024 0.015
3012   //   2 4: 0.008 0.005
3013   //   2 5: 0.075 0.024
3014   //   2 6: 0.019 0.010
3015   //   2 7: 0.007 0.004
3016   //   3 2: 0.029 0.009
3017   //   3 3: 0.032 0.012
3018   //   3 4: 0.005 0.003
3019   //   3 5: 0.086 0.013
3020   //   3 6: 0.015 0.006
3021   //   3 7: 0.008 0.003
3022
3023   Float_t uncertainty = 1;
3024   if (i == 0 && j == 1)
3025     uncertainty = 0.021;
3026   if (i == 1 && j == 1)
3027     uncertainty = 0.032;
3028   if (i == 1 && j == 2)
3029     uncertainty = 0.008;
3030   if (i == 2 && j == 1)
3031     uncertainty = 0.056;
3032   if (i == 2 && j == 2)
3033     uncertainty = 0.014;
3034   if (i == 2 && j == 3)
3035     uncertainty = 0.006;
3036   
3037   if (histId == 0)
3038     uncertainty *= 2;
3039   
3040   TH1* systUnc = (TH1*) hist->Clone(Form("%s_syst", hist->GetName()));
3041   for (Int_t n=1; n<=systUnc->GetNbinsX(); n++)
3042     systUnc->SetBinError(n, uncertainty);
3043   
3044   systUnc->SetFillColor(hist->GetLineColor());
3045   systUnc->SetFillStyle((etaPhi == 0) ? 3004 : 3005);
3046   systUnc->SetMarkerStyle(0);
3047   systUnc->SetLineColor(0);
3048   
3049   return systUnc;
3050 }
3051
3052 Bool_t disableUncertainties = kTRUE;
3053
3054 void DrawExample(const char* histFileName, Int_t i, Int_t j, Int_t histId, TH1** projPhi, TH1** projEta, TH1** projPhiSyst, TH1** projEtaSyst, Bool_t Ratio = kFALSE)
3055 {
3056   Float_t etaLimit = kEtaLimit;
3057   Float_t outerLimit = kOuterLimit;
3058   Float_t projectLimit = 0.8;
3059   Int_t maxAssocPt = 6;
3060   Int_t graphID = i * (maxAssocPt - 1) + j - 1;
3061
3062   TFile::Open(histFileName);
3063   
3064   TCanvas* c = new TCanvas(Form("ex_%d_%d_%d", i, j, histId), Form("ex_%d_%d_%d", i, j, histId), 1800, 600);
3065   c->Divide(3, 1);
3066
3067   TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
3068   if (!hist)
3069     return;
3070   // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
3071   hist->Scale(1.0 / hist->GetYaxis()->GetBinWidth(1));
3072   hist->GetZaxis()->SetTitle(kCorrFuncTitle);
3073   hist->GetZaxis()->SetTitleOffset(1.9);
3074   
3075   if (graphID >= 15)
3076   {
3077     etaLimit = 0.5;
3078     outerLimit = 0.99;
3079   }
3080 //   hist->Rebin2D(2, 2); hist->Scale(0.25);
3081   
3082   TString label(hist->GetTitle());
3083   label.ReplaceAll(".00", " GeV/c");
3084   label.ReplaceAll(".0", " GeV/c");
3085   TObjArray* objArray = label.Tokenize("-");
3086   TPaveText* paveText = new TPaveText(0.52, 0.72, 0.95, 0.95, "BRNDC");
3087   paveText->SetTextSize(0.035);
3088   paveText->SetFillColor(0);
3089   paveText->SetShadowColor(0);
3090   paveText->AddText(objArray->At(0)->GetName());
3091   paveText->AddText(objArray->At(1)->GetName());
3092   if (objArray->GetEntries() == 4)
3093     paveText->AddText(Form("Pb-Pb 2.76 TeV %s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()));
3094   else
3095     paveText->AddText(Form("%s 2.76 TeV", objArray->At(2)->GetName()));
3096   paveText->AddText("|#eta| < 0.9");
3097   
3098   c->cd(1);
3099   gPad->SetLeftMargin(0.15);
3100   hist->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
3101   hist->GetXaxis()->SetTitleOffset(1.5);
3102   hist->GetYaxis()->SetTitleOffset(1.7);
3103   hist->SetStats(0);
3104   hist->SetTitle("a) Correlation");
3105   TH2* clone = (TH2*) hist->Clone(Form("%s_clone", hist->GetName()));
3106   clone->GetXaxis()->SetRangeUser(-TMath::Pi() / 2, TMath::Pi() / 2);
3107   clone->Draw("SURF1");
3108   paveText->Draw();
3109   
3110   SubtractEtaGap(hist, etaLimit, outerLimit, kFALSE, kFALSE);
3111   c->cd(2);
3112   gPad->SetLeftMargin(0.15);
3113   hist->SetTitle("b) #eta-gap subtracted");
3114   hist->GetXaxis()->SetRangeUser(-TMath::Pi() / 2, TMath::Pi() / 2);
3115   hist->DrawCopy("SURF1");
3116   if (!disableUncertainties)
3117     DrawALICELogo(kTRUE, 0.7, 0.7, 0.9, 0.9);
3118     
3119   c->cd(3);
3120   gPad->SetLeftMargin(0.13);
3121   TH1* proj = hist->ProjectionX(Form("%s_proj1", hist->GetName()), hist->GetYaxis()->FindBin(-projectLimit+0.01), hist->GetYaxis()->FindBin(projectLimit-0.01));
3122   TH1* proj2 = hist->ProjectionY(Form("%s_proj2b", hist->GetName()), hist->GetXaxis()->FindBin(-projectLimit+0.01), hist->GetXaxis()->FindBin(projectLimit-0.01));
3123   // normalization
3124   proj->Scale(hist->GetYaxis()->GetBinWidth(1));
3125   proj2->Scale(hist->GetXaxis()->GetBinWidth(1));
3126
3127   proj->SetStats(0);
3128   proj->SetTitle("c) Projections");
3129   proj->GetXaxis()->SetRangeUser(-1.49, 1.49);
3130   proj->GetYaxis()->SetTitle(kProjYieldTitlePhi);
3131   proj->GetYaxis()->SetTitleOffset(1.3);
3132   proj->GetXaxis()->SetTitleOffset(1);
3133   proj->GetYaxis()->SetRangeUser(proj->GetMinimum(), proj->GetMaximum() * 1.6);
3134   TH1* systUncPhi = GetSystUnc(proj, i, j, histId, 0);
3135   TH1* copy = systUncPhi->DrawCopy("E2 ][");
3136   copy->GetXaxis()->SetTitle(Form("%s , %s", proj->GetXaxis()->GetTitle(), proj2->GetXaxis()->GetTitle()));
3137   copy->GetYaxis()->SetTitle(kProjYieldTitlePhiOrEta);
3138     
3139   proj2->SetLineColor(2);
3140   proj2->SetStats(0);
3141   proj2->GetYaxis()->SetTitle(kProjYieldTitleEta);
3142   proj2->GetYaxis()->SetTitleOffset(1.2);
3143   proj2->GetXaxis()->SetTitleOffset(1);
3144   proj2->GetYaxis()->SetRangeUser(proj2->GetMinimum(), proj2->GetMaximum() * 1.6);
3145   TH1* systUncEta = GetSystUnc(proj2, i, j, histId, 1);
3146   systUncEta->Draw("E2 ][ SAME");
3147   proj->DrawCopy((disableUncertainties) ? "" : "SAME");
3148   proj2->DrawCopy("SAME");
3149   if (!disableUncertainties)
3150   {
3151     if (histId == 0)
3152       DrawLatex(0.3, 0.85, 1, "Scale uncertainty: 20%", 0.04);
3153     else
3154       DrawLatex(0.3, 0.85, 1, "Scale uncertainty: 10%", 0.04);
3155   }
3156   DrawLatex(0.3, 0.80, 1, Form("#Delta#varphi projection in |#Delta#eta| < %.2f", hist->GetYaxis()->GetBinUpEdge(hist->GetYaxis()->FindBin(projectLimit-0.01))), 0.04);
3157   DrawLatex(0.3, 0.75, 2, Form("#Delta#eta projection in |#Delta#varphi| < %.2f", hist->GetXaxis()->GetBinUpEdge(hist->GetXaxis()->FindBin(projectLimit-0.01))), 0.04);
3158   
3159   proj->SetTitle(label);
3160   proj2->SetTitle(label);
3161   systUncPhi->SetTitle(label);
3162   systUncEta->SetTitle(label);
3163   
3164   TH1* etaRatio = (TH1*) proj2->Clone("etaRatio");
3165   etaRatio->GetXaxis()->SetRangeUser(0,2);
3166   etaRatio->GetYaxis()->SetRangeUser(etaRatio->GetMinimum(), etaRatio->GetMaximum() * 1.6);
3167   TH1* phiRatio = (TH1*) proj->Clone("phiRatio");
3168   phiRatio->GetXaxis()->SetRangeUser(0,2);
3169   phiRatio->GetYaxis()->SetRangeUser(phiRatio->GetMinimum(), phiRatio->GetMaximum() * 1.6);
3170
3171   for (Int_t x=1; x<=proj->GetNbinsX(); x++)
3172     if (proj->GetBinCenter(x)>=0)
3173     {
3174       phiRatio->SetBinContent(x,proj->GetBinContent(proj->FindBin(proj->GetBinCenter(x)))+proj->GetBinContent(proj->FindBin(-1*proj->GetBinCenter(x))));
3175       phiRatio->SetBinError(x,sqrt((proj->GetBinError(proj->FindBin(proj->GetBinCenter(x)))*proj->GetBinError(proj->FindBin(proj->GetBinCenter(x)))+proj->GetBinError(proj->FindBin(-1*proj->GetBinCenter(x)))*proj->GetBinError(proj->FindBin(-1*proj->GetBinCenter(x))))));
3176     }
3177   for (Int_t x=1; x<=proj2->GetNbinsX(); x++)
3178     if (proj2->GetBinCenter(x)>=0)
3179     {
3180       etaRatio->SetBinContent(x,proj2->GetBinContent(proj2->FindBin(proj2->GetBinCenter(x)))+proj2->GetBinContent(proj2->FindBin(-1*proj2->GetBinCenter(x))));
3181       etaRatio->SetBinError(x,sqrt((proj2->GetBinError(proj2->FindBin(proj2->GetBinCenter(x)))*proj2->GetBinError(proj2->FindBin(proj2->GetBinCenter(x)))+proj2->GetBinError(proj2->FindBin(-1*proj2->GetBinCenter(x)))*proj2->GetBinError(proj2->FindBin(-1*proj2->GetBinCenter(x))))));
3182     }
3183
3184   if (Ratio)
3185   {
3186     *projEta = etaRatio;
3187     *projPhi = phiRatio;
3188   }
3189   else
3190   {
3191     *projPhi = proj;
3192     *projEta = proj2;
3193   }
3194   *projPhiSyst = systUncPhi;
3195   *projEtaSyst = systUncEta;
3196
3197   c->SaveAs(Form("ex/%s.png", c->GetName()));
3198   c->SaveAs(Form("ex/%s.eps", c->GetName()));
3199
3200   if (1)
3201   {
3202     c = new TCanvas(Form("ex_%d_%d_%d_a", i, j, histId), Form("ex_%d_%d_%d_a", i, j, histId), 800, 800);
3203     gPad->SetLeftMargin(0.15);
3204     hist->SetTitle("");
3205     hist->DrawCopy("SURF1");
3206     paveText->Draw();
3207     DrawALICELogo(kTRUE, 0.2, 0.75, 0.4, 0.95);
3208     c->SaveAs(Form("ex/%s.png", c->GetName()));
3209     c->SaveAs(Form("ex/%s.eps", c->GetName()));
3210
3211     c = new TCanvas(Form("ex_%d_%d_%d_b", i, j, histId), Form("ex_%d_%d_%d_b", i, j, histId), 800, 800);
3212     gPad->SetLeftMargin(0.15);
3213     clone->SetTitle("");
3214     clone->DrawCopy("SURF1");
3215     paveText->Draw();
3216     DrawALICELogo(kTRUE, 0.2, 0.75, 0.4, 0.95);
3217     c->SaveAs(Form("ex/%s.png", c->GetName()));
3218     c->SaveAs(Form("ex/%s.eps", c->GetName()));
3219   }
3220 }
3221
3222 void DrawExampleAll(const char* histFileName, const char* graphFileName = 0)
3223 {
3224   if (graphFileName)
3225     ReadGraphs(graphFileName);
3226   
3227   Int_t colors[] = { 1, 2, 4 };
3228   
3229   //Float_t exampleI[] = { 0, 2, 2 };
3230   //Float_t exampleJ[] = { 1, 1, 2 };
3231   //Int_t exampleI[] = { 0, 1, 2 };
3232   //Int_t exampleJ[] = { 1, 1, 1 };
3233   Float_t exampleI[] = { 3, 3, 3 };
3234   Float_t exampleJ[] = { 3, 4, 5 };
3235   
3236   TH1* projectionsPhi[9];
3237   TH1* projectionsEta[9];
3238   TH1* projectionsPhiSyst[9];
3239   TH1* projectionsEtaSyst[9];
3240
3241   Bool_t Ratio = kFALSE;
3242   //if Ratio is true it divides by the 60-70% centrality bin the other two 
3243   //the label on the Y axis is wrong
3244
3245   disableUncertainties = 1;
3246   
3247   Int_t count = 0;
3248   for (Int_t i=0; i<3; i++)
3249   {
3250     for (Int_t histId = 0; histId<3; histId++)
3251     {
3252       //it is using the 0-10%, the 20-40% and the 60-70% centrality bins
3253       DrawExample(histFileName, exampleI[i], exampleJ[i], (histId==2) ? 3 : histId, &projectionsPhi[count], &projectionsEta[count], &projectionsPhiSyst[count], &projectionsEtaSyst[count], Ratio);
3254       count++;
3255     }
3256     
3257     TCanvas* c = new TCanvas(Form("centralities_%d", i), Form("centralities_%d", i), 1200, 600);
3258     c->Divide(2, 1);
3259     
3260     TLegend* legend = new TLegend(0.62, 0.7, 0.88, 0.88);
3261     legend->SetFillColor(0);
3262
3263     // syst first
3264     for (Int_t histId = 0; histId<3; histId++)
3265     {
3266       if (Ratio)
3267       {
3268         projectionsPhi[count-3+histId]->Divide(projectionsPhi[count-3+1]);
3269         projectionsEta[count-3+histId]->Divide(projectionsEta[count-3+1]);
3270
3271       }
3272       if(!Ratio || (Ratio && histId!=1))
3273       {
3274       c->cd(1);
3275       TH1* clone;
3276       gPad->SetLeftMargin(0.12);
3277       projectionsPhiSyst[count-3+histId]->SetFillStyle(3003+histId);
3278       if (!disableUncertainties) clone = projectionsPhiSyst[count-3+histId]->DrawCopy((histId > 0) ? "E2 ][ SAME" : "E2 ][");
3279       else clone = projectionsPhi[count-3+histId]->DrawCopy((histId > 0) ? "E ][ SAME" : "E ][");
3280       clone->SetFillColor(colors[histId]);
3281       TString label(clone->GetTitle());
3282       TObjArray* objArray = label.Tokenize("-");
3283       clone->SetTitle(Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()));
3284     
3285       c->cd(2);
3286       gPad->SetLeftMargin(0.12);
3287       projectionsEtaSyst[count-3+histId]->SetFillStyle(3003+histId);
3288       if (!disableUncertainties) clone = projectionsEtaSyst[count-3+histId]->DrawCopy((histId > 0) ? "E2 ][ SAME" : "E2 ][");
3289       else clone = projectionsEta[count-3+histId]->DrawCopy((histId > 0) ? "E ][ SAME" : "E ][");
3290       clone->SetFillColor(colors[histId]);
3291       clone->SetTitle(Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()));
3292       }
3293     }
3294     
3295     for (Int_t histId = 0; histId<3; histId++)
3296     {
3297       if(!Ratio || (Ratio && histId!=1))
3298       {
3299       c->cd(1);
3300       TH1* clone = projectionsPhi[count-3+histId]->DrawCopy("SAME");
3301       clone->SetLineColor(colors[histId]);
3302       TString label(clone->GetTitle());
3303       TObjArray* objArray = label.Tokenize("-");
3304       legend->AddEntry(clone, (objArray->GetEntries() == 4) ? Form("%s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()) : objArray->At(2)->GetName(), "L");
3305       
3306       c->cd(2);
3307       clone = projectionsEta[count-3+histId]->DrawCopy("SAME");
3308       clone->SetLineColor(colors[histId]);
3309       }
3310     }
3311     
3312     c->cd(1);
3313     legend->Draw();
3314     DrawLatex(0.15, 0.86,  1, "Pb-Pb #sqrt{s_{NN}} = 2.76 TeV", 0.03);
3315     DrawLatex(0.15, 0.82,  1, "pp #sqrt{s} = 2.76 TeV", 0.03);
3316     DrawLatex(0.15, 0.78, 1, "|#eta| < 0.9", 0.03);
3317     DrawLatex(0.15, 0.74, 1, "Projected within |#Delta#eta| < 0.80", 0.03);
3318
3319     c->cd(2);
3320     DrawLatex(0.15, 0.86, 1, "Scale uncertainty: 20% (for 0-10%) / 10% (other bins)", 0.03);
3321     DrawLatex(0.15, 0.82, 1, "Projected within |#Delta#varphi| < 0.87", 0.03);
3322     if (!disableUncertainties) DrawALICELogo(kTRUE, 0.75, 0.62, 0.95, 0.78);
3323     
3324     c->SaveAs(Form("ex/%s.png", c->GetName()));
3325     c->SaveAs(Form("ex/%s.eps", c->GetName()));
3326  
3327 //     return;
3328
3329     if (graphFileName)
3330     {
3331       // fit curves
3332       Int_t graphID = exampleI[i] * (6 - 1) + exampleJ[i] - 1;
3333       for (Int_t histId = 0; histId<3; histId++)
3334       {
3335         c->cd(1);
3336         // integral over y
3337 //      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());
3338         // see /home/jgrosseo/limited2dgaussianintegration.nb
3339         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());
3340
3341         func->SetParameter(0, 0);
3342 /*      if (histId == 0) func->SetParameter(0, -0.00129472);
3343         if (histId == 1) func->SetParameter(0, -0.00025859);
3344         if (histId == 2) func->SetParameter(0, -0.000833301);
3345 */
3346         for (Int_t k=0; k<6; k++)
3347           func->SetParameter(k+1, graphs[k][graphID]->GetY()[(histId==2) ? 3 : histId]);
3348         // scale by bin width (to compare with projection)
3349         // NOTE this scaling looks inconsistent, but is OK. See normalization in DrawExample
3350         func->SetParameter(1, func->GetParameter(1) / projectionsEta[count-3+histId]->GetXaxis()->GetBinWidth(1));
3351         func->SetParameter(7, 0.8); // project Limit
3352         func->SetLineColor(colors[histId]);
3353         func->SetLineWidth(1);
3354         func->DrawCopy("SAME");
3355         
3356         // draw contributions
3357         Float_t scale = func->GetParameter(1);
3358         Float_t weighting = func->GetParameter(4);
3359         func->SetParameter(1, scale * weighting);
3360         func->SetParameter(4, 1);
3361         func->SetLineStyle(2);
3362         if (weighting > 0.01)
3363           func->DrawCopy("SAME");
3364
3365         func->SetParameter(1, scale * (1.0-weighting));
3366         func->SetParameter(4, 0);
3367         func->SetLineStyle(3);
3368         if (1.0-weighting > 0.01)
3369           func->DrawCopy("SAME");
3370         
3371         c->cd(2);
3372         // integral over x
3373         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);
3374         func->SetParameter(0, 0);
3375         for (Int_t k=0; k<6; k++)
3376           func->SetParameter(k+1, graphs[k][graphID]->GetY()[(histId==2) ? 3 : histId]);
3377         // scale by bin width (to compare with projection)
3378         // NOTE this scaling looks inconsistent, but is OK. See normalization in DrawExample
3379         func->SetParameter(1, func->GetParameter(1) / projectionsEta[count-3+histId]->GetXaxis()->GetBinWidth(1));
3380         func->SetParameter(7, 0.87); // project Limit
3381         func->SetLineColor(colors[histId]);
3382         func->SetLineWidth(1);
3383         func->DrawCopy("SAME");
3384
3385         // draw contributions
3386         scale = func->GetParameter(1);
3387         weighting = func->GetParameter(4);
3388         func->SetParameter(1, scale * weighting);
3389         func->SetParameter(4, 1);
3390         func->SetLineStyle(2);
3391         if (weighting > 0.01)
3392           func->DrawCopy("SAME");
3393
3394         func->SetParameter(1, scale * (1.0-weighting));
3395         func->SetParameter(4, 0);
3396         func->SetLineStyle(3);
3397         if (1.0-weighting > 0.01)
3398           func->DrawCopy("SAME");
3399       }
3400       
3401       c->SaveAs(Form("ex/%s_fits.png", c->GetName()));
3402       c->SaveAs(Form("ex/%s_fits.eps", c->GetName()));
3403       
3404 //       return;
3405     }
3406   }
3407
3408   return;
3409   // TODO implement syst uncertainty plotting
3410   
3411   for (Int_t histId = 0; histId<3; histId++)
3412   {
3413     TCanvas* c = new TCanvas(Form("pt_%d", histId), Form("pt_%d", histId), 1200, 600);
3414     c->Divide(2, 1);
3415     
3416     TLegend* legend = new TLegend(0.15, 0.7, 0.88, 0.88);
3417     legend->SetFillColor(0);
3418     
3419     for (Int_t i=2; i>=0; i--)
3420     {
3421       c->cd(1);
3422       TH1* clone = projectionsPhi[i*3+histId]->DrawCopy((i < 2) ? "SAME" : "");
3423       clone->SetLineColor(colors[i]);
3424       clone->GetYaxis()->SetRangeUser(clone->GetMinimum(), clone->GetMaximum() * 1.2);
3425       
3426       TString label(clone->GetTitle());
3427       TObjArray* objArray = label.Tokenize("-");
3428       clone->SetTitle((objArray->GetEntries() == 4) ? Form("%s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()) : objArray->At(2)->GetName());
3429       
3430       legend->GetListOfPrimitives()->AddFirst(new TLegendEntry(clone, Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()), "L"));
3431       
3432       c->cd(2);
3433       clone = projectionsEta[i*3+histId]->DrawCopy((i < 2) ? "SAME" : "");
3434       clone->SetTitle((objArray->GetEntries() == 4) ? Form("%s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()) : objArray->At(2)->GetName());
3435       clone->SetLineColor(colors[i]);
3436 //       clone->GetYaxis()->SetRangeUser(clone->GetMinimum(), clone->GetMaximum() * 1.2);
3437
3438       if (!disableUncertainties) DrawALICELogo(kTRUE, 0.65, 0.65, 0.85, 0.85);
3439     }
3440     
3441     c->cd(1);
3442     legend->Draw();
3443
3444     c->SaveAs(Form("ex/%s.png", c->GetName()));
3445     c->SaveAs(Form("ex/%s.eps", c->GetName()));
3446   }
3447 }
3448
3449 void DrawDoubleHump(const char* histFileName)
3450 {
3451   Int_t exampleI[] = { 0, 1, 2, 0, 1, 2};
3452   Int_t exampleJ[] = { 1, 1, 1, 2, 2, 2};
3453   
3454   TH1* projectionsPhi[9];
3455   TH1* projectionsEta[9];
3456   TH1* projectionsPhiSyst[9];
3457   TH1* projectionsEtaSyst[9];
3458   
3459   TCanvas* c = new TCanvas("DrawDoubleHump", "DrawDoubleHump", 1200, 800);
3460   c->Divide(3, 2);
3461   
3462   TLegend* legend = new TLegend(0.15, 0.7, 0.88, 0.88);
3463   legend->SetFillColor(0);
3464   
3465   for (Int_t i=0; i<6; i++)
3466   {
3467     if (i == 3)
3468       continue;
3469     
3470     DrawExample(histFileName, exampleI[i], exampleJ[i], 0, &projectionsPhi[i], &projectionsEta[i], &projectionsPhiSyst[i], &projectionsEtaSyst[i]);
3471
3472     c->cd(i+1);
3473     gPad->SetLeftMargin(0.15);
3474     projectionsPhiSyst[i]->GetXaxis()->SetTitle(Form("%s , %s", "#Delta#varphi (rad.)", "#Delta#eta"));
3475     projectionsPhiSyst[i]->GetYaxis()->SetTitle(kProjYieldTitlePhiOrEta);
3476     projectionsPhiSyst[i]->GetYaxis()->SetTitleOffset(1.8);
3477 //     projectionsPhiSyst[i]->SetTitleSize(0.035);
3478     
3479     TString label(projectionsPhiSyst[i]->GetTitle());
3480     TObjArray* objArray = label.Tokenize("-");
3481     projectionsPhiSyst[i]->SetTitle(Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()));
3482     
3483     projectionsPhiSyst[i]->Draw("E2 ][");
3484     projectionsEtaSyst[i]->Draw("E2 ][ SAME");
3485     TH1* clone = projectionsPhi[i]->DrawCopy((disableUncertainties) ? "" : "SAME");
3486     clone->SetLineColor(1);
3487     clone->GetYaxis()->SetRangeUser(clone->GetMinimum(), clone->GetMaximum() * 1.2);
3488 //     clone->GetXaxis()->SetTitle(Form("%s , %s", clone->GetXaxis()->GetTitle(), "#Delta#eta"));
3489     clone->GetXaxis()->SetTitle(Form("%s , %s", "#Delta#varphi (rad.)", "#Delta#eta"));
3490     clone->GetYaxis()->SetTitle(kProjYieldTitlePhiOrEta);
3491     
3492     clone = projectionsEta[i]->DrawCopy("SAME");
3493     clone->SetLineColor(2);
3494
3495     if (!disableUncertainties)
3496       DrawLatex(0.3, 0.85, 1, "Scale uncertainty: 20%", 0.04);
3497     
3498     DrawLatex(0.3, 0.81, 1, Form("#Delta#varphi projection in |#Delta#eta| < %.2f", 0.8), 0.04);
3499     DrawLatex(0.3, 0.77, 2, Form("#Delta#eta projection in |#Delta#varphi| < %.2f", 0.87), 0.04);
3500 //     DrawLatex(0.3, 0.75, 1, Form("%s-%s centrality", objArray->At(2)->GetName(), objArray->At(3)->GetName()), 0.04);
3501     DrawLatex(0.3, 0.73, 1, "0-10% centrality", 0.04);
3502     
3503 //     return;
3504   }
3505   
3506   if (!disableUncertainties)
3507     DrawALICELogo(kTRUE, 0.65, 0.5, 0.85, 0.7);
3508   
3509   c->SaveAs(Form("ex/%s.png", c->GetName()));
3510   c->SaveAs(Form("ex/%s.eps", c->GetName()));
3511 }
3512
3513 void DrawFullCentralityDependence(const char* histFileName)
3514 {
3515   Int_t colors[] = { 1, 2, 5, 4, 3, 6 };
3516   
3517   TH1* projectionsPhi[9];
3518   TH1* projectionsEta[9];
3519   TH1* projectionsPhiSyst[9];
3520   TH1* projectionsEtaSyst[9];
3521   
3522   Int_t count = 0;
3523   for (Int_t histId = 0; histId<6; histId++)
3524   {
3525     projectionsPhi[count] = 0;
3526     DrawExample(histFileName, 0, 1, histId, &projectionsPhi[count], &projectionsEta[count], &projectionsPhiSyst[count], &projectionsEtaSyst[count]);
3527     count++;
3528   }
3529   
3530   TCanvas* c = new TCanvas(Form("centralities_%d", 0), Form("centralities_%d", 0), 800, 400);
3531   c->Divide(2, 1);
3532   
3533   TLegend* legend = new TLegend(0.62, 0.7, 0.88, 0.88);
3534   legend->SetFillColor(0);
3535   
3536   for (Int_t histId = 0; histId<6; histId++)
3537   {
3538     if (histId == 2)
3539       continue;
3540     
3541     if (!projectionsPhi[count-6+histId])
3542       continue;
3543     
3544     c->cd(1);
3545     TH1* clone = projectionsPhi[count-6+histId]->DrawCopy((histId > 0) ? "SAME" : "");
3546     clone->SetLineColor(colors[histId]);
3547     
3548     TString label(clone->GetTitle());
3549     TObjArray* objArray = label.Tokenize("-");
3550     clone->SetTitle(Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()));
3551     
3552     legend->AddEntry(clone, (objArray->GetEntries() == 4) ? Form("%s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()) : objArray->At(2)->GetName(), "L");
3553
3554     c->cd(2);
3555     clone = projectionsEta[count-6+histId]->DrawCopy((histId > 0) ? "SAME" : "");
3556     clone->SetTitle(Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName()));
3557     clone->SetLineColor(colors[histId]);
3558   }
3559     
3560   c->cd(1);
3561   legend->Draw();
3562   c->SaveAs(Form("ex/%s.png", c->GetName()));
3563   c->SaveAs(Form("ex/%s.eps", c->GetName()));
3564 }
3565
3566 void CompareExamples(const char* histFileName1, const char* histFileName2, Int_t i, Int_t j, Int_t histId)
3567 {
3568   Int_t exColors[] = { 1, 2, 4, 6 };
3569   
3570   TCanvas* c = new TCanvas(Form("c_%d", i), Form("c_%d", i), 1000, 600);
3571   c->Divide(2, 1);
3572   
3573   TCanvas* c2 = new TCanvas(Form("c_%d_ratio", i), Form("c_%d_ratio", i), 1000, 600);
3574   c2->Divide(2, 1);
3575   
3576   TH1* projFirst[2];
3577   
3578   for (Int_t k=0; k<2; k++)
3579   {
3580     if (k == 0)
3581       TFile::Open(histFileName1);
3582     else
3583       TFile::Open(histFileName2);
3584   
3585     TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
3586     if (!hist)
3587       continue;
3588 //     if (k == 1)
3589 //       hist->Scale(2);
3590 //     hist->Rebin2D(2, 2);
3591
3592 //     new TCanvas; hist->DrawCopy("COLZ"); return;
3593
3594     
3595     SubtractEtaGap(hist, kEtaLimit, kOuterLimit, kFALSE);
3596     
3597     // new TCanvas; hist->DrawCopy("COLZ");
3598
3599     TH1* proj1 = 0;
3600     TH1* proj2 = 0;
3601     GetProjections(hist, &proj1, &proj2, k);
3602     
3603     c->cd(1);
3604     proj1->SetLineColor(exColors[k]);
3605     proj1->GetXaxis()->SetRangeUser(-1.49, 1.49);
3606     proj1->GetYaxis()->SetRangeUser(proj1->GetMinimum() * 1.2, proj1->GetMaximum() * 1.5);
3607     proj1->DrawCopy((k == 0) ? "" : "SAME");
3608   
3609     if (k == 1)
3610     {
3611       c2->cd(1);
3612       proj1->Divide(projFirst[0]);
3613       proj1->Draw();
3614       proj1->GetYaxis()->SetRangeUser(0.5, 1.5);
3615     }
3616     else
3617       projFirst[0] = proj1;
3618
3619     c->cd(2);
3620     proj2->SetLineColor(exColors[k]);
3621     proj2->GetXaxis()->SetRangeUser(-1.49, 1.49);
3622     proj2->GetYaxis()->SetRangeUser(proj2->GetMinimum() * 1.2, proj2->GetMaximum() * 2);
3623     
3624     if (0 && k == 1)
3625     {
3626       Printf("Scaling!");
3627       proj2->Scale(projFirst[1]->Integral(projFirst[1]->FindBin(-0.5), projFirst[1]->FindBin(0.5)) / proj2->Integral(projFirst[1]->FindBin(-0.5), projFirst[1]->FindBin(0.5)));
3628     }
3629     
3630     proj2->DrawCopy((k == 0) ? "" : "SAME");
3631
3632     if (k == 1)
3633     {
3634       c2->cd(2);
3635       proj2->Divide(projFirst[1]);
3636       proj2->Draw();
3637       proj2->GetYaxis()->SetRangeUser(0.5, 1.5);
3638     }
3639     else
3640       projFirst[1] = proj2;
3641   }
3642 }
3643
3644 void TestMomentCode()
3645 {
3646   Float_t momentLimit = 1;
3647   TH1* proj = new TH1F("hist", "", 100, -momentLimit, momentLimit);
3648   TF1* func = new TF1("func", "gaus(0)", -momentLimit, momentLimit);
3649   func->SetParameters(1, 0, 0.2);
3650   proj->FillRandom("func", 100000);
3651   proj->Sumw2();
3652   proj->Draw();
3653   
3654   for (Int_t n=2; n <= 4; n++)
3655   {
3656     Float_t moment = 0;
3657     Float_t sum = 0;
3658     for (Int_t bin = proj->GetXaxis()->FindBin(-momentLimit); bin <= proj->GetXaxis()->FindBin(momentLimit); bin++)
3659     {
3660       moment += proj->GetBinContent(bin) * TMath::Power(proj->GetMean() - proj->GetXaxis()->GetBinCenter(bin), n);
3661       sum += proj->GetBinContent(bin);
3662     }
3663     moment /= sum;
3664     
3665     Float_t error = 0;
3666     for (Int_t bin = proj->GetXaxis()->FindBin(-momentLimit); bin <= proj->GetXaxis()->FindBin(momentLimit); bin++)
3667     {
3668       error += proj->GetBinError(bin) * proj->GetBinError(bin) * 
3669         TMath::Power(TMath::Power(proj->GetMean() - proj->GetXaxis()->GetBinCenter(bin), n) / sum 
3670           - moment / sum, 2);
3671     }
3672     
3673     Printf("%d %f +- %f <-> %f +- %f", n, moment, TMath::Sqrt(error), proj->GetRMS() * proj->GetRMS(), 2 * proj->GetRMSError() / proj->GetRMS() * proj->GetRMS() * proj->GetRMS());
3674   }
3675 }
3676
3677 void CompareGraph(const char* fileName1, const char* fileName2, Int_t graph, Int_t centrality)
3678 {
3679   ReadGraphs(fileName1);
3680   TGraphErrors* graph1 = (TGraphErrors*) graphs[graph][centrality]->Clone("graph1");
3681   
3682   ReadGraphs(fileName2);
3683   TGraphErrors* graph2 = (TGraphErrors*) graphs[graph][centrality]->Clone("graph2");
3684   
3685   graph1->Sort();
3686   graph2->Sort();
3687   
3688   TCanvas* c = new TCanvas(Form("%s_%s", fileName1, fileName2), Form("%s_%s", fileName1, fileName2), 800, 800);
3689   c->Divide(1, 2);
3690   
3691   c->cd(1);
3692   gPad->SetGridx();
3693   gPad->SetGridy();
3694   graph1->SetMarkerStyle(24);
3695   graph1->GetXaxis()->SetRangeUser(7, 38);
3696   graph1->DrawClone("AP");
3697   graph2->SetLineColor(2);
3698   graph2->SetMarkerColor(2);
3699   graph2->SetMarkerStyle(25);
3700   graph2->DrawClone("PSAME");
3701   
3702   c->cd(2);
3703   gPad->SetGridx();
3704   gPad->SetGridy();
3705   DivideGraphs(graph1, graph2);
3706   graph1->DrawClone("AP");
3707 }
3708
3709 void Compare2D(const char* fileName1, const char* fileName2, const char* histName)
3710 {
3711   TFile::Open(fileName1);
3712   TH1* hist1 = (TH1*) gFile->Get(histName);
3713
3714   TFile::Open(fileName2);
3715   TH1* hist2 = (TH1*) gFile->Get(histName);
3716   
3717   hist1->Divide(hist2);
3718   hist1->GetZaxis()->SetRangeUser(0.5, 1.5);
3719   hist1->Draw("COLZ");
3720 }
3721
3722 void TestTwoGaussian()
3723 {
3724 //   TF1* func = new TF1("func", "[0]/TMath::Sqrt(TMath::TwoPi())/[2]*exp(-0.5*((x-[1])/[2])**2) + [3]/TMath::Sqrt(TMath::TwoPi())/[5]*exp(-0.5*((x-[4])/[5])**2)", -2, 2);
3725   TF1* func = new TF1("func", "[0]*exp(-0.5*((x-[1])/[2])**2) + [3]*exp(-0.5*((x-[4])/[5])**2)", -2, 2);
3726   func->SetParameters(0.25, 0, 0.5, 0.75, 0, 0.2);
3727   
3728   TH1* hist = new TH1F("hist", "", 100, -2, 2);
3729   hist->FillRandom("func", 10000000);
3730   
3731   hist->Draw();
3732   
3733   TH1* proj = hist;
3734   
3735   Float_t momentLimit = 1.99;
3736   for (Int_t n=2; n <= 4; n++)
3737   {
3738     Float_t moment = 0;
3739     Float_t sum = 0;
3740     for (Int_t bin = proj->GetXaxis()->FindBin(-momentLimit); bin <= proj->GetXaxis()->FindBin(momentLimit); bin++)
3741     {
3742       moment += proj->GetBinContent(bin) * TMath::Power(proj->GetMean() - proj->GetXaxis()->GetBinCenter(bin), n);
3743       sum += proj->GetBinContent(bin);
3744     }
3745     moment /= sum;
3746     
3747     Float_t error = 0;
3748     for (Int_t bin = proj->GetXaxis()->FindBin(-momentLimit); bin <= proj->GetXaxis()->FindBin(momentLimit); bin++)
3749     {
3750       error += proj->GetBinError(bin) * proj->GetBinError(bin) * 
3751         TMath::Power(TMath::Power(proj->GetMean() - proj->GetXaxis()->GetBinCenter(bin), n) / sum 
3752           - moment / sum, 2);
3753     }
3754     
3755     Printf("%d %f %f", n, moment, TMath::Sqrt(error));
3756   }
3757 }
3758
3759 void DrawEtaGapExample(const char* histFileName, Int_t i = 2, Int_t j = 2, Int_t histId = 0)
3760 {
3761   Float_t etaLimit = kEtaLimit;
3762   Float_t outerLimit = kOuterLimit;
3763   Float_t projectLimit = 0.8;
3764
3765   TFile::Open(histFileName);
3766
3767   TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, histId));
3768   if (!hist)
3769     return;
3770   hist->GetZaxis()->SetTitle(kCorrFuncTitle);
3771   hist->GetZaxis()->SetTitleOffset(1.8);
3772   // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
3773   hist->Scale(1.0 / hist->GetYaxis()->GetBinWidth(1));
3774   
3775   TCanvas* c = new TCanvas("c",&