--- /dev/null
+/**************************************************************************
+ * 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 <TH1F.h>
+#include <TH2F.h>
+#include <TVirtualFitter.h>
+#include <TMath.h>
+#include <TCanvas.h>
+#include <TF1.h>
+
+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; i<fgMaxInput; ++i)
+ {
+ (*fgCurrentESDVector)[i] = measured->GetBinContent(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 <measured> using response from <correlation> and effiency <efficiency>
+ // output is in <result>
+ // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
+ // if <check> 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; i<fgMaxParams; i++)
+ (*fgEntropyAPriori)[i] = 1;
+
+ if (!initialConditions) {
+ initialConditions = measured;
+ } else {
+ Printf("AliUnfolding::UnfoldWithMinuit: Using different initial conditions...");
+ //new TCanvas; initialConditions->DrawCopy();
+ if (fgNormalizeInput)
+ initialConditions->Scale(1.0 / initialConditions->Integral());
+ }
+
+ // set initial conditions as a-priori distribution for MRX regularization
+ for (Int_t i=0; i<fgMaxParams; i++)
+ if (initialConditions->GetBinContent(i+1) > 0)
+ (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
+
+ Double_t* results = new Double_t[fgMaxParams+1];
+ for (Int_t i=0; i<fgMaxParams; ++i)
+ {
+ results[i] = initialConditions->GetBinContent(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; i<fgMaxParams; ++i)
+ {
+ results[i] = minuit->GetParameter(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; i<kMaxT; i++)
+ {
+ response[i] = new Double_t[kMaxM];
+ inverseResponse[i] = new Double_t[kMaxM];
+ }
+
+ // for normalization
+ Float_t measuredIntegral = measured->Integral();
+ for (Int_t m=0; m<kMaxM; m++)
+ {
+ measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
+ measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
+
+ for (Int_t t=0; t<kMaxT; t++)
+ {
+ response[t][m] = correlation->GetBinContent(t+1, m+1);
+ inverseResponse[t][m] = 0;
+ }
+ }
+
+ for (Int_t t=0; t<kMaxT; t++)
+ {
+ if (efficiency)
+ {
+ efficiency[t] = aEfficiency->GetBinContent(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; i<kMaxT; i++)
+ prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
+ }
+
+ //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
+
+ // unfold...
+ for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations< 0; i++)
+ {
+ if (fgDebug)
+ Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i);
+
+ // calculate IR from Beyes theorem
+ // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
+
+ Double_t chi2Measured = 0;
+ for (Int_t m=0; m<kMaxM; m++)
+ {
+ Float_t norm = 0;
+ for (Int_t t = kStartBin; t<kMaxT; t++)
+ norm += response[t][m] * prior[t];
+
+ // calc. chi2: (measured - response * prior) / error
+ if (measuredError[m] > 0)
+ {
+ Double_t value = (measuredCopy[m] - norm) / measuredError[m];
+ chi2Measured += value * value;
+ }
+
+ if (norm > 0)
+ {
+ for (Int_t t = kStartBin; t<kMaxT; t++)
+ inverseResponse[t][m] = response[t][m] * prior[t] / norm;
+ }
+ else
+ {
+ for (Int_t t = kStartBin; t<kMaxT; t++)
+ inverseResponse[t][m] = 0;
+ }
+ }
+ //Printf("chi2Measured of the last prior is %e", chi2Measured);
+
+ for (Int_t t = kStartBin; t<kMaxT; t++)
+ {
+ Float_t value = 0;
+ for (Int_t m=0; m<kMaxM; m++)
+ value += inverseResponse[t][m] * measuredCopy[m];
+
+ if (efficiency[t] > 0)
+ result[t] = value / efficiency[t];
+ else
+ result[t] = 0;
+ }
+
+ // draw intermediate result
+ /*
+ for (Int_t t=0; t<kMaxT; t++)
+ aResult->SetBinContent(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<kMaxT; t++)
+ {
+ Float_t newValue = 0;
+
+ // 0 bin excluded from smoothing
+ if (t > kStartBin+2 && t<kMaxT-1)
+ {
+ Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
+
+ // weight the average with the regularization parameter
+ newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average;
+ }
+ else
+ newValue = result[t];
+
+ // calculate chi2 (change from last iteration)
+ if (prior[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; t<kMaxT; t++)
+ aResult->SetBinContent(t+1, result[t]);
+
+ delete[] measuredCopy;
+ delete[] measuredError;
+ delete[] prior;
+ delete[] result;
+ delete[] efficiency;
+
+ for (Int_t i=0; i<kMaxT; i++)
+ {
+ delete[] response[i];
+ delete[] inverseResponse[i];
+ }
+ delete[] response;
+ delete[] inverseResponse;
+
+ 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");
+
+ // 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; k<xBins; k++)
+ {
+ printf("In Cov1: %d\n", k);
+ for (Int_t l=0; l<yBins; l++)
+ {
+ cov1[k][l] = 0;
+
+ // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
+ for (Int_t j=0; j<yBins; j++)
+ cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
+ * (1.0 - vectorn[j] / nEvents);
+
+ // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
+ for (Int_t i=0; i<yBins; i++)
+ for (Int_t j=0; j<yBins; j++)
+ {
+ if (i == j)
+ continue;
+ cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
+ * vectorn[j] / nEvents;
+ }
+ }
+ }
+
+ printf("Cov1 finished\n");
+
+ TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("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; r<yBins; r++)
+ for (Int_t u=0; u<xBins; u++)
+ for (Int_t s=0; s<yBins; s++)
+ {
+ if (r == s)
+ covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(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; k<xBins; k++)
+ for (Int_t l=0; l<yBins; l++)
+ {
+ cov2[k][l] = 0;
+ printf("In Cov2: %d %d\n", k, l);
+ for (Int_t i=0; i<yBins; i++)
+ for (Int_t j=0; j<yBins; j++)
+ {
+ //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
+ // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
+ Float_t tmpCov = 0;
+ for (Int_t r=0; r<yBins; r++)
+ for (Int_t u=0; u<xBins; u++)
+ for (Int_t s=0; s<yBins; s++)
+ {
+ if (hResponse->GetBinContent(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<fgMaxParams; ++i)
+ {
+ Double_t right = params[i];
+ Double_t left = params[i-1];
+
+ if (right != 0)
+ {
+ Double_t diff = 1 - left / right;
+ chi2 += diff * diff;
+ }
+ }
+
+ return chi2 / 100.0;
+}
+
+//____________________________________________________________________
+Double_t AliUnfolding::RegularizationPol1(TVectorD& params)
+{
+ // homogenity term for minuit fitting
+ // pure function of the parameters
+ // prefers linear function (pol1)
+
+ Double_t chi2 = 0;
+
+ for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
+ {
+ if (params[i-1] == 0)
+ continue;
+
+ Double_t right = params[i];
+ Double_t middle = params[i-1];
+ Double_t left = params[i-2];
+
+ Double_t der1 = (right - middle);
+ Double_t der2 = (middle - left);
+
+ Double_t diff = (der1 - der2) / middle;
+
+ chi2 += diff * diff;
+ }
+
+ return chi2;
+}
+
+//____________________________________________________________________
+Double_t AliUnfolding::RegularizationLog(TVectorD& params)
+{
+ // homogenity term for minuit fitting
+ // pure function of the parameters
+ // prefers linear function (pol1)
+
+ Double_t chi2 = 0;
+
+ for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
+ {
+ if (params[i-1] == 0)
+ continue;
+
+ Double_t right = log(params[i]);
+ Double_t middle = log(params[i-1]);
+ Double_t left = log(params[i-2]);
+
+ Double_t der1 = (right - middle);
+ Double_t der2 = (middle - left);
+
+ Double_t diff = (der1 - der2) / middle;
+
+ chi2 += diff * diff;
+ }
+
+ return chi2 * 100;
+}
+
+//____________________________________________________________________
+Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params)
+{
+ // homogenity term for minuit fitting
+ // pure function of the parameters
+ // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
+ // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
+
+ Double_t chi2 = 0;
+
+ for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
+ {
+ Double_t right = params[i];
+ Double_t middle = params[i-1];
+ Double_t left = params[i-2];
+
+ Double_t der1 = (right - middle);
+ Double_t der2 = (middle - left);
+
+ Double_t diff = (der1 - der2);
+
+ chi2 += diff * diff;
+ }
+
+ return chi2 * 1e4;
+}
+
+//____________________________________________________________________
+Double_t AliUnfolding::RegularizationEntropy(TVectorD& params)
+{
+ // homogenity term for minuit fitting
+ // pure function of the parameters
+ // calculates entropy, from
+ // The method of reduced cross-entropy (M. Schmelling 1993)
+
+ Double_t paramSum = 0;
+
+ for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
+ paramSum += params[i];
+
+ Double_t chi2 = 0;
+ for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
+ {
+ //Double_t tmp = params[i] / paramSum;
+ Double_t tmp = params[i];
+ if (tmp > 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<fgMaxParams; ++i)
+ paramsVector[i] = params[i] * params[i];
+
+ // calculate penalty factor
+ Double_t penaltyVal = 0;
+ switch (fgRegularizationType)
+ {
+ case kNone: break;
+ case kPol0: penaltyVal = RegularizationPol0(paramsVector); break;
+ case kPol1: penaltyVal = RegularizationPol1(paramsVector); break;
+ case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break;
+ case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break;
+ case kLog: penaltyVal = RegularizationLog(paramsVector); break;
+ }
+
+ //if (penaltyVal > 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<fgMaxInput; ++i)
+ measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
+
+ //measGuessVector.Print();
+
+ // (Ad - m) W (Ad - m)
+ // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
+ // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
+ Double_t chi2FromFit = measGuessVector * copy * 1e6;
+
+ chi2 = chi2FromFit + penaltyVal;
+
+ static Int_t callCount = 0;
+ if ((callCount++ % 10000) == 0)
+ Printf("AliUnfolding::Chi2Function: Iteration %d: %f %f --> %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; i<fgMaxParams; ++i)
+ paramsVector[i] = params[i] * params[i];
+
+ // Ad
+ TVectorD measGuessVector(fgMaxInput);
+ measGuessVector = (*fgCorrelationMatrix) * paramsVector;
+
+ // Ad - m
+ measGuessVector -= (*fgCurrentESDVector);
+
+ TVectorD copy(measGuessVector);
+ //copy.Print();
+
+ // (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<fgMaxInput; ++i)
+ measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
+
+ // (Ad - m) W (Ad - m)
+ // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
+ // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
+ //Double_t chi2FromFit = measGuessVector * copy * 1e6;
+
+ // draw residuals
+ TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
+ for (Int_t i=0; i<fgMaxInput; ++i)
+ residuals->SetBinContent(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; i<fgMaxParams; ++i)
+ {
+ if (params[i-1] == 0)
+ continue;
+
+ Double_t right = params[i];
+ Double_t middle = params[i-1];
+ Double_t left = params[i-2];
+
+ Double_t der1 = (right - middle);
+ Double_t der2 = (middle - left);
+
+ Double_t diff = (der1 - der2) / middle;
+
+ penalty->SetBinContent(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; i<fgFitFunction->GetNpar(); i++)
+ fgFitFunction->SetParameter(i, params[i]);
+
+ Double_t* params2 = new Double_t[fgMaxParams];
+
+ for (Int_t i=0; i<fgMaxParams; ++i)
+ params2[i] = fgFitFunction->Eval(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; i<fgFitFunction->GetNpar(); 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; i<fgFitFunction->GetNpar(); i++)
+ fgFitFunction->SetParameter(i, minuit->GetParameter(i));
+
+ for (Int_t i=0; i<fgMaxParams; ++i)
+ {
+ Double_t value = fgFitFunction->Eval(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");
+ }
+}
//
// 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)
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&);
#pragma link C++ class AliPWG0Helper+;
#pragma link C++ class AliCorrection+;
+#pragma link C++ class AliUnfolding+;
#pragma link C++ class AliMultiplicityCorrection+;
AliCorrectionMatrix3D.h \
AliCorrection.h \
AliPWG0Helper.h \
+ AliUnfolding.h \
multiplicity/AliMultiplicityCorrection.h
SRCS = $(HDRS:.h=.cxx)
#include <TH2F.h>
#include <TH3F.h>
#include <TDirectory.h>
-#include <TVirtualFitter.h>
#include <TCanvas.h>
#include <TString.h>
#include <TF1.h>
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!)
}
TH1::AddDirectory(oldStatus);
+
+ AliUnfolding::SetNbins(120, 120);
+ AliUnfolding::SetSkipBinsBegin(1);
+ AliUnfolding::SetNormalizeInput(kTRUE);
}
//____________________________________________________________________
}
//____________________________________________________________________
-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<fgNParamsMinuit; ++i)
- {
- Double_t right = params[i];
- Double_t left = params[i-1];
-
- if (right != 0)
- {
- Double_t diff = 1 - left / right;
- chi2 += diff * diff;
- }
- }
-
- return chi2 / 100.0;
-}
-
-//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params)
-{
- // homogenity term for minuit fitting
- // pure function of the parameters
- // prefers linear function (pol1)
-
- Double_t chi2 = 0;
-
- // ignore the first bin here. on purpose...
- for (Int_t i=2+1; i<fgNParamsMinuit; ++i)
- {
- if (params[i-1] == 0)
- continue;
-
- Double_t right = params[i];
- Double_t middle = params[i-1];
- Double_t left = params[i-2];
-
- Double_t der1 = (right - middle);
- Double_t der2 = (middle - left);
-
- Double_t diff = (der1 - der2) / middle;
-
- chi2 += diff * diff;
- }
-
- return chi2;
-}
-
-//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationLog(TVectorD& params)
-{
- // homogenity term for minuit fitting
- // pure function of the parameters
- // prefers linear function (pol1)
-
- Double_t chi2 = 0;
-
- /*Float_t der[fgkMaxParams];
-
- for (Int_t i=0; i<fgkMaxParams-1; ++i)
- der[i] = params[i+1] - params[i];
-
- for (Int_t i=0; i<fgkMaxParams-6; ++i)
- {
- for (Int_t j=1; j<=5; ++j)
- {
- Double_t diff = der[i+j] - der[i];
- chi2 += diff * diff;
- }
- }*/
-
- // ignore the first bin here. on purpose...
- for (Int_t i=2+1; i<fgNParamsMinuit; ++i)
- {
- if (params[i-1] == 0)
- continue;
-
- Double_t right = log(params[i]);
- Double_t middle = log(params[i-1]);
- Double_t left = log(params[i-2]);
-
- Double_t der1 = (right - middle);
- Double_t der2 = (middle - left);
-
- Double_t diff = (der1 - der2) / middle;
-
- chi2 += diff * diff;
- }
-
- return chi2 * 100;
-}
-
-//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& params)
-{
- // homogenity term for minuit fitting
- // pure function of the parameters
- // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
- // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
-
- Double_t chi2 = 0;
-
- // ignore the first bin here. on purpose...
- for (Int_t i=3; i<fgNParamsMinuit; ++i)
- {
- Double_t right = params[i];
- Double_t middle = params[i-1];
- Double_t left = params[i-2];
-
- Double_t der1 = (right - middle);
- Double_t der2 = (middle - left);
-
- Double_t diff = (der1 - der2);
-
- chi2 += diff * diff;
- }
-
- return chi2 * 1e4;
-}
-
-//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params)
-{
- // homogenity term for minuit fitting
- // pure function of the parameters
- // calculates entropy, from
- // The method of reduced cross-entropy (M. Schmelling 1993)
-
- Double_t paramSum = 0;
- // ignore the first bin here. on purpose...
- for (Int_t i=1; i<fgNParamsMinuit; ++i)
- paramSum += params[i];
-
- Double_t chi2 = 0;
- for (Int_t i=1; i<fgNParamsMinuit; ++i)
- {
- //Double_t tmp = params[i] / paramSum;
- Double_t tmp = params[i];
- if (tmp > 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; i<fgkMaxParams; ++i)
- params2[i] = fgNBD->Eval(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<fgNParamsMinuit; ++i)
- paramsVector[i] = params[i] * params[i];
-
- // calculate penalty factor
- Double_t penaltyVal = 0;
- switch (fgRegularizationType)
- {
- case kNone: break;
- case kPol0: penaltyVal = RegularizationPol0(paramsVector); break;
- case kPol1: penaltyVal = RegularizationPol1(paramsVector); break;
- case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break;
- case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break;
- case kLog: penaltyVal = RegularizationLog(paramsVector); break;
- }
-
- //if (penaltyVal > 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<fgkMaxInput; ++i)
- measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
-
- //measGuessVector.Print();
-
- // reduce influence of first bin
- //copy(0) *= 0.01;
-
- // (Ad - m) W (Ad - m)
- // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
- // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
- Double_t chi2FromFit = measGuessVector * copy * 1e6;
-
- chi2 = chi2FromFit + penaltyVal;
-
- static Int_t callCount = 0;
- if ((callCount++ % 10000) == 0)
- printf("%d) %f %f --> %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; i<fgNParamsMinuit; ++i)
- paramsVector[i] = params[i] * params[i];
-
- // Ad
- TVectorD measGuessVector(fgkMaxInput);
- measGuessVector = (*fgCorrelationMatrix) * paramsVector;
-
- // Ad - m
- measGuessVector -= (*fgCurrentESDVector);
-
- TVectorD copy(measGuessVector);
- //copy.Print();
-
- // (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<fgkMaxInput; ++i)
- measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
-
- // (Ad - m) W (Ad - m)
- // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
- // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
- //Double_t chi2FromFit = measGuessVector * copy * 1e6;
-
- // draw residuals
- TH1* residuals = new TH1F("residuals", "residuals", fgkMaxInput, -0.5, fgkMaxInput - 0.5);
- for (Int_t i=0; i<fgkMaxInput; ++i)
- residuals->SetBinContent(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; i<fgNParamsMinuit; ++i)
- {
- if (params[i-1] == 0)
- continue;
-
- Double_t right = params[i];
- Double_t middle = params[i-1];
- Double_t left = params[i-2];
-
- Double_t der1 = (right - middle);
- Double_t der2 = (middle - left);
-
- Double_t diff = (der1 - der2) / middle;
-
- penalty->SetBinContent(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();
}
}
- 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;
}
//____________________________________________________________________
-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 <measured> using response from <correlation> and effiency <aEfficiency>
- // output is in <result>
- // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
- // if <check> 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; i<fgNParamsMinuit; i++)
- (*fgEntropyAPriori)[i] = 1;
-
- if (!initialConditions) {
- initialConditions = measured;
- } else {
- printf("Using different starting conditions...\n");
- //new TCanvas; initialConditions->DrawCopy();
+ initialConditions = (TH1*) fCurrentESD->Clone("initialConditions");
initialConditions->Scale(1.0 / initialConditions->Integral());
- }
-
- for (Int_t i=0; i<fgNParamsMinuit; i++)
- if (initialConditions->GetBinContent(i+1) > 0)
- (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
-
- Double_t* results = new Double_t[fgNParamsMinuit+1];
- for (Int_t i=0; i<fgNParamsMinuit; ++i)
- {
- results[i] = initialConditions->GetBinContent(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; i<fgkMaxInput; ++i)
- {
- (*fgCurrentESDVector)[i] = measured->GetBinContent(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; i<fgNParamsMinuit; ++i)
- {
- results[i] = minuit->GetParameter(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; i<fgkMaxParams; ++i)
- {
- (*fgCurrentESDVector)[i] = fCurrentESD->GetBinContent(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; i<fgkMaxParams; ++i)
- {
- printf("%d : %f\n", i, fgNBD->Eval(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]);
}
//____________________________________________________________________
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)
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);
//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)
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)
}
//____________________________________________________________________
-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
// initialize seed with current time
gRandom->SetSeed(0);
+
+ if (methodType == AliUnfolding::kChi2Minimization)
+ AliUnfolding::SetCreateOverflowBin(5);
+ AliUnfolding::SetUnfoldingMethod(methodType);
const Int_t kErrorIterations = 150;
{
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");
}
// 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)
{
{
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;
// 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)
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)
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());
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; m<kMaxM; m++)
- {
- measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
- measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
-
- for (Int_t t=0; t<kMaxT; t++)
- {
- response[t][m] = correlation->GetBinContent(t+1, m+1);
- inverseResponse[t][m] = 0;
- }
- }
-
- for (Int_t t=0; t<kMaxT; t++)
- {
- efficiency[t] = aEfficiency->GetBinContent(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; i<kMaxT; i++)
- prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
- }
-
- //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
-
- // unfold...
- for (Int_t i=0; i<nIterations || nIterations < 0; i++)
- {
- //printf(" iteration %i \n", i);
-
- // calculate IR from Beyes theorem
- // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
-
- Double_t chi2Measured = 0;
- for (Int_t m=0; m<kMaxM; m++)
- {
- Float_t norm = 0;
- for (Int_t t = kStartBin; t<kMaxT; t++)
- norm += response[t][m] * prior[t];
-
- // calc. chi2: (measured - response * prior) / error
- if (measuredError[m] > 0)
- {
- Double_t value = (measuredCopy[m] - norm) / measuredError[m];
- chi2Measured += value * value;
- }
-
- if (norm > 0)
- {
- for (Int_t t = kStartBin; t<kMaxT; t++)
- inverseResponse[t][m] = response[t][m] * prior[t] / norm;
- }
- else
- {
- for (Int_t t = kStartBin; t<kMaxT; t++)
- inverseResponse[t][m] = 0;
- }
- }
- //Printf("chi2Measured of the last prior is %e", chi2Measured);
-
- for (Int_t t = kStartBin; t<kMaxT; t++)
- {
- Float_t value = 0;
- for (Int_t m=0; m<kMaxM; m++)
- value += inverseResponse[t][m] * measuredCopy[m];
-
- if (efficiency[t] > 0)
- result[t] = value / efficiency[t];
- else
- result[t] = 0;
- }
-
- // draw intermediate result
- /*
- for (Int_t t=0; t<kMaxT; t++)
- aResult->SetBinContent(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<kMaxT; t++)
- {
- Float_t newValue = 0;
-
- // 0 bin excluded from smoothing
- if (t > kStartBin+2 && t<kMaxT-1)
- {
- Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
-
- // weight the average with the regularization parameter
- newValue = (1 - regPar) * result[t] + regPar * average;
- }
- else
- newValue = result[t];
-
- // calculate chi2 (change from last iteration)
- if (prior[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; t<kMaxT; t++)
- aResult->SetBinContent(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; k<xBins; k++)
- {
- printf("In Cov1: %d\n", k);
- for (Int_t l=0; l<yBins; l++)
- {
- cov1[k][l] = 0;
-
- // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
- for (Int_t j=0; j<yBins; j++)
- cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
- * (1.0 - vectorn[j] / nEvents);
-
- // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
- for (Int_t i=0; i<yBins; i++)
- for (Int_t j=0; j<yBins; j++)
- {
- if (i == j)
- continue;
- cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
- * vectorn[j] / nEvents;
- }
- }
- }
-
- printf("Cov1 finished\n");
-
- TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("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; r<yBins; r++)
- for (Int_t u=0; u<xBins; u++)
- for (Int_t s=0; s<yBins; s++)
- {
- if (r == s)
- covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(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; k<xBins; k++)
- for (Int_t l=0; l<yBins; l++)
- {
- cov2[k][l] = 0;
- printf("In Cov2: %d %d\n", k, l);
- for (Int_t i=0; i<yBins; i++)
- for (Int_t j=0; j<yBins; j++)
- {
- //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
- // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
- Float_t tmpCov = 0;
- for (Int_t r=0; r<yBins; r++)
- for (Int_t u=0; u<xBins; u++)
- for (Int_t s=0; s<yBins; s++)
- {
- if (hResponse->GetBinContent(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)
{
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());
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());
#include <TMatrixD.h>
#include <TVectorD.h>
#include <AliPWG0Helper.h>
+#include <AliUnfolding.h>
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();
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);
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)
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);
//mult->SetMultiplicityESDCorrected(histID, (TH1F*) mcCompare);
mult->ApplyMinuitFit(histID, fullPhaseSpace, eventType, kFALSE); //hist2->ProjectionY("mymchist"));
+ //mult->ApplyNBDFit(histID, fullPhaseSpace, eventType);
}
else
{
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff