PWG2/FORWARD/analysis2/DrawdNdeta.C

   1 #include <TH1.h>
   2 #include <THStack.h>
   3 #include <TGraphErrors.h>
   4 #include <TGraphAsymmErrors.h>
   5 #include <TMultiGraph.h>
   6 #include <TFile.h>
   7 #include <TList.h>
   8 #include <TString.h>
   9 #include <TError.h>
  10 #include <TSystem.h>
  11 #include <TROOT.h>
  12 #include <TMath.h>
  13 #include <TCanvas.h>
  14 #include <TPad.h>
  15 #include <TStyle.h>
  16 #include <TLegend.h>
  17 #include <TLegendEntry.h>
  18 #include <TLatex.h>
  19 #include <TImage.h>
  20
  21 struct dNdetaDrawer
  22 {
  23   struct RangeParam
  24   {
  25     TAxis*       fMasterAxis; // Master axis
  26     TAxis*       fSlave1Axis; // First slave axis
  27     TVirtualPad* fSlave1Pad;  // First slave pad
  28     TAxis*       fSlave2Axis; // Second slave axis
  29     TVirtualPad* fSlave2Pad;  // Second slave pad
  30   };
  31   //__________________________________________________________________
  32   dNdetaDrawer()
  33     : fShowOthers(false),       // Bool_t
  34       fShowAlice(false),        // Bool_t
  35       fShowRatios(false),       // Bool_t
  36       fShowLeftRight(false),    // Bool_t
  37       fRebin(5),                // UShort_t
  38       fCutEdges(false),         // Bool_t
  39       fTitle(""),               // TString
  40       fHHDFile(""),             // TString
  41       fTrigString(0),           // TNamed*
  42       fSNNString(0),            // TNamed*
  43       fSysString(0),            // TNamed*
  44       fVtxAxis(0),              // TAxis*
  45       fForward(0),              // TH1*
  46       fForwardMC(0),            // TH1*
  47       fForwardHHD(0),           // TH1*
  48       fTruth(0),                // TH1*
  49       fCentral(0),              // TH1*
  50       fForwardSym(0),           // TH1*
  51       fForwardMCSym(0),         // TH1*
  52       fForwardHHDSym(0),        // TH1*
  53       fTriggers(0),             // TH1*
  54       fRangeParam(0)
  55
  56   {
  57     fRangeParam = new RangeParam;
  58     fRangeParam->fMasterAxis = 0;
  59     fRangeParam->fSlave1Axis = 0;
  60     fRangeParam->fSlave1Pad  = 0;
  61     fRangeParam->fSlave2Axis = 0;
  62     fRangeParam->fSlave2Pad  = 0;
  63   }
  64   void SetShowOthers(Bool_t x)    { fShowOthers = x; }
  65   void SetShowAlice(Bool_t x)     { fShowAlice = x; }
  66   void SetShowRatios(Bool_t x)    { fShowRatios = x; }
  67   void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; }
  68   void SetRebin(UShort_t x)       { fRebin = x; }
  69   void SetCutEdges(Bool_t x)      { fCutEdges = x; }
  70   void SetTitle(TString x)        { fTitle = x; }
  71   void SetHHDFile(const char* fn) { fHHDFile = fn; }
  72
  73   //__________________________________________________________________
  74   void Run(const char* filename)
  75   {
  76     if (!Open(filename)) return;
  77
  78     Double_t max = 0;
  79
  80     // Create our stack of results
  81     THStack* results = StackResults(max);
  82
  83     // Create our stack of other results
  84     TMultiGraph* other = 0;
  85     if (fShowOthers || fShowAlice) other = StackOther(max);
  86
  87     Double_t smax = 0;
  88     THStack* ratios = 0;
  89     if (fShowRatios) ratios = StackRatios(other, smax);
  90
  91     Double_t amax = 0;
  92     THStack* leftright = 0;
  93     if (fShowLeftRight) leftright = StackLeftRight(amax);
  94
  95     Plot(results, other, max, ratios, smax, leftright, amax);
  96   }
  97
  98   //__________________________________________________________________
  99   Bool_t Open(const char* filename)
 100   {
 101     TFile* file = TFile::Open(filename, "READ");
 102     if (!file) {
 103       Error("Open", "Cannot open %s", filename);
 104       return false;
 105     }
 106
 107     TList* results = static_cast<TList*>(file->Get("ForwardResults"));
 108     if (!results) {
 109       Error("Open", "Couldn't find list ForwardResults");
 110       return false;
 111     }
 112
 113     fForward   = GetResult(results, "dndetaForward");
 114     fForwardMC = GetResult(results, "dndetaForwardMC");
 115     fTruth     = GetResult(results, "dndetaTruth");
 116     fCentral   = GetResult(results, "dndetaCentral");
 117
 118     fTrigString = static_cast<TNamed*>(results->FindObject("trigString"));
 119     fSNNString  = static_cast<TNamed*>(results->FindObject("sNN"));
 120     fSysString  = static_cast<TNamed*>(results->FindObject("sys"));
 121     fVtxAxis    = static_cast<TAxis*>(results->FindObject("vtxAxis"));
 122
 123     TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
 124     if (sums)
 125       fTriggers = GetResult(sums, "triggers");
 126
 127     if (!fForward) {
 128       Error("Open", "Couldn't find the result of the forward analysis");
 129       return false;
 130     }
 131     file->Close();
 132
 133
 134     fForwardHHD = GetHHD();
 135
 136     return true;
 137   }
 138   //__________________________________________________________________
 139   TH1* GetResult(TList* list, const char* name) const
 140   {
 141     if (!list) return 0;
 142
 143     TH1* ret = static_cast<TH1*>(list->FindObject(name));
 144     if (!ret)
 145       Warning("GetResult", "Histogram %s not found", name);
 146
 147     return ret;
 148   }
 149   //__________________________________________________________________
 150   /**
 151    * Get the result from previous analysis code
 152    *
 153    * @param fn  File to open
 154    * @param nsd Whether this is NSD
 155    *
 156    * @return null or result of previous analysis code
 157    */
 158   TH1* GetHHD()
 159   {
 160     if (fHHDFile.IsNull()) return 0;
 161     const char* fn = fHHDFile.Data();
 162     Bool_t nsd = (fTrigString ? fTrigString->GetUniqueID() & 0x4 : false);
 163     TDirectory* savdir = gDirectory;
 164     if (gSystem->AccessPathName(fn)) {
 165       Warning("GetHHD", "Output of HHD analysis (%s) not available", fn);
 166       return 0;
 167     }
 168     TFile* file = TFile::Open(fn, "READ");
 169     if (!file) {
 170       Warning("GetHHD", "couldn't open HHD file %s", fn);
 171       return 0;
 172     }
 173     TString hist(Form("dNdeta_dNdeta%s", (nsd ? "NSD" : "")));
 174     TH1* h = static_cast<TH1*>(file->Get(hist.Data()));
 175     if (!h) {
 176       Warning("GetHHD", "Couldn't find HHD histogram %s in %s",
 177               hist.Data(), fn);
 178       file->Close();
 179       savdir->cd();
 180       return 0;
 181     }
 182     TH1* r = static_cast<TH1*>(h->Clone("dndeta_hhd"));
 183     r->SetTitle("ALICE Forward (HHD)");
 184     r->SetFillStyle(0);
 185     r->SetFillColor(0);
 186     r->SetMarkerStyle(21);
 187     r->SetMarkerColor(kPink+1);
 188     r->SetDirectory(0);
 189
 190     file->Close();
 191     savdir->cd();
 192     return r;
 193   }
 194   //__________________________________________________________________
 195   THStack* StackResults(Double_t& max)
 196   {
 197     THStack* stack = new THStack("results", "Stack of Results");
 198     max = TMath::Max(max, AddHistogram(stack, fTruth,      "e5 p"));
 199     max = TMath::Max(max, AddHistogram(stack, fForwardHHD, "", fForwardHHDSym));
 200     max = TMath::Max(max, AddHistogram(stack, fForwardMC,  "", fForwardMCSym));
 201     max = TMath::Max(max, AddHistogram(stack, fCentral,    ""));
 202     max = TMath::Max(max, AddHistogram(stack, fForward,    "", fForwardSym));
 203     return stack;
 204   }
 205   //__________________________________________________________________
 206   TMultiGraph* StackOther(Double_t& max) const
 207   {
 208     gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/OtherData.C");
 209     Int_t    error = 0;
 210     Bool_t   onlya = (fShowOthers ? false : true);
 211     Int_t    trg   = (fTrigString ? fTrigString->GetUniqueID() : 0);
 212     UShort_t snn   = (fSNNString  ? fSNNString->GetUniqueID() : 0);
 213     Long_t   ret   = gROOT->ProcessLine(Form("GetData(%d,%d,%d);",
 214                                              snn,trg,onlya));
 215     if (error) {
 216       Error("StackOther", "Failed to execute GetData(%d,%d,%d)",
 217             snn, trg, onlya);
 218       return 0;
 219     }
 220     if (!ret) {
 221       Warning("StackOther", "No other data found for sNN=%d, trigger=%d",
 222               snn, trg);
 223       return 0;
 224     }
 225     TMultiGraph* other = reinterpret_cast<TMultiGraph*>(ret);
 226
 227     TGraphAsymmErrors* o      = 0;
 228     TIter              next(other->GetListOfGraphs());
 229     while ((o = static_cast<TGraphAsymmErrors*>(next())))
 230       max = TMath::Max(max, TMath::MaxElement(o->GetN(), o->GetY()));
 231
 232     return other;
 233   }
 234   //__________________________________________________________________
 235   THStack* StackRatios(TMultiGraph* others, Double_t& max)
 236   {
 237     THStack* ratios = new THStack("ratios", "Ratios");
 238
 239     if (others) {
 240       TGraphAsymmErrors* ua5_1  = 0;
 241       TGraphAsymmErrors* ua5_2  = 0;
 242       TGraphAsymmErrors* alice  = 0;
 243       TGraphAsymmErrors* cms    = 0;
 244       TGraphAsymmErrors* o      = 0;
 245       TIter              nextG(others->GetListOfGraphs());
 246       while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
 247         ratios->Add(Ratio(fForward,          o, max));
 248         ratios->Add(Ratio(fForwardSym,       o, max));
 249         ratios->Add(Ratio(fForwardHHD,       o, max));
 250         ratios->Add(Ratio(fForwardHHDSym,    o, max));
 251         ratios->Add(Ratio(fCentral,          o, max));
 252         TString oName(o->GetName());
 253         oName.ToLower();
 254         if (oName.Contains("ua5"))  { if (ua5_1) ua5_2 = o; else ua5_1 = o; }
 255         if (oName.Contains("alice")) alice = o;
 256         if (oName.Contains("cms"))   cms = o;
 257       }
 258       if (ua5_1 && alice) ratios->Add(Ratio(alice, ua5_1, max));
 259       if (ua5_2 && alice) ratios->Add(Ratio(alice, ua5_2, max));
 260       if (cms   && alice) ratios->Add(Ratio(alice, cms,   max));
 261     }
 262
 263     // If we have data from HHD's analysis, then do the ratio of
 264     // our result to that data.
 265     if (fForwardHHD) {
 266       ratios->Add(Ratio(fForward,    fForwardHHD,    max));
 267       ratios->Add(Ratio(fForwardSym, fForwardHHDSym, max));
 268     }
 269
 270     // Do comparison to MC
 271     if (fForwardMC) {
 272       ratios->Add(Ratio(fForward,    fForwardMC,    max));
 273       ratios->Add(Ratio(fForwardSym, fForwardMCSym, max));
 274     }
 275
 276     // Check if we have ratios
 277     if (!ratios->GetHists() ||
 278         (ratios->GetHists()->GetEntries() <= 0)) {
 279       delete ratios;
 280       ratios = 0;
 281     }
 282     return ratios;
 283   }
 284   //__________________________________________________________________
 285   THStack* StackLeftRight(Double_t& max)
 286   {
 287     THStack* ret = new THStack("leftright", "Left-right asymmetry");
 288     ret->Add(Asymmetry(fForward,    max));
 289     ret->Add(Asymmetry(fForwardHHD, max));
 290     ret->Add(Asymmetry(fForwardMC,  max));
 291
 292     if (!ret->GetHists() ||
 293         (ret->GetHists()->GetEntries() <= 0)) {
 294       delete ret;
 295       ret = 0;
 296     }
 297     return ret;
 298   }
 299   //__________________________________________________________________
 300   TAxis* FindXAxis(TVirtualPad* p, const char* name)
 301   {
 302     TObject* o = p->GetListOfPrimitives()->FindObject(name);
 303     if (!o) {
 304       Warning("FindXAxis", "%s not found in pad", name);
 305       return 0;
 306     }
 307     THStack* stack = dynamic_cast<THStack*>(o);
 308     if (!stack) {
 309       Warning("FindXAxis", "%s is not a THStack", name);
 310       return 0;
 311     }
 312     if (!stack->GetHistogram()) {
 313       Warning("FindXAxis", "%s has no histogram", name);
 314       return 0;
 315     }
 316     TAxis* ret = stack->GetHistogram()->GetXaxis();
 317     return ret;
 318   }
 319   //__________________________________________________________________
 320   void Plot(THStack*     results,
 321             TMultiGraph* others,
 322             Double_t     max,
 323             THStack*     ratios,
 324             Double_t     rmax,
 325             THStack*     leftright,
 326             Double_t     amax)
 327   {
 328     gStyle->SetOptTitle(0);
 329     gStyle->SetTitleFont(132, "xyz");
 330     gStyle->SetLabelFont(132, "xyz");
 331
 332     Int_t    h = 800;
 333     Int_t    w = 800; // h / TMath::Sqrt(2);
 334     if (!ratios) w *= 1.4;
 335     if (!leftright) w *= 1.4;
 336     TCanvas* c = new TCanvas("Results", "Results", w, h);
 337     c->SetFillColor(0);
 338     c->SetBorderSize(0);
 339     c->SetBorderMode(0);
 340
 341     Double_t y1 = 0;
 342     Double_t y2 = 0;
 343     Double_t y3 = 0;
 344     if (ratios)    y1 = 0.3;
 345     if (leftright) {
 346       if (y1 > 0.0001) {
 347         y2 = 0.2;
 348         y1 = 0.4;
 349       }
 350       else {
 351         y1 = 0.2;
 352         y2 = y1;
 353       }
 354     }
 355     PlotResults(results, others, max, y1);
 356     c->cd();
 357
 358     PlotRatios(ratios, rmax, y2, y1);
 359     c->cd( );
 360
 361     PlotLeftRight(leftright, amax, y3, y2);
 362     c->cd();
 363   }
 364   //__________________________________________________________________
 365   void PlotResults(THStack* results, TMultiGraph* others,
 366                    Double_t max, Double_t yd)
 367   {
 368     // Make a sub-pad for the result itself
 369     TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0);
 370     p1->SetTopMargin(0.05);
 371     p1->SetBorderSize(0);
 372     p1->SetBorderMode(0);
 373     p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1);
 374     p1->SetRightMargin(0.05);
 375     p1->SetGridx();
 376     p1->SetTicks(1,1);
 377     p1->SetNumber(1);
 378     p1->Draw();
 379     p1->cd();
 380
 381     results->SetMaximum(1.15*max);
 382     results->SetMinimum(yd > 0.00001 ? -0.1 : 0);
 383
 384     FixAxis(results, 1/(1-yd)/1.7, "#frac{1}{N} #frac{dN_{ch}}{d#eta}");
 385
 386     p1->Clear();
 387     results->DrawClone("nostack e1");
 388
 389     fRangeParam->fSlave1Axis = results->GetXaxis();
 390     fRangeParam->fSlave1Pad  = p1;
 391
 392     // Draw other data
 393     if (others) {
 394       TGraphAsymmErrors* o      = 0;
 395       TIter              nextG(others->GetListOfGraphs());
 396       while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
 397         o->DrawClone("same p");
 398     }
 399
 400     // Make a legend in the main result pad
 401     TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
 402     l->SetNColumns(2);
 403     l->SetFillColor(0);
 404     l->SetFillStyle(0);
 405     l->SetBorderSize(0);
 406     l->SetTextFont(132);
 407
 408     // Put a title on top
 409     TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
 410     tit->SetNDC();
 411     tit->SetTextFont(132);
 412     tit->SetTextSize(0.05);
 413     tit->Draw();
 414
 415     // Put a nice label in the plot
 416     TString     eS;
 417     UShort_t    snn = (fSNNString ? fSNNString->GetUniqueID() : 0);
 418     const char* sys = (fSysString ? fSysString->GetTitle() : "?");
 419     if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
 420     else            eS = Form("%3dGeV", snn);
 421     TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s}=%s, %s",
 422                                            sys,
 423                                            eS.Data(),
 424                                            fTrigString ?
 425                                            fTrigString->GetTitle() : "?"));
 426     tt->SetNDC();
 427     tt->SetTextFont(132);
 428     tt->SetTextAlign(33);
 429     tt->Draw();
 430
 431     // Put number of accepted events on the plot
 432     Int_t nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
 433     TLatex* et = new TLatex(.93, .83, Form("%d events", nev));
 434     et->SetNDC();
 435     et->SetTextFont(132);
 436     et->SetTextAlign(33);
 437     et->Draw();
 438
 439     // Put number of accepted events on the plot
 440     if (fVtxAxis) {
 441       TLatex* vt = new TLatex(.93, .88, fVtxAxis->GetTitle());
 442       vt->SetNDC();
 443       vt->SetTextFont(132);
 444       vt->SetTextAlign(33);
 445       vt->Draw();
 446     }
 447     // results->Draw("nostack e1 same");
 448
 449     fRangeParam->fSlave1Axis = FindXAxis(p1, results->GetName());
 450     fRangeParam->fSlave1Pad  = p1;
 451
 452
 453     // Mark the plot as preliminary
 454     TLatex* pt = new TLatex(.12, .93, "Preliminary");
 455     pt->SetNDC();
 456     pt->SetTextFont(22);
 457     pt->SetTextSize(0.07);
 458     pt->SetTextColor(kRed+1);
 459     pt->SetTextAlign(13);
 460     pt->Draw();
 461
 462     if (!gSystem->AccessPathName("ALICE.png")) {
 463       TPad* logo = new TPad("logo", "logo", .12, .7, .32, .90, 0, 0, 0);
 464       logo->SetFillStyle(0);
 465       logo->Draw();
 466       logo->cd();
 467       TImage* i = TImage::Create();
 468       i->ReadImage("ALICE.png");
 469       i->Draw();
 470     }
 471     p1->cd();
 472   }
 473   //__________________________________________________________________
 474   void PlotRatios(THStack* ratios, Double_t max, Double_t y1, Double_t y2)
 475   {
 476     if (!ratios) return;
 477     bool isBottom = (y1 < 0.0001);
 478     Double_t yd = y2 - y1;
 479     // Make a sub-pad for the result itself
 480     TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0);
 481     p2->SetTopMargin(0.001);
 482     p2->SetRightMargin(0.05);
 483     p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
 484     p2->SetGridx();
 485     p2->SetTicks(1,1);
 486     p2->SetNumber(2);
 487     p2->Draw();
 488     p2->cd();
 489
 490     // Fix up axis
 491     FixAxis(ratios, 1/yd/1.7, "Ratios", 7);
 492
 493     ratios->SetMaximum(1+TMath::Max(.22,1.05*max));
 494     ratios->SetMinimum(1-TMath::Max(.32,1.05*max));
 495     p2->Clear();
 496     ratios->DrawClone("nostack e1");
 497
 498
 499     // Make a legend
 500     TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,
 501                                   isBottom ? .6 : .4);
 502     l2->SetNColumns(2);
 503     l2->SetFillColor(0);
 504     l2->SetFillStyle(0);
 505     l2->SetBorderSize(0);
 506     l2->SetTextFont(132);
 507
 508     // Make a nice band from 0.9 to 1.1
 509     TGraphErrors* band = new TGraphErrors(2);
 510     band->SetPoint(0, fForwardSym->GetXaxis()->GetXmin(), 1);
 511     band->SetPoint(1, fForward->GetXaxis()->GetXmax(), 1);
 512     band->SetPointError(0, 0, .1);
 513     band->SetPointError(1, 0, .1);
 514     band->SetFillColor(kYellow+2);
 515     band->SetFillStyle(3002);
 516     band->SetLineStyle(2);
 517     band->SetLineWidth(1);
 518     band->Draw("3 same");
 519     band->DrawClone("X L same");
 520
 521     // Replot the ratios on top
 522     ratios->DrawClone("nostack e1 same");
 523
 524     if (isBottom) {
 525       fRangeParam->fMasterAxis = FindXAxis(p2, ratios->GetName());
 526       p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
 527                                 fRangeParam));
 528     }
 529     else {
 530       fRangeParam->fSlave2Axis = FindXAxis(p2, ratios->GetName());
 531       fRangeParam->fSlave2Pad  = p2;
 532     }
 533   }
 534   //__________________________________________________________________
 535   void PlotLeftRight(THStack* leftright, Double_t max,
 536                      Double_t y1, Double_t y2)
 537   {
 538     if (!leftright) return;
 539     bool isBottom = (y1 < 0.0001);
 540     Double_t yd = y2 - y1;
 541     // Make a sub-pad for the result itself
 542     TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0);
 543     p3->SetTopMargin(0.001);
 544     p3->SetRightMargin(0.05);
 545     p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
 546     p3->SetGridx();
 547     p3->SetTicks(1,1);
 548     p3->SetNumber(2);
 549     p3->Draw();
 550     p3->cd();
 551
 552     TH1* dummy = 0;
 553     if (leftright->GetHists()->GetEntries() == 1) {
 554       // Add dummy histogram
 555       dummy = new TH1F("dummy","", 10, -6, 6);
 556       dummy->SetLineColor(0);
 557       dummy->SetFillColor(0);
 558       dummy->SetMarkerColor(0);
 559       leftright->Add(dummy);
 560     }
 561
 562     // Fix up axis
 563     FixAxis(leftright, 1/yd/1.7, "Right/Left", 4);
 564
 565     leftright->SetMaximum(1+TMath::Max(.12,1.05*max));
 566     leftright->SetMinimum(1-TMath::Max(.15,1.05*max));
 567     p3->Clear();
 568     leftright->DrawClone("nostack e1");
 569
 570
 571     // Make a legend
 572     TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
 573     l2->SetNColumns(2);
 574     l2->SetFillColor(0);
 575     l2->SetFillStyle(0);
 576     l2->SetBorderSize(0);
 577     l2->SetTextFont(132);
 578 #ifndef __CINT__
 579     if (dummy) {
 580       TList* prims = l2->GetListOfPrimitives();
 581       TIter next(prims);
 582       TLegendEntry* o = 0;
 583       while ((o = static_cast<TLegendEntry*>(next()))) {
 584         TString lbl(o->GetLabel());
 585         if (lbl != "dummy") continue;
 586         prims->Remove(o);
 587         break;
 588       }
 589     }
 590 #endif
 591     // Make a nice band from 0.9 to 1.1
 592     TGraphErrors* band = new TGraphErrors(2);
 593     band->SetPoint(0, fForwardSym->GetXaxis()->GetXmin(), 1);
 594     band->SetPoint(1, fForward->GetXaxis()->GetXmax(), 1);
 595     band->SetPointError(0, 0, .05);
 596     band->SetPointError(1, 0, .05);
 597     band->SetFillColor(kYellow+2);
 598     band->SetFillStyle(3002);
 599     band->SetLineStyle(2);
 600     band->SetLineWidth(1);
 601     band->Draw("3 same");
 602     band->DrawClone("X L same");
 603
 604     leftright->DrawClone("nostack e1 same");
 605     if (isBottom) {
 606       fRangeParam->fMasterAxis = FindXAxis(p3, leftright->GetName());
 607       p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
 608                                 fRangeParam));
 609     }
 610     else {
 611       fRangeParam->fSlave2Axis = FindXAxis(p3, leftright->GetName());
 612       fRangeParam->fSlave2Pad  = p3;
 613     }
 614   }
 615
 616   //__________________________________________________________________
 617   /**
 618    * Fix the apperance of the axis in a stack
 619    *
 620    * @param stack  stack of histogram
 621    * @param s      Scaling factor
 622    * @param ytitle Y axis title
 623    * @param force  Whether to draw the stack first or not
 624    * @param ynDiv  Divisions on Y axis
 625    */
 626   void FixAxis(THStack* stack, Double_t s, const char* ytitle,
 627                Int_t ynDiv=210, Bool_t force=true)
 628   {
 629     if (force) stack->Draw("nostack e1");
 630
 631     TH1* h = stack->GetHistogram();
 632     if (!h) return;
 633
 634     h->SetXTitle("#eta");
 635     h->SetYTitle(ytitle);
 636     TAxis* xa = h->GetXaxis();
 637     TAxis* ya = h->GetYaxis();
 638     if (xa) {
 639       xa->SetTitle("#eta");
 640       // xa->SetTicks("+-");
 641       xa->SetTitleSize(s*xa->GetTitleSize());
 642       xa->SetLabelSize(s*xa->GetLabelSize());
 643       xa->SetTickLength(s*xa->GetTickLength());
 644     }
 645     if (ya) {
 646       ya->SetTitle(ytitle);
 647       ya->SetDecimals();
 648       // ya->SetTicks("+-");
 649       ya->SetNdivisions(ynDiv);
 650       ya->SetTitleSize(s*ya->GetTitleSize());
 651       ya->SetTitleOffset(ya->GetTitleOffset()/s);
 652       ya->SetLabelSize(s*ya->GetLabelSize());
 653     }
 654   }
 655
 656   //__________________________________________________________________
 657   Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option) const
 658   {
 659     // Check if we have input
 660     if (!hist) return 0;
 661
 662     // Rebin if needed
 663     Rebin(hist);
 664
 665     stack->Add(hist, option);
 666     return hist->GetMaximum();
 667   }
 668   //__________________________________________________________________
 669   Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option,
 670                         TH1*& sym) const
 671   {
 672     // Check if we have input
 673     if (!hist) return 0;
 674
 675     // Rebin if needed
 676     Rebin(hist);
 677     stack->Add(hist, option);
 678
 679     // Now symmetrice the histogram
 680     sym = Symmetrice(hist);
 681     stack->Add(sym, option);
 682
 683     return hist->GetMaximum();
 684   }
 685
 686   //__________________________________________________________________
 687   /**
 688    * Rebin a histogram
 689    *
 690    * @param h     Histogram to rebin
 691    * @param rebin Rebinning factor
 692    *
 693    * @return
 694    */
 695   virtual void Rebin(TH1* h) const
 696   {
 697     if (fRebin <= 1) return;
 698
 699     Int_t nBins = h->GetNbinsX();
 700     if(nBins % fRebin != 0) {
 701       Warning("Rebin", "Rebin factor %d is not a devisor of current number "
 702               "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
 703       return;
 704     }
 705
 706     // Make a copy
 707     TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
 708     tmp->Rebin(fRebin);
 709     tmp->SetDirectory(0);
 710
 711     // The new number of bins
 712     Int_t nBinsNew = nBins / fRebin;
 713     for(Int_t i = 1;i<= nBinsNew; i++) {
 714       Double_t content = 0;
 715       Double_t sumw    = 0;
 716       Double_t wsum    = 0;
 717       Int_t    nbins   = 0;
 718       for(Int_t j = 1; j<=fRebin;j++) {
 719         Int_t    bin = (i-1)*fRebin + j;
 720         Double_t c   =  h->GetBinContent(bin);
 721
 722         if (c <= 0) continue;
 723
 724         if (fCutEdges) {
 725           if (h->GetBinContent(bin+1)<=0 ||
 726               h->GetBinContent(bin-1)) {
 727             Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
 728                     bin, c, h->GetName(),
 729                     bin+1, h->GetBinContent(bin+1),
 730                     bin-1, h->GetBinContent(bin-1));
 731             continue;
 732           }
 733         }
 734         Double_t e =  h->GetBinError(bin);
 735         Double_t w =  1 / (e*e); // 1/c/c
 736         content    += c;
 737         sumw       += w;
 738         wsum       += w * c;
 739         nbins++;
 740       }
 741
 742       if(content > 0 && nbins > 0) {
 743         tmp->SetBinContent(i, wsum / sumw);
 744         tmp->SetBinError(i,1./TMath::Sqrt(sumw));
 745       }
 746     }
 747
 748     // Finally, rebin the histogram, and set new content
 749     h->Rebin(fRebin);
 750     h->Reset();
 751     for(Int_t i = 1; i<= nBinsNew; i++) {
 752       h->SetBinContent(i,tmp->GetBinContent(i));
 753       h->SetBinError(i,  tmp->GetBinError(i));
 754     }
 755
 756     delete tmp;
 757   }
 758   //__________________________________________________________________
 759   /**
 760    * Make an extension of @a h to make it symmetric about 0
 761    *
 762    * @param h Histogram to symmertrice
 763    *
 764    * @return Symmetric extension of @a h
 765    */
 766   TH1* Symmetrice(const TH1* h) const
 767   {
 768     Int_t nBins = h->GetNbinsX();
 769     TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
 770     s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
 771     s->Reset();
 772     s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
 773     s->SetMarkerColor(h->GetMarkerColor());
 774     s->SetMarkerSize(h->GetMarkerSize());
 775     s->SetMarkerStyle(h->GetMarkerStyle()+4);
 776     s->SetFillColor(h->GetFillColor());
 777     s->SetFillStyle(h->GetFillStyle());
 778     s->SetDirectory(0);
 779
 780     // Find the first and last bin with data
 781     Int_t first = nBins+1;
 782     Int_t last  = 0;
 783     for (Int_t i = 1; i <= nBins; i++) {
 784       if (h->GetBinContent(i) <= 0) continue;
 785       first = TMath::Min(first, i);
 786       last  = TMath::Max(last,  i);
 787     }
 788
 789     Double_t xfirst = h->GetBinCenter(first-1);
 790     Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
 791     Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
 792     for (Int_t i = f1, j=l2; i <= last; i++,j--) {
 793       s->SetBinContent(j, h->GetBinContent(i));
 794       s->SetBinError(j, h->GetBinError(i));
 795     }
 796     // Fill in overlap bin
 797     s->SetBinContent(l2+1, h->GetBinContent(first));
 798     s->SetBinError(l2+1, h->GetBinError(first));
 799     return s;
 800   }
 801   //__________________________________________________________________
 802   Double_t RatioMax(TH1* h) const
 803   {
 804     Int_t    nBins = h->GetNbinsX();
 805     Double_t ret   = 0;
 806     for (Int_t i = 1; i <= nBins; i++) {
 807       Double_t c = h->GetBinContent(i);
 808       if (c == 0) continue;
 809       Double_t e = h->GetBinError(i);
 810       Double_t d = TMath::Abs(1-c-e);
 811       ret        = TMath::Max(d, ret);
 812     }
 813     return ret;
 814   }
 815   //__________________________________________________________________
 816   /**
 817    * Compute the ratio of @a h to @a g.  @a g is evaluated at the bin
 818    * centers of @a h
 819    *
 820    * @param h  Numerator
 821    * @param g  Divisor
 822    *
 823    * @return h/g
 824    */
 825   TH1* Ratio(const TH1* h, const TGraph* g, Double_t& max) const
 826   {
 827     if (!h || !g) return 0;
 828
 829     TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
 830     ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
 831     ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
 832     ret->Reset();
 833     ret->SetMarkerStyle(g->GetMarkerStyle());
 834     ret->SetMarkerColor(h->GetMarkerColor());
 835     ret->SetMarkerSize(0.9*g->GetMarkerSize());
 836     Double_t xlow  = g->GetX()[0];
 837     Double_t xhigh = g->GetX()[g->GetN()-1];
 838     if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
 839
 840     for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
 841       Double_t c = h->GetBinContent(i);
 842       if (c <= 0) continue;
 843
 844       Double_t x = h->GetBinCenter(i);
 845       if (x < xlow || x > xhigh) continue;
 846
 847       Double_t f = g->Eval(x);
 848       if (f <= 0) continue;
 849
 850       ret->SetBinContent(i, c / f);
 851       ret->SetBinError(i, h->GetBinError(i) / f);
 852     }
 853     if (ret->GetEntries() <= 0) {
 854       delete ret;
 855       ret = 0;
 856     }
 857     else
 858       max = TMath::Max(RatioMax(ret), max);
 859     return ret;
 860   }
 861   //__________________________________________________________________
 862   /**
 863    * Make the ratio of h1 to h2
 864    *
 865    * @param h1 First histogram (numerator)
 866    * @param h2 Second histogram (denominator)
 867    *
 868    * @return h1 / h2
 869    */
 870   TH1* Ratio(const TH1* h1, const TH1* h2, Double_t& max) const
 871   {
 872     if (!h1 || !h2) return 0;
 873     TH1* t1 = static_cast<TH1*>(h1->Clone(Form("%s_%s",
 874                                                h1->GetName(),
 875                                                h2->GetName())));
 876     t1->SetTitle(Form("%s / %s", h1->GetTitle(), h2->GetTitle()));
 877     t1->Divide(h2);
 878     t1->SetMarkerColor(h1->GetMarkerColor());
 879     t1->SetMarkerStyle(h2->GetMarkerStyle());
 880     max = TMath::Max(RatioMax(t1), max);
 881     return t1;
 882   }
 883   //__________________________________________________________________
 884   /**
 885    * Calculate the ratio of two graphs - g1 / g2
 886    *
 887    * @param g1 Numerator
 888    * @param g2 Denominator
 889    *
 890    * @return g1 / g2 in a histogram
 891    */
 892   TH1* Ratio(const TGraphAsymmErrors* g1,
 893              const TGraphAsymmErrors* g2, Double_t& max) const
 894   {
 895     Int_t    nBins = g1->GetN();
 896     TArrayF  bins(nBins+1);
 897     Double_t dx   = 0;
 898     for (Int_t i = 0; i < nBins; i++) {
 899       Double_t x   = g1->GetX()[i];
 900       Double_t exl = g1->GetEXlow()[i];
 901       Double_t exh = g1->GetEXhigh()[i];
 902       bins.fArray[i]   = x-exl;
 903       bins.fArray[i+1] = x+exh;
 904       Double_t dxi = exh+exl;
 905       if (i == 0) dx  = dxi;
 906       else if (dxi != dx) dx = 0;
 907     }
 908     TString name(Form("%s_%s", g1->GetName(), g2->GetName()));
 909     TString title(Form("%s / %s", g1->GetTitle(), g2->GetTitle()));
 910     TH1* h = 0;
 911     if (dx != 0) {
 912       h = new TH1F(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
 913     }
 914     else {
 915       h = new TH1F(name.Data(), title.Data(), nBins, bins.fArray);
 916     }
 917     h->SetMarkerStyle(g2->GetMarkerStyle());
 918     h->SetMarkerColor(g1->GetMarkerColor());
 919     h->SetMarkerSize(0.9*g2->GetMarkerSize());
 920
 921     Double_t low  = g2->GetX()[0];
 922     Double_t high = g2->GetX()[g2->GetN()-1];
 923     if (low > high) { Double_t t = low; low = high; high = t; }
 924     for (Int_t i = 0; i < nBins; i++) {
 925       Double_t x  = g1->GetX()[i];
 926       if (x < low-dx || x > high+dx) continue;
 927       Double_t c1 = g1->GetY()[i];
 928       Double_t e1 = g1->GetErrorY(i);
 929       Double_t c2 = g2->Eval(x);
 930
 931       h->SetBinContent(i+1, c1 / c2);
 932       h->SetBinError(i+1, e1 / c2);
 933     }
 934     max = TMath::Max(RatioMax(h), max);
 935     return h;
 936   }
 937
 938   //__________________________________________________________________
 939   TH1* Graph2Hist(const TGraphAsymmErrors* g) const
 940   {
 941     Int_t    nBins = g->GetN();
 942     TArrayF  bins(nBins+1);
 943     Double_t dx = 0;
 944     for (Int_t i = 0; i < nBins; i++) {
 945       Double_t x   = g->GetX()[i];
 946       Double_t exl = g->GetEXlow()[i];
 947       Double_t exh = g->GetEXhigh()[i];
 948       bins.fArray[i]   = x-exl;
 949       bins.fArray[i+1] = x+exh;
 950       Double_t dxi = exh+exl;
 951       if (i == 0) dx  = dxi;
 952       else if (dxi != dx) dx = 0;
 953     }
 954     TString name(g->GetName());
 955     TString title(g->GetTitle());
 956     TH1D* h = 0;
 957     if (dx != 0) {
 958       h = new TH1D(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
 959     }
 960     else {
 961       h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
 962     }
 963     h->SetMarkerStyle(g->GetMarkerStyle());
 964     h->SetMarkerColor(g->GetMarkerColor());
 965     h->SetMarkerSize(g->GetMarkerSize());
 966
 967     return h;
 968   }
 969   //__________________________________________________________________
 970   TH1* Asymmetry(TH1* h, Double_t& max)
 971   {
 972     if (!h) return 0;
 973
 974     TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName())));
 975     // Int_t    oBins = h->GetNbinsX();
 976     // Double_t high  = h->GetXaxis()->GetXmax();
 977     // Double_t low   = h->GetXaxis()->GetXmin();
 978     // Double_t dBin  = (high - low) / oBins;
 979     // Int_t    tBins = Int_t(2*high/dBin+.5);
 980     // ret->SetBins(tBins, -high, high);
 981     ret->Reset();
 982     ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
 983     ret->SetYTitle("Right/Left");
 984     Int_t nBins = h->GetNbinsX();
 985     for (Int_t i = 1; i <= nBins; i++) {
 986       Double_t x = h->GetBinCenter(i);
 987       if (x > 0) break;
 988
 989       Double_t c1 = h->GetBinContent(i);
 990       Double_t e1 = h->GetBinError(i);
 991       if (c1 <= 0) continue;
 992
 993       Int_t    j  = h->FindBin(-x);
 994       if (j <= 0 || j > nBins) continue;
 995
 996       Double_t c2 = h->GetBinContent(j);
 997       Double_t e2 = h->GetBinError(j);
 998
 999       Double_t c12 = c1*c1;
1000       Double_t e   = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12));
1001
1002       Int_t    k   = ret->FindBin(x);
1003       ret->SetBinContent(k, c2/c1);
1004       ret->SetBinError(k, e);
1005     }
1006     max = TMath::Max(max, RatioMax(ret));
1007
1008     return ret;
1009   }
1010
1011
1012   //__________________________________________________________________
1013   Bool_t      fShowOthers;
1014   Bool_t      fShowAlice;
1015   Bool_t      fShowRatios;
1016   Bool_t      fShowLeftRight;
1017   UShort_t    fRebin;
1018   Bool_t      fCutEdges;
1019   TString     fTitle;
1020   TString     fHHDFile;
1021   TNamed*     fTrigString;
1022   TNamed*     fSNNString;
1023   TNamed*     fSysString;
1024   TAxis*      fVtxAxis;
1025   TH1*        fForward;
1026   TH1*        fForwardMC;
1027   TH1*        fForwardHHD;
1028   TH1*        fTruth;
1029   TH1*        fCentral;
1030   TH1*        fForwardSym;
1031   TH1*        fForwardMCSym;
1032   TH1*        fForwardHHDSym;
1033   TH1*        fTriggers;
1034   RangeParam* fRangeParam;
1035
1036 };
1037
1038 //=== Stuff for auto zooming =========================================
1039 void UpdateRange(dNdetaDrawer::RangeParam* p)
1040 {
1041   if (!p) {
1042     Warning("UpdateRange", "No parameters %p", p);
1043     return;
1044   }
1045   if (!p->fMasterAxis) {
1046     Warning("UpdateRange", "No master axis %p", p->fMasterAxis);
1047     return;
1048   }
1049   Int_t    first = p->fMasterAxis->GetFirst();
1050   Int_t    last  = p->fMasterAxis->GetLast();
1051   Double_t x1    = p->fMasterAxis->GetBinCenter(first);
1052   Double_t x2    = p->fMasterAxis->GetBinCenter(last);
1053   //Info("UpdateRange", "Range set to [%3d,%3d]->[%f,%f]", first, last, x1,x2);
1054
1055   if (p->fSlave1Axis) {
1056     Int_t i1 = p->fSlave1Axis->FindBin(x1);
1057     Int_t i2 = p->fSlave1Axis->FindBin(x2);
1058     p->fSlave1Axis->SetRange(i1, i2);
1059     p->fSlave1Pad->Modified();
1060     p->fSlave1Pad->Update();
1061   }
1062   if (p->fSlave2Axis) {
1063     Int_t i1 = p->fSlave2Axis->FindBin(x1);
1064     Int_t i2 = p->fSlave2Axis->FindBin(x2);
1065     p->fSlave2Axis->SetRange(i1, i2);
1066     p->fSlave2Pad->Modified();
1067     p->fSlave2Pad->Update();
1068   }
1069   TCanvas*  c = gPad->GetCanvas();
1070   c->cd();
1071 }
1072
1073 //____________________________________________________________________
1074 void RangeExec(dNdetaDrawer::RangeParam* p)
1075 {
1076   // Event types:
1077   //  51:   Mouse move
1078   //  53:
1079   //  1:    Button down
1080   //  21:   Mouse drag
1081   //  11:   Mouse release
1082   // dNdetaDrawer::RangeParam* p =
1083   //   reinterpret_cast<dNdetaDrawer::RangeParam*>(addr);
1084   Int_t event     = gPad->GetEvent();
1085   TObject *select = gPad->GetSelected();
1086   if (event == 53) {
1087     UpdateRange(p);
1088     return;
1089   }
1090   if (event != 11 || !select || select != p->fMasterAxis) return;
1091   UpdateRange(p);
1092 }
1093
1094 //=== Steering function ==============================================
1095 void
1096 DrawdNdeta(const char* filename="forward_dndeta.root",
1097            Int_t       flags=0xf,
1098            const char* title=""
1099            UShort_t    rebin=5,
1100            Bool_t      cutEdges=false)
1101 {
1102   dNdetaDrawer d;
1103   d.SetRebin(rebin);
1104   d.SetCutEdges(cutEdges);
1105   d.SetTitle(title);
1106   d.SetHHDFile("");
1107   d.SetShowOthers(flags & 0x1);
1108   d.SetShowAlice(flags & 0x2);
1109   d.SetShowRatios(flags & 0x4);
1110   d.SetShowLeftRight(flags & 0x8);
1111   d.Run(filename);
1112 }