PWG2/FORWARD/analysis2/AnalyseAOD.C

   1 /**
   2  * @file
   3  *
   4  * @ingroup pwg2_forward_scripts
   5  */
   6 #include <TH1D.h>
   7 #include <TH2D.h>
   8 #include <TCanvas.h>
   9 #include <TPad.h>
  10 #include <TFile.h>
  11 #include <TTree.h>
  12 #include <TError.h>
  13 #include <TStyle.h>
  14 #include <THStack.h>
  15 #include <TLegend.h>
  16 #include <TMath.h>
  17 #include <TMultiGraph.h>
  18 #include <TGraph.h>
  19 #include <TGraphErrors.h>
  20 #include <TLatex.h>
  21 #include <TSystem.h>
  22 #include "AliAODForwardMult.h"
  23 #include "OtherData.C"
  24
  25 /**
  26  * Draw the data stored in the AOD
  27  *
  28  * To use this, do
  29  *
  30  * @code
  31  *   Root> .L $ALICE_ROOT/PWG2/FORWARD/analysis2/Compile.C
  32  *   Root> Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/AnalyseAOD.C")
  33  *   Root> AnalyseAOD dr
  34  *   Root> dr.Run("AliAODs.root",-10,10,5,AliAODForwardMult::kInel,900)
  35  * @endcode
  36  *
  37  * The output is stored in a ROOT file
  38  *
  39  * See also the script Pass2.C
  40  *
  41  * @ingroup pwg2_forward_scripts
  42  */
  43 struct AnalyseAOD
  44 {
  45 public:
  46   /** AOD tree */
  47   TTree*             fTree;
  48   /** AOD object */
  49   AliAODForwardMult* fAOD;
  50   /** AOD object */
  51   AliAODForwardMult* fMCAOD;
  52   /** Output file */
  53   TFile*             fOut;
  54   /** Summed histogram */
  55   TH2D*              fSum;
  56   /** Summed histogram */
  57   TH2D*              fMCSum;
  58   /** Primary information */
  59   TH2D*              fPrimary;
  60   /** Sum primary information */
  61   TH2D*              fSumPrimary;
  62   /** Vertex efficiency */
  63   Double_t           fVtxEff;
  64   /** Title to put on the plot */
  65   TString            fTitle;
  66   /** Do HHD comparison */
  67   Bool_t             fDoHHD;
  68   /** Number of events with a trigger */
  69   Int_t              fNTriggered;
  70   /** Number of events with a vertex */
  71   Int_t              fNWithVertex;
  72   /** Number of events accepted */
  73   Int_t              fNAccepted;
  74   /** Number of events accepted */
  75   Int_t              fNAll;
  76
  77   //__________________________________________________________________
  78   /**
  79    * Constructor
  80    *
  81    */
  82   AnalyseAOD()
  83     : fTree(0),
  84       fAOD(0),
  85       fMCAOD(0),
  86       fOut(0),
  87       fSum(0),
  88       fMCSum(0),
  89       fPrimary(0),
  90       fSumPrimary(0),
  91       fVtxEff(0),
  92       fTitle(""),
  93       fDoHHD(kTRUE)
  94   {}
  95   //__________________________________________________________________
  96   /**
  97    * Reset internal structures
  98    *
  99    */
 100   void Clear(Option_t* )
 101   {
 102     if (fTree && fTree->GetCurrentFile()) {
 103       fTree->GetCurrentFile()->Close();
 104       delete fTree;
 105     }
 106     if (fOut) {
 107       fOut->Close();
 108       delete fOut;
 109     }
 110     if (fSum)        delete fSum;
 111     if (fMCSum)      delete fMCSum;
 112     if (fSumPrimary) delete fSumPrimary;
 113     fTree       = 0;
 114     fOut        = 0;
 115     fSum        = 0;
 116     fMCSum      = 0;
 117     fSumPrimary = 0;
 118     fVtxEff     = 0;
 119
 120   }
 121   //__________________________________________________________________
 122   /**
 123    * Run
 124    *
 125    * @param file   Input file with AOD tree
 126    * @param vzMin  Minimum interaction point z coordinate
 127    * @param vzMax  Maximum interaction point z coordinate
 128    * @param rebin  How many times to re-bin the @f$ dN_{ch}/d\eta@f$
 129    * @param mask   Trigger mask
 130    * @param energy Collision energy @f$ \sqrt{s}@f$ or @f$\sqrt{s_{NN}}@f$
 131    * @param title  Title to put on the plot
 132    * @param doHHD  Whether to do HHD comparison
 133    * @param doComp Whether to do comparisons
 134    *
 135    * @return True on success, false otherwise
 136    */
 137   Bool_t Run(const char* file="AliAODs.root",
 138              Double_t vzMin=-10, Double_t vzMax=10, Int_t rebin=1,
 139              Int_t mask=AliAODForwardMult::kInel, Int_t energy=900,
 140              const char* title="", Bool_t doHHD=false, Bool_t doComp=false)
 141   {
 142     TString trgName;
 143     if    (mask & AliAODForwardMult::kInel)    trgName.Append("inel-");
 144     if    (mask & AliAODForwardMult::kInelGt0) trgName.Append("inelgt0-");
 145     if    (mask & AliAODForwardMult::kNSD)     trgName.Append("nsd-");
 146     if (trgName.EndsWith("-")) trgName.Remove(trgName.Length()-1);
 147     if (trgName.IsNull()) trgName = "unknown";
 148     TString outName =
 149       TString::Format("dndeta_%04dGeV_%c%02d-%c%02dcm_rb%02d_%s.root",
 150                       energy,
 151                       vzMin < 0 ? 'm' : 'p', Int_t(TMath::Abs(vzMin)+.5),
 152                       vzMax < 0 ? 'm' : 'p', Int_t(TMath::Abs(vzMax)+.5),
 153                       rebin, trgName.Data());
 154     fTitle = title;
 155     if (!Open(file, outName)) return kFALSE;
 156     if (!Process(vzMin,vzMax,mask)) return kFALSE;
 157     if (!Finish(rebin, mask, energy,doHHD,doComp)) return kFALSE;
 158
 159     return kTRUE;
 160   }
 161   //__________________________________________________________________
 162   /**
 163    * Open the files @a fname containg the tree with AliAODEvent objects,
 164    * and also open the output file @a outname for writting.
 165    *
 166    * @param fname   Input file name
 167    * @param outname Output file name
 168    *
 169    * @return true on success, false otherwise
 170    */
 171   Bool_t Open(const char* fname, const char* outname)
 172   {
 173     Clear("");
 174
 175     TFile* file = TFile::Open(fname, "READ");
 176     if (!file) {
 177       Error("Open", "Cannot open file %s", fname);
 178       return kFALSE;
 179     }
 180
 181     // Get the AOD tree
 182     fTree = static_cast<TTree*>(file->Get("aodTree"));
 183     if (!fTree) {
 184       Error("Init", "Couldn't get the tree");
 185       return kFALSE;
 186     }
 187
 188     // Set the branch pointer
 189     fTree->SetBranchAddress("Forward", &fAOD);
 190
 191     // Set the branch pointer
 192     fTree->SetBranchAddress("ForwardMC", &fMCAOD);
 193
 194     // Set the branch pointer
 195     fTree->SetBranchAddress("primary", &fPrimary);
 196
 197     fOut = TFile::Open(outname, "RECREATE");
 198     if (!fOut) {
 199       Error("Open", "Couldn't open %s", outname);
 200       return kFALSE;
 201     }
 202     return kTRUE;
 203   }
 204   //__________________________________________________________________
 205   /**
 206    * Process the events
 207    *
 208    * @param vzMin  Minimum interaction point z coordinate
 209    * @param vzMax  Maximum interaction point z coordinate
 210    * @param mask   Trigger mask
 211    *
 212    * @return true on success, false otherwise
 213    */
 214   Bool_t Process(Double_t vzMin, Double_t vzMax, Int_t mask)
 215   {
 216     fNTriggered       = 0;                    // # of triggered ev.
 217     fNWithVertex      = 0;                    // # of ev. w/vertex
 218     fNAccepted        = 0;                    // # of ev. used
 219     Int_t nAvailable  = fTree->GetEntries();  // How many entries
 220
 221     for (int i = 0; i < nAvailable; i++) {
 222       fTree->GetEntry(i);
 223       if (((i+1) % 100) == 0) {
 224         fprintf(stdout,"Event # %9d of %9d, %9d accepted so far\r",
 225                i+1, nAvailable, fNAccepted);
 226         fflush(stdout);
 227       }
 228
 229       // Create sum histogram on first event - to match binning to input
 230       if (!fSum) {
 231         fSum = static_cast<TH2D*>(fAOD->GetHistogram().Clone("d2ndetadphi"));
 232         Info("Process", "Created sum histogram (%d,%f,%f)x(%d,%f,%f)",
 233              fSum->GetNbinsX(),
 234              fSum->GetXaxis()->GetXmin(),
 235              fSum->GetXaxis()->GetXmax(),
 236              fSum->GetNbinsY(),
 237              fSum->GetYaxis()->GetXmin(),
 238              fSum->GetYaxis()->GetXmax());
 239       }
 240       if (!fMCSum && fTree->GetBranch("ForwardMC")) {
 241         fMCSum =
 242           static_cast<TH2D*>(fMCAOD->GetHistogram().Clone("d2ndetadphiMC"));
 243         Info("Process", "Created MC sum histogram (%d,%f,%f)x(%d,%f,%f)",
 244              fMCSum->GetNbinsX(),
 245              fMCSum->GetXaxis()->GetXmin(),
 246              fMCSum->GetXaxis()->GetXmax(),
 247              fMCSum->GetNbinsY(),
 248              fMCSum->GetYaxis()->GetXmin(),
 249              fMCSum->GetYaxis()->GetXmax());
 250       }
 251       if (!fSumPrimary && fTree->GetBranch("primary")) {
 252         fSumPrimary =
 253           static_cast<TH2D*>(fPrimary->Clone("primarySum"));
 254         Info("Process", "Created MC truth histogram (%d,%f,%f)x(%d,%f,%f)",
 255              fMCSum->GetNbinsX(),
 256              fMCSum->GetXaxis()->GetXmin(),
 257              fMCSum->GetXaxis()->GetXmax(),
 258              fMCSum->GetNbinsY(),
 259              fMCSum->GetYaxis()->GetXmin(),
 260              fMCSum->GetYaxis()->GetXmax());
 261       }
 262
 263       // Add contribution from this event
 264       if (fSumPrimary) fSumPrimary->Add(fPrimary);
 265
 266       fNAll++;
 267
 268       // fAOD->Print();
 269       // Other trigger/event requirements could be defined
 270       if (!fAOD->IsTriggerBits(mask)) continue;
 271       fNTriggered++;
 272
 273       // Check if there's a vertex
 274       if (!fAOD->HasIpZ()) continue;
 275       fNWithVertex++;
 276
 277       // Select vertex range (in centimeters)
 278       if (!fAOD->InRange(vzMin, vzMax)) continue;
 279       fNAccepted++;
 280
 281       // Add contribution from this event
 282       fSum->Add(&(fAOD->GetHistogram()));
 283
 284       // Add contribution from this event
 285       if (fMCSum) fMCSum->Add(&(fMCAOD->GetHistogram()));
 286
 287
 288     }
 289     printf("\n");
 290     fVtxEff = Double_t(fNWithVertex)/fNTriggered;
 291
 292     Info("Process", "Total of %9d events\n"
 293          "               %9d has a trigger\n"
 294          "               %9d has a vertex\n"
 295          "               %9d was used\n",
 296          nAvailable, fNTriggered, fNWithVertex, fNAccepted);
 297
 298     return kTRUE;
 299   }
 300   //__________________________________________________________________
 301   /**
 302    * Finish the stuff and draw
 303    *
 304    * @param rebin  How many times to rebin
 305    * @param energy Collision energy
 306    * @param mask   Trigger mask
 307    * @param doHHD  Whether to do HHD comparison
 308    * @param doComp Whether to do comparisons
 309    *
 310    * @return true on success, false otherwise
 311    */
 312   Bool_t Finish(Int_t rebin, Int_t mask, Int_t energy,
 313                 Bool_t doHHD, Bool_t doComp)
 314   {
 315     fOut->cd();
 316
 317     // Get acceptance normalisation from underflow bins
 318     TH1D* norm   = fSum->ProjectionX("norm", 0, 1, "");
 319     // Project onto eta axis - _ignoring_underflow_bins_!
 320     TH1D* dndeta = fSum->ProjectionX("dndeta", 1, -1, "e");
 321     dndeta->SetTitle("ALICE Forward");
 322     // Normalize to the acceptance
 323     dndeta->Divide(norm);
 324     // Scale by the vertex efficiency
 325     dndeta->Scale(fVtxEff, "width");
 326     dndeta->SetMarkerColor(kRed+1);
 327     dndeta->SetMarkerStyle(20);
 328     dndeta->SetMarkerSize(1);
 329     dndeta->SetFillStyle(0);
 330     Rebin(dndeta, rebin);
 331
 332     TH1D* dndetaMC = 0;
 333     if (fMCSum) {
 334       // Get acceptance normalisation from underflow bins
 335       norm = fMCSum->ProjectionX("norm", 0, 1, "");
 336       // Project onto eta axis - _ignoring_underflow_bins_!
 337       dndetaMC = fMCSum->ProjectionX("dndetaMC", 1, -1, "e");
 338       dndetaMC->SetTitle("ALICE Forward (MC)");
 339       // Normalize to the acceptance
 340       dndetaMC->Divide(norm);
 341       // Scale by the vertex efficiency
 342       dndetaMC->Scale(fVtxEff, "width");
 343       dndetaMC->SetMarkerColor(kRed+3);
 344       dndetaMC->SetMarkerStyle(21);
 345       dndetaMC->SetMarkerSize(1);
 346       dndetaMC->SetFillStyle(0);
 347       Rebin(dndetaMC, rebin);
 348     }
 349
 350     TH1D* dndetaTruth = 0;
 351     if (fSumPrimary) {
 352       dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth", -1, -1, "e");
 353       //dndetaTruth->Scale(1./fNTriggered, "width");
 354       dndetaTruth->Scale(1./fNAll, "width");
 355       dndetaTruth->SetMarkerColor(kGray+3);
 356       dndetaTruth->SetMarkerStyle(22);
 357       dndetaTruth->SetMarkerSize(1);
 358       dndetaTruth->SetFillStyle(0);
 359       Rebin(dndetaTruth, rebin);
 360     }
 361
 362     DrawIt(dndeta, dndetaMC, dndetaTruth, mask, energy, doHHD, doComp);
 363
 364     return kTRUE;
 365   }
 366   //__________________________________________________________________
 367   /**
 368    */
 369   void DrawIt(TH1* dndeta, TH1* dndetaMC, TH1* dndetaTruth,
 370               Int_t mask, Int_t energy, Bool_t doHHD,
 371               Bool_t doComp)
 372   {
 373     // --- 1st part - prepare data -----------------------------------
 374     TH1* dndetaSym = Symmetrice(dndeta);
 375
 376     Double_t max = dndeta->GetMaximum();
 377
 378     // Make our histogram stack
 379     THStack* stack = new THStack("results", "Results");
 380
 381     TH1* dndetaTruthSym = 0;
 382     if (dndetaTruth) {
 383       dndetaTruth->SetFillColor(kGray);
 384       dndetaTruth->SetFillStyle(3001);
 385       dndetaTruthSym = Symmetrice(dndetaTruth);
 386       stack->Add(dndetaTruthSym, "e5 p");
 387       stack->Add(dndetaTruth, "e5 p");
 388       Info("DrawIt", "Maximum of MC dndeta (truth): %f, was %f",
 389            dndetaTruth->GetMaximum(),dndeta->GetMaximum());
 390       max = TMath::Max(dndetaTruth->GetMaximum(),max);
 391     }
 392
 393     // Get the data from HHD's analysis - if available
 394     TH1* dndetaHHD    = 0;
 395     TH1* dndetaHHDSym = 0;
 396     // Info("DrawIt", "Will %sdraw HHD results", (doHHD ? "" : "not "));
 397     if (doHHD) {
 398       TString hhdName(fOut->GetName());
 399       hhdName.ReplaceAll("dndeta", "hhd");
 400       dndetaHHD    = GetHHD(hhdName.Data(), mask & AliAODForwardMult::kNSD);
 401       dndetaHHDSym = 0;
 402       if (dndetaHHD) {
 403         // Symmetrice and add to stack
 404         dndetaHHD->SetTitle("ALICE Forward (HHD)");
 405         dndetaHHDSym = Symmetrice(dndetaHHD);
 406         stack->Add(dndetaHHDSym);
 407         stack->Add(dndetaHHD);
 408         max = TMath::Max(dndetaHHD->GetMaximum(),max);
 409       }
 410       else
 411         Warning("DrawIt", "No HHD data found");
 412       fOut->cd();
 413     }
 414
 415     // If we have MC analysed data, plot it
 416     if (dndetaMC) {
 417       TH1* dndetaMCSym = Symmetrice(dndetaMC);
 418       stack->Add(dndetaMCSym);
 419       stack->Add(dndetaMC);
 420       max = TMath::Max(dndetaMC->GetMaximum(),max);
 421     }
 422
 423     // Add the analysis results to the list
 424     stack->Add(dndetaSym);
 425     stack->Add(dndeta);
 426
 427     // Get graph of 'other' data - e.g., UA5, CMS, ... - and check if
 428     // there's any graphs.  Set the pad division based on that.
 429     // Info("DrawIt", "Will %sdraw other results", (doComp ? "" : "not "));
 430     TMultiGraph* other    = (doComp ? GetData(energy, mask) : 0);
 431     THStack*     ratios   = MakeRatios(dndeta,      dndetaSym,
 432                                        dndetaHHD,   dndetaHHDSym,
 433                                        dndetaTruth, dndetaTruthSym,
 434                                        other);
 435
 436     // Check if we have ratios
 437
 438     // --- 2nd part - Draw in canvas ---------------------------------
 439     // Make a canvas
 440     gStyle->SetOptTitle(0);
 441     gStyle->SetTitleFont(132, "xyz");
 442     gStyle->SetLabelFont(132, "xyz");
 443     TCanvas* c = new TCanvas("c", "C", 900, 800);
 444     c->SetFillColor(0);
 445     c->SetBorderMode(0);
 446     c->SetBorderSize(0);
 447     c->cd();
 448
 449     Double_t yd = (ratios ? 0.3 : 0);
 450
 451     // Make a sub-pad for the result itself
 452     TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0);
 453     p1->SetTopMargin(0.05);
 454     p1->SetBottomMargin(ratios ? 0.001 : 0.1);
 455     p1->SetRightMargin(0.05);
 456     p1->SetGridx();
 457     p1->SetTicks(1,1);
 458     p1->Draw();
 459     p1->cd();
 460
 461     // Fix the apperance of the stack and redraw.
 462     if (other) {
 463       TGraphAsymmErrors* o      = 0;
 464       TIter              nextG(other->GetListOfGraphs());
 465       while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
 466         Double_t gmax = TMath::MaxElement(o->GetN(), o->GetY());
 467         //Info("DrawIt", "Maximum of %s is %f, was %f", o->GetName(),gmax,max);
 468         max = TMath::Max(max, gmax);
 469       }
 470     }
 471     max *= 1.1;
 472     // Info("DrawIt", "Setting maximum to %f", max);
 473     stack->SetMinimum(ratios ? -0.1 : 0);
 474     stack->SetMaximum(max);
 475     FixAxis(stack, 1/(1-yd)/1.5, "#frac{1}{N} #frac{dN_{ch}}{#eta}");
 476     p1->Clear();
 477     stack->DrawClone("nostack e1");
 478
 479     // Draw other data
 480     if (other) {
 481       TGraphAsymmErrors* o      = 0;
 482       TIter              nextG(other->GetListOfGraphs());
 483       while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
 484         o->Draw("same p");
 485     }
 486
 487     // Make a legend in the main result pad
 488     TString trigs(AliAODForwardMult::GetTriggerString(mask));
 489     TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
 490     l->SetNColumns(2);
 491     l->SetFillColor(0);
 492     l->SetFillStyle(0);
 493     l->SetBorderSize(0);
 494     l->SetTextFont(132);
 495     p1->cd();
 496
 497     // Put a title on top
 498     TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
 499     tit->SetNDC();
 500     tit->SetTextFont(132);
 501     tit->SetTextSize(0.05);
 502     tit->Draw();
 503
 504     // Put a nice label in the plot
 505     TLatex* tt = new TLatex(.93, .93,
 506                             Form("#sqrt{s}=%dGeV, %s", energy,
 507                                  AliAODForwardMult::GetTriggerString(mask)));
 508     tt->SetNDC();
 509     tt->SetTextFont(132);
 510     tt->SetTextAlign(33);
 511     tt->Draw();
 512
 513     // Mark the plot as preliminary
 514     TLatex* pt = new TLatex(.93, .88, "Preliminary");
 515     pt->SetNDC();
 516     pt->SetTextFont(22);
 517     pt->SetTextSize(0.07);
 518     pt->SetTextColor(kRed+1);
 519     pt->SetTextAlign(33);
 520     pt->Draw();
 521     c->cd();
 522
 523     // If we do not have the ratios, fix up the display
 524     // p1->SetPad(0, 0, 1, 1);
 525     // p1->SetBottomMargin(0.1);
 526     // l->SetY1(0.11);
 527     // stack->SetMinimum(0);
 528     // FixAxis(stack, (1-yd)/1,  "#frac{1}{N} #frac{dN_{ch}}{#eta}",10,false);
 529     if (ratios) {
 530       // If we do have the ratios, then make a new pad and draw the
 531       // ratios there
 532       c->cd();
 533       TPad* p2 = new TPad("p2", "p2", 0, 0.0, 1.0, yd, 0, 0);
 534       p2->SetTopMargin(0.001);
 535       p2->SetRightMargin(0.05);
 536       p2->SetBottomMargin(1/yd * 0.07);
 537       p2->SetGridx();
 538       p2->SetTicks(1,1);
 539       p2->Draw();
 540       p2->cd();
 541
 542       // Fix up axis
 543       FixAxis(ratios, 1/yd/1.5, "Ratios", 5);
 544
 545       // Fix up y range and redraw
 546       ratios->SetMinimum(.58);
 547       ratios->SetMaximum(1.22);
 548       p2->Clear();
 549       ratios->DrawClone("nostack e1");
 550
 551       // Make a legend
 552       TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,.6);
 553       l2->SetNColumns(2);
 554       l2->SetFillColor(0);
 555       l2->SetFillStyle(0);
 556       l2->SetBorderSize(0);
 557       l2->SetTextFont(132);
 558
 559       // Make a nice band from 0.9 to 1.1
 560       TGraphErrors* band = new TGraphErrors(2);
 561       band->SetPoint(0, dndetaSym->GetXaxis()->GetXmin(), 1);
 562       band->SetPoint(1, dndeta->GetXaxis()->GetXmax(), 1);
 563       band->SetPointError(0, 0, .1);
 564       band->SetPointError(1, 0, .1);
 565       band->SetFillColor(kYellow+2);
 566       band->SetFillStyle(3002);
 567       band->Draw("3 same");
 568
 569       // Replot the ratios on top
 570       ratios->DrawClone("nostack e1 same");
 571
 572       c->cd();
 573     }
 574
 575     // Plot to disk
 576     TString imgName(fOut->GetName());
 577     imgName.ReplaceAll(".root", ".png");
 578     c->SaveAs(imgName.Data());
 579     imgName.ReplaceAll(".png", ".C");
 580     c->SaveAs(imgName.Data());
 581
 582     fOut->cd();
 583     stack->Write();
 584     if (other)  other->Write();
 585     if (ratios) ratios->Write();
 586
 587     // Close our file
 588     fOut->Close();
 589   }
 590   //__________________________________________________________________
 591   /**
 592    * Get the result from previous analysis code
 593    *
 594    * @param fn  File to open
 595    * @param nsd Whether this is NSD
 596    *
 597    * @return null or result of previous analysis code
 598    */
 599   TH1* GetHHD(const char* fn="fmd_dNdeta_mult.root", bool nsd=false)
 600   {
 601     TDirectory* savdir = gDirectory;
 602     if (gSystem->AccessPathName(fn)) {
 603       Warning("GetHHD", "Output of HHD analysis (%s) not available", fn);
 604       return 0;
 605     }
 606     TFile* file = TFile::Open(fn, "READ");
 607     if (!file) {
 608       Warning("GetHHD", "couldn't open HHD file %s", fn);
 609       savdir->cd();
 610       return 0;
 611     }
 612     TString hist(Form("dNdeta_dNdeta%s", (nsd ? "NSD" : "")));
 613     TH1* h = static_cast<TH1*>(file->Get(hist.Data()));
 614     if (!h) {
 615       Warning("GetHHD", "Couldn't find HHD histogram %s in %s",
 616               hist.Data(), fn);
 617       file->Close();
 618       savdir->cd();
 619       return 0;
 620     }
 621     TH1* r = static_cast<TH1*>(h->Clone("dndeta_hhd"));
 622     r->SetTitle("1/N dN_{ch}/d#eta (HHD)");
 623     r->SetFillStyle(0);
 624     r->SetFillColor(0);
 625     r->SetMarkerStyle(21);
 626     r->SetMarkerColor(kPink+1);
 627     r->SetDirectory(savdir);
 628
 629     file->Close();
 630     savdir->cd();
 631     return r;
 632   }
 633   //__________________________________________________________________
 634   /**
 635    */
 636   THStack* MakeRatios(const TH1* dndeta, const TH1* sym,
 637                       const TH1* hhd,    const TH1* hhdsym,
 638                       const TH1* mc,     const TH1* mcsym,
 639                       TMultiGraph* other) const
 640   {
 641     // If we have 'other' data, then do the ratio of the results to that
 642     Bool_t hasOther = (other && other->GetListOfGraphs() &&
 643                        other->GetListOfGraphs()->GetEntries() > 0);
 644     Bool_t hasHhd   = (hhd && hhdsym);
 645     if (!hasOther && !hasHhd && !mc && !mcsym) return 0;
 646
 647     THStack* ratios = new THStack("ratios", "Ratios");
 648     if (hasOther) {
 649       TGraphAsymmErrors* o      = 0;
 650       TIter              nextG(other->GetListOfGraphs());
 651       while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
 652         ratios->Add(Ratio(dndeta, o));
 653         ratios->Add(Ratio(sym,    o));
 654         ratios->Add(Ratio(hhd,    o));
 655         ratios->Add(Ratio(hhdsym, o));
 656       }
 657     }
 658
 659     // If we have data from HHD's analysis, then do the ratio of
 660     // our result to that data.
 661     if (hasHhd) {
 662       TH1F* t1 = static_cast<TH1F*>(dndeta->Clone(Form("%s_%s",
 663                                                        dndeta->GetName(),
 664                                                        hhd->GetName())));
 665       TH1F* t2 = static_cast<TH1F*>(sym->Clone(Form("%s_%s",
 666                                                     sym->GetName(),
 667                                                     hhdsym->GetName())));
 668       t1->SetTitle(Form("%s / %s", dndeta->GetTitle(), hhd->GetTitle()));
 669       t2->SetTitle(Form("%s / %s", sym->GetTitle(), hhdsym->GetTitle()));
 670       t1->Divide(hhd);
 671       t2->Divide(hhdsym);
 672       t1->SetMarkerColor(hhd->GetMarkerColor());
 673       t2->SetMarkerColor(hhdsym->GetMarkerColor());
 674       ratios->Add(t1);
 675       ratios->Add(t2);
 676     }
 677
 678     // Do comparison to MC
 679     if (mc) {
 680       TH1D* t1 = static_cast<TH1D*>(dndeta->Clone(Form("%s_%s",
 681                                                        dndeta->GetName(),
 682                                                        mc->GetName())));
 683       TH1D* t2 = static_cast<TH1D*>(sym->Clone(Form("%s_%s",
 684                                                     sym->GetName(),
 685                                                     mcsym->GetName())));
 686       t1->SetTitle(Form("%s / %s", dndeta->GetTitle(), mc->GetTitle()));
 687       t2->SetTitle(Form("%s / %s", sym->GetTitle(), mcsym->GetTitle()));
 688       t1->Divide(mc);
 689       t2->Divide(mcsym);
 690       t1->SetMarkerColor(mc->GetMarkerColor());
 691       t2->SetMarkerColor(mcsym->GetMarkerColor());
 692       ratios->Add(t1);
 693       ratios->Add(t2);
 694     }
 695
 696     // Check if we have ratios
 697     Bool_t   hasRatios = (ratios->GetHists() &&
 698                           (ratios->GetHists()->GetEntries() > 0));
 699 #if 0
 700     Info("MakeRatios", "Got a total of %d ratios", !hasRatios ? 0 :
 701          ratios->GetHists()->GetEntries());
 702 #endif
 703
 704     if (!hasRatios) { delete ratios; ratios = 0; }
 705     return ratios;
 706   }
 707
 708   //__________________________________________________________________
 709   /**
 710    * Fix the apperance of the axis in a stack
 711    *
 712    * @param stack  stack of histogram
 713    * @param s      Scaling factor
 714    * @param ytitle Y axis title
 715    * @param force  Whether to draw the stack first or not
 716    * @param ynDiv  Divisions on Y axis
 717    */
 718   void FixAxis(THStack* stack, Double_t s, const char* ytitle,
 719                Int_t ynDiv=10, Bool_t force=true)
 720   {
 721     if (force) stack->Draw("nostack e1");
 722
 723     TH1* h = stack->GetHistogram();
 724     if (!h) return;
 725
 726     h->SetXTitle("#eta");
 727     h->SetYTitle(ytitle);
 728     TAxis* xa = h->GetXaxis();
 729     TAxis* ya = h->GetYaxis();
 730     if (xa) {
 731       xa->SetTitle("#eta");
 732       // xa->SetTicks("+-");
 733       xa->SetTitleSize(s*xa->GetTitleSize());
 734       xa->SetLabelSize(s*xa->GetLabelSize());
 735       xa->SetTickLength(s*xa->GetTickLength());
 736     }
 737     if (ya) {
 738       ya->SetTitle(ytitle);
 739       ya->SetDecimals();
 740       // ya->SetTicks("+-");
 741       ya->SetNdivisions(ynDiv);
 742       ya->SetTitleSize(s*ya->GetTitleSize());
 743       ya->SetLabelSize(s*ya->GetLabelSize());
 744     }
 745   }
 746   //__________________________________________________________________
 747   /**
 748    * Compute the ratio of @a h to @a g.  @a g is evaluated at the bin
 749    * centers of @a h
 750    *
 751    * @param h  Numerator
 752    * @param g  Divisor
 753    *
 754    * @return h/g
 755    */
 756   TH1* Ratio(const TH1* h, const TGraph* g) const
 757   {
 758     if (!h || !g) return 0;
 759
 760     TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
 761     ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
 762     ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
 763     ret->Reset();
 764     ret->SetMarkerStyle(h->GetMarkerStyle());
 765     ret->SetMarkerColor(g->GetMarkerColor());
 766     ret->SetMarkerSize(0.9*g->GetMarkerSize());
 767     Double_t xlow  = g->GetX()[0];
 768     Double_t xhigh = g->GetX()[g->GetN()-1];
 769     if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
 770
 771     for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
 772       Double_t c = h->GetBinContent(i);
 773       if (c <= 0) continue;
 774
 775       Double_t x = h->GetBinCenter(i);
 776       if (x < xlow || x > xhigh) continue;
 777
 778       Double_t f = g->Eval(x);
 779       if (f <= 0) continue;
 780
 781       ret->SetBinContent(i, c / f);
 782       ret->SetBinError(i, h->GetBinError(i) / f);
 783     }
 784     if (ret->GetEntries() <= 0) { delete ret; ret = 0; }
 785     return ret;
 786   }
 787
 788   //__________________________________________________________________
 789   /**
 790    * Make an extension of @a h to make it symmetric about 0
 791    *
 792    * @param h Histogram to symmertrice
 793    *
 794    * @return Symmetric extension of @a h
 795    */
 796   TH1* Symmetrice(const TH1* h) const
 797   {
 798     fOut->cd();
 799
 800     Int_t nBins = h->GetNbinsX();
 801     TH1*  s     = new TH1D(Form("%s_mirror", h->GetName()),
 802                            Form("%s (mirrored)", h->GetTitle()),
 803                            nBins,
 804                            -h->GetXaxis()->GetXmax(),
 805                            -h->GetXaxis()->GetXmin());
 806     s->SetMarkerColor(h->GetMarkerColor());
 807     s->SetMarkerSize(h->GetMarkerSize());
 808     s->SetMarkerStyle(h->GetMarkerStyle()+4);
 809     s->SetFillColor(h->GetFillColor());
 810     s->SetFillStyle(h->GetFillStyle());
 811     // s->SetDirectory(0);
 812
 813     // Find the first and last bin with data
 814     Int_t first = nBins+1;
 815     Int_t last  = 0;
 816     for (Int_t i = 1; i <= nBins; i++) {
 817       if (h->GetBinContent(i) <= 0) continue;
 818       first = TMath::Min(first, i);
 819       last  = TMath::Max(last,  i);
 820     }
 821
 822     Double_t xfirst = h->GetBinCenter(first-1);
 823     Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
 824     Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
 825     for (Int_t i = f1, j=l2; i <= last; i++,j--) {
 826       s->SetBinContent(j, h->GetBinContent(i));
 827       s->SetBinError(j, h->GetBinError(i));
 828     }
 829     // Fill in overlap bin
 830     s->SetBinContent(l2+1, h->GetBinContent(first));
 831     s->SetBinError(l2+1, h->GetBinError(first));
 832     return s;
 833   }
 834   //__________________________________________________________________
 835   /**
 836    * Rebin a histogram
 837    *
 838    * @param h     Histogram to rebin
 839    * @param rebin Rebinning factor
 840    *
 841    * @return
 842    */
 843   virtual void Rebin(TH1* h, Int_t rebin) const
 844   {
 845     if (rebin <= 1) return;
 846
 847     Int_t nBins = h->GetNbinsX();
 848     if(nBins % rebin != 0) {
 849       Warning("Rebin", "Rebin factor %d is not a devisor of current number "
 850               "of bins %d in the histogram %s", rebin, nBins, h->GetName());
 851       return;
 852     }
 853
 854     // Make a copy
 855     TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
 856     tmp->Rebin(rebin);
 857     tmp->SetDirectory(0);
 858
 859     // The new number of bins
 860     Int_t nBinsNew = nBins / rebin;
 861     for(Int_t i = 1;i<= nBinsNew; i++) {
 862       Double_t content = 0;
 863       Double_t sumw    = 0;
 864       Double_t wsum    = 0;
 865       Int_t    nbins   = 0;
 866       for(Int_t j = 1; j<=rebin;j++) {
 867         Int_t bin = (i-1)*rebin + j;
 868         if(h->GetBinContent(bin) <= 0) continue;
 869         Double_t c =  h->GetBinContent(bin);
 870         Double_t w = 1 / TMath::Power(c,2);
 871         content    += c;
 872         sumw       += w;
 873         wsum       += w * c;
 874         nbins++;
 875       }
 876
 877       if(content > 0 ) {
 878         tmp->SetBinContent(i, wsum / sumw);
 879         tmp->SetBinError(i,TMath::Sqrt(sumw));
 880       }
 881     }
 882
 883     // Finally, rebin the histogram, and set new content
 884     h->Rebin(rebin);
 885     for(Int_t i =1;i<=nBinsNew; i++) {
 886       h->SetBinContent(i,tmp->GetBinContent(i));
 887       // h->SetBinError(i,tmp->GetBinError(i));
 888     }
 889
 890     delete tmp;
 891   }
 892 };
 893
 894 //____________________________________________________________________
 895 //
 896 // EOF
 897 //
 898