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