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