Int_t AliUnfolding::fgSkipBinsBegin = 0;
Float_t AliUnfolding::fgMinuitStepSize = 0.1; // (usually not needed to be changed) step size in minimization
Float_t AliUnfolding::fgMinuitPrecision = 1e-6; // minuit precision
-Int_t AliUnfolding::fgMinuitMaxIterations = 1000000; // minuit maximum number of iterations\r
+Int_t AliUnfolding::fgMinuitMaxIterations = 1000000; // minuit maximum number of iterations
Double_t AliUnfolding::fgMinuitStrategy = 1.; // minuit strategy
Bool_t AliUnfolding::fgMinimumInitialValue = kFALSE; // set all initial values at least to the smallest value among the initial values
Float_t AliUnfolding::fgMinimumInitialValueFix = -1;
fgBayesianSmoothing = smoothing;
fgBayesianIterations = nIterations;
- Printf("AliUnfolding::SetBayesianParameters --> Parameters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing);\r
+ Printf("AliUnfolding::SetBayesianParameters --> Parameters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing);
}
//____________________________________________________________________
fgEntropyAPriori = new TVectorD(fgMaxParams);
if (!fgEfficiency)
fgEfficiency = new TVectorD(fgMaxParams);
- if (!fgUnfoldedAxis)
+ if (fgUnfoldedAxis)
+ {
delete fgUnfoldedAxis;
- fgUnfoldedAxis = new TAxis(*(correlation->GetXaxis()));
- if (!fgMeasuredAxis)
+ fgUnfoldedAxis = 0;
+ }
+ if (!fgUnfoldedAxis)
+ fgUnfoldedAxis = new TAxis(*(correlation->GetXaxis()));
+ if (fgMeasuredAxis)
+ {
delete fgMeasuredAxis;
- fgMeasuredAxis = new TAxis(*(correlation->GetYaxis()));
+ fgMeasuredAxis = 0;
+ }
+ if (!fgMeasuredAxis)
+ fgMeasuredAxis = new TAxis(*(correlation->GetYaxis()));
fgCorrelationMatrix->Zero();
fgCorrelationCovarianceMatrix->Zero();
if (sum <= 0)
continue;
Float_t maxValue = 0;
- Int_t maxBin = -1;
+ // 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;
+ // maxBin = j;
}
// npart sum to 1
Double_t* result = new Double_t[kMaxT];
Double_t* efficiency = new Double_t[kMaxT];
Double_t* binWidths = new Double_t[kMaxT];
- Double_t* normResponse = new Double_t[kMaxT]; \r
+ Double_t* normResponse = new Double_t[kMaxT];
Double_t** response = new Double_t*[kMaxT];
Double_t** inverseResponse = new Double_t*[kMaxT];
{
response[i] = new Double_t[kMaxM];
inverseResponse[i] = new Double_t[kMaxM];
- normResponse[i] = 0;\r
+ normResponse[i] = 0;
}
// for normalization
{
response[t][m] = correlation->GetBinContent(t+1, m+1);
inverseResponse[t][m] = 0;
- normResponse[t] += correlation->GetBinContent(t+1, m+1);\r
+ normResponse[t] += correlation->GetBinContent(t+1, m+1);
}
}
prior[t] = measuredCopy[t];
result[t] = 0;
binWidths[t] = aResult->GetXaxis()->GetBinWidth(t+1);
-\r
- for (Int_t m=0; m<kMaxM; m++) { // Normalise response matrix\r
- if (normResponse[t] != 0) \r
- response[t][m] /= normResponse[t];\r
- else\r
- Printf("AliUnfolding::UnfoldWithBayesian: Empty row,column in response matrix, for truth bin %d",t);\r
- }\r
+
+ for (Int_t m=0; m<kMaxM; m++) { // Normalise response matrix
+ if (normResponse[t] != 0)
+ response[t][m] /= normResponse[t];
+ else
+ Printf("AliUnfolding::UnfoldWithBayesian: Empty row,column in response matrix, for truth bin %d",t);
+ }
}
// pick prior distribution
{
Float_t norm = 0;
for (Int_t t = kStartBin; t<kMaxT; t++)
- norm += response[t][m] * prior[t] * efficiency[t];\r
+ norm += response[t][m] * prior[t] * efficiency[t];
// calc. chi2: (measured - response * prior) / error
if (measuredError[m] > 0)
if (norm > 0)
{
for (Int_t t = kStartBin; t<kMaxT; t++)
- inverseResponse[t][m] = response[t][m] * prior[t] * efficiency[t] / norm;\r
+ inverseResponse[t][m] = response[t][m] * prior[t] * efficiency[t] / norm;
}
else
{
delete[] result;
delete[] efficiency;
delete[] binWidths;
- delete[] normResponse;\r
+ delete[] normResponse;
for (Int_t i=0; i<kMaxT; i++)
{
middle = TMath::Power(params[i-1] / fgUnfoldedAxis->GetBinWidth(i),fgPowern);
if(middle>0) {
- right = TMath::Power(params[i] / fgUnfoldedAxis->GetBinWidth(i),fgPowern)/middle;
+ right = TMath::Power(params[i] / fgUnfoldedAxis->GetBinWidth(i+1),fgPowern)/middle;
left = TMath::Power(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1),fgPowern)/middle;
middle = 1.;
Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
- Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i-2)) / 2);
+ Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
- chi2 += (der1 - der2) * (der1 - der2)/ ( fgUnfoldedAxis->GetBinWidth(i)/2. + fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i-2)/2.)/( fgUnfoldedAxis->GetBinWidth(i)/2. + fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i-2)/2.);// / error;
+ chi2 += (der1 - der2) * (der1 - der2)/ ( fgUnfoldedAxis->GetBinWidth(i+1)/2. + fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)/2.)/( fgUnfoldedAxis->GetBinWidth(i+1)/2. + fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)/2.);// / error;
// printf("i: %d chi2 = %f\n",i,chi2);
}
if ((*fgEfficiency)(i-1) == 0 || (*fgEfficiency)(i) == 0 || (*fgEfficiency)(i-2) == 0)
continue;
-
- Double_t right = log(params[i] / (*fgEfficiency)(i) / fgUnfoldedAxis->GetBinWidth(i));
- Double_t middle = log(params[i-1] / (*fgEfficiency)(i-1) / fgUnfoldedAxis->GetBinWidth(i-1));
- Double_t left = log(params[i-2] / (*fgEfficiency)(i-2) / fgUnfoldedAxis->GetBinWidth(i-2));
+ Double_t right = log(params[i] / (*fgEfficiency)(i) / fgUnfoldedAxis->GetBinWidth(i+1));
+ Double_t middle = log(params[i-1] / (*fgEfficiency)(i-1) / fgUnfoldedAxis->GetBinWidth(i));
+ Double_t left = log(params[i-2] / (*fgEfficiency)(i-2) / fgUnfoldedAxis->GetBinWidth(i-1));
Double_t der1 = (right - middle) / ( log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i)) );
Double_t der2 = (middle - left) /( log(fgUnfoldedAxis->GetBinCenter(i)) - log(fgUnfoldedAxis->GetBinCenter(i-1)) );
-
+
double tmp = (log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i-1)))/2.;
Double_t dder = (der1-der2) / tmp;
{
params[i] = initialConditions->GetBinContent(i+1) * efficiency->GetBinContent(i+1);
- Bool_t fix = kFALSE;
if (params[i] < 0)
- {
- fix = kTRUE;
params[i] = -params[i];
- }
params[i]=TMath::Sqrt(params[i]);
-
- //cout << "params[" << i << "] " << params[i] << endl;
-
}
Double_t chi2;