#include <TCollection.h>
#include <TLegend.h>
#include <TLine.h>
+#include <TRandom.h>
+#include <TProfile.h>
+#include <TProfile2D.h>
ClassImp(AliMultiplicityCorrection)
Float_t AliMultiplicityCorrection::fRegularizationWeight = 5000;
TF1* AliMultiplicityCorrection::fNBD = 0;
+// 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!)
+
+// SPD: peak area - flat area - low stat area
+Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {4, 60, 180};
+Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {20, 140, 210};
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy)
+{
+ //
+ // sets the quality region definition to TPC or SPD
+ //
+
+ if (SPDStudy)
+ {
+ // SPD: peak area - flat area - low stat area
+ fgQualityRegionsB[0] = 4;
+ fgQualityRegionsE[0] = 20;
+
+ fgQualityRegionsB[1] = 60;
+ fgQualityRegionsE[1] = 140;
+
+ fgQualityRegionsB[2] = 180;
+ fgQualityRegionsE[2] = 210;
+ }
+ else
+ {
+ // TPC: peak area - flat area - low stat area
+ fgQualityRegionsB[0] = 4;
+ fgQualityRegionsE[0] = 12;
+
+ fgQualityRegionsB[1] = 25;
+ fgQualityRegionsE[1] = 55;
+
+ fgQualityRegionsB[2] = 88;
+ fgQualityRegionsE[2] = 108;
+ }
+}
+
//____________________________________________________________________
AliMultiplicityCorrection::AliMultiplicityCorrection() :
- TNamed(), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0)
+ TNamed(), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
{
//
// default constructor
//____________________________________________________________________
AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
- TNamed(name, title), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0)
+ TNamed(name, title), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
{
//
// named constructor
}
//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationTest(TVectorD& params)
+Double_t AliMultiplicityCorrection::RegularizationLog(TVectorD& params)
{
// homogenity term for minuit fitting
// pure function of the parameters
Double_t chi2 = 0;
- Float_t der[fgMaxParams];
+ /*Float_t der[fgMaxParams];
for (Int_t i=0; i<fgMaxParams-1; ++i)
der[i] = params[i+1] - params[i];
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<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; // 10000
+ return chi2 * 100;
}
//____________________________________________________________________
Double_t chi2 = 0;
for (Int_t i=1; i<fgMaxParams; ++i)
{
- Double_t tmp = params[i] / paramSum;
+ //Double_t tmp = params[i] / paramSum;
+ Double_t tmp = params[i];
if (tmp > 0 && (*fEntropyAPriori)[i] > 0)
- chi2 += tmp * log(tmp / (*fEntropyAPriori)[i]);
+ chi2 += tmp * TMath::Log(tmp / (*fEntropyAPriori)[i]);
}
return 10.0 + chi2;
//
// d
- TVectorD paramsVector(fgMaxParams);
+ static TVectorD paramsVector(fgMaxParams);
for (Int_t i=0; i<fgMaxParams; ++i)
paramsVector[i] = params[i] * params[i];
case kPol1: penaltyVal = RegularizationPol1(paramsVector); break;
case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break;
case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break;
- case kTest: penaltyVal = RegularizationTest(paramsVector); break;
+ case kLog: penaltyVal = RegularizationLog(paramsVector); break;
}
//if (penaltyVal > 1e10)
TVectorD copy(measGuessVector);
// (Ad - m) W
- measGuessVector *= (*fCorrelationCovarianceMatrix);
+ // this step can be optimized because currently only the diagonal elements of fCorrelationCovarianceMatrix are used
+ // normal way is like this:
+ // measGuessVector *= (*fCorrelationCovarianceMatrix);
+ // optimized way like this:
+ for (Int_t i=0; i<fgMaxParams; ++i)
+ measGuessVector[i] *= (*fCorrelationCovarianceMatrix)(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 ((*fCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
Double_t chi2FromFit = measGuessVector * copy * 1e6;
chi2 = chi2FromFit + penaltyVal;
// resets fMultiplicityESDCorrected
//
+ Bool_t display = kFALSE;
+
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
+
fMultiplicityESDCorrected[correlationID]->Reset();
+ fMultiplicityESDCorrected[correlationID]->Sumw2();
- fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY();
+ fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD"); //"fCurrentESD", -1, -1, "e");
fCurrentESD->Sumw2();
// empty under/overflow bins in x, otherwise Project3D takes them into account
if (maxBin > 0)
{
- TCanvas* canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
- canvas->Divide(2, 2);
-
- canvas->cd(1);
- fCurrentESD->DrawCopy();
- gPad->SetLogy();
+ TCanvas* canvas = 0;
- canvas->cd(2);
- fCurrentCorrelation->DrawCopy("COLZ");
+ 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", maxBin);
fCurrentESD->SetBinContent(maxBin, fCurrentESD->Integral(maxBin, fCurrentESD->GetNbinsX()));
printf("Adjusted correlation matrix!\n");
- canvas->cd(3);
- fCurrentESD->DrawCopy();
- gPad->SetLogy();
+ if (canvas)
+ {
+ canvas->cd(3);
+ fCurrentESD->DrawCopy();
+ gPad->SetLogy();
- canvas->cd(4);
- fCurrentCorrelation->DrawCopy("COLZ");
+ canvas->cd(4);
+ fCurrentCorrelation->DrawCopy("COLZ");
+ }
}
}
- //normalize ESD
- fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
+#if 0 // does not help
+ // null bins with one entry
+ Int_t nNulledBins = 0;
+ for (Int_t x=1; x<=fCurrentCorrelation->GetXaxis()->GetNbins(); ++x)
+ for (Int_t y=1; y<=fCurrentCorrelation->GetYaxis()->GetNbins(); ++y)
+ {
+ if (fCurrentCorrelation->GetBinContent(x, y) == 1)
+ {
+ fCurrentCorrelation->SetBinContent(x, y, 0);
+ fCurrentCorrelation->SetBinError(x, y, 0);
- fCurrentEfficiency = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency");
- TH2* divisor = 0;
+ ++nNulledBins;
+ }
+ }
+ Printf("Nulled %d bins", nNulledBins);
+#endif
+
+ fCurrentEfficiency = GetEfficiency(inputRange, eventType);
+ //fCurrentEfficiency->Rebin(2);
+ //fCurrentEfficiency->Scale(0.5);
+}
+
+//____________________________________________________________________
+TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventType)
+{
+ //
+ // calculates efficiency for given event type
+ //
+
+ TH1* divisor = 0;
switch (eventType)
{
- case kTrVtx : divisor = fMultiplicityVtx[inputRange]; break;
- case kMB: divisor = fMultiplicityMB[inputRange]; break;
- case kINEL: divisor = fMultiplicityINEL[inputRange]; break;
+ case kTrVtx : divisor = fMultiplicityVtx[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
+ case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
+ case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
}
- fCurrentEfficiency->Divide(divisor->ProjectionY());
- //fCurrentEfficiency->Rebin(2);
- //fCurrentEfficiency->Scale(0.5);
+ TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", -1, -1, "e");
+ eff->Divide(eff, divisor, 1, 1, "B");
+ return eff;
}
//____________________________________________________________________
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
- SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE); // TODO FAKE kTRUE
+ SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kTRUE);
+ //normalize ESD
+ fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
fCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
fCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
// 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 out why...
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<fgMaxParams; i++)
//proj->Rebin(2);
proj->Scale(1.0 / proj->Integral());
- Double_t results[fgMaxParams];
+ Double_t results[fgMaxParams+1];
for (Int_t i=0; i<fgMaxParams; ++i)
{
results[i] = inputDist->GetBinContent(i+1);
minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0);
}
+ //minuit->SetParameter(fgMaxParams, "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;
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
+ //normalize ESD
+ fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
fCorrelationMatrix = new TMatrixD(fgMaxParams, fgMaxParams);
fCorrelationCovarianceMatrix = new TMatrixD(fgMaxParams, fgMaxParams);
delete ratio;
mcHist->Scale(scalingFactor);*/
+ // find bin with <= 100 entries. this is used as limit for the chi2 calculation
+ Int_t mcBinLimit = 0;
+ for (Int_t i=20; i<mcHist->GetNbinsX(); ++i)
+ {
+ if (mcHist->GetBinContent(i) > 50)
+ {
+ mcBinLimit = i;
+ }
+ else
+ break;
+ }
+ Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit);
+
+ // scale to 1
+ mcHist->Sumw2();
mcHist->Scale(1.0 / mcHist->Integral());
// calculate residual
residual->Add(fCurrentESD, -1);
//residual->Divide(residual, fCurrentESD, 1, 1, "B");
- TH1* residualHist = new TH1F("residualHist", "residualHist", 50, -5, 5);
+ TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
// TODO fix errors
Float_t chi2 = 0;
fLastChi2Residuals = chi2;
}
- residualHist->Fit("gaus", "N");
- delete residualHist;
+ new TCanvas; residualHist->DrawCopy();
+
+ //residualHist->Fit("gaus", "N");
+ //delete residualHist;
printf("Difference (Residuals) is %f for bin 2-200\n", chi2);
legend->SetFillColor(0);
legend->Draw();
- TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
- diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
+ //TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
+ //diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
/*Float_t chi2 = 0;
Float_t chi = 0;
}
}*/
- chi2 = 0;
+ /*chi2 = 0;
Float_t chi = 0;
Int_t limit = 0;
for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
printf("limits %d %d\n", fLastChi2MCLimit, limit);
fLastChi2MCLimit = limit;
- printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);
+ printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);*/
if (!simple)
{
- canvas1->cd(3);
+ /*canvas1->cd(3);
- diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e");
+ diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(unfolded)");
//diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
diffMCUnfolded->DrawCopy("HIST");
for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
- new TCanvas;
- fluctuation->Draw();
+ //new TCanvas; fluctuation->DrawCopy();
+ delete fluctuation;*/
/*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
ratio->Draw("HIST");
Double_t ratioChi2 = 0;
+ fRatioAverage = 0;
fLastChi2MCLimit = 0;
for (Int_t i=2; i<=150; ++i)
{
value = -1e8;
ratioChi2 += value * value;
+ fRatioAverage += TMath::Abs(value);
if (ratioChi2 < 0.1)
fLastChi2MCLimit = i;
}
+ fRatioAverage /= 149;
- printf("Sum over (ratio-1)^2 (2..150) is %f\n", ratioChi2);
+ printf("Sum over (ratio-1)^2 (2..150) is %f; average of |ratio-1| is %f\n", ratioChi2, fRatioAverage);
// TODO FAKE
fLastChi2MC = ratioChi2;
+
+ // FFT of ratio
+ canvas1->cd(6);
+ const Int_t kFFT = 128;
+ Double_t fftReal[kFFT];
+ Double_t fftImag[kFFT];
+
+ for (Int_t i=0; i<kFFT; ++i)
+ {
+ fftReal[i] = ratio->GetBinContent(i+1+10);
+ // test: ;-)
+ //fftReal[i] = cos(TMath::Pi() * 5 * 2 * i / 128);
+ fftImag[i] = 0;
+ }
+
+ FFT(-1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
+
+ TH1* fftResult = (TH1*) ratio->Clone("fftResult");
+ fftResult->SetTitle("FFT;true multiplicity;coeff. (10...137)");
+ TH1* fftResult2 = (TH1*) ratio->Clone("fftResult2");
+ TH1* fftResult3 = (TH1*) ratio->Clone("fftResult3");
+ fftResult->Reset();
+ fftResult2->Reset();
+ fftResult3->Reset();
+ fftResult->GetYaxis()->UnZoom();
+ fftResult2->GetYaxis()->UnZoom();
+ fftResult3->GetYaxis()->UnZoom();
+
+ Printf("AFTER FFT");
+ for (Int_t i=0; i<kFFT; ++i)
+ {
+ Printf("%d: %f", i, fftReal[i]);
+ fftResult->SetBinContent(i+1, fftReal[i]);
+ /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5)
+ {
+ Printf("Nulled %d", i);
+ fftReal[i] = 0;
+ }*/
+ }
+
+ fftResult->SetLineColor(1);
+ fftResult->DrawCopy();
+
+
+ Printf("IMAG");
+ for (Int_t i=0; i<kFFT; ++i)
+ {
+ Printf("%d: %f", i, fftImag[i]);
+ fftResult2->SetBinContent(i+1, fftImag[i]);
+
+ fftResult3->SetBinContent(i+1, TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]));
+ }
+
+ fftResult2->SetLineColor(2);
+ fftResult2->DrawCopy("SAME");
+
+ fftResult3->SetLineColor(4);
+ fftResult3->DrawCopy("SAME");
+
+ for (Int_t i=1; i<kFFT - 1; ++i)
+ {
+ if (TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]) > 3)
+ {
+ fftReal[i-1] = 0;
+ fftReal[i] = 0;
+ fftReal[i+1] = 0;
+ fftImag[i-1] = 0;
+ fftImag[i] = 0;
+ fftImag[i+1] = 0;
+ //fftReal[i] = (fftReal[i-1] + fftReal[i+1]) / 2;
+ //fftImag[i] = (fftImag[i-1] + fftImag[i+1]) / 2;
+ Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]);
+ }
+ }
+
+ FFT(1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
+
+ TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4");
+ fftResult4->Reset();
+
+ Printf("AFTER BACK-TRAFO");
+ for (Int_t i=0; i<kFFT; ++i)
+ {
+ Printf("%d: %f", i, fftReal[i]);
+ fftResult4->SetBinContent(i+1+10, fftReal[i]);
+ }
+
+ canvas1->cd(5);
+ fftResult4->SetLineColor(4);
+ fftResult4->DrawCopy("SAME");
+
+ // plot (MC - Unfolded) / error (MC)
+ canvas1->cd(3);
+
+ TH1* diffMCUnfolded2 = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded2"));
+ diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1);
+
+ Int_t ndfQual[kQualityRegions];
+ for (Int_t region=0; region<kQualityRegions; ++region)
+ {
+ fQuality[region] = 0;
+ ndfQual[region] = 0;
+ }
+
+
+ Double_t newChi2 = 0;
+ Double_t newChi2_150 = 0;
+ Int_t ndf = 0;
+ for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgMaxParams+1); ++i)
+ {
+ Double_t value = 0;
+ if (proj->GetBinError(i) > 0)
+ {
+ value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
+ newChi2 += value * value;
+ if (i > 1 && i <= mcBinLimit)
+ newChi2_150 += value * value;
+ ++ndf;
+
+ for (Int_t region=0; region<kQualityRegions; ++region)
+ if (diffMCUnfolded2->GetXaxis()->GetBinCenter(i) >= fgQualityRegionsB[region] - 0.1 && diffMCUnfolded2->GetXaxis()->GetBinCenter(i) <= fgQualityRegionsE[region] + 0.1) // 0.1 to avoid e.g. 3.9999 < 4 problem
+ {
+ fQuality[region] += TMath::Abs(value);
+ ++ndfQual[region];
+ }
+ }
+
+ diffMCUnfolded2->SetBinContent(i, value);
+ }
+
+ // normalize region to the number of entries
+ for (Int_t region=0; region<kQualityRegions; ++region)
+ {
+ if (ndfQual[region] > 0)
+ fQuality[region] /= ndfQual[region];
+ Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]);
+ }
+
+ if (mcBinLimit > 1)
+ {
+ // TODO another FAKE
+ fLastChi2MC = newChi2_150 / (mcBinLimit - 1);
+ Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2_150, mcBinLimit - 1, fLastChi2MC);
+ }
+ else
+ fLastChi2MC = -1;
+
+ Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, newChi2 / ndf);
+
+ diffMCUnfolded2->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(MC)");
+ //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
+ diffMCUnfolded2->GetXaxis()->SetRangeUser(0, 200);
+ diffMCUnfolded2->DrawCopy("HIST");
}
canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
}
//____________________________________________________________________
-void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals)
+void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y)
+{
+/*-------------------------------------------------------------------------
+ This computes an in-place complex-to-complex FFT
+ x and y are the real and imaginary arrays of 2^m points.
+ dir = 1 gives forward transform
+ dir = -1 gives reverse transform
+
+ Formula: forward
+ N-1
+ ---
+ 1 \ - j k 2 pi n / N
+ X(n) = --- > x(k) e = forward transform
+ N / n=0..N-1
+ ---
+ k=0
+
+ Formula: reverse
+ N-1
+ ---
+ \ j k 2 pi n / N
+ X(n) = > x(k) e = forward transform
+ / n=0..N-1
+ ---
+ k=0
+*/
+
+ Long_t nn, i, i1, j, k, i2, l, l1, l2;
+ Double_t c1, c2, tx, ty, t1, t2, u1, u2, z;
+
+ /* Calculate the number of points */
+ nn = 1;
+ for (i = 0; i < m; i++)
+ nn *= 2;
+
+ /* Do the bit reversal */
+ i2 = nn >> 1;
+ j = 0;
+ for (i= 0; i < nn - 1; i++) {
+ if (i < j) {
+ tx = x[i];
+ ty = y[i];
+ x[i] = x[j];
+ y[i] = y[j];
+ x[j] = tx;
+ y[j] = ty;
+ }
+ k = i2;
+ while (k <= j) {
+ j -= k;
+ k >>= 1;
+ }
+ j += k;
+ }
+
+ /* Compute the FFT */
+ c1 = -1.0;
+ c2 = 0.0;
+ l2 = 1;
+ for (l = 0; l < m; l++) {
+ l1 = l2;
+ l2 <<= 1;
+ u1 = 1.0;
+ u2 = 0.0;
+ for (j = 0;j < l1; j++) {
+ for (i = j; i < nn; i += l2) {
+ i1 = i + l1;
+ t1 = u1 * x[i1] - u2 * y[i1];
+ t2 = u1 * y[i1] + u2 * x[i1];
+ x[i1] = x[i] - t1;
+ y[i1] = y[i] - t2;
+ x[i] += t1;
+ y[i] += t2;
+ }
+ z = u1 * c1 - u2 * c2;
+ u2 = u1 * c2 + u2 * c1;
+ u1 = z;
+ }
+ c2 = TMath::Sqrt((1.0 - c1) / 2.0);
+ if (dir == 1)
+ c2 = -c2;
+ c1 = TMath::Sqrt((1.0 + c1) / 2.0);
+ }
+
+ /* Scaling for forward transform */
+ if (dir == 1) {
+ for (i=0;i<nn;i++) {
+ x[i] /= (Double_t)nn;
+ y[i] /= (Double_t)nn;
+ }
+ }
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals, Float_t* ratioAverage)
{
// Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
// These values are computed during DrawComparison, thus this function picks up the
*mcLimit = fLastChi2MCLimit;
if (residuals)
*residuals = fLastChi2Residuals;
+ if (ratioAverage)
+ *ratioAverage = fRatioAverage;
}
-
//____________________________________________________________________
TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
{
}
//____________________________________________________________________
-void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* inputDist)
+TH1* AliMultiplicityCorrection::BayesianStatisticsEffect(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, Float_t regPar, Int_t nIterations, TH1* compareTo)
{
//
- // correct spectrum using bayesian method
+ // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
+ // the function unfolds the spectrum using the default response matrix and several modified ones
+ // the modified ones are created by randomizing each cell using poisson statistics with the mean = bin value
+ // these unfolded results are compared to the first result gained with the default response OR to the histogram given
+ // in <compareTo> (optional)
//
+ // returns the error assigned to the measurement
+ //
+
+ // initialize seed with current time
+ gRandom->SetSeed(0);
+
+ const Int_t kErrorIterations = 150;
+
+ TH1* maxError = 0;
+ TH1* firstResult = 0;
+
+ TH1** results = new TH1*[kErrorIterations];
+
+ for (Int_t n=0; n<kErrorIterations; ++n)
+ {
+ Printf("Iteration %d of %d...", n, kErrorIterations);
+
+ SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
+
+ TH1* measured = (TH1*) fCurrentESD->Clone("measured");
+
+ if (n > 0)
+ {
+ if (randomizeResponse)
+ {
+ // randomize response matrix
+ for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
+ for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
+ fCurrentCorrelation->SetBinContent(i, j, gRandom->Poisson(fCurrentCorrelation->GetBinContent(i, j)));
+ }
+
+ if (randomizeMeasured)
+ {
+ // randomize measured spectrum
+ for (Int_t x=1; x<=measured->GetNbinsX(); x++) // mult. axis
+ {
+ Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
+ measured->SetBinContent(x, randomValue);
+ measured->SetBinError(x, TMath::Sqrt(randomValue));
+ }
+ }
+ }
+
+ for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
+ {
+ // with this it is normalized to 1
+ Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
+
+ // with this normalized to the given efficiency
+ if (fCurrentEfficiency->GetBinContent(i) > 0)
+ sum /= fCurrentEfficiency->GetBinContent(i);
+ else
+ sum = 0;
+
+ for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
+ {
+ if (sum > 0)
+ {
+ fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
+ fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
+ }
+ else
+ {
+ fCurrentCorrelation->SetBinContent(i, j, 0);
+ fCurrentCorrelation->SetBinError(i, j, 0);
+ }
+ }
+ }
+
+ TH1* result = 0;
+ if (n == 0 && compareTo)
+ {
+ // in this case we just store the histogram we want to compare to
+ result = (TH1*) compareTo->Clone("compareTo");
+ result->Sumw2();
+ }
+ else
+ result = UnfoldWithBayesian(measured, regPar, nIterations, 0);
+
+ if (!result)
+ return 0;
+
+ // normalize
+ result->Scale(1.0 / result->Integral());
+
+ if (n == 0)
+ {
+ firstResult = (TH1*) result->Clone("firstResult");
+
+ maxError = (TH1*) result->Clone("maxError");
+ maxError->Reset();
+ }
+ else
+ {
+ // calculate ratio
+ TH1* ratio = (TH1*) firstResult->Clone("ratio");
+ ratio->Divide(result);
+
+ // find max. deviation
+ for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
+ maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
+
+ delete ratio;
+ }
+
+ results[n] = result;
+ }
+
+ // find covariance matrix
+ // results[n] is X_x
+ // cov. matrix is M_xy = E ( (X_x - E(X_x)) * (X_y - E(X_y))), with E() = expectation value
+
+ Int_t nBins = results[0]->GetNbinsX();
+ Float_t lowEdge = results[0]->GetXaxis()->GetBinLowEdge(1);
+ Float_t upEdge = results[0]->GetXaxis()->GetBinUpEdge(nBins);
+
+ // find average, E(X)
+ TProfile* average = new TProfile("average", "average", nBins, lowEdge, upEdge);
+ for (Int_t n=1; n<kErrorIterations; ++n)
+ for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
+ average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x));
+ //new TCanvas; average->DrawClone();
+
+ // find cov. matrix
+ TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge);
+
+ for (Int_t n=1; n<kErrorIterations; ++n)
+ for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
+ for (Int_t y=1; y<=results[n]->GetNbinsX(); y++)
+ {
+ // (X_x - E(X_x)) * (X_y - E(X_y)
+ Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y));
+ covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov);
+ }
+ new TCanvas; covMatrix->DrawClone("COLZ");
+
+ // fill 2D histogram that contains deviation from first
+ TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
+ for (Int_t n=1; n<kErrorIterations; ++n)
+ for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
+ deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x));
+ new TCanvas; deviations->DrawCopy("COLZ");
+
+ TH1* standardDeviation = (TH1*) maxError->Clone("standardDeviation");
+ standardDeviation->Reset();
+ // get standard deviation "by hand"
+ for (Int_t x=1; x<=nBins; x++)
+ {
+ if (results[0]->GetBinContent(x) > 0)
+ {
+ TH1* proj = deviations->ProjectionY("projRMS", x, x, "e");
+ Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact
+ standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
+ }
+ }
+
+ // compare maxError to RMS of profile (variable name: average)
+ // first: calculate rms in percent of value
+ TH1* rmsError = (TH1*) maxError->Clone("rmsError");
+ rmsError->Reset();
+
+ // enable error to be standard deviation (see http://root.cern.ch/root/html/TProfile.html#TProfile:SetErrorOption)
+ average->SetErrorOption("s");
+ for (Int_t x=1; x<=rmsError->GetNbinsX(); x++)
+ if (average->GetBinContent(x) > 0)
+ rmsError->SetBinContent(x, average->GetBinError(x) / average->GetBinContent(x));
+
+ // find maxError deviation from average (not from "first result"/mc as above)
+ TH1* maxError2 = (TH1*) maxError->Clone("maxError2");
+ maxError2->Reset();
+ for (Int_t n=1; n<kErrorIterations; ++n)
+ for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
+ if (average->GetBinContent(x) > 0)
+ maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x))));
+
+ new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME");
+
+ // plot difference between average and MC/first result
+ TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio");
+ averageFirstRatio->Reset();
+ averageFirstRatio->Divide(results[0], average);
+
+ new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
+ new TCanvas; averageFirstRatio->DrawCopy();
+
+ // clean up
+ for (Int_t n=0; n<kErrorIterations; ++n)
+ delete results[n];
+ delete[] results;
+
+ // fill into result histogram
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
- SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
+ for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
+ fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
- // smooth efficiency
- //TH1* tmp = (TH1*) fCurrentEfficiency->Clone("eff_clone");
- //for (Int_t i=2; i<fCurrentEfficiency->GetNbinsX(); ++i)
- // fCurrentEfficiency->SetBinContent(i, (tmp->GetBinContent(i-1) + tmp->GetBinContent(i) + tmp->GetBinContent(i+1)) / 3);
+ for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
+ fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
- // set efficiency to 1 above 150
- // FAKE TEST
- //for (Int_t i=150; i<=fCurrentEfficiency->GetNbinsX(); ++i)
- // fCurrentEfficiency->SetBinContent(i, 1);
+ return standardDeviation;
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* inputDist, Bool_t determineError)
+{
+ //
+ // correct spectrum using bayesian method
+ //
+
+ // initialize seed with current time
+ gRandom->SetSeed(0);
+
+ SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
// normalize correction for given nPart
for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
}
}
- // normalize nTrack
- /*for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
+ Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
+
+ TH1* result = UnfoldWithBayesian(fCurrentESD, regPar, nIterations, inputDist);
+ if (!result)
+ return;
+
+ for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
+ fMultiplicityESDCorrected[correlationID]->SetBinContent(i, result->GetBinContent(i));
+
+ if (!determineError)
{
- // with this it is normalized to 1
- Double_t sum = fCurrentCorrelation->Integral(1, fCurrentCorrelation->GetNbinsX(), j, j);
+ Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
+ return;
+ }
- for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
+ // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
+ // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
+ // spectrum calculated. This is performed N times and the maximum difference is taken as the systematic
+ // error of the unfolding method itself.
+
+ TH1* maxError = (TH1*) result->Clone("maxError");
+ maxError->Reset();
+
+ TH1* normalizedResult = (TH1*) result->Clone("normalizedResult");
+ normalizedResult->Scale(1.0 / normalizedResult->Integral());
+
+ const Int_t kErrorIterations = 20;
+
+ printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
+
+ TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
+ for (Int_t n=0; n<kErrorIterations; ++n)
+ {
+ // randomize the content of clone following a poisson with the mean = the value of that bin
+ for (Int_t x=1; x<=randomized->GetNbinsX(); x++) // mult. axis
{
- if (sum > 0)
- {
- fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
- fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
- }
- else
- {
- fCurrentCorrelation->SetBinContent(i, j, 0);
- fCurrentCorrelation->SetBinError(i, j, 0);
- }
+ Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
+ //printf("%d --> %d\n", fCurrentESD->GetBinContent(x), randomValue);
+ randomized->SetBinContent(x, randomValue);
+ randomized->SetBinError(x, TMath::Sqrt(randomValue));
}
- }*/
- // smooth input spectrum
- /*
- TH1* esdClone = (TH1*) fCurrentESD->Clone("esdClone");
+ TH1* result2 = UnfoldWithBayesian(randomized, regPar, nIterations, inputDist);
+ if (!result2)
+ return;
- for (Int_t i=2; i<fCurrentESD->GetNbinsX(); ++i)
- if (esdClone->GetBinContent(i) < 1e-3)
- fCurrentESD->SetBinContent(i, (esdClone->GetBinContent(i-1) + esdClone->GetBinContent(i) + esdClone->GetBinContent(i+1)) / 3);
+ result2->Scale(1.0 / result2->Integral());
- delete esdClone;
+ // calculate ratio
+ TH1* ratio = (TH1*) result2->Clone("ratio");
+ ratio->Divide(normalizedResult);
- // rescale to 1
- fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
- */
+ //new TCanvas; ratio->DrawCopy("HIST");
- /*new TCanvas;
- fCurrentEfficiency->Draw();
+ // find max. deviation
+ for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
+ maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
- new TCanvas;
- fCurrentCorrelation->DrawCopy("COLZ");
+ delete ratio;
+ delete result2;
+ }
- new TCanvas;
- ((TH2*) fCurrentCorrelation)->ProjectionX()->DrawCopy();
+ for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
+ fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
- new TCanvas;
- ((TH2*) fCurrentCorrelation)->ProjectionY()->DrawCopy();*/
+ delete maxError;
+ delete normalizedResult;
+}
- // pick prior distribution
- TH1* hPrior = 0;
- if (inputDist)
+//____________________________________________________________________
+TH1* AliMultiplicityCorrection::UnfoldWithBayesian(TH1* measured, Float_t regPar, Int_t nIterations, TH1* inputDist)
+{
+ //
+ // unfolds a spectrum
+ //
+
+ if (measured->Integral() <= 0)
{
- printf("Using different starting conditions...\n");
- hPrior = (TH1*)inputDist->Clone("prior");
+ Printf("AliMultiplicityCorrection::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
+ return 0;
}
- else
- hPrior = (TH1*)fCurrentESD->Clone("prior");
- Float_t norm = hPrior->Integral();
- for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
- hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
- // define temp hist
- TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
- hTemp->Reset();
+ const Int_t kStartBin = 0;
- // just a shortcut
- TH2F* hResponse = (TH2F*) fCurrentCorrelation;
+ const Int_t kMaxM = fgMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
+ const Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
- // histogram to store the inverse response calculated from bayes theorem (from prior and response) IR
- TH2F* hInverseResponseBayes = (TH2F*)hResponse->Clone("pji");
- hInverseResponseBayes->Reset();
+ // store information in arrays, to increase processing speed (~ factor 5)
+ Double_t measuredCopy[kMaxM];
+ Double_t prior[kMaxT];
+ Double_t result[kMaxT];
+ Double_t efficiency[kMaxT];
- TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5);
+ Double_t response[kMaxT][kMaxM];
+ Double_t inverseResponse[kMaxT][kMaxM];
- const Int_t kStartBin = 1;
+ // for normalization
+ Float_t measuredIntegral = measured->Integral();
+ for (Int_t m=0; m<kMaxM; m++)
+ {
+ measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
+
+ for (Int_t t=0; t<kMaxT; t++)
+ {
+ response[t][m] = fCurrentCorrelation->GetBinContent(t+1, m+1);
+ inverseResponse[t][m] = 0;
+ }
+ }
+
+ for (Int_t t=0; t<kMaxT; t++)
+ {
+ efficiency[t] = fCurrentEfficiency->GetBinContent(t+1);
+ prior[t] = measuredCopy[t];
+ result[t] = 0;
+ }
+
+ // pick prior distribution
+ if (inputDist)
+ {
+ printf("Using different starting conditions...\n");
+ // for normalization
+ Float_t inputDistIntegral = inputDist->Integral();
+ for (Int_t i=0; i<kMaxT; i++)
+ prior[i] = inputDist->GetBinContent(i+1) / inputDistIntegral;
+ }
+
+ //TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5);
// unfold...
for (Int_t i=0; i<nIterations; i++)
// calculate IR from Beyes theorem
// IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
- for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
+ for (Int_t m=0; m<kMaxM; m++)
{
Float_t norm = 0;
- for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++)
- norm += hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t);
- for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++)
+ for (Int_t t = kStartBin; t<kMaxT; t++)
+ norm += response[t][m] * prior[t];
+
+ if (norm > 0)
{
- if (norm==0)
- hInverseResponseBayes->SetBinContent(t,m,0);
- else
- hInverseResponseBayes->SetBinContent(t,m, hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t)/norm);
+ 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;
}
}
- for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++)
+ for (Int_t t = kStartBin; t<kMaxT; t++)
{
Float_t value = 0;
- for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
- value += fCurrentESD->GetBinContent(m) * hInverseResponseBayes->GetBinContent(t,m);
+ for (Int_t m=0; m<kMaxM; m++)
+ value += inverseResponse[t][m] * measuredCopy[m];
- if (fCurrentEfficiency->GetBinContent(t) > 0)
- hTemp->SetBinContent(t, value / fCurrentEfficiency->GetBinContent(t));
+ if (efficiency[t] > 0)
+ result[t] = value / efficiency[t];
else
- hTemp->SetBinContent(t, 0);
+ result[t] = 0;
}
- // this is the last guess, fill before (!) smoothing
- for (Int_t t=kStartBin; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
- {
- //as bin error put the difference to the last iteration
- //fMultiplicityESDCorrected[correlationID]->SetBinError(t, hTemp->GetBinContent(t) - fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
- fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
- fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t));
-
- //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
- }
-
- /*new TCanvas;
- fMultiplicityESDCorrected[correlationID]->DrawCopy();
- gPad->SetLogy();*/
-
// regularization (simple smoothing)
- TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
-
- for (Int_t t=kStartBin+2; t<hTrueSmoothed->GetNbinsX(); t++)
+ for (Int_t t=kStartBin; t<kMaxT; t++)
{
- Float_t average = (hTemp->GetBinContent(t-1)
- + hTemp->GetBinContent(t)
- + hTemp->GetBinContent(t+1)
- ) / 3.;
+ Float_t newValue = 0;
+ // 0 bin excluded from smoothing
+ if (t > kStartBin+1 && t<kMaxT-1)
+ {
+ Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
- // weight the average with the regularization parameter
- hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
+ // weight the average with the regularization parameter
+ newValue = (1 - regPar) * result[t] + regPar * average;
+ }
+ else
+ newValue = result[t];
+
+ prior[t] = newValue;
}
// calculate chi2 (change from last iteration)
- Double_t chi2 = 0;
-
- // use smoothed true (normalized) as new prior
- Float_t norm = 1;
- //for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
- // norm = norm + hTrueSmoothed->GetBinContent(t);
-
- for (Int_t t=kStartBin; t<=hTrueSmoothed->GetNbinsX(); t++)
+ //Double_t chi2 = 0;
+ /*for (Int_t t=kStartBin; t<=hTrueSmoothed->GetNbinsX(); t++)
{
- Float_t newValue = hTrueSmoothed->GetBinContent(t)/norm;
- Float_t diff = hPrior->GetBinContent(t) - newValue;
- chi2 += (Double_t) diff * diff;
+ Float_t newValue = hTrueSmoothed->GetBinContent(t);
+ //Float_t diff = hPrior->GetBinContent(t) - newValue;
+ //chi2 += (Double_t) diff * diff;
hPrior->SetBinContent(t, newValue);
}
-
//printf("Chi2 of %d iteration = %.10f\n", i, chi2);
-
- convergence->Fill(i+1, chi2);
+ //convergence->Fill(i+1, chi2);
//if (i % 5 == 0)
// DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
- delete hTrueSmoothed;
+ delete hTrueSmoothed;*/
} // end of iterations
//new TCanvas;
//convergence->DrawCopy();
//gPad->SetLogy();
- delete convergence;
+ //delete convergence;
+
+ // TODO this does not work when the number of bins is differnt
+ TH1* resultHist = (TH1*) measured->Clone("resultHist");
+ resultHist->Reset();
+ for (Int_t t=0; t<kMaxT; t++)
+ resultHist->SetBinContent(t+1, result[t]);
- return;
+ return resultHist;
// ********
// Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
}
//____________________________________________________________________
-Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse,
- TH1* fCurrentEfficiency, Int_t k, Int_t i, Int_t r, Int_t u)
+Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u)
{
//
// helper function for the covariance matrix of the bayesian method
Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
+ //normalize ESD
+ fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
// TODO should be taken from correlation map
//TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
+ //normalize ESD
+ fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
NormalizeToBinWidth((TH2*) fCurrentCorrelation);
for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
{
- mc->SetBinContent(mc->GetNbinsX() / 2, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
- mc->SetBinError(mc->GetNbinsX() / 2, i, 0);
+ mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
+ mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, 0);
}
//new TCanvas;
class AliMultiplicityCorrection : public TNamed {
public:
enum EventType { kTrVtx = 0, kMB, kINEL };
- enum RegularizationType { kNone = 0, kPol0, kPol1, kEntropy, kCurvature, kTest };
+ enum RegularizationType { kNone = 0, kPol0, kPol1, kLog, kEntropy, kCurvature };
+ enum { kESDHists = 4, kMCHists = 5, kCorrHists = 8, kQualityRegions = 3 };
AliMultiplicityCorrection();
AliMultiplicityCorrection(const Char_t* name, const Char_t* title);
void FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20);
- Bool_t LoadHistograms(const Char_t* dir);
+ Bool_t LoadHistograms(const Char_t* dir = 0);
void SaveHistograms();
void DrawHistograms();
void DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple = kFALSE);
void ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace);
- void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 0.1, Int_t nIterations = 15, TH1* inputDist = 0);
+ void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 1, Int_t nIterations = 100, TH1* inputDist = 0, Bool_t determineError = kTRUE);
+ TH1* BayesianStatisticsEffect(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, Float_t regPar = 1, Int_t nIterations = 100, TH1* compareTo = 0);
void ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace);
TH3F* GetCorrelation(Int_t i) { return fCorrelation[i]; }
TH1F* GetMultiplicityESDCorrected(Int_t i) { return fMultiplicityESDCorrected[i]; }
- void SetMultiplicityESD(Int_t i, TH2F* hist) { fMultiplicityESD[i] = hist; }
- void SetMultiplicityVtx(Int_t i, TH2F* hist) { fMultiplicityVtx[i] = hist; }
- void SetMultiplicityMB(Int_t i, TH2F* hist) { fMultiplicityMB[i] = hist; }
+ void SetMultiplicityESD(Int_t i, TH2F* hist) { fMultiplicityESD[i] = hist; }
+ void SetMultiplicityVtx(Int_t i, TH2F* hist) { fMultiplicityVtx[i] = hist; }
+ void SetMultiplicityMB(Int_t i, TH2F* hist) { fMultiplicityMB[i] = hist; }
void SetMultiplicityINEL(Int_t i, TH2F* hist) { fMultiplicityINEL[i] = hist; }
void SetCorrelation(Int_t i, TH3F* hist) { fCorrelation[i] = hist; }
void SetMultiplicityESDCorrected(Int_t i, TH1F* hist) { fMultiplicityESDCorrected[i] = hist; }
static void NormalizeToBinWidth(TH1* hist);
static void NormalizeToBinWidth(TH2* hist);
- void GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals);
+ void GetComparisonResults(Float_t* mc = 0, Int_t* mcLimit = 0, Float_t* residuals = 0, Float_t* ratioAverage = 0);
- protected:
- enum { kESDHists = 4, kMCHists = 5, kCorrHists = 8 };
+ TH1* GetEfficiency(Int_t inputRange, EventType eventType);
+
+ static void SetQualityRegions(Bool_t SPDStudy);
+ Float_t GetQuality(Int_t region) { return fQuality[region]; }
+
+ void FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y);
+ protected:
static const Int_t fgMaxParams; // bins in unfolded histogram = number of fit params
static const Int_t fgMaxInput; // bins in measured histogram
static Double_t RegularizationPol1(TVectorD& params);
static Double_t RegularizationTotalCurvature(TVectorD& params);
static Double_t RegularizationEntropy(TVectorD& params);
- static Double_t RegularizationTest(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 inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin);
- Float_t BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, TH1* fCurrentEfficiency, Int_t k, Int_t i, Int_t r, Int_t u);
+ Float_t BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u);
+ TH1* UnfoldWithBayesian(TH1* measured, Float_t regPar, Int_t nIterations, TH1* inputDist);
static TH1* fCurrentESD; // static variable to be accessed by MINUIT
static TH1* fCurrentCorrelation; // static variable to be accessed by MINUIT
Float_t fLastChi2MC; // last Chi2 between MC and unfolded ESD (calculated in DrawComparison)
Int_t fLastChi2MCLimit; // bin where the last chi2 breached a certain threshold, used to evaluate the multiplicity reach (calc. in DrawComparison)
Float_t fLastChi2Residuals; // last Chi2 of the ESD and the folded unfolded ESD (calculated in DrawComparison)
+ Float_t fRatioAverage; // last average of |ratio-1| where ratio = unfolded / mc (bin 2..150)
+
+ static Int_t fgQualityRegionsB[kQualityRegions]; // begin, given in multiplicity units
+ static Int_t fgQualityRegionsE[kQualityRegions]; // end
+ Float_t fQuality[kQualityRegions]; // stores the quality of the last comparison (calculated in DrawComparison). Contains 3 values that are averages of (MC - unfolded) / e(MC) in 3 regions, these are defined in fQualityRegionB,E
private:
AliMultiplicityCorrection(const AliMultiplicityCorrection&);
--- /dev/null
+//
+// plots for the note about multplicity measurements
+//
+
+const char* correctionFile = "multiplicityMC_2M.root";
+const char* measuredFile = "multiplicityMC_1M_3.root";
+Int_t etaRange = 3;
+
+const char* correctionFileTPC = "multiplicityMC_TPC_1.4M.root";
+const char* measuredFileTPC = "multiplicityMC_TPC_0.6M.root";
+Int_t etaRangeTPC = 1;
+
+void SetTPC()
+{
+ correctionFile = correctionFileTPC;
+ measuredFile = measuredFileTPC;
+ etaRange = etaRangeTPC;
+}
+
+void responseMatrixPlot()
+{
+ gSystem->Load("libPWG0base");
+
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+ TFile::Open(correctionFile);
+ mult->LoadHistograms("Multiplicity");
+
+ TH1* hist = mult->GetCorrelation(etaRange)->Project3D("zy");
+ hist->SetStats(kFALSE);
+
+ hist->SetTitle(";true multiplicity;measured multiplicity;Entries");
+ hist->GetXaxis()->SetRangeUser(0, 200);
+ hist->GetYaxis()->SetRangeUser(0, 200);
+
+ TCanvas* canvas = new TCanvas("c1", "c1", 800, 600);
+ canvas->SetRightMargin(0.15);
+ canvas->SetTopMargin(0.05);
+
+ gPad->SetLogz();
+ hist->Draw("COLZ");
+
+ canvas->SaveAs("responsematrix.eps");
+}
+
+TCanvas* DrawResultRatio(TH1* mcHist, TH1* result, TString epsName)
+{
+ // normalize unfolded result to mc hist
+ result->Scale(1.0 / result->Integral(2, 200));
+ result->Scale(mcHist->Integral(2, 200));
+
+ TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
+ canvas->Range(0, 0, 1, 1);
+
+ TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
+ pad1->Draw();
+
+ TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
+ pad2->Draw();
+
+ pad1->SetRightMargin(0.05);
+ pad2->SetRightMargin(0.05);
+
+ // no border between them
+ pad1->SetBottomMargin(0);
+ pad2->SetTopMargin(0);
+
+ pad1->cd();
+
+ mcHist->GetXaxis()->SetLabelSize(0.06);
+ mcHist->GetYaxis()->SetLabelSize(0.06);
+ mcHist->GetXaxis()->SetTitleSize(0.06);
+ mcHist->GetYaxis()->SetTitleSize(0.06);
+ mcHist->GetYaxis()->SetTitleOffset(0.6);
+
+ mcHist->GetXaxis()->SetRangeUser(0, 200);
+
+ mcHist->SetTitle(";true multiplicity;Entries");
+ mcHist->SetStats(kFALSE);
+
+ mcHist->DrawCopy("HIST E");
+ gPad->SetLogy();
+
+ result->SetLineColor(2);
+ result->DrawCopy("SAME HISTE");
+
+ TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
+ legend->AddEntry(mcHist, "true distribution");
+ legend->AddEntry(result, "unfolded distribution");
+ legend->SetFillColor(0);
+ legend->Draw();
+
+ pad2->cd();
+ pad2->SetBottomMargin(0.15);
+
+ // calculate ratio
+ mcHist->Sumw2();
+ TH1* ratio = (TH1*) mcHist->Clone("ratio");
+ result->Sumw2();
+ ratio->Divide(ratio, result, 1, 1, "");
+ ratio->GetYaxis()->SetTitle("Ratio (true / unfolded)");
+ ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
+
+ ratio->DrawCopy();
+
+ // get average of ratio
+ Float_t sum = 0;
+ for (Int_t i=2; i<=150; ++i)
+ {
+ sum += TMath::Abs(ratio->GetBinContent(i) - 1);
+ }
+ sum /= 149;
+
+ printf("Average (2..150) of |ratio - 1| is %f\n", sum);
+
+ TLine* line = new TLine(0, 1, 200, 1);
+ line->SetLineWidth(2);
+ line->Draw();
+
+ line = new TLine(0, 1.1, 200, 1.1);
+ line->SetLineWidth(2);
+ line->SetLineStyle(2);
+ line->Draw();
+ line = new TLine(0, 0.9, 200, 0.9);
+ line->SetLineWidth(2);
+ line->SetLineStyle(2);
+ line->Draw();
+
+ canvas->Modified();
+
+ canvas->SaveAs(epsName);
+
+ return canvas;
+}
+
+TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsName)
+{
+ // draws the 3 plots in the upper plot
+ // draws the ratio between result1 and result2 in the lower plot
+
+ // normalize unfolded result to mc hist
+ result1->Scale(1.0 / result1->Integral(2, 200));
+ result1->Scale(mcHist->Integral(2, 200));
+ result2->Scale(1.0 / result2->Integral(2, 200));
+ result2->Scale(mcHist->Integral(2, 200));
+
+ TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
+ canvas->Range(0, 0, 1, 1);
+
+ TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
+ pad1->Draw();
+
+ TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
+ pad2->Draw();
+
+ pad1->SetRightMargin(0.05);
+ pad2->SetRightMargin(0.05);
+
+ // no border between them
+ pad1->SetBottomMargin(0);
+ pad2->SetTopMargin(0);
+
+ pad1->cd();
+
+ mcHist->GetXaxis()->SetLabelSize(0.06);
+ mcHist->GetYaxis()->SetLabelSize(0.06);
+ mcHist->GetXaxis()->SetTitleSize(0.06);
+ mcHist->GetYaxis()->SetTitleSize(0.06);
+ mcHist->GetYaxis()->SetTitleOffset(0.6);
+
+ mcHist->GetXaxis()->SetRangeUser(0, 150);
+
+ mcHist->SetTitle(";true multiplicity;Entries");
+ mcHist->SetStats(kFALSE);
+
+ mcHist->DrawCopy("HIST E");
+ gPad->SetLogy();
+
+ result1->SetLineColor(2);
+ result1->DrawCopy("SAME HISTE");
+
+ result2->SetLineColor(4);
+ result2->DrawCopy("SAME HISTE");
+
+ TLegend* legend = new TLegend(0.55, 0.6, 0.95, 0.9);
+ legend->AddEntry(mcHist, "true distribution");
+ legend->AddEntry(result1, "unfolded distribution (syst)");
+ legend->AddEntry(result2, "unfolded distribution (normal)");
+ legend->SetFillColor(0);
+ legend->Draw();
+
+ pad2->cd();
+ pad2->SetBottomMargin(0.15);
+
+ result1->GetXaxis()->SetLabelSize(0.06);
+ result1->GetYaxis()->SetLabelSize(0.06);
+ result1->GetXaxis()->SetTitleSize(0.06);
+ result1->GetYaxis()->SetTitleSize(0.06);
+ result1->GetYaxis()->SetTitleOffset(0.6);
+
+ result1->GetXaxis()->SetRangeUser(0, 150);
+
+ result1->SetTitle(";true multiplicity;Entries");
+ result1->SetStats(kFALSE);
+
+ // calculate ratio
+ result1->Sumw2();
+ TH1* ratio = (TH1*) result1->Clone("ratio");
+ result2->Sumw2();
+ ratio->Divide(ratio, result2, 1, 1, "");
+ ratio->GetYaxis()->SetTitle("Ratio (syst / normal)");
+ ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
+
+ ratio->DrawCopy();
+
+ // get average of ratio
+ Float_t sum = 0;
+ for (Int_t i=2; i<=150; ++i)
+ {
+ sum += TMath::Abs(ratio->GetBinContent(i) - 1);
+ }
+ sum /= 149;
+
+ printf("Average (2..150) of |ratio - 1| is %f\n", sum);
+
+ TLine* line = new TLine(0, 1, 150, 1);
+ line->SetLineWidth(2);
+ line->Draw();
+
+ line = new TLine(0, 1.1, 150, 1.1);
+ line->SetLineWidth(2);
+ line->SetLineStyle(2);
+ line->Draw();
+ line = new TLine(0, 0.9, 150, 0.9);
+ line->SetLineWidth(2);
+ line->SetLineStyle(2);
+ line->Draw();
+
+ canvas->Modified();
+
+ canvas->SaveAs(epsName);
+
+ return canvas;
+}
+
+TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE)
+{
+ // compares n results with first results. E.g. one gained with the default response, another with a changed one to study
+ // a systematic effect
+
+ // normalize results
+ result->Scale(1.0 / result->Integral(2, 200));
+
+ TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
+
+ result->GetXaxis()->SetRangeUser(0, 150);
+ result->SetStats(kFALSE);
+
+ for (Int_t n=0; n<nResultSyst; ++n)
+ {
+ resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
+
+ // calculate ratio
+ TH1* ratio = (TH1*) result->Clone("ratio");
+ ratio->Divide(ratio, resultSyst[n], 1, 1, "B");
+ ratio->SetTitle(";true multiplicity;Ratio");
+ ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
+
+ if (firstMarker)
+ ratio->SetMarkerStyle(5);
+
+ ratio->SetLineColor(n+1);
+
+ ratio->DrawCopy((n == 0) ? ((firstMarker) ? "P" : "HIST") : "SAME HIST");
+
+ // get average of ratio
+ Float_t sum = 0;
+ for (Int_t i=2; i<=150; ++i)
+ sum += TMath::Abs(ratio->GetBinContent(i) - 1);
+ sum /= 149;
+
+ printf("%d) Average (2..150) of |ratio - 1| is %f\n", n, sum);
+ }
+
+ TLine* line = new TLine(0, 1, 150, 1);
+ line->SetLineWidth(2);
+ line->Draw();
+
+ line = new TLine(0, 1.1, 150, 1.1);
+ line->SetLineWidth(2);
+ line->SetLineStyle(2);
+ line->Draw();
+ line = new TLine(0, 0.9, 150, 0.9);
+ line->SetLineWidth(2);
+ line->SetLineStyle(2);
+ line->Draw();
+
+ canvas->Modified();
+
+ canvas->SaveAs(epsName);
+ canvas->SaveAs(Form("%s.gif", epsName.Data()));
+
+ return canvas;
+}
+
+TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
+{
+ // draws the ratios of each mc to the corresponding result
+
+ TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
+
+ for (Int_t n=0; n<nResultSyst; ++n)
+ {
+ // normalize
+ result[n]->Scale(1.0 / result[n]->Integral(2, 200));
+ mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
+
+ result[n]->GetXaxis()->SetRangeUser(0, 150);
+ result[n]->SetStats(kFALSE);
+
+ // calculate ratio
+ TH1* ratio = (TH1*) result[n]->Clone("ratio");
+ ratio->Divide(mc[n], ratio, 1, 1, "B");
+ ratio->SetTitle(";true multiplicity;Ratio (true / unfolded)");
+ ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
+
+ ratio->SetLineColor(n+1);
+
+ ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
+
+ // get average of ratio
+ Float_t sum = 0;
+ for (Int_t i=2; i<=150; ++i)
+ sum += TMath::Abs(ratio->GetBinContent(i) - 1);
+ sum /= 149;
+
+ printf("%d) Average (2..150) of |ratio - 1| is %f\n", n, sum);
+ }
+
+ TLine* line = new TLine(0, 1, 150, 1);
+ line->SetLineWidth(2);
+ line->Draw();
+
+ line = new TLine(0, 1.1, 150, 1.1);
+ line->SetLineWidth(2);
+ line->SetLineStyle(2);
+ line->Draw();
+ line = new TLine(0, 0.9, 150, 0.9);
+ line->SetLineWidth(2);
+ line->SetLineStyle(2);
+ line->Draw();
+
+ canvas->Modified();
+
+ canvas->SaveAs(epsName);
+ canvas->SaveAs(Form("%s.gif", epsName.Data()));
+
+ return canvas;
+}
+
+TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
+{
+ // draws the ratios of each mc to the corresponding result
+ // deducts from each ratio the ratio of mcBase / resultBase
+
+ // normalize
+ resultBase->Scale(1.0 / resultBase->Integral(2, 200));
+ mcBase->Scale(1.0 / mcBase->Integral(2, 200));
+
+ // calculate ratio
+ TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
+ ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
+
+ TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
+ canvas->SetRightMargin(0.05);
+ canvas->SetTopMargin(0.05);
+
+ for (Int_t n=0; n<nResultSyst; ++n)
+ {
+ // normalize
+ result[n]->Scale(1.0 / result[n]->Integral(2, 200));
+ mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
+
+ result[n]->GetXaxis()->SetRangeUser(0, 150);
+ result[n]->SetStats(kFALSE);
+
+ // calculate ratio
+ TH1* ratio = (TH1*) result[n]->Clone("ratio");
+ ratio->Divide(mc[n], ratio, 1, 1, "B");
+ ratio->Add(ratioBase, -1);
+
+ ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u)");
+ ratio->GetYaxis()->SetRangeUser(-1, 1);
+ ratio->SetLineColor(n+1);
+ ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
+
+ // get average of ratio
+ Float_t sum = 0;
+ for (Int_t i=2; i<=150; ++i)
+ sum += TMath::Abs(ratio->GetBinContent(i));
+ sum /= 149;
+
+ printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
+ }
+
+ TLine* line = new TLine(0, 0, 150, 0);
+ line->SetLineWidth(2);
+ line->Draw();
+
+ line = new TLine(0, 0.1, 150, 0.1);
+ line->SetLineWidth(2);
+ line->SetLineStyle(2);
+ line->Draw();
+ line = new TLine(0, -0.1, 150, -0.1);
+ line->SetLineWidth(2);
+ line->SetLineStyle(2);
+ line->Draw();
+
+ canvas->Modified();
+
+ canvas->SaveAs(epsName);
+ canvas->SaveAs(Form("%s.gif", epsName.Data()));
+
+ return canvas;
+}
+
+void Smooth(TH1* hist, Int_t windowWidth = 20)
+{
+ TH1* clone = (TH1*) hist->Clone("clone");
+ for (Int_t bin=3; bin<=clone->GetNbinsX(); ++bin)
+ {
+ Int_t min = TMath::Max(3, bin-windowWidth);
+ Int_t max = TMath::Min(clone->GetNbinsX(), bin+windowWidth);
+ Float_t average = clone->Integral(min, max) / (max - min + 1);
+
+ hist->SetBinContent(bin, average);
+ hist->SetBinError(bin, 0);
+ }
+
+ delete clone;
+}
+
+TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
+{
+ // draws the ratios of each mc to the corresponding result
+ // deducts from each ratio the ratio of mcBase / resultBase
+ // smoothens the ratios by a sliding window
+
+ // normalize
+ resultBase->Scale(1.0 / resultBase->Integral(2, 200));
+ mcBase->Scale(1.0 / mcBase->Integral(2, 200));
+
+ // calculate ratio
+ TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
+ ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
+
+ TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
+ canvas->SetRightMargin(0.05);
+ canvas->SetTopMargin(0.05);
+
+ for (Int_t n=0; n<nResultSyst; ++n)
+ {
+ // normalize
+ result[n]->Scale(1.0 / result[n]->Integral(2, 200));
+ mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
+
+ result[n]->GetXaxis()->SetRangeUser(0, 150);
+ result[n]->SetStats(kFALSE);
+
+ // calculate ratio
+ TH1* ratio = (TH1*) result[n]->Clone("ratio");
+ ratio->Divide(mc[n], ratio, 1, 1, "B");
+ ratio->Add(ratioBase, -1);
+
+ new TCanvas; ratio->DrawCopy();
+ Smooth(ratio);
+ ratio->SetLineColor(1);
+ ratio->DrawCopy("SAME");
+
+ canvas->cd();
+ ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u)");
+ ratio->GetYaxis()->SetRangeUser(-1, 1);
+ ratio->SetLineColor(n+1);
+ ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
+
+ // get average of ratio
+ Float_t sum = 0;
+ for (Int_t i=2; i<=150; ++i)
+ sum += TMath::Abs(ratio->GetBinContent(i));
+ sum /= 149;
+
+ printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
+ }
+
+ TLine* line = new TLine(0, 0, 150, 0);
+ line->SetLineWidth(2);
+ line->Draw();
+
+ line = new TLine(0, 0.1, 150, 0.1);
+ line->SetLineWidth(2);
+ line->SetLineStyle(2);
+ line->Draw();
+ line = new TLine(0, -0.1, 150, -0.1);
+ line->SetLineWidth(2);
+ line->SetLineStyle(2);
+ line->Draw();
+
+ canvas->Modified();
+
+ canvas->SaveAs(epsName);
+ canvas->SaveAs(Form("%s.gif", epsName.Data()));
+
+ return canvas;
+}
+
+void DrawResiduals(TH1* measured, TH1* unfoldedFolded, const char* epsName)
+{
+ // normalize
+ unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, 200));
+ unfoldedFolded->Scale(measured->Integral(2, 200));
+
+ TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
+ canvas->Range(0, 0, 1, 1);
+
+ TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
+ pad1->Draw();
+ pad1->SetGridx();
+ pad1->SetGridy();
+
+ TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
+ pad2->Draw();
+ pad2->SetGridx();
+ pad2->SetGridy();
+
+ TPad* pad3 = new TPad(Form("%s_pad3", epsName), "", 0.15, 0.5, 0.35, 0.75);
+ pad3->SetGridx();
+ pad3->SetGridy();
+ pad3->SetRightMargin(0.05);
+ pad3->SetTopMargin(0.05);
+ pad3->Draw();
+
+ pad1->SetRightMargin(0.05);
+ pad2->SetRightMargin(0.05);
+
+ // no border between them
+ pad1->SetBottomMargin(0);
+ pad2->SetTopMargin(0);
+
+ pad1->cd();
+
+ measured->GetXaxis()->SetLabelSize(0.06);
+ measured->GetYaxis()->SetLabelSize(0.06);
+ measured->GetXaxis()->SetTitleSize(0.06);
+ measured->GetYaxis()->SetTitleSize(0.06);
+ measured->GetYaxis()->SetTitleOffset(0.6);
+
+ measured->GetXaxis()->SetRangeUser(0, 150);
+
+ measured->SetTitle(";measured multiplicity;Entries");
+ measured->SetStats(kFALSE);
+
+ measured->DrawCopy("HIST");
+ gPad->SetLogy();
+
+ unfoldedFolded->SetMarkerStyle(5);
+ unfoldedFolded->SetMarkerColor(2);
+ unfoldedFolded->SetLineColor(0);
+ unfoldedFolded->DrawCopy("SAME P");
+
+ TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
+ legend->AddEntry(measured, "measured distribution");
+ legend->AddEntry(unfoldedFolded, "R #otimes unfolded distribution");
+ legend->SetFillColor(0);
+ legend->Draw();
+
+ pad2->cd();
+ pad2->SetBottomMargin(0.15);
+
+ // calculate ratio
+ measured->Sumw2();
+ TH1* residual = (TH1*) measured->Clone("residual");
+ unfoldedFolded->Sumw2();
+
+ residual->Add(unfoldedFolded, -1);
+
+ // projection
+ TH1* residualHist = new TH1F("residualHist", ";", 15, -3, 3);
+
+ for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
+ {
+ if (measured->GetBinError(i) > 0)
+ {
+ residual->SetBinContent(i, residual->GetBinContent(i) / measured->GetBinError(i));
+ residual->SetBinError(i, 1);
+
+ residualHist->Fill(residual->GetBinContent(i));
+ }
+ else
+ {
+ residual->SetBinContent(i, 0);
+ residual->SetBinError(i, 0);
+ }
+ }
+
+ residual->GetYaxis()->SetTitle("Residuals 1/e (M - R #otimes U)");
+ residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
+ residual->DrawCopy();
+
+ TLine* line = new TLine(-0.5, 0, 150.5, 0);
+ line->SetLineWidth(2);
+ line->Draw();
+
+ pad3->cd();
+ residualHist->SetStats(kFALSE);
+ residualHist->GetXaxis()->SetLabelSize(0.08);
+ residualHist->Fit("gaus");
+ residualHist->Draw();
+
+ canvas->Modified();
+ canvas->SaveAs(canvas->GetName());
+
+ //const char* epsName2 = "proj.eps";
+ //TCanvas* canvas = new TCanvas(epsName2, epsName2, 800, 600);
+ //canvas->SetGridx();
+ //canvas->SetGridy();
+
+ //canvas->SaveAs(canvas->GetName());
+}
+
+void bayesianExample()
+{
+ TStopwatch watch;
+ watch.Start();
+
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(correctionFile);
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ TFile::Open(measuredFile);
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+ mult2->LoadHistograms("Multiplicity");
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+
+ mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
+
+ TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
+ TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
+
+ //mult->DrawComparison("bayesianExample", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
+ DrawResultRatio(mcHist, result, "bayesianExample.eps");
+
+ // draw residual plot
+
+ // TODO take out efficiency correction if other than AliMultiplicityCorrection::kTrVtx
+ TH2* convoluted = mult->CalculateMultiplicityESD(result, etaRange);
+ TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
+
+ TH1* measured = mult2->GetMultiplicityESD(etaRange)->ProjectionY("measured");
+
+ DrawResiduals(measured, convolutedProj, "bayesianResiduals.eps");
+
+ watch.Stop();
+ watch.Print();
+}
+
+void chi2FluctuationResult()
+{
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(correctionFile);
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ TFile::Open(measuredFile);
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+ mult2->LoadHistograms("Multiplicity");
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0);
+ mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+
+ TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
+ TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
+
+ mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
+
+ TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_3");
+ canvas->SaveAs("chi2FluctuationResult.eps");
+}
+
+void chi2Example()
+{
+ TStopwatch watch;
+ watch.Start();
+
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(correctionFile);
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ TFile::Open(measuredFile);
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+ mult2->LoadHistograms("Multiplicity");
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+ mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+
+ TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
+ TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
+
+ DrawResultRatio(mcHist, result, "chi2Example.eps");
+
+ watch.Stop();
+ watch.Print();
+}
+
+void chi2ExampleTPC()
+{
+ TStopwatch watch;
+ watch.Start();
+
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(correctionFileTPC);
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ TFile::Open(measuredFileTPC);
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+ mult2->LoadHistograms("Multiplicity");
+
+ mult->SetMultiplicityESD(etaRangeTPC, mult2->GetMultiplicityESD(etaRangeTPC));
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+ mult->ApplyMinuitFit(etaRangeTPC, kFALSE, AliMultiplicityCorrection::kTrVtx);
+
+ TH1* mcHist = mult2->GetMultiplicityVtx(etaRangeTPC)->ProjectionY("mymc");
+ TH1* result = mult->GetMultiplicityESDCorrected(etaRangeTPC);
+
+ DrawResultRatio(mcHist, result, "chi2ExampleTPC.eps");
+
+ watch.Stop();
+ watch.Print();
+}
+
+void bayesianNBD()
+{
+ gSystem->Load("libPWG0base");
+ TFile::Open("multiplicityMC_3M.root");
+
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ TFile::Open("multiplicityMC_3M_NBD.root");
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+ mult2->LoadHistograms("Multiplicity");
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+ mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1.0, 100);
+
+ TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
+
+ mcHist->Sumw2();
+ TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
+
+ //mult->DrawComparison("bayesianNBD", etaRange, kFALSE, kTRUE, mcHist);
+ DrawResultRatio(mcHist, result, "bayesianNBD.eps");
+}
+
+void minimizationNBD()
+{
+ gSystem->Load("libPWG0base");
+ TFile::Open("multiplicityMC_3M.root");
+
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ TFile::Open("multiplicityMC_3M_NBD.root");
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+ mult2->LoadHistograms("Multiplicity");
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+ mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+
+ TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
+
+ mcHist->Sumw2();
+ TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
+
+ //mult->DrawComparison("minimizationNBD", etaRange, kFALSE, kTRUE, mcHist);
+ DrawResultRatio(mcHist, result, "minimizationNBD.eps");
+}
+
+void minimizationInfluenceAlpha()
+{
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(measuredFile);
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+ mult2->LoadHistograms("Multiplicity");
+
+ TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
+ mcHist->Scale(1.0 / mcHist->Integral());
+ mcHist->GetXaxis()->SetRangeUser(0, 200);
+ mcHist->SetStats(kFALSE);
+ mcHist->SetTitle(";true multiplicity;P_{N}");
+
+ TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 300);
+ canvas->Divide(3, 1);
+
+ TFile::Open("eval-2M-1M/EvaluateChi2MethodDetail.root");
+
+ TH1* hist1 = gFile->Get("MinuitChi2_00_2_100.000000");
+ TH1* hist2 = gFile->Get("MinuitChi2_03_2_100000.000000");
+ TH1* hist3 = gFile->Get("MinuitChi2_06_2_100000000.000000");
+
+ mcHist->Rebin(2); mcHist->Scale(0.5);
+ hist1->Rebin(2); hist1->Scale(0.5);
+ hist2->Rebin(2); hist2->Scale(0.5);
+ hist3->Rebin(2); hist3->Scale(0.5);
+
+ mcHist->GetXaxis()->SetRangeUser(0, 200);
+
+ canvas->cd(1);
+ gPad->SetLogy();
+ mcHist->Draw();
+ hist1->SetMarkerStyle(5);
+ hist1->SetMarkerColor(2);
+ hist1->Draw("SAME PE");
+
+ canvas->cd(2);
+ gPad->SetLogy();
+ mcHist->Draw();
+ hist2->SetMarkerStyle(5);
+ hist2->SetMarkerColor(2);
+ hist2->Draw("SAME PE");
+
+ canvas->cd(3);
+ gPad->SetLogy();
+ mcHist->Draw();
+ hist3->SetMarkerStyle(5);
+ hist3->SetMarkerColor(2);
+ hist3->Draw("SAME PE");
+
+ canvas->SaveAs("minimizationInfluenceAlpha.eps");
+}
+
+void NBDFit()
+{
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(correctionFile);
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ TH1* fCurrentESD = mult->GetMultiplicityVtx(etaRange)->ProjectionY();
+ fCurrentESD->Sumw2();
+ fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
+
+ TF1* 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])");
+ 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);
+ func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
+
+ TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
+ lognormal->SetParNames("scaling", "mean", "sigma");
+ lognormal->SetParameters(1, 1, 1);
+ lognormal->SetParLimits(0, 0, 10);
+ lognormal->SetParLimits(1, 0, 100);
+ lognormal->SetParLimits(2, 1e-3, 10);
+
+ TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
+ fCurrentESD->SetStats(kFALSE);
+ fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
+ fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
+ fCurrentESD->Draw("HIST");
+ fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
+ fCurrentESD->Fit(func, "W0", "", 0, 50);
+ func->SetRange(0, 100);
+ func->Draw("SAME");
+ printf("chi2 = %f\n", func->GetChisquare());
+
+ fCurrentESD->Fit(lognormal, "W0", "", 0.01, 100);
+ lognormal->SetLineColor(2);
+ lognormal->SetLineStyle(2);
+ lognormal->SetRange(0, 100);
+ lognormal->Draw("SAME");
+
+ canvas->SaveAs("NBDFit.eps");
+}
+
+void DifferentSamples()
+{
+ // data generated by runMultiplicitySelector.C DifferentSamples
+
+ const char* name = "DifferentSamples";
+
+ TFile* file = TFile::Open(Form("%s.root", name));
+
+ TCanvas* canvas = new TCanvas(name, name, 800, 600);
+ canvas->Divide(2, 2);
+
+ for (Int_t i=0; i<4; ++i)
+ {
+ canvas->cd(i+1);
+ gPad->SetTopMargin(0.05);
+ gPad->SetRightMargin(0.05);
+ TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", i));
+ TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", i));
+ TH1* mc = (TH1*) file->Get(Form("mc_%d", i));
+ mc->Sumw2();
+
+ chi2Result->Divide(chi2Result, mc, 1, 1, "");
+ bayesResult->Divide(bayesResult, mc, 1, 1, "");
+
+ chi2Result->SetTitle(";true multiplicity;unfolded measured/MC");
+ chi2Result->GetXaxis()->SetRangeUser(0, 150);
+ chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
+ chi2Result->GetYaxis()->SetTitleOffset(1.2);
+ chi2Result->SetLineColor(1);
+ chi2Result->SetStats(kFALSE);
+
+ bayesResult->SetStats(kFALSE);
+ bayesResult->SetLineColor(2);
+
+ chi2Result->DrawCopy("HIST");
+ bayesResult->DrawCopy("SAME HIST");
+
+ TLine* line = new TLine(0, 1, 150, 1);
+ line->Draw();
+ }
+
+ canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+}
+
+void StartingConditions()
+{
+ // data generated by runMultiplicitySelector.C StartingConditions
+
+ const char* name = "StartingConditions";
+
+ TFile* file = TFile::Open(Form("%s.root", name));
+
+ TCanvas* canvas = new TCanvas(name, name, 800, 400);
+ canvas->Divide(2, 1);
+
+ TH1* mc = (TH1*) file->Get("mc");
+ mc->Sumw2();
+ mc->Scale(1.0 / mc->Integral());
+
+ Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
+
+ for (Int_t i=0; i<6; ++i)
+ {
+ Int_t id = i;
+ if (id > 2)
+ id += 2;
+
+ TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", id));
+ TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", id));
+
+ chi2Result->Divide(chi2Result, mc, 1, 1, "");
+ bayesResult->Divide(bayesResult, mc, 1, 1, "");
+
+ chi2Result->SetTitle("a) #chi^{2} minimization;true multiplicity;unfolded / MC");
+ chi2Result->GetXaxis()->SetRangeUser(0, 150);
+ chi2Result->GetYaxis()->SetRangeUser(0.8, 1.2);
+ chi2Result->GetYaxis()->SetTitleOffset(1.5);
+ //chi2Result->SetMarkerStyle(marker[i]);
+ chi2Result->SetLineColor(i+1);
+ chi2Result->SetStats(kFALSE);
+
+ bayesResult->SetTitle("b) Bayesian method;true multiplicity;unfolded / MC");
+ bayesResult->GetXaxis()->SetRangeUser(0, 150);
+ bayesResult->GetYaxis()->SetRangeUser(0.8, 1.2);
+ bayesResult->GetYaxis()->SetTitleOffset(1.5);
+ bayesResult->SetStats(kFALSE);
+ bayesResult->SetLineColor(2);
+ bayesResult->SetLineColor(i+1);
+
+ canvas->cd(1);
+ gPad->SetLeftMargin(0.12);
+ chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
+
+ canvas->cd(2);
+ gPad->SetLeftMargin(0.12);
+ bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
+
+ //TLine* line = new TLine(0, 1, 150, 1);
+ //line->Draw();
+ }
+
+ canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+}
+
+void StatisticsPlot()
+{
+ const char* name = "StatisticsPlot";
+
+ TFile* file = TFile::Open(Form("%s.root", name));
+
+ TCanvas* canvas = new TCanvas(name, name, 600, 400);
+
+ TGraph* fitResultsChi2 = file->Get("fitResultsChi2");
+ fitResultsChi2->SetTitle(";number of measured events;P_{1}");
+ fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
+ fitResultsChi2->Draw("AP");
+
+ TF1* f = new TF1("f", "[0]/x", 1, 1e4);
+ fitResultsChi2->Fit(f);
+
+ canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+
+ TH1* mc[5];
+ TH1* result[5];
+
+ const char* plotname = "chi2Result";
+
+ const char* name = "StatisticsPlotRatios";
+
+ TCanvas* canvas = new TCanvas(name, name, 600, 400);
+
+ for (Int_t i=0; i<5; ++i)
+ {
+ mc[i] = (TH1*) file->Get(Form("mc_%d", i));
+ result[i] = (TH1*) file->Get(Form("%s_%d", plotname, i));
+
+ result[i]->SetLineColor(i+1);
+ result[i]->Draw(((i == 0) ? "" : "SAME"));
+ }
+}
+
+void SystematicLowEfficiency()
+{
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(correctionFile);
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ // calculate result with systematic effect
+ TFile::Open("multiplicityMC_100k_1_loweff.root");
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+ mult2->LoadHistograms("Multiplicity");
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+ mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+
+ TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
+ TH1* result1 = mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
+
+ DrawResultRatio(mcHist, result1, "SystematicLowEfficiencyLow.eps");
+
+ // calculate normal result
+ TFile::Open("multiplicityMC_100k_1.root");
+ mult2->LoadHistograms("Multiplicity");
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+ mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+
+ TH1* result2 = mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
+
+ DrawResultRatio(mcHist, result2, "SystematicLowEfficiencyOK.eps");
+
+ Draw2ResultRatio(mcHist, result1, result2, "SystematicLowEfficiency.eps");
+}
+
+void SystematicMisalignment()
+{
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(correctionFile);
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ TFile::Open("multiplicityMC_100k_fullmis.root");
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+ mult2->LoadHistograms("Multiplicity");
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+ mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+
+ TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
+ TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
+
+ DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
+}
+
+void EfficiencySpecies(const char* fileName = "multiplicityMC_400k_syst.root")
+{
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(fileName);
+ AliCorrection* correction[4];
+
+ TCanvas* canvas = new TCanvas("EfficiencySpeciesSPD", "EfficiencySpeciesSPD", 600, 400);
+ canvas->SetGridx();
+ canvas->SetGridy();
+ canvas->SetRightMargin(0.05);
+ canvas->SetTopMargin(0.05);
+
+ TLegend* legend = new TLegend(0.7, 0.4, 0.85, 0.6);
+ legend->SetFillColor(0);
+
+ Int_t marker[] = {24, 25, 26};
+ Int_t color[] = {1, 2, 4};
+
+ Float_t below = 0;
+ Float_t total = 0;
+
+ for (Int_t i=0; i<3; ++i)
+ {
+ Printf("correction %d", i);
+
+ TString name; name.Form("correction_%d", i);
+ correction[i] = new AliCorrection(name, name);
+ correction[i]->LoadHistograms();
+
+ TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
+ TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
+
+ // limit vtx axis
+ gene->GetXaxis()->SetRangeUser(-3.9, 3.9);
+ meas->GetXaxis()->SetRangeUser(-3.9, 3.9);
+
+ // empty over/underflow bin in eta, setting range to +-2 is not enough because this is the maximum range, Project3D takes them into account then (might be a bug)
+ /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
+ for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
+ {
+ gene->SetBinContent(x, 0, z, 0);
+ gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
+ meas->SetBinContent(x, 0, z, 0);
+ meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
+ }*/
+
+ // limit eta axis
+ gene->GetYaxis()->SetRangeUser(-0.49, 0.49);
+ meas->GetYaxis()->SetRangeUser(-0.49, 0.49);
+
+ TH1* genePt = gene->Project3D(Form("z_%d", i));
+ TH1* measPt = meas->Project3D(Form("z_%d", i));
+
+ genePt->Sumw2();
+ measPt->Sumw2();
+
+ effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
+ effPt->Reset();
+ effPt->Divide(measPt, genePt, 1, 1, "B");
+
+ Int_t bin = 0;
+ for (bin=20; bin>=1; bin--)
+ {
+ if (effPt->GetBinContent(bin) < 0.5)
+ break;
+ }
+
+ Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
+
+ Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
+ Printf("%.4f of the particles are below that momentum", fraction);
+
+ below += genePt->Integral(1, bin);
+ total += genePt->Integral();
+
+ effPt->SetLineColor(color[i]);
+ effPt->SetMarkerColor(color[i]);
+ effPt->SetMarkerStyle(marker[i]);
+
+ effPt->GetXaxis()->SetRangeUser(0, 1);
+ effPt->GetYaxis()->SetRangeUser(0, 1);
+
+ effPt->SetStats(kFALSE);
+ effPt->SetTitle("");
+ effPt->GetYaxis()->SetTitle("Efficiency");
+
+ effPt->DrawCopy((i == 0) ? "" : "SAME");
+
+ legend->AddEntry(effPt, ((i == 0) ? "#pi" : ((i == 1) ? "K" : "p")));
+ }
+
+ Printf("In total %.4f of the particles are below their pt cut off", (Float_t) below / total);
+
+ legend->Draw();
+
+ canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+}
+
+void ParticleSpeciesComparison1()
+{
+ gSystem->Load("libPWG0base");
+
+ const char* fileNameMC = "multiplicityMC_400k_syst_species.root";
+ const char* fileNameESD = "multiplicityMC_100k_syst.root";
+ Bool_t chi2 = 0;
+
+ TFile::Open(fileNameESD);
+ TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", etaRange));
+ TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", etaRange));
+
+ TH1* results[10];
+
+ // loop over cases (normal, enhanced/reduced ratios)
+ Int_t nMax = 7;
+ for (Int_t i = 0; i<nMax; ++i)
+ {
+ TString folder;
+ folder.Form("Multiplicity_%d", i);
+
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
+
+ TFile::Open(fileNameMC);
+ mult->LoadHistograms();
+
+ mult->SetMultiplicityESD(etaRange, hist);
+
+ if (chi2)
+ {
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
+ mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
+ //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
+ }
+ else
+ {
+ mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
+ //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
+ }
+
+ //Float_t averageRatio = 0;
+ //mult->GetComparisonResults(0, 0, 0, &averageRatio);
+
+ results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
+
+ //Printf("Case %d. Average ratio is %f", i, averageRatio);
+ }
+
+ DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[0], "ParticleSpeciesComparison1_1.eps");
+
+ TH1* mc = hist2->ProjectionY("mymchist2", -1, -1, "e");
+
+ for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
+ {
+ results[0]->SetBinError(i, 0);
+ mc->SetBinError(i, 0);
+ }
+
+ TCanvas* canvas = DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps");
+
+ TFile::Open("bayesianUncertainty_400k_100k_syst.root");
+ TH1* errorHist = (TH1*) gFile->Get("errorBoth");
+ errorHist->SetLineColor(1);
+ errorHist->SetLineWidth(2);
+ TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
+ for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
+ {
+ errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
+ errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
+ }
+
+ errorHist->DrawCopy("SAME");
+ errorHist2->DrawCopy("SAME");
+
+ canvas->SaveAs(canvas->GetName());
+
+ TCanvas* canvas2 = DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE);
+
+ errorHist->DrawCopy("SAME");
+ errorHist2->DrawCopy("SAME");
+
+ canvas2->SaveAs(canvas2->GetName());
+}
+
+void ParticleSpeciesComparison2()
+{
+ gSystem->Load("libPWG0base");
+
+ const char* fileNameMC = "multiplicityMC_400k_syst.root";
+ const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
+ Bool_t chi2 = 0;
+
+ TFile::Open(fileNameMC);
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms();
+
+ TH1* mc[10];
+ TH1* results[10];
+
+ // loop over cases (normal, enhanced/reduced ratios)
+ Int_t nMax = 7;
+ for (Int_t i = 0; i<nMax; ++i)
+ {
+ TString folder;
+ folder.Form("Multiplicity_%d", i);
+
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
+
+ TFile::Open(fileNameESD);
+ mult2->LoadHistograms();
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+
+ if (chi2)
+ {
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
+ mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
+ //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
+ }
+ else
+ {
+ mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
+ //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
+ }
+
+ //Float_t averageRatio = 0;
+ //mult->GetComparisonResults(0, 0, 0, &averageRatio);
+
+ results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
+
+ TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
+ mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
+
+ //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
+ //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
+
+ //Printf("Case %d. Average ratio is %f", i, averageRatio);
+ }
+
+ DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
+}
+
+TH1* Invert(TH1* eff)
+{
+ // calculate corr = 1 / eff
+
+ TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
+ corr->Reset();
+
+ for (Int_t i=1; i<=eff->GetNbinsX(); i++)
+ {
+ if (eff->GetBinContent(i) > 0)
+ {
+ corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
+ corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
+ }
+ }
+
+ return corr;
+}
+
+void TriggerVertexCorrection()
+{
+ //
+ // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
+ //
+
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(correctionFile);
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
+ TH1* corrMB = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
+
+ TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 600);
+
+ corrINEL->SetStats(kFALSE);
+ corrINEL->GetXaxis()->SetRangeUser(0, 30);
+ corrINEL->GetYaxis()->SetRangeUser(0, 5);
+ corrINEL->SetTitle(";true multiplicity;correction factor");
+ corrINEL->SetMarkerStyle(22);
+ corrINEL->Draw("PE");
+
+ corrMB->SetStats(kFALSE);
+ corrMB->SetLineColor(2);
+ corrMB->SetMarkerStyle(25);
+ corrMB->SetMarkerColor(2);
+ corrMB->Draw("SAME PE");
+
+ TLegend* legend = new TLegend(0.3, 0.5, 0.85, 0.65);
+ legend->SetFillColor(0);
+ legend->AddEntry(corrINEL, "correction to inelastic sample");
+ legend->AddEntry(corrMB, "correction to minimum bias sample");
+
+ legend->Draw();
+
+ canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+}
+
+void bayesianUncertainty(Bool_t mc = kFALSE)
+{
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(correctionFile);
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ TFile::Open(measuredFile);
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+ mult2->LoadHistograms("Multiplicity");
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+
+ TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
+
+ TH1* errorResponse = mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, 1, 100, ((mc) ? mcHist : 0))->Clone("errorResponse");
+ TH1* errorMeasured = mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, 1, 100, ((mc) ? mcHist : 0))->Clone("errorMeasured");
+ TH1* errorBoth = mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, 1, 100, ((mc) ? mcHist : 0))->Clone("errorBoth");
+
+ if (!mc)
+ {
+ TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
+ DrawResultRatio(mcHist, result, "bayesianUncertainty2.eps");
+ }
+
+ TCanvas* canvas = new TCanvas("bayesianUncertainty", "bayesianUncertainty", 600, 400);
+ canvas->SetGridx();
+ canvas->SetGridy();
+ canvas->SetRightMargin(0.05);
+ canvas->SetTopMargin(0.05);
+
+ errorResponse->SetLineColor(1);
+ errorResponse->GetXaxis()->SetRangeUser(0, 200);
+ errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
+ errorResponse->SetStats(kFALSE);
+ errorResponse->SetTitle(";true multiplicity;Uncertainty");
+
+ errorResponse->Draw();
+
+ errorMeasured->SetLineColor(2);
+ errorMeasured->Draw("SAME");
+
+ errorBoth->SetLineColor(3);
+ errorBoth->Draw("SAME");
+
+ Printf("Average errorResponse: %f", errorResponse->Integral(2, 150) / 149);
+ Printf("Average errorMeasured: %f", errorMeasured->Integral(2, 150) / 149);
+ Printf("Average errorBoth: %f", errorBoth->Integral(2, 150) / 149);
+
+ canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+
+ TFile* file = new TFile(Form("%s.root", canvas->GetName()), "RECREATE");
+ errorResponse->Write();
+ errorMeasured->Write();
+ errorBoth->Write();
+ file->Close();
+}
+
+void EfficiencyComparison(Int_t eventType = 2)
+{
+ const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root", "multiplicityMC_400k_syst_xsection.root" };
+
+ gSystem->Load("libPWG0base");
+
+ TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 500);
+ canvas->SetGridx();
+ canvas->SetGridy();
+ canvas->SetRightMargin(0.05);
+ canvas->SetTopMargin(0.05);
+
+ AliMultiplicityCorrection* data[4];
+ TH1* effArray[4];
+
+ Int_t markers[] = { 24, 25, 26, 5 };
+ Int_t colors[] = { 1, 2, 3, 4 };
+
+ TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
+ legend->SetFillColor(0);
+
+ TH1* effError = 0;
+
+ for (Int_t i=0; i<4; ++i)
+ {
+ TString name;
+ name.Form("Multiplicity_%d", i);
+
+ TFile::Open(files[i]);
+ data[i] = new AliMultiplicityCorrection(name, name);
+
+ if (i < 3)
+ {
+ data[i]->LoadHistograms("Multiplicity");
+ }
+ else
+ data[i]->LoadHistograms("Multiplicity_0");
+
+ TH1* eff = (TH1*) data[i]->GetEfficiency(3, eventType)->Clone(Form("eff_%d", i));
+ effArray[i] = eff;
+
+ eff->GetXaxis()->SetRangeUser(0, 15);
+ eff->GetYaxis()->SetRangeUser(0, 1.1);
+ eff->SetStats(kFALSE);
+ eff->SetTitle(";true multiplicity;Efficiency");
+ eff->SetLineColor(colors[i]);
+ eff->SetMarkerColor(colors[i]);
+ eff->SetMarkerStyle(markers[i]);
+
+ if (i == 3)
+ {
+ for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
+ eff->SetBinError(bin, 0);
+
+ // loop over cross section combinations
+ for (Int_t j=1; j<7; ++j)
+ {
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
+ mult->LoadHistograms(Form("Multiplicity_%d", j));
+
+ TH1* eff2 = mult->GetEfficiency(3, eventType);
+
+ for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
+ {
+ // TODO we could also do assymetric errors here
+ Float_t deviation = 10.0 * TMath::Abs(eff->GetBinContent(bin) - eff2->GetBinContent(bin));
+
+ eff->SetBinError(bin, TMath::Max(eff->GetBinError(bin), deviation));
+ }
+ }
+
+ for (Int_t bin=1; bin<=20; bin++)
+ if (eff->GetBinContent(bin) > 0)
+ Printf("Bin %d: Error: %.2f", bin, 100.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
+
+ effError = (TH1*) eff2->Clone("effError");
+ effError->Reset();
+
+ for (Int_t bin=2; bin<=eff->GetNbinsX(); bin++)
+ if (eff->GetBinContent(bin) > 0)
+ effError->SetBinContent(bin, eff->GetBinError(bin) / eff->GetBinContent(bin));
+
+ effError->DrawCopy("SAME HIST");
+ }
+
+ eff->SetBinContent(1, 0);
+ eff->SetBinError(1, 0);
+
+ canvas->cd();
+ if (i == 0)
+ {
+ eff->DrawCopy("P");
+ }
+ else
+ eff->DrawCopy("SAME P");
+
+ legend->AddEntry(eff, (((i == 0) ? "non diffractive" : ((i == 1) ? "single diffractive" : ((i == 2) ? "double diffractive" : "Pythia combined")))));
+ }
+
+ legend->AddEntry(effError, "relative syst. uncertainty #times 10");
+
+ legend->Draw();
+
+ canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+}
+
+void ModelDependencyPlot()
+{
+ gSystem->Load("libPWG0base");
+
+ TFile::Open("multiplicityMC_3M.root");
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ TH2* proj = (TH2*) mult->GetCorrelation(3)->Project3D("zy");
+
+ TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 800, 400);
+ canvas->SetGridx();
+ canvas->SetGridy();
+ //canvas->SetRightMargin(0.05);
+ //canvas->SetTopMargin(0.05);
+
+ canvas->Divide(2, 1);
+
+ canvas->cd(2);
+ gPad->SetLogy();
+
+ Int_t selectedMult = 30;
+ Int_t yMax = 200000;
+
+ TH1* full = proj->ProjectionX("full");
+ TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult));
+
+ full->SetStats(kFALSE);
+ full->GetXaxis()->SetRangeUser(0, 200);
+ full->GetYaxis()->SetRangeUser(5, yMax);
+ full->SetTitle(";multiplicity");
+
+ selected->SetLineColor(0);
+ selected->SetMarkerColor(2);
+ selected->SetMarkerStyle(7);
+
+ full->Draw();
+ selected->Draw("SAME P");
+
+ TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
+ legend->SetFillColor(0);
+ legend->AddEntry(full, "true");
+ legend->AddEntry(selected, "measured");
+ legend->Draw();
+
+ TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
+ line->SetLineWidth(2);
+ line->Draw();
+
+ canvas->cd(1);
+ gPad->SetLogy();
+
+ TH1* full = proj->ProjectionY("full2");
+ TH1* selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
+
+ full->SetStats(kFALSE);
+ full->GetXaxis()->SetRangeUser(0, 200);
+ full->GetYaxis()->SetRangeUser(5, yMax);
+ full->SetTitle(";multiplicity");
+
+ full->SetLineColor(0);
+ full->SetMarkerColor(2);
+ full->SetMarkerStyle(7);
+
+ full->Draw("P");
+ selected->Draw("SAME");
+
+ legend->Draw();
+
+ TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
+ line->SetLineWidth(2);
+ line->Draw();
+
+ canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+}
+
+void SystematicpTSpectrum()
+{
+ gSystem->Load("libPWG0base");
+
+ TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ TFile::Open("multiplicityMC_100k_syst.root");
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+ mult2->LoadHistograms("Multiplicity");
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+ mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+
+ TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
+ TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
+
+ DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
+}
+
+// to be deleted
+/*void covMatrix(Bool_t mc = kTRUE)
+{
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(correctionFile);
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ TFile::Open(measuredFile);
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+ mult2->LoadHistograms("Multiplicity");
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+
+ TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
+
+ mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, 1, 100, ((mc) ? mcHist : 0));
+}*/
+
+Double_t FitPtFunc(Double_t *x, Double_t *par)
+{
+ Double_t xx = x[0];
+
+ Float_t val1 = par[1] + par[2] * xx + par[3] * xx * xx;
+ Float_t val2 = TMath::Exp(par[4] + par[5] * xx);
+
+ const Float_t kTransitionWidth = 0;
+
+ // power law part
+ if (xx < par[0] - kTransitionWidth)
+ {
+ return val1;
+ }
+ /*else if (xx < par[0] + kTransitionWidth)
+ {
+ // smooth transition
+ Float_t factor = (xx - par[0] + kTransitionWidth) / kTransitionWidth / 2;
+ return (1 - factor) * val1 + factor * val2;
+ }*/
+ else
+ {
+ return val2;
+ }
+}
+
+void FitPt(const char* fileName = "firstplots100k_truept.root")
+{
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(fileName);
+
+ /*
+ // merge corrections
+ AliCorrection* correction[4];
+ TList list;
+
+ for (Int_t i=0; i<4; ++i)
+ {
+ Printf("correction %d", i);
+
+ TString name; name.Form("correction_%d", i);
+ correction[i] = new AliCorrection(name, name);
+ correction[i]->LoadHistograms();
+
+ if (i > 0)
+ list.Add(correction[i]);
+ }
+
+ correction[0]->Merge(&list);
+
+ TH3* gene = correction[0]->GetTrackCorrection()->GetGeneratedHistogram();
+
+ // limit vtx, eta axis
+ gene->GetXaxis()->SetRangeUser(-5.9, 5.9);
+ gene->GetYaxis()->SetRangeUser(-1.99, 0.99);
+
+ TH1* genePt = gene->Project3D("z");*/
+ TH1* genePt = (TH1*) gFile->Get("fdNdpTTrue");
+ genePt->Sumw2();
+
+ //genePt->Scale(1.0 / genePt->Integral());
+
+ // normalize by bin width
+ for (Int_t x=1; x<genePt->GetNbinsX(); x++)
+ genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
+
+ /// genePt->GetXaxis()->GetBinCenter(x));
+
+ genePt->GetXaxis()->SetRangeUser(0, 7.9);
+ genePt->GetYaxis()->SetTitle("a.u.");
+
+ TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x)+[2]/(1+(x*[4])**[3])", 0.001, 100);
+ //func->SetLineColor(2);
+ func->SetParameters(1, -1, 1, 1, 1);
+ func->SetParLimits(3, 1, 10);
+ func->SetParLimits(4, 0, 10);
+
+ //TF1* func = new TF1("func", "[1]*x**[0]", 0.001, 100);
+
+ //TF1* func = new TF1("func", &FitPtFunc, 0, 2, 6);
+ //func->SetParameters(0.3, -2.34909e-01, 1.54394e+01, -3.04134e+01, 1.41912e+00, -2.79284e+00);
+ //func->FixParameter(0, 0.314);
+ //func->SetParLimits(0, 0.1, 0.3);
+
+ genePt->SetMarkerStyle(25);
+ genePt->SetTitle("");
+ genePt->SetStats(kFALSE);
+ genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
+ //func->Draw("SAME");
+
+ // fit only exp. part
+ func->SetParameters(1, -1);
+ func->FixParameter(2, 0);
+ func->FixParameter(3, 1);
+ func->FixParameter(4, 1);
+ genePt->Fit(func, "0", "", 0.2, 1);
+
+ new TCanvas;
+ genePt->DrawCopy("P");
+ func->SetRange(0.02, 8);
+ func->DrawCopy("SAME");
+ gPad->SetLogy();
+
+ // now fix exp. parameters and fit second part
+ Double_t param0 = func->GetParameter(0);
+ Double_t param1 = func->GetParameter(1);
+ func->SetParameters(0, -1, 1, 1, 1);
+ func->FixParameter(0, 0);
+ func->FixParameter(1, -1);
+ func->ReleaseParameter(2);
+ func->ReleaseParameter(3);
+ func->ReleaseParameter(4);
+ func->SetParLimits(3, 1, 10);
+ func->SetParLimits(4, 0, 10);
+
+ genePt->Fit(func, "0", "", 1.5, 4);
+
+ new TCanvas;
+ genePt->DrawCopy("P");
+ func->SetRange(0.02, 8);
+ func->DrawCopy("SAME");
+ gPad->SetLogy();
+
+ // fit both
+ func->SetParameter(0, param0);
+ func->SetParameter(1, param1);
+ func->ReleaseParameter(0);
+ func->ReleaseParameter(1);
+
+ new TCanvas;
+ genePt->DrawCopy("P");
+ func->SetRange(0.02, 8);
+ func->DrawCopy("SAME");
+ gPad->SetLogy();
+
+ genePt->Fit(func, "0", "", 0.2, 4);
+
+ TCanvas* canvas = new TCanvas("FitPt", "FitPt", 600, 400);
+ canvas->SetGridx();
+ canvas->SetGridy();
+ canvas->SetRightMargin(0.05);
+ canvas->SetTopMargin(0.05);
+
+ genePt->DrawCopy("P");
+ func->SetRange(0.02, 8);
+ func->DrawCopy("SAME");
+ gPad->SetLogy();
+
+ TH1* first = (TH1*) func->GetHistogram()->Clone("first");
+
+ TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
+
+ TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
+
+ for (Int_t param=0; param<5; param++)
+ {
+ for (Int_t sign=0; sign<2; sign++)
+ {
+ TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign)); // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
+ func2->SetParameters(func->GetParameters());
+ //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
+
+ Float_t factor = ((sign == 0) ? 0.9 : 1.1);
+ func2->SetParameter(param, func2->GetParameter(param) * factor);
+ //func2->Print();
+
+ canvas->cd();
+ func2->SetLineWidth(1);
+ func2->SetLineColor(param + 2);
+ func2->DrawCopy("SAME");
+
+ canvas2->cd();
+ TH1* second = func2->GetHistogram();
+ second->Divide(first);
+ second->SetLineColor(param + 1);
+ second->GetYaxis()->SetRangeUser(0, 2);
+ second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
+ second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
+ }
+ }
+
+ canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+
+ file->Close();
+}
+
+void SystematicpT()
+{
+ gSystem->Load("libPWG0base");
+
+ TFile::Open("ptspectrum400.root");
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult->LoadHistograms("Multiplicity");
+
+ AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+
+ TH1* mcHist[12];
+ TH1* result[12];
+
+ Int_t nParams = 1;
+
+ for (Int_t param=0; param<nParams; param++)
+ {
+ for (Int_t sign=0; sign<2; sign++)
+ {
+ // calculate result with systematic effect
+ TFile* file = TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
+ mult2->LoadHistograms("Multiplicity");
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+
+ mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+ //mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+
+ Int_t id = param * 2 + sign;
+
+ mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
+ result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
+
+ TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
+ DrawResultRatio(mcHist[id], result[id], tmp);
+ }
+ }
+
+ // calculate normal result
+ TFile::Open("ptspectrum100_1.root");
+ mult2->LoadHistograms("Multiplicity");
+ TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc2");
+
+ mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+
+ mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
+ //mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+
+ TH1* result2 = mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
+
+ DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
+
+ DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
+
+ TCanvas* canvas = DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps");
+
+ TFile::Open("bayesianUncertainty_pt_400_100.root");
+ TH1* errorHist = (TH1*) gFile->Get("errorBoth");
+ errorHist->SetLineColor(1);
+ errorHist->SetLineWidth(2);
+ TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
+ for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
+ {
+ errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
+ errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
+ }
+
+ errorHist->DrawCopy("SAME");
+ errorHist2->DrawCopy("SAME");
+
+ canvas->SaveAs(canvas->GetName());
+
+ DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
+
+ // does not make sense: mc is different
+ //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
+}
void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, const char* option = "")
{
if (aProof)
- connectProof("proof40@lxb6046");
+ connectProof("lxb6046");
TString libraries("libEG;libGeom;libESD;libPWG0base");
TString packages("PWG0base");
}
}
-void draw(const char* fileName = "multiplicityMC.root")
+void SetTPC()
{
gSystem->Load("libPWG0base");
+ AliMultiplicityCorrection::SetQualityRegions(kFALSE);
+}
- AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+void draw(const char* fileName = "multiplicityMC.root", const char* folder = "Multiplicity")
+{
+ gSystem->Load("libPWG0base");
- TFile::Open(fileName);
- mult->LoadHistograms("Multiplicity");
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
+ TFile::Open(fileName);
+ mult->LoadHistograms();
mult->DrawHistograms();
+ return;
+
TH2* hist = (TH2*) gROOT->FindObject("fCorrelation3_zy");
canvas = new TCanvas("c1", "c1", 600, 500);
hist->SetStats(kFALSE);
canvas->SaveAs("Plot_Correlation.pdf");
}
-void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t eventType = 0)
-{
- gSystem->Load("libPWG0base");
-
- AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-
- TFile::Open(fileName);
- mult->LoadHistograms("Multiplicity");
-
- //mult->ApplyLaszloMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
-
- //return;
-
-
- /*mult->ApplyBayesianMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
- mult->DrawComparison("Bayesian", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());
- return;*/
-
- TStopwatch timer;
-
- timer.Start();
-
- mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.2);
- mult->ApplyMinuitFit(hist, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
- mult->DrawComparison("MinuitChi2", hist, kFALSE, kFALSE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());
-
- timer.Stop();
- timer.Print();
-
- return 0;
-
- //mult->ApplyGaussianMethod(hist, kFALSE);
-
- mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0);
- mult->ApplyNBDFit(hist, kFALSE);
- mult->DrawComparison("NBDChi2Fit", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, eventType)->ProjectionY());
-
- return mult;
-}
-
-void fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fileNameESD = "multiplicityMC_3M.root", Bool_t chi2 = kTRUE, Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE)
+void fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* folder = "Multiplicity", const char* fileNameESD = "multiplicityMC_3M.root", Bool_t chi2 = kTRUE, Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE)
{
gSystem->Load("libPWG0base");
- AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
TFile::Open(fileNameMC);
- mult->LoadHistograms("Multiplicity");
+ mult->LoadHistograms();
TFile::Open(fileNameESD);
TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
if (chi2)
{
mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
+ //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e3);
+ //mult->SetRegularizationParameters(AliMultiplicityCorrection::kLog, 1e5);
mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE); //, hist2->ProjectionY("mymchist"));
mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
}
case AliMultiplicityCorrection::kPol1: return "Pol1"; break;
case AliMultiplicityCorrection::kCurvature: return "TotalCurvature"; break;
case AliMultiplicityCorrection::kEntropy: return "Reduced cross-entropy"; break;
- case AliMultiplicityCorrection::kTest : return "Test"; break;
+ case AliMultiplicityCorrection::kLog : return "Log"; break;
}
return 0;
}
return 0;
}
-void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
+/*void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
{
gSystem->mkdir(targetDir);
Float_t chi2MC = 0;
Float_t residuals = 0;
- mult->ApplyBayesianMethod(histID, kFALSE, type, weight);
+ mult->ApplyBayesianMethod(histID, kFALSE, type, weight, 100, 0, kFALSE);
mult->DrawComparison(Form("%s/Bayesian_%02d_%d_%f", targetDir, count++, type, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY());
mult->GetComparisonResults(&chi2MC, 0, &residuals);
canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
-}
+}*/
void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
{
TFile* graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir), "RECREATE");
- for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
+ Int_t colors[3] = {1, 2, 4};
+ Int_t markers[12] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3};
+
+ for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type)
+ //for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
{
TString tmp;
tmp.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
TCanvas* canvas = new TCanvas(tmp, tmp, 800, 600);
- for (Int_t i = 1; i <= 7; i++)
+ for (Int_t i = 1; i <= 4; i++)
{
- Int_t iter = i * 20;
- TGraph* fitResultsMC = new TGraph;
- fitResultsMC->SetTitle(Form("%d iterations", iter));
- fitResultsMC->GetXaxis()->SetTitle("smoothing parameter");
- fitResultsMC->GetYaxis()->SetTitle("P_{1} (2 <= t <= 150)");
- fitResultsMC->SetName(Form("%s_MC_%d", tmp.Data(), i));
+ Int_t iterArray[4] = {5, 20, 50, 100};
+ //Int_t iter = i * 40 - 20;
+ Int_t iter = iterArray[i-1];
+
+ TGraph* fitResultsMC[3];
+ for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
+ {
+ fitResultsMC[region] = new TGraph;
+ fitResultsMC[region]->SetTitle(Form("%d iter. - reg. %d", iter, region+1));
+ fitResultsMC[region]->GetXaxis()->SetTitle("smoothing parameter #alpha");
+ fitResultsMC[region]->GetYaxis()->SetTitle(Form("P_{1} in region %d", region));
+ fitResultsMC[region]->SetName(Form("%s_MC_%d", tmp.Data(), i * AliMultiplicityCorrection::kQualityRegions + region - 2));
+ fitResultsMC[region]->SetFillColor(0);
+ fitResultsMC[region]->SetMarkerStyle(markers[(i-1) * AliMultiplicityCorrection::kQualityRegions + region]);
+ fitResultsMC[region]->SetLineColor(colors[region]);
+ }
+
TGraph* fitResultsRes = new TGraph;
fitResultsRes->SetTitle(Form("%d iterations", iter));
fitResultsRes->SetName(Form("%s_Res_%d", tmp.Data(), i));
fitResultsRes->GetXaxis()->SetTitle("smoothing parameter");
fitResultsRes->GetYaxis()->SetTitle("P_{2}");
- fitResultsMC->SetFillColor(0);
fitResultsRes->SetFillColor(0);
- fitResultsMC->SetMarkerStyle(19+i);
fitResultsRes->SetMarkerStyle(19+i);
- fitResultsRes->SetMarkerColor(kRed);
- fitResultsRes->SetLineColor(kRed);
+ fitResultsRes->SetMarkerColor(1);
+ fitResultsRes->SetLineColor(1);
for (Float_t weight = 0.0; weight < 1.01; weight += 0.2)
{
Float_t chi2MC = 0;
Float_t residuals = 0;
- mult->ApplyBayesianMethod(histID, kFALSE, type, weight, iter);
+ mult->ApplyBayesianMethod(histID, kFALSE, type, weight, iter, 0, kFALSE);
mult->DrawComparison(Form("%s/BayesianIterSmooth_%03d_%d_%d_%f", targetDir, count++, type, iter, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY());
mult->GetComparisonResults(&chi2MC, 0, &residuals);
- fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
+ for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
+ fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region));
+
fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
}
- fitResultsMC->Write();
+ graphFile->cd();
+ for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
+ fitResultsMC[region]->Write();
+
fitResultsRes->Write();
}
}
+
+ graphFile->Close();
}
void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = kTRUE)
TCanvas* canvas = new TCanvas(name, name, 800, 500);
if (type == -1)
+ {
canvas->SetLogx();
- canvas->SetLogy();
+ canvas->SetLogy();
+ }
+ canvas->SetTopMargin(0.05);
+ canvas->SetGridx();
+ canvas->SetGridy();
- TLegend* legend = new TLegend(0.6, 0.75, 0.88, 0.88);
+ TLegend* legend = new TLegend(0.8, 0.15, 0.98, 0.98);
legend->SetFillColor(0);
Int_t count = 1;
Float_t yMin = 1e20;
Float_t yMax = 0;
+ Float_t yMinRegion[3];
+ for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
+ yMinRegion[i] = 1e20;
+
TString xaxis, yaxis;
while (1)
TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
- if (!mc || !res)
+ if (!mc)
break;
xaxis = mc->GetXaxis()->GetTitle();
yaxis = mc->GetYaxis()->GetTitle();
mc->Print();
- res->Print();
- count++;
+ if (res)
+ res->Print();
- if (plotRes)
- {
- xMin = TMath::Min(TMath::Min(xMin, mc->GetXaxis()->GetXmin()), res->GetXaxis()->GetXmin());
- yMin = TMath::Min(TMath::Min(yMin, mc->GetYaxis()->GetXmin()), res->GetYaxis()->GetXmin());
+ xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin());
+ yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin());
- xMax = TMath::Max(TMath::Max(xMax, mc->GetXaxis()->GetXmax()), res->GetXaxis()->GetXmax());
- yMax = TMath::Max(TMath::Max(yMax, mc->GetYaxis()->GetXmax()), res->GetYaxis()->GetXmax());
- }
- else
+ xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax());
+ yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax());
+
+ if (plotRes && res)
{
- xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin());
- yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin());
+ xMin = TMath::Min(xMin, res->GetXaxis()->GetXmin());
+ yMin = TMath::Min(yMin, res->GetYaxis()->GetXmin());
- xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax());
- yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax());
+ xMax = TMath::Max(xMax, res->GetXaxis()->GetXmax());
+ yMax = TMath::Max(yMax, res->GetYaxis()->GetXmax());
}
+
+ for (Int_t i=0; i<mc->GetN(); ++i)
+ yMinRegion[(count-1) % 3] = TMath::Min(yMinRegion[(count-1) % 3], mc->GetY()[i]);
+
+ count++;
}
+ for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
+ Printf("Minimum for region %d is %f", i, yMinRegion[i]);
+
if (type >= 0)
{
xaxis = "smoothing parameter";
else if (type == -1)
{
xaxis = "weight parameter";
- //yaxis = "P_{3} (2 <= t <= 150)";
+ xMax *= 5;
}
- yaxis = "P_{1} (2 <= t <= 150)";
+ //yaxis = "P_{1} (2 <= t <= 150)";
printf("%f %f %f %f\n", xMin, xMax, yMin, yMax);
dummy->SetMarkerColor(0);
dummy->Draw("AP");
+ dummy->GetYaxis()->SetMoreLogLabels(1);
count = 1;
TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
- printf("%s_MC_%d %p %p\n", name.Data(), count, mc, res);
+ //printf("%s_MC_%d %p %p\n", name.Data(), count, mc, res);
- if (!mc || !res)
+ if (!mc)
break;
printf("Loaded %d sucessful.\n", count);
- if (type == -1)
+ /*if (type == -1)
{
- legend->AddEntry(mc, Form("Eq. (%d)", count+9));
+ legend->AddEntry(mc, Form("Eq. (%d)", (count-1)/3+11));
}
- else
+ else*/
legend->AddEntry(mc);
mc->Draw("SAME PC");
- if (plotRes)
+ if (plotRes && res)
{
legend->AddEntry(res);
res->Draw("SAME PC");
canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName()));
}
-void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
+void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC_2M.root", const char* fileNameESD = "multiplicityMC_1M_3.root", const char* targetDir, Int_t histID = 3)
{
gSystem->mkdir(targetDir);
mult->SetMultiplicityESD(histID, hist);
Int_t count = 0; // just to order the saved images...
+ Int_t colors[3] = {1, 2, 4};
+ Int_t markers[12] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3};
- TGraph* fitResultsMC = 0;
TGraph* fitResultsRes = 0;
TFile* graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir), "RECREATE");
- for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type)
-// for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type)
+ for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kLog; ++type)
+// for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kEntropy; type <= AliMultiplicityCorrection::kEntropy; ++type)
+// for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kPol0; ++type)
{
- fitResultsMC = new TGraph;
- fitResultsMC->SetTitle(Form("%s", GetRegName(type)));
- fitResultsMC->GetXaxis()->SetTitle("Weight Parameter");
- fitResultsMC->SetName(Form("EvaluateChi2Method_MC_%d", type));
+ TGraph* fitResultsMC[3];
+ for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
+ {
+ fitResultsMC[region] = new TGraph;
+ fitResultsMC[region]->SetTitle(Form("Eq. (%d) - reg. %d", type+10, region+1));
+ fitResultsMC[region]->GetXaxis()->SetTitle("weight parameter #alpha");
+ fitResultsMC[region]->GetYaxis()->SetTitle(Form("P_{1} in region %d", region));
+ fitResultsMC[region]->SetName(Form("EvaluateChi2Method_MC_%d", type * AliMultiplicityCorrection::kQualityRegions + region - 2));
+ fitResultsMC[region]->SetFillColor(0);
+ fitResultsMC[region]->SetMarkerStyle(markers[(type-1) * AliMultiplicityCorrection::kQualityRegions + region]);
+ fitResultsMC[region]->SetLineColor(colors[region]);
+ }
+
fitResultsRes = new TGraph;
fitResultsRes->SetTitle(Form("%s residual chi2", GetRegName(type)));
fitResultsRes->SetName(Form("EvaluateChi2Method_Res_%d", type));
fitResultsRes->GetXaxis()->SetTitle("Weight Parameter");
- fitResultsMC->SetFillColor(0);
fitResultsRes->SetFillColor(0);
- fitResultsMC->SetMarkerStyle(19+type);
fitResultsRes->SetMarkerStyle(23+type);
fitResultsRes->SetMarkerColor(kRed);
fitResultsRes->SetLineColor(kRed);
- for (Int_t i=0; i<5; ++i)
+ for (Int_t i=0; i<7; ++i)
{
Float_t weight = TMath::Power(TMath::Sqrt(10), i+6);
//Float_t weight = TMath::Power(10, i+2);
result->SetName(runName);
result->Write();
- fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
+ for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
+ fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region));
+
fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
}
graphFile->cd();
- fitResultsMC->Write();
+ for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
+ fitResultsMC[region]->Write();
fitResultsRes->Write();
}
const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_200k.root", "multiplicityMC_400k.root", "multiplicityMC_600k.root", "multiplicityMC_800k.root" };
- TGraph* fitResultsChi2 = new TGraph; fitResultsChi2->SetTitle(";Nevents;Chi2");
- TGraph* fitResultsBayes = new TGraph; fitResultsBayes->SetTitle(";Nevents;Chi2");
+ TGraph* fitResultsChi2[3];
+ TGraph* fitResultsBayes[3];
+
+ for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
+ {
+ fitResultsChi2[region] = new TGraph;
+ fitResultsChi2[region]->SetTitle(";Nevents;Chi2");
+ fitResultsChi2[region]->SetName(Form("fitResultsChi2_%d", region));
+ fitResultsChi2[region]->SetMarkerStyle(20+region);
+
+ fitResultsBayes[region] = new TGraph;
+ fitResultsBayes[region]->SetTitle(";Nevents;Chi2");
+ fitResultsBayes[region]->SetName(Form("fitResultsBayes_%d", region));
+ fitResultsBayes[region]->SetMarkerStyle(20+region);
+ fitResultsBayes[region]->SetMarkerColor(2);
+ }
+
TGraph* fitResultsChi2Limit = new TGraph; fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach");
TGraph* fitResultsBayesLimit = new TGraph; fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach");
+ TGraph* fitResultsChi2Res = new TGraph; fitResultsChi2Res->SetTitle(";Nevents;Chi2");
+ TGraph* fitResultsBayesRes = new TGraph; fitResultsBayesRes->SetTitle(";Nevents;Chi2");
- fitResultsChi2->SetName("fitResultsChi2");
- fitResultsBayes->SetName("fitResultsBayes");
fitResultsChi2Limit->SetName("fitResultsChi2Limit");
fitResultsBayesLimit->SetName("fitResultsBayesLimit");
+ fitResultsChi2Res->SetName("fitResultsChi2Res");
+ fitResultsBayesRes->SetName("fitResultsBayesRes");
TCanvas* canvas = new TCanvas("StatisticsPlot", "StatisticsPlot", 1200, 600);
canvas->Divide(5, 2);
mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
mult->DrawComparison(Form("StatisticsPlot_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
- Float_t chi2MC = 0;
Int_t chi2MCLimit = 0;
- mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
- fitResultsChi2->SetPoint(fitResultsChi2->GetN(), nEntries, chi2MC);
+ Float_t chi2Residuals = 0;
+ mult->GetComparisonResults(0, &chi2MCLimit, &chi2Residuals);
+ for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
+ {
+ fitResultsChi2[region]->SetPoint(fitResultsChi2[region]->GetN(), nEntries, mult->GetQuality(region));
+ min = TMath::Min(min, mult->GetQuality(region));
+ max = TMath::Max(max, mult->GetQuality(region));
+ }
fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), nEntries, chi2MCLimit);
- min = TMath::Min(min, chi2MC);
- max = TMath::Max(max, chi2MC);
+ fitResultsChi2Res->SetPoint(fitResultsChi2Res->GetN(), nEntries, chi2Residuals);
TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
- mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.07, 10);
+ mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0, kFALSE);
mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
- mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
- fitResultsBayes->SetPoint(fitResultsBayes->GetN(), nEntries, chi2MC);
+ mult->GetComparisonResults(0, &chi2MCLimit, &chi2Residuals);
+ for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
+ {
+ fitResultsBayes[region]->SetPoint(fitResultsBayes[region]->GetN(), nEntries, mult->GetQuality(region));
+ min = TMath::Min(min, mult->GetQuality(region));
+ max = TMath::Max(max, mult->GetQuality(region));
+ }
fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit);
+ fitResultsBayesRes->SetPoint(fitResultsBayesRes->GetN(), nEntries, chi2Residuals);
- min = TMath::Min(min, chi2MC);
- max = TMath::Max(max, chi2MC);
mc->GetXaxis()->SetRangeUser(0, 150);
chi2Result->GetXaxis()->SetRangeUser(0, 150);
canvas2->Divide(2, 1);
canvas2->cd(1);
- fitResultsChi2->SetMarkerStyle(20);
- fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
- fitResultsChi2->Draw("AP");
- fitResultsBayes->SetMarkerStyle(3);
- fitResultsBayes->SetMarkerColor(2);
- fitResultsBayes->Draw("P SAME");
+ for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
+ {
+ fitResultsChi2[region]->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
+ fitResultsChi2[region]->Draw(((region == 0) ? "AP" : "P SAME"));
+
+ fitResultsBayes[region]->Draw("P SAME");
+ }
gPad->SetLogy();
canvas2->SaveAs(Form("%s.C", canvas2->GetName()));
TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE");
- fitResultsChi2->Write();
- fitResultsBayes->Write();
+
+ for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
+ {
+ fitResultsChi2[region]->Write();
+ fitResultsBayes[region]->Write();
+ }
fitResultsChi2Limit->Write();
fitResultsBayesLimit->Write();
+ fitResultsChi2Res->Write();
+ fitResultsBayesRes->Write();
file->Close();
}
if (caseNo >= 4)
{
- 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, 250);
+ 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, 500);
func->SetParNames("scaling", "averagen", "k");
}
switch (caseNo)
{
- case 0: func = new TF1("flat", "1"); break;
+ case 0: func = new TF1("flat", "1000"); break;
case 1: func = new TF1("flat", "501-x"); break;
case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break;
case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break;
case 4: func->SetParameters(1e7, 10, 2); break;
case 5: func->SetParameters(1, 13, 7); break;
case 6: func->SetParameters(1e7, 30, 4); break;
- case 7: func->SetParameters(1e7, 30, 2); break;
+ case 7: func->SetParameters(1e7, 30, 2); break; // ***
case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break;
default: return;
TFile::Open("out.root", "RECREATE");
mult->SaveHistograms();
+ new TCanvas; mult->GetMultiplicityESD(3)->ProjectionY()->DrawCopy();
+ new TCanvas; mult->GetMultiplicityVtx(3)->ProjectionY()->DrawCopy();
+
//mult->ApplyBayesianMethod(2, kFALSE);
//mult->ApplyMinuitFit(2, kFALSE);
//mult->ApplyGaussianMethod(2, kFALSE);
proj2->Draw("SAME");
}
+void buildCorrelationMap(const char* fileName = "multiplicityMC_2M.root", Int_t corrMatrix = 3)
+{
+ gSystem->Load("libPWG0base");
+
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+ TFile::Open(fileName);
+ mult->LoadHistograms("Multiplicity");
+
+ TH3F* new = mult->GetCorrelation(corrMatrix);
+ new->Reset();
+
+ TF1* func = new TF1("func", "gaus(0)");
+
+ Int_t vtxBin = new->GetNbinsX() / 2;
+ if (vtxBin == 0)
+ vtxBin = 1;
+
+ Float_t sigma = 2;
+ for (Int_t i=1; i<=new->GetYaxis()->GetNbins(); i+=1)
+ {
+ Float_t x = new->GetYaxis()->GetBinCenter(i);
+ func->SetParameters(1, x * 0.8, sigma);
+ //func->SetParameters(1, x, sigma);
+
+ for (Int_t j=1; j<=new->GetZaxis()->GetNbins(); j+=1)
+ {
+ Float_t y = new->GetYaxis()->GetBinCenter(j);
+
+ // cut at 1 sigma
+ if (TMath::Abs(y-x*0.8) < sigma)
+ new->SetBinContent(vtxBin, i, j, func->Eval(y));
+
+ // test only bin 40 has smearing
+ //if (x != 40)
+ // new->SetBinContent(vtxBin, i, j, (i == j));
+ }
+ }
+
+ new TCanvas;
+ new->Project3D("zy")->DrawCopy("COLZ");
+
+ TFile* file = TFile::Open("out.root", "RECREATE");
+ mult->SetCorrelation(corrMatrix, new);
+ mult->SaveHistograms();
+ file->Close();
+}
+
void GetCrossSections(const char* fileName)
{
gSystem->Load("libPWG0base");
xSection15->Write();
gFile->Close();
}
+
+void BuildResponseFromTree(const char* fileName, const char* target)
+{
+ //
+ // builds several response matrices with different particle ratios (systematic study)
+ //
+
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(fileName);
+ TNtuple* fParticleSpecies = (TNtuple*) gFile->Get("fParticleSpecies");
+
+ TFile* file = TFile::Open(target, "RECREATE");
+ file->Close();
+
+ Int_t tracks = 0; // control variables
+ Int_t noLabel = 0;
+ Int_t doubleCount = 0;
+
+ for (Int_t num = 0; num < 7; num++)
+ {
+ AliMultiplicityCorrection* fMultiplicity = new AliMultiplicityCorrection(Form("Multiplicity_%d", num), Form("Multiplicity_%d", num));
+
+ Float_t ratio[4]; // pi, K, p, other
+ for (Int_t i = 0; i < 4; i++)
+ ratio[i] = 1;
+
+ switch (num)
+ {
+ case 1 : ratio[1] = 0.5; break;
+ case 2 : ratio[2] = 0.5; break;
+ case 3 : ratio[1] = 1.5; break;
+ case 4 : ratio[2] = 1.5; break;
+ case 5 : ratio[1] = 0.5; ratio[2] = 0.5; break;
+ case 6 : ratio[1] = 1.5; ratio[2] = 1.5; break;
+ }
+
+ for (Int_t i=0; i<fParticleSpecies->GetEntries(); i++)
+ {
+ fParticleSpecies->GetEvent(i);
+
+ Float_t* f = fParticleSpecies->GetArgs();
+
+ Float_t gene = 0;
+ Float_t meas = 0;
+
+ for (Int_t j = 0; j < 4; j++)
+ {
+ gene += ratio[j] * f[j+1];
+ meas += ratio[j] * f[j+1+4];
+ tracks += f[j+1+4];
+ }
+
+ // add the ones w/o label
+ tracks += f[9];
+ noLabel += f[9];
+
+ // add the double counted
+ tracks += f[10];
+ doubleCount += f[10];
+
+ // TODO fake, add the ones w/o label and the double counted to check the method
+ meas += f[9];
+ meas += f[10];
+
+ //Printf("%.f %.f %.f %.f %.f", f[5], f[6], f[7], f[8], f[9]);
+
+ fMultiplicity->FillCorrection(f[0], 0, 0, 0, gene, 0, 0, 0, 0, meas);
+ fMultiplicity->FillGenerated(f[0], kTRUE, kTRUE, 0, 0, 0, gene, 0);
+ fMultiplicity->FillMeasured(f[0], 0, 0, 0, meas);
+ }
+
+ //fMultiplicity->DrawHistograms();
+
+ TFile* file = TFile::Open(target, "UPDATE");
+ fMultiplicity->SaveHistograms();
+ file->Close();
+
+ if (num == 0)
+ Printf("%d total tracks, %d w/o label = %.2f %%, %d double counted = %.2f %% ", tracks, noLabel, 100.0 * noLabel / tracks, doubleCount, 100.0 * doubleCount / tracks);
+ }
+}
+
+void MergeModifyCrossSection(const char* output)
+{
+ const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root" };
+
+ gSystem->Load("libPWG0base");
+
+ TFile::Open(output, "RECREATE");
+ gFile->Close();
+
+ for (Int_t num=0; num<7; ++num)
+ {
+ AliMultiplicityCorrection* data[3];
+ TList list;
+
+ Float_t ratio[3];
+ switch (num)
+ {
+ case 0: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 1.0; break;
+ case 1: ratio[0] = 1.0; ratio[1] = 1.5; ratio[2] = 1.0; break;
+ case 2: ratio[0] = 1.0; ratio[1] = 0.5; ratio[2] = 1.0; break;
+ case 3: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 1.5; break;
+ case 4: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 0.5; break;
+ case 5: ratio[0] = 1.0; ratio[1] = 1.5; ratio[2] = 1.5; break;
+ case 6: ratio[0] = 1.0; ratio[1] = 0.5; ratio[2] = 0.5; break;
+ default: return;
+ }
+
+ for (Int_t i=0; i<3; ++i)
+ {
+ TString name;
+ name.Form("Multiplicity_%d", num);
+ if (i > 0)
+ name.Form("Multiplicity_%d_%d", num, i);
+
+ TFile::Open(files[i]);
+ data[i] = new AliMultiplicityCorrection(name, name);
+ data[i]->LoadHistograms("Multiplicity");
+
+ // modify x-section
+ for (Int_t j=0; j<AliMultiplicityCorrection::kMCHists; j++)
+ {
+ data[i]->GetMultiplicityVtx(j)->Scale(ratio[i]);
+ data[i]->GetMultiplicityMB(j)->Scale(ratio[i]);
+ data[i]->GetMultiplicityINEL(j)->Scale(ratio[i]);
+ }
+
+ for (Int_t j=0; j<AliMultiplicityCorrection::kESDHists; j++)
+ data[i]->GetMultiplicityESD(j)->Scale(ratio[i]);
+
+ for (Int_t j=0; j<AliMultiplicityCorrection::kCorrHists; j++)
+ data[i]->GetCorrelation(j)->Scale(ratio[i]);
+
+ if (i > 0)
+ list.Add(data[i]);
+ }
+
+ printf("Case %d, %s: Entries in response matrix 3: ND: %.2f SD: %.2f DD: %.2f", num, data[0]->GetName(), data[0]->GetCorrelation(3)->Integral(), data[1]->GetCorrelation(3)->Integral(), data[2]->GetCorrelation(3)->Integral());
+
+ data[0]->Merge(&list);
+
+ Printf(" Total: %.2f", data[0]->GetCorrelation(3)->Integral());
+
+ TFile::Open(output, "UPDATE");
+ data[0]->SaveHistograms();
+ gFile->Close();
+
+ list.Clear();
+
+ for (Int_t i=0; i<3; ++i)
+ delete data[i];
+ }
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff