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