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