]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/PHOSTasks/PHOS_PbPb/macros/production/DrawProduction.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPb / macros / production / DrawProduction.C
CommitLineData
f4b109c4 1#include "TCanvas.h"
2#include <TH1.h>
3#include <TF1.h>
4#include <TFile.h>
5#include <Rtypes.h>
6#include <TLegend.h>
7#include <TString.h>
6c3d03d7 8#include <TAttMarker.h>
9#include <RtypesCint.h>
10#include <TNamed.h>
11#include "TDirectoryFile.h"
12#include <TPRegexp.h>
97527d1f 13#include <TGraphAsymmErrors.h>
42bd76a1 14#include <cfloat>
f4b109c4 15
16namespace RawProduction {
17 class Output;
18}
19
6c3d03d7 20enum EffVersion { LHC10h_1234_Apr_10 };
21EffVersion effVersion = LHC10h_1234_Apr_10;
22bool canvHalfWidth=false;
42bd76a1 23int maxFailedCombined = 0;
24
25TH1* MakeCombinedMethodeProduction(const RawProduction::Output& rawOutput, const TString& trigger="kMB", int fromCent=0, int toCent=10, const TString& pid="All", const TString& graphName="yr1", Color_t color=kBlack, Style_t style=kFullDotSmall);
26TH1* MakeCombinedProduction(const RawProduction::Output& rawOutput, const TString& trigger="kMB", int fromCent=0, int toCent=10, const TString& pid="All", const TString& graphName="yr1", Color_t color=kBlack, Style_t style=kFullDotSmall);
6c3d03d7 27
97527d1f 28TF1* GetEfficiency(const TString& trigger, int fromCent, int toCent, const TString& pid, const TString& methode)
f4b109c4 29{
6c3d03d7 30 // Dmitri LHC10h Apr. 10 Efficiencies
31 // avalible pids: Allcore Disp2core CPVcore Both2core Disp2 All CPV Both2
32 if( effVersion == LHC10h_1234_Apr_10 ) {
33 if( ! methode.Contains("yr1") ) {
97527d1f 34 Printf("ERROR:GetEfficiency: only pol1 efficiancies avalible");
6c3d03d7 35 return 0x0;
36 }
37
38 // Cent bin defintion as extracted from PHOS_embedding/
39 int cent = -1;
40 if( 0 == fromCent && 5 == toCent ) cent = 0;
41 if( 5 == fromCent && 10 == toCent ) cent = 1;
42 if( 10 == fromCent && 20 == toCent ) cent = 2;
43 if( 20 == fromCent && 40 == toCent ) cent = 3;
44 if( 40 == fromCent && 60 == toCent ) cent = 4;
45 if( 60 == fromCent && 80 == toCent ) cent = 5;
46 if( cent < 0 ) {
97527d1f 47 Printf("ERROR:GetEfficiency not avalible for centrality [%i,%i)", fromCent, toCent);
6c3d03d7 48 return 0x0;
49 }
50
51 // determine name and return efficiancy function
52 TDirectory* pastDir = gDirectory;
53 TFile* file = TFile::Open("PHOS_eff_Full_PbPb_1234.root", "READ");
54 char funcName[256];
55 if ( methode.Contains("int") )
56 sprintf(funcName, "eff_int_Pi0_Gaus_PbPb_%s_cen%i", pid.Data(), cent);
57 else
58 sprintf(funcName, "eff_int_Pi0_Gaus_PbPb_%s_cen%i", pid.Data(), cent);
59 TF1* func = dynamic_cast<TF1*> ( file->Get(funcName) );
60 pastDir->cd();
61
62 if( ! func )
97527d1f 63 Printf("ERROR:GetEfficiency: efficiancy function %s does not exist", funcName);
6c3d03d7 64 return func;
65 }
66 return 0x0; // should not be reached
f4b109c4 67}
68
6c3d03d7 69TH1* GetRawProduction(const RawProduction::Output& rawOutput, const TString& trigger="kMB", int fromCent=0, int toCent=10,
70 const TString& pid="All", const TString& graphName="yr1", Color_t color=kBlack, Style_t style=kFullDotSmall)
f4b109c4 71{
1c471965 72 TString newName = Form("raw_%s_%02i-%02i_%s_%s", trigger.Data(), fromCent, toCent, pid.Data(), graphName.Data());
73 TH1* hist = dynamic_cast<TH1*> ( gDirectory->Get(newName.Data()) );
74
75 if( ! hist ) {
76 TString oldName(Form("%s/c%02i-%02i/%s/%s", trigger.Data(), fromCent, toCent, pid.Data(), graphName.Data()));
77 TH1* hist = rawOutput.GetHistogram(oldName.Data());
78 hist->SetName(newName.Data());
79 hist->SetTitle(Form("%s, %02i-%02i%%, %s", trigger.Data(), fromCent, toCent, graphName.Data()));
80 hist->GetXaxis()->SetTitle("p_{T}");
81 }
6c3d03d7 82 hist->SetLineColor(color);
83 hist->SetMarkerColor(color);
84 hist->SetMarkerStyle(style);
6c3d03d7 85 return hist;
86}
f4b109c4 87
42bd76a1 88
89
6c3d03d7 90TH1* MakeProduction(const RawProduction::Output& rawOutput, const TString& trigger="kMB", int fromCent=0, int toCent=10,
42bd76a1 91 const TString& pid="All", const TString& methode="yr1", Color_t color=kBlack, Style_t style=kFullDotSmall)
92{
6c3d03d7 93 // First, check if production histogram allready exist in cd.
42bd76a1 94 TString name = Form("prod_%s_%02i-%02i_%s_%s", trigger.Data(), fromCent, toCent, pid.Data(), methode.Data());
6c3d03d7 95 TH1* hist = dynamic_cast<TH1*> ( gDirectory->Get(name.Data()) );
42bd76a1 96
97527d1f 97 // Check if combined, The order corresponds to combining methodes then PID.
42bd76a1 98 if( ! hist ) hist = MakeCombinedProduction(rawOutput, trigger, fromCent, toCent, pid, methode, color, style);
99 if( ! hist ) hist = MakeCombinedMethodeProduction(rawOutput, trigger, fromCent, toCent, pid, methode, color, style);
100
f4b109c4 101
42bd76a1 102 // else clone raw and correct for efficiancy
103 if( ! hist ) {
104 const TH1* rawHist = GetRawProduction(rawOutput, trigger, fromCent, toCent, pid, methode, color, style);
1c471965 105 hist = (TH1*) rawHist->Clone(name.Data());
97527d1f 106 TF1* efficiency = GetEfficiency(trigger, fromCent, toCent, pid, methode);
107 hist->Divide(efficiency, 2*TMath::Pi());
108 hist->GetYaxis()->SetTitle("#frac{d^{2}N_{#pi^{0}}}{2#pi p_{T}dp_{T}dy N_{ev}}");
1c471965 109 }
110 hist->SetLineColor(color);
111 hist->SetMarkerColor(color);
112 hist->SetMarkerStyle(style);
6c3d03d7 113 return hist;
114}
f4b109c4 115
42bd76a1 116TH1* CombinePID(const RawProduction::Output& rawOutput, const TString& trigger, int fromCent, int toCent, const TString& combinedPidName, const TString& pids, const TString& methode, Color_t color, Style_t style)
117{
118 // First, check if production histogram allready exist in cd.
119 TString name = Form("prod_%s_%02i-%02i_%s_%s", trigger.Data(), fromCent, toCent, combinedPidName.Data(), methode.Data());
120 TH1* hist = dynamic_cast<TH1*> ( gDirectory->Get(name.Data()) );
121 if( hist ) return hist;
122
123 const int capacity = 16;
124 const int nHists = TMath::Min( pids.CountChar(' ')+1, capacity);
125 TH1* hists[capacity];
126 TStringToken pidst(pids, " ");
127 pidst.NextToken();
128 const TString firstPID = pidst;
129 hists[0] = MakeProduction(rawOutput, trigger, fromCent, toCent, firstPID, methode, color, style);
130 int index = 1;
131 while( pidst.NextToken() && index < capacity ) {
132 hists[index] = MakeProduction(rawOutput, trigger, fromCent, toCent, pidst, methode, color, style);
133 ++index;
134 }
135 TH1* combHist = hists[0]->Clone(name.Data());
136 combHist->SetTitle(TString(combHist->GetTitle()).ReplaceAll(firstPID.Data(), combinedPidName.Data()));
137
138 for(int ptBin=0; ptBin<combHist->GetNbinsX(); ++ptBin){
139 double wmeansum = 0.;
140 double sumw = 0.;
141 double mins = DBL_MAX;
142 int nFailed = 0;
143 for(int i=0; i<nHists; ++i) {
144 const double x = hists[i]->GetBinContent(ptBin);
145 const double s = hists[i]->GetBinError(ptBin);
146 if( 0.==x || 0.==s) {
147 nFailed++;
148 continue;
149 }
150 const double w = 1./(s*s);
151 wmeansum += w*x;
152 sumw +=w;
153 if( mins > s )
154 mins = s;
155 }
156 if(nFailed > maxFailedCombined) {
157 combHist->SetBinContent(ptBin, 0.);
158 combHist->SetBinError(ptBin, 0.);
159 }
160 else {
161 const double wmean = wmeansum/sumw;
162 combHist->SetBinContent(ptBin, wmean);
163 combHist->SetBinError(ptBin, mins);
164 }
165 }
166
167 return combHist;
168}
169
170
171TH1* MakeCombinedProduction(const RawProduction::Output& rawOutput, const TString& trigger, int fromCent, int toCent, const TString& combinedPidName, const TString& methode, Color_t color, Style_t style)
172{
173 TString pids;
174 if( combinedPidName.Contains("Combcore") )
175 pids = "Allcore CPVcore Disp2core Both2core";
176 else
177 return 0x0;
178
179 // First, check if production histogram allready exist in cd.
180 TString name = Form("prod_%s_%02i-%02i_%s_%s", trigger.Data(), fromCent, toCent, combinedPidName.Data(), methode.Data());
181 TH1* hist = dynamic_cast<TH1*> ( gDirectory->Get(name.Data()) );
182 if( hist )
183 return hist;
184
185 return CombinePID(rawOutput, trigger, fromCent, toCent, combinedPidName, pids, methode, color, style);
186}
187
188TH1* MergeMethodes(const RawProduction::Output& rawOutput, const TString& trigger, int fromCent, int toCent, const TString& pid, const TString& combMethodeName, const TString& methodes, Color_t color, Style_t style)
189{
190 // First, check if production histogram allready exist in cd.
191 TString name = Form("prod_%s_%02i-%02i_%s_%s", trigger.Data(), fromCent, toCent, pid.Data(), combMethodeName.Data());
192 TH1* hist = dynamic_cast<TH1*> ( gDirectory->Get(name.Data()) );
193 if( hist ) return hist;
194
195 const int capacity = 16;
196 const int nHists = TMath::Min( methodes.CountChar(' ')+1, capacity);
197 TH1* hists[capacity];
198 TStringToken mst(methodes, " ");
199 mst.NextToken();
200 const TString firstMethode = mst;
201 hists[0] = MakeProduction(rawOutput, trigger, fromCent, toCent, pid, firstMethode, color, style);
202 int index = 1;
203 while( mst.NextToken() && index < capacity ) {
204 hists[index] = MakeProduction(rawOutput, trigger, fromCent, toCent, pid, mst, color, style);
205 ++index;
206 }
207 TH1* combHist = hists[0]->Clone(name.Data());
208 combHist->SetTitle(TString(combHist->GetTitle()).ReplaceAll(firstMethode.Data(),combMethodeName.Data()));
209
210 for(int ptBin=0; ptBin<combHist->GetNbinsX(); ++ptBin){
211 double wmeansum = 0.;
212 double sumw = 0.;
213 double mins = DBL_MAX;
214 int nFailed = 0;
215 for(int i=0; i<nHists; ++i) {
216 const double x = hists[i]->GetBinContent(ptBin);
217 const double s = hists[i]->GetBinError(ptBin);
218 if( 0.==x || 0.==s) {
219 nFailed++;
220 continue;
221 }
222 const double w = 1./(s*s);
223 wmeansum += w*x;
224 sumw +=w;
225 if( mins > s )
226 mins = s;
227 }
228 if(nFailed > maxFailedCombined) {
229 combHist->SetBinContent(ptBin, 0.);
230 combHist->SetBinError(ptBin, 0.);
231 }
232 else {
233 const double wmean = wmeansum/sumw;
234 combHist->SetBinContent(ptBin, wmean);
235 combHist->SetBinError(ptBin, mins);
236 }
237 }
238
239 return combHist;
240}
241
242TH1* MakeCombinedMethodeProduction(const RawProduction::Output& rawOutput, const TString& trigger, int fromCent, int toCent, const TString& pid, const TString& combMethodeName, Color_t color, Style_t style)
243{
244 TString methodes;
245 if( combMethodeName.EqualTo("yr1comb") )
246 methodes = "yr1 yr1int";
247 else
248 return 0x0;
249
250 // First, check if production histogram allready exist in cd.
251 TString name = Form("prod_%s_%02i-%02i_%s_%s", trigger.Data(), fromCent, toCent, pid.Data(), combMethodeName.Data());
252 TH1* hist = dynamic_cast<TH1*> ( gDirectory->Get(name.Data()) );
253 if( hist )
254 return hist;
255
256 return MergeMethodes(rawOutput, trigger, fromCent, toCent, pid, combMethodeName, methodes, color, style);
257}
258
6c3d03d7 259TH1* MakeRatio(TH1* h1, TH1* h2, const TString& title ="")
260{
261 TString name = Form("%s_%s", h1->GetName(), h2->GetName());
262 TH1* hist = dynamic_cast<TH1*> ( gDirectory->Get(name.Data()) );
263 if( hist )
264 return hist;
f4b109c4 265
6c3d03d7 266 hist = (TH1*)h1->Clone(name.Data());
42bd76a1 267 hist->Divide(h1, h2, 1, 1, "B");
6c3d03d7 268 hist->SetTitle(title.Data());
269 hist->GetYaxis()->SetTitle("Ratio");
270 return hist;
271}
f4b109c4 272
f4b109c4 273
6c3d03d7 274TCanvas* DrawPIDProductionWithRatios(const RawProduction::Output& rawOutput, const TString& trigger,
275 int fromCent, int toCent, const TString& methode="yr1",
276 const TString& pids = TString("Allcore CPVcore Disp2core Both2core"), bool raw = false)
277{
278 const int capacity = 8;
279 const int nHists = TMath::Min( pids.CountChar(' ')+1, capacity);
280 TStringToken pidst = TStringToken(pids, " ");
281 char pidsa[capacity][64] ={""};
f4b109c4 282
6c3d03d7 283 TH1* hists[capacity] = {0x0};
284 const Style_t markers[capacity] = {22, 22, 23, 33, 24, 26, 32, 27};
285 const Color_t colors[capacity] = {kBlack, kRed, kBlue, kGreen, kGray, kMagenta, kCyan, kOrange};
f4b109c4 286
6c3d03d7 287 int index = 0;
288 while(pidst.NextToken() && index < capacity) {
289 sprintf(pidsa[index], "%s", pidst.Data());
290 if(raw)
291 hists[index] = GetRawProduction(rawOutput, trigger, fromCent, toCent, pidst, methode, colors[index], markers[index]);
292 else
293 hists[index] = MakeProduction(rawOutput, trigger, fromCent, toCent, pidst, methode, colors[index], markers[index]);
294 index++;
295 }
f4b109c4 296
6c3d03d7 297 //TString pids_underscore = TString(pids); pids_underscore.ReplaceAll(" ", "_");
298 TString key = Form("PIDRatios_%s_c%02i-%02i_%s_%s_raw%i_h%i", trigger.Data(), fromCent, toCent, methode.Data(), TString(pids).ReplaceAll(" ", "_").Data(), raw, canvHalfWidth);
299 //TString key = Form("PIDRatios_%s_c%02i-%02i_%s_raw%i", trigger.Data(), fromCent, toCent, methode.Data(), raw);
f4b109c4 300
6c3d03d7 301 TCanvas* canv = 0x0;
302 if( canvHalfWidth)
303 canv = new TCanvas(key.Data(), key.Data(), 1024/2, 768);
304 else
305 canv = new TCanvas(key.Data(), key.Data(), 1024, 768);
306
307 // Direct
308 canv->cd();
309 TPad *pad1 = new TPad("pad1","pad1",0,0.4,1,1);
310 pad1->SetBottomMargin(0);
311 pad1->SetLogy();
312 pad1->Draw();
313 pad1->cd();
314 if( canvHalfWidth ) {
315 hists[0]->GetYaxis()->SetLabelFont(63); //font in pixels
316 hists[0]->GetYaxis()->SetLabelSize(10); //in pixels
317 hists[0]->GetYaxis()->SetTitleSize(0.03);
318 hists[0]->GetYaxis()->SetTitleOffset(1.2);
319 } else {
320 hists[0]->GetYaxis()->SetLabelFont(63); //font in pixels
321 hists[0]->GetYaxis()->SetLabelSize(20); //in pixels
322 hists[0]->GetYaxis()->SetTitleSize(0.055);
323 hists[0]->GetYaxis()->SetTitleOffset(0.7);
f4b109c4 324 }
6c3d03d7 325 if( raw )
326 hists[0]->GetYaxis()->SetRangeUser(1.e-7, 1.e1);
327 else
97527d1f 328 hists[0]->GetYaxis()->SetRangeUser(1.e-6, 1.e3);
6c3d03d7 329 // Draw
330 hists[0]->DrawCopy();
331 for(int i=1;i<nHists;++i)
332 hists[i]->DrawCopy("same");
333 // Legend
334 TLegend* leg1 = new TLegend(0.7,0.6,0.85,0.88);
335 for(int i=0;i<nHists;++i)
336 leg1->AddEntry(hists[i], pidsa[i] , "lep");
337 leg1->Draw();
338
339 // Ratios
340 canv->cd();
341 TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.4);
342 pad2->SetTopMargin(0);
343 pad2->Draw();
344 pad2->cd();
345 pad2->SetGridy();
346 TH1* firstRatio = MakeRatio(hists[1], hists[0]);
347 if ( raw )
348 firstRatio->GetYaxis()->SetRangeUser(0., 2.);
349 else
350 firstRatio->GetYaxis()->SetRangeUser(0.6, 1.4);
351 firstRatio->GetYaxis()->SetTitle("");
352 firstRatio->GetYaxis()->SetLabelFont(63);
353 firstRatio->GetYaxis()->SetLabelSize(25);
354 firstRatio->GetXaxis()->SetLabelFont(63);
355 firstRatio->GetXaxis()->SetLabelSize(20);
356 firstRatio->SetTitle("Ratio");
357 // Draw
358 firstRatio->DrawCopy("AXIS");
359 firstRatio->DrawCopy("AXIGsame");
360 firstRatio->DrawCopy("same");
361 for(int i=2;i<nHists;++i)
362 MakeRatio(hists[i],hists[0])->DrawCopy("same");
363 // Ratios
364 TLegend* leg2 = new TLegend(0.7,0.63,0.85,0.98);
365 for(int i=1;i<nHists;++i)
366 leg2->AddEntry(MakeRatio(hists[i],hists[0]), Form("%s/%s", pidsa[i], pidsa[0]), "lep");
367 leg2->Draw();
368
369 canv->SaveAs(Form("imgs/%s.pdf", key.Data()));
370 canv->SaveAs(Form("imgs/%s.png", key.Data()));
371
372 return canv;
373 //delete canv;
f4b109c4 374}
375
376
42bd76a1 377TCanvas* DrawMethodeProductionWithRatios(const RawProduction::Output& rawOutput, const TString& trigger,
378 int fromCent, int toCent, const TString& methodes="yr1comb yr1 yr1int",
379 const TString& pid = "Allcore", bool raw = false)
380{
381 const int capacity = 8;
382 const int nHists = TMath::Min( methodes.CountChar(' ')+1, capacity);
383 TStringToken methodest = TStringToken(methodes, " ");
384 char methodesa[capacity][64] ={""};
385
386 TH1* hists[capacity] = {0x0};
387 const Style_t markers[capacity] = {22, 22, 23, 33, 24, 26, 32, 27};
388 const Color_t colors[capacity] = {kBlack, kRed, kBlue, kGreen, kGray, kMagenta, kCyan, kOrange};
389
390 int index = 0;
391 while(methodest.NextToken() && index < capacity) {
392 sprintf(methodesa[index], "%s", methodest.Data());
393 if(raw)
394 hists[index] = GetRawProduction(rawOutput, trigger, fromCent, toCent, pid, methodest, colors[index], markers[index]);
395 else
396 hists[index] = MakeProduction(rawOutput, trigger, fromCent, toCent, pid, methodest, colors[index], markers[index]);
397 index++;
398 }
399
400 TString key = Form("MethodeRatios_%s_c%02i-%02i_%s_%s_raw%i_h%i", trigger.Data(), fromCent, toCent, TString(methodes).ReplaceAll(" ", "_").Data(), pid.Data(), raw, canvHalfWidth);
401
402 TCanvas* canv = 0x0;
403 if( canvHalfWidth)
404 canv = new TCanvas(key.Data(), key.Data(), 1024/2, 768);
405 else
406 canv = new TCanvas(key.Data(), key.Data(), 1024, 768);
407
408 // Direct
409 canv->cd();
410 TPad *pad1 = new TPad("pad1","pad1",0,0.4,1,1);
411 pad1->SetBottomMargin(0);
412 pad1->SetLogy();
413 pad1->Draw();
414 pad1->cd();
415 if( canvHalfWidth ) {
416 hists[0]->GetYaxis()->SetLabelFont(63); //font in pixels
417 hists[0]->GetYaxis()->SetLabelSize(10); //in pixels
418 hists[0]->GetYaxis()->SetTitleSize(0.03);
419 hists[0]->GetYaxis()->SetTitleOffset(1.2);
420 } else {
421 hists[0]->GetYaxis()->SetLabelFont(63); //font in pixels
422 hists[0]->GetYaxis()->SetLabelSize(20); //in pixels
423 hists[0]->GetYaxis()->SetTitleSize(0.055);
424 hists[0]->GetYaxis()->SetTitleOffset(0.7);
425 }
426 if( raw )
427 hists[0]->GetYaxis()->SetRangeUser(1.e-7, 1.e1);
428 else
97527d1f 429 hists[0]->GetYaxis()->SetRangeUser(1.e-6, 1.e3);
42bd76a1 430 // Draw
431 hists[0]->DrawCopy();
432 for(int i=1;i<nHists;++i)
433 hists[i]->DrawCopy("same");
434 // Legend
435 TLegend* leg1 = new TLegend(0.7,0.6,0.85,0.88);
436 for(int i=0;i<nHists;++i)
437 leg1->AddEntry(hists[i], methodesa[i] , "lep");
438 leg1->Draw();
439
440 // Ratios
441 canv->cd();
442 TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.4);
443 pad2->SetTopMargin(0);
444 pad2->Draw();
445 pad2->cd();
446 pad2->SetGridy();
447 TH1* firstRatio = MakeRatio(hists[1], hists[0]);
448 if ( raw )
449 firstRatio->GetYaxis()->SetRangeUser(0., 2.);
450 else
451 firstRatio->GetYaxis()->SetRangeUser(0.6, 1.4);
452 firstRatio->GetYaxis()->SetTitle("");
453 firstRatio->GetYaxis()->SetLabelFont(63);
454 firstRatio->GetYaxis()->SetLabelSize(25);
455 firstRatio->GetXaxis()->SetLabelFont(63);
456 firstRatio->GetXaxis()->SetLabelSize(20);
457 firstRatio->SetTitle("Ratio");
458 // Draw
459 firstRatio->DrawCopy("AXIS");
460 firstRatio->DrawCopy("AXIGsame");
461 firstRatio->DrawCopy("same");
462 for(int i=2;i<nHists;++i)
463 MakeRatio(hists[i],hists[0])->DrawCopy("same");
464 // Ratios
465 TLegend* leg2 = new TLegend(0.7,0.63,0.85,0.98);
466 for(int i=1;i<nHists;++i)
467 leg2->AddEntry(MakeRatio(hists[i],hists[0]), Form("%s/%s", methodesa[i], methodesa[0]), "lep");
468 leg2->Draw();
469
470 canv->SaveAs(Form("imgs/%s.pdf", key.Data()));
471 canv->SaveAs(Form("imgs/%s.png", key.Data()));
472
473 return canv;
474 //delete canv;
475}
476
97527d1f 477TGraphAsymmErrors* GetPWWGACombinedYield(int fromCent, int toCent, const TString& detector = "")
478{
479 // detector may be "", "PHOS, "PCM"
480
481 // InvYieldPbPbStatErr_4060
482 // InvYieldPbPbPHOSStatErr_0010
483 const TString graphName = Form("InvYieldPbPb%sStatErr_%02i%02i", detector.Data(), fromCent, toCent);
484
485 const TString fileName = "CombinedResultsPbPb.root";
486 TFile* file = TFile::Open(fileName.Data());
487 TGraphAsymmErrors* graph = dynamic_cast<TGraphAsymmErrors*> ( file->Get(graphName.Data()) );
488
489 if( ! graph )
490 Printf("ERROR: GetPWWGACombinedYield(%i,%i,%s): %s was not found", fromCent, toCent, detector.Data(), graphName.Data());
491
492 return graph;
493}
494
495
496void CompareToPWGGA(const RawProduction::Output& rawOutput, const TString& trigger, int fromCent, int toCent,
497 const TString& methode="yr1comb", const TString& pid = "Combcore", const TString& detector="")
498{
499 TGraphAsymmErrors* pwgGraph = GetPWWGACombinedYield(fromCent, toCent, detector);
500 TH1* production = MakeProduction(rawOutput, trigger, fromCent, toCent, pid, methode, kRed, kFullDotMedium);
501 TString name = Form("pwggaRatio%s_%s", detector.Data(), production->GetName());
502 TH1* hRatio = production->Clone(name.Data());
503 TString detStr = detector;
504 if( detStr.EqualTo("") ) detStr = "Combined PWGGA (PHOS+PCM)";
505 hRatio->SetTitle(Form("%s vs %s", hRatio->GetTitle(), detStr.Data()));
506 hRatio->GetYaxis()->SetTitle(Form("PHOS 11h (10h eff.) / %s", detStr.Data()));
507
508
509 TString canvName = Form("PWGGAprod%s_%s", detector.Data(), production->GetName());
510 TCanvas* canv = new TCanvas(canvName.Data(), canvName.Data(), 1024, 768);
511 production->GetYaxis()->SetRangeUser(1.e-6, 1.e3);
512 canv->SetLogy();
513 production->Draw();
514 pwgGraph->Draw();
515 TLegend* leg = new TLegend(0.4,0.7,0.85,0.88);
516 leg->AddEntry(production, "11h (10h eff.)", "lep");
517 leg->AddEntry(pwgGraph, Form("10h %s", detStr.Data()), "lep");
518 leg->Draw();
519 canv->SaveAs(Form("imgs/%s.png", canvName.Data()));
520 canv->SaveAs(Form("imgs/%s.pdf", canvName.Data()));
521
522
523 for(int ptBin=1; ptBin <production->GetNbinsX(); ptBin++) {
524 int graphBin = ptBin-3;
525 if( graphBin >= 0 && pwgGraph->GetN() > ptBin ) {
526 //Printf("%f %f %f", production->GetBinLowEdge(ptBin), pwgGraph->GetX()[graphBin], production->GetBinLowEdge(ptBin+1));
527 double p = production->GetBinContent(ptBin);
528 double p_e = production->GetBinError(ptBin);
529 double pwgP = pwgGraph->GetY()[graphBin];
530 double pwgP_e = (pwgGraph->GetEYlow()[graphBin] + pwgGraph->GetEYhigh()[graphBin])/2;
531 if( 0. != p && 0. != pwgP) {
532 double ratio = p/pwgP;
533 double ratio_e = TMath::Sqrt((p_e*p_e + pwgP_e*pwgP_e*p/pwgP)/pwgP);
534 //Printf("%f %f - %f", p, pwgP, ratio);
535 hRatio->SetBinContent(ptBin, ratio);
536 hRatio->SetBinError(ptBin, ratio_e);
537 }
538 else {
539 hRatio->SetBinContent(ptBin, 0.);
540 hRatio->SetBinError(ptBin, 0.);
541 }
542 }
543 else {
544 hRatio->SetBinContent(ptBin, 0.);
545 hRatio->SetBinError(ptBin, 0.);
546 }
547 }
548 TString canvRName = Form("PWGGAratio%s_%s", detector.Data(), production->GetName());
549 TCanvas* canvRatio = new TCanvas(canvRName.Data(), canvRName.Data(), 1024, 768);;
550 canvRatio->SetLogy();
551 hRatio->GetYaxis()->SetRangeUser(0.1, 5.);
552 hRatio->Draw();
553 canvRatio->SaveAs(Form("imgs/%s.png", canvRName.Data()));
554 canvRatio->SaveAs(Form("imgs/%s.pdf", canvRName.Data()));
555}
42bd76a1 556
6c3d03d7 557void DrawPIDRatios(const RawProduction::Output& rawOutput)
f4b109c4 558{
6c3d03d7 559 const int nCent = 2;
560 int centBins[nCent][2] = {{0,5}, {5,10}/*, {10,20}, {20,40}, {40,60}, {60,80}*/};
561 TStringToken methodes("yr1 yr1int", " ");
562 while( methodes.NextToken() ) {
563 for(int ic=0; ic<nCent; ++ic) {
1c471965 564 DrawPIDProductionWithRatios(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], methodes.Data(), "Allcore CPVcore Disp2core Both2core All CPV Disp2 Both2", true );
565 DrawPIDProductionWithRatios(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], methodes.Data(), "Allcore CPVcore Disp2core Both2core All CPV Disp2 Both2", false );
6c3d03d7 566 DrawPIDProductionWithRatios(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], methodes.Data(), "Allcore CPVcore Disp2core Both2core", true );
567 DrawPIDProductionWithRatios(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], methodes.Data(), "Allcore CPVcore Disp2core Both2core", false );
1c471965 568 DrawPIDProductionWithRatios(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], methodes.Data(), "All Disp2 CPV Both2 Allcore", true );
569 DrawPIDProductionWithRatios(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], methodes.Data(), "All Disp2 CPV Both2 Allcore", false );
6c3d03d7 570 }
571 }
f4b109c4 572}
573
42bd76a1 574void DrawPIDRatiosCombined(const RawProduction::Output& rawOutput)
575{
576 const int nCent = 2;
577 int centBins[nCent][2] = {{0,5}, {5,10}/*, {10,20}, {20,40}, {40,60}, {60,80}*/};
578 TStringToken methodes("yr1comb", " ");
579 while( methodes.NextToken() ) {
580 for(int ic=0; ic<nCent; ++ic) {
581 DrawPIDProductionWithRatios(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], methodes.Data(), "Allcore CPVcore Disp2core Both2core All CPV Disp2 Both2", false );
582 DrawPIDProductionWithRatios(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], methodes.Data(), "Allcore CPVcore Disp2core Both2core", false );
583 DrawPIDProductionWithRatios(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], methodes.Data(), "Combcore Allcore CPVcore Disp2core Both2core", false );
584 DrawPIDProductionWithRatios(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], methodes.Data(), "All Disp2 CPV Both2 Allcore", false );
585 }
586 }
587}
588
589void DrawMethodeRatios(const RawProduction::Output& rawOutput)
590{
591 const int nCent = 2;
592 int centBins[nCent][2] = {{0,5}, {5,10}/*, {10,20}, {20,40}, {40,60}, {60,80}*/};
593 TStringToken pids("Allcore CPVcore Disp2core Both2core", " ");
594 while( pids.NextToken() ) {
595 for(int ic=0; ic<nCent; ++ic) {
596 DrawMethodeProductionWithRatios(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], "yr1comb yr1 yr1int", pids, false );
597 }
598 }
599}
600
97527d1f 601void ComparePWGGA(const RawProduction::Output& rawOutput)
602{
603 const int nCent = 2;
604 int centBins[nCent][2] = {{0,5}, {5,10}/*, {10,20}, {20,40}, {40,60}, {60,80}*/};
605 for(int ic=0; ic<nCent; ++ic) {// TODO: return to "<nCent"
606 CompareToPWGGA(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], "yr1comb", "Combcore");
607 CompareToPWGGA(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], "yr1comb", "Combcore", "PHOS");
608 CompareToPWGGA(rawOutput, "kCentral", centBins[ic][0], centBins[ic][1], "yr1comb", "Combcore", "PCM");
609 }
610}
611
612
f4b109c4 613void DrawProduction()
614{
615 gROOT->LoadMacro("MakeRawProduction.C+g");
6c3d03d7 616 RawProduction::Output rawOutput;
f4b109c4 617 gStyle->SetOptStat(0);
618
97527d1f 619 ComparePWGGA(rawOutput);
620
42bd76a1 621 DrawMethodeRatios(rawOutput);
622 DrawPIDRatiosCombined(rawOutput);
623
624 canvHalfWidth = true;
6c3d03d7 625 DrawPIDRatios(rawOutput);
f4b109c4 626}