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