]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliHFEBeautySpectrum.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEBeautySpectrum.cxx
CommitLineData
4bd4fc13 1
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//
26// Author:
27// Raphaelle Bailhache <R.Bailhache@gsi.de>
28// Markus Fasel <M.Fasel@gsi.de>
29//
30#include <TArrayD.h>
31#include <TH1.h>
32#include <TList.h>
33#include <TObjArray.h>
34#include <TROOT.h>
35#include <TCanvas.h>
36#include <TLegend.h>
37#include <TStyle.h>
38#include <TMath.h>
39#include <TAxis.h>
40#include <TGraphErrors.h>
41#include <TFile.h>
42#include <TPad.h>
43#include <TH2D.h>
44#include <TF1.h>
45
46#include "AliPID.h"
47#include "AliCFContainer.h"
48#include "AliCFDataGrid.h"
49#include "AliCFEffGrid.h"
50#include "AliCFGridSparse.h"
51#include "AliCFUnfolding.h"
52#include "AliLog.h"
53
54#include "AliHFEBeautySpectrum.h"
55#include "AliHFEBeautySpectrumQA.h"
56#include "AliHFEcuts.h"
57#include "AliHFEcontainer.h"
58#include "AliHFEtools.h"
59
60ClassImp(AliHFEBeautySpectrum)
61
62//____________________________________________________________
63AliHFEBeautySpectrum::AliHFEBeautySpectrum(const char *name):
64AliHFECorrectSpectrumBase(name),
65 fQA(NULL),
66 fTemporaryObjects(NULL),
67 fBackground(NULL),
68 fEfficiencyTOFPIDD(NULL),
69 fEfficiencyesdTOFPIDD(NULL),
70 fEfficiencyIPCharmD(NULL),
71 fEfficiencyIPBeautyD(NULL),
72 fEfficiencyIPBeautyesdD(NULL),
73 fEfficiencyIPConversionD(NULL),
74 fEfficiencyIPNonhfeD(NULL),
75 fNonHFEbg(NULL),
76 fWeightCharm(NULL),
77 fInclusiveSpectrum(kFALSE),
78 fDumpToFile(kFALSE),
79 fUnSetCorrelatedErrors(kTRUE),
80 fIPanaHadronBgSubtract(kFALSE),
81 fIPanaCharmBgSubtract(kFALSE),
82 fIPanaNonHFEBgSubtract(kFALSE),
83 fIPParameterizedEff(kFALSE),
84 fNonHFEsyst(0),
85 fBeauty2ndMethod(kFALSE),
86 fIPEffCombinedSamples(kTRUE),
87 fNMCEvents(0),
88 fNMCbgEvents(0),
89 fNCentralityBinAtTheEnd(0),
90 fFillMoreCorrelationMatrix(kFALSE),
91 fHadronEffbyIPcut(NULL),
92 fEfficiencyCharmSigD(NULL),
93 fEfficiencyBeautySigD(NULL),
94 fEfficiencyBeautySigesdD(NULL),
95 fConversionEff(NULL),
96 fNonHFEEff(NULL),
97 fCharmEff(NULL),
98 fBeautyEff(NULL),
99 fConversionEffbgc(NULL),
100 fNonHFEEffbgc(NULL),
101 fBSpectrum2ndMethod(NULL),
102 fkBeauty2ndMethodfilename(""),
103 fBeamType(0),
104 fEtaSyst(kTRUE),
105 fDebugLevel(0),
106 fWriteToFile(kFALSE),
107 fUnfoldBG(kFALSE)
108{
109 //
110 // Default constructor
111 //
112
113 fQA = new AliHFEBeautySpectrumQA("AliHFEBeautySpectrumQA");
114
115 for(Int_t k = 0; k < 20; k++){
116 fLowBoundaryCentralityBinAtTheEnd[k] = 0;
117 fHighBoundaryCentralityBinAtTheEnd[k] = 0;
118 }
119 fNMCbgEvents = 0;
120 fEfficiencyTOFPIDD = 0;
121 fEfficiencyesdTOFPIDD = 0;
122 fEfficiencyIPCharmD = 0;
123 fEfficiencyIPBeautyD = 0;
124 fEfficiencyIPBeautyesdD = 0;
125 fEfficiencyIPConversionD = 0;
126 fEfficiencyIPNonhfeD = 0;
127
128 fConversionEff = 0;
129 fNonHFEEff = 0;
130 fCharmEff = 0;
131 fBeautyEff = 0;
132 fEfficiencyCharmSigD = 0;
133 fEfficiencyBeautySigD = 0;
134 fEfficiencyBeautySigesdD = 0;
135
136
137 memset(fConvSourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
138 memset(fNonHFESourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
139}
140//____________________________________________________________
141AliHFEBeautySpectrum::AliHFEBeautySpectrum(const AliHFEBeautySpectrum &ref):
142 AliHFECorrectSpectrumBase(ref),
143 fQA(ref.fQA),
144 fTemporaryObjects(ref.fTemporaryObjects),
145 fBackground(ref.fBackground),
146 fEfficiencyTOFPIDD(ref.fEfficiencyTOFPIDD),
147 fEfficiencyesdTOFPIDD(ref.fEfficiencyesdTOFPIDD),
148 fEfficiencyIPCharmD(ref.fEfficiencyIPCharmD),
149 fEfficiencyIPBeautyD(ref.fEfficiencyIPBeautyD),
150 fEfficiencyIPBeautyesdD(ref.fEfficiencyIPBeautyesdD),
151 fEfficiencyIPConversionD(ref.fEfficiencyIPConversionD),
152 fEfficiencyIPNonhfeD(ref.fEfficiencyIPNonhfeD),
153 fNonHFEbg(ref.fNonHFEbg),
154 fWeightCharm(ref.fWeightCharm),
155 fInclusiveSpectrum(ref.fInclusiveSpectrum),
156 fDumpToFile(ref.fDumpToFile),
157 fUnSetCorrelatedErrors(ref.fUnSetCorrelatedErrors),
158 fIPanaHadronBgSubtract(ref.fIPanaHadronBgSubtract),
159 fIPanaCharmBgSubtract(ref.fIPanaCharmBgSubtract),
160 fIPanaNonHFEBgSubtract(ref.fIPanaNonHFEBgSubtract),
161 fIPParameterizedEff(ref.fIPParameterizedEff),
162 fNonHFEsyst(ref.fNonHFEsyst),
163 fBeauty2ndMethod(ref.fBeauty2ndMethod),
164 fIPEffCombinedSamples(ref.fIPEffCombinedSamples),
165 fNMCEvents(ref.fNMCEvents),
166 fNMCbgEvents(ref.fNMCbgEvents),
167 fNCentralityBinAtTheEnd(ref.fNCentralityBinAtTheEnd),
168 fFillMoreCorrelationMatrix(ref.fFillMoreCorrelationMatrix),
169 fHadronEffbyIPcut(ref.fHadronEffbyIPcut),
170 fEfficiencyCharmSigD(ref.fEfficiencyCharmSigD),
171 fEfficiencyBeautySigD(ref.fEfficiencyBeautySigD),
172 fEfficiencyBeautySigesdD(ref.fEfficiencyBeautySigesdD),
173 fConversionEff(ref.fConversionEff),
174 fNonHFEEff(ref.fNonHFEEff),
175 fCharmEff(ref.fCharmEff),
176 fBeautyEff(ref.fBeautyEff),
177 fConversionEffbgc(ref.fConversionEffbgc),
178 fNonHFEEffbgc(ref.fNonHFEEffbgc),
179 fBSpectrum2ndMethod(ref.fBSpectrum2ndMethod),
180 fkBeauty2ndMethodfilename(ref.fkBeauty2ndMethodfilename),
181 fBeamType(ref.fBeamType),
182 fEtaSyst(ref.fEtaSyst),
183 fDebugLevel(ref.fDebugLevel),
184 fWriteToFile(ref.fWriteToFile),
185 fUnfoldBG(ref.fUnfoldBG)
186{
187 //
188 // Copy constructor
189 //
190
191 ref.Copy(*this);
192}
193//____________________________________________________________
194AliHFEBeautySpectrum &AliHFEBeautySpectrum::operator=(const AliHFEBeautySpectrum &ref){
195 //
196 // Assignment operator
197 //
198 if(this == &ref)
199 ref.Copy(*this);
200 return *this;
201}
202//____________________________________________________________
203void AliHFEBeautySpectrum::Copy(TObject &o) const {
204 //
205 // Copy into object o
206 //
207 AliHFEBeautySpectrum &target = dynamic_cast<AliHFEBeautySpectrum &>(o);
208 target.fQA = fQA;
209
210 target.fQA = fQA;
211 target.fTemporaryObjects = fTemporaryObjects;
212 target.fBackground = fBackground;
213 target.fWeightCharm = fWeightCharm;
214 target.fDumpToFile = fDumpToFile;
215 target.fUnSetCorrelatedErrors = fUnSetCorrelatedErrors;
216 target.fIPanaHadronBgSubtract = fIPanaHadronBgSubtract;
217 target.fIPanaCharmBgSubtract = fIPanaCharmBgSubtract;
218 target.fIPanaNonHFEBgSubtract = fIPanaNonHFEBgSubtract;
219 target.fIPParameterizedEff = fIPParameterizedEff;
220 target.fNonHFEsyst = fNonHFEsyst;
221 target.fBeauty2ndMethod = fBeauty2ndMethod;
222 target.fIPEffCombinedSamples = fIPEffCombinedSamples;
223 target.fNMCEvents = fNMCEvents;
224 target.fNCentralityBinAtTheEnd = fNCentralityBinAtTheEnd;
225 target.fFillMoreCorrelationMatrix = fFillMoreCorrelationMatrix;
226 target.fHadronEffbyIPcut = fHadronEffbyIPcut;
227 target.fConversionEffbgc = fConversionEffbgc;
228 target.fNonHFEEffbgc = fNonHFEEffbgc;
229 target.fBSpectrum2ndMethod = fBSpectrum2ndMethod;
230 target.fkBeauty2ndMethodfilename = fkBeauty2ndMethodfilename;
231 target.fBeamType = fBeamType;
232 target.fEtaSyst = fEtaSyst;
233 target.fDebugLevel = fDebugLevel;
234 target.fWriteToFile = fWriteToFile;
235 target.fUnfoldBG = fUnfoldBG;
236}
237//____________________________________________________________
238AliHFEBeautySpectrum::~AliHFEBeautySpectrum(){
239 //
240 // Destructor
241 //
242 if(fTemporaryObjects){
243 fTemporaryObjects->Clear();
244 delete fTemporaryObjects;
245 }
246 if(fQA) delete fQA;
247}
248//____________________________________________________________
249Bool_t AliHFEBeautySpectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *bghfecontainer, const AliHFEcontainer */*v0hfecontainer*/,AliCFContainer */*photoniccontainerD*/){
250 //
251 // Init what we need for the correction:
252 //
253 // Raw spectrum, hadron contamination
254 // MC efficiency maps, correlation matrix
255 //
256 // This for a given dimension.
257 // If no inclusive spectrum, then take only efficiency map for beauty electron
258 // and the appropriate correlation matrix
259 //
260
261
262 if(fBeauty2ndMethod) CallInputFileForBeauty2ndMethod();
263
264 Int_t kNdim = 3;
265
266 Int_t dims[kNdim];
267
268 switch(fNbDimensions){
269 case 1: dims[0] = 0;
270 break;
271 case 2: for(Int_t i = 0; i < 2; i++) dims[i] = i;
272 break;
273 case 3: for(Int_t i = 0; i < 3; i++) dims[i] = i;
274 break;
275 default:
276 AliError("Container with this number of dimensions not foreseen (yet)");
277 return kFALSE;
278 };
279
280 //
281 // Data container: raw spectrum + hadron contamination
282 //
283 AliCFContainer *datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
284 AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
285 if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
286 AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
287 AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
288 if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
289 SetContainer(datacontainerD,AliHFECorrectSpectrumBase::kDataContainer);
290 SetContainer(contaminationcontainerD,AliHFECorrectSpectrumBase::kBackgroundData);
291
292 // QA : see with MJ if it is necessary for Beauty
293 Int_t dimqa = datacontainer->GetNVar();
294 Int_t dimsqa[dimqa];
295 for(Int_t i = 0; i < dimqa; i++) dimsqa[i] = i;
296 AliCFContainer *datacontainerDQA = GetSlicedContainer(datacontainer, dimqa, dimsqa, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
297 fQA->AddResultAt(datacontainerDQA,AliHFEBeautySpectrumQA::kDataProjection);
298
299 //
300 // MC container: ESD/MC efficiency maps + MC/MC efficiency maps
301 //
302 AliCFContainer *mccontaineresd = 0x0;
303 AliCFContainer *mccontaineresdbg = 0x0;
304 AliCFContainer *mccontainermc = 0x0;
305 AliCFContainer *mccontainermcbg = 0x0;
306 AliCFContainer *nonHFEweightedContainer = 0x0;
307 AliCFContainer *convweightedContainer = 0x0;
308 AliCFContainer *nonHFEweightedContainerSig = 0x0;//mjmj
309 AliCFContainer *convweightedContainerSig = 0x0;//mjmj
310 AliCFContainer *nonHFEtempContainer = 0x0;//temporary container to be sliced for the fnonHFESourceContainers
311 AliCFContainer *convtempContainer = 0x0;//temporary container to be sliced for the fConvSourceContainers
312
313 mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
314 mccontaineresdbg = bghfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
315 mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
316 mccontainermcbg = bghfecontainer->MakeMergedCFContainer("summcbg","summcbg","MCTrackCont:recTrackContMC:recTrackContDEMC");
317
318 if(fNonHFEsyst){
319 const Char_t *sourceName[kElecBgSources]={"Pion","Eta","Omega","Phi","EtaPrime","Rho","K0sSec","OtherSecPi0","Else"};
320 const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
321 for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
322 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
323 if(!(nonHFEtempContainer = bghfecontainer->GetCFContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]))))continue;
324 convtempContainer = bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
325 if(!(fConvSourceContainer[iSource][iLevel] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh)))continue;
326 if(!(fNonHFESourceContainer[iSource][iLevel] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh)))continue;
327 // if((!fConvSourceContainer[iSource][iLevel][icentr])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;
328 if(fBeamType == 1)break;//no systematics yet for PbPb non-HFE backgrounds
329 }
330 }
331 }
332 // else{
333 nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
334 convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
335 if((!convweightedContainer)||(!nonHFEweightedContainer)) return kFALSE;
336 nonHFEweightedContainerSig = mchfecontainer->GetCFContainer("mesonElecs");//mjmj
337 convweightedContainerSig = mchfecontainer->GetCFContainer("conversionElecs");//mjmj
338 if((!convweightedContainerSig)||(!nonHFEweightedContainerSig)) return kFALSE;
339 //}
340
341 if((!mccontaineresd) || (!mccontainermc)) return kFALSE;
342
343 Int_t source = 1; //beauty
344 AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
345 AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
346 if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
347 SetContainer(mccontainermcD,AliHFECorrectSpectrumBase::kMCContainerMC);
348 SetContainer(mccontaineresdD,AliHFECorrectSpectrumBase::kMCContainerESD);
349
350 // set charm, nonHFE container to estimate BG
351 source = 0; //charm
352 mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh); //mjmjmj
353 //mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
354 SetContainer(mccontainermcD,AliHFECorrectSpectrumBase::kMCContainerCharmMC);
355
356 //if(!fNonHFEsyst){
357 AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
358 SetContainer(nonHFEweightedContainerD,AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESD);
359 AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
360 SetContainer(convweightedContainerD,AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESD);
361
362 //mjmj
363 AliCFContainer *nonHFEweightedContainerSigD = GetSlicedContainer(nonHFEweightedContainerSig, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
364 SetContainer(nonHFEweightedContainerSigD,AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig);
365 if(!GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig)){printf("merged non-HFE container not found!\n");return kFALSE;};
366 AliCFContainer *convweightedContainerSigD = GetSlicedContainer(convweightedContainerSig, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
367 SetContainer(convweightedContainerSigD,AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig);
368 if(!GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig)){printf(" merged conversion container not found!\n");return kFALSE;};
369 //}
370
371 SetParameterizedEff(mccontainermc, mccontainermcbg, mccontaineresd, mccontaineresdbg, dims);
372
373 //
374 // MC container: correlation matrix
375 //
376 THnSparseF *mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
377 //mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
378 if(!mccorrelation) return kFALSE;
379 Int_t testCentralityLow = fTestCentralityLow;
380 Int_t testCentralityHigh = fTestCentralityHigh;
381 if(fFillMoreCorrelationMatrix) {
382 testCentralityLow = fTestCentralityLow-1;
383 testCentralityHigh = fTestCentralityHigh+1;
384 }
385 THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims, fChargeChoosen, testCentralityLow,testCentralityHigh);
386 if(!mccorrelationD) {
387 printf("No correlation\n");
388 return kFALSE;
389 }
390 SetCorrelation(mccorrelationD);
391
392 // QA
393 fQA->AddResultAt(mccorrelation,AliHFEBeautySpectrumQA::kCMProjection);
394 //
395 fQA->DrawProjections();
396
397
398 return kTRUE;
399}
400
401
402//____________________________________________________________
403void AliHFEBeautySpectrum::CallInputFileForBeauty2ndMethod(){
404 //
405 // get spectrum for beauty 2nd method
406 //
407 //
408 TFile *inputfile2ndmethod=TFile::Open(fkBeauty2ndMethodfilename);
409 fBSpectrum2ndMethod = new TH1D(*static_cast<TH1D*>(inputfile2ndmethod->Get("BSpectrum")));
410}
411
412//____________________________________________________________
413Bool_t AliHFEBeautySpectrum::Correct(Bool_t subtractcontamination, Bool_t /*subtractphotonic*/){
414 //
415 // Correct the spectrum for efficiency and unfolding for beauty analysis
416 // with both method and compare
417 //
418
419 gStyle->SetPalette(1);
420 gStyle->SetOptStat(1111);
421 gStyle->SetPadBorderMode(0);
422 gStyle->SetCanvasColor(10);
423 gStyle->SetPadLeftMargin(0.13);
424 gStyle->SetPadRightMargin(0.13);
425
426 printf("Steps are: stepdata %d, stepMC %d, steptrue %d\n",fStepData,fStepMC,fStepTrue);
427
428 ///////////////////////////
429 // Check initialization
430 ///////////////////////////
431
432 if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
433 AliInfo("You have to init before");
434 return kFALSE;
435 }
436
437 if((fStepTrue <= 0) && (fStepMC <= 0) && (fStepData <= 0)) {
438 AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
439 return kFALSE;
440 }
441
442 SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
443
444 AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
445 //////////////////////////////////
446 // Subtract hadron background
447 /////////////////////////////////
448 AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
449 AliCFDataGrid *unnormalizedRawSpectrum = 0x0;
450 TGraphErrors *gNormalizedRawSpectrum = 0x0;
451 if(subtractcontamination) {
452 if(!fBeauty2ndMethod) dataspectrumaftersubstraction = SubtractBackground(kTRUE);
453 else dataspectrumaftersubstraction = GetRawBspectra2ndMethod();
454 unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone();
455 dataGridAfterFirstSteps = dataspectrumaftersubstraction;
456 gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum);
457 }
458
459 /////////////////////////////////////////////////////////////////////////////////////////
460 // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
461 /////////////////////////////////////////////////////////////////////////////////////////
462
463 AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
464
465 if(fEfficiencyFunction){
466 dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
467 dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
468 }
469
470 ///////////////
471 // Unfold
472 //////////////
473 THnSparse *correctedspectrum = Unfold(dataGridAfterFirstSteps);
474 if(!correctedspectrum){
475 AliError("No corrected spectrum\n");
476 return kFALSE;
477 }
478
479 /////////////////////
480 // Simply correct
481 ////////////////////
482
483 AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
484
485
486 //////////
487 // Plot
488 //////////
489
490 // QA final results
491 TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
492 TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
493 fQA->AddResultAt(correctedspectrumD,AliHFEBeautySpectrumQA::kFinalResultUnfolded);
494 fQA->AddResultAt(alltogetherspectrumD,AliHFEBeautySpectrumQA::kFinalResultDirectEfficiency);
495 fQA->AddResultAt(correctedspectrum,AliHFEBeautySpectrumQA::kFinalResultUnfSparse);
496 fQA->AddResultAt(alltogetherCorrection,AliHFEBeautySpectrumQA::kFinalResultDirectEffSparse);
497 fQA->DrawResult();
498
499
500 if(fBeamType == 0){
501 if(fNonHFEsyst){
502 CalculateNonHFEsyst();
503 }
504 }
505
506 // Dump to file if needed
507
508 if(fDumpToFile) {
509 TFile *out;
510 out = new TFile("finalSpectrum.root","update");
511 out->cd();
512 // to do centrality dependent
513 correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
514 correctedspectrum->Write();
515 alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
516 alltogetherCorrection->Write();
517 //
518 if(unnormalizedRawSpectrum) {
519 unnormalizedRawSpectrum->SetName("beautyAfterIP");
520 unnormalizedRawSpectrum->Write();
521 }
522
523 if(gNormalizedRawSpectrum){
524 gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP");
525 gNormalizedRawSpectrum->Write();
526 }
527
528 fEfficiencyCharmSigD->SetTitle("IPEfficiencyForCharmSig");
529 fEfficiencyCharmSigD->SetName("IPEfficiencyForCharmSig");
530 fEfficiencyCharmSigD->Write();
531 fEfficiencyBeautySigD->SetTitle("IPEfficiencyForBeautySig");
532 fEfficiencyBeautySigD->SetName("IPEfficiencyForBeautySig");
533 fEfficiencyBeautySigD->Write();
534 fCharmEff->SetTitle("IPEfficiencyForCharm");
535 fCharmEff->SetName("IPEfficiencyForCharm");
536 fCharmEff->Write();
537 fBeautyEff->SetTitle("IPEfficiencyForBeauty");
538 fBeautyEff->SetName("IPEfficiencyForBeauty");
539 fBeautyEff->Write();
540 fConversionEff->SetTitle("IPEfficiencyForConversion");
541 fConversionEff->SetName("IPEfficiencyForConversion");
542 fConversionEff->Write();
543 fNonHFEEff->SetTitle("IPEfficiencyForNonHFE");
544 fNonHFEEff->SetName("IPEfficiencyForNonHFE");
545 fNonHFEEff->Write();
546
547
548 out->Close();
549 delete out;
550 }
551 return kTRUE;
552}
553//____________________________________________________________
554AliCFDataGrid* AliHFEBeautySpectrum::SubtractBackground(Bool_t setBackground){
555 //
556 // Apply background subtraction for IP analysis
557 //
558
559 Int_t nbins=1;
560 printf("subtracting background\n");
561 // Raw spectrum
562 AliCFContainer *dataContainer = GetContainer(kDataContainer);
563 if(!dataContainer){
564 AliError("Data Container not available");
565 return NULL;
566 }
567 printf("Step data: %d\n",fStepData);
568 AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
569
570 AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
571 dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
572
573 // Background Estimate
574 AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
575 if(!backgroundContainer){
576 AliError("MC background container not found");
577 return NULL;
578 }
579
580 Int_t stepbackground = 1;
581 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
582
583 const Int_t bgPlots = 8;
584 TH1D *incElec = 0x0;
585 TH1D *hadron = 0x0;
586 TH1D *charm = 0x0;
587 TH1D *nonHFE[bgPlots];
588 TH1D *subtracted = 0x0;
589
590 THnSparseF* sparseIncElec = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
591 incElec = (TH1D *) sparseIncElec->Projection(0);
afb48e1d 592 AliHFEtools::NormaliseBinWidth(incElec);
4bd4fc13 593
594 TH1D* htemp;
595 Int_t* bins=new Int_t[2];
596
597 if(fIPanaHadronBgSubtract){
598 // Hadron background
599 printf("Hadron background for IP analysis subtracted!\n");
600 htemp = (TH1D *) fHadronEffbyIPcut->Projection(0);
601 bins[0]=htemp->GetNbinsX();
602
603 AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",nbins,bins);
604 hbgContainer->SetGrid(fHadronEffbyIPcut);
605 backgroundGrid->Multiply(hbgContainer,1);
606 // subtract hadron contamination
607 spectrumSubtracted->Add(backgroundGrid,-1.0);
608 // draw raw hadron bg spectra
609 THnSparseF* sparseHbg = (THnSparseF *) hbgContainer->GetGrid();
610 hadron = (TH1D *) sparseHbg->Projection(0);
afb48e1d 611 AliHFEtools::NormaliseBinWidth(hadron);
4bd4fc13 612 }
613
614 if(fIPanaCharmBgSubtract){
615 // Charm background
616 printf("Charm background for IP analysis subtracted!\n");
617 AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
618 spectrumSubtracted->Add(charmbgContainer,-1.0);
619 THnSparseF *sparseCharmElec = (THnSparseF *) charmbgContainer->GetGrid();
620 charm = (TH1D *) sparseCharmElec->Projection(0);
afb48e1d 621 AliHFEtools::NormaliseBinWidth(charm);
4bd4fc13 622 }
623
624 const Char_t *sourceName[bgPlots]={"Conversion","Dalitz","K0s secondary pion Dalitz","K0s secondary pion conversions","Other secondary pion Dalitz","Other secondary pion Dalitz conversions","Other Dalitz", "Other conversions"};
625 const Int_t style[2] = {21,22};
626 const Int_t color[3] = {7,12,8};
627 if(fIPanaNonHFEBgSubtract){
628 Int_t iPlot = 0;
629 for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
630 if((iSource>0)&&(iSource<6))continue;//standard sources are summed up in case of iSource == 0; not treated separately in plots
631 for(Int_t decay = 0; decay < 2; decay++){
632 AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground(decay,iSource);//conversions/Dalitz decays from different sources
633 if(iSource == 0)
634 spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
635 THnSparseF* sparseNonHFEelecs = (THnSparseF *) nonHFEbgContainer->GetGrid();
636 nonHFE[iPlot] = (TH1D *) sparseNonHFEelecs->Projection(0);
afb48e1d 637 AliHFEtools::NormaliseBinWidth(nonHFE[iPlot]);
4bd4fc13 638 iPlot++;
639 }
640 if(!(fNonHFEsyst && (fBeamType == 1)))break;
641 }
642 }
643
644 TLegend *lRaw = new TLegend(0.55,0.55,0.85,0.85);
645 TCanvas *cRaw = new TCanvas("cRaw","cRaw",1000,800);
646 cRaw->cd();
647 gPad->SetLogx();
648 gPad->SetLogy();
649 incElec->GetXaxis()->SetRangeUser(0.4,8.);
650 incElec->Draw("p");
651 incElec->SetMarkerColor(1);
652 incElec->SetMarkerStyle(20);
653 lRaw->AddEntry(incElec,"inclusive electron spectrum");
654
655 if(fIPanaCharmBgSubtract){
656 charm->Draw("samep");
657 charm->SetMarkerColor(3);
658 charm->SetMarkerStyle(20);
659 charm->GetXaxis()->SetRangeUser(0.4,8.);
660 lRaw->AddEntry(charm,"charm elecs");
661 }
662
663 if(fIPanaNonHFEBgSubtract){
664 Int_t iPlot = 0;
665 for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
666 if((iSource>0)&&(iSource<6))continue;//standard sources are summed up in case of iSource == 0; not treated separately in plots
667 for(Int_t decay = 0; decay < 2; decay++){
668 nonHFE[iPlot]->Draw("samep");
669 if(iSource == 0){
670 nonHFE[iPlot]->SetMarkerStyle(20);
671 if(decay == 0)
672 nonHFE[iPlot]->SetMarkerColor(4);
673 else
674 nonHFE[iPlot]->SetMarkerColor(6);
675 }
676 else{
677 nonHFE[iPlot]->SetMarkerColor(color[iSource-6]);
678 nonHFE[iPlot]->SetMarkerStyle(style[decay]);
679 }
680 nonHFE[iPlot]->GetXaxis()->SetRangeUser(0.4,8.);
681 lRaw->AddEntry(nonHFE[iPlot],sourceName[iPlot]);
682 iPlot++;
683 }
684 if(!(fNonHFEsyst && (fBeamType == 1)))break;//only draw standard sources
685 }
686 }
687
688 THnSparseF* sparseSubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
689 subtracted = (TH1D *) sparseSubtracted->Projection(0);
afb48e1d 690 AliHFEtools::NormaliseBinWidth(subtracted);
4bd4fc13 691 subtracted->Draw("samep");
692 subtracted->SetMarkerStyle(24);
693 lRaw->AddEntry(subtracted,"subtracted electron spectrum");
694 lRaw->Draw("SAME");
695
696 delete[] bins;
697
698 TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(0);
afb48e1d 699 AliHFEtools::NormaliseBinWidth(measuredTH1background);
4bd4fc13 700
701 if(setBackground){
702 if(fBackground) delete fBackground;
703 fBackground = backgroundGrid;
704 } else delete backgroundGrid;
705
706 fQA->AddResultAt(subtracted,AliHFEBeautySpectrumQA::kAfterSC);
707 fQA->AddResultAt(incElec,AliHFEBeautySpectrumQA::kBeforeSC);
708 fQA->AddResultAt(measuredTH1background, AliHFEBeautySpectrumQA::kMeasBG);
709 fQA->DrawSubtractContamination();
710
711 return spectrumSubtracted;
712}
713//____________________________________________________________________________________________
714AliCFDataGrid* AliHFEBeautySpectrum::GetNonHFEBackground(Int_t decay, Int_t source){
715 //
716 // calculate non-HF electron background
717 // Arguments: decay = 0 for conversions, decay = 1 for meson decays to electrons
718 // source: gives mother either of the electron (for meson->electron) or of the gamma (conversions);
719 // source numbers as given by array in AliAnalysisTaskHFE:
720 // {"Pion","Eta","Omega","Phi","EtaPrime","Rho","K0sSec","OtherSecPi0","Else"};
721 //
722
723 Double_t evtnorm[1] = {1.};
724 if(fNMCbgEvents) evtnorm[0]= double(fNEvents)/double(fNMCbgEvents);
725
726 Int_t stepbackground = 1;//for non-standard background, step after IP cut is taken directly
727 if(source == 0) stepbackground = 3;//for standard sources, take step before PID -> correction of PID x IP efficiencies from enhanced samples
728 AliCFContainer *backgroundContainer = 0x0;
729 if(fNonHFEsyst){
730 if(decay == 0){//conversions
731 backgroundContainer = (AliCFContainer*)fConvSourceContainer[source][0]->Clone();
732 if(source == 0)//standard conversion containers
733 for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){//add all standard sources
734 backgroundContainer->Add(fConvSourceContainer[iSource][0]);
735 }
736 }
737 else if(decay == 1){//Dalitz decays, same procedure as above
738 backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[source][0]->Clone();
739 if(source == 0)
740 for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
afb48e1d 741 if(iSource == 1)
742 backgroundContainer->Add(fNonHFESourceContainer[iSource][0],1.41);//correction for the eta Dalitz decay branching ratio in PYTHIA
743 else
744 backgroundContainer->Add(fNonHFESourceContainer[iSource][0]);
4bd4fc13 745 }
746 }
747 }
748 else{//careful: these containers include the non-standard electron sources (K0s, etc.). If necessary, use SetRange in species axis to fix it?
749 if(source == 0){
750 // Background Estimate
751 if(decay == 0)
752 backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
753 else if(decay == 1)
754 backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
755 }
756 }
757 if(!backgroundContainer){
758 AliError("MC background container not found");
759 return NULL;
760 }
761 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("backgroundGrid","backgroundGrid",*backgroundContainer,stepbackground);
762
763 Double_t rangelow = 1.;
764 Double_t rangeup = 6.;
765 if(decay == 1) rangelow = 0.9;
766
767
768 const Char_t *dmode[2]={"Conversions","Dalitz decays"};
769 TF1 *fitHagedorn = new TF1("fitHagedorn", "[0]/TMath::Power(([1]+x/[2]), [3])", rangelow, rangeup);
770 fNonHFEbg = 0x0;
771 TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
772 TAxis *axis = h1->GetXaxis();
afb48e1d 773 AliHFEtools::NormaliseBinWidth(h1);
4bd4fc13 774 if(source == 0){
775 fitHagedorn->SetParameter(0, 0.15);
776 fitHagedorn->SetParameter(1, 0.09);
777 fitHagedorn->SetParameter(2, 8.4);
778 fitHagedorn->SetParameter(3, 6.3);
779 // TCanvas *ch1conversion = new TCanvas("ch1conversion","ch1conversion",500,400);
780 // ch1conversion->cd();
781 //fHagedorn->SetLineColor(2);
782 h1->Fit(fitHagedorn,"R");
783 // h1->Draw();
784 if(!(fNonHFEbg = h1->GetFunction("fitHagedorn")))printf("electron background fit failed for %s\n",dmode[decay]);
785 }
786
787 Int_t *nBinpp=new Int_t[1];
788 Int_t *binspp=new Int_t[1];
789 binspp[0]=kSignalPtBins;// number of pt bins
790
791 Int_t looppt=binspp[0];
792
793 for(Long_t iBin=1; iBin<= looppt;iBin++){
794 Double_t iipt= h1->GetBinCenter(iBin);
795 Double_t iiwidth= axis->GetBinWidth(iBin);
796 nBinpp[0]=iBin;
797 Double_t fitFactor = backgroundGrid->GetElement(nBinpp);//if no fit available, just take bin-by-bin information
798 if(fNonHFEbg)fitFactor = fNonHFEbg->Eval(iipt)*iiwidth;
799 backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
800 backgroundGrid->SetElement(nBinpp,fitFactor*evtnorm[0]);
801 }
4bd4fc13 802 //end of workaround for statistical errors
803 if(source == 0){
804 AliCFDataGrid *weightedBGContainer = 0x0;
805 weightedBGContainer = new AliCFDataGrid("weightedBGContainer","weightedBGContainer",1,binspp);
806 if(decay == 0)
807 weightedBGContainer->SetGrid(GetPIDxIPEff(2));
808 else if(decay == 1)
809 weightedBGContainer->SetGrid(GetPIDxIPEff(3));
810 if(stepbackground == 3){
811 backgroundGrid->Multiply(weightedBGContainer,1.0);
812 }
32c709ed 813 }
814 delete[] nBinpp;
815 delete[] binspp;
4bd4fc13 816 return backgroundGrid;
817}
818//____________________________________________________________
819AliCFDataGrid* AliHFEBeautySpectrum::GetCharmBackground(){
820 //
821 // calculate charm background; when using centrality binning 2: centBin = 0 for 0-20%, centBin = 5 for 40-80%
822 //
823
824 Int_t nDim = 1;
825
826 Double_t evtnorm=0;
827 if(fNMCEvents) evtnorm= double(fNEvents)/double(fNMCEvents); //mjmjmj
828 //if(fNMCbgEvents) evtnorm= double(fNEvents)/double(fNMCbgEvents);
829 printf("events for data: %d",fNEvents);
830 printf("events for MC: %d",fNMCEvents);
831 printf("normalization factor for charm: %f",evtnorm);
832
833 AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
834 if(!mcContainer){
835 AliError("MC Container not available");
836 return NULL;
837 }
838
839 if(!fCorrelation){
840 AliError("No Correlation map available");
841 return NULL;
842 }
843
844 AliCFDataGrid *charmBackgroundGrid= 0x0;
845 charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
846 TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
847 Int_t* bins=new Int_t[2];
848 bins[0]=charmbgaftertofpid->GetNbinsX();
849
850 AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins);
851 ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
852 TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(0);
853
854 charmBackgroundGrid->Multiply(ipWeightedCharmContainer,1.);
855 Int_t *nBinpp=new Int_t[1];
856 Int_t* binspp=new Int_t[1];
857 binspp[0]=charmbgaftertofpid->GetNbinsX(); // number of pt bins
858
859 Int_t looppt=binspp[0];
860
861 for(Long_t iBin=1; iBin<= looppt;iBin++){
862 nBinpp[0]=iBin;
863 charmBackgroundGrid->SetElementError(nBinpp, charmBackgroundGrid->GetElementError(nBinpp)*evtnorm);
864 charmBackgroundGrid->SetElement(nBinpp,charmBackgroundGrid->GetElement(nBinpp)*evtnorm);
865 }
866
867 TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(0);
868
869 AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
870 weightedCharmContainer->SetGrid(GetCharmWeights(fTestCentralityLow)); // get charm weighting factors
871 TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(0);
872 charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
873 TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(0);
874
875 // Efficiency (set efficiency to 1 for only folding)
876 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
877 efficiencyD->CalculateEfficiency(0,0);
878
879 // Folding
880 if(fBeamType==0)nDim = 1;
881 AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
882 folding.SetMaxNumberOfIterations(1);
883 folding.Unfold();
884
885 // Results
886 THnSparse* result1= folding.GetEstMeasured(); // folded spectra
887 THnSparse* result=(THnSparse*)result1->Clone();
888 charmBackgroundGrid->SetGrid(result);
889 TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(0);
890
891 //Charm background evaluation plots
892
893 TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
894 cCharmBgEval->Divide(3,1);
895
896 cCharmBgEval->cd(1);
897
898 if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm);
899
afb48e1d 900 AliHFEtools::NormaliseBinWidth(charmbgaftertofpid);
4bd4fc13 901 charmbgaftertofpid->SetMarkerStyle(25);
902 charmbgaftertofpid->Draw("p");
903 charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
904 charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
905 gPad->SetLogy();
906
afb48e1d 907 AliHFEtools::NormaliseBinWidth(charmbgafteripcut);
4bd4fc13 908 charmbgafteripcut->SetMarkerStyle(24);
909 charmbgafteripcut->Draw("samep");
910
afb48e1d 911 AliHFEtools::NormaliseBinWidth(charmbgafterweight);
4bd4fc13 912 charmbgafterweight->SetMarkerStyle(24);
913 charmbgafterweight->SetMarkerColor(4);
914 charmbgafterweight->Draw("samep");
915
afb48e1d 916 AliHFEtools::NormaliseBinWidth(charmbgafterfolding);
4bd4fc13 917 charmbgafterfolding->SetMarkerStyle(24);
918 charmbgafterfolding->SetMarkerColor(2);
919 charmbgafterfolding->Draw("samep");
920
921 cCharmBgEval->cd(2);
922 parametrizedcharmpidipeff->SetMarkerStyle(24);
923 parametrizedcharmpidipeff->Draw("p");
924 parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
925
926 cCharmBgEval->cd(3);
927 charmweightingfc->SetMarkerStyle(24);
928 charmweightingfc->Draw("p");
929 charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
930 charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
931
932 cCharmBgEval->cd(1);
933 TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
934 legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
935 legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
936 legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
937 legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
938 legcharmbg->Draw("same");
939
940 cCharmBgEval->cd(2);
941 TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
942 legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
943 legcharmbg2->Draw("same");
944
945 CorrectStatErr(charmBackgroundGrid);
946 if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
947
948 delete[] bins;
949 delete[] nBinpp;
950 delete[] binspp;
951
952 if(fUnfoldBG) UnfoldBG(charmBackgroundGrid);
953
954 return charmBackgroundGrid;
955}
956//____________________________________________________________
957AliCFDataGrid *AliHFEBeautySpectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
958 //
959 // Apply TPC pid efficiency correction from parametrisation
960 //
961
962 // Data in the right format
963 AliCFDataGrid *dataGrid = 0x0;
964 if(bgsubpectrum) {
965 dataGrid = bgsubpectrum;
966 }
967 else {
968
969 AliCFContainer *dataContainer = GetContainer(kDataContainer);
970 if(!dataContainer){
971 AliError("Data Container not available");
972 return NULL;
973 }
974 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
975 }
976 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
977 result->SetName("ParametrizedEfficiencyBefore");
978 THnSparse *h = result->GetGrid();
979 Int_t nbdimensions = h->GetNdimensions();
980 //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
981 AliCFContainer *dataContainer = GetContainer(kDataContainer);
982 if(!dataContainer){
983 AliError("Data Container not available");
984 return NULL;
985 }
986 AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
987 dataContainerbis->Add(dataContainerbis,-1.0);
988
989 Int_t* coord = new Int_t[nbdimensions];
990 memset(coord, 0, sizeof(Int_t) * nbdimensions);
991 Double_t* points = new Double_t[nbdimensions];
992
993 ULong64_t nEntries = h->GetNbins();
994 for (ULong64_t i = 0; i < nEntries; ++i) {
995 Double_t value = h->GetBinContent(i, coord);
996 //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
997 //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
998
999 // Get the bin co-ordinates given an coord
1000 for (Int_t j = 0; j < nbdimensions; ++j)
1001 points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1002
1003 if (!fEfficiencyFunction->IsInside(points))
1004 continue;
1005 TF1::RejectPoint(kFALSE);
1006
1007 // Evaulate function at points
1008 Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1009 //printf("Value efficiency is %f\n",valueEfficiency);
1010
1011 if(valueEfficiency > 0.0) {
1012 h->SetBinContent(coord,value/valueEfficiency);
1013 dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1014 }
1015 Double_t error = h->GetBinError(i);
1016 h->SetBinError(coord,error/valueEfficiency);
1017 dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1018 }
1019
1020 delete[] coord;
1021 delete[] points;
1022
1023 AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1024
1025// QA
1026 TH1D *afterE = (TH1D *) resultt->Project(0);
afb48e1d 1027 AliHFEtools::NormaliseBinWidth(afterE);
4bd4fc13 1028 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
afb48e1d 1029 AliHFEtools::NormaliseBinWidth(beforeE);
4bd4fc13 1030 fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterPE);
1031 fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforePE);
1032 fQA->AddResultAt(fEfficiencyFunction,AliHFEBeautySpectrumQA::kPEfficiency);
1033 fQA->DrawCorrectWithEfficiency(AliHFEBeautySpectrumQA::kParametrized);
1034
1035 return resultt;
1036}
1037//____________________________________________________________
1038THnSparse *AliHFEBeautySpectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
1039
1040 //
1041 // Unfold and eventually correct for efficiency the bgsubspectrum
1042 //
1043
1044 AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
1045 if(!mcContainer){
1046 AliError("MC Container not available");
1047 return NULL;
1048 }
1049
1050 if(!fCorrelation){
1051 AliError("No Correlation map available");
1052 return NULL;
1053 }
1054
1055 // Data
1056 AliCFDataGrid *dataGrid = 0x0;
1057 if(bgsubpectrum) {
1058 dataGrid = bgsubpectrum;
1059 }
1060 else {
1061
1062 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1063 if(!dataContainer){
1064 AliError("Data Container not available");
1065 return NULL;
1066 }
1067
1068 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1069 }
1070
1071 // Guessed
1072 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
1073 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
1074
1075 // Efficiency
1076 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1077 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1078
1079 if(!fBeauty2ndMethod)
1080 {
1081 // Consider parameterized IP cut efficiency
1082 Int_t* bins=new Int_t[1];
1083 bins[0]=kSignalPtBins;
1084
1085 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1086 beffContainer->SetGrid(GetBeautyIPEff(kTRUE));
1087 efficiencyD->Multiply(beffContainer,1);
1088 }
1089
1090
1091 // Unfold
1092
1093 AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse,1.e-06,0,fNumberOfIterations);//check with MinJung if last arguments are correct (taken from inclusive analysis...
1094 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
1095 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
1096 if(fSetSmoothing) unfolding.UseSmoothing();
1097 unfolding.Unfold();
1098
1099 // Results
1100 THnSparse* result = unfolding.GetUnfolded();
1101 THnSparse* residual = unfolding.GetEstMeasured();
1102
1103 // QA
1104 TH1D *residualh = (TH1D *) residual->Projection(0);
1105 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1106 TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
1107 TH1D *afterE = (TH1D *) result->Projection(0);
1108
afb48e1d 1109 AliHFEtools::NormaliseBinWidth(residualh);
1110 AliHFEtools::NormaliseBinWidth(beforeE);
1111 AliHFEtools::NormaliseBinWidth(afterE);
4bd4fc13 1112 fQA->AddResultAt(residualh,AliHFEBeautySpectrumQA::kResidualU);
1113 fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterU);
1114 fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforeU);
1115 fQA->AddResultAt(efficiencyDproj,AliHFEBeautySpectrumQA::kUEfficiency);
1116 fQA->DrawUnfolding();
1117
1118 return (THnSparse *) result->Clone();
1119}
1120
1121//____________________________________________________________
1122AliCFDataGrid *AliHFEBeautySpectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
1123
1124 //
1125 // Apply unfolding and efficiency correction together to bgsubspectrum
1126 //
1127
1128 AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
1129 if(!mcContainer){
1130 AliError("MC Container not available");
1131 return NULL;
1132 }
1133
1134 // Efficiency
1135 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1136 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1137
1138
1139 if(!fBeauty2ndMethod)
1140 {
1141 // Consider parameterized IP cut efficiency
1142 Int_t* bins=new Int_t[1];
1143 bins[0]=kSignalPtBins;
1144
1145 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1146 beffContainer->SetGrid(GetBeautyIPEff(kFALSE));
1147 efficiencyD->Multiply(beffContainer,1);
1148 }
1149
1150 // Data in the right format
1151 AliCFDataGrid *dataGrid = 0x0;
1152 if(bgsubpectrum) {
1153 dataGrid = bgsubpectrum;
1154 }
1155 else {
1156
1157 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1158 if(!dataContainer){
1159 AliError("Data Container not available");
1160 return NULL;
1161 }
1162
1163 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1164 }
1165
1166 // Correct
1167 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1168 result->ApplyEffCorrection(*efficiencyD);
1169
1170 // QA
1171 TH1D *afterE = (TH1D *) result->Project(0);
afb48e1d 1172 AliHFEtools::NormaliseBinWidth(afterE);
4bd4fc13 1173 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
afb48e1d 1174 AliHFEtools::NormaliseBinWidth(beforeE);
4bd4fc13 1175 TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
1176 fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterMCE);
1177 fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforeMCE);
1178 fQA->AddResultAt(efficiencyDproj,AliHFEBeautySpectrumQA::kMCEfficiency);
1179 fQA->DrawCorrectWithEfficiency(AliHFEBeautySpectrumQA::kMC);
1180
1181 return result;
1182
1183}
1184//____________________________________________________________
1185void AliHFEBeautySpectrum::AddTemporaryObject(TObject *o){
1186 //
1187 // Emulate garbage collection on class level
1188 //
1189 if(!fTemporaryObjects) {
1190 fTemporaryObjects= new TList;
1191 fTemporaryObjects->SetOwner();
1192 }
1193 fTemporaryObjects->Add(o);
1194}
1195
1196//____________________________________________________________
1197void AliHFEBeautySpectrum::ClearObject(TObject *o){
1198 //
1199 // Do a safe deletion
1200 //
1201 if(fTemporaryObjects){
1202 if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
1203 delete o;
1204 } else{
1205 // just delete
1206 delete o;
1207 }
1208}
1209//_________________________________________________________________________
1210THnSparse* AliHFEBeautySpectrum::GetCharmWeights(Int_t centBin){
1211
1212 //
1213 // Measured D->e based weighting factors
1214 //
1215
1216 const Int_t nDimpp=1;
1217 Int_t nBinpp[nDimpp] = {kSignalPtBins};
1218 Double_t ptbinning1[kSignalPtBins+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.};
1219 Int_t looppt=nBinpp[0];
1220
1221 fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
1222 fWeightCharm->SetBinEdges(0, ptbinning1);
1223
1224 // //if(fBeamType == 0){// Weighting factor for pp
1225 //Double_t weightpp[kSignalPtBins]={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};
1226//TPC+TOF standard cut 4800
1227 Double_t weightpp[kSignalPtBins]={0.050326, 0.045826, 0.042043, 0.039641, 0.039872, 0.041796, 0.044320, 0.046273, 0.047376, 0.047657, 0.047973, 0.048307, 0.045325, 0.044067, 0.043367, 0.040417, 0.037048, 0.033695, 0.032192, 0.029270, 0.027270, 0.024451, 0.020846, 0.019032, 0.018210, 0.017554, 0.015604, 0.015194, 0.013542, 0.013447, 0.015160, 0.015507, 0.014989, 0.012533, 0.025603};
1228 //}
1229 //else{
1230 //if(centBin == 0){
1231 // Weighting factor for PbPb (0-20%)
1232 Double_t weightPbPb1[kSignalPtBins]={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};
1233 //}
1234 //else{
1235 // Weighting factor for PbPb (40-80%)
1236 Double_t weightPbPb2[kSignalPtBins]={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};
1237 // }
1238 //}
1239 Double_t weight[kSignalPtBins];
1240 for(Int_t ipt = 0; ipt < kSignalPtBins; ipt++){
1241 if(fBeamType == 0)weight[ipt] = weightpp[ipt];
1242 else if(centBin == 0)weight[ipt] = weightPbPb1[ipt];
1243 else weight[ipt] = weightPbPb2[ipt];
1244 }
1245
1246 //points
1247 Double_t pt[1];
1248 Double_t contents[2];
1249
1250 for(int i=0; i<looppt; i++){
1251 pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
1252 contents[0]=pt[0];
1253 fWeightCharm->Fill(contents,weight[i]);
1254 }
1255
1256
1257 Int_t nDimSparse = fWeightCharm->GetNdimensions();
1258 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1259 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1260
1261 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1262 binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
1263 nBins *= binsvar[iVar];
1264 }
1265
1266 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1267 // loop that sets 0 error in each bin
1268 for (Long_t iBin=0; iBin<nBins; iBin++) {
1269 Long_t bintmp = iBin ;
1270 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1271 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1272 bintmp /= binsvar[iVar] ;
1273 }
1274 fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
1275 }
1276 delete[] binsvar;
1277 delete[] binfill;
1278
1279 return fWeightCharm;
1280}
1281//____________________________________________________________________________
1282void AliHFEBeautySpectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
1283
1284 // TOF PID efficiencies
1285
1286 TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1287 TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,6.);
1288 TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.0);
1289
1290 if(fBeamType == 1){//not in use yet - adapt as soon as possible!
1291 fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1292 fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1293 fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
1294 }
1295
1296 TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600);
1297 cefficiencyParamtof->cd();
1298
1299 AliCFContainer *mccontainermcD = 0x0;
1300 AliCFContainer *mccontaineresdD = 0x0;
1301 TH1D* efficiencysigTOFPIDD;
1302 TH1D* efficiencyTOFPIDD;
1303 TH1D* efficiencysigesdTOFPIDD;
1304 TH1D* efficiencyesdTOFPIDD;
1305 Int_t source = -1; //get parameterized TOF PID efficiencies
1306
1307 // signal sample
1308 mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1309 AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
1310 efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1311
1312 // mb sample for double check
1313 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1314 AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
1315 efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1316
1317 // mb sample with reconstructed variables
1318 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1319 AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD);
1320 efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1321
1322 // mb sample with reconstructed variables
1323 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1324 AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
1325 efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1326
1327 //fill histo
1328 efficiencysigTOFPIDD = (TH1D *) efficiencymcsigParamTOFPID->Project(0);
1329 efficiencyTOFPIDD = (TH1D *) efficiencymcParamTOFPID->Project(0);
1330 efficiencysigesdTOFPIDD = (TH1D *) efficiencysigesdParamTOFPID->Project(0);
1331 efficiencyesdTOFPIDD = (TH1D *) efficiencyesdParamTOFPID->Project(0);
1332 efficiencysigTOFPIDD->SetName("efficiencysigTOFPIDD");
1333 efficiencyTOFPIDD->SetName("efficiencyTOFPIDD");
1334 efficiencysigesdTOFPIDD->SetName("efficiencysigesdTOFPIDD");
1335 efficiencyesdTOFPIDD->SetName("efficiencyesdTOFPIDD");
1336
1337 //fit (mc enhenced sample)
1338 fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1339 efficiencysigTOFPIDD->Fit(fittofpid,"R");
1340 efficiencysigTOFPIDD->GetYaxis()->SetTitle("Efficiency");
1341 fEfficiencyTOFPIDD = efficiencysigTOFPIDD->GetFunction("fittofpid");
1342
1343 //fit (esd enhenced sample)
1344 efficiencysigesdTOFPIDD->Fit(fittofpid,"R");
1345 efficiencysigesdTOFPIDD->GetYaxis()->SetTitle("Efficiency");
1346 fEfficiencyesdTOFPIDD = efficiencysigesdTOFPIDD->GetFunction("fittofpid");
1347
1348 // draw (for PbPb, only 1st bin)
1349 //sig mc
1350 efficiencysigTOFPIDD->SetTitle("");
1351 efficiencysigTOFPIDD->SetStats(0);
1352 efficiencysigTOFPIDD->SetMarkerStyle(25);
1353 efficiencysigTOFPIDD->SetMarkerColor(2);
1354 efficiencysigTOFPIDD->SetLineColor(2);
1355 efficiencysigTOFPIDD->Draw();
1356
1357 //mb mc
1358 efficiencyTOFPIDD->SetTitle("");
1359 efficiencyTOFPIDD->SetStats(0);
1360 efficiencyTOFPIDD->SetMarkerStyle(24);
1361 efficiencyTOFPIDD->SetMarkerColor(4);
1362 efficiencyTOFPIDD->SetLineColor(4);
1363 efficiencyTOFPIDD->Draw("same");
1364
1365 //sig esd
1366 efficiencysigesdTOFPIDD->SetTitle("");
1367 efficiencysigesdTOFPIDD->SetStats(0);
1368 efficiencysigesdTOFPIDD->SetMarkerStyle(25);
1369 efficiencysigesdTOFPIDD->SetMarkerColor(3);
1370 efficiencysigesdTOFPIDD->SetLineColor(3);
1371 efficiencysigesdTOFPIDD->Draw("same");
1372
1373 //mb esd
1374 efficiencyesdTOFPIDD->SetTitle("");
1375 efficiencyesdTOFPIDD->SetStats(0);
1376 efficiencyesdTOFPIDD->SetMarkerStyle(25);
1377 efficiencyesdTOFPIDD->SetMarkerColor(1);
1378 efficiencyesdTOFPIDD->SetLineColor(1);
1379 efficiencyesdTOFPIDD->Draw("same");
1380
1381 //signal mc fit
1382 if(fEfficiencyTOFPIDD){
1383 fEfficiencyTOFPIDD->SetLineColor(2);
1384 fEfficiencyTOFPIDD->Draw("same");
1385 }
1386 //mb esd fit
1387 if(fEfficiencyesdTOFPIDD){
1388 fEfficiencyesdTOFPIDD->SetLineColor(3);
1389 fEfficiencyesdTOFPIDD->Draw("same");
1390 }
1391
1392 TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
1393 legtofeff->AddEntry(efficiencysigTOFPIDD,"TOF PID Step Efficiency","");
1394 legtofeff->AddEntry(efficiencysigTOFPIDD,"vs MC p_{t} for enhenced samples","p");
1395 legtofeff->AddEntry(efficiencyTOFPIDD,"vs MC p_{t} for mb samples","p");
1396 legtofeff->AddEntry(efficiencysigesdTOFPIDD,"vs esd p_{t} for enhenced samples","p");
1397 legtofeff->AddEntry(efficiencyesdTOFPIDD,"vs esd p_{t} for mb samples","p");
1398 legtofeff->Draw("same");
1399
1400
1401 TCanvas * cefficiencyParamIP = new TCanvas("efficiencyParamIP","efficiencyParamIP",500,500);
1402 cefficiencyParamIP->cd();
1403 gStyle->SetOptStat(0);
1404
1405 // IP cut efficiencies
1406 AliCFContainer *charmCombined = 0x0;
1407 AliCFContainer *beautyCombined = 0x0;
1408 AliCFContainer *beautyCombinedesd = 0x0;
1409
1410 source = 0; //charm enhenced
1411 mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1412 AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
1413 efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1414
1415 charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");
1416
1417 source = 1; //beauty enhenced
1418 mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1419 AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
1420 efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1421
1422 beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined");
1423
1424 mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1425 AliCFEffGrid* efficiencyBeautySigesd = new AliCFEffGrid("efficiencyBeautySigesd","",*mccontaineresdD);
1426 efficiencyBeautySigesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1427
1428 beautyCombinedesd = (AliCFContainer*)mccontaineresdD->Clone("beautyCombinedesd");
1429
1430 source = 0; //charm mb
1431 mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1432 AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
1433 efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1434
1435 charmCombined->Add(mccontainermcD);
1436 AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined);
1437 efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
1438
1439 source = 1; //beauty mb
1440 mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1441 AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
1442 efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1443
1444 beautyCombined->Add(mccontainermcD);
1445 AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined);
1446 efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
1447
1448 mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1449 AliCFEffGrid* efficiencyBeautyesd = new AliCFEffGrid("efficiencyBeautyesd","",*mccontaineresdD);
1450 efficiencyBeautyesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1451
1452 beautyCombinedesd->Add(mccontaineresdD);
1453 AliCFEffGrid* efficiencyBeautyCombinedesd = new AliCFEffGrid("efficiencyBeautyCombinedesd","",*beautyCombinedesd);
1454 efficiencyBeautyCombinedesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
1455
1456 source = 2; //conversion mb
1457 mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1458 AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
1459 efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1460
1461 source = 3; //non HFE except for the conversion mb
1462 mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1463 AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
1464 efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1465
1466 if(fIPEffCombinedSamples){printf("combined samples taken for beauty and charm\n");
1467 fEfficiencyCharmSigD = (TH1D*)efficiencyCharmCombined->Project(0); //signal enhenced + mb
1468 fEfficiencyBeautySigD = (TH1D*)efficiencyBeautyCombined->Project(0); //signal enhenced + mb
1469 fEfficiencyBeautySigesdD = (TH1D*)efficiencyBeautyCombinedesd->Project(0); //signal enhenced + mb
1470 }
1471 else{printf("signal enhanced samples taken for beauty and charm\n");
1472 fEfficiencyCharmSigD = (TH1D*)efficiencyCharmSig->Project(0); //signal enhenced only
1473 fEfficiencyBeautySigD = (TH1D*)efficiencyBeautySig->Project(0); //signal enhenced only
1474 fEfficiencyBeautySigesdD = (TH1D*)efficiencyBeautySigesd->Project(0); //signal enhenced only
1475 }
1476 fCharmEff = (TH1D*)efficiencyCharm->Project(0); //mb only
1477 fBeautyEff = (TH1D*)efficiencyBeauty->Project(0); //mb only
1478 fConversionEff = (TH1D*)efficiencyConv->Project(0); //mb only
1479 fNonHFEEff = (TH1D*)efficiencyNonhfe->Project(0); //mb only
1480
1481 if(fBeamType==0){
1482 //AliCFEffGrid *nonHFEEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
32c709ed 1483 AliCFContainer *cfcontainer = GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig);
1484 if(!cfcontainer) return;
1485 AliCFEffGrid *nonHFEEffGrid = (AliCFEffGrid*) GetEfficiency(cfcontainer,1,0); //mjmj
4bd4fc13 1486 fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
1487
1488 //AliCFEffGrid *conversionEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
32c709ed 1489 AliCFContainer *cfcontainerr = GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig);
1490 if(!cfcontainerr) return;
1491 AliCFEffGrid *conversionEffGrid = (AliCFEffGrid*) GetEfficiency(cfcontainerr,1,0); //mjmj
4bd4fc13 1492 fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
1493
1494 //MHMH
1495 if(fBeamType == 0){
1496 fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1497 fNonHFEEffbgc->Fit(fipfitnonhfe,"R");
1498 fEfficiencyIPNonhfeD = fNonHFEEffbgc->GetFunction("fipfitnonhfe");
1499
1500 fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
1501 fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1502 fConversionEffbgc->Fit(fipfitnonhfe,"R");
1503 fEfficiencyIPConversionD = fConversionEffbgc->GetFunction("fipfitnonhfe");
1504 }
1505 //MHMH
1506 }
1507
1508 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1509 fipfit->SetLineColor(2);
1510 fEfficiencyBeautySigD->Fit(fipfit,"R");
1511 fEfficiencyBeautySigD->GetYaxis()->SetTitle("Efficiency");
1512 fEfficiencyIPBeautyD = fEfficiencyBeautySigD->GetFunction("fipfit");
1513
1514 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1515 fipfit->SetLineColor(6);
1516 fEfficiencyBeautySigesdD->Fit(fipfit,"R");
1517 fEfficiencyBeautySigesdD->GetYaxis()->SetTitle("Efficiency");
1518 fEfficiencyIPBeautyesdD = fEfficiencyBeautySigesdD->GetFunction("fipfit");
1519
1520 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1521 fipfit->SetLineColor(1);
1522 fEfficiencyCharmSigD->Fit(fipfit,"R");
1523 fEfficiencyCharmSigD->GetYaxis()->SetTitle("Efficiency");
1524 fEfficiencyIPCharmD = fEfficiencyCharmSigD->GetFunction("fipfit");
1525
1526 //if(0){
1527 if(fIPParameterizedEff){
1528 // fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1529 fipfitnonhfe->SetLineColor(3);
1530 if(fBeamType==1){
1531 fConversionEff->Fit(fipfitnonhfe,"R");
1532 fConversionEff->GetYaxis()->SetTitle("Efficiency");
1533 fEfficiencyIPConversionD = fConversionEff->GetFunction("fipfitnonhfe");
1534 }
1535 else{
1536 fConversionEffbgc->Fit(fipfitnonhfe,"R");
1537 fConversionEffbgc->GetYaxis()->SetTitle("Efficiency");
1538 fEfficiencyIPConversionD = fConversionEffbgc->GetFunction("fipfitnonhfe");
1539 }
1540 // fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1541 fipfitnonhfe->SetLineColor(4);
1542 if(fBeamType==1){
1543 fNonHFEEff->Fit(fipfitnonhfe,"R");
1544 fNonHFEEff->GetYaxis()->SetTitle("Efficiency");
1545 fEfficiencyIPNonhfeD = fNonHFEEff->GetFunction("fipfitnonhfe");
1546 }
1547 else{
1548 fNonHFEEffbgc->Fit(fipfitnonhfe,"R");
1549 fNonHFEEffbgc->GetYaxis()->SetTitle("Efficiency");
1550 fEfficiencyIPNonhfeD = fNonHFEEffbgc->GetFunction("fipfitnonhfe");
1551 }
1552 }
1553
1554 // draw
1555 fEfficiencyCharmSigD->SetMarkerStyle(21);
1556 fEfficiencyCharmSigD->SetMarkerColor(1);
1557 fEfficiencyCharmSigD->SetLineColor(1);
1558 fEfficiencyBeautySigD->SetMarkerStyle(21);
1559 fEfficiencyBeautySigD->SetMarkerColor(2);
1560 fEfficiencyBeautySigD->SetLineColor(2);
1561 fEfficiencyBeautySigesdD->SetStats(0);
1562 fEfficiencyBeautySigesdD->SetMarkerStyle(21);
1563 fEfficiencyBeautySigesdD->SetMarkerColor(6);
1564 fEfficiencyBeautySigesdD->SetLineColor(6);
1565 fCharmEff->SetMarkerStyle(24);
1566 fCharmEff->SetMarkerColor(1);
1567 fCharmEff->SetLineColor(1);
1568 fBeautyEff->SetMarkerStyle(24);
1569 fBeautyEff->SetMarkerColor(2);
1570 fBeautyEff->SetLineColor(2);
1571 fConversionEff->SetMarkerStyle(24);
1572 fConversionEff->SetMarkerColor(3);
1573 fConversionEff->SetLineColor(3);
1574 fNonHFEEff->SetMarkerStyle(24);
1575 fNonHFEEff->SetMarkerColor(4);
1576 fNonHFEEff->SetLineColor(4);
1577
1578 fEfficiencyCharmSigD->Draw();
1579 fEfficiencyCharmSigD->GetXaxis()->SetRangeUser(0.0,7.9);
1580 fEfficiencyCharmSigD->GetYaxis()->SetRangeUser(0.0,0.5);
1581
1582 fEfficiencyBeautySigD->Draw("same");
1583 fEfficiencyBeautySigesdD->Draw("same");
1584 if(fBeamType == 1){
1585 fNonHFEEff->Draw("same");
1586 fConversionEff->Draw("same");
1587 }
1588 //fCharmEff->Draw("same");
1589 //fBeautyEff->Draw("same");
1590
1591 if(fBeamType==0){
1592 fConversionEffbgc->SetMarkerStyle(25);
1593 fConversionEffbgc->SetMarkerColor(3);
1594 fConversionEffbgc->SetLineColor(3);
1595 fNonHFEEffbgc->SetMarkerStyle(25);
1596 fNonHFEEffbgc->SetMarkerColor(4);
1597 fNonHFEEffbgc->SetLineColor(4);
1598 fConversionEffbgc->Draw("same");
1599 fNonHFEEffbgc->Draw("same");
1600 }
1601
1602 if(fEfficiencyIPBeautyD)
1603 fEfficiencyIPBeautyD->Draw("same");
1604 if(fEfficiencyIPBeautyesdD)
1605 fEfficiencyIPBeautyesdD->Draw("same");
1606 if( fEfficiencyIPCharmD)
1607 fEfficiencyIPCharmD->Draw("same");
1608 if(fIPParameterizedEff){
1609 if(fEfficiencyIPConversionD)
1610 fEfficiencyIPConversionD->Draw("same");
1611 if(fEfficiencyIPNonhfeD)
1612 fEfficiencyIPNonhfeD->Draw("same");
1613 }
1614 TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39);
1615 legipeff->AddEntry(fEfficiencyBeautySigD,"IP Step Efficiency","");
1616 legipeff->AddEntry(fEfficiencyBeautySigD,"beauty e","p");
1617 legipeff->AddEntry(fEfficiencyBeautySigesdD,"beauty e(esd pt)","p");
1618 legipeff->AddEntry(fEfficiencyCharmSigD,"charm e","p");
1619 if(fBeamType == 0){
1620 legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p");
1621 legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p");
1622 }
1623 else{
1624 legipeff->AddEntry(fConversionEff,"conversion e","p");
1625 legipeff->AddEntry(fNonHFEEff,"Dalitz e","p");
1626 }
1627 legipeff->Draw("same");
1628 gPad->SetGrid();
1629 //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps");
1630}
1631
1632//____________________________________________________________________________
1633THnSparse* AliHFEBeautySpectrum::GetBeautyIPEff(Bool_t isMCpt){
1634 //
1635 // Return beauty electron IP cut efficiency
1636 //
1637
1638 const Int_t nPtbinning1 = kSignalPtBins;//number of pt bins, according to new binning
1639 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
1640
1641 Int_t nDim=1; //dimensions of the efficiency weighting grid
1642
1643 Int_t nBin[1] = {nPtbinning1};
1644
1645 THnSparseF *ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
1646
1647 ipcut->SetBinEdges(0, kPtRange);
1648
1649 Double_t pt[1];
1650 Double_t weight;
1651 Double_t weightErr;
1652 Double_t contents[2];
1653
1654 weight = 1.0;
1655 weightErr = 1.0;
1656
1657 Int_t looppt=nBin[0];
1658
1659 Int_t ibin[2];
1660
1661 for(int i=0; i<looppt; i++)
1662 {
1663 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1664 if(isMCpt){
1665 if(fEfficiencyIPBeautyD){
1666 weight=fEfficiencyIPBeautyD->Eval(pt[0]);
1667 weightErr = 0;
1668 }
1669 else{
1670 printf("Fit failed on beauty IP cut efficiency. Contents in histo used!\n");
1671 weight = fEfficiencyBeautySigD->GetBinContent(i+1);
1672 weightErr = fEfficiencyBeautySigD->GetBinError(i+1);
1673 }
1674 }
1675 else{
1676 if(fEfficiencyIPBeautyesdD){
1677 weight=fEfficiencyIPBeautyesdD->Eval(pt[0]);
1678 weightErr = 0;
1679 }
1680 else{
1681 printf("Fit failed on beauty IP cut efficiency. Contents in histo used!\n");
1682 weight = fEfficiencyBeautySigesdD->GetBinContent(i+1);
1683 weightErr = fEfficiencyBeautySigD->GetBinError(i+1);
1684 }
1685 }
1686
1687 contents[0]=pt[0];
1688 ibin[0]=i+1;
1689
1690 ipcut->Fill(contents,weight);
1691 ipcut->SetBinError(ibin,weightErr);
1692 }
1693
1694 Int_t nDimSparse = ipcut->GetNdimensions();
1695 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1696 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1697
1698 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1699 binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
1700 nBins *= binsvar[iVar];
1701 }
1702
1703 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1704 // loop that sets 0 error in each bin
1705 for (Long_t iBin=0; iBin<nBins; iBin++) {
1706 Long_t bintmp = iBin ;
1707 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1708 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1709 bintmp /= binsvar[iVar] ;
1710 }
1711 //ipcut->SetBinError(binfill,0.); // put 0 everywhere
1712 }
1713
1714 delete[] binsvar;
1715 delete[] binfill;
1716
1717 return ipcut;
1718}
1719
1720//____________________________________________________________________________
1721THnSparse* AliHFEBeautySpectrum::GetPIDxIPEff(Int_t source){
1722 //
1723 // Return PID x IP cut efficiency
1724 //
1725 const Int_t nPtbinning1 = kSignalPtBins;//number of pt bins, according to new binning
1726 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
1727 Int_t nDim=1; //dimensions of the efficiency weighting grid
1728
1729 Int_t nBin[1] = {nPtbinning1};
1730
1731 THnSparseF *pideff;
1732 pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
1733 pideff->SetBinEdges(0, kPtRange);
1734
1735 Double_t pt[1];
1736 Double_t weight;
1737 Double_t weightErr;
1738 Double_t contents[2];
1739
1740 weight = 1.0;
1741 weightErr = 1.0;
1742
1743 Int_t looppt=nBin[0];
1744 Int_t ibin[2];
1745
1746 Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
1747 for(int i=0; i<looppt; i++)
1748 {
1749 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1750
1751 Double_t tofpideff = 1.;
1752 Double_t tofpideffesd = 1.;
1753 if(fEfficiencyTOFPIDD)
1754 tofpideff = fEfficiencyTOFPIDD->Eval(pt[0]);
1755 else{
1756 printf("TOF PID fit failed on conversion. The result is wrong!\n");
1757 }
1758 if(fEfficiencyesdTOFPIDD)
1759 tofpideffesd = fEfficiencyesdTOFPIDD->Eval(pt[0]);
1760 else{
1761 printf("TOF esd PID fit failed on conversion. The result is wrong!\n");
1762 }
1763
1764 //tof pid eff x tpc pid eff x ip cut eff
1765 if(fIPParameterizedEff){
1766 if(source==0) {
1767 if(fEfficiencyIPCharmD){
1768 weight = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD->Eval(pt[0]);
1769 weightErr = 0;
1770 }
1771 else{
1772 printf("Fit failed on charm IP cut efficiency\n");
1773 weight = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinContent(i+1);
1774 weightErr = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinError(i+1);
1775 }
1776 }
1777 else if(source==2) {
1778 if(fEfficiencyIPConversionD){
1779 weight = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD->Eval(pt[0]);
1780 weightErr = 0;
1781 }
1782 else{
1783 printf("Fit failed on conversion IP cut efficiency\n");
1784 weight = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1);
1785 weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinError(i+1);
1786 }
1787 }
1788 else if(source==3) {
1789 if(fEfficiencyIPNonhfeD){
1790 weight = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD->Eval(pt[0]);
1791 weightErr = 0;
1792 }
1793 else{
1794 printf("Fit failed on dalitz IP cut efficiency\n");
1795 weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1);
1796 weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinError(i+1);
1797 }
1798 }
1799 }
1800 else{
1801 if(source==0){
1802 if(fEfficiencyIPCharmD){
1803 weight = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD->Eval(pt[0]);
1804 weightErr = 0;
1805 }
1806 else{
1807 printf("Fit failed on charm IP cut efficiency\n");
1808 weight = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinContent(i+1);
1809 weightErr = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinError(i+1);
1810 }
1811 }
1812 else if(source==2){
1813 if(fBeamType==0){
1814 weight = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
1815 weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
1816 }
1817 else{
1818 weight = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1); // conversion
1819 weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinError(i+1);
1820 }
1821 }
1822 else if(source==3){
1823 if(fBeamType==0){
1824 weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
1825 weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
1826 }
1827 else{
1828 weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1); // Dalitz
1829 weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinError(i+1);
1830 }
1831 }
1832 }
1833
1834 contents[0]=pt[0];
1835 ibin[0]=i+1;
1836
1837 pideff->Fill(contents,weight);
1838 pideff->SetBinError(ibin,weightErr);
1839 }
1840
1841 Int_t nDimSparse = pideff->GetNdimensions();
1842 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1843 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1844
1845 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1846 binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
1847 nBins *= binsvar[iVar];
1848 }
1849
1850 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1851 // loop that sets 0 error in each bin
1852 for (Long_t iBin=0; iBin<nBins; iBin++) {
1853 Long_t bintmp = iBin ;
1854 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1855 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1856 bintmp /= binsvar[iVar] ;
1857 }
1858 }
1859
1860 delete[] binsvar;
1861 delete[] binfill;
1862
1863
1864 return pideff;
1865}
1866
1867//__________________________________________________________________________
1868AliCFDataGrid *AliHFEBeautySpectrum::GetRawBspectra2ndMethod(){
1869 //
1870 // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method
1871 //
1872 Int_t nDim = 1;
1873
1874 const Int_t nPtbinning1 = 18;//number of pt bins, according to new binning
1875 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};
1876
1877 Int_t nBin[1] = {nPtbinning1};
1878
1879 AliCFDataGrid *rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
1880
1881 THnSparseF *brawspectra;
1882 brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin);
1883
1884 brawspectra->SetBinEdges(0, kPtRange);
1885
1886 Double_t pt[1];
1887 Double_t yields= 0.;
1888 Double_t valuesb[2];
1889
1890 for(int i=0; i<fBSpectrum2ndMethod->GetNbinsX(); i++){
1891 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1892
1893 yields = fBSpectrum2ndMethod->GetBinContent(i+1);
1894
1895 valuesb[0]=pt[0];
1896 brawspectra->Fill(valuesb,yields);
1897 }
1898
1899
1900
1901 Int_t nDimSparse = brawspectra->GetNdimensions();
1902 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1903 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1904
1905 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1906 binsvar[iVar] = brawspectra->GetAxis(iVar)->GetNbins();
1907 nBins *= binsvar[iVar];
1908 }
1909
1910 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1911 // loop that sets 0 error in each bin
1912 for (Long_t iBin=0; iBin<nBins; iBin++) {
1913 Long_t bintmp = iBin ;
1914 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1915 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1916 bintmp /= binsvar[iVar] ;
1917 }
1918 brawspectra->SetBinError(binfill,0.); // put 0 everywhere
1919 }
1920
1921
1922 rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency
1923 TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(0);
1924
1925 new TCanvas;
1926 fBSpectrum2ndMethod->SetMarkerStyle(24);
1927 fBSpectrum2ndMethod->Draw("p");
1928 hRawBeautySpectra->SetMarkerStyle(25);
1929 hRawBeautySpectra->Draw("samep");
1930
1931 delete[] binfill;
1932 delete[] binsvar;
1933
1934 return rawBeautyContainer;
1935}
1936//__________________________________________________________________________
1937void AliHFEBeautySpectrum::CalculateNonHFEsyst(){
1938 //
1939 // Calculate non HFE sys
1940 //
1941 //
1942
1943 if(!fNonHFEsyst)
1944 return;
1945
1946 Double_t evtnorm[1] = {0.0};
1947 if(fNMCbgEvents>0) evtnorm[0]= double(fNEvents)/double(fNMCbgEvents);
1948
1949 AliCFDataGrid *convSourceGrid[kElecBgSources-3][kBgLevels];
1950 AliCFDataGrid *nonHFESourceGrid[kElecBgSources-3][kBgLevels];
1951
1952 AliCFDataGrid *bgLevelGrid[2][kBgLevels];//for pi0 and eta based errors
1953 AliCFDataGrid *bgNonHFEGrid[kBgLevels];
1954 AliCFDataGrid *bgConvGrid[kBgLevels];
1955
1956 Int_t stepbackground = 3;
1957 Int_t* bins=new Int_t[1];
1958 const Char_t *bgBase[2] = {"pi0","eta"};
1959
1960 bins[0]=kSignalPtBins;//fConversionEff[centrality]->GetNbinsX();
1961
1962 AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
1963 AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
1964
1965 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1966 for(Int_t iSource = 0; iSource < kElecBgSources-3; iSource++){
1967 convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d",iSource,iLevel),Form("convGrid_%d_%d",iSource,iLevel),*fConvSourceContainer[iSource][iLevel],stepbackground);
1968 weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
1969 convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
1970
1971 nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d",iSource,iLevel),Form("nonHFEGrid_%d_%d",iSource,iLevel),*fNonHFESourceContainer[iSource][iLevel],stepbackground);
1972 weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1973 nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
1974 }
1975
1976 bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
1977 for(Int_t iSource = 2; iSource < kElecBgSources-3; iSource++){
1978 bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
1979 }
1980 if(!fEtaSyst)
1981 bgConvGrid[iLevel]->Add(convSourceGrid[1][iLevel]);
1982
1983 bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone();
1984 for(Int_t iSource = 2; iSource < kElecBgSources-3; iSource++){//add other sources to pi0, to get overall background from all meson decays, exception: eta (independent error calculation)
1985 bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
1986 }
1987 if(!fEtaSyst)
1988 bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[1][iLevel]);
1989
1990 bgLevelGrid[0][iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
1991 bgLevelGrid[0][iLevel]->Add(bgNonHFEGrid[iLevel]);
1992 if(fEtaSyst){
1993 bgLevelGrid[1][iLevel] = (AliCFDataGrid*)nonHFESourceGrid[1][iLevel]->Clone();//background for eta source
1994 bgLevelGrid[1][iLevel]->Add(convSourceGrid[1][iLevel]);
1995 }
1996 }
1997
1998 //Now subtract the mean from upper, and lower from mean container to get the error based on the pion yield uncertainty (-> this error sums linearly, since its contribution to all meson yields is correlated; exception: eta errors in pp 7 TeV sum with others the gaussian way, as they are independent from pi0)
1999 AliCFDataGrid *bgErrorGrid[2][2];//for pions/eta error base, for lower/upper
2000 TH1D* hBaseErrors[2][2];//pi0/eta and lower/upper
2001 for(Int_t iErr = 0; iErr < 2; iErr++){//errors for pi0 and eta base
2002 bgErrorGrid[iErr][0] = (AliCFDataGrid*)bgLevelGrid[iErr][1]->Clone();
2003 bgErrorGrid[iErr][0]->Add(bgLevelGrid[iErr][0],-1.);
2004 bgErrorGrid[iErr][1] = (AliCFDataGrid*)bgLevelGrid[iErr][2]->Clone();
2005 bgErrorGrid[iErr][1]->Add(bgLevelGrid[iErr][0],-1.);
2006
2007 //plot absolute differences between limit yields (upper/lower limit, based on pi0 and eta errors) and best estimate
2008
2009 hBaseErrors[iErr][0] = (TH1D*)bgErrorGrid[iErr][0]->Project(0);
2010 hBaseErrors[iErr][0]->Scale(-1.);
2011 hBaseErrors[iErr][0]->SetTitle(Form("Absolute %s-based systematic errors from non-HF meson decays and conversions",bgBase[iErr]));
2012 hBaseErrors[iErr][1] = (TH1D*)bgErrorGrid[iErr][1]->Project(0);
2013 if(!fEtaSyst)break;
2014 }
2015
2016 //Calculate the scaling errors for electrons from all mesons except for pions (and in pp 7 TeV case eta): square sum of (0.3 * best yield estimate), where 0.3 is the error generally assumed for m_t scaling
2017 TH1D *hSpeciesErrors[kElecBgSources-4];
2018 for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
2019 if(fEtaSyst && (iSource == 1))continue;
2020 hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0);
2021 TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0);
2022 hSpeciesErrors[iSource-1]->Add(hNonHFEtemp);
2023 hSpeciesErrors[iSource-1]->Scale(0.3);
2024 }
2025
2026 //Int_t firstBgSource = 0;//if eta systematics are not from scaling
2027 //if(fEtaSyst){firstBgSource = 1;}//source 0 histograms are not filled if eta errors are independently determined!
2028 TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[1]->Clone();
2029 TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[1]->Clone();
2030 TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[1]->Clone();
2031
2032 TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
2033 TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
2034
2035 for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
2036 Double_t pi0basedErrLow,pi0basedErrUp,etaErrLow,etaErrUp;
2037 pi0basedErrLow = hBaseErrors[0][0]->GetBinContent(iBin);
2038 pi0basedErrUp = hBaseErrors[0][1]->GetBinContent(iBin);
2039 if(fEtaSyst){
2040 etaErrLow = hBaseErrors[1][0]->GetBinContent(iBin);
2041 etaErrUp = hBaseErrors[1][1]->GetBinContent(iBin);
2042 }
2043 else{ etaErrLow = etaErrUp = 0.;}
2044
2045 Double_t sqrsumErrs= 0;
2046 for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
2047 if(fEtaSyst && (iSource == 1))continue;
2048 Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
2049 sqrsumErrs+=(scalingErr*scalingErr);
2050 }
2051 for(Int_t iErr = 0; iErr < 2; iErr++){
2052 for(Int_t iLevel = 0; iLevel < 2; iLevel++){
2053 hBaseErrors[iErr][iLevel]->SetBinContent(iBin,hBaseErrors[iErr][iLevel]->GetBinContent(iBin)/hBaseErrors[iErr][iLevel]->GetBinWidth(iBin));
2054 }
2055 if(!fEtaSyst)break;
2056 }
2057 hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+(etaErrUp*etaErrUp)+sqrsumErrs));
2058 hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+(etaErrLow*etaErrLow)+sqrsumErrs));
2059 hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin));
2060
2061 hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin));
2062 hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin));
2063 }
2064
2065 TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
2066 cPiErrors->cd();
2067 cPiErrors->SetLogx();
2068 cPiErrors->SetLogy();
2069 hBaseErrors[0][0]->Draw();
2070 //hBaseErrors[0][1]->SetMarkerColor(kBlack);
2071 //hBaseErrors[0][1]->SetLineColor(kBlack);
2072 //hBaseErrors[0][1]->Draw("SAME");
2073 if(fEtaSyst){
2074 hBaseErrors[1][0]->Draw("SAME");
2075 hBaseErrors[1][0]->SetMarkerColor(kBlack);
2076 hBaseErrors[1][0]->SetLineColor(kBlack);
2077 //hBaseErrors[1][1]->SetMarkerColor(13);
2078 //hBaseErrors[1][1]->SetLineColor(13);
2079 //hBaseErrors[1][1]->Draw("SAME");
2080 }
2081 //hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
2082 //hOverallBinScaledErrsUp->SetLineColor(kBlue);
2083 //hOverallBinScaledErrsUp->Draw("SAME");
2084 hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
2085 hOverallBinScaledErrsLow->SetLineColor(kGreen);
2086 hOverallBinScaledErrsLow->Draw("SAME");
2087 hScalingErrors->SetLineColor(kBlue);
2088 hScalingErrors->Draw("SAME");
2089
2090 TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
2091 lPiErr->AddEntry(hBaseErrors[0][0],"Lower error from pion error");
2092 //lPiErr->AddEntry(hBaseErrors[0][1],"Upper error from pion error");
2093 if(fEtaSyst){
2094 lPiErr->AddEntry(hBaseErrors[1][0],"Lower error from eta error");
2095 //lPiErr->AddEntry(hBaseErrors[1][1],"Upper error from eta error");
2096 }
2097 lPiErr->AddEntry(hScalingErrors, "scaling error");
2098 lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics");
2099 //lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
2100 lPiErr->Draw("SAME");
2101
2102 //Normalize errors
2103 TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone();
2104 TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone();
2105 hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!)
2106 hLowSystScaled->Scale(evtnorm[0]);
2107 TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone();
2108 TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone();
2109 //histograms to be normalized to TGraphErrors
afb48e1d 2110 AliHFEtools::NormaliseBinWidth(hNormAllSystErrUp);
2111 AliHFEtools::NormaliseBinWidth(hNormAllSystErrLow);
4bd4fc13 2112
2113 TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs");
2114 cNormOvErrs->cd();
2115 cNormOvErrs->SetLogx();
2116 cNormOvErrs->SetLogy();
2117
2118 TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp);
2119 TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow);
2120 gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors");
2121 gOverallSystErrUp->SetMarkerColor(kBlack);
2122 gOverallSystErrUp->SetLineColor(kBlack);
2123 gOverallSystErrLow->SetMarkerColor(kRed);
2124 gOverallSystErrLow->SetLineColor(kRed);
2125 gOverallSystErrUp->Draw("AP");
2126 gOverallSystErrLow->Draw("PSAME");
2127 TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89);
2128 lAllSys->AddEntry(gOverallSystErrLow,"lower","p");
2129 lAllSys->AddEntry(gOverallSystErrUp,"upper","p");
2130 lAllSys->Draw("same");
2131
2132
2133 AliCFDataGrid *bgYieldGrid;
2134 if(fEtaSyst){
2135 bgLevelGrid[0][0]->Add(bgLevelGrid[1][0]);//Addition of the eta background best estimate to the rest. Needed to be separated for error treatment - now overall background necessary! If no separate eta systematics exist, the corresponding grid has already been added before.
2136 }
2137 bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0][0]->Clone();
2138
2139 TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0);
2140 TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone();
2141 hRelErrUp->Divide(hBgYield);
2142 TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone();
2143 hRelErrLow->Divide(hBgYield);
2144
2145 TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs");
2146 cRelErrs->cd();
2147 cRelErrs->SetLogx();
2148 hRelErrUp->SetTitle("Relative error of non-HFE background yield");
2149 hRelErrUp->Draw();
2150 hRelErrLow->SetLineColor(kBlack);
2151 hRelErrLow->Draw("SAME");
2152
2153 TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95);
2154 lRel->AddEntry(hRelErrUp, "upper");
2155 lRel->AddEntry(hRelErrLow, "lower");
2156 lRel->Draw("SAME");
2157
afb48e1d 2158 //AliHFEtools::NormaliseBinWidth(hBgYield);
4bd4fc13 2159 //hBgYield->Scale(evtnorm[0]);
2160
2161
2162 //write histograms/TGraphs to file
2163 TFile *output = new TFile("systHists.root","recreate");
2164
2165 hBgYield->SetName("hBgYield");
2166 hBgYield->Write();
2167 hRelErrUp->SetName("hRelErrUp");
2168 hRelErrUp->Write();
2169 hRelErrLow->SetName("hRelErrLow");
2170 hRelErrLow->Write();
2171 hUpSystScaled->SetName("hOverallSystErrUp");
2172 hUpSystScaled->Write();
2173 hLowSystScaled->SetName("hOverallSystErrLow");
2174 hLowSystScaled->Write();
2175 gOverallSystErrUp->SetName("gOverallSystErrUp");
2176 gOverallSystErrUp->Write();
2177 gOverallSystErrLow->SetName("gOverallSystErrLow");
2178 gOverallSystErrLow->Write();
2179
2180 output->Close();
2181 delete output;
2182}
2183//____________________________________________________________
2184void AliHFEBeautySpectrum::UnfoldBG(AliCFDataGrid* const bgsubpectrum){
2185 //
2186 // Unfold backgrounds to check its sanity
2187 //
2188
2189 AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
2190 //AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
2191 if(!mcContainer){
2192 AliError("MC Container not available");
2193 return;
2194 }
2195
2196 if(!fCorrelation){
2197 AliError("No Correlation map available");
2198 return;
2199 }
2200
2201 // Data
2202 AliCFDataGrid *dataGrid = 0x0;
2203 dataGrid = bgsubpectrum;
2204
2205 // Guessed
2206 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
2207 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
2208
2209 // Efficiency
2210 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
2211 efficiencyD->CalculateEfficiency(fStepMC+2,fStepTrue);
2212
2213 // Unfold background spectra
2214 Int_t nDim=1;
2215 AliCFUnfolding unfolding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
2216 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
2217 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
2218 if(fSetSmoothing) unfolding.UseSmoothing();
2219 unfolding.Unfold();
2220
2221 // Results
2222 THnSparse* result = unfolding.GetUnfolded();
2223 TCanvas *ctest = new TCanvas("yvonnetest","yvonnetest",1000,600);
2224 if(fBeamType==1)
2225 {
2226 ctest->Divide(2);
2227 ctest->cd(1);
2228 TH1D* htest1=(TH1D*)result->Projection(0);
2229 htest1->Draw();
2230 ctest->cd(2);
2231 TH1D* htest2=(TH1D*)result->Projection(1);
2232 htest2->Draw();
2233 }
2234
2235
2236 TGraphErrors* unfoldedbgspectrumD = Normalize(result);
2237 if(!unfoldedbgspectrumD) {
2238 AliError("Unfolded background spectrum doesn't exist");
2239 }
2240 else{
2241 TFile *file = TFile::Open("unfoldedbgspectrum.root","recreate");
2242 if(fBeamType==0)unfoldedbgspectrumD->Write("unfoldedbgspectrum");
2243
2244 file->Close();
2245 }
2246}