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