Mior fixes
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / DrawdNdeta.C
CommitLineData
e1f47419 1/**
ffca499d 2 * @file DrawdNdeta.C
3 * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
4 * @date Wed Mar 23 14:07:10 2011
5 *
5bb5d1f6 6 * @brief Script to visualise the dN/deta for pp and PbPb
e1f47419 7 *
ffca499d 8 * This script is independent of any AliROOT code - and should stay
9 * that way.
10 *
5bb5d1f6 11 * The script is <emph>very</emph> long - sigh - the joy of drawing
12 * things nicely in ROOT
ffca499d 13 *
14 * @ingroup pwg2_forward_dndeta
e1f47419 15 */
b2e7f2d6 16#include <TH1.h>
f709afcf 17#include <TColor.h>
b2e7f2d6 18#include <THStack.h>
19#include <TGraphErrors.h>
20#include <TGraphAsymmErrors.h>
21#include <TMultiGraph.h>
22#include <TFile.h>
23#include <TList.h>
24#include <TString.h>
25#include <TError.h>
26#include <TSystem.h>
27#include <TROOT.h>
28#include <TMath.h>
29#include <TCanvas.h>
30#include <TPad.h>
31#include <TStyle.h>
32#include <TLegend.h>
33#include <TLegendEntry.h>
34#include <TLatex.h>
35#include <TImage.h>
33ab4a11 36#include <TRandom.h>
37#include <fstream>
37079f20 38#define SYSERR_COLOR kBlue-10
39#define SYSERR_STYLE 1001
b2e7f2d6 40
797161e8 41Double_t myFunc(Double_t* xp, Double_t* pp);
42
2f92ff0e 43/**
44 * Class to draw dN/deta results
45 *
ffca499d 46 * @ingroup pwg2_forward_tasks_dndeta
47 * @ingroup pwg2_forward_dndeta
2f92ff0e 48 */
b2e7f2d6 49struct dNdetaDrawer
50{
2f92ff0e 51 /**
52 * POD of data for range zooming
53 */
b2e7f2d6 54 struct RangeParam
55 {
56 TAxis* fMasterAxis; // Master axis
57 TAxis* fSlave1Axis; // First slave axis
58 TVirtualPad* fSlave1Pad; // First slave pad
59 TAxis* fSlave2Axis; // Second slave axis
60 TVirtualPad* fSlave2Pad; // Second slave pad
61 };
62 //__________________________________________________________________
2f92ff0e 63 /**
64 * Constructor
65 *
66 */
b2e7f2d6 67 dNdetaDrawer()
33ab4a11 68 : // Options
69 fShowRatios(false), // Show ratios
70 fShowLeftRight(false), // Show asymmetry
71 fShowRings(false), // Show rings too
72 fExport(false), // Export data to script
73 fCutEdges(false), // Whether to cut edges
74 fRemoveOuters(false), // Whether to remove outers
75 fShowOthers(0), // Show other data
76 // Settings
77 fRebin(0), // Rebinning factor
78 fFwdSysErr(0.076), // Systematic error in forward range
79 fCenSysErr(0), // Systematic error in central range
80 fTitle(""), // Title on plot
81 fClusterScale(""), // Scaling of clusters to tracklets
82 // Read (or set) information
83 fTrigString(0), // Trigger string (read, or set)
84 fNormString(0), // Normalisation string (read, or set)
85 fSNNString(0), // Energy string (read, or set)
86 fSysString(0), // Collision system string (read or set)
87 fVtxAxis(0), // Vertex cuts (read or set)
88 fCentAxis(0), // Centrality axis
89 // Resulting plots
90 fResults(0), // Stack of results
91 fRatios(0), // Stack of ratios
92 fLeftRight(0), // Left-right asymmetry
93 fOthers(0), // Older data
94 fTriggers(0), // Number of triggers
95 fTruth(0), // Pointer to truth
96 fRangeParam(0) // Parameter object for range zoom
b2e7f2d6 97 {
98 fRangeParam = new RangeParam;
99 fRangeParam->fMasterAxis = 0;
100 fRangeParam->fSlave1Axis = 0;
101 fRangeParam->fSlave1Pad = 0;
102 fRangeParam->fSlave2Axis = 0;
103 fRangeParam->fSlave2Pad = 0;
104 }
33ab4a11 105 dNdetaDrawer(const dNdetaDrawer&) {}
106 dNdetaDrawer& operator=(const dNdetaDrawer&) { return *this; }
107
ffca499d 108 //__________________________________________________________________
109 virtual ~dNdetaDrawer()
110 {
111 if (fRatios && fRatios->GetHists()) fRatios->GetHists()->Delete();
112 if (fResults && fResults->GetHists()) fResults->GetHists()->Delete();
113
114 if (fTrigString) { delete fTrigString; fTrigString = 0; }
115 if (fSNNString) { delete fSNNString; fSNNString = 0; }
116 if (fSysString) { delete fSysString; fSysString = 0; }
117 if (fVtxAxis) { delete fVtxAxis; fVtxAxis = 0; }
118 if (fCentAxis) { delete fCentAxis; fCentAxis = 0; }
119 if (fResults) { delete fResults; fResults = 0; }
120 if (fRatios) { delete fRatios; fRatios = 0; }
121 if (fOthers) { delete fOthers; fOthers = 0; }
122 if (fTriggers) { delete fTriggers; fTriggers = 0; }
123 fRangeParam = 0;
124 }
125
2f92ff0e 126 //==================================================================
127 /**
128 * @{
129 * @name Set parameters
130 */
131 /**
132 * Show other (UA5, CMS, ...) data
133 *
134 * @param x Whether to show or not
135 */
797161e8 136 void SetShowOthers(UShort_t x) { fShowOthers = x; }
2f92ff0e 137 //__________________________________________________________________
2f92ff0e 138 /**
139 * Whether to show ratios or not. If there's nothing to compare to,
140 * the ratio panel will be implicitly disabled
141 *
142 * @param x Whether to show or not
143 */
b2e7f2d6 144 void SetShowRatios(Bool_t x) { fShowRatios = x; }
2f92ff0e 145 //__________________________________________________________________
146 /**
147 *
148 * Whether to show the left/right asymmetry
149 *
150 * @param x To show or not
151 */
b2e7f2d6 152 void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; }
797161e8 153 //__________________________________________________________________
154 /**
155 * Whether to show rings
156 *
157 * @param x To show or not
158 */
f709afcf 159 void SetShowRings(Bool_t x) { fShowRings = x; }
2f92ff0e 160 //__________________________________________________________________
33ab4a11 161 /**
162 * Whether to export results to a script
163 *
164 * @param x Wheter to export results to a script
165 */
166 void SetExport(Bool_t x) { fExport = x; }
167 //__________________________________________________________________
2f92ff0e 168 /**
169 * Set the rebinning factor
170 *
171 * @param x Rebinning factor (must be a divisor in the number of bins)
172 */
b2e7f2d6 173 void SetRebin(UShort_t x) { fRebin = x; }
2f92ff0e 174 //__________________________________________________________________
175 /**
176 * Wheter to cut away the edges
177 *
178 * @param x Whether or not to cut away edges
179 */
b2e7f2d6 180 void SetCutEdges(Bool_t x) { fCutEdges = x; }
2f92ff0e 181 //__________________________________________________________________
182 /**
183 * Set the title of the plot
184 *
185 * @param x Title
186 */
b2e7f2d6 187 void SetTitle(TString x) { fTitle = x; }
797161e8 188 //__________________________________________________________________
189 /**
190 * Set the systematic error in the forward region
191 *
192 * @param e Systematic error in the forward region
193 */
194 void SetForwardSysError(Double_t e=0) { fFwdSysErr = e; }
195 //__________________________________________________________________
196 /**
197 * Set the systematic error in the forward region
198 *
199 * @param e Systematic error in the forward region
200 */
201 void SetCentralSysError(Double_t e=0) { fCenSysErr = e; }
2f92ff0e 202 /* @} */
203 //==================================================================
204 /**
205 * @{
206 * @name Override settings from input
207 */
208 /**
209 * Override setting from file
210 *
211 * @param sNN Center of mass energy per nucleon pair (GeV)
212 */
213 void SetSNN(UShort_t sNN)
214 {
cdb9e20f 215 fSNNString = new TNamed("sNN", Form("%04dGeV", sNN));
2f92ff0e 216 fSNNString->SetUniqueID(sNN);
217 }
218 //__________________________________________________________________
219 /**
220 * Set the collision system
221 * - 1: pp
222 * - 2: PbPb
223 *
224 * @param sys collision system
225 */
226 void SetSys(UShort_t sys)
227 {
fe52e455 228 fSysString = new TNamed("sys", (sys == 1 ? "pp" :
2f92ff0e 229 sys == 2 ? "PbPb" : "unknown"));
fe52e455 230 fSysString->SetUniqueID(sys);
2f92ff0e 231 }
232 //__________________________________________________________________
233 /**
234 * Set the vertex range in centimeters
235 *
236 * @param vzMin Min @f$ v_z@f$
237 * @param vzMax Max @f$ v_z@f$
238 */
239 void SetVertexRange(Double_t vzMin, Double_t vzMax)
240 {
241 fVtxAxis = new TAxis(10, vzMin, vzMax);
242 fVtxAxis->SetName("vtxAxis");
243 fVtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", vzMin, vzMax));
244 }
b2e7f2d6 245 //__________________________________________________________________
2f92ff0e 246 void SetTrigger(UShort_t trig)
247 {
248 fTrigString = new TNamed("trigString", (trig & 0x1 ? "INEL" :
249 trig & 0x2 ? "INEL>0" :
250 trig & 0x4 ? "NSD" :
251 "unknown"));
252 fTrigString->SetUniqueID(trig);
253 }
254
255
256 //==================================================================
257 /**
258 * @{
259 * @name Main steering functions
260 */
261 /**
262 * Run the code to produce the final result.
263 *
264 * @param filename File containing the data
265 */
e308a636 266 void Run(const char* filename="forward_dndeta.root")
b2e7f2d6 267 {
e308a636 268 Double_t max = 0, rmax=0, amax=0;
b2e7f2d6 269
e308a636 270 gStyle->SetPalette(1);
e1f47419 271
e308a636 272 // --- Open input file -------------------------------------------
273 TFile* file = TFile::Open(filename, "READ");
274 if (!file) {
275 Error("Open", "Cannot open %s", filename);
276 return;
277 }
278 // --- Get forward list ------------------------------------------
279 TList* forward = static_cast<TList*>(file->Get("ForwardResults"));
280 if (!forward) {
281 Error("Open", "Couldn't find list ForwardResults");
282 return;
283 }
284 // --- Get information on the run --------------------------------
285 FetchInformation(forward);
286 // --- Set the macro pathand load other data script --------------
287 gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
288 gROOT->GetMacroPath()));
289 gROOT->LoadMacro("OtherData.C");
290
291 // --- Get the central results -----------------------------------
292 TList* clusters = static_cast<TList*>(file->Get("CentralResults"));
293 if (!clusters) Warning("Open", "Couldn't find list CentralResults");
b2e7f2d6 294
9ecab72f 295 // --- Get the central results -----------------------------------
296 TList* mcTruth = static_cast<TList*>(file->Get("MCTruthResults"));
297 if (!mcTruth) Warning("Open", "Couldn't find list MCTruthResults");
298
e308a636 299 // --- Make our containtes ---------------------------------------
300 fResults = new THStack("results", "Results");
301 fRatios = new THStack("ratios", "Ratios");
302 fLeftRight = new THStack("asymmetry", "Left-right asymmetry");
303 fOthers = new TMultiGraph();
b2e7f2d6 304
e308a636 305 // --- Loop over input data --------------------------------------
797161e8 306 /*TObjArray* mcA =*/ FetchResults(mcTruth, "MCTruth", max, rmax, amax);
307 TObjArray* fwdA = FetchResults(forward, "Forward", max, rmax, amax);
308 TObjArray* cenA = FetchResults(clusters, "Central", max, rmax, amax);
b2e7f2d6 309
e308a636 310 // --- Get trigger information -----------------------------------
311 TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
312 if (sums) {
313 TList* all = static_cast<TList*>(sums->FindObject("all"));
314 if (all) {
315 fTriggers = FetchResult(all, "triggers");
316 if (!fTriggers) all->ls();
317 }
318 else {
319 Warning("Run", "List all not found in ForwardSums");
320 sums->ls();
321 }
322 }
323 else {
324 Warning("Run", "No ForwardSums directory found in %s", file->GetName());
325 file->ls();
326 }
327
328 // --- Check our stacks ------------------------------------------
329 if (!fResults->GetHists() ||
330 fResults->GetHists()->GetEntries() <= 0) {
331 Error("Run", "No histograms in result stack!");
332 return;
333 }
334 if (!fOthers->GetListOfGraphs() ||
335 fOthers->GetListOfGraphs()->GetEntries() <= 0) {
336 Warning("Run", "No other data found - disabling that");
797161e8 337 fShowOthers = 0;
e308a636 338 }
339 if (!fRatios->GetHists() ||
340 fRatios->GetHists()->GetEntries() <= 0) {
341 Warning("Run", "No ratio data found - disabling that");
342 // fRatios->ls();
343 fShowRatios = false;
344 }
345 if (!fLeftRight->GetHists() ||
346 fLeftRight->GetHists()->GetEntries() <= 0) {
347 Warning("Run", "No left/right data found - disabling that");
348 // fLeftRight->ls();
349 fShowLeftRight = false;
350 }
797161e8 351 if (fFwdSysErr > 0) {
352 if (fCenSysErr <= 0) fCenSysErr = fFwdSysErr;
353 for (Int_t i = 0; i < fwdA->GetEntriesFast(); i++) {
354 TH1* fwd = static_cast<TH1*>(fwdA->At(i));
355 TH1* cen = static_cast<TH1*>(cenA->At(i));
356 CorrectForward(fwd);
357 CorrectCentral(cen);
358 Double_t low, high;
359 TH1* tmp = Merge(cen, fwd, low, high);
360 TF1* f = FitMerged(tmp, low, high);
361 MakeSysError(tmp, cen, fwd, f);
362 delete f;
363 Info("", "Adding systematic error histogram %s",
364 tmp->GetName());
365 fResults->GetHists()->AddFirst(tmp, "e5");
366 TH1* tmp2 = Symmetrice(tmp);
367 tmp2->SetFillColor(tmp->GetFillColor());
368 tmp2->SetFillStyle(tmp->GetFillStyle());
369 tmp2->SetMarkerStyle(tmp->GetMarkerStyle());
370 tmp2->SetLineWidth(tmp->GetLineWidth());
371 fResults->GetHists()->AddFirst(tmp2, "e5");
372 fResults->Modified();
373 }
374 }
375 delete fwdA;
376 delete cenA;
e308a636 377
378 // --- Close the input file --------------------------------------
379 file->Close();
380
381
b2e7f2d6 382
e308a636 383 // --- Plot results ----------------------------------------------
384 Plot(max, rmax, amax);
b2e7f2d6 385 }
e1f47419 386
b2e7f2d6 387 //__________________________________________________________________
2f92ff0e 388 /**
e1f47419 389 * Fetch the information on the run from the results list
2f92ff0e 390 *
e1f47419 391 * @param results Results list
2f92ff0e 392 */
e1f47419 393 void FetchInformation(const TList* results)
b2e7f2d6 394 {
2f92ff0e 395 if (!fTrigString)
0be6c8cd 396 fTrigString = static_cast<TNamed*>(results->FindObject("trigger"));
397 if (!fNormString)
398 fNormString = static_cast<TNamed*>(results->FindObject("scheme"));
2f92ff0e 399 if (!fSNNString)
400 fSNNString = static_cast<TNamed*>(results->FindObject("sNN"));
401 if (!fSysString)
402 fSysString = static_cast<TNamed*>(results->FindObject("sys"));
403 if (!fVtxAxis)
404 fVtxAxis = static_cast<TAxis*>(results->FindObject("vtxAxis"));
e308a636 405 if (!fCentAxis)
406 fCentAxis = static_cast<TAxis*>(results->FindObject("centAxis"));
c25b5e1b 407
408 TNamed* options = static_cast<TAxis*>(results->FindObject("options"));
0be6c8cd 409 if (!fTrigString) fTrigString = new TNamed("trigger", "unknown");
410 if (!fNormString) fNormString = new TNamed("scheme", "unknown");
bad9a3c1 411 if (!fSNNString) fSNNString = new TNamed("sNN", "unknown");
412 if (!fSysString) fSysString = new TNamed("sys", "unknown");
413 if (!fVtxAxis) {
414 fVtxAxis = new TAxis(1,0,0);
415 fVtxAxis->SetName("vtxAxis");
416 fVtxAxis->SetTitle("v_{z} range unspecified");
417 }
b2e7f2d6 418
0be6c8cd 419 TString centTxt("none");
e308a636 420 if (fCentAxis) {
421 Int_t nCent = fCentAxis->GetNbins();
0be6c8cd 422 centTxt = Form("%d bins", nCent);
e308a636 423 for (Int_t i = 0; i <= nCent; i++)
0be6c8cd 424 centTxt.Append(Form("%c%d", i == 0 ? ' ' : '-',
ffca499d 425 int(fCentAxis->GetXbins()->At(i))));
e308a636 426 }
e1f47419 427 Info("FetchInformation",
fe52e455 428 "Initialized for\n"
0be6c8cd 429 " Trigger: %-30s (%d)\n"
430 " sqrt(sNN): %-30s (%dGeV)\n"
431 " System: %-30s (%d)\n"
432 " Vz range: %-30s (%f,%f)\n"
433 " Normalization: %-30s (%d)\n"
c25b5e1b 434 " Centrality: %s\n"
435 " Options: %s",
fe52e455 436 fTrigString->GetTitle(), fTrigString->GetUniqueID(),
437 fSNNString->GetTitle(), fSNNString->GetUniqueID(),
438 fSysString->GetTitle(), fSysString->GetUniqueID(),
e308a636 439 fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax(),
0be6c8cd 440 fNormString->GetTitle(), fNormString->GetUniqueID(),
c25b5e1b 441 centTxt.Data(), (options ? options->GetTitle() : "none"));
e1f47419 442 }
e308a636 443 //__________________________________________________________________
444 TMultiGraph* FetchOthers(UShort_t centLow, UShort_t centHigh)
445 {
446 TMultiGraph* thisOther = 0;
797161e8 447 if (fShowOthers == 0) return 0;
e308a636 448
e308a636 449 UShort_t sys = (fSysString ? fSysString->GetUniqueID() : 0);
450 UShort_t trg = (fTrigString ? fTrigString->GetUniqueID() : 0);
451 UShort_t snn = (fSNNString ? fSNNString->GetUniqueID() : 0);
797161e8 452 Long_t ret = gROOT->ProcessLine(Form("GetData(%d,%d,%d,%d,%d,%d);",
e308a636 453 sys,snn,trg,
797161e8 454 centLow,centHigh,
455 fShowOthers));
456 if (!ret) return 0;
457
458 thisOther = reinterpret_cast<TMultiGraph*>(ret);
e308a636 459 return thisOther;
460 }
e1f47419 461 //__________________________________________________________________
462 /**
e308a636 463 * Get the results from the top-level list
e1f47419 464 *
e308a636 465 * @param list List
466 * @param name name
467 * @param max On return, maximum of data
468 * @param rmax On return, maximum of ratios
469 * @param amax On return, maximum of left-right comparisons
e1f47419 470 */
797161e8 471 TObjArray*
472 FetchResults(const TList* list,
473 const char* name,
474 Double_t& max,
475 Double_t& rmax,
476 Double_t& amax)
e1f47419 477 {
797161e8 478 UShort_t n = fCentAxis ? fCentAxis->GetNbins() : 0;
e308a636 479 if (n == 0) {
480 TList* all = static_cast<TList*>(list->FindObject("all"));
797161e8 481 if (!all) {
e308a636 482 Error("FetchResults", "Couldn't find list 'all' in %s",
483 list->GetName());
797161e8 484 return 0;
485 }
486 TObjArray* a = new TObjArray;
487 TH1* h = FetchResults(all, name, FetchOthers(0,0),
488 -1000, 0, max, rmax, amax);
489 a->AddAt(h, 0);
490 return a;
b2e7f2d6 491 }
e1f47419 492
797161e8 493 TObjArray* a = new TObjArray;
e308a636 494 for (UShort_t i = 0; i < n; i++) {
5bb5d1f6 495 UShort_t centLow = fCentAxis->GetBinLowEdge(i+1);
496 UShort_t centHigh = fCentAxis->GetBinUpEdge(i+1);
e308a636 497 TString lname = Form("cent%03d_%03d", centLow, centHigh);
498 TList* thisCent = static_cast<TList*>(list->FindObject(lname));
5bb5d1f6 499 Int_t col = GetCentralityColor(i+1);
e308a636 500
501 TString centTxt = Form("%3d%%-%3d%% central", centLow, centHigh);
797161e8 502 if (!thisCent) {
e308a636 503 Error("FetchResults", "Couldn't find list '%s' in %s",
504 lname.Data(), list->GetName());
797161e8 505 continue;
506 }
507 TH1* h = FetchResults(thisCent, name, FetchOthers(centLow, centHigh),
508 col, centTxt.Data(), max, rmax, amax);
509 a->AddAt(h, i);
e1f47419 510 }
797161e8 511 return a;
e308a636 512 }
513 //__________________________________________________________________
5bb5d1f6 514 Int_t GetCentralityColor(Int_t bin) const
515 {
516 UShort_t centLow = fCentAxis->GetBinLowEdge(bin);
517 UShort_t centHigh = fCentAxis->GetBinUpEdge(bin);
518 Float_t fc = (centLow+double(centHigh-centLow)/2) / 100;
519 Int_t nCol = gStyle->GetNumberOfColors();
520 Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5));
521 Int_t col = gStyle->GetColorPalette(icol);
522 //Info("GetCentralityColor","%3d: %3d-%3d -> %3d",bin,centLow,centHigh,col);
523 return col;
524 }
525 //__________________________________________________________________
e308a636 526 void SetAttributes(TH1* h, Int_t color)
527 {
528 if (!h) return;
529 if (color < 0) return;
ffca499d 530 // h->SetLineColor(color);
e308a636 531 h->SetMarkerColor(color);
ffca499d 532 // h->SetFillColor(color);
e308a636 533 }
534 //__________________________________________________________________
535 void SetAttributes(TGraph* g, Int_t color)
536 {
537 if (!g) return;
538 if (color < 0) return;
ffca499d 539 // g->SetLineColor(color);
e308a636 540 g->SetMarkerColor(color);
ffca499d 541 // g->SetFillColor(color);
b2e7f2d6 542 }
543 //__________________________________________________________________
9ecab72f 544 void ModifyTitle(TNamed* /*h*/, const char* /*centTxt*/)
e308a636 545 {
5bb5d1f6 546 return;
547 // if (!centTxt || !h) return;
548 // h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
e308a636 549 }
550
551 //__________________________________________________________________
552 /**
553 * Fetch results for a particular centrality bin
554 *
555 * @param list List
556 * @param name Name
ffca499d 557 * @param thisOther Other graphs
558 * @param color Color
559 * @param centTxt Centrality text
e308a636 560 * @param max On return, data maximum
561 * @param rmax On return, ratio maximum
562 * @param amax On return, left-right maximum
563 */
797161e8 564 TH1* FetchResults(const TList* list,
e308a636 565 const char* name,
566 TMultiGraph* thisOther,
567 Int_t color,
568 const char* centTxt,
569 Double_t& max,
570 Double_t& rmax,
571 Double_t& amax)
b2e7f2d6 572 {
e1f47419 573 TH1* dndeta = FetchResult(list, Form("dndeta%s", name));
574 TH1* dndetaMC = FetchResult(list, Form("dndeta%sMC", name));
575 TH1* dndetaTruth = FetchResult(list, "dndetaTruth");
576 TH1* dndetaSym = 0;
577 TH1* dndetaMCSym = 0;
e308a636 578 SetAttributes(dndeta, color);
5bb5d1f6 579 SetAttributes(dndetaMC, fCentAxis ? color : color+2);
e308a636 580 SetAttributes(dndetaTruth,color);
581 SetAttributes(dndetaSym, color);
5bb5d1f6 582 SetAttributes(dndetaMCSym,fCentAxis ? color : color+2);
5bb5d1f6 583 if (dndetaMC && fCentAxis)
584 dndetaMC->SetMarkerStyle(dndetaMC->GetMarkerStyle()+2);
585 if (dndetaMCSym && fCentAxis)
586 dndetaMCSym->SetMarkerStyle(dndetaMCSym->GetMarkerStyle()+2);
9ecab72f 587 if (dndetaTruth && fCentAxis) {
5bb5d1f6 588 dndetaTruth->SetMarkerStyle(34);
9ecab72f 589 dndetaTruth->SetMarkerColor(kYellow-1);
590 }
e308a636 591 ModifyTitle(dndeta, centTxt);
592 ModifyTitle(dndetaMC, centTxt);
593 ModifyTitle(dndetaTruth,centTxt);
594 ModifyTitle(dndetaSym, centTxt);
595 ModifyTitle(dndetaMCSym,centTxt);
e308a636 596
e1f47419 597
e308a636 598 max = TMath::Max(max, AddHistogram(fResults, dndetaTruth, "e5 p"));
599 max = TMath::Max(max, AddHistogram(fResults, dndetaMC, dndetaMCSym));
600 max = TMath::Max(max, AddHistogram(fResults, dndeta, dndetaSym));
5bb5d1f6 601
9ecab72f 602 if (dndetaTruth) {
603 fTruth = dndetaTruth;
5bb5d1f6 604 }
9ecab72f 605 else {
606 if (fShowRings) {
607 THStack* rings = static_cast<THStack*>(list->FindObject("dndetaRings"));
608 if (rings) {
609 TIter next(rings->GetHists());
610 TH1* hist = 0;
611 while ((hist = static_cast<TH1*>(next())))
612 max = TMath::Max(max, AddHistogram(fResults, hist));
ffca499d 613 }
b2e7f2d6 614 }
9ecab72f 615 // Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f",
616 // dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
617
618 if (fShowLeftRight) {
619 fLeftRight->Add(Asymmetry(dndeta, amax));
620 fLeftRight->Add(Asymmetry(dndetaMC, amax));
9ecab72f 621 }
622
623 if (thisOther) {
624 TIter next(thisOther->GetListOfGraphs());
625 TGraph* g = 0;
626 while ((g = static_cast<TGraph*>(next()))) {
627 fRatios->Add(Ratio(dndeta, g, rmax));
628 fRatios->Add(Ratio(dndetaSym, g, rmax));
629 SetAttributes(g, color);
630 ModifyTitle(g, centTxt);
631 if (!fOthers->GetListOfGraphs() ||
632 !fOthers->GetListOfGraphs()->FindObject(g->GetName())) {
633 max = TMath::Max(max,TMath::MaxElement(g->GetN(), g->GetY()));
634 fOthers->Add(g);
635 }
636 }
637 // fOthers->Add(thisOther);
638 }
b2e7f2d6 639 }
5bb5d1f6 640 if (dndetaMC) {
641 fRatios->Add(Ratio(dndeta, dndetaMC, rmax));
642 fRatios->Add(Ratio(dndetaSym, dndetaMCSym, rmax));
643 }
9ecab72f 644 if (fTruth) {
645 fRatios->Add(Ratio(dndeta, fTruth, rmax));
646 fRatios->Add(Ratio(dndetaSym, fTruth, rmax));
b2e7f2d6 647 }
797161e8 648 return dndeta;
b2e7f2d6 649 }
650 //__________________________________________________________________
2f92ff0e 651 /**
652 * Plot the results
2f92ff0e 653 * @param max Max value
2f92ff0e 654 * @param rmax Maximum diviation from 1 of ratios
2f92ff0e 655 * @param amax Maximum diviation from 1 of asymmetries
656 */
e308a636 657 void Plot(Double_t max,
b2e7f2d6 658 Double_t rmax,
b2e7f2d6 659 Double_t amax)
660 {
661 gStyle->SetOptTitle(0);
662 gStyle->SetTitleFont(132, "xyz");
663 gStyle->SetLabelFont(132, "xyz");
664
e308a636 665 Int_t h = 800;
666 Int_t w = 800; // h / TMath::Sqrt(2);
667 Double_t y1 = 0;
668 Double_t y2 = 0;
669 Double_t y3 = 0;
797161e8 670 if (!fShowRatios) w *= 1.3;
e308a636 671 else y1 = 0.3;
797161e8 672 if (!fShowLeftRight) w *= 1.3;
e308a636 673 else {
674 Double_t y11 = y1;
675 y1 = (y11 > 0.0001 ? 0.4 : 0.2);
ffca499d 676 y2 = (y11 > 0.0001 ? 0.2 : 0.3);
e308a636 677 }
b2e7f2d6 678 TCanvas* c = new TCanvas("Results", "Results", w, h);
679 c->SetFillColor(0);
680 c->SetBorderSize(0);
681 c->SetBorderMode(0);
682
e308a636 683 PlotResults(max, y1);
b2e7f2d6 684 c->cd();
685
e308a636 686 PlotRatios(rmax, y2, y1);
b2e7f2d6 687 c->cd( );
688
e308a636 689 PlotLeftRight(amax, y3, y2);
b2e7f2d6 690 c->cd();
bad9a3c1 691
692
693 Int_t vMin = fVtxAxis->GetXmin();
694 Int_t vMax = fVtxAxis->GetXmax();
695 TString trg(fTrigString->GetTitle());
e308a636 696 Int_t nev = 0;
ffca499d 697 if (fTriggers) nev = fTriggers->GetBinContent(1);
bad9a3c1 698 trg = trg.Strip(TString::kBoth);
f709afcf 699 trg.ReplaceAll(" ", "_");
700 trg.ReplaceAll(">", "Gt");
bad9a3c1 701 TString base(Form("dndeta_%s_%s_%s_%c%02d%c%02dcm_%09dev",
702 fSysString->GetTitle(),
703 fSNNString->GetTitle(),
704 trg.Data(),
705 vMin < 0 ? 'm' : 'p', TMath::Abs(vMin),
706 vMax < 0 ? 'm' : 'p', TMath::Abs(vMax),
707 nev));
708 c->SaveAs(Form("%s.png", base.Data()));
709 c->SaveAs(Form("%s.root", base.Data()));
710 c->SaveAs(Form("%s.C", base.Data()));
33ab4a11 711 base.ReplaceAll("dndeta", "export");
712 Export(base);
b2e7f2d6 713 }
714 //__________________________________________________________________
e308a636 715 /**
716 * Build main legend
717 *
ffca499d 718 * @param stack Stack to include
719 * @param mg (optional) Multi graph to include
720 * @param x1 Lower X coordinate in the range [0,1]
721 * @param y1 Lower Y coordinate in the range [0,1]
722 * @param x2 Upper X coordinate in the range [0,1]
723 * @param y2 Upper Y coordinate in the range [0,1]
e308a636 724 */
ffca499d 725 void BuildLegend(THStack* stack, TMultiGraph* mg,
726 Double_t x1, Double_t y1, Double_t x2, Double_t y2)
e308a636 727 {
728 TLegend* l = new TLegend(x1,y1,x2,y2);
5bb5d1f6 729 l->SetNColumns(fCentAxis ? 1 : 2);
e308a636 730 l->SetFillColor(0);
731 l->SetFillStyle(0);
732 l->SetBorderSize(0);
733 l->SetTextFont(132);
734
797161e8 735 // Loop over items in stack and get unique items, while ignoring
736 // mirrored data and systematic error bands
ffca499d 737 TIter next(stack->GetHists());
5bb5d1f6 738 TH1* hist = 0;
739 TObjArray unique;
740 unique.SetOwner();
797161e8 741 Bool_t sysErrSeen = false;
5bb5d1f6 742 while ((hist = static_cast<TH1*>(next()))) {
797161e8 743 TString t(hist->GetTitle());
744 TString n(hist->GetName());
745 n.ToLower();
746 if (t.Contains("mirrored")) continue;
747 if (n.Contains("syserror")) { sysErrSeen = true; continue; }
748 if (unique.FindObject(t.Data())) continue;
749 TObjString* s = new TObjString(hist->GetTitle());
5bb5d1f6 750 s->SetUniqueID(((hist->GetMarkerStyle() & 0xFFFF) << 16) |
751 ((hist->GetMarkerColor() & 0xFFFF) << 0));
752 unique.Add(s);
753 // l->AddEntry(hist, hist->GetTitle(), "pl");
e308a636 754 }
ffca499d 755 if (mg) {
797161e8 756 // If we have other stuff, scan for unique names
ffca499d 757 TIter nexto(mg->GetListOfGraphs());
5bb5d1f6 758 TGraph* g = 0;
759 while ((g = static_cast<TGraph*>(nexto()))) {
760 TString n(g->GetTitle());
ffca499d 761 if (n.Contains("mirrored")) continue;
5bb5d1f6 762 if (unique.FindObject(n.Data())) continue;
763 TObjString* s = new TObjString(n);
764 s->SetUniqueID(((g->GetMarkerStyle() & 0xFFFF) << 16) |
765 ((g->GetMarkerColor() & 0xFFFF) << 0));
766 unique.Add(s);
767 // l->AddEntry(hist, hist->GetTitle(), "pl");
ffca499d 768 }
e308a636 769 }
797161e8 770
771 // Add legend entries for unique items only
5bb5d1f6 772 TIter nextu(&unique);
773 TObject* s = 0;
774 Int_t i = 0;
775 while ((s = nextu())) {
776 TLegendEntry* dd = l->AddEntry(Form("data%2d", i++),
777 s->GetName(), "lp");
778 Int_t style = (s->GetUniqueID() >> 16) & 0xFFFF;
779 Int_t color = (s->GetUniqueID() >> 0) & 0xFFFF;
780 dd->SetLineColor(kBlack);
781 if (fCentAxis) dd->SetMarkerColor(kBlack);
782 else dd->SetMarkerColor(color);
783 dd->SetMarkerStyle(style);
784 }
797161e8 785 if (sysErrSeen) {
786 // Add entry for systematic errors
787 TLegendEntry* d0 = l->AddEntry("d0", Form("%4.1f%% Systematic error",
788 100*fFwdSysErr), "f");
e42bc740 789 d0->SetLineColor(SYSERR_COLOR);
790 d0->SetMarkerColor(SYSERR_COLOR);
791 d0->SetFillColor(SYSERR_COLOR);
792 d0->SetFillStyle(SYSERR_STYLE);
797161e8 793 d0->SetMarkerStyle(0);
794 d0->SetLineWidth(0);
795 i++;
796 }
797 if (!fCentAxis && i % 2 == 1) {
798 // To make sure the 'data' and 'mirrored' entries are on a line
799 // by themselves
800 TLegendEntry* dd = l->AddEntry("dd", " ", "");
801 dd->SetTextSize(0);
802 dd->SetFillStyle(0);
803 dd->SetFillColor(0);
804 dd->SetLineWidth(0);
805 dd->SetLineColor(0);
806 dd->SetMarkerSize(0);
807 }
808 // Add entry for 'data'
e308a636 809 TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp");
810 d1->SetLineColor(kBlack);
811 d1->SetMarkerColor(kBlack);
812 d1->SetMarkerStyle(20);
797161e8 813 // Add entry for 'mirrored data'
e308a636 814 TLegendEntry* d2 = l->AddEntry("d2", "Mirrored data", "lp");
815 d2->SetLineColor(kBlack);
816 d2->SetMarkerColor(kBlack);
817 d2->SetMarkerStyle(24);
818
819 l->Draw();
820 }
821 //__________________________________________________________________
5bb5d1f6 822 /**
823 * Build centrality legend
824 *
825 * @param axis Stack to include
826 * @param x1 Lower X coordinate in the range [0,1]
827 * @param y1 Lower Y coordinate in the range [0,1]
828 * @param x2 Upper X coordinate in the range [0,1]
829 * @param y2 Upper Y coordinate in the range [0,1]
830 */
831 void BuildCentLegend(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
832 {
833 if (!fCentAxis) return;
834
835 TLegend* l = new TLegend(x1,y1,x2,y2);
836 l->SetNColumns(1);
837 l->SetFillColor(0);
838 l->SetFillStyle(0);
839 l->SetBorderSize(0);
840 l->SetTextFont(132);
841
842 Int_t n = fCentAxis->GetNbins();
843 for (Int_t i = 1; i <= n; i++) {
844 Double_t low = fCentAxis->GetBinLowEdge(i);
845 Double_t upp = fCentAxis->GetBinUpEdge(i);
846 TLegendEntry* e = l->AddEntry(Form("dummy%02d", i),
847 Form("%3d%% - %3d%%",
848 int(low), int(upp)), "pl");
849 e->SetMarkerColor(GetCentralityColor(i));
850 }
851 l->Draw();
852 }
853 //__________________________________________________________________
2f92ff0e 854 /**
855 * Plot the results
856 *
2f92ff0e 857 * @param max Maximum
858 * @param yd Bottom position of pad
859 */
e308a636 860 void PlotResults(Double_t max, Double_t yd)
b2e7f2d6 861 {
862 // Make a sub-pad for the result itself
863 TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0);
864 p1->SetTopMargin(0.05);
865 p1->SetBorderSize(0);
866 p1->SetBorderMode(0);
867 p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1);
797161e8 868 p1->SetRightMargin(0.03);
869 if (fShowLeftRight || fShowRatios) p1->SetGridx();
b2e7f2d6 870 p1->SetTicks(1,1);
871 p1->SetNumber(1);
872 p1->Draw();
873 p1->cd();
33ab4a11 874
e308a636 875 // Info("PlotResults", "Plotting results with max=%f", max);
876 fResults->SetMaximum(1.15*max);
877 fResults->SetMinimum(yd > 0.00001 ? -0.1 : 0);
b2e7f2d6 878
797161e8 879 FixAxis(fResults, (1-yd)*(yd > .001 ? 1 : .9 / 1.2),
880 "#font[12]{#frac{1}{N} "
881 "#frac{dN_{#font[132]{ch}}}{d#font[152]{#eta}}}");
b2e7f2d6 882
883 p1->Clear();
e308a636 884 fResults->DrawClone("nostack e1");
b2e7f2d6 885
e308a636 886 fRangeParam->fSlave1Axis = fResults->GetXaxis();
b2e7f2d6 887 fRangeParam->fSlave1Pad = p1;
888
889 // Draw other data
797161e8 890 if (fShowOthers != 0) {
b2e7f2d6 891 TGraphAsymmErrors* o = 0;
e308a636 892 TIter nextG(fOthers->GetListOfGraphs());
b2e7f2d6 893 while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
894 o->DrawClone("same p");
895 }
896
897 // Make a legend in the main result pad
5bb5d1f6 898 BuildCentLegend(.12, 1-p1->GetTopMargin()-.01-.5,
899 .35, 1-p1->GetTopMargin()-.01-.1);
900 Double_t x1 = (fCentAxis ? .7 : .15);
901 Double_t x2 = (fCentAxis ? 1-p1->GetRightMargin()-.01: .90);
902 Double_t y1 = (fCentAxis ? .5: p1->GetBottomMargin()+.01);
903 Double_t y2 = (fCentAxis ? 1-p1->GetTopMargin()-.01-.15 : .35);
904
905 BuildLegend(fResults, fOthers, x1, y1, x2, y2);
b2e7f2d6 906
907 // Put a title on top
ffca499d 908 fTitle.ReplaceAll("@", " ");
b2e7f2d6 909 TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
910 tit->SetNDC();
911 tit->SetTextFont(132);
912 tit->SetTextSize(0.05);
913 tit->Draw();
914
f709afcf 915 Int_t aliceBlue = TColor::GetColor(41,73,156);
b2e7f2d6 916 // Put a nice label in the plot
917 TString eS;
bad9a3c1 918 UShort_t snn = fSNNString->GetUniqueID();
919 const char* sys = fSysString->GetTitle();
5bb5d1f6 920 if (snn == 2750) snn = 2760;
b2e7f2d6 921 if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
922 else eS = Form("%3dGeV", snn);
5bb5d1f6 923 TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s%s}=%s, %s",
b2e7f2d6 924 sys,
5bb5d1f6 925 (fCentAxis ? "_{NN}" : ""),
b2e7f2d6 926 eS.Data(),
5bb5d1f6 927 fCentAxis ? "by centrality" :
bad9a3c1 928 fTrigString->GetTitle()));
f709afcf 929 tt->SetTextColor(aliceBlue);
b2e7f2d6 930 tt->SetNDC();
931 tt->SetTextFont(132);
932 tt->SetTextAlign(33);
933 tt->Draw();
934
935 // Put number of accepted events on the plot
e308a636 936 Int_t nev = 0;
ffca499d 937 if (fTriggers) nev = fTriggers->GetBinContent(1);
b2e7f2d6 938 TLatex* et = new TLatex(.93, .83, Form("%d events", nev));
f709afcf 939 et->SetTextColor(aliceBlue);
b2e7f2d6 940 et->SetNDC();
941 et->SetTextFont(132);
942 et->SetTextAlign(33);
943 et->Draw();
944
945 // Put number of accepted events on the plot
946 if (fVtxAxis) {
947 TLatex* vt = new TLatex(.93, .88, fVtxAxis->GetTitle());
948 vt->SetNDC();
949 vt->SetTextFont(132);
950 vt->SetTextAlign(33);
f709afcf 951 vt->SetTextColor(aliceBlue);
b2e7f2d6 952 vt->Draw();
953 }
954 // results->Draw("nostack e1 same");
955
e308a636 956 fRangeParam->fSlave1Axis = FindXAxis(p1, fResults->GetName());
b2e7f2d6 957 fRangeParam->fSlave1Pad = p1;
958
959
960 // Mark the plot as preliminary
f709afcf 961 TLatex* pt = new TLatex(.12, .93, "Work in progress");
b2e7f2d6 962 pt->SetNDC();
963 pt->SetTextFont(22);
f709afcf 964 // pt->SetTextSize();
965 pt->SetTextColor(TColor::GetColor(234,26,46));
b2e7f2d6 966 pt->SetTextAlign(13);
967 pt->Draw();
968
f709afcf 969 const char* logos[] = { "ALICE.png", "FMD.png", 0 };
970 const char** logo = logos;
971 while (*logo) {
972 if (gSystem->AccessPathName(*logo)) {
973 logo++;
974 continue;
975 }
797161e8 976 TPad* pad = new TPad("logo", "logo", .12, .7, .25, .9, 0, 0, 0);
f709afcf 977 pad->SetFillStyle(0);
978 pad->Draw();
979 pad->cd();
b2e7f2d6 980 TImage* i = TImage::Create();
f709afcf 981 i->ReadImage(*logo);
b2e7f2d6 982 i->Draw();
f709afcf 983 break;
b2e7f2d6 984 }
985 p1->cd();
986 }
987 //__________________________________________________________________
2f92ff0e 988 /**
989 * Plot the ratios
990 *
2f92ff0e 991 * @param max Maximum diviation from 1
992 * @param y1 Lower y coordinate of pad
993 * @param y2 Upper y coordinate of pad
994 */
e308a636 995 void PlotRatios(Double_t max, Double_t y1, Double_t y2)
b2e7f2d6 996 {
797161e8 997 if (fShowRatios == 0) return;
e308a636 998
b2e7f2d6 999 bool isBottom = (y1 < 0.0001);
1000 Double_t yd = y2 - y1;
1001 // Make a sub-pad for the result itself
1002 TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0);
1003 p2->SetTopMargin(0.001);
797161e8 1004 p2->SetRightMargin(0.03);
b2e7f2d6 1005 p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
1006 p2->SetGridx();
1007 p2->SetTicks(1,1);
1008 p2->SetNumber(2);
1009 p2->Draw();
1010 p2->cd();
1011
1012 // Fix up axis
797161e8 1013 FixAxis(fRatios, yd, "Ratios", 7);
b2e7f2d6 1014
e308a636 1015 fRatios->SetMaximum(1+TMath::Max(.22,1.05*max));
1016 fRatios->SetMinimum(1-TMath::Max(.32,1.05*max));
b2e7f2d6 1017 p2->Clear();
e308a636 1018 fRatios->DrawClone("nostack e1");
b2e7f2d6 1019
1020
1021 // Make a legend
ffca499d 1022 BuildLegend(fRatios, 0, .15,p2->GetBottomMargin()+.01,.9,
1023 isBottom ? .6 : .4);
1024#if 0
b2e7f2d6 1025 TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,
1026 isBottom ? .6 : .4);
1027 l2->SetNColumns(2);
1028 l2->SetFillColor(0);
1029 l2->SetFillStyle(0);
1030 l2->SetBorderSize(0);
1031 l2->SetTextFont(132);
ffca499d 1032#endif
b2e7f2d6 1033 // Make a nice band from 0.9 to 1.1
1034 TGraphErrors* band = new TGraphErrors(2);
e1f47419 1035 band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
1036 band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
b2e7f2d6 1037 band->SetPointError(0, 0, .1);
1038 band->SetPointError(1, 0, .1);
1039 band->SetFillColor(kYellow+2);
1040 band->SetFillStyle(3002);
1041 band->SetLineStyle(2);
1042 band->SetLineWidth(1);
1043 band->Draw("3 same");
1044 band->DrawClone("X L same");
1045
1046 // Replot the ratios on top
e308a636 1047 fRatios->DrawClone("nostack e1 same");
b2e7f2d6 1048
1049 if (isBottom) {
e308a636 1050 fRangeParam->fMasterAxis = FindXAxis(p2, fRatios->GetName());
b2e7f2d6 1051 p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
1052 fRangeParam));
1053 }
1054 else {
e308a636 1055 fRangeParam->fSlave2Axis = FindXAxis(p2, fRatios->GetName());
b2e7f2d6 1056 fRangeParam->fSlave2Pad = p2;
1057 }
1058 }
1059 //__________________________________________________________________
2f92ff0e 1060 /**
1061 * Plot the asymmetries
1062 *
2f92ff0e 1063 * @param max Maximum diviation from 1
1064 * @param y1 Lower y coordinate of pad
1065 * @param y2 Upper y coordinate of pad
1066 */
e308a636 1067 void PlotLeftRight(Double_t max, Double_t y1, Double_t y2)
b2e7f2d6 1068 {
e308a636 1069 if (!fShowLeftRight) return;
1070
b2e7f2d6 1071 bool isBottom = (y1 < 0.0001);
1072 Double_t yd = y2 - y1;
1073 // Make a sub-pad for the result itself
1074 TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0);
1075 p3->SetTopMargin(0.001);
797161e8 1076 p3->SetRightMargin(0.03);
b2e7f2d6 1077 p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
1078 p3->SetGridx();
1079 p3->SetTicks(1,1);
1080 p3->SetNumber(2);
1081 p3->Draw();
1082 p3->cd();
1083
b2e7f2d6 1084 // Fix up axis
797161e8 1085 FixAxis(fLeftRight, yd, "Right/Left", 4);
b2e7f2d6 1086
e308a636 1087 fLeftRight->SetMaximum(1+TMath::Max(.12,1.05*max));
1088 fLeftRight->SetMinimum(1-TMath::Max(.15,1.05*max));
b2e7f2d6 1089 p3->Clear();
e308a636 1090 fLeftRight->DrawClone("nostack e1");
b2e7f2d6 1091
1092
1093 // Make a legend
5bb5d1f6 1094 Double_t xx1 = (fCentAxis ? .7 : .15);
1095 Double_t xx2 = (fCentAxis ? 1-p3->GetRightMargin()-.01 : .90);
1096 Double_t yy1 = p3->GetBottomMargin()+.01;
1097 Double_t yy2 = (fCentAxis ? 1-p3->GetTopMargin()-.01-.15 : .5);
1098 BuildLegend(fLeftRight, 0, xx1, yy1, xx2, yy2);
1099 // TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
1100 // l2->SetNColumns(2);
1101 // l2->SetFillColor(0);
1102 // l2->SetFillStyle(0);
1103 // l2->SetBorderSize(0);
1104 // l2->SetTextFont(132);
ffca499d 1105
b2e7f2d6 1106 // Make a nice band from 0.9 to 1.1
1107 TGraphErrors* band = new TGraphErrors(2);
e1f47419 1108 band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
1109 band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
b2e7f2d6 1110 band->SetPointError(0, 0, .05);
1111 band->SetPointError(1, 0, .05);
1112 band->SetFillColor(kYellow+2);
1113 band->SetFillStyle(3002);
1114 band->SetLineStyle(2);
1115 band->SetLineWidth(1);
1116 band->Draw("3 same");
1117 band->DrawClone("X L same");
1118
e308a636 1119 fLeftRight->DrawClone("nostack e1 same");
b2e7f2d6 1120 if (isBottom) {
e308a636 1121 fRangeParam->fMasterAxis = FindXAxis(p3, fLeftRight->GetName());
b2e7f2d6 1122 p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
1123 fRangeParam));
1124 }
1125 else {
e308a636 1126 fRangeParam->fSlave2Axis = FindXAxis(p3, fLeftRight->GetName());
b2e7f2d6 1127 fRangeParam->fSlave2Pad = p3;
1128 }
1129 }
2f92ff0e 1130 /** @} */
1131 //==================================================================
1132 /**
1133 * @{
1134 * @name Data utility functions
1135 */
1136 /**
1137 * Get a result from the passed list
1138 *
1139 * @param list List to search
1140 * @param name Object name to search for
1141 *
1142 * @return
1143 */
e1f47419 1144 TH1* FetchResult(const TList* list, const char* name) const
2f92ff0e 1145 {
1146 if (!list) return 0;
1147
e308a636 1148 TH1* ret = static_cast<TH1*>(list->FindObject(name));
5bb5d1f6 1149#if 0
e1f47419 1150 if (!ret) {
1151 // all->ls();
1152 Warning("GetResult", "Histogram %s not found", name);
2f92ff0e 1153 }
5bb5d1f6 1154#endif
b2e7f2d6 1155
e1f47419 1156 return ret;
2f92ff0e 1157 }
b2e7f2d6 1158 //__________________________________________________________________
2f92ff0e 1159 /**
1160 * Add a histogram to the stack after possibly rebinning it
1161 *
1162 * @param stack Stack to add to
1163 * @param hist histogram
1164 * @param option Draw options
1165 *
1166 * @return Maximum of histogram
1167 */
e308a636 1168 Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option="") const
b2e7f2d6 1169 {
1170 // Check if we have input
1171 if (!hist) return 0;
1172
1173 // Rebin if needed
1174 Rebin(hist);
1175
1176 stack->Add(hist, option);
1177 return hist->GetMaximum();
1178 }
1179 //__________________________________________________________________
2f92ff0e 1180 /**
1181 * Add a histogram to the stack after possibly rebinning it
1182 *
1183 * @param stack Stack to add to
1184 * @param hist histogram
1185 * @param option Draw options
1186 * @param sym On return, the data symmetriced (added to stack)
1187 *
1188 * @return Maximum of histogram
1189 */
e308a636 1190 Double_t AddHistogram(THStack* stack, TH1* hist, TH1*& sym,
1191 Option_t* option="") const
b2e7f2d6 1192 {
1193 // Check if we have input
1194 if (!hist) return 0;
1195
1196 // Rebin if needed
1197 Rebin(hist);
1198 stack->Add(hist, option);
1199
1200 // Now symmetrice the histogram
1201 sym = Symmetrice(hist);
1202 stack->Add(sym, option);
1203
1204 return hist->GetMaximum();
1205 }
1206
1207 //__________________________________________________________________
1208 /**
1209 * Rebin a histogram
1210 *
1211 * @param h Histogram to rebin
b2e7f2d6 1212 *
1213 * @return
1214 */
1215 virtual void Rebin(TH1* h) const
1216 {
1217 if (fRebin <= 1) return;
1218
1219 Int_t nBins = h->GetNbinsX();
1220 if(nBins % fRebin != 0) {
1221 Warning("Rebin", "Rebin factor %d is not a devisor of current number "
1222 "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
1223 return;
1224 }
1225
1226 // Make a copy
1227 TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
1228 tmp->Rebin(fRebin);
1229 tmp->SetDirectory(0);
e28f5fc5 1230 tmp->Reset();
b2e7f2d6 1231 // The new number of bins
1232 Int_t nBinsNew = nBins / fRebin;
1233 for(Int_t i = 1;i<= nBinsNew; i++) {
1234 Double_t content = 0;
1235 Double_t sumw = 0;
1236 Double_t wsum = 0;
1237 Int_t nbins = 0;
1238 for(Int_t j = 1; j<=fRebin;j++) {
1239 Int_t bin = (i-1)*fRebin + j;
1240 Double_t c = h->GetBinContent(bin);
1241
1242 if (c <= 0) continue;
1243
1244 if (fCutEdges) {
1245 if (h->GetBinContent(bin+1)<=0 ||
1246 h->GetBinContent(bin-1)) {
1247 Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
1248 bin, c, h->GetName(),
1249 bin+1, h->GetBinContent(bin+1),
1250 bin-1, h->GetBinContent(bin-1));
1251 continue;
1252 }
1253 }
1254 Double_t e = h->GetBinError(bin);
1255 Double_t w = 1 / (e*e); // 1/c/c
1256 content += c;
1257 sumw += w;
1258 wsum += w * c;
1259 nbins++;
1260 }
1261
e28f5fc5 1262 if(content > 0 && nbins > 1 ) {
b2e7f2d6 1263 tmp->SetBinContent(i, wsum / sumw);
1264 tmp->SetBinError(i,1./TMath::Sqrt(sumw));
1265 }
1266 }
1267
1268 // Finally, rebin the histogram, and set new content
1269 h->Rebin(fRebin);
1270 h->Reset();
1271 for(Int_t i = 1; i<= nBinsNew; i++) {
1272 h->SetBinContent(i,tmp->GetBinContent(i));
1273 h->SetBinError(i, tmp->GetBinError(i));
1274 }
1275
1276 delete tmp;
1277 }
1278 //__________________________________________________________________
1279 /**
1280 * Make an extension of @a h to make it symmetric about 0
1281 *
1282 * @param h Histogram to symmertrice
1283 *
1284 * @return Symmetric extension of @a h
1285 */
1286 TH1* Symmetrice(const TH1* h) const
1287 {
1288 Int_t nBins = h->GetNbinsX();
1289 TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
1290 s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
1291 s->Reset();
1292 s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
1293 s->SetMarkerColor(h->GetMarkerColor());
1294 s->SetMarkerSize(h->GetMarkerSize());
1295 s->SetMarkerStyle(h->GetMarkerStyle()+4);
1296 s->SetFillColor(h->GetFillColor());
1297 s->SetFillStyle(h->GetFillStyle());
1298 s->SetDirectory(0);
1299
1300 // Find the first and last bin with data
1301 Int_t first = nBins+1;
1302 Int_t last = 0;
1303 for (Int_t i = 1; i <= nBins; i++) {
1304 if (h->GetBinContent(i) <= 0) continue;
1305 first = TMath::Min(first, i);
1306 last = TMath::Max(last, i);
1307 }
1308
1309 Double_t xfirst = h->GetBinCenter(first-1);
1310 Int_t f1 = h->GetXaxis()->FindBin(-xfirst);
1311 Int_t l2 = s->GetXaxis()->FindBin(xfirst);
1312 for (Int_t i = f1, j=l2; i <= last; i++,j--) {
1313 s->SetBinContent(j, h->GetBinContent(i));
1314 s->SetBinError(j, h->GetBinError(i));
1315 }
1316 // Fill in overlap bin
1317 s->SetBinContent(l2+1, h->GetBinContent(first));
1318 s->SetBinError(l2+1, h->GetBinError(first));
1319 return s;
1320 }
1321 //__________________________________________________________________
2f92ff0e 1322 /**
1323 * Calculate the left-right asymmetry of input histogram
1324 *
1325 * @param h Input histogram
1326 * @param max On return, the maximum distance from 1 of the histogram
1327 *
1328 * @return Asymmetry
1329 */
1330 TH1* Asymmetry(TH1* h, Double_t& max)
1331 {
1332 if (!h) return 0;
1333
1334 TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName())));
1335 // Int_t oBins = h->GetNbinsX();
1336 // Double_t high = h->GetXaxis()->GetXmax();
1337 // Double_t low = h->GetXaxis()->GetXmin();
1338 // Double_t dBin = (high - low) / oBins;
1339 // Int_t tBins = Int_t(2*high/dBin+.5);
1340 // ret->SetBins(tBins, -high, high);
e308a636 1341 ret->SetDirectory(0);
2f92ff0e 1342 ret->Reset();
1343 ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
1344 ret->SetYTitle("Right/Left");
1345 Int_t nBins = h->GetNbinsX();
1346 for (Int_t i = 1; i <= nBins; i++) {
1347 Double_t x = h->GetBinCenter(i);
1348 if (x > 0) break;
1349
1350 Double_t c1 = h->GetBinContent(i);
1351 Double_t e1 = h->GetBinError(i);
1352 if (c1 <= 0) continue;
1353
1354 Int_t j = h->FindBin(-x);
1355 if (j <= 0 || j > nBins) continue;
1356
1357 Double_t c2 = h->GetBinContent(j);
1358 Double_t e2 = h->GetBinError(j);
1359
1360 Double_t c12 = c1*c1;
1361 Double_t e = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12));
1362
1363 Int_t k = ret->FindBin(x);
1364 ret->SetBinContent(k, c2/c1);
1365 ret->SetBinError(k, e);
1366 }
1367 max = TMath::Max(max, RatioMax(ret));
1368
1369 return ret;
1370 }
1371 //__________________________________________________________________
1372 /**
1373 * Transform a graph into a histogram
1374 *
1375 * @param g
1376 *
1377 * @return
1378 */
1379 TH1* Graph2Hist(const TGraphAsymmErrors* g) const
1380 {
1381 Int_t nBins = g->GetN();
1382 TArrayF bins(nBins+1);
1383 Double_t dx = 0;
1384 for (Int_t i = 0; i < nBins; i++) {
1385 Double_t x = g->GetX()[i];
1386 Double_t exl = g->GetEXlow()[i];
1387 Double_t exh = g->GetEXhigh()[i];
1388 bins.fArray[i] = x-exl;
1389 bins.fArray[i+1] = x+exh;
1390 Double_t dxi = exh+exl;
1391 if (i == 0) dx = dxi;
1392 else if (dxi != dx) dx = 0;
1393 }
1394 TString name(g->GetName());
1395 TString title(g->GetTitle());
1396 TH1D* h = 0;
1397 if (dx != 0) {
1398 h = new TH1D(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
1399 }
1400 else {
1401 h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
1402 }
1403 h->SetMarkerStyle(g->GetMarkerStyle());
1404 h->SetMarkerColor(g->GetMarkerColor());
1405 h->SetMarkerSize(g->GetMarkerSize());
1406
1407 return h;
1408 }
1409 /* @} */
1410 //==================================================================
1411 /**
1412 * @{
1413 * @name Ratio utility functions
1414 */
1415 /**
1416 * Get the maximum diviation from 1 in the passed ratio
1417 *
1418 * @param h Ratio histogram
1419 *
1420 * @return Max diviation from 1
1421 */
b2e7f2d6 1422 Double_t RatioMax(TH1* h) const
1423 {
1424 Int_t nBins = h->GetNbinsX();
1425 Double_t ret = 0;
1426 for (Int_t i = 1; i <= nBins; i++) {
1427 Double_t c = h->GetBinContent(i);
1428 if (c == 0) continue;
1429 Double_t e = h->GetBinError(i);
1430 Double_t d = TMath::Abs(1-c-e);
1431 ret = TMath::Max(d, ret);
1432 }
1433 return ret;
1434 }
1435 //__________________________________________________________________
e1f47419 1436 /**
1437 * Make the ratio of h1 to h2
1438 *
ffca499d 1439 * @param o1 First object (numerator)
1440 * @param o2 Second object (denominator)
1441 * @param max Maximum diviation from 1
e1f47419 1442 *
ffca499d 1443 * @return o1 / o2
e1f47419 1444 */
1445 TH1* Ratio(const TObject* o1, const TObject* o2, Double_t& max) const
1446 {
b30dee70 1447 if (!o1 || !o2) return 0;
1448
5bb5d1f6 1449 TH1* r = 0;
1450 const TAttMarker* m1 = 0;
1451 const TAttMarker* m2 = 0;
e1f47419 1452 const TH1* h1 = dynamic_cast<const TH1*>(o1);
1453 if (h1) {
5bb5d1f6 1454 m1 = h1;
e1f47419 1455 const TH1* h2 = dynamic_cast<const TH1*>(o2);
ffca499d 1456 if (h2) {
5bb5d1f6 1457 m2 = h2;
1458 r = RatioHH(h1,h2);
ffca499d 1459 }
1460 else {
1461 const TGraph* g2 = dynamic_cast<const TGraph*>(o2);
1462 if (g2) {
5bb5d1f6 1463 m2 = g2;
1464 r = RatioHG(h1,g2);
ffca499d 1465 }
1466 }
e1f47419 1467 }
ffca499d 1468 else {
1469 const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(o1);
1470 if (g1) {
5bb5d1f6 1471 m1 = g1;
ffca499d 1472 const TGraphAsymmErrors* g2 =
1473 dynamic_cast<const TGraphAsymmErrors*>(o2);
1474 if (g2) {
5bb5d1f6 1475 m2 = g2;
1476 r = RatioGG(g1, g2);
ffca499d 1477 }
1478 }
e1f47419 1479 }
ffca499d 1480 if (!r) {
1481 Warning("Ratio", "Don't know how to divide a %s (%s) with a %s (%s)",
1482 o1->ClassName(), o1->GetName(), o2->ClassName(), o2->GetName());
1483 return 0;
1484 }
1485 // Check that the histogram isn't empty
1486 if (r->GetEntries() <= 0) {
1487 delete r;
1488 r = 0;
1489 }
5bb5d1f6 1490 if (r) {
9ecab72f 1491 r->SetMarkerStyle(m2->GetMarkerStyle());
1492 r->SetMarkerColor(m1->GetMarkerColor());
5bb5d1f6 1493 r->SetMarkerSize(0.9*m1->GetMarkerSize());
1494 r->SetName(Form("%s_over_%s", o1->GetName(), o2->GetName()));
1495 r->SetTitle(Form("%s / %s", o1->GetTitle(), o2->GetTitle()));
1496 r->SetDirectory(0);
1497 max = TMath::Max(RatioMax(r), max);
1498 }
1499
ffca499d 1500 return r;
e1f47419 1501 }
1502 //__________________________________________________________________
b2e7f2d6 1503 /**
1504 * Compute the ratio of @a h to @a g. @a g is evaluated at the bin
1505 * centers of @a h
1506 *
1507 * @param h Numerator
1508 * @param g Divisor
1509 *
1510 * @return h/g
1511 */
5bb5d1f6 1512 TH1* RatioHG(const TH1* h, const TGraph* g) const
b2e7f2d6 1513 {
1514 if (!h || !g) return 0;
1515
1516 TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
b2e7f2d6 1517 ret->Reset();
b2e7f2d6 1518 Double_t xlow = g->GetX()[0];
1519 Double_t xhigh = g->GetX()[g->GetN()-1];
1520 if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
1521
1522 for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1523 Double_t c = h->GetBinContent(i);
1524 if (c <= 0) continue;
1525
1526 Double_t x = h->GetBinCenter(i);
1527 if (x < xlow || x > xhigh) continue;
1528
1529 Double_t f = g->Eval(x);
1530 if (f <= 0) continue;
1531
1532 ret->SetBinContent(i, c / f);
1533 ret->SetBinError(i, h->GetBinError(i) / f);
1534 }
b2e7f2d6 1535 return ret;
1536 }
1537 //__________________________________________________________________
1538 /**
1539 * Make the ratio of h1 to h2
1540 *
1541 * @param h1 First histogram (numerator)
1542 * @param h2 Second histogram (denominator)
1543 *
1544 * @return h1 / h2
1545 */
5bb5d1f6 1546 TH1* RatioHH(const TH1* h1, const TH1* h2) const
b2e7f2d6 1547 {
1548 if (!h1 || !h2) return 0;
5bb5d1f6 1549 TH1* t1 = static_cast<TH1*>(h1->Clone("tmp"));
b2e7f2d6 1550 t1->Divide(h2);
b2e7f2d6 1551 return t1;
1552 }
1553 //__________________________________________________________________
1554 /**
1555 * Calculate the ratio of two graphs - g1 / g2
1556 *
1557 * @param g1 Numerator
1558 * @param g2 Denominator
1559 *
1560 * @return g1 / g2 in a histogram
1561 */
ffca499d 1562 TH1* RatioGG(const TGraphAsymmErrors* g1,
5bb5d1f6 1563 const TGraphAsymmErrors* g2) const
b2e7f2d6 1564 {
1565 Int_t nBins = g1->GetN();
1566 TArrayF bins(nBins+1);
1567 Double_t dx = 0;
1568 for (Int_t i = 0; i < nBins; i++) {
1569 Double_t x = g1->GetX()[i];
1570 Double_t exl = g1->GetEXlow()[i];
1571 Double_t exh = g1->GetEXhigh()[i];
1572 bins.fArray[i] = x-exl;
1573 bins.fArray[i+1] = x+exh;
1574 Double_t dxi = exh+exl;
1575 if (i == 0) dx = dxi;
1576 else if (dxi != dx) dx = 0;
1577 }
b2e7f2d6 1578 TH1* h = 0;
1579 if (dx != 0) {
5bb5d1f6 1580 h = new TH1F("tmp", "tmp", nBins, bins[0], bins[nBins]);
b2e7f2d6 1581 }
1582 else {
5bb5d1f6 1583 h = new TH1F("tmp", "tmp", nBins, bins.fArray);
b2e7f2d6 1584 }
b2e7f2d6 1585
1586 Double_t low = g2->GetX()[0];
1587 Double_t high = g2->GetX()[g2->GetN()-1];
1588 if (low > high) { Double_t t = low; low = high; high = t; }
1589 for (Int_t i = 0; i < nBins; i++) {
1590 Double_t x = g1->GetX()[i];
1591 if (x < low-dx || x > high+dx) continue;
1592 Double_t c1 = g1->GetY()[i];
1593 Double_t e1 = g1->GetErrorY(i);
1594 Double_t c2 = g2->Eval(x);
1595
1596 h->SetBinContent(i+1, c1 / c2);
1597 h->SetBinError(i+1, e1 / c2);
1598 }
b2e7f2d6 1599 return h;
2f92ff0e 1600 }
1601 /* @} */
1602 //==================================================================
1603 /**
1604 * @{
1605 * @name Graphics utility functions
1606 */
1607 /**
1608 * Find an X axis in a pad
1609 *
1610 * @param p Pad
1611 * @param name Histogram to find axis for
1612 *
1613 * @return Found axis or null
1614 */
1615 TAxis* FindXAxis(TVirtualPad* p, const char* name)
b2e7f2d6 1616 {
2f92ff0e 1617 TObject* o = p->GetListOfPrimitives()->FindObject(name);
1618 if (!o) {
1619 Warning("FindXAxis", "%s not found in pad", name);
1620 return 0;
b2e7f2d6 1621 }
2f92ff0e 1622 THStack* stack = dynamic_cast<THStack*>(o);
1623 if (!stack) {
1624 Warning("FindXAxis", "%s is not a THStack", name);
1625 return 0;
b2e7f2d6 1626 }
2f92ff0e 1627 if (!stack->GetHistogram()) {
1628 Warning("FindXAxis", "%s has no histogram", name);
1629 return 0;
b2e7f2d6 1630 }
2f92ff0e 1631 TAxis* ret = stack->GetHistogram()->GetXaxis();
1632 return ret;
b2e7f2d6 1633 }
2f92ff0e 1634
b2e7f2d6 1635 //__________________________________________________________________
2f92ff0e 1636 /**
1637 * Fix the apperance of the axis in a stack
1638 *
1639 * @param stack stack of histogram
1640 * @param s Scaling factor
1641 * @param ytitle Y axis title
1642 * @param force Whether to draw the stack first or not
1643 * @param ynDiv Divisions on Y axis
1644 */
797161e8 1645 void FixAxis(THStack* stack, Double_t yd, const char* ytitle,
2f92ff0e 1646 Int_t ynDiv=210, Bool_t force=true)
b2e7f2d6 1647 {
e308a636 1648 if (!stack) {
1649 Warning("FixAxis", "No stack passed for %s!", ytitle);
1650 return;
1651 }
2f92ff0e 1652 if (force) stack->Draw("nostack e1");
b2e7f2d6 1653
2f92ff0e 1654 TH1* h = stack->GetHistogram();
e308a636 1655 if (!h) {
1656 Warning("FixAxis", "Stack %s has no histogram", stack->GetName());
1657 return;
1658 }
797161e8 1659 Double_t s = 1/yd/1.2;
1660 // Info("FixAxis", "for %s, s=1/%f=%f", stack->GetName(), yd, s);
1661
1662 h->SetXTitle("#font[152]{#eta}");
2f92ff0e 1663 h->SetYTitle(ytitle);
1664 TAxis* xa = h->GetXaxis();
1665 TAxis* ya = h->GetYaxis();
797161e8 1666
1667 // Int_t npixels = 20
1668 // Float_t dy = gPad->AbsPixeltoY(0) - gPad->AbsPixeltoY(npixels);
1669 // Float_t ts = dy/(gPad->GetY2() - gPad->GetY1());
1670
2f92ff0e 1671 if (xa) {
797161e8 1672 // xa->SetTitle(h->GetXTitle());
2f92ff0e 1673 // xa->SetTicks("+-");
1674 xa->SetTitleSize(s*xa->GetTitleSize());
1675 xa->SetLabelSize(s*xa->GetLabelSize());
1676 xa->SetTickLength(s*xa->GetTickLength());
797161e8 1677 // xa->SetTitleOffset(xa->GetTitleOffset()/s);
e308a636 1678
1679 if (stack != fResults) {
1680 TAxis* rxa = fResults->GetXaxis();
1681 xa->Set(rxa->GetNbins(), rxa->GetXmin(), rxa->GetXmax());
1682 }
2f92ff0e 1683 }
1684 if (ya) {
797161e8 1685 // ya->SetTitle(h->GetYTitle());
2f92ff0e 1686 ya->SetDecimals();
1687 // ya->SetTicks("+-");
1688 ya->SetNdivisions(ynDiv);
1689 ya->SetTitleSize(s*ya->GetTitleSize());
1690 ya->SetTitleOffset(ya->GetTitleOffset()/s);
1691 ya->SetLabelSize(s*ya->GetLabelSize());
b2e7f2d6 1692 }
b2e7f2d6 1693 }
797161e8 1694 //__________________________________________________________________
1695 /**
1696 * Merge two histograms into one
1697 *
1698 * @param cen Central part
1699 * @param fwd Forward part
1700 * @param xlow On return, lower eta bound
1701 * @param xhigh On return, upper eta bound
1702 *
1703 * @return Newly allocated histogram or null
1704 */
1705 TH1*
1706 Merge(const TH1* cen, const TH1* fwd, Double_t& xlow, Double_t& xhigh)
1707 {
1708 TH1* tmp = static_cast<TH1*>(fwd->Clone("tmp"));
1709 TString name(fwd->GetName());
1710 name.ReplaceAll("Forward", "Merged");
1711 tmp->SetName(name);
1712
1713 // tmp->SetMarkerStyle(28);
1714 // tmp->SetMarkerColor(kBlack);
1715 tmp->SetDirectory(0);
1716 xlow = 100;
1717 xhigh = -100;
1718 for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
1719 Double_t cc = cen->GetBinContent(i);
1720 Double_t cf = fwd->GetBinContent(i);
1721 Double_t ec = cen->GetBinError(i);
1722 Double_t ef = fwd->GetBinError(i);
1723 Double_t nc = cf;
1724 Double_t ne = ef;
1725 if (cc < 0.001 && cf < 0.01) continue;
1726 xlow = TMath::Min(tmp->GetXaxis()->GetBinLowEdge(i),xlow);
1727 xhigh = TMath::Max(tmp->GetXaxis()->GetBinUpEdge(i),xhigh);
1728 if (cc > 0.001) {
1729 nc = cc;
1730 ne = ec;
1731 if (cf > 0.001) {
1732 nc = (cf + cc) / 2;
1733 ne = TMath::Sqrt(ec*ec + ef*ef);
1734 }
1735 }
1736 tmp->SetBinContent(i, nc);
1737 tmp->SetBinError(i, ne);
1738 }
1739 return tmp;
1740 }
1741 //____________________________________________________________________
1742 /**
1743 * Fit @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$ to histogram data
1744 *
1745 * @param tmp Histogram
1746 * @param xlow Lower x bound
1747 * @param xhigh Upper x bound
1748 *
1749 * @return Fitted function
1750 */
1751 TF1*
1752 FitMerged(TH1* tmp, Double_t xlow, Double_t xhigh)
1753 {
1754 TF1* tmpf = new TF1("tmpf", "gaus", xlow, xhigh);
1755 tmp->Fit(tmpf, "NQ", "");
1756 tmp->SetDirectory(0);
1757
1758 TF1* fit = new TF1("f", myFunc, xlow, xhigh, 4);
1759 fit->SetParNames("a_{1}", "a_{2}", "#sigma_{1}", "#sigma_{2}");
1760 fit->SetParameters(tmpf->GetParameter(0),
1761 .2,
1762 tmpf->GetParameter(2),
1763 tmpf->GetParameter(2)/4);
1764 fit->SetParLimits(3, 0, 100);
1765 fit->SetParLimits(4, 0, 100);
1766 tmp->Fit(fit,"0W","");
1767
1768 delete tmpf;
1769 return fit;
1770 }
1771 //____________________________________________________________________
1772 /**
1773 * Make band of systematic errors
1774 *
1775 * @param tmp Histogram
1776 * @param fit Fit
1777 */
1778 void
1779 MakeSysError(TH1* tmp, TH1* cen, TH1* fwd, TF1* fit)
1780 {
1781 for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
1782 Double_t tc = tmp->GetBinContent(i);
1783 if (tc < 0.01) continue;
1784 Double_t fc = fwd->GetBinContent(i);
1785 Double_t cc = cen->GetBinContent(i);
1786 Double_t sysErr = fFwdSysErr;
1787 if (cc > .01 && fc > 0.01)
1788 sysErr = (fFwdSysErr+fCenSysErr) / 2;
1789 else if (cc > .01)
1790 sysErr = fCenSysErr;
1791 Double_t x = tmp->GetXaxis()->GetBinCenter(i);
1792 Double_t y = fit->Eval(x);
1793 tmp->SetBinContent(i, y);
1794 tmp->SetBinError(i,sysErr*y);
1795 }
1796 TString name(tmp->GetName());
1797 name.ReplaceAll("Merged", "SysError");
1798 tmp->SetName(name);
e42bc740 1799 tmp->SetMarkerColor(SYSERR_COLOR);
1800 tmp->SetLineColor(SYSERR_COLOR);
1801 tmp->SetFillColor(SYSERR_COLOR);
1802 tmp->SetFillStyle(SYSERR_STYLE);
797161e8 1803 tmp->SetMarkerStyle(0);
1804 tmp->SetLineWidth(0);
1805 }
1806 void CorrectForward(TH1* h) const
1807 {
1808 if (!fRemoveOuters) return;
1809
1810 for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1811 Double_t eta = h->GetBinCenter(i);
1812 if (TMath::Abs(eta) < 2.3) {
1813 h->SetBinContent(i, 0);
1814 h->SetBinError(i, 0);
1815 }
1816 }
1817 }
1818 void CorrectCentral(TH1* h) const
1819 {
1820 if (fClusterScale.IsNull()) return;
1821 TString t(h->GetTitle());
1822 Info("CorrectCentral", "Replacing Central with Tracklets in %s", t.Data());
1823 t.ReplaceAll("Central", "Tracklets");
1824 h->SetTitle(t);
1825
1826 TF1* cf = new TF1("clusterScale", fClusterScale);
1827 for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1828 Double_t eta = h->GetBinCenter(i);
1829 Double_t f = cf->Eval(eta);
1830 Double_t c = h->GetBinContent(i);
1831 if (f < .1) f = 1;
1832 h->SetBinContent(i, c / f);
1833 }
1834 delete cf;
1835 }
33ab4a11 1836 //____________________________________________________________________
1837 void Export(const char* basename)
1838 {
1839 TString bname(basename);
1840 bname.ReplaceAll(" ", "_");
1841 bname.ReplaceAll("-", "_");
1842 TString fname(Form("%s.C", bname.Data()));
1843
1844 std::ofstream out(fname.Data());
1845 if (!out) {
1846 Error("Export", "Failed to open output file %s", fname.Data());
1847 return;
1848 }
1849 out << "// Create by dNdetaDrawer\n"
1850 << "void " << bname << "(THStack* stack, TLegend* l, Int_t m)\n"
1851 << "{"
1852 << " Int_t ma[] = { 24, 25, 26, 32,\n"
1853 << " 20, 21, 22, 33,\n"
1854 << " 34, 30, 29, 0, \n"
1855 << " 23, 27 };\n"
1856 << " Int_t mm = ((m < 20 || m > 34) ? 0 : ma[m-20]);\n\n";
1857 TList* hists = fResults->GetHists();
1858 TIter next(hists);
1859 TH1* hist = 0;
1860 while ((hist = static_cast<TH1*>(next()))) {
1861 TString hname = hist->GetName();
1862 hname.Append(Form("_%04x", (gRandom->Integer(0xffff) & 0xffff)));
1863 hist->SetName(hname);
1864 hist->GetListOfFunctions()->Clear();
1865 hist->SavePrimitive(out, "nodraw");
1866 bool mirror = hname.Contains("mirror");
1867 bool syserr = hname.Contains("SysError");
1868 if (!syserr)
1869 out << " " << hname << "->SetMarkerStyle("
1870 << (mirror ? "mm" : "m") << ");\n";
1871 else
1872 out << " " << hname << "->SetMarkerStyle(1);\n";
1873 out << " stack->Add(" << hname
1874 << (syserr ? ",\"e5\"" : "") << ");\n\n";
1875 }
1876 UShort_t snn = fSNNString->GetUniqueID();
1877 // const char* sys = fSysString->GetTitle();
1878 TString eS;
1879 if (snn == 2750) snn = 2760;
1880 if (snn < 1000) eS = Form("%3dGeV", snn);
1881 else if (snn % 1000 == 0) eS = Form("%dTeV", snn/1000);
1882 else eS = Form("%4.2fTeV", float(snn)/1000);
1883 out << " if (l) {\n"
1884 << " TLegendEntry* e = l->AddEntry(\"\",\"" << eS << "\",\"pl\");\n"
1885 << " e->SetMarkerStyle(m);\n"
1886 << " e->SetMarkerColor(kBlack);\n"
1887 << " }\n"
1888 << "}\n" << std::endl;
1889 }
797161e8 1890
2f92ff0e 1891 /* @} */
1892
b2e7f2d6 1893
1894
1895 //__________________________________________________________________
33ab4a11 1896 /**
1897 * @{
1898 * @name Options
1899 */
e308a636 1900 Bool_t fShowRatios; // Show ratios
1901 Bool_t fShowLeftRight;// Show asymmetry
f709afcf 1902 Bool_t fShowRings; // Show rings too
33ab4a11 1903 Bool_t fExport; // Export results to file
e308a636 1904 Bool_t fCutEdges; // Whether to cut edges
33ab4a11 1905 Bool_t fRemoveOuters; // Whether to remove outers
1906 UShort_t fShowOthers; // Show other data
1907 /* @} */
1908 /**
1909 * @{
1910 * @name Settings
1911 */
1912 UShort_t fRebin; // Rebinning factor
1913 Double_t fFwdSysErr; // Systematic error in forward range
1914 Double_t fCenSysErr; // Systematic error in central range
e308a636 1915 TString fTitle; // Title on plot
33ab4a11 1916 TString fClusterScale; // Scaling of clusters to tracklets
1917 /* @} */
1918 /**
1919 * @{
1920 * @name Read (or set) information
1921 */
e308a636 1922 TNamed* fTrigString; // Trigger string (read, or set)
0be6c8cd 1923 TNamed* fNormString; // Normalisation string (read, or set)
e308a636 1924 TNamed* fSNNString; // Energy string (read, or set)
1925 TNamed* fSysString; // Collision system string (read or set)
1926 TAxis* fVtxAxis; // Vertex cuts (read or set)
1927 TAxis* fCentAxis; // Centrality axis
33ab4a11 1928 /* @} */
1929 /**
1930 * @{
1931 * @name Resulting plots
1932 */
e308a636 1933 THStack* fResults; // Stack of results
1934 THStack* fRatios; // Stack of ratios
1935 THStack* fLeftRight; // Left-right asymmetry
1936 TMultiGraph* fOthers; // Older data
1937 TH1* fTriggers; // Number of triggers
9ecab72f 1938 TH1* fTruth; // Pointer to truth
33ab4a11 1939 /* @} */
e308a636 1940 RangeParam* fRangeParam; // Parameter object for range zoom
b2e7f2d6 1941};
1942
797161e8 1943//____________________________________________________________________
1944/**
1945 * Function to calculate
1946 * @f[
1947 * g(x;A_1,A_2,\sigma_1,\sigma_2) =
1948 * A_1\left(\frac{1}{2\pi\sigma_1}e^{(x/\sigma_1)^2} -
1949 * A_2\frac{1}{2\pi\sigma_2}e^{(x/\sigma_2)^2}\right)
1950 * @f]
1951 *
1952 * @param xp Pointer to x array
1953 * @param pp Pointer to parameter array
1954 *
1955 * @return @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$
1956 */
1957Double_t myFunc(Double_t* xp, Double_t* pp)
1958{
1959 Double_t x = xp[0];
1960 Double_t a1 = pp[0];
1961 Double_t a2 = pp[1];
1962 Double_t s1 = pp[2];
1963 Double_t s2 = pp[3];
1964 return a1*(TMath::Gaus(x, 0, s1) - a2 * TMath::Gaus(x, 0, s2));
1965}
1966
b2e7f2d6 1967//=== Stuff for auto zooming =========================================
1968void UpdateRange(dNdetaDrawer::RangeParam* p)
1969{
1970 if (!p) {
1971 Warning("UpdateRange", "No parameters %p", p);
1972 return;
1973 }
1974 if (!p->fMasterAxis) {
1975 Warning("UpdateRange", "No master axis %p", p->fMasterAxis);
1976 return;
1977 }
1978 Int_t first = p->fMasterAxis->GetFirst();
1979 Int_t last = p->fMasterAxis->GetLast();
1980 Double_t x1 = p->fMasterAxis->GetBinCenter(first);
1981 Double_t x2 = p->fMasterAxis->GetBinCenter(last);
b2e7f2d6 1982
1983 if (p->fSlave1Axis) {
1984 Int_t i1 = p->fSlave1Axis->FindBin(x1);
1985 Int_t i2 = p->fSlave1Axis->FindBin(x2);
1986 p->fSlave1Axis->SetRange(i1, i2);
1987 p->fSlave1Pad->Modified();
1988 p->fSlave1Pad->Update();
1989 }
1990 if (p->fSlave2Axis) {
1991 Int_t i1 = p->fSlave2Axis->FindBin(x1);
1992 Int_t i2 = p->fSlave2Axis->FindBin(x2);
1993 p->fSlave2Axis->SetRange(i1, i2);
1994 p->fSlave2Pad->Modified();
1995 p->fSlave2Pad->Update();
1996 }
1997 TCanvas* c = gPad->GetCanvas();
1998 c->cd();
1999}
2000
2001//____________________________________________________________________
2002void RangeExec(dNdetaDrawer::RangeParam* p)
2003{
2004 // Event types:
2005 // 51: Mouse move
2006 // 53:
2007 // 1: Button down
2008 // 21: Mouse drag
2009 // 11: Mouse release
2010 // dNdetaDrawer::RangeParam* p =
2011 // reinterpret_cast<dNdetaDrawer::RangeParam*>(addr);
2012 Int_t event = gPad->GetEvent();
2013 TObject *select = gPad->GetSelected();
2014 if (event == 53) {
2015 UpdateRange(p);
2016 return;
2017 }
2018 if (event != 11 || !select || select != p->fMasterAxis) return;
2019 UpdateRange(p);
2020}
2021
33ab4a11 2022//=== Steering functions
2023//==============================================
2024void
2025Usage()
2026{
2027 Info("DrawdNdeta", "Usage: DrawdNdeta(FILE,TITLE,REBIN,OTHERS,FLAGS)\n\n"
2028 " const char* FILE File name to open (\"forward_root\")\n"
2029 " const char* TITLE Title to put on plot (\"\")\n"
2030 " UShort_t REBIN Rebinning factor (1)\n"
2031 " UShort_t OTHERS Other data to draw - more below (0x7)\n"
2032 " UShort_t FLAGS Visualisation flags - more below (0x7)\n\n"
2033 " OTHERS is a bit mask of\n\n"
2034 " 0x1 Show UA5 data (INEL,NSD, ppbar, 900GeV)\n"
2035 " 0x2 Show CMS data (NSD, pp)\n"
2036 " 0x4 Show published ALICE data (INEL,INEL>0,NSD, pp)\n"
2037 " 0x8 Show event genertor data\n\n"
2038 " FLAGS is a bit mask of\n\n"
2039 " 0x1 Show ratios of data to other data and possibly MC\n"
2040 " 0x2 Show left-right asymmetry\n"
2041 " 0x4 Show systematic error band\n"
2042 " 0x8 Show individual ring results (INEL only)\n"
2043 " 0x10 Cut edges when rebinning\n"
2044 " 0x20 Remove FMDxO points\n"
2045 " 0x40 Do not make our own canvas\n"
2046 );
2047}
2048
2049//____________________________________________________________________
ffca499d 2050/**
2051 * Draw @f$ dN/d\eta@f$
2052 *
2053 * @param filename File name
2054 * @param flags Flags
2055 * @param title Title
2056 * @param rebin Rebinning factor
2057 * @param cutEdges Whether to cut edges when rebinning
2058 * @param sNN (optional) Collision energy [GeV]
2059 * @param sys (optional) Collision system (1: pp, 2: PbPb)
2060 * @param trg (optional) Trigger (1: INEL, 2: INEL>0, 4: NSD)
2061 * @param vzMin Least @f$ v_z@f$
2062 * @param vzMax Largest @f$ v_z@f$
2063 *
2064 * @ingroup pwg2_forward_dndeta
2065 */
b2e7f2d6 2066void
2067DrawdNdeta(const char* filename="forward_dndeta.root",
bad9a3c1 2068 const char* title="",
b2e7f2d6 2069 UShort_t rebin=5,
797161e8 2070 UShort_t others=0x7,
2071 UShort_t flags=0x7,
fe52e455 2072 UShort_t sNN=0,
2073 UShort_t sys=0,
ffca499d 2074 UShort_t trg=0,
fe52e455 2075 Float_t vzMin=999,
2076 Float_t vzMax=-999)
b2e7f2d6 2077{
33ab4a11 2078 TString fname(filename);
2079 fname.ToLower();
2080 if (fname.CompareTo("help") == 0 ||
2081 fname.CompareTo("--help") == 0) {
2082 Usage();
2083 return;
2084 }
ffca499d 2085 dNdetaDrawer* pd = new dNdetaDrawer;
2086 dNdetaDrawer& d = *pd;
b2e7f2d6 2087 d.SetRebin(rebin);
b2e7f2d6 2088 d.SetTitle(title);
797161e8 2089 d.SetShowOthers(others);
2090 d.SetShowRatios(flags & 0x1);
2091 d.SetShowLeftRight(flags & 0x2);
2092 d.SetForwardSysError(flags & 0x4 ? 0.076 : 0);
2093 d.SetShowRings(flags & 0x8);
2094 d.SetCutEdges(flags & 0x10);
33ab4a11 2095 d.fRemoveOuters = (flags & 0x20);
2096 d.SetExport(flags & 0x40);
797161e8 2097 // d.fClusterScale = "1.06 -0.003*x +0.0119*x*x";
2f92ff0e 2098 // Do the below if your input data does not contain these settings
fe52e455 2099 if (sNN > 0) d.SetSNN(sNN); // Collision energy per nucleon pair (GeV)
2100 if (sys > 0) d.SetSys(sys); // Collision system (1:pp, 2:PbPB)
2101 if (trg > 0) d.SetTrigger(trg); // Collision trigger (1:INEL, 2:INEL>0, 4:NSD)
2102 if (vzMin < 999 && vzMax > -999)
2103 d.SetVertexRange(vzMin,vzMax); // Collision vertex range (cm)
b2e7f2d6 2104 d.Run(filename);
2105}
2f92ff0e 2106//____________________________________________________________________
2107//
2108// EOF
2109//
2110