//
//  plot_theta13_proj_wreactor.C
////////////////////////////////////////////////////
//
//  To plot chisq as function of theta13 from GLoBES output
//
///////////////////////////////////////////////////
//  Created by S. Cao, cvson@utexas.edu
//  Modified by Ngoc Tran April, 2018
{
   auto c1 = new TCanvas();
   //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);
   //TH2 to fill the plots
   //this is a little tricky. Think about it.

   //This is no of lines in the dat file

   //############## for 3% constraint on sin^2(2th13)  #################//
   Int_t TotalLine = (th13_steps+1);
   double *achisq = new double[TotalLine];
   double *achisq2 = new double[TotalLine];
   double *atheta13 = new double[TotalLine];

   // fprintf(outfile, "%g %g %g\n", this_th13*180.0/M_PI, this_delta*180.0/M_PI, res);
   ifstream cfile("theta13_proj_reactor_6.dat");
   double MinNLL=1e6;//to get the minimum value of chisq
   double MinNLL2=1e6;
   float th13_val, chisq_val, chisq_val2;
   for(int iline=0; iline<TotalLine; iline++){
      cfile >> th13_val >> chisq_val >> chisq_val2;
      //if(iline%10==0) cout<<" at th13 "<<th13_val<<" dcp "<<dcp_val<<" chisq "<<chisq_val<<endl;
      if(chisq_val2<MinNLL2) MinNLL2 = chisq_val2;
      if(chisq_val<MinNLL) MinNLL = chisq_val;
      //to fill 2D plots of chisq
      atheta13[iline] = th13_val;
      achisq[iline] = chisq_val;
      achisq2[iline] = chisq_val2;
   }
   cout <<"The mean value of chisq is "<<MinNLL<<" and2 "<<MinNLL2<<endl; //should be close to zero

    //subtracting zero
    for(int iline=0; iline<TotalLine; iline++){
        achisq[iline]  -= MinNLL;
        achisq2[iline] -= MinNLL2;
    }

   //TGraph
   //TGraph *pchisq = new TGraph(TotalLine,atheta13,achisq);
   TGraph *pchisq2 = new TGraph(TotalLine,atheta13,achisq2);

   new TCanvas;
   pchisq2->Draw("AL");
   pchisq2->SetTitle("");
   pchisq2->SetLineWidth(2);
   pchisq2->GetXaxis()->SetTitle("sin^{2}2#theta_{13}");
   pchisq2->GetYaxis()->SetTitle("#chi^{2}");
   pchisq2->SetMaximum(2);
   pchisq2->GetXaxis()->SetRangeUser(0.09, 0.11);
   pchisq2->SetLineColor(4);
   pchisq2->GetXaxis()->CenterTitle();
   pchisq2->GetYaxis()->CenterTitle();

   /* Draw horizontal lines   */
   c1->Update();
   TLine *l1 = new TLine(c1->GetUxmin(), 1, c1->GetUxmax(), 1);

   l1->SetLineStyle(2);
   l1->Draw();

   double XL = atheta13[0];
   double XH = atheta13[TotalLine-1];
   printf("%g %g\n",XL,XH);
   double dx = (XH-XL)/100000.;
   double y  = 1.;
   double YL, YH;
   TLine *l = new TLine();
   l->SetLineStyle(2);
   for (double x=XL; x<XH; x = x+dx) {
      YL = pchisq2->Eval(x);
      YH = pchisq2->Eval(x+dx);
      if (YL<=y) {
         if (YH>=y) {
            printf("x = %g\n", x);
            l->DrawLine(x,0.,x,y);
         }
      }
      if (YL>=y) {
         if (YH<=y) {
            printf("x = %g\n", x);
            l->DrawLine(x,0.,x,y);
         }
      }
   }

   TLegend * pleg2 = new TLegend(0.4,0.7,0.6,0.85);
   pleg2 -> AddEntry( pchisq2 , "6\% precision" ,"l");

   pleg2 -> SetFillColor(0);
   pleg2 -> SetBorderSize(0);
   pleg2 ->SetTextSize(16);
   pleg2 ->SetTextFont(43);
   pleg2 -> Draw();

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

   TLegend * pleg1 = new TLegend(0.45,0.45,0.55,0.49);
   pleg1 -> AddEntry( l1 , "1#sigma C.L." ,"");
   pleg1 -> SetFillColor(0);
   pleg1 -> SetBorderSize(0);
   pleg1 ->SetTextSize(16);
   pleg1 ->SetTextFont(43);
   pleg1 -> Draw();

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

}