#include <TF1.h>
TMatrixD* AliUnfolding::fgCorrelationMatrix = 0;
+TMatrixD* AliUnfolding::fgCorrelationMatrixSquared = 0;
TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0;
TVectorD* AliUnfolding::fgCurrentESDVector = 0;
TVectorD* AliUnfolding::fgEntropyAPriori = 0;
+TVectorD* AliUnfolding::fgEfficiency = 0;
+TVectorD* AliUnfolding::fgBinWidths = 0;
TF1* AliUnfolding::fgFitFunction = 0;
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::fgMinimumInitialValue = kFALSE; // set all initial values at least to the smallest value among the initial values
+Float_t AliUnfolding::fgMinimumInitialValueFix = -1;
Bool_t AliUnfolding::fgNormalizeInput = kFALSE; // normalize input spectrum
+Float_t AliUnfolding::fgNotFoundEvents = 0;
+Bool_t AliUnfolding::fgSkipBin0InChi2 = kFALSE;
Float_t AliUnfolding::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing)
-Int_t AliUnfolding::fgBayesianIterations = 10; // number of iterations in Bayesian method
+Int_t AliUnfolding::fgBayesianIterations = 10; // number of iterations in Bayesian method
Bool_t AliUnfolding::fgDebug = kFALSE;
+Int_t AliUnfolding::fgCallCount = 0;
+
ClassImp(AliUnfolding)
//____________________________________________________________________
fgMaxInput = nMeasured;
fgMaxParams = nUnfolded;
+ if (fgCorrelationMatrix)
+ {
+ delete fgCorrelationMatrix;
+ fgCorrelationMatrix = 0;
+ }
+ if (fgCorrelationMatrixSquared)
+ {
+ fgCorrelationMatrixSquared = 0;
+ delete fgCorrelationMatrixSquared;
+ }
+ if (fgCorrelationCovarianceMatrix)
+ {
+ delete fgCorrelationCovarianceMatrix;
+ fgCorrelationCovarianceMatrix = 0;
+ }
+ if (fgCurrentESDVector)
+ {
+ delete fgCurrentESDVector;
+ fgCurrentESDVector = 0;
+ }
+ if (fgEntropyAPriori)
+ {
+ delete fgEntropyAPriori;
+ fgEntropyAPriori = 0;
+ }
+ if (fgEfficiency)
+ {
+ delete fgEfficiency;
+ fgEfficiency = 0;
+ }
+ if (fgBinWidths)
+ {
+ delete fgBinWidths;
+ fgBinWidths = 0;
+ }
+
Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded);
}
// 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
+ //
+ // return code: see UnfoldWithMinuit/UnfoldWithBayesian/UnfoldWithFunction
if (fgMaxInput == -1)
{
}
//____________________________________________________________________
-void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured)
+void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured, TH1* efficiency)
{
// fill static variables needed for minuit fit
if (!fgCorrelationMatrix)
fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
+ if (!fgCorrelationMatrixSquared)
+ fgCorrelationMatrixSquared = new TMatrixD(fgMaxInput, fgMaxParams);
if (!fgCorrelationCovarianceMatrix)
fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
if (!fgCurrentESDVector)
fgCurrentESDVector = new TVectorD(fgMaxInput);
if (!fgEntropyAPriori)
fgEntropyAPriori = new TVectorD(fgMaxParams);
-
+ if (!fgEfficiency)
+ fgEfficiency = new TVectorD(fgMaxParams);
+ if (!fgBinWidths)
+ fgBinWidths = new TVectorD(fgMaxParams);
+
fgCorrelationMatrix->Zero();
fgCorrelationCovarianceMatrix->Zero();
fgCurrentESDVector->Zero();
}
// npart sum to 1
- correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
+ correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);// * correlation->GetXaxis()->GetBinWidth(i));
correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
if (i <= fgMaxParams && j <= fgMaxInput)
+ {
(*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
+ (*fgCorrelationMatrixSquared)(j-1, i-1) = correlation->GetBinContent(i, j) * correlation->GetBinContent(i, j);
+ }
}
//printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
}
//normalize measured
+ Float_t smallestError = 1;
if (fgNormalizeInput)
- measured->Scale(1.0 / measured->Integral());
+ {
+ Float_t sumMeasured = measured->Integral();
+ measured->Scale(1.0 / sumMeasured);
+ smallestError /= sumMeasured;
+ }
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);
+ }
+ else // in this case put error of 1, otherwise 0 bins are not added to the chi2...
+ (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / smallestError / smallestError;
if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
(*fgCorrelationCovarianceMatrix)(i, i) = 0;
//Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i));
}
+
+ // efficiency is expected to match bin width of result
+ for (Int_t i=0; i<fgMaxParams; ++i)
+ {
+ (*fgEfficiency)(i) = efficiency->GetBinContent(i+1);
+ (*fgBinWidths)(i) = efficiency->GetXaxis()->GetBinWidth(i+1);
+ }
}
//____________________________________________________________________
// 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
+ // negative values in initialConditions mean that the given parameter is fixed to the absolute of the value
// 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);
+ SetStaticVariables(correlation, measured, efficiency);
// Initialize TMinuit via generic fitter interface
- TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
+ Int_t params = fgMaxParams;
+ if (fgNotFoundEvents > 0)
+ params++;
+
+ TVirtualFitter *minuit = TVirtualFitter::Fitter(0, params);
Double_t arglist[100];
// disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
for (Int_t i=0; i<fgMaxParams; i++)
(*fgEntropyAPriori)[i] = 1;
+ // set initial conditions as a-priori distribution for MRX regularization
+ /*
+ for (Int_t i=0; i<fgMaxParams; i++)
+ if (initialConditions && initialConditions->GetBinContent(i+1) > 0)
+ (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
+ */
+
if (!initialConditions) {
initialConditions = measured;
} else {
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);
-
+ // extract minimum value from initial conditions (if we set a value to 0 it will stay 0)
+ Float_t minValue = 1e100;
+ if (fgMinimumInitialValueFix < 0)
+ {
+ for (Int_t i=0; i<fgMaxParams; ++i)
+ {
+ Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
+ if (initialConditions->GetBinContent(bin) > 0)
+ minValue = TMath::Min(minValue, (Float_t) initialConditions->GetBinContent(bin));
+ }
+ }
+ else
+ minValue = fgMinimumInitialValueFix;
+
Double_t* results = new Double_t[fgMaxParams+1];
for (Int_t i=0; i<fgMaxParams; ++i)
{
- results[i] = initialConditions->GetBinContent(i+1);
+ Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
+ results[i] = initialConditions->GetBinContent(bin);
+ Bool_t fix = kFALSE;
+ if (results[i] < 0)
+ {
+ fix = kTRUE;
+ results[i] = -results[i];
+ }
+
+ if (!fix && fgMinimumInitialValue && results[i] < minValue)
+ results[i] = minValue;
+
// 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);
+ minuit->SetParameter(i, Form("param%d", i), results[i], (fix) ? 0 : fgMinuitStepSize, 0, 0);
+ }
+ if (fgNotFoundEvents > 0)
+ {
+ results[fgMaxParams] = efficiency->GetBinContent(1);
+ minuit->SetParameter(fgMaxParams, "vtx0", results[fgMaxParams], fgMinuitStepSize / 100, 0.01, 0.80);
}
Int_t dummy = 0;
//printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
//minuit->ExecuteCommand("MINOS", arglist, 0);
+ if (fgNotFoundEvents > 0)
+ {
+ results[fgMaxParams] = minuit->GetParameter(fgMaxParams);
+ Printf("Efficiency for bin 0 changed from %f to %f", efficiency->GetBinContent(1), results[fgMaxParams]);
+ efficiency->SetBinContent(1, results[fgMaxParams]);
+ }
+
for (Int_t i=0; i<fgMaxParams; ++i)
{
results[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);
+ value /= efficiency->GetBinContent(i+1);
+ error /= efficiency->GetBinContent(i+1);
}
else
{
- value = 0;
- error = 0;
+ value = 0;
+ error = 0;
}
}
result->SetBinError(i+1, error);
}
+ fgCallCount = 0;
+ Chi2Function(dummy, 0, chi2, results, 0);
+ printf("AliUnfolding::UnfoldWithMinuit: Chi2 of final parameters is = %f\n", chi2);
+
if (fgDebug)
DrawGuess(results);
Double_t* prior = new Double_t[kMaxT];
Double_t* result = new Double_t[kMaxT];
Double_t* efficiency = new Double_t[kMaxT];
+ Double_t* binWidths = new Double_t[kMaxT];
Double_t** response = new Double_t*[kMaxT];
Double_t** inverseResponse = new Double_t*[kMaxT];
for (Int_t t=0; t<kMaxT; t++)
{
- if (efficiency)
+ if (aEfficiency)
{
efficiency[t] = aEfficiency->GetBinContent(t+1);
}
prior[t] = measuredCopy[t];
result[t] = 0;
+ binWidths[t] = aResult->GetXaxis()->GetBinWidth(t+1);
}
// pick prior distribution
//TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
+ //new TCanvas;
// unfold...
- for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations< 0; i++)
+ for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations < 0; i++)
{
if (fgDebug)
Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i);
- // calculate IR from Beyes theorem
+ // calculate IR from Bayes theorem
// IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
Double_t chi2Measured = 0;
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->SetMarkerStyle(24+i);
aResult->SetMarkerColor(2);
- aResult->DrawCopy("P SAME HIST");
+ aResult->DrawCopy((i == 0) ? "P" : "PSAME");
*/
-
+
Double_t chi2LastIter = 0;
// regularization (simple smoothing)
for (Int_t t=kStartBin; t<kMaxT; t++)
// 0 bin excluded from smoothing
if (t > kStartBin+2 && t<kMaxT-1)
{
- Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
+ Float_t average = (result[t-1] / binWidths[t-1] + result[t] / binWidths[t] + result[t+1] / binWidths[t+1]) / 3 * binWidths[t];
// weight the average with the regularization parameter
newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average;
//new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
//delete convergence;
+ Float_t factor = 1;
+ if (!fgNormalizeInput)
+ factor = measuredIntegral;
for (Int_t t=0; t<kMaxT; t++)
- aResult->SetBinContent(t+1, result[t]);
+ aResult->SetBinContent(t+1, result[t] * factor);
delete[] measuredCopy;
delete[] measuredError;
delete[] prior;
delete[] result;
delete[] efficiency;
+ delete[] binWidths;
for (Int_t i=0; i<kMaxT; i++)
{
for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
{
- Double_t right = params[i];
- Double_t left = params[i-1];
+ Double_t right = params[i] / (*fgBinWidths)(i);
+ Double_t left = params[i-1] / (*fgBinWidths)(i-1);
- if (right != 0)
+ if (left != 0)
{
- Double_t diff = 1 - left / right;
- chi2 += diff * diff;
+ Double_t diff = (right - left);
+ chi2 += diff * diff / left / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
}
}
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 right = params[i] / (*fgBinWidths)(i);
+ Double_t middle = params[i-1] / (*fgBinWidths)(i-1);
+ Double_t left = params[i-2] / (*fgBinWidths)(i-2);
- Double_t diff = (der1 - der2) / middle;
+ Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
+ Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 2);
- chi2 += diff * diff;
+ //Double_t diff = (der1 - der2) / middle;
+ //chi2 += diff * diff;
+ chi2 += (der1 - der2) * (der1 - der2) / middle * (*fgBinWidths)(i-1);
}
return chi2;
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);
+ if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0)
+ continue;
- Double_t diff = (der1 - der2) / middle;
+ Double_t right = log(params[i] / (*fgBinWidths)(i));
+ Double_t middle = log(params[i-1] / (*fgBinWidths)(i-1));
+ Double_t left = log(params[i-2] / (*fgBinWidths)(i-2));
+
+ Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
+ Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 2);
+
+ //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
- chi2 += diff * diff;
+ //if (fgCallCount == 0)
+ // Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
+ chi2 += (der1 - der2) * (der1 - der2);// / error;
}
- return chi2 * 100;
+ return chi2;
}
//____________________________________________________________________
Double_t chi2 = 0;
for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
{
- //Double_t tmp = params[i] / paramSum;
- Double_t tmp = params[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]);
}
+ else
+ chi2 += 100;
+ }
+
+ return -chi2;
+}
+
+//____________________________________________________________________
+Double_t AliUnfolding::RegularizationRatio(TVectorD& params)
+{
+ // homogenity term for minuit fitting
+ // pure function of the parameters
+
+ Double_t chi2 = 0;
+
+ for (Int_t i=5+fgSkipBinsBegin; i<fgMaxParams; ++i)
+ {
+ if (params[i-1] == 0 || params[i] == 0)
+ continue;
+
+ Double_t right = params[i] / (*fgBinWidths)(i);
+ Double_t middle = params[i-1] / (*fgBinWidths)(i-1);
+ Double_t left = params[i-2] / (*fgBinWidths)(i-2);
+ Double_t left2 = params[i-3] / (*fgBinWidths)(i-3);
+ Double_t left3 = params[i-4] / (*fgBinWidths)(i-4);
+ Double_t left4 = params[i-5] / (*fgBinWidths)(i-5);
+
+ //Double_t diff = left / middle - middle / right;
+ //Double_t diff = 2 * left / middle - middle / right - left2 / left;
+ Double_t diff = 4 * left2 / left - middle / right - left / middle - left3 / left2 - left4 / left3;
+
+ chi2 += diff * diff;// / middle;
}
- return 100.0 + chi2;
+ return chi2;
}
//____________________________________________________________________
// fit function for minuit
// does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
//
+
+ // TODO use static members for the variables here to speed up processing (no construction/deconstruction)
// d
- static TVectorD paramsVector(fgMaxParams);
+ TVectorD paramsVector(fgMaxParams);
for (Int_t i=0; i<fgMaxParams; ++i)
paramsVector[i] = params[i] * params[i];
case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break;
case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break;
case kLog: penaltyVal = RegularizationLog(paramsVector); break;
+ case kRatio: penaltyVal = RegularizationRatio(paramsVector); break;
}
//if (penaltyVal > 1e10)
// Ad - m
measGuessVector -= (*fgCurrentESDVector);
+
+#if 0
+ // new error calcuation using error on the guess
+
+ // error from guess
+ TVectorD errorGuessVector(fgMaxInput);
+ //errorGuessVector = (*fgCorrelationMatrixSquared) * paramsVector;
+ errorGuessVector = (*fgCorrelationMatrix) * paramsVector;
+ Double_t chi2FromFit = 0;
+ for (Int_t i=0; i<fgMaxInput; ++i)
+ {
+ Float_t error = 1;
+ if (errorGuessVector(i) > 0)
+ error = errorGuessVector(i);
+ chi2FromFit += measGuessVector(i) * measGuessVector(i) / error;
+ }
//measGuessVector.Print();
+#else
+ // old error calcuation using the error on the measured
TVectorD copy(measGuessVector);
// (Ad - m) W
//measGuessVector.Print();
+ if (fgSkipBin0InChi2)
+ measGuessVector[0] = 0;
+
// (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)
+ // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see UnfoldWithMinuit)
Double_t chi2FromFit = measGuessVector * copy * 1e6;
+#endif
- chi2 = chi2FromFit + penaltyVal;
+ Double_t notFoundEventsConstraint = 0;
+ Double_t currentNotFoundEvents = 0;
+ Double_t errorNotFoundEvents = 0;
+
+ if (fgNotFoundEvents > 0)
+ {
+ for (Int_t i=0; i<fgMaxParams; ++i)
+ {
+ Double_t eff = (1.0 / (*fgEfficiency)(i) - 1);
+ if (i == 0)
+ eff = (1.0 / params[fgMaxParams] - 1);
+ currentNotFoundEvents += eff * paramsVector(i);
+ errorNotFoundEvents += eff * eff * paramsVector(i); // error due to guess (paramsVector)
+ if (i <= 3)
+ errorNotFoundEvents += (eff * 0.03) * (eff * 0.03) * paramsVector(i) * paramsVector(i); // error on eff
+ //if ((fgCallCount % 10000) == 0)
+ // Printf("%d %f %f %f", i, (*fgEfficiency)(i), paramsVector(i), currentNotFoundEvents);
+ }
+ errorNotFoundEvents += fgNotFoundEvents;
+ // TODO add error on background, background estimate
+
+ notFoundEventsConstraint = (currentNotFoundEvents - fgNotFoundEvents) * (currentNotFoundEvents - fgNotFoundEvents) / errorNotFoundEvents;
+ }
+
+ Float_t avg = 0;
+ Float_t sum = 0;
+ Float_t currentMult = 0;
+ // try to match dn/deta
+ for (Int_t i=0; i<fgMaxParams; ++i)
+ {
+ avg += paramsVector(i) * currentMult;
+ sum += paramsVector(i);
+ currentMult += (*fgBinWidths)(i);
+ }
+ avg /= sum;
+ Float_t chi2Avg = 0; //(avg - 3.73) * (avg - 3.73) * 100;
- static Int_t callCount = 0;
- if ((callCount++ % 10000) == 0)
- Printf("AliUnfolding::Chi2Function: Iteration %d: %f %f --> %f", callCount, chi2FromFit, penaltyVal, chi2);
+ chi2 = chi2FromFit + penaltyVal + notFoundEventsConstraint + chi2Avg;
+
+ if ((fgCallCount++ % 1000) == 0)
+ {
+ Printf("AliUnfolding::Chi2Function: Iteration %d (ev %d %d +- %f) (%f) (%f): %f %f %f %f --> %f", fgCallCount-1, (Int_t) fgNotFoundEvents, (Int_t) currentNotFoundEvents, TMath::Sqrt(errorNotFoundEvents), params[fgMaxParams], avg, chi2FromFit, penaltyVal, notFoundEventsConstraint, chi2Avg, chi2);
+ //for (Int_t i=0; i<fgMaxInput; ++i)
+ // Printf("%d: %f", i, measGuessVector(i) * copy(i) * 1e6);
+ }
}
//____________________________________________________________________
//
// d
- static TVectorD paramsVector(fgMaxParams);
+ TVectorD paramsVector(fgMaxParams);
for (Int_t i=0; i<fgMaxParams; ++i)
paramsVector[i] = params[i] * params[i];
new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
// draw penalty
- if (fgRegularizationType != kPol1) {
- Printf("Drawing not supported");
- return;
- }
+ TH1* penalty = GetPenaltyPlot(params);
+
+ new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
+}
- TH1* penalty = new TH1F("penalty", "penalty", fgMaxParams, -0.5, fgMaxParams - 0.5);
+//____________________________________________________________________
+TH1* AliUnfolding::GetPenaltyPlot(TH1* corrected)
+{
+ // draws the penalty factors as function of multiplicity of the current selected regularization
- for (Int_t i=2+1; i<fgMaxParams; ++i)
+ Double_t* params = new Double_t[fgMaxParams];
+ for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++)
+ params[i] = corrected->GetBinContent(i+1);
+
+ TH1* penalty = GetPenaltyPlot(params);
+
+ delete[] params;
+
+ return penalty;
+}
+
+//____________________________________________________________________
+TH1* AliUnfolding::GetPenaltyPlot(Double_t* params)
+{
+ // draws the penalty factors as function of multiplicity of the current selected regularization
+
+ TH1* penalty = new TH1F("penalty", ";unfolded multiplicity;penalty factor", fgMaxParams, -0.5, fgMaxParams - 0.5);
+
+ for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
{
- if (params[i-1] == 0)
- continue;
+ Double_t diff = 0;
+ if (fgRegularizationType == kPol0)
+ {
+ Double_t right = params[i] / (*fgBinWidths)(i);
+ Double_t left = params[i-1] / (*fgBinWidths)(i-1);
- Double_t right = params[i];
- Double_t middle = params[i-1];
- Double_t left = params[i-2];
+ if (left != 0)
+ {
+ Double_t diffTmp = (right - left);
+ diff = diffTmp * diffTmp / left / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2) / 100;
+ }
+ }
+ if (fgRegularizationType == kPol1 && i > 1+fgSkipBinsBegin)
+ {
+ if (params[i-1] == 0)
+ continue;
- Double_t der1 = (right - middle);
- Double_t der2 = (middle - left);
+ Double_t right = params[i] / (*fgBinWidths)(i);
+ Double_t middle = params[i-1] / (*fgBinWidths)(i-1);
+ Double_t left = params[i-2] / (*fgBinWidths)(i-2);
- Double_t diff = (der1 - der2) / middle;
+ Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
+ Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 2);
+
+ diff = (der1 - der2) * (der1 - der2) / middle;
+ }
+ if (fgRegularizationType == kLog && i > 1+fgSkipBinsBegin)
+ {
+ 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 error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
+ //Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
+
+ diff = (der1 - der2) * (der1 - der2);// / error;
+ }
- penalty->SetBinContent(i-1, diff * diff);
+ penalty->SetBinContent(i, diff);
//Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
}
- new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
+
+ return penalty;
}
//____________________________________________________________________
SetChi2Regularization(kNone, 0);
- SetStaticVariables(correlation, measured);
+ SetStaticVariables(correlation, measured, efficiency);
// Initialize TMinuit via generic fitter interface
TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar());
correlation->DrawCopy("COLZ");
}
}
+
+Int_t AliUnfolding::UnfoldGetBias(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result)
+{
+ // unfolds and assigns bias as errors with Eq. 19 of Cowan, "a survey of unfolding methods for particle physics"
+ // b_i = sum_j dmu_i/dn_j (nu_j - n_j) with nu_j as folded guess, n_j as data
+ // dmu_i/dn_j is found numerically by changing the bin content and re-unfolding
+ //
+ // return code: 0 (success) -1 (error: from Unfold(...) )
+
+ if (Unfold(correlation, efficiency, measured, initialConditions, result) != 0)
+ return -1;
+
+ TMatrixD matrix(fgMaxInput, fgMaxParams);
+
+ TH1* newResult[4];
+ for (Int_t loop=0; loop<4; loop++)
+ newResult[loop] = (TH1*) result->Clone(Form("newresult_%d", loop));
+
+ // change bin-by-bin and built matrix of effects
+ for (Int_t m=0; m<fgMaxInput; m++)
+ {
+ if (measured->GetBinContent(m+1) < 1)
+ continue;
+
+ for (Int_t loop=0; loop<4; loop++)
+ {
+ //Printf("%d %d", i, loop);
+ Float_t factor = 1;
+ switch (loop)
+ {
+ case 0: factor = 0.5; break;
+ case 1: factor = -0.5; break;
+ case 2: factor = 1; break;
+ case 3: factor = -1; break;
+ default: return -1;
+ }
+
+ TH1* measuredClone = (TH1*) measured->Clone("measuredClone");
+
+ measuredClone->SetBinContent(m+1, measured->GetBinContent(m+1) + factor * TMath::Sqrt(measured->GetBinContent(m+1)));
+ //new TCanvas; measuredClone->Draw("COLZ");
+
+ newResult[loop]->Reset();
+ Unfold(correlation, efficiency, measuredClone, measuredClone, newResult[loop]);
+ // WARNING if we leave here, then nothing is calculated
+ //return -1;
+
+ delete measuredClone;
+ }
+
+ for (Int_t t=0; t<fgMaxParams; t++)
+ {
+ // one value
+ //matrix(i, j) = (result->GetBinContent(j+1) - newResult->GetBinContent(j+1)) / TMath::Sqrt(changedHist->GetBinContent(1, i+1));
+
+ // four values from the derivate (procedure taken from ROOT -- suggestion by Ruben)
+ // = 1/2D * [ 8 (f(D/2) - f(-D/2)) - (f(D)-f(-D)) ]/3
+
+ /*
+ // check formula
+ measured->SetBinContent(1, m+1, 100);
+ newResult[0]->SetBinContent(t+1, 5 * 0.5 + 10);
+ newResult[1]->SetBinContent(t+1, 5 * -0.5 + 10);
+ newResult[2]->SetBinContent(t+1, 5 * 1 + 10);
+ newResult[3]->SetBinContent(t+1, 5 * -1 + 10);
+ */
+
+ matrix(m, t) = 0.5 / TMath::Sqrt(measured->GetBinContent(m+1)) *
+ ( 8. * (newResult[0]->GetBinContent(t+1) - newResult[1]->GetBinContent(t+1)) -
+ (newResult[2]->GetBinContent(t+1) - newResult[3]->GetBinContent(t+1))
+ ) / 3;
+ }
+ }
+
+ for (Int_t loop=0; loop<4; loop++)
+ delete newResult[loop];
+
+ //matrix.Print();
+ //TH2* matrixHist = new TH2D(matrix); new TCanvas; matrixHist->Draw("BOX");
+
+ // ...calculate folded guess
+ TH1* convoluted = (TH1*) measured->Clone("convoluted");
+ convoluted->Reset();
+ for (Int_t m=0; m<fgMaxInput; m++)
+ {
+ Float_t value = 0;
+ for (Int_t t = 0; t<fgMaxParams; t++)
+ {
+ Float_t tmp = correlation->GetBinContent(t+1, m+1) * result->GetBinContent(t+1);
+ if (efficiency)
+ tmp *= efficiency->GetBinContent(t+1);
+ value += tmp;
+ }
+ convoluted->SetBinContent(m+1, value);
+ }
+
+ // rescale
+ convoluted->Scale(measured->Integral() / convoluted->Integral());
+
+ //new TCanvas; convoluted->Draw(); measured->Draw("SAME"); measured->SetLineColor(2);
+ //return;
+
+ // difference
+ convoluted->Add(measured, -1);
+
+ // ...and the bias
+ for (Int_t t = 0; t<fgMaxParams; t++)
+ {
+ Double_t error = 0;
+ for (Int_t m=0; m<fgMaxInput; m++)
+ error += matrix(m, t) * convoluted->GetBinContent(m+1);
+ result->SetBinError(t+1, error);
+ }
+
+ //new TCanvas; result->Draw(); gPad->SetLogy();
+
+ return 0;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff