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