Fix compilation errors
[u/mrichter/AliRoot.git] / PWG0 / AliUnfolding.cxx
CommitLineData
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
31TMatrixD* AliUnfolding::fgCorrelationMatrix = 0;
95e970ca 32TMatrixD* AliUnfolding::fgCorrelationMatrixSquared = 0;
19442b86 33TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0;
34TVectorD* AliUnfolding::fgCurrentESDVector = 0;
35TVectorD* AliUnfolding::fgEntropyAPriori = 0;
95e970ca 36TVectorD* AliUnfolding::fgEfficiency = 0;
37TVectorD* AliUnfolding::fgBinWidths = 0;
19442b86 38
39TF1* AliUnfolding::fgFitFunction = 0;
40
41AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid;
42Int_t AliUnfolding::fgMaxInput = -1; // bins in measured histogram
43Int_t AliUnfolding::fgMaxParams = -1; // bins in unfolded histogram = number of fit params
44Float_t AliUnfolding::fgOverflowBinLimit = -1;
45
46AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1;
47Float_t AliUnfolding::fgRegularizationWeight = 10000;
48Int_t AliUnfolding::fgSkipBinsBegin = 0;
49Float_t AliUnfolding::fgMinuitStepSize = 0.1; // (usually not needed to be changed) step size in minimization
95e970ca 50Bool_t AliUnfolding::fgMinimumInitialValue = kFALSE; // set all initial values at least to the smallest value among the initial values
51Float_t AliUnfolding::fgMinimumInitialValueFix = -1;
19442b86 52Bool_t AliUnfolding::fgNormalizeInput = kFALSE; // normalize input spectrum
95e970ca 53Float_t AliUnfolding::fgNotFoundEvents = 0;
54Bool_t AliUnfolding::fgSkipBin0InChi2 = kFALSE;
19442b86 55
56Float_t AliUnfolding::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing)
651e2088 57Int_t AliUnfolding::fgBayesianIterations = 10; // number of iterations in Bayesian method
19442b86 58
59Bool_t AliUnfolding::fgDebug = kFALSE;
60
95e970ca 61Int_t AliUnfolding::fgCallCount = 0;
62
19442b86 63ClassImp(AliUnfolding)
64
65//____________________________________________________________________
66void 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//____________________________________________________________________
83void 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//____________________________________________________________________
93void 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//____________________________________________________________________
103void 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//____________________________________________________________________
149void 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//____________________________________________________________________
162void 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//____________________________________________________________________
175void 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//____________________________________________________________________
185Int_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 231void 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//____________________________________________________________________
319Int_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)
deaac8b1 373 Float_t minValue = 1e35;
95e970ca 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//____________________________________________________________________
483Int_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//____________________________________________________________________
809Double_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//____________________________________________________________________
833Double_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//____________________________________________________________________
862Double_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//____________________________________________________________________
893Double_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//____________________________________________________________________
920Double_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//____________________________________________________________________
949Double_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//____________________________________________________________________
979void 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//____________________________________________________________________
1107void 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//____________________________________________________________________
1154TH1* 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//____________________________________________________________________
1170TH1* 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//____________________________________________________________________
1231void 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//____________________________________________________________________
1255Int_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//____________________________________________________________________
1314void 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
1391Int_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}