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