3 * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
4 * @date Wed Mar 23 14:07:10 2011
6 * @brief Script to visualise the dN/deta for pp and PbPb
8 * This script is independent of any AliROOT code - and should stay
11 * The script is <emph>very</emph> long - sigh - the joy of drawing
12 * things nicely in ROOT
14 * @ingroup pwg2_forward_dndeta
19 #include <TGraphErrors.h>
20 #include <TGraphAsymmErrors.h>
21 #include <TMultiGraph.h>
33 #include <TLegendEntry.h>
37 Double_t myFunc(Double_t* xp, Double_t* pp);
40 * Class to draw dN/deta results
42 * @ingroup pwg2_forward_tasks_dndeta
43 * @ingroup pwg2_forward_dndeta
48 * POD of data for range zooming
52 TAxis* fMasterAxis; // Master axis
53 TAxis* fSlave1Axis; // First slave axis
54 TVirtualPad* fSlave1Pad; // First slave pad
55 TAxis* fSlave2Axis; // Second slave axis
56 TVirtualPad* fSlave2Pad; // Second slave pad
58 //__________________________________________________________________
64 : fShowOthers(0), // Bool_t
65 fShowAlice(false), // Bool_t
66 fShowRatios(false), // Bool_t
67 fShowLeftRight(false), // Bool_t
68 fShowRings(false), // Bool_t
69 fRebin(5), // UShort_t
70 fCutEdges(false), // Bool_t
71 fTitle(""), // TString
72 fTrigString(0), // TNamed*
73 fNormString(0), // TNamed*
74 fSNNString(0), // TNamed*
75 fSysString(0), // TNamed*
76 fVtxAxis(0), // TAxis*
77 fCentAxis(0), // TAxis*
86 fRangeParam = new RangeParam;
87 fRangeParam->fMasterAxis = 0;
88 fRangeParam->fSlave1Axis = 0;
89 fRangeParam->fSlave1Pad = 0;
90 fRangeParam->fSlave2Axis = 0;
91 fRangeParam->fSlave2Pad = 0;
93 //__________________________________________________________________
94 virtual ~dNdetaDrawer()
96 if (fRatios && fRatios->GetHists()) fRatios->GetHists()->Delete();
97 if (fResults && fResults->GetHists()) fResults->GetHists()->Delete();
99 if (fTrigString) { delete fTrigString; fTrigString = 0; }
100 if (fSNNString) { delete fSNNString; fSNNString = 0; }
101 if (fSysString) { delete fSysString; fSysString = 0; }
102 if (fVtxAxis) { delete fVtxAxis; fVtxAxis = 0; }
103 if (fCentAxis) { delete fCentAxis; fCentAxis = 0; }
104 if (fResults) { delete fResults; fResults = 0; }
105 if (fRatios) { delete fRatios; fRatios = 0; }
106 if (fOthers) { delete fOthers; fOthers = 0; }
107 if (fTriggers) { delete fTriggers; fTriggers = 0; }
111 //==================================================================
114 * @name Set parameters
117 * Show other (UA5, CMS, ...) data
119 * @param x Whether to show or not
121 void SetShowOthers(UShort_t x) { fShowOthers = x; }
122 //__________________________________________________________________
124 * Show ALICE published data
126 * @param x Wheter to show or not
128 void SetShowAlice(Bool_t x) { fShowAlice = x; }
129 //__________________________________________________________________
131 * Whether to show ratios or not. If there's nothing to compare to,
132 * the ratio panel will be implicitly disabled
134 * @param x Whether to show or not
136 void SetShowRatios(Bool_t x) { fShowRatios = x; }
137 //__________________________________________________________________
140 * Whether to show the left/right asymmetry
142 * @param x To show or not
144 void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; }
145 //__________________________________________________________________
147 * Whether to show rings
149 * @param x To show or not
151 void SetShowRings(Bool_t x) { fShowRings = x; }
152 //__________________________________________________________________
154 * Set the rebinning factor
156 * @param x Rebinning factor (must be a divisor in the number of bins)
158 void SetRebin(UShort_t x) { fRebin = x; }
159 //__________________________________________________________________
161 * Wheter to cut away the edges
163 * @param x Whether or not to cut away edges
165 void SetCutEdges(Bool_t x) { fCutEdges = x; }
166 //__________________________________________________________________
168 * Set the title of the plot
172 void SetTitle(TString x) { fTitle = x; }
173 //__________________________________________________________________
175 * Set the systematic error in the forward region
177 * @param e Systematic error in the forward region
179 void SetForwardSysError(Double_t e=0) { fFwdSysErr = e; }
180 //__________________________________________________________________
182 * Set the systematic error in the forward region
184 * @param e Systematic error in the forward region
186 void SetCentralSysError(Double_t e=0) { fCenSysErr = e; }
188 //==================================================================
191 * @name Override settings from input
194 * Override setting from file
196 * @param sNN Center of mass energy per nucleon pair (GeV)
198 void SetSNN(UShort_t sNN)
200 fSNNString = new TNamed("sNN", Form("%04dGeV", sNN));
201 fSNNString->SetUniqueID(sNN);
203 //__________________________________________________________________
205 * Set the collision system
209 * @param sys collision system
211 void SetSys(UShort_t sys)
213 fSysString = new TNamed("sys", (sys == 1 ? "pp" :
214 sys == 2 ? "PbPb" : "unknown"));
215 fSysString->SetUniqueID(sys);
217 //__________________________________________________________________
219 * Set the vertex range in centimeters
221 * @param vzMin Min @f$ v_z@f$
222 * @param vzMax Max @f$ v_z@f$
224 void SetVertexRange(Double_t vzMin, Double_t vzMax)
226 fVtxAxis = new TAxis(10, vzMin, vzMax);
227 fVtxAxis->SetName("vtxAxis");
228 fVtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", vzMin, vzMax));
230 //__________________________________________________________________
231 void SetTrigger(UShort_t trig)
233 fTrigString = new TNamed("trigString", (trig & 0x1 ? "INEL" :
234 trig & 0x2 ? "INEL>0" :
237 fTrigString->SetUniqueID(trig);
241 //==================================================================
244 * @name Main steering functions
247 * Run the code to produce the final result.
249 * @param filename File containing the data
251 void Run(const char* filename="forward_dndeta.root")
253 Double_t max = 0, rmax=0, amax=0;
255 gStyle->SetPalette(1);
257 // --- Open input file -------------------------------------------
258 TFile* file = TFile::Open(filename, "READ");
260 Error("Open", "Cannot open %s", filename);
263 // --- Get forward list ------------------------------------------
264 TList* forward = static_cast<TList*>(file->Get("ForwardResults"));
266 Error("Open", "Couldn't find list ForwardResults");
269 // --- Get information on the run --------------------------------
270 FetchInformation(forward);
271 // --- Set the macro pathand load other data script --------------
272 gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
273 gROOT->GetMacroPath()));
274 gROOT->LoadMacro("OtherData.C");
276 // --- Get the central results -----------------------------------
277 TList* clusters = static_cast<TList*>(file->Get("CentralResults"));
278 if (!clusters) Warning("Open", "Couldn't find list CentralResults");
280 // --- Get the central results -----------------------------------
281 TList* mcTruth = static_cast<TList*>(file->Get("MCTruthResults"));
282 if (!mcTruth) Warning("Open", "Couldn't find list MCTruthResults");
284 // --- Make our containtes ---------------------------------------
285 fResults = new THStack("results", "Results");
286 fRatios = new THStack("ratios", "Ratios");
287 fLeftRight = new THStack("asymmetry", "Left-right asymmetry");
288 fOthers = new TMultiGraph();
290 // --- Loop over input data --------------------------------------
291 /*TObjArray* mcA =*/ FetchResults(mcTruth, "MCTruth", max, rmax, amax);
292 TObjArray* fwdA = FetchResults(forward, "Forward", max, rmax, amax);
293 TObjArray* cenA = FetchResults(clusters, "Central", max, rmax, amax);
295 // --- Get trigger information -----------------------------------
296 TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
298 TList* all = static_cast<TList*>(sums->FindObject("all"));
300 fTriggers = FetchResult(all, "triggers");
301 if (!fTriggers) all->ls();
304 Warning("Run", "List all not found in ForwardSums");
309 Warning("Run", "No ForwardSums directory found in %s", file->GetName());
313 // --- Check our stacks ------------------------------------------
314 if (!fResults->GetHists() ||
315 fResults->GetHists()->GetEntries() <= 0) {
316 Error("Run", "No histograms in result stack!");
319 if (!fOthers->GetListOfGraphs() ||
320 fOthers->GetListOfGraphs()->GetEntries() <= 0) {
321 Warning("Run", "No other data found - disabling that");
324 if (!fRatios->GetHists() ||
325 fRatios->GetHists()->GetEntries() <= 0) {
326 Warning("Run", "No ratio data found - disabling that");
330 if (!fLeftRight->GetHists() ||
331 fLeftRight->GetHists()->GetEntries() <= 0) {
332 Warning("Run", "No left/right data found - disabling that");
334 fShowLeftRight = false;
336 if (fFwdSysErr > 0) {
337 if (fCenSysErr <= 0) fCenSysErr = fFwdSysErr;
338 for (Int_t i = 0; i < fwdA->GetEntriesFast(); i++) {
339 TH1* fwd = static_cast<TH1*>(fwdA->At(i));
340 TH1* cen = static_cast<TH1*>(cenA->At(i));
344 TH1* tmp = Merge(cen, fwd, low, high);
345 TF1* f = FitMerged(tmp, low, high);
346 MakeSysError(tmp, cen, fwd, f);
348 Info("", "Adding systematic error histogram %s",
350 fResults->GetHists()->AddFirst(tmp, "e5");
351 TH1* tmp2 = Symmetrice(tmp);
352 tmp2->SetFillColor(tmp->GetFillColor());
353 tmp2->SetFillStyle(tmp->GetFillStyle());
354 tmp2->SetMarkerStyle(tmp->GetMarkerStyle());
355 tmp2->SetLineWidth(tmp->GetLineWidth());
356 fResults->GetHists()->AddFirst(tmp2, "e5");
357 fResults->Modified();
363 // --- Close the input file --------------------------------------
368 // --- Plot results ----------------------------------------------
369 Plot(max, rmax, amax);
372 //__________________________________________________________________
374 * Fetch the information on the run from the results list
376 * @param results Results list
378 void FetchInformation(const TList* results)
381 fTrigString = static_cast<TNamed*>(results->FindObject("trigger"));
383 fNormString = static_cast<TNamed*>(results->FindObject("scheme"));
385 fSNNString = static_cast<TNamed*>(results->FindObject("sNN"));
387 fSysString = static_cast<TNamed*>(results->FindObject("sys"));
389 fVtxAxis = static_cast<TAxis*>(results->FindObject("vtxAxis"));
391 fCentAxis = static_cast<TAxis*>(results->FindObject("centAxis"));
393 TNamed* options = static_cast<TAxis*>(results->FindObject("options"));
394 if (!fTrigString) fTrigString = new TNamed("trigger", "unknown");
395 if (!fNormString) fNormString = new TNamed("scheme", "unknown");
396 if (!fSNNString) fSNNString = new TNamed("sNN", "unknown");
397 if (!fSysString) fSysString = new TNamed("sys", "unknown");
399 fVtxAxis = new TAxis(1,0,0);
400 fVtxAxis->SetName("vtxAxis");
401 fVtxAxis->SetTitle("v_{z} range unspecified");
404 TString centTxt("none");
406 Int_t nCent = fCentAxis->GetNbins();
407 centTxt = Form("%d bins", nCent);
408 for (Int_t i = 0; i <= nCent; i++)
409 centTxt.Append(Form("%c%d", i == 0 ? ' ' : '-',
410 int(fCentAxis->GetXbins()->At(i))));
412 Info("FetchInformation",
414 " Trigger: %-30s (%d)\n"
415 " sqrt(sNN): %-30s (%dGeV)\n"
416 " System: %-30s (%d)\n"
417 " Vz range: %-30s (%f,%f)\n"
418 " Normalization: %-30s (%d)\n"
421 fTrigString->GetTitle(), fTrigString->GetUniqueID(),
422 fSNNString->GetTitle(), fSNNString->GetUniqueID(),
423 fSysString->GetTitle(), fSysString->GetUniqueID(),
424 fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax(),
425 fNormString->GetTitle(), fNormString->GetUniqueID(),
426 centTxt.Data(), (options ? options->GetTitle() : "none"));
428 //__________________________________________________________________
429 TMultiGraph* FetchOthers(UShort_t centLow, UShort_t centHigh)
431 TMultiGraph* thisOther = 0;
432 if (fShowOthers == 0) return 0;
434 UShort_t sys = (fSysString ? fSysString->GetUniqueID() : 0);
435 UShort_t trg = (fTrigString ? fTrigString->GetUniqueID() : 0);
436 UShort_t snn = (fSNNString ? fSNNString->GetUniqueID() : 0);
437 Long_t ret = gROOT->ProcessLine(Form("GetData(%d,%d,%d,%d,%d,%d);",
443 thisOther = reinterpret_cast<TMultiGraph*>(ret);
446 //__________________________________________________________________
448 * Get the results from the top-level list
452 * @param max On return, maximum of data
453 * @param rmax On return, maximum of ratios
454 * @param amax On return, maximum of left-right comparisons
457 FetchResults(const TList* list,
463 UShort_t n = fCentAxis ? fCentAxis->GetNbins() : 0;
465 TList* all = static_cast<TList*>(list->FindObject("all"));
467 Error("FetchResults", "Couldn't find list 'all' in %s",
471 TObjArray* a = new TObjArray;
472 TH1* h = FetchResults(all, name, FetchOthers(0,0),
473 -1000, 0, max, rmax, amax);
478 TObjArray* a = new TObjArray;
479 for (UShort_t i = 0; i < n; i++) {
480 UShort_t centLow = fCentAxis->GetBinLowEdge(i+1);
481 UShort_t centHigh = fCentAxis->GetBinUpEdge(i+1);
482 TString lname = Form("cent%03d_%03d", centLow, centHigh);
483 TList* thisCent = static_cast<TList*>(list->FindObject(lname));
484 Int_t col = GetCentralityColor(i+1);
486 TString centTxt = Form("%3d%%-%3d%% central", centLow, centHigh);
488 Error("FetchResults", "Couldn't find list '%s' in %s",
489 lname.Data(), list->GetName());
492 TH1* h = FetchResults(thisCent, name, FetchOthers(centLow, centHigh),
493 col, centTxt.Data(), max, rmax, amax);
498 //__________________________________________________________________
499 Int_t GetCentralityColor(Int_t bin) const
501 UShort_t centLow = fCentAxis->GetBinLowEdge(bin);
502 UShort_t centHigh = fCentAxis->GetBinUpEdge(bin);
503 Float_t fc = (centLow+double(centHigh-centLow)/2) / 100;
504 Int_t nCol = gStyle->GetNumberOfColors();
505 Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5));
506 Int_t col = gStyle->GetColorPalette(icol);
507 //Info("GetCentralityColor","%3d: %3d-%3d -> %3d",bin,centLow,centHigh,col);
510 //__________________________________________________________________
511 void SetAttributes(TH1* h, Int_t color)
514 if (color < 0) return;
515 // h->SetLineColor(color);
516 h->SetMarkerColor(color);
517 // h->SetFillColor(color);
519 //__________________________________________________________________
520 void SetAttributes(TGraph* g, Int_t color)
523 if (color < 0) return;
524 // g->SetLineColor(color);
525 g->SetMarkerColor(color);
526 // g->SetFillColor(color);
528 //__________________________________________________________________
529 void ModifyTitle(TNamed* /*h*/, const char* /*centTxt*/)
532 // if (!centTxt || !h) return;
533 // h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
536 //__________________________________________________________________
538 * Fetch results for a particular centrality bin
542 * @param thisOther Other graphs
544 * @param centTxt Centrality text
545 * @param max On return, data maximum
546 * @param rmax On return, ratio maximum
547 * @param amax On return, left-right maximum
549 TH1* FetchResults(const TList* list,
551 TMultiGraph* thisOther,
558 TH1* dndeta = FetchResult(list, Form("dndeta%s", name));
559 TH1* dndetaMC = FetchResult(list, Form("dndeta%sMC", name));
560 TH1* dndetaTruth = FetchResult(list, "dndetaTruth");
562 TH1* dndetaMCSym = 0;
563 SetAttributes(dndeta, color);
564 SetAttributes(dndetaMC, fCentAxis ? color : color+2);
565 SetAttributes(dndetaTruth,color);
566 SetAttributes(dndetaSym, color);
567 SetAttributes(dndetaMCSym,fCentAxis ? color : color+2);
568 if (dndetaMC && fCentAxis)
569 dndetaMC->SetMarkerStyle(dndetaMC->GetMarkerStyle()+2);
570 if (dndetaMCSym && fCentAxis)
571 dndetaMCSym->SetMarkerStyle(dndetaMCSym->GetMarkerStyle()+2);
572 if (dndetaTruth && fCentAxis) {
573 dndetaTruth->SetMarkerStyle(34);
574 dndetaTruth->SetMarkerColor(kYellow-1);
576 ModifyTitle(dndeta, centTxt);
577 ModifyTitle(dndetaMC, centTxt);
578 ModifyTitle(dndetaTruth,centTxt);
579 ModifyTitle(dndetaSym, centTxt);
580 ModifyTitle(dndetaMCSym,centTxt);
583 max = TMath::Max(max, AddHistogram(fResults, dndetaTruth, "e5 p"));
584 max = TMath::Max(max, AddHistogram(fResults, dndetaMC, dndetaMCSym));
585 max = TMath::Max(max, AddHistogram(fResults, dndeta, dndetaSym));
588 fTruth = dndetaTruth;
592 THStack* rings = static_cast<THStack*>(list->FindObject("dndetaRings"));
594 TIter next(rings->GetHists());
596 while ((hist = static_cast<TH1*>(next())))
597 max = TMath::Max(max, AddHistogram(fResults, hist));
600 // Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f",
601 // dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
603 if (fShowLeftRight) {
604 fLeftRight->Add(Asymmetry(dndeta, amax));
605 fLeftRight->Add(Asymmetry(dndetaMC, amax));
609 TIter next(thisOther->GetListOfGraphs());
611 while ((g = static_cast<TGraph*>(next()))) {
612 fRatios->Add(Ratio(dndeta, g, rmax));
613 fRatios->Add(Ratio(dndetaSym, g, rmax));
614 SetAttributes(g, color);
615 ModifyTitle(g, centTxt);
616 if (!fOthers->GetListOfGraphs() ||
617 !fOthers->GetListOfGraphs()->FindObject(g->GetName())) {
618 max = TMath::Max(max,TMath::MaxElement(g->GetN(), g->GetY()));
622 // fOthers->Add(thisOther);
626 fRatios->Add(Ratio(dndeta, dndetaMC, rmax));
627 fRatios->Add(Ratio(dndetaSym, dndetaMCSym, rmax));
630 fRatios->Add(Ratio(dndeta, fTruth, rmax));
631 fRatios->Add(Ratio(dndetaSym, fTruth, rmax));
635 //__________________________________________________________________
638 * @param max Max value
639 * @param rmax Maximum diviation from 1 of ratios
640 * @param amax Maximum diviation from 1 of asymmetries
642 void Plot(Double_t max,
646 gStyle->SetOptTitle(0);
647 gStyle->SetTitleFont(132, "xyz");
648 gStyle->SetLabelFont(132, "xyz");
651 Int_t w = 800; // h / TMath::Sqrt(2);
655 if (!fShowRatios) w *= 1.3;
657 if (!fShowLeftRight) w *= 1.3;
660 y1 = (y11 > 0.0001 ? 0.4 : 0.2);
661 y2 = (y11 > 0.0001 ? 0.2 : 0.3);
663 TCanvas* c = new TCanvas("Results", "Results", w, h);
668 PlotResults(max, y1);
671 PlotRatios(rmax, y2, y1);
674 PlotLeftRight(amax, y3, y2);
678 Int_t vMin = fVtxAxis->GetXmin();
679 Int_t vMax = fVtxAxis->GetXmax();
680 TString trg(fTrigString->GetTitle());
682 if (fTriggers) nev = fTriggers->GetBinContent(1);
683 trg = trg.Strip(TString::kBoth);
684 trg.ReplaceAll(" ", "_");
685 trg.ReplaceAll(">", "Gt");
686 TString base(Form("dndeta_%s_%s_%s_%c%02d%c%02dcm_%09dev",
687 fSysString->GetTitle(),
688 fSNNString->GetTitle(),
690 vMin < 0 ? 'm' : 'p', TMath::Abs(vMin),
691 vMax < 0 ? 'm' : 'p', TMath::Abs(vMax),
693 c->SaveAs(Form("%s.png", base.Data()));
694 c->SaveAs(Form("%s.root", base.Data()));
695 c->SaveAs(Form("%s.C", base.Data()));
697 //__________________________________________________________________
701 * @param stack Stack to include
702 * @param mg (optional) Multi graph to include
703 * @param x1 Lower X coordinate in the range [0,1]
704 * @param y1 Lower Y coordinate in the range [0,1]
705 * @param x2 Upper X coordinate in the range [0,1]
706 * @param y2 Upper Y coordinate in the range [0,1]
708 void BuildLegend(THStack* stack, TMultiGraph* mg,
709 Double_t x1, Double_t y1, Double_t x2, Double_t y2)
711 TLegend* l = new TLegend(x1,y1,x2,y2);
712 l->SetNColumns(fCentAxis ? 1 : 2);
718 // Loop over items in stack and get unique items, while ignoring
719 // mirrored data and systematic error bands
720 TIter next(stack->GetHists());
724 Bool_t sysErrSeen = false;
725 while ((hist = static_cast<TH1*>(next()))) {
726 TString t(hist->GetTitle());
727 TString n(hist->GetName());
729 if (t.Contains("mirrored")) continue;
730 if (n.Contains("syserror")) { sysErrSeen = true; continue; }
731 if (unique.FindObject(t.Data())) continue;
732 TObjString* s = new TObjString(hist->GetTitle());
733 s->SetUniqueID(((hist->GetMarkerStyle() & 0xFFFF) << 16) |
734 ((hist->GetMarkerColor() & 0xFFFF) << 0));
736 // l->AddEntry(hist, hist->GetTitle(), "pl");
739 // If we have other stuff, scan for unique names
740 TIter nexto(mg->GetListOfGraphs());
742 while ((g = static_cast<TGraph*>(nexto()))) {
743 TString n(g->GetTitle());
744 if (n.Contains("mirrored")) continue;
745 if (unique.FindObject(n.Data())) continue;
746 TObjString* s = new TObjString(n);
747 s->SetUniqueID(((g->GetMarkerStyle() & 0xFFFF) << 16) |
748 ((g->GetMarkerColor() & 0xFFFF) << 0));
750 // l->AddEntry(hist, hist->GetTitle(), "pl");
754 // Add legend entries for unique items only
755 TIter nextu(&unique);
758 while ((s = nextu())) {
759 TLegendEntry* dd = l->AddEntry(Form("data%2d", i++),
761 Int_t style = (s->GetUniqueID() >> 16) & 0xFFFF;
762 Int_t color = (s->GetUniqueID() >> 0) & 0xFFFF;
763 dd->SetLineColor(kBlack);
764 if (fCentAxis) dd->SetMarkerColor(kBlack);
765 else dd->SetMarkerColor(color);
766 dd->SetMarkerStyle(style);
769 // Add entry for systematic errors
770 TLegendEntry* d0 = l->AddEntry("d0", Form("%4.1f%% Systematic error",
771 100*fFwdSysErr), "f");
772 d0->SetLineColor(kGray);
773 d0->SetMarkerColor(kGray);
774 d0->SetFillColor(kGray);
775 d0->SetFillStyle(3001);
776 d0->SetMarkerStyle(0);
780 if (!fCentAxis && i % 2 == 1) {
781 // To make sure the 'data' and 'mirrored' entries are on a line
783 TLegendEntry* dd = l->AddEntry("dd", " ", "");
789 dd->SetMarkerSize(0);
791 // Add entry for 'data'
792 TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp");
793 d1->SetLineColor(kBlack);
794 d1->SetMarkerColor(kBlack);
795 d1->SetMarkerStyle(20);
796 // Add entry for 'mirrored data'
797 TLegendEntry* d2 = l->AddEntry("d2", "Mirrored data", "lp");
798 d2->SetLineColor(kBlack);
799 d2->SetMarkerColor(kBlack);
800 d2->SetMarkerStyle(24);
804 //__________________________________________________________________
806 * Build centrality legend
808 * @param axis Stack to include
809 * @param x1 Lower X coordinate in the range [0,1]
810 * @param y1 Lower Y coordinate in the range [0,1]
811 * @param x2 Upper X coordinate in the range [0,1]
812 * @param y2 Upper Y coordinate in the range [0,1]
814 void BuildCentLegend(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
816 if (!fCentAxis) return;
818 TLegend* l = new TLegend(x1,y1,x2,y2);
825 Int_t n = fCentAxis->GetNbins();
826 for (Int_t i = 1; i <= n; i++) {
827 Double_t low = fCentAxis->GetBinLowEdge(i);
828 Double_t upp = fCentAxis->GetBinUpEdge(i);
829 TLegendEntry* e = l->AddEntry(Form("dummy%02d", i),
830 Form("%3d%% - %3d%%",
831 int(low), int(upp)), "pl");
832 e->SetMarkerColor(GetCentralityColor(i));
836 //__________________________________________________________________
841 * @param yd Bottom position of pad
843 void PlotResults(Double_t max, Double_t yd)
845 // Make a sub-pad for the result itself
846 TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0);
847 p1->SetTopMargin(0.05);
848 p1->SetBorderSize(0);
849 p1->SetBorderMode(0);
850 p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1);
851 p1->SetRightMargin(0.03);
852 if (fShowLeftRight || fShowRatios) p1->SetGridx();
858 // Info("PlotResults", "Plotting results with max=%f", max);
859 fResults->SetMaximum(1.15*max);
860 fResults->SetMinimum(yd > 0.00001 ? -0.1 : 0);
862 FixAxis(fResults, (1-yd)*(yd > .001 ? 1 : .9 / 1.2),
863 "#font[12]{#frac{1}{N} "
864 "#frac{dN_{#font[132]{ch}}}{d#font[152]{#eta}}}");
867 fResults->DrawClone("nostack e1");
869 fRangeParam->fSlave1Axis = fResults->GetXaxis();
870 fRangeParam->fSlave1Pad = p1;
873 if (fShowOthers != 0) {
874 TGraphAsymmErrors* o = 0;
875 TIter nextG(fOthers->GetListOfGraphs());
876 while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
877 o->DrawClone("same p");
880 // Make a legend in the main result pad
881 BuildCentLegend(.12, 1-p1->GetTopMargin()-.01-.5,
882 .35, 1-p1->GetTopMargin()-.01-.1);
883 Double_t x1 = (fCentAxis ? .7 : .15);
884 Double_t x2 = (fCentAxis ? 1-p1->GetRightMargin()-.01: .90);
885 Double_t y1 = (fCentAxis ? .5: p1->GetBottomMargin()+.01);
886 Double_t y2 = (fCentAxis ? 1-p1->GetTopMargin()-.01-.15 : .35);
888 BuildLegend(fResults, fOthers, x1, y1, x2, y2);
890 // Put a title on top
891 fTitle.ReplaceAll("@", " ");
892 TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
894 tit->SetTextFont(132);
895 tit->SetTextSize(0.05);
898 Int_t aliceBlue = TColor::GetColor(41,73,156);
899 // Put a nice label in the plot
901 UShort_t snn = fSNNString->GetUniqueID();
902 const char* sys = fSysString->GetTitle();
903 if (snn == 2750) snn = 2760;
904 if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
905 else eS = Form("%3dGeV", snn);
906 TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s%s}=%s, %s",
908 (fCentAxis ? "_{NN}" : ""),
910 fCentAxis ? "by centrality" :
911 fTrigString->GetTitle()));
912 tt->SetTextColor(aliceBlue);
914 tt->SetTextFont(132);
915 tt->SetTextAlign(33);
918 // Put number of accepted events on the plot
920 if (fTriggers) nev = fTriggers->GetBinContent(1);
921 TLatex* et = new TLatex(.93, .83, Form("%d events", nev));
922 et->SetTextColor(aliceBlue);
924 et->SetTextFont(132);
925 et->SetTextAlign(33);
928 // Put number of accepted events on the plot
930 TLatex* vt = new TLatex(.93, .88, fVtxAxis->GetTitle());
932 vt->SetTextFont(132);
933 vt->SetTextAlign(33);
934 vt->SetTextColor(aliceBlue);
937 // results->Draw("nostack e1 same");
939 fRangeParam->fSlave1Axis = FindXAxis(p1, fResults->GetName());
940 fRangeParam->fSlave1Pad = p1;
943 // Mark the plot as preliminary
944 TLatex* pt = new TLatex(.12, .93, "Work in progress");
947 // pt->SetTextSize();
948 pt->SetTextColor(TColor::GetColor(234,26,46));
949 pt->SetTextAlign(13);
952 const char* logos[] = { "ALICE.png", "FMD.png", 0 };
953 const char** logo = logos;
955 if (gSystem->AccessPathName(*logo)) {
959 TPad* pad = new TPad("logo", "logo", .12, .7, .25, .9, 0, 0, 0);
960 pad->SetFillStyle(0);
963 TImage* i = TImage::Create();
970 //__________________________________________________________________
974 * @param max Maximum diviation from 1
975 * @param y1 Lower y coordinate of pad
976 * @param y2 Upper y coordinate of pad
978 void PlotRatios(Double_t max, Double_t y1, Double_t y2)
980 if (fShowRatios == 0) return;
982 bool isBottom = (y1 < 0.0001);
983 Double_t yd = y2 - y1;
984 // Make a sub-pad for the result itself
985 TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0);
986 p2->SetTopMargin(0.001);
987 p2->SetRightMargin(0.03);
988 p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
996 FixAxis(fRatios, yd, "Ratios", 7);
998 fRatios->SetMaximum(1+TMath::Max(.22,1.05*max));
999 fRatios->SetMinimum(1-TMath::Max(.32,1.05*max));
1001 fRatios->DrawClone("nostack e1");
1005 BuildLegend(fRatios, 0, .15,p2->GetBottomMargin()+.01,.9,
1006 isBottom ? .6 : .4);
1008 TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,
1009 isBottom ? .6 : .4);
1011 l2->SetFillColor(0);
1012 l2->SetFillStyle(0);
1013 l2->SetBorderSize(0);
1014 l2->SetTextFont(132);
1016 // Make a nice band from 0.9 to 1.1
1017 TGraphErrors* band = new TGraphErrors(2);
1018 band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
1019 band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
1020 band->SetPointError(0, 0, .1);
1021 band->SetPointError(1, 0, .1);
1022 band->SetFillColor(kYellow+2);
1023 band->SetFillStyle(3002);
1024 band->SetLineStyle(2);
1025 band->SetLineWidth(1);
1026 band->Draw("3 same");
1027 band->DrawClone("X L same");
1029 // Replot the ratios on top
1030 fRatios->DrawClone("nostack e1 same");
1033 fRangeParam->fMasterAxis = FindXAxis(p2, fRatios->GetName());
1034 p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
1038 fRangeParam->fSlave2Axis = FindXAxis(p2, fRatios->GetName());
1039 fRangeParam->fSlave2Pad = p2;
1042 //__________________________________________________________________
1044 * Plot the asymmetries
1046 * @param max Maximum diviation from 1
1047 * @param y1 Lower y coordinate of pad
1048 * @param y2 Upper y coordinate of pad
1050 void PlotLeftRight(Double_t max, Double_t y1, Double_t y2)
1052 if (!fShowLeftRight) return;
1054 bool isBottom = (y1 < 0.0001);
1055 Double_t yd = y2 - y1;
1056 // Make a sub-pad for the result itself
1057 TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0);
1058 p3->SetTopMargin(0.001);
1059 p3->SetRightMargin(0.03);
1060 p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
1068 FixAxis(fLeftRight, yd, "Right/Left", 4);
1070 fLeftRight->SetMaximum(1+TMath::Max(.12,1.05*max));
1071 fLeftRight->SetMinimum(1-TMath::Max(.15,1.05*max));
1073 fLeftRight->DrawClone("nostack e1");
1077 Double_t xx1 = (fCentAxis ? .7 : .15);
1078 Double_t xx2 = (fCentAxis ? 1-p3->GetRightMargin()-.01 : .90);
1079 Double_t yy1 = p3->GetBottomMargin()+.01;
1080 Double_t yy2 = (fCentAxis ? 1-p3->GetTopMargin()-.01-.15 : .5);
1081 BuildLegend(fLeftRight, 0, xx1, yy1, xx2, yy2);
1082 // TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
1083 // l2->SetNColumns(2);
1084 // l2->SetFillColor(0);
1085 // l2->SetFillStyle(0);
1086 // l2->SetBorderSize(0);
1087 // l2->SetTextFont(132);
1089 // Make a nice band from 0.9 to 1.1
1090 TGraphErrors* band = new TGraphErrors(2);
1091 band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
1092 band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
1093 band->SetPointError(0, 0, .05);
1094 band->SetPointError(1, 0, .05);
1095 band->SetFillColor(kYellow+2);
1096 band->SetFillStyle(3002);
1097 band->SetLineStyle(2);
1098 band->SetLineWidth(1);
1099 band->Draw("3 same");
1100 band->DrawClone("X L same");
1102 fLeftRight->DrawClone("nostack e1 same");
1104 fRangeParam->fMasterAxis = FindXAxis(p3, fLeftRight->GetName());
1105 p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
1109 fRangeParam->fSlave2Axis = FindXAxis(p3, fLeftRight->GetName());
1110 fRangeParam->fSlave2Pad = p3;
1114 //==================================================================
1117 * @name Data utility functions
1120 * Get a result from the passed list
1122 * @param list List to search
1123 * @param name Object name to search for
1127 TH1* FetchResult(const TList* list, const char* name) const
1129 if (!list) return 0;
1131 TH1* ret = static_cast<TH1*>(list->FindObject(name));
1135 Warning("GetResult", "Histogram %s not found", name);
1141 //__________________________________________________________________
1143 * Add a histogram to the stack after possibly rebinning it
1145 * @param stack Stack to add to
1146 * @param hist histogram
1147 * @param option Draw options
1149 * @return Maximum of histogram
1151 Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option="") const
1153 // Check if we have input
1154 if (!hist) return 0;
1159 stack->Add(hist, option);
1160 return hist->GetMaximum();
1162 //__________________________________________________________________
1164 * Add a histogram to the stack after possibly rebinning it
1166 * @param stack Stack to add to
1167 * @param hist histogram
1168 * @param option Draw options
1169 * @param sym On return, the data symmetriced (added to stack)
1171 * @return Maximum of histogram
1173 Double_t AddHistogram(THStack* stack, TH1* hist, TH1*& sym,
1174 Option_t* option="") const
1176 // Check if we have input
1177 if (!hist) return 0;
1181 stack->Add(hist, option);
1183 // Now symmetrice the histogram
1184 sym = Symmetrice(hist);
1185 stack->Add(sym, option);
1187 return hist->GetMaximum();
1190 //__________________________________________________________________
1194 * @param h Histogram to rebin
1198 virtual void Rebin(TH1* h) const
1200 if (fRebin <= 1) return;
1202 Int_t nBins = h->GetNbinsX();
1203 if(nBins % fRebin != 0) {
1204 Warning("Rebin", "Rebin factor %d is not a devisor of current number "
1205 "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
1210 TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
1212 tmp->SetDirectory(0);
1214 // The new number of bins
1215 Int_t nBinsNew = nBins / fRebin;
1216 for(Int_t i = 1;i<= nBinsNew; i++) {
1217 Double_t content = 0;
1221 for(Int_t j = 1; j<=fRebin;j++) {
1222 Int_t bin = (i-1)*fRebin + j;
1223 Double_t c = h->GetBinContent(bin);
1225 if (c <= 0) continue;
1228 if (h->GetBinContent(bin+1)<=0 ||
1229 h->GetBinContent(bin-1)) {
1230 Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
1231 bin, c, h->GetName(),
1232 bin+1, h->GetBinContent(bin+1),
1233 bin-1, h->GetBinContent(bin-1));
1237 Double_t e = h->GetBinError(bin);
1238 Double_t w = 1 / (e*e); // 1/c/c
1245 if(content > 0 && nbins > 1 ) {
1246 tmp->SetBinContent(i, wsum / sumw);
1247 tmp->SetBinError(i,1./TMath::Sqrt(sumw));
1251 // Finally, rebin the histogram, and set new content
1254 for(Int_t i = 1; i<= nBinsNew; i++) {
1255 h->SetBinContent(i,tmp->GetBinContent(i));
1256 h->SetBinError(i, tmp->GetBinError(i));
1261 //__________________________________________________________________
1263 * Make an extension of @a h to make it symmetric about 0
1265 * @param h Histogram to symmertrice
1267 * @return Symmetric extension of @a h
1269 TH1* Symmetrice(const TH1* h) const
1271 Int_t nBins = h->GetNbinsX();
1272 TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
1273 s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
1275 s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
1276 s->SetMarkerColor(h->GetMarkerColor());
1277 s->SetMarkerSize(h->GetMarkerSize());
1278 s->SetMarkerStyle(h->GetMarkerStyle()+4);
1279 s->SetFillColor(h->GetFillColor());
1280 s->SetFillStyle(h->GetFillStyle());
1283 // Find the first and last bin with data
1284 Int_t first = nBins+1;
1286 for (Int_t i = 1; i <= nBins; i++) {
1287 if (h->GetBinContent(i) <= 0) continue;
1288 first = TMath::Min(first, i);
1289 last = TMath::Max(last, i);
1292 Double_t xfirst = h->GetBinCenter(first-1);
1293 Int_t f1 = h->GetXaxis()->FindBin(-xfirst);
1294 Int_t l2 = s->GetXaxis()->FindBin(xfirst);
1295 for (Int_t i = f1, j=l2; i <= last; i++,j--) {
1296 s->SetBinContent(j, h->GetBinContent(i));
1297 s->SetBinError(j, h->GetBinError(i));
1299 // Fill in overlap bin
1300 s->SetBinContent(l2+1, h->GetBinContent(first));
1301 s->SetBinError(l2+1, h->GetBinError(first));
1304 //__________________________________________________________________
1306 * Calculate the left-right asymmetry of input histogram
1308 * @param h Input histogram
1309 * @param max On return, the maximum distance from 1 of the histogram
1313 TH1* Asymmetry(TH1* h, Double_t& max)
1317 TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName())));
1318 // Int_t oBins = h->GetNbinsX();
1319 // Double_t high = h->GetXaxis()->GetXmax();
1320 // Double_t low = h->GetXaxis()->GetXmin();
1321 // Double_t dBin = (high - low) / oBins;
1322 // Int_t tBins = Int_t(2*high/dBin+.5);
1323 // ret->SetBins(tBins, -high, high);
1324 ret->SetDirectory(0);
1326 ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
1327 ret->SetYTitle("Right/Left");
1328 Int_t nBins = h->GetNbinsX();
1329 for (Int_t i = 1; i <= nBins; i++) {
1330 Double_t x = h->GetBinCenter(i);
1333 Double_t c1 = h->GetBinContent(i);
1334 Double_t e1 = h->GetBinError(i);
1335 if (c1 <= 0) continue;
1337 Int_t j = h->FindBin(-x);
1338 if (j <= 0 || j > nBins) continue;
1340 Double_t c2 = h->GetBinContent(j);
1341 Double_t e2 = h->GetBinError(j);
1343 Double_t c12 = c1*c1;
1344 Double_t e = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12));
1346 Int_t k = ret->FindBin(x);
1347 ret->SetBinContent(k, c2/c1);
1348 ret->SetBinError(k, e);
1350 max = TMath::Max(max, RatioMax(ret));
1354 //__________________________________________________________________
1356 * Transform a graph into a histogram
1362 TH1* Graph2Hist(const TGraphAsymmErrors* g) const
1364 Int_t nBins = g->GetN();
1365 TArrayF bins(nBins+1);
1367 for (Int_t i = 0; i < nBins; i++) {
1368 Double_t x = g->GetX()[i];
1369 Double_t exl = g->GetEXlow()[i];
1370 Double_t exh = g->GetEXhigh()[i];
1371 bins.fArray[i] = x-exl;
1372 bins.fArray[i+1] = x+exh;
1373 Double_t dxi = exh+exl;
1374 if (i == 0) dx = dxi;
1375 else if (dxi != dx) dx = 0;
1377 TString name(g->GetName());
1378 TString title(g->GetTitle());
1381 h = new TH1D(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
1384 h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
1386 h->SetMarkerStyle(g->GetMarkerStyle());
1387 h->SetMarkerColor(g->GetMarkerColor());
1388 h->SetMarkerSize(g->GetMarkerSize());
1393 //==================================================================
1396 * @name Ratio utility functions
1399 * Get the maximum diviation from 1 in the passed ratio
1401 * @param h Ratio histogram
1403 * @return Max diviation from 1
1405 Double_t RatioMax(TH1* h) const
1407 Int_t nBins = h->GetNbinsX();
1409 for (Int_t i = 1; i <= nBins; i++) {
1410 Double_t c = h->GetBinContent(i);
1411 if (c == 0) continue;
1412 Double_t e = h->GetBinError(i);
1413 Double_t d = TMath::Abs(1-c-e);
1414 ret = TMath::Max(d, ret);
1418 //__________________________________________________________________
1420 * Make the ratio of h1 to h2
1422 * @param o1 First object (numerator)
1423 * @param o2 Second object (denominator)
1424 * @param max Maximum diviation from 1
1428 TH1* Ratio(const TObject* o1, const TObject* o2, Double_t& max) const
1430 if (!o1 || !o2) return 0;
1433 const TAttMarker* m1 = 0;
1434 const TAttMarker* m2 = 0;
1435 const TH1* h1 = dynamic_cast<const TH1*>(o1);
1438 const TH1* h2 = dynamic_cast<const TH1*>(o2);
1444 const TGraph* g2 = dynamic_cast<const TGraph*>(o2);
1452 const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(o1);
1455 const TGraphAsymmErrors* g2 =
1456 dynamic_cast<const TGraphAsymmErrors*>(o2);
1459 r = RatioGG(g1, g2);
1464 Warning("Ratio", "Don't know how to divide a %s (%s) with a %s (%s)",
1465 o1->ClassName(), o1->GetName(), o2->ClassName(), o2->GetName());
1468 // Check that the histogram isn't empty
1469 if (r->GetEntries() <= 0) {
1474 r->SetMarkerStyle(m2->GetMarkerStyle());
1475 r->SetMarkerColor(m1->GetMarkerColor());
1476 r->SetMarkerSize(0.9*m1->GetMarkerSize());
1477 r->SetName(Form("%s_over_%s", o1->GetName(), o2->GetName()));
1478 r->SetTitle(Form("%s / %s", o1->GetTitle(), o2->GetTitle()));
1480 max = TMath::Max(RatioMax(r), max);
1485 //__________________________________________________________________
1487 * Compute the ratio of @a h to @a g. @a g is evaluated at the bin
1490 * @param h Numerator
1495 TH1* RatioHG(const TH1* h, const TGraph* g) const
1497 if (!h || !g) return 0;
1499 TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
1501 Double_t xlow = g->GetX()[0];
1502 Double_t xhigh = g->GetX()[g->GetN()-1];
1503 if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
1505 for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1506 Double_t c = h->GetBinContent(i);
1507 if (c <= 0) continue;
1509 Double_t x = h->GetBinCenter(i);
1510 if (x < xlow || x > xhigh) continue;
1512 Double_t f = g->Eval(x);
1513 if (f <= 0) continue;
1515 ret->SetBinContent(i, c / f);
1516 ret->SetBinError(i, h->GetBinError(i) / f);
1520 //__________________________________________________________________
1522 * Make the ratio of h1 to h2
1524 * @param h1 First histogram (numerator)
1525 * @param h2 Second histogram (denominator)
1529 TH1* RatioHH(const TH1* h1, const TH1* h2) const
1531 if (!h1 || !h2) return 0;
1532 TH1* t1 = static_cast<TH1*>(h1->Clone("tmp"));
1536 //__________________________________________________________________
1538 * Calculate the ratio of two graphs - g1 / g2
1540 * @param g1 Numerator
1541 * @param g2 Denominator
1543 * @return g1 / g2 in a histogram
1545 TH1* RatioGG(const TGraphAsymmErrors* g1,
1546 const TGraphAsymmErrors* g2) const
1548 Int_t nBins = g1->GetN();
1549 TArrayF bins(nBins+1);
1551 for (Int_t i = 0; i < nBins; i++) {
1552 Double_t x = g1->GetX()[i];
1553 Double_t exl = g1->GetEXlow()[i];
1554 Double_t exh = g1->GetEXhigh()[i];
1555 bins.fArray[i] = x-exl;
1556 bins.fArray[i+1] = x+exh;
1557 Double_t dxi = exh+exl;
1558 if (i == 0) dx = dxi;
1559 else if (dxi != dx) dx = 0;
1563 h = new TH1F("tmp", "tmp", nBins, bins[0], bins[nBins]);
1566 h = new TH1F("tmp", "tmp", nBins, bins.fArray);
1569 Double_t low = g2->GetX()[0];
1570 Double_t high = g2->GetX()[g2->GetN()-1];
1571 if (low > high) { Double_t t = low; low = high; high = t; }
1572 for (Int_t i = 0; i < nBins; i++) {
1573 Double_t x = g1->GetX()[i];
1574 if (x < low-dx || x > high+dx) continue;
1575 Double_t c1 = g1->GetY()[i];
1576 Double_t e1 = g1->GetErrorY(i);
1577 Double_t c2 = g2->Eval(x);
1579 h->SetBinContent(i+1, c1 / c2);
1580 h->SetBinError(i+1, e1 / c2);
1585 //==================================================================
1588 * @name Graphics utility functions
1591 * Find an X axis in a pad
1594 * @param name Histogram to find axis for
1596 * @return Found axis or null
1598 TAxis* FindXAxis(TVirtualPad* p, const char* name)
1600 TObject* o = p->GetListOfPrimitives()->FindObject(name);
1602 Warning("FindXAxis", "%s not found in pad", name);
1605 THStack* stack = dynamic_cast<THStack*>(o);
1607 Warning("FindXAxis", "%s is not a THStack", name);
1610 if (!stack->GetHistogram()) {
1611 Warning("FindXAxis", "%s has no histogram", name);
1614 TAxis* ret = stack->GetHistogram()->GetXaxis();
1618 //__________________________________________________________________
1620 * Fix the apperance of the axis in a stack
1622 * @param stack stack of histogram
1623 * @param s Scaling factor
1624 * @param ytitle Y axis title
1625 * @param force Whether to draw the stack first or not
1626 * @param ynDiv Divisions on Y axis
1628 void FixAxis(THStack* stack, Double_t yd, const char* ytitle,
1629 Int_t ynDiv=210, Bool_t force=true)
1632 Warning("FixAxis", "No stack passed for %s!", ytitle);
1635 if (force) stack->Draw("nostack e1");
1637 TH1* h = stack->GetHistogram();
1639 Warning("FixAxis", "Stack %s has no histogram", stack->GetName());
1642 Double_t s = 1/yd/1.2;
1643 // Info("FixAxis", "for %s, s=1/%f=%f", stack->GetName(), yd, s);
1645 h->SetXTitle("#font[152]{#eta}");
1646 h->SetYTitle(ytitle);
1647 TAxis* xa = h->GetXaxis();
1648 TAxis* ya = h->GetYaxis();
1650 // Int_t npixels = 20
1651 // Float_t dy = gPad->AbsPixeltoY(0) - gPad->AbsPixeltoY(npixels);
1652 // Float_t ts = dy/(gPad->GetY2() - gPad->GetY1());
1655 // xa->SetTitle(h->GetXTitle());
1656 // xa->SetTicks("+-");
1657 xa->SetTitleSize(s*xa->GetTitleSize());
1658 xa->SetLabelSize(s*xa->GetLabelSize());
1659 xa->SetTickLength(s*xa->GetTickLength());
1660 // xa->SetTitleOffset(xa->GetTitleOffset()/s);
1662 if (stack != fResults) {
1663 TAxis* rxa = fResults->GetXaxis();
1664 xa->Set(rxa->GetNbins(), rxa->GetXmin(), rxa->GetXmax());
1668 // ya->SetTitle(h->GetYTitle());
1670 // ya->SetTicks("+-");
1671 ya->SetNdivisions(ynDiv);
1672 ya->SetTitleSize(s*ya->GetTitleSize());
1673 ya->SetTitleOffset(ya->GetTitleOffset()/s);
1674 ya->SetLabelSize(s*ya->GetLabelSize());
1677 //__________________________________________________________________
1679 * Merge two histograms into one
1681 * @param cen Central part
1682 * @param fwd Forward part
1683 * @param xlow On return, lower eta bound
1684 * @param xhigh On return, upper eta bound
1686 * @return Newly allocated histogram or null
1689 Merge(const TH1* cen, const TH1* fwd, Double_t& xlow, Double_t& xhigh)
1691 TH1* tmp = static_cast<TH1*>(fwd->Clone("tmp"));
1692 TString name(fwd->GetName());
1693 name.ReplaceAll("Forward", "Merged");
1696 // tmp->SetMarkerStyle(28);
1697 // tmp->SetMarkerColor(kBlack);
1698 tmp->SetDirectory(0);
1701 for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
1702 Double_t cc = cen->GetBinContent(i);
1703 Double_t cf = fwd->GetBinContent(i);
1704 Double_t ec = cen->GetBinError(i);
1705 Double_t ef = fwd->GetBinError(i);
1708 if (cc < 0.001 && cf < 0.01) continue;
1709 xlow = TMath::Min(tmp->GetXaxis()->GetBinLowEdge(i),xlow);
1710 xhigh = TMath::Max(tmp->GetXaxis()->GetBinUpEdge(i),xhigh);
1716 ne = TMath::Sqrt(ec*ec + ef*ef);
1719 tmp->SetBinContent(i, nc);
1720 tmp->SetBinError(i, ne);
1724 //____________________________________________________________________
1726 * Fit @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$ to histogram data
1728 * @param tmp Histogram
1729 * @param xlow Lower x bound
1730 * @param xhigh Upper x bound
1732 * @return Fitted function
1735 FitMerged(TH1* tmp, Double_t xlow, Double_t xhigh)
1737 TF1* tmpf = new TF1("tmpf", "gaus", xlow, xhigh);
1738 tmp->Fit(tmpf, "NQ", "");
1739 tmp->SetDirectory(0);
1741 TF1* fit = new TF1("f", myFunc, xlow, xhigh, 4);
1742 fit->SetParNames("a_{1}", "a_{2}", "#sigma_{1}", "#sigma_{2}");
1743 fit->SetParameters(tmpf->GetParameter(0),
1745 tmpf->GetParameter(2),
1746 tmpf->GetParameter(2)/4);
1747 fit->SetParLimits(3, 0, 100);
1748 fit->SetParLimits(4, 0, 100);
1749 tmp->Fit(fit,"0W","");
1754 //____________________________________________________________________
1756 * Make band of systematic errors
1758 * @param tmp Histogram
1762 MakeSysError(TH1* tmp, TH1* cen, TH1* fwd, TF1* fit)
1764 for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
1765 Double_t tc = tmp->GetBinContent(i);
1766 if (tc < 0.01) continue;
1767 Double_t fc = fwd->GetBinContent(i);
1768 Double_t cc = cen->GetBinContent(i);
1769 Double_t sysErr = fFwdSysErr;
1770 if (cc > .01 && fc > 0.01)
1771 sysErr = (fFwdSysErr+fCenSysErr) / 2;
1773 sysErr = fCenSysErr;
1774 Double_t x = tmp->GetXaxis()->GetBinCenter(i);
1775 Double_t y = fit->Eval(x);
1776 tmp->SetBinContent(i, y);
1777 tmp->SetBinError(i,sysErr*y);
1779 TString name(tmp->GetName());
1780 name.ReplaceAll("Merged", "SysError");
1782 tmp->SetFillColor(kGray);
1783 tmp->SetFillStyle(3001);
1784 tmp->SetMarkerStyle(0);
1785 tmp->SetLineWidth(0);
1787 void CorrectForward(TH1* h) const
1789 if (!fRemoveOuters) return;
1791 for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1792 Double_t eta = h->GetBinCenter(i);
1793 if (TMath::Abs(eta) < 2.3) {
1794 h->SetBinContent(i, 0);
1795 h->SetBinError(i, 0);
1799 void CorrectCentral(TH1* h) const
1801 if (fClusterScale.IsNull()) return;
1802 TString t(h->GetTitle());
1803 Info("CorrectCentral", "Replacing Central with Tracklets in %s", t.Data());
1804 t.ReplaceAll("Central", "Tracklets");
1807 TF1* cf = new TF1("clusterScale", fClusterScale);
1808 for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1809 Double_t eta = h->GetBinCenter(i);
1810 Double_t f = cf->Eval(eta);
1811 Double_t c = h->GetBinContent(i);
1813 h->SetBinContent(i, c / f);
1823 //__________________________________________________________________
1824 UShort_t fShowOthers; // Show other data
1825 Bool_t fShowAlice; // Show ALICE published data
1826 Bool_t fShowRatios; // Show ratios
1827 Bool_t fShowLeftRight;// Show asymmetry
1828 Bool_t fShowRings; // Show rings too
1829 UShort_t fRebin; // Rebinning factor
1830 Bool_t fCutEdges; // Whether to cut edges
1831 TString fTitle; // Title on plot
1832 TNamed* fTrigString; // Trigger string (read, or set)
1833 TNamed* fNormString; // Normalisation string (read, or set)
1834 TNamed* fSNNString; // Energy string (read, or set)
1835 TNamed* fSysString; // Collision system string (read or set)
1836 TAxis* fVtxAxis; // Vertex cuts (read or set)
1837 TAxis* fCentAxis; // Centrality axis
1838 THStack* fResults; // Stack of results
1839 THStack* fRatios; // Stack of ratios
1840 THStack* fLeftRight; // Left-right asymmetry
1841 TMultiGraph* fOthers; // Older data
1842 TH1* fTriggers; // Number of triggers
1843 TH1* fTruth; // Pointer to truth
1844 RangeParam* fRangeParam; // Parameter object for range zoom
1845 Double_t fFwdSysErr; // Systematic error in forward range
1846 Double_t fCenSysErr; // Systematic error in central range
1847 Bool_t fRemoveOuters; // Whether to remove outers
1848 TString fClusterScale; // Scaling of clusters to tracklets
1851 //____________________________________________________________________
1853 * Function to calculate
1855 * g(x;A_1,A_2,\sigma_1,\sigma_2) =
1856 * A_1\left(\frac{1}{2\pi\sigma_1}e^{(x/\sigma_1)^2} -
1857 * A_2\frac{1}{2\pi\sigma_2}e^{(x/\sigma_2)^2}\right)
1860 * @param xp Pointer to x array
1861 * @param pp Pointer to parameter array
1863 * @return @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$
1865 Double_t myFunc(Double_t* xp, Double_t* pp)
1868 Double_t a1 = pp[0];
1869 Double_t a2 = pp[1];
1870 Double_t s1 = pp[2];
1871 Double_t s2 = pp[3];
1872 return a1*(TMath::Gaus(x, 0, s1) - a2 * TMath::Gaus(x, 0, s2));
1875 //=== Stuff for auto zooming =========================================
1876 void UpdateRange(dNdetaDrawer::RangeParam* p)
1879 Warning("UpdateRange", "No parameters %p", p);
1882 if (!p->fMasterAxis) {
1883 Warning("UpdateRange", "No master axis %p", p->fMasterAxis);
1886 Int_t first = p->fMasterAxis->GetFirst();
1887 Int_t last = p->fMasterAxis->GetLast();
1888 Double_t x1 = p->fMasterAxis->GetBinCenter(first);
1889 Double_t x2 = p->fMasterAxis->GetBinCenter(last);
1891 if (p->fSlave1Axis) {
1892 Int_t i1 = p->fSlave1Axis->FindBin(x1);
1893 Int_t i2 = p->fSlave1Axis->FindBin(x2);
1894 p->fSlave1Axis->SetRange(i1, i2);
1895 p->fSlave1Pad->Modified();
1896 p->fSlave1Pad->Update();
1898 if (p->fSlave2Axis) {
1899 Int_t i1 = p->fSlave2Axis->FindBin(x1);
1900 Int_t i2 = p->fSlave2Axis->FindBin(x2);
1901 p->fSlave2Axis->SetRange(i1, i2);
1902 p->fSlave2Pad->Modified();
1903 p->fSlave2Pad->Update();
1905 TCanvas* c = gPad->GetCanvas();
1909 //____________________________________________________________________
1910 void RangeExec(dNdetaDrawer::RangeParam* p)
1917 // 11: Mouse release
1918 // dNdetaDrawer::RangeParam* p =
1919 // reinterpret_cast<dNdetaDrawer::RangeParam*>(addr);
1920 Int_t event = gPad->GetEvent();
1921 TObject *select = gPad->GetSelected();
1926 if (event != 11 || !select || select != p->fMasterAxis) return;
1930 //=== Steering function ==============================================
1932 * Draw @f$ dN/d\eta@f$
1934 * @param filename File name
1935 * @param flags Flags
1936 * @param title Title
1937 * @param rebin Rebinning factor
1938 * @param cutEdges Whether to cut edges when rebinning
1939 * @param sNN (optional) Collision energy [GeV]
1940 * @param sys (optional) Collision system (1: pp, 2: PbPb)
1941 * @param trg (optional) Trigger (1: INEL, 2: INEL>0, 4: NSD)
1942 * @param vzMin Least @f$ v_z@f$
1943 * @param vzMax Largest @f$ v_z@f$
1945 * @ingroup pwg2_forward_dndeta
1948 DrawdNdeta(const char* filename="forward_dndeta.root",
1949 const char* title="",
1951 UShort_t others=0x7,
1959 dNdetaDrawer* pd = new dNdetaDrawer;
1960 dNdetaDrawer& d = *pd;
1963 d.SetShowOthers(others);
1964 d.SetShowRatios(flags & 0x1);
1965 d.SetShowLeftRight(flags & 0x2);
1966 d.SetForwardSysError(flags & 0x4 ? 0.076 : 0);
1967 d.SetShowRings(flags & 0x8);
1968 d.SetCutEdges(flags & 0x10);
1969 // d.fRemoveOuters = (flags & 0x20);
1970 // d.fClusterScale = "1.06 -0.003*x +0.0119*x*x";
1971 // Do the below if your input data does not contain these settings
1972 if (sNN > 0) d.SetSNN(sNN); // Collision energy per nucleon pair (GeV)
1973 if (sys > 0) d.SetSys(sys); // Collision system (1:pp, 2:PbPB)
1974 if (trg > 0) d.SetTrigger(trg); // Collision trigger (1:INEL, 2:INEL>0, 4:NSD)
1975 if (vzMin < 999 && vzMax > -999)
1976 d.SetVertexRange(vzMin,vzMax); // Collision vertex range (cm)
1979 //____________________________________________________________________