#include<iostream>
#include<fstream>
using namespace std;

const int sets_num=1000;

//world data from PDG

const double dC1_world=0.0068;
const double dC2_world=0.06;
const double corr12_world=-0.22;//-0.22;
const double dC1_ofs_world=0.;  // -0.0030;
const double dC2_ofs_world=0.;  // -0.0346;

/*more /group/eic/users/aemmert/2022/polTest/pol85/11_12400_ciq_fitSummary.txt
==================================================
11GeV,12400,Ciq,b_unfixed,stat+full_syst+pdf
==================================================
-0.712943 0.0335298
-0.140684 0.0404602
0.00123607 0.00447485
-0.0215336 0.0373782
==================================================
0.00112424 -0.00133772 4.31014e-05 -0.00018349 
-0.00133772 0.00163703 -6.29858e-05 0.000368755 
4.31014e-05 -6.29858e-05 2.00243e-05 -0.00015065 
-0.00018349 0.000368755 -0.00015065 0.00139713 
==================================================
*/

const double dC1_solid=0.03353;
const double dC2_solid=0.04046;
const double corr12_solid=-0.986;
const double dC1_ofs_solid=0.;
const double dC2_ofs_solid=0.;

/* more /group/eic/users/aemmert/2022/polTest/pol85/11_12400_ciq_fitSummary.txt
==================================================
11GeV,12400,Ciq(P2),b_unfixed,stat+full_syst+pdf
==================================================
-0.714194 0.00235418
-0.139196 0.00729006
0.00118811 0.00428716
-0.0213335 0.0369768
==================================================
5.54215e-06 -6.59418e-06 2.12287e-07 -9.02248e-07 
-6.59418e-06 5.3145e-05 -1.19544e-05 0.000151488 
2.12287e-07 -1.19544e-05 1.83797e-05 -0.00014365 
-9.02248e-07 0.000151488 -0.00014365 0.00136729 
==================================================
*/

const double dC1_solid_p2=0.00235;
const double dC2_solid_p2=0.00729;
const double corr12_solid_p2=-0.38; // <- from the 2x2 sub matrix above
const double dC1_ofs_solid_p2=0;
const double dC2_ofs_solid_p2=0;


//edge
//const double dC1_ofs=0.005;
//const double dC2_ofs=-0.01;

//const double theta0=+acos(-1)/4.*3./4.;
//const double theta0=0.;
const double PI=acos(-1.);
const double constF=0.24622*sqrt(8.*PI); // in TeV

