]> 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
Merge branch 'master' into TPCdev
[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++){
ff8249bd 1688 if(iSource == 1)
1689 backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0],1.41);//correction for the eta Dalitz decay branching ratio in PYTHIA
1690 else
1691 backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0]);
e17c1f86 1692 }
1693 }
1694 else{
1695 // Background Estimate
1696 backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
1697 }
8c1c76e9 1698 if(!backgroundContainer){
1699 AliError("MC background container not found");
1700 return NULL;
c2690925 1701 }
e17c1f86 1702
1703
8c1c76e9 1704 Int_t stepbackground = 3;
a8ef1999 1705
8c1c76e9 1706 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
a8ef1999 1707 Int_t *nBinpp=new Int_t[1];
1708 Int_t* binspp=new Int_t[1];
1709 binspp[0]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1710
1711 Int_t *nBinPbPb=new Int_t[2];
1712 Int_t* binsPbPb=new Int_t[2];
1713 binsPbPb[1]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1714 binsPbPb[0]=12;
1715
1716 Int_t looppt=binspp[0];
1717 if(fBeamType==1) looppt=binsPbPb[1];
1718
1719
1720 for(Long_t iBin=1; iBin<= looppt;iBin++){
1721 if(fBeamType==0)
1722 {
1723 nBinpp[0]=iBin;
1724 backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
1725 backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
1726 }
1727 if(fBeamType==1)
1728 {
1729 for(Long_t iiBin=1; iiBin<=binsPbPb[0];iiBin++){
1730 nBinPbPb[0]=iiBin;
1731 nBinPbPb[1]=iBin;
1732 Double_t evtnormPbPb=0;
11ff28c5 1733 if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
a8ef1999 1734 backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1735 backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1736 }
1737 }
e17c1f86 1738 }
1739 //end of workaround for statistical errors
a8ef1999 1740 AliCFDataGrid *weightedNonHFEContainer;
1741 if(fBeamType==0) weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,binspp);
1742 else weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",2,binsPbPb);
8c1c76e9 1743 weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1744 backgroundGrid->Multiply(weightedNonHFEContainer,1.0);
c2690925 1745
11ff28c5 1746 delete[] nBinpp;
dcef324e 1747 delete[] binspp;
1748 delete[] nBinPbPb;
11ff28c5 1749 delete[] binsPbPb;
dcef324e 1750
8c1c76e9 1751 return backgroundGrid;
c2690925 1752}
1753
1754//____________________________________________________________
1755AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
1756
1757 //
1758 // Apply TPC pid efficiency correction from parametrisation
1759 //
1760
1761 // Data in the right format
1762 AliCFDataGrid *dataGrid = 0x0;
1763 if(bgsubpectrum) {
1764 dataGrid = bgsubpectrum;
1765 }
1766 else {
1767
1768 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1769 if(!dataContainer){
1770 AliError("Data Container not available");
1771 return NULL;
1772 }
c2690925 1773 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1774 }
c2690925 1775 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1776 result->SetName("ParametrizedEfficiencyBefore");
1777 THnSparse *h = result->GetGrid();
1778 Int_t nbdimensions = h->GetNdimensions();
1779 //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
c2690925 1780 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1781 if(!dataContainer){
1782 AliError("Data Container not available");
1783 return NULL;
1784 }
1785 AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
1786 dataContainerbis->Add(dataContainerbis,-1.0);
1787
1788
1789 Int_t* coord = new Int_t[nbdimensions];
1790 memset(coord, 0, sizeof(Int_t) * nbdimensions);
1791 Double_t* points = new Double_t[nbdimensions];
1792
c2690925 1793 ULong64_t nEntries = h->GetNbins();
1794 for (ULong64_t i = 0; i < nEntries; ++i) {
1795
1796 Double_t value = h->GetBinContent(i, coord);
1797 //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
1798 //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
1799
1800 // Get the bin co-ordinates given an coord
1801 for (Int_t j = 0; j < nbdimensions; ++j)
1802 points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1803
1804 if (!fEfficiencyFunction->IsInside(points))
1805 continue;
1806 TF1::RejectPoint(kFALSE);
1807
1808 // Evaulate function at points
1809 Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1810 //printf("Value efficiency is %f\n",valueEfficiency);
1811
1812 if(valueEfficiency > 0.0) {
1813 h->SetBinContent(coord,value/valueEfficiency);
1814 dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1815 }
1816 Double_t error = h->GetBinError(i);
1817 h->SetBinError(coord,error/valueEfficiency);
1818 dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1819
1820
1821 }
1822
6c70d827 1823 delete[] coord;
1824 delete[] points;
1825
c2690925 1826 AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1827
1828 if(fDebugLevel > 0) {
a8ef1999 1829
c2690925 1830 TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700);
1831 cEfficiencyParametrized->Divide(2,1);
1832 cEfficiencyParametrized->cd(1);
1833 TH1D *afterE = (TH1D *) resultt->Project(0);
1834 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1835 CorrectFromTheWidth(afterE);
1836 CorrectFromTheWidth(beforeE);
1837 afterE->SetStats(0);
1838 afterE->SetTitle("");
1839 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1840 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1841 afterE->SetMarkerStyle(25);
1842 afterE->SetMarkerColor(kBlack);
1843 afterE->SetLineColor(kBlack);
1844 beforeE->SetStats(0);
1845 beforeE->SetTitle("");
1846 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1847 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1848 beforeE->SetMarkerStyle(24);
1849 beforeE->SetMarkerColor(kBlue);
1850 beforeE->SetLineColor(kBlue);
1851 gPad->SetLogy();
1852 afterE->Draw();
1853 beforeE->Draw("same");
1854 TLegend *legefficiencyparametrized = new TLegend(0.4,0.6,0.89,0.89);
1855 legefficiencyparametrized->AddEntry(beforeE,"Before Efficiency correction","p");
1856 legefficiencyparametrized->AddEntry(afterE,"After Efficiency correction","p");
1857 legefficiencyparametrized->Draw("same");
1858 cEfficiencyParametrized->cd(2);
1859 fEfficiencyFunction->Draw();
1860 //cEfficiencyParametrized->cd(3);
1861 //TH1D *ratioefficiency = (TH1D *) beforeE->Clone();
1862 //ratioefficiency->Divide(afterE);
1863 //ratioefficiency->Draw();
1864
e17c1f86 1865 if(fWriteToFile) cEfficiencyParametrized->SaveAs("efficiency.eps");
a8ef1999 1866
c2690925 1867 }
1868
c2690925 1869 return resultt;
1870
1871}
c04c80e6 1872//____________________________________________________________
3a72645a 1873AliCFDataGrid *AliHFEspectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
1874
c04c80e6 1875 //
3a72645a 1876 // Apply TPC pid efficiency correction from V0
c04c80e6 1877 //
1878
3a72645a 1879 AliCFContainer *v0Container = GetContainer(kDataContainerV0);
1880 if(!v0Container){
1881 AliError("V0 Container not available");
c04c80e6 1882 return NULL;
1883 }
3a72645a 1884
1885 // Efficiency
1886 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
1887 efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
1888
1889 // Data in the right format
1890 AliCFDataGrid *dataGrid = 0x0;
1891 if(bgsubpectrum) {
1892 dataGrid = bgsubpectrum;
c04c80e6 1893 }
3a72645a 1894 else {
1895
1896 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1897 if(!dataContainer){
1898 AliError("Data Container not available");
1899 return NULL;
1900 }
c04c80e6 1901
3a72645a 1902 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1903 }
c04c80e6 1904
3a72645a 1905 // Correct
1906 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1907 result->ApplyEffCorrection(*efficiencyD);
c04c80e6 1908
3a72645a 1909 if(fDebugLevel > 0) {
c2690925 1910
1911 Int_t ptpr;
1912 if(fBeamType==0) ptpr=0;
1913 if(fBeamType==1) ptpr=1;
3a72645a 1914
1915 TCanvas * cV0Efficiency = new TCanvas("V0Efficiency","V0Efficiency",1000,700);
1916 cV0Efficiency->Divide(2,1);
1917 cV0Efficiency->cd(1);
c2690925 1918 TH1D *afterE = (TH1D *) result->Project(ptpr);
1919 TH1D *beforeE = (TH1D *) dataGrid->Project(ptpr);
3a72645a 1920 afterE->SetStats(0);
1921 afterE->SetTitle("");
1922 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1923 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1924 afterE->SetMarkerStyle(25);
1925 afterE->SetMarkerColor(kBlack);
1926 afterE->SetLineColor(kBlack);
1927 beforeE->SetStats(0);
1928 beforeE->SetTitle("");
1929 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1930 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1931 beforeE->SetMarkerStyle(24);
1932 beforeE->SetMarkerColor(kBlue);
1933 beforeE->SetLineColor(kBlue);
1934 afterE->Draw();
1935 beforeE->Draw("same");
1936 TLegend *legV0efficiency = new TLegend(0.4,0.6,0.89,0.89);
1937 legV0efficiency->AddEntry(beforeE,"Before Efficiency correction","p");
1938 legV0efficiency->AddEntry(afterE,"After Efficiency correction","p");
1939 legV0efficiency->Draw("same");
1940 cV0Efficiency->cd(2);
c2690925 1941 TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(ptpr);
3a72645a 1942 efficiencyDproj->SetTitle("");
1943 efficiencyDproj->SetStats(0);
1944 efficiencyDproj->SetMarkerStyle(25);
1945 efficiencyDproj->Draw();
c04c80e6 1946
c2690925 1947 if(fBeamType==1) {
1948
1949 TCanvas * cV0Efficiencye = new TCanvas("V0Efficiency_allspectra","V0Efficiency_allspectra",1000,700);
1950 cV0Efficiencye->Divide(2,1);
1951 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1952 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1953
1954 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1955 TAxis *cenaxisa = sparseafter->GetAxis(0);
1956 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1957 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1958 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1959 TAxis *cenaxisc = efficiencya->GetAxis(0);
1960 Int_t nbbin = cenaxisb->GetNbins();
1961 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1962 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1963 for(Int_t binc = 0; binc < nbbin; binc++){
e17c1f86 1964 TString titlee("V0Efficiency_centrality_bin_");
1965 titlee += binc;
1966 TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1967 ccV0Efficiency->Divide(2,1);
1968 ccV0Efficiency->cd(1);
1969 gPad->SetLogy();
1970 cenaxisa->SetRange(binc+1,binc+1);
1971 cenaxisb->SetRange(binc+1,binc+1);
1972 cenaxisc->SetRange(binc+1,binc+1);
1973 TH1D *aftere = (TH1D *) sparseafter->Projection(1);
1974 TH1D *beforee = (TH1D *) sparsebefore->Projection(1);
1975 CorrectFromTheWidth(aftere);
1976 CorrectFromTheWidth(beforee);
1977 aftere->SetStats(0);
1978 aftere->SetTitle((const char*)titlee);
1979 aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1980 aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1981 aftere->SetMarkerStyle(25);
1982 aftere->SetMarkerColor(kBlack);
1983 aftere->SetLineColor(kBlack);
1984 beforee->SetStats(0);
1985 beforee->SetTitle((const char*)titlee);
1986 beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1987 beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1988 beforee->SetMarkerStyle(24);
1989 beforee->SetMarkerColor(kBlue);
1990 beforee->SetLineColor(kBlue);
1991 aftere->Draw();
1992 beforee->Draw("same");
1993 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1994 lega->AddEntry(beforee,"Before correction","p");
1995 lega->AddEntry(aftere,"After correction","p");
1996 lega->Draw("same");
1997 cV0Efficiencye->cd(1);
1998 gPad->SetLogy();
1999 TH1D *afteree = (TH1D *) aftere->Clone();
2000 afteree->SetMarkerStyle(stylee[binc]);
2001 afteree->SetMarkerColor(colorr[binc]);
2002 if(binc==0) afteree->Draw();
2003 else afteree->Draw("same");
2004 legtotal->AddEntry(afteree,(const char*) titlee,"p");
2005 cV0Efficiencye->cd(2);
2006 gPad->SetLogy();
2007 TH1D *aftereeu = (TH1D *) aftere->Clone();
2008 aftereeu->SetMarkerStyle(stylee[binc]);
2009 aftereeu->SetMarkerColor(colorr[binc]);
2010 if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]);
2011 if(binc==0) aftereeu->Draw();
2012 else aftereeu->Draw("same");
2013 legtotalg->AddEntry(aftereeu,(const char*) titlee,"p");
2014 ccV0Efficiency->cd(2);
2015 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1);
2016 efficiencyDDproj->SetTitle("");
2017 efficiencyDDproj->SetStats(0);
2018 efficiencyDDproj->SetMarkerStyle(25);
2019 efficiencyDDproj->Draw();
c2690925 2020 }
2021
2022 cV0Efficiencye->cd(1);
2023 legtotal->Draw("same");
2024 cV0Efficiencye->cd(2);
2025 legtotalg->Draw("same");
2026
2027 cenaxisa->SetRange(0,nbbin);
2028 cenaxisb->SetRange(0,nbbin);
2029 cenaxisc->SetRange(0,nbbin);
e17c1f86 2030
2031 if(fWriteToFile) cV0Efficiencye->SaveAs("V0efficiency.eps");
c2690925 2032 }
3a72645a 2033
2034 }
2035
2036
2037 return result;
2038
2039}
c04c80e6 2040//____________________________________________________________
3a72645a 2041TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
c04c80e6 2042
2043 //
2044 // Unfold and eventually correct for efficiency the bgsubspectrum
2045 //
2046
3a72645a 2047 AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
c04c80e6 2048 if(!mcContainer){
2049 AliError("MC Container not available");
2050 return NULL;
2051 }
2052
2053 if(!fCorrelation){
2054 AliError("No Correlation map available");
2055 return NULL;
2056 }
2057
3a72645a 2058 // Data
c04c80e6 2059 AliCFDataGrid *dataGrid = 0x0;
2060 if(bgsubpectrum) {
c04c80e6 2061 dataGrid = bgsubpectrum;
c04c80e6 2062 }
2063 else {
2064
2065 AliCFContainer *dataContainer = GetContainer(kDataContainer);
2066 if(!dataContainer){
2067 AliError("Data Container not available");
2068 return NULL;
2069 }
2070
3a72645a 2071 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
c04c80e6 2072 }
2073
c04c80e6 2074 // Guessed
3a72645a 2075 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
c04c80e6 2076 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
2077
2078 // Efficiency
3a72645a 2079 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
c04c80e6 2080 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
8c1c76e9 2081
a8ef1999 2082 if(!fBeauty2ndMethod)
2083 {
2084 // Consider parameterized IP cut efficiency
2085 if(!fInclusiveSpectrum){
2086 Int_t* bins=new Int_t[1];
2087 bins[0]=35;
2088
2089 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
0e30407a 2090 beffContainer->SetGrid(GetBeautyIPEff(kTRUE));
a8ef1999 2091 efficiencyD->Multiply(beffContainer,1);
2092 }
8c1c76e9 2093 }
2094
c04c80e6 2095
2096 // Unfold
2097
3a72645a 2098 AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
c2690925 2099 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
c04c80e6 2100 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
e156c3bb 2101 if(fSetSmoothing) unfolding.UseSmoothing();
c04c80e6 2102 unfolding.Unfold();
2103
2104 // Results
2105 THnSparse* result = unfolding.GetUnfolded();
2106 THnSparse* residual = unfolding.GetEstMeasured();
2107 TList *listofresults = new TList;
2108 listofresults->SetOwner();
2109 listofresults->AddAt((THnSparse*)result->Clone(),0);
2110 listofresults->AddAt((THnSparse*)residual->Clone(),1);
2111
3a72645a 2112 if(fDebugLevel > 0) {
c2690925 2113
2114 Int_t ptpr;
2115 if(fBeamType==0) ptpr=0;
2116 if(fBeamType==1) ptpr=1;
3a72645a 2117
2118 TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
2119 cresidual->Divide(2,1);
2120 cresidual->cd(1);
2121 gPad->SetLogy();
2122 TGraphErrors* residualspectrumD = Normalize(residual);
2123 if(!residualspectrumD) {
2124 AliError("Number of Events not set for the normalization");
3ccf8f4c 2125 return NULL;
3a72645a 2126 }
2127 residualspectrumD->SetTitle("");
2128 residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
2129 residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
2130 residualspectrumD->SetMarkerStyle(26);
2131 residualspectrumD->SetMarkerColor(kBlue);
2132 residualspectrumD->SetLineColor(kBlue);
2133 residualspectrumD->Draw("AP");
2134 AliCFDataGrid *dataGridBis = (AliCFDataGrid *) (dataGrid->Clone());
2135 dataGridBis->SetName("dataGridBis");
2136 TGraphErrors* measuredspectrumD = Normalize(dataGridBis);
2137 measuredspectrumD->SetTitle("");
2138 measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
2139 measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
2140 measuredspectrumD->SetMarkerStyle(25);
2141 measuredspectrumD->SetMarkerColor(kBlack);
2142 measuredspectrumD->SetLineColor(kBlack);
2143 measuredspectrumD->Draw("P");
2144 TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
2145 legres->AddEntry(residualspectrumD,"Residual","p");
2146 legres->AddEntry(measuredspectrumD,"Measured","p");
2147 legres->Draw("same");
2148 cresidual->cd(2);
c2690925 2149 TH1D *residualTH1D = residual->Projection(ptpr);
2150 TH1D *measuredTH1D = (TH1D *) dataGridBis->Project(ptpr);
3a72645a 2151 TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
2152 ratioresidual->SetName("ratioresidual");
2153 ratioresidual->SetTitle("");
2154 ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
2155 ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2156 ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
2157 ratioresidual->SetStats(0);
2158 ratioresidual->Draw();
e17c1f86 2159
2160 if(fWriteToFile) cresidual->SaveAs("Unfolding.eps");
3a72645a 2161 }
2162
c04c80e6 2163 return listofresults;
2164
2165}
2166
2167//____________________________________________________________
3a72645a 2168AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
c04c80e6 2169
2170 //
2171 // Apply unfolding and efficiency correction together to bgsubspectrum
2172 //
2173
3a72645a 2174 AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
c04c80e6 2175 if(!mcContainer){
2176 AliError("MC Container not available");
2177 return NULL;
2178 }
2179
c04c80e6 2180 // Efficiency
3a72645a 2181 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
c04c80e6 2182 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
2183
a8ef1999 2184
2185 if(!fBeauty2ndMethod)
2186 {
2187 // Consider parameterized IP cut efficiency
2188 if(!fInclusiveSpectrum){
2189 Int_t* bins=new Int_t[1];
2190 bins[0]=35;
2191
2192 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
0e30407a 2193 beffContainer->SetGrid(GetBeautyIPEff(kFALSE));
a8ef1999 2194 efficiencyD->Multiply(beffContainer,1);
2195 }
8c1c76e9 2196 }
2197
c04c80e6 2198 // Data in the right format
2199 AliCFDataGrid *dataGrid = 0x0;
2200 if(bgsubpectrum) {
c04c80e6 2201 dataGrid = bgsubpectrum;
c04c80e6 2202 }
2203 else {
3a72645a 2204
c04c80e6 2205 AliCFContainer *dataContainer = GetContainer(kDataContainer);
2206 if(!dataContainer){
2207 AliError("Data Container not available");
2208 return NULL;
2209 }
2210
3a72645a 2211 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
c04c80e6 2212 }
2213
2214 // Correct
2215 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
2216 result->ApplyEffCorrection(*efficiencyD);
c04c80e6 2217
3a72645a 2218 if(fDebugLevel > 0) {
c2690925 2219
2220 Int_t ptpr;
2221 if(fBeamType==0) ptpr=0;
2222 if(fBeamType==1) ptpr=1;
3a72645a 2223
bf892a6a 2224 printf("Step MC: %d\n",fStepMC);
2225 printf("Step tracking: %d\n",AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
2226 printf("Step MC true: %d\n",fStepTrue);
3a72645a 2227 AliCFEffGrid *efficiencymcPID = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
2228 AliCFEffGrid *efficiencymctrackinggeo = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack,fStepTrue);
2229 AliCFEffGrid *efficiencymcall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,fStepTrue);
2230
2231 AliCFEffGrid *efficiencyesdall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerESD),fStepMC,fStepTrue);
2232
2233 TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
2234 cefficiency->cd(1);
c2690925 2235 TH1D* efficiencymcPIDD = (TH1D *) efficiencymcPID->Project(ptpr);
3a72645a 2236 efficiencymcPIDD->SetTitle("");
2237 efficiencymcPIDD->SetStats(0);
2238 efficiencymcPIDD->SetMarkerStyle(25);
2239 efficiencymcPIDD->Draw();
c2690925 2240 TH1D* efficiencymctrackinggeoD = (TH1D *) efficiencymctrackinggeo->Project(ptpr);
3a72645a 2241 efficiencymctrackinggeoD->SetTitle("");
2242 efficiencymctrackinggeoD->SetStats(0);
2243 efficiencymctrackinggeoD->SetMarkerStyle(26);
2244 efficiencymctrackinggeoD->Draw("same");
c2690925 2245 TH1D* efficiencymcallD = (TH1D *) efficiencymcall->Project(ptpr);
3a72645a 2246 efficiencymcallD->SetTitle("");
2247 efficiencymcallD->SetStats(0);
2248 efficiencymcallD->SetMarkerStyle(27);
2249 efficiencymcallD->Draw("same");
c2690925 2250 TH1D* efficiencyesdallD = (TH1D *) efficiencyesdall->Project(ptpr);
3a72645a 2251 efficiencyesdallD->SetTitle("");
2252 efficiencyesdallD->SetStats(0);
2253 efficiencyesdallD->SetMarkerStyle(24);
2254 efficiencyesdallD->Draw("same");
2255 TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
2256 legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
2257 legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
2258 legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
2259 legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
2260 legeff->Draw("same");
c2690925 2261
2262 if(fBeamType==1) {
2263
2264 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
2265 TAxis *cenaxisa = sparseafter->GetAxis(0);
2266 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
2267 TAxis *cenaxisb = sparsebefore->GetAxis(0);
2268 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
2269 TAxis *cenaxisc = efficiencya->GetAxis(0);
2270 //Int_t nbbin = cenaxisb->GetNbins();
2271 //Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
2272 //Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
2273 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
e17c1f86 2274 TString titlee("Efficiency_centrality_bin_");
2275 titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
2276 titlee += "_";
2277 titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
2278 TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
2279 cefficiencye->Divide(2,1);
2280 cefficiencye->cd(1);
2281 gPad->SetLogy();
2282 cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
2283 cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
2284 TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr);
2285 TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr);
2286 CorrectFromTheWidth(afterefficiency);
2287 CorrectFromTheWidth(beforeefficiency);
2288 afterefficiency->SetStats(0);
2289 afterefficiency->SetTitle((const char*)titlee);
2290 afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2291 afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2292 afterefficiency->SetMarkerStyle(25);
2293 afterefficiency->SetMarkerColor(kBlack);
2294 afterefficiency->SetLineColor(kBlack);
2295 beforeefficiency->SetStats(0);
2296 beforeefficiency->SetTitle((const char*)titlee);
2297 beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2298 beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2299 beforeefficiency->SetMarkerStyle(24);
2300 beforeefficiency->SetMarkerColor(kBlue);
2301 beforeefficiency->SetLineColor(kBlue);
2302 afterefficiency->Draw();
2303 beforeefficiency->Draw("same");
2304 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
2305 lega->AddEntry(beforeefficiency,"Before efficiency correction","p");
2306 lega->AddEntry(afterefficiency,"After efficiency correction","p");
2307 lega->Draw("same");
2308 cefficiencye->cd(2);
2309 cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1);
2310 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr);
2311 efficiencyDDproj->SetTitle("");
2312 efficiencyDDproj->SetStats(0);
2313 efficiencyDDproj->SetMarkerStyle(25);
2314 efficiencyDDproj->SetMarkerColor(2);
2315 efficiencyDDproj->Draw();
2316 cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]);
2317 TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr);
2318 efficiencyDDproja->SetTitle("");
2319 efficiencyDDproja->SetStats(0);
2320 efficiencyDDproja->SetMarkerStyle(26);
2321 efficiencyDDproja->SetMarkerColor(4);
2322 efficiencyDDproja->Draw("same");
c2690925 2323 }
2324 }
2325
e17c1f86 2326 if(fWriteToFile) cefficiency->SaveAs("efficiencyCorrected.eps");
3a72645a 2327 }
2328
c04c80e6 2329 return result;
2330
2331}
3a72645a 2332
c04c80e6 2333//__________________________________________________________________________________
c2690925 2334TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum,Int_t i) const {
c04c80e6 2335 //
2336 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2337 // Give the final pt spectrum to be compared
2338 //
cedf0381 2339
c2690925 2340 if(fNEvents[i] > 0) {
2341
a199006c 2342 Int_t ptpr = 0;
c2690925 2343 if(fBeamType==0) ptpr=0;
2344 if(fBeamType==1) ptpr=1;
c04c80e6 2345
c2690925 2346 TH1D* projection = spectrum->Projection(ptpr);
c04c80e6 2347 CorrectFromTheWidth(projection);
c2690925 2348 TGraphErrors *graphError = NormalizeTH1(projection,i);
c04c80e6 2349 return graphError;
2350
2351 }
2352
2353 return 0x0;
2354
2355
2356}
2357//__________________________________________________________________________________
c2690925 2358TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i) const {
c04c80e6 2359 //
2360 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2361 // Give the final pt spectrum to be compared
2362 //
c2690925 2363 if(fNEvents[i] > 0) {
c04c80e6 2364
a199006c 2365 Int_t ptpr=0;
c2690925 2366 if(fBeamType==0) ptpr=0;
2367 if(fBeamType==1) ptpr=1;
2368
2369 TH1D* projection = (TH1D *) spectrum->Project(ptpr);
c04c80e6 2370 CorrectFromTheWidth(projection);
c2690925 2371 TGraphErrors *graphError = NormalizeTH1(projection,i);
a8ef1999 2372
c04c80e6 2373 return graphError;
2374
2375 }
2376
2377 return 0x0;
2378
2379
2380}
2381//__________________________________________________________________________________
c2690925 2382TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
2383 //
2384 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2385 // Give the final pt spectrum to be compared
2386 //
8c1c76e9 2387 Double_t chargecoefficient = 0.5;
e17c1f86 2388 if(fChargeChoosen != 0) chargecoefficient = 1.0;
8c1c76e9 2389
e17c1f86 2390 Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
8c1c76e9 2391 printf("Normalizing Eta Range %f\n", etarange);
c2690925 2392 if(fNEvents[i] > 0) {
2393
2394 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
2395 Double_t p = 0, dp = 0; Int_t point = 1;
2396 Double_t n = 0, dN = 0;
2397 Double_t nCorr = 0, dNcorr = 0;
7483e78e 2398 //Double_t errdN = 0, errdp = 0;
2399 Double_t errdN = 0;
c2690925 2400 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
2401 point = ibin - input->GetXaxis()->GetFirst();
2402 p = input->GetXaxis()->GetBinCenter(ibin);
7483e78e 2403 //dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
c2690925 2404 n = input->GetBinContent(ibin);
2405 dN = input->GetBinError(ibin);
c2690925 2406 // New point
8c1c76e9 2407 nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
c2690925 2408 errdN = 1./(2. * TMath::Pi() * p);
7483e78e 2409 //errdp = 1./(2. * TMath::Pi() * p*p) * n;
2410 //dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
2411 dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN);
c2690925 2412
2413 spectrumNormalized->SetPoint(point, p, nCorr);
2414 spectrumNormalized->SetPointError(point, dp, dNcorr);
2415 }
2416 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2417 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2418 spectrumNormalized->SetMarkerStyle(22);
2419 spectrumNormalized->SetMarkerColor(kBlue);
2420 spectrumNormalized->SetLineColor(kBlue);
a8ef1999 2421
c2690925 2422 return spectrumNormalized;
2423
2424 }
2425
2426 return 0x0;
2427
2428
2429}
2430//__________________________________________________________________________________
2431TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) const {
c04c80e6 2432 //
2433 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2434 // Give the final pt spectrum to be compared
2435 //
8c1c76e9 2436 Double_t chargecoefficient = 0.5;
e17c1f86 2437 if(fChargeChoosen != kAllCharge) chargecoefficient = 1.0;
8c1c76e9 2438
e17c1f86 2439 Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
8c1c76e9 2440 printf("Normalizing Eta Range %f\n", etarange);
c2690925 2441 if(normalization > 0) {
c04c80e6 2442
2443 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
2444 Double_t p = 0, dp = 0; Int_t point = 1;
2445 Double_t n = 0, dN = 0;
2446 Double_t nCorr = 0, dNcorr = 0;
7483e78e 2447 //Double_t errdN = 0, errdp = 0;
2448 Double_t errdN = 0;
c04c80e6 2449 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
2450 point = ibin - input->GetXaxis()->GetFirst();
2451 p = input->GetXaxis()->GetBinCenter(ibin);
7483e78e 2452 //dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
c04c80e6 2453 n = input->GetBinContent(ibin);
2454 dN = input->GetBinError(ibin);
c04c80e6 2455 // New point
8c1c76e9 2456 nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
c04c80e6 2457 errdN = 1./(2. * TMath::Pi() * p);
7483e78e 2458 //errdp = 1./(2. * TMath::Pi() * p*p) * n;
2459 //dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
2460 dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN);
c04c80e6 2461
2462 spectrumNormalized->SetPoint(point, p, nCorr);
2463 spectrumNormalized->SetPointError(point, dp, dNcorr);
2464 }
2465 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2466 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2467 spectrumNormalized->SetMarkerStyle(22);
2468 spectrumNormalized->SetMarkerColor(kBlue);
2469 spectrumNormalized->SetLineColor(kBlue);
2470
2471 return spectrumNormalized;
2472
2473 }
2474
2475 return 0x0;
2476
2477
2478}
2479//____________________________________________________________
2480void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
2481 //
2482 // Set the container for a given step to the
2483 //
e17c1f86 2484 if(!fCFContainers) fCFContainers = new TObjArray(kDataContainerV0+1);
c04c80e6 2485 fCFContainers->AddAt(cont, type);
2486}
2487
2488//____________________________________________________________
2489AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
2490 //
2491 // Get Correction Framework Container for given type
2492 //
2493 if(!fCFContainers) return NULL;
2494 return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
2495}
c04c80e6 2496//____________________________________________________________
0e30407a 2497AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centralitylow, Int_t centralityhigh) {
c04c80e6 2498 //
3a72645a 2499 // Slice bin for a given source of electron
c2690925 2500 // nDim is the number of dimension the corrections are done
2501 // dimensions are the definition of the dimensions
2502 // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
2503 // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
a8ef1999 2504 // centrality (-1 means we do not cut on centrality)
c04c80e6 2505 //
2506
2507 Double_t *varMin = new Double_t[container->GetNVar()],
2508 *varMax = new Double_t[container->GetNVar()];
2509
2510 Double_t *binLimits;
2511 for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
2512
2513 binLimits = new Double_t[container->GetNBins(ivar)+1];
2514 container->GetBinLimits(ivar,binLimits);
c2690925 2515 varMin[ivar] = binLimits[0];
2516 varMax[ivar] = binLimits[container->GetNBins(ivar)];
2517 // source
2518 if(ivar == 4){
2519 if((source>= 0) && (source<container->GetNBins(ivar))) {
ff8249bd 2520 varMin[ivar] = container->GetAxis(4,0)->GetBinLowEdge(container->GetAxis(4,0)->FindBin(binLimits[source]));
2521 varMax[ivar] = container->GetAxis(4,0)->GetBinUpEdge(container->GetAxis(4,0)->FindBin(binLimits[source]));
c2690925 2522 }
3a72645a 2523 }
c2690925 2524 // charge
2525 if(ivar == 3) {
ff8249bd 2526 if(charge != kAllCharge){
2527 varMin[ivar] = container->GetAxis(3,0)->GetBinLowEdge(container->GetAxis(3,0)->FindBin(charge));
2528 varMax[ivar] = container->GetAxis(3,0)->GetBinUpEdge(container->GetAxis(3,0)->FindBin(charge));
2529 }
3a72645a 2530 }
8c1c76e9 2531 // eta
2532 if(ivar == 1){
e17c1f86 2533 for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++)
2534 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 2535 if(fEtaSelected){
2536 varMin[ivar] = fEtaRange[0];
2537 varMax[ivar] = fEtaRange[1];
2538 }
2539 }
e17c1f86 2540 if(fEtaSelected){
2541 fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0]));
2542 fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1]));
2543 AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[0]));
2544 }
a8ef1999 2545 if(ivar == 5){
0e30407a 2546 if((centralitylow>= 0) && (centralitylow<container->GetNBins(ivar)) && (centralityhigh>= 0) && (centralityhigh<container->GetNBins(ivar))) {
2547 varMin[ivar] = binLimits[centralitylow];
2548 varMax[ivar] = binLimits[centralityhigh];
2549
2550 TAxis *axistest = container->GetAxis(5,0);
2551 printf("GetSlicedContainer: Number of bin in centrality direction %d\n",axistest->GetNbins());
2552 printf("Project from %f to %f\n",binLimits[centralitylow],binLimits[centralityhigh]);
2553 Double_t lowcentrality = axistest->GetBinLowEdge(axistest->FindBin(binLimits[centralitylow]));
2554 Double_t highcentrality = axistest->GetBinUpEdge(axistest->FindBin(binLimits[centralityhigh]));
2555 printf("GetSlicedContainer: Low centrality %f and high centrality %f\n",lowcentrality,highcentrality);
2556
a8ef1999 2557 }
0e30407a 2558
a8ef1999 2559 }
3a72645a 2560
11ff28c5 2561
c04c80e6 2562 delete[] binLimits;
2563 }
2564
2565 AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
2566 AddTemporaryObject(k);
2567 delete[] varMin; delete[] varMax;
2568
2569 return k;
2570
2571}
2572
2573//_________________________________________________________________________
0e30407a 2574THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions, Int_t centralitylow, Int_t centralityhigh) const {
c04c80e6 2575 //
3a72645a 2576 // Slice correlation
c04c80e6 2577 //
2578
3a72645a 2579 Int_t ndimensions = correlationmatrix->GetNdimensions();
c2690925 2580 //printf("Number of dimension %d correlation map\n",ndimensions);
c04c80e6 2581 if(ndimensions < (2*nDim)) {
2582 AliError("Problem in the dimensions");
2583 return NULL;
2584 }
0e30407a 2585
2586 // Cut in centrality is centrality > -1
2587 if((centralitylow >=0) && (centralityhigh >=0)) {
2588
2589 TAxis *axiscentrality0 = correlationmatrix->GetAxis(5);
2590 TAxis *axiscentrality1 = correlationmatrix->GetAxis(5+((Int_t)(ndimensions/2.)));
2591
2592 Int_t bins0 = axiscentrality0->GetNbins();
2593 Int_t bins1 = axiscentrality1->GetNbins();
2594
2595 printf("Number of centrality bins: %d and %d\n",bins0,bins1);
2596 if(bins0 != bins1) {
2597 AliError("Problem in the dimensions");
2598 return NULL;
2599 }
2600
2601 if((centralitylow>= 0) && (centralitylow<bins0) && (centralityhigh>= 0) && (centralityhigh<bins0)) {
2602 axiscentrality0->SetRangeUser(centralitylow,centralityhigh);
2603 axiscentrality1->SetRangeUser(centralitylow,centralityhigh);
2604
2605 Double_t lowcentrality0 = axiscentrality0->GetBinLowEdge(axiscentrality0->FindBin(centralitylow));
2606 Double_t highcentrality0 = axiscentrality0->GetBinUpEdge(axiscentrality0->FindBin(centralityhigh));
2607 Double_t lowcentrality1 = axiscentrality1->GetBinLowEdge(axiscentrality1->FindBin(centralitylow));
2608 Double_t highcentrality1 = axiscentrality1->GetBinUpEdge(axiscentrality1->FindBin(centralityhigh));
2609 printf("GetCorrelation0: Low centrality %f and high centrality %f\n",lowcentrality0,highcentrality0);
2610 printf("GetCorrelation1: Low centrality %f and high centrality %f\n",lowcentrality1,highcentrality1);
2611
2612 TCanvas * ctestcorrelation = new TCanvas("testcorrelationprojection","testcorrelationprojection",1000,700);
2613 ctestcorrelation->cd(1);
2614 TH2D* projection = (TH2D *) correlationmatrix->Projection(5,5+((Int_t)(ndimensions/2.)));
2615 projection->Draw("colz");
2616
2617 }
2618
2619 }
2620
2621
c04c80e6 2622 Int_t ndimensionsContainer = (Int_t) ndimensions/2;
c2690925 2623 //printf("Number of dimension %d container\n",ndimensionsContainer);
c04c80e6 2624
2625 Int_t *dim = new Int_t[nDim*2];
2626 for(Int_t iter=0; iter < nDim; iter++){
2627 dim[iter] = dimensions[iter];
2628 dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
c2690925 2629 //printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
c04c80e6 2630 }
2631
3a72645a 2632 THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
c04c80e6 2633
2634 delete[] dim;
2635 return k;
2636
2637}
2638//___________________________________________________________________________
2639void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
2640 //
2641 // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
2642 //
2643
2644 TAxis *axis = h1->GetXaxis();
2645 Int_t nbinX = h1->GetNbinsX();
2646
2647 for(Int_t i = 1; i <= nbinX; i++) {
2648
2649 Double_t width = axis->GetBinWidth(i);
2650 Double_t content = h1->GetBinContent(i);
2651 Double_t error = h1->GetBinError(i);
2652 h1->SetBinContent(i,content/width);
2653 h1->SetBinError(i,error/width);
2654 }
2655
2656}
8c1c76e9 2657
2658//___________________________________________________________________________
2659void AliHFEspectrum::CorrectStatErr(AliCFDataGrid *backgroundGrid) const {
2660 //
2661 // Correct statistical error
2662 //
2663
2664 TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
2665 Int_t nbinX = h1->GetNbinsX();
2666 Int_t bins[1];
2667 for(Long_t i = 1; i <= nbinX; i++) {
2668 bins[0] = i;
2669 Float_t content = h1->GetBinContent(i);
2670 if(content>0){
2671 Float_t error = TMath::Sqrt(content);
2672 backgroundGrid->SetElementError(bins, error);
2673 }
2674 }
2675}
2676
c04c80e6 2677//____________________________________________________________
2678void AliHFEspectrum::AddTemporaryObject(TObject *o){
2679 //
2680 // Emulate garbage collection on class level
2681 //
2682 if(!fTemporaryObjects) {
2683 fTemporaryObjects= new TList;
2684 fTemporaryObjects->SetOwner();
2685 }
2686 fTemporaryObjects->Add(o);
2687}
2688
2689//____________________________________________________________
2690void AliHFEspectrum::ClearObject(TObject *o){
2691 //
2692 // Do a safe deletion
2693 //
2694 if(fTemporaryObjects){
2695 if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
2696 delete o;
2697 } else{
2698 // just delete
2699 delete o;
2700 }
2701}
2702//_________________________________________________________________________
c2690925 2703TObject* AliHFEspectrum::GetSpectrum(const AliCFContainer * const c, Int_t step) {
245877d0 2704 AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
c04c80e6 2705 return data;
2706}
2707//_________________________________________________________________________
c2690925 2708TObject* AliHFEspectrum::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
c04c80e6 2709 //
2710 // Create efficiency grid and calculate efficiency
2711 // of step to step0
2712 //
2713 TString name("eff");
2714 name += step;
2715 name+= step0;
2716 AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
2717 eff->CalculateEfficiency(step,step0);
2718 return eff;
2719}
c2690925 2720
2721//____________________________________________________________________________
8c1c76e9 2722THnSparse* AliHFEspectrum::GetCharmWeights(){
c2690925 2723
2724 //
2725 // Measured D->e based weighting factors
2726 //
2727
a8ef1999 2728 const Int_t nDimpp=1;
2729 Int_t nBinpp[nDimpp] = {35};
e17c1f86 2730 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 2731 const Int_t nDimPbPb=2;
2732 Int_t nBinPbPb[nDimPbPb] = {11,35};
2733 Double_t kCentralityRange[12] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
2734 Int_t loopcentr=1;
2735 Int_t looppt=nBinpp[0];
2736 if(fBeamType==0)
2737 {
2738 fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
2739 fWeightCharm->SetBinEdges(0, ptbinning1);
2740 }
2741 if(fBeamType==1)
2742 {
2743 fWeightCharm = new THnSparseF("weightHisto", "weighting factor; centrality bin; pt[GeV/c]", nDimPbPb, nBinPbPb);
2744 fWeightCharm->SetBinEdges(1, ptbinning1);
2745 fWeightCharm->SetBinEdges(0, kCentralityRange);
2746 loopcentr=nBinPbPb[0];
2747 }
2748
0e30407a 2749 // Weighting factor for pp
11ff28c5 2750 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};
2751
0e30407a 2752 // Weighting factor for PbPb (0-20%)
2753 //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};
2754
2755 // Weighting factor for PbPb (40-80%)
2756 //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 2757
c2690925 2758 //points
8c1c76e9 2759 Double_t pt[1];
a8ef1999 2760 Double_t contents[2];
2761
2762 for(int icentr=0; icentr<loopcentr; icentr++)
2763 {
2764 for(int i=0; i<looppt; i++){
2765 pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
2766 if(fBeamType==1)
2767 {
2768 contents[0]=icentr;
2769 contents[1]=pt[0];
2770 }
2771 if(fBeamType==0)
2772 {
2773 contents[0]=pt[0];
2774 }
2775
2776 fWeightCharm->Fill(contents,weight[i]);
2777 }
2778 }
2779
2780 Int_t nDimSparse = fWeightCharm->GetNdimensions();
2781 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
2782 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
2783
2784 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2785 binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
2786 nBins *= binsvar[iVar];
2787 }
2788
2789 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
2790 // loop that sets 0 error in each bin
2791 for (Long_t iBin=0; iBin<nBins; iBin++) {
2792 Long_t bintmp = iBin ;
2793 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2794 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
2795 bintmp /= binsvar[iVar] ;
2796 }
2797 fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
c2690925 2798 }
c2690925 2799
dcef324e 2800 delete[] binsvar;
2801 delete[] binfill;
2802
c2690925 2803 return fWeightCharm;
2804}
8c1c76e9 2805
2806//____________________________________________________________________________
11ff28c5 2807void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
a8ef1999 2808
2809 // TOF PID efficiencies
2810 Int_t ptpr;
2811 if(fBeamType==0) ptpr=0;
2812 if(fBeamType==1) ptpr=1;
2813
2814 Int_t loopcentr=1;
2815 const Int_t nCentralitybinning=11; //number of centrality bins
2816 if(fBeamType==1)
2817 {
2818 loopcentr=nCentralitybinning;
2819 }
8c1c76e9 2820
cedf0381 2821 TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
2822 TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
0e30407a 2823 TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.3,10.0);
a8ef1999 2824
0e30407a 2825 TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600);
2826 cefficiencyParamtof->cd();
a8ef1999 2827
2828 AliCFContainer *mccontainermcD = 0x0;
0e30407a 2829 AliCFContainer *mccontaineresdD = 0x0;
a8ef1999 2830 TH1D* efficiencysigTOFPIDD[nCentralitybinning];
2831 TH1D* efficiencyTOFPIDD[nCentralitybinning];
11ff28c5 2832 TH1D* efficiencysigesdTOFPIDD[nCentralitybinning];
2833 TH1D* efficiencyesdTOFPIDD[nCentralitybinning];
a8ef1999 2834 Int_t source = -1; //get parameterized TOF PID efficiencies
2835
2836 for(int icentr=0; icentr<loopcentr; icentr++) {
2837 // signal sample
2838 if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 2839 else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
a8ef1999 2840 AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
2841 efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2842
2843 // mb sample for double check
2844 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 2845 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
a8ef1999 2846 AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
2847 efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2848
11ff28c5 2849 // mb sample with reconstructed variables
2850 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
2851 else mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2852 AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD);
2853 efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2854
2855 // mb sample with reconstructed variables
2856 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
2857 else mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2858 AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
2859 efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2860
2861 //fill histo
a8ef1999 2862 efficiencysigTOFPIDD[icentr] = (TH1D *) efficiencymcsigParamTOFPID->Project(ptpr);
2863 efficiencyTOFPIDD[icentr] = (TH1D *) efficiencymcParamTOFPID->Project(ptpr);
11ff28c5 2864 efficiencysigesdTOFPIDD[icentr] = (TH1D *) efficiencysigesdParamTOFPID->Project(ptpr);
2865 efficiencyesdTOFPIDD[icentr] = (TH1D *) efficiencyesdParamTOFPID->Project(ptpr);
a8ef1999 2866 efficiencysigTOFPIDD[icentr]->SetName(Form("efficiencysigTOFPIDD%d",icentr));
2867 efficiencyTOFPIDD[icentr]->SetName(Form("efficiencyTOFPIDD%d",icentr));
11ff28c5 2868 efficiencysigesdTOFPIDD[icentr]->SetName(Form("efficiencysigesdTOFPIDD%d",icentr));
2869 efficiencyesdTOFPIDD[icentr]->SetName(Form("efficiencyesdTOFPIDD%d",icentr));
a8ef1999 2870
11ff28c5 2871 //fit (mc enhenced sample)
a8ef1999 2872 fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
a8ef1999 2873 efficiencysigTOFPIDD[icentr]->Fit(fittofpid,"R");
2874 efficiencysigTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
2875 fEfficiencyTOFPIDD[icentr] = efficiencysigTOFPIDD[icentr]->GetFunction("fittofpid");
11ff28c5 2876
2877 //fit (esd enhenced sample)
2878 efficiencysigesdTOFPIDD[icentr]->Fit(fittofpid,"R");
2879 efficiencysigesdTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
2880 fEfficiencyesdTOFPIDD[icentr] = efficiencysigesdTOFPIDD[icentr]->GetFunction("fittofpid");
2881
a8ef1999 2882 }
8c1c76e9 2883
a8ef1999 2884 // draw (for PbPb, only 1st bin)
11ff28c5 2885 //sig mc
a8ef1999 2886 efficiencysigTOFPIDD[0]->SetTitle("");
2887 efficiencysigTOFPIDD[0]->SetStats(0);
2888 efficiencysigTOFPIDD[0]->SetMarkerStyle(25);
2889 efficiencysigTOFPIDD[0]->SetMarkerColor(2);
2890 efficiencysigTOFPIDD[0]->SetLineColor(2);
2891 efficiencysigTOFPIDD[0]->Draw();
8c1c76e9 2892
11ff28c5 2893 //mb mc
a8ef1999 2894 efficiencyTOFPIDD[0]->SetTitle("");
2895 efficiencyTOFPIDD[0]->SetStats(0);
2896 efficiencyTOFPIDD[0]->SetMarkerStyle(24);
11ff28c5 2897 efficiencyTOFPIDD[0]->SetMarkerColor(4);
2898 efficiencyTOFPIDD[0]->SetLineColor(4);
a8ef1999 2899 efficiencyTOFPIDD[0]->Draw("same");
2900
11ff28c5 2901 //sig esd
2902 efficiencysigesdTOFPIDD[0]->SetTitle("");
2903 efficiencysigesdTOFPIDD[0]->SetStats(0);
2904 efficiencysigesdTOFPIDD[0]->SetMarkerStyle(25);
2905 efficiencysigesdTOFPIDD[0]->SetMarkerColor(3);
2906 efficiencysigesdTOFPIDD[0]->SetLineColor(3);
2907 efficiencysigesdTOFPIDD[0]->Draw("same");
2908
2909 //mb esd
2910 efficiencyesdTOFPIDD[0]->SetTitle("");
2911 efficiencyesdTOFPIDD[0]->SetStats(0);
2912 efficiencyesdTOFPIDD[0]->SetMarkerStyle(25);
2913 efficiencyesdTOFPIDD[0]->SetMarkerColor(1);
2914 efficiencyesdTOFPIDD[0]->SetLineColor(1);
2915 efficiencyesdTOFPIDD[0]->Draw("same");
2916
2917 //signal mc fit
cedf0381 2918 if(fEfficiencyTOFPIDD[0]){
2919 fEfficiencyTOFPIDD[0]->SetLineColor(2);
2920 fEfficiencyTOFPIDD[0]->Draw("same");
2921 }
11ff28c5 2922 //mb esd fit
cedf0381 2923 if(fEfficiencyesdTOFPIDD[0]){
2924 fEfficiencyesdTOFPIDD[0]->SetLineColor(3);
2925 fEfficiencyesdTOFPIDD[0]->Draw("same");
2926 }
a8ef1999 2927
11ff28c5 2928 TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
2929 legtofeff->AddEntry(efficiencysigTOFPIDD[0],"TOF PID Step Efficiency","");
2930 legtofeff->AddEntry(efficiencysigTOFPIDD[0],"vs MC p_{t} for enhenced samples","p");
2931 legtofeff->AddEntry(efficiencyTOFPIDD[0],"vs MC p_{t} for mb samples","p");
2932 legtofeff->AddEntry(efficiencysigesdTOFPIDD[0],"vs esd p_{t} for enhenced samples","p");
2933 legtofeff->AddEntry(efficiencyesdTOFPIDD[0],"vs esd p_{t} for mb samples","p");
2934 legtofeff->Draw("same");
a8ef1999 2935
11ff28c5 2936
0e30407a 2937 TCanvas * cefficiencyParamIP = new TCanvas("efficiencyParamIP","efficiencyParamIP",500,500);
2938 cefficiencyParamIP->cd();
2939 gStyle->SetOptStat(0);
2940
a8ef1999 2941 // IP cut efficiencies
a8ef1999 2942 for(int icentr=0; icentr<loopcentr; icentr++) {
11ff28c5 2943
2944 AliCFContainer *charmCombined = 0x0;
2945 AliCFContainer *beautyCombined = 0x0;
0e30407a 2946 AliCFContainer *beautyCombinedesd = 0x0;
11ff28c5 2947
2948 printf("centrality printing %i \n",icentr);
2949
2950 source = 0; //charm enhenced
a8ef1999 2951 if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 2952 else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
a8ef1999 2953 AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
2954 efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2955
11ff28c5 2956 charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");
2957
2958 source = 1; //beauty enhenced
a8ef1999 2959 if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 2960 else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
a8ef1999 2961 AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
2962 efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2963
11ff28c5 2964 beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined");
2965
0e30407a 2966 if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
2967 else mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2968 AliCFEffGrid* efficiencyBeautySigesd = new AliCFEffGrid("efficiencyBeautySigesd","",*mccontaineresdD);
2969 efficiencyBeautySigesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2970
2971 beautyCombinedesd = (AliCFContainer*)mccontaineresdD->Clone("beautyCombinedesd");
2972
11ff28c5 2973 source = 0; //charm 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* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
2977 efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2978
11ff28c5 2979 charmCombined->Add(mccontainermcD);
2980 AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined);
2981 efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
2982
2983 source = 1; //beauty mb
a8ef1999 2984 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 2985 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
a8ef1999 2986 AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
2987 efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2988
11ff28c5 2989 beautyCombined->Add(mccontainermcD);
2990 AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined);
2991 efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
2992
0e30407a 2993 if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
2994 else mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2995 AliCFEffGrid* efficiencyBeautyesd = new AliCFEffGrid("efficiencyBeautyesd","",*mccontaineresdD);
2996 efficiencyBeautyesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2997
2998 beautyCombinedesd->Add(mccontaineresdD);
2999 AliCFEffGrid* efficiencyBeautyCombinedesd = new AliCFEffGrid("efficiencyBeautyCombinedesd","",*beautyCombinedesd);
3000 efficiencyBeautyCombinedesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
3001
11ff28c5 3002 source = 2; //conversion mb
a8ef1999 3003 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 3004 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
a8ef1999 3005 AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
3006 efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
3007
11ff28c5 3008 source = 3; //non HFE except for the conversion mb
a8ef1999 3009 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
dcef324e 3010 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
a8ef1999 3011 AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
3012 efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
3013
11ff28c5 3014 if(fIPEffCombinedSamples){
3015 fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmCombined->Project(ptpr); //signal enhenced + mb
3016 fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautyCombined->Project(ptpr); //signal enhenced + mb
0e30407a 3017 fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautyCombinedesd->Project(ptpr); //signal enhenced + mb
11ff28c5 3018 }
3019 else{
3020 fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmSig->Project(ptpr); //signal enhenced only
3021 fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautySig->Project(ptpr); //signal enhenced only
0e30407a 3022 fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautySigesd->Project(ptpr); //signal enhenced only
11ff28c5 3023 }
3024 fCharmEff[icentr] = (TH1D*)efficiencyCharm->Project(ptpr); //mb only
3025 fBeautyEff[icentr] = (TH1D*)efficiencyBeauty->Project(ptpr); //mb only
3026 fConversionEff[icentr] = (TH1D*)efficiencyConv->Project(ptpr); //mb only
3027 fNonHFEEff[icentr] = (TH1D*)efficiencyNonhfe->Project(ptpr); //mb only
3028
a8ef1999 3029 }
3030
0e30407a 3031 if(fBeamType==0){
3032 AliCFEffGrid *nonHFEEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
3033 fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
3034
3035 AliCFEffGrid *conversionEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
3036 fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
3037 }
3038
a8ef1999 3039 for(int icentr=0; icentr<loopcentr; icentr++) {
3040 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3041 fipfit->SetLineColor(2);
3042 fEfficiencyBeautySigD[icentr]->Fit(fipfit,"R");
3043 fEfficiencyBeautySigD[icentr]->GetYaxis()->SetTitle("Efficiency");
11ff28c5 3044 if(fBeauty2ndMethod)fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[0]->GetFunction("fipfit"); //why do we need this line?
a8ef1999 3045 else fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[icentr]->GetFunction("fipfit");
3046
0e30407a 3047 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3048 fipfit->SetLineColor(6);
3049 fEfficiencyBeautySigesdD[icentr]->Fit(fipfit,"R");
3050 fEfficiencyBeautySigesdD[icentr]->GetYaxis()->SetTitle("Efficiency");
3051 if(fBeauty2ndMethod)fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[0]->GetFunction("fipfit"); //why do we need this line?
3052 else fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[icentr]->GetFunction("fipfit");
3053
11ff28c5 3054 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3055 fipfit->SetLineColor(1);
3056 fEfficiencyCharmSigD[icentr]->Fit(fipfit,"R");
3057 fEfficiencyCharmSigD[icentr]->GetYaxis()->SetTitle("Efficiency");
3058 fEfficiencyIPCharmD[icentr] = fEfficiencyCharmSigD[icentr]->GetFunction("fipfit");
cedf0381 3059
11ff28c5 3060 if(fIPParameterizedEff){
0e30407a 3061 fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3062 fipfitnonhfe->SetLineColor(3);
3063 fConversionEff[icentr]->Fit(fipfitnonhfe,"R");
a8ef1999 3064 fConversionEff[icentr]->GetYaxis()->SetTitle("Efficiency");
0e30407a 3065 fEfficiencyIPConversionD[icentr] = fConversionEff[icentr]->GetFunction("fipfitnonhfe");
a8ef1999 3066
0e30407a 3067 fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3068 fipfitnonhfe->SetLineColor(4);
3069 fNonHFEEff[icentr]->Fit(fipfitnonhfe,"R");
a8ef1999 3070 fNonHFEEff[icentr]->GetYaxis()->SetTitle("Efficiency");
0e30407a 3071 fEfficiencyIPNonhfeD[icentr] = fNonHFEEff[icentr]->GetFunction("fipfitnonhfe");
a8ef1999 3072 }
3073 }
3074
3075 // draw (for PbPb, only 1st bin)
a8ef1999 3076 fEfficiencyCharmSigD[0]->SetMarkerStyle(21);
3077 fEfficiencyCharmSigD[0]->SetMarkerColor(1);
3078 fEfficiencyCharmSigD[0]->SetLineColor(1);
3079 fEfficiencyBeautySigD[0]->SetMarkerStyle(21);
3080 fEfficiencyBeautySigD[0]->SetMarkerColor(2);
3081 fEfficiencyBeautySigD[0]->SetLineColor(2);
0e30407a 3082 fEfficiencyBeautySigesdD[0]->SetStats(0);
3083 fEfficiencyBeautySigesdD[0]->SetMarkerStyle(21);
3084 fEfficiencyBeautySigesdD[0]->SetMarkerColor(6);
3085 fEfficiencyBeautySigesdD[0]->SetLineColor(6);
a8ef1999 3086 fCharmEff[0]->SetMarkerStyle(24);
3087 fCharmEff[0]->SetMarkerColor(1);
3088 fCharmEff[0]->SetLineColor(1);
3089 fBeautyEff[0]->SetMarkerStyle(24);
3090 fBeautyEff[0]->SetMarkerColor(2);
3091 fBeautyEff[0]->SetLineColor(2);
3092 fConversionEff[0]->SetMarkerStyle(24);
3093 fConversionEff[0]->SetMarkerColor(3);
3094 fConversionEff[0]->SetLineColor(3);
3095 fNonHFEEff[0]->SetMarkerStyle(24);
3096 fNonHFEEff[0]->SetMarkerColor(4);
3097 fNonHFEEff[0]->SetLineColor(4);
3098
3099 fEfficiencyCharmSigD[0]->Draw();
0e30407a 3100 fEfficiencyCharmSigD[0]->GetXaxis()->SetRangeUser(0.0,7.9);
3101 fEfficiencyCharmSigD[0]->GetYaxis()->SetRangeUser(0.0,0.5);
3102
a8ef1999 3103 fEfficiencyBeautySigD[0]->Draw("same");
0e30407a 3104 fEfficiencyBeautySigesdD[0]->Draw("same");
11ff28c5 3105 //fCharmEff[0]->Draw("same");
3106 //fBeautyEff[0]->Draw("same");
0e30407a 3107
3108 if(fBeamType==0){
3109 fConversionEffbgc->SetMarkerStyle(25);
3110 fConversionEffbgc->SetMarkerColor(3);
3111 fConversionEffbgc->SetLineColor(3);
3112 fNonHFEEffbgc->SetMarkerStyle(25);
3113 fNonHFEEffbgc->SetMarkerColor(4);
3114 fNonHFEEffbgc->SetLineColor(4);
3115 fConversionEffbgc->Draw("same");
3116 fNonHFEEffbgc->Draw("same");
3117 }
3118 else{
3119 fConversionEff[0]->Draw("same");
3120 fNonHFEEff[0]->Draw("same");
3121 }
cedf0381 3122 if(fEfficiencyIPBeautyD[0])
3123 fEfficiencyIPBeautyD[0]->Draw("same");
3124 if(fEfficiencyIPBeautyesdD[0])
3125 fEfficiencyIPBeautyesdD[0]->Draw("same");
3126 if( fEfficiencyIPCharmD[0])
3127 fEfficiencyIPCharmD[0]->Draw("same");
a8ef1999 3128 if(fIPParameterizedEff){
a8ef1999 3129 fEfficiencyIPConversionD[0]->Draw("same");
3130 fEfficiencyIPNonhfeD[0]->Draw("same");
3131 }
0e30407a 3132 TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39);
11ff28c5 3133 legipeff->AddEntry(fEfficiencyBeautySigD[0],"IP Step Efficiency","");
3134 legipeff->AddEntry(fEfficiencyBeautySigD[0],"beauty e","p");
0e30407a 3135 legipeff->AddEntry(fEfficiencyBeautySigesdD[0],"beauty e(esd pt)","p");
11ff28c5 3136 legipeff->AddEntry(fEfficiencyCharmSigD[0],"charm e","p");
0e30407a 3137 legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p");
3138 legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p");
3139 //legipeff->AddEntry(fConversionEff[0],"conversion e","p");
3140 //legipeff->AddEntry(fNonHFEEff[0],"Dalitz e","p");
11ff28c5 3141 legipeff->Draw("same");
0e30407a 3142 gPad->SetGrid();
3143 //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps");
8c1c76e9 3144}
3145
3146//____________________________________________________________________________
0e30407a 3147THnSparse* AliHFEspectrum::GetBeautyIPEff(Bool_t isMCpt){
8c1c76e9 3148 //
a8ef1999 3149 // Return beauty electron IP cut efficiency
8c1c76e9 3150 //
3151
a8ef1999 3152 const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
3153 const Int_t nCentralitybinning=11;//number of centrality bins
3154 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
3155 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3156 Int_t ptpr = 0;
3157 Int_t nDim=1; //dimensions of the efficiency weighting grid
3158 if(fBeamType==1)
3159 {
3160 ptpr=1;
3161 nDim=2; //dimensions of the efficiency weighting grid
3162 }
3163 Int_t nBin[1] = {nPtbinning1};
3164 Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
8c1c76e9 3165
8c1c76e9 3166
a8ef1999 3167 THnSparseF *ipcut;
3168 if(fBeamType==0) ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
3169 else ipcut = new THnSparseF("beff", "b IP efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
3170
8c1c76e9 3171 for(Int_t idim = 0; idim < nDim; idim++)
a8ef1999 3172 {
3173 if(nDim==1) ipcut->SetBinEdges(idim, kPtRange);
3174 if(nDim==2)
3175 {
3176 ipcut->SetBinEdges(0, kCentralityRange);
3177 ipcut->SetBinEdges(1, kPtRange);
3178 }
3179 }
8c1c76e9 3180 Double_t pt[1];
a8ef1999 3181 Double_t weight[nCentralitybinning];
0e30407a 3182 Double_t weightErr[nCentralitybinning];
a8ef1999 3183 Double_t contents[2];
3184
3185 for(Int_t a=0;a<11;a++)
3186 {
3187 weight[a] = 1.0;
0e30407a 3188 weightErr[a] = 1.0;
a8ef1999 3189 }
3190
3191
3192 Int_t looppt=nBin[0];
3193 Int_t loopcentr=1;
0e30407a 3194 Int_t ibin[2];
a8ef1999 3195 if(fBeamType==1)
3196 {
3197 loopcentr=nBinPbPb[0];
3198 }
3199
3200
3201 for(int icentr=0; icentr<loopcentr; icentr++)
3202 {
3203 for(int i=0; i<looppt; i++)
3204 {
3205 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
0e30407a 3206 if(isMCpt){
3207 if(fEfficiencyIPBeautyD[icentr]){
3208 weight[icentr]=fEfficiencyIPBeautyD[icentr]->Eval(pt[0]);
3209 weightErr[icentr] = 0;
3210 }
3211 else{
3212 printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
3213 weight[icentr] = fEfficiencyBeautySigD[icentr]->GetBinContent(i+1);
3214 weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
3215 }
3216 }
a8ef1999 3217 else{
0e30407a 3218 if(fEfficiencyIPBeautyesdD[icentr]){
3219 weight[icentr]=fEfficiencyIPBeautyesdD[icentr]->Eval(pt[0]);
3220 weightErr[icentr] = 0;
3221 }
3222 else{
3223 printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
3224 weight[icentr] = fEfficiencyBeautySigesdD[icentr]->GetBinContent(i+1);
3225 weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
3226 }
a8ef1999 3227 }
3228
3229 if(fBeamType==1){
3230 contents[0]=icentr;
3231 contents[1]=pt[0];
0e30407a 3232 ibin[0]=icentr;
3233 ibin[1]=i+1;
a8ef1999 3234 }
3235 if(fBeamType==0){
3236 contents[0]=pt[0];
0e30407a 3237 ibin[0]=i+1;
a8ef1999 3238 }
3239 ipcut->Fill(contents,weight[icentr]);
0e30407a 3240 ipcut->SetBinError(ibin,weightErr[icentr]);
a8ef1999 3241 }
3242 }
3243
a8ef1999 3244 Int_t nDimSparse = ipcut->GetNdimensions();
3245 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3246 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3247
3248 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3249 binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
3250 nBins *= binsvar[iVar];
3251 }
3252
3253 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3254 // loop that sets 0 error in each bin
3255 for (Long_t iBin=0; iBin<nBins; iBin++) {
3256 Long_t bintmp = iBin ;
3257 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3258 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3259 bintmp /= binsvar[iVar] ;
3260 }
0e30407a 3261 //ipcut->SetBinError(binfill,0.); // put 0 everywhere
a8ef1999 3262 }
8c1c76e9 3263
dcef324e 3264 delete[] binsvar;
3265 delete[] binfill;
3266
8c1c76e9 3267 return ipcut;
3268}
3269
3270//____________________________________________________________________________
3271THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
3272 //
3273 // Return PID x IP cut efficiency
3274 //
a8ef1999 3275 const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
3276 const Int_t nCentralitybinning=11;//number of centrality bins
3277 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
3278 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3279 Int_t ptpr = 0;
3280 Int_t nDim=1; //dimensions of the efficiency weighting grid
3281 if(fBeamType==1)
3282 {
3283 ptpr=1;
3284 nDim=2; //dimensions of the efficiency weighting grid
3285 }
3286 Int_t nBin[1] = {nPtbinning1};
3287 Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3288
3289 THnSparseF *pideff;
3290 if(fBeamType==0) pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
3291 else pideff = new THnSparseF("pideff", "PID efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
3292 for(Int_t idim = 0; idim < nDim; idim++)
3293 {
3294
3295 if(nDim==1) pideff->SetBinEdges(idim, kPtRange);
3296 if(nDim==2)
3297 {
3298 pideff->SetBinEdges(0, kCentralityRange);
3299 pideff->SetBinEdges(1, kPtRange);
3300 }
3301 }
8c1c76e9 3302
a8ef1999 3303 Double_t pt[1];
3304 Double_t weight[nCentralitybinning];
11ff28c5 3305 Double_t weightErr[nCentralitybinning];
a8ef1999 3306 Double_t contents[2];
3307
3024f297 3308 for(Int_t a=0;a<nCentralitybinning;a++)
a8ef1999 3309 {
3310 weight[a] = 1.0;
3024f297 3311 weightErr[a] = 1.0;
a8ef1999 3312 }
3313
3314 Int_t looppt=nBin[0];
3315 Int_t loopcentr=1;
11ff28c5 3316 Int_t ibin[2];
a8ef1999 3317 if(fBeamType==1)
3318 {
3319 loopcentr=nBinPbPb[0];
3320 }
3321
3322 for(int icentr=0; icentr<loopcentr; icentr++)
3323 {
3324 Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
3325 for(int i=0; i<looppt; i++)
3326 {
3327 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3328
3329 Double_t tofpideff = 0.;
11ff28c5 3330 Double_t tofpideffesd = 0.;
a8ef1999 3331 if(fEfficiencyTOFPIDD[icentr])
11ff28c5 3332 tofpideff = fEfficiencyTOFPIDD[icentr]->Eval(pt[0]);
a8ef1999 3333 else{
3334 printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
3335 }
11ff28c5 3336 if(fEfficiencyesdTOFPIDD[icentr])
3337 tofpideffesd = fEfficiencyesdTOFPIDD[icentr]->Eval(pt[0]);
3338 else{
3339 printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
3340 }
a8ef1999 3341
3342 //tof pid eff x tpc pid eff x ip cut eff
3343 if(fIPParameterizedEff){
3344 if(source==0) {
11ff28c5 3345 if(fEfficiencyIPCharmD[icentr]){
a8ef1999 3346 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
11ff28c5 3347 weightErr[icentr] = 0;
3348 }
a8ef1999 3349 else{
3350 printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
3351 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
11ff28c5 3352 weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
a8ef1999 3353 }
3354 }
3355 else if(source==2) {
11ff28c5 3356 if(fEfficiencyIPConversionD[icentr]){
3357 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD[icentr]->Eval(pt[0]);
3358 weightErr[icentr] = 0;
3359 }
a8ef1999 3360 else{
3361 printf("Fit failed on conversion IP cut efficiency for centrality %d\n",icentr);
11ff28c5 3362 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1);
3363 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
a8ef1999 3364 }
3365 }
3366 else if(source==3) {
11ff28c5 3367 if(fEfficiencyIPNonhfeD[icentr]){
3368 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD[icentr]->Eval(pt[0]);
3369 weightErr[icentr] = 0;
3370 }
a8ef1999 3371 else{
3372 printf("Fit failed on dalitz IP cut efficiency for centrality %d\n",icentr);
11ff28c5 3373 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1);
3374 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
a8ef1999 3375 }
3376 }
3377 }
3378 else{
11ff28c5 3379 if(source==0){
3380 if(fEfficiencyIPCharmD[icentr]){
3381 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
3382 weightErr[icentr] = 0;
3383 }
3384 else{
3385 printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
3386 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
3387 weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
3388 }
3389 }
3390 else if(source==2){
3391 if(fBeamType==0){
3392 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
3393 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
3394 }
3395 else{
3396 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); // conversion
3397 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
3398 }
3399 }
3400 else if(source==3){
3401 if(fBeamType==0){
3402 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
3403 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
3404 }
3405 else{
3406 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); // Dalitz
3407 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
3408 }
3409 }
a8ef1999 3410 }
3411
3412 if(fBeamType==1){
3413 contents[0]=icentr;
3414 contents[1]=pt[0];
11ff28c5 3415 ibin[0]=icentr;
3416 ibin[1]=i+1;
a8ef1999 3417 }
3418 if(fBeamType==0){
3419 contents[0]=pt[0];
11ff28c5 3420 ibin[0]=i+1;
a8ef1999 3421 }
3422
3423 pideff->Fill(contents,weight[icentr]);
11ff28c5 3424 pideff->SetBinError(ibin,weightErr[icentr]);
a8ef1999 3425 }
3426 }
3427
3428 Int_t nDimSparse = pideff->GetNdimensions();
3429 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3430 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3431
3432 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3433 binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
3434 nBins *= binsvar[iVar];
3435 }
3436
3437 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3438 // loop that sets 0 error in each bin
3439 for (Long_t iBin=0; iBin<nBins; iBin++) {
3440 Long_t bintmp = iBin ;
3441 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3442 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3443 bintmp /= binsvar[iVar] ;
3444 }
a8ef1999 3445 }
8c1c76e9 3446
dcef324e 3447 delete[] binsvar;
3448 delete[] binfill;
8c1c76e9 3449
11ff28c5 3450
8c1c76e9 3451 return pideff;
3452}
a8ef1999 3453
3454//__________________________________________________________________________
3455AliCFDataGrid *AliHFEspectrum::GetRawBspectra2ndMethod(){
3456 //
3457 // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method
3458 //
3459 Int_t ptpr = 0;
3460 Int_t nDim = 1;
3461 if(fBeamType==0)
3462 {
3463 ptpr=0;
3464 }
3465 if(fBeamType==1)
3466 {
3467 ptpr=1;
3468 nDim=2;
3469 }
3470
11ff28c5 3471 const Int_t nPtbinning1 = 18;//number of pt bins, according to new binning
a8ef1999 3472 const Int_t nCentralitybinning=11;//number of centrality bins
11ff28c5 3473 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};
3474
a8ef1999 3475 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3476 Int_t nBin[1] = {nPtbinning1};
3477 Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
a8ef1999 3478
3479 AliCFDataGrid *rawBeautyContainer;
3480 if(fBeamType==0) rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
3481 else rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBinPbPb);
3482 // printf("number of bins= %d\n",bins[0]);
3483
3484
3485
3486
3487 THnSparseF *brawspectra;
3488 if(fBeamType==0) brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin);
3489 else brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBinPbPb);
3490 if(fBeamType==0) brawspectra->SetBinEdges(0, kPtRange);
3491 if(fBeamType==1)
3492 {
3493 // brawspectra->SetBinEdges(0, centralityBins);
3494 brawspectra->SetBinEdges(0, kCentralityRange);
3495 brawspectra->SetBinEdges(1, kPtRange);
3496 }
3497
3498 Double_t pt[1];
3499 Double_t yields= 0.;
3500 Double_t valuesb[2];
3501
3502 //Int_t looppt=nBin[0];
3503 Int_t loopcentr=1;
3504 if(fBeamType==1)
3505 {
3506 loopcentr=nBinPbPb[0];
3507 }
3508
3509 for(int icentr=0; icentr<loopcentr; icentr++)
3510 {
3511
3512 for(int i=0; i<fBSpectrum2ndMethod->GetNbinsX(); i++){
3513 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3514
3515 yields = fBSpectrum2ndMethod->GetBinContent(i+1);
3516
3517 if(fBeamType==1)
3518 {
3519 valuesb[0]=icentr;
3520 valuesb[1]=pt[0];
3521 }
3522 if(fBeamType==0) valuesb[0]=pt[0];
3523 brawspectra->Fill(valuesb,yields);
3524 }
3525 }
3526
3527
3528
3529 Int_t nDimSparse = brawspectra->GetNdimensions();
3530 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3531 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3532
3533 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3534 binsvar[iVar] = brawspectra->GetAxis(iVar)->GetNbins();
3535 nBins *= binsvar[iVar];
3536 }
3537
3538 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3539 // loop that sets 0 error in each bin
3540 for (Long_t iBin=0; iBin<nBins; iBin++) {
3541 Long_t bintmp = iBin ;
3542 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3543 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3544 bintmp /= binsvar[iVar] ;
3545 }
3546 brawspectra->SetBinError(binfill,0.); // put 0 everywhere
3547 }
3548
3549
3550 rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency
3551 TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(ptpr);
3552
3553 new TCanvas;
3554 fBSpectrum2ndMethod->SetMarkerStyle(24);
3555 fBSpectrum2ndMethod->Draw("p");
3556 hRawBeautySpectra->SetMarkerStyle(25);
3557 hRawBeautySpectra->Draw("samep");
dcef324e 3558
11ff28c5 3559 delete[] binfill;
dcef324e 3560 delete[] binsvar;
11ff28c5 3561
a8ef1999 3562 return rawBeautyContainer;
3563}
3564
e17c1f86 3565//__________________________________________________________________________
a8ef1999 3566void AliHFEspectrum::CalculateNonHFEsyst(Int_t centrality){
3567 //
3568 // Calculate non HFE sys
e17c1f86 3569 //
e17c1f86 3570 //
a8ef1999 3571
e17c1f86 3572 if(!fNonHFEsyst)
3573 return;
3574
3575 Double_t evtnorm[1] = {0.0};
dcef324e 3576 if(fNMCbgEvents[0]>0) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
e17c1f86 3577
3578 AliCFDataGrid *convSourceGrid[kElecBgSources][kBgLevels];
3579 AliCFDataGrid *nonHFESourceGrid[kElecBgSources][kBgLevels];
3580
cedf0381 3581 AliCFDataGrid *bgLevelGrid[2][kBgLevels];//for pi0 and eta based errors
e17c1f86 3582 AliCFDataGrid *bgNonHFEGrid[kBgLevels];
3583 AliCFDataGrid *bgConvGrid[kBgLevels];
3584
3585 Int_t stepbackground = 3;
3586 Int_t* bins=new Int_t[1];
cedf0381 3587 const Char_t *bgBase[2] = {"pi0","eta"};
a8ef1999 3588
3589 bins[0]=fConversionEff[centrality]->GetNbinsX();
3590
e17c1f86 3591 AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
3592 AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
3593
a8ef1999 3594 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
e17c1f86 3595 for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
a8ef1999 3596 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 3597 weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
3598 convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
a8ef1999 3599
3600 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 3601 weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
3602 nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
3603 }
a8ef1999 3604
e17c1f86 3605 bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
cedf0381 3606 for(Int_t iSource = 2; iSource < kElecBgSources; iSource++){
e17c1f86 3607 bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
3608 }
cedf0381 3609 if(!fEtaSyst)
3610 bgConvGrid[iLevel]->Add(convSourceGrid[1][iLevel]);
e17c1f86 3611
3612 bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone();
cedf0381 3613 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 3614 bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
3615 }
cedf0381 3616 if(!fEtaSyst)
3617 bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[1][iLevel]);
a8ef1999 3618
cedf0381 3619 bgLevelGrid[0][iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
3620 bgLevelGrid[0][iLevel]->Add(bgNonHFEGrid[iLevel]);
3621 if(fEtaSyst){
3622 bgLevelGrid[1][iLevel] = (AliCFDataGrid*)nonHFESourceGrid[1][iLevel]->Clone();//background for eta source
3623 bgLevelGrid[1][iLevel]->Add(convSourceGrid[1][iLevel]);
3624 }
e17c1f86 3625 }
cedf0381 3626
e17c1f86 3627
cedf0381 3628 //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)
3629 AliCFDataGrid *bgErrorGrid[2][2];//for pions/eta error base, for lower/upper
3630 TH1D* hBaseErrors[2][2];//pi0/eta and lower/upper
3631 for(Int_t iErr = 0; iErr < 2; iErr++){//errors for pi0 and eta base
3632 bgErrorGrid[iErr][0] = (AliCFDataGrid*)bgLevelGrid[iErr][1]->Clone();
3633 bgErrorGrid[iErr][0]->Add(bgLevelGrid[iErr][0],-1.);
3634 bgErrorGrid[iErr][1] = (AliCFDataGrid*)bgLevelGrid[iErr][2]->Clone();
3635 bgErrorGrid[iErr][1]->Add(bgLevelGrid[iErr][0],-1.);
3636
3637 //plot absolute differences between limit yields (upper/lower limit, based on pi0 and eta errors) and best estimate
e17c1f86 3638
cedf0381 3639 hBaseErrors[iErr][0] = (TH1D*)bgErrorGrid[iErr][0]->Project(0);
3640 hBaseErrors[iErr][0]->Scale(-1.);
3641 hBaseErrors[iErr][0]->SetTitle(Form("Absolute %s-based systematic errors from non-HF meson decays and conversions",bgBase[iErr]));
3642 hBaseErrors[iErr][1] = (TH1D*)bgErrorGrid[iErr][1]->Project(0);
3643 if(!fEtaSyst)break;
3644 }
e17c1f86 3645
cedf0381 3646
3647 //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 3648 TH1D *hSpeciesErrors[kElecBgSources-1];
3649 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
cedf0381 3650 if(fEtaSyst && (iSource == 1))continue;
e17c1f86 3651 hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0);
3652 TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0);
3653 hSpeciesErrors[iSource-1]->Add(hNonHFEtemp);
3654 hSpeciesErrors[iSource-1]->Scale(0.3);
3655 }
cedf0381 3656
3657 //Int_t firstBgSource = 0;//if eta systematics are not from scaling
3658 //if(fEtaSyst){firstBgSource = 1;}//source 0 histograms are not filled if eta errors are independently determined!
3659 TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[1]->Clone();
3660 TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[1]->Clone();
3661 TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[1]->Clone();
e17c1f86 3662
3663 TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
3664 TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
3665
3666 for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
cedf0381 3667 Double_t pi0basedErrLow,pi0basedErrUp,etaErrLow,etaErrUp;
3668 pi0basedErrLow = hBaseErrors[0][0]->GetBinContent(iBin);
3669 pi0basedErrUp = hBaseErrors[0][1]->GetBinContent(iBin);
3670 if(fEtaSyst){
3671 etaErrLow = hBaseErrors[1][0]->GetBinContent(iBin);
3672 etaErrUp = hBaseErrors[1][1]->GetBinContent(iBin);
3673 }
3674 else{ etaErrLow = etaErrUp = 0.;}
3675
e17c1f86 3676 Double_t sqrsumErrs= 0;
3677 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
cedf0381 3678 if(fEtaSyst && (iSource == 1))continue;
e17c1f86 3679 Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
3680 sqrsumErrs+=(scalingErr*scalingErr);
3681 }
3682 for(Int_t iErr = 0; iErr < 2; iErr++){
cedf0381 3683 for(Int_t iLevel = 0; iLevel < 2; iLevel++){
3684 hBaseErrors[iErr][iLevel]->SetBinContent(iBin,hBaseErrors[iErr][iLevel]->GetBinContent(iBin)/hBaseErrors[iErr][iLevel]->GetBinWidth(iBin));
3685 }
3686 if(!fEtaSyst)break;
e17c1f86 3687 }
cedf0381 3688 hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+(etaErrUp*etaErrUp)+sqrsumErrs));
3689 hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+(etaErrLow*etaErrLow)+sqrsumErrs));
e17c1f86 3690 hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin));
3691
3692 hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin));
3693 hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin));
3694 }
3695
e17c1f86 3696
3697 TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
3698 cPiErrors->cd();
3699 cPiErrors->SetLogx();
3700 cPiErrors->SetLogy();
cedf0381 3701 hBaseErrors[0][0]->Draw();
3702 //hBaseErrors[0][1]->SetMarkerColor(kBlack);
3703 //hBaseErrors[0][1]->SetLineColor(kBlack);
3704 //hBaseErrors[0][1]->Draw("SAME");
3705 if(fEtaSyst){
3706 hBaseErrors[1][0]->Draw("SAME");
3707 hBaseErrors[1][0]->SetMarkerColor(kBlack);
3708 hBaseErrors[1][0]->SetLineColor(kBlack);
3709 //hBaseErrors[1][1]->SetMarkerColor(13);
3710 //hBaseErrors[1][1]->SetLineColor(13);
3711 //hBaseErrors[1][1]->Draw("SAME");
3712 }
3713 //hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
3714 //hOverallBinScaledErrsUp->SetLineColor(kBlue);
3715 //hOverallBinScaledErrsUp->Draw("SAME");
e17c1f86 3716 hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
3717 hOverallBinScaledErrsLow->SetLineColor(kGreen);
3718 hOverallBinScaledErrsLow->Draw("SAME");
cedf0381 3719 hScalingErrors->SetLineColor(kBlue);
e17c1f86 3720 hScalingErrors->Draw("SAME");
3721
3722 TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
cedf0381 3723 lPiErr->AddEntry(hBaseErrors[0][0],"Lower error from pion error");
3724 //lPiErr->AddEntry(hBaseErrors[0][1],"Upper error from pion error");
3725 if(fEtaSyst){
3726 lPiErr->AddEntry(hBaseErrors[1][0],"Lower error from eta error");
3727 //lPiErr->AddEntry(hBaseErrors[1][1],"Upper error from eta error");
3728 }
e17c1f86 3729 lPiErr->AddEntry(hScalingErrors, "scaling error");
3730 lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics");
cedf0381 3731 //lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
e17c1f86 3732 lPiErr->Draw("SAME");
3733
3734 //Normalize errors
3735 TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone();
3736 TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone();
3737 hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!)
3738 hLowSystScaled->Scale(evtnorm[0]);
3739 TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone();
3740 TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone();
3741 //histograms to be normalized to TGraphErrors
3742 CorrectFromTheWidth(hNormAllSystErrUp);
3743 CorrectFromTheWidth(hNormAllSystErrLow);
3744
3745 TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs");
3746 cNormOvErrs->cd();
3747 cNormOvErrs->SetLogx();
3748 cNormOvErrs->SetLogy();
3749
3750 TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp);
3751 TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow);
3752 gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors");
3753 gOverallSystErrUp->SetMarkerColor(kBlack);
3754 gOverallSystErrUp->SetLineColor(kBlack);
3755 gOverallSystErrLow->SetMarkerColor(kRed);
3756 gOverallSystErrLow->SetLineColor(kRed);
3757 gOverallSystErrUp->Draw("AP");
3758 gOverallSystErrLow->Draw("PSAME");
3759 TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89);
3760 lAllSys->AddEntry(gOverallSystErrLow,"lower","p");
3761 lAllSys->AddEntry(gOverallSystErrUp,"upper","p");
3762 lAllSys->Draw("same");
3763
3764
3765 AliCFDataGrid *bgYieldGrid;
cedf0381 3766 if(fEtaSyst){
3767 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.
3768 }
3769 bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0][0]->Clone();
e17c1f86 3770
3771 TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0);
3772 TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone();
3773 hRelErrUp->Divide(hBgYield);
3774 TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone();
3775 hRelErrLow->Divide(hBgYield);
3776
3777 TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs");
3778 cRelErrs->cd();
3779 cRelErrs->SetLogx();
3780 hRelErrUp->SetTitle("Relative error of non-HFE background yield");
3781 hRelErrUp->Draw();
3782 hRelErrLow->SetLineColor(kBlack);
3783 hRelErrLow->Draw("SAME");
3784
3785 TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95);
3786 lRel->AddEntry(hRelErrUp, "upper");
3787 lRel->AddEntry(hRelErrLow, "lower");
3788 lRel->Draw("SAME");
3789
3790 //CorrectFromTheWidth(hBgYield);
3791 //hBgYield->Scale(evtnorm[0]);
3792
3793
3794 //write histograms/TGraphs to file
3795 TFile *output = new TFile("systHists.root","recreate");
3796
3797 hBgYield->SetName("hBgYield");
3798 hBgYield->Write();
3799 hRelErrUp->SetName("hRelErrUp");
3800 hRelErrUp->Write();
3801 hRelErrLow->SetName("hRelErrLow");
3802 hRelErrLow->Write();
3803 hUpSystScaled->SetName("hOverallSystErrUp");
3804 hUpSystScaled->Write();
3805 hLowSystScaled->SetName("hOverallSystErrLow");
3806 hLowSystScaled->Write();
3807 gOverallSystErrUp->SetName("gOverallSystErrUp");
3808 gOverallSystErrUp->Write();
3809 gOverallSystErrLow->SetName("gOverallSystErrLow");
3810 gOverallSystErrLow->Write();
3811
3812 output->Close();
a8ef1999 3813 delete output;
3814
e17c1f86 3815}
a8ef1999 3816
cedf0381 3817//____________________________________________________________
3818void AliHFEspectrum::UnfoldBG(AliCFDataGrid* const bgsubpectrum){
3819
3820 //
3821 // Unfold backgrounds to check its sanity
3822 //
3823
3824 AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
3825 //AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
3826 if(!mcContainer){
3827 AliError("MC Container not available");
91e50e2b 3828 return;
cedf0381 3829 }
3830
3831 if(!fCorrelation){
3832 AliError("No Correlation map available");
91e50e2b 3833 return;
cedf0381 3834 }
3835
3836 // Data
3837 AliCFDataGrid *dataGrid = 0x0;
3838 dataGrid = bgsubpectrum;
3839
3840 // Guessed
3841 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
3842 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
3843
3844 // Efficiency
3845 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
3846 efficiencyD->CalculateEfficiency(fStepMC+2,fStepTrue);
3847
3848 // Unfold background spectra
3849 Int_t nDim=1;
3850 if(fBeamType==0)nDim = 1;
3851 if(fBeamType==1)nDim = 2;
3852 AliCFUnfolding unfolding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
3853 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
3854 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
3855 if(fSetSmoothing) unfolding.UseSmoothing();
3856 unfolding.Unfold();
3857
3858 // Results
3859 THnSparse* result = unfolding.GetUnfolded();
3860 TCanvas *ctest = new TCanvas("yvonnetest","yvonnetest",1000,600);
3861 if(fBeamType==1)
3862 {
3863 ctest->Divide(2,2);
3864 ctest->cd(1);
3865 result->GetAxis(0)->SetRange(1,1);
3866 TH1D* htest1=(TH1D*)result->Projection(0);
3867 htest1->Draw();
3868 ctest->cd(2);
3869 result->GetAxis(0)->SetRange(1,1);
3870 TH1D* htest2=(TH1D*)result->Projection(1);
3871 htest2->Draw();
3872 ctest->cd(3);
3873 result->GetAxis(0)->SetRange(6,6);
3874 TH1D* htest3=(TH1D*)result->Projection(0);
3875 htest3->Draw();
3876 ctest->cd(4);
3877 result->GetAxis(0)->SetRange(6,6);
3878 TH1D* htest4=(TH1D*)result->Projection(1);
3879 htest4->Draw();
3880
3881 }
3882
3883
3884
3885
3886
3887 TGraphErrors* unfoldedbgspectrumD = Normalize(result);
3888 if(!unfoldedbgspectrumD) {
3889 AliError("Unfolded background spectrum doesn't exist");
3890 }
3891 else{
3892 TFile *file = TFile::Open("unfoldedbgspectrum.root","recreate");
3893 if(fBeamType==0)unfoldedbgspectrumD->Write("unfoldedbgspectrum");
3894
3895 if(fBeamType==1)
3896 {
3897 Int_t centr=1;
3898 result->GetAxis(0)->SetRange(centr,centr);
3899 unfoldedbgspectrumD = Normalize(result,centr-1);
3900 unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr0_20");
3901 centr=6;
3902 result->GetAxis(0)->SetRange(centr,centr);
3903 unfoldedbgspectrumD = Normalize(result,centr-1);
3904 unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr40_80");
3905 }
3906
3907 file->Close();
3908 }
3909}