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