void plot_masslimit1(){

	double const_sin2thw=0.23119;
	double dC1,dC2,corr12,dC1_ofs,dC2_ofs;
	double S11,S12,S22,axisa,axisb,theta0;
	  
	double C1[sets_num],C2[sets_num];
	double Lx[sets_num], Ly[sets_num];
	double Lx_alt[sets_num], Ly_alt[sets_num];
	double C1_solid[sets_num],C2_solid[sets_num];
	double Lx_solid[sets_num], Ly_solid[sets_num];
	double Lx_alt_solid[sets_num], Ly_alt_solid[sets_num];
	double C1_solid_p2[sets_num],C2_solid_p2[sets_num];
	double Lx_solid_p2[sets_num], Ly_solid_p2[sets_num];
	double Lx_alt_solid_p2[sets_num], Ly_alt_solid_p2[sets_num];
	double tmp, tmp_th, d, theta;

	//------------------existing world data from PDG------------------
	dC1=dC1_world;
	dC2=dC2_world;
	corr12=corr12_world;
	dC1_ofs=dC1_ofs_world;
	dC2_ofs=dC2_ofs_world;

	S11=dC1*dC1; S12=corr12*dC1*dC2; S22=dC2*dC2;

	// following https://cookierobotics.com/007/
	
	axisa=sqrt((S11+S22)/2+sqrt(((S11-S22)*(S11-S22)/4)+S12*S12));
	axisb=sqrt((S11+S22)/2-sqrt(((S11-S22)*(S11-S22)/4)+S12*S12));
	theta0=atan2(axisa*axisa-S11,S12);
	  
	for(int i=0;i<sets_num;i++){
	  theta=2.*PI/sets_num * i;
	  // generate the ellipse and rotate to desired orientation
	  C1[i]=(axisa*cos(theta) * cos(theta0) - axisb*sin(theta) * sin(theta0)) + dC1_ofs;
	  C2[i]=(axisa*cos(theta) * sin(theta0) + axisb*sin(theta) * cos(theta0)) + dC2_ofs;
	  d=sqrt(C1[i]*C1[i]+C2[i]*C2[i]);
	  tmp=constF/sqrt(d);
	  // theta+theta0 angle is the direction of BSM physics
	  Lx[i]=tmp*cos(theta+theta0);
	  Ly[i]=tmp*sin(theta+theta0);

	  d=sqrt(C1[i]*C1[i]*cos(theta)*cos(theta) +C2[i]*C2[i]*sin(theta)*sin(theta));
	  //	  d=sqrt(C1[i]*C1[i]*cos(theta)*cos(theta) +C2[i]*C2[i]*sin(theta)*sin(theta)+2 * C1[i] * C2[i]*sin(theta)*cos(theta)*corr12);
	  tmp=constF/sqrt(d);
	  Lx_alt[i]=tmp*cos(theta);
	  Ly_alt[i]=tmp*sin(theta);
	  
	}
	
	//------------------SoLID alone------------------
	dC1=dC1_solid;
	dC2=dC2_solid;
	corr12=corr12_solid;
	dC1_ofs=dC1_ofs_solid;
	dC2_ofs=dC2_ofs_solid;


	S11=dC1*dC1; S12=corr12*dC1*dC2; S22=dC2*dC2;

	// following https://cookierobotics.com/007/
	
	axisa=sqrt((S11+S22)/2+sqrt(((S11-S22)*(S11-S22)/4)+S12*S12));
	axisb=sqrt((S11+S22)/2-sqrt(((S11-S22)*(S11-S22)/4)+S12*S12));
	theta0=atan2(axisa*axisa-S11,S12);
	  
	for(int i=0;i<sets_num;i++){
	  theta=2.*PI/sets_num * i;
	  // generate the ellipse and rotate to desired orientation
	  // theta angle is the direction of BSM physics, irrelevant of the uncertainty of C1,C2
	  C1_solid[i]=(axisa*cos(theta) * cos(theta0) - axisb*sin(theta) * sin(theta0)) + dC1_ofs;
	  C2_solid[i]=(axisa*cos(theta) * sin(theta0) + axisb*sin(theta) * cos(theta0)) + dC2_ofs;
	  d=sqrt(C1_solid[i]*C1_solid[i]+C2_solid[i]*C2_solid[i]);
	  tmp=constF/sqrt(d);
	  // theta+theta0 angle is the direction of BSM physics
	  Lx_solid[i]=tmp*cos(theta+theta0);
	  Ly_solid[i]=tmp*sin(theta+theta0);

	  d=sqrt(C1_solid[i]*C1_solid[i]*cos(theta)*cos(theta) +C2_solid[i]*C2_solid[i]*sin(theta)*sin(theta));
	  //	  d=sqrt(C1_solid[i]*C1_solid[i]*cos(theta)*cos(theta) +C2_solid[i]*C2_solid[i]*sin(theta)*sin(theta)+2 * C1_solid[i] * C2_solid[i]*sin(theta)*cos(theta)*corr12);
	  tmp=constF/sqrt(d);
	  Lx_alt_solid[i]=tmp*cos(theta);
	  Ly_alt_solid[i]=tmp*sin(theta);
	  
	}
	  
	//------------------SoLID and P2, no other PVDIS data ------------------
	dC1=dC1_solid_p2;
	dC2=dC2_solid_p2;
	corr12=corr12_solid_p2;
	dC1_ofs=dC1_ofs_solid_p2;
	dC2_ofs=dC2_ofs_solid_p2;


	S11=dC1*dC1; S12=corr12*dC1*dC2; S22=dC2*dC2;

	// following https://cookierobotics.com/007/
	
	axisa=sqrt((S11+S22)/2+sqrt(((S11-S22)*(S11-S22)/4)+S12*S12));
	axisb=sqrt((S11+S22)/2-sqrt(((S11-S22)*(S11-S22)/4)+S12*S12));
	theta0=atan2(axisa*axisa-S11,S12);
	  
	for(int i=0;i<sets_num;i++){
	  theta=2.*PI/sets_num * i;
	  // generate the ellipse and rotate to desired orientation
	  // theta angle is the direction of BSM physics, irrelevant of the uncertainty of C1,C2
	  C1_solid_p2[i]=(axisa*cos(theta) * cos(theta0) - axisb*sin(theta) * sin(theta0)) + dC1_ofs;
	  C2_solid_p2[i]=(axisa*cos(theta) * sin(theta0) + axisb*sin(theta) * cos(theta0)) + dC2_ofs;
	  d=sqrt(C1_solid_p2[i]*C1_solid_p2[i]+C2_solid_p2[i]*C2_solid_p2[i]);
	  tmp=constF/sqrt(d);
	  // theta+theta0 angle is the direction of BSM physics
	  Lx_solid_p2[i]=tmp*cos(theta+theta0);
	  Ly_solid_p2[i]=tmp*sin(theta+theta0);

	  d=sqrt(C1_solid_p2[i]*C1_solid_p2[i]*cos(theta)*cos(theta) + C2_solid_p2[i]*C2_solid_p2[i]*sin(theta)*sin(theta));
	  //	  d=sqrt(C1_solid_p2[i]*C1_solid_p2[i]*cos(theta)*cos(theta) +C2_solid_p2[i]*C2_solid_p2[i]*sin(theta)*sin(theta)+2 * C1_solid_p2[i] * C2_solid_p2[i]*sin(theta)*cos(theta)*corr12);
	  tmp=constF/sqrt(d);
	  Lx_alt_solid_p2[i]=tmp*cos(theta);
	  Ly_alt_solid_p2[i]=tmp*sin(theta);
	  
	}
	  
//------------------abs errors

	TGraph *g_C1C2=new TGraph(sets_num,C1,C2);
	g_C1C2->SetTitle("C2 vs. C1 (SM at 0,0)");
	g_C1C2->GetXaxis()->SetTitle("C1");
	g_C1C2->GetYaxis()->SetTitle("C2");
	TGraph *g_C1C2_solid=new TGraph(sets_num,C1_solid,C2_solid);
	TGraph *g_C1C2_solid_p2=new TGraph(sets_num,C1_solid_p2,C2_solid_p2);
	
	TGraph *g_lambda12=new TGraph(sets_num,Lx,Ly);
	g_lambda12->SetTitle("mass limit (TeV)");
	g_lambda12->GetXaxis()->SetTitle("Lambda1");
	g_lambda12->GetYaxis()->SetTitle("Lambda2");
	TGraph *g_lambda12_solid=new TGraph(sets_num,Lx_solid,Ly_solid);
	g_lambda12_solid->SetTitle("mass limit (TeV)");
	g_lambda12_solid->GetXaxis()->SetTitle("Lambda(C1)");
	g_lambda12_solid->GetYaxis()->SetTitle("Lambda(C2)");
	TGraph *g_lambda12_solid_p2=new TGraph(sets_num,Lx_solid_p2,Ly_solid_p2);
	g_lambda12_solid_p2->SetTitle("mass limit (TeV)");
	g_lambda12_solid_p2->GetXaxis()->SetTitle("Lambda(C1)");
	g_lambda12_solid_p2->GetYaxis()->SetTitle("Lambda(C2)");
	
	TGraph *g_lambda12_alt=new TGraph(sets_num,Lx_alt,Ly_alt);
	TGraph *g_lambda12_alt_solid=new TGraph(sets_num,Lx_alt_solid,Ly_alt_solid);
	TGraph *g_lambda12_alt_solid_p2=new TGraph(sets_num,Lx_alt_solid_p2,Ly_alt_solid_p2);


//---abs
	TCanvas *can1=new TCanvas("can1","can1",700,1200);
	can1->Divide(1,2);
	can1->cd(1);
	g_C1C2->SetMarkerColor(1);
	g_C1C2->SetMarkerStyle(1);
	g_C1C2->SetMarkerSize(2);
	g_C1C2->SetLineColor(1);
	g_C1C2->GetXaxis()->SetLimits(-0.08,0.08);
	g_C1C2->GetYaxis()->SetRangeUser(-0.08,0.08);
	g_C1C2->Draw("APL");
	can1->Update();

	g_C1C2_solid->SetMarkerColor(4);
	g_C1C2_solid->SetMarkerStyle(1);
	g_C1C2_solid->SetMarkerSize(2);
	g_C1C2_solid->SetLineColor(4);
	g_C1C2_solid->Draw("same");

	g_C1C2_solid_p2->SetMarkerColor(2);
	g_C1C2_solid_p2->SetMarkerStyle(1);
	g_C1C2_solid_p2->SetMarkerSize(2);
	g_C1C2_solid_p2->SetLineColor(2);
	g_C1C2_solid_p2->Draw("same");
        gPad->SetGridx(1);
        gPad->SetGridy(1);

	can1->Update();
	
	can1->cd(2);
	
	g_lambda12_solid_p2->SetMarkerColor(2);
	g_lambda12_solid_p2->SetMarkerStyle(1);
	g_lambda12_solid_p2->SetMarkerSize(2);
	g_lambda12_solid_p2->SetLineColor(2);
	g_lambda12_solid_p2->Draw("APL");
	g_lambda12_solid_p2->GetXaxis()->SetLimits(-40,40);
	g_lambda12_solid_p2->GetYaxis()->SetRangeUser(-40,40);
	
	g_lambda12_solid->SetMarkerColor(4);
	g_lambda12_solid->SetMarkerStyle(1);
	g_lambda12_solid->SetMarkerSize(2);
	g_lambda12_solid->SetLineColor(4);
	g_lambda12_solid->Draw("same");
	g_lambda12_solid->GetXaxis()->SetLimits(-40,40);
	g_lambda12_solid->GetYaxis()->SetRangeUser(-40,40);
	
	
	g_lambda12->SetMarkerColor(1);
	g_lambda12->SetMarkerStyle(1);
	g_lambda12->SetMarkerSize(2);
	g_lambda12->SetLineColor(1);
	g_lambda12->Draw("same"); // same
	
        gPad->SetGridx(1);
        gPad->SetGridy(1);
	can1->Update();
	
	can1->Print("masslimit1.pdf");

	/*
	TCanvas *can2=new TCanvas("can2","can2",600,600);
	can2->cd(1);
	
	g_lambda12_alt_solid->SetMarkerColor(2);
	g_lambda12_alt_solid->SetMarkerStyle(1);
	g_lambda12_alt_solid->SetMarkerSize(2);
	g_lambda12_alt_solid->SetLineColor(2);
	g_lambda12_alt_solid->Draw("APL");
	g_lambda12_alt_solid->GetXaxis()->SetLimits(-40,40);
	g_lambda12_alt_solid->GetYaxis()->SetRangeUser(-40,40);
	
	
	g_lambda12_alt->SetMarkerColor(1);
	g_lambda12_alt->SetMarkerStyle(1);
	g_lambda12_alt->SetMarkerSize(2);
	g_lambda12_alt->SetLineColor(1);
	g_lambda12_alt->Draw("same"); // same
	
        gPad->SetGridx(1);
        gPad->SetGridy(1);
	can2->Update();
	
	can2->Print("masslimit1_alt.pdf");
	*/
	
}