/* $Id$ */ // // plots for the note about multplicity measurements // #if !defined(__CINT__) || defined(__MAKECINT__) #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "AliMultiplicityCorrection.h" #include "AliCorrection.h" #include "AliCorrectionMatrix3D.h" #endif const char* correctionFile = "multiplicityMC.root"; const char* measuredFile = "multiplicityESD.root"; Int_t etaRange = 1; Int_t displayRange = 80; // axis range Int_t ratioRange = 151; // range to calculate difference Int_t longDisplayRange = 120; const char* correctionFileTPC = "multiplicityMC_TPC_1.4M.root"; const char* measuredFileTPC = "multiplicityMC_TPC_0.6M.root"; Int_t etaRangeTPC = 1; void loadlibs() { gSystem->Load("libANALYSIS"); gSystem->Load("libANALYSISalice"); gSystem->Load("libPWG0base"); } void SetTPC() { correctionFile = correctionFileTPC; measuredFile = measuredFileTPC; etaRange = etaRangeTPC; displayRange = 100; ratioRange = 76; longDisplayRange = 100; } const char* GetMultLabel(Int_t etaR = -1, Bool_t trueM = kTRUE) { if (etaR == -1) etaR = etaRange; TString tmpStr((trueM) ? "True " : "Measured "); tmpStr += "multiplicity"; //return Form("%s", tmpStr.Data()); if (etaR == 4) { tmpStr += " (full phase space)"; } else if (etaR == 2) { tmpStr += " in |#eta| < 1.3"; } else tmpStr += Form(" in |#eta| < %.1f", (etaR+1)* 0.5); return Form("%s", tmpStr.Data()); } void Smooth(TH1* hist, Int_t windowWidth = 20) { TH1* clone = (TH1*) hist->Clone("clone"); for (Int_t bin=2; bin<=clone->GetNbinsX(); ++bin) { Int_t min = TMath::Max(2, 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; } TH1* responseMatrixPlot(const char* fileName = 0, const char* folder = "Multiplicity") { loadlibs(); AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder); if (fileName == 0) fileName = correctionFile; TFile::Open(fileName); mult->LoadHistograms(); // empty under/overflow bins in x, otherwise Project3D takes them into account TH2* hist = (TH2*) mult->GetCorrelation(etaRange); for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y) { for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z) { hist->SetBinContent(0, y, z, 0); hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); } } hist = (TH2*) (((TH3*) mult->GetCorrelation(etaRange))->Project3D("zy")); hist->SetStats(kFALSE); if (0) { // normalize // normalize correction for given nPart for (Int_t i=1; i<=hist->GetNbinsX(); ++i) { Double_t sum = hist->Integral(i, i, 1, hist->GetNbinsY()); if (sum <= 0) continue; for (Int_t j=1; j<=hist->GetNbinsY(); ++j) { // npart sum to 1 hist->SetBinContent(i, j, hist->GetBinContent(i, j) / sum);// * correlation->GetXaxis()->GetBinWidth(i)); hist->SetBinError(i, j, hist->GetBinError(i, j) / sum); } } } hist->SetTitle(Form(";%s;%s;Entries", GetMultLabel(), GetMultLabel(etaRange, kFALSE))); hist->GetXaxis()->SetRangeUser(0, longDisplayRange); hist->GetYaxis()->SetRangeUser(0, longDisplayRange); hist->GetYaxis()->SetTitleOffset(1.3); hist->GetZaxis()->SetTitleOffset(1.2); TCanvas* canvas = new TCanvas("c1", "c1", 600, 600); canvas->SetRightMargin(0.15); canvas->SetTopMargin(0.05); gPad->SetLogz(); hist->Draw("COLZ"); canvas->SaveAs("responsematrix.eps"); return hist; } void multPythiaPhojet() { loadlibs(); TFile::Open("LHC08c11_10TeV_0.5T/mb1/spd/multiplicity.root"); AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); mult->LoadHistograms("Multiplicity"); hist1 = mult->GetMultiplicityINEL(1)->ProjectionY(); hist1->Sumw2(); TFile::Open("LHC08c15_10TeV_0.5T_Phojet/mb1/spd/multiplicity.root"); AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); mult->LoadHistograms("Multiplicity"); hist2 = mult->GetMultiplicityINEL(1)->ProjectionY(); hist2->Sumw2(); legend = new TLegend(0.6, 0.7, 0.9, 0.9); legend->SetFillColor(0); legend->SetTextSize(0.04); legend->AddEntry(hist1, "Pythia", "L"); legend->AddEntry(hist2, "Phojet", "L"); c1 = new TCanvas("c", "c", 600, 600); c1->SetTopMargin(0.05); c1->SetRightMargin(0.05); c1->SetLeftMargin(0.12); c1->SetGridx(); c1->SetGridy(); c1->SetLogy(); //hist1->SetMarkerStyle(20); //hist2->SetMarkerStyle(24); //hist2->SetMarkerColor(2); hist1->SetLineWidth(2); hist2->SetLineWidth(2); hist2->SetLineStyle(2); hist2->SetLineColor(2); hist1->Scale(1.0 / hist1->Integral()); hist2->Scale(1.0 / hist2->Integral()); hist1->SetStats(0); hist1->GetYaxis()->SetTitleOffset(1.3); hist1->GetXaxis()->SetRangeUser(0, 100); hist1->SetTitle(";N_{ch};P(N_{ch})"); hist1->Draw(""); hist2->Draw("SAME"); legend->Draw(); c1->SaveAs("mult_pythia_phojet.eps"); } TCanvas* DrawResultRatio(TH1* mcHist, TH1* result, TString epsName) { // normalize unfolded result to mc hist result->Scale(1.0 / result->Integral(2, displayRange)); result->Scale(mcHist->Integral(2, displayRange)); TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600); canvas->Range(0, 0, 1, 1); TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98); pad1->Draw(); TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 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(); pad1->SetGridx(); pad1->SetGridy(); 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, displayRange); mcHist->SetTitle(Form(";%s;Entries", GetMultLabel())); mcHist->SetStats(kFALSE); mcHist->DrawCopy("HIST E"); gPad->SetLogy(); result->SetLineColor(2); result->SetMarkerColor(2); result->SetMarkerStyle(5); result->DrawCopy("SAME PE"); TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9); legend->AddEntry(mcHist, "True distribution"); legend->AddEntry(result, "Unfolded distribution", "P"); 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<=ratioRange; ++i) { sum += TMath::Abs(ratio->GetBinContent(i) - 1); } sum /= ratioRange-1; printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum); TLine* line = new TLine(0, 1, displayRange, 1); line->SetLineWidth(2); line->Draw(); line = new TLine(0, 1.1, displayRange, 1.1); line->SetLineWidth(2); line->SetLineStyle(2); line->Draw(); line = new TLine(0, 0.9, displayRange, 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, displayRange)); result1->Scale(mcHist->Integral(2, displayRange)); result2->Scale(1.0 / result2->Integral(2, displayRange)); result2->Scale(mcHist->Integral(2, displayRange)); TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600); canvas->Range(0, 0, 1, 1); TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98); pad1->Draw(); TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 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(); gPad->SetGridx(); gPad->SetGridy(); 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, displayRange); mcHist->SetTitle(";true multiplicity;Entries"); mcHist->SetStats(kFALSE); mcHist->DrawCopy("HIST E"); gPad->SetLogy(); result1->SetLineColor(2); result1->SetMarkerStyle(24); result1->DrawCopy("SAME HISTE"); result2->SetLineColor(4); result2->SetMarkerColor(4); result2->DrawCopy("SAME HISTE"); TLegend* legend = new TLegend(0.5, 0.6, 0.95, 0.9); legend->AddEntry(mcHist, "True distribution"); legend->AddEntry(result1, "Unfolded distribution (syst)", "P"); legend->AddEntry(result2, "Unfolded distribution (normal)", "P"); legend->SetFillColor(0); legend->SetTextSize(0.06); legend->Draw(); pad2->cd(); pad2->SetBottomMargin(0.15); //gPad->SetGridx(); //gPad->SetGridy(); 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, displayRange); result1->SetTitle(Form(";%s;Entries", GetMultLabel())); result1->SetStats(kFALSE); // calculate ratio result1->Sumw2(); TH1* ratio = (TH1*) result1->Clone("ratio"); result2->Sumw2(); ratio->Divide(ratio, result2, 1, 1, ""); ratio->SetLineColor(1); ratio->SetMarkerColor(1); ratio->SetMarkerStyle(0); 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<=ratioRange; ++i) { sum += TMath::Abs(ratio->GetBinContent(i) - 1); } sum /= ratioRange-1; printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum); TLine* line = new TLine(0, 1, displayRange, 1); line->SetLineWidth(2); line->Draw(); line = new TLine(0, 1.1, displayRange, 1.1); line->SetLineWidth(2); line->SetLineStyle(2); line->Draw(); line = new TLine(0, 0.9, displayRange, 0.9); line->SetLineWidth(2); line->SetLineStyle(2); line->Draw(); canvas->Modified(); canvas->SaveAs(epsName); return canvas; } TCanvas* Draw2ResultRatios(TH1* mcHist, TH1* result1, TH1* result2, TH1* ratio2, TH1* ratio3, TString epsName) { // draws the 3 plots in the upper plot // draws the ratio between result1 and result2 in the lower plot // also draws ratio2 and ratio3 in the lower plot, uses their name for the legend // 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.Data()), "", 0, 0.5, 0.98, 0.98); pad1->Draw(); TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 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(); gPad->SetGridx(); gPad->SetGridy(); 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, displayRange); mcHist->SetTitle(";True multiplicity;Entries"); mcHist->SetStats(kFALSE); mcHist->DrawCopy("HIST E"); gPad->SetLogy(); result1->SetLineColor(2); result1->SetMarkerColor(2); result1->SetMarkerStyle(24); result1->DrawCopy("SAME HISTE"); result2->SetLineColor(4); result2->SetMarkerColor(4); result2->SetMarkerStyle(5); 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 (5%)", "P"); legend->AddEntry(result2, "Unfolded distribution (normal)", "P"); legend->SetFillColor(0); legend->SetTextSize(0.06); legend->Draw(); pad2->cd(); pad2->SetBottomMargin(0.15); //gPad->SetGridx(); //gPad->SetGridy(); 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, displayRange); result1->SetTitle(Form(";%s;Entries", GetMultLabel())); result1->SetStats(kFALSE); // calculate ratio result1->Sumw2(); TH1* ratio = (TH1*) result1->Clone("ratio"); result2->Sumw2(); ratio->Divide(ratio, result2, 1, 1, ""); ratio->SetLineColor(1); ratio->SetMarkerColor(1); ratio->SetMarkerStyle(24); ratio->GetYaxis()->SetTitle("Ratio (change / normal)"); ratio->GetYaxis()->SetRangeUser(0.55, 1.45); ratio2->SetLineColor(2); ratio2->SetMarkerColor(2); ratio2->SetMarkerStyle(25); ratio3->SetLineColor(4); ratio3->SetMarkerColor(4); ratio3->SetMarkerStyle(26); ratio->DrawCopy(); ratio2->DrawCopy("SAME"); ratio3->DrawCopy("SAME"); legend2 = new TLegend(0.3, 0.8, 0.8, 0.95); legend2->SetNColumns(3); legend2->SetFillColor(0); legend2->SetTextSize(0.06); legend2->AddEntry(ratio, "5% change", "P"); legend2->AddEntry(ratio2, "2% change", "P"); legend2->AddEntry(ratio3, "1% change", "P"); legend2->Draw(); // get average of ratio Float_t sum = 0; for (Int_t i=2; i<=ratioRange; ++i) { sum += TMath::Abs(ratio->GetBinContent(i) - 1); } sum /= ratioRange-1; printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum); TLine* line = new TLine(0, 1, displayRange, 1); line->SetLineWidth(2); line->Draw(); line = new TLine(0, 1.1, displayRange, 1.1); line->SetLineWidth(2); line->SetLineStyle(2); line->Draw(); line = new TLine(0, 0.9, displayRange, 0.9); line->SetLineWidth(2); line->SetLineStyle(2); line->Draw(); canvas->Modified(); canvas->SaveAs(epsName); return canvas; } void DrawEfficiencyChange() { loadlibs(); const char* fileName = "chi2a_inel.root"; base = AliMultiplicityCorrection::Open(Form("spd/%s", fileName)); low = AliMultiplicityCorrection::Open(Form("spd-loweff/%s", fileName)); high = AliMultiplicityCorrection::Open(Form("spd-higheff/%s", fileName)); const char* legendStrings[] = { "low efficiency", "high efficiency" }; file = TFile::Open("systunc_detectorefficiency.root", "RECREATE"); for (Int_t etaR = 1; etaR < 2; etaR++) { base->GetMultiplicityESDCorrected(etaR)->Scale(1.0 / base->GetMultiplicityESDCorrected(etaR)->Integral(2, base->GetMultiplicityESDCorrected(etaR)->GetNbinsX())); TH1* hists[2]; hists[0] = low->GetMultiplicityESDCorrected(etaR); hists[0]->Scale(1.0 / hists[0]->Integral(2, hists[0]->GetNbinsX())); if (high) { hists[1] = high->GetMultiplicityESDCorrected(etaR); hists[1]->Scale(1.0 / hists[1]->Integral(2, hists[1]->GetNbinsX())); } largestError = (TH1*) hists[0]->Clone(Form("detectorefficiency_%d", etaR)); largestError->Reset(); DrawRatio(base->GetMultiplicityESDCorrected(etaR), (high) ? 2 : 1, hists, Form("eff_%d.png", etaR), kFALSE, legendStrings, kFALSE, largestError); largestError->Write(); new TCanvas; largestError->DrawCopy(); } file->Close(); } TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0, Bool_t errors = kFALSE, TH1* largestErrorLow = 0, TH1* largestErrorHigh = 0) { // 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(1, 200)); TCanvas* canvas = new TCanvas(epsName, epsName, 800, 500); canvas->SetTopMargin(0.05); canvas->SetRightMargin(0.05); canvas->SetGridx(); canvas->SetGridy(); result->GetXaxis()->SetRangeUser(0, displayRange); result->GetYaxis()->SetRangeUser(0.55, 1.45); result->SetStats(kFALSE); // to get the axis how we want it TH1* dummy = (TH1*) result->Clone("dummy"); dummy->Reset(); dummy->SetTitle(Form(";%s;Ratio", GetMultLabel())); dummy->DrawCopy(); delete dummy; Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10}; TLegend* legend = new TLegend(0.2, 0.7, 0.7, 0.93); legend->SetFillColor(0); //legend->SetTextSize(0.04); if (nResultSyst > 6) legend->SetNColumns(2); for (Int_t n=0; nScale(1.0 / resultSyst[n]->Integral(1, 200)); // calculate ratio TH1* ratio = (TH1*) result->Clone("ratio"); ratio->Divide(ratio, resultSyst[n], 1, 1, ""); ratio->GetXaxis()->SetRangeUser(0, displayRange); if (firstMarker) ratio->SetMarkerStyle(5); ratio->SetLineColor(colors[n / 2]); if ((n % 2)) ratio->SetLineStyle(2); TString drawStr("SAME HIST"); if (n == 0 && firstMarker) drawStr = "SAME P"; if (errors) drawStr += " E"; ratio->DrawCopy(drawStr); if (legendStrings && legendStrings[n]) legend->AddEntry(ratio, legendStrings[n], "L"); /* s = new TSpline3(ratio); s->SetNpx(5); s->SetLineColor(kRed); s->Draw("same"); */ if (largestErrorLow && largestErrorHigh) for (Int_t i=1; i<=ratio->GetNbinsX(); ++i) { if (ratio->GetBinContent(i) - 1 > 0) largestErrorHigh->SetBinContent(i, TMath::Max(ratio->GetBinContent(i) - 1, largestErrorHigh->GetBinContent(i))); if (ratio->GetBinContent(i) - 1 < 0) largestErrorLow->SetBinContent(i, TMath::Min(ratio->GetBinContent(i) - 1, largestErrorLow->GetBinContent(i))); } if (largestErrorLow && !largestErrorHigh) for (Int_t i=1; i<=ratio->GetNbinsX(); ++i) largestErrorLow->SetBinContent(i, TMath::Max(TMath::Abs(ratio->GetBinContent(i) - 1), largestErrorLow->GetBinContent(i))); // get average of ratio Float_t sum = 0; for (Int_t i=2; i<=ratioRange; ++i) sum += TMath::Abs(ratio->GetBinContent(i) - 1); sum /= ratioRange-1; printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum); } if (legendStrings) legend->Draw(); TLine* line = new TLine(-0.5, 1, displayRange, 1); line->SetLineWidth(2); line->Draw(); line = new TLine(-0.5, 1.1, displayRange, 1.1); line->SetLineWidth(2); line->SetLineStyle(2); line->Draw(); line = new TLine(-0.5, 0.9, displayRange, 0.9); line->SetLineWidth(2); line->SetLineStyle(2); line->Draw(); canvas->SaveAs(epsName); return canvas; } TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName, Bool_t smooth = kFALSE, Bool_t dashed = kFALSE) { // draws the ratios of each mc to the corresponding result 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, displayRange); result[n]->SetStats(kFALSE); // calculate ratio TH1* ratio = (TH1*) result[n]->Clone("ratio"); ratio->Divide(mc[n], ratio, 1, 1, "B"); // SetRangeUser(1, ...) would be the same, but the 0 should be still on the axis... ratio->SetBinContent(1, 1); ratio->SetBinError(1, 0); if (smooth) Smooth(ratio); ratio->SetTitle(Form(";true multiplicity;Ratio (true / unfolded)%s", ((smooth) ? " (smoothed)" : ""))); ratio->GetYaxis()->SetRangeUser(0.55, 1.45); if (dashed) { ratio->SetLineColor((n/2)+1); ratio->SetLineStyle((n%2)+1); } else 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<=ratioRange; ++i) sum += TMath::Abs(ratio->GetBinContent(i) - 1); sum /= ratioRange-1; printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum); } TLine* line = new TLine(0, 1, displayRange, 1); line->SetLineWidth(2); line->Draw(); line = new TLine(0, 1.1, displayRange, 1.1); line->SetLineWidth(2); line->SetLineStyle(2); line->Draw(); line = new TLine(0, 0.9, displayRange, 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, displayRange); 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<=ratioRange; ++i) sum += TMath::Abs(ratio->GetBinContent(i)); sum /= ratioRange-1; printf("%d) Average (2..%d) of |ratio - ratioBase| is %f\n", n, ratioRange, sum); } TLine* line = new TLine(0, 0, displayRange, 0); line->SetLineWidth(2); line->Draw(); line = new TLine(0, 0.1, displayRange, 0.1); line->SetLineWidth(2); line->SetLineStyle(2); line->Draw(); line = new TLine(0, -0.1, displayRange, -0.1); line->SetLineWidth(2); line->SetLineStyle(2); line->Draw(); canvas->Modified(); canvas->SaveAs(epsName); canvas->SaveAs(Form("%s.gif", epsName.Data())); return canvas; } 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, displayRange); 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(); // clear 0 bin ratio->SetBinContent(1, 0); ratio->SetBinError(1, 0); Smooth(ratio); //ratio->SetLineColor(1); ratio->DrawCopy("SAME"); canvas->cd(); ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u) (smoothed)"); ratio->GetYaxis()->SetRangeUser(-0.3, 0.3); ratio->SetLineColor((n / 2)+1); ratio->SetLineStyle((n % 2)+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, displayRange, 0); line->SetLineWidth(2); line->Draw(); line = new TLine(0, 0.1, displayRange, 0.1); line->SetLineWidth(2); line->SetLineStyle(2); line->Draw(); line = new TLine(0, -0.1, displayRange, -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(const char* fileName, const char* epsName) { loadlibs(); AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileName); TH1* measured = mult->GetMultiplicityESD(etaRange)->ProjectionY("myesd", 1, 1); if (0) { // test how far we are from a normal exponential in the unfolded corrected = mult->GetMultiplicityESDCorrected(etaRange); new TCanvas; corrected->DrawCopy(); gPad->SetLogy(); expo = new TF1("exp", "expo(0)", 0, 50); //expo->SetParameters(10, -0.18); expo->SetParameters(9.96, -0.176); expo->Draw("SAME"); for (Int_t i=21; i<=50; i++) corrected->SetBinContent(i, expo->Eval(corrected->GetXaxis()->GetBinCenter(i))); corrected->SetLineColor(2); corrected->Draw("SAME"); } TH1* unfoldedFolded = mult->GetConvoluted(etaRange, AliMultiplicityCorrection::kINEL); //mult->CalculateMultiplicityESD(mult->GetMultiplicityESDCorrected(etaRange), etaRange)->ProjectionY("myfolded", 1, 1); // normalize //unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, displayRange+1)); //unfoldedFolded->Scale(measured->Integral(2, displayRange+1)); 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, 1, 1); pad1->Draw(); pad1->SetGridx(); pad1->SetGridy(); TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 1, 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, displayRange); measured->SetTitle(Form(";%s;Entries", GetMultLabel(etaRange, kFALSE))); measured->SetStats(kFALSE); measured->DrawCopy("HIST"); gPad->SetLogy(); unfoldedFolded->SetMarkerStyle(5); unfoldedFolded->SetMarkerColor(2); unfoldedFolded->SetLineColor(2); unfoldedFolded->DrawCopy("SAME PHIST"); TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9); legend->AddEntry(measured, "Measured distribution", "L"); legend->AddEntry(unfoldedFolded, "R #otimes unfolded distribution", "P"); legend->SetFillColor(0); legend->SetTextSize(0.06); 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", ";", 16, -3, 3); residualHist->Sumw2(); Float_t chi2 = 0; Float_t chi2_hump = 0; for (Int_t i=1; i<=displayRange+1; ++i) { if (measured->GetBinError(i) > 0) { residual->SetBinContent(i, residual->GetBinContent(i) / measured->GetBinError(i)); residual->SetBinError(i, 1); residualHist->Fill(residual->GetBinContent(i)); if (i >= 15 && i <= 23) chi2_hump += residual->GetBinContent(i) * residual->GetBinContent(i); chi2 += residual->GetBinContent(i) * residual->GetBinContent(i); } else { residual->SetBinContent(i, 0); residual->SetBinError(i, 0); } } Printf("chi2 / ndf = %f / %d = %f", chi2, displayRange+1, chi2 / (displayRange+1)); Printf("chi2 (hump) / ndf = %f / %d = %f", chi2_hump, 23-15+1, chi2_hump / (23-15+1)); 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, displayRange + 0.5, 0); line->SetLineWidth(2); line->Draw(); pad3->cd(); residualHist->SetStats(kFALSE); residualHist->SetLabelSize(0.08, "xy"); residualHist->Fit("gaus"); residualHist->Draw("HIST"); residualHist->FindObject("gaus")->Draw("SAME"); 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 chi2FluctuationResult() { loadlibs(); TFile::Open("chi2_noregularization.root"); AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); mult->LoadHistograms("Multiplicity"); TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1); mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE); TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_1"); ((TH1*) canvas->FindObject("proj"))->GetXaxis()->SetRangeUser(0, displayRange); ((TH1*) canvas->FindObject("fCurrentESD"))->GetXaxis()->SetRangeUser(0, displayRange); //((TH1*) canvas->FindObject("proj"))->GetXaxis()->SetTitle(GetMultTitle()); //((TH1*) canvas->FindObject("fCurrentESD"))->GetXaxis()->->SetTitle(GetMultTitle(etaRange, kFALSE)); canvas->SaveAs("chi2FluctuationResult.eps"); } void DrawUnfolded(const char* fileName, const char* eps) { loadlibs(); TFile::Open(fileName); AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); mult->LoadHistograms("Multiplicity"); TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, mult->GetMultiplicityVtx(etaRange)->GetNbinsX()); TH1* result = mult->GetMultiplicityESDCorrected(etaRange); DrawResultRatio(mcHist, result, eps); } void minimizationInfluenceAlpha() { loadlibs(); TFile::Open(measuredFile); AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2"); mult2->LoadHistograms("Multiplicity"); TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, mult2->GetMultiplicityVtx(etaRange)->GetNbinsX()); mcHist->Scale(1.0 / mcHist->Integral()); mcHist->SetStats(kFALSE); mcHist->SetTitle(";True multiplicity n in |#eta| < 1.0;P(n)"); TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 350); canvas->Divide(3, 1, 0.005); TFile::Open("chi2compare-influencealpha/EvaluateChi2Method.root"); TH1* hist1 = (TH1*) gFile->Get("MinuitChi2_00_2_3.162278"); TH1* hist2 = (TH1*) gFile->Get("MinuitChi2_07_2_10000.000000"); TH1* hist3 = (TH1*) gFile->Get("MinuitChi2_13_2_10000000.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, displayRange); mcHist->SetLabelSize(0.06, "xy"); mcHist->SetTitleSize(0.06, "xy"); mcHist->GetYaxis()->SetTitleOffset(1.5); canvas->cd(1); gPad->SetLogy(); gPad->SetRightMargin(0.03); gPad->SetLeftMargin(0.19); gPad->SetTopMargin(0.05); gPad->SetBottomMargin(0.13); gPad->SetGridx(); gPad->SetGridy(); mcHist->Draw(); hist1->SetMarkerStyle(5); hist1->SetMarkerColor(2); hist1->Draw("SAME PE"); canvas->cd(2); gPad->SetRightMargin(0.03); gPad->SetLeftMargin(0.19); gPad->SetTopMargin(0.05); gPad->SetBottomMargin(0.13); gPad->SetLogy(); gPad->SetGridx(); gPad->SetGridy(); mcHist->Draw(); hist2->SetMarkerStyle(5); hist2->SetMarkerColor(2); hist2->Draw("SAME PE"); canvas->cd(3); gPad->SetRightMargin(0.03); gPad->SetLeftMargin(0.19); gPad->SetTopMargin(0.05); gPad->SetBottomMargin(0.13); gPad->SetLogy(); gPad->SetGridx(); gPad->SetGridy(); 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 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, 1200, 600); canvas->Divide(2, 1); TH1* mc = (TH1*) file->Get("mc"); mc->Sumw2(); mc->Scale(1.0 / mc->Integral(2, displayRange)); //Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5}; TLegend* legend = new TLegend(0.6, 0.7, 0.99, 0.99); legend->SetFillColor(0); legend->SetTextSize(0.04); const char* names[] = { "True", "Measured 1", "Measured 2", "Measured 3", "NBD", "Flat" }; 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->Scale(1.0 / chi2Result->Integral(2, displayRange)); bayesResult->Scale(1.0 / bayesResult->Integral(2, displayRange)); chi2Result->Divide(mc, chi2Result, 1, 1, ""); bayesResult->Divide(mc, bayesResult, 1, 1, ""); chi2Result->SetTitle(Form("a) #chi^{2}-minimization;%s;MC / unfolded", GetMultLabel())); chi2Result->GetXaxis()->SetRangeUser(0, displayRange); chi2Result->GetYaxis()->SetRangeUser(0.7, 1.3); chi2Result->GetYaxis()->SetTitleOffset(1.7); //chi2Result->SetMarkerStyle(marker[i]); chi2Result->SetLineColor(i+1); chi2Result->SetMarkerColor(i+1); chi2Result->SetStats(kFALSE); bayesResult->SetTitle(Form("b) Bayesian unfolding;%s;MC / unfolded", GetMultLabel())); bayesResult->GetXaxis()->SetRangeUser(0, displayRange); bayesResult->GetYaxis()->SetRangeUser(0.7, 1.3); bayesResult->GetYaxis()->SetTitleOffset(1.7); bayesResult->SetStats(kFALSE); //bayesResult->SetLineColor(2); bayesResult->SetLineColor(i+1); canvas->cd(1); gPad->SetRightMargin(0.05); gPad->SetLeftMargin(0.12); gPad->SetGridx(); gPad->SetGridy(); chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME"); canvas->cd(2); gPad->SetRightMargin(0.05); gPad->SetLeftMargin(0.12); gPad->SetGridx(); gPad->SetGridy(); bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME"); //TLine* line = new TLine(0, 1, 150, 1); //line->Draw(); legend->AddEntry(chi2Result, names[i]); } canvas->cd(1); legend->Draw(); canvas->cd(2); legend->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 = (TGraph*) 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"; name = "StatisticsPlotRatios"; 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 Draw2Unfolded(const char* file1, const char* file2, const char* output) { loadlibs(); TFile::Open(file1); AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); mult->LoadHistograms("Multiplicity"); // result with systematic effect TFile::Open(file2); AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2"); mult2->LoadHistograms("Multiplicity"); TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1); TH1* result1 = (TH1*) mult2->GetMultiplicityESDCorrected(etaRange)->Clone("result1"); // from file2 (with syst) TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2"); // from file1 (without syst) DrawResultRatio(mcHist, result1, "tmp1.eps"); DrawResultRatio(mcHist, result2, "tmp2.eps"); Draw2ResultRatio(mcHist, result1, result2, output); } void PythiaPhojet() { loadlibs(); displayRange = 55; Draw2Unfolded("self.root", "pythia.root", "test.eps"); canvas = (TCanvas*) gROOT->GetListOfCanvases()->Last(); pad1 = (TCanvas*)canvas->GetListOfPrimitives()->First(); pad2 = (TCanvas*)canvas->GetListOfPrimitives()->Last(); legend = (TLegend*)pad1->GetListOfPrimitives()->Last(); ((TH1*)pad2->GetListOfPrimitives()->At(1))->GetYaxis()->SetTitle("Ratio (Pythia / Phojet)"); ((TLegendEntry*)legend->GetListOfPrimitives()->At(1))->SetLabel("Unfolded distribution (Pythia)"); ((TLegendEntry*)legend->GetListOfPrimitives()->At(2))->SetLabel("Unfolded distribution (Phojet)"); canvas->SaveAs("PythiaPhojet.eps"); } void Misalignment() { loadlibs(); Draw2Unfolded("chi2_ideal.root", "chi2_misaligned.root", "test.eps"); canvas = (TCanvas*) gROOT->GetListOfCanvases()->Last(); pad1 = (TCanvas*)canvas->GetListOfPrimitives()->First(); pad2 = (TCanvas*)canvas->GetListOfPrimitives()->Last(); legend = (TLegend*)pad1->GetListOfPrimitives()->Last(); ((TH1*)pad2->GetListOfPrimitives()->At(1))->GetYaxis()->SetTitle("Ratio (misaligned / realigned)"); ((TLegendEntry*)legend->GetListOfPrimitives()->At(1))->SetLabel("Unfolded distribution (misaligned)"); ((TLegendEntry*)legend->GetListOfPrimitives()->At(2))->SetLabel("Unfolded distribution (realigned)"); canvas->SaveAs("SystematicMisalignment.eps"); } void SystematicLowEfficiency() { Draw2Unfolded("chi2.root", "chi2_loweff_5.root", "SystematicLowEfficiency.eps"); } void SystematicLowEfficiency2() { loadlibs(); AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("chi2.root"); TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2"); TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1); result2->Scale(1.0 / result2->Integral(2, displayRange)); result2->Scale(mcHist->Integral(2, displayRange)); mult = AliMultiplicityCorrection::Open("chi2_loweff_5.root"); TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1"); mult = AliMultiplicityCorrection::Open("chi2_loweff_2.root"); TH1* ratio2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("ratio2"); ratio2->Scale(1.0 / ratio2->Integral(2, displayRange)); ratio2->Scale(mcHist->Integral(2, displayRange)); ratio2->Divide(result2); mult = AliMultiplicityCorrection::Open("chi2_loweff_1.root"); TH1* ratio3 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("ratio3"); ratio3->Scale(1.0 / ratio3->Integral(2, displayRange)); ratio3->Scale(mcHist->Integral(2, displayRange)); ratio3->Divide(result2); Draw2ResultRatios(mcHist, result1, result2, ratio2, ratio3, "SystematicLowEfficiency2.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 SystematicMisalignmentTPC() { gSystem->Load("libPWG0base"); SetTPC(); TFile::Open(correctionFile); AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); mult->LoadHistograms("Multiplicity"); TFile::Open("multiplicityMC_TPC_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, "SystematicMisalignmentTPC.eps"); } void LowMomentumEffectSPD() { // this function increases/reduces the correction as function of pt between 0 and 0.2 gev by +-50% to 0% and checks the effect on the overall correction factor // only a normal acceptance region is considered to not get a bias by the edges loadlibs(); TFile::Open("multiplicity.root"); AliCorrection* correction[8]; Float_t values[3]; for (Int_t loop=0; loop<3; loop++) { Float_t sumGen = 0; Float_t sumMeas = 0; Printf("loop %d", loop); 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(); TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram(); TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram(); Float_t vtxRange = 5.9; gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange); meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange); Float_t etaRange = 0.99; gene->GetYaxis()->SetRangeUser(-etaRange, etaRange); meas->GetYaxis()->SetRangeUser(-etaRange, etaRange); TH1* genePt = gene->Project3D(Form("z_%d", i)); TH1* measPt = meas->Project3D(Form("z_%d", i)); if (loop > 0) { for (Int_t x=1; x<=genePt->GetNbinsX(); x++) { Float_t pt = genePt->GetXaxis()->GetBinCenter(x); //Printf("%f", pt); if (pt < 0.2) { Float_t factor = 1; if (loop == 1) factor = 1.5 - pt / 0.2 * 0.5; if (loop == 2) factor = 0.5 + pt / 0.2 * 0.5; //Printf("%f", factor); genePt->SetBinContent(x, genePt->GetBinContent(x) * factor); measPt->SetBinContent(x, measPt->GetBinContent(x) * factor); } } } //new TCanvas; genePt->DrawCopy(); measPt->DrawCopy("SAME"); sumGen += genePt->Integral(); sumMeas += measPt->Integral(); Float_t average = measPt->Integral() / genePt->Integral(); Printf("The average efficiency of this correction is %f", average); } Float_t average = sumMeas / sumGen; Printf("The average efficiency of all corrections is %f", average); values[loop] = average; } Printf("relative is %f and %f", values[1] / values[0], values[2] / values[0]); } void EfficiencySpecies(Bool_t addDecayStopped = kFALSE) { loadlibs(); Int_t marker[] = {24, 25, 26, 27}; Int_t color[] = {1, 2, 4, 3}; // SPD TPC //const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" }; //const char* fileName[] = { "LHC09b9_0.9TeV_0T/mb1/spd/multiplicity.root", "LHC09b8_0.9TeV_0.5T/mb1/tpc/multiplicity.root" }; //const char* fileName[] = { "spd/multiplicity.root", "tpc/multiplicity.root" }; const char* fileName[] = { "multiplicity.root", "multiplicity.root" }; Float_t etaRangeArr[] = {0.49, 0.9}; const char* titles[] = { "SPD Tracklets", "TPC Tracks" }; TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500); canvas->Divide(2, 1); TCanvas* canvas3 = new TCanvas("EfficiencySpecies_comb", "EfficiencySpecies_comb", 600, 600); gPad->SetGridx(); gPad->SetGridy(); gPad->SetRightMargin(0.05); gPad->SetTopMargin(0.05); TLegend* legends[2]; TH1* effGlobal = 0; for (Int_t loop=0; loop<1; ++loop) { Printf("%s", fileName[loop]); TCanvas* canvas2 = new TCanvas(Form("EfficiencySpecies_%d", loop), Form("EfficiencySpecies_%d", loop), 600, 600); gPad->SetGridx(); gPad->SetGridy(); gPad->SetRightMargin(0.05); gPad->SetTopMargin(0.05); AliCorrection* correction[8]; canvas->cd(loop+1); gPad->SetGridx(); gPad->SetGridy(); gPad->SetRightMargin(0.05); TLegend* legend = new TLegend(0.6, 0.4, 0.85, 0.6); legend->SetFillColor(0); legend->SetEntrySeparation(0.2); legend->SetTextSize(gStyle->GetTextSize()); legends[loop] = new TLegend(0.4+loop*0.3, 0.2, 0.6+loop*0.3, 0.5); legends[loop]->SetFillColor(0); legends[loop]->SetEntrySeparation(0.2); legends[loop]->SetTextSize(gStyle->GetTextSize()); legends[loop]->SetHeader((loop == 0) ? "SPD" : "TPC"); Float_t below = 0; Float_t total = 0; TFile* file = TFile::Open(fileName[loop]); if (!file) { Printf("Could not open %s", fileName[loop]); return; } Float_t sumGen = 0; Float_t sumMeas = 0; Float_t sumGenAbove02 = 0; Float_t sumMeasAbove02 = 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 Float_t vtxRange = 3.9; gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange); meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange); // 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 Float_t etaBegin = -etaRangeArr[loop]; Float_t etaEnd = etaRangeArr[loop]; //etaBegin = 0.01; //etaEnd = -0.01; gene->GetYaxis()->SetRangeUser(etaBegin, etaEnd); meas->GetYaxis()->SetRangeUser(etaBegin, etaEnd); TH1* genePt = gene->Project3D(Form("z_%d", i)); TH1* measPt = meas->Project3D(Form("z_%d", i)); // if (i == 2) // { // Printf("WARNING: Rebinning for protons!"); // genePt->Rebin(2); // measPt->Rebin(2); // } genePt->Sumw2(); measPt->Sumw2(); for (Int_t x=0; x<=genePt->GetNbinsX()+1; x++) { genePt->SetBinError(x, TMath::Sqrt(genePt->GetBinContent(x))); measPt->SetBinError(x, TMath::Sqrt(measPt->GetBinContent(x))); } sumGen += genePt->Integral(); sumMeas += measPt->Integral(); sumGenAbove02 += genePt->Integral(genePt->GetXaxis()->FindBin(0.21), genePt->GetNbinsX()); sumMeasAbove02 += measPt->Integral(genePt->GetXaxis()->FindBin(0.21), genePt->GetNbinsX()); TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i)); effPt->Reset(); effPt->Divide(measPt, genePt, 1, 1, "B"); if (i == 0) effGlobal = effPt; 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->SetMarkerSize(2); effPt->GetXaxis()->SetRangeUser(0, 1); effPt->GetYaxis()->SetRangeUser(0.001, 1); effPt->GetXaxis()->SetTitleOffset(1.1); effPt->GetYaxis()->SetTitleOffset(1.2); effPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); effPt->SetStats(kFALSE); effPt->SetTitle(titles[loop]); effPt->GetYaxis()->SetTitle("Efficiency"); canvas->cd(loop+1); effPt->DrawCopy((i == 0) ? "" : "SAME"); canvas2->cd(); effPt->SetTitle(""); effPt->DrawCopy((i == 0) ? "" : "SAME"); canvas3->cd(); effPtClone = (TH1*) effPt->Clone("effPtClone"); effPtClone->SetMarkerStyle(marker[i]-4*loop); effPtClone->DrawCopy((i == 0 && loop == 0) ? "" : "SAME"); legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")), "P"); legends[loop]->AddEntry(effPtClone, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")), "P"); //legend2->AddEntry(effPt, Form("%s %s", (loop == 0) ? "SPD" : "TPC", ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}"))), "P"); if (addDecayStopped) { name.Form("correction_%d", i+4); corr = new AliCorrection(name, name); corr->LoadHistograms(); TH3* gene = corr->GetTrackCorrection()->GetGeneratedHistogram(); TH3* meas = corr->GetTrackCorrection()->GetMeasuredHistogram(); // limit axes gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange); meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange); gene->GetYaxis()->SetRangeUser(-etaRangeArr[loop], etaRangeArr[loop]); meas->GetYaxis()->SetRangeUser(-etaRangeArr[loop], etaRangeArr[loop]); TH1* decayed = gene->Project3D(Form("z_%d", i+4)); TH1* stopped = meas->Project3D(Form("z_%d", i+4)); Printf("%d: %d decayed, %d stopped, out of %d", i, (Int_t) decayed->Integral(), (Int_t) stopped->Integral(), (Int_t) genePt->Integral()); decayed->Divide(decayed, genePt, 1, 1, "B"); stopped->Divide(stopped, genePt, 1, 1, "B"); decayed->SetMarkerStyle(20); stopped->SetMarkerStyle(21); stopped->SetMarkerColor(2); new TCanvas(Form("all_%d_%d", loop, i), Form("all_%d_%d", loop, i), 600, 600); effPt->DrawCopy(); decayed->DrawCopy("SAME"); stopped->DrawCopy("SAME"); decayed->Add(stopped); decayed->Add(effPt); decayed->SetMarkerStyle(22); decayed->SetMarkerColor(4); decayed->DrawCopy("SAME"); } } Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total); Printf("%f measured, %f generated, effiency: %f", sumGen, sumMeas, sumMeas / sumGen); Printf("Above 0.2 GeV/c: %f measured, %f generated, effiency: %f", sumGenAbove02, sumMeasAbove02, sumMeasAbove02 / sumGenAbove02); canvas->cd(loop+1); legend->Draw(); canvas2->cd(); legend->Draw(); canvas2->SaveAs(Form("%s.eps", canvas2->GetName())); } TH1* pt = 0; for (Int_t i=0; i<3; ++i) { TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram(); // limit axes gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange); gene->GetYaxis()->SetRangeUser(-0.49, 0.49); if (!pt) pt = gene->Project3D(Form("z_pt_%d", i)); else pt->Add(gene->Project3D(Form("z_pt_%d", i))); } new TCanvas; Float_t sum = 0; for (Int_t i=1; i<30; i++) { Float_t tmpCorr = 100.0 * pt->GetBinContent(i) / pt->Integral() * (0.532348 - effGlobal->GetBinContent(i)) / 0.532348; Printf("Yield in %d bin: %.1f%%, efficiency %.1f%%, correction: %.2f%%", i, 100.0 * pt->GetBinContent(i) / pt->Integral(), 100.0 * effGlobal->GetBinContent(i), tmpCorr); sum += tmpCorr; } Printf("%f", sum); AliPWG0Helper::NormalizeToBinWidth(pt); pt->Draw(); return; canvas->cd(); canvas->SaveAs(Form("%s.eps", canvas->GetName())); canvas3->cd(); legends[0]->Draw(); legends[1]->Draw(); canvas3->SaveAs(Form("%s.eps", canvas3->GetName())); } void DrawpTCutOff() { /* aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt_ref.root", "Multiplicity", "multiplicityESD.root", kTRUE)' mv unfolded.root chi2_ptref.root aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt0.root", "Multiplicity", "multiplicityESD.root", kTRUE)' mv unfolded.root chi2_pt0.root aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt1.root", "Multiplicity", "multiplicityESD.root", kTRUE)' mv unfolded.root chi2_pt1.root aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt0_25.root", "Multiplicity", "multiplicityESD.root", kTRUE)' mv unfolded.root chi2_pt0_25.root aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt1_25.root", "Multiplicity", "multiplicityESD.root", kTRUE)' mv unfolded.root chi2_pt1_25.root */ loadlibs(); TH1* results[10]; const char* files[] = { "chi2_ptref.root", "chi2_pt0.root", "chi2_pt1.root", "chi2_pt0_25.root", "chi2_pt1_25.root"}; Int_t nMax = 5; for (Int_t i = 0; iGetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i)); } const char* legendStrings[] = { "Reduced 50%", "Enhanced 50%", "Reduced 25%", "Enhanced 25%" }; DrawRatio(results[0], nMax-1, results+1, "LowMomentumSyst.eps", kFALSE, legendStrings); } void ParticleSpeciesComparison() { loadlibs(); TH1* results[10]; TH1* mc = 0; // loop over cases (normal, enhanced/reduced ratios) Int_t nMax = 7; for (Int_t i = 0; iGetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i)); } for (Int_t i=1; i<=results[0]->GetNbinsX(); i++) results[0]->SetBinError(i, 0); const char* legendStrings[] = { "K #times 0.7", "K #times 1.3", "p #times 0.7", "p #times 1.3", "others #times 0.7", "others #times 1.3", }; DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings); } /*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 CompareVertexRecoEfficiencies() { loadlibs(); const char* file1 = "LHC10a12_run10482X"; const char* file2 = "LHC10a14_run10482X_Phojet"; const char* suffix = "/all/spd/multiplicityMC_xsection.root"; hist1 = AliMultiplicityCorrection::Open(Form("%s%s", file1, suffix))->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB); hist2 = AliMultiplicityCorrection::Open(Form("%s%s", file2, suffix))->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB); ratio = (TH1*) hist1->Clone("ratio"); ratio->Divide(hist2); new TCanvas; hist1->Draw(); hist2->SetLineColor(2); hist2->Draw("SAME"); new TCanvas; ratio->Draw(); } void TriggerVertexCorrection() { // // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample // loadlibs(); TFile::Open(correctionFile); AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); mult->LoadHistograms("Multiplicity"); TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL)); TH1* corrNSD = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kNSD)); TH1* corrMB = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB)); TH1* triggerEff = Invert(mult->GetTriggerEfficiency(etaRange, AliMultiplicityCorrection::kNSD)); TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 500); gPad->SetGridx(); gPad->SetGridy(); gPad->SetTopMargin(0.05); gPad->SetRightMargin(0.05); corrINEL->SetStats(kFALSE); corrINEL->GetXaxis()->SetRangeUser(0, 12); corrINEL->GetYaxis()->SetRangeUser(0, 8); corrINEL->SetTitle(Form(";%s;Correction factor", GetMultLabel())); corrINEL->SetMarkerStyle(22); corrINEL->Draw("PE"); corrMB->SetStats(kFALSE); corrMB->SetLineColor(2); corrMB->SetMarkerStyle(25); corrMB->SetMarkerColor(2); corrMB->Draw("SAME PE"); corrNSD->SetLineColor(4); corrNSD->SetMarkerStyle(24); corrNSD->SetMarkerColor(4); corrNSD->Draw("SAME PE"); triggerEff->SetLineColor(6); triggerEff->SetMarkerStyle(25); triggerEff->SetMarkerColor(6); //triggerEff->Draw("SAME PE"); Printf(" MB INEL NSD TRIGINEL"); Printf("bin 0: %.2f %.2f %.2f %.2f", corrMB->GetBinContent(1), corrINEL->GetBinContent(1), corrNSD->GetBinContent(1), triggerEff->GetBinContent(1)); Printf("bin 1: %.2f %.2f %.2f %.2f", corrMB->GetBinContent(2), corrINEL->GetBinContent(2), corrNSD->GetBinContent(2), triggerEff->GetBinContent(2)); TLegend* legend = new TLegend(0.3, 0.6, 0.85, 0.85); legend->SetFillColor(0); legend->AddEntry(corrINEL, "Correction to inelastic sample"); legend->AddEntry(corrNSD, "Correction to NSD sample"); legend->AddEntry(corrMB, "Correction to triggered sample"); legend->SetTextSize(0.04); legend->Draw(); canvas->SaveAs(Form("%s.eps", canvas->GetName())); } void DrawStatisticalUncertainty(const char* fileName = "StatisticalUncertainty.root") { TFile::Open(fileName); errorResponse = (TH1*) gFile->Get("errorResponse"); errorMeasured = (TH1*) gFile->Get("errorMeasured"); errorBoth = (TH1*) gFile->Get("errorBoth"); TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400); canvas->SetGridx(); canvas->SetGridy(); canvas->SetRightMargin(0.05); canvas->SetTopMargin(0.05); errorResponse->SetLineColor(1); errorResponse->GetXaxis()->SetRangeUser(0, longDisplayRange); errorResponse->GetYaxis()->SetRangeUser(0, 1); errorResponse->SetStats(kFALSE); errorResponse->SetTitle(";true multiplicity;Uncertainty"); errorResponse->Draw(); errorMeasured->SetLineColor(2); errorMeasured->Draw("SAME"); errorBoth->SetLineColor(4); errorBoth->Draw("SAME"); errorResponse->Scale(1.0 / sqrt(2)); errorMeasured->Scale(1.0 / sqrt(2)); errorBoth->Scale(1.0 / sqrt(2)); Printf("Average errorResponse: %f", errorResponse->Integral(2, displayRange) / (displayRange - 1)); Printf("Average errorMeasured: %f", errorMeasured->Integral(2, displayRange) / (displayRange - 1)); Printf("Average errorBoth: %f", errorBoth->Integral(2, displayRange) / (displayRange - 1)); canvas->SaveAs(Form("%s.eps", canvas->GetName())); } void StatisticalUncertaintyCompare(const char* det = "SPD") { TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det)); TH1* errorResponse = (TH1*) file1->Get("errorResponse"); TH1* errorMeasured = (TH1*) file1->Get("errorMeasured"); TH1* errorBoth = (TH1*) file1->Get("errorBoth"); TString str; str.Form("StatisticalUncertaintyCompare%s", det); TCanvas* canvas = new TCanvas(str, str, 800, 500); canvas->SetGridx(); canvas->SetGridy(); canvas->SetRightMargin(0.05); canvas->SetTopMargin(0.05); errorResponse->Scale(1.0 / sqrt(2)); errorMeasured->Scale(1.0 / sqrt(2)); errorBoth->Scale(1.0 / sqrt(2)); errorResponse->SetLineColor(1); errorResponse->GetXaxis()->SetRangeUser(0, displayRange); errorResponse->GetYaxis()->SetRangeUser(0, 0.18); errorResponse->SetStats(kFALSE); errorResponse->GetYaxis()->SetTitleOffset(1.2); errorResponse->SetTitle(Form(";%s;#sqrt{2}^{-1} #sigma(unfolded - unfolded_{0}) / unfolded_{0}", GetMultLabel())); errorResponse->Draw(); errorMeasured->SetLineColor(2); errorMeasured->Draw("SAME"); errorBoth->SetLineColor(4); errorBoth->Draw("SAME"); TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det)); TH1* errorResponse2 = (TH1*) file2->Get("errorResponse"); TH1* errorMeasured2 = (TH1*) file2->Get("errorMeasured"); TH1* errorBoth2 = (TH1*) file2->Get("errorBoth"); errorResponse2->Scale(1.0 / sqrt(2)); errorMeasured2->Scale(1.0 / sqrt(2)); errorBoth2->Scale(1.0 / sqrt(2)); errorResponse2->SetLineStyle(2); errorResponse2->Draw("SAME"); errorMeasured2->SetLineColor(2); errorMeasured2->SetLineStyle(2); errorMeasured2->Draw("SAME"); errorBoth2->SetLineColor(4); errorBoth2->SetLineStyle(2); errorBoth2->Draw("SAME"); TLegend* legend = new TLegend(0.2, 0.5, 0.8, 0.9); legend->SetFillColor(0); legend->SetTextSize(0.04); legend->AddEntry(errorBoth, "Both (Bayesian unfolding)"); legend->AddEntry(errorMeasured, "Measured (Bayesian unfolding)"); legend->AddEntry(errorResponse, "Response matrix (Bayesian unfolding)"); legend->AddEntry(errorBoth2, "Both (#chi^{2}-minimization)"); legend->AddEntry(errorMeasured2, "Measured (#chi^{2}-minimization)"); legend->AddEntry(errorResponse2, "Response matrix (#chi^{2}-minimization)"); legend->Draw(); canvas->SaveAs(Form("%s.eps", canvas->GetName())); } void EfficiencyComparison(Int_t eventType = 2, Bool_t uncertainty = kTRUE) { const char* files[] = { "multiplicityMC_nd.root", "multiplicityMC_sd.root", "multiplicityMC_dd.root", "multiplicityMC_xsection.root" }; loadlibs(); TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 600); canvas->SetGridx(); canvas->SetGridy(); canvas->SetRightMargin(0.05); canvas->SetTopMargin(0.05); AliMultiplicityCorrection* data[4]; TH1* effArray[4]; TH1* effErrorArray[2]; Int_t markers[] = { 24, 25, 26, 5 }; //Int_t markers[] = { 2, 25, 24, 5 }; Int_t colors[] = { 1, 2, 4, 6 }; //Int_t colors[] = { 1, 1, 1, 1 }; //TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7); TLegend* legend = new TLegend(0.3, 0.3, 0.9, 0.6); legend->SetTextSize(0.04); legend->SetFillColor(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 = 0; if (eventType == -1) { eff = (TH1*) data[i]->GetTriggerEfficiency(etaRange)->Clone(Form("eff_%d", i)); } else eff = (TH1*) data[i]->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i)); effArray[i] = eff; eff->GetXaxis()->SetRangeUser(0, 15); eff->GetYaxis()->SetRangeUser(0, 1.19); eff->SetStats(kFALSE); eff->GetXaxis()->SetTitle(GetMultLabel()); eff->GetYaxis()->SetTitle("Efficiency"); eff->SetTitle(""); eff->SetLineColor(colors[i]); eff->SetMarkerColor(colors[i]); eff->SetMarkerStyle(markers[i]); if (i == 3) { // once for INEL, once for NSD for (AliMultiplicityCorrection::EventType eventType2 = AliMultiplicityCorrection::kINEL; eventType2 <= AliMultiplicityCorrection::kNSD; eventType2++) { effDiff = (TH1*) data[i]->GetEfficiency(etaRange, eventType2)->Clone(Form("effDiff_%d", i)); for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++) effDiff->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(etaRange, eventType2); for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++) { // TODO we could also do asymmetric errors here Float_t deviation = TMath::Abs(effDiff->GetBinContent(bin) - eff2->GetBinContent(bin)); effDiff->SetBinError(bin, TMath::Max(effDiff->GetBinError(bin), (Double_t) deviation)); } } for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++) { //if (eventType2 == AliMultiplicityCorrection::kINEL) //eff->SetBinError(bin, 0); //eff->SetBinError(bin, effDiff->GetBinError(bin)); if (bin < 20 && effDiff->GetBinContent(bin) > 0) Printf("Bin %d: Error: %.2f", bin, 100.0 * effDiff->GetBinError(bin) / effDiff->GetBinContent(bin)); } if (uncertainty) { TH1* effError = (TH1*) effDiff->Clone(Form("effError_%s", (eventType2 == AliMultiplicityCorrection::kINEL) ? "INEL" : "NSD")); effErrorArray[eventType2 - AliMultiplicityCorrection::kINEL] = effError; effError->Reset(); for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++) if (effDiff->GetBinContent(bin) > 0) effError->SetBinContent(bin, 1.0 * effDiff->GetBinError(bin) / effDiff->GetBinContent(bin)); effError->SetLineColor(1); effError->SetMarkerStyle(1); if (eventType2 == AliMultiplicityCorrection::kNSD) effError->SetLineStyle(2); effError->DrawCopy("SAME HIST"); } } } 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"))))); } if (uncertainty) { legend->AddEntry(effErrorArray[0], "Relative syst. uncertainty: inelastic"); legend->AddEntry(effErrorArray[1], "Relative syst. uncertainty: NSD"); file = TFile::Open("uncertainty_xsection.root", "RECREATE"); effErrorArray[0]->Write(); effErrorArray[1]->Write(); file->Close(); } legend->Draw(); canvas->SaveAs(Form("%s.eps", canvas->GetName())); } void ModelDependencyPlot() { loadlibs(); TFile::Open("multiplicityMC.root"); AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); mult->LoadHistograms("Multiplicity"); hist = mult->GetCorrelation(etaRange); for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y) { for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z) { hist->SetBinContent(0, y, z, 0); hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); } } TH2* proj = (TH2*) hist->Project3D("zy"); TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 1200, 600); canvas->Divide(2, 1); canvas->cd(2); gPad->SetLogy(); gPad->SetGridx(); gPad->SetGridy(); gPad->SetTopMargin(0.05); gPad->SetRightMargin(0.05); Int_t selectedMult = 30; Int_t yMax = 9e4; 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, displayRange); full->GetYaxis()->SetRangeUser(5, yMax); full->SetTitle(";Multiplicity;Entries"); selected->SetLineColor(0); selected->SetMarkerColor(2); selected->SetMarkerStyle(5); 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(); gPad->SetGridx(); gPad->SetGridy(); gPad->SetTopMargin(0.05); gPad->SetRightMargin(0.05); full = proj->ProjectionY("full2"); selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult)); full->SetStats(kFALSE); full->GetXaxis()->SetRangeUser(0, displayRange); full->GetYaxis()->SetRangeUser(5, yMax); full->SetTitle(";Multiplicity;Entries"); full->SetLineColor(0); full->SetMarkerColor(2); full->SetMarkerStyle(5); full->Draw("P"); selected->Draw("SAME"); legend->Draw(); 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"); } void FitPt(const char* fileName) { // needs a MC file from the dNdEta analysis TFile::Open(fileName); TH1* genePt = (TH1*) gFile->Get("fdNdpT"); genePt->SetTitle(";p_{T} (GeV/c);1/p_{T} dN_{ch}/dp_{T} (GeV/c)^{-2}"); // number of events genePt->Scale(1.0 / 287800); // bin width genePt->Scale(1.0 / genePt->GetXaxis()->GetBinWidth(1)); genePt->GetXaxis()->SetRangeUser(0, 0.4); TF1* func = new TF1("func", "[1]*x*exp(x*[0])"); func->SetParameters(-1, 1); genePt->SetMarkerStyle(25); genePt->SetTitle(""); genePt->SetStats(kFALSE); new TCanvas; genePt->DrawCopy("P"); func->DrawCopy("SAME"); gPad->SetLogy(); genePt->Fit(func, "0", "", 0, 0.25); genePt->Fit(func, "0", "", 0, 0.25); TCanvas* canvas = new TCanvas("FitPt", "FitPt", 600, 600); gPad->SetGridx(); gPad->SetGridy(); gPad->SetLeftMargin(0.13); gPad->SetRightMargin(0.05); gPad->SetTopMargin(0.05); //genePt->GetXaxis()->SetRangeUser(0, 1); genePt->GetYaxis()->SetRangeUser(2, 200); genePt->GetYaxis()->SetTitleOffset(1.4); genePt->GetXaxis()->SetTitleOffset(1.1); genePt->DrawCopy("P"); //func->SetRange(0, 0.3); func->DrawCopy("SAME"); gPad->SetLogy(); canvas->SaveAs(Form("%s.eps", canvas->GetName())); 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<2; 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.75 : 1.25); func2->SetParameter(param, func2->GetParameter(param) * factor); //func2->Print(); canvas->cd(); func2->SetLineWidth(2); func2->SetLineColor(2); func2->SetLineStyle(2); func2->DrawCopy("SAME"); canvas2->cd(); TH1* second = func2->GetHistogram(); second->Divide(first); second->SetLineColor(param + 1); // set to 1 above 0.2 GeV for (Int_t bin=second->GetXaxis()->FindBin(0.20001); bin<=second->GetNbinsX(); bin++) second->SetBinContent(bin, 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())); canvas2->SaveAs(Form("%s.eps", canvas2->GetName())); file->Close(); } void DrawSystematicpT() { TFile* file = TFile::Open("SystematicpT.root"); TH1* mcHist2 = (TH1*) file->Get("mymc_unity"); TH1* result2 = (TH1*) file->Get("result_unity"); TH1* mcHist[12]; TH1* result[12]; Int_t nParams = 5; for (Int_t id=0; idGet(Form("mymc_%d_%d.root", id / 2, id % 2)); result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2)); } DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps"); //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps"); DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE); //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps"); // does not make sense: mc is different //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps"); } void SystematicpT(Bool_t chi2 = 1) { gSystem->Load("libPWG0base"); TFile::Open("ptspectrum900.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 = 5; for (Int_t param=0; paramLoadHistograms("Multiplicity"); mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange)); if (chi2) { mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx); } else mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0); 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("mymc_unity"); mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange)); if (chi2) { mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx); } else mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100); TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity"); TFile* file = TFile::Open("SystematicpT.root", "RECREATE"); mcHist2->Write(); result2->Write(); for (Int_t id=0; idWrite(); result[id]->Write(); } file->Close(); DrawSystematicpT(); } void DrawSystematicpT2() { //displayRange = 200; // read from file TFile* file = TFile::Open("SystematicpT2.root"); TH1* mcHist = (TH1*) file->Get("mymc"); TH1* result[12]; result[0] = (TH1*) file->Get("result_unity"); Int_t nParams = 5; for (Int_t id=0; idGet(Form("result_%d_%d", id / 2, id % 2)); DrawResultRatio((TH1*) mcHist->Clone(), (TH1*) result[0]->Clone(), "SystematicpT_OK.eps"); DrawRatio(mcHist, nParams*2+1, result, "SystematicpT_Ratios_MC.eps", kTRUE); DrawRatio(result[0], nParams*2, result+1, "SystematicpT_Ratios.eps"); } void SystematicpT2(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE) { gSystem->Load("libPWG0base"); if (tpc) { SetTPC(); TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root"); } else TFile::Open("ptspectrum100_1.root"); AliMultiplicityCorrection* measured = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); measured->LoadHistograms("Multiplicity"); TH1* mcHist = measured->GetMultiplicityVtx(etaRange)->ProjectionY("mymc"); TH1* result[12]; Int_t nParams = 5; // -1 = unity change, 0...4 parameters for (Int_t id=-1; idLoadHistograms("Multiplicity"); response->SetMultiplicityESD(etaRange, measured->GetMultiplicityESD(etaRange)); if (chi2) { response->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); response->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx); } else response->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0); result[id+1] = (TH1*) response->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%s", idStr.Data())); TString tmp; tmp.Form("SystematicpT_%s.eps", idStr.Data()); DrawResultRatio(mcHist, result[id+1], tmp); } TFile* file = TFile::Open("SystematicpT2.root", "RECREATE"); mcHist->Write(); for (Int_t id=0; idWrite(); file->Close(); DrawSystematicpT2(); } void SystematicpTCutOff(Bool_t chi2 = kTRUE) { // only needed for TPC SetTPC(); gSystem->Load("libPWG0base"); TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root"); AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); mult->LoadHistograms("Multiplicity"); TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root"); AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2"); mult2->LoadHistograms("Multiplicity"); // "normal" result mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange)); if (chi2) { mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx); } else mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100); TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc"); TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1"); TH1* syst[2]; // change of pt spectrum (down) TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root"); AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3"); mult3->LoadHistograms("Multiplicity"); mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange)); if (chi2) { mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx); } else mult3->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100); syst[0] = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2"); // change of pt spectrum (up) TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_inc.root"); AliMultiplicityCorrection* mult4 = new AliMultiplicityCorrection("Multiplicity4", "Multiplicity4"); mult4->LoadHistograms("Multiplicity"); mult4->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange)); if (chi2) { mult4->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); mult4->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx); } else mult4->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100); syst[1] = (TH1*) mult4->GetMultiplicityESDCorrected(etaRange)->Clone("result3"); DrawRatio(result1, 2, syst, "SystematicpTCutOff.eps", kFALSE, 0, kTRUE); Draw2ResultRatio(mcHist, result1, syst[0], "SystematicpTCutOff1.eps"); Draw2ResultRatio(mcHist, result1, syst[1], "SystematicpTCutOff2.eps"); } TH1* SystematicsSummary(Bool_t tpc = 0, Bool_t nsd = kTRUE) { Int_t nEffects = 7; TH1* effects[10]; const char** names = 0; Int_t colors[] = { 1, 2, 4, 1, 2, 4 }; Int_t styles[] = { 1, 2, 3, 1, 2, 3 }; Int_t widths[] = { 1, 1, 1, 2, 2, 2 }; Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26 }; TH1* dummy = new TH2F("dummy", Form(";%s;Uncertainty", GetMultLabel()), 202, -1.5, 200.5, 100, 0, 0.4); dummy->SetStats(0); for (Int_t i=0; iGet("errorBoth"); // smooth a bit, but skip 0 bin effects[0]->SetBinContent(2, hist->GetBinContent(2)); for (Int_t i=3; i<=200; ++i) effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2); // relative x-section effects[1]->SetBinContent(2, 0.005); effects[1]->SetBinContent(3, 0.0025); effects[1]->SetBinContent(4, 0.0025); // particle composition for (Int_t i=2; i<=101; ++i) { if (i < 41) { effects[2]->SetBinContent(i, 0.01); } else if (i < 76) { effects[2]->SetBinContent(i, 0.02); } else effects[2]->SetBinContent(i, 0.02 + 0.08 / 25 * (i - 76)); } // pt cut off (only tpc) for (Int_t i=2; i<=101; ++i) { if (i < 11) { effects[3]->SetBinContent(i, 0.05); } else if (i < 51) { effects[3]->SetBinContent(i, 0.01); } else effects[3]->SetBinContent(i, 0.01 + 0.1 / 30 * (i - 51)); } // track selection (only tpc) for (Int_t i=2; i<=101; ++i) effects[4]->SetBinContent(i, 0.03); // secondaries for (Int_t i=2; i<=101; ++i) effects[5]->SetBinContent(i, 0.01); // pt spectrum for (Int_t i=2; i<=101; ++i) { if (i < 21) { effects[6]->SetBinContent(i, 0.05); } else if (i < 51) { effects[6]->SetBinContent(i, 0.02); } else effects[6]->SetBinContent(i, 0.02 + 0.13 / 25 * (i - 51)); } } else { nEffects = 5; //const char* namesSPD[] = { "Particle composition", "p_{t} cut-off", "Unfolding Method (#chi^{2})", "Relative cross-sections (INEL)", "Relative cross-sections (NSD)"}; const char* namesSPD[] = { "Unfolding Method (#chi^{2})", "Relative cross-sections (INEL)", "Relative cross-sections (NSD)", "Particle composition", "p_{t} cut-off"}; names = namesSPD; currentEffect = 0; // method TFile* file = TFile::Open("StatisticalUncertaintySPDChi2.root"); TH1* hist = (TH1*) file->Get("errorBoth"); hist->Scale(1.0 / sqrt(2)); // smooth a bit, but skip 0 bin /*effects[currentEffect]->SetBinContent(1, hist->GetBinContent(1)); for (Int_t i=2; i<=201; ++i) effects[currentEffect]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);*/ effects[currentEffect] = hist; currentEffect++; // relative x-section file = TFile::Open("uncertainty_xsection.root"); effects[currentEffect++] = (TH1*) file->Get("effError_INEL"); effects[currentEffect] = (TH1*) file->Get("effError_NSD"); effects[currentEffect]->SetLineStyle(1); //effects[2]->SetBinContent(1, 0.20); //effects[2]->SetBinContent(2, 0.01); //effects[2]->SetBinContent(3, 0.002); currentEffect++; // particle composition effects[currentEffect]->SetBinContent(1, 0.16); for (Int_t i=2; i<=81; ++i) { effects[currentEffect]->SetBinContent(i, 0.01 + 0.05 * i / 81); } currentEffect++; // pt spectrum effects[currentEffect]->SetBinContent(1, 0.06); effects[currentEffect]->SetBinContent(2, 0.03); for (Int_t i=3; i<=81; ++i) { if (i <= 61) { effects[currentEffect]->SetBinContent(i, 0.01); } else if (i <= 81) { effects[currentEffect]->SetBinContent(i, 0.01 + 0.05 * (i - 61) / 20); } } // currentEffect++; // // // material budget // for (Int_t i=1; i<=81; ++i) // { // if (i < 5) // effects[currentEffect]->SetBinContent(i, 0.05 - 0.01 * i); // if (i > 51) // effects[currentEffect]->SetBinContent(i, 0.05 * (i - 50) / 30); // } // currentEffect++; } TCanvas* canvas = new TCanvas("SystematicsSummary.eps", "SystematicsSummary.eps", 800, 500); canvas->SetRightMargin(0.05); canvas->SetTopMargin(0.05); //canvas->SetGridx(); canvas->SetGridy(); TLegend* legend = new TLegend(0.2, 0.4, 0.7, 0.4 + 0.5 * nEffects / 7); legend->SetFillColor(0); legend->SetTextSize(0.04); dummy->Draw(); dummy->GetXaxis()->SetRangeUser(0, displayRange); for (Int_t i=0; iClone(Form("current_%d", i)); /*current->Reset(); for (Int_t j=0; jAdd(effects[j]);*/ current->SetLineColor(colors[i]); current->SetLineStyle(styles[i]); current->SetLineWidth(widths[i]); //current->SetFillColor(colors[i]); current->SetMarkerColor(colors[i]); //current->SetMarkerStyle(markers[i]); current->SetStats(kFALSE); current->GetYaxis()->SetRangeUser(0, 0.4); current->DrawCopy("SAME"); legend->AddEntry(current, names[i]); //TLatex* text = new TLatex(displayRange+2, current->GetBinContent(displayRange+1), names[i]); TLatex* text = new TLatex(displayRange+2, 0.1 - i * 0.02, names[i]); text->SetTextSize(0.04); text->SetTextColor(colors[i]); //text->Draw(); } // add total in square TH1* totalINEL = (TH1*) effects[0]->Clone("totalINEL"); totalINEL->Reset(); TH1* totalNSD = (TH1*) totalINEL->Clone("totalNSD"); for (Int_t i=0; iGetBinContent(20)); effects[i]->Multiply(effects[i]); if (i != 2) totalINEL->Add(effects[i]); if (i != 1) totalNSD->Add(effects[i]); } for (Int_t i=1; i<=totalINEL->GetNbinsX(); ++i) { totalINEL->SetBinError(i, 0); if (totalINEL->GetBinContent(i) > 0) totalINEL->SetBinContent(i, TMath::Min(sqrt(totalINEL->GetBinContent(i)), 1.0)); totalNSD->SetBinError(i, 0); if (totalNSD->GetBinContent(i) > 0) totalNSD->SetBinContent(i, TMath::Min(sqrt(totalNSD->GetBinContent(i)), 1.0)); } //Printf("%f", total->GetBinContent(20)); totalINEL->SetMarkerStyle(5); totalINEL->SetMarkerColor(1); legend->AddEntry(totalINEL, "Total (INEL)", "P"); totalNSD->SetMarkerStyle(24); totalNSD->SetMarkerColor(2); legend->AddEntry(totalNSD, "Total (NSD)", "P"); Printf("total in bin 0 is INEL: %f NSD: %f", totalINEL->GetBinContent(1), totalNSD->GetBinContent(1)); totalINEL->DrawCopy("SAME P"); //->SetBinContent(1, 0); totalNSD->DrawCopy("SAME P"); //->SetBinContent(1, 0); legend->Draw(); canvas->SaveAs(canvas->GetName()); file = TFile::Open(Form("%s_syst_error.root", (tpc) ? "tpc" : "spd"), "RECREATE"); totalINEL->Write("inel_1"); totalNSD->Write("nsd_1"); file->Close(); return (nsd) ? totalNSD : totalINEL; } void finalPlot(Bool_t tpc = 0, Bool_t small = kFALSE) { loadlibs(); if (tpc) SetTPC(); //TH1* errorNSD = SystematicsSummary(tpc, 1); TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", 800, 500); canvas->SetRightMargin(0.05); canvas->SetTopMargin(0.05); canvas->SetGridx(); canvas->SetGridy(); legend = new TLegend(0.5, 0.6, 0.9, 0.8); legend->SetFillColor(0); legend->SetTextSize(0.04); for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kINEL; eventType <= AliMultiplicityCorrection::kNSD; eventType++) { AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open((eventType == AliMultiplicityCorrection::kINEL) ? "chi2_inel.root" : "chi2_nsd.root"); TH1* mcHist = mult->GetMultiplicityMC(etaRange, eventType)->ProjectionY("mymc"); TH1* result = mult->GetMultiplicityESDCorrected(etaRange); DrawResultRatio(mcHist, result, Form("finalPlotCheck_%d.eps", (Int_t) eventType)); // normalize result //result->Scale(1.0 / result->Integral(2, displayRange)); result->GetXaxis()->SetRangeUser(0, displayRange); //result->SetBinContent(1, 0); result->SetBinError(1, 0); result->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (etaRange+1) * 0.5)); result->SetMarkerStyle(0); result->SetLineColor(1); result->SetStats(kFALSE); // systematic error TH1* error = SystematicsSummary(tpc, (eventType == AliMultiplicityCorrection::kNSD)); TH1* systError = (TH1*) result->Clone("systError"); for (Int_t i=1; i<=systError->GetNbinsX(); ++i) systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(i)); // change error drawing style systError->SetFillColor(15); if (eventType == AliMultiplicityCorrection::kNSD) { result->SetLineColor(2); result->SetMarkerColor(2); result->SetMarkerStyle(5); } canvas->cd(); systError->DrawCopy(Form("E2 ][ %s", (eventType == AliMultiplicityCorrection::kINEL) ? "" : "SAME")); result->DrawCopy("SAME E ]["); canvas->SetLogy(); legend->AddEntry(result, (eventType == AliMultiplicityCorrection::kINEL) ? "Inelastic cross-section" : "NSD cross-section", (eventType == AliMultiplicityCorrection::kINEL) ? "L" : "P"); } legend->Draw(); /* //TPaveText* text = new TPaveText(10, 1e-3, 50, 1e-4, "B"); TPaveText* text = new TPaveText(0.15, 0.2, 0.5, 0.4, "B NDC"); text->SetFillColor(0); text->SetTextAlign(12); text->AddText("Systematic errors summed quadratically"); text->AddText("0.6 million minimum bias events"); text->AddText("corrected to inelastic events"); text->Draw("B"); TPaveText* text2 = new TPaveText(0.4, 0.7, 0.6, 0.85, "B NDC"); text2->SetFillColor(0); text2->SetTextAlign(12); text2->AddText("#sqrt{s} = 14 TeV"); if (tpc) { text2->AddText("|#eta| < 0.9"); } else text2->AddText("|#eta| < 2.0"); text2->AddText("simulated data (PYTHIA)"); text2->Draw("B"); if (tpc) { TText* text3 = new TText(0.75, 0.6, "TPC - full tracking"); text3->SetNDC(); text3->Draw(); } else { TText* text3 = new TText(0.75, 0.6, "SPD - Tracklets"); text3->SetNDC(); text3->Draw(); } // alice logo TPad* pad = new TPad("pad", "pad", 0.8, 0.7, 0.9, 0.9); pad->Draw(); pad->cd(); TImage* img = TImage::Open("$HOME/alice.png"); img->SetImageQuality(TAttImage::kImgBest); img->Draw(); canvas->Modified(); */ /* TText* text = new TText(10, 1e-4, "Systematic errors summed quadratically"); text->SetTextSize(0.04); text->DrawText(10, 5e-5, "0.6 #cdot 10^{6} minimum bias events"); text->DrawText(10, 3e-5, "TPC tracks in |#eta| < 0.9"); text->DrawText(10, 1e-5, "corrected to ineleastic events in |#eta| < 0.9"); text->Draw();*/ canvas->SaveAs(canvas->GetName()); } void finalPlot2(Bool_t tpc = 0) { loadlibs(); if (tpc) SetTPC(); //TH1* errorNSD = SystematicsSummary(tpc, 1); TCanvas* canvas = new TCanvas("finalPlot2.eps", "finalPlot2.eps", 800, 600); canvas->SetRightMargin(0.05); canvas->SetTopMargin(0.05); canvas->SetGridx(); canvas->SetGridy(); canvas->SetLogy(); legend = new TLegend(0.5, 0.7, 0.9, 0.9); legend->SetFillColor(0); legend->SetTextSize(0.04); Int_t displayRanges[] = { 30, 45, 65 }; TH2* dummy = new TH2F("dummy", ";True multiplicity N_{ch};P(N_{ch})", 100, -0.5, displayRanges[2]+10, 1000, 5e-6, 5); dummy->SetStats(0); dummy->Draw(); TList list; // get syst error file = TFile::Open(Form("%s_syst_error.root", (tpc) ? "tpc" : "spd")); TH1* totalINEL = (TH1*) file->Get("inel_1"); TH1* totalNSD = (TH1*) file->Get("nsd_1"); Int_t count = 0; for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kINEL; eventType <= AliMultiplicityCorrection::kNSD; eventType++) { AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open((eventType == AliMultiplicityCorrection::kINEL) ? "chi2_inel.root" : "chi2_nsd.root"); //TH1* mcHist = mult->GetMultiplicityMC(etaRange, eventType)->ProjectionY("mymc"); for (Int_t etaR = 0; etaR <= 2; etaR++) { TH1* result = mult->GetMultiplicityESDCorrected(etaR); //DrawResultRatio(mcHist, result, Form("finalPlotCheck_%d.eps", (Int_t) eventType)); // normalize result result->Scale(TMath::Power(5, etaR) / result->Integral(1, displayRanges[etaR])); result->GetXaxis()->SetRangeUser(0, displayRanges[etaR]); //result->SetBinContent(1, 0); result->SetBinError(1, 0); //result->SetTitle(Form("")); result->SetMarkerStyle(0); result->SetLineColor(1); result->SetLineWidth(2); //result->SetMarkerStyle(4); //result->SetStats(kFALSE); // systematic error //TH1* error = SystematicsSummary(tpc, (eventType == AliMultiplicityCorrection::kNSD)); TH1* error = (eventType == AliMultiplicityCorrection::kNSD) ? totalNSD : totalINEL; TH1* systError = (TH1*) result->Clone("systError"); // small hack until we have syst. errors for all eta ranges Float_t factor = 80.0 / displayRanges[etaR]; for (Int_t i=1; i<=systError->GetNbinsX(); ++i) { systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(factor * i)); } // change error drawing style systError->SetFillColor(15); if (eventType == AliMultiplicityCorrection::kNSD) { result->SetLineColor(2); result->SetMarkerColor(2); result->SetMarkerStyle(5); } canvas->cd(); systError->GetXaxis()->SetRangeUser(0, displayRanges[etaR]); systError->DrawCopy("E2 ][ SAME"); list.Add(result); if (eventType == AliMultiplicityCorrection::kINEL) { TLatex* Tl = new TLatex; Tl->SetTextSize(0.04); //Tl->SetBit(TLatex::kTextNDC); Tl->SetTextAlign(32); Float_t etaRangeArr[] = { 0.5, 1.0, 1.4 }; TString tmpStr; tmpStr.Form("|#eta| < %.1f (x %d)", etaRangeArr[etaR], (Int_t) TMath::Power(5, etaR)); if (etaR == 0) Tl->DrawLatex(15, result->GetBinContent(20), tmpStr); if (etaR == 1) { Tl->SetTextAlign(12); Tl->DrawLatex(40, result->GetBinContent(40), tmpStr); } if (etaR == 2) { Tl->SetTextAlign(12); Tl->DrawLatex(60, result->GetBinContent(50), tmpStr); } } if (etaR == 0) legend->AddEntry(result, (eventType == AliMultiplicityCorrection::kINEL) ? "Inelastic events" : "NSD events", (eventType == AliMultiplicityCorrection::kINEL) ? "L" : "P"); count++; } } for (Int_t i=0; iDrawCopy("SAME E ]["); legend->Draw(); canvas->SaveAs(canvas->GetName()); } TMatrixD* NonInvertable() { const Int_t kSize = 5; TMatrixD matrix(kSize, kSize); for (Int_t x=0; xClone("unfolded"); TH1F* measuredIdealDist = inputDist->Clone("measuredIdealDist"); measuredIdealDist->SetTitle(";m;Entries"); TH1F* measuredDist = measuredIdealDist->Clone("measuredDist"); TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize); // norm: 1/(sqrt(2pi)sigma) gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8); //gaus->Print(); for (Int_t x=1; x<=inputDist->GetNbinsX(); x++) { Float_t value = gaus->Eval(inputDist->GetBinCenter(x)); inputDist->SetBinContent(x, value); inputDistVector(x-1) = value; } TVectorD measuredDistIdealVector = matrix * inputDistVector; for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++) measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1)); measuredDist->FillRandom(measuredIdealDist, 10000); // fill error matrix before scaling TMatrixD covarianceMatrixMeasured(kSize, kSize); for (Int_t x=1; x<=unfolded->GetNbinsX(); x++) covarianceMatrixMeasured(x-1, x-1) = TMath::Sqrt(measuredDist->GetBinContent(x)); TMatrixD covarianceMatrix = inverted * covarianceMatrixMeasured * inverted; //covarianceMatrix.Print(); TVectorD measuredDistVector(kSize); for (Int_t x=1; x<=measuredDist->GetNbinsX(); x++) measuredDistVector(x-1) = measuredDist->GetBinContent(x); TVectorD unfoldedVector = inverted * measuredDistVector; for (Int_t x=1; x<=unfolded->GetNbinsX(); x++) unfolded->SetBinContent(x, unfoldedVector(x-1)); TCanvas* canvas = new TCanvas("BlobelUnfoldingExample", "BlobelUnfoldingExample", 1200, 600); canvas->SetTopMargin(0.05); canvas->Divide(2, 1); canvas->cd(1); canvas->cd(1)->SetLeftMargin(0.15); canvas->cd(1)->SetRightMargin(0.05); canvas->cd(1)->SetTopMargin(0.05); gPad->SetGridx(); gPad->SetGridy(); measuredDist->GetYaxis()->SetRangeUser(-600, 2799); measuredDist->GetYaxis()->SetTitleOffset(1.9); measuredDist->SetStats(0); measuredDist->DrawCopy(); gaus->Draw("SAME"); canvas->cd(2); canvas->cd(2)->SetLeftMargin(0.15); canvas->cd(2)->SetRightMargin(0.05); canvas->cd(2)->SetTopMargin(0.05); gPad->SetGridx(); gPad->SetGridy(); unfolded->GetYaxis()->SetRangeUser(-600, 2799); unfolded->GetYaxis()->SetTitleOffset(1.9); unfolded->SetStats(0); unfolded->DrawCopy(); gaus->Draw("SAME"); canvas->SaveAs("BlobelUnfoldingExample.eps"); return; // now unfold this with Bayesian loadlibs(); // fill a multiplicity object mult = new AliMultiplicityCorrection("mult", "mult"); for (Int_t x=0; xGetMultiplicityVtx(0)->SetBinContent(1, x+1, inputDistVector(x)); mult->GetMultiplicityESD(0)->SetBinContent(1, x+1, measuredDistVector(x)*10000); for (Int_t y=0; yGetCorrelation(0)->SetBinContent(1, x+1, y+1, matrix(x, y)); } //mult->DrawHistograms(); mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, 0); //mult->SetCreateBigBin(kFALSE); mult->ApplyMinuitFit(0, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE); //hist2->ProjectionY("mymchist")); //mult->ApplyBayesianMethod(0, kFALSE, AliMultiplicityCorrection::kTrVtx, 0, -1, 0, kFALSE); mult->DrawComparison("BlobelExample", 0, kFALSE, kTRUE, mult->GetMultiplicityVtx(0)->ProjectionY("mcmchist", 1, mult->GetMultiplicityVtx(0)->GetNbinsX())); } void TestWithDiagonalMatrix() { const Int_t kSize = 20; TMatrixD matrix(kSize, kSize); for (Int_t x=0; xClone("measuredIdealDist"); measuredIdealDist->SetTitle(";m;Entries"); TH1F* measuredDist = measuredIdealDist->Clone("measuredDist"); TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize); // norm: 1/(sqrt(2pi)sigma) gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8); //gaus->Print(); for (Int_t x=1; x<=inputDist->GetNbinsX(); x++) { Float_t value = gaus->Eval(inputDist->GetBinCenter(x)); inputDist->SetBinContent(x, value); inputDistVector(x-1) = value; } TVectorD measuredDistIdealVector = matrix * inputDistVector; for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++) measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1)); // randomize //measuredDist->FillRandom(measuredIdealDist, 10000); measuredDist = measuredIdealDist; measuredDist->Sumw2(); for (Int_t x=0; xGetBinContent(x+1), measuredDist->GetBinError(x+1)); // now unfold this loadlibs(); // fill a multiplicity object mult = new AliMultiplicityCorrection("mult", "mult"); for (Int_t x=0; xGetMultiplicityVtx(0)->SetBinContent(1, x+1, inputDistVector(x)); mult->GetMultiplicityESD(0)->SetBinContent(1, x+1, measuredDist->GetBinContent(x+1)); for (Int_t y=0; yGetCorrelation(0)->SetBinContent(1, x+1, y+1, matrix(x, y)); } //mult->DrawHistograms(); AliUnfolding::SetNbins(20, 20); AliUnfolding::SetChi2Regularization(AliUnfolding::kNone, 0); AliUnfolding::SetChi2Regularization(AliUnfolding::kPol1, 1e-14); //mult->SetCreateBigBin(kFALSE); mult->ApplyMinuitFit(0, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE); //hist2->ProjectionY("mymchist")); //mult->ApplyBayesianMethod(0, kFALSE, AliMultiplicityCorrection::kTrVtx, 0, -1, 0, kFALSE); mult->DrawComparison("TestWithDiagonalMatrix", 0, kFALSE, kTRUE, mult->GetMultiplicityVtx(0)->ProjectionY("mcmchist", 1, mult->GetMultiplicityVtx(0)->GetNbinsX())); } void E735Fit() { TH1* fCurrentESD = new TH1F("mult", "mult", 501, -0.5, 500.5); fCurrentESD->Sumw2(); // Open the input stream ifstream in; in.open("e735data.txt"); while(in.good()) { Float_t x, y, ye; in >> x >> y >> ye; //Printf("%f %f %f", x, y, ye); fCurrentESD->SetBinContent(fCurrentESD->FindBin(x), y); fCurrentESD->SetBinError(fCurrentESD->FindBin(x), ye); } in.close(); //new TCanvas; fCurrentESD->DrawCopy(); gPad->SetLogy(); 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->SetMarkerStyle(0); fCurrentESD->SetLineColor(1); fCurrentESD->GetYaxis()->SetTitleOffset(1.3); fCurrentESD->SetTitle(";true multiplicity (N);P_{N}"); fCurrentESD->Draw(""); fCurrentESD->GetXaxis()->SetRangeUser(0, 250); fCurrentESD->Fit(func, "0", "", 0, 150); func->SetRange(0, 250); func->Draw("SAME"); printf("chi2 = %f\n", func->GetChisquare()); fCurrentESD->Fit(lognormal, "0", "", 0.01, 150); lognormal->SetLineColor(2); lognormal->SetLineStyle(2); lognormal->SetRange(0, 250); lognormal->Draw("SAME"); gPad->SetLogy(); canvas->SaveAs("E735Fit.eps"); } void DifferentSamples() { loadlibs(); Int_t n = 2; const char* filesChi2[] = { "chi2_100k_1.root", "chi2_100k_2.root" }; const char* filesBayesian[] = { "bayesian_100k_1.root", "bayesian_100k_2.root" }; TCanvas* canvas = new TCanvas("DifferentSamples", "DifferentSamples", 1200, 600); canvas->Divide(2, 1); legend = new TLegend(0.15, 0.7, 0.65, 0.9); legend->SetFillColor(0); legend->SetTextSize(0.04); for (Int_t i=0; iLoadHistograms("Multiplicity"); AliMultiplicityCorrection* bayesian = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); TFile::Open(filesBayesian[i]); bayesian->LoadHistograms("Multiplicity"); chi2Hist = chi2->GetMultiplicityESDCorrected(etaRange); bayesianHist = bayesian->GetMultiplicityESDCorrected(etaRange); mc = chi2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1); // normalize and divide chi2Hist->Scale(1.0 / chi2Hist->Integral(2, displayRange+1) * mc->Integral(2, displayRange)); bayesianHist->Scale(1.0 / bayesianHist->Integral(2, displayRange+1) * mc->Integral(2, displayRange)); chi2Hist->Divide(mc, chi2Hist); bayesianHist->Divide(mc, bayesianHist); canvas->cd(i+1); gPad->SetTopMargin(0.05); gPad->SetRightMargin(0.05); //gPad->SetLeftMargin(0.12); gPad->SetGridx(); gPad->SetGridy(); chi2Hist->GetXaxis()->SetRangeUser(0, displayRange); chi2Hist->GetYaxis()->SetTitleOffset(1.3); chi2Hist->SetStats(0); chi2Hist->SetTitle(Form(";%s;MC / unfolded", GetMultLabel())); chi2Hist->GetYaxis()->SetRangeUser(0.2, 1.8); chi2Hist->Draw("HIST"); for (Int_t x=1; x<=bayesianHist->GetNbinsX(); x++) bayesianHist->SetBinError(x, 1e-6); bayesianHist->SetLineColor(2); bayesianHist->SetMarkerColor(2); bayesianHist->SetMarkerStyle(5); bayesianHist->Draw("HIST E SAME"); if (i == 0) { legend->AddEntry(chi2Hist, "#chi^{2}-minimization", "L"); legend->AddEntry(bayesianHist, "Bayesian unfolding", "LP"); } legend->Draw(); } canvas->SaveAs("DifferentSamples.eps"); } void PileUp() { loadlibs(); AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicityMC.root"); hist2d = mult->GetMultiplicityMC(etaRange, AliMultiplicityCorrection::kINEL); mult1 = hist2d->ProjectionY("mult1", 1, hist2d->GetNbinsX()); conv = (TH1*) mult1->Clone("conv"); conv->Reset(); mult1->Scale(1.0 / mult1->Integral()); for (Int_t i=1; i<=mult1->GetNbinsX(); i++) for (Int_t j=1; j<=mult1->GetNbinsX(); j++) conv->Fill(mult1->GetBinCenter(i)+mult1->GetBinCenter(j), mult1->GetBinContent(i) * mult1->GetBinContent(j)); conv->Scale(1.0 / conv->Integral()); c = new TCanvas("c", "c", 800, 500); gPad->SetLogy(); gPad->SetTopMargin(0.05); gPad->SetRightMargin(0.05); mult1->SetTitle(Form(";%s;Probability", GetMultLabel())); mult1->SetStats(0); gPad->SetGridx(); gPad->SetGridy(); mult1->Draw(); mult1->GetYaxis()->SetRangeUser(1e-7, 2 * mult1->GetMaximum()); mult1->GetXaxis()->SetRangeUser(0, displayRange); mult1->GetXaxis()->SetTitleOffset(1.15); conv->SetLineColor(2); conv->SetMarkerColor(2); conv->SetMarkerStyle(5); conv->DrawCopy("SAME P"); conv->Scale(0.00058); conv->DrawCopy("SAME P"); legend = new TLegend(0.73, 0.73, 0.93, 0.93); legend->SetFillColor(0); legend->SetTextSize(0.04); legend->AddEntry(mult1, "1 collision"); legend->AddEntry(conv, "2 collisions", "P"); legend->Draw(); c->SaveAs("pileup.eps"); new TCanvas; conv->Divide(mult1); conv->Draw(); } void TestErrorDetermination(Int_t nRandomizations) { TF1* func = new TF1("nbd", "exp(log([0]) + TMath::LnGamma([2]+x) - TMath::LnGamma([2]) - TMath::LnGamma(x+1) + log([1] / ([1]+[2])) * x + log(1.0 + [1]/[2]) * -[2])", 0, 100); func->SetParNames("scaling", "averagen", "k"); func->SetParameters(1, 15, 2); TF1* func2 = new TF1("nbd2", "exp(log([0]) + TMath::LnGamma([2]+x) - TMath::LnGamma([2]) - TMath::LnGamma(x+1) + log([1] / ([1]+[2])) * x + log(1.0 + [1]/[2]) * -[2])", 0, 100); func2->SetParNames("scaling", "averagen", "k"); func2->SetParLimits(0, 0.5, 2); func2->SetParLimits(1, 1, 50); func2->SetParLimits(2, 1, 10); func2->SetParameters(1, 15, 2); //func2->FixParameter(0, 1); //new TCanvas; func->Draw("L"); hist1 = new TH1F("hist1", "", 100, 0.5, 100.5); hist2 = new TH1F("hist2", "", 100, 0.5, 100.5); hist1->Sumw2(); TH1* params[3]; params[0] = new TH1F("param_0", Form("param_%d", 0), 100, 0.95, 1.05); params[1] = new TH1F("param_1", Form("param_%d", 1), 100, 14, 16); params[2] = new TH1F("param_2", Form("param_%d", 2), 100, 1.8, 2.2); const Int_t nTries = 1000; for (Int_t i=0; iReset(); if (nRandomizations == 1) { hist1->FillRandom("nbd", 10000); } else if (nRandomizations == 2) { hist2->Reset(); hist2->FillRandom("nbd", 10000); hist1->FillRandom(hist2, 10000); } else if (nRandomizations == 3) { hist2->Reset(); hist1->FillRandom("nbd", 10000); hist2->FillRandom(hist1, 10000); hist1->Reset(); hist1->FillRandom(hist2, 10000); } else return; //new TCanvas; hist1->Draw(); hist1->Scale(1.0 / hist1->Integral()); hist1->Fit(func2, "NQ"); hist1->Fit(func2, "NQ"); for (Int_t j=0; j<3; j++) params[j]->Fill(func2->GetParameter(j)); } for (Int_t j=0; j<3; j++) { new TCanvas; params[j]->Draw(); params[j]->Fit("gaus"); Printf("sigma of param %d if %f", j, ((TF1*) params[j]->FindObject("gaus"))->GetParameter(2)); } } void DrawRawDistributions(const char* fileName = "multiplicityESD.root") { loadlibs(); AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileName); c = new TCanvas("c", "c", 600, 600); dummy = new TH2F("dummy", ";measured multiplicity", 100, -0.5, 149.5, 100, 0.5, 4e4); dummy->SetStats(0); dummy->Draw(); gPad->SetGridx(); gPad->SetGridy(); gPad->SetLogy(); Int_t colors[] = { 1, 2, 4 }; for (Int_t i=2; i>=0; i--) { hist = mult->GetMultiplicityESD(i)->ProjectionY(); hist->SetLineColor(colors[i]); hist->DrawCopy("SAME"); } } void FindUnfoldedLimit() { loadlibs(); AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicityMC.root"); AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open("multiplicityESD.root"); TH1* hist = mult->GetCorrelation(etaRange); for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y) { for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z) { hist->SetBinContent(0, y, z, 0); hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); } } TH2* corr = (TH2*) ((TH3*) mult->GetCorrelation(etaRange))->Project3D("zy"); TH1* esd_proj = esd->GetMultiplicityESD(etaRange)->ProjectionY("esd_proj", 1, esd->GetMultiplicityESD(etaRange)->GetNbinsX()); //esd_proj = corr->ProjectionY("esd_proj", 1, corr->GetNbinsX()); new TCanvas; corr->Draw("COLZ"); new TCanvas; esd_proj->DrawCopy(); TH1* percentage = (TH1*) (esd_proj->Clone("percentage")); percentage->SetTitle("percentage"); percentage->Reset(); for (Int_t i=1; i<=esd_proj->GetNbinsX(); i++) if (esd_proj->GetBinContent(i) > 0) esd_proj->SetBinContent(i, 1); Int_t limit = -1; for (Int_t i=1; i<=percentage->GetNbinsX(); i++) { TH1* binResponse = corr->ProjectionY("proj", i, i); if (binResponse->Integral() <= 0) continue; binResponse->Scale(1.0 / binResponse->Integral()); binResponse->Multiply(esd_proj); //new TCanvas; binResponse->Draw(); Float_t value = binResponse->Integral(); percentage->SetBinContent(i, value); if (limit == -1 && value < 0.9) limit = percentage->GetXaxis()->GetBinCenter(i-1); //return; } Printf("Limit is %d", limit); new TCanvas; percentage->Draw(); new TCanvas; mc = esd->GetMultiplicityVtx(etaRange)->ProjectionY("mc", 1, esd->GetMultiplicityVtx(etaRange)->GetNbinsX()); mc->SetLineColor(2); mc->Draw(""); } void FindUnfoldedLimitAll() { loadlibs(); for (etaRange = 0; etaRange <= 2; etaRange++) FindUnfoldedLimit(); } void CompareUnfoldedWithMC() { loadlibs(); AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("unfolded.root"); //mult->GetMultiplicityESDCorrected(etaRange)->Scale(1.0 / mult->GetMultiplicityESDCorrected(etaRange)->Integral()); mult->GetMultiplicityESDCorrected(etaRange)->SetTitle(";multiplicity;events"); mult->GetMultiplicityESDCorrected(etaRange)->SetStats(0); mult->GetMultiplicityESDCorrected(etaRange)->Draw(); //mcHist->Scale(1.0 / mcHist->Integral()); mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1); mcHist->SetLineColor(2); mcHist->Draw("SAME"); gPad->SetLogy(); } void PaperNumbers() { const char* label[] = { "SD", "DD", "ND" }; const char* files[] = { "multiplicity_syst_sd.root", "multiplicity_syst_dd.root", "multiplicity_syst_nd.root" }; loadlibs(); Printf("vertex reco"); Float_t oneParticle[3]; for (Int_t i=0; i<3; i++) { AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(files[i]); eff = mult->GetEfficiency(2, AliMultiplicityCorrection::kMB); //eff = mult->GetTriggerEfficiency(2); oneParticle[i] = 100.0 * eff->GetBinContent(2); Printf("%s: %.2f", label[i], oneParticle[i]); } Float_t ua5_SD = 0.153; Float_t ua5_DD = 0.080; Float_t ua5_ND = 0.767; Float_t vtxINELUA5 = ua5_SD * oneParticle[0] + ua5_DD * oneParticle[1] + ua5_ND * oneParticle[2]; Float_t vtxNSDUA5 = (ua5_DD * oneParticle[1] + ua5_ND * oneParticle[2]) / (ua5_DD + ua5_ND); Printf("INEL (UA5) = %.1f", vtxINELUA5); Printf("NSD (UA5) = %.1f", vtxNSDUA5); Printf("total for >= 2 charged tracks"); AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicity.root"); vtx = mult->GetMultiplicityVtx(2)->ProjectionY("vtx", 1, mult->GetMultiplicityVtx(2)->GetNbinsX(), "e"); all = mult->GetMultiplicityINEL(2)->ProjectionY("all", 1, mult->GetMultiplicityINEL(2)->GetNbinsX(), "e"); Int_t begin = vtx->GetXaxis()->FindBin(2); Int_t end = vtx->GetNbinsX(); Float_t above2 = vtx->Integral(begin, end) / all->Integral(begin, end); Printf("%.1f", 100.0 * above2); } void DrawRawEta() { TFile::Open(measuredFile); Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10}; for (Int_t i=0; i<3; i++) { hist = dynamic_cast (gFile->Get(Form("fEta_%d", i))); hist->GetYaxis()->SetRangeUser(0, 8); hist->SetLineColor(colors[i]); hist->SetStats(kFALSE); hist->Draw((i == 0) ? "HIST" : "HIST SAME"); } gPad->SetGridx(); gPad->SetGridy(); } void CrossSectionUncertainties(Int_t xsectionID, Int_t eventType) { const char* files[] = { "multiplicitySD.root", "multiplicityDD.root", "multiplicityND.root" }; loadlibs(); AliMultiplicityCorrection* data[3]; Float_t ref_SD = -1; Float_t ref_DD = -1; Float_t ref_ND = -1; Float_t error_SD = -1; Float_t error_DD = -1; Float_t error_ND = -1; gROOT->ProcessLine(gSystem->ExpandPathName(".L $ALICE_ROOT/PWG0/dNdEta/drawSystematics.C")); GetRelativeFractions(xsectionID, ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND); TH1* unc[9*3]; TH1* unc2[9*3]; const char* legendStrings[8]; for (Int_t i=0; i<9; i++) { Float_t factorSD = 0; Float_t factorDD = 0; if (i > 0 && i < 3) factorSD = (i % 2 == 0) ? 1 : -1; else if (i >= 3 && i < 5) factorDD = (i % 2 == 0) ? 1 : -1; else if (i >= 5 && i < 9) { factorSD = ((i % 2 == 0) ? 1.0 : -1.0) / TMath::Sqrt(2); if (i == 5 || i == 6) factorDD = factorSD; else factorDD = -factorSD; } Float_t scalesSD = ref_SD + factorSD * error_SD; Float_t scalesDD = ref_DD + factorDD * error_DD; Float_t scalesND = 1.0 - scalesDD[i] - scalesSD[i]; str = new TString; str->Form("Case %d: SD: %.2f DD: %.2f ND: %.2f", i, scalesSD, scalesDD, scalesND); Printf("%s", str->Data()); if (i > 0) legendStrings[i-1] = str->Data(); for (Int_t j=0; j<3; ++j) { TString name; name.Form("Multiplicity_%d", j); TFile::Open(files[j]); data[j] = new AliMultiplicityCorrection(name, name); data[j]->LoadHistograms("Multiplicity"); } // calculating relative Float_t sd = data[0]->GetMultiplicityINEL(0)->Integral(0, data[0]->GetMultiplicityINEL(0)->GetNbinsX()+1); Float_t dd = data[1]->GetMultiplicityINEL(0)->Integral(0, data[1]->GetMultiplicityINEL(0)->GetNbinsX()+1); Float_t nd = data[2]->GetMultiplicityINEL(0)->Integral(0, data[2]->GetMultiplicityINEL(0)->GetNbinsX()+1); Float_t total = nd + dd + sd; nd /= total; sd /= total; dd /= total; Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd); Float_t ratio[3]; ratio[0] = scalesSD / sd; ratio[1] = scalesDD / dd; ratio[2] = scalesND / nd; Printf("SD=%.2f, DD=%.2f, ND=%.2f", ratio[0], ratio[1], ratio[2]); TList list; for (Int_t k=0; k<3; ++k) { // modify x-section for (Int_t j=0; jGetMultiplicityVtx(j)->Scale(ratio[k]); data[k]->GetMultiplicityMB(j)->Scale(ratio[k]); data[k]->GetMultiplicityINEL(j)->Scale(ratio[k]); data[k]->GetMultiplicityNSD(j)->Scale(ratio[k]); } if (k > 0) list.Add(data[k]); } printf("Case %d, %s: Entries in response matrix: SD: %.2f DD: %.2f ND: %.2f", xsectionID, data[0]->GetName(), data[0]->GetCorrelation(0)->Integral(), data[1]->GetCorrelation(0)->Integral(), data[2]->GetCorrelation(0)->Integral()); data[0]->Merge(&list); data[0]->FixTriggerEfficiencies(20); Printf(" Total: %.2f", data[0]->GetCorrelation(0)->Integral()); for (Int_t etaR = 0; etaR < 3; etaR++) { unc[i+(etaR*9)] = (TH1*) data[0]->GetTriggerEfficiency(etaR, eventType)->Clone(Form("eff_%d_%d", etaR, i)); unc2[i+(etaR*9)] = (TH1*) data[0]->GetEfficiency(etaR, eventType)->Clone(Form("eff_%d_%d", etaR, i)); } } file = TFile::Open(Form("systunc_fractions_%s.root", (eventType == 2) ? "inel" : "nsd"), "RECREATE"); for (Int_t etaR = 0; etaR < 3; etaR++) { largestError1 = (TH1*) unc[0]->Clone(Form("fractions_trig_%d", etaR)); largestError1->Reset(); largestError2 = (TH1*) unc[0]->Clone(Form("fractions_trigvtx_%d", etaR)); largestError2->Reset(); etaRange = etaR; // correct labels in DrawRatio DrawRatio(unc[etaR*9], 8, unc+1+etaR*9, Form("xsections_trig%d.png", etaR), kFALSE, legendStrings, kFALSE, largestError1); DrawRatio(unc2[etaR*9], 8, unc2+1+etaR*9, Form("xsections_trigvtx%d.png", etaR), kFALSE, legendStrings, kFALSE, largestError2); largestError2->SetBinContent(1, largestError1->GetBinContent(1)); largestError2->GetYaxis()->UnZoom(); largestError2->Write(); new TCanvas; largestError2->DrawCopy(); } file->Close(); } void PythiaPhojetUncertainties() { PythiaPhojetUncertainties(900, 2); PythiaPhojetUncertainties(900, 3); PythiaPhojetUncertainties(2360, 2); PythiaPhojetUncertainties(2360, 3); } void PythiaPhojetUncertainties(Int_t energy, Int_t eventType) { loadlibs(); AliMultiplicityCorrection* mult[2]; const char* trigger = "all"; if (eventType == 3) trigger = "v0and"; if (energy == 2360) trigger = "spdgfobits"; if (energy == 7000) trigger = "all-onepart"; if (energy == 900) TFile::Open(Form("LHC10a12_run10482X/%s/spd/multiplicityMC_xsection.root", trigger)); else if (energy == 2360) TFile::Open(Form("LHC10a10_run105054_7/%s/spd/multiplicityMC_xsection.root", trigger)); else TFile::Open(Form("LHC10b2_run114931/%s/spd/multiplicity.root", trigger)); mult[0] = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); mult[0]->LoadHistograms("Multiplicity"); if (energy == 900) TFile::Open(Form("LHC10a14_run10482X_Phojet/%s/spd/multiplicityMC_xsection.root", trigger)); else if (energy == 2360) TFile::Open(Form("LHC10a15_run10505X_Phojet/%s/spd/multiplicityMC_xsection.root", trigger)); else TFile::Open(Form("LHC10b1_run114931/%s/spd/multiplicity.root", trigger)); mult[1] = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); mult[1]->LoadHistograms("Multiplicity"); TH1* unc[6]; TH1* unc2[6]; for (Int_t i=0; i<2; i++) { mult[i]->FixTriggerEfficiencies(20); for (Int_t etaR = 0; etaR < 3; etaR++) { unc[i+(etaR*2)] = (TH1*) mult[i]->GetTriggerEfficiency(etaR, eventType)->Clone(Form("eff_%d_%d", etaR, i)); unc2[i+(etaR*2)] = (TH1*) mult[i]->GetEfficiency(etaR, eventType)->Clone(Form("eff_%d_%d", etaR, i)); } } file = TFile::Open(Form("systunc_pythiaphojet_%d_%s.root", energy, (eventType == 2) ? "inel" : "nsd"), "RECREATE"); for (Int_t etaR = 0; etaR < 3; etaR++) { largestError1Low = (TH1*) unc[0]->Clone(Form("pythiaphojet_trig_%d_low", etaR)); largestError1Low->Reset(); largestError1High = (TH1*) largestError1Low->Clone(Form("pythiaphojet_trig_%d_high", etaR)); largestError2Low = (TH1*) largestError1Low->Clone(Form("pythiaphojet_trigvtx_%d_low", etaR)); largestError2High = (TH1*) largestError1Low->Clone(Form("pythiaphojet_trigvtx_%d_high", etaR)); DrawRatio(unc[etaR*2], 1, unc+1+etaR*2, Form("pythiaphojet_trig%d.png", etaR), kFALSE, 0, kFALSE, largestError1Low, largestError1High); DrawRatio(unc2[etaR*2], 1, unc2+1+etaR*2, Form("pythiaphojet_trigvtx%d.png", etaR), kFALSE, 0, kFALSE, largestError2Low, largestError2High); largestError2Low->SetBinContent(1, largestError1Low->GetBinContent(1)); largestError2High->SetBinContent(1, largestError1High->GetBinContent(1)); largestError2Low->GetYaxis()->UnZoom(); largestError2High->GetYaxis()->UnZoom(); largestError2Low->Write(); largestError2High->Write(); new TCanvas; largestError2Low->DrawCopy(); new TCanvas; largestError2High->DrawCopy(); } file->Close(); } void CombinatoricsAsFunctionOfN(Int_t etaR) { loadlibs(); mult = AliMultiplicityCorrection::Open("multiplicityESD.root"); //Float_t integratedCombinatorics[] = { 0.00062, 0.00074, 0.001 }; // 900 GeV Float_t integratedCombinatorics[] = { 7.5e-4, 1.1e-3, 1.2e-3 }; // 2.36 TeV multHist = mult->GetMultiplicityESD(etaR)->ProjectionY(); multHist->Scale(1.0 / multHist->Integral()); Float_t integral = 0; Float_t tracks = 0; for (Int_t i=1; i<= multHist->GetNbinsX(); i++) { Float_t x = multHist->GetXaxis()->GetBinCenter(i); integral += x*x * multHist->GetBinContent(i+1); tracks += x * multHist->GetBinContent(i+1); Printf("%f %f %f", x, multHist->GetBinContent(i+1), integral); } Printf("Integral is %f; %f", integral, integral / tracks); Printf("alpha is %e", integratedCombinatorics[etaR] / (integral / tracks)); } void TryCombinatoricsCorrection(Float_t etaR, Float_t alpha) { loadlibs(); mult = AliMultiplicityCorrection::Open("multiplicity.root"); multHist = mult->GetMultiplicityESD(etaR)->ProjectionY(); target = (TH1*) multHist->Clone("target"); target->Reset(); Int_t totalTracks1 = 0; Int_t totalTracks2 = 0; for (Int_t i=1; i<=multHist->GetNbinsX(); i++) { for (Int_t j=0; jGetBinContent(i); j++) { Int_t origTracks = (Int_t) multHist->GetXaxis()->GetBinCenter(i); tracks = origTracks; //tracks -= gRandom->Poisson(1.3e-4 * origTracks * origTracks); if (gRandom->Uniform() < alpha * origTracks * origTracks) tracks--; target->Fill(tracks); totalTracks1 += origTracks; totalTracks2 += tracks; } } new TCanvas; multHist->Draw(); target->Draw("SAME"); target->SetLineColor(2); Printf("%d %d %f %f", totalTracks1, totalTracks2, 1. - (Float_t) totalTracks2 / totalTracks1, 1. - target->GetMean() / multHist->GetMean()); } void ApplyCombinatoricsCorrection() { loadlibs(); mult = AliMultiplicityCorrection::Open("multiplicity.root"); fractionMC = new TF1("fractionMC", "[0]+[1]*x", 0, 200); fractionMC->SetParameters(0.0004832, 5.82e-5); // mariella's presentation 16.04.10 fractionData = new TF1("fractionData", "[0]+[1]*x", 0, 200); fractionData->SetParameters(0.00052, 0.0001241); // mariella's presentation 16.04.10 Int_t etaR = 1; esd = mult->GetMultiplicityESD(etaR); for (Int_t x=1; x<=esd->GetNbinsX(); x++) { for (Int_t y=2; y<=esd->GetNbinsY(); y++) { printf("Before: %f %f ", esd->GetBinContent(x, y-1), esd->GetBinContent(x, y)); Float_t n = esd->GetYaxis()->GetBinCenter(y); Int_t addComb = (fractionData->Eval(n) - fractionMC->Eval(n)) * n * esd->GetBinContent(x, y); // shift addComb down one bin esd->SetBinContent(x, y-1, esd->GetBinContent(x, y-1) + addComb); esd->SetBinContent(x, y, esd->GetBinContent(x, y) - addComb); // recalc errors esd->SetBinError(x, y-1, TMath::Sqrt(esd->GetBinContent(x, y-1))); esd->SetBinError(x, y, TMath::Sqrt(esd->GetBinContent(x, y))); Printf("comb: %d; after: %f +- %f %f +- %f", addComb, esd->GetBinContent(x, y-1), esd->GetBinError(x, y-1), esd->GetBinContent(x, y), esd->GetBinError(x, y)); } } TFile::Open("out.root", "RECREATE"); mult->SaveHistograms(); }