{
    //copy from th13delta.c, change th13_steps and delta_steps to Int_t
    /* Scan the th13-delta plane */
    double this_th13, this_delta;
    double th13_lower  = 0.05;
    double th13_upper  = 0.12;
    //you should increase No of steps to have finer histogram, change in the th13delta.c, recompile and run to produce th13delta.dat again
    Int_t th13_steps  = 360;//15 is too small

    double th13_halfstep = (th13_upper-th13_lower)/(th13_steps*2);

    Int_t TotalLineb = (th13_steps+1);
    double *achisqb = new double[TotalLineb];
    double *achisq2b = new double[TotalLineb];
    double *atheta13b = new double[TotalLineb];

    ifstream cfile("theta13_proj_reactor_current.dat");
    double MinNLLb  = 1e6;//to get the minimum value of chisq
    double MinNLL2b = 1e6;
    float th13_valb, chisq_valb, chisq_val2b;
    for (int ilineb=0; ilineb<TotalLineb; ilineb++){
        cfile >> th13_valb >> chisq_valb >> chisq_val2b;
        if(chisq_val2b<MinNLL2b) MinNLL2b = chisq_val2b;
        if(chisq_valb<MinNLLb) MinNLLb = chisq_valb;
        atheta13b[ilineb] = th13_valb;
        achisqb[ilineb]   = chisq_valb;
        achisq2b[ilineb]  = chisq_val2b;
    }
    cout <<"The mean value of chisq is "<<MinNLLb<<" and2 "<<MinNLL2b<<endl; //should be close to zero

    for(int ilineb=0; ilineb<TotalLineb; ilineb++){
        achisqb[ilineb]  -= MinNLLb;
        achisq2b[ilineb] -= MinNLL2b;
    }

    TGraph *pchisq2  = new TGraph(TotalLineb,atheta13b,achisq2b);
    TGraph *pchisq2b = new TGraph(TotalLineb,atheta13b,achisq2b);

    pchisq2b->Draw("AL");
    pchisq2b->SetTitle("");
    pchisq2b->SetLineWidth(2);
    pchisq2b->GetXaxis()->SetTitle("sin^{2}2#theta_{13}");
    pchisq2b->GetYaxis()->SetTitle("#chi^{2}");
    pchisq2b->SetMaximum(2);
    pchisq2b->GetXaxis()->SetRangeUser(0.075, 0.095);
    pchisq2b->SetLineColor(2);
    pchisq2b->GetXaxis()->CenterTitle();
    pchisq2b->GetYaxis()->CenterTitle();

    pchisq2->SetLineWidth(2);
    pchisq2->SetLineColor(4);
    pchisq2->Draw("*");

    TLegend * pleg1 = new TLegend(0.75,0.62,0.8,0.75);
    pleg1 -> AddEntry( "" , "sin^{2}#theta_{23} = 0.5" ,"");
    pleg1 -> AddEntry("" , "#delta_{CP} = -#pi/2" ,"");
    pleg1 -> SetFillColor(0);
    pleg1 -> SetBorderSize(0);
    pleg1 ->SetTextSize(16);
    pleg1 ->SetTextFont(43);
    pleg1 -> Draw();

    gPad->Print("2020_theta13_proj_reactor.pdf");
}