]>
Commit | Line | Data |
---|---|---|
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: AliUnfolding.cxx 31168 2009-02-23 15:18:45Z jgrosseo $ */ | |
17 | ||
18 | // This class allows 1-dimensional unfolding. | |
19 | // Methods that are implemented are chi2 minimization and bayesian unfolding. | |
20 | // | |
21 | // Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch | |
22 | ||
23 | #include "AliUnfolding.h" | |
24 | #include <TH1F.h> | |
25 | #include <TH2F.h> | |
26 | #include <TVirtualFitter.h> | |
27 | #include <TMath.h> | |
28 | #include <TCanvas.h> | |
29 | #include <TF1.h> | |
30 | ||
31 | TMatrixD* AliUnfolding::fgCorrelationMatrix = 0; | |
32 | TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0; | |
33 | TVectorD* AliUnfolding::fgCurrentESDVector = 0; | |
34 | TVectorD* AliUnfolding::fgEntropyAPriori = 0; | |
35 | ||
36 | TF1* AliUnfolding::fgFitFunction = 0; | |
37 | ||
38 | AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid; | |
39 | Int_t AliUnfolding::fgMaxInput = -1; // bins in measured histogram | |
40 | Int_t AliUnfolding::fgMaxParams = -1; // bins in unfolded histogram = number of fit params | |
41 | Float_t AliUnfolding::fgOverflowBinLimit = -1; | |
42 | ||
43 | AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1; | |
44 | Float_t AliUnfolding::fgRegularizationWeight = 10000; | |
45 | Int_t AliUnfolding::fgSkipBinsBegin = 0; | |
46 | Float_t AliUnfolding::fgMinuitStepSize = 0.1; // (usually not needed to be changed) step size in minimization | |
47 | Bool_t AliUnfolding::fgNormalizeInput = kFALSE; // normalize input spectrum | |
48 | ||
49 | Float_t AliUnfolding::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing) | |
50 | Int_t AliUnfolding::fgBayesianIterations = 10; // number of iterations in Bayesian method | |
51 | ||
52 | Bool_t AliUnfolding::fgDebug = kFALSE; | |
53 | ||
54 | ClassImp(AliUnfolding) | |
55 | ||
56 | //____________________________________________________________________ | |
57 | void AliUnfolding::SetUnfoldingMethod(MethodType methodType) | |
58 | { | |
59 | // set unfolding method | |
60 | fgMethodType = methodType; | |
61 | ||
62 | const char* name = 0; | |
63 | switch (methodType) | |
64 | { | |
65 | case kInvalid: name = "INVALID"; break; | |
66 | case kChi2Minimization: name = "Chi2 Minimization"; break; | |
67 | case kBayesian: name = "Bayesian unfolding"; break; | |
68 | case kFunction: name = "Functional fit"; break; | |
69 | } | |
70 | Printf("AliUnfolding::SetUnfoldingMethod: %s enabled.", name); | |
71 | } | |
72 | ||
73 | //____________________________________________________________________ | |
74 | void AliUnfolding::SetCreateOverflowBin(Float_t overflowBinLimit) | |
75 | { | |
76 | // enable the creation of a overflow bin that includes all statistics below the given limit | |
77 | ||
78 | fgOverflowBinLimit = overflowBinLimit; | |
79 | ||
80 | Printf("AliUnfolding::SetCreateOverflowBin: overflow bin limit set to %f", overflowBinLimit); | |
81 | } | |
82 | ||
83 | //____________________________________________________________________ | |
84 | void AliUnfolding::SetSkipBinsBegin(Int_t nBins) | |
85 | { | |
86 | // set number of skipped bins in regularization | |
87 | ||
88 | fgSkipBinsBegin = nBins; | |
89 | ||
90 | Printf("AliUnfolding::SetSkipBinsBegin: skipping %d bins at the beginning of the spectrum in the regularization.", fgSkipBinsBegin); | |
91 | } | |
92 | ||
93 | //____________________________________________________________________ | |
94 | void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded) | |
95 | { | |
96 | // set number of bins in the input (measured) distribution and in the unfolded distribution | |
97 | fgMaxInput = nMeasured; | |
98 | fgMaxParams = nUnfolded; | |
99 | ||
100 | Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded); | |
101 | } | |
102 | ||
103 | //____________________________________________________________________ | |
104 | void AliUnfolding::SetChi2Regularization(RegularizationType type, Float_t weight) | |
105 | { | |
106 | // | |
107 | // sets the parameters for chi2 minimization | |
108 | // | |
109 | ||
110 | fgRegularizationType = type; | |
111 | fgRegularizationWeight = weight; | |
112 | ||
113 | Printf("AliUnfolding::SetChi2Regularization --> Regularization set to %d with weight %f", (Int_t) type, weight); | |
114 | } | |
115 | ||
116 | //____________________________________________________________________ | |
117 | void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations) | |
118 | { | |
119 | // | |
120 | // sets the parameters for Bayesian unfolding | |
121 | // | |
122 | ||
123 | fgBayesianSmoothing = smoothing; | |
124 | fgBayesianIterations = nIterations; | |
125 | ||
126 | Printf("AliUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing); | |
127 | } | |
128 | ||
129 | //____________________________________________________________________ | |
130 | void AliUnfolding::SetFunction(TF1* function) | |
131 | { | |
132 | // set function for unfolding with a fit function | |
133 | ||
134 | fgFitFunction = function; | |
135 | ||
136 | Printf("AliUnfolding::SetFunction: Set fit function with %d parameters.", function->GetNpar()); | |
137 | } | |
138 | ||
139 | //____________________________________________________________________ | |
140 | Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check) | |
141 | { | |
142 | // unfolds with unfolding method fgMethodType | |
143 | // | |
144 | // parameters: | |
145 | // correlation: response matrix as measured vs. generated | |
146 | // efficiency: (optional) efficiency that is applied on the unfolded spectrum, i.e. it has to be in unfolded variables. If 0 no efficiency is applied. | |
147 | // measured: the measured spectrum | |
148 | // initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions. | |
149 | // result: target for the unfolded result | |
150 | // check: depends on the unfolding method, see comments in specific functions | |
151 | ||
152 | if (fgMaxInput == -1) | |
153 | { | |
154 | Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution"); | |
155 | fgMaxInput = measured->GetNbinsX(); | |
156 | } | |
157 | if (fgMaxParams == -1) | |
158 | { | |
159 | Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution"); | |
160 | fgMaxParams = measured->GetNbinsX(); | |
161 | } | |
162 | ||
163 | if (fgOverflowBinLimit > 0) | |
164 | CreateOverflowBin(correlation, measured); | |
165 | ||
166 | switch (fgMethodType) | |
167 | { | |
168 | case kInvalid: | |
169 | { | |
170 | Printf("AliUnfolding::Unfold: ERROR: Unfolding method not set. Use SetUnfoldingMethod. Exiting..."); | |
171 | return -1; | |
172 | } | |
173 | case kChi2Minimization: | |
174 | return UnfoldWithMinuit(correlation, efficiency, measured, initialConditions, result, check); | |
175 | case kBayesian: | |
176 | return UnfoldWithBayesian(correlation, efficiency, measured, initialConditions, result); | |
177 | case kFunction: | |
178 | return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result); | |
179 | } | |
180 | return -1; | |
181 | } | |
182 | ||
183 | //____________________________________________________________________ | |
184 | void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured) | |
185 | { | |
186 | // fill static variables needed for minuit fit | |
187 | ||
188 | if (!fgCorrelationMatrix) | |
189 | fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams); | |
190 | if (!fgCorrelationCovarianceMatrix) | |
191 | fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput); | |
192 | if (!fgCurrentESDVector) | |
193 | fgCurrentESDVector = new TVectorD(fgMaxInput); | |
194 | if (!fgEntropyAPriori) | |
195 | fgEntropyAPriori = new TVectorD(fgMaxParams); | |
196 | ||
197 | fgCorrelationMatrix->Zero(); | |
198 | fgCorrelationCovarianceMatrix->Zero(); | |
199 | fgCurrentESDVector->Zero(); | |
200 | fgEntropyAPriori->Zero(); | |
201 | ||
202 | // normalize correction for given nPart | |
203 | for (Int_t i=1; i<=correlation->GetNbinsX(); ++i) | |
204 | { | |
205 | Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY()); | |
206 | if (sum <= 0) | |
207 | continue; | |
208 | Float_t maxValue = 0; | |
209 | Int_t maxBin = -1; | |
210 | for (Int_t j=1; j<=correlation->GetNbinsY(); ++j) | |
211 | { | |
212 | // find most probably value | |
213 | if (maxValue < correlation->GetBinContent(i, j)) | |
214 | { | |
215 | maxValue = correlation->GetBinContent(i, j); | |
216 | maxBin = j; | |
217 | } | |
218 | ||
219 | // npart sum to 1 | |
220 | correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum); | |
221 | correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum); | |
222 | ||
223 | if (i <= fgMaxParams && j <= fgMaxInput) | |
224 | (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j); | |
225 | } | |
226 | ||
227 | //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin)); | |
228 | } | |
229 | ||
230 | //normalize measured | |
231 | if (fgNormalizeInput) | |
232 | measured->Scale(1.0 / measured->Integral()); | |
233 | ||
234 | for (Int_t i=0; i<fgMaxInput; ++i) | |
235 | { | |
236 | (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1); | |
237 | if (measured->GetBinError(i+1) > 0) | |
238 | (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1); | |
239 | ||
240 | if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7) | |
241 | (*fgCorrelationCovarianceMatrix)(i, i) = 0; | |
242 | //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i)); | |
243 | } | |
244 | } | |
245 | ||
246 | //____________________________________________________________________ | |
247 | Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check) | |
248 | { | |
249 | // | |
250 | // implementation of unfolding (internal function) | |
251 | // | |
252 | // unfolds <measured> using response from <correlation> and effiency <efficiency> | |
253 | // output is in <result> | |
254 | // <initialConditions> set the initial values for the minimization, if 0 <measured> is used | |
255 | // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed | |
256 | // | |
257 | // returns minuit status (0 = success), (-1 when check was set) | |
258 | // | |
259 | ||
260 | SetStaticVariables(correlation, measured); | |
261 | ||
262 | // Initialize TMinuit via generic fitter interface | |
263 | TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams); | |
264 | Double_t arglist[100]; | |
265 | ||
266 | // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way... | |
267 | arglist[0] = 0; | |
268 | minuit->ExecuteCommand("SET PRINT", arglist, 1); | |
269 | ||
270 | // however, enable warnings | |
271 | //minuit->ExecuteCommand("SET WAR", arglist, 0); | |
272 | ||
273 | // set minimization function | |
274 | minuit->SetFCN(Chi2Function); | |
275 | ||
276 | for (Int_t i=0; i<fgMaxParams; i++) | |
277 | (*fgEntropyAPriori)[i] = 1; | |
278 | ||
279 | if (!initialConditions) { | |
280 | initialConditions = measured; | |
281 | } else { | |
282 | Printf("AliUnfolding::UnfoldWithMinuit: Using different initial conditions..."); | |
283 | //new TCanvas; initialConditions->DrawCopy(); | |
284 | if (fgNormalizeInput) | |
285 | initialConditions->Scale(1.0 / initialConditions->Integral()); | |
286 | } | |
287 | ||
288 | // set initial conditions as a-priori distribution for MRX regularization | |
289 | for (Int_t i=0; i<fgMaxParams; i++) | |
290 | if (initialConditions->GetBinContent(i+1) > 0) | |
291 | (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1); | |
292 | ||
293 | Double_t* results = new Double_t[fgMaxParams+1]; | |
294 | for (Int_t i=0; i<fgMaxParams; ++i) | |
295 | { | |
296 | results[i] = initialConditions->GetBinContent(i+1); | |
297 | ||
298 | // minuit sees squared values to prevent it from going negative... | |
299 | results[i] = TMath::Sqrt(results[i]); | |
300 | ||
301 | minuit->SetParameter(i, Form("param%d", i), results[i], fgMinuitStepSize, 0, 0); | |
302 | } | |
303 | ||
304 | Int_t dummy = 0; | |
305 | Double_t chi2 = 0; | |
306 | Chi2Function(dummy, 0, chi2, results, 0); | |
307 | printf("AliUnfolding::UnfoldWithMinuit: Chi2 of initial parameters is = %f\n", chi2); | |
308 | ||
309 | if (check) | |
310 | { | |
311 | DrawGuess(results); | |
312 | delete[] results; | |
313 | return -1; | |
314 | } | |
315 | ||
316 | // first param is number of iterations, second is precision.... | |
317 | arglist[0] = 1e6; | |
318 | //arglist[1] = 1e-5; | |
319 | //minuit->ExecuteCommand("SCAN", arglist, 0); | |
320 | Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1); | |
321 | Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status); | |
322 | //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!"); | |
323 | //minuit->ExecuteCommand("MINOS", arglist, 0); | |
324 | ||
325 | for (Int_t i=0; i<fgMaxParams; ++i) | |
326 | { | |
327 | results[i] = minuit->GetParameter(i); | |
328 | Double_t value = results[i] * results[i]; | |
329 | // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i)) | |
330 | Double_t error = minuit->GetParError(i) * results[i]; | |
331 | ||
332 | if (efficiency) | |
333 | { | |
334 | if (efficiency->GetBinContent(i+1) > 0) | |
335 | { | |
336 | value /= efficiency->GetBinContent(i+1); | |
337 | error /= efficiency->GetBinContent(i+1); | |
338 | } | |
339 | else | |
340 | { | |
341 | value = 0; | |
342 | error = 0; | |
343 | } | |
344 | } | |
345 | ||
346 | result->SetBinContent(i+1, value); | |
347 | result->SetBinError(i+1, error); | |
348 | } | |
349 | ||
350 | if (fgDebug) | |
351 | DrawGuess(results); | |
352 | ||
353 | delete[] results; | |
354 | ||
355 | return status; | |
356 | } | |
357 | ||
358 | //____________________________________________________________________ | |
359 | Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult) | |
360 | { | |
361 | // | |
362 | // unfolds a spectrum using the Bayesian method | |
363 | // | |
364 | ||
365 | if (measured->Integral() <= 0) | |
366 | { | |
367 | Printf("AliUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty"); | |
368 | return -1; | |
369 | } | |
370 | ||
371 | const Int_t kStartBin = 0; | |
372 | ||
373 | Int_t kMaxM = fgMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis | |
374 | Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis | |
375 | ||
376 | // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4) | |
377 | const Double_t kConvergenceLimit = kMaxT * 1e-6; | |
378 | ||
379 | // store information in arrays, to increase processing speed (~ factor 5) | |
380 | Double_t* measuredCopy = new Double_t[kMaxM]; | |
381 | Double_t* measuredError = new Double_t[kMaxM]; | |
382 | Double_t* prior = new Double_t[kMaxT]; | |
383 | Double_t* result = new Double_t[kMaxT]; | |
384 | Double_t* efficiency = new Double_t[kMaxT]; | |
385 | ||
386 | Double_t** response = new Double_t*[kMaxT]; | |
387 | Double_t** inverseResponse = new Double_t*[kMaxT]; | |
388 | for (Int_t i=0; i<kMaxT; i++) | |
389 | { | |
390 | response[i] = new Double_t[kMaxM]; | |
391 | inverseResponse[i] = new Double_t[kMaxM]; | |
392 | } | |
393 | ||
394 | // for normalization | |
395 | Float_t measuredIntegral = measured->Integral(); | |
396 | for (Int_t m=0; m<kMaxM; m++) | |
397 | { | |
398 | measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral; | |
399 | measuredError[m] = measured->GetBinError(m+1) / measuredIntegral; | |
400 | ||
401 | for (Int_t t=0; t<kMaxT; t++) | |
402 | { | |
403 | response[t][m] = correlation->GetBinContent(t+1, m+1); | |
404 | inverseResponse[t][m] = 0; | |
405 | } | |
406 | } | |
407 | ||
408 | for (Int_t t=0; t<kMaxT; t++) | |
409 | { | |
410 | if (efficiency) | |
411 | { | |
412 | efficiency[t] = aEfficiency->GetBinContent(t+1); | |
413 | } | |
414 | else | |
415 | efficiency[t] = 1; | |
416 | ||
417 | prior[t] = measuredCopy[t]; | |
418 | result[t] = 0; | |
419 | } | |
420 | ||
421 | // pick prior distribution | |
422 | if (initialConditions) | |
423 | { | |
424 | printf("Using different starting conditions...\n"); | |
425 | // for normalization | |
426 | Float_t inputDistIntegral = initialConditions->Integral(); | |
427 | for (Int_t i=0; i<kMaxT; i++) | |
428 | prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral; | |
429 | } | |
430 | ||
431 | //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5); | |
432 | ||
433 | // unfold... | |
434 | for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations< 0; i++) | |
435 | { | |
436 | if (fgDebug) | |
437 | Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i); | |
438 | ||
439 | // calculate IR from Beyes theorem | |
440 | // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k) | |
441 | ||
442 | Double_t chi2Measured = 0; | |
443 | for (Int_t m=0; m<kMaxM; m++) | |
444 | { | |
445 | Float_t norm = 0; | |
446 | for (Int_t t = kStartBin; t<kMaxT; t++) | |
447 | norm += response[t][m] * prior[t]; | |
448 | ||
449 | // calc. chi2: (measured - response * prior) / error | |
450 | if (measuredError[m] > 0) | |
451 | { | |
452 | Double_t value = (measuredCopy[m] - norm) / measuredError[m]; | |
453 | chi2Measured += value * value; | |
454 | } | |
455 | ||
456 | if (norm > 0) | |
457 | { | |
458 | for (Int_t t = kStartBin; t<kMaxT; t++) | |
459 | inverseResponse[t][m] = response[t][m] * prior[t] / norm; | |
460 | } | |
461 | else | |
462 | { | |
463 | for (Int_t t = kStartBin; t<kMaxT; t++) | |
464 | inverseResponse[t][m] = 0; | |
465 | } | |
466 | } | |
467 | //Printf("chi2Measured of the last prior is %e", chi2Measured); | |
468 | ||
469 | for (Int_t t = kStartBin; t<kMaxT; t++) | |
470 | { | |
471 | Float_t value = 0; | |
472 | for (Int_t m=0; m<kMaxM; m++) | |
473 | value += inverseResponse[t][m] * measuredCopy[m]; | |
474 | ||
475 | if (efficiency[t] > 0) | |
476 | result[t] = value / efficiency[t]; | |
477 | else | |
478 | result[t] = 0; | |
479 | } | |
480 | ||
481 | // draw intermediate result | |
482 | /* | |
483 | for (Int_t t=0; t<kMaxT; t++) | |
484 | aResult->SetBinContent(t+1, result[t]); | |
485 | aResult->SetMarkerStyle(20+i); | |
486 | aResult->SetMarkerColor(2); | |
487 | aResult->DrawCopy("P SAME HIST"); | |
488 | */ | |
489 | ||
490 | Double_t chi2LastIter = 0; | |
491 | // regularization (simple smoothing) | |
492 | for (Int_t t=kStartBin; t<kMaxT; t++) | |
493 | { | |
494 | Float_t newValue = 0; | |
495 | ||
496 | // 0 bin excluded from smoothing | |
497 | if (t > kStartBin+2 && t<kMaxT-1) | |
498 | { | |
499 | Float_t average = (result[t-1] + result[t] + result[t+1]) / 3; | |
500 | ||
501 | // weight the average with the regularization parameter | |
502 | newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average; | |
503 | } | |
504 | else | |
505 | newValue = result[t]; | |
506 | ||
507 | // calculate chi2 (change from last iteration) | |
508 | if (prior[t] > 1e-5) | |
509 | { | |
510 | Double_t diff = (prior[t] - newValue) / prior[t]; | |
511 | chi2LastIter += diff * diff; | |
512 | } | |
513 | ||
514 | prior[t] = newValue; | |
515 | } | |
516 | //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter); | |
517 | //convergence->Fill(i+1, chi2LastIter); | |
518 | ||
519 | if (fgBayesianIterations < 0 && chi2LastIter < kConvergenceLimit) | |
520 | { | |
521 | Printf("AliUnfolding::UnfoldWithBayesian: Stopped Bayesian unfolding after %d iterations at chi2(change since last iteration) of %e; chi2Measured of the last prior is %e", i, chi2LastIter, chi2Measured); | |
522 | break; | |
523 | } | |
524 | } // end of iterations | |
525 | ||
526 | //new TCanvas; convergence->DrawCopy(); gPad->SetLogy(); | |
527 | //delete convergence; | |
528 | ||
529 | for (Int_t t=0; t<kMaxT; t++) | |
530 | aResult->SetBinContent(t+1, result[t]); | |
531 | ||
532 | delete[] measuredCopy; | |
533 | delete[] measuredError; | |
534 | delete[] prior; | |
535 | delete[] result; | |
536 | delete[] efficiency; | |
537 | ||
538 | for (Int_t i=0; i<kMaxT; i++) | |
539 | { | |
540 | delete[] response[i]; | |
541 | delete[] inverseResponse[i]; | |
542 | } | |
543 | delete[] response; | |
544 | delete[] inverseResponse; | |
545 | ||
546 | return 0; | |
547 | ||
548 | // ******** | |
549 | // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995 | |
550 | ||
551 | /*printf("Calculating covariance matrix. This may take some time...\n"); | |
552 | ||
553 | // check if this is the right one... | |
554 | TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX()); | |
555 | ||
556 | Int_t xBins = hInverseResponseBayes->GetNbinsX(); | |
557 | Int_t yBins = hInverseResponseBayes->GetNbinsY(); | |
558 | ||
559 | // calculate "unfolding matrix" Mij | |
560 | Float_t matrixM[251][251]; | |
561 | for (Int_t i=1; i<=xBins; i++) | |
562 | { | |
563 | for (Int_t j=1; j<=yBins; j++) | |
564 | { | |
565 | if (fCurrentEfficiency->GetBinContent(i) > 0) | |
566 | matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i); | |
567 | else | |
568 | matrixM[i-1][j-1] = 0; | |
569 | } | |
570 | } | |
571 | ||
572 | Float_t* vectorn = new Float_t[yBins]; | |
573 | for (Int_t j=1; j<=yBins; j++) | |
574 | vectorn[j-1] = fCurrentESD->GetBinContent(j); | |
575 | ||
576 | // first part of covariance matrix, depends on input distribution n(E) | |
577 | Float_t cov1[251][251]; | |
578 | ||
579 | Float_t nEvents = fCurrentESD->Integral(); // N | |
580 | ||
581 | xBins = 20; | |
582 | yBins = 20; | |
583 | ||
584 | for (Int_t k=0; k<xBins; k++) | |
585 | { | |
586 | printf("In Cov1: %d\n", k); | |
587 | for (Int_t l=0; l<yBins; l++) | |
588 | { | |
589 | cov1[k][l] = 0; | |
590 | ||
591 | // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N) | |
592 | for (Int_t j=0; j<yBins; j++) | |
593 | cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j] | |
594 | * (1.0 - vectorn[j] / nEvents); | |
595 | ||
596 | // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N | |
597 | for (Int_t i=0; i<yBins; i++) | |
598 | for (Int_t j=0; j<yBins; j++) | |
599 | { | |
600 | if (i == j) | |
601 | continue; | |
602 | cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i] | |
603 | * vectorn[j] / nEvents; | |
604 | } | |
605 | } | |
606 | } | |
607 | ||
608 | printf("Cov1 finished\n"); | |
609 | ||
610 | TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov"); | |
611 | cov->Reset(); | |
612 | ||
613 | for (Int_t i=1; i<=xBins; i++) | |
614 | for (Int_t j=1; j<=yBins; j++) | |
615 | cov->SetBinContent(i, j, cov1[i-1][j-1]); | |
616 | ||
617 | new TCanvas; | |
618 | cov->Draw("COLZ"); | |
619 | ||
620 | // second part of covariance matrix, depends on response matrix | |
621 | Float_t cov2[251][251]; | |
622 | ||
623 | // Cov[P(Er|Cu), P(Es|Cu)] term | |
624 | Float_t covTerm[100][100][100]; | |
625 | for (Int_t r=0; r<yBins; r++) | |
626 | for (Int_t u=0; u<xBins; u++) | |
627 | for (Int_t s=0; s<yBins; s++) | |
628 | { | |
629 | if (r == s) | |
630 | covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) | |
631 | * (1.0 - hResponse->GetBinContent(u+1, r+1)); | |
632 | else | |
633 | covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) | |
634 | * hResponse->GetBinContent(u+1, s+1); | |
635 | } | |
636 | ||
637 | for (Int_t k=0; k<xBins; k++) | |
638 | for (Int_t l=0; l<yBins; l++) | |
639 | { | |
640 | cov2[k][l] = 0; | |
641 | printf("In Cov2: %d %d\n", k, l); | |
642 | for (Int_t i=0; i<yBins; i++) | |
643 | for (Int_t j=0; j<yBins; j++) | |
644 | { | |
645 | //printf("In Cov2: %d %d %d %d\n", k, l, i, j); | |
646 | // calculate Cov(Mki, Mlj) = sum{ru},{su} ... | |
647 | Float_t tmpCov = 0; | |
648 | for (Int_t r=0; r<yBins; r++) | |
649 | for (Int_t u=0; u<xBins; u++) | |
650 | for (Int_t s=0; s<yBins; s++) | |
651 | { | |
652 | if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0 | |
653 | || hResponse->GetBinContent(u+1, i+1) == 0) | |
654 | continue; | |
655 | ||
656 | tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u) | |
657 | * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u) | |
658 | * covTerm[r][u][s]; | |
659 | } | |
660 | ||
661 | cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov; | |
662 | } | |
663 | } | |
664 | ||
665 | printf("Cov2 finished\n"); | |
666 | ||
667 | for (Int_t i=1; i<=xBins; i++) | |
668 | for (Int_t j=1; j<=yBins; j++) | |
669 | cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]); | |
670 | ||
671 | new TCanvas; | |
672 | cov->Draw("COLZ");*/ | |
673 | } | |
674 | ||
675 | //____________________________________________________________________ | |
676 | Double_t AliUnfolding::RegularizationPol0(TVectorD& params) | |
677 | { | |
678 | // homogenity term for minuit fitting | |
679 | // pure function of the parameters | |
680 | // prefers constant function (pol0) | |
681 | ||
682 | Double_t chi2 = 0; | |
683 | ||
684 | for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
685 | { | |
686 | Double_t right = params[i]; | |
687 | Double_t left = params[i-1]; | |
688 | ||
689 | if (right != 0) | |
690 | { | |
691 | Double_t diff = 1 - left / right; | |
692 | chi2 += diff * diff; | |
693 | } | |
694 | } | |
695 | ||
696 | return chi2 / 100.0; | |
697 | } | |
698 | ||
699 | //____________________________________________________________________ | |
700 | Double_t AliUnfolding::RegularizationPol1(TVectorD& params) | |
701 | { | |
702 | // homogenity term for minuit fitting | |
703 | // pure function of the parameters | |
704 | // prefers linear function (pol1) | |
705 | ||
706 | Double_t chi2 = 0; | |
707 | ||
708 | for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
709 | { | |
710 | if (params[i-1] == 0) | |
711 | continue; | |
712 | ||
713 | Double_t right = params[i]; | |
714 | Double_t middle = params[i-1]; | |
715 | Double_t left = params[i-2]; | |
716 | ||
717 | Double_t der1 = (right - middle); | |
718 | Double_t der2 = (middle - left); | |
719 | ||
720 | Double_t diff = (der1 - der2) / middle; | |
721 | ||
722 | chi2 += diff * diff; | |
723 | } | |
724 | ||
725 | return chi2; | |
726 | } | |
727 | ||
728 | //____________________________________________________________________ | |
729 | Double_t AliUnfolding::RegularizationLog(TVectorD& params) | |
730 | { | |
731 | // homogenity term for minuit fitting | |
732 | // pure function of the parameters | |
733 | // prefers linear function (pol1) | |
734 | ||
735 | Double_t chi2 = 0; | |
736 | ||
737 | for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
738 | { | |
739 | if (params[i-1] == 0) | |
740 | continue; | |
741 | ||
742 | Double_t right = log(params[i]); | |
743 | Double_t middle = log(params[i-1]); | |
744 | Double_t left = log(params[i-2]); | |
745 | ||
746 | Double_t der1 = (right - middle); | |
747 | Double_t der2 = (middle - left); | |
748 | ||
749 | Double_t diff = (der1 - der2) / middle; | |
750 | ||
751 | chi2 += diff * diff; | |
752 | } | |
753 | ||
754 | return chi2 * 100; | |
755 | } | |
756 | ||
757 | //____________________________________________________________________ | |
758 | Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params) | |
759 | { | |
760 | // homogenity term for minuit fitting | |
761 | // pure function of the parameters | |
762 | // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments, | |
763 | // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp. | |
764 | ||
765 | Double_t chi2 = 0; | |
766 | ||
767 | for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
768 | { | |
769 | Double_t right = params[i]; | |
770 | Double_t middle = params[i-1]; | |
771 | Double_t left = params[i-2]; | |
772 | ||
773 | Double_t der1 = (right - middle); | |
774 | Double_t der2 = (middle - left); | |
775 | ||
776 | Double_t diff = (der1 - der2); | |
777 | ||
778 | chi2 += diff * diff; | |
779 | } | |
780 | ||
781 | return chi2 * 1e4; | |
782 | } | |
783 | ||
784 | //____________________________________________________________________ | |
785 | Double_t AliUnfolding::RegularizationEntropy(TVectorD& params) | |
786 | { | |
787 | // homogenity term for minuit fitting | |
788 | // pure function of the parameters | |
789 | // calculates entropy, from | |
790 | // The method of reduced cross-entropy (M. Schmelling 1993) | |
791 | ||
792 | Double_t paramSum = 0; | |
793 | ||
794 | for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i) | |
795 | paramSum += params[i]; | |
796 | ||
797 | Double_t chi2 = 0; | |
798 | for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i) | |
799 | { | |
800 | //Double_t tmp = params[i] / paramSum; | |
801 | Double_t tmp = params[i]; | |
802 | if (tmp > 0 && (*fgEntropyAPriori)[i] > 0) | |
803 | { | |
804 | chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]); | |
805 | } | |
806 | } | |
807 | ||
808 | return 100.0 + chi2; | |
809 | } | |
810 | ||
811 | //____________________________________________________________________ | |
812 | void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t) | |
813 | { | |
814 | // | |
815 | // fit function for minuit | |
816 | // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix | |
817 | // | |
818 | ||
819 | // d | |
820 | static TVectorD paramsVector(fgMaxParams); | |
821 | for (Int_t i=0; i<fgMaxParams; ++i) | |
822 | paramsVector[i] = params[i] * params[i]; | |
823 | ||
824 | // calculate penalty factor | |
825 | Double_t penaltyVal = 0; | |
826 | switch (fgRegularizationType) | |
827 | { | |
828 | case kNone: break; | |
829 | case kPol0: penaltyVal = RegularizationPol0(paramsVector); break; | |
830 | case kPol1: penaltyVal = RegularizationPol1(paramsVector); break; | |
831 | case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break; | |
832 | case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break; | |
833 | case kLog: penaltyVal = RegularizationLog(paramsVector); break; | |
834 | } | |
835 | ||
836 | //if (penaltyVal > 1e10) | |
837 | // paramsVector2.Print(); | |
838 | ||
839 | penaltyVal *= fgRegularizationWeight; | |
840 | ||
841 | // Ad | |
842 | TVectorD measGuessVector(fgMaxInput); | |
843 | measGuessVector = (*fgCorrelationMatrix) * paramsVector; | |
844 | ||
845 | // Ad - m | |
846 | measGuessVector -= (*fgCurrentESDVector); | |
847 | ||
848 | //measGuessVector.Print(); | |
849 | ||
850 | TVectorD copy(measGuessVector); | |
851 | ||
852 | // (Ad - m) W | |
853 | // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used | |
854 | // normal way is like this: | |
855 | // measGuessVector *= (*fgCorrelationCovarianceMatrix); | |
856 | // optimized way like this: | |
857 | for (Int_t i=0; i<fgMaxInput; ++i) | |
858 | measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i); | |
859 | ||
860 | //measGuessVector.Print(); | |
861 | ||
862 | // (Ad - m) W (Ad - m) | |
863 | // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very | |
864 | // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit) | |
865 | Double_t chi2FromFit = measGuessVector * copy * 1e6; | |
866 | ||
867 | chi2 = chi2FromFit + penaltyVal; | |
868 | ||
869 | static Int_t callCount = 0; | |
870 | if ((callCount++ % 10000) == 0) | |
871 | Printf("AliUnfolding::Chi2Function: Iteration %d: %f %f --> %f", callCount, chi2FromFit, penaltyVal, chi2); | |
872 | } | |
873 | ||
874 | //____________________________________________________________________ | |
875 | void AliUnfolding::DrawGuess(Double_t *params) | |
876 | { | |
877 | // | |
878 | // draws residuals of solution suggested by params and effect of regularization | |
879 | // | |
880 | ||
881 | // d | |
882 | static TVectorD paramsVector(fgMaxParams); | |
883 | for (Int_t i=0; i<fgMaxParams; ++i) | |
884 | paramsVector[i] = params[i] * params[i]; | |
885 | ||
886 | // Ad | |
887 | TVectorD measGuessVector(fgMaxInput); | |
888 | measGuessVector = (*fgCorrelationMatrix) * paramsVector; | |
889 | ||
890 | // Ad - m | |
891 | measGuessVector -= (*fgCurrentESDVector); | |
892 | ||
893 | TVectorD copy(measGuessVector); | |
894 | //copy.Print(); | |
895 | ||
896 | // (Ad - m) W | |
897 | // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used | |
898 | // normal way is like this: | |
899 | // measGuessVector *= (*fgCorrelationCovarianceMatrix); | |
900 | // optimized way like this: | |
901 | for (Int_t i=0; i<fgMaxInput; ++i) | |
902 | measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i); | |
903 | ||
904 | // (Ad - m) W (Ad - m) | |
905 | // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very | |
906 | // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit) | |
907 | //Double_t chi2FromFit = measGuessVector * copy * 1e6; | |
908 | ||
909 | // draw residuals | |
910 | TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5); | |
911 | for (Int_t i=0; i<fgMaxInput; ++i) | |
912 | residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6); | |
913 | new TCanvas; residuals->DrawCopy(); gPad->SetLogy(); | |
914 | ||
915 | // draw penalty | |
916 | if (fgRegularizationType != kPol1) { | |
917 | Printf("Drawing not supported"); | |
918 | return; | |
919 | } | |
920 | ||
921 | TH1* penalty = new TH1F("penalty", "penalty", fgMaxParams, -0.5, fgMaxParams - 0.5); | |
922 | ||
923 | for (Int_t i=2+1; i<fgMaxParams; ++i) | |
924 | { | |
925 | if (params[i-1] == 0) | |
926 | continue; | |
927 | ||
928 | Double_t right = params[i]; | |
929 | Double_t middle = params[i-1]; | |
930 | Double_t left = params[i-2]; | |
931 | ||
932 | Double_t der1 = (right - middle); | |
933 | Double_t der2 = (middle - left); | |
934 | ||
935 | Double_t diff = (der1 - der2) / middle; | |
936 | ||
937 | penalty->SetBinContent(i-1, diff * diff); | |
938 | ||
939 | //Printf("%d %f %f %f %f", i-1, left, middle, right, diff); | |
940 | } | |
941 | new TCanvas; penalty->DrawCopy(); gPad->SetLogy(); | |
942 | } | |
943 | ||
944 | //____________________________________________________________________ | |
945 | void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3) | |
946 | { | |
947 | // | |
948 | // fit function for minuit | |
949 | // uses the TF1 stored in fgFitFunction | |
950 | // | |
951 | ||
952 | for (Int_t i=0; i<fgFitFunction->GetNpar(); i++) | |
953 | fgFitFunction->SetParameter(i, params[i]); | |
954 | ||
955 | Double_t* params2 = new Double_t[fgMaxParams]; | |
956 | ||
957 | for (Int_t i=0; i<fgMaxParams; ++i) | |
958 | params2[i] = fgFitFunction->Eval(i); | |
959 | ||
960 | Chi2Function(unused1, unused2, chi2, params2, unused3); | |
961 | ||
962 | delete[] params2; | |
963 | ||
964 | if (fgDebug) | |
965 | Printf("%f", chi2); | |
966 | } | |
967 | ||
968 | //____________________________________________________________________ | |
969 | Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult) | |
970 | { | |
971 | // | |
972 | // correct spectrum using minuit chi2 method applying a functional fit | |
973 | // | |
974 | ||
975 | if (!fgFitFunction) | |
976 | { | |
977 | Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting."); | |
978 | return -1; | |
979 | } | |
980 | ||
981 | SetChi2Regularization(kNone, 0); | |
982 | ||
983 | SetStaticVariables(correlation, measured); | |
984 | ||
985 | // Initialize TMinuit via generic fitter interface | |
986 | TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar()); | |
987 | ||
988 | minuit->SetFCN(TF1Function); | |
989 | for (Int_t i=0; i<fgFitFunction->GetNpar(); i++) | |
990 | { | |
991 | Double_t lower, upper; | |
992 | fgFitFunction->GetParLimits(i, lower, upper); | |
993 | minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper); | |
994 | } | |
995 | ||
996 | Double_t arglist[100]; | |
997 | arglist[0] = 0; | |
998 | minuit->ExecuteCommand("SET PRINT", arglist, 1); | |
999 | minuit->ExecuteCommand("MIGRAD", arglist, 0); | |
1000 | //minuit->ExecuteCommand("MINOS", arglist, 0); | |
1001 | ||
1002 | for (Int_t i=0; i<fgFitFunction->GetNpar(); i++) | |
1003 | fgFitFunction->SetParameter(i, minuit->GetParameter(i)); | |
1004 | ||
1005 | for (Int_t i=0; i<fgMaxParams; ++i) | |
1006 | { | |
1007 | Double_t value = fgFitFunction->Eval(i); | |
1008 | if (fgDebug) | |
1009 | Printf("%d : %f", i, value); | |
1010 | ||
1011 | if (efficiency) | |
1012 | { | |
1013 | if (efficiency->GetBinContent(i+1) > 0) | |
1014 | { | |
1015 | value /= efficiency->GetBinContent(i+1); | |
1016 | } | |
1017 | else | |
1018 | value = 0; | |
1019 | } | |
1020 | aResult->SetBinContent(i+1, value); | |
1021 | aResult->SetBinError(i+1, 0); | |
1022 | } | |
1023 | ||
1024 | return 0; | |
1025 | } | |
1026 | ||
1027 | //____________________________________________________________________ | |
1028 | void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured) | |
1029 | { | |
1030 | // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins | |
1031 | // The same limit is applied to the correlation | |
1032 | ||
1033 | Int_t lastBin = 0; | |
1034 | for (Int_t i=1; i<=measured->GetNbinsX(); ++i) | |
1035 | { | |
1036 | if (measured->GetBinContent(i) <= fgOverflowBinLimit) | |
1037 | { | |
1038 | lastBin = i; | |
1039 | break; | |
1040 | } | |
1041 | } | |
1042 | ||
1043 | if (lastBin == 0) | |
1044 | { | |
1045 | Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit); | |
1046 | return; | |
1047 | } | |
1048 | ||
1049 | Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin); | |
1050 | ||
1051 | TCanvas* canvas = 0; | |
1052 | ||
1053 | if (fgDebug) | |
1054 | { | |
1055 | canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800); | |
1056 | canvas->Divide(2, 2); | |
1057 | ||
1058 | canvas->cd(1); | |
1059 | measured->SetStats(kFALSE); | |
1060 | measured->DrawCopy(); | |
1061 | gPad->SetLogy(); | |
1062 | ||
1063 | canvas->cd(2); | |
1064 | correlation->SetStats(kFALSE); | |
1065 | correlation->DrawCopy("COLZ"); | |
1066 | } | |
1067 | ||
1068 | measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX())); | |
1069 | for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i) | |
1070 | { | |
1071 | measured->SetBinContent(i, 0); | |
1072 | measured->SetBinError(i, 0); | |
1073 | } | |
1074 | // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? | |
1075 | measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin))); | |
1076 | ||
1077 | Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin)); | |
1078 | ||
1079 | for (Int_t i=1; i<=correlation->GetNbinsX(); ++i) | |
1080 | { | |
1081 | correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY())); | |
1082 | // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? | |
1083 | correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin))); | |
1084 | ||
1085 | for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j) | |
1086 | { | |
1087 | correlation->SetBinContent(i, j, 0); | |
1088 | correlation->SetBinError(i, j, 0); | |
1089 | } | |
1090 | } | |
1091 | ||
1092 | Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!"); | |
1093 | ||
1094 | if (canvas) | |
1095 | { | |
1096 | canvas->cd(3); | |
1097 | measured->DrawCopy(); | |
1098 | gPad->SetLogy(); | |
1099 | ||
1100 | canvas->cd(4); | |
1101 | correlation->DrawCopy("COLZ"); | |
1102 | } | |
1103 | } |