void Higgs2ggDiscovery()
{
    double mmin = 110;
    double mmax = 150;
    double step = 0.5;
    double sigma = 1.5;   // detector resolution (GeV)
    bool useLikelihood = true;

    TFile *f = new TFile("mgg.root");
    TH1F *h = (TH1F*) f->Get("h_mgg");

    h->Rebin(5);
    h->SetMarkerSize(1);
    h->SetMarkerStyle(20);

    if (!h) {
        std::cerr << "Error: histogram is null!" << std::endl;
        return;
    }

    // ----------------------------
    // Background-only model
    // ----------------------------
    TF1* fb = new TF1("fb", "expo(0)", mmin, mmax);

    // ----------------------------
    // Signal + background model
    // params: [0]=amp, [1]=mean, [2]=sigma, [3]=exp norm, [4]=exp slope
    // ----------------------------
    TF1* fsb = new TF1("fsb", "gaus(0) + expo(3)", mmin, mmax);

    std::vector<double> masses;
    std::vector<double> p0_vals;
    std::vector<double> Z_vals;

    // Fit options
    TString fitOpt = "RQ0";
    if (useLikelihood) fitOpt = "RLQ0";  // likelihood fit

    // ----------------------------
    // Scan loop ove mass
    // ----------------------------
    for (double mH = mmin; mH <= mmax; mH += step) {

        // ---- Background-only fit
        h->Fit(fb, fitOpt);
        double chi2_b = fb->GetChisquare();

        // ---- Signal + background fit
        fsb->SetParameters(10, mH, sigma, 1, -0.05); // reasonable seeds
        fsb->FixParameter(1, mH);    // fix mean
        fsb->FixParameter(2, sigma); // fix width

        h->Fit(fsb, fitOpt);
        double chi2_sb = fsb->GetChisquare();

        // ---- Test statistic
        double q0 = chi2_b - chi2_sb;
        if (q0 < 0) q0 = 0;

        // ---- significance and local p-value
        double Z = sqrt(q0);
        double p0 = 1 -  ROOT::Math::normal_cdf(Z);

        masses.push_back(mH);
        p0_vals.push_back(p0);
        Z_vals.push_back(Z);

        // Printout
        std::cout << "mH = " << mH
                  << "  q0 = " << q0
                  << "  Z = " << Z
                  << "  p0 = " << p0 << std::endl;
    }

    // ----------------------------
    // Build graph
    // ----------------------------
    TGraph* g_p0 = new TGraph(masses.size());
    for (size_t i = 0; i < masses.size(); i++) {
        g_p0->SetPoint(i, masses[i], p0_vals[i]);
    }

    // ----------------------------
    // Plot
    // ----------------------------
    TCanvas* c = new TCanvas("c_p0", "Local p0 scan", 800, 600);

    g_p0->SetTitle(";m_{#gamma#gamma} [GeV];Local p_{0}");
    g_p0->SetLineWidth(2);
    g_p0->SetMarkerStyle(20);
    g_p0->Draw("ALP");

    gPad->SetLogy();

    // ----------------------------
    // Reference lines (significance)
    // ----------------------------
    std::vector<double> sigmas = {1, 2, 3, 4, 5};

    for (double s : sigmas) {
        double p = 0.5 * TMath::Erfc(s / sqrt(2));

        TLine* line = new TLine(mmin, p, mmax, p);
        line->SetLineStyle(2);
        line->Draw("same");

        TLatex* latex = new TLatex(mmax + 0.5, p, Form("%g#sigma", s));
        latex->SetTextSize(0.03);
        latex->Draw();
    }

    c->Update();

    // ----------------------------
    // Final fitted mgg plot
    // ----------------------------
    TCanvas* c2 = new TCanvas("c_2", "mass", 800, 600);
    fsb->SetParameter(1, 125); // fix mean
    fsb->SetParameter(2, 1.5); // fix width
    h->Fit(fsb);
    h->Draw("e");
}