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