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