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