2 * Script to visualise the dN/deta
4 * This script is independent of any AliROOT code - and should stay that way.
8 #include <TGraphErrors.h>
9 #include <TGraphAsymmErrors.h>
10 #include <TMultiGraph.h>
22 #include <TLegendEntry.h>
27 * Class to draw dN/deta results
33 * POD of data for range zooming
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
43 //__________________________________________________________________
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
56 fTrigString(0), // TNamed*
57 fSNNString(0), // TNamed*
58 fSysString(0), // TNamed*
59 fVtxAxis(0), // TAxis*
60 // fForward(0), // TH1*
61 // fForwardMC(0), // TH1*
62 // fForwardHHD(0), // TH1*
64 // fCentral(0), // TH1*
65 // fForwardSym(0), // TH1*
66 // fForwardMCSym(0), // TH1*
67 // fForwardHHDSym(0), // TH1*
72 fRangeParam = new RangeParam;
73 fRangeParam->fMasterAxis = 0;
74 fRangeParam->fSlave1Axis = 0;
75 fRangeParam->fSlave1Pad = 0;
76 fRangeParam->fSlave2Axis = 0;
77 fRangeParam->fSlave2Pad = 0;
79 //==================================================================
82 * @name Set parameters
85 * Show other (UA5, CMS, ...) data
87 * @param x Whether to show or not
89 void SetShowOthers(Bool_t x) { fShowOthers = x; }
90 //__________________________________________________________________
92 * Show ALICE published data
94 * @param x Wheter to show or not
96 void SetShowAlice(Bool_t x) { fShowAlice = x; }
97 //__________________________________________________________________
99 * Whether to show ratios or not. If there's nothing to compare to,
100 * the ratio panel will be implicitly disabled
102 * @param x Whether to show or not
104 void SetShowRatios(Bool_t x) { fShowRatios = x; }
105 //__________________________________________________________________
108 * Whether to show the left/right asymmetry
110 * @param x To show or not
112 void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; }
113 //__________________________________________________________________
115 * Set the rebinning factor
117 * @param x Rebinning factor (must be a divisor in the number of bins)
119 void SetRebin(UShort_t x) { fRebin = x; }
120 //__________________________________________________________________
122 * Wheter to cut away the edges
124 * @param x Whether or not to cut away edges
126 void SetCutEdges(Bool_t x) { fCutEdges = x; }
127 //__________________________________________________________________
129 * Set the title of the plot
133 void SetTitle(TString x) { fTitle = x; }
135 //==================================================================
138 * @name Override settings from input
141 * Override setting from file
143 * @param sNN Center of mass energy per nucleon pair (GeV)
145 void SetSNN(UShort_t sNN)
147 fSNNString = new TNamed("sNN", Form("%04dGeV", sNN));
148 fSNNString->SetUniqueID(sNN);
150 //__________________________________________________________________
152 * Set the collision system
156 * @param sys collision system
158 void SetSys(UShort_t sys)
160 fSysString = new TNamed("sys", (sys == 1 ? "pp" :
161 sys == 2 ? "PbPb" : "unknown"));
162 fSysString->SetUniqueID(sys);
164 //__________________________________________________________________
166 * Set the vertex range in centimeters
168 * @param vzMin Min @f$ v_z@f$
169 * @param vzMax Max @f$ v_z@f$
171 void SetVertexRange(Double_t vzMin, Double_t vzMax)
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));
177 //__________________________________________________________________
178 void SetTrigger(UShort_t trig)
180 fTrigString = new TNamed("trigString", (trig & 0x1 ? "INEL" :
181 trig & 0x2 ? "INEL>0" :
184 fTrigString->SetUniqueID(trig);
188 //==================================================================
191 * @name Main steering functions
194 * Run the code to produce the final result.
196 * @param filename File containing the data
198 void Run(const char* filename="forward_dndeta.root")
200 Double_t max = 0, rmax=0, amax=0;
202 gStyle->SetPalette(1);
204 // --- Open input file -------------------------------------------
205 TFile* file = TFile::Open(filename, "READ");
207 Error("Open", "Cannot open %s", filename);
210 // --- Get forward list ------------------------------------------
211 TList* forward = static_cast<TList*>(file->Get("ForwardResults"));
213 Error("Open", "Couldn't find list ForwardResults");
216 // --- Get information on the run --------------------------------
217 FetchInformation(forward);
218 // --- Set the macro pathand load other data script --------------
219 gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
220 gROOT->GetMacroPath()));
221 gROOT->LoadMacro("OtherData.C");
223 // --- Get the central results -----------------------------------
224 TList* clusters = static_cast<TList*>(file->Get("CentralResults"));
225 if (!clusters) Warning("Open", "Couldn't find list CentralResults");
227 // --- Make our containtes ---------------------------------------
228 fResults = new THStack("results", "Results");
229 fRatios = new THStack("ratios", "Ratios");
230 fLeftRight = new THStack("asymmetry", "Left-right asymmetry");
231 fOthers = new TMultiGraph();
233 // --- Loop over input data --------------------------------------
234 FetchResults(forward, "Forward", max, rmax, amax);
235 FetchResults(clusters, "Central", max, rmax, amax);
237 // --- Get trigger information -----------------------------------
238 TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
240 TList* all = static_cast<TList*>(sums->FindObject("all"));
242 fTriggers = FetchResult(all, "triggers");
243 if (!fTriggers) all->ls();
246 Warning("Run", "List all not found in ForwardSums");
251 Warning("Run", "No ForwardSums directory found in %s", file->GetName());
255 // --- Check our stacks ------------------------------------------
256 if (!fResults->GetHists() ||
257 fResults->GetHists()->GetEntries() <= 0) {
258 Error("Run", "No histograms in result stack!");
261 if (!fOthers->GetListOfGraphs() ||
262 fOthers->GetListOfGraphs()->GetEntries() <= 0) {
263 Warning("Run", "No other data found - disabling that");
266 if (!fRatios->GetHists() ||
267 fRatios->GetHists()->GetEntries() <= 0) {
268 Warning("Run", "No ratio data found - disabling that");
272 if (!fLeftRight->GetHists() ||
273 fLeftRight->GetHists()->GetEntries() <= 0) {
274 Warning("Run", "No left/right data found - disabling that");
276 fShowLeftRight = false;
279 // --- Close the input file --------------------------------------
284 // --- Plot results ----------------------------------------------
285 Plot(max, rmax, amax);
288 //__________________________________________________________________
290 * Fetch the information on the run from the results list
292 * @param results Results list
294 void FetchInformation(const TList* results)
297 fTrigString = static_cast<TNamed*>(results->FindObject("trigString"));
299 fSNNString = static_cast<TNamed*>(results->FindObject("sNN"));
301 fSysString = static_cast<TNamed*>(results->FindObject("sys"));
303 fVtxAxis = static_cast<TAxis*>(results->FindObject("vtxAxis"));
305 fCentAxis = static_cast<TAxis*>(results->FindObject("centAxis"));
306 if (!fTrigString) fTrigString = new TNamed("trigString", "unknown");
307 if (!fSNNString) fSNNString = new TNamed("sNN", "unknown");
308 if (!fSysString) fSysString = new TNamed("sys", "unknown");
310 fVtxAxis = new TAxis(1,0,0);
311 fVtxAxis->SetName("vtxAxis");
312 fVtxAxis->SetTitle("v_{z} range unspecified");
317 Int_t nCent = fCentAxis->GetNbins();
318 centTxt = Form("\n Centrality: %d bins", nCent);
319 for (Int_t i = 0; i <= nCent; i++)
320 centTxt.Append(Form("%c%d", i == 0 ? ' ' : ',',
321 fCentAxis->GetXbins()->At(i)));
323 Info("FetchInformation",
325 " Trigger: %s (%d)\n"
326 " sqrt(sNN): %s (%dGeV)\n"
328 " Vz range: %s (%f,%f)%s",
329 fTrigString->GetTitle(), fTrigString->GetUniqueID(),
330 fSNNString->GetTitle(), fSNNString->GetUniqueID(),
331 fSysString->GetTitle(), fSysString->GetUniqueID(),
332 fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax(),
335 //__________________________________________________________________
336 TMultiGraph* FetchOthers(UShort_t centLow, UShort_t centHigh)
338 TMultiGraph* thisOther = 0;
339 if (!fShowOthers && !fShowAlice) return 0;
341 Bool_t onlya = (fShowOthers ? false : true);
342 UShort_t sys = (fSysString ? fSysString->GetUniqueID() : 0);
343 UShort_t trg = (fTrigString ? fTrigString->GetUniqueID() : 0);
344 UShort_t snn = (fSNNString ? fSNNString->GetUniqueID() : 0);
345 Long_t ret = gROOT->ProcessLine(Form("GetData(%d,%d,%d,%d,%d,%d);",
347 centLow,centHigh,onlya));
349 Warning("FetchResults", "No other data found for sys=%d, sNN=%d, "
350 "trigger=%d %d%%-%d%% central %s",
351 sys, snn, trg, centLow, centHigh,
352 onlya ? " (ALICE results)" : "all");
355 thisOther = reinterpret_cast<TMultiGraph*>(ret);
359 //__________________________________________________________________
361 * Get the results from the top-level list
365 * @param max On return, maximum of data
366 * @param rmax On return, maximum of ratios
367 * @param amax On return, maximum of left-right comparisons
369 void FetchResults(const TList* list,
375 UShort_t n = fCentAxis ? fCentAxis->GetNbins() : 0;
377 TList* all = static_cast<TList*>(list->FindObject("all"));
379 Error("FetchResults", "Couldn't find list 'all' in %s",
382 FetchResults(all, name, FetchOthers(0,0), -1, 0, max, rmax, amax);
386 Int_t nCol = gStyle->GetNumberOfColors();
387 for (UShort_t i = 0; i < n; i++) {
388 UShort_t centLow = fCentAxis->GetXbins()->At(i);
389 UShort_t centHigh = fCentAxis->GetXbins()->At(i+1);
390 TString lname = Form("cent%03d_%03d", centLow, centHigh);
391 TList* thisCent = static_cast<TList*>(list->FindObject(lname));
393 Float_t fc = (centLow+double(centHigh-centLow)/2) / 100;
394 Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5));
395 Int_t col = gStyle->GetColorPalette(icol);
396 Info("FetchResults","Centrality %d-%d color index %d (=%d*%f) -> %d",
397 centLow, centHigh, icol, nCol, fc, col);
399 TString centTxt = Form("%3d%%-%3d%% central", centLow, centHigh);
401 Error("FetchResults", "Couldn't find list '%s' in %s",
402 lname.Data(), list->GetName());
404 FetchResults(thisCent, name, FetchOthers(centLow, centHigh),
405 col, centTxt.Data(), max, rmax, amax);
408 //__________________________________________________________________
409 void SetAttributes(TH1* h, Int_t color)
412 if (color < 0) return;
413 h->SetLineColor(color);
414 h->SetMarkerColor(color);
415 h->SetFillColor(color);
417 //__________________________________________________________________
418 void SetAttributes(TGraph* g, Int_t color)
421 if (color < 0) return;
422 g->SetLineColor(color);
423 g->SetMarkerColor(color);
424 g->SetFillColor(color);
426 //__________________________________________________________________
427 void ModifyTitle(TNamed* h, const char* centTxt)
429 if (!centTxt || !h) return;
430 h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
433 //__________________________________________________________________
435 * Fetch results for a particular centrality bin
439 * @param centLow Low cut on centrality
440 * @param centHigh high cut on centrality
441 * @param max On return, data maximum
442 * @param rmax On return, ratio maximum
443 * @param amax On return, left-right maximum
445 void FetchResults(const TList* list,
447 TMultiGraph* thisOther,
454 TH1* dndeta = FetchResult(list, Form("dndeta%s", name));
455 TH1* dndetaMC = FetchResult(list, Form("dndeta%sMC", name));
456 TH1* dndetaTruth = FetchResult(list, "dndetaTruth");
458 TH1* dndetaMCSym = 0;
459 TH1* tracks = FetchResult(list, "tracks");
460 SetAttributes(dndeta, color);
461 SetAttributes(dndetaMC, color+2);
462 SetAttributes(dndetaTruth,color);
463 SetAttributes(dndetaSym, color);
464 SetAttributes(dndetaMCSym,color+2);
465 SetAttributes(tracks, color+3);
466 ModifyTitle(dndeta, centTxt);
467 ModifyTitle(dndetaMC, centTxt);
468 ModifyTitle(dndetaTruth,centTxt);
469 ModifyTitle(dndetaSym, centTxt);
470 ModifyTitle(dndetaMCSym,centTxt);
471 ModifyTitle(tracks, centTxt);
474 max = TMath::Max(max, AddHistogram(fResults, dndetaTruth, "e5 p"));
475 max = TMath::Max(max, AddHistogram(fResults, dndetaMC, dndetaMCSym));
476 max = TMath::Max(max, AddHistogram(fResults, dndeta, dndetaSym));
477 max = TMath::Max(max, AddHistogram(fResults, tracks));
479 // Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f",
480 // dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
482 if (fShowLeftRight) {
483 fLeftRight->Add(Asymmetry(dndeta, amax));
484 fLeftRight->Add(Asymmetry(dndetaMC, amax));
485 fLeftRight->Add(Asymmetry(tracks, amax));
489 TIter next(thisOther->GetListOfGraphs());
491 while ((g = static_cast<TGraph*>(next()))) {
492 fRatios->Add(Ratio(dndeta, g, rmax));
493 fRatios->Add(Ratio(dndetaSym, g, rmax));
494 SetAttributes(g, color);
495 ModifyTitle(g, centTxt);
496 if (!fOthers->GetListOfGraphs() ||
497 !fOthers->GetListOfGraphs()->FindObject(g->GetName()))
500 // fOthers->Add(thisOther);
503 if (!fRatios->GetHists()->FindObject(tracks->GetName()))
504 fRatios->Add(Ratio(dndeta, tracks, rmax));
508 fRatios->Add(Ratio(dndeta, dndetaTruth, rmax));
509 fRatios->Add(Ratio(dndetaSym, dndetaTruth, rmax));
510 fRatios->Add(Ratio(dndetaMC, dndetaTruth, rmax));
511 fRatios->Add(Ratio(dndetaMCSym, dndetaTruth, rmax));
514 //__________________________________________________________________
517 * @param max Max value
518 * @param rmax Maximum diviation from 1 of ratios
519 * @param amax Maximum diviation from 1 of asymmetries
521 void Plot(Double_t max,
525 gStyle->SetOptTitle(0);
526 gStyle->SetTitleFont(132, "xyz");
527 gStyle->SetLabelFont(132, "xyz");
530 Int_t w = 800; // h / TMath::Sqrt(2);
534 if (!fShowRatios) w *= 1.4;
536 if (!fShowLeftRight) w *= 1.4;
539 y1 = (y11 > 0.0001 ? 0.4 : 0.2);
540 y2 = (y11 > 0.0001 ? y11 : 0.2);
542 TCanvas* c = new TCanvas("Results", "Results", w, h);
548 Info("Plot", "y1=%f, y2=%f, y3=%f extra: %s %s",
549 y1, y2, y2, (fShowRatios ? "ratios" : ""),
550 (fShowLeftRight ? "right/left" : ""));
552 PlotResults(max, y1);
555 PlotRatios(rmax, y2, y1);
558 PlotLeftRight(amax, y3, y2);
562 Int_t vMin = fVtxAxis->GetXmin();
563 Int_t vMax = fVtxAxis->GetXmax();
564 TString trg(fTrigString->GetTitle());
566 if (fTriggers) nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
567 trg = trg.Strip(TString::kBoth);
568 TString base(Form("dndeta_%s_%s_%s_%c%02d%c%02dcm_%09dev",
569 fSysString->GetTitle(),
570 fSNNString->GetTitle(),
572 vMin < 0 ? 'm' : 'p', TMath::Abs(vMin),
573 vMax < 0 ? 'm' : 'p', TMath::Abs(vMax),
575 c->SaveAs(Form("%s.png", base.Data()));
576 c->SaveAs(Form("%s.root", base.Data()));
577 c->SaveAs(Form("%s.C", base.Data()));
579 //__________________________________________________________________
588 void BuildLegend(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
590 TLegend* l = new TLegend(x1,y1,x2,y2);
597 TIter next(fResults->GetHists());
599 while ((hist = next())) {
600 TString n(hist->GetTitle());
601 if (n.Contains("mirrored")) continue;
602 l->AddEntry(hist, hist->GetTitle(), "pl");
604 TIter nexto(fOthers->GetListOfGraphs());
605 while ((hist = nexto())) {
606 TString n(hist->GetTitle());
607 if (n.Contains("mirrored")) continue;
608 l->AddEntry(hist, hist->GetTitle(), "pl");
610 TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp");
611 d1->SetLineColor(kBlack);
612 d1->SetMarkerColor(kBlack);
613 d1->SetMarkerStyle(20);
614 TLegendEntry* d2 = l->AddEntry("d2", "Mirrored data", "lp");
615 d2->SetLineColor(kBlack);
616 d2->SetMarkerColor(kBlack);
617 d2->SetMarkerStyle(24);
621 //__________________________________________________________________
626 * @param yd Bottom position of pad
628 void PlotResults(Double_t max, Double_t yd)
630 // Make a sub-pad for the result itself
631 TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0);
632 p1->SetTopMargin(0.05);
633 p1->SetBorderSize(0);
634 p1->SetBorderMode(0);
635 p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1);
636 p1->SetRightMargin(0.05);
643 // Info("PlotResults", "Plotting results with max=%f", max);
644 fResults->SetMaximum(1.15*max);
645 fResults->SetMinimum(yd > 0.00001 ? -0.1 : 0);
647 FixAxis(fResults, 1/(1-yd)/1.7, "#frac{1}{N} #frac{dN_{ch}}{d#eta}");
650 fResults->DrawClone("nostack e1");
652 fRangeParam->fSlave1Axis = fResults->GetXaxis();
653 fRangeParam->fSlave1Pad = p1;
656 if (fShowOthers || fShowAlice) {
657 TGraphAsymmErrors* o = 0;
658 TIter nextG(fOthers->GetListOfGraphs());
659 while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
660 o->DrawClone("same p");
663 // Make a legend in the main result pad
664 BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
666 TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
674 // Put a title on top
675 TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
677 tit->SetTextFont(132);
678 tit->SetTextSize(0.05);
681 // Put a nice label in the plot
683 UShort_t snn = fSNNString->GetUniqueID();
684 const char* sys = fSysString->GetTitle();
685 if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
686 else eS = Form("%3dGeV", snn);
687 TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s}=%s, %s",
690 fTrigString->GetTitle()));
692 tt->SetTextFont(132);
693 tt->SetTextAlign(33);
696 // Put number of accepted events on the plot
698 if (fTriggers) nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
699 TLatex* et = new TLatex(.93, .83, Form("%d events", nev));
701 et->SetTextFont(132);
702 et->SetTextAlign(33);
705 // Put number of accepted events on the plot
707 TLatex* vt = new TLatex(.93, .88, fVtxAxis->GetTitle());
709 vt->SetTextFont(132);
710 vt->SetTextAlign(33);
713 // results->Draw("nostack e1 same");
715 fRangeParam->fSlave1Axis = FindXAxis(p1, fResults->GetName());
716 fRangeParam->fSlave1Pad = p1;
719 // Mark the plot as preliminary
720 TLatex* pt = new TLatex(.12, .93, "Preliminary");
723 pt->SetTextSize(0.07);
724 pt->SetTextColor(kRed+1);
725 pt->SetTextAlign(13);
728 if (!gSystem->AccessPathName("ALICE.png")) {
729 TPad* logo = new TPad("logo", "logo", .12, .65, .25, .85, 0, 0, 0);
730 logo->SetFillStyle(0);
733 TImage* i = TImage::Create();
734 i->ReadImage("ALICE.png");
739 //__________________________________________________________________
743 * @param max Maximum diviation from 1
744 * @param y1 Lower y coordinate of pad
745 * @param y2 Upper y coordinate of pad
747 void PlotRatios(Double_t max, Double_t y1, Double_t y2)
749 if (!fShowRatios) return;
751 bool isBottom = (y1 < 0.0001);
752 Double_t yd = y2 - y1;
753 // Make a sub-pad for the result itself
754 TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0);
755 p2->SetTopMargin(0.001);
756 p2->SetRightMargin(0.05);
757 p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
765 FixAxis(fRatios, 1/yd/1.7, "Ratios", 7);
767 fRatios->SetMaximum(1+TMath::Max(.22,1.05*max));
768 fRatios->SetMinimum(1-TMath::Max(.32,1.05*max));
770 fRatios->DrawClone("nostack e1");
774 TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,
779 l2->SetBorderSize(0);
780 l2->SetTextFont(132);
782 // Make a nice band from 0.9 to 1.1
783 TGraphErrors* band = new TGraphErrors(2);
784 band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
785 band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
786 band->SetPointError(0, 0, .1);
787 band->SetPointError(1, 0, .1);
788 band->SetFillColor(kYellow+2);
789 band->SetFillStyle(3002);
790 band->SetLineStyle(2);
791 band->SetLineWidth(1);
792 band->Draw("3 same");
793 band->DrawClone("X L same");
795 // Replot the ratios on top
796 fRatios->DrawClone("nostack e1 same");
799 fRangeParam->fMasterAxis = FindXAxis(p2, fRatios->GetName());
800 p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
804 fRangeParam->fSlave2Axis = FindXAxis(p2, fRatios->GetName());
805 fRangeParam->fSlave2Pad = p2;
808 //__________________________________________________________________
810 * Plot the asymmetries
812 * @param max Maximum diviation from 1
813 * @param y1 Lower y coordinate of pad
814 * @param y2 Upper y coordinate of pad
816 void PlotLeftRight(Double_t max, Double_t y1, Double_t y2)
818 if (!fShowLeftRight) return;
820 bool isBottom = (y1 < 0.0001);
821 Double_t yd = y2 - y1;
822 // Make a sub-pad for the result itself
823 TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0);
824 p3->SetTopMargin(0.001);
825 p3->SetRightMargin(0.05);
826 p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
834 if (fLeftRight->GetHists()->GetEntries() == 1) {
835 // Add dummy histogram
836 dummy = new TH1F("dummy","", 10, -6, 6);
837 dummy->SetLineColor(0);
838 dummy->SetFillColor(0);
839 dummy->SetMarkerColor(0);
840 fLeftRight->Add(dummy);
844 FixAxis(fLeftRight, 1/yd/1.7, "Right/Left", 4);
846 fLeftRight->SetMaximum(1+TMath::Max(.12,1.05*max));
847 fLeftRight->SetMinimum(1-TMath::Max(.15,1.05*max));
849 fLeftRight->DrawClone("nostack e1");
853 TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
857 l2->SetBorderSize(0);
858 l2->SetTextFont(132);
861 TList* prims = l2->GetListOfPrimitives();
864 while ((o = static_cast<TLegendEntry*>(next()))) {
865 TString lbl(o->GetLabel());
866 if (lbl != "dummy") continue;
872 // Make a nice band from 0.9 to 1.1
873 TGraphErrors* band = new TGraphErrors(2);
874 band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
875 band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
876 band->SetPointError(0, 0, .05);
877 band->SetPointError(1, 0, .05);
878 band->SetFillColor(kYellow+2);
879 band->SetFillStyle(3002);
880 band->SetLineStyle(2);
881 band->SetLineWidth(1);
882 band->Draw("3 same");
883 band->DrawClone("X L same");
885 fLeftRight->DrawClone("nostack e1 same");
887 fRangeParam->fMasterAxis = FindXAxis(p3, fLeftRight->GetName());
888 p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
892 fRangeParam->fSlave2Axis = FindXAxis(p3, fLeftRight->GetName());
893 fRangeParam->fSlave2Pad = p3;
897 //==================================================================
900 * @name Data utility functions
903 * Get a result from the passed list
905 * @param list List to search
906 * @param name Object name to search for
910 TH1* FetchResult(const TList* list, const char* name) const
914 TH1* ret = static_cast<TH1*>(list->FindObject(name));
917 Warning("GetResult", "Histogram %s not found", name);
922 //__________________________________________________________________
924 * Add a histogram to the stack after possibly rebinning it
926 * @param stack Stack to add to
927 * @param hist histogram
928 * @param option Draw options
930 * @return Maximum of histogram
932 Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option="") const
934 // Check if we have input
940 // Info("AddHistogram", "Adding %s to %s",
941 // hist->GetName(), stack->GetName());
942 stack->Add(hist, option);
944 return hist->GetMaximum();
946 //__________________________________________________________________
948 * Add a histogram to the stack after possibly rebinning it
950 * @param stack Stack to add to
951 * @param hist histogram
952 * @param option Draw options
953 * @param sym On return, the data symmetriced (added to stack)
955 * @return Maximum of histogram
957 Double_t AddHistogram(THStack* stack, TH1* hist, TH1*& sym,
958 Option_t* option="") const
960 // Check if we have input
965 stack->Add(hist, option);
967 // Now symmetrice the histogram
968 sym = Symmetrice(hist);
969 stack->Add(sym, option);
971 // Info("AddHistogram", "Adding %s and %s to %s",
972 // hist->GetName(), sym->GetName(), stack->GetName());
973 return hist->GetMaximum();
976 //__________________________________________________________________
980 * @param h Histogram to rebin
981 * @param rebin Rebinning factor
985 virtual void Rebin(TH1* h) const
987 if (fRebin <= 1) return;
989 Int_t nBins = h->GetNbinsX();
990 if(nBins % fRebin != 0) {
991 Warning("Rebin", "Rebin factor %d is not a devisor of current number "
992 "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
997 TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
999 tmp->SetDirectory(0);
1001 // The new number of bins
1002 Int_t nBinsNew = nBins / fRebin;
1003 for(Int_t i = 1;i<= nBinsNew; i++) {
1004 Double_t content = 0;
1008 for(Int_t j = 1; j<=fRebin;j++) {
1009 Int_t bin = (i-1)*fRebin + j;
1010 Double_t c = h->GetBinContent(bin);
1012 if (c <= 0) continue;
1015 if (h->GetBinContent(bin+1)<=0 ||
1016 h->GetBinContent(bin-1)) {
1017 Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
1018 bin, c, h->GetName(),
1019 bin+1, h->GetBinContent(bin+1),
1020 bin-1, h->GetBinContent(bin-1));
1024 Double_t e = h->GetBinError(bin);
1025 Double_t w = 1 / (e*e); // 1/c/c
1032 if(content > 0 && nbins > 1 ) {
1033 tmp->SetBinContent(i, wsum / sumw);
1034 tmp->SetBinError(i,1./TMath::Sqrt(sumw));
1038 // Finally, rebin the histogram, and set new content
1041 for(Int_t i = 1; i<= nBinsNew; i++) {
1042 h->SetBinContent(i,tmp->GetBinContent(i));
1043 h->SetBinError(i, tmp->GetBinError(i));
1048 //__________________________________________________________________
1050 * Make an extension of @a h to make it symmetric about 0
1052 * @param h Histogram to symmertrice
1054 * @return Symmetric extension of @a h
1056 TH1* Symmetrice(const TH1* h) const
1058 Int_t nBins = h->GetNbinsX();
1059 TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
1060 s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
1062 s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
1063 s->SetMarkerColor(h->GetMarkerColor());
1064 s->SetMarkerSize(h->GetMarkerSize());
1065 s->SetMarkerStyle(h->GetMarkerStyle()+4);
1066 s->SetFillColor(h->GetFillColor());
1067 s->SetFillStyle(h->GetFillStyle());
1070 // Find the first and last bin with data
1071 Int_t first = nBins+1;
1073 for (Int_t i = 1; i <= nBins; i++) {
1074 if (h->GetBinContent(i) <= 0) continue;
1075 first = TMath::Min(first, i);
1076 last = TMath::Max(last, i);
1079 Double_t xfirst = h->GetBinCenter(first-1);
1080 Int_t f1 = h->GetXaxis()->FindBin(-xfirst);
1081 Int_t l2 = s->GetXaxis()->FindBin(xfirst);
1082 for (Int_t i = f1, j=l2; i <= last; i++,j--) {
1083 s->SetBinContent(j, h->GetBinContent(i));
1084 s->SetBinError(j, h->GetBinError(i));
1086 // Fill in overlap bin
1087 s->SetBinContent(l2+1, h->GetBinContent(first));
1088 s->SetBinError(l2+1, h->GetBinError(first));
1091 //__________________________________________________________________
1093 * Calculate the left-right asymmetry of input histogram
1095 * @param h Input histogram
1096 * @param max On return, the maximum distance from 1 of the histogram
1100 TH1* Asymmetry(TH1* h, Double_t& max)
1104 TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName())));
1105 // Int_t oBins = h->GetNbinsX();
1106 // Double_t high = h->GetXaxis()->GetXmax();
1107 // Double_t low = h->GetXaxis()->GetXmin();
1108 // Double_t dBin = (high - low) / oBins;
1109 // Int_t tBins = Int_t(2*high/dBin+.5);
1110 // ret->SetBins(tBins, -high, high);
1111 ret->SetDirectory(0);
1113 ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
1114 ret->SetYTitle("Right/Left");
1115 Int_t nBins = h->GetNbinsX();
1116 for (Int_t i = 1; i <= nBins; i++) {
1117 Double_t x = h->GetBinCenter(i);
1120 Double_t c1 = h->GetBinContent(i);
1121 Double_t e1 = h->GetBinError(i);
1122 if (c1 <= 0) continue;
1124 Int_t j = h->FindBin(-x);
1125 if (j <= 0 || j > nBins) continue;
1127 Double_t c2 = h->GetBinContent(j);
1128 Double_t e2 = h->GetBinError(j);
1130 Double_t c12 = c1*c1;
1131 Double_t e = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12));
1133 Int_t k = ret->FindBin(x);
1134 ret->SetBinContent(k, c2/c1);
1135 ret->SetBinError(k, e);
1137 max = TMath::Max(max, RatioMax(ret));
1141 //__________________________________________________________________
1143 * Transform a graph into a histogram
1149 TH1* Graph2Hist(const TGraphAsymmErrors* g) const
1151 Int_t nBins = g->GetN();
1152 TArrayF bins(nBins+1);
1154 for (Int_t i = 0; i < nBins; i++) {
1155 Double_t x = g->GetX()[i];
1156 Double_t exl = g->GetEXlow()[i];
1157 Double_t exh = g->GetEXhigh()[i];
1158 bins.fArray[i] = x-exl;
1159 bins.fArray[i+1] = x+exh;
1160 Double_t dxi = exh+exl;
1161 if (i == 0) dx = dxi;
1162 else if (dxi != dx) dx = 0;
1164 TString name(g->GetName());
1165 TString title(g->GetTitle());
1168 h = new TH1D(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
1171 h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
1173 h->SetMarkerStyle(g->GetMarkerStyle());
1174 h->SetMarkerColor(g->GetMarkerColor());
1175 h->SetMarkerSize(g->GetMarkerSize());
1180 //==================================================================
1183 * @name Ratio utility functions
1186 * Get the maximum diviation from 1 in the passed ratio
1188 * @param h Ratio histogram
1190 * @return Max diviation from 1
1192 Double_t RatioMax(TH1* h) const
1194 Int_t nBins = h->GetNbinsX();
1196 for (Int_t i = 1; i <= nBins; i++) {
1197 Double_t c = h->GetBinContent(i);
1198 if (c == 0) continue;
1199 Double_t e = h->GetBinError(i);
1200 Double_t d = TMath::Abs(1-c-e);
1201 ret = TMath::Max(d, ret);
1205 //__________________________________________________________________
1207 * Make the ratio of h1 to h2
1209 * @param h1 First histogram (numerator)
1210 * @param h2 Second histogram (denominator)
1214 TH1* Ratio(const TObject* o1, const TObject* o2, Double_t& max) const
1216 const TH1* h1 = dynamic_cast<const TH1*>(o1);
1218 const TH1* h2 = dynamic_cast<const TH1*>(o2);
1219 if (h2) return Ratio(h1,h2,max);
1221 const TGraph* g2 = dynamic_cast<const TGraph*>(o2);
1222 if (g2) return Ratio(h1,g2,max);
1225 const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(o1);
1227 const TGraphAsymmErrors* g2 = dynamic_cast<const TGraphAsymmErrors*>(o2);
1228 if (g2) return Ratio(g1, g2, max);
1231 Warning("Ratio", "Don't know how to divide a %s (%s) with a %s (%s)",
1232 o1->ClassName(), o1->GetName(), o2->ClassName(), o2->GetName());
1235 //__________________________________________________________________
1237 * Compute the ratio of @a h to @a g. @a g is evaluated at the bin
1240 * @param h Numerator
1245 TH1* Ratio(const TH1* h, const TGraph* g, Double_t& max) const
1247 if (!h || !g) return 0;
1249 TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
1250 ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
1251 ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
1253 ret->SetMarkerStyle(g->GetMarkerStyle());
1254 ret->SetMarkerColor(h->GetMarkerColor());
1255 ret->SetMarkerSize(0.9*g->GetMarkerSize());
1256 ret->SetDirectory(0);
1257 Double_t xlow = g->GetX()[0];
1258 Double_t xhigh = g->GetX()[g->GetN()-1];
1259 if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
1261 for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1262 Double_t c = h->GetBinContent(i);
1263 if (c <= 0) continue;
1265 Double_t x = h->GetBinCenter(i);
1266 if (x < xlow || x > xhigh) continue;
1268 Double_t f = g->Eval(x);
1269 if (f <= 0) continue;
1271 ret->SetBinContent(i, c / f);
1272 ret->SetBinError(i, h->GetBinError(i) / f);
1274 if (ret->GetEntries() <= 0) {
1279 max = TMath::Max(RatioMax(ret), max);
1282 //__________________________________________________________________
1284 * Make the ratio of h1 to h2
1286 * @param h1 First histogram (numerator)
1287 * @param h2 Second histogram (denominator)
1291 TH1* Ratio(const TH1* h1, const TH1* h2, Double_t& max) const
1293 if (!h1 || !h2) return 0;
1294 TH1* t1 = static_cast<TH1*>(h1->Clone(Form("%s_%s",
1297 t1->SetTitle(Form("%s / %s", h1->GetTitle(), h2->GetTitle()));
1299 t1->SetMarkerColor(h1->GetMarkerColor());
1300 t1->SetMarkerStyle(h2->GetMarkerStyle());
1301 t1->SetDirectory(0);
1302 max = TMath::Max(RatioMax(t1), max);
1305 //__________________________________________________________________
1307 * Calculate the ratio of two graphs - g1 / g2
1309 * @param g1 Numerator
1310 * @param g2 Denominator
1312 * @return g1 / g2 in a histogram
1314 TH1* Ratio(const TGraphAsymmErrors* g1,
1315 const TGraphAsymmErrors* g2, Double_t& max) const
1317 Int_t nBins = g1->GetN();
1318 TArrayF bins(nBins+1);
1320 for (Int_t i = 0; i < nBins; i++) {
1321 Double_t x = g1->GetX()[i];
1322 Double_t exl = g1->GetEXlow()[i];
1323 Double_t exh = g1->GetEXhigh()[i];
1324 bins.fArray[i] = x-exl;
1325 bins.fArray[i+1] = x+exh;
1326 Double_t dxi = exh+exl;
1327 if (i == 0) dx = dxi;
1328 else if (dxi != dx) dx = 0;
1330 TString name(Form("%s_%s", g1->GetName(), g2->GetName()));
1331 TString title(Form("%s / %s", g1->GetTitle(), g2->GetTitle()));
1334 h = new TH1F(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
1337 h = new TH1F(name.Data(), title.Data(), nBins, bins.fArray);
1339 h->SetMarkerStyle(g2->GetMarkerStyle());
1340 h->SetMarkerColor(g1->GetMarkerColor());
1341 h->SetMarkerSize(0.9*g2->GetMarkerSize());
1344 Double_t low = g2->GetX()[0];
1345 Double_t high = g2->GetX()[g2->GetN()-1];
1346 if (low > high) { Double_t t = low; low = high; high = t; }
1347 for (Int_t i = 0; i < nBins; i++) {
1348 Double_t x = g1->GetX()[i];
1349 if (x < low-dx || x > high+dx) continue;
1350 Double_t c1 = g1->GetY()[i];
1351 Double_t e1 = g1->GetErrorY(i);
1352 Double_t c2 = g2->Eval(x);
1354 h->SetBinContent(i+1, c1 / c2);
1355 h->SetBinError(i+1, e1 / c2);
1357 max = TMath::Max(RatioMax(h), max);
1361 //==================================================================
1364 * @name Graphics utility functions
1367 * Find an X axis in a pad
1370 * @param name Histogram to find axis for
1372 * @return Found axis or null
1374 TAxis* FindXAxis(TVirtualPad* p, const char* name)
1376 TObject* o = p->GetListOfPrimitives()->FindObject(name);
1378 Warning("FindXAxis", "%s not found in pad", name);
1381 THStack* stack = dynamic_cast<THStack*>(o);
1383 Warning("FindXAxis", "%s is not a THStack", name);
1386 if (!stack->GetHistogram()) {
1387 Warning("FindXAxis", "%s has no histogram", name);
1390 TAxis* ret = stack->GetHistogram()->GetXaxis();
1394 //__________________________________________________________________
1396 * Fix the apperance of the axis in a stack
1398 * @param stack stack of histogram
1399 * @param s Scaling factor
1400 * @param ytitle Y axis title
1401 * @param force Whether to draw the stack first or not
1402 * @param ynDiv Divisions on Y axis
1404 void FixAxis(THStack* stack, Double_t s, const char* ytitle,
1405 Int_t ynDiv=210, Bool_t force=true)
1408 Warning("FixAxis", "No stack passed for %s!", ytitle);
1411 if (force) stack->Draw("nostack e1");
1413 TH1* h = stack->GetHistogram();
1415 Warning("FixAxis", "Stack %s has no histogram", stack->GetName());
1419 h->SetXTitle("#eta");
1420 h->SetYTitle(ytitle);
1421 TAxis* xa = h->GetXaxis();
1422 TAxis* ya = h->GetYaxis();
1424 xa->SetTitle("#eta");
1425 // xa->SetTicks("+-");
1426 xa->SetTitleSize(s*xa->GetTitleSize());
1427 xa->SetLabelSize(s*xa->GetLabelSize());
1428 xa->SetTickLength(s*xa->GetTickLength());
1430 if (stack != fResults) {
1431 TAxis* rxa = fResults->GetXaxis();
1432 xa->Set(rxa->GetNbins(), rxa->GetXmin(), rxa->GetXmax());
1436 ya->SetTitle(ytitle);
1438 // ya->SetTicks("+-");
1439 ya->SetNdivisions(ynDiv);
1440 ya->SetTitleSize(s*ya->GetTitleSize());
1441 ya->SetTitleOffset(ya->GetTitleOffset()/s);
1442 ya->SetLabelSize(s*ya->GetLabelSize());
1449 //__________________________________________________________________
1450 Bool_t fShowOthers; // Show other data
1451 Bool_t fShowAlice; // Show ALICE published data
1452 Bool_t fShowRatios; // Show ratios
1453 Bool_t fShowLeftRight;// Show asymmetry
1454 UShort_t fRebin; // Rebinning factor
1455 Bool_t fCutEdges; // Whether to cut edges
1456 TString fTitle; // Title on plot
1457 TNamed* fTrigString; // Trigger string (read, or set)
1458 TNamed* fSNNString; // Energy string (read, or set)
1459 TNamed* fSysString; // Collision system string (read or set)
1460 TAxis* fVtxAxis; // Vertex cuts (read or set)
1461 TAxis* fCentAxis; // Centrality axis
1462 THStack* fResults; // Stack of results
1463 THStack* fRatios; // Stack of ratios
1464 THStack* fLeftRight; // Left-right asymmetry
1465 TMultiGraph* fOthers; // Older data
1466 TH1* fTriggers; // Number of triggers
1467 RangeParam* fRangeParam; // Parameter object for range zoom
1471 //=== Stuff for auto zooming =========================================
1472 void UpdateRange(dNdetaDrawer::RangeParam* p)
1475 Warning("UpdateRange", "No parameters %p", p);
1478 if (!p->fMasterAxis) {
1479 Warning("UpdateRange", "No master axis %p", p->fMasterAxis);
1482 Int_t first = p->fMasterAxis->GetFirst();
1483 Int_t last = p->fMasterAxis->GetLast();
1484 Double_t x1 = p->fMasterAxis->GetBinCenter(first);
1485 Double_t x2 = p->fMasterAxis->GetBinCenter(last);
1486 //Info("UpdateRange", "Range set to [%3d,%3d]->[%f,%f]", first, last, x1,x2);
1488 if (p->fSlave1Axis) {
1489 Int_t i1 = p->fSlave1Axis->FindBin(x1);
1490 Int_t i2 = p->fSlave1Axis->FindBin(x2);
1491 p->fSlave1Axis->SetRange(i1, i2);
1492 p->fSlave1Pad->Modified();
1493 p->fSlave1Pad->Update();
1495 if (p->fSlave2Axis) {
1496 Int_t i1 = p->fSlave2Axis->FindBin(x1);
1497 Int_t i2 = p->fSlave2Axis->FindBin(x2);
1498 p->fSlave2Axis->SetRange(i1, i2);
1499 p->fSlave2Pad->Modified();
1500 p->fSlave2Pad->Update();
1502 TCanvas* c = gPad->GetCanvas();
1506 //____________________________________________________________________
1507 void RangeExec(dNdetaDrawer::RangeParam* p)
1514 // 11: Mouse release
1515 // dNdetaDrawer::RangeParam* p =
1516 // reinterpret_cast<dNdetaDrawer::RangeParam*>(addr);
1517 Int_t event = gPad->GetEvent();
1518 TObject *select = gPad->GetSelected();
1523 if (event != 11 || !select || select != p->fMasterAxis) return;
1527 //=== Steering function ==============================================
1529 DrawdNdeta(const char* filename="forward_dndeta.root",
1531 const char* title="",
1533 Bool_t cutEdges=false,
1542 d.SetCutEdges(cutEdges);
1544 d.SetShowOthers(flags & 0x1);
1545 d.SetShowAlice(flags & 0x2);
1546 d.SetShowRatios(flags & 0x4);
1547 d.SetShowLeftRight(flags & 0x8);
1548 // Do the below if your input data does not contain these settings
1549 if (sNN > 0) d.SetSNN(sNN); // Collision energy per nucleon pair (GeV)
1550 if (sys > 0) d.SetSys(sys); // Collision system (1:pp, 2:PbPB)
1551 if (trg > 0) d.SetTrigger(trg); // Collision trigger (1:INEL, 2:INEL>0, 4:NSD)
1552 if (vzMin < 999 && vzMax > -999)
1553 d.SetVertexRange(vzMin,vzMax); // Collision vertex range (cm)
1556 //____________________________________________________________________