{
    const int N_reactor = 6;
    const int N_detector = 6;
    // reactor core coordinate (x, y, z) in meters
    double x_reactor[N_reactor] = { 43.0, -44.6, 856.0, 792.3,  1143.6, 1076.5 };
    double y_reactor[N_reactor] = { -7.0, 6.9, 830.9, 767.9,  1206.1, 1138.5 };
    double z_reactor[N_reactor] = { -12.0, -12.0, -12.0, -12.0, -12.0, -12.0 };

    // detector coordinate (x, y, z) in meters
    double x_detector[N_detector] = { 94.5, 97.8, 584.1, -254.3, -259.5, -257.3 };
    double y_detector[N_detector] = { 350.2, 345.2, 1216.2, 1892.6, 1889.6, 1897.8 };
    double z_detector[N_detector] = { -20.0, -20.0, -16.6, -15.4, -15.4, -15.4 };

    // Measured signal rate per day (after subtracting background and correcting for efficiency)
    double rate_measured[N_detector] = {721.7, 705.2, 498.8, 63.2, 66.0, 65.6 };
    double rate_error[N_detector] = { 13.0,  12.9, 11.1, 3.7, 3.8, 3.8}; // squared root of the total candidate, normalized by the live time and efficiency

    // Calculate expected flux based on reactor 1/r^2 inverse power law, with arbitrary normalization
    double flux_expected[N_detector] = {0} ;
    double flux_error[N_detector] = {0};
    double distance_eff[N_detector] = {0};

    for (int i=0; i<N_detector; i++) {
        for (int j=0; j<N_reactor; j++) {
            double distance = sqrt( pow((x_detector[i]-x_reactor[j]),2) + pow((y_detector[i]-y_reactor[j]),2) + pow((z_detector[i]-z_reactor[j]),2) );
            if (j!=3) { // excluding the L2 reactor, which was powered off during that period
                flux_expected[i] += 1/distance/distance;
            }
        }
        distance_eff[i] = sqrt(1/flux_expected[i]*(N_reactor-1)); // weighted distance 1/L_eff^2 = sum(1/L_j^2)/sum(j)
    }

    TGraphErrors *g = new TGraphErrors(N_detector, flux_expected, rate_measured, flux_error, rate_error);
    g->Draw("Ap");
    g->SetMarkerStyle(24);
    g->SetMarkerSize(0.6);
    g->SetTitle("");
    g->GetXaxis()->SetTitle("Expected flux (arbitrary unit)");
    g->GetYaxis()->SetTitle("Measured rate (events/day)");
    g->GetYaxis()->SetRangeUser(0, 800);

    // only fit the first 3 data points from EH1 and EH2, because we want to use the fit to predict EH3 rate
    TGraphErrors *g2 = new TGraphErrors(3, flux_expected, rate_measured, flux_error, rate_error); // only use the near detector to predict the far detector
    TF1 *f = new TF1("f", "[0] * x "); 
    g2->Fit(f);
    f->Draw("same");


    // plot the predicted ratio vs distance 
    TCanvas *c2 = new TCanvas;
    double ratio_fit[N_detector] = {0};
    double ratio_error[N_detector] = {0};
    for (int i=0; i<N_detector; i++) {
        double e = f->Eval(flux_expected[i]);
        ratio_fit[i] = rate_measured[i] / e;
        ratio_error[i] = rate_error[i] / e;
    }
    TGraphErrors *g3 = new TGraphErrors(N_detector, distance_eff, ratio_fit, flux_error, ratio_error);
    g3->Draw("Ap");
    g3->SetMarkerStyle(24);
    g3->SetMarkerSize(0.6);
    g3->SetTitle("");
    g3->GetXaxis()->SetTitle("Weighted Distance (m)");
    g3->GetYaxis()->SetTitle("Measured / Expected");
    g3->GetYaxis()->SetRangeUser(0.7, 1.07);
    g3->GetXaxis()->SetRangeUser(0, 2000);
    c2->SetGridx();
    c2->SetGridy();

    cout << "average EH3 ratio: " << (ratio_fit[3] + ratio_fit[4] + ratio_fit[5]) / 3 
     << " +- " << sqrt(pow(ratio_error[3],2) + pow(ratio_error[4],2) + pow(ratio_error[5],2)) / 3 << endl;
    
    double dm2 = 2.4e-3;  // eV^2
    double L = 1.66e3; // m
    double E = 3.5; // MeV
    cout << "theat13 = " << (1- (ratio_fit[3] + ratio_fit[4] + ratio_fit[5]) / 3) / pow(sin(1.27*dm2*L/E),2) << endl;

}