From dd70110949cdcdc8c7c1cba590a9a2ea7e6de638 Mon Sep 17 00:00:00 2001 From: jgrosseo Date: Wed, 11 Apr 2007 13:18:30 +0000 Subject: [PATCH] o) Performance studies added (mainly in runMultiplicitySelector) that evaluate the free parameters of the chi2 and bayesian method. o) Cleaned up AliMultiplicityCorrection o) Implemented selector to read the spectrum from the ESD --- PWG0/dNdEta/AliMultiplicityCorrection.cxx | 575 +++++++------ PWG0/dNdEta/AliMultiplicityCorrection.h | 20 +- PWG0/dNdEta/AliMultiplicityESDSelector.cxx | 168 +++- PWG0/dNdEta/AliMultiplicityESDSelector.h | 19 +- PWG0/dNdEta/AliMultiplicityMCSelector.cxx | 7 +- PWG0/dNdEta/AliMultiplicityMCSelector.h | 2 - PWG0/dNdEta/runMultiplicitySelector.C | 915 ++++++++++++++++++++- 7 files changed, 1360 insertions(+), 346 deletions(-) diff --git a/PWG0/dNdEta/AliMultiplicityCorrection.cxx b/PWG0/dNdEta/AliMultiplicityCorrection.cxx index 3f28db6f288..bdb827f0d06 100644 --- a/PWG0/dNdEta/AliMultiplicityCorrection.cxx +++ b/PWG0/dNdEta/AliMultiplicityCorrection.cxx @@ -46,10 +46,11 @@ TMatrixF* AliMultiplicityCorrection::fCorrelationCovarianceMatrix = 0; TVectorF* AliMultiplicityCorrection::fCurrentESDVector = 0; AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fRegularizationType = AliMultiplicityCorrection::kPol1; Float_t AliMultiplicityCorrection::fRegularizationWeight = 1e4; +TF1* AliMultiplicityCorrection::fNBD = 0; //____________________________________________________________________ AliMultiplicityCorrection::AliMultiplicityCorrection() : - TNamed(), fLastChi2MC(0), fLastChi2Residuals(0) + TNamed(), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0) { // // default constructor @@ -74,7 +75,7 @@ AliMultiplicityCorrection::AliMultiplicityCorrection() : //____________________________________________________________________ AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) : - TNamed(name, title), fLastChi2MC(0), fLastChi2Residuals(0) + TNamed(name, title), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0) { // // named constructor @@ -361,17 +362,14 @@ Double_t AliMultiplicityCorrection::RegularizationPol0(Double_t *params) for (Int_t i=1; iGetBinWidth(i+1); - Double_t left = params[i-1]; // / fCurrentESD->GetBinWidth(i); + Double_t right = params[i]; + Double_t left = params[i-1]; - Double_t diff = (right - left) / right; + Double_t diff = (right - left); chi2 += diff * diff; } - return chi2; + return chi2 * 100; } //____________________________________________________________________ @@ -385,25 +383,48 @@ Double_t AliMultiplicityCorrection::RegularizationPol1(Double_t *params) for (Int_t i=2; iGetBinWidth(i+1); - Double_t middle = params[i-1]; // / fCurrentESD->GetBinWidth(i); - Double_t left = params[i-2]; // / fCurrentESD->GetBinWidth(i-1); + Double_t right = params[i]; + Double_t middle = params[i-1]; + Double_t left = params[i-2]; + + Double_t der1 = (right - middle); + Double_t der2 = (middle - left); + + Double_t diff = (der1 - der2) / middle; + + chi2 += diff * diff; + } - Double_t der1 = (right - middle); // / fCurrentESD->GetBinWidth(i); - Double_t der2 = (middle - left); // / fCurrentESD->GetBinWidth(i-1); + return chi2; +} - Double_t diff = (der1 - der2) / middle; /// fCurrentESD->GetBinWidth(i); +//____________________________________________________________________ +Double_t AliMultiplicityCorrection::RegularizationTest(Double_t *params) +{ + // homogenity term for minuit fitting + // pure function of the parameters + // prefers linear function (pol1) - // TODO give additonal weight to big bins - chi2 += diff * diff; // * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1); + Double_t chi2 = 0; - //printf("%d --> %f\n", i, diff); + Float_t der[fgMaxParams]; + + for (Int_t i=0; iGetBinWidth(i+1); - Double_t middle = params[i-1]; // / fCurrentESD->GetBinWidth(i); - Double_t left = params[i-2]; // / fCurrentESD->GetBinWidth(i-1); + Double_t der1 = (right - middle); + Double_t der2 = (middle - left); - Double_t der1 = (right - middle); // / fCurrentESD->GetBinWidth(i); - Double_t der2 = (middle - left); // / fCurrentESD->GetBinWidth(i-1); - - Double_t secDer = (der1 - der2); // / fCurrentESD->GetBinWidth(i); + Double_t secDer = (der1 - der2); // square and weight with the bin width - chi2 += secDer * secDer; // * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1); - - //printf("%d --> %f\n", i, secDer); + chi2 += secDer * secDer; } - return chi2; // 10 + return chi2 * 100; } //____________________________________________________________________ @@ -447,31 +463,41 @@ Double_t AliMultiplicityCorrection::RegularizationEntropy(Double_t *params) // calculates entropy, from // The method of reduced cross-entropy (M. Schmelling 1993) - //static Int_t callCount = 0; - Double_t paramSum = 0; for (Int_t i=0; iGetBinWidth(i+1); + paramSum += params[i]; Double_t chi2 = 0; for (Int_t i=0; iGetBinWidth(i+1); + Double_t tmp = params[i] / paramSum; if (tmp > 0) - { - chi2 += tmp * log(tmp); /* * fCurrentESD->GetBinWidth(i+1); - // /\ - // ------------------------------------ - // TODO WARNING the absolute fitting values seem to depend on that value | - // this should be impossible... - //if ((callCount % 100) == 0) - // printf("%f => %f\n", params[i], tmp * log(tmp)); */ - } + chi2 += tmp * log(tmp); } - //if ((callCount++ % 100) == 0) - // printf("\n"); - return chi2; // 1000 + return chi2; +} + +//____________________________________________________________________ +void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3) +{ + // + // fit function for minuit + // does: nbd + // 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, 50); + // func->SetParNames("scaling", "averagen", "k"); + // + + fNBD->SetParameters(params[0], params[1], params[2]); + + Double_t params2[fgMaxParams]; + + for (Int_t i=0; iEval(i); + + MinuitFitFunction(unused1, unused2, chi2, params2, unused3); + + printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2); } //____________________________________________________________________ @@ -505,68 +531,16 @@ void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& c // (Ad - m) W (Ad - m) chi2FromFit = paramsVector * copy; - /*printf("chi2FromFit = %f\n", chi2FromFit); - chi2FromFit = 0; - - // loop over multiplicity (x axis of fMultiplicityESD) - for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i) - { - if (i > fCurrentCorrelation->GetNbinsY() || i > fgMaxParams) - break; - - Double_t sum = 0; - //Double_t error = 0; - // loop over generated particles (x axis of fCorrelation) that resulted in reconstruction of i tracks - for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsX(); ++j) - { - if (j > fgMaxParams) - break; - - sum += fCurrentCorrelation->GetBinContent(j, i) * params[j-1]; - - //if (params[j-1] > 0) - // error += fCurrentCorrelation->GetBinError(j, i) * fCurrentCorrelation->GetBinError(j, i) * params[j-1]; - } - - //printf("%f\n", sum); - - Double_t diff = fCurrentESD->GetBinContent(i) - sum; - - // include error - if (fCurrentESD->GetBinError(i) > 0) - diff /= fCurrentESD->GetBinError(i); - - //if (fCurrentESD->GetBinContent(i) > 0) - // diff /= fCurrentESD->GetBinContent(i); - //else - // diff /= fCurrentESD->Integral(); - - // weight with bin width - //diff *= fCurrentESD->GetBinWidth(i); - - //error = TMath::Sqrt(error) + fCurrentESD->GetBinError(i); - //if (error <= 0) - // error = 1; - - //Double_t tmp = diff / error; - //chi2 += tmp * tmp; - chi2FromFit += diff * diff; - - //printf("\nExpected sum = %f; Diff for bin %d is %f\n**********************************\n", fCurrentESD->GetBinContent(i), i, diff); - //printf("Diff for bin %d is %f\n", i, diff); - } - - printf("chi2FromFit = %f\n", chi2FromFit);*/ - Double_t penaltyVal = 0; switch (fRegularizationType) { - case kNone: break; - case kPol0: penaltyVal = RegularizationPol0(params); break; - case kPol1: penaltyVal = RegularizationPol1(params); break; - case kCurvature: penaltyVal = RegularizationTotalCurvature(params); break; - case kEntropy: penaltyVal = RegularizationEntropy(params); break; + case kNone: break; + case kPol0: penaltyVal = RegularizationPol0(params); break; + case kPol1: penaltyVal = RegularizationPol1(params); break; + case kCurvature: penaltyVal = RegularizationTotalCurvature(params); break; + case kEntropy: penaltyVal = RegularizationEntropy(params); break; + case kTest: penaltyVal = RegularizationTest(params); break; } penaltyVal *= fRegularizationWeight; @@ -617,7 +591,7 @@ void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullP } //____________________________________________________________________ -void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, Bool_t check) +void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* inputDist) { // // correct spectrum using minuit chi2 method @@ -629,7 +603,7 @@ void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhas Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); - SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx); + SetupCurrentHists(inputRange, fullPhaseSpace, eventType); fCorrelationMatrix = new TMatrixF(fgMaxParams, fgMaxParams); fCorrelationCovarianceMatrix = new TMatrixF(fgMaxParams, fgMaxParams); @@ -671,45 +645,51 @@ void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhas fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i-1, j)); }*/ - /*TCanvas* canvas1 = new TCanvas("ApplyMinuitFit", "ApplyMinuitFit", 800, 400); - canvas1->Divide(2, 1); - canvas1->cd(1); - fCurrentESD->DrawCopy(); - canvas1->cd(2); - fCurrentCorrelation->DrawCopy("COLZ");*/ - // Initialize TMinuit via generic fitter interface TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams); minuit->SetFCN(MinuitFitFunction); static Double_t results[fgMaxParams]; - //printf("%x\n", (void*) results); - TH1* proj = fMultiplicityVtx[mcTarget]->ProjectionY(); + if (inputDist) + { + printf("Using different starting conditions...\n"); + new TCanvas; + inputDist->DrawCopy(); + } + else + inputDist = fCurrentESD; + + // normalize ESD + fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); + inputDist->Scale(1.0 / inputDist->Integral()); + + Float_t minStartValue = 1e-3; + + TH1* proj = fMultiplicityVtx[mcTarget]->ProjectionY("check-proj"); + proj->Scale(1.0 / proj->Integral()); for (Int_t i=0; iGetBinContent(i+1); - if (results[i] < 1e-2) - results[i] = 1e-2; + results[i] = inputDist->GetBinContent(i+1); if (check) results[i] = proj->GetBinContent(i+1); - minuit->SetParameter(i, Form("param%d", i), results[i], results[i] / 10, 0.001, fCurrentESD->GetMaximum() * 10); + + printf("%f %f %f\n", inputDist->GetBinContent(i+1), TMath::Sqrt(inputDist->GetBinContent(i+1)), inputDist->GetBinError(i+1)); (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1); if (fCurrentESD->GetBinError(i+1) > 0) - (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1)); + (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / fCurrentESD->GetBinError(i+1); - //minuit->SetParameter(i, Form("param%d", i), fCurrentESD->GetBinWidth(i+1), 0.01, 0.001, 50); + minuit->SetParameter(i, Form("param%d", i), ((results[i] > minStartValue) ? results[i] : minStartValue), 0.1, 0, 1); } - minuit->SetParameter(0, "param0", results[1], ((results[1] > 1) ? (results[1] / 10) : 0), 0.001, fCurrentESD->GetMaximum() * 10); - //minuit->SetParameter(0, "param0", 1, 0, 1, 1); + // 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, 1); Int_t dummy; Double_t chi2; MinuitFitFunction(dummy, 0, chi2, results, 0); - printf("Chi2 of right solution is = %f\n", chi2); + printf("Chi2 of initial parameters is = %f\n", chi2); if (check) return; @@ -718,16 +698,95 @@ void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhas arglist[0] = 0; minuit->ExecuteCommand("SET PRINT", arglist, 1); minuit->ExecuteCommand("MIGRAD", arglist, 1); + //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!"); //minuit->ExecuteCommand("MINOS", arglist, 0); for (Int_t i=0; iGetParameter(i); - fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, results[i]); - fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0); + if (fCurrentEfficiency->GetBinContent(i+1) > 0) + { + fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, minuit->GetParameter(i) / fCurrentEfficiency->GetBinContent(i+1)); + fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, minuit->GetParError(i) / fCurrentEfficiency->GetBinContent(i+1)); + } + else + { + fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, 0); + fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0); + } + } +} + +//____________________________________________________________________ +void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace) +{ + // + // correct spectrum using minuit chi2 method applying a NBD fit + // + + Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); + + SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx); + + fCorrelationMatrix = new TMatrixF(fgMaxParams, fgMaxParams); + fCorrelationCovarianceMatrix = new TMatrixF(fgMaxParams, fgMaxParams); + fCurrentESDVector = new TVectorF(fgMaxParams); + + // normalize correction for given nPart + for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) + { + Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); + if (sum <= 0) + continue; + for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) + { + // npart sum to 1 + fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); + fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); + + if (i <= fgMaxParams && j <= fgMaxParams) + (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j); + } + } + + for (Int_t i=0; iGetBinContent(i+1); + if (fCurrentESD->GetBinError(i+1) > 0) + (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1)); + } + + // Create NBD function + if (!fNBD) + { + fNBD = 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); + fNBD->SetParNames("scaling", "averagen", "k"); } - //printf("Penalty is %f\n", RegularizationTotalCurvature(results)); + // Initialize TMinuit via generic fitter interface + TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3); + + minuit->SetFCN(MinuitNBD); + minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000); + minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000); + minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000); + + Double_t arglist[100]; + arglist[0] = 0; + minuit->ExecuteCommand("SET PRINT", arglist, 1); + minuit->ExecuteCommand("MIGRAD", arglist, 1); + //minuit->ExecuteCommand("MINOS", arglist, 0); + + fNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2)); + + for (Int_t i=0; iEval(i)); + if (fNBD->Eval(i) > 0) + { + fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fNBD->Eval(i)); + fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0); + } + } } //____________________________________________________________________ @@ -760,20 +819,6 @@ void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist) } } -//____________________________________________________________________ -void AliMultiplicityCorrection::ApplyMinuitFitAll() -{ - // - // fit all eta ranges - // - - for (Int_t i=0; iClone("tmpESDEfficiencyRecorrected"); + tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency); + TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId, !normalizeESD); TH1* proj2 = convoluted->ProjectionY("proj2", -1, -1, "e"); NormalizeToBinWidth(proj2); - TH1* residual = convoluted->ProjectionY("residual", -1, -1, "e"); - residual->SetTitle("Residuals"); + + TH1* residual = (TH1*) proj2->Clone("residual"); + residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / measured"); residual->Add(fCurrentESD, -1); residual->Divide(residual, fCurrentESD, 1, 1, "B"); @@ -846,6 +894,16 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang residual->SetBinError(i, 0); } + // normalize mc to 1 + //mcHist->Scale(1.0 / mcHist->Integral(1, 200)); + + // normalize unfolded esd to 1 + //fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral(1, 200)); + + // normalize unfolded result to mc hist + fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral()); + fMultiplicityESDCorrected[esdCorrId]->Scale(mcHist->Integral(1, 200)); + TCanvas* canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 800); canvas1->Divide(2, 2); @@ -856,7 +914,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang if (normalizeESD) NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]); - proj->GetXaxis()->SetRangeUser(0, 150); + proj->GetXaxis()->SetRangeUser(0, 200); proj->DrawCopy("HIST"); gPad->SetLogy(); @@ -865,23 +923,29 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST"); Float_t chi2 = 0; - for (Int_t i=1; i<=100; ++i) + fLastChi2MCLimit = 0; + for (Int_t i=2; i<=200; ++i) { - if (fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i) != 0) + if (proj->GetBinContent(i) != 0) { - Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i); + Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i); chi2 += value * value; + + if (chi2 < 0.1) + fLastChi2MCLimit = i; + + if (i == 100) + fLastChi2MC = chi2; } } - printf("Difference (from MC) is %f for bin 1-100\n", chi2); - fLastChi2MC = chi2; + printf("Difference (from MC) is %f for bin 2-100. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit); canvas1->cd(2); NormalizeToBinWidth(fCurrentESD); gPad->SetLogy(); - fCurrentESD->GetXaxis()->SetRangeUser(0, 150); + fCurrentESD->GetXaxis()->SetRangeUser(0, 200); //fCurrentESD->SetLineColor(2); fCurrentESD->DrawCopy("HIST"); @@ -891,7 +955,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang proj2->DrawCopy("HIST SAME"); chi2 = 0; - for (Int_t i=1; i<=100; ++i) + for (Int_t i=2; i<=100; ++i) { if (fCurrentESD->GetBinContent(i) != 0) { @@ -900,15 +964,15 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang } } - printf("Difference (Residuals) is %f for bin 1-100\n", chi2); + printf("Difference (Residuals) is %f for bin 2-100\n", chi2); fLastChi2Residuals = chi2; canvas1->cd(3); TH1* clone = dynamic_cast (proj->Clone("clone")); clone->Divide(fMultiplicityESDCorrected[esdCorrId]); - clone->SetTitle("Ratio;Npart;MC/ESD"); + clone->SetTitle("Ratio;Npart;MC/Unfolded Measured"); clone->GetYaxis()->SetRangeUser(0.8, 1.2); - clone->GetXaxis()->SetRangeUser(0, 150); + clone->GetXaxis()->SetRangeUser(0, 200); clone->DrawCopy("HIST"); /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85); @@ -918,23 +982,26 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang legend->Draw();*/ canvas1->cd(4); - residual->GetYaxis()->SetRangeUser(-0.2, 0.2); - residual->GetXaxis()->SetRangeUser(0, 150); + residual->GetXaxis()->SetRangeUser(0, 200); residual->DrawCopy(); canvas1->SaveAs(Form("%s.gif", canvas1->GetName())); } //____________________________________________________________________ -void AliMultiplicityCorrection::GetComparisonResults(Float_t& mc, Float_t& residuals) +void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals) { // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD // These values are computed during DrawComparison, thus this function picks up the // last calculation - mc = fLastChi2MC; - residuals = fLastChi2Residuals; + if (mc) + *mc = fLastChi2MC; + if (mcLimit) + *mcLimit = fLastChi2MCLimit; + if (residuals) + *residuals = fLastChi2Residuals; } @@ -956,7 +1023,7 @@ TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType) } //____________________________________________________________________ -void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar) +void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* inputDist) { // // correct spectrum using bayesian method @@ -966,52 +1033,59 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful SetupCurrentHists(inputRange, fullPhaseSpace, eventType); - // TODO should be taken from correlation map - TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX()); - - //new TCanvas; - //sumHist->Draw(); + // normalize ESD + fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); // normalize correction for given nPart for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i) { - //Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY()); - Double_t sum = sumHist->GetBinContent(i); - if (sum <= 0) - continue; + // 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) { - // npart sum to 1 - fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); - fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); + 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); + } } } - //new TCanvas; - //fCurrentCorrelation->Draw("COLZ"); - //return; + // smooth input spectrum + /* + TH1* esdClone = (TH1*) fCurrentESD->Clone("esdClone"); - // normalize correction for given nTrack - /*for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j) - { - Double_t sum = fCurrentCorrelation->Integral(1, fCurrentCorrelation->GetNbinsX(), j, j); - if (sum <= 0) - continue; - for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsY(); ++i) - { - // ntrack sum to 1 - fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum); - fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum); - } - }*/ + 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); - // FAKE - fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff"); - //new TCanvas; - //fCurrentEfficiency->Draw(); + delete esdClone; + + // rescale to 1 + fCurrentESD->Scale(1.0 / fCurrentESD->Integral()); + */ // pick prior distribution - TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior"); + TH1* hPrior = 0; + if (inputDist) + { + printf("Using different starting conditions...\n"); + hPrior = (TH1*)inputDist->Clone("prior"); + } + 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); @@ -1027,9 +1101,10 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful TH2F* hInverseResponseBayes = (TH2F*)hResponse->Clone("pji"); hInverseResponseBayes->Reset(); + TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5); + // unfold... - Int_t iterations = 20; - for (Int_t i=0; iGetNbinsY(); m++) value += fCurrentESD->GetBinContent(m) * hInverseResponseBayes->GetBinContent(t,m); - /*if (eventType == kTrVtx) - { - hTemp->SetBinContent(t, value); - } - else*/ - { - if (fCurrentEfficiency->GetBinContent(t) > 0) - hTemp->SetBinContent(t, value / fCurrentEfficiency->GetBinContent(t)); - else - hTemp->SetBinContent(t, 0); - } + if (fCurrentEfficiency->GetBinContent(t) > 0) + hTemp->SetBinContent(t, value / fCurrentEfficiency->GetBinContent(t)); + else + hTemp->SetBinContent(t, 0); } // this is the last guess, fill before (!) smoothing @@ -1083,10 +1151,10 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful for (Int_t t=2; tGetNbinsX(); t++) { - Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1) - + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t) - + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.; - average *= hTrueSmoothed->GetBinWidth(t); + Float_t average = (hTemp->GetBinContent(t-1) + + hTemp->GetBinContent(t) + + hTemp->GetBinContent(t+1) + ) / 3.; // weight the average with the regularization parameter hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average); @@ -1111,19 +1179,29 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful //printf("Chi2 of %d iteration = %.10f\n", i, chi2); + convergence->Fill(i+1, chi2); + //if (i % 5 == 0) // DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType); delete hTrueSmoothed; } // end of iterations + //new TCanvas; + //convergence->DrawCopy(); + //gPad->SetLogy(); + + delete convergence; return; // ******** // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995 - printf("Calculating covariance matrix. This may take some time...\n"); + /*printf("Calculating covariance matrix. This may take some time...\n"); + + // TODO check if this is the right one... + TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX()); Int_t xBins = hInverseResponseBayes->GetNbinsX(); Int_t yBins = hInverseResponseBayes->GetNbinsY(); @@ -1241,7 +1319,31 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]); new TCanvas; - cov->Draw("COLZ"); + cov->Draw("COLZ");*/ +} + +//____________________________________________________________________ +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) +{ + // + // helper function for the covariance matrix of the bayesian method + // + + Float_t result = 0; + + if (k == u && r == i) + result += 1.0 / hResponse->GetBinContent(u+1, r+1); + + if (k == u) + result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1); + + if (r == i) + result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1); + + result *= matrixM[k][i]; + + return result; } //____________________________________________________________________ @@ -1251,7 +1353,7 @@ void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullP // correct spectrum using bayesian method // - Float_t regPar = 0.05; + Float_t regPar = 0; Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); @@ -1299,7 +1401,7 @@ void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullP TH2F* hResponse = (TH2F*) fCurrentCorrelation; // unfold... - Int_t iterations = 20; + Int_t iterations = 25; for (Int_t i=0; iGetNbinsX(); t++) { - Float_t newValue = hTrueSmoothed->GetBinContent(t)/norm; + Float_t newValue = hTemp->GetBinContent(t)/norm; Float_t diff = hPrior->GetBinContent(t) - newValue; chi2 += (Double_t) diff * diff; @@ -1357,7 +1459,7 @@ void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullP printf("Chi2 of %d iteration = %.10f\n", i, chi2); - if (i % 5 == 0) + //if (i % 5 == 0) DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY()); delete hTrueSmoothed; @@ -1366,30 +1468,6 @@ void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullP DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY()); } -//____________________________________________________________________ -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 result = 0; - - if (k == u && r == i) - result += 1.0 / hResponse->GetBinContent(u+1, r+1); - - if (k == u) - result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1); - - if (r == i) - result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1); - - result *= matrixM[k][i]; - - return result; -} - //____________________________________________________________________ void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace) { @@ -1556,7 +1634,7 @@ void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id) if (id < 0 || id >= kESDHists) return; - TH2F* mc = fMultiplicityINEL[id]; + TH2F* mc = fMultiplicityVtx[id]; mc->Reset(); @@ -1571,5 +1649,8 @@ void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id) TH1* proj = mc->ProjectionY(); + TString tmp(fMultiplicityESD[id]->GetName()); + delete fMultiplicityESD[id]; fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id); + fMultiplicityESD[id]->SetName(tmp); } diff --git a/PWG0/dNdEta/AliMultiplicityCorrection.h b/PWG0/dNdEta/AliMultiplicityCorrection.h index 5ce87473dec..45c94f5fe52 100644 --- a/PWG0/dNdEta/AliMultiplicityCorrection.h +++ b/PWG0/dNdEta/AliMultiplicityCorrection.h @@ -24,7 +24,7 @@ class TCollection; class AliMultiplicityCorrection : public TNamed { public: enum EventType { kTrVtx = 0, kMB, kINEL }; - enum RegularizationType { kNone = 0, kPol0, kPol1, kCurvature, kEntropy }; + enum RegularizationType { kNone = 0, kPol0, kPol1, kCurvature, kEntropy, kTest }; AliMultiplicityCorrection(); AliMultiplicityCorrection(const Char_t* name, const Char_t* title); @@ -42,11 +42,12 @@ class AliMultiplicityCorrection : public TNamed { void DrawHistograms(); void DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist); - void ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, Bool_t check = kFALSE); - void ApplyMinuitFitAll(); + void ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* inputDist = 0); void SetRegularizationParameters(RegularizationType type, Float_t weight) { fRegularizationType = type; fRegularizationWeight = weight; }; - void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 0.07); + void ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace); + + void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 0.07, Int_t nIterations = 30, TH1* inputDist = 0); void ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace); @@ -58,22 +59,22 @@ class AliMultiplicityCorrection : public TNamed { TH2F* GetMultiplicityINEL(Int_t i) { return fMultiplicityINEL[i]; } TH2F* GetMultiplicityMC(Int_t i, EventType eventType); 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 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; } void SetGenMeasFromFunc(TF1* inputMC, Int_t id); TH2F* CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap, Bool_t normalized = kFALSE); - TH1F* GetMultiplicityESDCorrected(Int_t i) { return fMultiplicityESDCorrected[i]; } - static void NormalizeToBinWidth(TH1* hist); static void NormalizeToBinWidth(TH2* hist); - void GetComparisonResults(Float_t& mc, Float_t& residuals); + void GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals); protected: enum { kESDHists = 4, kMCHists = 5, kCorrHists = 8 }; @@ -84,8 +85,10 @@ class AliMultiplicityCorrection : public TNamed { static Double_t RegularizationPol1(Double_t *params); static Double_t RegularizationTotalCurvature(Double_t *params); static Double_t RegularizationEntropy(Double_t *params); + static Double_t RegularizationTest(Double_t *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); @@ -99,6 +102,8 @@ class AliMultiplicityCorrection : public TNamed { static TMatrixF* fCorrelationCovarianceMatrix; // contains the errors of fCurrentESD static TVectorF* fCurrentESDVector; // contains fCurrentESD + static TF1* fNBD; // negative binomial distribution + static RegularizationType fRegularizationType; // regularization that is used during Chi2 method static Float_t fRegularizationWeight; // factor for regularization term @@ -111,6 +116,7 @@ class AliMultiplicityCorrection : public TNamed { TH1F* fMultiplicityESDCorrected[kCorrHists]; // corrected histograms 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) private: diff --git a/PWG0/dNdEta/AliMultiplicityESDSelector.cxx b/PWG0/dNdEta/AliMultiplicityESDSelector.cxx index 4a57794db8d..7f26a66e5c2 100644 --- a/PWG0/dNdEta/AliMultiplicityESDSelector.cxx +++ b/PWG0/dNdEta/AliMultiplicityESDSelector.cxx @@ -2,23 +2,22 @@ #include "AliMultiplicityESDSelector.h" -#include -#include -#include #include -#include #include -#include +#include +#include +#include #include #include +#include #include "esdTrackCuts/AliESDtrackCuts.h" #include "AliPWG0Helper.h" +#include "dNdEta/AliMultiplicityCorrection.h" -#ifdef ALISELECTOR_USEMONALISA - #include -#endif +//#define TPCMEASUREMENT +#define ITSMEASUREMENT ClassImp(AliMultiplicityESDSelector) @@ -26,9 +25,6 @@ AliMultiplicityESDSelector::AliMultiplicityESDSelector() : AliSelector(), fMultiplicity(0), fEsdTrackCuts(0) -#ifdef ALISELECTOR_USEMONALISA - ,fMonaLisaWriter(0) -#endif { // // Constructor. Initialization of pointers @@ -78,25 +74,30 @@ void AliMultiplicityESDSelector::SlaveBegin(TTree* tree) ReadUserObjects(tree); - fMultiplicity = new TH1F("fMultiplicity", "multiplicity", 201, 0.5, 200.5); + fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); +} - #ifdef ALISELECTOR_USEMONALISA - TNamed *nm = 0; - if (fInput) - nm = dynamic_cast (fInput->FindObject("PROOF_QueryTag")); - if (!nm) - { - AliDebug(AliLog::kError, "Query tag not found. Cannot enable monitoring"); - return; - } +void AliMultiplicityESDSelector::Init(TTree* tree) +{ + // read the user objects - TString option = GetOption(); - option.ReplaceAll("#+", ""); + AliSelector::Init(tree); - TString id; - id.Form("%s_%s%d", gSystem->HostName(), nm->GetTitle(), gSystem->GetPid()); - fMonaLisaWriter = new TMonaLisaWriter(option, id, "CAF", "aliendb6.cern.ch"); - #endif + // enable only the needed branches + if (tree) + { + tree->SetBranchStatus("*", 0); + tree->SetBranchStatus("fTriggerMask", 1); + tree->SetBranchStatus("fSPDVertex*", 1); + + #ifdef ITSMEASUREMENT + tree->SetBranchStatus("fSPDMult*", 1); + #endif + + #ifdef TPCMEASUREMENT + AliESDtrackCuts::EnableNeededBranches(tree); + #endif + } } Bool_t AliMultiplicityESDSelector::Process(Long64_t entry) @@ -129,22 +130,103 @@ Bool_t AliMultiplicityESDSelector::Process(Long64_t entry) return kFALSE; } + Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD); + Bool_t eventVertex = AliPWG0Helper::IsVertexReconstructed(fESD); + + if (!eventTriggered || !eventVertex) + return kTRUE; + + // get the ESD vertex + const AliESDVertex* vtxESD = fESD->GetVertex(); + Double_t vtx[3]; + vtxESD->GetXYZ(vtx); + + Int_t nESDTracks05 = 0; + Int_t nESDTracks10 = 0; + Int_t nESDTracks15 = 0; + Int_t nESDTracks20 = 0; + +#ifdef ITSMEASUREMENT + // get tracklets + const AliMultiplicity* mult = fESD->GetMultiplicity(); + if (!mult) + { + AliDebug(AliLog::kError, "AliMultiplicity not available"); + return kFALSE; + } + + // get multiplicity from ITS tracklets + for (Int_t i=0; iGetNumberOfTracklets(); ++i) + { + //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i)); + + // this removes non-tracklets. Very bad solution. SPD guys are working on better solution... + if (mult->GetDeltaPhi(i) < -1000) + continue; + + Float_t theta = mult->GetTheta(i); + Float_t eta = -TMath::Log(TMath::Tan(theta/2.)); + + if (TMath::Abs(eta) < 0.5) + nESDTracks05++; + + if (TMath::Abs(eta) < 1.0) + nESDTracks10++; + + if (TMath::Abs(eta) < 1.5) + nESDTracks15++; + + if (TMath::Abs(eta) < 2.0) + nESDTracks20++; + } +#endif + +#ifdef TPCMEASUREMENT if (!fEsdTrackCuts) { AliDebug(AliLog::kError, "fESDTrackCuts not available"); return kFALSE; } - if (AliPWG0Helper::IsEventTriggered(fESD) == kFALSE) - return kTRUE; + // get multiplicity from ESD tracks + TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD); + Int_t nGoodTracks = list->GetEntries(); + // loop over esd tracks + for (Int_t i=0; i (list->At(i)); + if (!esdTrack) + { + AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i)); + continue; + } - if (AliPWG0Helper::IsVertexReconstructed(fESD) == kFALSE) - return kTRUE; + Double_t p[3]; + esdTrack->GetConstrainedPxPyPz(p); // ### TODO should be okay because we have a vertex, however GetInnerPxPyPy / GetOuterPxPyPy also exist + TVector3 vector(p); + + Float_t theta = vector.Theta(); + Float_t eta = -TMath::Log(TMath::Tan(theta/2.)); + Float_t pt = vector.Pt(); + + //if (pt < kPtCut) + // continue; - // get number of "good" tracks - Int_t nGoodTracks = fEsdTrackCuts->CountAcceptedTracks(fESD); + if (TMath::Abs(eta) < 0.5) + nESDTracks05++; - fMultiplicity->Fill(nGoodTracks); + if (TMath::Abs(eta) < 1.0) + nESDTracks10++; + + if (TMath::Abs(eta) < 1.5) + nESDTracks15++; + + if (TMath::Abs(eta) < 2.0) + nESDTracks20++; + } +#endif + + fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20); return kTRUE; } @@ -157,14 +239,6 @@ void AliMultiplicityESDSelector::SlaveTerminate() AliSelector::SlaveTerminate(); - #ifdef ALISELECTOR_USEMONALISA - if (fMonaLisaWriter) - { - delete fMonaLisaWriter; - fMonaLisaWriter = 0; - } - #endif - // Add the histograms to the output on each slave server if (!fOutput) { @@ -183,15 +257,19 @@ void AliMultiplicityESDSelector::Terminate() AliSelector::Terminate(); - fMultiplicity = dynamic_cast (fOutput->FindObject("fMultiplicity")); + fMultiplicity = dynamic_cast (fOutput->FindObject("Multiplicity")); if (!fMultiplicity) { - AliDebug(AliLog::kError, Form("ERROR: Histogram not available %p", (void*) fMultiplicity)); + AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity)); return; } TFile* file = TFile::Open("multiplicityESD.root", "RECREATE"); - fMultiplicity->Write(); + + fMultiplicity->SaveHistograms(); + file->Close(); + + fMultiplicity->DrawHistograms(); } diff --git a/PWG0/dNdEta/AliMultiplicityESDSelector.h b/PWG0/dNdEta/AliMultiplicityESDSelector.h index f818ad3c976..6487fefd674 100644 --- a/PWG0/dNdEta/AliMultiplicityESDSelector.h +++ b/PWG0/dNdEta/AliMultiplicityESDSelector.h @@ -5,15 +5,8 @@ #include "AliSelector.h" -// uncomment this to enable mona lisa monitoring -//#define ALISELECTOR_USEMONALISA - class AliESDtrackCuts; -class TH1F; - -#ifdef ALISELECTOR_USEMONALISA - class TMonaLisaWriter; -#endif +class AliMultiplicityCorrection; class AliMultiplicityESDSelector : public AliSelector { public: @@ -22,6 +15,7 @@ class AliMultiplicityESDSelector : public AliSelector { virtual void Begin(TTree* tree); virtual void SlaveBegin(TTree *tree); + virtual void Init(TTree *tree); virtual Bool_t Process(Long64_t entry); virtual void SlaveTerminate(); virtual void Terminate(); @@ -29,18 +23,13 @@ class AliMultiplicityESDSelector : public AliSelector { protected: void ReadUserObjects(TTree* tree); - TH1F* fMultiplicity; // multiplicity histogram - - AliESDtrackCuts* fEsdTrackCuts; // Object containing the parameters of the esd track cuts + AliMultiplicityCorrection* fMultiplicity; // object containing the extracted data + AliESDtrackCuts* fEsdTrackCuts; // Object containing the parameters of the esd track cuts private: AliMultiplicityESDSelector(const AliMultiplicityESDSelector&); AliMultiplicityESDSelector& operator=(const AliMultiplicityESDSelector&); -#ifdef ALISELECTOR_USEMONALISA - TMonaLisaWriter* fMonaLisaWriter; //! ML instance for monitoring -#endif - ClassDef(AliMultiplicityESDSelector, 0); }; diff --git a/PWG0/dNdEta/AliMultiplicityMCSelector.cxx b/PWG0/dNdEta/AliMultiplicityMCSelector.cxx index aae710d0327..dbecc097f52 100644 --- a/PWG0/dNdEta/AliMultiplicityMCSelector.cxx +++ b/PWG0/dNdEta/AliMultiplicityMCSelector.cxx @@ -96,7 +96,10 @@ void AliMultiplicityMCSelector::Init(TTree* tree) tree->SetBranchStatus("*", 0); tree->SetBranchStatus("fTriggerMask", 1); tree->SetBranchStatus("fSPDVertex*", 1); - tree->SetBranchStatus("fSPDMult*", 1); + + #ifdef ITSMEASUREMENT + tree->SetBranchStatus("fSPDMult*", 1); + #endif #ifdef TPCMEASUREMENT AliESDtrackCuts::EnableNeededBranches(tree); @@ -298,7 +301,7 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry) } #endif - fMultiplicity->FillMeasured(vtxMC[2], nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20); + fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20); // TODO should this be vtx[2] or vtxMC[2] ? fMultiplicity->FillCorrection(vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20); diff --git a/PWG0/dNdEta/AliMultiplicityMCSelector.h b/PWG0/dNdEta/AliMultiplicityMCSelector.h index 675643b8cae..3dfec893c9b 100644 --- a/PWG0/dNdEta/AliMultiplicityMCSelector.h +++ b/PWG0/dNdEta/AliMultiplicityMCSelector.h @@ -6,8 +6,6 @@ #include "AliSelectorRL.h" class AliESDtrackCuts; -class TH2F; -class TH3F; class AliMultiplicityCorrection; class AliMultiplicityMCSelector : public AliSelectorRL { diff --git a/PWG0/dNdEta/runMultiplicitySelector.C b/PWG0/dNdEta/runMultiplicitySelector.C index 33109a48df5..9dd5e60fb08 100644 --- a/PWG0/dNdEta/runMultiplicitySelector.C +++ b/PWG0/dNdEta/runMultiplicitySelector.C @@ -24,10 +24,6 @@ void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_ if (!prepareQuery(libraries, packages, kTRUE)) return; - gProof->SetParallel(0); - gProof->SetLogLevel(4); - gProof->SetParallel(9999); - gROOT->ProcessLine(".L CreateCuts.C"); gROOT->ProcessLine(".L drawPlots.C"); @@ -82,18 +78,39 @@ void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t ev TFile::Open(fileName); mult->LoadHistograms("Multiplicity"); - mult->ApplyBayesianMethod(hist, kFALSE, eventType); - mult->DrawComparison("Bayesian", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, eventType)->ProjectionY()); + //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; - //mult->ApplyMinuitFit(hist, kFALSE); - //mult->DrawComparison("MinuitChi2", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY()); + timer.Start(); + + mult->SetRegularizationParameters(AliMultiplicityCorrection::kTest, 1); + mult->ApplyMinuitFit(hist, kFALSE, AliMultiplicityCorrection::kTrVtx); + 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.root", const char* fileNameESD = "multiplicityESD.root", Int_t histID = 2) +void* fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fileNameESD = "multiplicityMC_3M.root", Int_t histID = 2) { gSystem->Load("libPWG0base"); @@ -105,16 +122,32 @@ void* fitOther(const char* fileNameMC = "multiplicityMC.root", const char* fileN TFile::Open(fileNameESD); TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID)); TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID)); + //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityMB%d", histID)); mult->SetMultiplicityESD(histID, hist); + mult->SetMultiplicityVtx(histID, hist2); + mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); + return; + + mult->SetRegularizationParameters(AliMultiplicityCorrection::kTest, 1.1); + mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); + mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY()); + + return; + //mult->ApplyGaussianMethod(histID, kFALSE); - //for (Float_t f=0; f<0.1; f+=0.01) - //mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); - mult->ApplyMinuitFit(histID, kFALSE); - mult->DrawComparison("MinuitChi2", hist, kFALSE, kTRUE, hist2->ProjectionY()); - //mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); + for (Float_t f=0.1; f<=0.11; f+=0.05) + { + mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, f); + mult->DrawComparison(Form("Bayesian_%f", f), histID, kFALSE, kTRUE, hist2->ProjectionY()); + } + + //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e7); + //mult->ApplyMinuitFit(histID, kFALSE); + //mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY()); + return mult; } @@ -128,12 +161,188 @@ 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; + } + return 0; +} + +const char* GetEventTypeName(Int_t type) +{ + switch (type) + { + case AliMultiplicityCorrection::kTrVtx: return "trigger, vertex"; break; + case AliMultiplicityCorrection::kMB: return "minimum bias"; break; + case AliMultiplicityCorrection::kINEL: return "inelastic"; break; } return 0; } -void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", Int_t histID = 2) +void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2) { + gSystem->mkdir(targetDir); + + gSystem->Load("libPWG0base"); + + AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); + TFile::Open(fileNameMC); + mult->LoadHistograms("Multiplicity"); + + AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); + TFile::Open(fileNameESD); + multESD->LoadHistograms("Multiplicity"); + mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); + + TCanvas* canvas = new TCanvas("EvaluateBayesianMethod", "EvaluateBayesianMethod", 800, 600); + TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98); + legend->SetFillColor(0); + + Float_t min = 1e10; + Float_t max = 0; + + TGraph* first = 0; + Int_t count = 0; // just to order the saved images... + + for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type) + { + TGraph* fitResultsMC = new TGraph; + fitResultsMC->SetTitle(";Weight Parameter"); + TGraph* fitResultsRes = new TGraph; + fitResultsRes->SetTitle(";Weight Parameter"); + + fitResultsMC->SetFillColor(0); + fitResultsRes->SetFillColor(0); + fitResultsMC->SetMarkerStyle(20+type); + fitResultsRes->SetMarkerStyle(24+type); + fitResultsRes->SetMarkerColor(kRed); + fitResultsRes->SetLineColor(kRed); + + legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type))); + legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type))); + + for (Float_t weight = 0; weight < 0.301; weight += 0.02) + { + Float_t chi2MC = 0; + Float_t residuals = 0; + + mult->ApplyBayesianMethod(histID, kFALSE, type, weight); + 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); + + fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC); + fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); + + min = TMath::Min(min, TMath::Min(chi2MC, residuals)); + max = TMath::Max(max, TMath::Max(chi2MC, residuals)); + } + + fitResultsMC->Print(); + fitResultsRes->Print(); + + canvas->cd(); + fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME")); + fitResultsRes->Draw("SAME CP"); + + if (first == 0) + first = fitResultsMC; + } + + gPad->SetLogy(); + printf("min = %f, max = %f\n", min, max); + if (min <= 0) + min = 1e-5; + first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5); + + legend->Draw(); + + canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); + canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); +} + +void EvaluateBayesianMethodIterations(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2) +{ + gSystem->mkdir(targetDir); + + gSystem->Load("libPWG0base"); + + AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); + TFile::Open(fileNameMC); + mult->LoadHistograms("Multiplicity"); + + AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); + TFile::Open(fileNameESD); + multESD->LoadHistograms("Multiplicity"); + mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); + + TCanvas* canvas = new TCanvas("EvaluateBayesianMethodIterations", "EvaluateBayesianMethodIterations", 800, 600); + TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98); + legend->SetFillColor(0); + + Float_t min = 1e10; + Float_t max = 0; + + TGraph* first = 0; + Int_t count = 0; // just to order the saved images... + + for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type) + { + TGraph* fitResultsMC = new TGraph; + fitResultsMC->SetTitle(";Iterations"); + TGraph* fitResultsRes = new TGraph; + fitResultsRes->SetTitle(";Iterations"); + + fitResultsMC->SetFillColor(0); + fitResultsRes->SetFillColor(0); + fitResultsMC->SetMarkerStyle(20+type); + fitResultsRes->SetMarkerStyle(24+type); + fitResultsRes->SetMarkerColor(kRed); + fitResultsRes->SetLineColor(kRed); + + legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type))); + legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type))); + + for (Int_t iter = 5; iter <= 50; iter += 5) + { + Float_t chi2MC = 0; + Float_t residuals = 0; + + mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1, iter); + mult->DrawComparison(Form("%s/BayesianIter_%02d_%d_%d", targetDir, count++, type, iter), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY()); + mult->GetComparisonResults(&chi2MC, 0, &residuals); + + fitResultsMC->SetPoint(fitResultsMC->GetN(), iter, chi2MC); + fitResultsRes->SetPoint(fitResultsRes->GetN(), iter, residuals); + + min = TMath::Min(min, TMath::Min(chi2MC, residuals)); + max = TMath::Max(max, TMath::Max(chi2MC, residuals)); + } + + fitResultsMC->Print(); + fitResultsRes->Print(); + + canvas->cd(); + fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME")); + fitResultsRes->Draw("SAME CP"); + + if (first == 0) + first = fitResultsMC; + } + + gPad->SetLogy(); + printf("min = %f, max = %f\n", min, max); + if (min <= 0) + min = 1e-5; + first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5); + + legend->Draw(); + + canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); + canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); +} + +void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2) +{ + gSystem->mkdir(targetDir); + gSystem->Load("libPWG0base"); AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); @@ -148,15 +357,19 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch mult->SetMultiplicityESD(histID, hist); TCanvas* canvas = new TCanvas("EvaluateChi2Method", "EvaluateChi2Method", 800, 600); - TLegend* legend = new TLegend(0.6, 0.1, 0.98, 0.4); + TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98); legend->SetFillColor(0); Float_t min = 1e10; Float_t max = 0; TGraph* first = 0; + Int_t count = 0; // just to order the saved images... + + Bool_t firstLoop = kTRUE; for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type) + //for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type) { TGraph* fitResultsMC = new TGraph; fitResultsMC->SetTitle(";Weight Parameter"); @@ -166,7 +379,7 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch fitResultsMC->SetFillColor(0); fitResultsRes->SetFillColor(0); fitResultsMC->SetMarkerStyle(19+type); - fitResultsRes->SetMarkerStyle(19+type); + fitResultsRes->SetMarkerStyle(23+type); fitResultsRes->SetMarkerColor(kRed); fitResultsRes->SetLineColor(kRed); @@ -176,15 +389,16 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch if (first == 0) first = fitResultsMC; - for (Float_t weight = 1e-2; weight < 2e4; weight *= 1e2) + for (Float_t weight = 1e-4; weight < 2e4; weight *= 10) + //for (Float_t weight = 0.1; weight < 10; weight *= TMath::Sqrt(TMath::Sqrt(10))) { Float_t chi2MC = 0; Float_t residuals = 0; mult->SetRegularizationParameters(type, weight); - mult->ApplyMinuitFit(histID, kFALSE); - mult->DrawComparison(Form("MinuitChi2_%d_%f", type, weight), histID, kFALSE, kTRUE, hist2->ProjectionY()); - mult->GetComparisonResults(chi2MC, residuals); + mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); + mult->DrawComparison(Form("%s/MinuitChi2_%02d_%d_%f", targetDir, count++, type, weight), histID, kFALSE, kTRUE, hist2->ProjectionY()); + mult->GetComparisonResults(&chi2MC, 0, &residuals); fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC); fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); @@ -197,8 +411,10 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch fitResultsRes->Print(); canvas->cd(); - fitResultsMC->Draw(Form("%s CP", (type == AliMultiplicityCorrection::kPol0) ? "A" : "SAME")); + fitResultsMC->Draw(Form("%s CP", (firstLoop) ? "A" : "SAME")); fitResultsRes->Draw("SAME CP"); + + firstLoop = kFALSE; } gPad->SetLogx(); @@ -210,7 +426,416 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch legend->Draw(); + canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); + canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); +} + +void EvaluateChi2MethodAll() +{ + EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M"); + EvaluateChi2Method("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M"); + EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD"); + EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M"); + EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD"); +} + +void EvaluateBayesianMethodAll() +{ + EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M"); + EvaluateBayesianMethod("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M"); + EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD"); + EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M"); + EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD"); +} + +void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2) +{ + gSystem->mkdir(targetDir); + + gSystem->Load("libPWG0base"); + + AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); + + TFile::Open(fileNameMC); + mult->LoadHistograms("Multiplicity"); + + TFile::Open(fileNameESD); + AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); + multESD->LoadHistograms("Multiplicity"); + + mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); + + TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 800); + canvas->Divide(3, 2); + + Int_t count = 0; + + for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type) + { + TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY(); + + mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3); + mult->ApplyMinuitFit(histID, kFALSE, type); + mult->DrawComparison(Form("%s/CompareMethods_%d_MinuitChi2", targetDir, type), histID, kFALSE, kTRUE, mc); + TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result"); + + mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1); + mult->DrawComparison(Form("%s/CompareMethods_%d_Bayesian", targetDir, type), histID, kFALSE, kTRUE, mc); + TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult"); + + mc->GetXaxis()->SetRangeUser(0, 150); + chi2Result->GetXaxis()->SetRangeUser(0, 150); + + // skip errors for now + for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i) + { + chi2Result->SetBinError(i, 0); + bayesResult->SetBinError(i, 0); + } + + canvas->cd(++count); + mc->SetFillColor(kYellow); + mc->DrawCopy(); + chi2Result->SetLineColor(kRed); + chi2Result->DrawCopy("SAME"); + bayesResult->SetLineColor(kBlue); + bayesResult->DrawCopy("SAME"); + gPad->SetLogy(); + + canvas->cd(count + 3); + chi2Result->Divide(chi2Result, mc, 1, 1, "B"); + bayesResult->Divide(bayesResult, mc, 1, 1, "B"); + + // skip errors for now + for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i) + { + chi2Result->SetBinError(i, 0); + bayesResult->SetBinError(i, 0); + } + + chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); + + chi2Result->DrawCopy(""); + bayesResult->DrawCopy("SAME"); + } + + canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); + canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); +} + +void StatisticsPlot(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 2) +{ + gSystem->Load("libPWG0base"); + + AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); + + TFile::Open(fileNameMC); + mult->LoadHistograms("Multiplicity"); + + const char* files[] = { "multiplicityMC_0.5M.root", "multiplicityMC_0.75M.root", "multiplicityMC_1M_3.root", "multiplicityMC_1.25M.root", "multiplicityMC_1.5M.root" }; + + TGraph* fitResultsChi2 = new TGraph; + fitResultsChi2->SetTitle(";Nevents;Chi2"); + TGraph* fitResultsBayes = new TGraph; + fitResultsBayes->SetTitle(";Nevents;Chi2"); + TGraph* fitResultsChi2Limit = new TGraph; + fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach"); + TGraph* fitResultsBayesLimit = new TGraph; + fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach"); + + TCanvas* canvas = new TCanvas("StatisticsPlot", "StatisticsPlot", 1200, 600); + canvas->Divide(5, 2); + + Float_t min = 1e10; + Float_t max = 0; + + for (Int_t i=0; i<5; ++i) + { + TFile::Open(files[i]); + AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); + multESD->LoadHistograms("Multiplicity"); + + mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); + Int_t nEntries = multESD->GetMultiplicityESD(histID)->GetEntries(); + TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(); + + mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3); + 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); + fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), nEntries, chi2MCLimit); + min = TMath::Min(min, chi2MC); + max = TMath::Max(max, chi2MC); + + TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result"); + + mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.1); + mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc); + TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult"); + mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); + fitResultsBayes->SetPoint(fitResultsBayes->GetN(), nEntries, chi2MC); + fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit); + + min = TMath::Min(min, chi2MC); + max = TMath::Max(max, chi2MC); + mc->GetXaxis()->SetRangeUser(0, 150); + chi2Result->GetXaxis()->SetRangeUser(0, 150); + + // skip errors for now + for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) + { + chi2Result->SetBinError(j, 0); + bayesResult->SetBinError(j, 0); + } + + canvas->cd(i+1); + mc->SetFillColor(kYellow); + mc->DrawCopy(); + chi2Result->SetLineColor(kRed); + chi2Result->DrawCopy("SAME"); + bayesResult->SetLineColor(kBlue); + bayesResult->DrawCopy("SAME"); + gPad->SetLogy(); + + canvas->cd(i+6); + chi2Result->Divide(chi2Result, mc, 1, 1, "B"); + bayesResult->Divide(bayesResult, mc, 1, 1, "B"); + + // skip errors for now + for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) + { + chi2Result->SetBinError(j, 0); + bayesResult->SetBinError(j, 0); + } + + chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC"); + chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); + + chi2Result->DrawCopy(""); + bayesResult->DrawCopy("SAME"); + } + + canvas->SaveAs(Form("%s.gif", canvas->GetName())); + canvas->SaveAs(Form("%s.C", canvas->GetName())); + + TCanvas* canvas2 = new TCanvas("StatisticsPlot2", "StatisticsPlot2", 800, 400); + 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"); + + gPad->SetLogy(); + + canvas2->cd(2); + fitResultsChi2Limit->SetMarkerStyle(20); + fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax())); + fitResultsChi2Limit->Draw("AP"); + + fitResultsBayesLimit->SetMarkerStyle(3); + fitResultsBayesLimit->SetMarkerColor(2); + fitResultsBayesLimit->Draw("P SAME"); + + canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); + canvas2->SaveAs(Form("%s.C", canvas2->GetName())); +} + +void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 2) +{ + gSystem->Load("libPWG0base"); + + AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); + + TFile::Open(fileNameMC); + mult->LoadHistograms("Multiplicity"); + + const char* files[] = { "multiplicityMC_0.5M.root", "multiplicityMC_0.75M.root", "multiplicityMC_1M_3.root", "multiplicityMC_1.25M.root", "multiplicityMC_1.5M.root" }; + + // this one we try to unfold + TFile::Open(files[0]); + AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); + multESD->LoadHistograms("Multiplicity"); + mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); + TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(); + + TGraph* fitResultsChi2 = new TGraph; + fitResultsChi2->SetTitle(";Input Dist ID;Chi2"); + TGraph* fitResultsBayes = new TGraph; + fitResultsBayes->SetTitle(";Input Dist ID;Chi2"); + TGraph* fitResultsChi2Limit = new TGraph; + fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach"); + TGraph* fitResultsBayesLimit = new TGraph; + fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach"); + + TCanvas* canvas = new TCanvas("StartingConditions", "StartingConditions", 1200, 600); + canvas->Divide(8, 2); + + TCanvas* canvas3 = new TCanvas("StartingConditions3", "StartingConditions3", 1000, 400); + canvas3->Divide(2, 1); + + Float_t min = 1e10; + Float_t max = 0; + + TH1* firstChi = 0; + TH1* firstBayesian = 0; + TH1* startCond = multESD->GetMultiplicityESD(histID)->ProjectionY("startCond"); + + TLegend* legend = new TLegend(0.7, 0.7, 1, 1); + + for (Int_t i=0; i<8; ++i) + { + if (i == 0) + { + startCond = (TH1*) mc->Clone("startCond2"); + } + else if (i < 6) + { + TFile::Open(files[i-1]); + AliMultiplicityCorrection* multESD2 = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2"); + multESD2->LoadHistograms("Multiplicity"); + startCond = multESD2->GetMultiplicityESD(histID)->ProjectionY("startCond"); + } + else if (i == 6) + { + 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, 50); + func->SetParNames("scaling", "averagen", "k"); + func->SetParLimits(0, 1e-3, 1e10); + func->SetParLimits(1, 0.001, 1000); + func->SetParLimits(2, 0.001, 1000); + + func->SetParameters(1, 10, 2); + for (Int_t j=2; j<=startCond->GetNbinsX(); j++) + startCond->SetBinContent(j, func->Eval(j-1)); + } + else if (i == 7) + { + for (Int_t j=1; j<=startCond->GetNbinsX(); j++) + startCond->SetBinContent(j, 1); + } + + mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3); + mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, startCond); + mult->DrawComparison(Form("StartingConditions_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc); + + Float_t chi2MC = 0; + Int_t chi2MCLimit = 0; + mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); + fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC); + fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit); + min = TMath::Min(min, chi2MC); + max = TMath::Max(max, chi2MC); + + TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result"); + if (!firstChi) + firstChi = (TH1*) chi2Result->Clone("firstChi"); + + mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.1, 30, startCond); + mult->DrawComparison(Form("StartingConditions_%d_Bayesian", i), histID, kFALSE, kTRUE, mc); + TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult"); + if (!firstBayesian) + firstBayesian = (TH1*) bayesResult->Clone("firstBayesian"); + + mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); + fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC); + fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit); + + min = TMath::Min(min, chi2MC); + max = TMath::Max(max, chi2MC); + mc->GetXaxis()->SetRangeUser(0, 150); + chi2Result->GetXaxis()->SetRangeUser(0, 150); + + // skip errors for now + for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) + { + chi2Result->SetBinError(j, 0); + bayesResult->SetBinError(j, 0); + } + + canvas3->cd(1); + TH1* tmp = (TH1*) chi2Result->Clone("tmp"); + tmp->SetTitle("Difference to best initial conditions;Npart;Ratio"); + tmp->Divide(firstChi); + tmp->GetYaxis()->SetRangeUser(0.5, 1.5); + tmp->GetXaxis()->SetRangeUser(0, 200); + tmp->SetLineColor(i+1); + legend->AddEntry(tmp, Form("%d", i)); + tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); + + canvas3->cd(2); + tmp = (TH1*) bayesResult->Clone("tmp"); + tmp->SetTitle("Difference to best initial conditions;Npart;Ratio"); + tmp->Divide(firstBayesian); + tmp->SetLineColor(i+1); + tmp->GetYaxis()->SetRangeUser(0.5, 1.5); + tmp->GetXaxis()->SetRangeUser(0, 200); + tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); + + canvas->cd(i+1); + mc->SetFillColor(kYellow); + mc->DrawCopy(); + chi2Result->SetLineColor(kRed); + chi2Result->DrawCopy("SAME"); + bayesResult->SetLineColor(kBlue); + bayesResult->DrawCopy("SAME"); + gPad->SetLogy(); + + canvas->cd(i+9); + chi2Result->Divide(chi2Result, mc, 1, 1, "B"); + bayesResult->Divide(bayesResult, mc, 1, 1, "B"); + + // skip errors for now + for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) + { + chi2Result->SetBinError(j, 0); + bayesResult->SetBinError(j, 0); + } + + chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC"); + chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); + + chi2Result->DrawCopy(""); + bayesResult->DrawCopy("SAME"); + } + + canvas3->cd(1); + legend->Draw(); + canvas->SaveAs(Form("%s.gif", canvas->GetName())); + + TCanvas* canvas2 = new TCanvas("StartingConditions2", "StartingConditions2", 800, 400); + 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"); + + gPad->SetLogy(); + + canvas2->cd(2); + fitResultsChi2Limit->SetMarkerStyle(20); + fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax())); + fitResultsChi2Limit->Draw("AP"); + + fitResultsBayesLimit->SetMarkerStyle(3); + fitResultsBayesLimit->SetMarkerColor(2); + fitResultsBayesLimit->Draw("P SAME"); + + canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); + canvas3->SaveAs(Form("%s.gif", canvas3->GetName())); } void Merge(Int_t n, const char** files, const char* output) @@ -264,17 +889,251 @@ void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root") 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(1000, 10, 2); break; - case 5: func->SetParameters(1000, 20, 3); break; - case 6: func->SetParameters(1000, 30, 4); break; - case 7: func->SetParameters(1000, 40, 5); break; + case 4: func->SetParameters(1e7, 10, 2); break; + case 5: func->SetParameters(1e7, 20, 3); break; + case 6: func->SetParameters(1e7, 30, 4); break; + case 7: func->SetParameters(1e7, 70, 2); break; + case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break; default: return; } mult->SetGenMeasFromFunc(func, 2); - mult->ApplyBayesianMethod(2, kFALSE); - mult->ApplyMinuitFit(2, kFALSE); + TFile::Open("out.root", "RECREATE"); + mult->SaveHistograms(); + + //mult->ApplyBayesianMethod(2, kFALSE); + //mult->ApplyMinuitFit(2, kFALSE); //mult->ApplyGaussianMethod(2, kFALSE); + mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx); +} + +void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t corrMatrix = 2) +{ + gSystem->Load("libPWG0base"); + + AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); + + TFile::Open(fileName); + mult->LoadHistograms("Multiplicity"); + + // empty under/overflow bins in x, otherwise Project3D takes them into account + TH3* corr = mult->GetCorrelation(corrMatrix); + for (Int_t j=1; j<=corr->GetYaxis()->GetNbins(); ++j) + { + for (Int_t k=1; k<=corr->GetZaxis()->GetNbins(); ++k) + { + corr->SetBinContent(0, j, k, 0); + corr->SetBinContent(corr->GetXaxis()->GetNbins()+1, j, k, 0); + } + } + + TH2* proj = (TH2*) corr->Project3D("zy"); + + // normalize correction for given nPart + for (Int_t i=1; i<=proj->GetNbinsX(); ++i) + { + Double_t sum = proj->Integral(i, i, 1, proj->GetNbinsY()); + if (sum <= 0) + continue; + + for (Int_t j=1; j<=proj->GetNbinsY(); ++j) + { + // npart sum to 1 + proj->SetBinContent(i, j, proj->GetBinContent(i, j) / sum); + proj->SetBinError(i, j, proj->GetBinError(i, j) / sum); + } + } + + new TCanvas; + proj->Draw("COLZ"); + + TH1* scaling = proj->ProjectionY("scaling", 1, 1); + scaling->Reset(); + scaling->SetMarkerStyle(3); + //scaling->GetXaxis()->SetRangeUser(0, 50); + TH1* mean = (TH1F*) scaling->Clone("mean"); + TH1* width = (TH1F*) scaling->Clone("width"); + + 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, 1, 1); + lognormal->SetParLimits(1, 0, 100); + lognormal->SetParLimits(2, 1e-3, 1); + + TF1* nbd = 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, 50); + nbd->SetParNames("scaling", "averagen", "k"); + nbd->SetParameters(1, 13, 5); + nbd->SetParLimits(0, 1, 1); + nbd->SetParLimits(1, 1, 100); + nbd->SetParLimits(2, 1, 1e8); + + TF1* poisson = new TF1("poisson", "[0] * exp(-(x+[2])) * (x+[2])**[1] / TMath::Factorial([1])", 0.01, 50); + poisson->SetParNames("scaling", "k", "deltax"); + poisson->SetParameters(1, 1, 0); + poisson->SetParLimits(0, 0, 10); + poisson->SetParLimits(1, 0.01, 100); + poisson->SetParLimits(2, 0, 10); + + TF1* mygaus = new TF1("mygaus", "[0] * exp(-(x-[1])**2 / 2 / [2] - [3] * log(x + [4]) / [5])", 0.01, 50); + mygaus->SetParNames("scaling", "mean", "width", "scale2log", "logmean", "logwidth"); + mygaus->SetParameters(1, 0, 1, 1, 0, 1); + mygaus->SetParLimits(2, 1e-5, 10); + mygaus->SetParLimits(4, 1, 1); + mygaus->SetParLimits(5, 1e-5, 10); + + //TF1* sqrt = new TF1("sqrt", "[0] + [1] * sqrt((x + [3]) * [2])", 0, 50); + TF1* sqrt = new TF1("sqrt", "[0] + (x + [1])**[2]", 0, 50); + sqrt->SetParNames("ydelta", "exp", "xdelta"); + sqrt->SetParameters(0, 0, 1); + sqrt->SetParLimits(1, 0, 10); + + const char* fitWith = "gaus"; + + for (Int_t i=1; i<=150; ++i) + { + printf("Fitting %d...\n", i); + + TH1* hist = proj->ProjectionY(Form("proj%d", i), i, i, "e"); + //hist->GetXaxis()->SetRangeUser(0, 50); + //lognormal->SetParameter(0, hist->GetMaximum()); + hist->Fit(fitWith, "0 M", ""); + + TF1* func = hist->GetListOfFunctions()->FindObject(fitWith); + + if (((i-1) % 15 == 0) || ((i % 5 == 0) && i < 30)) + { + new TCanvas; + hist->Draw(); + func->Clone()->Draw("SAME"); + gPad->SetLogy(); + } + + scaling->Fill(i, func->GetParameter(0)); + mean->Fill(i, func->GetParameter(1)); + width->Fill(i, func->GetParameter(2)); + } + + TF1* log = new TF1("log", "[0] + [1] * log([2] * x)", 0.01, 500); + log->SetParameters(0, 1, 1); + log->SetParLimits(1, 0, 100); + log->SetParLimits(2, 1e-3, 10); + + TF1* over = new TF1("over", "[0] + [1] / (x+[2])", 0.01, 500); + over->SetParameters(0, 1, 0); + //over->SetParLimits(0, 0, 100); + over->SetParLimits(1, 1e-3, 10); + over->SetParLimits(2, 0, 100); + + c1 = new TCanvas("fitparams", "fitparams", 1200, 400); + c1->Divide(3, 1); + + c1->cd(1); + scaling->Draw("P"); + + //TF1* scalingFit = new TF1("mypol0", "[0]"); + TF1* scalingFit = over; + scaling->Fit(scalingFit, "", "", 3, 100); + + c1->cd(2); + mean->Draw("P"); + + //TF1* meanFit = log; + TF1* meanFit = new TF1("mypol1", "[0]+[1]*x"); + mean->Fit(meanFit, "", "", 3, 100); + + c1->cd(3); + width->Draw("P"); + + //TF1* widthFit = over; + TF1* widthFit = new TF1("mypol2", "[0]+[1]*x+[2]*x*x"); + width->Fit(widthFit, "", "", 5, 100); + + // build new correction matrix + TH2* new = (TH2*) proj->Clone("new"); + new->Reset(); + Float_t x, y; + for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1) + { + TF1* func = (TF1*) gROOT->FindObject(fitWith); + x = new->GetXaxis()->GetBinCenter(i); + //if (i == 1) + // x = 0.1; + x++; + func->SetParameters(scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x)); + printf("%f %f %f %f\n", x, scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x)); + + for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1) + { + if (i < 21) + { + // leave bins 1..20 untouched + new->SetBinContent(i, j, corr->Integral(1, corr->GetNbinsX(), i, i, j, j)); + } + else + { + y = new->GetYaxis()->GetBinCenter(j); + if (j == 1) + y = 0.1; + if (func->Eval(y) > 1e-4) + new->SetBinContent(i, j, func->Eval(y)); + } + } + } + + // fill 0 multiplicity bins, this cannot be done with the function because it does not accept 0 + // we take the values from the old response matrix + //for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1) + // new->SetBinContent(i, 1, proj->GetBinContent(i, 1)); + + //for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1) + // new->SetBinContent(1, j, proj->GetBinContent(1, j)); + + // normalize correction for given nPart + for (Int_t i=1; i<=new->GetNbinsX(); ++i) + { + Double_t sum = new->Integral(i, i, 1, proj->GetNbinsY()); + if (sum <= 0) + continue; + + for (Int_t j=1; j<=new->GetNbinsY(); ++j) + { + // npart sum to 1 + new->SetBinContent(i, j, new->GetBinContent(i, j) / sum); + new->SetBinError(i, j, new->GetBinError(i, j) / sum); + } + } + + new TCanvas; + new->Draw("COLZ"); + + TH2* diff = (TH2*) new->Clone("diff"); + diff->Add(proj, -1); + + new TCanvas; + diff->Draw("COLZ"); + diff->SetMinimum(-0.05); + diff->SetMaximum(0.05); + + corr->Reset(); + + for (Int_t i=1; i<=new->GetNbinsX(); ++i) + for (Int_t j=1; j<=new->GetNbinsY(); ++j) + corr->SetBinContent(corr->GetXaxis()->GetNbins() / 2, i, j, new->GetBinContent(i, j)); + + new TCanvas; + corr->Project3D("zy")->Draw("COLZ"); + + TFile::Open("out.root", "RECREATE"); + mult->SaveHistograms(); + + TH1* proj1 = proj->ProjectionY("proj1", 36, 36); + TH1* proj2 = new->ProjectionY("proj2", 36, 36); + proj2->SetLineColor(2); + + new TCanvas; + proj1->Draw(); + proj2->Draw("SAME"); } -- 2.39.3