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