]>
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> |
0a173978 | 36 | |
37 | ClassImp(AliMultiplicityCorrection) | |
38 | ||
9ca6be09 | 39 | const Int_t AliMultiplicityCorrection::fgMaxParams = 90; |
40 | TH1* AliMultiplicityCorrection::fCurrentESD = 0; | |
41 | TH1* AliMultiplicityCorrection::fCurrentCorrelation = 0; | |
0a173978 | 42 | |
43 | //____________________________________________________________________ | |
44 | AliMultiplicityCorrection::AliMultiplicityCorrection() : | |
45 | TNamed() | |
46 | { | |
47 | // | |
48 | // default constructor | |
49 | // | |
50 | ||
51 | for (Int_t i = 0; i < kESDHists; ++i) | |
52 | fMultiplicityESD[i] = 0; | |
53 | ||
54 | for (Int_t i = 0; i < kMCHists; ++i) | |
55 | fMultiplicityMC[i] = 0; | |
56 | ||
57 | for (Int_t i = 0; i < kCorrHists; ++i) | |
58 | { | |
59 | fCorrelation[i] = 0; | |
60 | fMultiplicityESDCorrected[i] = 0; | |
61 | } | |
62 | } | |
63 | ||
64 | //____________________________________________________________________ | |
65 | AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) : | |
66 | TNamed(name, title) | |
67 | { | |
68 | // | |
69 | // named constructor | |
70 | // | |
71 | ||
72 | // do not add this hists to the directory | |
73 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
74 | TH1::AddDirectory(kFALSE); | |
75 | ||
76 | Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10}; | |
77 | 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, | |
78 | 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, | |
79 | 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, | |
80 | 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, | |
81 | 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, | |
82 | 50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5, | |
83 | 100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5, | |
84 | 150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5, | |
85 | 250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5, | |
9ca6be09 | 86 | 500.5 }; |
87 | //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5, | |
88 | //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5, | |
89 | //1000.5 }; | |
0a173978 | 90 | |
9ca6be09 | 91 | #define BINNING fgMaxParams, binLimitsN |
0a173978 | 92 | |
93 | for (Int_t i = 0; i < kESDHists; ++i) | |
94 | fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", 10, binLimitsVtx, BINNING); | |
95 | ||
96 | for (Int_t i = 0; i < kMCHists; ++i) | |
97 | { | |
98 | fMultiplicityMC[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityMC%d", i))); | |
99 | fMultiplicityMC[i]->SetTitle("fMultiplicityMC"); | |
100 | } | |
101 | ||
102 | for (Int_t i = 0; i < kCorrHists; ++i) | |
103 | { | |
104 | fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", 10, binLimitsVtx, BINNING, BINNING); | |
105 | fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", BINNING); | |
106 | } | |
107 | ||
108 | TH1::AddDirectory(oldStatus); | |
109 | } | |
110 | ||
111 | //____________________________________________________________________ | |
112 | AliMultiplicityCorrection::~AliMultiplicityCorrection() | |
113 | { | |
114 | // | |
115 | // Destructor | |
116 | // | |
117 | } | |
118 | ||
119 | //____________________________________________________________________ | |
120 | Long64_t AliMultiplicityCorrection::Merge(TCollection* list) | |
121 | { | |
122 | // Merge a list of AliMultiplicityCorrection objects with this (needed for | |
123 | // PROOF). | |
124 | // Returns the number of merged objects (including this). | |
125 | ||
126 | if (!list) | |
127 | return 0; | |
128 | ||
129 | if (list->IsEmpty()) | |
130 | return 1; | |
131 | ||
132 | TIterator* iter = list->MakeIterator(); | |
133 | TObject* obj; | |
134 | ||
135 | // collections of all histograms | |
136 | TList collections[kESDHists+kMCHists+kCorrHists*2]; | |
137 | ||
138 | Int_t count = 0; | |
139 | while ((obj = iter->Next())) { | |
140 | ||
141 | AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj); | |
142 | if (entry == 0) | |
143 | continue; | |
144 | ||
145 | for (Int_t i = 0; i < kESDHists; ++i) | |
146 | collections[i].Add(entry->fMultiplicityESD[i]); | |
147 | ||
148 | for (Int_t i = 0; i < kMCHists; ++i) | |
149 | collections[kESDHists+i].Add(entry->fMultiplicityMC[i]); | |
150 | ||
151 | for (Int_t i = 0; i < kCorrHists; ++i) | |
152 | collections[kESDHists+kMCHists+i].Add(entry->fCorrelation[i]); | |
153 | ||
154 | for (Int_t i = 0; i < kCorrHists; ++i) | |
155 | collections[kESDHists+kMCHists+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]); | |
156 | ||
157 | count++; | |
158 | } | |
159 | ||
160 | for (Int_t i = 0; i < kESDHists; ++i) | |
161 | fMultiplicityESD[i]->Merge(&collections[i]); | |
162 | ||
163 | for (Int_t i = 0; i < kMCHists; ++i) | |
164 | fMultiplicityMC[i]->Merge(&collections[kESDHists+i]); | |
165 | ||
166 | for (Int_t i = 0; i < kCorrHists; ++i) | |
167 | fCorrelation[i]->Merge(&collections[kESDHists+kMCHists+i]); | |
168 | ||
169 | for (Int_t i = 0; i < kCorrHists; ++i) | |
170 | fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists+kCorrHists+i]); | |
171 | ||
172 | delete iter; | |
173 | ||
174 | return count+1; | |
175 | } | |
176 | ||
177 | //____________________________________________________________________ | |
178 | Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir) | |
179 | { | |
180 | // | |
181 | // loads the histograms from a file | |
182 | // if dir is empty a directory with the name of this object is taken (like in SaveHistogram) | |
183 | // | |
184 | ||
185 | if (!dir) | |
186 | dir = GetName(); | |
187 | ||
188 | if (!gDirectory->cd(dir)) | |
189 | return kFALSE; | |
190 | ||
191 | // TODO memory leak. old histograms needs to be deleted. | |
192 | ||
193 | Bool_t success = kTRUE; | |
194 | ||
195 | for (Int_t i = 0; i < kESDHists; ++i) | |
196 | { | |
197 | fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName())); | |
198 | if (!fMultiplicityESD[i]) | |
199 | success = kFALSE; | |
200 | } | |
201 | ||
202 | for (Int_t i = 0; i < kMCHists; ++i) | |
203 | { | |
204 | fMultiplicityMC[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMC[i]->GetName())); | |
205 | if (!fMultiplicityMC[i]) | |
206 | success = kFALSE; | |
207 | } | |
208 | ||
209 | for (Int_t i = 0; i < kCorrHists; ++i) | |
210 | { | |
211 | fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName())); | |
212 | if (!fCorrelation[i]) | |
213 | success = kFALSE; | |
214 | fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName())); | |
215 | if (!fMultiplicityESDCorrected[i]) | |
216 | success = kFALSE; | |
217 | } | |
218 | ||
219 | gDirectory->cd(".."); | |
220 | ||
221 | return success; | |
222 | } | |
223 | ||
224 | //____________________________________________________________________ | |
225 | void AliMultiplicityCorrection::SaveHistograms() | |
226 | { | |
227 | // | |
228 | // saves the histograms | |
229 | // | |
230 | ||
231 | gDirectory->mkdir(GetName()); | |
232 | gDirectory->cd(GetName()); | |
233 | ||
234 | for (Int_t i = 0; i < kESDHists; ++i) | |
235 | if (fMultiplicityESD[i]) | |
236 | fMultiplicityESD[i]->Write(); | |
237 | ||
238 | for (Int_t i = 0; i < kMCHists; ++i) | |
239 | if (fMultiplicityMC[i]) | |
240 | fMultiplicityMC[i]->Write(); | |
241 | ||
242 | for (Int_t i = 0; i < kCorrHists; ++i) | |
243 | { | |
244 | if (fCorrelation[i]) | |
245 | fCorrelation[i]->Write(); | |
246 | if (fMultiplicityESDCorrected[i]) | |
247 | fMultiplicityESDCorrected[i]->Write(); | |
248 | } | |
249 | ||
250 | gDirectory->cd(".."); | |
251 | } | |
252 | ||
253 | //____________________________________________________________________ | |
254 | void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll) | |
255 | { | |
256 | // | |
257 | // Fills an event from MC | |
258 | // | |
259 | ||
260 | fMultiplicityMC[0]->Fill(vtx, generated05); | |
261 | fMultiplicityMC[1]->Fill(vtx, generated10); | |
262 | fMultiplicityMC[2]->Fill(vtx, generated15); | |
263 | fMultiplicityMC[3]->Fill(vtx, generated20); | |
264 | fMultiplicityMC[4]->Fill(vtx, generatedAll); | |
265 | } | |
266 | ||
267 | //____________________________________________________________________ | |
268 | void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20) | |
269 | { | |
270 | // | |
271 | // Fills an event from ESD | |
272 | // | |
273 | ||
274 | fMultiplicityESD[0]->Fill(vtx, measured05); | |
275 | fMultiplicityESD[1]->Fill(vtx, measured10); | |
276 | fMultiplicityESD[2]->Fill(vtx, measured15); | |
277 | fMultiplicityESD[3]->Fill(vtx, measured20); | |
278 | } | |
279 | ||
280 | //____________________________________________________________________ | |
281 | 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) | |
282 | { | |
283 | // | |
284 | // Fills an event into the correlation map with the information from MC and ESD | |
285 | // | |
286 | ||
287 | fCorrelation[0]->Fill(vtx, generated05, measured05); | |
288 | fCorrelation[1]->Fill(vtx, generated10, measured10); | |
289 | fCorrelation[2]->Fill(vtx, generated15, measured15); | |
290 | fCorrelation[3]->Fill(vtx, generated20, measured20); | |
291 | ||
292 | fCorrelation[4]->Fill(vtx, generatedAll, measured05); | |
293 | fCorrelation[5]->Fill(vtx, generatedAll, measured10); | |
294 | fCorrelation[6]->Fill(vtx, generatedAll, measured15); | |
295 | fCorrelation[7]->Fill(vtx, generatedAll, measured20); | |
296 | } | |
297 | ||
298 | //____________________________________________________________________ | |
9ca6be09 | 299 | Double_t AliMultiplicityCorrection::RegularizationPol0(Double_t *params) |
0a173978 | 300 | { |
301 | // homogenity term for minuit fitting | |
302 | // pure function of the parameters | |
303 | // prefers constant function (pol0) | |
304 | ||
305 | Double_t chi2 = 0; | |
306 | ||
307 | for (Int_t i=1; i<fgMaxParams; ++i) | |
308 | { | |
309 | if (params[i] == 0) | |
310 | continue; | |
311 | ||
9ca6be09 | 312 | Double_t right = params[i] / fCurrentESD->GetBinWidth(i+1); |
313 | Double_t left = params[i-1] / fCurrentESD->GetBinWidth(i); | |
0a173978 | 314 | |
315 | Double_t diff = (right - left) / right; | |
316 | chi2 += diff * diff; | |
317 | } | |
318 | ||
319 | return chi2; | |
320 | } | |
321 | ||
322 | //____________________________________________________________________ | |
9ca6be09 | 323 | Double_t AliMultiplicityCorrection::RegularizationPol1(Double_t *params) |
0a173978 | 324 | { |
325 | // homogenity term for minuit fitting | |
326 | // pure function of the parameters | |
327 | // prefers linear function (pol1) | |
328 | ||
329 | Double_t chi2 = 0; | |
330 | ||
331 | for (Int_t i=2; i<fgMaxParams; ++i) | |
332 | { | |
333 | if (params[i] == 0 || params[i-1] == 0) | |
334 | continue; | |
335 | ||
9ca6be09 | 336 | Double_t right = params[i] / fCurrentESD->GetBinWidth(i+1); |
337 | Double_t middle = params[i-1] / fCurrentESD->GetBinWidth(i); | |
338 | Double_t left = params[i-2] / fCurrentESD->GetBinWidth(i-1); | |
0a173978 | 339 | |
9ca6be09 | 340 | Double_t der1 = (right - middle) / fCurrentESD->GetBinWidth(i); |
341 | Double_t der2 = (middle - left) / fCurrentESD->GetBinWidth(i-1); | |
0a173978 | 342 | |
343 | Double_t diff = der1 - der2; | |
344 | ||
345 | // TODO give additonal weight to big bins | |
9ca6be09 | 346 | chi2 += diff * diff * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1); |
0a173978 | 347 | |
348 | //printf("%d --> %f\n", i, diff); | |
349 | } | |
350 | ||
351 | return chi2 / 1e5 / 2; | |
352 | } | |
353 | ||
9ca6be09 | 354 | //____________________________________________________________________ |
355 | Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(Double_t *params) | |
356 | { | |
357 | // homogenity term for minuit fitting | |
358 | // pure function of the parameters | |
359 | // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments, | |
360 | // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp. | |
361 | ||
362 | Double_t chi2 = 0; | |
363 | ||
364 | for (Int_t i=2; i<fgMaxParams; ++i) | |
365 | { | |
366 | if (params[i] == 0 || params[i-1] == 0) | |
367 | continue; | |
368 | ||
369 | Double_t right = params[i] / fCurrentESD->GetBinWidth(i+1); | |
370 | Double_t middle = params[i-1] / fCurrentESD->GetBinWidth(i); | |
371 | Double_t left = params[i-2] / fCurrentESD->GetBinWidth(i-1); | |
372 | ||
373 | Double_t der1 = (right - middle) / fCurrentESD->GetBinWidth(i); | |
374 | Double_t der2 = (middle - left) / fCurrentESD->GetBinWidth(i-1); | |
375 | ||
376 | Double_t secDer = (der1 - der2) / fCurrentESD->GetBinWidth(i); | |
377 | ||
378 | // square and weight with the bin width | |
379 | chi2 += secDer * secDer * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1); | |
380 | ||
381 | //printf("%d --> %f\n", i, secDer); | |
382 | } | |
383 | ||
384 | return chi2 / 1e5; | |
385 | } | |
386 | ||
0a173978 | 387 | //____________________________________________________________________ |
388 | void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t) | |
389 | { | |
390 | // | |
391 | // fit function for minuit | |
392 | // | |
393 | ||
394 | // TODO take errors into account | |
395 | ||
396 | static Int_t callCount = 0; | |
397 | ||
398 | Double_t chi2FromFit = 0; | |
399 | ||
400 | // loop over multiplicity (x axis of fMultiplicityESD) | |
9ca6be09 | 401 | for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i) |
0a173978 | 402 | { |
9ca6be09 | 403 | if (i > fCurrentCorrelation->GetNbinsY()) |
0a173978 | 404 | break; |
405 | ||
406 | Double_t sum = 0; | |
407 | //Double_t error = 0; | |
408 | // loop over generated particles (x axis of fCorrelation) that resulted in reconstruction of i tracks | |
9ca6be09 | 409 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsX(); ++j) |
0a173978 | 410 | { |
411 | if (j > fgMaxParams) | |
412 | break; | |
413 | ||
9ca6be09 | 414 | sum += fCurrentCorrelation->GetBinContent(j, i) * params[j-1]; |
0a173978 | 415 | |
416 | //if (params[j-1] > 0) | |
9ca6be09 | 417 | // error += fCurrentCorrelation->GetBinError(j, i) * fCurrentCorrelation->GetBinError(j, i) * params[j-1]; |
0a173978 | 418 | //printf("%f ", sum); |
419 | } | |
420 | ||
9ca6be09 | 421 | Double_t diff = fCurrentESD->GetBinContent(i) - sum; |
422 | if (fCurrentESD->GetBinContent(i) > 0) | |
423 | diff /= fCurrentESD->GetBinContent(i); | |
0a173978 | 424 | else |
9ca6be09 | 425 | diff /= fCurrentESD->Integral(); |
426 | ||
427 | // weight with bin width | |
428 | //diff *= fCurrentESD->GetBinWidth(i); | |
429 | ||
430 | //error = TMath::Sqrt(error) + fCurrentESD->GetBinError(i); | |
0a173978 | 431 | //if (error <= 0) |
432 | // error = 1; | |
433 | ||
434 | //Double_t tmp = diff / error; | |
435 | //chi2 += tmp * tmp; | |
436 | chi2FromFit += diff * diff; | |
437 | ||
9ca6be09 | 438 | //printf("\nExpected sum = %f; Diff for bin %d is %f\n**********************************\n", fCurrentESD->GetBinContent(i), i, diff); |
0a173978 | 439 | //printf("Diff for bin %d is %f\n", i, diff); |
440 | } | |
441 | ||
9ca6be09 | 442 | Double_t penaltyVal = RegularizationTotalCurvature(params); |
0a173978 | 443 | |
444 | chi2 = chi2FromFit * chi2FromFit + penaltyVal * penaltyVal; | |
445 | ||
446 | if ((callCount++ % 100) == 0) | |
447 | printf("%f %f --> %f\n", chi2FromFit, penaltyVal, chi2); | |
448 | } | |
449 | ||
450 | //____________________________________________________________________ | |
9ca6be09 | 451 | void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace) |
0a173978 | 452 | { |
453 | // | |
9ca6be09 | 454 | // fills fCurrentESD, fCurrentCorrelation |
455 | // resets fMultiplicityESDCorrected | |
0a173978 | 456 | // |
457 | ||
458 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
9ca6be09 | 459 | fMultiplicityESDCorrected[correlationID]->Reset(); |
0a173978 | 460 | |
9ca6be09 | 461 | fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY(); |
462 | fCurrentESD->Sumw2(); | |
0a173978 | 463 | |
464 | // empty under/overflow bins in x, otherwise Project3D takes them into account | |
465 | TH3* hist = fCorrelation[correlationID]; | |
466 | for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y) | |
467 | { | |
468 | for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z) | |
469 | { | |
470 | hist->SetBinContent(0, y, z, 0); | |
471 | hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); | |
472 | } | |
473 | } | |
474 | ||
9ca6be09 | 475 | fCurrentCorrelation = hist->Project3D("zy"); |
476 | fCurrentCorrelation->Sumw2(); | |
477 | } | |
478 | ||
479 | //____________________________________________________________________ | |
480 | void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, Bool_t check) | |
481 | { | |
482 | // | |
483 | // correct spectrum using minuit chi2 method | |
484 | // | |
485 | // if check is kTRUE the input MC solution (by definition the right solution) is used | |
486 | // no fit is made and just the chi2 is calculated | |
487 | // | |
488 | ||
489 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
490 | Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); | |
491 | ||
492 | SetupCurrentHists(inputRange, fullPhaseSpace); | |
0a173978 | 493 | |
494 | // normalize correction for given nPart | |
9ca6be09 | 495 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) |
0a173978 | 496 | { |
9ca6be09 | 497 | Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); |
0a173978 | 498 | if (sum <= 0) |
499 | continue; | |
9ca6be09 | 500 | Float_t maxValue = 0; |
501 | Int_t maxBin = -1; | |
502 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) | |
0a173978 | 503 | { |
9ca6be09 | 504 | // find most probably value |
505 | if (maxValue < fCurrentCorrelation->GetBinContent(i, j)) | |
506 | { | |
507 | maxValue = fCurrentCorrelation->GetBinContent(i, j); | |
508 | maxBin = j; | |
509 | } | |
510 | ||
0a173978 | 511 | // npart sum to 1 |
9ca6be09 | 512 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); |
513 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
0a173978 | 514 | } |
9ca6be09 | 515 | |
516 | printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin)); | |
0a173978 | 517 | } |
518 | ||
519 | // small hack to get around charge conservation for full phase space ;-) | |
520 | /*if (fullPhaseSpace) | |
521 | { | |
522 | for (Int_t i=2; i<=50; i+=2) | |
9ca6be09 | 523 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) |
524 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i-1, j)); | |
0a173978 | 525 | }*/ |
526 | ||
527 | TCanvas* canvas1 = new TCanvas("ApplyMinuitFit", "ApplyMinuitFit", 800, 400); | |
528 | canvas1->Divide(2, 1); | |
529 | canvas1->cd(1); | |
9ca6be09 | 530 | fCurrentESD->DrawCopy(); |
0a173978 | 531 | canvas1->cd(2); |
9ca6be09 | 532 | fCurrentCorrelation->DrawCopy("COLZ"); |
0a173978 | 533 | |
534 | // Initialize TMinuit via generic fitter interface | |
535 | TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams); | |
536 | ||
537 | minuit->SetFCN(MinuitFitFunction); | |
538 | ||
539 | Double_t results[fgMaxParams]; | |
540 | ||
9ca6be09 | 541 | TH1* proj = fMultiplicityMC[mcTarget]->ProjectionY(); |
0a173978 | 542 | for (Int_t i=0; i<fgMaxParams; ++i) |
543 | { | |
544 | //results[i] = 1; | |
9ca6be09 | 545 | results[i] = fCurrentESD->GetBinContent(i+1); |
546 | if (check) | |
547 | results[i] = proj->GetBinContent(i+1); | |
548 | minuit->SetParameter(i, Form("param%d", i), results[i], ((results[i] > 1) ? (results[i] / 10) : 0), 0, fCurrentESD->GetMaximum() * 100); | |
0a173978 | 549 | } |
550 | ||
551 | Int_t dummy; | |
552 | Double_t chi2; | |
553 | MinuitFitFunction(dummy, 0, chi2, results, 0); | |
554 | printf("Chi2 of right solution is = %f\n", chi2); | |
555 | ||
9ca6be09 | 556 | if (check) |
557 | return; | |
0a173978 | 558 | |
559 | Double_t arglist[100]; | |
560 | arglist[0] = 100000; | |
561 | //minuit->ExecuteCommand("SET PRINT", arglist, 1); | |
562 | minuit->ExecuteCommand("MIGRAD", arglist, 1); | |
563 | //minuit->ExecuteCommand("MINOS", arglist, 0); | |
564 | ||
565 | for (Int_t i=0; i<fgMaxParams; ++i) | |
566 | { | |
567 | results[i] = minuit->GetParameter(i); | |
568 | fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, results[i]); | |
569 | fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0); | |
570 | } | |
571 | ||
9ca6be09 | 572 | printf("Penalty is %f\n", RegularizationTotalCurvature(results)); |
0a173978 | 573 | |
9ca6be09 | 574 | DrawComparison("MinuitChi2", mcTarget, correlationID); |
0a173978 | 575 | } |
576 | ||
577 | //____________________________________________________________________ | |
578 | void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist) | |
579 | { | |
580 | // | |
581 | // normalizes a 1-d histogram to its bin width | |
582 | // | |
583 | ||
584 | for (Int_t i=1; i<=hist->GetNbinsX(); ++i) | |
585 | { | |
586 | hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i)); | |
587 | hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i)); | |
588 | } | |
589 | } | |
590 | ||
591 | //____________________________________________________________________ | |
592 | void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist) | |
593 | { | |
594 | // | |
595 | // normalizes a 2-d histogram to its bin width (x width * y width) | |
596 | // | |
597 | ||
598 | for (Int_t i=1; i<=hist->GetNbinsX(); ++i) | |
599 | for (Int_t j=1; j<=hist->GetNbinsY(); ++j) | |
600 | { | |
601 | Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j); | |
602 | hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor); | |
603 | hist->SetBinError(i, j, hist->GetBinError(i, j) / factor); | |
604 | } | |
605 | } | |
606 | ||
607 | //____________________________________________________________________ | |
608 | void AliMultiplicityCorrection::ApplyMinuitFitAll() | |
609 | { | |
610 | // | |
611 | // fit all eta ranges | |
612 | // | |
613 | ||
614 | for (Int_t i=0; i<kESDHists; ++i) | |
615 | { | |
616 | ApplyMinuitFit(i, kFALSE); | |
617 | ApplyMinuitFit(i, kTRUE); | |
618 | } | |
619 | } | |
620 | ||
621 | //____________________________________________________________________ | |
622 | void AliMultiplicityCorrection::DrawHistograms() | |
623 | { | |
624 | printf("ESD:\n"); | |
625 | ||
626 | TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600); | |
627 | canvas1->Divide(3, 2); | |
628 | for (Int_t i = 0; i < kESDHists; ++i) | |
629 | { | |
630 | canvas1->cd(i+1); | |
631 | fMultiplicityESD[i]->DrawCopy("COLZ"); | |
632 | printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean()); | |
633 | } | |
634 | ||
635 | printf("MC:\n"); | |
636 | ||
637 | TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600); | |
638 | canvas2->Divide(3, 2); | |
639 | for (Int_t i = 0; i < kMCHists; ++i) | |
640 | { | |
641 | canvas2->cd(i+1); | |
642 | fMultiplicityMC[i]->DrawCopy("COLZ"); | |
643 | printf("%d --> %f\n", i, fMultiplicityMC[i]->ProjectionY()->GetMean()); | |
644 | } | |
645 | ||
646 | TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900); | |
647 | canvas3->Divide(3, 3); | |
648 | for (Int_t i = 0; i < kCorrHists; ++i) | |
649 | { | |
650 | canvas3->cd(i+1); | |
651 | TH1* proj = fCorrelation[i]->Project3D("zy"); | |
652 | proj->DrawCopy("COLZ"); | |
653 | } | |
654 | } | |
655 | ||
656 | //____________________________________________________________________ | |
9ca6be09 | 657 | void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t mcID, Int_t esdCorrId, Bool_t normalizeESD) |
0a173978 | 658 | { |
9ca6be09 | 659 | TString tmpStr; |
660 | tmpStr.Form("%s_DrawComparison_%d_%d", name, mcID, esdCorrId); | |
0a173978 | 661 | |
9ca6be09 | 662 | TCanvas* canvas1 = new TCanvas(tmpStr, tmpStr, 900, 300); |
0a173978 | 663 | canvas1->Divide(3, 1); |
664 | ||
665 | canvas1->cd(1); | |
666 | TH1* proj = fMultiplicityMC[mcID]->ProjectionY(); | |
667 | NormalizeToBinWidth(proj); | |
9ca6be09 | 668 | |
669 | if (normalizeESD) | |
670 | NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]); | |
671 | ||
0a173978 | 672 | proj->GetXaxis()->SetRangeUser(0, 200); |
673 | proj->DrawCopy("HIST"); | |
674 | gPad->SetLogy(); | |
675 | ||
9ca6be09 | 676 | NormalizeToBinWidth(fCurrentESD); |
677 | fCurrentESD->SetLineColor(2); | |
678 | fCurrentESD->DrawCopy("HISTSAME"); | |
679 | ||
0a173978 | 680 | fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3); |
681 | fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAMEP"); | |
682 | ||
683 | canvas1->cd(2); | |
9ca6be09 | 684 | fMultiplicityESDCorrected[esdCorrId]->GetXaxis()->SetRangeUser(0, 200); |
685 | fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST"); | |
0a173978 | 686 | |
687 | canvas1->cd(3); | |
688 | TH1* clone = dynamic_cast<TH1*> (proj->Clone("clone")); | |
689 | clone->Divide(fMultiplicityESDCorrected[esdCorrId]); | |
690 | clone->SetTitle("Ratio;Npart;MC/ESD"); | |
691 | clone->GetYaxis()->SetRangeUser(0.8, 1.2); | |
692 | clone->GetXaxis()->SetRangeUser(0, 200); | |
693 | clone->DrawCopy("HIST"); | |
694 | ||
695 | /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85); | |
696 | legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected"); | |
697 | legend->AddEntry(fMultiplicityMC, "MC"); | |
698 | legend->AddEntry(fMultiplicityESD, "ESD"); | |
699 | legend->Draw();*/ | |
700 | } | |
701 | ||
702 | //____________________________________________________________________ | |
703 | void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace) | |
704 | { | |
705 | // | |
706 | // correct spectrum using bayesian method | |
707 | // | |
708 | ||
709 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
710 | Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); | |
711 | ||
9ca6be09 | 712 | SetupCurrentHists(inputRange, fullPhaseSpace); |
0a173978 | 713 | |
714 | // normalize correction for given nPart | |
9ca6be09 | 715 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) |
0a173978 | 716 | { |
9ca6be09 | 717 | Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); |
0a173978 | 718 | if (sum <= 0) |
719 | continue; | |
9ca6be09 | 720 | for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) |
0a173978 | 721 | { |
722 | // npart sum to 1 | |
9ca6be09 | 723 | fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); |
724 | fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); | |
0a173978 | 725 | } |
726 | } | |
727 | ||
6127aca6 | 728 | Float_t regPar = 0.01; |
729 | ||
730 | Float_t norm = 0; | |
731 | // pick prior distribution | |
9ca6be09 | 732 | TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior"); |
6127aca6 | 733 | for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) |
734 | norm = norm + hPrior->GetBinContent(t); | |
735 | for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) { | |
736 | hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm); | |
737 | ||
738 | printf(" bin %d content %f \n", t, hPrior->GetBinContent(t)); | |
739 | ||
740 | } | |
741 | ||
742 | // define hist to store guess of true | |
9ca6be09 | 743 | TH1F* hTrue = (TH1F*)fCurrentESD->Clone("prior"); |
6127aca6 | 744 | // hTrue->Reset(); |
745 | ||
746 | // clone response R | |
9ca6be09 | 747 | TH2F* hResponse = (TH2F*)fCurrentCorrelation->Clone("pij"); |
6127aca6 | 748 | |
9ca6be09 | 749 | // histogram to store the inverse response calculated from bayes theorem (from prior and response) |
6127aca6 | 750 | // IR |
9ca6be09 | 751 | //TAxis* xAxis = hResponse->GetYaxis(); |
752 | //TAxis* yAxis = hResponse->GetXaxis(); | |
6127aca6 | 753 | |
754 | TH2F* hInverseResponseBayes = (TH2F*)hResponse->Clone("pji"); | |
755 | //new TH2F("pji","pji", | |
756 | // xAxis->GetNbins(),xAxis->GetXbins()->GetArray(), | |
757 | // yAxis->GetNbins(),yAxis->GetXbins()->GetArray()); | |
758 | hInverseResponseBayes->Reset(); | |
759 | ||
760 | // unfold... | |
9ca6be09 | 761 | Int_t iterations = 50; |
6127aca6 | 762 | for (Int_t i=0; i<iterations; i++) { |
763 | printf(" iteration %i \n", i); | |
764 | ||
765 | // calculate IR from Beyes theorem | |
766 | // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k) | |
9ca6be09 | 767 | for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++) { |
6127aca6 | 768 | for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) { |
769 | ||
770 | norm = 0; | |
771 | for(Int_t t2 = 1; t2<=hResponse->GetNbinsX(); t2++) | |
9ca6be09 | 772 | norm += (hResponse->GetBinContent(t2,m) * hPrior->GetBinContent(t2)); |
773 | ||
6127aca6 | 774 | if (norm==0) |
775 | hInverseResponseBayes->SetBinContent(t,m,0); | |
776 | else | |
777 | hInverseResponseBayes->SetBinContent(t,m, hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t)/norm); | |
778 | } | |
779 | } | |
9ca6be09 | 780 | // calculate true assuming IR is the correct inverse response |
781 | for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++) { | |
6127aca6 | 782 | hTrue ->SetBinContent(t, 0); |
783 | for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) | |
9ca6be09 | 784 | hTrue->SetBinContent(t, hTrue->GetBinContent(t) + fCurrentESD->GetBinContent(m)*hInverseResponseBayes->GetBinContent(t,m)); |
6127aca6 | 785 | } |
9ca6be09 | 786 | |
6127aca6 | 787 | // regularization (simple smoothing) NB : not good bin width!!! |
788 | TH1F* hTrueSmoothed = (TH1F*)hTrue->Clone("truesmoothed"); | |
9ca6be09 | 789 | |
6127aca6 | 790 | for (Int_t t=2; t<hTrue->GetNbinsX()-1; t++) { |
9ca6be09 | 791 | Float_t average = (hTrue->GetBinContent(t-1) / hTrue->GetBinWidth(t-1) |
792 | + hTrue->GetBinContent(t) / hTrue->GetBinWidth(t) | |
793 | + hTrue->GetBinContent(t+1) / hTrue->GetBinWidth(t+1)) / 3.; | |
794 | average *= hTrueSmoothed->GetBinWidth(t); | |
795 | ||
6127aca6 | 796 | // weight the average with the regularization parameter |
797 | hTrueSmoothed->SetBinContent(t, (1-regPar)*hTrue->GetBinContent(t) + (regPar*average)); | |
798 | } | |
9ca6be09 | 799 | |
800 | // use smoothed true (normalized) as new prior | |
6127aca6 | 801 | norm = 0; |
802 | for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) | |
803 | norm = norm + hTrueSmoothed->GetBinContent(t); | |
804 | for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) { | |
805 | hPrior->SetBinContent(t, hTrueSmoothed->GetBinContent(t)/norm); | |
806 | hTrue->SetBinContent(t, hTrueSmoothed->GetBinContent(t)); | |
807 | } | |
9ca6be09 | 808 | |
809 | delete hTrueSmoothed; | |
6127aca6 | 810 | } // end of iterations |
811 | ||
9ca6be09 | 812 | for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++) { |
813 | fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTrue->GetBinContent(t)); | |
814 | fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05*hTrue->GetBinContent(t)); | |
6127aca6 | 815 | |
9ca6be09 | 816 | printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t)); |
6127aca6 | 817 | } |
818 | ||
9ca6be09 | 819 | DrawComparison("Bayesian", mcTarget, correlationID); |
820 | } | |
821 | ||
822 | //____________________________________________________________________ | |
823 | void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace) | |
824 | { | |
825 | // | |
826 | // correct spectrum using a simple Gaussian approach, that is model-dependent | |
827 | // | |
828 | ||
829 | Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); | |
830 | Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); | |
831 | ||
832 | SetupCurrentHists(inputRange, fullPhaseSpace); | |
833 | ||
834 | NormalizeToBinWidth((TH2*) fCurrentCorrelation); | |
835 | ||
836 | TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean")); | |
837 | correction->SetTitle("GaussianMean"); | |
838 | ||
839 | TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth")); | |
840 | correction->SetTitle("GaussianWidth"); | |
841 | ||
842 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
843 | { | |
844 | TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1); | |
845 | proj->Fit("gaus", "0Q"); | |
846 | correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1)); | |
847 | correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2)); | |
848 | ||
849 | continue; | |
850 | ||
851 | // draw for debugging | |
852 | new TCanvas; | |
853 | proj->DrawCopy(); | |
854 | proj->GetFunction("gaus")->DrawCopy("SAME"); | |
855 | } | |
856 | ||
857 | TH1* target = fMultiplicityESDCorrected[correlationID]; | |
858 | ||
859 | Int_t nBins = target->GetNbinsX()*10+1; | |
860 | Float_t* binning = new Float_t[nBins]; | |
861 | for (Int_t i=1; i<=target->GetNbinsX(); ++i) | |
862 | for (Int_t j=0; j<10; ++j) | |
863 | binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j; | |
864 | ||
865 | binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX()); | |
866 | ||
867 | TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning); | |
868 | ||
869 | for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) | |
870 | { | |
871 | Float_t mean = correction->GetBinContent(i); | |
872 | Float_t width = correctionWidth->GetBinContent(i); | |
873 | ||
874 | Int_t fillBegin = fineBinned->FindBin(mean - width * 3); | |
875 | Int_t fillEnd = fineBinned->FindBin(mean + width * 3); | |
876 | printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd); | |
877 | ||
878 | for (Int_t j=fillBegin; j <= fillEnd; ++j) | |
879 | { | |
880 | fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i)); | |
881 | } | |
882 | } | |
883 | ||
884 | for (Int_t i=1; i<=target->GetNbinsX(); ++i) | |
885 | { | |
886 | Float_t sum = 0; | |
887 | for (Int_t j=1; j<=10; ++j) | |
888 | sum += fineBinned->GetBinContent((i-1)*10 + j); | |
889 | target->SetBinContent(i, sum / 10); | |
890 | } | |
891 | ||
892 | delete[] binning; | |
893 | ||
894 | DrawComparison("Gaussian", mcTarget, correlationID, kFALSE); | |
0a173978 | 895 | } |