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