]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliHFEspectrum.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added pass1 and pass2 directories
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEspectrum.cxx
CommitLineData
c04c80e6 1/**************************************************************************
2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
7* Permission to use, copy, modify and distribute this software and its *
8* documentation strictly for non-commercial purposes is hereby granted *
9* without fee, provided that the above copyright notice appears in all *
10* copies and that both the copyright notice and this permission notice *
11* appear in the supporting documentation. The authors make no claims *
12* about the suitability of this software for any purpose. It is *
13* provided "as is" without express or implied warranty. *
14**************************************************************************/
27de2dfb 15
16/* $Id$ */
17
c04c80e6 18//
19// Class for spectrum correction
20// Subtraction of hadronic background, Unfolding of the data and
21// Renormalization done here
22// The following containers have to be set:
23// - Correction framework container for real data
24// - Correction framework container for MC (Efficiency Map)
25// - Correction framework container for background coming from data
26// - Correction framework container for background coming from MC
27//
245877d0 28// Author:
29// Raphaelle Bailhache <R.Bailhache@gsi.de>
30// Markus Fasel <M.Fasel@gsi.de>
c04c80e6 31//
32
33#include <TArrayD.h>
34#include <TH1.h>
35#include <TList.h>
36#include <TObjArray.h>
37#include <TROOT.h>
38#include <TCanvas.h>
39#include <TLegend.h>
40#include <TStyle.h>
41#include <TMath.h>
42#include <TAxis.h>
43#include <TGraphErrors.h>
44#include <TFile.h>
45#include <TPad.h>
67fe7bd0 46#include <TH2D.h>
c04c80e6 47
3a72645a 48#include "AliPID.h"
c04c80e6 49#include "AliCFContainer.h"
50#include "AliCFDataGrid.h"
51#include "AliCFEffGrid.h"
52#include "AliCFGridSparse.h"
53#include "AliCFUnfolding.h"
54#include "AliLog.h"
55
56#include "AliHFEspectrum.h"
3a72645a 57#include "AliHFEcuts.h"
58#include "AliHFEcontainer.h"
c04c80e6 59
60ClassImp(AliHFEspectrum)
61
62//____________________________________________________________
63AliHFEspectrum::AliHFEspectrum(const char *name):
64 TNamed(name, ""),
65 fCFContainers(NULL),
66 fTemporaryObjects(NULL),
67 fCorrelation(NULL),
68 fBackground(NULL),
3a72645a 69 fInclusiveSpectrum(kTRUE),
c04c80e6 70 fDumpToFile(kFALSE),
3a72645a 71 fNbDimensions(1),
c04c80e6 72 fNEvents(0),
73 fStepMC(-1),
74 fStepTrue(-1),
75 fStepData(-1),
3a72645a 76 fStepBeforeCutsV0(-1),
77 fStepAfterCutsV0(-1),
c04c80e6 78 fStepGuessedUnfolding(-1),
3a72645a 79 fNumberOfIterations(5),
80 fDebugLevel(0)
c04c80e6 81{
82 //
83 // Default constructor
84 //
85}
86
87//____________________________________________________________
88AliHFEspectrum::~AliHFEspectrum(){
89 //
90 // Destructor
91 //
92 if(fCFContainers) delete fCFContainers;
93 if(fTemporaryObjects){
94 fTemporaryObjects->Clear();
95 delete fTemporaryObjects;
96 }
97}
3a72645a 98//____________________________________________________________
6555e2ad 99Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *v0hfecontainer){
3a72645a 100 //
101 // Init what we need for the correction:
102 //
103 // Raw spectrum, hadron contamination
104 // MC efficiency maps, correlation matrix
105 // V0 efficiency if wanted
106 //
107 // This for a given dimension.
108 // If no inclusive spectrum, then take only efficiency map for beauty electron
109 // and the appropriate correlation matrix
110 //
111
112 // Get the requested format
113 Int_t dims[3];
114 switch(fNbDimensions){
115 case 1: dims[0] = 0;
116 break;
117 case 2: for(Int_t i = 0; i < 2; i++) dims[i] = i;
118 break;
119 case 3: for(Int_t i = 0; i < 3; i++) dims[i] = i;
120 break;
121 default:
122 AliError("Container with this number of dimensions not foreseen (yet)");
123 return kFALSE;
124 };
125
126 // Data container: raw spectrum + hadron contamination
127 AliCFContainer *datacontainer = datahfecontainer->GetCFContainer("recTrackContReco");
128 AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
129 if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
130
131 AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims);
132 AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims);
133 if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
134 SetContainer(datacontainerD,AliHFEspectrum::kDataContainer);
135 SetContainer(contaminationcontainerD,AliHFEspectrum::kBackgroundData);
136
137 // MC container: ESD/MC efficiency maps + MC/MC efficiency maps
138 AliCFContainer *mccontaineresd = 0x0;
139 AliCFContainer *mccontainermc = 0x0;
140 if(fInclusiveSpectrum) {
141 mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco");
142 mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC");
143 }
144 else {
145 mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
146 mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
147 }
148 if((!mccontaineresd) || (!mccontainermc)) return kFALSE;
149
150 Int_t source = -1;
151 if(!fInclusiveSpectrum) source = 1;
152 AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source);
153 AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source);
154 if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
155 SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerMC);
156 SetContainer(mccontaineresdD,AliHFEspectrum::kMCContainerESD);
157
158 // MC container: correlation matrix
159 THnSparseF *mccorrelation = 0x0;
bf892a6a 160 if(fInclusiveSpectrum) {
161 if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
162 if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterTOF");
163 if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
164 }
3a72645a 165 else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
166 if(!mccorrelation) return kFALSE;
167 THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims);
168 if(!mccorrelationD) {
169 printf("No correlation\n");
170 return kFALSE;
171 }
172 SetCorrelation(mccorrelationD);
173
174 // V0 container Electron, pt eta phi
175 if(v0hfecontainer) {
176 AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
177 if(!containerV0) return kFALSE;
178 AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron+1);
179 if(!containerV0Electron) return kFALSE;
180 SetContainer(containerV0Electron,AliHFEspectrum::kDataContainerV0);
181 }
182
183 if(fDebugLevel>0){
184
185 AliCFDataGrid *contaminationspectrum = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(contaminationcontainer,1))->Clone();
186 contaminationspectrum->SetName("contaminationspectrum");
187 TCanvas * ccontaminationspectrum = new TCanvas("contaminationspectrum","contaminationspectrum",1000,700);
188 ccontaminationspectrum->Divide(3,1);
189 ccontaminationspectrum->cd(1);
6555e2ad 190 TH2D * contaminationspectrum2dpteta = (TH2D *) contaminationspectrum->Project(1,0);
191 TH2D * contaminationspectrum2dptphi = (TH2D *) contaminationspectrum->Project(2,0);
192 TH2D * contaminationspectrum2detaphi = (TH2D *) contaminationspectrum->Project(1,2);
3a72645a 193 contaminationspectrum2dpteta->SetStats(0);
194 contaminationspectrum2dpteta->SetTitle("");
195 contaminationspectrum2dpteta->GetXaxis()->SetTitle("#eta");
196 contaminationspectrum2dpteta->GetYaxis()->SetTitle("p_{T} [GeV/c]");
197 contaminationspectrum2dptphi->SetStats(0);
198 contaminationspectrum2dptphi->SetTitle("");
199 contaminationspectrum2dptphi->GetXaxis()->SetTitle("#phi [rad]");
200 contaminationspectrum2dptphi->GetYaxis()->SetTitle("p_{T} [GeV/c]");
201 contaminationspectrum2detaphi->SetStats(0);
202 contaminationspectrum2detaphi->SetTitle("");
203 contaminationspectrum2detaphi->GetXaxis()->SetTitle("#eta");
204 contaminationspectrum2detaphi->GetYaxis()->SetTitle("#phi [rad]");
205 contaminationspectrum2dptphi->Draw("colz");
206 ccontaminationspectrum->cd(2);
207 contaminationspectrum2dpteta->Draw("colz");
208 ccontaminationspectrum->cd(3);
209 contaminationspectrum2detaphi->Draw("colz");
210
211 }
212
213
214 return kTRUE;
215
216
217}
c04c80e6 218
219//____________________________________________________________
3a72645a 220Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
c04c80e6 221 //
222 // Correct the spectrum for efficiency and unfolding
223 // with both method and compare
224 //
225
226 gStyle->SetPalette(1);
227 gStyle->SetOptStat(1111);
228 gStyle->SetPadBorderMode(0);
229 gStyle->SetCanvasColor(10);
230 gStyle->SetPadLeftMargin(0.13);
231 gStyle->SetPadRightMargin(0.13);
67fe7bd0 232
3a72645a 233 ///////////////////////////
234 // Check initialization
235 ///////////////////////////
c04c80e6 236
3a72645a 237 if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
238 AliInfo("You have to init before");
239 return kFALSE;
240 }
241
242 if((fStepTrue == 0) && (fStepMC == 0) && (fStepData == 0)) {
243 AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
244 return kFALSE;
245 }
246
247 SetNumberOfIteration(50);
248 SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
249
250 AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
251 //////////////////////////////////
252 // Subtract hadron background
253 /////////////////////////////////
67fe7bd0 254 AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
3a72645a 255 if(subtractcontamination) {
256 dataspectrumaftersubstraction = SubtractBackground(kTRUE);
257 dataGridAfterFirstSteps = dataspectrumaftersubstraction;
258 }
259
260 ////////////////////////////////////////////////
261 // Correct for TPC efficiency from V0
262 ///////////////////////////////////////////////
263 AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
264 AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
265 if(dataContainerV0){
266 dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
267 dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
67fe7bd0 268 }
c04c80e6 269
3a72645a 270 ///////////////
c04c80e6 271 // Unfold
3a72645a 272 //////////////
273 TList *listunfolded = Unfold(dataGridAfterFirstSteps);
c04c80e6 274 if(!listunfolded){
275 printf("Unfolded failed\n");
3a72645a 276 return kFALSE;
c04c80e6 277 }
278 THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
279 THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
280 if(!correctedspectrum){
281 AliError("No corrected spectrum\n");
3a72645a 282 return kFALSE;
c04c80e6 283 }
67fe7bd0 284 if(!residualspectrum){
c04c80e6 285 AliError("No residul spectrum\n");
3a72645a 286 return kFALSE;
c04c80e6 287 }
67fe7bd0 288
3a72645a 289 /////////////////////
c04c80e6 290 // Simply correct
3a72645a 291 ////////////////////
292 AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
67fe7bd0 293
3a72645a 294
67fe7bd0 295 //////////
c04c80e6 296 // Plot
297 //////////
3a72645a 298 if(fDebugLevel > 0.0) {
299
300 TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
301 ccorrected->Divide(2,1);
302 ccorrected->cd(1);
303 gPad->SetLogy();
304 TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
305 correctedspectrumD->SetTitle("");
306 correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
307 correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
308 correctedspectrumD->SetMarkerStyle(26);
309 correctedspectrumD->SetMarkerColor(kBlue);
310 correctedspectrumD->SetLineColor(kBlue);
311 correctedspectrumD->Draw("AP");
312 TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
313 alltogetherspectrumD->SetTitle("");
314 alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
315 alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
316 alltogetherspectrumD->SetMarkerStyle(25);
317 alltogetherspectrumD->SetMarkerColor(kBlack);
318 alltogetherspectrumD->SetLineColor(kBlack);
319 alltogetherspectrumD->Draw("P");
320 TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
321 legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
322 legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
323 legcorrected->Draw("same");
324 ccorrected->cd(2);
325 TH1D *correctedTH1D = correctedspectrum->Projection(0);
6555e2ad 326 TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(0);
3a72645a 327 TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
328 ratiocorrected->SetName("ratiocorrected");
329 ratiocorrected->SetTitle("");
330 ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
331 ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
332 ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
333 ratiocorrected->SetStats(0);
334 ratiocorrected->Draw();
335
336
337 // Dump to file if needed
338
339 if(fDumpToFile) {
340
341 TFile *out = new TFile("finalSpectrum.root","recreate");
342 out->cd();
343 //
344 correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
345 correctedspectrumD->Write();
346 alltogetherspectrumD->SetName("AlltogetherSpectrum");
347 alltogetherspectrumD->Write();
348 ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
349 ratiocorrected->Write();
350 //
351 correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
352 correctedspectrum->Write();
353 alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
354 alltogetherCorrection->Write();
355 //
356 out->Close(); delete out;
357 }
358
359
360 }
361
362
363
364
365 return kTRUE;
366}
367//____________________________________________________________
368AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
369 //
370 // Apply background subtraction
371 //
372
373 // Raw spectrum
374 AliCFContainer *dataContainer = GetContainer(kDataContainer);
375 if(!dataContainer){
376 AliError("Data Container not available");
377 return NULL;
378 }
379 AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
380
381 AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
382 dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
383
384 // Background Estimate
385 AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
386 if(!backgroundContainer){
387 AliError("MC background container not found");
388 return NULL;
389 }
390
391 Int_t stepbackground = 1;
392 if(!fInclusiveSpectrum) stepbackground = 2;
393 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
394
395 // Subtract
396 spectrumSubtracted->Add(backgroundGrid,-1.0);
397 if(setBackground){
398 if(fBackground) delete fBackground;
399 fBackground = backgroundGrid;
400 } else delete backgroundGrid;
401
402
403 if(fDebugLevel > 0) {
67fe7bd0 404
405 TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
3a72645a 406 cbackgroundsubtraction->Divide(3,1);
67fe7bd0 407 cbackgroundsubtraction->cd(1);
3a72645a 408 gPad->SetLogy();
6555e2ad 409 TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(0);
410 TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(0);
67fe7bd0 411 CorrectFromTheWidth(measuredTH1Daftersubstraction);
412 CorrectFromTheWidth(measuredTH1Dbeforesubstraction);
413 measuredTH1Daftersubstraction->SetStats(0);
414 measuredTH1Daftersubstraction->SetTitle("");
415 measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
416 measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
417 measuredTH1Daftersubstraction->SetMarkerStyle(25);
418 measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
419 measuredTH1Daftersubstraction->SetLineColor(kBlack);
420 measuredTH1Dbeforesubstraction->SetStats(0);
421 measuredTH1Dbeforesubstraction->SetTitle("");
422 measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
423 measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
424 measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
425 measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
426 measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
427 measuredTH1Daftersubstraction->Draw();
428 measuredTH1Dbeforesubstraction->Draw("same");
429 TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
3a72645a 430 legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
431 legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
67fe7bd0 432 legsubstraction->Draw("same");
433 cbackgroundsubtraction->cd(2);
3a72645a 434 gPad->SetLogy();
67fe7bd0 435 TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
436 ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
437 ratiomeasuredcontamination->SetTitle("");
438 ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
439 ratiomeasuredcontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
440 ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
441 ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
442 ratiomeasuredcontamination->SetStats(0);
443 ratiomeasuredcontamination->SetMarkerStyle(26);
444 ratiomeasuredcontamination->SetMarkerColor(kBlack);
445 ratiomeasuredcontamination->SetLineColor(kBlack);
446 ratiomeasuredcontamination->Draw();
3a72645a 447 cbackgroundsubtraction->cd(3);
6555e2ad 448 TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(0);
3a72645a 449 CorrectFromTheWidth(measuredTH1background);
450 measuredTH1background->SetStats(0);
451 measuredTH1background->SetTitle("");
452 measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
453 measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]");
454 measuredTH1background->SetMarkerStyle(26);
455 measuredTH1background->SetMarkerColor(kBlack);
456 measuredTH1background->SetLineColor(kBlack);
457 measuredTH1background->Draw();
458
67fe7bd0 459 }
460
c04c80e6 461
3a72645a 462 return spectrumSubtracted;
c04c80e6 463}
c04c80e6 464//____________________________________________________________
3a72645a 465AliCFDataGrid *AliHFEspectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
466
c04c80e6 467 //
3a72645a 468 // Apply TPC pid efficiency correction from V0
c04c80e6 469 //
470
3a72645a 471 AliCFContainer *v0Container = GetContainer(kDataContainerV0);
472 if(!v0Container){
473 AliError("V0 Container not available");
c04c80e6 474 return NULL;
475 }
3a72645a 476
477 // Efficiency
478 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
479 efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
480
481 // Data in the right format
482 AliCFDataGrid *dataGrid = 0x0;
483 if(bgsubpectrum) {
484 dataGrid = bgsubpectrum;
c04c80e6 485 }
3a72645a 486 else {
487
488 AliCFContainer *dataContainer = GetContainer(kDataContainer);
489 if(!dataContainer){
490 AliError("Data Container not available");
491 return NULL;
492 }
c04c80e6 493
3a72645a 494 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
495 }
c04c80e6 496
3a72645a 497 // Correct
498 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
499 result->ApplyEffCorrection(*efficiencyD);
c04c80e6 500
3a72645a 501 if(fDebugLevel > 0) {
502
503 TCanvas * cV0Efficiency = new TCanvas("V0Efficiency","V0Efficiency",1000,700);
504 cV0Efficiency->Divide(2,1);
505 cV0Efficiency->cd(1);
6555e2ad 506 TH1D *afterE = (TH1D *) result->Project(0);
507 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
3a72645a 508 afterE->SetStats(0);
509 afterE->SetTitle("");
510 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
511 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
512 afterE->SetMarkerStyle(25);
513 afterE->SetMarkerColor(kBlack);
514 afterE->SetLineColor(kBlack);
515 beforeE->SetStats(0);
516 beforeE->SetTitle("");
517 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
518 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
519 beforeE->SetMarkerStyle(24);
520 beforeE->SetMarkerColor(kBlue);
521 beforeE->SetLineColor(kBlue);
522 afterE->Draw();
523 beforeE->Draw("same");
524 TLegend *legV0efficiency = new TLegend(0.4,0.6,0.89,0.89);
525 legV0efficiency->AddEntry(beforeE,"Before Efficiency correction","p");
526 legV0efficiency->AddEntry(afterE,"After Efficiency correction","p");
527 legV0efficiency->Draw("same");
528 cV0Efficiency->cd(2);
6555e2ad 529 TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
3a72645a 530 efficiencyDproj->SetTitle("");
531 efficiencyDproj->SetStats(0);
532 efficiencyDproj->SetMarkerStyle(25);
533 efficiencyDproj->Draw();
c04c80e6 534
3a72645a 535
536 }
537
538
539 return result;
540
541}
c04c80e6 542//____________________________________________________________
3a72645a 543TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
c04c80e6 544
545 //
546 // Unfold and eventually correct for efficiency the bgsubspectrum
547 //
548
3a72645a 549 AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
c04c80e6 550 if(!mcContainer){
551 AliError("MC Container not available");
552 return NULL;
553 }
554
555 if(!fCorrelation){
556 AliError("No Correlation map available");
557 return NULL;
558 }
559
3a72645a 560 // Data
c04c80e6 561 AliCFDataGrid *dataGrid = 0x0;
562 if(bgsubpectrum) {
c04c80e6 563 dataGrid = bgsubpectrum;
c04c80e6 564 }
565 else {
566
567 AliCFContainer *dataContainer = GetContainer(kDataContainer);
568 if(!dataContainer){
569 AliError("Data Container not available");
570 return NULL;
571 }
572
3a72645a 573 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
c04c80e6 574 }
575
c04c80e6 576 // Guessed
3a72645a 577 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
c04c80e6 578 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
579
580 // Efficiency
3a72645a 581 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
c04c80e6 582 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
583
584 // Unfold
585
3a72645a 586 AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
c04c80e6 587 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
588 unfolding.UseSmoothing();
589 unfolding.Unfold();
590
591 // Results
592 THnSparse* result = unfolding.GetUnfolded();
593 THnSparse* residual = unfolding.GetEstMeasured();
594 TList *listofresults = new TList;
595 listofresults->SetOwner();
596 listofresults->AddAt((THnSparse*)result->Clone(),0);
597 listofresults->AddAt((THnSparse*)residual->Clone(),1);
598
3a72645a 599 if(fDebugLevel > 0) {
600
601 TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
602 cresidual->Divide(2,1);
603 cresidual->cd(1);
604 gPad->SetLogy();
605 TGraphErrors* residualspectrumD = Normalize(residual);
606 if(!residualspectrumD) {
607 AliError("Number of Events not set for the normalization");
608 return kFALSE;
609 }
610 residualspectrumD->SetTitle("");
611 residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
612 residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
613 residualspectrumD->SetMarkerStyle(26);
614 residualspectrumD->SetMarkerColor(kBlue);
615 residualspectrumD->SetLineColor(kBlue);
616 residualspectrumD->Draw("AP");
617 AliCFDataGrid *dataGridBis = (AliCFDataGrid *) (dataGrid->Clone());
618 dataGridBis->SetName("dataGridBis");
619 TGraphErrors* measuredspectrumD = Normalize(dataGridBis);
620 measuredspectrumD->SetTitle("");
621 measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
622 measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
623 measuredspectrumD->SetMarkerStyle(25);
624 measuredspectrumD->SetMarkerColor(kBlack);
625 measuredspectrumD->SetLineColor(kBlack);
626 measuredspectrumD->Draw("P");
627 TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
628 legres->AddEntry(residualspectrumD,"Residual","p");
629 legres->AddEntry(measuredspectrumD,"Measured","p");
630 legres->Draw("same");
631 cresidual->cd(2);
632 TH1D *residualTH1D = residual->Projection(0);
6555e2ad 633 TH1D *measuredTH1D = (TH1D *) dataGridBis->Project(0);
3a72645a 634 TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
635 ratioresidual->SetName("ratioresidual");
636 ratioresidual->SetTitle("");
637 ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
638 ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
639 ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
640 ratioresidual->SetStats(0);
641 ratioresidual->Draw();
642
643 }
644
c04c80e6 645 return listofresults;
646
647}
648
649//____________________________________________________________
3a72645a 650AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
c04c80e6 651
652 //
653 // Apply unfolding and efficiency correction together to bgsubspectrum
654 //
655
3a72645a 656 AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
c04c80e6 657 if(!mcContainer){
658 AliError("MC Container not available");
659 return NULL;
660 }
661
c04c80e6 662 // Efficiency
3a72645a 663 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
c04c80e6 664 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
665
666 // Data in the right format
667 AliCFDataGrid *dataGrid = 0x0;
668 if(bgsubpectrum) {
c04c80e6 669 dataGrid = bgsubpectrum;
c04c80e6 670 }
671 else {
3a72645a 672
c04c80e6 673 AliCFContainer *dataContainer = GetContainer(kDataContainer);
674 if(!dataContainer){
675 AliError("Data Container not available");
676 return NULL;
677 }
678
3a72645a 679 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
c04c80e6 680 }
681
682 // Correct
683 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
684 result->ApplyEffCorrection(*efficiencyD);
c04c80e6 685
3a72645a 686 if(fDebugLevel > 0) {
687
bf892a6a 688 printf("Step MC: %d\n",fStepMC);
689 printf("Step tracking: %d\n",AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
690 printf("Step MC true: %d\n",fStepTrue);
3a72645a 691 AliCFEffGrid *efficiencymcPID = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
692 AliCFEffGrid *efficiencymctrackinggeo = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack,fStepTrue);
693 AliCFEffGrid *efficiencymcall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,fStepTrue);
694
695 AliCFEffGrid *efficiencyesdall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerESD),fStepMC,fStepTrue);
696
697 TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
698 cefficiency->cd(1);
6555e2ad 699 TH1D* efficiencymcPIDD = (TH1D *) efficiencymcPID->Project(0);
3a72645a 700 efficiencymcPIDD->SetTitle("");
701 efficiencymcPIDD->SetStats(0);
702 efficiencymcPIDD->SetMarkerStyle(25);
703 efficiencymcPIDD->Draw();
6555e2ad 704 TH1D* efficiencymctrackinggeoD = (TH1D *) efficiencymctrackinggeo->Project(0);
3a72645a 705 efficiencymctrackinggeoD->SetTitle("");
706 efficiencymctrackinggeoD->SetStats(0);
707 efficiencymctrackinggeoD->SetMarkerStyle(26);
708 efficiencymctrackinggeoD->Draw("same");
6555e2ad 709 TH1D* efficiencymcallD = (TH1D *) efficiencymcall->Project(0);
3a72645a 710 efficiencymcallD->SetTitle("");
711 efficiencymcallD->SetStats(0);
712 efficiencymcallD->SetMarkerStyle(27);
713 efficiencymcallD->Draw("same");
6555e2ad 714 TH1D* efficiencyesdallD = (TH1D *) efficiencyesdall->Project(0);
3a72645a 715 efficiencyesdallD->SetTitle("");
716 efficiencyesdallD->SetStats(0);
717 efficiencyesdallD->SetMarkerStyle(24);
718 efficiencyesdallD->Draw("same");
719 TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
720 legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
721 legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
722 legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
723 legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
724 legeff->Draw("same");
725 }
726
c04c80e6 727 return result;
728
729}
3a72645a 730
c04c80e6 731//__________________________________________________________________________________
732TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum) const {
733 //
734 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
735 // Give the final pt spectrum to be compared
736 //
737
738 if(fNEvents > 0) {
739
740 TH1D* projection = spectrum->Projection(0);
741 CorrectFromTheWidth(projection);
742 TGraphErrors *graphError = NormalizeTH1(projection);
743 return graphError;
744
745 }
746
747 return 0x0;
748
749
750}
751//__________________________________________________________________________________
752TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum) const {
753 //
754 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
755 // Give the final pt spectrum to be compared
756 //
757 if(fNEvents > 0) {
758
6555e2ad 759 TH1D* projection = (TH1D *) spectrum->Project(0);
c04c80e6 760 CorrectFromTheWidth(projection);
761 TGraphErrors *graphError = NormalizeTH1(projection);
762 return graphError;
763
764 }
765
766 return 0x0;
767
768
769}
770//__________________________________________________________________________________
771TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input) const {
772 //
773 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
774 // Give the final pt spectrum to be compared
775 //
776 if(fNEvents > 0) {
777
778 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
779 Double_t p = 0, dp = 0; Int_t point = 1;
780 Double_t n = 0, dN = 0;
781 Double_t nCorr = 0, dNcorr = 0;
782 Double_t errdN = 0, errdp = 0;
783 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
784 point = ibin - input->GetXaxis()->GetFirst();
785 p = input->GetXaxis()->GetBinCenter(ibin);
786 dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
787 n = input->GetBinContent(ibin);
788 dN = input->GetBinError(ibin);
789
790 // New point
791 nCorr = 0.5 * 1/1.6 * 1./(Double_t)(fNEvents) * 1/(2. * TMath::Pi() * p) * n;
792 errdN = 1./(2. * TMath::Pi() * p);
793 errdp = 1./(2. * TMath::Pi() * p*p) * n;
794 dNcorr = 0.5 * 1/1.6 * 1./(Double_t)(fNEvents) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
795
796 spectrumNormalized->SetPoint(point, p, nCorr);
797 spectrumNormalized->SetPointError(point, dp, dNcorr);
798 }
799 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
800 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
801 spectrumNormalized->SetMarkerStyle(22);
802 spectrumNormalized->SetMarkerColor(kBlue);
803 spectrumNormalized->SetLineColor(kBlue);
804
805 return spectrumNormalized;
806
807 }
808
809 return 0x0;
810
811
812}
813//____________________________________________________________
814void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
815 //
816 // Set the container for a given step to the
817 //
818 if(!fCFContainers) fCFContainers = new TList;
819 fCFContainers->AddAt(cont, type);
820}
821
822//____________________________________________________________
823AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
824 //
825 // Get Correction Framework Container for given type
826 //
827 if(!fCFContainers) return NULL;
828 return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
829}
c04c80e6 830//____________________________________________________________
3a72645a 831AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source) {
c04c80e6 832 //
3a72645a 833 // Slice bin for a given source of electron
c04c80e6 834 //
835
836 Double_t *varMin = new Double_t[container->GetNVar()],
837 *varMax = new Double_t[container->GetNVar()];
838
839 Double_t *binLimits;
840 for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
841
842 binLimits = new Double_t[container->GetNBins(ivar)+1];
843 container->GetBinLimits(ivar,binLimits);
3a72645a 844 if((ivar == 4) && ((source>= 0) && (source<container->GetNBins(ivar)))) {
845 varMin[ivar] = binLimits[source];
846 varMax[ivar] = binLimits[source];
847 }
848 else {
849 varMin[ivar] = binLimits[0];
850 varMax[ivar] = binLimits[container->GetNBins(ivar)];
851 }
852
c04c80e6 853 delete[] binLimits;
854 }
855
856 AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
857 AddTemporaryObject(k);
858 delete[] varMin; delete[] varMax;
859
860 return k;
861
862}
863
864//_________________________________________________________________________
3a72645a 865THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions) const {
c04c80e6 866 //
3a72645a 867 // Slice correlation
c04c80e6 868 //
869
3a72645a 870 Int_t ndimensions = correlationmatrix->GetNdimensions();
c04c80e6 871 printf("Number of dimension %d correlation map\n",ndimensions);
872 if(ndimensions < (2*nDim)) {
873 AliError("Problem in the dimensions");
874 return NULL;
875 }
876 Int_t ndimensionsContainer = (Int_t) ndimensions/2;
877 printf("Number of dimension %d container\n",ndimensionsContainer);
878
879 Int_t *dim = new Int_t[nDim*2];
880 for(Int_t iter=0; iter < nDim; iter++){
881 dim[iter] = dimensions[iter];
882 dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
883 printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
884 }
885
3a72645a 886 THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
c04c80e6 887
888 delete[] dim;
889 return k;
890
891}
892//___________________________________________________________________________
893void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
894 //
895 // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
896 //
897
898 TAxis *axis = h1->GetXaxis();
899 Int_t nbinX = h1->GetNbinsX();
900
901 for(Int_t i = 1; i <= nbinX; i++) {
902
903 Double_t width = axis->GetBinWidth(i);
904 Double_t content = h1->GetBinContent(i);
905 Double_t error = h1->GetBinError(i);
906 h1->SetBinContent(i,content/width);
907 h1->SetBinError(i,error/width);
908 }
909
910}
911//____________________________________________________________
912void AliHFEspectrum::AddTemporaryObject(TObject *o){
913 //
914 // Emulate garbage collection on class level
915 //
916 if(!fTemporaryObjects) {
917 fTemporaryObjects= new TList;
918 fTemporaryObjects->SetOwner();
919 }
920 fTemporaryObjects->Add(o);
921}
922
923//____________________________________________________________
924void AliHFEspectrum::ClearObject(TObject *o){
925 //
926 // Do a safe deletion
927 //
928 if(fTemporaryObjects){
929 if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
930 delete o;
931 } else{
932 // just delete
933 delete o;
934 }
935}
936//_________________________________________________________________________
937TObject* AliHFEspectrum::GetSpectrum(AliCFContainer * const c, Int_t step) {
245877d0 938 AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
c04c80e6 939 return data;
940}
941//_________________________________________________________________________
3a72645a 942TObject* AliHFEspectrum::GetEfficiency(AliCFContainer * const c, Int_t step, Int_t step0){
c04c80e6 943 //
944 // Create efficiency grid and calculate efficiency
945 // of step to step0
946 //
947 TString name("eff");
948 name += step;
949 name+= step0;
950 AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
951 eff->CalculateEfficiency(step,step0);
952 return eff;
953}