]>
Commit | Line | Data |
---|---|---|
19442b86 | 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; | |
95e970ca | 32 | TMatrixD* AliUnfolding::fgCorrelationMatrixSquared = 0; |
19442b86 | 33 | TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0; |
34 | TVectorD* AliUnfolding::fgCurrentESDVector = 0; | |
35 | TVectorD* AliUnfolding::fgEntropyAPriori = 0; | |
95e970ca | 36 | TVectorD* AliUnfolding::fgEfficiency = 0; |
37 | TVectorD* AliUnfolding::fgBinWidths = 0; | |
19442b86 | 38 | |
39 | TF1* AliUnfolding::fgFitFunction = 0; | |
40 | ||
41 | AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid; | |
42 | Int_t AliUnfolding::fgMaxInput = -1; // bins in measured histogram | |
43 | Int_t AliUnfolding::fgMaxParams = -1; // bins in unfolded histogram = number of fit params | |
44 | Float_t AliUnfolding::fgOverflowBinLimit = -1; | |
45 | ||
46 | AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1; | |
47 | Float_t AliUnfolding::fgRegularizationWeight = 10000; | |
48 | Int_t AliUnfolding::fgSkipBinsBegin = 0; | |
49 | Float_t AliUnfolding::fgMinuitStepSize = 0.1; // (usually not needed to be changed) step size in minimization | |
95e970ca | 50 | Bool_t AliUnfolding::fgMinimumInitialValue = kFALSE; // set all initial values at least to the smallest value among the initial values |
51 | Float_t AliUnfolding::fgMinimumInitialValueFix = -1; | |
19442b86 | 52 | Bool_t AliUnfolding::fgNormalizeInput = kFALSE; // normalize input spectrum |
95e970ca | 53 | Float_t AliUnfolding::fgNotFoundEvents = 0; |
54 | Bool_t AliUnfolding::fgSkipBin0InChi2 = kFALSE; | |
19442b86 | 55 | |
56 | Float_t AliUnfolding::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing) | |
651e2088 | 57 | Int_t AliUnfolding::fgBayesianIterations = 10; // number of iterations in Bayesian method |
19442b86 | 58 | |
59 | Bool_t AliUnfolding::fgDebug = kFALSE; | |
60 | ||
95e970ca | 61 | Int_t AliUnfolding::fgCallCount = 0; |
62 | ||
19442b86 | 63 | ClassImp(AliUnfolding) |
64 | ||
65 | //____________________________________________________________________ | |
66 | void AliUnfolding::SetUnfoldingMethod(MethodType methodType) | |
67 | { | |
68 | // set unfolding method | |
69 | fgMethodType = methodType; | |
70 | ||
71 | const char* name = 0; | |
72 | switch (methodType) | |
73 | { | |
74 | case kInvalid: name = "INVALID"; break; | |
75 | case kChi2Minimization: name = "Chi2 Minimization"; break; | |
76 | case kBayesian: name = "Bayesian unfolding"; break; | |
77 | case kFunction: name = "Functional fit"; break; | |
78 | } | |
79 | Printf("AliUnfolding::SetUnfoldingMethod: %s enabled.", name); | |
80 | } | |
81 | ||
82 | //____________________________________________________________________ | |
83 | void AliUnfolding::SetCreateOverflowBin(Float_t overflowBinLimit) | |
84 | { | |
85 | // enable the creation of a overflow bin that includes all statistics below the given limit | |
86 | ||
87 | fgOverflowBinLimit = overflowBinLimit; | |
88 | ||
89 | Printf("AliUnfolding::SetCreateOverflowBin: overflow bin limit set to %f", overflowBinLimit); | |
90 | } | |
91 | ||
92 | //____________________________________________________________________ | |
93 | void AliUnfolding::SetSkipBinsBegin(Int_t nBins) | |
94 | { | |
95 | // set number of skipped bins in regularization | |
96 | ||
97 | fgSkipBinsBegin = nBins; | |
98 | ||
99 | Printf("AliUnfolding::SetSkipBinsBegin: skipping %d bins at the beginning of the spectrum in the regularization.", fgSkipBinsBegin); | |
100 | } | |
101 | ||
102 | //____________________________________________________________________ | |
103 | void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded) | |
104 | { | |
105 | // set number of bins in the input (measured) distribution and in the unfolded distribution | |
106 | fgMaxInput = nMeasured; | |
107 | fgMaxParams = nUnfolded; | |
108 | ||
95e970ca | 109 | if (fgCorrelationMatrix) |
110 | { | |
111 | delete fgCorrelationMatrix; | |
112 | fgCorrelationMatrix = 0; | |
113 | } | |
114 | if (fgCorrelationMatrixSquared) | |
115 | { | |
116 | fgCorrelationMatrixSquared = 0; | |
117 | delete fgCorrelationMatrixSquared; | |
118 | } | |
119 | if (fgCorrelationCovarianceMatrix) | |
120 | { | |
121 | delete fgCorrelationCovarianceMatrix; | |
122 | fgCorrelationCovarianceMatrix = 0; | |
123 | } | |
124 | if (fgCurrentESDVector) | |
125 | { | |
126 | delete fgCurrentESDVector; | |
127 | fgCurrentESDVector = 0; | |
128 | } | |
129 | if (fgEntropyAPriori) | |
130 | { | |
131 | delete fgEntropyAPriori; | |
132 | fgEntropyAPriori = 0; | |
133 | } | |
134 | if (fgEfficiency) | |
135 | { | |
136 | delete fgEfficiency; | |
137 | fgEfficiency = 0; | |
138 | } | |
139 | if (fgBinWidths) | |
140 | { | |
141 | delete fgBinWidths; | |
142 | fgBinWidths = 0; | |
143 | } | |
144 | ||
19442b86 | 145 | Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded); |
146 | } | |
147 | ||
148 | //____________________________________________________________________ | |
149 | void AliUnfolding::SetChi2Regularization(RegularizationType type, Float_t weight) | |
150 | { | |
151 | // | |
152 | // sets the parameters for chi2 minimization | |
153 | // | |
154 | ||
155 | fgRegularizationType = type; | |
156 | fgRegularizationWeight = weight; | |
157 | ||
158 | Printf("AliUnfolding::SetChi2Regularization --> Regularization set to %d with weight %f", (Int_t) type, weight); | |
159 | } | |
160 | ||
161 | //____________________________________________________________________ | |
162 | void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations) | |
163 | { | |
164 | // | |
165 | // sets the parameters for Bayesian unfolding | |
166 | // | |
167 | ||
168 | fgBayesianSmoothing = smoothing; | |
169 | fgBayesianIterations = nIterations; | |
170 | ||
171 | Printf("AliUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing); | |
172 | } | |
173 | ||
174 | //____________________________________________________________________ | |
175 | void AliUnfolding::SetFunction(TF1* function) | |
176 | { | |
177 | // set function for unfolding with a fit function | |
178 | ||
179 | fgFitFunction = function; | |
180 | ||
181 | Printf("AliUnfolding::SetFunction: Set fit function with %d parameters.", function->GetNpar()); | |
182 | } | |
183 | ||
184 | //____________________________________________________________________ | |
185 | Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check) | |
186 | { | |
187 | // unfolds with unfolding method fgMethodType | |
188 | // | |
189 | // parameters: | |
190 | // correlation: response matrix as measured vs. generated | |
191 | // 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. | |
192 | // measured: the measured spectrum | |
193 | // initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions. | |
194 | // result: target for the unfolded result | |
195 | // check: depends on the unfolding method, see comments in specific functions | |
95e970ca | 196 | // |
197 | // return code: see UnfoldWithMinuit/UnfoldWithBayesian/UnfoldWithFunction | |
19442b86 | 198 | |
199 | if (fgMaxInput == -1) | |
200 | { | |
201 | Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution"); | |
202 | fgMaxInput = measured->GetNbinsX(); | |
203 | } | |
204 | if (fgMaxParams == -1) | |
205 | { | |
206 | Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution"); | |
207 | fgMaxParams = measured->GetNbinsX(); | |
208 | } | |
209 | ||
210 | if (fgOverflowBinLimit > 0) | |
211 | CreateOverflowBin(correlation, measured); | |
212 | ||
213 | switch (fgMethodType) | |
214 | { | |
215 | case kInvalid: | |
216 | { | |
217 | Printf("AliUnfolding::Unfold: ERROR: Unfolding method not set. Use SetUnfoldingMethod. Exiting..."); | |
218 | return -1; | |
219 | } | |
220 | case kChi2Minimization: | |
221 | return UnfoldWithMinuit(correlation, efficiency, measured, initialConditions, result, check); | |
222 | case kBayesian: | |
223 | return UnfoldWithBayesian(correlation, efficiency, measured, initialConditions, result); | |
224 | case kFunction: | |
225 | return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result); | |
226 | } | |
227 | return -1; | |
228 | } | |
229 | ||
230 | //____________________________________________________________________ | |
95e970ca | 231 | void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured, TH1* efficiency) |
19442b86 | 232 | { |
233 | // fill static variables needed for minuit fit | |
234 | ||
235 | if (!fgCorrelationMatrix) | |
236 | fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams); | |
95e970ca | 237 | if (!fgCorrelationMatrixSquared) |
238 | fgCorrelationMatrixSquared = new TMatrixD(fgMaxInput, fgMaxParams); | |
19442b86 | 239 | if (!fgCorrelationCovarianceMatrix) |
240 | fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput); | |
241 | if (!fgCurrentESDVector) | |
242 | fgCurrentESDVector = new TVectorD(fgMaxInput); | |
243 | if (!fgEntropyAPriori) | |
244 | fgEntropyAPriori = new TVectorD(fgMaxParams); | |
95e970ca | 245 | if (!fgEfficiency) |
246 | fgEfficiency = new TVectorD(fgMaxParams); | |
247 | if (!fgBinWidths) | |
248 | fgBinWidths = new TVectorD(fgMaxParams); | |
249 | ||
19442b86 | 250 | fgCorrelationMatrix->Zero(); |
251 | fgCorrelationCovarianceMatrix->Zero(); | |
252 | fgCurrentESDVector->Zero(); | |
253 | fgEntropyAPriori->Zero(); | |
254 | ||
255 | // normalize correction for given nPart | |
256 | for (Int_t i=1; i<=correlation->GetNbinsX(); ++i) | |
257 | { | |
258 | Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY()); | |
259 | if (sum <= 0) | |
260 | continue; | |
261 | Float_t maxValue = 0; | |
262 | Int_t maxBin = -1; | |
263 | for (Int_t j=1; j<=correlation->GetNbinsY(); ++j) | |
264 | { | |
265 | // find most probably value | |
266 | if (maxValue < correlation->GetBinContent(i, j)) | |
267 | { | |
268 | maxValue = correlation->GetBinContent(i, j); | |
269 | maxBin = j; | |
270 | } | |
271 | ||
272 | // npart sum to 1 | |
95e970ca | 273 | correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);// * correlation->GetXaxis()->GetBinWidth(i)); |
19442b86 | 274 | correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum); |
275 | ||
276 | if (i <= fgMaxParams && j <= fgMaxInput) | |
95e970ca | 277 | { |
19442b86 | 278 | (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j); |
95e970ca | 279 | (*fgCorrelationMatrixSquared)(j-1, i-1) = correlation->GetBinContent(i, j) * correlation->GetBinContent(i, j); |
280 | } | |
19442b86 | 281 | } |
282 | ||
283 | //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin)); | |
284 | } | |
285 | ||
286 | //normalize measured | |
95e970ca | 287 | Float_t smallestError = 1; |
19442b86 | 288 | if (fgNormalizeInput) |
95e970ca | 289 | { |
290 | Float_t sumMeasured = measured->Integral(); | |
291 | measured->Scale(1.0 / sumMeasured); | |
292 | smallestError /= sumMeasured; | |
293 | } | |
19442b86 | 294 | |
295 | for (Int_t i=0; i<fgMaxInput; ++i) | |
296 | { | |
297 | (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1); | |
298 | if (measured->GetBinError(i+1) > 0) | |
95e970ca | 299 | { |
19442b86 | 300 | (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1); |
95e970ca | 301 | } |
302 | else // in this case put error of 1, otherwise 0 bins are not added to the chi2... | |
303 | (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / smallestError / smallestError; | |
19442b86 | 304 | |
305 | if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7) | |
306 | (*fgCorrelationCovarianceMatrix)(i, i) = 0; | |
307 | //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i)); | |
308 | } | |
95e970ca | 309 | |
310 | // efficiency is expected to match bin width of result | |
311 | for (Int_t i=0; i<fgMaxParams; ++i) | |
312 | { | |
313 | (*fgEfficiency)(i) = efficiency->GetBinContent(i+1); | |
314 | (*fgBinWidths)(i) = efficiency->GetXaxis()->GetBinWidth(i+1); | |
315 | } | |
19442b86 | 316 | } |
317 | ||
318 | //____________________________________________________________________ | |
319 | Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check) | |
320 | { | |
321 | // | |
322 | // implementation of unfolding (internal function) | |
323 | // | |
324 | // unfolds <measured> using response from <correlation> and effiency <efficiency> | |
325 | // output is in <result> | |
326 | // <initialConditions> set the initial values for the minimization, if 0 <measured> is used | |
95e970ca | 327 | // negative values in initialConditions mean that the given parameter is fixed to the absolute of the value |
19442b86 | 328 | // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed |
329 | // | |
330 | // returns minuit status (0 = success), (-1 when check was set) | |
331 | // | |
332 | ||
95e970ca | 333 | SetStaticVariables(correlation, measured, efficiency); |
19442b86 | 334 | |
335 | // Initialize TMinuit via generic fitter interface | |
95e970ca | 336 | Int_t params = fgMaxParams; |
337 | if (fgNotFoundEvents > 0) | |
338 | params++; | |
339 | ||
340 | TVirtualFitter *minuit = TVirtualFitter::Fitter(0, params); | |
19442b86 | 341 | Double_t arglist[100]; |
342 | ||
343 | // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way... | |
344 | arglist[0] = 0; | |
345 | minuit->ExecuteCommand("SET PRINT", arglist, 1); | |
346 | ||
347 | // however, enable warnings | |
348 | //minuit->ExecuteCommand("SET WAR", arglist, 0); | |
349 | ||
350 | // set minimization function | |
351 | minuit->SetFCN(Chi2Function); | |
352 | ||
353 | for (Int_t i=0; i<fgMaxParams; i++) | |
354 | (*fgEntropyAPriori)[i] = 1; | |
355 | ||
95e970ca | 356 | // set initial conditions as a-priori distribution for MRX regularization |
357 | /* | |
358 | for (Int_t i=0; i<fgMaxParams; i++) | |
359 | if (initialConditions && initialConditions->GetBinContent(i+1) > 0) | |
360 | (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1); | |
361 | */ | |
362 | ||
19442b86 | 363 | if (!initialConditions) { |
364 | initialConditions = measured; | |
365 | } else { | |
366 | Printf("AliUnfolding::UnfoldWithMinuit: Using different initial conditions..."); | |
367 | //new TCanvas; initialConditions->DrawCopy(); | |
368 | if (fgNormalizeInput) | |
369 | initialConditions->Scale(1.0 / initialConditions->Integral()); | |
370 | } | |
371 | ||
95e970ca | 372 | // extract minimum value from initial conditions (if we set a value to 0 it will stay 0) |
373 | Float_t minValue = 1e100; | |
374 | if (fgMinimumInitialValueFix < 0) | |
375 | { | |
376 | for (Int_t i=0; i<fgMaxParams; ++i) | |
377 | { | |
378 | Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1)); | |
379 | if (initialConditions->GetBinContent(bin) > 0) | |
380 | minValue = TMath::Min(minValue, (Float_t) initialConditions->GetBinContent(bin)); | |
381 | } | |
382 | } | |
383 | else | |
384 | minValue = fgMinimumInitialValueFix; | |
385 | ||
19442b86 | 386 | Double_t* results = new Double_t[fgMaxParams+1]; |
387 | for (Int_t i=0; i<fgMaxParams; ++i) | |
388 | { | |
95e970ca | 389 | Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1)); |
390 | results[i] = initialConditions->GetBinContent(bin); | |
19442b86 | 391 | |
95e970ca | 392 | Bool_t fix = kFALSE; |
393 | if (results[i] < 0) | |
394 | { | |
395 | fix = kTRUE; | |
396 | results[i] = -results[i]; | |
397 | } | |
398 | ||
399 | if (!fix && fgMinimumInitialValue && results[i] < minValue) | |
400 | results[i] = minValue; | |
401 | ||
19442b86 | 402 | // minuit sees squared values to prevent it from going negative... |
403 | results[i] = TMath::Sqrt(results[i]); | |
404 | ||
95e970ca | 405 | minuit->SetParameter(i, Form("param%d", i), results[i], (fix) ? 0 : fgMinuitStepSize, 0, 0); |
406 | } | |
407 | if (fgNotFoundEvents > 0) | |
408 | { | |
409 | results[fgMaxParams] = efficiency->GetBinContent(1); | |
410 | minuit->SetParameter(fgMaxParams, "vtx0", results[fgMaxParams], fgMinuitStepSize / 100, 0.01, 0.80); | |
19442b86 | 411 | } |
412 | ||
413 | Int_t dummy = 0; | |
414 | Double_t chi2 = 0; | |
415 | Chi2Function(dummy, 0, chi2, results, 0); | |
416 | printf("AliUnfolding::UnfoldWithMinuit: Chi2 of initial parameters is = %f\n", chi2); | |
417 | ||
418 | if (check) | |
419 | { | |
420 | DrawGuess(results); | |
421 | delete[] results; | |
422 | return -1; | |
423 | } | |
424 | ||
425 | // first param is number of iterations, second is precision.... | |
426 | arglist[0] = 1e6; | |
427 | //arglist[1] = 1e-5; | |
428 | //minuit->ExecuteCommand("SCAN", arglist, 0); | |
429 | Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1); | |
430 | Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status); | |
431 | //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!"); | |
432 | //minuit->ExecuteCommand("MINOS", arglist, 0); | |
433 | ||
95e970ca | 434 | if (fgNotFoundEvents > 0) |
435 | { | |
436 | results[fgMaxParams] = minuit->GetParameter(fgMaxParams); | |
437 | Printf("Efficiency for bin 0 changed from %f to %f", efficiency->GetBinContent(1), results[fgMaxParams]); | |
438 | efficiency->SetBinContent(1, results[fgMaxParams]); | |
439 | } | |
440 | ||
19442b86 | 441 | for (Int_t i=0; i<fgMaxParams; ++i) |
442 | { | |
443 | results[i] = minuit->GetParameter(i); | |
444 | Double_t value = results[i] * results[i]; | |
445 | // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i)) | |
b203a518 | 446 | Double_t error = 0; |
447 | if (TMath::IsNaN(minuit->GetParError(i))) | |
448 | Printf("WARNING: Parameter %d error is nan", i); | |
449 | else | |
450 | error = minuit->GetParError(i) * results[i]; | |
19442b86 | 451 | |
452 | if (efficiency) | |
95e970ca | 453 | { |
19442b86 | 454 | if (efficiency->GetBinContent(i+1) > 0) |
455 | { | |
95e970ca | 456 | value /= efficiency->GetBinContent(i+1); |
457 | error /= efficiency->GetBinContent(i+1); | |
19442b86 | 458 | } |
459 | else | |
460 | { | |
95e970ca | 461 | value = 0; |
462 | error = 0; | |
19442b86 | 463 | } |
464 | } | |
465 | ||
466 | result->SetBinContent(i+1, value); | |
467 | result->SetBinError(i+1, error); | |
468 | } | |
469 | ||
95e970ca | 470 | fgCallCount = 0; |
471 | Chi2Function(dummy, 0, chi2, results, 0); | |
b203a518 | 472 | Printf("AliUnfolding::UnfoldWithMinuit: Chi2 of final parameters is = %f", chi2); |
95e970ca | 473 | |
19442b86 | 474 | if (fgDebug) |
475 | DrawGuess(results); | |
476 | ||
477 | delete[] results; | |
478 | ||
479 | return status; | |
480 | } | |
481 | ||
482 | //____________________________________________________________________ | |
483 | Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult) | |
484 | { | |
485 | // | |
486 | // unfolds a spectrum using the Bayesian method | |
487 | // | |
488 | ||
489 | if (measured->Integral() <= 0) | |
490 | { | |
491 | Printf("AliUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty"); | |
492 | return -1; | |
493 | } | |
494 | ||
495 | const Int_t kStartBin = 0; | |
496 | ||
497 | Int_t kMaxM = fgMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis | |
498 | Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis | |
499 | ||
500 | // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4) | |
501 | const Double_t kConvergenceLimit = kMaxT * 1e-6; | |
502 | ||
503 | // store information in arrays, to increase processing speed (~ factor 5) | |
504 | Double_t* measuredCopy = new Double_t[kMaxM]; | |
505 | Double_t* measuredError = new Double_t[kMaxM]; | |
506 | Double_t* prior = new Double_t[kMaxT]; | |
507 | Double_t* result = new Double_t[kMaxT]; | |
508 | Double_t* efficiency = new Double_t[kMaxT]; | |
95e970ca | 509 | Double_t* binWidths = new Double_t[kMaxT]; |
19442b86 | 510 | |
511 | Double_t** response = new Double_t*[kMaxT]; | |
512 | Double_t** inverseResponse = new Double_t*[kMaxT]; | |
513 | for (Int_t i=0; i<kMaxT; i++) | |
514 | { | |
515 | response[i] = new Double_t[kMaxM]; | |
516 | inverseResponse[i] = new Double_t[kMaxM]; | |
517 | } | |
518 | ||
519 | // for normalization | |
520 | Float_t measuredIntegral = measured->Integral(); | |
521 | for (Int_t m=0; m<kMaxM; m++) | |
522 | { | |
523 | measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral; | |
524 | measuredError[m] = measured->GetBinError(m+1) / measuredIntegral; | |
525 | ||
526 | for (Int_t t=0; t<kMaxT; t++) | |
527 | { | |
528 | response[t][m] = correlation->GetBinContent(t+1, m+1); | |
529 | inverseResponse[t][m] = 0; | |
530 | } | |
531 | } | |
532 | ||
533 | for (Int_t t=0; t<kMaxT; t++) | |
534 | { | |
651e2088 | 535 | if (aEfficiency) |
19442b86 | 536 | { |
537 | efficiency[t] = aEfficiency->GetBinContent(t+1); | |
538 | } | |
539 | else | |
540 | efficiency[t] = 1; | |
541 | ||
542 | prior[t] = measuredCopy[t]; | |
543 | result[t] = 0; | |
95e970ca | 544 | binWidths[t] = aResult->GetXaxis()->GetBinWidth(t+1); |
19442b86 | 545 | } |
546 | ||
547 | // pick prior distribution | |
548 | if (initialConditions) | |
549 | { | |
550 | printf("Using different starting conditions...\n"); | |
551 | // for normalization | |
552 | Float_t inputDistIntegral = initialConditions->Integral(); | |
553 | for (Int_t i=0; i<kMaxT; i++) | |
554 | prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral; | |
555 | } | |
556 | ||
557 | //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5); | |
558 | ||
95e970ca | 559 | //new TCanvas; |
19442b86 | 560 | // unfold... |
95e970ca | 561 | for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations < 0; i++) |
19442b86 | 562 | { |
563 | if (fgDebug) | |
564 | Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i); | |
565 | ||
95e970ca | 566 | // calculate IR from Bayes theorem |
19442b86 | 567 | // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k) |
568 | ||
569 | Double_t chi2Measured = 0; | |
570 | for (Int_t m=0; m<kMaxM; m++) | |
571 | { | |
572 | Float_t norm = 0; | |
573 | for (Int_t t = kStartBin; t<kMaxT; t++) | |
574 | norm += response[t][m] * prior[t]; | |
575 | ||
576 | // calc. chi2: (measured - response * prior) / error | |
577 | if (measuredError[m] > 0) | |
578 | { | |
579 | Double_t value = (measuredCopy[m] - norm) / measuredError[m]; | |
580 | chi2Measured += value * value; | |
581 | } | |
582 | ||
583 | if (norm > 0) | |
584 | { | |
585 | for (Int_t t = kStartBin; t<kMaxT; t++) | |
586 | inverseResponse[t][m] = response[t][m] * prior[t] / norm; | |
587 | } | |
588 | else | |
589 | { | |
590 | for (Int_t t = kStartBin; t<kMaxT; t++) | |
591 | inverseResponse[t][m] = 0; | |
592 | } | |
593 | } | |
594 | //Printf("chi2Measured of the last prior is %e", chi2Measured); | |
595 | ||
596 | for (Int_t t = kStartBin; t<kMaxT; t++) | |
597 | { | |
598 | Float_t value = 0; | |
599 | for (Int_t m=0; m<kMaxM; m++) | |
600 | value += inverseResponse[t][m] * measuredCopy[m]; | |
601 | ||
602 | if (efficiency[t] > 0) | |
603 | result[t] = value / efficiency[t]; | |
604 | else | |
605 | result[t] = 0; | |
606 | } | |
607 | ||
95e970ca | 608 | /* |
19442b86 | 609 | // draw intermediate result |
19442b86 | 610 | for (Int_t t=0; t<kMaxT; t++) |
95e970ca | 611 | { |
19442b86 | 612 | aResult->SetBinContent(t+1, result[t]); |
95e970ca | 613 | } |
614 | aResult->SetMarkerStyle(24+i); | |
19442b86 | 615 | aResult->SetMarkerColor(2); |
95e970ca | 616 | aResult->DrawCopy((i == 0) ? "P" : "PSAME"); |
19442b86 | 617 | */ |
95e970ca | 618 | |
19442b86 | 619 | Double_t chi2LastIter = 0; |
620 | // regularization (simple smoothing) | |
621 | for (Int_t t=kStartBin; t<kMaxT; t++) | |
622 | { | |
623 | Float_t newValue = 0; | |
624 | ||
625 | // 0 bin excluded from smoothing | |
626 | if (t > kStartBin+2 && t<kMaxT-1) | |
627 | { | |
95e970ca | 628 | Float_t average = (result[t-1] / binWidths[t-1] + result[t] / binWidths[t] + result[t+1] / binWidths[t+1]) / 3 * binWidths[t]; |
19442b86 | 629 | |
630 | // weight the average with the regularization parameter | |
631 | newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average; | |
632 | } | |
633 | else | |
634 | newValue = result[t]; | |
635 | ||
636 | // calculate chi2 (change from last iteration) | |
637 | if (prior[t] > 1e-5) | |
638 | { | |
639 | Double_t diff = (prior[t] - newValue) / prior[t]; | |
640 | chi2LastIter += diff * diff; | |
641 | } | |
642 | ||
643 | prior[t] = newValue; | |
644 | } | |
645 | //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter); | |
646 | //convergence->Fill(i+1, chi2LastIter); | |
647 | ||
648 | if (fgBayesianIterations < 0 && chi2LastIter < kConvergenceLimit) | |
649 | { | |
650 | 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); | |
651 | break; | |
652 | } | |
653 | } // end of iterations | |
654 | ||
655 | //new TCanvas; convergence->DrawCopy(); gPad->SetLogy(); | |
656 | //delete convergence; | |
657 | ||
95e970ca | 658 | Float_t factor = 1; |
659 | if (!fgNormalizeInput) | |
660 | factor = measuredIntegral; | |
19442b86 | 661 | for (Int_t t=0; t<kMaxT; t++) |
95e970ca | 662 | aResult->SetBinContent(t+1, result[t] * factor); |
19442b86 | 663 | |
664 | delete[] measuredCopy; | |
665 | delete[] measuredError; | |
666 | delete[] prior; | |
667 | delete[] result; | |
668 | delete[] efficiency; | |
95e970ca | 669 | delete[] binWidths; |
19442b86 | 670 | |
671 | for (Int_t i=0; i<kMaxT; i++) | |
672 | { | |
673 | delete[] response[i]; | |
674 | delete[] inverseResponse[i]; | |
675 | } | |
676 | delete[] response; | |
677 | delete[] inverseResponse; | |
678 | ||
679 | return 0; | |
680 | ||
681 | // ******** | |
682 | // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995 | |
683 | ||
684 | /*printf("Calculating covariance matrix. This may take some time...\n"); | |
685 | ||
686 | // check if this is the right one... | |
687 | TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX()); | |
688 | ||
689 | Int_t xBins = hInverseResponseBayes->GetNbinsX(); | |
690 | Int_t yBins = hInverseResponseBayes->GetNbinsY(); | |
691 | ||
692 | // calculate "unfolding matrix" Mij | |
693 | Float_t matrixM[251][251]; | |
694 | for (Int_t i=1; i<=xBins; i++) | |
695 | { | |
696 | for (Int_t j=1; j<=yBins; j++) | |
697 | { | |
698 | if (fCurrentEfficiency->GetBinContent(i) > 0) | |
699 | matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i); | |
700 | else | |
701 | matrixM[i-1][j-1] = 0; | |
702 | } | |
703 | } | |
704 | ||
705 | Float_t* vectorn = new Float_t[yBins]; | |
706 | for (Int_t j=1; j<=yBins; j++) | |
707 | vectorn[j-1] = fCurrentESD->GetBinContent(j); | |
708 | ||
709 | // first part of covariance matrix, depends on input distribution n(E) | |
710 | Float_t cov1[251][251]; | |
711 | ||
712 | Float_t nEvents = fCurrentESD->Integral(); // N | |
713 | ||
714 | xBins = 20; | |
715 | yBins = 20; | |
716 | ||
717 | for (Int_t k=0; k<xBins; k++) | |
718 | { | |
719 | printf("In Cov1: %d\n", k); | |
720 | for (Int_t l=0; l<yBins; l++) | |
721 | { | |
722 | cov1[k][l] = 0; | |
723 | ||
724 | // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N) | |
725 | for (Int_t j=0; j<yBins; j++) | |
726 | cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j] | |
727 | * (1.0 - vectorn[j] / nEvents); | |
728 | ||
729 | // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N | |
730 | for (Int_t i=0; i<yBins; i++) | |
731 | for (Int_t j=0; j<yBins; j++) | |
732 | { | |
733 | if (i == j) | |
734 | continue; | |
735 | cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i] | |
736 | * vectorn[j] / nEvents; | |
737 | } | |
738 | } | |
739 | } | |
740 | ||
741 | printf("Cov1 finished\n"); | |
742 | ||
743 | TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov"); | |
744 | cov->Reset(); | |
745 | ||
746 | for (Int_t i=1; i<=xBins; i++) | |
747 | for (Int_t j=1; j<=yBins; j++) | |
748 | cov->SetBinContent(i, j, cov1[i-1][j-1]); | |
749 | ||
750 | new TCanvas; | |
751 | cov->Draw("COLZ"); | |
752 | ||
753 | // second part of covariance matrix, depends on response matrix | |
754 | Float_t cov2[251][251]; | |
755 | ||
756 | // Cov[P(Er|Cu), P(Es|Cu)] term | |
757 | Float_t covTerm[100][100][100]; | |
758 | for (Int_t r=0; r<yBins; r++) | |
759 | for (Int_t u=0; u<xBins; u++) | |
760 | for (Int_t s=0; s<yBins; s++) | |
761 | { | |
762 | if (r == s) | |
763 | covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) | |
764 | * (1.0 - hResponse->GetBinContent(u+1, r+1)); | |
765 | else | |
766 | covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) | |
767 | * hResponse->GetBinContent(u+1, s+1); | |
768 | } | |
769 | ||
770 | for (Int_t k=0; k<xBins; k++) | |
771 | for (Int_t l=0; l<yBins; l++) | |
772 | { | |
773 | cov2[k][l] = 0; | |
774 | printf("In Cov2: %d %d\n", k, l); | |
775 | for (Int_t i=0; i<yBins; i++) | |
776 | for (Int_t j=0; j<yBins; j++) | |
777 | { | |
778 | //printf("In Cov2: %d %d %d %d\n", k, l, i, j); | |
779 | // calculate Cov(Mki, Mlj) = sum{ru},{su} ... | |
780 | Float_t tmpCov = 0; | |
781 | for (Int_t r=0; r<yBins; r++) | |
782 | for (Int_t u=0; u<xBins; u++) | |
783 | for (Int_t s=0; s<yBins; s++) | |
784 | { | |
785 | if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0 | |
786 | || hResponse->GetBinContent(u+1, i+1) == 0) | |
787 | continue; | |
788 | ||
789 | tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u) | |
790 | * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u) | |
791 | * covTerm[r][u][s]; | |
792 | } | |
793 | ||
794 | cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov; | |
795 | } | |
796 | } | |
797 | ||
798 | printf("Cov2 finished\n"); | |
799 | ||
800 | for (Int_t i=1; i<=xBins; i++) | |
801 | for (Int_t j=1; j<=yBins; j++) | |
802 | cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]); | |
803 | ||
804 | new TCanvas; | |
805 | cov->Draw("COLZ");*/ | |
806 | } | |
807 | ||
808 | //____________________________________________________________________ | |
809 | Double_t AliUnfolding::RegularizationPol0(TVectorD& params) | |
810 | { | |
811 | // homogenity term for minuit fitting | |
812 | // pure function of the parameters | |
813 | // prefers constant function (pol0) | |
814 | ||
815 | Double_t chi2 = 0; | |
816 | ||
817 | for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
818 | { | |
95e970ca | 819 | Double_t right = params[i] / (*fgBinWidths)(i); |
820 | Double_t left = params[i-1] / (*fgBinWidths)(i-1); | |
19442b86 | 821 | |
95e970ca | 822 | if (left != 0) |
19442b86 | 823 | { |
95e970ca | 824 | Double_t diff = (right - left); |
825 | chi2 += diff * diff / left / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2); | |
19442b86 | 826 | } |
827 | } | |
828 | ||
829 | return chi2 / 100.0; | |
830 | } | |
831 | ||
832 | //____________________________________________________________________ | |
833 | Double_t AliUnfolding::RegularizationPol1(TVectorD& params) | |
834 | { | |
835 | // homogenity term for minuit fitting | |
836 | // pure function of the parameters | |
837 | // prefers linear function (pol1) | |
838 | ||
839 | Double_t chi2 = 0; | |
840 | ||
841 | for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
842 | { | |
843 | if (params[i-1] == 0) | |
844 | continue; | |
845 | ||
95e970ca | 846 | Double_t right = params[i] / (*fgBinWidths)(i); |
847 | Double_t middle = params[i-1] / (*fgBinWidths)(i-1); | |
848 | Double_t left = params[i-2] / (*fgBinWidths)(i-2); | |
19442b86 | 849 | |
95e970ca | 850 | Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2); |
851 | Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 2); | |
19442b86 | 852 | |
95e970ca | 853 | //Double_t diff = (der1 - der2) / middle; |
854 | //chi2 += diff * diff; | |
855 | chi2 += (der1 - der2) * (der1 - der2) / middle * (*fgBinWidths)(i-1); | |
19442b86 | 856 | } |
857 | ||
858 | return chi2; | |
859 | } | |
860 | ||
861 | //____________________________________________________________________ | |
862 | Double_t AliUnfolding::RegularizationLog(TVectorD& params) | |
863 | { | |
864 | // homogenity term for minuit fitting | |
865 | // pure function of the parameters | |
866 | // prefers linear function (pol1) | |
867 | ||
868 | Double_t chi2 = 0; | |
869 | ||
870 | for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
871 | { | |
95e970ca | 872 | if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0) |
873 | continue; | |
19442b86 | 874 | |
95e970ca | 875 | Double_t right = log(params[i] / (*fgBinWidths)(i)); |
876 | Double_t middle = log(params[i-1] / (*fgBinWidths)(i-1)); | |
877 | Double_t left = log(params[i-2] / (*fgBinWidths)(i-2)); | |
878 | ||
879 | Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2); | |
880 | Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 2); | |
881 | ||
882 | //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2]; | |
19442b86 | 883 | |
95e970ca | 884 | //if (fgCallCount == 0) |
885 | // Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error); | |
886 | chi2 += (der1 - der2) * (der1 - der2);// / error; | |
19442b86 | 887 | } |
888 | ||
95e970ca | 889 | return chi2; |
19442b86 | 890 | } |
891 | ||
892 | //____________________________________________________________________ | |
893 | Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params) | |
894 | { | |
895 | // homogenity term for minuit fitting | |
896 | // pure function of the parameters | |
897 | // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments, | |
898 | // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp. | |
899 | ||
900 | Double_t chi2 = 0; | |
901 | ||
902 | for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
903 | { | |
904 | Double_t right = params[i]; | |
905 | Double_t middle = params[i-1]; | |
906 | Double_t left = params[i-2]; | |
907 | ||
908 | Double_t der1 = (right - middle); | |
909 | Double_t der2 = (middle - left); | |
910 | ||
911 | Double_t diff = (der1 - der2); | |
912 | ||
913 | chi2 += diff * diff; | |
914 | } | |
915 | ||
916 | return chi2 * 1e4; | |
917 | } | |
918 | ||
919 | //____________________________________________________________________ | |
920 | Double_t AliUnfolding::RegularizationEntropy(TVectorD& params) | |
921 | { | |
922 | // homogenity term for minuit fitting | |
923 | // pure function of the parameters | |
924 | // calculates entropy, from | |
925 | // The method of reduced cross-entropy (M. Schmelling 1993) | |
926 | ||
927 | Double_t paramSum = 0; | |
928 | ||
929 | for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i) | |
930 | paramSum += params[i]; | |
931 | ||
932 | Double_t chi2 = 0; | |
933 | for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i) | |
934 | { | |
95e970ca | 935 | Double_t tmp = params[i] / paramSum; |
936 | //Double_t tmp = params[i]; | |
19442b86 | 937 | if (tmp > 0 && (*fgEntropyAPriori)[i] > 0) |
938 | { | |
939 | chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]); | |
940 | } | |
95e970ca | 941 | else |
942 | chi2 += 100; | |
943 | } | |
944 | ||
945 | return -chi2; | |
946 | } | |
947 | ||
948 | //____________________________________________________________________ | |
949 | Double_t AliUnfolding::RegularizationRatio(TVectorD& params) | |
950 | { | |
951 | // homogenity term for minuit fitting | |
952 | // pure function of the parameters | |
953 | ||
954 | Double_t chi2 = 0; | |
955 | ||
956 | for (Int_t i=5+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
957 | { | |
958 | if (params[i-1] == 0 || params[i] == 0) | |
959 | continue; | |
960 | ||
961 | Double_t right = params[i] / (*fgBinWidths)(i); | |
962 | Double_t middle = params[i-1] / (*fgBinWidths)(i-1); | |
963 | Double_t left = params[i-2] / (*fgBinWidths)(i-2); | |
964 | Double_t left2 = params[i-3] / (*fgBinWidths)(i-3); | |
965 | Double_t left3 = params[i-4] / (*fgBinWidths)(i-4); | |
966 | Double_t left4 = params[i-5] / (*fgBinWidths)(i-5); | |
967 | ||
968 | //Double_t diff = left / middle - middle / right; | |
969 | //Double_t diff = 2 * left / middle - middle / right - left2 / left; | |
970 | Double_t diff = 4 * left2 / left - middle / right - left / middle - left3 / left2 - left4 / left3; | |
971 | ||
972 | chi2 += diff * diff;// / middle; | |
19442b86 | 973 | } |
974 | ||
95e970ca | 975 | return chi2; |
19442b86 | 976 | } |
977 | ||
978 | //____________________________________________________________________ | |
979 | void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t) | |
980 | { | |
981 | // | |
982 | // fit function for minuit | |
983 | // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix | |
984 | // | |
95e970ca | 985 | |
986 | // TODO use static members for the variables here to speed up processing (no construction/deconstruction) | |
19442b86 | 987 | |
988 | // d | |
95e970ca | 989 | TVectorD paramsVector(fgMaxParams); |
19442b86 | 990 | for (Int_t i=0; i<fgMaxParams; ++i) |
991 | paramsVector[i] = params[i] * params[i]; | |
992 | ||
993 | // calculate penalty factor | |
994 | Double_t penaltyVal = 0; | |
995 | switch (fgRegularizationType) | |
996 | { | |
997 | case kNone: break; | |
998 | case kPol0: penaltyVal = RegularizationPol0(paramsVector); break; | |
999 | case kPol1: penaltyVal = RegularizationPol1(paramsVector); break; | |
1000 | case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break; | |
1001 | case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break; | |
1002 | case kLog: penaltyVal = RegularizationLog(paramsVector); break; | |
95e970ca | 1003 | case kRatio: penaltyVal = RegularizationRatio(paramsVector); break; |
19442b86 | 1004 | } |
1005 | ||
1006 | //if (penaltyVal > 1e10) | |
1007 | // paramsVector2.Print(); | |
1008 | ||
1009 | penaltyVal *= fgRegularizationWeight; | |
1010 | ||
1011 | // Ad | |
1012 | TVectorD measGuessVector(fgMaxInput); | |
1013 | measGuessVector = (*fgCorrelationMatrix) * paramsVector; | |
1014 | ||
1015 | // Ad - m | |
1016 | measGuessVector -= (*fgCurrentESDVector); | |
95e970ca | 1017 | |
1018 | #if 0 | |
1019 | // new error calcuation using error on the guess | |
1020 | ||
1021 | // error from guess | |
1022 | TVectorD errorGuessVector(fgMaxInput); | |
1023 | //errorGuessVector = (*fgCorrelationMatrixSquared) * paramsVector; | |
1024 | errorGuessVector = (*fgCorrelationMatrix) * paramsVector; | |
19442b86 | 1025 | |
95e970ca | 1026 | Double_t chi2FromFit = 0; |
1027 | for (Int_t i=0; i<fgMaxInput; ++i) | |
1028 | { | |
1029 | Float_t error = 1; | |
1030 | if (errorGuessVector(i) > 0) | |
1031 | error = errorGuessVector(i); | |
1032 | chi2FromFit += measGuessVector(i) * measGuessVector(i) / error; | |
1033 | } | |
19442b86 | 1034 | //measGuessVector.Print(); |
1035 | ||
95e970ca | 1036 | #else |
1037 | // old error calcuation using the error on the measured | |
19442b86 | 1038 | TVectorD copy(measGuessVector); |
1039 | ||
1040 | // (Ad - m) W | |
1041 | // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used | |
1042 | // normal way is like this: | |
1043 | // measGuessVector *= (*fgCorrelationCovarianceMatrix); | |
1044 | // optimized way like this: | |
1045 | for (Int_t i=0; i<fgMaxInput; ++i) | |
1046 | measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i); | |
1047 | ||
1048 | //measGuessVector.Print(); | |
1049 | ||
95e970ca | 1050 | if (fgSkipBin0InChi2) |
1051 | measGuessVector[0] = 0; | |
1052 | ||
19442b86 | 1053 | // (Ad - m) W (Ad - m) |
1054 | // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very | |
95e970ca | 1055 | // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see UnfoldWithMinuit) |
19442b86 | 1056 | Double_t chi2FromFit = measGuessVector * copy * 1e6; |
95e970ca | 1057 | #endif |
19442b86 | 1058 | |
95e970ca | 1059 | Double_t notFoundEventsConstraint = 0; |
1060 | Double_t currentNotFoundEvents = 0; | |
1061 | Double_t errorNotFoundEvents = 0; | |
1062 | ||
1063 | if (fgNotFoundEvents > 0) | |
1064 | { | |
1065 | for (Int_t i=0; i<fgMaxParams; ++i) | |
1066 | { | |
1067 | Double_t eff = (1.0 / (*fgEfficiency)(i) - 1); | |
1068 | if (i == 0) | |
1069 | eff = (1.0 / params[fgMaxParams] - 1); | |
1070 | currentNotFoundEvents += eff * paramsVector(i); | |
1071 | errorNotFoundEvents += eff * eff * paramsVector(i); // error due to guess (paramsVector) | |
1072 | if (i <= 3) | |
1073 | errorNotFoundEvents += (eff * 0.03) * (eff * 0.03) * paramsVector(i) * paramsVector(i); // error on eff | |
1074 | //if ((fgCallCount % 10000) == 0) | |
1075 | // Printf("%d %f %f %f", i, (*fgEfficiency)(i), paramsVector(i), currentNotFoundEvents); | |
1076 | } | |
1077 | errorNotFoundEvents += fgNotFoundEvents; | |
1078 | // TODO add error on background, background estimate | |
1079 | ||
1080 | notFoundEventsConstraint = (currentNotFoundEvents - fgNotFoundEvents) * (currentNotFoundEvents - fgNotFoundEvents) / errorNotFoundEvents; | |
1081 | } | |
1082 | ||
1083 | Float_t avg = 0; | |
1084 | Float_t sum = 0; | |
1085 | Float_t currentMult = 0; | |
1086 | // try to match dn/deta | |
1087 | for (Int_t i=0; i<fgMaxParams; ++i) | |
1088 | { | |
1089 | avg += paramsVector(i) * currentMult; | |
1090 | sum += paramsVector(i); | |
1091 | currentMult += (*fgBinWidths)(i); | |
1092 | } | |
1093 | avg /= sum; | |
1094 | Float_t chi2Avg = 0; //(avg - 3.73) * (avg - 3.73) * 100; | |
19442b86 | 1095 | |
95e970ca | 1096 | chi2 = chi2FromFit + penaltyVal + notFoundEventsConstraint + chi2Avg; |
1097 | ||
1098 | if ((fgCallCount++ % 1000) == 0) | |
1099 | { | |
1100 | Printf("AliUnfolding::Chi2Function: Iteration %d (ev %d %d +- %f) (%f) (%f): %f %f %f %f --> %f", fgCallCount-1, (Int_t) fgNotFoundEvents, (Int_t) currentNotFoundEvents, TMath::Sqrt(errorNotFoundEvents), params[fgMaxParams], avg, chi2FromFit, penaltyVal, notFoundEventsConstraint, chi2Avg, chi2); | |
1101 | //for (Int_t i=0; i<fgMaxInput; ++i) | |
1102 | // Printf("%d: %f", i, measGuessVector(i) * copy(i) * 1e6); | |
1103 | } | |
19442b86 | 1104 | } |
1105 | ||
1106 | //____________________________________________________________________ | |
1107 | void AliUnfolding::DrawGuess(Double_t *params) | |
1108 | { | |
1109 | // | |
1110 | // draws residuals of solution suggested by params and effect of regularization | |
1111 | // | |
1112 | ||
1113 | // d | |
95e970ca | 1114 | TVectorD paramsVector(fgMaxParams); |
19442b86 | 1115 | for (Int_t i=0; i<fgMaxParams; ++i) |
1116 | paramsVector[i] = params[i] * params[i]; | |
1117 | ||
1118 | // Ad | |
1119 | TVectorD measGuessVector(fgMaxInput); | |
1120 | measGuessVector = (*fgCorrelationMatrix) * paramsVector; | |
1121 | ||
1122 | // Ad - m | |
1123 | measGuessVector -= (*fgCurrentESDVector); | |
1124 | ||
1125 | TVectorD copy(measGuessVector); | |
1126 | //copy.Print(); | |
1127 | ||
1128 | // (Ad - m) W | |
1129 | // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used | |
1130 | // normal way is like this: | |
1131 | // measGuessVector *= (*fgCorrelationCovarianceMatrix); | |
1132 | // optimized way like this: | |
1133 | for (Int_t i=0; i<fgMaxInput; ++i) | |
1134 | measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i); | |
1135 | ||
1136 | // (Ad - m) W (Ad - m) | |
1137 | // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very | |
1138 | // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit) | |
1139 | //Double_t chi2FromFit = measGuessVector * copy * 1e6; | |
1140 | ||
1141 | // draw residuals | |
1142 | TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5); | |
1143 | for (Int_t i=0; i<fgMaxInput; ++i) | |
1144 | residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6); | |
1145 | new TCanvas; residuals->DrawCopy(); gPad->SetLogy(); | |
1146 | ||
1147 | // draw penalty | |
95e970ca | 1148 | TH1* penalty = GetPenaltyPlot(params); |
1149 | ||
1150 | new TCanvas; penalty->DrawCopy(); gPad->SetLogy(); | |
1151 | } | |
19442b86 | 1152 | |
95e970ca | 1153 | //____________________________________________________________________ |
1154 | TH1* AliUnfolding::GetPenaltyPlot(TH1* corrected) | |
1155 | { | |
1156 | // draws the penalty factors as function of multiplicity of the current selected regularization | |
19442b86 | 1157 | |
95e970ca | 1158 | Double_t* params = new Double_t[fgMaxParams]; |
1159 | for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++) | |
1160 | params[i] = corrected->GetBinContent(i+1); | |
1161 | ||
1162 | TH1* penalty = GetPenaltyPlot(params); | |
1163 | ||
1164 | delete[] params; | |
1165 | ||
1166 | return penalty; | |
1167 | } | |
1168 | ||
1169 | //____________________________________________________________________ | |
1170 | TH1* AliUnfolding::GetPenaltyPlot(Double_t* params) | |
1171 | { | |
1172 | // draws the penalty factors as function of multiplicity of the current selected regularization | |
1173 | ||
1174 | TH1* penalty = new TH1F("penalty", ";unfolded multiplicity;penalty factor", fgMaxParams, -0.5, fgMaxParams - 0.5); | |
1175 | ||
1176 | for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
19442b86 | 1177 | { |
95e970ca | 1178 | Double_t diff = 0; |
1179 | if (fgRegularizationType == kPol0) | |
1180 | { | |
1181 | Double_t right = params[i] / (*fgBinWidths)(i); | |
1182 | Double_t left = params[i-1] / (*fgBinWidths)(i-1); | |
19442b86 | 1183 | |
95e970ca | 1184 | if (left != 0) |
1185 | { | |
1186 | Double_t diffTmp = (right - left); | |
1187 | diff = diffTmp * diffTmp / left / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2) / 100; | |
1188 | } | |
1189 | } | |
1190 | if (fgRegularizationType == kPol1 && i > 1+fgSkipBinsBegin) | |
1191 | { | |
1192 | if (params[i-1] == 0) | |
1193 | continue; | |
19442b86 | 1194 | |
95e970ca | 1195 | Double_t right = params[i] / (*fgBinWidths)(i); |
1196 | Double_t middle = params[i-1] / (*fgBinWidths)(i-1); | |
1197 | Double_t left = params[i-2] / (*fgBinWidths)(i-2); | |
19442b86 | 1198 | |
95e970ca | 1199 | Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2); |
1200 | Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 2); | |
1201 | ||
1202 | diff = (der1 - der2) * (der1 - der2) / middle; | |
1203 | } | |
1204 | if (fgRegularizationType == kLog && i > 1+fgSkipBinsBegin) | |
1205 | { | |
1206 | if (params[i-1] == 0) | |
1207 | continue; | |
1208 | ||
1209 | Double_t right = log(params[i]); | |
1210 | Double_t middle = log(params[i-1]); | |
1211 | Double_t left = log(params[i-2]); | |
1212 | ||
1213 | Double_t der1 = (right - middle); | |
1214 | Double_t der2 = (middle - left); | |
1215 | ||
1216 | //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2]; | |
1217 | //Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error); | |
1218 | ||
1219 | diff = (der1 - der2) * (der1 - der2);// / error; | |
1220 | } | |
19442b86 | 1221 | |
95e970ca | 1222 | penalty->SetBinContent(i, diff); |
19442b86 | 1223 | |
1224 | //Printf("%d %f %f %f %f", i-1, left, middle, right, diff); | |
1225 | } | |
95e970ca | 1226 | |
1227 | return penalty; | |
19442b86 | 1228 | } |
1229 | ||
1230 | //____________________________________________________________________ | |
1231 | void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3) | |
1232 | { | |
1233 | // | |
1234 | // fit function for minuit | |
1235 | // uses the TF1 stored in fgFitFunction | |
1236 | // | |
1237 | ||
1238 | for (Int_t i=0; i<fgFitFunction->GetNpar(); i++) | |
1239 | fgFitFunction->SetParameter(i, params[i]); | |
1240 | ||
1241 | Double_t* params2 = new Double_t[fgMaxParams]; | |
1242 | ||
1243 | for (Int_t i=0; i<fgMaxParams; ++i) | |
1244 | params2[i] = fgFitFunction->Eval(i); | |
1245 | ||
1246 | Chi2Function(unused1, unused2, chi2, params2, unused3); | |
1247 | ||
1248 | delete[] params2; | |
1249 | ||
1250 | if (fgDebug) | |
1251 | Printf("%f", chi2); | |
1252 | } | |
1253 | ||
1254 | //____________________________________________________________________ | |
1255 | Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult) | |
1256 | { | |
1257 | // | |
1258 | // correct spectrum using minuit chi2 method applying a functional fit | |
1259 | // | |
1260 | ||
1261 | if (!fgFitFunction) | |
1262 | { | |
1263 | Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting."); | |
1264 | return -1; | |
1265 | } | |
1266 | ||
1267 | SetChi2Regularization(kNone, 0); | |
1268 | ||
95e970ca | 1269 | SetStaticVariables(correlation, measured, efficiency); |
19442b86 | 1270 | |
1271 | // Initialize TMinuit via generic fitter interface | |
1272 | TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar()); | |
1273 | ||
1274 | minuit->SetFCN(TF1Function); | |
1275 | for (Int_t i=0; i<fgFitFunction->GetNpar(); i++) | |
1276 | { | |
1277 | Double_t lower, upper; | |
1278 | fgFitFunction->GetParLimits(i, lower, upper); | |
1279 | minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper); | |
1280 | } | |
1281 | ||
1282 | Double_t arglist[100]; | |
1283 | arglist[0] = 0; | |
1284 | minuit->ExecuteCommand("SET PRINT", arglist, 1); | |
1285 | minuit->ExecuteCommand("MIGRAD", arglist, 0); | |
1286 | //minuit->ExecuteCommand("MINOS", arglist, 0); | |
1287 | ||
1288 | for (Int_t i=0; i<fgFitFunction->GetNpar(); i++) | |
1289 | fgFitFunction->SetParameter(i, minuit->GetParameter(i)); | |
1290 | ||
1291 | for (Int_t i=0; i<fgMaxParams; ++i) | |
1292 | { | |
1293 | Double_t value = fgFitFunction->Eval(i); | |
1294 | if (fgDebug) | |
1295 | Printf("%d : %f", i, value); | |
1296 | ||
1297 | if (efficiency) | |
1298 | { | |
1299 | if (efficiency->GetBinContent(i+1) > 0) | |
1300 | { | |
1301 | value /= efficiency->GetBinContent(i+1); | |
1302 | } | |
1303 | else | |
1304 | value = 0; | |
1305 | } | |
1306 | aResult->SetBinContent(i+1, value); | |
1307 | aResult->SetBinError(i+1, 0); | |
1308 | } | |
1309 | ||
1310 | return 0; | |
1311 | } | |
1312 | ||
1313 | //____________________________________________________________________ | |
1314 | void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured) | |
1315 | { | |
1316 | // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins | |
1317 | // The same limit is applied to the correlation | |
1318 | ||
1319 | Int_t lastBin = 0; | |
1320 | for (Int_t i=1; i<=measured->GetNbinsX(); ++i) | |
1321 | { | |
1322 | if (measured->GetBinContent(i) <= fgOverflowBinLimit) | |
1323 | { | |
1324 | lastBin = i; | |
1325 | break; | |
1326 | } | |
1327 | } | |
1328 | ||
1329 | if (lastBin == 0) | |
1330 | { | |
1331 | Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit); | |
1332 | return; | |
1333 | } | |
1334 | ||
1335 | Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin); | |
1336 | ||
1337 | TCanvas* canvas = 0; | |
1338 | ||
1339 | if (fgDebug) | |
1340 | { | |
1341 | canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800); | |
1342 | canvas->Divide(2, 2); | |
1343 | ||
1344 | canvas->cd(1); | |
1345 | measured->SetStats(kFALSE); | |
1346 | measured->DrawCopy(); | |
1347 | gPad->SetLogy(); | |
1348 | ||
1349 | canvas->cd(2); | |
1350 | correlation->SetStats(kFALSE); | |
1351 | correlation->DrawCopy("COLZ"); | |
1352 | } | |
1353 | ||
1354 | measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX())); | |
1355 | for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i) | |
1356 | { | |
1357 | measured->SetBinContent(i, 0); | |
1358 | measured->SetBinError(i, 0); | |
1359 | } | |
1360 | // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? | |
1361 | measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin))); | |
1362 | ||
1363 | Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin)); | |
1364 | ||
1365 | for (Int_t i=1; i<=correlation->GetNbinsX(); ++i) | |
1366 | { | |
1367 | correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY())); | |
1368 | // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? | |
1369 | correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin))); | |
1370 | ||
1371 | for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j) | |
1372 | { | |
1373 | correlation->SetBinContent(i, j, 0); | |
1374 | correlation->SetBinError(i, j, 0); | |
1375 | } | |
1376 | } | |
1377 | ||
1378 | Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!"); | |
1379 | ||
1380 | if (canvas) | |
1381 | { | |
1382 | canvas->cd(3); | |
1383 | measured->DrawCopy(); | |
1384 | gPad->SetLogy(); | |
1385 | ||
1386 | canvas->cd(4); | |
1387 | correlation->DrawCopy("COLZ"); | |
1388 | } | |
1389 | } | |
95e970ca | 1390 | |
1391 | Int_t AliUnfolding::UnfoldGetBias(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result) | |
1392 | { | |
1393 | // unfolds and assigns bias as errors with Eq. 19 of Cowan, "a survey of unfolding methods for particle physics" | |
1394 | // b_i = sum_j dmu_i/dn_j (nu_j - n_j) with nu_j as folded guess, n_j as data | |
1395 | // dmu_i/dn_j is found numerically by changing the bin content and re-unfolding | |
1396 | // | |
1397 | // return code: 0 (success) -1 (error: from Unfold(...) ) | |
1398 | ||
1399 | if (Unfold(correlation, efficiency, measured, initialConditions, result) != 0) | |
1400 | return -1; | |
1401 | ||
1402 | TMatrixD matrix(fgMaxInput, fgMaxParams); | |
1403 | ||
1404 | TH1* newResult[4]; | |
1405 | for (Int_t loop=0; loop<4; loop++) | |
1406 | newResult[loop] = (TH1*) result->Clone(Form("newresult_%d", loop)); | |
1407 | ||
1408 | // change bin-by-bin and built matrix of effects | |
1409 | for (Int_t m=0; m<fgMaxInput; m++) | |
1410 | { | |
1411 | if (measured->GetBinContent(m+1) < 1) | |
1412 | continue; | |
1413 | ||
1414 | for (Int_t loop=0; loop<4; loop++) | |
1415 | { | |
1416 | //Printf("%d %d", i, loop); | |
1417 | Float_t factor = 1; | |
1418 | switch (loop) | |
1419 | { | |
1420 | case 0: factor = 0.5; break; | |
1421 | case 1: factor = -0.5; break; | |
1422 | case 2: factor = 1; break; | |
1423 | case 3: factor = -1; break; | |
1424 | default: return -1; | |
1425 | } | |
1426 | ||
1427 | TH1* measuredClone = (TH1*) measured->Clone("measuredClone"); | |
1428 | ||
1429 | measuredClone->SetBinContent(m+1, measured->GetBinContent(m+1) + factor * TMath::Sqrt(measured->GetBinContent(m+1))); | |
1430 | //new TCanvas; measuredClone->Draw("COLZ"); | |
1431 | ||
1432 | newResult[loop]->Reset(); | |
1433 | Unfold(correlation, efficiency, measuredClone, measuredClone, newResult[loop]); | |
1434 | // WARNING if we leave here, then nothing is calculated | |
1435 | //return -1; | |
1436 | ||
1437 | delete measuredClone; | |
1438 | } | |
1439 | ||
1440 | for (Int_t t=0; t<fgMaxParams; t++) | |
1441 | { | |
1442 | // one value | |
1443 | //matrix(i, j) = (result->GetBinContent(j+1) - newResult->GetBinContent(j+1)) / TMath::Sqrt(changedHist->GetBinContent(1, i+1)); | |
1444 | ||
1445 | // four values from the derivate (procedure taken from ROOT -- suggestion by Ruben) | |
1446 | // = 1/2D * [ 8 (f(D/2) - f(-D/2)) - (f(D)-f(-D)) ]/3 | |
1447 | ||
1448 | /* | |
1449 | // check formula | |
1450 | measured->SetBinContent(1, m+1, 100); | |
1451 | newResult[0]->SetBinContent(t+1, 5 * 0.5 + 10); | |
1452 | newResult[1]->SetBinContent(t+1, 5 * -0.5 + 10); | |
1453 | newResult[2]->SetBinContent(t+1, 5 * 1 + 10); | |
1454 | newResult[3]->SetBinContent(t+1, 5 * -1 + 10); | |
1455 | */ | |
1456 | ||
1457 | matrix(m, t) = 0.5 / TMath::Sqrt(measured->GetBinContent(m+1)) * | |
1458 | ( 8. * (newResult[0]->GetBinContent(t+1) - newResult[1]->GetBinContent(t+1)) - | |
1459 | (newResult[2]->GetBinContent(t+1) - newResult[3]->GetBinContent(t+1)) | |
1460 | ) / 3; | |
1461 | } | |
1462 | } | |
1463 | ||
1464 | for (Int_t loop=0; loop<4; loop++) | |
1465 | delete newResult[loop]; | |
1466 | ||
1467 | //matrix.Print(); | |
1468 | //TH2* matrixHist = new TH2D(matrix); new TCanvas; matrixHist->Draw("BOX"); | |
1469 | ||
1470 | // ...calculate folded guess | |
1471 | TH1* convoluted = (TH1*) measured->Clone("convoluted"); | |
1472 | convoluted->Reset(); | |
1473 | for (Int_t m=0; m<fgMaxInput; m++) | |
1474 | { | |
1475 | Float_t value = 0; | |
1476 | for (Int_t t = 0; t<fgMaxParams; t++) | |
1477 | { | |
1478 | Float_t tmp = correlation->GetBinContent(t+1, m+1) * result->GetBinContent(t+1); | |
1479 | if (efficiency) | |
1480 | tmp *= efficiency->GetBinContent(t+1); | |
1481 | value += tmp; | |
1482 | } | |
1483 | convoluted->SetBinContent(m+1, value); | |
1484 | } | |
1485 | ||
1486 | // rescale | |
1487 | convoluted->Scale(measured->Integral() / convoluted->Integral()); | |
1488 | ||
1489 | //new TCanvas; convoluted->Draw(); measured->Draw("SAME"); measured->SetLineColor(2); | |
1490 | //return; | |
1491 | ||
1492 | // difference | |
1493 | convoluted->Add(measured, -1); | |
1494 | ||
1495 | // ...and the bias | |
1496 | for (Int_t t = 0; t<fgMaxParams; t++) | |
1497 | { | |
1498 | Double_t error = 0; | |
1499 | for (Int_t m=0; m<fgMaxInput; m++) | |
1500 | error += matrix(m, t) * convoluted->GetBinContent(m+1); | |
1501 | result->SetBinError(t+1, error); | |
1502 | } | |
1503 | ||
1504 | //new TCanvas; result->Draw(); gPad->SetLogy(); | |
1505 | ||
1506 | return 0; | |
1507 | } |