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