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