Initialisation corrected.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / DrawdNdeta.C
CommitLineData
b2e7f2d6 1#include <TH1.h>
2#include <THStack.h>
3#include <TGraphErrors.h>
4#include <TGraphAsymmErrors.h>
5#include <TMultiGraph.h>
6#include <TFile.h>
7#include <TList.h>
8#include <TString.h>
9#include <TError.h>
10#include <TSystem.h>
11#include <TROOT.h>
12#include <TMath.h>
13#include <TCanvas.h>
14#include <TPad.h>
15#include <TStyle.h>
16#include <TLegend.h>
17#include <TLegendEntry.h>
18#include <TLatex.h>
19#include <TImage.h>
20
21struct dNdetaDrawer
22{
23 struct RangeParam
24 {
25 TAxis* fMasterAxis; // Master axis
26 TAxis* fSlave1Axis; // First slave axis
27 TVirtualPad* fSlave1Pad; // First slave pad
28 TAxis* fSlave2Axis; // Second slave axis
29 TVirtualPad* fSlave2Pad; // Second slave pad
30 };
31 //__________________________________________________________________
32 dNdetaDrawer()
33 : fShowOthers(false), // Bool_t
34 fShowAlice(false), // Bool_t
35 fShowRatios(false), // Bool_t
36 fShowLeftRight(false), // Bool_t
37 fRebin(5), // UShort_t
38 fCutEdges(false), // Bool_t
39 fTitle(""), // TString
40 fHHDFile(""), // TString
41 fTrigString(0), // TNamed*
42 fSNNString(0), // TNamed*
43 fSysString(0), // TNamed*
44 fVtxAxis(0), // TAxis*
45 fForward(0), // TH1*
46 fForwardMC(0), // TH1*
47 fForwardHHD(0), // TH1*
48 fTruth(0), // TH1*
49 fCentral(0), // TH1*
50 fForwardSym(0), // TH1*
51 fForwardMCSym(0), // TH1*
52 fForwardHHDSym(0), // TH1*
53 fTriggers(0), // TH1*
54 fRangeParam(0)
55
56 {
57 fRangeParam = new RangeParam;
58 fRangeParam->fMasterAxis = 0;
59 fRangeParam->fSlave1Axis = 0;
60 fRangeParam->fSlave1Pad = 0;
61 fRangeParam->fSlave2Axis = 0;
62 fRangeParam->fSlave2Pad = 0;
63 }
64 void SetShowOthers(Bool_t x) { fShowOthers = x; }
65 void SetShowAlice(Bool_t x) { fShowAlice = x; }
66 void SetShowRatios(Bool_t x) { fShowRatios = x; }
67 void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; }
68 void SetRebin(UShort_t x) { fRebin = x; }
69 void SetCutEdges(Bool_t x) { fCutEdges = x; }
70 void SetTitle(TString x) { fTitle = x; }
71 void SetHHDFile(const char* fn) { fHHDFile = fn; }
72
73 //__________________________________________________________________
74 void Run(const char* filename)
75 {
76 if (!Open(filename)) return;
77
78 Double_t max = 0;
79
80 // Create our stack of results
81 THStack* results = StackResults(max);
82
83 // Create our stack of other results
84 TMultiGraph* other = 0;
85 if (fShowOthers || fShowAlice) other = StackOther(max);
86
87 Double_t smax = 0;
88 THStack* ratios = 0;
89 if (fShowRatios) ratios = StackRatios(other, smax);
90
91 Double_t amax = 0;
92 THStack* leftright = 0;
93 if (fShowLeftRight) leftright = StackLeftRight(amax);
94
95 Plot(results, other, max, ratios, smax, leftright, amax);
96 }
97
98 //__________________________________________________________________
99 Bool_t Open(const char* filename)
100 {
101 TFile* file = TFile::Open(filename, "READ");
102 if (!file) {
103 Error("Open", "Cannot open %s", filename);
104 return false;
105 }
106
107 TList* results = static_cast<TList*>(file->Get("ForwardResults"));
108 if (!results) {
109 Error("Open", "Couldn't find list ForwardResults");
110 return false;
111 }
112
113 fForward = GetResult(results, "dndetaForward");
114 fForwardMC = GetResult(results, "dndetaForwardMC");
115 fTruth = GetResult(results, "dndetaTruth");
116 fCentral = GetResult(results, "dndetaCentral");
117
118 fTrigString = static_cast<TNamed*>(results->FindObject("trigString"));
119 fSNNString = static_cast<TNamed*>(results->FindObject("sNN"));
120 fSysString = static_cast<TNamed*>(results->FindObject("sys"));
121 fVtxAxis = static_cast<TAxis*>(results->FindObject("vtxAxis"));
bad9a3c1 122
123 if (!fTrigString) fTrigString = new TNamed("trigString", "unknown");
124 if (!fSNNString) fSNNString = new TNamed("sNN", "unknown");
125 if (!fSysString) fSysString = new TNamed("sys", "unknown");
126 if (!fVtxAxis) {
127 fVtxAxis = new TAxis(1,0,0);
128 fVtxAxis->SetName("vtxAxis");
129 fVtxAxis->SetTitle("v_{z} range unspecified");
130 }
b2e7f2d6 131
132 TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
133 if (sums)
134 fTriggers = GetResult(sums, "triggers");
135
136 if (!fForward) {
137 Error("Open", "Couldn't find the result of the forward analysis");
138 return false;
139 }
140 file->Close();
141
142
143 fForwardHHD = GetHHD();
144
145 return true;
146 }
147 //__________________________________________________________________
148 TH1* GetResult(TList* list, const char* name) const
149 {
150 if (!list) return 0;
151
152 TH1* ret = static_cast<TH1*>(list->FindObject(name));
153 if (!ret)
154 Warning("GetResult", "Histogram %s not found", name);
155
156 return ret;
157 }
158 //__________________________________________________________________
159 /**
160 * Get the result from previous analysis code
161 *
162 * @param fn File to open
163 * @param nsd Whether this is NSD
164 *
165 * @return null or result of previous analysis code
166 */
167 TH1* GetHHD()
168 {
169 if (fHHDFile.IsNull()) return 0;
170 const char* fn = fHHDFile.Data();
171 Bool_t nsd = (fTrigString ? fTrigString->GetUniqueID() & 0x4 : false);
172 TDirectory* savdir = gDirectory;
173 if (gSystem->AccessPathName(fn)) {
174 Warning("GetHHD", "Output of HHD analysis (%s) not available", fn);
175 return 0;
176 }
177 TFile* file = TFile::Open(fn, "READ");
178 if (!file) {
179 Warning("GetHHD", "couldn't open HHD file %s", fn);
180 return 0;
181 }
182 TString hist(Form("dNdeta_dNdeta%s", (nsd ? "NSD" : "")));
183 TH1* h = static_cast<TH1*>(file->Get(hist.Data()));
184 if (!h) {
185 Warning("GetHHD", "Couldn't find HHD histogram %s in %s",
186 hist.Data(), fn);
187 file->Close();
188 savdir->cd();
189 return 0;
190 }
191 TH1* r = static_cast<TH1*>(h->Clone("dndeta_hhd"));
192 r->SetTitle("ALICE Forward (HHD)");
193 r->SetFillStyle(0);
194 r->SetFillColor(0);
195 r->SetMarkerStyle(21);
196 r->SetMarkerColor(kPink+1);
197 r->SetDirectory(0);
198
199 file->Close();
200 savdir->cd();
201 return r;
202 }
203 //__________________________________________________________________
204 THStack* StackResults(Double_t& max)
205 {
206 THStack* stack = new THStack("results", "Stack of Results");
207 max = TMath::Max(max, AddHistogram(stack, fTruth, "e5 p"));
208 max = TMath::Max(max, AddHistogram(stack, fForwardHHD, "", fForwardHHDSym));
209 max = TMath::Max(max, AddHistogram(stack, fForwardMC, "", fForwardMCSym));
210 max = TMath::Max(max, AddHistogram(stack, fCentral, ""));
211 max = TMath::Max(max, AddHistogram(stack, fForward, "", fForwardSym));
212 return stack;
213 }
214 //__________________________________________________________________
215 TMultiGraph* StackOther(Double_t& max) const
216 {
217 gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/OtherData.C");
218 Int_t error = 0;
219 Bool_t onlya = (fShowOthers ? false : true);
220 Int_t trg = (fTrigString ? fTrigString->GetUniqueID() : 0);
221 UShort_t snn = (fSNNString ? fSNNString->GetUniqueID() : 0);
222 Long_t ret = gROOT->ProcessLine(Form("GetData(%d,%d,%d);",
223 snn,trg,onlya));
224 if (error) {
225 Error("StackOther", "Failed to execute GetData(%d,%d,%d)",
226 snn, trg, onlya);
227 return 0;
228 }
229 if (!ret) {
230 Warning("StackOther", "No other data found for sNN=%d, trigger=%d",
231 snn, trg);
232 return 0;
233 }
234 TMultiGraph* other = reinterpret_cast<TMultiGraph*>(ret);
235
236 TGraphAsymmErrors* o = 0;
237 TIter next(other->GetListOfGraphs());
238 while ((o = static_cast<TGraphAsymmErrors*>(next())))
239 max = TMath::Max(max, TMath::MaxElement(o->GetN(), o->GetY()));
240
241 return other;
242 }
243 //__________________________________________________________________
244 THStack* StackRatios(TMultiGraph* others, Double_t& max)
245 {
246 THStack* ratios = new THStack("ratios", "Ratios");
247
248 if (others) {
249 TGraphAsymmErrors* ua5_1 = 0;
250 TGraphAsymmErrors* ua5_2 = 0;
251 TGraphAsymmErrors* alice = 0;
252 TGraphAsymmErrors* cms = 0;
253 TGraphAsymmErrors* o = 0;
254 TIter nextG(others->GetListOfGraphs());
255 while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
256 ratios->Add(Ratio(fForward, o, max));
257 ratios->Add(Ratio(fForwardSym, o, max));
258 ratios->Add(Ratio(fForwardHHD, o, max));
259 ratios->Add(Ratio(fForwardHHDSym, o, max));
260 ratios->Add(Ratio(fCentral, o, max));
261 TString oName(o->GetName());
262 oName.ToLower();
263 if (oName.Contains("ua5")) { if (ua5_1) ua5_2 = o; else ua5_1 = o; }
264 if (oName.Contains("alice")) alice = o;
265 if (oName.Contains("cms")) cms = o;
266 }
267 if (ua5_1 && alice) ratios->Add(Ratio(alice, ua5_1, max));
268 if (ua5_2 && alice) ratios->Add(Ratio(alice, ua5_2, max));
269 if (cms && alice) ratios->Add(Ratio(alice, cms, max));
270 }
271
272 // If we have data from HHD's analysis, then do the ratio of
273 // our result to that data.
274 if (fForwardHHD) {
275 ratios->Add(Ratio(fForward, fForwardHHD, max));
276 ratios->Add(Ratio(fForwardSym, fForwardHHDSym, max));
277 }
278
279 // Do comparison to MC
280 if (fForwardMC) {
281 ratios->Add(Ratio(fForward, fForwardMC, max));
282 ratios->Add(Ratio(fForwardSym, fForwardMCSym, max));
283 }
284
285 // Check if we have ratios
286 if (!ratios->GetHists() ||
287 (ratios->GetHists()->GetEntries() <= 0)) {
288 delete ratios;
289 ratios = 0;
290 }
291 return ratios;
292 }
293 //__________________________________________________________________
294 THStack* StackLeftRight(Double_t& max)
295 {
296 THStack* ret = new THStack("leftright", "Left-right asymmetry");
297 ret->Add(Asymmetry(fForward, max));
298 ret->Add(Asymmetry(fForwardHHD, max));
299 ret->Add(Asymmetry(fForwardMC, max));
300
301 if (!ret->GetHists() ||
302 (ret->GetHists()->GetEntries() <= 0)) {
303 delete ret;
304 ret = 0;
305 }
306 return ret;
307 }
308 //__________________________________________________________________
309 TAxis* FindXAxis(TVirtualPad* p, const char* name)
310 {
311 TObject* o = p->GetListOfPrimitives()->FindObject(name);
312 if (!o) {
313 Warning("FindXAxis", "%s not found in pad", name);
314 return 0;
315 }
316 THStack* stack = dynamic_cast<THStack*>(o);
317 if (!stack) {
318 Warning("FindXAxis", "%s is not a THStack", name);
319 return 0;
320 }
321 if (!stack->GetHistogram()) {
322 Warning("FindXAxis", "%s has no histogram", name);
323 return 0;
324 }
325 TAxis* ret = stack->GetHistogram()->GetXaxis();
326 return ret;
327 }
328 //__________________________________________________________________
329 void Plot(THStack* results,
330 TMultiGraph* others,
331 Double_t max,
332 THStack* ratios,
333 Double_t rmax,
334 THStack* leftright,
335 Double_t amax)
336 {
337 gStyle->SetOptTitle(0);
338 gStyle->SetTitleFont(132, "xyz");
339 gStyle->SetLabelFont(132, "xyz");
340
341 Int_t h = 800;
342 Int_t w = 800; // h / TMath::Sqrt(2);
343 if (!ratios) w *= 1.4;
344 if (!leftright) w *= 1.4;
345 TCanvas* c = new TCanvas("Results", "Results", w, h);
346 c->SetFillColor(0);
347 c->SetBorderSize(0);
348 c->SetBorderMode(0);
349
350 Double_t y1 = 0;
351 Double_t y2 = 0;
352 Double_t y3 = 0;
353 if (ratios) y1 = 0.3;
354 if (leftright) {
355 if (y1 > 0.0001) {
356 y2 = 0.2;
357 y1 = 0.4;
358 }
359 else {
360 y1 = 0.2;
361 y2 = y1;
362 }
363 }
364 PlotResults(results, others, max, y1);
365 c->cd();
366
367 PlotRatios(ratios, rmax, y2, y1);
368 c->cd( );
369
370 PlotLeftRight(leftright, amax, y3, y2);
371 c->cd();
bad9a3c1 372
373
374 Int_t vMin = fVtxAxis->GetXmin();
375 Int_t vMax = fVtxAxis->GetXmax();
376 TString trg(fTrigString->GetTitle());
377 Int_t nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
378 trg = trg.Strip(TString::kBoth);
379 TString base(Form("dndeta_%s_%s_%s_%c%02d%c%02dcm_%09dev",
380 fSysString->GetTitle(),
381 fSNNString->GetTitle(),
382 trg.Data(),
383 vMin < 0 ? 'm' : 'p', TMath::Abs(vMin),
384 vMax < 0 ? 'm' : 'p', TMath::Abs(vMax),
385 nev));
386 c->SaveAs(Form("%s.png", base.Data()));
387 c->SaveAs(Form("%s.root", base.Data()));
388 c->SaveAs(Form("%s.C", base.Data()));
b2e7f2d6 389 }
390 //__________________________________________________________________
391 void PlotResults(THStack* results, TMultiGraph* others,
392 Double_t max, Double_t yd)
393 {
394 // Make a sub-pad for the result itself
395 TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0);
396 p1->SetTopMargin(0.05);
397 p1->SetBorderSize(0);
398 p1->SetBorderMode(0);
399 p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1);
400 p1->SetRightMargin(0.05);
401 p1->SetGridx();
402 p1->SetTicks(1,1);
403 p1->SetNumber(1);
404 p1->Draw();
405 p1->cd();
406
407 results->SetMaximum(1.15*max);
408 results->SetMinimum(yd > 0.00001 ? -0.1 : 0);
409
410 FixAxis(results, 1/(1-yd)/1.7, "#frac{1}{N} #frac{dN_{ch}}{d#eta}");
411
412 p1->Clear();
413 results->DrawClone("nostack e1");
414
415 fRangeParam->fSlave1Axis = results->GetXaxis();
416 fRangeParam->fSlave1Pad = p1;
417
418 // Draw other data
419 if (others) {
420 TGraphAsymmErrors* o = 0;
421 TIter nextG(others->GetListOfGraphs());
422 while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
423 o->DrawClone("same p");
424 }
425
426 // Make a legend in the main result pad
427 TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
428 l->SetNColumns(2);
429 l->SetFillColor(0);
430 l->SetFillStyle(0);
431 l->SetBorderSize(0);
432 l->SetTextFont(132);
433
434 // Put a title on top
435 TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
436 tit->SetNDC();
437 tit->SetTextFont(132);
438 tit->SetTextSize(0.05);
439 tit->Draw();
440
441 // Put a nice label in the plot
442 TString eS;
bad9a3c1 443 UShort_t snn = fSNNString->GetUniqueID();
444 const char* sys = fSysString->GetTitle();
b2e7f2d6 445 if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
446 else eS = Form("%3dGeV", snn);
447 TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s}=%s, %s",
448 sys,
449 eS.Data(),
bad9a3c1 450 fTrigString->GetTitle()));
b2e7f2d6 451 tt->SetNDC();
452 tt->SetTextFont(132);
453 tt->SetTextAlign(33);
454 tt->Draw();
455
456 // Put number of accepted events on the plot
457 Int_t nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
458 TLatex* et = new TLatex(.93, .83, Form("%d events", nev));
459 et->SetNDC();
460 et->SetTextFont(132);
461 et->SetTextAlign(33);
462 et->Draw();
463
464 // Put number of accepted events on the plot
465 if (fVtxAxis) {
466 TLatex* vt = new TLatex(.93, .88, fVtxAxis->GetTitle());
467 vt->SetNDC();
468 vt->SetTextFont(132);
469 vt->SetTextAlign(33);
470 vt->Draw();
471 }
472 // results->Draw("nostack e1 same");
473
474 fRangeParam->fSlave1Axis = FindXAxis(p1, results->GetName());
475 fRangeParam->fSlave1Pad = p1;
476
477
478 // Mark the plot as preliminary
479 TLatex* pt = new TLatex(.12, .93, "Preliminary");
480 pt->SetNDC();
481 pt->SetTextFont(22);
482 pt->SetTextSize(0.07);
483 pt->SetTextColor(kRed+1);
484 pt->SetTextAlign(13);
485 pt->Draw();
486
487 if (!gSystem->AccessPathName("ALICE.png")) {
bad9a3c1 488 TPad* logo = new TPad("logo", "logo", .12, .65, .25, .85, 0, 0, 0);
b2e7f2d6 489 logo->SetFillStyle(0);
490 logo->Draw();
491 logo->cd();
492 TImage* i = TImage::Create();
493 i->ReadImage("ALICE.png");
494 i->Draw();
495 }
496 p1->cd();
497 }
498 //__________________________________________________________________
499 void PlotRatios(THStack* ratios, Double_t max, Double_t y1, Double_t y2)
500 {
501 if (!ratios) return;
502 bool isBottom = (y1 < 0.0001);
503 Double_t yd = y2 - y1;
504 // Make a sub-pad for the result itself
505 TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0);
506 p2->SetTopMargin(0.001);
507 p2->SetRightMargin(0.05);
508 p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
509 p2->SetGridx();
510 p2->SetTicks(1,1);
511 p2->SetNumber(2);
512 p2->Draw();
513 p2->cd();
514
515 // Fix up axis
516 FixAxis(ratios, 1/yd/1.7, "Ratios", 7);
517
518 ratios->SetMaximum(1+TMath::Max(.22,1.05*max));
519 ratios->SetMinimum(1-TMath::Max(.32,1.05*max));
520 p2->Clear();
521 ratios->DrawClone("nostack e1");
522
523
524 // Make a legend
525 TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,
526 isBottom ? .6 : .4);
527 l2->SetNColumns(2);
528 l2->SetFillColor(0);
529 l2->SetFillStyle(0);
530 l2->SetBorderSize(0);
531 l2->SetTextFont(132);
532
533 // Make a nice band from 0.9 to 1.1
534 TGraphErrors* band = new TGraphErrors(2);
535 band->SetPoint(0, fForwardSym->GetXaxis()->GetXmin(), 1);
536 band->SetPoint(1, fForward->GetXaxis()->GetXmax(), 1);
537 band->SetPointError(0, 0, .1);
538 band->SetPointError(1, 0, .1);
539 band->SetFillColor(kYellow+2);
540 band->SetFillStyle(3002);
541 band->SetLineStyle(2);
542 band->SetLineWidth(1);
543 band->Draw("3 same");
544 band->DrawClone("X L same");
545
546 // Replot the ratios on top
547 ratios->DrawClone("nostack e1 same");
548
549 if (isBottom) {
550 fRangeParam->fMasterAxis = FindXAxis(p2, ratios->GetName());
551 p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
552 fRangeParam));
553 }
554 else {
555 fRangeParam->fSlave2Axis = FindXAxis(p2, ratios->GetName());
556 fRangeParam->fSlave2Pad = p2;
557 }
558 }
559 //__________________________________________________________________
560 void PlotLeftRight(THStack* leftright, Double_t max,
561 Double_t y1, Double_t y2)
562 {
563 if (!leftright) return;
564 bool isBottom = (y1 < 0.0001);
565 Double_t yd = y2 - y1;
566 // Make a sub-pad for the result itself
567 TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0);
568 p3->SetTopMargin(0.001);
569 p3->SetRightMargin(0.05);
570 p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
571 p3->SetGridx();
572 p3->SetTicks(1,1);
573 p3->SetNumber(2);
574 p3->Draw();
575 p3->cd();
576
577 TH1* dummy = 0;
578 if (leftright->GetHists()->GetEntries() == 1) {
579 // Add dummy histogram
580 dummy = new TH1F("dummy","", 10, -6, 6);
581 dummy->SetLineColor(0);
582 dummy->SetFillColor(0);
583 dummy->SetMarkerColor(0);
584 leftright->Add(dummy);
585 }
586
587 // Fix up axis
588 FixAxis(leftright, 1/yd/1.7, "Right/Left", 4);
589
590 leftright->SetMaximum(1+TMath::Max(.12,1.05*max));
591 leftright->SetMinimum(1-TMath::Max(.15,1.05*max));
592 p3->Clear();
593 leftright->DrawClone("nostack e1");
594
595
596 // Make a legend
597 TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
598 l2->SetNColumns(2);
599 l2->SetFillColor(0);
600 l2->SetFillStyle(0);
601 l2->SetBorderSize(0);
602 l2->SetTextFont(132);
603#ifndef __CINT__
604 if (dummy) {
605 TList* prims = l2->GetListOfPrimitives();
606 TIter next(prims);
607 TLegendEntry* o = 0;
608 while ((o = static_cast<TLegendEntry*>(next()))) {
609 TString lbl(o->GetLabel());
610 if (lbl != "dummy") continue;
611 prims->Remove(o);
612 break;
613 }
614 }
615#endif
616 // Make a nice band from 0.9 to 1.1
617 TGraphErrors* band = new TGraphErrors(2);
618 band->SetPoint(0, fForwardSym->GetXaxis()->GetXmin(), 1);
619 band->SetPoint(1, fForward->GetXaxis()->GetXmax(), 1);
620 band->SetPointError(0, 0, .05);
621 band->SetPointError(1, 0, .05);
622 band->SetFillColor(kYellow+2);
623 band->SetFillStyle(3002);
624 band->SetLineStyle(2);
625 band->SetLineWidth(1);
626 band->Draw("3 same");
627 band->DrawClone("X L same");
628
629 leftright->DrawClone("nostack e1 same");
630 if (isBottom) {
631 fRangeParam->fMasterAxis = FindXAxis(p3, leftright->GetName());
632 p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
633 fRangeParam));
634 }
635 else {
636 fRangeParam->fSlave2Axis = FindXAxis(p3, leftright->GetName());
637 fRangeParam->fSlave2Pad = p3;
638 }
639 }
640
641 //__________________________________________________________________
642 /**
643 * Fix the apperance of the axis in a stack
644 *
645 * @param stack stack of histogram
646 * @param s Scaling factor
647 * @param ytitle Y axis title
648 * @param force Whether to draw the stack first or not
649 * @param ynDiv Divisions on Y axis
650 */
651 void FixAxis(THStack* stack, Double_t s, const char* ytitle,
652 Int_t ynDiv=210, Bool_t force=true)
653 {
654 if (force) stack->Draw("nostack e1");
655
656 TH1* h = stack->GetHistogram();
657 if (!h) return;
658
659 h->SetXTitle("#eta");
660 h->SetYTitle(ytitle);
661 TAxis* xa = h->GetXaxis();
662 TAxis* ya = h->GetYaxis();
663 if (xa) {
664 xa->SetTitle("#eta");
665 // xa->SetTicks("+-");
666 xa->SetTitleSize(s*xa->GetTitleSize());
667 xa->SetLabelSize(s*xa->GetLabelSize());
668 xa->SetTickLength(s*xa->GetTickLength());
669 }
670 if (ya) {
671 ya->SetTitle(ytitle);
672 ya->SetDecimals();
673 // ya->SetTicks("+-");
674 ya->SetNdivisions(ynDiv);
675 ya->SetTitleSize(s*ya->GetTitleSize());
676 ya->SetTitleOffset(ya->GetTitleOffset()/s);
677 ya->SetLabelSize(s*ya->GetLabelSize());
678 }
679 }
680
681 //__________________________________________________________________
682 Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option) const
683 {
684 // Check if we have input
685 if (!hist) return 0;
686
687 // Rebin if needed
688 Rebin(hist);
689
690 stack->Add(hist, option);
691 return hist->GetMaximum();
692 }
693 //__________________________________________________________________
694 Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option,
695 TH1*& sym) const
696 {
697 // Check if we have input
698 if (!hist) return 0;
699
700 // Rebin if needed
701 Rebin(hist);
702 stack->Add(hist, option);
703
704 // Now symmetrice the histogram
705 sym = Symmetrice(hist);
706 stack->Add(sym, option);
707
708 return hist->GetMaximum();
709 }
710
711 //__________________________________________________________________
712 /**
713 * Rebin a histogram
714 *
715 * @param h Histogram to rebin
716 * @param rebin Rebinning factor
717 *
718 * @return
719 */
720 virtual void Rebin(TH1* h) const
721 {
722 if (fRebin <= 1) return;
723
724 Int_t nBins = h->GetNbinsX();
725 if(nBins % fRebin != 0) {
726 Warning("Rebin", "Rebin factor %d is not a devisor of current number "
727 "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
728 return;
729 }
730
731 // Make a copy
732 TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
733 tmp->Rebin(fRebin);
734 tmp->SetDirectory(0);
735
736 // The new number of bins
737 Int_t nBinsNew = nBins / fRebin;
738 for(Int_t i = 1;i<= nBinsNew; i++) {
739 Double_t content = 0;
740 Double_t sumw = 0;
741 Double_t wsum = 0;
742 Int_t nbins = 0;
743 for(Int_t j = 1; j<=fRebin;j++) {
744 Int_t bin = (i-1)*fRebin + j;
745 Double_t c = h->GetBinContent(bin);
746
747 if (c <= 0) continue;
748
749 if (fCutEdges) {
750 if (h->GetBinContent(bin+1)<=0 ||
751 h->GetBinContent(bin-1)) {
752 Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
753 bin, c, h->GetName(),
754 bin+1, h->GetBinContent(bin+1),
755 bin-1, h->GetBinContent(bin-1));
756 continue;
757 }
758 }
759 Double_t e = h->GetBinError(bin);
760 Double_t w = 1 / (e*e); // 1/c/c
761 content += c;
762 sumw += w;
763 wsum += w * c;
764 nbins++;
765 }
766
767 if(content > 0 && nbins > 0) {
768 tmp->SetBinContent(i, wsum / sumw);
769 tmp->SetBinError(i,1./TMath::Sqrt(sumw));
770 }
771 }
772
773 // Finally, rebin the histogram, and set new content
774 h->Rebin(fRebin);
775 h->Reset();
776 for(Int_t i = 1; i<= nBinsNew; i++) {
777 h->SetBinContent(i,tmp->GetBinContent(i));
778 h->SetBinError(i, tmp->GetBinError(i));
779 }
780
781 delete tmp;
782 }
783 //__________________________________________________________________
784 /**
785 * Make an extension of @a h to make it symmetric about 0
786 *
787 * @param h Histogram to symmertrice
788 *
789 * @return Symmetric extension of @a h
790 */
791 TH1* Symmetrice(const TH1* h) const
792 {
793 Int_t nBins = h->GetNbinsX();
794 TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
795 s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
796 s->Reset();
797 s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
798 s->SetMarkerColor(h->GetMarkerColor());
799 s->SetMarkerSize(h->GetMarkerSize());
800 s->SetMarkerStyle(h->GetMarkerStyle()+4);
801 s->SetFillColor(h->GetFillColor());
802 s->SetFillStyle(h->GetFillStyle());
803 s->SetDirectory(0);
804
805 // Find the first and last bin with data
806 Int_t first = nBins+1;
807 Int_t last = 0;
808 for (Int_t i = 1; i <= nBins; i++) {
809 if (h->GetBinContent(i) <= 0) continue;
810 first = TMath::Min(first, i);
811 last = TMath::Max(last, i);
812 }
813
814 Double_t xfirst = h->GetBinCenter(first-1);
815 Int_t f1 = h->GetXaxis()->FindBin(-xfirst);
816 Int_t l2 = s->GetXaxis()->FindBin(xfirst);
817 for (Int_t i = f1, j=l2; i <= last; i++,j--) {
818 s->SetBinContent(j, h->GetBinContent(i));
819 s->SetBinError(j, h->GetBinError(i));
820 }
821 // Fill in overlap bin
822 s->SetBinContent(l2+1, h->GetBinContent(first));
823 s->SetBinError(l2+1, h->GetBinError(first));
824 return s;
825 }
826 //__________________________________________________________________
827 Double_t RatioMax(TH1* h) const
828 {
829 Int_t nBins = h->GetNbinsX();
830 Double_t ret = 0;
831 for (Int_t i = 1; i <= nBins; i++) {
832 Double_t c = h->GetBinContent(i);
833 if (c == 0) continue;
834 Double_t e = h->GetBinError(i);
835 Double_t d = TMath::Abs(1-c-e);
836 ret = TMath::Max(d, ret);
837 }
838 return ret;
839 }
840 //__________________________________________________________________
841 /**
842 * Compute the ratio of @a h to @a g. @a g is evaluated at the bin
843 * centers of @a h
844 *
845 * @param h Numerator
846 * @param g Divisor
847 *
848 * @return h/g
849 */
850 TH1* Ratio(const TH1* h, const TGraph* g, Double_t& max) const
851 {
852 if (!h || !g) return 0;
853
854 TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
855 ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
856 ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
857 ret->Reset();
858 ret->SetMarkerStyle(g->GetMarkerStyle());
859 ret->SetMarkerColor(h->GetMarkerColor());
860 ret->SetMarkerSize(0.9*g->GetMarkerSize());
861 Double_t xlow = g->GetX()[0];
862 Double_t xhigh = g->GetX()[g->GetN()-1];
863 if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
864
865 for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
866 Double_t c = h->GetBinContent(i);
867 if (c <= 0) continue;
868
869 Double_t x = h->GetBinCenter(i);
870 if (x < xlow || x > xhigh) continue;
871
872 Double_t f = g->Eval(x);
873 if (f <= 0) continue;
874
875 ret->SetBinContent(i, c / f);
876 ret->SetBinError(i, h->GetBinError(i) / f);
877 }
878 if (ret->GetEntries() <= 0) {
879 delete ret;
880 ret = 0;
881 }
882 else
883 max = TMath::Max(RatioMax(ret), max);
884 return ret;
885 }
886 //__________________________________________________________________
887 /**
888 * Make the ratio of h1 to h2
889 *
890 * @param h1 First histogram (numerator)
891 * @param h2 Second histogram (denominator)
892 *
893 * @return h1 / h2
894 */
895 TH1* Ratio(const TH1* h1, const TH1* h2, Double_t& max) const
896 {
897 if (!h1 || !h2) return 0;
898 TH1* t1 = static_cast<TH1*>(h1->Clone(Form("%s_%s",
899 h1->GetName(),
900 h2->GetName())));
901 t1->SetTitle(Form("%s / %s", h1->GetTitle(), h2->GetTitle()));
902 t1->Divide(h2);
903 t1->SetMarkerColor(h1->GetMarkerColor());
904 t1->SetMarkerStyle(h2->GetMarkerStyle());
905 max = TMath::Max(RatioMax(t1), max);
906 return t1;
907 }
908 //__________________________________________________________________
909 /**
910 * Calculate the ratio of two graphs - g1 / g2
911 *
912 * @param g1 Numerator
913 * @param g2 Denominator
914 *
915 * @return g1 / g2 in a histogram
916 */
917 TH1* Ratio(const TGraphAsymmErrors* g1,
918 const TGraphAsymmErrors* g2, Double_t& max) const
919 {
920 Int_t nBins = g1->GetN();
921 TArrayF bins(nBins+1);
922 Double_t dx = 0;
923 for (Int_t i = 0; i < nBins; i++) {
924 Double_t x = g1->GetX()[i];
925 Double_t exl = g1->GetEXlow()[i];
926 Double_t exh = g1->GetEXhigh()[i];
927 bins.fArray[i] = x-exl;
928 bins.fArray[i+1] = x+exh;
929 Double_t dxi = exh+exl;
930 if (i == 0) dx = dxi;
931 else if (dxi != dx) dx = 0;
932 }
933 TString name(Form("%s_%s", g1->GetName(), g2->GetName()));
934 TString title(Form("%s / %s", g1->GetTitle(), g2->GetTitle()));
935 TH1* h = 0;
936 if (dx != 0) {
937 h = new TH1F(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
938 }
939 else {
940 h = new TH1F(name.Data(), title.Data(), nBins, bins.fArray);
941 }
942 h->SetMarkerStyle(g2->GetMarkerStyle());
943 h->SetMarkerColor(g1->GetMarkerColor());
944 h->SetMarkerSize(0.9*g2->GetMarkerSize());
945
946 Double_t low = g2->GetX()[0];
947 Double_t high = g2->GetX()[g2->GetN()-1];
948 if (low > high) { Double_t t = low; low = high; high = t; }
949 for (Int_t i = 0; i < nBins; i++) {
950 Double_t x = g1->GetX()[i];
951 if (x < low-dx || x > high+dx) continue;
952 Double_t c1 = g1->GetY()[i];
953 Double_t e1 = g1->GetErrorY(i);
954 Double_t c2 = g2->Eval(x);
955
956 h->SetBinContent(i+1, c1 / c2);
957 h->SetBinError(i+1, e1 / c2);
958 }
959 max = TMath::Max(RatioMax(h), max);
960 return h;
961 }
962
963 //__________________________________________________________________
964 TH1* Graph2Hist(const TGraphAsymmErrors* g) const
965 {
966 Int_t nBins = g->GetN();
967 TArrayF bins(nBins+1);
968 Double_t dx = 0;
969 for (Int_t i = 0; i < nBins; i++) {
970 Double_t x = g->GetX()[i];
971 Double_t exl = g->GetEXlow()[i];
972 Double_t exh = g->GetEXhigh()[i];
973 bins.fArray[i] = x-exl;
974 bins.fArray[i+1] = x+exh;
975 Double_t dxi = exh+exl;
976 if (i == 0) dx = dxi;
977 else if (dxi != dx) dx = 0;
978 }
979 TString name(g->GetName());
980 TString title(g->GetTitle());
981 TH1D* h = 0;
982 if (dx != 0) {
983 h = new TH1D(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
984 }
985 else {
986 h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
987 }
988 h->SetMarkerStyle(g->GetMarkerStyle());
989 h->SetMarkerColor(g->GetMarkerColor());
990 h->SetMarkerSize(g->GetMarkerSize());
991
992 return h;
993 }
994 //__________________________________________________________________
995 TH1* Asymmetry(TH1* h, Double_t& max)
996 {
997 if (!h) return 0;
998
999 TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName())));
1000 // Int_t oBins = h->GetNbinsX();
1001 // Double_t high = h->GetXaxis()->GetXmax();
1002 // Double_t low = h->GetXaxis()->GetXmin();
1003 // Double_t dBin = (high - low) / oBins;
1004 // Int_t tBins = Int_t(2*high/dBin+.5);
1005 // ret->SetBins(tBins, -high, high);
1006 ret->Reset();
1007 ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
1008 ret->SetYTitle("Right/Left");
1009 Int_t nBins = h->GetNbinsX();
1010 for (Int_t i = 1; i <= nBins; i++) {
1011 Double_t x = h->GetBinCenter(i);
1012 if (x > 0) break;
1013
1014 Double_t c1 = h->GetBinContent(i);
1015 Double_t e1 = h->GetBinError(i);
1016 if (c1 <= 0) continue;
1017
1018 Int_t j = h->FindBin(-x);
1019 if (j <= 0 || j > nBins) continue;
1020
1021 Double_t c2 = h->GetBinContent(j);
1022 Double_t e2 = h->GetBinError(j);
1023
1024 Double_t c12 = c1*c1;
1025 Double_t e = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12));
1026
1027 Int_t k = ret->FindBin(x);
1028 ret->SetBinContent(k, c2/c1);
1029 ret->SetBinError(k, e);
1030 }
1031 max = TMath::Max(max, RatioMax(ret));
1032
1033 return ret;
1034 }
1035
1036
1037 //__________________________________________________________________
1038 Bool_t fShowOthers;
1039 Bool_t fShowAlice;
1040 Bool_t fShowRatios;
1041 Bool_t fShowLeftRight;
1042 UShort_t fRebin;
1043 Bool_t fCutEdges;
1044 TString fTitle;
1045 TString fHHDFile;
1046 TNamed* fTrigString;
1047 TNamed* fSNNString;
1048 TNamed* fSysString;
1049 TAxis* fVtxAxis;
1050 TH1* fForward;
1051 TH1* fForwardMC;
1052 TH1* fForwardHHD;
1053 TH1* fTruth;
1054 TH1* fCentral;
1055 TH1* fForwardSym;
1056 TH1* fForwardMCSym;
1057 TH1* fForwardHHDSym;
1058 TH1* fTriggers;
1059 RangeParam* fRangeParam;
1060
1061};
1062
1063//=== Stuff for auto zooming =========================================
1064void UpdateRange(dNdetaDrawer::RangeParam* p)
1065{
1066 if (!p) {
1067 Warning("UpdateRange", "No parameters %p", p);
1068 return;
1069 }
1070 if (!p->fMasterAxis) {
1071 Warning("UpdateRange", "No master axis %p", p->fMasterAxis);
1072 return;
1073 }
1074 Int_t first = p->fMasterAxis->GetFirst();
1075 Int_t last = p->fMasterAxis->GetLast();
1076 Double_t x1 = p->fMasterAxis->GetBinCenter(first);
1077 Double_t x2 = p->fMasterAxis->GetBinCenter(last);
1078 //Info("UpdateRange", "Range set to [%3d,%3d]->[%f,%f]", first, last, x1,x2);
1079
1080 if (p->fSlave1Axis) {
1081 Int_t i1 = p->fSlave1Axis->FindBin(x1);
1082 Int_t i2 = p->fSlave1Axis->FindBin(x2);
1083 p->fSlave1Axis->SetRange(i1, i2);
1084 p->fSlave1Pad->Modified();
1085 p->fSlave1Pad->Update();
1086 }
1087 if (p->fSlave2Axis) {
1088 Int_t i1 = p->fSlave2Axis->FindBin(x1);
1089 Int_t i2 = p->fSlave2Axis->FindBin(x2);
1090 p->fSlave2Axis->SetRange(i1, i2);
1091 p->fSlave2Pad->Modified();
1092 p->fSlave2Pad->Update();
1093 }
1094 TCanvas* c = gPad->GetCanvas();
1095 c->cd();
1096}
1097
1098//____________________________________________________________________
1099void RangeExec(dNdetaDrawer::RangeParam* p)
1100{
1101 // Event types:
1102 // 51: Mouse move
1103 // 53:
1104 // 1: Button down
1105 // 21: Mouse drag
1106 // 11: Mouse release
1107 // dNdetaDrawer::RangeParam* p =
1108 // reinterpret_cast<dNdetaDrawer::RangeParam*>(addr);
1109 Int_t event = gPad->GetEvent();
1110 TObject *select = gPad->GetSelected();
1111 if (event == 53) {
1112 UpdateRange(p);
1113 return;
1114 }
1115 if (event != 11 || !select || select != p->fMasterAxis) return;
1116 UpdateRange(p);
1117}
1118
1119//=== Steering function ==============================================
1120void
1121DrawdNdeta(const char* filename="forward_dndeta.root",
1122 Int_t flags=0xf,
bad9a3c1 1123 const char* title="",
b2e7f2d6 1124 UShort_t rebin=5,
1125 Bool_t cutEdges=false)
1126{
1127 dNdetaDrawer d;
1128 d.SetRebin(rebin);
1129 d.SetCutEdges(cutEdges);
1130 d.SetTitle(title);
1131 d.SetHHDFile("");
1132 d.SetShowOthers(flags & 0x1);
1133 d.SetShowAlice(flags & 0x2);
1134 d.SetShowRatios(flags & 0x4);
1135 d.SetShowLeftRight(flags & 0x8);
1136 d.Run(filename);
1137}