]>
Commit | Line | Data |
---|---|---|
1 | ||
2 | #include "TROOT.h" | |
3 | #include "TF1.h" | |
4 | #include "TH1.h" | |
5 | #include "TH2.h" | |
6 | #include "TNtuple.h" | |
7 | #include "TFile.h" | |
8 | #include "TKey.h" | |
9 | #include "TRandom.h" | |
10 | #include "TMath.h" | |
11 | #include "TGraph.h" | |
12 | #include "TCanvas.h" | |
13 | #include "TLatex.h" | |
14 | #include "TBox.h" | |
15 | #include "TGraphErrors.h" | |
16 | #include "TStyle.h" | |
17 | #include "TPolyLine.h" | |
18 | #include "TLegend.h" | |
19 | #include "TLine.h" | |
20 | #include "TSystem.h" | |
21 | #include "TObjArray.h" | |
22 | #include "TGraph2D.h" | |
23 | #include "TGaxis.h" | |
24 | #include "TArrow.h" | |
25 | #include "TASImage.h" | |
26 | #include "TFitResult.h" | |
27 | #include "TStopwatch.h" | |
28 | #include "TProfile2D.h" | |
29 | #include "TRegexp.h" | |
30 | #include <iostream> | |
31 | #include <fstream> | |
32 | #include <vector> | |
33 | ||
34 | // Plot utilities class - loaded by rootlogon.C + .rootrc | |
35 | #include "Utils.h" | |
36 | PlotUtils utils; | |
37 | const char* inputFileName = ""; // Assign in macro | |
38 | const char* cfDef = "cA"; // s, m, cA or cB. Only for phi projections. | |
39 | double minRidgeDEta = 0.8; | |
40 | double maxRidgeDEta = 1.799; | |
41 | const int gNMax = 12; // Global const -- largest n in VnDelta | |
42 | ||
43 | // Don't do global fit to points outside these pt bin | |
44 | // indices. Globally scoped for availability in vntvna(). Reassigned | |
45 | // within individual functions. | |
46 | double minFitPtBin = 0, maxFitPtBin = 999; | |
47 | ||
48 | // Switches | |
49 | //============================================================ | |
50 | bool usePtBinCenter = 1; | |
51 | bool ispp = 0; | |
52 | ||
53 | // Functions | |
54 | //============================================================ | |
55 | void SetProtonProtonFlag(int val) { ispp = val; } | |
56 | void SetMinRidgeDEta(int val) { minRidgeDEta = val; } | |
57 | bool IsProtonProton() { return ispp; } | |
58 | void SetCFDef(char* def) { cfDef = def; } | |
59 | void SaveGraphs(const char* outputFileName, TString opt = ""); | |
60 | TF1* DoGlobalFit(int n, int k, int includedPtBin1 = 0, int includedPtBin2 = 999, | |
61 | const char* region = "RIDGE", const char* cfDef = "cA", TString opt = ""); | |
62 | void VnDelta(int n, const TH1 &h, double &vnd, double &vnd_err, TString opt=""); | |
63 | void VnDelta(int n, const TF1 &f, double &vnd, double &vnd_err, TString opt=""); | |
64 | void MakeVnDVsQGraphs(int n1, int n2, int k, const char* region, const char* corrtype, TString opt=""); | |
65 | TF1* Harmonic(TH1* h, int n, TString opt); | |
66 | TF1* HarmonicSum(TH1* h, int n1, int n2, TString opt=""); | |
67 | TH1* Hist(const char* region, const char* type, int i, int j, int k, TString opt=""); | |
68 | TH2* EtaPhiHist(int cent, int ptt, int pta); | |
69 | TH2* PttPtaHist(); | |
70 | TH2F* Agreement2DHist(int k, int n); | |
71 | void ijkFromHistName(TH1* h, int& i, int& j, int& k); | |
72 | ||
73 | TGraphErrors* VnDeltaVsN(int i, int j, int k, int nMax, TString opt); | |
74 | TGraphErrors* VnDeltaNFVsN(int i, int j, int k, int nMax, TString opt); | |
75 | TGraphErrors* GlobalvnVsN(int i, int k, int nMax, TString opt); | |
76 | TGraphErrors* VnDeltaVsPtAssoc(int i, int k, int n, TString opt=""); | |
77 | TGraphErrors* VnDeltaVsPtTrig(int j, int k, int n, TString opt=""); | |
78 | TGraphErrors* AgreementVsPtSum(int k, int n); | |
79 | TGraphErrors* VnDVsQ(int n, int k, const char* region = "RIDGE", const char* cfDef = "cA", TString opt=""); | |
80 | TGraphErrors* GlobalFitVsPtAssoc(int i, int k, int n, TString opt=""); // v_n(ptt)v_n(pta) vs. assoc pt | |
81 | TGraphErrors* vnGFvsPt(int n, int k, int ipt1, int ipt2, | |
82 | const char* region = "RIDGE", | |
83 | const char* corrtype = "cA", TString opt = ""); | |
84 | TF1* GlobalFitVsQ(int n, int k, int ipt1=0, int ipt2=999, const char* region="RIDGE", | |
85 | const char* corrtype="cA", TString opt=""); // The real fit function | |
86 | //TF1* GlobalFitVsQSubset(int n, int k, int i, TString opt = ""); | |
87 | TGraph* Luzumv1(int cent, const char* pid); | |
88 | double MomCorrection(int i, int j, int k); | |
89 | ||
90 | // Global fit function - give to TF1 constructor | |
91 | double vntvna(double *x, double *par); | |
92 | double Chi2(TH1* h, TF1* f, TString opt = ""); // "", "ndf" or "prob" | |
93 | double ReducedChi2(int i, int k, int n); | |
94 | double SineError(int i, int j, int k); // RMS of <sin n D phi> vs. n about zero | |
95 | double vnGF(int n, int k, int ptBin, int includedPtBin1 = 0, int includedPtBin2 = 999, | |
96 | const char* region = "RIDGE", const char* cfDef = "cA", TString opt = ""); | |
97 | ||
98 | // Show +/- solutions in chi2 space by varying one fixed parameter (par) | |
99 | TCanvas* FitStudy(int icent, int n, int par); | |
100 | TCanvas* DrawQ(int k, int n, TString opt = ""); | |
101 | TCanvas* SingleDraw(int cent, int ptt, int pta, TString opt=""); | |
102 | TCanvas* SingleDrawEta(int cent, int ptt, int pta, TString opt); | |
103 | TCanvas* SingleDrawPlain(int cent, int ptt, int pta, TString opt); | |
104 | TCanvas* DrawVnVsPtTrig(int k, int npta, int ptabins[], TString opt = ""); | |
105 | TCanvas* SingleDraw2D(int cent, int ptt, int pta, TString opt=""); | |
106 | TCanvas* MultiDraw(int cent, TString opt); | |
107 | TCanvas* DrawVnFromGlobalFit(int n, int ipt1=0, int ipt2=999, int ncb = 999, | |
108 | int centbins[] = 0, TString opt=""); // little v_n | |
109 | TCanvas* TwoPanelVert(double ydiv, const char* canvName, const char* canvTitle); | |
110 | TCanvas* DrawChi2vsPtTrig(int npta, int ptabins[], int ncb, int centbins[], TString opt); | |
111 | TCanvas* DrawAgreement(int icent, TString opt = ""); | |
112 | TCanvas* Drawv1to5(int ncb, int centbins[], int ipt1=0, int ipt2=999, TString opt=""); | |
113 | TCanvas* Drawv2to5(int ncb, int centbins[], int ipt1=0, int ipt2=999, TString opt=""); | |
114 | TCanvas* DrawGlobalvnVsN(int npt, int ptbins[], int ncb, int centbins[], TString opt=""); | |
115 | ||
116 | // functions needed for above routine | |
117 | void initialize(const char* ifName = inputFileName); | |
118 | void AddPrelimStampToCurrentPad(float x1, float y1, float x2, float y2, TString opt=""); | |
119 | ||
120 | void SaveCanvases(TObjArray* canvases, const char* fileName); | |
121 | void SaveCanvasesFromFile(const char* rootFile, const char* targetDir, const char* tag, const char* fileType); | |
122 | TObjArray* GetObjectsFromFile(TFile& file, TString clname, TString dir=""); | |
123 | ||
124 | ||
125 | // Global variables | |
126 | //============================================================ | |
127 | TFile *fin=0; | |
128 | int maxcent; | |
129 | int centlow[100]; // low and high percentile for each selection | |
130 | int centhigh[100]; | |
131 | int maxpt; | |
132 | double ptmean[100]; | |
133 | double ptlow[100]; | |
134 | double pthigh[100]; | |
135 | ||
136 | int colorsPale[] = {kOrange-3, kGreen-6, kAzure-9, kGray+1, kRed-9, kRed-7, kBlue-9, kViolet-9}; | |
137 | int colorsDark[] = {kOrange-8, kGreen+1, kAzure+2, kBlack, kCyan+1, kRed+0, kBlue+1, kViolet+1}; | |
138 | int centColors[] = {kBlack, kRed, kOrange-3, kGreen+1, kCyan+2, kBlue-1, kViolet, kGreen+1, kOrange-3, kRed+0, kOrange-8, kRed+1, kBlue+1}; | |
139 | int centColorsPale[] = {kGray, kRed-9, kOrange-4, kGreen-7, kCyan-7, kBlue-9, kViolet-9, kGreen, kOrange, kRed-4, kOrange-9, kRed-1, kBlue-1}; | |
140 | ||
141 | TLatex ltx; // NDC | |
142 | TLatex ltx2; // absolute | |
143 | ||
144 | // Global container for graphs, TF1s, etc. | |
145 | TObjArray* gList = new TObjArray(); | |
146 | TObjArray* gFitList = new TObjArray(); | |
147 | TF1* axisFn = 0; | |
148 | ||
149 | TGraphErrors *gi, *gj, *gk; | |
150 | char* trigLabel(int i); | |
151 | char* asscLabel(int j); | |
152 | char* centLabel(int k); | |
153 | double MeanPt(int i, int j, int k, TString t_or_a = "t" /* "t" or "a" */, | |
154 | TString opt="" /* "" or "squared" */); | |
155 | int CentBin(double cntLo, double cntHi); | |
156 | int PtBin(double ptLo, double ptHi); | |
157 | int CentColor(double cen1, double cen2); | |
158 | int CentColorPale(double cen1, double cen2); | |
159 | int PtColor(double p1, double p2); | |
160 | int PtColorPale(double p1, double p2); | |
161 | TStopwatch timer; | |
162 | ||
163 | // The business | |
164 | //============================================================ | |
165 | ||
166 | void SaveGraphs(const char* outputFileName, TString opt) | |
167 | { | |
168 | const char* mode = "recreate"; | |
169 | if (opt.Contains("update")) | |
170 | mode = "update"; | |
171 | ||
172 | TFile* outputFile = new TFile(outputFileName, mode); | |
173 | gList->Write(); | |
174 | gFitList->Write(); | |
175 | outputFile->Close(); | |
176 | return; | |
177 | } | |
178 | ||
179 | void initialize(const char* ifName) { | |
180 | ||
181 | // TH1F files with histograms for: | |
182 | // NSJET_y_PTTRIG_PTASSOC_CENT (tracks with |delta-eta| < 0.8) | |
183 | // RIDGE_y_PTTRIG_PTASSOC_CENT (tracks with |delta-eta| = 0.8 - 1.6) | |
184 | ||
185 | if (fin) // don't keep re-initializing | |
186 | return; | |
187 | ||
188 | fin = new TFile(ifName); | |
189 | if (!fin) { | |
190 | Error("FourierPlus - initialize()", "%s not found", ifName); | |
191 | gSystem->Exit(1); | |
192 | } | |
193 | ||
194 | TList* info = (TList*)fin->Get("FileInfo"); | |
195 | if (info) | |
196 | info->Print(); | |
197 | ||
198 | gi = (TGraphErrors*) fin->Get("TrigPtBins"); | |
199 | gj = (TGraphErrors*) fin->Get("AsscPtBins"); | |
200 | gk = (TGraphErrors*) fin->Get("EvCentBins"); | |
201 | ||
202 | maxpt = gi->GetN(); | |
203 | maxcent = gk->GetN(); | |
204 | gList->Add(gi); gList->Add(gj); gList->Add(gk); | |
205 | ||
206 | if (gj->GetN() != maxpt) | |
207 | Warning("initialize()", | |
208 | "Trigger and assc. binning are not symmetric: %d vs. %d", | |
209 | maxpt, gj->GetN() ); | |
210 | ||
211 | for (int k=0; k<maxcent; k++) { | |
212 | centlow[k] = gk->GetX()[k] - gk->GetEX()[k]; | |
213 | centhigh[k] = gk->GetX()[k] + gk->GetEX()[k]; | |
214 | // cout << centlow[k] << " " << centhigh[k] << endl; | |
215 | } | |
216 | ||
217 | for (int i=0; i<maxpt; i++) { | |
218 | ptlow[i] = gi->GetX()[i] - gi->GetEX()[i]; | |
219 | pthigh[i] = gi->GetX()[i] + gi->GetEX()[i]; | |
220 | // cout << ptlow[i] << " " << pthigh[i] << endl; | |
221 | } | |
222 | ||
223 | // OK, leave hard-coded for now.... | |
224 | ptmean[0] = 0.66; // 0.5-1 GeV | |
225 | ptmean[1] = 1.30; // 1-2 GeV | |
226 | ptmean[2] = 2.30; // 2-3 GeV | |
227 | ptmean[3] = 3.30; // 3-4 GeV | |
228 | ptmean[4] = 4.60; // 4-6 GeV | |
229 | ptmean[5] = 6.60; // 6-8 GeV | |
230 | ptmean[6] = 11.7; // 8-15 GeV | |
231 | ||
232 | if (usePtBinCenter) | |
233 | for (int i=0; i<maxpt; i++) { | |
234 | ptmean[i] = gi->GetX()[i]; | |
235 | } | |
236 | ||
237 | gROOT->SetStyle("Plain"); | |
238 | gStyle->SetOptTitle(0); | |
239 | gStyle->SetOptStat(0); | |
240 | ltx.SetNDC(); | |
241 | } | |
242 | ||
243 | TH1* Hist(const char* region, const char* type, int i, int j, int k, TString opt) | |
244 | { | |
245 | // Options: | |
246 | // "sys": return measured pts. w/systematic errors. | |
247 | // "hi_sys": meas. pts + sys, unmodified stat. errors. | |
248 | // "lo_sys": meas. pts - sys, unmodified stat. errors. | |
249 | ||
250 | TH1* h_in=0, *h_out=0; | |
251 | const char* name = Form("PbPb/%s_%s_%d_%d_%d",region,type,i,j,k); | |
252 | h_in = (TH1*)fin->Get(name); | |
253 | if (!h_in) { | |
254 | Error("Hist()","%s not found, exiting.", name); | |
255 | gSystem->Exit(-1); | |
256 | } | |
257 | ||
258 | h_out = (TH1*)h_in->Clone(Form("%s_%s", name, opt.Data())); | |
259 | ||
260 | // Stat --> sys error bars. Error estimated as difference in Fourier | |
261 | // series (0<n<6) between same-event dist. and CF. | |
262 | if (opt.Contains("sys")) { | |
263 | const char* s = Form("PbPb/%s_s_%d_%d_%d",region,i,j,k); | |
264 | TH1* htmp = (TH1*)fin->Get(s); | |
265 | TH1* hs = 0; | |
266 | if (htmp) | |
267 | hs = (TH1*)htmp->Clone(); | |
268 | else | |
269 | Error("Hist()", "No hist %s", s); | |
270 | ||
271 | double c = double(hs->GetNbinsX())/hs->Integral(); | |
272 | hs->Scale(c); | |
273 | ||
274 | TF1* fs = HarmonicSum(hs, 1,5); | |
275 | TF1* fc = HarmonicSum(h_in,1,5); | |
276 | ||
277 | for (int n=1; n<=hs->GetNbinsX(); n++) { | |
278 | double x = hs->GetBinCenter(n); | |
279 | double err = TMath::Abs(fs->Eval(x) - fc->Eval(x)); | |
280 | double bc = h_in->GetBinContent(n); | |
281 | ||
282 | if (opt.Contains("hi_sys")) | |
283 | h_out->SetBinContent(n, bc+err); | |
284 | else if (opt.Contains("lo_sys")) | |
285 | h_out->SetBinContent(n, bc-err); | |
286 | else { | |
287 | h_out->SetBinError(n, err); | |
288 | h_out->SetFillColor(kGray); | |
289 | } | |
290 | } | |
291 | ||
292 | delete htmp; | |
293 | delete hs; | |
294 | delete fs; | |
295 | delete fc; | |
296 | } | |
297 | ||
298 | return h_out; | |
299 | } | |
300 | ||
301 | void ijkFromHistName(TH1* h, int& i, int& j, int& k) | |
302 | { | |
303 | TString tok, name(h->GetName()); | |
304 | Ssiz_t from = 0; | |
305 | int counter = 0; | |
306 | while (name.Tokenize(tok, from, "_")) { | |
307 | if (counter==2) i = tok.Atoi(); | |
308 | if (counter==3) j = tok.Atoi(); | |
309 | if (counter==4) k = tok.Atoi(); | |
310 | counter++; | |
311 | } | |
312 | if (0) | |
313 | cout << i<<j<<k << endl; | |
314 | ||
315 | return; | |
316 | } | |
317 | ||
318 | double PtBinCenterSum(int ptt, int pta) | |
319 | { | |
320 | if ((ptt<0 || ptt>=maxpt) || (pta>ptt)) | |
321 | Error("PtBinCenterSum()","Bad arg(s): ptt %d pta %d", ptt, pta); | |
322 | ||
323 | double pti = gi->GetX()[ptt]; | |
324 | double ptj = gj->GetX()[pta]; | |
325 | return pti+ptj; | |
326 | } | |
327 | ||
328 | int GlobalIndex(int ptt, int pta) | |
329 | { | |
330 | // Find the global index from the trig and assc. pt bins. | |
331 | // It is useful that TMath::Binomial(n,k) returns 0 when n < k. | |
332 | if (pta > ptt) return -1; | |
333 | return TMath::Binomial(ptt+1, 2) + pta; | |
334 | } | |
335 | ||
336 | void ijFromGlobalIndex(int q, int& ptt, int& pta) | |
337 | { | |
338 | // Assign ptt and pta from the global index q. This function is | |
339 | // effectively the inverse of GlobalIndex(i,j). | |
340 | for (int i=0; i<maxpt; i++) { | |
341 | for (int j=0; j<=i; j++) { | |
342 | if (GlobalIndex(i,j)==q) { | |
343 | ptt = i; pta = j; | |
344 | return; | |
345 | } | |
346 | } | |
347 | } | |
348 | return; | |
349 | } | |
350 | ||
351 | // Needs to be updated! Has hard-coded numbers like 27.9. -AMA 11/15/2011 | |
352 | TCanvas* FitStudy(int icent, int n, int par) | |
353 | { | |
354 | // First get the VnDelta vs q graph... | |
355 | TGraphErrors* gVn = VnDVsQ(n, icent); | |
356 | ||
357 | // Then setup fit function | |
358 | TF1 *fn = new TF1(Form("fn_cent%dto%d_n%d", centlow[icent],centhigh[icent], n), | |
359 | vntvna, 0.0, 27.9, maxpt); | |
360 | fn->SetNpx(1000); | |
361 | fn->SetLineColor(kRed); | |
362 | ||
363 | TCanvas* c = new TCanvas(Form("chi2_cent%dto%d_n%d", centlow[icent],centhigh[icent], n), | |
364 | Form("chi2 test for %s, varying v_{%d}(%.1f)", | |
365 | centLabel(icent), n, ptmean[par]), | |
366 | 700, 1200); | |
367 | c->Divide(1, 3, 0.001, 0.001); | |
368 | ||
369 | TGraphErrors *g = new TGraphErrors(); | |
370 | TGraphErrors *chi2 = new TGraphErrors(); | |
371 | chi2->SetTitle(Form("Total #chi^{2} for n=%d;v_{%d}(p_{T}=%.2f);total #chi^{2}", n,n, ptmean[par])); | |
372 | utils.set_tgraph_props(g, kGray, kGray, kFullCircle, 1.3); | |
373 | utils.set_tgraph_props(chi2, kBlue, kBlue, kFullCircle, 0.8); | |
374 | g->SetLineWidth(3); | |
375 | ||
376 | double bestvn = vnGF(n,icent,par,0,999,"RIDGE","cA",""); | |
377 | double bestvn_err = vnGF(n,icent,par,0,999,"RIDGE","cA","err"); | |
378 | ||
379 | int nsteps = 500; | |
380 | ||
381 | double parmax = TMath::Abs(2*(bestvn + bestvn_err)); | |
382 | double parmin = -parmax; | |
383 | ||
384 | // First panel: step thru one fixed vn | |
385 | c->cd(1); | |
386 | TH1F* hf1 = gPad->DrawFrame(0.0, 1.2*parmin, 12.5, 1.2*parmax); | |
387 | hf1->SetTitle(Form("title;p_{T};v_{%d}(p_{T})", n)); | |
388 | ||
389 | for (int step=0; step<=nsteps; step++) { | |
390 | double stepsize = (parmax-parmin)/nsteps; | |
391 | double val = parmin + step*stepsize; | |
392 | ||
393 | for (int ip=0;ip<maxpt;ip++) fn->SetParameter(ip,0.0); | |
394 | fn->FixParameter(par, val); | |
395 | ||
396 | gVn->Fit(fn, "QR"); | |
397 | for (int i=0;i<maxpt;i++) { | |
398 | g->SetPoint(i, ptmean[i], fn->GetParameter(i)); | |
399 | g->SetPointError(i, 0.0, fn->GetParError(i)); | |
400 | } | |
401 | if (step%50==0) | |
402 | g->DrawClone("elpsame"); | |
403 | ||
404 | // store chi2 at this val | |
405 | chi2->SetPoint(step, val, fn->GetChisquare()); | |
406 | } | |
407 | // Finally, draw once with no fixed parameters... | |
408 | utils.set_tgraph_props(g, kBlack, kBlack, kFullCircle, 1.3); | |
409 | for (int ip=0;ip<maxpt;ip++) fn->ReleaseParameter(ip); | |
410 | gVn->Fit(fn,"QR"); | |
411 | for (int i=0;i<maxpt;i++) { | |
412 | g->SetPoint(i, ptmean[i], fn->GetParameter(i)); | |
413 | g->SetPointError(i, 0.0, fn->GetParError(i)); | |
414 | } | |
415 | g->DrawClone("elpsame"); | |
416 | ltx.SetTextSize(0.08); | |
417 | ltx.DrawLatex(0.6, 0.8, Form("v_{%d}, %s", n, centLabel(icent))); | |
418 | ||
419 | // Plot chi2 over full range | |
420 | c->cd(2); | |
421 | c->Update(); | |
422 | chi2->DrawClone("alp"); | |
423 | ||
424 | // and zoomed to see 1 sigma interval | |
425 | c->cd(3); | |
426 | c->Update(); | |
427 | double ymin = 1e9;//chi2->GetHistogram()->GetMinimum(); | |
428 | for (int j=0; j<chi2->GetN(); j++) | |
429 | if (chi2->GetY()[j] < ymin) | |
430 | ymin = chi2->GetY()[j]; | |
431 | double nsigma = n==2? 3 : 3; | |
432 | TH1F* hf3 = gPad->DrawFrame(parmin, ymin - 1, parmax, ymin + nsigma); | |
433 | hf3->SetTitle(Form("Total #chi^{2} for n=%d;v_{%d}(p_{T}=%.2f);total #chi^{2}", n,n, ptmean[par])); | |
434 | chi2->DrawClone("lpsame"); | |
435 | TLine best, onesig; | |
436 | best.SetLineWidth(3); | |
437 | onesig.SetLineWidth(3); | |
438 | onesig.SetLineColor(kGray); | |
439 | best.DrawLine(parmin, ymin, parmax, ymin); | |
440 | onesig.DrawLine(parmin, ymin+1, parmax, ymin+1); | |
441 | TLatex lbl; | |
442 | lbl.DrawLatex(parmin + 0.05*(parmax-parmin), ymin+0.1, Form("Best v_{%d}(%.2f)",n,ptmean[par])); | |
443 | lbl.DrawLatex(parmin + 0.05*(parmax-parmin), ymin+1.1, "+ 1#sigma"); | |
444 | ||
445 | return c; | |
446 | } | |
447 | ||
448 | TCanvas* MultiDraw(int cent, TString opt) | |
449 | { | |
450 | // TODO: make this function much more flexible. Pass in an array of | |
451 | // bins or histos and m,n ints to plot in an m x n grid. | |
452 | TCanvas* c = new TCanvas(Form("cfmulti_%s_%d",opt.Data(),cent), | |
453 | Form("Corr. Functions - %s", centLabel(cent)), | |
454 | 1200, 900); | |
455 | c->Divide(maxpt, maxpt-2, 0, 0); | |
456 | ||
457 | int ipad=1; | |
458 | ||
459 | // Combine first 3 rows (i=0,1,2) to save space | |
460 | // Then draw the rest as usual | |
461 | for (int i=0; i<maxpt; i++) { | |
462 | for (int j=0; j<maxpt; j++) { | |
463 | if(j<=i) { | |
464 | TCanvas* cs = 0; | |
465 | if (opt.Contains("2D")) | |
466 | cs = SingleDraw2D(cent, i, j, ""); | |
467 | else if (opt.Contains("deta")) | |
468 | cs = SingleDrawEta(cent, i, j, ""); | |
469 | else | |
470 | cs = SingleDraw(cent, i, j, ""); | |
471 | c->cd(ipad); | |
472 | cs->DrawClonePad(); | |
473 | ipad++; | |
474 | if (ipad==maxpt) ipad++; | |
475 | } | |
476 | else if (i > 2) ipad++; | |
477 | } | |
478 | } | |
479 | ||
480 | // Add some labeling in empty pads | |
481 | ltx.SetTextSize(0.15); | |
482 | c->cd(1*maxpt); | |
483 | ltx.DrawLatex(0.1, 0.8, "p_{T}^{trig}: 0.5-1,"); | |
484 | ltx.DrawLatex(0.1, 0.6, "1-2, 2-3 GeV/c"); | |
485 | c->cd(2*maxpt); | |
486 | ltx.DrawLatex(0.1, 0.8, "p_{T}^{trig}:"); | |
487 | ltx.DrawLatex(0.1, 0.6, "3-4 GeV/c"); | |
488 | c->cd(3*maxpt); | |
489 | ltx.DrawLatex(0.1, 0.8, "p_{T}^{trig}:"); | |
490 | ltx.DrawLatex(0.1, 0.6, "4-6 GeV/c"); | |
491 | c->cd(4*maxpt); | |
492 | ltx.DrawLatex(0.1, 0.8, "p_{T}^{trig}:"); | |
493 | ltx.DrawLatex(0.1, 0.6, "6-8 GeV/c"); | |
494 | ||
495 | c->cd(2*maxpt-2); | |
496 | ltx.SetTextSize(0.2); | |
497 | ltx.DrawLatex(0.1, 0.7, "Pb+Pb"); | |
498 | ltx.DrawLatex(0.1, 0.4, "LHC10h"); | |
499 | ltx.DrawLatex(0.1, 0.1, centLabel(cent)); | |
500 | ||
501 | return c; | |
502 | } | |
503 | ||
504 | TH2* EtaPhiHist(int cent, int ptt, int pta) | |
505 | { | |
506 | initialize(); | |
507 | TH2* htmp = (TH2*)fin->Get(Form("PbPb/ETAPHI_c_%1d_%1d_%1d",ptt,pta,cent)); | |
508 | if (!htmp) { | |
509 | Error("EtaPhiHist()", | |
510 | "Problem getting pointer to histogram. cent %d ptt %d pta %d", | |
511 | cent, ptt, pta); | |
512 | return 0; | |
513 | } | |
514 | TH2* h = (TH2*)htmp->Clone(Form("EtaPhi_c_%1d_%1d_%1d",ptt,pta,cent)); | |
515 | TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis(), *az = h->GetZaxis(); | |
516 | ax->SetTitle("#Delta#phi"); | |
517 | ay->SetTitle("#Delta#eta"); | |
518 | az->SetTitle("C(#Delta#phi, #Delta#eta)"); | |
519 | ay->SetRangeUser(-maxRidgeDEta, maxRidgeDEta); | |
520 | ax->SetTitleOffset(1.3); | |
521 | ay->SetTitleOffset(1.3); | |
522 | az->SetTitleOffset(1.4); | |
523 | ||
524 | return h; | |
525 | } | |
526 | ||
527 | TCanvas* SingleDraw2D(int cent, int ptt, int pta, TString opt) | |
528 | { | |
529 | initialize(); | |
530 | TCanvas* c = 0; // The canvas to be returned | |
531 | double labelScale = 1.3; | |
532 | double nDiv = 5; | |
533 | ||
534 | // See if we have already made this once | |
535 | const char* name; | |
536 | if (opt=="") | |
537 | name = Form("etaphi_cent%dto%d_%d_%d", | |
538 | centlow[cent],centhigh[cent],ptt,pta ); | |
539 | else | |
540 | name = Form("etaphi_cent%dto%d_%d_%d_%s", | |
541 | centlow[cent],centhigh[cent],ptt,pta, opt.Data() ); | |
542 | ||
543 | const char* title = Form("dphideta_cent%dto%d_pTtrig%.2gto%.2g_pTassoc%.2gto%.2g_%s", | |
544 | centlow[cent],centhigh[cent], | |
545 | ptlow[ptt], pthigh[ptt], | |
546 | ptlow[pta], pthigh[pta], | |
547 | opt.Data()); | |
548 | ||
549 | if (gDirectory->FindObject(name)) { | |
550 | c = (TCanvas*) gDirectory->FindObject(name); | |
551 | return c; | |
552 | } | |
553 | else | |
554 | c = new TCanvas(name, title, 1); | |
555 | ||
556 | TH2* h = EtaPhiHist(cent, ptt,pta); | |
557 | if (!h) { | |
558 | Error("SingleDraw2D()", | |
559 | "Problem getting pointer to histogram. cent %d ptt %d pta %d", | |
560 | cent, ptt, pta); | |
561 | return 0; | |
562 | } | |
563 | utils.make_nice_axes(c, h, labelScale, nDiv, nDiv); | |
564 | TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis(), *az = h->GetZaxis(); | |
565 | ax->SetTitle("#Delta#phi [rad]"); | |
566 | ay->SetTitle("#Delta#eta"); | |
567 | az->SetTitle("C(#Delta#phi, #Delta#eta)"); | |
568 | ay->SetRangeUser(-maxRidgeDEta, maxRidgeDEta); | |
569 | ax->SetTitleOffset(1.3); | |
570 | ay->SetTitleOffset(1.3); | |
571 | az->SetTitleOffset(1.4); | |
572 | ||
573 | double czlo=0, czhi=5.0; | |
574 | if (opt.Contains("zoom")) | |
575 | az->SetRangeUser(czlo, czhi); | |
576 | ||
577 | if (opt.Contains("colz")) | |
578 | h->Draw("colz"); | |
579 | else | |
580 | h->Draw("surf1"); | |
581 | ||
582 | ltx.SetTextSize(0.05); | |
583 | if (opt.Contains("colz")) { | |
584 | if (opt.Contains("zoom")) | |
585 | ltx.DrawLatex(0.6, 0.9, Form("C(#Delta#phi) (C(#Delta#phi) < %.3g)", czhi)); | |
586 | else | |
587 | ltx.DrawLatex(0.6, 0.9, Form("C(#Delta#phi) (full range)")); | |
588 | } | |
589 | ||
590 | // Label pt, centrality, energy | |
591 | TString ptStr(Form("#splitline" | |
592 | "{%.3g < p_{T}^{t} < %.3g GeV/c}" | |
593 | "{%.3g < p_{T}^{a} < %.3g GeV/c}", | |
594 | ptlow[ptt], pthigh[ptt], | |
595 | ptlow[pta], pthigh[pta])); | |
596 | ltx.SetTextSize(0.07); | |
597 | ltx.DrawLatex(0.01, 0.86, ptStr.Data()); | |
598 | if (ispp) | |
599 | ltx.DrawLatex(0.65, 0.9, Form("#splitline{pp 2.76 TeV}{}")); | |
600 | else | |
601 | ltx.DrawLatex(0.65, 0.9, Form("#splitline{Pb-Pb 2.76 TeV}{%s}", centLabel(cent))); | |
602 | ||
603 | return c; | |
604 | } | |
605 | ||
606 | TCanvas* SingleDrawEta(int cent, int ptt, int pta, TString opt) | |
607 | { | |
608 | TCanvas* c = 0; | |
609 | int lineWidth = opt.Contains("small") ? 1 : 4; | |
610 | double mrkSize = opt.Contains("small") ? 0.5 : 1.5; | |
611 | ||
612 | const char* name = Form("deta_cent%dto%d_%d_%d", | |
613 | centlow[cent],centhigh[cent],ptt,pta); | |
614 | ||
615 | // See if we have already made this once | |
616 | if (gDirectory->FindObject(name)) { | |
617 | c = (TCanvas*) gDirectory->FindObject(name); | |
618 | Info("SingleDrawEta()", "%s already created", c->GetName()); | |
619 | return c; | |
620 | } | |
621 | else | |
622 | c = new TCanvas(name, name, 600, 600); | |
623 | ||
624 | TH1 *hNS = 0, *hAS = 0; // divide-then-project | |
625 | TH1 *hNS2 = 0, *hAS2 = 0; // project-then-divide | |
626 | ||
627 | bool includeProjThenDiv = 0; | |
628 | ||
629 | if (1) { // Recreate CFs using divide-then-project method | |
630 | ||
631 | TH2* hc = (TH2*)fin->Get(Form("PbPb/ETAPHI_c_%1d_%1d_%1d",ptt,pta,cent)); | |
632 | if (!hc) { | |
633 | Error("SingleDrawEta()", "Problem getting hc\n"); | |
634 | gSystem->Exit(-1); | |
635 | } | |
636 | ||
637 | int bin1 = hc->GetXaxis()->FindBin(-TMath::PiOver2()+0.001); | |
638 | int bin2 = hc->GetXaxis()->FindBin(TMath::PiOver2()-0.001); | |
639 | int bin3 = hc->GetXaxis()->FindBin(TMath::PiOver2()+0.001); | |
640 | int bin4 = hc->GetXaxis()->FindBin(3*TMath::PiOver2()-0.001); | |
641 | ||
642 | bool smallEtaRange = false; | |
643 | if (smallEtaRange) { | |
644 | bin1 = hc->GetXaxis()->FindBin(-TMath::PiOver2()+0.001); | |
645 | bin2 = hc->GetXaxis()->FindBin(TMath::PiOver2()-0.001); | |
646 | bin3 = hc->GetXaxis()->FindBin(TMath::Pi()-0.5); | |
647 | bin4 = hc->GetXaxis()->FindBin(TMath::Pi()+0.5); | |
648 | } | |
649 | ||
650 | const char* ns_name = Form("ns_deta_cent%dto%d_%d_%d", | |
651 | centlow[cent],centhigh[cent],ptt,pta); | |
652 | const char* as_name = Form("as_deta_cent%dto%d_%d_%d", | |
653 | centlow[cent],centhigh[cent],ptt,pta); | |
654 | ||
655 | hNS = hc->ProjectionY(ns_name, bin1, bin2, "e"); | |
656 | hAS = hc->ProjectionY(as_name, bin3, bin4, "e"); | |
657 | ||
658 | hNS->Scale(1./(bin2-bin1+1)); | |
659 | hAS->Scale(1./(bin4-bin3+1)); | |
660 | ||
661 | } | |
662 | ||
663 | if (includeProjThenDiv) { // As created by Project.C (project-then-divide) | |
664 | const char* ns = Form("PbPb/ETA_NS_c_%1d_%1d_%1d",ptt,pta,cent); | |
665 | const char* as = Form("PbPb/ETA_AS_c_%1d_%1d_%1d",ptt,pta,cent); | |
666 | hNS2 = (TH1*)fin->Get(ns); | |
667 | hAS2 = (TH1*)fin->Get(as); | |
668 | if (!hNS2) { | |
669 | Error("SingleDrawEta()", "Problem getting %s\n", ns); | |
670 | gSystem->Exit(-1); | |
671 | } | |
672 | if (!hAS2) { | |
673 | Error("SingleDrawEta()", "Problem getting %s\n", as); | |
674 | gSystem->Exit(-1); | |
675 | } | |
676 | } | |
677 | ||
678 | utils.set_hist_props(hNS, kBlack, kNone, kBlack, kFullDotLarge, mrkSize); | |
679 | utils.set_hist_props(hAS, kRed, kNone, kRed, kFullDotLarge, mrkSize); | |
680 | utils.set_ylimits(hNS, hAS); | |
681 | hNS->SetLineWidth(lineWidth); | |
682 | hAS->SetLineWidth(lineWidth); | |
683 | ||
684 | // hNS->Draw("hist"); | |
685 | // hAS->Draw("histsame"); | |
686 | hNS->DrawClone("ep"); | |
687 | hAS->DrawClone("epsame"); | |
688 | ||
689 | if (includeProjThenDiv) { | |
690 | utils.set_hist_props(hNS2, kGray+2, kNone, kGray+2, kFullDotLarge, mrkSize); | |
691 | utils.set_hist_props(hAS2, kCyan, kNone, kCyan, kFullDotLarge, mrkSize); | |
692 | hNS2->SetLineWidth(lineWidth); | |
693 | hAS2->SetLineWidth(lineWidth); | |
694 | ||
695 | hNS2->DrawClone("epsame"); | |
696 | hAS2->DrawClone("epsame"); | |
697 | } | |
698 | ||
699 | return c; | |
700 | } | |
701 | ||
702 | TCanvas* SingleDrawPlain(int cent, int ptt, int pta, TString opt) | |
703 | { | |
704 | initialize(); | |
705 | const char* region = "RIDGE"; | |
706 | const char* yTitle = Form("%.2g < |#Delta#eta| < %.2g", | |
707 | minRidgeDEta, maxRidgeDEta); | |
708 | if (opt.Contains("ALL")) { | |
709 | region = "ALL"; | |
710 | yTitle = Form("%.2g < |#Delta#eta| < %.2g", 0.0, maxRidgeDEta); | |
711 | } | |
712 | if (opt.Contains("NSJET")) { | |
713 | region = "NSJET"; | |
714 | yTitle = Form("%.2g < |#Delta#eta| < %.2g", 0.0, minRidgeDEta); | |
715 | } | |
716 | ||
717 | const char* name = Form("%s_cent%dto%d_%d_%d_%s", region, | |
718 | centlow[cent],centhigh[cent],ptt,pta, opt.Data()); | |
719 | const char* title = Form("dphiPlain_%s_cent%dto%d_pTtrig%.2gto%.2g_pTassoc%.2gto%.2g_%s", | |
720 | region, | |
721 | centlow[cent],centhigh[cent], | |
722 | ptlow[ptt], pthigh[ptt], | |
723 | ptlow[pta], pthigh[pta], | |
724 | opt.Data()); | |
725 | TCanvas* c = 0; | |
726 | if (gDirectory->FindObject(name)) { | |
727 | ||
728 | c = (TCanvas*) gDirectory->FindObject(name); | |
729 | Info("SingleDrawPlain()", "%s already created", c->GetName()); | |
730 | return c; | |
731 | } | |
732 | else | |
733 | c = new TCanvas(name, title, 700, 700); | |
734 | c->cd(); | |
735 | ||
736 | TH1* h_sys = Hist(region, cfDef, ptt,pta,cent, "sys"); | |
737 | TH1* h = Hist(region, cfDef, ptt,pta,cent); | |
738 | h->SetName("h"); | |
739 | ||
740 | TAxis* ax = h->GetXaxis(); | |
741 | TAxis* ay = h->GetYaxis(); | |
742 | ||
743 | h->SetTitle(Form(";#Delta#phi [rad];C(#Delta#phi), %s", yTitle)); | |
744 | h->SetLineWidth(2); | |
745 | ax->SetNdivisions(208); | |
746 | ay->SetNdivisions(208); | |
747 | ax->SetLabelSize(0.04); | |
748 | ay->SetLabelSize(0.04); | |
749 | ax->SetTitleSize(0.05); | |
750 | ay->SetTitleSize(0.05); | |
751 | ax->SetTitleOffset(1.5); | |
752 | ay->SetTitleOffset(-1.9); | |
753 | c->SetFrameLineWidth(2); | |
754 | c->SetLeftMargin(0.2); | |
755 | c->SetBottomMargin(0.15); | |
756 | c->SetTopMargin(0.01); | |
757 | c->SetRightMargin(0.01); | |
758 | ax->CenterTitle(); | |
759 | ay->CenterTitle(); | |
760 | ax->SetTicks("+-"); | |
761 | ay->SetTicks("+-"); | |
762 | ||
763 | utils.set_ylimits(h, h, 0.05); | |
764 | utils.set_hist_props(h, kBlack, kNone, kBlack, kFullDotLarge, 1.4); | |
765 | h->Draw(); | |
766 | h_sys->Draw("e2psame"); | |
767 | h->Draw("same"); | |
768 | if (opt.Contains("harmonics")) | |
769 | for (int n=1; n<6; n++) | |
770 | Harmonic(h, n, "draw"); | |
771 | if (opt.Contains("sum")) | |
772 | HarmonicSum(h,1,5,"draw"); | |
773 | ||
774 | ltx.SetTextSize(0.045); | |
775 | double top = 0.92, gap = 0.08; | |
776 | double x1 = (TString(region)=="RIDGE" && (ptlow[ptt] > 4 || ispp))? 0.25 : 0.5; | |
777 | ltx.DrawLatex(x1, top-0*gap, Form("p_{T}^{t} %s GeV/c", trigLabel(ptt))); | |
778 | ltx.DrawLatex(x1, top-1*gap, Form("p_{T}^{a} %s GeV/c", asscLabel(pta))); | |
779 | if (ispp) | |
780 | ltx.DrawLatex(x1, top-2*gap, Form("pp 2.76 TeV")); | |
781 | else | |
782 | ltx.DrawLatex(x1, top-2*gap, Form("Pb-Pb 2.76 TeV, %s", centLabel(cent))); | |
783 | if (0) | |
784 | ltx.DrawLatex(0.22, 0.92, Form("|#Delta#eta| > %.1f", minRidgeDEta)); | |
785 | ||
786 | // ltx.SetTextColor(kBlack); | |
787 | // ltx.SetTextSize(0.03); | |
788 | // ltx.DrawLatex(0.7, 0.18, "statistical error only"); | |
789 | ||
790 | return c; | |
791 | } | |
792 | ||
793 | TF1* Harmonic(TH1* h, int n, TString opt) | |
794 | { | |
795 | TF1* f = 0; | |
796 | double vndelta=0, vndelta_err=0; | |
797 | double phiMin = -0.5*TMath::Pi(); | |
798 | double phiMax = +1.5*TMath::Pi(); | |
799 | double b0 = h->Integral() / h->GetNbinsX(); | |
800 | int color = opt.Contains("gray")? kGray : colorsPale[n-1]; | |
801 | ||
802 | VnDelta(n, *h, vndelta, vndelta_err); | |
803 | ||
804 | f = new TF1(Form("%scos%ddphi", h->GetName(), n), | |
805 | "[0]*(1. + 2.*[1]*cos([2]*x))", | |
806 | phiMin, phiMax); | |
807 | f->SetParameters(b0, vndelta, n); | |
808 | f->SetLineColor(color); | |
809 | f->SetLineWidth(3); | |
810 | if (opt.Contains("draw")) | |
811 | f->Draw("same"); | |
812 | ||
813 | return f; | |
814 | } | |
815 | ||
816 | TF1* HarmonicSum(TH1* h, int n1, int n2, TString opt) | |
817 | { | |
818 | if (n1<1) | |
819 | Warning("HarmonicSum()", "n1=%d but should be > 0", n1); | |
820 | ||
821 | double phiMin = -0.5*TMath::Pi(); | |
822 | double phiMax = +1.5*TMath::Pi(); | |
823 | ||
824 | // Function definition string | |
825 | TString fdef("[0]*(1.0 + 2.0*("); | |
826 | for (int n=n1; n<=n2; n++) { | |
827 | fdef.Append(Form(" [%d]*cos(%d*x) ", n, n)); | |
828 | fdef.Append( n<n2 ? "+" : "))"); | |
829 | } | |
830 | ||
831 | // Normalization | |
832 | TF1 *ff = new TF1("ff", fdef.Data(), phiMin,phiMax); | |
833 | ff->SetParameter(0, h->Integral()/h->GetNbinsX()); | |
834 | ||
835 | // Set parameters | |
836 | for (int n=n1; n<=n2; n++) { | |
837 | ||
838 | // Sum of vnt x vna from global fit | |
839 | if (opt.Contains("global") ) { | |
840 | int i=-1, j=-1, k=-1; | |
841 | ijkFromHistName(h,i,j,k); | |
842 | double vnt = vnGF(n,k,i,0,999,"RIDGE","cA",""); | |
843 | double vna = vnGF(n,k,j,0,999,"RIDGE","cA",""); | |
844 | ff->SetParameter(n, vnt*vna); // TODO: consider including errors | |
845 | ff->SetLineColor(kRed+1); | |
846 | } | |
847 | else { // sum of directly-extracted pair moments | |
848 | double Vn=0, VnErr=0; | |
849 | VnDelta(n, *h, Vn, VnErr); | |
850 | ff->SetParameter(n, Vn); | |
851 | ff->SetLineColor(kBlack); // ff->SetLineColor(kAzure); | |
852 | // ff->SetLineWidth(4); | |
853 | ff->SetLineStyle(kDashed); | |
854 | if (n2==10) { | |
855 | ff->SetLineColor(kGreen); | |
856 | } | |
857 | } | |
858 | } | |
859 | ||
860 | if (opt.Contains("draw")) | |
861 | ff->Draw("same"); | |
862 | ||
863 | return ff; | |
864 | } | |
865 | ||
866 | TCanvas* SingleDraw(int cent, int ptt, int pta, TString opt) | |
867 | { | |
868 | int lineWidth = opt.Contains("small") ? 1 : 4; | |
869 | double mrkSize = opt.Contains("small") ? 0.5 : 1.5; | |
870 | TString dEtaLabel = ""; | |
871 | ||
872 | const char* region = "RIDGE"; | |
873 | if (opt.Contains("ALL")) { | |
874 | region = "ALL"; | |
875 | dEtaLabel = Form("%.2g < |#Delta#eta| < %.2g", 0.0, maxRidgeDEta); | |
876 | } | |
877 | if (opt.Contains("NSJET")) { | |
878 | region = "NSJET"; | |
879 | dEtaLabel = Form("%.2g < |#Delta#eta| < %.2g", 0.0, minRidgeDEta); | |
880 | } | |
881 | if (TString(region)=="RIDGE") | |
882 | dEtaLabel = Form("%.2g < |#Delta#eta| < %.2g", minRidgeDEta, maxRidgeDEta); | |
883 | ||
884 | initialize(); | |
885 | ||
886 | TCanvas* c1 = 0; // The canvas to be returned | |
887 | ||
888 | // See if we have already made this once | |
889 | TString name = Form("dphi_etamin%02d_cent%dto%d_%d_%d", (int)(10*minRidgeDEta), | |
890 | centlow[cent],centhigh[cent],ptt,pta); | |
891 | TString title = Form("dphiFourier_%s_cent%dto%d_pTtrig%.2gto%.2g_pTassoc%.2gto%.2g", | |
892 | region, | |
893 | centlow[cent],centhigh[cent], | |
894 | ptlow[ptt], pthigh[ptt], | |
895 | ptlow[pta], pthigh[pta]); | |
896 | ||
897 | if (!opt.IsNull()) { | |
898 | name.Append("_"); name.Append(opt); | |
899 | title.Append("_"); title.Append(opt); | |
900 | } | |
901 | ||
902 | if (gDirectory->FindObject(name.Data())) { | |
903 | ||
904 | c1 = (TCanvas*) gDirectory->FindObject(name.Data()); | |
905 | Info("SingleDraw()", "%s already created", c1->GetName()); | |
906 | return c1; | |
907 | } | |
908 | else | |
909 | c1 = TwoPanelVert(0.3, name.Data(), title.Data()); | |
910 | ||
911 | c1->cd(1); | |
912 | gPad->SetTopMargin(0.1); | |
913 | c1->cd(); | |
914 | ||
915 | TH1* h = Hist(region, cfDef, ptt,pta,cent); | |
916 | TH1* h_sys = Hist(region, cfDef, ptt,pta,cent, "sys"); | |
917 | utils.set_hist_props(h, kBlack, kNone, kBlack, kFullDotLarge, mrkSize); | |
918 | h->SetLineWidth(opt.Contains("small") ? 1 : 2); | |
919 | h->GetXaxis()->SetTitle("#Delta#phi [rad]"); | |
920 | h->GetYaxis()->SetTitle("C(#Delta#phi)"); | |
921 | ||
922 | TF1* ff5 = HarmonicSum(h, 1, 5, ""); | |
923 | ff5->SetLineWidth(lineWidth); | |
924 | TF1* fflowg = HarmonicSum(h, 1, 5, "global"); | |
925 | fflowg->SetLineWidth(lineWidth); | |
926 | TF1* ff10 = HarmonicSum(h, 1, 10, ""); | |
927 | ff10->SetLineWidth(lineWidth); | |
928 | TF1* ff12 = HarmonicSum(h, 1, 12, ""); | |
929 | ff12->SetLineWidth(lineWidth); | |
930 | ||
931 | c1->cd(1); | |
932 | ||
933 | int nXticks = 404, nYticks = 210; | |
934 | double pad1Scale = 1.6; | |
935 | double pad2Scale = 2.3; | |
936 | utils.make_nice_axes(c1, h, pad1Scale, nXticks, nYticks, 0.05, 0.01); | |
937 | h->GetXaxis()->SetTitleOffset(1.2); | |
938 | h->GetYaxis()->SetTitleOffset(1.5); | |
939 | utils.set_ylimits(h,h); | |
940 | ||
941 | h->Draw("ep"); | |
942 | h_sys->Draw("e2psame"); | |
943 | h->Draw("epsame"); | |
944 | ||
945 | if (opt.Contains("gray")) | |
946 | for (int n=1; n<6; n++) | |
947 | Harmonic(h, n, "gray, draw"); | |
948 | else if (opt.Contains("color")) | |
949 | for (int n=1; n<6; n++) | |
950 | Harmonic(h, n, "draw"); | |
951 | ||
952 | if (opt.Contains("upto10")) { | |
953 | ff10->SetLineWidth(lineWidth+2); | |
954 | ff10->SetLineColor(kBlack); | |
955 | ff10->DrawClone("same"); | |
956 | ff10->SetLineWidth(lineWidth); | |
957 | ff10->SetLineColor(kGreen); | |
958 | ff10->SetLineStyle(kDashed); | |
959 | ff10->Draw("same"); | |
960 | } | |
961 | ||
962 | // ff10->Draw("same"); | |
963 | ff5->Draw("same"); | |
964 | h->Draw("p,e,same"); | |
965 | ||
966 | if (opt.Contains("global")) { | |
967 | fflowg->Draw("same"); | |
968 | } | |
969 | ||
970 | ltx.SetTextSize(0.06); | |
971 | double top = 0.82, gap = 0.08; | |
972 | double x1 = (TString(region)=="RIDGE" && ptlow[ptt] > 6)? 0.25 : 0.6; | |
973 | if (ispp) | |
974 | ltx.DrawLatex(0.2, 0.92, Form("pp 2.76 TeV")); | |
975 | else | |
976 | ltx.DrawLatex(0.2, 0.92, Form("Pb-Pb 2.76 TeV, %s central", centLabel(cent))); | |
977 | ltx.DrawLatex(x1, top-0*gap, Form("%.3g < p_{T}^{t} < %.3g GeV/c", | |
978 | ptlow[ptt], pthigh[ptt])); | |
979 | ltx.DrawLatex(x1, top-1*gap, Form("%.3g < p_{T}^{a} < %.3g GeV/c", | |
980 | ptlow[pta], pthigh[pta])); | |
981 | // ltx.DrawLatex(x1, top-0*gap, Form("p_{T}^{trig} %s GeV/c", trigLabel(ptt))); | |
982 | // ltx.DrawLatex(x1, top-1*gap, Form("p_{T}^{assoc} %s GeV/c", asscLabel(pta))); | |
983 | ltx.DrawLatex(x1, top-2*gap, dEtaLabel.Data()); | |
984 | ||
985 | ||
986 | c1->cd(2); | |
987 | ||
988 | TH1F *hsub = static_cast <TH1F *> (h->Clone()); | |
989 | TH1F *hsub10 = static_cast <TH1F *> (h->Clone()); | |
990 | TH1F *hsubg = (TH1F*) (h->Clone()); | |
991 | ||
992 | utils.set_hist_props(hsub10, kBlack, kNone, kBlack, kOpenCircle, mrkSize); | |
993 | utils.set_hist_props(hsub, kAzure, kNone, kAzure, kFullSquare, mrkSize); | |
994 | utils.set_hist_props(hsubg, kRed+1, kNone, kRed+1, kFullSquare, mrkSize); | |
995 | ||
996 | for (int i=1;i<=hsub->GetNbinsX();i++) { | |
997 | double val, err; | |
998 | val = h->GetBinContent(i) * h->GetBinWidth(i); | |
999 | err = h->GetBinError(i) * h->GetBinWidth(i); | |
1000 | double dphi1 = h->GetXaxis()->GetBinLowEdge(i); | |
1001 | double dphi2 = h->GetXaxis()->GetBinUpEdge(i); | |
1002 | double fSum5, fSum10, fSumGlobal; | |
1003 | fSum5 = ff5->Integral(dphi1, dphi2); | |
1004 | fSum10 = ff10->Integral(dphi1, dphi2); | |
1005 | fSumGlobal = fflowg->Integral(dphi1, dphi2); | |
1006 | hsub->SetBinContent(i, val/fSum5); | |
1007 | hsub->SetBinError(i, err/fSum5); | |
1008 | hsub10->SetBinContent(i, val/fSum10); | |
1009 | hsub10->SetBinError(i, err/fSum10); | |
1010 | hsubg->SetBinContent(i, val/fSumGlobal); | |
1011 | hsubg->SetBinError(i, err/fSumGlobal); | |
1012 | } | |
1013 | utils.make_nice_axes(c1, hsub, pad2Scale, nXticks, 5, 0.05, 0.01); | |
1014 | if (opt.Contains("global")) | |
1015 | utils.set_ylimits(hsub, hsubg, 0.2); | |
1016 | else | |
1017 | utils.set_ylimits(hsub, hsub10, 0.2); | |
1018 | ||
1019 | hsub->SetTitle(Form(";#Delta#phi [rad];ratio")); | |
1020 | // hsub->SetTitle(Form(";#Delta#phi [rad];data/#Sigma_{n=1}^{%s}", opt.Contains("global")?"5":"N")); | |
1021 | hsub->GetXaxis()->SetTitleOffset(1.25); | |
1022 | hsub->GetYaxis()->SetTitleOffset(0.7); | |
1023 | hsub->GetXaxis()->SetTicks("+-"); hsub->GetXaxis()->SetTickLength(0.08); | |
1024 | hsub10->SetTitle(";#Delta#phi [rad];data/#Sigma_{n=1}^{N}"); | |
1025 | hsub10->GetXaxis()->SetTitleOffset(1.25); | |
1026 | hsub10->GetYaxis()->SetTitleOffset(0.7); | |
1027 | hsubg->SetTitle(";#Delta#phi [rad];data/calc"); | |
1028 | hsubg->GetXaxis()->SetTitleOffset(1.25); | |
1029 | hsubg->GetYaxis()->SetTitleOffset(0.7); | |
1030 | ||
1031 | TLegend* leg = new TLegend(0.2, 0.01, 0.38, 0.38); | |
1032 | leg->SetFillColor(kNone); | |
1033 | leg->SetBorderSize(1); | |
1034 | if (opt.Contains("global")) { | |
1035 | leg->AddEntry(hsub, "#LTcos(n#Delta#phi)#GT", "epl"); | |
1036 | leg->SetBorderSize(0); | |
1037 | } | |
1038 | else | |
1039 | leg->AddEntry(hsub, "1#leqn#leq5", "epl"); | |
1040 | ||
1041 | TLine unity; | |
1042 | hsub->Draw("ep"); | |
1043 | ||
1044 | if (opt.Contains("global")) { | |
1045 | hsubg->Draw("epsame"); | |
1046 | leg->AddEntry(hsubg, "v(p_{Tt})v(p_{Ta})", "epl"); | |
1047 | } | |
1048 | else if (opt.Contains("upto10")) { | |
1049 | TH1* hsub10green = (TH1*)hsub10->Clone(); | |
1050 | utils.set_hist_props(hsub10green, kBlack, kNone, kGreen, kFullCircle, mrkSize); | |
1051 | hsub10green->Draw("epsame"); | |
1052 | hsub10->Draw("epsame"); | |
1053 | leg->AddEntry(hsub10green, "1#leqn#leq10", "epl"); | |
1054 | } | |
1055 | ||
1056 | unity.DrawLine(-TMath::PiOver2(), 1.0, 3*TMath::PiOver2(), 1.0); | |
1057 | ||
1058 | c1->cd(2); | |
1059 | if (opt.Contains("global")) | |
1060 | leg->Draw(); | |
1061 | ||
1062 | if (1) { | |
1063 | // calculate chi^2 | |
1064 | TF1* flat = new TF1("flat", "1.0", -TMath::PiOver2(), 3*TMath::PiOver2()); | |
1065 | double chi2 = Chi2(hsub, flat); | |
1066 | double chi2_10 = Chi2(hsub10, flat); | |
1067 | int ndof = hsub->GetNbinsX(); | |
1068 | delete flat; | |
1069 | c1->cd(1); | |
1070 | ltx.DrawLatex(0.65, 0.05, Form("#chi^{2}/ndf = %.3g / %d", chi2, ndof-1)); | |
1071 | if (0) | |
1072 | cout << "Chi2[5] = " << chi2 << " Chi2[10] = " << chi2_10 | |
1073 | << " ndof = " << ndof-1 << " (Prob = " << TMath::Prob(chi2,ndof-6) | |
1074 | << " , " << TMath::Prob(chi2_10,ndof-11) << " )" << endl; | |
1075 | } | |
1076 | return c1; | |
1077 | } | |
1078 | ||
1079 | double Chi2(TH1* h, TF1* f, TString opt) { | |
1080 | double chi2 = 0; | |
1081 | double ndf = 0; | |
1082 | for (int i=1;i<=h->GetNbinsX();i++) { | |
1083 | double fy = f->Eval(h->GetBinCenter(i)); | |
1084 | chi2 += pow( (h->GetBinContent(i) - fy)/h->GetBinError(i) , 2.0); | |
1085 | ndf++; | |
1086 | } | |
1087 | ndf -= 1; | |
1088 | ||
1089 | if (opt.Contains("ndf") && ndf > 0) | |
1090 | chi2 /= ndf; | |
1091 | ||
1092 | if (opt.Contains("prob") ) | |
1093 | return TMath::Prob(chi2, ndf); | |
1094 | ||
1095 | return chi2; | |
1096 | } | |
1097 | ||
1098 | double vntvna(double *x, double *par) | |
1099 | { | |
1100 | // This function is called many times during a fit. For a 1-d fit, | |
1101 | // x[0] is the x-position of the nth point, | |
1102 | // i.e. graph->GetX()[n]. For a single-valued TGraph, it can be used | |
1103 | // to indicate which point the fitter is currently working on. The | |
1104 | // par array contains maxpt free parameters that are continually | |
1105 | // updated during the fit process. | |
1106 | ||
1107 | // q is the global index | |
1108 | int q = (int)x[0]; // Assume points are at q + 0.5. | |
1109 | int i = 999; | |
1110 | int j = 999; | |
1111 | int nPoints = (maxpt*(maxpt+1))/2; | |
1112 | TF1::fgRejectPoint = kFALSE; | |
1113 | ||
1114 | // Assign i and j from q | |
1115 | ijFromGlobalIndex(q, i, j); | |
1116 | ||
1117 | if (q < 0 || q >= nPoints) { | |
1118 | Error("vntvna()", | |
1119 | "Problem finding global index. x[0]=%f, q=%d, i=%d, j=%d", | |
1120 | x[0], q, i, j); | |
1121 | ||
1122 | return 999; | |
1123 | } | |
1124 | ||
1125 | // Exclude points from fit. Warning: TF1::fgRejectPoint is a global | |
1126 | // flag checked by all TF1s, beware of side-effects. | |
1127 | if (j < minFitPtBin || j > maxFitPtBin) { | |
1128 | TF1::RejectPoint(); | |
1129 | if (0) | |
1130 | cout << Form("Rejecting x[0],q,i,j = %f %d %d %d", | |
1131 | x[0], q, i, j) | |
1132 | << endl; | |
1133 | return 0; | |
1134 | } | |
1135 | ||
1136 | double v2t = par[i]; | |
1137 | double v2a = par[j]; | |
1138 | int n = par[maxpt]; // Fourier coeff. index | |
1139 | int k = par[maxpt+1]; // cent. bin index | |
1140 | double fitval = v2t * v2a; | |
1141 | ||
1142 | // Add extra fit term for momentum conservation | |
1143 | // It is 1/<\sum pt2> * <ptt><pta> | |
1144 | // From ML & JYO, PRL 106, 102301 (2011) eqs 3-7 | |
1145 | double mptt=0, mpta=0, v1fac=0; | |
1146 | if (n==1) { | |
1147 | mptt = MeanPt(i, j, k, "t"); | |
1148 | mpta = MeanPt(i, j, k, "a"); | |
1149 | v1fac = par[maxpt+2]; | |
1150 | ||
1151 | fitval += v1fac * mptt * mpta; | |
1152 | } | |
1153 | ||
1154 | // Odd vnd's < 0 at high pt. // can I replace the below by saying if | |
1155 | // (vndelta < 0) then make fitval < 0? | |
1156 | if (j >= PtBin(5., 6.) && (n%2) ) | |
1157 | fitval = -fitval; | |
1158 | ||
1159 | return fitval; | |
1160 | } | |
1161 | ||
1162 | void VnDelta(int n, const TH1 &h, double &vnd, double &vnd_err, TString opt) | |
1163 | { | |
1164 | // Extract fourier coefficient n directly from h. If opt is "hi" or | |
1165 | // "lo", the coefficient is calculated for the points + or - their | |
1166 | // errors. "hi" and "lo" are meant for an h with systematic (not | |
1167 | // statistical) errors, so vnd_err is not calculated. | |
1168 | double VN = 0.0; | |
1169 | double norm = 0.0; | |
1170 | ||
1171 | // Calculate and set coeff. vnd | |
1172 | for (int ib=1; ib<=h.GetNbinsX(); ib++) { | |
1173 | double deltaphi = h.GetBinCenter(ib); | |
1174 | double weight = h.GetBinContent(ib); | |
1175 | ||
1176 | if (opt.Contains("hi")) | |
1177 | weight += h.GetBinError(ib); | |
1178 | else if (opt.Contains("lo")) | |
1179 | weight -= h.GetBinError(ib); | |
1180 | ||
1181 | if (opt.Contains("sin") ) | |
1182 | VN += sin(n * deltaphi ) * weight; | |
1183 | else | |
1184 | VN += cos(n * deltaphi ) * weight; | |
1185 | ||
1186 | norm += weight; | |
1187 | } | |
1188 | ||
1189 | if (norm==0) | |
1190 | ::Error("VnDelta()", "Div/0 error"); | |
1191 | ||
1192 | VN = VN / norm; | |
1193 | ||
1194 | // // Mom. cons correction????????????????????????????????? | |
1195 | // if (n==1) VN += 0.002*ptmean[i]*ptmean[j]; | |
1196 | ||
1197 | vnd = VN; | |
1198 | if (opt.Contains("hi") || opt.Contains("lo")) { | |
1199 | vnd_err = 0; | |
1200 | return; | |
1201 | } | |
1202 | ||
1203 | // Statistical uncertainty vnd_err | |
1204 | double quad_sum_uncertainty = 0.0; | |
1205 | for (int ib=1; ib<=h.GetNbinsX(); ib++) { | |
1206 | double deltaphi = h.GetBinCenter(ib); | |
1207 | double dfdwi = cos(n * deltaphi)/norm - VN/norm; | |
1208 | ||
1209 | if (opt.Contains("sin") ) | |
1210 | dfdwi = sin(n * deltaphi)/norm - VN/norm; | |
1211 | ||
1212 | double sigma = h.GetBinError(ib); | |
1213 | quad_sum_uncertainty += dfdwi*dfdwi * sigma*sigma; | |
1214 | } | |
1215 | ||
1216 | vnd_err = TMath::Sqrt(quad_sum_uncertainty); | |
1217 | return; | |
1218 | } | |
1219 | ||
1220 | void VnDelta(int n, const TF1 &f, double &vnd, double &vnd_err, TString opt) | |
1221 | { | |
1222 | // Extract fourier coefficient n from f. | |
1223 | double VN = 0.0; | |
1224 | double norm = 0.0; | |
1225 | int nSteps = 100; | |
1226 | double x1=0, x2=0; | |
1227 | f.GetRange(x1, x2); | |
1228 | double dx = (x2-x1)/nSteps; | |
1229 | ||
1230 | if (!opt.IsNull() ) cout<<opt.Data()<<endl; | |
1231 | ||
1232 | // Calculate and set coeff. vnd | |
1233 | for (int ib=0; ib<nSteps; ib++) { | |
1234 | ||
1235 | double deltaphi = x1 + ib*dx; | |
1236 | double weight = f.Eval(deltaphi); | |
1237 | ||
1238 | VN += cos(n * deltaphi ) * weight; | |
1239 | norm += weight; | |
1240 | } | |
1241 | ||
1242 | if (norm==0) | |
1243 | Error("VnDelta()", "Div/0 error"); | |
1244 | ||
1245 | VN = VN / norm; | |
1246 | vnd = VN; | |
1247 | vnd_err = 0; // for now | |
1248 | return; | |
1249 | } | |
1250 | ||
1251 | ||
1252 | ||
1253 | // The argument ydiv is the division between pads as a fraction of the | |
1254 | // height. | |
1255 | TCanvas* TwoPanelVert(double ydiv, const char* canvName, const char* canvTitle) | |
1256 | { | |
1257 | TCanvas* c = new TCanvas(canvName, canvTitle, 800, 800); | |
1258 | // Extents of pads inside canvas | |
1259 | double x1=0.01, x2=0.99, y1=0.01, y2=0.99; | |
1260 | double lm=0.2, rm=0.02; | |
1261 | ||
1262 | // Division between top and bottom pads | |
1263 | double ysplit = y1 + ydiv*(y2-y1); | |
1264 | ||
1265 | // Divide first, adjust afterward. | |
1266 | c->Divide(1,2); | |
1267 | TPad* ptop = (TPad*)(c->GetPad(1)); | |
1268 | TPad* pbot = (TPad*)(c->GetPad(2)); | |
1269 | ||
1270 | ptop->SetPad(x1, ysplit, x2, y2); | |
1271 | ptop->SetMargin(lm, rm, 0.01, 0.01); // left, right, bottom, top. | |
1272 | ptop->SetFrameLineWidth(2); | |
1273 | ||
1274 | pbot->SetPad(x1, y1, x2, ysplit); | |
1275 | pbot->SetMargin(lm, rm, 0.4, 0.01); // left, right, bottom, top. | |
1276 | pbot->SetFrameLineWidth(2); | |
1277 | pbot->SetTicks(); | |
1278 | ||
1279 | c->cd(); | |
1280 | return c; | |
1281 | } | |
1282 | ||
1283 | char* trigLabel(int i) | |
1284 | { | |
1285 | if (i<0 || i>maxpt) | |
1286 | return Form("Error: no trig pt bin %d", i); | |
1287 | ||
1288 | double x1 = gi->GetX()[i] - gi->GetEX()[i]; | |
1289 | double x2 = gi->GetX()[i] + gi->GetEX()[i]; | |
1290 | ||
1291 | // int prec1 = 0, prec2 = 0; | |
1292 | // if (x1 - (int)x1 > 0) | |
1293 | // prec1 = 1; | |
1294 | // if (x2 - (int)x2 > 0) | |
1295 | // prec2 = 1; | |
1296 | ||
1297 | return Form("%.3g-%.3g", x1, x2); | |
1298 | } | |
1299 | ||
1300 | char* asscLabel(int j) | |
1301 | { | |
1302 | if (j<0 || j>maxpt) | |
1303 | return Form("Error: no assc pt bin %d", j); | |
1304 | ||
1305 | double x1 = gj->GetX()[j] - gj->GetEX()[j]; | |
1306 | double x2 = gj->GetX()[j] + gj->GetEX()[j]; | |
1307 | ||
1308 | /* Stupid coding... | |
1309 | int prec1 = 0, prec2 = 0; | |
1310 | ||
1311 | if (x1 - (int)x1 > 0) | |
1312 | prec1 = 1; | |
1313 | if (x2 - (int)x2 > 0) | |
1314 | prec2 = 1; | |
1315 | */ | |
1316 | ||
1317 | return Form("%.3g-%.3g", x1, x2); | |
1318 | } | |
1319 | ||
1320 | char* centLabel(int k) | |
1321 | { | |
1322 | if (k<0 || k>maxcent) | |
1323 | return Form("Error: no cent bin %d", k); | |
1324 | ||
1325 | double c1 = gk->GetX()[k] - gk->GetEX()[k]; | |
1326 | double c2 = gk->GetX()[k] + gk->GetEX()[k]; | |
1327 | ||
1328 | return Form("%.0f-%.0f%%", c1, c2); | |
1329 | } | |
1330 | ||
1331 | ||
1332 | void VnLimits(double &minVn_k, double &maxVn_k, int k, int imin, int imax) | |
1333 | { | |
1334 | // Find limits for each centrality bin k between imin (inclusive) | |
1335 | // and imax (also inclusive) to plot a common y axis. | |
1336 | // for (int k=0; k < maxcent; k++) { | |
1337 | for (int n=1; n<=5; n++) { | |
1338 | for (int i=imin; i<=imax; i++) { | |
1339 | TGraphErrors* gc = VnDeltaVsPtAssoc(i, k, n); | |
1340 | if (!gc) { | |
1341 | Error("VnLimits()", "No graph i %d k %d n %d", i,k,n ); | |
1342 | continue; | |
1343 | } | |
1344 | for (int j=0; j<gc->GetN(); j++) { | |
1345 | double vn = gc->GetY()[j]; | |
1346 | double evn = gc->GetEY()[j]; | |
1347 | if (vn+evn > maxVn_k) maxVn_k = vn+evn; | |
1348 | if (vn-evn < minVn_k) minVn_k = vn-evn; | |
1349 | } | |
1350 | } | |
1351 | } | |
1352 | ||
1353 | return; | |
1354 | } | |
1355 | ||
1356 | TCanvas* DrawQ(int k, int n, TString opt) | |
1357 | { | |
1358 | TStopwatch ts1, ts2; | |
1359 | int fitCurveColor = kRed-4; | |
1360 | int divColor = kBlue; | |
1361 | const char* region = "RIDGE"; | |
1362 | if (opt.Contains("ALL")) | |
1363 | region = "ALL"; | |
1364 | ||
1365 | TString title = Form("globalfit_%d_cent%dto%d", | |
1366 | n, centlow[k],centhigh[k]); | |
1367 | TString name = Form("V%dvsQ_etamin%02d_cent%dto%d", | |
1368 | n, (int)(10*minRidgeDEta), | |
1369 | centlow[k],centhigh[k]); | |
1370 | ||
1371 | if (!opt.IsNull()) { | |
1372 | name.Append("_"); | |
1373 | name.Append(opt.Data()); | |
1374 | title.Append("_"); | |
1375 | title.Append(opt.Data()); | |
1376 | } | |
1377 | ||
1378 | TCanvas* c = TwoPanelVert(0.4, name.Data(), title.Data()); | |
1379 | int cwidth = opt.Contains("highptfit") ? 500 : 1200; | |
1380 | c->SetCanvasSize(cwidth, 450); | |
1381 | c->SetWindowSize(cwidth+50, 500); | |
1382 | ||
1383 | if (opt.Contains("highptfit")) { | |
1384 | c->GetPad(1)->SetLeftMargin(0.2); | |
1385 | c->GetPad(2)->SetLeftMargin(0.2); | |
1386 | c->GetPad(1)->SetRightMargin(0.05); | |
1387 | c->GetPad(2)->SetRightMargin(0.05); | |
1388 | ||
1389 | } | |
1390 | ||
1391 | TLegend* lq = new TLegend(0.14, 0.56, 0.25, 0.83); | |
1392 | lq->SetFillColor(kNone); | |
1393 | lq->SetBorderSize(0); | |
1394 | ||
1395 | // g: global VnDelta vs global index (stat errors). | |
1396 | // gs: global VnDelta vs global index (sys errors). | |
1397 | // fn: global fit TF1. | |
1398 | // r: ratio of g/fit. | |
1399 | TString opt1 = ""; | |
1400 | if (opt.Contains("ALL")) | |
1401 | opt1.Append("_ALL"); | |
1402 | if (opt.Contains("ptcons")) | |
1403 | opt1.Append("_ptcons"); | |
1404 | ||
1405 | TGraphErrors* g = VnDVsQ(n, k, region, cfDef, opt1); | |
1406 | opt1.Append("_sys"); | |
1407 | TGraphErrors* gs = VnDVsQ(n, k, region, cfDef, opt1); | |
1408 | ||
1409 | int ptaLo = 0, ptaHi = 999; | |
1410 | if (opt.Contains("lowptfit")) { | |
1411 | ptaHi = PtBin(3., 4.); | |
1412 | } | |
1413 | if (opt.Contains("highptfit")) { | |
1414 | ptaLo = PtBin(5., 6.); | |
1415 | } | |
1416 | // If custom pta range was passed in, e.g. | |
1417 | // Form("ptabin%dto%d", PtBin(1.5, 2.0), PtBin(8, 15)) | |
1418 | if (opt.Contains("ptabin")) { | |
1419 | TRegexp re("ptabin[0-9]+to[0-9]+"); | |
1420 | TRegexp re1("ptabin[0-9]+"); | |
1421 | TRegexp re2("to[0-9]+"); | |
1422 | TString rx = opt(re), pta1s = rx(re1), pta2s = rx(re2); | |
1423 | pta1s.ReplaceAll("ptabin",""); | |
1424 | pta2s.ReplaceAll("to",""); | |
1425 | ptaLo = pta1s.Atoi(); | |
1426 | ptaHi = pta2s.Atoi(); | |
1427 | } | |
1428 | ||
1429 | TF1* fn = GlobalFitVsQ(n, k, ptaLo, ptaHi, region, cfDef, ""); | |
1430 | TF1* fn_hi = GlobalFitVsQ(n, k, ptaLo, ptaHi, region, cfDef, "hi_sys"); | |
1431 | TF1* fn_lo = GlobalFitVsQ(n, k, ptaLo, ptaHi, region, cfDef, "lo_sys"); | |
1432 | ||
1433 | // Systematic band on ratio line at unity | |
1434 | TH1F *r_sys = new TH1F(Form("r_sys_%s", title.Data()), "r_sys", | |
1435 | g->GetN()-0.5, 0, g->GetN()-0.5); | |
1436 | utils.set_hist_props(r_sys, kBlack, kGray, kBlack, kDot, 1.0); | |
1437 | ||
1438 | TGraphErrors* r = new TGraphErrors(); | |
1439 | TString rname(g->GetName()); | |
1440 | rname.Append("_ratio"); | |
1441 | r->SetName(rname.Data()); | |
1442 | ||
1443 | if (!fn) { | |
1444 | Error("DrawQ()","Problem getting/creating global fit function"); | |
1445 | gSystem->Exit(-1); | |
1446 | } | |
1447 | ||
1448 | // Set visual properties | |
1449 | utils.set_tgraph_props(g, kBlack, kBlack, kFullSquare, 1.0); | |
1450 | utils.set_tgraph_props(gs, kBlack, kBlack, kDot, 0.5); | |
1451 | gs->SetFillColor(kGray); | |
1452 | utils.set_tgraph_props(r, kBlack, kBlack, kFullCircle, 1.0); | |
1453 | g->SetLineWidth(1); | |
1454 | r->SetLineWidth(1); | |
1455 | fn->SetLineWidth(1); | |
1456 | fn->SetLineColor(fitCurveColor); | |
1457 | fn_hi->SetLineStyle(kDashed); | |
1458 | fn_lo->SetLineStyle(kDashed); | |
1459 | fn_hi->SetLineWidth(1); | |
1460 | fn_lo->SetLineWidth(1); | |
1461 | ||
1462 | // Compute ratio for r graph | |
1463 | for (int q=0; q<g->GetN(); q++) { | |
1464 | int ii=0, jj=0; | |
1465 | ijFromGlobalIndex(q, ii, jj); | |
1466 | ||
1467 | double fy = fn->Eval(q+0.5); | |
1468 | double y = fy==0? 1e9 : g->GetY()[q]/fy; | |
1469 | double rerr = fn_hi->Eval(q+0.5) / fy - 1; | |
1470 | ||
1471 | if (opt.Contains("highptfit") && n%2 && g->GetY()[q]!=0) | |
1472 | y = fy/g->GetY()[q]; | |
1473 | ||
1474 | // No error band when no fit point! | |
1475 | if (jj < ptaLo || jj > ptaHi) | |
1476 | rerr = 0; | |
1477 | ||
1478 | r->SetPoint(q, g->GetX()[q], y); | |
1479 | r->SetPointError(q, 0.0, TMath::Abs(g->GetEY()[q]/fy)); | |
1480 | ||
1481 | r_sys->SetBinContent(q+1, 1.0); | |
1482 | r_sys->SetBinError(q+1, rerr); | |
1483 | } | |
1484 | ||
1485 | // Set up plot ranges for top panel (g and f) | |
1486 | int nx = g->GetN(); | |
1487 | double x1,y1,x2,y2; | |
1488 | double sf = 1; // scale factor | |
1489 | fn->GetRange(x1,x2); | |
1490 | if (opt.Contains("highptfit")) { | |
1491 | x1 = 45; | |
1492 | } | |
1493 | y1 = sf*fn->GetMinimum(x1,x2); | |
1494 | y2 = sf*fn->GetMaximum(x1,x2); | |
1495 | double marg = 0.4*(y2-y1); | |
1496 | // if (opt.Contains("lowptfit") || opt.Contains("highptfit") ) marg = 2*(y2-y1); | |
1497 | x1 -= 0.5; | |
1498 | x2 += 0.5; | |
1499 | y1 -= marg; | |
1500 | y2 += marg; | |
1501 | ||
1502 | // Set up plot ranges for lower panel (r) | |
1503 | double ry1=0.21, ry2=1.79; | |
1504 | if (opt.Contains("highptfit")) { | |
1505 | ry1 = 0.01; | |
1506 | ry2 = 1.99; | |
1507 | } | |
1508 | ||
1509 | // Get / make frame histos for top and bottom | |
1510 | TH2F* hf1 = 0; // (TH2F*)gDirectory->Get("qhf1"); | |
1511 | TH2F* hf2 = 0; // (TH2F*)gDirectory->Get("qhf2"); | |
1512 | TString name1 = Form("qhf1_%s", title.Data()); | |
1513 | TString name2 = Form("qhf2_%s", title.Data()); | |
1514 | // TString sv = n==1 ? "#LTcos(#Delta#phi)#GT" : | |
1515 | // Form("#LTcos(%d#Delta#phi)#GT", n); | |
1516 | TString sv = Form("V_{%d#Delta}", n); | |
1517 | TString sgf = Form("(v_{%d}^{t}v_{%d}^{a})_{fit}", n, n); | |
1518 | ||
1519 | if (!hf1) | |
1520 | hf1 = new TH2F(name1.Data(), "hf1", nx+1, x1, x2, 1000, y1, y2); | |
1521 | if (!hf2) | |
1522 | hf2 = new TH2F(name2.Data(), "hf2", nx+1, x1, x2, 1000, ry1, ry2); | |
1523 | ||
1524 | // Axes | |
1525 | TAxis* ax1 = hf1->GetXaxis(); | |
1526 | TAxis* ay1 = hf1->GetYaxis(); | |
1527 | TAxis* ax2 = hf2->GetXaxis(); | |
1528 | TAxis* ay2 = hf2->GetYaxis(); | |
1529 | hf1->SetBins(nx+1, x1, x2, 100, y1, y2); | |
1530 | hf2->SetBins(nx+1, x1, x2, 100, ry1, ry2); | |
1531 | hf1->SetAxisRange(x1, x2, "X"); | |
1532 | hf1->SetAxisRange(y1, y2, "Y"); | |
1533 | hf2->SetAxisRange(x1, x2, "X"); | |
1534 | hf2->SetAxisRange(ry1, ry2, "Y"); | |
1535 | ax1->SetLabelOffset(0.02); | |
1536 | ax1->SetLabelSize(0.04); ay1->SetLabelSize(0.08); | |
1537 | ax2->SetLabelSize(0.18); ay2->SetLabelSize(0.12); | |
1538 | ax1->SetNdivisions(1); ax2->SetNdivisions(1); | |
1539 | ay1->SetNdivisions(208); ay2->SetNdivisions(opt.Contains("highptfit")?104:210); | |
1540 | ax1->SetTickLength(0.02); | |
1541 | ||
1542 | //ay1->SetTitle(Form("%s and #color[%d]{%s}",sv.Data(), fitCurveColor, sgf.Data())); | |
1543 | ay1->SetTitle(Form("%s", sv.Data())); | |
1544 | ay2->SetTitle(Form("#frac{%s}{fit}", sv.Data())); | |
1545 | // ay2->SetTitle(Form("#frac{%s}{%s}", sv.Data(), sgf.Data())); | |
1546 | // ax2->SetTitle(Form("#splitline{associated p_{T}^{a}}{#color[%d]{ trigger p_{T}^{t}}} [GeV/c]", divColor)); | |
1547 | ax2->SetTitle(Form("#color[%d]{p_{T}^{t}}, p_{T}^{a} [GeV/c]", divColor)); | |
1548 | ay1->SetTitleOffset(opt.Contains("highptfit") ? 1. : 0.45); | |
1549 | ay1->SetTitleSize(0.09); | |
1550 | ay1->CenterTitle(); | |
1551 | ay2->SetTitleOffset(opt.Contains("highptfit") ? 0.6 : 0.26); | |
1552 | ay2->SetTitleSize(0.14); | |
1553 | ay2->CenterTitle(); | |
1554 | ax2->SetLabelOffset(2.1); | |
1555 | ax2->SetTitleOffset(1.3); | |
1556 | ax2->SetTitleSize(0.14); | |
1557 | ax2->CenterTitle(); | |
1558 | ||
1559 | double leftMarg = opt.Contains("highptfit") ? 0.2 : 0.08; | |
1560 | c->cd(1); | |
1561 | gPad->SetFrameLineWidth(1); | |
1562 | gPad->SetLeftMargin(leftMarg); | |
1563 | hf1->Draw(); | |
1564 | TLine zero; | |
1565 | zero.DrawLine(x1, 0., x2, 0.); | |
1566 | ||
1567 | ||
1568 | // NEW======================================== | |
1569 | for (int i=0; i<gi->GetN(); i++) { | |
1570 | int q1 = GlobalIndex(i, 1), q2 = GlobalIndex(i, i); | |
1571 | TF1* ifn = 0, *ifn_hi = 0, *ifn_lo = 0; | |
1572 | // TF1* _ifn = 0, *_ifn_hi = 0, *_ifn_lo = 0; // temp, unscaled. | |
1573 | ||
1574 | // This works, but has lines between \ptt groups :( | |
1575 | // TH1F *hfn=0, *hfn_hi=0, *hfn_lo=0; | |
1576 | // int nq = g->GetN(); | |
1577 | // hfn = new TH1F(Form("%s_%d", fn->GetName(), i), Form("%s_%d", fn->GetName(), i), | |
1578 | // nq, 0, nq); | |
1579 | // hfn->SetLineColor(kMagenta); | |
1580 | // for (int q=0; q<g->GetN(); q++) { | |
1581 | // hfn->SetBinContent(q+1, sf*fn->Eval(q+0.5)); | |
1582 | // } | |
1583 | ||
1584 | ifn = (TF1*) fn->Clone(Form("%s_%d", fn->GetName(), i)); | |
1585 | ifn_hi = (TF1*) fn_hi->Clone(Form("%s_%d", fn_hi->GetName(), i)); | |
1586 | ifn_lo = (TF1*) fn_lo->Clone(Form("%s_%d", fn_lo->GetName(), i)); | |
1587 | ||
1588 | ifn->SetName(Form("%s_%d", fn->GetName(), i)); | |
1589 | ifn_hi->SetName(Form("%s_%d", fn_hi->GetName(), i)); | |
1590 | ifn_lo->SetName(Form("%s_%d", fn_lo->GetName(), i)); | |
1591 | ||
1592 | double r1 = q1-0.9, r2 = q2+0.92; | |
1593 | ifn->SetRange(r1, r2); | |
1594 | ifn_hi->SetRange(r1, r2); | |
1595 | ifn_lo->SetRange(r1, r2); | |
1596 | ||
1597 | ||
1598 | // for (int ip=0; ip<ifn->GetNpar(); ip++) { | |
1599 | // double par = ifn->GetParameter(ip); | |
1600 | // ifn->SetParameter(ip, 1000*par); | |
1601 | // cout<<par<< " " << ifn->GetParameter(ip) << endl; | |
1602 | // } | |
1603 | ||
1604 | ||
1605 | // ifn = new TF1(Form("%s_%d", fn->GetName(), i), Form("100*%s", _ifn->GetName())); | |
1606 | // ifn_hi = new TF1(Form("%s_%d", fn_hi->GetName(), i), Form("100*%s", _ifn_hi->GetName())); | |
1607 | // ifn_lo = new TF1(Form("%s_%d", fn_lo->GetName(), i), Form("100*%s", _ifn_lo->GetName())); | |
1608 | ||
1609 | ifn_hi->Draw("same"); | |
1610 | ifn_lo->Draw("same"); | |
1611 | ifn->Draw("same"); | |
1612 | // hfn->Draw("same"); | |
1613 | } | |
1614 | // NEW======================================== | |
1615 | ||
1616 | for (int i=0; i<gi->GetN(); i++) { | |
1617 | char* newName = Form("%s_ptt%d", g->GetName(), i); | |
1618 | TGraphErrors* gcol = (TGraphErrors*)g->Clone(newName); | |
1619 | int col = PtColor(ptlow[i], pthigh[i]); | |
1620 | utils.set_tgraph_props(gcol, kBlack, col, kFullSquare, 1.0); | |
1621 | gcol->SetLineWidth(1); | |
1622 | ||
1623 | newName = Form("%s_ptt%d", gs->GetName(), i); | |
1624 | TGraphErrors* g_sy = (TGraphErrors*)gs->Clone(newName); | |
1625 | int col_sy = PtColorPale(ptlow[i], pthigh[i]); | |
1626 | utils.set_tgraph_props(g_sy, kBlack, col_sy, kDot, 0.5); | |
1627 | g_sy->SetFillColor(col_sy); | |
1628 | ||
1629 | for (int q=0; q<gcol->GetN(); q++) { | |
1630 | int i2=0, j=0; | |
1631 | ijFromGlobalIndex(q,i2,j); | |
1632 | if (i2!=i) { | |
1633 | gcol->SetPoint(q, gcol->GetX()[q], -99.); | |
1634 | g_sy->SetPoint(q, gcol->GetX()[q], -99.); | |
1635 | } | |
1636 | ||
1637 | // Scale points by sf | |
1638 | gcol->SetPoint(q, gcol->GetX()[q], sf*gcol->GetY()[q]); | |
1639 | g_sy->SetPoint(q, g_sy->GetX()[q], sf*g_sy->GetY()[q]); | |
1640 | gcol->SetPointError(q, gcol->GetEX()[q], sf*gcol->GetEY()[q]); | |
1641 | g_sy->SetPointError(q, g_sy->GetEX()[q], sf*g_sy->GetEY()[q]); | |
1642 | ||
1643 | } | |
1644 | g_sy->Draw("e2psame"); | |
1645 | gcol->Draw("ezpsame"); | |
1646 | } | |
1647 | ||
1648 | c->cd(2); | |
1649 | gPad->SetFrameLineWidth(1); | |
1650 | gPad->SetLeftMargin(leftMarg); | |
1651 | hf2->Draw(); | |
1652 | r_sys->Draw("e2psame"); | |
1653 | ||
1654 | // Overlay lines, arrows, text, etc. | |
1655 | TLine div, div_a; // division ticks for trigger and assoc | |
1656 | TArrow a; | |
1657 | double tipSize = 0.008; // size of arrowhead | |
1658 | a.SetLineWidth(1); | |
1659 | div.SetLineWidth(2); | |
1660 | div.SetLineColor(fitCurveColor); | |
1661 | div.DrawLine(x1, 1.0, x2, 1.0); | |
1662 | r->Draw("ezpsame"); | |
1663 | gList->Add(r); | |
1664 | ||
1665 | for (int i=0; i<gi->GetN(); i++) { | |
1666 | char* newName = Form("%s_ptt%d", r->GetName(), i); | |
1667 | TGraphErrors* rcol = (TGraphErrors*)r->Clone(newName); | |
1668 | int col = PtColor(ptlow[i], pthigh[i]); | |
1669 | utils.set_tgraph_props(rcol, kBlack, col, kFullCircle, 1.0); | |
1670 | for (int q=0; q<rcol->GetN(); q++) { | |
1671 | int i2=0, j=0; | |
1672 | ijFromGlobalIndex(q,i2,j); | |
1673 | if (i2!=i) rcol->SetPoint(q, rcol->GetX()[q], -99.); | |
1674 | } | |
1675 | rcol->SetLineWidth(1); | |
1676 | rcol->Draw("ezpsame"); | |
1677 | } | |
1678 | ||
1679 | div.SetLineColor(divColor); | |
1680 | div.SetLineStyle(kSolid); | |
1681 | ||
1682 | // Label centrality and fourier index | |
1683 | c->cd(1); | |
1684 | ltx.SetTextSize(0.14); | |
1685 | ltx.DrawLatex(opt.Contains("highptfit")?0.25:0.14, 0.85, Form("n = %d", n)); | |
1686 | ||
1687 | if (!opt.Contains("highptfit")) { | |
1688 | if (n<=2) { | |
1689 | ltx.SetTextSize(0.1); | |
1690 | if (ispp) | |
1691 | ltx.DrawLatex(0.14, 0.70, Form("pp 2.76 TeV")); | |
1692 | else { | |
1693 | ltx.DrawLatex(0.14, 0.70, Form("Pb-Pb 2.76 TeV")); | |
1694 | ltx.DrawLatex(0.14, 0.61, Form("%d-%d%%", centlow[k],centhigh[k])); | |
1695 | } | |
1696 | if (opt.Contains("ptcons")) { | |
1697 | ltx.DrawLatex(0.14, 0.58, Form("#LT#sum p_{T}^{2}#GT = 4.76 GeV^{2}")); | |
1698 | } | |
1699 | ||
1700 | } | |
1701 | if (n==3) { | |
1702 | lq->AddEntry(g, Form("V_{n#Delta}"), "epl"); | |
1703 | lq->AddEntry(fn, Form("fit"), "l"); | |
1704 | lq->Draw(); | |
1705 | } | |
1706 | // Always show | |
1707 | ltx.SetTextSize(0.07); | |
1708 | double ptaFitLow = ptlow[ptaLo]; | |
1709 | double ptaFitHigh = ptaHi<999? pthigh[ptaHi] : pthigh[maxpt-1]; | |
1710 | // if (!opt.Contains("highptfit")) | |
1711 | ltx.DrawLatex(0.14, 0.09, | |
1712 | Form("Global fit range: %3.2g < p_{T}^{a} < min(p_{T}^{t}, %3.2g GeV/c)", | |
1713 | ptaFitLow, ptaFitHigh)); | |
1714 | } | |
1715 | ||
1716 | for (int q=0; q<g->GetN(); q++) { | |
1717 | int i=0, j=0; | |
1718 | ijFromGlobalIndex(q,i,j); | |
1719 | double drop = 0.35; | |
1720 | double ptt = gi->GetX()[i] + gi->GetEX()[i]; | |
1721 | double pta = gj->GetX()[j] + gj->GetEX()[j]; | |
1722 | double tickLength = 0.03*(y2-y1); | |
1723 | ||
1724 | // Draw divisions between trigger pt intervals... | |
1725 | bool skipTicks = opt.Contains("highptfit") && (pta < 5.); | |
1726 | if (i==j && !skipTicks) { | |
1727 | // if (ptt > 5) { | |
1728 | c->cd(1); | |
1729 | div.DrawLine(q+1, y1, q+1, y1+tickLength); | |
1730 | c->cd(2); | |
1731 | tickLength = drop-0.1; | |
1732 | div.DrawLine(q+1, ry2-tickLength, q+1, ry2+tickLength); | |
1733 | div.DrawLine(q+1, ry1-tickLength, q+1, ry1+tickLength); | |
1734 | // } | |
1735 | ||
1736 | // Alignment: align = 10*HorizontalAlign + VerticalAlign | |
1737 | // For horizontal alignment the following convention applies: | |
1738 | // 1=left adjusted, 2=centered, 3=right adjusted | |
1739 | // For vertical alignment the following convention applies: | |
1740 | // 1=bottom adjusted, 2=centered, 3=top adjusted | |
1741 | ltx2.SetTextAlign(23); | |
1742 | ltx2.SetTextSize(0.12); | |
1743 | ltx2.SetTextColor(divColor); | |
1744 | ||
1745 | if (ptt != 0.75) // too cluttered otherwise | |
1746 | ltx2.DrawLatex(q+1, ry1-drop, Form("%.2g", ptt) ); // pt_trig labels | |
1747 | } | |
1748 | ||
1749 | if (pta < ptt && (pta == (int)pta || pta == 0.5) && !skipTicks) { | |
1750 | div_a.DrawLine(q+1, ry1-drop/8, q+1, ry1+drop/4); | |
1751 | ltx2.SetTextColor(kBlack); | |
1752 | ltx2.SetTextFont(52); // italic | |
1753 | ltx2.DrawLatex(q+1, ry1-drop/4, Form("%.2g", pta) ); // pt_assoc labels | |
1754 | ltx2.SetTextColor(divColor); | |
1755 | ltx2.SetTextFont(62); // 62 is default | |
1756 | } | |
1757 | ||
1758 | // Draw arrows to off-scale points | |
1759 | c->cd(1); | |
1760 | double y = g->GetY()[q]; | |
1761 | double ymax = hf1->GetYaxis()->GetXmax(); | |
1762 | double ymin = hf1->GetYaxis()->GetXmin(); | |
1763 | double len = 0.06*(ymax-ymin); | |
1764 | double gap = 0.01*(ymax-ymin); | |
1765 | ||
1766 | if (y > ymax) | |
1767 | a.DrawArrow(q+0.5, ymax-len-gap, q+0.5, ymax-gap, tipSize, ">"); | |
1768 | if (y < ymin) | |
1769 | a.DrawArrow(q+0.5, ymin+len+gap, q+0.5, ymin+gap, tipSize, ">"); | |
1770 | ||
1771 | // Do the same for the ratio. | |
1772 | // If fit has not been applied to some points, don't draw the arrow. | |
1773 | c->cd(2); | |
1774 | // gPad->SetTickx(); | |
1775 | y = r->GetY()[q]; | |
1776 | ymax = hf2->GetYaxis()->GetXmax(); | |
1777 | ymin = hf2->GetYaxis()->GetXmin(); | |
1778 | len = 0.06/0.4*(ymax-ymin); // 0.4 from TwoPanelVert ctor | |
1779 | gap = 0.02*(ymax-ymin); | |
1780 | ||
1781 | if (0) { // Draw at upper or lower edge of panel | |
1782 | if (y > ymax) | |
1783 | a.DrawArrow(q+0.5, ymax-len-gap, q+0.5, ymax-gap, tipSize, ">"); | |
1784 | if (y < ymin) | |
1785 | a.DrawArrow(q+0.5, ymin+len+gap, q+0.5, ymin+gap, tipSize, ">"); | |
1786 | } | |
1787 | // Draw with arrow based from ref. line. | |
1788 | // No arrow when no fit point! | |
1789 | if (j < ptaLo || j > ptaHi) | |
1790 | continue; | |
1791 | if (y > ymax) | |
1792 | a.DrawArrow(q+0.5, 1.0, q+0.5, 1.0+len, tipSize, ">"); | |
1793 | if (y < ymin) | |
1794 | a.DrawArrow(q+0.5, 1.0, q+0.5, 1.0-len, tipSize, ">"); | |
1795 | } // end q loop | |
1796 | ||
1797 | // ax2->LabelsOption("h"); | |
1798 | c->cd(); | |
1799 | c->Update(); | |
1800 | return c; | |
1801 | } | |
1802 | ||
1803 | TCanvas* DrawVn(int k, int ptt1, int ptt2) | |
1804 | { | |
1805 | // Draw Fourier coefficients on one canvas for each centrality bin k | |
1806 | // for trigger pt bins from ptt1 to (and including) ptt2. | |
1807 | // The x-axis is associated pt. | |
1808 | const char* name = Form("vndelta_etamin%02d_cent%dto%d_trig%dto%d", | |
1809 | (int)(10*minRidgeDEta), centlow[k],centhigh[k], ptt1, ptt2); | |
1810 | TCanvas* cv = new TCanvas(name, name, 1200, 500); | |
1811 | TLegend* lv = new TLegend(0.09, 0.72, 0.27, 0.88); // w.r.t. cv | |
1812 | lv->SetFillColor(kNone); | |
1813 | lv->SetBorderSize(0); | |
1814 | ||
1815 | double lv2bottom = ptt2<5? 0.65 : 0.8; | |
1816 | TLegend* lv2 = new TLegend(0.82, lv2bottom, 0.94, 0.95, "p_{T}^{t} (GeV/c)"); | |
1817 | lv2->SetFillColor(kNone); | |
1818 | lv2->SetBorderSize(0); | |
1819 | utils.padsetup(cv, 5, 1, "", 0.12, 0.01, 0.03, 0.2); | |
1820 | ||
1821 | // Find limits | |
1822 | double minVn_k=1000, maxVn_k=-1000; | |
1823 | VnLimits(minVn_k, maxVn_k, k, ptt1, ptt2); | |
1824 | ||
1825 | for (int n=1; n<=5; n++) { | |
1826 | cv->cd(n); | |
1827 | double vertMargin = 0.1*(maxVn_k - minVn_k); | |
1828 | double xmax = gj->GetX()[ptt2] + gj->GetEX()[ptt2]; | |
1829 | double xmin = -0.43; | |
1830 | TH1F* hf = gPad->DrawFrame(xmin, minVn_k - vertMargin, | |
1831 | xmax, maxVn_k + vertMargin); | |
1832 | int nXticks = 404; | |
1833 | int nYticks = 210; | |
1834 | if (n==1) { | |
1835 | utils.make_nice_axes(cv, hf, 2., nXticks, nYticks, -0.02, 0.02); | |
1836 | hf->SetLabelOffset(-0.01, "X"); | |
1837 | } | |
1838 | else { | |
1839 | utils.make_nice_axes(cv, hf, 3.4, nXticks, nYticks, -0.075, 0.0); | |
1840 | hf->SetLabelOffset(-0.057, "X"); | |
1841 | hf->SetTickLength(0.02, "X"); | |
1842 | hf->SetTickLength(0.05, "Y"); | |
1843 | } | |
1844 | hf->GetXaxis()->SetTicks("+-"); | |
1845 | hf->GetYaxis()->SetTicks("+-"); | |
1846 | ||
1847 | // Loop backwards to improve visibility of low-pt points | |
1848 | for (int i=ptt2; i>=ptt1; i--) { | |
1849 | ||
1850 | TGraphErrors *gc=0; | |
1851 | gc = VnDeltaVsPtAssoc(i, k, n); | |
1852 | if (!gc) { | |
1853 | Error("DrawVn()", "Problem finding graph n%d, i%d, k%d" ,n,i,k); | |
1854 | continue; | |
1855 | } | |
1856 | ||
1857 | utils.set_tgraph_props(gc, colorsDark[i], colorsDark[i], kFullSquare, 1.6); | |
1858 | ||
1859 | TGraphErrors* gfit = GlobalFitVsPtAssoc(i, k, n); | |
1860 | ||
1861 | gfit->Draw("plsame"); | |
1862 | gc->Draw("epsame"); | |
1863 | ||
1864 | if (n==1) | |
1865 | lv2->AddEntry(gc,trigLabel(i), "p"); | |
1866 | ||
1867 | } | |
1868 | } | |
1869 | ||
1870 | cv->cd(); | |
1871 | TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1); | |
1872 | overlay->SetFillStyle(4000); // transparent | |
1873 | overlay->Draw(); | |
1874 | overlay->cd(); | |
1875 | ||
1876 | ltx.SetTextSize(0.07); | |
1877 | lv2->Draw(); | |
1878 | ltx.DrawLatex(0.14, 0.76, Form("%s", centLabel(k))); | |
1879 | ||
1880 | for (int n=1; n<=4; n++) { | |
1881 | cv->cd(n); | |
1882 | ltx.SetTextSize(n==1? 0.1 : 0.15); | |
1883 | ltx.DrawLatex((n==1)? 0.45 : 0.1, 0.86, Form("V_{%d#Delta}", n)); | |
1884 | } | |
1885 | cv->cd(5); | |
1886 | ltx.DrawLatex(0.72, 0.86, Form("V_{%d#Delta}", 5)); | |
1887 | overlay->cd(); | |
1888 | ltx.SetTextSize(0.075); | |
1889 | ltx.DrawLatex(0.4, 0.05, "associated p_{T} [GeV/c]"); | |
1890 | ltx.SetTextAngle(90); | |
1891 | ltx.DrawLatex(0.03, 0.4, "V_{n#Delta}(p_{T}^{assoc.})"); | |
1892 | ltx.SetTextAngle(0); | |
1893 | ||
1894 | return cv; | |
1895 | } | |
1896 | ||
1897 | TCanvas* DrawVnVsPtTrig(int k, int npta, int ptabins[], TString opt) | |
1898 | { | |
1899 | TString title = Form("VnDeltaVsPtTrig_cent%dto%d_pTassoc%.2gto%.2g", | |
1900 | centlow[k],centhigh[k], | |
1901 | ptlow[ptabins[0]], pthigh[ptabins[npta-1]]); | |
1902 | if (!opt.IsNull()) { | |
1903 | title.Append("_"); | |
1904 | title.Append(opt); | |
1905 | } | |
1906 | TCanvas* cv = new TCanvas(title.Data(), title.Data(), 1400, 500); | |
1907 | double lv2bot = 0.9 - 0.07*npta; | |
1908 | TLegend* lv2 = new TLegend(0.82, lv2bot, 0.98, 0.95, "p_{T}^{a} (GeV/c)"); | |
1909 | lv2->SetFillColor(kNone); | |
1910 | lv2->SetBorderSize(0); | |
1911 | utils.padsetup(cv, 5, 1, "", 0.12, 0.01, 0.03, 0.2); | |
1912 | ||
1913 | // Start without knowing y-limits, adjust below. | |
1914 | double minVn_k=-1., maxVn_k=1.; | |
1915 | TObjArray* graphs = new TObjArray(); | |
1916 | ||
1917 | for (int n=1; n<=5; n++) { | |
1918 | cv->cd(n); | |
1919 | double vertMargin = 0.1*(maxVn_k - minVn_k); | |
1920 | double xmax = 11.1; | |
1921 | double xmin = -0.43; | |
1922 | ||
1923 | TH1F* hf = gPad->DrawFrame(xmin, minVn_k - 4*vertMargin, | |
1924 | xmax, maxVn_k + vertMargin); | |
1925 | TLine l; | |
1926 | l.DrawLine(xmin, 0, xmax, 0); | |
1927 | ||
1928 | int nXticks = 404; | |
1929 | int nYticks = 210; | |
1930 | if (n==1) { | |
1931 | utils.make_nice_axes(cv, hf, 2., nXticks, nYticks, -0.02, 0.02); | |
1932 | hf->SetLabelOffset(-0.01, "X"); | |
1933 | } | |
1934 | else { | |
1935 | utils.make_nice_axes(cv, hf, 3.4, nXticks, nYticks, -0.075, 0.0); | |
1936 | hf->SetLabelOffset(-0.057, "X"); | |
1937 | hf->SetTickLength(0.02, "X"); | |
1938 | hf->SetTickLength(0.05, "Y"); | |
1939 | } | |
1940 | hf->GetXaxis()->SetTicks("+-"); | |
1941 | hf->GetYaxis()->SetTicks("+-"); | |
1942 | ||
1943 | for (int jpta=0; jpta<npta; jpta++) { | |
1944 | int j = ptabins[jpta]; | |
1945 | TGraphErrors *gc = VnDeltaVsPtTrig(j, k, n); | |
1946 | if (!gc) { | |
1947 | Error("DrawVnVsPtTrig()", "Problem finding gc n%d, i%d, k%d" ,n,j,k); | |
1948 | continue; | |
1949 | } | |
1950 | TGraphErrors *gs = VnDeltaVsPtTrig(j, k, n, "sys"); | |
1951 | if (!gs) { | |
1952 | Error("DrawVnVsPtTrig()", "Problem finding gs n%d, i%d, k%d" ,n,j,k); | |
1953 | continue; | |
1954 | } | |
1955 | ||
1956 | int col = PtColor(ptlow[j], pthigh[j]); | |
1957 | utils.set_tgraph_props(gc, col, col, kFullSquare, 1.2); | |
1958 | gs->SetFillColor(PtColorPale(ptlow[j], pthigh[j])); | |
1959 | TGraphErrors* open = (TGraphErrors*)gc->Clone(); | |
1960 | utils.set_tgraph_props(open, kBlack, kBlack, kOpenSquare, 1.2); | |
1961 | ||
1962 | gs->Draw("e2psame"); | |
1963 | gc->Draw("epsame"); | |
1964 | open->Draw("epsame"); | |
1965 | ||
1966 | if (n==1) | |
1967 | lv2->AddEntry(gc,asscLabel(j), "elp"); | |
1968 | ||
1969 | graphs->Add(gc); // store graph array to find y limits for plot | |
1970 | } | |
1971 | } | |
1972 | ||
1973 | // Adjust y limits | |
1974 | double newYmin = 0, newYmax = 0, newSpace = 0; | |
1975 | for (int m=0; m<graphs->GetEntries(); m++) { | |
1976 | TGraphErrors* tg = (TGraphErrors*)graphs->At(m); | |
1977 | double yhi = TMath::MaxElement(tg->GetN()-1, tg->GetY()); | |
1978 | double ylo = TMath::MinElement(tg->GetN()-1, tg->GetY()); | |
1979 | double ey = TMath::MaxElement(tg->GetN()-1, tg->GetEY()); | |
1980 | ylo -= ey; | |
1981 | if (yhi > newYmax) newYmax = yhi; | |
1982 | if (ylo < newYmin) newYmin = ylo; | |
1983 | } | |
1984 | newSpace = 0.1*(newYmax - newYmin); | |
1985 | for (int n=1; n<=5; n++) { | |
1986 | cv->cd(n); | |
1987 | TH1F* hf = (TH1F*)gPad->FindObject("hframe"); | |
1988 | hf->GetYaxis()->SetRangeUser(newYmin-3*newSpace, newYmax+newSpace); | |
1989 | gPad->Modified(); // Signal to redraw | |
1990 | } | |
1991 | ||
1992 | cv->cd(); | |
1993 | TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1); | |
1994 | overlay->SetFillStyle(4000); // transparent | |
1995 | overlay->Draw(); | |
1996 | overlay->cd(); | |
1997 | lv2->Draw(); | |
1998 | ltx.SetTextSize(0.05); | |
1999 | ltx.DrawLatex(0.14, 0.88, Form("%s 2.76 TeV", ispp?"pp":"Pb-Pb")); | |
2000 | ltx.SetTextSize(0.07); | |
2001 | if (!ispp) | |
2002 | ltx.DrawLatex(0.14, 0.8, Form("%s", centLabel(k))); | |
2003 | ||
2004 | for (int n=1; n<=5; n++) { | |
2005 | cv->cd(n); | |
2006 | ltx.SetTextSize(n==1? 0.1 : 0.15); | |
2007 | if (n==1) | |
2008 | ltx.DrawLatex(0.44, 0.24, Form("n=%d", n)); | |
2009 | else | |
2010 | ltx.DrawLatex(0.04, 0.24, Form("n=%d", n)); | |
2011 | ||
2012 | } | |
2013 | overlay->cd(); | |
2014 | ltx.SetTextSize(0.075); | |
2015 | ltx.DrawLatex(0.48, 0.05, "trigger p_{T}^{t} [GeV/c]"); | |
2016 | ltx.SetTextAngle(90); | |
2017 | ltx.DrawLatex(0.06, 0.48, "V_{n#Delta} [10^{-2}]"); | |
2018 | ltx.SetTextAngle(0); | |
2019 | return cv; | |
2020 | } | |
2021 | ||
2022 | TGraphErrors* AgreementVsPtSum(int k, int n) | |
2023 | { | |
2024 | TGraphErrors* gVn = VnDVsQ(n, k); | |
2025 | TF1* fVn = GlobalFitVsQ(n,k); | |
2026 | TGraphErrors* g = new TGraphErrors(); | |
2027 | int N = gVn->GetN(); // And f better have the same N. | |
2028 | int i = 0, j = 0; | |
2029 | if(N<=0) { | |
2030 | Error("AgreementVsPtSum()", "Invalid number of points: %d", N); | |
2031 | gSystem->Exit(1); | |
2032 | } | |
2033 | ||
2034 | for (int q=0; q<N; q++) { | |
2035 | ijFromGlobalIndex(q,i,j); | |
2036 | double gy = gVn->GetY()[q]; | |
2037 | double fy = fVn->Eval(q+0.5); | |
2038 | double r = TMath::Abs((gy-fy)/gy); | |
2039 | g->SetPoint(q, PtBinCenterSum(i, j), r); | |
2040 | } | |
2041 | g->SetMarkerStyle(kFullCircle); | |
2042 | return g; | |
2043 | } | |
2044 | ||
2045 | TGraph2D* Agreement2D(int k, int n) | |
2046 | { | |
2047 | TGraphErrors* gVn = VnDVsQ(n, k); | |
2048 | TF1* fVn = GlobalFitVsQ(n,k); | |
2049 | TGraph2D* g = new TGraph2D(); | |
2050 | int N = gVn->GetN(); // And f better have the same N. | |
2051 | int i = 0, j = 0; | |
2052 | if(N<=0) { | |
2053 | Error("AgreementVsPtSum()", "Invalid number of points: %d", N); | |
2054 | gSystem->Exit(1); | |
2055 | } | |
2056 | ||
2057 | for (int q=0; q<N; q++) { | |
2058 | ijFromGlobalIndex(q,i,j); | |
2059 | double gy = gVn->GetY()[q]; | |
2060 | double fy = fVn->Eval(q+0.5); | |
2061 | double r = TMath::Abs((gy-fy)/fy); | |
2062 | g->SetPoint(q, ptmean[j], ptmean[i], r); | |
2063 | } | |
2064 | return g; | |
2065 | } | |
2066 | ||
2067 | TH2F* Agreement2DHist(int k, int n) | |
2068 | { | |
2069 | TGraphErrors* gVn = VnDVsQ(n, k); | |
2070 | TF1* fVn = GlobalFitVsQ(n,k); | |
2071 | ||
2072 | // Pt binning array | |
2073 | double ptarr[100]; | |
2074 | for (int i=0; i<maxpt; i++) ptarr[i] = gi->GetX()[i] - gi->GetEX()[i]; | |
2075 | ptarr[maxpt] = gi->GetX()[maxpt-1] + gi->GetEX()[maxpt-1]; | |
2076 | ||
2077 | int N = gVn->GetN(); // And f better have the same N. | |
2078 | int i = 0, j = 0; | |
2079 | if(N<=0) { | |
2080 | Error("AgreementVsPtSum()", "Invalid number of points: %d", N); | |
2081 | gSystem->Exit(1); | |
2082 | } | |
2083 | const char* name = Form("agmt_cent%dto%d_n%d", | |
2084 | centlow[k],centhigh[k], n); | |
2085 | TH2F* h = new TH2F(name, name, maxpt, ptarr, maxpt, ptarr); | |
2086 | h->SetTitle("title;p_{T}^{a};p_{T}^{t};"); | |
2087 | ||
2088 | for (int q=0; q<N; q++) { | |
2089 | ijFromGlobalIndex(q,i,j); | |
2090 | double gy = gVn->GetY()[q]; | |
2091 | double fy = fVn->Eval(q+0.5); | |
2092 | ||
2093 | // double r = fy==0? 0 : TMath::Abs(gy-fy) / fy; | |
2094 | double r = gy==0? 0 : TMath::Abs((gy-fy) / gy); | |
2095 | // double r = fy==0? 0 : TMath::Abs(gy / fy); | |
2096 | h->SetBinContent(j+1, i+1, r); | |
2097 | } | |
2098 | return h; | |
2099 | } | |
2100 | ||
2101 | TH2* PttPtaHist() | |
2102 | { | |
2103 | // Pt binning array | |
2104 | double ptarr[100]; | |
2105 | for (int i=0; i<maxpt; i++) ptarr[i] = gi->GetX()[i] - gi->GetEX()[i]; | |
2106 | ptarr[maxpt] = gi->GetX()[maxpt-1] + gi->GetEX()[maxpt-1]; | |
2107 | ||
2108 | int i = 0, j = 0; | |
2109 | const char* name = Form("PttPtaHist"); | |
2110 | TH2F* h = new TH2F(name, name, maxpt, ptarr, maxpt, ptarr); | |
2111 | h->SetTitle("title;p_{T}^{a} (GeV/c);p_{T}^{t} (GeV/c);"); | |
2112 | h->GetXaxis()->CenterTitle(); | |
2113 | h->GetYaxis()->CenterTitle(); | |
2114 | int N = maxpt * (maxpt + 1) / 2; | |
2115 | for (int q=0; q<N; q++) { | |
2116 | ijFromGlobalIndex(q,i,j); | |
2117 | h->SetBinContent(j+1, i+1, 1.); | |
2118 | } | |
2119 | return h; | |
2120 | } | |
2121 | ||
2122 | TGraphErrors* Chi2VsTrigPt(int k, int n) | |
2123 | { | |
2124 | TGraphErrors* g = new TGraphErrors(); | |
2125 | for (int i=0; i<maxpt; i++) { | |
2126 | g->SetPoint(i, ptmean[i], ReducedChi2(i, k, n)); | |
2127 | } | |
2128 | g->SetMarkerStyle(kFullCircle); | |
2129 | return g; | |
2130 | } | |
2131 | ||
2132 | TGraphErrors* | |
2133 | CorrFnChi2VsTrigPt(int j, int k, int n1, int n2, TString opt) | |
2134 | { | |
2135 | TGraphErrors* g = new TGraphErrors(); | |
2136 | int ip=0; | |
2137 | for (int i=j; i<maxpt; i++) { | |
2138 | TH1* h = Hist("RIDGE", "cA", i, j, k); | |
2139 | double x = MeanPt(i, j, k, "t"); | |
2140 | double y = Chi2(h, HarmonicSum(h, n1, n2, "global"), opt); | |
2141 | g->SetPoint(ip, x, y); | |
2142 | ip++; | |
2143 | } | |
2144 | g->SetMarkerStyle(kFullCircle); | |
2145 | return g; | |
2146 | } | |
2147 | ||
2148 | double ReducedChi2(int i, int k, int n) | |
2149 | { | |
2150 | TGraphErrors* g = VnDeltaVsPtAssoc(i, k, n); | |
2151 | TGraphErrors* f = GlobalFitVsPtAssoc(i, k, n); | |
2152 | int N = g->GetN(); // And f better have the same N. | |
2153 | double chi2 = 0; | |
2154 | ||
2155 | // Value (or one-sigma error if opt=="err") of v_n | |
2156 | // double vnerr = Bestvn(k,n,i,"err"); | |
2157 | double vnerr = vnGF(n,k,i,0,999,"RIDGE","cA","err"); | |
2158 | ||
2159 | if(N<=0 || N!= f->GetN()) { | |
2160 | Error("ReducedChi2()", "Invalid number of points: %d", N); | |
2161 | gSystem->Exit(1); | |
2162 | } | |
2163 | ||
2164 | for (int j=0; j<N; j++) { | |
2165 | double errj = g->GetEY()[j]; | |
2166 | ||
2167 | chi2 += | |
2168 | (g->GetY()[j] - f->GetY()[j]) * (g->GetY()[j] - f->GetY()[j]) / | |
2169 | pow(errj + vnerr, 2); | |
2170 | } | |
2171 | chi2 /= N; | |
2172 | return chi2; | |
2173 | } | |
2174 | ||
2175 | TGraphErrors* VnDeltaVsPtAssoc(int i, int k, int n, TString opt) | |
2176 | { | |
2177 | // Create a TGraph from VnDelta points | |
2178 | TGraphErrors* gALL = VnDVsQ(n, k, "RIDGE", "cA", opt); | |
2179 | TGraphErrors* g = new TGraphErrors(); | |
2180 | TString gname(Form("VnDeltaVsPtAssoc_ptt%d_cent%d_n%d",i,k,n)); | |
2181 | gname.Append(opt); | |
2182 | g->SetName(gname.Data()); | |
2183 | ||
2184 | for (int j=0; j<=i; j++) { | |
2185 | int q = GlobalIndex(i, j); | |
2186 | g->SetPoint(j, gj->GetX()[j], gALL->GetY()[q] ); | |
2187 | g->SetPointError(j, 0, gALL->GetEY()[q] ); | |
2188 | utils.set_tgraph_props(g, colorsDark[i], colorsDark[i], kFullSquare, 1.0); | |
2189 | g->SetLineWidth(2); | |
2190 | } | |
2191 | ||
2192 | return g; | |
2193 | } | |
2194 | ||
2195 | TGraphErrors* VnDeltaVsPtTrig(int j, int k, int n, TString opt) | |
2196 | { | |
2197 | // Create a TGraph from VnDelta points | |
2198 | TGraphErrors* gALL = VnDVsQ(n, k, "RIDGE", "cA", opt); | |
2199 | TGraphErrors* g = new TGraphErrors(); | |
2200 | TString gname(Form("VnDeltaVsPtTrig_pta%d_cent%dto%d_n%d", | |
2201 | j,centlow[k],centhigh[k],n)); | |
2202 | gname.Append(opt); | |
2203 | g->SetName(gname.Data()); | |
2204 | ||
2205 | int ip=0; | |
2206 | for (int i=j; i<maxpt; i++) { | |
2207 | int q = GlobalIndex(i, j); | |
2208 | ||
2209 | double xerr = 0.; | |
2210 | if (opt.Contains("sys")) { | |
2211 | if (gi->GetEX()[i] < 0.126) | |
2212 | xerr = gi->GetEX()[i]; | |
2213 | else | |
2214 | xerr = 0.250; | |
2215 | } | |
2216 | double mpt = MeanPt(i, j, k, "t"); | |
2217 | double sf = 100; // Scale y-value; make sure to label accordingly!! | |
2218 | g->SetPoint(ip, mpt, sf*gALL->GetY()[q] ); | |
2219 | g->SetPointError(ip, xerr, sf*gALL->GetEY()[q] ); | |
2220 | ||
2221 | utils.set_tgraph_props(g, colorsDark[j], colorsDark[j], kFullSquare, 1.0); | |
2222 | g->SetLineWidth(2); | |
2223 | ip++; | |
2224 | } | |
2225 | ||
2226 | return g; | |
2227 | } | |
2228 | ||
2229 | TGraphErrors* VnDeltaVsN(int i, int j, int k, int nMax, TString opt) | |
2230 | { | |
2231 | // Create a TGraph from VnDelta points | |
2232 | // Stat errors by default. Call with opt "sys" to get sys errors. | |
2233 | // Call with opt "sine" to get residual <sin n Delta phi> for sys. est. | |
2234 | int q = GlobalIndex(i, j); | |
2235 | TString gname(Form("VnDeltaVsN_ptt%.2gto%.2g_pta%.2gto%.2g_cent%dto%d", | |
2236 | ptlow[i],pthigh[i],ptlow[j],pthigh[j], | |
2237 | centlow[k],centhigh[k])); | |
2238 | if (!opt.IsNull()) { | |
2239 | gname.Append("_"); | |
2240 | gname.Append(opt); | |
2241 | } | |
2242 | TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data()); | |
2243 | if (g) { | |
2244 | return g; | |
2245 | } | |
2246 | else { | |
2247 | g = new TGraphErrors(); | |
2248 | g->SetName(gname.Data()); | |
2249 | } | |
2250 | ||
2251 | int ip=0; | |
2252 | gStyle->SetEndErrorSize(5); | |
2253 | for (int n=1; n<=nMax; n++) { | |
2254 | ||
2255 | TGraphErrors* gALL = VnDVsQ(n, k, "RIDGE", "cA", opt); | |
2256 | double xerr = opt.Contains("sys") ? 0. : 0.0; // 0.2 | |
2257 | g->SetPoint(ip, n, 100*gALL->GetY()[q] ); | |
2258 | g->SetPointError(ip, xerr, 100*gALL->GetEY()[q] ); | |
2259 | if (0) | |
2260 | cout<<gALL->GetName()<< " " << g->GetY()[ip] << endl; | |
2261 | ||
2262 | ip++; | |
2263 | } | |
2264 | ||
2265 | utils.set_tgraph_props(g, kBlack, kBlack, kFullCircle, 1.2); | |
2266 | g->SetLineWidth(2); | |
2267 | gList->Add(g); | |
2268 | return g; | |
2269 | } | |
2270 | ||
2271 | TGraphErrors* GlobalvnVsN(int i, int k, int nMax, TString opt) | |
2272 | { | |
2273 | // Create a TGraph from vn{GF} points | |
2274 | // Stat errors by default. Call with opt "sys" to get sys errors. | |
2275 | ||
2276 | TString gname(Form("GlobalvnVsN_pt%.2gto%.2g_cent%dto%d", | |
2277 | ptlow[i],pthigh[i], centlow[k],centhigh[k])); | |
2278 | if (!opt.IsNull()) { | |
2279 | gname.Append("_"); | |
2280 | gname.Append(opt); | |
2281 | } | |
2282 | TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data()); | |
2283 | if (g) { | |
2284 | return g; | |
2285 | } | |
2286 | else { | |
2287 | g = new TGraphErrors(); | |
2288 | g->SetName(gname.Data()); | |
2289 | } | |
2290 | g->SetName(gname.Data()); | |
2291 | ||
2292 | int ip=0; | |
2293 | for (int n=1; n<=nMax; n++) { | |
2294 | ||
2295 | // point, stat err, sys err. | |
2296 | double vnp = vnGF(n, k, i, 0, 999, "RIDGE", "cA", ""); | |
2297 | double vne = vnGF(n, k, i, 0, 999, "RIDGE", "cA", "err"); | |
2298 | double vns = vnGF(n, k, i, 0, 999, "RIDGE", "cA", "sys"); | |
2299 | double xerr = opt.Contains("sys") ? 0. : 0.0; // 0.2 | |
2300 | double yerr = opt.Contains("sys") ? vns : vne; | |
2301 | ||
2302 | g->SetPoint(ip, n, vnp ); | |
2303 | g->SetPointError(ip, xerr, yerr); | |
2304 | ip++; | |
2305 | } | |
2306 | ||
2307 | utils.set_tgraph_props(g, kBlack, kBlack, kFullCircle, 1.2); | |
2308 | g->SetLineWidth(2); | |
2309 | ||
2310 | return g; | |
2311 | } | |
2312 | ||
2313 | TGraphErrors* VnDeltaNFVsN(int i, int j, int k, int nMax, TString opt) | |
2314 | { | |
2315 | ||
2316 | TGraphErrors* VnD = VnDeltaVsN(i,j,k,nMax,opt); | |
2317 | TGraphErrors* vnt = GlobalvnVsN(i,k,nMax,opt); | |
2318 | TGraphErrors* vna = GlobalvnVsN(j,k,nMax,opt); | |
2319 | TString gname(Form("VnDeltaNFVsN_ptt%.2gto%.2g_pta%.2gto%.2g_cent%dto%d", | |
2320 | ptlow[i],pthigh[i],ptlow[j],pthigh[j], | |
2321 | centlow[k],centhigh[k])); | |
2322 | if (!opt.IsNull()) { | |
2323 | gname.Append("_"); | |
2324 | gname.Append(opt); | |
2325 | } | |
2326 | TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data()); | |
2327 | if (g) { | |
2328 | return g; | |
2329 | } | |
2330 | else { | |
2331 | g = new TGraphErrors(); | |
2332 | g->SetName(gname.Data()); | |
2333 | } | |
2334 | g->SetName(gname.Data()); | |
2335 | ||
2336 | int ip=0; | |
2337 | for (int n=1; n<=nMax; n++) { | |
2338 | ||
2339 | // Uncertainty vne is stat or sys depending on opt. | |
2340 | double vnp = VnD->GetY()[ip] - 100*vnt->GetY()[ip]*vna->GetY()[ip]; | |
2341 | double vne = VnD->GetEY()[ip] - 100*vnt->GetEY()[ip]*vna->GetEY()[ip]; | |
2342 | g->SetPoint(ip, n, vnp ); | |
2343 | g->SetPointError(ip, 0.0, vne); | |
2344 | ip++; | |
2345 | } | |
2346 | ||
2347 | utils.set_tgraph_props(g, kBlack, kBlack, kOpenSquare, 1.2); | |
2348 | g->SetLineWidth(2); | |
2349 | return g; | |
2350 | } | |
2351 | ||
2352 | ||
2353 | TF1* GlobalFitVsQ(int n, int k, int ipt1, int ipt2, const char* region, | |
2354 | const char* corrtype, TString opt) | |
2355 | { | |
2356 | ||
2357 | // vn{GF}(trig) x vn{GF}(assc) vs. q, where the points between (and | |
2358 | // including) ipt1 and ipt2 were used in the global fit. | |
2359 | // Note: gfname must be identical format to fnName in GlobalFitVsQ(). | |
2360 | TF1* fitfn = 0; | |
2361 | TString errType("meas"); | |
2362 | if (opt.Contains("hi_sys")) | |
2363 | errType = "hi"; | |
2364 | if (opt.Contains("lo_sys")) | |
2365 | errType = "lo"; | |
2366 | TString gfname = | |
2367 | Form("globalFitFn_cent%dto%d_n%d_ptbins%dto%d_%s_corrtype_%s_%s", | |
2368 | centlow[k],centhigh[k], n, ipt1, ipt2, | |
2369 | region, corrtype, errType.Data()); | |
2370 | fitfn = (TF1*)gFitList->FindObject(gfname.Data()); | |
2371 | if (fitfn) { | |
2372 | if (0) | |
2373 | Info("GlobalFitVsQ()", "Returning cached TF1 %s", gfname.Data()); | |
2374 | return fitfn; | |
2375 | } | |
2376 | ||
2377 | // If nonexistent, create (& store) a new one | |
2378 | fitfn = DoGlobalFit(n,k,ipt1,ipt2,region,corrtype,opt); | |
2379 | // fitfn = (TF1*)gFitList->FindObject(gfname.Data()); | |
2380 | ||
2381 | if (!fitfn) | |
2382 | Warning("GlobalFitVsQ()", | |
2383 | "no fitfn %s,\neven after calling DoGlobalFit()", | |
2384 | gfname.Data()); | |
2385 | ||
2386 | return fitfn; | |
2387 | } | |
2388 | ||
2389 | TF1* DoGlobalFit(int n, int k, int ipt1, int ipt2, | |
2390 | const char* region, const char* corrtype, TString opt) | |
2391 | { | |
2392 | // fnName must be identical format to gfname in GlobalFitVsQ(). | |
2393 | TF1 *globalFitFn = 0; | |
2394 | TGraphErrors* g = 0; | |
2395 | TString errType("meas"); | |
2396 | if (opt.Contains("hi_sys")) | |
2397 | errType = "hi"; | |
2398 | if (opt.Contains("lo_sys")) | |
2399 | errType = "lo"; | |
2400 | TString fnName = | |
2401 | Form("globalFitFn_cent%dto%d_n%d_ptbins%dto%d_%s_corrtype_%s_%s", | |
2402 | centlow[k],centhigh[k], n, ipt1, ipt2, | |
2403 | region, corrtype, errType.Data()); | |
2404 | ||
2405 | g = VnDVsQ(n, k, region, corrtype, opt); | |
2406 | ||
2407 | if (!g) { | |
2408 | Error("DoGlobalFit()", "No graph found"); | |
2409 | return 0; | |
2410 | } | |
2411 | ||
2412 | // Assign global variables that will get checked in the vntvna | |
2413 | // function during the fit. | |
2414 | minFitPtBin = ipt1; | |
2415 | maxFitPtBin = ipt2; | |
2416 | ||
2417 | // Set up the global fit function... | |
2418 | double qmax = g->GetN() - 0.4999; | |
2419 | int nPars = maxpt + 3; | |
2420 | globalFitFn = new TF1(fnName.Data(), vntvna, 0.0, qmax, nPars); | |
2421 | ||
2422 | /* Not needed now, I just flip the v1 points later...AMA 11/15/2011 */ | |
2423 | // Each parameter has an equally good minimum at + or - v_nbest. The | |
2424 | // starting value is chosen to steer the fit to the positive "root". | |
2425 | // v_1 is special because it is the only coefficient that crosses | |
2426 | // zero at low pt, so it gets a negative starting point. Pars are | |
2427 | // vn(ipt), then last one is n. | |
2428 | // double startVal = (n==1)? -0.02 : 0.01; | |
2429 | ||
2430 | double startVal = 0.01; | |
2431 | for (int ip=0; ip<maxpt; ip++) { | |
2432 | globalFitFn->SetParameter(ip, startVal); | |
2433 | globalFitFn->SetParName(ip, Form("v_%d(%.3g)",n,ptmean[ip])); | |
2434 | } | |
2435 | ||
2436 | // Pass n into vntvna(). Access as par[maxpt]. | |
2437 | globalFitFn->FixParameter(maxpt, double(n)); | |
2438 | ||
2439 | // Pass k into vntvna(). Access as par[maxpt+1]. | |
2440 | globalFitFn->FixParameter(maxpt+1, double(k)); | |
2441 | ||
2442 | // Mom. conservation term prefactor. Access as par[maxpt+2]. | |
2443 | if (n==1) { | |
2444 | globalFitFn->FixParameter(maxpt+2, 0.); | |
2445 | // globalFitFn->SetParLimits(maxpt+2, 0., 1e-3); | |
2446 | } | |
2447 | else | |
2448 | globalFitFn->FixParameter(maxpt+2, 0.); | |
2449 | ||
2450 | globalFitFn->SetNpx(1000); | |
2451 | ||
2452 | TFitResultPtr result = g->Fit(globalFitFn,"QRNS"); // fast, takes ~15 ms | |
2453 | ||
2454 | if (0) { | |
2455 | cout<<globalFitFn->GetName(); | |
2456 | result->Print(); | |
2457 | } | |
2458 | // TMinuit status: | |
2459 | // migradResult + 10*minosResult + 100*hesseResult + 1000*improveResult | |
2460 | // 0:ok, <0:user error, >0: see above. | |
2461 | int fitStatus = result; | |
2462 | if (fitStatus) | |
2463 | Warning("DoGlobalFit()", "TMinuit status %d", fitStatus); | |
2464 | ||
2465 | gFitList->Add(globalFitFn); | |
2466 | ||
2467 | return globalFitFn; | |
2468 | } | |
2469 | ||
2470 | double vnGF(int n, int k, int ptBin, int ipt1, int ipt2, | |
2471 | const char* region, const char* corrtype, TString opt) | |
2472 | { | |
2473 | ||
2474 | double vn, err; | |
2475 | TF1* gfn = GlobalFitVsQ(n,k,ipt1,ipt2,region,corrtype,opt); | |
2476 | if (!gfn) | |
2477 | Error("vnGF()", "!gfn"); | |
2478 | vn = gfn->GetParameter(ptBin); // ptBin is the q index | |
2479 | err = gfn->GetParError(ptBin); | |
2480 | ||
2481 | if (opt.Contains("sys")) { | |
2482 | TF1* hi = GlobalFitVsQ(n,k,ipt1,ipt2,region,corrtype,"hi_sys"); | |
2483 | TF1* lo = GlobalFitVsQ(n,k,ipt1,ipt2,region,corrtype,"lo_sys"); | |
2484 | double vhi = hi->GetParameter(ptBin); | |
2485 | double vlo = lo->GetParameter(ptBin); | |
2486 | return 0.5*TMath::Abs(vhi-vlo); | |
2487 | } | |
2488 | ||
2489 | // v1 can be negative, otherwise ensure the positive solution is chosen | |
2490 | // if (n>1 && vn < 0) | |
2491 | // vn *= -1; | |
2492 | ||
2493 | // if (ptlow[ipt1] > 4.) // !!!!!!!!!!!!!!!!!!!!!!!!! | |
2494 | // vn = TMath::Abs(vn); | |
2495 | ||
2496 | // if (n==1) | |
2497 | // vn += MomCorrection(i, j, k); | |
2498 | ||
2499 | return opt.Contains("err") ? err : vn; | |
2500 | } | |
2501 | ||
2502 | TGraphErrors* vnGFvsPt(int n, int k, int ipt1, int ipt2, | |
2503 | const char* region, const char* corrtype, TString opt) | |
2504 | { | |
2505 | ||
2506 | // vn{GF} vs. pt, where the points between (and including) ipt1 and | |
2507 | // ipt2 were used in the global fit. | |
2508 | bool isHighPtFit = ipt1 > 0 && ipt2==999; | |
2509 | bool flipSign = n==1; // Flip v1? | |
2510 | if (isHighPtFit) | |
2511 | flipSign = 0; | |
2512 | if (opt.Contains("keepsign")) | |
2513 | flipSign = 0; | |
2514 | TString name(Form("v%dGFvsPt_fitpt%dto%d_cent%dto%d", | |
2515 | n, ipt1, ipt2, centlow[k], centhigh[k])); | |
2516 | if (opt.Contains("sys")) | |
2517 | name.Append("_sys"); | |
2518 | if (opt.Contains("ptcons")) | |
2519 | name.Append("_ptcons"); | |
2520 | ||
2521 | TGraphErrors* g = (TGraphErrors*)gList->FindObject(name.Data()); | |
2522 | if (g) | |
2523 | return g; | |
2524 | ||
2525 | g = new TGraphErrors(); | |
2526 | g->SetName(name); | |
2527 | ||
2528 | for (int i=0; i<maxpt; i++) { | |
2529 | TString opt1 = ""; | |
2530 | if (opt.Contains("ptcons")) | |
2531 | opt1.Append("_ptcons"); | |
2532 | double vn = vnGF(n,k,i,ipt1,ipt2,region,corrtype,opt1); | |
2533 | ||
2534 | if (flipSign) | |
2535 | vn = -vn; | |
2536 | ||
2537 | // if (isHighPtFit && i<ipt1) | |
2538 | if (i<ipt1) | |
2539 | vn = 0.; | |
2540 | ||
2541 | double mpt = MeanPt(i, 0, k, "t"); | |
2542 | if (i<ipt1) | |
2543 | mpt = -10.; | |
2544 | g->SetPoint(i, mpt, vn); | |
2545 | ||
2546 | if (opt.Contains("sys")) { | |
2547 | opt1.Append("_sys"); | |
2548 | double ex = (gi->GetEX()[i] < 0.126) ? gi->GetEX()[i] : 0.250; | |
2549 | double vns = vnGF(n,k,i,ipt1,ipt2,region,corrtype, opt1); | |
2550 | g->SetPointError(i, ex, vns); | |
2551 | // g->SetPointError(i, gi->GetEX()[i], vns); // too wide | |
2552 | } | |
2553 | else { | |
2554 | opt1.Append("_err"); | |
2555 | double vne = vnGF(n,k,i,ipt1,ipt2,region,corrtype,opt1); | |
2556 | g->SetPointError(i, 0.0, vne); | |
2557 | } | |
2558 | } | |
2559 | ||
2560 | int col = CentColor(centlow[k], centhigh[k]); | |
2561 | int col_pale = CentColorPale(centlow[k], centhigh[k]); | |
2562 | g->SetLineWidth(2); | |
2563 | g->SetFillColor(col_pale); // for sys | |
2564 | utils.set_tgraph_props(g, col, col, kFullCircle, 1.2); | |
2565 | ||
2566 | if (0) cout<<opt.Data()<<endl; | |
2567 | gList->Add(g); | |
2568 | return g; | |
2569 | } | |
2570 | ||
2571 | TGraph* Luzumv1(int cent, const char* pid) | |
2572 | { | |
2573 | // v1 for pions vs pt from viscous hydro calc by Matt Luzum | |
2574 | // valid cent args 0-4 - 10% bins. | |
2575 | // pid = "pi", "K", "p". | |
2576 | const char* dat = Form("v1LHCGlauber%ssv08.dat", pid); | |
2577 | TString fmt("%lg "); // pt | |
2578 | for (int ky=0; ky<5; ky++) | |
2579 | fmt.Append( ky==cent? "%lg " : "%*lg "); // centrality - skip if *lg | |
2580 | return new TGraph(dat, fmt.Data()); | |
2581 | } | |
2582 | ||
2583 | TGraphErrors* GlobalFitVsPtAssoc(int i, int k, int n, TString opt) | |
2584 | { | |
2585 | // Create a TGraph of vnt{gf} x vna{gf} vs. pt_assoc from gfit | |
2586 | TF1* gfit = GlobalFitVsQ(n, k, 0, 999, "RIDGE", cfDef, opt); | |
2587 | TGraphErrors* g = new TGraphErrors(); | |
2588 | TString gname(Form("GlobalFitVsPtAssoc_ptt%d_cent%d_n%d",i,k,n)); | |
2589 | if (!opt.IsNull()) | |
2590 | gname.Append(opt); | |
2591 | g->SetName(gname); | |
2592 | ||
2593 | for (int j=0; j<=i; j++) { | |
2594 | int q = GlobalIndex(i, j); | |
2595 | g->SetPoint(j, gj->GetX()[j], gfit->Eval(q + 0.5) ); | |
2596 | g->SetPointError(j, 0, 0 ); | |
2597 | utils.set_tgraph_props(g, colorsPale[i], colorsPale[i], kFullCircle, 1.0); | |
2598 | g->SetLineWidth(4); | |
2599 | } | |
2600 | return g; | |
2601 | } | |
2602 | ||
2603 | double SineError(int i, int j, int k) | |
2604 | { | |
2605 | double rms = 0; | |
2606 | double resid_sin=0, rs_err=0; | |
2607 | TH1* hst = Hist("RIDGE", "cA", i, j, k, ""); | |
2608 | double arr[100] = {0}; | |
2609 | // arr[0] = 0; | |
2610 | for (int n=1; n<=gNMax; n++) { | |
2611 | VnDelta(n, *hst, resid_sin, rs_err, "sin"); | |
2612 | arr[n-1] = resid_sin; | |
2613 | } | |
2614 | ||
2615 | rms = TMath::RMS(gNMax, arr) / 2; | |
2616 | if (0) | |
2617 | cout << " SineError: " << rms << endl; | |
2618 | ||
2619 | return rms; | |
2620 | } | |
2621 | ||
2622 | TCanvas* DrawVnDeltaVsN(int nptt, int pttbins[], int npta, | |
2623 | int ptabins[], int ncb, int centbins[], TString opt) | |
2624 | { | |
2625 | ||
2626 | double t1=ptlow[pttbins[0]], t2=pthigh[pttbins[nptt-1]]; | |
2627 | double a1=ptlow[ptabins[0]], a2=pthigh[ptabins[npta-1]]; | |
2628 | TString cname = Form("VnDeltaVsN_ptt%.2gto%.2g_pta%.2gto%.2g",t1,t2,a1,a2); | |
2629 | if (!opt.IsNull()) { | |
2630 | cname.Append("_"); | |
2631 | cname.Append(opt); | |
2632 | } | |
2633 | TCanvas* c = new TCanvas(cname.Data(), cname.Data(), 650,700); | |
2634 | double legbottom = 0.9 - 0.07*ncb; | |
2635 | TLegend* leg = new TLegend(0.62, legbottom, 0.98, 0.98, "Centrality"); | |
2636 | leg->SetFillColor(kNone); | |
2637 | leg->SetBorderSize(0); | |
2638 | leg->SetName("leg"); | |
2639 | ||
2640 | TObjArray* graphs = new TObjArray(); | |
2641 | double xmin = 0.4; | |
2642 | double xmax = opt.Contains("ext") ? gNMax + 0.5 : 8.5; | |
2643 | int nXticks = 10, nYticks = 210; | |
2644 | TLine l; | |
2645 | TH1F* hf = gPad->DrawFrame(xmin, 0.0001, xmax, 0.26); | |
2646 | // bool scale1000 = 0; | |
2647 | ||
2648 | hf->SetTitle(Form(";n;V_{n#Delta} [10^{-2}]")); | |
2649 | if (opt.Contains("sine")) | |
2650 | hf->SetTitle(Form(";n;#LTsin(n#Delta#phi)#GT [10^{-2}]")); | |
2651 | ||
2652 | // if (opt.Contains("ext") ) | |
2653 | // if (!opt.Contains("sine")) { | |
2654 | // scale1000 = 0; | |
2655 | // hf->SetTitle(Form(";n;V_{n#Delta}")); | |
2656 | // } | |
2657 | ||
2658 | l.SetLineStyle(kDashed); | |
2659 | ||
2660 | utils.make_nice_axes(c, hf, 1.3, nXticks, nYticks, 1.4,-1.4); | |
2661 | hf->SetLabelOffset(0.0, "X"); | |
2662 | hf->SetTitleOffset(-1.96, "Y"); | |
2663 | ||
2664 | hf->GetXaxis()->SetTicks("+-"); | |
2665 | hf->GetYaxis()->SetTicks("+-"); | |
2666 | l.DrawLine(hf->GetXaxis()->GetXmin(), 0, hf->GetXaxis()->GetXmax(),0); | |
2667 | gPad->SetLeftMargin(0.22); | |
2668 | gPad->SetTopMargin(0.01); | |
2669 | gPad->SetRightMargin(0.01); | |
2670 | ||
2671 | for (int i=0; i<nptt; i++) { | |
2672 | for (int j=0; j<npta; j++) { | |
2673 | for(int k=0; k<ncb; k++) { | |
2674 | ||
2675 | // Stat errors on g, sys errors on gs. | |
2676 | TString opt1 = opt.Contains("sine") ? "sine" : ""; | |
2677 | TGraphErrors* g = | |
2678 | VnDeltaVsN(pttbins[i], ptabins[j], centbins[k], gNMax, opt1); | |
2679 | TGraphErrors* gs = | |
2680 | VnDeltaVsN(pttbins[i], ptabins[j], centbins[k], gNMax, "sys"); | |
2681 | TGraphErrors* gc = (TGraphErrors*)g->Clone(); | |
2682 | int dark = CentColor(centlow[centbins[k]], centhigh[centbins[k]]); | |
2683 | int pale = CentColorPale(centlow[centbins[k]], centhigh[centbins[k]]); | |
2684 | utils.set_tgraph_props(g, dark, dark, kFullCircle, 1.5); | |
2685 | utils.set_tgraph_props(gs, dark, dark, kDot, 2.); | |
2686 | utils.set_tgraph_props(gc, kBlack, kBlack, kOpenCircle, 1.5); | |
2687 | g->SetFillColor(pale); | |
2688 | gs->SetFillColor(dark); | |
2689 | ||
2690 | TGraphErrors* gnf = 0; | |
2691 | TGraphErrors* gnf_sys = 0; | |
2692 | TGraphErrors* gnf_solid = 0; | |
2693 | if (opt.Contains("vnf")) { // Non-factorizing part | |
2694 | gnf = | |
2695 | VnDeltaNFVsN(pttbins[i], ptabins[j], centbins[k], gNMax, ""); | |
2696 | gnf_sys = | |
2697 | VnDeltaNFVsN(pttbins[i], ptabins[j], centbins[k], gNMax, "sys"); | |
2698 | ||
2699 | // Just for appearance... | |
2700 | gnf_solid = (TGraphErrors*) gnf->Clone(Form("%s_solid", gnf->GetName())); | |
2701 | utils.set_tgraph_props(gnf_solid, pale, pale, kFullSquare, gnf->GetMarkerSize()); | |
2702 | } | |
2703 | // if (scale1000) { | |
2704 | // for (int ip=0; ip<g->GetN(); ip++) { | |
2705 | // g->SetPoint(ip, g->GetX()[ip], 1000* g->GetY()[ip]); | |
2706 | // gs->SetPoint(ip, gs->GetX()[ip], 1000*gs->GetY()[ip]); | |
2707 | // gc->SetPoint(ip, gc->GetX()[ip], 1000*gc->GetY()[ip]); | |
2708 | ||
2709 | // g->SetPointError(ip, g->GetEX()[ip], 1000* g->GetEY()[ip]); | |
2710 | // gs->SetPointError(ip, gs->GetEX()[ip], 1000*gs->GetEY()[ip]); | |
2711 | // gc->SetPointError(ip, gc->GetEX()[ip], 1000*gc->GetEY()[ip]); | |
2712 | // } | |
2713 | // } | |
2714 | ||
2715 | if (opt.Contains("sine")) { | |
2716 | g->Fit("pol0", "Q"); | |
2717 | TF1* gfit = g->GetFunction("pol0"); | |
2718 | gfit->SetLineColor(dark); | |
2719 | } | |
2720 | ||
2721 | else if (!opt.Contains("nobars")) | |
2722 | g->Draw("Bsame"); // bars | |
2723 | ||
2724 | if (opt.Contains("vnf")) { // Non-factorizing part | |
2725 | utils.draw_errorbox(gnf_sys, pale, 0.3); | |
2726 | gnf_solid->Draw("epsame"); | |
2727 | gnf->Draw("epsame"); | |
2728 | } | |
2729 | ||
2730 | g->Draw("psame"); // colored markers | |
2731 | gc->Draw("epsame"); // open markers w/black stat error bars | |
2732 | ||
2733 | if (!opt.Contains("sine")) | |
2734 | utils.draw_errorbox(gs, dark, 0.3); | |
2735 | ||
2736 | graphs->Add(g); | |
2737 | leg->AddEntry(g,Form("%d-%d%%",centlow[centbins[k]],centhigh[centbins[k]]),"l,p"); | |
2738 | ||
2739 | if (opt.Contains("vnf")) { // Non-factorizing part | |
2740 | leg->AddEntry(gnf_solid,"NF V_{n#Delta}","l,p"); | |
2741 | // leg->AddEntry(gnf_solid,"[v_{n}(p_{T}^{t}) #times v_{n}(p_{T}^{a})]_{GF}","l,p"); | |
2742 | leg->SetHeader(""); | |
2743 | } | |
2744 | ||
2745 | } | |
2746 | } | |
2747 | } | |
2748 | ||
2749 | // Adjust y limits | |
2750 | double newYmin = 0, newYmax = 0, newSpace = 0; | |
2751 | for (int m=0; m<graphs->GetEntries(); m++) { | |
2752 | TGraphErrors* tg = (TGraphErrors*)graphs->At(m); | |
2753 | double yhi = TMath::MaxElement(tg->GetN()-1, tg->GetY()); | |
2754 | double ylo = TMath::MinElement(tg->GetN()-1, tg->GetY()); | |
2755 | if (yhi > newYmax) newYmax = yhi; | |
2756 | if (ylo < newYmin) newYmin = ylo; | |
2757 | } | |
2758 | ||
2759 | newSpace = 0.1*(newYmax - newYmin); | |
2760 | ||
2761 | if (opt.Contains("ext")) | |
2762 | hf->GetYaxis()->SetRangeUser(newYmin-newSpace, newYmax+4*newSpace); | |
2763 | else | |
2764 | hf->GetYaxis()->SetRangeUser(newYmin-newSpace, newYmax+newSpace); | |
2765 | ||
2766 | // if (opt.Contains("sine")) { | |
2767 | // hf->GetYaxis()->SetRangeUser(0.0001, newYmax+4*newSpace); | |
2768 | // gPad->SetLogy(); | |
2769 | // } | |
2770 | ||
2771 | if (!opt.Contains("sine")) | |
2772 | leg->Draw(); | |
2773 | ||
2774 | // Draw text | |
2775 | ltx.SetTextSize(0.04); // 0.05 | |
2776 | ||
2777 | if ( opt.Contains("sine") ) | |
2778 | legbottom = 0.99; | |
2779 | else if (opt.Contains("ext") && !opt.Contains("vnf") ) | |
2780 | legbottom -= 0.4; | |
2781 | ||
2782 | double xstart = 0.6; | |
2783 | ltx.DrawLatex(xstart, legbottom-1*0.07, Form("%.2g < p_{T}^{t} < %.2g GeV/c", t1, t2)); | |
2784 | ltx.DrawLatex(xstart, legbottom-2*0.07, Form("%.2g < p_{T}^{a} < %.2g GeV/c", a1, a2)); | |
2785 | ||
2786 | gPad->Update(); | |
2787 | gPad->Modified(); | |
2788 | return c; | |
2789 | } | |
2790 | ||
2791 | TCanvas* DrawGlobalvnVsN(int npt, int ptbins[], int ncb, int centbins[], TString opt) | |
2792 | { | |
2793 | double t1=ptlow[ptbins[0]], t2=pthigh[ptbins[npt-1]]; | |
2794 | TString cname = Form("GlobalvnVsN_pt%.2gto%.2g",t1,t2); | |
2795 | if (!opt.IsNull()) { | |
2796 | cname.Append("_"); | |
2797 | cname.Append(opt); | |
2798 | } | |
2799 | TCanvas* c = new TCanvas(cname.Data(), cname.Data(), 600,700); | |
2800 | double legbottom = 0.9 - 0.07*ncb; | |
2801 | TLegend* leg = new TLegend(0.62, legbottom, 0.98, 0.98, "Centrality"); | |
2802 | leg->SetFillColor(kNone); | |
2803 | leg->SetBorderSize(0); | |
2804 | leg->SetName("leg"); | |
2805 | ||
2806 | TObjArray* graphs = new TObjArray(); | |
2807 | double xmin = opt.Contains("zoom") ? 2.6 : 0.4; | |
2808 | double xmax = 8.5; | |
2809 | int nXticks = 10, nYticks = 210; | |
2810 | TLine l; | |
2811 | TH1F* hf = gPad->DrawFrame(xmin, -0.06, xmax, 0.26); | |
2812 | hf->SetTitle(Form(";n;v_{n}{GF}")); | |
2813 | l.SetLineStyle(kDashed); | |
2814 | ||
2815 | utils.make_nice_axes(c, hf, 1.3, nXticks, nYticks, 1.4,-1.4); | |
2816 | hf->SetLabelOffset(0.0, "X"); | |
2817 | if (opt.Contains("zoom")) | |
2818 | hf->SetTitleOffset(-2.3, "Y"); | |
2819 | else | |
2820 | hf->SetTitleOffset(-1.9, "Y"); | |
2821 | ||
2822 | hf->GetXaxis()->SetTicks("+-"); | |
2823 | hf->GetYaxis()->SetTicks("+-"); | |
2824 | l.DrawLine(hf->GetXaxis()->GetXmin(), 0, hf->GetXaxis()->GetXmax(),0); | |
2825 | if (opt.Contains("zoom")) | |
2826 | gPad->SetLeftMargin(0.25); | |
2827 | else | |
2828 | gPad->SetLeftMargin(0.2); | |
2829 | gPad->SetTopMargin(0.01); | |
2830 | gPad->SetRightMargin(0.01); | |
2831 | ||
2832 | for (int i=0; i<npt; i++) { | |
2833 | for(int k=0; k<ncb; k++) { | |
2834 | ||
2835 | // Stat errors on g, sys errors on gs. | |
2836 | TGraphErrors* g = GlobalvnVsN(ptbins[i], centbins[k], gNMax, ""); | |
2837 | TGraphErrors* gs = GlobalvnVsN(ptbins[i], centbins[k], gNMax, "sys"); | |
2838 | TGraphErrors* gc = (TGraphErrors*)g->Clone(); | |
2839 | int dark = CentColor(centlow[centbins[k]], centhigh[centbins[k]]); | |
2840 | int pale = CentColorPale(centlow[centbins[k]], centhigh[centbins[k]]); | |
2841 | utils.set_tgraph_props(g, dark, dark, kFullCircle, 1.5); | |
2842 | // utils.set_tgraph_props(gs, dark, dark, kDot, 2.); | |
2843 | utils.set_tgraph_props(gc, kBlack, kBlack, kOpenCircle, 1.5); | |
2844 | g->SetFillColor(pale); | |
2845 | // gs->SetFillColor(dark); | |
2846 | g->Draw("Bsame"); // colored bars | |
2847 | utils.draw_errorbox(gs, dark, 0.3); | |
2848 | g->Draw("psame"); // colored markers | |
2849 | // gs->Draw("e[]psame"); // sys error filled rectangles | |
2850 | gc->Draw("epsame"); // open markers w/black stat error bars | |
2851 | graphs->Add(g); | |
2852 | leg->AddEntry(g,Form("%d-%d%%",centlow[centbins[k]],centhigh[centbins[k]]),"l,p"); | |
2853 | } | |
2854 | } | |
2855 | ||
2856 | // Adjust y limits | |
2857 | double newYmin = 0, newYmax = 0, newSpace = 0; | |
2858 | for (int m=0; m<graphs->GetEntries(); m++) { | |
2859 | TGraphErrors* tg = (TGraphErrors*)graphs->At(m); | |
2860 | double yhi = TMath::MaxElement(tg->GetN()-1, tg->GetY()); | |
2861 | double ylo = TMath::MinElement(tg->GetN()-1, tg->GetY()); | |
2862 | if (yhi > newYmax) newYmax = yhi; | |
2863 | if (ylo < newYmin) newYmin = ylo; | |
2864 | } | |
2865 | ||
2866 | newSpace = 0.1*(newYmax - newYmin); | |
2867 | hf->GetYaxis()->SetRangeUser(newYmin-newSpace, newYmax+newSpace); | |
2868 | if (opt.Contains("zoom")) | |
2869 | hf->GetYaxis()->SetRangeUser(-0.0008, 0.0085); | |
2870 | ||
2871 | leg->Draw(); | |
2872 | c->Update(); | |
2873 | ||
2874 | // Draw text | |
2875 | // if (opt.Contains("zoom")) | |
2876 | // ltx.SetTextSize(0.03); // 0.04 | |
2877 | // else | |
2878 | ||
2879 | // ltx.SetTextSize(0.04); // 0.05 | |
2880 | ||
2881 | double xstart = (opt.Contains("zoom")) ? 0.7 : 0.6; | |
2882 | ltx.DrawLatex(xstart, legbottom-1*0.07, Form("%.2g < p_{T} < %.2g GeV/c", t1, t2)); | |
2883 | // ltx.DrawLatex(xstart, legbottom-2*0.07, Form("%.2g < p_{T}^{assoc} < %.2g GeV/c", a1, a2)); | |
2884 | ||
2885 | return c; | |
2886 | } | |
2887 | ||
2888 | ||
2889 | TCanvas* DrawVnFromGlobalFit(int n, int ipt1, int ipt2, int ncb, | |
2890 | int centbins[], TString opt) | |
2891 | { | |
2892 | // Draw vn for one n on one canvas - many centralities | |
2893 | initialize(); | |
2894 | double maxPtFit = ipt2==999 ? pthigh[maxpt-1] : pthigh[ipt2]; | |
2895 | const char* cname = Form("v%d_global%s_ptfit%.2gto%.2g", | |
2896 | n, opt.Data(), ptlow[ipt1], maxPtFit); | |
2897 | const char* title = Form("v%d_global%s_ptfit%.2gto%.2g", | |
2898 | n, opt.Data(), ptlow[ipt1], maxPtFit); | |
2899 | double lvleft = opt.Contains("multi")? 0.74 : 0.82; | |
2900 | double lvbot = opt.Contains("multi")? 0.25 : 0.56; | |
2901 | int cvHeight = opt.Contains("multi")? 500 : 500; | |
2902 | TCanvas* cv = new TCanvas(cname, title, 700, cvHeight); | |
2903 | TLegend* lv = new TLegend(lvleft, lvbot, 0.99, 0.97); | |
2904 | lv->SetFillColor(kNone); lv->SetName("lv"); | |
2905 | lv->SetMargin(0.3); | |
2906 | lv->SetBorderSize(1); | |
2907 | lv->SetHeader("#splitline{Fit range}{(GeV/c):}"); | |
2908 | utils.padsetup(cv, 1, 1, "", 0.12, 0.01, 0.03, 0.2); // ###### | |
2909 | // gPad->SetRightMargin(0.15); | |
2910 | ||
2911 | ||
2912 | TPad* inset = 0; | |
2913 | ||
2914 | TObjArray* graphs = new TObjArray(); | |
2915 | double xmin = 0., xmax = 14.8; // 8.3; | |
2916 | double ymin = -0.06, ymax = 0.26; | |
2917 | if (opt.Contains("multi")) { | |
2918 | ymin = -1.19; | |
2919 | ymax = 0.52; | |
2920 | } | |
2921 | TH1F* hf = gPad->DrawFrame(xmin, ymin, xmax, ymax); | |
2922 | hf->SetTitle(Form(";p_{T}^{t} (GeV/c);v_{%d}{GF}",n)); | |
2923 | int nXticks = 210; | |
2924 | int nYticks = 210; | |
2925 | TLine l; | |
2926 | l.SetLineStyle(kDashed); | |
2927 | if (n==1) { | |
2928 | cv->cd(); | |
2929 | inset = new TPad("inset", "inset", 0.18, 0.18, 0.45, 0.59, kNone, 1, 0); | |
2930 | inset->Draw(); | |
2931 | inset->cd(); | |
2932 | TH1F* hfi = gPad->DrawFrame(-0.05, -0.04, 4.2, 0.16); | |
2933 | TAxis* hx = hfi->GetXaxis(); | |
2934 | TAxis* hy = hfi->GetYaxis(); | |
2935 | hx->SetLabelSize(0.1); | |
2936 | hy->SetLabelSize(0.1); | |
2937 | hx->SetNdivisions(104); | |
2938 | hy->SetNdivisions(104); | |
2939 | gPad->SetBottomMargin(0.2); | |
2940 | gPad->SetLeftMargin(0.2); | |
2941 | gPad->SetRightMargin(0.01); | |
2942 | gPad->SetTopMargin(0.01); | |
2943 | // utils.make_nice_axes(inset, hfi, 3, nXticks, nYticks, 1.4,-1.4); | |
2944 | l.DrawLine(-0.05, 0, 4.2, 0); | |
2945 | cv->cd(); | |
2946 | } | |
2947 | ||
2948 | utils.make_nice_axes(cv, hf, 1.3, nXticks, nYticks, 1.4,-1.4); | |
2949 | hf->SetLabelOffset(0.0, "X"); | |
2950 | hf->SetTitleOffset(1.5, "Y"); | |
2951 | hf->GetXaxis()->SetTicks("+"); | |
2952 | // hf->GetYaxis()->SetTicks("+-"); | |
2953 | ||
2954 | TGraph* v1h[5][3]; | |
2955 | if (opt.Contains("luzum")) { | |
2956 | int lSty[] = {7, kSolid, kSolid}; | |
2957 | int lCol[] = {kBlack, kCyan, kBlue, kGreen+1, kRed}; | |
2958 | TString pidString[] = {"pi", "K", "p"}; | |
2959 | for (int m=0; m<5; m++) { | |
2960 | for (int pid=0; pid<3; pid++) { | |
2961 | v1h[m][pid] = Luzumv1(m, pidString[pid].Data()); | |
2962 | TGraph* g1 = v1h[m][pid]; | |
2963 | g1->SetLineStyle(lSty[pid]); | |
2964 | g1->SetLineColor(lCol[m]); | |
2965 | g1->SetLineWidth(2); | |
2966 | } | |
2967 | } | |
2968 | hf->GetXaxis()->SetRangeUser(0, 4.8); | |
2969 | l.DrawLine(0, 0, 4.8, 0); | |
2970 | v1h[0][0]->Draw("plsame"); | |
2971 | v1h[0][2]->Draw("plsame"); | |
2972 | v1h[4][0]->Draw("plsame"); | |
2973 | v1h[4][2]->Draw("plsame"); | |
2974 | } | |
2975 | else | |
2976 | l.DrawLine(hf->GetXaxis()->GetXmin(), 0, hf->GetXaxis()->GetXmax(),0); | |
2977 | ||
2978 | for (int cb=0; cb<ncb; cb++) { | |
2979 | int k = centbins[cb]; | |
2980 | ||
2981 | if (!opt.Contains("multi")) { | |
2982 | // last 2 args are included pt bins | |
2983 | TGraphErrors* gtmp = vnGFvsPt(n, k, ipt1, ipt2); | |
2984 | ||
2985 | // clone to modify without ripples | |
2986 | TGraphErrors* gc = (TGraphErrors*) gtmp->Clone(); | |
2987 | if (!gc) { | |
2988 | Error("DrawVnFromGlobalFit()", "Problem making graph" ); | |
2989 | continue; | |
2990 | } | |
2991 | TGraphErrors* gtmp_sys = vnGFvsPt(n, k, ipt1, ipt2, "RIDGE", cfDef, "sys"); | |
2992 | TGraphErrors* gs = (TGraphErrors*) gtmp_sys->Clone(); | |
2993 | ||
2994 | if (opt.Contains("luzum")) { | |
2995 | if (centlow[k]==20) | |
2996 | utils.set_tgraph_props(gc, kBlue, kBlue, kFullCircle, 1.2); | |
2997 | if (centlow[k]==40) | |
2998 | utils.set_tgraph_props(gc, kRed, kRed, kFullCircle, 1.2); | |
2999 | } | |
3000 | ||
3001 | utils.draw_errorbox(gs, gs->GetFillColor()); | |
3002 | gc->Draw("epsame"); | |
3003 | graphs->Add(gc); | |
3004 | ||
3005 | TGraphErrors* open = (TGraphErrors*) gc->Clone(Form("%s_open", gc->GetName())); | |
3006 | utils.set_tgraph_props(open, kBlack, kBlack, kOpenCircle, 1.2); | |
3007 | open->Draw("epsame"); | |
3008 | ||
3009 | lv->AddEntry(gc,Form("%d-%d%%",centlow[k],centhigh[k]),"l,p"); | |
3010 | } | |
3011 | // NEW | |
3012 | if (opt.Contains("multi")) { | |
3013 | int b1,b2,b3,b4; | |
3014 | b1=PtBin(0.25,0.5); | |
3015 | b2=PtBin(0.75,1.0); | |
3016 | b3=PtBin(2.0, 2.5); | |
3017 | b4=PtBin(3.0, 4.0); | |
3018 | TGraphErrors* gg1 = vnGFvsPt(n, k, b1, b2, "RIDGE", cfDef, ""); | |
3019 | TGraphErrors* gg2 = vnGFvsPt(n, k, b3, b4, "RIDGE", cfDef, "keepsign"); | |
3020 | TGraphErrors* gs1 = vnGFvsPt(n, k, b1, b2, "RIDGE", cfDef, "sys"); | |
3021 | TGraphErrors* gs2 = vnGFvsPt(n, k, b3, b4, "RIDGE", cfDef, "sys_keepsign"); | |
3022 | utils.set_tgraph_props(gg1, gg1->GetLineColor(), gg1->GetMarkerColor(), kOpenCircle, 1.0); | |
3023 | utils.set_tgraph_props(gg2, gg2->GetLineColor(), gg2->GetMarkerColor(), kFullSquare, 1.0); | |
3024 | ||
3025 | utils.draw_errorbox(gs1, gs1->GetFillColor()); | |
3026 | utils.draw_errorbox(gs2, gs2->GetFillColor()); | |
3027 | gg1->Draw("epsame"); | |
3028 | gg2->Draw("epsame"); | |
3029 | graphs->Add(gs1); | |
3030 | graphs->Add(gs2); | |
3031 | const char* s1 = Form("#splitline{%.2g < p_{T}^{a} < %.2g }{%d-%d%%}", | |
3032 | ptlow[b1], pthigh[b2], centlow[k],centhigh[k]); | |
3033 | const char* s2 = Form("#splitline{%.2g < p_{T}^{a} < %.2g }{%d-%d%%}", | |
3034 | ptlow[b3], pthigh[b4], centlow[k],centhigh[k]); | |
3035 | lv->AddEntry(gg1, s1, "l,p"); | |
3036 | lv->AddEntry(gg2, s2, "l,p"); | |
3037 | ||
3038 | ||
3039 | if (n==1) { | |
3040 | inset->cd(); | |
3041 | gg1->Draw("epsame"); | |
3042 | gg2->Draw("epsame"); | |
3043 | cv->cd(); | |
3044 | } | |
3045 | ||
3046 | } | |
3047 | ||
3048 | } // centrality bin loop | |
3049 | ||
3050 | if (opt.Contains("luzum")) { | |
3051 | lv->AddEntry(v1h[0][0], Form("v.h. (#pi)"),"l,p"); | |
3052 | lv->AddEntry(v1h[0][2], Form("v.h. (p)"), "l,p"); | |
3053 | } | |
3054 | lv->Draw(); | |
3055 | ltx.SetTextSize(0.05); | |
3056 | ||
3057 | if (opt.Contains("multi")) { | |
3058 | gPad->SetTopMargin(0.02); | |
3059 | if (ispp) | |
3060 | ltx.DrawLatex(0.19, 0.92, Form("%s", "#splitline{pp}{2.76 TeV}")); | |
3061 | else | |
3062 | ltx.DrawLatex(0.19, 0.92, Form("%s", "#splitline{Pb-Pb}{2.76 TeV}")); | |
3063 | ltx.DrawLatex(0.49, 0.94, Form("%.2g < |#Delta#eta| < 1.8", minRidgeDEta) ); | |
3064 | ||
3065 | } | |
3066 | else { | |
3067 | ltx.DrawLatex(0.17, 0.88, Form("%.2g < |#Delta#eta| < 1.8", minRidgeDEta) ); | |
3068 | ltx.DrawLatex(0.17, 0.94, Form("%.2g < fit p_{T}^{a} < %.2g GeV/c", | |
3069 | ptlow[ipt1], maxPtFit)); | |
3070 | ltx.DrawLatex(0.19, 0.78, Form("%s", "#splitline{Pb-Pb}{2.76 TeV}")); | |
3071 | } | |
3072 | double newYmin = 0, newYmax = 0, newSpace = 0; | |
3073 | for (int m=0; m<graphs->GetEntries(); m++) { | |
3074 | TGraphErrors* tg = (TGraphErrors*)graphs->At(m); | |
3075 | double yhi = TMath::MaxElement(tg->GetN()-1, tg->GetY()); | |
3076 | double ylo = TMath::MinElement(tg->GetN()-1, tg->GetY()); | |
3077 | if (yhi > newYmax) newYmax = yhi; | |
3078 | if (ylo < newYmin) newYmin = ylo; | |
3079 | } | |
3080 | ||
3081 | if (opt.Contains("multi")) { | |
3082 | if(n==1) | |
3083 | newYmin *= 2.2; | |
3084 | if(n==2) | |
3085 | newYmin = 0.02; | |
3086 | } | |
3087 | ||
3088 | newSpace = 0.2*(newYmax - newYmin); | |
3089 | // if (opt.Contains("multi")) newSpace /= 2; | |
3090 | if (n==5) newSpace *= 2; | |
3091 | hf->GetYaxis()->SetRangeUser(newYmin-newSpace, newYmax+newSpace); | |
3092 | cv->Modified(); | |
3093 | cv->Update(); | |
3094 | cv->Draw(); | |
3095 | return cv; | |
3096 | } | |
3097 | ||
3098 | TCanvas* Drawv1to5(int ncb, int centbins[], int ipt1, int ipt2, TString opt) | |
3099 | { | |
3100 | if (!opt.IsNull()) | |
3101 | cout<<opt.Data()<<endl; | |
3102 | ||
3103 | initialize(); | |
3104 | int cent1 = centlow[0]; | |
3105 | int cent2 = centhigh[centbins[ncb-1]]; | |
3106 | double pmax = 15.; | |
3107 | if (ipt2 < 999.) | |
3108 | pmax = pthigh[ipt2]; | |
3109 | const char* cname = Form("vn_etamin%02d_cent%dto%d_fitptbin%dto%d", | |
3110 | (int)(10*minRidgeDEta), cent1, cent2, ipt1, ipt2); | |
3111 | const char* title = Form("v1to5gf_etamin%02d_cent%dto%d_fitpt%.2gto%.2g", | |
3112 | (int)(10*minRidgeDEta), cent1, cent2, ptlow[ipt1], pmax); | |
3113 | TCanvas* cv = new TCanvas(cname, title, 1400, 500); | |
3114 | TLegend* lv = new TLegend(0.09, 0.72, 0.27, 0.88); | |
3115 | lv->SetFillColor(kNone); | |
3116 | lv->SetBorderSize(0); | |
3117 | TLine zero; | |
3118 | double lv2bottom = 0.62; | |
3119 | TLegend* lv2 = new TLegend(0.82, lv2bottom, 0.93, 0.95, "Centrality"); | |
3120 | lv2->SetFillColor(kNone); | |
3121 | lv2->SetBorderSize(0); | |
3122 | utils.padsetup(cv, 5, 1, "", 0.12, 0.01, 0.03, 0.2); | |
3123 | ||
3124 | for (int n=1; n<=5; n++) { | |
3125 | Printf("Drawv1to5() n = %d / 5", n); | |
3126 | cv->cd(n); | |
3127 | double xmin = -0.43, xmax = 8.3; | |
3128 | double ymin = -0.09, ymax = 0.26; | |
3129 | if (ptlow[ipt1] > 4.) { | |
3130 | xmin = 4.9; | |
3131 | ymin = -0.1; | |
3132 | } | |
3133 | if (pthigh[ipt2] > 6.) { | |
3134 | xmax = 11.8; | |
3135 | ymax = 1; | |
3136 | } | |
3137 | if (opt.Contains("split")) { | |
3138 | xmin = -0.43; | |
3139 | xmax = 11.8; | |
3140 | ymin = -0.2; | |
3141 | ymax = 1; | |
3142 | } | |
3143 | ||
3144 | TH1F* hf = gPad->DrawFrame(xmin, ymin, xmax, ymax); | |
3145 | int nXticks = 206; | |
3146 | int nYticks = 208; | |
3147 | if (n==1) { | |
3148 | utils.make_nice_axes(cv, hf, 2., nXticks, nYticks, -0.02, 0.02); | |
3149 | hf->SetLabelOffset(-0.01, "X"); | |
3150 | } | |
3151 | else { | |
3152 | utils.make_nice_axes(cv, hf, 3.4, nXticks, nYticks, -0.075, 0.0); | |
3153 | hf->SetLabelOffset(-0.057, "X"); | |
3154 | hf->SetTickLength(0.02, "X"); | |
3155 | hf->SetTickLength(0.05, "Y"); | |
3156 | } | |
3157 | hf->GetXaxis()->SetTicks("+-"); | |
3158 | hf->GetYaxis()->SetTicks("+-"); | |
3159 | zero.SetLineStyle(kDashed); | |
3160 | zero.DrawLine(xmin, 0., xmax, 0.); | |
3161 | ||
3162 | // Loop thru every cent bin. Skip non-requested ones. | |
3163 | // for (int k=0; k<maxcent; k++) { | |
3164 | // bool centOK = 1; | |
3165 | // if (ncb != 999) { | |
3166 | // centOK = 0; | |
3167 | // for (int cb=0; cb<ncb; cb++) | |
3168 | // if(k==centbins[cb]) | |
3169 | // centOK = 1; | |
3170 | // } | |
3171 | // if (!centOK) | |
3172 | // continue; | |
3173 | ||
3174 | for (int cb=0; cb<ncb; cb++) { | |
3175 | int k = centbins[cb]; | |
3176 | TGraphErrors* gc = vnGFvsPt(n, k, ipt1, ipt2); | |
3177 | TGraphErrors* gs = vnGFvsPt(n, k, ipt1, ipt2, "RIDGE", cfDef, "sys"); | |
3178 | ||
3179 | if (!gc) { | |
3180 | Error("Drawv1to5()", "Problem finding graph"); | |
3181 | continue; | |
3182 | } | |
3183 | ||
3184 | TGraphErrors* gc_hipt=0; | |
3185 | TGraphErrors* gs_hipt=0; | |
3186 | if (opt.Contains("highptfit") || opt.Contains("split")) { | |
3187 | gc_hipt = vnGFvsPt(n, k, ipt2+1, 999, "highptfit"); | |
3188 | gs_hipt = vnGFvsPt(n, k, ipt2+1, 999, "RIDGE", cfDef, "sys"); | |
3189 | gc_hipt->SetMarkerStyle(kOpenCircle); | |
3190 | } | |
3191 | ||
3192 | int pale = gs->GetFillColor(); // it is CentColor(centlow[k], centhigh[k]). | |
3193 | utils.draw_errorbox(gs, pale); | |
3194 | gc->Draw("epsame"); | |
3195 | ||
3196 | if (opt.Contains("split")) { | |
3197 | utils.draw_errorbox(gs_hipt, pale); | |
3198 | gc_hipt->Draw("epsame"); | |
3199 | } | |
3200 | ||
3201 | TGraphErrors* open = (TGraphErrors*) gc->Clone(Form("%s_open", gc->GetName())); | |
3202 | utils.set_tgraph_props(open, kBlack, kBlack, kOpenCircle, 1.2); | |
3203 | open->Draw("epsame"); | |
3204 | ||
3205 | if (n==1) | |
3206 | lv2->AddEntry(gc,Form("%d-%d%%",centlow[k],centhigh[k]),"l,p"); | |
3207 | ||
3208 | } // cent bin k | |
3209 | } // Fourier moment n | |
3210 | ||
3211 | cv->cd(); | |
3212 | TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1); | |
3213 | overlay->SetFillStyle(4000); // transparent | |
3214 | overlay->Draw(); | |
3215 | overlay->cd(); | |
3216 | ||
3217 | ltx.SetTextSize(0.07); | |
3218 | lv2->Draw(); | |
3219 | ||
3220 | for (int n=1; n<=5; n++) { | |
3221 | cv->cd(n); | |
3222 | ltx.SetTextSize(n==1? 0.1 : 0.15); | |
3223 | ltx.DrawLatex((n==1)? 0.85 : 0.7, 0.86, Form("v_{%d}", n)); | |
3224 | ||
3225 | if (0 && n==1) | |
3226 | AddPrelimStampToCurrentPad(0.45,0.7,0.75,0.95, ""); | |
3227 | } | |
3228 | overlay->cd(); | |
3229 | ltx.SetTextSize(0.075); | |
3230 | ltx.DrawLatex(0.5, 0.05, "p_{T} [GeV/c]"); | |
3231 | ltx.SetTextAngle(90); | |
3232 | ltx.DrawLatex(0.05, 0.5, "v_{n }(p_{T})"); | |
3233 | ltx.SetTextAngle(0); | |
3234 | ltx.SetTextSize(0.05); | |
3235 | ltx.DrawLatex(0.14, 0.88, Form("%s", "#splitline{Pb-Pb}{2.76 TeV}")); | |
3236 | ltx.SetTextSize(0.05); | |
3237 | ltx.DrawLatex(0.13, 0.24, Form("%.2g < |#Delta#eta| < 1.8", minRidgeDEta) ); | |
3238 | ||
3239 | const char* frange = Form("%.2g < fit p_{T}^{a} < %.2g GeV", | |
3240 | ptlow[ipt1], pmax); | |
3241 | const char* frange1 = Form("#splitline{solid:}{%.2g < fit p_{T}^{a} < %.2g GeV}", | |
3242 | ptlow[ipt1], pthigh[ipt2]); | |
3243 | const char* frange2 = Form("#splitline{open:}{%.2g < fit p_{T}^{a} < %.2g GeV}", | |
3244 | ptlow[ipt2+1], pthigh[maxpt-1]); | |
3245 | ||
3246 | if (opt.Contains("split")) { | |
3247 | ltx.SetTextSize(0.04); | |
3248 | ltx.DrawLatex(0.31, 0.8, frange1); | |
3249 | ltx.DrawLatex(0.31, 0.7, frange2); | |
3250 | } | |
3251 | else | |
3252 | ltx.DrawLatex(0.3, 0.27, frange); | |
3253 | ||
3254 | return cv; | |
3255 | } | |
3256 | ||
3257 | TCanvas* Drawv2to5(int ncb, int centbins[], int ipt1, int ipt2, TString opt) | |
3258 | { | |
3259 | if (!opt.IsNull()) | |
3260 | cout<<opt.Data()<<endl; | |
3261 | ||
3262 | initialize(); | |
3263 | int cent1 = centlow[0]; | |
3264 | int cent2 = centhigh[centbins[ncb-1]]; | |
3265 | double pmax = 15.; | |
3266 | if (ipt2 < 999.) | |
3267 | pmax = pthigh[ipt2]; | |
3268 | const char* cname = Form("v2to5_etamin%02d_cent%dto%d_fitptbin%dto%d", | |
3269 | (int)(10*minRidgeDEta), cent1, cent2, ipt1, ipt2); | |
3270 | const char* title = Form("v2to5gf_etamin%02d_cent%dto%d_fitpt%.2gto%.2g", | |
3271 | (int)(10*minRidgeDEta), cent1, cent2, ptlow[ipt1], pmax); | |
3272 | TCanvas* cv = new TCanvas(cname, title, 1400, 500); | |
3273 | TLine zero; | |
3274 | double lv2bottom = 0.58; // 0.5; | |
3275 | TLegend* lv2 = new TLegend(0.77, lv2bottom, 0.9, 0.95, "Centrality"); | |
3276 | lv2->SetFillColor(kNone); | |
3277 | lv2->SetBorderSize(0); | |
3278 | utils.padsetup(cv, 4, 1, "", 0.08, 0.01, 0.03, 0.2); | |
3279 | ||
3280 | for (int n=2; n<=5; n++) { | |
3281 | Printf("Drawv2to5() n = %d", n); | |
3282 | cv->cd(n-1); | |
3283 | double xmin = -0.2; | |
3284 | double xmax = 7.6; | |
3285 | TH1F* hf = gPad->DrawFrame(xmin, -0.12, xmax, 0.26); //was y = -0.03, 0.26 | |
3286 | int nXticks = 206; | |
3287 | int nYticks = 208; | |
3288 | if (n==2) { | |
3289 | utils.make_nice_axes(cv, hf, 2., nXticks, nYticks, -0.02, 0.02); | |
3290 | hf->SetLabelOffset(0.0, "X"); // hf->SetLabelOffset(-0.01, "X"); | |
3291 | } | |
3292 | else { | |
3293 | utils.make_nice_axes(cv, hf, 2.6, nXticks, nYticks, -0.075, 0.0); | |
3294 | hf->SetLabelOffset(-0.02, "X"); | |
3295 | hf->SetTickLength(0.02, "X"); | |
3296 | hf->SetTickLength(0.05, "Y"); | |
3297 | } | |
3298 | hf->GetXaxis()->SetTicks("+-"); | |
3299 | hf->GetYaxis()->SetTicks("+-"); | |
3300 | zero.SetLineStyle(kDashed); | |
3301 | zero.DrawLine(xmin, 0., xmax, 0.); | |
3302 | ||
3303 | for (int cb=0; cb<ncb; cb++) { | |
3304 | int k = centbins[cb]; | |
3305 | ||
3306 | TGraphErrors* gc = vnGFvsPt(n, k, ipt1, ipt2); | |
3307 | TGraphErrors* gs = vnGFvsPt(n, k, ipt1, ipt2, "RIDGE", cfDef, "sys"); | |
3308 | ||
3309 | if (!gc) { | |
3310 | Error("Drawv2to5()", "Problem finding graph"); | |
3311 | continue; | |
3312 | } | |
3313 | ||
3314 | int pale = gs->GetFillColor(); // it is CentColor(centlow[k], centhigh[k]). | |
3315 | utils.draw_errorbox(gs, pale); | |
3316 | // gs->Draw("e2psame"); | |
3317 | gc->Draw("epsame"); | |
3318 | TGraphErrors* open = (TGraphErrors*) gc->Clone(Form("%s_open", gc->GetName())); | |
3319 | utils.set_tgraph_props(open, kBlack, kBlack, kOpenCircle, 1.2); | |
3320 | open->Draw("epsame"); | |
3321 | ||
3322 | if (n==2) | |
3323 | lv2->AddEntry(gc,Form("%d-%d%%",centlow[k],centhigh[k]),"l,p"); | |
3324 | ||
3325 | } // cent bin k | |
3326 | } // Fourier moment n | |
3327 | ||
3328 | cv->cd(); | |
3329 | TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1); | |
3330 | overlay->SetFillStyle(4000); // transparent | |
3331 | overlay->Draw(); | |
3332 | overlay->cd(); | |
3333 | ||
3334 | ltx.SetTextSize(0.07); | |
3335 | lv2->Draw(); | |
3336 | ||
3337 | for (int n=2; n<=5; n++) { | |
3338 | cv->cd(n-1); | |
3339 | ltx.SetTextSize(n==2? 0.1 : 0.15); | |
3340 | ltx.DrawLatex((n==2)? 0.85 : 0.7, 0.86, Form("v_{%d}", n)); | |
3341 | } | |
3342 | overlay->cd(); | |
3343 | ltx.SetTextSize(0.075); | |
3344 | ltx.DrawLatex(0.48, 0.04, "p_{T} [GeV/c]"); | |
3345 | ltx.SetTextAngle(90); | |
3346 | ltx.DrawLatex(0.02, 0.5, "v_{n }{GF}"); | |
3347 | ltx.SetTextAngle(0); | |
3348 | ltx.SetTextSize(0.05); | |
3349 | ltx.DrawLatex(0.33, 0.88, Form("%s", "#splitline{Pb-Pb}{2.76 TeV}")); | |
3350 | ltx.SetTextSize(0.05); | |
3351 | ltx.DrawLatex(0.33, 0.75, Form("%.2g < |#Delta#eta| < 1.8", minRidgeDEta) ); | |
3352 | ||
3353 | const char* frange = Form("#splitline{Global fit range}{%.2g < p_{T} < %.2g GeV}", | |
3354 | ptlow[ipt1], pmax); | |
3355 | ltx.DrawLatex(0.56, 0.75, frange); | |
3356 | ||
3357 | return cv; | |
3358 | } | |
3359 | ||
3360 | TCanvas* DrawChi2vsPtTrig(int npta, int ptabins[], int ncb, int centbins[], TString opt) | |
3361 | { | |
3362 | ||
3363 | // TODO: options-- global or 1to5, chi2/N or prob | |
3364 | double a1=ptlow[ptabins[0]], a2=pthigh[ptabins[npta-1]]; | |
3365 | int c1=centlow[centbins[0]], c2=centhigh[centbins[ncb-1]]; | |
3366 | TString cname = Form("chi2_etamin%02d_pta%.2gto%.2g_cent%dto%d", | |
3367 | (int)(10*minRidgeDEta), a1, a2, c1, c2 ); | |
3368 | ||
3369 | if (!opt.IsNull() ) { | |
3370 | cname.Append("_"); | |
3371 | cname.Append(opt.Data()); | |
3372 | } | |
3373 | ||
3374 | TCanvas* cv = new TCanvas(cname, cname, 1); | |
3375 | TLegend* lv = new TLegend(0.09, 0.72, 0.27, 0.88); // w.r.t. cv | |
3376 | lv->SetFillColor(kNone); | |
3377 | lv->SetBorderSize(0); | |
3378 | ||
3379 | double lv2bottom = 0.62; | |
3380 | TLegend* lv2 = new TLegend(0.82, lv2bottom, 0.93, 0.95, "Centrality"); // w.r.t. cv | |
3381 | lv2->SetFillColor(kNone); | |
3382 | lv2->SetBorderSize(0); | |
3383 | utils.padsetup(cv, 1, 1, "", 0.12, 0.01, 0.03, 0.2); | |
3384 | ||
3385 | double xmin = -0.43; | |
3386 | double xmax = 12.5; | |
3387 | TH1F* hf = gPad->DrawFrame(xmin, -0.09, xmax, 26); | |
3388 | // int nXticks = 210; | |
3389 | // int nYticks = 210; | |
3390 | hf->GetXaxis()->SetTicks("+-"); | |
3391 | hf->GetYaxis()->SetTicks("+-"); | |
3392 | ||
3393 | // IN PROGRESS--------------- | |
3394 | for (int j=0; j<npta; j++) { | |
3395 | for(int k=0; k<ncb; k++) { | |
3396 | ||
3397 | TGraphErrors *gc = CorrFnChi2VsTrigPt(ptabins[j], centbins[k], 1, 5, "ndf"); | |
3398 | if (!gc) { | |
3399 | Error("DrawChi2()", "!gc"); | |
3400 | continue; | |
3401 | } | |
3402 | int col = CentColor(centlow[centbins[k]], centhigh[centbins[k]]); | |
3403 | utils.set_tgraph_props(gc, col, col, kFullCircle, 1.6); | |
3404 | gc->Draw("epsame"); | |
3405 | lv2->AddEntry(gc, centLabel(centbins[k]), "l,p"); | |
3406 | } | |
3407 | } | |
3408 | ||
3409 | // cv->cd(); | |
3410 | // TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1); | |
3411 | // overlay->SetFillStyle(4000); // transparent | |
3412 | // overlay->Draw(); | |
3413 | // overlay->cd(); | |
3414 | ||
3415 | // ltx.SetTextSize(0.07); | |
3416 | // lv2->Draw(); | |
3417 | // overlay->cd(); | |
3418 | // ltx.SetTextSize(0.075); | |
3419 | // ltx.DrawLatex(0.5, 0.05, "trigger p_{T} [GeV/c]"); | |
3420 | // ltx.SetTextAngle(90); | |
3421 | // ltx.DrawLatex(0.05, 0.5, "#chi^{2} / NDF"); | |
3422 | // ltx.SetTextAngle(0); | |
3423 | ||
3424 | return cv; | |
3425 | } | |
3426 | ||
3427 | ||
3428 | TCanvas* DrawAgreement(int icent, TString opt) | |
3429 | { | |
3430 | // Draw Chi2 on one canvas for each | |
3431 | // centrality bin k. | |
3432 | ||
3433 | initialize(); | |
3434 | int cent1 = centlow[icent]; | |
3435 | int cent2 = centhigh[icent]; | |
3436 | ||
3437 | if (icent==999) { | |
3438 | cent1 = 0; | |
3439 | cent2 = 90; | |
3440 | } | |
3441 | int k = icent; | |
3442 | TString cname = Form("agmt_etamin%02d_cent%dto%d_%s", | |
3443 | (int)(10*minRidgeDEta), cent1, cent2, opt.Data()); | |
3444 | if (!opt.IsNull()) { | |
3445 | cname.Append("_"); | |
3446 | cname.Append(opt.Data()); | |
3447 | } | |
3448 | ||
3449 | TCanvas* cv = new TCanvas(cname.Data(), cname.Data(), 1400, 500); | |
3450 | TLegend* lv = new TLegend(0.09, 0.72, 0.27, 0.88); | |
3451 | lv->SetFillColor(kNone); | |
3452 | lv->SetBorderSize(0); | |
3453 | ||
3454 | double lv2bottom = 0.62; | |
3455 | TLegend* lv2 = new TLegend(0.82, lv2bottom, 0.93, 0.95, "Centrality"); | |
3456 | lv2->SetFillColor(kNone); | |
3457 | lv2->SetBorderSize(0); | |
3458 | cv->Divide(5,1); | |
3459 | ||
3460 | for (int n=1; n<=5; n++) { | |
3461 | cv->cd(n); | |
3462 | ||
3463 | TH2F *gc = Agreement2DHist(k,n); | |
3464 | if (!gc) { | |
3465 | Error("DrawAgreement()", "Problem finding graph" ); | |
3466 | continue; | |
3467 | } | |
3468 | ||
3469 | // gc->GetZaxis()->SetRangeUser(1.e-1, 1000); | |
3470 | ||
3471 | gPad->SetLogx(); gPad->SetLogy(); gPad->SetLogz(); | |
3472 | ||
3473 | if (opt.Contains("lego")) { | |
3474 | gc->SetTitleOffset(0.8, "X"); | |
3475 | gc->GetXaxis()->SetLabelSize(0.07); | |
3476 | gc->GetYaxis()->SetLabelSize(0.07); | |
3477 | gc->GetZaxis()->SetLabelSize(0.07); | |
3478 | gc->GetXaxis()->SetTitleSize(0.1); | |
3479 | gc->GetYaxis()->SetTitleSize(0.1); | |
3480 | gc->GetZaxis()->SetTitleSize(0.1); | |
3481 | gc->GetXaxis()->SetNdivisions(8); | |
3482 | gc->GetYaxis()->SetNdivisions(8); | |
3483 | ||
3484 | gPad->SetTopMargin(0.2); | |
3485 | gPad->SetLeftMargin(0.12); | |
3486 | gPad->SetRightMargin(0.01); | |
3487 | gc->Draw("lego2"); | |
3488 | } | |
3489 | else if (opt.Contains("contour")) | |
3490 | gc->Draw("colz"); | |
3491 | ||
3492 | if (n==1) { | |
3493 | lv2->AddEntry(gc, Form("%d-%d%%",centlow[k],centhigh[k]),"l,p"); | |
3494 | } | |
3495 | ||
3496 | } // Fourier moment n | |
3497 | ||
3498 | cv->cd(); | |
3499 | TPad* overlay = new TPad("overlay", "overlay", 0,0,1,1); | |
3500 | overlay->SetFillStyle(4000); // transparent | |
3501 | overlay->Draw(); | |
3502 | overlay->cd(); | |
3503 | ||
3504 | ltx.SetTextSize(0.07); | |
3505 | const char* str = Form("#left|V_{n#Delta} / (v_{n}^{t} v_{n}^{a})#right|"); | |
3506 | ltx.DrawLatex(0.42, 0.92, str); | |
3507 | ltx.DrawLatex(0.91, 0.92, Form("%s", centLabel(k))); | |
3508 | ||
3509 | for (int n=1; n<=5; n++) { | |
3510 | cv->cd(n); | |
3511 | ltx.SetTextSize(0.10); | |
3512 | ltx.DrawLatex(0.05, 0.76, Form("n = %d", n)); | |
3513 | } | |
3514 | return cv; | |
3515 | } | |
3516 | ||
3517 | void AddPrelimStampToCurrentPad(float x1, float y1, float x2, float y2, TString opt) | |
3518 | { | |
3519 | TASImage* logo = new TASImage("../commonfigs/alice-prelim.png"); // 287x323 | |
3520 | TPad* logopad = new TPad("logopad", "logo", x1,y1,x2,y2); | |
3521 | if (opt.Contains("border")) { | |
3522 | logopad->SetFillColor(kRed); | |
3523 | logopad->SetBorderSize(2); | |
3524 | } | |
3525 | logopad->Draw(); | |
3526 | logopad->cd(); | |
3527 | logo->Draw(""); | |
3528 | return; | |
3529 | } | |
3530 | ||
3531 | int CentBin(double cntLo, double cntHi) | |
3532 | { | |
3533 | int bin = -123; | |
3534 | for (int k=0; k<maxcent; k++) { | |
3535 | if (centlow[k]==cntLo && centhigh[k]==cntHi) | |
3536 | bin = k; | |
3537 | } | |
3538 | if (bin < 0 && 0) | |
3539 | Warning("CentBin","No bin for %.2g - %.2g", cntLo, cntHi); | |
3540 | ||
3541 | return bin; | |
3542 | } | |
3543 | ||
3544 | int PtBin(double ptLo, double ptHi) | |
3545 | { | |
3546 | int bin = -123; | |
3547 | for (int i=0; i<maxpt; i++) { | |
3548 | if (ptlow[i]==ptLo && pthigh[i]==ptHi) | |
3549 | bin = i; | |
3550 | } | |
3551 | if (bin < 0) | |
3552 | Warning("PtBin","No bin for %.2g - %.2g", ptLo, ptHi); | |
3553 | ||
3554 | return bin; | |
3555 | } | |
3556 | ||
3557 | double MomCorrection(int i, int j, int k) | |
3558 | { | |
3559 | if(0)cout<<i<<j<<k<<endl; | |
3560 | // For momentum (non)conservation in v1. | |
3561 | // Multiplicitave weight: w = w_t*w_a. | |
3562 | // w_t = ptt - <pt^2>/<pt> | |
3563 | // w_a = pta - <pt^2>/<pt> | |
3564 | ||
3565 | // double c = 1./3 * -0.00185; | |
3566 | ||
3567 | // (ptmean[i]-)* ptmean[j] | |
3568 | ||
3569 | double correction = 1./3 * -0.00185 * ptmean[i]* ptmean[j]; | |
3570 | return correction; | |
3571 | } | |
3572 | ||
3573 | void MakeVnDVsQGraphs(int n1, int n2, int k, const char* region, const char* corrtype, TString opt) | |
3574 | { | |
3575 | // The first two errType constants have points set at the central measured value of | |
3576 | // VnDelta, with statistical and systematic error bars | |
3577 | // respectively. The last two have points set at the upper and lower | |
3578 | // systematic value, with statistical error bars set to be the same | |
3579 | // as for meas_stat. | |
3580 | const int nErrTypes = 4; | |
3581 | enum errType {MEAS_STAT, MEAS_SYS, HI, LO}; | |
3582 | // const char* errTypeStr[] = {"meas_stat","meas_sys","hi","lo"}; | |
3583 | ||
3584 | // Fill a huge array with all the coefficients, then use it to set the tgraphs. | |
3585 | const int nN = gNMax + 1; // 1-10, don't use 0. | |
3586 | const int nI = maxpt+1; | |
3587 | const int nJ = maxpt+1; | |
3588 | const int nT = nErrTypes; | |
3589 | double vVal[nN][nI][nJ][nT]; | |
3590 | double vErr[nN][nI][nJ][nT]; | |
3591 | double vMix[nN][nI][nJ][nT]; // Just the acceptance correction unc. | |
3592 | ||
3593 | double vSin[nN][nI][nJ][nT]; // <sin n delta phi> for sys. estimation | |
3594 | double vSinErr[nN][nI][nJ][nT]; | |
3595 | ||
3596 | for (int i=0; i<maxpt; i++) { | |
3597 | for (int j=0; j<=i; j++) { | |
3598 | ||
3599 | double esin = 0; // residual sine error | |
3600 | double ebw = 0; // bin width error \propto n / (# dphi bins) | |
3601 | double ecnt = 0; // centrality determination | |
3602 | double emix = 0; // acc. correction error | |
3603 | double etrk = 0; // tracking resolution | |
3604 | double ev1 = 0; // v1 | |
3605 | double esys = 0; // quadrature sum | |
3606 | ||
3607 | TH1* hst=0, *hsys=0; // stat and sys error, but same points | |
3608 | ||
3609 | hst = Hist(region, corrtype, i, j, k, ""); | |
3610 | hsys = Hist(region, corrtype, i, j, k, "sys"); | |
3611 | ||
3612 | if (!hst) Error("MakeVnDVsQGraphs()","!hst"); | |
3613 | if (!hsys) Error("MakeVnDVsQGraphs()","!hsys"); | |
3614 | ||
3615 | esin = SineError(i,j,k); // resid_sin / 2; | |
3616 | ||
3617 | for (int n=n1; n<=n2; n++) { // Order of coeff. | |
3618 | ||
3619 | // VnDelta on measured points, with statistical errors. | |
3620 | double val=0, err=0, resid_sin=0, rs_err=0; | |
3621 | VnDelta(n, *hst, val, err, ""); | |
3622 | VnDelta(n, *hst, resid_sin, rs_err, "sin"); | |
3623 | ||
3624 | // VnDelta on high and low sys values. Stat err. not used. | |
3625 | double val_hi=0, xxx=0; | |
3626 | VnDelta(n, *hsys, val_hi, xxx, "hi"); | |
3627 | double val_lo=0, yyy=0; | |
3628 | VnDelta(n, *hsys, val_lo, yyy, "lo"); | |
3629 | ||
3630 | // Momentum conservation correction | |
3631 | if (n==1 && opt.Contains("ptcons")) { // eqs. 3-7, PRL 106, 102301 (2011) | |
3632 | double factor = ispp ? 0.212 : 0; // 7e-5; <------ trial & error: 7e-5 must be close for 0-20% | |
3633 | val += factor*ptmean[i]*ptmean[j]; | |
3634 | val_hi += factor*ptmean[i]*ptmean[j]; | |
3635 | val_lo += factor*ptmean[i]*ptmean[j]; | |
3636 | } | |
3637 | ||
3638 | int nbins = hst->GetNbinsX(); | |
3639 | double mean = (val_hi + val_lo + val) / 3.; | |
3640 | double ms = | |
3641 | (val_hi-mean)*(val_hi-mean) + | |
3642 | (val_lo-mean)*(val_lo-mean) + | |
3643 | (val -mean)*(val -mean); | |
3644 | double rms = TMath::Sqrt(ms/2); // sqrt(1/(N-1) sum (x_i - <x>)^2 ) | |
3645 | ||
3646 | ||
3647 | ebw = n/nbins/TMath::Sqrt(12)*val; // This is about 0.008 * n * val | |
3648 | etrk = 0.01*ptmean[j]*val; // see comments in comparisons/VnDcomparison.C | |
3649 | ecnt = 0.01*val; | |
3650 | emix = rms/2; // take half not to double-count stat err. // 0.25*TMath::Abs(val_hi - val_lo); | |
3651 | ev1 = (n==1) ? 0.001*ptmean[i]*ptmean[j]*val : 0.; | |
3652 | ||
3653 | esys = TMath::Sqrt(ecnt*ecnt + emix*emix + ebw*ebw + etrk*etrk + ev1*ev1 + esin*esin); | |
3654 | ||
3655 | double pointVal = val; | |
3656 | double pointErr = err; | |
3657 | ||
3658 | vVal[n][i][j][MEAS_STAT] = pointVal; | |
3659 | vErr[n][i][j][MEAS_STAT] = pointErr; | |
3660 | vVal[n][i][j][MEAS_SYS] = pointVal; | |
3661 | vErr[n][i][j][MEAS_SYS] = esys; | |
3662 | vVal[n][i][j][HI] = pointVal + esys; | |
3663 | vErr[n][i][j][HI] = pointErr; | |
3664 | vVal[n][i][j][LO] = pointVal - esys; | |
3665 | vErr[n][i][j][LO] = pointErr; | |
3666 | vSin[n][i][j][MEAS_STAT] = resid_sin; | |
3667 | vSinErr[n][i][j][MEAS_STAT] = rs_err; | |
3668 | vMix[n][i][j][MEAS_SYS] = emix; | |
3669 | ||
3670 | } // n | |
3671 | } // j | |
3672 | } // i | |
3673 | ||
3674 | // Now set the graphs from the arrays. | |
3675 | for (int n=n1; n<=n2; n++) { | |
3676 | for (int t=0; t<nErrTypes; t++) { | |
3677 | ||
3678 | TString gname = Form("V%dDeltaVsGlobalIndex_cent%dto%d", | |
3679 | n,centlow[k],centhigh[k]); | |
3680 | if (t==HI) | |
3681 | gname.Append("hi"); | |
3682 | else if (t==LO) | |
3683 | gname.Append("lo"); | |
3684 | else if (t==MEAS_SYS) | |
3685 | gname.Append("meas_sys"); | |
3686 | else | |
3687 | gname.Append("meas_stat"); | |
3688 | ||
3689 | if (opt.Contains("sine")) | |
3690 | gname.Append("_sine"); | |
3691 | ||
3692 | if (opt.Contains("ALL")) | |
3693 | gname.Append("_ALL"); | |
3694 | ||
3695 | TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data()); | |
3696 | if (!g) { | |
3697 | g = new TGraphErrors(); | |
3698 | g->SetName(gname.Data()); | |
3699 | } | |
3700 | ||
3701 | TGraphErrors* gmix = 0; | |
3702 | if (t==MEAS_SYS) { | |
3703 | gmix = new TGraphErrors(); | |
3704 | gname.Append("_mix"); | |
3705 | gmix->SetName(gname.Data()); | |
3706 | } | |
3707 | ||
3708 | double xer = 0.4; // was 0.5 | |
3709 | for (int i=0; i<maxpt; i++) { | |
3710 | for (int j=0; j<=i; j++) { | |
3711 | int q = GlobalIndex(i, j); | |
3712 | g->SetPoint(q, q+0.5, vVal[n][i][j][t]); | |
3713 | g->SetPointError(q, xer, vErr[n][i][j][t]); | |
3714 | ||
3715 | if (opt.Contains("sine")) { | |
3716 | g->SetPoint(q, q+0.5, vSin[n][i][j][MEAS_STAT]); | |
3717 | g->SetPointError(q, xer, vSinErr[n][i][j][MEAS_STAT]); | |
3718 | } | |
3719 | ||
3720 | if (t==MEAS_SYS) { | |
3721 | gmix->SetPoint(q, q+0.5, vVal[n][i][j][t]); | |
3722 | gmix->SetPointError(q, xer, vMix[n][i][j][t]); | |
3723 | } | |
3724 | ||
3725 | } | |
3726 | } | |
3727 | ||
3728 | gList->Add(g); | |
3729 | if (gmix) gList->Add(gmix); | |
3730 | } | |
3731 | } | |
3732 | ||
3733 | return; | |
3734 | } | |
3735 | ||
3736 | TGraphErrors* VnDVsQ(int n, int k, const char* region, const char* corrtype, TString opt) | |
3737 | { | |
3738 | ||
3739 | // Graph of pair fourier coefficients VnDelta vs. the linear | |
3740 | // global index q. Arguments: | |
3741 | // | |
3742 | // n, k: fourier index, centrality bin | |
3743 | // region: "RIDGE", "NSJET" or "ALL" | |
3744 | // corrtype: "s", "m", "cA", "cB" | |
3745 | // | |
3746 | // Select opt "hi_sys" or "lo_sys" for hi and low pts | |
3747 | // with central stat errors and "sys" for the central points | |
3748 | // with sys. errors. Just like Hist(). | |
3749 | ||
3750 | TString gname = Form("V%dDeltaVsGlobalIndex_cent%dto%d", | |
3751 | n,centlow[k],centhigh[k]); | |
3752 | ||
3753 | if (opt.Contains("hi")) | |
3754 | gname.Append("hi"); | |
3755 | else if (opt.Contains("lo")) | |
3756 | gname.Append("lo"); | |
3757 | else if (opt.Contains("sys")) | |
3758 | gname.Append("meas_sys"); | |
3759 | else | |
3760 | gname.Append("meas_stat"); | |
3761 | ||
3762 | if (opt.Contains("sine")) | |
3763 | gname.Append("_sine"); | |
3764 | ||
3765 | if (opt.Contains("ALL")) | |
3766 | gname.Append("_ALL"); | |
3767 | ||
3768 | TGraphErrors* g = (TGraphErrors*) gList->FindObject(gname.Data()); | |
3769 | if (g) { | |
3770 | if (0) | |
3771 | Info("VnDVsQ()", "Returning preexisting graph"); | |
3772 | return g; | |
3773 | } | |
3774 | ||
3775 | TString opt1 = (opt.Contains("sine")) ? "sine" : ""; | |
3776 | if (opt.Contains("ALL")) | |
3777 | opt1.Append("_ALL"); | |
3778 | if (opt.Contains("ptcons")) | |
3779 | opt1.Append("_ptcons"); | |
3780 | ||
3781 | MakeVnDVsQGraphs(1, gNMax, k, region, corrtype, opt1); | |
3782 | g = (TGraphErrors*) gList->FindObject(gname.Data()); | |
3783 | if (!g) { | |
3784 | Error("VnDVsQ()", "No graph %s", gname.Data()); | |
3785 | gSystem->Exit(-1); | |
3786 | } | |
3787 | return g; | |
3788 | } | |
3789 | ||
3790 | double MeanPt(int i, int j, int k, TString t_or_a, TString opt) | |
3791 | { | |
3792 | // For combined centralities, map k's to the closest available bin | |
3793 | // double cb1[] = {0, 0, 0,2,2, 1,3,0, 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 60 }; | |
3794 | // double cb2[] = {10,20,2,5,10,3,5,5, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 90 }; | |
3795 | // OBJ: TProfile2D hMeanPtTrgCen1 cen=1 (0.5) : 0 at: 0x102af06e0 | |
3796 | // OBJ: TProfile2D hMeanPtAssCen1 cen=1 (0.5) : 0 at: 0x108308e60 | |
3797 | // OBJ: TProfile2D hMean2PtTrgCen1 cen=1 (0.5) : 0 at: 0x108309480 | |
3798 | // OBJ: TProfile2D hMean2PtAssCen1 cen=1 (0.5) : 0 at: 0x108309ad0 | |
3799 | // OBJ: TProfile2D hMeanPtTrgCen2 cen=2 (1.5) : 0 at: 0x10830a120 | |
3800 | // OBJ: TProfile2D hMeanPtAssCen2 cen=2 (1.5) : 0 at: 0x10830a770 | |
3801 | // OBJ: TProfile2D hMean2PtTrgCen2 cen=2 (1.5) : 0 at: 0x10830adc0 | |
3802 | // OBJ: TProfile2D hMean2PtAssCen2 cen=2 (1.5) : 0 at: 0x10830b410 | |
3803 | // OBJ: TProfile2D hMeanPtTrgCen3 cen=3 (2.5) : 0 at: 0x10830ba60 | |
3804 | // OBJ: TProfile2D hMeanPtAssCen3 cen=3 (2.5) : 0 at: 0x10830c0b0 | |
3805 | // OBJ: TProfile2D hMean2PtTrgCen3 cen=3 (2.5) : 0 at: 0x10830c700 | |
3806 | // OBJ: TProfile2D hMean2PtAssCen3 cen=3 (2.5) : 0 at: 0x10830cd50 | |
3807 | // OBJ: TProfile2D hMeanPtTrgCen4 cen=4 (3.5) : 0 at: 0x10830d3a0 | |
3808 | // OBJ: TProfile2D hMeanPtAssCen4 cen=4 (3.5) : 0 at: 0x10830d9f0 | |
3809 | // OBJ: TProfile2D hMean2PtTrgCen4 cen=4 (3.5) : 0 at: 0x10830e040 | |
3810 | // OBJ: TProfile2D hMean2PtAssCen4 cen=4 (3.5) : 0 at: 0x10830e690 | |
3811 | // OBJ: TProfile2D hMeanPtTrgCen5 cen=5 (4.5) : 0 at: 0x10830ece0 | |
3812 | // OBJ: TProfile2D hMeanPtAssCen5 cen=5 (4.5) : 0 at: 0x10830f330 | |
3813 | // OBJ: TProfile2D hMean2PtTrgCen5 cen=5 (4.5) : 0 at: 0x10830f980 | |
3814 | // OBJ: TProfile2D hMean2PtAssCen5 cen=5 (4.5) : 0 at: 0x10830ffd0 | |
3815 | // OBJ: TProfile2D hMeanPtTrgCen6 cen=6 (7.5) : 0 at: 0x108310620 | |
3816 | // OBJ: TProfile2D hMeanPtAssCen6 cen=6 (7.5) : 0 at: 0x108310c70 | |
3817 | // OBJ: TProfile2D hMean2PtTrgCen6 cen=6 (7.5) : 0 at: 0x1083112c0 | |
3818 | // OBJ: TProfile2D hMean2PtAssCen6 cen=6 (7.5) : 0 at: 0x108311910 | |
3819 | // OBJ: TProfile2D hMeanPtTrgCen7 cen=7 (15.0) : 0 at: 0x108311f60 | |
3820 | // OBJ: TProfile2D hMeanPtAssCen7 cen=7 (15.0) : 0 at: 0x1083125b0 | |
3821 | // OBJ: TProfile2D hMean2PtTrgCen7 cen=7 (15.0) : 0 at: 0x108312c00 | |
3822 | // OBJ: TProfile2D hMean2PtAssCen7 cen=7 (15.0) : 0 at: 0x108313250 | |
3823 | // OBJ: TProfile2D hMeanPtTrgCen8 cen=8 (25.0) : 0 at: 0x1083138a0 | |
3824 | // OBJ: TProfile2D hMeanPtAssCen8 cen=8 (25.0) : 0 at: 0x108313ef0 | |
3825 | // OBJ: TProfile2D hMean2PtTrgCen8 cen=8 (25.0) : 0 at: 0x108314540 | |
3826 | // OBJ: TProfile2D hMean2PtAssCen8 cen=8 (25.0) : 0 at: 0x108314b90 | |
3827 | // OBJ: TProfile2D hMeanPtTrgCen9 cen=9 (35.0) : 0 at: 0x1083151e0 | |
3828 | // OBJ: TProfile2D hMeanPtAssCen9 cen=9 (35.0) : 0 at: 0x108315830 | |
3829 | // OBJ: TProfile2D hMean2PtTrgCen9 cen=9 (35.0) : 0 at: 0x108315e80 | |
3830 | // OBJ: TProfile2D hMean2PtAssCen9 cen=9 (35.0) : 0 at: 0x1083164d0 | |
3831 | // OBJ: TProfile2D hMeanPtTrgCen10 cen=10 (45.0) : 0 at: 0x108316b20 | |
3832 | // OBJ: TProfile2D hMeanPtAssCen10 cen=10 (45.0) : 0 at: 0x108317170 | |
3833 | // OBJ: TProfile2D hMean2PtTrgCen10 cen=10 (45.0) : 0 at: 0x1083177c0 | |
3834 | // OBJ: TProfile2D hMean2PtAssCen10 cen=10 (45.0) : 0 at: 0x108317e10 | |
3835 | // OBJ: TProfile2D hMeanPtTrgCen11 cen=11 (55.0) : 0 at: 0x108318460 | |
3836 | // OBJ: TProfile2D hMeanPtAssCen11 cen=11 (55.0) : 0 at: 0x108318ab0 | |
3837 | // OBJ: TProfile2D hMean2PtTrgCen11 cen=11 (55.0) : 0 at: 0x108319100 | |
3838 | // OBJ: TProfile2D hMean2PtAssCen11 cen=11 (55.0) : 0 at: 0x108319750 | |
3839 | // OBJ: TProfile2D hMeanPtTrgCen12 cen=12 (75.0) : 0 at: 0x108319da0 | |
3840 | // OBJ: TProfile2D hMeanPtAssCen12 cen=12 (75.0) : 0 at: 0x10831a3f0 | |
3841 | // OBJ: TProfile2D hMean2PtTrgCen12 cen=12 (75.0) : 0 at: 0x10831aa40 | |
3842 | // OBJ: TProfile2D hMean2PtAssCen12 cen=12 (75.0) : 0 at: 0x10831b090 | |
3843 | ||
3844 | int k2 = -1; | |
3845 | if (k == CentBin(0, 1)) k2 = 1; | |
3846 | else if (k == CentBin(1, 2)) k2 = 2; | |
3847 | else if (k == CentBin(2, 3)) k2 = 3; | |
3848 | else if (k == CentBin(3, 4)) k2 = 4; | |
3849 | else if (k == CentBin(4, 5)) k2 = 5; | |
3850 | else if (k == CentBin(5, 10)) k2 = 6; | |
3851 | else if (k == CentBin(10,20)) k2 = 7; | |
3852 | else if (k == CentBin(20,30)) k2 = 8; | |
3853 | else if (k == CentBin(30,40)) k2 = 9; | |
3854 | else if (k == CentBin(40,50)) k2 = 10; | |
3855 | else if (k == CentBin(50,60)) k2 = 11; | |
3856 | else if (k == CentBin(60,90)) k2 = 12; | |
3857 | else if (k == CentBin(0, 10)) k2 = 5; //CentBin(4, 5); | |
3858 | else if (k == CentBin(0, 20)) k2 = 6; //CentBin(5, 10); | |
3859 | else if (k == CentBin(0, 2)) k2 = 1; //CentBin(0, 1); | |
3860 | else if (k == CentBin(2, 5)) k2 = 4; //CentBin(3, 4); | |
3861 | else if (k == CentBin(2, 10)) k2 = 6; //CentBin(5, 10); | |
3862 | else if (k == CentBin(1, 3)) k2 = 2; //CentBin(1, 2); | |
3863 | else if (k == CentBin(3, 5)) k2 = 4; //CentBin(3, 4); | |
3864 | else if (k == CentBin(0, 5)) k2 = 3; //CentBin(2, 3); | |
3865 | ||
3866 | else if (k == CentBin(0, 0)) k2 = 1; //pp | |
3867 | else { | |
3868 | Error("MeanPt()", "cent bin %d (%d-%d%%) not assigned", k, centlow[k], centhigh[k]); | |
3869 | } | |
3870 | ||
3871 | const char* m = opt.Contains("squared") ? "Mean2" : "Mean"; | |
3872 | ||
3873 | TString prName_t = Form("h%sPtTrgCen%d", m, k2); | |
3874 | TString prName_a = Form("h%sPtAssCen%d", m, k2); | |
3875 | TProfile2D* hmpt_t = (TProfile2D*)fin->Get(prName_t.Data()); | |
3876 | TProfile2D* hmpt_a = (TProfile2D*)fin->Get(prName_a.Data()); | |
3877 | ||
3878 | if (!hmpt_t) { | |
3879 | Error("MeanPt", "Problem finding TProfile %s. Input %d", prName_t.Data(), k); | |
3880 | gSystem->Exit(-1); | |
3881 | } | |
3882 | // cout << hmpt_t->GetTitle() << " " << centlow[k] << " " << centhigh[k] << endl; | |
3883 | ||
3884 | double mptt = hmpt_t->GetBinContent(i+1, j+1); | |
3885 | double mpta = hmpt_a->GetBinContent(i+1, j+1); | |
3886 | ||
3887 | if (0) { | |
3888 | // x axis is trig pt, y axis assc pt | |
3889 | TAxis* ax = hmpt_t->GetXaxis(); | |
3890 | TAxis* ay = hmpt_t->GetYaxis(); | |
3891 | Printf("trig %.3g - %.3g, assc: %.3g - %.3g, <pt_t,a> %.3g, %.3g", | |
3892 | ax->GetBinLowEdge(i+1), | |
3893 | ax->GetBinUpEdge(i+1), | |
3894 | ay->GetBinLowEdge(j+1), | |
3895 | ay->GetBinUpEdge(j+1), | |
3896 | mptt, mpta ); | |
3897 | } | |
3898 | ||
3899 | if (t_or_a.Contains("t")) | |
3900 | return mptt; | |
3901 | if (t_or_a.Contains("a")) | |
3902 | return mpta; | |
3903 | else | |
3904 | return -1; | |
3905 | } | |
3906 | ||
3907 | ||
3908 | int CentColor(double cen1, double cen2) | |
3909 | { | |
3910 | if (cen1==0 && cen2==10) | |
3911 | return kGray+3; | |
3912 | if (cen1==0 && cen2==20) | |
3913 | return kMagenta+3; | |
3914 | if (cen1==0 && cen2==2) | |
3915 | return kBlack; | |
3916 | if (cen1==2 && cen2==10) | |
3917 | return kRed; | |
3918 | if (cen1==10 && cen2==20) | |
3919 | return kOrange-3; | |
3920 | if (cen1==20 && cen2==30) | |
3921 | return kGreen+1; | |
3922 | if (cen1==30 && cen2==40) | |
3923 | return kYellow+2; | |
3924 | if (cen1==40 && cen2==50) | |
3925 | return kAzure; | |
3926 | if (cen1==50 && cen2==60) | |
3927 | return kBlue-1; | |
3928 | if (cen1==60 && cen2==90) | |
3929 | return kViolet; | |
3930 | ||
3931 | return kBlack; | |
3932 | } | |
3933 | ||
3934 | int CentColorPale(double cen1, double cen2) | |
3935 | { | |
3936 | if (cen1==0 && cen2==10) | |
3937 | return kGray+1; | |
3938 | if (cen1==0 && cen2==20) | |
3939 | return kPink-1; // kPink-6; | |
3940 | if (cen1==0 && cen2==2) | |
3941 | return kGray+1; | |
3942 | if (cen1==2 && cen2==10) | |
3943 | return kRed-9; | |
3944 | if (cen1==10 && cen2==20) | |
3945 | return kOrange-4; | |
3946 | if (cen1==20 && cen2==30) | |
3947 | return kGreen-7; | |
3948 | if (cen1==30 && cen2==40) | |
3949 | return kYellow-8; | |
3950 | if (cen1==40 && cen2==50) | |
3951 | return kAzure-4; | |
3952 | if (cen1==50 && cen2==60) | |
3953 | return kBlue-9; | |
3954 | if (cen1==60 && cen2==90) | |
3955 | return kViolet-9; | |
3956 | ||
3957 | return kBlack; | |
3958 | } | |
3959 | ||
3960 | int PtColor(double p1, double p2) | |
3961 | { | |
3962 | if (p1==0.25 && p2==0.5) | |
3963 | return kGray+1; | |
3964 | if (p1==0.5 && p2==0.75) | |
3965 | return kYellow+3; | |
3966 | if (p1==0.75 && p2==1.0) | |
3967 | return kPink-6; | |
3968 | if (p1==1.0 && p2==1.5) | |
3969 | return kBlack; | |
3970 | if (p1==1.5 && p2==2.0) | |
3971 | return kRed; | |
3972 | if (p1==2.0 && p2==2.5) | |
3973 | return kOrange-3; | |
3974 | if (p1==2.5 && p2==3) | |
3975 | return kGreen+2; | |
3976 | if (p1==3 && p2==4) | |
3977 | return kBlue; | |
3978 | if (p1==4 && p2==5) | |
3979 | return kViolet + 2; | |
3980 | if (p1==5 && p2==6) | |
3981 | return kViolet+10; | |
3982 | if (p1==4 && p2==6) | |
3983 | return kViolet + 2; | |
3984 | if (p1==6 && p2==8) | |
3985 | return kOrange-7; | |
3986 | if (p1==8 && p2==15) | |
3987 | return kOrange+9; | |
3988 | ||
3989 | return kBlack; | |
3990 | } | |
3991 | ||
3992 | int PtColorPale(double p1, double p2) | |
3993 | { | |
3994 | if (p1==0.25 && p2==0.5) | |
3995 | return kGray; | |
3996 | if (p1==0.5 && p2==0.75) | |
3997 | return kYellow+1; | |
3998 | if (p1==0.75 && p2==1.0) | |
3999 | return kPink-4; | |
4000 | if (p1==1.0 && p2==1.5) | |
4001 | return kGray+1; | |
4002 | if (p1==1.5 && p2==2.0) | |
4003 | return kRed-7; | |
4004 | if (p1==2.0 && p2==2.5) | |
4005 | return kOrange-4; | |
4006 | if (p1==2.5 && p2==3) | |
4007 | return kSpring-4; | |
4008 | if (p1==3 && p2==4) | |
4009 | return kAzure+6; | |
4010 | if (p1==4 && p2==5) | |
4011 | return kViolet+1; | |
4012 | if (p1==5 && p2==6) | |
4013 | return kViolet+8; | |
4014 | if (p1==4 && p2==6) | |
4015 | return kViolet+1; | |
4016 | if (p1==6 && p2==8) | |
4017 | return kOrange-9; | |
4018 | if (p1==8 && p2==15) | |
4019 | return kOrange+6; | |
4020 | ||
4021 | return kGray; | |
4022 | } | |
4023 | ||
4024 | int PtaColor(double p1, double p2) | |
4025 | { | |
4026 | if (p1==0.25 && p2==0.5) | |
4027 | return kGray+1; | |
4028 | if (p1==0.5 && p2==0.75) | |
4029 | return kBlack; | |
4030 | if (p1==0.75 && p2==1.0) | |
4031 | return kRed+1; | |
4032 | if (p1==1.0 && p2==1.5) | |
4033 | return kOrange-3;//kBlack; | |
4034 | if (p1==1.5 && p2==2.0) | |
4035 | return kGreen+2;//kRed; | |
4036 | if (p1==2.0 && p2==2.5) | |
4037 | return kAzure-3;//kOrange-3; | |
4038 | if (p1==2.5 && p2==3) | |
4039 | return kViolet + 2;//kBlue; | |
4040 | if (p1==3 && p2==4) | |
4041 | return kViolet + 10;//kBlue; | |
4042 | if (p1==4 && p2==5) | |
4043 | return kMagenta+3; //kViolet + 2; | |
4044 | if (p1==5 && p2==6) | |
4045 | return kViolet; | |
4046 | if (p1==4 && p2==6) | |
4047 | return kMagenta-5; | |
4048 | if (p1==6 && p2==8) | |
4049 | return kGray+2; | |
4050 | if (p1==8 && p2==15) | |
4051 | return kGray+3; | |
4052 | ||
4053 | return kBlack; | |
4054 | } | |
4055 | ||
4056 | int PtaColorPale(double p1, double p2) | |
4057 | { | |
4058 | if (p1==0.25 && p2==0.5) | |
4059 | return kGray; | |
4060 | if (p1==0.5 && p2==0.75) | |
4061 | return kGray+2; | |
4062 | if (p1==0.75 && p2==1.0) | |
4063 | return kPink-2; | |
4064 | if (p1==1.0 && p2==1.5) | |
4065 | return kOrange-4; | |
4066 | if (p1==1.5 && p2==2.0) | |
4067 | return kGreen-3; | |
4068 | if (p1==2.0 && p2==2.5) | |
4069 | return kAzure-4; | |
4070 | if (p1==2.5 && p2==3) | |
4071 | return kViolet; | |
4072 | if (p1==3 && p2==4) | |
4073 | return kViolet + 2; | |
4074 | if (p1==4 && p2==5) | |
4075 | return kMagenta+3; | |
4076 | if (p1==5 && p2==6) | |
4077 | return kViolet+10; | |
4078 | if (p1==4 && p2==6) | |
4079 | return kMagenta-5; | |
4080 | if (p1==6 && p2==8) | |
4081 | return kGray+2; | |
4082 | if (p1==8 && p2==15) | |
4083 | return kGray+3; | |
4084 | ||
4085 | return kGray; | |
4086 | } | |
4087 | ||
4088 | void SaveCanvases(TObjArray* canvases, const char* fileName) | |
4089 | { | |
4090 | TFile* f = new TFile(fileName, "recreate"); | |
4091 | ||
4092 | for (int n=0; n<canvases->GetEntries(); n++) { | |
4093 | TCanvas* c = (TCanvas*)canvases->At(n); | |
4094 | c->Write(c->GetTitle()); | |
4095 | } | |
4096 | f->Close(); | |
4097 | return; | |
4098 | } | |
4099 | ||
4100 | void SaveCanvasesFromFile(const char* rootFile, const char* targetDir, const char* tag, const char* fileType) | |
4101 | { | |
4102 | // Get a list of canvases from rootFile into array, then save each | |
4103 | // to its own file in targetDir/. fileType = "eps", "pdf", "C", | |
4104 | // "png", etc. Not all formats have been tested. | |
4105 | gStyle->SetOptTitle(0); | |
4106 | gStyle->SetOptStat(0); | |
4107 | TString name = ""; | |
4108 | TString base(targetDir); | |
4109 | TFile *cFile = new TFile(rootFile, "read"); | |
4110 | cout << cFile->GetName() << endl; | |
4111 | TObjArray* cList = GetObjectsFromFile(*cFile, "TCanvas"); | |
4112 | ||
4113 | for (int n=0; n<cList->GetEntries(); n++) { | |
4114 | TCanvas* c = (TCanvas*)cList->At(n); | |
4115 | if (c) { | |
4116 | name = ""; | |
4117 | name = base; | |
4118 | name += TString("/"); | |
4119 | name += TString(fileType); | |
4120 | name += TString("/"); | |
4121 | name += TString(c->GetTitle()); | |
4122 | name += TString("."); | |
4123 | name += TString(fileType); | |
4124 | cout<<name.Data()<<endl; | |
4125 | ||
4126 | c->Draw(); | |
4127 | c->Modified(); | |
4128 | c->Update(); | |
4129 | c->SaveAs(name.Data()); | |
4130 | } | |
4131 | else | |
4132 | Error("SaveCanvasesFromFile()", "!c"); | |
4133 | } | |
4134 | ||
4135 | if (1) { | |
4136 | utils.print_pdf(cList, Form("%s/all-figs%s", targetDir, tag), "pdf"); | |
4137 | } | |
4138 | ||
4139 | return; | |
4140 | } | |
4141 | ||
4142 | //void SaveCanvases(TObjArray* canvases, const char* fileName) | |
4143 | //void SaveCanvasesFromFile(const char* rootFile, const char* targetDir, const char* fileType) | |
4144 | TObjArray* GetObjectsFromFile(TFile& file, TString clname, TString dir) | |
4145 | { | |
4146 | file.cd(dir.Data()); | |
4147 | ||
4148 | TObjArray* objList = new TObjArray(); | |
4149 | TIter next(gDirectory->GetListOfKeys()); | |
4150 | TKey *key; | |
4151 | ||
4152 | while ((key=(TKey*)next())) { | |
4153 | TString className(key->GetClassName()); | |
4154 | TString keyName(key->GetName()); | |
4155 | if (1) | |
4156 | printf("%10s %20s\n", className.Data(), keyName.Data()); | |
4157 | ||
4158 | if (className.Contains(clname)) { | |
4159 | objList->Add(gDirectory->Get(keyName.Data())); | |
4160 | } | |
4161 | } | |
4162 | ||
4163 | cout << objList->GetEntries() << " objects retrieved from " | |
4164 | << file.GetName() << "/" << gDirectory->GetName() | |
4165 | << endl; | |
4166 | ||
4167 | return objList; | |
4168 | } |