]> 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
Updates for TRD HFE analysis
[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);
592 CorrectFromTheWidth(incElec);
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);
611 CorrectFromTheWidth(hadron);
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);
621 CorrectFromTheWidth(charm);
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);
637 CorrectFromTheWidth(nonHFE[iPlot]);
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);
690 CorrectFromTheWidth(subtracted);
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);
699 CorrectFromTheWidth(measuredTH1background);
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++){
741 backgroundContainer->Add(fNonHFESourceContainer[iSource][0]);
742 }
743 }
744 }
745 else{//careful: these containers include the non-standard electron sources (K0s, etc.). If necessary, use SetRange in species axis to fix it?
746 if(source == 0){
747 // Background Estimate
748 if(decay == 0)
749 backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
750 else if(decay == 1)
751 backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
752 }
753 }
754 if(!backgroundContainer){
755 AliError("MC background container not found");
756 return NULL;
757 }
758 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("backgroundGrid","backgroundGrid",*backgroundContainer,stepbackground);
759
760 Double_t rangelow = 1.;
761 Double_t rangeup = 6.;
762 if(decay == 1) rangelow = 0.9;
763
764
765 const Char_t *dmode[2]={"Conversions","Dalitz decays"};
766 TF1 *fitHagedorn = new TF1("fitHagedorn", "[0]/TMath::Power(([1]+x/[2]), [3])", rangelow, rangeup);
767 fNonHFEbg = 0x0;
768 TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
769 TAxis *axis = h1->GetXaxis();
770 CorrectFromTheWidth(h1);
771 if(source == 0){
772 fitHagedorn->SetParameter(0, 0.15);
773 fitHagedorn->SetParameter(1, 0.09);
774 fitHagedorn->SetParameter(2, 8.4);
775 fitHagedorn->SetParameter(3, 6.3);
776 // TCanvas *ch1conversion = new TCanvas("ch1conversion","ch1conversion",500,400);
777 // ch1conversion->cd();
778 //fHagedorn->SetLineColor(2);
779 h1->Fit(fitHagedorn,"R");
780 // h1->Draw();
781 if(!(fNonHFEbg = h1->GetFunction("fitHagedorn")))printf("electron background fit failed for %s\n",dmode[decay]);
782 }
783
784 Int_t *nBinpp=new Int_t[1];
785 Int_t *binspp=new Int_t[1];
786 binspp[0]=kSignalPtBins;// number of pt bins
787
788 Int_t looppt=binspp[0];
789
790 for(Long_t iBin=1; iBin<= looppt;iBin++){
791 Double_t iipt= h1->GetBinCenter(iBin);
792 Double_t iiwidth= axis->GetBinWidth(iBin);
793 nBinpp[0]=iBin;
794 Double_t fitFactor = backgroundGrid->GetElement(nBinpp);//if no fit available, just take bin-by-bin information
795 if(fNonHFEbg)fitFactor = fNonHFEbg->Eval(iipt)*iiwidth;
796 backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
797 backgroundGrid->SetElement(nBinpp,fitFactor*evtnorm[0]);
798 }
4bd4fc13 799 //end of workaround for statistical errors
800 if(source == 0){
801 AliCFDataGrid *weightedBGContainer = 0x0;
802 weightedBGContainer = new AliCFDataGrid("weightedBGContainer","weightedBGContainer",1,binspp);
803 if(decay == 0)
804 weightedBGContainer->SetGrid(GetPIDxIPEff(2));
805 else if(decay == 1)
806 weightedBGContainer->SetGrid(GetPIDxIPEff(3));
807 if(stepbackground == 3){
808 backgroundGrid->Multiply(weightedBGContainer,1.0);
809 }
32c709ed 810 }
811 delete[] nBinpp;
812 delete[] binspp;
4bd4fc13 813 return backgroundGrid;
814}
815//____________________________________________________________
816AliCFDataGrid* AliHFEBeautySpectrum::GetCharmBackground(){
817 //
818 // calculate charm background; when using centrality binning 2: centBin = 0 for 0-20%, centBin = 5 for 40-80%
819 //
820
821 Int_t nDim = 1;
822
823 Double_t evtnorm=0;
824 if(fNMCEvents) evtnorm= double(fNEvents)/double(fNMCEvents); //mjmjmj
825 //if(fNMCbgEvents) evtnorm= double(fNEvents)/double(fNMCbgEvents);
826 printf("events for data: %d",fNEvents);
827 printf("events for MC: %d",fNMCEvents);
828 printf("normalization factor for charm: %f",evtnorm);
829
830 AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
831 if(!mcContainer){
832 AliError("MC Container not available");
833 return NULL;
834 }
835
836 if(!fCorrelation){
837 AliError("No Correlation map available");
838 return NULL;
839 }
840
841 AliCFDataGrid *charmBackgroundGrid= 0x0;
842 charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
843 TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
844 Int_t* bins=new Int_t[2];
845 bins[0]=charmbgaftertofpid->GetNbinsX();
846
847 AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins);
848 ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
849 TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(0);
850
851 charmBackgroundGrid->Multiply(ipWeightedCharmContainer,1.);
852 Int_t *nBinpp=new Int_t[1];
853 Int_t* binspp=new Int_t[1];
854 binspp[0]=charmbgaftertofpid->GetNbinsX(); // number of pt bins
855
856 Int_t looppt=binspp[0];
857
858 for(Long_t iBin=1; iBin<= looppt;iBin++){
859 nBinpp[0]=iBin;
860 charmBackgroundGrid->SetElementError(nBinpp, charmBackgroundGrid->GetElementError(nBinpp)*evtnorm);
861 charmBackgroundGrid->SetElement(nBinpp,charmBackgroundGrid->GetElement(nBinpp)*evtnorm);
862 }
863
864 TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(0);
865
866 AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
867 weightedCharmContainer->SetGrid(GetCharmWeights(fTestCentralityLow)); // get charm weighting factors
868 TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(0);
869 charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
870 TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(0);
871
872 // Efficiency (set efficiency to 1 for only folding)
873 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
874 efficiencyD->CalculateEfficiency(0,0);
875
876 // Folding
877 if(fBeamType==0)nDim = 1;
878 AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
879 folding.SetMaxNumberOfIterations(1);
880 folding.Unfold();
881
882 // Results
883 THnSparse* result1= folding.GetEstMeasured(); // folded spectra
884 THnSparse* result=(THnSparse*)result1->Clone();
885 charmBackgroundGrid->SetGrid(result);
886 TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(0);
887
888 //Charm background evaluation plots
889
890 TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
891 cCharmBgEval->Divide(3,1);
892
893 cCharmBgEval->cd(1);
894
895 if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm);
896
897 CorrectFromTheWidth(charmbgaftertofpid);
898 charmbgaftertofpid->SetMarkerStyle(25);
899 charmbgaftertofpid->Draw("p");
900 charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
901 charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
902 gPad->SetLogy();
903
904 CorrectFromTheWidth(charmbgafteripcut);
905 charmbgafteripcut->SetMarkerStyle(24);
906 charmbgafteripcut->Draw("samep");
907
908 CorrectFromTheWidth(charmbgafterweight);
909 charmbgafterweight->SetMarkerStyle(24);
910 charmbgafterweight->SetMarkerColor(4);
911 charmbgafterweight->Draw("samep");
912
913 CorrectFromTheWidth(charmbgafterfolding);
914 charmbgafterfolding->SetMarkerStyle(24);
915 charmbgafterfolding->SetMarkerColor(2);
916 charmbgafterfolding->Draw("samep");
917
918 cCharmBgEval->cd(2);
919 parametrizedcharmpidipeff->SetMarkerStyle(24);
920 parametrizedcharmpidipeff->Draw("p");
921 parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
922
923 cCharmBgEval->cd(3);
924 charmweightingfc->SetMarkerStyle(24);
925 charmweightingfc->Draw("p");
926 charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
927 charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
928
929 cCharmBgEval->cd(1);
930 TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
931 legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
932 legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
933 legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
934 legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
935 legcharmbg->Draw("same");
936
937 cCharmBgEval->cd(2);
938 TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
939 legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
940 legcharmbg2->Draw("same");
941
942 CorrectStatErr(charmBackgroundGrid);
943 if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
944
945 delete[] bins;
946 delete[] nBinpp;
947 delete[] binspp;
948
949 if(fUnfoldBG) UnfoldBG(charmBackgroundGrid);
950
951 return charmBackgroundGrid;
952}
953//____________________________________________________________
954AliCFDataGrid *AliHFEBeautySpectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
955 //
956 // Apply TPC pid efficiency correction from parametrisation
957 //
958
959 // Data in the right format
960 AliCFDataGrid *dataGrid = 0x0;
961 if(bgsubpectrum) {
962 dataGrid = bgsubpectrum;
963 }
964 else {
965
966 AliCFContainer *dataContainer = GetContainer(kDataContainer);
967 if(!dataContainer){
968 AliError("Data Container not available");
969 return NULL;
970 }
971 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
972 }
973 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
974 result->SetName("ParametrizedEfficiencyBefore");
975 THnSparse *h = result->GetGrid();
976 Int_t nbdimensions = h->GetNdimensions();
977 //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
978 AliCFContainer *dataContainer = GetContainer(kDataContainer);
979 if(!dataContainer){
980 AliError("Data Container not available");
981 return NULL;
982 }
983 AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
984 dataContainerbis->Add(dataContainerbis,-1.0);
985
986 Int_t* coord = new Int_t[nbdimensions];
987 memset(coord, 0, sizeof(Int_t) * nbdimensions);
988 Double_t* points = new Double_t[nbdimensions];
989
990 ULong64_t nEntries = h->GetNbins();
991 for (ULong64_t i = 0; i < nEntries; ++i) {
992 Double_t value = h->GetBinContent(i, coord);
993 //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
994 //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
995
996 // Get the bin co-ordinates given an coord
997 for (Int_t j = 0; j < nbdimensions; ++j)
998 points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
999
1000 if (!fEfficiencyFunction->IsInside(points))
1001 continue;
1002 TF1::RejectPoint(kFALSE);
1003
1004 // Evaulate function at points
1005 Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1006 //printf("Value efficiency is %f\n",valueEfficiency);
1007
1008 if(valueEfficiency > 0.0) {
1009 h->SetBinContent(coord,value/valueEfficiency);
1010 dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1011 }
1012 Double_t error = h->GetBinError(i);
1013 h->SetBinError(coord,error/valueEfficiency);
1014 dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1015 }
1016
1017 delete[] coord;
1018 delete[] points;
1019
1020 AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1021
1022// QA
1023 TH1D *afterE = (TH1D *) resultt->Project(0);
1024 CorrectFromTheWidth(afterE);
1025 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1026 CorrectFromTheWidth(beforeE);
1027 fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterPE);
1028 fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforePE);
1029 fQA->AddResultAt(fEfficiencyFunction,AliHFEBeautySpectrumQA::kPEfficiency);
1030 fQA->DrawCorrectWithEfficiency(AliHFEBeautySpectrumQA::kParametrized);
1031
1032 return resultt;
1033}
1034//____________________________________________________________
1035THnSparse *AliHFEBeautySpectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
1036
1037 //
1038 // Unfold and eventually correct for efficiency the bgsubspectrum
1039 //
1040
1041 AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
1042 if(!mcContainer){
1043 AliError("MC Container not available");
1044 return NULL;
1045 }
1046
1047 if(!fCorrelation){
1048 AliError("No Correlation map available");
1049 return NULL;
1050 }
1051
1052 // Data
1053 AliCFDataGrid *dataGrid = 0x0;
1054 if(bgsubpectrum) {
1055 dataGrid = bgsubpectrum;
1056 }
1057 else {
1058
1059 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1060 if(!dataContainer){
1061 AliError("Data Container not available");
1062 return NULL;
1063 }
1064
1065 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1066 }
1067
1068 // Guessed
1069 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
1070 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
1071
1072 // Efficiency
1073 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1074 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1075
1076 if(!fBeauty2ndMethod)
1077 {
1078 // Consider parameterized IP cut efficiency
1079 Int_t* bins=new Int_t[1];
1080 bins[0]=kSignalPtBins;
1081
1082 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1083 beffContainer->SetGrid(GetBeautyIPEff(kTRUE));
1084 efficiencyD->Multiply(beffContainer,1);
1085 }
1086
1087
1088 // Unfold
1089
1090 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...
1091 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
1092 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
1093 if(fSetSmoothing) unfolding.UseSmoothing();
1094 unfolding.Unfold();
1095
1096 // Results
1097 THnSparse* result = unfolding.GetUnfolded();
1098 THnSparse* residual = unfolding.GetEstMeasured();
1099
1100 // QA
1101 TH1D *residualh = (TH1D *) residual->Projection(0);
1102 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1103 TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
1104 TH1D *afterE = (TH1D *) result->Projection(0);
1105
1106 CorrectFromTheWidth(residualh);
1107 CorrectFromTheWidth(beforeE);
1108 CorrectFromTheWidth(afterE);
1109 fQA->AddResultAt(residualh,AliHFEBeautySpectrumQA::kResidualU);
1110 fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterU);
1111 fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforeU);
1112 fQA->AddResultAt(efficiencyDproj,AliHFEBeautySpectrumQA::kUEfficiency);
1113 fQA->DrawUnfolding();
1114
1115 return (THnSparse *) result->Clone();
1116}
1117
1118//____________________________________________________________
1119AliCFDataGrid *AliHFEBeautySpectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
1120
1121 //
1122 // Apply unfolding and efficiency correction together to bgsubspectrum
1123 //
1124
1125 AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
1126 if(!mcContainer){
1127 AliError("MC Container not available");
1128 return NULL;
1129 }
1130
1131 // Efficiency
1132 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1133 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1134
1135
1136 if(!fBeauty2ndMethod)
1137 {
1138 // Consider parameterized IP cut efficiency
1139 Int_t* bins=new Int_t[1];
1140 bins[0]=kSignalPtBins;
1141
1142 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1143 beffContainer->SetGrid(GetBeautyIPEff(kFALSE));
1144 efficiencyD->Multiply(beffContainer,1);
1145 }
1146
1147 // Data in the right format
1148 AliCFDataGrid *dataGrid = 0x0;
1149 if(bgsubpectrum) {
1150 dataGrid = bgsubpectrum;
1151 }
1152 else {
1153
1154 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1155 if(!dataContainer){
1156 AliError("Data Container not available");
1157 return NULL;
1158 }
1159
1160 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1161 }
1162
1163 // Correct
1164 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1165 result->ApplyEffCorrection(*efficiencyD);
1166
1167 // QA
1168 TH1D *afterE = (TH1D *) result->Project(0);
1169 CorrectFromTheWidth(afterE);
1170 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1171 CorrectFromTheWidth(beforeE);
1172 TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
1173 fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterMCE);
1174 fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforeMCE);
1175 fQA->AddResultAt(efficiencyDproj,AliHFEBeautySpectrumQA::kMCEfficiency);
1176 fQA->DrawCorrectWithEfficiency(AliHFEBeautySpectrumQA::kMC);
1177
1178 return result;
1179
1180}
1181//____________________________________________________________
1182void AliHFEBeautySpectrum::AddTemporaryObject(TObject *o){
1183 //
1184 // Emulate garbage collection on class level
1185 //
1186 if(!fTemporaryObjects) {
1187 fTemporaryObjects= new TList;
1188 fTemporaryObjects->SetOwner();
1189 }
1190 fTemporaryObjects->Add(o);
1191}
1192
1193//____________________________________________________________
1194void AliHFEBeautySpectrum::ClearObject(TObject *o){
1195 //
1196 // Do a safe deletion
1197 //
1198 if(fTemporaryObjects){
1199 if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
1200 delete o;
1201 } else{
1202 // just delete
1203 delete o;
1204 }
1205}
1206//_________________________________________________________________________
1207THnSparse* AliHFEBeautySpectrum::GetCharmWeights(Int_t centBin){
1208
1209 //
1210 // Measured D->e based weighting factors
1211 //
1212
1213 const Int_t nDimpp=1;
1214 Int_t nBinpp[nDimpp] = {kSignalPtBins};
1215 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.};
1216 Int_t looppt=nBinpp[0];
1217
1218 fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
1219 fWeightCharm->SetBinEdges(0, ptbinning1);
1220
1221 // //if(fBeamType == 0){// Weighting factor for pp
1222 //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};
1223//TPC+TOF standard cut 4800
1224 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};
1225 //}
1226 //else{
1227 //if(centBin == 0){
1228 // Weighting factor for PbPb (0-20%)
1229 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};
1230 //}
1231 //else{
1232 // Weighting factor for PbPb (40-80%)
1233 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};
1234 // }
1235 //}
1236 Double_t weight[kSignalPtBins];
1237 for(Int_t ipt = 0; ipt < kSignalPtBins; ipt++){
1238 if(fBeamType == 0)weight[ipt] = weightpp[ipt];
1239 else if(centBin == 0)weight[ipt] = weightPbPb1[ipt];
1240 else weight[ipt] = weightPbPb2[ipt];
1241 }
1242
1243 //points
1244 Double_t pt[1];
1245 Double_t contents[2];
1246
1247 for(int i=0; i<looppt; i++){
1248 pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
1249 contents[0]=pt[0];
1250 fWeightCharm->Fill(contents,weight[i]);
1251 }
1252
1253
1254 Int_t nDimSparse = fWeightCharm->GetNdimensions();
1255 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1256 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1257
1258 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1259 binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
1260 nBins *= binsvar[iVar];
1261 }
1262
1263 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1264 // loop that sets 0 error in each bin
1265 for (Long_t iBin=0; iBin<nBins; iBin++) {
1266 Long_t bintmp = iBin ;
1267 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1268 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1269 bintmp /= binsvar[iVar] ;
1270 }
1271 fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
1272 }
1273 delete[] binsvar;
1274 delete[] binfill;
1275
1276 return fWeightCharm;
1277}
1278//____________________________________________________________________________
1279void AliHFEBeautySpectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
1280
1281 // TOF PID efficiencies
1282
1283 TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1284 TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,6.);
1285 TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.0);
1286
1287 if(fBeamType == 1){//not in use yet - adapt as soon as possible!
1288 fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1289 fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1290 fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
1291 }
1292
1293 TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600);
1294 cefficiencyParamtof->cd();
1295
1296 AliCFContainer *mccontainermcD = 0x0;
1297 AliCFContainer *mccontaineresdD = 0x0;
1298 TH1D* efficiencysigTOFPIDD;
1299 TH1D* efficiencyTOFPIDD;
1300 TH1D* efficiencysigesdTOFPIDD;
1301 TH1D* efficiencyesdTOFPIDD;
1302 Int_t source = -1; //get parameterized TOF PID efficiencies
1303
1304 // signal sample
1305 mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1306 AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
1307 efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1308
1309 // mb sample for double check
1310 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1311 AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
1312 efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1313
1314 // mb sample with reconstructed variables
1315 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1316 AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD);
1317 efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1318
1319 // mb sample with reconstructed variables
1320 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1321 AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
1322 efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1323
1324 //fill histo
1325 efficiencysigTOFPIDD = (TH1D *) efficiencymcsigParamTOFPID->Project(0);
1326 efficiencyTOFPIDD = (TH1D *) efficiencymcParamTOFPID->Project(0);
1327 efficiencysigesdTOFPIDD = (TH1D *) efficiencysigesdParamTOFPID->Project(0);
1328 efficiencyesdTOFPIDD = (TH1D *) efficiencyesdParamTOFPID->Project(0);
1329 efficiencysigTOFPIDD->SetName("efficiencysigTOFPIDD");
1330 efficiencyTOFPIDD->SetName("efficiencyTOFPIDD");
1331 efficiencysigesdTOFPIDD->SetName("efficiencysigesdTOFPIDD");
1332 efficiencyesdTOFPIDD->SetName("efficiencyesdTOFPIDD");
1333
1334 //fit (mc enhenced sample)
1335 fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1336 efficiencysigTOFPIDD->Fit(fittofpid,"R");
1337 efficiencysigTOFPIDD->GetYaxis()->SetTitle("Efficiency");
1338 fEfficiencyTOFPIDD = efficiencysigTOFPIDD->GetFunction("fittofpid");
1339
1340 //fit (esd enhenced sample)
1341 efficiencysigesdTOFPIDD->Fit(fittofpid,"R");
1342 efficiencysigesdTOFPIDD->GetYaxis()->SetTitle("Efficiency");
1343 fEfficiencyesdTOFPIDD = efficiencysigesdTOFPIDD->GetFunction("fittofpid");
1344
1345 // draw (for PbPb, only 1st bin)
1346 //sig mc
1347 efficiencysigTOFPIDD->SetTitle("");
1348 efficiencysigTOFPIDD->SetStats(0);
1349 efficiencysigTOFPIDD->SetMarkerStyle(25);
1350 efficiencysigTOFPIDD->SetMarkerColor(2);
1351 efficiencysigTOFPIDD->SetLineColor(2);
1352 efficiencysigTOFPIDD->Draw();
1353
1354 //mb mc
1355 efficiencyTOFPIDD->SetTitle("");
1356 efficiencyTOFPIDD->SetStats(0);
1357 efficiencyTOFPIDD->SetMarkerStyle(24);
1358 efficiencyTOFPIDD->SetMarkerColor(4);
1359 efficiencyTOFPIDD->SetLineColor(4);
1360 efficiencyTOFPIDD->Draw("same");
1361
1362 //sig esd
1363 efficiencysigesdTOFPIDD->SetTitle("");
1364 efficiencysigesdTOFPIDD->SetStats(0);
1365 efficiencysigesdTOFPIDD->SetMarkerStyle(25);
1366 efficiencysigesdTOFPIDD->SetMarkerColor(3);
1367 efficiencysigesdTOFPIDD->SetLineColor(3);
1368 efficiencysigesdTOFPIDD->Draw("same");
1369
1370 //mb esd
1371 efficiencyesdTOFPIDD->SetTitle("");
1372 efficiencyesdTOFPIDD->SetStats(0);
1373 efficiencyesdTOFPIDD->SetMarkerStyle(25);
1374 efficiencyesdTOFPIDD->SetMarkerColor(1);
1375 efficiencyesdTOFPIDD->SetLineColor(1);
1376 efficiencyesdTOFPIDD->Draw("same");
1377
1378 //signal mc fit
1379 if(fEfficiencyTOFPIDD){
1380 fEfficiencyTOFPIDD->SetLineColor(2);
1381 fEfficiencyTOFPIDD->Draw("same");
1382 }
1383 //mb esd fit
1384 if(fEfficiencyesdTOFPIDD){
1385 fEfficiencyesdTOFPIDD->SetLineColor(3);
1386 fEfficiencyesdTOFPIDD->Draw("same");
1387 }
1388
1389 TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
1390 legtofeff->AddEntry(efficiencysigTOFPIDD,"TOF PID Step Efficiency","");
1391 legtofeff->AddEntry(efficiencysigTOFPIDD,"vs MC p_{t} for enhenced samples","p");
1392 legtofeff->AddEntry(efficiencyTOFPIDD,"vs MC p_{t} for mb samples","p");
1393 legtofeff->AddEntry(efficiencysigesdTOFPIDD,"vs esd p_{t} for enhenced samples","p");
1394 legtofeff->AddEntry(efficiencyesdTOFPIDD,"vs esd p_{t} for mb samples","p");
1395 legtofeff->Draw("same");
1396
1397
1398 TCanvas * cefficiencyParamIP = new TCanvas("efficiencyParamIP","efficiencyParamIP",500,500);
1399 cefficiencyParamIP->cd();
1400 gStyle->SetOptStat(0);
1401
1402 // IP cut efficiencies
1403 AliCFContainer *charmCombined = 0x0;
1404 AliCFContainer *beautyCombined = 0x0;
1405 AliCFContainer *beautyCombinedesd = 0x0;
1406
1407 source = 0; //charm enhenced
1408 mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1409 AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
1410 efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1411
1412 charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");
1413
1414 source = 1; //beauty enhenced
1415 mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1416 AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
1417 efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1418
1419 beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined");
1420
1421 mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1422 AliCFEffGrid* efficiencyBeautySigesd = new AliCFEffGrid("efficiencyBeautySigesd","",*mccontaineresdD);
1423 efficiencyBeautySigesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1424
1425 beautyCombinedesd = (AliCFContainer*)mccontaineresdD->Clone("beautyCombinedesd");
1426
1427 source = 0; //charm mb
1428 mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1429 AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
1430 efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1431
1432 charmCombined->Add(mccontainermcD);
1433 AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined);
1434 efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
1435
1436 source = 1; //beauty mb
1437 mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1438 AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
1439 efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1440
1441 beautyCombined->Add(mccontainermcD);
1442 AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined);
1443 efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
1444
1445 mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1446 AliCFEffGrid* efficiencyBeautyesd = new AliCFEffGrid("efficiencyBeautyesd","",*mccontaineresdD);
1447 efficiencyBeautyesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1448
1449 beautyCombinedesd->Add(mccontaineresdD);
1450 AliCFEffGrid* efficiencyBeautyCombinedesd = new AliCFEffGrid("efficiencyBeautyCombinedesd","",*beautyCombinedesd);
1451 efficiencyBeautyCombinedesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
1452
1453 source = 2; //conversion mb
1454 mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1455 AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
1456 efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1457
1458 source = 3; //non HFE except for the conversion mb
1459 mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1460 AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
1461 efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1462
1463 if(fIPEffCombinedSamples){printf("combined samples taken for beauty and charm\n");
1464 fEfficiencyCharmSigD = (TH1D*)efficiencyCharmCombined->Project(0); //signal enhenced + mb
1465 fEfficiencyBeautySigD = (TH1D*)efficiencyBeautyCombined->Project(0); //signal enhenced + mb
1466 fEfficiencyBeautySigesdD = (TH1D*)efficiencyBeautyCombinedesd->Project(0); //signal enhenced + mb
1467 }
1468 else{printf("signal enhanced samples taken for beauty and charm\n");
1469 fEfficiencyCharmSigD = (TH1D*)efficiencyCharmSig->Project(0); //signal enhenced only
1470 fEfficiencyBeautySigD = (TH1D*)efficiencyBeautySig->Project(0); //signal enhenced only
1471 fEfficiencyBeautySigesdD = (TH1D*)efficiencyBeautySigesd->Project(0); //signal enhenced only
1472 }
1473 fCharmEff = (TH1D*)efficiencyCharm->Project(0); //mb only
1474 fBeautyEff = (TH1D*)efficiencyBeauty->Project(0); //mb only
1475 fConversionEff = (TH1D*)efficiencyConv->Project(0); //mb only
1476 fNonHFEEff = (TH1D*)efficiencyNonhfe->Project(0); //mb only
1477
1478 if(fBeamType==0){
1479 //AliCFEffGrid *nonHFEEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
32c709ed 1480 AliCFContainer *cfcontainer = GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig);
1481 if(!cfcontainer) return;
1482 AliCFEffGrid *nonHFEEffGrid = (AliCFEffGrid*) GetEfficiency(cfcontainer,1,0); //mjmj
4bd4fc13 1483 fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
1484
1485 //AliCFEffGrid *conversionEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
32c709ed 1486 AliCFContainer *cfcontainerr = GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig);
1487 if(!cfcontainerr) return;
1488 AliCFEffGrid *conversionEffGrid = (AliCFEffGrid*) GetEfficiency(cfcontainerr,1,0); //mjmj
4bd4fc13 1489 fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
1490
1491 //MHMH
1492 if(fBeamType == 0){
1493 fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1494 fNonHFEEffbgc->Fit(fipfitnonhfe,"R");
1495 fEfficiencyIPNonhfeD = fNonHFEEffbgc->GetFunction("fipfitnonhfe");
1496
1497 fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
1498 fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1499 fConversionEffbgc->Fit(fipfitnonhfe,"R");
1500 fEfficiencyIPConversionD = fConversionEffbgc->GetFunction("fipfitnonhfe");
1501 }
1502 //MHMH
1503 }
1504
1505 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1506 fipfit->SetLineColor(2);
1507 fEfficiencyBeautySigD->Fit(fipfit,"R");
1508 fEfficiencyBeautySigD->GetYaxis()->SetTitle("Efficiency");
1509 fEfficiencyIPBeautyD = fEfficiencyBeautySigD->GetFunction("fipfit");
1510
1511 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1512 fipfit->SetLineColor(6);
1513 fEfficiencyBeautySigesdD->Fit(fipfit,"R");
1514 fEfficiencyBeautySigesdD->GetYaxis()->SetTitle("Efficiency");
1515 fEfficiencyIPBeautyesdD = fEfficiencyBeautySigesdD->GetFunction("fipfit");
1516
1517 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1518 fipfit->SetLineColor(1);
1519 fEfficiencyCharmSigD->Fit(fipfit,"R");
1520 fEfficiencyCharmSigD->GetYaxis()->SetTitle("Efficiency");
1521 fEfficiencyIPCharmD = fEfficiencyCharmSigD->GetFunction("fipfit");
1522
1523 //if(0){
1524 if(fIPParameterizedEff){
1525 // fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1526 fipfitnonhfe->SetLineColor(3);
1527 if(fBeamType==1){
1528 fConversionEff->Fit(fipfitnonhfe,"R");
1529 fConversionEff->GetYaxis()->SetTitle("Efficiency");
1530 fEfficiencyIPConversionD = fConversionEff->GetFunction("fipfitnonhfe");
1531 }
1532 else{
1533 fConversionEffbgc->Fit(fipfitnonhfe,"R");
1534 fConversionEffbgc->GetYaxis()->SetTitle("Efficiency");
1535 fEfficiencyIPConversionD = fConversionEffbgc->GetFunction("fipfitnonhfe");
1536 }
1537 // fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1538 fipfitnonhfe->SetLineColor(4);
1539 if(fBeamType==1){
1540 fNonHFEEff->Fit(fipfitnonhfe,"R");
1541 fNonHFEEff->GetYaxis()->SetTitle("Efficiency");
1542 fEfficiencyIPNonhfeD = fNonHFEEff->GetFunction("fipfitnonhfe");
1543 }
1544 else{
1545 fNonHFEEffbgc->Fit(fipfitnonhfe,"R");
1546 fNonHFEEffbgc->GetYaxis()->SetTitle("Efficiency");
1547 fEfficiencyIPNonhfeD = fNonHFEEffbgc->GetFunction("fipfitnonhfe");
1548 }
1549 }
1550
1551 // draw
1552 fEfficiencyCharmSigD->SetMarkerStyle(21);
1553 fEfficiencyCharmSigD->SetMarkerColor(1);
1554 fEfficiencyCharmSigD->SetLineColor(1);
1555 fEfficiencyBeautySigD->SetMarkerStyle(21);
1556 fEfficiencyBeautySigD->SetMarkerColor(2);
1557 fEfficiencyBeautySigD->SetLineColor(2);
1558 fEfficiencyBeautySigesdD->SetStats(0);
1559 fEfficiencyBeautySigesdD->SetMarkerStyle(21);
1560 fEfficiencyBeautySigesdD->SetMarkerColor(6);
1561 fEfficiencyBeautySigesdD->SetLineColor(6);
1562 fCharmEff->SetMarkerStyle(24);
1563 fCharmEff->SetMarkerColor(1);
1564 fCharmEff->SetLineColor(1);
1565 fBeautyEff->SetMarkerStyle(24);
1566 fBeautyEff->SetMarkerColor(2);
1567 fBeautyEff->SetLineColor(2);
1568 fConversionEff->SetMarkerStyle(24);
1569 fConversionEff->SetMarkerColor(3);
1570 fConversionEff->SetLineColor(3);
1571 fNonHFEEff->SetMarkerStyle(24);
1572 fNonHFEEff->SetMarkerColor(4);
1573 fNonHFEEff->SetLineColor(4);
1574
1575 fEfficiencyCharmSigD->Draw();
1576 fEfficiencyCharmSigD->GetXaxis()->SetRangeUser(0.0,7.9);
1577 fEfficiencyCharmSigD->GetYaxis()->SetRangeUser(0.0,0.5);
1578
1579 fEfficiencyBeautySigD->Draw("same");
1580 fEfficiencyBeautySigesdD->Draw("same");
1581 if(fBeamType == 1){
1582 fNonHFEEff->Draw("same");
1583 fConversionEff->Draw("same");
1584 }
1585 //fCharmEff->Draw("same");
1586 //fBeautyEff->Draw("same");
1587
1588 if(fBeamType==0){
1589 fConversionEffbgc->SetMarkerStyle(25);
1590 fConversionEffbgc->SetMarkerColor(3);
1591 fConversionEffbgc->SetLineColor(3);
1592 fNonHFEEffbgc->SetMarkerStyle(25);
1593 fNonHFEEffbgc->SetMarkerColor(4);
1594 fNonHFEEffbgc->SetLineColor(4);
1595 fConversionEffbgc->Draw("same");
1596 fNonHFEEffbgc->Draw("same");
1597 }
1598
1599 if(fEfficiencyIPBeautyD)
1600 fEfficiencyIPBeautyD->Draw("same");
1601 if(fEfficiencyIPBeautyesdD)
1602 fEfficiencyIPBeautyesdD->Draw("same");
1603 if( fEfficiencyIPCharmD)
1604 fEfficiencyIPCharmD->Draw("same");
1605 if(fIPParameterizedEff){
1606 if(fEfficiencyIPConversionD)
1607 fEfficiencyIPConversionD->Draw("same");
1608 if(fEfficiencyIPNonhfeD)
1609 fEfficiencyIPNonhfeD->Draw("same");
1610 }
1611 TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39);
1612 legipeff->AddEntry(fEfficiencyBeautySigD,"IP Step Efficiency","");
1613 legipeff->AddEntry(fEfficiencyBeautySigD,"beauty e","p");
1614 legipeff->AddEntry(fEfficiencyBeautySigesdD,"beauty e(esd pt)","p");
1615 legipeff->AddEntry(fEfficiencyCharmSigD,"charm e","p");
1616 if(fBeamType == 0){
1617 legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p");
1618 legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p");
1619 }
1620 else{
1621 legipeff->AddEntry(fConversionEff,"conversion e","p");
1622 legipeff->AddEntry(fNonHFEEff,"Dalitz e","p");
1623 }
1624 legipeff->Draw("same");
1625 gPad->SetGrid();
1626 //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps");
1627}
1628
1629//____________________________________________________________________________
1630THnSparse* AliHFEBeautySpectrum::GetBeautyIPEff(Bool_t isMCpt){
1631 //
1632 // Return beauty electron IP cut efficiency
1633 //
1634
1635 const Int_t nPtbinning1 = kSignalPtBins;//number of pt bins, according to new binning
1636 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
1637
1638 Int_t nDim=1; //dimensions of the efficiency weighting grid
1639
1640 Int_t nBin[1] = {nPtbinning1};
1641
1642 THnSparseF *ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
1643
1644 ipcut->SetBinEdges(0, kPtRange);
1645
1646 Double_t pt[1];
1647 Double_t weight;
1648 Double_t weightErr;
1649 Double_t contents[2];
1650
1651 weight = 1.0;
1652 weightErr = 1.0;
1653
1654 Int_t looppt=nBin[0];
1655
1656 Int_t ibin[2];
1657
1658 for(int i=0; i<looppt; i++)
1659 {
1660 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1661 if(isMCpt){
1662 if(fEfficiencyIPBeautyD){
1663 weight=fEfficiencyIPBeautyD->Eval(pt[0]);
1664 weightErr = 0;
1665 }
1666 else{
1667 printf("Fit failed on beauty IP cut efficiency. Contents in histo used!\n");
1668 weight = fEfficiencyBeautySigD->GetBinContent(i+1);
1669 weightErr = fEfficiencyBeautySigD->GetBinError(i+1);
1670 }
1671 }
1672 else{
1673 if(fEfficiencyIPBeautyesdD){
1674 weight=fEfficiencyIPBeautyesdD->Eval(pt[0]);
1675 weightErr = 0;
1676 }
1677 else{
1678 printf("Fit failed on beauty IP cut efficiency. Contents in histo used!\n");
1679 weight = fEfficiencyBeautySigesdD->GetBinContent(i+1);
1680 weightErr = fEfficiencyBeautySigD->GetBinError(i+1);
1681 }
1682 }
1683
1684 contents[0]=pt[0];
1685 ibin[0]=i+1;
1686
1687 ipcut->Fill(contents,weight);
1688 ipcut->SetBinError(ibin,weightErr);
1689 }
1690
1691 Int_t nDimSparse = ipcut->GetNdimensions();
1692 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1693 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1694
1695 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1696 binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
1697 nBins *= binsvar[iVar];
1698 }
1699
1700 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1701 // loop that sets 0 error in each bin
1702 for (Long_t iBin=0; iBin<nBins; iBin++) {
1703 Long_t bintmp = iBin ;
1704 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1705 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1706 bintmp /= binsvar[iVar] ;
1707 }
1708 //ipcut->SetBinError(binfill,0.); // put 0 everywhere
1709 }
1710
1711 delete[] binsvar;
1712 delete[] binfill;
1713
1714 return ipcut;
1715}
1716
1717//____________________________________________________________________________
1718THnSparse* AliHFEBeautySpectrum::GetPIDxIPEff(Int_t source){
1719 //
1720 // Return PID x IP cut efficiency
1721 //
1722 const Int_t nPtbinning1 = kSignalPtBins;//number of pt bins, according to new binning
1723 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
1724 Int_t nDim=1; //dimensions of the efficiency weighting grid
1725
1726 Int_t nBin[1] = {nPtbinning1};
1727
1728 THnSparseF *pideff;
1729 pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
1730 pideff->SetBinEdges(0, kPtRange);
1731
1732 Double_t pt[1];
1733 Double_t weight;
1734 Double_t weightErr;
1735 Double_t contents[2];
1736
1737 weight = 1.0;
1738 weightErr = 1.0;
1739
1740 Int_t looppt=nBin[0];
1741 Int_t ibin[2];
1742
1743 Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
1744 for(int i=0; i<looppt; i++)
1745 {
1746 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1747
1748 Double_t tofpideff = 1.;
1749 Double_t tofpideffesd = 1.;
1750 if(fEfficiencyTOFPIDD)
1751 tofpideff = fEfficiencyTOFPIDD->Eval(pt[0]);
1752 else{
1753 printf("TOF PID fit failed on conversion. The result is wrong!\n");
1754 }
1755 if(fEfficiencyesdTOFPIDD)
1756 tofpideffesd = fEfficiencyesdTOFPIDD->Eval(pt[0]);
1757 else{
1758 printf("TOF esd PID fit failed on conversion. The result is wrong!\n");
1759 }
1760
1761 //tof pid eff x tpc pid eff x ip cut eff
1762 if(fIPParameterizedEff){
1763 if(source==0) {
1764 if(fEfficiencyIPCharmD){
1765 weight = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD->Eval(pt[0]);
1766 weightErr = 0;
1767 }
1768 else{
1769 printf("Fit failed on charm IP cut efficiency\n");
1770 weight = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinContent(i+1);
1771 weightErr = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinError(i+1);
1772 }
1773 }
1774 else if(source==2) {
1775 if(fEfficiencyIPConversionD){
1776 weight = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD->Eval(pt[0]);
1777 weightErr = 0;
1778 }
1779 else{
1780 printf("Fit failed on conversion IP cut efficiency\n");
1781 weight = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1);
1782 weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinError(i+1);
1783 }
1784 }
1785 else if(source==3) {
1786 if(fEfficiencyIPNonhfeD){
1787 weight = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD->Eval(pt[0]);
1788 weightErr = 0;
1789 }
1790 else{
1791 printf("Fit failed on dalitz IP cut efficiency\n");
1792 weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1);
1793 weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinError(i+1);
1794 }
1795 }
1796 }
1797 else{
1798 if(source==0){
1799 if(fEfficiencyIPCharmD){
1800 weight = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD->Eval(pt[0]);
1801 weightErr = 0;
1802 }
1803 else{
1804 printf("Fit failed on charm IP cut efficiency\n");
1805 weight = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinContent(i+1);
1806 weightErr = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinError(i+1);
1807 }
1808 }
1809 else if(source==2){
1810 if(fBeamType==0){
1811 weight = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
1812 weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
1813 }
1814 else{
1815 weight = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1); // conversion
1816 weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinError(i+1);
1817 }
1818 }
1819 else if(source==3){
1820 if(fBeamType==0){
1821 weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
1822 weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
1823 }
1824 else{
1825 weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1); // Dalitz
1826 weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinError(i+1);
1827 }
1828 }
1829 }
1830
1831 contents[0]=pt[0];
1832 ibin[0]=i+1;
1833
1834 pideff->Fill(contents,weight);
1835 pideff->SetBinError(ibin,weightErr);
1836 }
1837
1838 Int_t nDimSparse = pideff->GetNdimensions();
1839 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1840 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1841
1842 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1843 binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
1844 nBins *= binsvar[iVar];
1845 }
1846
1847 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1848 // loop that sets 0 error in each bin
1849 for (Long_t iBin=0; iBin<nBins; iBin++) {
1850 Long_t bintmp = iBin ;
1851 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1852 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1853 bintmp /= binsvar[iVar] ;
1854 }
1855 }
1856
1857 delete[] binsvar;
1858 delete[] binfill;
1859
1860
1861 return pideff;
1862}
1863
1864//__________________________________________________________________________
1865AliCFDataGrid *AliHFEBeautySpectrum::GetRawBspectra2ndMethod(){
1866 //
1867 // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method
1868 //
1869 Int_t nDim = 1;
1870
1871 const Int_t nPtbinning1 = 18;//number of pt bins, according to new binning
1872 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};
1873
1874 Int_t nBin[1] = {nPtbinning1};
1875
1876 AliCFDataGrid *rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
1877
1878 THnSparseF *brawspectra;
1879 brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin);
1880
1881 brawspectra->SetBinEdges(0, kPtRange);
1882
1883 Double_t pt[1];
1884 Double_t yields= 0.;
1885 Double_t valuesb[2];
1886
1887 for(int i=0; i<fBSpectrum2ndMethod->GetNbinsX(); i++){
1888 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1889
1890 yields = fBSpectrum2ndMethod->GetBinContent(i+1);
1891
1892 valuesb[0]=pt[0];
1893 brawspectra->Fill(valuesb,yields);
1894 }
1895
1896
1897
1898 Int_t nDimSparse = brawspectra->GetNdimensions();
1899 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1900 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1901
1902 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1903 binsvar[iVar] = brawspectra->GetAxis(iVar)->GetNbins();
1904 nBins *= binsvar[iVar];
1905 }
1906
1907 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1908 // loop that sets 0 error in each bin
1909 for (Long_t iBin=0; iBin<nBins; iBin++) {
1910 Long_t bintmp = iBin ;
1911 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1912 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1913 bintmp /= binsvar[iVar] ;
1914 }
1915 brawspectra->SetBinError(binfill,0.); // put 0 everywhere
1916 }
1917
1918
1919 rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency
1920 TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(0);
1921
1922 new TCanvas;
1923 fBSpectrum2ndMethod->SetMarkerStyle(24);
1924 fBSpectrum2ndMethod->Draw("p");
1925 hRawBeautySpectra->SetMarkerStyle(25);
1926 hRawBeautySpectra->Draw("samep");
1927
1928 delete[] binfill;
1929 delete[] binsvar;
1930
1931 return rawBeautyContainer;
1932}
1933//__________________________________________________________________________
1934void AliHFEBeautySpectrum::CalculateNonHFEsyst(){
1935 //
1936 // Calculate non HFE sys
1937 //
1938 //
1939
1940 if(!fNonHFEsyst)
1941 return;
1942
1943 Double_t evtnorm[1] = {0.0};
1944 if(fNMCbgEvents>0) evtnorm[0]= double(fNEvents)/double(fNMCbgEvents);
1945
1946 AliCFDataGrid *convSourceGrid[kElecBgSources-3][kBgLevels];
1947 AliCFDataGrid *nonHFESourceGrid[kElecBgSources-3][kBgLevels];
1948
1949 AliCFDataGrid *bgLevelGrid[2][kBgLevels];//for pi0 and eta based errors
1950 AliCFDataGrid *bgNonHFEGrid[kBgLevels];
1951 AliCFDataGrid *bgConvGrid[kBgLevels];
1952
1953 Int_t stepbackground = 3;
1954 Int_t* bins=new Int_t[1];
1955 const Char_t *bgBase[2] = {"pi0","eta"};
1956
1957 bins[0]=kSignalPtBins;//fConversionEff[centrality]->GetNbinsX();
1958
1959 AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
1960 AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
1961
1962 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1963 for(Int_t iSource = 0; iSource < kElecBgSources-3; iSource++){
1964 convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d",iSource,iLevel),Form("convGrid_%d_%d",iSource,iLevel),*fConvSourceContainer[iSource][iLevel],stepbackground);
1965 weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
1966 convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
1967
1968 nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d",iSource,iLevel),Form("nonHFEGrid_%d_%d",iSource,iLevel),*fNonHFESourceContainer[iSource][iLevel],stepbackground);
1969 weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1970 nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
1971 }
1972
1973 bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
1974 for(Int_t iSource = 2; iSource < kElecBgSources-3; iSource++){
1975 bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
1976 }
1977 if(!fEtaSyst)
1978 bgConvGrid[iLevel]->Add(convSourceGrid[1][iLevel]);
1979
1980 bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone();
1981 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)
1982 bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
1983 }
1984 if(!fEtaSyst)
1985 bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[1][iLevel]);
1986
1987 bgLevelGrid[0][iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
1988 bgLevelGrid[0][iLevel]->Add(bgNonHFEGrid[iLevel]);
1989 if(fEtaSyst){
1990 bgLevelGrid[1][iLevel] = (AliCFDataGrid*)nonHFESourceGrid[1][iLevel]->Clone();//background for eta source
1991 bgLevelGrid[1][iLevel]->Add(convSourceGrid[1][iLevel]);
1992 }
1993 }
1994
1995 //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)
1996 AliCFDataGrid *bgErrorGrid[2][2];//for pions/eta error base, for lower/upper
1997 TH1D* hBaseErrors[2][2];//pi0/eta and lower/upper
1998 for(Int_t iErr = 0; iErr < 2; iErr++){//errors for pi0 and eta base
1999 bgErrorGrid[iErr][0] = (AliCFDataGrid*)bgLevelGrid[iErr][1]->Clone();
2000 bgErrorGrid[iErr][0]->Add(bgLevelGrid[iErr][0],-1.);
2001 bgErrorGrid[iErr][1] = (AliCFDataGrid*)bgLevelGrid[iErr][2]->Clone();
2002 bgErrorGrid[iErr][1]->Add(bgLevelGrid[iErr][0],-1.);
2003
2004 //plot absolute differences between limit yields (upper/lower limit, based on pi0 and eta errors) and best estimate
2005
2006 hBaseErrors[iErr][0] = (TH1D*)bgErrorGrid[iErr][0]->Project(0);
2007 hBaseErrors[iErr][0]->Scale(-1.);
2008 hBaseErrors[iErr][0]->SetTitle(Form("Absolute %s-based systematic errors from non-HF meson decays and conversions",bgBase[iErr]));
2009 hBaseErrors[iErr][1] = (TH1D*)bgErrorGrid[iErr][1]->Project(0);
2010 if(!fEtaSyst)break;
2011 }
2012
2013 //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
2014 TH1D *hSpeciesErrors[kElecBgSources-4];
2015 for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
2016 if(fEtaSyst && (iSource == 1))continue;
2017 hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0);
2018 TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0);
2019 hSpeciesErrors[iSource-1]->Add(hNonHFEtemp);
2020 hSpeciesErrors[iSource-1]->Scale(0.3);
2021 }
2022
2023 //Int_t firstBgSource = 0;//if eta systematics are not from scaling
2024 //if(fEtaSyst){firstBgSource = 1;}//source 0 histograms are not filled if eta errors are independently determined!
2025 TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[1]->Clone();
2026 TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[1]->Clone();
2027 TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[1]->Clone();
2028
2029 TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
2030 TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
2031
2032 for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
2033 Double_t pi0basedErrLow,pi0basedErrUp,etaErrLow,etaErrUp;
2034 pi0basedErrLow = hBaseErrors[0][0]->GetBinContent(iBin);
2035 pi0basedErrUp = hBaseErrors[0][1]->GetBinContent(iBin);
2036 if(fEtaSyst){
2037 etaErrLow = hBaseErrors[1][0]->GetBinContent(iBin);
2038 etaErrUp = hBaseErrors[1][1]->GetBinContent(iBin);
2039 }
2040 else{ etaErrLow = etaErrUp = 0.;}
2041
2042 Double_t sqrsumErrs= 0;
2043 for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
2044 if(fEtaSyst && (iSource == 1))continue;
2045 Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
2046 sqrsumErrs+=(scalingErr*scalingErr);
2047 }
2048 for(Int_t iErr = 0; iErr < 2; iErr++){
2049 for(Int_t iLevel = 0; iLevel < 2; iLevel++){
2050 hBaseErrors[iErr][iLevel]->SetBinContent(iBin,hBaseErrors[iErr][iLevel]->GetBinContent(iBin)/hBaseErrors[iErr][iLevel]->GetBinWidth(iBin));
2051 }
2052 if(!fEtaSyst)break;
2053 }
2054 hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+(etaErrUp*etaErrUp)+sqrsumErrs));
2055 hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+(etaErrLow*etaErrLow)+sqrsumErrs));
2056 hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin));
2057
2058 hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin));
2059 hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin));
2060 }
2061
2062 TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
2063 cPiErrors->cd();
2064 cPiErrors->SetLogx();
2065 cPiErrors->SetLogy();
2066 hBaseErrors[0][0]->Draw();
2067 //hBaseErrors[0][1]->SetMarkerColor(kBlack);
2068 //hBaseErrors[0][1]->SetLineColor(kBlack);
2069 //hBaseErrors[0][1]->Draw("SAME");
2070 if(fEtaSyst){
2071 hBaseErrors[1][0]->Draw("SAME");
2072 hBaseErrors[1][0]->SetMarkerColor(kBlack);
2073 hBaseErrors[1][0]->SetLineColor(kBlack);
2074 //hBaseErrors[1][1]->SetMarkerColor(13);
2075 //hBaseErrors[1][1]->SetLineColor(13);
2076 //hBaseErrors[1][1]->Draw("SAME");
2077 }
2078 //hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
2079 //hOverallBinScaledErrsUp->SetLineColor(kBlue);
2080 //hOverallBinScaledErrsUp->Draw("SAME");
2081 hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
2082 hOverallBinScaledErrsLow->SetLineColor(kGreen);
2083 hOverallBinScaledErrsLow->Draw("SAME");
2084 hScalingErrors->SetLineColor(kBlue);
2085 hScalingErrors->Draw("SAME");
2086
2087 TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
2088 lPiErr->AddEntry(hBaseErrors[0][0],"Lower error from pion error");
2089 //lPiErr->AddEntry(hBaseErrors[0][1],"Upper error from pion error");
2090 if(fEtaSyst){
2091 lPiErr->AddEntry(hBaseErrors[1][0],"Lower error from eta error");
2092 //lPiErr->AddEntry(hBaseErrors[1][1],"Upper error from eta error");
2093 }
2094 lPiErr->AddEntry(hScalingErrors, "scaling error");
2095 lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics");
2096 //lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
2097 lPiErr->Draw("SAME");
2098
2099 //Normalize errors
2100 TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone();
2101 TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone();
2102 hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!)
2103 hLowSystScaled->Scale(evtnorm[0]);
2104 TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone();
2105 TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone();
2106 //histograms to be normalized to TGraphErrors
2107 CorrectFromTheWidth(hNormAllSystErrUp);
2108 CorrectFromTheWidth(hNormAllSystErrLow);
2109
2110 TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs");
2111 cNormOvErrs->cd();
2112 cNormOvErrs->SetLogx();
2113 cNormOvErrs->SetLogy();
2114
2115 TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp);
2116 TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow);
2117 gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors");
2118 gOverallSystErrUp->SetMarkerColor(kBlack);
2119 gOverallSystErrUp->SetLineColor(kBlack);
2120 gOverallSystErrLow->SetMarkerColor(kRed);
2121 gOverallSystErrLow->SetLineColor(kRed);
2122 gOverallSystErrUp->Draw("AP");
2123 gOverallSystErrLow->Draw("PSAME");
2124 TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89);
2125 lAllSys->AddEntry(gOverallSystErrLow,"lower","p");
2126 lAllSys->AddEntry(gOverallSystErrUp,"upper","p");
2127 lAllSys->Draw("same");
2128
2129
2130 AliCFDataGrid *bgYieldGrid;
2131 if(fEtaSyst){
2132 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.
2133 }
2134 bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0][0]->Clone();
2135
2136 TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0);
2137 TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone();
2138 hRelErrUp->Divide(hBgYield);
2139 TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone();
2140 hRelErrLow->Divide(hBgYield);
2141
2142 TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs");
2143 cRelErrs->cd();
2144 cRelErrs->SetLogx();
2145 hRelErrUp->SetTitle("Relative error of non-HFE background yield");
2146 hRelErrUp->Draw();
2147 hRelErrLow->SetLineColor(kBlack);
2148 hRelErrLow->Draw("SAME");
2149
2150 TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95);
2151 lRel->AddEntry(hRelErrUp, "upper");
2152 lRel->AddEntry(hRelErrLow, "lower");
2153 lRel->Draw("SAME");
2154
2155 //CorrectFromTheWidth(hBgYield);
2156 //hBgYield->Scale(evtnorm[0]);
2157
2158
2159 //write histograms/TGraphs to file
2160 TFile *output = new TFile("systHists.root","recreate");
2161
2162 hBgYield->SetName("hBgYield");
2163 hBgYield->Write();
2164 hRelErrUp->SetName("hRelErrUp");
2165 hRelErrUp->Write();
2166 hRelErrLow->SetName("hRelErrLow");
2167 hRelErrLow->Write();
2168 hUpSystScaled->SetName("hOverallSystErrUp");
2169 hUpSystScaled->Write();
2170 hLowSystScaled->SetName("hOverallSystErrLow");
2171 hLowSystScaled->Write();
2172 gOverallSystErrUp->SetName("gOverallSystErrUp");
2173 gOverallSystErrUp->Write();
2174 gOverallSystErrLow->SetName("gOverallSystErrLow");
2175 gOverallSystErrLow->Write();
2176
2177 output->Close();
2178 delete output;
2179}
2180//____________________________________________________________
2181void AliHFEBeautySpectrum::UnfoldBG(AliCFDataGrid* const bgsubpectrum){
2182 //
2183 // Unfold backgrounds to check its sanity
2184 //
2185
2186 AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
2187 //AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
2188 if(!mcContainer){
2189 AliError("MC Container not available");
2190 return;
2191 }
2192
2193 if(!fCorrelation){
2194 AliError("No Correlation map available");
2195 return;
2196 }
2197
2198 // Data
2199 AliCFDataGrid *dataGrid = 0x0;
2200 dataGrid = bgsubpectrum;
2201
2202 // Guessed
2203 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
2204 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
2205
2206 // Efficiency
2207 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
2208 efficiencyD->CalculateEfficiency(fStepMC+2,fStepTrue);
2209
2210 // Unfold background spectra
2211 Int_t nDim=1;
2212 AliCFUnfolding unfolding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
2213 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
2214 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
2215 if(fSetSmoothing) unfolding.UseSmoothing();
2216 unfolding.Unfold();
2217
2218 // Results
2219 THnSparse* result = unfolding.GetUnfolded();
2220 TCanvas *ctest = new TCanvas("yvonnetest","yvonnetest",1000,600);
2221 if(fBeamType==1)
2222 {
2223 ctest->Divide(2);
2224 ctest->cd(1);
2225 TH1D* htest1=(TH1D*)result->Projection(0);
2226 htest1->Draw();
2227 ctest->cd(2);
2228 TH1D* htest2=(TH1D*)result->Projection(1);
2229 htest2->Draw();
2230 }
2231
2232
2233 TGraphErrors* unfoldedbgspectrumD = Normalize(result);
2234 if(!unfoldedbgspectrumD) {
2235 AliError("Unfolded background spectrum doesn't exist");
2236 }
2237 else{
2238 TFile *file = TFile::Open("unfoldedbgspectrum.root","recreate");
2239 if(fBeamType==0)unfoldedbgspectrumD->Write("unfoldedbgspectrum");
2240
2241 file->Close();
2242 }
2243}