2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
27 // Raphaelle Bailhache <R.Bailhache@gsi.de>
33 #include <TObjArray.h>
41 #include <TGraphErrors.h>
48 #include "AliCFContainer.h"
49 #include "AliCFDataGrid.h"
50 #include "AliCFEffGrid.h"
51 #include "AliCFGridSparse.h"
52 #include "AliCFUnfolding.h"
55 #include "AliHFEInclusiveSpectrumQA.h"
56 #include "AliHFECorrectSpectrumBase.h"
57 #include "AliHFEcuts.h"
58 #include "AliHFEcontainer.h"
59 #include "AliHFEtools.h"
61 ClassImp(AliHFEInclusiveSpectrumQA)
63 const Char_t *AliHFEInclusiveSpectrumQA::fgkNameCanvas[AliHFEInclusiveSpectrumQA::kNTypeEfficiency] = {
66 "ParametrizedEfficiency"
69 //____________________________________________________________
70 AliHFEInclusiveSpectrumQA::AliHFEInclusiveSpectrumQA():
77 // Default constructor
80 fListOfResult = new TObjArray(kNResults);
81 fListOfResult->SetName("ListOfResults");
85 //____________________________________________________________
86 AliHFEInclusiveSpectrumQA::AliHFEInclusiveSpectrumQA(const char *name):
93 // Default constructor
96 fListOfResult = new TObjArray(kNResults);
97 fListOfResult->SetName("ListOfResults");
102 //____________________________________________________________
103 AliHFEInclusiveSpectrumQA::~AliHFEInclusiveSpectrumQA(){
107 if(fListOfResult) delete fListOfResult;
110 //____________________________________________________________
111 void AliHFEInclusiveSpectrumQA::AddResultAt(TObject *obj,Int_t index)
114 // Init what we need for the correction:
117 if(fListOfResult) fListOfResult->AddAt(obj,index);
120 //____________________________________________________________
121 TObject *AliHFEInclusiveSpectrumQA::GetResult(Int_t index)
127 if(fListOfResult) return fListOfResult->UncheckedAt(index);
131 //____________________________________________________________
132 void AliHFEInclusiveSpectrumQA::DrawProjections() const
135 // get spectrum for beauty 2nd method
138 AliCFContainer *data = (AliCFContainer *) fListOfResult->UncheckedAt(kDataProjection);
139 THnSparseF *correlation = (THnSparseF *) fListOfResult->UncheckedAt(kCMProjection);
140 if(!data || !correlation) return;
142 Int_t ndimcont = data->GetNVar();
143 Int_t ndimcor = correlation->GetNdimensions();
145 Int_t centrality = 5;
148 TCanvas * canvas = new TCanvas("Projections","Projections",1000,700);
150 if(charge < ndimcont) n++;
151 if(centrality < ndimcont) n++;
152 if(eta < ndimcont) n++;
156 if(charge < ndimcont) {
159 TH1 *checkcharge = (TH1 *) data->Project(data->GetNStep()-1,charge);
163 TH2F* projectioncharge = (TH2F *) correlation->Projection(charge,charge+((Int_t)(ndimcor/2.)));
164 projectioncharge->Draw("colz");
169 if(centrality < ndimcont) {
171 TH1 *checkcentrality = (TH1 *) data->Project(data->GetNStep()-1,centrality);
172 checkcentrality->Draw();
175 TH2F *projectioncentrality = (TH2F *) correlation->Projection(centrality,centrality+((Int_t)(ndimcor/2.)));
176 projectioncentrality->Draw("colz");
182 TH1 *checketa = (TH1 *) data->Project(data->GetNStep()-1,eta);
186 TH2D* projectioneta = (TH2D *) correlation->Projection(eta,eta+((Int_t)(ndimcor/2.)));
187 projectioneta->Draw("colz");
192 //____________________________________________________________
193 void AliHFEInclusiveSpectrumQA::DrawSubtractContamination() const
196 // get spectrum for beauty 2nd method
199 TH1D *measuredTH1Daftersubstraction = (TH1D *) fListOfResult->UncheckedAt(kAfterSC);
200 TH1D *measuredTH1Dbeforesubstraction = (TH1D *) fListOfResult->UncheckedAt(kBeforeSC);
201 if(!measuredTH1Daftersubstraction || !measuredTH1Dbeforesubstraction) return;
205 TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
206 cbackgroundsubtraction->Divide(2,1);
207 cbackgroundsubtraction->cd(1);
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);
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);
257 ratiomeasuredcontamination->Draw("P");
258 if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.png");
263 //____________________________________________________________
264 void AliHFEInclusiveSpectrumQA::DrawCorrectWithEfficiency(Int_t typeeff) const
267 // Correct the spectrum for efficiency and unfolding
268 // with both method and compare
273 TH1D *efficiencyDproj = 0x0;
274 TF1 *efficiencyparametrized = 0x0;
277 afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterV0);
278 beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforeV0);
279 efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kV0Efficiency);
282 afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
283 beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforeMCE);
284 efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kMCEfficiency);
286 if(typeeff== kParametrized) {
287 afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterPE);
288 beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforePE);
289 efficiencyparametrized = (TF1 *) fListOfResult->UncheckedAt(kPEfficiency);
292 if((typeeff==kV0 || typeeff==kMC) && (!afterE || !beforeE || !efficiencyDproj)) return;
293 if(typeeff==kParametrized && (!afterE || !beforeE || !efficiencyparametrized)) return;
297 TCanvas * cEfficiency = new TCanvas(AliHFEInclusiveSpectrumQA::fgkNameCanvas[typeeff],AliHFEInclusiveSpectrumQA::fgkNameCanvas[typeeff],1000,700);
298 cEfficiency->Divide(2,1);
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);
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");
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();
341 if(typeeff==kParametrized) {
342 efficiencyparametrized->Draw();
346 if(typeeff==kV0) cEfficiency->SaveAs("EfficiencyV0.png");
347 if(typeeff==kMC) cEfficiency->SaveAs("EfficiencyMC.png");
348 if(typeeff==kParametrized) cEfficiency->SaveAs("EfficiencyParametrized.png");
353 //____________________________________________________________
354 void AliHFEInclusiveSpectrumQA::DrawUnfolding() const
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);
364 if(!measuredspectrumD || !residualspectrumD || !efficiencyDproj || !correlation) return;
366 Int_t ndimcor = (Int_t) correlation->GetNdimensions()/2.;
370 TCanvas * cunfolding = new TCanvas("unfolding","unfolding",1000,700);
371 cunfolding->Divide(2,2);
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");
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();
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");
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");
436 cunfolding->SaveAs("Unfolding.png");
440 //____________________________________________________________
441 void AliHFEInclusiveSpectrumQA::DrawResult()
446 TGraphErrors* correctedspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt(kFinalResultUnfolded);
447 TGraphErrors* alltogetherspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt( kFinalResultDirectEfficiency);
448 if(!correctedspectrumD || !alltogetherspectrumD) return;
452 TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
453 ccorrected->Divide(2,1);
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");
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");
494 //____________________________________________________________
495 TH1D *AliHFEInclusiveSpectrumQA::DivideSpectra(TGraphErrors *ga, TGraphErrors *gb)
501 TH1D *afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
502 if(!afterE) return 0x0;
504 TH1D *histoB = (TH1D*) afterE->Clone();
506 histoB->SetName("ratio");
507 TH1D *histoa = (TH1D*) afterE->Clone();
509 histoa->SetName("a");
510 TH1D *histob = (TH1D*) afterE->Clone();
512 histob->SetName("b");
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");
521 for(Int_t k = 0; k < npointsa; k++){
522 ga->GetPoint(k,xa,ya);
523 gb->GetPoint(k,xb,yb);
525 Double_t centerhisto = histoa->GetBinCenter(k+1);
527 if((TMath::Abs(xa-xb) > 0.0001) || (TMath::Abs(xa-centerhisto) > 0.0001)) {
528 printf("Problem not the same x axis\n");
531 histoa->SetBinContent(k+1,ya);
532 histob->SetBinContent(k+1,yb);
534 eya = ga->GetErrorY(k);
535 eyb = gb->GetErrorY(k);
537 histoa->SetBinError(k+1,eya);
538 histob->SetBinError(k+1,eyb);
543 histoB->Divide(histoa,histob,1.0,1.0,"B");
548 //__________________________________________
549 void AliHFEInclusiveSpectrumQA::SetStyle() const
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);