]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliHFEspectrum.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updated treatment of D0/D0bar mass assumption (Carlos)
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEspectrum.cxx
CommitLineData
a8ef1999 1
c04c80e6 2/**************************************************************************
3* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4* *
5* Author: The ALICE Off-line Project. *
6* Contributors are mentioned in the code where appropriate. *
7* *
8* Permission to use, copy, modify and distribute this software and its *
9* documentation strictly for non-commercial purposes is hereby granted *
10* without fee, provided that the above copyright notice appears in all *
11* copies and that both the copyright notice and this permission notice *
12* appear in the supporting documentation. The authors make no claims *
13* about the suitability of this software for any purpose. It is *
14* provided "as is" without express or implied warranty. *
15**************************************************************************/
16//
17// Class for spectrum correction
18// Subtraction of hadronic background, Unfolding of the data and
19// Renormalization done here
20// The following containers have to be set:
21// - Correction framework container for real data
22// - Correction framework container for MC (Efficiency Map)
23// - Correction framework container for background coming from data
24// - Correction framework container for background coming from MC
25//
245877d0 26// Author:
27// Raphaelle Bailhache <R.Bailhache@gsi.de>
28// Markus Fasel <M.Fasel@gsi.de>
c04c80e6 29//
30
31#include <TArrayD.h>
32#include <TH1.h>
33#include <TList.h>
34#include <TObjArray.h>
35#include <TROOT.h>
36#include <TCanvas.h>
37#include <TLegend.h>
38#include <TStyle.h>
39#include <TMath.h>
40#include <TAxis.h>
41#include <TGraphErrors.h>
42#include <TFile.h>
43#include <TPad.h>
67fe7bd0 44#include <TH2D.h>
c2690925 45#include <TF1.h>
c04c80e6 46
3a72645a 47#include "AliPID.h"
c04c80e6 48#include "AliCFContainer.h"
49#include "AliCFDataGrid.h"
50#include "AliCFEffGrid.h"
51#include "AliCFGridSparse.h"
52#include "AliCFUnfolding.h"
53#include "AliLog.h"
54
55#include "AliHFEspectrum.h"
3a72645a 56#include "AliHFEcuts.h"
57#include "AliHFEcontainer.h"
c2690925 58#include "AliHFEtools.h"
c04c80e6 59
60ClassImp(AliHFEspectrum)
61
62//____________________________________________________________
63AliHFEspectrum::AliHFEspectrum(const char *name):
64 TNamed(name, ""),
e17c1f86 65 fCFContainers(new TObjArray(kDataContainerV0+1)),
c04c80e6 66 fTemporaryObjects(NULL),
67 fCorrelation(NULL),
68 fBackground(NULL),
c2690925 69 fEfficiencyFunction(NULL),
70 fWeightCharm(NULL),
3a72645a 71 fInclusiveSpectrum(kTRUE),
c04c80e6 72 fDumpToFile(kFALSE),
8c1c76e9 73 fEtaSelected(kFALSE),
c2690925 74 fUnSetCorrelatedErrors(kTRUE),
e156c3bb 75 fSetSmoothing(kFALSE),
c2690925 76 fIPanaHadronBgSubtract(kFALSE),
77 fIPanaCharmBgSubtract(kFALSE),
78 fIPanaConversionBgSubtract(kFALSE),
79 fIPanaNonHFEBgSubtract(kFALSE),
a8ef1999 80 fIPParameterizedEff(kFALSE),
e17c1f86 81 fNonHFEsyst(0),
a8ef1999 82 fBeauty2ndMethod(kFALSE),
11ff28c5 83 fIPEffCombinedSamples(kTRUE),
3a72645a 84 fNbDimensions(1),
c2690925 85 fNMCEvents(0),
c04c80e6 86 fStepMC(-1),
87 fStepTrue(-1),
88 fStepData(-1),
3a72645a 89 fStepBeforeCutsV0(-1),
90 fStepAfterCutsV0(-1),
c04c80e6 91 fStepGuessedUnfolding(-1),
3a72645a 92 fNumberOfIterations(5),
e17c1f86 93 fChargeChoosen(kAllCharge),
c2690925 94 fNCentralityBinAtTheEnd(0),
95 fHadronEffbyIPcut(NULL),
11ff28c5 96 fConversionEffbgc(NULL),
97 fNonHFEEffbgc(NULL),
a8ef1999 98 fBSpectrum2ndMethod(NULL),
99 fkBeauty2ndMethodfilename(""),
c2690925 100 fBeamType(0),
e17c1f86 101 fDebugLevel(0),
102 fWriteToFile(kFALSE)
c04c80e6 103{
104 //
105 // Default constructor
106 //
c2690925 107
108 for(Int_t k = 0; k < 20; k++){
a8ef1999 109 fNEvents[k] = 0;
110 fNMCbgEvents[k] = 0;
111 fLowBoundaryCentralityBinAtTheEnd[k] = 0;
112 fHighBoundaryCentralityBinAtTheEnd[k] = 0;
113 if(k<kCentrality)
114 {
115 fEfficiencyTOFPIDD[k] = 0;
11ff28c5 116 fEfficiencyesdTOFPIDD[k] = 0;
a8ef1999 117 fEfficiencyIPCharmD[k] = 0;
118 fEfficiencyIPBeautyD[k] = 0;
119 fEfficiencyIPConversionD[k] = 0;
120 fEfficiencyIPNonhfeD[k] = 0;
121
122 fConversionEff[k] = 0;
123 fNonHFEEff[k] = 0;
124 fCharmEff[k] = 0;
125 fBeautyEff[k] = 0;
126 fEfficiencyCharmSigD[k] = 0;
127 fEfficiencyBeautySigD[k] = 0;
128 }
c2690925 129 }
8c1c76e9 130 memset(fEtaRange, 0, sizeof(Double_t) * 2);
e17c1f86 131 memset(fEtaRangeNorm, 0, sizeof(Double_t) * 2);
132 memset(fConvSourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
133 memset(fNonHFESourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
c04c80e6 134}
135
136//____________________________________________________________
137AliHFEspectrum::~AliHFEspectrum(){
138 //
139 // Destructor
140 //
141 if(fCFContainers) delete fCFContainers;
142 if(fTemporaryObjects){
143 fTemporaryObjects->Clear();
144 delete fTemporaryObjects;
145 }
146}
3a72645a 147//____________________________________________________________
c2690925 148Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *v0hfecontainer, const AliHFEcontainer *bghfecontainer){
3a72645a 149 //
150 // Init what we need for the correction:
151 //
152 // Raw spectrum, hadron contamination
153 // MC efficiency maps, correlation matrix
154 // V0 efficiency if wanted
155 //
156 // This for a given dimension.
157 // If no inclusive spectrum, then take only efficiency map for beauty electron
158 // and the appropriate correlation matrix
159 //
a8ef1999 160
161
162 if(fBeauty2ndMethod) CallInputFileForBeauty2ndMethod();
163
e17c1f86 164 Int_t kNdim = 3;
a8ef1999 165 Int_t kNcentr =1;
166 Int_t ptpr =0;
e17c1f86 167 if(fBeamType==0) kNdim=3;
a8ef1999 168 if(fBeamType==1)
169 {
170 kNdim=4;
171 kNcentr=11;
172 ptpr=1;
173 }
e17c1f86 174
175 Int_t dims[kNdim];
176 // Get the requested format
177 if(fBeamType==0){
178 // pp analysis
179 switch(fNbDimensions){
180 case 1: dims[0] = 0;
181 break;
182 case 2: for(Int_t i = 0; i < 2; i++) dims[i] = i;
183 break;
184 case 3: for(Int_t i = 0; i < 3; i++) dims[i] = i;
185 break;
186 default:
187 AliError("Container with this number of dimensions not foreseen (yet)");
188 return kFALSE;
189 };
190 }
191
192 if(fBeamType==1){
193 // PbPb analysis; centrality as first dimension
194 Int_t nbDimensions = fNbDimensions;
195 fNbDimensions = fNbDimensions + 1;
196 switch(nbDimensions){
197 case 1: dims[0] = 5;
198 dims[1] = 0;
199 break;
200 case 2: dims[0] = 5;
201 for(Int_t i = 0; i < 2; i++) dims[(i+1)] = i;
202 break;
203 case 3: dims[0] = 5;
204 for(Int_t i = 0; i < 3; i++) dims[(i+1)] = i;
205 break;
206 default:
207 AliError("Container with this number of dimensions not foreseen (yet)");
208 return kFALSE;
209 };
210 }
3a72645a 211
3a72645a 212 // Data container: raw spectrum + hadron contamination
c2690925 213 AliCFContainer *datacontainer = 0x0;
214 if(fInclusiveSpectrum) {
215 datacontainer = datahfecontainer->GetCFContainer("recTrackContReco");
216 }
217 else{
a8ef1999 218
219 datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
c2690925 220 }
3a72645a 221 AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
222 if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
223
c2690925 224 AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen);
225 AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen);
3a72645a 226 if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
c2690925 227
3a72645a 228 SetContainer(datacontainerD,AliHFEspectrum::kDataContainer);
229 SetContainer(contaminationcontainerD,AliHFEspectrum::kBackgroundData);
230
231 // MC container: ESD/MC efficiency maps + MC/MC efficiency maps
232 AliCFContainer *mccontaineresd = 0x0;
11ff28c5 233 AliCFContainer *mccontaineresdbg = 0x0;
3a72645a 234 AliCFContainer *mccontainermc = 0x0;
8c1c76e9 235 AliCFContainer *mccontainermcbg = 0x0;
236 AliCFContainer *nonHFEweightedContainer = 0x0;
237 AliCFContainer *convweightedContainer = 0x0;
e17c1f86 238 AliCFContainer *nonHFEtempContainer = 0x0;//temporary container to be sliced for the fnonHFESourceContainers
239 AliCFContainer *convtempContainer = 0x0;//temporary container to be sliced for the fConvSourceContainers
8c1c76e9 240
3a72645a 241 if(fInclusiveSpectrum) {
242 mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco");
243 mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC");
244 }
245 else {
246 mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
11ff28c5 247 mccontaineresdbg = bghfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
3a72645a 248 mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
8c1c76e9 249 mccontainermcbg = bghfecontainer->MakeMergedCFContainer("summcbg","summcbg","MCTrackCont:recTrackContMC:recTrackContDEMC");
e17c1f86 250
251 if(fNonHFEsyst){
252 const Char_t *sourceName[kElecBgSources]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
253 const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
254 for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
255 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
256 nonHFEtempContainer = bghfecontainer->GetCFContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]));
a8ef1999 257 convtempContainer = bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
258 for(Int_t icentr=0;icentr<kNcentr;icentr++)
259 {
260 if(fBeamType==0)
261 {
262 fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
263 fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
264 }
265 if(fBeamType==1)
266 {
11ff28c5 267
268 fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr+1);
269 fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr+1);
a8ef1999 270 }
11ff28c5 271// if((!fConvSourceContainer[iSource][iLevel][icentr])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;
a8ef1999 272 }
e17c1f86 273 }
274 }
275 }
276 else{
277 nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
278 convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
279 if((!convweightedContainer)||(!nonHFEweightedContainer)) return kFALSE;
280 }
3a72645a 281 }
282 if((!mccontaineresd) || (!mccontainermc)) return kFALSE;
e17c1f86 283
3a72645a 284 Int_t source = -1;
8c1c76e9 285 if(!fInclusiveSpectrum) source = 1; //beauty
c2690925 286 AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen);
287 AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen);
3a72645a 288 if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
289 SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerMC);
290 SetContainer(mccontaineresdD,AliHFEspectrum::kMCContainerESD);
291
c2690925 292 // set charm, nonHFE container to estimate BG
293 if(!fInclusiveSpectrum){
294 source = 0; //charm
8c1c76e9 295 mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
c2690925 296 SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerCharmMC);
8c1c76e9 297
11ff28c5 298 SetParameterizedEff(mccontainermc, mccontainermcbg, mccontaineresd, mccontaineresdbg, dims);
8c1c76e9 299
e17c1f86 300 if(!fNonHFEsyst){
301 AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
302 SetContainer(nonHFEweightedContainerD,AliHFEspectrum::kMCWeightedContainerNonHFEESD);
303 AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
304 SetContainer(convweightedContainerD,AliHFEspectrum::kMCWeightedContainerConversionESD);
11ff28c5 305
306 if(fBeamType==0){
307 AliCFEffGrid *nonHFEEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
308 fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
309
310 AliCFEffGrid *conversionEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
311 fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
312 }
313
e17c1f86 314 }
c2690925 315 }
e17c1f86 316 // MC container: correlation matrix
3a72645a 317 THnSparseF *mccorrelation = 0x0;
bf892a6a 318 if(fInclusiveSpectrum) {
3ccf8f4c 319 if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
320 else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
321 else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
322 else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD - 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
6c70d827 323 else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
e156c3bb 324
325 if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
bf892a6a 326 }
c2690925 327 else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
328 //else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
3a72645a 329 if(!mccorrelation) return kFALSE;
330 THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims);
331 if(!mccorrelationD) {
332 printf("No correlation\n");
333 return kFALSE;
334 }
335 SetCorrelation(mccorrelationD);
336
337 // V0 container Electron, pt eta phi
338 if(v0hfecontainer) {
339 AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
340 if(!containerV0) return kFALSE;
e17c1f86 341 AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron);
3a72645a 342 if(!containerV0Electron) return kFALSE;
343 SetContainer(containerV0Electron,AliHFEspectrum::kDataContainerV0);
8c1c76e9 344 }
345
3a72645a 346
347 if(fDebugLevel>0){
348
349 AliCFDataGrid *contaminationspectrum = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(contaminationcontainer,1))->Clone();
350 contaminationspectrum->SetName("contaminationspectrum");
351 TCanvas * ccontaminationspectrum = new TCanvas("contaminationspectrum","contaminationspectrum",1000,700);
352 ccontaminationspectrum->Divide(3,1);
353 ccontaminationspectrum->cd(1);
6555e2ad 354 TH2D * contaminationspectrum2dpteta = (TH2D *) contaminationspectrum->Project(1,0);
355 TH2D * contaminationspectrum2dptphi = (TH2D *) contaminationspectrum->Project(2,0);
356 TH2D * contaminationspectrum2detaphi = (TH2D *) contaminationspectrum->Project(1,2);
3a72645a 357 contaminationspectrum2dpteta->SetStats(0);
358 contaminationspectrum2dpteta->SetTitle("");
359 contaminationspectrum2dpteta->GetXaxis()->SetTitle("#eta");
360 contaminationspectrum2dpteta->GetYaxis()->SetTitle("p_{T} [GeV/c]");
361 contaminationspectrum2dptphi->SetStats(0);
362 contaminationspectrum2dptphi->SetTitle("");
363 contaminationspectrum2dptphi->GetXaxis()->SetTitle("#phi [rad]");
364 contaminationspectrum2dptphi->GetYaxis()->SetTitle("p_{T} [GeV/c]");
365 contaminationspectrum2detaphi->SetStats(0);
366 contaminationspectrum2detaphi->SetTitle("");
367 contaminationspectrum2detaphi->GetXaxis()->SetTitle("#eta");
368 contaminationspectrum2detaphi->GetYaxis()->SetTitle("#phi [rad]");
369 contaminationspectrum2dptphi->Draw("colz");
370 ccontaminationspectrum->cd(2);
371 contaminationspectrum2dpteta->Draw("colz");
372 ccontaminationspectrum->cd(3);
373 contaminationspectrum2detaphi->Draw("colz");
374
c2690925 375 TCanvas * ccorrelation = new TCanvas("ccorrelation","ccorrelation",1000,700);
376 ccorrelation->cd(1);
377 if(mccorrelationD) {
378 if(fBeamType==0){
e17c1f86 379 TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions,0);
380 ptcorrelation->Draw("colz");
c2690925 381 }
382 if(fBeamType==1){
e17c1f86 383 TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions+1,1);
384 ptcorrelation->Draw("colz");
c2690925 385 }
386 }
e17c1f86 387 if(fWriteToFile) ccontaminationspectrum->SaveAs("contaminationspectrum.eps");
3a72645a 388 }
389
390
391 return kTRUE;
3a72645a 392}
c04c80e6 393
c2690925 394
a8ef1999 395//____________________________________________________________
396void AliHFEspectrum::CallInputFileForBeauty2ndMethod(){
397 //
398 // get spectrum for beauty 2nd method
399 //
400 //
401 TFile *inputfile2ndmethod=TFile::Open(fkBeauty2ndMethodfilename);
402 fBSpectrum2ndMethod = new TH1D(*static_cast<TH1D*>(inputfile2ndmethod->Get("BSpectrum")));
403}
404
405
c04c80e6 406//____________________________________________________________
3a72645a 407Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
c04c80e6 408 //
409 // Correct the spectrum for efficiency and unfolding
410 // with both method and compare
411 //
412
413 gStyle->SetPalette(1);
414 gStyle->SetOptStat(1111);
415 gStyle->SetPadBorderMode(0);
416 gStyle->SetCanvasColor(10);
417 gStyle->SetPadLeftMargin(0.13);
418 gStyle->SetPadRightMargin(0.13);
67fe7bd0 419
c2690925 420 printf("Steps are: stepdata %d, stepMC %d, steptrue %d, stepV0after %d, stepV0before %d\n",fStepData,fStepMC,fStepTrue,fStepAfterCutsV0,fStepBeforeCutsV0);
421
422 ///////////////////////////
423 // Check initialization
424 ///////////////////////////
425
426 if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
427 AliInfo("You have to init before");
428 return kFALSE;
429 }
430
a199006c 431 if((fStepTrue < 0) && (fStepMC < 0) && (fStepData < 0)) {
c2690925 432 AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
433 return kFALSE;
434 }
435
436 SetNumberOfIteration(10);
437 SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
438
439 AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
440 //////////////////////////////////
441 // Subtract hadron background
442 /////////////////////////////////
443 AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
444 if(subtractcontamination) {
445 dataspectrumaftersubstraction = SubtractBackground(kTRUE);
446 dataGridAfterFirstSteps = dataspectrumaftersubstraction;
447 }
448
449 ////////////////////////////////////////////////
450 // Correct for TPC efficiency from V0
451 ///////////////////////////////////////////////
452 AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
453 AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
8c1c76e9 454 if(dataContainerV0){printf("Got the V0 container\n");
c2690925 455 dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
456 dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
457 }
458
459 //////////////////////////////////////////////////////////////////////////////
460 // Correct for efficiency parametrized (if TPC efficiency is parametrized)
461 /////////////////////////////////////////////////////////////////////////////
462 AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
463 if(fEfficiencyFunction){
464 dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
465 dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
466 }
467
468 ///////////////
469 // Unfold
470 //////////////
471 TList *listunfolded = Unfold(dataGridAfterFirstSteps);
472 if(!listunfolded){
473 printf("Unfolded failed\n");
474 return kFALSE;
475 }
476 THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
477 THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
478 if(!correctedspectrum){
479 AliError("No corrected spectrum\n");
480 return kFALSE;
481 }
482 if(!residualspectrum){
483 AliError("No residul spectrum\n");
484 return kFALSE;
485 }
486
487 /////////////////////
488 // Simply correct
489 ////////////////////
490 AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
491
492
493 //////////
494 // Plot
495 //////////
496 if(fDebugLevel > 0.0) {
497
a199006c 498 Int_t ptpr = 0;
c2690925 499 if(fBeamType==0) ptpr=0;
500 if(fBeamType==1) ptpr=1;
501
502 TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
503 ccorrected->Divide(2,1);
504 ccorrected->cd(1);
505 gPad->SetLogy();
506 TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
507 correctedspectrumD->SetTitle("");
508 correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
509 correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
510 correctedspectrumD->SetMarkerStyle(26);
511 correctedspectrumD->SetMarkerColor(kBlue);
512 correctedspectrumD->SetLineColor(kBlue);
513 correctedspectrumD->Draw("AP");
514 TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
515 alltogetherspectrumD->SetTitle("");
516 alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
517 alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
518 alltogetherspectrumD->SetMarkerStyle(25);
519 alltogetherspectrumD->SetMarkerColor(kBlack);
520 alltogetherspectrumD->SetLineColor(kBlack);
521 alltogetherspectrumD->Draw("P");
522 TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
523 legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
524 legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
525 legcorrected->Draw("same");
526 ccorrected->cd(2);
527 TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
528 TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
529 TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
530 ratiocorrected->SetName("ratiocorrected");
531 ratiocorrected->SetTitle("");
532 ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
533 ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
534 ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
535 ratiocorrected->SetStats(0);
536 ratiocorrected->Draw();
e17c1f86 537 if(fWriteToFile)ccorrected->SaveAs("CorrectedPbPb.eps");
c2690925 538
a199006c 539 //TH1D unfoldingspectrac[fNCentralityBinAtTheEnd];
540 //TGraphErrors unfoldingspectracn[fNCentralityBinAtTheEnd];
541 //TH1D correctedspectrac[fNCentralityBinAtTheEnd];
542 //TGraphErrors correctedspectracn[fNCentralityBinAtTheEnd];
543
544 TH1D *unfoldingspectrac = new TH1D[fNCentralityBinAtTheEnd];
545 TGraphErrors *unfoldingspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
546 TH1D *correctedspectrac = new TH1D[fNCentralityBinAtTheEnd];
547 TGraphErrors *correctedspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
c2690925 548
c2690925 549 if(fBeamType==1) {
550
e156c3bb 551 TCanvas * ccorrectedallspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700);
552 ccorrectedallspectra->Divide(2,1);
c2690925 553 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
554 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
555
556 THnSparseF* sparsesu = (THnSparseF *) correctedspectrum;
557 TAxis *cenaxisa = sparsesu->GetAxis(0);
558 THnSparseF* sparsed = (THnSparseF *) alltogetherCorrection->GetGrid();
559 TAxis *cenaxisb = sparsed->GetAxis(0);
560 Int_t nbbin = cenaxisb->GetNbins();
561 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
562 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
563 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
e17c1f86 564 TString titlee("corrected_centrality_bin_");
565 titlee += "[";
566 titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
567 titlee += "_";
568 titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
569 titlee += "[";
570 TString titleec("corrected_check_projection_bin_");
571 titleec += "[";
572 titleec += fLowBoundaryCentralityBinAtTheEnd[binc];
573 titleec += "_";
574 titleec += fHighBoundaryCentralityBinAtTheEnd[binc];
575 titleec += "[";
576 TString titleenameunotnormalized("Unfolded_Notnormalized_centrality_bin_");
577 titleenameunotnormalized += "[";
578 titleenameunotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
579 titleenameunotnormalized += "_";
580 titleenameunotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
581 titleenameunotnormalized += "[";
582 TString titleenameunormalized("Unfolded_normalized_centrality_bin_");
583 titleenameunormalized += "[";
584 titleenameunormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
585 titleenameunormalized += "_";
586 titleenameunormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
587 titleenameunormalized += "[";
588 TString titleenamednotnormalized("Dirrectcorrected_Notnormalized_centrality_bin_");
589 titleenamednotnormalized += "[";
590 titleenamednotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
591 titleenamednotnormalized += "_";
592 titleenamednotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
593 titleenamednotnormalized += "[";
594 TString titleenamednormalized("Dirrectedcorrected_normalized_centrality_bin_");
595 titleenamednormalized += "[";
596 titleenamednormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
597 titleenamednormalized += "_";
598 titleenamednormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
599 titleenamednormalized += "[";
600 Int_t nbEvents = 0;
601 for(Int_t k = fLowBoundaryCentralityBinAtTheEnd[binc]; k < fHighBoundaryCentralityBinAtTheEnd[binc]; k++) {
602 printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
603 nbEvents += fNEvents[k];
604 }
605 Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
606 Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
607 printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
608 Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
609 Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
610 printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
611 cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
612 cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
613 TCanvas * ccorrectedcheck = new TCanvas((const char*) titleec,(const char*) titleec,1000,700);
614 ccorrectedcheck->cd(1);
615 TH1D *aftersuc = (TH1D *) sparsesu->Projection(0);
616 TH1D *aftersdc = (TH1D *) sparsed->Projection(0);
617 aftersuc->Draw();
618 aftersdc->Draw("same");
619 TCanvas * ccorrectede = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
620 ccorrectede->Divide(2,1);
621 ccorrectede->cd(1);
622 gPad->SetLogy();
623 TH1D *aftersu = (TH1D *) sparsesu->Projection(1);
624 CorrectFromTheWidth(aftersu);
625 aftersu->SetName((const char*)titleenameunotnormalized);
626 unfoldingspectrac[binc] = *aftersu;
627 ccorrectede->cd(1);
628 TGraphErrors* aftersun = NormalizeTH1N(aftersu,nbEvents);
629 aftersun->SetTitle("");
630 aftersun->GetYaxis()->SetTitleOffset(1.5);
631 aftersun->GetYaxis()->SetRangeUser(0.000000001,1.0);
632 aftersun->SetMarkerStyle(26);
633 aftersun->SetMarkerColor(kBlue);
634 aftersun->SetLineColor(kBlue);
635 aftersun->Draw("AP");
636 aftersun->SetName((const char*)titleenameunormalized);
637 unfoldingspectracn[binc] = *aftersun;
638 ccorrectede->cd(1);
639 TH1D *aftersd = (TH1D *) sparsed->Projection(1);
640 CorrectFromTheWidth(aftersd);
641 aftersd->SetName((const char*)titleenamednotnormalized);
642 correctedspectrac[binc] = *aftersd;
643 ccorrectede->cd(1);
644 TGraphErrors* aftersdn = NormalizeTH1N(aftersd,nbEvents);
645 aftersdn->SetTitle("");
646 aftersdn->GetYaxis()->SetTitleOffset(1.5);
647 aftersdn->GetYaxis()->SetRangeUser(0.000000001,1.0);
648 aftersdn->SetMarkerStyle(25);
649 aftersdn->SetMarkerColor(kBlack);
650 aftersdn->SetLineColor(kBlack);
651 aftersdn->Draw("P");
652 aftersdn->SetName((const char*)titleenamednormalized);
653 correctedspectracn[binc] = *aftersdn;
654 TLegend *legcorrectedud = new TLegend(0.4,0.6,0.89,0.89);
655 legcorrectedud->AddEntry(aftersun,"Corrected","p");
656 legcorrectedud->AddEntry(aftersdn,"Alltogether","p");
657 legcorrectedud->Draw("same");
658 ccorrectedallspectra->cd(1);
659 gPad->SetLogy();
660 TH1D *aftersunn = (TH1D *) aftersun->Clone();
661 aftersunn->SetMarkerStyle(stylee[binc]);
662 aftersunn->SetMarkerColor(colorr[binc]);
663 if(binc==0) aftersunn->Draw("AP");
664 else aftersunn->Draw("P");
665 legtotal->AddEntry(aftersunn,(const char*) titlee,"p");
666 ccorrectedallspectra->cd(2);
667 gPad->SetLogy();
668 TH1D *aftersdnn = (TH1D *) aftersdn->Clone();
669 aftersdnn->SetMarkerStyle(stylee[binc]);
670 aftersdnn->SetMarkerColor(colorr[binc]);
671 if(binc==0) aftersdnn->Draw("AP");
672 else aftersdnn->Draw("P");
673 legtotalg->AddEntry(aftersdnn,(const char*) titlee,"p");
674 ccorrectede->cd(2);
675 TH1D* ratiocorrectedbinc = (TH1D*)aftersu->Clone();
676 TString titleee("ratiocorrected_bin_");
677 titleee += binc;
678 ratiocorrectedbinc->SetName((const char*) titleee);
679 ratiocorrectedbinc->SetTitle("");
680 ratiocorrectedbinc->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
681 ratiocorrectedbinc->GetXaxis()->SetTitle("p_{T} [GeV/c]");
682 ratiocorrectedbinc->Divide(aftersu,aftersd,1,1);
683 ratiocorrectedbinc->SetStats(0);
684 ratiocorrectedbinc->Draw();
c2690925 685 }
686
e156c3bb 687 ccorrectedallspectra->cd(1);
c2690925 688 legtotal->Draw("same");
e156c3bb 689 ccorrectedallspectra->cd(2);
c2690925 690 legtotalg->Draw("same");
691
692 cenaxisa->SetRange(0,nbbin);
693 cenaxisb->SetRange(0,nbbin);
e17c1f86 694 if(fWriteToFile) ccorrectedallspectra->SaveAs("CorrectedPbPb.eps");
c2690925 695 }
696
697 // Dump to file if needed
698 if(fDumpToFile) {
699 TFile *out = new TFile("finalSpectrum.root","recreate");
700 correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
701 correctedspectrumD->Write();
702 alltogetherspectrumD->SetName("AlltogetherSpectrum");
703 alltogetherspectrumD->Write();
704 ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
705 ratiocorrected->Write();
706 correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
707 correctedspectrum->Write();
708 alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
709 alltogetherCorrection->Write();
710 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
e17c1f86 711 unfoldingspectrac[binc].Write();
712 unfoldingspectracn[binc].Write();
713 correctedspectrac[binc].Write();
714 correctedspectracn[binc].Write();
c2690925 715 }
716 out->Close(); delete out;
717 }
718
a199006c 719 if (unfoldingspectrac) delete[] unfoldingspectrac;
720 if (unfoldingspectracn) delete[] unfoldingspectracn;
721 if (correctedspectrac) delete[] correctedspectrac;
722 if (correctedspectracn) delete[] correctedspectracn;
723
017dcb19 724 }
c2690925 725
726 return kTRUE;
727}
728
729//____________________________________________________________
730Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
731 //
732 // Correct the spectrum for efficiency and unfolding for beauty analysis
733 // with both method and compare
734 //
735
736 gStyle->SetPalette(1);
737 gStyle->SetOptStat(1111);
738 gStyle->SetPadBorderMode(0);
739 gStyle->SetCanvasColor(10);
740 gStyle->SetPadLeftMargin(0.13);
741 gStyle->SetPadRightMargin(0.13);
742
3a72645a 743 ///////////////////////////
744 // Check initialization
745 ///////////////////////////
c04c80e6 746
3a72645a 747 if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
748 AliInfo("You have to init before");
749 return kFALSE;
750 }
751
752 if((fStepTrue == 0) && (fStepMC == 0) && (fStepData == 0)) {
753 AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
754 return kFALSE;
755 }
756
c2690925 757 SetNumberOfIteration(10);
3a72645a 758 SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
759
760 AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
761 //////////////////////////////////
762 // Subtract hadron background
763 /////////////////////////////////
67fe7bd0 764 AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
e17c1f86 765 AliCFDataGrid *unnormalizedRawSpectrum = 0x0;
766 TGraphErrors *gNormalizedRawSpectrum = 0x0;
3a72645a 767 if(subtractcontamination) {
a8ef1999 768 if(!fBeauty2ndMethod) dataspectrumaftersubstraction = SubtractBackground(kTRUE);
769 else dataspectrumaftersubstraction = GetRawBspectra2ndMethod();
770 unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone();
771 dataGridAfterFirstSteps = dataspectrumaftersubstraction;
772 gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum);
3a72645a 773 }
774
a8ef1999 775 printf("after normalize getting IP \n");
776
c2690925 777 /////////////////////////////////////////////////////////////////////////////////////////
778 // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
779 /////////////////////////////////////////////////////////////////////////////////////////
780
8c1c76e9 781 AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
3a72645a 782 AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
783 AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
8c1c76e9 784
785 if(fEfficiencyFunction){
786 dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
a8ef1999 787 dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
8c1c76e9 788 }
789 else if(dataContainerV0){
3a72645a 790 dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
791 dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
8c1c76e9 792 }
793
794
a8ef1999 795
3a72645a 796 ///////////////
c04c80e6 797 // Unfold
3a72645a 798 //////////////
799 TList *listunfolded = Unfold(dataGridAfterFirstSteps);
c04c80e6 800 if(!listunfolded){
801 printf("Unfolded failed\n");
3a72645a 802 return kFALSE;
c04c80e6 803 }
804 THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
805 THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
806 if(!correctedspectrum){
807 AliError("No corrected spectrum\n");
3a72645a 808 return kFALSE;
c04c80e6 809 }
67fe7bd0 810 if(!residualspectrum){
8c1c76e9 811 AliError("No residual spectrum\n");
3a72645a 812 return kFALSE;
c04c80e6 813 }
67fe7bd0 814
3a72645a 815 /////////////////////
c04c80e6 816 // Simply correct
3a72645a 817 ////////////////////
a8ef1999 818
3a72645a 819 AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
67fe7bd0 820
3a72645a 821
67fe7bd0 822 //////////
c04c80e6 823 // Plot
824 //////////
a8ef1999 825
3a72645a 826 if(fDebugLevel > 0.0) {
a8ef1999 827
828 Int_t ptpr = 0;
829 if(fBeamType==0) ptpr=0;
830 if(fBeamType==1) ptpr=1;
3a72645a 831
832 TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
833 ccorrected->Divide(2,1);
834 ccorrected->cd(1);
835 gPad->SetLogy();
836 TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
837 correctedspectrumD->SetTitle("");
838 correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
839 correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
840 correctedspectrumD->SetMarkerStyle(26);
841 correctedspectrumD->SetMarkerColor(kBlue);
842 correctedspectrumD->SetLineColor(kBlue);
843 correctedspectrumD->Draw("AP");
844 TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
845 alltogetherspectrumD->SetTitle("");
846 alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
847 alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
848 alltogetherspectrumD->SetMarkerStyle(25);
849 alltogetherspectrumD->SetMarkerColor(kBlack);
850 alltogetherspectrumD->SetLineColor(kBlack);
851 alltogetherspectrumD->Draw("P");
852 TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
853 legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
854 legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
855 legcorrected->Draw("same");
856 ccorrected->cd(2);
a8ef1999 857 TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
858 TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
3a72645a 859 TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
860 ratiocorrected->SetName("ratiocorrected");
861 ratiocorrected->SetTitle("");
862 ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
863 ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
864 ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
865 ratiocorrected->SetStats(0);
866 ratiocorrected->Draw();
e17c1f86 867 if(fWriteToFile) ccorrected->SaveAs("CorrectedBeauty.eps");
868
a8ef1999 869 if(fBeamType == 0){
870 if(fNonHFEsyst){
871 CalculateNonHFEsyst(0);
872 }
e17c1f86 873 }
3a72645a 874
3a72645a 875 // Dump to file if needed
876
877 if(fDumpToFile) {
a8ef1999 878 // to do centrality dependent
879
8c1c76e9 880 TFile *out;
e17c1f86 881 out = new TFile("finalSpectrum.root","recreate");
3a72645a 882 out->cd();
883 //
884 correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
885 correctedspectrumD->Write();
886 alltogetherspectrumD->SetName("AlltogetherSpectrum");
887 alltogetherspectrumD->Write();
888 ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
889 ratiocorrected->Write();
890 //
891 correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
892 correctedspectrum->Write();
893 alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
894 alltogetherCorrection->Write();
895 //
dcef324e 896 if(unnormalizedRawSpectrum) {
11ff28c5 897 unnormalizedRawSpectrum->SetName("beautyAfterIP");
898 unnormalizedRawSpectrum->Write();
dcef324e 899 }
11ff28c5 900
e17c1f86 901 if(gNormalizedRawSpectrum){
902 gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP");
903 gNormalizedRawSpectrum->Write();
904 }
905
11ff28c5 906 if(fBeamType==0) {
907 Int_t countpp=0;
908 fEfficiencyCharmSigD[countpp]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",countpp));
909 fEfficiencyCharmSigD[countpp]->SetName(Form("IPEfficiencyForCharmSigCent%i",countpp));
910 fEfficiencyCharmSigD[countpp]->Write();
911 fEfficiencyBeautySigD[countpp]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",countpp));
912 fEfficiencyBeautySigD[countpp]->SetName(Form("IPEfficiencyForBeautySigCent%i",countpp));
913 fEfficiencyBeautySigD[countpp]->Write();
914 fCharmEff[countpp]->SetTitle(Form("IPEfficiencyForCharmCent%i",countpp));
915 fCharmEff[countpp]->SetName(Form("IPEfficiencyForCharmCent%i",countpp));
916 fCharmEff[countpp]->Write();
917 fBeautyEff[countpp]->SetTitle(Form("IPEfficiencyForBeautyCent%i",countpp));
918 fBeautyEff[countpp]->SetName(Form("IPEfficiencyForBeautyCent%i",countpp));
919 fBeautyEff[countpp]->Write();
920 fConversionEff[countpp]->SetTitle(Form("IPEfficiencyForConversionCent%i",countpp));
921 fConversionEff[countpp]->SetName(Form("IPEfficiencyForConversionCent%i",countpp));
922 fConversionEff[countpp]->Write();
923 fNonHFEEff[countpp]->SetTitle(Form("IPEfficiencyForNonHFECent%i",countpp));
924 fNonHFEEff[countpp]->SetName(Form("IPEfficiencyForNonHFECent%i",countpp));
925 fNonHFEEff[countpp]->Write();
926 }
927
a8ef1999 928 if(fBeamType==1) {
929
930 TGraphErrors* correctedspectrumDc[kCentrality];
931 TGraphErrors* alltogetherspectrumDc[kCentrality];
11ff28c5 932 for(Int_t i=0;i<kCentrality-2;i++)
a8ef1999 933 {
934 correctedspectrum->GetAxis(0)->SetRange(i+1,i+1);
935 correctedspectrumDc[i] = Normalize(correctedspectrum,i);
936 correctedspectrumDc[i]->SetTitle(Form("UnfoldingCorrectedSpectrum_%i",i));
937 correctedspectrumDc[i]->SetName(Form("UnfoldingCorrectedSpectrum_%i",i));
938 alltogetherCorrection->GetAxis(0)->SetRange(i+1,i+1);
939 alltogetherspectrumDc[i] = Normalize(alltogetherCorrection,i);
940 alltogetherspectrumDc[i]->SetTitle(Form("AlltogetherSpectrum_%i",i));
941 alltogetherspectrumDc[i]->SetName(Form("AlltogetherSpectrum_%i",i));
942 correctedspectrumDc[i]->Write();
943 alltogetherspectrumDc[i]->Write();
944
945
946 TH1D *centrcrosscheck = correctedspectrum->Projection(0);
947 centrcrosscheck->SetTitle(Form("centrality_%i",i));
948 centrcrosscheck->SetName(Form("centrality_%i",i));
949 centrcrosscheck->Write();
950
951 TH1D *correctedTH1Dc = correctedspectrum->Projection(ptpr);
952 TH1D *alltogetherTH1Dc = (TH1D *) alltogetherCorrection->Project(ptpr);
953
954 TH1D *centrcrosscheck2 = (TH1D *) alltogetherCorrection->Project(0);
955 centrcrosscheck2->SetTitle(Form("centrality2_%i",i));
956 centrcrosscheck2->SetName(Form("centrality2_%i",i));
957 centrcrosscheck2->Write();
958
959 TH1D* ratiocorrectedc = (TH1D*)correctedTH1D->Clone();
960 ratiocorrectedc->Divide(correctedTH1Dc,alltogetherTH1Dc,1,1);
961 ratiocorrectedc->SetTitle(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
962 ratiocorrectedc->SetName(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
963 ratiocorrectedc->Write();
964
11ff28c5 965 fEfficiencyCharmSigD[i]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",i));
966 fEfficiencyCharmSigD[i]->SetName(Form("IPEfficiencyForCharmSigCent%i",i));
a8ef1999 967 fEfficiencyCharmSigD[i]->Write();
11ff28c5 968 fEfficiencyBeautySigD[i]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",i));
969 fEfficiencyBeautySigD[i]->SetName(Form("IPEfficiencyForBeautySigCent%i",i));
a8ef1999 970 fEfficiencyBeautySigD[i]->Write();
971 fCharmEff[i]->SetTitle(Form("IPEfficiencyForCharmCent%i",i));
972 fCharmEff[i]->SetName(Form("IPEfficiencyForCharmCent%i",i));
973 fCharmEff[i]->Write();
974 fBeautyEff[i]->SetTitle(Form("IPEfficiencyForBeautyCent%i",i));
975 fBeautyEff[i]->SetName(Form("IPEfficiencyForBeautyCent%i",i));
976 fBeautyEff[i]->Write();
977 fConversionEff[i]->SetTitle(Form("IPEfficiencyForConversionCent%i",i));
978 fConversionEff[i]->SetName(Form("IPEfficiencyForConversionCent%i",i));
979 fConversionEff[i]->Write();
980 fNonHFEEff[i]->SetTitle(Form("IPEfficiencyForNonHFECent%i",i));
981 fNonHFEEff[i]->SetName(Form("IPEfficiencyForNonHFECent%i",i));
982 fNonHFEEff[i]->Write();
983 }
984
985 }
986
8c1c76e9 987 out->Close();
988 delete out;
3a72645a 989 }
3a72645a 990 }
991
3a72645a 992 return kTRUE;
993}
c2690925 994
3a72645a 995//____________________________________________________________
996AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
997 //
998 // Apply background subtraction
999 //
a8ef1999 1000
1001 Int_t ptpr = 0;
1002 Int_t nbins=1;
1003 if(fBeamType==0)
1004 {
1005 ptpr=0;
1006 nbins=1;
1007 }
1008 if(fBeamType==1)
1009 {
1010 ptpr=1;
1011 nbins=2;
1012 }
1013
3a72645a 1014 // Raw spectrum
1015 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1016 if(!dataContainer){
1017 AliError("Data Container not available");
1018 return NULL;
1019 }
c2690925 1020 printf("Step data: %d\n",fStepData);
3a72645a 1021 AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
1022
1023 AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
1024 dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
1025
1026 // Background Estimate
1027 AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
1028 if(!backgroundContainer){
1029 AliError("MC background container not found");
1030 return NULL;
1031 }
1032
c2690925 1033 Int_t stepbackground = 1; // 2 for !fInclusiveSpectrum analysis(old method)
3a72645a 1034 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
1035
c2690925 1036 if(!fInclusiveSpectrum){
1037 //Background subtraction for IP analysis
a8ef1999 1038 TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
8c1c76e9 1039 CorrectFromTheWidth(measuredTH1Draw);
1040 TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",600,500);
1041 rawspectra->cd();
e17c1f86 1042 rawspectra->SetLogy();
1043 TLegend *lRaw = new TLegend(0.65,0.65,0.95,0.95);
8c1c76e9 1044 measuredTH1Draw->SetMarkerStyle(20);
1045 measuredTH1Draw->Draw();
e17c1f86 1046 lRaw->AddEntry(measuredTH1Draw,"measured raw spectrum");
a8ef1999 1047 TH1D* htemp;
dcef324e 1048 Int_t* bins=new Int_t[2];
c2690925 1049 if(fIPanaHadronBgSubtract){
1050 // Hadron background
a8ef1999 1051 printf("Hadron background for IP analysis subtracted!\n");
1052 if(fBeamType==0)
1053 {
1054
1055 htemp = (TH1D *) fHadronEffbyIPcut->Projection(0);
1056 bins[0]=htemp->GetNbinsX();
1057 }
1058 if(fBeamType==1)
1059 {
a8ef1999 1060 htemp = (TH1D *) fHadronEffbyIPcut->Projection(0);
1061 bins[0]=htemp->GetNbinsX();
1062 htemp = (TH1D *) fHadronEffbyIPcut->Projection(1);
1063 bins[1]=htemp->GetNbinsX();
1064 }
1065 AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",nbins,bins);
c2690925 1066 hbgContainer->SetGrid(fHadronEffbyIPcut);
1067 backgroundGrid->Multiply(hbgContainer,1);
8c1c76e9 1068 // draw raw hadron bg spectra
a8ef1999 1069 TH1D *hadronbg= (TH1D *) backgroundGrid->Project(ptpr);
8c1c76e9 1070 CorrectFromTheWidth(hadronbg);
1071 hadronbg->SetMarkerColor(7);
1072 hadronbg->SetMarkerStyle(20);
1073 rawspectra->cd();
1074 hadronbg->Draw("samep");
e17c1f86 1075 lRaw->AddEntry(hadronbg,"hadrons");
8c1c76e9 1076 // subtract hadron contamination
c2690925 1077 spectrumSubtracted->Add(backgroundGrid,-1.0);
1078 }
1079 if(fIPanaCharmBgSubtract){
1080 // Charm background
1081 printf("Charm background for IP analysis subtracted!\n");
1082 AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
8c1c76e9 1083 // draw charm bg spectra
a8ef1999 1084 TH1D *charmbg= (TH1D *) charmbgContainer->Project(ptpr);
8c1c76e9 1085 CorrectFromTheWidth(charmbg);
1086 charmbg->SetMarkerColor(3);
1087 charmbg->SetMarkerStyle(20);
1088 rawspectra->cd();
1089 charmbg->Draw("samep");
e17c1f86 1090 lRaw->AddEntry(charmbg,"charm elecs");
8c1c76e9 1091 // subtract charm background
c2690925 1092 spectrumSubtracted->Add(charmbgContainer,-1.0);
1093 }
1094 if(fIPanaConversionBgSubtract){
1095 // Conversion background
a8ef1999 1096 AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground();
8c1c76e9 1097 // draw conversion bg spectra
a8ef1999 1098 TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(ptpr);
8c1c76e9 1099 CorrectFromTheWidth(conversionbg);
1100 conversionbg->SetMarkerColor(4);
1101 conversionbg->SetMarkerStyle(20);
1102 rawspectra->cd();
1103 conversionbg->Draw("samep");
e17c1f86 1104 lRaw->AddEntry(conversionbg,"conversion elecs");
8c1c76e9 1105 // subtract conversion background
c2690925 1106 spectrumSubtracted->Add(conversionbgContainer,-1.0);
1107 printf("Conversion background subtraction is preliminary!\n");
1108 }
1109 if(fIPanaNonHFEBgSubtract){
1110 // NonHFE background
1111 AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground();
8c1c76e9 1112 // draw Dalitz/dielectron bg spectra
a8ef1999 1113 TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(ptpr);
8c1c76e9 1114 CorrectFromTheWidth(nonhfebg);
1115 nonhfebg->SetMarkerColor(6);
1116 nonhfebg->SetMarkerStyle(20);
1117 rawspectra->cd();
1118 nonhfebg->Draw("samep");
e17c1f86 1119 lRaw->AddEntry(nonhfebg,"non-HF elecs");
8c1c76e9 1120 // subtract Dalitz/dielectron background
c2690925 1121 spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
1122 printf("Non HFE background subtraction is preliminary!\n");
1123 }
a8ef1999 1124 TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(ptpr);
8c1c76e9 1125 CorrectFromTheWidth(rawbgsubtracted);
1126 rawbgsubtracted->SetMarkerStyle(24);
1127 rawspectra->cd();
e17c1f86 1128 lRaw->AddEntry(rawbgsubtracted,"subtracted raw spectrum");
8c1c76e9 1129 rawbgsubtracted->Draw("samep");
e17c1f86 1130 lRaw->Draw("SAME");
dcef324e 1131
1132 delete[] bins;
11ff28c5 1133
c2690925 1134 }
1135 else{
1136 // Subtract
1137 spectrumSubtracted->Add(backgroundGrid,-1.0);
1138 }
1139
3a72645a 1140 if(setBackground){
1141 if(fBackground) delete fBackground;
1142 fBackground = backgroundGrid;
1143 } else delete backgroundGrid;
1144
1145
1146 if(fDebugLevel > 0) {
c2690925 1147
a8ef1999 1148 Int_t ptprd;
1149 if(fBeamType==0) ptprd=0;
1150 if(fBeamType==1) ptprd=1;
67fe7bd0 1151
1152 TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
3a72645a 1153 cbackgroundsubtraction->Divide(3,1);
67fe7bd0 1154 cbackgroundsubtraction->cd(1);
3a72645a 1155 gPad->SetLogy();
a8ef1999 1156 TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptprd);
1157 TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptprd);
67fe7bd0 1158 CorrectFromTheWidth(measuredTH1Daftersubstraction);
1159 CorrectFromTheWidth(measuredTH1Dbeforesubstraction);
1160 measuredTH1Daftersubstraction->SetStats(0);
1161 measuredTH1Daftersubstraction->SetTitle("");
1162 measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1163 measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1164 measuredTH1Daftersubstraction->SetMarkerStyle(25);
1165 measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
1166 measuredTH1Daftersubstraction->SetLineColor(kBlack);
1167 measuredTH1Dbeforesubstraction->SetStats(0);
1168 measuredTH1Dbeforesubstraction->SetTitle("");
1169 measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1170 measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1171 measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
1172 measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
1173 measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
1174 measuredTH1Daftersubstraction->Draw();
1175 measuredTH1Dbeforesubstraction->Draw("same");
1176 TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
3a72645a 1177 legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
1178 legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
67fe7bd0 1179 legsubstraction->Draw("same");
1180 cbackgroundsubtraction->cd(2);
3a72645a 1181 gPad->SetLogy();
67fe7bd0 1182 TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
1183 ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
1184 ratiomeasuredcontamination->SetTitle("");
1185 ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1186 ratiomeasuredcontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
c2690925 1187 ratiomeasuredcontamination->Sumw2();
67fe7bd0 1188 ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
1189 ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
1190 ratiomeasuredcontamination->SetStats(0);
1191 ratiomeasuredcontamination->SetMarkerStyle(26);
1192 ratiomeasuredcontamination->SetMarkerColor(kBlack);
1193 ratiomeasuredcontamination->SetLineColor(kBlack);
c2690925 1194 for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
1195 ratiomeasuredcontamination->SetBinError(k+1,0.0);
1196 }
1197 ratiomeasuredcontamination->Draw("P");
3a72645a 1198 cbackgroundsubtraction->cd(3);
a8ef1999 1199 TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptprd);
3a72645a 1200 CorrectFromTheWidth(measuredTH1background);
1201 measuredTH1background->SetStats(0);
1202 measuredTH1background->SetTitle("");
1203 measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1204 measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1205 measuredTH1background->SetMarkerStyle(26);
1206 measuredTH1background->SetMarkerColor(kBlack);
1207 measuredTH1background->SetLineColor(kBlack);
1208 measuredTH1background->Draw();
e17c1f86 1209 if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.eps");
c2690925 1210
1211 if(fBeamType==1) {
1212
1213 TCanvas * cbackgrounde = new TCanvas("BackgroundSubtraction_allspectra","BackgroundSubtraction_allspectra",1000,700);
1214 cbackgrounde->Divide(2,1);
1215 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1216 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1217
1218 THnSparseF* sparsesubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
1219 TAxis *cenaxisa = sparsesubtracted->GetAxis(0);
1220 THnSparseF* sparsebefore = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
1221 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1222 Int_t nbbin = cenaxisb->GetNbins();
1223 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1224 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1225 for(Int_t binc = 0; binc < nbbin; binc++){
e17c1f86 1226 TString titlee("BackgroundSubtraction_centrality_bin_");
1227 titlee += binc;
1228 TCanvas * cbackground = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1229 cbackground->Divide(2,1);
1230 cbackground->cd(1);
1231 gPad->SetLogy();
1232 cenaxisa->SetRange(binc+1,binc+1);
1233 cenaxisb->SetRange(binc+1,binc+1);
1234 TH1D *aftersubstraction = (TH1D *) sparsesubtracted->Projection(1);
1235 TH1D *beforesubstraction = (TH1D *) sparsebefore->Projection(1);
1236 CorrectFromTheWidth(aftersubstraction);
1237 CorrectFromTheWidth(beforesubstraction);
1238 aftersubstraction->SetStats(0);
1239 aftersubstraction->SetTitle((const char*)titlee);
1240 aftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1241 aftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1242 aftersubstraction->SetMarkerStyle(25);
1243 aftersubstraction->SetMarkerColor(kBlack);
1244 aftersubstraction->SetLineColor(kBlack);
1245 beforesubstraction->SetStats(0);
1246 beforesubstraction->SetTitle((const char*)titlee);
1247 beforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1248 beforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1249 beforesubstraction->SetMarkerStyle(24);
1250 beforesubstraction->SetMarkerColor(kBlue);
1251 beforesubstraction->SetLineColor(kBlue);
1252 aftersubstraction->Draw();
1253 beforesubstraction->Draw("same");
1254 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1255 lega->AddEntry(beforesubstraction,"With hadron contamination","p");
1256 lega->AddEntry(aftersubstraction,"Without hadron contamination ","p");
1257 lega->Draw("same");
1258 cbackgrounde->cd(1);
1259 gPad->SetLogy();
1260 TH1D *aftersubtractionn = (TH1D *) aftersubstraction->Clone();
1261 aftersubtractionn->SetMarkerStyle(stylee[binc]);
1262 aftersubtractionn->SetMarkerColor(colorr[binc]);
1263 if(binc==0) aftersubtractionn->Draw();
1264 else aftersubtractionn->Draw("same");
1265 legtotal->AddEntry(aftersubtractionn,(const char*) titlee,"p");
1266 cbackgrounde->cd(2);
1267 gPad->SetLogy();
1268 TH1D *aftersubtractionng = (TH1D *) aftersubstraction->Clone();
1269 aftersubtractionng->SetMarkerStyle(stylee[binc]);
1270 aftersubtractionng->SetMarkerColor(colorr[binc]);
1271 if(fNEvents[binc] > 0.0) aftersubtractionng->Scale(1/(Double_t)fNEvents[binc]);
1272 if(binc==0) aftersubtractionng->Draw();
1273 else aftersubtractionng->Draw("same");
1274 legtotalg->AddEntry(aftersubtractionng,(const char*) titlee,"p");
1275 cbackground->cd(2);
1276 TH1D* ratiocontamination = (TH1D*)beforesubstraction->Clone();
1277 ratiocontamination->SetName("ratiocontamination");
1278 ratiocontamination->SetTitle((const char*)titlee);
1279 ratiocontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1280 ratiocontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1281 ratiocontamination->Add(aftersubstraction,-1.0);
1282 ratiocontamination->Divide(beforesubstraction);
1283 Int_t totalbin = ratiocontamination->GetXaxis()->GetNbins();
1284 for(Int_t nbinpt = 0; nbinpt < totalbin; nbinpt++) {
1285 ratiocontamination->SetBinError(nbinpt+1,0.0);
1286 }
1287 ratiocontamination->SetStats(0);
1288 ratiocontamination->SetMarkerStyle(26);
1289 ratiocontamination->SetMarkerColor(kBlack);
1290 ratiocontamination->SetLineColor(kBlack);
1291 ratiocontamination->Draw("P");
c2690925 1292 }
1293
1294 cbackgrounde->cd(1);
1295 legtotal->Draw("same");
1296 cbackgrounde->cd(2);
1297 legtotalg->Draw("same");
1298
1299 cenaxisa->SetRange(0,nbbin);
1300 cenaxisb->SetRange(0,nbbin);
e17c1f86 1301 if(fWriteToFile) cbackgrounde->SaveAs("BackgroundSubtractedPbPb.eps");
c2690925 1302 }
67fe7bd0 1303 }
1304
3a72645a 1305 return spectrumSubtracted;
c04c80e6 1306}
c2690925 1307
1308//____________________________________________________________
1309AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
1310 //
1311 // calculate charm background
1312 //
a8ef1999 1313 Int_t ptpr = 0;
1314 Int_t nDim = 1;
1315 if(fBeamType==0)
1316 {
1317 ptpr=0;
1318 }
1319 if(fBeamType==1)
1320 {
1321 ptpr=1;
1322 nDim=2;
1323 }
c2690925 1324
1325 Double_t evtnorm=0;
a8ef1999 1326 if(fNMCbgEvents[0]) evtnorm= double(fNEvents[0])/double(fNMCbgEvents[0]);
c2690925 1327
1328 AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
1329 if(!mcContainer){
1330 AliError("MC Container not available");
1331 return NULL;
1332 }
1333
1334 if(!fCorrelation){
1335 AliError("No Correlation map available");
1336 return NULL;
1337 }
1338
c2690925 1339 AliCFDataGrid *charmBackgroundGrid= 0x0;
8c1c76e9 1340 charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
8c1c76e9 1341
a8ef1999 1342 TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
dcef324e 1343 Int_t* bins=new Int_t[2];
8c1c76e9 1344 bins[0]=charmbgaftertofpid->GetNbinsX();
8c1c76e9 1345
a8ef1999 1346 if(fBeamType==1)
1347 {
3024f297 1348 //charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
11ff28c5 1349 bins[0]=12;
a8ef1999 1350 charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(1);
1351 bins[1]=charmbgaftertofpid->GetNbinsX();
8c1c76e9 1352
a8ef1999 1353 }
1354
1355 AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins);
1356 ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
1357 TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(ptpr);
1358
1359 if(fBeamType==0)charmBackgroundGrid->Multiply(ipWeightedCharmContainer,evtnorm);
1360 Double_t contents[2];
1361 AliCFDataGrid *eventTemp = new AliCFDataGrid("eventTemp","eventTemp",nDim,bins);
1362 if(fBeamType==1)
1363 {
1364 for(Int_t kCentr=0;kCentr<bins[0];kCentr++)
1365 {
1366 for(Int_t kpt=0;kpt<bins[1];kpt++)
1367 {
1368 Double_t evtnormPbPb=0;
1369 if(fNMCbgEvents[kCentr]) evtnormPbPb= double(fNEvents[kCentr])/double(fNMCbgEvents[kCentr]);
1370 contents[0]=kCentr;
1371 contents[1]=kpt;
1372 eventTemp->Fill(contents,evtnormPbPb);
1373 }
1374 }
1375 charmBackgroundGrid->Multiply(eventTemp,1);
1376 }
1377 TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(ptpr);
1378
1379 AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
8c1c76e9 1380 weightedCharmContainer->SetGrid(GetCharmWeights()); // get charm weighting factors
a8ef1999 1381 TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(ptpr);
8c1c76e9 1382 charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
a8ef1999 1383 TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(ptpr);
c2690925 1384
1385 // Efficiency (set efficiency to 1 for only folding)
1386 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1387 efficiencyD->CalculateEfficiency(0,0);
1388
a8ef1999 1389 // Folding
1390 if(fBeamType==0)nDim = 1;
1391 if(fBeamType==1)nDim = 2;
c2690925 1392 AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
1393 folding.SetMaxNumberOfIterations(1);
1394 folding.Unfold();
1395
1396 // Results
1397 THnSparse* result1= folding.GetEstMeasured(); // folded spectra
1398 THnSparse* result=(THnSparse*)result1->Clone();
1399 charmBackgroundGrid->SetGrid(result);
a8ef1999 1400 TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(ptpr);
8c1c76e9 1401
1402 //Charm background evaluation plots
1403
1404 TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
1405 cCharmBgEval->Divide(3,1);
1406
1407 cCharmBgEval->cd(1);
a8ef1999 1408
1409 if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm);
1410 if(fBeamType==1)
1411 {
1412 Double_t evtnormPbPb=0;
1413 for(Int_t kCentr=0;kCentr<bins[0];kCentr++)
1414 {
1415 if(fNMCbgEvents[kCentr]) evtnormPbPb= evtnormPbPb+double(fNEvents[kCentr])/double(fNMCbgEvents[kCentr]);
1416 }
1417 charmbgaftertofpid->Scale(evtnormPbPb);
1418 }
1419
8c1c76e9 1420 CorrectFromTheWidth(charmbgaftertofpid);
1421 charmbgaftertofpid->SetMarkerStyle(25);
1422 charmbgaftertofpid->Draw("p");
1423 charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
1424 charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1425 gPad->SetLogy();
1426
1427 CorrectFromTheWidth(charmbgafteripcut);
1428 charmbgafteripcut->SetMarkerStyle(24);
1429 charmbgafteripcut->Draw("samep");
1430
1431 CorrectFromTheWidth(charmbgafterweight);
1432 charmbgafterweight->SetMarkerStyle(24);
1433 charmbgafterweight->SetMarkerColor(4);
1434 charmbgafterweight->Draw("samep");
1435
1436 CorrectFromTheWidth(charmbgafterfolding);
1437 charmbgafterfolding->SetMarkerStyle(24);
1438 charmbgafterfolding->SetMarkerColor(2);
1439 charmbgafterfolding->Draw("samep");
1440
1441 cCharmBgEval->cd(2);
1442 parametrizedcharmpidipeff->SetMarkerStyle(24);
1443 parametrizedcharmpidipeff->Draw("p");
1444 parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1445
1446 cCharmBgEval->cd(3);
1447 charmweightingfc->SetMarkerStyle(24);
1448 charmweightingfc->Draw("p");
1449 charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
1450 charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1451
1452 cCharmBgEval->cd(1);
1453 TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
1454 legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
1455 legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
1456 legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
1457 legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
1458 legcharmbg->Draw("same");
1459
1460 cCharmBgEval->cd(2);
1461 TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
1462 legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
1463 legcharmbg2->Draw("same");
1464
1465 CorrectStatErr(charmBackgroundGrid);
e17c1f86 1466 if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
c2690925 1467
dcef324e 1468 delete[] bins;
1469
c2690925 1470 return charmBackgroundGrid;
1471}
1472
1473//____________________________________________________________
1474AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
1475 //
1476 // calculate conversion background
1477 //
e17c1f86 1478
a199006c 1479 Double_t evtnorm[1] = {0.0};
a8ef1999 1480 if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
c2690925 1481 printf("check event!!! %lf \n",evtnorm[0]);
e17c1f86 1482
1483 AliCFContainer *backgroundContainer = 0x0;
1484
1485 if(fNonHFEsyst){
a8ef1999 1486 backgroundContainer = (AliCFContainer*)fConvSourceContainer[0][0][0]->Clone();
e17c1f86 1487 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
a8ef1999 1488 backgroundContainer->Add(fConvSourceContainer[iSource][0][0]); // make centrality dependent
e17c1f86 1489 }
1490 }
1491 else{
1492 // Background Estimate
1493 backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
1494 }
c2690925 1495 if(!backgroundContainer){
1496 AliError("MC background container not found");
1497 return NULL;
1498 }
e17c1f86 1499
8c1c76e9 1500 Int_t stepbackground = 3;
a8ef1999 1501 if(TMath::Abs(fEtaRange[0]) < 0.5) stepbackground = 4;
1502
c2690925 1503 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground);
a8ef1999 1504 Int_t *nBinpp=new Int_t[1];
1505 Int_t* binspp=new Int_t[1];
1506 binspp[0]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1507
1508 Int_t *nBinPbPb=new Int_t[2];
1509 Int_t* binsPbPb=new Int_t[2];
1510 binsPbPb[1]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1511 binsPbPb[0]=12;
1512
1513 Int_t looppt=binspp[0];
1514 if(fBeamType==1) looppt=binsPbPb[1];
1515
1516 for(Long_t iBin=1; iBin<= looppt;iBin++){
1517 if(fBeamType==0)
1518 {
1519 nBinpp[0]=iBin;
1520 backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
1521 backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
1522 }
1523 if(fBeamType==1)
1524 {
1525 // loop over centrality
1526 for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){
1527 nBinPbPb[0]=iiBin;
1528 nBinPbPb[1]=iBin;
1529 Double_t evtnormPbPb=0;
11ff28c5 1530 if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
a8ef1999 1531 backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1532 backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1533 }
1534 }
e17c1f86 1535 }
1536 //end of workaround for statistical errors
a8ef1999 1537
1538 AliCFDataGrid *weightedConversionContainer;
1539 if(fBeamType==0) weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,binspp);
1540 else weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",2,binsPbPb);
8c1c76e9 1541 weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
1542 backgroundGrid->Multiply(weightedConversionContainer,1.0);
11ff28c5 1543
dcef324e 1544 delete[] nBinpp;
1545 delete[] binspp;
1546 delete[] nBinPbPb;
11ff28c5 1547 delete[] binsPbPb;
dcef324e 1548
c2690925 1549 return backgroundGrid;
1550}
1551
1552
1553//____________________________________________________________
1554AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
1555 //
8c1c76e9 1556 // calculate non-HFE background
c2690925 1557 //
1558
a199006c 1559 Double_t evtnorm[1] = {0.0};
a8ef1999 1560 if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
8c1c76e9 1561 printf("check event!!! %lf \n",evtnorm[0]);
e17c1f86 1562
1563 AliCFContainer *backgroundContainer = 0x0;
1564 if(fNonHFEsyst){
a8ef1999 1565 backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[0][0][0]->Clone();
e17c1f86 1566 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
a8ef1999 1567 backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0]);
e17c1f86 1568 }
1569 }
1570 else{
1571 // Background Estimate
1572 backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
1573 }
8c1c76e9 1574 if(!backgroundContainer){
1575 AliError("MC background container not found");
1576 return NULL;
c2690925 1577 }
e17c1f86 1578
1579
8c1c76e9 1580 Int_t stepbackground = 3;
a8ef1999 1581 if(TMath::Abs(fEtaRange[0]) < 0.5) stepbackground = 4;
1582
8c1c76e9 1583 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
a8ef1999 1584 Int_t *nBinpp=new Int_t[1];
1585 Int_t* binspp=new Int_t[1];
1586 binspp[0]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1587
1588 Int_t *nBinPbPb=new Int_t[2];
1589 Int_t* binsPbPb=new Int_t[2];
1590 binsPbPb[1]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1591 binsPbPb[0]=12;
1592
1593 Int_t looppt=binspp[0];
1594 if(fBeamType==1) looppt=binsPbPb[1];
1595
1596
1597 for(Long_t iBin=1; iBin<= looppt;iBin++){
1598 if(fBeamType==0)
1599 {
1600 nBinpp[0]=iBin;
1601 backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
1602 backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
1603 }
1604 if(fBeamType==1)
1605 {
1606 for(Long_t iiBin=1; iiBin<=binsPbPb[0];iiBin++){
1607 nBinPbPb[0]=iiBin;
1608 nBinPbPb[1]=iBin;
1609 Double_t evtnormPbPb=0;
11ff28c5 1610 if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
a8ef1999 1611 backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1612 backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1613 }
1614 }
e17c1f86 1615 }
1616 //end of workaround for statistical errors
a8ef1999 1617 AliCFDataGrid *weightedNonHFEContainer;
1618 if(fBeamType==0) weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,binspp);
1619 else weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",2,binsPbPb);
8c1c76e9 1620 weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1621 backgroundGrid->Multiply(weightedNonHFEContainer,1.0);
c2690925 1622
11ff28c5 1623 delete[] nBinpp;
dcef324e 1624 delete[] binspp;
1625 delete[] nBinPbPb;
11ff28c5 1626 delete[] binsPbPb;
dcef324e 1627
8c1c76e9 1628 return backgroundGrid;
c2690925 1629}
1630
1631//____________________________________________________________
1632AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
1633
1634 //
1635 // Apply TPC pid efficiency correction from parametrisation
1636 //
1637
1638 // Data in the right format
1639 AliCFDataGrid *dataGrid = 0x0;
1640 if(bgsubpectrum) {
1641 dataGrid = bgsubpectrum;
1642 }
1643 else {
1644
1645 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1646 if(!dataContainer){
1647 AliError("Data Container not available");
1648 return NULL;
1649 }
c2690925 1650 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1651 }
c2690925 1652 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1653 result->SetName("ParametrizedEfficiencyBefore");
1654 THnSparse *h = result->GetGrid();
1655 Int_t nbdimensions = h->GetNdimensions();
1656 //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
c2690925 1657 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1658 if(!dataContainer){
1659 AliError("Data Container not available");
1660 return NULL;
1661 }
1662 AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
1663 dataContainerbis->Add(dataContainerbis,-1.0);
1664
1665
1666 Int_t* coord = new Int_t[nbdimensions];
1667 memset(coord, 0, sizeof(Int_t) * nbdimensions);
1668 Double_t* points = new Double_t[nbdimensions];
1669
c2690925 1670 ULong64_t nEntries = h->GetNbins();
1671 for (ULong64_t i = 0; i < nEntries; ++i) {
1672
1673 Double_t value = h->GetBinContent(i, coord);
1674 //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
1675 //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
1676
1677 // Get the bin co-ordinates given an coord
1678 for (Int_t j = 0; j < nbdimensions; ++j)
1679 points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1680
1681 if (!fEfficiencyFunction->IsInside(points))
1682 continue;
1683 TF1::RejectPoint(kFALSE);
1684
1685 // Evaulate function at points
1686 Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1687 //printf("Value efficiency is %f\n",valueEfficiency);
1688
1689 if(valueEfficiency > 0.0) {
1690 h->SetBinContent(coord,value/valueEfficiency);
1691 dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1692 }
1693 Double_t error = h->GetBinError(i);
1694 h->SetBinError(coord,error/valueEfficiency);
1695 dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1696
1697
1698 }
1699
6c70d827 1700 delete[] coord;
1701 delete[] points;
1702
c2690925 1703 AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1704
1705 if(fDebugLevel > 0) {
a8ef1999 1706
c2690925 1707 TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700);
1708 cEfficiencyParametrized->Divide(2,1);
1709 cEfficiencyParametrized->cd(1);
1710 TH1D *afterE = (TH1D *) resultt->Project(0);
1711 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1712 CorrectFromTheWidth(afterE);
1713 CorrectFromTheWidth(beforeE);
1714 afterE->SetStats(0);
1715 afterE->SetTitle("");
1716 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1717 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1718 afterE->SetMarkerStyle(25);
1719 afterE->SetMarkerColor(kBlack);
1720 afterE->SetLineColor(kBlack);
1721 beforeE->SetStats(0);
1722 beforeE->SetTitle("");
1723 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1724 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1725 beforeE->SetMarkerStyle(24);
1726 beforeE->SetMarkerColor(kBlue);
1727 beforeE->SetLineColor(kBlue);
1728 gPad->SetLogy();
1729 afterE->Draw();
1730 beforeE->Draw("same");
1731 TLegend *legefficiencyparametrized = new TLegend(0.4,0.6,0.89,0.89);
1732 legefficiencyparametrized->AddEntry(beforeE,"Before Efficiency correction","p");
1733 legefficiencyparametrized->AddEntry(afterE,"After Efficiency correction","p");
1734 legefficiencyparametrized->Draw("same");
1735 cEfficiencyParametrized->cd(2);
1736 fEfficiencyFunction->Draw();
1737 //cEfficiencyParametrized->cd(3);
1738 //TH1D *ratioefficiency = (TH1D *) beforeE->Clone();
1739 //ratioefficiency->Divide(afterE);
1740 //ratioefficiency->Draw();
1741
e17c1f86 1742 if(fWriteToFile) cEfficiencyParametrized->SaveAs("efficiency.eps");
a8ef1999 1743
c2690925 1744 }
1745
c2690925 1746 return resultt;
1747
1748}
c04c80e6 1749//____________________________________________________________
3a72645a 1750AliCFDataGrid *AliHFEspectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
1751
c04c80e6 1752 //
3a72645a 1753 // Apply TPC pid efficiency correction from V0
c04c80e6 1754 //
1755
3a72645a 1756 AliCFContainer *v0Container = GetContainer(kDataContainerV0);
1757 if(!v0Container){
1758 AliError("V0 Container not available");
c04c80e6 1759 return NULL;
1760 }
3a72645a 1761
1762 // Efficiency
1763 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
1764 efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
1765
1766 // Data in the right format
1767 AliCFDataGrid *dataGrid = 0x0;
1768 if(bgsubpectrum) {
1769 dataGrid = bgsubpectrum;
c04c80e6 1770 }
3a72645a 1771 else {
1772
1773 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1774 if(!dataContainer){
1775 AliError("Data Container not available");
1776 return NULL;
1777 }
c04c80e6 1778
3a72645a 1779 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1780 }
c04c80e6 1781
3a72645a 1782 // Correct
1783 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1784 result->ApplyEffCorrection(*efficiencyD);
c04c80e6 1785
3a72645a 1786 if(fDebugLevel > 0) {
c2690925 1787
1788 Int_t ptpr;
1789 if(fBeamType==0) ptpr=0;
1790 if(fBeamType==1) ptpr=1;
3a72645a 1791
1792 TCanvas * cV0Efficiency = new TCanvas("V0Efficiency","V0Efficiency",1000,700);
1793 cV0Efficiency->Divide(2,1);
1794 cV0Efficiency->cd(1);
c2690925 1795 TH1D *afterE = (TH1D *) result->Project(ptpr);
1796 TH1D *beforeE = (TH1D *) dataGrid->Project(ptpr);
3a72645a 1797 afterE->SetStats(0);
1798 afterE->SetTitle("");
1799 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1800 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1801 afterE->SetMarkerStyle(25);
1802 afterE->SetMarkerColor(kBlack);
1803 afterE->SetLineColor(kBlack);
1804 beforeE->SetStats(0);
1805 beforeE->SetTitle("");
1806 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1807 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1808 beforeE->SetMarkerStyle(24);
1809 beforeE->SetMarkerColor(kBlue);
1810 beforeE->SetLineColor(kBlue);
1811 afterE->Draw();
1812 beforeE->Draw("same");
1813 TLegend *legV0efficiency = new TLegend(0.4,0.6,0.89,0.89);
1814 legV0efficiency->AddEntry(beforeE,"Before Efficiency correction","p");
1815 legV0efficiency->AddEntry(afterE,"After Efficiency correction","p");
1816 legV0efficiency->Draw("same");
1817 cV0Efficiency->cd(2);
c2690925 1818 TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(ptpr);
3a72645a 1819 efficiencyDproj->SetTitle("");
1820 efficiencyDproj->SetStats(0);
1821 efficiencyDproj->SetMarkerStyle(25);
1822 efficiencyDproj->Draw();
c04c80e6 1823
c2690925 1824 if(fBeamType==1) {
1825
1826 TCanvas * cV0Efficiencye = new TCanvas("V0Efficiency_allspectra","V0Efficiency_allspectra",1000,700);
1827 cV0Efficiencye->Divide(2,1);
1828 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1829 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1830
1831 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1832 TAxis *cenaxisa = sparseafter->GetAxis(0);
1833 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1834 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1835 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1836 TAxis *cenaxisc = efficiencya->GetAxis(0);
1837 Int_t nbbin = cenaxisb->GetNbins();
1838 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1839 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1840 for(Int_t binc = 0; binc < nbbin; binc++){
e17c1f86 1841 TString titlee("V0Efficiency_centrality_bin_");
1842 titlee += binc;
1843 TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1844 ccV0Efficiency->Divide(2,1);
1845 ccV0Efficiency->cd(1);
1846 gPad->SetLogy();
1847 cenaxisa->SetRange(binc+1,binc+1);
1848 cenaxisb->SetRange(binc+1,binc+1);
1849 cenaxisc->SetRange(binc+1,binc+1);
1850 TH1D *aftere = (TH1D *) sparseafter->Projection(1);
1851 TH1D *beforee = (TH1D *) sparsebefore->Projection(1);
1852 CorrectFromTheWidth(aftere);
1853 CorrectFromTheWidth(beforee);
1854 aftere->SetStats(0);
1855 aftere->SetTitle((const char*)titlee);
1856 aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1857 aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1858 aftere->SetMarkerStyle(25);
1859 aftere->SetMarkerColor(kBlack);
1860 aftere->SetLineColor(kBlack);
1861 beforee->SetStats(0);
1862 beforee->SetTitle((const char*)titlee);
1863 beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1864 beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1865 beforee->SetMarkerStyle(24);
1866 beforee->SetMarkerColor(kBlue);
1867 beforee->SetLineColor(kBlue);
1868 aftere->Draw();
1869 beforee->Draw("same");
1870 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1871 lega->AddEntry(beforee,"Before correction","p");
1872 lega->AddEntry(aftere,"After correction","p");
1873 lega->Draw("same");
1874 cV0Efficiencye->cd(1);
1875 gPad->SetLogy();
1876 TH1D *afteree = (TH1D *) aftere->Clone();
1877 afteree->SetMarkerStyle(stylee[binc]);
1878 afteree->SetMarkerColor(colorr[binc]);
1879 if(binc==0) afteree->Draw();
1880 else afteree->Draw("same");
1881 legtotal->AddEntry(afteree,(const char*) titlee,"p");
1882 cV0Efficiencye->cd(2);
1883 gPad->SetLogy();
1884 TH1D *aftereeu = (TH1D *) aftere->Clone();
1885 aftereeu->SetMarkerStyle(stylee[binc]);
1886 aftereeu->SetMarkerColor(colorr[binc]);
1887 if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]);
1888 if(binc==0) aftereeu->Draw();
1889 else aftereeu->Draw("same");
1890 legtotalg->AddEntry(aftereeu,(const char*) titlee,"p");
1891 ccV0Efficiency->cd(2);
1892 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1);
1893 efficiencyDDproj->SetTitle("");
1894 efficiencyDDproj->SetStats(0);
1895 efficiencyDDproj->SetMarkerStyle(25);
1896 efficiencyDDproj->Draw();
c2690925 1897 }
1898
1899 cV0Efficiencye->cd(1);
1900 legtotal->Draw("same");
1901 cV0Efficiencye->cd(2);
1902 legtotalg->Draw("same");
1903
1904 cenaxisa->SetRange(0,nbbin);
1905 cenaxisb->SetRange(0,nbbin);
1906 cenaxisc->SetRange(0,nbbin);
e17c1f86 1907
1908 if(fWriteToFile) cV0Efficiencye->SaveAs("V0efficiency.eps");
c2690925 1909 }
3a72645a 1910
1911 }
1912
1913
1914 return result;
1915
1916}
c04c80e6 1917//____________________________________________________________
3a72645a 1918TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
c04c80e6 1919
1920 //
1921 // Unfold and eventually correct for efficiency the bgsubspectrum
1922 //
1923
3a72645a 1924 AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
c04c80e6 1925 if(!mcContainer){
1926 AliError("MC Container not available");
1927 return NULL;
1928 }
1929
1930 if(!fCorrelation){
1931 AliError("No Correlation map available");
1932 return NULL;
1933 }
1934
3a72645a 1935 // Data
c04c80e6 1936 AliCFDataGrid *dataGrid = 0x0;
1937 if(bgsubpectrum) {
c04c80e6 1938 dataGrid = bgsubpectrum;
c04c80e6 1939 }
1940 else {
1941
1942 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1943 if(!dataContainer){
1944 AliError("Data Container not available");
1945 return NULL;
1946 }
1947
3a72645a 1948 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
c04c80e6 1949 }
1950
c04c80e6 1951 // Guessed
3a72645a 1952 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
c04c80e6 1953 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
1954
1955 // Efficiency
3a72645a 1956 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
c04c80e6 1957 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
8c1c76e9 1958
a8ef1999 1959 if(!fBeauty2ndMethod)
1960 {
1961 // Consider parameterized IP cut efficiency
1962 if(!fInclusiveSpectrum){
1963 Int_t* bins=new Int_t[1];
1964 bins[0]=35;
1965
1966 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1967 beffContainer->SetGrid(GetBeautyIPEff());
1968 efficiencyD->Multiply(beffContainer,1);
1969 }
8c1c76e9 1970 }
1971
c04c80e6 1972
1973 // Unfold
1974
3a72645a 1975 AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
c2690925 1976 //unfolding.SetUseCorrelatedErrors();
1977 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
c04c80e6 1978 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
e156c3bb 1979 if(fSetSmoothing) unfolding.UseSmoothing();
c04c80e6 1980 unfolding.Unfold();
1981
1982 // Results
1983 THnSparse* result = unfolding.GetUnfolded();
1984 THnSparse* residual = unfolding.GetEstMeasured();
1985 TList *listofresults = new TList;
1986 listofresults->SetOwner();
1987 listofresults->AddAt((THnSparse*)result->Clone(),0);
1988 listofresults->AddAt((THnSparse*)residual->Clone(),1);
1989
3a72645a 1990 if(fDebugLevel > 0) {
c2690925 1991
1992 Int_t ptpr;
1993 if(fBeamType==0) ptpr=0;
1994 if(fBeamType==1) ptpr=1;
3a72645a 1995
1996 TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
1997 cresidual->Divide(2,1);
1998 cresidual->cd(1);
1999 gPad->SetLogy();
2000 TGraphErrors* residualspectrumD = Normalize(residual);
2001 if(!residualspectrumD) {
2002 AliError("Number of Events not set for the normalization");
3ccf8f4c 2003 return NULL;
3a72645a 2004 }
2005 residualspectrumD->SetTitle("");
2006 residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
2007 residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
2008 residualspectrumD->SetMarkerStyle(26);
2009 residualspectrumD->SetMarkerColor(kBlue);
2010 residualspectrumD->SetLineColor(kBlue);
2011 residualspectrumD->Draw("AP");
2012 AliCFDataGrid *dataGridBis = (AliCFDataGrid *) (dataGrid->Clone());
2013 dataGridBis->SetName("dataGridBis");
2014 TGraphErrors* measuredspectrumD = Normalize(dataGridBis);
2015 measuredspectrumD->SetTitle("");
2016 measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
2017 measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
2018 measuredspectrumD->SetMarkerStyle(25);
2019 measuredspectrumD->SetMarkerColor(kBlack);
2020 measuredspectrumD->SetLineColor(kBlack);
2021 measuredspectrumD->Draw("P");
2022 TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
2023 legres->AddEntry(residualspectrumD,"Residual","p");
2024 legres->AddEntry(measuredspectrumD,"Measured","p");
2025 legres->Draw("same");
2026 cresidual->cd(2);
c2690925 2027 TH1D *residualTH1D = residual->Projection(ptpr);
2028 TH1D *measuredTH1D = (TH1D *) dataGridBis->Project(ptpr);
3a72645a 2029 TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
2030 ratioresidual->SetName("ratioresidual");
2031 ratioresidual->SetTitle("");
2032 ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
2033 ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2034 ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
2035 ratioresidual->SetStats(0);
2036 ratioresidual->Draw();
e17c1f86 2037
2038 if(fWriteToFile) cresidual->SaveAs("Unfolding.eps");
3a72645a 2039 }
2040
c04c80e6 2041 return listofresults;
2042
2043}
2044
2045//____________________________________________________________
3a72645a 2046AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
c04c80e6 2047
2048 //
2049 // Apply unfolding and efficiency correction together to bgsubspectrum
2050 //
2051
3a72645a 2052 AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
c04c80e6 2053 if(!mcContainer){
2054 AliError("MC Container not available");
2055 return NULL;
2056 }
2057
c04c80e6 2058 // Efficiency
3a72645a 2059 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
c04c80e6 2060 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
2061
a8ef1999 2062
2063 if(!fBeauty2ndMethod)
2064 {
2065 // Consider parameterized IP cut efficiency
2066 if(!fInclusiveSpectrum){
2067 Int_t* bins=new Int_t[1];
2068 bins[0]=35;
2069
2070 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
2071 beffContainer->SetGrid(GetBeautyIPEff());
2072 efficiencyD->Multiply(beffContainer,1);
2073 }
8c1c76e9 2074 }
2075
c04c80e6 2076 // Data in the right format
2077 AliCFDataGrid *dataGrid = 0x0;
2078 if(bgsubpectrum) {
c04c80e6 2079 dataGrid = bgsubpectrum;
c04c80e6 2080 }
2081 else {
3a72645a 2082
c04c80e6 2083 AliCFContainer *dataContainer = GetContainer(kDataContainer);
2084 if(!dataContainer){
2085 AliError("Data Container not available");
2086 return NULL;
2087 }
2088
3a72645a 2089 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
c04c80e6 2090 }
2091
2092 // Correct
2093 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
2094 result->ApplyEffCorrection(*efficiencyD);
c04c80e6 2095
3a72645a 2096 if(fDebugLevel > 0) {
c2690925 2097
2098 Int_t ptpr;
2099 if(fBeamType==0) ptpr=0;
2100 if(fBeamType==1) ptpr=1;
3a72645a 2101
bf892a6a 2102 printf("Step MC: %d\n",fStepMC);
2103 printf("Step tracking: %d\n",AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
2104 printf("Step MC true: %d\n",fStepTrue);
3a72645a 2105 AliCFEffGrid *efficiencymcPID = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
2106 AliCFEffGrid *efficiencymctrackinggeo = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack,fStepTrue);
2107 AliCFEffGrid *efficiencymcall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,fStepTrue);
2108
2109 AliCFEffGrid *efficiencyesdall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerESD),fStepMC,fStepTrue);
2110
2111 TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
2112 cefficiency->cd(1);
c2690925 2113 TH1D* efficiencymcPIDD = (TH1D *) efficiencymcPID->Project(ptpr);
3a72645a 2114 efficiencymcPIDD->SetTitle("");
2115 efficiencymcPIDD->SetStats(0);
2116 efficiencymcPIDD->SetMarkerStyle(25);
2117 efficiencymcPIDD->Draw();
c2690925 2118 TH1D* efficiencymctrackinggeoD = (TH1D *) efficiencymctrackinggeo->Project(ptpr);
3a72645a 2119 efficiencymctrackinggeoD->SetTitle("");
2120 efficiencymctrackinggeoD->SetStats(0);
2121 efficiencymctrackinggeoD->SetMarkerStyle(26);
2122 efficiencymctrackinggeoD->Draw("same");
c2690925 2123 TH1D* efficiencymcallD = (TH1D *) efficiencymcall->Project(ptpr);
3a72645a 2124 efficiencymcallD->SetTitle("");
2125 efficiencymcallD->SetStats(0);
2126 efficiencymcallD->SetMarkerStyle(27);
2127 efficiencymcallD->Draw("same");
c2690925 2128 TH1D* efficiencyesdallD = (TH1D *) efficiencyesdall->Project(ptpr);
3a72645a 2129 efficiencyesdallD->SetTitle("");
2130 efficiencyesdallD->SetStats(0);
2131 efficiencyesdallD->SetMarkerStyle(24);
2132 efficiencyesdallD->Draw("same");
2133 TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
2134 legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
2135 legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
2136 legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
2137 legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
2138 legeff->Draw("same");
c2690925 2139
2140 if(fBeamType==1) {
2141
2142 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
2143 TAxis *cenaxisa = sparseafter->GetAxis(0);
2144 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
2145 TAxis *cenaxisb = sparsebefore->GetAxis(0);
2146 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
2147 TAxis *cenaxisc = efficiencya->GetAxis(0);
2148 //Int_t nbbin = cenaxisb->GetNbins();
2149 //Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
2150 //Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
2151 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
e17c1f86 2152 TString titlee("Efficiency_centrality_bin_");
2153 titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
2154 titlee += "_";
2155 titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
2156 TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
2157 cefficiencye->Divide(2,1);
2158 cefficiencye->cd(1);
2159 gPad->SetLogy();
2160 cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
2161 cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
2162 TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr);
2163 TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr);
2164 CorrectFromTheWidth(afterefficiency);
2165 CorrectFromTheWidth(beforeefficiency);
2166 afterefficiency->SetStats(0);
2167 afterefficiency->SetTitle((const char*)titlee);
2168 afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2169 afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2170 afterefficiency->SetMarkerStyle(25);
2171 afterefficiency->SetMarkerColor(kBlack);
2172 afterefficiency->SetLineColor(kBlack);
2173 beforeefficiency->SetStats(0);
2174 beforeefficiency->SetTitle((const char*)titlee);
2175 beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2176 beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2177 beforeefficiency->SetMarkerStyle(24);
2178 beforeefficiency->SetMarkerColor(kBlue);
2179 beforeefficiency->SetLineColor(kBlue);
2180 afterefficiency->Draw();
2181 beforeefficiency->Draw("same");
2182 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
2183 lega->AddEntry(beforeefficiency,"Before efficiency correction","p");
2184 lega->AddEntry(afterefficiency,"After efficiency correction","p");
2185 lega->Draw("same");
2186 cefficiencye->cd(2);
2187 cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1);
2188 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr);
2189 efficiencyDDproj->SetTitle("");
2190 efficiencyDDproj->SetStats(0);
2191 efficiencyDDproj->SetMarkerStyle(25);
2192 efficiencyDDproj->SetMarkerColor(2);
2193 efficiencyDDproj->Draw();
2194 cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]);
2195 TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr);
2196 efficiencyDDproja->SetTitle("");
2197 efficiencyDDproja->SetStats(0);
2198 efficiencyDDproja->SetMarkerStyle(26);
2199 efficiencyDDproja->SetMarkerColor(4);
2200 efficiencyDDproja->Draw("same");
c2690925 2201 }
2202 }
2203
e17c1f86 2204 if(fWriteToFile) cefficiency->SaveAs("efficiencyCorrected.eps");
3a72645a 2205 }
2206
c04c80e6 2207 return result;
2208
2209}
3a72645a 2210
c04c80e6 2211//__________________________________________________________________________________
c2690925 2212TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum,Int_t i) const {
c04c80e6 2213 //
2214 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2215 // Give the final pt spectrum to be compared
2216 //
2217
c2690925 2218 if(fNEvents[i] > 0) {
2219
a199006c 2220 Int_t ptpr = 0;
c2690925 2221 if(fBeamType==0) ptpr=0;
2222 if(fBeamType==1) ptpr=1;
c04c80e6 2223
c2690925 2224 TH1D* projection = spectrum->Projection(ptpr);
c04c80e6 2225 CorrectFromTheWidth(projection);
c2690925 2226 TGraphErrors *graphError = NormalizeTH1(projection,i);
c04c80e6 2227 return graphError;
2228
2229 }
2230
2231 return 0x0;
2232
2233
2234}
2235//__________________________________________________________________________________
c2690925 2236TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i) const {
c04c80e6 2237 //
2238 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2239 // Give the final pt spectrum to be compared
2240 //
c2690925 2241 if(fNEvents[i] > 0) {
c04c80e6 2242
a199006c 2243 Int_t ptpr=0;
c2690925 2244 if(fBeamType==0) ptpr=0;
2245 if(fBeamType==1) ptpr=1;
2246
2247 TH1D* projection = (TH1D *) spectrum->Project(ptpr);
c04c80e6 2248 CorrectFromTheWidth(projection);
c2690925 2249 TGraphErrors *graphError = NormalizeTH1(projection,i);
a8ef1999 2250
c04c80e6 2251 return graphError;
2252
2253 }
2254
2255 return 0x0;
2256
2257
2258}
2259//__________________________________________________________________________________
c2690925 2260TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
2261 //
2262 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2263 // Give the final pt spectrum to be compared
2264 //
8c1c76e9 2265 Double_t chargecoefficient = 0.5;
e17c1f86 2266 if(fChargeChoosen != 0) chargecoefficient = 1.0;
8c1c76e9 2267
e17c1f86 2268 Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
8c1c76e9 2269 printf("Normalizing Eta Range %f\n", etarange);
c2690925 2270 if(fNEvents[i] > 0) {
2271
2272 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
2273 Double_t p = 0, dp = 0; Int_t point = 1;
2274 Double_t n = 0, dN = 0;
2275 Double_t nCorr = 0, dNcorr = 0;
2276 Double_t errdN = 0, errdp = 0;
2277 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
2278 point = ibin - input->GetXaxis()->GetFirst();
2279 p = input->GetXaxis()->GetBinCenter(ibin);
2280 dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
2281 n = input->GetBinContent(ibin);
2282 dN = input->GetBinError(ibin);
c2690925 2283 // New point
8c1c76e9 2284 nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
c2690925 2285 errdN = 1./(2. * TMath::Pi() * p);
2286 errdp = 1./(2. * TMath::Pi() * p*p) * n;
8c1c76e9 2287 dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
c2690925 2288
2289 spectrumNormalized->SetPoint(point, p, nCorr);
2290 spectrumNormalized->SetPointError(point, dp, dNcorr);
2291 }
2292 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2293 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2294 spectrumNormalized->SetMarkerStyle(22);
2295 spectrumNormalized->SetMarkerColor(kBlue);
2296 spectrumNormalized->SetLineColor(kBlue);
a8ef1999 2297
c2690925 2298 return spectrumNormalized;
2299
2300 }
2301
2302 return 0x0;
2303
2304
2305}
2306//__________________________________________________________________________________
2307TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) const {
c04c80e6 2308 //
2309 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2310 // Give the final pt spectrum to be compared
2311 //
8c1c76e9 2312 Double_t chargecoefficient = 0.5;
e17c1f86 2313 if(fChargeChoosen != kAllCharge) chargecoefficient = 1.0;
8c1c76e9 2314
e17c1f86 2315 Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
8c1c76e9 2316 printf("Normalizing Eta Range %f\n", etarange);
c2690925 2317 if(normalization > 0) {
c04c80e6 2318
2319 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
2320 Double_t p = 0, dp = 0; Int_t point = 1;
2321 Double_t n = 0, dN = 0;
2322 Double_t nCorr = 0, dNcorr = 0;
2323 Double_t errdN = 0, errdp = 0;
2324 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
2325 point = ibin - input->GetXaxis()->GetFirst();
2326 p = input->GetXaxis()->GetBinCenter(ibin);
2327 dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
2328 n = input->GetBinContent(ibin);
2329 dN = input->GetBinError(ibin);
c04c80e6 2330 // New point
8c1c76e9 2331 nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
c04c80e6 2332 errdN = 1./(2. * TMath::Pi() * p);
2333 errdp = 1./(2. * TMath::Pi() * p*p) * n;
8c1c76e9 2334 dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
c04c80e6 2335
2336 spectrumNormalized->SetPoint(point, p, nCorr);
2337 spectrumNormalized->SetPointError(point, dp, dNcorr);
2338 }
2339 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2340 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2341 spectrumNormalized->SetMarkerStyle(22);
2342 spectrumNormalized->SetMarkerColor(kBlue);
2343 spectrumNormalized->SetLineColor(kBlue);
2344
2345 return spectrumNormalized;
2346
2347 }
2348
2349 return 0x0;
2350
2351
2352}
2353//____________________________________________________________
2354void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
2355 //
2356 // Set the container for a given step to the
2357 //
e17c1f86 2358 if(!fCFContainers) fCFContainers = new TObjArray(kDataContainerV0+1);
c04c80e6 2359 fCFContainers->AddAt(cont, type);
2360}
2361
2362//____________________________________________________________
2363AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
2364 //
2365 // Get Correction Framework Container for given type
2366 //
2367 if(!fCFContainers) return NULL;
2368 return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
2369}
c04c80e6 2370//____________________________________________________________
a8ef1999 2371AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centrality) {
c04c80e6 2372 //
3a72645a 2373 // Slice bin for a given source of electron
c2690925 2374 // nDim is the number of dimension the corrections are done
2375 // dimensions are the definition of the dimensions
2376 // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
2377 // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
a8ef1999 2378 // centrality (-1 means we do not cut on centrality)
c04c80e6 2379 //
2380
2381 Double_t *varMin = new Double_t[container->GetNVar()],
2382 *varMax = new Double_t[container->GetNVar()];
2383
2384 Double_t *binLimits;
2385 for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
2386
2387 binLimits = new Double_t[container->GetNBins(ivar)+1];
2388 container->GetBinLimits(ivar,binLimits);
c2690925 2389 varMin[ivar] = binLimits[0];
2390 varMax[ivar] = binLimits[container->GetNBins(ivar)];
2391 // source
2392 if(ivar == 4){
2393 if((source>= 0) && (source<container->GetNBins(ivar))) {
e17c1f86 2394 varMin[ivar] = binLimits[source];
2395 varMax[ivar] = binLimits[source];
c2690925 2396 }
3a72645a 2397 }
c2690925 2398 // charge
2399 if(ivar == 3) {
e17c1f86 2400 if(charge != kAllCharge) varMin[ivar] = varMax[ivar] = charge;
3a72645a 2401 }
8c1c76e9 2402 // eta
2403 if(ivar == 1){
e17c1f86 2404 for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++)
2405 AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic)));
8c1c76e9 2406 if(fEtaSelected){
2407 varMin[ivar] = fEtaRange[0];
2408 varMax[ivar] = fEtaRange[1];
2409 }
2410 }
e17c1f86 2411 if(fEtaSelected){
2412 fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0]));
2413 fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1]));
2414 AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[0]));
2415 }
a8ef1999 2416 if(ivar == 5){
2417 if((centrality>= 0) && (centrality<container->GetNBins(ivar))) {
2418 varMin[ivar] = binLimits[centrality];
2419 varMax[ivar] = binLimits[centrality];
2420 }
2421 }
3a72645a 2422
11ff28c5 2423
c04c80e6 2424 delete[] binLimits;
2425 }
2426
2427 AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
2428 AddTemporaryObject(k);
2429 delete[] varMin; delete[] varMax;
2430
2431 return k;
2432
2433}
2434
2435//_________________________________________________________________________
3a72645a 2436THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions) const {
c04c80e6 2437 //
3a72645a 2438 // Slice correlation
c04c80e6 2439 //
2440
3a72645a 2441 Int_t ndimensions = correlationmatrix->GetNdimensions();
c2690925 2442 //printf("Number of dimension %d correlation map\n",ndimensions);
c04c80e6 2443 if(ndimensions < (2*nDim)) {
2444 AliError("Problem in the dimensions");
2445 return NULL;
2446 }
2447 Int_t ndimensionsContainer = (Int_t) ndimensions/2;
c2690925 2448 //printf("Number of dimension %d container\n",ndimensionsContainer);
c04c80e6 2449
2450 Int_t *dim = new Int_t[nDim*2];
2451 for(Int_t iter=0; iter < nDim; iter++){
2452 dim[iter] = dimensions[iter];
2453 dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
c2690925 2454 //printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
c04c80e6 2455 }
2456
3a72645a 2457 THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
c04c80e6 2458
2459 delete[] dim;
2460 return k;
2461
2462}
2463//___________________________________________________________________________
2464void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
2465 //
2466 // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
2467 //
2468
2469 TAxis *axis = h1->GetXaxis();
2470 Int_t nbinX = h1->GetNbinsX();
2471
2472 for(Int_t i = 1; i <= nbinX; i++) {
2473
2474 Double_t width = axis->GetBinWidth(i);
2475 Double_t content = h1->GetBinContent(i);
2476 Double_t error = h1->GetBinError(i);
2477 h1->SetBinContent(i,content/width);
2478 h1->SetBinError(i,error/width);
2479 }
2480
2481}
8c1c76e9 2482
2483//___________________________________________________________________________
2484void AliHFEspectrum::CorrectStatErr(AliCFDataGrid *backgroundGrid) const {
2485 //
2486 // Correct statistical error
2487 //
2488
2489 TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
2490 Int_t nbinX = h1->GetNbinsX();
2491 Int_t bins[1];
2492 for(Long_t i = 1; i <= nbinX; i++) {
2493 bins[0] = i;
2494 Float_t content = h1->GetBinContent(i);
2495 if(content>0){
2496 Float_t error = TMath::Sqrt(content);
2497 backgroundGrid->SetElementError(bins, error);
2498 }
2499 }
2500}
2501
c04c80e6 2502//____________________________________________________________
2503void AliHFEspectrum::AddTemporaryObject(TObject *o){
2504 //
2505 // Emulate garbage collection on class level
2506 //
2507 if(!fTemporaryObjects) {
2508 fTemporaryObjects= new TList;
2509 fTemporaryObjects->SetOwner();
2510 }
2511 fTemporaryObjects->Add(o);
2512}
2513
2514//____________________________________________________________
2515void AliHFEspectrum::ClearObject(TObject *o){
2516 //
2517 // Do a safe deletion
2518 //
2519 if(fTemporaryObjects){
2520 if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
2521 delete o;
2522 } else{
2523 // just delete
2524 delete o;
2525 }
2526}
2527//_________________________________________________________________________
c2690925 2528TObject* AliHFEspectrum::GetSpectrum(const AliCFContainer * const c, Int_t step) {
245877d0 2529 AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
c04c80e6 2530 return data;
2531}
2532//_________________________________________________________________________
c2690925 2533TObject* AliHFEspectrum::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
c04c80e6 2534 //
2535 // Create efficiency grid and calculate efficiency
2536 // of step to step0
2537 //
2538 TString name("eff");
2539 name += step;
2540 name+= step0;
2541 AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
2542 eff->CalculateEfficiency(step,step0);
2543 return eff;
2544}
c2690925 2545
2546//____________________________________________________________________________
8c1c76e9 2547THnSparse* AliHFEspectrum::GetCharmWeights(){
c2690925 2548
2549 //
2550 // Measured D->e based weighting factors
2551 //
2552
a8ef1999 2553 const Int_t nDimpp=1;
2554 Int_t nBinpp[nDimpp] = {35};
e17c1f86 2555 Double_t ptbinning1[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
a8ef1999 2556 const Int_t nDimPbPb=2;
2557 Int_t nBinPbPb[nDimPbPb] = {11,35};
2558 Double_t kCentralityRange[12] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
2559 Int_t loopcentr=1;
2560 Int_t looppt=nBinpp[0];
2561 if(fBeamType==0)
2562 {
2563 fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
2564 fWeightCharm->SetBinEdges(0, ptbinning1);
2565 }
2566 if(fBeamType==1)
2567 {
2568 fWeightCharm = new THnSparseF("weightHisto", "weighting factor; centrality bin; pt[GeV/c]", nDimPbPb, nBinPbPb);
2569 fWeightCharm->SetBinEdges(1, ptbinning1);
2570 fWeightCharm->SetBinEdges(0, kCentralityRange);
2571 loopcentr=nBinPbPb[0];
2572 }
2573
11ff28c5 2574 Double_t weight[35]={0.859260, 0.872552, 0.847475, 0.823631, 0.839386, 0.874024, 0.916755, 0.942801, 0.965856, 0.933905, 0.933414, 0.931936, 0.847826, 0.810902, 0.796608, 0.727002, 0.659227, 0.583610, 0.549956, 0.512633, 0.472254, 0.412364, 0.353191, 0.319145, 0.305412, 0.290334, 0.269863, 0.254646, 0.230245, 0.200859, 0.275953, 0.276271, 0.227332, 0.197004, 0.474385};
2575
2576 // Weighting factor for PbPb
2577 //Double_t weight[35]={0.645016, 0.643397, 0.617149, 0.641176, 0.752285, 0.838097, 0.996226, 1.152757, 1.243257, 1.441421, 1.526796, 1.755632, 1.731234, 1.900269, 1.859628, 1.945138, 1.943707, 1.739451, 1.640120, 1.318674, 1.041654, 0.826493, 0.704605, 0.662040, 0.572361, 0.644030, 0.479850, 0.559513, 0.504044, 0.449495, 0.227605, 0.186836, 0.038681, 0.000000, 0.000000};
c2690925 2578
c2690925 2579 //points
8c1c76e9 2580 Double_t pt[1];
a8ef1999 2581 Double_t contents[2];
2582
2583 for(int icentr=0; icentr<loopcentr; icentr++)
2584 {
2585 for(int i=0; i<looppt; i++){
2586 pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
2587 if(fBeamType==1)
2588 {
2589 contents[0]=icentr;
2590 contents[1]=pt[0];
2591 }
2592 if(fBeamType==0)
2593 {
2594 contents[0]=pt[0];
2595 }
2596
2597 fWeightCharm->Fill(contents,weight[i]);
2598 }
2599 }
2600
2601 Int_t nDimSparse = fWeightCharm->GetNdimensions();
2602 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
2603 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
2604
2605 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2606 binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
2607 nBins *= binsvar[iVar];
2608 }
2609
2610 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
2611 // loop that sets 0 error in each bin
2612 for (Long_t iBin=0; iBin<nBins; iBin++) {
2613 Long_t bintmp = iBin ;
2614 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2615 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
2616 bintmp /= binsvar[iVar] ;
2617 }
2618 fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
c2690925 2619 }
c2690925 2620
dcef324e 2621 delete[] binsvar;
2622 delete[] binfill;
2623
c2690925 2624 return fWeightCharm;
2625}
8c1c76e9 2626
2627//____________________________________________________________________________
11ff28c5 2628void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
a8ef1999 2629
2630 // TOF PID efficiencies
2631 Int_t ptpr;
2632 if(fBeamType==0) ptpr=0;
2633 if(fBeamType==1) ptpr=1;
2634
2635 Int_t loopcentr=1;
2636 const Int_t nCentralitybinning=11; //number of centrality bins
2637 if(fBeamType==1)
2638 {
2639 loopcentr=nCentralitybinning;
2640 }
8c1c76e9 2641
11ff28c5 2642 TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
2643 TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
a8ef1999 2644 TF1 *fipfit2 = new TF1("fipfit2","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,10.0);
a8ef1999 2645
2646 TCanvas * cefficiencyParam = new TCanvas("efficiencyParam","efficiencyParam",1200,600);
2647 cefficiencyParam->Divide(2,1);
2648 cefficiencyParam->cd(1);
2649
2650 AliCFContainer *mccontainermcD = 0x0;
2651 TH1D* efficiencysigTOFPIDD[nCentralitybinning];
2652 TH1D* efficiencyTOFPIDD[nCentralitybinning];
11ff28c5 2653 TH1D* efficiencysigesdTOFPIDD[nCentralitybinning];
2654 TH1D* efficiencyesdTOFPIDD[nCentralitybinning];
a8ef1999 2655 Int_t source = -1; //get parameterized TOF PID efficiencies
2656
2657 for(int icentr=0; icentr<loopcentr; icentr++) {
2658 // signal sample
2659 if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 2660 else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
a8ef1999 2661 AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
2662 efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2663
2664 // mb sample for double check
2665 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 2666 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
a8ef1999 2667 AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
2668 efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2669
11ff28c5 2670 // mb sample with reconstructed variables
2671 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
2672 else mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2673 AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD);
2674 efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2675
2676 // mb sample with reconstructed variables
2677 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
2678 else mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2679 AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
2680 efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2681
2682 //fill histo
a8ef1999 2683 efficiencysigTOFPIDD[icentr] = (TH1D *) efficiencymcsigParamTOFPID->Project(ptpr);
2684 efficiencyTOFPIDD[icentr] = (TH1D *) efficiencymcParamTOFPID->Project(ptpr);
11ff28c5 2685 efficiencysigesdTOFPIDD[icentr] = (TH1D *) efficiencysigesdParamTOFPID->Project(ptpr);
2686 efficiencyesdTOFPIDD[icentr] = (TH1D *) efficiencyesdParamTOFPID->Project(ptpr);
a8ef1999 2687 efficiencysigTOFPIDD[icentr]->SetName(Form("efficiencysigTOFPIDD%d",icentr));
2688 efficiencyTOFPIDD[icentr]->SetName(Form("efficiencyTOFPIDD%d",icentr));
11ff28c5 2689 efficiencysigesdTOFPIDD[icentr]->SetName(Form("efficiencysigesdTOFPIDD%d",icentr));
2690 efficiencyesdTOFPIDD[icentr]->SetName(Form("efficiencyesdTOFPIDD%d",icentr));
a8ef1999 2691
11ff28c5 2692 //fit (mc enhenced sample)
a8ef1999 2693 fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
a8ef1999 2694 efficiencysigTOFPIDD[icentr]->Fit(fittofpid,"R");
2695 efficiencysigTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
2696 fEfficiencyTOFPIDD[icentr] = efficiencysigTOFPIDD[icentr]->GetFunction("fittofpid");
11ff28c5 2697
2698 //fit (esd enhenced sample)
2699 efficiencysigesdTOFPIDD[icentr]->Fit(fittofpid,"R");
2700 efficiencysigesdTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
2701 fEfficiencyesdTOFPIDD[icentr] = efficiencysigesdTOFPIDD[icentr]->GetFunction("fittofpid");
2702
a8ef1999 2703 }
8c1c76e9 2704
a8ef1999 2705 // draw (for PbPb, only 1st bin)
11ff28c5 2706 //sig mc
a8ef1999 2707 efficiencysigTOFPIDD[0]->SetTitle("");
2708 efficiencysigTOFPIDD[0]->SetStats(0);
2709 efficiencysigTOFPIDD[0]->SetMarkerStyle(25);
2710 efficiencysigTOFPIDD[0]->SetMarkerColor(2);
2711 efficiencysigTOFPIDD[0]->SetLineColor(2);
2712 efficiencysigTOFPIDD[0]->Draw();
8c1c76e9 2713
11ff28c5 2714 //mb mc
a8ef1999 2715 efficiencyTOFPIDD[0]->SetTitle("");
2716 efficiencyTOFPIDD[0]->SetStats(0);
2717 efficiencyTOFPIDD[0]->SetMarkerStyle(24);
11ff28c5 2718 efficiencyTOFPIDD[0]->SetMarkerColor(4);
2719 efficiencyTOFPIDD[0]->SetLineColor(4);
a8ef1999 2720 efficiencyTOFPIDD[0]->Draw("same");
2721
11ff28c5 2722 //sig esd
2723 efficiencysigesdTOFPIDD[0]->SetTitle("");
2724 efficiencysigesdTOFPIDD[0]->SetStats(0);
2725 efficiencysigesdTOFPIDD[0]->SetMarkerStyle(25);
2726 efficiencysigesdTOFPIDD[0]->SetMarkerColor(3);
2727 efficiencysigesdTOFPIDD[0]->SetLineColor(3);
2728 efficiencysigesdTOFPIDD[0]->Draw("same");
2729
2730 //mb esd
2731 efficiencyesdTOFPIDD[0]->SetTitle("");
2732 efficiencyesdTOFPIDD[0]->SetStats(0);
2733 efficiencyesdTOFPIDD[0]->SetMarkerStyle(25);
2734 efficiencyesdTOFPIDD[0]->SetMarkerColor(1);
2735 efficiencyesdTOFPIDD[0]->SetLineColor(1);
2736 efficiencyesdTOFPIDD[0]->Draw("same");
2737
2738 //signal mc fit
2739 fEfficiencyTOFPIDD[0]->SetLineColor(2);
2740 fEfficiencyTOFPIDD[0]->Draw("same");
2741 //mb esd fit
2742 fEfficiencyesdTOFPIDD[0]->SetLineColor(3);
2743 fEfficiencyesdTOFPIDD[0]->Draw("same");
a8ef1999 2744
11ff28c5 2745 TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
2746 legtofeff->AddEntry(efficiencysigTOFPIDD[0],"TOF PID Step Efficiency","");
2747 legtofeff->AddEntry(efficiencysigTOFPIDD[0],"vs MC p_{t} for enhenced samples","p");
2748 legtofeff->AddEntry(efficiencyTOFPIDD[0],"vs MC p_{t} for mb samples","p");
2749 legtofeff->AddEntry(efficiencysigesdTOFPIDD[0],"vs esd p_{t} for enhenced samples","p");
2750 legtofeff->AddEntry(efficiencyesdTOFPIDD[0],"vs esd p_{t} for mb samples","p");
2751 legtofeff->Draw("same");
a8ef1999 2752
11ff28c5 2753
2754 cefficiencyParam->cd(2);
a8ef1999 2755 // IP cut efficiencies
2756
2757 for(int icentr=0; icentr<loopcentr; icentr++) {
11ff28c5 2758
2759 AliCFContainer *charmCombined = 0x0;
2760 AliCFContainer *beautyCombined = 0x0;
2761
2762 printf("centrality printing %i \n",icentr);
2763
2764 source = 0; //charm enhenced
a8ef1999 2765 if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 2766 else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
a8ef1999 2767 AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
2768 efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2769
11ff28c5 2770 charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");
2771
2772 source = 1; //beauty enhenced
a8ef1999 2773 if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 2774 else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
a8ef1999 2775 AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
2776 efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2777
11ff28c5 2778 beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined");
2779
2780 source = 0; //charm mb
a8ef1999 2781 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 2782 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
a8ef1999 2783 AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
2784 efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2785
11ff28c5 2786 charmCombined->Add(mccontainermcD);
2787 AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined);
2788 efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
2789
2790 source = 1; //beauty mb
a8ef1999 2791 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 2792 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
a8ef1999 2793 AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
2794 efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2795
11ff28c5 2796 beautyCombined->Add(mccontainermcD);
2797 AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined);
2798 efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
2799
2800 source = 2; //conversion mb
a8ef1999 2801 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 2802 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
a8ef1999 2803 AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
2804 efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2805
11ff28c5 2806 source = 3; //non HFE except for the conversion mb
a8ef1999 2807 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 2808 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
a8ef1999 2809 AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
2810 efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2811
11ff28c5 2812 if(fIPEffCombinedSamples){
2813 fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmCombined->Project(ptpr); //signal enhenced + mb
2814 fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautyCombined->Project(ptpr); //signal enhenced + mb
2815 }
2816 else{
2817 fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmSig->Project(ptpr); //signal enhenced only
2818 fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautySig->Project(ptpr); //signal enhenced only
2819 }
2820 fCharmEff[icentr] = (TH1D*)efficiencyCharm->Project(ptpr); //mb only
2821 fBeautyEff[icentr] = (TH1D*)efficiencyBeauty->Project(ptpr); //mb only
2822 fConversionEff[icentr] = (TH1D*)efficiencyConv->Project(ptpr); //mb only
2823 fNonHFEEff[icentr] = (TH1D*)efficiencyNonhfe->Project(ptpr); //mb only
2824
a8ef1999 2825 }
2826
2827 for(int icentr=0; icentr<loopcentr; icentr++) {
2828 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
2829 fipfit->SetLineColor(2);
2830 fEfficiencyBeautySigD[icentr]->Fit(fipfit,"R");
2831 fEfficiencyBeautySigD[icentr]->GetYaxis()->SetTitle("Efficiency");
11ff28c5 2832 if(fBeauty2ndMethod)fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[0]->GetFunction("fipfit"); //why do we need this line?
a8ef1999 2833 else fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[icentr]->GetFunction("fipfit");
2834
11ff28c5 2835 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
2836 fipfit->SetLineColor(1);
2837 fEfficiencyCharmSigD[icentr]->Fit(fipfit,"R");
2838 fEfficiencyCharmSigD[icentr]->GetYaxis()->SetTitle("Efficiency");
2839 fEfficiencyIPCharmD[icentr] = fEfficiencyCharmSigD[icentr]->GetFunction("fipfit");
a8ef1999 2840
11ff28c5 2841 if(fIPParameterizedEff){
a8ef1999 2842 fipfit2->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
2843 fipfit2->SetLineColor(3);
2844 fConversionEff[icentr]->Fit(fipfit2,"R");
2845 fConversionEff[icentr]->GetYaxis()->SetTitle("Efficiency");
2846 fEfficiencyIPConversionD[icentr] = fConversionEff[icentr]->GetFunction("fipfit2");
2847
2848 fipfit2->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
2849 fipfit2->SetLineColor(4);
2850 fNonHFEEff[icentr]->Fit(fipfit2,"R");
2851 fNonHFEEff[icentr]->GetYaxis()->SetTitle("Efficiency");
2852 fEfficiencyIPNonhfeD[icentr] = fNonHFEEff[icentr]->GetFunction("fipfit2");
2853 }
2854 }
2855
2856 // draw (for PbPb, only 1st bin)
a8ef1999 2857 fEfficiencyCharmSigD[0]->SetMarkerStyle(21);
2858 fEfficiencyCharmSigD[0]->SetMarkerColor(1);
2859 fEfficiencyCharmSigD[0]->SetLineColor(1);
2860 fEfficiencyBeautySigD[0]->SetMarkerStyle(21);
2861 fEfficiencyBeautySigD[0]->SetMarkerColor(2);
2862 fEfficiencyBeautySigD[0]->SetLineColor(2);
2863 fCharmEff[0]->SetMarkerStyle(24);
2864 fCharmEff[0]->SetMarkerColor(1);
2865 fCharmEff[0]->SetLineColor(1);
2866 fBeautyEff[0]->SetMarkerStyle(24);
2867 fBeautyEff[0]->SetMarkerColor(2);
2868 fBeautyEff[0]->SetLineColor(2);
2869 fConversionEff[0]->SetMarkerStyle(24);
2870 fConversionEff[0]->SetMarkerColor(3);
2871 fConversionEff[0]->SetLineColor(3);
2872 fNonHFEEff[0]->SetMarkerStyle(24);
2873 fNonHFEEff[0]->SetMarkerColor(4);
2874 fNonHFEEff[0]->SetLineColor(4);
2875
2876 fEfficiencyCharmSigD[0]->Draw();
2877 fEfficiencyBeautySigD[0]->Draw("same");
11ff28c5 2878 //fCharmEff[0]->Draw("same");
2879 //fBeautyEff[0]->Draw("same");
a8ef1999 2880 fConversionEff[0]->Draw("same");
2881 fNonHFEEff[0]->Draw("same");
2882
2883 fEfficiencyIPBeautyD[0]->Draw("same");
11ff28c5 2884 fEfficiencyIPCharmD[0]->Draw("same");
a8ef1999 2885 if(fIPParameterizedEff){
a8ef1999 2886 fEfficiencyIPConversionD[0]->Draw("same");
2887 fEfficiencyIPNonhfeD[0]->Draw("same");
2888 }
11ff28c5 2889 TLegend *legipeff = new TLegend(0.55,0.2,0.85,0.39);
2890 legipeff->AddEntry(fEfficiencyBeautySigD[0],"IP Step Efficiency","");
2891 legipeff->AddEntry(fEfficiencyBeautySigD[0],"beauty e","p");
2892 legipeff->AddEntry(fEfficiencyCharmSigD[0],"charm e","p");
2893 legipeff->AddEntry(fConversionEff[0],"conversion e","p");
2894 legipeff->AddEntry(fNonHFEEff[0],"Dalitz e","p");
2895 legipeff->Draw("same");
8c1c76e9 2896}
2897
2898//____________________________________________________________________________
a8ef1999 2899THnSparse* AliHFEspectrum::GetBeautyIPEff(){
8c1c76e9 2900 //
a8ef1999 2901 // Return beauty electron IP cut efficiency
8c1c76e9 2902 //
2903
a8ef1999 2904 const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
2905 const Int_t nCentralitybinning=11;//number of centrality bins
2906 Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
2907 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
2908 Int_t ptpr = 0;
2909 Int_t nDim=1; //dimensions of the efficiency weighting grid
2910 if(fBeamType==1)
2911 {
2912 ptpr=1;
2913 nDim=2; //dimensions of the efficiency weighting grid
2914 }
2915 Int_t nBin[1] = {nPtbinning1};
2916 Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
8c1c76e9 2917
8c1c76e9 2918
a8ef1999 2919 THnSparseF *ipcut;
2920 if(fBeamType==0) ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
2921 else ipcut = new THnSparseF("beff", "b IP efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
2922
8c1c76e9 2923 for(Int_t idim = 0; idim < nDim; idim++)
a8ef1999 2924 {
2925 if(nDim==1) ipcut->SetBinEdges(idim, kPtRange);
2926 if(nDim==2)
2927 {
2928 ipcut->SetBinEdges(0, kCentralityRange);
2929 ipcut->SetBinEdges(1, kPtRange);
2930 }
2931 }
8c1c76e9 2932 Double_t pt[1];
a8ef1999 2933 Double_t weight[nCentralitybinning];
2934 Double_t contents[2];
2935
2936 for(Int_t a=0;a<11;a++)
2937 {
2938 weight[a] = 1.0;
2939 }
2940
2941
2942 Int_t looppt=nBin[0];
2943 Int_t loopcentr=1;
2944 if(fBeamType==1)
2945 {
2946 loopcentr=nBinPbPb[0];
2947 }
2948
2949
2950 for(int icentr=0; icentr<loopcentr; icentr++)
2951 {
2952 for(int i=0; i<looppt; i++)
2953 {
2954 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
2955 if(fEfficiencyIPBeautyD[icentr])
2956 weight[icentr]=fEfficiencyIPBeautyD[icentr]->Eval(pt[0]);
2957 else{
2958 printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
2959 weight[icentr] = fEfficiencyBeautySigD[icentr]->GetBinContent(i+1);
2960 }
2961
2962 if(fBeamType==1){
2963 contents[0]=icentr;
2964 contents[1]=pt[0];
2965 }
2966 if(fBeamType==0){
2967 contents[0]=pt[0];
2968 }
2969 ipcut->Fill(contents,weight[icentr]);
2970 }
2971 }
2972
2973
2974 Int_t nDimSparse = ipcut->GetNdimensions();
2975 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
2976 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
2977
2978 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2979 binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
2980 nBins *= binsvar[iVar];
2981 }
2982
2983 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
2984 // loop that sets 0 error in each bin
2985 for (Long_t iBin=0; iBin<nBins; iBin++) {
2986 Long_t bintmp = iBin ;
2987 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2988 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
2989 bintmp /= binsvar[iVar] ;
2990 }
2991 ipcut->SetBinError(binfill,0.); // put 0 everywhere
2992 }
8c1c76e9 2993
dcef324e 2994 delete[] binsvar;
2995 delete[] binfill;
2996
8c1c76e9 2997 return ipcut;
2998}
2999
3000//____________________________________________________________________________
3001THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
3002 //
3003 // Return PID x IP cut efficiency
3004 //
a8ef1999 3005 const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
3006 const Int_t nCentralitybinning=11;//number of centrality bins
3007 Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
3008 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3009 Int_t ptpr = 0;
3010 Int_t nDim=1; //dimensions of the efficiency weighting grid
3011 if(fBeamType==1)
3012 {
3013 ptpr=1;
3014 nDim=2; //dimensions of the efficiency weighting grid
3015 }
3016 Int_t nBin[1] = {nPtbinning1};
3017 Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3018
3019 THnSparseF *pideff;
3020 if(fBeamType==0) pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
3021 else pideff = new THnSparseF("pideff", "PID efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
3022 for(Int_t idim = 0; idim < nDim; idim++)
3023 {
3024
3025 if(nDim==1) pideff->SetBinEdges(idim, kPtRange);
3026 if(nDim==2)
3027 {
3028 pideff->SetBinEdges(0, kCentralityRange);
3029 pideff->SetBinEdges(1, kPtRange);
3030 }
3031 }
8c1c76e9 3032
a8ef1999 3033 Double_t pt[1];
3034 Double_t weight[nCentralitybinning];
11ff28c5 3035 Double_t weightErr[nCentralitybinning];
a8ef1999 3036 Double_t contents[2];
3037
3024f297 3038 for(Int_t a=0;a<nCentralitybinning;a++)
a8ef1999 3039 {
3040 weight[a] = 1.0;
3024f297 3041 weightErr[a] = 1.0;
3042
a8ef1999 3043 }
3044
3045 Int_t looppt=nBin[0];
3046 Int_t loopcentr=1;
11ff28c5 3047 Int_t ibin[2];
a8ef1999 3048 if(fBeamType==1)
3049 {
3050 loopcentr=nBinPbPb[0];
3051 }
3052
3053 for(int icentr=0; icentr<loopcentr; icentr++)
3054 {
3055 Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
3056 for(int i=0; i<looppt; i++)
3057 {
3058 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3059
3060 Double_t tofpideff = 0.;
11ff28c5 3061 Double_t tofpideffesd = 0.;
a8ef1999 3062 if(fEfficiencyTOFPIDD[icentr])
11ff28c5 3063 tofpideff = fEfficiencyTOFPIDD[icentr]->Eval(pt[0]);
a8ef1999 3064 else{
3065 printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
3066 }
11ff28c5 3067 if(fEfficiencyesdTOFPIDD[icentr])
3068 tofpideffesd = fEfficiencyesdTOFPIDD[icentr]->Eval(pt[0]);
3069 else{
3070 printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
3071 }
a8ef1999 3072
3073 //tof pid eff x tpc pid eff x ip cut eff
3074 if(fIPParameterizedEff){
3075 if(source==0) {
11ff28c5 3076 if(fEfficiencyIPCharmD[icentr]){
a8ef1999 3077 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
11ff28c5 3078 weightErr[icentr] = 0;
3079 }
a8ef1999 3080 else{
3081 printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
3082 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
11ff28c5 3083 weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
a8ef1999 3084 }
3085 }
3086 else if(source==2) {
11ff28c5 3087 if(fEfficiencyIPConversionD[icentr]){
3088 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD[icentr]->Eval(pt[0]);
3089 weightErr[icentr] = 0;
3090 }
a8ef1999 3091 else{
3092 printf("Fit failed on conversion IP cut efficiency for centrality %d\n",icentr);
11ff28c5 3093 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1);
3094 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
a8ef1999 3095 }
3096 }
3097 else if(source==3) {
11ff28c5 3098 if(fEfficiencyIPNonhfeD[icentr]){
3099 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD[icentr]->Eval(pt[0]);
3100 weightErr[icentr] = 0;
3101 }
a8ef1999 3102 else{
3103 printf("Fit failed on dalitz IP cut efficiency for centrality %d\n",icentr);
11ff28c5 3104 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1);
3105 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
a8ef1999 3106 }
3107 }
3108 }
3109 else{
11ff28c5 3110 if(source==0){
3111 if(fEfficiencyIPCharmD[icentr]){
3112 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
3113 weightErr[icentr] = 0;
3114 }
3115 else{
3116 printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
3117 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
3118 weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
3119 }
3120 }
3121 else if(source==2){
3122 if(fBeamType==0){
3123 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
3124 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
3125 }
3126 else{
3127 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); // conversion
3128 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
3129 }
3130 }
3131 else if(source==3){
3132 if(fBeamType==0){
3133 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
3134 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
3135 }
3136 else{
3137 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); // Dalitz
3138 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
3139 }
3140 }
a8ef1999 3141 }
3142
3143 if(fBeamType==1){
3144 contents[0]=icentr;
3145 contents[1]=pt[0];
11ff28c5 3146 ibin[0]=icentr;
3147 ibin[1]=i+1;
a8ef1999 3148 }
3149 if(fBeamType==0){
3150 contents[0]=pt[0];
11ff28c5 3151 ibin[0]=i+1;
a8ef1999 3152 }
3153
3154 pideff->Fill(contents,weight[icentr]);
11ff28c5 3155 pideff->SetBinError(ibin,weightErr[icentr]);
a8ef1999 3156 }
3157 }
3158
3159 Int_t nDimSparse = pideff->GetNdimensions();
3160 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3161 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3162
3163 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3164 binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
3165 nBins *= binsvar[iVar];
3166 }
3167
3168 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3169 // loop that sets 0 error in each bin
3170 for (Long_t iBin=0; iBin<nBins; iBin++) {
3171 Long_t bintmp = iBin ;
3172 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3173 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3174 bintmp /= binsvar[iVar] ;
3175 }
a8ef1999 3176 }
8c1c76e9 3177
dcef324e 3178 delete[] binsvar;
3179 delete[] binfill;
8c1c76e9 3180
11ff28c5 3181
8c1c76e9 3182 return pideff;
3183}
a8ef1999 3184
3185//__________________________________________________________________________
3186AliCFDataGrid *AliHFEspectrum::GetRawBspectra2ndMethod(){
3187 //
3188 // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method
3189 //
3190 Int_t ptpr = 0;
3191 Int_t nDim = 1;
3192 if(fBeamType==0)
3193 {
3194 ptpr=0;
3195 }
3196 if(fBeamType==1)
3197 {
3198 ptpr=1;
3199 nDim=2;
3200 }
3201
11ff28c5 3202 // const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
3203 const Int_t nPtbinning1 = 18;//number of pt bins, according to new binning
a8ef1999 3204 const Int_t nCentralitybinning=11;//number of centrality bins
11ff28c5 3205 // Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
3206 Double_t kPtRange[19] = {0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 2, 2.5, 3, 4, 5, 6, 8, 12, 16, 20};
3207
a8ef1999 3208 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3209 Int_t nBin[1] = {nPtbinning1};
3210 Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3211 //fBSpectrum2ndMethod->GetNbinsX()
3212
3213 AliCFDataGrid *rawBeautyContainer;
3214 if(fBeamType==0) rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
3215 else rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBinPbPb);
3216 // printf("number of bins= %d\n",bins[0]);
3217
3218
3219
3220
3221 THnSparseF *brawspectra;
3222 if(fBeamType==0) brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin);
3223 else brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBinPbPb);
3224 if(fBeamType==0) brawspectra->SetBinEdges(0, kPtRange);
3225 if(fBeamType==1)
3226 {
3227 // brawspectra->SetBinEdges(0, centralityBins);
3228 brawspectra->SetBinEdges(0, kCentralityRange);
3229 brawspectra->SetBinEdges(1, kPtRange);
3230 }
3231
3232 Double_t pt[1];
3233 Double_t yields= 0.;
3234 Double_t valuesb[2];
3235
3236 //Int_t looppt=nBin[0];
3237 Int_t loopcentr=1;
3238 if(fBeamType==1)
3239 {
3240 loopcentr=nBinPbPb[0];
3241 }
3242
3243 for(int icentr=0; icentr<loopcentr; icentr++)
3244 {
3245
3246 for(int i=0; i<fBSpectrum2ndMethod->GetNbinsX(); i++){
3247 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3248
3249 yields = fBSpectrum2ndMethod->GetBinContent(i+1);
3250
3251 if(fBeamType==1)
3252 {
3253 valuesb[0]=icentr;
3254 valuesb[1]=pt[0];
3255 }
3256 if(fBeamType==0) valuesb[0]=pt[0];
3257 brawspectra->Fill(valuesb,yields);
3258 }
3259 }
3260
3261
3262
3263 Int_t nDimSparse = brawspectra->GetNdimensions();
3264 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3265 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3266
3267 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3268 binsvar[iVar] = brawspectra->GetAxis(iVar)->GetNbins();
3269 nBins *= binsvar[iVar];
3270 }
3271
3272 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3273 // loop that sets 0 error in each bin
3274 for (Long_t iBin=0; iBin<nBins; iBin++) {
3275 Long_t bintmp = iBin ;
3276 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3277 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3278 bintmp /= binsvar[iVar] ;
3279 }
3280 brawspectra->SetBinError(binfill,0.); // put 0 everywhere
3281 }
3282
3283
3284 rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency
3285 TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(ptpr);
3286
3287 new TCanvas;
3288 fBSpectrum2ndMethod->SetMarkerStyle(24);
3289 fBSpectrum2ndMethod->Draw("p");
3290 hRawBeautySpectra->SetMarkerStyle(25);
3291 hRawBeautySpectra->Draw("samep");
dcef324e 3292
11ff28c5 3293 delete[] binfill;
dcef324e 3294 delete[] binsvar;
11ff28c5 3295
a8ef1999 3296 return rawBeautyContainer;
3297}
3298
e17c1f86 3299//__________________________________________________________________________
a8ef1999 3300void AliHFEspectrum::CalculateNonHFEsyst(Int_t centrality){
3301 //
3302 // Calculate non HFE sys
e17c1f86 3303 //
e17c1f86 3304 //
a8ef1999 3305
e17c1f86 3306 if(!fNonHFEsyst)
3307 return;
3308
3309 Double_t evtnorm[1] = {0.0};
dcef324e 3310 if(fNMCbgEvents[0]>0) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
e17c1f86 3311
3312 AliCFDataGrid *convSourceGrid[kElecBgSources][kBgLevels];
3313 AliCFDataGrid *nonHFESourceGrid[kElecBgSources][kBgLevels];
3314
3315 AliCFDataGrid *bgLevelGrid[kBgLevels];
3316 AliCFDataGrid *bgNonHFEGrid[kBgLevels];
3317 AliCFDataGrid *bgConvGrid[kBgLevels];
3318
3319 Int_t stepbackground = 3;
3320 Int_t* bins=new Int_t[1];
a8ef1999 3321
3322 bins[0]=fConversionEff[centrality]->GetNbinsX();
3323
e17c1f86 3324 AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
3325 AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
3326
a8ef1999 3327 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
e17c1f86 3328 for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
a8ef1999 3329 convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d_%d",iSource,iLevel,centrality),Form("convGrid_%d_%d_%d",iSource,iLevel,centrality),*fConvSourceContainer[iSource][iLevel][centrality],stepbackground);
e17c1f86 3330 weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
3331 convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
a8ef1999 3332
3333 nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d_%d",iSource,iLevel,centrality),Form("nonHFEGrid_%d_%d_%d",iSource,iLevel,centrality),*fNonHFESourceContainer[iSource][iLevel][centrality],stepbackground);
e17c1f86 3334 weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
3335 nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
3336 }
a8ef1999 3337
e17c1f86 3338 bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
3339 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
3340 bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
3341 }
3342
3343 bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone();
3344 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){//add other sources to get overall background from all meson decays
3345 bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
3346 }
a8ef1999 3347
e17c1f86 3348 bgLevelGrid[iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
3349 bgLevelGrid[iLevel]->Add(bgNonHFEGrid[iLevel]);
3350 }
3351
a8ef1999 3352
e17c1f86 3353 //Now subtract the mean from upper, and lower from mean container to get the error based on the pion yield uncertainty (-> this error sums linearly, since its contribution to all meson yields is correlated)
3354 AliCFDataGrid *bgErrorGrid[2];
3355 bgErrorGrid[0] = (AliCFDataGrid*)bgLevelGrid[1]->Clone();
3356 bgErrorGrid[1] = (AliCFDataGrid*)bgLevelGrid[2]->Clone();
3357 bgErrorGrid[0]->Add(bgLevelGrid[0],-1.);
3358 bgErrorGrid[1]->Add(bgLevelGrid[0],-1.);
3359
3360 //plot absolute differences between limit yields (upper/lower limit, based on pi0 errors) and best estimate
3361 TH1D* hpiErrors[2];
3362 hpiErrors[0] = (TH1D*)bgErrorGrid[0]->Project(0);
3363 hpiErrors[0]->Scale(-1.);
3364 hpiErrors[0]->SetTitle("Absolute systematic errors from non-HF meson decays and conversions");
3365 hpiErrors[1] = (TH1D*)bgErrorGrid[1]->Project(0);
3366
3367
3368
3369 //Calculate the scaling errors for electrons from all mesons except for pions: square sum of (0.3 * best yield estimate), where 0.3 is the error generally assumed for m_t scaling
3370 TH1D *hSpeciesErrors[kElecBgSources-1];
3371 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
3372 hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0);
3373 TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0);
3374 hSpeciesErrors[iSource-1]->Add(hNonHFEtemp);
3375 hSpeciesErrors[iSource-1]->Scale(0.3);
3376 }
3377
3378 TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[0]->Clone();
3379 TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[0]->Clone();
3380 TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[0]->Clone();
3381
3382 TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
3383 TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
3384
3385 for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
3386 Double_t pi0basedErrLow = hpiErrors[0]->GetBinContent(iBin);
3387 Double_t pi0basedErrUp = hpiErrors[1]->GetBinContent(iBin);
3388
3389 Double_t sqrsumErrs= 0;
3390 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
3391 Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
3392 sqrsumErrs+=(scalingErr*scalingErr);
3393 }
3394 for(Int_t iErr = 0; iErr < 2; iErr++){
3395 hpiErrors[iErr]->SetBinContent(iBin,hpiErrors[iErr]->GetBinContent(iBin)/hpiErrors[iErr]->GetBinWidth(iBin));
3396 }
3397 hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+sqrsumErrs));
3398 hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+sqrsumErrs));
3399 hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin));
3400
3401 hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin));
3402 hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin));
3403 }
3404
3405
3406 // /hOverallSystErrUp->GetBinWidth(iBin))
3407
3408 TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
3409 cPiErrors->cd();
3410 cPiErrors->SetLogx();
3411 cPiErrors->SetLogy();
3412 hpiErrors[0]->Draw();
3413 hpiErrors[1]->SetMarkerColor(kBlack);
3414 hpiErrors[1]->SetLineColor(kBlack);
3415 hpiErrors[1]->Draw("SAME");
3416 hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
3417 hOverallBinScaledErrsUp->SetLineColor(kBlue);
3418 hOverallBinScaledErrsUp->Draw("SAME");
3419 hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
3420 hOverallBinScaledErrsLow->SetLineColor(kGreen);
3421 hOverallBinScaledErrsLow->Draw("SAME");
3422 hScalingErrors->SetLineColor(11);
3423 hScalingErrors->Draw("SAME");
3424
3425 TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
3426 lPiErr->AddEntry(hpiErrors[0],"Lower error from pion error");
3427 lPiErr->AddEntry(hpiErrors[1],"Upper error from pion error");
3428 lPiErr->AddEntry(hScalingErrors, "scaling error");
3429 lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics");
3430 lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
3431 lPiErr->Draw("SAME");
3432
3433 //Normalize errors
3434 TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone();
3435 TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone();
3436 hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!)
3437 hLowSystScaled->Scale(evtnorm[0]);
3438 TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone();
3439 TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone();
3440 //histograms to be normalized to TGraphErrors
3441 CorrectFromTheWidth(hNormAllSystErrUp);
3442 CorrectFromTheWidth(hNormAllSystErrLow);
3443
3444 TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs");
3445 cNormOvErrs->cd();
3446 cNormOvErrs->SetLogx();
3447 cNormOvErrs->SetLogy();
3448
3449 TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp);
3450 TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow);
3451 gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors");
3452 gOverallSystErrUp->SetMarkerColor(kBlack);
3453 gOverallSystErrUp->SetLineColor(kBlack);
3454 gOverallSystErrLow->SetMarkerColor(kRed);
3455 gOverallSystErrLow->SetLineColor(kRed);
3456 gOverallSystErrUp->Draw("AP");
3457 gOverallSystErrLow->Draw("PSAME");
3458 TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89);
3459 lAllSys->AddEntry(gOverallSystErrLow,"lower","p");
3460 lAllSys->AddEntry(gOverallSystErrUp,"upper","p");
3461 lAllSys->Draw("same");
3462
3463
3464 AliCFDataGrid *bgYieldGrid;
3465 bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0]->Clone();
3466
3467 TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0);
3468 TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone();
3469 hRelErrUp->Divide(hBgYield);
3470 TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone();
3471 hRelErrLow->Divide(hBgYield);
3472
3473 TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs");
3474 cRelErrs->cd();
3475 cRelErrs->SetLogx();
3476 hRelErrUp->SetTitle("Relative error of non-HFE background yield");
3477 hRelErrUp->Draw();
3478 hRelErrLow->SetLineColor(kBlack);
3479 hRelErrLow->Draw("SAME");
3480
3481 TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95);
3482 lRel->AddEntry(hRelErrUp, "upper");
3483 lRel->AddEntry(hRelErrLow, "lower");
3484 lRel->Draw("SAME");
3485
3486 //CorrectFromTheWidth(hBgYield);
3487 //hBgYield->Scale(evtnorm[0]);
3488
3489
3490 //write histograms/TGraphs to file
3491 TFile *output = new TFile("systHists.root","recreate");
3492
3493 hBgYield->SetName("hBgYield");
3494 hBgYield->Write();
3495 hRelErrUp->SetName("hRelErrUp");
3496 hRelErrUp->Write();
3497 hRelErrLow->SetName("hRelErrLow");
3498 hRelErrLow->Write();
3499 hUpSystScaled->SetName("hOverallSystErrUp");
3500 hUpSystScaled->Write();
3501 hLowSystScaled->SetName("hOverallSystErrLow");
3502 hLowSystScaled->Write();
3503 gOverallSystErrUp->SetName("gOverallSystErrUp");
3504 gOverallSystErrUp->Write();
3505 gOverallSystErrLow->SetName("gOverallSystErrLow");
3506 gOverallSystErrLow->Write();
3507
3508 output->Close();
a8ef1999 3509 delete output;
3510
e17c1f86 3511}
a8ef1999 3512