#include "AliCFUnfolding.h"
#include "TMath.h"
#include "TAxis.h"
#include "TF1.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "TRandom3.h"
ClassImp(AliCFUnfolding)
AliCFUnfolding::AliCFUnfolding() :
TNamed(),
fResponseOrig(0x0),
fPriorOrig(0x0),
fEfficiencyOrig(0x0),
fMeasuredOrig(0x0),
fMaxNumIterations(0),
fNVariables(0),
fUseSmoothing(kFALSE),
fSmoothFunction(0x0),
fSmoothOption("iremn"),
fMaxConvergence(0),
fNRandomIterations(0),
fResponse(0x0),
fPrior(0x0),
fEfficiency(0x0),
fMeasured(0x0),
fInverseResponse(0x0),
fMeasuredEstimate(0x0),
fConditional(0x0),
fUnfolded(0x0),
fUnfoldedFinal(0x0),
fCoordinates2N(0x0),
fCoordinatesN_M(0x0),
fCoordinatesN_T(0x0),
fRandomResponse(0x0),
fRandomEfficiency(0x0),
fRandomMeasured(0x0),
fRandom3(0x0),
fDeltaUnfoldedP(0x0),
fDeltaUnfoldedN(0x0),
fNCalcCorrErrors(0),
fRandomSeed(0)
{
}
AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar,
const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior ,
Double_t maxConvergencePerDOF, UInt_t randomSeed, Int_t maxNumIterations
) :
TNamed(name,title),
fResponseOrig((THnSparse*)response->Clone()),
fPriorOrig(0x0),
fEfficiencyOrig((THnSparse*)efficiency->Clone()),
fMeasuredOrig((THnSparse*)measured->Clone()),
fMaxNumIterations(maxNumIterations),
fNVariables(nVar),
fUseSmoothing(kFALSE),
fSmoothFunction(0x0),
fSmoothOption("iremn"),
fMaxConvergence(0),
fNRandomIterations(maxNumIterations),
fResponse((THnSparse*)response->Clone()),
fPrior(0x0),
fEfficiency((THnSparse*)efficiency->Clone()),
fMeasured((THnSparse*)measured->Clone()),
fInverseResponse(0x0),
fMeasuredEstimate(0x0),
fConditional(0x0),
fUnfolded(0x0),
fUnfoldedFinal(0x0),
fCoordinates2N(0x0),
fCoordinatesN_M(0x0),
fCoordinatesN_T(0x0),
fRandomResponse((THnSparse*)response->Clone()),
fRandomEfficiency((THnSparse*)efficiency->Clone()),
fRandomMeasured((THnSparse*)measured->Clone()),
fRandom3(0x0),
fDeltaUnfoldedP(0x0),
fDeltaUnfoldedN(0x0),
fNCalcCorrErrors(0),
fRandomSeed(randomSeed)
{
AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions()));
if (!prior) CreateFlatPrior();
else {
fPrior = (THnSparse*) prior->Clone();
fPriorOrig = (THnSparse*)fPrior->Clone();
if (fPrior->GetNdimensions() != fNVariables)
AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
}
if (fEfficiency->GetNdimensions() != fNVariables)
AliFatal(Form("The efficiency matrix should have %d dimensions, and it has actually %d",fNVariables,fEfficiency->GetNdimensions()));
if (fMeasured->GetNdimensions() != fNVariables)
AliFatal(Form("The measured matrix should have %d dimensions, and it has actually %d",fNVariables,fMeasured->GetNdimensions()));
if (fResponse->GetNdimensions() != 2*fNVariables)
AliFatal(Form("The response matrix should have %d dimensions, and it has actually %d",2*fNVariables,fResponse->GetNdimensions()));
for (Int_t iVar=0; iVar<fNVariables; iVar++) {
AliInfo(Form("prior matrix has %d bins in dimension %d",fPrior ->GetAxis(iVar)->GetNbins(),iVar));
AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar));
AliInfo(Form("measured matrix has %d bins in dimension %d",fMeasured ->GetAxis(iVar)->GetNbins(),iVar));
}
fRandomResponse ->SetTitle("Randomized response matrix");
fRandomEfficiency->SetTitle("Randomized efficiency");
fRandomMeasured ->SetTitle("Randomized measured");
SetMaxConvergencePerDOF(maxConvergencePerDOF) ;
Init();
}
AliCFUnfolding::~AliCFUnfolding() {
if (fResponse) delete fResponse;
if (fPrior) delete fPrior;
if (fEfficiency) delete fEfficiency;
if (fEfficiencyOrig) delete fEfficiencyOrig;
if (fMeasured) delete fMeasured;
if (fMeasuredOrig) delete fMeasuredOrig;
if (fSmoothFunction) delete fSmoothFunction;
if (fPriorOrig) delete fPriorOrig;
if (fInverseResponse) delete fInverseResponse;
if (fMeasuredEstimate) delete fMeasuredEstimate;
if (fConditional) delete fConditional;
if (fUnfolded) delete fUnfolded;
if (fUnfoldedFinal) delete fUnfoldedFinal;
if (fCoordinates2N) delete [] fCoordinates2N;
if (fCoordinatesN_M) delete [] fCoordinatesN_M;
if (fCoordinatesN_T) delete [] fCoordinatesN_T;
if (fRandomResponse) delete fRandomResponse;
if (fRandomEfficiency) delete fRandomEfficiency;
if (fRandomMeasured) delete fRandomMeasured;
if (fRandom3) delete fRandom3;
if (fDeltaUnfoldedP) delete fDeltaUnfoldedP;
if (fDeltaUnfoldedN) delete fDeltaUnfoldedN;
}
void AliCFUnfolding::Init() {
fRandom3 = new TRandom3(fRandomSeed);
fCoordinates2N = new Int_t[2*fNVariables];
fCoordinatesN_M = new Int_t[fNVariables];
fCoordinatesN_T = new Int_t[fNVariables];
CreateConditional();
fInverseResponse = (THnSparse*) fResponse->Clone();
fUnfolded = (THnSparse*) fPrior->Clone();
fUnfolded->SetTitle("Unfolded");
fMeasuredEstimate = (THnSparse*) fMeasured->Clone();
fDeltaUnfoldedP = (THnSparse*)fPrior->Clone();
fDeltaUnfoldedP->SetTitle("#Delta unfolded");
fDeltaUnfoldedP->Reset();
fDeltaUnfoldedN = (THnSparse*)fPrior->Clone();
fDeltaUnfoldedN->SetTitle("");
fDeltaUnfoldedN->Reset();
}
void AliCFUnfolding::CreateEstMeasured() {
fMeasuredEstimate->Reset();
THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
priorTimesEff->Multiply(fEfficiency);
for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
GetCoordinates();
Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T);
Double_t fill = conditionalValue * priorTimesEffValue ;
if (fill>0.) {
fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill);
fMeasuredEstimate->SetBinError(fCoordinatesN_M,0.);
}
}
delete priorTimesEff ;
}
void AliCFUnfolding::CreateInvResponse() {
THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
priorTimesEff->Multiply(fEfficiency);
for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
GetCoordinates();
Double_t estMeasuredValue = fMeasuredEstimate->GetBinContent(fCoordinatesN_M);
Double_t priorTimesEffValue = priorTimesEff ->GetBinContent(fCoordinatesN_T);
Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ;
if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) {
fInverseResponse->SetBinContent(fCoordinates2N,fill);
fInverseResponse->SetBinError (fCoordinates2N,0.);
}
}
delete priorTimesEff ;
}
void AliCFUnfolding::Unfold() {
Int_t iIterBayes = 0 ;
Double_t convergence = 0.;
for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) {
CreateEstMeasured();
CreateInvResponse();
CreateUnfolded();
convergence = GetConvergence();
AliDebug(0,Form("convergence at iteration %d is %e",iIterBayes,convergence));
if (fMaxConvergence>0. && convergence<fMaxConvergence && fNCalcCorrErrors == 0) {
fNRandomIterations = iIterBayes;
AliDebug(0,Form("convergence is met at iteration %d",iIterBayes));
break;
}
if (fUseSmoothing) {
if (Smooth()) {
AliError("Couldn't smooth the unfolded spectrum!!");
if (fNCalcCorrErrors>0) {
AliInfo(Form("=======================\nUnfold of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
}
else {
AliInfo(Form("\n\n=======================\nFinish at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
}
return;
}
}
if (fPrior) delete fPrior ;
fPrior = (THnSparse*)fUnfolded->Clone() ;
fPrior->SetTitle("Prior");
}
if (fNCalcCorrErrors==0) fUnfoldedFinal = (THnSparse*) fUnfolded->Clone() ;
if (fNCalcCorrErrors == 0) {
AliInfo("\n================================================\nFinished bayes iteration, now calculating errors...\n================================================\n");
fNCalcCorrErrors = 1;
CalculateCorrelatedErrors();
}
if (fNCalcCorrErrors >1 ) {
AliInfo(Form("\n\n=======================\nFinished at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
}
else if(fNCalcCorrErrors>0) {
AliInfo(Form("=======================\nUnfolding of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
}
}
void AliCFUnfolding::CreateUnfolded() {
fUnfolded->Reset();
for (Long_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) {
Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N);
GetCoordinates();
Double_t effValue = fEfficiency->GetBinContent(fCoordinatesN_T);
Double_t measuredValue = fMeasured ->GetBinContent(fCoordinatesN_M);
Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ;
if (fill>0.) {
Double_t err = 0.;
fUnfolded->SetBinError (fCoordinatesN_T,err);
fUnfolded->AddBinContent(fCoordinatesN_T,fill);
}
}
}
void AliCFUnfolding::CalculateCorrelatedErrors() {
for (int i=0; i<fNRandomIterations; i++) {
if (fPrior) delete fPrior ;
fPrior = (THnSparse*) fPriorOrig->Clone();
CreateRandomizedDist();
if (fResponse) delete fResponse ;
fResponse = (THnSparse*) fRandomResponse->Clone();
fResponse->SetTitle("Response");
if (fEfficiency) delete fEfficiency ;
fEfficiency = (THnSparse*) fRandomEfficiency->Clone();
fEfficiency->SetTitle("Efficiency");
if (fMeasured) delete fMeasured ;
fMeasured = (THnSparse*) fRandomMeasured->Clone();
fMeasured->SetTitle("Measured");
Unfold();
FillDeltaUnfoldedProfile();
}
Double_t meanx2 = 0.;
Double_t mean = 0.;
Double_t checksigma = 0.;
Double_t entriesInBin = 0.;
for (Long_t iBin=0; iBin<fUnfoldedFinal->GetNbins(); iBin++) {
fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M);
mean = fDeltaUnfoldedP->GetBinContent(fCoordinatesN_M);
meanx2 = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M);
entriesInBin = fDeltaUnfoldedN->GetBinContent(fCoordinatesN_M);
if(entriesInBin > 1.) checksigma = TMath::Sqrt((entriesInBin/(entriesInBin-1.))*TMath::Abs(meanx2-mean*mean));
fUnfoldedFinal->SetBinError(fCoordinatesN_M,checksigma);
}
fNCalcCorrErrors = 2;
}
void AliCFUnfolding::CreateRandomizedDist() {
for (Long_t iBin=0; iBin<fResponseOrig->GetNbins(); iBin++) {
Double_t val = fResponseOrig->GetBinContent(iBin,fCoordinatesN_M);
Double_t err = fResponseOrig->GetBinError(fCoordinatesN_M);
Double_t ran = fRandom3->Gaus(val,err);
fRandomResponse->SetBinContent(iBin,ran);
}
for (Long_t iBin=0; iBin<fEfficiencyOrig->GetNbins(); iBin++) {
Double_t val = fEfficiencyOrig->GetBinContent(iBin,fCoordinatesN_M);
Double_t err = fEfficiencyOrig->GetBinError(fCoordinatesN_M);
Double_t ran = fRandom3->Gaus(val,err);
fRandomEfficiency->SetBinContent(iBin,ran);
}
for (Long_t iBin=0; iBin<fMeasuredOrig->GetNbins(); iBin++) {
Double_t val = fMeasuredOrig->GetBinContent(iBin,fCoordinatesN_M);
Double_t err = fMeasuredOrig->GetBinError(fCoordinatesN_M);
Double_t ran = fRandom3->Gaus(val,err);
fRandomMeasured->SetBinContent(iBin,ran);
}
}
void AliCFUnfolding::FillDeltaUnfoldedProfile() {
for (Long_t iBin=0; iBin<fUnfoldedFinal->GetNbins(); iBin++) {
Double_t deltaInBin = fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M) - fUnfolded->GetBinContent(fCoordinatesN_M);
Double_t entriesInBin = fDeltaUnfoldedN->GetBinContent(fCoordinatesN_M);
Double_t mean_n = fDeltaUnfoldedP->GetBinContent(fCoordinatesN_M) ;
Double_t mean_nplus1 = mean_n ;
mean_nplus1 *= entriesInBin ;
mean_nplus1 += deltaInBin ;
mean_nplus1 /= (entriesInBin+1) ;
Double_t meanx2_n = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M) ;
Double_t meanx2_nplus1 = meanx2_n ;
meanx2_nplus1 *= entriesInBin ;
meanx2_nplus1 += (deltaInBin*deltaInBin) ;
meanx2_nplus1 /= (entriesInBin+1) ;
fDeltaUnfoldedP->SetBinError(fCoordinatesN_M,meanx2_nplus1) ;
fDeltaUnfoldedP->SetBinContent(fCoordinatesN_M,mean_nplus1) ;
fDeltaUnfoldedN->SetBinContent(fCoordinatesN_M,entriesInBin+1);
}
}
void AliCFUnfolding::GetCoordinates() {
for (Int_t i = 0; i<fNVariables ; i++) {
fCoordinatesN_M[i] = fCoordinates2N[i];
fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables];
}
}
void AliCFUnfolding::CreateConditional() {
fConditional = (THnSparse*) fResponse->Clone();
Int_t* dim = new Int_t [fNVariables];
for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ;
THnSparse* responseInT = fConditional->Projection(fNVariables,dim,"E");
delete [] dim;
for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
GetCoordinates();
Double_t projValue = responseInT->GetBinContent(fCoordinatesN_T);
Double_t fill = responseValue / projValue ;
if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) {
fConditional->SetBinContent(fCoordinates2N,fill);
Double_t err = 0.;
fConditional->SetBinError (fCoordinates2N,err);
}
}
delete responseInT ;
}
Int_t AliCFUnfolding::GetDOF() {
Int_t nDOF = 1 ;
for (Int_t iDim=0; iDim<fNVariables; iDim++) {
nDOF *= fPrior->GetAxis(iDim)->GetNbins();
}
AliDebug(0,Form("Number of degrees of freedom = %d",nDOF));
return nDOF;
}
Double_t AliCFUnfolding::GetChi2() {
Double_t chi2 = 0. ;
Double_t error_unf = 0.;
for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
Double_t priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
error_unf = fUnfolded->GetBinError(fCoordinatesN_T);
chi2 += (error_unf > 0. ? TMath::Power((fUnfolded->GetBinContent(fCoordinatesN_T) - priorValue)/error_unf,2) / priorValue : 0.) ;
}
return chi2;
}
Double_t AliCFUnfolding::GetConvergence() {
Double_t convergence = 0.;
Double_t priorValue = 0.;
Double_t currentValue = 0.;
for (Long_t iBin=0; iBin < fPrior->GetNbins(); iBin++) {
priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
currentValue = fUnfolded->GetBinContent(fCoordinatesN_T);
if (priorValue > 0.)
convergence += ((priorValue-currentValue)/priorValue)*((priorValue-currentValue)/priorValue);
else
AliWarning(Form("priorValue = %f. Adding 0 to convergence criterion.",priorValue));
}
return convergence;
}
void AliCFUnfolding::SetMaxConvergencePerDOF(Double_t val) {
Int_t nDOF = GetDOF() ;
fMaxConvergence = val * nDOF ;
AliInfo(Form("MaxConvergence = %e. Number of degrees of freedom = %d",fMaxConvergence,nDOF));
}
Short_t AliCFUnfolding::Smooth() {
if (fSmoothFunction) {
AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction));
return SmoothUsingFunction();
}
else return SmoothUsingNeighbours(fUnfolded);
}
Short_t AliCFUnfolding::SmoothUsingNeighbours(THnSparse* hist) {
Int_t const nDimensions = hist->GetNdimensions() ;
Int_t* coordinates = new Int_t[nDimensions];
Int_t* numBins = new Int_t[nDimensions];
for (Int_t iVar=0; iVar<nDimensions; iVar++) numBins[iVar] = hist->GetAxis(iVar)->GetNbins();
THnSparse* copy = (THnSparse*)hist->Clone();
for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) {
Double_t content = copy->GetBinContent(iBin,coordinates);
Double_t error2 = TMath::Power(copy->GetBinError(iBin),2);
Bool_t isOutside = kFALSE ;
for (Int_t iVar=0; iVar<nDimensions; iVar++) {
if (coordinates[iVar]<1 || coordinates[iVar]>numBins[iVar]) {
isOutside=kTRUE;
break;
}
}
if (isOutside) continue;
Int_t neighbours = 0;
for (Int_t iVar=0; iVar<nDimensions; iVar++) {
if (coordinates[iVar] > 1) {
coordinates[iVar]-- ;
content += copy->GetBinContent(coordinates);
error2 += TMath::Power(copy->GetBinError(coordinates),2);
neighbours++;
coordinates[iVar]++ ;
}
if (coordinates[iVar] < numBins[iVar]) {
coordinates[iVar]++ ;
content += copy->GetBinContent(coordinates);
error2 += TMath::Power(copy->GetBinError(coordinates),2);
neighbours++;
coordinates[iVar]-- ;
}
}
hist->SetBinContent(coordinates,content/(1.+neighbours));
hist->SetBinError (coordinates,TMath::Sqrt(error2)/(1.+neighbours));
}
delete [] numBins;
delete [] coordinates ;
delete copy;
return 0;
}
Short_t AliCFUnfolding::SmoothUsingFunction() {
AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar()));
for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar)));
Int_t fitResult = 0;
switch (fNVariables) {
case 1 : fitResult = fUnfolded->Projection(0) ->Fit(fSmoothFunction,fSmoothOption); break;
case 2 : fitResult = fUnfolded->Projection(1,0) ->Fit(fSmoothFunction,fSmoothOption); break;
case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break;
default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1;
}
if (fitResult != 0) {
AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult));
return 1;
}
Int_t nDim = fNVariables;
Int_t* bins = new Int_t[nDim];
Long_t nBins = 1;
for (Int_t iVar=0; iVar<nDim; iVar++) {
bins[iVar] = fUnfolded->GetAxis(iVar)->GetNbins();
nBins *= bins[iVar];
}
Int_t *bin = new Int_t[nDim];
Double_t x[3] = {0,0,0} ;
for (Long_t iBin=0; iBin<nBins; iBin++) {
Long_t bin_tmp = iBin ;
for (Int_t iVar=0; iVar<nDim; iVar++) {
bin[iVar] = 1 + bin_tmp % bins[iVar] ;
bin_tmp /= bins[iVar] ;
x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]);
}
Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ;
fUnfolded->SetBinError (bin,fUnfolded->GetBinError(bin)*functionValue/fUnfolded->GetBinContent(bin));
fUnfolded->SetBinContent(bin,functionValue);
}
delete [] bins;
delete [] bin ;
return 0;
}
void AliCFUnfolding::CreateFlatPrior() {
AliInfo("Creating a flat a priori distribution");
fPrior = (THnSparse*) fEfficiency->Clone();
fPrior->SetTitle("Prior");
if (fNVariables != fPrior->GetNdimensions())
AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
Int_t nDim = fNVariables;
Int_t* bins = new Int_t[nDim];
Long_t nBins = 1;
for (Int_t iVar=0; iVar<nDim; iVar++) {
bins[iVar] = fPrior->GetAxis(iVar)->GetNbins();
nBins *= bins[iVar];
}
Int_t *bin = new Int_t[nDim];
for (Long_t iBin=0; iBin<nBins; iBin++) {
Long_t bin_tmp = iBin ;
for (Int_t iVar=0; iVar<nDim; iVar++) {
bin[iVar] = 1 + bin_tmp % bins[iVar] ;
bin_tmp /= bins[iVar] ;
}
fPrior->SetBinContent(bin,1.);
fPrior->SetBinError (bin,0.);
}
fPriorOrig = (THnSparse*)fPrior->Clone();
delete [] bin;
delete [] bins;
}