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