]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGHF/hfe/AliHFEInclusiveSpectrumQA.cxx
Bug fix
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEInclusiveSpectrumQA.cxx
1
2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 *                                                                        *
5 * Author: The ALICE Off-line Project.                                    *
6 * Contributors are mentioned in the code where appropriate.              *
7 *                                                                        *
8 * Permission to use, copy, modify and distribute this software and its   *
9 * documentation strictly for non-commercial purposes is hereby granted   *
10 * without fee, provided that the above copyright notice appears in all   *
11 * copies and that both the copyright notice and this permission notice   *
12 * appear in the supporting documentation. The authors make no claims     *
13 * about the suitability of this software for any purpose. It is          *
14 * provided "as is" without express or implied warranty.                  *
15 **************************************************************************/
16 //
17 // Class for spectrum correction
18 // Subtraction of hadronic background, Unfolding of the data and
19 // Renormalization done here
20 // The following containers have to be set:
21 //  - Correction framework container for real data
22 //  - Correction framework container for MC (Efficiency Map)
23 //  - Correction framework container for background coming from data
24 //  - Correction framework container for background coming from MC
25 //
26 //  Author: 
27 //            Raphaelle Bailhache <R.Bailhache@gsi.de>
28 //
29
30 #include <TArrayD.h>
31 #include <TH1.h>
32 #include <TList.h>
33 #include <TObjArray.h>
34 #include <TObject.h>
35 #include <TROOT.h>
36 #include <TCanvas.h>
37 #include <TLegend.h>
38 #include <TStyle.h>
39 #include <TMath.h>
40 #include <TAxis.h>
41 #include <TGraphErrors.h>
42 #include <TFile.h>
43 #include <TPad.h>
44 #include <TH2D.h>
45 #include <TF1.h>
46
47 #include "AliPID.h"
48 #include "AliCFContainer.h"
49 #include "AliCFDataGrid.h"
50 #include "AliCFEffGrid.h"
51 #include "AliCFGridSparse.h"
52 #include "AliCFUnfolding.h"
53 #include "AliLog.h"
54
55 #include "AliHFEInclusiveSpectrumQA.h"
56 #include "AliHFECorrectSpectrumBase.h"
57 #include "AliHFEcuts.h"
58 #include "AliHFEcontainer.h"
59 #include "AliHFEtools.h"
60
61 ClassImp(AliHFEInclusiveSpectrumQA)
62
63 const Char_t *AliHFEInclusiveSpectrumQA::fgkNameCanvas[AliHFEInclusiveSpectrumQA::kNTypeEfficiency] = {
64   "V0Efficiency",
65   "MCEfficiency",
66   "ParametrizedEfficiency"
67 };
68
69 //____________________________________________________________
70 AliHFEInclusiveSpectrumQA::AliHFEInclusiveSpectrumQA():
71   TNamed(),
72   fPtMax(7.0),
73   fListOfResult(),
74   fWriteToFile(kTRUE)
75 {
76   //
77   // Default constructor
78   //
79
80   fListOfResult = new TObjArray(kNResults);
81   fListOfResult->SetName("ListOfResults");
82  
83
84 }
85 //____________________________________________________________
86 AliHFEInclusiveSpectrumQA::AliHFEInclusiveSpectrumQA(const char *name):
87   TNamed(name, ""),
88   fPtMax(7.0),
89   fListOfResult(),
90   fWriteToFile(kTRUE)
91 {
92   //
93   // Default constructor
94   //
95
96   fListOfResult = new TObjArray(kNResults);
97   fListOfResult->SetName("ListOfResults");
98  
99
100 }
101
102 //____________________________________________________________
103 AliHFEInclusiveSpectrumQA::~AliHFEInclusiveSpectrumQA(){
104   //
105   // Destructor
106   //
107   if(fListOfResult) delete fListOfResult;
108  
109 }
110 //____________________________________________________________
111 void AliHFEInclusiveSpectrumQA::AddResultAt(TObject *obj,Int_t index)
112 {
113   //
114   // Init what we need for the correction:
115   //
116
117   if(fListOfResult) fListOfResult->AddAt(obj,index);
118
119 }
120 //____________________________________________________________
121 TObject *AliHFEInclusiveSpectrumQA::GetResult(Int_t index)
122 {
123   //
124   // Get result
125   //
126
127   if(fListOfResult) return fListOfResult->UncheckedAt(index);
128   else return 0x0;
129
130 }
131 //____________________________________________________________
132 void AliHFEInclusiveSpectrumQA::DrawProjections() const
133 {
134   //
135   // get spectrum for beauty 2nd method
136   //
137   //
138   AliCFContainer *data = (AliCFContainer *) fListOfResult->UncheckedAt(kDataProjection);
139   THnSparseF *correlation = (THnSparseF *) fListOfResult->UncheckedAt(kCMProjection);
140   if(!data || !correlation) return;
141
142   Int_t ndimcont = data->GetNVar();
143   Int_t ndimcor = correlation->GetNdimensions();
144   Int_t charge = 3;
145   Int_t centrality = 5;
146   Int_t eta = 1;
147
148   TCanvas * canvas = new TCanvas("Projections","Projections",1000,700);
149   Int_t n = 0;
150   if(charge < ndimcont) n++;
151   if(centrality < ndimcont) n++;
152   if(eta < ndimcont) n++;
153   canvas->Divide(2,n);
154   Int_t counter = 1;
155
156   if(charge < ndimcont) {
157    
158     canvas->cd(counter);
159     TH1 *checkcharge = (TH1 *) data->Project(data->GetNStep()-1,charge);
160     checkcharge->Draw();
161     counter++;
162     canvas->cd(counter);
163     TH2F* projectioncharge = (TH2F *) correlation->Projection(charge,charge+((Int_t)(ndimcor/2.)));
164     projectioncharge->Draw("colz");
165     counter++;  
166
167   }
168
169   if(centrality < ndimcont) {
170     canvas->cd(counter);
171     TH1 *checkcentrality = (TH1 *) data->Project(data->GetNStep()-1,centrality);
172     checkcentrality->Draw();
173     counter++;
174     canvas->cd(counter);
175     TH2F *projectioncentrality = (TH2F *) correlation->Projection(centrality,centrality+((Int_t)(ndimcor/2.)));
176     projectioncentrality->Draw("colz");
177     counter++;  
178   }
179
180   if(eta < ndimcont) {
181     canvas->cd(counter);
182     TH1 *checketa = (TH1 *) data->Project(data->GetNStep()-1,eta);
183     checketa->Draw();
184     counter++;
185     canvas->cd(counter);
186     TH2D* projectioneta = (TH2D *) correlation->Projection(eta,eta+((Int_t)(ndimcor/2.)));
187     projectioneta->Draw("colz");
188   }
189
190
191 }
192 //____________________________________________________________
193 void AliHFEInclusiveSpectrumQA::DrawSubtractContamination() const
194 {
195   //
196   // get spectrum for beauty 2nd method
197   //
198   //
199   TH1D *measuredTH1Daftersubstraction = (TH1D *) fListOfResult->UncheckedAt(kAfterSC);
200   TH1D *measuredTH1Dbeforesubstraction = (TH1D *) fListOfResult->UncheckedAt(kBeforeSC);
201   if(!measuredTH1Daftersubstraction || !measuredTH1Dbeforesubstraction) return;
202
203   SetStyle();
204
205   TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
206   cbackgroundsubtraction->Divide(2,1);
207   cbackgroundsubtraction->cd(1);
208   gPad->SetLogy();
209   gPad->SetTicks();
210   measuredTH1Daftersubstraction->SetStats(0);
211   measuredTH1Daftersubstraction->SetTitle("");
212   measuredTH1Daftersubstraction->GetYaxis()->SetTitleOffset(1.5);
213   measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
214   measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
215   measuredTH1Daftersubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
216   measuredTH1Daftersubstraction->SetMarkerStyle(25);
217   measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
218   measuredTH1Daftersubstraction->SetLineColor(kBlack);
219   measuredTH1Dbeforesubstraction->SetStats(0);
220   measuredTH1Dbeforesubstraction->SetTitle("");
221   measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
222   measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
223   measuredTH1Dbeforesubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
224   measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
225   measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
226   measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
227   measuredTH1Daftersubstraction->Draw();
228   measuredTH1Dbeforesubstraction->Draw("same");
229   TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
230   legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
231   legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
232   legsubstraction->SetFillStyle(0);
233   legsubstraction->SetLineStyle(0);
234   legsubstraction->SetLineColor(0);
235   legsubstraction->Draw("same");
236   cbackgroundsubtraction->cd(2);
237   gPad->SetLogy();
238   gPad->SetTicks();
239   TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
240   ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
241   ratiomeasuredcontamination->SetTitle("");
242   ratiomeasuredcontamination->GetYaxis()->SetTitleOffset(1.5);
243   ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
244   ratiomeasuredcontamination->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
245   ratiomeasuredcontamination->GetYaxis()->SetRangeUser(0.8,1.2);
246   ratiomeasuredcontamination->GetXaxis()->SetRangeUser(0.0,fPtMax);
247   ratiomeasuredcontamination->Sumw2();
248   ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
249   ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
250   ratiomeasuredcontamination->SetStats(0);
251   ratiomeasuredcontamination->SetMarkerStyle(26);
252   ratiomeasuredcontamination->SetMarkerColor(kBlack);
253   ratiomeasuredcontamination->SetLineColor(kBlack);
254   for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
255     ratiomeasuredcontamination->SetBinError(k+1,0.0);
256   }
257   ratiomeasuredcontamination->Draw("P");
258   if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.png");
259
260 }
261
262
263 //____________________________________________________________
264 void AliHFEInclusiveSpectrumQA::DrawCorrectWithEfficiency(Int_t typeeff) const
265 {
266   //
267   // Correct the spectrum for efficiency and unfolding
268   // with both method and compare
269   //
270   
271   TH1D *afterE = 0x0;
272   TH1D *beforeE = 0x0;
273   TH1D *efficiencyDproj = 0x0;
274   TF1 *efficiencyparametrized = 0x0;
275
276   if(typeeff== kV0) {
277     afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterV0);
278     beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforeV0);
279     efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kV0Efficiency);
280   }
281   if(typeeff== kMC) {
282     afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
283     beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforeMCE);
284     efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kMCEfficiency);
285   }
286  if(typeeff== kParametrized) {
287     afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterPE);
288     beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforePE);
289     efficiencyparametrized = (TF1 *) fListOfResult->UncheckedAt(kPEfficiency);
290   }
291
292  if((typeeff==kV0 || typeeff==kMC) && (!afterE || !beforeE || !efficiencyDproj)) return;
293  if(typeeff==kParametrized && (!afterE || !beforeE || !efficiencyparametrized)) return;
294
295   SetStyle();
296
297   TCanvas * cEfficiency = new TCanvas(AliHFEInclusiveSpectrumQA::fgkNameCanvas[typeeff],AliHFEInclusiveSpectrumQA::fgkNameCanvas[typeeff],1000,700);
298   cEfficiency->Divide(2,1);
299   cEfficiency->cd(1);
300   gPad->SetLogy();
301   gPad->SetTicks();
302   afterE->SetStats(0);
303   afterE->SetTitle("");
304   afterE->GetYaxis()->SetTitleOffset(1.5);
305   afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
306   afterE->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
307   afterE->GetXaxis()->SetRangeUser(0.0,fPtMax);
308   afterE->SetMarkerStyle(25);
309   afterE->SetMarkerColor(kBlack);
310   afterE->SetLineColor(kBlack);
311   beforeE->SetStats(0);
312   beforeE->SetTitle("");
313   beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
314   beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
315   beforeE->GetXaxis()->SetRangeUser(0.0,fPtMax);
316   beforeE->SetMarkerStyle(24);
317   beforeE->SetMarkerColor(kBlue);
318   beforeE->SetLineColor(kBlue);
319   afterE->Draw();
320   beforeE->Draw("same");
321   TLegend *legefficiency = new TLegend(0.4,0.6,0.89,0.89);
322   legefficiency->AddEntry(beforeE,"Before Efficiency correction","p");
323   legefficiency->AddEntry(afterE,"After Efficiency correction","p");
324   legefficiency->SetFillStyle(0);
325   legefficiency->SetLineStyle(0);
326   legefficiency->SetLineColor(0);
327   legefficiency->Draw("same");
328   cEfficiency->cd(2);
329   gPad->SetTicks();
330   if((typeeff==kV0 || typeeff==kMC)) {
331     efficiencyDproj->SetTitle("");
332     efficiencyDproj->SetStats(0);
333     efficiencyDproj->GetYaxis()->SetTitleOffset(1.5);
334     efficiencyDproj->GetYaxis()->SetRangeUser(0.0,1.0);
335     efficiencyDproj->GetYaxis()->SetTitle("Efficiency");
336     efficiencyDproj->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
337     efficiencyDproj->GetXaxis()->SetRangeUser(0.0,fPtMax);
338     efficiencyDproj->SetMarkerStyle(25);
339     efficiencyDproj->Draw();
340   }
341   if(typeeff==kParametrized) {
342     efficiencyparametrized->Draw();
343   }
344   
345   if(fWriteToFile) {
346     if(typeeff==kV0) cEfficiency->SaveAs("EfficiencyV0.png");
347     if(typeeff==kMC) cEfficiency->SaveAs("EfficiencyMC.png");
348     if(typeeff==kParametrized) cEfficiency->SaveAs("EfficiencyParametrized.png");
349   }
350
351 }
352
353 //____________________________________________________________
354 void AliHFEInclusiveSpectrumQA::DrawUnfolding() const
355 {
356   //
357   // Draw unfolding
358   //
359   TH1D *measuredspectrumD = (TH1D *) fListOfResult->UncheckedAt(kBeforeU);
360   TH1D *residualspectrumD = (TH1D *) fListOfResult->UncheckedAt(kResidualU);
361   TH1D *efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kUEfficiency);
362   THnSparseF *correlation = (THnSparseF *) fListOfResult->UncheckedAt(kCMProjection);
363   
364   if(!measuredspectrumD || !residualspectrumD || !efficiencyDproj || !correlation) return;
365   
366   Int_t ndimcor = (Int_t) correlation->GetNdimensions()/2.;
367
368   SetStyle();
369
370   TCanvas * cunfolding = new TCanvas("unfolding","unfolding",1000,700);
371   cunfolding->Divide(2,2);
372   cunfolding->cd(1);
373   gPad->SetLogy();
374   gPad->SetTicks();
375   residualspectrumD->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
376   residualspectrumD->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
377   residualspectrumD->GetXaxis()->SetRangeUser(0.0,fPtMax);
378   residualspectrumD->SetStats(0);
379   residualspectrumD->SetTitle("");
380   residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
381   residualspectrumD->SetMarkerStyle(26);
382   residualspectrumD->SetMarkerColor(kBlue);
383   residualspectrumD->SetLineColor(kBlue);
384   residualspectrumD->Sumw2();
385   residualspectrumD->Draw("P");
386   measuredspectrumD->SetStats(0);
387   measuredspectrumD->SetTitle("");  
388   measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
389   measuredspectrumD->SetMarkerStyle(25);
390   measuredspectrumD->SetMarkerColor(kBlack);
391   measuredspectrumD->SetLineColor(kBlack);
392   measuredspectrumD->Draw("same");
393   TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
394   legres->AddEntry(residualspectrumD,"Residual","p");
395   legres->AddEntry(measuredspectrumD,"Measured","p");
396   legres->SetFillStyle(0);
397   legres->SetLineStyle(0);
398   legres->SetLineColor(0);
399   legres->Draw("same");
400   cunfolding->cd(2);
401   gPad->SetTicks();
402   TH1D* ratioresidual = (TH1D*)residualspectrumD->Clone();
403   ratioresidual->SetName("ratioresidual");
404   ratioresidual->SetTitle("");
405   ratioresidual->GetYaxis()->SetRangeUser(0.6,1.4);
406   ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
407   ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
408   ratioresidual->Divide(measuredspectrumD);
409   ratioresidual->SetStats(0);
410   ratioresidual->Draw();
411   cunfolding->cd(3);
412   gPad->SetTicks();
413   efficiencyDproj->GetYaxis()->SetTitle("Efficiency");
414   efficiencyDproj->GetXaxis()->SetTitle("p^{MC}_{T} [GeV/c]");
415   efficiencyDproj->GetXaxis()->SetRangeUser(0.0,fPtMax);
416   efficiencyDproj->GetYaxis()->SetRangeUser(0.0,1.0);
417   efficiencyDproj->SetStats(0);
418   efficiencyDproj->SetTitle("");
419   efficiencyDproj->GetYaxis()->SetTitleOffset(1.5);
420   efficiencyDproj->SetMarkerStyle(26);
421   efficiencyDproj->SetMarkerColor(kBlue);
422   efficiencyDproj->SetLineColor(kBlue);
423   efficiencyDproj->Sumw2();
424   efficiencyDproj->Draw("P");
425   cunfolding->cd(4);
426   TH2F *projectioncorr = (TH2F *) correlation->Projection(0,ndimcor);
427   projectioncorr->GetYaxis()->SetTitle("p^{ESD}_{T} [GeV/c]");
428   projectioncorr->GetXaxis()->SetTitle("p^{MC}_{T} [GeV/c]");
429   projectioncorr->GetXaxis()->SetRangeUser(0.0,fPtMax);
430   projectioncorr->GetYaxis()->SetRangeUser(0.0,fPtMax);
431   projectioncorr->SetStats(0);
432   projectioncorr->SetTitle("");
433   projectioncorr->Draw("colz");
434
435   if(fWriteToFile){
436     cunfolding->SaveAs("Unfolding.png");
437   }
438   
439 }
440 //____________________________________________________________
441 void AliHFEInclusiveSpectrumQA::DrawResult()
442 {
443   //
444   // Draw Results
445   //
446   TGraphErrors* correctedspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt(kFinalResultUnfolded);
447   TGraphErrors* alltogetherspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt( kFinalResultDirectEfficiency);
448   if(!correctedspectrumD || !alltogetherspectrumD) return;
449   
450   SetStyle();
451
452   TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
453   ccorrected->Divide(2,1);
454   ccorrected->cd(1);
455   gPad->SetLogy();
456   gPad->SetTicks();
457   correctedspectrumD->SetTitle("");
458   correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
459   correctedspectrumD->GetYaxis()->SetRangeUser(0.000001,100.0);
460   correctedspectrumD->GetXaxis()->SetRangeUser(0.0,fPtMax);
461   correctedspectrumD->SetMarkerStyle(26);
462   correctedspectrumD->SetMarkerColor(kBlue);
463   correctedspectrumD->SetLineColor(kBlue);
464   correctedspectrumD->Draw("AP");
465   alltogetherspectrumD->SetTitle("");
466   alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
467   alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
468   alltogetherspectrumD->SetMarkerStyle(25);
469   alltogetherspectrumD->SetMarkerColor(kBlack);
470   alltogetherspectrumD->SetLineColor(kBlack);
471   alltogetherspectrumD->Draw("P");
472   TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
473   legcorrected->AddEntry(correctedspectrumD,"Unfolded","p");
474   legcorrected->AddEntry(alltogetherspectrumD,"Direct corrected","p");
475   legcorrected->SetFillStyle(0);
476   legcorrected->SetLineStyle(0);
477   legcorrected->SetLineColor(0);
478   legcorrected->Draw("same");
479   ccorrected->cd(2);
480   gPad->SetTicks();
481   TH1D* ratiocorrected = DivideSpectra(correctedspectrumD,alltogetherspectrumD);
482   ratiocorrected->SetName("ratiocorrected");
483   ratiocorrected->SetTitle("");
484   ratiocorrected->GetYaxis()->SetTitleOffset(1.5);
485   ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
486   ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
487   ratiocorrected->GetXaxis()->SetRangeUser(0.0,fPtMax);
488   ratiocorrected->GetYaxis()->SetRangeUser(0.4,1.4);
489   ratiocorrected->SetStats(0);
490   ratiocorrected->Draw();
491   if(fWriteToFile)ccorrected->SaveAs("CorrectedResults.png");
492
493 }
494 //____________________________________________________________
495 TH1D *AliHFEInclusiveSpectrumQA::DivideSpectra(TGraphErrors *ga, TGraphErrors *gb) 
496 {
497   //
498   // Divide Spectra
499   //
500
501   TH1D *afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
502   if(!afterE) return 0x0;
503
504   TH1D *histoB = (TH1D*) afterE->Clone();
505   histoB->Sumw2();
506   histoB->SetName("ratio");
507   TH1D *histoa = (TH1D*) afterE->Clone();
508   histoa->Sumw2();
509   histoa->SetName("a");
510   TH1D *histob = (TH1D*) afterE->Clone();
511   histob->Sumw2();
512   histob->SetName("b");
513   
514   double xa,ya,xb,yb,eya,eyb;
515   Int_t npointsa = ga->GetN();
516   Int_t npointsb = gb->GetN();
517   if(npointsa != npointsb) {
518     printf("Problem the two spectra have not the same number of points\n");
519     return 0x0;
520   }
521   for(Int_t k = 0; k < npointsa; k++){
522     ga->GetPoint(k,xa,ya);
523     gb->GetPoint(k,xb,yb);
524     //
525     Double_t centerhisto = histoa->GetBinCenter(k+1);
526     //
527     if((TMath::Abs(xa-xb) > 0.0001) || (TMath::Abs(xa-centerhisto) > 0.0001)) {
528       printf("Problem not the same x axis\n");
529       return 0x0;
530     }
531     histoa->SetBinContent(k+1,ya);
532     histob->SetBinContent(k+1,yb);
533     //
534     eya = ga->GetErrorY(k);
535     eyb = gb->GetErrorY(k);
536     //
537     histoa->SetBinError(k+1,eya);
538     histob->SetBinError(k+1,eyb);
539    
540   }
541   
542   histoB->Sumw2();
543   histoB->Divide(histoa,histob,1.0,1.0,"B");
544
545   return histoB;  
546
547 }
548 //__________________________________________
549 void AliHFEInclusiveSpectrumQA::SetStyle() const
550 {
551   //
552   // Set style
553   //
554
555   gStyle->SetPalette(1);
556   gStyle->SetOptStat(1111);
557   gStyle->SetPadBorderMode(0);
558   gStyle->SetCanvasColor(10);
559   gStyle->SetPadLeftMargin(0.13);
560   gStyle->SetPadRightMargin(0.13);
561
562 }