analysis_INTT/plot/combine_CMSHIN10001.C

0001 #include "./util_combine.h"
0002
0003 std::vector<double> 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_Strangeness", "h1D_error_Generator", "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_strangeness", "hM_maxreldiff_eventgen", "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", "Secondaries", "Event generator", "Geometric offset"};
0008 std::vector<std::string> v_plotname = {"dRcut", "clusAdcCut", "clusPhiSizeCut", "segment", "statUnc", "Secondaries", "EventGenerator", "GeoOffset"};
0009 std::vector<std::string> v_color = {"#f2777a", "#6699cc", "#9999cc", "#99cc99", "#e99960", "#FFC0CB", "#ffcc66", "#7FE9DE"};
0010 std::vector<double> v_alpha = {0.7, 0.7, 0.6, 0.5, 0.8, 0.7, 0.7, 0.7};
0011 std::vector<int> v_marker = {20, 21, 33, 34, 47, 45, 43, 41};
0012 std::vector<double> v_corr = {}; // the same order as v_phobos_histname/v_cms_histname
0013 int Nbins_dNdeta = 27;
0014 double xbinlow_dNdeta = -2.7;
0015 double xbinhigh_dNdeta = 2.7;
0016
0017 std::vector<float> v_weightratio = {};
0018
0019 void setbinattr()
0020 {
0021     TFile *f = new TFile("./systematics/Centrality0to3_Zvtxm10p0to10p0_noasel/finalhists_systematics_Centrality0to3_Zvtxm10p0to10p0_noasel.root", "READ");
0022     TH1 *h = (TH1 *)f->Get("hM_final");
0023     h->SetDirectory(0);
0024     Nbins_dNdeta = h->GetNbinsX();
0025     xbinlow_dNdeta = h->GetXaxis()->GetXmin();
0026     xbinhigh_dNdeta = h->GetXaxis()->GetXmax();
0027     f->Close();
0028 }
0029
0030 TH1 *gethist(std::string fname, std::string histname)
0031 {
0032     TFile *f = new TFile(fname.c_str(), "READ");
0033     TH1 *h;
0034     if (!f->Get(histname.c_str()))
0035     {
0036         std::cout << "Warning: cannot find histogram " << histname << " in file " << fname << std::endl;
0037         std::cout << "Set a hisgoram with all bins = 0" << std::endl;
0038         h = new TH1D(histname.c_str(), histname.c_str(), Nbins_dNdeta, xbinlow_dNdeta, xbinhigh_dNdeta);
0039         for (int i = 1; i <= Nbins_dNdeta; i++)
0040         {
0041             h->SetBinContent(i, 0);
0042         }
0043     }
0044     else
0045     {
0046         h = (TH1 *)f->Get(histname.c_str());
0047     }
0048     h->SetDirectory(0);
0049     f->Close();
0050     return h;
0051 }
0052
0053 void drawTGraph(TGraph *g, std::string xtitle, std::string ytitle, std::vector<std::string> plotinfo, std::string plotname)
0054 {
0055     TCanvas *c = new TCanvas("c", "c", 800, 700);
0056     gPad->SetRightMargin(0.07);
0057     gPad->SetTopMargin(0.07);
0058     gPad->SetLeftMargin(0.17);
0059     c->cd();
0060     g->SetMarkerStyle(20);
0061     g->SetMarkerSize(1);
0062     g->GetXaxis()->SetTitle(xtitle.c_str());
0063     g->GetYaxis()->SetTitle(ytitle.c_str());
0064     g->GetYaxis()->SetTitleOffset(1.7);
0065     // g->GetXaxis()->SetLimits(g->GetXaxis()->GetXmin() * 0.9, g->GetXaxis()->GetXmax() * 1.1);
0066     g->GetXaxis()->SetLimits(0, g->GetXaxis()->GetXmax() * 1.1);
0067     g->GetHistogram()->SetMaximum(g->GetYaxis()->GetXmax() * 1.1);
0068     g->GetHistogram()->SetMinimum(0);
0069     g->Draw("AP");
0070     TLatex *t = new TLatex();
0071     t->SetTextSize(0.035);
0072     for (size_t i = 0; i < plotinfo.size(); ++i)
0073     {
0074         t->DrawLatexNDC(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.08 - i * 0.045, plotinfo[i].c_str());
0075     }
0076     // covariance and correlation factors
0077     t->DrawLatexNDC(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.08 - plotinfo.size() * 0.045, Form("Covariance: %.3g", g->GetCovariance()));
0078     t->DrawLatexNDC(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.08 - (plotinfo.size() + 1) * 0.045, Form("Correlation: %.3g", g->GetCorrelationFactor()));
0079     c->SaveAs(Form("%s.pdf", plotname.c_str()));
0080     c->SaveAs(Form("%s.png", plotname.c_str()));
0081 }
0082
0083 void drawTMultiGraph(std::vector<TGraph *> v_graph, std::string xtitle, std::string ytitle, std::string plotname)
0084 {
0085     TCanvas *c = new TCanvas("c", "c", 800, 700);
0086     c->cd();
0087     TMultiGraph *mg = new TMultiGraph();
0088     for (size_t i = 0; i < v_graph.size(); i++)
0089     {
0090         v_graph[i]->SetMarkerColorAlpha(TColor::GetColor(v_color[i].c_str()), v_alpha[i]);
0091         v_graph[i]->SetMarkerStyle(v_marker[i]);
0092         v_graph[i]->SetMarkerSize(0.8);
0093         mg->Add(v_graph[i]);
0094     }
0095     mg->Draw("AP");
0096     mg->GetXaxis()->SetTitle(xtitle.c_str());
0097     mg->GetYaxis()->SetTitle(ytitle.c_str());
0098     mg->GetYaxis()->SetTitleOffset(1.5);
0099     gPad->Modified();
0100     float axismax = (mg->GetXaxis()->GetXmax() > mg->GetYaxis()->GetXmax()) ? mg->GetXaxis()->GetXmax() : mg->GetYaxis()->GetXmax();
0101     mg->GetXaxis()->SetLimits(0, axismax * 1.1);
0102     mg->GetHistogram()->SetMaximum(axismax * 1.1);
0103     mg->GetHistogram()->SetMinimum(0);
0104     TLegend *leg = new TLegend(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.05, gPad->GetLeftMargin() + 0.25, 1 - gPad->GetTopMargin() - 0.05 - 0.05 * v_graph.size());
0105     leg->SetTextAlign(kHAlignLeft + kVAlignTop);
0106     leg->SetTextSize(0.03);
0107     for (size_t i = 0; i < v_graph.size(); i++)
0108     {
0109         leg->AddEntry(v_graph[i], Form("%s (corr. coefficient = %.3g)", v_unclabel[i].c_str(), v_graph[i]->GetCorrelationFactor()), "p");
0110     }
0111     leg->SetBorderSize(0);
0112     leg->SetFillStyle(0);
0113     leg->Draw();
0114     gPad->Modified();
0115     gPad->Update();
0116     c->SaveAs(Form("%s.pdf", plotname.c_str()));
0117     c->SaveAs(Form("%s.png", plotname.c_str()));
0118 }
0119
0120 // calculate the correlation in the systematic uncertainties between the PHOBOS and CMS approaches
0121 void sysunccorr()
0122 {
0123     system("mkdir -p ./systematics/correlation/");
0124
0125     std::vector<TGraph *> v_g_unc;
0126
0127     for (size_t j = 0; j < v_phobos_histname.size(); j++)
0128     {
0129         std::vector<double> v_unc_phobos, v_unc_cms;
0130         v_unc_phobos.clear();
0131         v_unc_cms.clear();
0132
0133         for (size_t i = 0; i < v_centralitybin.size() - 1; i++)
0134         {
0135             std::string centstr = Form("Centrality%dto%d", (int)v_centralitybin[i], (int)v_centralitybin[i + 1]);
0136             std::string fphobos = 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", GetMbinNum(centstr), GetMbinNum(centstr), GetMbinNum(centstr));
0137             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]);
0138
0139             TH1 *h_phobos = gethist(fphobos, v_phobos_histname[j]);
0140             TH1 *h_cms = gethist(fcms, v_cms_histname[j]);
0141
0142             // now loop over the bins with non-zero content
0143             for (int k = 1; k <= h_phobos->GetNbinsX(); k++)
0144             {
0145                 // if the bin center is outside |eta|<= 1.1, skip
0146                 double bincenter = h_phobos->GetXaxis()->GetBinCenter(k);
0147                 if (std::abs(bincenter) > 1.1)
0148                     continue;
0149
0150                 double unc_phobos = h_phobos->GetBinContent(k);
0151                 double unc_cms = h_cms->GetBinContent(k);
0152                 if (unc_cms == 0)
0153                     continue;
0154                 v_unc_phobos.push_back(unc_phobos * 100);
0155                 v_unc_cms.push_back(unc_cms * 100);
0156             }
0157         }
0158
0159         TGraph *g_unc = new TGraph(v_unc_phobos.size(), v_unc_phobos.data(), v_unc_cms.data());
0160         v_corr.push_back(g_unc->GetCorrelationFactor());
0161         v_g_unc.push_back(g_unc);
0162         drawTGraph(g_unc, "The combinatoric method [%]", "The closest-match method [%]", {"Uncertainty: " + v_unclabel[j]}, Form("./systematics/correlation/correlation_%s", v_plotname[j].c_str()));
0163     }
0164
0165     drawTMultiGraph(v_g_unc, "The combinatoric method [%]", "The closest-match method [%]", "./systematics/correlation/correlation_all");
0166 }
0167
0168 // combination method according to the method in CMS HIN-10-001
0169 // return <the combined value, uncorrelated uncertainty, correlated uncertainty, total uncertainty>
0170 std::tuple<double, double, double, double> combined_measurement(double X_phobos,                //
0171                                                                 double X_cms,                   //
0172                                                                 std::vector<double> unc_phobos, //
0173                                                                 std::vector<double> unc_cms,    //
0174                                                                 std::vector<double> corr        //
0175 )
0176 {
0177     int num_sources = unc_phobos.size();
0178     double unc_corr = 0., unc_uncorr = 0., unc_uncorr_cms = 0., unc_uncorr_phobos = 0., unc_total = 0.;
0179     for (int i = 0; i < num_sources; i++)
0180     {
0181         // convert the number to only include 3 digit after the decimal point
0182         corr[i] = std::round(corr[i] * 100.) / 100.;
0183         if (fabs(corr[i]) <= 0.10 && (unc_phobos[i] != 0 && unc_cms[i] != 0)) // treated as uncorrelated uncertainty -> quadrature sum
0184         {
0185             std::cout << "  Uncertainty (uncorrelated) - " << v_unclabel[i] << " (phobos, cms, correlation): " << unc_phobos[i] << ", " << unc_cms[i] << ", " << corr[i] << std::endl;
0186             // unc_uncorr += unc_phobos[i] * unc_phobos[i] + unc_cms[i] * unc_cms[i];
0187             unc_uncorr_phobos += unc_phobos[i] * unc_phobos[i];
0188             unc_uncorr_cms += unc_cms[i] * unc_cms[i];
0189         }
0190         // else // treated as fully correlated uncertainty -> sum quadrature
0191         // {
0192         //     std::cout << "  Uncertainty (correlated) - " << v_unclabel[i] << " (phobos, cms, correlation): " << unc_phobos[i] << ", " << unc_cms[i] << ", " << corr[i] << std::endl;
0193         //     unc_corr += std::pow(unc_phobos[i] + unc_cms[i], 2);
0194         // }
0195     }
0196     unc_uncorr_phobos = std::sqrt(unc_uncorr_phobos);
0197     unc_uncorr_cms = std::sqrt(unc_uncorr_cms);
0198     cout << "----- unc_uncorr_phobos = " << unc_uncorr_phobos << ", unc_uncorr_cms = " << unc_uncorr_cms << endl;
0199
0200     // the combined result is the weighted average of the two measurements based on the uncorrelated unccertainty, with the weights being the inverse square of the total uncertainty w_i = 1/(sigma_i^2)
0201     double w_phobos = 1. / (unc_uncorr_phobos * unc_uncorr_phobos);
0202     double w_cms = 1. / (unc_uncorr_cms * unc_uncorr_cms);
0203     double w_total = w_phobos + w_cms;
0204     double X_combined = (w_phobos * X_phobos + w_cms * X_cms) / w_total;
0205     // now combine the uncorrelated uncertainties: sum the weights and take the inverse square root
0206     unc_uncorr = 1. / std::sqrt(w_total);
0207
0208     cout << "----- w_phobos = " << w_phobos << ", w_cms = " << w_cms << ", w_total = " << w_total << endl;
0209     if (unc_uncorr_phobos != 0 && unc_uncorr_cms != 0 && X_phobos > 0)
0210     {
0211         v_weightratio.push_back(w_cms / w_total);
0212         std::cout << "Weight ratio (w_cms/w_total) = " << w_cms / w_total << std::endl;
0213     }
0214
0215     // now combine the correlated uncertainties: sum the squares of the correlated uncertainties and take the square root
0216     for (int i = 0; i < num_sources; i++)
0217     {
0218         // convert the number to only include 3 digit after the decimal point
0219         corr[i] = std::round(corr[i] * 100.) / 100.;
0220         if (fabs(corr[i]) > 0.10 && (unc_phobos[i] != 0 && unc_cms[i] != 0)) // treated as fully correlated uncertainty -> sum quadrature
0221         {
0222             std::cout << "  Uncertainty (correlated) - " << v_unclabel[i] << " (phobos, cms, correlation): " << unc_phobos[i] << ", " << unc_cms[i] << ", " << corr[i] << std::endl;
0223             unc_corr += std::pow((w_phobos / w_total) * unc_phobos[i] + (w_cms / w_total) * unc_cms[i], 2);
0224         }
0225     }
0226
0227     unc_corr = std::sqrt(unc_corr);
0228     unc_total = std::sqrt(unc_corr * unc_corr + unc_uncorr * unc_uncorr);
0229
0230     std::cout << "  Uncorrelated uncertainty (phobos, cms): " << unc_uncorr_phobos << ", " << unc_uncorr_cms << " -> (w_phobos, w_cms)=(" << w_phobos << ", " << w_cms << ") -> w_total=" << w_total << std::endl;
0231     std::cout << "  Total (uncorrelated, correlated) uncertainty = (" << unc_uncorr << ", " << unc_corr << ") -> Total uncertainty = " << unc_total << std::endl;
0232     std::cout << "  Measurement: " << X_combined << " +/- " << unc_total << std::endl;
0233
0234     return std::make_tuple(X_combined, unc_uncorr, unc_corr, unc_total);
0235 }
0236
0237 void combine_CMSHIN10001()
0238 {
0239     setbinattr();
0240
0241     // First, calculate the correlation between the systematic uncertainties in the PHOBOS and CMS approaches
0242     sysunccorr();
0243
0244     system("mkdir -p ./systematics/combined/");
0245
0246     std::cout << "Correlation in the uncertainty between PHOBOS and CMS approaches:" << std::endl;
0247     for (size_t i = 0; i < v_corr.size(); i++)
0248     {
0249         std::cout << v_unclabel[i] << ": " << v_corr[i] << std::endl;
0250     }
0251
0252     std::vector<TGraphAsymmErrors *> v_tgae_combined;
0253     v_tgae_combined.clear();
0254     std::vector<TGraphAsymmErrors *> v_tgae_uncorr;
0255     v_tgae_uncorr.clear();
0256     std::vector<TGraphAsymmErrors *> v_tgae_corr;
0257     v_tgae_corr.clear();
0258     std::vector<TH1 *> v_combined;
0259     v_combined.clear();
0260
0261     float uncorr_unc_min = 1E6, uncorr_unc_max = -1E6, corr_unc_min = 1E6, corr_unc_max = -1E6, total_unc_min = 1E6, total_unc_max = -1E6; // in percent
0262
0263     // start the loop over the centrality bins
0264     for (size_t i = 0; i < v_centralitybin.size() - 1; i++)
0265     // for (size_t i = 0; i < 1; i++)
0266     {
0267         std::cout << "Combine the PHOBOS and CMS measurements for centrality bin " << v_centralitybin[i] << " to " << v_centralitybin[i + 1] << std::endl;
0268
0269         std::string centstr = Form("Centrality%dto%d", (int)v_centralitybin[i], (int)v_centralitybin[i + 1]);
0270         std::string fphobos = 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", GetMbinNum(centstr), GetMbinNum(centstr), GetMbinNum(centstr));
0271         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]);
0272
0273         // histograms for the central values
0274         TH1 *h_phobos_X = gethist(fphobos, "h1D_dNdEta_reco");
0275         TH1 *h_cms_X = gethist(fcms, "hM_final");
0276         // histograms for the total uncertainties, for checking
0277         TH1 *h_phobos_totlunc = gethist(fphobos, "h1D_error_Final");
0278         TH1 *h_cms_totlunc = gethist(fcms, "hM_TotalRelUnc");
0279         // histograms for the uncertainties
0280         std::vector<TH1 *> v_h_phobos_unc, v_h_cms_unc;
0281         v_h_phobos_unc.clear();
0282         v_h_cms_unc.clear();
0283         for (size_t j = 0; j < v_phobos_histname.size(); j++)
0284         {
0285             TH1F *hp = (TH1F *)gethist(fphobos, v_phobos_histname[j]);
0286             TH1F *hc = (TH1F *)gethist(fcms, v_cms_histname[j]);
0287             v_h_phobos_unc.push_back(hp);
0288             v_h_cms_unc.push_back(hc);
0289         }
0290
0291         // histogram for the combined measurement
0292         // 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());
0293
0294         std::vector<double> v_combined_x, v_combined_y, v_combined_xerr, v_combined_yerrup, v_combined_yerrdown, v_combined_uncorrup, v_combined_uncorrdown, v_combined_corrup, v_combined_corrdown;
0295         // loop over the bins
0296         for (int k = 1; k <= h_phobos_X->GetNbinsX(); k++)
0297         {
0298             double X_phobos = h_phobos_X->GetBinContent(k);
0299             double X_cms = h_cms_X->GetBinContent(k);
0300             for (int b = 0; b < h_cms_X->GetNbinsX(); b++)
0301             {
0302                 if (h_cms_X->GetBinContent(b) < 0.5)
0303                 {
0304                     h_cms_X->SetBinContent(b, 0);
0305                     h_cms_X->SetBinError(b, 0);
0306                 }
0307             }
0308
0309             std::vector<double> unc_phobos, unc_cms;
0310             unc_phobos.clear();
0311             unc_cms.clear();
0312             for (size_t j = 0; j < v_phobos_histname.size(); j++)
0313             {
0314                 unc_phobos.push_back(v_h_phobos_unc[j]->GetBinContent(k) * X_phobos);
0315                 unc_cms.push_back(v_h_cms_unc[j]->GetBinContent(k) * X_cms);
0316             }
0317
0318             std::cout << "----- Bin " << k << ", eta = " << h_phobos_X->GetBinCenter(k) << " -----" << std::endl;
0319             std::cout << "PHOBOS measurement: " << X_phobos << " +/- " << h_phobos_totlunc->GetBinContent(k) * X_phobos << " (" << h_phobos_totlunc->GetBinContent(k) * 100. << "%)" << std::endl;
0320             std::cout << "CMS measurement: " << X_cms << " +/- " << h_cms_totlunc->GetBinContent(k) * X_cms << " (" << h_cms_totlunc->GetBinContent(k) * 100. << "%)" << std::endl;
0321             std::tuple<double, double, double, double> combined = combined_measurement(X_phobos, X_cms, unc_phobos, unc_cms, v_corr); // return <the combined value, uncorrelated uncertainty, correlated uncertainty, total uncertainty>
0322
0323             // print out the details
0324             std::cout << "Combined measurement: " << std::get<0>(combined) << " +/- " << std::get<3>(combined) << " (" << std::get<3>(combined) / std::get<0>(combined) * 100. << "%)" << std::endl;
0325             std::cout << "---------------------------------" << std::endl;
0326
0327             // only when both measurements are non-zero or NAN (not a number) store the combined measurement
0328             if (X_phobos == 0 || X_cms == 0 || std::isnan(X_phobos) || std::isnan(X_cms) || std::isnan(std::get<0>(combined)) || std::isnan(std::get<3>(combined)))
0329                 continue;
0330
0331             // also don"t store the measurement with |eta| >= 1.1
0332             if (std::abs(h_phobos_X->GetBinCenter(k)) >= 1.1)
0333                 continue;
0334
0335             // assign the combined measurement to the histogram
0336             v_combined_x.push_back(h_phobos_X->GetBinCenter(k));
0337             v_combined_y.push_back(std::get<0>(combined));
0338             v_combined_xerr.push_back(h_phobos_X->GetBinWidth(k) * 0.49);
0339             v_combined_yerrup.push_back(std::get<3>(combined));
0340             v_combined_yerrdown.push_back(std::get<3>(combined));
0341             v_combined_uncorrup.push_back(std::get<1>(combined));
0342             v_combined_uncorrdown.push_back(std::get<1>(combined));
0343             v_combined_corrup.push_back(std::get<2>(combined));
0344             v_combined_corrdown.push_back(std::get<2>(combined));
0345
0346             if ((std::get<1>(combined) / std::get<0>(combined) * 100.) < uncorr_unc_min)
0347                 uncorr_unc_min = std::get<1>(combined) / std::get<0>(combined) * 100.;
0348             if ((std::get<1>(combined) / std::get<0>(combined) * 100.) > uncorr_unc_max)
0349                 uncorr_unc_max = std::get<1>(combined) / std::get<0>(combined) * 100.;
0350             if ((std::get<2>(combined) / std::get<0>(combined) * 100.) < corr_unc_min)
0351                 corr_unc_min = std::get<2>(combined) / std::get<0>(combined) * 100.;
0352             if ((std::get<2>(combined) / std::get<0>(combined) * 100.) > corr_unc_max)
0353                 corr_unc_max = std::get<2>(combined) / std::get<0>(combined) * 100.;
0354             if ((std::get<3>(combined) / std::get<0>(combined) * 100.) < total_unc_min)
0355                 total_unc_min = std::get<3>(combined) / std::get<0>(combined) * 100.;
0356             if ((std::get<3>(combined) / std::get<0>(combined) * 100.) > total_unc_max)
0357                 total_unc_max = std::get<3>(combined) / std::get<0>(combined) * 100.;
0358         }
0359
0360         // set up the TGraphAssymErrors for the combined measurement
0361         TGraphAsymmErrors *g_combined = new TGraphAsymmErrors(v_combined_x.size(), v_combined_x.data(), v_combined_y.data(), v_combined_xerr.data(), v_combined_xerr.data(), v_combined_yerrdown.data(), v_combined_yerrup.data());
0362         g_combined->SetName(Form("tgae_combine_%s", centstr.c_str()));
0363         // make histogram from the TGraphAsymmErrors
0364         TH1 *h_combined = new TH1D(Form("h_combined_%s", centstr.c_str()), Form("Combined dN/d#eta for %s", centstr.c_str()), g_combined->GetN(), g_combined->GetX()[0] - g_combined->GetErrorXlow(0), g_combined->GetX()[g_combined->GetN() - 1] + g_combined->GetErrorXhigh(g_combined->GetN() - 1));
0365         for (int k = 0; k < g_combined->GetN(); k++)
0366         {
0367             h_combined->SetBinContent(h_combined->FindBin(g_combined->GetX()[k]), g_combined->GetY()[k]);
0368             h_combined->SetBinError(h_combined->FindBin(g_combined->GetX()[k]), g_combined->GetErrorY(k));
0369         }
0370         v_tgae_combined.push_back(g_combined);
0371         v_combined.push_back(h_combined);
0372         // set up the TGraphAssymErrors for the uncorrelated uncertainty
0373         TGraphAsymmErrors *g_uncorr = new TGraphAsymmErrors(v_combined_x.size(), v_combined_x.data(), v_combined_y.data(), v_combined_xerr.data(), v_combined_xerr.data(), v_combined_uncorrdown.data(), v_combined_uncorrup.data());
0374         g_uncorr->SetName(Form("tgae_uncorr_%s", centstr.c_str()));
0375         v_tgae_uncorr.push_back(g_uncorr);
0376         // set up the TGraphAssymErrors for the correlated uncertainty
0377         TGraphAsymmErrors *g_corr = new TGraphAsymmErrors(v_combined_x.size(), v_combined_x.data(), v_combined_y.data(), v_combined_xerr.data(), v_combined_xerr.data(), v_combined_corrdown.data(), v_combined_corrup.data());
0378         g_corr->SetName(Form("tgae_corr_%s", centstr.c_str()));
0379         v_tgae_corr.push_back(g_corr);
0380
0381         // draw the combined measurement
0382         TCanvas *c = new TCanvas("c", "c", 800, 700);
0383         c->cd();
0384         gPad->SetRightMargin(0.07);
0385         gPad->SetTopMargin(0.07);
0386         gPad->SetLeftMargin(0.17);
0387         h_phobos_X->SetLineColor(kRed);
0388         h_phobos_X->SetLineWidth(0);
0389         h_phobos_X->SetMarkerColor(kRed);
0390         h_phobos_X->SetMarkerStyle(20);
0391         h_phobos_X->SetMarkerSize(1);
0392         h_phobos_X->SetFillColorAlpha(kRed, 0.1);
0393         h_phobos_X->GetXaxis()->SetTitle("#eta");
0394         h_phobos_X->GetYaxis()->SetTitle("dN_{ch}/d#eta");
0395         h_phobos_X->GetYaxis()->SetTitleOffset(1.7);
0396         h_phobos_X->GetYaxis()->SetRangeUser(h_cms_X->GetMinimum(0) * 0.7, h_cms_X->GetMaximum() * 1.3);
0397         h_phobos_X->Draw("E2P");
0398         h_cms_X->SetLineColor(kBlue);
0399         h_cms_X->SetLineWidth(0);
0400         h_cms_X->SetMarkerColor(kBlue);
0401         h_cms_X->SetMarkerStyle(21);
0402         h_cms_X->SetMarkerSize(1);
0403         h_cms_X->SetFillColorAlpha(kBlue, 0.1);
0404         h_cms_X->Draw("E2P SAME");
0405         g_combined->SetLineColor(kBlack);
0406         g_combined->SetLineWidth(1);
0407         g_combined->SetMarkerColor(kBlack);
0408         g_combined->SetMarkerStyle(33);
0409         g_combined->SetMarkerSize(1);
0410         g_combined->SetFillColorAlpha(kBlack, 0.4);
0411         g_combined->Draw("5E P SAME");
0412         TLegend *leg = new TLegend(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.05, gPad->GetLeftMargin() + 0.25, 1 - gPad->GetTopMargin() - 0.25);
0413         leg->SetHeader(Form("Centrality %d-%d%%", (int)v_centralitybin[i], (int)v_centralitybin[i + 1]));
0414         leg->SetTextAlign(kHAlignLeft + kVAlignCenter);
0415         leg->SetTextSize(0.03);
0416         leg->AddEntry(h_phobos_X, "The combinatoric method", "pf");
0417         leg->AddEntry(h_cms_X, "The closest-match method", "pf");
0418         leg->AddEntry(g_combined, "Combined (adapted from JHEP 08 (2011) 141) ", "pf");
0419         leg->SetBorderSize(0);
0420         leg->SetFillStyle(0);
0421         leg->Draw();
0422         gPad->Modified();
0423         gPad->Update();
0424         c->SaveAs(Form("./systematics/combined/combined_%s_CMSHIN10001.pdf", centstr.c_str()));
0425         c->SaveAs(Form("./systematics/combined/combined_%s_CMSHIN10001.png", centstr.c_str()));
0426     }
0427
0428     // TFile to store the combined results
0429     TFile *fout = new TFile("./systematics/combined/combined.root", "RECREATE");
0430     for (size_t i = 0; i < v_tgae_combined.size(); i++)
0431     {
0432         // print all points
0433         // for (int j = 0; j < v_tgae_combined[i]->GetN(); j++)
0434         // {
0435         //     std::cout << "Centrality: " << v_centralitybin[i] << " to " << v_centralitybin[i + 1] << ", eta = " << v_tgae_combined[i]->GetX()[j] << ", dN/deta = " << v_tgae_combined[i]->GetY()[j] << " +/- " << v_tgae_combined[i]->GetErrorY(j) << std::endl;
0436         // }
0437
0438         v_tgae_combined[i]->Write();
0439         v_combined[i]->Write();
0440         v_tgae_uncorr[i]->Write();
0441         v_tgae_corr[i]->Write();
0442     }
0443     fout->Close();
0444
0445     std::cout << "Uncorrelated uncertainty range: " << uncorr_unc_min << "%-" << uncorr_unc_max << "%" << std::endl;
0446     std::cout << "Correlated uncertainty range: " << corr_unc_min << "%-" << corr_unc_max << "%" << std::endl;
0447     std::cout << "Total uncertainty range: " << total_unc_min << "%-" << total_unc_max << "%" << std::endl;
0448
0449     // print out the weight ratio, minimum and maximum
0450     std::cout << "Weight ratio (w_cms/w_total) range: " << *std::min_element(v_weightratio.begin(), v_weightratio.end()) << "-" << *std::max_element(v_weightratio.begin(), v_weightratio.end()) << std::endl;
0451     std::cout << "Weight ratio (w_cms/w_total) average: " << std::accumulate(v_weightratio.begin(), v_weightratio.end(), 0.0) / v_weightratio.size() << std::endl;
0452     std::cout << "Weight ratio (w_cms/w_total) median: " << v_weightratio[v_weightratio.size() / 2] << std::endl;
0453     std::cout << "Weight ratio (w_cms/w_total) standard deviation: " << std::sqrt(std::accumulate(v_weightratio.begin(), v_weightratio.end(), 0.0, [](double sum, double val) { return sum + (val - (std::accumulate(v_weightratio.begin(), v_weightratio.end(), 0.0) / v_weightratio.size())) * (val - (std::accumulate(v_weightratio.begin(), v_weightratio.end(), 0.0) / v_weightratio.size())); }) / v_weightratio.size()) << std::endl;
0454     std::cout << "Done!" << std::endl;
0455 }