]>
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> | |
32 | #include <TVirtualFitter.h> | |
33 | #include <TCanvas.h> | |
34 | #include <TString.h> | |
9ca6be09 | 35 | #include <TF1.h> |
d560b581 | 36 | #include <TMath.h> |
3602328d | 37 | #include <TCollection.h> |
447c325d | 38 | #include <TLegend.h> |
39 | #include <TLine.h> | |
0b4bfd98 | 40 | #include <TRandom.h> |
41 | #include <TProfile.h> | |
42 | #include <TProfile2D.h> | |
0a173978 | 43 | |
44 | ClassImp(AliMultiplicityCorrection) | |
45 | ||
6d81c2de | 46 | const Int_t AliMultiplicityCorrection::fgkMaxInput = 250; // bins in measured histogram |
47 | const Int_t AliMultiplicityCorrection::fgkMaxParams = 250; // bins in unfolded histogram = number of fit params | |
447c325d | 48 | |
6d81c2de | 49 | TMatrixD* AliMultiplicityCorrection::fgCorrelationMatrix = 0; |
50 | TMatrixD* AliMultiplicityCorrection::fgCorrelationCovarianceMatrix = 0; | |
51 | TVectorD* AliMultiplicityCorrection::fgCurrentESDVector = 0; | |
52 | TVectorD* AliMultiplicityCorrection::fgEntropyAPriori = 0; | |
53 | ||
44466df2 | 54 | Bool_t AliMultiplicityCorrection::fgCreateBigBin = kTRUE; |
6d81c2de | 55 | AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fgRegularizationType = AliMultiplicityCorrection::kPol1; |
56 | Float_t AliMultiplicityCorrection::fgRegularizationWeight = 10000; | |
44466df2 | 57 | Int_t AliMultiplicityCorrection::fgNParamsMinuit = AliMultiplicityCorrection::fgkMaxParams; |
6d81c2de | 58 | |
59 | TF1* AliMultiplicityCorrection::fgNBD = 0; | |
60 | ||
61 | Float_t AliMultiplicityCorrection::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing) | |
62 | Int_t AliMultiplicityCorrection::fgBayesianIterations = 100; // number of iterations in Bayesian method | |
0a173978 | 63 | |
0b4bfd98 | 64 | // These are the areas where the quality of the unfolding results are evaluated |
65 | // Default defined here, call SetQualityRegions to change them | |
66 | // unit is in multiplicity (not in bin!) | |
67 | ||
68 | // SPD: peak area - flat area - low stat area | |
69 | Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {4, 60, 180}; | |
70 | Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {20, 140, 210}; | |
71 | ||
72 | //____________________________________________________________________ | |
73 | void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy) | |
74 | { | |
75 | // | |
76 | // sets the quality region definition to TPC or SPD | |
77 | // | |
78 | ||
79 | if (SPDStudy) | |
80 | { | |
81 | // SPD: peak area - flat area - low stat area | |
82 | fgQualityRegionsB[0] = 4; | |
83 | fgQualityRegionsE[0] = 20; | |
84 | ||
85 | fgQualityRegionsB[1] = 60; | |
86 | fgQualityRegionsE[1] = 140; | |
87 | ||
88 | fgQualityRegionsB[2] = 180; | |
89 | fgQualityRegionsE[2] = 210; | |
6d81c2de | 90 | |
91 | Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for SPD"); | |
0b4bfd98 | 92 | } |
93 | else | |
94 | { | |
95 | // TPC: peak area - flat area - low stat area | |
96 | fgQualityRegionsB[0] = 4; | |
97 | fgQualityRegionsE[0] = 12; | |
98 | ||
99 | fgQualityRegionsB[1] = 25; | |
100 | fgQualityRegionsE[1] = 55; | |
101 | ||
102 | fgQualityRegionsB[2] = 88; | |
103 | fgQualityRegionsE[2] = 108; | |
6d81c2de | 104 | |
105 | Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for TPC"); | |
0b4bfd98 | 106 | } |
107 | } | |
108 | ||
0a173978 | 109 | //____________________________________________________________________ |
110 | AliMultiplicityCorrection::AliMultiplicityCorrection() : | |
0f67a57c | 111 | TNamed(), fCurrentESD(0), fCurrentCorrelation(0), fCurrentEfficiency(0), fLastBinLimit(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0) |
0a173978 | 112 | { |
113 | // | |
114 | // default constructor | |
115 | // | |
116 | ||
117 | for (Int_t i = 0; i < kESDHists; ++i) | |
118 | fMultiplicityESD[i] = 0; | |
119 | ||
120 | for (Int_t i = 0; i < kMCHists; ++i) | |
cfc19dd5 | 121 | { |
122 | fMultiplicityVtx[i] = 0; | |
123 | fMultiplicityMB[i] = 0; | |
124 | fMultiplicityINEL[i] = 0; | |
125 | } | |
0a173978 | 126 | |
127 | for (Int_t i = 0; i < kCorrHists; ++i) | |
128 | { | |
129 | fCorrelation[i] = 0; | |
130 | fMultiplicityESDCorrected[i] = 0; | |
131 | } | |
6d81c2de | 132 | |
133 | for (Int_t i = 0; i < kQualityRegions; ++i) | |
134 | fQuality[i] = 0; | |
0a173978 | 135 | } |
136 | ||
137 | //____________________________________________________________________ | |
138 | AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) : | |
7307d52c | 139 | TNamed(name, title), |
140 | fCurrentESD(0), | |
141 | fCurrentCorrelation(0), | |
142 | fCurrentEfficiency(0), | |
143 | fLastBinLimit(0), | |
144 | fLastChi2MC(0), | |
145 | fLastChi2MCLimit(0), | |
146 | fLastChi2Residuals(0), | |
147 | fRatioAverage(0) | |
0a173978 | 148 | { |
149 | // | |
150 | // named constructor | |
151 | // | |
152 | ||
153 | // do not add this hists to the directory | |
154 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
155 | TH1::AddDirectory(kFALSE); | |
156 | ||
b477f4f2 | 157 | /*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 | 158 | 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, |
159 | 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, | |
160 | 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, | |
161 | 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, | |
162 | 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, | |
163 | 50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5, | |
164 | 100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5, | |
165 | 150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5, | |
166 | 250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5, | |
9ca6be09 | 167 | 500.5 }; |
168 | //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5, | |
169 | //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5, | |
170 | //1000.5 }; | |
0a173978 | 171 | |
b477f4f2 | 172 | #define VTXBINNING 10, binLimitsVtx |
6d81c2de | 173 | #define NBINNING fgkMaxParams, binLimitsN*/ |
b477f4f2 | 174 | |
447c325d | 175 | #define NBINNING 501, -0.5, 500.5 |
176 | #define VTXBINNING 1, -10, 10 | |
0a173978 | 177 | |
178 | for (Int_t i = 0; i < kESDHists; ++i) | |
b477f4f2 | 179 | fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING); |
0a173978 | 180 | |
181 | for (Int_t i = 0; i < kMCHists; ++i) | |
182 | { | |
cfc19dd5 | 183 | fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i))); |
184 | fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart"); | |
185 | ||
186 | fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i))); | |
187 | fMultiplicityMB[i]->SetTitle("fMultiplicityMB"); | |
188 | ||
189 | fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i))); | |
190 | fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL"); | |
0a173978 | 191 | } |
192 | ||
193 | for (Int_t i = 0; i < kCorrHists; ++i) | |
194 | { | |
b477f4f2 | 195 | fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING); |
196 | fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING); | |
0a173978 | 197 | } |
198 | ||
199 | TH1::AddDirectory(oldStatus); | |
200 | } | |
201 | ||
202 | //____________________________________________________________________ | |
203 | AliMultiplicityCorrection::~AliMultiplicityCorrection() | |
204 | { | |
205 | // | |
206 | // Destructor | |
207 | // | |
6d81c2de | 208 | |
0f67a57c | 209 | Printf("AliMultiplicityCorrection::~AliMultiplicityCorrection called"); |
210 | ||
6d81c2de | 211 | for (Int_t i = 0; i < kESDHists; ++i) |
212 | { | |
213 | if (fMultiplicityESD[i]) | |
214 | delete fMultiplicityESD[i]; | |
215 | fMultiplicityESD[i] = 0; | |
216 | } | |
217 | ||
218 | for (Int_t i = 0; i < kMCHists; ++i) | |
219 | { | |
220 | if (fMultiplicityVtx[i]) | |
221 | delete fMultiplicityVtx[i]; | |
222 | fMultiplicityVtx[i] = 0; | |
223 | ||
224 | if (fMultiplicityMB[i]) | |
225 | delete fMultiplicityMB[i]; | |
226 | fMultiplicityMB[i] = 0; | |
227 | ||
228 | if (fMultiplicityINEL[i]) | |
229 | delete fMultiplicityINEL[i]; | |
230 | fMultiplicityINEL[i] = 0; | |
231 | } | |
232 | ||
233 | for (Int_t i = 0; i < kCorrHists; ++i) | |
234 | { | |
235 | if (fCorrelation[i]) | |
236 | delete fCorrelation[i]; | |
237 | fCorrelation[i] = 0; | |
238 | ||
239 | if (fMultiplicityESDCorrected[i]) | |
240 | delete fMultiplicityESDCorrected[i]; | |
241 | fMultiplicityESDCorrected[i] = 0; | |
242 | } | |
0a173978 | 243 | } |
244 | ||
245 | //____________________________________________________________________ | |
246 | Long64_t AliMultiplicityCorrection::Merge(TCollection* list) | |
247 | { | |
248 | // Merge a list of AliMultiplicityCorrection objects with this (needed for | |
249 | // PROOF). | |
250 | // Returns the number of merged objects (including this). | |
251 | ||
252 | if (!list) | |
253 | return 0; | |
254 | ||
255 | if (list->IsEmpty()) | |
256 | return 1; | |
257 | ||
258 | TIterator* iter = list->MakeIterator(); | |
259 | TObject* obj; | |
260 | ||
261 | // collections of all histograms | |
cfc19dd5 | 262 | TList collections[kESDHists+kMCHists*3+kCorrHists*2]; |
0a173978 | 263 | |
264 | Int_t count = 0; | |
265 | while ((obj = iter->Next())) { | |
266 | ||
267 | AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj); | |
268 | if (entry == 0) | |
269 | continue; | |
270 | ||
271 | for (Int_t i = 0; i < kESDHists; ++i) | |
272 | collections[i].Add(entry->fMultiplicityESD[i]); | |
273 | ||
274 | for (Int_t i = 0; i < kMCHists; ++i) | |
cfc19dd5 | 275 | { |
276 | collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]); | |
277 | collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]); | |
278 | collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]); | |
279 | } | |
0a173978 | 280 | |
281 | for (Int_t i = 0; i < kCorrHists; ++i) | |
cfc19dd5 | 282 | collections[kESDHists+kMCHists*3+i].Add(entry->fCorrelation[i]); |
0a173978 | 283 | |
284 | for (Int_t i = 0; i < kCorrHists; ++i) | |
cfc19dd5 | 285 | collections[kESDHists+kMCHists*3+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]); |
0a173978 | 286 | |
287 | count++; | |
288 | } | |
289 | ||
290 | for (Int_t i = 0; i < kESDHists; ++i) | |
291 | fMultiplicityESD[i]->Merge(&collections[i]); | |
292 | ||
293 | for (Int_t i = 0; i < kMCHists; ++i) | |
cfc19dd5 | 294 | { |
295 | fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]); | |
296 | fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]); | |
297 | fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]); | |
298 | } | |
0a173978 | 299 | |
300 | for (Int_t i = 0; i < kCorrHists; ++i) | |
cfc19dd5 | 301 | fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*3+i]); |
0a173978 | 302 | |
303 | for (Int_t i = 0; i < kCorrHists; ++i) | |
cfc19dd5 | 304 | fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*3+kCorrHists+i]); |
0a173978 | 305 | |
306 | delete iter; | |
307 | ||
308 | return count+1; | |
309 | } | |
310 | ||
311 | //____________________________________________________________________ | |
312 | Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir) | |
313 | { | |
314 | // | |
315 | // loads the histograms from a file | |
316 | // if dir is empty a directory with the name of this object is taken (like in SaveHistogram) | |
317 | // | |
318 | ||
319 | if (!dir) | |
320 | dir = GetName(); | |
321 | ||
322 | if (!gDirectory->cd(dir)) | |
323 | return kFALSE; | |
324 | ||
144ff489 | 325 | // store old hists to delete them later |
326 | TList oldObjects; | |
327 | oldObjects.SetOwner(1); | |
328 | for (Int_t i = 0; i < kESDHists; ++i) | |
329 | if (fMultiplicityESD[i]) | |
330 | oldObjects.Add(fMultiplicityESD[i]); | |
331 | ||
332 | for (Int_t i = 0; i < kMCHists; ++i) | |
333 | { | |
334 | if (fMultiplicityVtx[i]) | |
335 | oldObjects.Add(fMultiplicityVtx[i]); | |
336 | if (fMultiplicityMB[i]) | |
337 | oldObjects.Add(fMultiplicityMB[i]); | |
338 | if (fMultiplicityINEL[i]) | |
339 | oldObjects.Add(fMultiplicityINEL[i]); | |
340 | } | |
341 | ||
342 | for (Int_t i = 0; i < kCorrHists; ++i) | |
343 | if (fCorrelation[i]) | |
344 | oldObjects.Add(fCorrelation[i]); | |
345 | ||
346 | // load histograms | |
0a173978 | 347 | |
348 | Bool_t success = kTRUE; | |
349 | ||
350 | for (Int_t i = 0; i < kESDHists; ++i) | |
351 | { | |
352 | fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName())); | |
353 | if (!fMultiplicityESD[i]) | |
354 | success = kFALSE; | |
355 | } | |
356 | ||
357 | for (Int_t i = 0; i < kMCHists; ++i) | |
358 | { | |
cfc19dd5 | 359 | fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName())); |
360 | fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName())); | |
361 | fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName())); | |
362 | if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i]) | |
0a173978 | 363 | success = kFALSE; |
364 | } | |
365 | ||
366 | for (Int_t i = 0; i < kCorrHists; ++i) | |
367 | { | |
368 | fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName())); | |
369 | if (!fCorrelation[i]) | |
370 | success = kFALSE; | |
371 | fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName())); | |
372 | if (!fMultiplicityESDCorrected[i]) | |
373 | success = kFALSE; | |
374 | } | |
375 | ||
376 | gDirectory->cd(".."); | |
377 | ||
144ff489 | 378 | // delete old hists |
379 | oldObjects.Delete(); | |
380 | ||
0a173978 | 381 | return success; |
382 | } | |
383 | ||
384 | //____________________________________________________________________ | |
144ff489 | 385 | void AliMultiplicityCorrection::SaveHistograms(const char* dir) |
0a173978 | 386 | { |
387 | // | |
388 | // saves the histograms | |
389 | // | |
390 | ||
144ff489 | 391 | if (!dir) |
392 | dir = GetName(); | |
393 | ||
394 | gDirectory->mkdir(dir); | |
395 | gDirectory->cd(dir); | |
0a173978 | 396 | |
397 | for (Int_t i = 0; i < kESDHists; ++i) | |
398 | if (fMultiplicityESD[i]) | |
399 | fMultiplicityESD[i]->Write(); | |
400 | ||
401 | for (Int_t i = 0; i < kMCHists; ++i) | |
cfc19dd5 | 402 | { |
403 | if (fMultiplicityVtx[i]) | |
404 | fMultiplicityVtx[i]->Write(); | |
405 | if (fMultiplicityMB[i]) | |
406 | fMultiplicityMB[i]->Write(); | |
407 | if (fMultiplicityINEL[i]) | |
408 | fMultiplicityINEL[i]->Write(); | |
409 | } | |
0a173978 | 410 | |
411 | for (Int_t i = 0; i < kCorrHists; ++i) | |
412 | { | |
413 | if (fCorrelation[i]) | |
414 | fCorrelation[i]->Write(); | |
415 | if (fMultiplicityESDCorrected[i]) | |
416 | fMultiplicityESDCorrected[i]->Write(); | |
417 | } | |
418 | ||
419 | gDirectory->cd(".."); | |
420 | } | |
421 | ||
422 | //____________________________________________________________________ | |
cfc19dd5 | 423 | void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll) |
0a173978 | 424 | { |
425 | // | |
426 | // Fills an event from MC | |
427 | // | |
428 | ||
cfc19dd5 | 429 | if (triggered) |
430 | { | |
431 | fMultiplicityMB[0]->Fill(vtx, generated05); | |
432 | fMultiplicityMB[1]->Fill(vtx, generated10); | |
433 | fMultiplicityMB[2]->Fill(vtx, generated15); | |
434 | fMultiplicityMB[3]->Fill(vtx, generated20); | |
435 | fMultiplicityMB[4]->Fill(vtx, generatedAll); | |
436 | ||
437 | if (vertex) | |
438 | { | |
439 | fMultiplicityVtx[0]->Fill(vtx, generated05); | |
440 | fMultiplicityVtx[1]->Fill(vtx, generated10); | |
441 | fMultiplicityVtx[2]->Fill(vtx, generated15); | |
442 | fMultiplicityVtx[3]->Fill(vtx, generated20); | |
443 | fMultiplicityVtx[4]->Fill(vtx, generatedAll); | |
444 | } | |
445 | } | |
446 | ||
447 | fMultiplicityINEL[0]->Fill(vtx, generated05); | |
448 | fMultiplicityINEL[1]->Fill(vtx, generated10); | |
449 | fMultiplicityINEL[2]->Fill(vtx, generated15); | |
450 | fMultiplicityINEL[3]->Fill(vtx, generated20); | |
451 | fMultiplicityINEL[4]->Fill(vtx, generatedAll); | |
0a173978 | 452 | } |
453 | ||
454 | //____________________________________________________________________ | |
455 | void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20) | |
456 | { | |
457 | // | |
458 | // Fills an event from ESD | |
459 | // | |
460 | ||
461 | fMultiplicityESD[0]->Fill(vtx, measured05); | |
462 | fMultiplicityESD[1]->Fill(vtx, measured10); | |
463 | fMultiplicityESD[2]->Fill(vtx, measured15); | |
464 | fMultiplicityESD[3]->Fill(vtx, measured20); | |
465 | } | |
466 | ||
467 | //____________________________________________________________________ | |
468 | void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20) | |
469 | { | |
470 | // | |
471 | // Fills an event into the correlation map with the information from MC and ESD | |
472 | // | |
473 | ||
474 | fCorrelation[0]->Fill(vtx, generated05, measured05); | |
475 | fCorrelation[1]->Fill(vtx, generated10, measured10); | |
476 | fCorrelation[2]->Fill(vtx, generated15, measured15); | |
477 | fCorrelation[3]->Fill(vtx, generated20, measured20); | |
478 | ||
479 | fCorrelation[4]->Fill(vtx, generatedAll, measured05); | |
480 | fCorrelation[5]->Fill(vtx, generatedAll, measured10); | |
481 | fCorrelation[6]->Fill(vtx, generatedAll, measured15); | |
482 | fCorrelation[7]->Fill(vtx, generatedAll, measured20); | |
483 | } | |
484 | ||
485 | //____________________________________________________________________ | |
447c325d | 486 | Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params) |
0a173978 | 487 | { |
488 | // homogenity term for minuit fitting | |
489 | // pure function of the parameters | |
490 | // prefers constant function (pol0) | |
491 | ||
492 | Double_t chi2 = 0; | |
493 | ||
447c325d | 494 | // ignore the first bin here. on purpose... |
44466df2 | 495 | for (Int_t i=2; i<fgNParamsMinuit; ++i) |
0a173978 | 496 | { |
dd701109 | 497 | Double_t right = params[i]; |
498 | Double_t left = params[i-1]; | |
0a173978 | 499 | |
447c325d | 500 | if (right != 0) |
501 | { | |
502 | Double_t diff = 1 - left / right; | |
503 | chi2 += diff * diff; | |
504 | } | |
0a173978 | 505 | } |
506 | ||
447c325d | 507 | return chi2 / 100.0; |
0a173978 | 508 | } |
509 | ||
510 | //____________________________________________________________________ | |
447c325d | 511 | Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params) |
0a173978 | 512 | { |
513 | // homogenity term for minuit fitting | |
514 | // pure function of the parameters | |
515 | // prefers linear function (pol1) | |
516 | ||
517 | Double_t chi2 = 0; | |
518 | ||
447c325d | 519 | // ignore the first bin here. on purpose... |
44466df2 | 520 | for (Int_t i=2+1; i<fgNParamsMinuit; ++i) |
0a173978 | 521 | { |
dd701109 | 522 | if (params[i-1] == 0) |
0a173978 | 523 | continue; |
524 | ||
dd701109 | 525 | Double_t right = params[i]; |
526 | Double_t middle = params[i-1]; | |
527 | Double_t left = params[i-2]; | |
528 | ||
529 | Double_t der1 = (right - middle); | |
530 | Double_t der2 = (middle - left); | |
531 | ||
532 | Double_t diff = (der1 - der2) / middle; | |
533 | ||
534 | chi2 += diff * diff; | |
535 | } | |
0a173978 | 536 | |
dd701109 | 537 | return chi2; |
538 | } | |
0a173978 | 539 | |
dd701109 | 540 | //____________________________________________________________________ |
0b4bfd98 | 541 | Double_t AliMultiplicityCorrection::RegularizationLog(TVectorD& params) |
dd701109 | 542 | { |
543 | // homogenity term for minuit fitting | |
544 | // pure function of the parameters | |
545 | // prefers linear function (pol1) | |
0a173978 | 546 | |
dd701109 | 547 | Double_t chi2 = 0; |
0a173978 | 548 | |
6d81c2de | 549 | /*Float_t der[fgkMaxParams]; |
dd701109 | 550 | |
6d81c2de | 551 | for (Int_t i=0; i<fgkMaxParams-1; ++i) |
dd701109 | 552 | der[i] = params[i+1] - params[i]; |
553 | ||
6d81c2de | 554 | for (Int_t i=0; i<fgkMaxParams-6; ++i) |
dd701109 | 555 | { |
556 | for (Int_t j=1; j<=5; ++j) | |
557 | { | |
558 | Double_t diff = der[i+j] - der[i]; | |
559 | chi2 += diff * diff; | |
560 | } | |
0b4bfd98 | 561 | }*/ |
562 | ||
563 | // ignore the first bin here. on purpose... | |
44466df2 | 564 | for (Int_t i=2+1; i<fgNParamsMinuit; ++i) |
0b4bfd98 | 565 | { |
566 | if (params[i-1] == 0) | |
567 | continue; | |
568 | ||
569 | Double_t right = log(params[i]); | |
570 | Double_t middle = log(params[i-1]); | |
571 | Double_t left = log(params[i-2]); | |
572 | ||
573 | Double_t der1 = (right - middle); | |
574 | Double_t der2 = (middle - left); | |
575 | ||
576 | Double_t diff = (der1 - der2) / middle; | |
577 | ||
578 | chi2 += diff * diff; | |
0a173978 | 579 | } |
580 | ||
0b4bfd98 | 581 | return chi2 * 100; |
0a173978 | 582 | } |
583 | ||
9ca6be09 | 584 | //____________________________________________________________________ |
447c325d | 585 | Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& params) |
9ca6be09 | 586 | { |
587 | // homogenity term for minuit fitting | |
588 | // pure function of the parameters | |
589 | // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments, | |
590 | // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp. | |
591 | ||
592 | Double_t chi2 = 0; | |
593 | ||
447c325d | 594 | // ignore the first bin here. on purpose... |
44466df2 | 595 | for (Int_t i=3; i<fgNParamsMinuit; ++i) |
9ca6be09 | 596 | { |
dd701109 | 597 | Double_t right = params[i]; |
598 | Double_t middle = params[i-1]; | |
599 | Double_t left = params[i-2]; | |
9ca6be09 | 600 | |
dd701109 | 601 | Double_t der1 = (right - middle); |
602 | Double_t der2 = (middle - left); | |
9ca6be09 | 603 | |
447c325d | 604 | Double_t diff = (der1 - der2); |
9ca6be09 | 605 | |
447c325d | 606 | chi2 += diff * diff; |
9ca6be09 | 607 | } |
608 | ||
447c325d | 609 | return chi2 * 1e4; |
3602328d | 610 | } |
611 | ||
612 | //____________________________________________________________________ | |
447c325d | 613 | Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params) |
3602328d | 614 | { |
615 | // homogenity term for minuit fitting | |
616 | // pure function of the parameters | |
617 | // calculates entropy, from | |
618 | // The method of reduced cross-entropy (M. Schmelling 1993) | |
619 | ||
3602328d | 620 | Double_t paramSum = 0; |
447c325d | 621 | // ignore the first bin here. on purpose... |
44466df2 | 622 | for (Int_t i=1; i<fgNParamsMinuit; ++i) |
dd701109 | 623 | paramSum += params[i]; |
3602328d | 624 | |
625 | Double_t chi2 = 0; | |
44466df2 | 626 | for (Int_t i=1; i<fgNParamsMinuit; ++i) |
3602328d | 627 | { |
0b4bfd98 | 628 | //Double_t tmp = params[i] / paramSum; |
629 | Double_t tmp = params[i]; | |
6d81c2de | 630 | if (tmp > 0 && (*fgEntropyAPriori)[i] > 0) |
631 | chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]); | |
3602328d | 632 | } |
3602328d | 633 | |
0f67a57c | 634 | return 100.0 + chi2; |
dd701109 | 635 | } |
636 | ||
637 | //____________________________________________________________________ | |
638 | void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3) | |
639 | { | |
640 | // | |
641 | // fit function for minuit | |
642 | // does: nbd | |
643 | // func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50); | |
644 | // func->SetParNames("scaling", "averagen", "k"); | |
447c325d | 645 | // func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000); |
646 | // func->SetParLimits(1, 0.001, 1000); | |
647 | // func->SetParLimits(2, 0.001, 1000); | |
dd701109 | 648 | // |
649 | ||
6d81c2de | 650 | fgNBD->SetParameters(params[0], params[1], params[2]); |
dd701109 | 651 | |
6d81c2de | 652 | Double_t params2[fgkMaxParams]; |
dd701109 | 653 | |
6d81c2de | 654 | for (Int_t i=0; i<fgkMaxParams; ++i) |
655 | params2[i] = fgNBD->Eval(i); | |
dd701109 | 656 | |
657 | MinuitFitFunction(unused1, unused2, chi2, params2, unused3); | |
658 | ||
659 | printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2); | |
9ca6be09 | 660 | } |
661 | ||
0a173978 | 662 | //____________________________________________________________________ |
663 | void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t) | |
664 | { | |
665 | // | |
666 | // fit function for minuit | |
cfc19dd5 | 667 | // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix |
0a173978 | 668 | // |
669 | ||
cfc19dd5 | 670 | // d |
44466df2 | 671 | static TVectorD paramsVector(fgNParamsMinuit); |
672 | for (Int_t i=0; i<fgNParamsMinuit; ++i) | |
447c325d | 673 | paramsVector[i] = params[i] * params[i]; |
674 | ||
675 | // calculate penalty factor | |
676 | Double_t penaltyVal = 0; | |
6d81c2de | 677 | switch (fgRegularizationType) |
447c325d | 678 | { |
679 | case kNone: break; | |
680 | case kPol0: penaltyVal = RegularizationPol0(paramsVector); break; | |
681 | case kPol1: penaltyVal = RegularizationPol1(paramsVector); break; | |
682 | case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break; | |
683 | case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break; | |
0b4bfd98 | 684 | case kLog: penaltyVal = RegularizationLog(paramsVector); break; |
447c325d | 685 | } |
686 | ||
687 | //if (penaltyVal > 1e10) | |
688 | // paramsVector2.Print(); | |
689 | ||
6d81c2de | 690 | penaltyVal *= fgRegularizationWeight; |
cfc19dd5 | 691 | |
692 | // Ad | |
6d81c2de | 693 | TVectorD measGuessVector(fgkMaxInput); |
694 | measGuessVector = (*fgCorrelationMatrix) * paramsVector; | |
cfc19dd5 | 695 | |
696 | // Ad - m | |
6d81c2de | 697 | measGuessVector -= (*fgCurrentESDVector); |
cfc19dd5 | 698 | |
44466df2 | 699 | //measGuessVector.Print(); |
700 | ||
447c325d | 701 | TVectorD copy(measGuessVector); |
cfc19dd5 | 702 | |
703 | // (Ad - m) W | |
6d81c2de | 704 | // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used |
0b4bfd98 | 705 | // normal way is like this: |
6d81c2de | 706 | // measGuessVector *= (*fgCorrelationCovarianceMatrix); |
0b4bfd98 | 707 | // optimized way like this: |
0f67a57c | 708 | for (Int_t i=0; i<fgkMaxInput; ++i) |
6d81c2de | 709 | measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i); |
cfc19dd5 | 710 | |
447c325d | 711 | //measGuessVector.Print(); |
cfc19dd5 | 712 | |
0f67a57c | 713 | // reduce influence of first bin |
714 | copy(0) *= 0.01; | |
715 | ||
447c325d | 716 | // (Ad - m) W (Ad - m) |
0b4bfd98 | 717 | // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very |
6d81c2de | 718 | // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit) |
447c325d | 719 | Double_t chi2FromFit = measGuessVector * copy * 1e6; |
0a173978 | 720 | |
3602328d | 721 | chi2 = chi2FromFit + penaltyVal; |
0a173978 | 722 | |
447c325d | 723 | static Int_t callCount = 0; |
724 | if ((callCount++ % 10000) == 0) | |
725 | printf("%d) %f %f --> %f\n", callCount, chi2FromFit, penaltyVal, chi2); | |
0a173978 | 726 | } |
727 | ||
0f67a57c | 728 | //____________________________________________________________________ |
729 | void AliMultiplicityCorrection::DrawGuess(Double_t *params) | |
730 | { | |
731 | // | |
732 | // draws residuals of solution suggested by params and effect of regularization | |
733 | // | |
734 | ||
735 | // d | |
736 | static TVectorD paramsVector(fgNParamsMinuit); | |
737 | for (Int_t i=0; i<fgNParamsMinuit; ++i) | |
738 | paramsVector[i] = params[i] * params[i]; | |
739 | ||
740 | // Ad | |
741 | TVectorD measGuessVector(fgkMaxInput); | |
742 | measGuessVector = (*fgCorrelationMatrix) * paramsVector; | |
743 | ||
744 | // Ad - m | |
745 | measGuessVector -= (*fgCurrentESDVector); | |
746 | ||
747 | TVectorD copy(measGuessVector); | |
748 | ||
749 | // (Ad - m) W | |
750 | // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used | |
751 | // normal way is like this: | |
752 | // measGuessVector *= (*fgCorrelationCovarianceMatrix); | |
753 | // optimized way like this: | |
754 | for (Int_t i=0; i<fgkMaxInput; ++i) | |
755 | measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i); | |
756 | ||
757 | // (Ad - m) W (Ad - m) | |
758 | // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very | |
759 | // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit) | |
760 | //Double_t chi2FromFit = measGuessVector * copy * 1e6; | |
761 | ||
762 | // draw residuals | |
763 | TH1* residuals = new TH1F("residuals", "residuals", fgkMaxInput, -0.5, fgkMaxInput - 0.5); | |
764 | for (Int_t i=0; i<fgkMaxInput; ++i) | |
765 | residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6); | |
766 | new TCanvas; residuals->DrawCopy(); gPad->SetLogy(); | |
767 | ||
768 | // draw penalty | |
769 | if (fgRegularizationType != kPol1) { | |
770 | Printf("Drawing not supported"); | |
771 | return; | |
772 | } | |
773 | ||
774 | TH1* penalty = new TH1F("penalty", "penalty", fgNParamsMinuit, -0.5, fgNParamsMinuit - 0.5); | |
775 | ||
776 | for (Int_t i=2; i<fgNParamsMinuit; ++i) | |
777 | { | |
778 | if (params[i-1] == 0) | |
779 | continue; | |
780 | ||
781 | Double_t right = params[i]; | |
782 | Double_t middle = params[i-1]; | |
783 | Double_t left = params[i-2]; | |
784 | ||
785 | Double_t der1 = (right - middle); | |
786 | Double_t der2 = (middle - left); | |
787 | ||
788 | Double_t diff = (der1 - der2) / middle; | |
789 | ||
790 | penalty->SetBinContent(i, diff * diff); | |
791 | } | |
792 | new TCanvas; penalty->DrawCopy(); gPad->SetLogy(); | |
793 | } | |
794 | ||
0a173978 | 795 | //____________________________________________________________________ |
447c325d | 796 | void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin) |
0a173978 | 797 | { |
798 | // | |
9ca6be09 | 799 | // fills fCurrentESD, fCurrentCorrelation |
800 | // resets fMultiplicityESDCorrected | |
0a173978 | 801 | // |
802 | ||
0b4bfd98 | 803 | Bool_t display = kFALSE; |
804 | ||
0a173978 | 805 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); |
0b4bfd98 | 806 | |
9ca6be09 | 807 | fMultiplicityESDCorrected[correlationID]->Reset(); |
0b4bfd98 | 808 | fMultiplicityESDCorrected[correlationID]->Sumw2(); |
0a173978 | 809 | |
039e265e | 810 | // project without under/overflow bins |
811 | fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD", 1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins()); | |
9ca6be09 | 812 | fCurrentESD->Sumw2(); |
0a173978 | 813 | |
814 | // empty under/overflow bins in x, otherwise Project3D takes them into account | |
815 | TH3* hist = fCorrelation[correlationID]; | |
039e265e | 816 | for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y) |
0a173978 | 817 | { |
039e265e | 818 | for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z) |
0a173978 | 819 | { |
820 | hist->SetBinContent(0, y, z, 0); | |
821 | hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); | |
822 | } | |
823 | } | |
824 | ||
9ca6be09 | 825 | fCurrentCorrelation = hist->Project3D("zy"); |
447c325d | 826 | //((TH2*) fCurrentCorrelation)->Rebin2D(2, 1); |
827 | //fMultiplicityESDCorrected[correlationID]->Rebin(2); | |
9ca6be09 | 828 | fCurrentCorrelation->Sumw2(); |
cfc19dd5 | 829 | |
44466df2 | 830 | Printf("AliMultiplicityCorrection::SetupCurrentHists: Statistics information: %.f entries in correlation map; %.f entries in measured spectrum", fCurrentCorrelation->Integral(), fCurrentESD->Integral()); |
831 | ||
447c325d | 832 | if (createBigBin) |
833 | { | |
0f67a57c | 834 | fLastBinLimit = 0; |
447c325d | 835 | for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i) |
836 | { | |
837 | if (fCurrentESD->GetBinContent(i) <= 5) | |
838 | { | |
0f67a57c | 839 | fLastBinLimit = i; |
447c325d | 840 | break; |
841 | } | |
842 | } | |
843 | ||
0f67a57c | 844 | if (fLastBinLimit > 0) |
447c325d | 845 | { |
0b4bfd98 | 846 | TCanvas* canvas = 0; |
447c325d | 847 | |
0b4bfd98 | 848 | if (display) |
849 | { | |
850 | canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800); | |
851 | canvas->Divide(2, 2); | |
852 | ||
853 | canvas->cd(1); | |
854 | fCurrentESD->GetXaxis()->SetRangeUser(0, 200); | |
855 | fCurrentESD->SetStats(kFALSE); | |
856 | fCurrentESD->SetTitle(";measured multiplicity"); | |
857 | fCurrentESD->DrawCopy(); | |
858 | gPad->SetLogy(); | |
859 | ||
860 | canvas->cd(2); | |
861 | fCurrentCorrelation->GetXaxis()->SetRangeUser(0, 250); | |
862 | fCurrentCorrelation->GetYaxis()->SetRangeUser(0, 250); | |
863 | fCurrentCorrelation->SetTitle(";true multiplicity;measured multiplicity"); | |
864 | fCurrentCorrelation->SetStats(kFALSE); | |
865 | fCurrentCorrelation->DrawCopy("COLZ"); | |
866 | } | |
447c325d | 867 | |
0f67a57c | 868 | printf("Bin limit in measured spectrum is %d.\n", fLastBinLimit); |
869 | fCurrentESD->SetBinContent(fLastBinLimit, fCurrentESD->Integral(fLastBinLimit, fCurrentESD->GetNbinsX())); | |
870 | for (Int_t i=fLastBinLimit+1; i<=fCurrentESD->GetNbinsX(); ++i) | |
447c325d | 871 | { |
872 | fCurrentESD->SetBinContent(i, 0); | |
873 | fCurrentESD->SetBinError(i, 0); | |
874 | } | |
875 | // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? | |
0f67a57c | 876 | fCurrentESD->SetBinError(fLastBinLimit, TMath::Sqrt(fCurrentESD->GetBinContent(fLastBinLimit))); |
447c325d | 877 | |
0f67a57c | 878 | printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(fLastBinLimit), fCurrentESD->GetBinError(fLastBinLimit)); |
447c325d | 879 | |
880 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
881 | { | |
0f67a57c | 882 | fCurrentCorrelation->SetBinContent(i, fLastBinLimit, fCurrentCorrelation->Integral(i, i, fLastBinLimit, fCurrentCorrelation->GetNbinsY())); |
447c325d | 883 | // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? |
0f67a57c | 884 | fCurrentCorrelation->SetBinError(i, fLastBinLimit, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, fLastBinLimit))); |
447c325d | 885 | |
0f67a57c | 886 | for (Int_t j=fLastBinLimit+1; j<=fCurrentCorrelation->GetNbinsY(); ++j) |
447c325d | 887 | { |
888 | fCurrentCorrelation->SetBinContent(i, j, 0); | |
889 | fCurrentCorrelation->SetBinError(i, j, 0); | |
890 | } | |
891 | } | |
892 | ||
893 | printf("Adjusted correlation matrix!\n"); | |
894 | ||
0b4bfd98 | 895 | if (canvas) |
896 | { | |
897 | canvas->cd(3); | |
898 | fCurrentESD->DrawCopy(); | |
899 | gPad->SetLogy(); | |
447c325d | 900 | |
0b4bfd98 | 901 | canvas->cd(4); |
902 | fCurrentCorrelation->DrawCopy("COLZ"); | |
903 | } | |
447c325d | 904 | } |
905 | } | |
906 | ||
0b4bfd98 | 907 | #if 0 // does not help |
908 | // null bins with one entry | |
909 | Int_t nNulledBins = 0; | |
910 | for (Int_t x=1; x<=fCurrentCorrelation->GetXaxis()->GetNbins(); ++x) | |
911 | for (Int_t y=1; y<=fCurrentCorrelation->GetYaxis()->GetNbins(); ++y) | |
912 | { | |
913 | if (fCurrentCorrelation->GetBinContent(x, y) == 1) | |
914 | { | |
915 | fCurrentCorrelation->SetBinContent(x, y, 0); | |
916 | fCurrentCorrelation->SetBinError(x, y, 0); | |
447c325d | 917 | |
0b4bfd98 | 918 | ++nNulledBins; |
919 | } | |
920 | } | |
921 | Printf("Nulled %d bins", nNulledBins); | |
922 | #endif | |
923 | ||
924 | fCurrentEfficiency = GetEfficiency(inputRange, eventType); | |
925 | //fCurrentEfficiency->Rebin(2); | |
926 | //fCurrentEfficiency->Scale(0.5); | |
927 | } | |
928 | ||
929 | //____________________________________________________________________ | |
930 | TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventType) | |
931 | { | |
932 | // | |
933 | // calculates efficiency for given event type | |
934 | // | |
935 | ||
936 | TH1* divisor = 0; | |
cfc19dd5 | 937 | switch (eventType) |
938 | { | |
0b4bfd98 | 939 | case kTrVtx : divisor = fMultiplicityVtx[inputRange]->ProjectionY("divisor", -1, -1, "e"); break; |
940 | case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", -1, -1, "e"); break; | |
941 | case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", -1, -1, "e"); break; | |
cfc19dd5 | 942 | } |
0b4bfd98 | 943 | TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", -1, -1, "e"); |
944 | eff->Divide(eff, divisor, 1, 1, "B"); | |
945 | return eff; | |
9ca6be09 | 946 | } |
947 | ||
948 | //____________________________________________________________________ | |
44466df2 | 949 | void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight, Int_t minuitParams) |
447c325d | 950 | { |
6d81c2de | 951 | // |
952 | // sets the parameters for chi2 minimization | |
953 | // | |
954 | ||
955 | fgRegularizationType = type; | |
956 | fgRegularizationWeight = weight; | |
447c325d | 957 | |
958 | printf("AliMultiplicityCorrection::SetRegularizationParameters --> Regularization set to %d with weight %f\n", (Int_t) type, weight); | |
44466df2 | 959 | |
960 | if (minuitParams > 0) | |
961 | { | |
962 | fgNParamsMinuit = minuitParams; | |
963 | printf("AliMultiplicityCorrection::SetRegularizationParameters --> Number of MINUIT iterations set to %d\n", minuitParams); | |
964 | } | |
447c325d | 965 | } |
966 | ||
967 | //____________________________________________________________________ | |
6d81c2de | 968 | void AliMultiplicityCorrection::SetBayesianParameters(Float_t smoothing, Int_t nIterations) |
9ca6be09 | 969 | { |
970 | // | |
6d81c2de | 971 | // sets the parameters for Bayesian unfolding |
9ca6be09 | 972 | // |
973 | ||
6d81c2de | 974 | fgBayesianSmoothing = smoothing; |
975 | fgBayesianIterations = nIterations; | |
9ca6be09 | 976 | |
6d81c2de | 977 | printf("AliMultiplicityCorrection::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing); |
978 | } | |
cfc19dd5 | 979 | |
6d81c2de | 980 | //____________________________________________________________________ |
981 | Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check) | |
982 | { | |
983 | // | |
984 | // implemenation of unfolding (internal function) | |
985 | // | |
986 | // unfolds <measured> using response from <correlation> and effiency <aEfficiency> | |
987 | // output is in <result> | |
988 | // <initialConditions> set the initial values for the minimization, if 0 <measured> is used | |
989 | // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed | |
990 | // | |
991 | // returns minuit status (0 = success), (-1 when check was set) | |
992 | // | |
447c325d | 993 | |
6d81c2de | 994 | //normalize measured |
995 | measured->Scale(1.0 / measured->Integral()); | |
996 | ||
997 | if (!fgCorrelationMatrix) | |
44466df2 | 998 | fgCorrelationMatrix = new TMatrixD(fgkMaxInput, fgNParamsMinuit); |
6d81c2de | 999 | if (!fgCorrelationCovarianceMatrix) |
1000 | fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxInput, fgkMaxInput); | |
1001 | if (!fgCurrentESDVector) | |
1002 | fgCurrentESDVector = new TVectorD(fgkMaxInput); | |
1003 | if (!fgEntropyAPriori) | |
44466df2 | 1004 | fgEntropyAPriori = new TVectorD(fgNParamsMinuit); |
6d81c2de | 1005 | |
1006 | fgCorrelationMatrix->Zero(); | |
1007 | fgCorrelationCovarianceMatrix->Zero(); | |
1008 | fgCurrentESDVector->Zero(); | |
1009 | fgEntropyAPriori->Zero(); | |
0a173978 | 1010 | |
1011 | // normalize correction for given nPart | |
6d81c2de | 1012 | for (Int_t i=1; i<=correlation->GetNbinsX(); ++i) |
0a173978 | 1013 | { |
6d81c2de | 1014 | Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY()); |
0a173978 | 1015 | if (sum <= 0) |
1016 | continue; | |
9ca6be09 | 1017 | Float_t maxValue = 0; |
1018 | Int_t maxBin = -1; | |
6d81c2de | 1019 | for (Int_t j=1; j<=correlation->GetNbinsY(); ++j) |
0a173978 | 1020 | { |
9ca6be09 | 1021 | // find most probably value |
6d81c2de | 1022 | if (maxValue < correlation->GetBinContent(i, j)) |
9ca6be09 | 1023 | { |
6d81c2de | 1024 | maxValue = correlation->GetBinContent(i, j); |
9ca6be09 | 1025 | maxBin = j; |
1026 | } | |
1027 | ||
0a173978 | 1028 | // npart sum to 1 |
6d81c2de | 1029 | correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum); |
1030 | correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum); | |
cfc19dd5 | 1031 | |
44466df2 | 1032 | if (i <= fgNParamsMinuit && j <= fgkMaxInput) |
6d81c2de | 1033 | (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j); |
0a173978 | 1034 | } |
9ca6be09 | 1035 | |
cfc19dd5 | 1036 | //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin)); |
0a173978 | 1037 | } |
1038 | ||
0a173978 | 1039 | // Initialize TMinuit via generic fitter interface |
44466df2 | 1040 | TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgNParamsMinuit); |
447c325d | 1041 | Double_t arglist[100]; |
0b4bfd98 | 1042 | |
6d81c2de | 1043 | // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way... |
447c325d | 1044 | arglist[0] = 0; |
1045 | minuit->ExecuteCommand("SET PRINT", arglist, 1); | |
0a173978 | 1046 | |
0b4bfd98 | 1047 | // however, enable warnings |
1048 | //minuit->ExecuteCommand("SET WAR", arglist, 0); | |
1049 | ||
1050 | // set minimization function | |
0a173978 | 1051 | minuit->SetFCN(MinuitFitFunction); |
1052 | ||
44466df2 | 1053 | for (Int_t i=0; i<fgNParamsMinuit; i++) |
6d81c2de | 1054 | (*fgEntropyAPriori)[i] = 1; |
0a173978 | 1055 | |
0f67a57c | 1056 | if (!initialConditions) { |
1057 | initialConditions = measured; | |
1058 | } else { | |
dd701109 | 1059 | printf("Using different starting conditions...\n"); |
6d81c2de | 1060 | //new TCanvas; initialConditions->DrawCopy(); |
6d81c2de | 1061 | initialConditions->Scale(1.0 / initialConditions->Integral()); |
dd701109 | 1062 | } |
0f67a57c | 1063 | |
1064 | for (Int_t i=0; i<fgNParamsMinuit; i++) | |
1065 | if (initialConditions->GetBinContent(i+1) > 0) | |
1066 | (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1); | |
447c325d | 1067 | |
44466df2 | 1068 | Double_t* results = new Double_t[fgNParamsMinuit+1]; |
1069 | for (Int_t i=0; i<fgNParamsMinuit; ++i) | |
0a173978 | 1070 | { |
6d81c2de | 1071 | results[i] = initialConditions->GetBinContent(i+1); |
dd701109 | 1072 | |
447c325d | 1073 | // minimum value |
44466df2 | 1074 | if (!check) |
1075 | results[i] = TMath::Max(results[i], 1e-3); | |
447c325d | 1076 | |
1077 | Float_t stepSize = 0.1; | |
1078 | ||
1079 | // minuit sees squared values to prevent it from going negative... | |
1080 | results[i] = TMath::Sqrt(results[i]); | |
1081 | //stepSize /= results[i]; // keep relative step size | |
1082 | ||
1083 | minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0); | |
1084 | } | |
6d81c2de | 1085 | //minuit->SetParameter(fgkMaxParams, "alpha", 1e4, 1, 0, 0); |
447c325d | 1086 | // bin 0 is filled with value from bin 1 (otherwise it's 0) |
1087 | //minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 0); | |
1088 | //results[0] = 0; | |
1089 | //minuit->SetParameter(0, "param0", 0, 0, 0, 0); | |
1090 | ||
6d81c2de | 1091 | for (Int_t i=0; i<fgkMaxInput; ++i) |
447c325d | 1092 | { |
6d81c2de | 1093 | (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1); |
1094 | if (measured->GetBinError(i+1) > 0) | |
1095 | (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1); | |
cfc19dd5 | 1096 | |
6d81c2de | 1097 | if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7) |
1098 | (*fgCorrelationCovarianceMatrix)(i, i) = 0; | |
0a173978 | 1099 | } |
1100 | ||
0f67a57c | 1101 | Int_t dummy = 0; |
1102 | Double_t chi2 = 0; | |
0a173978 | 1103 | MinuitFitFunction(dummy, 0, chi2, results, 0); |
0f67a57c | 1104 | DrawGuess(results); |
dd701109 | 1105 | printf("Chi2 of initial parameters is = %f\n", chi2); |
0a173978 | 1106 | |
9ca6be09 | 1107 | if (check) |
0f67a57c | 1108 | { |
1109 | delete[] results; | |
447c325d | 1110 | return -1; |
0f67a57c | 1111 | } |
447c325d | 1112 | |
1113 | // first param is number of iterations, second is precision.... | |
1114 | arglist[0] = 1e6; | |
1115 | //arglist[1] = 1e-5; | |
1116 | //minuit->ExecuteCommand("SCAN", arglist, 0); | |
1117 | Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1); | |
1118 | printf("MINUIT status is %d\n", status); | |
1119 | //minuit->ExecuteCommand("MIGRAD", arglist, 1); | |
1120 | //minuit->ExecuteCommand("MIGRAD", arglist, 1); | |
dd701109 | 1121 | //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!"); |
0a173978 | 1122 | //minuit->ExecuteCommand("MINOS", arglist, 0); |
1123 | ||
44466df2 | 1124 | for (Int_t i=0; i<fgNParamsMinuit; ++i) |
0a173978 | 1125 | { |
6d81c2de | 1126 | if (aEfficiency->GetBinContent(i+1) > 0) |
dd701109 | 1127 | { |
0f67a57c | 1128 | results[i] = minuit->GetParameter(i); |
1129 | result->SetBinContent(i+1, results[i] * results[i] / aEfficiency->GetBinContent(i+1)); | |
447c325d | 1130 | // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i)) |
0f67a57c | 1131 | result->SetBinError(i+1, minuit->GetParError(i) * results[i] / aEfficiency->GetBinContent(i+1)); |
dd701109 | 1132 | } |
1133 | else | |
1134 | { | |
6d81c2de | 1135 | result->SetBinContent(i+1, 0); |
1136 | result->SetBinError(i+1, 0); | |
dd701109 | 1137 | } |
1138 | } | |
0f67a57c | 1139 | DrawGuess(results); |
1140 | ||
1141 | delete[] results; | |
447c325d | 1142 | |
1143 | return status; | |
dd701109 | 1144 | } |
1145 | ||
6d81c2de | 1146 | //____________________________________________________________________ |
1147 | Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions) | |
1148 | { | |
1149 | // | |
1150 | // correct spectrum using minuit chi2 method | |
1151 | // | |
1152 | // for description of parameters, see UnfoldWithMinuit | |
1153 | // | |
1154 | ||
1155 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
44466df2 | 1156 | SetupCurrentHists(inputRange, fullPhaseSpace, eventType, fgCreateBigBin); |
6d81c2de | 1157 | |
1158 | return UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check); | |
1159 | } | |
1160 | ||
dd701109 | 1161 | //____________________________________________________________________ |
1162 | void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace) | |
1163 | { | |
1164 | // | |
1165 | // correct spectrum using minuit chi2 method applying a NBD fit | |
1166 | // | |
1167 | ||
1168 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
1169 | ||
447c325d | 1170 | SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE); |
0b4bfd98 | 1171 | //normalize ESD |
1172 | fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); | |
dd701109 | 1173 | |
6d81c2de | 1174 | fgCorrelationMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams); |
1175 | fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams); | |
1176 | fgCurrentESDVector = new TVectorD(fgkMaxParams); | |
dd701109 | 1177 | |
1178 | // normalize correction for given nPart | |
1179 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
1180 | { | |
1181 | Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); | |
1182 | if (sum <= 0) | |
1183 | continue; | |
1184 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) | |
1185 | { | |
1186 | // npart sum to 1 | |
1187 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); | |
1188 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
1189 | ||
6d81c2de | 1190 | if (i <= fgkMaxParams && j <= fgkMaxParams) |
1191 | (*fgCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j); | |
dd701109 | 1192 | } |
1193 | } | |
1194 | ||
6d81c2de | 1195 | for (Int_t i=0; i<fgkMaxParams; ++i) |
dd701109 | 1196 | { |
6d81c2de | 1197 | (*fgCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1); |
dd701109 | 1198 | if (fCurrentESD->GetBinError(i+1) > 0) |
6d81c2de | 1199 | (*fgCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1)); |
dd701109 | 1200 | } |
1201 | ||
1202 | // Create NBD function | |
6d81c2de | 1203 | if (!fgNBD) |
dd701109 | 1204 | { |
6d81c2de | 1205 | fgNBD = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 250); |
1206 | fgNBD->SetParNames("scaling", "averagen", "k"); | |
0a173978 | 1207 | } |
1208 | ||
dd701109 | 1209 | // Initialize TMinuit via generic fitter interface |
1210 | TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3); | |
1211 | ||
1212 | minuit->SetFCN(MinuitNBD); | |
1213 | minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000); | |
1214 | minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000); | |
1215 | minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000); | |
1216 | ||
1217 | Double_t arglist[100]; | |
1218 | arglist[0] = 0; | |
1219 | minuit->ExecuteCommand("SET PRINT", arglist, 1); | |
447c325d | 1220 | minuit->ExecuteCommand("MIGRAD", arglist, 0); |
dd701109 | 1221 | //minuit->ExecuteCommand("MINOS", arglist, 0); |
1222 | ||
6d81c2de | 1223 | fgNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2)); |
dd701109 | 1224 | |
6d81c2de | 1225 | for (Int_t i=0; i<fgkMaxParams; ++i) |
dd701109 | 1226 | { |
6d81c2de | 1227 | printf("%d : %f\n", i, fgNBD->Eval(i)); |
1228 | if (fgNBD->Eval(i) > 0) | |
dd701109 | 1229 | { |
6d81c2de | 1230 | fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fgNBD->Eval(i)); |
dd701109 | 1231 | fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0); |
1232 | } | |
1233 | } | |
0a173978 | 1234 | } |
1235 | ||
0a173978 | 1236 | //____________________________________________________________________ |
1237 | void AliMultiplicityCorrection::DrawHistograms() | |
1238 | { | |
6d81c2de | 1239 | // |
1240 | // draws the histograms of this class | |
1241 | // | |
1242 | ||
0a173978 | 1243 | printf("ESD:\n"); |
1244 | ||
1245 | TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600); | |
1246 | canvas1->Divide(3, 2); | |
1247 | for (Int_t i = 0; i < kESDHists; ++i) | |
1248 | { | |
1249 | canvas1->cd(i+1); | |
1250 | fMultiplicityESD[i]->DrawCopy("COLZ"); | |
1251 | printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean()); | |
1252 | } | |
1253 | ||
cfc19dd5 | 1254 | printf("Vtx:\n"); |
0a173978 | 1255 | |
1256 | TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600); | |
1257 | canvas2->Divide(3, 2); | |
1258 | for (Int_t i = 0; i < kMCHists; ++i) | |
1259 | { | |
1260 | canvas2->cd(i+1); | |
cfc19dd5 | 1261 | fMultiplicityVtx[i]->DrawCopy("COLZ"); |
1262 | printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean()); | |
0a173978 | 1263 | } |
1264 | ||
1265 | TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900); | |
1266 | canvas3->Divide(3, 3); | |
1267 | for (Int_t i = 0; i < kCorrHists; ++i) | |
1268 | { | |
1269 | canvas3->cd(i+1); | |
b477f4f2 | 1270 | TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone()); |
1271 | for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y) | |
1272 | { | |
1273 | for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z) | |
1274 | { | |
1275 | hist->SetBinContent(0, y, z, 0); | |
1276 | hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); | |
1277 | } | |
1278 | } | |
1279 | TH1* proj = hist->Project3D("zy"); | |
0a173978 | 1280 | proj->DrawCopy("COLZ"); |
1281 | } | |
1282 | } | |
1283 | ||
1284 | //____________________________________________________________________ | |
1c91c693 | 1285 | void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t /*normalizeESD*/, TH1* mcHist, Bool_t simple) |
0a173978 | 1286 | { |
039e265e | 1287 | // draw comparison plots |
1288 | ||
1289 | ||
447c325d | 1290 | //mcHist->Rebin(2); |
1291 | ||
cfc19dd5 | 1292 | Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); |
1293 | ||
9ca6be09 | 1294 | TString tmpStr; |
cfc19dd5 | 1295 | tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId); |
0a173978 | 1296 | |
447c325d | 1297 | if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0) |
1298 | { | |
1299 | printf("ERROR. Unfolded histogram is empty\n"); | |
1300 | return; | |
1301 | } | |
1302 | ||
1303 | //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists | |
039e265e | 1304 | fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY("fCurrentESD", 1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins()); |
447c325d | 1305 | fCurrentESD->Sumw2(); |
1306 | fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); | |
1307 | ||
1308 | // normalize unfolded result to 1 | |
1309 | fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral()); | |
1310 | ||
1311 | //fCurrentESD->Scale(mcHist->Integral(2, 200)); | |
1312 | ||
1313 | //new TCanvas; | |
1314 | /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio"); | |
1315 | ratio->Divide(mcHist); | |
1316 | ratio->Draw("HIST"); | |
1317 | ratio->Fit("pol0", "W0", "", 20, 120); | |
1318 | Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0); | |
1319 | delete ratio; | |
1320 | mcHist->Scale(scalingFactor);*/ | |
1321 | ||
0f67a57c | 1322 | // find bin with <= 50 entries. this is used as limit for the chi2 calculation |
0b4bfd98 | 1323 | Int_t mcBinLimit = 0; |
1324 | for (Int_t i=20; i<mcHist->GetNbinsX(); ++i) | |
1325 | { | |
1326 | if (mcHist->GetBinContent(i) > 50) | |
1327 | { | |
1328 | mcBinLimit = i; | |
1329 | } | |
1330 | else | |
1331 | break; | |
1332 | } | |
1333 | Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit); | |
1334 | ||
1335 | // scale to 1 | |
1336 | mcHist->Sumw2(); | |
144ff489 | 1337 | if (mcHist->Integral() > 0) |
1338 | mcHist->Scale(1.0 / mcHist->Integral()); | |
447c325d | 1339 | |
b477f4f2 | 1340 | // calculate residual |
1341 | ||
1342 | // for that we convolute the response matrix with the gathered result | |
1343 | // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD | |
dd701109 | 1344 | TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected"); |
1345 | tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency); | |
447c325d | 1346 | TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId); |
1347 | TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e"); | |
1348 | if (convolutedProj->Integral() > 0) | |
1349 | { | |
1350 | convolutedProj->Scale(1.0 / convolutedProj->Integral()); | |
1351 | } | |
1352 | else | |
1353 | printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n"); | |
1354 | ||
447c325d | 1355 | TH1* residual = (TH1*) convolutedProj->Clone("residual"); |
1356 | residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e"); | |
b477f4f2 | 1357 | |
1358 | residual->Add(fCurrentESD, -1); | |
447c325d | 1359 | //residual->Divide(residual, fCurrentESD, 1, 1, "B"); |
1360 | ||
0b4bfd98 | 1361 | TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5); |
b477f4f2 | 1362 | |
1363 | // TODO fix errors | |
447c325d | 1364 | Float_t chi2 = 0; |
1365 | for (Int_t i=1; i<=residual->GetNbinsX(); ++i) | |
b477f4f2 | 1366 | { |
447c325d | 1367 | if (fCurrentESD->GetBinError(i) > 0) |
1368 | { | |
1369 | Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i); | |
0f67a57c | 1370 | if (i > 1 && i <= fLastBinLimit) |
447c325d | 1371 | chi2 += value * value; |
0f67a57c | 1372 | Printf("%d --> %f (%f)", i, value * value, chi2); |
447c325d | 1373 | residual->SetBinContent(i, value); |
1374 | residualHist->Fill(value); | |
1375 | } | |
1376 | else | |
1377 | { | |
1378 | //printf("Residual bin %d set to 0\n", i); | |
1379 | residual->SetBinContent(i, 0); | |
1380 | } | |
1381 | convolutedProj->SetBinError(i, 0); | |
b477f4f2 | 1382 | residual->SetBinError(i, 0); |
447c325d | 1383 | } |
0f67a57c | 1384 | fLastChi2Residuals = chi2; |
dd701109 | 1385 | |
0b4bfd98 | 1386 | new TCanvas; residualHist->DrawCopy(); |
1387 | ||
1388 | //residualHist->Fit("gaus", "N"); | |
1389 | //delete residualHist; | |
dd701109 | 1390 | |
0f67a57c | 1391 | printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, fLastBinLimit); |
dd701109 | 1392 | |
447c325d | 1393 | TCanvas* canvas1 = 0; |
1394 | if (simple) | |
1395 | { | |
1396 | canvas1 = new TCanvas(tmpStr, tmpStr, 900, 400); | |
1397 | canvas1->Divide(2, 1); | |
1398 | } | |
1399 | else | |
1400 | { | |
1401 | canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200); | |
1402 | canvas1->Divide(2, 3); | |
1403 | } | |
0a173978 | 1404 | |
1405 | canvas1->cd(1); | |
cfc19dd5 | 1406 | TH1* proj = (TH1*) mcHist->Clone("proj"); |
9ca6be09 | 1407 | |
039e265e | 1408 | // normalize without 0 bin |
1409 | proj->Scale(1.0 / proj->Integral(2, proj->GetNbinsX())); | |
1410 | Printf("Normalized without 0 bin!"); | |
dd701109 | 1411 | proj->GetXaxis()->SetRangeUser(0, 200); |
447c325d | 1412 | proj->SetTitle(";true multiplicity;Entries"); |
1413 | proj->SetStats(kFALSE); | |
0a173978 | 1414 | |
cfc19dd5 | 1415 | fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2); |
039e265e | 1416 | // normalize without 0 bin |
1417 | fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral(2, fMultiplicityESDCorrected[esdCorrId]->GetNbinsX())); | |
1418 | Printf("Normalized without 0 bin!"); | |
0f67a57c | 1419 | |
1420 | if (proj->GetEntries() > 0) { | |
1421 | proj->DrawCopy("HIST"); | |
1422 | fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E"); | |
1423 | } | |
1424 | else | |
1425 | fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST E"); | |
1426 | ||
1427 | gPad->SetLogy(); | |
cfc19dd5 | 1428 | |
447c325d | 1429 | TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93); |
1430 | legend->AddEntry(proj, "true distribution"); | |
1431 | legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "unfolded distribution"); | |
1432 | legend->SetFillColor(0); | |
1433 | legend->Draw(); | |
1434 | // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14); | |
1435 | ||
1436 | canvas1->cd(2); | |
1437 | ||
1438 | gPad->SetLogy(); | |
1439 | fCurrentESD->GetXaxis()->SetRangeUser(0, 200); | |
1440 | //fCurrentESD->SetLineColor(2); | |
1441 | fCurrentESD->SetTitle(";measured multiplicity;Entries"); | |
1442 | fCurrentESD->SetStats(kFALSE); | |
1443 | fCurrentESD->DrawCopy("HIST E"); | |
1444 | ||
1445 | convolutedProj->SetLineColor(2); | |
1446 | //proj2->SetMarkerColor(2); | |
1447 | //proj2->SetMarkerStyle(5); | |
1448 | convolutedProj->DrawCopy("HIST SAME"); | |
1449 | ||
1450 | legend = new TLegend(0.3, 0.8, 0.93, 0.93); | |
1451 | legend->AddEntry(fCurrentESD, "measured distribution"); | |
1452 | legend->AddEntry(convolutedProj, "R #otimes unfolded distribution"); | |
1453 | legend->SetFillColor(0); | |
1454 | legend->Draw(); | |
1455 | ||
0b4bfd98 | 1456 | //TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded")); |
1457 | //diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1); | |
447c325d | 1458 | |
1459 | /*Float_t chi2 = 0; | |
1460 | Float_t chi = 0; | |
dd701109 | 1461 | fLastChi2MCLimit = 0; |
447c325d | 1462 | Int_t limit = 0; |
dd701109 | 1463 | for (Int_t i=2; i<=200; ++i) |
cfc19dd5 | 1464 | { |
dd701109 | 1465 | if (proj->GetBinContent(i) != 0) |
cfc19dd5 | 1466 | { |
dd701109 | 1467 | Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i); |
cfc19dd5 | 1468 | chi2 += value * value; |
447c325d | 1469 | chi += TMath::Abs(value); |
dd701109 | 1470 | |
447c325d | 1471 | //printf("%d %f\n", i, chi); |
cfc19dd5 | 1472 | |
447c325d | 1473 | if (chi2 < 0.2) |
1474 | fLastChi2MCLimit = i; | |
cfc19dd5 | 1475 | |
447c325d | 1476 | if (chi < 3) |
1477 | limit = i; | |
cfc19dd5 | 1478 | |
447c325d | 1479 | } |
1480 | }*/ | |
9ca6be09 | 1481 | |
0b4bfd98 | 1482 | /*chi2 = 0; |
447c325d | 1483 | Float_t chi = 0; |
1484 | Int_t limit = 0; | |
1485 | for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i) | |
cfc19dd5 | 1486 | { |
447c325d | 1487 | if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0) |
cfc19dd5 | 1488 | { |
447c325d | 1489 | Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i); |
1490 | if (value > 1e8) | |
1491 | value = 1e8; //prevent arithmetic exception | |
1492 | else if (value < -1e8) | |
1493 | value = -1e8; | |
1494 | if (i > 1) | |
1495 | { | |
1496 | chi2 += value * value; | |
1497 | chi += TMath::Abs(value); | |
1498 | } | |
1499 | diffMCUnfolded->SetBinContent(i, value); | |
cfc19dd5 | 1500 | } |
447c325d | 1501 | else |
1502 | { | |
1503 | //printf("diffMCUnfolded bin %d set to 0\n", i); | |
1504 | diffMCUnfolded->SetBinContent(i, 0); | |
1505 | } | |
1506 | if (chi2 < 1000) | |
1507 | fLastChi2MCLimit = i; | |
1508 | if (chi < 1000) | |
1509 | limit = i; | |
1510 | if (i == 150) | |
1511 | fLastChi2MC = chi2; | |
cfc19dd5 | 1512 | } |
0a173978 | 1513 | |
447c325d | 1514 | printf("limits %d %d\n", fLastChi2MCLimit, limit); |
1515 | fLastChi2MCLimit = limit; | |
1516 | ||
0b4bfd98 | 1517 | printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);*/ |
0a173978 | 1518 | |
447c325d | 1519 | if (!simple) |
1520 | { | |
0b4bfd98 | 1521 | /*canvas1->cd(3); |
447c325d | 1522 | |
0b4bfd98 | 1523 | diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(unfolded)"); |
447c325d | 1524 | //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20); |
1525 | diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200); | |
1526 | diffMCUnfolded->DrawCopy("HIST"); | |
0a173978 | 1527 | |
447c325d | 1528 | TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5); |
1529 | for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i) | |
1530 | fluctuation->Fill(diffMCUnfolded->GetBinContent(i)); | |
3602328d | 1531 | |
0b4bfd98 | 1532 | //new TCanvas; fluctuation->DrawCopy(); |
1533 | delete fluctuation;*/ | |
447c325d | 1534 | |
1535 | /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85); | |
1536 | legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected"); | |
1537 | legend->AddEntry(fMultiplicityMC, "MC"); | |
1538 | legend->AddEntry(fMultiplicityESD, "ESD"); | |
1539 | legend->Draw();*/ | |
1540 | ||
1541 | canvas1->cd(4); | |
039e265e | 1542 | residual->GetYaxis()->SetRangeUser(-5, 5); |
447c325d | 1543 | residual->GetXaxis()->SetRangeUser(0, 200); |
1544 | residual->DrawCopy(); | |
1545 | ||
1546 | canvas1->cd(5); | |
1547 | ||
1548 | TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio"); | |
1549 | ratio->Divide(mcHist); | |
1550 | ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC"); | |
1551 | ratio->GetYaxis()->SetRangeUser(0.5, 1.5); | |
1552 | ratio->GetXaxis()->SetRangeUser(0, 200); | |
1553 | ratio->SetStats(kFALSE); | |
1554 | ratio->Draw("HIST"); | |
1555 | ||
1556 | Double_t ratioChi2 = 0; | |
0b4bfd98 | 1557 | fRatioAverage = 0; |
447c325d | 1558 | fLastChi2MCLimit = 0; |
1559 | for (Int_t i=2; i<=150; ++i) | |
1560 | { | |
1561 | Float_t value = ratio->GetBinContent(i) - 1; | |
1562 | if (value > 1e8) | |
1563 | value = 1e8; //prevent arithmetic exception | |
1564 | else if (value < -1e8) | |
1565 | value = -1e8; | |
1566 | ||
1567 | ratioChi2 += value * value; | |
0b4bfd98 | 1568 | fRatioAverage += TMath::Abs(value); |
447c325d | 1569 | |
1570 | if (ratioChi2 < 0.1) | |
1571 | fLastChi2MCLimit = i; | |
1572 | } | |
0b4bfd98 | 1573 | fRatioAverage /= 149; |
447c325d | 1574 | |
0b4bfd98 | 1575 | printf("Sum over (ratio-1)^2 (2..150) is %f; average of |ratio-1| is %f\n", ratioChi2, fRatioAverage); |
447c325d | 1576 | // TODO FAKE |
1577 | fLastChi2MC = ratioChi2; | |
0b4bfd98 | 1578 | |
1579 | // FFT of ratio | |
1580 | canvas1->cd(6); | |
1581 | const Int_t kFFT = 128; | |
1582 | Double_t fftReal[kFFT]; | |
1583 | Double_t fftImag[kFFT]; | |
1584 | ||
1585 | for (Int_t i=0; i<kFFT; ++i) | |
1586 | { | |
1587 | fftReal[i] = ratio->GetBinContent(i+1+10); | |
1588 | // test: ;-) | |
1589 | //fftReal[i] = cos(TMath::Pi() * 5 * 2 * i / 128); | |
1590 | fftImag[i] = 0; | |
1591 | } | |
1592 | ||
04054a84 | 1593 | FFT(-1, TMath::Nint(TMath::Log(kFFT) / TMath::Log(2)), fftReal, fftImag); |
0b4bfd98 | 1594 | |
1595 | TH1* fftResult = (TH1*) ratio->Clone("fftResult"); | |
1596 | fftResult->SetTitle("FFT;true multiplicity;coeff. (10...137)"); | |
1597 | TH1* fftResult2 = (TH1*) ratio->Clone("fftResult2"); | |
1598 | TH1* fftResult3 = (TH1*) ratio->Clone("fftResult3"); | |
1599 | fftResult->Reset(); | |
1600 | fftResult2->Reset(); | |
1601 | fftResult3->Reset(); | |
1602 | fftResult->GetYaxis()->UnZoom(); | |
1603 | fftResult2->GetYaxis()->UnZoom(); | |
1604 | fftResult3->GetYaxis()->UnZoom(); | |
1605 | ||
0f67a57c | 1606 | //Printf("AFTER FFT"); |
0b4bfd98 | 1607 | for (Int_t i=0; i<kFFT; ++i) |
1608 | { | |
0f67a57c | 1609 | //Printf("%d: %f", i, fftReal[i]); |
0b4bfd98 | 1610 | fftResult->SetBinContent(i+1, fftReal[i]); |
1611 | /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5) | |
1612 | { | |
1613 | Printf("Nulled %d", i); | |
1614 | fftReal[i] = 0; | |
1615 | }*/ | |
1616 | } | |
1617 | ||
1618 | fftResult->SetLineColor(1); | |
1619 | fftResult->DrawCopy(); | |
1620 | ||
1621 | ||
0f67a57c | 1622 | //Printf("IMAG"); |
0b4bfd98 | 1623 | for (Int_t i=0; i<kFFT; ++i) |
1624 | { | |
0f67a57c | 1625 | //Printf("%d: %f", i, fftImag[i]); |
0b4bfd98 | 1626 | fftResult2->SetBinContent(i+1, fftImag[i]); |
1627 | ||
1628 | fftResult3->SetBinContent(i+1, TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i])); | |
1629 | } | |
1630 | ||
1631 | fftResult2->SetLineColor(2); | |
1632 | fftResult2->DrawCopy("SAME"); | |
1633 | ||
1634 | fftResult3->SetLineColor(4); | |
1635 | fftResult3->DrawCopy("SAME"); | |
1636 | ||
1637 | for (Int_t i=1; i<kFFT - 1; ++i) | |
1638 | { | |
1639 | if (TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]) > 3) | |
1640 | { | |
1641 | fftReal[i-1] = 0; | |
1642 | fftReal[i] = 0; | |
1643 | fftReal[i+1] = 0; | |
1644 | fftImag[i-1] = 0; | |
1645 | fftImag[i] = 0; | |
1646 | fftImag[i+1] = 0; | |
1647 | //fftReal[i] = (fftReal[i-1] + fftReal[i+1]) / 2; | |
1648 | //fftImag[i] = (fftImag[i-1] + fftImag[i+1]) / 2; | |
0f67a57c | 1649 | //Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]); |
0b4bfd98 | 1650 | } |
1651 | } | |
1652 | ||
04054a84 | 1653 | FFT(1, TMath::Nint(TMath::Log(kFFT) / TMath::Log(2)), fftReal, fftImag); |
0b4bfd98 | 1654 | |
1655 | TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4"); | |
1656 | fftResult4->Reset(); | |
1657 | ||
0f67a57c | 1658 | //Printf("AFTER BACK-TRAFO"); |
0b4bfd98 | 1659 | for (Int_t i=0; i<kFFT; ++i) |
1660 | { | |
0f67a57c | 1661 | //Printf("%d: %f", i, fftReal[i]); |
0b4bfd98 | 1662 | fftResult4->SetBinContent(i+1+10, fftReal[i]); |
1663 | } | |
1664 | ||
1665 | canvas1->cd(5); | |
1666 | fftResult4->SetLineColor(4); | |
1667 | fftResult4->DrawCopy("SAME"); | |
1668 | ||
1669 | // plot (MC - Unfolded) / error (MC) | |
1670 | canvas1->cd(3); | |
1671 | ||
1672 | TH1* diffMCUnfolded2 = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded2")); | |
1673 | diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1); | |
1674 | ||
1675 | Int_t ndfQual[kQualityRegions]; | |
1676 | for (Int_t region=0; region<kQualityRegions; ++region) | |
1677 | { | |
1678 | fQuality[region] = 0; | |
1679 | ndfQual[region] = 0; | |
1680 | } | |
1681 | ||
1682 | ||
1683 | Double_t newChi2 = 0; | |
6d81c2de | 1684 | Double_t newChi2Limit150 = 0; |
0b4bfd98 | 1685 | Int_t ndf = 0; |
6d81c2de | 1686 | for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgkMaxParams+1); ++i) |
0b4bfd98 | 1687 | { |
1688 | Double_t value = 0; | |
1689 | if (proj->GetBinError(i) > 0) | |
1690 | { | |
1691 | value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i); | |
1692 | newChi2 += value * value; | |
1693 | if (i > 1 && i <= mcBinLimit) | |
6d81c2de | 1694 | newChi2Limit150 += value * value; |
0b4bfd98 | 1695 | ++ndf; |
1696 | ||
1697 | for (Int_t region=0; region<kQualityRegions; ++region) | |
1698 | 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 | |
1699 | { | |
1700 | fQuality[region] += TMath::Abs(value); | |
1701 | ++ndfQual[region]; | |
1702 | } | |
1703 | } | |
1704 | ||
1705 | diffMCUnfolded2->SetBinContent(i, value); | |
1706 | } | |
1707 | ||
1708 | // normalize region to the number of entries | |
1709 | for (Int_t region=0; region<kQualityRegions; ++region) | |
1710 | { | |
1711 | if (ndfQual[region] > 0) | |
1712 | fQuality[region] /= ndfQual[region]; | |
1713 | Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]); | |
1714 | } | |
1715 | ||
1716 | if (mcBinLimit > 1) | |
1717 | { | |
1718 | // TODO another FAKE | |
6d81c2de | 1719 | fLastChi2MC = newChi2Limit150 / (mcBinLimit - 1); |
1720 | Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2Limit150, mcBinLimit - 1, fLastChi2MC); | |
0b4bfd98 | 1721 | } |
1722 | else | |
1723 | fLastChi2MC = -1; | |
1724 | ||
144ff489 | 1725 | 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 | 1726 | |
1727 | diffMCUnfolded2->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(MC)"); | |
039e265e | 1728 | diffMCUnfolded2->GetYaxis()->SetRangeUser(-5, 5); |
0b4bfd98 | 1729 | diffMCUnfolded2->GetXaxis()->SetRangeUser(0, 200); |
1730 | diffMCUnfolded2->DrawCopy("HIST"); | |
447c325d | 1731 | } |
3602328d | 1732 | |
b477f4f2 | 1733 | canvas1->SaveAs(Form("%s.gif", canvas1->GetName())); |
0a173978 | 1734 | } |
1735 | ||
1736 | //____________________________________________________________________ | |
0b4bfd98 | 1737 | void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y) |
1738 | { | |
1739 | /*------------------------------------------------------------------------- | |
1740 | This computes an in-place complex-to-complex FFT | |
1741 | x and y are the real and imaginary arrays of 2^m points. | |
1742 | dir = 1 gives forward transform | |
1743 | dir = -1 gives reverse transform | |
1744 | ||
1745 | Formula: forward | |
1746 | N-1 | |
1747 | --- | |
1748 | 1 \ - j k 2 pi n / N | |
1749 | X(n) = --- > x(k) e = forward transform | |
1750 | N / n=0..N-1 | |
1751 | --- | |
1752 | k=0 | |
1753 | ||
1754 | Formula: reverse | |
1755 | N-1 | |
1756 | --- | |
1757 | \ j k 2 pi n / N | |
1758 | X(n) = > x(k) e = forward transform | |
1759 | / n=0..N-1 | |
1760 | --- | |
1761 | k=0 | |
1762 | */ | |
1763 | ||
1764 | Long_t nn, i, i1, j, k, i2, l, l1, l2; | |
1765 | Double_t c1, c2, tx, ty, t1, t2, u1, u2, z; | |
1766 | ||
1767 | /* Calculate the number of points */ | |
1768 | nn = 1; | |
1769 | for (i = 0; i < m; i++) | |
1770 | nn *= 2; | |
1771 | ||
1772 | /* Do the bit reversal */ | |
1773 | i2 = nn >> 1; | |
1774 | j = 0; | |
1775 | for (i= 0; i < nn - 1; i++) { | |
1776 | if (i < j) { | |
1777 | tx = x[i]; | |
1778 | ty = y[i]; | |
1779 | x[i] = x[j]; | |
1780 | y[i] = y[j]; | |
1781 | x[j] = tx; | |
1782 | y[j] = ty; | |
1783 | } | |
1784 | k = i2; | |
1785 | while (k <= j) { | |
1786 | j -= k; | |
1787 | k >>= 1; | |
1788 | } | |
1789 | j += k; | |
1790 | } | |
1791 | ||
1792 | /* Compute the FFT */ | |
1793 | c1 = -1.0; | |
1794 | c2 = 0.0; | |
1795 | l2 = 1; | |
1796 | for (l = 0; l < m; l++) { | |
1797 | l1 = l2; | |
1798 | l2 <<= 1; | |
1799 | u1 = 1.0; | |
1800 | u2 = 0.0; | |
1801 | for (j = 0;j < l1; j++) { | |
1802 | for (i = j; i < nn; i += l2) { | |
1803 | i1 = i + l1; | |
1804 | t1 = u1 * x[i1] - u2 * y[i1]; | |
1805 | t2 = u1 * y[i1] + u2 * x[i1]; | |
1806 | x[i1] = x[i] - t1; | |
1807 | y[i1] = y[i] - t2; | |
1808 | x[i] += t1; | |
1809 | y[i] += t2; | |
1810 | } | |
1811 | z = u1 * c1 - u2 * c2; | |
1812 | u2 = u1 * c2 + u2 * c1; | |
1813 | u1 = z; | |
1814 | } | |
1815 | c2 = TMath::Sqrt((1.0 - c1) / 2.0); | |
1816 | if (dir == 1) | |
1817 | c2 = -c2; | |
1818 | c1 = TMath::Sqrt((1.0 + c1) / 2.0); | |
1819 | } | |
1820 | ||
1821 | /* Scaling for forward transform */ | |
1822 | if (dir == 1) { | |
1823 | for (i=0;i<nn;i++) { | |
1824 | x[i] /= (Double_t)nn; | |
1825 | y[i] /= (Double_t)nn; | |
1826 | } | |
1827 | } | |
1828 | } | |
1829 | ||
1830 | //____________________________________________________________________ | |
6d81c2de | 1831 | void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals, Float_t* ratioAverage) const |
cfc19dd5 | 1832 | { |
1833 | // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD | |
1834 | // These values are computed during DrawComparison, thus this function picks up the | |
1835 | // last calculation | |
1836 | ||
dd701109 | 1837 | if (mc) |
1838 | *mc = fLastChi2MC; | |
1839 | if (mcLimit) | |
1840 | *mcLimit = fLastChi2MCLimit; | |
1841 | if (residuals) | |
1842 | *residuals = fLastChi2Residuals; | |
0b4bfd98 | 1843 | if (ratioAverage) |
1844 | *ratioAverage = fRatioAverage; | |
cfc19dd5 | 1845 | } |
1846 | ||
cfc19dd5 | 1847 | //____________________________________________________________________ |
1848 | TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType) | |
1849 | { | |
1850 | // | |
1851 | // returns the corresponding MC spectrum | |
1852 | // | |
1853 | ||
1854 | switch (eventType) | |
1855 | { | |
1856 | case kTrVtx : return fMultiplicityVtx[i]; break; | |
1857 | case kMB : return fMultiplicityMB[i]; break; | |
1858 | case kINEL : return fMultiplicityINEL[i]; break; | |
1859 | } | |
1860 | ||
1861 | return 0; | |
1862 | } | |
1863 | ||
0f67a57c | 1864 | //____________________________________________________________________ |
1865 | TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max) | |
1866 | { | |
1867 | // calculate standard deviation of (results[0] - results[k]) k=1...max-1 | |
1868 | // per bin one gets: sigma(r[0] - r[n]) / r[0] | |
1869 | ||
1870 | TH1* standardDeviation = (TH1*) results[0]->Clone("standardDeviation"); | |
1871 | standardDeviation->Reset(); | |
1872 | ||
1873 | for (Int_t x=1; x<=results[0]->GetNbinsX(); x++) | |
1874 | { | |
1875 | if (results[0]->GetBinContent(x) > 0) | |
1876 | { | |
1877 | Double_t average = 0; | |
1878 | for (Int_t n=1; n<max; ++n) | |
1879 | average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x); | |
1880 | average /= max-1; | |
1881 | ||
1882 | Double_t variance = 0; | |
1883 | for (Int_t n=1; n<max; ++n) | |
1884 | { | |
1885 | Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average; | |
1886 | variance += value * value; | |
1887 | } | |
1888 | variance /= max-1; | |
1889 | ||
1890 | Double_t standardDev = TMath::Sqrt(variance); | |
1891 | standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x)); | |
1892 | //Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x)); | |
1893 | } | |
1894 | } | |
1895 | ||
1896 | return standardDeviation; | |
1897 | } | |
1898 | ||
cfc19dd5 | 1899 | //____________________________________________________________________ |
6d81c2de | 1900 | TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo) |
0a173978 | 1901 | { |
1902 | // | |
0b4bfd98 | 1903 | // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix |
1904 | // the function unfolds the spectrum using the default response matrix and several modified ones | |
1905 | // the modified ones are created by randomizing each cell using poisson statistics with the mean = bin value | |
1906 | // these unfolded results are compared to the first result gained with the default response OR to the histogram given | |
1907 | // in <compareTo> (optional) | |
0a173978 | 1908 | // |
0b4bfd98 | 1909 | // returns the error assigned to the measurement |
1910 | // | |
1911 | ||
6d81c2de | 1912 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); |
1913 | ||
0b4bfd98 | 1914 | // initialize seed with current time |
1915 | gRandom->SetSeed(0); | |
1916 | ||
1917 | const Int_t kErrorIterations = 150; | |
1918 | ||
1919 | TH1* maxError = 0; | |
1920 | TH1* firstResult = 0; | |
1921 | ||
1922 | TH1** results = new TH1*[kErrorIterations]; | |
1923 | ||
1924 | for (Int_t n=0; n<kErrorIterations; ++n) | |
1925 | { | |
1926 | Printf("Iteration %d of %d...", n, kErrorIterations); | |
1927 | ||
6d81c2de | 1928 | // only done for chi2 minimization |
1929 | Bool_t createBigBin = kFALSE; | |
1930 | if (methodType == kChi2Minimization) | |
1931 | createBigBin = kTRUE; | |
1932 | ||
1933 | SetupCurrentHists(inputRange, fullPhaseSpace, eventType, createBigBin); | |
0b4bfd98 | 1934 | |
1935 | TH1* measured = (TH1*) fCurrentESD->Clone("measured"); | |
1936 | ||
1937 | if (n > 0) | |
1938 | { | |
1939 | if (randomizeResponse) | |
1940 | { | |
1941 | // randomize response matrix | |
1942 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
1943 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) | |
1944 | fCurrentCorrelation->SetBinContent(i, j, gRandom->Poisson(fCurrentCorrelation->GetBinContent(i, j))); | |
1945 | } | |
1946 | ||
1947 | if (randomizeMeasured) | |
1948 | { | |
1949 | // randomize measured spectrum | |
1950 | for (Int_t x=1; x<=measured->GetNbinsX(); x++) // mult. axis | |
1951 | { | |
1952 | Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x)); | |
1953 | measured->SetBinContent(x, randomValue); | |
1954 | measured->SetBinError(x, TMath::Sqrt(randomValue)); | |
1955 | } | |
1956 | } | |
1957 | } | |
1958 | ||
6d81c2de | 1959 | // only for bayesian method we have to do it before the call to Unfold... |
1960 | if (methodType == kBayesian) | |
0b4bfd98 | 1961 | { |
6d81c2de | 1962 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) |
0b4bfd98 | 1963 | { |
6d81c2de | 1964 | // with this it is normalized to 1 |
1965 | Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); | |
1966 | ||
1967 | // with this normalized to the given efficiency | |
1968 | if (fCurrentEfficiency->GetBinContent(i) > 0) | |
1969 | sum /= fCurrentEfficiency->GetBinContent(i); | |
0b4bfd98 | 1970 | else |
6d81c2de | 1971 | sum = 0; |
1972 | ||
1973 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) | |
0b4bfd98 | 1974 | { |
6d81c2de | 1975 | if (sum > 0) |
1976 | { | |
1977 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); | |
1978 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
1979 | } | |
1980 | else | |
1981 | { | |
1982 | fCurrentCorrelation->SetBinContent(i, j, 0); | |
1983 | fCurrentCorrelation->SetBinError(i, j, 0); | |
1984 | } | |
0b4bfd98 | 1985 | } |
1986 | } | |
1987 | } | |
1988 | ||
1989 | TH1* result = 0; | |
1990 | if (n == 0 && compareTo) | |
1991 | { | |
1992 | // in this case we just store the histogram we want to compare to | |
1993 | result = (TH1*) compareTo->Clone("compareTo"); | |
1994 | result->Sumw2(); | |
1995 | } | |
1996 | else | |
6d81c2de | 1997 | { |
1998 | result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n)); | |
0b4bfd98 | 1999 | |
6d81c2de | 2000 | Int_t returnCode = -1; |
2001 | if (methodType == kChi2Minimization) | |
2002 | { | |
2003 | returnCode = UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, kFALSE); | |
2004 | } | |
2005 | else if (methodType == kBayesian) | |
2006 | returnCode = UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, fgBayesianSmoothing, fgBayesianIterations); | |
2007 | ||
2008 | if (returnCode != 0) | |
2009 | return 0; | |
2010 | } | |
0b4bfd98 | 2011 | |
2012 | // normalize | |
2013 | result->Scale(1.0 / result->Integral()); | |
2014 | ||
2015 | if (n == 0) | |
2016 | { | |
2017 | firstResult = (TH1*) result->Clone("firstResult"); | |
2018 | ||
2019 | maxError = (TH1*) result->Clone("maxError"); | |
2020 | maxError->Reset(); | |
2021 | } | |
2022 | else | |
2023 | { | |
2024 | // calculate ratio | |
2025 | TH1* ratio = (TH1*) firstResult->Clone("ratio"); | |
2026 | ratio->Divide(result); | |
2027 | ||
2028 | // find max. deviation | |
2029 | for (Int_t x=1; x<=ratio->GetNbinsX(); x++) | |
2030 | maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x)))); | |
2031 | ||
2032 | delete ratio; | |
2033 | } | |
2034 | ||
2035 | results[n] = result; | |
2036 | } | |
2037 | ||
2038 | // find covariance matrix | |
2039 | // results[n] is X_x | |
2040 | // cov. matrix is M_xy = E ( (X_x - E(X_x)) * (X_y - E(X_y))), with E() = expectation value | |
2041 | ||
2042 | Int_t nBins = results[0]->GetNbinsX(); | |
2043 | Float_t lowEdge = results[0]->GetXaxis()->GetBinLowEdge(1); | |
2044 | Float_t upEdge = results[0]->GetXaxis()->GetBinUpEdge(nBins); | |
2045 | ||
2046 | // find average, E(X) | |
2047 | TProfile* average = new TProfile("average", "average", nBins, lowEdge, upEdge); | |
2048 | for (Int_t n=1; n<kErrorIterations; ++n) | |
2049 | for (Int_t x=1; x<=results[n]->GetNbinsX(); x++) | |
2050 | average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x)); | |
2051 | //new TCanvas; average->DrawClone(); | |
2052 | ||
2053 | // find cov. matrix | |
2054 | TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge); | |
2055 | ||
2056 | for (Int_t n=1; n<kErrorIterations; ++n) | |
2057 | for (Int_t x=1; x<=results[n]->GetNbinsX(); x++) | |
2058 | for (Int_t y=1; y<=results[n]->GetNbinsX(); y++) | |
2059 | { | |
2060 | // (X_x - E(X_x)) * (X_y - E(X_y) | |
2061 | Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y)); | |
2062 | covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov); | |
2063 | } | |
2064 | new TCanvas; covMatrix->DrawClone("COLZ"); | |
2065 | ||
6d81c2de | 2066 | // // fill 2D histogram that contains deviation from first |
2067 | // TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01); | |
2068 | // for (Int_t n=1; n<kErrorIterations; ++n) | |
2069 | // for (Int_t x=1; x<=results[n]->GetNbinsX(); x++) | |
2070 | // deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x)); | |
2071 | // //new TCanvas; deviations->DrawCopy("COLZ"); | |
2072 | // | |
2073 | // // get standard deviation "by hand" | |
2074 | // for (Int_t x=1; x<=nBins; x++) | |
2075 | // { | |
2076 | // if (results[0]->GetBinContent(x) > 0) | |
2077 | // { | |
2078 | // TH1* proj = deviations->ProjectionY("projRMS", x, x, "e"); | |
2079 | // Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact | |
2080 | // //standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x)); | |
2081 | // Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x)); | |
2082 | // } | |
2083 | // } | |
2084 | ||
0f67a57c | 2085 | TH1* standardDeviation = CalculateStdDev(results, kErrorIterations); |
0b4bfd98 | 2086 | |
2087 | // compare maxError to RMS of profile (variable name: average) | |
2088 | // first: calculate rms in percent of value | |
2089 | TH1* rmsError = (TH1*) maxError->Clone("rmsError"); | |
2090 | rmsError->Reset(); | |
2091 | ||
2092 | // enable error to be standard deviation (see http://root.cern.ch/root/html/TProfile.html#TProfile:SetErrorOption) | |
2093 | average->SetErrorOption("s"); | |
2094 | for (Int_t x=1; x<=rmsError->GetNbinsX(); x++) | |
2095 | if (average->GetBinContent(x) > 0) | |
2096 | rmsError->SetBinContent(x, average->GetBinError(x) / average->GetBinContent(x)); | |
2097 | ||
2098 | // find maxError deviation from average (not from "first result"/mc as above) | |
2099 | TH1* maxError2 = (TH1*) maxError->Clone("maxError2"); | |
2100 | maxError2->Reset(); | |
2101 | for (Int_t n=1; n<kErrorIterations; ++n) | |
2102 | for (Int_t x=1; x<=results[n]->GetNbinsX(); x++) | |
2103 | if (average->GetBinContent(x) > 0) | |
2104 | maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x)))); | |
2105 | ||
6d81c2de | 2106 | //new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME"); |
0b4bfd98 | 2107 | |
2108 | // plot difference between average and MC/first result | |
2109 | TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio"); | |
2110 | averageFirstRatio->Reset(); | |
2111 | averageFirstRatio->Divide(results[0], average); | |
2112 | ||
6d81c2de | 2113 | //new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME"); |
2114 | //new TCanvas; averageFirstRatio->DrawCopy(); | |
0b4bfd98 | 2115 | |
2116 | // clean up | |
2117 | for (Int_t n=0; n<kErrorIterations; ++n) | |
2118 | delete results[n]; | |
2119 | delete[] results; | |
2120 | ||
2121 | // fill into result histogram | |
0b4bfd98 | 2122 | for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i) |
2123 | fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i)); | |
cfc19dd5 | 2124 | |
0b4bfd98 | 2125 | for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i) |
2126 | fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i)); | |
447c325d | 2127 | |
0b4bfd98 | 2128 | return standardDeviation; |
2129 | } | |
2130 | ||
2131 | //____________________________________________________________________ | |
6d81c2de | 2132 | void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Bool_t determineError) |
0b4bfd98 | 2133 | { |
2134 | // | |
2135 | // correct spectrum using bayesian method | |
2136 | // | |
2137 | ||
2138 | // initialize seed with current time | |
2139 | gRandom->SetSeed(0); | |
2140 | ||
2141 | SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE); | |
0a173978 | 2142 | |
2143 | // normalize correction for given nPart | |
9ca6be09 | 2144 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) |
0a173978 | 2145 | { |
dd701109 | 2146 | // with this it is normalized to 1 |
2147 | Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); | |
2148 | ||
2149 | // with this normalized to the given efficiency | |
2150 | if (fCurrentEfficiency->GetBinContent(i) > 0) | |
2151 | sum /= fCurrentEfficiency->GetBinContent(i); | |
2152 | else | |
2153 | sum = 0; | |
2154 | ||
9ca6be09 | 2155 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) |
0a173978 | 2156 | { |
dd701109 | 2157 | if (sum > 0) |
2158 | { | |
2159 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); | |
2160 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
2161 | } | |
2162 | else | |
2163 | { | |
2164 | fCurrentCorrelation->SetBinContent(i, j, 0); | |
2165 | fCurrentCorrelation->SetBinError(i, j, 0); | |
2166 | } | |
0a173978 | 2167 | } |
2168 | } | |
2169 | ||
0b4bfd98 | 2170 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); |
2171 | ||
6d81c2de | 2172 | if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], regPar, nIterations) != 0) |
0b4bfd98 | 2173 | return; |
2174 | ||
0b4bfd98 | 2175 | if (!determineError) |
447c325d | 2176 | { |
0b4bfd98 | 2177 | Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated."); |
2178 | return; | |
2179 | } | |
447c325d | 2180 | |
0b4bfd98 | 2181 | // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics |
2182 | // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured | |
0f67a57c | 2183 | // spectrum calculated. This is performed N times and the sigma is taken as the statistical |
0b4bfd98 | 2184 | // error of the unfolding method itself. |
2185 | ||
0b4bfd98 | 2186 | const Int_t kErrorIterations = 20; |
2187 | ||
2188 | printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations); | |
2189 | ||
2190 | TH1* randomized = (TH1*) fCurrentESD->Clone("randomized"); | |
0f67a57c | 2191 | TH1* resultArray[kErrorIterations+1]; |
0b4bfd98 | 2192 | for (Int_t n=0; n<kErrorIterations; ++n) |
2193 | { | |
2194 | // randomize the content of clone following a poisson with the mean = the value of that bin | |
2195 | for (Int_t x=1; x<=randomized->GetNbinsX(); x++) // mult. axis | |
447c325d | 2196 | { |
0b4bfd98 | 2197 | Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x)); |
2198 | //printf("%d --> %d\n", fCurrentESD->GetBinContent(x), randomValue); | |
2199 | randomized->SetBinContent(x, randomValue); | |
2200 | randomized->SetBinError(x, TMath::Sqrt(randomValue)); | |
447c325d | 2201 | } |
447c325d | 2202 | |
0f67a57c | 2203 | TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2"); |
6d81c2de | 2204 | result2->Reset(); |
2205 | if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2, regPar, nIterations) != 0) | |
0b4bfd98 | 2206 | return; |
cfc19dd5 | 2207 | |
0b4bfd98 | 2208 | result2->Scale(1.0 / result2->Integral()); |
cfc19dd5 | 2209 | |
0f67a57c | 2210 | resultArray[n+1] = result2; |
0b4bfd98 | 2211 | } |
6d81c2de | 2212 | delete randomized; |
0f67a57c | 2213 | |
2214 | resultArray[0] = fMultiplicityESDCorrected[correlationID]; | |
2215 | TH1* error = CalculateStdDev(resultArray, kErrorIterations+1); | |
2216 | ||
2217 | for (Int_t n=0; n<kErrorIterations; ++n) | |
2218 | delete resultArray[n+1]; | |
447c325d | 2219 | |
0b4bfd98 | 2220 | for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i) |
0f67a57c | 2221 | fMultiplicityESDCorrected[correlationID]->SetBinError(i, error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i)); |
447c325d | 2222 | |
0f67a57c | 2223 | delete error; |
0b4bfd98 | 2224 | } |
447c325d | 2225 | |
0b4bfd98 | 2226 | //____________________________________________________________________ |
6d81c2de | 2227 | Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations) |
0b4bfd98 | 2228 | { |
2229 | // | |
2230 | // unfolds a spectrum | |
2231 | // | |
2232 | ||
2233 | if (measured->Integral() <= 0) | |
dd701109 | 2234 | { |
0b4bfd98 | 2235 | Printf("AliMultiplicityCorrection::UnfoldWithBayesian: ERROR: The measured spectrum is empty"); |
6d81c2de | 2236 | return -1; |
dd701109 | 2237 | } |
6127aca6 | 2238 | |
0b4bfd98 | 2239 | const Int_t kStartBin = 0; |
6127aca6 | 2240 | |
0f67a57c | 2241 | const Int_t kMaxM = fgkMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis |
6d81c2de | 2242 | const Int_t kMaxT = fgkMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis |
6127aca6 | 2243 | |
0f67a57c | 2244 | // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4) |
2245 | const Double_t kConvergenceLimit = kMaxT * 1e-6; | |
2246 | ||
0b4bfd98 | 2247 | // store information in arrays, to increase processing speed (~ factor 5) |
2248 | Double_t measuredCopy[kMaxM]; | |
0f67a57c | 2249 | Double_t measuredError[kMaxM]; |
0b4bfd98 | 2250 | Double_t prior[kMaxT]; |
2251 | Double_t result[kMaxT]; | |
2252 | Double_t efficiency[kMaxT]; | |
b477f4f2 | 2253 | |
0b4bfd98 | 2254 | Double_t response[kMaxT][kMaxM]; |
2255 | Double_t inverseResponse[kMaxT][kMaxM]; | |
dd701109 | 2256 | |
0b4bfd98 | 2257 | // for normalization |
2258 | Float_t measuredIntegral = measured->Integral(); | |
2259 | for (Int_t m=0; m<kMaxM; m++) | |
2260 | { | |
2261 | measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral; | |
0f67a57c | 2262 | measuredError[m] = measured->GetBinError(m+1) / measuredIntegral; |
0b4bfd98 | 2263 | |
2264 | for (Int_t t=0; t<kMaxT; t++) | |
2265 | { | |
6d81c2de | 2266 | response[t][m] = correlation->GetBinContent(t+1, m+1); |
0b4bfd98 | 2267 | inverseResponse[t][m] = 0; |
2268 | } | |
2269 | } | |
2270 | ||
2271 | for (Int_t t=0; t<kMaxT; t++) | |
2272 | { | |
6d81c2de | 2273 | efficiency[t] = aEfficiency->GetBinContent(t+1); |
0b4bfd98 | 2274 | prior[t] = measuredCopy[t]; |
2275 | result[t] = 0; | |
2276 | } | |
2277 | ||
2278 | // pick prior distribution | |
6d81c2de | 2279 | if (initialConditions) |
0b4bfd98 | 2280 | { |
2281 | printf("Using different starting conditions...\n"); | |
2282 | // for normalization | |
6d81c2de | 2283 | Float_t inputDistIntegral = initialConditions->Integral(); |
0b4bfd98 | 2284 | for (Int_t i=0; i<kMaxT; i++) |
6d81c2de | 2285 | prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral; |
0b4bfd98 | 2286 | } |
2287 | ||
0f67a57c | 2288 | //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5); |
447c325d | 2289 | |
6127aca6 | 2290 | // unfold... |
0f67a57c | 2291 | for (Int_t i=0; i<nIterations || nIterations < 0; i++) |
cfc19dd5 | 2292 | { |
2293 | //printf(" iteration %i \n", i); | |
2294 | ||
6127aca6 | 2295 | // calculate IR from Beyes theorem |
2296 | // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k) | |
9ca6be09 | 2297 | |
0f67a57c | 2298 | Double_t chi2Measured = 0; |
0b4bfd98 | 2299 | for (Int_t m=0; m<kMaxM; m++) |
cfc19dd5 | 2300 | { |
2301 | Float_t norm = 0; | |
0b4bfd98 | 2302 | for (Int_t t = kStartBin; t<kMaxT; t++) |
2303 | norm += response[t][m] * prior[t]; | |
2304 | ||
0f67a57c | 2305 | // calc. chi2: (measured - response * prior) / error |
2306 | if (measuredError[m] > 0) | |
2307 | { | |
2308 | Double_t value = (measuredCopy[m] - norm) / measuredError[m]; | |
2309 | chi2Measured += value * value; | |
2310 | } | |
2311 | ||
0b4bfd98 | 2312 | if (norm > 0) |
cfc19dd5 | 2313 | { |
0b4bfd98 | 2314 | for (Int_t t = kStartBin; t<kMaxT; t++) |
2315 | inverseResponse[t][m] = response[t][m] * prior[t] / norm; | |
2316 | } | |
2317 | else | |
2318 | { | |
2319 | for (Int_t t = kStartBin; t<kMaxT; t++) | |
2320 | inverseResponse[t][m] = 0; | |
6127aca6 | 2321 | } |
2322 | } | |
0f67a57c | 2323 | //Printf("chi2Measured of the last prior is %e", chi2Measured); |
cfc19dd5 | 2324 | |
0b4bfd98 | 2325 | for (Int_t t = kStartBin; t<kMaxT; t++) |
cfc19dd5 | 2326 | { |
2327 | Float_t value = 0; | |
0b4bfd98 | 2328 | for (Int_t m=0; m<kMaxM; m++) |
2329 | value += inverseResponse[t][m] * measuredCopy[m]; | |
cfc19dd5 | 2330 | |
0b4bfd98 | 2331 | if (efficiency[t] > 0) |
2332 | result[t] = value / efficiency[t]; | |
dd701109 | 2333 | else |
0b4bfd98 | 2334 | result[t] = 0; |
cfc19dd5 | 2335 | } |
2336 | ||
0f67a57c | 2337 | Double_t chi2LastIter = 0; |
b477f4f2 | 2338 | // regularization (simple smoothing) |
0b4bfd98 | 2339 | for (Int_t t=kStartBin; t<kMaxT; t++) |
cfc19dd5 | 2340 | { |
0b4bfd98 | 2341 | Float_t newValue = 0; |
2342 | // 0 bin excluded from smoothing | |
2343 | if (t > kStartBin+1 && t<kMaxT-1) | |
2344 | { | |
2345 | Float_t average = (result[t-1] + result[t] + result[t+1]) / 3; | |
9ca6be09 | 2346 | |
0b4bfd98 | 2347 | // weight the average with the regularization parameter |
2348 | newValue = (1 - regPar) * result[t] + regPar * average; | |
2349 | } | |
2350 | else | |
2351 | newValue = result[t]; | |
2352 | ||
0f67a57c | 2353 | // calculate chi2 (change from last iteration) |
2354 | if (prior[t] > 1e-5) | |
2355 | { | |
2356 | Double_t diff = (prior[t] - newValue) / prior[t]; | |
2357 | chi2LastIter += diff * diff; | |
2358 | } | |
2359 | ||
0b4bfd98 | 2360 | prior[t] = newValue; |
6127aca6 | 2361 | } |
0f67a57c | 2362 | //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter); |
2363 | //convergence->Fill(i+1, chi2LastIter); | |
9ca6be09 | 2364 | |
0f67a57c | 2365 | if (nIterations < 0 && chi2LastIter < kConvergenceLimit) |
cfc19dd5 | 2366 | { |
0f67a57c | 2367 | Printf("Stopped Bayesian unfolding after %d iterations at chi2(change since last iteration) of %e; chi2Measured of the last prior is %e", i, chi2LastIter, chi2Measured); |
2368 | break; | |
6127aca6 | 2369 | } |
dd701109 | 2370 | |
cfc19dd5 | 2371 | //if (i % 5 == 0) |
2372 | // DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType); | |
6127aca6 | 2373 | } // end of iterations |
2374 | ||
0f67a57c | 2375 | //new TCanvas; convergence->DrawCopy(); gPad->SetLogy(); |
0b4bfd98 | 2376 | //delete convergence; |
2377 | ||
0b4bfd98 | 2378 | for (Int_t t=0; t<kMaxT; t++) |
6d81c2de | 2379 | aResult->SetBinContent(t+1, result[t]); |
6127aca6 | 2380 | |
6d81c2de | 2381 | return 0; |
cfc19dd5 | 2382 | |
2383 | // ******** | |
2384 | // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995 | |
2385 | ||
dd701109 | 2386 | /*printf("Calculating covariance matrix. This may take some time...\n"); |
2387 | ||
2388 | // TODO check if this is the right one... | |
2389 | TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX()); | |
cfc19dd5 | 2390 | |
2391 | Int_t xBins = hInverseResponseBayes->GetNbinsX(); | |
2392 | Int_t yBins = hInverseResponseBayes->GetNbinsY(); | |
2393 | ||
2394 | // calculate "unfolding matrix" Mij | |
2395 | Float_t matrixM[251][251]; | |
2396 | for (Int_t i=1; i<=xBins; i++) | |
2397 | { | |
2398 | for (Int_t j=1; j<=yBins; j++) | |
2399 | { | |
2400 | if (fCurrentEfficiency->GetBinContent(i) > 0) | |
2401 | matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i); | |
2402 | else | |
2403 | matrixM[i-1][j-1] = 0; | |
2404 | } | |
2405 | } | |
2406 | ||
2407 | Float_t* vectorn = new Float_t[yBins]; | |
2408 | for (Int_t j=1; j<=yBins; j++) | |
2409 | vectorn[j-1] = fCurrentESD->GetBinContent(j); | |
2410 | ||
2411 | // first part of covariance matrix, depends on input distribution n(E) | |
2412 | Float_t cov1[251][251]; | |
2413 | ||
2414 | Float_t nEvents = fCurrentESD->Integral(); // N | |
2415 | ||
2416 | xBins = 20; | |
2417 | yBins = 20; | |
2418 | ||
2419 | for (Int_t k=0; k<xBins; k++) | |
2420 | { | |
2421 | printf("In Cov1: %d\n", k); | |
2422 | for (Int_t l=0; l<yBins; l++) | |
2423 | { | |
2424 | cov1[k][l] = 0; | |
2425 | ||
2426 | // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N) | |
2427 | for (Int_t j=0; j<yBins; j++) | |
2428 | cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j] | |
2429 | * (1.0 - vectorn[j] / nEvents); | |
2430 | ||
2431 | // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N | |
2432 | for (Int_t i=0; i<yBins; i++) | |
2433 | for (Int_t j=0; j<yBins; j++) | |
2434 | { | |
2435 | if (i == j) | |
2436 | continue; | |
2437 | cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i] | |
2438 | * vectorn[j] / nEvents; | |
2439 | } | |
2440 | } | |
2441 | } | |
2442 | ||
2443 | printf("Cov1 finished\n"); | |
2444 | ||
2445 | TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov"); | |
2446 | cov->Reset(); | |
2447 | ||
2448 | for (Int_t i=1; i<=xBins; i++) | |
2449 | for (Int_t j=1; j<=yBins; j++) | |
2450 | cov->SetBinContent(i, j, cov1[i-1][j-1]); | |
2451 | ||
2452 | new TCanvas; | |
2453 | cov->Draw("COLZ"); | |
2454 | ||
2455 | // second part of covariance matrix, depends on response matrix | |
2456 | Float_t cov2[251][251]; | |
2457 | ||
2458 | // Cov[P(Er|Cu), P(Es|Cu)] term | |
2459 | Float_t covTerm[100][100][100]; | |
2460 | for (Int_t r=0; r<yBins; r++) | |
2461 | for (Int_t u=0; u<xBins; u++) | |
2462 | for (Int_t s=0; s<yBins; s++) | |
2463 | { | |
2464 | if (r == s) | |
2465 | covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) | |
2466 | * (1.0 - hResponse->GetBinContent(u+1, r+1)); | |
2467 | else | |
2468 | covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) | |
2469 | * hResponse->GetBinContent(u+1, s+1); | |
2470 | } | |
2471 | ||
2472 | for (Int_t k=0; k<xBins; k++) | |
2473 | for (Int_t l=0; l<yBins; l++) | |
2474 | { | |
2475 | cov2[k][l] = 0; | |
2476 | printf("In Cov2: %d %d\n", k, l); | |
2477 | for (Int_t i=0; i<yBins; i++) | |
2478 | for (Int_t j=0; j<yBins; j++) | |
2479 | { | |
2480 | //printf("In Cov2: %d %d %d %d\n", k, l, i, j); | |
2481 | // calculate Cov(Mki, Mlj) = sum{ru},{su} ... | |
2482 | Float_t tmpCov = 0; | |
2483 | for (Int_t r=0; r<yBins; r++) | |
2484 | for (Int_t u=0; u<xBins; u++) | |
2485 | for (Int_t s=0; s<yBins; s++) | |
2486 | { | |
2487 | if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0 | |
2488 | || hResponse->GetBinContent(u+1, i+1) == 0) | |
2489 | continue; | |
2490 | ||
2491 | tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u) | |
2492 | * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u) | |
2493 | * covTerm[r][u][s]; | |
2494 | } | |
2495 | ||
2496 | cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov; | |
2497 | } | |
2498 | } | |
2499 | ||
2500 | printf("Cov2 finished\n"); | |
2501 | ||
2502 | for (Int_t i=1; i<=xBins; i++) | |
2503 | for (Int_t j=1; j<=yBins; j++) | |
2504 | cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]); | |
2505 | ||
2506 | new TCanvas; | |
dd701109 | 2507 | cov->Draw("COLZ");*/ |
2508 | } | |
2509 | ||
2510 | //____________________________________________________________________ | |
0b4bfd98 | 2511 | Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u) |
dd701109 | 2512 | { |
2513 | // | |
2514 | // helper function for the covariance matrix of the bayesian method | |
2515 | // | |
2516 | ||
2517 | Float_t result = 0; | |
2518 | ||
2519 | if (k == u && r == i) | |
2520 | result += 1.0 / hResponse->GetBinContent(u+1, r+1); | |
2521 | ||
2522 | if (k == u) | |
2523 | result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1); | |
2524 | ||
2525 | if (r == i) | |
2526 | result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1); | |
2527 | ||
2528 | result *= matrixM[k][i]; | |
2529 | ||
2530 | return result; | |
cfc19dd5 | 2531 | } |
2532 | ||
2533 | //____________________________________________________________________ | |
2534 | void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType) | |
2535 | { | |
2536 | // | |
2537 | // correct spectrum using bayesian method | |
2538 | // | |
2539 | ||
dd701109 | 2540 | Float_t regPar = 0; |
cfc19dd5 | 2541 | |
2542 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
2543 | Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); | |
2544 | ||
447c325d | 2545 | SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE); |
0b4bfd98 | 2546 | //normalize ESD |
2547 | fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); | |
cfc19dd5 | 2548 | |
2549 | // TODO should be taken from correlation map | |
2550 | //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX()); | |
2551 | ||
2552 | // normalize correction for given nPart | |
2553 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
2554 | { | |
2555 | Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); | |
2556 | //Double_t sum = sumHist->GetBinContent(i); | |
2557 | if (sum <= 0) | |
2558 | continue; | |
2559 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) | |
2560 | { | |
2561 | // npart sum to 1 | |
2562 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); | |
2563 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
2564 | } | |
6127aca6 | 2565 | } |
cfc19dd5 | 2566 | |
2567 | new TCanvas; | |
2568 | fCurrentCorrelation->Draw("COLZ"); | |
2569 | ||
2570 | // FAKE | |
2571 | fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff"); | |
2572 | ||
2573 | // pick prior distribution | |
2574 | TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior"); | |
2575 | Float_t norm = 1; //hPrior->Integral(); | |
2576 | for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) | |
2577 | hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm); | |
2578 | ||
2579 | // zero distribution | |
2580 | TH1F* zero = (TH1F*)hPrior->Clone("zero"); | |
2581 | ||
2582 | // define temp hist | |
2583 | TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp"); | |
2584 | hTemp->Reset(); | |
2585 | ||
2586 | // just a shortcut | |
2587 | TH2F* hResponse = (TH2F*) fCurrentCorrelation; | |
2588 | ||
2589 | // unfold... | |
dd701109 | 2590 | Int_t iterations = 25; |
cfc19dd5 | 2591 | for (Int_t i=0; i<iterations; i++) |
2592 | { | |
2593 | //printf(" iteration %i \n", i); | |
2594 | ||
2595 | for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) | |
2596 | { | |
2597 | Float_t value = 0; | |
2598 | for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++) | |
2599 | value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t); | |
2600 | hTemp->SetBinContent(m, value); | |
2601 | //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value); | |
2602 | } | |
2603 | ||
2604 | // regularization (simple smoothing) | |
2605 | TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed"); | |
2606 | ||
2607 | for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++) | |
2608 | { | |
2609 | Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1) | |
2610 | + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t) | |
2611 | + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.; | |
2612 | average *= hTrueSmoothed->GetBinWidth(t); | |
2613 | ||
2614 | // weight the average with the regularization parameter | |
2615 | hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average); | |
2616 | } | |
2617 | ||
2618 | for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) | |
2619 | hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m)); | |
2620 | ||
2621 | // fill guess | |
2622 | for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++) | |
2623 | { | |
2624 | fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t)); | |
2625 | fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO | |
2626 | ||
2627 | //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t)); | |
2628 | } | |
2629 | ||
2630 | ||
2631 | // calculate chi2 (change from last iteration) | |
2632 | Double_t chi2 = 0; | |
2633 | ||
2634 | // use smoothed true (normalized) as new prior | |
745d6088 | 2635 | norm = 1; //hTrueSmoothed->Integral(); |
cfc19dd5 | 2636 | |
2637 | for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++) | |
2638 | { | |
dd701109 | 2639 | Float_t newValue = hTemp->GetBinContent(t)/norm; |
cfc19dd5 | 2640 | Float_t diff = hPrior->GetBinContent(t) - newValue; |
2641 | chi2 += (Double_t) diff * diff; | |
2642 | ||
2643 | hPrior->SetBinContent(t, newValue); | |
2644 | } | |
2645 | ||
2646 | printf("Chi2 of %d iteration = %.10f\n", i, chi2); | |
2647 | ||
dd701109 | 2648 | //if (i % 5 == 0) |
cfc19dd5 | 2649 | DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY()); |
2650 | ||
2651 | delete hTrueSmoothed; | |
2652 | } // end of iterations | |
2653 | ||
2654 | DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY()); | |
2655 | } | |
2656 | ||
9ca6be09 | 2657 | //____________________________________________________________________ |
2658 | void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace) | |
2659 | { | |
2660 | // | |
2661 | // correct spectrum using a simple Gaussian approach, that is model-dependent | |
2662 | // | |
2663 | ||
2664 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
2665 | Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); | |
2666 | ||
447c325d | 2667 | SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE); |
0b4bfd98 | 2668 | //normalize ESD |
2669 | fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); | |
9ca6be09 | 2670 | |
9ca6be09 | 2671 | TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean")); |
2672 | correction->SetTitle("GaussianMean"); | |
2673 | ||
2674 | TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth")); | |
2675 | correction->SetTitle("GaussianWidth"); | |
2676 | ||
2677 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
2678 | { | |
2679 | TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1); | |
2680 | proj->Fit("gaus", "0Q"); | |
2681 | correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1)); | |
2682 | correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2)); | |
2683 | ||
2684 | continue; | |
2685 | ||
2686 | // draw for debugging | |
2687 | new TCanvas; | |
2688 | proj->DrawCopy(); | |
2689 | proj->GetFunction("gaus")->DrawCopy("SAME"); | |
2690 | } | |
2691 | ||
2692 | TH1* target = fMultiplicityESDCorrected[correlationID]; | |
2693 | ||
2694 | Int_t nBins = target->GetNbinsX()*10+1; | |
2695 | Float_t* binning = new Float_t[nBins]; | |
2696 | for (Int_t i=1; i<=target->GetNbinsX(); ++i) | |
2697 | for (Int_t j=0; j<10; ++j) | |
2698 | binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j; | |
2699 | ||
2700 | binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX()); | |
2701 | ||
2702 | TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning); | |
2703 | ||
2704 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
2705 | { | |
2706 | Float_t mean = correction->GetBinContent(i); | |
2707 | Float_t width = correctionWidth->GetBinContent(i); | |
2708 | ||
3602328d | 2709 | Int_t fillBegin = fineBinned->FindBin(mean - width * 5); |
2710 | Int_t fillEnd = fineBinned->FindBin(mean + width * 5); | |
2711 | //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd); | |
9ca6be09 | 2712 | |
2713 | for (Int_t j=fillBegin; j <= fillEnd; ++j) | |
2714 | { | |
2715 | fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i)); | |
2716 | } | |
2717 | } | |
2718 | ||
2719 | for (Int_t i=1; i<=target->GetNbinsX(); ++i) | |
2720 | { | |
2721 | Float_t sum = 0; | |
2722 | for (Int_t j=1; j<=10; ++j) | |
2723 | sum += fineBinned->GetBinContent((i-1)*10 + j); | |
2724 | target->SetBinContent(i, sum / 10); | |
2725 | } | |
2726 | ||
2727 | delete[] binning; | |
2728 | ||
cfc19dd5 | 2729 | DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY()); |
0a173978 | 2730 | } |
3602328d | 2731 | |
2732 | //____________________________________________________________________ | |
447c325d | 2733 | TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap) |
3602328d | 2734 | { |
2735 | // runs the distribution given in inputMC through the response matrix identified by | |
2736 | // correlationMap and produces a measured distribution | |
2737 | // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex | |
b477f4f2 | 2738 | // if normalized is set, inputMC is expected to be normalized to the bin width |
3602328d | 2739 | |
2740 | if (!inputMC) | |
2741 | return 0; | |
2742 | ||
2743 | if (correlationMap < 0 || correlationMap >= kCorrHists) | |
2744 | return 0; | |
2745 | ||
2746 | // empty under/overflow bins in x, otherwise Project3D takes them into account | |
2747 | TH3* hist = fCorrelation[correlationMap]; | |
039e265e | 2748 | for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y) |
3602328d | 2749 | { |
039e265e | 2750 | for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z) |
3602328d | 2751 | { |
2752 | hist->SetBinContent(0, y, z, 0); | |
2753 | hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); | |
2754 | } | |
2755 | } | |
2756 | ||
447c325d | 2757 | TH2* corr = (TH2*) hist->Project3D("zy"); |
2758 | //corr->Rebin2D(2, 1); | |
3602328d | 2759 | corr->Sumw2(); |
039e265e | 2760 | Printf("Correction histogram used for convolution has %f entries", corr->Integral()); |
3602328d | 2761 | |
2762 | // normalize correction for given nPart | |
2763 | for (Int_t i=1; i<=corr->GetNbinsX(); ++i) | |
2764 | { | |
2765 | Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY()); | |
2766 | if (sum <= 0) | |
2767 | continue; | |
2768 | ||
2769 | for (Int_t j=1; j<=corr->GetNbinsY(); ++j) | |
2770 | { | |
2771 | // npart sum to 1 | |
2772 | corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum); | |
2773 | corr->SetBinError(i, j, corr->GetBinError(i, j) / sum); | |
2774 | } | |
2775 | } | |
2776 | ||
2777 | TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName()))); | |
2778 | target->Reset(); | |
2779 | ||
2780 | for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas) | |
2781 | { | |
2782 | Float_t measured = 0; | |
2783 | Float_t error = 0; | |
2784 | ||
447c325d | 2785 | for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen) |
3602328d | 2786 | { |
447c325d | 2787 | Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen)); |
b477f4f2 | 2788 | |
447c325d | 2789 | measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas); |
2790 | error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas); | |
3602328d | 2791 | } |
2792 | ||
cfc19dd5 | 2793 | //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error); |
3602328d | 2794 | |
447c325d | 2795 | target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured); |
2796 | target->SetBinError(1 + target->GetNbinsX() / 2, meas, error); | |
3602328d | 2797 | } |
2798 | ||
2799 | return target; | |
2800 | } | |
2801 | ||
2802 | //____________________________________________________________________ | |
2803 | void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id) | |
2804 | { | |
2805 | // uses the given function to fill the input MC histogram and generates from that | |
2806 | // the measured histogram by applying the response matrix | |
2807 | // this can be used to evaluate if the methods work indepedently of the input | |
2808 | // distribution | |
2809 | // WARNING does not respect the vertex distribution, just fills central vertex bin | |
2810 | ||
2811 | if (!inputMC) | |
2812 | return; | |
2813 | ||
2814 | if (id < 0 || id >= kESDHists) | |
2815 | return; | |
2816 | ||
dd701109 | 2817 | TH2F* mc = fMultiplicityVtx[id]; |
3602328d | 2818 | |
2819 | mc->Reset(); | |
2820 | ||
2821 | for (Int_t i=1; i<=mc->GetNbinsY(); ++i) | |
b477f4f2 | 2822 | { |
0b4bfd98 | 2823 | mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i)); |
2824 | mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, 0); | |
b477f4f2 | 2825 | } |
3602328d | 2826 | |
2827 | //new TCanvas; | |
2828 | //mc->Draw("COLZ"); | |
2829 | ||
2830 | TH1* proj = mc->ProjectionY(); | |
2831 | ||
dd701109 | 2832 | TString tmp(fMultiplicityESD[id]->GetName()); |
2833 | delete fMultiplicityESD[id]; | |
3602328d | 2834 | fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id); |
dd701109 | 2835 | fMultiplicityESD[id]->SetName(tmp); |
3602328d | 2836 | } |