From 0b4bfd988bb8310a2839ddb2b0001abdb2c7f24f Mon Sep 17 00:00:00 2001 From: jgrosseo Date: Thu, 19 Jul 2007 07:20:28 +0000 Subject: [PATCH] added study for bayesian uncertainty optimized unfolding methods in speed by factor 20 --- PWG0/dNdEta/AliMultiplicityCorrection.cxx | 992 +++++++++-- PWG0/dNdEta/AliMultiplicityCorrection.h | 35 +- PWG0/dNdEta/plotsMultiplicity.C | 1945 +++++++++++++++++++++ PWG0/dNdEta/runMultiplicitySelector.C | 515 ++++-- 4 files changed, 3172 insertions(+), 315 deletions(-) create mode 100644 PWG0/dNdEta/plotsMultiplicity.C diff --git a/PWG0/dNdEta/AliMultiplicityCorrection.cxx b/PWG0/dNdEta/AliMultiplicityCorrection.cxx index 45bc2413466..e25996e265f 100644 --- a/PWG0/dNdEta/AliMultiplicityCorrection.cxx +++ b/PWG0/dNdEta/AliMultiplicityCorrection.cxx @@ -36,6 +36,9 @@ #include #include #include +#include +#include +#include ClassImp(AliMultiplicityCorrection) @@ -53,9 +56,50 @@ AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fRegula 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 @@ -80,7 +124,7 @@ AliMultiplicityCorrection::AliMultiplicityCorrection() : //____________________________________________________________________ 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 @@ -412,7 +456,7 @@ Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params) } //____________________________________________________________________ -Double_t AliMultiplicityCorrection::RegularizationTest(TVectorD& params) +Double_t AliMultiplicityCorrection::RegularizationLog(TVectorD& params) { // homogenity term for minuit fitting // pure function of the parameters @@ -420,7 +464,7 @@ Double_t AliMultiplicityCorrection::RegularizationTest(TVectorD& params) Double_t chi2 = 0; - Float_t der[fgMaxParams]; + /*Float_t der[fgMaxParams]; for (Int_t i=0; i 0 && (*fEntropyAPriori)[i] > 0) - chi2 += tmp * log(tmp / (*fEntropyAPriori)[i]); + chi2 += tmp * TMath::Log(tmp / (*fEntropyAPriori)[i]); } return 10.0 + chi2; @@ -523,7 +586,7 @@ void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& c // // d - TVectorD paramsVector(fgMaxParams); + static TVectorD paramsVector(fgMaxParams); for (Int_t i=0; i 1e10) @@ -554,11 +617,18 @@ void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& c 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; iReset(); + 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 @@ -612,15 +686,27 @@ void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullP 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())); @@ -649,29 +735,57 @@ void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullP 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; } //____________________________________________________________________ @@ -696,7 +810,9 @@ Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPha 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); @@ -741,9 +857,15 @@ Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPha // 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; iRebin(2); proj->Scale(1.0 / proj->Integral()); - Double_t results[fgMaxParams]; + Double_t results[fgMaxParams+1]; for (Int_t i=0; iGetBinContent(i+1); @@ -790,6 +912,7 @@ Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPha 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; @@ -856,6 +979,8 @@ void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSp 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); @@ -1028,6 +1153,21 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang 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; iGetNbinsX(); ++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 @@ -1053,7 +1193,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang 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; @@ -1079,8 +1219,10 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang 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); @@ -1140,8 +1282,8 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang legend->SetFillColor(0); legend->Draw(); - TH1* diffMCUnfolded = dynamic_cast (proj->Clone("diffMCUnfolded")); - diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1); + //TH1* diffMCUnfolded = dynamic_cast (proj->Clone("diffMCUnfolded")); + //diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1); /*Float_t chi2 = 0; Float_t chi = 0; @@ -1166,7 +1308,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang } }*/ - chi2 = 0; + /*chi2 = 0; Float_t chi = 0; Int_t limit = 0; for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i) @@ -1201,13 +1343,13 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang 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"); @@ -1216,8 +1358,8 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang 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"); @@ -1241,6 +1383,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang ratio->Draw("HIST"); Double_t ratioChi2 = 0; + fRatioAverage = 0; fLastChi2MCLimit = 0; for (Int_t i=2; i<=150; ++i) { @@ -1251,21 +1394,270 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang 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; iGetBinContent(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; iSetBinContent(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; iSetBinContent(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 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; iSetBinContent(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 (proj->Clone("diffMCUnfolded2")); + diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1); + + Int_t ndfQual[kQualityRegions]; + for (Int_t region=0; regionGetNbinsX(), 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; regionGetXaxis()->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 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 (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; nClone("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; nGetNbinsX(); 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; nGetNbinsX(); 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; nGetNbinsX(); 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; nGetNbinsX(); 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; nGetNbinsX(); ++i) + fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i)); - // smooth efficiency - //TH1* tmp = (TH1*) fCurrentEfficiency->Clone("eff_clone"); - //for (Int_t i=2; iGetNbinsX(); ++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) @@ -1345,80 +1942,133 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful } } - // 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; nGetNbinsX(); 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; iGetNbinsX(); ++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; mGetBinContent(m+1) / measuredIntegral; + + for (Int_t t=0; tGetBinContent(t+1, m+1); + inverseResponse[t][m] = 0; + } + } + + for (Int_t t=0; tGetBinContent(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; iGetBinContent(i+1) / inputDistIntegral; + } + + //TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5); // unfold... for (Int_t i=0; iGetNbinsY(); m++) + for (Int_t m=0; mGetNbinsX(); t++) - norm += hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t); - for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++) + for (Int_t t = kStartBin; t 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; tGetNbinsX(); t++) + for (Int_t t = kStartBin; tGetNbinsY(); m++) - value += fCurrentESD->GetBinContent(m) * hInverseResponseBayes->GetBinContent(t,m); + for (Int_t m=0; mGetBinContent(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; tGetNbinsX(); t++) + for (Int_t t=kStartBin; tGetBinContent(t-1) - + hTemp->GetBinContent(t) - + hTemp->GetBinContent(t+1) - ) / 3.; + Float_t newValue = 0; + // 0 bin excluded from smoothing + if (t > kStartBin+1 && tSetBinContent(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; tSetBinContent(t+1, result[t]); - return; + return resultHist; // ******** // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995 @@ -1646,8 +2287,7 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful } //____________________________________________________________________ -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 @@ -1682,6 +2322,8 @@ void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullP 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()); @@ -1802,6 +2444,8 @@ void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t ful Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE); + //normalize ESD + fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); NormalizeToBinWidth((TH2*) fCurrentCorrelation); @@ -1956,8 +2600,8 @@ void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id) 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; diff --git a/PWG0/dNdEta/AliMultiplicityCorrection.h b/PWG0/dNdEta/AliMultiplicityCorrection.h index a04f580d2f4..3da2d88fd3b 100644 --- a/PWG0/dNdEta/AliMultiplicityCorrection.h +++ b/PWG0/dNdEta/AliMultiplicityCorrection.h @@ -24,7 +24,8 @@ class TCollection; 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); @@ -37,7 +38,7 @@ class AliMultiplicityCorrection : public TNamed { 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); @@ -47,7 +48,8 @@ class AliMultiplicityCorrection : public TNamed { 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); @@ -61,9 +63,9 @@ class AliMultiplicityCorrection : public TNamed { 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; } @@ -74,11 +76,16 @@ class AliMultiplicityCorrection : public TNamed { 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 @@ -86,14 +93,15 @@ class AliMultiplicityCorrection : public TNamed { 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 @@ -120,6 +128,11 @@ class AliMultiplicityCorrection : public TNamed { 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&); diff --git a/PWG0/dNdEta/plotsMultiplicity.C b/PWG0/dNdEta/plotsMultiplicity.C new file mode 100644 index 00000000000..6d418e76a73 --- /dev/null +++ b/PWG0/dNdEta/plotsMultiplicity.C @@ -0,0 +1,1945 @@ +// +// 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; nScale(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; nScale(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; nScale(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; nScale(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; iLoadHistograms(); + + 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; iLoadHistograms(); + + 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; xGetNbinsX(); 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; paramLoadHistograms("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"); +} diff --git a/PWG0/dNdEta/runMultiplicitySelector.C b/PWG0/dNdEta/runMultiplicitySelector.C index 28364d69833..0a26fcee07d 100644 --- a/PWG0/dNdEta/runMultiplicitySelector.C +++ b/PWG0/dNdEta/runMultiplicitySelector.C @@ -10,7 +10,7 @@ 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"); @@ -58,17 +58,24 @@ void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_ } } -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); @@ -81,54 +88,14 @@ void draw(const char* fileName = "multiplicityMC.root") 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)); @@ -157,6 +124,8 @@ void fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fil 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")); } @@ -223,7 +192,7 @@ const char* GetRegName(Int_t type) 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; } @@ -239,7 +208,7 @@ const char* GetEventTypeName(Int_t type) 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); @@ -286,7 +255,7 @@ void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", cons 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); @@ -399,7 +368,7 @@ void EvaluateBayesianMethodIterations(const char* fileNameMC = "multiplicityMC.r 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) { @@ -420,51 +389,71 @@ void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multipl 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; regionSetTitle(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; regionSetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region)); + fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); } - fitResultsMC->Write(); + graphFile->cd(); + for (Int_t region=0; regionWrite(); + fitResultsRes->Write(); } } + + graphFile->Close(); } void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = kTRUE) @@ -486,10 +475,15 @@ void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = 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; @@ -500,6 +494,10 @@ void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = Float_t yMin = 1e20; Float_t yMax = 0; + Float_t yMinRegion[3]; + for (Int_t i=0; iGet(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; iGetN(); ++i) + yMinRegion[(count-1) % 3] = TMath::Min(yMinRegion[(count-1) % 3], mc->GetY()[i]); + + count++; } + for (Int_t i=0; i= 0) { xaxis = "smoothing parameter"; @@ -543,9 +547,9 @@ void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = 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); @@ -556,6 +560,7 @@ void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = dummy->SetMarkerColor(0); dummy->Draw("AP"); + dummy->GetYaxis()->SetMoreLogLabels(1); count = 1; @@ -564,23 +569,23 @@ void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = 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"); @@ -595,7 +600,7 @@ void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = 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); @@ -613,32 +618,41 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch 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; regionSetTitle(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); @@ -666,12 +680,15 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch result->SetName(runName); result->Write(); - fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC); + for (Int_t region=0; regionSetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region)); + fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); } graphFile->cd(); - fitResultsMC->Write(); + for (Int_t region=0; regionWrite(); fitResultsRes->Write(); } @@ -783,15 +800,32 @@ void StatisticsPlot(const char* fileNameMC = "multiplicityMC_2M.root", Int_t his 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; regionSetTitle(";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); @@ -816,25 +850,33 @@ void StatisticsPlot(const char* fileNameMC = "multiplicityMC_2M.root", Int_t his 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; regionSetPoint(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; regionSetPoint(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); @@ -885,13 +927,14 @@ void StatisticsPlot(const char* fileNameMC = "multiplicityMC_2M.root", Int_t his 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; regionGetYaxis()->SetRangeUser(0.5 * min, 1.5 * max); + fitResultsChi2[region]->Draw(((region == 0) ? "AP" : "P SAME")); + + fitResultsBayes[region]->Draw("P SAME"); + } gPad->SetLogy(); @@ -908,10 +951,16 @@ void StatisticsPlot(const char* fileNameMC = "multiplicityMC_2M.root", Int_t his canvas2->SaveAs(Form("%s.C", canvas2->GetName())); TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE"); - fitResultsChi2->Write(); - fitResultsBayes->Write(); + + for (Int_t region=0; regionWrite(); + fitResultsBayes[region]->Write(); + } fitResultsChi2Limit->Write(); fitResultsBayesLimit->Write(); + fitResultsChi2Res->Write(); + fitResultsBayesRes->Write(); file->Close(); } @@ -1363,20 +1412,20 @@ void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root") 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; @@ -1390,6 +1439,9 @@ void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root") 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); @@ -1633,6 +1685,54 @@ void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t co 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"); @@ -1655,3 +1755,158 @@ void GetCrossSections(const char* fileName) 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; iGetEntries(); 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; jGetMultiplicityVtx(j)->Scale(ratio[i]); + data[i]->GetMultiplicityMB(j)->Scale(ratio[i]); + data[i]->GetMultiplicityINEL(j)->Scale(ratio[i]); + } + + for (Int_t j=0; jGetMultiplicityESD(j)->Scale(ratio[i]); + + for (Int_t j=0; jGetCorrelation(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]; + } +} -- 2.39.3