analysis_INTT/plot/ratio_dNdEta.C

0001 #include "./util_combine.h"
0002
0003 std::vector<float> v_sphenix_centralitybin = {0, 3, 6, 10, 15, 20, 25, 30, 35, 40, 45, 50};
0004 int nsphenix = v_sphenix_centralitybin.size() - 1;
0005 float rho = 2537.15 / (sqrt(4476.51) * sqrt(1974.34));
0006 bool drawall = false;
0007 int cent_low = 20, cent_high = 25;
0008
0009 float LeftMargin = 0.17;
0010 float RightMargin = 0.05;
0011 float TopMargin = 0.08;
0012 float BottomMargin = 0.15;
0013
0014 void ratio_dNdEta()
0015 {
0016     // gStyle->SetPalette(kLightTemperature);
0017
0018     std::string plotdir = "./dNdEtaFinal";
0019     system(Form("mkdir -p %s", plotdir.c_str()));
0020
0021     std::vector<TH1 *> v_sphenix_dNdeta_ratio; // ratio of the two approaches
0022     for (int i = 0; i < nsphenix; i++)
0023     {
0024         if (!drawall && !(v_sphenix_centralitybin[i] == cent_low && v_sphenix_centralitybin[i + 1] == cent_high))
0025             continue;
0026
0027         int Mbin = GetMbinNum(Form("Centrality%dto%d", (int)v_sphenix_centralitybin[i], (int)v_sphenix_centralitybin[i + 1]));
0028
0029         // CMS-inspired approach
0030         TFile *f = new TFile(Form("./systematics/Centrality%dto%d_Zvtxm10p0to10p0_noasel/finalhists_systematics_Centrality%dto%d_Zvtxm10p0to10p0_noasel.root", (int)v_sphenix_centralitybin[i], (int)v_sphenix_centralitybin[i + 1], (int)v_sphenix_centralitybin[i], (int)v_sphenix_centralitybin[i + 1]), "READ");
0031         TH1D *h = (TH1D *)f->Get("hM_final");
0032         h->SetDirectory(0);
0033         f->Close();
0034
0035         // PHOBOS-inspired approach
0036         TFile *fphobos = new TFile(Form("/sphenix/user/ChengWei/sPH_dNdeta/Run24AuAuMC/Sim_HIJING_MDC2_ana472_20250307/Run7/EvtVtxZ/FinalResult_10cm_Pol2BkgFit_DeltaPhi0p026/completed/vtxZ_-10_10cm_MBin%d/Final_Mbin%d_00054280/Final_Mbin%d_00054280.root", Mbin, Mbin, Mbin), "READ");
0037         TH1D *hphobos = (TH1D *)fphobos->Get("h1D_dNdEta_reco");
0038         hphobos->SetDirectory(0);
0039         fphobos->Close();
0040
0041         // ratio
0042         TH1 *hratio = (TH1 *)h->Clone(Form("hratio_centraliy%dto%d", (int)v_sphenix_centralitybin[i], (int)v_sphenix_centralitybin[i + 1]));
0043         // loop over bins between |eta|<=1.1
0044         for (int j = 1; j <= hratio->GetNbinsX(); j++)
0045         {
0046             if (fabs(hratio->GetBinCenter(j)) > 1.1)
0047                 continue;
0048
0049             float dNdeta_cms = h->GetBinContent(j);
0050             float dNdeta_phobos = hphobos->GetBinContent(j);
0051             float dNdeta_ratio = dNdeta_cms / dNdeta_phobos;
0052             // error propagation for the ratio taking into account the correlation between the two approaches
0053             float dNdeta_cms_err = h->GetBinError(j);
0054             float dNdeta_phobos_err = hphobos->GetBinError(j);
0055             float dNdeta_ratio_err = dNdeta_ratio * sqrt(pow(dNdeta_cms_err / dNdeta_cms, 2) + pow(dNdeta_phobos_err / dNdeta_phobos, 2) - 2 * rho * dNdeta_cms_err * dNdeta_phobos_err / (dNdeta_cms * dNdeta_phobos));
0056             std::cout << "rho = " << rho << std::endl;
0057             std::cout << "dNdeta_cms = " << dNdeta_cms << " +/- " << dNdeta_cms_err << std::endl;
0058             std::cout << "dNdeta_phobos = " << dNdeta_phobos << " +/- " << dNdeta_phobos_err << std::endl;
0059             std::cout << "dNdeta_ratio = " << dNdeta_ratio << " +/- " << dNdeta_ratio_err << std::endl;
0060             hratio->SetBinContent(j, dNdeta_ratio);
0061             hratio->SetBinError(j, dNdeta_ratio_err);
0062         }
0063         v_sphenix_dNdeta_ratio.push_back(hratio);
0064     }
0065
0066     TCanvas *c = new TCanvas("c", "c", 800, 700);
0067     gPad->SetTopMargin(TopMargin);
0068     gPad->SetLeftMargin(LeftMargin);
0069     gPad->SetRightMargin(RightMargin);
0070     gPad->SetBottomMargin(BottomMargin);
0071     c->cd();
0072     for (size_t i = 0; i < v_sphenix_dNdeta_ratio.size(); i++)
0073     {
0074         if (i == 0)
0075         {
0076             v_sphenix_dNdeta_ratio[i]->GetXaxis()->SetTitle("Pseudorapidity #eta");
0077             v_sphenix_dNdeta_ratio[i]->GetXaxis()->SetRangeUser(-1.15, 1.15);
0078             v_sphenix_dNdeta_ratio[i]->GetYaxis()->SetTitle("The closest-match method / The combinatoric method");
0079             v_sphenix_dNdeta_ratio[i]->GetYaxis()->SetTitleSize(0.035);
0080             v_sphenix_dNdeta_ratio[i]->GetYaxis()->SetTitleOffset(2.1);
0081             v_sphenix_dNdeta_ratio[i]->GetYaxis()->SetRangeUser(0.7, 1.3);
0082             v_sphenix_dNdeta_ratio[i]->GetXaxis()->SetNdivisions(505);
0083             // v_sphenix_dNdeta_ratio[i]->GetYaxis()->SetNdivisions(505);
0084             if (drawall)
0085                 v_sphenix_dNdeta_ratio[i]->Draw("PLC PMC");
0086             else
0087             {
0088                 v_sphenix_dNdeta_ratio[i]->SetLineColor(kBlack);
0089                 v_sphenix_dNdeta_ratio[i]->SetMarkerColor(kBlack);
0090                 v_sphenix_dNdeta_ratio[i]->Draw("PE");
0091             }
0092         }
0093         else
0094         {
0095             if (drawall)
0096                 v_sphenix_dNdeta_ratio[i]->Draw("PLC PMC same");
0097         }
0098     }
0099
0100     c->RedrawAxis();
0101
0102     TLegend *leg = new TLegend(gPad->GetLeftMargin() + 0.05,    //
0103                                (drawall) ? 1 - gPad->GetTopMargin() - 0.27 : 1 - gPad->GetTopMargin() - 0.14, //
0104                                (drawall) ? gPad->GetLeftMargin() + 0.8 : gPad->GetLeftMargin() + 0.3,     //
0105                                1 - gPad->GetTopMargin() - 0.04  //
0106     );
0107     leg->SetHeader("Centrality interval");
0108     if (drawall)
0109         leg->SetNColumns(2);
0110     leg->SetBorderSize(0);
0111     leg->SetFillStyle(0);
0112     leg->SetTextSize(0.03);
0113     if (drawall)
0114     {
0115         for (size_t i = 0; i < v_sphenix_dNdeta_ratio.size(); i++)
0116         {
0117             leg->AddEntry(v_sphenix_dNdeta_ratio[i], Form("%d-%d%%", (int)v_sphenix_centralitybin[i], (int)v_sphenix_centralitybin[i + 1]), "l");
0118         }
0119     }
0120     else
0121         leg->AddEntry(v_sphenix_dNdeta_ratio[0], Form("%d-%d%%", cent_low, cent_high), "pl");
0122
0123     leg->Draw();
0124     c->SaveAs(Form("%s/ratio_dNdEta.pdf", plotdir.c_str()));
0125     c->SaveAs(Form("%s/ratio_dNdEta.png", plotdir.c_str()));
0126 }