]> 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
Some more clang warnings corrected
[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**************************************************************************/
15//
16// Class for spectrum correction
17// Subtraction of hadronic background, Unfolding of the data and
18// Renormalization done here
19// The following containers have to be set:
20// - Correction framework container for real data
21// - Correction framework container for MC (Efficiency Map)
22// - Correction framework container for background coming from data
23// - Correction framework container for background coming from MC
24//
245877d0 25// Author:
26// Raphaelle Bailhache <R.Bailhache@gsi.de>
27// Markus Fasel <M.Fasel@gsi.de>
c04c80e6 28//
29
30#include <TArrayD.h>
31#include <TH1.h>
32#include <TList.h>
33#include <TObjArray.h>
34#include <TROOT.h>
35#include <TCanvas.h>
36#include <TLegend.h>
37#include <TStyle.h>
38#include <TMath.h>
39#include <TAxis.h>
40#include <TGraphErrors.h>
41#include <TFile.h>
42#include <TPad.h>
67fe7bd0 43#include <TH2D.h>
c2690925 44#include <TF1.h>
c04c80e6 45
3a72645a 46#include "AliPID.h"
c04c80e6 47#include "AliCFContainer.h"
48#include "AliCFDataGrid.h"
49#include "AliCFEffGrid.h"
50#include "AliCFGridSparse.h"
51#include "AliCFUnfolding.h"
52#include "AliLog.h"
53
54#include "AliHFEspectrum.h"
3a72645a 55#include "AliHFEcuts.h"
56#include "AliHFEcontainer.h"
c2690925 57#include "AliHFEtools.h"
c04c80e6 58
59ClassImp(AliHFEspectrum)
60
61//____________________________________________________________
62AliHFEspectrum::AliHFEspectrum(const char *name):
63 TNamed(name, ""),
64 fCFContainers(NULL),
65 fTemporaryObjects(NULL),
66 fCorrelation(NULL),
67 fBackground(NULL),
c2690925 68 fEfficiencyFunction(NULL),
69 fWeightCharm(NULL),
3a72645a 70 fInclusiveSpectrum(kTRUE),
c04c80e6 71 fDumpToFile(kFALSE),
c2690925 72 fUnSetCorrelatedErrors(kTRUE),
e156c3bb 73 fSetSmoothing(kFALSE),
c2690925 74 fIPanaHadronBgSubtract(kFALSE),
75 fIPanaCharmBgSubtract(kFALSE),
76 fIPanaConversionBgSubtract(kFALSE),
77 fIPanaNonHFEBgSubtract(kFALSE),
78 fNonHFEbgMethod2(kFALSE),
3a72645a 79 fNbDimensions(1),
c2690925 80 fNMCEvents(0),
81 fNMCbgEvents(0),
c04c80e6 82 fStepMC(-1),
83 fStepTrue(-1),
84 fStepData(-1),
3a72645a 85 fStepBeforeCutsV0(-1),
86 fStepAfterCutsV0(-1),
c04c80e6 87 fStepGuessedUnfolding(-1),
3a72645a 88 fNumberOfIterations(5),
c2690925 89 fChargeChoosen(-1),
90 fNCentralityBinAtTheEnd(0),
91 fHadronEffbyIPcut(NULL),
92 fBeamType(0),
3a72645a 93 fDebugLevel(0)
c04c80e6 94{
95 //
96 // Default constructor
97 //
c2690925 98
99 for(Int_t k = 0; k < 20; k++){
100 fNEvents[k] = 0;
101 fLowBoundaryCentralityBinAtTheEnd[k] = 0;
102 fHighBoundaryCentralityBinAtTheEnd[k] = 0;
103 }
c04c80e6 104}
105
106//____________________________________________________________
107AliHFEspectrum::~AliHFEspectrum(){
108 //
109 // Destructor
110 //
111 if(fCFContainers) delete fCFContainers;
112 if(fTemporaryObjects){
113 fTemporaryObjects->Clear();
114 delete fTemporaryObjects;
115 }
116}
3a72645a 117//____________________________________________________________
c2690925 118Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *v0hfecontainer, const AliHFEcontainer *bghfecontainer){
3a72645a 119 //
120 // Init what we need for the correction:
121 //
122 // Raw spectrum, hadron contamination
123 // MC efficiency maps, correlation matrix
124 // V0 efficiency if wanted
125 //
126 // This for a given dimension.
127 // If no inclusive spectrum, then take only efficiency map for beauty electron
128 // and the appropriate correlation matrix
129 //
c2690925 130 Int_t kNdim;
131 if(fBeamType==0) kNdim=3;
132 if(fBeamType==1) kNdim=4;
133
134 Int_t dims[kNdim];
135 // Get the requested format
136 if(fBeamType==0)
137 {
138 // pp analysis
139 switch(fNbDimensions){
140 case 1: dims[0] = 0;
141 break;
142 case 2: for(Int_t i = 0; i < 2; i++) dims[i] = i;
143 break;
144 case 3: for(Int_t i = 0; i < 3; i++) dims[i] = i;
145 break;
146 default:
147 AliError("Container with this number of dimensions not foreseen (yet)");
148 return kFALSE;
149 };
150 }
151
152 if(fBeamType==1)
153 {
154 // PbPb analysis; centrality as first dimension
155 Int_t nbDimensions = fNbDimensions;
156 fNbDimensions = fNbDimensions + 1;
157 switch(nbDimensions){
158 case 1:
6c70d827 159 dims[0] = 5;
160 dims[1] = 0;
161 break;
c2690925 162 case 2:
6c70d827 163 dims[0] = 5;
164 for(Int_t i = 0; i < 2; i++) dims[(i+1)] = i;
165 break;
166 case 3:
167 dims[0] = 5;
168 for(Int_t i = 0; i < 3; i++) dims[(i+1)] = i;
169 break;
c2690925 170 default:
171 AliError("Container with this number of dimensions not foreseen (yet)");
172 return kFALSE;
173 };
174 }
175
176
3a72645a 177
3a72645a 178 // Data container: raw spectrum + hadron contamination
c2690925 179 AliCFContainer *datacontainer = 0x0;
180 if(fInclusiveSpectrum) {
181 datacontainer = datahfecontainer->GetCFContainer("recTrackContReco");
182 }
183 else{
184 datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
185 }
3a72645a 186 AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
187 if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
188
c2690925 189 AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen);
190 AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen);
3a72645a 191 if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
c2690925 192
3a72645a 193 SetContainer(datacontainerD,AliHFEspectrum::kDataContainer);
194 SetContainer(contaminationcontainerD,AliHFEspectrum::kBackgroundData);
195
c2690925 196 if(bghfecontainer){
197 // set nonHFE backgrounds using temporary stored in hadronicBackground, will be replaced with proper container
198 AliCFContainer *bgcontainer = bghfecontainer->GetCFContainer("hadronicBackground");
199 if(!bgcontainer) return kFALSE;
200 AliCFContainer *bgcontainerD = GetSlicedContainer(bgcontainer, fNbDimensions, dims, -1, fChargeChoosen);
201 if(!bgcontainerD) return kFALSE;
202 SetContainer(bgcontainerD,AliHFEspectrum::kNonHFEBackgroundData);
203 }
204
3a72645a 205 // MC container: ESD/MC efficiency maps + MC/MC efficiency maps
206 AliCFContainer *mccontaineresd = 0x0;
207 AliCFContainer *mccontainermc = 0x0;
208 if(fInclusiveSpectrum) {
209 mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco");
210 mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC");
211 }
212 else {
213 mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
214 mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
215 }
216 if((!mccontaineresd) || (!mccontainermc)) return kFALSE;
217
218 Int_t source = -1;
219 if(!fInclusiveSpectrum) source = 1;
c2690925 220 AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen);
221 AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen);
3a72645a 222 if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
223 SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerMC);
224 SetContainer(mccontaineresdD,AliHFEspectrum::kMCContainerESD);
225
c2690925 226 // set charm, nonHFE container to estimate BG
227 if(!fInclusiveSpectrum){
228 source = 0; //charm
229 mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen);
230 SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerCharmMC);
231 source = 2; //conversion
232 mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen);
233 SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerConversionMC);
234 source = 3; //non HFE except for the conversion
235 mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen);
236 SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerNonHFEMC);
237 }
238
3a72645a 239 // MC container: correlation matrix
240 THnSparseF *mccorrelation = 0x0;
bf892a6a 241 if(fInclusiveSpectrum) {
242 if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
6c70d827 243 else if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
244 else if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
245 else if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD - 1))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
246 else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
e156c3bb 247
248 if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
bf892a6a 249 }
c2690925 250 else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
251 //else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
3a72645a 252 if(!mccorrelation) return kFALSE;
253 THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims);
254 if(!mccorrelationD) {
255 printf("No correlation\n");
256 return kFALSE;
257 }
258 SetCorrelation(mccorrelationD);
259
260 // V0 container Electron, pt eta phi
261 if(v0hfecontainer) {
262 AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
263 if(!containerV0) return kFALSE;
264 AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron+1);
265 if(!containerV0Electron) return kFALSE;
266 SetContainer(containerV0Electron,AliHFEspectrum::kDataContainerV0);
267 }
268
269 if(fDebugLevel>0){
270
271 AliCFDataGrid *contaminationspectrum = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(contaminationcontainer,1))->Clone();
272 contaminationspectrum->SetName("contaminationspectrum");
273 TCanvas * ccontaminationspectrum = new TCanvas("contaminationspectrum","contaminationspectrum",1000,700);
274 ccontaminationspectrum->Divide(3,1);
275 ccontaminationspectrum->cd(1);
6555e2ad 276 TH2D * contaminationspectrum2dpteta = (TH2D *) contaminationspectrum->Project(1,0);
277 TH2D * contaminationspectrum2dptphi = (TH2D *) contaminationspectrum->Project(2,0);
278 TH2D * contaminationspectrum2detaphi = (TH2D *) contaminationspectrum->Project(1,2);
3a72645a 279 contaminationspectrum2dpteta->SetStats(0);
280 contaminationspectrum2dpteta->SetTitle("");
281 contaminationspectrum2dpteta->GetXaxis()->SetTitle("#eta");
282 contaminationspectrum2dpteta->GetYaxis()->SetTitle("p_{T} [GeV/c]");
283 contaminationspectrum2dptphi->SetStats(0);
284 contaminationspectrum2dptphi->SetTitle("");
285 contaminationspectrum2dptphi->GetXaxis()->SetTitle("#phi [rad]");
286 contaminationspectrum2dptphi->GetYaxis()->SetTitle("p_{T} [GeV/c]");
287 contaminationspectrum2detaphi->SetStats(0);
288 contaminationspectrum2detaphi->SetTitle("");
289 contaminationspectrum2detaphi->GetXaxis()->SetTitle("#eta");
290 contaminationspectrum2detaphi->GetYaxis()->SetTitle("#phi [rad]");
291 contaminationspectrum2dptphi->Draw("colz");
292 ccontaminationspectrum->cd(2);
293 contaminationspectrum2dpteta->Draw("colz");
294 ccontaminationspectrum->cd(3);
295 contaminationspectrum2detaphi->Draw("colz");
296
c2690925 297 TCanvas * ccorrelation = new TCanvas("ccorrelation","ccorrelation",1000,700);
298 ccorrelation->cd(1);
299 if(mccorrelationD) {
300 if(fBeamType==0){
301 TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions,0);
302 ptcorrelation->Draw("colz");
303 }
304 if(fBeamType==1){
305 TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions+1,1);
306 ptcorrelation->Draw("colz");
307 }
308 }
3a72645a 309 }
310
311
312 return kTRUE;
313
314
315}
c04c80e6 316
c2690925 317
c04c80e6 318//____________________________________________________________
3a72645a 319Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
c04c80e6 320 //
321 // Correct the spectrum for efficiency and unfolding
322 // with both method and compare
323 //
324
325 gStyle->SetPalette(1);
326 gStyle->SetOptStat(1111);
327 gStyle->SetPadBorderMode(0);
328 gStyle->SetCanvasColor(10);
329 gStyle->SetPadLeftMargin(0.13);
330 gStyle->SetPadRightMargin(0.13);
67fe7bd0 331
c2690925 332 printf("Steps are: stepdata %d, stepMC %d, steptrue %d, stepV0after %d, stepV0before %d\n",fStepData,fStepMC,fStepTrue,fStepAfterCutsV0,fStepBeforeCutsV0);
333
334 ///////////////////////////
335 // Check initialization
336 ///////////////////////////
337
338 if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
339 AliInfo("You have to init before");
340 return kFALSE;
341 }
342
343 if((fStepTrue == -1) && (fStepMC == -1) && (fStepData == -1)) {
344 AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
345 return kFALSE;
346 }
347
348 SetNumberOfIteration(10);
349 SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
350
351 AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
352 //////////////////////////////////
353 // Subtract hadron background
354 /////////////////////////////////
355 AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
356 if(subtractcontamination) {
357 dataspectrumaftersubstraction = SubtractBackground(kTRUE);
358 dataGridAfterFirstSteps = dataspectrumaftersubstraction;
359 }
360
361 ////////////////////////////////////////////////
362 // Correct for TPC efficiency from V0
363 ///////////////////////////////////////////////
364 AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
365 AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
366 if(dataContainerV0){
367 dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
368 dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
369 }
370
371 //////////////////////////////////////////////////////////////////////////////
372 // Correct for efficiency parametrized (if TPC efficiency is parametrized)
373 /////////////////////////////////////////////////////////////////////////////
374 AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
375 if(fEfficiencyFunction){
376 dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
377 dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
378 }
379
380 ///////////////
381 // Unfold
382 //////////////
383 TList *listunfolded = Unfold(dataGridAfterFirstSteps);
384 if(!listunfolded){
385 printf("Unfolded failed\n");
386 return kFALSE;
387 }
388 THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
389 THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
390 if(!correctedspectrum){
391 AliError("No corrected spectrum\n");
392 return kFALSE;
393 }
394 if(!residualspectrum){
395 AliError("No residul spectrum\n");
396 return kFALSE;
397 }
398
399 /////////////////////
400 // Simply correct
401 ////////////////////
402 AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
403
404
405 //////////
406 // Plot
407 //////////
408 if(fDebugLevel > 0.0) {
409
410
411
412
413 Int_t ptpr;
414 if(fBeamType==0) ptpr=0;
415 if(fBeamType==1) ptpr=1;
416
417 TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
418 ccorrected->Divide(2,1);
419 ccorrected->cd(1);
420 gPad->SetLogy();
421 TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
422 correctedspectrumD->SetTitle("");
423 correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
424 correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
425 correctedspectrumD->SetMarkerStyle(26);
426 correctedspectrumD->SetMarkerColor(kBlue);
427 correctedspectrumD->SetLineColor(kBlue);
428 correctedspectrumD->Draw("AP");
429 TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
430 alltogetherspectrumD->SetTitle("");
431 alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
432 alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
433 alltogetherspectrumD->SetMarkerStyle(25);
434 alltogetherspectrumD->SetMarkerColor(kBlack);
435 alltogetherspectrumD->SetLineColor(kBlack);
436 alltogetherspectrumD->Draw("P");
437 TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
438 legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
439 legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
440 legcorrected->Draw("same");
441 ccorrected->cd(2);
442 TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
443 TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
444 TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
445 ratiocorrected->SetName("ratiocorrected");
446 ratiocorrected->SetTitle("");
447 ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
448 ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
449 ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
450 ratiocorrected->SetStats(0);
451 ratiocorrected->Draw();
452
017dcb19 453 TH1D *unfoldingspectrac = new TH1D[fNCentralityBinAtTheEnd];
454 TGraphErrors *unfoldingspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
455 TH1D *correctedspectrac = new TH1D[fNCentralityBinAtTheEnd];
456 TGraphErrors *correctedspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
457
c2690925 458
c2690925 459
460 if(fBeamType==1) {
461
e156c3bb 462 TCanvas * ccorrectedallspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700);
463 ccorrectedallspectra->Divide(2,1);
c2690925 464 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
465 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
466
467 THnSparseF* sparsesu = (THnSparseF *) correctedspectrum;
468 TAxis *cenaxisa = sparsesu->GetAxis(0);
469 THnSparseF* sparsed = (THnSparseF *) alltogetherCorrection->GetGrid();
470 TAxis *cenaxisb = sparsed->GetAxis(0);
471 Int_t nbbin = cenaxisb->GetNbins();
472 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
473 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
474 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
475 TString titlee("corrected_centrality_bin_");
476 titlee += "[";
477 titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
478 titlee += "_";
479 titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
480 titlee += "[";
481 TString titleec("corrected_check_projection_bin_");
482 titleec += "[";
483 titleec += fLowBoundaryCentralityBinAtTheEnd[binc];
484 titleec += "_";
485 titleec += fHighBoundaryCentralityBinAtTheEnd[binc];
486 titleec += "[";
487 TString titleenameunotnormalized("Unfolded_Notnormalized_centrality_bin_");
488 titleenameunotnormalized += "[";
489 titleenameunotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
490 titleenameunotnormalized += "_";
491 titleenameunotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
492 titleenameunotnormalized += "[";
493 TString titleenameunormalized("Unfolded_normalized_centrality_bin_");
494 titleenameunormalized += "[";
495 titleenameunormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
496 titleenameunormalized += "_";
497 titleenameunormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
498 titleenameunormalized += "[";
499 TString titleenamednotnormalized("Dirrectcorrected_Notnormalized_centrality_bin_");
500 titleenamednotnormalized += "[";
501 titleenamednotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
502 titleenamednotnormalized += "_";
503 titleenamednotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
504 titleenamednotnormalized += "[";
505 TString titleenamednormalized("Dirrectedcorrected_normalized_centrality_bin_");
506 titleenamednormalized += "[";
507 titleenamednormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
508 titleenamednormalized += "_";
509 titleenamednormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
510 titleenamednormalized += "[";
511 Int_t nbEvents = 0;
512 for(Int_t k = fLowBoundaryCentralityBinAtTheEnd[binc]; k < fHighBoundaryCentralityBinAtTheEnd[binc]; k++) {
e156c3bb 513 printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
c2690925 514 nbEvents += fNEvents[k];
515 }
e156c3bb 516 Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
517 Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
518 printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
519 Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
520 Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
521 printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
c2690925 522 cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
523 cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
524 TCanvas * ccorrectedcheck = new TCanvas((const char*) titleec,(const char*) titleec,1000,700);
525 ccorrectedcheck->cd(1);
526 TH1D *aftersuc = (TH1D *) sparsesu->Projection(0);
527 TH1D *aftersdc = (TH1D *) sparsed->Projection(0);
528 aftersuc->Draw();
529 aftersdc->Draw("same");
530 TCanvas * ccorrectede = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
531 ccorrectede->Divide(2,1);
532 ccorrectede->cd(1);
533 gPad->SetLogy();
534 TH1D *aftersu = (TH1D *) sparsesu->Projection(1);
535 CorrectFromTheWidth(aftersu);
536 aftersu->SetName((const char*)titleenameunotnormalized);
6c70d827 537 unfoldingspectrac[binc] = *aftersu;
c2690925 538 ccorrectede->cd(1);
539 TGraphErrors* aftersun = NormalizeTH1N(aftersu,nbEvents);
540 aftersun->SetTitle("");
541 aftersun->GetYaxis()->SetTitleOffset(1.5);
542 aftersun->GetYaxis()->SetRangeUser(0.000000001,1.0);
543 aftersun->SetMarkerStyle(26);
544 aftersun->SetMarkerColor(kBlue);
545 aftersun->SetLineColor(kBlue);
546 aftersun->Draw("AP");
547 aftersun->SetName((const char*)titleenameunormalized);
6c70d827 548 unfoldingspectracn[binc] = *aftersun;
c2690925 549 ccorrectede->cd(1);
550 TH1D *aftersd = (TH1D *) sparsed->Projection(1);
551 CorrectFromTheWidth(aftersd);
552 aftersd->SetName((const char*)titleenamednotnormalized);
6c70d827 553 correctedspectrac[binc] = *aftersd;
c2690925 554 ccorrectede->cd(1);
555 TGraphErrors* aftersdn = NormalizeTH1N(aftersd,nbEvents);
556 aftersdn->SetTitle("");
557 aftersdn->GetYaxis()->SetTitleOffset(1.5);
558 aftersdn->GetYaxis()->SetRangeUser(0.000000001,1.0);
559 aftersdn->SetMarkerStyle(25);
560 aftersdn->SetMarkerColor(kBlack);
561 aftersdn->SetLineColor(kBlack);
562 aftersdn->Draw("P");
563 aftersdn->SetName((const char*)titleenamednormalized);
6c70d827 564 correctedspectracn[binc] = *aftersdn;
c2690925 565 TLegend *legcorrectedud = new TLegend(0.4,0.6,0.89,0.89);
566 legcorrectedud->AddEntry(aftersun,"Corrected","p");
567 legcorrectedud->AddEntry(aftersdn,"Alltogether","p");
568 legcorrectedud->Draw("same");
e156c3bb 569 ccorrectedallspectra->cd(1);
c2690925 570 gPad->SetLogy();
571 TH1D *aftersunn = (TH1D *) aftersun->Clone();
572 aftersunn->SetMarkerStyle(stylee[binc]);
573 aftersunn->SetMarkerColor(colorr[binc]);
574 if(binc==0) aftersunn->Draw("AP");
575 else aftersunn->Draw("P");
576 legtotal->AddEntry(aftersunn,(const char*) titlee,"p");
e156c3bb 577 ccorrectedallspectra->cd(2);
c2690925 578 gPad->SetLogy();
579 TH1D *aftersdnn = (TH1D *) aftersdn->Clone();
580 aftersdnn->SetMarkerStyle(stylee[binc]);
581 aftersdnn->SetMarkerColor(colorr[binc]);
582 if(binc==0) aftersdnn->Draw("AP");
583 else aftersdnn->Draw("P");
584 legtotalg->AddEntry(aftersdnn,(const char*) titlee,"p");
585 ccorrectede->cd(2);
586 TH1D* ratiocorrectedbinc = (TH1D*)aftersu->Clone();
587 TString titleee("ratiocorrected_bin_");
588 titleee += binc;
589 ratiocorrectedbinc->SetName((const char*) titleee);
590 ratiocorrectedbinc->SetTitle("");
591 ratiocorrectedbinc->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
592 ratiocorrectedbinc->GetXaxis()->SetTitle("p_{T} [GeV/c]");
593 ratiocorrectedbinc->Divide(aftersu,aftersd,1,1);
594 ratiocorrectedbinc->SetStats(0);
595 ratiocorrectedbinc->Draw();
596 }
597
e156c3bb 598 ccorrectedallspectra->cd(1);
c2690925 599 legtotal->Draw("same");
e156c3bb 600 ccorrectedallspectra->cd(2);
c2690925 601 legtotalg->Draw("same");
602
603 cenaxisa->SetRange(0,nbbin);
604 cenaxisb->SetRange(0,nbbin);
605
606 }
607
608 // Dump to file if needed
609 if(fDumpToFile) {
610 TFile *out = new TFile("finalSpectrum.root","recreate");
611 correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
612 correctedspectrumD->Write();
613 alltogetherspectrumD->SetName("AlltogetherSpectrum");
614 alltogetherspectrumD->Write();
615 ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
616 ratiocorrected->Write();
617 correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
618 correctedspectrum->Write();
619 alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
620 alltogetherCorrection->Write();
621 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
6c70d827 622 unfoldingspectrac[binc].Write();
623 unfoldingspectracn[binc].Write();
624 correctedspectrac[binc].Write();
625 correctedspectracn[binc].Write();
c2690925 626 }
627 out->Close(); delete out;
628 }
629
017dcb19 630 delete [] unfoldingspectrac;
631 delete [] unfoldingspectracn;
632 delete [] correctedspectrac;
633 delete [] correctedspectracn;
c2690925 634
017dcb19 635 }
c2690925 636
637 return kTRUE;
638}
639
640//____________________________________________________________
641Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
642 //
643 // Correct the spectrum for efficiency and unfolding for beauty analysis
644 // with both method and compare
645 //
646
647 gStyle->SetPalette(1);
648 gStyle->SetOptStat(1111);
649 gStyle->SetPadBorderMode(0);
650 gStyle->SetCanvasColor(10);
651 gStyle->SetPadLeftMargin(0.13);
652 gStyle->SetPadRightMargin(0.13);
653
3a72645a 654 ///////////////////////////
655 // Check initialization
656 ///////////////////////////
c04c80e6 657
3a72645a 658 if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
659 AliInfo("You have to init before");
660 return kFALSE;
661 }
662
663 if((fStepTrue == 0) && (fStepMC == 0) && (fStepData == 0)) {
664 AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
665 return kFALSE;
666 }
667
c2690925 668 SetNumberOfIteration(10);
3a72645a 669 SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
670
671 AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
672 //////////////////////////////////
673 // Subtract hadron background
674 /////////////////////////////////
67fe7bd0 675 AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
3a72645a 676 if(subtractcontamination) {
677 dataspectrumaftersubstraction = SubtractBackground(kTRUE);
678 dataGridAfterFirstSteps = dataspectrumaftersubstraction;
679 }
680
c2690925 681 /////////////////////////////////////////////////////////////////////////////////////////
682 // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
683 /////////////////////////////////////////////////////////////////////////////////////////
684
685
3a72645a 686 ////////////////////////////////////////////////
687 // Correct for TPC efficiency from V0
688 ///////////////////////////////////////////////
689 AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
690 AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
691 if(dataContainerV0){
692 dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
693 dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
67fe7bd0 694 }
c04c80e6 695
3a72645a 696 ///////////////
c04c80e6 697 // Unfold
3a72645a 698 //////////////
699 TList *listunfolded = Unfold(dataGridAfterFirstSteps);
c04c80e6 700 if(!listunfolded){
701 printf("Unfolded failed\n");
3a72645a 702 return kFALSE;
c04c80e6 703 }
704 THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
705 THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
706 if(!correctedspectrum){
707 AliError("No corrected spectrum\n");
3a72645a 708 return kFALSE;
c04c80e6 709 }
67fe7bd0 710 if(!residualspectrum){
c04c80e6 711 AliError("No residul spectrum\n");
3a72645a 712 return kFALSE;
c04c80e6 713 }
67fe7bd0 714
3a72645a 715 /////////////////////
c04c80e6 716 // Simply correct
3a72645a 717 ////////////////////
718 AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
67fe7bd0 719
3a72645a 720
67fe7bd0 721 //////////
c04c80e6 722 // Plot
723 //////////
3a72645a 724 if(fDebugLevel > 0.0) {
725
726 TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
727 ccorrected->Divide(2,1);
728 ccorrected->cd(1);
729 gPad->SetLogy();
730 TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
731 correctedspectrumD->SetTitle("");
732 correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
733 correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
734 correctedspectrumD->SetMarkerStyle(26);
735 correctedspectrumD->SetMarkerColor(kBlue);
736 correctedspectrumD->SetLineColor(kBlue);
737 correctedspectrumD->Draw("AP");
738 TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
739 alltogetherspectrumD->SetTitle("");
740 alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
741 alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
742 alltogetherspectrumD->SetMarkerStyle(25);
743 alltogetherspectrumD->SetMarkerColor(kBlack);
744 alltogetherspectrumD->SetLineColor(kBlack);
745 alltogetherspectrumD->Draw("P");
746 TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
747 legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
748 legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
749 legcorrected->Draw("same");
750 ccorrected->cd(2);
751 TH1D *correctedTH1D = correctedspectrum->Projection(0);
6555e2ad 752 TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(0);
3a72645a 753 TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
754 ratiocorrected->SetName("ratiocorrected");
755 ratiocorrected->SetTitle("");
756 ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
757 ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
758 ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
759 ratiocorrected->SetStats(0);
760 ratiocorrected->Draw();
761
762
763 // Dump to file if needed
764
765 if(fDumpToFile) {
766
767 TFile *out = new TFile("finalSpectrum.root","recreate");
768 out->cd();
769 //
770 correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
771 correctedspectrumD->Write();
772 alltogetherspectrumD->SetName("AlltogetherSpectrum");
773 alltogetherspectrumD->Write();
774 ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
775 ratiocorrected->Write();
776 //
777 correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
778 correctedspectrum->Write();
779 alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
780 alltogetherCorrection->Write();
781 //
782 out->Close(); delete out;
783 }
784
785
786 }
787
788
789
790
791 return kTRUE;
792}
c2690925 793
3a72645a 794//____________________________________________________________
795AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
796 //
797 // Apply background subtraction
798 //
799
800 // Raw spectrum
801 AliCFContainer *dataContainer = GetContainer(kDataContainer);
802 if(!dataContainer){
803 AliError("Data Container not available");
804 return NULL;
805 }
c2690925 806 printf("Step data: %d\n",fStepData);
3a72645a 807 AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
808
809 AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
810 dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
811
812 // Background Estimate
813 AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
814 if(!backgroundContainer){
815 AliError("MC background container not found");
816 return NULL;
817 }
818
c2690925 819 Int_t stepbackground = 1; // 2 for !fInclusiveSpectrum analysis(old method)
3a72645a 820 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
821
c2690925 822 if(!fInclusiveSpectrum){
823 //Background subtraction for IP analysis
824 if(fIPanaHadronBgSubtract){
825 // Hadron background
826 printf("Hadron background for IP analysis subtracted!\n");
827 Int_t* bins=new Int_t[1];
828 TH1D* htemp = (TH1D *) fHadronEffbyIPcut->Projection(0);
829 bins[0]=htemp->GetNbinsX();
830 AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",1,bins);
831 hbgContainer->SetGrid(fHadronEffbyIPcut);
832 backgroundGrid->Multiply(hbgContainer,1);
833 spectrumSubtracted->Add(backgroundGrid,-1.0);
834 }
835 if(fIPanaCharmBgSubtract){
836 // Charm background
837 printf("Charm background for IP analysis subtracted!\n");
838 AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
839 spectrumSubtracted->Add(charmbgContainer,-1.0);
840 }
841 if(fIPanaConversionBgSubtract){
842 // Conversion background
843 AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground();
844 spectrumSubtracted->Add(conversionbgContainer,-1.0);
845 printf("Conversion background subtraction is preliminary!\n");
846 }
847 if(fIPanaNonHFEBgSubtract){
848 // NonHFE background
849 AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground();
850 spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
851 printf("Non HFE background subtraction is preliminary!\n");
852 }
853 }
854 else{
855 // Subtract
856 spectrumSubtracted->Add(backgroundGrid,-1.0);
857 }
858
3a72645a 859 if(setBackground){
860 if(fBackground) delete fBackground;
861 fBackground = backgroundGrid;
862 } else delete backgroundGrid;
863
864
865 if(fDebugLevel > 0) {
c2690925 866
867 Int_t ptpr;
868 if(fBeamType==0) ptpr=0;
869 if(fBeamType==1) ptpr=1;
67fe7bd0 870
871 TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
3a72645a 872 cbackgroundsubtraction->Divide(3,1);
67fe7bd0 873 cbackgroundsubtraction->cd(1);
3a72645a 874 gPad->SetLogy();
c2690925 875 TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptpr);
876 TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
67fe7bd0 877 CorrectFromTheWidth(measuredTH1Daftersubstraction);
878 CorrectFromTheWidth(measuredTH1Dbeforesubstraction);
879 measuredTH1Daftersubstraction->SetStats(0);
880 measuredTH1Daftersubstraction->SetTitle("");
881 measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
882 measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
883 measuredTH1Daftersubstraction->SetMarkerStyle(25);
884 measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
885 measuredTH1Daftersubstraction->SetLineColor(kBlack);
886 measuredTH1Dbeforesubstraction->SetStats(0);
887 measuredTH1Dbeforesubstraction->SetTitle("");
888 measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
889 measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
890 measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
891 measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
892 measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
893 measuredTH1Daftersubstraction->Draw();
894 measuredTH1Dbeforesubstraction->Draw("same");
895 TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
3a72645a 896 legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
897 legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
67fe7bd0 898 legsubstraction->Draw("same");
899 cbackgroundsubtraction->cd(2);
3a72645a 900 gPad->SetLogy();
67fe7bd0 901 TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
902 ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
903 ratiomeasuredcontamination->SetTitle("");
904 ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
905 ratiomeasuredcontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
c2690925 906 ratiomeasuredcontamination->Sumw2();
67fe7bd0 907 ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
908 ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
909 ratiomeasuredcontamination->SetStats(0);
910 ratiomeasuredcontamination->SetMarkerStyle(26);
911 ratiomeasuredcontamination->SetMarkerColor(kBlack);
912 ratiomeasuredcontamination->SetLineColor(kBlack);
c2690925 913 for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
914 ratiomeasuredcontamination->SetBinError(k+1,0.0);
915 }
916 ratiomeasuredcontamination->Draw("P");
3a72645a 917 cbackgroundsubtraction->cd(3);
c2690925 918 TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptpr);
3a72645a 919 CorrectFromTheWidth(measuredTH1background);
920 measuredTH1background->SetStats(0);
921 measuredTH1background->SetTitle("");
922 measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
923 measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]");
924 measuredTH1background->SetMarkerStyle(26);
925 measuredTH1background->SetMarkerColor(kBlack);
926 measuredTH1background->SetLineColor(kBlack);
927 measuredTH1background->Draw();
c2690925 928
929 if(fBeamType==1) {
930
931 TCanvas * cbackgrounde = new TCanvas("BackgroundSubtraction_allspectra","BackgroundSubtraction_allspectra",1000,700);
932 cbackgrounde->Divide(2,1);
933 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
934 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
935
936 THnSparseF* sparsesubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
937 TAxis *cenaxisa = sparsesubtracted->GetAxis(0);
938 THnSparseF* sparsebefore = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
939 TAxis *cenaxisb = sparsebefore->GetAxis(0);
940 Int_t nbbin = cenaxisb->GetNbins();
941 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
942 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
943 for(Int_t binc = 0; binc < nbbin; binc++){
944 TString titlee("BackgroundSubtraction_centrality_bin_");
945 titlee += binc;
946 TCanvas * cbackground = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
947 cbackground->Divide(2,1);
948 cbackground->cd(1);
949 gPad->SetLogy();
950 cenaxisa->SetRange(binc+1,binc+1);
951 cenaxisb->SetRange(binc+1,binc+1);
952 TH1D *aftersubstraction = (TH1D *) sparsesubtracted->Projection(1);
953 TH1D *beforesubstraction = (TH1D *) sparsebefore->Projection(1);
954 CorrectFromTheWidth(aftersubstraction);
955 CorrectFromTheWidth(beforesubstraction);
956 aftersubstraction->SetStats(0);
957 aftersubstraction->SetTitle((const char*)titlee);
958 aftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
959 aftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
960 aftersubstraction->SetMarkerStyle(25);
961 aftersubstraction->SetMarkerColor(kBlack);
962 aftersubstraction->SetLineColor(kBlack);
963 beforesubstraction->SetStats(0);
964 beforesubstraction->SetTitle((const char*)titlee);
965 beforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
966 beforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
967 beforesubstraction->SetMarkerStyle(24);
968 beforesubstraction->SetMarkerColor(kBlue);
969 beforesubstraction->SetLineColor(kBlue);
970 aftersubstraction->Draw();
971 beforesubstraction->Draw("same");
972 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
973 lega->AddEntry(beforesubstraction,"With hadron contamination","p");
974 lega->AddEntry(aftersubstraction,"Without hadron contamination ","p");
975 lega->Draw("same");
976 cbackgrounde->cd(1);
977 gPad->SetLogy();
978 TH1D *aftersubtractionn = (TH1D *) aftersubstraction->Clone();
979 aftersubtractionn->SetMarkerStyle(stylee[binc]);
980 aftersubtractionn->SetMarkerColor(colorr[binc]);
981 if(binc==0) aftersubtractionn->Draw();
982 else aftersubtractionn->Draw("same");
983 legtotal->AddEntry(aftersubtractionn,(const char*) titlee,"p");
984 cbackgrounde->cd(2);
985 gPad->SetLogy();
986 TH1D *aftersubtractionng = (TH1D *) aftersubstraction->Clone();
987 aftersubtractionng->SetMarkerStyle(stylee[binc]);
988 aftersubtractionng->SetMarkerColor(colorr[binc]);
989 if(fNEvents[binc] > 0.0) aftersubtractionng->Scale(1/(Double_t)fNEvents[binc]);
990 if(binc==0) aftersubtractionng->Draw();
991 else aftersubtractionng->Draw("same");
992 legtotalg->AddEntry(aftersubtractionng,(const char*) titlee,"p");
993 cbackground->cd(2);
994 TH1D* ratiocontamination = (TH1D*)beforesubstraction->Clone();
995 ratiocontamination->SetName("ratiocontamination");
996 ratiocontamination->SetTitle((const char*)titlee);
997 ratiocontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
998 ratiocontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
999 ratiocontamination->Add(aftersubstraction,-1.0);
1000 ratiocontamination->Divide(beforesubstraction);
1001 Int_t totalbin = ratiocontamination->GetXaxis()->GetNbins();
1002 for(Int_t nbinpt = 0; nbinpt < totalbin; nbinpt++) {
1003 ratiocontamination->SetBinError(nbinpt+1,0.0);
1004 }
1005 ratiocontamination->SetStats(0);
1006 ratiocontamination->SetMarkerStyle(26);
1007 ratiocontamination->SetMarkerColor(kBlack);
1008 ratiocontamination->SetLineColor(kBlack);
1009 ratiocontamination->Draw("P");
1010
1011 }
1012
1013 cbackgrounde->cd(1);
1014 legtotal->Draw("same");
1015 cbackgrounde->cd(2);
1016 legtotalg->Draw("same");
1017
1018 cenaxisa->SetRange(0,nbbin);
1019 cenaxisb->SetRange(0,nbbin);
1020
1021 }
1022
3a72645a 1023
67fe7bd0 1024 }
1025
c04c80e6 1026
3a72645a 1027 return spectrumSubtracted;
c04c80e6 1028}
c2690925 1029
1030//____________________________________________________________
1031AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
1032 //
1033 // calculate charm background
1034 //
1035
1036 Double_t evtnorm=0;
1037 if(fNMCEvents) evtnorm= double(fNEvents[0])/double(fNMCEvents);
1038
1039 AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
1040 if(!mcContainer){
1041 AliError("MC Container not available");
1042 return NULL;
1043 }
1044
1045 if(!fCorrelation){
1046 AliError("No Correlation map available");
1047 return NULL;
1048 }
1049
1050 Int_t nDim = 1;
1051 AliCFDataGrid *charmBackgroundGrid= 0x0;
1052 charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC+2);
1053 TH1D *tmphisto = (TH1D *) charmBackgroundGrid->Project(0);
1054 Int_t* bins=new Int_t[1];
1055 bins[0]=tmphisto->GetNbinsX();
1056 AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",1,bins);
1057 weightedCharmContainer->SetGrid(GetWeights());
1058 charmBackgroundGrid->Multiply(weightedCharmContainer,evtnorm);
1059
1060 // Efficiency (set efficiency to 1 for only folding)
1061 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1062 efficiencyD->CalculateEfficiency(0,0);
1063
1064 // Folding
1065 AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
1066 folding.SetMaxNumberOfIterations(1);
1067 folding.Unfold();
1068
1069 // Results
1070 THnSparse* result1= folding.GetEstMeasured(); // folded spectra
1071 THnSparse* result=(THnSparse*)result1->Clone();
1072 charmBackgroundGrid->SetGrid(result);
1073
1074 return charmBackgroundGrid;
1075}
1076
1077//____________________________________________________________
1078AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
1079 //
1080 // calculate conversion background
1081 //
1082
1083 Double_t evtnorm[1];
1084 if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
1085 printf("check event!!! %lf \n",evtnorm[0]);
1086
1087 // Background Estimate
1088 AliCFContainer *backgroundContainer = GetContainer(kNonHFEBackgroundData);
1089 if(!backgroundContainer){
1090 AliError("MC background container not found");
1091 return NULL;
1092 }
1093
1094 Int_t stepbackground = 2;
1095 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground);
1096 backgroundGrid->Scale(evtnorm);
1097
1098 new TCanvas;
1099 TH1D *histodraw111 = (TH1D *) backgroundGrid->Project(0);
1100 CorrectFromTheWidth(histodraw111);
1101 histodraw111->SetMarkerStyle(20);
1102 histodraw111->SetLineColor(4);
1103 histodraw111->SetMarkerColor(4);
1104 histodraw111->Draw("p");
1105
1106
1107 return backgroundGrid;
1108}
1109
1110
1111//____________________________________________________________
1112AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
1113 //
1114 // calculate non HFE background
1115 //
1116
1117 if(fNonHFEbgMethod2){
1118 Double_t evtnorm=0;
1119 if(fNMCEvents) evtnorm= double(fNEvents[0])/double(fNMCEvents);
1120
1121 AliCFContainer *mcContainer = GetContainer(kMCContainerNonHFEMC);
1122 if(!mcContainer){
1123 AliError("MC Container not available");
1124 return NULL;
1125 }
1126
1127 if(!fCorrelation){
1128 AliError("No Correlation map available");
1129 return NULL;
1130 }
1131
1132 Int_t nDim = 1;
1133 AliCFDataGrid *nonHFEBackgroundGrid= 0x0;
1134 nonHFEBackgroundGrid = new AliCFDataGrid("nonHFEBackgroundGrid","nonHFEBackgroundGrid",*mcContainer, fStepMC+2);
1135 TH1D *tmphisto = (TH1D *) nonHFEBackgroundGrid->Project(0);
1136 Int_t* bins=new Int_t[1];
1137 bins[0]=tmphisto->GetNbinsX();
1138 AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
1139 weightedNonHFEContainer->SetGrid(GetWeights()); // we have to use different weithing function
1140 nonHFEBackgroundGrid->Multiply(weightedNonHFEContainer,evtnorm);
1141
1142 // Efficiency (set efficiency to 1 for only folding)
1143 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1144 efficiencyD->CalculateEfficiency(0,0);
1145
1146 // Folding
1147 AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),nonHFEBackgroundGrid->GetGrid(),nonHFEBackgroundGrid->GetGrid());
1148 folding.SetMaxNumberOfIterations(1);
1149 folding.Unfold();
1150
1151 // Results
1152 THnSparse* result1= folding.GetEstMeasured(); // folded spectra
1153 THnSparse* result=(THnSparse*)result1->Clone();
1154 nonHFEBackgroundGrid->SetGrid(result);
1155
1156 return nonHFEBackgroundGrid;
1157 }
1158 else{
1159 Double_t evtnorm[1];
1160 if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
1161
1162 // Background Estimate
1163 AliCFContainer *backgroundContainer = GetContainer(kNonHFEBackgroundData);
1164 if(!backgroundContainer){
1165 AliError("MC background container not found");
1166 return NULL;
1167 }
1168
1169 Int_t stepbackground = 3;
1170 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
1171 backgroundGrid->Scale(evtnorm);
1172
1173 new TCanvas;
1174 TH1D *histodraw222 = (TH1D *) backgroundGrid->Project(0);
1175 CorrectFromTheWidth(histodraw222);
1176 histodraw222->SetMarkerStyle(20);
1177 histodraw222->SetLineColor(4);
1178 histodraw222->SetMarkerColor(4);
1179 histodraw222->Draw("p");
1180
1181 return backgroundGrid;
1182 }
1183
1184}
1185
1186//____________________________________________________________
1187AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
1188
1189 //
1190 // Apply TPC pid efficiency correction from parametrisation
1191 //
1192
1193 // Data in the right format
1194 AliCFDataGrid *dataGrid = 0x0;
1195 if(bgsubpectrum) {
1196 dataGrid = bgsubpectrum;
1197 }
1198 else {
1199
1200 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1201 if(!dataContainer){
1202 AliError("Data Container not available");
1203 return NULL;
1204 }
1205
1206 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1207 }
1208
1209 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1210 result->SetName("ParametrizedEfficiencyBefore");
1211 THnSparse *h = result->GetGrid();
1212 Int_t nbdimensions = h->GetNdimensions();
1213 //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
1214
1215 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1216 if(!dataContainer){
1217 AliError("Data Container not available");
1218 return NULL;
1219 }
1220 AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
1221 dataContainerbis->Add(dataContainerbis,-1.0);
1222
1223
1224 Int_t* coord = new Int_t[nbdimensions];
1225 memset(coord, 0, sizeof(Int_t) * nbdimensions);
1226 Double_t* points = new Double_t[nbdimensions];
1227
1228
1229 ULong64_t nEntries = h->GetNbins();
1230 for (ULong64_t i = 0; i < nEntries; ++i) {
1231
1232 Double_t value = h->GetBinContent(i, coord);
1233 //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
1234 //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
1235
1236 // Get the bin co-ordinates given an coord
1237 for (Int_t j = 0; j < nbdimensions; ++j)
1238 points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1239
1240 if (!fEfficiencyFunction->IsInside(points))
1241 continue;
1242 TF1::RejectPoint(kFALSE);
1243
1244 // Evaulate function at points
1245 Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1246 //printf("Value efficiency is %f\n",valueEfficiency);
1247
1248 if(valueEfficiency > 0.0) {
1249 h->SetBinContent(coord,value/valueEfficiency);
1250 dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1251 }
1252 Double_t error = h->GetBinError(i);
1253 h->SetBinError(coord,error/valueEfficiency);
1254 dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1255
1256
1257 }
1258
6c70d827 1259 delete[] coord;
1260 delete[] points;
1261
c2690925 1262 AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1263
1264 if(fDebugLevel > 0) {
1265
1266 TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700);
1267 cEfficiencyParametrized->Divide(2,1);
1268 cEfficiencyParametrized->cd(1);
1269 TH1D *afterE = (TH1D *) resultt->Project(0);
1270 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1271 CorrectFromTheWidth(afterE);
1272 CorrectFromTheWidth(beforeE);
1273 afterE->SetStats(0);
1274 afterE->SetTitle("");
1275 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1276 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1277 afterE->SetMarkerStyle(25);
1278 afterE->SetMarkerColor(kBlack);
1279 afterE->SetLineColor(kBlack);
1280 beforeE->SetStats(0);
1281 beforeE->SetTitle("");
1282 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1283 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1284 beforeE->SetMarkerStyle(24);
1285 beforeE->SetMarkerColor(kBlue);
1286 beforeE->SetLineColor(kBlue);
1287 gPad->SetLogy();
1288 afterE->Draw();
1289 beforeE->Draw("same");
1290 TLegend *legefficiencyparametrized = new TLegend(0.4,0.6,0.89,0.89);
1291 legefficiencyparametrized->AddEntry(beforeE,"Before Efficiency correction","p");
1292 legefficiencyparametrized->AddEntry(afterE,"After Efficiency correction","p");
1293 legefficiencyparametrized->Draw("same");
1294 cEfficiencyParametrized->cd(2);
1295 fEfficiencyFunction->Draw();
1296 //cEfficiencyParametrized->cd(3);
1297 //TH1D *ratioefficiency = (TH1D *) beforeE->Clone();
1298 //ratioefficiency->Divide(afterE);
1299 //ratioefficiency->Draw();
1300
1301
1302 }
1303
1304
1305 return resultt;
1306
1307}
c04c80e6 1308//____________________________________________________________
3a72645a 1309AliCFDataGrid *AliHFEspectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
1310
c04c80e6 1311 //
3a72645a 1312 // Apply TPC pid efficiency correction from V0
c04c80e6 1313 //
1314
3a72645a 1315 AliCFContainer *v0Container = GetContainer(kDataContainerV0);
1316 if(!v0Container){
1317 AliError("V0 Container not available");
c04c80e6 1318 return NULL;
1319 }
3a72645a 1320
1321 // Efficiency
1322 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
1323 efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
1324
1325 // Data in the right format
1326 AliCFDataGrid *dataGrid = 0x0;
1327 if(bgsubpectrum) {
1328 dataGrid = bgsubpectrum;
c04c80e6 1329 }
3a72645a 1330 else {
1331
1332 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1333 if(!dataContainer){
1334 AliError("Data Container not available");
1335 return NULL;
1336 }
c04c80e6 1337
3a72645a 1338 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1339 }
c04c80e6 1340
3a72645a 1341 // Correct
1342 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1343 result->ApplyEffCorrection(*efficiencyD);
c04c80e6 1344
3a72645a 1345 if(fDebugLevel > 0) {
c2690925 1346
1347 Int_t ptpr;
1348 if(fBeamType==0) ptpr=0;
1349 if(fBeamType==1) ptpr=1;
3a72645a 1350
1351 TCanvas * cV0Efficiency = new TCanvas("V0Efficiency","V0Efficiency",1000,700);
1352 cV0Efficiency->Divide(2,1);
1353 cV0Efficiency->cd(1);
c2690925 1354 TH1D *afterE = (TH1D *) result->Project(ptpr);
1355 TH1D *beforeE = (TH1D *) dataGrid->Project(ptpr);
3a72645a 1356 afterE->SetStats(0);
1357 afterE->SetTitle("");
1358 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1359 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1360 afterE->SetMarkerStyle(25);
1361 afterE->SetMarkerColor(kBlack);
1362 afterE->SetLineColor(kBlack);
1363 beforeE->SetStats(0);
1364 beforeE->SetTitle("");
1365 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1366 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1367 beforeE->SetMarkerStyle(24);
1368 beforeE->SetMarkerColor(kBlue);
1369 beforeE->SetLineColor(kBlue);
1370 afterE->Draw();
1371 beforeE->Draw("same");
1372 TLegend *legV0efficiency = new TLegend(0.4,0.6,0.89,0.89);
1373 legV0efficiency->AddEntry(beforeE,"Before Efficiency correction","p");
1374 legV0efficiency->AddEntry(afterE,"After Efficiency correction","p");
1375 legV0efficiency->Draw("same");
1376 cV0Efficiency->cd(2);
c2690925 1377 TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(ptpr);
3a72645a 1378 efficiencyDproj->SetTitle("");
1379 efficiencyDproj->SetStats(0);
1380 efficiencyDproj->SetMarkerStyle(25);
1381 efficiencyDproj->Draw();
c04c80e6 1382
c2690925 1383 if(fBeamType==1) {
1384
1385 TCanvas * cV0Efficiencye = new TCanvas("V0Efficiency_allspectra","V0Efficiency_allspectra",1000,700);
1386 cV0Efficiencye->Divide(2,1);
1387 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1388 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1389
1390 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1391 TAxis *cenaxisa = sparseafter->GetAxis(0);
1392 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1393 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1394 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1395 TAxis *cenaxisc = efficiencya->GetAxis(0);
1396 Int_t nbbin = cenaxisb->GetNbins();
1397 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1398 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1399 for(Int_t binc = 0; binc < nbbin; binc++){
1400 TString titlee("V0Efficiency_centrality_bin_");
1401 titlee += binc;
1402 TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1403 ccV0Efficiency->Divide(2,1);
1404 ccV0Efficiency->cd(1);
1405 gPad->SetLogy();
1406 cenaxisa->SetRange(binc+1,binc+1);
1407 cenaxisb->SetRange(binc+1,binc+1);
1408 cenaxisc->SetRange(binc+1,binc+1);
1409 TH1D *aftere = (TH1D *) sparseafter->Projection(1);
1410 TH1D *beforee = (TH1D *) sparsebefore->Projection(1);
1411 CorrectFromTheWidth(aftere);
1412 CorrectFromTheWidth(beforee);
1413 aftere->SetStats(0);
1414 aftere->SetTitle((const char*)titlee);
1415 aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1416 aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1417 aftere->SetMarkerStyle(25);
1418 aftere->SetMarkerColor(kBlack);
1419 aftere->SetLineColor(kBlack);
1420 beforee->SetStats(0);
1421 beforee->SetTitle((const char*)titlee);
1422 beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1423 beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1424 beforee->SetMarkerStyle(24);
1425 beforee->SetMarkerColor(kBlue);
1426 beforee->SetLineColor(kBlue);
1427 aftere->Draw();
1428 beforee->Draw("same");
1429 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1430 lega->AddEntry(beforee,"Before correction","p");
1431 lega->AddEntry(aftere,"After correction","p");
1432 lega->Draw("same");
1433 cV0Efficiencye->cd(1);
1434 gPad->SetLogy();
1435 TH1D *afteree = (TH1D *) aftere->Clone();
1436 afteree->SetMarkerStyle(stylee[binc]);
1437 afteree->SetMarkerColor(colorr[binc]);
1438 if(binc==0) afteree->Draw();
1439 else afteree->Draw("same");
1440 legtotal->AddEntry(afteree,(const char*) titlee,"p");
1441 cV0Efficiencye->cd(2);
1442 gPad->SetLogy();
1443 TH1D *aftereeu = (TH1D *) aftere->Clone();
1444 aftereeu->SetMarkerStyle(stylee[binc]);
1445 aftereeu->SetMarkerColor(colorr[binc]);
1446 if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]);
1447 if(binc==0) aftereeu->Draw();
1448 else aftereeu->Draw("same");
1449 legtotalg->AddEntry(aftereeu,(const char*) titlee,"p");
1450 ccV0Efficiency->cd(2);
1451 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1);
1452 efficiencyDDproj->SetTitle("");
1453 efficiencyDDproj->SetStats(0);
1454 efficiencyDDproj->SetMarkerStyle(25);
1455 efficiencyDDproj->Draw();
1456
1457 }
1458
1459 cV0Efficiencye->cd(1);
1460 legtotal->Draw("same");
1461 cV0Efficiencye->cd(2);
1462 legtotalg->Draw("same");
1463
1464 cenaxisa->SetRange(0,nbbin);
1465 cenaxisb->SetRange(0,nbbin);
1466 cenaxisc->SetRange(0,nbbin);
1467
1468 }
3a72645a 1469
1470 }
1471
1472
1473 return result;
1474
1475}
c04c80e6 1476//____________________________________________________________
3a72645a 1477TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
c04c80e6 1478
1479 //
1480 // Unfold and eventually correct for efficiency the bgsubspectrum
1481 //
1482
3a72645a 1483 AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
c04c80e6 1484 if(!mcContainer){
1485 AliError("MC Container not available");
1486 return NULL;
1487 }
1488
1489 if(!fCorrelation){
1490 AliError("No Correlation map available");
1491 return NULL;
1492 }
1493
3a72645a 1494 // Data
c04c80e6 1495 AliCFDataGrid *dataGrid = 0x0;
1496 if(bgsubpectrum) {
c04c80e6 1497 dataGrid = bgsubpectrum;
c04c80e6 1498 }
1499 else {
1500
1501 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1502 if(!dataContainer){
1503 AliError("Data Container not available");
1504 return NULL;
1505 }
1506
3a72645a 1507 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
c04c80e6 1508 }
1509
c04c80e6 1510 // Guessed
3a72645a 1511 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
c04c80e6 1512 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
1513
1514 // Efficiency
3a72645a 1515 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
c04c80e6 1516 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1517
1518 // Unfold
1519
3a72645a 1520 AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
c2690925 1521 //unfolding.SetUseCorrelatedErrors();
1522 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
c04c80e6 1523 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
e156c3bb 1524 if(fSetSmoothing) unfolding.UseSmoothing();
c04c80e6 1525 unfolding.Unfold();
1526
1527 // Results
1528 THnSparse* result = unfolding.GetUnfolded();
1529 THnSparse* residual = unfolding.GetEstMeasured();
1530 TList *listofresults = new TList;
1531 listofresults->SetOwner();
1532 listofresults->AddAt((THnSparse*)result->Clone(),0);
1533 listofresults->AddAt((THnSparse*)residual->Clone(),1);
1534
3a72645a 1535 if(fDebugLevel > 0) {
c2690925 1536
1537 Int_t ptpr;
1538 if(fBeamType==0) ptpr=0;
1539 if(fBeamType==1) ptpr=1;
3a72645a 1540
1541 TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
1542 cresidual->Divide(2,1);
1543 cresidual->cd(1);
1544 gPad->SetLogy();
1545 TGraphErrors* residualspectrumD = Normalize(residual);
1546 if(!residualspectrumD) {
1547 AliError("Number of Events not set for the normalization");
1548 return kFALSE;
1549 }
1550 residualspectrumD->SetTitle("");
1551 residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
1552 residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
1553 residualspectrumD->SetMarkerStyle(26);
1554 residualspectrumD->SetMarkerColor(kBlue);
1555 residualspectrumD->SetLineColor(kBlue);
1556 residualspectrumD->Draw("AP");
1557 AliCFDataGrid *dataGridBis = (AliCFDataGrid *) (dataGrid->Clone());
1558 dataGridBis->SetName("dataGridBis");
1559 TGraphErrors* measuredspectrumD = Normalize(dataGridBis);
1560 measuredspectrumD->SetTitle("");
1561 measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
1562 measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
1563 measuredspectrumD->SetMarkerStyle(25);
1564 measuredspectrumD->SetMarkerColor(kBlack);
1565 measuredspectrumD->SetLineColor(kBlack);
1566 measuredspectrumD->Draw("P");
1567 TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
1568 legres->AddEntry(residualspectrumD,"Residual","p");
1569 legres->AddEntry(measuredspectrumD,"Measured","p");
1570 legres->Draw("same");
1571 cresidual->cd(2);
c2690925 1572 TH1D *residualTH1D = residual->Projection(ptpr);
1573 TH1D *measuredTH1D = (TH1D *) dataGridBis->Project(ptpr);
3a72645a 1574 TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
1575 ratioresidual->SetName("ratioresidual");
1576 ratioresidual->SetTitle("");
1577 ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
1578 ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1579 ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
1580 ratioresidual->SetStats(0);
1581 ratioresidual->Draw();
1582
1583 }
1584
c04c80e6 1585 return listofresults;
1586
1587}
1588
1589//____________________________________________________________
3a72645a 1590AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
c04c80e6 1591
1592 //
1593 // Apply unfolding and efficiency correction together to bgsubspectrum
1594 //
1595
3a72645a 1596 AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
c04c80e6 1597 if(!mcContainer){
1598 AliError("MC Container not available");
1599 return NULL;
1600 }
1601
c04c80e6 1602 // Efficiency
3a72645a 1603 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
c04c80e6 1604 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1605
1606 // Data in the right format
1607 AliCFDataGrid *dataGrid = 0x0;
1608 if(bgsubpectrum) {
c04c80e6 1609 dataGrid = bgsubpectrum;
c04c80e6 1610 }
1611 else {
3a72645a 1612
c04c80e6 1613 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1614 if(!dataContainer){
1615 AliError("Data Container not available");
1616 return NULL;
1617 }
1618
3a72645a 1619 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
c04c80e6 1620 }
1621
1622 // Correct
1623 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1624 result->ApplyEffCorrection(*efficiencyD);
c04c80e6 1625
3a72645a 1626 if(fDebugLevel > 0) {
c2690925 1627
1628 Int_t ptpr;
1629 if(fBeamType==0) ptpr=0;
1630 if(fBeamType==1) ptpr=1;
3a72645a 1631
bf892a6a 1632 printf("Step MC: %d\n",fStepMC);
1633 printf("Step tracking: %d\n",AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
1634 printf("Step MC true: %d\n",fStepTrue);
3a72645a 1635 AliCFEffGrid *efficiencymcPID = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
1636 AliCFEffGrid *efficiencymctrackinggeo = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack,fStepTrue);
1637 AliCFEffGrid *efficiencymcall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,fStepTrue);
1638
1639 AliCFEffGrid *efficiencyesdall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerESD),fStepMC,fStepTrue);
1640
1641 TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
1642 cefficiency->cd(1);
c2690925 1643 TH1D* efficiencymcPIDD = (TH1D *) efficiencymcPID->Project(ptpr);
3a72645a 1644 efficiencymcPIDD->SetTitle("");
1645 efficiencymcPIDD->SetStats(0);
1646 efficiencymcPIDD->SetMarkerStyle(25);
1647 efficiencymcPIDD->Draw();
c2690925 1648 TH1D* efficiencymctrackinggeoD = (TH1D *) efficiencymctrackinggeo->Project(ptpr);
3a72645a 1649 efficiencymctrackinggeoD->SetTitle("");
1650 efficiencymctrackinggeoD->SetStats(0);
1651 efficiencymctrackinggeoD->SetMarkerStyle(26);
1652 efficiencymctrackinggeoD->Draw("same");
c2690925 1653 TH1D* efficiencymcallD = (TH1D *) efficiencymcall->Project(ptpr);
3a72645a 1654 efficiencymcallD->SetTitle("");
1655 efficiencymcallD->SetStats(0);
1656 efficiencymcallD->SetMarkerStyle(27);
1657 efficiencymcallD->Draw("same");
c2690925 1658 TH1D* efficiencyesdallD = (TH1D *) efficiencyesdall->Project(ptpr);
3a72645a 1659 efficiencyesdallD->SetTitle("");
1660 efficiencyesdallD->SetStats(0);
1661 efficiencyesdallD->SetMarkerStyle(24);
1662 efficiencyesdallD->Draw("same");
1663 TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
1664 legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
1665 legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
1666 legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
1667 legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
1668 legeff->Draw("same");
c2690925 1669
1670 if(fBeamType==1) {
1671
1672 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1673 TAxis *cenaxisa = sparseafter->GetAxis(0);
1674 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1675 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1676 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1677 TAxis *cenaxisc = efficiencya->GetAxis(0);
1678 //Int_t nbbin = cenaxisb->GetNbins();
1679 //Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1680 //Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1681 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
1682 TString titlee("Efficiency_centrality_bin_");
1683 titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
1684 titlee += "_";
1685 titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
1686 TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1687 cefficiencye->Divide(2,1);
1688 cefficiencye->cd(1);
1689 gPad->SetLogy();
1690 cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
1691 cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
1692 TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr);
1693 TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr);
1694 CorrectFromTheWidth(afterefficiency);
1695 CorrectFromTheWidth(beforeefficiency);
1696 afterefficiency->SetStats(0);
1697 afterefficiency->SetTitle((const char*)titlee);
1698 afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1699 afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1700 afterefficiency->SetMarkerStyle(25);
1701 afterefficiency->SetMarkerColor(kBlack);
1702 afterefficiency->SetLineColor(kBlack);
1703 beforeefficiency->SetStats(0);
1704 beforeefficiency->SetTitle((const char*)titlee);
1705 beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1706 beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1707 beforeefficiency->SetMarkerStyle(24);
1708 beforeefficiency->SetMarkerColor(kBlue);
1709 beforeefficiency->SetLineColor(kBlue);
1710 afterefficiency->Draw();
1711 beforeefficiency->Draw("same");
1712 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1713 lega->AddEntry(beforeefficiency,"Before efficiency correction","p");
1714 lega->AddEntry(afterefficiency,"After efficiency correction","p");
1715 lega->Draw("same");
1716 cefficiencye->cd(2);
1717 cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1);
1718 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr);
1719 efficiencyDDproj->SetTitle("");
1720 efficiencyDDproj->SetStats(0);
1721 efficiencyDDproj->SetMarkerStyle(25);
1722 efficiencyDDproj->SetMarkerColor(2);
1723 efficiencyDDproj->Draw();
1724 cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]);
1725 TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr);
1726 efficiencyDDproja->SetTitle("");
1727 efficiencyDDproja->SetStats(0);
1728 efficiencyDDproja->SetMarkerStyle(26);
1729 efficiencyDDproja->SetMarkerColor(4);
1730 efficiencyDDproja->Draw("same");
1731
1732 }
1733 }
1734
1735
3a72645a 1736 }
1737
c04c80e6 1738 return result;
1739
1740}
3a72645a 1741
c04c80e6 1742//__________________________________________________________________________________
c2690925 1743TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum,Int_t i) const {
c04c80e6 1744 //
1745 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1746 // Give the final pt spectrum to be compared
1747 //
1748
c2690925 1749 if(fNEvents[i] > 0) {
1750
1751 Int_t ptpr;
1752 if(fBeamType==0) ptpr=0;
1753 if(fBeamType==1) ptpr=1;
c04c80e6 1754
c2690925 1755 TH1D* projection = spectrum->Projection(ptpr);
c04c80e6 1756 CorrectFromTheWidth(projection);
c2690925 1757 TGraphErrors *graphError = NormalizeTH1(projection,i);
c04c80e6 1758 return graphError;
1759
1760 }
1761
1762 return 0x0;
1763
1764
1765}
1766//__________________________________________________________________________________
c2690925 1767TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i) const {
c04c80e6 1768 //
1769 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1770 // Give the final pt spectrum to be compared
1771 //
c2690925 1772 if(fNEvents[i] > 0) {
c04c80e6 1773
c2690925 1774 Int_t ptpr;
1775 if(fBeamType==0) ptpr=0;
1776 if(fBeamType==1) ptpr=1;
1777
1778 TH1D* projection = (TH1D *) spectrum->Project(ptpr);
c04c80e6 1779 CorrectFromTheWidth(projection);
c2690925 1780 TGraphErrors *graphError = NormalizeTH1(projection,i);
c04c80e6 1781 return graphError;
1782
1783 }
1784
1785 return 0x0;
1786
1787
1788}
1789//__________________________________________________________________________________
c2690925 1790TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
1791 //
1792 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1793 // Give the final pt spectrum to be compared
1794 //
1795 if(fNEvents[i] > 0) {
1796
1797 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
1798 Double_t p = 0, dp = 0; Int_t point = 1;
1799 Double_t n = 0, dN = 0;
1800 Double_t nCorr = 0, dNcorr = 0;
1801 Double_t errdN = 0, errdp = 0;
1802 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
1803 point = ibin - input->GetXaxis()->GetFirst();
1804 p = input->GetXaxis()->GetBinCenter(ibin);
1805 dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
1806 n = input->GetBinContent(ibin);
1807 dN = input->GetBinError(ibin);
1808
1809 // New point
e156c3bb 1810 nCorr = 0.5 * 1./1.6 * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
c2690925 1811 errdN = 1./(2. * TMath::Pi() * p);
1812 errdp = 1./(2. * TMath::Pi() * p*p) * n;
e156c3bb 1813 dNcorr = 0.5 * 1./1.6 * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
c2690925 1814
1815 spectrumNormalized->SetPoint(point, p, nCorr);
1816 spectrumNormalized->SetPointError(point, dp, dNcorr);
1817 }
1818 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1819 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
1820 spectrumNormalized->SetMarkerStyle(22);
1821 spectrumNormalized->SetMarkerColor(kBlue);
1822 spectrumNormalized->SetLineColor(kBlue);
1823
1824 return spectrumNormalized;
1825
1826 }
1827
1828 return 0x0;
1829
1830
1831}
1832//__________________________________________________________________________________
1833TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) const {
c04c80e6 1834 //
1835 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1836 // Give the final pt spectrum to be compared
1837 //
c2690925 1838 if(normalization > 0) {
c04c80e6 1839
1840 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
1841 Double_t p = 0, dp = 0; Int_t point = 1;
1842 Double_t n = 0, dN = 0;
1843 Double_t nCorr = 0, dNcorr = 0;
1844 Double_t errdN = 0, errdp = 0;
1845 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
1846 point = ibin - input->GetXaxis()->GetFirst();
1847 p = input->GetXaxis()->GetBinCenter(ibin);
1848 dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
1849 n = input->GetBinContent(ibin);
1850 dN = input->GetBinError(ibin);
1851
1852 // New point
e156c3bb 1853 nCorr = 0.5 * 1./1.6 * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
c04c80e6 1854 errdN = 1./(2. * TMath::Pi() * p);
1855 errdp = 1./(2. * TMath::Pi() * p*p) * n;
e156c3bb 1856 dNcorr = 0.5 * 1./1.6 * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
c04c80e6 1857
1858 spectrumNormalized->SetPoint(point, p, nCorr);
1859 spectrumNormalized->SetPointError(point, dp, dNcorr);
1860 }
1861 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1862 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
1863 spectrumNormalized->SetMarkerStyle(22);
1864 spectrumNormalized->SetMarkerColor(kBlue);
1865 spectrumNormalized->SetLineColor(kBlue);
1866
1867 return spectrumNormalized;
1868
1869 }
1870
1871 return 0x0;
1872
1873
1874}
1875//____________________________________________________________
1876void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
1877 //
1878 // Set the container for a given step to the
1879 //
1880 if(!fCFContainers) fCFContainers = new TList;
1881 fCFContainers->AddAt(cont, type);
1882}
1883
1884//____________________________________________________________
1885AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
1886 //
1887 // Get Correction Framework Container for given type
1888 //
1889 if(!fCFContainers) return NULL;
1890 return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
1891}
c04c80e6 1892//____________________________________________________________
c2690925 1893AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Int_t positivenegative) {
c04c80e6 1894 //
3a72645a 1895 // Slice bin for a given source of electron
c2690925 1896 // nDim is the number of dimension the corrections are done
1897 // dimensions are the definition of the dimensions
1898 // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
1899 // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
c04c80e6 1900 //
1901
1902 Double_t *varMin = new Double_t[container->GetNVar()],
1903 *varMax = new Double_t[container->GetNVar()];
1904
1905 Double_t *binLimits;
1906 for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
1907
1908 binLimits = new Double_t[container->GetNBins(ivar)+1];
1909 container->GetBinLimits(ivar,binLimits);
c2690925 1910 varMin[ivar] = binLimits[0];
1911 varMax[ivar] = binLimits[container->GetNBins(ivar)];
1912 // source
1913 if(ivar == 4){
1914 if((source>= 0) && (source<container->GetNBins(ivar))) {
1915 varMin[ivar] = binLimits[source];
1916 varMax[ivar] = binLimits[source];
1917 }
3a72645a 1918 }
c2690925 1919 // charge
1920 if(ivar == 3) {
1921 if((positivenegative>= 0) && (positivenegative<container->GetNBins(ivar))) {
1922 varMin[ivar] = binLimits[positivenegative];
1923 varMax[ivar] = binLimits[positivenegative];
1924 }
3a72645a 1925 }
c2690925 1926
3a72645a 1927
c04c80e6 1928 delete[] binLimits;
1929 }
1930
1931 AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
1932 AddTemporaryObject(k);
1933 delete[] varMin; delete[] varMax;
1934
1935 return k;
1936
1937}
1938
1939//_________________________________________________________________________
3a72645a 1940THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions) const {
c04c80e6 1941 //
3a72645a 1942 // Slice correlation
c04c80e6 1943 //
1944
3a72645a 1945 Int_t ndimensions = correlationmatrix->GetNdimensions();
c2690925 1946 //printf("Number of dimension %d correlation map\n",ndimensions);
c04c80e6 1947 if(ndimensions < (2*nDim)) {
1948 AliError("Problem in the dimensions");
1949 return NULL;
1950 }
1951 Int_t ndimensionsContainer = (Int_t) ndimensions/2;
c2690925 1952 //printf("Number of dimension %d container\n",ndimensionsContainer);
c04c80e6 1953
1954 Int_t *dim = new Int_t[nDim*2];
1955 for(Int_t iter=0; iter < nDim; iter++){
1956 dim[iter] = dimensions[iter];
1957 dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
c2690925 1958 //printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
c04c80e6 1959 }
1960
3a72645a 1961 THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
c04c80e6 1962
1963 delete[] dim;
1964 return k;
1965
1966}
1967//___________________________________________________________________________
1968void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
1969 //
1970 // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
1971 //
1972
1973 TAxis *axis = h1->GetXaxis();
1974 Int_t nbinX = h1->GetNbinsX();
1975
1976 for(Int_t i = 1; i <= nbinX; i++) {
1977
1978 Double_t width = axis->GetBinWidth(i);
1979 Double_t content = h1->GetBinContent(i);
1980 Double_t error = h1->GetBinError(i);
1981 h1->SetBinContent(i,content/width);
1982 h1->SetBinError(i,error/width);
1983 }
1984
1985}
1986//____________________________________________________________
1987void AliHFEspectrum::AddTemporaryObject(TObject *o){
1988 //
1989 // Emulate garbage collection on class level
1990 //
1991 if(!fTemporaryObjects) {
1992 fTemporaryObjects= new TList;
1993 fTemporaryObjects->SetOwner();
1994 }
1995 fTemporaryObjects->Add(o);
1996}
1997
1998//____________________________________________________________
1999void AliHFEspectrum::ClearObject(TObject *o){
2000 //
2001 // Do a safe deletion
2002 //
2003 if(fTemporaryObjects){
2004 if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
2005 delete o;
2006 } else{
2007 // just delete
2008 delete o;
2009 }
2010}
2011//_________________________________________________________________________
c2690925 2012TObject* AliHFEspectrum::GetSpectrum(const AliCFContainer * const c, Int_t step) {
245877d0 2013 AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
c04c80e6 2014 return data;
2015}
2016//_________________________________________________________________________
c2690925 2017TObject* AliHFEspectrum::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
c04c80e6 2018 //
2019 // Create efficiency grid and calculate efficiency
2020 // of step to step0
2021 //
2022 TString name("eff");
2023 name += step;
2024 name+= step0;
2025 AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
2026 eff->CalculateEfficiency(step,step0);
2027 return eff;
2028}
c2690925 2029
2030//____________________________________________________________________________
2031THnSparse* AliHFEspectrum::GetWeights(){
2032
2033 //
2034 // Measured D->e based weighting factors
2035 //
2036
2037 const Int_t nDim=1;
2038 Int_t nBin[nDim] = {44};
2039 const Double_t kPtbound[2] = {0.1, 20.};
2040
2041 Double_t* binEdges[nDim];
2042 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
2043
2044 fWeightCharm = new THnSparseF("weightHisto", "weiting factor; pt[GeV/c]", nDim, nBin);
2045 for(Int_t idim = 0; idim < nDim; idim++)
2046 fWeightCharm->SetBinEdges(idim, binEdges[idim]);
2047 /* //fit
2048 Double_t weight[44]={ 2.20939,2.1742,2.13569,2.09369,2.04806,1.99869,1.94552,1.88854,1.82782,1.76349,1.69577,1.62497,1.55149,1.4758,1.39848,1.32016,1.24153,1.16328,1.08615,1.01081,0.937916,0.868033,0.801642,0.739117,0.680724,0.626621,0.576859,0.5314,0.490125,0.452854,0.419361,0.389388,0.362661,0.338898,0.317822,0.299165,0.282674,0.268115,0.255271,0.243946,0.233964,0.225166,0.217413,0.210578};*/
2049 //points
2050 Double_t weight[44]={ 0.443336,0.437069,0.424299,0.41365,0.410649,0.385012,0.364252,0.360181,0.347073,0.354445,0.369578,0.358253,0.381115,0.412863,0.451792,0.47061,0.544857,0.566261,0.656843,0.669739,0.725247,0.779511,0.750385,0.777815,0.673391,0.694613,0.598251,0.513847,0.45729,0.468098,0.379422,0.358718,0.403793,0.392902,0.349981,0.27205,0.302447,0.279313,0.258986,0.14738,0.414487,0.0552302,0.0184101,0};
2051 Double_t pt[44]={0.106398,0.120014,0.135371,0.152694,0.172234,0.194274,0.219135,0.247177,0.278807,0.314485,0.354729,0.400122,0.451324,0.509078,0.574223,0.647704,0.730589,0.824079,0.929534,1.04848,1.18265,1.33399,1.5047,1.69725,1.91444,2.15943,2.43576,2.74745,3.09904,3.49561,3.94293,4.44749,5.01662,5.65858,6.38268,7.19945,8.12074,9.15992,10.3321,11.6542,13.1456,14.8278,16.7252,18.8655};
2052 Double_t dataE[1];
2053 for(int i=0; i<44; i++){
2054 dataE[0]=pt[i];
2055 fWeightCharm->Fill(dataE,weight[i]);
2056 }
2057 Int_t* ibins[nDim];
2058 for(Int_t ivar = 0; ivar < nDim; ivar++)
2059 ibins[ivar] = new Int_t[nBin[ivar] + 1];
2060 fWeightCharm->SetBinError(ibins[0],0);
2061
2062 return fWeightCharm;
2063}