#include "TCanvas.h"
#include "TGraph.h"
#include "TLine.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TPad.h"
#include "TMath.h"
#include "TAxis.h"

#include <cmath>
#include <iostream>
 
 
  TCanvas *c1 = new TCanvas("c1", "", 200, 10, 700, 500);
    double myEne;
    double myDmsq23 = 2.5e-3;//in eV^2
    double myDmsq21 = 7.6e-5;//in eV^2
    double myBaseL;// T2K 295km, Nova 
    double myTheta12 = asin(sqrt(0.8704))/2; //in radian 
    double myTheta23 = asin(sqrt(0.5)); //in radian
    double myTheta13 = asin(sqrt(0.085))/2; //in radian
    double myDelta = -TMath::Pi()/2; //in radian
    
// ... first .... all helper functions ...
  double prob_T2K(double Ene, double dmsq_23, double dmsq_21, double baseL, double theta_12, double theta_23, double theta_13, double delta) {
  
  double probMu2ea =  
-16.*std::sin(theta_12)*std::sin(theta_13)*std::sin(theta_23)*std::cos(theta_12)*std::pow(std::cos(theta_13),2)*std::cos(theta_23)*std::sin(delta)*std::pow(std::sin(1.27*dmsq_23*baseL/Ene),2)*std::sin(1.27*dmsq_21*baseL/Ene);

  return probMu2ea;
}
 double prob_Nova(double Ene, double dmsq_23, double dmsq_21, double baseL, double theta_12, double theta_23, double theta_13, double delta) {
  
  double probMu2eb =  
-16.*std::sin(theta_12)*std::sin(theta_13)*std::sin(theta_23)*std::cos(theta_12)*std::pow(std::cos(theta_13),2)*std::cos(theta_23)*std::sin(delta)*std::pow(std::sin(1.27*dmsq_23*baseL/Ene),2)*std::sin(1.27*dmsq_21*baseL/Ene);

  return probMu2eb;
}
// ... then .... the main macro function (with the same name as the file)  ...
  void CPviolation_vac_compare(){

  if (myEne == 0.6 and myBaseL == 295) {
  //set number of steps you want to calculate
  const Int_t nEnStep = 200;
  Double_t pEne[nEnStep], pProb[nEnStep];
  
  // set maximal energy to plot
  Double_t EnergyMax = 2.0;
  double myProbMu2ea = prob_T2K(myEne, myDmsq23, myDmsq21, myBaseL, myTheta12, myTheta23, myTheta13, myDelta);
 
  //Graph of oscillation prob. as function of energy
  //Note:formula is not work at zero energy
  
  for (Int_t istep=1; istep<=nEnStep; ++istep) {
    //energy value at each step
    pEne[istep-1] = EnergyMax * ((Double_t)istep) / ((Double_t)nEnStep);
    //calculated osc. prob. at each step
    pProb[istep-1] = prob_T2K(pEne[istep-1], myDmsq23, myDmsq21, myBaseL, myTheta12, myTheta23, myTheta13, myDelta);
  }
  //create new Graph based on the array
   TGraph *pGraphMu2ea = new TGraph(nEnStep, pEne, pProb);
   pGraphMu2ea->SetTitle(""); 
   pGraphMu2ea->GetXaxis()->SetTitle("Neutrino Energy [GeV]");
   pGraphMu2ea->GetYaxis()->SetTitle("P(#nu_{#mu}#rightarrow #nu_{e}) - P(#bar{#nu}_{#mu}#rightarrow #bar{#nu}_{e})");
   pGraphMu2ea->GetXaxis()->SetRangeUser(0,4.);
   pGraphMu2ea->GetYaxis()->SetRangeUser(0, 0.1);
   pGraphMu2ea->GetXaxis()->CenterTitle();
   pGraphMu2ea->GetYaxis()->CenterTitle();
   pGraphMu2ea->GetYaxis()->SetTitleOffset(1.2);
   pGraphMu2ea->SetLineWidth(2);
   pGraphMu2ea->SetLineColor(2);
   pGraphMu2ea->Draw("AL");
   }
   if ( myEne == 2 and myBaseL == 810){
   //set number of steps you want to calculate
  const Int_t nEnStep = 200;
  Double_t pEne[nEnStep], pProb[nEnStep];
  
  // set maximal energy to plot
  Double_t EnergyMax = 4.0;
  
  double myProbMu2eb = prob_Nova(myEne, myDmsq23, myDmsq21, myBaseL, myTheta12, myTheta23, myTheta13, myDelta);
 
  //Graph of oscillation prob. as function of energy
  //Note:formula is not work at zero energy
  
  for (Int_t istep=1; istep<=nEnStep; ++istep) {
    //energy value at each step
    pEne[istep-1] = EnergyMax * ((Double_t)istep) / ((Double_t)nEnStep);
    //calculated osc. prob. at each step
    pProb[istep-1] = prob_Nova(pEne[istep-1], myDmsq23, myDmsq21, myBaseL, myTheta12, myTheta23, myTheta13, myDelta);
  }
  //create new Graph based on the array
   TGraph *pGraphMu2eb = new TGraph(nEnStep, pEne, pProb);
   pGraphMu2eb->SetTitle(""); 
   pGraphMu2eb->GetXaxis()->SetTitle("Neutrino Energy [GeV]");
   pGraphMu2eb->GetYaxis()->SetTitle("P(#nu_{#mu}#rightarrow #nu_{e}) - P(#bar{#nu}_{#mu}#rightarrow #bar{#nu}_{e})");
   pGraphMu2eb->GetXaxis()->SetRangeUser(0,4.);
   pGraphMu2eb->GetYaxis()->SetRangeUser(0, 0.1);
   pGraphMu2eb->GetXaxis()->CenterTitle();
   pGraphMu2eb->GetYaxis()->CenterTitle();
   pGraphMu2eb->GetYaxis()->SetTitleOffset(1.2);
   pGraphMu2eb->SetLineWidth(2);
   pGraphMu2eb->SetLineColor(4);
   pGraphMu2eb->Draw("AL");
   }
   
  /* Draw horizontal lines   */
    c1->Update();
    TLine *l1 = new TLine(c1->GetUxmin(), 0.0, c1->GetUxmax(), 0.0);
    l1->SetLineColor(1);
    l1->Draw("L SAME");
   TLegend * pleg = new TLegend(0.7,0.65,0.85,0.87);
    pleg -> AddEntry( pGraphMu2ea , "T2K - CP violation" ,"l");
    pleg -> AddEntry( pGraphMu2eb, "Nova -CP violation" ,"l");
    
    pleg -> SetFillColor(0);
    pleg -> SetBorderSize(0);
    pleg ->SetTextSize(14);
    pleg ->SetTextFont(43);
    pleg -> Draw();

    gStyle->SetOptStat(0);
    TLegend* leg = new TLegend(0.7,0.4,0.85,0.5);
    leg->SetTextSize(0.03);
    leg->SetTextColor(1);
    leg->SetFillColor(0);
    leg->SetBorderSize(0);
    leg ->SetTextSize(14);
    leg ->SetTextFont(43);
    leg->AddEntry(pGraphMu2ea,"sin^{2}2#theta_{12} = 0.8704","");
    leg->AddEntry(pGraphMu2ea,"sin^{2}2#theta_{13} = 0.085","");
    leg->AddEntry(pGraphMu2ea,"sin^{2} #theta_{23} = 0.5","");
    leg->AddEntry(pGraphMu2ea,"|#Deltam^{2}_{21}| = 7.6x10^{-5} eV^{2}/c^{4}","");
    leg->AddEntry(pGraphMu2ea,"|#Deltam^{2}_{32}| = 2.5x10^{-3} eV^{2}/c^{4}","");
    leg->AddEntry(pGraphMu2ea,"#delta_{CP} = -#pi/2","");
    leg->Draw();//

    
  //save the graph
    gPad->Print("CPviolation_Nova_vac");//pdf format

}