PWG2/FORWARD/analysis2/DrawRes.C

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