]>
Commit | Line | Data |
---|---|---|
e1f47419 | 1 | /** |
2 | * Script to visualise the dN/deta | |
3 | * | |
4 | * This script is independent of any AliROOT code - and should stay that way. | |
5 | */ | |
b2e7f2d6 | 6 | #include <TH1.h> |
7 | #include <THStack.h> | |
8 | #include <TGraphErrors.h> | |
9 | #include <TGraphAsymmErrors.h> | |
10 | #include <TMultiGraph.h> | |
11 | #include <TFile.h> | |
12 | #include <TList.h> | |
13 | #include <TString.h> | |
14 | #include <TError.h> | |
15 | #include <TSystem.h> | |
16 | #include <TROOT.h> | |
17 | #include <TMath.h> | |
18 | #include <TCanvas.h> | |
19 | #include <TPad.h> | |
20 | #include <TStyle.h> | |
21 | #include <TLegend.h> | |
22 | #include <TLegendEntry.h> | |
23 | #include <TLatex.h> | |
24 | #include <TImage.h> | |
25 | ||
2f92ff0e | 26 | /** |
27 | * Class to draw dN/deta results | |
28 | * | |
29 | */ | |
b2e7f2d6 | 30 | struct dNdetaDrawer |
31 | { | |
2f92ff0e | 32 | /** |
33 | * POD of data for range zooming | |
34 | */ | |
b2e7f2d6 | 35 | struct RangeParam |
36 | { | |
37 | TAxis* fMasterAxis; // Master axis | |
38 | TAxis* fSlave1Axis; // First slave axis | |
39 | TVirtualPad* fSlave1Pad; // First slave pad | |
40 | TAxis* fSlave2Axis; // Second slave axis | |
41 | TVirtualPad* fSlave2Pad; // Second slave pad | |
42 | }; | |
43 | //__________________________________________________________________ | |
2f92ff0e | 44 | /** |
45 | * Constructor | |
46 | * | |
47 | */ | |
b2e7f2d6 | 48 | dNdetaDrawer() |
49 | : fShowOthers(false), // Bool_t | |
50 | fShowAlice(false), // Bool_t | |
51 | fShowRatios(false), // Bool_t | |
52 | fShowLeftRight(false), // Bool_t | |
53 | fRebin(5), // UShort_t | |
54 | fCutEdges(false), // Bool_t | |
55 | fTitle(""), // TString | |
b2e7f2d6 | 56 | fTrigString(0), // TNamed* |
57 | fSNNString(0), // TNamed* | |
58 | fSysString(0), // TNamed* | |
59 | fVtxAxis(0), // TAxis* | |
e1f47419 | 60 | // fForward(0), // TH1* |
61 | // fForwardMC(0), // TH1* | |
62 | // fForwardHHD(0), // TH1* | |
63 | // fTruth(0), // TH1* | |
64 | // fCentral(0), // TH1* | |
65 | // fForwardSym(0), // TH1* | |
66 | // fForwardMCSym(0), // TH1* | |
67 | // fForwardHHDSym(0), // TH1* | |
b2e7f2d6 | 68 | fTriggers(0), // TH1* |
69 | fRangeParam(0) | |
70 | ||
71 | { | |
72 | fRangeParam = new RangeParam; | |
73 | fRangeParam->fMasterAxis = 0; | |
74 | fRangeParam->fSlave1Axis = 0; | |
75 | fRangeParam->fSlave1Pad = 0; | |
76 | fRangeParam->fSlave2Axis = 0; | |
77 | fRangeParam->fSlave2Pad = 0; | |
78 | } | |
2f92ff0e | 79 | //================================================================== |
80 | /** | |
81 | * @{ | |
82 | * @name Set parameters | |
83 | */ | |
84 | /** | |
85 | * Show other (UA5, CMS, ...) data | |
86 | * | |
87 | * @param x Whether to show or not | |
88 | */ | |
b2e7f2d6 | 89 | void SetShowOthers(Bool_t x) { fShowOthers = x; } |
2f92ff0e | 90 | //__________________________________________________________________ |
91 | /** | |
92 | * Show ALICE published data | |
93 | * | |
94 | * @param x Wheter to show or not | |
95 | */ | |
b2e7f2d6 | 96 | void SetShowAlice(Bool_t x) { fShowAlice = x; } |
2f92ff0e | 97 | //__________________________________________________________________ |
98 | /** | |
99 | * Whether to show ratios or not. If there's nothing to compare to, | |
100 | * the ratio panel will be implicitly disabled | |
101 | * | |
102 | * @param x Whether to show or not | |
103 | */ | |
b2e7f2d6 | 104 | void SetShowRatios(Bool_t x) { fShowRatios = x; } |
2f92ff0e | 105 | //__________________________________________________________________ |
106 | /** | |
107 | * | |
108 | * Whether to show the left/right asymmetry | |
109 | * | |
110 | * @param x To show or not | |
111 | */ | |
b2e7f2d6 | 112 | void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; } |
2f92ff0e | 113 | //__________________________________________________________________ |
114 | /** | |
115 | * Set the rebinning factor | |
116 | * | |
117 | * @param x Rebinning factor (must be a divisor in the number of bins) | |
118 | */ | |
b2e7f2d6 | 119 | void SetRebin(UShort_t x) { fRebin = x; } |
2f92ff0e | 120 | //__________________________________________________________________ |
121 | /** | |
122 | * Wheter to cut away the edges | |
123 | * | |
124 | * @param x Whether or not to cut away edges | |
125 | */ | |
b2e7f2d6 | 126 | void SetCutEdges(Bool_t x) { fCutEdges = x; } |
2f92ff0e | 127 | //__________________________________________________________________ |
128 | /** | |
129 | * Set the title of the plot | |
130 | * | |
131 | * @param x Title | |
132 | */ | |
b2e7f2d6 | 133 | void SetTitle(TString x) { fTitle = x; } |
2f92ff0e | 134 | /* @} */ |
135 | //================================================================== | |
136 | /** | |
137 | * @{ | |
138 | * @name Override settings from input | |
139 | */ | |
140 | /** | |
141 | * Override setting from file | |
142 | * | |
143 | * @param sNN Center of mass energy per nucleon pair (GeV) | |
144 | */ | |
145 | void SetSNN(UShort_t sNN) | |
146 | { | |
cdb9e20f | 147 | fSNNString = new TNamed("sNN", Form("%04dGeV", sNN)); |
2f92ff0e | 148 | fSNNString->SetUniqueID(sNN); |
149 | } | |
150 | //__________________________________________________________________ | |
151 | /** | |
152 | * Set the collision system | |
153 | * - 1: pp | |
154 | * - 2: PbPb | |
155 | * | |
156 | * @param sys collision system | |
157 | */ | |
158 | void SetSys(UShort_t sys) | |
159 | { | |
fe52e455 | 160 | fSysString = new TNamed("sys", (sys == 1 ? "pp" : |
2f92ff0e | 161 | sys == 2 ? "PbPb" : "unknown")); |
fe52e455 | 162 | fSysString->SetUniqueID(sys); |
2f92ff0e | 163 | } |
164 | //__________________________________________________________________ | |
165 | /** | |
166 | * Set the vertex range in centimeters | |
167 | * | |
168 | * @param vzMin Min @f$ v_z@f$ | |
169 | * @param vzMax Max @f$ v_z@f$ | |
170 | */ | |
171 | void SetVertexRange(Double_t vzMin, Double_t vzMax) | |
172 | { | |
173 | fVtxAxis = new TAxis(10, vzMin, vzMax); | |
174 | fVtxAxis->SetName("vtxAxis"); | |
175 | fVtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", vzMin, vzMax)); | |
176 | } | |
b2e7f2d6 | 177 | //__________________________________________________________________ |
2f92ff0e | 178 | void SetTrigger(UShort_t trig) |
179 | { | |
180 | fTrigString = new TNamed("trigString", (trig & 0x1 ? "INEL" : | |
181 | trig & 0x2 ? "INEL>0" : | |
182 | trig & 0x4 ? "NSD" : | |
183 | "unknown")); | |
184 | fTrigString->SetUniqueID(trig); | |
185 | } | |
186 | ||
187 | ||
188 | //================================================================== | |
189 | /** | |
190 | * @{ | |
191 | * @name Main steering functions | |
192 | */ | |
193 | /** | |
194 | * Run the code to produce the final result. | |
195 | * | |
196 | * @param filename File containing the data | |
197 | */ | |
198 | void Run(const char* filename="forward_dndeta.C") | |
b2e7f2d6 | 199 | { |
b2e7f2d6 | 200 | Double_t max = 0; |
201 | ||
e1f47419 | 202 | if (!Open(filename, max)) return; |
203 | Info("Run", "Got data from %s with maximum %f", filename, max); | |
204 | ||
b2e7f2d6 | 205 | // Create our stack of results |
e1f47419 | 206 | THStack* results = fResults; // StackResults(max); |
b2e7f2d6 | 207 | |
208 | // Create our stack of other results | |
209 | TMultiGraph* other = 0; | |
210 | if (fShowOthers || fShowAlice) other = StackOther(max); | |
211 | ||
212 | Double_t smax = 0; | |
213 | THStack* ratios = 0; | |
214 | if (fShowRatios) ratios = StackRatios(other, smax); | |
215 | ||
216 | Double_t amax = 0; | |
217 | THStack* leftright = 0; | |
218 | if (fShowLeftRight) leftright = StackLeftRight(amax); | |
219 | ||
220 | Plot(results, other, max, ratios, smax, leftright, amax); | |
221 | } | |
e1f47419 | 222 | |
b2e7f2d6 | 223 | //__________________________________________________________________ |
2f92ff0e | 224 | /** |
e1f47419 | 225 | * Fetch the information on the run from the results list |
2f92ff0e | 226 | * |
e1f47419 | 227 | * @param results Results list |
2f92ff0e | 228 | */ |
e1f47419 | 229 | void FetchInformation(const TList* results) |
b2e7f2d6 | 230 | { |
2f92ff0e | 231 | if (!fTrigString) |
232 | fTrigString = static_cast<TNamed*>(results->FindObject("trigString")); | |
233 | if (!fSNNString) | |
234 | fSNNString = static_cast<TNamed*>(results->FindObject("sNN")); | |
235 | if (!fSysString) | |
236 | fSysString = static_cast<TNamed*>(results->FindObject("sys")); | |
237 | if (!fVtxAxis) | |
238 | fVtxAxis = static_cast<TAxis*>(results->FindObject("vtxAxis")); | |
bad9a3c1 | 239 | if (!fTrigString) fTrigString = new TNamed("trigString", "unknown"); |
240 | if (!fSNNString) fSNNString = new TNamed("sNN", "unknown"); | |
241 | if (!fSysString) fSysString = new TNamed("sys", "unknown"); | |
242 | if (!fVtxAxis) { | |
243 | fVtxAxis = new TAxis(1,0,0); | |
244 | fVtxAxis->SetName("vtxAxis"); | |
245 | fVtxAxis->SetTitle("v_{z} range unspecified"); | |
246 | } | |
b2e7f2d6 | 247 | |
e1f47419 | 248 | Info("FetchInformation", |
fe52e455 | 249 | "Initialized for\n" |
250 | " Trigger: %s (%d)\n" | |
251 | " sqrt(sNN): %s (%dGeV)\n" | |
252 | " System: %s (%d)\n" | |
253 | " Vz range: %s (%f,%f)", | |
254 | fTrigString->GetTitle(), fTrigString->GetUniqueID(), | |
255 | fSNNString->GetTitle(), fSNNString->GetUniqueID(), | |
256 | fSysString->GetTitle(), fSysString->GetUniqueID(), | |
257 | fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax()); | |
e1f47419 | 258 | } |
259 | ||
260 | //__________________________________________________________________ | |
261 | /** | |
262 | * Open input file, and find data | |
263 | * | |
264 | * @param filename File name | |
265 | * | |
266 | * @return true on success | |
267 | */ | |
268 | Bool_t Open(const char* filename, Double_t& max) | |
269 | { | |
270 | TFile* file = TFile::Open(filename, "READ"); | |
271 | if (!file) { | |
272 | Error("Open", "Cannot open %s", filename); | |
b2e7f2d6 | 273 | return false; |
274 | } | |
e1f47419 | 275 | |
276 | TList* results = static_cast<TList*>(file->Get("ForwardResults")); | |
277 | if (!results) { | |
278 | Error("Open", "Couldn't find list ForwardResults"); | |
279 | return false; | |
280 | } | |
281 | ||
282 | // Get information on the run | |
283 | FetchInformation(results); | |
284 | ||
285 | TList* clusters = static_cast<TList*>(file->Get("CentralResults")); | |
286 | if (!clusters) Warning("Open", "Couldn't find list CentralResults"); | |
b2e7f2d6 | 287 | |
e1f47419 | 288 | fResults = new THStack("results", "Results"); |
289 | FetchResults(fResults, results, "Forward", max); | |
290 | FetchResults(fResults, clusters, "Central", max); | |
291 | ||
292 | TList* sums = static_cast<TList*>(file->Get("ForwardSums")); | |
293 | if (sums) fTriggers = FetchResult(sums, "triggers"); | |
b2e7f2d6 | 294 | |
e1f47419 | 295 | file->Close(); |
b2e7f2d6 | 296 | |
297 | return true; | |
298 | } | |
299 | //__________________________________________________________________ | |
e1f47419 | 300 | void FetchResults(THStack* stack, const TList* list, const char* name, |
301 | Double_t& max) | |
b2e7f2d6 | 302 | { |
e1f47419 | 303 | TH1* dndeta = FetchResult(list, Form("dndeta%s", name)); |
304 | TH1* dndetaMC = FetchResult(list, Form("dndeta%sMC", name)); | |
305 | TH1* dndetaTruth = FetchResult(list, "dndetaTruth"); | |
306 | TH1* dndetaSym = 0; | |
307 | TH1* dndetaMCSym = 0; | |
308 | ||
309 | max = TMath::Max(max, AddHistogram(stack, dndetaTruth, "e5 p")); | |
310 | max = TMath::Max(max, AddHistogram(stack, dndetaMC, "", dndetaMCSym)); | |
311 | max = TMath::Max(max, AddHistogram(stack, dndeta, "", dndetaSym)); | |
312 | ||
313 | Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f", | |
314 | dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max); | |
315 | ||
316 | if (dndetaMCSym) fNumerators.Add(dndetaMCSym); | |
317 | if (dndetaMC) fNumerators.Add(dndetaMC); | |
318 | if (dndetaSym) fNumerators.Add(dndetaSym); | |
319 | if (dndeta) fNumerators.Add(dndeta); | |
320 | if (dndetaTruth) fDenominators.Add(dndetaTruth); | |
b2e7f2d6 | 321 | } |
322 | //__________________________________________________________________ | |
2f92ff0e | 323 | /** |
324 | * Make a histogram stack of results | |
325 | * | |
326 | * @param max On return, the maximum value in the stack | |
327 | * | |
328 | * @return Newly allocated stack | |
329 | */ | |
e1f47419 | 330 | TMultiGraph* StackOther(Double_t& max) |
b2e7f2d6 | 331 | { |
332 | gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/OtherData.C"); | |
333 | Int_t error = 0; | |
334 | Bool_t onlya = (fShowOthers ? false : true); | |
335 | Int_t trg = (fTrigString ? fTrigString->GetUniqueID() : 0); | |
336 | UShort_t snn = (fSNNString ? fSNNString->GetUniqueID() : 0); | |
337 | Long_t ret = gROOT->ProcessLine(Form("GetData(%d,%d,%d);", | |
338 | snn,trg,onlya)); | |
339 | if (error) { | |
340 | Error("StackOther", "Failed to execute GetData(%d,%d,%d)", | |
341 | snn, trg, onlya); | |
342 | return 0; | |
343 | } | |
344 | if (!ret) { | |
345 | Warning("StackOther", "No other data found for sNN=%d, trigger=%d", | |
346 | snn, trg); | |
347 | return 0; | |
348 | } | |
349 | TMultiGraph* other = reinterpret_cast<TMultiGraph*>(ret); | |
350 | ||
351 | TGraphAsymmErrors* o = 0; | |
352 | TIter next(other->GetListOfGraphs()); | |
e1f47419 | 353 | while ((o = static_cast<TGraphAsymmErrors*>(next()))) { |
b2e7f2d6 | 354 | max = TMath::Max(max, TMath::MaxElement(o->GetN(), o->GetY())); |
e1f47419 | 355 | fDenominators.Add(o); |
356 | } | |
b2e7f2d6 | 357 | |
358 | return other; | |
359 | } | |
360 | //__________________________________________________________________ | |
2f92ff0e | 361 | /** |
362 | * Make a histogram stack of ratios of results to other data | |
363 | * | |
364 | * @param max On return, the maximum value in the stack | |
365 | * | |
366 | * @return Newly allocated stack | |
367 | */ | |
b2e7f2d6 | 368 | THStack* StackRatios(TMultiGraph* others, Double_t& max) |
369 | { | |
370 | THStack* ratios = new THStack("ratios", "Ratios"); | |
371 | ||
372 | if (others) { | |
e1f47419 | 373 | TObject* ua5_1 = 0; |
374 | TObject* ua5_2 = 0; | |
375 | TObject* alice = 0; | |
376 | TObject* cms = 0; | |
377 | TObject* denom = 0; | |
378 | TIter nextD(&fDenominators); | |
379 | while ((denom = nextD())) { | |
380 | TIter nextN(&fNumerators); | |
381 | TObject* numer = 0; | |
382 | while ((numer = nextN())) { | |
383 | ratios->Add(Ratio(numer, denom, max)); | |
384 | } | |
385 | TString oName(denom->GetName()); | |
b2e7f2d6 | 386 | oName.ToLower(); |
e1f47419 | 387 | if (oName.Contains("ua5")) { |
388 | if (ua5_1) ua5_2 = denom; | |
389 | else ua5_1 = denom; | |
390 | } | |
391 | if (oName.Contains("alice")) alice = denom; | |
392 | if (oName.Contains("cms")) cms = denom; | |
b2e7f2d6 | 393 | } |
394 | if (ua5_1 && alice) ratios->Add(Ratio(alice, ua5_1, max)); | |
395 | if (ua5_2 && alice) ratios->Add(Ratio(alice, ua5_2, max)); | |
396 | if (cms && alice) ratios->Add(Ratio(alice, cms, max)); | |
397 | } | |
b2e7f2d6 | 398 | // Check if we have ratios |
399 | if (!ratios->GetHists() || | |
400 | (ratios->GetHists()->GetEntries() <= 0)) { | |
401 | delete ratios; | |
402 | ratios = 0; | |
403 | } | |
404 | return ratios; | |
405 | } | |
406 | //__________________________________________________________________ | |
2f92ff0e | 407 | /** |
408 | * Make a histogram stack of the left-right asymmetry | |
409 | * | |
410 | * @param max On return, the maximum value in the stack | |
411 | * | |
412 | * @return Newly allocated stack | |
413 | */ | |
b2e7f2d6 | 414 | THStack* StackLeftRight(Double_t& max) |
415 | { | |
416 | THStack* ret = new THStack("leftright", "Left-right asymmetry"); | |
e1f47419 | 417 | // ret->Add(Asymmetry(fForward, max)); |
418 | // ret->Add(Asymmetry(fForwardMC, max)); | |
b2e7f2d6 | 419 | |
420 | if (!ret->GetHists() || | |
421 | (ret->GetHists()->GetEntries() <= 0)) { | |
422 | delete ret; | |
423 | ret = 0; | |
424 | } | |
425 | return ret; | |
426 | } | |
427 | //__________________________________________________________________ | |
2f92ff0e | 428 | /** |
429 | * Plot the results | |
430 | * | |
431 | * @param results Results | |
432 | * @param others Other data | |
433 | * @param max Max value | |
434 | * @param ratios Stack of ratios (optional) | |
435 | * @param rmax Maximum diviation from 1 of ratios | |
436 | * @param leftright Stack of left-right asymmetry (optional) | |
437 | * @param amax Maximum diviation from 1 of asymmetries | |
438 | */ | |
b2e7f2d6 | 439 | void Plot(THStack* results, |
440 | TMultiGraph* others, | |
441 | Double_t max, | |
442 | THStack* ratios, | |
443 | Double_t rmax, | |
444 | THStack* leftright, | |
445 | Double_t amax) | |
446 | { | |
447 | gStyle->SetOptTitle(0); | |
448 | gStyle->SetTitleFont(132, "xyz"); | |
449 | gStyle->SetLabelFont(132, "xyz"); | |
450 | ||
451 | Int_t h = 800; | |
452 | Int_t w = 800; // h / TMath::Sqrt(2); | |
453 | if (!ratios) w *= 1.4; | |
454 | if (!leftright) w *= 1.4; | |
455 | TCanvas* c = new TCanvas("Results", "Results", w, h); | |
456 | c->SetFillColor(0); | |
457 | c->SetBorderSize(0); | |
458 | c->SetBorderMode(0); | |
459 | ||
460 | Double_t y1 = 0; | |
461 | Double_t y2 = 0; | |
462 | Double_t y3 = 0; | |
463 | if (ratios) y1 = 0.3; | |
464 | if (leftright) { | |
465 | if (y1 > 0.0001) { | |
466 | y2 = 0.2; | |
467 | y1 = 0.4; | |
468 | } | |
469 | else { | |
470 | y1 = 0.2; | |
471 | y2 = y1; | |
472 | } | |
473 | } | |
474 | PlotResults(results, others, max, y1); | |
475 | c->cd(); | |
476 | ||
477 | PlotRatios(ratios, rmax, y2, y1); | |
478 | c->cd( ); | |
479 | ||
480 | PlotLeftRight(leftright, amax, y3, y2); | |
481 | c->cd(); | |
bad9a3c1 | 482 | |
483 | ||
484 | Int_t vMin = fVtxAxis->GetXmin(); | |
485 | Int_t vMax = fVtxAxis->GetXmax(); | |
486 | TString trg(fTrigString->GetTitle()); | |
487 | Int_t nev = fTriggers->GetBinContent(fTriggers->GetNbinsX()); | |
488 | trg = trg.Strip(TString::kBoth); | |
489 | TString base(Form("dndeta_%s_%s_%s_%c%02d%c%02dcm_%09dev", | |
490 | fSysString->GetTitle(), | |
491 | fSNNString->GetTitle(), | |
492 | trg.Data(), | |
493 | vMin < 0 ? 'm' : 'p', TMath::Abs(vMin), | |
494 | vMax < 0 ? 'm' : 'p', TMath::Abs(vMax), | |
495 | nev)); | |
496 | c->SaveAs(Form("%s.png", base.Data())); | |
497 | c->SaveAs(Form("%s.root", base.Data())); | |
498 | c->SaveAs(Form("%s.C", base.Data())); | |
b2e7f2d6 | 499 | } |
500 | //__________________________________________________________________ | |
2f92ff0e | 501 | /** |
502 | * Plot the results | |
503 | * | |
504 | * @param results Results | |
505 | * @param others Other data | |
506 | * @param max Maximum | |
507 | * @param yd Bottom position of pad | |
508 | */ | |
b2e7f2d6 | 509 | void PlotResults(THStack* results, TMultiGraph* others, |
510 | Double_t max, Double_t yd) | |
511 | { | |
512 | // Make a sub-pad for the result itself | |
513 | TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0); | |
514 | p1->SetTopMargin(0.05); | |
515 | p1->SetBorderSize(0); | |
516 | p1->SetBorderMode(0); | |
517 | p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1); | |
518 | p1->SetRightMargin(0.05); | |
519 | p1->SetGridx(); | |
520 | p1->SetTicks(1,1); | |
521 | p1->SetNumber(1); | |
522 | p1->Draw(); | |
523 | p1->cd(); | |
524 | ||
e1f47419 | 525 | Info("PlotResults", "Plotting results with max=%f", max); |
526 | results->ls(); | |
b2e7f2d6 | 527 | results->SetMaximum(1.15*max); |
528 | results->SetMinimum(yd > 0.00001 ? -0.1 : 0); | |
529 | ||
530 | FixAxis(results, 1/(1-yd)/1.7, "#frac{1}{N} #frac{dN_{ch}}{d#eta}"); | |
531 | ||
532 | p1->Clear(); | |
533 | results->DrawClone("nostack e1"); | |
534 | ||
535 | fRangeParam->fSlave1Axis = results->GetXaxis(); | |
536 | fRangeParam->fSlave1Pad = p1; | |
537 | ||
538 | // Draw other data | |
539 | if (others) { | |
540 | TGraphAsymmErrors* o = 0; | |
541 | TIter nextG(others->GetListOfGraphs()); | |
542 | while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) | |
543 | o->DrawClone("same p"); | |
544 | } | |
545 | ||
546 | // Make a legend in the main result pad | |
547 | TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35); | |
548 | l->SetNColumns(2); | |
549 | l->SetFillColor(0); | |
550 | l->SetFillStyle(0); | |
551 | l->SetBorderSize(0); | |
552 | l->SetTextFont(132); | |
553 | ||
554 | // Put a title on top | |
555 | TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data()); | |
556 | tit->SetNDC(); | |
557 | tit->SetTextFont(132); | |
558 | tit->SetTextSize(0.05); | |
559 | tit->Draw(); | |
560 | ||
561 | // Put a nice label in the plot | |
562 | TString eS; | |
bad9a3c1 | 563 | UShort_t snn = fSNNString->GetUniqueID(); |
564 | const char* sys = fSysString->GetTitle(); | |
b2e7f2d6 | 565 | if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000); |
566 | else eS = Form("%3dGeV", snn); | |
567 | TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s}=%s, %s", | |
568 | sys, | |
569 | eS.Data(), | |
bad9a3c1 | 570 | fTrigString->GetTitle())); |
b2e7f2d6 | 571 | tt->SetNDC(); |
572 | tt->SetTextFont(132); | |
573 | tt->SetTextAlign(33); | |
574 | tt->Draw(); | |
575 | ||
576 | // Put number of accepted events on the plot | |
577 | Int_t nev = fTriggers->GetBinContent(fTriggers->GetNbinsX()); | |
578 | TLatex* et = new TLatex(.93, .83, Form("%d events", nev)); | |
579 | et->SetNDC(); | |
580 | et->SetTextFont(132); | |
581 | et->SetTextAlign(33); | |
582 | et->Draw(); | |
583 | ||
584 | // Put number of accepted events on the plot | |
585 | if (fVtxAxis) { | |
586 | TLatex* vt = new TLatex(.93, .88, fVtxAxis->GetTitle()); | |
587 | vt->SetNDC(); | |
588 | vt->SetTextFont(132); | |
589 | vt->SetTextAlign(33); | |
590 | vt->Draw(); | |
591 | } | |
592 | // results->Draw("nostack e1 same"); | |
593 | ||
594 | fRangeParam->fSlave1Axis = FindXAxis(p1, results->GetName()); | |
595 | fRangeParam->fSlave1Pad = p1; | |
596 | ||
597 | ||
598 | // Mark the plot as preliminary | |
599 | TLatex* pt = new TLatex(.12, .93, "Preliminary"); | |
600 | pt->SetNDC(); | |
601 | pt->SetTextFont(22); | |
602 | pt->SetTextSize(0.07); | |
603 | pt->SetTextColor(kRed+1); | |
604 | pt->SetTextAlign(13); | |
605 | pt->Draw(); | |
606 | ||
607 | if (!gSystem->AccessPathName("ALICE.png")) { | |
bad9a3c1 | 608 | TPad* logo = new TPad("logo", "logo", .12, .65, .25, .85, 0, 0, 0); |
b2e7f2d6 | 609 | logo->SetFillStyle(0); |
610 | logo->Draw(); | |
611 | logo->cd(); | |
612 | TImage* i = TImage::Create(); | |
613 | i->ReadImage("ALICE.png"); | |
614 | i->Draw(); | |
615 | } | |
616 | p1->cd(); | |
617 | } | |
618 | //__________________________________________________________________ | |
2f92ff0e | 619 | /** |
620 | * Plot the ratios | |
621 | * | |
622 | * @param ratios Ratios to plot (if any) | |
623 | * @param max Maximum diviation from 1 | |
624 | * @param y1 Lower y coordinate of pad | |
625 | * @param y2 Upper y coordinate of pad | |
626 | */ | |
b2e7f2d6 | 627 | void PlotRatios(THStack* ratios, Double_t max, Double_t y1, Double_t y2) |
628 | { | |
629 | if (!ratios) return; | |
630 | bool isBottom = (y1 < 0.0001); | |
631 | Double_t yd = y2 - y1; | |
632 | // Make a sub-pad for the result itself | |
633 | TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0); | |
634 | p2->SetTopMargin(0.001); | |
635 | p2->SetRightMargin(0.05); | |
636 | p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001); | |
637 | p2->SetGridx(); | |
638 | p2->SetTicks(1,1); | |
639 | p2->SetNumber(2); | |
640 | p2->Draw(); | |
641 | p2->cd(); | |
642 | ||
643 | // Fix up axis | |
644 | FixAxis(ratios, 1/yd/1.7, "Ratios", 7); | |
645 | ||
646 | ratios->SetMaximum(1+TMath::Max(.22,1.05*max)); | |
647 | ratios->SetMinimum(1-TMath::Max(.32,1.05*max)); | |
648 | p2->Clear(); | |
649 | ratios->DrawClone("nostack e1"); | |
650 | ||
651 | ||
652 | // Make a legend | |
653 | TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9, | |
654 | isBottom ? .6 : .4); | |
655 | l2->SetNColumns(2); | |
656 | l2->SetFillColor(0); | |
657 | l2->SetFillStyle(0); | |
658 | l2->SetBorderSize(0); | |
659 | l2->SetTextFont(132); | |
660 | ||
661 | // Make a nice band from 0.9 to 1.1 | |
662 | TGraphErrors* band = new TGraphErrors(2); | |
e1f47419 | 663 | band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1); |
664 | band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1); | |
b2e7f2d6 | 665 | band->SetPointError(0, 0, .1); |
666 | band->SetPointError(1, 0, .1); | |
667 | band->SetFillColor(kYellow+2); | |
668 | band->SetFillStyle(3002); | |
669 | band->SetLineStyle(2); | |
670 | band->SetLineWidth(1); | |
671 | band->Draw("3 same"); | |
672 | band->DrawClone("X L same"); | |
673 | ||
674 | // Replot the ratios on top | |
675 | ratios->DrawClone("nostack e1 same"); | |
676 | ||
677 | if (isBottom) { | |
678 | fRangeParam->fMasterAxis = FindXAxis(p2, ratios->GetName()); | |
679 | p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)", | |
680 | fRangeParam)); | |
681 | } | |
682 | else { | |
683 | fRangeParam->fSlave2Axis = FindXAxis(p2, ratios->GetName()); | |
684 | fRangeParam->fSlave2Pad = p2; | |
685 | } | |
686 | } | |
687 | //__________________________________________________________________ | |
2f92ff0e | 688 | /** |
689 | * Plot the asymmetries | |
690 | * | |
691 | * @param ratios Asymmetries to plot (if any) | |
692 | * @param max Maximum diviation from 1 | |
693 | * @param y1 Lower y coordinate of pad | |
694 | * @param y2 Upper y coordinate of pad | |
695 | */ | |
b2e7f2d6 | 696 | void PlotLeftRight(THStack* leftright, Double_t max, |
697 | Double_t y1, Double_t y2) | |
698 | { | |
699 | if (!leftright) return; | |
700 | bool isBottom = (y1 < 0.0001); | |
701 | Double_t yd = y2 - y1; | |
702 | // Make a sub-pad for the result itself | |
703 | TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0); | |
704 | p3->SetTopMargin(0.001); | |
705 | p3->SetRightMargin(0.05); | |
706 | p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001); | |
707 | p3->SetGridx(); | |
708 | p3->SetTicks(1,1); | |
709 | p3->SetNumber(2); | |
710 | p3->Draw(); | |
711 | p3->cd(); | |
712 | ||
713 | TH1* dummy = 0; | |
714 | if (leftright->GetHists()->GetEntries() == 1) { | |
715 | // Add dummy histogram | |
716 | dummy = new TH1F("dummy","", 10, -6, 6); | |
717 | dummy->SetLineColor(0); | |
718 | dummy->SetFillColor(0); | |
719 | dummy->SetMarkerColor(0); | |
720 | leftright->Add(dummy); | |
721 | } | |
722 | ||
723 | // Fix up axis | |
724 | FixAxis(leftright, 1/yd/1.7, "Right/Left", 4); | |
725 | ||
726 | leftright->SetMaximum(1+TMath::Max(.12,1.05*max)); | |
727 | leftright->SetMinimum(1-TMath::Max(.15,1.05*max)); | |
728 | p3->Clear(); | |
729 | leftright->DrawClone("nostack e1"); | |
730 | ||
731 | ||
732 | // Make a legend | |
733 | TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5); | |
734 | l2->SetNColumns(2); | |
735 | l2->SetFillColor(0); | |
736 | l2->SetFillStyle(0); | |
737 | l2->SetBorderSize(0); | |
738 | l2->SetTextFont(132); | |
739 | #ifndef __CINT__ | |
740 | if (dummy) { | |
741 | TList* prims = l2->GetListOfPrimitives(); | |
742 | TIter next(prims); | |
743 | TLegendEntry* o = 0; | |
744 | while ((o = static_cast<TLegendEntry*>(next()))) { | |
745 | TString lbl(o->GetLabel()); | |
746 | if (lbl != "dummy") continue; | |
747 | prims->Remove(o); | |
748 | break; | |
749 | } | |
750 | } | |
751 | #endif | |
752 | // Make a nice band from 0.9 to 1.1 | |
753 | TGraphErrors* band = new TGraphErrors(2); | |
e1f47419 | 754 | band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1); |
755 | band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1); | |
b2e7f2d6 | 756 | band->SetPointError(0, 0, .05); |
757 | band->SetPointError(1, 0, .05); | |
758 | band->SetFillColor(kYellow+2); | |
759 | band->SetFillStyle(3002); | |
760 | band->SetLineStyle(2); | |
761 | band->SetLineWidth(1); | |
762 | band->Draw("3 same"); | |
763 | band->DrawClone("X L same"); | |
764 | ||
765 | leftright->DrawClone("nostack e1 same"); | |
766 | if (isBottom) { | |
767 | fRangeParam->fMasterAxis = FindXAxis(p3, leftright->GetName()); | |
768 | p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)", | |
769 | fRangeParam)); | |
770 | } | |
771 | else { | |
772 | fRangeParam->fSlave2Axis = FindXAxis(p3, leftright->GetName()); | |
773 | fRangeParam->fSlave2Pad = p3; | |
774 | } | |
775 | } | |
2f92ff0e | 776 | /** @} */ |
777 | //================================================================== | |
778 | /** | |
779 | * @{ | |
780 | * @name Data utility functions | |
781 | */ | |
782 | /** | |
783 | * Get a result from the passed list | |
784 | * | |
785 | * @param list List to search | |
786 | * @param name Object name to search for | |
787 | * | |
788 | * @return | |
789 | */ | |
e1f47419 | 790 | TH1* FetchResult(const TList* list, const char* name) const |
2f92ff0e | 791 | { |
792 | if (!list) return 0; | |
793 | ||
e1f47419 | 794 | TList* all = static_cast<TList*>(list->FindObject("all")); |
795 | if (!all) { | |
796 | Warning("GetResult", "No 'all' list find in %s", list->GetName()); | |
797 | // list->ls(); | |
2f92ff0e | 798 | return 0; |
b2e7f2d6 | 799 | } |
e1f47419 | 800 | |
801 | TH1* ret = static_cast<TH1*>(all->FindObject(name)); | |
802 | if (!ret) { | |
803 | // all->ls(); | |
804 | Warning("GetResult", "Histogram %s not found", name); | |
2f92ff0e | 805 | } |
b2e7f2d6 | 806 | |
e1f47419 | 807 | return ret; |
2f92ff0e | 808 | } |
b2e7f2d6 | 809 | //__________________________________________________________________ |
2f92ff0e | 810 | /** |
811 | * Add a histogram to the stack after possibly rebinning it | |
812 | * | |
813 | * @param stack Stack to add to | |
814 | * @param hist histogram | |
815 | * @param option Draw options | |
816 | * | |
817 | * @return Maximum of histogram | |
818 | */ | |
b2e7f2d6 | 819 | Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option) const |
820 | { | |
821 | // Check if we have input | |
822 | if (!hist) return 0; | |
823 | ||
824 | // Rebin if needed | |
825 | Rebin(hist); | |
826 | ||
827 | stack->Add(hist, option); | |
828 | return hist->GetMaximum(); | |
829 | } | |
830 | //__________________________________________________________________ | |
2f92ff0e | 831 | /** |
832 | * Add a histogram to the stack after possibly rebinning it | |
833 | * | |
834 | * @param stack Stack to add to | |
835 | * @param hist histogram | |
836 | * @param option Draw options | |
837 | * @param sym On return, the data symmetriced (added to stack) | |
838 | * | |
839 | * @return Maximum of histogram | |
840 | */ | |
b2e7f2d6 | 841 | Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option, |
842 | TH1*& sym) const | |
843 | { | |
844 | // Check if we have input | |
845 | if (!hist) return 0; | |
846 | ||
847 | // Rebin if needed | |
848 | Rebin(hist); | |
849 | stack->Add(hist, option); | |
850 | ||
851 | // Now symmetrice the histogram | |
852 | sym = Symmetrice(hist); | |
853 | stack->Add(sym, option); | |
854 | ||
855 | return hist->GetMaximum(); | |
856 | } | |
857 | ||
858 | //__________________________________________________________________ | |
859 | /** | |
860 | * Rebin a histogram | |
861 | * | |
862 | * @param h Histogram to rebin | |
863 | * @param rebin Rebinning factor | |
864 | * | |
865 | * @return | |
866 | */ | |
867 | virtual void Rebin(TH1* h) const | |
868 | { | |
869 | if (fRebin <= 1) return; | |
870 | ||
871 | Int_t nBins = h->GetNbinsX(); | |
872 | if(nBins % fRebin != 0) { | |
873 | Warning("Rebin", "Rebin factor %d is not a devisor of current number " | |
874 | "of bins %d in the histogram %s", fRebin, nBins, h->GetName()); | |
875 | return; | |
876 | } | |
877 | ||
878 | // Make a copy | |
879 | TH1* tmp = static_cast<TH1*>(h->Clone("tmp")); | |
880 | tmp->Rebin(fRebin); | |
881 | tmp->SetDirectory(0); | |
e28f5fc5 | 882 | tmp->Reset(); |
b2e7f2d6 | 883 | // The new number of bins |
884 | Int_t nBinsNew = nBins / fRebin; | |
885 | for(Int_t i = 1;i<= nBinsNew; i++) { | |
886 | Double_t content = 0; | |
887 | Double_t sumw = 0; | |
888 | Double_t wsum = 0; | |
889 | Int_t nbins = 0; | |
890 | for(Int_t j = 1; j<=fRebin;j++) { | |
891 | Int_t bin = (i-1)*fRebin + j; | |
892 | Double_t c = h->GetBinContent(bin); | |
893 | ||
894 | if (c <= 0) continue; | |
895 | ||
896 | if (fCutEdges) { | |
897 | if (h->GetBinContent(bin+1)<=0 || | |
898 | h->GetBinContent(bin-1)) { | |
899 | Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)", | |
900 | bin, c, h->GetName(), | |
901 | bin+1, h->GetBinContent(bin+1), | |
902 | bin-1, h->GetBinContent(bin-1)); | |
903 | continue; | |
904 | } | |
905 | } | |
906 | Double_t e = h->GetBinError(bin); | |
907 | Double_t w = 1 / (e*e); // 1/c/c | |
908 | content += c; | |
909 | sumw += w; | |
910 | wsum += w * c; | |
911 | nbins++; | |
912 | } | |
913 | ||
e28f5fc5 | 914 | if(content > 0 && nbins > 1 ) { |
b2e7f2d6 | 915 | tmp->SetBinContent(i, wsum / sumw); |
916 | tmp->SetBinError(i,1./TMath::Sqrt(sumw)); | |
917 | } | |
918 | } | |
919 | ||
920 | // Finally, rebin the histogram, and set new content | |
921 | h->Rebin(fRebin); | |
922 | h->Reset(); | |
923 | for(Int_t i = 1; i<= nBinsNew; i++) { | |
924 | h->SetBinContent(i,tmp->GetBinContent(i)); | |
925 | h->SetBinError(i, tmp->GetBinError(i)); | |
926 | } | |
927 | ||
928 | delete tmp; | |
929 | } | |
930 | //__________________________________________________________________ | |
931 | /** | |
932 | * Make an extension of @a h to make it symmetric about 0 | |
933 | * | |
934 | * @param h Histogram to symmertrice | |
935 | * | |
936 | * @return Symmetric extension of @a h | |
937 | */ | |
938 | TH1* Symmetrice(const TH1* h) const | |
939 | { | |
940 | Int_t nBins = h->GetNbinsX(); | |
941 | TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName()))); | |
942 | s->SetTitle(Form("%s (mirrored)", h->GetTitle())); | |
943 | s->Reset(); | |
944 | s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin()); | |
945 | s->SetMarkerColor(h->GetMarkerColor()); | |
946 | s->SetMarkerSize(h->GetMarkerSize()); | |
947 | s->SetMarkerStyle(h->GetMarkerStyle()+4); | |
948 | s->SetFillColor(h->GetFillColor()); | |
949 | s->SetFillStyle(h->GetFillStyle()); | |
950 | s->SetDirectory(0); | |
951 | ||
952 | // Find the first and last bin with data | |
953 | Int_t first = nBins+1; | |
954 | Int_t last = 0; | |
955 | for (Int_t i = 1; i <= nBins; i++) { | |
956 | if (h->GetBinContent(i) <= 0) continue; | |
957 | first = TMath::Min(first, i); | |
958 | last = TMath::Max(last, i); | |
959 | } | |
960 | ||
961 | Double_t xfirst = h->GetBinCenter(first-1); | |
962 | Int_t f1 = h->GetXaxis()->FindBin(-xfirst); | |
963 | Int_t l2 = s->GetXaxis()->FindBin(xfirst); | |
964 | for (Int_t i = f1, j=l2; i <= last; i++,j--) { | |
965 | s->SetBinContent(j, h->GetBinContent(i)); | |
966 | s->SetBinError(j, h->GetBinError(i)); | |
967 | } | |
968 | // Fill in overlap bin | |
969 | s->SetBinContent(l2+1, h->GetBinContent(first)); | |
970 | s->SetBinError(l2+1, h->GetBinError(first)); | |
971 | return s; | |
972 | } | |
973 | //__________________________________________________________________ | |
2f92ff0e | 974 | /** |
975 | * Calculate the left-right asymmetry of input histogram | |
976 | * | |
977 | * @param h Input histogram | |
978 | * @param max On return, the maximum distance from 1 of the histogram | |
979 | * | |
980 | * @return Asymmetry | |
981 | */ | |
982 | TH1* Asymmetry(TH1* h, Double_t& max) | |
983 | { | |
984 | if (!h) return 0; | |
985 | ||
986 | TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName()))); | |
987 | // Int_t oBins = h->GetNbinsX(); | |
988 | // Double_t high = h->GetXaxis()->GetXmax(); | |
989 | // Double_t low = h->GetXaxis()->GetXmin(); | |
990 | // Double_t dBin = (high - low) / oBins; | |
991 | // Int_t tBins = Int_t(2*high/dBin+.5); | |
992 | // ret->SetBins(tBins, -high, high); | |
993 | ret->Reset(); | |
994 | ret->SetTitle(Form("%s (+/-)", h->GetTitle())); | |
995 | ret->SetYTitle("Right/Left"); | |
996 | Int_t nBins = h->GetNbinsX(); | |
997 | for (Int_t i = 1; i <= nBins; i++) { | |
998 | Double_t x = h->GetBinCenter(i); | |
999 | if (x > 0) break; | |
1000 | ||
1001 | Double_t c1 = h->GetBinContent(i); | |
1002 | Double_t e1 = h->GetBinError(i); | |
1003 | if (c1 <= 0) continue; | |
1004 | ||
1005 | Int_t j = h->FindBin(-x); | |
1006 | if (j <= 0 || j > nBins) continue; | |
1007 | ||
1008 | Double_t c2 = h->GetBinContent(j); | |
1009 | Double_t e2 = h->GetBinError(j); | |
1010 | ||
1011 | Double_t c12 = c1*c1; | |
1012 | Double_t e = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12)); | |
1013 | ||
1014 | Int_t k = ret->FindBin(x); | |
1015 | ret->SetBinContent(k, c2/c1); | |
1016 | ret->SetBinError(k, e); | |
1017 | } | |
1018 | max = TMath::Max(max, RatioMax(ret)); | |
1019 | ||
1020 | return ret; | |
1021 | } | |
1022 | //__________________________________________________________________ | |
1023 | /** | |
1024 | * Transform a graph into a histogram | |
1025 | * | |
1026 | * @param g | |
1027 | * | |
1028 | * @return | |
1029 | */ | |
1030 | TH1* Graph2Hist(const TGraphAsymmErrors* g) const | |
1031 | { | |
1032 | Int_t nBins = g->GetN(); | |
1033 | TArrayF bins(nBins+1); | |
1034 | Double_t dx = 0; | |
1035 | for (Int_t i = 0; i < nBins; i++) { | |
1036 | Double_t x = g->GetX()[i]; | |
1037 | Double_t exl = g->GetEXlow()[i]; | |
1038 | Double_t exh = g->GetEXhigh()[i]; | |
1039 | bins.fArray[i] = x-exl; | |
1040 | bins.fArray[i+1] = x+exh; | |
1041 | Double_t dxi = exh+exl; | |
1042 | if (i == 0) dx = dxi; | |
1043 | else if (dxi != dx) dx = 0; | |
1044 | } | |
1045 | TString name(g->GetName()); | |
1046 | TString title(g->GetTitle()); | |
1047 | TH1D* h = 0; | |
1048 | if (dx != 0) { | |
1049 | h = new TH1D(name.Data(), title.Data(), nBins, bins[0], bins[nBins]); | |
1050 | } | |
1051 | else { | |
1052 | h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray); | |
1053 | } | |
1054 | h->SetMarkerStyle(g->GetMarkerStyle()); | |
1055 | h->SetMarkerColor(g->GetMarkerColor()); | |
1056 | h->SetMarkerSize(g->GetMarkerSize()); | |
1057 | ||
1058 | return h; | |
1059 | } | |
1060 | /* @} */ | |
1061 | //================================================================== | |
1062 | /** | |
1063 | * @{ | |
1064 | * @name Ratio utility functions | |
1065 | */ | |
1066 | /** | |
1067 | * Get the maximum diviation from 1 in the passed ratio | |
1068 | * | |
1069 | * @param h Ratio histogram | |
1070 | * | |
1071 | * @return Max diviation from 1 | |
1072 | */ | |
b2e7f2d6 | 1073 | Double_t RatioMax(TH1* h) const |
1074 | { | |
1075 | Int_t nBins = h->GetNbinsX(); | |
1076 | Double_t ret = 0; | |
1077 | for (Int_t i = 1; i <= nBins; i++) { | |
1078 | Double_t c = h->GetBinContent(i); | |
1079 | if (c == 0) continue; | |
1080 | Double_t e = h->GetBinError(i); | |
1081 | Double_t d = TMath::Abs(1-c-e); | |
1082 | ret = TMath::Max(d, ret); | |
1083 | } | |
1084 | return ret; | |
1085 | } | |
1086 | //__________________________________________________________________ | |
e1f47419 | 1087 | /** |
1088 | * Make the ratio of h1 to h2 | |
1089 | * | |
1090 | * @param h1 First histogram (numerator) | |
1091 | * @param h2 Second histogram (denominator) | |
1092 | * | |
1093 | * @return h1 / h2 | |
1094 | */ | |
1095 | TH1* Ratio(const TObject* o1, const TObject* o2, Double_t& max) const | |
1096 | { | |
1097 | const TH1* h1 = dynamic_cast<const TH1*>(o1); | |
1098 | if (h1) { | |
1099 | const TH1* h2 = dynamic_cast<const TH1*>(o2); | |
1100 | if (h2) return Ratio(h1,h2,max); | |
1101 | ||
1102 | const TGraph* g2 = dynamic_cast<const TGraph*>(o2); | |
1103 | if (g2) return Ratio(h1,g2,max); | |
1104 | } | |
1105 | ||
1106 | const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(o1); | |
1107 | if (g1) { | |
1108 | const TGraphAsymmErrors* g2 = dynamic_cast<const TGraphAsymmErrors*>(o2); | |
1109 | if (g2) return Ratio(g1, g2, max); | |
1110 | } | |
1111 | ||
1112 | Warning("Ratio", "Don't know how to divide a %s (%s) with a %s (%s)", | |
1113 | o1->ClassName(), o1->GetName(), o2->ClassName(), o2->GetName()); | |
1114 | return 0; | |
1115 | } | |
1116 | //__________________________________________________________________ | |
b2e7f2d6 | 1117 | /** |
1118 | * Compute the ratio of @a h to @a g. @a g is evaluated at the bin | |
1119 | * centers of @a h | |
1120 | * | |
1121 | * @param h Numerator | |
1122 | * @param g Divisor | |
1123 | * | |
1124 | * @return h/g | |
1125 | */ | |
1126 | TH1* Ratio(const TH1* h, const TGraph* g, Double_t& max) const | |
1127 | { | |
1128 | if (!h || !g) return 0; | |
1129 | ||
1130 | TH1* ret = static_cast<TH1*>(h->Clone("tmp")); | |
1131 | ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName())); | |
1132 | ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle())); | |
1133 | ret->Reset(); | |
1134 | ret->SetMarkerStyle(g->GetMarkerStyle()); | |
1135 | ret->SetMarkerColor(h->GetMarkerColor()); | |
1136 | ret->SetMarkerSize(0.9*g->GetMarkerSize()); | |
1137 | Double_t xlow = g->GetX()[0]; | |
1138 | Double_t xhigh = g->GetX()[g->GetN()-1]; | |
1139 | if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; } | |
1140 | ||
1141 | for (Int_t i = 1; i <= h->GetNbinsX(); i++) { | |
1142 | Double_t c = h->GetBinContent(i); | |
1143 | if (c <= 0) continue; | |
1144 | ||
1145 | Double_t x = h->GetBinCenter(i); | |
1146 | if (x < xlow || x > xhigh) continue; | |
1147 | ||
1148 | Double_t f = g->Eval(x); | |
1149 | if (f <= 0) continue; | |
1150 | ||
1151 | ret->SetBinContent(i, c / f); | |
1152 | ret->SetBinError(i, h->GetBinError(i) / f); | |
1153 | } | |
1154 | if (ret->GetEntries() <= 0) { | |
1155 | delete ret; | |
1156 | ret = 0; | |
1157 | } | |
1158 | else | |
1159 | max = TMath::Max(RatioMax(ret), max); | |
1160 | return ret; | |
1161 | } | |
1162 | //__________________________________________________________________ | |
1163 | /** | |
1164 | * Make the ratio of h1 to h2 | |
1165 | * | |
1166 | * @param h1 First histogram (numerator) | |
1167 | * @param h2 Second histogram (denominator) | |
1168 | * | |
1169 | * @return h1 / h2 | |
1170 | */ | |
1171 | TH1* Ratio(const TH1* h1, const TH1* h2, Double_t& max) const | |
1172 | { | |
1173 | if (!h1 || !h2) return 0; | |
1174 | TH1* t1 = static_cast<TH1*>(h1->Clone(Form("%s_%s", | |
1175 | h1->GetName(), | |
1176 | h2->GetName()))); | |
1177 | t1->SetTitle(Form("%s / %s", h1->GetTitle(), h2->GetTitle())); | |
1178 | t1->Divide(h2); | |
1179 | t1->SetMarkerColor(h1->GetMarkerColor()); | |
1180 | t1->SetMarkerStyle(h2->GetMarkerStyle()); | |
1181 | max = TMath::Max(RatioMax(t1), max); | |
1182 | return t1; | |
1183 | } | |
1184 | //__________________________________________________________________ | |
1185 | /** | |
1186 | * Calculate the ratio of two graphs - g1 / g2 | |
1187 | * | |
1188 | * @param g1 Numerator | |
1189 | * @param g2 Denominator | |
1190 | * | |
1191 | * @return g1 / g2 in a histogram | |
1192 | */ | |
1193 | TH1* Ratio(const TGraphAsymmErrors* g1, | |
1194 | const TGraphAsymmErrors* g2, Double_t& max) const | |
1195 | { | |
1196 | Int_t nBins = g1->GetN(); | |
1197 | TArrayF bins(nBins+1); | |
1198 | Double_t dx = 0; | |
1199 | for (Int_t i = 0; i < nBins; i++) { | |
1200 | Double_t x = g1->GetX()[i]; | |
1201 | Double_t exl = g1->GetEXlow()[i]; | |
1202 | Double_t exh = g1->GetEXhigh()[i]; | |
1203 | bins.fArray[i] = x-exl; | |
1204 | bins.fArray[i+1] = x+exh; | |
1205 | Double_t dxi = exh+exl; | |
1206 | if (i == 0) dx = dxi; | |
1207 | else if (dxi != dx) dx = 0; | |
1208 | } | |
1209 | TString name(Form("%s_%s", g1->GetName(), g2->GetName())); | |
1210 | TString title(Form("%s / %s", g1->GetTitle(), g2->GetTitle())); | |
1211 | TH1* h = 0; | |
1212 | if (dx != 0) { | |
1213 | h = new TH1F(name.Data(), title.Data(), nBins, bins[0], bins[nBins]); | |
1214 | } | |
1215 | else { | |
1216 | h = new TH1F(name.Data(), title.Data(), nBins, bins.fArray); | |
1217 | } | |
1218 | h->SetMarkerStyle(g2->GetMarkerStyle()); | |
1219 | h->SetMarkerColor(g1->GetMarkerColor()); | |
1220 | h->SetMarkerSize(0.9*g2->GetMarkerSize()); | |
1221 | ||
1222 | Double_t low = g2->GetX()[0]; | |
1223 | Double_t high = g2->GetX()[g2->GetN()-1]; | |
1224 | if (low > high) { Double_t t = low; low = high; high = t; } | |
1225 | for (Int_t i = 0; i < nBins; i++) { | |
1226 | Double_t x = g1->GetX()[i]; | |
1227 | if (x < low-dx || x > high+dx) continue; | |
1228 | Double_t c1 = g1->GetY()[i]; | |
1229 | Double_t e1 = g1->GetErrorY(i); | |
1230 | Double_t c2 = g2->Eval(x); | |
1231 | ||
1232 | h->SetBinContent(i+1, c1 / c2); | |
1233 | h->SetBinError(i+1, e1 / c2); | |
1234 | } | |
1235 | max = TMath::Max(RatioMax(h), max); | |
1236 | return h; | |
2f92ff0e | 1237 | } |
1238 | /* @} */ | |
1239 | //================================================================== | |
1240 | /** | |
1241 | * @{ | |
1242 | * @name Graphics utility functions | |
1243 | */ | |
1244 | /** | |
1245 | * Find an X axis in a pad | |
1246 | * | |
1247 | * @param p Pad | |
1248 | * @param name Histogram to find axis for | |
1249 | * | |
1250 | * @return Found axis or null | |
1251 | */ | |
1252 | TAxis* FindXAxis(TVirtualPad* p, const char* name) | |
b2e7f2d6 | 1253 | { |
2f92ff0e | 1254 | TObject* o = p->GetListOfPrimitives()->FindObject(name); |
1255 | if (!o) { | |
1256 | Warning("FindXAxis", "%s not found in pad", name); | |
1257 | return 0; | |
b2e7f2d6 | 1258 | } |
2f92ff0e | 1259 | THStack* stack = dynamic_cast<THStack*>(o); |
1260 | if (!stack) { | |
1261 | Warning("FindXAxis", "%s is not a THStack", name); | |
1262 | return 0; | |
b2e7f2d6 | 1263 | } |
2f92ff0e | 1264 | if (!stack->GetHistogram()) { |
1265 | Warning("FindXAxis", "%s has no histogram", name); | |
1266 | return 0; | |
b2e7f2d6 | 1267 | } |
2f92ff0e | 1268 | TAxis* ret = stack->GetHistogram()->GetXaxis(); |
1269 | return ret; | |
b2e7f2d6 | 1270 | } |
2f92ff0e | 1271 | |
b2e7f2d6 | 1272 | //__________________________________________________________________ |
2f92ff0e | 1273 | /** |
1274 | * Fix the apperance of the axis in a stack | |
1275 | * | |
1276 | * @param stack stack of histogram | |
1277 | * @param s Scaling factor | |
1278 | * @param ytitle Y axis title | |
1279 | * @param force Whether to draw the stack first or not | |
1280 | * @param ynDiv Divisions on Y axis | |
1281 | */ | |
1282 | void FixAxis(THStack* stack, Double_t s, const char* ytitle, | |
1283 | Int_t ynDiv=210, Bool_t force=true) | |
b2e7f2d6 | 1284 | { |
2f92ff0e | 1285 | if (force) stack->Draw("nostack e1"); |
b2e7f2d6 | 1286 | |
2f92ff0e | 1287 | TH1* h = stack->GetHistogram(); |
1288 | if (!h) return; | |
b2e7f2d6 | 1289 | |
2f92ff0e | 1290 | h->SetXTitle("#eta"); |
1291 | h->SetYTitle(ytitle); | |
1292 | TAxis* xa = h->GetXaxis(); | |
1293 | TAxis* ya = h->GetYaxis(); | |
1294 | if (xa) { | |
1295 | xa->SetTitle("#eta"); | |
1296 | // xa->SetTicks("+-"); | |
1297 | xa->SetTitleSize(s*xa->GetTitleSize()); | |
1298 | xa->SetLabelSize(s*xa->GetLabelSize()); | |
1299 | xa->SetTickLength(s*xa->GetTickLength()); | |
1300 | } | |
1301 | if (ya) { | |
1302 | ya->SetTitle(ytitle); | |
1303 | ya->SetDecimals(); | |
1304 | // ya->SetTicks("+-"); | |
1305 | ya->SetNdivisions(ynDiv); | |
1306 | ya->SetTitleSize(s*ya->GetTitleSize()); | |
1307 | ya->SetTitleOffset(ya->GetTitleOffset()/s); | |
1308 | ya->SetLabelSize(s*ya->GetLabelSize()); | |
b2e7f2d6 | 1309 | } |
b2e7f2d6 | 1310 | } |
2f92ff0e | 1311 | /* @} */ |
1312 | ||
b2e7f2d6 | 1313 | |
1314 | ||
1315 | //__________________________________________________________________ | |
2f92ff0e | 1316 | Bool_t fShowOthers; // Show other data |
1317 | Bool_t fShowAlice; // Show ALICE published data | |
1318 | Bool_t fShowRatios; // Show ratios | |
1319 | Bool_t fShowLeftRight;// Show asymmetry | |
1320 | UShort_t fRebin; // Rebinning factor | |
1321 | Bool_t fCutEdges; // Whether to cut edges | |
1322 | TString fTitle; // Title on plot | |
2f92ff0e | 1323 | TNamed* fTrigString; // Trigger string (read, or set) |
1324 | TNamed* fSNNString; // Energy string (read, or set) | |
1325 | TNamed* fSysString; // Collision system string (read or set) | |
1326 | TAxis* fVtxAxis; // Vertex cuts (read or set) | |
e1f47419 | 1327 | TList fNumerators; // Numerators in ratios |
1328 | TList fDenominators; // Denominators in ratios | |
1329 | THStack* fResults; // Stack of results | |
1330 | THStack* fRatios; // Stack of ratios | |
1331 | // TH1* fForward; // Results | |
1332 | // TH1* fForwardMC; // MC results | |
1333 | // TH1* fForwardHHD; // Old results | |
1334 | // TH1* fTruth; // MC truth | |
1335 | // TH1* fCentral; // Central region data | |
1336 | // TH1* fForwardSym; // Symmetric extension | |
1337 | // TH1* fForwardMCSym; // Symmetric extension | |
1338 | // TH1* fForwardHHDSym;// Symmetric extension | |
2f92ff0e | 1339 | TH1* fTriggers; // Number of triggers |
1340 | RangeParam* fRangeParam; // Parameter object for range zoom | |
b2e7f2d6 | 1341 | |
1342 | }; | |
1343 | ||
1344 | //=== Stuff for auto zooming ========================================= | |
1345 | void UpdateRange(dNdetaDrawer::RangeParam* p) | |
1346 | { | |
1347 | if (!p) { | |
1348 | Warning("UpdateRange", "No parameters %p", p); | |
1349 | return; | |
1350 | } | |
1351 | if (!p->fMasterAxis) { | |
1352 | Warning("UpdateRange", "No master axis %p", p->fMasterAxis); | |
1353 | return; | |
1354 | } | |
1355 | Int_t first = p->fMasterAxis->GetFirst(); | |
1356 | Int_t last = p->fMasterAxis->GetLast(); | |
1357 | Double_t x1 = p->fMasterAxis->GetBinCenter(first); | |
1358 | Double_t x2 = p->fMasterAxis->GetBinCenter(last); | |
1359 | //Info("UpdateRange", "Range set to [%3d,%3d]->[%f,%f]", first, last, x1,x2); | |
1360 | ||
1361 | if (p->fSlave1Axis) { | |
1362 | Int_t i1 = p->fSlave1Axis->FindBin(x1); | |
1363 | Int_t i2 = p->fSlave1Axis->FindBin(x2); | |
1364 | p->fSlave1Axis->SetRange(i1, i2); | |
1365 | p->fSlave1Pad->Modified(); | |
1366 | p->fSlave1Pad->Update(); | |
1367 | } | |
1368 | if (p->fSlave2Axis) { | |
1369 | Int_t i1 = p->fSlave2Axis->FindBin(x1); | |
1370 | Int_t i2 = p->fSlave2Axis->FindBin(x2); | |
1371 | p->fSlave2Axis->SetRange(i1, i2); | |
1372 | p->fSlave2Pad->Modified(); | |
1373 | p->fSlave2Pad->Update(); | |
1374 | } | |
1375 | TCanvas* c = gPad->GetCanvas(); | |
1376 | c->cd(); | |
1377 | } | |
1378 | ||
1379 | //____________________________________________________________________ | |
1380 | void RangeExec(dNdetaDrawer::RangeParam* p) | |
1381 | { | |
1382 | // Event types: | |
1383 | // 51: Mouse move | |
1384 | // 53: | |
1385 | // 1: Button down | |
1386 | // 21: Mouse drag | |
1387 | // 11: Mouse release | |
1388 | // dNdetaDrawer::RangeParam* p = | |
1389 | // reinterpret_cast<dNdetaDrawer::RangeParam*>(addr); | |
1390 | Int_t event = gPad->GetEvent(); | |
1391 | TObject *select = gPad->GetSelected(); | |
1392 | if (event == 53) { | |
1393 | UpdateRange(p); | |
1394 | return; | |
1395 | } | |
1396 | if (event != 11 || !select || select != p->fMasterAxis) return; | |
1397 | UpdateRange(p); | |
1398 | } | |
1399 | ||
1400 | //=== Steering function ============================================== | |
1401 | void | |
1402 | DrawdNdeta(const char* filename="forward_dndeta.root", | |
1403 | Int_t flags=0xf, | |
bad9a3c1 | 1404 | const char* title="", |
b2e7f2d6 | 1405 | UShort_t rebin=5, |
fe52e455 | 1406 | Bool_t cutEdges=false, |
1407 | UShort_t sNN=0, | |
1408 | UShort_t sys=0, | |
1409 | UShort_t trg=1, | |
1410 | Float_t vzMin=999, | |
1411 | Float_t vzMax=-999) | |
b2e7f2d6 | 1412 | { |
1413 | dNdetaDrawer d; | |
1414 | d.SetRebin(rebin); | |
1415 | d.SetCutEdges(cutEdges); | |
1416 | d.SetTitle(title); | |
b2e7f2d6 | 1417 | d.SetShowOthers(flags & 0x1); |
1418 | d.SetShowAlice(flags & 0x2); | |
1419 | d.SetShowRatios(flags & 0x4); | |
1420 | d.SetShowLeftRight(flags & 0x8); | |
2f92ff0e | 1421 | // Do the below if your input data does not contain these settings |
fe52e455 | 1422 | if (sNN > 0) d.SetSNN(sNN); // Collision energy per nucleon pair (GeV) |
1423 | if (sys > 0) d.SetSys(sys); // Collision system (1:pp, 2:PbPB) | |
1424 | if (trg > 0) d.SetTrigger(trg); // Collision trigger (1:INEL, 2:INEL>0, 4:NSD) | |
1425 | if (vzMin < 999 && vzMax > -999) | |
1426 | d.SetVertexRange(vzMin,vzMax); // Collision vertex range (cm) | |
b2e7f2d6 | 1427 | d.Run(filename); |
1428 | } | |
2f92ff0e | 1429 | //____________________________________________________________________ |
1430 | // | |
1431 | // EOF | |
1432 | // | |
1433 |