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