]> 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
Moving to ansi isfinite
[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
6c70d827 453 TH1D unfoldingspectrac[fNCentralityBinAtTheEnd];
454 TGraphErrors unfoldingspectracn[fNCentralityBinAtTheEnd];
455 TH1D correctedspectrac[fNCentralityBinAtTheEnd];
456 TGraphErrors correctedspectracn[fNCentralityBinAtTheEnd];
c2690925 457
458
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
630 }
631
632
633
634
635 return kTRUE;
636}
637
638//____________________________________________________________
639Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
640 //
641 // Correct the spectrum for efficiency and unfolding for beauty analysis
642 // with both method and compare
643 //
644
645 gStyle->SetPalette(1);
646 gStyle->SetOptStat(1111);
647 gStyle->SetPadBorderMode(0);
648 gStyle->SetCanvasColor(10);
649 gStyle->SetPadLeftMargin(0.13);
650 gStyle->SetPadRightMargin(0.13);
651
3a72645a 652 ///////////////////////////
653 // Check initialization
654 ///////////////////////////
c04c80e6 655
3a72645a 656 if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
657 AliInfo("You have to init before");
658 return kFALSE;
659 }
660
661 if((fStepTrue == 0) && (fStepMC == 0) && (fStepData == 0)) {
662 AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
663 return kFALSE;
664 }
665
c2690925 666 SetNumberOfIteration(10);
3a72645a 667 SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
668
669 AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
670 //////////////////////////////////
671 // Subtract hadron background
672 /////////////////////////////////
67fe7bd0 673 AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
3a72645a 674 if(subtractcontamination) {
675 dataspectrumaftersubstraction = SubtractBackground(kTRUE);
676 dataGridAfterFirstSteps = dataspectrumaftersubstraction;
677 }
678
c2690925 679 /////////////////////////////////////////////////////////////////////////////////////////
680 // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
681 /////////////////////////////////////////////////////////////////////////////////////////
682
683
3a72645a 684 ////////////////////////////////////////////////
685 // Correct for TPC efficiency from V0
686 ///////////////////////////////////////////////
687 AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
688 AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
689 if(dataContainerV0){
690 dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
691 dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
67fe7bd0 692 }
c04c80e6 693
3a72645a 694 ///////////////
c04c80e6 695 // Unfold
3a72645a 696 //////////////
697 TList *listunfolded = Unfold(dataGridAfterFirstSteps);
c04c80e6 698 if(!listunfolded){
699 printf("Unfolded failed\n");
3a72645a 700 return kFALSE;
c04c80e6 701 }
702 THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
703 THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
704 if(!correctedspectrum){
705 AliError("No corrected spectrum\n");
3a72645a 706 return kFALSE;
c04c80e6 707 }
67fe7bd0 708 if(!residualspectrum){
c04c80e6 709 AliError("No residul spectrum\n");
3a72645a 710 return kFALSE;
c04c80e6 711 }
67fe7bd0 712
3a72645a 713 /////////////////////
c04c80e6 714 // Simply correct
3a72645a 715 ////////////////////
716 AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
67fe7bd0 717
3a72645a 718
67fe7bd0 719 //////////
c04c80e6 720 // Plot
721 //////////
3a72645a 722 if(fDebugLevel > 0.0) {
723
724 TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
725 ccorrected->Divide(2,1);
726 ccorrected->cd(1);
727 gPad->SetLogy();
728 TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
729 correctedspectrumD->SetTitle("");
730 correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
731 correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
732 correctedspectrumD->SetMarkerStyle(26);
733 correctedspectrumD->SetMarkerColor(kBlue);
734 correctedspectrumD->SetLineColor(kBlue);
735 correctedspectrumD->Draw("AP");
736 TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
737 alltogetherspectrumD->SetTitle("");
738 alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
739 alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
740 alltogetherspectrumD->SetMarkerStyle(25);
741 alltogetherspectrumD->SetMarkerColor(kBlack);
742 alltogetherspectrumD->SetLineColor(kBlack);
743 alltogetherspectrumD->Draw("P");
744 TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
745 legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
746 legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
747 legcorrected->Draw("same");
748 ccorrected->cd(2);
749 TH1D *correctedTH1D = correctedspectrum->Projection(0);
6555e2ad 750 TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(0);
3a72645a 751 TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
752 ratiocorrected->SetName("ratiocorrected");
753 ratiocorrected->SetTitle("");
754 ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
755 ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
756 ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
757 ratiocorrected->SetStats(0);
758 ratiocorrected->Draw();
759
760
761 // Dump to file if needed
762
763 if(fDumpToFile) {
764
765 TFile *out = new TFile("finalSpectrum.root","recreate");
766 out->cd();
767 //
768 correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
769 correctedspectrumD->Write();
770 alltogetherspectrumD->SetName("AlltogetherSpectrum");
771 alltogetherspectrumD->Write();
772 ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
773 ratiocorrected->Write();
774 //
775 correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
776 correctedspectrum->Write();
777 alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
778 alltogetherCorrection->Write();
779 //
780 out->Close(); delete out;
781 }
782
783
784 }
785
786
787
788
789 return kTRUE;
790}
c2690925 791
3a72645a 792//____________________________________________________________
793AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
794 //
795 // Apply background subtraction
796 //
797
798 // Raw spectrum
799 AliCFContainer *dataContainer = GetContainer(kDataContainer);
800 if(!dataContainer){
801 AliError("Data Container not available");
802 return NULL;
803 }
c2690925 804 printf("Step data: %d\n",fStepData);
3a72645a 805 AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
806
807 AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
808 dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
809
810 // Background Estimate
811 AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
812 if(!backgroundContainer){
813 AliError("MC background container not found");
814 return NULL;
815 }
816
c2690925 817 Int_t stepbackground = 1; // 2 for !fInclusiveSpectrum analysis(old method)
3a72645a 818 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
819
c2690925 820 if(!fInclusiveSpectrum){
821 //Background subtraction for IP analysis
822 if(fIPanaHadronBgSubtract){
823 // Hadron background
824 printf("Hadron background for IP analysis subtracted!\n");
825 Int_t* bins=new Int_t[1];
826 TH1D* htemp = (TH1D *) fHadronEffbyIPcut->Projection(0);
827 bins[0]=htemp->GetNbinsX();
828 AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",1,bins);
829 hbgContainer->SetGrid(fHadronEffbyIPcut);
830 backgroundGrid->Multiply(hbgContainer,1);
831 spectrumSubtracted->Add(backgroundGrid,-1.0);
832 }
833 if(fIPanaCharmBgSubtract){
834 // Charm background
835 printf("Charm background for IP analysis subtracted!\n");
836 AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
837 spectrumSubtracted->Add(charmbgContainer,-1.0);
838 }
839 if(fIPanaConversionBgSubtract){
840 // Conversion background
841 AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground();
842 spectrumSubtracted->Add(conversionbgContainer,-1.0);
843 printf("Conversion background subtraction is preliminary!\n");
844 }
845 if(fIPanaNonHFEBgSubtract){
846 // NonHFE background
847 AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground();
848 spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
849 printf("Non HFE background subtraction is preliminary!\n");
850 }
851 }
852 else{
853 // Subtract
854 spectrumSubtracted->Add(backgroundGrid,-1.0);
855 }
856
3a72645a 857 if(setBackground){
858 if(fBackground) delete fBackground;
859 fBackground = backgroundGrid;
860 } else delete backgroundGrid;
861
862
863 if(fDebugLevel > 0) {
c2690925 864
865 Int_t ptpr;
866 if(fBeamType==0) ptpr=0;
867 if(fBeamType==1) ptpr=1;
67fe7bd0 868
869 TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
3a72645a 870 cbackgroundsubtraction->Divide(3,1);
67fe7bd0 871 cbackgroundsubtraction->cd(1);
3a72645a 872 gPad->SetLogy();
c2690925 873 TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptpr);
874 TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
67fe7bd0 875 CorrectFromTheWidth(measuredTH1Daftersubstraction);
876 CorrectFromTheWidth(measuredTH1Dbeforesubstraction);
877 measuredTH1Daftersubstraction->SetStats(0);
878 measuredTH1Daftersubstraction->SetTitle("");
879 measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
880 measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
881 measuredTH1Daftersubstraction->SetMarkerStyle(25);
882 measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
883 measuredTH1Daftersubstraction->SetLineColor(kBlack);
884 measuredTH1Dbeforesubstraction->SetStats(0);
885 measuredTH1Dbeforesubstraction->SetTitle("");
886 measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
887 measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
888 measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
889 measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
890 measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
891 measuredTH1Daftersubstraction->Draw();
892 measuredTH1Dbeforesubstraction->Draw("same");
893 TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
3a72645a 894 legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
895 legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
67fe7bd0 896 legsubstraction->Draw("same");
897 cbackgroundsubtraction->cd(2);
3a72645a 898 gPad->SetLogy();
67fe7bd0 899 TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
900 ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
901 ratiomeasuredcontamination->SetTitle("");
902 ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
903 ratiomeasuredcontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
c2690925 904 ratiomeasuredcontamination->Sumw2();
67fe7bd0 905 ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
906 ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
907 ratiomeasuredcontamination->SetStats(0);
908 ratiomeasuredcontamination->SetMarkerStyle(26);
909 ratiomeasuredcontamination->SetMarkerColor(kBlack);
910 ratiomeasuredcontamination->SetLineColor(kBlack);
c2690925 911 for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
912 ratiomeasuredcontamination->SetBinError(k+1,0.0);
913 }
914 ratiomeasuredcontamination->Draw("P");
3a72645a 915 cbackgroundsubtraction->cd(3);
c2690925 916 TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptpr);
3a72645a 917 CorrectFromTheWidth(measuredTH1background);
918 measuredTH1background->SetStats(0);
919 measuredTH1background->SetTitle("");
920 measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
921 measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]");
922 measuredTH1background->SetMarkerStyle(26);
923 measuredTH1background->SetMarkerColor(kBlack);
924 measuredTH1background->SetLineColor(kBlack);
925 measuredTH1background->Draw();
c2690925 926
927 if(fBeamType==1) {
928
929 TCanvas * cbackgrounde = new TCanvas("BackgroundSubtraction_allspectra","BackgroundSubtraction_allspectra",1000,700);
930 cbackgrounde->Divide(2,1);
931 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
932 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
933
934 THnSparseF* sparsesubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
935 TAxis *cenaxisa = sparsesubtracted->GetAxis(0);
936 THnSparseF* sparsebefore = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
937 TAxis *cenaxisb = sparsebefore->GetAxis(0);
938 Int_t nbbin = cenaxisb->GetNbins();
939 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
940 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
941 for(Int_t binc = 0; binc < nbbin; binc++){
942 TString titlee("BackgroundSubtraction_centrality_bin_");
943 titlee += binc;
944 TCanvas * cbackground = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
945 cbackground->Divide(2,1);
946 cbackground->cd(1);
947 gPad->SetLogy();
948 cenaxisa->SetRange(binc+1,binc+1);
949 cenaxisb->SetRange(binc+1,binc+1);
950 TH1D *aftersubstraction = (TH1D *) sparsesubtracted->Projection(1);
951 TH1D *beforesubstraction = (TH1D *) sparsebefore->Projection(1);
952 CorrectFromTheWidth(aftersubstraction);
953 CorrectFromTheWidth(beforesubstraction);
954 aftersubstraction->SetStats(0);
955 aftersubstraction->SetTitle((const char*)titlee);
956 aftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
957 aftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
958 aftersubstraction->SetMarkerStyle(25);
959 aftersubstraction->SetMarkerColor(kBlack);
960 aftersubstraction->SetLineColor(kBlack);
961 beforesubstraction->SetStats(0);
962 beforesubstraction->SetTitle((const char*)titlee);
963 beforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
964 beforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
965 beforesubstraction->SetMarkerStyle(24);
966 beforesubstraction->SetMarkerColor(kBlue);
967 beforesubstraction->SetLineColor(kBlue);
968 aftersubstraction->Draw();
969 beforesubstraction->Draw("same");
970 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
971 lega->AddEntry(beforesubstraction,"With hadron contamination","p");
972 lega->AddEntry(aftersubstraction,"Without hadron contamination ","p");
973 lega->Draw("same");
974 cbackgrounde->cd(1);
975 gPad->SetLogy();
976 TH1D *aftersubtractionn = (TH1D *) aftersubstraction->Clone();
977 aftersubtractionn->SetMarkerStyle(stylee[binc]);
978 aftersubtractionn->SetMarkerColor(colorr[binc]);
979 if(binc==0) aftersubtractionn->Draw();
980 else aftersubtractionn->Draw("same");
981 legtotal->AddEntry(aftersubtractionn,(const char*) titlee,"p");
982 cbackgrounde->cd(2);
983 gPad->SetLogy();
984 TH1D *aftersubtractionng = (TH1D *) aftersubstraction->Clone();
985 aftersubtractionng->SetMarkerStyle(stylee[binc]);
986 aftersubtractionng->SetMarkerColor(colorr[binc]);
987 if(fNEvents[binc] > 0.0) aftersubtractionng->Scale(1/(Double_t)fNEvents[binc]);
988 if(binc==0) aftersubtractionng->Draw();
989 else aftersubtractionng->Draw("same");
990 legtotalg->AddEntry(aftersubtractionng,(const char*) titlee,"p");
991 cbackground->cd(2);
992 TH1D* ratiocontamination = (TH1D*)beforesubstraction->Clone();
993 ratiocontamination->SetName("ratiocontamination");
994 ratiocontamination->SetTitle((const char*)titlee);
995 ratiocontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
996 ratiocontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
997 ratiocontamination->Add(aftersubstraction,-1.0);
998 ratiocontamination->Divide(beforesubstraction);
999 Int_t totalbin = ratiocontamination->GetXaxis()->GetNbins();
1000 for(Int_t nbinpt = 0; nbinpt < totalbin; nbinpt++) {
1001 ratiocontamination->SetBinError(nbinpt+1,0.0);
1002 }
1003 ratiocontamination->SetStats(0);
1004 ratiocontamination->SetMarkerStyle(26);
1005 ratiocontamination->SetMarkerColor(kBlack);
1006 ratiocontamination->SetLineColor(kBlack);
1007 ratiocontamination->Draw("P");
1008
1009 }
1010
1011 cbackgrounde->cd(1);
1012 legtotal->Draw("same");
1013 cbackgrounde->cd(2);
1014 legtotalg->Draw("same");
1015
1016 cenaxisa->SetRange(0,nbbin);
1017 cenaxisb->SetRange(0,nbbin);
1018
1019 }
1020
3a72645a 1021
67fe7bd0 1022 }
1023
c04c80e6 1024
3a72645a 1025 return spectrumSubtracted;
c04c80e6 1026}
c2690925 1027
1028//____________________________________________________________
1029AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
1030 //
1031 // calculate charm background
1032 //
1033
1034 Double_t evtnorm=0;
1035 if(fNMCEvents) evtnorm= double(fNEvents[0])/double(fNMCEvents);
1036
1037 AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
1038 if(!mcContainer){
1039 AliError("MC Container not available");
1040 return NULL;
1041 }
1042
1043 if(!fCorrelation){
1044 AliError("No Correlation map available");
1045 return NULL;
1046 }
1047
1048 Int_t nDim = 1;
1049 AliCFDataGrid *charmBackgroundGrid= 0x0;
1050 charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC+2);
1051 TH1D *tmphisto = (TH1D *) charmBackgroundGrid->Project(0);
1052 Int_t* bins=new Int_t[1];
1053 bins[0]=tmphisto->GetNbinsX();
1054 AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",1,bins);
1055 weightedCharmContainer->SetGrid(GetWeights());
1056 charmBackgroundGrid->Multiply(weightedCharmContainer,evtnorm);
1057
1058 // Efficiency (set efficiency to 1 for only folding)
1059 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1060 efficiencyD->CalculateEfficiency(0,0);
1061
1062 // Folding
1063 AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
1064 folding.SetMaxNumberOfIterations(1);
1065 folding.Unfold();
1066
1067 // Results
1068 THnSparse* result1= folding.GetEstMeasured(); // folded spectra
1069 THnSparse* result=(THnSparse*)result1->Clone();
1070 charmBackgroundGrid->SetGrid(result);
1071
1072 return charmBackgroundGrid;
1073}
1074
1075//____________________________________________________________
1076AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
1077 //
1078 // calculate conversion background
1079 //
1080
1081 Double_t evtnorm[1];
1082 if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
1083 printf("check event!!! %lf \n",evtnorm[0]);
1084
1085 // Background Estimate
1086 AliCFContainer *backgroundContainer = GetContainer(kNonHFEBackgroundData);
1087 if(!backgroundContainer){
1088 AliError("MC background container not found");
1089 return NULL;
1090 }
1091
1092 Int_t stepbackground = 2;
1093 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground);
1094 backgroundGrid->Scale(evtnorm);
1095
1096 new TCanvas;
1097 TH1D *histodraw111 = (TH1D *) backgroundGrid->Project(0);
1098 CorrectFromTheWidth(histodraw111);
1099 histodraw111->SetMarkerStyle(20);
1100 histodraw111->SetLineColor(4);
1101 histodraw111->SetMarkerColor(4);
1102 histodraw111->Draw("p");
1103
1104
1105 return backgroundGrid;
1106}
1107
1108
1109//____________________________________________________________
1110AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
1111 //
1112 // calculate non HFE background
1113 //
1114
1115 if(fNonHFEbgMethod2){
1116 Double_t evtnorm=0;
1117 if(fNMCEvents) evtnorm= double(fNEvents[0])/double(fNMCEvents);
1118
1119 AliCFContainer *mcContainer = GetContainer(kMCContainerNonHFEMC);
1120 if(!mcContainer){
1121 AliError("MC Container not available");
1122 return NULL;
1123 }
1124
1125 if(!fCorrelation){
1126 AliError("No Correlation map available");
1127 return NULL;
1128 }
1129
1130 Int_t nDim = 1;
1131 AliCFDataGrid *nonHFEBackgroundGrid= 0x0;
1132 nonHFEBackgroundGrid = new AliCFDataGrid("nonHFEBackgroundGrid","nonHFEBackgroundGrid",*mcContainer, fStepMC+2);
1133 TH1D *tmphisto = (TH1D *) nonHFEBackgroundGrid->Project(0);
1134 Int_t* bins=new Int_t[1];
1135 bins[0]=tmphisto->GetNbinsX();
1136 AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
1137 weightedNonHFEContainer->SetGrid(GetWeights()); // we have to use different weithing function
1138 nonHFEBackgroundGrid->Multiply(weightedNonHFEContainer,evtnorm);
1139
1140 // Efficiency (set efficiency to 1 for only folding)
1141 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1142 efficiencyD->CalculateEfficiency(0,0);
1143
1144 // Folding
1145 AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),nonHFEBackgroundGrid->GetGrid(),nonHFEBackgroundGrid->GetGrid());
1146 folding.SetMaxNumberOfIterations(1);
1147 folding.Unfold();
1148
1149 // Results
1150 THnSparse* result1= folding.GetEstMeasured(); // folded spectra
1151 THnSparse* result=(THnSparse*)result1->Clone();
1152 nonHFEBackgroundGrid->SetGrid(result);
1153
1154 return nonHFEBackgroundGrid;
1155 }
1156 else{
1157 Double_t evtnorm[1];
1158 if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
1159
1160 // Background Estimate
1161 AliCFContainer *backgroundContainer = GetContainer(kNonHFEBackgroundData);
1162 if(!backgroundContainer){
1163 AliError("MC background container not found");
1164 return NULL;
1165 }
1166
1167 Int_t stepbackground = 3;
1168 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
1169 backgroundGrid->Scale(evtnorm);
1170
1171 new TCanvas;
1172 TH1D *histodraw222 = (TH1D *) backgroundGrid->Project(0);
1173 CorrectFromTheWidth(histodraw222);
1174 histodraw222->SetMarkerStyle(20);
1175 histodraw222->SetLineColor(4);
1176 histodraw222->SetMarkerColor(4);
1177 histodraw222->Draw("p");
1178
1179 return backgroundGrid;
1180 }
1181
1182}
1183
1184//____________________________________________________________
1185AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
1186
1187 //
1188 // Apply TPC pid efficiency correction from parametrisation
1189 //
1190
1191 // Data in the right format
1192 AliCFDataGrid *dataGrid = 0x0;
1193 if(bgsubpectrum) {
1194 dataGrid = bgsubpectrum;
1195 }
1196 else {
1197
1198 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1199 if(!dataContainer){
1200 AliError("Data Container not available");
1201 return NULL;
1202 }
1203
1204 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1205 }
1206
1207 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1208 result->SetName("ParametrizedEfficiencyBefore");
1209 THnSparse *h = result->GetGrid();
1210 Int_t nbdimensions = h->GetNdimensions();
1211 //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
1212
1213 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1214 if(!dataContainer){
1215 AliError("Data Container not available");
1216 return NULL;
1217 }
1218 AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
1219 dataContainerbis->Add(dataContainerbis,-1.0);
1220
1221
1222 Int_t* coord = new Int_t[nbdimensions];
1223 memset(coord, 0, sizeof(Int_t) * nbdimensions);
1224 Double_t* points = new Double_t[nbdimensions];
1225
1226
1227 ULong64_t nEntries = h->GetNbins();
1228 for (ULong64_t i = 0; i < nEntries; ++i) {
1229
1230 Double_t value = h->GetBinContent(i, coord);
1231 //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
1232 //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
1233
1234 // Get the bin co-ordinates given an coord
1235 for (Int_t j = 0; j < nbdimensions; ++j)
1236 points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1237
1238 if (!fEfficiencyFunction->IsInside(points))
1239 continue;
1240 TF1::RejectPoint(kFALSE);
1241
1242 // Evaulate function at points
1243 Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1244 //printf("Value efficiency is %f\n",valueEfficiency);
1245
1246 if(valueEfficiency > 0.0) {
1247 h->SetBinContent(coord,value/valueEfficiency);
1248 dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1249 }
1250 Double_t error = h->GetBinError(i);
1251 h->SetBinError(coord,error/valueEfficiency);
1252 dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1253
1254
1255 }
1256
6c70d827 1257 delete[] coord;
1258 delete[] points;
1259
c2690925 1260 AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1261
1262 if(fDebugLevel > 0) {
1263
1264 TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700);
1265 cEfficiencyParametrized->Divide(2,1);
1266 cEfficiencyParametrized->cd(1);
1267 TH1D *afterE = (TH1D *) resultt->Project(0);
1268 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1269 CorrectFromTheWidth(afterE);
1270 CorrectFromTheWidth(beforeE);
1271 afterE->SetStats(0);
1272 afterE->SetTitle("");
1273 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1274 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1275 afterE->SetMarkerStyle(25);
1276 afterE->SetMarkerColor(kBlack);
1277 afterE->SetLineColor(kBlack);
1278 beforeE->SetStats(0);
1279 beforeE->SetTitle("");
1280 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1281 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1282 beforeE->SetMarkerStyle(24);
1283 beforeE->SetMarkerColor(kBlue);
1284 beforeE->SetLineColor(kBlue);
1285 gPad->SetLogy();
1286 afterE->Draw();
1287 beforeE->Draw("same");
1288 TLegend *legefficiencyparametrized = new TLegend(0.4,0.6,0.89,0.89);
1289 legefficiencyparametrized->AddEntry(beforeE,"Before Efficiency correction","p");
1290 legefficiencyparametrized->AddEntry(afterE,"After Efficiency correction","p");
1291 legefficiencyparametrized->Draw("same");
1292 cEfficiencyParametrized->cd(2);
1293 fEfficiencyFunction->Draw();
1294 //cEfficiencyParametrized->cd(3);
1295 //TH1D *ratioefficiency = (TH1D *) beforeE->Clone();
1296 //ratioefficiency->Divide(afterE);
1297 //ratioefficiency->Draw();
1298
1299
1300 }
1301
1302
1303 return resultt;
1304
1305}
c04c80e6 1306//____________________________________________________________
3a72645a 1307AliCFDataGrid *AliHFEspectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
1308
c04c80e6 1309 //
3a72645a 1310 // Apply TPC pid efficiency correction from V0
c04c80e6 1311 //
1312
3a72645a 1313 AliCFContainer *v0Container = GetContainer(kDataContainerV0);
1314 if(!v0Container){
1315 AliError("V0 Container not available");
c04c80e6 1316 return NULL;
1317 }
3a72645a 1318
1319 // Efficiency
1320 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
1321 efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
1322
1323 // Data in the right format
1324 AliCFDataGrid *dataGrid = 0x0;
1325 if(bgsubpectrum) {
1326 dataGrid = bgsubpectrum;
c04c80e6 1327 }
3a72645a 1328 else {
1329
1330 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1331 if(!dataContainer){
1332 AliError("Data Container not available");
1333 return NULL;
1334 }
c04c80e6 1335
3a72645a 1336 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1337 }
c04c80e6 1338
3a72645a 1339 // Correct
1340 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1341 result->ApplyEffCorrection(*efficiencyD);
c04c80e6 1342
3a72645a 1343 if(fDebugLevel > 0) {
c2690925 1344
1345 Int_t ptpr;
1346 if(fBeamType==0) ptpr=0;
1347 if(fBeamType==1) ptpr=1;
3a72645a 1348
1349 TCanvas * cV0Efficiency = new TCanvas("V0Efficiency","V0Efficiency",1000,700);
1350 cV0Efficiency->Divide(2,1);
1351 cV0Efficiency->cd(1);
c2690925 1352 TH1D *afterE = (TH1D *) result->Project(ptpr);
1353 TH1D *beforeE = (TH1D *) dataGrid->Project(ptpr);
3a72645a 1354 afterE->SetStats(0);
1355 afterE->SetTitle("");
1356 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1357 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1358 afterE->SetMarkerStyle(25);
1359 afterE->SetMarkerColor(kBlack);
1360 afterE->SetLineColor(kBlack);
1361 beforeE->SetStats(0);
1362 beforeE->SetTitle("");
1363 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1364 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1365 beforeE->SetMarkerStyle(24);
1366 beforeE->SetMarkerColor(kBlue);
1367 beforeE->SetLineColor(kBlue);
1368 afterE->Draw();
1369 beforeE->Draw("same");
1370 TLegend *legV0efficiency = new TLegend(0.4,0.6,0.89,0.89);
1371 legV0efficiency->AddEntry(beforeE,"Before Efficiency correction","p");
1372 legV0efficiency->AddEntry(afterE,"After Efficiency correction","p");
1373 legV0efficiency->Draw("same");
1374 cV0Efficiency->cd(2);
c2690925 1375 TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(ptpr);
3a72645a 1376 efficiencyDproj->SetTitle("");
1377 efficiencyDproj->SetStats(0);
1378 efficiencyDproj->SetMarkerStyle(25);
1379 efficiencyDproj->Draw();
c04c80e6 1380
c2690925 1381 if(fBeamType==1) {
1382
1383 TCanvas * cV0Efficiencye = new TCanvas("V0Efficiency_allspectra","V0Efficiency_allspectra",1000,700);
1384 cV0Efficiencye->Divide(2,1);
1385 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1386 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1387
1388 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1389 TAxis *cenaxisa = sparseafter->GetAxis(0);
1390 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1391 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1392 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1393 TAxis *cenaxisc = efficiencya->GetAxis(0);
1394 Int_t nbbin = cenaxisb->GetNbins();
1395 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1396 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1397 for(Int_t binc = 0; binc < nbbin; binc++){
1398 TString titlee("V0Efficiency_centrality_bin_");
1399 titlee += binc;
1400 TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1401 ccV0Efficiency->Divide(2,1);
1402 ccV0Efficiency->cd(1);
1403 gPad->SetLogy();
1404 cenaxisa->SetRange(binc+1,binc+1);
1405 cenaxisb->SetRange(binc+1,binc+1);
1406 cenaxisc->SetRange(binc+1,binc+1);
1407 TH1D *aftere = (TH1D *) sparseafter->Projection(1);
1408 TH1D *beforee = (TH1D *) sparsebefore->Projection(1);
1409 CorrectFromTheWidth(aftere);
1410 CorrectFromTheWidth(beforee);
1411 aftere->SetStats(0);
1412 aftere->SetTitle((const char*)titlee);
1413 aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1414 aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1415 aftere->SetMarkerStyle(25);
1416 aftere->SetMarkerColor(kBlack);
1417 aftere->SetLineColor(kBlack);
1418 beforee->SetStats(0);
1419 beforee->SetTitle((const char*)titlee);
1420 beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1421 beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1422 beforee->SetMarkerStyle(24);
1423 beforee->SetMarkerColor(kBlue);
1424 beforee->SetLineColor(kBlue);
1425 aftere->Draw();
1426 beforee->Draw("same");
1427 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1428 lega->AddEntry(beforee,"Before correction","p");
1429 lega->AddEntry(aftere,"After correction","p");
1430 lega->Draw("same");
1431 cV0Efficiencye->cd(1);
1432 gPad->SetLogy();
1433 TH1D *afteree = (TH1D *) aftere->Clone();
1434 afteree->SetMarkerStyle(stylee[binc]);
1435 afteree->SetMarkerColor(colorr[binc]);
1436 if(binc==0) afteree->Draw();
1437 else afteree->Draw("same");
1438 legtotal->AddEntry(afteree,(const char*) titlee,"p");
1439 cV0Efficiencye->cd(2);
1440 gPad->SetLogy();
1441 TH1D *aftereeu = (TH1D *) aftere->Clone();
1442 aftereeu->SetMarkerStyle(stylee[binc]);
1443 aftereeu->SetMarkerColor(colorr[binc]);
1444 if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]);
1445 if(binc==0) aftereeu->Draw();
1446 else aftereeu->Draw("same");
1447 legtotalg->AddEntry(aftereeu,(const char*) titlee,"p");
1448 ccV0Efficiency->cd(2);
1449 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1);
1450 efficiencyDDproj->SetTitle("");
1451 efficiencyDDproj->SetStats(0);
1452 efficiencyDDproj->SetMarkerStyle(25);
1453 efficiencyDDproj->Draw();
1454
1455 }
1456
1457 cV0Efficiencye->cd(1);
1458 legtotal->Draw("same");
1459 cV0Efficiencye->cd(2);
1460 legtotalg->Draw("same");
1461
1462 cenaxisa->SetRange(0,nbbin);
1463 cenaxisb->SetRange(0,nbbin);
1464 cenaxisc->SetRange(0,nbbin);
1465
1466 }
3a72645a 1467
1468 }
1469
1470
1471 return result;
1472
1473}
c04c80e6 1474//____________________________________________________________
3a72645a 1475TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
c04c80e6 1476
1477 //
1478 // Unfold and eventually correct for efficiency the bgsubspectrum
1479 //
1480
3a72645a 1481 AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
c04c80e6 1482 if(!mcContainer){
1483 AliError("MC Container not available");
1484 return NULL;
1485 }
1486
1487 if(!fCorrelation){
1488 AliError("No Correlation map available");
1489 return NULL;
1490 }
1491
3a72645a 1492 // Data
c04c80e6 1493 AliCFDataGrid *dataGrid = 0x0;
1494 if(bgsubpectrum) {
c04c80e6 1495 dataGrid = bgsubpectrum;
c04c80e6 1496 }
1497 else {
1498
1499 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1500 if(!dataContainer){
1501 AliError("Data Container not available");
1502 return NULL;
1503 }
1504
3a72645a 1505 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
c04c80e6 1506 }
1507
c04c80e6 1508 // Guessed
3a72645a 1509 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
c04c80e6 1510 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
1511
1512 // Efficiency
3a72645a 1513 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
c04c80e6 1514 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1515
1516 // Unfold
1517
3a72645a 1518 AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
c2690925 1519 //unfolding.SetUseCorrelatedErrors();
1520 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
c04c80e6 1521 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
e156c3bb 1522 if(fSetSmoothing) unfolding.UseSmoothing();
c04c80e6 1523 unfolding.Unfold();
1524
1525 // Results
1526 THnSparse* result = unfolding.GetUnfolded();
1527 THnSparse* residual = unfolding.GetEstMeasured();
1528 TList *listofresults = new TList;
1529 listofresults->SetOwner();
1530 listofresults->AddAt((THnSparse*)result->Clone(),0);
1531 listofresults->AddAt((THnSparse*)residual->Clone(),1);
1532
3a72645a 1533 if(fDebugLevel > 0) {
c2690925 1534
1535 Int_t ptpr;
1536 if(fBeamType==0) ptpr=0;
1537 if(fBeamType==1) ptpr=1;
3a72645a 1538
1539 TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
1540 cresidual->Divide(2,1);
1541 cresidual->cd(1);
1542 gPad->SetLogy();
1543 TGraphErrors* residualspectrumD = Normalize(residual);
1544 if(!residualspectrumD) {
1545 AliError("Number of Events not set for the normalization");
1546 return kFALSE;
1547 }
1548 residualspectrumD->SetTitle("");
1549 residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
1550 residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
1551 residualspectrumD->SetMarkerStyle(26);
1552 residualspectrumD->SetMarkerColor(kBlue);
1553 residualspectrumD->SetLineColor(kBlue);
1554 residualspectrumD->Draw("AP");
1555 AliCFDataGrid *dataGridBis = (AliCFDataGrid *) (dataGrid->Clone());
1556 dataGridBis->SetName("dataGridBis");
1557 TGraphErrors* measuredspectrumD = Normalize(dataGridBis);
1558 measuredspectrumD->SetTitle("");
1559 measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
1560 measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
1561 measuredspectrumD->SetMarkerStyle(25);
1562 measuredspectrumD->SetMarkerColor(kBlack);
1563 measuredspectrumD->SetLineColor(kBlack);
1564 measuredspectrumD->Draw("P");
1565 TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
1566 legres->AddEntry(residualspectrumD,"Residual","p");
1567 legres->AddEntry(measuredspectrumD,"Measured","p");
1568 legres->Draw("same");
1569 cresidual->cd(2);
c2690925 1570 TH1D *residualTH1D = residual->Projection(ptpr);
1571 TH1D *measuredTH1D = (TH1D *) dataGridBis->Project(ptpr);
3a72645a 1572 TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
1573 ratioresidual->SetName("ratioresidual");
1574 ratioresidual->SetTitle("");
1575 ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
1576 ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1577 ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
1578 ratioresidual->SetStats(0);
1579 ratioresidual->Draw();
1580
1581 }
1582
c04c80e6 1583 return listofresults;
1584
1585}
1586
1587//____________________________________________________________
3a72645a 1588AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
c04c80e6 1589
1590 //
1591 // Apply unfolding and efficiency correction together to bgsubspectrum
1592 //
1593
3a72645a 1594 AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
c04c80e6 1595 if(!mcContainer){
1596 AliError("MC Container not available");
1597 return NULL;
1598 }
1599
c04c80e6 1600 // Efficiency
3a72645a 1601 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
c04c80e6 1602 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1603
1604 // Data in the right format
1605 AliCFDataGrid *dataGrid = 0x0;
1606 if(bgsubpectrum) {
c04c80e6 1607 dataGrid = bgsubpectrum;
c04c80e6 1608 }
1609 else {
3a72645a 1610
c04c80e6 1611 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1612 if(!dataContainer){
1613 AliError("Data Container not available");
1614 return NULL;
1615 }
1616
3a72645a 1617 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
c04c80e6 1618 }
1619
1620 // Correct
1621 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1622 result->ApplyEffCorrection(*efficiencyD);
c04c80e6 1623
3a72645a 1624 if(fDebugLevel > 0) {
c2690925 1625
1626 Int_t ptpr;
1627 if(fBeamType==0) ptpr=0;
1628 if(fBeamType==1) ptpr=1;
3a72645a 1629
bf892a6a 1630 printf("Step MC: %d\n",fStepMC);
1631 printf("Step tracking: %d\n",AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
1632 printf("Step MC true: %d\n",fStepTrue);
3a72645a 1633 AliCFEffGrid *efficiencymcPID = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
1634 AliCFEffGrid *efficiencymctrackinggeo = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack,fStepTrue);
1635 AliCFEffGrid *efficiencymcall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,fStepTrue);
1636
1637 AliCFEffGrid *efficiencyesdall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerESD),fStepMC,fStepTrue);
1638
1639 TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
1640 cefficiency->cd(1);
c2690925 1641 TH1D* efficiencymcPIDD = (TH1D *) efficiencymcPID->Project(ptpr);
3a72645a 1642 efficiencymcPIDD->SetTitle("");
1643 efficiencymcPIDD->SetStats(0);
1644 efficiencymcPIDD->SetMarkerStyle(25);
1645 efficiencymcPIDD->Draw();
c2690925 1646 TH1D* efficiencymctrackinggeoD = (TH1D *) efficiencymctrackinggeo->Project(ptpr);
3a72645a 1647 efficiencymctrackinggeoD->SetTitle("");
1648 efficiencymctrackinggeoD->SetStats(0);
1649 efficiencymctrackinggeoD->SetMarkerStyle(26);
1650 efficiencymctrackinggeoD->Draw("same");
c2690925 1651 TH1D* efficiencymcallD = (TH1D *) efficiencymcall->Project(ptpr);
3a72645a 1652 efficiencymcallD->SetTitle("");
1653 efficiencymcallD->SetStats(0);
1654 efficiencymcallD->SetMarkerStyle(27);
1655 efficiencymcallD->Draw("same");
c2690925 1656 TH1D* efficiencyesdallD = (TH1D *) efficiencyesdall->Project(ptpr);
3a72645a 1657 efficiencyesdallD->SetTitle("");
1658 efficiencyesdallD->SetStats(0);
1659 efficiencyesdallD->SetMarkerStyle(24);
1660 efficiencyesdallD->Draw("same");
1661 TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
1662 legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
1663 legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
1664 legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
1665 legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
1666 legeff->Draw("same");
c2690925 1667
1668 if(fBeamType==1) {
1669
1670 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1671 TAxis *cenaxisa = sparseafter->GetAxis(0);
1672 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1673 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1674 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1675 TAxis *cenaxisc = efficiencya->GetAxis(0);
1676 //Int_t nbbin = cenaxisb->GetNbins();
1677 //Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1678 //Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1679 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
1680 TString titlee("Efficiency_centrality_bin_");
1681 titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
1682 titlee += "_";
1683 titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
1684 TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1685 cefficiencye->Divide(2,1);
1686 cefficiencye->cd(1);
1687 gPad->SetLogy();
1688 cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
1689 cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
1690 TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr);
1691 TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr);
1692 CorrectFromTheWidth(afterefficiency);
1693 CorrectFromTheWidth(beforeefficiency);
1694 afterefficiency->SetStats(0);
1695 afterefficiency->SetTitle((const char*)titlee);
1696 afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1697 afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1698 afterefficiency->SetMarkerStyle(25);
1699 afterefficiency->SetMarkerColor(kBlack);
1700 afterefficiency->SetLineColor(kBlack);
1701 beforeefficiency->SetStats(0);
1702 beforeefficiency->SetTitle((const char*)titlee);
1703 beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1704 beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1705 beforeefficiency->SetMarkerStyle(24);
1706 beforeefficiency->SetMarkerColor(kBlue);
1707 beforeefficiency->SetLineColor(kBlue);
1708 afterefficiency->Draw();
1709 beforeefficiency->Draw("same");
1710 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1711 lega->AddEntry(beforeefficiency,"Before efficiency correction","p");
1712 lega->AddEntry(afterefficiency,"After efficiency correction","p");
1713 lega->Draw("same");
1714 cefficiencye->cd(2);
1715 cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1);
1716 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr);
1717 efficiencyDDproj->SetTitle("");
1718 efficiencyDDproj->SetStats(0);
1719 efficiencyDDproj->SetMarkerStyle(25);
1720 efficiencyDDproj->SetMarkerColor(2);
1721 efficiencyDDproj->Draw();
1722 cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]);
1723 TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr);
1724 efficiencyDDproja->SetTitle("");
1725 efficiencyDDproja->SetStats(0);
1726 efficiencyDDproja->SetMarkerStyle(26);
1727 efficiencyDDproja->SetMarkerColor(4);
1728 efficiencyDDproja->Draw("same");
1729
1730 }
1731 }
1732
1733
3a72645a 1734 }
1735
c04c80e6 1736 return result;
1737
1738}
3a72645a 1739
c04c80e6 1740//__________________________________________________________________________________
c2690925 1741TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum,Int_t i) const {
c04c80e6 1742 //
1743 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1744 // Give the final pt spectrum to be compared
1745 //
1746
c2690925 1747 if(fNEvents[i] > 0) {
1748
1749 Int_t ptpr;
1750 if(fBeamType==0) ptpr=0;
1751 if(fBeamType==1) ptpr=1;
c04c80e6 1752
c2690925 1753 TH1D* projection = spectrum->Projection(ptpr);
c04c80e6 1754 CorrectFromTheWidth(projection);
c2690925 1755 TGraphErrors *graphError = NormalizeTH1(projection,i);
c04c80e6 1756 return graphError;
1757
1758 }
1759
1760 return 0x0;
1761
1762
1763}
1764//__________________________________________________________________________________
c2690925 1765TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i) const {
c04c80e6 1766 //
1767 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1768 // Give the final pt spectrum to be compared
1769 //
c2690925 1770 if(fNEvents[i] > 0) {
c04c80e6 1771
c2690925 1772 Int_t ptpr;
1773 if(fBeamType==0) ptpr=0;
1774 if(fBeamType==1) ptpr=1;
1775
1776 TH1D* projection = (TH1D *) spectrum->Project(ptpr);
c04c80e6 1777 CorrectFromTheWidth(projection);
c2690925 1778 TGraphErrors *graphError = NormalizeTH1(projection,i);
c04c80e6 1779 return graphError;
1780
1781 }
1782
1783 return 0x0;
1784
1785
1786}
1787//__________________________________________________________________________________
c2690925 1788TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
1789 //
1790 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1791 // Give the final pt spectrum to be compared
1792 //
1793 if(fNEvents[i] > 0) {
1794
1795 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
1796 Double_t p = 0, dp = 0; Int_t point = 1;
1797 Double_t n = 0, dN = 0;
1798 Double_t nCorr = 0, dNcorr = 0;
1799 Double_t errdN = 0, errdp = 0;
1800 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
1801 point = ibin - input->GetXaxis()->GetFirst();
1802 p = input->GetXaxis()->GetBinCenter(ibin);
1803 dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
1804 n = input->GetBinContent(ibin);
1805 dN = input->GetBinError(ibin);
1806
1807 // New point
e156c3bb 1808 nCorr = 0.5 * 1./1.6 * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
c2690925 1809 errdN = 1./(2. * TMath::Pi() * p);
1810 errdp = 1./(2. * TMath::Pi() * p*p) * n;
e156c3bb 1811 dNcorr = 0.5 * 1./1.6 * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
c2690925 1812
1813 spectrumNormalized->SetPoint(point, p, nCorr);
1814 spectrumNormalized->SetPointError(point, dp, dNcorr);
1815 }
1816 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1817 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
1818 spectrumNormalized->SetMarkerStyle(22);
1819 spectrumNormalized->SetMarkerColor(kBlue);
1820 spectrumNormalized->SetLineColor(kBlue);
1821
1822 return spectrumNormalized;
1823
1824 }
1825
1826 return 0x0;
1827
1828
1829}
1830//__________________________________________________________________________________
1831TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) const {
c04c80e6 1832 //
1833 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1834 // Give the final pt spectrum to be compared
1835 //
c2690925 1836 if(normalization > 0) {
c04c80e6 1837
1838 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
1839 Double_t p = 0, dp = 0; Int_t point = 1;
1840 Double_t n = 0, dN = 0;
1841 Double_t nCorr = 0, dNcorr = 0;
1842 Double_t errdN = 0, errdp = 0;
1843 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
1844 point = ibin - input->GetXaxis()->GetFirst();
1845 p = input->GetXaxis()->GetBinCenter(ibin);
1846 dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
1847 n = input->GetBinContent(ibin);
1848 dN = input->GetBinError(ibin);
1849
1850 // New point
e156c3bb 1851 nCorr = 0.5 * 1./1.6 * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
c04c80e6 1852 errdN = 1./(2. * TMath::Pi() * p);
1853 errdp = 1./(2. * TMath::Pi() * p*p) * n;
e156c3bb 1854 dNcorr = 0.5 * 1./1.6 * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
c04c80e6 1855
1856 spectrumNormalized->SetPoint(point, p, nCorr);
1857 spectrumNormalized->SetPointError(point, dp, dNcorr);
1858 }
1859 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1860 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
1861 spectrumNormalized->SetMarkerStyle(22);
1862 spectrumNormalized->SetMarkerColor(kBlue);
1863 spectrumNormalized->SetLineColor(kBlue);
1864
1865 return spectrumNormalized;
1866
1867 }
1868
1869 return 0x0;
1870
1871
1872}
1873//____________________________________________________________
1874void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
1875 //
1876 // Set the container for a given step to the
1877 //
1878 if(!fCFContainers) fCFContainers = new TList;
1879 fCFContainers->AddAt(cont, type);
1880}
1881
1882//____________________________________________________________
1883AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
1884 //
1885 // Get Correction Framework Container for given type
1886 //
1887 if(!fCFContainers) return NULL;
1888 return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
1889}
c04c80e6 1890//____________________________________________________________
c2690925 1891AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Int_t positivenegative) {
c04c80e6 1892 //
3a72645a 1893 // Slice bin for a given source of electron
c2690925 1894 // nDim is the number of dimension the corrections are done
1895 // dimensions are the definition of the dimensions
1896 // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
1897 // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
c04c80e6 1898 //
1899
1900 Double_t *varMin = new Double_t[container->GetNVar()],
1901 *varMax = new Double_t[container->GetNVar()];
1902
1903 Double_t *binLimits;
1904 for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
1905
1906 binLimits = new Double_t[container->GetNBins(ivar)+1];
1907 container->GetBinLimits(ivar,binLimits);
c2690925 1908 varMin[ivar] = binLimits[0];
1909 varMax[ivar] = binLimits[container->GetNBins(ivar)];
1910 // source
1911 if(ivar == 4){
1912 if((source>= 0) && (source<container->GetNBins(ivar))) {
1913 varMin[ivar] = binLimits[source];
1914 varMax[ivar] = binLimits[source];
1915 }
3a72645a 1916 }
c2690925 1917 // charge
1918 if(ivar == 3) {
1919 if((positivenegative>= 0) && (positivenegative<container->GetNBins(ivar))) {
1920 varMin[ivar] = binLimits[positivenegative];
1921 varMax[ivar] = binLimits[positivenegative];
1922 }
3a72645a 1923 }
c2690925 1924
3a72645a 1925
c04c80e6 1926 delete[] binLimits;
1927 }
1928
1929 AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
1930 AddTemporaryObject(k);
1931 delete[] varMin; delete[] varMax;
1932
1933 return k;
1934
1935}
1936
1937//_________________________________________________________________________
3a72645a 1938THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions) const {
c04c80e6 1939 //
3a72645a 1940 // Slice correlation
c04c80e6 1941 //
1942
3a72645a 1943 Int_t ndimensions = correlationmatrix->GetNdimensions();
c2690925 1944 //printf("Number of dimension %d correlation map\n",ndimensions);
c04c80e6 1945 if(ndimensions < (2*nDim)) {
1946 AliError("Problem in the dimensions");
1947 return NULL;
1948 }
1949 Int_t ndimensionsContainer = (Int_t) ndimensions/2;
c2690925 1950 //printf("Number of dimension %d container\n",ndimensionsContainer);
c04c80e6 1951
1952 Int_t *dim = new Int_t[nDim*2];
1953 for(Int_t iter=0; iter < nDim; iter++){
1954 dim[iter] = dimensions[iter];
1955 dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
c2690925 1956 //printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
c04c80e6 1957 }
1958
3a72645a 1959 THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
c04c80e6 1960
1961 delete[] dim;
1962 return k;
1963
1964}
1965//___________________________________________________________________________
1966void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
1967 //
1968 // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
1969 //
1970
1971 TAxis *axis = h1->GetXaxis();
1972 Int_t nbinX = h1->GetNbinsX();
1973
1974 for(Int_t i = 1; i <= nbinX; i++) {
1975
1976 Double_t width = axis->GetBinWidth(i);
1977 Double_t content = h1->GetBinContent(i);
1978 Double_t error = h1->GetBinError(i);
1979 h1->SetBinContent(i,content/width);
1980 h1->SetBinError(i,error/width);
1981 }
1982
1983}
1984//____________________________________________________________
1985void AliHFEspectrum::AddTemporaryObject(TObject *o){
1986 //
1987 // Emulate garbage collection on class level
1988 //
1989 if(!fTemporaryObjects) {
1990 fTemporaryObjects= new TList;
1991 fTemporaryObjects->SetOwner();
1992 }
1993 fTemporaryObjects->Add(o);
1994}
1995
1996//____________________________________________________________
1997void AliHFEspectrum::ClearObject(TObject *o){
1998 //
1999 // Do a safe deletion
2000 //
2001 if(fTemporaryObjects){
2002 if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
2003 delete o;
2004 } else{
2005 // just delete
2006 delete o;
2007 }
2008}
2009//_________________________________________________________________________
c2690925 2010TObject* AliHFEspectrum::GetSpectrum(const AliCFContainer * const c, Int_t step) {
245877d0 2011 AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
c04c80e6 2012 return data;
2013}
2014//_________________________________________________________________________
c2690925 2015TObject* AliHFEspectrum::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
c04c80e6 2016 //
2017 // Create efficiency grid and calculate efficiency
2018 // of step to step0
2019 //
2020 TString name("eff");
2021 name += step;
2022 name+= step0;
2023 AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
2024 eff->CalculateEfficiency(step,step0);
2025 return eff;
2026}
c2690925 2027
2028//____________________________________________________________________________
2029THnSparse* AliHFEspectrum::GetWeights(){
2030
2031 //
2032 // Measured D->e based weighting factors
2033 //
2034
2035 const Int_t nDim=1;
2036 Int_t nBin[nDim] = {44};
2037 const Double_t kPtbound[2] = {0.1, 20.};
2038
2039 Double_t* binEdges[nDim];
2040 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
2041
2042 fWeightCharm = new THnSparseF("weightHisto", "weiting factor; pt[GeV/c]", nDim, nBin);
2043 for(Int_t idim = 0; idim < nDim; idim++)
2044 fWeightCharm->SetBinEdges(idim, binEdges[idim]);
2045 /* //fit
2046 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};*/
2047 //points
2048 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};
2049 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};
2050 Double_t dataE[1];
2051 for(int i=0; i<44; i++){
2052 dataE[0]=pt[i];
2053 fWeightCharm->Fill(dataE,weight[i]);
2054 }
2055 Int_t* ibins[nDim];
2056 for(Int_t ivar = 0; ivar < nDim; ivar++)
2057 ibins[ivar] = new Int_t[nBin[ivar] + 1];
2058 fWeightCharm->SetBinError(ibins[0],0);
2059
2060 return fWeightCharm;
2061}