]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/FORWARD/analysis2/DrawdNdeta.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Better layout of documentation pages
[u/mrichter/AliRoot.git] / PWGLF / 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 *
970b1a8a 11 * The script is <i>very</i> long - sigh - the joy of drawing
5bb5d1f6 12 * things nicely in ROOT
ffca499d 13 *
bd6f5206 14 * @ingroup pwglf_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>
bfab35d9 37#include <TParameter.h>
33ab4a11 38#include <fstream>
298a7264 39#include <iostream>
290052e7 40/** Systematic error color */
37079f20 41#define SYSERR_COLOR kBlue-10
290052e7 42/** Systematic error style */
37079f20 43#define SYSERR_STYLE 1001
b2e7f2d6 44
797161e8 45Double_t myFunc(Double_t* xp, Double_t* pp);
46
2f92ff0e 47/**
48 * Class to draw dN/deta results
49 *
bd6f5206 50 * @ingroup pwglf_forward_tasks_dndeta
51 * @ingroup pwglf_forward_dndeta
2f92ff0e 52 */
b2e7f2d6 53struct dNdetaDrawer
54{
2f92ff0e 55 /**
56 * POD of data for range zooming
57 */
b2e7f2d6 58 struct RangeParam
59 {
60 TAxis* fMasterAxis; // Master axis
61 TAxis* fSlave1Axis; // First slave axis
62 TVirtualPad* fSlave1Pad; // First slave pad
63 TAxis* fSlave2Axis; // Second slave axis
64 TVirtualPad* fSlave2Pad; // Second slave pad
65 };
66 //__________________________________________________________________
2f92ff0e 67 /**
68 * Constructor
69 *
70 */
b2e7f2d6 71 dNdetaDrawer()
33ab4a11 72 : // Options
73 fShowRatios(false), // Show ratios
74 fShowLeftRight(false), // Show asymmetry
75 fShowRings(false), // Show rings too
76 fExport(false), // Export data to script
77 fCutEdges(false), // Whether to cut edges
78 fRemoveOuters(false), // Whether to remove outers
79 fShowOthers(0), // Show other data
298a7264 80 fMirror(false),
81 fForceMB(false),
73b32206 82 fAddExec(false),
d43c6cd0 83 fSaveAsPDF(false),
84 fSaveAsPNG(false),
85 fSaveAsROOT(false),
86 fSaveAsScript(false),
33ab4a11 87 // Settings
88 fRebin(0), // Rebinning factor
89 fFwdSysErr(0.076), // Systematic error in forward range
90 fCenSysErr(0), // Systematic error in central range
91 fTitle(""), // Title on plot
d43c6cd0 92 fBase(""), // Optional base name of output files
33ab4a11 93 fClusterScale(""), // Scaling of clusters to tracklets
94 // Read (or set) information
95 fTrigString(0), // Trigger string (read, or set)
96 fNormString(0), // Normalisation string (read, or set)
97 fSNNString(0), // Energy string (read, or set)
98 fSysString(0), // Collision system string (read or set)
99 fVtxAxis(0), // Vertex cuts (read or set)
100 fCentAxis(0), // Centrality axis
bfab35d9 101 fTriggerEff(1), // Trigger efficency
d43c6cd0 102 fExtTriggerEff(false), // True if fTriggerEff was read
33ab4a11 103 // Resulting plots
104 fResults(0), // Stack of results
105 fRatios(0), // Stack of ratios
106 fLeftRight(0), // Left-right asymmetry
107 fOthers(0), // Older data
108 fTriggers(0), // Number of triggers
109 fTruth(0), // Pointer to truth
110 fRangeParam(0) // Parameter object for range zoom
b2e7f2d6 111 {
112 fRangeParam = new RangeParam;
113 fRangeParam->fMasterAxis = 0;
114 fRangeParam->fSlave1Axis = 0;
115 fRangeParam->fSlave1Pad = 0;
116 fRangeParam->fSlave2Axis = 0;
117 fRangeParam->fSlave2Pad = 0;
118 }
290052e7 119 /**
120 * Cpoy constructor
121 */
33ab4a11 122 dNdetaDrawer(const dNdetaDrawer&) {}
290052e7 123 /**
124 * Assignment operator
125 *
126 *
127 * @return Reference to this object
128 */
33ab4a11 129 dNdetaDrawer& operator=(const dNdetaDrawer&) { return *this; }
130
ffca499d 131 //__________________________________________________________________
290052e7 132 /**
133 * Destructor
134 */
ffca499d 135 virtual ~dNdetaDrawer()
136 {
137 if (fRatios && fRatios->GetHists()) fRatios->GetHists()->Delete();
138 if (fResults && fResults->GetHists()) fResults->GetHists()->Delete();
139
140 if (fTrigString) { delete fTrigString; fTrigString = 0; }
141 if (fSNNString) { delete fSNNString; fSNNString = 0; }
142 if (fSysString) { delete fSysString; fSysString = 0; }
143 if (fVtxAxis) { delete fVtxAxis; fVtxAxis = 0; }
144 if (fCentAxis) { delete fCentAxis; fCentAxis = 0; }
145 if (fResults) { delete fResults; fResults = 0; }
146 if (fRatios) { delete fRatios; fRatios = 0; }
147 if (fOthers) { delete fOthers; fOthers = 0; }
148 if (fTriggers) { delete fTriggers; fTriggers = 0; }
149 fRangeParam = 0;
150 }
151
2f92ff0e 152 //==================================================================
153 /**
154 * @{
155 * @name Set parameters
156 */
2e6b5465 157 void SetOld(Bool_t use=true) { fOld = use; }
2f92ff0e 158 /**
159 * Show other (UA5, CMS, ...) data
160 *
161 * @param x Whether to show or not
162 */
797161e8 163 void SetShowOthers(UShort_t x) { fShowOthers = x; }
2f92ff0e 164 //__________________________________________________________________
2f92ff0e 165 /**
166 * Whether to show ratios or not. If there's nothing to compare to,
167 * the ratio panel will be implicitly disabled
168 *
169 * @param x Whether to show or not
170 */
b2e7f2d6 171 void SetShowRatios(Bool_t x) { fShowRatios = x; }
2f92ff0e 172 //__________________________________________________________________
173 /**
174 *
175 * Whether to show the left/right asymmetry
176 *
177 * @param x To show or not
178 */
b2e7f2d6 179 void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; }
797161e8 180 //__________________________________________________________________
181 /**
182 * Whether to show rings
183 *
184 * @param x To show or not
185 */
f709afcf 186 void SetShowRings(Bool_t x) { fShowRings = x; }
2f92ff0e 187 //__________________________________________________________________
33ab4a11 188 /**
189 * Whether to export results to a script
190 *
191 * @param x Wheter to export results to a script
192 */
193 void SetExport(Bool_t x) { fExport = x; }
194 //__________________________________________________________________
2f92ff0e 195 /**
196 * Set the rebinning factor
197 *
198 * @param x Rebinning factor (must be a divisor in the number of bins)
199 */
b2e7f2d6 200 void SetRebin(UShort_t x) { fRebin = x; }
2f92ff0e 201 //__________________________________________________________________
202 /**
203 * Wheter to cut away the edges
204 *
205 * @param x Whether or not to cut away edges
206 */
b2e7f2d6 207 void SetCutEdges(Bool_t x) { fCutEdges = x; }
2f92ff0e 208 //__________________________________________________________________
209 /**
210 * Set the title of the plot
211 *
212 * @param x Title
213 */
b2e7f2d6 214 void SetTitle(TString x) { fTitle = x; }
797161e8 215 //__________________________________________________________________
d43c6cd0 216 /**
217 * Set the base name of the output files
218 *
219 * @param x Base name
220 */
221 void SetBase(TString x) { fBase = x; }
222 void SetSaveAsPNG(Bool_t x) { fSaveAsPNG = x; }
223 void SetSaveAsPDF(Bool_t x) { fSaveAsPDF = x; }
224 void SetSaveAsROOT(Bool_t x) { fSaveAsROOT = x; }
225 void SetSaveAsScript(Bool_t x) { fSaveAsScript = x; }
226 //__________________________________________________________________
797161e8 227 /**
228 * Set the systematic error in the forward region
229 *
230 * @param e Systematic error in the forward region
231 */
232 void SetForwardSysError(Double_t e=0) { fFwdSysErr = e; }
233 //__________________________________________________________________
234 /**
235 * Set the systematic error in the forward region
236 *
237 * @param e Systematic error in the forward region
238 */
239 void SetCentralSysError(Double_t e=0) { fCenSysErr = e; }
73b32206 240 /**
241 * Force the plot of minimum bias, even if centrality dependent data
242 * is present
243 *
244 * @param force if true, force minimum bias
245 */
298a7264 246 void SetForceMB(Bool_t force=true) { fForceMB = force; }
73b32206 247 /**
248 * Force the plot of minimum bias, even if centrality dependent data
249 * is present
250 *
33438b4c 251 * @param add if true, force minimum bias
73b32206 252 */
253 void SetAddExec(Bool_t add=true) { fAddExec = add; }
254 /**
255 * Mirror data to regions with no coverage (@f$-5.0<\eta<-3.5@f$)
256 *
257 * @param mirror If true, mirror data
258 */
298a7264 259 void SetMirror(Bool_t mirror=true) { fMirror = mirror; }
73b32206 260 /**
261 * Set the 'Final MC' correction file. This is needed if the
262 * secondary maps where produced using the old code
263 *
264 * @param file Filename
265 */
b615e211 266 void SetFinalMC(const TString& file) { fFinalMC = file; }
73b32206 267 /**
268 * Set the file that contains the empirical correction. This is
269 * needed when the secondary maps was generated with an in-accurate
270 * geometry, and when we're analysing data at nominal interaction
271 * points
272 *
273 * @param file Filename
274 */
b615e211 275 void SetEmpirical(const TString& file) { fEmpirical = file; }
2f92ff0e 276 /* @} */
277 //==================================================================
278 /**
279 * @{
280 * @name Override settings from input
281 */
282 /**
283 * Override setting from file
284 *
285 * @param sNN Center of mass energy per nucleon pair (GeV)
286 */
287 void SetSNN(UShort_t sNN)
288 {
cdb9e20f 289 fSNNString = new TNamed("sNN", Form("%04dGeV", sNN));
2f92ff0e 290 fSNNString->SetUniqueID(sNN);
291 }
292 //__________________________________________________________________
293 /**
294 * Set the collision system
295 * - 1: pp
296 * - 2: PbPb
297 *
298 * @param sys collision system
299 */
300 void SetSys(UShort_t sys)
301 {
fe52e455 302 fSysString = new TNamed("sys", (sys == 1 ? "pp" :
4d97b4bb 303 sys == 2 ? "PbPb" :
304 sys == 3 ? "pPb" :
305 "unknown"));
fe52e455 306 fSysString->SetUniqueID(sys);
2f92ff0e 307 }
308 //__________________________________________________________________
309 /**
310 * Set the vertex range in centimeters
311 *
312 * @param vzMin Min @f$ v_z@f$
313 * @param vzMax Max @f$ v_z@f$
314 */
315 void SetVertexRange(Double_t vzMin, Double_t vzMax)
316 {
317 fVtxAxis = new TAxis(10, vzMin, vzMax);
318 fVtxAxis->SetName("vtxAxis");
319 fVtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", vzMin, vzMax));
320 }
b2e7f2d6 321 //__________________________________________________________________
bfab35d9 322 /**
323 * Set the trigger mask (overrides what's in the file)
324 *
325 * @param trig Trigger mask (0x1: INEL, 0x2: INEL>0, 0x4: NSD)
326 */
2f92ff0e 327 void SetTrigger(UShort_t trig)
328 {
329 fTrigString = new TNamed("trigString", (trig & 0x1 ? "INEL" :
330 trig & 0x2 ? "INEL>0" :
331 trig & 0x4 ? "NSD" :
332 "unknown"));
333 fTrigString->SetUniqueID(trig);
334 }
bfab35d9 335 //__________________________________________________________________
336 /**
337 * Set the trigger efficiency - if set, then scale result histograms
338 * by this factor
339 *
340 * @param eff @f$\varepsilon_{T}@f$
341 */
d43c6cd0 342 void SetTriggerEfficiency(Float_t eff)
343 {
344 fTriggerEff = eff;
345 fExtTriggerEff = false;
346 }
2f92ff0e 347
348 //==================================================================
349 /**
350 * @{
351 * @name Main steering functions
352 */
353 /**
354 * Run the code to produce the final result.
355 *
356 * @param filename File containing the data
357 */
e308a636 358 void Run(const char* filename="forward_dndeta.root")
b2e7f2d6 359 {
e308a636 360 Double_t max = 0, rmax=0, amax=0;
b2e7f2d6 361
e308a636 362 gStyle->SetPalette(1);
e1f47419 363
e308a636 364 // --- Open input file -------------------------------------------
365 TFile* file = TFile::Open(filename, "READ");
366 if (!file) {
73b32206 367 Error("Run", "Cannot open %s", filename);
e308a636 368 return;
369 }
73b32206 370 Info("Run", "Drawing results from %s", file->GetName());
29befb1a 371
e308a636 372 // --- Get forward list ------------------------------------------
c8b1a7db 373 TList* forward = static_cast<TList*>(file->Get("ForwarddNdetaResults"));
e308a636 374 if (!forward) {
c8b1a7db 375 Error("Run", "Couldn't find list ForwarddNdetaResults");
e308a636 376 return;
377 }
c8b1a7db 378 TList* sums = static_cast<TList*>(file->Get("ForwarddNdetaSums"));
bfab35d9 379 if (!sums) {
c8b1a7db 380 Error("Run", "Couldn't find list ForwarddNdetaSums");
bfab35d9 381 return;
382 }
383 TParameter<bool>* p =
384 static_cast<TParameter<bool>*>(sums->FindObject("empirical"));
385 if (p && p->GetVal() && !fEmpirical.IsNull()) {
386 Warning("Run", "Empirical correction already applied");
387 fEmpirical = "__task__";
388 }
e308a636 389 // --- Get information on the run --------------------------------
390 FetchInformation(forward);
73b32206 391
392 // --- Print settings --------------------------------------------
393 Info("Run", "Settings for the drawer:\n"
394 " Show ratios: %5s\n"
395 " Show Left/right: %5s\n"
396 " Show rings: %5s\n"
397 " Export to file: %5s\n"
398 " Cut edges when rebinning: %5s\n"
399 " Remove outer rings: %5s\n"
400 " Mirror to un-covered regions: %5s\n"
401 " Force minimum bias: %5s\n"
402 " Show other results: 0x%03x\n"
403 " Rebinning factor: %5d\n"
404 " Forward systematic error: %5.1f%%\n"
405 " Central systematic error: %5.1f%%\n"
406 " Title on plot: %s\n"
407 " Scaling of clusters to tracklets: %s\n"
408 " Final MC correction file: %s\n"
d43c6cd0 409 " Empirical correction file: %s",
73b32206 410 (fShowRatios ? "yes" : "no"),
411 (fShowLeftRight ? "yes" : "no"),
412 (fShowRings ? "yes" : "no"),
413 (fExport ? "yes" : "no"),
414 (fCutEdges ? "yes" : "no"),
415 (fRemoveOuters ? "yes" : "no"),
416 (fMirror ? "yes" : "no"),
417 (fForceMB ? "yes" : "no"),
418 fShowOthers, fRebin, (100*fFwdSysErr), (100*fCenSysErr),
419 fTitle.Data(), fClusterScale.Data(), fFinalMC.Data(),
d43c6cd0 420 fEmpirical.Data());
73b32206 421
e308a636 422 // --- Set the macro pathand load other data script --------------
d43c6cd0 423 TString savPath(gROOT->GetMacroPath());
bd6f5206 424 gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2",
e308a636 425 gROOT->GetMacroPath()));
bfab35d9 426 // Always recompile
d43c6cd0 427 if (!gROOT->GetClass("RefData"))
428 gROOT->LoadMacro("OtherData.C++");
429 gROOT->SetMacroPath(savPath);
e308a636 430
431 // --- Get the central results -----------------------------------
c8b1a7db 432 TList* clusters = static_cast<TList*>(file->Get("CentraldNdetaResults"));
433 if (!clusters) Warning("Run", "Couldn't find list CentraldNdetaResults");
b2e7f2d6 434
9ecab72f 435 // --- Get the central results -----------------------------------
c8b1a7db 436 TList* mcTruth = static_cast<TList*>(file->Get("MCTruthdNdetaResults"));
437 if (!mcTruth) Warning("Run", "Couldn't find list MCTruthdNdetaResults");
9ecab72f 438
e308a636 439 // --- Make our containtes ---------------------------------------
440 fResults = new THStack("results", "Results");
441 fRatios = new THStack("ratios", "Ratios");
442 fLeftRight = new THStack("asymmetry", "Left-right asymmetry");
443 fOthers = new TMultiGraph();
b615e211 444
445 // --- Try to open the final MC file, and find relevant lists ----
446 TList* forwardMC = 0;
29befb1a 447 // TList* centralMC = 0;
b615e211 448 if (!fFinalMC.IsNull()) {
449 TFile* finalMC = TFile::Open(fFinalMC, "READ");
450 if (!finalMC) {
451 Warning("Run", "Failed to open file %s for final MC corrections",
452 fFinalMC.Data());
453 }
454 else {
c8b1a7db 455 forwardMC = static_cast<TList*>(finalMC->Get("ForwarddNdetaResults"));
b615e211 456 if (!forwardMC)
c8b1a7db 457 Warning("Run","Couldn't find list ForwarddNdetaResults for final MC");
29befb1a 458#if 0
c8b1a7db 459 centralMC = static_cast<TList*>(finalMC->Get("CentradNdetalResults"));
29befb1a 460 if (!centralMC)
c8b1a7db 461 Warning("Run","Couldn't find list CentraldNdetaResults for final MC");
29befb1a 462#endif
b615e211 463 }
464 }
465 if (!forwardMC) fFinalMC = "";
466
467 // --- Try to get the emperical correction -----------------------
468 TGraphErrors* empCorr = 0;
bfab35d9 469 if (!fEmpirical.IsNull() && !fEmpirical.EqualTo("__task__")) {
5934a3e3 470 if (gSystem->AccessPathName(fEmpirical.Data())) { // Not found here
471 fEmpirical =
bfab35d9 472 gSystem->ExpandPathName(Form("$ALICE_ROOT/PWGLF/FORWARD/"
473 "corrections/Empirical/%s",
474 fEmpirical.Data()));
5934a3e3 475 if (gSystem->AccessPathName(fEmpirical.Data())) { // Not found here
73b32206 476 Warning("Run", "Couldn't get empirical correction file");
5934a3e3 477 fEmpirical = "";
478 }
b615e211 479 }
5934a3e3 480 if (!fEmpirical.IsNull()) {
481 TFile* empirical = TFile::Open(fEmpirical, "READ");
482 if (!empirical) {
73b32206 483 Warning("Run", "couldn't open empirical correction file: %s",
5934a3e3 484 fEmpirical.Data());
485 fEmpirical = "";
486 }
487 const char* empPath = "fmdfull/average";
488 empCorr = static_cast<TGraphErrors*>(empirical->Get(empPath));
489 if (!empCorr) {
73b32206 490 Warning("Run", "Didn't find the graph %s in %s",
5934a3e3 491 empPath, fEmpirical.Data());
492 fEmpirical = "";
493 }
b615e211 494 }
495 }
bfab35d9 496 if (!empCorr && !fEmpirical.EqualTo("__task__")) fEmpirical = "";
b615e211 497
e308a636 498 // --- Loop over input data --------------------------------------
cd3b3781 499 TObjArray truths;
b615e211 500 FetchResults(mcTruth, 0, 0, "MCTruth", max, rmax, amax,truths);
501 TObjArray* fwdA = FetchResults(forward, forwardMC, empCorr, "Forward",
502 max, rmax, amax,truths);
503 TObjArray* cenA = FetchResults(clusters, 0, 0, "Central",
504 max, rmax, amax,truths);
b2e7f2d6 505
e308a636 506 // --- Get trigger information -----------------------------------
bfab35d9 507 // TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
e308a636 508 if (sums) {
2e6b5465 509 TList* all = 0;
510 if (fOld) all = sums;
511 else all = static_cast<TList*>(sums->FindObject("all"));
e308a636 512 if (all) {
513 fTriggers = FetchResult(all, "triggers");
514 if (!fTriggers) all->ls();
515 }
516 else {
517 Warning("Run", "List all not found in ForwardSums");
518 sums->ls();
519 }
520 }
521 else {
522 Warning("Run", "No ForwardSums directory found in %s", file->GetName());
523 file->ls();
524 }
525
526 // --- Check our stacks ------------------------------------------
527 if (!fResults->GetHists() ||
528 fResults->GetHists()->GetEntries() <= 0) {
529 Error("Run", "No histograms in result stack!");
530 return;
531 }
532 if (!fOthers->GetListOfGraphs() ||
533 fOthers->GetListOfGraphs()->GetEntries() <= 0) {
534 Warning("Run", "No other data found - disabling that");
797161e8 535 fShowOthers = 0;
e308a636 536 }
537 if (!fRatios->GetHists() ||
538 fRatios->GetHists()->GetEntries() <= 0) {
539 Warning("Run", "No ratio data found - disabling that");
540 // fRatios->ls();
541 fShowRatios = false;
542 }
543 if (!fLeftRight->GetHists() ||
544 fLeftRight->GetHists()->GetEntries() <= 0) {
545 Warning("Run", "No left/right data found - disabling that");
546 // fLeftRight->ls();
547 fShowLeftRight = false;
548 }
797161e8 549 if (fFwdSysErr > 0) {
550 if (fCenSysErr <= 0) fCenSysErr = fFwdSysErr;
551 for (Int_t i = 0; i < fwdA->GetEntriesFast(); i++) {
552 TH1* fwd = static_cast<TH1*>(fwdA->At(i));
553 TH1* cen = static_cast<TH1*>(cenA->At(i));
554 CorrectForward(fwd);
555 CorrectCentral(cen);
556 Double_t low, high;
557 TH1* tmp = Merge(cen, fwd, low, high);
558 TF1* f = FitMerged(tmp, low, high);
559 MakeSysError(tmp, cen, fwd, f);
560 delete f;
561 Info("", "Adding systematic error histogram %s",
562 tmp->GetName());
563 fResults->GetHists()->AddFirst(tmp, "e5");
298a7264 564
565 if (!fMirror) continue;
566
797161e8 567 TH1* tmp2 = Symmetrice(tmp);
568 tmp2->SetFillColor(tmp->GetFillColor());
569 tmp2->SetFillStyle(tmp->GetFillStyle());
570 tmp2->SetMarkerStyle(tmp->GetMarkerStyle());
571 tmp2->SetLineWidth(tmp->GetLineWidth());
572 fResults->GetHists()->AddFirst(tmp2, "e5");
573 fResults->Modified();
574 }
575 }
576 delete fwdA;
577 delete cenA;
e308a636 578
579 // --- Close the input file --------------------------------------
580 file->Close();
581
582
b2e7f2d6 583
e308a636 584 // --- Plot results ----------------------------------------------
585 Plot(max, rmax, amax);
b2e7f2d6 586 }
e1f47419 587
b2e7f2d6 588 //__________________________________________________________________
2f92ff0e 589 /**
e1f47419 590 * Fetch the information on the run from the results list
2f92ff0e 591 *
e1f47419 592 * @param results Results list
2f92ff0e 593 */
e1f47419 594 void FetchInformation(const TList* results)
b2e7f2d6 595 {
2f92ff0e 596 if (!fTrigString)
0be6c8cd 597 fTrigString = static_cast<TNamed*>(results->FindObject("trigger"));
598 if (!fNormString)
599 fNormString = static_cast<TNamed*>(results->FindObject("scheme"));
2f92ff0e 600 if (!fSNNString)
601 fSNNString = static_cast<TNamed*>(results->FindObject("sNN"));
602 if (!fSysString)
603 fSysString = static_cast<TNamed*>(results->FindObject("sys"));
604 if (!fVtxAxis)
605 fVtxAxis = static_cast<TAxis*>(results->FindObject("vtxAxis"));
e308a636 606 if (!fCentAxis)
607 fCentAxis = static_cast<TAxis*>(results->FindObject("centAxis"));
d43c6cd0 608 if (fTriggerEff <= 0 || (1-fTriggerEff)<1e-6) {
609 TParameter<double>* eff =
610 static_cast<TParameter<double>*>(results->FindObject("triggerEff"));
611 if (eff) {
612 fTriggerEff = eff->GetVal();
613 fExtTriggerEff = true;
614 }
615 if (fTriggerEff <= 0) SetTriggerEfficiency(1);
616 }
c25b5e1b 617
618 TNamed* options = static_cast<TAxis*>(results->FindObject("options"));
0be6c8cd 619 if (!fTrigString) fTrigString = new TNamed("trigger", "unknown");
620 if (!fNormString) fNormString = new TNamed("scheme", "unknown");
bad9a3c1 621 if (!fSNNString) fSNNString = new TNamed("sNN", "unknown");
622 if (!fSysString) fSysString = new TNamed("sys", "unknown");
623 if (!fVtxAxis) {
624 fVtxAxis = new TAxis(1,0,0);
625 fVtxAxis->SetName("vtxAxis");
626 fVtxAxis->SetTitle("v_{z} range unspecified");
627 }
b2e7f2d6 628
0be6c8cd 629 TString centTxt("none");
e308a636 630 if (fCentAxis) {
631 Int_t nCent = fCentAxis->GetNbins();
0be6c8cd 632 centTxt = Form("%d bins", nCent);
e308a636 633 for (Int_t i = 0; i <= nCent; i++)
0be6c8cd 634 centTxt.Append(Form("%c%d", i == 0 ? ' ' : '-',
ffca499d 635 int(fCentAxis->GetXbins()->At(i))));
e308a636 636 }
e1f47419 637 Info("FetchInformation",
fe52e455 638 "Initialized for\n"
29befb1a 639 " Trigger: %-30s (0x%x)\n"
d43c6cd0 640 " Efficiency: %-6.4f\n"
0be6c8cd 641 " sqrt(sNN): %-30s (%dGeV)\n"
642 " System: %-30s (%d)\n"
643 " Vz range: %-30s (%f,%f)\n"
644 " Normalization: %-30s (%d)\n"
c25b5e1b 645 " Centrality: %s\n"
646 " Options: %s",
fe52e455 647 fTrigString->GetTitle(), fTrigString->GetUniqueID(),
d43c6cd0 648 fTriggerEff,
fe52e455 649 fSNNString->GetTitle(), fSNNString->GetUniqueID(),
650 fSysString->GetTitle(), fSysString->GetUniqueID(),
e308a636 651 fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax(),
0be6c8cd 652 fNormString->GetTitle(), fNormString->GetUniqueID(),
c25b5e1b 653 centTxt.Data(), (options ? options->GetTitle() : "none"));
29befb1a 654 if (fSysString->GetUniqueID() == 3) {
655 Info("FetchResults", "Left/Right assymmetry, mirror, and systematic "
656 "errors explicitly disabled for pPb");
657 fShowLeftRight = false;
658 fMirror = false;
659 fFwdSysErr = 0;
660 fCenSysErr = 0;
661 }
e1f47419 662 }
e308a636 663 //__________________________________________________________________
664 TMultiGraph* FetchOthers(UShort_t centLow, UShort_t centHigh)
665 {
666 TMultiGraph* thisOther = 0;
797161e8 667 if (fShowOthers == 0) return 0;
e308a636 668
e308a636 669 UShort_t sys = (fSysString ? fSysString->GetUniqueID() : 0);
670 UShort_t trg = (fTrigString ? fTrigString->GetUniqueID() : 0);
671 UShort_t snn = (fSNNString ? fSNNString->GetUniqueID() : 0);
d43c6cd0 672 Long_t ret =
673 gROOT->ProcessLine(Form("RefData::GetData(%d,%d,%d,%d,%d,%d);",
674 sys,snn,trg,centLow,centHigh,fShowOthers));
797161e8 675 if (!ret) return 0;
676
677 thisOther = reinterpret_cast<TMultiGraph*>(ret);
e308a636 678 return thisOther;
679 }
e1f47419 680 //__________________________________________________________________
681 /**
e308a636 682 * Get the results from the top-level list
e1f47419 683 *
66cf95f2 684 * @param list List
685 * @param mcList List of histograms from MC
686 * @param empCorr Emperical correction if any
687 * @param name name
688 * @param max On return, maximum of data
689 * @param rmax On return, maximum of ratios
690 * @param amax On return, maximum of left-right comparisons
691 * @param truths List of MC truths to compare to.
290052e7 692 *
693 * @return Array of results
e1f47419 694 */
797161e8 695 TObjArray*
b615e211 696 FetchResults(const TList* list,
697 const TList* mcList,
698 TGraphErrors* empCorr,
699 const char* name,
700 Double_t& max,
701 Double_t& rmax,
702 Double_t& amax,
703 TObjArray& truths)
e1f47419 704 {
3d532276 705 if (!list) return 0;
298a7264 706 UShort_t n = HasCent() ? fCentAxis->GetNbins() : 0;
cd3b3781 707 // Info("FetchResults","got %d centrality bins", n);
e308a636 708 if (n == 0) {
b615e211 709 TH1* h = FetchOne(list, mcList, empCorr, name, "all",
710 FetchOthers(0,0), -1000, 0,
711 max, rmax, amax, fTruth);
712 if (!h) return 0;
797161e8 713 TObjArray* a = new TObjArray;
cd3b3781 714 // Info("FetchResults", "Adding %s to result stack", h->GetName());
797161e8 715 a->AddAt(h, 0);
716 return a;
b2e7f2d6 717 }
e1f47419 718
797161e8 719 TObjArray* a = new TObjArray;
cd3b3781 720 truths.Expand(n);
e308a636 721 for (UShort_t i = 0; i < n; i++) {
5bb5d1f6 722 UShort_t centLow = fCentAxis->GetBinLowEdge(i+1);
723 UShort_t centHigh = fCentAxis->GetBinUpEdge(i+1);
e308a636 724 TString lname = Form("cent%03d_%03d", centLow, centHigh);
5bb5d1f6 725 Int_t col = GetCentralityColor(i+1);
b615e211 726 TString centTxt = Form("%3d%%-%3d%% central", centLow, centHigh);
e308a636 727
cd3b3781 728 TH1* tt = static_cast<TH1*>(truths.At(i));
729 TH1* ot = tt;
b615e211 730 TH1* h = FetchOne(list, mcList, empCorr, name, lname,
731 FetchOthers(centLow,centHigh), col,
732 centTxt.Data(), max, rmax, amax, fTruth);
733 if (!h) continue;
cd3b3781 734 if (ot != tt) {
735 //Info("FetchResults", "old truth=%p new truth=%p (%s)", ot, tt, name);
736 truths.AddAt(tt, i);
737 }
738 // Info("FetchResults", "Adding %p to result stack", h);
797161e8 739 a->AddAt(h, i);
e1f47419 740 }
797161e8 741 return a;
e308a636 742 }
743 //__________________________________________________________________
290052e7 744 /**
745 * Get the color for a centrality bin
746 *
747 * @param bin Centrality bin
748 *
749 * @return Color
750 */
5bb5d1f6 751 Int_t GetCentralityColor(Int_t bin) const
752 {
3a8ad972 753 if (fCentAxis->GetNbins() < 6) {
754 switch (bin) {
755 case 1: return kRed+2;
756 case 2: return kGreen+2;
757 case 3: return kBlue+1;
758 case 4: return kCyan+1;
759 case 5: return kMagenta+1;
760 case 6: return kYellow+2;
761 }
762 }
5bb5d1f6 763 UShort_t centLow = fCentAxis->GetBinLowEdge(bin);
764 UShort_t centHigh = fCentAxis->GetBinUpEdge(bin);
765 Float_t fc = (centLow+double(centHigh-centLow)/2) / 100;
766 Int_t nCol = gStyle->GetNumberOfColors();
767 Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5));
768 Int_t col = gStyle->GetColorPalette(icol);
769 //Info("GetCentralityColor","%3d: %3d-%3d -> %3d",bin,centLow,centHigh,col);
770 return col;
771 }
772 //__________________________________________________________________
290052e7 773 /**
774 * Set attributed on a histogram
775 *
776 * @param h Histogram
777 * @param color Color
778 */
e308a636 779 void SetAttributes(TH1* h, Int_t color)
780 {
781 if (!h) return;
782 if (color < 0) return;
ffca499d 783 // h->SetLineColor(color);
e308a636 784 h->SetMarkerColor(color);
ffca499d 785 // h->SetFillColor(color);
e308a636 786 }
787 //__________________________________________________________________
290052e7 788 /**
789 * Set attributed on a graph
790 *
791 * @param g Graph
792 * @param color Color
793 */
e308a636 794 void SetAttributes(TGraph* g, Int_t color)
795 {
796 if (!g) return;
797 if (color < 0) return;
ffca499d 798 // g->SetLineColor(color);
e308a636 799 g->SetMarkerColor(color);
ffca499d 800 // g->SetFillColor(color);
b2e7f2d6 801 }
802 //__________________________________________________________________
290052e7 803 /**
804 * Modify the title
805 *
806 */
3d532276 807 void ModifyTitle(TNamed* h, const char* /*centTxt*/)
e308a636 808 {
3d532276 809 if (!h) return;
810
811 TString title(h->GetTitle());
812 title.ReplaceAll("ALICE ","");
813 if (title.Contains("Central"))
c8b1a7db 814 title.ReplaceAll("CentraldNdeta", "SPD clusters");
3d532276 815 if (title.Contains("Forward"))
c8b1a7db 816 title.ReplaceAll("ForwarddNdeta", "FMD");
3d532276 817 h->SetTitle(title);
818
819
5bb5d1f6 820 return;
821 // if (!centTxt || !h) return;
822 // h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
e308a636 823 }
b615e211 824 //__________________________________________________________________
825 TH1* FetchOne(const TList* list,
826 const TList* mcList,
827 TGraphErrors* empCorr,
828 const char* name,
829 const char* folderName,
830 TMultiGraph* others,
831 Int_t col,
832 const char* txt,
833 Double_t& max,
834 Double_t& rmax,
835 Double_t& amax,
836 TH1*& truth)
837 {
2e6b5465 838 TList* folder = 0;
839 if (fOld) folder = const_cast<TList*>(list);
840 else folder = static_cast<TList*>(list->FindObject(folderName));
b615e211 841 if (!folder) {
842 Error("FetchResults", "Couldn't find list '%s' in %s",
843 folderName, list->GetName());
844 return 0;
845 }
846 TList* mcFolder = 0;
847 if (mcList) {
848 mcFolder = static_cast<TList*>(mcList->FindObject(folderName));
849 if (!mcFolder)
850 Warning("FetchResults",
851 "Didn't find the list '%s' in %s for final MC correction",
852 folderName, mcList->GetName());
853 }
29befb1a 854 TObject* normCalc = folder->FindObject("normCalc");
855 if (normCalc) Info("FetchOne", "%s:\n%s", folderName, normCalc->GetTitle());
b615e211 856 TH1* h = FetchResults(folder, mcFolder, empCorr, name,
857 others, col, txt, max, rmax, amax, truth);
858 return h;
859 }
e308a636 860 //__________________________________________________________________
861 /**
862 * Fetch results for a particular centrality bin
863 *
864 * @param list List
66cf95f2 865 * @param mcList List of MC results
866 * @param empCorr Emperical correction if any
e308a636 867 * @param name Name
ffca499d 868 * @param thisOther Other graphs
869 * @param color Color
870 * @param centTxt Centrality text
e308a636 871 * @param max On return, data maximum
872 * @param rmax On return, ratio maximum
873 * @param amax On return, left-right maximum
66cf95f2 874 * @param truth MC truth to compare to or possibly update
290052e7 875 *
876 * @return Histogram of results
e308a636 877 */
b615e211 878 TH1* FetchResults(const TList* list,
879 const TList* mcList,
880 TGraphErrors* empCorr,
881 const char* name,
882 TMultiGraph* thisOther,
883 Int_t color,
884 const char* centTxt,
885 Double_t& max,
886 Double_t& rmax,
887 Double_t& amax,
888 TH1*& truth)
b2e7f2d6 889 {
cd3b3781 890
e1f47419 891 TH1* dndeta = FetchResult(list, Form("dndeta%s", name));
892 TH1* dndetaMC = FetchResult(list, Form("dndeta%sMC", name));
893 TH1* dndetaTruth = FetchResult(list, "dndetaTruth");
b615e211 894
895 if (mcList && FetchResult(mcList, "finalMCCorr"))
896 Warning("FetchResults", "dNdeta already corrected for final MC");
897 else
898 CorrectFinalMC(dndeta, mcList);
899
900 CorrectEmpirical(dndeta, empCorr);
bfab35d9 901 CorrectTriggerEff(dndeta);
902 CorrectTriggerEff(dndetaMC);
b615e211 903
e1f47419 904 TH1* dndetaSym = 0;
905 TH1* dndetaMCSym = 0;
e308a636 906 SetAttributes(dndeta, color);
298a7264 907 SetAttributes(dndetaMC, HasCent() ? color : color+2);
e308a636 908 SetAttributes(dndetaTruth,color);
909 SetAttributes(dndetaSym, color);
298a7264 910 SetAttributes(dndetaMCSym,HasCent() ? color : color+2);
911 if (dndetaMC && HasCent())
5bb5d1f6 912 dndetaMC->SetMarkerStyle(dndetaMC->GetMarkerStyle()+2);
298a7264 913 if (dndetaMCSym && HasCent())
5bb5d1f6 914 dndetaMCSym->SetMarkerStyle(dndetaMCSym->GetMarkerStyle()+2);
298a7264 915 if (dndetaTruth && HasCent()) {
5bb5d1f6 916 dndetaTruth->SetMarkerStyle(34);
9ecab72f 917 dndetaTruth->SetMarkerColor(kYellow-1);
918 }
298a7264 919 if (dndetaTruth) {
920 dndetaTruth->SetLineColor(kBlack);
cd3b3781 921 dndetaTruth->SetFillColor(kBlack);
922 dndetaTruth->SetFillStyle(3002);
923 // dndetaTruth->SetLineColor(kBlack);
298a7264 924 }
e308a636 925 ModifyTitle(dndeta, centTxt);
926 ModifyTitle(dndetaMC, centTxt);
927 ModifyTitle(dndetaTruth,centTxt);
928 ModifyTitle(dndetaSym, centTxt);
929 ModifyTitle(dndetaMCSym,centTxt);
b615e211 930
e1f47419 931
cd3b3781 932 max = TMath::Max(max, AddHistogram(fResults, dndetaTruth, "e5"));
e308a636 933 max = TMath::Max(max, AddHistogram(fResults, dndetaMC, dndetaMCSym));
934 max = TMath::Max(max, AddHistogram(fResults, dndeta, dndetaSym));
5bb5d1f6 935
9ecab72f 936 if (dndetaTruth) {
cd3b3781 937 truth = dndetaTruth;
5bb5d1f6 938 }
9ecab72f 939 else {
940 if (fShowRings) {
941 THStack* rings = static_cast<THStack*>(list->FindObject("dndetaRings"));
942 if (rings) {
943 TIter next(rings->GetHists());
944 TH1* hist = 0;
945 while ((hist = static_cast<TH1*>(next())))
946 max = TMath::Max(max, AddHistogram(fResults, hist));
ffca499d 947 }
b2e7f2d6 948 }
9ecab72f 949 // Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f",
950 // dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
951
952 if (fShowLeftRight) {
953 fLeftRight->Add(Asymmetry(dndeta, amax));
954 fLeftRight->Add(Asymmetry(dndetaMC, amax));
9ecab72f 955 }
956
957 if (thisOther) {
958 TIter next(thisOther->GetListOfGraphs());
959 TGraph* g = 0;
960 while ((g = static_cast<TGraph*>(next()))) {
961 fRatios->Add(Ratio(dndeta, g, rmax));
962 fRatios->Add(Ratio(dndetaSym, g, rmax));
963 SetAttributes(g, color);
964 ModifyTitle(g, centTxt);
965 if (!fOthers->GetListOfGraphs() ||
966 !fOthers->GetListOfGraphs()->FindObject(g->GetName())) {
967 max = TMath::Max(max,TMath::MaxElement(g->GetN(), g->GetY()));
968 fOthers->Add(g);
969 }
970 }
971 // fOthers->Add(thisOther);
972 }
b2e7f2d6 973 }
5bb5d1f6 974 if (dndetaMC) {
975 fRatios->Add(Ratio(dndeta, dndetaMC, rmax));
976 fRatios->Add(Ratio(dndetaSym, dndetaMCSym, rmax));
977 }
cd3b3781 978 if (truth) {
979 fRatios->Add(Ratio(dndeta, truth, rmax));
980 fRatios->Add(Ratio(dndetaSym, truth, rmax));
b2e7f2d6 981 }
797161e8 982 return dndeta;
b2e7f2d6 983 }
b615e211 984 //__________________________________________________________________
985 void CorrectFinalMC(TH1* dndeta, const TList* mcList)
986 {
987 if (!dndeta) return;
988 if (!mcList) return;
989
990 TH1* dndetaMC = FetchResult(mcList, dndeta->GetName());
991 TH1* dndetaTruth = FetchResult(mcList, "dndetaTruth");
992 if (!dndetaMC || !dndetaTruth) return;
993
994 TH1* corr = static_cast<TH1*>(dndetaMC->Clone("finalMCCorr"));
995 corr->Divide(dndetaTruth);
996
997 Info("CorrectFinalMC", "Correcting dN/deta with final MC correction");
998 dndeta->Divide(corr);
999 }
1000 //__________________________________________________________________
1001 void CorrectEmpirical(TH1* dndeta, const TGraphErrors* empCorr)
1002 {
1003 if (!dndeta) return;
1004 if (!empCorr) return;
1005
1006 Info("CorrectEmpirical", "Doing empirical correction of dN/deta");
1007 TAxis* xAxis = dndeta->GetXaxis();
1008 for (Int_t i = 1; i <= xAxis->GetNbins(); i++) {
1009 Double_t x = xAxis->GetBinCenter(i);
1010 Double_t y = dndeta->GetBinContent(i);
1011 Double_t c = empCorr->Eval(x);
1012 dndeta->SetBinContent(i, y / c);
1013 }
1014 }
bfab35d9 1015 //__________________________________________________________________
1016 void CorrectTriggerEff(TH1* dndeta)
1017 {
d43c6cd0 1018 if (fExtTriggerEff) return;
bfab35d9 1019 if (!dndeta) return;
1020 if (fTriggerEff <= 0 || fTriggerEff >= 1) return;
1021 dndeta->Scale(fTriggerEff);
1022 }
b2e7f2d6 1023 //__________________________________________________________________
2f92ff0e 1024 /**
1025 * Plot the results
2f92ff0e 1026 * @param max Max value
2f92ff0e 1027 * @param rmax Maximum diviation from 1 of ratios
2f92ff0e 1028 * @param amax Maximum diviation from 1 of asymmetries
1029 */
e308a636 1030 void Plot(Double_t max,
b2e7f2d6 1031 Double_t rmax,
b2e7f2d6 1032 Double_t amax)
1033 {
1034 gStyle->SetOptTitle(0);
d43c6cd0 1035 gStyle->SetTitleFont(kFont, "xyz");
1036 gStyle->SetLabelFont(kFont, "xyz");
b2e7f2d6 1037
e308a636 1038 Int_t h = 800;
1039 Int_t w = 800; // h / TMath::Sqrt(2);
1040 Double_t y1 = 0;
1041 Double_t y2 = 0;
1042 Double_t y3 = 0;
797161e8 1043 if (!fShowRatios) w *= 1.3;
e308a636 1044 else y1 = 0.3;
797161e8 1045 if (!fShowLeftRight) w *= 1.3;
e308a636 1046 else {
1047 Double_t y11 = y1;
1048 y1 = (y11 > 0.0001 ? 0.4 : 0.2);
ffca499d 1049 y2 = (y11 > 0.0001 ? 0.2 : 0.3);
e308a636 1050 }
b2e7f2d6 1051 TCanvas* c = new TCanvas("Results", "Results", w, h);
1052 c->SetFillColor(0);
1053 c->SetBorderSize(0);
1054 c->SetBorderMode(0);
1055
e308a636 1056 PlotResults(max, y1);
b2e7f2d6 1057 c->cd();
1058
e308a636 1059 PlotRatios(rmax, y2, y1);
b2e7f2d6 1060 c->cd( );
1061
e308a636 1062 PlotLeftRight(amax, y3, y2);
b2e7f2d6 1063 c->cd();
bad9a3c1 1064
1065
1066 Int_t vMin = fVtxAxis->GetXmin();
1067 Int_t vMax = fVtxAxis->GetXmax();
1068 TString trg(fTrigString->GetTitle());
e308a636 1069 Int_t nev = 0;
ffca499d 1070 if (fTriggers) nev = fTriggers->GetBinContent(1);
bad9a3c1 1071 trg = trg.Strip(TString::kBoth);
f709afcf 1072 trg.ReplaceAll(" ", "_");
1073 trg.ReplaceAll(">", "Gt");
29befb1a 1074 trg.ReplaceAll("&", "AND");
1075 trg.ReplaceAll("|", "OR");
d43c6cd0 1076 if (fBase.IsNull())
1077 fBase = "dndeta_<sys>_<snn>_<trig>_<ipmin><ipmax>cm_<nev>ev";
1078 fBase.ReplaceAll("<sys>", fSysString->GetTitle());
1079 fBase.ReplaceAll("<snn>", fSNNString->GetTitle());
1080 fBase.ReplaceAll("<trig>", trg.Data());
1081 fBase.ReplaceAll("<ipmin>", Form("%c%02d",vMin<0?'m':'p',TMath::Abs(vMin)));
1082 fBase.ReplaceAll("<ipmax>", Form("%c%02d",vMax<0?'m':'p',TMath::Abs(vMax)));
1083 fBase.ReplaceAll("<nev>", Form("%09d", nev));
1084 if (fSaveAsPNG) c->SaveAs(Form("%s.png", fBase.Data()));
1085 if (fSaveAsROOT) c->SaveAs(Form("%s.root", fBase.Data()));
1086 if (fSaveAsScript)c->SaveAs(Form("%s.C", fBase.Data()));
1087 if (fSaveAsPDF) c->SaveAs(Form("%s.pdf", fBase.Data()));
1088 TString exp(fBase);
1089 exp.ReplaceAll("dndeta", "export");
1090 Export(exp);
b2e7f2d6 1091 }
1092 //__________________________________________________________________
e308a636 1093 /**
1094 * Build main legend
1095 *
66cf95f2 1096 * @param stack Stack to include
1097 * @param mg (optional) Multi graph to include
1098 * @param x1 Lower X coordinate in the range [0,1]
1099 * @param y1 Lower Y coordinate in the range [0,1]
1100 * @param x2 Upper X coordinate in the range [0,1]
1101 * @param y2 Upper Y coordinate in the range [0,1]
1102 * @param forceCol If non-zero, force this many columns
e308a636 1103 */
ffca499d 1104 void BuildLegend(THStack* stack, TMultiGraph* mg,
cd3b3781 1105 Double_t x1, Double_t y1, Double_t x2, Double_t y2,
1106 Int_t forceCol=0)
e308a636 1107 {
1108 TLegend* l = new TLegend(x1,y1,x2,y2);
cd3b3781 1109 Int_t nCol = forceCol;
1110 if (nCol <= 0) nCol = HasCent() ? 1 : 2;
1111 l->SetNColumns(nCol);
e308a636 1112 l->SetFillColor(0);
1113 l->SetFillStyle(0);
1114 l->SetBorderSize(0);
d43c6cd0 1115 l->SetTextFont(kFont);
1116 l->SetTextColor(kAliceBlue);
e308a636 1117
797161e8 1118 // Loop over items in stack and get unique items, while ignoring
1119 // mirrored data and systematic error bands
ffca499d 1120 TIter next(stack->GetHists());
5bb5d1f6 1121 TH1* hist = 0;
1122 TObjArray unique;
1123 unique.SetOwner();
797161e8 1124 Bool_t sysErrSeen = false;
5bb5d1f6 1125 while ((hist = static_cast<TH1*>(next()))) {
797161e8 1126 TString t(hist->GetTitle());
1127 TString n(hist->GetName());
1128 n.ToLower();
1129 if (t.Contains("mirrored")) continue;
1130 if (n.Contains("syserror")) { sysErrSeen = true; continue; }
1131 if (unique.FindObject(t.Data())) continue;
298a7264 1132 TObjString* s1 = new TObjString(hist->GetTitle());
1133 s1->SetUniqueID(((hist->GetMarkerStyle() & 0xFFFF) << 16) |
1134 ((hist->GetMarkerColor() & 0xFFFF) << 0));
1135 unique.Add(s1);
5bb5d1f6 1136 // l->AddEntry(hist, hist->GetTitle(), "pl");
e308a636 1137 }
ffca499d 1138 if (mg) {
797161e8 1139 // If we have other stuff, scan for unique names
ffca499d 1140 TIter nexto(mg->GetListOfGraphs());
5bb5d1f6 1141 TGraph* g = 0;
1142 while ((g = static_cast<TGraph*>(nexto()))) {
1143 TString n(g->GetTitle());
ffca499d 1144 if (n.Contains("mirrored")) continue;
5bb5d1f6 1145 if (unique.FindObject(n.Data())) continue;
298a7264 1146 TObjString* s2 = new TObjString(n);
1147 s2->SetUniqueID(((g->GetMarkerStyle() & 0xFFFF) << 16) |
1148 ((g->GetMarkerColor() & 0xFFFF) << 0));
1149 unique.Add(s2);
5bb5d1f6 1150 // l->AddEntry(hist, hist->GetTitle(), "pl");
ffca499d 1151 }
e308a636 1152 }
797161e8 1153
1154 // Add legend entries for unique items only
5bb5d1f6 1155 TIter nextu(&unique);
1156 TObject* s = 0;
1157 Int_t i = 0;
1158 while ((s = nextu())) {
1159 TLegendEntry* dd = l->AddEntry(Form("data%2d", i++),
1160 s->GetName(), "lp");
1161 Int_t style = (s->GetUniqueID() >> 16) & 0xFFFF;
1162 Int_t color = (s->GetUniqueID() >> 0) & 0xFFFF;
1163 dd->SetLineColor(kBlack);
298a7264 1164 if (HasCent()) dd->SetMarkerColor(kBlack);
5bb5d1f6 1165 else dd->SetMarkerColor(color);
1166 dd->SetMarkerStyle(style);
1167 }
797161e8 1168 if (sysErrSeen) {
1169 // Add entry for systematic errors
1170 TLegendEntry* d0 = l->AddEntry("d0", Form("%4.1f%% Systematic error",
1171 100*fFwdSysErr), "f");
e42bc740 1172 d0->SetLineColor(SYSERR_COLOR);
1173 d0->SetMarkerColor(SYSERR_COLOR);
1174 d0->SetFillColor(SYSERR_COLOR);
1175 d0->SetFillStyle(SYSERR_STYLE);
797161e8 1176 d0->SetMarkerStyle(0);
1177 d0->SetLineWidth(0);
1178 i++;
1179 }
cd3b3781 1180 if (nCol == 2 && i % 2 == 1) {
797161e8 1181 // To make sure the 'data' and 'mirrored' entries are on a line
1182 // by themselves
1183 TLegendEntry* dd = l->AddEntry("dd", " ", "");
1184 dd->SetTextSize(0);
1185 dd->SetFillStyle(0);
1186 dd->SetFillColor(0);
1187 dd->SetLineWidth(0);
1188 dd->SetLineColor(0);
1189 dd->SetMarkerSize(0);
1190 }
298a7264 1191 if (fMirror) {
cd3b3781 1192 // Add entry for 'data'
1193 TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp");
1194 d1->SetLineColor(kBlack);
1195 d1->SetMarkerColor(kBlack);
1196 d1->SetMarkerStyle(20);
1197
1198 // Add entry for 'mirrored data'
298a7264 1199 TLegendEntry* d2 = l->AddEntry("d2", "Mirrored data", "lp");
1200 d2->SetLineColor(kBlack);
1201 d2->SetMarkerColor(kBlack);
1202 d2->SetMarkerStyle(24);
1203 }
e308a636 1204
1205 l->Draw();
1206 }
1207 //__________________________________________________________________
5bb5d1f6 1208 /**
1209 * Build centrality legend
1210 *
5bb5d1f6 1211 * @param x1 Lower X coordinate in the range [0,1]
1212 * @param y1 Lower Y coordinate in the range [0,1]
1213 * @param x2 Upper X coordinate in the range [0,1]
1214 * @param y2 Upper Y coordinate in the range [0,1]
1215 */
1216 void BuildCentLegend(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
1217 {
298a7264 1218 if (!HasCent()) return;
5bb5d1f6 1219
1220 TLegend* l = new TLegend(x1,y1,x2,y2);
1221 l->SetNColumns(1);
1222 l->SetFillColor(0);
1223 l->SetFillStyle(0);
1224 l->SetBorderSize(0);
d43c6cd0 1225 l->SetTextFont(kFont);
1226 l->SetTextColor(kAliceBlue);
5bb5d1f6 1227
1228 Int_t n = fCentAxis->GetNbins();
1229 for (Int_t i = 1; i <= n; i++) {
1230 Double_t low = fCentAxis->GetBinLowEdge(i);
1231 Double_t upp = fCentAxis->GetBinUpEdge(i);
1232 TLegendEntry* e = l->AddEntry(Form("dummy%02d", i),
1233 Form("%3d%% - %3d%%",
1234 int(low), int(upp)), "pl");
1235 e->SetMarkerColor(GetCentralityColor(i));
1236 }
1237 l->Draw();
1238 }
1239 //__________________________________________________________________
2f92ff0e 1240 /**
1241 * Plot the results
1242 *
2f92ff0e 1243 * @param max Maximum
1244 * @param yd Bottom position of pad
1245 */
e308a636 1246 void PlotResults(Double_t max, Double_t yd)
b2e7f2d6 1247 {
1248 // Make a sub-pad for the result itself
1249 TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0);
1250 p1->SetTopMargin(0.05);
1251 p1->SetBorderSize(0);
1252 p1->SetBorderMode(0);
1253 p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1);
d43c6cd0 1254 p1->SetRightMargin(kRightMargin);
797161e8 1255 if (fShowLeftRight || fShowRatios) p1->SetGridx();
b2e7f2d6 1256 p1->SetTicks(1,1);
1257 p1->SetNumber(1);
1258 p1->Draw();
1259 p1->cd();
33ab4a11 1260
73b32206 1261 Info("PlotResults", "Plotting results with max=%f", max);
e308a636 1262 fResults->SetMaximum(1.15*max);
cd3b3781 1263 fResults->SetMinimum(yd > 0.00001 ? -0.02*max : 0);
73b32206 1264 // fResults->SetMinimum(yd > 0.00001 ? -0.02*max : 0);
b2e7f2d6 1265
797161e8 1266 FixAxis(fResults, (1-yd)*(yd > .001 ? 1 : .9 / 1.2),
bfab35d9 1267 "#frac{1}{#it{N}}#kern[.1]{#frac{d#it{N}_{ch}}{d#it{#eta}}}");
b2e7f2d6 1268
1269 p1->Clear();
e308a636 1270 fResults->DrawClone("nostack e1");
b2e7f2d6 1271
e308a636 1272 fRangeParam->fSlave1Axis = fResults->GetXaxis();
b2e7f2d6 1273 fRangeParam->fSlave1Pad = p1;
1274
1275 // Draw other data
797161e8 1276 if (fShowOthers != 0) {
b2e7f2d6 1277 TGraphAsymmErrors* o = 0;
e308a636 1278 TIter nextG(fOthers->GetListOfGraphs());
b2e7f2d6 1279 while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
1280 o->DrawClone("same p");
1281 }
1282
1283 // Make a legend in the main result pad
5bb5d1f6 1284 BuildCentLegend(.12, 1-p1->GetTopMargin()-.01-.5,
1285 .35, 1-p1->GetTopMargin()-.01-.1);
298a7264 1286 Double_t x1 = (HasCent() ? .7 : .15);
1287 Double_t x2 = (HasCent() ? 1-p1->GetRightMargin()-.01: .90);
1288 Double_t y1 = (HasCent() ? .5: p1->GetBottomMargin()+.01);
1289 Double_t y2 = (HasCent() ? 1-p1->GetTopMargin()-.01-.15 : .35);
5bb5d1f6 1290
1291 BuildLegend(fResults, fOthers, x1, y1, x2, y2);
b2e7f2d6 1292
1293 // Put a title on top
ffca499d 1294 fTitle.ReplaceAll("@", " ");
d43c6cd0 1295 TLatex* tit = new TLatex(0.10, .99, fTitle.Data());
b2e7f2d6 1296 tit->SetNDC();
d43c6cd0 1297 tit->SetTextFont(kFont);
1298 tit->SetTextAlign(13);
1299 tit->SetTextSize(0.045);
1300 tit->SetTextColor(kAlicePurple);
b2e7f2d6 1301 tit->Draw();
1302
d43c6cd0 1303 // TColor::GetColor(41,73,156);
1304 Double_t x = .95;
bfab35d9 1305 Double_t y = .93;
b2e7f2d6 1306 // Put a nice label in the plot
1307 TString eS;
bad9a3c1 1308 UShort_t snn = fSNNString->GetUniqueID();
1309 const char* sys = fSysString->GetTitle();
5bb5d1f6 1310 if (snn == 2750) snn = 2760;
b2e7f2d6 1311 if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
1312 else eS = Form("%3dGeV", snn);
bfab35d9 1313 TLatex* tt = new TLatex(x, y, Form("%s #sqrt{s%s}=%s, %s",
1314 sys,
1315 (HasCent() ? "_{NN}" : ""),
1316 eS.Data(),
1317 HasCent() ? "by centrality" :
1318 fTrigString->GetTitle()));
d43c6cd0 1319 tt->SetTextColor(kAliceBlue);
b2e7f2d6 1320 tt->SetNDC();
d43c6cd0 1321 tt->SetTextFont(kFont);
b2e7f2d6 1322 tt->SetTextAlign(33);
1323 tt->Draw();
bfab35d9 1324 y -= tt->GetTextSize() + .01;
1325
b2e7f2d6 1326 // Put number of accepted events on the plot
e308a636 1327 Int_t nev = 0;
ffca499d 1328 if (fTriggers) nev = fTriggers->GetBinContent(1);
bfab35d9 1329 TLatex* et = new TLatex(x, y, Form("%d events", nev));
d43c6cd0 1330 et->SetTextColor(kAliceBlue);
b2e7f2d6 1331 et->SetNDC();
d43c6cd0 1332 et->SetTextFont(kFont);
b2e7f2d6 1333 et->SetTextAlign(33);
1334 et->Draw();
bfab35d9 1335 y -= et->GetTextSize() + .01;
b2e7f2d6 1336
1337 // Put number of accepted events on the plot
1338 if (fVtxAxis) {
bfab35d9 1339 TLatex* vt = new TLatex(x, y, fVtxAxis->GetTitle());
b2e7f2d6 1340 vt->SetNDC();
d43c6cd0 1341 vt->SetTextFont(kFont);
b2e7f2d6 1342 vt->SetTextAlign(33);
d43c6cd0 1343 vt->SetTextColor(kAliceBlue);
b2e7f2d6 1344 vt->Draw();
bfab35d9 1345 y -= vt->GetTextSize() + .01;
b2e7f2d6 1346 }
1347 // results->Draw("nostack e1 same");
1348
b615e211 1349 TString corrs;
1350 if (!fEmpirical.IsNull()) corrs.Append("Emperical");
1351 if (!fFinalMC.IsNull()) {
1352 if (!corrs.IsNull()) corrs.Append("+");
1353 corrs.Append("Final MC");
1354 }
1355
1356 if (!corrs.IsNull()) {
1357 corrs.Append(" correction");
1358 if (corrs.Index("+") != kNPOS) corrs.Append("s");
bfab35d9 1359 TLatex* em = new TLatex(x, y, corrs);
b615e211 1360 em->SetNDC();
d43c6cd0 1361 em->SetTextFont(kFont);
b615e211 1362 em->SetTextAlign(33);
d43c6cd0 1363 em->SetTextColor(kAliceBlue);
b615e211 1364 em->Draw();
bfab35d9 1365 y -= em->GetTextSize() + .01;
b615e211 1366 }
1367
bfab35d9 1368 if (fTriggerEff > 0 && fTriggerEff <= 1 && !HasCent()) {
1369 TLatex* ef = new TLatex(x, y, Form("#varepsilon_{%s} = %5.3f",
1370 fTrigString->GetTitle(),
1371 fTriggerEff));
1372 ef->SetNDC();
d43c6cd0 1373 ef->SetTextFont(kFont);
bfab35d9 1374 ef->SetTextAlign(33);
d43c6cd0 1375 ef->SetTextColor(kAliceBlue);
bfab35d9 1376 ef->Draw();
1377 y -= ef->GetTextSize() + .01;
1378 }
1379
e308a636 1380 fRangeParam->fSlave1Axis = FindXAxis(p1, fResults->GetName());
b2e7f2d6 1381 fRangeParam->fSlave1Pad = p1;
1382
1383
1384 // Mark the plot as preliminary
f709afcf 1385 TLatex* pt = new TLatex(.12, .93, "Work in progress");
b2e7f2d6 1386 pt->SetNDC();
d43c6cd0 1387 pt->SetTextFont(62);
f709afcf 1388 // pt->SetTextSize();
d43c6cd0 1389 pt->SetTextColor(kAliceRed);
b2e7f2d6 1390 pt->SetTextAlign(13);
1391 pt->Draw();
d43c6cd0 1392 TDatime now;
1393 TLatex* dt = new TLatex(.12, .88, now.AsSQLString());
1394 dt->SetNDC();
1395 dt->SetTextFont(42);
1396 dt->SetTextSize(0.04);
1397 dt->SetTextColor(kAliceBlue); // kAliceRed);
1398 dt->SetTextAlign(13);
1399 dt->Draw();
1400
1401 TString savPath(gROOT->GetMacroPath());
1402 gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2/scripts",
1403 gROOT->GetMacroPath()));
1404 // Always recompile
1405 if (!gROOT->GetClass("AliceLogo"))
1406 gROOT->LoadMacro("AliceLogo.C++");
1407 gROOT->SetMacroPath(savPath);
b2e7f2d6 1408
d43c6cd0 1409 if (gROOT->GetClass("AliceLogo")) {
1410 p1->cd();
1411 p1->Range(0,0,1,1);
1412 gROOT->ProcessLine("AliceLogo* al = new AliceLogo();");
1413 gROOT->ProcessLine("al->Draw(0,.88,.5,.2, 0, 0);");
b2e7f2d6 1414 }
1415 p1->cd();
1416 }
1417 //__________________________________________________________________
2f92ff0e 1418 /**
1419 * Plot the ratios
1420 *
2f92ff0e 1421 * @param max Maximum diviation from 1
1422 * @param y1 Lower y coordinate of pad
1423 * @param y2 Upper y coordinate of pad
1424 */
e308a636 1425 void PlotRatios(Double_t max, Double_t y1, Double_t y2)
b2e7f2d6 1426 {
797161e8 1427 if (fShowRatios == 0) return;
e308a636 1428
b2e7f2d6 1429 bool isBottom = (y1 < 0.0001);
1430 Double_t yd = y2 - y1;
1431 // Make a sub-pad for the result itself
1432 TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0);
1433 p2->SetTopMargin(0.001);
d43c6cd0 1434 p2->SetRightMargin(kRightMargin);
b2e7f2d6 1435 p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
1436 p2->SetGridx();
1437 p2->SetTicks(1,1);
1438 p2->SetNumber(2);
1439 p2->Draw();
1440 p2->cd();
1441
1442 // Fix up axis
797161e8 1443 FixAxis(fRatios, yd, "Ratios", 7);
b2e7f2d6 1444
e308a636 1445 fRatios->SetMaximum(1+TMath::Max(.22,1.05*max));
1446 fRatios->SetMinimum(1-TMath::Max(.32,1.05*max));
b2e7f2d6 1447 p2->Clear();
e308a636 1448 fRatios->DrawClone("nostack e1");
b2e7f2d6 1449
1450
1451 // Make a legend
ffca499d 1452 BuildLegend(fRatios, 0, .15,p2->GetBottomMargin()+.01,.9,
cd3b3781 1453 isBottom ? .6 : .4, 2);
ffca499d 1454#if 0
b2e7f2d6 1455 TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,
1456 isBottom ? .6 : .4);
1457 l2->SetNColumns(2);
1458 l2->SetFillColor(0);
1459 l2->SetFillStyle(0);
1460 l2->SetBorderSize(0);
d43c6cd0 1461 l2->SetTextFont(kFont);
ffca499d 1462#endif
b2e7f2d6 1463 // Make a nice band from 0.9 to 1.1
1464 TGraphErrors* band = new TGraphErrors(2);
e1f47419 1465 band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
1466 band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
b2e7f2d6 1467 band->SetPointError(0, 0, .1);
1468 band->SetPointError(1, 0, .1);
1469 band->SetFillColor(kYellow+2);
1470 band->SetFillStyle(3002);
1471 band->SetLineStyle(2);
1472 band->SetLineWidth(1);
1473 band->Draw("3 same");
1474 band->DrawClone("X L same");
1475
1476 // Replot the ratios on top
e308a636 1477 fRatios->DrawClone("nostack e1 same");
b2e7f2d6 1478
73b32206 1479 if (fAddExec) {
1480 if (isBottom) {
1481 fRangeParam->fMasterAxis = FindXAxis(p2, fRatios->GetName());
1482 p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
1483 fRangeParam));
1484 }
1485 else {
1486 fRangeParam->fSlave2Axis = FindXAxis(p2, fRatios->GetName());
1487 fRangeParam->fSlave2Pad = p2;
1488 }
b2e7f2d6 1489 }
1490 }
1491 //__________________________________________________________________
2f92ff0e 1492 /**
1493 * Plot the asymmetries
1494 *
2f92ff0e 1495 * @param max Maximum diviation from 1
1496 * @param y1 Lower y coordinate of pad
1497 * @param y2 Upper y coordinate of pad
1498 */
e308a636 1499 void PlotLeftRight(Double_t max, Double_t y1, Double_t y2)
b2e7f2d6 1500 {
e308a636 1501 if (!fShowLeftRight) return;
1502
b2e7f2d6 1503 bool isBottom = (y1 < 0.0001);
1504 Double_t yd = y2 - y1;
1505 // Make a sub-pad for the result itself
1506 TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0);
1507 p3->SetTopMargin(0.001);
d43c6cd0 1508 p3->SetRightMargin(kRightMargin);
b2e7f2d6 1509 p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
1510 p3->SetGridx();
1511 p3->SetTicks(1,1);
1512 p3->SetNumber(2);
1513 p3->Draw();
1514 p3->cd();
1515
b2e7f2d6 1516 // Fix up axis
797161e8 1517 FixAxis(fLeftRight, yd, "Right/Left", 4);
b2e7f2d6 1518
e308a636 1519 fLeftRight->SetMaximum(1+TMath::Max(.12,1.05*max));
1520 fLeftRight->SetMinimum(1-TMath::Max(.15,1.05*max));
b2e7f2d6 1521 p3->Clear();
e308a636 1522 fLeftRight->DrawClone("nostack e1");
b2e7f2d6 1523
1524
1525 // Make a legend
298a7264 1526 Double_t xx1 = (HasCent() ? .7 : .15);
1527 Double_t xx2 = (HasCent() ? 1-p3->GetRightMargin()-.01 : .90);
5bb5d1f6 1528 Double_t yy1 = p3->GetBottomMargin()+.01;
298a7264 1529 Double_t yy2 = (HasCent() ? 1-p3->GetTopMargin()-.01-.15 : .5);
5bb5d1f6 1530 BuildLegend(fLeftRight, 0, xx1, yy1, xx2, yy2);
1531 // TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
1532 // l2->SetNColumns(2);
1533 // l2->SetFillColor(0);
1534 // l2->SetFillStyle(0);
1535 // l2->SetBorderSize(0);
1536 // l2->SetTextFont(132);
ffca499d 1537
b2e7f2d6 1538 // Make a nice band from 0.9 to 1.1
1539 TGraphErrors* band = new TGraphErrors(2);
e1f47419 1540 band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
1541 band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
b2e7f2d6 1542 band->SetPointError(0, 0, .05);
1543 band->SetPointError(1, 0, .05);
1544 band->SetFillColor(kYellow+2);
1545 band->SetFillStyle(3002);
1546 band->SetLineStyle(2);
1547 band->SetLineWidth(1);
1548 band->Draw("3 same");
1549 band->DrawClone("X L same");
1550
e308a636 1551 fLeftRight->DrawClone("nostack e1 same");
b2e7f2d6 1552 if (isBottom) {
e308a636 1553 fRangeParam->fMasterAxis = FindXAxis(p3, fLeftRight->GetName());
b2e7f2d6 1554 p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
1555 fRangeParam));
1556 }
1557 else {
e308a636 1558 fRangeParam->fSlave2Axis = FindXAxis(p3, fLeftRight->GetName());
b2e7f2d6 1559 fRangeParam->fSlave2Pad = p3;
1560 }
1561 }
2f92ff0e 1562 /** @} */
1563 //==================================================================
1564 /**
1565 * @{
1566 * @name Data utility functions
1567 */
1568 /**
1569 * Get a result from the passed list
1570 *
1571 * @param list List to search
1572 * @param name Object name to search for
1573 *
290052e7 1574 * @return Histogram
2f92ff0e 1575 */
e1f47419 1576 TH1* FetchResult(const TList* list, const char* name) const
2f92ff0e 1577 {
1578 if (!list) return 0;
1579
e308a636 1580 TH1* ret = static_cast<TH1*>(list->FindObject(name));
5bb5d1f6 1581#if 0
e1f47419 1582 if (!ret) {
1583 // all->ls();
1584 Warning("GetResult", "Histogram %s not found", name);
2f92ff0e 1585 }
5bb5d1f6 1586#endif
b2e7f2d6 1587
e1f47419 1588 return ret;
2f92ff0e 1589 }
b2e7f2d6 1590 //__________________________________________________________________
2f92ff0e 1591 /**
1592 * Add a histogram to the stack after possibly rebinning it
1593 *
1594 * @param stack Stack to add to
1595 * @param hist histogram
1596 * @param option Draw options
1597 *
1598 * @return Maximum of histogram
1599 */
e308a636 1600 Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option="") const
b2e7f2d6 1601 {
1602 // Check if we have input
1603 if (!hist) return 0;
1604
1605 // Rebin if needed
1606 Rebin(hist);
1607
1608 stack->Add(hist, option);
1609 return hist->GetMaximum();
1610 }
1611 //__________________________________________________________________
2f92ff0e 1612 /**
1613 * Add a histogram to the stack after possibly rebinning it
1614 *
1615 * @param stack Stack to add to
1616 * @param hist histogram
1617 * @param option Draw options
1618 * @param sym On return, the data symmetriced (added to stack)
1619 *
1620 * @return Maximum of histogram
1621 */
e308a636 1622 Double_t AddHistogram(THStack* stack, TH1* hist, TH1*& sym,
1623 Option_t* option="") const
b2e7f2d6 1624 {
1625 // Check if we have input
1626 if (!hist) return 0;
1627
1628 // Rebin if needed
1629 Rebin(hist);
1630 stack->Add(hist, option);
1631
1632 // Now symmetrice the histogram
298a7264 1633 if (fMirror) {
1634 sym = Symmetrice(hist);
1635 stack->Add(sym, option);
1636 }
b2e7f2d6 1637
1638 return hist->GetMaximum();
1639 }
1640
1641 //__________________________________________________________________
1642 /**
1643 * Rebin a histogram
1644 *
1645 * @param h Histogram to rebin
b2e7f2d6 1646 */
1647 virtual void Rebin(TH1* h) const
1648 {
1649 if (fRebin <= 1) return;
1650
1651 Int_t nBins = h->GetNbinsX();
1652 if(nBins % fRebin != 0) {
1653 Warning("Rebin", "Rebin factor %d is not a devisor of current number "
1654 "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
1655 return;
1656 }
1657
1658 // Make a copy
1659 TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
1660 tmp->Rebin(fRebin);
1661 tmp->SetDirectory(0);
e28f5fc5 1662 tmp->Reset();
b2e7f2d6 1663 // The new number of bins
1664 Int_t nBinsNew = nBins / fRebin;
1665 for(Int_t i = 1;i<= nBinsNew; i++) {
1666 Double_t content = 0;
1667 Double_t sumw = 0;
1668 Double_t wsum = 0;
1669 Int_t nbins = 0;
1670 for(Int_t j = 1; j<=fRebin;j++) {
1671 Int_t bin = (i-1)*fRebin + j;
1672 Double_t c = h->GetBinContent(bin);
1673
1674 if (c <= 0) continue;
1675
1676 if (fCutEdges) {
1677 if (h->GetBinContent(bin+1)<=0 ||
73b32206 1678 h->GetBinContent(bin-1)<=0) {
b2e7f2d6 1679 Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
1680 bin, c, h->GetName(),
1681 bin+1, h->GetBinContent(bin+1),
1682 bin-1, h->GetBinContent(bin-1));
1683 continue;
1684 }
1685 }
1686 Double_t e = h->GetBinError(bin);
1687 Double_t w = 1 / (e*e); // 1/c/c
1688 content += c;
1689 sumw += w;
1690 wsum += w * c;
1691 nbins++;
1692 }
1693
e28f5fc5 1694 if(content > 0 && nbins > 1 ) {
b2e7f2d6 1695 tmp->SetBinContent(i, wsum / sumw);
1696 tmp->SetBinError(i,1./TMath::Sqrt(sumw));
1697 }
1698 }
1699
1700 // Finally, rebin the histogram, and set new content
1701 h->Rebin(fRebin);
1702 h->Reset();
1703 for(Int_t i = 1; i<= nBinsNew; i++) {
1704 h->SetBinContent(i,tmp->GetBinContent(i));
1705 h->SetBinError(i, tmp->GetBinError(i));
1706 }
1707
1708 delete tmp;
1709 }
1710 //__________________________________________________________________
1711 /**
1712 * Make an extension of @a h to make it symmetric about 0
1713 *
1714 * @param h Histogram to symmertrice
1715 *
1716 * @return Symmetric extension of @a h
1717 */
1718 TH1* Symmetrice(const TH1* h) const
1719 {
1720 Int_t nBins = h->GetNbinsX();
1721 TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
1722 s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
1723 s->Reset();
1724 s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
1725 s->SetMarkerColor(h->GetMarkerColor());
1726 s->SetMarkerSize(h->GetMarkerSize());
1727 s->SetMarkerStyle(h->GetMarkerStyle()+4);
1728 s->SetFillColor(h->GetFillColor());
1729 s->SetFillStyle(h->GetFillStyle());
1730 s->SetDirectory(0);
1731
1732 // Find the first and last bin with data
1733 Int_t first = nBins+1;
1734 Int_t last = 0;
1735 for (Int_t i = 1; i <= nBins; i++) {
1736 if (h->GetBinContent(i) <= 0) continue;
1737 first = TMath::Min(first, i);
1738 last = TMath::Max(last, i);
1739 }
1740
1741 Double_t xfirst = h->GetBinCenter(first-1);
1742 Int_t f1 = h->GetXaxis()->FindBin(-xfirst);
1743 Int_t l2 = s->GetXaxis()->FindBin(xfirst);
1744 for (Int_t i = f1, j=l2; i <= last; i++,j--) {
1745 s->SetBinContent(j, h->GetBinContent(i));
1746 s->SetBinError(j, h->GetBinError(i));
1747 }
1748 // Fill in overlap bin
1749 s->SetBinContent(l2+1, h->GetBinContent(first));
1750 s->SetBinError(l2+1, h->GetBinError(first));
1751 return s;
1752 }
1753 //__________________________________________________________________
2f92ff0e 1754 /**
1755 * Calculate the left-right asymmetry of input histogram
1756 *
1757 * @param h Input histogram
1758 * @param max On return, the maximum distance from 1 of the histogram
1759 *
1760 * @return Asymmetry
1761 */
1762 TH1* Asymmetry(TH1* h, Double_t& max)
1763 {
1764 if (!h) return 0;
1765
1766 TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName())));
1767 // Int_t oBins = h->GetNbinsX();
1768 // Double_t high = h->GetXaxis()->GetXmax();
1769 // Double_t low = h->GetXaxis()->GetXmin();
1770 // Double_t dBin = (high - low) / oBins;
1771 // Int_t tBins = Int_t(2*high/dBin+.5);
1772 // ret->SetBins(tBins, -high, high);
e308a636 1773 ret->SetDirectory(0);
2f92ff0e 1774 ret->Reset();
1775 ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
1776 ret->SetYTitle("Right/Left");
1777 Int_t nBins = h->GetNbinsX();
1778 for (Int_t i = 1; i <= nBins; i++) {
1779 Double_t x = h->GetBinCenter(i);
1780 if (x > 0) break;
1781
1782 Double_t c1 = h->GetBinContent(i);
1783 Double_t e1 = h->GetBinError(i);
1784 if (c1 <= 0) continue;
1785
1786 Int_t j = h->FindBin(-x);
1787 if (j <= 0 || j > nBins) continue;
1788
1789 Double_t c2 = h->GetBinContent(j);
1790 Double_t e2 = h->GetBinError(j);
1791
1792 Double_t c12 = c1*c1;
1793 Double_t e = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12));
1794
1795 Int_t k = ret->FindBin(x);
1796 ret->SetBinContent(k, c2/c1);
1797 ret->SetBinError(k, e);
1798 }
1799 max = TMath::Max(max, RatioMax(ret));
1800
1801 return ret;
1802 }
1803 //__________________________________________________________________
1804 /**
1805 * Transform a graph into a histogram
1806 *
290052e7 1807 * @param g Graph
2f92ff0e 1808 *
290052e7 1809 * @return Histogram
2f92ff0e 1810 */
1811 TH1* Graph2Hist(const TGraphAsymmErrors* g) const
1812 {
1813 Int_t nBins = g->GetN();
1814 TArrayF bins(nBins+1);
1815 Double_t dx = 0;
1816 for (Int_t i = 0; i < nBins; i++) {
1817 Double_t x = g->GetX()[i];
1818 Double_t exl = g->GetEXlow()[i];
1819 Double_t exh = g->GetEXhigh()[i];
1820 bins.fArray[i] = x-exl;
1821 bins.fArray[i+1] = x+exh;
1822 Double_t dxi = exh+exl;
1823 if (i == 0) dx = dxi;
1824 else if (dxi != dx) dx = 0;
1825 }
1826 TString name(g->GetName());
1827 TString title(g->GetTitle());
1828 TH1D* h = 0;
1829 if (dx != 0) {
1830 h = new TH1D(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
1831 }
1832 else {
1833 h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
1834 }
1835 h->SetMarkerStyle(g->GetMarkerStyle());
1836 h->SetMarkerColor(g->GetMarkerColor());
1837 h->SetMarkerSize(g->GetMarkerSize());
1838
1839 return h;
1840 }
1841 /* @} */
1842 //==================================================================
1843 /**
1844 * @{
1845 * @name Ratio utility functions
1846 */
1847 /**
1848 * Get the maximum diviation from 1 in the passed ratio
1849 *
1850 * @param h Ratio histogram
1851 *
1852 * @return Max diviation from 1
1853 */
b2e7f2d6 1854 Double_t RatioMax(TH1* h) const
1855 {
1856 Int_t nBins = h->GetNbinsX();
1857 Double_t ret = 0;
1858 for (Int_t i = 1; i <= nBins; i++) {
1859 Double_t c = h->GetBinContent(i);
1860 if (c == 0) continue;
1861 Double_t e = h->GetBinError(i);
1862 Double_t d = TMath::Abs(1-c-e);
1863 ret = TMath::Max(d, ret);
1864 }
1865 return ret;
1866 }
1867 //__________________________________________________________________
e1f47419 1868 /**
1869 * Make the ratio of h1 to h2
1870 *
ffca499d 1871 * @param o1 First object (numerator)
1872 * @param o2 Second object (denominator)
1873 * @param max Maximum diviation from 1
e1f47419 1874 *
ffca499d 1875 * @return o1 / o2
e1f47419 1876 */
1877 TH1* Ratio(const TObject* o1, const TObject* o2, Double_t& max) const
1878 {
b30dee70 1879 if (!o1 || !o2) return 0;
1880
5bb5d1f6 1881 TH1* r = 0;
1882 const TAttMarker* m1 = 0;
1883 const TAttMarker* m2 = 0;
e1f47419 1884 const TH1* h1 = dynamic_cast<const TH1*>(o1);
1885 if (h1) {
5bb5d1f6 1886 m1 = h1;
e1f47419 1887 const TH1* h2 = dynamic_cast<const TH1*>(o2);
ffca499d 1888 if (h2) {
5bb5d1f6 1889 m2 = h2;
1890 r = RatioHH(h1,h2);
ffca499d 1891 }
1892 else {
1893 const TGraph* g2 = dynamic_cast<const TGraph*>(o2);
1894 if (g2) {
5bb5d1f6 1895 m2 = g2;
1896 r = RatioHG(h1,g2);
ffca499d 1897 }
1898 }
e1f47419 1899 }
ffca499d 1900 else {
1901 const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(o1);
1902 if (g1) {
5bb5d1f6 1903 m1 = g1;
ffca499d 1904 const TGraphAsymmErrors* g2 =
1905 dynamic_cast<const TGraphAsymmErrors*>(o2);
1906 if (g2) {
5bb5d1f6 1907 m2 = g2;
1908 r = RatioGG(g1, g2);
ffca499d 1909 }
1910 }
e1f47419 1911 }
ffca499d 1912 if (!r) {
c8b1a7db 1913 // Warning("Ratio", "Don't know how to divide a %s (%s) with a %s (%s)",
1914 // o1->ClassName(),o1->GetName(),o2->ClassName(),o2->GetName());
ffca499d 1915 return 0;
1916 }
1917 // Check that the histogram isn't empty
1918 if (r->GetEntries() <= 0) {
1919 delete r;
1920 r = 0;
1921 }
5bb5d1f6 1922 if (r) {
9ecab72f 1923 r->SetMarkerStyle(m2->GetMarkerStyle());
1924 r->SetMarkerColor(m1->GetMarkerColor());
cd3b3781 1925 if (TString(o2->GetName()).Contains("truth", TString::kIgnoreCase))
1926 r->SetMarkerStyle(m1->GetMarkerStyle());
5bb5d1f6 1927 r->SetMarkerSize(0.9*m1->GetMarkerSize());
1928 r->SetName(Form("%s_over_%s", o1->GetName(), o2->GetName()));
1929 r->SetTitle(Form("%s / %s", o1->GetTitle(), o2->GetTitle()));
1930 r->SetDirectory(0);
1931 max = TMath::Max(RatioMax(r), max);
1932 }
1933
ffca499d 1934 return r;
e1f47419 1935 }
1936 //__________________________________________________________________
b2e7f2d6 1937 /**
1938 * Compute the ratio of @a h to @a g. @a g is evaluated at the bin
1939 * centers of @a h
1940 *
1941 * @param h Numerator
1942 * @param g Divisor
1943 *
1944 * @return h/g
1945 */
5bb5d1f6 1946 TH1* RatioHG(const TH1* h, const TGraph* g) const
b2e7f2d6 1947 {
1948 if (!h || !g) return 0;
1949
1950 TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
b2e7f2d6 1951 ret->Reset();
b2e7f2d6 1952 Double_t xlow = g->GetX()[0];
1953 Double_t xhigh = g->GetX()[g->GetN()-1];
1954 if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
1955
1956 for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1957 Double_t c = h->GetBinContent(i);
1958 if (c <= 0) continue;
1959
1960 Double_t x = h->GetBinCenter(i);
1961 if (x < xlow || x > xhigh) continue;
1962
1963 Double_t f = g->Eval(x);
1964 if (f <= 0) continue;
1965
1966 ret->SetBinContent(i, c / f);
1967 ret->SetBinError(i, h->GetBinError(i) / f);
1968 }
b2e7f2d6 1969 return ret;
1970 }
1971 //__________________________________________________________________
1972 /**
1973 * Make the ratio of h1 to h2
1974 *
1975 * @param h1 First histogram (numerator)
1976 * @param h2 Second histogram (denominator)
1977 *
1978 * @return h1 / h2
1979 */
5bb5d1f6 1980 TH1* RatioHH(const TH1* h1, const TH1* h2) const
b2e7f2d6 1981 {
1982 if (!h1 || !h2) return 0;
c8b1a7db 1983 Bool_t bad = false;
1984 if (h1->GetNbinsX() != h2->GetNbinsX()) {
1985 Error("RatioHH", "They have differnet number of bins");
1986 bad = true;
1987 }
1988 for (Int_t i = 1; i <= h1->GetNbinsX(); i++) {
1989 if (h1->GetXaxis()->GetBinLowEdge(i) !=
1990 h2->GetXaxis()->GetBinLowEdge(i)) {
1991 // Error("RatioHH", "They have incompatible variable bins");
1992 bad = true;
1993 break;
1994 }
1995 }
1996 if (bad) return 0;
1997
5bb5d1f6 1998 TH1* t1 = static_cast<TH1*>(h1->Clone("tmp"));
c8b1a7db 1999 // Printf("Dividing %s with %s", h1->GetName(), h2->GetName());
b2e7f2d6 2000 t1->Divide(h2);
b2e7f2d6 2001 return t1;
2002 }
2003 //__________________________________________________________________
2004 /**
2005 * Calculate the ratio of two graphs - g1 / g2
2006 *
2007 * @param g1 Numerator
2008 * @param g2 Denominator
2009 *
2010 * @return g1 / g2 in a histogram
2011 */
ffca499d 2012 TH1* RatioGG(const TGraphAsymmErrors* g1,
5bb5d1f6 2013 const TGraphAsymmErrors* g2) const
b2e7f2d6 2014 {
2015 Int_t nBins = g1->GetN();
2016 TArrayF bins(nBins+1);
2017 Double_t dx = 0;
2018 for (Int_t i = 0; i < nBins; i++) {
2019 Double_t x = g1->GetX()[i];
2020 Double_t exl = g1->GetEXlow()[i];
2021 Double_t exh = g1->GetEXhigh()[i];
2022 bins.fArray[i] = x-exl;
2023 bins.fArray[i+1] = x+exh;
2024 Double_t dxi = exh+exl;
2025 if (i == 0) dx = dxi;
2026 else if (dxi != dx) dx = 0;
2027 }
b2e7f2d6 2028 TH1* h = 0;
2029 if (dx != 0) {
5bb5d1f6 2030 h = new TH1F("tmp", "tmp", nBins, bins[0], bins[nBins]);
b2e7f2d6 2031 }
2032 else {
5bb5d1f6 2033 h = new TH1F("tmp", "tmp", nBins, bins.fArray);
b2e7f2d6 2034 }
b2e7f2d6 2035
2036 Double_t low = g2->GetX()[0];
2037 Double_t high = g2->GetX()[g2->GetN()-1];
2038 if (low > high) { Double_t t = low; low = high; high = t; }
2039 for (Int_t i = 0; i < nBins; i++) {
2040 Double_t x = g1->GetX()[i];
2041 if (x < low-dx || x > high+dx) continue;
2042 Double_t c1 = g1->GetY()[i];
2043 Double_t e1 = g1->GetErrorY(i);
2044 Double_t c2 = g2->Eval(x);
2045
2046 h->SetBinContent(i+1, c1 / c2);
2047 h->SetBinError(i+1, e1 / c2);
2048 }
b2e7f2d6 2049 return h;
2f92ff0e 2050 }
2051 /* @} */
2052 //==================================================================
2053 /**
2054 * @{
2055 * @name Graphics utility functions
2056 */
2057 /**
2058 * Find an X axis in a pad
2059 *
2060 * @param p Pad
2061 * @param name Histogram to find axis for
2062 *
2063 * @return Found axis or null
2064 */
2065 TAxis* FindXAxis(TVirtualPad* p, const char* name)
b2e7f2d6 2066 {
2f92ff0e 2067 TObject* o = p->GetListOfPrimitives()->FindObject(name);
2068 if (!o) {
2069 Warning("FindXAxis", "%s not found in pad", name);
2070 return 0;
b2e7f2d6 2071 }
2f92ff0e 2072 THStack* stack = dynamic_cast<THStack*>(o);
2073 if (!stack) {
2074 Warning("FindXAxis", "%s is not a THStack", name);
2075 return 0;
b2e7f2d6 2076 }
2f92ff0e 2077 if (!stack->GetHistogram()) {
2078 Warning("FindXAxis", "%s has no histogram", name);
2079 return 0;
b2e7f2d6 2080 }
2f92ff0e 2081 TAxis* ret = stack->GetHistogram()->GetXaxis();
2082 return ret;
b2e7f2d6 2083 }
2f92ff0e 2084
b2e7f2d6 2085 //__________________________________________________________________
2f92ff0e 2086 /**
2087 * Fix the apperance of the axis in a stack
2088 *
2089 * @param stack stack of histogram
c6115ede 2090 * @param yd How the canvas is cut
2f92ff0e 2091 * @param ytitle Y axis title
2f92ff0e 2092 * @param ynDiv Divisions on Y axis
c6115ede 2093 * @param force Whether to draw the stack first or not
2f92ff0e 2094 */
797161e8 2095 void FixAxis(THStack* stack, Double_t yd, const char* ytitle,
2f92ff0e 2096 Int_t ynDiv=210, Bool_t force=true)
b2e7f2d6 2097 {
e308a636 2098 if (!stack) {
2099 Warning("FixAxis", "No stack passed for %s!", ytitle);
2100 return;
2101 }
2f92ff0e 2102 if (force) stack->Draw("nostack e1");
b2e7f2d6 2103
2f92ff0e 2104 TH1* h = stack->GetHistogram();
e308a636 2105 if (!h) {
2106 Warning("FixAxis", "Stack %s has no histogram", stack->GetName());
2107 return;
2108 }
797161e8 2109 Double_t s = 1/yd/1.2;
2110 // Info("FixAxis", "for %s, s=1/%f=%f", stack->GetName(), yd, s);
2111
d43c6cd0 2112 h->SetXTitle("#it{#eta}");
2f92ff0e 2113 h->SetYTitle(ytitle);
2114 TAxis* xa = h->GetXaxis();
2115 TAxis* ya = h->GetYaxis();
797161e8 2116
2117 // Int_t npixels = 20
2118 // Float_t dy = gPad->AbsPixeltoY(0) - gPad->AbsPixeltoY(npixels);
2119 // Float_t ts = dy/(gPad->GetY2() - gPad->GetY1());
2120
2f92ff0e 2121 if (xa) {
797161e8 2122 // xa->SetTitle(h->GetXTitle());
2f92ff0e 2123 // xa->SetTicks("+-");
2124 xa->SetTitleSize(s*xa->GetTitleSize());
2125 xa->SetLabelSize(s*xa->GetLabelSize());
2126 xa->SetTickLength(s*xa->GetTickLength());
797161e8 2127 // xa->SetTitleOffset(xa->GetTitleOffset()/s);
e308a636 2128
2129 if (stack != fResults) {
2130 TAxis* rxa = fResults->GetXaxis();
2131 xa->Set(rxa->GetNbins(), rxa->GetXmin(), rxa->GetXmax());
2132 }
2f92ff0e 2133 }
2134 if (ya) {
797161e8 2135 // ya->SetTitle(h->GetYTitle());
2f92ff0e 2136 ya->SetDecimals();
2137 // ya->SetTicks("+-");
2138 ya->SetNdivisions(ynDiv);
2139 ya->SetTitleSize(s*ya->GetTitleSize());
d43c6cd0 2140 ya->SetTitleOffset(1.15*ya->GetTitleOffset()/s);
2f92ff0e 2141 ya->SetLabelSize(s*ya->GetLabelSize());
b2e7f2d6 2142 }
b2e7f2d6 2143 }
797161e8 2144 //__________________________________________________________________
2145 /**
2146 * Merge two histograms into one
2147 *
2148 * @param cen Central part
2149 * @param fwd Forward part
2150 * @param xlow On return, lower eta bound
2151 * @param xhigh On return, upper eta bound
2152 *
2153 * @return Newly allocated histogram or null
2154 */
2155 TH1*
2156 Merge(const TH1* cen, const TH1* fwd, Double_t& xlow, Double_t& xhigh)
2157 {
2158 TH1* tmp = static_cast<TH1*>(fwd->Clone("tmp"));
2159 TString name(fwd->GetName());
2160 name.ReplaceAll("Forward", "Merged");
2161 tmp->SetName(name);
2162
2163 // tmp->SetMarkerStyle(28);
2164 // tmp->SetMarkerColor(kBlack);
2165 tmp->SetDirectory(0);
2166 xlow = 100;
2167 xhigh = -100;
2168 for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
2169 Double_t cc = cen->GetBinContent(i);
2170 Double_t cf = fwd->GetBinContent(i);
2171 Double_t ec = cen->GetBinError(i);
2172 Double_t ef = fwd->GetBinError(i);
2173 Double_t nc = cf;
2174 Double_t ne = ef;
2175 if (cc < 0.001 && cf < 0.01) continue;
2176 xlow = TMath::Min(tmp->GetXaxis()->GetBinLowEdge(i),xlow);
2177 xhigh = TMath::Max(tmp->GetXaxis()->GetBinUpEdge(i),xhigh);
2178 if (cc > 0.001) {
2179 nc = cc;
2180 ne = ec;
2181 if (cf > 0.001) {
2182 nc = (cf + cc) / 2;
2183 ne = TMath::Sqrt(ec*ec + ef*ef);
2184 }
2185 }
2186 tmp->SetBinContent(i, nc);
2187 tmp->SetBinError(i, ne);
2188 }
2189 return tmp;
2190 }
2191 //____________________________________________________________________
2192 /**
2193 * Fit @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$ to histogram data
2194 *
2195 * @param tmp Histogram
2196 * @param xlow Lower x bound
2197 * @param xhigh Upper x bound
2198 *
2199 * @return Fitted function
2200 */
2201 TF1*
2202 FitMerged(TH1* tmp, Double_t xlow, Double_t xhigh)
2203 {
2204 TF1* tmpf = new TF1("tmpf", "gaus", xlow, xhigh);
2205 tmp->Fit(tmpf, "NQ", "");
2206 tmp->SetDirectory(0);
2207
2208 TF1* fit = new TF1("f", myFunc, xlow, xhigh, 4);
2209 fit->SetParNames("a_{1}", "a_{2}", "#sigma_{1}", "#sigma_{2}");
2210 fit->SetParameters(tmpf->GetParameter(0),
2211 .2,
2212 tmpf->GetParameter(2),
2213 tmpf->GetParameter(2)/4);
2214 fit->SetParLimits(3, 0, 100);
2215 fit->SetParLimits(4, 0, 100);
2216 tmp->Fit(fit,"0W","");
2217
2218 delete tmpf;
2219 return fit;
2220 }
2221 //____________________________________________________________________
2222 /**
2223 * Make band of systematic errors
2224 *
2225 * @param tmp Histogram
c6115ede 2226 * @param cen Central
2227 * @param fwd Forward
797161e8 2228 * @param fit Fit
2229 */
2230 void
2231 MakeSysError(TH1* tmp, TH1* cen, TH1* fwd, TF1* fit)
2232 {
2233 for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
2234 Double_t tc = tmp->GetBinContent(i);
2235 if (tc < 0.01) continue;
2236 Double_t fc = fwd->GetBinContent(i);
2237 Double_t cc = cen->GetBinContent(i);
2238 Double_t sysErr = fFwdSysErr;
2239 if (cc > .01 && fc > 0.01)
2240 sysErr = (fFwdSysErr+fCenSysErr) / 2;
2241 else if (cc > .01)
2242 sysErr = fCenSysErr;
2243 Double_t x = tmp->GetXaxis()->GetBinCenter(i);
2244 Double_t y = fit->Eval(x);
2245 tmp->SetBinContent(i, y);
2246 tmp->SetBinError(i,sysErr*y);
2247 }
2248 TString name(tmp->GetName());
2249 name.ReplaceAll("Merged", "SysError");
2250 tmp->SetName(name);
e42bc740 2251 tmp->SetMarkerColor(SYSERR_COLOR);
2252 tmp->SetLineColor(SYSERR_COLOR);
2253 tmp->SetFillColor(SYSERR_COLOR);
2254 tmp->SetFillStyle(SYSERR_STYLE);
797161e8 2255 tmp->SetMarkerStyle(0);
2256 tmp->SetLineWidth(0);
2257 }
2258 void CorrectForward(TH1* h) const
2259 {
2260 if (!fRemoveOuters) return;
2261
2262 for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
2263 Double_t eta = h->GetBinCenter(i);
2264 if (TMath::Abs(eta) < 2.3) {
2265 h->SetBinContent(i, 0);
2266 h->SetBinError(i, 0);
2267 }
2268 }
2269 }
2270 void CorrectCentral(TH1* h) const
2271 {
2272 if (fClusterScale.IsNull()) return;
2273 TString t(h->GetTitle());
2274 Info("CorrectCentral", "Replacing Central with Tracklets in %s", t.Data());
2275 t.ReplaceAll("Central", "Tracklets");
2276 h->SetTitle(t);
2277
2278 TF1* cf = new TF1("clusterScale", fClusterScale);
2279 for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
2280 Double_t eta = h->GetBinCenter(i);
2281 Double_t f = cf->Eval(eta);
2282 Double_t c = h->GetBinContent(i);
2283 if (f < .1) f = 1;
2284 h->SetBinContent(i, c / f);
2285 }
2286 delete cf;
2287 }
33ab4a11 2288 //____________________________________________________________________
2289 void Export(const char* basename)
2290 {
2291 TString bname(basename);
2292 bname.ReplaceAll(" ", "_");
2293 bname.ReplaceAll("-", "_");
2294 TString fname(Form("%s.C", bname.Data()));
2295
298a7264 2296 std::ofstream outf(fname.Data());
2297 if (!outf) {
33ab4a11 2298 Error("Export", "Failed to open output file %s", fname.Data());
2299 return;
2300 }
298a7264 2301 outf << "// Create by dNdetaDrawer\n"
2302 << "void " << bname << "(THStack* stack, TLegend* l, Int_t m)\n"
2303 << "{"
2304 << " Int_t ma[] = { 24, 25, 26, 32,\n"
2305 << " 20, 21, 22, 33,\n"
2306 << " 34, 30, 29, 0, \n"
2307 << " 23, 27 };\n"
2308 << " Int_t mm = ((m < 20 || m > 34) ? 0 : ma[m-20]);\n\n";
33ab4a11 2309 TList* hists = fResults->GetHists();
2310 TIter next(hists);
2311 TH1* hist = 0;
2312 while ((hist = static_cast<TH1*>(next()))) {
2313 TString hname = hist->GetName();
2314 hname.Append(Form("_%04x", (gRandom->Integer(0xffff) & 0xffff)));
2315 hist->SetName(hname);
2316 hist->GetListOfFunctions()->Clear();
298a7264 2317 hist->SavePrimitive(outf, "nodraw");
33ab4a11 2318 bool mirror = hname.Contains("mirror");
2319 bool syserr = hname.Contains("SysError");
2320 if (!syserr)
298a7264 2321 outf << " " << hname << "->SetMarkerStyle("
2322 << (mirror ? "mm" : "m") << ");\n";
33ab4a11 2323 else
298a7264 2324 outf << " " << hname << "->SetMarkerStyle(1);\n";
2325 outf << " stack->Add(" << hname
2326 << (syserr ? ",\"e5\"" : "") << ");\n\n";
33ab4a11 2327 }
2328 UShort_t snn = fSNNString->GetUniqueID();
2329 // const char* sys = fSysString->GetTitle();
2330 TString eS;
2331 if (snn == 2750) snn = 2760;
2332 if (snn < 1000) eS = Form("%3dGeV", snn);
2333 else if (snn % 1000 == 0) eS = Form("%dTeV", snn/1000);
2334 else eS = Form("%4.2fTeV", float(snn)/1000);
298a7264 2335 outf << " if (l) {\n"
2336 << " TLegendEntry* e = l->AddEntry(\"\",\"" << eS << "\",\"pl\");\n"
2337 << " e->SetMarkerStyle(m);\n"
2338 << " e->SetMarkerColor(kBlack);\n"
2339 << " }\n"
2340 << "}\n" << std::endl;
33ab4a11 2341 }
09d8cf0a 2342 /* @} */
73b32206 2343 /**
2344 * Check if we have centrality dependent information, and we're not
2345 * forcing to use minimum bias
2346 *
2347 *
2348 * @return True if we should do centrality dependent ploting
2349 */
c8b1a7db 2350 Bool_t HasCent() const
2351 {
2352 return fCentAxis && fCentAxis->GetNbins() > 0 && !fForceMB;
2353 }
2f92ff0e 2354
b2e7f2d6 2355
2356
2357 //__________________________________________________________________
33ab4a11 2358 /**
2359 * @{
2360 * @name Options
2361 */
e308a636 2362 Bool_t fShowRatios; // Show ratios
2363 Bool_t fShowLeftRight;// Show asymmetry
f709afcf 2364 Bool_t fShowRings; // Show rings too
33ab4a11 2365 Bool_t fExport; // Export results to file
e308a636 2366 Bool_t fCutEdges; // Whether to cut edges
33ab4a11 2367 Bool_t fRemoveOuters; // Whether to remove outers
2368 UShort_t fShowOthers; // Show other data
298a7264 2369 Bool_t fMirror; // Whether to mirror
2370 Bool_t fForceMB; // Force min-bias
73b32206 2371 Bool_t fAddExec; // Add code to do combined zooms
d43c6cd0 2372 Bool_t fSaveAsPDF; // Flag whether to save as PDF
2373 Bool_t fSaveAsPNG; // Flag whether to save as PNG
2374 Bool_t fSaveAsROOT; // Flag whether to save as ROOT file
2375 Bool_t fSaveAsScript; // Flag whether to save as ROOT script
33ab4a11 2376 /* @} */
2377 /**
2378 * @{
2379 * @name Settings
2380 */
2381 UShort_t fRebin; // Rebinning factor
2382 Double_t fFwdSysErr; // Systematic error in forward range
2383 Double_t fCenSysErr; // Systematic error in central range
e308a636 2384 TString fTitle; // Title on plot
d43c6cd0 2385 TString fBase; // Base name of output
33ab4a11 2386 TString fClusterScale; // Scaling of clusters to tracklets
b615e211 2387 TString fFinalMC; // Final MC correction file name
2388 TString fEmpirical; // Empirical correction file name
33ab4a11 2389 /* @} */
2390 /**
2391 * @{
2392 * @name Read (or set) information
2393 */
e308a636 2394 TNamed* fTrigString; // Trigger string (read, or set)
0be6c8cd 2395 TNamed* fNormString; // Normalisation string (read, or set)
e308a636 2396 TNamed* fSNNString; // Energy string (read, or set)
2397 TNamed* fSysString; // Collision system string (read or set)
2398 TAxis* fVtxAxis; // Vertex cuts (read or set)
2399 TAxis* fCentAxis; // Centrality axis
bfab35d9 2400 Float_t fTriggerEff; // Trigger efficiency
d43c6cd0 2401 Bool_t fExtTriggerEff;// True if read externally
33ab4a11 2402 /* @} */
2403 /**
2404 * @{
2405 * @name Resulting plots
2406 */
e308a636 2407 THStack* fResults; // Stack of results
2408 THStack* fRatios; // Stack of ratios
2409 THStack* fLeftRight; // Left-right asymmetry
2410 TMultiGraph* fOthers; // Older data
2411 TH1* fTriggers; // Number of triggers
9ecab72f 2412 TH1* fTruth; // Pointer to truth
33ab4a11 2413 /* @} */
e308a636 2414 RangeParam* fRangeParam; // Parameter object for range zoom
2e6b5465 2415 Bool_t fOld;
d43c6cd0 2416
2417 static const Float_t kRightMargin;
2418 static const Int_t kFont;
2419 static const Int_t kAliceBlue;
2420 static const Int_t kAliceRed;
2421 static const Int_t kAlicePurple;
2422 static const Int_t kAliceYellow;
b2e7f2d6 2423};
2424
d43c6cd0 2425const Float_t dNdetaDrawer::kRightMargin = 0.02;
2426const Int_t dNdetaDrawer::kFont = 42; // 132 for serif
2427const Int_t dNdetaDrawer::kAliceBlue = TColor::GetColor(40, 58, 68);
2428const Int_t dNdetaDrawer::kAliceRed = TColor::GetColor(226, 0, 26);
2429const Int_t dNdetaDrawer::kAlicePurple = TColor::GetColor(202, 71, 67);
2430const Int_t dNdetaDrawer::kAliceYellow = TColor::GetColor(238, 125, 17);
2431
2432
797161e8 2433//____________________________________________________________________
2434/**
2435 * Function to calculate
2436 * @f[
2437 * g(x;A_1,A_2,\sigma_1,\sigma_2) =
2438 * A_1\left(\frac{1}{2\pi\sigma_1}e^{(x/\sigma_1)^2} -
2439 * A_2\frac{1}{2\pi\sigma_2}e^{(x/\sigma_2)^2}\right)
2440 * @f]
2441 *
2442 * @param xp Pointer to x array
2443 * @param pp Pointer to parameter array
2444 *
2445 * @return @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$
2446 */
2447Double_t myFunc(Double_t* xp, Double_t* pp)
2448{
2449 Double_t x = xp[0];
2450 Double_t a1 = pp[0];
2451 Double_t a2 = pp[1];
2452 Double_t s1 = pp[2];
2453 Double_t s2 = pp[3];
2454 return a1*(TMath::Gaus(x, 0, s1) - a2 * TMath::Gaus(x, 0, s2));
2455}
2456
b2e7f2d6 2457//=== Stuff for auto zooming =========================================
290052e7 2458/**
2459 * Update canvas range
2460 *
2461 * @param p Parameter
2462 */
b2e7f2d6 2463void UpdateRange(dNdetaDrawer::RangeParam* p)
2464{
2465 if (!p) {
2466 Warning("UpdateRange", "No parameters %p", p);
2467 return;
2468 }
2469 if (!p->fMasterAxis) {
2470 Warning("UpdateRange", "No master axis %p", p->fMasterAxis);
2471 return;
2472 }
2473 Int_t first = p->fMasterAxis->GetFirst();
2474 Int_t last = p->fMasterAxis->GetLast();
2475 Double_t x1 = p->fMasterAxis->GetBinCenter(first);
2476 Double_t x2 = p->fMasterAxis->GetBinCenter(last);
b2e7f2d6 2477
2478 if (p->fSlave1Axis) {
2479 Int_t i1 = p->fSlave1Axis->FindBin(x1);
2480 Int_t i2 = p->fSlave1Axis->FindBin(x2);
2481 p->fSlave1Axis->SetRange(i1, i2);
2482 p->fSlave1Pad->Modified();
2483 p->fSlave1Pad->Update();
2484 }
2485 if (p->fSlave2Axis) {
2486 Int_t i1 = p->fSlave2Axis->FindBin(x1);
2487 Int_t i2 = p->fSlave2Axis->FindBin(x2);
2488 p->fSlave2Axis->SetRange(i1, i2);
2489 p->fSlave2Pad->Modified();
2490 p->fSlave2Pad->Update();
2491 }
2492 TCanvas* c = gPad->GetCanvas();
2493 c->cd();
2494}
2495
2496//____________________________________________________________________
290052e7 2497/**
2498 * Called when user changes X range
2499 *
2500 * @param p Parameter
2501 */
b2e7f2d6 2502void RangeExec(dNdetaDrawer::RangeParam* p)
2503{
2504 // Event types:
2505 // 51: Mouse move
2506 // 53:
2507 // 1: Button down
2508 // 21: Mouse drag
2509 // 11: Mouse release
2510 // dNdetaDrawer::RangeParam* p =
2511 // reinterpret_cast<dNdetaDrawer::RangeParam*>(addr);
2512 Int_t event = gPad->GetEvent();
2513 TObject *select = gPad->GetSelected();
2514 if (event == 53) {
2515 UpdateRange(p);
2516 return;
2517 }
2518 if (event != 11 || !select || select != p->fMasterAxis) return;
2519 UpdateRange(p);
2520}
2521
33ab4a11 2522//=== Steering functions
2523//==============================================
290052e7 2524/**
2525 * Display usage information
2526 *
2527 */
33ab4a11 2528void
2529Usage()
2530{
29befb1a 2531 printf("Usage: DrawdNdeta(FILE,TITLE,REBIN,OTHERS,FLAGS,"
d43c6cd0 2532 "SNN,SYS,TRIG,IPZMIN,IPZMAX,BASE,FMT)\n"
29befb1a 2533 " const char* FILE File name to open (\"forward_dndeta.root\")\n"
2534 " const char* TITLE Title to put on plot (\"\")\n"
2535 " UShort_t REBIN Rebinning factor (1)\n"
2536 " UShort_t OTHERS Other data to draw - more below (0x7)\n"
2537 " UShort_t FLAGS Visualisation flags - more below (0x7)\n"
2538 " UShort_t SYS (optional) 1:pp, 2:PbPb, 3:pPb\n"
2539 " UShort_t SNN (optional) sqrt(s_NN) in GeV\n"
2540 " UShort_t TRIG (optional) 1: INEL, 2: INEL>0, 4: NSD, ...\n"
bfab35d9 2541 " Float_t EFF (optional) Trigger efficiency\n"
29befb1a 2542 " Float_t IPZMIN (optional) Least z coordinate of IP\n"
d43c6cd0 2543 " Float_t IPZMAX (optional) Largest z coordinate of IP\n"
2544 " const char* BASE (optional) base name of output files\n"
2545 " UShort_t FMT (optional) Output formats\n"
2546 "\n"
2547 " OTHERS is a bit mask of\n"
73b32206 2548 " 0x1 Show UA5 data (INEL,NSD, ppbar, 900GeV)\n"
2549 " 0x2 Show CMS data (NSD, pp)\n"
2550 " 0x4 Show published ALICE data (INEL,INEL>0,NSD, pp)\n"
d43c6cd0 2551 " 0x8 Show event genertor data\n"
2552 "\n"
2553 " FLAGS is a bit mask of\n"
73b32206 2554 " 0x1 Show ratios of data to other data and possibly MC\n"
2555 " 0x2 Show left-right asymmetry\n"
2556 " 0x4 Show systematic error band\n"
2557 " 0x8 Show individual ring results (INEL only)\n"
2558 " 0x10 Cut edges when rebinning\n"
2559 " 0x20 Remove FMDxO points\n"
2560 " 0x40 Do not make our own canvas\n"
2561 " 0x80 Force use of MB\n"
2562 " 0x100 Mirror data\n"
2563 " 0x200 Apply `final MC' correction\n"
2564 " 0x400 Apply `Emperical' correction\n"
2565 " 0x800 Export results to script\n"
2566 " 0x1000 Add code to do combined zooms on eta axis\n"
d43c6cd0 2567 " 0x2000 Assume old-style input\n"
2568 "\n"
2569 " FMD is a bit mask of\n"
2570 " 0x1 Make PNG output\n"
2571 " 0x2 Make PDF output\n"
2572 " 0x4 Make ROOT file output\n"
2573 " 0x8 Make ROOT script output\n"
2574 "\n"
29befb1a 2575 "0x200 requires the file forward_dndetamc.root\n"
2576 "0x400 requires the file EmpiricalCorrection.root\n"
2577 "To specify that you want ratios, force MB, apply empirical "
2578 "correction, and export to script, set flags to\n\n"
2579 " 0x1|0x80|0x400|0x800=0xC81\n\n"
2580 );
33ab4a11 2581}
2582
2583//____________________________________________________________________
ffca499d 2584/**
2585 * Draw @f$ dN/d\eta@f$
2586 *
2587 * @param filename File name
ffca499d 2588 * @param title Title
2589 * @param rebin Rebinning factor
c6115ede 2590 * @param others What other data to show
2591 * @param flags Flags
ffca499d 2592 * @param sNN (optional) Collision energy [GeV]
2593 * @param sys (optional) Collision system (1: pp, 2: PbPb)
2594 * @param trg (optional) Trigger (1: INEL, 2: INEL>0, 4: NSD)
c8b1a7db 2595 * @param eff (optional) Trigger efficiency
ffca499d 2596 * @param vzMin Least @f$ v_z@f$
2597 * @param vzMax Largest @f$ v_z@f$
c8b1a7db 2598 * @param base Base name
2599 * @param outflg Output flags
ffca499d 2600 *
bd6f5206 2601 * @ingroup pwglf_forward_dndeta
ffca499d 2602 */
b2e7f2d6 2603void
2604DrawdNdeta(const char* filename="forward_dndeta.root",
bad9a3c1 2605 const char* title="",
b2e7f2d6 2606 UShort_t rebin=5,
797161e8 2607 UShort_t others=0x7,
298a7264 2608 UShort_t flags=0x187,
fe52e455 2609 UShort_t sNN=0,
2610 UShort_t sys=0,
ffca499d 2611 UShort_t trg=0,
bfab35d9 2612 Float_t eff=0,
fe52e455 2613 Float_t vzMin=999,
d43c6cd0 2614 Float_t vzMax=-999,
2615 const char* base="",
2616 UShort_t outflg=0xf)
b2e7f2d6 2617{
33ab4a11 2618 TString fname(filename);
2619 fname.ToLower();
2620 if (fname.CompareTo("help") == 0 ||
2621 fname.CompareTo("--help") == 0) {
2622 Usage();
2623 return;
2624 }
ffca499d 2625 dNdetaDrawer* pd = new dNdetaDrawer;
2626 dNdetaDrawer& d = *pd;
b2e7f2d6 2627 d.SetRebin(rebin);
b2e7f2d6 2628 d.SetTitle(title);
797161e8 2629 d.SetShowOthers(others);
2630 d.SetShowRatios(flags & 0x1);
2631 d.SetShowLeftRight(flags & 0x2);
2632 d.SetForwardSysError(flags & 0x4 ? 0.076 : 0);
2633 d.SetShowRings(flags & 0x8);
2634 d.SetCutEdges(flags & 0x10);
33ab4a11 2635 d.fRemoveOuters = (flags & 0x20);
2636 d.SetExport(flags & 0x40);
298a7264 2637 d.SetForceMB(flags & 0x80);
2638 d.SetMirror(flags & 0x100);
b615e211 2639 d.SetFinalMC(flags & 0x200 ? "forward_dndetamc.root" : "");
2640 d.SetEmpirical(flags & 0x400 ? "EmpiricalCorrection.root" : "");
29befb1a 2641 d.SetExport(flags & 0x800);
73b32206 2642 d.SetAddExec(flags & 0x1000);
2643 d.SetOld(flags & 0x2000);
d43c6cd0 2644 d.SetBase(base);
2645 d.SetSaveAsPNG(outflg & 0x1);
2646 d.SetSaveAsPDF(outflg & 0x2);
2647 d.SetSaveAsROOT(outflg & 0x4);
2648 d.SetSaveAsScript(outflg & 0x8);
797161e8 2649 // d.fClusterScale = "1.06 -0.003*x +0.0119*x*x";
2f92ff0e 2650 // Do the below if your input data does not contain these settings
fe52e455 2651 if (sNN > 0) d.SetSNN(sNN); // Collision energy per nucleon pair (GeV)
2652 if (sys > 0) d.SetSys(sys); // Collision system (1:pp, 2:PbPB)
2653 if (trg > 0) d.SetTrigger(trg); // Collision trigger (1:INEL, 2:INEL>0, 4:NSD)
bfab35d9 2654 if (eff > 0) d.SetTriggerEfficiency(eff); // Trigger efficiency
fe52e455 2655 if (vzMin < 999 && vzMax > -999)
2656 d.SetVertexRange(vzMin,vzMax); // Collision vertex range (cm)
b2e7f2d6 2657 d.Run(filename);
2658}
2f92ff0e 2659//____________________________________________________________________
2660//
2661// EOF
2662//
2663