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