From 19442b862ec527d1f200935229c934c210596a0b Mon Sep 17 00:00:00 2001 From: jgrosseo Date: Wed, 25 Feb 2009 09:59:20 +0000 Subject: [PATCH] refactoring multiplicity analysis new class created that does the pure unfolding adapted the analysis accordingly --- PWG0/AliUnfolding.cxx | 1103 +++++++++++++++++ PWG0/AliUnfolding.h | 100 +- PWG0/PWG0baseLinkDef.h | 1 + PWG0/libPWG0base.pkg | 1 + .../AliMultiplicityCorrection.cxx | 1065 +--------------- PWG0/multiplicity/AliMultiplicityCorrection.h | 52 +- PWG0/multiplicity/correct.C | 3 +- 7 files changed, 1226 insertions(+), 1099 deletions(-) create mode 100644 PWG0/AliUnfolding.cxx diff --git a/PWG0/AliUnfolding.cxx b/PWG0/AliUnfolding.cxx new file mode 100644 index 00000000000..6915330d62f --- /dev/null +++ b/PWG0/AliUnfolding.cxx @@ -0,0 +1,1103 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id: AliUnfolding.cxx 31168 2009-02-23 15:18:45Z jgrosseo $ */ + +// This class allows 1-dimensional unfolding. +// Methods that are implemented are chi2 minimization and bayesian unfolding. +// +// Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch + +#include "AliUnfolding.h" +#include +#include +#include +#include +#include +#include + +TMatrixD* AliUnfolding::fgCorrelationMatrix = 0; +TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0; +TVectorD* AliUnfolding::fgCurrentESDVector = 0; +TVectorD* AliUnfolding::fgEntropyAPriori = 0; + +TF1* AliUnfolding::fgFitFunction = 0; + +AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid; +Int_t AliUnfolding::fgMaxInput = -1; // bins in measured histogram +Int_t AliUnfolding::fgMaxParams = -1; // bins in unfolded histogram = number of fit params +Float_t AliUnfolding::fgOverflowBinLimit = -1; + +AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1; +Float_t AliUnfolding::fgRegularizationWeight = 10000; +Int_t AliUnfolding::fgSkipBinsBegin = 0; +Float_t AliUnfolding::fgMinuitStepSize = 0.1; // (usually not needed to be changed) step size in minimization +Bool_t AliUnfolding::fgNormalizeInput = kFALSE; // normalize input spectrum + +Float_t AliUnfolding::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing) +Int_t AliUnfolding::fgBayesianIterations = 10; // number of iterations in Bayesian method + +Bool_t AliUnfolding::fgDebug = kFALSE; + +ClassImp(AliUnfolding) + +//____________________________________________________________________ +void AliUnfolding::SetUnfoldingMethod(MethodType methodType) +{ + // set unfolding method + fgMethodType = methodType; + + const char* name = 0; + switch (methodType) + { + case kInvalid: name = "INVALID"; break; + case kChi2Minimization: name = "Chi2 Minimization"; break; + case kBayesian: name = "Bayesian unfolding"; break; + case kFunction: name = "Functional fit"; break; + } + Printf("AliUnfolding::SetUnfoldingMethod: %s enabled.", name); +} + +//____________________________________________________________________ +void AliUnfolding::SetCreateOverflowBin(Float_t overflowBinLimit) +{ + // enable the creation of a overflow bin that includes all statistics below the given limit + + fgOverflowBinLimit = overflowBinLimit; + + Printf("AliUnfolding::SetCreateOverflowBin: overflow bin limit set to %f", overflowBinLimit); +} + +//____________________________________________________________________ +void AliUnfolding::SetSkipBinsBegin(Int_t nBins) +{ + // set number of skipped bins in regularization + + fgSkipBinsBegin = nBins; + + Printf("AliUnfolding::SetSkipBinsBegin: skipping %d bins at the beginning of the spectrum in the regularization.", fgSkipBinsBegin); +} + +//____________________________________________________________________ +void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded) +{ + // set number of bins in the input (measured) distribution and in the unfolded distribution + fgMaxInput = nMeasured; + fgMaxParams = nUnfolded; + + Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded); +} + +//____________________________________________________________________ +void AliUnfolding::SetChi2Regularization(RegularizationType type, Float_t weight) +{ + // + // sets the parameters for chi2 minimization + // + + fgRegularizationType = type; + fgRegularizationWeight = weight; + + Printf("AliUnfolding::SetChi2Regularization --> Regularization set to %d with weight %f", (Int_t) type, weight); +} + +//____________________________________________________________________ +void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations) +{ + // + // sets the parameters for Bayesian unfolding + // + + fgBayesianSmoothing = smoothing; + fgBayesianIterations = nIterations; + + Printf("AliUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing); +} + +//____________________________________________________________________ +void AliUnfolding::SetFunction(TF1* function) +{ + // set function for unfolding with a fit function + + fgFitFunction = function; + + Printf("AliUnfolding::SetFunction: Set fit function with %d parameters.", function->GetNpar()); +} + +//____________________________________________________________________ +Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check) +{ + // unfolds with unfolding method fgMethodType + // + // parameters: + // correlation: response matrix as measured vs. generated + // 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. + // measured: the measured spectrum + // initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions. + // result: target for the unfolded result + // check: depends on the unfolding method, see comments in specific functions + + if (fgMaxInput == -1) + { + Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution"); + fgMaxInput = measured->GetNbinsX(); + } + if (fgMaxParams == -1) + { + Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution"); + fgMaxParams = measured->GetNbinsX(); + } + + if (fgOverflowBinLimit > 0) + CreateOverflowBin(correlation, measured); + + switch (fgMethodType) + { + case kInvalid: + { + Printf("AliUnfolding::Unfold: ERROR: Unfolding method not set. Use SetUnfoldingMethod. Exiting..."); + return -1; + } + case kChi2Minimization: + return UnfoldWithMinuit(correlation, efficiency, measured, initialConditions, result, check); + case kBayesian: + return UnfoldWithBayesian(correlation, efficiency, measured, initialConditions, result); + case kFunction: + return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result); + } + return -1; +} + +//____________________________________________________________________ +void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured) +{ + // fill static variables needed for minuit fit + + if (!fgCorrelationMatrix) + fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams); + if (!fgCorrelationCovarianceMatrix) + fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput); + if (!fgCurrentESDVector) + fgCurrentESDVector = new TVectorD(fgMaxInput); + if (!fgEntropyAPriori) + fgEntropyAPriori = new TVectorD(fgMaxParams); + + fgCorrelationMatrix->Zero(); + fgCorrelationCovarianceMatrix->Zero(); + fgCurrentESDVector->Zero(); + fgEntropyAPriori->Zero(); + + // normalize correction for given nPart + for (Int_t i=1; i<=correlation->GetNbinsX(); ++i) + { + Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY()); + if (sum <= 0) + continue; + Float_t maxValue = 0; + Int_t maxBin = -1; + for (Int_t j=1; j<=correlation->GetNbinsY(); ++j) + { + // find most probably value + if (maxValue < correlation->GetBinContent(i, j)) + { + maxValue = correlation->GetBinContent(i, j); + maxBin = j; + } + + // npart sum to 1 + correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum); + correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum); + + if (i <= fgMaxParams && j <= fgMaxInput) + (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j); + } + + //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin)); + } + + //normalize measured + if (fgNormalizeInput) + measured->Scale(1.0 / measured->Integral()); + + for (Int_t i=0; iGetBinContent(i+1); + if (measured->GetBinError(i+1) > 0) + (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1); + + if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7) + (*fgCorrelationCovarianceMatrix)(i, i) = 0; + //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i)); + } +} + +//____________________________________________________________________ +Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check) +{ + // + // implementation of unfolding (internal function) + // + // unfolds using response from and effiency + // output is in + // set the initial values for the minimization, if 0 is used + // if is true no unfolding is made, instead only the chi2 without unfolding is printed + // + // returns minuit status (0 = success), (-1 when check was set) + // + + SetStaticVariables(correlation, measured); + + // Initialize TMinuit via generic fitter interface + TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams); + Double_t arglist[100]; + + // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way... + arglist[0] = 0; + minuit->ExecuteCommand("SET PRINT", arglist, 1); + + // however, enable warnings + //minuit->ExecuteCommand("SET WAR", arglist, 0); + + // set minimization function + minuit->SetFCN(Chi2Function); + + for (Int_t i=0; iDrawCopy(); + if (fgNormalizeInput) + initialConditions->Scale(1.0 / initialConditions->Integral()); + } + + // set initial conditions as a-priori distribution for MRX regularization + for (Int_t i=0; iGetBinContent(i+1) > 0) + (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1); + + Double_t* results = new Double_t[fgMaxParams+1]; + for (Int_t i=0; iGetBinContent(i+1); + + // minuit sees squared values to prevent it from going negative... + results[i] = TMath::Sqrt(results[i]); + + minuit->SetParameter(i, Form("param%d", i), results[i], fgMinuitStepSize, 0, 0); + } + + Int_t dummy = 0; + Double_t chi2 = 0; + Chi2Function(dummy, 0, chi2, results, 0); + printf("AliUnfolding::UnfoldWithMinuit: Chi2 of initial parameters is = %f\n", chi2); + + if (check) + { + DrawGuess(results); + delete[] results; + return -1; + } + + // first param is number of iterations, second is precision.... + arglist[0] = 1e6; + //arglist[1] = 1e-5; + //minuit->ExecuteCommand("SCAN", arglist, 0); + Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1); + Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status); + //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!"); + //minuit->ExecuteCommand("MINOS", arglist, 0); + + for (Int_t i=0; iGetParameter(i); + Double_t value = results[i] * results[i]; + // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i)) + Double_t error = minuit->GetParError(i) * results[i]; + + if (efficiency) + { + if (efficiency->GetBinContent(i+1) > 0) + { + value /= efficiency->GetBinContent(i+1); + error /= efficiency->GetBinContent(i+1); + } + else + { + value = 0; + error = 0; + } + } + + result->SetBinContent(i+1, value); + result->SetBinError(i+1, error); + } + + if (fgDebug) + DrawGuess(results); + + delete[] results; + + return status; +} + +//____________________________________________________________________ +Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult) +{ + // + // unfolds a spectrum using the Bayesian method + // + + if (measured->Integral() <= 0) + { + Printf("AliUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty"); + return -1; + } + + const Int_t kStartBin = 0; + + Int_t kMaxM = fgMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis + Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis + + // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4) + const Double_t kConvergenceLimit = kMaxT * 1e-6; + + // store information in arrays, to increase processing speed (~ factor 5) + Double_t* measuredCopy = new Double_t[kMaxM]; + Double_t* measuredError = new Double_t[kMaxM]; + Double_t* prior = new Double_t[kMaxT]; + Double_t* result = new Double_t[kMaxT]; + Double_t* efficiency = new Double_t[kMaxT]; + + Double_t** response = new Double_t*[kMaxT]; + Double_t** inverseResponse = new Double_t*[kMaxT]; + for (Int_t i=0; iIntegral(); + for (Int_t m=0; mGetBinContent(m+1) / measuredIntegral; + measuredError[m] = measured->GetBinError(m+1) / measuredIntegral; + + for (Int_t t=0; tGetBinContent(t+1, m+1); + inverseResponse[t][m] = 0; + } + } + + for (Int_t t=0; tGetBinContent(t+1); + } + else + efficiency[t] = 1; + + prior[t] = measuredCopy[t]; + result[t] = 0; + } + + // pick prior distribution + if (initialConditions) + { + printf("Using different starting conditions...\n"); + // for normalization + Float_t inputDistIntegral = initialConditions->Integral(); + for (Int_t i=0; iGetBinContent(i+1) / inputDistIntegral; + } + + //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5); + + // unfold... + for (Int_t i=0; i 0) + { + Double_t value = (measuredCopy[m] - norm) / measuredError[m]; + chi2Measured += value * value; + } + + if (norm > 0) + { + for (Int_t t = kStartBin; t 0) + result[t] = value / efficiency[t]; + else + result[t] = 0; + } + + // draw intermediate result + /* + for (Int_t t=0; tSetBinContent(t+1, result[t]); + aResult->SetMarkerStyle(20+i); + aResult->SetMarkerColor(2); + aResult->DrawCopy("P SAME HIST"); + */ + + Double_t chi2LastIter = 0; + // regularization (simple smoothing) + for (Int_t t=kStartBin; t kStartBin+2 && t 1e-5) + { + Double_t diff = (prior[t] - newValue) / prior[t]; + chi2LastIter += diff * diff; + } + + prior[t] = newValue; + } + //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter); + //convergence->Fill(i+1, chi2LastIter); + + if (fgBayesianIterations < 0 && chi2LastIter < kConvergenceLimit) + { + 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); + break; + } + } // end of iterations + + //new TCanvas; convergence->DrawCopy(); gPad->SetLogy(); + //delete convergence; + + for (Int_t t=0; tSetBinContent(t+1, result[t]); + + delete[] measuredCopy; + delete[] measuredError; + delete[] prior; + delete[] result; + delete[] efficiency; + + for (Int_t i=0; iProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX()); + + Int_t xBins = hInverseResponseBayes->GetNbinsX(); + Int_t yBins = hInverseResponseBayes->GetNbinsY(); + + // calculate "unfolding matrix" Mij + Float_t matrixM[251][251]; + for (Int_t i=1; i<=xBins; i++) + { + for (Int_t j=1; j<=yBins; j++) + { + if (fCurrentEfficiency->GetBinContent(i) > 0) + matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i); + else + matrixM[i-1][j-1] = 0; + } + } + + Float_t* vectorn = new Float_t[yBins]; + for (Int_t j=1; j<=yBins; j++) + vectorn[j-1] = fCurrentESD->GetBinContent(j); + + // first part of covariance matrix, depends on input distribution n(E) + Float_t cov1[251][251]; + + Float_t nEvents = fCurrentESD->Integral(); // N + + xBins = 20; + yBins = 20; + + for (Int_t k=0; kClone("cov"); + cov->Reset(); + + for (Int_t i=1; i<=xBins; i++) + for (Int_t j=1; j<=yBins; j++) + cov->SetBinContent(i, j, cov1[i-1][j-1]); + + new TCanvas; + cov->Draw("COLZ"); + + // second part of covariance matrix, depends on response matrix + Float_t cov2[251][251]; + + // Cov[P(Er|Cu), P(Es|Cu)] term + Float_t covTerm[100][100][100]; + for (Int_t r=0; rGetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) + * (1.0 - hResponse->GetBinContent(u+1, r+1)); + else + covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) + * hResponse->GetBinContent(u+1, s+1); + } + + for (Int_t k=0; kGetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0 + || hResponse->GetBinContent(u+1, i+1) == 0) + continue; + + tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u) + * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u) + * covTerm[r][u][s]; + } + + cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov; + } + } + + printf("Cov2 finished\n"); + + for (Int_t i=1; i<=xBins; i++) + for (Int_t j=1; j<=yBins; j++) + cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]); + + new TCanvas; + cov->Draw("COLZ");*/ +} + +//____________________________________________________________________ +Double_t AliUnfolding::RegularizationPol0(TVectorD& params) +{ + // homogenity term for minuit fitting + // pure function of the parameters + // prefers constant function (pol0) + + Double_t chi2 = 0; + + for (Int_t i=1+fgSkipBinsBegin; i 0 && (*fgEntropyAPriori)[i] > 0) + { + chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]); + } + } + + return 100.0 + chi2; +} + +//____________________________________________________________________ +void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t) +{ + // + // fit function for minuit + // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix + // + + // d + static TVectorD paramsVector(fgMaxParams); + for (Int_t i=0; i 1e10) + // paramsVector2.Print(); + + penaltyVal *= fgRegularizationWeight; + + // Ad + TVectorD measGuessVector(fgMaxInput); + measGuessVector = (*fgCorrelationMatrix) * paramsVector; + + // Ad - m + measGuessVector -= (*fgCurrentESDVector); + + //measGuessVector.Print(); + + TVectorD copy(measGuessVector); + + // (Ad - m) W + // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used + // normal way is like this: + // measGuessVector *= (*fgCorrelationCovarianceMatrix); + // optimized way like this: + for (Int_t i=0; i %f", callCount, chi2FromFit, penaltyVal, chi2); +} + +//____________________________________________________________________ +void AliUnfolding::DrawGuess(Double_t *params) +{ + // + // draws residuals of solution suggested by params and effect of regularization + // + + // d + static TVectorD paramsVector(fgMaxParams); + for (Int_t i=0; iSetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6); + new TCanvas; residuals->DrawCopy(); gPad->SetLogy(); + + // draw penalty + if (fgRegularizationType != kPol1) { + Printf("Drawing not supported"); + return; + } + + TH1* penalty = new TH1F("penalty", "penalty", fgMaxParams, -0.5, fgMaxParams - 0.5); + + for (Int_t i=2+1; iSetBinContent(i-1, diff * diff); + + //Printf("%d %f %f %f %f", i-1, left, middle, right, diff); + } + new TCanvas; penalty->DrawCopy(); gPad->SetLogy(); +} + +//____________________________________________________________________ +void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3) +{ + // + // fit function for minuit + // uses the TF1 stored in fgFitFunction + // + + for (Int_t i=0; iGetNpar(); i++) + fgFitFunction->SetParameter(i, params[i]); + + Double_t* params2 = new Double_t[fgMaxParams]; + + for (Int_t i=0; iEval(i); + + Chi2Function(unused1, unused2, chi2, params2, unused3); + + delete[] params2; + + if (fgDebug) + Printf("%f", chi2); +} + +//____________________________________________________________________ +Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult) +{ + // + // correct spectrum using minuit chi2 method applying a functional fit + // + + if (!fgFitFunction) + { + Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting."); + return -1; + } + + SetChi2Regularization(kNone, 0); + + SetStaticVariables(correlation, measured); + + // Initialize TMinuit via generic fitter interface + TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar()); + + minuit->SetFCN(TF1Function); + for (Int_t i=0; iGetNpar(); i++) + { + Double_t lower, upper; + fgFitFunction->GetParLimits(i, lower, upper); + minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper); + } + + Double_t arglist[100]; + arglist[0] = 0; + minuit->ExecuteCommand("SET PRINT", arglist, 1); + minuit->ExecuteCommand("MIGRAD", arglist, 0); + //minuit->ExecuteCommand("MINOS", arglist, 0); + + for (Int_t i=0; iGetNpar(); i++) + fgFitFunction->SetParameter(i, minuit->GetParameter(i)); + + for (Int_t i=0; iEval(i); + if (fgDebug) + Printf("%d : %f", i, value); + + if (efficiency) + { + if (efficiency->GetBinContent(i+1) > 0) + { + value /= efficiency->GetBinContent(i+1); + } + else + value = 0; + } + aResult->SetBinContent(i+1, value); + aResult->SetBinError(i+1, 0); + } + + return 0; +} + +//____________________________________________________________________ +void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured) +{ + // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins + // The same limit is applied to the correlation + + Int_t lastBin = 0; + for (Int_t i=1; i<=measured->GetNbinsX(); ++i) + { + if (measured->GetBinContent(i) <= fgOverflowBinLimit) + { + lastBin = i; + break; + } + } + + if (lastBin == 0) + { + Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit); + return; + } + + Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin); + + TCanvas* canvas = 0; + + if (fgDebug) + { + canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800); + canvas->Divide(2, 2); + + canvas->cd(1); + measured->SetStats(kFALSE); + measured->DrawCopy(); + gPad->SetLogy(); + + canvas->cd(2); + correlation->SetStats(kFALSE); + correlation->DrawCopy("COLZ"); + } + + measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX())); + for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i) + { + measured->SetBinContent(i, 0); + measured->SetBinError(i, 0); + } + // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? + measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin))); + + Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin)); + + for (Int_t i=1; i<=correlation->GetNbinsX(); ++i) + { + correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY())); + // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? + correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin))); + + for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j) + { + correlation->SetBinContent(i, j, 0); + correlation->SetBinError(i, j, 0); + } + } + + Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!"); + + if (canvas) + { + canvas->cd(3); + measured->DrawCopy(); + gPad->SetLogy(); + + canvas->cd(4); + correlation->DrawCopy("COLZ"); + } +} diff --git a/PWG0/AliUnfolding.h b/PWG0/AliUnfolding.h index f07e27b3fe3..d11de2237e6 100644 --- a/PWG0/AliUnfolding.h +++ b/PWG0/AliUnfolding.h @@ -5,7 +5,8 @@ // // class that implements several unfolding methods -// E.g. chi2 minimization and bayesian unfolding +// I.e. chi2 minimization and bayesian unfolding +// The whole class is static and not thread-safe (due to the fact that MINUIT unfolding is not thread-safe) // // TMatrixD, TVectorD defined here, because it does not seem possible to predeclare these (or i do not know how) @@ -19,83 +20,76 @@ class TH1; class TH2; -class TH1F; -class TH2F; -class TH3F; class TF1; -class TCollection; class AliUnfolding : public TObject { public: enum RegularizationType { kNone = 0, kPol0, kPol1, kLog, kEntropy, kCurvature }; - enum MethodType { kChi2Minimization = 0, kBayesian = 1 }; + enum MethodType { kInvalid = -1, kChi2Minimization = 0, kBayesian = 1, kFunction = 2}; - AliUnfolding(); - virtual ~AliUnfolding(); + virtual ~AliUnfolding() {}; - void SetInput(TH2* correlationMatrix, TH1* efficiency, TH1* measured) { fCurrentCorrelation = correlationMatrix; fCurrentEfficiency = efficiency; fCurrentESD = measured; } - void SetInitialConditions(TH1* initialConditions) { fInitialConditions = initialConditions; } - const TH1* GetResult() const { return fResult; } + static void SetUnfoldingMethod(MethodType methodType); + static void SetCreateOverflowBin(Float_t overflowBinLimit); + static void SetSkipBinsBegin(Int_t nBins); + static void SetNbins(Int_t nMeasured, Int_t nUnfolded); + static void SetChi2Regularization(RegularizationType type, Float_t weight); + static void SetMinuitStepSize(Float_t stepSize) { fgMinuitStepSize = stepSize; } + static void SetNormalizeInput(Bool_t flag) { fgNormalizeInput = flag; } + static void SetBayesianParameters(Float_t smoothing, Int_t nIterations); + static void SetFunction(TF1* function); + static void SetDebug(Bool_t flag) { fgDebug = flag; } - static void SetParameters(Int_t measuredBins, Int_t unfoldedBins, Bool_t bigbin) { fMaxInput = measuredBins; fMaxParams = unfoldedBins; fgCreateBigBin = bigbin; } - static void SetChi2MinimizationParameters(RegularizationType type, Float_t weight) { fgRegularizationType = type; fgRegularizationWeight = weight; } - static void SetRegularizationRange(Int_t start, Int_t end) { fgRegularizationRangeStart = start; fgRegularizationRangeEnd = end; } - static void SetBayesianParameters(Float_t smoothing, Int_t nIterations) { fgBayesianSmoothing = smoothing; fgBayesianIterations = nIterations; } + static Int_t Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check = kFALSE); - Int_t ApplyMinuitFit(Bool_t check = kFALSE); - Int_t ApplyBayesianMethod(Bool_t determineError = kTRUE); - Int_t ApplyNBDFit(); - Int_t ApplyLaszloMethod(); + protected: + AliUnfolding() {}; - TH1* StatisticalUncertainty(MethodType methodType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo = 0); + static Int_t UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check); + static Int_t UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult); + static Int_t UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult); + + static void DrawGuess(Double_t *params); + static void CreateOverflowBin(TH2* correlation, TH1* measured); + static void SetStaticVariables(TH2* correlation, TH1* measured); - protected: static Double_t RegularizationPol0(TVectorD& params); static Double_t RegularizationPol1(TVectorD& params); static Double_t RegularizationTotalCurvature(TVectorD& params); static Double_t RegularizationEntropy(TVectorD& params); static Double_t RegularizationLog(TVectorD& params); - static void MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t); - static void MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3); - - void SetupCurrentHists(); - - Int_t UnfoldWithBayesian(* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations); - Int_t UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check); - - Float_t BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u); - - TH1* fCurrentESD; //! measured spectrum - TH2* fCurrentCorrelation; //! correlation matrix - TH1* fCurrentEfficiency; //! efficiency - TH1* fInitialConditions; //! initial conditions - TH1* fResult; //! unfolding result + static void Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t); + static void TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3); // static variable to be accessed by MINUIT - static TMatrixD* fgCorrelationMatrix; //! contains fCurrentCorrelation in matrix form - static TMatrixD* fgCorrelationCovarianceMatrix; //! contains the errors of fCurrentESD - static TVectorD* fgCurrentESDVector; //! contains fCurrentESD - static TVectorD* fgEntropyAPriori; //! a-priori distribution for entropy regularization + static TMatrixD* fgCorrelationMatrix; // contains fCurrentCorrelation in matrix form + static TMatrixD* fgCorrelationCovarianceMatrix; // contains the errors of fCurrentESD + static TVectorD* fgCurrentESDVector; // contains fCurrentESD + static TVectorD* fgEntropyAPriori; // a-priori distribution for entropy regularization + + static TF1* fgFitFunction; // fit function - static TF1* fgNBD; //! negative binomial distribution + // --- configuration params follow --- + static MethodType fgMethodType; // unfolding method to be used + static Int_t fgMaxParams; // bins in unfolded histogram = number of fit params + static Int_t fgMaxInput; // bins in measured histogram + static Float_t fgOverflowBinLimit; // to fix fluctuations at high multiplicities, all entries above the limit are summarized in one bin - static Int_t fgMaxParams; //! bins in unfolded histogram = number of fit params - static Int_t fgMaxInput; //! bins in measured histogram + static RegularizationType fgRegularizationType; // regularization that is used during Chi2 method + static Float_t fgRegularizationWeight; // factor for regularization term + static Int_t fgSkipBinsBegin; // (optional) skip the given number of bins in the regularization + static Float_t fgMinuitStepSize; // (usually not needed to be changed) step size in minimization + static Bool_t fgNormalizeInput; // normalize input spectrum - // configuration params follow - static RegularizationType fgRegularizationType; //! regularization that is used during Chi2 method - static Float_t fgRegularizationWeight; //! factor for regularization term - static Int_t fgRegularizationRangeStart; //! first bin where regularization is applied - static Int_t fgRegularizationRangeEnd; //! last bin + 1 where regularization is applied - static Bool_t fgCreateBigBin; //! to fix fluctuations at high multiplicities, all entries above a certain limit are summarized in one bin + static Float_t fgBayesianSmoothing; // smoothing parameter (0 = no smoothing) + static Int_t fgBayesianIterations; // number of iterations in Bayesian method - static Float_t fgBayesianSmoothing; //! smoothing parameter (0 = no smoothing) - static Int_t fgBayesianIterations; //! number of iterations in Bayesian method - // end of configuration + static Bool_t fgDebug; // debug flag + // --- end of configuration --- - private: +private: AliUnfolding(const AliUnfolding&); AliUnfolding& operator=(const AliUnfolding&); diff --git a/PWG0/PWG0baseLinkDef.h b/PWG0/PWG0baseLinkDef.h index ef2767c54eb..9a61f4b9a5b 100644 --- a/PWG0/PWG0baseLinkDef.h +++ b/PWG0/PWG0baseLinkDef.h @@ -18,6 +18,7 @@ #pragma link C++ class AliPWG0Helper+; #pragma link C++ class AliCorrection+; +#pragma link C++ class AliUnfolding+; #pragma link C++ class AliMultiplicityCorrection+; diff --git a/PWG0/libPWG0base.pkg b/PWG0/libPWG0base.pkg index a23aedbc783..7a17ec14e17 100644 --- a/PWG0/libPWG0base.pkg +++ b/PWG0/libPWG0base.pkg @@ -9,6 +9,7 @@ HDRS = dNdEta/dNdEtaAnalysis.h \ AliCorrectionMatrix3D.h \ AliCorrection.h \ AliPWG0Helper.h \ + AliUnfolding.h \ multiplicity/AliMultiplicityCorrection.h SRCS = $(HDRS:.h=.cxx) diff --git a/PWG0/multiplicity/AliMultiplicityCorrection.cxx b/PWG0/multiplicity/AliMultiplicityCorrection.cxx index cc6f34aeaf3..45d42174fc5 100644 --- a/PWG0/multiplicity/AliMultiplicityCorrection.cxx +++ b/PWG0/multiplicity/AliMultiplicityCorrection.cxx @@ -29,7 +29,6 @@ #include #include #include -#include #include #include #include @@ -43,24 +42,6 @@ ClassImp(AliMultiplicityCorrection) -const Int_t AliMultiplicityCorrection::fgkMaxInput = 120; // bins in measured histogram -const Int_t AliMultiplicityCorrection::fgkMaxParams = 120; // bins in unfolded histogram = number of fit params - -TMatrixD* AliMultiplicityCorrection::fgCorrelationMatrix = 0; -TMatrixD* AliMultiplicityCorrection::fgCorrelationCovarianceMatrix = 0; -TVectorD* AliMultiplicityCorrection::fgCurrentESDVector = 0; -TVectorD* AliMultiplicityCorrection::fgEntropyAPriori = 0; - -Bool_t AliMultiplicityCorrection::fgCreateBigBin = kTRUE; -AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fgRegularizationType = AliMultiplicityCorrection::kPol1; -Float_t AliMultiplicityCorrection::fgRegularizationWeight = 10000; -Int_t AliMultiplicityCorrection::fgNParamsMinuit = AliMultiplicityCorrection::fgkMaxParams; - -TF1* AliMultiplicityCorrection::fgNBD = 0; - -Float_t AliMultiplicityCorrection::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing) -Int_t AliMultiplicityCorrection::fgBayesianIterations = 10; // number of iterations in Bayesian method - // These are the areas where the quality of the unfolding results are evaluated // Default defined here, call SetQualityRegions to change them // unit is in multiplicity (not in bin!) @@ -221,6 +202,10 @@ AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const C } TH1::AddDirectory(oldStatus); + + AliUnfolding::SetNbins(120, 120); + AliUnfolding::SetSkipBinsBegin(1); + AliUnfolding::SetNormalizeInput(kTRUE); } //____________________________________________________________________ @@ -542,330 +527,13 @@ void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, I } //____________________________________________________________________ -Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params) -{ - // homogenity term for minuit fitting - // pure function of the parameters - // prefers constant function (pol0) - - Double_t chi2 = 0; - - // ignore the first bin here. on purpose... - for (Int_t i=2; i 0 && (*fgEntropyAPriori)[i] > 0) - { - chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]); - } - } - - return 100.0 + chi2; -} - -//____________________________________________________________________ -void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3) -{ - // - // fit function for minuit - // does: nbd - // func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50); - // func->SetParNames("scaling", "averagen", "k"); - // func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000); - // func->SetParLimits(1, 0.001, 1000); - // func->SetParLimits(2, 0.001, 1000); - // - - fgNBD->SetParameters(params[0], params[1], params[2]); - - Double_t params2[fgkMaxParams]; - - for (Int_t i=0; iEval(i); - - MinuitFitFunction(unused1, unused2, chi2, params2, unused3); - - printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2); -} - -//____________________________________________________________________ -void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t) -{ - // - // fit function for minuit - // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix - // - - // d - static TVectorD paramsVector(fgNParamsMinuit); - for (Int_t i=0; i 1e10) - // paramsVector2.Print(); - - penaltyVal *= fgRegularizationWeight; - - // Ad - TVectorD measGuessVector(fgkMaxInput); - measGuessVector = (*fgCorrelationMatrix) * paramsVector; - - // Ad - m - measGuessVector -= (*fgCurrentESDVector); - - //measGuessVector.Print(); - - TVectorD copy(measGuessVector); - - // (Ad - m) W - // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used - // normal way is like this: - // measGuessVector *= (*fgCorrelationCovarianceMatrix); - // optimized way like this: - for (Int_t i=0; i %f\n", callCount, chi2FromFit, penaltyVal, chi2); -} - -//____________________________________________________________________ -void AliMultiplicityCorrection::DrawGuess(Double_t *params) -{ - // - // draws residuals of solution suggested by params and effect of regularization - // - - // d - static TVectorD paramsVector(fgNParamsMinuit); - for (Int_t i=0; iSetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6); - new TCanvas; residuals->DrawCopy(); gPad->SetLogy(); - - // draw penalty - if (fgRegularizationType != kPol1) { - Printf("Drawing not supported"); - return; - } - - TH1* penalty = new TH1F("penalty", "penalty", fgNParamsMinuit, -0.5, fgNParamsMinuit - 0.5); - - for (Int_t i=2+1; iSetBinContent(i-1, diff * diff); - - //Printf("%d %f %f %f %f", i-1, left, middle, right, diff); - } - new TCanvas; penalty->DrawCopy(); gPad->SetLogy(); -} - -//____________________________________________________________________ -void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin) +void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType) { // // fills fCurrentESD, fCurrentCorrelation // resets fMultiplicityESDCorrected // - Bool_t display = kFALSE; - Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); fMultiplicityESDCorrected[correlationID]->Reset(); @@ -886,88 +554,11 @@ void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullP } } - fCurrentCorrelation = hist->Project3D("zy"); + fCurrentCorrelation = (TH2*) hist->Project3D("zy"); fCurrentCorrelation->Sumw2(); Printf("AliMultiplicityCorrection::SetupCurrentHists: Statistics information: %.f entries in correlation map; %.f entries in measured spectrum", fCurrentCorrelation->Integral(), fCurrentESD->Integral()); - fLastBinLimit = 0; - for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i) - { - if (fCurrentESD->GetBinContent(i) <= 5) - { - fLastBinLimit = i; - break; - } - } - - Printf("AliMultiplicityCorrection::SetupCurrentHists: Bin limit in measured spectrum determined to be %d", fLastBinLimit); - - if (createBigBin) - { - if (fLastBinLimit > 0) - { - TCanvas* canvas = 0; - - if (display) - { - canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800); - canvas->Divide(2, 2); - - canvas->cd(1); - fCurrentESD->GetXaxis()->SetRangeUser(0, 200); - fCurrentESD->SetStats(kFALSE); - fCurrentESD->SetTitle(";measured multiplicity"); - fCurrentESD->DrawCopy(); - gPad->SetLogy(); - - canvas->cd(2); - fCurrentCorrelation->GetXaxis()->SetRangeUser(0, 250); - fCurrentCorrelation->GetYaxis()->SetRangeUser(0, 250); - fCurrentCorrelation->SetTitle(";true multiplicity;measured multiplicity"); - fCurrentCorrelation->SetStats(kFALSE); - fCurrentCorrelation->DrawCopy("COLZ"); - } - - printf("Bin limit in measured spectrum is %d.\n", fLastBinLimit); - fCurrentESD->SetBinContent(fLastBinLimit, fCurrentESD->Integral(fLastBinLimit, fCurrentESD->GetNbinsX())); - for (Int_t i=fLastBinLimit+1; i<=fCurrentESD->GetNbinsX(); ++i) - { - fCurrentESD->SetBinContent(i, 0); - fCurrentESD->SetBinError(i, 0); - } - // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? - fCurrentESD->SetBinError(fLastBinLimit, TMath::Sqrt(fCurrentESD->GetBinContent(fLastBinLimit))); - - printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(fLastBinLimit), fCurrentESD->GetBinError(fLastBinLimit)); - - for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) - { - fCurrentCorrelation->SetBinContent(i, fLastBinLimit, fCurrentCorrelation->Integral(i, i, fLastBinLimit, fCurrentCorrelation->GetNbinsY())); - // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? - fCurrentCorrelation->SetBinError(i, fLastBinLimit, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, fLastBinLimit))); - - for (Int_t j=fLastBinLimit+1; j<=fCurrentCorrelation->GetNbinsY(); ++j) - { - fCurrentCorrelation->SetBinContent(i, j, 0); - fCurrentCorrelation->SetBinError(i, j, 0); - } - } - - printf("Adjusted correlation matrix!\n"); - - if (canvas) - { - canvas->cd(3); - fCurrentESD->DrawCopy(); - gPad->SetLogy(); - - canvas->cd(4); - fCurrentCorrelation->DrawCopy("COLZ"); - } - } - } - #if 0 // does not help // null bins with one entry Int_t nNulledBins = 0; @@ -1032,294 +623,58 @@ TH1* AliMultiplicityCorrection::GetTriggerEfficiency(Int_t inputRange) } //____________________________________________________________________ -void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight, Int_t minuitParams) -{ - // - // sets the parameters for chi2 minimization - // - - fgRegularizationType = type; - fgRegularizationWeight = weight; - - printf("AliMultiplicityCorrection::SetRegularizationParameters --> Regularization set to %d with weight %f\n", (Int_t) type, weight); - - if (minuitParams > 0) - { - fgNParamsMinuit = minuitParams; - printf("AliMultiplicityCorrection::SetRegularizationParameters --> Number of MINUIT iterations set to %d\n", minuitParams); - } -} - -//____________________________________________________________________ -void AliMultiplicityCorrection::SetBayesianParameters(Float_t smoothing, Int_t nIterations) -{ - // - // sets the parameters for Bayesian unfolding - // - - fgBayesianSmoothing = smoothing; - fgBayesianIterations = nIterations; - - printf("AliMultiplicityCorrection::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing); -} - -//____________________________________________________________________ -Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check) +Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions) { // - // implementation of unfolding (internal function) - // - // unfolds using response from and effiency - // output is in - // set the initial values for the minimization, if 0 is used - // if is true no unfolding is made, instead only the chi2 without unfolding is printed + // correct spectrum using minuit chi2 method // - // returns minuit status (0 = success), (-1 when check was set) + // for description of parameters, see AliUnfolding::Unfold // - //normalize measured - measured->Scale(1.0 / measured->Integral()); - - if (!fgCorrelationMatrix) - fgCorrelationMatrix = new TMatrixD(fgkMaxInput, fgNParamsMinuit); - if (!fgCorrelationCovarianceMatrix) - fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxInput, fgkMaxInput); - if (!fgCurrentESDVector) - fgCurrentESDVector = new TVectorD(fgkMaxInput); - if (!fgEntropyAPriori) - fgEntropyAPriori = new TVectorD(fgNParamsMinuit); - - fgCorrelationMatrix->Zero(); - fgCorrelationCovarianceMatrix->Zero(); - fgCurrentESDVector->Zero(); - fgEntropyAPriori->Zero(); - - // normalize correction for given nPart - for (Int_t i=1; i<=correlation->GetNbinsX(); ++i) + Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); + + AliUnfolding::SetCreateOverflowBin(5); + AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization); + SetupCurrentHists(inputRange, fullPhaseSpace, eventType); + + if (!initialConditions) { - Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY()); - if (sum <= 0) - continue; - Float_t maxValue = 0; - Int_t maxBin = -1; - for (Int_t j=1; j<=correlation->GetNbinsY(); ++j) - { - // find most probably value - if (maxValue < correlation->GetBinContent(i, j)) - { - maxValue = correlation->GetBinContent(i, j); - maxBin = j; - } - - // npart sum to 1 - correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum); - correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum); - - if (i <= fgNParamsMinuit && j <= fgkMaxInput) - (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j); - } - - //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin)); - } - - // Initialize TMinuit via generic fitter interface - TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgNParamsMinuit); - Double_t arglist[100]; - - // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way... - arglist[0] = 0; - minuit->ExecuteCommand("SET PRINT", arglist, 1); - - // however, enable warnings - //minuit->ExecuteCommand("SET WAR", arglist, 0); - - // set minimization function - minuit->SetFCN(MinuitFitFunction); - - for (Int_t i=0; iDrawCopy(); + initialConditions = (TH1*) fCurrentESD->Clone("initialConditions"); initialConditions->Scale(1.0 / initialConditions->Integral()); - } - - for (Int_t i=0; iGetBinContent(i+1) > 0) - (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1); - - Double_t* results = new Double_t[fgNParamsMinuit+1]; - for (Int_t i=0; iGetBinContent(i+1); - - // minimum value if (!check) - results[i] = TMath::Max(results[i], 1e-3); - - Float_t stepSize = 0.1; - - // minuit sees squared values to prevent it from going negative... - results[i] = TMath::Sqrt(results[i]); - //stepSize /= results[i]; // keep relative step size - - minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0); - } - //minuit->SetParameter(fgkMaxParams, "alpha", 1e4, 1, 0, 0); - // bin 0 is filled with value from bin 1 (otherwise it's 0) - //minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 0); - //results[0] = 0; - //minuit->SetParameter(0, "param0", 0, 0, 0, 0); - - for (Int_t i=0; iGetBinContent(i+1); - if (measured->GetBinError(i+1) > 0) - (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1); - - if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7) - (*fgCorrelationCovarianceMatrix)(i, i) = 0; - //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i)); - } - - Int_t dummy = 0; - Double_t chi2 = 0; - MinuitFitFunction(dummy, 0, chi2, results, 0); - printf("Chi2 of initial parameters is = %f\n", chi2); - - if (check) - { - DrawGuess(results); - delete[] results; - return -1; - } - - // first param is number of iterations, second is precision.... - arglist[0] = 1e6; - //arglist[1] = 1e-5; - //minuit->ExecuteCommand("SCAN", arglist, 0); - Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1); - printf("MINUIT status is %d\n", status); - //minuit->ExecuteCommand("MIGRAD", arglist, 1); - //minuit->ExecuteCommand("MIGRAD", arglist, 1); - //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!"); - //minuit->ExecuteCommand("MINOS", arglist, 0); - - //new TCanvas; aEfficiency->Draw(); - - for (Int_t i=0; iGetParameter(i); - if (aEfficiency->GetBinContent(i+1) > 0) - { - result->SetBinContent(i+1, results[i] * results[i] / aEfficiency->GetBinContent(i+1)); - // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i)) - result->SetBinError(i+1, minuit->GetParError(i) * results[i] / aEfficiency->GetBinContent(i+1)); - } - else { - result->SetBinContent(i+1, 0); - result->SetBinError(i+1, 0); + // set minimum value to prevent MINUIT just staying in the small value + for (Int_t i=1; i<=initialConditions->GetNbinsX(); i++) + initialConditions->SetBinContent(i, TMath::Max(initialConditions->GetBinContent(i), 1e-3)); } } - //DrawGuess(results); - delete[] results; - - return status; + return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check); } //____________________________________________________________________ -Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions) +Int_t AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType) { // - // correct spectrum using minuit chi2 method - // - // for description of parameters, see UnfoldWithMinuit - // - - Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); - SetupCurrentHists(inputRange, fullPhaseSpace, eventType, fgCreateBigBin); - - return UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check); -} - -//____________________________________________________________________ -void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace) -{ + // correct spectrum using minuit chi2 method with a NBD function // - // correct spectrum using minuit chi2 method applying a NBD fit + // for description of parameters, see AliUnfolding::Unfold // Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); - - SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE); - //normalize ESD - fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); - - fgCorrelationMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams); - fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams); - fgCurrentESDVector = new TVectorD(fgkMaxParams); - - // normalize correction for given nPart - for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) - { - Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); - if (sum <= 0) - continue; - for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) - { - // npart sum to 1 - fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); - fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); - - if (i <= fgkMaxParams && j <= fgkMaxParams) - (*fgCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j); - } - } - - for (Int_t i=0; iGetBinContent(i+1); - if (fCurrentESD->GetBinError(i+1) > 0) - (*fgCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1)); - } - - // Create NBD function - if (!fgNBD) - { - fgNBD = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 250); - fgNBD->SetParNames("scaling", "averagen", "k"); - } - - // Initialize TMinuit via generic fitter interface - TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3); - - minuit->SetFCN(MinuitNBD); - minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000); - minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000); - minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000); - - Double_t arglist[100]; - arglist[0] = 0; - minuit->ExecuteCommand("SET PRINT", arglist, 1); - minuit->ExecuteCommand("MIGRAD", arglist, 0); - //minuit->ExecuteCommand("MINOS", arglist, 0); - - fgNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2)); - - for (Int_t i=0; iEval(i)); - if (fgNBD->Eval(i) > 0) - { - fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fgNBD->Eval(i)); - fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0); - } - } + + AliUnfolding::SetUnfoldingMethod(AliUnfolding::kFunction); + SetupCurrentHists(inputRange, fullPhaseSpace, eventType); + + TF1* func = new TF1("nbd", "[0] * TMath::Gamma([2]+x) / TMath::Gamma([2]) / TMath::Gamma(x+1) * pow([1] / ([1]+[2]), x) * pow(1.0 + [1]/[2], -[2])"); + func->SetParNames("scaling", "averagen", "k"); + func->SetParLimits(0, 0, 1000); + func->SetParLimits(1, 1, 50); + func->SetParLimits(2, 1, 10); + func->SetParameters(1, 10, 2); + AliUnfolding::SetFunction(func); + + return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, 0, fMultiplicityESDCorrected[correlationID]); } //____________________________________________________________________ @@ -1449,6 +804,17 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5); + // find bin limit + Int_t lastBin = 0; + for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i) + { + if (fCurrentESD->GetBinContent(i) <= 5) + { + lastBin = i; + break; + } + } + // TODO fix errors Float_t chi2 = 0; for (Int_t i=1; i<=residual->GetNbinsX(); ++i) @@ -1456,7 +822,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang if (fCurrentESD->GetBinError(i) > 0) { Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i); - if (i > 1 && i <= fLastBinLimit) + if (i > 1 && i <= lastBin) chi2 += value * value; Printf("%d --> %f (%f)", i, value * value, chi2); residual->SetBinContent(i, value); @@ -1477,7 +843,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang //residualHist->Fit("gaus", "N"); //delete residualHist; - printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, fLastBinLimit); + printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, lastBin); TCanvas* canvas1 = 0; if (simple) @@ -1796,7 +1162,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang Double_t newChi2 = 0; Double_t newChi2Limit150 = 0; Int_t ndf = 0; - for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgkMaxParams+1); ++i) + for (Int_t i=1; i<=diffMCUnfolded2->GetNbinsX(); ++i) { Double_t value = 0; if (proj->GetBinError(i) > 0) @@ -2027,7 +1393,7 @@ TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max) } //____________________________________________________________________ -TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo) +TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo) { // // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix @@ -2043,6 +1409,10 @@ TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, In // initialize seed with current time gRandom->SetSeed(0); + + if (methodType == AliUnfolding::kChi2Minimization) + AliUnfolding::SetCreateOverflowBin(5); + AliUnfolding::SetUnfoldingMethod(methodType); const Int_t kErrorIterations = 150; @@ -2055,12 +1425,7 @@ TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, In { Printf("Iteration %d of %d...", n, kErrorIterations); - // only done for chi2 minimization - Bool_t createBigBin = kFALSE; - if (methodType == kChi2Minimization) - createBigBin = kTRUE; - - SetupCurrentHists(inputRange, fullPhaseSpace, eventType, createBigBin); + SetupCurrentHists(inputRange, fullPhaseSpace, eventType); TH1* measured = (TH1*) fCurrentESD->Clone("measured"); @@ -2087,7 +1452,7 @@ TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, In } // only for bayesian method we have to do it before the call to Unfold... - if (methodType == kBayesian) + if (methodType == AliUnfolding::kBayesian) { for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) { @@ -2127,13 +1492,7 @@ TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, In { result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n)); - Int_t returnCode = -1; - if (methodType == kChi2Minimization) - { - returnCode = UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, kFALSE); - } - else if (methodType == kBayesian) - returnCode = UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, fgBayesianSmoothing, fgBayesianIterations); + Int_t returnCode = AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result); if (returnCode != 0) return 0; @@ -2286,7 +1645,7 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful // initialize seed with current time gRandom->SetSeed(0); - SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE); + SetupCurrentHists(inputRange, fullPhaseSpace, eventType); // normalize correction for given nPart for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) @@ -2317,7 +1676,9 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); - if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], regPar, nIterations) != 0) + AliUnfolding::SetBayesianParameters(regPar, nIterations); + AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian); + if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]) != 0) return; if (!determineError) @@ -2350,7 +1711,7 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2"); result2->Reset(); - if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2, regPar, nIterations) != 0) + if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2) != 0) return; result2->Scale(1.0 / result2->Integral()); @@ -2371,300 +1732,6 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful delete error; } -//____________________________________________________________________ -Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations) -{ - // - // unfolds a spectrum - // - - if (measured->Integral() <= 0) - { - Printf("AliMultiplicityCorrection::UnfoldWithBayesian: ERROR: The measured spectrum is empty"); - return -1; - } - - const Int_t kStartBin = 0; - - const Int_t kMaxM = fgkMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis - const Int_t kMaxT = fgkMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis - - // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4) - const Double_t kConvergenceLimit = kMaxT * 1e-6; - - // store information in arrays, to increase processing speed (~ factor 5) - Double_t measuredCopy[kMaxM]; - Double_t measuredError[kMaxM]; - Double_t prior[kMaxT]; - Double_t result[kMaxT]; - Double_t efficiency[kMaxT]; - - Double_t response[kMaxT][kMaxM]; - Double_t inverseResponse[kMaxT][kMaxM]; - - // for normalization - Float_t measuredIntegral = measured->Integral(); - for (Int_t m=0; mGetBinContent(m+1) / measuredIntegral; - measuredError[m] = measured->GetBinError(m+1) / measuredIntegral; - - for (Int_t t=0; tGetBinContent(t+1, m+1); - inverseResponse[t][m] = 0; - } - } - - for (Int_t t=0; tGetBinContent(t+1); - prior[t] = measuredCopy[t]; - result[t] = 0; - } - - // pick prior distribution - if (initialConditions) - { - printf("Using different starting conditions...\n"); - // for normalization - Float_t inputDistIntegral = initialConditions->Integral(); - for (Int_t i=0; iGetBinContent(i+1) / inputDistIntegral; - } - - //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5); - - // unfold... - for (Int_t i=0; i 0) - { - Double_t value = (measuredCopy[m] - norm) / measuredError[m]; - chi2Measured += value * value; - } - - if (norm > 0) - { - for (Int_t t = kStartBin; t 0) - result[t] = value / efficiency[t]; - else - result[t] = 0; - } - - // draw intermediate result - /* - for (Int_t t=0; tSetBinContent(t+1, result[t]); - aResult->SetMarkerStyle(20+i); - aResult->SetMarkerColor(2); - aResult->DrawCopy("P SAME HIST"); - */ - - Double_t chi2LastIter = 0; - // regularization (simple smoothing) - for (Int_t t=kStartBin; t kStartBin+2 && t 1e-5) - { - Double_t diff = (prior[t] - newValue) / prior[t]; - chi2LastIter += diff * diff; - } - - prior[t] = newValue; - } - //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter); - //convergence->Fill(i+1, chi2LastIter); - - if (nIterations < 0 && chi2LastIter < kConvergenceLimit) - { - Printf("Stopped Bayesian unfolding after %d iterations at chi2(change since last iteration) of %e; chi2Measured of the last prior is %e", i, chi2LastIter, chi2Measured); - break; - } - - //if (i % 5 == 0) - // DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType); - } // end of iterations - - //new TCanvas; convergence->DrawCopy(); gPad->SetLogy(); - //delete convergence; - - for (Int_t t=0; tSetBinContent(t+1, result[t]); - - return 0; - - // ******** - // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995 - - /*printf("Calculating covariance matrix. This may take some time...\n"); - - // TODO check if this is the right one... - TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX()); - - Int_t xBins = hInverseResponseBayes->GetNbinsX(); - Int_t yBins = hInverseResponseBayes->GetNbinsY(); - - // calculate "unfolding matrix" Mij - Float_t matrixM[251][251]; - for (Int_t i=1; i<=xBins; i++) - { - for (Int_t j=1; j<=yBins; j++) - { - if (fCurrentEfficiency->GetBinContent(i) > 0) - matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i); - else - matrixM[i-1][j-1] = 0; - } - } - - Float_t* vectorn = new Float_t[yBins]; - for (Int_t j=1; j<=yBins; j++) - vectorn[j-1] = fCurrentESD->GetBinContent(j); - - // first part of covariance matrix, depends on input distribution n(E) - Float_t cov1[251][251]; - - Float_t nEvents = fCurrentESD->Integral(); // N - - xBins = 20; - yBins = 20; - - for (Int_t k=0; kClone("cov"); - cov->Reset(); - - for (Int_t i=1; i<=xBins; i++) - for (Int_t j=1; j<=yBins; j++) - cov->SetBinContent(i, j, cov1[i-1][j-1]); - - new TCanvas; - cov->Draw("COLZ"); - - // second part of covariance matrix, depends on response matrix - Float_t cov2[251][251]; - - // Cov[P(Er|Cu), P(Es|Cu)] term - Float_t covTerm[100][100][100]; - for (Int_t r=0; rGetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) - * (1.0 - hResponse->GetBinContent(u+1, r+1)); - else - covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1) - * hResponse->GetBinContent(u+1, s+1); - } - - for (Int_t k=0; kGetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0 - || hResponse->GetBinContent(u+1, i+1) == 0) - continue; - - tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u) - * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u) - * covTerm[r][u][s]; - } - - cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov; - } - } - - printf("Cov2 finished\n"); - - for (Int_t i=1; i<=xBins; i++) - for (Int_t j=1; j<=yBins; j++) - cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]); - - new TCanvas; - cov->Draw("COLZ");*/ -} - //____________________________________________________________________ Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u) { @@ -2700,7 +1767,7 @@ void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullP Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); - SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE); + SetupCurrentHists(inputRange, fullPhaseSpace, eventType); //normalize ESD fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); @@ -2822,7 +1889,7 @@ void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t ful Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); - SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE); + SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx); //normalize ESD fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); diff --git a/PWG0/multiplicity/AliMultiplicityCorrection.h b/PWG0/multiplicity/AliMultiplicityCorrection.h index fede05b5cc9..63204d3013b 100644 --- a/PWG0/multiplicity/AliMultiplicityCorrection.h +++ b/PWG0/multiplicity/AliMultiplicityCorrection.h @@ -27,12 +27,11 @@ class TCollection; #include #include #include +#include class AliMultiplicityCorrection : public TNamed { public: enum EventType { kTrVtx = 0, kMB, kINEL, kNSD }; - enum RegularizationType { kNone = 0, kPol0, kPol1, kLog, kEntropy, kCurvature }; - enum MethodType { kChi2Minimization = 0, kBayesian = 1 }; enum { kESDHists = 4, kMCHists = 5, kCorrHists = 8, kQualityRegions = 3 }; AliMultiplicityCorrection(); @@ -55,17 +54,12 @@ class AliMultiplicityCorrection : public TNamed { Int_t ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* initialConditions = 0); - static void SetRegularizationParameters(RegularizationType type, Float_t weight, Int_t minuitParams = -1); - static void SetBayesianParameters(Float_t smoothing, Int_t nIterations); - static void SetCreateBigBin(Bool_t flag) { fgCreateBigBin = flag; } - - void ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace); - void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 1, Int_t nIterations = 100, TH1* initialConditions = 0, Bool_t determineError = kTRUE); static TH1* CalculateStdDev(TH1** results, Int_t max); - TH1* StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo = 0); + TH1* StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo = 0); + Int_t ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType); void ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace); void ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType); @@ -102,48 +96,14 @@ class AliMultiplicityCorrection : public TNamed { void FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y); protected: - static const Int_t fgkMaxParams; //! bins in unfolded histogram = number of fit params - static const Int_t fgkMaxInput; //! bins in measured histogram - - static Double_t RegularizationPol0(TVectorD& params); - static Double_t RegularizationPol1(TVectorD& params); - static Double_t RegularizationTotalCurvature(TVectorD& params); - static Double_t RegularizationEntropy(TVectorD& params); - static Double_t RegularizationLog(TVectorD& params); - - static void MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t); - static void MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3); - - static void DrawGuess(Double_t *params); - - void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin); + void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType); Float_t BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u); - static Int_t UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations); - static Int_t UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check); - + TH1* fCurrentESD; //! current input esd - TH1* fCurrentCorrelation; //! current correlation + TH2* fCurrentCorrelation; //! current correlation TH1* fCurrentEfficiency; //! current efficiency - // static variable to be accessed by MINUIT - static TMatrixD* fgCorrelationMatrix; //! contains fCurrentCorrelation in matrix form - static TMatrixD* fgCorrelationCovarianceMatrix; //! contains the errors of fCurrentESD - static TVectorD* fgCurrentESDVector; //! contains fCurrentESD - static TVectorD* fgEntropyAPriori; //! a-priori distribution for entropy regularization - - static TF1* fgNBD; //! negative binomial distribution - - // configuration params follow - static RegularizationType fgRegularizationType; //! regularization that is used during Chi2 method - static Float_t fgRegularizationWeight; //! factor for regularization term - static Bool_t fgCreateBigBin; //! to fix fluctuations at high multiplicities, all entries above a certain limit are summarized in one bin - static Int_t fgNParamsMinuit; //! number of parameters minuit uses for unfolding (todo: to be merged w/ fgkMaxParams that has to be const. for the moment) - - static Float_t fgBayesianSmoothing; //! smoothing parameter (0 = no smoothing) - static Int_t fgBayesianIterations; //! number of iterations in Bayesian method - // end of configuration - TH2F* fMultiplicityESD[kESDHists]; // multiplicity histogram: vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.5, 2 (0..3) TH2F* fMultiplicityVtx[kMCHists]; // multiplicity histogram of events that have a reconstructed vertex : vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.5, 2, inf (0..4) diff --git a/PWG0/multiplicity/correct.C b/PWG0/multiplicity/correct.C index c69255372b8..25fd74fee7a 100644 --- a/PWG0/multiplicity/correct.C +++ b/PWG0/multiplicity/correct.C @@ -88,7 +88,7 @@ void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder if (chi2) { - mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, beta); + AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, beta); //mult->SetCreateBigBin(kFALSE); //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0); //mult->SetCreateBigBin(kFALSE); //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125); mult->SetCreateBigBin(kFALSE); @@ -99,6 +99,7 @@ void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder //mult->SetMultiplicityESDCorrected(histID, (TH1F*) mcCompare); mult->ApplyMinuitFit(histID, fullPhaseSpace, eventType, kFALSE); //hist2->ProjectionY("mymchist")); + //mult->ApplyNBDFit(histID, fullPhaseSpace, eventType); } else { -- 2.39.3