1) UnloadFromCache does not only remove the entry (or the entries) from the cache...
[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];
1518 for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++)
3a44f1b9 1519 params[i] = TMath::Sqrt(TMath::Abs(corrected->GetBinContent(i+1)*(*fgEfficiency)(i)));
9e065ad2 1520
1521
2d99332c 1522 TH1 * plot = GetResidualsPlot(params);
1523 delete [] params;
1524 return plot;
9e065ad2 1525}
1526
1527//____________________________________________________________________
1528TH1* AliUnfolding::GetResidualsPlot(Double_t* params)
1529{
1530 //
1531 // fill residuals histogram of solution suggested by params and effect of regularization
1532 //
1533
1534 // d
1535 TVectorD paramsVector(fgMaxParams);
1536 for (Int_t i=0; i<fgMaxParams; ++i)
1537 paramsVector[i] = params[i] * params[i];
1538
1539 // Ad
1540 TVectorD measGuessVector(fgMaxInput);
1541 measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1542
1543 // Ad - m
1544 measGuessVector -= (*fgCurrentESDVector);
1545
1546 TVectorD copy(measGuessVector);
1547
1548 // (Ad - m) W
1549 // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1550 // normal way is like this:
1551 // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1552 // optimized way like this:
1553 for (Int_t i=0; i<fgMaxInput; ++i)
1554 measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1555
1556 // (Ad - m) W (Ad - m)
1557 // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1558 // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
1559 //Double_t chi2FromFit = measGuessVector * copy * 1e6;
1560
1561 // draw residuals
1562 TH1* residuals = 0;
1563 if (fgMeasuredAxis->GetXbins()->GetArray()) // variable bins
1564 residuals = new TH1F("residuals", "residuals;unfolded pos;residual",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXbins()->GetArray());
1565 else
1566 residuals = new TH1F("residuals", "residuals;unfolded pos;residual",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXmin(), fgMeasuredAxis->GetXmax());
1567 // TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
1568
1569 Double_t sumResiduals = 0.;
1570 for (Int_t i=0; i<fgMaxInput; ++i) {
1571 residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
1572 sumResiduals += measGuessVector(i) * copy(i) * 1e6;
1573 }
1574 fAvgResidual = sumResiduals/(double)fgMaxInput;
1575
1576 return residuals;
95e970ca 1577}
19442b86 1578
95e970ca 1579//____________________________________________________________________
1580TH1* AliUnfolding::GetPenaltyPlot(TH1* corrected)
1581{
1582 // draws the penalty factors as function of multiplicity of the current selected regularization
19442b86 1583
95e970ca 1584 Double_t* params = new Double_t[fgMaxParams];
1585 for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++)
9e065ad2 1586 params[i] = (*fgEfficiency)(i)*corrected->GetBinContent(i+1);
95e970ca 1587
1588 TH1* penalty = GetPenaltyPlot(params);
1589
1590 delete[] params;
1591
1592 return penalty;
1593}
1594
1595//____________________________________________________________________
1596TH1* AliUnfolding::GetPenaltyPlot(Double_t* params)
1597{
1598 // draws the penalty factors as function of multiplicity of the current selected regularization
1599
9e065ad2 1600 //TH1* penalty = new TH1F("penalty", ";unfolded multiplicity;penalty factor", fgMaxParams, -0.5, fgMaxParams - 0.5);
1601 // 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) );
1602
1603 TH1* penalty = 0;
1604 if (fgUnfoldedAxis->GetXbins()->GetArray())
1605 penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgUnfoldedAxis->GetNbins(),fgUnfoldedAxis->GetXbins()->GetArray());
1606 else
1607 penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgUnfoldedAxis->GetNbins(),fgUnfoldedAxis->GetXmin(),fgUnfoldedAxis->GetXmax());
95e970ca 1608
1609 for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
19442b86 1610 {
95e970ca 1611 Double_t diff = 0;
1612 if (fgRegularizationType == kPol0)
1613 {
9e065ad2 1614 Double_t right = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
1615 Double_t left = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
19442b86 1616
95e970ca 1617 if (left != 0)
1618 {
1619 Double_t diffTmp = (right - left);
9e065ad2 1620 diff = diffTmp * diffTmp / left / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2) / 100;
95e970ca 1621 }
1622 }
1623 if (fgRegularizationType == kPol1 && i > 1+fgSkipBinsBegin)
1624 {
1625 if (params[i-1] == 0)
1626 continue;
19442b86 1627
9e065ad2 1628 Double_t right = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
1629 Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
1630 Double_t left = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1);
19442b86 1631
9e065ad2 1632 Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
1633 Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
95e970ca 1634
1635 diff = (der1 - der2) * (der1 - der2) / middle;
1636 }
9e065ad2 1637
95e970ca 1638 if (fgRegularizationType == kLog && i > 1+fgSkipBinsBegin)
1639 {
1640 if (params[i-1] == 0)
1641 continue;
1642
1643 Double_t right = log(params[i]);
1644 Double_t middle = log(params[i-1]);
1645 Double_t left = log(params[i-2]);
1646
1647 Double_t der1 = (right - middle);
1648 Double_t der2 = (middle - left);
1649
1650 //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
1651 //Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
1652
1653 diff = (der1 - der2) * (der1 - der2);// / error;
1654 }
9e065ad2 1655 if (fgRegularizationType == kCurvature && i > 1+fgSkipBinsBegin)
1656 {
1657 Double_t right = params[i]; // params are sqrt
1658 Double_t middle = params[i-1];
1659 Double_t left = params[i-2];
1660
1661 Double_t der1 = (right - middle)/0.5/(fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i));
1662 Double_t der2 = (middle - left)/0.5/(fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i+1));
1663
1664 diff = (der1 - der2)/(fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1))*3.0;
1665 diff = 1e4*diff*diff;
1666 }
1667 if (fgRegularizationType == kPowerLaw && i > 1+fgSkipBinsBegin)
1668 {
19442b86 1669
9e065ad2 1670 if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
1671 continue;
1672
1673 if (fgUnfoldedAxis->GetBinWidth(i+1) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8)
1674 continue;
1675
1676 double middle = TMath::Power(params[i-1] / fgUnfoldedAxis->GetBinWidth(i),fgPowern);
1677
1678 if(middle>0) {
1679 double right = TMath::Power(params[i] / fgUnfoldedAxis->GetBinWidth(i+1),fgPowern)/middle;
1680
1681 double left = TMath::Power(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1),fgPowern)/middle;
1682
1683 middle = 1.;
1684
1685 Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
1686 Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
1687
1688 diff = (der1 - der2) * (der1 - der2);// / error;
1689 }
1690 }
1691
1692 if (fgRegularizationType == kLogLog && i > 1+fgSkipBinsBegin)
1693 {
1694
1695 if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
1696 continue;
1697
1698 Double_t right = log(params[i] / (*fgEfficiency)(i) / fgUnfoldedAxis->GetBinWidth(i+1));
1699 Double_t middle = log(params[i-1] / (*fgEfficiency)(i-1) / fgUnfoldedAxis->GetBinWidth(i));
1700 Double_t left = log(params[i-2] / (*fgEfficiency)(i-2) / fgUnfoldedAxis->GetBinWidth(i-1));
1701
1702 Double_t der1 = (right - middle) / ( log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i)) );
1703 Double_t der2 = (middle - left) /( log(fgUnfoldedAxis->GetBinCenter(i)) - log(fgUnfoldedAxis->GetBinCenter(i-1)) );
1704
1705 double tmp = (log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i-1)))/2.;
1706 Double_t dder = (der1-der2) / tmp;
1707
1708 diff = dder * dder;
1709 }
1710
1711 penalty->SetBinContent(i, diff*fgRegularizationWeight);
19442b86 1712
1713 //Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
1714 }
95e970ca 1715
1716 return penalty;
19442b86 1717}
1718
1719//____________________________________________________________________
1720void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
1721{
1722 //
1723 // fit function for minuit
1724 // uses the TF1 stored in fgFitFunction
1725 //
1726
1727 for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1728 fgFitFunction->SetParameter(i, params[i]);
1729
1730 Double_t* params2 = new Double_t[fgMaxParams];
1731
1732 for (Int_t i=0; i<fgMaxParams; ++i)
1733 params2[i] = fgFitFunction->Eval(i);
1734
1735 Chi2Function(unused1, unused2, chi2, params2, unused3);
1736
1737 delete[] params2;
1738
1739 if (fgDebug)
1740 Printf("%f", chi2);
1741}
1742
1743//____________________________________________________________________
1744Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult)
1745{
1746 //
1747 // correct spectrum using minuit chi2 method applying a functional fit
1748 //
1749
1750 if (!fgFitFunction)
1751 {
1752 Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting.");
1753 return -1;
1754 }
1755
1756 SetChi2Regularization(kNone, 0);
1757
95e970ca 1758 SetStaticVariables(correlation, measured, efficiency);
19442b86 1759
1760 // Initialize TMinuit via generic fitter interface
1761 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar());
1762
1763 minuit->SetFCN(TF1Function);
1764 for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1765 {
1766 Double_t lower, upper;
1767 fgFitFunction->GetParLimits(i, lower, upper);
1768 minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper);
1769 }
1770
1771 Double_t arglist[100];
1772 arglist[0] = 0;
1773 minuit->ExecuteCommand("SET PRINT", arglist, 1);
9e065ad2 1774 minuit->ExecuteCommand("SCAN", arglist, 0);
19442b86 1775 minuit->ExecuteCommand("MIGRAD", arglist, 0);
1776 //minuit->ExecuteCommand("MINOS", arglist, 0);
1777
1778 for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1779 fgFitFunction->SetParameter(i, minuit->GetParameter(i));
1780
1781 for (Int_t i=0; i<fgMaxParams; ++i)
1782 {
1783 Double_t value = fgFitFunction->Eval(i);
1784 if (fgDebug)
1785 Printf("%d : %f", i, value);
1786
1787 if (efficiency)
1788 {
1789 if (efficiency->GetBinContent(i+1) > 0)
1790 {
1791 value /= efficiency->GetBinContent(i+1);
1792 }
1793 else
1794 value = 0;
1795 }
1796 aResult->SetBinContent(i+1, value);
1797 aResult->SetBinError(i+1, 0);
1798 }
1799
1800 return 0;
1801}
1802
1803//____________________________________________________________________
1804void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured)
1805{
1806 // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins
1807 // The same limit is applied to the correlation
1808
1809 Int_t lastBin = 0;
1810 for (Int_t i=1; i<=measured->GetNbinsX(); ++i)
1811 {
1812 if (measured->GetBinContent(i) <= fgOverflowBinLimit)
1813 {
1814 lastBin = i;
1815 break;
1816 }
1817 }
1818
1819 if (lastBin == 0)
1820 {
1821 Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit);
1822 return;
1823 }
1824
1825 Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin);
1826
1827 TCanvas* canvas = 0;
1828
1829 if (fgDebug)
1830 {
1831 canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
1832 canvas->Divide(2, 2);
1833
1834 canvas->cd(1);
1835 measured->SetStats(kFALSE);
1836 measured->DrawCopy();
1837 gPad->SetLogy();
1838
1839 canvas->cd(2);
1840 correlation->SetStats(kFALSE);
1841 correlation->DrawCopy("COLZ");
1842 }
1843
1844 measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX()));
1845 for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i)
1846 {
1847 measured->SetBinContent(i, 0);
1848 measured->SetBinError(i, 0);
1849 }
1850 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1851 measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin)));
1852
1853 Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin));
1854
1855 for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
1856 {
1857 correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY()));
1858 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1859 correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin)));
1860
1861 for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j)
1862 {
1863 correlation->SetBinContent(i, j, 0);
1864 correlation->SetBinError(i, j, 0);
1865 }
1866 }
1867
1868 Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!");
1869
1870 if (canvas)
1871 {
1872 canvas->cd(3);
1873 measured->DrawCopy();
1874 gPad->SetLogy();
1875
1876 canvas->cd(4);
1877 correlation->DrawCopy("COLZ");
1878 }
1879}
95e970ca 1880
1881Int_t AliUnfolding::UnfoldGetBias(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result)
1882{
1883 // unfolds and assigns bias as errors with Eq. 19 of Cowan, "a survey of unfolding methods for particle physics"
1884 // b_i = sum_j dmu_i/dn_j (nu_j - n_j) with nu_j as folded guess, n_j as data
1885 // dmu_i/dn_j is found numerically by changing the bin content and re-unfolding
1886 //
1887 // return code: 0 (success) -1 (error: from Unfold(...) )
1888
1889 if (Unfold(correlation, efficiency, measured, initialConditions, result) != 0)
1890 return -1;
1891
1892 TMatrixD matrix(fgMaxInput, fgMaxParams);
1893
1894 TH1* newResult[4];
1895 for (Int_t loop=0; loop<4; loop++)
1896 newResult[loop] = (TH1*) result->Clone(Form("newresult_%d", loop));
1897
1898 // change bin-by-bin and built matrix of effects
1899 for (Int_t m=0; m<fgMaxInput; m++)
1900 {
1901 if (measured->GetBinContent(m+1) < 1)
1902 continue;
1903
1904 for (Int_t loop=0; loop<4; loop++)
1905 {
1906 //Printf("%d %d", i, loop);
1907 Float_t factor = 1;
1908 switch (loop)
1909 {
1910 case 0: factor = 0.5; break;
1911 case 1: factor = -0.5; break;
1912 case 2: factor = 1; break;
1913 case 3: factor = -1; break;
1914 default: return -1;
1915 }
1916
1917 TH1* measuredClone = (TH1*) measured->Clone("measuredClone");
1918
1919 measuredClone->SetBinContent(m+1, measured->GetBinContent(m+1) + factor * TMath::Sqrt(measured->GetBinContent(m+1)));
1920 //new TCanvas; measuredClone->Draw("COLZ");
1921
1922 newResult[loop]->Reset();
1923 Unfold(correlation, efficiency, measuredClone, measuredClone, newResult[loop]);
1924 // WARNING if we leave here, then nothing is calculated
1925 //return -1;
1926
1927 delete measuredClone;
1928 }
1929
1930 for (Int_t t=0; t<fgMaxParams; t++)
1931 {
1932 // one value
1933 //matrix(i, j) = (result->GetBinContent(j+1) - newResult->GetBinContent(j+1)) / TMath::Sqrt(changedHist->GetBinContent(1, i+1));
1934
1935 // four values from the derivate (procedure taken from ROOT -- suggestion by Ruben)
1936 // = 1/2D * [ 8 (f(D/2) - f(-D/2)) - (f(D)-f(-D)) ]/3
1937
1938 /*
1939 // check formula
1940 measured->SetBinContent(1, m+1, 100);
1941 newResult[0]->SetBinContent(t+1, 5 * 0.5 + 10);
1942 newResult[1]->SetBinContent(t+1, 5 * -0.5 + 10);
1943 newResult[2]->SetBinContent(t+1, 5 * 1 + 10);
1944 newResult[3]->SetBinContent(t+1, 5 * -1 + 10);
1945 */
1946
1947 matrix(m, t) = 0.5 / TMath::Sqrt(measured->GetBinContent(m+1)) *
1948 ( 8. * (newResult[0]->GetBinContent(t+1) - newResult[1]->GetBinContent(t+1)) -
1949 (newResult[2]->GetBinContent(t+1) - newResult[3]->GetBinContent(t+1))
1950 ) / 3;
1951 }
1952 }
1953
1954 for (Int_t loop=0; loop<4; loop++)
1955 delete newResult[loop];
1956
95e970ca 1957 // ...calculate folded guess
1958 TH1* convoluted = (TH1*) measured->Clone("convoluted");
1959 convoluted->Reset();
1960 for (Int_t m=0; m<fgMaxInput; m++)
1961 {
1962 Float_t value = 0;
1963 for (Int_t t = 0; t<fgMaxParams; t++)
1964 {
1965 Float_t tmp = correlation->GetBinContent(t+1, m+1) * result->GetBinContent(t+1);
1966 if (efficiency)
1967 tmp *= efficiency->GetBinContent(t+1);
1968 value += tmp;
1969 }
1970 convoluted->SetBinContent(m+1, value);
1971 }
1972
1973 // rescale
1974 convoluted->Scale(measured->Integral() / convoluted->Integral());
1975
1976 //new TCanvas; convoluted->Draw(); measured->Draw("SAME"); measured->SetLineColor(2);
1977 //return;
1978
1979 // difference
1980 convoluted->Add(measured, -1);
1981
1982 // ...and the bias
1983 for (Int_t t = 0; t<fgMaxParams; t++)
1984 {
1985 Double_t error = 0;
1986 for (Int_t m=0; m<fgMaxInput; m++)
1987 error += matrix(m, t) * convoluted->GetBinContent(m+1);
1988 result->SetBinError(t+1, error);
1989 }
1990
1991 //new TCanvas; result->Draw(); gPad->SetLogy();
1992
1993 return 0;
1994}