]>
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> | |
9e065ad2 | 30 | #include <TExec.h> |
31 | #include "Riostream.h" | |
32 | #include "TROOT.h" | |
33 | ||
34 | using namespace std; //required for resolving the 'cout' symbol | |
19442b86 | 35 | |
36 | TMatrixD* AliUnfolding::fgCorrelationMatrix = 0; | |
95e970ca | 37 | TMatrixD* AliUnfolding::fgCorrelationMatrixSquared = 0; |
19442b86 | 38 | TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0; |
39 | TVectorD* AliUnfolding::fgCurrentESDVector = 0; | |
40 | TVectorD* AliUnfolding::fgEntropyAPriori = 0; | |
95e970ca | 41 | TVectorD* AliUnfolding::fgEfficiency = 0; |
9e065ad2 | 42 | |
43 | TAxis* AliUnfolding::fgUnfoldedAxis = 0; | |
44 | TAxis* AliUnfolding::fgMeasuredAxis = 0; | |
19442b86 | 45 | |
46 | TF1* AliUnfolding::fgFitFunction = 0; | |
47 | ||
48 | AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid; | |
49 | Int_t AliUnfolding::fgMaxInput = -1; // bins in measured histogram | |
50 | Int_t AliUnfolding::fgMaxParams = -1; // bins in unfolded histogram = number of fit params | |
51 | Float_t AliUnfolding::fgOverflowBinLimit = -1; | |
52 | ||
53 | AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1; | |
54 | Float_t AliUnfolding::fgRegularizationWeight = 10000; | |
55 | Int_t AliUnfolding::fgSkipBinsBegin = 0; | |
56 | Float_t AliUnfolding::fgMinuitStepSize = 0.1; // (usually not needed to be changed) step size in minimization | |
9e065ad2 | 57 | Float_t AliUnfolding::fgMinuitPrecision = 1e-6; // minuit precision |
736125ad | 58 | Int_t AliUnfolding::fgMinuitMaxIterations = 1000000; // minuit maximum number of iterations |
9f070ef5 | 59 | Double_t AliUnfolding::fgMinuitStrategy = 1.; // minuit strategy |
95e970ca | 60 | Bool_t AliUnfolding::fgMinimumInitialValue = kFALSE; // set all initial values at least to the smallest value among the initial values |
61 | Float_t AliUnfolding::fgMinimumInitialValueFix = -1; | |
9e065ad2 | 62 | Bool_t AliUnfolding::fgNormalizeInput = kFALSE; // normalize input spectrum |
95e970ca | 63 | Float_t AliUnfolding::fgNotFoundEvents = 0; |
64 | Bool_t AliUnfolding::fgSkipBin0InChi2 = kFALSE; | |
19442b86 | 65 | |
66 | Float_t AliUnfolding::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing) | |
651e2088 | 67 | Int_t AliUnfolding::fgBayesianIterations = 10; // number of iterations in Bayesian method |
19442b86 | 68 | |
69 | Bool_t AliUnfolding::fgDebug = kFALSE; | |
70 | ||
95e970ca | 71 | Int_t AliUnfolding::fgCallCount = 0; |
72 | ||
9e065ad2 | 73 | Int_t AliUnfolding::fgPowern = 5; |
74 | ||
75 | Double_t AliUnfolding::fChi2FromFit = 0.; | |
76 | Double_t AliUnfolding::fPenaltyVal = 0.; | |
77 | Double_t AliUnfolding::fAvgResidual = 0.; | |
78 | ||
79 | Int_t AliUnfolding::fgPrintChi2Details = 0; | |
80 | ||
81 | TCanvas *AliUnfolding::fgCanvas = 0; | |
82 | TH1 *AliUnfolding::fghUnfolded = 0; | |
83 | TH2 *AliUnfolding::fghCorrelation = 0; | |
84 | TH1 *AliUnfolding::fghEfficiency = 0; | |
85 | TH1 *AliUnfolding::fghMeasured = 0; | |
86 | ||
19442b86 | 87 | ClassImp(AliUnfolding) |
88 | ||
89 | //____________________________________________________________________ | |
90 | void AliUnfolding::SetUnfoldingMethod(MethodType methodType) | |
91 | { | |
92 | // set unfolding method | |
93 | fgMethodType = methodType; | |
94 | ||
95 | const char* name = 0; | |
96 | switch (methodType) | |
97 | { | |
98 | case kInvalid: name = "INVALID"; break; | |
99 | case kChi2Minimization: name = "Chi2 Minimization"; break; | |
100 | case kBayesian: name = "Bayesian unfolding"; break; | |
101 | case kFunction: name = "Functional fit"; break; | |
102 | } | |
103 | Printf("AliUnfolding::SetUnfoldingMethod: %s enabled.", name); | |
104 | } | |
105 | ||
106 | //____________________________________________________________________ | |
107 | void AliUnfolding::SetCreateOverflowBin(Float_t overflowBinLimit) | |
108 | { | |
109 | // enable the creation of a overflow bin that includes all statistics below the given limit | |
110 | ||
111 | fgOverflowBinLimit = overflowBinLimit; | |
112 | ||
113 | Printf("AliUnfolding::SetCreateOverflowBin: overflow bin limit set to %f", overflowBinLimit); | |
114 | } | |
115 | ||
116 | //____________________________________________________________________ | |
117 | void AliUnfolding::SetSkipBinsBegin(Int_t nBins) | |
118 | { | |
119 | // set number of skipped bins in regularization | |
120 | ||
121 | fgSkipBinsBegin = nBins; | |
122 | ||
123 | Printf("AliUnfolding::SetSkipBinsBegin: skipping %d bins at the beginning of the spectrum in the regularization.", fgSkipBinsBegin); | |
124 | } | |
125 | ||
126 | //____________________________________________________________________ | |
127 | void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded) | |
128 | { | |
129 | // set number of bins in the input (measured) distribution and in the unfolded distribution | |
130 | fgMaxInput = nMeasured; | |
131 | fgMaxParams = nUnfolded; | |
132 | ||
95e970ca | 133 | if (fgCorrelationMatrix) |
134 | { | |
135 | delete fgCorrelationMatrix; | |
136 | fgCorrelationMatrix = 0; | |
137 | } | |
138 | if (fgCorrelationMatrixSquared) | |
139 | { | |
140 | fgCorrelationMatrixSquared = 0; | |
141 | delete fgCorrelationMatrixSquared; | |
142 | } | |
143 | if (fgCorrelationCovarianceMatrix) | |
144 | { | |
145 | delete fgCorrelationCovarianceMatrix; | |
146 | fgCorrelationCovarianceMatrix = 0; | |
147 | } | |
148 | if (fgCurrentESDVector) | |
149 | { | |
150 | delete fgCurrentESDVector; | |
151 | fgCurrentESDVector = 0; | |
152 | } | |
153 | if (fgEntropyAPriori) | |
154 | { | |
155 | delete fgEntropyAPriori; | |
156 | fgEntropyAPriori = 0; | |
157 | } | |
158 | if (fgEfficiency) | |
159 | { | |
160 | delete fgEfficiency; | |
161 | fgEfficiency = 0; | |
162 | } | |
9e065ad2 | 163 | if (fgUnfoldedAxis) |
164 | { | |
165 | delete fgUnfoldedAxis; | |
166 | fgUnfoldedAxis = 0; | |
167 | } | |
168 | if (fgMeasuredAxis) | |
95e970ca | 169 | { |
9e065ad2 | 170 | delete fgMeasuredAxis; |
171 | fgMeasuredAxis = 0; | |
95e970ca | 172 | } |
173 | ||
19442b86 | 174 | Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded); |
175 | } | |
176 | ||
177 | //____________________________________________________________________ | |
178 | void AliUnfolding::SetChi2Regularization(RegularizationType type, Float_t weight) | |
179 | { | |
180 | // | |
181 | // sets the parameters for chi2 minimization | |
182 | // | |
183 | ||
184 | fgRegularizationType = type; | |
185 | fgRegularizationWeight = weight; | |
186 | ||
187 | Printf("AliUnfolding::SetChi2Regularization --> Regularization set to %d with weight %f", (Int_t) type, weight); | |
188 | } | |
189 | ||
190 | //____________________________________________________________________ | |
191 | void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations) | |
192 | { | |
193 | // | |
194 | // sets the parameters for Bayesian unfolding | |
195 | // | |
196 | ||
197 | fgBayesianSmoothing = smoothing; | |
198 | fgBayesianIterations = nIterations; | |
199 | ||
736125ad | 200 | Printf("AliUnfolding::SetBayesianParameters --> Parameters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing); |
19442b86 | 201 | } |
202 | ||
203 | //____________________________________________________________________ | |
204 | void AliUnfolding::SetFunction(TF1* function) | |
205 | { | |
206 | // set function for unfolding with a fit function | |
207 | ||
208 | fgFitFunction = function; | |
209 | ||
210 | Printf("AliUnfolding::SetFunction: Set fit function with %d parameters.", function->GetNpar()); | |
211 | } | |
212 | ||
213 | //____________________________________________________________________ | |
214 | Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check) | |
215 | { | |
216 | // unfolds with unfolding method fgMethodType | |
217 | // | |
218 | // parameters: | |
219 | // correlation: response matrix as measured vs. generated | |
220 | // 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. | |
221 | // measured: the measured spectrum | |
222 | // initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions. | |
223 | // result: target for the unfolded result | |
224 | // check: depends on the unfolding method, see comments in specific functions | |
95e970ca | 225 | // |
226 | // return code: see UnfoldWithMinuit/UnfoldWithBayesian/UnfoldWithFunction | |
19442b86 | 227 | |
228 | if (fgMaxInput == -1) | |
229 | { | |
230 | Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution"); | |
231 | fgMaxInput = measured->GetNbinsX(); | |
232 | } | |
233 | if (fgMaxParams == -1) | |
234 | { | |
235 | Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution"); | |
236 | fgMaxParams = measured->GetNbinsX(); | |
237 | } | |
238 | ||
239 | if (fgOverflowBinLimit > 0) | |
240 | CreateOverflowBin(correlation, measured); | |
241 | ||
242 | switch (fgMethodType) | |
243 | { | |
244 | case kInvalid: | |
245 | { | |
246 | Printf("AliUnfolding::Unfold: ERROR: Unfolding method not set. Use SetUnfoldingMethod. Exiting..."); | |
247 | return -1; | |
248 | } | |
249 | case kChi2Minimization: | |
250 | return UnfoldWithMinuit(correlation, efficiency, measured, initialConditions, result, check); | |
251 | case kBayesian: | |
252 | return UnfoldWithBayesian(correlation, efficiency, measured, initialConditions, result); | |
253 | case kFunction: | |
254 | return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result); | |
255 | } | |
9e065ad2 | 256 | |
257 | ||
258 | ||
19442b86 | 259 | return -1; |
260 | } | |
261 | ||
262 | //____________________________________________________________________ | |
95e970ca | 263 | void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured, TH1* efficiency) |
19442b86 | 264 | { |
265 | // fill static variables needed for minuit fit | |
266 | ||
267 | if (!fgCorrelationMatrix) | |
268 | fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams); | |
95e970ca | 269 | if (!fgCorrelationMatrixSquared) |
270 | fgCorrelationMatrixSquared = new TMatrixD(fgMaxInput, fgMaxParams); | |
19442b86 | 271 | if (!fgCorrelationCovarianceMatrix) |
272 | fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput); | |
273 | if (!fgCurrentESDVector) | |
274 | fgCurrentESDVector = new TVectorD(fgMaxInput); | |
275 | if (!fgEntropyAPriori) | |
276 | fgEntropyAPriori = new TVectorD(fgMaxParams); | |
95e970ca | 277 | if (!fgEfficiency) |
278 | fgEfficiency = new TVectorD(fgMaxParams); | |
736125ad | 279 | if (fgUnfoldedAxis) |
280 | { | |
9e065ad2 | 281 | delete fgUnfoldedAxis; |
736125ad | 282 | fgUnfoldedAxis = 0; |
283 | } | |
284 | if (!fgUnfoldedAxis) | |
285 | fgUnfoldedAxis = new TAxis(*(correlation->GetXaxis())); | |
286 | if (fgMeasuredAxis) | |
287 | { | |
9e065ad2 | 288 | delete fgMeasuredAxis; |
736125ad | 289 | fgMeasuredAxis = 0; |
290 | } | |
291 | if (!fgMeasuredAxis) | |
292 | fgMeasuredAxis = new TAxis(*(correlation->GetYaxis())); | |
9e065ad2 | 293 | |
19442b86 | 294 | fgCorrelationMatrix->Zero(); |
295 | fgCorrelationCovarianceMatrix->Zero(); | |
296 | fgCurrentESDVector->Zero(); | |
297 | fgEntropyAPriori->Zero(); | |
298 | ||
299 | // normalize correction for given nPart | |
300 | for (Int_t i=1; i<=correlation->GetNbinsX(); ++i) | |
301 | { | |
302 | Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY()); | |
303 | if (sum <= 0) | |
304 | continue; | |
305 | Float_t maxValue = 0; | |
306 | Int_t maxBin = -1; | |
307 | for (Int_t j=1; j<=correlation->GetNbinsY(); ++j) | |
308 | { | |
309 | // find most probably value | |
310 | if (maxValue < correlation->GetBinContent(i, j)) | |
311 | { | |
312 | maxValue = correlation->GetBinContent(i, j); | |
313 | maxBin = j; | |
314 | } | |
315 | ||
316 | // npart sum to 1 | |
95e970ca | 317 | correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);// * correlation->GetXaxis()->GetBinWidth(i)); |
19442b86 | 318 | correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum); |
319 | ||
320 | if (i <= fgMaxParams && j <= fgMaxInput) | |
95e970ca | 321 | { |
19442b86 | 322 | (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j); |
95e970ca | 323 | (*fgCorrelationMatrixSquared)(j-1, i-1) = correlation->GetBinContent(i, j) * correlation->GetBinContent(i, j); |
324 | } | |
19442b86 | 325 | } |
326 | ||
327 | //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin)); | |
328 | } | |
329 | ||
330 | //normalize measured | |
95e970ca | 331 | Float_t smallestError = 1; |
19442b86 | 332 | if (fgNormalizeInput) |
95e970ca | 333 | { |
334 | Float_t sumMeasured = measured->Integral(); | |
335 | measured->Scale(1.0 / sumMeasured); | |
336 | smallestError /= sumMeasured; | |
337 | } | |
19442b86 | 338 | |
339 | for (Int_t i=0; i<fgMaxInput; ++i) | |
340 | { | |
341 | (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1); | |
342 | if (measured->GetBinError(i+1) > 0) | |
95e970ca | 343 | { |
19442b86 | 344 | (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1); |
95e970ca | 345 | } |
346 | else // in this case put error of 1, otherwise 0 bins are not added to the chi2... | |
347 | (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / smallestError / smallestError; | |
19442b86 | 348 | |
349 | if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7) | |
350 | (*fgCorrelationCovarianceMatrix)(i, i) = 0; | |
351 | //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i)); | |
352 | } | |
95e970ca | 353 | |
354 | // efficiency is expected to match bin width of result | |
355 | for (Int_t i=0; i<fgMaxParams; ++i) | |
356 | { | |
357 | (*fgEfficiency)(i) = efficiency->GetBinContent(i+1); | |
95e970ca | 358 | } |
9e065ad2 | 359 | |
360 | if (correlation->GetNbinsX() != fgMaxParams || correlation->GetNbinsY() != fgMaxInput) | |
361 | cout << "Response histo has incorrect dimensions; expect (" << fgMaxParams << ", " << fgMaxInput << "), got (" << correlation->GetNbinsX() << ", " << correlation->GetNbinsY() << ")" << endl; | |
362 | ||
19442b86 | 363 | } |
364 | ||
365 | //____________________________________________________________________ | |
366 | Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check) | |
367 | { | |
368 | // | |
369 | // implementation of unfolding (internal function) | |
370 | // | |
371 | // unfolds <measured> using response from <correlation> and effiency <efficiency> | |
372 | // output is in <result> | |
373 | // <initialConditions> set the initial values for the minimization, if 0 <measured> is used | |
95e970ca | 374 | // negative values in initialConditions mean that the given parameter is fixed to the absolute of the value |
19442b86 | 375 | // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed |
376 | // | |
377 | // returns minuit status (0 = success), (-1 when check was set) | |
378 | // | |
379 | ||
95e970ca | 380 | SetStaticVariables(correlation, measured, efficiency); |
19442b86 | 381 | |
382 | // Initialize TMinuit via generic fitter interface | |
95e970ca | 383 | Int_t params = fgMaxParams; |
384 | if (fgNotFoundEvents > 0) | |
385 | params++; | |
386 | ||
387 | TVirtualFitter *minuit = TVirtualFitter::Fitter(0, params); | |
19442b86 | 388 | Double_t arglist[100]; |
9e065ad2 | 389 | // minuit->SetDefaultFitter("Minuit2"); |
19442b86 | 390 | |
391 | // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way... | |
392 | arglist[0] = 0; | |
393 | minuit->ExecuteCommand("SET PRINT", arglist, 1); | |
394 | ||
395 | // however, enable warnings | |
396 | //minuit->ExecuteCommand("SET WAR", arglist, 0); | |
397 | ||
398 | // set minimization function | |
399 | minuit->SetFCN(Chi2Function); | |
400 | ||
9e065ad2 | 401 | // set precision |
402 | minuit->SetPrecision(fgMinuitPrecision); | |
403 | ||
d22beb7f | 404 | minuit->SetMaxIterations(fgMinuitMaxIterations); |
405 | ||
9f070ef5 | 406 | minuit->ExecuteCommand("SET STRATEGY",&fgMinuitStrategy,1); |
407 | ||
19442b86 | 408 | for (Int_t i=0; i<fgMaxParams; i++) |
409 | (*fgEntropyAPriori)[i] = 1; | |
410 | ||
95e970ca | 411 | // set initial conditions as a-priori distribution for MRX regularization |
412 | /* | |
413 | for (Int_t i=0; i<fgMaxParams; i++) | |
414 | if (initialConditions && initialConditions->GetBinContent(i+1) > 0) | |
415 | (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1); | |
416 | */ | |
417 | ||
19442b86 | 418 | if (!initialConditions) { |
419 | initialConditions = measured; | |
420 | } else { | |
421 | Printf("AliUnfolding::UnfoldWithMinuit: Using different initial conditions..."); | |
422 | //new TCanvas; initialConditions->DrawCopy(); | |
423 | if (fgNormalizeInput) | |
424 | initialConditions->Scale(1.0 / initialConditions->Integral()); | |
425 | } | |
426 | ||
95e970ca | 427 | // extract minimum value from initial conditions (if we set a value to 0 it will stay 0) |
deaac8b1 | 428 | Float_t minValue = 1e35; |
95e970ca | 429 | if (fgMinimumInitialValueFix < 0) |
430 | { | |
431 | for (Int_t i=0; i<fgMaxParams; ++i) | |
432 | { | |
433 | Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1)); | |
434 | if (initialConditions->GetBinContent(bin) > 0) | |
435 | minValue = TMath::Min(minValue, (Float_t) initialConditions->GetBinContent(bin)); | |
436 | } | |
437 | } | |
438 | else | |
439 | minValue = fgMinimumInitialValueFix; | |
440 | ||
19442b86 | 441 | Double_t* results = new Double_t[fgMaxParams+1]; |
442 | for (Int_t i=0; i<fgMaxParams; ++i) | |
443 | { | |
95e970ca | 444 | Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1)); |
445 | results[i] = initialConditions->GetBinContent(bin); | |
19442b86 | 446 | |
95e970ca | 447 | Bool_t fix = kFALSE; |
448 | if (results[i] < 0) | |
449 | { | |
450 | fix = kTRUE; | |
451 | results[i] = -results[i]; | |
452 | } | |
453 | ||
454 | if (!fix && fgMinimumInitialValue && results[i] < minValue) | |
455 | results[i] = minValue; | |
456 | ||
19442b86 | 457 | // minuit sees squared values to prevent it from going negative... |
458 | results[i] = TMath::Sqrt(results[i]); | |
459 | ||
95e970ca | 460 | minuit->SetParameter(i, Form("param%d", i), results[i], (fix) ? 0 : fgMinuitStepSize, 0, 0); |
461 | } | |
462 | if (fgNotFoundEvents > 0) | |
463 | { | |
464 | results[fgMaxParams] = efficiency->GetBinContent(1); | |
465 | minuit->SetParameter(fgMaxParams, "vtx0", results[fgMaxParams], fgMinuitStepSize / 100, 0.01, 0.80); | |
19442b86 | 466 | } |
467 | ||
468 | Int_t dummy = 0; | |
469 | Double_t chi2 = 0; | |
470 | Chi2Function(dummy, 0, chi2, results, 0); | |
471 | printf("AliUnfolding::UnfoldWithMinuit: Chi2 of initial parameters is = %f\n", chi2); | |
472 | ||
473 | if (check) | |
474 | { | |
475 | DrawGuess(results); | |
476 | delete[] results; | |
477 | return -1; | |
478 | } | |
479 | ||
480 | // first param is number of iterations, second is precision.... | |
03650114 | 481 | arglist[0] = (float)fgMinuitMaxIterations; |
482 | // arglist[1] = 1e-5; | |
9e065ad2 | 483 | // minuit->ExecuteCommand("SET PRINT", arglist, 3); |
484 | // minuit->ExecuteCommand("SCAN", arglist, 0); | |
19442b86 | 485 | Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1); |
486 | Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status); | |
487 | //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!"); | |
488 | //minuit->ExecuteCommand("MINOS", arglist, 0); | |
489 | ||
95e970ca | 490 | if (fgNotFoundEvents > 0) |
491 | { | |
492 | results[fgMaxParams] = minuit->GetParameter(fgMaxParams); | |
493 | Printf("Efficiency for bin 0 changed from %f to %f", efficiency->GetBinContent(1), results[fgMaxParams]); | |
494 | efficiency->SetBinContent(1, results[fgMaxParams]); | |
495 | } | |
496 | ||
19442b86 | 497 | for (Int_t i=0; i<fgMaxParams; ++i) |
498 | { | |
499 | results[i] = minuit->GetParameter(i); | |
500 | Double_t value = results[i] * results[i]; | |
a9edd311 | 501 | // error is : 2 * (relError on results[i]) * (value) = 2 * (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i)) |
b203a518 | 502 | Double_t error = 0; |
503 | if (TMath::IsNaN(minuit->GetParError(i))) | |
504 | Printf("WARNING: Parameter %d error is nan", i); | |
505 | else | |
a9edd311 | 506 | error = 2 * minuit->GetParError(i) * results[i]; |
19442b86 | 507 | |
508 | if (efficiency) | |
95e970ca | 509 | { |
9e065ad2 | 510 | //printf("value before efficiency correction: %f\n",value); |
19442b86 | 511 | if (efficiency->GetBinContent(i+1) > 0) |
512 | { | |
95e970ca | 513 | value /= efficiency->GetBinContent(i+1); |
514 | error /= efficiency->GetBinContent(i+1); | |
19442b86 | 515 | } |
516 | else | |
517 | { | |
95e970ca | 518 | value = 0; |
519 | error = 0; | |
19442b86 | 520 | } |
521 | } | |
9e065ad2 | 522 | //printf("value after efficiency correction: %f +/- %f\n",value,error); |
19442b86 | 523 | result->SetBinContent(i+1, value); |
524 | result->SetBinError(i+1, error); | |
525 | } | |
d22beb7f | 526 | |
527 | Int_t tmpCallCount = fgCallCount; | |
528 | fgCallCount = 0; // needs to be 0 so that the Chi2Function prints its output | |
95e970ca | 529 | Chi2Function(dummy, 0, chi2, results, 0); |
d22beb7f | 530 | |
531 | Printf("AliUnfolding::UnfoldWithMinuit: iterations %d. Chi2 of final parameters is = %f", tmpCallCount, chi2); | |
95e970ca | 532 | |
19442b86 | 533 | delete[] results; |
534 | ||
535 | return status; | |
536 | } | |
537 | ||
538 | //____________________________________________________________________ | |
539 | Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult) | |
540 | { | |
541 | // | |
542 | // unfolds a spectrum using the Bayesian method | |
543 | // | |
544 | ||
545 | if (measured->Integral() <= 0) | |
546 | { | |
547 | Printf("AliUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty"); | |
548 | return -1; | |
549 | } | |
550 | ||
551 | const Int_t kStartBin = 0; | |
552 | ||
553 | Int_t kMaxM = fgMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis | |
554 | Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis | |
555 | ||
556 | // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4) | |
557 | const Double_t kConvergenceLimit = kMaxT * 1e-6; | |
558 | ||
559 | // store information in arrays, to increase processing speed (~ factor 5) | |
560 | Double_t* measuredCopy = new Double_t[kMaxM]; | |
561 | Double_t* measuredError = new Double_t[kMaxM]; | |
562 | Double_t* prior = new Double_t[kMaxT]; | |
563 | Double_t* result = new Double_t[kMaxT]; | |
564 | Double_t* efficiency = new Double_t[kMaxT]; | |
95e970ca | 565 | Double_t* binWidths = new Double_t[kMaxT]; |
736125ad | 566 | Double_t* normResponse = new Double_t[kMaxT]; |
19442b86 | 567 | |
568 | Double_t** response = new Double_t*[kMaxT]; | |
569 | Double_t** inverseResponse = new Double_t*[kMaxT]; | |
570 | for (Int_t i=0; i<kMaxT; i++) | |
571 | { | |
572 | response[i] = new Double_t[kMaxM]; | |
573 | inverseResponse[i] = new Double_t[kMaxM]; | |
736125ad | 574 | normResponse[i] = 0; |
19442b86 | 575 | } |
576 | ||
577 | // for normalization | |
578 | Float_t measuredIntegral = measured->Integral(); | |
579 | for (Int_t m=0; m<kMaxM; m++) | |
580 | { | |
581 | measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral; | |
582 | measuredError[m] = measured->GetBinError(m+1) / measuredIntegral; | |
583 | ||
584 | for (Int_t t=0; t<kMaxT; t++) | |
585 | { | |
586 | response[t][m] = correlation->GetBinContent(t+1, m+1); | |
587 | inverseResponse[t][m] = 0; | |
736125ad | 588 | normResponse[t] += correlation->GetBinContent(t+1, m+1); |
19442b86 | 589 | } |
590 | } | |
591 | ||
592 | for (Int_t t=0; t<kMaxT; t++) | |
593 | { | |
651e2088 | 594 | if (aEfficiency) |
19442b86 | 595 | { |
596 | efficiency[t] = aEfficiency->GetBinContent(t+1); | |
597 | } | |
598 | else | |
599 | efficiency[t] = 1; | |
600 | ||
601 | prior[t] = measuredCopy[t]; | |
602 | result[t] = 0; | |
95e970ca | 603 | binWidths[t] = aResult->GetXaxis()->GetBinWidth(t+1); |
736125ad | 604 | |
605 | for (Int_t m=0; m<kMaxM; m++) { // Normalise response matrix | |
606 | if (normResponse[t] != 0) | |
607 | response[t][m] /= normResponse[t]; | |
608 | else | |
609 | Printf("AliUnfolding::UnfoldWithBayesian: Empty row,column in response matrix, for truth bin %d",t); | |
610 | } | |
19442b86 | 611 | } |
612 | ||
613 | // pick prior distribution | |
614 | if (initialConditions) | |
615 | { | |
616 | printf("Using different starting conditions...\n"); | |
617 | // for normalization | |
618 | Float_t inputDistIntegral = initialConditions->Integral(); | |
619 | for (Int_t i=0; i<kMaxT; i++) | |
620 | prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral; | |
621 | } | |
622 | ||
623 | //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5); | |
624 | ||
95e970ca | 625 | //new TCanvas; |
19442b86 | 626 | // unfold... |
95e970ca | 627 | for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations < 0; i++) |
19442b86 | 628 | { |
629 | if (fgDebug) | |
630 | Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i); | |
631 | ||
95e970ca | 632 | // calculate IR from Bayes theorem |
19442b86 | 633 | // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k) |
634 | ||
635 | Double_t chi2Measured = 0; | |
636 | for (Int_t m=0; m<kMaxM; m++) | |
637 | { | |
638 | Float_t norm = 0; | |
639 | for (Int_t t = kStartBin; t<kMaxT; t++) | |
736125ad | 640 | norm += response[t][m] * prior[t] * efficiency[t]; |
19442b86 | 641 | |
642 | // calc. chi2: (measured - response * prior) / error | |
643 | if (measuredError[m] > 0) | |
644 | { | |
645 | Double_t value = (measuredCopy[m] - norm) / measuredError[m]; | |
646 | chi2Measured += value * value; | |
647 | } | |
648 | ||
649 | if (norm > 0) | |
650 | { | |
651 | for (Int_t t = kStartBin; t<kMaxT; t++) | |
736125ad | 652 | inverseResponse[t][m] = response[t][m] * prior[t] * efficiency[t] / norm; |
19442b86 | 653 | } |
654 | else | |
655 | { | |
656 | for (Int_t t = kStartBin; t<kMaxT; t++) | |
657 | inverseResponse[t][m] = 0; | |
658 | } | |
659 | } | |
660 | //Printf("chi2Measured of the last prior is %e", chi2Measured); | |
661 | ||
662 | for (Int_t t = kStartBin; t<kMaxT; t++) | |
663 | { | |
664 | Float_t value = 0; | |
665 | for (Int_t m=0; m<kMaxM; m++) | |
666 | value += inverseResponse[t][m] * measuredCopy[m]; | |
667 | ||
668 | if (efficiency[t] > 0) | |
669 | result[t] = value / efficiency[t]; | |
670 | else | |
671 | result[t] = 0; | |
672 | } | |
673 | ||
95e970ca | 674 | /* |
19442b86 | 675 | // draw intermediate result |
19442b86 | 676 | for (Int_t t=0; t<kMaxT; t++) |
95e970ca | 677 | { |
19442b86 | 678 | aResult->SetBinContent(t+1, result[t]); |
95e970ca | 679 | } |
680 | aResult->SetMarkerStyle(24+i); | |
19442b86 | 681 | aResult->SetMarkerColor(2); |
95e970ca | 682 | aResult->DrawCopy((i == 0) ? "P" : "PSAME"); |
19442b86 | 683 | */ |
95e970ca | 684 | |
19442b86 | 685 | Double_t chi2LastIter = 0; |
686 | // regularization (simple smoothing) | |
687 | for (Int_t t=kStartBin; t<kMaxT; t++) | |
688 | { | |
689 | Float_t newValue = 0; | |
690 | ||
691 | // 0 bin excluded from smoothing | |
692 | if (t > kStartBin+2 && t<kMaxT-1) | |
693 | { | |
95e970ca | 694 | Float_t average = (result[t-1] / binWidths[t-1] + result[t] / binWidths[t] + result[t+1] / binWidths[t+1]) / 3 * binWidths[t]; |
19442b86 | 695 | |
696 | // weight the average with the regularization parameter | |
697 | newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average; | |
698 | } | |
699 | else | |
700 | newValue = result[t]; | |
701 | ||
702 | // calculate chi2 (change from last iteration) | |
703 | if (prior[t] > 1e-5) | |
704 | { | |
705 | Double_t diff = (prior[t] - newValue) / prior[t]; | |
706 | chi2LastIter += diff * diff; | |
707 | } | |
708 | ||
709 | prior[t] = newValue; | |
710 | } | |
711 | //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter); | |
712 | //convergence->Fill(i+1, chi2LastIter); | |
713 | ||
714 | if (fgBayesianIterations < 0 && chi2LastIter < kConvergenceLimit) | |
715 | { | |
716 | 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); | |
717 | break; | |
718 | } | |
719 | } // end of iterations | |
720 | ||
721 | //new TCanvas; convergence->DrawCopy(); gPad->SetLogy(); | |
722 | //delete convergence; | |
723 | ||
95e970ca | 724 | Float_t factor = 1; |
725 | if (!fgNormalizeInput) | |
726 | factor = measuredIntegral; | |
19442b86 | 727 | for (Int_t t=0; t<kMaxT; t++) |
95e970ca | 728 | aResult->SetBinContent(t+1, result[t] * factor); |
19442b86 | 729 | |
730 | delete[] measuredCopy; | |
731 | delete[] measuredError; | |
732 | delete[] prior; | |
733 | delete[] result; | |
734 | delete[] efficiency; | |
95e970ca | 735 | delete[] binWidths; |
736125ad | 736 | delete[] normResponse; |
19442b86 | 737 | |
738 | for (Int_t i=0; i<kMaxT; i++) | |
739 | { | |
740 | delete[] response[i]; | |
741 | delete[] inverseResponse[i]; | |
742 | } | |
743 | delete[] response; | |
744 | delete[] inverseResponse; | |
745 | ||
746 | return 0; | |
747 | ||
748 | // ******** | |
749 | // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995 | |
750 | ||
751 | /*printf("Calculating covariance matrix. This may take some time...\n"); | |
752 | ||
753 | // check if this is the right one... | |
754 | TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX()); | |
755 | ||
756 | Int_t xBins = hInverseResponseBayes->GetNbinsX(); | |
757 | Int_t yBins = hInverseResponseBayes->GetNbinsY(); | |
758 | ||
759 | // calculate "unfolding matrix" Mij | |
760 | Float_t matrixM[251][251]; | |
761 | for (Int_t i=1; i<=xBins; i++) | |
762 | { | |
763 | for (Int_t j=1; j<=yBins; j++) | |
764 | { | |
765 | if (fCurrentEfficiency->GetBinContent(i) > 0) | |
766 | matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i); | |
767 | else | |
768 | matrixM[i-1][j-1] = 0; | |
769 | } | |
770 | } | |
771 | ||
772 | Float_t* vectorn = new Float_t[yBins]; | |
773 | for (Int_t j=1; j<=yBins; j++) | |
774 | vectorn[j-1] = fCurrentESD->GetBinContent(j); | |
775 | ||
776 | // first part of covariance matrix, depends on input distribution n(E) | |
777 | Float_t cov1[251][251]; | |
778 | ||
779 | Float_t nEvents = fCurrentESD->Integral(); // N | |
780 | ||
781 | xBins = 20; | |
782 | yBins = 20; | |
783 | ||
784 | for (Int_t k=0; k<xBins; k++) | |
785 | { | |
786 | printf("In Cov1: %d\n", k); | |
787 | for (Int_t l=0; l<yBins; l++) | |
788 | { | |
789 | cov1[k][l] = 0; | |
790 | ||
791 | // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N) | |
792 | for (Int_t j=0; j<yBins; j++) | |
793 | cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j] | |
794 | * (1.0 - vectorn[j] / nEvents); | |
795 | ||
796 | // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N | |
797 | for (Int_t i=0; i<yBins; i++) | |
798 | for (Int_t j=0; j<yBins; j++) | |
799 | { | |
800 | if (i == j) | |
801 | continue; | |
802 | cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i] | |
803 | * vectorn[j] / nEvents; | |
804 | } | |
805 | } | |
806 | } | |
807 | ||
808 | printf("Cov1 finished\n"); | |
809 | ||
810 | TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov"); | |
811 | cov->Reset(); | |
812 | ||
813 | for (Int_t i=1; i<=xBins; i++) | |
814 | for (Int_t j=1; j<=yBins; j++) | |
815 | cov->SetBinContent(i, j, cov1[i-1][j-1]); | |
816 | ||
817 | new TCanvas; | |
818 | cov->Draw("COLZ"); | |
819 | ||
820 | // second part of covariance matrix, depends on response matrix | |
821 | Float_t cov2[251][251]; | |
822 | ||
823 | // Cov[P(Er|Cu), P(Es|Cu)] term | |
824 | Float_t covTerm[100][100][100]; | |
825 | for (Int_t r=0; r<yBins; r++) | |
826 | for (Int_t u=0; u<xBins; u++) | |
827 | for (Int_t s=0; s<yBins; s++) | |
828 | { | |
829 | if (r == s) | |
830 | covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) | |
831 | * (1.0 - hResponse->GetBinContent(u+1, r+1)); | |
832 | else | |
833 | covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) | |
834 | * hResponse->GetBinContent(u+1, s+1); | |
835 | } | |
836 | ||
837 | for (Int_t k=0; k<xBins; k++) | |
838 | for (Int_t l=0; l<yBins; l++) | |
839 | { | |
840 | cov2[k][l] = 0; | |
841 | printf("In Cov2: %d %d\n", k, l); | |
842 | for (Int_t i=0; i<yBins; i++) | |
843 | for (Int_t j=0; j<yBins; j++) | |
844 | { | |
845 | //printf("In Cov2: %d %d %d %d\n", k, l, i, j); | |
846 | // calculate Cov(Mki, Mlj) = sum{ru},{su} ... | |
847 | Float_t tmpCov = 0; | |
848 | for (Int_t r=0; r<yBins; r++) | |
849 | for (Int_t u=0; u<xBins; u++) | |
850 | for (Int_t s=0; s<yBins; s++) | |
851 | { | |
852 | if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0 | |
853 | || hResponse->GetBinContent(u+1, i+1) == 0) | |
854 | continue; | |
855 | ||
856 | tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u) | |
857 | * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u) | |
858 | * covTerm[r][u][s]; | |
859 | } | |
860 | ||
861 | cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov; | |
862 | } | |
863 | } | |
864 | ||
865 | printf("Cov2 finished\n"); | |
866 | ||
867 | for (Int_t i=1; i<=xBins; i++) | |
868 | for (Int_t j=1; j<=yBins; j++) | |
869 | cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]); | |
870 | ||
871 | new TCanvas; | |
872 | cov->Draw("COLZ");*/ | |
873 | } | |
874 | ||
875 | //____________________________________________________________________ | |
876 | Double_t AliUnfolding::RegularizationPol0(TVectorD& params) | |
877 | { | |
878 | // homogenity term for minuit fitting | |
879 | // pure function of the parameters | |
880 | // prefers constant function (pol0) | |
9e065ad2 | 881 | // |
882 | // Does not take into account efficiency | |
19442b86 | 883 | Double_t chi2 = 0; |
884 | ||
885 | for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
886 | { | |
9e065ad2 | 887 | Double_t right = params[i] / fgUnfoldedAxis->GetBinWidth(i+1); |
888 | Double_t left = params[i-1] / fgUnfoldedAxis->GetBinWidth(i); | |
19442b86 | 889 | |
95e970ca | 890 | if (left != 0) |
19442b86 | 891 | { |
95e970ca | 892 | Double_t diff = (right - left); |
9e065ad2 | 893 | chi2 += diff * diff / left / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2); |
19442b86 | 894 | } |
895 | } | |
896 | ||
897 | return chi2 / 100.0; | |
898 | } | |
899 | ||
900 | //____________________________________________________________________ | |
901 | Double_t AliUnfolding::RegularizationPol1(TVectorD& params) | |
902 | { | |
903 | // homogenity term for minuit fitting | |
904 | // pure function of the parameters | |
905 | // prefers linear function (pol1) | |
9e065ad2 | 906 | // |
907 | // Does not take into account efficiency | |
19442b86 | 908 | Double_t chi2 = 0; |
909 | ||
910 | for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
911 | { | |
912 | if (params[i-1] == 0) | |
913 | continue; | |
914 | ||
9e065ad2 | 915 | Double_t right = params[i] / fgUnfoldedAxis->GetBinWidth(i+1); |
916 | Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i); | |
917 | Double_t left = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1); | |
19442b86 | 918 | |
9e065ad2 | 919 | Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2); |
920 | Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2); | |
19442b86 | 921 | |
95e970ca | 922 | //Double_t diff = (der1 - der2) / middle; |
923 | //chi2 += diff * diff; | |
9e065ad2 | 924 | chi2 += (der1 - der2) * (der1 - der2) / middle * fgUnfoldedAxis->GetBinWidth(i); |
19442b86 | 925 | } |
926 | ||
927 | return chi2; | |
928 | } | |
929 | ||
930 | //____________________________________________________________________ | |
931 | Double_t AliUnfolding::RegularizationLog(TVectorD& params) | |
932 | { | |
933 | // homogenity term for minuit fitting | |
934 | // pure function of the parameters | |
9e065ad2 | 935 | // prefers logarithmic function (log) |
936 | // | |
937 | // Does not take into account efficiency | |
19442b86 | 938 | |
939 | Double_t chi2 = 0; | |
940 | ||
941 | for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
942 | { | |
95e970ca | 943 | if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0) |
944 | continue; | |
19442b86 | 945 | |
9e065ad2 | 946 | Double_t right = log(params[i] / fgUnfoldedAxis->GetBinWidth(i+1)); |
947 | Double_t middle = log(params[i-1] / fgUnfoldedAxis->GetBinWidth(i)); | |
948 | Double_t left = log(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1)); | |
95e970ca | 949 | |
9e065ad2 | 950 | Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2); |
951 | Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2); | |
95e970ca | 952 | |
953 | //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2]; | |
19442b86 | 954 | |
95e970ca | 955 | //if (fgCallCount == 0) |
956 | // Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error); | |
957 | chi2 += (der1 - der2) * (der1 - der2);// / error; | |
19442b86 | 958 | } |
959 | ||
95e970ca | 960 | return chi2; |
19442b86 | 961 | } |
962 | ||
963 | //____________________________________________________________________ | |
964 | Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params) | |
965 | { | |
966 | // homogenity term for minuit fitting | |
967 | // pure function of the parameters | |
968 | // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments, | |
969 | // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp. | |
9e065ad2 | 970 | // |
971 | // Does not take into account efficiency | |
19442b86 | 972 | |
973 | Double_t chi2 = 0; | |
974 | ||
975 | for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
976 | { | |
977 | Double_t right = params[i]; | |
978 | Double_t middle = params[i-1]; | |
979 | Double_t left = params[i-2]; | |
980 | ||
981 | Double_t der1 = (right - middle); | |
982 | Double_t der2 = (middle - left); | |
983 | ||
984 | Double_t diff = (der1 - der2); | |
985 | ||
986 | chi2 += diff * diff; | |
987 | } | |
988 | ||
989 | return chi2 * 1e4; | |
990 | } | |
991 | ||
992 | //____________________________________________________________________ | |
993 | Double_t AliUnfolding::RegularizationEntropy(TVectorD& params) | |
994 | { | |
995 | // homogenity term for minuit fitting | |
996 | // pure function of the parameters | |
997 | // calculates entropy, from | |
998 | // The method of reduced cross-entropy (M. Schmelling 1993) | |
9e065ad2 | 999 | // |
1000 | // Does not take into account efficiency | |
19442b86 | 1001 | |
1002 | Double_t paramSum = 0; | |
1003 | ||
1004 | for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i) | |
1005 | paramSum += params[i]; | |
1006 | ||
1007 | Double_t chi2 = 0; | |
1008 | for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i) | |
1009 | { | |
95e970ca | 1010 | Double_t tmp = params[i] / paramSum; |
1011 | //Double_t tmp = params[i]; | |
19442b86 | 1012 | if (tmp > 0 && (*fgEntropyAPriori)[i] > 0) |
1013 | { | |
1014 | chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]); | |
1015 | } | |
95e970ca | 1016 | else |
1017 | chi2 += 100; | |
1018 | } | |
1019 | ||
1020 | return -chi2; | |
1021 | } | |
1022 | ||
1023 | //____________________________________________________________________ | |
1024 | Double_t AliUnfolding::RegularizationRatio(TVectorD& params) | |
1025 | { | |
1026 | // homogenity term for minuit fitting | |
1027 | // pure function of the parameters | |
9e065ad2 | 1028 | // |
1029 | // Does not take into account efficiency | |
95e970ca | 1030 | |
1031 | Double_t chi2 = 0; | |
1032 | ||
1033 | for (Int_t i=5+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
1034 | { | |
1035 | if (params[i-1] == 0 || params[i] == 0) | |
1036 | continue; | |
1037 | ||
9e065ad2 | 1038 | Double_t right = params[i] / fgUnfoldedAxis->GetBinWidth(i+1); |
1039 | Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i); | |
1040 | Double_t left = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1); | |
1041 | Double_t left2 = params[i-3] / fgUnfoldedAxis->GetBinWidth(i-2); | |
1042 | Double_t left3 = params[i-4] / fgUnfoldedAxis->GetBinWidth(i-3); | |
1043 | Double_t left4 = params[i-5] / fgUnfoldedAxis->GetBinWidth(i-4); | |
95e970ca | 1044 | |
1045 | //Double_t diff = left / middle - middle / right; | |
1046 | //Double_t diff = 2 * left / middle - middle / right - left2 / left; | |
1047 | Double_t diff = 4 * left2 / left - middle / right - left / middle - left3 / left2 - left4 / left3; | |
1048 | ||
1049 | chi2 += diff * diff;// / middle; | |
19442b86 | 1050 | } |
1051 | ||
95e970ca | 1052 | return chi2; |
19442b86 | 1053 | } |
1054 | ||
9e065ad2 | 1055 | //____________________________________________________________________ |
1056 | Double_t AliUnfolding::RegularizationPowerLaw(TVectorD& params) | |
1057 | { | |
1058 | // homogenity term for minuit fitting | |
1059 | // pure function of the parameters | |
1060 | // prefers power law with n = -5 | |
1061 | // | |
1062 | // Does not take into account efficiency | |
1063 | ||
1064 | Double_t chi2 = 0; | |
1065 | ||
1066 | Double_t right = 0.; | |
1067 | Double_t middle = 0.; | |
1068 | Double_t left = 0.; | |
1069 | ||
1070 | for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
1071 | { | |
1072 | if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8) | |
1073 | continue; | |
1074 | ||
1075 | if (fgUnfoldedAxis->GetBinWidth(i+1) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i-1) < 1e-8) | |
1076 | continue; | |
1077 | ||
1078 | middle = TMath::Power(params[i-1] / fgUnfoldedAxis->GetBinWidth(i),fgPowern); | |
1079 | ||
1080 | if(middle>0) { | |
59076376 | 1081 | right = TMath::Power(params[i] / fgUnfoldedAxis->GetBinWidth(i+1),fgPowern)/middle; |
9e065ad2 | 1082 | |
1083 | left = TMath::Power(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1),fgPowern)/middle; | |
1084 | ||
1085 | middle = 1.; | |
1086 | ||
1087 | Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2); | |
59076376 | 1088 | Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2); |
9e065ad2 | 1089 | |
59076376 | 1090 | chi2 += (der1 - der2) * (der1 - der2)/ ( fgUnfoldedAxis->GetBinWidth(i+1)/2. + fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)/2.)/( fgUnfoldedAxis->GetBinWidth(i+1)/2. + fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)/2.);// / error; |
9e065ad2 | 1091 | // printf("i: %d chi2 = %f\n",i,chi2); |
1092 | } | |
1093 | ||
1094 | } | |
1095 | ||
1096 | return chi2; | |
1097 | } | |
1098 | ||
1099 | //____________________________________________________________________ | |
1100 | Double_t AliUnfolding::RegularizationLogLog(TVectorD& params) | |
1101 | { | |
1102 | // homogenity term for minuit fitting | |
1103 | // pure function of the parameters | |
1104 | // prefers a powerlaw (linear on a log-log scale) | |
1105 | // | |
1106 | // The calculation takes into account the efficiencies | |
1107 | ||
1108 | Double_t chi2 = 0; | |
1109 | ||
1110 | for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
1111 | { | |
1112 | if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0) | |
1113 | continue; | |
1114 | if ((*fgEfficiency)(i-1) == 0 || (*fgEfficiency)(i) == 0 || (*fgEfficiency)(i-2) == 0) | |
1115 | continue; | |
1116 | ||
59076376 | 1117 | Double_t right = log(params[i] / (*fgEfficiency)(i) / fgUnfoldedAxis->GetBinWidth(i+1)); |
1118 | Double_t middle = log(params[i-1] / (*fgEfficiency)(i-1) / fgUnfoldedAxis->GetBinWidth(i)); | |
1119 | Double_t left = log(params[i-2] / (*fgEfficiency)(i-2) / fgUnfoldedAxis->GetBinWidth(i-1)); | |
9e065ad2 | 1120 | |
1121 | Double_t der1 = (right - middle) / ( log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i)) ); | |
1122 | Double_t der2 = (middle - left) /( log(fgUnfoldedAxis->GetBinCenter(i)) - log(fgUnfoldedAxis->GetBinCenter(i-1)) ); | |
59076376 | 1123 | |
9e065ad2 | 1124 | double tmp = (log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i-1)))/2.; |
1125 | Double_t dder = (der1-der2) / tmp; | |
1126 | ||
1127 | chi2 += dder * dder; | |
1128 | } | |
1129 | ||
1130 | return chi2; | |
1131 | } | |
1132 | ||
1133 | ||
1134 | ||
19442b86 | 1135 | //____________________________________________________________________ |
1136 | void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t) | |
1137 | { | |
1138 | // | |
1139 | // fit function for minuit | |
1140 | // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix | |
1141 | // | |
95e970ca | 1142 | |
1143 | // TODO use static members for the variables here to speed up processing (no construction/deconstruction) | |
19442b86 | 1144 | |
9e065ad2 | 1145 | // d = guess |
95e970ca | 1146 | TVectorD paramsVector(fgMaxParams); |
19442b86 | 1147 | for (Int_t i=0; i<fgMaxParams; ++i) |
1148 | paramsVector[i] = params[i] * params[i]; | |
1149 | ||
1150 | // calculate penalty factor | |
1151 | Double_t penaltyVal = 0; | |
9e065ad2 | 1152 | |
19442b86 | 1153 | switch (fgRegularizationType) |
1154 | { | |
1155 | case kNone: break; | |
1156 | case kPol0: penaltyVal = RegularizationPol0(paramsVector); break; | |
1157 | case kPol1: penaltyVal = RegularizationPol1(paramsVector); break; | |
1158 | case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break; | |
1159 | case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break; | |
1160 | case kLog: penaltyVal = RegularizationLog(paramsVector); break; | |
95e970ca | 1161 | case kRatio: penaltyVal = RegularizationRatio(paramsVector); break; |
9e065ad2 | 1162 | case kPowerLaw: penaltyVal = RegularizationPowerLaw(paramsVector); break; |
1163 | case kLogLog: penaltyVal = RegularizationLogLog(paramsVector); break; | |
19442b86 | 1164 | } |
1165 | ||
9e065ad2 | 1166 | penaltyVal *= fgRegularizationWeight; //beta*PU |
19442b86 | 1167 | |
1168 | // Ad | |
1169 | TVectorD measGuessVector(fgMaxInput); | |
1170 | measGuessVector = (*fgCorrelationMatrix) * paramsVector; | |
1171 | ||
1172 | // Ad - m | |
1173 | measGuessVector -= (*fgCurrentESDVector); | |
95e970ca | 1174 | |
1175 | #if 0 | |
1176 | // new error calcuation using error on the guess | |
1177 | ||
1178 | // error from guess | |
1179 | TVectorD errorGuessVector(fgMaxInput); | |
1180 | //errorGuessVector = (*fgCorrelationMatrixSquared) * paramsVector; | |
1181 | errorGuessVector = (*fgCorrelationMatrix) * paramsVector; | |
19442b86 | 1182 | |
95e970ca | 1183 | Double_t chi2FromFit = 0; |
1184 | for (Int_t i=0; i<fgMaxInput; ++i) | |
1185 | { | |
1186 | Float_t error = 1; | |
1187 | if (errorGuessVector(i) > 0) | |
1188 | error = errorGuessVector(i); | |
1189 | chi2FromFit += measGuessVector(i) * measGuessVector(i) / error; | |
1190 | } | |
19442b86 | 1191 | |
95e970ca | 1192 | #else |
1193 | // old error calcuation using the error on the measured | |
19442b86 | 1194 | TVectorD copy(measGuessVector); |
1195 | ||
1196 | // (Ad - m) W | |
1197 | // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used | |
1198 | // normal way is like this: | |
1199 | // measGuessVector *= (*fgCorrelationCovarianceMatrix); | |
1200 | // optimized way like this: | |
1201 | for (Int_t i=0; i<fgMaxInput; ++i) | |
1202 | measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i); | |
1203 | ||
19442b86 | 1204 | |
95e970ca | 1205 | if (fgSkipBin0InChi2) |
1206 | measGuessVector[0] = 0; | |
1207 | ||
19442b86 | 1208 | // (Ad - m) W (Ad - m) |
1209 | // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very | |
95e970ca | 1210 | // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see UnfoldWithMinuit) |
19442b86 | 1211 | Double_t chi2FromFit = measGuessVector * copy * 1e6; |
95e970ca | 1212 | #endif |
19442b86 | 1213 | |
95e970ca | 1214 | Double_t notFoundEventsConstraint = 0; |
1215 | Double_t currentNotFoundEvents = 0; | |
1216 | Double_t errorNotFoundEvents = 0; | |
1217 | ||
1218 | if (fgNotFoundEvents > 0) | |
1219 | { | |
1220 | for (Int_t i=0; i<fgMaxParams; ++i) | |
1221 | { | |
1222 | Double_t eff = (1.0 / (*fgEfficiency)(i) - 1); | |
1223 | if (i == 0) | |
1224 | eff = (1.0 / params[fgMaxParams] - 1); | |
1225 | currentNotFoundEvents += eff * paramsVector(i); | |
1226 | errorNotFoundEvents += eff * eff * paramsVector(i); // error due to guess (paramsVector) | |
1227 | if (i <= 3) | |
1228 | errorNotFoundEvents += (eff * 0.03) * (eff * 0.03) * paramsVector(i) * paramsVector(i); // error on eff | |
9e065ad2 | 1229 | // if ((fgCallCount % 10000) == 0) |
1230 | //Printf("%d %f %f %f", i, (*fgEfficiency)(i), paramsVector(i), currentNotFoundEvents); | |
95e970ca | 1231 | } |
1232 | errorNotFoundEvents += fgNotFoundEvents; | |
1233 | // TODO add error on background, background estimate | |
1234 | ||
1235 | notFoundEventsConstraint = (currentNotFoundEvents - fgNotFoundEvents) * (currentNotFoundEvents - fgNotFoundEvents) / errorNotFoundEvents; | |
1236 | } | |
1237 | ||
1238 | Float_t avg = 0; | |
1239 | Float_t sum = 0; | |
1240 | Float_t currentMult = 0; | |
1241 | // try to match dn/deta | |
1242 | for (Int_t i=0; i<fgMaxParams; ++i) | |
1243 | { | |
1244 | avg += paramsVector(i) * currentMult; | |
1245 | sum += paramsVector(i); | |
9e065ad2 | 1246 | currentMult += fgUnfoldedAxis->GetBinWidth(i); |
95e970ca | 1247 | } |
1248 | avg /= sum; | |
1249 | Float_t chi2Avg = 0; //(avg - 3.73) * (avg - 3.73) * 100; | |
19442b86 | 1250 | |
95e970ca | 1251 | chi2 = chi2FromFit + penaltyVal + notFoundEventsConstraint + chi2Avg; |
9e065ad2 | 1252 | |
95e970ca | 1253 | if ((fgCallCount++ % 1000) == 0) |
1254 | { | |
9e065ad2 | 1255 | |
1256 | 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-1], avg, chi2FromFit, penaltyVal, notFoundEventsConstraint, chi2Avg, chi2); | |
1257 | ||
95e970ca | 1258 | //for (Int_t i=0; i<fgMaxInput; ++i) |
1259 | // Printf("%d: %f", i, measGuessVector(i) * copy(i) * 1e6); | |
1260 | } | |
9e065ad2 | 1261 | |
1262 | fChi2FromFit = chi2FromFit; | |
1263 | fPenaltyVal = penaltyVal; | |
19442b86 | 1264 | } |
1265 | ||
1266 | //____________________________________________________________________ | |
9e065ad2 | 1267 | void AliUnfolding::MakePads() { |
1268 | TPad *presult = new TPad("presult","result",0,0.4,1,1); | |
1269 | presult->SetNumber(1); | |
1270 | presult->Draw(); | |
1271 | presult->SetLogy(); | |
1272 | TPad *pres = new TPad("pres","residuals",0,0.2,1,0.4); | |
1273 | pres->SetNumber(2); | |
1274 | pres->Draw(); | |
1275 | TPad *ppen = new TPad("ppen","penalty",0,0.0,1,0.2); | |
1276 | ppen->SetNumber(3); | |
1277 | ppen->Draw(); | |
1278 | ||
1279 | } | |
1280 | //____________________________________________________________________ | |
1281 | void AliUnfolding::DrawResults(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TCanvas *canv, Int_t reuseHists,TH1 *unfolded) | |
1282 | { | |
1283 | // Draw histograms of | |
1284 | // - Result folded with response | |
1285 | // - Penalty factors | |
1286 | // - Chisquare contributions | |
1287 | // (Useful for debugging/sanity checks and the interactive unfolder) | |
1288 | // | |
1289 | // If a canvas pointer is given (canv != 0), it will be used for all | |
1290 | // plots; 3 pads are made if needed. | |
1291 | ||
1292 | ||
1293 | Int_t blankCanvas = 0; | |
1294 | TVirtualPad *presult = 0; | |
1295 | TVirtualPad *pres = 0; | |
1296 | TVirtualPad *ppen = 0; | |
1297 | ||
1298 | if (canv) { | |
1299 | canv->cd(); | |
1300 | presult = canv->GetPad(1); | |
1301 | pres = canv->GetPad(2); | |
1302 | ppen = canv->GetPad(3); | |
1303 | if (presult == 0 || pres == 0 || ppen == 0) { | |
1304 | canv->Clear(); | |
1305 | MakePads(); | |
1306 | blankCanvas = 1; | |
1307 | presult = canv->GetPad(1); | |
1308 | pres = canv->GetPad(2); | |
1309 | ppen = canv->GetPad(3); | |
1310 | } | |
1311 | } | |
1312 | else { | |
1313 | presult = new TCanvas; | |
1314 | pres = new TCanvas; | |
1315 | ppen = new TCanvas; | |
1316 | } | |
1317 | ||
1318 | ||
1319 | if (fgMaxInput == -1) | |
1320 | { | |
1321 | Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution"); | |
1322 | fgMaxInput = measured->GetNbinsX(); | |
1323 | } | |
1324 | if (fgMaxParams == -1) | |
1325 | { | |
1326 | Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution"); | |
1327 | fgMaxParams = initialConditions->GetNbinsX(); | |
1328 | } | |
1329 | ||
1330 | if (fgOverflowBinLimit > 0) | |
1331 | CreateOverflowBin(correlation, measured); | |
1332 | ||
1333 | // Uses Minuit methods | |
1334 | ||
1335 | SetStaticVariables(correlation, measured, efficiency); | |
1336 | ||
1337 | Double_t* params = new Double_t[fgMaxParams+1]; | |
1338 | for (Int_t i=0; i<fgMaxParams; ++i) | |
1339 | { | |
1340 | params[i] = initialConditions->GetBinContent(i+1) * efficiency->GetBinContent(i+1); | |
1341 | ||
1342 | Bool_t fix = kFALSE; | |
1343 | if (params[i] < 0) | |
1344 | { | |
1345 | fix = kTRUE; | |
1346 | params[i] = -params[i]; | |
1347 | } | |
1348 | params[i]=TMath::Sqrt(params[i]); | |
9e065ad2 | 1349 | } |
1350 | ||
1351 | Double_t chi2; | |
1352 | Int_t dummy; | |
1353 | ||
1354 | //fgPrintChi2Details = kTRUE; | |
1355 | fgCallCount = 0; // To make sure that Chi2Function prints the components | |
1356 | Chi2Function(dummy, 0, chi2, params, 0); | |
1357 | ||
1358 | presult->cd(); | |
03650114 | 1359 | TH1 *meas2 = (TH1*)measured->Clone("meas2"); |
1360 | meas2->SetTitle("meas2"); | |
1361 | meas2->SetName("meas2"); | |
9e065ad2 | 1362 | meas2->SetLineColor(1); |
1363 | meas2->SetLineWidth(2); | |
1364 | meas2->SetMarkerColor(meas2->GetLineColor()); | |
1365 | meas2->SetMarkerStyle(20); | |
1366 | Float_t scale = unfolded->GetXaxis()->GetBinWidth(1)/meas2->GetXaxis()->GetBinWidth(1); | |
1367 | meas2->Scale(scale); | |
1368 | if (blankCanvas) | |
1369 | meas2->DrawCopy(); | |
1370 | else | |
1371 | meas2->DrawCopy("same"); | |
1372 | //meas2->GetXaxis()->SetLimits(0,fgMaxInput); | |
1373 | meas2->SetBit(kCannotPick); | |
1374 | DrawGuess(params, presult, pres, ppen, reuseHists,unfolded); | |
2d99332c | 1375 | delete [] params; |
9e065ad2 | 1376 | } |
1377 | //____________________________________________________________________ | |
1378 | void AliUnfolding::RedrawInteractive() { | |
1379 | // | |
1380 | // Helper function for interactive unfolding | |
1381 | // | |
1382 | DrawResults(fghCorrelation,fghEfficiency,fghMeasured,fghUnfolded,fgCanvas,1,fghUnfolded); | |
1383 | } | |
1384 | //____________________________________________________________________ | |
1385 | void AliUnfolding::InteractiveUnfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions) | |
1386 | { | |
1387 | // | |
1388 | // Function to do interactive unfolding | |
1389 | // A canvas is drawn with the unfolding result | |
1390 | // Change the histogram with your mouse and all histograms | |
1391 | // will be updated automatically | |
1392 | ||
1393 | fgCanvas = new TCanvas("UnfoldingCanvas","Interactive unfolding",500,800); | |
1394 | fgCanvas->cd(); | |
1395 | ||
1396 | MakePads(); | |
1397 | ||
1398 | if (fghUnfolded) { | |
1399 | delete fghUnfolded; fghUnfolded = 0; | |
1400 | } | |
1401 | ||
1402 | gROOT->SetEditHistograms(kTRUE); | |
1403 | ||
03650114 | 1404 | fghUnfolded = new TH1F("AliUnfoldingfghUnfolded","Unfolded result (Interactive unfolder",efficiency->GetNbinsX(),efficiency->GetXaxis()->GetXmin(),efficiency->GetXaxis()->GetXmax()); |
9e065ad2 | 1405 | |
1406 | fghCorrelation = correlation; | |
1407 | fghEfficiency = efficiency; | |
1408 | fghMeasured = measured; | |
1409 | ||
03650114 | 1410 | if(fgMinuitMaxIterations>0) |
1411 | UnfoldWithMinuit(correlation,efficiency,measured,initialConditions,fghUnfolded,kFALSE); | |
1412 | else | |
1413 | fghUnfolded = initialConditions; | |
1414 | ||
1415 | fghUnfolded->SetLineColor(2); | |
1416 | fghUnfolded->SetMarkerColor(2); | |
1417 | fghUnfolded->SetLineWidth(2); | |
1418 | ||
9e065ad2 | 1419 | |
1420 | fgCanvas->cd(1); | |
1421 | fghUnfolded->Draw(); | |
1422 | DrawResults(correlation,efficiency,measured,fghUnfolded,fgCanvas,kFALSE,fghUnfolded); | |
1423 | ||
1424 | TExec *execRedraw = new TExec("redraw","AliUnfolding::RedrawInteractive()"); | |
1425 | fghUnfolded->GetListOfFunctions()->Add(execRedraw); | |
1426 | } | |
1427 | //____________________________________________________________________ | |
1428 | void AliUnfolding::DrawGuess(Double_t *params, TVirtualPad *pfolded, TVirtualPad *pres, TVirtualPad *ppen, Int_t reuseHists,TH1* unfolded) | |
19442b86 | 1429 | { |
1430 | // | |
1431 | // draws residuals of solution suggested by params and effect of regularization | |
1432 | // | |
1433 | ||
9e065ad2 | 1434 | if (pfolded == 0) |
1435 | pfolded = new TCanvas; | |
1436 | if (ppen == 0) | |
1437 | ppen = new TCanvas; | |
1438 | if (pres == 0) | |
1439 | pres = new TCanvas; | |
1440 | ||
19442b86 | 1441 | // d |
95e970ca | 1442 | TVectorD paramsVector(fgMaxParams); |
19442b86 | 1443 | for (Int_t i=0; i<fgMaxParams; ++i) |
1444 | paramsVector[i] = params[i] * params[i]; | |
1445 | ||
1446 | // Ad | |
1447 | TVectorD measGuessVector(fgMaxInput); | |
1448 | measGuessVector = (*fgCorrelationMatrix) * paramsVector; | |
1449 | ||
9e065ad2 | 1450 | TH1* folded = 0; |
1451 | if (reuseHists) | |
1452 | folded = dynamic_cast<TH1*>(gROOT->FindObject("__hfolded")); | |
1453 | if (!reuseHists || folded == 0) { | |
1454 | if (fgMeasuredAxis->GetXbins()->GetArray()) // variable bins | |
1455 | folded = new TH1F("__hfolded","Folded histo from AliUnfolding",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXbins()->GetArray()); | |
1456 | else | |
1457 | folded = new TH1F("__hfolded","Folded histo from AliUnfolding",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXmin(),fgMeasuredAxis->GetXmax()); | |
1458 | } | |
1459 | ||
1460 | folded->SetBit(kCannotPick); | |
1461 | folded->SetLineColor(4); | |
1462 | folded->SetLineWidth(2); | |
1463 | ||
1464 | for (Int_t ibin =0; ibin < fgMaxInput; ibin++) | |
1465 | folded->SetBinContent(ibin+1, measGuessVector[ibin]); | |
1466 | ||
1467 | Float_t scale = unfolded->GetXaxis()->GetBinWidth(1)/folded->GetXaxis()->GetBinWidth(1); | |
1468 | folded->Scale(scale); | |
1469 | ||
1470 | pfolded->cd(); | |
1471 | ||
1472 | folded->Draw("same"); | |
1473 | ||
19442b86 | 1474 | // Ad - m |
1475 | measGuessVector -= (*fgCurrentESDVector); | |
1476 | ||
1477 | TVectorD copy(measGuessVector); | |
19442b86 | 1478 | |
1479 | // (Ad - m) W | |
1480 | // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used | |
1481 | // normal way is like this: | |
1482 | // measGuessVector *= (*fgCorrelationCovarianceMatrix); | |
1483 | // optimized way like this: | |
1484 | for (Int_t i=0; i<fgMaxInput; ++i) | |
1485 | measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i); | |
1486 | ||
1487 | // (Ad - m) W (Ad - m) | |
1488 | // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very | |
1489 | // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit) | |
1490 | //Double_t chi2FromFit = measGuessVector * copy * 1e6; | |
1491 | ||
1492 | // draw residuals | |
9e065ad2 | 1493 | // Double_t pTarray[fgMaxParams+1]; |
1494 | // for(int i=0; i<fgMaxParams; i++) { | |
1495 | // pTarray[i] = fgUnfoldedAxis->GetBinCenter(i)-0.5*fgUnfoldedAxis->GetBinWidth(i); | |
1496 | // } | |
1497 | // pTarray[fgMaxParams] = fgUnfoldedAxis->GetBinCenter(fgMaxParams-1)+0.5*fgUnfoldedAxis->GetBinWidth(fgMaxParams-1); | |
1498 | // TH1* residuals = new TH1F("residuals", "residuals", fgMaxParams,pTarray); | |
1499 | // TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5); | |
1500 | // for (Int_t i=0; i<fgMaxInput; ++i) | |
1501 | // residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);7 | |
1502 | TH1* residuals = GetResidualsPlot(params); | |
1503 | residuals->GetXaxis()->SetTitleSize(0.06); | |
1504 | residuals->GetXaxis()->SetTitleOffset(0.7); | |
1505 | residuals->GetXaxis()->SetLabelSize(0.07); | |
1506 | residuals->GetYaxis()->SetTitleSize(0.08); | |
1507 | residuals->GetYaxis()->SetTitleOffset(0.5); | |
1508 | residuals->GetYaxis()->SetLabelSize(0.06); | |
1509 | pres->cd(); residuals->DrawCopy(); gPad->SetLogy(); | |
1510 | ||
19442b86 | 1511 | |
1512 | // draw penalty | |
95e970ca | 1513 | TH1* penalty = GetPenaltyPlot(params); |
9e065ad2 | 1514 | penalty->GetXaxis()->SetTitleSize(0.06); |
1515 | penalty->GetXaxis()->SetTitleOffset(0.7); | |
1516 | penalty->GetXaxis()->SetLabelSize(0.07); | |
1517 | penalty->GetYaxis()->SetTitleSize(0.08); | |
1518 | penalty->GetYaxis()->SetTitleOffset(0.5); | |
1519 | penalty->GetYaxis()->SetLabelSize(0.06); | |
1520 | ppen->cd(); penalty->DrawCopy(); gPad->SetLogy(); | |
1521 | } | |
1522 | ||
1523 | //____________________________________________________________________ | |
1524 | TH1* AliUnfolding::GetResidualsPlot(TH1* corrected) | |
1525 | { | |
1526 | // | |
1527 | // MvL: THIS MUST BE INCORRECT. | |
1528 | // Need heff to calculate params from TH1 'corrected' | |
1529 | // | |
1530 | ||
1531 | // | |
1532 | // fill residuals histogram of solution suggested by params and effect of regularization | |
1533 | // | |
1534 | ||
1535 | Double_t* params = new Double_t[fgMaxParams]; | |
b9b1eb8a | 1536 | for (Int_t i=0; i<fgMaxParams; i++) |
1537 | params[i] = 0; | |
1538 | ||
9e065ad2 | 1539 | for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++) |
3a44f1b9 | 1540 | params[i] = TMath::Sqrt(TMath::Abs(corrected->GetBinContent(i+1)*(*fgEfficiency)(i))); |
9e065ad2 | 1541 | |
2d99332c | 1542 | TH1 * plot = GetResidualsPlot(params); |
1543 | delete [] params; | |
1544 | return plot; | |
9e065ad2 | 1545 | } |
1546 | ||
1547 | //____________________________________________________________________ | |
1548 | TH1* AliUnfolding::GetResidualsPlot(Double_t* params) | |
1549 | { | |
1550 | // | |
1551 | // fill residuals histogram of solution suggested by params and effect of regularization | |
1552 | // | |
1553 | ||
1554 | // d | |
1555 | TVectorD paramsVector(fgMaxParams); | |
1556 | for (Int_t i=0; i<fgMaxParams; ++i) | |
1557 | paramsVector[i] = params[i] * params[i]; | |
1558 | ||
1559 | // Ad | |
1560 | TVectorD measGuessVector(fgMaxInput); | |
1561 | measGuessVector = (*fgCorrelationMatrix) * paramsVector; | |
1562 | ||
1563 | // Ad - m | |
1564 | measGuessVector -= (*fgCurrentESDVector); | |
1565 | ||
1566 | TVectorD copy(measGuessVector); | |
1567 | ||
1568 | // (Ad - m) W | |
1569 | // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used | |
1570 | // normal way is like this: | |
1571 | // measGuessVector *= (*fgCorrelationCovarianceMatrix); | |
1572 | // optimized way like this: | |
1573 | for (Int_t i=0; i<fgMaxInput; ++i) | |
1574 | measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i); | |
1575 | ||
1576 | // (Ad - m) W (Ad - m) | |
1577 | // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very | |
1578 | // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit) | |
1579 | //Double_t chi2FromFit = measGuessVector * copy * 1e6; | |
1580 | ||
1581 | // draw residuals | |
1582 | TH1* residuals = 0; | |
1583 | if (fgMeasuredAxis->GetXbins()->GetArray()) // variable bins | |
1584 | residuals = new TH1F("residuals", "residuals;unfolded pos;residual",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXbins()->GetArray()); | |
1585 | else | |
1586 | residuals = new TH1F("residuals", "residuals;unfolded pos;residual",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXmin(), fgMeasuredAxis->GetXmax()); | |
1587 | // TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5); | |
1588 | ||
1589 | Double_t sumResiduals = 0.; | |
1590 | for (Int_t i=0; i<fgMaxInput; ++i) { | |
1591 | residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6); | |
1592 | sumResiduals += measGuessVector(i) * copy(i) * 1e6; | |
1593 | } | |
1594 | fAvgResidual = sumResiduals/(double)fgMaxInput; | |
1595 | ||
1596 | return residuals; | |
95e970ca | 1597 | } |
19442b86 | 1598 | |
95e970ca | 1599 | //____________________________________________________________________ |
1600 | TH1* AliUnfolding::GetPenaltyPlot(TH1* corrected) | |
1601 | { | |
1602 | // draws the penalty factors as function of multiplicity of the current selected regularization | |
19442b86 | 1603 | |
95e970ca | 1604 | Double_t* params = new Double_t[fgMaxParams]; |
b9b1eb8a | 1605 | for (Int_t i=0; i<fgMaxParams; i++) |
1606 | params[i] = 0; | |
1607 | ||
95e970ca | 1608 | for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++) |
9e065ad2 | 1609 | params[i] = (*fgEfficiency)(i)*corrected->GetBinContent(i+1); |
95e970ca | 1610 | |
1611 | TH1* penalty = GetPenaltyPlot(params); | |
1612 | ||
1613 | delete[] params; | |
1614 | ||
1615 | return penalty; | |
1616 | } | |
1617 | ||
1618 | //____________________________________________________________________ | |
1619 | TH1* AliUnfolding::GetPenaltyPlot(Double_t* params) | |
1620 | { | |
1621 | // draws the penalty factors as function of multiplicity of the current selected regularization | |
1622 | ||
9e065ad2 | 1623 | //TH1* penalty = new TH1F("penalty", ";unfolded multiplicity;penalty factor", fgMaxParams, -0.5, fgMaxParams - 0.5); |
1624 | // TH1* penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgMaxParams, fgUnfoldedAxis->GetBinCenter(0)-0.5*fgUnfoldedAxis->GetBinWidth(0),fgUnfoldedAxis->GetBinCenter(fgMaxParams)+0.5*fgUnfoldedAxis->GetBinWidth(fgMaxParams) ); | |
1625 | ||
1626 | TH1* penalty = 0; | |
1627 | if (fgUnfoldedAxis->GetXbins()->GetArray()) | |
1628 | penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgUnfoldedAxis->GetNbins(),fgUnfoldedAxis->GetXbins()->GetArray()); | |
1629 | else | |
1630 | penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgUnfoldedAxis->GetNbins(),fgUnfoldedAxis->GetXmin(),fgUnfoldedAxis->GetXmax()); | |
95e970ca | 1631 | |
1632 | for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i) | |
19442b86 | 1633 | { |
95e970ca | 1634 | Double_t diff = 0; |
1635 | if (fgRegularizationType == kPol0) | |
1636 | { | |
9e065ad2 | 1637 | Double_t right = params[i] / fgUnfoldedAxis->GetBinWidth(i+1); |
1638 | Double_t left = params[i-1] / fgUnfoldedAxis->GetBinWidth(i); | |
19442b86 | 1639 | |
95e970ca | 1640 | if (left != 0) |
1641 | { | |
1642 | Double_t diffTmp = (right - left); | |
9e065ad2 | 1643 | diff = diffTmp * diffTmp / left / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2) / 100; |
95e970ca | 1644 | } |
1645 | } | |
1646 | if (fgRegularizationType == kPol1 && i > 1+fgSkipBinsBegin) | |
1647 | { | |
1648 | if (params[i-1] == 0) | |
1649 | continue; | |
19442b86 | 1650 | |
9e065ad2 | 1651 | Double_t right = params[i] / fgUnfoldedAxis->GetBinWidth(i+1); |
1652 | Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i); | |
1653 | Double_t left = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1); | |
19442b86 | 1654 | |
9e065ad2 | 1655 | Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2); |
1656 | Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2); | |
95e970ca | 1657 | |
1658 | diff = (der1 - der2) * (der1 - der2) / middle; | |
1659 | } | |
9e065ad2 | 1660 | |
95e970ca | 1661 | if (fgRegularizationType == kLog && i > 1+fgSkipBinsBegin) |
1662 | { | |
1663 | if (params[i-1] == 0) | |
1664 | continue; | |
1665 | ||
1666 | Double_t right = log(params[i]); | |
1667 | Double_t middle = log(params[i-1]); | |
1668 | Double_t left = log(params[i-2]); | |
1669 | ||
1670 | Double_t der1 = (right - middle); | |
1671 | Double_t der2 = (middle - left); | |
1672 | ||
1673 | //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2]; | |
1674 | //Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error); | |
1675 | ||
1676 | diff = (der1 - der2) * (der1 - der2);// / error; | |
1677 | } | |
9e065ad2 | 1678 | if (fgRegularizationType == kCurvature && i > 1+fgSkipBinsBegin) |
1679 | { | |
1680 | Double_t right = params[i]; // params are sqrt | |
1681 | Double_t middle = params[i-1]; | |
1682 | Double_t left = params[i-2]; | |
1683 | ||
1684 | Double_t der1 = (right - middle)/0.5/(fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i)); | |
1685 | Double_t der2 = (middle - left)/0.5/(fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i+1)); | |
1686 | ||
1687 | diff = (der1 - der2)/(fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1))*3.0; | |
1688 | diff = 1e4*diff*diff; | |
1689 | } | |
1690 | if (fgRegularizationType == kPowerLaw && i > 1+fgSkipBinsBegin) | |
1691 | { | |
19442b86 | 1692 | |
9e065ad2 | 1693 | if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8) |
1694 | continue; | |
1695 | ||
1696 | if (fgUnfoldedAxis->GetBinWidth(i+1) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8) | |
1697 | continue; | |
1698 | ||
1699 | double middle = TMath::Power(params[i-1] / fgUnfoldedAxis->GetBinWidth(i),fgPowern); | |
1700 | ||
1701 | if(middle>0) { | |
1702 | double right = TMath::Power(params[i] / fgUnfoldedAxis->GetBinWidth(i+1),fgPowern)/middle; | |
1703 | ||
1704 | double left = TMath::Power(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1),fgPowern)/middle; | |
1705 | ||
1706 | middle = 1.; | |
1707 | ||
1708 | Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2); | |
1709 | Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2); | |
1710 | ||
1711 | diff = (der1 - der2) * (der1 - der2);// / error; | |
1712 | } | |
1713 | } | |
1714 | ||
1715 | if (fgRegularizationType == kLogLog && i > 1+fgSkipBinsBegin) | |
1716 | { | |
1717 | ||
1718 | if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8) | |
1719 | continue; | |
1720 | ||
1721 | Double_t right = log(params[i] / (*fgEfficiency)(i) / fgUnfoldedAxis->GetBinWidth(i+1)); | |
1722 | Double_t middle = log(params[i-1] / (*fgEfficiency)(i-1) / fgUnfoldedAxis->GetBinWidth(i)); | |
1723 | Double_t left = log(params[i-2] / (*fgEfficiency)(i-2) / fgUnfoldedAxis->GetBinWidth(i-1)); | |
1724 | ||
1725 | Double_t der1 = (right - middle) / ( log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i)) ); | |
1726 | Double_t der2 = (middle - left) /( log(fgUnfoldedAxis->GetBinCenter(i)) - log(fgUnfoldedAxis->GetBinCenter(i-1)) ); | |
1727 | ||
1728 | double tmp = (log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i-1)))/2.; | |
1729 | Double_t dder = (der1-der2) / tmp; | |
1730 | ||
1731 | diff = dder * dder; | |
1732 | } | |
1733 | ||
1734 | penalty->SetBinContent(i, diff*fgRegularizationWeight); | |
19442b86 | 1735 | |
1736 | //Printf("%d %f %f %f %f", i-1, left, middle, right, diff); | |
1737 | } | |
95e970ca | 1738 | |
1739 | return penalty; | |
19442b86 | 1740 | } |
1741 | ||
1742 | //____________________________________________________________________ | |
1743 | void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3) | |
1744 | { | |
1745 | // | |
1746 | // fit function for minuit | |
1747 | // uses the TF1 stored in fgFitFunction | |
1748 | // | |
1749 | ||
1750 | for (Int_t i=0; i<fgFitFunction->GetNpar(); i++) | |
1751 | fgFitFunction->SetParameter(i, params[i]); | |
1752 | ||
1753 | Double_t* params2 = new Double_t[fgMaxParams]; | |
1754 | ||
1755 | for (Int_t i=0; i<fgMaxParams; ++i) | |
1756 | params2[i] = fgFitFunction->Eval(i); | |
1757 | ||
1758 | Chi2Function(unused1, unused2, chi2, params2, unused3); | |
1759 | ||
1760 | delete[] params2; | |
1761 | ||
1762 | if (fgDebug) | |
1763 | Printf("%f", chi2); | |
1764 | } | |
1765 | ||
1766 | //____________________________________________________________________ | |
1767 | Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult) | |
1768 | { | |
1769 | // | |
1770 | // correct spectrum using minuit chi2 method applying a functional fit | |
1771 | // | |
1772 | ||
1773 | if (!fgFitFunction) | |
1774 | { | |
1775 | Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting."); | |
1776 | return -1; | |
1777 | } | |
1778 | ||
1779 | SetChi2Regularization(kNone, 0); | |
1780 | ||
95e970ca | 1781 | SetStaticVariables(correlation, measured, efficiency); |
19442b86 | 1782 | |
1783 | // Initialize TMinuit via generic fitter interface | |
1784 | TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar()); | |
1785 | ||
1786 | minuit->SetFCN(TF1Function); | |
1787 | for (Int_t i=0; i<fgFitFunction->GetNpar(); i++) | |
1788 | { | |
1789 | Double_t lower, upper; | |
1790 | fgFitFunction->GetParLimits(i, lower, upper); | |
1791 | minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper); | |
1792 | } | |
1793 | ||
1794 | Double_t arglist[100]; | |
1795 | arglist[0] = 0; | |
1796 | minuit->ExecuteCommand("SET PRINT", arglist, 1); | |
9e065ad2 | 1797 | minuit->ExecuteCommand("SCAN", arglist, 0); |
19442b86 | 1798 | minuit->ExecuteCommand("MIGRAD", arglist, 0); |
1799 | //minuit->ExecuteCommand("MINOS", arglist, 0); | |
1800 | ||
1801 | for (Int_t i=0; i<fgFitFunction->GetNpar(); i++) | |
1802 | fgFitFunction->SetParameter(i, minuit->GetParameter(i)); | |
1803 | ||
1804 | for (Int_t i=0; i<fgMaxParams; ++i) | |
1805 | { | |
1806 | Double_t value = fgFitFunction->Eval(i); | |
1807 | if (fgDebug) | |
1808 | Printf("%d : %f", i, value); | |
1809 | ||
1810 | if (efficiency) | |
1811 | { | |
1812 | if (efficiency->GetBinContent(i+1) > 0) | |
1813 | { | |
1814 | value /= efficiency->GetBinContent(i+1); | |
1815 | } | |
1816 | else | |
1817 | value = 0; | |
1818 | } | |
1819 | aResult->SetBinContent(i+1, value); | |
1820 | aResult->SetBinError(i+1, 0); | |
1821 | } | |
1822 | ||
1823 | return 0; | |
1824 | } | |
1825 | ||
1826 | //____________________________________________________________________ | |
1827 | void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured) | |
1828 | { | |
1829 | // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins | |
1830 | // The same limit is applied to the correlation | |
1831 | ||
1832 | Int_t lastBin = 0; | |
1833 | for (Int_t i=1; i<=measured->GetNbinsX(); ++i) | |
1834 | { | |
1835 | if (measured->GetBinContent(i) <= fgOverflowBinLimit) | |
1836 | { | |
1837 | lastBin = i; | |
1838 | break; | |
1839 | } | |
1840 | } | |
1841 | ||
1842 | if (lastBin == 0) | |
1843 | { | |
1844 | Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit); | |
1845 | return; | |
1846 | } | |
1847 | ||
1848 | Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin); | |
1849 | ||
1850 | TCanvas* canvas = 0; | |
1851 | ||
1852 | if (fgDebug) | |
1853 | { | |
1854 | canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800); | |
1855 | canvas->Divide(2, 2); | |
1856 | ||
1857 | canvas->cd(1); | |
1858 | measured->SetStats(kFALSE); | |
1859 | measured->DrawCopy(); | |
1860 | gPad->SetLogy(); | |
1861 | ||
1862 | canvas->cd(2); | |
1863 | correlation->SetStats(kFALSE); | |
1864 | correlation->DrawCopy("COLZ"); | |
1865 | } | |
1866 | ||
1867 | measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX())); | |
1868 | for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i) | |
1869 | { | |
1870 | measured->SetBinContent(i, 0); | |
1871 | measured->SetBinError(i, 0); | |
1872 | } | |
1873 | // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? | |
1874 | measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin))); | |
1875 | ||
1876 | Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin)); | |
1877 | ||
1878 | for (Int_t i=1; i<=correlation->GetNbinsX(); ++i) | |
1879 | { | |
1880 | correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY())); | |
1881 | // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? | |
1882 | correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin))); | |
1883 | ||
1884 | for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j) | |
1885 | { | |
1886 | correlation->SetBinContent(i, j, 0); | |
1887 | correlation->SetBinError(i, j, 0); | |
1888 | } | |
1889 | } | |
1890 | ||
1891 | Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!"); | |
1892 | ||
1893 | if (canvas) | |
1894 | { | |
1895 | canvas->cd(3); | |
1896 | measured->DrawCopy(); | |
1897 | gPad->SetLogy(); | |
1898 | ||
1899 | canvas->cd(4); | |
1900 | correlation->DrawCopy("COLZ"); | |
1901 | } | |
1902 | } | |
95e970ca | 1903 | |
1904 | Int_t AliUnfolding::UnfoldGetBias(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result) | |
1905 | { | |
1906 | // unfolds and assigns bias as errors with Eq. 19 of Cowan, "a survey of unfolding methods for particle physics" | |
1907 | // b_i = sum_j dmu_i/dn_j (nu_j - n_j) with nu_j as folded guess, n_j as data | |
1908 | // dmu_i/dn_j is found numerically by changing the bin content and re-unfolding | |
1909 | // | |
1910 | // return code: 0 (success) -1 (error: from Unfold(...) ) | |
1911 | ||
1912 | if (Unfold(correlation, efficiency, measured, initialConditions, result) != 0) | |
1913 | return -1; | |
1914 | ||
1915 | TMatrixD matrix(fgMaxInput, fgMaxParams); | |
1916 | ||
1917 | TH1* newResult[4]; | |
1918 | for (Int_t loop=0; loop<4; loop++) | |
1919 | newResult[loop] = (TH1*) result->Clone(Form("newresult_%d", loop)); | |
1920 | ||
1921 | // change bin-by-bin and built matrix of effects | |
1922 | for (Int_t m=0; m<fgMaxInput; m++) | |
1923 | { | |
1924 | if (measured->GetBinContent(m+1) < 1) | |
1925 | continue; | |
1926 | ||
1927 | for (Int_t loop=0; loop<4; loop++) | |
1928 | { | |
1929 | //Printf("%d %d", i, loop); | |
1930 | Float_t factor = 1; | |
1931 | switch (loop) | |
1932 | { | |
1933 | case 0: factor = 0.5; break; | |
1934 | case 1: factor = -0.5; break; | |
1935 | case 2: factor = 1; break; | |
1936 | case 3: factor = -1; break; | |
1937 | default: return -1; | |
1938 | } | |
1939 | ||
1940 | TH1* measuredClone = (TH1*) measured->Clone("measuredClone"); | |
1941 | ||
1942 | measuredClone->SetBinContent(m+1, measured->GetBinContent(m+1) + factor * TMath::Sqrt(measured->GetBinContent(m+1))); | |
1943 | //new TCanvas; measuredClone->Draw("COLZ"); | |
1944 | ||
1945 | newResult[loop]->Reset(); | |
1946 | Unfold(correlation, efficiency, measuredClone, measuredClone, newResult[loop]); | |
1947 | // WARNING if we leave here, then nothing is calculated | |
1948 | //return -1; | |
1949 | ||
1950 | delete measuredClone; | |
1951 | } | |
1952 | ||
1953 | for (Int_t t=0; t<fgMaxParams; t++) | |
1954 | { | |
1955 | // one value | |
1956 | //matrix(i, j) = (result->GetBinContent(j+1) - newResult->GetBinContent(j+1)) / TMath::Sqrt(changedHist->GetBinContent(1, i+1)); | |
1957 | ||
1958 | // four values from the derivate (procedure taken from ROOT -- suggestion by Ruben) | |
1959 | // = 1/2D * [ 8 (f(D/2) - f(-D/2)) - (f(D)-f(-D)) ]/3 | |
1960 | ||
1961 | /* | |
1962 | // check formula | |
1963 | measured->SetBinContent(1, m+1, 100); | |
1964 | newResult[0]->SetBinContent(t+1, 5 * 0.5 + 10); | |
1965 | newResult[1]->SetBinContent(t+1, 5 * -0.5 + 10); | |
1966 | newResult[2]->SetBinContent(t+1, 5 * 1 + 10); | |
1967 | newResult[3]->SetBinContent(t+1, 5 * -1 + 10); | |
1968 | */ | |
1969 | ||
1970 | matrix(m, t) = 0.5 / TMath::Sqrt(measured->GetBinContent(m+1)) * | |
1971 | ( 8. * (newResult[0]->GetBinContent(t+1) - newResult[1]->GetBinContent(t+1)) - | |
1972 | (newResult[2]->GetBinContent(t+1) - newResult[3]->GetBinContent(t+1)) | |
1973 | ) / 3; | |
1974 | } | |
1975 | } | |
1976 | ||
1977 | for (Int_t loop=0; loop<4; loop++) | |
1978 | delete newResult[loop]; | |
1979 | ||
95e970ca | 1980 | // ...calculate folded guess |
1981 | TH1* convoluted = (TH1*) measured->Clone("convoluted"); | |
1982 | convoluted->Reset(); | |
1983 | for (Int_t m=0; m<fgMaxInput; m++) | |
1984 | { | |
1985 | Float_t value = 0; | |
1986 | for (Int_t t = 0; t<fgMaxParams; t++) | |
1987 | { | |
1988 | Float_t tmp = correlation->GetBinContent(t+1, m+1) * result->GetBinContent(t+1); | |
1989 | if (efficiency) | |
1990 | tmp *= efficiency->GetBinContent(t+1); | |
1991 | value += tmp; | |
1992 | } | |
1993 | convoluted->SetBinContent(m+1, value); | |
1994 | } | |
1995 | ||
1996 | // rescale | |
1997 | convoluted->Scale(measured->Integral() / convoluted->Integral()); | |
1998 | ||
1999 | //new TCanvas; convoluted->Draw(); measured->Draw("SAME"); measured->SetLineColor(2); | |
2000 | //return; | |
2001 | ||
2002 | // difference | |
2003 | convoluted->Add(measured, -1); | |
2004 | ||
2005 | // ...and the bias | |
2006 | for (Int_t t = 0; t<fgMaxParams; t++) | |
2007 | { | |
2008 | Double_t error = 0; | |
2009 | for (Int_t m=0; m<fgMaxInput; m++) | |
2010 | error += matrix(m, t) * convoluted->GetBinContent(m+1); | |
2011 | result->SetBinError(t+1, error); | |
2012 | } | |
2013 | ||
2014 | //new TCanvas; result->Draw(); gPad->SetLogy(); | |
2015 | ||
2016 | return 0; | |
2017 | } |