]>
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. | |
21 | // | |
22 | // Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch | |
23 | ||
24 | #include "AliMultiplicityCorrection.h" | |
25 | ||
26 | #include <TFile.h> | |
27 | #include <TH1F.h> | |
28 | #include <TH2F.h> | |
29 | #include <TH3F.h> | |
30 | #include <TDirectory.h> | |
31 | #include <TVirtualFitter.h> | |
32 | #include <TCanvas.h> | |
33 | #include <TString.h> | |
9ca6be09 | 34 | #include <TF1.h> |
d560b581 | 35 | #include <TMath.h> |
3602328d | 36 | #include <TCollection.h> |
447c325d | 37 | #include <TLegend.h> |
38 | #include <TLine.h> | |
0a173978 | 39 | |
40 | ClassImp(AliMultiplicityCorrection) | |
41 | ||
447c325d | 42 | const Int_t AliMultiplicityCorrection::fgMaxInput = 250; // bins in measured histogram |
43 | const Int_t AliMultiplicityCorrection::fgMaxParams = 250; // bins in unfolded histogram = number of fit params | |
44 | ||
9ca6be09 | 45 | TH1* AliMultiplicityCorrection::fCurrentESD = 0; |
46 | TH1* AliMultiplicityCorrection::fCurrentCorrelation = 0; | |
cfc19dd5 | 47 | TH1* AliMultiplicityCorrection::fCurrentEfficiency = 0; |
447c325d | 48 | TMatrixD* AliMultiplicityCorrection::fCorrelationMatrix = 0; |
49 | TMatrixD* AliMultiplicityCorrection::fCorrelationCovarianceMatrix = 0; | |
50 | TVectorD* AliMultiplicityCorrection::fCurrentESDVector = 0; | |
51 | TVectorD* AliMultiplicityCorrection::fEntropyAPriori = 0; | |
cfc19dd5 | 52 | AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fRegularizationType = AliMultiplicityCorrection::kPol1; |
447c325d | 53 | Float_t AliMultiplicityCorrection::fRegularizationWeight = 5000; |
dd701109 | 54 | TF1* AliMultiplicityCorrection::fNBD = 0; |
0a173978 | 55 | |
56 | //____________________________________________________________________ | |
57 | AliMultiplicityCorrection::AliMultiplicityCorrection() : | |
dd701109 | 58 | TNamed(), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0) |
0a173978 | 59 | { |
60 | // | |
61 | // default constructor | |
62 | // | |
63 | ||
64 | for (Int_t i = 0; i < kESDHists; ++i) | |
65 | fMultiplicityESD[i] = 0; | |
66 | ||
67 | for (Int_t i = 0; i < kMCHists; ++i) | |
cfc19dd5 | 68 | { |
69 | fMultiplicityVtx[i] = 0; | |
70 | fMultiplicityMB[i] = 0; | |
71 | fMultiplicityINEL[i] = 0; | |
72 | } | |
0a173978 | 73 | |
74 | for (Int_t i = 0; i < kCorrHists; ++i) | |
75 | { | |
76 | fCorrelation[i] = 0; | |
77 | fMultiplicityESDCorrected[i] = 0; | |
78 | } | |
79 | } | |
80 | ||
81 | //____________________________________________________________________ | |
82 | AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) : | |
dd701109 | 83 | TNamed(name, title), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0) |
0a173978 | 84 | { |
85 | // | |
86 | // named constructor | |
87 | // | |
88 | ||
89 | // do not add this hists to the directory | |
90 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
91 | TH1::AddDirectory(kFALSE); | |
92 | ||
b477f4f2 | 93 | /*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 | 94 | 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, |
95 | 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, | |
96 | 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, | |
97 | 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, | |
98 | 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, | |
99 | 50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5, | |
100 | 100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5, | |
101 | 150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5, | |
102 | 250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5, | |
9ca6be09 | 103 | 500.5 }; |
104 | //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5, | |
105 | //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5, | |
106 | //1000.5 }; | |
0a173978 | 107 | |
b477f4f2 | 108 | #define VTXBINNING 10, binLimitsVtx |
109 | #define NBINNING fgMaxParams, binLimitsN*/ | |
110 | ||
447c325d | 111 | #define NBINNING 501, -0.5, 500.5 |
112 | #define VTXBINNING 1, -10, 10 | |
0a173978 | 113 | |
114 | for (Int_t i = 0; i < kESDHists; ++i) | |
b477f4f2 | 115 | fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING); |
0a173978 | 116 | |
117 | for (Int_t i = 0; i < kMCHists; ++i) | |
118 | { | |
cfc19dd5 | 119 | fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i))); |
120 | fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart"); | |
121 | ||
122 | fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i))); | |
123 | fMultiplicityMB[i]->SetTitle("fMultiplicityMB"); | |
124 | ||
125 | fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i))); | |
126 | fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL"); | |
0a173978 | 127 | } |
128 | ||
129 | for (Int_t i = 0; i < kCorrHists; ++i) | |
130 | { | |
b477f4f2 | 131 | fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING); |
132 | fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING); | |
0a173978 | 133 | } |
134 | ||
135 | TH1::AddDirectory(oldStatus); | |
136 | } | |
137 | ||
138 | //____________________________________________________________________ | |
139 | AliMultiplicityCorrection::~AliMultiplicityCorrection() | |
140 | { | |
141 | // | |
142 | // Destructor | |
143 | // | |
144 | } | |
145 | ||
146 | //____________________________________________________________________ | |
147 | Long64_t AliMultiplicityCorrection::Merge(TCollection* list) | |
148 | { | |
149 | // Merge a list of AliMultiplicityCorrection objects with this (needed for | |
150 | // PROOF). | |
151 | // Returns the number of merged objects (including this). | |
152 | ||
153 | if (!list) | |
154 | return 0; | |
155 | ||
156 | if (list->IsEmpty()) | |
157 | return 1; | |
158 | ||
159 | TIterator* iter = list->MakeIterator(); | |
160 | TObject* obj; | |
161 | ||
162 | // collections of all histograms | |
cfc19dd5 | 163 | TList collections[kESDHists+kMCHists*3+kCorrHists*2]; |
0a173978 | 164 | |
165 | Int_t count = 0; | |
166 | while ((obj = iter->Next())) { | |
167 | ||
168 | AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj); | |
169 | if (entry == 0) | |
170 | continue; | |
171 | ||
172 | for (Int_t i = 0; i < kESDHists; ++i) | |
173 | collections[i].Add(entry->fMultiplicityESD[i]); | |
174 | ||
175 | for (Int_t i = 0; i < kMCHists; ++i) | |
cfc19dd5 | 176 | { |
177 | collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]); | |
178 | collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]); | |
179 | collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]); | |
180 | } | |
0a173978 | 181 | |
182 | for (Int_t i = 0; i < kCorrHists; ++i) | |
cfc19dd5 | 183 | collections[kESDHists+kMCHists*3+i].Add(entry->fCorrelation[i]); |
0a173978 | 184 | |
185 | for (Int_t i = 0; i < kCorrHists; ++i) | |
cfc19dd5 | 186 | collections[kESDHists+kMCHists*3+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]); |
0a173978 | 187 | |
188 | count++; | |
189 | } | |
190 | ||
191 | for (Int_t i = 0; i < kESDHists; ++i) | |
192 | fMultiplicityESD[i]->Merge(&collections[i]); | |
193 | ||
194 | for (Int_t i = 0; i < kMCHists; ++i) | |
cfc19dd5 | 195 | { |
196 | fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]); | |
197 | fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]); | |
198 | fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]); | |
199 | } | |
0a173978 | 200 | |
201 | for (Int_t i = 0; i < kCorrHists; ++i) | |
cfc19dd5 | 202 | fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*3+i]); |
0a173978 | 203 | |
204 | for (Int_t i = 0; i < kCorrHists; ++i) | |
cfc19dd5 | 205 | fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*3+kCorrHists+i]); |
0a173978 | 206 | |
207 | delete iter; | |
208 | ||
209 | return count+1; | |
210 | } | |
211 | ||
212 | //____________________________________________________________________ | |
213 | Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir) | |
214 | { | |
215 | // | |
216 | // loads the histograms from a file | |
217 | // if dir is empty a directory with the name of this object is taken (like in SaveHistogram) | |
218 | // | |
219 | ||
220 | if (!dir) | |
221 | dir = GetName(); | |
222 | ||
223 | if (!gDirectory->cd(dir)) | |
224 | return kFALSE; | |
225 | ||
226 | // TODO memory leak. old histograms needs to be deleted. | |
227 | ||
228 | Bool_t success = kTRUE; | |
229 | ||
230 | for (Int_t i = 0; i < kESDHists; ++i) | |
231 | { | |
232 | fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName())); | |
233 | if (!fMultiplicityESD[i]) | |
234 | success = kFALSE; | |
235 | } | |
236 | ||
237 | for (Int_t i = 0; i < kMCHists; ++i) | |
238 | { | |
cfc19dd5 | 239 | fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName())); |
240 | fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName())); | |
241 | fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName())); | |
242 | if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i]) | |
0a173978 | 243 | success = kFALSE; |
244 | } | |
245 | ||
246 | for (Int_t i = 0; i < kCorrHists; ++i) | |
247 | { | |
248 | fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName())); | |
249 | if (!fCorrelation[i]) | |
250 | success = kFALSE; | |
251 | fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName())); | |
252 | if (!fMultiplicityESDCorrected[i]) | |
253 | success = kFALSE; | |
254 | } | |
255 | ||
256 | gDirectory->cd(".."); | |
257 | ||
258 | return success; | |
259 | } | |
260 | ||
261 | //____________________________________________________________________ | |
262 | void AliMultiplicityCorrection::SaveHistograms() | |
263 | { | |
264 | // | |
265 | // saves the histograms | |
266 | // | |
267 | ||
268 | gDirectory->mkdir(GetName()); | |
269 | gDirectory->cd(GetName()); | |
270 | ||
271 | for (Int_t i = 0; i < kESDHists; ++i) | |
272 | if (fMultiplicityESD[i]) | |
273 | fMultiplicityESD[i]->Write(); | |
274 | ||
275 | for (Int_t i = 0; i < kMCHists; ++i) | |
cfc19dd5 | 276 | { |
277 | if (fMultiplicityVtx[i]) | |
278 | fMultiplicityVtx[i]->Write(); | |
279 | if (fMultiplicityMB[i]) | |
280 | fMultiplicityMB[i]->Write(); | |
281 | if (fMultiplicityINEL[i]) | |
282 | fMultiplicityINEL[i]->Write(); | |
283 | } | |
0a173978 | 284 | |
285 | for (Int_t i = 0; i < kCorrHists; ++i) | |
286 | { | |
287 | if (fCorrelation[i]) | |
288 | fCorrelation[i]->Write(); | |
289 | if (fMultiplicityESDCorrected[i]) | |
290 | fMultiplicityESDCorrected[i]->Write(); | |
291 | } | |
292 | ||
293 | gDirectory->cd(".."); | |
294 | } | |
295 | ||
296 | //____________________________________________________________________ | |
cfc19dd5 | 297 | 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 | 298 | { |
299 | // | |
300 | // Fills an event from MC | |
301 | // | |
302 | ||
cfc19dd5 | 303 | if (triggered) |
304 | { | |
305 | fMultiplicityMB[0]->Fill(vtx, generated05); | |
306 | fMultiplicityMB[1]->Fill(vtx, generated10); | |
307 | fMultiplicityMB[2]->Fill(vtx, generated15); | |
308 | fMultiplicityMB[3]->Fill(vtx, generated20); | |
309 | fMultiplicityMB[4]->Fill(vtx, generatedAll); | |
310 | ||
311 | if (vertex) | |
312 | { | |
313 | fMultiplicityVtx[0]->Fill(vtx, generated05); | |
314 | fMultiplicityVtx[1]->Fill(vtx, generated10); | |
315 | fMultiplicityVtx[2]->Fill(vtx, generated15); | |
316 | fMultiplicityVtx[3]->Fill(vtx, generated20); | |
317 | fMultiplicityVtx[4]->Fill(vtx, generatedAll); | |
318 | } | |
319 | } | |
320 | ||
321 | fMultiplicityINEL[0]->Fill(vtx, generated05); | |
322 | fMultiplicityINEL[1]->Fill(vtx, generated10); | |
323 | fMultiplicityINEL[2]->Fill(vtx, generated15); | |
324 | fMultiplicityINEL[3]->Fill(vtx, generated20); | |
325 | fMultiplicityINEL[4]->Fill(vtx, generatedAll); | |
0a173978 | 326 | } |
327 | ||
328 | //____________________________________________________________________ | |
329 | void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20) | |
330 | { | |
331 | // | |
332 | // Fills an event from ESD | |
333 | // | |
334 | ||
335 | fMultiplicityESD[0]->Fill(vtx, measured05); | |
336 | fMultiplicityESD[1]->Fill(vtx, measured10); | |
337 | fMultiplicityESD[2]->Fill(vtx, measured15); | |
338 | fMultiplicityESD[3]->Fill(vtx, measured20); | |
339 | } | |
340 | ||
341 | //____________________________________________________________________ | |
342 | 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) | |
343 | { | |
344 | // | |
345 | // Fills an event into the correlation map with the information from MC and ESD | |
346 | // | |
347 | ||
348 | fCorrelation[0]->Fill(vtx, generated05, measured05); | |
349 | fCorrelation[1]->Fill(vtx, generated10, measured10); | |
350 | fCorrelation[2]->Fill(vtx, generated15, measured15); | |
351 | fCorrelation[3]->Fill(vtx, generated20, measured20); | |
352 | ||
353 | fCorrelation[4]->Fill(vtx, generatedAll, measured05); | |
354 | fCorrelation[5]->Fill(vtx, generatedAll, measured10); | |
355 | fCorrelation[6]->Fill(vtx, generatedAll, measured15); | |
356 | fCorrelation[7]->Fill(vtx, generatedAll, measured20); | |
357 | } | |
358 | ||
359 | //____________________________________________________________________ | |
447c325d | 360 | Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params) |
0a173978 | 361 | { |
362 | // homogenity term for minuit fitting | |
363 | // pure function of the parameters | |
364 | // prefers constant function (pol0) | |
365 | ||
366 | Double_t chi2 = 0; | |
367 | ||
447c325d | 368 | // ignore the first bin here. on purpose... |
369 | for (Int_t i=2; i<fgMaxParams; ++i) | |
0a173978 | 370 | { |
dd701109 | 371 | Double_t right = params[i]; |
372 | Double_t left = params[i-1]; | |
0a173978 | 373 | |
447c325d | 374 | if (right != 0) |
375 | { | |
376 | Double_t diff = 1 - left / right; | |
377 | chi2 += diff * diff; | |
378 | } | |
0a173978 | 379 | } |
380 | ||
447c325d | 381 | return chi2 / 100.0; |
0a173978 | 382 | } |
383 | ||
384 | //____________________________________________________________________ | |
447c325d | 385 | Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params) |
0a173978 | 386 | { |
387 | // homogenity term for minuit fitting | |
388 | // pure function of the parameters | |
389 | // prefers linear function (pol1) | |
390 | ||
391 | Double_t chi2 = 0; | |
392 | ||
447c325d | 393 | // ignore the first bin here. on purpose... |
394 | for (Int_t i=2+1; i<fgMaxParams; ++i) | |
0a173978 | 395 | { |
dd701109 | 396 | if (params[i-1] == 0) |
0a173978 | 397 | continue; |
398 | ||
dd701109 | 399 | Double_t right = params[i]; |
400 | Double_t middle = params[i-1]; | |
401 | Double_t left = params[i-2]; | |
402 | ||
403 | Double_t der1 = (right - middle); | |
404 | Double_t der2 = (middle - left); | |
405 | ||
406 | Double_t diff = (der1 - der2) / middle; | |
407 | ||
408 | chi2 += diff * diff; | |
409 | } | |
0a173978 | 410 | |
dd701109 | 411 | return chi2; |
412 | } | |
0a173978 | 413 | |
dd701109 | 414 | //____________________________________________________________________ |
447c325d | 415 | Double_t AliMultiplicityCorrection::RegularizationTest(TVectorD& params) |
dd701109 | 416 | { |
417 | // homogenity term for minuit fitting | |
418 | // pure function of the parameters | |
419 | // prefers linear function (pol1) | |
0a173978 | 420 | |
dd701109 | 421 | Double_t chi2 = 0; |
0a173978 | 422 | |
dd701109 | 423 | Float_t der[fgMaxParams]; |
424 | ||
425 | for (Int_t i=0; i<fgMaxParams-1; ++i) | |
426 | der[i] = params[i+1] - params[i]; | |
427 | ||
428 | for (Int_t i=0; i<fgMaxParams-6; ++i) | |
429 | { | |
430 | for (Int_t j=1; j<=5; ++j) | |
431 | { | |
432 | Double_t diff = der[i+j] - der[i]; | |
433 | chi2 += diff * diff; | |
434 | } | |
0a173978 | 435 | } |
436 | ||
dd701109 | 437 | return chi2 * 100; // 10000 |
0a173978 | 438 | } |
439 | ||
9ca6be09 | 440 | //____________________________________________________________________ |
447c325d | 441 | Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& params) |
9ca6be09 | 442 | { |
443 | // homogenity term for minuit fitting | |
444 | // pure function of the parameters | |
445 | // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments, | |
446 | // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp. | |
447 | ||
448 | Double_t chi2 = 0; | |
449 | ||
447c325d | 450 | // ignore the first bin here. on purpose... |
451 | for (Int_t i=3; i<fgMaxParams; ++i) | |
9ca6be09 | 452 | { |
dd701109 | 453 | Double_t right = params[i]; |
454 | Double_t middle = params[i-1]; | |
455 | Double_t left = params[i-2]; | |
9ca6be09 | 456 | |
dd701109 | 457 | Double_t der1 = (right - middle); |
458 | Double_t der2 = (middle - left); | |
9ca6be09 | 459 | |
447c325d | 460 | Double_t diff = (der1 - der2); |
9ca6be09 | 461 | |
447c325d | 462 | chi2 += diff * diff; |
9ca6be09 | 463 | } |
464 | ||
447c325d | 465 | return chi2 * 1e4; |
3602328d | 466 | } |
467 | ||
468 | //____________________________________________________________________ | |
447c325d | 469 | Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params) |
3602328d | 470 | { |
471 | // homogenity term for minuit fitting | |
472 | // pure function of the parameters | |
473 | // calculates entropy, from | |
474 | // The method of reduced cross-entropy (M. Schmelling 1993) | |
475 | ||
3602328d | 476 | Double_t paramSum = 0; |
447c325d | 477 | // ignore the first bin here. on purpose... |
478 | for (Int_t i=1; i<fgMaxParams; ++i) | |
dd701109 | 479 | paramSum += params[i]; |
3602328d | 480 | |
481 | Double_t chi2 = 0; | |
447c325d | 482 | for (Int_t i=1; i<fgMaxParams; ++i) |
3602328d | 483 | { |
dd701109 | 484 | Double_t tmp = params[i] / paramSum; |
447c325d | 485 | if (tmp > 0 && (*fEntropyAPriori)[i] > 0) |
486 | chi2 += tmp * log(tmp / (*fEntropyAPriori)[i]); | |
3602328d | 487 | } |
3602328d | 488 | |
447c325d | 489 | return 10.0 + chi2; |
dd701109 | 490 | } |
491 | ||
492 | //____________________________________________________________________ | |
493 | void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3) | |
494 | { | |
495 | // | |
496 | // fit function for minuit | |
497 | // does: nbd | |
498 | // 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); | |
499 | // func->SetParNames("scaling", "averagen", "k"); | |
447c325d | 500 | // func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000); |
501 | // func->SetParLimits(1, 0.001, 1000); | |
502 | // func->SetParLimits(2, 0.001, 1000); | |
dd701109 | 503 | // |
504 | ||
505 | fNBD->SetParameters(params[0], params[1], params[2]); | |
506 | ||
507 | Double_t params2[fgMaxParams]; | |
508 | ||
509 | for (Int_t i=0; i<fgMaxParams; ++i) | |
510 | params2[i] = fNBD->Eval(i); | |
511 | ||
512 | MinuitFitFunction(unused1, unused2, chi2, params2, unused3); | |
513 | ||
514 | printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2); | |
9ca6be09 | 515 | } |
516 | ||
0a173978 | 517 | //____________________________________________________________________ |
518 | void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t) | |
519 | { | |
520 | // | |
521 | // fit function for minuit | |
cfc19dd5 | 522 | // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix |
0a173978 | 523 | // |
524 | ||
cfc19dd5 | 525 | // d |
447c325d | 526 | TVectorD paramsVector(fgMaxParams); |
cfc19dd5 | 527 | for (Int_t i=0; i<fgMaxParams; ++i) |
447c325d | 528 | paramsVector[i] = params[i] * params[i]; |
529 | ||
530 | // calculate penalty factor | |
531 | Double_t penaltyVal = 0; | |
532 | switch (fRegularizationType) | |
533 | { | |
534 | case kNone: break; | |
535 | case kPol0: penaltyVal = RegularizationPol0(paramsVector); break; | |
536 | case kPol1: penaltyVal = RegularizationPol1(paramsVector); break; | |
537 | case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break; | |
538 | case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break; | |
539 | case kTest: penaltyVal = RegularizationTest(paramsVector); break; | |
540 | } | |
541 | ||
542 | //if (penaltyVal > 1e10) | |
543 | // paramsVector2.Print(); | |
544 | ||
545 | penaltyVal *= fRegularizationWeight; | |
cfc19dd5 | 546 | |
547 | // Ad | |
447c325d | 548 | TVectorD measGuessVector(fgMaxInput); |
549 | measGuessVector = (*fCorrelationMatrix) * paramsVector; | |
cfc19dd5 | 550 | |
551 | // Ad - m | |
447c325d | 552 | measGuessVector -= (*fCurrentESDVector); |
cfc19dd5 | 553 | |
447c325d | 554 | TVectorD copy(measGuessVector); |
cfc19dd5 | 555 | |
556 | // (Ad - m) W | |
447c325d | 557 | measGuessVector *= (*fCorrelationCovarianceMatrix); |
cfc19dd5 | 558 | |
447c325d | 559 | //measGuessVector.Print(); |
cfc19dd5 | 560 | |
447c325d | 561 | // (Ad - m) W (Ad - m) |
562 | Double_t chi2FromFit = measGuessVector * copy * 1e6; | |
0a173978 | 563 | |
3602328d | 564 | chi2 = chi2FromFit + penaltyVal; |
0a173978 | 565 | |
447c325d | 566 | static Int_t callCount = 0; |
567 | if ((callCount++ % 10000) == 0) | |
568 | printf("%d) %f %f --> %f\n", callCount, chi2FromFit, penaltyVal, chi2); | |
0a173978 | 569 | } |
570 | ||
571 | //____________________________________________________________________ | |
447c325d | 572 | void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin) |
0a173978 | 573 | { |
574 | // | |
9ca6be09 | 575 | // fills fCurrentESD, fCurrentCorrelation |
576 | // resets fMultiplicityESDCorrected | |
0a173978 | 577 | // |
578 | ||
579 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
9ca6be09 | 580 | fMultiplicityESDCorrected[correlationID]->Reset(); |
0a173978 | 581 | |
9ca6be09 | 582 | fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY(); |
583 | fCurrentESD->Sumw2(); | |
0a173978 | 584 | |
585 | // empty under/overflow bins in x, otherwise Project3D takes them into account | |
586 | TH3* hist = fCorrelation[correlationID]; | |
587 | for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y) | |
588 | { | |
589 | for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z) | |
590 | { | |
591 | hist->SetBinContent(0, y, z, 0); | |
592 | hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); | |
593 | } | |
594 | } | |
595 | ||
9ca6be09 | 596 | fCurrentCorrelation = hist->Project3D("zy"); |
447c325d | 597 | //((TH2*) fCurrentCorrelation)->Rebin2D(2, 1); |
598 | //fMultiplicityESDCorrected[correlationID]->Rebin(2); | |
9ca6be09 | 599 | fCurrentCorrelation->Sumw2(); |
cfc19dd5 | 600 | |
447c325d | 601 | if (createBigBin) |
602 | { | |
603 | Int_t maxBin = 0; | |
604 | for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i) | |
605 | { | |
606 | if (fCurrentESD->GetBinContent(i) <= 5) | |
607 | { | |
608 | maxBin = i; | |
609 | break; | |
610 | } | |
611 | } | |
612 | ||
613 | if (maxBin > 0) | |
614 | { | |
615 | TCanvas* canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800); | |
616 | canvas->Divide(2, 2); | |
617 | ||
618 | canvas->cd(1); | |
619 | fCurrentESD->DrawCopy(); | |
620 | gPad->SetLogy(); | |
621 | ||
622 | canvas->cd(2); | |
623 | fCurrentCorrelation->DrawCopy("COLZ"); | |
624 | ||
625 | printf("Bin limit in measured spectrum is %d.\n", maxBin); | |
626 | fCurrentESD->SetBinContent(maxBin, fCurrentESD->Integral(maxBin, fCurrentESD->GetNbinsX())); | |
627 | for (Int_t i=maxBin+1; i<=fCurrentESD->GetNbinsX(); ++i) | |
628 | { | |
629 | fCurrentESD->SetBinContent(i, 0); | |
630 | fCurrentESD->SetBinError(i, 0); | |
631 | } | |
632 | // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? | |
633 | fCurrentESD->SetBinError(maxBin, TMath::Sqrt(fCurrentESD->GetBinContent(maxBin))); | |
634 | ||
635 | printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(maxBin), fCurrentESD->GetBinError(maxBin)); | |
636 | ||
637 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
638 | { | |
639 | fCurrentCorrelation->SetBinContent(i, maxBin, fCurrentCorrelation->Integral(i, i, maxBin, fCurrentCorrelation->GetNbinsY())); | |
640 | // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? | |
641 | fCurrentCorrelation->SetBinError(i, maxBin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, maxBin))); | |
642 | ||
643 | for (Int_t j=maxBin+1; j<=fCurrentCorrelation->GetNbinsY(); ++j) | |
644 | { | |
645 | fCurrentCorrelation->SetBinContent(i, j, 0); | |
646 | fCurrentCorrelation->SetBinError(i, j, 0); | |
647 | } | |
648 | } | |
649 | ||
650 | printf("Adjusted correlation matrix!\n"); | |
651 | ||
652 | canvas->cd(3); | |
653 | fCurrentESD->DrawCopy(); | |
654 | gPad->SetLogy(); | |
655 | ||
656 | canvas->cd(4); | |
657 | fCurrentCorrelation->DrawCopy("COLZ"); | |
658 | } | |
659 | } | |
660 | ||
661 | //normalize ESD | |
662 | fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); | |
663 | ||
cfc19dd5 | 664 | fCurrentEfficiency = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency"); |
665 | TH2* divisor = 0; | |
666 | switch (eventType) | |
667 | { | |
668 | case kTrVtx : divisor = fMultiplicityVtx[inputRange]; break; | |
669 | case kMB: divisor = fMultiplicityMB[inputRange]; break; | |
670 | case kINEL: divisor = fMultiplicityINEL[inputRange]; break; | |
671 | } | |
672 | fCurrentEfficiency->Divide(divisor->ProjectionY()); | |
447c325d | 673 | //fCurrentEfficiency->Rebin(2); |
674 | //fCurrentEfficiency->Scale(0.5); | |
9ca6be09 | 675 | } |
676 | ||
677 | //____________________________________________________________________ | |
447c325d | 678 | void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight) |
679 | { | |
680 | fRegularizationType = type; | |
681 | fRegularizationWeight = weight; | |
682 | ||
683 | printf("AliMultiplicityCorrection::SetRegularizationParameters --> Regularization set to %d with weight %f\n", (Int_t) type, weight); | |
684 | } | |
685 | ||
686 | //____________________________________________________________________ | |
687 | Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* inputDist) | |
9ca6be09 | 688 | { |
689 | // | |
690 | // correct spectrum using minuit chi2 method | |
691 | // | |
692 | // if check is kTRUE the input MC solution (by definition the right solution) is used | |
693 | // no fit is made and just the chi2 is calculated | |
694 | // | |
695 | ||
696 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
697 | Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); | |
698 | ||
447c325d | 699 | SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE); // TODO FAKE kTRUE |
cfc19dd5 | 700 | |
447c325d | 701 | fCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams); |
702 | fCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput); | |
703 | fCurrentESDVector = new TVectorD(fgMaxInput); | |
704 | fEntropyAPriori = new TVectorD(fgMaxParams); | |
705 | ||
706 | /*new TCanvas; fCurrentESD->DrawCopy(); | |
707 | fCurrentESD = ((TH2*) fCurrentCorrelation)->ProjectionY("check-proj2"); | |
708 | fCurrentESD->Sumw2(); | |
709 | fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); | |
710 | fCurrentESD->SetLineColor(2); | |
711 | fCurrentESD->DrawCopy("SAME");*/ | |
0a173978 | 712 | |
713 | // normalize correction for given nPart | |
9ca6be09 | 714 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) |
0a173978 | 715 | { |
9ca6be09 | 716 | Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); |
0a173978 | 717 | if (sum <= 0) |
718 | continue; | |
9ca6be09 | 719 | Float_t maxValue = 0; |
720 | Int_t maxBin = -1; | |
721 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) | |
0a173978 | 722 | { |
9ca6be09 | 723 | // find most probably value |
724 | if (maxValue < fCurrentCorrelation->GetBinContent(i, j)) | |
725 | { | |
726 | maxValue = fCurrentCorrelation->GetBinContent(i, j); | |
727 | maxBin = j; | |
728 | } | |
729 | ||
0a173978 | 730 | // npart sum to 1 |
9ca6be09 | 731 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); |
732 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
cfc19dd5 | 733 | |
447c325d | 734 | if (i <= fgMaxParams && j <= fgMaxInput) |
cfc19dd5 | 735 | (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j); |
0a173978 | 736 | } |
9ca6be09 | 737 | |
cfc19dd5 | 738 | //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin)); |
0a173978 | 739 | } |
740 | ||
0a173978 | 741 | // Initialize TMinuit via generic fitter interface |
742 | TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams); | |
447c325d | 743 | Double_t arglist[100]; |
744 | arglist[0] = 0; | |
745 | minuit->ExecuteCommand("SET PRINT", arglist, 1); | |
0a173978 | 746 | |
747 | minuit->SetFCN(MinuitFitFunction); | |
748 | ||
447c325d | 749 | for (Int_t i=0; i<fgMaxParams; i++) |
750 | (*fEntropyAPriori)[i] = 1; | |
0a173978 | 751 | |
dd701109 | 752 | if (inputDist) |
753 | { | |
754 | printf("Using different starting conditions...\n"); | |
755 | new TCanvas; | |
756 | inputDist->DrawCopy(); | |
447c325d | 757 | |
758 | inputDist->Scale(1.0 / inputDist->Integral()); | |
759 | for (Int_t i=0; i<fgMaxParams; i++) | |
760 | if (inputDist->GetBinContent(i+1) > 0) | |
761 | (*fEntropyAPriori)[i] = inputDist->GetBinContent(i+1); | |
dd701109 | 762 | } |
763 | else | |
764 | inputDist = fCurrentESD; | |
765 | ||
dd701109 | 766 | |
447c325d | 767 | //Float_t minStartValue = 0; //1e-3; |
dd701109 | 768 | |
447c325d | 769 | //new TCanvas; fMultiplicityVtx[mcTarget]->Draw("COLZ"); |
dd701109 | 770 | TH1* proj = fMultiplicityVtx[mcTarget]->ProjectionY("check-proj"); |
447c325d | 771 | //proj->Rebin(2); |
dd701109 | 772 | proj->Scale(1.0 / proj->Integral()); |
447c325d | 773 | |
774 | Double_t results[fgMaxParams]; | |
0a173978 | 775 | for (Int_t i=0; i<fgMaxParams; ++i) |
776 | { | |
dd701109 | 777 | results[i] = inputDist->GetBinContent(i+1); |
447c325d | 778 | |
9ca6be09 | 779 | if (check) |
780 | results[i] = proj->GetBinContent(i+1); | |
dd701109 | 781 | |
447c325d | 782 | // minimum value |
783 | results[i] = TMath::Max(results[i], 1e-3); | |
784 | ||
785 | Float_t stepSize = 0.1; | |
786 | ||
787 | // minuit sees squared values to prevent it from going negative... | |
788 | results[i] = TMath::Sqrt(results[i]); | |
789 | //stepSize /= results[i]; // keep relative step size | |
790 | ||
791 | minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0); | |
792 | } | |
793 | // bin 0 is filled with value from bin 1 (otherwise it's 0) | |
794 | //minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 0); | |
795 | //results[0] = 0; | |
796 | //minuit->SetParameter(0, "param0", 0, 0, 0, 0); | |
797 | ||
798 | for (Int_t i=0; i<fgMaxInput; ++i) | |
799 | { | |
800 | //printf("%f %f %f\n", inputDist->GetBinContent(i+1), TMath::Sqrt(inputDist->GetBinContent(i+1)), inputDist->GetBinError(i+1)); | |
cfc19dd5 | 801 | |
802 | (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1); | |
803 | if (fCurrentESD->GetBinError(i+1) > 0) | |
447c325d | 804 | (*fCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / fCurrentESD->GetBinError(i+1) / fCurrentESD->GetBinError(i+1); |
805 | ||
806 | if ((*fCorrelationCovarianceMatrix)(i, i) > 1e7) | |
807 | (*fCorrelationCovarianceMatrix)(i, i) = 0; | |
cfc19dd5 | 808 | |
447c325d | 809 | //printf("%d --> %e\n", i, (*fCorrelationCovarianceMatrix)(i, i)); |
0a173978 | 810 | } |
811 | ||
812 | Int_t dummy; | |
813 | Double_t chi2; | |
814 | MinuitFitFunction(dummy, 0, chi2, results, 0); | |
dd701109 | 815 | printf("Chi2 of initial parameters is = %f\n", chi2); |
0a173978 | 816 | |
9ca6be09 | 817 | if (check) |
447c325d | 818 | return -1; |
819 | ||
820 | // first param is number of iterations, second is precision.... | |
821 | arglist[0] = 1e6; | |
822 | //arglist[1] = 1e-5; | |
823 | //minuit->ExecuteCommand("SCAN", arglist, 0); | |
824 | Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1); | |
825 | printf("MINUIT status is %d\n", status); | |
826 | //minuit->ExecuteCommand("MIGRAD", arglist, 1); | |
827 | //minuit->ExecuteCommand("MIGRAD", arglist, 1); | |
dd701109 | 828 | //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!"); |
0a173978 | 829 | //minuit->ExecuteCommand("MINOS", arglist, 0); |
830 | ||
831 | for (Int_t i=0; i<fgMaxParams; ++i) | |
832 | { | |
dd701109 | 833 | if (fCurrentEfficiency->GetBinContent(i+1) > 0) |
834 | { | |
447c325d | 835 | fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, minuit->GetParameter(i) * minuit->GetParameter(i) / fCurrentEfficiency->GetBinContent(i+1)); |
836 | // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i)) | |
837 | fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, minuit->GetParError(i) * minuit->GetParameter(i) / fCurrentEfficiency->GetBinContent(i+1)); | |
dd701109 | 838 | } |
839 | else | |
840 | { | |
841 | fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, 0); | |
842 | fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0); | |
843 | } | |
844 | } | |
447c325d | 845 | |
846 | return status; | |
dd701109 | 847 | } |
848 | ||
849 | //____________________________________________________________________ | |
850 | void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace) | |
851 | { | |
852 | // | |
853 | // correct spectrum using minuit chi2 method applying a NBD fit | |
854 | // | |
855 | ||
856 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
857 | ||
447c325d | 858 | SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE); |
dd701109 | 859 | |
447c325d | 860 | fCorrelationMatrix = new TMatrixD(fgMaxParams, fgMaxParams); |
861 | fCorrelationCovarianceMatrix = new TMatrixD(fgMaxParams, fgMaxParams); | |
862 | fCurrentESDVector = new TVectorD(fgMaxParams); | |
dd701109 | 863 | |
864 | // normalize correction for given nPart | |
865 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
866 | { | |
867 | Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); | |
868 | if (sum <= 0) | |
869 | continue; | |
870 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) | |
871 | { | |
872 | // npart sum to 1 | |
873 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); | |
874 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
875 | ||
876 | if (i <= fgMaxParams && j <= fgMaxParams) | |
877 | (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j); | |
878 | } | |
879 | } | |
880 | ||
881 | for (Int_t i=0; i<fgMaxParams; ++i) | |
882 | { | |
883 | (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1); | |
884 | if (fCurrentESD->GetBinError(i+1) > 0) | |
885 | (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1)); | |
886 | } | |
887 | ||
888 | // Create NBD function | |
889 | if (!fNBD) | |
890 | { | |
891 | fNBD = 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); | |
892 | fNBD->SetParNames("scaling", "averagen", "k"); | |
0a173978 | 893 | } |
894 | ||
dd701109 | 895 | // Initialize TMinuit via generic fitter interface |
896 | TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3); | |
897 | ||
898 | minuit->SetFCN(MinuitNBD); | |
899 | minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000); | |
900 | minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000); | |
901 | minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000); | |
902 | ||
903 | Double_t arglist[100]; | |
904 | arglist[0] = 0; | |
905 | minuit->ExecuteCommand("SET PRINT", arglist, 1); | |
447c325d | 906 | minuit->ExecuteCommand("MIGRAD", arglist, 0); |
dd701109 | 907 | //minuit->ExecuteCommand("MINOS", arglist, 0); |
908 | ||
909 | fNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2)); | |
910 | ||
911 | for (Int_t i=0; i<fgMaxParams; ++i) | |
912 | { | |
913 | printf("%d : %f\n", i, fNBD->Eval(i)); | |
914 | if (fNBD->Eval(i) > 0) | |
915 | { | |
916 | fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fNBD->Eval(i)); | |
917 | fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0); | |
918 | } | |
919 | } | |
0a173978 | 920 | } |
921 | ||
922 | //____________________________________________________________________ | |
923 | void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist) | |
924 | { | |
925 | // | |
926 | // normalizes a 1-d histogram to its bin width | |
927 | // | |
928 | ||
929 | for (Int_t i=1; i<=hist->GetNbinsX(); ++i) | |
930 | { | |
931 | hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i)); | |
932 | hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i)); | |
933 | } | |
934 | } | |
935 | ||
936 | //____________________________________________________________________ | |
937 | void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist) | |
938 | { | |
939 | // | |
940 | // normalizes a 2-d histogram to its bin width (x width * y width) | |
941 | // | |
942 | ||
943 | for (Int_t i=1; i<=hist->GetNbinsX(); ++i) | |
944 | for (Int_t j=1; j<=hist->GetNbinsY(); ++j) | |
945 | { | |
946 | Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j); | |
947 | hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor); | |
948 | hist->SetBinError(i, j, hist->GetBinError(i, j) / factor); | |
949 | } | |
950 | } | |
951 | ||
0a173978 | 952 | //____________________________________________________________________ |
953 | void AliMultiplicityCorrection::DrawHistograms() | |
954 | { | |
955 | printf("ESD:\n"); | |
956 | ||
957 | TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600); | |
958 | canvas1->Divide(3, 2); | |
959 | for (Int_t i = 0; i < kESDHists; ++i) | |
960 | { | |
961 | canvas1->cd(i+1); | |
962 | fMultiplicityESD[i]->DrawCopy("COLZ"); | |
963 | printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean()); | |
964 | } | |
965 | ||
cfc19dd5 | 966 | printf("Vtx:\n"); |
0a173978 | 967 | |
968 | TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600); | |
969 | canvas2->Divide(3, 2); | |
970 | for (Int_t i = 0; i < kMCHists; ++i) | |
971 | { | |
972 | canvas2->cd(i+1); | |
cfc19dd5 | 973 | fMultiplicityVtx[i]->DrawCopy("COLZ"); |
974 | printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean()); | |
0a173978 | 975 | } |
976 | ||
977 | TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900); | |
978 | canvas3->Divide(3, 3); | |
979 | for (Int_t i = 0; i < kCorrHists; ++i) | |
980 | { | |
981 | canvas3->cd(i+1); | |
b477f4f2 | 982 | TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone()); |
983 | for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y) | |
984 | { | |
985 | for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z) | |
986 | { | |
987 | hist->SetBinContent(0, y, z, 0); | |
988 | hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); | |
989 | } | |
990 | } | |
991 | TH1* proj = hist->Project3D("zy"); | |
0a173978 | 992 | proj->DrawCopy("COLZ"); |
993 | } | |
994 | } | |
995 | ||
996 | //____________________________________________________________________ | |
447c325d | 997 | void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple) |
0a173978 | 998 | { |
447c325d | 999 | //mcHist->Rebin(2); |
1000 | ||
cfc19dd5 | 1001 | Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); |
1002 | ||
9ca6be09 | 1003 | TString tmpStr; |
cfc19dd5 | 1004 | tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId); |
0a173978 | 1005 | |
447c325d | 1006 | if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0) |
1007 | { | |
1008 | printf("ERROR. Unfolded histogram is empty\n"); | |
1009 | return; | |
1010 | } | |
1011 | ||
1012 | //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists | |
1013 | fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY(); | |
1014 | fCurrentESD->Sumw2(); | |
1015 | fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); | |
1016 | ||
1017 | // normalize unfolded result to 1 | |
1018 | fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral()); | |
1019 | ||
1020 | //fCurrentESD->Scale(mcHist->Integral(2, 200)); | |
1021 | ||
1022 | //new TCanvas; | |
1023 | /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio"); | |
1024 | ratio->Divide(mcHist); | |
1025 | ratio->Draw("HIST"); | |
1026 | ratio->Fit("pol0", "W0", "", 20, 120); | |
1027 | Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0); | |
1028 | delete ratio; | |
1029 | mcHist->Scale(scalingFactor);*/ | |
1030 | ||
1031 | mcHist->Scale(1.0 / mcHist->Integral()); | |
1032 | ||
b477f4f2 | 1033 | // calculate residual |
1034 | ||
1035 | // for that we convolute the response matrix with the gathered result | |
1036 | // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD | |
dd701109 | 1037 | TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected"); |
1038 | tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency); | |
447c325d | 1039 | TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId); |
1040 | TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e"); | |
1041 | if (convolutedProj->Integral() > 0) | |
1042 | { | |
1043 | convolutedProj->Scale(1.0 / convolutedProj->Integral()); | |
1044 | } | |
1045 | else | |
1046 | printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n"); | |
1047 | ||
1048 | //NormalizeToBinWidth(proj2); | |
dd701109 | 1049 | |
447c325d | 1050 | TH1* residual = (TH1*) convolutedProj->Clone("residual"); |
1051 | residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e"); | |
b477f4f2 | 1052 | |
1053 | residual->Add(fCurrentESD, -1); | |
447c325d | 1054 | //residual->Divide(residual, fCurrentESD, 1, 1, "B"); |
1055 | ||
1056 | TH1* residualHist = new TH1F("residualHist", "residualHist", 50, -5, 5); | |
b477f4f2 | 1057 | |
1058 | // TODO fix errors | |
447c325d | 1059 | Float_t chi2 = 0; |
1060 | for (Int_t i=1; i<=residual->GetNbinsX(); ++i) | |
b477f4f2 | 1061 | { |
447c325d | 1062 | if (fCurrentESD->GetBinError(i) > 0) |
1063 | { | |
1064 | Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i); | |
1065 | if (i > 1) | |
1066 | chi2 += value * value; | |
1067 | residual->SetBinContent(i, value); | |
1068 | residualHist->Fill(value); | |
1069 | } | |
1070 | else | |
1071 | { | |
1072 | //printf("Residual bin %d set to 0\n", i); | |
1073 | residual->SetBinContent(i, 0); | |
1074 | } | |
1075 | convolutedProj->SetBinError(i, 0); | |
b477f4f2 | 1076 | residual->SetBinError(i, 0); |
b477f4f2 | 1077 | |
447c325d | 1078 | if (i == 200) |
1079 | fLastChi2Residuals = chi2; | |
1080 | } | |
dd701109 | 1081 | |
447c325d | 1082 | residualHist->Fit("gaus", "N"); |
1083 | delete residualHist; | |
dd701109 | 1084 | |
447c325d | 1085 | printf("Difference (Residuals) is %f for bin 2-200\n", chi2); |
dd701109 | 1086 | |
447c325d | 1087 | TCanvas* canvas1 = 0; |
1088 | if (simple) | |
1089 | { | |
1090 | canvas1 = new TCanvas(tmpStr, tmpStr, 900, 400); | |
1091 | canvas1->Divide(2, 1); | |
1092 | } | |
1093 | else | |
1094 | { | |
1095 | canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200); | |
1096 | canvas1->Divide(2, 3); | |
1097 | } | |
0a173978 | 1098 | |
1099 | canvas1->cd(1); | |
cfc19dd5 | 1100 | TH1* proj = (TH1*) mcHist->Clone("proj"); |
0a173978 | 1101 | NormalizeToBinWidth(proj); |
9ca6be09 | 1102 | |
1103 | if (normalizeESD) | |
1104 | NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]); | |
1105 | ||
dd701109 | 1106 | proj->GetXaxis()->SetRangeUser(0, 200); |
447c325d | 1107 | proj->SetTitle(";true multiplicity;Entries"); |
1108 | proj->SetStats(kFALSE); | |
cfc19dd5 | 1109 | proj->DrawCopy("HIST"); |
0a173978 | 1110 | gPad->SetLogy(); |
1111 | ||
cfc19dd5 | 1112 | //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3); |
1113 | fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2); | |
447c325d | 1114 | fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E"); |
cfc19dd5 | 1115 | |
447c325d | 1116 | TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93); |
1117 | legend->AddEntry(proj, "true distribution"); | |
1118 | legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "unfolded distribution"); | |
1119 | legend->SetFillColor(0); | |
1120 | legend->Draw(); | |
1121 | // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14); | |
1122 | ||
1123 | canvas1->cd(2); | |
1124 | ||
1125 | gPad->SetLogy(); | |
1126 | fCurrentESD->GetXaxis()->SetRangeUser(0, 200); | |
1127 | //fCurrentESD->SetLineColor(2); | |
1128 | fCurrentESD->SetTitle(";measured multiplicity;Entries"); | |
1129 | fCurrentESD->SetStats(kFALSE); | |
1130 | fCurrentESD->DrawCopy("HIST E"); | |
1131 | ||
1132 | convolutedProj->SetLineColor(2); | |
1133 | //proj2->SetMarkerColor(2); | |
1134 | //proj2->SetMarkerStyle(5); | |
1135 | convolutedProj->DrawCopy("HIST SAME"); | |
1136 | ||
1137 | legend = new TLegend(0.3, 0.8, 0.93, 0.93); | |
1138 | legend->AddEntry(fCurrentESD, "measured distribution"); | |
1139 | legend->AddEntry(convolutedProj, "R #otimes unfolded distribution"); | |
1140 | legend->SetFillColor(0); | |
1141 | legend->Draw(); | |
1142 | ||
1143 | TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded")); | |
1144 | diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1); | |
1145 | ||
1146 | /*Float_t chi2 = 0; | |
1147 | Float_t chi = 0; | |
dd701109 | 1148 | fLastChi2MCLimit = 0; |
447c325d | 1149 | Int_t limit = 0; |
dd701109 | 1150 | for (Int_t i=2; i<=200; ++i) |
cfc19dd5 | 1151 | { |
dd701109 | 1152 | if (proj->GetBinContent(i) != 0) |
cfc19dd5 | 1153 | { |
dd701109 | 1154 | Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i); |
cfc19dd5 | 1155 | chi2 += value * value; |
447c325d | 1156 | chi += TMath::Abs(value); |
dd701109 | 1157 | |
447c325d | 1158 | //printf("%d %f\n", i, chi); |
cfc19dd5 | 1159 | |
447c325d | 1160 | if (chi2 < 0.2) |
1161 | fLastChi2MCLimit = i; | |
cfc19dd5 | 1162 | |
447c325d | 1163 | if (chi < 3) |
1164 | limit = i; | |
cfc19dd5 | 1165 | |
447c325d | 1166 | } |
1167 | }*/ | |
9ca6be09 | 1168 | |
cfc19dd5 | 1169 | chi2 = 0; |
447c325d | 1170 | Float_t chi = 0; |
1171 | Int_t limit = 0; | |
1172 | for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i) | |
cfc19dd5 | 1173 | { |
447c325d | 1174 | if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0) |
cfc19dd5 | 1175 | { |
447c325d | 1176 | Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i); |
1177 | if (value > 1e8) | |
1178 | value = 1e8; //prevent arithmetic exception | |
1179 | else if (value < -1e8) | |
1180 | value = -1e8; | |
1181 | if (i > 1) | |
1182 | { | |
1183 | chi2 += value * value; | |
1184 | chi += TMath::Abs(value); | |
1185 | } | |
1186 | diffMCUnfolded->SetBinContent(i, value); | |
cfc19dd5 | 1187 | } |
447c325d | 1188 | else |
1189 | { | |
1190 | //printf("diffMCUnfolded bin %d set to 0\n", i); | |
1191 | diffMCUnfolded->SetBinContent(i, 0); | |
1192 | } | |
1193 | if (chi2 < 1000) | |
1194 | fLastChi2MCLimit = i; | |
1195 | if (chi < 1000) | |
1196 | limit = i; | |
1197 | if (i == 150) | |
1198 | fLastChi2MC = chi2; | |
cfc19dd5 | 1199 | } |
0a173978 | 1200 | |
447c325d | 1201 | printf("limits %d %d\n", fLastChi2MCLimit, limit); |
1202 | fLastChi2MCLimit = limit; | |
1203 | ||
1204 | printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit); | |
0a173978 | 1205 | |
447c325d | 1206 | if (!simple) |
1207 | { | |
1208 | canvas1->cd(3); | |
1209 | ||
1210 | diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e"); | |
1211 | //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20); | |
1212 | diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200); | |
1213 | diffMCUnfolded->DrawCopy("HIST"); | |
0a173978 | 1214 | |
447c325d | 1215 | TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5); |
1216 | for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i) | |
1217 | fluctuation->Fill(diffMCUnfolded->GetBinContent(i)); | |
3602328d | 1218 | |
447c325d | 1219 | new TCanvas; |
1220 | fluctuation->Draw(); | |
1221 | ||
1222 | /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85); | |
1223 | legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected"); | |
1224 | legend->AddEntry(fMultiplicityMC, "MC"); | |
1225 | legend->AddEntry(fMultiplicityESD, "ESD"); | |
1226 | legend->Draw();*/ | |
1227 | ||
1228 | canvas1->cd(4); | |
1229 | //residual->GetYaxis()->SetRangeUser(-0.2, 0.2); | |
1230 | residual->GetXaxis()->SetRangeUser(0, 200); | |
1231 | residual->DrawCopy(); | |
1232 | ||
1233 | canvas1->cd(5); | |
1234 | ||
1235 | TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio"); | |
1236 | ratio->Divide(mcHist); | |
1237 | ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC"); | |
1238 | ratio->GetYaxis()->SetRangeUser(0.5, 1.5); | |
1239 | ratio->GetXaxis()->SetRangeUser(0, 200); | |
1240 | ratio->SetStats(kFALSE); | |
1241 | ratio->Draw("HIST"); | |
1242 | ||
1243 | Double_t ratioChi2 = 0; | |
1244 | fLastChi2MCLimit = 0; | |
1245 | for (Int_t i=2; i<=150; ++i) | |
1246 | { | |
1247 | Float_t value = ratio->GetBinContent(i) - 1; | |
1248 | if (value > 1e8) | |
1249 | value = 1e8; //prevent arithmetic exception | |
1250 | else if (value < -1e8) | |
1251 | value = -1e8; | |
1252 | ||
1253 | ratioChi2 += value * value; | |
1254 | ||
1255 | if (ratioChi2 < 0.1) | |
1256 | fLastChi2MCLimit = i; | |
1257 | } | |
1258 | ||
1259 | printf("Sum over (ratio-1)^2 (2..150) is %f\n", ratioChi2); | |
1260 | // TODO FAKE | |
1261 | fLastChi2MC = ratioChi2; | |
1262 | } | |
3602328d | 1263 | |
b477f4f2 | 1264 | canvas1->SaveAs(Form("%s.gif", canvas1->GetName())); |
0a173978 | 1265 | } |
1266 | ||
1267 | //____________________________________________________________________ | |
dd701109 | 1268 | void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals) |
cfc19dd5 | 1269 | { |
1270 | // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD | |
1271 | // These values are computed during DrawComparison, thus this function picks up the | |
1272 | // last calculation | |
1273 | ||
dd701109 | 1274 | if (mc) |
1275 | *mc = fLastChi2MC; | |
1276 | if (mcLimit) | |
1277 | *mcLimit = fLastChi2MCLimit; | |
1278 | if (residuals) | |
1279 | *residuals = fLastChi2Residuals; | |
cfc19dd5 | 1280 | } |
1281 | ||
1282 | ||
1283 | //____________________________________________________________________ | |
1284 | TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType) | |
1285 | { | |
1286 | // | |
1287 | // returns the corresponding MC spectrum | |
1288 | // | |
1289 | ||
1290 | switch (eventType) | |
1291 | { | |
1292 | case kTrVtx : return fMultiplicityVtx[i]; break; | |
1293 | case kMB : return fMultiplicityMB[i]; break; | |
1294 | case kINEL : return fMultiplicityINEL[i]; break; | |
1295 | } | |
1296 | ||
1297 | return 0; | |
1298 | } | |
1299 | ||
1300 | //____________________________________________________________________ | |
dd701109 | 1301 | void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* inputDist) |
0a173978 | 1302 | { |
1303 | // | |
1304 | // correct spectrum using bayesian method | |
1305 | // | |
1306 | ||
1307 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
0a173978 | 1308 | |
447c325d | 1309 | SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE); |
cfc19dd5 | 1310 | |
447c325d | 1311 | // smooth efficiency |
1312 | //TH1* tmp = (TH1*) fCurrentEfficiency->Clone("eff_clone"); | |
1313 | //for (Int_t i=2; i<fCurrentEfficiency->GetNbinsX(); ++i) | |
1314 | // fCurrentEfficiency->SetBinContent(i, (tmp->GetBinContent(i-1) + tmp->GetBinContent(i) + tmp->GetBinContent(i+1)) / 3); | |
1315 | ||
1316 | // set efficiency to 1 above 150 | |
1317 | // FAKE TEST | |
1318 | //for (Int_t i=150; i<=fCurrentEfficiency->GetNbinsX(); ++i) | |
1319 | // fCurrentEfficiency->SetBinContent(i, 1); | |
0a173978 | 1320 | |
1321 | // normalize correction for given nPart | |
9ca6be09 | 1322 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) |
0a173978 | 1323 | { |
dd701109 | 1324 | // with this it is normalized to 1 |
1325 | Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); | |
1326 | ||
1327 | // with this normalized to the given efficiency | |
1328 | if (fCurrentEfficiency->GetBinContent(i) > 0) | |
1329 | sum /= fCurrentEfficiency->GetBinContent(i); | |
1330 | else | |
1331 | sum = 0; | |
1332 | ||
9ca6be09 | 1333 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) |
0a173978 | 1334 | { |
dd701109 | 1335 | if (sum > 0) |
1336 | { | |
1337 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); | |
1338 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
1339 | } | |
1340 | else | |
1341 | { | |
1342 | fCurrentCorrelation->SetBinContent(i, j, 0); | |
1343 | fCurrentCorrelation->SetBinError(i, j, 0); | |
1344 | } | |
0a173978 | 1345 | } |
1346 | } | |
1347 | ||
447c325d | 1348 | // normalize nTrack |
1349 | /*for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) | |
1350 | { | |
1351 | // with this it is normalized to 1 | |
1352 | Double_t sum = fCurrentCorrelation->Integral(1, fCurrentCorrelation->GetNbinsX(), j, j); | |
1353 | ||
1354 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
1355 | { | |
1356 | if (sum > 0) | |
1357 | { | |
1358 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); | |
1359 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
1360 | } | |
1361 | else | |
1362 | { | |
1363 | fCurrentCorrelation->SetBinContent(i, j, 0); | |
1364 | fCurrentCorrelation->SetBinError(i, j, 0); | |
1365 | } | |
1366 | } | |
1367 | }*/ | |
1368 | ||
dd701109 | 1369 | // smooth input spectrum |
1370 | /* | |
1371 | TH1* esdClone = (TH1*) fCurrentESD->Clone("esdClone"); | |
cfc19dd5 | 1372 | |
dd701109 | 1373 | for (Int_t i=2; i<fCurrentESD->GetNbinsX(); ++i) |
1374 | if (esdClone->GetBinContent(i) < 1e-3) | |
1375 | fCurrentESD->SetBinContent(i, (esdClone->GetBinContent(i-1) + esdClone->GetBinContent(i) + esdClone->GetBinContent(i+1)) / 3); | |
cfc19dd5 | 1376 | |
dd701109 | 1377 | delete esdClone; |
1378 | ||
1379 | // rescale to 1 | |
1380 | fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); | |
1381 | */ | |
6127aca6 | 1382 | |
447c325d | 1383 | /*new TCanvas; |
1384 | fCurrentEfficiency->Draw(); | |
1385 | ||
1386 | new TCanvas; | |
1387 | fCurrentCorrelation->DrawCopy("COLZ"); | |
1388 | ||
1389 | new TCanvas; | |
1390 | ((TH2*) fCurrentCorrelation)->ProjectionX()->DrawCopy(); | |
1391 | ||
1392 | new TCanvas; | |
1393 | ((TH2*) fCurrentCorrelation)->ProjectionY()->DrawCopy();*/ | |
1394 | ||
6127aca6 | 1395 | // pick prior distribution |
dd701109 | 1396 | TH1* hPrior = 0; |
1397 | if (inputDist) | |
1398 | { | |
1399 | printf("Using different starting conditions...\n"); | |
1400 | hPrior = (TH1*)inputDist->Clone("prior"); | |
1401 | } | |
1402 | else | |
1403 | hPrior = (TH1*)fCurrentESD->Clone("prior"); | |
cfc19dd5 | 1404 | Float_t norm = hPrior->Integral(); |
6127aca6 | 1405 | for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) |
6127aca6 | 1406 | hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm); |
1407 | ||
cfc19dd5 | 1408 | // define temp hist |
1409 | TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp"); | |
1410 | hTemp->Reset(); | |
6127aca6 | 1411 | |
cfc19dd5 | 1412 | // just a shortcut |
1413 | TH2F* hResponse = (TH2F*) fCurrentCorrelation; | |
6127aca6 | 1414 | |
cfc19dd5 | 1415 | // histogram to store the inverse response calculated from bayes theorem (from prior and response) IR |
6127aca6 | 1416 | TH2F* hInverseResponseBayes = (TH2F*)hResponse->Clone("pji"); |
6127aca6 | 1417 | hInverseResponseBayes->Reset(); |
b477f4f2 | 1418 | |
dd701109 | 1419 | TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5); |
1420 | ||
447c325d | 1421 | const Int_t kStartBin = 1; |
1422 | ||
6127aca6 | 1423 | // unfold... |
dd701109 | 1424 | for (Int_t i=0; i<nIterations; i++) |
cfc19dd5 | 1425 | { |
1426 | //printf(" iteration %i \n", i); | |
1427 | ||
6127aca6 | 1428 | // calculate IR from Beyes theorem |
1429 | // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k) | |
9ca6be09 | 1430 | |
cfc19dd5 | 1431 | for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) |
1432 | { | |
1433 | Float_t norm = 0; | |
447c325d | 1434 | for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++) |
cfc19dd5 | 1435 | norm += hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t); |
447c325d | 1436 | for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++) |
cfc19dd5 | 1437 | { |
6127aca6 | 1438 | if (norm==0) |
cfc19dd5 | 1439 | hInverseResponseBayes->SetBinContent(t,m,0); |
6127aca6 | 1440 | else |
cfc19dd5 | 1441 | hInverseResponseBayes->SetBinContent(t,m, hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t)/norm); |
6127aca6 | 1442 | } |
1443 | } | |
cfc19dd5 | 1444 | |
447c325d | 1445 | for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++) |
cfc19dd5 | 1446 | { |
1447 | Float_t value = 0; | |
6127aca6 | 1448 | for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) |
cfc19dd5 | 1449 | value += fCurrentESD->GetBinContent(m) * hInverseResponseBayes->GetBinContent(t,m); |
1450 | ||
dd701109 | 1451 | if (fCurrentEfficiency->GetBinContent(t) > 0) |
1452 | hTemp->SetBinContent(t, value / fCurrentEfficiency->GetBinContent(t)); | |
1453 | else | |
1454 | hTemp->SetBinContent(t, 0); | |
cfc19dd5 | 1455 | } |
1456 | ||
1457 | // this is the last guess, fill before (!) smoothing | |
447c325d | 1458 | for (Int_t t=kStartBin; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++) |
cfc19dd5 | 1459 | { |
447c325d | 1460 | //as bin error put the difference to the last iteration |
1461 | //fMultiplicityESDCorrected[correlationID]->SetBinError(t, hTemp->GetBinContent(t) - fMultiplicityESDCorrected[correlationID]->GetBinContent(t)); | |
cfc19dd5 | 1462 | fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t)); |
447c325d | 1463 | fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); |
cfc19dd5 | 1464 | |
1465 | //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t)); | |
6127aca6 | 1466 | } |
9ca6be09 | 1467 | |
447c325d | 1468 | /*new TCanvas; |
1469 | fMultiplicityESDCorrected[correlationID]->DrawCopy(); | |
1470 | gPad->SetLogy();*/ | |
1471 | ||
b477f4f2 | 1472 | // regularization (simple smoothing) |
cfc19dd5 | 1473 | TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed"); |
9ca6be09 | 1474 | |
447c325d | 1475 | for (Int_t t=kStartBin+2; t<hTrueSmoothed->GetNbinsX(); t++) |
cfc19dd5 | 1476 | { |
dd701109 | 1477 | Float_t average = (hTemp->GetBinContent(t-1) |
1478 | + hTemp->GetBinContent(t) | |
1479 | + hTemp->GetBinContent(t+1) | |
1480 | ) / 3.; | |
9ca6be09 | 1481 | |
6127aca6 | 1482 | // weight the average with the regularization parameter |
cfc19dd5 | 1483 | hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average); |
6127aca6 | 1484 | } |
9ca6be09 | 1485 | |
cfc19dd5 | 1486 | // calculate chi2 (change from last iteration) |
1487 | Double_t chi2 = 0; | |
1488 | ||
9ca6be09 | 1489 | // use smoothed true (normalized) as new prior |
cfc19dd5 | 1490 | Float_t norm = 1; |
1491 | //for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) | |
1492 | // norm = norm + hTrueSmoothed->GetBinContent(t); | |
1493 | ||
447c325d | 1494 | for (Int_t t=kStartBin; t<=hTrueSmoothed->GetNbinsX(); t++) |
cfc19dd5 | 1495 | { |
1496 | Float_t newValue = hTrueSmoothed->GetBinContent(t)/norm; | |
1497 | Float_t diff = hPrior->GetBinContent(t) - newValue; | |
1498 | chi2 += (Double_t) diff * diff; | |
1499 | ||
1500 | hPrior->SetBinContent(t, newValue); | |
6127aca6 | 1501 | } |
9ca6be09 | 1502 | |
cfc19dd5 | 1503 | //printf("Chi2 of %d iteration = %.10f\n", i, chi2); |
1504 | ||
dd701109 | 1505 | convergence->Fill(i+1, chi2); |
1506 | ||
cfc19dd5 | 1507 | //if (i % 5 == 0) |
1508 | // DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType); | |
1509 | ||
9ca6be09 | 1510 | delete hTrueSmoothed; |
6127aca6 | 1511 | } // end of iterations |
1512 | ||
dd701109 | 1513 | //new TCanvas; |
1514 | //convergence->DrawCopy(); | |
1515 | //gPad->SetLogy(); | |
1516 | ||
1517 | delete convergence; | |
6127aca6 | 1518 | |
cfc19dd5 | 1519 | return; |
1520 | ||
1521 | // ******** | |
1522 | // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995 | |
1523 | ||
dd701109 | 1524 | /*printf("Calculating covariance matrix. This may take some time...\n"); |
1525 | ||
1526 | // TODO check if this is the right one... | |
1527 | TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX()); | |
cfc19dd5 | 1528 | |
1529 | Int_t xBins = hInverseResponseBayes->GetNbinsX(); | |
1530 | Int_t yBins = hInverseResponseBayes->GetNbinsY(); | |
1531 | ||
1532 | // calculate "unfolding matrix" Mij | |
1533 | Float_t matrixM[251][251]; | |
1534 | for (Int_t i=1; i<=xBins; i++) | |
1535 | { | |
1536 | for (Int_t j=1; j<=yBins; j++) | |
1537 | { | |
1538 | if (fCurrentEfficiency->GetBinContent(i) > 0) | |
1539 | matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i); | |
1540 | else | |
1541 | matrixM[i-1][j-1] = 0; | |
1542 | } | |
1543 | } | |
1544 | ||
1545 | Float_t* vectorn = new Float_t[yBins]; | |
1546 | for (Int_t j=1; j<=yBins; j++) | |
1547 | vectorn[j-1] = fCurrentESD->GetBinContent(j); | |
1548 | ||
1549 | // first part of covariance matrix, depends on input distribution n(E) | |
1550 | Float_t cov1[251][251]; | |
1551 | ||
1552 | Float_t nEvents = fCurrentESD->Integral(); // N | |
1553 | ||
1554 | xBins = 20; | |
1555 | yBins = 20; | |
1556 | ||
1557 | for (Int_t k=0; k<xBins; k++) | |
1558 | { | |
1559 | printf("In Cov1: %d\n", k); | |
1560 | for (Int_t l=0; l<yBins; l++) | |
1561 | { | |
1562 | cov1[k][l] = 0; | |
1563 | ||
1564 | // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N) | |
1565 | for (Int_t j=0; j<yBins; j++) | |
1566 | cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j] | |
1567 | * (1.0 - vectorn[j] / nEvents); | |
1568 | ||
1569 | // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N | |
1570 | for (Int_t i=0; i<yBins; i++) | |
1571 | for (Int_t j=0; j<yBins; j++) | |
1572 | { | |
1573 | if (i == j) | |
1574 | continue; | |
1575 | cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i] | |
1576 | * vectorn[j] / nEvents; | |
1577 | } | |
1578 | } | |
1579 | } | |
1580 | ||
1581 | printf("Cov1 finished\n"); | |
1582 | ||
1583 | TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov"); | |
1584 | cov->Reset(); | |
1585 | ||
1586 | for (Int_t i=1; i<=xBins; i++) | |
1587 | for (Int_t j=1; j<=yBins; j++) | |
1588 | cov->SetBinContent(i, j, cov1[i-1][j-1]); | |
1589 | ||
1590 | new TCanvas; | |
1591 | cov->Draw("COLZ"); | |
1592 | ||
1593 | // second part of covariance matrix, depends on response matrix | |
1594 | Float_t cov2[251][251]; | |
1595 | ||
1596 | // Cov[P(Er|Cu), P(Es|Cu)] term | |
1597 | Float_t covTerm[100][100][100]; | |
1598 | for (Int_t r=0; r<yBins; r++) | |
1599 | for (Int_t u=0; u<xBins; u++) | |
1600 | for (Int_t s=0; s<yBins; s++) | |
1601 | { | |
1602 | if (r == s) | |
1603 | covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) | |
1604 | * (1.0 - hResponse->GetBinContent(u+1, r+1)); | |
1605 | else | |
1606 | covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) | |
1607 | * hResponse->GetBinContent(u+1, s+1); | |
1608 | } | |
1609 | ||
1610 | for (Int_t k=0; k<xBins; k++) | |
1611 | for (Int_t l=0; l<yBins; l++) | |
1612 | { | |
1613 | cov2[k][l] = 0; | |
1614 | printf("In Cov2: %d %d\n", k, l); | |
1615 | for (Int_t i=0; i<yBins; i++) | |
1616 | for (Int_t j=0; j<yBins; j++) | |
1617 | { | |
1618 | //printf("In Cov2: %d %d %d %d\n", k, l, i, j); | |
1619 | // calculate Cov(Mki, Mlj) = sum{ru},{su} ... | |
1620 | Float_t tmpCov = 0; | |
1621 | for (Int_t r=0; r<yBins; r++) | |
1622 | for (Int_t u=0; u<xBins; u++) | |
1623 | for (Int_t s=0; s<yBins; s++) | |
1624 | { | |
1625 | if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0 | |
1626 | || hResponse->GetBinContent(u+1, i+1) == 0) | |
1627 | continue; | |
1628 | ||
1629 | tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u) | |
1630 | * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u) | |
1631 | * covTerm[r][u][s]; | |
1632 | } | |
1633 | ||
1634 | cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov; | |
1635 | } | |
1636 | } | |
1637 | ||
1638 | printf("Cov2 finished\n"); | |
1639 | ||
1640 | for (Int_t i=1; i<=xBins; i++) | |
1641 | for (Int_t j=1; j<=yBins; j++) | |
1642 | cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]); | |
1643 | ||
1644 | new TCanvas; | |
dd701109 | 1645 | cov->Draw("COLZ");*/ |
1646 | } | |
1647 | ||
1648 | //____________________________________________________________________ | |
1649 | Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, | |
1650 | TH1* fCurrentEfficiency, Int_t k, Int_t i, Int_t r, Int_t u) | |
1651 | { | |
1652 | // | |
1653 | // helper function for the covariance matrix of the bayesian method | |
1654 | // | |
1655 | ||
1656 | Float_t result = 0; | |
1657 | ||
1658 | if (k == u && r == i) | |
1659 | result += 1.0 / hResponse->GetBinContent(u+1, r+1); | |
1660 | ||
1661 | if (k == u) | |
1662 | result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1); | |
1663 | ||
1664 | if (r == i) | |
1665 | result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1); | |
1666 | ||
1667 | result *= matrixM[k][i]; | |
1668 | ||
1669 | return result; | |
cfc19dd5 | 1670 | } |
1671 | ||
1672 | //____________________________________________________________________ | |
1673 | void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType) | |
1674 | { | |
1675 | // | |
1676 | // correct spectrum using bayesian method | |
1677 | // | |
1678 | ||
dd701109 | 1679 | Float_t regPar = 0; |
cfc19dd5 | 1680 | |
1681 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
1682 | Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); | |
1683 | ||
447c325d | 1684 | SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE); |
cfc19dd5 | 1685 | |
1686 | // TODO should be taken from correlation map | |
1687 | //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX()); | |
1688 | ||
1689 | // normalize correction for given nPart | |
1690 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
1691 | { | |
1692 | Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); | |
1693 | //Double_t sum = sumHist->GetBinContent(i); | |
1694 | if (sum <= 0) | |
1695 | continue; | |
1696 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) | |
1697 | { | |
1698 | // npart sum to 1 | |
1699 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); | |
1700 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
1701 | } | |
6127aca6 | 1702 | } |
cfc19dd5 | 1703 | |
1704 | new TCanvas; | |
1705 | fCurrentCorrelation->Draw("COLZ"); | |
1706 | ||
1707 | // FAKE | |
1708 | fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff"); | |
1709 | ||
1710 | // pick prior distribution | |
1711 | TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior"); | |
1712 | Float_t norm = 1; //hPrior->Integral(); | |
1713 | for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) | |
1714 | hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm); | |
1715 | ||
1716 | // zero distribution | |
1717 | TH1F* zero = (TH1F*)hPrior->Clone("zero"); | |
1718 | ||
1719 | // define temp hist | |
1720 | TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp"); | |
1721 | hTemp->Reset(); | |
1722 | ||
1723 | // just a shortcut | |
1724 | TH2F* hResponse = (TH2F*) fCurrentCorrelation; | |
1725 | ||
1726 | // unfold... | |
dd701109 | 1727 | Int_t iterations = 25; |
cfc19dd5 | 1728 | for (Int_t i=0; i<iterations; i++) |
1729 | { | |
1730 | //printf(" iteration %i \n", i); | |
1731 | ||
1732 | for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) | |
1733 | { | |
1734 | Float_t value = 0; | |
1735 | for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++) | |
1736 | value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t); | |
1737 | hTemp->SetBinContent(m, value); | |
1738 | //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value); | |
1739 | } | |
1740 | ||
1741 | // regularization (simple smoothing) | |
1742 | TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed"); | |
1743 | ||
1744 | for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++) | |
1745 | { | |
1746 | Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1) | |
1747 | + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t) | |
1748 | + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.; | |
1749 | average *= hTrueSmoothed->GetBinWidth(t); | |
1750 | ||
1751 | // weight the average with the regularization parameter | |
1752 | hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average); | |
1753 | } | |
1754 | ||
1755 | for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) | |
1756 | hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m)); | |
1757 | ||
1758 | // fill guess | |
1759 | for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++) | |
1760 | { | |
1761 | fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t)); | |
1762 | fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO | |
1763 | ||
1764 | //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t)); | |
1765 | } | |
1766 | ||
1767 | ||
1768 | // calculate chi2 (change from last iteration) | |
1769 | Double_t chi2 = 0; | |
1770 | ||
1771 | // use smoothed true (normalized) as new prior | |
1772 | Float_t norm = 1; //hTrueSmoothed->Integral(); | |
1773 | ||
1774 | for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++) | |
1775 | { | |
dd701109 | 1776 | Float_t newValue = hTemp->GetBinContent(t)/norm; |
cfc19dd5 | 1777 | Float_t diff = hPrior->GetBinContent(t) - newValue; |
1778 | chi2 += (Double_t) diff * diff; | |
1779 | ||
1780 | hPrior->SetBinContent(t, newValue); | |
1781 | } | |
1782 | ||
1783 | printf("Chi2 of %d iteration = %.10f\n", i, chi2); | |
1784 | ||
dd701109 | 1785 | //if (i % 5 == 0) |
cfc19dd5 | 1786 | DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY()); |
1787 | ||
1788 | delete hTrueSmoothed; | |
1789 | } // end of iterations | |
1790 | ||
1791 | DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY()); | |
1792 | } | |
1793 | ||
9ca6be09 | 1794 | //____________________________________________________________________ |
1795 | void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace) | |
1796 | { | |
1797 | // | |
1798 | // correct spectrum using a simple Gaussian approach, that is model-dependent | |
1799 | // | |
1800 | ||
1801 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
1802 | Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); | |
1803 | ||
447c325d | 1804 | SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE); |
9ca6be09 | 1805 | |
1806 | NormalizeToBinWidth((TH2*) fCurrentCorrelation); | |
1807 | ||
1808 | TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean")); | |
1809 | correction->SetTitle("GaussianMean"); | |
1810 | ||
1811 | TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth")); | |
1812 | correction->SetTitle("GaussianWidth"); | |
1813 | ||
1814 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
1815 | { | |
1816 | TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1); | |
1817 | proj->Fit("gaus", "0Q"); | |
1818 | correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1)); | |
1819 | correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2)); | |
1820 | ||
1821 | continue; | |
1822 | ||
1823 | // draw for debugging | |
1824 | new TCanvas; | |
1825 | proj->DrawCopy(); | |
1826 | proj->GetFunction("gaus")->DrawCopy("SAME"); | |
1827 | } | |
1828 | ||
1829 | TH1* target = fMultiplicityESDCorrected[correlationID]; | |
1830 | ||
1831 | Int_t nBins = target->GetNbinsX()*10+1; | |
1832 | Float_t* binning = new Float_t[nBins]; | |
1833 | for (Int_t i=1; i<=target->GetNbinsX(); ++i) | |
1834 | for (Int_t j=0; j<10; ++j) | |
1835 | binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j; | |
1836 | ||
1837 | binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX()); | |
1838 | ||
1839 | TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning); | |
1840 | ||
1841 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
1842 | { | |
1843 | Float_t mean = correction->GetBinContent(i); | |
1844 | Float_t width = correctionWidth->GetBinContent(i); | |
1845 | ||
3602328d | 1846 | Int_t fillBegin = fineBinned->FindBin(mean - width * 5); |
1847 | Int_t fillEnd = fineBinned->FindBin(mean + width * 5); | |
1848 | //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd); | |
9ca6be09 | 1849 | |
1850 | for (Int_t j=fillBegin; j <= fillEnd; ++j) | |
1851 | { | |
1852 | fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i)); | |
1853 | } | |
1854 | } | |
1855 | ||
1856 | for (Int_t i=1; i<=target->GetNbinsX(); ++i) | |
1857 | { | |
1858 | Float_t sum = 0; | |
1859 | for (Int_t j=1; j<=10; ++j) | |
1860 | sum += fineBinned->GetBinContent((i-1)*10 + j); | |
1861 | target->SetBinContent(i, sum / 10); | |
1862 | } | |
1863 | ||
1864 | delete[] binning; | |
1865 | ||
cfc19dd5 | 1866 | DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY()); |
0a173978 | 1867 | } |
3602328d | 1868 | |
1869 | //____________________________________________________________________ | |
447c325d | 1870 | TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap) |
3602328d | 1871 | { |
1872 | // runs the distribution given in inputMC through the response matrix identified by | |
1873 | // correlationMap and produces a measured distribution | |
1874 | // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex | |
b477f4f2 | 1875 | // if normalized is set, inputMC is expected to be normalized to the bin width |
3602328d | 1876 | |
1877 | if (!inputMC) | |
1878 | return 0; | |
1879 | ||
1880 | if (correlationMap < 0 || correlationMap >= kCorrHists) | |
1881 | return 0; | |
1882 | ||
1883 | // empty under/overflow bins in x, otherwise Project3D takes them into account | |
1884 | TH3* hist = fCorrelation[correlationMap]; | |
1885 | for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y) | |
1886 | { | |
1887 | for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z) | |
1888 | { | |
1889 | hist->SetBinContent(0, y, z, 0); | |
1890 | hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); | |
1891 | } | |
1892 | } | |
1893 | ||
447c325d | 1894 | TH2* corr = (TH2*) hist->Project3D("zy"); |
1895 | //corr->Rebin2D(2, 1); | |
3602328d | 1896 | corr->Sumw2(); |
1897 | ||
1898 | // normalize correction for given nPart | |
1899 | for (Int_t i=1; i<=corr->GetNbinsX(); ++i) | |
1900 | { | |
1901 | Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY()); | |
1902 | if (sum <= 0) | |
1903 | continue; | |
1904 | ||
1905 | for (Int_t j=1; j<=corr->GetNbinsY(); ++j) | |
1906 | { | |
1907 | // npart sum to 1 | |
1908 | corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum); | |
1909 | corr->SetBinError(i, j, corr->GetBinError(i, j) / sum); | |
1910 | } | |
1911 | } | |
1912 | ||
1913 | TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName()))); | |
1914 | target->Reset(); | |
1915 | ||
1916 | for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas) | |
1917 | { | |
1918 | Float_t measured = 0; | |
1919 | Float_t error = 0; | |
1920 | ||
447c325d | 1921 | for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen) |
3602328d | 1922 | { |
447c325d | 1923 | Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen)); |
b477f4f2 | 1924 | |
447c325d | 1925 | measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas); |
1926 | error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas); | |
3602328d | 1927 | } |
1928 | ||
cfc19dd5 | 1929 | //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error); |
3602328d | 1930 | |
447c325d | 1931 | target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured); |
1932 | target->SetBinError(1 + target->GetNbinsX() / 2, meas, error); | |
3602328d | 1933 | } |
1934 | ||
1935 | return target; | |
1936 | } | |
1937 | ||
1938 | //____________________________________________________________________ | |
1939 | void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id) | |
1940 | { | |
1941 | // uses the given function to fill the input MC histogram and generates from that | |
1942 | // the measured histogram by applying the response matrix | |
1943 | // this can be used to evaluate if the methods work indepedently of the input | |
1944 | // distribution | |
1945 | // WARNING does not respect the vertex distribution, just fills central vertex bin | |
1946 | ||
1947 | if (!inputMC) | |
1948 | return; | |
1949 | ||
1950 | if (id < 0 || id >= kESDHists) | |
1951 | return; | |
1952 | ||
dd701109 | 1953 | TH2F* mc = fMultiplicityVtx[id]; |
3602328d | 1954 | |
1955 | mc->Reset(); | |
1956 | ||
1957 | for (Int_t i=1; i<=mc->GetNbinsY(); ++i) | |
b477f4f2 | 1958 | { |
3602328d | 1959 | mc->SetBinContent(mc->GetNbinsX() / 2, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i)); |
b477f4f2 | 1960 | mc->SetBinError(mc->GetNbinsX() / 2, i, 0); |
1961 | } | |
3602328d | 1962 | |
1963 | //new TCanvas; | |
1964 | //mc->Draw("COLZ"); | |
1965 | ||
1966 | TH1* proj = mc->ProjectionY(); | |
1967 | ||
dd701109 | 1968 | TString tmp(fMultiplicityESD[id]->GetName()); |
1969 | delete fMultiplicityESD[id]; | |
3602328d | 1970 | fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id); |
dd701109 | 1971 | fMultiplicityESD[id]->SetName(tmp); |
3602328d | 1972 | } |