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