1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliUnfolding.cxx 31168 2009-02-23 15:18:45Z jgrosseo $ */
18 // This class allows 1-dimensional unfolding.
19 // Methods that are implemented are chi2 minimization and bayesian unfolding.
21 // Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
23 #include "AliUnfolding.h"
26 #include <TVirtualFitter.h>
31 TMatrixD* AliUnfolding::fgCorrelationMatrix = 0;
32 TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0;
33 TVectorD* AliUnfolding::fgCurrentESDVector = 0;
34 TVectorD* AliUnfolding::fgEntropyAPriori = 0;
36 TF1* AliUnfolding::fgFitFunction = 0;
38 AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid;
39 Int_t AliUnfolding::fgMaxInput = -1; // bins in measured histogram
40 Int_t AliUnfolding::fgMaxParams = -1; // bins in unfolded histogram = number of fit params
41 Float_t AliUnfolding::fgOverflowBinLimit = -1;
43 AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1;
44 Float_t AliUnfolding::fgRegularizationWeight = 10000;
45 Int_t AliUnfolding::fgSkipBinsBegin = 0;
46 Float_t AliUnfolding::fgMinuitStepSize = 0.1; // (usually not needed to be changed) step size in minimization
47 Bool_t AliUnfolding::fgNormalizeInput = kFALSE; // normalize input spectrum
49 Float_t AliUnfolding::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing)
50 Int_t AliUnfolding::fgBayesianIterations = 10; // number of iterations in Bayesian method
52 Bool_t AliUnfolding::fgDebug = kFALSE;
54 ClassImp(AliUnfolding)
56 //____________________________________________________________________
57 void AliUnfolding::SetUnfoldingMethod(MethodType methodType)
59 // set unfolding method
60 fgMethodType = methodType;
65 case kInvalid: name = "INVALID"; break;
66 case kChi2Minimization: name = "Chi2 Minimization"; break;
67 case kBayesian: name = "Bayesian unfolding"; break;
68 case kFunction: name = "Functional fit"; break;
70 Printf("AliUnfolding::SetUnfoldingMethod: %s enabled.", name);
73 //____________________________________________________________________
74 void AliUnfolding::SetCreateOverflowBin(Float_t overflowBinLimit)
76 // enable the creation of a overflow bin that includes all statistics below the given limit
78 fgOverflowBinLimit = overflowBinLimit;
80 Printf("AliUnfolding::SetCreateOverflowBin: overflow bin limit set to %f", overflowBinLimit);
83 //____________________________________________________________________
84 void AliUnfolding::SetSkipBinsBegin(Int_t nBins)
86 // set number of skipped bins in regularization
88 fgSkipBinsBegin = nBins;
90 Printf("AliUnfolding::SetSkipBinsBegin: skipping %d bins at the beginning of the spectrum in the regularization.", fgSkipBinsBegin);
93 //____________________________________________________________________
94 void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded)
96 // set number of bins in the input (measured) distribution and in the unfolded distribution
97 fgMaxInput = nMeasured;
98 fgMaxParams = nUnfolded;
100 Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded);
103 //____________________________________________________________________
104 void AliUnfolding::SetChi2Regularization(RegularizationType type, Float_t weight)
107 // sets the parameters for chi2 minimization
110 fgRegularizationType = type;
111 fgRegularizationWeight = weight;
113 Printf("AliUnfolding::SetChi2Regularization --> Regularization set to %d with weight %f", (Int_t) type, weight);
116 //____________________________________________________________________
117 void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
120 // sets the parameters for Bayesian unfolding
123 fgBayesianSmoothing = smoothing;
124 fgBayesianIterations = nIterations;
126 Printf("AliUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing);
129 //____________________________________________________________________
130 void AliUnfolding::SetFunction(TF1* function)
132 // set function for unfolding with a fit function
134 fgFitFunction = function;
136 Printf("AliUnfolding::SetFunction: Set fit function with %d parameters.", function->GetNpar());
139 //____________________________________________________________________
140 Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
142 // unfolds with unfolding method fgMethodType
145 // correlation: response matrix as measured vs. generated
146 // 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.
147 // measured: the measured spectrum
148 // initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions.
149 // result: target for the unfolded result
150 // check: depends on the unfolding method, see comments in specific functions
152 if (fgMaxInput == -1)
154 Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution");
155 fgMaxInput = measured->GetNbinsX();
157 if (fgMaxParams == -1)
159 Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution");
160 fgMaxParams = measured->GetNbinsX();
163 if (fgOverflowBinLimit > 0)
164 CreateOverflowBin(correlation, measured);
166 switch (fgMethodType)
170 Printf("AliUnfolding::Unfold: ERROR: Unfolding method not set. Use SetUnfoldingMethod. Exiting...");
173 case kChi2Minimization:
174 return UnfoldWithMinuit(correlation, efficiency, measured, initialConditions, result, check);
176 return UnfoldWithBayesian(correlation, efficiency, measured, initialConditions, result);
178 return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result);
183 //____________________________________________________________________
184 void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured)
186 // fill static variables needed for minuit fit
188 if (!fgCorrelationMatrix)
189 fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
190 if (!fgCorrelationCovarianceMatrix)
191 fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
192 if (!fgCurrentESDVector)
193 fgCurrentESDVector = new TVectorD(fgMaxInput);
194 if (!fgEntropyAPriori)
195 fgEntropyAPriori = new TVectorD(fgMaxParams);
197 fgCorrelationMatrix->Zero();
198 fgCorrelationCovarianceMatrix->Zero();
199 fgCurrentESDVector->Zero();
200 fgEntropyAPriori->Zero();
202 // normalize correction for given nPart
203 for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
205 Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
208 Float_t maxValue = 0;
210 for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
212 // find most probably value
213 if (maxValue < correlation->GetBinContent(i, j))
215 maxValue = correlation->GetBinContent(i, j);
220 correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
221 correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
223 if (i <= fgMaxParams && j <= fgMaxInput)
224 (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
227 //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
231 if (fgNormalizeInput)
232 measured->Scale(1.0 / measured->Integral());
234 for (Int_t i=0; i<fgMaxInput; ++i)
236 (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
237 if (measured->GetBinError(i+1) > 0)
238 (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
240 if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
241 (*fgCorrelationCovarianceMatrix)(i, i) = 0;
242 //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i));
246 //____________________________________________________________________
247 Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
250 // implementation of unfolding (internal function)
252 // unfolds <measured> using response from <correlation> and effiency <efficiency>
253 // output is in <result>
254 // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
255 // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
257 // returns minuit status (0 = success), (-1 when check was set)
260 SetStaticVariables(correlation, measured);
262 // Initialize TMinuit via generic fitter interface
263 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
264 Double_t arglist[100];
266 // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
268 minuit->ExecuteCommand("SET PRINT", arglist, 1);
270 // however, enable warnings
271 //minuit->ExecuteCommand("SET WAR", arglist, 0);
273 // set minimization function
274 minuit->SetFCN(Chi2Function);
276 for (Int_t i=0; i<fgMaxParams; i++)
277 (*fgEntropyAPriori)[i] = 1;
279 if (!initialConditions) {
280 initialConditions = measured;
282 Printf("AliUnfolding::UnfoldWithMinuit: Using different initial conditions...");
283 //new TCanvas; initialConditions->DrawCopy();
284 if (fgNormalizeInput)
285 initialConditions->Scale(1.0 / initialConditions->Integral());
288 // set initial conditions as a-priori distribution for MRX regularization
289 for (Int_t i=0; i<fgMaxParams; i++)
290 if (initialConditions->GetBinContent(i+1) > 0)
291 (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
293 Double_t* results = new Double_t[fgMaxParams+1];
294 for (Int_t i=0; i<fgMaxParams; ++i)
296 results[i] = initialConditions->GetBinContent(i+1);
298 // minuit sees squared values to prevent it from going negative...
299 results[i] = TMath::Sqrt(results[i]);
301 minuit->SetParameter(i, Form("param%d", i), results[i], fgMinuitStepSize, 0, 0);
306 Chi2Function(dummy, 0, chi2, results, 0);
307 printf("AliUnfolding::UnfoldWithMinuit: Chi2 of initial parameters is = %f\n", chi2);
316 // first param is number of iterations, second is precision....
319 //minuit->ExecuteCommand("SCAN", arglist, 0);
320 Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
321 Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status);
322 //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
323 //minuit->ExecuteCommand("MINOS", arglist, 0);
325 for (Int_t i=0; i<fgMaxParams; ++i)
327 results[i] = minuit->GetParameter(i);
328 Double_t value = results[i] * results[i];
329 // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
330 Double_t error = minuit->GetParError(i) * results[i];
334 if (efficiency->GetBinContent(i+1) > 0)
336 value /= efficiency->GetBinContent(i+1);
337 error /= efficiency->GetBinContent(i+1);
346 result->SetBinContent(i+1, value);
347 result->SetBinError(i+1, error);
358 //____________________________________________________________________
359 Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult)
362 // unfolds a spectrum using the Bayesian method
365 if (measured->Integral() <= 0)
367 Printf("AliUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
371 const Int_t kStartBin = 0;
373 Int_t kMaxM = fgMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
374 Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
376 // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
377 const Double_t kConvergenceLimit = kMaxT * 1e-6;
379 // store information in arrays, to increase processing speed (~ factor 5)
380 Double_t* measuredCopy = new Double_t[kMaxM];
381 Double_t* measuredError = new Double_t[kMaxM];
382 Double_t* prior = new Double_t[kMaxT];
383 Double_t* result = new Double_t[kMaxT];
384 Double_t* efficiency = new Double_t[kMaxT];
386 Double_t** response = new Double_t*[kMaxT];
387 Double_t** inverseResponse = new Double_t*[kMaxT];
388 for (Int_t i=0; i<kMaxT; i++)
390 response[i] = new Double_t[kMaxM];
391 inverseResponse[i] = new Double_t[kMaxM];
395 Float_t measuredIntegral = measured->Integral();
396 for (Int_t m=0; m<kMaxM; m++)
398 measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
399 measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
401 for (Int_t t=0; t<kMaxT; t++)
403 response[t][m] = correlation->GetBinContent(t+1, m+1);
404 inverseResponse[t][m] = 0;
408 for (Int_t t=0; t<kMaxT; t++)
412 efficiency[t] = aEfficiency->GetBinContent(t+1);
417 prior[t] = measuredCopy[t];
421 // pick prior distribution
422 if (initialConditions)
424 printf("Using different starting conditions...\n");
426 Float_t inputDistIntegral = initialConditions->Integral();
427 for (Int_t i=0; i<kMaxT; i++)
428 prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
431 //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
434 for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations< 0; i++)
437 Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i);
439 // calculate IR from Beyes theorem
440 // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
442 Double_t chi2Measured = 0;
443 for (Int_t m=0; m<kMaxM; m++)
446 for (Int_t t = kStartBin; t<kMaxT; t++)
447 norm += response[t][m] * prior[t];
449 // calc. chi2: (measured - response * prior) / error
450 if (measuredError[m] > 0)
452 Double_t value = (measuredCopy[m] - norm) / measuredError[m];
453 chi2Measured += value * value;
458 for (Int_t t = kStartBin; t<kMaxT; t++)
459 inverseResponse[t][m] = response[t][m] * prior[t] / norm;
463 for (Int_t t = kStartBin; t<kMaxT; t++)
464 inverseResponse[t][m] = 0;
467 //Printf("chi2Measured of the last prior is %e", chi2Measured);
469 for (Int_t t = kStartBin; t<kMaxT; t++)
472 for (Int_t m=0; m<kMaxM; m++)
473 value += inverseResponse[t][m] * measuredCopy[m];
475 if (efficiency[t] > 0)
476 result[t] = value / efficiency[t];
481 // draw intermediate result
483 for (Int_t t=0; t<kMaxT; t++)
484 aResult->SetBinContent(t+1, result[t]);
485 aResult->SetMarkerStyle(20+i);
486 aResult->SetMarkerColor(2);
487 aResult->DrawCopy("P SAME HIST");
490 Double_t chi2LastIter = 0;
491 // regularization (simple smoothing)
492 for (Int_t t=kStartBin; t<kMaxT; t++)
494 Float_t newValue = 0;
496 // 0 bin excluded from smoothing
497 if (t > kStartBin+2 && t<kMaxT-1)
499 Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
501 // weight the average with the regularization parameter
502 newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average;
505 newValue = result[t];
507 // calculate chi2 (change from last iteration)
510 Double_t diff = (prior[t] - newValue) / prior[t];
511 chi2LastIter += diff * diff;
516 //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
517 //convergence->Fill(i+1, chi2LastIter);
519 if (fgBayesianIterations < 0 && chi2LastIter < kConvergenceLimit)
521 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);
524 } // end of iterations
526 //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
527 //delete convergence;
529 for (Int_t t=0; t<kMaxT; t++)
530 aResult->SetBinContent(t+1, result[t]);
532 delete[] measuredCopy;
533 delete[] measuredError;
538 for (Int_t i=0; i<kMaxT; i++)
540 delete[] response[i];
541 delete[] inverseResponse[i];
544 delete[] inverseResponse;
549 // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
551 /*printf("Calculating covariance matrix. This may take some time...\n");
553 // check if this is the right one...
554 TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
556 Int_t xBins = hInverseResponseBayes->GetNbinsX();
557 Int_t yBins = hInverseResponseBayes->GetNbinsY();
559 // calculate "unfolding matrix" Mij
560 Float_t matrixM[251][251];
561 for (Int_t i=1; i<=xBins; i++)
563 for (Int_t j=1; j<=yBins; j++)
565 if (fCurrentEfficiency->GetBinContent(i) > 0)
566 matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
568 matrixM[i-1][j-1] = 0;
572 Float_t* vectorn = new Float_t[yBins];
573 for (Int_t j=1; j<=yBins; j++)
574 vectorn[j-1] = fCurrentESD->GetBinContent(j);
576 // first part of covariance matrix, depends on input distribution n(E)
577 Float_t cov1[251][251];
579 Float_t nEvents = fCurrentESD->Integral(); // N
584 for (Int_t k=0; k<xBins; k++)
586 printf("In Cov1: %d\n", k);
587 for (Int_t l=0; l<yBins; l++)
591 // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
592 for (Int_t j=0; j<yBins; j++)
593 cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
594 * (1.0 - vectorn[j] / nEvents);
596 // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
597 for (Int_t i=0; i<yBins; i++)
598 for (Int_t j=0; j<yBins; j++)
602 cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
603 * vectorn[j] / nEvents;
608 printf("Cov1 finished\n");
610 TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
613 for (Int_t i=1; i<=xBins; i++)
614 for (Int_t j=1; j<=yBins; j++)
615 cov->SetBinContent(i, j, cov1[i-1][j-1]);
620 // second part of covariance matrix, depends on response matrix
621 Float_t cov2[251][251];
623 // Cov[P(Er|Cu), P(Es|Cu)] term
624 Float_t covTerm[100][100][100];
625 for (Int_t r=0; r<yBins; r++)
626 for (Int_t u=0; u<xBins; u++)
627 for (Int_t s=0; s<yBins; s++)
630 covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
631 * (1.0 - hResponse->GetBinContent(u+1, r+1));
633 covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
634 * hResponse->GetBinContent(u+1, s+1);
637 for (Int_t k=0; k<xBins; k++)
638 for (Int_t l=0; l<yBins; l++)
641 printf("In Cov2: %d %d\n", k, l);
642 for (Int_t i=0; i<yBins; i++)
643 for (Int_t j=0; j<yBins; j++)
645 //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
646 // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
648 for (Int_t r=0; r<yBins; r++)
649 for (Int_t u=0; u<xBins; u++)
650 for (Int_t s=0; s<yBins; s++)
652 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
653 || hResponse->GetBinContent(u+1, i+1) == 0)
656 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
657 * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
661 cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
665 printf("Cov2 finished\n");
667 for (Int_t i=1; i<=xBins; i++)
668 for (Int_t j=1; j<=yBins; j++)
669 cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
675 //____________________________________________________________________
676 Double_t AliUnfolding::RegularizationPol0(TVectorD& params)
678 // homogenity term for minuit fitting
679 // pure function of the parameters
680 // prefers constant function (pol0)
684 for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
686 Double_t right = params[i];
687 Double_t left = params[i-1];
691 Double_t diff = 1 - left / right;
699 //____________________________________________________________________
700 Double_t AliUnfolding::RegularizationPol1(TVectorD& params)
702 // homogenity term for minuit fitting
703 // pure function of the parameters
704 // prefers linear function (pol1)
708 for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
710 if (params[i-1] == 0)
713 Double_t right = params[i];
714 Double_t middle = params[i-1];
715 Double_t left = params[i-2];
717 Double_t der1 = (right - middle);
718 Double_t der2 = (middle - left);
720 Double_t diff = (der1 - der2) / middle;
728 //____________________________________________________________________
729 Double_t AliUnfolding::RegularizationLog(TVectorD& params)
731 // homogenity term for minuit fitting
732 // pure function of the parameters
733 // prefers linear function (pol1)
737 for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
739 if (params[i-1] == 0)
742 Double_t right = log(params[i]);
743 Double_t middle = log(params[i-1]);
744 Double_t left = log(params[i-2]);
746 Double_t der1 = (right - middle);
747 Double_t der2 = (middle - left);
749 Double_t diff = (der1 - der2) / middle;
757 //____________________________________________________________________
758 Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params)
760 // homogenity term for minuit fitting
761 // pure function of the parameters
762 // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
763 // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
767 for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
769 Double_t right = params[i];
770 Double_t middle = params[i-1];
771 Double_t left = params[i-2];
773 Double_t der1 = (right - middle);
774 Double_t der2 = (middle - left);
776 Double_t diff = (der1 - der2);
784 //____________________________________________________________________
785 Double_t AliUnfolding::RegularizationEntropy(TVectorD& params)
787 // homogenity term for minuit fitting
788 // pure function of the parameters
789 // calculates entropy, from
790 // The method of reduced cross-entropy (M. Schmelling 1993)
792 Double_t paramSum = 0;
794 for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
795 paramSum += params[i];
798 for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
800 //Double_t tmp = params[i] / paramSum;
801 Double_t tmp = params[i];
802 if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
804 chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
811 //____________________________________________________________________
812 void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
815 // fit function for minuit
816 // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
820 static TVectorD paramsVector(fgMaxParams);
821 for (Int_t i=0; i<fgMaxParams; ++i)
822 paramsVector[i] = params[i] * params[i];
824 // calculate penalty factor
825 Double_t penaltyVal = 0;
826 switch (fgRegularizationType)
829 case kPol0: penaltyVal = RegularizationPol0(paramsVector); break;
830 case kPol1: penaltyVal = RegularizationPol1(paramsVector); break;
831 case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break;
832 case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break;
833 case kLog: penaltyVal = RegularizationLog(paramsVector); break;
836 //if (penaltyVal > 1e10)
837 // paramsVector2.Print();
839 penaltyVal *= fgRegularizationWeight;
842 TVectorD measGuessVector(fgMaxInput);
843 measGuessVector = (*fgCorrelationMatrix) * paramsVector;
846 measGuessVector -= (*fgCurrentESDVector);
848 //measGuessVector.Print();
850 TVectorD copy(measGuessVector);
853 // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
854 // normal way is like this:
855 // measGuessVector *= (*fgCorrelationCovarianceMatrix);
856 // optimized way like this:
857 for (Int_t i=0; i<fgMaxInput; ++i)
858 measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
860 //measGuessVector.Print();
862 // (Ad - m) W (Ad - m)
863 // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
864 // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
865 Double_t chi2FromFit = measGuessVector * copy * 1e6;
867 chi2 = chi2FromFit + penaltyVal;
869 static Int_t callCount = 0;
870 if ((callCount++ % 10000) == 0)
871 Printf("AliUnfolding::Chi2Function: Iteration %d: %f %f --> %f", callCount, chi2FromFit, penaltyVal, chi2);
874 //____________________________________________________________________
875 void AliUnfolding::DrawGuess(Double_t *params)
878 // draws residuals of solution suggested by params and effect of regularization
882 static TVectorD paramsVector(fgMaxParams);
883 for (Int_t i=0; i<fgMaxParams; ++i)
884 paramsVector[i] = params[i] * params[i];
887 TVectorD measGuessVector(fgMaxInput);
888 measGuessVector = (*fgCorrelationMatrix) * paramsVector;
891 measGuessVector -= (*fgCurrentESDVector);
893 TVectorD copy(measGuessVector);
897 // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
898 // normal way is like this:
899 // measGuessVector *= (*fgCorrelationCovarianceMatrix);
900 // optimized way like this:
901 for (Int_t i=0; i<fgMaxInput; ++i)
902 measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
904 // (Ad - m) W (Ad - m)
905 // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
906 // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
907 //Double_t chi2FromFit = measGuessVector * copy * 1e6;
910 TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
911 for (Int_t i=0; i<fgMaxInput; ++i)
912 residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
913 new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
916 if (fgRegularizationType != kPol1) {
917 Printf("Drawing not supported");
921 TH1* penalty = new TH1F("penalty", "penalty", fgMaxParams, -0.5, fgMaxParams - 0.5);
923 for (Int_t i=2+1; i<fgMaxParams; ++i)
925 if (params[i-1] == 0)
928 Double_t right = params[i];
929 Double_t middle = params[i-1];
930 Double_t left = params[i-2];
932 Double_t der1 = (right - middle);
933 Double_t der2 = (middle - left);
935 Double_t diff = (der1 - der2) / middle;
937 penalty->SetBinContent(i-1, diff * diff);
939 //Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
941 new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
944 //____________________________________________________________________
945 void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
948 // fit function for minuit
949 // uses the TF1 stored in fgFitFunction
952 for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
953 fgFitFunction->SetParameter(i, params[i]);
955 Double_t* params2 = new Double_t[fgMaxParams];
957 for (Int_t i=0; i<fgMaxParams; ++i)
958 params2[i] = fgFitFunction->Eval(i);
960 Chi2Function(unused1, unused2, chi2, params2, unused3);
968 //____________________________________________________________________
969 Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult)
972 // correct spectrum using minuit chi2 method applying a functional fit
977 Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting.");
981 SetChi2Regularization(kNone, 0);
983 SetStaticVariables(correlation, measured);
985 // Initialize TMinuit via generic fitter interface
986 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar());
988 minuit->SetFCN(TF1Function);
989 for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
991 Double_t lower, upper;
992 fgFitFunction->GetParLimits(i, lower, upper);
993 minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper);
996 Double_t arglist[100];
998 minuit->ExecuteCommand("SET PRINT", arglist, 1);
999 minuit->ExecuteCommand("MIGRAD", arglist, 0);
1000 //minuit->ExecuteCommand("MINOS", arglist, 0);
1002 for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1003 fgFitFunction->SetParameter(i, minuit->GetParameter(i));
1005 for (Int_t i=0; i<fgMaxParams; ++i)
1007 Double_t value = fgFitFunction->Eval(i);
1009 Printf("%d : %f", i, value);
1013 if (efficiency->GetBinContent(i+1) > 0)
1015 value /= efficiency->GetBinContent(i+1);
1020 aResult->SetBinContent(i+1, value);
1021 aResult->SetBinError(i+1, 0);
1027 //____________________________________________________________________
1028 void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured)
1030 // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins
1031 // The same limit is applied to the correlation
1034 for (Int_t i=1; i<=measured->GetNbinsX(); ++i)
1036 if (measured->GetBinContent(i) <= fgOverflowBinLimit)
1045 Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit);
1049 Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin);
1051 TCanvas* canvas = 0;
1055 canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
1056 canvas->Divide(2, 2);
1059 measured->SetStats(kFALSE);
1060 measured->DrawCopy();
1064 correlation->SetStats(kFALSE);
1065 correlation->DrawCopy("COLZ");
1068 measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX()));
1069 for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i)
1071 measured->SetBinContent(i, 0);
1072 measured->SetBinError(i, 0);
1074 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1075 measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin)));
1077 Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin));
1079 for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
1081 correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY()));
1082 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1083 correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin)));
1085 for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j)
1087 correlation->SetBinContent(i, j, 0);
1088 correlation->SetBinError(i, j, 0);
1092 Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!");
1097 measured->DrawCopy();
1101 correlation->DrawCopy("COLZ");