]>
Commit | Line | Data |
---|---|---|
0a173978 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id$ */ | |
17 | ||
18 | // This class is used to store correction maps, raw input and results of the multiplicity | |
19 | // measurement with the ITS or TPC | |
20 | // It also contains functions to correct the spectrum using different methods. | |
6d81c2de | 21 | // e.g. chi2 minimization and bayesian unfolding |
0a173978 | 22 | // |
23 | // Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch | |
24 | ||
25 | #include "AliMultiplicityCorrection.h" | |
26 | ||
27 | #include <TFile.h> | |
28 | #include <TH1F.h> | |
29 | #include <TH2F.h> | |
30 | #include <TH3F.h> | |
31 | #include <TDirectory.h> | |
0a173978 | 32 | #include <TCanvas.h> |
33 | #include <TString.h> | |
9ca6be09 | 34 | #include <TF1.h> |
d560b581 | 35 | #include <TMath.h> |
3602328d | 36 | #include <TCollection.h> |
447c325d | 37 | #include <TLegend.h> |
38 | #include <TLine.h> | |
0b4bfd98 | 39 | #include <TRandom.h> |
40 | #include <TProfile.h> | |
41 | #include <TProfile2D.h> | |
eb9356d5 | 42 | #include <AliLog.h> |
0a173978 | 43 | |
44 | ClassImp(AliMultiplicityCorrection) | |
45 | ||
eb9356d5 | 46 | // Defined where the efficiency drops below 1/3 |
47 | // |eta| < 1.4 --> -0.3 ... 0.8 | |
48 | // |eta| < 1.3 --> -1.9 ... 2.4 | |
49 | // |eta| < 1.0 --> -5.6 ... 6.1 | |
50 | //Double_t AliMultiplicityCorrection::fgVtxRangeBegin[kESDHists] = { -10.0, -5.6, -1.9 }; | |
51 | //Double_t AliMultiplicityCorrection::fgVtxRangeEnd[kESDHists] = { 10.0, 6.1, 2.4 }; | |
52 | Double_t AliMultiplicityCorrection::fgVtxRangeBegin[kESDHists] = { -10.0, -5.5, -1.9 }; | |
53 | Double_t AliMultiplicityCorrection::fgVtxRangeEnd[kESDHists] = { 10.0, 5.5, 2.4 }; | |
54 | ||
0b4bfd98 | 55 | // These are the areas where the quality of the unfolding results are evaluated |
56 | // Default defined here, call SetQualityRegions to change them | |
57 | // unit is in multiplicity (not in bin!) | |
0b4bfd98 | 58 | // SPD: peak area - flat area - low stat area |
69b09e3b | 59 | Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {1, 20, 70}; |
60 | Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {10, 65, 80}; | |
0b4bfd98 | 61 | |
62 | //____________________________________________________________________ | |
63 | void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy) | |
64 | { | |
65 | // | |
66 | // sets the quality region definition to TPC or SPD | |
67 | // | |
68 | ||
69 | if (SPDStudy) | |
70 | { | |
71 | // SPD: peak area - flat area - low stat area | |
69b09e3b | 72 | fgQualityRegionsB[0] = 1; |
73 | fgQualityRegionsE[0] = 10; | |
0b4bfd98 | 74 | |
69b09e3b | 75 | fgQualityRegionsB[1] = 20; |
76 | fgQualityRegionsE[1] = 65; | |
0b4bfd98 | 77 | |
69b09e3b | 78 | fgQualityRegionsB[2] = 70; |
79 | fgQualityRegionsE[2] = 80; | |
6d81c2de | 80 | |
81 | Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for SPD"); | |
0b4bfd98 | 82 | } |
83 | else | |
84 | { | |
85 | // TPC: peak area - flat area - low stat area | |
86 | fgQualityRegionsB[0] = 4; | |
87 | fgQualityRegionsE[0] = 12; | |
88 | ||
89 | fgQualityRegionsB[1] = 25; | |
90 | fgQualityRegionsE[1] = 55; | |
91 | ||
92 | fgQualityRegionsB[2] = 88; | |
93 | fgQualityRegionsE[2] = 108; | |
6d81c2de | 94 | |
95 | Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for TPC"); | |
0b4bfd98 | 96 | } |
97 | } | |
98 | ||
0a173978 | 99 | //____________________________________________________________________ |
100 | AliMultiplicityCorrection::AliMultiplicityCorrection() : | |
d29b31c5 | 101 | TNamed(), fCurrentESD(0), fCurrentCorrelation(0), fCurrentEfficiency(0), fLastBinLimit(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0), fVtxBegin(0), fVtxEnd(0) |
0a173978 | 102 | { |
103 | // | |
104 | // default constructor | |
105 | // | |
106 | ||
107 | for (Int_t i = 0; i < kESDHists; ++i) | |
eb9356d5 | 108 | { |
0a173978 | 109 | fMultiplicityESD[i] = 0; |
eb9356d5 | 110 | fTriggeredEvents[i] = 0; |
111 | fNoVertexEvents[i] = 0; | |
112 | } | |
0a173978 | 113 | |
114 | for (Int_t i = 0; i < kMCHists; ++i) | |
cfc19dd5 | 115 | { |
116 | fMultiplicityVtx[i] = 0; | |
117 | fMultiplicityMB[i] = 0; | |
118 | fMultiplicityINEL[i] = 0; | |
69b09e3b | 119 | fMultiplicityNSD[i] = 0; |
cfc19dd5 | 120 | } |
0a173978 | 121 | |
122 | for (Int_t i = 0; i < kCorrHists; ++i) | |
123 | { | |
124 | fCorrelation[i] = 0; | |
125 | fMultiplicityESDCorrected[i] = 0; | |
126 | } | |
6d81c2de | 127 | |
128 | for (Int_t i = 0; i < kQualityRegions; ++i) | |
129 | fQuality[i] = 0; | |
0a173978 | 130 | } |
131 | ||
69b09e3b | 132 | //____________________________________________________________________ |
133 | AliMultiplicityCorrection* AliMultiplicityCorrection::Open(const char* fileName, const char* folderName) | |
134 | { | |
135 | // opens the given file, reads the multiplicity from the given folder and returns the object | |
136 | ||
137 | TFile* file = TFile::Open(fileName); | |
138 | if (!file) | |
139 | { | |
140 | Printf("ERROR: Could not open %s", fileName); | |
141 | return 0; | |
142 | } | |
143 | ||
eb9356d5 | 144 | Printf("AliMultiplicityCorrection::Open: Reading file %s", fileName); |
145 | ||
69b09e3b | 146 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folderName, folderName); |
147 | mult->LoadHistograms(); | |
148 | ||
149 | // TODO closing the file does not work here, because the histograms cannot be read anymore. LoadHistograms need to be adapted | |
150 | ||
151 | return mult; | |
152 | } | |
153 | ||
0a173978 | 154 | //____________________________________________________________________ |
155 | AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) : | |
7307d52c | 156 | TNamed(name, title), |
157 | fCurrentESD(0), | |
158 | fCurrentCorrelation(0), | |
159 | fCurrentEfficiency(0), | |
160 | fLastBinLimit(0), | |
161 | fLastChi2MC(0), | |
162 | fLastChi2MCLimit(0), | |
163 | fLastChi2Residuals(0), | |
eb9356d5 | 164 | fRatioAverage(0), |
165 | fVtxBegin(0), | |
166 | fVtxEnd(0) | |
0a173978 | 167 | { |
168 | // | |
169 | // named constructor | |
170 | // | |
eb9356d5 | 171 | |
0a173978 | 172 | // do not add this hists to the directory |
173 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
174 | TH1::AddDirectory(kFALSE); | |
175 | ||
b477f4f2 | 176 | /*Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10}; |
0a173978 | 177 | Float_t binLimitsN[] = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, |
178 | 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, | |
179 | 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, | |
180 | 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, | |
181 | 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, | |
182 | 50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5, | |
183 | 100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5, | |
184 | 150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5, | |
185 | 250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5, | |
9ca6be09 | 186 | 500.5 }; |
187 | //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5, | |
188 | //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5, | |
189 | //1000.5 }; | |
0a173978 | 190 | |
b477f4f2 | 191 | #define VTXBINNING 10, binLimitsVtx |
6d81c2de | 192 | #define NBINNING fgkMaxParams, binLimitsN*/ |
b477f4f2 | 193 | |
69b09e3b | 194 | #define NBINNING 201, -0.5, 200.5 |
b3b260d1 | 195 | |
0a173978 | 196 | for (Int_t i = 0; i < kESDHists; ++i) |
eb9356d5 | 197 | { |
198 | fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;measured multiplicity;Count", 1, fgVtxRangeBegin[i], fgVtxRangeEnd[i], NBINNING); | |
199 | fTriggeredEvents[i] = new TH1F(Form("fTriggeredEvents%d", i), "fTriggeredEvents;measured multiplicity;Count", NBINNING); | |
200 | fNoVertexEvents[i] = new TH1F(Form("fNoVertexEvents%d", i), "fNoVertexEvents;generated multiplicity;Count", NBINNING); | |
201 | } | |
0a173978 | 202 | |
203 | for (Int_t i = 0; i < kMCHists; ++i) | |
204 | { | |
b3b260d1 | 205 | fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[i%3]->Clone(Form("fMultiplicityVtx%d", i))); |
cfc19dd5 | 206 | fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart"); |
207 | ||
b3b260d1 | 208 | fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[i]->Clone(Form("fMultiplicityMB%d", i))); |
cfc19dd5 | 209 | fMultiplicityMB[i]->SetTitle("fMultiplicityMB"); |
210 | ||
b3b260d1 | 211 | fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[i]->Clone(Form("fMultiplicityINEL%d", i))); |
cfc19dd5 | 212 | fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL"); |
69b09e3b | 213 | |
b3b260d1 | 214 | fMultiplicityNSD[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[i]->Clone(Form("fMultiplicityNSD%d", i))); |
69b09e3b | 215 | fMultiplicityNSD[i]->SetTitle("fMultiplicityNSD"); |
0a173978 | 216 | } |
217 | ||
218 | for (Int_t i = 0; i < kCorrHists; ++i) | |
219 | { | |
eb9356d5 | 220 | fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;true multiplicity;measured multiplicity", 1, fgVtxRangeBegin[i%3], fgVtxRangeEnd[i%3], NBINNING, NBINNING); |
221 | fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;true multiplicity;Count", NBINNING); | |
0a173978 | 222 | } |
223 | ||
b9b1eb8a | 224 | for (Int_t i = 0; i < kQualityRegions; ++i) |
225 | fQuality[i] = 0; | |
226 | ||
0a173978 | 227 | TH1::AddDirectory(oldStatus); |
19442b86 | 228 | |
229 | AliUnfolding::SetNbins(120, 120); | |
230 | AliUnfolding::SetSkipBinsBegin(1); | |
eb9356d5 | 231 | //AliUnfolding::SetNormalizeInput(kTRUE); |
232 | } | |
233 | ||
234 | //____________________________________________________________________ | |
235 | void AliMultiplicityCorrection::Rebin2DY(TH2F*& hist, Int_t nBins, Double_t* newBins) const | |
236 | { | |
237 | // | |
238 | // rebins the y axis of a two-dimensional histogram giving variable size binning (missing in ROOT v5/25/02) | |
239 | // | |
240 | ||
241 | TH2F* temp = new TH2F(hist->GetName(), Form("%s;%s;%s", hist->GetTitle(), hist->GetXaxis()->GetTitle(), hist->GetYaxis()->GetTitle()), hist->GetNbinsX(), hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax(), nBins, newBins); | |
242 | ||
243 | for (Int_t x=0; x<=hist->GetNbinsX()+1; x++) | |
244 | for (Int_t y=0; y<=hist->GetNbinsY()+1; y++) | |
245 | temp->Fill(hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetBinContent(x, y)); | |
246 | ||
247 | for (Int_t x=0; x<=temp->GetNbinsX()+1; x++) | |
248 | for (Int_t y=0; y<=temp->GetNbinsY()+1; y++) | |
249 | temp->SetBinError(x, y, TMath::Sqrt(temp->GetBinContent(x, y))); | |
250 | ||
251 | delete hist; | |
252 | hist = temp; | |
253 | } | |
254 | ||
255 | //____________________________________________________________________ | |
256 | void AliMultiplicityCorrection::Rebin3DY(TH3F*& hist, Int_t nBins, Double_t* newBins) const | |
257 | { | |
258 | // | |
259 | // rebins the y axis of a three-dimensional histogram giving variable size binning (missing in ROOT v5/25/02) | |
260 | // this function is a mess - and it should have been Fons who should have gone through the pain of writing it! (JF) | |
261 | // | |
262 | ||
263 | // construct variable size arrays for fixed size binning axes because TH3 lacks some constructors | |
264 | Double_t* xBins = new Double_t[hist->GetNbinsX()+1]; | |
265 | Double_t* zBins = new Double_t[hist->GetNbinsZ()+1]; | |
266 | ||
267 | for (Int_t x=1; x<=hist->GetNbinsX()+1; x++) | |
268 | { | |
269 | xBins[x-1] = hist->GetXaxis()->GetBinLowEdge(x); | |
270 | //Printf("%d %f", x, xBins[x-1]); | |
271 | } | |
272 | ||
273 | for (Int_t z=1; z<=hist->GetNbinsZ()+1; z++) | |
274 | { | |
275 | zBins[z-1] = hist->GetZaxis()->GetBinLowEdge(z); | |
276 | //Printf("%d %f", y, yBins[y-1]); | |
277 | } | |
278 | ||
279 | TH3F* temp = new TH3F(hist->GetName(), Form("%s;%s;%s;%s", hist->GetTitle(), hist->GetXaxis()->GetTitle(), hist->GetYaxis()->GetTitle(), hist->GetZaxis()->GetTitle()), hist->GetNbinsX(), xBins, nBins, newBins, hist->GetNbinsZ(), zBins); | |
280 | ||
281 | for (Int_t x=0; x<=hist->GetNbinsX()+1; x++) | |
282 | for (Int_t y=0; y<=hist->GetNbinsY()+1; y++) | |
283 | for (Int_t z=0; z<=hist->GetNbinsZ()+1; z++) | |
284 | temp->Fill(hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z), hist->GetBinContent(x, y, z)); | |
285 | ||
286 | for (Int_t x=0; x<=temp->GetNbinsX()+1; x++) | |
287 | for (Int_t y=0; y<=temp->GetNbinsY()+1; y++) | |
288 | for (Int_t z=0; z<=hist->GetNbinsZ()+1; z++) | |
289 | temp->SetBinError(x, y, z, TMath::Sqrt(temp->GetBinContent(x, y, z))); | |
290 | ||
291 | delete[] xBins; | |
292 | delete[] zBins; | |
293 | ||
294 | delete hist; | |
295 | hist = temp; | |
296 | } | |
297 | ||
298 | //____________________________________________________________________ | |
299 | void AliMultiplicityCorrection::RebinGenerated(Int_t nBins05, Double_t* newBins05, Int_t nBins10, Double_t* newBins10, Int_t nBins13, Double_t* newBins13) | |
300 | { | |
301 | // | |
302 | // Rebins the (and only the) generated multiplicity axis | |
303 | // | |
304 | ||
305 | Printf("Rebinning generated-multiplicity axis..."); | |
306 | ||
307 | // do not add this hists to the directory | |
308 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
309 | TH1::AddDirectory(kFALSE); | |
310 | ||
311 | if (kESDHists != 3) | |
312 | AliFatal("This function only works for three ESD hists!"); | |
313 | ||
314 | for (Int_t i = 0; i < kESDHists; ++i) | |
315 | { | |
316 | Int_t nBins = -1; | |
317 | Double_t* newBins = 0; | |
318 | ||
319 | switch (i) | |
320 | { | |
321 | case 0: | |
322 | nBins = nBins05; | |
323 | newBins = newBins05; | |
324 | break; | |
325 | case 1: | |
326 | nBins = nBins10; | |
327 | newBins = newBins10; | |
328 | break; | |
329 | case 2: | |
330 | nBins = nBins13; | |
331 | newBins = newBins13; | |
332 | break; | |
333 | } | |
334 | ||
335 | // 1D | |
336 | // TODO mem leak | |
337 | fNoVertexEvents[i] = (TH1F*) fNoVertexEvents[i]->Rebin(nBins, fNoVertexEvents[i]->GetName(), newBins); | |
338 | fMultiplicityESDCorrected[i] = (TH1F*) fMultiplicityESDCorrected[i]->Rebin(nBins, fMultiplicityESDCorrected[i]->GetName(), newBins); | |
339 | ||
340 | // 2D | |
341 | Rebin2DY(fMultiplicityVtx[i], nBins, newBins); | |
342 | Rebin2DY(fMultiplicityMB[i], nBins, newBins); | |
343 | Rebin2DY(fMultiplicityINEL[i], nBins, newBins); | |
344 | Rebin2DY(fMultiplicityNSD[i], nBins, newBins); | |
345 | ||
346 | // 3D | |
347 | Rebin3DY(fCorrelation[i], nBins, newBins); | |
348 | } | |
349 | ||
350 | TH1::AddDirectory(oldStatus); | |
0a173978 | 351 | } |
352 | ||
353 | //____________________________________________________________________ | |
354 | AliMultiplicityCorrection::~AliMultiplicityCorrection() | |
355 | { | |
356 | // | |
357 | // Destructor | |
358 | // | |
6d81c2de | 359 | |
0f67a57c | 360 | Printf("AliMultiplicityCorrection::~AliMultiplicityCorrection called"); |
361 | ||
6d81c2de | 362 | for (Int_t i = 0; i < kESDHists; ++i) |
363 | { | |
364 | if (fMultiplicityESD[i]) | |
365 | delete fMultiplicityESD[i]; | |
366 | fMultiplicityESD[i] = 0; | |
eb9356d5 | 367 | |
368 | if (fTriggeredEvents[i]) | |
369 | delete fTriggeredEvents[i]; | |
370 | fTriggeredEvents[i]= 0; | |
371 | ||
372 | if (fNoVertexEvents[i]) | |
373 | delete fNoVertexEvents[i]; | |
374 | fNoVertexEvents[i]= 0; | |
6d81c2de | 375 | } |
376 | ||
377 | for (Int_t i = 0; i < kMCHists; ++i) | |
378 | { | |
379 | if (fMultiplicityVtx[i]) | |
380 | delete fMultiplicityVtx[i]; | |
381 | fMultiplicityVtx[i] = 0; | |
382 | ||
383 | if (fMultiplicityMB[i]) | |
384 | delete fMultiplicityMB[i]; | |
385 | fMultiplicityMB[i] = 0; | |
386 | ||
387 | if (fMultiplicityINEL[i]) | |
388 | delete fMultiplicityINEL[i]; | |
389 | fMultiplicityINEL[i] = 0; | |
69b09e3b | 390 | |
391 | if (fMultiplicityNSD[i]) | |
392 | delete fMultiplicityNSD[i]; | |
393 | fMultiplicityNSD[i] = 0; | |
394 | } | |
6d81c2de | 395 | |
396 | for (Int_t i = 0; i < kCorrHists; ++i) | |
397 | { | |
398 | if (fCorrelation[i]) | |
399 | delete fCorrelation[i]; | |
400 | fCorrelation[i] = 0; | |
401 | ||
402 | if (fMultiplicityESDCorrected[i]) | |
403 | delete fMultiplicityESDCorrected[i]; | |
404 | fMultiplicityESDCorrected[i] = 0; | |
405 | } | |
0a173978 | 406 | } |
407 | ||
408 | //____________________________________________________________________ | |
0427591c | 409 | Long64_t AliMultiplicityCorrection::Merge(const TCollection* list) |
0a173978 | 410 | { |
411 | // Merge a list of AliMultiplicityCorrection objects with this (needed for | |
412 | // PROOF). | |
413 | // Returns the number of merged objects (including this). | |
414 | ||
415 | if (!list) | |
416 | return 0; | |
417 | ||
418 | if (list->IsEmpty()) | |
419 | return 1; | |
420 | ||
421 | TIterator* iter = list->MakeIterator(); | |
422 | TObject* obj; | |
423 | ||
424 | // collections of all histograms | |
eb9356d5 | 425 | TList collections[3*kESDHists+kMCHists*4+kCorrHists*2]; |
0a173978 | 426 | |
427 | Int_t count = 0; | |
428 | while ((obj = iter->Next())) { | |
429 | ||
430 | AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj); | |
431 | if (entry == 0) | |
432 | continue; | |
433 | ||
434 | for (Int_t i = 0; i < kESDHists; ++i) | |
eb9356d5 | 435 | { |
0a173978 | 436 | collections[i].Add(entry->fMultiplicityESD[i]); |
eb9356d5 | 437 | collections[kESDHists+i].Add(entry->fTriggeredEvents[i]); |
438 | collections[kESDHists*2+i].Add(entry->fNoVertexEvents[i]); | |
439 | } | |
0a173978 | 440 | |
441 | for (Int_t i = 0; i < kMCHists; ++i) | |
cfc19dd5 | 442 | { |
eb9356d5 | 443 | collections[3*kESDHists+i].Add(entry->fMultiplicityVtx[i]); |
444 | collections[3*kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]); | |
445 | collections[3*kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]); | |
446 | collections[3*kESDHists+kMCHists*3+i].Add(entry->fMultiplicityNSD[i]); | |
cfc19dd5 | 447 | } |
0a173978 | 448 | |
449 | for (Int_t i = 0; i < kCorrHists; ++i) | |
eb9356d5 | 450 | collections[3*kESDHists+kMCHists*4+i].Add(entry->fCorrelation[i]); |
0a173978 | 451 | |
452 | for (Int_t i = 0; i < kCorrHists; ++i) | |
eb9356d5 | 453 | collections[3*kESDHists+kMCHists*4+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]); |
0a173978 | 454 | |
455 | count++; | |
456 | } | |
457 | ||
458 | for (Int_t i = 0; i < kESDHists; ++i) | |
eb9356d5 | 459 | { |
0a173978 | 460 | fMultiplicityESD[i]->Merge(&collections[i]); |
eb9356d5 | 461 | fTriggeredEvents[i]->Merge(&collections[kESDHists+i]); |
462 | fNoVertexEvents[i]->Merge(&collections[2*kESDHists+i]); | |
463 | } | |
464 | ||
0a173978 | 465 | for (Int_t i = 0; i < kMCHists; ++i) |
cfc19dd5 | 466 | { |
eb9356d5 | 467 | fMultiplicityVtx[i]->Merge(&collections[3*kESDHists+i]); |
468 | fMultiplicityMB[i]->Merge(&collections[3*kESDHists+kMCHists+i]); | |
469 | fMultiplicityINEL[i]->Merge(&collections[3*kESDHists+kMCHists*2+i]); | |
470 | fMultiplicityNSD[i]->Merge(&collections[3*kESDHists+kMCHists*3+i]); | |
cfc19dd5 | 471 | } |
0a173978 | 472 | |
473 | for (Int_t i = 0; i < kCorrHists; ++i) | |
eb9356d5 | 474 | fCorrelation[i]->Merge(&collections[3*kESDHists+kMCHists*4+i]); |
0a173978 | 475 | |
476 | for (Int_t i = 0; i < kCorrHists; ++i) | |
eb9356d5 | 477 | fMultiplicityESDCorrected[i]->Merge(&collections[3*kESDHists+kMCHists*4+kCorrHists+i]); |
0a173978 | 478 | |
479 | delete iter; | |
480 | ||
481 | return count+1; | |
482 | } | |
483 | ||
484 | //____________________________________________________________________ | |
485 | Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir) | |
486 | { | |
487 | // | |
488 | // loads the histograms from a file | |
489 | // if dir is empty a directory with the name of this object is taken (like in SaveHistogram) | |
490 | // | |
491 | ||
492 | if (!dir) | |
493 | dir = GetName(); | |
494 | ||
495 | if (!gDirectory->cd(dir)) | |
496 | return kFALSE; | |
497 | ||
144ff489 | 498 | // store old hists to delete them later |
499 | TList oldObjects; | |
500 | oldObjects.SetOwner(1); | |
501 | for (Int_t i = 0; i < kESDHists; ++i) | |
eb9356d5 | 502 | { |
144ff489 | 503 | if (fMultiplicityESD[i]) |
504 | oldObjects.Add(fMultiplicityESD[i]); | |
eb9356d5 | 505 | if (fTriggeredEvents[i]) |
506 | oldObjects.Add(fTriggeredEvents[i]); | |
507 | if (fNoVertexEvents[i]) | |
508 | oldObjects.Add(fNoVertexEvents[i]); | |
509 | } | |
510 | ||
144ff489 | 511 | for (Int_t i = 0; i < kMCHists; ++i) |
512 | { | |
513 | if (fMultiplicityVtx[i]) | |
514 | oldObjects.Add(fMultiplicityVtx[i]); | |
515 | if (fMultiplicityMB[i]) | |
516 | oldObjects.Add(fMultiplicityMB[i]); | |
517 | if (fMultiplicityINEL[i]) | |
518 | oldObjects.Add(fMultiplicityINEL[i]); | |
69b09e3b | 519 | if (fMultiplicityNSD[i]) |
520 | oldObjects.Add(fMultiplicityNSD[i]); | |
144ff489 | 521 | } |
522 | ||
523 | for (Int_t i = 0; i < kCorrHists; ++i) | |
524 | if (fCorrelation[i]) | |
525 | oldObjects.Add(fCorrelation[i]); | |
526 | ||
527 | // load histograms | |
0a173978 | 528 | |
529 | Bool_t success = kTRUE; | |
530 | ||
531 | for (Int_t i = 0; i < kESDHists; ++i) | |
532 | { | |
533 | fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName())); | |
eb9356d5 | 534 | fTriggeredEvents[i] = dynamic_cast<TH1F*> (gDirectory->Get(fTriggeredEvents[i]->GetName())); |
535 | fNoVertexEvents[i] = dynamic_cast<TH1F*> (gDirectory->Get(fNoVertexEvents[i]->GetName())); | |
536 | if (!fMultiplicityESD[i] || !fTriggeredEvents[i] || !fNoVertexEvents[i]) | |
0a173978 | 537 | success = kFALSE; |
538 | } | |
539 | ||
540 | for (Int_t i = 0; i < kMCHists; ++i) | |
541 | { | |
cfc19dd5 | 542 | fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName())); |
543 | fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName())); | |
544 | fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName())); | |
69b09e3b | 545 | fMultiplicityNSD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityNSD[i]->GetName())); |
cfc19dd5 | 546 | if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i]) |
0a173978 | 547 | success = kFALSE; |
548 | } | |
549 | ||
550 | for (Int_t i = 0; i < kCorrHists; ++i) | |
551 | { | |
552 | fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName())); | |
553 | if (!fCorrelation[i]) | |
554 | success = kFALSE; | |
555 | fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName())); | |
556 | if (!fMultiplicityESDCorrected[i]) | |
557 | success = kFALSE; | |
558 | } | |
559 | ||
560 | gDirectory->cd(".."); | |
561 | ||
144ff489 | 562 | // delete old hists |
563 | oldObjects.Delete(); | |
564 | ||
0a173978 | 565 | return success; |
566 | } | |
567 | ||
568 | //____________________________________________________________________ | |
144ff489 | 569 | void AliMultiplicityCorrection::SaveHistograms(const char* dir) |
0a173978 | 570 | { |
571 | // | |
572 | // saves the histograms | |
573 | // | |
574 | ||
144ff489 | 575 | if (!dir) |
576 | dir = GetName(); | |
577 | ||
578 | gDirectory->mkdir(dir); | |
579 | gDirectory->cd(dir); | |
0a173978 | 580 | |
581 | for (Int_t i = 0; i < kESDHists; ++i) | |
eb9356d5 | 582 | { |
0a173978 | 583 | if (fMultiplicityESD[i]) |
69b09e3b | 584 | { |
0a173978 | 585 | fMultiplicityESD[i]->Write(); |
69b09e3b | 586 | fMultiplicityESD[i]->ProjectionY(Form("%s_px", fMultiplicityESD[i]->GetName()), 1, fMultiplicityESD[i]->GetNbinsX())->Write(); |
587 | } | |
eb9356d5 | 588 | if (fTriggeredEvents[i]) |
589 | fTriggeredEvents[i]->Write(); | |
590 | if (fNoVertexEvents[i]) | |
591 | fNoVertexEvents[i]->Write(); | |
592 | } | |
0a173978 | 593 | |
594 | for (Int_t i = 0; i < kMCHists; ++i) | |
cfc19dd5 | 595 | { |
596 | if (fMultiplicityVtx[i]) | |
69b09e3b | 597 | { |
cfc19dd5 | 598 | fMultiplicityVtx[i]->Write(); |
69b09e3b | 599 | fMultiplicityVtx[i]->ProjectionY(Form("%s_px", fMultiplicityVtx[i]->GetName()), 1, fMultiplicityVtx[i]->GetNbinsX())->Write(); |
600 | } | |
cfc19dd5 | 601 | if (fMultiplicityMB[i]) |
69b09e3b | 602 | { |
cfc19dd5 | 603 | fMultiplicityMB[i]->Write(); |
69b09e3b | 604 | fMultiplicityMB[i]->ProjectionY(Form("%s_px", fMultiplicityMB[i]->GetName()), 1, fMultiplicityMB[i]->GetNbinsX())->Write(); |
605 | } | |
cfc19dd5 | 606 | if (fMultiplicityINEL[i]) |
69b09e3b | 607 | { |
cfc19dd5 | 608 | fMultiplicityINEL[i]->Write(); |
69b09e3b | 609 | fMultiplicityINEL[i]->ProjectionY(Form("%s_px", fMultiplicityINEL[i]->GetName()), 1, fMultiplicityINEL[i]->GetNbinsX())->Write(); |
610 | } | |
611 | if (fMultiplicityNSD[i]) | |
612 | { | |
613 | fMultiplicityNSD[i]->Write(); | |
614 | fMultiplicityNSD[i]->ProjectionY(Form("%s_px", fMultiplicityNSD[i]->GetName()), 1, fMultiplicityNSD[i]->GetNbinsX())->Write(); | |
615 | } | |
cfc19dd5 | 616 | } |
0a173978 | 617 | |
618 | for (Int_t i = 0; i < kCorrHists; ++i) | |
619 | { | |
620 | if (fCorrelation[i]) | |
621 | fCorrelation[i]->Write(); | |
622 | if (fMultiplicityESDCorrected[i]) | |
623 | fMultiplicityESDCorrected[i]->Write(); | |
624 | } | |
625 | ||
626 | gDirectory->cd(".."); | |
627 | } | |
628 | ||
629 | //____________________________________________________________________ | |
b3b260d1 | 630 | void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, AliPWG0Helper::MCProcessType processType, Int_t generated05, Int_t generated10, Int_t generated14, Int_t generatedAll) |
0a173978 | 631 | { |
632 | // | |
633 | // Fills an event from MC | |
634 | // | |
635 | ||
cfc19dd5 | 636 | if (triggered) |
637 | { | |
638 | fMultiplicityMB[0]->Fill(vtx, generated05); | |
639 | fMultiplicityMB[1]->Fill(vtx, generated10); | |
b3b260d1 | 640 | fMultiplicityMB[2]->Fill(vtx, generated14); |
641 | fMultiplicityMB[3]->Fill(vtx, generatedAll); | |
cfc19dd5 | 642 | |
643 | if (vertex) | |
644 | { | |
645 | fMultiplicityVtx[0]->Fill(vtx, generated05); | |
646 | fMultiplicityVtx[1]->Fill(vtx, generated10); | |
b3b260d1 | 647 | fMultiplicityVtx[2]->Fill(vtx, generated14); |
648 | fMultiplicityVtx[3]->Fill(vtx, generatedAll); | |
cfc19dd5 | 649 | } |
650 | } | |
651 | ||
652 | fMultiplicityINEL[0]->Fill(vtx, generated05); | |
653 | fMultiplicityINEL[1]->Fill(vtx, generated10); | |
b3b260d1 | 654 | fMultiplicityINEL[2]->Fill(vtx, generated14); |
655 | fMultiplicityINEL[3]->Fill(vtx, generatedAll); | |
69b09e3b | 656 | |
657 | if (processType != AliPWG0Helper::kSD) | |
658 | { | |
659 | fMultiplicityNSD[0]->Fill(vtx, generated05); | |
660 | fMultiplicityNSD[1]->Fill(vtx, generated10); | |
b3b260d1 | 661 | fMultiplicityNSD[2]->Fill(vtx, generated14); |
662 | fMultiplicityNSD[3]->Fill(vtx, generatedAll); | |
69b09e3b | 663 | } |
0a173978 | 664 | } |
665 | ||
666 | //____________________________________________________________________ | |
b3b260d1 | 667 | void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured14) |
0a173978 | 668 | { |
669 | // | |
670 | // Fills an event from ESD | |
671 | // | |
672 | ||
673 | fMultiplicityESD[0]->Fill(vtx, measured05); | |
674 | fMultiplicityESD[1]->Fill(vtx, measured10); | |
b3b260d1 | 675 | fMultiplicityESD[2]->Fill(vtx, measured14); |
0a173978 | 676 | } |
677 | ||
eb9356d5 | 678 | //____________________________________________________________________ |
679 | void AliMultiplicityCorrection::FillTriggeredEvent(Int_t measured05, Int_t measured10, Int_t measured14) | |
680 | { | |
681 | // | |
682 | // fills raw distribution of triggered events | |
683 | // | |
684 | ||
685 | fTriggeredEvents[0]->Fill(measured05); | |
686 | fTriggeredEvents[1]->Fill(measured10); | |
687 | fTriggeredEvents[2]->Fill(measured14); | |
688 | } | |
689 | ||
690 | //____________________________________________________________________ | |
691 | void AliMultiplicityCorrection::FillNoVertexEvent(Float_t vtx, Bool_t vertexReconstructed, Int_t generated05, Int_t generated10, Int_t generated14, Int_t measured05, Int_t measured10, Int_t measured14) | |
692 | { | |
693 | // | |
694 | // fills raw distribution of triggered events | |
695 | // | |
696 | ||
697 | if (vtx > fgVtxRangeBegin[0] && vtx < fgVtxRangeEnd[0] && (!vertexReconstructed || measured05 == 0)) | |
698 | fNoVertexEvents[0]->Fill(generated05); | |
699 | ||
700 | if (vtx > fgVtxRangeBegin[1] && vtx < fgVtxRangeEnd[1] && (!vertexReconstructed || measured10 == 0)) | |
701 | fNoVertexEvents[1]->Fill(generated10); | |
702 | ||
703 | if (vtx > fgVtxRangeBegin[2] && vtx < fgVtxRangeEnd[2] && (!vertexReconstructed || measured14 == 0)) | |
704 | fNoVertexEvents[2]->Fill(generated14); | |
705 | } | |
706 | ||
0a173978 | 707 | //____________________________________________________________________ |
b3b260d1 | 708 | void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated14, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured14) |
0a173978 | 709 | { |
710 | // | |
711 | // Fills an event into the correlation map with the information from MC and ESD | |
712 | // | |
713 | ||
714 | fCorrelation[0]->Fill(vtx, generated05, measured05); | |
715 | fCorrelation[1]->Fill(vtx, generated10, measured10); | |
b3b260d1 | 716 | fCorrelation[2]->Fill(vtx, generated14, measured14); |
0a173978 | 717 | |
b3b260d1 | 718 | fCorrelation[3]->Fill(vtx, generatedAll, measured05); |
719 | fCorrelation[4]->Fill(vtx, generatedAll, measured10); | |
720 | fCorrelation[5]->Fill(vtx, generatedAll, measured14); | |
0a173978 | 721 | } |
722 | ||
723 | //____________________________________________________________________ | |
19442b86 | 724 | void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType) |
0a173978 | 725 | { |
726 | // | |
9ca6be09 | 727 | // fills fCurrentESD, fCurrentCorrelation |
728 | // resets fMultiplicityESDCorrected | |
0a173978 | 729 | // |
730 | ||
731 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
0b4bfd98 | 732 | |
9ca6be09 | 733 | fMultiplicityESDCorrected[correlationID]->Reset(); |
0b4bfd98 | 734 | fMultiplicityESDCorrected[correlationID]->Sumw2(); |
0a173978 | 735 | |
039e265e | 736 | // project without under/overflow bins |
eb9356d5 | 737 | Int_t begin = 1; |
738 | Int_t end = fMultiplicityESD[inputRange]->GetXaxis()->GetNbins(); | |
739 | if (fVtxEnd > fVtxBegin) | |
740 | { | |
741 | begin = fVtxBegin; | |
742 | end = fVtxEnd; | |
743 | } | |
744 | fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD", begin, end); | |
9ca6be09 | 745 | fCurrentESD->Sumw2(); |
0a173978 | 746 | |
747 | // empty under/overflow bins in x, otherwise Project3D takes them into account | |
748 | TH3* hist = fCorrelation[correlationID]; | |
039e265e | 749 | for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y) |
0a173978 | 750 | { |
039e265e | 751 | for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z) |
0a173978 | 752 | { |
753 | hist->SetBinContent(0, y, z, 0); | |
754 | hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); | |
755 | } | |
756 | } | |
757 | ||
eb9356d5 | 758 | if (fVtxEnd > fVtxBegin) |
759 | hist->GetXaxis()->SetRange(fVtxBegin, fVtxEnd); | |
760 | ||
19442b86 | 761 | fCurrentCorrelation = (TH2*) hist->Project3D("zy"); |
9ca6be09 | 762 | fCurrentCorrelation->Sumw2(); |
69b09e3b | 763 | |
44466df2 | 764 | Printf("AliMultiplicityCorrection::SetupCurrentHists: Statistics information: %.f entries in correlation map; %.f entries in measured spectrum", fCurrentCorrelation->Integral(), fCurrentESD->Integral()); |
765 | ||
0b4bfd98 | 766 | #if 0 // does not help |
767 | // null bins with one entry | |
768 | Int_t nNulledBins = 0; | |
769 | for (Int_t x=1; x<=fCurrentCorrelation->GetXaxis()->GetNbins(); ++x) | |
770 | for (Int_t y=1; y<=fCurrentCorrelation->GetYaxis()->GetNbins(); ++y) | |
771 | { | |
772 | if (fCurrentCorrelation->GetBinContent(x, y) == 1) | |
773 | { | |
774 | fCurrentCorrelation->SetBinContent(x, y, 0); | |
775 | fCurrentCorrelation->SetBinError(x, y, 0); | |
447c325d | 776 | |
0b4bfd98 | 777 | ++nNulledBins; |
778 | } | |
779 | } | |
780 | Printf("Nulled %d bins", nNulledBins); | |
781 | #endif | |
782 | ||
783 | fCurrentEfficiency = GetEfficiency(inputRange, eventType); | |
784 | //fCurrentEfficiency->Rebin(2); | |
785 | //fCurrentEfficiency->Scale(0.5); | |
786 | } | |
787 | ||
788 | //____________________________________________________________________ | |
789 | TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventType) | |
790 | { | |
791 | // | |
792 | // calculates efficiency for given event type | |
793 | // | |
69b09e3b | 794 | |
eb9356d5 | 795 | TString name1; |
796 | name1.Form("divisor%d", inputRange); | |
797 | ||
798 | TString name2; | |
799 | name2.Form("CurrentEfficiency%d", inputRange); | |
800 | ||
0b4bfd98 | 801 | TH1* divisor = 0; |
cfc19dd5 | 802 | switch (eventType) |
803 | { | |
69b09e3b | 804 | case kTrVtx : break; |
eb9356d5 | 805 | case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY(name1, 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); break; |
806 | case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY(name1, 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e"); break; | |
807 | case kNSD: divisor = fMultiplicityNSD[inputRange]->ProjectionY(name1, 1, fMultiplicityNSD[inputRange]->GetNbinsX(), "e"); break; | |
69b09e3b | 808 | } |
eb9356d5 | 809 | TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY(name2, 1, fMultiplicityVtx[inputRange]->GetNbinsX(), "e"); |
69b09e3b | 810 | |
811 | if (eventType == kTrVtx) | |
812 | { | |
813 | for (Int_t i=0; i<= eff->GetNbinsX()+1; i++) | |
814 | eff->SetBinContent(i, 1); | |
cfc19dd5 | 815 | } |
69b09e3b | 816 | else |
817 | eff->Divide(eff, divisor, 1, 1, "B"); | |
818 | ||
819 | return eff; | |
820 | } | |
821 | ||
822 | //____________________________________________________________________ | |
eb9356d5 | 823 | TH1* AliMultiplicityCorrection::GetTriggerEfficiency(Int_t inputRange, EventType eventType) |
69b09e3b | 824 | { |
825 | // | |
826 | // calculates efficiency for given event type | |
827 | // | |
828 | ||
eb9356d5 | 829 | TString name1; |
830 | name1.Form("divisor%d", inputRange); | |
831 | ||
832 | TString name2; | |
833 | name2.Form("CurrentEfficiency%d", inputRange); | |
834 | ||
835 | TH1* divisor = 0; | |
836 | switch (eventType) | |
837 | { | |
838 | case kTrVtx : AliFatal("Not supported!"); break; | |
839 | case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY (name1, 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); break; | |
840 | case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY(name1, 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e"); break; | |
841 | case kNSD: divisor = fMultiplicityNSD[inputRange]->ProjectionY (name1, 1, fMultiplicityNSD[inputRange]->GetNbinsX(), "e"); break; | |
842 | } | |
843 | TH1* eff = fMultiplicityMB[inputRange]->ProjectionY(name2, 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); | |
844 | ||
0b4bfd98 | 845 | eff->Divide(eff, divisor, 1, 1, "B"); |
846 | return eff; | |
9ca6be09 | 847 | } |
848 | ||
849 | //____________________________________________________________________ | |
eb9356d5 | 850 | Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Int_t zeroBinEvents, Bool_t check, TH1* initialConditions, Bool_t errorAsBias) |
6d81c2de | 851 | { |
852 | // | |
19442b86 | 853 | // correct spectrum using minuit chi2 method |
6d81c2de | 854 | // |
19442b86 | 855 | // for description of parameters, see AliUnfolding::Unfold |
6d81c2de | 856 | // |
447c325d | 857 | |
19442b86 | 858 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); |
859 | ||
eb9356d5 | 860 | //AliUnfolding::SetCreateOverflowBin(5); |
19442b86 | 861 | AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization); |
eb9356d5 | 862 | AliUnfolding::SetMinimumInitialValue(kTRUE, 0.1); |
863 | ||
864 | // use here only vtx efficiency (to MB sample) which is always needed if we use the 0 bin | |
865 | SetupCurrentHists(inputRange, fullPhaseSpace, (eventType == kTrVtx) ? kTrVtx : kMB); | |
866 | ||
867 | // TODO set errors on measured with 0.5 * TMath::ChisquareQuantile(0.1, 20000) | |
868 | // see PDG: Statistics / Poission or binomial data / Eq. 32.49a/b in 2004 edition | |
869 | ||
870 | Calculate0Bin(inputRange, eventType, zeroBinEvents); | |
871 | ||
872 | Int_t resultCode = -1; | |
873 | if (errorAsBias == kFALSE) | |
874 | { | |
875 | resultCode = AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check); | |
876 | } | |
877 | else | |
878 | { | |
879 | resultCode = AliUnfolding::UnfoldGetBias(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]); | |
880 | } | |
881 | ||
882 | // HACK store new vertex reco efficiency for bin 0, changing number of events with trigger and vertex in MC map | |
883 | if (zeroBinEvents > 0) | |
884 | { | |
885 | Printf("WARNING: Stored vertex reco efficiency from unfolding for bin 0."); | |
886 | fMultiplicityVtx[inputRange]->SetBinContent(1, 1, fMultiplicityMB[inputRange]->GetBinContent(1, 1) * fCurrentEfficiency->GetBinContent(1)); | |
887 | } | |
888 | ||
889 | // correct for the trigger bias if requested | |
890 | if (eventType > kMB) | |
891 | { | |
892 | Printf("Applying trigger efficiency"); | |
893 | TH1* eff = GetTriggerEfficiency(inputRange, eventType); | |
894 | for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); i++) | |
895 | { | |
896 | fMultiplicityESDCorrected[correlationID]->SetBinContent(i, fMultiplicityESDCorrected[correlationID]->GetBinContent(i) / eff->GetBinContent(i)); | |
897 | fMultiplicityESDCorrected[correlationID]->SetBinError(i, fMultiplicityESDCorrected[correlationID]->GetBinError(i) / eff->GetBinContent(i)); | |
898 | } | |
899 | } | |
900 | ||
901 | return resultCode; | |
902 | } | |
903 | ||
904 | //____________________________________________________________________ | |
905 | void AliMultiplicityCorrection::Calculate0Bin(Int_t inputRange, EventType eventType, Int_t zeroBinEvents) | |
906 | { | |
907 | // fills the 0 bin | |
19442b86 | 908 | |
eb9356d5 | 909 | if (eventType == kTrVtx) |
910 | return; | |
911 | ||
912 | Double_t fractionEventsInVertexRange = fMultiplicityESD[inputRange]->Integral(1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins()) / fMultiplicityESD[inputRange]->Integral(0, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins() + 1); | |
913 | ||
914 | // difference of fraction that is inside the considered range between triggered events and events with vertex | |
915 | // Extension to NSD not needed, INEL and NSD vertex distributions are nature-given and unbiased! | |
916 | Double_t differenceVtxDist = (fMultiplicityINEL[inputRange]->Integral(1, fMultiplicityINEL[inputRange]->GetXaxis()->GetNbins()) / fMultiplicityINEL[inputRange]->Integral(0, fMultiplicityINEL[inputRange]->GetXaxis()->GetNbins() + 1)) / (fMultiplicityVtx[inputRange]->Integral(1, fMultiplicityVtx[inputRange]->GetXaxis()->GetNbins()) / fMultiplicityVtx[inputRange]->Integral(0, fMultiplicityVtx[inputRange]->GetXaxis()->GetNbins() + 1)); | |
917 | ||
918 | Printf("Enabling 0 bin estimate for triggered events without vertex."); | |
919 | Printf(" Events in 0 bin: %d", zeroBinEvents); | |
920 | Printf(" Fraction in range: %.1f%%", fractionEventsInVertexRange * 100); | |
921 | Printf(" Difference Vtx Dist: %f", differenceVtxDist); | |
922 | ||
923 | AliUnfolding::SetNotFoundEvents(zeroBinEvents * fractionEventsInVertexRange * differenceVtxDist); | |
924 | } | |
925 | ||
926 | //____________________________________________________________________ | |
927 | void AliMultiplicityCorrection::FixTriggerEfficiencies(Int_t start) | |
928 | { | |
929 | // | |
930 | // sets trigger and vertex efficiencies to 1 for large multiplicities where no event was simulated | |
931 | // | |
932 | ||
933 | for (Int_t etaRange = 0; etaRange < kMCHists; etaRange++) | |
0a173978 | 934 | { |
eb9356d5 | 935 | for (Int_t i = 1; i <= fMultiplicityINEL[etaRange]->GetNbinsY(); i++) |
dd701109 | 936 | { |
eb9356d5 | 937 | if (fMultiplicityINEL[etaRange]->GetYaxis()->GetBinCenter(i) < start) |
938 | continue; | |
939 | ||
940 | if (fMultiplicityINEL[etaRange]->GetBinContent(1, i) > 0) | |
941 | continue; | |
942 | ||
943 | fMultiplicityINEL[etaRange]->SetBinContent(1, i, 1); | |
944 | fMultiplicityNSD[etaRange]->SetBinContent(1, i, 1); | |
945 | fMultiplicityMB[etaRange]->SetBinContent(1, i, 1); | |
946 | fMultiplicityVtx[etaRange]->SetBinContent(1, i, 1); | |
dd701109 | 947 | } |
948 | } | |
eb9356d5 | 949 | } |
0f67a57c | 950 | |
eb9356d5 | 951 | //____________________________________________________________________ |
952 | Float_t AliMultiplicityCorrection::GetFraction0Generated(Int_t inputRange) | |
953 | { | |
954 | // | |
955 | // returns the fraction of events that have 0 generated particles in the given range among all events without vertex OR 0 tracklets | |
956 | // | |
957 | ||
958 | TH1* multMB = GetNoVertexEvents(inputRange); | |
959 | return multMB->GetBinContent(1) / multMB->Integral(); | |
dd701109 | 960 | } |
961 | ||
6d81c2de | 962 | //____________________________________________________________________ |
19442b86 | 963 | Int_t AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType) |
6d81c2de | 964 | { |
965 | // | |
19442b86 | 966 | // correct spectrum using minuit chi2 method with a NBD function |
dd701109 | 967 | // |
19442b86 | 968 | // for description of parameters, see AliUnfolding::Unfold |
dd701109 | 969 | // |
970 | ||
971 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
19442b86 | 972 | |
973 | AliUnfolding::SetUnfoldingMethod(AliUnfolding::kFunction); | |
974 | SetupCurrentHists(inputRange, fullPhaseSpace, eventType); | |
975 | ||
976 | TF1* func = new TF1("nbd", "[0] * TMath::Gamma([2]+x) / TMath::Gamma([2]) / TMath::Gamma(x+1) * pow([1] / ([1]+[2]), x) * pow(1.0 + [1]/[2], -[2])"); | |
977 | func->SetParNames("scaling", "averagen", "k"); | |
978 | func->SetParLimits(0, 0, 1000); | |
979 | func->SetParLimits(1, 1, 50); | |
980 | func->SetParLimits(2, 1, 10); | |
981 | func->SetParameters(1, 10, 2); | |
982 | AliUnfolding::SetFunction(func); | |
983 | ||
984 | return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, 0, fMultiplicityESDCorrected[correlationID]); | |
0a173978 | 985 | } |
986 | ||
0a173978 | 987 | //____________________________________________________________________ |
988 | void AliMultiplicityCorrection::DrawHistograms() | |
989 | { | |
6d81c2de | 990 | // |
991 | // draws the histograms of this class | |
992 | // | |
993 | ||
0a173978 | 994 | printf("ESD:\n"); |
995 | ||
996 | TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600); | |
997 | canvas1->Divide(3, 2); | |
998 | for (Int_t i = 0; i < kESDHists; ++i) | |
999 | { | |
1000 | canvas1->cd(i+1); | |
1001 | fMultiplicityESD[i]->DrawCopy("COLZ"); | |
1002 | printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean()); | |
1003 | } | |
1004 | ||
cfc19dd5 | 1005 | printf("Vtx:\n"); |
0a173978 | 1006 | |
1007 | TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600); | |
1008 | canvas2->Divide(3, 2); | |
1009 | for (Int_t i = 0; i < kMCHists; ++i) | |
1010 | { | |
1011 | canvas2->cd(i+1); | |
cfc19dd5 | 1012 | fMultiplicityVtx[i]->DrawCopy("COLZ"); |
1013 | printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean()); | |
0a173978 | 1014 | } |
1015 | ||
1016 | TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900); | |
1017 | canvas3->Divide(3, 3); | |
1018 | for (Int_t i = 0; i < kCorrHists; ++i) | |
1019 | { | |
1020 | canvas3->cd(i+1); | |
335a5e1b | 1021 | TH3* hist = static_cast<TH3*> (fCorrelation[i]->Clone()); |
b477f4f2 | 1022 | for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y) |
1023 | { | |
1024 | for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z) | |
1025 | { | |
1026 | hist->SetBinContent(0, y, z, 0); | |
1027 | hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); | |
1028 | } | |
1029 | } | |
1030 | TH1* proj = hist->Project3D("zy"); | |
0a173978 | 1031 | proj->DrawCopy("COLZ"); |
1032 | } | |
1033 | } | |
1034 | ||
1035 | //____________________________________________________________________ | |
7dcf959e | 1036 | void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t /*normalizeESD*/, TH1* mcHist, Bool_t simple, EventType eventType) |
0a173978 | 1037 | { |
039e265e | 1038 | // draw comparison plots |
eb9356d5 | 1039 | |
cfc19dd5 | 1040 | Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); |
eb9356d5 | 1041 | |
9ca6be09 | 1042 | TString tmpStr; |
cfc19dd5 | 1043 | tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId); |
0a173978 | 1044 | |
447c325d | 1045 | if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0) |
1046 | { | |
1047 | printf("ERROR. Unfolded histogram is empty\n"); | |
1048 | return; | |
1049 | } | |
eb9356d5 | 1050 | |
1051 | Int_t begin = 1; | |
1052 | Int_t end = fMultiplicityESD[inputRange]->GetXaxis()->GetNbins(); | |
1053 | if (fVtxEnd > fVtxBegin) | |
0b4bfd98 | 1054 | { |
eb9356d5 | 1055 | begin = fVtxBegin; |
1056 | end = fVtxEnd; | |
0b4bfd98 | 1057 | } |
eb9356d5 | 1058 | fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY("fCurrentESD", begin, end); |
1059 | fCurrentESD->Sumw2(); | |
b477f4f2 | 1060 | |
eb9356d5 | 1061 | mcHist->Sumw2(); |
1062 | Int_t mcMax = 0; | |
1063 | for (Int_t i=5; i<=mcHist->GetNbinsX(); ++i) | |
447c325d | 1064 | { |
eb9356d5 | 1065 | if (mcHist->GetBinContent(i) > 0) |
ad0385e9 | 1066 | mcMax = (Int_t) mcHist->GetXaxis()->GetBinCenter(i) + 2; |
447c325d | 1067 | } |
eb9356d5 | 1068 | if (mcMax == 0) |
1069 | { | |
1070 | for (Int_t i=5; i<=fMultiplicityESDCorrected[esdCorrId]->GetNbinsX(); ++i) | |
1071 | if (fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i) > 1) | |
ad0385e9 | 1072 | mcMax = (Int_t) fMultiplicityESDCorrected[esdCorrId]->GetXaxis()->GetBinCenter(i) + 2; |
eb9356d5 | 1073 | } |
1074 | Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcMax); | |
1075 | // calculate residual | |
1076 | Float_t tmp; | |
1077 | TH1* convolutedProj = (TH1*) GetConvoluted(esdCorrId, eventType)->Clone("convolutedProj"); | |
1078 | TH1* residual = GetResiduals(esdCorrId, eventType, tmp); | |
447c325d | 1079 | |
0b4bfd98 | 1080 | TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5); |
b477f4f2 | 1081 | |
447c325d | 1082 | Float_t chi2 = 0; |
eb9356d5 | 1083 | for (Int_t i=1; i<=TMath::Min(residual->GetNbinsX(), 75); ++i) |
b477f4f2 | 1084 | { |
eb9356d5 | 1085 | Float_t value = residual->GetBinContent(i); |
1086 | // TODO has to get a parameter (used in Chi2Function, GetResiduals, and here) | |
1087 | if (i > 1) | |
1088 | chi2 += value * value; | |
1089 | Printf("%d --> %f (%f)", i, value * value, chi2); | |
1090 | residualHist->Fill(value); | |
447c325d | 1091 | convolutedProj->SetBinError(i, 0); |
447c325d | 1092 | } |
0f67a57c | 1093 | fLastChi2Residuals = chi2; |
dd701109 | 1094 | |
eb9356d5 | 1095 | //new TCanvas; residualHist->DrawCopy(); |
dd701109 | 1096 | |
eb9356d5 | 1097 | printf("Difference (Residuals) is %f\n", fLastChi2Residuals); |
dd701109 | 1098 | |
447c325d | 1099 | TCanvas* canvas1 = 0; |
1100 | if (simple) | |
1101 | { | |
69b09e3b | 1102 | canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 600); |
447c325d | 1103 | canvas1->Divide(2, 1); |
1104 | } | |
1105 | else | |
1106 | { | |
1107 | canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200); | |
1108 | canvas1->Divide(2, 3); | |
1109 | } | |
0a173978 | 1110 | |
1111 | canvas1->cd(1); | |
69b09e3b | 1112 | canvas1->cd(1)->SetGridx(); |
1113 | canvas1->cd(1)->SetGridy(); | |
1114 | canvas1->cd(1)->SetTopMargin(0.05); | |
1115 | canvas1->cd(1)->SetRightMargin(0.05); | |
1116 | canvas1->cd(1)->SetLeftMargin(0.12); | |
1117 | canvas1->cd(1)->SetBottomMargin(0.12); | |
cfc19dd5 | 1118 | TH1* proj = (TH1*) mcHist->Clone("proj"); |
eb9356d5 | 1119 | if (proj->GetEntries() > 0) |
1120 | AliPWG0Helper::NormalizeToBinWidth(proj); | |
9ca6be09 | 1121 | |
eb9356d5 | 1122 | proj->GetXaxis()->SetRangeUser(0, mcMax); |
69b09e3b | 1123 | proj->GetYaxis()->SetTitleOffset(1.4); |
69b09e3b | 1124 | proj->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5)); |
447c325d | 1125 | proj->SetStats(kFALSE); |
0a173978 | 1126 | |
eb9356d5 | 1127 | fMultiplicityESDCorrected[esdCorrId]->GetXaxis()->SetRangeUser(0, mcMax); |
1128 | fMultiplicityESDCorrected[esdCorrId]->GetYaxis()->SetRangeUser(0.1, fMultiplicityESDCorrected[esdCorrId]->GetMaximum() * 1.5); | |
1129 | fMultiplicityESDCorrected[esdCorrId]->GetYaxis()->SetTitleOffset(1.4); | |
1130 | fMultiplicityESDCorrected[esdCorrId]->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5)); | |
1131 | fMultiplicityESDCorrected[esdCorrId]->SetStats(kFALSE); | |
1132 | ||
cfc19dd5 | 1133 | fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2); |
69b09e3b | 1134 | fMultiplicityESDCorrected[esdCorrId]->SetMarkerColor(2); |
0f67a57c | 1135 | |
eb9356d5 | 1136 | TH1* esdCorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("esdCorrected"); |
1137 | AliPWG0Helper::NormalizeToBinWidth(esdCorrected); | |
1138 | ||
0f67a57c | 1139 | if (proj->GetEntries() > 0) { |
1140 | proj->DrawCopy("HIST"); | |
eb9356d5 | 1141 | esdCorrected->DrawCopy("SAME HIST E"); |
0f67a57c | 1142 | } |
1143 | else | |
eb9356d5 | 1144 | esdCorrected->DrawCopy("HIST E"); |
1145 | ||
0f67a57c | 1146 | gPad->SetLogy(); |
cfc19dd5 | 1147 | |
447c325d | 1148 | TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93); |
69b09e3b | 1149 | legend->AddEntry(proj, "True distribution"); |
1150 | legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "Unfolded distribution"); | |
447c325d | 1151 | legend->SetFillColor(0); |
69b09e3b | 1152 | legend->SetTextSize(0.04); |
447c325d | 1153 | legend->Draw(); |
447c325d | 1154 | |
1155 | canvas1->cd(2); | |
69b09e3b | 1156 | canvas1->cd(2)->SetGridx(); |
1157 | canvas1->cd(2)->SetGridy(); | |
1158 | canvas1->cd(2)->SetTopMargin(0.05); | |
1159 | canvas1->cd(2)->SetRightMargin(0.05); | |
1160 | canvas1->cd(2)->SetLeftMargin(0.12); | |
1161 | canvas1->cd(2)->SetBottomMargin(0.12); | |
447c325d | 1162 | |
1163 | gPad->SetLogy(); | |
eb9356d5 | 1164 | fCurrentESD->GetXaxis()->SetRangeUser(0, mcMax); |
69b09e3b | 1165 | fCurrentESD->SetTitle(Form(";Measured multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5)); |
447c325d | 1166 | fCurrentESD->SetStats(kFALSE); |
69b09e3b | 1167 | fCurrentESD->GetYaxis()->SetTitleOffset(1.4); |
447c325d | 1168 | fCurrentESD->DrawCopy("HIST E"); |
1169 | ||
1170 | convolutedProj->SetLineColor(2); | |
69b09e3b | 1171 | convolutedProj->SetMarkerColor(2); |
1172 | convolutedProj->SetMarkerStyle(5); | |
69b09e3b | 1173 | convolutedProj->DrawCopy("HIST SAME P"); |
447c325d | 1174 | |
1175 | legend = new TLegend(0.3, 0.8, 0.93, 0.93); | |
69b09e3b | 1176 | legend->AddEntry(fCurrentESD, "Measured distribution"); |
1177 | legend->AddEntry(convolutedProj, "R #otimes unfolded distribution", "P"); | |
447c325d | 1178 | legend->SetFillColor(0); |
69b09e3b | 1179 | legend->SetTextSize(0.04); |
447c325d | 1180 | legend->Draw(); |
1181 | ||
447c325d | 1182 | if (!simple) |
1183 | { | |
447c325d | 1184 | canvas1->cd(4); |
039e265e | 1185 | residual->GetYaxis()->SetRangeUser(-5, 5); |
eb9356d5 | 1186 | residual->GetXaxis()->SetRangeUser(0, mcMax); |
1187 | residual->SetStats(kFALSE); | |
447c325d | 1188 | residual->DrawCopy(); |
1189 | ||
1190 | canvas1->cd(5); | |
447c325d | 1191 | TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio"); |
1192 | ratio->Divide(mcHist); | |
1193 | ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC"); | |
1194 | ratio->GetYaxis()->SetRangeUser(0.5, 1.5); | |
eb9356d5 | 1195 | ratio->GetXaxis()->SetRangeUser(0, mcMax); |
447c325d | 1196 | ratio->SetStats(kFALSE); |
1197 | ratio->Draw("HIST"); | |
1198 | ||
0b4bfd98 | 1199 | // plot (MC - Unfolded) / error (MC) |
1200 | canvas1->cd(3); | |
1201 | ||
335a5e1b | 1202 | TH1* diffMCUnfolded2 = static_cast<TH1*> (proj->Clone("diffMCUnfolded2")); |
eb9356d5 | 1203 | diffMCUnfolded2->Add(esdCorrected, -1); |
0b4bfd98 | 1204 | |
1205 | Int_t ndfQual[kQualityRegions]; | |
1206 | for (Int_t region=0; region<kQualityRegions; ++region) | |
1207 | { | |
1208 | fQuality[region] = 0; | |
1209 | ndfQual[region] = 0; | |
1210 | } | |
1211 | ||
0b4bfd98 | 1212 | Double_t newChi2 = 0; |
6d81c2de | 1213 | Double_t newChi2Limit150 = 0; |
0b4bfd98 | 1214 | Int_t ndf = 0; |
19442b86 | 1215 | for (Int_t i=1; i<=diffMCUnfolded2->GetNbinsX(); ++i) |
0b4bfd98 | 1216 | { |
1217 | Double_t value = 0; | |
1218 | if (proj->GetBinError(i) > 0) | |
1219 | { | |
1220 | value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i); | |
1221 | newChi2 += value * value; | |
eb9356d5 | 1222 | if (i > 1 && i <= mcMax) |
6d81c2de | 1223 | newChi2Limit150 += value * value; |
0b4bfd98 | 1224 | ++ndf; |
1225 | ||
1226 | for (Int_t region=0; region<kQualityRegions; ++region) | |
1227 | if (diffMCUnfolded2->GetXaxis()->GetBinCenter(i) >= fgQualityRegionsB[region] - 0.1 && diffMCUnfolded2->GetXaxis()->GetBinCenter(i) <= fgQualityRegionsE[region] + 0.1) // 0.1 to avoid e.g. 3.9999 < 4 problem | |
1228 | { | |
1229 | fQuality[region] += TMath::Abs(value); | |
1230 | ++ndfQual[region]; | |
1231 | } | |
1232 | } | |
1233 | ||
1234 | diffMCUnfolded2->SetBinContent(i, value); | |
1235 | } | |
1236 | ||
1237 | // normalize region to the number of entries | |
1238 | for (Int_t region=0; region<kQualityRegions; ++region) | |
1239 | { | |
1240 | if (ndfQual[region] > 0) | |
1241 | fQuality[region] /= ndfQual[region]; | |
1242 | Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]); | |
1243 | } | |
1244 | ||
eb9356d5 | 1245 | if (mcMax > 1) |
0b4bfd98 | 1246 | { |
eb9356d5 | 1247 | fLastChi2MC = newChi2Limit150 / (mcMax - 1); |
1248 | Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcMax, newChi2Limit150, mcMax - 1, fLastChi2MC); | |
0b4bfd98 | 1249 | } |
1250 | else | |
1251 | fLastChi2MC = -1; | |
1252 | ||
144ff489 | 1253 | Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, ((ndf > 0) ? newChi2 / ndf : -1)); |
0b4bfd98 | 1254 | |
eb9356d5 | 1255 | diffMCUnfolded2->SetTitle("#chi^{2};true multiplicity;(MC - Unfolded) / e(MC)"); |
039e265e | 1256 | diffMCUnfolded2->GetYaxis()->SetRangeUser(-5, 5); |
eb9356d5 | 1257 | diffMCUnfolded2->GetXaxis()->SetRangeUser(0, mcMax); |
0b4bfd98 | 1258 | diffMCUnfolded2->DrawCopy("HIST"); |
eb9356d5 | 1259 | |
1260 | canvas1->cd(6); | |
1261 | // draw penalty factor | |
1262 | ||
1263 | TH1* penalty = AliUnfolding::GetPenaltyPlot(fMultiplicityESDCorrected[esdCorrId]); | |
1264 | penalty->SetStats(0); | |
1265 | penalty->GetXaxis()->SetRangeUser(0, mcMax); | |
1266 | penalty->DrawCopy("HIST"); | |
447c325d | 1267 | } |
0a173978 | 1268 | } |
1269 | ||
1270 | //____________________________________________________________________ | |
0427591c | 1271 | void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y) const |
0b4bfd98 | 1272 | { |
1273 | /*------------------------------------------------------------------------- | |
1274 | This computes an in-place complex-to-complex FFT | |
1275 | x and y are the real and imaginary arrays of 2^m points. | |
1276 | dir = 1 gives forward transform | |
1277 | dir = -1 gives reverse transform | |
1278 | ||
1279 | Formula: forward | |
1280 | N-1 | |
1281 | --- | |
1282 | 1 \ - j k 2 pi n / N | |
1283 | X(n) = --- > x(k) e = forward transform | |
1284 | N / n=0..N-1 | |
1285 | --- | |
1286 | k=0 | |
1287 | ||
1288 | Formula: reverse | |
1289 | N-1 | |
1290 | --- | |
1291 | \ j k 2 pi n / N | |
1292 | X(n) = > x(k) e = forward transform | |
1293 | / n=0..N-1 | |
1294 | --- | |
1295 | k=0 | |
1296 | */ | |
1297 | ||
1298 | Long_t nn, i, i1, j, k, i2, l, l1, l2; | |
1299 | Double_t c1, c2, tx, ty, t1, t2, u1, u2, z; | |
1300 | ||
1301 | /* Calculate the number of points */ | |
1302 | nn = 1; | |
1303 | for (i = 0; i < m; i++) | |
1304 | nn *= 2; | |
1305 | ||
1306 | /* Do the bit reversal */ | |
1307 | i2 = nn >> 1; | |
1308 | j = 0; | |
1309 | for (i= 0; i < nn - 1; i++) { | |
1310 | if (i < j) { | |
1311 | tx = x[i]; | |
1312 | ty = y[i]; | |
1313 | x[i] = x[j]; | |
1314 | y[i] = y[j]; | |
1315 | x[j] = tx; | |
1316 | y[j] = ty; | |
1317 | } | |
1318 | k = i2; | |
1319 | while (k <= j) { | |
1320 | j -= k; | |
1321 | k >>= 1; | |
1322 | } | |
1323 | j += k; | |
1324 | } | |
1325 | ||
1326 | /* Compute the FFT */ | |
1327 | c1 = -1.0; | |
1328 | c2 = 0.0; | |
1329 | l2 = 1; | |
1330 | for (l = 0; l < m; l++) { | |
1331 | l1 = l2; | |
1332 | l2 <<= 1; | |
1333 | u1 = 1.0; | |
1334 | u2 = 0.0; | |
1335 | for (j = 0;j < l1; j++) { | |
1336 | for (i = j; i < nn; i += l2) { | |
1337 | i1 = i + l1; | |
1338 | t1 = u1 * x[i1] - u2 * y[i1]; | |
1339 | t2 = u1 * y[i1] + u2 * x[i1]; | |
1340 | x[i1] = x[i] - t1; | |
1341 | y[i1] = y[i] - t2; | |
1342 | x[i] += t1; | |
1343 | y[i] += t2; | |
1344 | } | |
1345 | z = u1 * c1 - u2 * c2; | |
1346 | u2 = u1 * c2 + u2 * c1; | |
1347 | u1 = z; | |
1348 | } | |
1349 | c2 = TMath::Sqrt((1.0 - c1) / 2.0); | |
1350 | if (dir == 1) | |
1351 | c2 = -c2; | |
1352 | c1 = TMath::Sqrt((1.0 + c1) / 2.0); | |
1353 | } | |
1354 | ||
1355 | /* Scaling for forward transform */ | |
1356 | if (dir == 1) { | |
1357 | for (i=0;i<nn;i++) { | |
1358 | x[i] /= (Double_t)nn; | |
1359 | y[i] /= (Double_t)nn; | |
1360 | } | |
1361 | } | |
1362 | } | |
1363 | ||
1364 | //____________________________________________________________________ | |
9e190068 | 1365 | void AliMultiplicityCorrection::GetComparisonResults(Float_t* const mc, Int_t* const mcLimit, Float_t* const residuals, Float_t* const ratioAverage) const |
cfc19dd5 | 1366 | { |
1367 | // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD | |
1368 | // These values are computed during DrawComparison, thus this function picks up the | |
1369 | // last calculation | |
1370 | ||
dd701109 | 1371 | if (mc) |
1372 | *mc = fLastChi2MC; | |
1373 | if (mcLimit) | |
1374 | *mcLimit = fLastChi2MCLimit; | |
1375 | if (residuals) | |
1376 | *residuals = fLastChi2Residuals; | |
0b4bfd98 | 1377 | if (ratioAverage) |
1378 | *ratioAverage = fRatioAverage; | |
cfc19dd5 | 1379 | } |
1380 | ||
cfc19dd5 | 1381 | //____________________________________________________________________ |
0427591c | 1382 | TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType) const |
cfc19dd5 | 1383 | { |
1384 | // | |
1385 | // returns the corresponding MC spectrum | |
1386 | // | |
1387 | ||
1388 | switch (eventType) | |
1389 | { | |
1390 | case kTrVtx : return fMultiplicityVtx[i]; break; | |
1391 | case kMB : return fMultiplicityMB[i]; break; | |
1392 | case kINEL : return fMultiplicityINEL[i]; break; | |
69b09e3b | 1393 | case kNSD : return fMultiplicityNSD[i]; break; |
cfc19dd5 | 1394 | } |
1395 | ||
1396 | return 0; | |
1397 | } | |
1398 | ||
69b09e3b | 1399 | //____________________________________________________________________ |
9e190068 | 1400 | void AliMultiplicityCorrection::SetMultiplicityMC(Int_t i, EventType eventType, TH2F* const hist) |
69b09e3b | 1401 | { |
1402 | // | |
1403 | // returns the corresponding MC spectrum | |
1404 | // | |
1405 | ||
1406 | switch (eventType) | |
1407 | { | |
1408 | case kTrVtx : fMultiplicityVtx[i] = hist; break; | |
1409 | case kMB : fMultiplicityMB[i] = hist; break; | |
1410 | case kINEL : fMultiplicityINEL[i] = hist; break; | |
1411 | case kNSD : fMultiplicityNSD[i] = hist; break; | |
1412 | } | |
1413 | } | |
1414 | ||
0f67a57c | 1415 | //____________________________________________________________________ |
1416 | TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max) | |
1417 | { | |
1418 | // calculate standard deviation of (results[0] - results[k]) k=1...max-1 | |
1419 | // per bin one gets: sigma(r[0] - r[n]) / r[0] | |
1420 | ||
1421 | TH1* standardDeviation = (TH1*) results[0]->Clone("standardDeviation"); | |
1422 | standardDeviation->Reset(); | |
1423 | ||
1424 | for (Int_t x=1; x<=results[0]->GetNbinsX(); x++) | |
1425 | { | |
1426 | if (results[0]->GetBinContent(x) > 0) | |
1427 | { | |
1428 | Double_t average = 0; | |
1429 | for (Int_t n=1; n<max; ++n) | |
1430 | average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x); | |
1431 | average /= max-1; | |
1432 | ||
1433 | Double_t variance = 0; | |
1434 | for (Int_t n=1; n<max; ++n) | |
1435 | { | |
1436 | Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average; | |
1437 | variance += value * value; | |
1438 | } | |
1439 | variance /= max-1; | |
1440 | ||
1441 | Double_t standardDev = TMath::Sqrt(variance); | |
1442 | standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x)); | |
1443 | //Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x)); | |
1444 | } | |
1445 | } | |
1446 | ||
1447 | return standardDeviation; | |
1448 | } | |
1449 | ||
cfc19dd5 | 1450 | //____________________________________________________________________ |
eb9356d5 | 1451 | TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Int_t zeroBinEvents, Bool_t randomizeMeasured, Bool_t randomizeResponse, const TH1* compareTo) |
0a173978 | 1452 | { |
1453 | // | |
0b4bfd98 | 1454 | // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix |
1455 | // the function unfolds the spectrum using the default response matrix and several modified ones | |
1456 | // the modified ones are created by randomizing each cell using poisson statistics with the mean = bin value | |
1457 | // these unfolded results are compared to the first result gained with the default response OR to the histogram given | |
1458 | // in <compareTo> (optional) | |
0a173978 | 1459 | // |
0b4bfd98 | 1460 | // returns the error assigned to the measurement |
1461 | // | |
1462 | ||
6d81c2de | 1463 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); |
1464 | ||
0b4bfd98 | 1465 | // initialize seed with current time |
1466 | gRandom->SetSeed(0); | |
19442b86 | 1467 | |
1468 | if (methodType == AliUnfolding::kChi2Minimization) | |
eb9356d5 | 1469 | { |
1470 | Calculate0Bin(inputRange, eventType, zeroBinEvents); | |
1471 | AliUnfolding::SetMinimumInitialValue(kTRUE, 0.1); | |
1472 | } | |
1473 | ||
19442b86 | 1474 | AliUnfolding::SetUnfoldingMethod(methodType); |
0b4bfd98 | 1475 | |
eb9356d5 | 1476 | const Int_t kErrorIterations = 20; |
0b4bfd98 | 1477 | |
1478 | TH1* maxError = 0; | |
1479 | TH1* firstResult = 0; | |
1480 | ||
1481 | TH1** results = new TH1*[kErrorIterations]; | |
1482 | ||
1483 | for (Int_t n=0; n<kErrorIterations; ++n) | |
1484 | { | |
1485 | Printf("Iteration %d of %d...", n, kErrorIterations); | |
1486 | ||
19442b86 | 1487 | SetupCurrentHists(inputRange, fullPhaseSpace, eventType); |
0b4bfd98 | 1488 | |
1489 | TH1* measured = (TH1*) fCurrentESD->Clone("measured"); | |
1490 | ||
1491 | if (n > 0) | |
1492 | { | |
1493 | if (randomizeResponse) | |
1494 | { | |
1495 | // randomize response matrix | |
1496 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
1497 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) | |
1498 | fCurrentCorrelation->SetBinContent(i, j, gRandom->Poisson(fCurrentCorrelation->GetBinContent(i, j))); | |
1499 | } | |
1500 | ||
1501 | if (randomizeMeasured) | |
1502 | { | |
1503 | // randomize measured spectrum | |
1504 | for (Int_t x=1; x<=measured->GetNbinsX(); x++) // mult. axis | |
1505 | { | |
1506 | Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x)); | |
1507 | measured->SetBinContent(x, randomValue); | |
1508 | measured->SetBinError(x, TMath::Sqrt(randomValue)); | |
1509 | } | |
1510 | } | |
1511 | } | |
1512 | ||
6d81c2de | 1513 | // only for bayesian method we have to do it before the call to Unfold... |
19442b86 | 1514 | if (methodType == AliUnfolding::kBayesian) |
0b4bfd98 | 1515 | { |
6d81c2de | 1516 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) |
0b4bfd98 | 1517 | { |
6d81c2de | 1518 | // with this it is normalized to 1 |
1519 | Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); | |
1520 | ||
1521 | // with this normalized to the given efficiency | |
1522 | if (fCurrentEfficiency->GetBinContent(i) > 0) | |
1523 | sum /= fCurrentEfficiency->GetBinContent(i); | |
0b4bfd98 | 1524 | else |
6d81c2de | 1525 | sum = 0; |
1526 | ||
1527 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) | |
0b4bfd98 | 1528 | { |
6d81c2de | 1529 | if (sum > 0) |
1530 | { | |
1531 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); | |
1532 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
1533 | } | |
1534 | else | |
1535 | { | |
1536 | fCurrentCorrelation->SetBinContent(i, j, 0); | |
1537 | fCurrentCorrelation->SetBinError(i, j, 0); | |
1538 | } | |
0b4bfd98 | 1539 | } |
1540 | } | |
1541 | } | |
1542 | ||
1543 | TH1* result = 0; | |
1544 | if (n == 0 && compareTo) | |
1545 | { | |
1546 | // in this case we just store the histogram we want to compare to | |
1547 | result = (TH1*) compareTo->Clone("compareTo"); | |
1548 | result->Sumw2(); | |
1549 | } | |
1550 | else | |
6d81c2de | 1551 | { |
1552 | result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n)); | |
0b4bfd98 | 1553 | |
19442b86 | 1554 | Int_t returnCode = AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result); |
6d81c2de | 1555 | |
1556 | if (returnCode != 0) | |
eb9356d5 | 1557 | { |
1558 | n--; | |
1559 | continue; | |
1560 | } | |
6d81c2de | 1561 | } |
0b4bfd98 | 1562 | |
1563 | // normalize | |
1564 | result->Scale(1.0 / result->Integral()); | |
1565 | ||
1566 | if (n == 0) | |
1567 | { | |
1568 | firstResult = (TH1*) result->Clone("firstResult"); | |
1569 | ||
1570 | maxError = (TH1*) result->Clone("maxError"); | |
1571 | maxError->Reset(); | |
1572 | } | |
1573 | else | |
1574 | { | |
1575 | // calculate ratio | |
1576 | TH1* ratio = (TH1*) firstResult->Clone("ratio"); | |
1577 | ratio->Divide(result); | |
1578 | ||
1579 | // find max. deviation | |
1580 | for (Int_t x=1; x<=ratio->GetNbinsX(); x++) | |
1581 | maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x)))); | |
1582 | ||
1583 | delete ratio; | |
1584 | } | |
1585 | ||
1586 | results[n] = result; | |
1587 | } | |
1588 | ||
1589 | // find covariance matrix | |
1590 | // results[n] is X_x | |
1591 | // cov. matrix is M_xy = E ( (X_x - E(X_x)) * (X_y - E(X_y))), with E() = expectation value | |
1592 | ||
1593 | Int_t nBins = results[0]->GetNbinsX(); | |
1594 | Float_t lowEdge = results[0]->GetXaxis()->GetBinLowEdge(1); | |
1595 | Float_t upEdge = results[0]->GetXaxis()->GetBinUpEdge(nBins); | |
1596 | ||
1597 | // find average, E(X) | |
1598 | TProfile* average = new TProfile("average", "average", nBins, lowEdge, upEdge); | |
1599 | for (Int_t n=1; n<kErrorIterations; ++n) | |
1600 | for (Int_t x=1; x<=results[n]->GetNbinsX(); x++) | |
1601 | average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x)); | |
1602 | //new TCanvas; average->DrawClone(); | |
69b09e3b | 1603 | |
0b4bfd98 | 1604 | // find cov. matrix |
1605 | TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge); | |
1606 | ||
1607 | for (Int_t n=1; n<kErrorIterations; ++n) | |
1608 | for (Int_t x=1; x<=results[n]->GetNbinsX(); x++) | |
1609 | for (Int_t y=1; y<=results[n]->GetNbinsX(); y++) | |
1610 | { | |
1611 | // (X_x - E(X_x)) * (X_y - E(X_y) | |
1612 | Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y)); | |
1613 | covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov); | |
1614 | } | |
69b09e3b | 1615 | TCanvas* c = new TCanvas; c->cd(); covMatrix->DrawCopy("COLZ"); |
0b4bfd98 | 1616 | |
6d81c2de | 1617 | // // fill 2D histogram that contains deviation from first |
1618 | // TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01); | |
1619 | // for (Int_t n=1; n<kErrorIterations; ++n) | |
1620 | // for (Int_t x=1; x<=results[n]->GetNbinsX(); x++) | |
1621 | // deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x)); | |
1622 | // //new TCanvas; deviations->DrawCopy("COLZ"); | |
1623 | // | |
1624 | // // get standard deviation "by hand" | |
1625 | // for (Int_t x=1; x<=nBins; x++) | |
1626 | // { | |
1627 | // if (results[0]->GetBinContent(x) > 0) | |
1628 | // { | |
1629 | // TH1* proj = deviations->ProjectionY("projRMS", x, x, "e"); | |
1630 | // Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact | |
1631 | // //standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x)); | |
1632 | // Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x)); | |
1633 | // } | |
1634 | // } | |
1635 | ||
0f67a57c | 1636 | TH1* standardDeviation = CalculateStdDev(results, kErrorIterations); |
0b4bfd98 | 1637 | |
1638 | // compare maxError to RMS of profile (variable name: average) | |
1639 | // first: calculate rms in percent of value | |
1640 | TH1* rmsError = (TH1*) maxError->Clone("rmsError"); | |
1641 | rmsError->Reset(); | |
1642 | ||
1643 | // enable error to be standard deviation (see http://root.cern.ch/root/html/TProfile.html#TProfile:SetErrorOption) | |
1644 | average->SetErrorOption("s"); | |
1645 | for (Int_t x=1; x<=rmsError->GetNbinsX(); x++) | |
1646 | if (average->GetBinContent(x) > 0) | |
1647 | rmsError->SetBinContent(x, average->GetBinError(x) / average->GetBinContent(x)); | |
1648 | ||
1649 | // find maxError deviation from average (not from "first result"/mc as above) | |
1650 | TH1* maxError2 = (TH1*) maxError->Clone("maxError2"); | |
1651 | maxError2->Reset(); | |
1652 | for (Int_t n=1; n<kErrorIterations; ++n) | |
1653 | for (Int_t x=1; x<=results[n]->GetNbinsX(); x++) | |
1654 | if (average->GetBinContent(x) > 0) | |
1655 | maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x)))); | |
1656 | ||
6d81c2de | 1657 | //new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME"); |
0b4bfd98 | 1658 | |
1659 | // plot difference between average and MC/first result | |
1660 | TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio"); | |
1661 | averageFirstRatio->Reset(); | |
1662 | averageFirstRatio->Divide(results[0], average); | |
1663 | ||
6d81c2de | 1664 | //new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME"); |
1665 | //new TCanvas; averageFirstRatio->DrawCopy(); | |
0b4bfd98 | 1666 | |
69b09e3b | 1667 | static TH1* temp = 0; |
1668 | if (!temp) | |
1669 | { | |
1670 | temp = (TH1*) standardDeviation->Clone("temp"); | |
1671 | for (Int_t x=1; x<=results[0]->GetNbinsX(); x++) | |
1672 | temp->SetBinContent(x, temp->GetBinContent(x) * results[0]->GetBinContent(x)); | |
1673 | } | |
1674 | else | |
1675 | { | |
1676 | // find difference from result[0] as TH2 | |
1677 | TH2F* pulls = new TH2F("pulls", "pulls;multiplicity;difference", nBins, lowEdge, upEdge, 1000, -10, 10); | |
1678 | for (Int_t n=1; n<kErrorIterations; ++n) | |
1679 | for (Int_t x=1; x<=results[n]->GetNbinsX(); x++) | |
1680 | if (temp->GetBinContent(x) > 0) | |
1681 | pulls->Fill(results[n]->GetXaxis()->GetBinCenter(x), (results[0]->GetBinContent(x) - results[n]->GetBinContent(x)) / temp->GetBinContent(x)); | |
1682 | new TCanvas("pulls", "pulls", 800, 600); pulls->DrawCopy(); pulls->FitSlicesY(); | |
1683 | } | |
1684 | ||
0b4bfd98 | 1685 | // clean up |
1686 | for (Int_t n=0; n<kErrorIterations; ++n) | |
1687 | delete results[n]; | |
1688 | delete[] results; | |
1689 | ||
1690 | // fill into result histogram | |
0b4bfd98 | 1691 | for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i) |
1692 | fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i)); | |
cfc19dd5 | 1693 | |
0b4bfd98 | 1694 | for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i) |
eb9356d5 | 1695 | fMultiplicityESDCorrected[correlationID]->SetBinError(i, standardDeviation->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i)); |
447c325d | 1696 | |
0b4bfd98 | 1697 | return standardDeviation; |
1698 | } | |
1699 | ||
1700 | //____________________________________________________________________ | |
eb9356d5 | 1701 | void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Int_t determineError) |
0b4bfd98 | 1702 | { |
1703 | // | |
1704 | // correct spectrum using bayesian method | |
1705 | // | |
eb9356d5 | 1706 | // determineError: |
1707 | // 0 = no errors | |
1708 | // 1 = from randomizing | |
1709 | // 2 = with UnfoldGetBias | |
0b4bfd98 | 1710 | |
1711 | // initialize seed with current time | |
1712 | gRandom->SetSeed(0); | |
1713 | ||
19442b86 | 1714 | SetupCurrentHists(inputRange, fullPhaseSpace, eventType); |
0a173978 | 1715 | |
1716 | // normalize correction for given nPart | |
9ca6be09 | 1717 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) |
0a173978 | 1718 | { |
dd701109 | 1719 | // with this it is normalized to 1 |
1720 | Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); | |
1721 | ||
1722 | // with this normalized to the given efficiency | |
1723 | if (fCurrentEfficiency->GetBinContent(i) > 0) | |
1724 | sum /= fCurrentEfficiency->GetBinContent(i); | |
1725 | else | |
1726 | sum = 0; | |
1727 | ||
9ca6be09 | 1728 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) |
0a173978 | 1729 | { |
dd701109 | 1730 | if (sum > 0) |
1731 | { | |
1732 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); | |
1733 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
1734 | } | |
1735 | else | |
1736 | { | |
1737 | fCurrentCorrelation->SetBinContent(i, j, 0); | |
1738 | fCurrentCorrelation->SetBinError(i, j, 0); | |
1739 | } | |
0a173978 | 1740 | } |
1741 | } | |
1742 | ||
0b4bfd98 | 1743 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); |
1744 | ||
19442b86 | 1745 | AliUnfolding::SetBayesianParameters(regPar, nIterations); |
1746 | AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian); | |
eb9356d5 | 1747 | |
1748 | if (determineError <= 1) | |
1749 | { | |
1750 | if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]) != 0) | |
1751 | return; | |
1752 | } | |
1753 | else if (determineError == 2) | |
1754 | { | |
1755 | AliUnfolding::UnfoldGetBias(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]); | |
0b4bfd98 | 1756 | return; |
eb9356d5 | 1757 | } |
0b4bfd98 | 1758 | |
eb9356d5 | 1759 | if (determineError == 0) |
447c325d | 1760 | { |
0b4bfd98 | 1761 | Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated."); |
1762 | return; | |
1763 | } | |
447c325d | 1764 | |
0b4bfd98 | 1765 | // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics |
1766 | // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured | |
0f67a57c | 1767 | // spectrum calculated. This is performed N times and the sigma is taken as the statistical |
0b4bfd98 | 1768 | // error of the unfolding method itself. |
1769 | ||
0b4bfd98 | 1770 | const Int_t kErrorIterations = 20; |
1771 | ||
eb9356d5 | 1772 | Printf("Spectrum unfolded. Determining error (%d iterations)...", kErrorIterations); |
0b4bfd98 | 1773 | |
1774 | TH1* randomized = (TH1*) fCurrentESD->Clone("randomized"); | |
0f67a57c | 1775 | TH1* resultArray[kErrorIterations+1]; |
0b4bfd98 | 1776 | for (Int_t n=0; n<kErrorIterations; ++n) |
1777 | { | |
1778 | // randomize the content of clone following a poisson with the mean = the value of that bin | |
1779 | for (Int_t x=1; x<=randomized->GetNbinsX(); x++) // mult. axis | |
447c325d | 1780 | { |
0b4bfd98 | 1781 | Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x)); |
1782 | //printf("%d --> %d\n", fCurrentESD->GetBinContent(x), randomValue); | |
1783 | randomized->SetBinContent(x, randomValue); | |
1784 | randomized->SetBinError(x, TMath::Sqrt(randomValue)); | |
447c325d | 1785 | } |
447c325d | 1786 | |
0f67a57c | 1787 | TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2"); |
6d81c2de | 1788 | result2->Reset(); |
19442b86 | 1789 | if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2) != 0) |
eb9356d5 | 1790 | { |
1791 | n--; | |
1792 | continue; | |
1793 | } | |
cfc19dd5 | 1794 | |
0f67a57c | 1795 | resultArray[n+1] = result2; |
0b4bfd98 | 1796 | } |
6d81c2de | 1797 | delete randomized; |
0f67a57c | 1798 | |
1799 | resultArray[0] = fMultiplicityESDCorrected[correlationID]; | |
1800 | TH1* error = CalculateStdDev(resultArray, kErrorIterations+1); | |
1801 | ||
1802 | for (Int_t n=0; n<kErrorIterations; ++n) | |
1803 | delete resultArray[n+1]; | |
447c325d | 1804 | |
eb9356d5 | 1805 | Printf("Comparing bias and error:"); |
0b4bfd98 | 1806 | for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i) |
eb9356d5 | 1807 | { |
1808 | Printf("Bin %d: Content: %f Error: %f Bias: %f", i, fMultiplicityESDCorrected[correlationID]->GetBinContent(i), error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i), fMultiplicityESDCorrected[correlationID]->GetBinError(i)); | |
0f67a57c | 1809 | fMultiplicityESDCorrected[correlationID]->SetBinError(i, error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i)); |
eb9356d5 | 1810 | } |
447c325d | 1811 | |
0f67a57c | 1812 | delete error; |
0b4bfd98 | 1813 | } |
447c325d | 1814 | |
dd701109 | 1815 | //____________________________________________________________________ |
0427591c | 1816 | Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], const TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u) |
dd701109 | 1817 | { |
1818 | // | |
1819 | // helper function for the covariance matrix of the bayesian method | |
1820 | // | |
1821 | ||
1822 | Float_t result = 0; | |
1823 | ||
1824 | if (k == u && r == i) | |
1825 | result += 1.0 / hResponse->GetBinContent(u+1, r+1); | |
1826 | ||
1827 | if (k == u) | |
1828 | result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1); | |
1829 | ||
1830 | if (r == i) | |
1831 | result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1); | |
1832 | ||
1833 | result *= matrixM[k][i]; | |
1834 | ||
1835 | return result; | |
cfc19dd5 | 1836 | } |
1837 | ||
1838 | //____________________________________________________________________ | |
1839 | void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType) | |
1840 | { | |
1841 | // | |
1842 | // correct spectrum using bayesian method | |
1843 | // | |
1844 | ||
dd701109 | 1845 | Float_t regPar = 0; |
cfc19dd5 | 1846 | |
335a5e1b | 1847 | Int_t correlationID = inputRange; // + ((fullPhaseSpace == kFALSE) ? 0 : 4); |
1848 | Int_t mcTarget = inputRange; //((fullPhaseSpace == kFALSE) ? inputRange : 4); | |
cfc19dd5 | 1849 | |
19442b86 | 1850 | SetupCurrentHists(inputRange, fullPhaseSpace, eventType); |
0b4bfd98 | 1851 | //normalize ESD |
1852 | fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); | |
cfc19dd5 | 1853 | |
1854 | // TODO should be taken from correlation map | |
1855 | //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX()); | |
1856 | ||
1857 | // normalize correction for given nPart | |
1858 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
1859 | { | |
1860 | Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); | |
1861 | //Double_t sum = sumHist->GetBinContent(i); | |
1862 | if (sum <= 0) | |
1863 | continue; | |
1864 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) | |
1865 | { | |
1866 | // npart sum to 1 | |
1867 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); | |
1868 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
1869 | } | |
6127aca6 | 1870 | } |
cfc19dd5 | 1871 | |
1872 | new TCanvas; | |
1873 | fCurrentCorrelation->Draw("COLZ"); | |
1874 | ||
1875 | // FAKE | |
1876 | fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff"); | |
1877 | ||
1878 | // pick prior distribution | |
1879 | TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior"); | |
1880 | Float_t norm = 1; //hPrior->Integral(); | |
1881 | for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) | |
1882 | hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm); | |
1883 | ||
1884 | // zero distribution | |
1885 | TH1F* zero = (TH1F*)hPrior->Clone("zero"); | |
1886 | ||
1887 | // define temp hist | |
1888 | TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp"); | |
1889 | hTemp->Reset(); | |
1890 | ||
1891 | // just a shortcut | |
1892 | TH2F* hResponse = (TH2F*) fCurrentCorrelation; | |
1893 | ||
1894 | // unfold... | |
dd701109 | 1895 | Int_t iterations = 25; |
cfc19dd5 | 1896 | for (Int_t i=0; i<iterations; i++) |
1897 | { | |
1898 | //printf(" iteration %i \n", i); | |
1899 | ||
1900 | for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) | |
1901 | { | |
1902 | Float_t value = 0; | |
1903 | for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++) | |
1904 | value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t); | |
1905 | hTemp->SetBinContent(m, value); | |
1906 | //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value); | |
1907 | } | |
1908 | ||
1909 | // regularization (simple smoothing) | |
1910 | TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed"); | |
1911 | ||
1912 | for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++) | |
1913 | { | |
1914 | Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1) | |
1915 | + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t) | |
1916 | + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.; | |
1917 | average *= hTrueSmoothed->GetBinWidth(t); | |
1918 | ||
1919 | // weight the average with the regularization parameter | |
1920 | hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average); | |
1921 | } | |
1922 | ||
1923 | for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) | |
1924 | hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m)); | |
1925 | ||
1926 | // fill guess | |
1927 | for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++) | |
1928 | { | |
1929 | fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t)); | |
1930 | fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO | |
1931 | ||
1932 | //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t)); | |
1933 | } | |
1934 | ||
1935 | ||
1936 | // calculate chi2 (change from last iteration) | |
1937 | Double_t chi2 = 0; | |
1938 | ||
1939 | // use smoothed true (normalized) as new prior | |
745d6088 | 1940 | norm = 1; //hTrueSmoothed->Integral(); |
cfc19dd5 | 1941 | |
1942 | for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++) | |
1943 | { | |
dd701109 | 1944 | Float_t newValue = hTemp->GetBinContent(t)/norm; |
cfc19dd5 | 1945 | Float_t diff = hPrior->GetBinContent(t) - newValue; |
1946 | chi2 += (Double_t) diff * diff; | |
1947 | ||
1948 | hPrior->SetBinContent(t, newValue); | |
1949 | } | |
1950 | ||
1951 | printf("Chi2 of %d iteration = %.10f\n", i, chi2); | |
1952 | ||
dd701109 | 1953 | //if (i % 5 == 0) |
cfc19dd5 | 1954 | DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY()); |
1955 | ||
1956 | delete hTrueSmoothed; | |
1957 | } // end of iterations | |
1958 | ||
1959 | DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY()); | |
1960 | } | |
1961 | ||
9ca6be09 | 1962 | //____________________________________________________________________ |
1963 | void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace) | |
1964 | { | |
1965 | // | |
1966 | // correct spectrum using a simple Gaussian approach, that is model-dependent | |
1967 | // | |
1968 | ||
335a5e1b | 1969 | Int_t correlationID = inputRange; // + ((fullPhaseSpace == kFALSE) ? 0 : 4); |
1970 | Int_t mcTarget = inputRange; //((fullPhaseSpace == kFALSE) ? inputRange : 4); | |
9ca6be09 | 1971 | |
19442b86 | 1972 | SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx); |
0b4bfd98 | 1973 | //normalize ESD |
1974 | fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); | |
9ca6be09 | 1975 | |
9ca6be09 | 1976 | TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean")); |
1977 | correction->SetTitle("GaussianMean"); | |
1978 | ||
1979 | TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth")); | |
1980 | correction->SetTitle("GaussianWidth"); | |
1981 | ||
1982 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
1983 | { | |
1984 | TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1); | |
1985 | proj->Fit("gaus", "0Q"); | |
1986 | correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1)); | |
1987 | correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2)); | |
ff458f69 | 1988 | /* |
9ca6be09 | 1989 | // draw for debugging |
1990 | new TCanvas; | |
1991 | proj->DrawCopy(); | |
1992 | proj->GetFunction("gaus")->DrawCopy("SAME"); | |
ff458f69 | 1993 | */ |
9ca6be09 | 1994 | } |
1995 | ||
1996 | TH1* target = fMultiplicityESDCorrected[correlationID]; | |
1997 | ||
1998 | Int_t nBins = target->GetNbinsX()*10+1; | |
1999 | Float_t* binning = new Float_t[nBins]; | |
2000 | for (Int_t i=1; i<=target->GetNbinsX(); ++i) | |
2001 | for (Int_t j=0; j<10; ++j) | |
2002 | binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j; | |
2003 | ||
2004 | binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX()); | |
2005 | ||
2006 | TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning); | |
2007 | ||
2008 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
2009 | { | |
2010 | Float_t mean = correction->GetBinContent(i); | |
2011 | Float_t width = correctionWidth->GetBinContent(i); | |
2012 | ||
3602328d | 2013 | Int_t fillBegin = fineBinned->FindBin(mean - width * 5); |
2014 | Int_t fillEnd = fineBinned->FindBin(mean + width * 5); | |
2015 | //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd); | |
9ca6be09 | 2016 | |
2017 | for (Int_t j=fillBegin; j <= fillEnd; ++j) | |
2018 | { | |
2019 | fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i)); | |
2020 | } | |
2021 | } | |
2022 | ||
2023 | for (Int_t i=1; i<=target->GetNbinsX(); ++i) | |
2024 | { | |
2025 | Float_t sum = 0; | |
2026 | for (Int_t j=1; j<=10; ++j) | |
2027 | sum += fineBinned->GetBinContent((i-1)*10 + j); | |
2028 | target->SetBinContent(i, sum / 10); | |
2029 | } | |
2030 | ||
2031 | delete[] binning; | |
2032 | ||
cfc19dd5 | 2033 | DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY()); |
0a173978 | 2034 | } |
3602328d | 2035 | |
eb9356d5 | 2036 | //____________________________________________________________________ |
2037 | TH1* AliMultiplicityCorrection::GetConvoluted(Int_t i, EventType eventType) | |
2038 | { | |
2039 | // convolutes the corrected histogram i with the response matrix | |
2040 | ||
2041 | TH1* corrected = (TH1*) fMultiplicityESDCorrected[i]->Clone("corrected"); | |
2042 | ||
2043 | // undo efficiency correction (hist must not be deleted, is reused) | |
2044 | TH1* efficiency = GetEfficiency(i, eventType); | |
2045 | //new TCanvas; efficiency->DrawCopy(); | |
2046 | corrected->Multiply(efficiency); | |
2047 | ||
2048 | TH2* convoluted = CalculateMultiplicityESD(corrected, i); | |
2049 | TH1* convolutedProj = convoluted->ProjectionY("GetConvoluted_convolutedProj", 1, convoluted->GetNbinsX()); | |
2050 | ||
2051 | delete convoluted; | |
2052 | delete corrected; | |
2053 | ||
2054 | return convolutedProj; | |
2055 | } | |
2056 | ||
2057 | //____________________________________________________________________ | |
2058 | TH1* AliMultiplicityCorrection::GetResiduals(Int_t i, EventType eventType, Float_t& residualSum) | |
2059 | { | |
2060 | // creates the residual histogram from the corrected histogram i corresponding to an eventType event sample using the corresponding correlation matrix | |
2061 | // residual is : M - UT / eM | |
2062 | // residualSum contains the squared sum of the residuals | |
2063 | ||
2064 | TH1* corrected = (TH1*) fMultiplicityESDCorrected[i]->Clone("corrected"); | |
2065 | TH1* convolutedProj = GetConvoluted(i, eventType); | |
2066 | ||
2067 | Int_t begin = 1; | |
2068 | Int_t end = fMultiplicityESD[i]->GetNbinsX(); | |
2069 | if (fVtxEnd > fVtxBegin) | |
2070 | { | |
2071 | begin = fVtxBegin; | |
2072 | end = fVtxEnd; | |
2073 | } | |
2074 | TH1* measuredProj = fMultiplicityESD[i]->ProjectionY("measuredProj", begin, end); | |
2075 | ||
2076 | TH1* residuals = (TH1*) measuredProj->Clone("GetResiduals_residuals"); | |
2077 | residuals->SetTitle(";measured multiplicity;residuals (M-Ut)/e"); | |
2078 | residuals->Add(convolutedProj, -1); | |
2079 | ||
2080 | residualSum = 0; | |
d29b31c5 | 2081 | for (Int_t j=1; j<=residuals->GetNbinsX(); j++) |
eb9356d5 | 2082 | { |
d29b31c5 | 2083 | if (measuredProj->GetBinContent(j) > 0) |
2084 | residuals->SetBinContent(j, residuals->GetBinContent(j) / TMath::Sqrt(measuredProj->GetBinContent(j))); | |
2085 | residuals->SetBinError(j, 0); | |
eb9356d5 | 2086 | |
d29b31c5 | 2087 | if (j > 1) |
2088 | residualSum += residuals->GetBinContent(j) * residuals->GetBinContent(j); | |
eb9356d5 | 2089 | } |
2090 | ||
2091 | delete corrected; | |
2092 | delete convolutedProj; | |
2093 | delete measuredProj; | |
2094 | ||
2095 | return residuals; | |
2096 | } | |
2097 | ||
3602328d | 2098 | //____________________________________________________________________ |
447c325d | 2099 | TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap) |
3602328d | 2100 | { |
2101 | // runs the distribution given in inputMC through the response matrix identified by | |
2102 | // correlationMap and produces a measured distribution | |
2103 | // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex | |
2104 | ||
2105 | if (!inputMC) | |
2106 | return 0; | |
2107 | ||
2108 | if (correlationMap < 0 || correlationMap >= kCorrHists) | |
2109 | return 0; | |
2110 | ||
2111 | // empty under/overflow bins in x, otherwise Project3D takes them into account | |
2112 | TH3* hist = fCorrelation[correlationMap]; | |
039e265e | 2113 | for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y) |
3602328d | 2114 | { |
039e265e | 2115 | for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z) |
3602328d | 2116 | { |
2117 | hist->SetBinContent(0, y, z, 0); | |
2118 | hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); | |
2119 | } | |
2120 | } | |
2121 | ||
eb9356d5 | 2122 | if (fVtxEnd > fVtxBegin) |
2123 | hist->GetXaxis()->SetRange(fVtxBegin, fVtxEnd); | |
2124 | ||
447c325d | 2125 | TH2* corr = (TH2*) hist->Project3D("zy"); |
2126 | //corr->Rebin2D(2, 1); | |
3602328d | 2127 | corr->Sumw2(); |
039e265e | 2128 | Printf("Correction histogram used for convolution has %f entries", corr->Integral()); |
3602328d | 2129 | |
2130 | // normalize correction for given nPart | |
2131 | for (Int_t i=1; i<=corr->GetNbinsX(); ++i) | |
2132 | { | |
2133 | Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY()); | |
2134 | if (sum <= 0) | |
2135 | continue; | |
2136 | ||
2137 | for (Int_t j=1; j<=corr->GetNbinsY(); ++j) | |
2138 | { | |
2139 | // npart sum to 1 | |
2140 | corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum); | |
2141 | corr->SetBinError(i, j, corr->GetBinError(i, j) / sum); | |
2142 | } | |
2143 | } | |
2144 | ||
335a5e1b | 2145 | TH2F* target = static_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName()))); |
3602328d | 2146 | target->Reset(); |
2147 | ||
2148 | for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas) | |
2149 | { | |
2150 | Float_t measured = 0; | |
2151 | Float_t error = 0; | |
2152 | ||
447c325d | 2153 | for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen) |
3602328d | 2154 | { |
447c325d | 2155 | Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen)); |
b477f4f2 | 2156 | |
447c325d | 2157 | measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas); |
eb9356d5 | 2158 | // TODO fix error |
447c325d | 2159 | error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas); |
3602328d | 2160 | } |
2161 | ||
cfc19dd5 | 2162 | //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error); |
3602328d | 2163 | |
447c325d | 2164 | target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured); |
2165 | target->SetBinError(1 + target->GetNbinsX() / 2, meas, error); | |
3602328d | 2166 | } |
2167 | ||
2168 | return target; | |
2169 | } | |
2170 | ||
2171 | //____________________________________________________________________ | |
0427591c | 2172 | void AliMultiplicityCorrection::SetGenMeasFromFunc(const TF1* inputMC, Int_t id) |
3602328d | 2173 | { |
2174 | // uses the given function to fill the input MC histogram and generates from that | |
2175 | // the measured histogram by applying the response matrix | |
2176 | // this can be used to evaluate if the methods work indepedently of the input | |
2177 | // distribution | |
2178 | // WARNING does not respect the vertex distribution, just fills central vertex bin | |
2179 | ||
2180 | if (!inputMC) | |
2181 | return; | |
2182 | ||
2183 | if (id < 0 || id >= kESDHists) | |
2184 | return; | |
2185 | ||
69b09e3b | 2186 | // fill histogram used for random generation |
2187 | TH1* tmp = fMultiplicityVtx[id]->ProjectionY("tmp"); | |
2188 | tmp->Reset(); | |
2189 | ||
2190 | for (Int_t i=1; i<=tmp->GetNbinsX(); ++i) | |
2191 | tmp->SetBinContent(i, inputMC->Eval(tmp->GetXaxis()->GetBinCenter(i)) * tmp->GetXaxis()->GetBinWidth(i)); | |
2192 | ||
2193 | TH1* mcRnd = fMultiplicityVtx[id]->ProjectionY("mcRnd"); | |
2194 | mcRnd->Reset(); | |
ad0385e9 | 2195 | mcRnd->FillRandom(tmp, (Int_t) tmp->Integral()); |
69b09e3b | 2196 | |
2197 | //new TCanvas; tmp->Draw(); | |
2198 | //new TCanvas; mcRnd->Draw(); | |
2199 | ||
2200 | // and move into 2d histogram | |
2201 | TH1* mc = fMultiplicityVtx[id]; | |
3602328d | 2202 | mc->Reset(); |
3602328d | 2203 | for (Int_t i=1; i<=mc->GetNbinsY(); ++i) |
b477f4f2 | 2204 | { |
69b09e3b | 2205 | mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, mcRnd->GetBinContent(i)); |
2206 | mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, TMath::Sqrt(mcRnd->GetBinContent(i))); | |
2207 | } | |
2208 | ||
2209 | //new TCanvas; mc->Draw("COLZ"); | |
2210 | ||
2211 | // now randomize the measured histogram; funcMeasured is used as pilot function to generated the measured entries | |
2212 | TH1* funcMeasured = CalculateMultiplicityESD(tmp, id)->ProjectionY("funcMeasured"); | |
2213 | ||
2214 | //new TCanvas; funcMeasured->Draw(); | |
2215 | ||
2216 | fMultiplicityESD[id]->Reset(); | |
2217 | ||
2218 | TH1* measRnd = fMultiplicityESD[id]->ProjectionY("measRnd"); | |
ad0385e9 | 2219 | measRnd->FillRandom(funcMeasured, (Int_t) tmp->Integral()); |
69b09e3b | 2220 | |
2221 | //new TCanvas; measRnd->Draw(); | |
2222 | ||
2223 | fMultiplicityESD[id]->Reset(); | |
2224 | for (Int_t i=1; i<=fMultiplicityESD[id]->GetNbinsY(); ++i) | |
2225 | { | |
2226 | fMultiplicityESD[id]->SetBinContent(fMultiplicityESD[id]->GetNbinsX() / 2 + 1, i, measRnd->GetBinContent(i)); | |
2227 | fMultiplicityESD[id]->SetBinError(fMultiplicityESD[id]->GetNbinsX() / 2 + 1, i, TMath::Sqrt(measRnd->GetBinContent(i))); | |
b477f4f2 | 2228 | } |
3602328d | 2229 | } |