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