6b3ad8654f02f68faa0eb1d7cdf6123099e253cb
[u/mrichter/AliRoot.git] / PWGCF / Correlations / macros / dphicorrelations / pA.C
1 // 0       1       2       3       4  5  6              7     8             9                            10         11
2 // nsyield,asyield,nswidth,aswidth,v2,v3,nsasyieldratio,v3/v2,remainingpeak,remainingjetyield/ridgeyield,chi2(v2v3),baseline
3 const char* graphTitles[] = { "NS Ridge Yield", "AS Ridge Yield", "NS Width", "AS Width", "v2", "v3", "NS Yield / AS Yield", "v3 / v2", "remaining peak in %", "remaining jet / NS yield in %", "chi2/ndf", "baseline" };
4 const Int_t NGraphs = 28; // pt index
5 const Int_t NHists = 12; 
6 TGraphErrors*** graphs = 0;
7
8 const char* kCorrFuncTitle = "#frac{1}{#it{N}_{trig}} #frac{d^{2}#it{N}_{assoc}}{d#Delta#etad#Delta#varphi} (rad^{-1})";
9 // const char* kCorrFuncTitle = "1/N_{trig} dN_{assoc}/d#Delta#etad#Delta#varphi (1/rad)";
10 const char* kProjYieldTitlePhi = "1/#it{N}_{trig} d#it{N}_{assoc}/d#Delta#varphi per #Delta#eta (rad^{-1})";
11 const char* kProjYieldTitleEta = "1/#it{N}_{trig} d#it{N}_{assoc}/d#Delta#eta per #Delta#varphi (rad^{-1})";
12 // const char* kProjYieldTitlePhiOrEta = "1/N_{trig} dN_{assoc}/d#Delta#varphi (1/rad) , dN_{assoc}/d#Delta#eta";
13 const char* kProjYieldTitlePhiOrEta = "#frac{1}{#it{N}_{trig}} #frac{d#it{N}_{assoc}}{d#Delta#varphi} (rad^{-1}) , d#it{N}_{assoc}/d#Delta#eta";
14
15 Float_t etaMax = 1.6;
16
17 void CreateGraphStructure()
18 {
19   graphs = new TGraphErrors**[NGraphs];
20   for (Int_t i=0; i<NGraphs; i++)
21   {
22     graphs[i] = new TGraphErrors*[NHists];
23     for (Int_t j=0; j<NHists; j++)
24     {
25       graphs[i][j] = new TGraphErrors;
26       graphs[i][j]->GetYaxis()->SetTitle(graphTitles[j]);
27     }
28   }
29 }
30
31 void WriteGraphs(const char* outputFileName = "graphs.root")
32 {
33   TFile::Open(outputFileName, "RECREATE");
34   for (Int_t i=0; i<NGraphs; i++)
35   {
36     if (!graphs[i])
37       continue;
38     for (Int_t j=0; j<NHists; j++)
39     {
40       if (!graphs[i][j])
41         continue;
42       graphs[i][j]->GetYaxis()->SetTitle(graphTitles[j]);
43       graphs[i][j]->Write(Form("graph_%d_%d", i, j));
44     }
45   }
46
47   gFile->Close();
48 }
49
50 void ReadGraphs(const char* fileName = "graphs.root")
51 {
52   CreateGraphStructure();
53   TFile* file = TFile::Open(fileName);
54   for (Int_t i=0; i<NGraphs; i++)
55     for (Int_t j=0; j<NHists; j++)
56       graphs[i][j] = (TGraphErrors*) file->Get(Form("graph_%d_%d", i, j));
57 }
58
59 void AddPoint(TGraphErrors* graph, Float_t x, Float_t y, Float_t xe, Float_t ye)
60 {
61         graph->SetPoint(graph->GetN(), x, y);
62         graph->SetPointError(graph->GetN() - 1, xe, ye);
63 }
64
65 void RemovePointsBelowX(TGraphErrors* graph, Float_t minX)
66 {
67   Int_t i=0;
68   while (i<graph->GetN())
69   {
70     if (graph->GetX()[i] < minX)
71       graph->RemovePoint(i);
72     else
73       i++;
74   }
75 }
76
77 void RemovePointsAboveX(TGraphErrors* graph, Float_t maxX)
78 {
79   Int_t i=0;
80   while (i<graph->GetN())
81   {
82     if (graph->GetX()[i] > maxX)
83       graph->RemovePoint(i);
84     else
85       i++;
86   }
87 }
88
89 void DrawLatex(Float_t x, Float_t y, Int_t color, const char* text, Float_t textSize = 0.06)
90 {
91   TLatex* latex = new TLatex(x, y, text);
92   latex->SetNDC();
93   latex->SetTextSize(textSize);
94   latex->SetTextColor(color);
95   latex->Draw();
96 }
97
98 void PadFor2DCorr()
99 {
100   gPad->SetPad(0, 0, 1, 1);
101   gPad->SetLeftMargin(0.2);
102   gPad->SetTopMargin(0.05);
103   gPad->SetRightMargin(0.05);
104   gPad->SetBottomMargin(0.05);
105   gPad->SetTheta(61.62587);
106   gPad->SetPhi(45);
107 }
108
109 void DivideGraphs(TGraphErrors* graph1, TGraphErrors* graph2)
110 {
111 //   graph1->Print();
112 //   graph2->Print();
113
114   graph1->Sort();
115   graph2->Sort();
116
117   for (Int_t bin1 = 0; bin1 < graph1->GetN(); bin1++)
118   {
119     Float_t x = graph1->GetX()[bin1];
120
121     Int_t bin2 = 0;
122     for (bin2 = 0; bin2<graph2->GetN(); bin2++)
123       if (graph2->GetX()[bin2] >= x)
124         break;
125
126     if (bin2 == graph2->GetN())
127             bin2--;
128
129     if (bin2 > 0)
130       if (TMath::Abs(graph2->GetX()[bin2-1] - x) < TMath::Abs(graph2->GetX()[bin2] - x))
131         bin2--;
132
133     if (graph1->GetY()[bin1] == 0 || graph2->GetY()[bin2] == 0 || bin2 == graph2->GetN())
134     {
135       Printf("%d %d removed", bin1, bin2);
136       graph1->RemovePoint(bin1--);
137       continue;
138     }
139
140     Float_t graph2Extrapolated = graph2->GetY()[bin2];
141     if (TMath::Abs(x - graph2->GetX()[bin2]) > 0.0001)
142     {
143       Printf("%f %f %d %d not found", x, graph2->GetX()[bin2], bin1, bin2);
144       graph1->RemovePoint(bin1--);
145       continue;
146     }
147
148     Float_t value = graph1->GetY()[bin1] / graph2Extrapolated;
149     Float_t error = value * TMath::Sqrt(TMath::Power(graph1->GetEY()[bin1] / graph1->GetY()[bin1], 2) + TMath::Power(graph2->GetEY()[bin2] / graph2->GetY()[bin2], 2));
150
151     graph1->GetY()[bin1] = value;
152     graph1->GetEY()[bin1] = error;
153
154 //     Printf("%d %d %f %f %f %f", bin1, bin2, x, graph2Extrapolated, value, error);
155   }
156 }
157
158 void GraphShiftX(TGraphErrors* graph, Float_t offset)
159 {
160   for (Int_t i=0; i<graph->GetN(); i++)
161     graph->GetX()[i] += offset;
162 }
163
164 TGraphErrors* ReadHepdata(const char* fileName, Bool_t errorsAreAdded = kFALSE, Int_t skipYErrors = 0, Int_t skipXerrors = 1)
165 {
166   // expected format: x [x2] y [ye] [ye2] [xe]
167   //
168   // skipYErrors:   0 --> ye present
169   //                1 --> no errors ye
170   //                2 --> y and ye are lower and upper error, i.e. y' = (y + ye) / 2 and ye = (ye - y) / 2
171   //                3 --> ye and ye2 are stat and syst error, will be added in quadrature
172   // 
173   // skipXerrors:   0 --> xe present
174   //                1 --> no errors xe
175   //                2 --> x2 present, xe not present and is calculated from x2 - x
176   
177   ifstream fin(fileName);
178
179   graph = new TGraphErrors(0);
180
181   Double_t sum = 0;
182
183   while (fin.good())
184   {
185     char buffer[2000];
186     if (fin.peek() == '#')
187     {
188       fin.getline(buffer, 2000);
189       continue;
190     }
191   
192     Double_t x = -1;
193     Double_t x2 = -1;
194     Double_t y = -1;
195     Double_t ye = 0;
196     Double_t xe = 0;
197
198     fin >> x;
199     
200     if (skipXerrors == 2)
201     {
202       fin >> x2;
203       xe = (x2 - x + 1) / 2;
204       x = x + (x2 - x) / 2;
205     }
206     
207     fin >> y;
208
209     if (y == -1)
210       continue;
211
212     if (skipYErrors == 0)
213     {
214       ye = -1;
215       fin >> ye;
216       if (ye == -1)
217         continue;
218     }
219     else if (skipYErrors == 2)
220     {
221       ye = -1;
222       fin >> ye;
223       if (ye == -1)
224         continue;
225       
226       Double_t newy = (y + ye) / 2;
227       ye = (ye - y) / 2;
228       y = newy;
229     }
230     else if (skipYErrors == 3)
231     {
232       ye = -1;
233       fin >> ye;
234       if (ye == -1)
235         continue;
236       
237       Double_t ye2 = -1;
238       fin >> ye2;
239       if (ye2 == -1)
240         continue;
241
242       ye = TMath::Sqrt(ye*ye + ye2*ye2);
243     }
244
245     if (skipXerrors == 0)
246     {
247       xe = -1;
248       fin >> xe;
249       if (xe == -1)
250         continue;
251     }
252
253     //Printf("%f %f %f %f", x, y, xe, ye);
254
255     if (errorsAreAdded)
256       ye -= y;
257
258     graph->SetPoint(graph->GetN(), x, y);
259     graph->SetPointError(graph->GetN()-1, xe, ye);
260
261     sum += y;
262     
263     // read rest until end of line...
264     fin.getline(buffer, 2000);
265   }
266   fin.close();
267
268   Printf("%s: %f", fileName, sum);
269
270   return graph;
271 }
272
273 TH2* SubtractEtaGapNS(TH2* hist, Float_t etaLimit, Float_t outerLimit, Bool_t drawEtaGapDist = kFALSE)
274 {
275   TString histName(hist->GetName());
276   Int_t etaBins = 0;
277
278   TH1D* etaGap = hist->ProjectionX(histName + "_1", TMath::Max(1, hist->GetYaxis()->FindBin(-outerLimit + 0.01)), hist->GetYaxis()->FindBin(-etaLimit - 0.01));
279 //   Printf("%f", etaGap->GetEntries());
280   if (etaGap->GetEntries() > 0)
281     etaBins += hist->GetYaxis()->FindBin(-etaLimit - 0.01) - TMath::Max(1, hist->GetYaxis()->FindBin(-outerLimit + 0.01)) + 1;
282
283   TH1D* tracksTmp = hist->ProjectionX(histName + "_2", hist->GetYaxis()->FindBin(etaLimit + 0.01), TMath::Min(hist->GetYaxis()->GetNbins(), hist->GetYaxis()->FindBin(outerLimit - 0.01)));
284 //   Printf("%f", tracksTmp->GetEntries());
285   if (tracksTmp->GetEntries() > 0)
286     etaBins += TMath::Min(hist->GetYaxis()->GetNbins(), hist->GetYaxis()->FindBin(outerLimit - 0.01)) - hist->GetYaxis()->FindBin(etaLimit + 0.01) + 1;
287   
288   etaGap->Add(tracksTmp);
289
290   // get per bin result
291   if (etaBins > 0)
292     etaGap->Scale(1.0 / etaBins);
293  
294   for (Int_t i=1; i<=etaGap->GetNbinsX()/2; i++)
295   {
296 //     Printf("%d -> %d", i, etaGap->GetNbinsX()+1-i);
297     etaGap->SetBinContent(etaGap->GetNbinsX()+1-i, etaGap->GetBinContent(i));
298     etaGap->SetBinError(etaGap->GetNbinsX()+1-i, etaGap->GetBinError(i));
299   }
300   
301   if (drawEtaGapDist)
302   {
303     TH1D* centralRegion = hist->ProjectionX(histName + "_3", hist->GetYaxis()->FindBin(-etaLimit + 0.01), hist->GetYaxis()->FindBin(etaLimit - 0.01));
304     
305 //    centralRegion->Scale(1.0 / (hist->GetYaxis()->FindBin(etaLimit - 0.01) - hist->GetYaxis()->FindBin(-etaLimit + 0.01) + 1));
306     centralRegion->Scale(hist->GetXaxis()->GetBinWidth(1));
307
308     TCanvas* c = new TCanvas("SubtractEtaGap", "SubtractEtaGap", 800, 800);
309     gPad->SetLeftMargin(0.13);
310     centralRegion->SetStats(0);
311     TString label(centralRegion->GetTitle());
312     label.ReplaceAll(".00", " GeV/#it{c}");
313     label.ReplaceAll(".0", " GeV/#it{c}");
314     centralRegion->SetTitle(label);
315     centralRegion->SetLineColor(3);
316     centralRegion->Draw();
317 //     centralRegion->GetYaxis()->SetTitle(kProjYieldTitlePhi);
318     centralRegion->GetYaxis()->SetTitleOffset(1.6);
319     TH1* copy = etaGap->DrawCopy("SAME");
320     copy->Scale(hist->GetXaxis()->GetBinWidth(1));
321     copy->Scale((hist->GetYaxis()->FindBin(etaLimit - 0.01) - hist->GetYaxis()->FindBin(-etaLimit + 0.01) + 1));
322     copy->SetLineColor(2);
323     TLegend* legend = new TLegend(0.41, 0.73, 0.69, 0.85);
324     legend->SetFillColor(0);
325     legend->AddEntry(centralRegion, Form("|#Delta#eta| < %.1f", etaLimit), "L");
326     legend->AddEntry(copy, Form("%.1f < |#Delta#eta| < %.1f (scaled)", etaLimit, outerLimit), "L");
327     legend->Draw();
328     
329 //     DrawLatex(0.705, 0.62, 1, "Pb-Pb 2.76 TeV", 0.025);
330 //     DrawLatex(0.705, 0.58, 1, "Stat. unc. only", 0.025);
331   }
332   
333 //   new TCanvas; etaGap->DrawCopy();
334   
335   TH2* histTmp2D = (TH2*) hist->Clone("histTmp2D");
336   histTmp2D->Reset();
337   
338   for (Int_t xbin=1; xbin<=histTmp2D->GetNbinsX(); xbin++)
339     for (Int_t y=1; y<=histTmp2D->GetNbinsY(); y++)
340       histTmp2D->SetBinContent(xbin, y, etaGap->GetBinContent(xbin));
341     
342   hist->Add(histTmp2D, -1);  
343   return histTmp2D;
344 }
345
346 TH1* GetProjections(Int_t i, Int_t j, Int_t centr, char** label, Float_t etaBegin = 1.0, TH1** etaProj = 0)
347 {
348   TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr));
349   if (!hist1)
350     return 0;
351   hist1 = (TH2*) hist1->Clone(Form("%s_%.1f", hist1->GetName(), etaBegin));
352
353   // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
354   hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
355   
356   hist1->Rebin2D(2, 1); hist1->Scale(0.5);
357   
358 //   new TCanvas; hist1->Draw("surf1");
359   
360   if (etaBegin > 0)
361     SubtractEtaGapNS(hist1, etaBegin, etaMax, kTRUE);
362   
363   tokens = TString(hist1->GetTitle()).Tokenize("-");
364   centralityStr = new TString;
365   if (tokens->GetEntries() > 2)
366     *centralityStr = tokens->At(2)->GetName();
367   if (tokens->GetEntries() > 3)
368     *centralityStr = *centralityStr + "-" + tokens->At(3)->GetName();
369   *label = centralityStr->Data();
370 //   Printf("%s", label);
371   
372   proj1x = ((TH2*) hist1)->ProjectionX(Form("proj1x_%d_%d_%d_%.1f", i, j, centr, etaBegin), hist1->GetYaxis()->FindBin(-etaMax+0.01), hist1->GetYaxis()->FindBin(etaMax-0.01));
373   proj1x->GetXaxis()->SetTitleOffset(1);
374   proj1x->Scale(hist1->GetYaxis()->GetBinWidth(1));
375
376   proj1y = ((TH2*) hist1)->ProjectionY(Form("proj1y_%d_%d_%d_%.1f", i, j, centr, etaBegin), hist1->GetXaxis()->FindBin(-0.5), hist1->GetXaxis()->FindBin(0.5));
377   proj1y->Scale(hist1->GetXaxis()->GetBinWidth(1));
378   proj1y->GetXaxis()->SetTitleOffset(1);
379   proj1y->GetXaxis()->SetRangeUser(-1.99, 1.99);
380   Float_t etaPhiScale = 1.0 * (hist1->GetXaxis()->FindBin(0.5) - hist1->GetXaxis()->FindBin(-0.5) + 1) / (hist1->GetYaxis()->FindBin(etaMax-0.01) - hist1->GetYaxis()->FindBin(-etaMax+0.01) + 1);
381   
382   if (gStudySystematic == 20)
383   {
384     Printf(">>>>>>>>>>>> Applying non-closure systematics <<<<<<<<<<<<");
385     file2 = TFile::Open("non_closure.root");
386     non_closure = (TH1*) file2->Get(Form("non_closure_all_%d_%d_%d", i, j, 0));
387     for (Int_t bin=1; bin<=non_closure->GetNbinsX(); bin++)
388       non_closure->SetBinError(bin, 0);
389     
390     proj1x->Multiply(non_closure);  
391   }
392   
393   Float_t zyam = 0;
394   if (0)
395   {  
396     clone = (TH1*) proj1x->Clone();
397     clone->Fit("pol0", "0", "", TMath::Pi()/2 - 0.2, TMath::Pi()/2);
398     zyam = clone->GetFunction("pol0")->GetParameter(0);
399   }
400   else
401   {
402 //     zyam = (proj1x->GetBinContent(proj1x->FindBin(TMath::Pi()/2)) + proj1x->GetBinContent(proj1x->FindBin(-TMath::Pi()/2))) / 2;
403 //     zyam = proj1x->GetBinContent(proj1x->FindBin(TMath::Pi()/2));
404     zyam = proj1x->GetBinContent(proj1x->FindBin(1.3));
405 //     zyam = proj1x->GetMinimum();
406   }
407     
408   proj1x->Add(new TF1("func", "-1", -100, 100), zyam);
409   proj1y->Add(new TF1("func", "-1", -100, 100), zyam * etaPhiScale);
410   
411   proj1x->Scale(1.0 / (2.0 * etaMax));
412   
413   if (etaProj != 0)
414     *etaProj = proj1y;
415   return proj1x;
416 }
417
418 TH1* GetProjectionsNew(Int_t i, Int_t j, Int_t centr, char** label, Float_t etaBegin = 1.0, TH1** etaProj = 0)
419 {
420   TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr));
421   if (!hist1)
422     return 0;
423   hist1 = (TH2*) hist1->Clone(Form("%s_%.1f", hist1->GetName(), etaBegin));
424
425   // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
426   hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
427   
428   hist1->Rebin2D(2, 1); hist1->Scale(0.5);
429   
430 //   new TCanvas; hist1->Draw("surf1");
431   
432   tokens = TString(hist1->GetTitle()).Tokenize("-");
433   centralityStr = new TString;
434   if (tokens->GetEntries() > 2)
435     *centralityStr = tokens->At(2)->GetName();
436   if (tokens->GetEntries() > 3)
437     *centralityStr = *centralityStr + "-" + tokens->At(3)->GetName();
438   *label = centralityStr->Data();
439 //   Printf("%s", label);
440   
441   if (etaBegin > 0)
442   {
443     proj1x = ((TH2*) hist1)->ProjectionX(Form("proj1x_%d_%d_%d_%.1f", i, j, centr, etaBegin), hist1->GetYaxis()->FindBin(-etaMax+0.01), hist1->GetYaxis()->FindBin(etaMax-0.01));
444     proj1x->GetXaxis()->SetTitleOffset(1);
445     proj1x->Scale(hist1->GetYaxis()->GetBinWidth(1));
446
447     proj1xR1 = ((TH2*) hist1)->ProjectionX(Form("proj2x_%d_%d_%d_%.1f", i, j, centr, etaBegin), hist1->GetYaxis()->FindBin(-etaMax+0.01), hist1->GetYaxis()->FindBin(-etaBegin-0.01));
448     proj1xR2 = ((TH2*) hist1)->ProjectionX(Form("proj3x_%d_%d_%d_%.1f", i, j, centr, etaBegin), hist1->GetYaxis()->FindBin(etaBegin+0.01), hist1->GetYaxis()->FindBin(etaMax-0.01));
449     proj1xR1->Add(proj1xR2);
450     
451     proj1xR1->GetXaxis()->SetTitleOffset(1);
452     proj1xR1->Scale(hist1->GetYaxis()->GetBinWidth(1));
453     
454     proj1xR1->Scale((1.0 * hist1->GetYaxis()->FindBin(etaMax-0.01) - hist1->GetYaxis()->FindBin(-etaMax+0.01) + 1) / (hist1->GetYaxis()->FindBin(-etaBegin-0.01) - hist1->GetYaxis()->FindBin(-etaMax+0.01) + 1 + hist1->GetYaxis()->FindBin(etaMax-0.01) - hist1->GetYaxis()->FindBin(etaBegin+0.01) + 1));
455     
456     // mirror
457     for (Int_t i=1; i<=proj1xR1->GetNbinsX()/2; i++)
458     {
459   //     Printf("%d -> %d", i, etaGap->GetNbinsX()+1-i);
460       proj1xR1->SetBinContent(proj1xR1->GetNbinsX()+1-i, proj1xR1->GetBinContent(i));
461       proj1xR1->SetBinError(proj1xR1->GetNbinsX()+1-i, proj1xR1->GetBinError(i));
462     }
463     
464     proj1x->Add(proj1xR1, -1);
465     
466     proj1x->Scale(1.0 / (2.0 * etaMax));
467   
468 //     new TCanvas; proj1xR1->DrawCopy();
469   }
470   else
471   {
472     proj1x = ((TH2*) hist1)->ProjectionX(Form("proj1x_%d_%d_%d_%.1f", i, j, centr, etaBegin), hist1->GetYaxis()->FindBin(-etaMax+0.01), hist1->GetYaxis()->FindBin(etaMax-0.01));
473     proj1x->GetXaxis()->SetTitleOffset(1);
474     proj1x->Scale(hist1->GetYaxis()->GetBinWidth(1));
475
476     proj1x->Scale(1.0 / (2.0 * etaMax));
477   }
478
479   if (gStudySystematic == 20)
480   {
481     Printf(">>>>>>>>>>>> Applying non-closure systematics <<<<<<<<<<<<");
482     file2 = TFile::Open("non_closure.root");
483     non_closure = (TH1*) file2->Get(Form("non_closure_all_%d_%d_%d", i, j, 0));
484     for (Int_t bin=1; bin<=non_closure->GetNbinsX(); bin++)
485       non_closure->SetBinError(bin, 0);
486     
487     proj1x->Multiply(non_closure);  
488   }
489   
490   Float_t zyam = 0;
491   if (0)
492   {  
493     clone = (TH1*) proj1x->Clone();
494     clone->Fit("pol0", "0", "", TMath::Pi()/2 - 0.2, TMath::Pi()/2);
495     zyam = clone->GetFunction("pol0")->GetParameter(0);
496   }
497   else
498   {
499 //     zyam = (proj1x->GetBinContent(proj1x->FindBin(TMath::Pi()/2)) + proj1x->GetBinContent(proj1x->FindBin(-TMath::Pi()/2))) / 2;
500 //     zyam = proj1x->GetBinContent(proj1x->FindBin(TMath::Pi()/2));
501     zyam = proj1x->GetBinContent(proj1x->FindBin(1.3));
502 //     zyam = proj1x->GetMinimum();
503   }
504     
505   proj1x->Add(new TF1("func", "-1", -100, 100), zyam);
506   
507   return proj1x;
508 }
509
510 void DrawEtaDep(const char* fileName, Int_t i, Int_t j, Int_t centr)
511 {
512   c = new TCanvas;
513   gPad->SetGridx();
514   gPad->SetGridy();
515   
516   TFile::Open(fileName);
517
518   legend = new TLegend(0.8, 0.6, 0.99, 0.99);
519   legend->SetFillColor(0);
520
521   for (Int_t n=0; n<4; n++)
522   {
523     Float_t eta = 0.8 + n * 0.2;
524     const char* label = 0;
525     TH1* hist = GetProjections(i, j, centr, &label, eta);
526     hist->SetStats(0);  
527     hist->SetLineColor(n+1);
528     hist->Draw((n == 0) ? "" : "SAME");
529     legend->AddEntry(hist, Form("|#Delta#eta| > %.1f", eta));
530   }
531   legend->Draw();
532 }
533
534 void DrawProjectionsTim(const char* fileName, Int_t i, Int_t j, Float_t eta = 1.0, Bool_t etaPhi = kFALSE)
535 {
536     gStyle->SetErrorX(0.0);
537     c = new TCanvas;
538 //    gPad->SetGridx();
539 //    gPad->SetGridy();
540     gPad->SetTopMargin(0.025);
541     gPad->SetRightMargin(0.01);
542     gPad->SetBottomMargin(0.15);
543     gPad->SetLeftMargin(0.15);
544
545     TFile::Open(fileName);
546     
547 //    TLegend *legend = new TLegend(0.45, 0.60, 0.65, 0.90);
548     TLegend *legend = new TLegend(0.45, 0.45, 0.65, 0.90);
549     legend->SetFillColor(0);
550     legend->SetBorderSize(0);
551     legend->SetTextFont(62);
552     legend->SetTextSize(0.04);
553     legend->SetFillStyle(0);
554     
555 //    TLegend *legend2 = new TLegend(0.65, 0.75, 0.85, 0.90);
556 //    legend2->SetFillColor(0);
557 //    legend2->SetBorderSize(0);
558 //    legend2->SetTextFont(62);
559 //    legend2->SetTextSize(0.04);
560 //    legend2->SetFillStyle(0);
561     
562 //     Int_t colors[] = { 1, 2, 1, 4, 6, 2 };
563     Int_t colors[] = { kRed+1, kOrange+7, kBlack, kGreen+2, kAzure+2, kBlack };
564     Int_t markers[] = { 20, 21, 1, 22, 23, 1 };
565     
566     TH1* first = 0;
567     
568     Float_t min = 100;
569     Float_t max = -100;
570     
571     Int_t centSeq[] = { 0, 1, 3, 4, 2, 5 };
572     
573     for (Int_t otcentr=0; otcentr<6; otcentr++)
574     {
575         Int_t centr = centSeq[otcentr];
576         /*    if (centr >= 5)
577          continue;*/
578         const char* label = 0;
579         TH1* etaHist = 0;
580         TH1* hist = GetProjectionsNew(i, j, centr, &label, (centr == 2 || centr == 5 || centr == 4) ? -1 : eta, &etaHist);
581         if (!hist)
582             continue;
583         if (etaPhi)
584             hist = etaHist;
585         hist->SetStats(0);
586         hist->GetXaxis()->CenterTitle();
587         hist->GetXaxis()->SetLabelSize(0.05);
588         hist->GetXaxis()->SetTitleSize(0.05);
589         hist->GetXaxis()->SetTitleOffset(1.1);
590         hist->GetXaxis()->SetTitle("#Delta#varphi (rad)");
591         hist->GetYaxis()->SetNdivisions(506);
592         hist->GetYaxis()->CenterTitle();
593         hist->GetYaxis()->SetLabelSize(0.05);
594         hist->GetYaxis()->SetTitleSize(0.05);
595         hist->GetYaxis()->SetTitleOffset(1.35);
596         
597         if (centr == 2 || centr == 5)
598         {
599           hist->SetLineWidth(2);
600           if (centr == 2)
601             hist->SetLineStyle(2);
602         }
603           
604 //         hist->GetYaxis()->SetTitle("#frac{1}{#it{N}_{trig}}#frac{d#it{N}_{assoc}}{d#Delta#varphi} (rad^{-1})");
605         hist->GetYaxis()->SetTitle("1/#it{N}_{trig} d#it{N}_{assoc}/d#Delta#varphi per #Delta#eta - const (rad^{-1})");
606
607         hist->SetLineColor(colors[centr]);
608         hist->SetMarkerColor(colors[centr]);
609         hist->SetMarkerStyle(markers[centr]);
610         if (etaPhi)
611             hist->GetXaxis()->SetRangeUser(-1.79, 1.79);
612         c->cd();
613         
614         // -----
615         tokens = TString(hist->GetTitle()).Tokenize("-");
616         sPtTRange = new TString;
617         *sPtTRange = tokens->At(0)->GetName();
618         *sPtTRange = *sPtTRange + "GeV/#it{c}";
619         sPtTRange->ReplaceAll(".0", "");
620         sPtARange = new TString;
621         *sPtARange = tokens->At(1)->GetName();
622         *sPtARange = *sPtARange + "GeV/#it{c}";
623         sPtARange->ReplaceAll(".00", "");
624         sPtARange->ReplaceAll(".0", "");
625         sPtARange->ReplaceAll(" 1", "1");
626         sPtTRange->ReplaceAll("p_", "#it{p}_");
627         sPtARange->ReplaceAll("p_", "#it{p}_");
628         cout << sPtTRange->Data() << endl;
629         cout << sPtARange->Data() << endl;
630         hist->SetTitle("");
631         // -----
632         
633         if (centr == 0)
634           hist->Draw("");
635         else if (centr == 2 || centr == 5)
636           hist->Draw("HISTE SAME");
637         else
638           hist->Draw("SAME");
639         
640         min = TMath::Min(min, hist->GetMinimum());
641         max = TMath::Max(max, hist->GetMaximum());
642         if (!first)
643             first = hist;
644         if (centr==2) legend->AddEntry(hist,"pp 2.76 TeV","l");
645         else if (centr==5) legend->AddEntry(hist,"pp 7 TeV","l");
646         else legend->AddEntry(hist, label, "p");
647         //     break;
648     }
649 //    first->GetYaxis()->SetRangeUser(min * 1.1, max * 1.1);
650 //     first->GetYaxis()->SetRangeUser(min * 1.1, 0.67/3.6);
651     first->GetYaxis()->SetRangeUser(-0.009, 0.2);
652     cout << max*1.1 << endl;
653     cout << 0.925*first->GetMaximum() << endl;
654     legend->Draw();
655 //    legend2->Draw();
656     TLine * li = new TLine(first->GetXaxis()->GetXmin(),0.0,first->GetXaxis()->GetXmax(),0.0);
657     li->SetLineStyle(kDashed);
658     li->SetLineColor(kBlack);
659     li->Draw("same");
660
661     TLatex * tex_Pbp = new TLatex(0.65,0.945*first->GetMaximum(),"p-Pb #sqrt{#it{s}_{_{NN}}} = 5.02 TeV");
662     tex_Pbp->SetTextFont(62);
663     tex_Pbp->SetTextSize(0.05);
664     tex_Pbp->Draw();
665     
666     TLatex * tex_statu = new TLatex(2.3,0.65*first->GetMaximum(),"stat. uncertainties only");
667     tex_statu->SetTextFont(62);
668     tex_statu->SetTextSize(0.04);
669     tex_statu->Draw();
670     
671     if (eta > 0)
672     {
673       TLatex * tex_statu = new TLatex(2.3,0.56*first->GetMaximum(),"ridge subtracted");
674       tex_statu->SetTextFont(62);
675       tex_statu->SetTextSize(0.04);
676       tex_statu->Draw();
677     }
678
679     TLatex * tex_PtT = new TLatex(2.3,0.85*first->GetMaximum(),sPtTRange->Data());
680     tex_PtT->SetTextFont(62);
681     tex_PtT->SetTextSize(0.04);
682     tex_PtT->Draw();
683     
684     TLatex * tex_PtA = new TLatex(2.3,0.75*first->GetMaximum(),sPtARange->Data());
685     tex_PtA->SetTextFont(62);
686     tex_PtA->SetTextSize(0.04);
687     tex_PtA->Draw();
688     
689 //     box = new TBox(-1.4, 0.185, -1.3, 0.195);
690 //     box->SetFillColor(kGray);
691 //     box->Draw();
692
693     if (i == 2 && j == 2 && eta < 0)
694     {
695       c->SaveAs("fig2.eps");
696       c->SaveAs("fig2.png");
697     }
698     else if (i == 2 && j == 2 && eta > 0)
699     {
700       c->SaveAs("fig5.eps");
701       c->SaveAs("fig5.png");
702     }
703     
704 //     c->SaveAs(Form("%s_%d_%d.png", (etaPhi) ? "eta" : "phi", i, j));
705 //     c->SaveAs(Form("%s_%d_%d.eps", (etaPhi) ? "eta" : "phi", i, j));
706     
707     return;
708     
709     c = new TCanvas;
710 //    gPad->SetGridx();
711 //    gPad->SetGridy();
712     
713     ppHist = GetProjections(i, j, 2, &label);
714     
715     for (Int_t centr=0; centr<6; centr++)
716     {
717         if (centr >= 4 || centr == 2)
718             continue;
719         const char* label = 0;
720         TH1* hist = GetProjections(i, j, centr, &label);
721         hist->SetStats(0);
722         hist->SetLineColor(centr+1);
723         hist->Divide(ppHist);
724         hist->GetYaxis()->SetRangeUser(1, 3);
725         hist->Draw((centr == 0) ? "" : "SAME");
726     }
727     legend->Draw();
728     
729     c->SaveAs(Form("phi_%d_%d_ratio.png", i, j));
730     c->SaveAs(Form("phi_%d_%d_ratio.eps", i, j));
731 }
732
733
734 void DrawProjections(const char* fileName, Int_t i, Int_t j, Float_t eta = 1.0, Bool_t etaPhi = kFALSE)
735 {
736   c = new TCanvas;
737   gPad->SetGridx();
738   gPad->SetGridy();
739   
740   TFile::Open(fileName);
741
742   legend = new TLegend(0.8, 0.6, 0.99, 0.99);
743   legend->SetFillColor(0);
744   
745   Int_t colors[] = { 1, 2, 1, 4, 6, 2 };
746   Int_t markers[] = { 1, 1, 5, 1, 1, 5 };
747   
748   TH1* first = 0;
749   
750   Float_t min = 100;
751   Float_t max = -100;
752   
753   for (Int_t centr=0; centr<6; centr++)
754   {
755 /*    if (centr >= 5)
756       continue;*/
757     const char* label = 0;
758     TH1* etaHist = 0;
759     TH1* hist = GetProjectionsNew(i, j, centr, &label, eta, &etaHist);
760     if (!hist)
761       continue;
762     if (etaPhi)
763       hist = etaHist;
764     hist->SetStats(0);
765     hist->SetLineColor(colors[centr]);
766     hist->SetMarkerColor(colors[centr]);
767     hist->SetMarkerStyle(markers[centr]);
768     if (etaPhi)
769       hist->GetXaxis()->SetRangeUser(-1.79, 1.79);
770     c->cd();
771     hist->Draw((centr == 0) ? "" : "SAME");
772     min = TMath::Min(min, hist->GetMinimum());
773     max = TMath::Max(max, hist->GetMaximum());
774     if (!first)
775       first = hist;
776     legend->AddEntry(hist, label);
777 //     break;
778   }
779   first->GetYaxis()->SetRangeUser(min * 1.1, max * 1.1);
780   legend->Draw();
781   
782   c->SaveAs(Form("%s_%d_%d.png", (etaPhi) ? "eta" : "phi", i, j));
783   
784   return;
785   
786   c = new TCanvas;
787   gPad->SetGridx();
788   gPad->SetGridy();
789   
790   ppHist = GetProjections(i, j, 2, &label);
791   
792   for (Int_t centr=0; centr<6; centr++)
793   {
794     if (centr >= 4 || centr == 2)
795       continue;
796     const char* label = 0;
797     TH1* hist = GetProjections(i, j, centr, &label);
798     hist->SetStats(0);
799     hist->SetLineColor(centr+1);
800     hist->Divide(ppHist);
801     hist->GetYaxis()->SetRangeUser(1, 3);
802     hist->Draw((centr == 0) ? "" : "SAME");
803   }
804   legend->Draw();
805   
806   c->SaveAs(Form("phi_%d_%d_ratio.png", i, j));
807 }
808
809 void CompareProjections(const char* fileName, Int_t i, Int_t j, Int_t centr)
810 {
811   TFile::Open(fileName);
812
813   const char* label = 0;
814   TH1* etaHist = 0;
815   TH1* hist = GetProjections(i, j, centr, &label, -1, &etaHist);
816   if (!hist)
817     continue;
818   TH1* hist2 = GetProjections(i, j, centr, &label, 1.2, &etaHist);
819
820   c = new TCanvas;
821   gPad->SetGridx();
822   gPad->SetGridy();
823   
824   hist->SetStats(0);
825   hist->SetMarkerStyle(20);
826   hist->Draw("");
827
828   hist2->SetMarkerStyle(21);
829   hist2->SetMarkerColor(2);
830   hist2->SetLineColor(2);
831   hist2->Draw("SAME");
832 }
833
834 void DrawProjectionsAll(const char* fileName, Float_t eta = 1, Bool_t etaPhi = kFALSE)
835 {
836   DrawProjections(fileName, 0, 1, eta, etaPhi);
837   DrawProjections(fileName, 0, 2, eta, etaPhi);
838   DrawProjections(fileName, 1, 1, eta, etaPhi);
839   DrawProjections(fileName, 1, 2, eta, etaPhi);
840   DrawProjections(fileName, 1, 3, eta, etaPhi);
841   DrawProjections(fileName, 2, 3, eta, etaPhi);
842 }
843
844 void DrawProjectionsRidge(const char* fileName, Int_t i, Int_t j, Int_t centr1, Int_t centr2)
845 {
846   Float_t etaLimit = 1.0;
847   Float_t outerLimit = 1.6;
848   
849   TFile::Open(fileName);
850   
851   TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr1));
852   TH2* hist2 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr2));
853   
854   // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
855   hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
856   hist2->Scale(1.0 / hist2->GetYaxis()->GetBinWidth(1));
857
858   SubtractEtaGapNS(hist1, etaLimit, outerLimit, kFALSE);
859   
860   proj1y = ((TH2*) hist1)->ProjectionY("proj1y", hist1->GetXaxis()->FindBin(-0.5), hist1->GetXaxis()->FindBin(0.5));
861   proj1x = ((TH2*) hist1)->ProjectionX("proj1x", hist1->GetYaxis()->FindBin(-1.79), hist1->GetYaxis()->FindBin(1.79));
862   
863   Float_t etaPhiScale = 1.0 * (hist1->GetXaxis()->FindBin(0.5) - hist1->GetXaxis()->FindBin(-0.5) + 1) / (hist1->GetYaxis()->FindBin(1.79) - hist1->GetYaxis()->FindBin(-1.79) + 1);
864   Printf("%f", etaPhiScale);
865   
866   proj1y->GetXaxis()->SetTitleOffset(1);
867   proj1x->GetXaxis()->SetTitleOffset(1);
868   
869   clone = (TH1*) proj1x->Clone();
870   clone->Fit("pol0", "0", "", 1.2, 1.8);
871   Float_t zyam = clone->GetFunction("pol0")->GetParameter(0);
872
873   proj1x->Add(new TF1("func", "-1", -100, 100), zyam);
874   proj1y->Add(new TF1("func", "-1", -100, 100), zyam * etaPhiScale);
875   
876   new TCanvas("c", "c", 800, 800);
877   gPad->SetLeftMargin(0.15);
878 //   hist1->SetTitle("");
879
880   hist1->GetYaxis()->SetRangeUser(-1.79, 1.79);
881   hist1->GetXaxis()->SetTitleOffset(1.5);
882   hist1->GetYaxis()->SetTitleOffset(2);
883   hist1->SetStats(kFALSE);
884   hist1->DrawCopy("SURF1");
885   
886   proj2y = ((TH2*) hist2)->ProjectionY("proj2y", hist1->GetXaxis()->FindBin(-0.5), hist1->GetXaxis()->FindBin(0.5));
887   proj2x = ((TH2*) hist2)->ProjectionX("proj2x", hist1->GetYaxis()->FindBin(-1.79), hist1->GetYaxis()->FindBin(1.79));
888
889 //   proj2y->Scale(1.0 / (hist1->GetXaxis()->FindBin(0.5) - hist1->GetXaxis()->FindBin(-0.5) + 1));
890 //   proj2x->Scale(1.0 / (hist1->GetYaxis()->FindBin(1.79) - hist1->GetYaxis()->FindBin(-1.79) + 1));
891  
892   clone = (TH1*) proj2x->Clone();
893   clone->Fit("pol0", "0", "", 1.2, 1.8);
894   zyam = clone->GetFunction("pol0")->GetParameter(0);
895
896   proj2x->Add(new TF1("func", "-1", -100, 100), zyam);
897   proj2y->Add(new TF1("func", "-1", -100, 100), zyam * etaPhiScale);
898
899   proj2y->SetLineColor(2); proj2x->SetLineColor(2);
900   
901   new TCanvas; proj1y->Draw(); proj2y->Draw("SAME"); gPad->SetGridx(); gPad->SetGridy();
902   new TCanvas; proj1x->Draw(); proj2x->Draw("SAME"); gPad->SetGridx(); gPad->SetGridy();
903   
904   new TCanvas("c2", "c2", 800, 800);
905   gPad->SetLeftMargin(0.15);
906 //   hist2->SetTitle("");
907   hist2->GetYaxis()->SetRangeUser(-1.79, 1.79);
908   hist2->GetXaxis()->SetTitleOffset(1.5);
909   hist2->GetYaxis()->SetTitleOffset(2);
910   hist2->SetStats(kFALSE);
911   hist2->DrawCopy("SURF1");
912 }
913
914 void CalculateIAA(Int_t i, Int_t j, Int_t centr, Float_t etaBegin, Double_t& nsPeak, Double_t& nsPeakE, Double_t& asPeak, Double_t& asPeakE, Double_t& nsRidge, Double_t& nsRidgeE, char** label, Bool_t subtractRidge)
915 {
916   TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr));
917   if (!hist)
918   {
919     nsPeak = 0;
920     nsPeakE = 0;
921     asPeak = 0;
922     asPeakE = 0;
923     nsRidge = 0;
924     nsRidgeE = 0;
925     return;
926   }
927   
928 //   new TCanvas; hist->Draw("COLZ");
929   
930   TH2* hist1 = (TH2*) hist->Clone();
931   hist1->Reset();
932   
933   // copy to one quadrant
934   for (Int_t x=1; x<=hist->GetNbinsX(); x++)
935     for (Int_t y=1; y<=hist->GetNbinsY(); y++)
936     {
937       Int_t xTarget = x;
938       if (hist->GetXaxis()->GetBinCenter(x) < 0)
939         xTarget = hist->GetXaxis()->FindBin(-1.0 * hist->GetXaxis()->GetBinCenter(x));
940       else if (hist->GetXaxis()->GetBinCenter(x) > TMath::Pi())
941         xTarget = hist->GetXaxis()->FindBin(TMath::TwoPi() - hist->GetXaxis()->GetBinCenter(x));
942       
943       Int_t yTarget = y;
944       if (hist->GetYaxis()->GetBinCenter(y) < 0)
945         yTarget = hist->GetYaxis()->FindBin(-1.0 * hist->GetYaxis()->GetBinCenter(y));
946       
947 //       Printf("%d %d --> %d %d", x, y, xTarget, yTarget);
948       
949       Float_t value = 0;
950       Float_t error = 0;
951       value += hist->GetBinContent(x, y);
952       error += hist->GetBinError(x, y) * hist->GetBinError(x, y);
953
954       value += hist1->GetBinContent(xTarget, yTarget);
955       error += hist1->GetBinError(xTarget, yTarget) * hist1->GetBinError(xTarget, yTarget);
956
957       error = TMath::Sqrt(error);
958       
959       hist1->SetBinContent(xTarget, yTarget, value);
960       hist1->SetBinError(xTarget, yTarget, error);
961     }
962   
963   // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
964   hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
965   new TCanvas; hist1->Draw("COLZ");
966 //   new TCanvas; hist1->Draw("SURF1");
967   
968   tokens = TString(hist1->GetTitle()).Tokenize("-");
969   centralityStr = new TString;
970   if (tokens->GetEntries() > 1)
971   {
972     *centralityStr = tokens->At(0)->GetName();
973     *centralityStr = *centralityStr + "-" + tokens->At(1)->GetName();
974   }
975   *label = centralityStr->Data();
976
977   Int_t phi1 = hist1->GetXaxis()->FindBin(0.0001);
978   Int_t phi3 = hist1->GetXaxis()->FindBin(TMath::Pi()/2 - 0.3);
979   Int_t phi2 = phi3 - 1;
980   Int_t phi4 = hist1->GetXaxis()->FindBin(TMath::Pi()/2 - 0.1);
981   Int_t phi5 = phi4 + 1;
982   Int_t phi6 = hist1->GetXaxis()->FindBin(TMath::Pi() - 0.0001);
983 //   Printf("phi = %d %d %d %d %d %d", phi1, phi2, phi3, phi4, phi5, phi6);
984   
985   Int_t eta1 = hist1->GetYaxis()->FindBin(0.0001);
986   Int_t eta2 = hist1->GetYaxis()->FindBin(etaBegin - 0.0001);
987   Int_t eta3 = eta2 + 1;
988   Int_t eta4 = hist1->GetYaxis()->FindBin(etaMax - 0.0001);
989 //   Printf("eta = %d %d %d %d", eta1, eta2, eta3, eta4);
990   
991   Double_t zyamYieldE, zyamYieldE1, zyamYieldE2;
992   nsPeak              = hist1->IntegralAndError(phi1, phi2, eta1, eta2, nsPeakE, "width");
993   nsRidge             = hist1->IntegralAndError(phi1, phi2, eta3, eta4, nsRidgeE, "width");
994   asPeak              = hist1->IntegralAndError(phi5, phi6, eta1, eta4, asPeakE, "width");
995   Double_t zyamYield  = hist1->IntegralAndError(phi3, phi4, eta1, eta4, zyamYieldE, "width");
996   Double_t zyamYield1 = hist1->IntegralAndError(phi3, phi4, eta1, eta2, zyamYieldE1, "width");
997   Double_t zyamYield2 = hist1->IntegralAndError(phi3, phi4, eta3, eta4, zyamYieldE2, "width");
998   
999   // factor 4 from folding to one quadrant above
1000   Double_t zyamYield2Density = zyamYield2 / (hist1->GetXaxis()->GetBinUpEdge(phi4) - hist1->GetXaxis()->GetBinLowEdge(phi3)) / (hist1->GetYaxis()->GetBinUpEdge(eta4) - hist1->GetYaxis()->GetBinLowEdge(eta3)) / 4;
1001   Double_t zyamYieldE2Density = zyamYieldE2 / (hist1->GetXaxis()->GetBinUpEdge(phi4) - hist1->GetXaxis()->GetBinLowEdge(phi3)) / (hist1->GetYaxis()->GetBinUpEdge(eta4) - hist1->GetYaxis()->GetBinLowEdge(eta3)) / 4;
1002   
1003   Double_t nsZyamScaling = 1.0 * (phi2 - phi1 + 1) / (phi4 - phi3 + 1);
1004   Double_t asZyamScaling = 1.0 * (phi6 - phi5 + 1) / (phi4 - phi3 + 1);
1005 //   Printf("%f %f", nsZyamScaling, asZyamScaling);
1006   
1007   nsPeak -= zyamYield1 * nsZyamScaling;
1008   nsPeakE = TMath::Sqrt(zyamYieldE1 * zyamYieldE1 * nsZyamScaling * nsZyamScaling + nsPeakE * nsPeakE);
1009   
1010   nsRidge -= zyamYield2 * nsZyamScaling;
1011   nsRidgeE = TMath::Sqrt(zyamYieldE2 * zyamYieldE2 * nsZyamScaling * nsZyamScaling + nsRidgeE * nsRidgeE);
1012
1013   asPeak -= zyamYield * asZyamScaling;
1014   asPeakE = TMath::Sqrt(zyamYieldE * zyamYieldE * asZyamScaling * asZyamScaling + asPeakE * asPeakE);
1015   
1016   if (subtractRidge)
1017   {
1018     Double_t nsRidgeScaling = 1.0 * (eta2 - eta1 + 1) / (eta4 - eta3 + 1);
1019     Double_t asRidgeScaling = 1.0 * (eta4 - eta1 + 1) / (eta4 - eta3 + 1);
1020     
1021     nsPeak -= nsRidge * nsRidgeScaling;
1022     asPeak -= nsRidge * asRidgeScaling;
1023     nsPeakE = TMath::Sqrt(nsRidgeE * nsRidgeE * nsRidgeScaling * nsRidgeScaling + nsPeakE * nsPeakE);
1024     asPeakE = TMath::Sqrt(nsRidgeE * nsRidgeE * asRidgeScaling * asRidgeScaling + asPeakE * asPeakE);
1025   }
1026
1027   nsRidge /= (etaMax - etaBegin) * 2;
1028   nsRidgeE /= (etaMax - etaBegin) * 2;
1029   
1030   Printf("Peak yields (%d %d %d): %f +- %f; %f +- %f; %f +- %f; %f +- %f", i, j, centr, nsPeak, nsPeakE, asPeak, asPeakE, nsRidge, nsRidgeE, zyamYield2Density, zyamYieldE2Density);
1031
1032 //   Printf("");
1033 }
1034
1035 void PlotIAA(const char* fileName)
1036 {
1037   Int_t colors[] = { 1, 2, 3, 4, 6, 7 };
1038   Int_t markers[] = { 20, 21, 22, 23, 24, 25 };
1039   
1040   if (0)
1041   {
1042     Int_t n = 6;
1043     Int_t is[] = { 0, 0, 1, 1, 1, 2 };
1044     Int_t js[] = { 1, 2, 1, 2, 3, 3 };
1045   
1046     Int_t centralityBins = 6;
1047     Float_t centralityX[] = { 10, 30, 110, 50, 80, 120 };
1048     Float_t centralityEX[] = { 10, 10, 0, 10, 20, 0 };
1049   }
1050   else if (0)
1051   {
1052     Int_t n = 3;
1053     Int_t is[] = { 0, 1, 2, 3 };
1054     Int_t js[] = { 1, 2, 3, 4 };
1055
1056     Int_t centralityBins = 4;
1057     Float_t centralityX[] = { 1.5, 6.5, 30, 75 };
1058     Float_t centralityEX[] = { 1.5, 2.5, 20, 25 };
1059   }
1060   else
1061   {
1062     Int_t n = 1;
1063     Int_t is[] = { 0 };
1064     Int_t js[] = { 1 };
1065   
1066     Int_t centralityBins = 6;
1067     Float_t centralityX[] = { 10, 30, 110, 50, 80, 120 };
1068     Float_t centralityEX[] = { 10, 10, 0, 10, 20, 0 };
1069   }
1070   
1071   const char* graphTitles[] = { "NS Yield", "AS Yield", "NS Ridge" };
1072
1073   TCanvas* canvas[3];
1074   for (Int_t ci=0; ci<3; ci++)
1075   {
1076     canvas[ci] = new TCanvas;
1077     gPad->SetGridx();
1078     gPad->SetGridy();
1079   }
1080   
1081   TFile::Open(fileName);
1082
1083   legend = new TLegend(0.4, 0.6, 0.99, 0.99);
1084   legend->SetFillColor(0);
1085
1086   for (Int_t i=0; i<n; i++)
1087   {
1088     TGraphErrors* graph[5];
1089     for (Int_t ci=0; ci<5; ci++)
1090       graph[ci] = new TGraphErrors;
1091
1092     char* label = 0;
1093     for (Int_t c=0; c<centralityBins; c++)
1094     {
1095       Double_t nsYield, nsError, asYield, asError, nsRidge, nsRidgeE;
1096       CalculateIAA(is[i], js[i], c, 1.2, nsYield, nsError, asYield, asError, nsRidge, nsRidgeE, &label, kTRUE);
1097
1098       AddPoint(graph[0], centralityX[c], nsYield, centralityEX[c], nsError);
1099       AddPoint(graph[1], centralityX[c], asYield, centralityEX[c], asError);
1100       AddPoint(graph[2], centralityX[c], nsRidge, centralityEX[c], nsRidgeE);
1101
1102 //       if (c != 2 && c != 5)
1103       {
1104         CalculateIAA(is[i], js[i], c, 1.2, nsYield, nsError, asYield, asError, nsRidge, nsRidgeE, &label, kFALSE);
1105         AddPoint(graph[3], centralityX[c], nsYield, 0, nsError);
1106         AddPoint(graph[4], centralityX[c], asYield, 0, asError);
1107       }
1108     }
1109
1110     for (Int_t ci=0; ci<3; ci++)
1111     {
1112       canvas[ci]->cd();
1113       graph[ci]->SetMarkerStyle(markers[i]);
1114       graph[ci]->SetMarkerColor(colors[i]);
1115       graph[ci]->SetLineColor(colors[i]);
1116       graph[ci]->Draw((i == 0) ? "AP" : "PSAME");
1117       graph[ci]->GetYaxis()->SetTitle(graphTitles[ci]);
1118       graph[ci]->GetYaxis()->SetRangeUser(0, 1);
1119     }
1120     for (Int_t ci=3; ci<5; ci++)
1121     {
1122       canvas[ci-3]->cd();
1123       graph[ci]->SetLineColor(colors[i]);
1124       graph[ci]->Sort();
1125       graph[ci]->Draw("LSAME");
1126     }
1127     
1128     legend->AddEntry(graph[0], label, "P");
1129   }
1130   
1131   canvas[0]->cd();
1132   legend->Draw();
1133   canvas[0]->SaveAs("ns_yield.png");
1134
1135   canvas[1]->cd();
1136   legend->Draw();
1137   canvas[1]->SaveAs("as_yield.png");
1138
1139   canvas[2]->cd();
1140   legend->Draw();
1141   canvas[2]->SaveAs("ns_ridge.png");
1142 }
1143
1144 void Draw2D(const char* fileName, Int_t i, Int_t j, Int_t centr)
1145 {
1146   TFile::Open(fileName);
1147   
1148   TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr));
1149   if (!hist1)
1150     return 0;
1151   // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
1152   hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
1153   
1154 //   hist1->Rebin2D(2, 2); hist1->Scale(0.25);
1155   
1156   hist1->GetYaxis()->SetRangeUser(-1.79, 1.79);
1157   new TCanvas; 
1158   hist1->DrawCopy("SURF1");
1159 }
1160
1161 void CorrelationSubtraction(const char* fileName, Int_t i, Int_t j, Int_t centr)
1162 {
1163   CreateGraphStructure();
1164
1165   CorrelationSubtraction(fileName, i, j, centr, graphs[0], 0);
1166 }
1167
1168 void CorrelationSubtractionAll(const char* fileName, Int_t centr = 0)
1169 {
1170   Int_t n = 6;
1171   Int_t is[] = { 0, 1, 1, 2, 2, 2, 3 };
1172   Int_t js[] = { 1, 1, 2, 1, 2, 3, 3 };
1173
1174   CreateGraphStructure();
1175
1176   for (Int_t i=0; i<n; i++)
1177     CorrelationSubtraction(fileName, is[i], js[i], centr, graphs[0], 0);
1178 }
1179
1180 Int_t gStudySystematic = 0; // 10 = exclusion zone to 0.5; 11 = exclusion off; 12 = exclusion 0.8 and mirror to AS; 13 = scale peripheral; 14 = exclusion 1.2; 20 = non closure; 30 = baseline; 40 = track cuts; 50/51 = cms comparison; 60 = other peripheral bin
1181
1182 void CorrelationSubtractionHistogram(const char* fileName, Bool_t centralityAxis = kTRUE, const char* outputFileName = "graphs.root")
1183 {
1184 //   Int_t n = 10; //15;
1185 //   // sorted by pT,T
1186 //   Int_t is[] =    { 0, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4 };
1187 //   Int_t js[] =    { 1, 1, 2, 1, 2, 3, 1, 2, 3, 4, 1, 2, 3, 4, 5 };
1188 //   Bool_t symm[] = { 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1 };
1189
1190   // for extraction when trigger pT is from whole range
1191   Int_t n = 5;
1192   Int_t is[] =    { 0, 0, 0, 0, 0 };
1193   Int_t js[] =    { 1, 2, 3, 4, 6 };
1194   Bool_t symm[] = { 1, 1, 1, 1, 1 };
1195
1196   // sorted by pT,a
1197 //   Int_t n = 10; // 15
1198 /*  Int_t is[] =    { 0, 1, 2, 3, 1, 2, 3, 2, 3, 3 };
1199   Int_t js[] =    { 1, 1, 1, 1, 2, 2, 2, 3, 3, 4 };
1200   Bool_t symm[] = { 1, 0, 0, 0, 1, 0, 0, 1, 0, 1 };*/
1201 /*  Int_t is[] =    { 0, 1, 2, 3, 4, 1, 2, 3, 4, 2, 3, 4, 3, 4, 4 };
1202   Int_t js[] =    { 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 5 };
1203   Bool_t symm[] = { 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1 };*/
1204
1205   //   Int_t n = 6;
1206 //   Int_t is[] = { 0, 0, 1, 1, 1, 2 };
1207 //   Int_t js[] = { 1, 2, 1, 2, 3, 3 };
1208
1209   Int_t centralityBins = 4;
1210
1211   Int_t colors[] = { 1, 2, 3, 4, 6, 7 };
1212   Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34 };
1213
1214   CreateGraphStructure();
1215   
1216   Float_t yMax[] = { 0.1, 2, 0.4, 4, 1.5, 50, 50, 4 };
1217
1218   TCanvas* canvas[8];
1219   for (Int_t ci=0; ci<8; ci++)
1220   {
1221     canvas[ci] = new TCanvas;
1222     gPad->SetGridx();
1223     gPad->SetGridy();
1224     dummy = new TH2F("dummy", Form(";%s;%s", (centralityAxis) ? "Event class" : "<mult>", graphTitles[(ci < 3) ? ci*2 : ci+3]), 100, 0, 60, 100, 0, yMax[ci]);
1225     if (ci == 0)
1226       dummy->GetYaxis()->SetTitle("Yield");
1227     if (ci == 1)
1228       dummy->GetYaxis()->SetTitle("Width");
1229     if (ci == 2)
1230       dummy->GetYaxis()->SetTitle("v2 , v3");
1231     dummy->SetStats(0);
1232     dummy->Draw();
1233   }
1234   
1235   legend = new TLegend(0.4, 0.6, 0.99, 0.99);
1236   legend->SetFillColor(0);
1237
1238   legend2 = new TLegend(0.4, 0.6, 0.99, 0.99);
1239   legend2->SetFillColor(0);
1240
1241   for (Int_t i=0; i<n; i++)
1242   {
1243     for (Int_t c=0; c<centralityBins; c++)
1244     {
1245       if (c == 2)
1246         continue;
1247       
1248       Double_t nsYield, nsYieldE, asYield, asYieldE, nsWidth, nsWidthE, asWidth, asWidthE, v2, v2E, v3, v3E, centrality;
1249       CorrelationSubtraction(fileName, is[i], js[i], c, graphs[i], 1.5 * (-n/2 + i), kTRUE, symm[i], centralityAxis);
1250     }
1251
1252     for (Int_t ci=0; ci<6; ci++)
1253     {
1254       canvas[ci/2]->cd();
1255
1256       graphs[i][ci]->SetMarkerStyle(markers[i]);
1257       graphs[i][ci]->SetMarkerColor((ci % 2 == 0) ? 1 : 2);
1258       graphs[i][ci]->SetLineColor((ci % 2 == 0) ? 1 : 2);
1259       graphs[i][ci]->GetXaxis()->SetTitle(dummy->GetXaxis()->GetTitle());
1260       graphs[i][ci]->Draw("PSAME");
1261     }
1262     
1263     for (Int_t ci=6; ci<11; ci++)
1264     {
1265       canvas[ci-3]->cd();
1266
1267       graphs[i][ci]->SetMarkerStyle(markers[i]);
1268       graphs[i][ci]->SetMarkerColor(1);
1269       graphs[i][ci]->SetLineColor(1);
1270       graphs[i][ci]->GetXaxis()->SetTitle(dummy->GetXaxis()->GetTitle());
1271       graphs[i][ci]->Draw("PSAME");
1272     }
1273
1274     legend->AddEntry(graphs[i][0], graphs[i][0]->GetTitle(), "P");
1275     if (symm[i])
1276       legend2->AddEntry(graphs[i][0], graphs[i][0]->GetTitle(), "P");
1277   }
1278   
1279   const char* fileNames[] = { "ridge_yield.png", "ridge_width.png", "v2_v3.png", "ridge_yield_ratio.png", "v3_over_v2.png", "remaining_peak.png", "remaining_jet.png", "chi2ndf.png" };
1280   
1281   for (Int_t ci=0; ci<8; ci++)
1282   {
1283     canvas[ci]->cd();
1284     if (ci == 2 || ci == 4)
1285       legend2->Draw();
1286     else
1287       legend->Draw();
1288     canvas[ci]->SaveAs(fileNames[ci]);
1289   }
1290   
1291   if (gStudySystematic != 0)
1292     Printf(">>>>>>>>>>>> WARNING: gStudySystematic set to %d", gStudySystematic);
1293
1294   WriteGraphs(outputFileName);
1295 }
1296
1297 void CorrelationSubtraction(const char* fileName, Int_t i, Int_t j, Int_t centr, TGraphErrors** graph, Float_t axisOffset, Bool_t silent = kFALSE, Bool_t symmetricpT = kFALSE, Bool_t centralityAxis = kTRUE)
1298 {
1299   TFile::Open(fileName);
1300   
1301   TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr));
1302   TH2* hist2 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, ((gStudySystematic == 60) ? 6 : 4)));
1303   TH1* refMult = (TH1*) gFile->Get("refMult");
1304
1305   hist1 = (TH2*) hist1->Clone();
1306   hist2 = (TH2*) hist2->Clone();
1307   
1308   if (!hist1 || !hist2)
1309     return 0;
1310   // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
1311   hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
1312   hist2->Scale(1.0 / hist2->GetYaxis()->GetBinWidth(1));
1313   
1314 //   hist1->Scale(1.0 / 283550); hist2->Scale(1.0 / 190000); // for 1, 2, 0
1315 //   hist1->Scale(1.0 / 283550); hist2->Scale(1.0 / 75000); // for 2, 3, 0
1316 //   hist1->Scale(1.0 / 283550); hist2->Scale(1.0 / 100000); // for 2, 2, 0
1317   
1318   tokens = TString(hist1->GetTitle()).Tokenize("-");
1319   if (tokens->GetEntries() > 1)
1320   {
1321     TString centralityStr;
1322     centralityStr = tokens->At(0)->GetName();
1323     centralityStr = centralityStr + "-" + tokens->At(1)->GetName();
1324 //     Printf("%s", centralityStr.Data());
1325     for (Int_t gID=0; gID<NHists; gID++)
1326     {
1327       TString xTitle = graph[gID]->GetXaxis()->GetTitle();
1328       TString yTitle = graph[gID]->GetYaxis()->GetTitle();
1329 //       Printf("%s", centralityStr.Data());
1330       graph[gID]->SetTitle(centralityStr);
1331       graph[gID]->GetXaxis()->SetTitle(xTitle);
1332       graph[gID]->GetYaxis()->SetTitle(yTitle);
1333 //       Printf("%s %s", yTitle.Data(), graph[gID]->GetYaxis()->GetTitle());
1334     }
1335     Double_t centrality = -1;
1336     if (tokens->GetEntries() == 4)
1337       centrality = 0.5 * (TString(tokens->At(3)->GetName()).Atoi() + TString(tokens->At(2)->GetName()).Atoi());
1338   }
1339   
1340   hist1->GetYaxis()->SetRangeUser(-1.99, 1.99); hist2->GetYaxis()->SetRangeUser(-1.99, 1.99);
1341   if (!silent)
1342   {
1343     new TCanvas; 
1344     copy = (TH2*) hist1->DrawCopy("SURF1");
1345     copy->Rebin2D(2, 2);
1346     copy->Scale(0.25);
1347     copy->GetYaxis()->SetRangeUser(-1.99, 1.99);
1348     
1349     new TCanvas; 
1350     copy = (TH2*) hist2->DrawCopy("SURF1");
1351     copy->Rebin2D(2, 2);
1352     copy->Scale(0.25);
1353     copy->GetYaxis()->SetRangeUser(-1.99, 1.99);
1354   }
1355   
1356   if (gStudySystematic == 13)
1357   {
1358     Float_t factor[10][10][10];
1359     
1360     if (1)
1361     {
1362       // VZERO, from dphi_corr_pA_121108_hybrid_corrected.root
1363       factor[0][1][0] = 0.259;
1364       factor[0][1][1] = 0.212;
1365       factor[0][1][3] = 0.150;
1366       factor[1][1][0] = 0.105;
1367       factor[1][1][1] = 0.105;
1368       factor[1][1][3] = 0.081;
1369       factor[1][2][0] = 0.209;
1370       factor[1][2][1] = 0.176;
1371       factor[1][2][3] = 0.121;
1372       factor[2][1][0] = 0.050;
1373       factor[2][1][1] = -0.005;
1374       factor[2][1][3] = 0.031;
1375       factor[2][2][0] = 0.042;
1376       factor[2][2][1] = 0.061;
1377       factor[2][2][3] = 0.022;
1378       factor[2][3][0] = 0.081;
1379       factor[2][3][1] = 0.073;
1380       factor[2][3][3] = 0.058;
1381     }
1382     else if (0)
1383     {
1384       // ZNA, from dphi_corr_pA_121029_zna.root
1385       factor[0][1][0] = 0.125;
1386       factor[0][1][1] = 0.114;
1387       factor[0][1][3] = 0.052;
1388       factor[1][1][0] = 0.052;
1389       factor[1][1][1] = 0.065;
1390       factor[1][1][3] = 0.024;
1391       factor[1][2][0] = 0.113;
1392       factor[1][2][1] = 0.102;
1393       factor[1][2][3] = 0.060;
1394       factor[2][1][0] = 0.014;
1395       factor[2][1][1] = 0.023;
1396       factor[2][1][3] = -0.017;
1397       factor[2][2][0] = 0.018;
1398       factor[2][2][1] = 0.037;
1399       factor[2][2][3] = -0.005;
1400       factor[2][3][0] = 0.068;
1401       factor[2][3][1] = 0.066;
1402       factor[2][3][3] = 0.035;
1403     }
1404     else if (0) //V0, cmsmethod
1405     {
1406       factor[0][1][0] = 0.470 +  0.221 + 0.103;
1407       factor[0][1][1] = 0.432 +  0.186 + 0.080;
1408       factor[0][1][3] = 0.339 +  0.115 + 0.039;
1409       factor[1][1][0] = 0.164 +  0.027 + 0.005;
1410       factor[1][1][1] = 0.160 +  0.025 + 0.004;
1411       factor[1][1][3] = 0.111 +  0.013 + 0.001;
1412       factor[1][2][0] = 0.458 +  0.210 + 0.097;
1413       factor[1][2][1] = 0.385 +  0.148 + 0.057;
1414       factor[1][2][3] = 0.280 +  0.078 + 0.022;
1415       factor[2][1][0] = 0.026 +  0.000 + 0.000;
1416       factor[2][1][1] = -0.018 + 0.001 + -0.000;
1417       factor[2][1][3] = -0.005 + 0.000 + -0.000;
1418       factor[2][2][0] = 0.047 +  0.003 + -0.000;
1419       factor[2][2][1] = 0.044 +  0.002 + 0.000;
1420       factor[2][2][3] = 0.013 +  0.000 + 0.000;
1421       factor[2][3][0] = 0.258 +  0.067 + 0.017;
1422       factor[2][3][1] = 0.176 +  0.031 + 0.005;
1423       factor[2][3][3] = 0.118 +  0.014 + 0.002;
1424     }
1425     else
1426     {
1427       factor[1][2][0] = 0.3;
1428     }
1429     
1430     hist2->Scale(1.0 + factor[i][j][centr]);
1431   }
1432   
1433   const Float_t etaFlat = 1.2;
1434   const Float_t centralEta = 0.5;
1435
1436   projCentral = hist1->ProjectionY(Form("%s_proj3y", hist1->GetName()), hist1->GetXaxis()->FindBin(-0.7), hist1->GetXaxis()->FindBin(0.7));
1437   projCentral->Scale(hist1->GetXaxis()->GetBinWidth(1));
1438   projCentral2 = hist2->ProjectionY(Form("%s_proj4y", hist2->GetName()), hist2->GetXaxis()->FindBin(-0.7), hist2->GetXaxis()->FindBin(0.7));
1439   projCentral2->Scale(hist2->GetXaxis()->GetBinWidth(1));
1440
1441   Float_t yEta1 = projCentral->FindBin(-etaMax+0.01);
1442   Float_t yEta2 = projCentral->FindBin(-etaFlat-0.01);
1443   Float_t yEta3 = projCentral->FindBin(-centralEta+0.01);
1444   Float_t yEta4 = projCentral->FindBin(centralEta-0.01);
1445   Float_t yEta5 = projCentral->FindBin(etaFlat+0.01);
1446   Float_t yEta6 = projCentral->FindBin(etaMax-0.01);
1447   Float_t yBaseline = (projCentral->Integral(yEta1, yEta2) + projCentral->Integral(yEta5, yEta6)) / (yEta2 - yEta1 + 1 + yEta6 - yEta5 + 1);
1448   Float_t yPeak = projCentral->Integral(yEta3, yEta4, "width") - yBaseline * (yEta4 - yEta3 + 1) * projCentral->GetBinWidth(1);
1449   Printf("Peak (unsubtracted): %f %f", yBaseline, yPeak);
1450
1451   if (!silent)
1452   {
1453     Float_t yBaseline = (projCentral->Integral(yEta1, yEta2) + projCentral->Integral(yEta5, yEta6)) / (yEta2 - yEta1 + 1 + yEta6 - yEta5 + 1);
1454     Float_t yBaseline2 = (projCentral2->Integral(yEta1, yEta2) + projCentral2->Integral(yEta5, yEta6)) / (yEta2 - yEta1 + 1 + yEta6 - yEta5 + 1);
1455     Printf("baselines: %f %f", yBaseline, yBaseline2);
1456     
1457     projCentral2->Add(new TF1("flat", "1", -5, 5), yBaseline - yBaseline2);
1458     
1459     new TCanvas; 
1460     projCentral->DrawCopy();
1461     projCentral2->SetLineColor(2);
1462     projCentral2->DrawCopy("SAME");
1463   }
1464   
1465   Float_t xValue1 = -1;
1466   Float_t xValue2 = -1;
1467   Float_t xValue1vn = -1;
1468   Float_t xValue2vn = -1;
1469   if (centralityAxis)
1470   {
1471     xValue1 = centrality + axisOffset;
1472     xValue2 = xValue1 + 0.5;
1473     xValue1vn = centrality + (Int_t) axisOffset;
1474     xValue2vn = xValue1vn + 0.5;
1475   }
1476   else
1477   {
1478     if (centrality > 0)
1479     {
1480       xValue1 = refMult->GetBinContent(refMult->GetXaxis()->FindBin(centrality));
1481       xValue2 = xValue1;
1482       xValue1vn = xValue1;
1483       xValue2vn = xValue1;
1484 //       Printf("%f %f %f", centrality, xValue1, xValue2);
1485     }
1486   }
1487   
1488   Float_t exclusion = 0.8;
1489   if (gStudySystematic == 10)
1490     exclusion = 0.5;
1491   else if (gStudySystematic == 11)
1492     exclusion = 0.0;
1493   else if (gStudySystematic == 14 || gStudySystematic == 51)
1494     exclusion = 1.2;
1495   
1496   Int_t eta1 = hist1->GetYaxis()->FindBin(-etaMax + 0.001);
1497   Int_t eta2 = hist1->GetYaxis()->FindBin(-exclusion - 0.001);
1498   Int_t eta3 = hist1->GetYaxis()->FindBin(exclusion + 0.001);
1499   Int_t eta6 = hist1->GetYaxis()->FindBin(etaMax - 0.001);
1500   Int_t phi1Z = hist1->GetXaxis()->FindBin(TMath::Pi()/2 - 0.2);
1501   Int_t phi2Z = hist1->GetXaxis()->FindBin(TMath::Pi()/2 + 0.2);
1502   
1503 //   new TCanvas; tmpProj = hist1->ProjectionX("tmp", eta1, eta6); tmpProj->Scale(hist1->GetYaxis()->GetBinWidth(1)); tmpProj->Scale(1.0 / (etaMax * 2)); tmpProj->Draw();
1504   
1505   Double_t baseLineE;
1506   Float_t baseLine = hist1->IntegralAndError(phi1Z, phi2Z, eta1, eta6, baseLineE);
1507   baseLine /= (phi2Z - phi1Z + 1) * (eta6 - eta1 + 1);
1508   baseLineE /= (phi2Z - phi1Z + 1) * (eta6 - eta1 + 1);
1509   
1510   hist1->Add(hist2, -1);
1511   
1512   // phi projection
1513   // NS
1514   proj = hist1->ProjectionX(Form("%s_proj1x", hist1->GetName()), eta1, eta2);
1515   proj2 = hist1->ProjectionX(Form("%s_proj2x", hist1->GetName()), eta3, eta6);
1516   proj->Add(proj2, 1);
1517   
1518   // AS
1519   projAS = hist1->ProjectionX(Form("%s_proj3x", hist1->GetName()), eta1, eta6);
1520
1521   // match NS and AS yield
1522   proj->Scale(1.0 * (eta6 - eta1 + 1) / (eta6 - eta3 + 1 + eta2 - eta1 + 1));
1523   
1524   // copy AS
1525   for (Int_t bin=proj->FindBin(TMath::Pi()/2+0.001); bin<=proj->GetNbinsX(); bin++)
1526   {
1527 //     Printf("%d %f", bin, projAS->GetBinContent(bin));
1528     proj->SetBinContent(bin, projAS->GetBinContent(bin));
1529     proj->SetBinError(bin, projAS->GetBinError(bin));
1530   }
1531   
1532   if (gStudySystematic == 12)
1533   {
1534     // get difference between exclusion 0.8 and 0.0, shift to AS, remove from AS
1535     
1536     projAll = hist1->ProjectionX(Form("%s_proj4x", hist1->GetName()), eta1, eta6);
1537     projAll->Add(proj, -1);
1538     if (!silent)
1539     {
1540       new TCanvas; projAll->DrawCopy();
1541     }
1542
1543     // From Tim, 13.11.12
1544     Float_t ratioNSAS[10][10][10];
1545     ratioNSAS[0][1][4] = 2.060797;
1546     ratioNSAS[1][1][4] = 1.433904;
1547     ratioNSAS[1][2][4] = 1.390280;
1548     ratioNSAS[2][1][4] = 1.013925;
1549     ratioNSAS[2][2][4] = 1.070422;
1550     ratioNSAS[2][3][4] = 1.087238;
1551     ratioNSAS[3][3][4] = 1.023756;
1552     
1553     Printf("Using NS/AS ratio %f", ratioNSAS[i][j][4]);
1554
1555     for (Int_t k=0; k<=projAll->GetNbinsX()/2; k++)
1556     {
1557   //     Printf("%d -> %d", i, projAll->GetNbinsX()+1-i);
1558       projAll->SetBinContent(projAll->GetNbinsX()+1-k, projAll->GetBinContent(k) / ratioNSAS[i][j][4]);
1559       projAll->SetBinError(projAll->GetNbinsX()+1-k, projAll->GetBinError(k) / ratioNSAS[i][j][4]);
1560       
1561       projAll->SetBinContent(k, 0);
1562       projAll->SetBinError(k, 0);
1563     }
1564     if (!silent)
1565     {
1566       new TCanvas; projAll->DrawCopy();
1567     }
1568     
1569 //     new TCanvas; proj->DrawCopy();
1570     proj->Add(projAll, -1);
1571 //     new TCanvas; proj->DrawCopy();
1572   }
1573   
1574   proj->Scale(hist1->GetYaxis()->GetBinWidth(1));
1575   proj->Scale(1.0 / (etaMax * 2));
1576   proj->Rebin(2); proj->Scale(0.5);
1577
1578   if (gStudySystematic == 20)
1579   {
1580     Printf(">>>>>>>>>>>> Applying non-closure systematics <<<<<<<<<<<<");
1581     file2 = TFile::Open("non_closure.root");
1582     non_closure = (TH1*) file2->Get(Form("non_closure_%d_%d_%d", i, j, 0));
1583     for (Int_t bin=1; bin<=non_closure->GetNbinsX(); bin++)
1584       non_closure->SetBinError(bin, 0);
1585     
1586     proj->Multiply(non_closure);
1587   }
1588   
1589 //   proj = hist1->ProjectionX("proj1x", hist1->GetYaxis()->FindBin(0.5), eta6);
1590
1591   // eta projection
1592   Int_t binEtaProj1 = 1;
1593   Int_t binEtaProj2 = hist1->GetXaxis()->FindBin(-TMath::Pi()/3-0.001);
1594   Int_t binEtaProj3 = binEtaProj2 + 1;
1595   Int_t binEtaProj4 = hist1->GetXaxis()->FindBin(TMath::Pi()/3-0.001);
1596   Int_t binEtaProj5 = binEtaProj4+1;
1597   Int_t binEtaProj6 = hist1->GetXaxis()->FindBin(TMath::Pi()-TMath::Pi()/3-0.001);
1598   Int_t binEtaProj7 = binEtaProj6+1;
1599   Int_t binEtaProj8 = hist1->GetXaxis()->FindBin(TMath::Pi()+TMath::Pi()/3-0.001);
1600   Int_t binEtaProj9 = binEtaProj8+1;
1601   Int_t binEtaProj10 = hist1->GetNbinsX();
1602   Printf("%d %d %d %d %d %d %d %d %d %d", binEtaProj1, binEtaProj2, binEtaProj3, binEtaProj4, binEtaProj5, binEtaProj6, binEtaProj7, binEtaProj8, binEtaProj9, binEtaProj10);
1603   
1604   etaProj = hist1->ProjectionY(Form("%s_proj1y", hist1->GetName()), binEtaProj3, binEtaProj4);
1605   etaProj->Scale(1.0 / (binEtaProj4 - binEtaProj3 + 1));
1606   etaProj2 = hist1->ProjectionY(Form("%s_proj2y", hist1->GetName()), binEtaProj7, binEtaProj8);
1607   etaProj2->Scale(1.0 / (binEtaProj8 - binEtaProj7 + 1));
1608   etaProj3 = hist1->ProjectionY(Form("%s_proj5y", hist1->GetName()), binEtaProj1, binEtaProj2);
1609   etaProj3b = hist1->ProjectionY(Form("%s_proj5by", hist1->GetName()), binEtaProj5, binEtaProj6);
1610   etaProj3c = hist1->ProjectionY(Form("%s_proj5cy", hist1->GetName()), binEtaProj9, binEtaProj10);
1611   etaProj3->Add(etaProj3b);
1612   etaProj3->Add(etaProj3c);
1613   etaProj3->Scale(1.0 / (binEtaProj2 - binEtaProj1 + 1 + binEtaProj6 - binEtaProj5 + 1 + binEtaProj10 - binEtaProj9 + 1));
1614
1615   Float_t ySubBaseline = (etaProj->Integral(yEta1, yEta2) + etaProj->Integral(yEta5, yEta6)) / (yEta2 - yEta1 + 1 + yEta6 - yEta5 + 1);
1616   Float_t ySubPeak = etaProj->Integral(yEta3, yEta4, "width") - ySubBaseline * (yEta4 - yEta3 + 1) * etaProj->GetBinWidth(1);
1617   Printf("Peak (subtracted): %f %f (%.2f%% of unsubtracted peak)", ySubBaseline, ySubPeak, ySubPeak / yPeak * 100);
1618   Printf("    factor[%d][%d][%d] = %.3f;", i, j, centr, ySubPeak / yPeak);
1619   AddPoint(graph[8], xValue1, 100.0 * ySubPeak / yPeak, 0, 0);
1620   
1621   hist1->Rebin2D(4, 4); hist1->Scale(1.0 / 16);
1622   hist1->GetYaxis()->SetRangeUser(-1.99, 1.99);
1623   
1624   TString fitOption = "0I";
1625   if (!silent)
1626   {
1627     c = new TCanvas("c", "c", 600, 600);
1628
1629     PadFor2DCorr();
1630     
1631     hist1->GetYaxis()->SetRangeUser(-1.99, 1.99);
1632     hist1->GetXaxis()->SetTitleOffset(1.7);
1633     hist1->GetYaxis()->SetTitleOffset(2);
1634     hist1->GetZaxis()->SetTitle(kCorrFuncTitle);
1635     hist1->GetZaxis()->SetTitleOffset(2);
1636     hist1->GetZaxis()->SetNdivisions(504);
1637     hist1->SetStats(kFALSE);
1638     hist1->GetZaxis()->CenterTitle(kTRUE);
1639     hist1->GetYaxis()->CenterTitle(kTRUE);
1640     hist1->GetXaxis()->CenterTitle(kTRUE);
1641     hist1->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1642     hist1->GetYaxis()->SetNdivisions(505);
1643     
1644     TString label(hist1->GetTitle());
1645     hist1->SetTitle("");
1646     label.ReplaceAll(".00", "");
1647     label.ReplaceAll(".0", "");
1648     TObjArray* objArray = label.Tokenize("-");
1649     TPaveText* paveText = new TPaveText(0.03, 0.86, 0.44, 0.97, "BRNDC");
1650     paveText->SetTextSize(0.035);
1651     paveText->SetFillColor(0);
1652     paveText->SetShadowColor(0);
1653     paveText->SetBorderSize(0);
1654     paveText->SetFillStyle(0);
1655     TString tmpStr(objArray->At(0)->GetName());
1656     tmpStr.ReplaceAll("p_", "#it{p}_");
1657     paveText->AddText(tmpStr + "GeV/#it{c}");
1658
1659     TString tmpStr(objArray->At(1)->GetName());
1660     tmpStr.ReplaceAll("p_", "#it{p}_");
1661     paveText->AddText(tmpStr + "GeV/#it{c}");
1662   
1663     TString label2(hist2->GetTitle());
1664     TObjArray* objArray2 = label2.Tokenize("-");
1665
1666     TPaveText* paveText2 = new TPaveText(0.65, 0.86, 0.98, 0.97, "BRNDC");
1667     paveText2->SetTextSize(0.035);
1668     paveText2->SetBorderSize(0);
1669     paveText2->SetFillColor(0);
1670     paveText2->SetFillStyle(0);
1671     paveText2->SetShadowColor(0);
1672     paveText2->AddText("p-Pb #sqrt{s_{NN}} = 5.02 TeV");
1673     if (objArray->GetEntries() > 3 && objArray2->GetEntries() > 3)
1674     {
1675       TString centrBeginStr;
1676       centrBeginStr = objArray->At(2)->GetName();
1677       centrBeginStr.ReplaceAll(" ", "");
1678       TString centrEndStr;
1679       centrEndStr = objArray2->At(2)->GetName();
1680       centrEndStr.ReplaceAll(" ", "");
1681       paveText2->AddText(Form("(%s-%s) - (%s-%s)", centrBeginStr.Data(), objArray->At(3)->GetName(), centrEndStr.Data(), objArray2->At(3)->GetName()));
1682     }
1683
1684     hist1->DrawCopy("SURF1");
1685     paveText->Draw();
1686     paveText2->Draw();
1687     gPad->GetCanvas()->SaveAs(Form("ridge_%d_%d.png", i, j));
1688     gPad->GetCanvas()->SaveAs(Form("ridge_%d_%d.eps", i, j));
1689     gPad->GetCanvas()->SaveAs("fig3a.eps");
1690     
1691     Float_t fontSize = 0.05;
1692
1693     c3 = new TCanvas("c3", "c3", 600, 400);
1694     gPad->SetLeftMargin(0.12);
1695     gPad->SetRightMargin(0.01);
1696     gPad->SetBottomMargin(0.12);
1697     gPad->SetTopMargin(0.01);
1698     etaProj->SetTitle();
1699 //     etaProj->GetYaxis()->SetNdivisions(505);
1700     etaProj->GetYaxis()->SetLabelSize(fontSize);
1701     etaProj->GetXaxis()->SetLabelSize(fontSize);
1702     etaProj->GetXaxis()->SetTitleSize(fontSize);
1703     etaProj->GetYaxis()->SetTitleSize(fontSize);
1704     etaProj->GetYaxis()->SetTitle(kProjYieldTitleEta);
1705     etaProj->GetYaxis()->SetTitleOffset(1.1);
1706     etaProj->SetStats(0);
1707     etaProj->GetXaxis()->SetNdivisions(505);
1708     etaProj->Rebin(2); etaProj->Scale(0.5);
1709     etaProj2->Rebin(2); etaProj2->Scale(0.5);
1710     etaProj3->Rebin(2); etaProj3->Scale(0.5);
1711     etaProj->Rebin(2); etaProj->Scale(0.5);
1712     etaProj2->Rebin(2); etaProj2->Scale(0.5);
1713     etaProj3->Rebin(2); etaProj3->Scale(0.5);
1714     etaProj->GetXaxis()->SetRangeUser(-1.99, 1.99);
1715 //     etaProj->GetYaxis()->SetRangeUser(TMath::Min(etaProj->GetMinimum(), etaProj3->GetMinimum()) * 0.8, TMath::Max(etaProj->GetMaximum(), etaProj3->GetMaximum()) * 1.2);
1716     etaProj->GetYaxis()->SetRangeUser(proj->GetMinimum() * 0.98, proj->GetMaximum() * 1.065);
1717     etaProj2->SetLineColor(2); etaProj2->SetMarkerColor(2);
1718     etaProj3->SetLineColor(4); etaProj3->SetMarkerColor(4);
1719     etaProj->SetMarkerStyle(24);
1720     etaProj2->SetMarkerStyle(25);
1721     etaProj3->SetMarkerStyle(26);
1722     etaProj->DrawCopy("E0 X0");
1723     etaProj2->DrawCopy("E0 X0 SAME");
1724     etaProj3->DrawCopy("E0 X0 SAME");
1725     
1726     legend5 = new TLegend(0.48, 0.74, 0.92, 0.97);
1727     legend5->SetBorderSize(0);
1728     legend5->SetTextSize(fontSize);
1729     legend5->SetFillColor(0);
1730     legend5->AddEntry(etaProj,  "|#Delta#varphi| < #pi/3", "P");
1731     legend5->AddEntry(etaProj2, "|#Delta#varphi - #pi| < #pi/3", "P");
1732 //     legend5->AddEntry(etaProj3, "#pi/3 < |#Delta#varphi| < 2#pi/3, #Delta#varphi > 4#pi/3", "P");
1733     legend5->AddEntry(etaProj3, "Remaining #Delta#varphi", "P");
1734     legend5->Draw();
1735     
1736     c = new TCanvas("c2", "c2", 600, 400);
1737     gPad->SetLeftMargin(0.12);
1738     gPad->SetRightMargin(0.01);
1739     gPad->SetBottomMargin(0.12);
1740     gPad->SetTopMargin(0.01);
1741     
1742     proj->SetStats(0);
1743 //     proj->GetYaxis()->SetNdivisions(505);
1744     proj->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1745     proj->GetYaxis()->SetTitle(kProjYieldTitlePhi);
1746     proj->GetYaxis()->SetTitleOffset(1.1);
1747     proj->GetYaxis()->SetLabelSize(fontSize);
1748     proj->GetXaxis()->SetLabelSize(fontSize);
1749     proj->GetXaxis()->SetTitleSize(fontSize);
1750     proj->GetYaxis()->SetTitleSize(fontSize);
1751     proj->SetTitle();
1752     proj->SetMarkerStyle(21);
1753     proj->SetMarkerSize(0.7);
1754     proj->Draw("E0 X0");
1755     proj->SetMinimum(proj->GetMinimum() * 0.98); proj->SetMaximum(proj->GetMaximum() * 1.065);
1756 //     proj->SetMinimum(proj->GetMinimum() * 1.2); proj->SetMaximum(proj->GetMaximum() * 0.6);
1757     fitOption = "I";
1758     
1759     if (0)
1760     {
1761       Printf("\nPer-trigger yield per unit of delta eta (y) as function of delta phi (x) in the bin %s", centralityStr.Data());
1762       Printf("Systematic uncertainties are mostly correlated and affect the baseline. Uncorrelated uncertainties are less than 1%% and not indicated explicitly in the table below.");
1763
1764       for (Int_t k=1; k<=proj->GetNbinsX(); k++)
1765         Printf("x = %.2f, y = %.4f +- %.4f (stat)", proj->GetXaxis()->GetBinCenter(k), proj->GetBinContent(k), proj->GetBinError(k));
1766     }
1767   }
1768   
1769   fileProj = TFile::Open("dphi_proj.root", "UPDATE");
1770   proj->Write(Form("proj_%d_%d_%d_%d", i, j, centr, gStudySystematic));
1771   fileProj->Close();
1772
1773   TF1* v2 = new TF1("func", "[0]+2*[1]*cos(2*x)", -5, 5);
1774   v2->SetLineStyle(2);
1775 //   v2->SetLineWidth(1);
1776   
1777   TF1* v2v3 = new TF1("func", "[0]+2*[1]*cos(2*x)+2*[2]*cos(3*x)", -5, 5);
1778   v2v3->SetLineColor(2);
1779 //   v2v3->FixParameter(2, 0);
1780   proj->Fit(v2, fitOption);
1781 //   return;
1782     
1783   fitOption += "+";
1784   proj->Fit(v2v3, fitOption, "E0 X0");
1785
1786   Float_t min = v2v3->GetMinimum();
1787   Float_t diffMinParam0 = v2v3->GetParameter(0) - min;
1788   Printf("Chi2: %f ndf: %d chi2/ndf %f:Min: %f %f", v2v3->GetChisquare(), v2v3->GetNDF(), v2v3->GetChisquare() / v2v3->GetNDF(), min, diffMinParam0);
1789   AddPoint(graph[10], xValue1, v2v3->GetChisquare() / v2v3->GetNDF(), 0, 0);
1790
1791   AddPoint(graph[11], xValue1, baseLine + diffMinParam0, 0, baseLineE);
1792 //   AddPoint(graph[11], xValue1, baseLine, 0, baseLineE);
1793   
1794   Float_t v2value = 0, v2E = 0, v3 = 0, v3E = 0;
1795   if (v2v3->GetParameter(1) > 0)
1796   {
1797     v2value = TMath::Sqrt(v2v3->GetParameter(1) / (baseLine + diffMinParam0));
1798     v2E = 0.5 * v2value * TMath::Sqrt(v2v3->GetParError(1) * v2v3->GetParError(1) / v2v3->GetParameter(1) / v2v3->GetParameter(1) + baseLineE * baseLineE / baseLine / baseLine);
1799     if (symmetricpT)
1800       AddPoint(graph[4], xValue1vn, v2value, 0, v2E);
1801   }
1802   if (v2v3->GetParameter(2) > 0)
1803   {
1804     v3 = TMath::Sqrt(v2v3->GetParameter(2) / (baseLine + diffMinParam0));
1805     v3E = 0.5 * v3 * TMath::Sqrt(v2v3->GetParError(2) * v2v3->GetParError(2) / v2v3->GetParameter(2) / v2v3->GetParameter(2) + baseLineE * baseLineE / baseLine / baseLine);
1806     if (symmetricpT)
1807       AddPoint(graph[5], xValue2vn, v3, 0, v3E);
1808   }
1809   if (v2v3->GetParameter(1) > 0 && v2v3->GetParameter(2) > 0 && symmetricpT)
1810   {
1811     Float_t v3v2Ratio = TMath::Sqrt(v2v3->GetParameter(2) / v2v3->GetParameter(1));
1812     Float_t v3v2RatioE = 0.5 * v3v2Ratio * TMath::Sqrt(v2v3->GetParError(1) * v2v3->GetParError(1) / v2v3->GetParameter(1) / v2v3->GetParameter(1) + v2v3->GetParError(2) * v2v3->GetParError(2) / v2v3->GetParameter(2) / v2v3->GetParameter(2));
1813     AddPoint(graph[7], xValue1vn, v3v2Ratio, 0, v3v2RatioE);
1814   }
1815   
1816   Printf("Baseline: %f +- %f; v2 = %f +- %f; v3 = %f +- %f", baseLine, baseLineE, v2value, v2E, v3, v3E);
1817   
1818   if (gStudySystematic == 30)
1819   {
1820     // alternative way for baseline
1821     
1822     parabola = new TF1("parabola", "[0] + [1]*(x - [2])**2",  0, TMath::Pi());
1823     parabola->SetLineColor(3);
1824     parabola->SetParameters(min, -0.1, TMath::Pi() / 2);
1825 //     parabola->SetParLimits(1, 0, 1);
1826     proj->Fit(parabola, fitOption, "", TMath::Pi() / 2 - 1, TMath::Pi() / 2 + 1);
1827 //     proj->Fit(parabola, fitOption, "", 0.1, 2);
1828     
1829     min = parabola->GetParameter(0);
1830     Printf("New baseline: %f", min);
1831   }
1832
1833   fitOption += "R";
1834   
1835   if (0)
1836   {
1837     // single gauss fit
1838     
1839     TF1* gaus1 = new TF1("gaus1", "[0]+gaus(1)", -0.5 * TMath::Pi(), 0.5 * TMath::Pi());
1840     gaus1->SetParameters(min, 1, 0, 1);
1841   //   gaus1->FixParameter(0, min);
1842     gaus1->FixParameter(2, 0);
1843     gaus1->SetParLimits(1, 0.001, 10);
1844     gaus1->SetParLimits(3, 0.1, 2);
1845     gaus1->SetLineColor(3);
1846     proj->Fit(gaus1, fitOption);
1847     
1848     TF1* gaus2 = new TF1("gaus2", "[0]+gaus(1)", 0.5 * TMath::Pi(), 1.5 * TMath::Pi());
1849     gaus2->SetParameters(min, 1, TMath::Pi(), 1);
1850     gaus2->FixParameter(2, TMath::Pi());
1851     gaus2->SetParLimits(1, 0.001, 10);
1852     gaus2->SetParLimits(3, 0.1, 2);
1853     gaus2->SetLineColor(3);
1854     proj->Fit(gaus2, fitOption);
1855   
1856     AddPoint(graph[2], xValue1, gaus1->GetParameter(3), 0, gaus1->GetParError(3));
1857     AddPoint(graph[3], xValue2, gaus2->GetParameter(3), 0, gaus2->GetParError(3));
1858   }
1859   else if (0)
1860   {
1861     // combined gauss fit
1862     
1863     TF1* gausBoth = new TF1("gaus2", "[0]+[1]*(exp(-0.5*(x/[2])**2)+exp(-0.5*((x-TMath::TwoPi())/[2])**2))+[3]*(exp(-0.5*((x-TMath::Pi())/[4])**2)+exp(-0.5*((x+TMath::Pi())/[4])**2))", -0.5 * TMath::Pi(), 1.5 * TMath::Pi());
1864     gausBoth->SetParameters(min, 1, 1, 1, 1);
1865   //   gausBoth->FixParameter(0, min);
1866     gausBoth->SetParLimits(1, 0.0001, 10);
1867     gausBoth->SetParLimits(3, 0.0001, 10);
1868     gausBoth->SetParLimits(2, 0.1, 4);
1869     gausBoth->SetParLimits(4, 0.1, 4);
1870     gausBoth->SetLineColor(6);
1871     
1872     proj->Fit(gausBoth, fitOption);
1873   
1874     AddPoint(graph[2], xValue1, gausBoth->GetParameter(2), 0, gausBoth->GetParError(3));
1875     AddPoint(graph[3], xValue2, gausBoth->GetParameter(4), 0, gausBoth->GetParError(4));
1876   }
1877   else
1878   {
1879     // RMS without baseline
1880     projRMS = (TH1*) proj->Clone();
1881     projRMS->Add(new TF1("f", "1", -10, 10), -min);
1882     projRMS->GetXaxis()->SetRangeUser(-TMath::Pi()/2 + 0.01, TMath::Pi()/2 - 0.01);
1883     AddPoint(graph[2], xValue1, projRMS->GetRMS(), 0, projRMS->GetRMSError());
1884     projRMS->GetXaxis()->SetRangeUser(TMath::Pi()/2 + 0.01, 3 * TMath::Pi()/2 - 0.01);
1885     AddPoint(graph[3], xValue2, projRMS->GetRMS(), 0, projRMS->GetRMSError());
1886 //     new TCanvas; projRMS->Draw(); return;
1887   }
1888   
1889 //   new TCanvas; gausBoth->Draw(); return;
1890
1891   if (!silent)
1892   {
1893     TF1* v2v3_v2 = new TF1("func", "[0]+2*[1]*cos(2*x)", -5, 5);
1894     v2v3_v2->SetParameters(v2v3->GetParameter(0), v2v3->GetParameter(1));
1895     v2v3_v2->SetLineStyle(3);
1896     v2v3_v2->SetLineColor(2);
1897 //     v2v3_v2->Draw("SAME");
1898     TF1* v2v3_v3 = new TF1("func", "[0]+2*[1]*cos(3*x)", -5, 5);
1899     v2v3_v3->SetParameters(v2v3->GetParameter(0), v2v3->GetParameter(2));
1900     v2v3_v3->SetLineStyle(4);
1901     v2v3_v3->SetLineColor(2);
1902 //     v2v3_v3->Draw("SAME");
1903     
1904     line = new TLine(-0.5 * TMath::Pi(), min, 1.5 * TMath::Pi(), min);
1905     line->SetLineWidth(2);
1906     line->SetLineColor(4);
1907     line->Draw();
1908
1909     legend = new TLegend(0.47, 0.65, 0.88, 0.95);
1910     legend->SetFillColor(0);
1911     legend->SetBorderSize(0);
1912
1913     legend->AddEntry(proj, "Data", "P");
1914     legend->AddEntry(v2v3, "a_{0} + a_{2} cos(2#Delta#varphi) + a_{3} cos(3#Delta#varphi)", "L");
1915     legend->AddEntry(v2, "a_{0} + a_{2} cos(2#Delta#varphi)", "L");
1916     legend->AddEntry(line, "Baseline for yield extraction", "L");
1917
1918     if (1)
1919     {
1920       hijingFile = TFile::Open("dphi_proj_hijing.root");
1921       if (hijingFile)
1922       {
1923         hijing = (TH1*) hijingFile->Get(Form("proj_%d_%d_%d_0", i, j, centr));
1924         if (hijing)
1925         {
1926           hijing->Add(new TF1("flat", "1", -5, 5), min - ((i == 1) ? 0.192078 : 0.373710));
1927           hijing->SetMarkerColor(2);
1928           hijing->SetLineColor(2);
1929           hijing->SetMarkerStyle(25);
1930           hijing->Draw("SAME E0X0");
1931           legend->AddEntry(hijing, "HIJING shifted", "P");
1932         }
1933       }
1934     }
1935     
1936     legend->SetTextSize(fontSize);
1937     legend->Draw();
1938     
1939 //     paveText3 = (TPaveText*) paveText2->Clone();
1940 //     paveText3->SetX1NDC(0.16);
1941 //     paveText3->SetY1NDC(0.77);
1942 //     paveText3->SetX2NDC(0.42);
1943 //     paveText3->SetY2NDC(0.94);
1944 //     paveText3->Draw();
1945     paveText4 = (TPaveText*) paveText2->Clone();
1946     paveText4->SetTextSize(fontSize);
1947     paveText4->SetX1NDC(0.16);
1948     paveText4->SetY1NDC(0.68);
1949     paveText4->SetX2NDC(0.42);
1950     paveText4->SetY2NDC(0.96);
1951
1952     TString tmpStr(objArray->At(0)->GetName());
1953     tmpStr.ReplaceAll("p_", "#it{p}_");
1954     paveText4->AddText(tmpStr + "GeV/#it{c}");
1955
1956     TString tmpStr(objArray->At(1)->GetName());
1957     tmpStr.ReplaceAll("p_", "#it{p}_");
1958     paveText4->AddText(tmpStr + "GeV/#it{c}");
1959     
1960     paveText4->Draw();
1961     
1962 //     DrawLatex(0.27, 0.19, 1, Form("%sGeV/#it{c}       %sGeV/#it{c}", objArray->At(0)->GetName(), objArray->At(1)->GetName()), fontSize);
1963     
1964     gPad->GetCanvas()->SaveAs(Form("ridge_fit_%d_%d.png", i, j));
1965     gPad->GetCanvas()->SaveAs(Form("ridge_fit_%d_%d.eps", i, j));
1966     gPad->GetCanvas()->SaveAs("fig3b.eps");
1967     
1968     c3->cd();
1969     paveText3 = (TPaveText*) paveText4->Clone();
1970     paveText3->Draw();
1971     gPad->GetCanvas()->SaveAs(Form("ridge_eta_%d_%d.png", i, j));
1972     gPad->GetCanvas()->SaveAs(Form("ridge_eta_%d_%d.eps", i, j));
1973     gPad->GetCanvas()->SaveAs("fig3c.eps");
1974   }
1975   
1976   if (gStudySystematic != 50 && gStudySystematic != 51)
1977   {
1978     Int_t phi1 = 1;
1979     Int_t phi2 = proj->FindBin(TMath::Pi()/2-0.001);
1980     Int_t phi3 = phi2 + 1;
1981     Int_t phi4 = proj->GetNbinsX();
1982   }
1983   else
1984   {
1985     Printf(">>>> Using |dphi| < 1.2 <<<<");
1986     Int_t phi1 = proj->FindBin(-1.2);;
1987     Int_t phi2 = proj->FindBin(1.2);
1988     Int_t phi3 = proj->FindBin(TMath::Pi() - 1.2);
1989     Int_t phi4 = proj->FindBin(TMath::Pi() + 1.2);
1990     Printf("%d %d %d %d", phi1, phi2, phi3, phi4);
1991   }
1992   
1993   Double_t nsYield, asYield, nsYieldE, asYieldE;
1994   
1995   nsYield = proj->IntegralAndError(phi1, phi2, nsYieldE, "width");
1996   nsYield -= min * proj->GetBinWidth(1) * (phi2 - phi1 + 1);
1997   
1998   asYield = proj->IntegralAndError(phi3, phi4, asYieldE, "width");
1999   asYield -= min * proj->GetBinWidth(1) * (phi4 - phi3 + 1);
2000   
2001   AddPoint(graph[0], xValue1, nsYield, 0, nsYieldE);
2002   AddPoint(graph[1], xValue2, asYield, 0, asYieldE);
2003   
2004   if (nsYieldE > 0 && asYieldE > 0)
2005     AddPoint(graph[6], xValue1, nsYield / asYield, 0, nsYield / asYield * TMath::Sqrt(nsYieldE * nsYieldE / nsYield / nsYield + asYieldE * asYieldE / asYield / asYield));
2006   
2007   AddPoint(graph[9], xValue1, 100.0 * ySubPeak / (etaMax * 2) / nsYield, 0, 0);
2008   
2009   Printf("Yields: %f +- %f ; %f +- %f", nsYield, nsYieldE, asYield, asYieldE);
2010 }
2011
2012 Int_t colors[] = { 1, 2, 3, 4, kGray+1, 6, 7, 8 };
2013
2014 void DrawSeveral(Int_t n, const char** graphFiles, Int_t id)
2015 {
2016   ReadGraphs(graphFiles[0]);
2017   TGraphErrors*** base = graphs;
2018
2019   Float_t yMax[] = { 0.1, 0.1, 2, 2, 0.25, 0.25, 4, 1.5, 50, 50, 4, 4 };
2020   Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34 };
2021
2022   TString baseName(graphFiles[0]);
2023   baseName.ReplaceAll(".", "_");
2024   TCanvas* canvas = new TCanvas(Form("%s_%d", baseName.Data(), id), Form("%s_%d", baseName.Data(), id), 800, 600);
2025 //   Printf("%p", canvas);
2026   gPad->SetGridx();
2027   gPad->SetGridy();
2028   dummy = new TH2F(Form("hist_%s_%d", graphFiles[0], id), Form(";%s;%s", graphs[0][id]->GetXaxis()->GetTitle(), graphs[0][id]->GetYaxis()->GetTitle()), 100, 0, 100, 100, 0, yMax[id]);
2029   dummy->SetStats(0);
2030   dummy->DrawCopy();
2031   
2032   TCanvas* canvas2 = new TCanvas(Form("%s_%d_ratio", baseName.Data(), id), Form("%s_%d", baseName.Data(), id), 800, 600);
2033   gPad->SetGridx();
2034   gPad->SetGridy();
2035   dummy = new TH2F(Form("hist_%s_%d_ratio", graphFiles[0], id), Form(";%s;%s ratio", graphs[0][id]->GetXaxis()->GetTitle(), graphs[0][id]->GetYaxis()->GetTitle()), 100, 0, 100, 100, 0, 2);
2036   dummy->SetStats(0);
2037   dummy->DrawCopy();
2038
2039   legend = new TLegend(0.4, 0.6, 0.99, 0.99);
2040   legend->SetFillColor(0);
2041   
2042   for (Int_t fc = 0; fc<n; fc++)
2043   {
2044     ReadGraphs(graphFiles[fc]);
2045     
2046     for (Int_t i=0; i<NGraphs; i++)
2047     {
2048       if (TString(graphFiles[0]).Contains("cms") && i != 2 && i != 5) continue;
2049     
2050       if (!graphs[i][id])
2051         continue;
2052       
2053       graphs[i][id]->Print();
2054     
2055       canvas->cd();
2056 //       Printf("%p", canvas);
2057       graphs[i][id]->SetMarkerStyle(markers[i%16]);
2058       graphs[i][id]->SetMarkerColor(colors[fc]);
2059       graphs[i][id]->SetLineColor(colors[fc]);
2060       GraphShiftX((TGraphErrors*) graphs[i][id]->DrawClone("PSAME"), 0.5 / n * fc);
2061
2062       if (fc == 0 && graphs[i][id]->GetN() > 0)
2063         legend->AddEntry(graphs[i][id], graphs[i][id]->GetTitle(), "P");
2064       
2065       if (fc > 0)
2066       {
2067         canvas2->cd();
2068         DivideGraphs(graphs[i][id], base[i][id]);
2069         GraphShiftX((TGraphErrors*) graphs[i][id]->DrawClone("PSAME"), 0.5 / n * fc);
2070       }
2071     }
2072   }
2073
2074   canvas->cd();
2075   legend->Draw();
2076   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2077   canvas->SaveAs(Form("%s.png", canvas->GetName()));
2078
2079   canvas2->cd();
2080   legend->Draw();
2081   canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
2082   canvas2->SaveAs(Form("%s.png", canvas2->GetName()));
2083 }
2084
2085 void CMSRidge()
2086 {
2087   CreateGraphStructure();
2088   
2089   // scaling done for zyam level difference (note that the values in fig2 have to be multiplied by 2 due to the |dphi| instead of dphi
2090   
2091   // Fig3a
2092   // 0.334516912727  0.0104728088208
2093   // 0.721360493366  0.0263651653896
2094   // 1.1900579331  0.0301999626238
2095   // 1.68965302436  0.024186548724
2096   // 2.19286114745  0.0123441101352
2097   // 2.69021366723  0.00994871155963
2098   // 3.38790257273  0.00434394401877
2099
2100 //   AddPoint(graphs[0][0], 55.3, 0.0263651653896, 0, 0);
2101   AddPoint(graphs[2][0], 55.383968, 0.0301999626238 * 0.72 * 2 / 1.336, 0, 0);
2102   AddPoint(graphs[5][0], 55.383968, 0.0123441101352 * 0.154 * 2 / 0.229, 0, 0);
2103   
2104   // Fig3b
2105   // 128.344370861  0.0254622516556
2106   // 98.1456953642  0.0154940397351
2107   // 56.6225165563  0.00493774834437
2108   // 17.6158940397  0.00037880794702
2109   // Ref multiplicity for centrality 0.000000 to 3.000000: 55.284322
2110   // Ref multiplicity for centrality 3.000000 to 10.000000: 41.124550
2111   // Ref multiplicity for centrality 10.000000 to 50.000000: 24.334303
2112   // Ref multiplicity for centrality 50.000000 to 100.000000: 8.008294
2113
2114 //   AddPoint(graphs[2][0], 55.3, 0.0254622516556, 0, 0);
2115   AddPoint(graphs[2][0], 40.649288, 0.0154940397351 * 0.5 * 2 / 0.985, 0, 0);
2116   AddPoint(graphs[2][0], 23.783224, 0.00493774834437 * 0.29 * 2 / 0.506, 0, 0);
2117
2118   WriteGraphs("graphs_cms.root");
2119 }
2120
2121 void ALICECMSMethod()
2122 {
2123   CreateGraphStructure();
2124   
2125   AddPoint(graphs[2][0], 10, 0.027, 0, 0);
2126   AddPoint(graphs[4][0], 10, 0.115 / 1.34 * 1, 0, 0);
2127   AddPoint(graphs[5][0], 10, 0.019 / 0.152 * 0.126, 0, 0);
2128   
2129   AddPoint(graphs[2][4], 10, 0.09, 0, 0);
2130 //   AddPoint(graphs[4][4], 10, 0.115, 0, 0);
2131   AddPoint(graphs[5][4], 10, 0.137, 0, 0);
2132
2133   AddPoint(graphs[2][5], 10, 0.04, 0, 0);
2134 //   AddPoint(graphs[4][5], 10, 0.046, 0, 0);
2135   AddPoint(graphs[5][5], 10, 0.06, 0, 0);
2136
2137   WriteGraphs("graphs_alice_cmsmethod.root");
2138 }
2139
2140 void GetSystematic()
2141 {
2142   fileProj = TFile::Open("dphi_proj.root", "RECREATE");
2143   fileProj->Close();
2144
2145   gROOT->SetBatch(kTRUE);
2146   
2147   if (1)
2148   {
2149     const char* baseFile = "dphi_corr_pA_121119_3.root";
2150 //     const char* baseFile = "dphi_corr_pA_121108_hybrid_corrected.root";
2151     TString fileTag = "graphs_121119";
2152   }
2153   else if (0)
2154   {
2155     const char* baseFile = "dphi_corr_pA_121029_zna.root";
2156     TString fileTag = "graphs_121113_zna_uncorrected";
2157   }
2158   else if (0)
2159   {
2160     // needs dphi setting...
2161     const char* baseFile = "dphi_corr_pA_121119_cms_2.root";
2162     TString fileTag = "graphs_cms_121119";
2163   }
2164   else if (1)
2165   {
2166     const char* baseFile = "dphi_corr_pA_121121_cmsmethod.root";
2167     TString fileTag = "graphs_cmmethod_121121";
2168   }
2169   
2170   CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + ".root");
2171   
2172   gStudySystematic = 10;
2173   CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_exclusion05.root");
2174   
2175   gStudySystematic = 11;
2176   CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_exclusion00.root");
2177
2178   gStudySystematic = 12;
2179   CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_exclusionAS.root");
2180
2181   gStudySystematic = 13;
2182   CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_exclusionScale.root");
2183
2184   gStudySystematic = 14;
2185   CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_exclusion12.root");
2186   
2187 //   return;
2188   
2189   gStudySystematic = 20;
2190   CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_nonclosure.root");
2191
2192   gStudySystematic = 30;
2193   CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_baseline.root");
2194   
2195   gStudySystematic = 40;
2196   CorrelationSubtractionHistogram("dphi_corr_pA_121116_global.root", kTRUE, fileTag + "_trackcuts.root");
2197
2198   gStudySystematic = 60;
2199   CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_otherperipheral.root");
2200 }
2201
2202 void DrawSystematics(TString fileTag = "graphs_121119", const char* suffix = "")
2203 {
2204   if (0)
2205   {
2206     const Int_t n = 6;
2207     const char* filesBase[] = { "", "_exclusion05", "_exclusion12", "_exclusion00", "_exclusionAS", "_exclusionScale" };
2208   }
2209   else if (0)
2210   {
2211     const Int_t n = 5;
2212     const char* filesBase[] = { "", "_trackcuts", "_nonclosure", "_baseline", "_otherperipheral" };
2213   }
2214   else
2215   {
2216     const Int_t n = 3;
2217     const char* filesBase[] = { "", "_analytical", "_asgap" };
2218   }
2219   
2220   const char* files[n];
2221   for (Int_t i=0; i<n; i++)
2222   {
2223     str = new TString;
2224     str->Form("%s%s%s.root", fileTag.Data(), filesBase[i], suffix);
2225     files[i] = str->Data();
2226   }
2227   
2228 //   const Int_t plots = 6;
2229 //   Int_t ids[7] = { 0, 1, 2, 3, 4, 5, 7 };
2230   
2231   const Int_t plots = 1;
2232   Int_t ids[] = { 4 };
2233
2234   for (Int_t i=0; i<plots; i++)
2235   {
2236     DrawSeveral(n, files, ids[i]);
2237 //     break;
2238   }
2239   
2240   new TCanvas;
2241   for (Int_t i=0; i<n; i++)
2242     DrawLatex(0.2, 0.9 - 0.1 * i, colors[i], files[i]);
2243 }
2244
2245 void DrawSeveral(const char* graphFile1, const char* graphFile2, const char* graphFile3, Int_t id)
2246 {
2247   const char* files[3] = { graphFile1, graphFile2, graphFile3 };
2248   Int_t n = 1;
2249   if (graphFile2)
2250     n++;
2251   if (graphFile3)
2252     n++;
2253   
2254   DrawSeveral(n, files, id);
2255 }
2256   
2257 void DrawYield(const char* graphFile)
2258 {
2259   DrawGraph(graphFile, 0, 1, "Ridge yield per #Delta#eta", "fig4b.eps", kTRUE);
2260 }
2261
2262 void DrawRMS(const char* graphFile)
2263 {
2264   DrawGraph(graphFile, 2, 3, "#sigma", 0);
2265 }
2266
2267 void Drawv2v3(const char* graphFile)
2268 {
2269   DrawGraph(graphFile, 4, 5, "#it{v}_{2} , #it{v}_{3}", "fig4a.eps");
2270 }
2271
2272 void AddSystUnc(TGraphErrors* graph, Float_t syst020, Float_t syst2060)
2273 {
2274   for (Int_t j=0; j<graph->GetN(); j++)
2275   {
2276     Float_t syst = syst2060;
2277     if (graph->GetX()[j] < 20)
2278       syst = syst020;
2279   
2280 //     Printf("%d %f", j, syst);
2281     graph->GetEY()[j] = TMath::Sqrt(graph->GetEY()[j] * graph->GetEY()[j] + syst * syst * graph->GetY()[j] * graph->GetY()[j]);
2282   }
2283 }
2284
2285 void SetSystUnc(TGraphErrors* graph, Float_t syst020, Float_t syst2060)
2286 {
2287   for (Int_t j=0; j<graph->GetN(); j++)
2288   {
2289     Float_t syst = syst2060;
2290     if (graph->GetX()[j] < 20)
2291       syst = syst020;
2292   
2293 //     Printf("%d %f", j, syst);
2294     graph->GetEY()[j] = syst * graph->GetY()[j];
2295   }
2296 }
2297
2298 void SetXError(TGraphErrors* graph, Float_t value)
2299 {
2300   for (Int_t j=0; j<graph->GetN(); j++)
2301     graph->GetEX()[j] = value;
2302 }
2303
2304 void SetYError(TGraphErrors* graph, Float_t value)
2305 {
2306   for (Int_t j=0; j<graph->GetN(); j++)
2307     graph->GetEY()[j] = value;
2308 }
2309
2310 void DrawGraph(const char* graphFile, Int_t id1, Int_t id2, const char* yLabel = 0, const char* outputFileName, Bool_t corrGraph = kFALSE)
2311 {
2312   ReadGraphs(graphFile);
2313   TGraphErrors*** graphsSyst = graphs;
2314   ReadGraphs(graphFile);
2315   TGraphErrors*** graphsCombinedError = graphs;
2316   ReadGraphs(graphFile);
2317   
2318   if (yLabel == 0)
2319     yLabel = graphs[0][id1]->GetYaxis()->GetTitle();
2320   
2321   Float_t yMax[] = { 0.12, 0.12, 2, 2, 0.2, 0.25, 4, 1.5, 50, 50, 4 };
2322   Int_t markers[] = { 20, 21, 22, 23, 29, 33, 34, 2, 3, 5, 2, 3, 5, 2, 3, 5 };
2323   Int_t markers2[] = { 24, 25, 26, 32, 30, 27, 28, 2, 3, 5, 2, 3, 5, 2, 3, 5 };
2324   Float_t markerSize[] = { 1.7, 1.7, 1.7, 1.7, 2.0, 2.0, 1.7, 1.7, 1.7, 1.7, 1.7, 1.7, 1.7, 1.7, 1.7 };
2325
2326   if (1)
2327   {
2328     // default systs (for all pT, centr)
2329     Float_t syst020Array[] = { 0.16, 0.18, 0.15, 0.15, 0.14, 0.23 };
2330     Float_t syst2060Array[] = { 0.23, 0.25, 0.23, 0.23, 0.18, 0.42 };
2331     
2332     // 0.5<1 0.5<1
2333     Float_t syst020ArrayGr0[] = { 0.23, 0.25, 0.15, 0.15, 0.14, 0.23 };
2334     Float_t syst2060ArrayGr0[] = { 0.42, 0.43, 0.23, 0.23, 0.18, 0.42 };
2335
2336     for (Int_t i=0; i<NGraphs; i++)
2337     {
2338       Float_t syst020 = syst020Array[id1];
2339       Float_t syst2060 = syst2060Array[id1];
2340       
2341       if (i == 0)
2342       {
2343         syst020 = syst020ArrayGr0[id1];
2344         syst2060 = syst2060ArrayGr0[id1];
2345       }
2346       
2347 //       Printf("%d %d %f %f", i, id1, syst020, syst2060);
2348       SetSystUnc(graphsSyst[i][id1], syst020, syst2060);
2349       AddSystUnc(graphsCombinedError[i][id1], syst020, syst2060);
2350       SetXError(graphsSyst[i][id1], 0.3);
2351
2352       Float_t syst020 = syst020Array[id2];
2353       Float_t syst2060 = syst2060Array[id2];
2354       
2355       if (i == 0)
2356       {
2357         syst020 = syst020ArrayGr0[id2];
2358         syst2060 = syst2060ArrayGr0[id2];
2359       }
2360
2361       SetSystUnc(graphsSyst[i][id2], syst020, syst2060);
2362       AddSystUnc(graphsCombinedError[i][id2], syst020, syst2060);
2363       SetXError(graphsSyst[i][id2], 0.3);
2364     }
2365     
2366 //     graphs[0][id1]->Print();
2367   }
2368   
2369   const char* eventClass[] = { "0-20%", "20-40%", "40-60%" };
2370   Printf("\n%s:", graphTitles[id1]);
2371   for (Int_t i=0; i<3; i++)
2372   {
2373     Printf("Event class: %s minus 60-100%%", eventClass[i]);
2374     
2375     for (Int_t j=0; j<NGraphs; j++)
2376       if (graphs[j][id1]->GetN() > i)
2377         Printf("%s: %.4f +- %.4f (stat) +- %.4f (syst)", graphs[j][id1]->GetTitle(), graphs[j][id1]->GetY()[i], graphs[j][id1]->GetEY()[i], graphsSyst[j][id1]->GetEY()[i]);
2378   }
2379     
2380   Printf("\n%s:", graphTitles[id2]);
2381   for (Int_t i=0; i<3; i++)
2382   {
2383     Printf("Event class: %s minus 60-100%%", eventClass[i]);
2384     
2385     for (Int_t j=0; j<NGraphs; j++)
2386       if (graphs[j][id2]->GetN() > i)
2387         Printf("%s: %.4f +- %.4f (stat) +- %.4f (syst)", graphs[j][id2]->GetTitle(), graphs[j][id2]->GetY()[i], graphs[j][id2]->GetEY()[i], graphsSyst[j][id2]->GetEY()[i]);
2388   }
2389
2390   TCanvas* canvas = new TCanvas;
2391   gPad->SetTopMargin(0.03);
2392   gPad->SetRightMargin(0.01);
2393   gPad->SetLeftMargin(0.14);
2394   gPad->SetBottomMargin(0.13);
2395 //   gPad->SetGridx();
2396 //   gPad->SetGridy();
2397   TH2F* dummy = new TH2F("dummy", Form(";%s;%s", graphs[0][id1]->GetXaxis()->GetTitle(), yLabel), 3, 0, 60, 100, 0, yMax[id1]);
2398   dummy->GetYaxis()->SetTitleOffset(1.1);
2399   dummy->SetStats(0);
2400   dummy->GetYaxis()->SetNdivisions(505);
2401   dummy->GetXaxis()->SetLabelSize(0.06);
2402   dummy->GetYaxis()->SetLabelSize(0.06);
2403   dummy->GetXaxis()->SetTitleSize(0.06);
2404   dummy->GetYaxis()->SetTitleSize(0.06);
2405   dummy->Draw();
2406   
2407   if (strcmp(graphs[0][id1]->GetXaxis()->GetTitle(), "Event class") == 0)
2408   {
2409     dummy->GetXaxis()->SetBinLabel(1, "0-20%");
2410     dummy->GetXaxis()->SetBinLabel(2, "20-40%");
2411     dummy->GetXaxis()->SetBinLabel(3, "40-60%");
2412     dummy->GetXaxis()->SetLabelSize(0.09);
2413     dummy->GetXaxis()->SetTitleOffset(1);
2414   }
2415   
2416   legend = new TLegend((id1 == 4) ? 0.47 : 0.33, (id1 == 4) ? 0.70 : 0.55, 0.95, 0.92);
2417   legend->SetNColumns(2);
2418   if (id1 == 0 || id1 == 2)
2419     legend->SetHeader("Near side   Away side");
2420   else if (id1 == 4)
2421     legend->SetHeader("  #it{v}_{2}       #it{v}_{3}");
2422     
2423   legend->SetFillColor(0);
2424   legend->SetBorderSize(0);
2425
2426   if (0)
2427   {
2428     Int_t fillStyle[11] = { 3001, 3002, 3003, 3004, 3005, 3006, 3007, 3008, 3009, 3010, 3011 };
2429     
2430     for (Int_t i=0; i<NGraphs; i++)
2431     {
2432       graphsSyst[i][id1]->SetMarkerStyle(1);
2433       graphsSyst[i][id1]->SetMarkerColor(1);
2434       graphsSyst[i][id1]->SetLineColor(1);
2435       graphsSyst[i][id1]->SetFillColor(1);
2436       graphsSyst[i][id1]->SetFillStyle(3001);
2437       graphsSyst[i][id1]->Draw("2SAME");
2438
2439       graphsSyst[i][id2]->SetMarkerStyle(1);
2440       graphsSyst[i][id2]->SetMarkerColor(2);
2441       graphsSyst[i][id2]->SetLineColor(2);
2442       graphsSyst[i][id2]->SetFillColor(2);
2443       graphsSyst[i][id2]->SetFillStyle(3001);
2444       graphsSyst[i][id2]->Draw("2SAME");
2445     }
2446   }
2447   else
2448     Printf(">>>>>>>>>>>> SKIPPING SYST");
2449
2450   
2451   for (Int_t i=0; i<NGraphs; i++)
2452   {
2453     if (graphs[i][id1]->GetN() == 0)
2454       continue;
2455     
2456     graphs[i][id1]->SetMarkerStyle(markers[i]);
2457     graphs[i][id1]->SetMarkerSize(markerSize[i]);
2458     graphs[i][id1]->SetMarkerColor(1);
2459     graphs[i][id1]->SetLineColor(1);
2460     graphs[i][id1]->Draw("PSAME");
2461
2462     graphs[i][id2]->SetMarkerStyle(markers2[i]);
2463     graphs[i][id2]->SetMarkerSize(markerSize[i]);
2464     graphs[i][id2]->SetMarkerColor(2);
2465     graphs[i][id2]->SetLineColor(2);
2466     graphs[i][id2]->Draw("PSAME");
2467
2468     TString label(graphs[i][id1]->GetTitle());
2469     label.ReplaceAll(".00", ".0");
2470     label.ReplaceAll(".50", ".5");
2471     
2472     
2473 /*    label.ReplaceAll(".0", " GeV/#it{c}");
2474     label.ReplaceAll(".5", ".5 GeV/#it{c}");*/
2475     TObjArray* objArray = label.Tokenize("-");
2476     label.Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName());
2477     label.ReplaceAll("-", ";");
2478     
2479     if (id1 == 4)
2480     {
2481       // reduce label
2482       label.ReplaceAll(" ", "");
2483       label.ReplaceAll(";", "<");
2484       tokens = label.Tokenize("<");
2485       label.Form("%s < %s < %s < %s ", tokens->At(0)->GetName(), tokens->At(1)->GetName(), tokens->At(4)->GetName(), tokens->At(5)->GetName());
2486     }
2487     
2488     label.ReplaceAll("p_", "#it{p}_");
2489     label += "GeV/#it{c}";
2490
2491     if (graphs[i][id1]->GetN() > 0)
2492     {
2493       legend->AddEntry(graphs[i][id1], (id1 == 4) ? " " : "    ", "P");
2494       legend->AddEntry(graphs[i][id2], label, "P");
2495     }
2496   }
2497   
2498   legend->Draw();
2499   DrawLatex(0.7, 0.92, 1, "p-Pb #sqrt{s_{NN}} = 5.02 TeV", 0.04);
2500   
2501   canvas->SaveAs(outputFileName);
2502   
2503   if (!corrGraph)
2504     return;
2505   
2506   corr = new TGraphErrors;
2507   for (Int_t i=0; i<NGraphs; i++)
2508   {
2509     for (Int_t j=0; j<graphs[i][id1]->GetN(); j++)
2510       AddPoint(corr, graphsCombinedError[i][id1]->GetY()[j], graphsCombinedError[i][id2]->GetY()[j], graphsCombinedError[i][id1]->GetEY()[j], graphsCombinedError[i][id2]->GetEY()[j]);
2511   }
2512
2513   new TCanvas;
2514   gPad->SetTopMargin(0.03);
2515   gPad->SetRightMargin(0.02);
2516   gPad->SetLeftMargin(0.14);
2517   gPad->SetBottomMargin(0.13);
2518   TString titleString;
2519   titleString.Form(";%s;%s", graphs[0][id1]->GetYaxis()->GetTitle(), graphs[0][id2]->GetYaxis()->GetTitle());
2520   titleString.ReplaceAll("NS", "Near-side");
2521   titleString.ReplaceAll("AS", "Away-side");
2522   titleString.ReplaceAll("Ridge Yield", "ridge yield per #Delta#eta");
2523   dummy = new TH2F("dummy2", titleString, 100, 0, 0.095, 100, 1e-5, 0.095);
2524   dummy->GetYaxis()->SetTitleOffset(1.1);
2525   dummy->SetStats(0);
2526   dummy->GetXaxis()->SetNdivisions(505);
2527   dummy->GetYaxis()->SetNdivisions(505);
2528   dummy->GetXaxis()->SetLabelSize(0.06);
2529   dummy->GetYaxis()->SetLabelSize(0.06);
2530   dummy->GetXaxis()->SetTitleSize(0.06);
2531   dummy->GetYaxis()->SetTitleSize(0.06);
2532   dummy->Draw();
2533
2534   corr->Draw("PSAME");
2535   
2536   line = new TLine(0, 0, 0.095, 0.095);
2537   line->SetLineWidth(2);
2538   line->SetLineStyle(2);
2539   line->Draw();
2540
2541   DrawLatex(0.18, 0.92, 1, "p-Pb #sqrt{s_{NN}} = 5.02 TeV", 0.04);
2542 }
2543
2544 void PaperCorrFuncAll(const char* fileName)
2545 {
2546   Int_t n = 6;
2547   Int_t is[] = { 0, 1, 1, 2, 2, 2, 3 };
2548   Int_t js[] = { 1, 1, 2, 1, 2, 3, 3 };
2549   Int_t centr[] = { 0, 1, 3, 4 };
2550
2551   for (Int_t i=0; i<n; i++)
2552     for (Int_t j=0; j<4; j++)
2553       PaperCorrFunc(fileName, is[i], js[i], centr[j]);
2554 }
2555
2556 void PaperCorrFunc(const char* fileName, Int_t i = 2, Int_t j = 2, Int_t centr = 0)
2557 {
2558   TFile::Open(fileName);
2559   
2560   TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr));
2561
2562   if (!hist1)
2563     return 0;
2564   // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
2565   hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
2566
2567   hist1 = (TH2*) hist1->Clone();
2568   hist1->Rebin2D(2, 2); hist1->Scale(0.25);
2569 //   hist1->Rebin2D(2, 1); hist1->Scale(0.5);
2570   
2571   hist1->GetYaxis()->SetRangeUser(-1.99, 1.99);
2572   hist1->GetXaxis()->SetTitleOffset(1.5);
2573   hist1->GetYaxis()->SetTitleOffset(2);
2574   hist1->GetZaxis()->SetTitle(kCorrFuncTitle);
2575   hist1->SetStats(kFALSE);
2576   hist1->GetZaxis()->SetTitleOffset(2);
2577   hist1->GetZaxis()->SetNdivisions(504);
2578   hist1->SetStats(kFALSE);
2579   hist1->GetZaxis()->CenterTitle(kTRUE);
2580   hist1->GetYaxis()->CenterTitle(kTRUE);
2581   hist1->GetXaxis()->CenterTitle(kTRUE);
2582   hist1->GetYaxis()->SetNdivisions(505);
2583   hist1->GetXaxis()->SetTitle("#Delta#varphi (rad)");
2584
2585   c = new TCanvas("c", "c", 600, 600);
2586   gPad->SetPad(0, 0, 1, 1);
2587   gPad->SetLeftMargin(0.2);
2588
2589   TString label(hist1->GetTitle());
2590   hist1->SetTitle("");
2591   label.ReplaceAll(".00", "");
2592   label.ReplaceAll(".0", "");
2593 //   label.ReplaceAll(".00", " GeV/#it{c}");
2594 //   label.ReplaceAll(".0", " GeV/#it{c}");
2595   TObjArray* objArray = label.Tokenize("-");
2596   TPaveText* paveText = new TPaveText(0.03, 0.86, 0.44, 0.97, "BRNDC");
2597   paveText->SetTextSize(0.035);
2598   paveText->SetFillColor(0);
2599   paveText->SetShadowColor(0);
2600   paveText->SetBorderSize(0);
2601   paveText->SetFillStyle(0);
2602   
2603   TString tmpStr(objArray->At(0)->GetName());
2604   tmpStr.ReplaceAll("p_", "#it{p}_");
2605   paveText->AddText(tmpStr + "GeV/#it{c}");
2606
2607   TString tmpStr(objArray->At(1)->GetName());
2608   tmpStr.ReplaceAll("p_", "#it{p}_");
2609   paveText->AddText(tmpStr + "GeV/#it{c}");
2610 //   paveText->AddText(objArray->At(1)->GetName());
2611
2612   TPaveText* paveText2 = new TPaveText(0.63, 0.86, 0.97, 0.97, "BRNDC");
2613   paveText2->SetTextSize(0.035);
2614   paveText2->SetBorderSize(0);
2615   paveText2->SetFillColor(0);
2616   paveText2->SetFillStyle(0);
2617   paveText2->SetShadowColor(0);
2618   if (objArray->GetEntries() == 4)
2619   {
2620     paveText2->AddText(Form("p-Pb #sqrt{s_{NN}} = 5.02 TeV"));
2621     paveText2->AddText(Form("%s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()));
2622   }
2623   else
2624     paveText2->AddText(Form("%s #sqrt{s} = 2.76 TeV", objArray->At(2)->GetName()));
2625
2626   hist1->Draw("SURF1");
2627   paveText->Draw();
2628   paveText2->Draw();
2629   
2630   if (i == 2 && j == 2 && centr == 0)
2631     c->SaveAs("fig1b.eps");
2632   else if (i == 2 && j == 2 && centr == 4)
2633     c->SaveAs("fig1a.eps");
2634   
2635   c->SaveAs(Form("corr_%d_%d_%d.eps", i, j, centr));
2636 }
2637
2638 void ExtractSystematics(const char* fileName = "dphi_proj.root")
2639 {
2640   Int_t i = 2; Int_t j = 2;
2641 //   Int_t i = 0; Int_t j = 1;
2642   Int_t centr = 0;
2643   
2644   Int_t nSyst = 8;
2645   Int_t syst[] = { 0, 10, 11, 12, 13, 20, 30, 40 };
2646   
2647   TFile::Open(fileName);
2648   base = (TH1*) gFile->Get(Form("proj_%d_%d_%d_%d", i, j, centr, syst[0]));
2649   
2650   c = new TCanvas;
2651 //   base->Rebin(2); base->Scale(0.5);
2652   base->Draw();
2653   
2654   Float_t baseValue = base->Integral() / base->GetNbinsX();
2655 //   (base->GetBinContent(1) + base->GetBinContent(base->GetNbinsX()/2+1)) / 2;
2656   
2657   c2 = new TCanvas;
2658   c3 = new TCanvas;
2659   
2660   Int_t color = 2;
2661
2662   for (Int_t n = 1; n<nSyst; n++)
2663   {
2664     hist = (TH1*) gFile->Get(Form("proj_%d_%d_%d_%d", i, j, centr, syst[n]));
2665     if (!hist)
2666       continue;
2667     
2668 //     hist->Rebin(2); hist->Scale(0.5);
2669     c->cd();
2670     hist->SetLineColor(color++);
2671     Float_t baseValue2 = hist->Integral() / hist->GetNbinsX();
2672     Printf("%f %f", baseValue, baseValue2);
2673 //     hist->Add(new TF1("func", "1", -5, 5), baseValue - baseValue2);
2674     hist->Scale(baseValue / baseValue2);
2675     hist->DrawCopy("SAME");
2676     
2677     c2->cd();
2678     hist->DrawCopy((n == 1) ? "HIST" : "HIST SAME")->Divide(base);
2679
2680     c3->cd();
2681     hist->Add(base, -1);
2682     hist->DrawCopy((n == 1) ? "HIST" : "HIST SAME")->GetYaxis()->SetRangeUser(-0.1, 0.1);
2683   }
2684 }
2685
2686 void GetProjectionSystematics(Float_t eta)
2687 {
2688   fileProj = TFile::Open("dphi_proj2.root", "RECREATE");
2689   fileProj->Close();
2690
2691   Int_t i = 2; Int_t j = 2;
2692 //   Int_t i = 0; Int_t j = 1;
2693
2694   Int_t n = 5;
2695   const char* files[] = { "dphi_corr_pA_121119_3.root", "dphi_corr_pA_121116_global.root", "dphi_corr_pA_121119_3.root", "dphi_corr_pA_121116_hybrid.root" /*pair cuts*/, "dphi_corr_pA_121119_hijing.root" };
2696   
2697   for (Int_t centr=0; centr<6; centr++)
2698   {
2699     for (Int_t k=0; k<n; k++)
2700     {
2701       gStudySystematic = 0;
2702       if (k == 2)
2703         gStudySystematic = 20;
2704       
2705       TFile::Open(files[k]);
2706       const char* label = 0;
2707       TH1* hist = GetProjections(i, j, centr, &label, eta);
2708       if (!hist)
2709         continue;
2710       
2711       fileProj = TFile::Open("dphi_proj2.root", "UPDATE");
2712       hist->Write(Form("proj_%d_%d_%d_%d", i, j, centr, k*10));
2713       fileProj->Close();
2714     }
2715   }
2716 }
2717
2718 void CMSPlot()
2719 {
2720   Int_t centrArr[] = { 0, 1, 3, 4 };
2721   Int_t nCentr = 4;
2722   
2723 //   Int_t i = 1; Int_t j = 2;
2724   Int_t i = 2; Int_t j = 2;
2725
2726   if (1)
2727   {
2728     const char* labels[] = { "Data", "HIJING", "DPMJET" };
2729   //   const char* files[] = { "dphi_corr_pA_121119_3.root", "dphi_corr_121121_hijing_uncorrected.root", "dphi_corr_121121_dpmjet_uncorrected.root" };
2730 //     const char* files[] = { "dphi_corr_pA_121119_3.root", "dphi_corr_121121_hijing_step0.root", "dphi_corr_121121_dpmjet_step0.root" };
2731     const char* files[] = { "dphi_corr_pA_121121_cmsmethod.root", "dphi_corr_121122_cmsmethod_hijing_step0.root", "dphi_corr_121122_cmsmethod_dpmjet_step0.root" };
2732 //     const char* files[] = { "dphi_corr_pA_121122_cnd_cmsmethod.root", "dphi_corr_121122_cmsmethod_hijing_cnd_step0.root", 0 };
2733     Int_t nFiles = 3;
2734   }
2735   else
2736   {
2737     const char* labels[] = { "ALICE method", "CMS method" };
2738 //     const char* files[] = { "dphi_corr_pA_121119_3.root", "dphi_corr_pA_121121_cmsmethod.root" };
2739     const char* files[] = { "dphi_corr_121121_hijing_step0.root", "dphi_corr_121122_cmsmethod_hijing_step0.root" };
2740     Int_t nFiles = 2;
2741   }
2742   
2743   Int_t colors[] = { 1, 2, 4 };
2744   
2745   c = new TCanvas("c", "c", 800, 800);
2746   c->Divide(2, 2);
2747   
2748   Float_t maxValue = -1;
2749   for (Int_t centr = 0; centr < nCentr; centr++)
2750   {
2751     c->cd(centr+1);
2752     
2753     legend = new TLegend(0.15, 0.55, 0.46, 0.85);
2754     legend->SetFillColor(0);
2755     
2756     for (Int_t fileId = 0; fileId < nFiles; fileId++)
2757     {
2758       if (files[fileId] == 0)
2759         continue;
2760       
2761       TFile::Open(files[fileId]);
2762   
2763       TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centrArr[centr]));
2764       if (!hist1)
2765         return 0;
2766       hist1 = (TH2*) hist1->Clone(Form("%s_%.1f", hist1->GetName(), 0.0));
2767
2768       // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
2769       hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
2770       
2771       tokens = TString(hist1->GetTitle()).Tokenize("-");
2772       centralityStr = new TString;
2773       if (tokens->GetEntries() > 2)
2774         *centralityStr = tokens->At(2)->GetName();
2775       if (tokens->GetEntries() > 3)
2776         *centralityStr = *centralityStr + "-" + tokens->At(3)->GetName();
2777       
2778       Float_t etaMin = 1.0;
2779       
2780       proj1x = hist1->ProjectionX(Form("proj1x_%d_%d_%d_%d", i, j, centr, fileId), hist1->GetYaxis()->FindBin(-etaMax+0.01), hist1->GetYaxis()->FindBin(-etaMin-0.01));
2781       proj2x = hist1->ProjectionX(Form("proj2x_%d_%d_%d_%d", i, j, centr, fileId), hist1->GetYaxis()->FindBin(etaMin+0.01), hist1->GetYaxis()->FindBin(etaMax-0.01));
2782       proj1x->Add(proj2x);
2783       
2784       proj1x->Scale(hist1->GetYaxis()->GetBinWidth(1));
2785       proj1x->Scale(1.0 / (etaMax - etaMin) / 2);
2786
2787       proj1x->Rebin(2); proj1x->Scale(0.5);
2788       proj1x->Rebin(2); proj1x->Scale(0.5);
2789       
2790       // zyam
2791       Float_t zyam = proj1x->Integral(proj1x->FindBin(TMath::Pi()/2 - 0.5), proj1x->FindBin(TMath::Pi()/2 + 0.5)) / (proj1x->FindBin(TMath::Pi()/2 + 0.5) - proj1x->FindBin(TMath::Pi()/2 - 0.5) + 1);
2792       proj1x->Add(new TF1("f", "1", -5, 5), -zyam);
2793
2794       proj1x->SetStats(kFALSE);
2795       proj1x->GetXaxis()->SetTitleOffset(1);
2796       
2797       proj1x->SetLineColor(colors[fileId]);
2798       if (maxValue < 0)
2799         maxValue = proj1x->GetMaximum() * 2;
2800       proj1x->SetMaximum(maxValue);
2801 //       proj1x->SetMinimum(-0.01);
2802       proj1x->Draw((fileId == 0) ? "" : "SAME");
2803       
2804       legend->AddEntry(proj1x, labels[fileId]);
2805     }
2806     legend->Draw();
2807   }
2808 }
2809
2810 void FourierFactorization(const char* fileName)
2811 {
2812   ReadGraphs(fileName);
2813   
2814   Int_t n = 6;
2815   Int_t is[] =    { 0, 1, 1, 2, 2, 2, 3 };
2816   Int_t js[] =    { 1, 1, 2, 1, 2, 3, 3 };
2817   Bool_t symm[] = { 1, 0, 1, 0, 0, 1, 0 };
2818
2819   Int_t graphID = 5;
2820   
2821   TGraphErrors** symmGraphs[] = { graphs[0], graphs[2], graphs[5] };
2822   
2823   for (Int_t i=0; i<symmGraphs[0][graphID]->GetN(); i++)
2824   {
2825     Printf("%d", i);
2826     graphs[1][graphID]->GetY()[i] = TMath::Sqrt(symmGraphs[0][graphID]->GetY()[i] * symmGraphs[1][graphID]->GetY()[i]);
2827     graphs[3][graphID]->GetY()[i] = TMath::Sqrt(symmGraphs[2][graphID]->GetY()[i] * symmGraphs[1][graphID]->GetY()[i]);
2828     graphs[4][graphID]->GetY()[i] = TMath::Sqrt(symmGraphs[2][graphID]->GetY()[i] * symmGraphs[2][graphID]->GetY()[i]);
2829   }
2830   
2831 //   graphs[1][graphID]->Draw("A*");
2832   
2833   WriteGraphs();
2834 }
2835
2836 void CompareATLAS()
2837 {
2838   atlas = ReadHepdata("/home/jgrosseo/Dropbox/alice-paper-paridge/comparison_atlas/atlas_figaux10.txt", kFALSE, 1, 1);
2839   atlas->SetMarkerStyle(20);
2840   atlas->Draw("PA");
2841   atlas->SetTitle("");
2842   atlas->GetXaxis()->SetTitle("p_{T} (GeV/c)");
2843   atlas->GetYaxis()->SetTitle("s_{2} / v_{2}");
2844   atlas->GetYaxis()->SetRangeUser(0.02, 0.18);
2845
2846   alice = new TGraphErrors;
2847   AddPoint(alice, 0.75, 0.0583503, 0.25, 0.00831755);
2848   AddPoint(alice, 1.5, 0.0953562, 0.5, 0.0134597);
2849   AddPoint(alice, 3, 0.128309, 1, 0.0189634);
2850   
2851   alice->SetMarkerStyle(21);
2852   alice->SetLineColor(2);
2853   alice->SetMarkerColor(2);
2854   alice->Draw("PSAME");
2855 }
2856
2857 TGraphErrors* GetGraphAsFunctionOfPt(const char* graphFile, Int_t graphID, Int_t centrID, Int_t maxpTIndex)
2858 {
2859   ReadGraphs(graphFile);
2860   
2861   result = new TGraphErrors;
2862   result->GetXaxis()->SetTitle("p_{T} (GeV/c)");
2863   
2864 //   Printf("%s", graphFile);
2865   
2866   for (Int_t i=0; i<maxpTIndex; i++)
2867   {
2868     if (!graphs[i][graphID])
2869       continue;
2870     if (graphs[i][graphID]->GetN() == 0)
2871       continue;
2872     
2873     TObjArray* tokens = TString(graphs[i][graphID]->GetTitle()).Tokenize("-");
2874 //     tokens->Print();
2875     tokens = TString(tokens->At(1)->GetName()).Tokenize("<");
2876     Float_t pt = (TString(tokens->At(0)->GetName()).Atof() + TString(tokens->At(2)->GetName()).Atof()) / 2;
2877     Float_t width = (TString(tokens->At(2)->GetName()).Atof() - TString(tokens->At(0)->GetName()).Atof()) / 2;
2878     
2879     AddPoint(result, pt, graphs[i][graphID]->GetY()[centrID], width, graphs[i][graphID]->GetEY()[centrID]);
2880     
2881     result->GetYaxis()->SetTitle(graphs[i][graphID]->GetYaxis()->GetTitle());
2882   }
2883   
2884   return result;
2885 }
2886
2887 void DrawAsFunctionOfPt(const char* graphFile1, Int_t graphID, Int_t centrID, Int_t maxpTIndex)
2888 {
2889   graph = GetGraphAsFunctionOfPt(graphFile1, graphID, centrID, maxpTIndex);
2890   graph->Draw("AP");
2891 }
2892
2893 TGraphErrors* DrawAsFunctionOfPt(const char* baseFile, Int_t n, const char** graphFiles, Int_t graphID, Int_t centrID, Int_t* maxpTIndex, const char** titles = 0)
2894 {
2895   c = new TCanvas;
2896   if (n > 4)
2897     c2 = new TCanvas;
2898   legend = new TLegend(0.125, 0.71, 0.62, 0.98);
2899   legend->SetFillColor(0);
2900   
2901   TGraphErrors* allGraphs[20];
2902   
2903   for (Int_t i=0; i<n; i++)
2904   {
2905     graph = GetGraphAsFunctionOfPt(TString(baseFile) + graphFiles[i], graphID, centrID, maxpTIndex[i]);
2906     graph->SetLineColor((i%4)+1);
2907     graph->SetMarkerColor((i%4)+1);
2908     graph->SetMarkerStyle(20+i);
2909     RemovePointsBelowX(graph, 0.3);
2910     if (TString(graphFiles[i]).Contains("Kaons"))
2911     {
2912       RemovePointsAboveX(graph, 3.0);
2913       RemovePointsBelowX(graph, 0.3);
2914     }
2915     if (TString(graphFiles[i]).Contains("Protons"))
2916     {
2917       RemovePointsBelowX(graph, 0.5);
2918     }
2919 //     graph->Print();
2920     
2921     c->cd();
2922     if (1 || i < 4)
2923       graph->Draw((i == 0) ? "AP" : "PSAME");
2924     else
2925     {
2926       graph2 = (TGraphErrors*) graph->Clone();
2927       SetXError(graph2, 0);
2928       GraphShiftX(graph2, 0.05);
2929       graph2->Draw("LSAME");
2930     }
2931       
2932     legend->AddEntry(graph, (titles) ? titles[i] : graphFiles[i], "P");
2933     
2934     allGraphs[i] = graph;
2935     
2936     if (i >= 4)
2937     {
2938       graph = (TGraphErrors*) graph->Clone();
2939 //       SetYError(graph, 0);
2940       DivideGraphs(graph, allGraphs[i%4]);
2941       c2->cd();
2942       graph->Draw((i == 4) ? "AP" : "PSAME");
2943       graph->GetYaxis()->SetRangeUser(0.5, 1.5);
2944     }
2945   }
2946   
2947   c->cd();
2948   legend->Draw();
2949   
2950   return allGraphs[0];
2951 }
2952
2953 void DrawAsFunctionOfPt()
2954 {
2955   const char* baseFile = "graphs_130513b"; Float_t maxY = 0.25;
2956 //   const char* baseFile = "graphs_130503_nosub"; Float_t maxY = 0.4;
2957   const char* graphFiles[] = { ".root", "_pions.root", "_kaons.root", "_protons.root", 
2958   "_wingremoved.root", "_wingremoved_pions.root", "_wingremoved_kaons.root", "_wingremoved_protons.root",
2959   "_unfolded.root", "_pions_unfolded.root", "_kaons_unfolded.root", "_protons_unfolded.root", 
2960   "_wingremoved_unfolded.root", "_wingremoved_pions_unfolded.root", "_wingremoved_kaons_unfolded.root", "_wingremoved_protons_unfolded.root" 
2961   };
2962    const char* titles[] = { 
2963       "h-h", "h-#pi", "h-K", "h-p", 
2964       "h-h wingremoved", "h-#pi wingremoved", "h-K wingremoved", "h-p wingremoved", 
2965       "h", "#pi", "K", "p", 
2966       "h wingremoved", "#pi wingremoved", "K wingremoved", "p wingremoved" };
2967   Int_t max[] = { NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, 7, 7, 7, 7 };
2968   
2969   Int_t graphID = 4;
2970   Int_t centrID = 1;
2971   
2972   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles, graphID, centrID, max, titles);
2973   first->GetYaxis()->SetRangeUser(0, maxY);
2974   
2975 //   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+4, graphID, centrID, max, titles+4);
2976 //   first->GetYaxis()->SetRangeUser(0, maxY);
2977
2978   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+8, graphID, centrID, max, titles+8);
2979   first->GetYaxis()->SetRangeUser(0, maxY);
2980   
2981   return;
2982
2983   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+12, graphID, centrID, max, titles+12);
2984   first->GetYaxis()->SetRangeUser(0, maxY);
2985
2986   first = DrawAsFunctionOfPt(baseFile, 8, graphFiles, graphID, centrID, max, titles);
2987   first->GetYaxis()->SetRangeUser(0, maxY);
2988
2989   first = DrawAsFunctionOfPt(baseFile, 8, graphFiles+8, graphID, centrID, max, titles+8);
2990   first->GetYaxis()->SetRangeUser(0, maxY);
2991 }
2992
2993 void DrawAsFunctionOfPt2()
2994 {
2995   const char* baseFile = "graphs_130423_nosub";
2996   const char* graphFiles[] = { ".root", "_nomixed.root", "_allpt_unfolded.root", "_nomixed_allpt_unfolded.root", "_nogap.root", "_nomixed_nogap.root",  "_allpt_nogap_unfolded.root", "_nomixed_allpt_nogap_unfolded.root" };
2997   Int_t max[] = { NGraphs, NGraphs, 4, 4, NGraphs, NGraphs, 4, 4 };
2998   
2999   new TCanvas;
3000   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles, 4, 0, max);
3001   first->GetYaxis()->SetRangeUser(0, 0.4);
3002
3003   new TCanvas;
3004   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+4, 4, 0, max+4);
3005   first->GetYaxis()->SetRangeUser(0, 0.4);
3006 }
3007
3008 void DrawAsFunctionOfPt3()
3009 {
3010 //   const char* baseFile = "graphs_130515"; Float_t maxY = 0.25;
3011   const char* baseFile = "graphs_130515_nosub"; Float_t maxY = 0.4;
3012   const char* graphFiles[] = { ".root", "_pions.root", "_kaons.root", "_protons.root", "_unfolded.root", "_pions_unfolded.root", "_kaons_unfolded.root", "_protons_unfolded.root", "_allpt_unfolded.root", "_allpt_pions_unfolded.root", "_allpt_kaons_unfolded.root", "_allpt_protons_unfolded.root" };
3013   const char* titles[] = { "h-h", "h-#pi", "h-K", "h-p", "h", "#pi", "K", "p", "h allpt", "#pi allpt", "K allpt", "p allpt" };
3014   Int_t max[] = { NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, 7, 7, 7, 7 };
3015   
3016   Int_t graphID = 4;
3017   
3018   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles, graphID, 0, max, titles);
3019   first->GetYaxis()->SetRangeUser(0, maxY);
3020   
3021   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+4, graphID, 0, max+4, titles+4);
3022   first->GetYaxis()->SetRangeUser(0, maxY);
3023
3024   first = DrawAsFunctionOfPt(baseFile, 8, graphFiles, graphID, 0, max, titles);
3025   first->GetYaxis()->SetRangeUser(0, maxY);
3026   
3027   return;
3028
3029   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+8, graphID, 0, max+8, titles+8);
3030   first->GetYaxis()->SetRangeUser(0, maxY);
3031
3032   first = DrawAsFunctionOfPt(baseFile, 8, graphFiles+4, graphID, 0, max+4, titles+4);
3033   first->GetYaxis()->SetRangeUser(0, maxY);
3034 }
3035
3036 void DrawAsFunctionOfPt4()
3037 {
3038 //   const char* baseFile = "graphs_130503";
3039   const char* baseFile = "graphs_130503_nosub";
3040   const char* graphFiles[] = { "_allpt_unfolded.root", "_allpt_pions_unfolded.root", "_allpt_kaons_unfolded.root", "_allpt_protons_unfolded.root", "_allpt.root", "_allpt_pions.root", "_allpt_kaons.root", "_allpt_protons.root" };
3041
3042   Int_t max[8] = { 8, 8, 8, 8, 8, 8, 8, 8 };
3043   
3044   Int_t centrID = 0;
3045   
3046   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles, 4, centrID, max);
3047   first->GetYaxis()->SetRangeUser(0, 0.25);
3048   
3049   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+4, 4, centrID, max+4);
3050   first->GetYaxis()->SetRangeUser(0, 0.25);
3051
3052   first = DrawAsFunctionOfPt(baseFile, 8, graphFiles, 4, centrID, max);
3053   first->GetYaxis()->SetRangeUser(0, 0.25);
3054 }
3055
3056 void DrawAsFunctionOfPt5()
3057 {
3058   const char* baseFile = "~/Dropbox/pid-ridge/pPb/graphs_LHC13bc_20130510"; Float_t maxY = 0.4;
3059   const char* graphFiles[] = { 
3060     "_HadronsNoSub.root", "_PionsNoSub.root", "_KaonsNoSub.root", "_ProtonsNoSub.root",
3061     "_HadronsNoSubAnalytical.root", "_PionsNoSubAnalytical.root", "_KaonsNoSubAnalytical.root", "_ProtonsNoSubAnalytical.root"
3062   };
3063   
3064   Int_t max[] = { NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs };
3065   
3066   Int_t graphID = 4;
3067   Int_t centrID = 0;
3068   
3069   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles, graphID, centrID, max, 0);
3070   first->GetYaxis()->SetRangeUser(0, maxY); 
3071
3072   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+4, graphID, centrID, max, 0);
3073   first->GetYaxis()->SetRangeUser(0, maxY); 
3074
3075   first = DrawAsFunctionOfPt(baseFile, 8, graphFiles, graphID, centrID, max, 0);
3076   first->GetYaxis()->SetRangeUser(0, maxY); 
3077 }
3078
3079 void DrawAsFunctionOfPt6()
3080 {
3081   const char* baseFile = "graphs_130515b_nosub"; Float_t maxY = 0.4;
3082   const char* graphFiles[] = { 
3083     ".root", "_pions.root", "_kaons.root", "_protons.root",
3084     "_analytical.root", "_analytical_pions.root", "_analytical_kaons.root", "_analytical_protons.root"
3085   };
3086   
3087   Int_t max[] = { 6, 6, 6, 6, 6, 6, 6, 6 };
3088   
3089   Int_t graphID = 4;
3090   Int_t centrID = 0;
3091   
3092   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles, graphID, centrID, max, 0);
3093   first->GetYaxis()->SetRangeUser(0, maxY); 
3094
3095   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+4, graphID, centrID, max, 0);
3096   first->GetYaxis()->SetRangeUser(0, maxY); 
3097
3098   first = DrawAsFunctionOfPt(baseFile, 8, graphFiles, graphID, centrID, max, 0);
3099   first->GetYaxis()->SetRangeUser(0, maxY); 
3100 }
3101
3102 void CompareCentralityEstimators()
3103 {
3104   const char* baseFile = "~/Dropbox/pid-ridge/pPb/graphs_LHC13bc_"; Float_t maxY = 0.25;
3105   const char* graphFiles[] = { 
3106     "20130510_Hadrons.root", "20130510_Pions.root", "20130510_Kaons.root", "20130510_Protons.root",
3107     "20130517_Hadrons_ZNA.root", "20130517_Pions_ZNA.root", "20130517_Kaons_ZNA.root", "20130517_Protons_ZNA.root"
3108   };
3109   
3110   Int_t max[] = { NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs };
3111   
3112   const char* titles[] = { 
3113       "h V0A", "#pi V0A", "K V0A", "p V0A", 
3114       "h ZNA", "#pi ZNA", "K ZNA", "p ZNA"
3115    };
3116       
3117   Int_t graphID = 4;
3118   Int_t centrID = 0;
3119   
3120   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles, graphID, centrID, max, titles);
3121   first->GetYaxis()->SetRangeUser(0, maxY); 
3122
3123   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+4, graphID, centrID, max, titles+4);
3124   first->GetYaxis()->SetRangeUser(0, maxY); 
3125
3126   first = DrawAsFunctionOfPt(baseFile, 8, graphFiles, graphID, centrID, max, titles);
3127   first->GetYaxis()->SetRangeUser(0, maxY); 
3128 }
3129
3130 void DrawLikeUnlikeSign()
3131 {
3132   const char* baseFile = "~/Dropbox/pid-ridge/pPb/graphs_LHC13bc_20130510"; Float_t maxY = 0.35;
3133   const char* graphFiles[] = { 
3134     "_Hadrons_LikeSign.root", "_Pions_LikeSign.root", "_Kaons_LikeSign.root", "_Protons_LikeSign.root", 
3135     "_Hadrons_UnlikeSign.root", "_Pions_UnlikeSign.root", "_Kaons_UnlikeSign.root", "_Protons_UnlikeSign.root"
3136   };
3137   
3138    const char* titles[] = { 
3139       "h LS", "#pi LS", "K LS", "p LS", 
3140       "h US", "#pi US", "K US", "p US" };
3141   Int_t max[] = { NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs };
3142   
3143   Int_t graphID = 4;
3144   Int_t centrID = 0;
3145   
3146   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles, graphID, centrID, max, titles);
3147   first->GetYaxis()->SetRangeUser(0, maxY);
3148   
3149   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+4, graphID, centrID, max, titles+4);
3150   first->GetYaxis()->SetRangeUser(0, maxY);
3151   
3152   first = DrawAsFunctionOfPt(baseFile, 8, graphFiles, graphID, centrID, max, titles);
3153   first->GetYaxis()->SetRangeUser(0, maxY);
3154 }
3155
3156 void DrawRelativeError(Int_t id = 0)
3157 {
3158   if (id == 0)
3159   {
3160     Int_t centrID = 0;
3161     const char* file1 = "~/Dropbox/pid-ridge/pPb/graphs_LHC13bc_20130510_Kaons.root";
3162     const char* file2 = "~/Dropbox/pid-ridge/syst020WithSub.root";
3163     const char* label = "0-20% with subtraction";
3164   }
3165   else if (id == 1)
3166   {
3167     Int_t centrID = 0;
3168     const char* file1 = "~/Dropbox/pid-ridge/pPb/graphs_LHC13bc_20130510_KaonsNoSub.root";
3169     const char* file2 = "~/Dropbox/pid-ridge/syst020WithoutSub.root";
3170     const char* label = "0-20% without subtraction";
3171   }
3172   else if (id == 2)
3173   {
3174     Int_t centrID = 3;
3175     const char* file1 = "~/Dropbox/pid-ridge/pPb/graphs_LHC13bc_20130510_KaonsNoSub.root";
3176     const char* file2 = "~/Dropbox/pid-ridge/syst60100WithoutSub.root";
3177     const char* label = "60-100% without subtraction";
3178   }
3179   
3180   graph = GetGraphAsFunctionOfPt(file1, 4, 0, NGraphs);
3181   RemovePointsAboveX(graph, 3.0);
3182   RemovePointsBelowX(graph, 0.3);
3183   
3184    new TCanvas;
3185   graph->SetMarkerStyle(20);
3186   graph->DrawClone("AP");
3187   
3188   for (Int_t i=0; i<graph->GetN(); i++)
3189   {
3190     if (graph->GetEY()[i] <= 0 || graph->GetY()[i] <= 0)
3191       continue;
3192     
3193     graph->GetY()[i] = graph->GetEY()[i] / graph->GetY()[i];
3194     graph->GetEY()[i] = 0;
3195   }
3196   
3197   legend = new TLegend(0.6, 0.7, 0.9, 0.9);
3198   new TCanvas;
3199   clone = (TGraphErrors*) graph->DrawClone("AP");
3200   clone->GetYaxis()->SetRangeUser(0, 0.5);
3201   legend->AddEntry(clone, "stat", "P");
3202   
3203   TFile::Open(file2);
3204   hist = (TH1*) gFile->Get("hsyst2");
3205   hist->SetLineColor(2);
3206   hist->Draw("SAME");
3207   legend->AddEntry(hist, "syst", "L");
3208   
3209   for (Int_t i=0; i<graph->GetN(); i++)
3210   {
3211     if (graph->GetY()[i] <= 0)
3212       continue;
3213   
3214     Float_t syst = hist->GetBinContent(hist->FindBin(graph->GetX()[i]));
3215     graph->GetY()[i] = TMath::Sqrt(graph->GetY()[i] * graph->GetY()[i] + syst * syst);
3216   }
3217   
3218   graph->SetLineColor(4);
3219   graph->SetMarkerStyle(21);
3220   graph->SetMarkerColor(4);
3221   graph->DrawClone("PSAME");
3222   legend->AddEntry(graph->DrawClone("PSAME"), "comb", "P");
3223   
3224   line = new TLine(0, 0.1, 4.0, 0.1);
3225   line->Draw();
3226   legend->AddEntry(line, "K0-Kch unc", "L");
3227   legend->SetHeader(label);  
3228   
3229   legend->Draw();
3230 }