]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGHF/hfe/AliHFEBeautySpectrumQA.cxx
Add Classes (Thanks Constantin)
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEBeautySpectrumQA.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 "AliHFEBeautySpectrumQA.h"
56 #include "AliHFECorrectSpectrumBase.h"
57 #include "AliHFEcuts.h"
58 #include "AliHFEcontainer.h"
59 #include "AliHFEtools.h"
60
61 ClassImp(AliHFEBeautySpectrumQA)
62
63 const Char_t *AliHFEBeautySpectrumQA::fgkNameCanvas[AliHFEBeautySpectrumQA::kNTypeEfficiency] = {
64   "MCEfficiency",
65   "ParametrizedEfficiency"
66 };
67
68 //____________________________________________________________
69 AliHFEBeautySpectrumQA::AliHFEBeautySpectrumQA():
70   TNamed(),
71   fPtMax(8.0),
72   fListOfResult(),
73   fWriteToFile(kTRUE)
74   {
75   //
76   // Default constructor
77   //
78
79   fListOfResult = new TObjArray(kNResults);
80   fListOfResult->SetName("ListOfResults");
81  
82
83 }
84 //____________________________________________________________
85 AliHFEBeautySpectrumQA::AliHFEBeautySpectrumQA(const char *name):
86   TNamed(name, ""),
87   fPtMax(8.0),
88   fListOfResult(),
89   fWriteToFile(kTRUE)
90   {
91   //
92   // Default constructor
93   //
94
95   fListOfResult = new TObjArray(kNResults);
96   fListOfResult->SetName("ListOfResults");
97  
98
99 }
100
101 //____________________________________________________________
102 AliHFEBeautySpectrumQA::~AliHFEBeautySpectrumQA(){
103   //
104   // Destructor
105   //
106   if(fListOfResult) delete fListOfResult;
107  
108 }
109 //____________________________________________________________
110 void AliHFEBeautySpectrumQA::AddResultAt(TObject *obj,Int_t index)
111 {
112   //
113   // Init what we need for the correction:
114   //
115
116   if(fListOfResult) fListOfResult->AddAt(obj,index);
117
118 }
119 //____________________________________________________________
120 TObject *AliHFEBeautySpectrumQA::GetResult(Int_t index)
121 {
122   //
123   // Get result
124   //
125
126   if(fListOfResult) return fListOfResult->UncheckedAt(index);
127   else return 0x0;
128
129 }
130 //____________________________________________________________
131 void AliHFEBeautySpectrumQA::DrawProjections() const
132 {
133   //
134   // get spectrum for beauty 2nd method
135   //
136   //
137   AliCFContainer *data = (AliCFContainer *) fListOfResult->UncheckedAt(kDataProjection);
138   THnSparseF *correlation = (THnSparseF *) fListOfResult->UncheckedAt(kCMProjection);
139   if(!data || !correlation) return;
140
141   Int_t ndimcont = data->GetNVar();
142   Int_t ndimcor = correlation->GetNdimensions();
143   Int_t charge = 3;
144   Int_t centrality = 5;
145   Int_t eta = 1;
146
147   TCanvas * canvas = new TCanvas("Projections","Projections",1000,700);
148   Int_t n = 0;
149   if(charge < ndimcont) n++;
150   if(centrality < ndimcont) n++;
151   if(eta < ndimcont) n++;
152   canvas->Divide(2,n);
153   Int_t counter = 1;
154
155   if(charge < ndimcont) {
156    
157     canvas->cd(counter);
158     TH1 *checkcharge = (TH1 *) data->Project(data->GetNStep()-1,charge);
159     checkcharge->Draw();
160     counter++;
161     canvas->cd(counter);
162     TH2F* projectioncharge = (TH2F *) correlation->Projection(charge,charge+((Int_t)(ndimcor/2.)));
163     projectioncharge->Draw("colz");
164     counter++;  
165
166   }
167
168   if(centrality < ndimcont) {
169     canvas->cd(counter);
170     TH1 *checkcentrality = (TH1 *) data->Project(data->GetNStep()-1,centrality);
171     checkcentrality->Draw();
172     counter++;
173     canvas->cd(counter);
174     TH2F *projectioncentrality = (TH2F *) correlation->Projection(centrality,centrality+((Int_t)(ndimcor/2.)));
175     projectioncentrality->Draw("colz");
176     counter++;  
177   }
178
179   if(eta < ndimcont) {
180     canvas->cd(counter);
181     TH1 *checketa = (TH1 *) data->Project(data->GetNStep()-1,eta);
182     checketa->Draw();
183     counter++;
184     canvas->cd(counter);
185     TH2D* projectioneta = (TH2D *) correlation->Projection(eta,eta+((Int_t)(ndimcor/2.)));
186     projectioneta->Draw("colz");
187   }
188
189
190 }
191 //____________________________________________________________
192 void AliHFEBeautySpectrumQA::DrawSubtractContamination() const
193 {
194   //
195   // get spectrum for beauty 2nd method
196   //
197   //
198   TH1D *measuredTH1Daftersubstraction = (TH1D *) fListOfResult->UncheckedAt(kAfterSC);
199   TH1D *measuredTH1Dbeforesubstraction = (TH1D *) fListOfResult->UncheckedAt(kBeforeSC);
200   TH1D *measuredTH1background = (TH1D *) fListOfResult->UncheckedAt(kMeasBG);
201   if(!measuredTH1Daftersubstraction || !measuredTH1Dbeforesubstraction) return;
202
203   SetStyle();
204
205   TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
206   cbackgroundsubtraction->Divide(3,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   cbackgroundsubtraction->cd(3);
259   measuredTH1background->SetStats(0);
260   measuredTH1background->SetTitle("");
261   measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
262   measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]");
263   measuredTH1background->SetMarkerStyle(26);
264   measuredTH1background->SetMarkerColor(kBlack);
265   measuredTH1background->SetLineColor(kBlack);
266   measuredTH1background->Draw();
267   if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.png");
268
269 }
270
271
272 //____________________________________________________________
273 void AliHFEBeautySpectrumQA::DrawCorrectWithEfficiency(Int_t typeeff) const
274 {
275   //
276   // Correct the spectrum for efficiency and unfolding
277   // with both method and compare
278   //
279   
280   TH1D *afterE = 0x0;
281   TH1D *beforeE = 0x0;
282   TH1D *efficiencyDproj = 0x0;
283   TF1 *efficiencyparametrized = 0x0;
284   
285   if(typeeff== kMC) {
286     afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
287     beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforeMCE);
288     efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kMCEfficiency);
289   }
290   if(typeeff== kParametrized) {
291     afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterPE);
292     beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforePE);
293     efficiencyparametrized = (TF1 *) fListOfResult->UncheckedAt(kPEfficiency);
294   }
295   
296   if(typeeff==kMC && (!afterE || !beforeE || !efficiencyDproj)) return;
297   if(typeeff==kParametrized && (!afterE || !beforeE || !efficiencyparametrized)) return;
298   
299   SetStyle();
300   
301   TCanvas * cEfficiency = new TCanvas(AliHFEBeautySpectrumQA::fgkNameCanvas[typeeff],AliHFEBeautySpectrumQA::fgkNameCanvas[typeeff],1000,700);
302   cEfficiency->Divide(2,1);
303   cEfficiency->cd(1);
304   gPad->SetLogy();
305   gPad->SetTicks();
306   afterE->SetStats(0);
307   afterE->SetTitle("");
308   afterE->GetYaxis()->SetTitleOffset(1.5);
309   afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
310   afterE->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
311   afterE->GetXaxis()->SetRangeUser(0.0,fPtMax);
312   afterE->SetMarkerStyle(25);
313   afterE->SetMarkerColor(kBlack);
314   afterE->SetLineColor(kBlack);
315   beforeE->SetStats(0);
316   beforeE->SetTitle("");
317   beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
318   beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
319   beforeE->GetXaxis()->SetRangeUser(0.0,fPtMax);
320   beforeE->SetMarkerStyle(24);
321   beforeE->SetMarkerColor(kBlue);
322   beforeE->SetLineColor(kBlue);
323   afterE->Draw();
324   beforeE->Draw("same");
325   TLegend *legefficiency = new TLegend(0.4,0.6,0.89,0.89);
326   legefficiency->AddEntry(beforeE,"Before Efficiency correction","p");
327   legefficiency->AddEntry(afterE,"After Efficiency correction","p");
328   legefficiency->SetFillStyle(0);
329   legefficiency->SetLineStyle(0);
330   legefficiency->SetLineColor(0);
331   legefficiency->Draw("same");
332   cEfficiency->cd(2);
333   gPad->SetTicks();
334   if(typeeff==kMC) {
335     if(efficiencyDproj) {
336       efficiencyDproj->SetTitle("");
337       efficiencyDproj->SetStats(0);
338       efficiencyDproj->GetYaxis()->SetTitleOffset(1.5);
339       efficiencyDproj->GetYaxis()->SetRangeUser(0.0,1.0);
340       efficiencyDproj->GetYaxis()->SetTitle("Efficiency");
341       efficiencyDproj->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
342       efficiencyDproj->GetXaxis()->SetRangeUser(0.0,fPtMax);
343       efficiencyDproj->SetMarkerStyle(25);
344       efficiencyDproj->Draw();
345     }
346   }
347   if(typeeff==kParametrized) {
348     if(efficiencyparametrized) efficiencyparametrized->Draw();
349   }
350   
351   if(fWriteToFile) {
352     if(typeeff==kMC) cEfficiency->SaveAs("EfficiencyMC.png");
353     if(typeeff==kParametrized) cEfficiency->SaveAs("EfficiencyParametrized.png");
354   }
355
356 }
357
358 //____________________________________________________________
359 void AliHFEBeautySpectrumQA::DrawUnfolding() const
360 {
361   //
362   // Draw unfolding
363   //
364   TH1D *measuredspectrumD = (TH1D *) fListOfResult->UncheckedAt(kBeforeU);
365   TH1D *residualspectrumD = (TH1D *) fListOfResult->UncheckedAt(kResidualU);
366   TH1D *efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kUEfficiency);
367   THnSparseF *correlation = (THnSparseF *) fListOfResult->UncheckedAt(kCMProjection);
368   
369   if(!measuredspectrumD || !residualspectrumD || !efficiencyDproj || !correlation) return;
370   
371   Int_t ndimcor = (Int_t) correlation->GetNdimensions()/2.;
372
373   SetStyle();
374
375   TCanvas * cunfolding = new TCanvas("unfolding","unfolding",1000,700);
376   cunfolding->Divide(2,2);
377   cunfolding->cd(1);
378   gPad->SetLogy();
379   gPad->SetTicks();
380   residualspectrumD->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
381   residualspectrumD->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
382   residualspectrumD->GetXaxis()->SetRangeUser(0.0,fPtMax);
383   residualspectrumD->SetStats(0);
384   residualspectrumD->SetTitle("");
385   residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
386   residualspectrumD->SetMarkerStyle(26);
387   residualspectrumD->SetMarkerColor(kBlue);
388   residualspectrumD->SetLineColor(kBlue);
389   residualspectrumD->Sumw2();
390   residualspectrumD->Draw("P");
391   measuredspectrumD->SetStats(0);
392   measuredspectrumD->SetTitle("");  
393   measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
394   measuredspectrumD->SetMarkerStyle(25);
395   measuredspectrumD->SetMarkerColor(kBlack);
396   measuredspectrumD->SetLineColor(kBlack);
397   measuredspectrumD->Draw("same");
398   TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
399   legres->AddEntry(residualspectrumD,"Residual","p");
400   legres->AddEntry(measuredspectrumD,"Measured","p");
401   legres->SetFillStyle(0);
402   legres->SetLineStyle(0);
403   legres->SetLineColor(0);
404   legres->Draw("same");
405   cunfolding->cd(2);
406   gPad->SetTicks();
407   TH1D* ratioresidual = (TH1D*)residualspectrumD->Clone();
408   ratioresidual->SetName("ratioresidual");
409   ratioresidual->SetTitle("");
410   ratioresidual->GetYaxis()->SetRangeUser(0.6,1.4);
411   ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
412   ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
413   ratioresidual->Divide(measuredspectrumD);
414   ratioresidual->SetStats(0);
415   ratioresidual->Draw();
416   cunfolding->cd(3);
417   gPad->SetTicks();
418   efficiencyDproj->GetYaxis()->SetTitle("Efficiency");
419   efficiencyDproj->GetXaxis()->SetTitle("p^{MC}_{T} [GeV/c]");
420   efficiencyDproj->GetXaxis()->SetRangeUser(0.0,fPtMax);
421   efficiencyDproj->GetYaxis()->SetRangeUser(0.0,1.0);
422   efficiencyDproj->SetStats(0);
423   efficiencyDproj->SetTitle("");
424   efficiencyDproj->GetYaxis()->SetTitleOffset(1.5);
425   efficiencyDproj->SetMarkerStyle(26);
426   efficiencyDproj->SetMarkerColor(kBlue);
427   efficiencyDproj->SetLineColor(kBlue);
428   efficiencyDproj->Sumw2();
429   efficiencyDproj->Draw("P");
430   cunfolding->cd(4);
431   TH2F *projectioncorr = (TH2F *) correlation->Projection(0,ndimcor);
432   projectioncorr->GetYaxis()->SetTitle("p^{ESD}_{T} [GeV/c]");
433   projectioncorr->GetXaxis()->SetTitle("p^{MC}_{T} [GeV/c]");
434   projectioncorr->GetXaxis()->SetRangeUser(0.0,fPtMax);
435   projectioncorr->GetYaxis()->SetRangeUser(0.0,fPtMax);
436   projectioncorr->SetStats(0);
437   projectioncorr->SetTitle("");
438   projectioncorr->Draw("colz");
439
440   if(fWriteToFile){
441     cunfolding->SaveAs("Unfolding.png");
442   }
443   
444 }
445 //____________________________________________________________
446 void AliHFEBeautySpectrumQA::DrawResult()
447 {
448   //
449   // Draw Results
450   //
451   TGraphErrors* correctedspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt(kFinalResultUnfolded);
452   TGraphErrors* alltogetherspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt(kFinalResultDirectEfficiency);
453   THnSparse* correctedspectrum = (THnSparse *) fListOfResult->UncheckedAt(kFinalResultUnfSparse);
454   AliCFDataGrid* alltogetherCorrection = (AliCFDataGrid *) fListOfResult->UncheckedAt( kFinalResultDirectEffSparse);
455   if(!correctedspectrumD || !alltogetherspectrumD) return;
456   
457   SetStyle();
458
459   TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
460   ccorrected->Divide(2,1);
461   ccorrected->cd(1);
462   gPad->SetLogy();
463   correctedspectrumD->SetTitle("");
464   correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
465   correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
466   correctedspectrumD->SetMarkerStyle(26);
467   correctedspectrumD->SetMarkerColor(kBlue);
468   correctedspectrumD->SetLineColor(kBlue);
469   correctedspectrumD->Draw("AP");
470   alltogetherspectrumD->SetTitle("");
471   alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
472   alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
473   alltogetherspectrumD->SetMarkerStyle(25);
474   alltogetherspectrumD->SetMarkerColor(kBlack);
475   alltogetherspectrumD->SetLineColor(kBlack);
476   alltogetherspectrumD->Draw("P");
477   TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
478   legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
479   legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
480   legcorrected->Draw("same");
481   ccorrected->cd(2);
482   TH1D *correctedTH1D = correctedspectrum->Projection(0);
483   TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(0);
484   TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
485   ratiocorrected->SetName("ratiocorrected");
486   ratiocorrected->SetTitle("");
487   ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
488   ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
489   ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
490   ratiocorrected->SetStats(0);
491   ratiocorrected->Draw();
492   if(fWriteToFile){ 
493     ccorrected->SaveAs("CorrectedBeauty.eps");    
494     TFile *out;
495     out = new TFile("finalSpectrum.root","recreate");
496     out->cd();
497     //
498     correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
499     correctedspectrumD->Write();
500     alltogetherspectrumD->SetName("AlltogetherSpectrum");
501     alltogetherspectrumD->Write();
502     ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
503     ratiocorrected->Write();
504     //
505  
506     out->Close(); 
507     delete out;
508   }
509 }
510 //____________________________________________________________
511 TH1D *AliHFEBeautySpectrumQA::DivideSpectra(TGraphErrors *ga, TGraphErrors *gb) 
512 {
513   //
514   // Divide Spectra
515   //
516
517   TH1D *afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
518   if(!afterE) return 0x0;
519
520   TH1D *histoB = (TH1D*) afterE->Clone();
521   histoB->Sumw2();
522   histoB->SetName("ratio");
523   TH1D *histoa = (TH1D*) afterE->Clone();
524   histoa->Sumw2();
525   histoa->SetName("a");
526   TH1D *histob = (TH1D*) afterE->Clone();
527   histob->Sumw2();
528   histob->SetName("b");
529   
530   double xa,ya,xb,yb,eya,eyb;
531   Int_t npointsa = ga->GetN();
532   Int_t npointsb = gb->GetN();
533   if(npointsa != npointsb) {
534     printf("Problem the two spectra have not the same number of points\n");
535     return 0x0;
536   }
537   for(Int_t k = 0; k < npointsa; k++){
538     ga->GetPoint(k,xa,ya);
539     gb->GetPoint(k,xb,yb);
540     //
541     Double_t centerhisto = histoa->GetBinCenter(k+1);
542     //
543     if((TMath::Abs(xa-xb) > 0.0001) || (TMath::Abs(xa-centerhisto) > 0.0001)) {
544       printf("Problem not the same x axis\n");
545       return 0x0;
546     }
547     histoa->SetBinContent(k+1,ya);
548     histob->SetBinContent(k+1,yb);
549     //
550     eya = ga->GetErrorY(k);
551     eyb = gb->GetErrorY(k);
552     //
553     histoa->SetBinError(k+1,eya);
554     histob->SetBinError(k+1,eyb);
555    
556   }
557   
558   histoB->Sumw2();
559   histoB->Divide(histoa,histob,1.0,1.0,"B");
560
561   return histoB;  
562
563 }
564 //__________________________________________
565 void AliHFEBeautySpectrumQA::SetStyle() const
566 {
567   //
568   // Set style
569   //
570
571   gStyle->SetPalette(1);
572   gStyle->SetOptStat(1111);
573   gStyle->SetPadBorderMode(0);
574   gStyle->SetCanvasColor(10);
575   gStyle->SetPadLeftMargin(0.13);
576   gStyle->SetPadRightMargin(0.13);
577
578 }