analysis_INTT/plot/combine_profilelikelihood.C

0001 #include "./plotutil.h"
0002 #include "./util_combine.h"
0003
0004 std::vector<float> v_centralitybin = {0, 3, 6, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70};
0005
0006 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"};
0007 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"};
0008 std::vector<std::string> v_unclabel = {"Tracklet reconstruction selection", "Cluster ADC cut", "Cluster #phi-size cut", "Run segments", "Statistical uncertainty in corrections", "Geometric offset"};
0009 std::vector<std::string> v_plotname = {"dRcut", "clusAdcCut", "clusPhiSizeCut", "segment", "statUnc", "GeoOffset"};
0010 std::vector<std::string> v_color = {"#f2777a", "#6699cc", "#9999cc", "#99cc99", "#888888", "#ffcc66"};
0011 std::vector<float> v_alpha = {0.7, 0.7, 0.6, 0.5, 0.8, 0.7};
0012 std::vector<int> v_marker = {20, 21, 33, 34, 47, 45};
0013 std::vector<double> v_corr = {}; // the same order as v_phobos_histname/v_cms_histname
0014 int Nbins_dNdeta = 27;
0015 float xbinlow_dNdeta = -2.7;
0016 float xbinhigh_dNdeta = 2.7;
0017 double desiredCL = 0.99;
0018 double step = 1E-5;        // Fine resolution for scanning the profile likelihood
0019 bool set_phobosunc = true; // for testing
0020
0021 void setbinattr()
0022 {
0023     TFile *f = new TFile("./systematics/Centrality0to3_Zvtxm10p0to10p0_noasel/finalhists_systematics_Centrality0to3_Zvtxm10p0to10p0_noasel.root", "READ");
0024     TH1 *h = (TH1 *)f->Get("hM_final");
0025     h->SetDirectory(0);
0026     Nbins_dNdeta = h->GetNbinsX();
0027     xbinlow_dNdeta = h->GetXaxis()->GetXmin();
0028     xbinhigh_dNdeta = h->GetXaxis()->GetXmax();
0029     f->Close();
0030 }
0031
0032 TH1 *gethist(std::string fname, std::string histname)
0033 {
0034     TFile *f = new TFile(fname.c_str(), "READ");
0035     TH1 *h;
0036     if (!f->Get(histname.c_str()))
0037     {
0038         std::cout << "Warning: cannot find histogram " << histname << " in file " << fname << std::endl;
0039         std::cout << "Set a hisgoram with all bins = 0" << std::endl;
0040         h = new TH1D(histname.c_str(), histname.c_str(), Nbins_dNdeta, xbinlow_dNdeta, xbinhigh_dNdeta);
0041         for (int i = 1; i <= Nbins_dNdeta; i++)
0042         {
0043             h->SetBinContent(i, 0);
0044         }
0045     }
0046     else
0047     {
0048         h = (TH1 *)f->Get(histname.c_str());
0049     }
0050     h->SetDirectory(0);
0051     f->Close();
0052     return h;
0053 }
0054
0055 void sysunccorr()
0056 {
0057     system("mkdir -p ./systematics/correlation/");
0058
0059     std::vector<TGraph *> v_g_unc;
0060
0061     for (size_t j = 0; j < v_phobos_histname.size(); j++)
0062     {
0063         std::vector<float> v_unc_phobos, v_unc_cms;
0064         v_unc_phobos.clear();
0065         v_unc_cms.clear();
0066
0067         for (size_t i = 0; i < v_centralitybin.size() - 1; i++)
0068         {
0069             std::string centstr = Form("Centrality%dto%d", (int)v_centralitybin[i], (int)v_centralitybin[i + 1]);
0070             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));
0071             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]);
0072
0073             TH1 *h_phobos = gethist(fphobos, v_phobos_histname[j]);
0074             TH1 *h_cms = gethist(fcms, v_cms_histname[j]);
0075
0076             // now loop over the bins with non-zero content
0077             for (int k = 1; k <= h_phobos->GetNbinsX(); k++)
0078             {
0079                 float unc_phobos = h_phobos->GetBinContent(k);
0080                 float unc_cms = h_cms->GetBinContent(k);
0081                 if (unc_cms == 0)
0082                     continue;
0083                 v_unc_phobos.push_back(unc_phobos * 100);
0084                 v_unc_cms.push_back(unc_cms * 100);
0085             }
0086         }
0087
0088         TGraph *g_unc = new TGraph(v_unc_phobos.size(), v_unc_phobos.data(), v_unc_cms.data());
0089         v_corr.push_back(g_unc->GetCorrelationFactor());
0090         v_g_unc.push_back(g_unc);
0091         drawTGraph(g_unc, "PHOBOS approach [%]", "CMS approach [%]", {"Uncertainty: " + v_unclabel[j]}, Form("./systematics/correlation/correlation_%s", v_plotname[j].c_str()));
0092     }
0093
0094     drawTMultiGraph(v_g_unc, "PHOBOS approach [%]", "CMS approach [%]", "./systematics/correlation/correlation_all");
0095 }
0096
0097 //------------------------------------------------------
0098 // For a chi-square distribution with 1 degree of freedom,
0099 // the cumulative distribution function can be expressed as:
0100 //    F(x;1) = erf( sqrt(x/2) )
0101 // where erf is the error function provided in <cmath>.
0102 //------------------------------------------------------
0103 double chi_square_cdf_1dof(double x) { return erf(sqrt(x / 2.0)); }
0104
0105 //------------------------------------------------------
0106 // Inverse chi-square function for 1 degree of freedom using binary search.
0107 // For a given probability p (0 < p < 1), it finds x such that
0108 //    chi_square_cdf_1dof(x) = p.
0109 // 'tol' sets the tolerance for convergence.
0110 //------------------------------------------------------
0111 double chi_square_inv_1dof(double p, double tol = 1e-6)
0112 {
0113     double low = 0.0;
0114     double high = 10.0;
0115     // Expand the upper bound until the CDF is at least p.
0116     while (chi_square_cdf_1dof(high) < p)
0117     {
0118         high *= 2;
0119     }
0120     // Binary search for the value x such that CDF(x) ~ p.
0121     while (high - low > tol)
0122     {
0123         double mid = (low + high) / 2.0;
0124         double cdf_mid = chi_square_cdf_1dof(mid);
0125         if (cdf_mid < p)
0126             low = mid;
0127         else
0128             high = mid;
0129     }
0130     return (low + high) / 2.0;
0131 }
0132
0133 //------------------------------------------------------
0134 // get_threshold_for_confidence()
0135 // For a given desired confidence level (e.g., 0.95 for 95%),
0136 // this function returns the corresponding chi-square threshold
0137 // for one degree of freedom using the inverse CDF function.
0138 //------------------------------------------------------
0139 double get_threshold_for_confidence(double desiredCL)
0140 {
0141     if (desiredCL <= 0.0 || desiredCL >= 1.0)
0142     {
0143         cout << "Invalid confidence level. Please enter a value between 0 and 1.\n";
0144         return -1;
0145     }
0146     return chi_square_inv_1dof(desiredCL);
0147 }
0148
0149 //------------------------------------------------------
0150 // Gaussian likelihood function.
0151 // Computes the likelihood of observing X given a predicted value theta
0152 // and a Gaussian uncertainty sigma.
0153 //------------------------------------------------------
0154 double gaussian_likelihood(double X, double theta, double sigma)
0155 {
0156     double exponent = -0.5 * pow((X - theta) / sigma, 2);
0157     return exp(exponent) / (sqrt(2 * M_PI) * sigma);
0158 }
0159
0160 //------------------------------------------------------
0161 // Profile likelihood function (systematics only).
0162 // This function loops over the available systematic uncertainty sources
0163 // (each element in the vectors represents one source) and computes the likelihood
0164 // for a given theta. In this simplified example, we take the maximum likelihood
0165 // over the systematic configurations.
0166 // Note: There is no statistical uncertainty in this version.
0167 //------------------------------------------------------
0168 double profile_likelihood(double theta,                     //
0169                           double X1,                        //
0170                           double X2,                        //
0171                           const vector<double> &sigma1_sys, //
0172                           const vector<double> &sigma2_sys, //
0173                           const vector<double> &correlation //
0174 )
0175 {
0176     // Combine systematic uncertainties in quadrature.
0177     double V11 = 0.0, V22 = 0.0, V12 = 0.0;
0178     int n = sigma1_sys.size();
0179     for (int i = 0; i < n; i++)
0180     {
0181         V11 += sigma1_sys[i] * sigma1_sys[i];
0182         V22 += sigma2_sys[i] * sigma2_sys[i];
0183         V12 += correlation[i] * sigma1_sys[i] * sigma2_sys[i];
0184     }
0185     // Form the covariance matrix:
0186     // V = [ V11, V12 ]
0187     //     [ V12, V22 ]
0188     double detV = V11 * V22 - V12 * V12;
0189     if (detV <= 0)
0190     {
0191         // In case of numerical issues, return a very small likelihood.
0192         return 1e-300;
0193     }
0194     // Inverse covariance matrix elements:
0195     double invV11 = V22 / detV;
0196     double invV12 = -V12 / detV;
0197     double invV22 = V11 / detV;
0198
0199     // The difference vector: d = [X1 - theta, X2 - theta]^T.
0200     double d1 = X1 - theta;
0201     double d2 = X2 - theta;
0202
0203     // Exponent term: -0.5 * d^T * V^-1 * d.
0204     double exponent = -0.5 * (d1 * d1 * invV11 + 2.0 * d1 * d2 * invV12 + d2 * d2 * invV22);
0205
0206     // Joint likelihood:
0207     double L = exp(exponent) / (2.0 * M_PI * sqrt(detV));
0208     return L;
0209 }
0210
0211 //------------------------------------------------------
0212 // Find the best estimate of theta by scanning a range.
0213 // Returns the value of theta that maximizes the profile likelihood.
0214 //------------------------------------------------------
0215 double find_best_estimate(double X1,                         //
0216                           double X2,                         //
0217                           const vector<double> &sigma1_sys,  //
0218                           const vector<double> &sigma2_sys,  //
0219                           const vector<double> &correlation, //
0220                           double theta_min,                  //
0221                           double theta_max,                  //
0222                           double step                        //
0223 )
0224 {
0225     double best_theta = X1; // initial guess
0226     double max_likelihood = -numeric_limits<double>::infinity();
0227
0228     for (double theta = theta_min; theta <= theta_max; theta += step)
0229     {
0230         double L = profile_likelihood(theta, X1, X2, sigma1_sys, sigma2_sys, correlation);
0231         if (L > max_likelihood)
0232         {
0233             max_likelihood = L;
0234             best_theta = theta;
0235         }
0236     }
0237     return best_theta;
0238 }
0239
0240 //------------------------------------------------------
0241 // Compute the confidence interval using the likelihood ratio method.
0242 // For a given threshold (obtained from the desired confidence level),
0243 // this function scans from the best estimate toward lower and higher theta values
0244 // to find where the likelihood ratio -2ln(L/L_max) crosses the threshold.
0245 // Linear interpolation is used for better accuracy.
0246 //------------------------------------------------------
0247 pair<double, double> compute_confidence_interval(double X1,                         //
0248                                                  double X2,                         //
0249                                                  const vector<double> &sigma1_sys,  //
0250                                                  const vector<double> &sigma2_sys,  //
0251                                                  const vector<double> &correlation, //
0252                                                  double best_theta,                 //
0253                                                  double L_max,                      //
0254                                                  double theta_min,                  //
0255                                                  double theta_max,                  //
0256                                                  double step,                       //
0257                                                  double threshold                   //
0258 )
0259 {
0260     double lower_bound = best_theta;
0261     double upper_bound = best_theta;
0262     bool found_lower = false, found_upper = false;
0263
0264     // Scan downward from best_theta to find the lower bound.
0265     double prev_theta = best_theta;
0266     double prev_delta = 0.0;
0267     for (double theta = best_theta; theta >= theta_min; theta -= step)
0268     {
0269         double L = profile_likelihood(theta, X1, X2, sigma1_sys, sigma2_sys, correlation);
0270         double delta = -2.0 * log(L / L_max);
0271         if (theta == best_theta)
0272         {
0273             prev_delta = delta;
0274             prev_theta = theta;
0275             continue;
0276         }
0277         if (delta >= threshold && prev_delta < threshold)
0278         {
0279             double t = (threshold - prev_delta) / (delta - prev_delta);
0280             lower_bound = prev_theta + t * (theta - prev_theta);
0281             found_lower = true;
0282             break;
0283         }
0284         prev_delta = delta;
0285         prev_theta = theta;
0286     }
0287
0288     // Scan upward from best_theta to find the upper bound.
0289     prev_theta = best_theta;
0290     prev_delta = 0.0;
0291     for (double theta = best_theta; theta <= theta_max; theta += step)
0292     {
0293         double L = profile_likelihood(theta, X1, X2, sigma1_sys, sigma2_sys, correlation);
0294         double delta = -2.0 * log(L / L_max);
0295         if (theta == best_theta)
0296         {
0297             prev_delta = delta;
0298             prev_theta = theta;
0299             continue;
0300         }
0301         if (delta >= threshold && prev_delta < threshold)
0302         {
0303             double t = (threshold - prev_delta) / (delta - prev_delta);
0304             upper_bound = prev_theta + t * (theta - prev_theta);
0305             found_upper = true;
0306             break;
0307         }
0308         prev_delta = delta;
0309         prev_theta = theta;
0310     }
0311
0312     if (!found_lower)
0313         lower_bound = theta_min;
0314     if (!found_upper)
0315         upper_bound = theta_max;
0316
0317     return make_pair(lower_bound, upper_bound);
0318 }
0319
0320 void combine_profilelikelihood()
0321 {
0322     setbinattr();
0323
0324     // First, calculate the correlation between the systematic uncertainties in the PHOBOS and CMS approaches
0325     sysunccorr();
0326
0327     system("mkdir -p ./systematics/combined/");
0328
0329     std::cout << "Correlation in the uncertainty between PHOBOS and CMS approaches:" << std::endl;
0330     for (size_t i = 0; i < v_corr.size(); i++)
0331     {
0332         std::cout << v_unclabel[i] << ": " << v_corr[i] << std::endl;
0333     }
0334
0335     // start the loop over the centrality bins
0336     for (size_t i = 0; i < v_centralitybin.size() - 1; i++)
0337     // for (size_t i = 0; i < 1; i++)
0338     {
0339         std::cout << "Combine the PHOBOS and CMS measurements for centrality bin " << v_centralitybin[i] << " to " << v_centralitybin[i + 1] << std::endl;
0340
0341         std::string centstr = Form("Centrality%dto%d", (int)v_centralitybin[i], (int)v_centralitybin[i + 1]);
0342         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));
0343         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]);
0344
0345         // histograms for the central values
0346         TH1 *h_phobos_X = gethist(fphobos, "h1D_dNdEta_reco");
0347         TH1 *h_cms_X = gethist(fcms, "hM_final");
0348         // histograms for the total uncertainties, for checking
0349         TH1 *h_phobos_totlunc = gethist(fphobos, "h1D_error_Final");
0350         TH1 *h_cms_totlunc = gethist(fcms, "hM_TotalRelUnc");
0351         for (int b = 0; b < h_cms_X->GetNbinsX(); b++)
0352         {
0353             if (h_cms_X->GetBinContent(b) < 0.5)
0354             {
0355                 h_cms_X->SetBinContent(b, 0);
0356                 h_cms_X->SetBinError(b, 0);
0357             }
0358             // set the uncertainty of the CMS measurement to be the same as phobos
0359             if (set_phobosunc)
0360                 h_cms_X->SetBinError(b, h_phobos_totlunc->GetBinContent(b) * h_cms_X->GetBinContent(b));
0361         }
0362
0363         // histograms for the uncertainties
0364         std::vector<TH1 *> v_h_phobos_unc, v_h_cms_unc;
0365         v_h_phobos_unc.clear();
0366         v_h_cms_unc.clear();
0367         for (size_t j = 0; j < v_phobos_histname.size(); j++)
0368         {
0369             v_h_phobos_unc.push_back(gethist(fphobos, v_phobos_histname[j]));
0370             v_h_cms_unc.push_back(gethist(fcms, v_cms_histname[j]));
0371         }
0372
0373         std::vector<double> v_combined_x, v_combined_y, v_combined_xerr, v_combined_yerrup, v_combined_yerrdown;
0374         // loop over the bins
0375         for (int k = 1; k <= h_phobos_X->GetNbinsX(); k++)
0376         {
0377             double X_phobos = h_phobos_X->GetBinContent(k);
0378             double X_cms = h_cms_X->GetBinContent(k);
0379             if (X_phobos < 5E-1 || X_cms < 5E-1)
0380                 continue;
0381
0382             std::vector<double> unc_phobos, unc_cms;
0383             unc_phobos.clear();
0384             unc_cms.clear();
0385             for (size_t j = 0; j < v_phobos_histname.size(); j++)
0386             {
0387                 if (set_phobosunc)
0388                     v_h_cms_unc[j]->SetBinContent(k, v_h_phobos_unc[j]->GetBinContent(k));
0389
0390                 unc_phobos.push_back(v_h_phobos_unc[j]->GetBinContent(k) * X_phobos);
0391                 unc_cms.push_back(v_h_cms_unc[j]->GetBinContent(k) * X_cms);
0392             }
0393
0394             std::cout << "----- Bin " << k << ", eta = " << h_phobos_X->GetBinCenter(k) << " -----" << std::endl;
0395             std::cout << "PHOBOS measurement: " << X_phobos << " +/- " << h_phobos_totlunc->GetBinContent(k) * X_phobos << std::endl;
0396             std::cout << "CMS measurement: " << X_cms << " +/- " << h_cms_totlunc->GetBinContent(k) * X_cms << std::endl;
0397             double threshold = get_threshold_for_confidence(desiredCL);
0398             if (threshold < 0)
0399             {
0400                 std::cout << "Invalid confidence level. Please enter a value between 0 and 1." << std::endl;
0401                 exit(1);
0402             }
0403             // Define the scanning range and step size for theta (parameter of interest, POI)
0404             double theta_min = min(X_phobos, X_cms) - 1.;
0405             double theta_max = max(X_phobos, X_cms) + 1.;
0406
0407             // Find the best estimate by maximizing the profile likelihood
0408             double best_theta = find_best_estimate(X_phobos, X_cms, unc_phobos, unc_cms, v_corr, theta_min, theta_max, step);
0409
0410             // Compute the maximum likelihood at the best estimate
0411             double L_max = profile_likelihood(best_theta, X_phobos, X_cms, unc_phobos, unc_cms, v_corr);
0412
0413             // Compute the confidence interval using the likelihood ratio method
0414             pair<double, double> conf_interval = compute_confidence_interval(X_phobos, X_cms, unc_phobos, unc_cms, v_corr, best_theta, L_max, theta_min, theta_max, step, threshold);
0415
0416             cout << "Combined measurement, best estimate = " << best_theta << " with confidence interval (" << desiredCL * 100 << "%): [" << conf_interval.first << ", " << conf_interval.second << "]" << endl;
0417
0418             // assign the combined measurement to the histogram
0419             v_combined_x.push_back(h_phobos_X->GetBinCenter(k));
0420             v_combined_y.push_back(best_theta);
0421             v_combined_xerr.push_back(h_phobos_X->GetBinWidth(k) * 0.49);
0422             v_combined_yerrup.push_back(conf_interval.second - best_theta);
0423             v_combined_yerrdown.push_back(best_theta - conf_interval.first);
0424         }
0425
0426         // set up the TGraphAssymErrors for the combined measurement
0427         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());
0428
0429         // draw the combined measurement
0430         TCanvas *c = new TCanvas("c", "c", 800, 700);
0431         c->cd();
0432         gPad->SetRightMargin(0.07);
0433         gPad->SetTopMargin(0.07);
0434         gPad->SetLeftMargin(0.17);
0435         h_phobos_X->SetLineColor(kRed);
0436         h_phobos_X->SetLineWidth(0);
0437         h_phobos_X->SetMarkerColor(kRed);
0438         h_phobos_X->SetMarkerStyle(20);
0439         h_phobos_X->SetMarkerSize(1);
0440         h_phobos_X->SetFillColorAlpha(kRed, 0.1);
0441         h_phobos_X->GetXaxis()->SetTitle("#eta");
0442         h_phobos_X->GetYaxis()->SetTitle("dN_{ch}/d#eta");
0443         h_phobos_X->GetYaxis()->SetTitleOffset(1.7);
0444         h_phobos_X->GetYaxis()->SetRangeUser(h_cms_X->GetMinimum(0) * 0.7, h_cms_X->GetMaximum() * 1.3);
0445         h_phobos_X->Draw("E2P");
0446         h_cms_X->SetLineColor(kBlue);
0447         h_cms_X->SetLineWidth(0);
0448         h_cms_X->SetMarkerColor(kBlue);
0449         h_cms_X->SetMarkerStyle(21);
0450         h_cms_X->SetMarkerSize(1);
0451         h_cms_X->SetFillColorAlpha(kBlue, 0.1);
0452         h_cms_X->Draw("E2P SAME");
0453         g_combined->SetLineColor(kBlack);
0454         g_combined->SetLineWidth(1);
0455         g_combined->SetMarkerColor(kBlack);
0456         g_combined->SetMarkerStyle(33);
0457         g_combined->SetMarkerSize(1);
0458         g_combined->SetFillColorAlpha(kBlack, 0.4);
0459         g_combined->Draw("5E P SAME");
0460         TLegend *leg = new TLegend(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.05, gPad->GetLeftMargin() + 0.25, 1 - gPad->GetTopMargin() - 0.25);
0461         leg->SetHeader(Form("Centrality %d-%d%%", (int)v_centralitybin[i], (int)v_centralitybin[i + 1]));
0462         leg->SetTextAlign(kHAlignLeft + kVAlignCenter);
0463         leg->SetTextSize(0.03);
0464         leg->AddEntry(h_phobos_X, "PHOBOS approach", "pf");
0465         leg->AddEntry(h_cms_X, "CMS approach", "pf");
0466         leg->AddEntry(g_combined, Form("Combined (Profile likelihood, C.L=%.3g%%) ", desiredCL * 100), "pf");
0467         leg->SetBorderSize(0);
0468         leg->SetFillStyle(0);
0469         leg->Draw();
0470         gPad->Modified();
0471         gPad->Update();
0472         if (set_phobosunc)
0473         {
0474             c->SaveAs(Form("./systematics/combined/combined_%s_profilelikelihood_setphobosunc.pdf", centstr.c_str()));
0475             c->SaveAs(Form("./systematics/combined/combined_%s_profilelikelihood_setphobosunc.png", centstr.c_str()));
0476         }
0477         else
0478         {
0479             c->SaveAs(Form("./systematics/combined/combined_%s_profilelikelihood.pdf", centstr.c_str()));
0480             c->SaveAs(Form("./systematics/combined/combined_%s_profilelikelihood.png", centstr.c_str()));
0481         }
0482     }
0483 }