analysis_INTT/plot/combine_leastsquare.C

0001 #include "./util_combine.h"
0002
0003 std::vector<float> v_centralitybin = {0, 3, 6, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70};
0004
0005 std::vector<std::string> v_phobos_histname = {"h1D_error_DeltaPhi", "h1D_error_ClusAdc", "h1D_error_ClusPhiSize", "h1D_error_Run_segmentation", "h1D_error_statistic", "h1D_error_GeoOffset"};
0006 std::vector<std::string> v_cms_histname = {"hM_maxreldiff_dRcut", "hM_maxreldiff_clusAdcCut", "hM_maxreldiff_clusPhiSizeCut", "hM_maxreldiff_segment", "hM_statunc_nominal", "hM_maxreldiff_GeoOffset"};
0007 std::vector<std::string> v_unclabel = {"Tracklet reconstruction selection", "Cluster ADC cut", "Cluster #phi-size cut", "Run segments", "Statistical uncertainty in corrections", "Geometric offset"};
0008 std::vector<std::string> v_plotname = {"dRcut", "clusAdcCut", "clusPhiSizeCut", "segment", "statUnc", "GeoOffset"};
0009 std::vector<std::string> v_color = {"#f2777a", "#6699cc", "#9999cc", "#99cc99", "#888888", "#ffcc66"};
0010 std::vector<float> v_alpha = {0.7, 0.7, 0.6, 0.5, 0.8, 0.7};
0011 std::vector<int> v_marker = {20, 21, 33, 34, 47, 45};
0012 std::vector<float> v_corr = {}; // the same order as v_phobos_histname/v_cms_histname
0013 int Nbins_dNdeta = 27;
0014 float xbinlow_dNdeta = -2.7;
0015 float xbinhigh_dNdeta = 2.7;
0016
0017 void setbinattr()
0018 {
0019     TFile *f = new TFile("./systematics/Centrality0to3_Zvtxm10p0to10p0_noasel/finalhists_systematics_Centrality0to3_Zvtxm10p0to10p0_noasel.root", "READ");
0020     TH1 *h = (TH1 *)f->Get("hM_final");
0021     h->SetDirectory(0);
0022     Nbins_dNdeta = h->GetNbinsX();
0023     xbinlow_dNdeta = h->GetXaxis()->GetXmin();
0024     xbinhigh_dNdeta = h->GetXaxis()->GetXmax();
0025     f->Close();
0026 }
0027
0028 TH1 *gethist(std::string fname, std::string histname)
0029 {
0030     TFile *f = new TFile(fname.c_str(), "READ");
0031     TH1 *h;
0032     if (!f->Get(histname.c_str()))
0033     {
0034         std::cout << "Warning: cannot find histogram " << histname << " in file " << fname << std::endl;
0035         std::cout << "Set a hisgoram with all bins = 0" << std::endl;
0036         h = new TH1D(histname.c_str(), histname.c_str(), Nbins_dNdeta, xbinlow_dNdeta, xbinhigh_dNdeta);
0037         for (int i = 1; i <= Nbins_dNdeta; i++)
0038         {
0039             h->SetBinContent(i, 0);
0040         }
0041     }
0042     else
0043     {
0044         h = (TH1 *)f->Get(histname.c_str());
0045     }
0046     h->SetDirectory(0);
0047     f->Close();
0048     return h;
0049 }
0050
0051 void sysunccorr()
0052 {
0053     system("mkdir -p ./systematics/correlation/");
0054
0055     std::vector<TGraph *> v_g_unc;
0056
0057     for (size_t j = 0; j < v_phobos_histname.size(); j++)
0058     {
0059         std::vector<float> v_unc_phobos, v_unc_cms;
0060         v_unc_phobos.clear();
0061         v_unc_cms.clear();
0062
0063         for (size_t i = 0; i < v_centralitybin.size() - 1; i++)
0064         {
0065             std::string centstr = Form("Centrality%dto%d", (int)v_centralitybin[i], (int)v_centralitybin[i + 1]);
0066             std::string fphobos = Form("/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280_HR_Jan172025/Run4/EvtVtxZ/FinalResult/completed/vtxZ_-10_10cm_MBin%d/Final_Mbin%d_00054280/Final_Mbin%d_00054280.root", GetMbinNum(centstr), GetMbinNum(centstr), GetMbinNum(centstr));
0067             std::string fcms = Form("./systematics/Centrality%dto%d_Zvtxm10p0to10p0_noasel/finalhists_systematics_Centrality%dto%d_Zvtxm10p0to10p0_noasel.root", (int)v_centralitybin[i], (int)v_centralitybin[i + 1], (int)v_centralitybin[i], (int)v_centralitybin[i + 1]);
0068
0069             TH1 *h_phobos = gethist(fphobos, v_phobos_histname[j]);
0070             TH1 *h_cms = gethist(fcms, v_cms_histname[j]);
0071
0072             // now loop over the bins with non-zero content
0073             for (int k = 1; k <= h_phobos->GetNbinsX(); k++)
0074             {
0075                 float unc_phobos = h_phobos->GetBinContent(k);
0076                 float unc_cms = h_cms->GetBinContent(k);
0077                 if (unc_cms == 0)
0078                     continue;
0079                 v_unc_phobos.push_back(unc_phobos * 100);
0080                 v_unc_cms.push_back(unc_cms * 100);
0081             }
0082         }
0083
0084         TGraph *g_unc = new TGraph(v_unc_phobos.size(), v_unc_phobos.data(), v_unc_cms.data());
0085         v_corr.push_back(g_unc->GetCorrelationFactor());
0086         v_g_unc.push_back(g_unc);
0087         drawTGraph(g_unc, "PHOBOS approach [%]", "CMS approach [%]", {"Uncertainty: " + v_unclabel[j]}, Form("./systematics/correlation/correlation_%s", v_plotname[j].c_str()));
0088     }
0089
0090     drawTMultiGraph(v_g_unc, "PHOBOS approach [%]", "CMS approach [%]", "./systematics/correlation/correlation_all");
0091 }
0092
0093 std::tuple<float, float> combined_measurement(float X_phobos,                //
0094                                               float X_cms,                   //
0095                                               std::vector<float> unc_phobos, //
0096                                               std::vector<float> unc_cms,    //
0097                                               std::vector<float> corr        //
0098 )
0099 {
0100     int num_sources = unc_phobos.size();
0101     float sigma1_sq = 0., sigma2_sq = 0., covariance = 0.;
0102     for (int i = 0; i < num_sources; i++)
0103     {
0104         std::cout << "  Uncertainty - " << v_unclabel[i] << " (phobos, cms, correlation): " << unc_phobos[i] << ", " << unc_cms[i] << ", " << corr[i] << std::endl;
0105         sigma1_sq += unc_phobos[i] * unc_phobos[i];
0106         sigma2_sq += unc_cms[i] * unc_cms[i];
0107         covariance += unc_phobos[i] * unc_cms[i] * corr[i];
0108     }
0109     // Covariance matrix V = [[V_11, V_12], [V_21, V_22]] = [[sigma1_sq, covariance], [covariance, sigma2_sq]]
0110     float V_11 = sigma1_sq;
0111     float V_12 = covariance; // = V_21
0112     float V_22 = sigma2_sq;
0113     // weights
0114     float denominator = V_11 + V_22 - 2 * covariance;
0115     float w1 = (V_22 - covariance) / denominator;
0116     float w2 = (V_11 - covariance) / denominator;
0117
0118     std::cout << "  V_11 (sigma1_sq) = " << V_11 << ", V_12 (covariance) = " << V_12 << ", V_22 (sigma2_sq) = " << V_22 << std::endl                                                                          //
0119               << "  denominator (V_11 + V_22 - 2 * covariance) = " << denominator << ", w1 ((V_22 - covariance) / denominator) = " << w1 << ", w2 ((V_11 - covariance) / denominator) = " << w2 << std::endl; //
0120     float X_combined = (X_phobos * sigma2_sq + X_cms * sigma1_sq - covariance * (X_phobos + X_cms)) / denominator;
0121     float sigma_combined = std::sqrt((V_11 * V_22 - V_12 * V_12) / (V_11 + V_22 - 2 * V_12));
0122
0123     return std::make_tuple(X_combined, sigma_combined);
0124 }
0125
0126 void combine_leastsquare()
0127 {
0128     setbinattr();
0129
0130     // First, calculate the correlation between the systematic uncertainties in the PHOBOS and CMS approaches
0131     sysunccorr();
0132
0133     system("mkdir -p ./systematics/combined/");
0134
0135     std::cout << "Correlation in the uncertainty between PHOBOS and CMS approaches:" << std::endl;
0136     for (size_t i = 0; i < v_corr.size(); i++)
0137     {
0138         std::cout << v_unclabel[i] << ": " << v_corr[i] << std::endl;
0139     }
0140
0141     // start the loop over the centrality bins
0142     for (size_t i = 0; i < v_centralitybin.size() - 1; i++)
0143     // for (size_t i = 0; i < 1; i++)
0144     {
0145         std::cout << "Combine the PHOBOS and CMS measurements for centrality bin " << v_centralitybin[i] << " to " << v_centralitybin[i + 1] << std::endl;
0146
0147         std::string centstr = Form("Centrality%dto%d", (int)v_centralitybin[i], (int)v_centralitybin[i + 1]);
0148         std::string fphobos = Form("/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280_HR_Jan172025/Run4/EvtVtxZ/FinalResult/completed/vtxZ_-10_10cm_MBin%d/Final_Mbin%d_00054280/Final_Mbin%d_00054280.root", GetMbinNum(centstr), GetMbinNum(centstr), GetMbinNum(centstr));
0149         std::string fcms = Form("./systematics/Centrality%dto%d_Zvtxm10p0to10p0_noasel/finalhists_systematics_Centrality%dto%d_Zvtxm10p0to10p0_noasel.root", (int)v_centralitybin[i], (int)v_centralitybin[i + 1], (int)v_centralitybin[i], (int)v_centralitybin[i + 1]);
0150
0151         // histograms for the central values
0152         TH1 *h_phobos_X = gethist(fphobos, "h1D_dNdEta_reco");
0153         TH1 *h_cms_X = gethist(fcms, "hM_final");
0154         // histograms for the total uncertainties, for checking
0155         TH1 *h_phobos_totlunc = gethist(fphobos, "h1D_error_Final");
0156         TH1 *h_cms_totlunc = gethist(fcms, "hM_TotalRelUnc");
0157         // histograms for the uncertainties
0158         std::vector<TH1 *> v_h_phobos_unc, v_h_cms_unc;
0159         v_h_phobos_unc.clear();
0160         v_h_cms_unc.clear();
0161         for (size_t j = 0; j < v_phobos_histname.size(); j++)
0162         {
0163             v_h_phobos_unc.push_back(gethist(fphobos, v_phobos_histname[j]));
0164             v_h_cms_unc.push_back(gethist(fcms, v_cms_histname[j]));
0165         }
0166
0167         // histogram for the combined measurement
0168         TH1 *h_combined = new TH1D(Form("h_combined_%s", centstr.c_str()), Form("Combined dN/d#eta for %s", centstr.c_str()), h_phobos_X->GetNbinsX(), h_phobos_X->GetXaxis()->GetXmin(), h_phobos_X->GetXaxis()->GetXmax());
0169
0170         // loop over the bins
0171         for (int k = 1; k <= h_phobos_X->GetNbinsX(); k++)
0172         {
0173             float X_phobos = h_phobos_X->GetBinContent(k);
0174             float X_cms = h_cms_X->GetBinContent(k);
0175             for (int b = 0; b < h_cms_X->GetNbinsX(); b++)
0176             {
0177                 if (h_cms_X->GetBinContent(b) < 0.5)
0178                 {
0179                     h_cms_X->SetBinContent(b, 0);
0180                     h_cms_X->SetBinError(b, 0);
0181                 }
0182             }
0183
0184             std::vector<float> unc_phobos, unc_cms;
0185             unc_phobos.clear();
0186             unc_cms.clear();
0187             for (size_t j = 0; j < v_phobos_histname.size(); j++)
0188             {
0189                 unc_phobos.push_back(v_h_phobos_unc[j]->GetBinContent(k) * X_phobos);
0190                 unc_cms.push_back(v_h_cms_unc[j]->GetBinContent(k) * X_cms);
0191             }
0192
0193             std::cout << "----- Bin " << k << ", eta = " << h_phobos_X->GetBinCenter(k) << " -----" << std::endl;
0194             std::cout << "PHOBOS measurement: " << X_phobos << " +/- " << h_phobos_totlunc->GetBinContent(k) * X_phobos << std::endl;
0195             std::cout << "CMS measurement: " << X_cms << " +/- " << h_cms_totlunc->GetBinContent(k) * X_cms << std::endl;
0196             std::tuple<float, float> combined = combined_measurement(X_phobos, X_cms, unc_phobos, unc_cms, v_corr);
0197             h_combined->SetBinContent(k, std::get<0>(combined));
0198             h_combined->SetBinError(k, std::get<1>(combined));
0199
0200             // print out the details
0201             std::cout << "Combined measurement: " << std::get<0>(combined) << " +/- " << std::get<1>(combined) << std::endl;
0202             std::cout << "---------------------------------" << std::endl;
0203         }
0204
0205         // draw the combined measurement
0206         TCanvas *c = new TCanvas("c", "c", 800, 700);
0207         c->cd();
0208         gPad->SetRightMargin(0.07);
0209         gPad->SetTopMargin(0.07);
0210         gPad->SetLeftMargin(0.17);
0211         h_phobos_X->SetLineColor(kRed);
0212         h_phobos_X->SetLineWidth(0);
0213         h_phobos_X->SetMarkerColor(kRed);
0214         h_phobos_X->SetMarkerStyle(20);
0215         h_phobos_X->SetMarkerSize(1);
0216         h_phobos_X->SetFillColorAlpha(kRed, 0.1);
0217         h_phobos_X->GetXaxis()->SetTitle("#eta");
0218         h_phobos_X->GetYaxis()->SetTitle("dN_{ch}/d#eta");
0219         h_phobos_X->GetYaxis()->SetTitleOffset(1.7);
0220         h_phobos_X->GetYaxis()->SetRangeUser(h_cms_X->GetMinimum(0) * 0.7, h_cms_X->GetMaximum() * 1.3);
0221         h_phobos_X->Draw("E2P");
0222         h_cms_X->SetLineColor(kBlue);
0223         h_cms_X->SetLineWidth(0);
0224         h_cms_X->SetMarkerColor(kBlue);
0225         h_cms_X->SetMarkerStyle(21);
0226         h_cms_X->SetMarkerSize(1);
0227         h_cms_X->SetFillColorAlpha(kBlue, 0.1);
0228         h_cms_X->Draw("E2P SAME");
0229         h_combined->SetLineColor(kBlack);
0230         h_combined->SetLineWidth(1);
0231         h_combined->SetMarkerColor(kBlack);
0232         h_combined->SetMarkerStyle(33);
0233         h_combined->SetMarkerSize(1);
0234         h_combined->SetFillColorAlpha(kBlack, 0.4);
0235         h_combined->Draw("E2P SAME");
0236         TLegend *leg = new TLegend(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.05, gPad->GetLeftMargin() + 0.25, 1 - gPad->GetTopMargin() - 0.25);
0237         leg->SetHeader(Form("Centrality %d-%d%%", (int)v_centralitybin[i], (int)v_centralitybin[i + 1]));
0238         leg->SetTextAlign(kHAlignLeft + kVAlignCenter);
0239         leg->SetTextSize(0.03);
0240         leg->AddEntry(h_phobos_X, "PHOBOS approach", "pf");
0241         leg->AddEntry(h_cms_X, "CMS approach", "pf");
0242         leg->AddEntry(h_combined, "Combined measurement", "pf");
0243         leg->SetBorderSize(0);
0244         leg->SetFillStyle(0);
0245         leg->Draw();
0246         gPad->Modified();
0247         gPad->Update();
0248         c->SaveAs(Form("./systematics/combined/combined_%s_covmatrix.pdf", centstr.c_str()));
0249         c->SaveAs(Form("./systematics/combined/combined_%s_covmatrix.png", centstr.c_str()));
0250     }
0251 }