PWG2/FORWARD/analysis2/DrawRes.C

   1 #include <TH1D.h>
   2 #include <TH2D.h>
   3 #include <TCanvas.h>
   4 #include <TPad.h>
   5 #include <TFile.h>
   6 #include <TTree.h>
   7 #include <TError.h>
   8 #include <TStyle.h>
   9 #include <THStack.h>
  10 #include <TLegend.h>
  11 #include <TMath.h>
  12 #include <TMultiGraph.h>
  13 #include <TGraph.h>
  14 #include <TGraphErrors.h>
  15 #include <TLatex.h>
  16 #include "AliAODForwardMult.h"
  17 #include "OtherData.C"
  18
  19 /**
  20  * Draw the data stored in the AOD
  21  *
  22  * To use this, do
  23  *
  24  * @code
  25  *   Root> .L $ALICE_ROOT/PWG2/FORWARD/analysis2/Compile.C
  26  *   Root> Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/DrawRes.C")
  27  *   Root> DrawRes dr
  28  *   Root> dr.Run("AliAODs.root",-10,10,5,AliAODForwardMult::kInel,900)
  29  * @endcode
  30  *
  31  * The output is stored in a ROOT file
  32  *
  33  * See also the script Pass2.C
  34  *
  35  * @ingroup pwg2_forward_analysis_scripts
  36  */
  37 struct DrawRes
  38 {
  39 public:
  40   /** AOD tree */
  41   TTree*             fTree;
  42   /** AOD object */
  43   AliAODForwardMult* fAOD;
  44   /** Output file */
  45   TFile*             fOut;
  46   /** Summed histogram */
  47   TH2D*              fSum;
  48   /** Vertex efficiency */
  49   Double_t           fVtxEff;
  50   /** Title to put on the plot */
  51   TString            fTitle;
  52
  53   //__________________________________________________________________
  54   /**
  55    * Constructor
  56    *
  57    */
  58   DrawRes()
  59     : fTree(0),
  60       fAOD(0),
  61       fOut(0),
  62       fSum(0),
  63       fVtxEff(0),
  64       fTitle("")
  65   {}
  66   //__________________________________________________________________
  67   /**
  68    * Reset internal structures
  69    *
  70    */
  71   void Clear(Option_t* )
  72   {
  73     if (fTree && fTree->GetCurrentFile()) {
  74       fTree->GetCurrentFile()->Close();
  75       delete fTree;
  76     }
  77     if (fOut) {
  78       fOut->Close();
  79       delete fOut;
  80     }
  81     if (fSum) delete fSum;
  82     fTree   = 0;
  83     fOut    = 0;
  84     fSum    = 0;
  85     fVtxEff = 0;
  86
  87   }
  88   //__________________________________________________________________
  89   /**
  90    * Run
  91    *
  92    * @param file   Input file with AOD tree
  93    * @param vzMin  Minimum interaction point z coordinate
  94    * @param vzMax  Maximum interaction point z coordinate
  95    * @param rebin  How many times to re-bin the @f$ dN_{ch}/d\eta@f$
  96    * @param mask   Trigger mask
  97    * @param energy Collision energy @f$ \sqrt{s}@f$ or @f$\sqrt{s_{NN}}@f$
  98    * @param title  Title to put on the plot
  99    *
 100    * @return True on success, false otherwise
 101    */
 102   Bool_t Run(const char* file="AliAODs.root",
 103              Double_t vzMin=-10, Double_t vzMax=10, Int_t rebin=1,
 104              Int_t mask=AliAODForwardMult::kInel, Int_t energy=900,
 105              const char* title="")
 106   {
 107     TString trgName;
 108     if    (mask & AliAODForwardMult::kInel)    trgName.Append("inel-");
 109     if    (mask & AliAODForwardMult::kInelGt0) trgName.Append("inelgt0-");
 110     if    (mask & AliAODForwardMult::kNSD)     trgName.Append("nsd-");
 111     if (trgName.EndsWith("-")) trgName.Remove(trgName.Length()-1);
 112     if (trgName.IsNull()) trgName = "unknown";
 113     TString outName =
 114       TString::Format("dndeta_%04dGeV_%c%02d-%c%02dcm_rb%02d_%s.root",
 115                       energy,
 116                       vzMin < 0 ? 'm' : 'p', Int_t(TMath::Abs(vzMin)+.5),
 117                       vzMax < 0 ? 'm' : 'p', Int_t(TMath::Abs(vzMax)+.5),
 118                       rebin, trgName.Data());
 119     fTitle = title;
 120     if (!Open(file, outName)) return kFALSE;
 121     if (!Process(vzMin,vzMax,mask)) return kFALSE;
 122     if (!Finish(rebin, mask, energy)) return kFALSE;
 123
 124     return kTRUE;
 125   }
 126   //__________________________________________________________________
 127   /**
 128    * Open the files @a fname containg the tree with AliAODEvent objects,
 129    * and also open the output file @a outname for writting.
 130    *
 131    * @param fname   Input file name
 132    * @param outname Output file name
 133    *
 134    * @return true on success, false otherwise
 135    */
 136   Bool_t Open(const char* fname, const char* outname)
 137   {
 138     Clear("");
 139
 140     TFile* file = TFile::Open(fname, "READ");
 141     if (!file) {
 142       Error("Open", "Cannot open file %s", fname);
 143       return kFALSE;
 144     }
 145
 146     // Get the AOD tree
 147     fTree = static_cast<TTree*>(file->Get("aodTree"));
 148     if (!fTree) {
 149       Error("Init", "Couldn't get the tree");
 150       return kFALSE;
 151     }
 152
 153     // Set the branch pointer
 154     fTree->SetBranchAddress("Forward", &fAOD);
 155
 156     fOut = TFile::Open(outname, "RECREATE");
 157     if (!fOut) {
 158       Error("Open", "Couldn't open %s", outname);
 159       return kFALSE;
 160     }
 161     return kTRUE;
 162   }
 163   //__________________________________________________________________
 164   /**
 165    * Process the events
 166    *
 167    * @param vzMin  Minimum interaction point z coordinate
 168    * @param vzMax  Maximum interaction point z coordinate
 169    * @param mask   Trigger mask
 170    *
 171    * @return true on success, false otherwise
 172    */
 173   Bool_t Process(Double_t vzMin, Double_t vzMax, Int_t mask)
 174   {
 175     Int_t nTriggered  = 0;                    // # of triggered ev.
 176     Int_t nWithVertex = 0;                    // # of ev. w/vertex
 177     Int_t nAccepted   = 0;                    // # of ev. used
 178     Int_t nAvailable  = fTree->GetEntries();  // How many entries
 179
 180     for (int i = 0; i < nAvailable; i++) {
 181       fTree->GetEntry(i);
 182
 183       // Create sum histogram on first event - to match binning to input
 184       if (!fSum) {
 185         fSum = static_cast<TH2D*>(fAOD->GetHistogram().Clone("d2ndetadphi"));
 186         Info("Process", "Created sum histogram (%d,%f,%f)x(%d,%f,%f)",
 187              fSum->GetNbinsX(),
 188              fSum->GetXaxis()->GetXmin(),
 189              fSum->GetXaxis()->GetXmax(),
 190              fSum->GetNbinsY(),
 191              fSum->GetYaxis()->GetXmin(),
 192              fSum->GetYaxis()->GetXmax());
 193       }
 194
 195       // fAOD->Print();
 196       // Other trigger/event requirements could be defined
 197       if (!fAOD->IsTriggerBits(mask)) continue;
 198       nTriggered++;
 199
 200       // Check if there's a vertex
 201       if (!fAOD->HasIpZ()) continue;
 202       nWithVertex++;
 203
 204       // Select vertex range (in centimeters)
 205       if (!fAOD->InRange(vzMin, vzMax)) continue;
 206       nAccepted++;
 207
 208       // Add contribution from this event
 209       fSum->Add(&(fAOD->GetHistogram()));
 210     }
 211     fVtxEff = Double_t(nWithVertex)/nTriggered;
 212
 213     Info("Process", "Total of %9d events\n"
 214          "               %9d has a trigger\n"
 215          "               %9d has a vertex\n"
 216          "               %9d was used\n",
 217          nAvailable, nTriggered, nWithVertex, nAccepted);
 218
 219     return kTRUE;
 220   }
 221   //__________________________________________________________________
 222   /**
 223    * Finish the stuff and draw
 224    *
 225    * @param rebin  How many times to rebin
 226    * @param energy Collision energy
 227    * @param mask   Trigger mask
 228    *
 229    * @return true on success, false otherwise
 230    */
 231   Bool_t Finish(Int_t rebin, Int_t mask, Int_t energy)
 232   {
 233     fOut->cd();
 234
 235     // Get acceptance normalisation from underflow bins
 236     TH1D* norm   = fSum->ProjectionX("norm", 0, 1, "");
 237     // Project onto eta axis - _ignoring_underflow_bins_!
 238     TH1D* dndeta = fSum->ProjectionX("dndeta", 1, -1, "e");
 239     dndeta->SetTitle("1/N dN_{ch}/d#eta|_{forward}");
 240     // Normalize to the acceptance
 241     dndeta->Divide(norm);
 242     // Scale by the vertex efficiency
 243     dndeta->Scale(fVtxEff, "width");
 244     dndeta->SetMarkerColor(kRed+1);
 245     dndeta->SetMarkerStyle(20);
 246     dndeta->SetMarkerSize(1);
 247     dndeta->SetFillStyle(0);
 248
 249     TH1* sym = Symmetrice(dndeta);
 250
 251     Rebin(dndeta, rebin);
 252     Rebin(sym, rebin);
 253
 254     THStack* stack = new THStack("results", "Results");
 255     stack->Add(dndeta);
 256     stack->Add(sym);
 257
 258     TString hhdName(fOut->GetName());
 259     hhdName.ReplaceAll("dndeta", "hhd");
 260     TH1* hhd    = GetHHD(hhdName.Data(), mask & AliAODForwardMult::kNSD);
 261     TH1* hhdsym = 0;
 262     if (hhd) {
 263       hhdsym = Symmetrice(hhd);
 264       stack->Add(hhd);
 265       stack->Add(hhdsym);
 266     }
 267
 268     TMultiGraph* mg     = GetData(energy, mask);
 269
 270     gStyle->SetOptTitle(0);
 271     Double_t yd = (mg->GetListOfGraphs()->GetEntries() ? 0.3 : 0);
 272     TCanvas* c = new TCanvas("c", "C", 900, 800);
 273     c->SetFillColor(0);
 274     c->SetBorderMode(0);
 275     c->SetBorderSize(0);
 276     c->cd();
 277
 278     TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0);
 279     p1->SetTopMargin(0.05);
 280     p1->SetBottomMargin(0.001);
 281     p1->SetRightMargin(0.05);
 282     p1->SetGridx();
 283     p1->SetTicks(1,1);
 284     p1->Draw();
 285     p1->cd();
 286     stack->SetMinimum(-0.1);
 287     FixAxis(stack, 1/(1-yd)/1.5, "#frac{1}{N} #frac{dN_{ch}}{#eta}");
 288     p1->Clear();
 289     stack->DrawClone("nostack e1");
 290
 291
 292     THStack* ratios = new THStack("ratios", "Ratios");
 293     TGraphAsymmErrors* o = 0;
 294     TIter nextG(mg->GetListOfGraphs());
 295     while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
 296       o->Draw("same p");
 297       ratios->Add(Ratio(dndeta, o));
 298       ratios->Add(Ratio(sym,    o));
 299       ratios->Add(Ratio(hhd,    o));
 300       ratios->Add(Ratio(hhdsym, o));
 301     }
 302     if (hhd && hhdsym) {
 303       TH1F* t1 = static_cast<TH1F*>(dndeta->Clone(Form("%s_%s",
 304                                                        dndeta->GetName(),
 305                                                        hhd->GetName())));
 306       t1->SetTitle(Form("%s / %s", dndeta->GetTitle(), hhd->GetTitle()));
 307       TH1F* t2 = static_cast<TH1F*>(sym->Clone(Form("%s_%s",
 308                                                     sym->GetName(),
 309                                                     hhdsym->GetName())));
 310       t2->SetTitle(Form("%s / %s", sym->GetTitle(), hhdsym->GetTitle()));
 311       t1->Divide(hhd);
 312       t2->Divide(hhdsym);
 313       ratios->Add(t1);
 314       ratios->Add(t2);
 315     }
 316
 317     // Make a legend
 318     TString trigs(AliAODForwardMult::GetTriggerString(mask));
 319     TLegend* l = p1->BuildLegend(.3,p1->GetBottomMargin()+.01,.7,.4,
 320                                  Form("#sqrt{s}=%dGeV, %s",
 321                                       energy, trigs.Data()));
 322     l->SetFillColor(0);
 323     l->SetFillStyle(0);
 324     l->SetBorderSize(0);
 325     p1->cd();
 326     c->cd();
 327
 328     if (!ratios->GetHists() || ratios->GetHists()->GetEntries() <= 0) {
 329       p1->SetPad(0, 0, 1, 1);
 330       p1->SetBottomMargin(0.1);
 331       l->SetY1(0.11);
 332       FixAxis(stack, (1-yd)/1,  "#frac{1}{N} #frac{dN_{ch}}{#eta}",false);
 333       p1->cd();
 334       c->cd();
 335     }
 336     else {
 337       TPad* p2 = new TPad("p2", "p2", 0, 0.0, 1.0, yd, 0, 0);
 338       p2->SetTopMargin(0.001);
 339       p2->SetRightMargin(0.05);
 340       p2->SetBottomMargin(1/yd * 0.07);
 341       p2->SetGridx();
 342       p2->SetTicks(1,1);
 343       p2->Draw();
 344       p2->cd();
 345       FixAxis(ratios, 1/yd/1.5, "Ratios");
 346       p2->Clear();
 347       ratios->DrawClone("nostack e1");
 348
 349       TLegend* l2 = p2->BuildLegend(.3,p2->GetBottomMargin()+.01,.7,.6);
 350       l2->SetFillColor(0);
 351       l2->SetFillStyle(0);
 352       l2->SetBorderSize(0);
 353
 354       p2->cd();
 355       TGraphErrors* band = new TGraphErrors(2);
 356       band->SetPoint(0, sym->GetXaxis()->GetXmin(), 1);
 357       band->SetPoint(1, dndeta->GetXaxis()->GetXmax(), 1);
 358       band->SetPointError(0, 0, .1);
 359       band->SetPointError(1, 0, .1);
 360       band->SetFillColor(kYellow+2);
 361       band->SetFillStyle(3002);
 362       band->Draw("3 same");
 363       ratios->DrawClone("nostack e1 same");
 364
 365       c->cd();
 366     }
 367     p1->cd();
 368     TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
 369     tit->SetNDC();
 370     tit->SetTextSize(0.05);
 371     tit->Draw();
 372
 373     TString imgName(fOut->GetName());
 374     imgName.ReplaceAll(".root", ".png");
 375     c->SaveAs(imgName.Data());
 376
 377     stack->Write();
 378     mg->Write();
 379     ratios->Write();
 380
 381     fOut->Close();
 382
 383     return kTRUE;
 384   }
 385   //__________________________________________________________________
 386   /**
 387    * Get the result from previous analysis code
 388    *
 389    * @param fn File to open
 390    *
 391    * @return null or result of previous analysis code
 392    */
 393   TH1* GetHHD(const char* fn="fmd_dNdeta_mult.root", bool nsd=false)
 394   {
 395     TDirectory* savdir = gDirectory;
 396     TFile* file = TFile::Open(fn, "READ");
 397     if (!file) {
 398       Warning("GetHHD", "couldn't open HHD file %s", fn);
 399       savdir->cd();
 400       return 0;
 401     }
 402     TString hist(Form("dNdeta_dNdeta%s", (nsd ? "NSD" : "")));
 403     TH1* h = static_cast<TH1*>(file->Get(hist.Data()));
 404     if (!h) {
 405       Warning("GetHHD", "Couldn't find HHD histogram %s in %s",
 406               hist.Data(), fn);
 407       file->ls();
 408       file->Close();
 409       savdir->cd();
 410       return 0;
 411     }
 412     TH1* r = static_cast<TH1*>(h->Clone("dndeta_hhd"));
 413     r->SetTitle("1/N dN_{ch}/d#eta (HHD)");
 414     r->SetFillStyle(0);
 415     r->SetFillColor(0);
 416     r->SetMarkerStyle(21);
 417     r->SetMarkerColor(kPink+1);
 418     r->SetDirectory(savdir);
 419
 420     file->Close();
 421     savdir->cd();
 422     return r;
 423   }
 424   //__________________________________________________________________
 425   /**
 426    * Fix the apperance of the axis in a stack
 427    *
 428    * @param stack  stack of histogram
 429    * @param s      Scaling factor
 430    * @param ytitle Y axis title
 431    * @param force  Whether to draw the stack first or not
 432    */
 433   void FixAxis(THStack* stack, Double_t s, const char* ytitle,
 434                Bool_t force=true)
 435   {
 436     if (force) stack->Draw("nostack e1");
 437
 438     TH1* h = stack->GetHistogram();
 439     if (!h) return;
 440
 441     h->SetXTitle("#eta");
 442     h->SetYTitle(ytitle);
 443     TAxis* xa = h->GetXaxis();
 444     TAxis* ya = h->GetYaxis();
 445     if (xa) {
 446       xa->SetTitle("#eta");
 447       // xa->SetTicks("+-");
 448       xa->SetTitleSize(s*xa->GetTitleSize());
 449       xa->SetLabelSize(s*xa->GetLabelSize());
 450       xa->SetTickLength(s*xa->GetTickLength());
 451     }
 452     if (ya) {
 453       ya->SetTitle(ytitle);
 454       // ya->SetTicks("+-");
 455       ya->SetNdivisions(10);
 456       ya->SetTitleSize(s*ya->GetTitleSize());
 457       ya->SetLabelSize(s*ya->GetLabelSize());
 458     }
 459   }
 460   //__________________________________________________________________
 461   /**
 462    * Compute the ratio of @a h to @a g.  @a g is evaluated at the bin
 463    * centers of @a h
 464    *
 465    * @param h  Numerator
 466    * @param g  Divisor
 467    *
 468    * @return h/g
 469    */
 470   TH1* Ratio(const TH1* h, const TGraph* g) const
 471   {
 472     if (!h || !g) return 0;
 473
 474     TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
 475     ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
 476     ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
 477     ret->Reset();
 478     ret->SetMarkerStyle(h->GetMarkerStyle());
 479     ret->SetMarkerColor(g->GetMarkerColor());
 480     ret->SetMarkerSize(0.7*g->GetMarkerSize());
 481     Double_t xlow  = g->GetX()[0];
 482     Double_t xhigh = g->GetX()[g->GetN()-1];
 483     if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
 484
 485     for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
 486       Double_t c = h->GetBinContent(i);
 487       if (c <= 0) continue;
 488
 489       Double_t x = h->GetBinCenter(i);
 490       if (x < xlow || x > xhigh) continue;
 491
 492       Double_t f = g->Eval(x);
 493       if (f <= 0) continue;
 494
 495       ret->SetBinContent(i, c / f);
 496       ret->SetBinError(i, h->GetBinError(i) / f);
 497     }
 498     if (ret->GetEntries() <= 0) { delete ret; ret = 0; }
 499     return ret;
 500   }
 501
 502   //__________________________________________________________________
 503   /**
 504    * Make an extension of @a h to make it symmetric about 0
 505    *
 506    * @param h Histogram to symmertrice
 507    *
 508    * @return Symmetric extension of @a h
 509    */
 510   TH1* Symmetrice(const TH1* h) const
 511   {
 512     Int_t nBins = h->GetNbinsX();
 513     TH1*  s     = new TH1D(Form("%s_mirror", h->GetName()),
 514                            Form("%s (mirrored)", h->GetTitle()),
 515                            nBins,
 516                            -h->GetXaxis()->GetXmax(),
 517                            -h->GetXaxis()->GetXmin());
 518     s->SetMarkerColor(h->GetMarkerColor());
 519     s->SetMarkerSize(h->GetMarkerSize());
 520     s->SetMarkerStyle(h->GetMarkerStyle()+4);
 521     s->SetFillColor(h->GetFillColor());
 522     s->SetFillStyle(h->GetFillStyle());
 523
 524     // Find the first and last bin with data
 525     Int_t first = nBins+1;
 526     Int_t last  = 0;
 527     for (Int_t i = 1; i <= nBins; i++) {
 528       if (h->GetBinContent(i) <= 0) continue;
 529       first = TMath::Min(first, i);
 530       last  = TMath::Max(last,  i);
 531     }
 532
 533     Double_t xfirst = h->GetBinCenter(first-1);
 534     Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
 535     Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
 536     for (Int_t i = f1, j=l2; i <= last; i++,j--) {
 537       s->SetBinContent(j, h->GetBinContent(i));
 538       s->SetBinError(j, h->GetBinError(i));
 539     }
 540     return s;
 541   }
 542   //__________________________________________________________________
 543   /**
 544    * Rebin a histogram
 545    *
 546    * @param h     Histogram to rebin
 547    * @param rebin Rebinning factor
 548    *
 549    * @return
 550    */
 551   virtual void Rebin(TH1* h, Int_t rebin) const
 552   {
 553     if (rebin <= 1) return;
 554
 555     Int_t nBins = h->GetNbinsX();
 556     if(nBins % rebin != 0) {
 557       Warning("Rebin", "Rebin factor %d is not a devisor of current number "
 558               "of bins %d in the histogram %s", rebin, nBins, h->GetName());
 559       return;
 560     }
 561
 562     // Make a copy
 563     TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
 564     tmp->Rebin(rebin);
 565     tmp->SetDirectory(0);
 566
 567     // The new number of bins
 568     Int_t nBinsNew = nBins / rebin;
 569     for(Int_t i = 1;i<= nBinsNew; i++) {
 570       Double_t content = 0;
 571       Double_t sumw    = 0;
 572       Double_t wsum    = 0;
 573       Int_t    nbins   = 0;
 574       for(Int_t j = 1; j<=rebin;j++) {
 575         Int_t bin = (i-1)*rebin + j;
 576         if(h->GetBinContent(bin) <= 0) continue;
 577         Double_t c =  h->GetBinContent(bin);
 578         Double_t w = 1 / TMath::Power(c,2);
 579         content    += c;
 580         sumw       += w;
 581         wsum       += w * c;
 582         nbins++;
 583       }
 584
 585       if(content > 0 ) {
 586         tmp->SetBinContent(i, wsum / sumw);
 587         tmp->SetBinError(i,TMath::Sqrt(sumw));
 588       }
 589     }
 590
 591     // Finally, rebin the histogram, and set new content
 592     h->Rebin(rebin);
 593     for(Int_t i =1;i<=nBinsNew; i++) {
 594       h->SetBinContent(i,tmp->GetBinContent(i));
 595       // h->SetBinError(i,tmp->GetBinError(i));
 596     }
 597
 598     delete tmp;
 599   }
 600 };
 601
 602 //____________________________________________________________________
 603 //
 604 // EOF
 605 //
 606