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