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