Minor fixes
[u/mrichter/AliRoot.git] / 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