#include <iomanip>

void asymmetry(TString A180_file, TString A080_file) {
  gROOT->Reset();
  gStyle->SetOptStat(0);
  using namespace std;

  const Float_t d2r = TMath::Pi()/180;

  Float_t Wmin = 2.04;//1.05;//2.01;//1.050;//1.155;//2.01;//1.05;//2.0;//0.90;//1.05; //To apply dilution factor in [0.0,3.0] with 200 bins, Wmin should be divided by 0.015 GeV.
  Int_t Wbin = 9;//20;//11;//9;//11;//20;//34;//40;//22;//33;//75; // For 15MeV resoultion deltaW=1.05 for Wbin=70, 30MeV for Wbin=35
  Int_t bin_option = 2; // 1 for 15MeV, 2 for 30MeV
  Float_t Wbinsize = 0.030; //0.010; // 0.015 for 15MeV, 0.030 for 30MeV
  Float_t Wmax = Wmin + Wbin*Wbinsize;

  TFile *A180 = new TFile(A180_file,"READ");
  TFile *A080 = new TFile(A080_file,"READ");

  TH1F *asym_180 = (TH1F*)A180->Get("W_asym_CPbPtdf_all");
  TH1F *asym_080 = (TH1F*)A080->Get("W_asym_CPbPtdf_all");
  TH1F *asym_090 = new TH1F("W_asym_90","A_{-90} = - #frac{(A_{180}cos80 + A_{80})}{sin80}",Wbin,Wmin,Wmax);

  TH1F *asym_180_N = new TH1F("W_asym_180_N","A_{180}/C_{N}",Wbin,Wmin,Wmax);
  TH1F *asym_080_N = new TH1F("W_asym_080_N","A_{80}/C_{N}",Wbin,Wmin,Wmax);
  TH1F *asym_090_N = new TH1F("W_asym_090_N","A_{-90}/C_{N}",Wbin,Wmin,Wmax);

  //asym_090->Add(asym_180,asym_080,-0.17651445,-1.01650622); // with phi = 2.987 and 3.204 -> average 3.0955
  //asym_090->Add(asym_180,asym_080,-0.17193178,-0.99011569); // with phi = 2.856 and 3.220 -> average 3.038
  asym_090->Add(asym_180,asym_080,-1/TMath::Tan(80*d2r),-1/TMath::Sin(80*d2r)); 


  //Loading the nitrogen correction and interpolating it
  Long_t *id,*size,*flags,*mt;
  
  if (gSystem->GetPathInfo("ncorr_bot_para_RSS.dat",id,size,flags,mt)==0) {

    TH1F *C_N = new TH1F("C_N","Nitrogen correction factor",Wbin,Wmin,Wmax);
    
    ifstream C_N_text;
    C_N_text.open("ncorr_bot_para_RSS.dat");
    
    Float_t W_C_N[100];
    Float_t C_N_value[100];
    TString line;
    Int_t i = 0;
    while(!C_N_text.eof()) {
      
      C_N_text >> W_C_N[i] >> C_N_value[i];
      
      //TObjArray *sublist = line.Tokenize("\t ");
      //TObjString *sub = sublist->At(0);
      //W_C_N[i] = atof(sub->GetString().Data());
      
      //TObjString *sub2 = sublist->At(1);
      //C_N_value[i] = atof(sub2->GetString().Data());
      
      i++;
    }
    
    for(Int_t ibin=1;ibin<=Wbin;ibin++) {
      
      Float_t wcen = (C_N->GetXaxis()->GetBinCenter(ibin))*1e3;
      for(Int_t j=0;j<i;j++) {
	if(wcen>=W_C_N[j] && wcen<W_C_N[j+1]) {
	  cout << W_C_N[j] << " " << wcen << " " << W_C_N[j+1] << endl;
	  Float_t slope = (C_N_value[j+1]-C_N_value[j])/(W_C_N[j+1]-W_C_N[j]);
	  Float_t interpol = slope*(wcen-W_C_N[j])+C_N_value[j];
	  C_N->SetBinContent(ibin,interpol);
	  C_N->SetBinError(ibin,0);
	}
      }
    }

    C_N_text.close();
    
    TCanvas *C_N_c = new TCanvas("C_N_c","C_N canvas");
    C_N_c->cd(1); C_N->Draw();

    for(Int_t ibin = 1; ibin<=Wbin; ibin++) {
      Float_t A180_w = asym_180->GetBinContent(ibin);
      Float_t E180_w = asym_180->GetBinError(ibin);
      Float_t A080_w = asym_080->GetBinContent(ibin);
      Float_t E080_w = asym_080->GetBinError(ibin);
      Float_t C_N_w = 1.;//C_N->GetBinContent(ibin); excluding Nitrogen correction
      
      asym_180_N->SetBinContent(ibin,A180_w/C_N_w);
      asym_180_N->SetBinError(ibin,E180_w/C_N_w);
      asym_080_N->SetBinContent(ibin,A080_w/C_N_w);
      asym_080_N->SetBinError(ibin,E080_w/C_N_w);
    }
    asym_090_N->Add(asym_180_N,asym_080_N,-1/TMath::Tan(80*d2r),-1/TMath::Sin(80*d2r));
  }
  else {
    cout << "No nitrogen correction exists!" << endl;
  }

  
  asym_180->SetMaximum(1.0);
  asym_180->SetMinimum(-1.0);
  asym_080->SetMaximum(1.0);
  asym_080->SetMinimum(-1.0);
  asym_090->SetMaximum(1.0);
  asym_090->SetMinimum(-1.0);

  asym_180_N->SetMaximum(1.0);
  asym_180_N->SetMinimum(-1.0);
  asym_080_N->SetMaximum(1.0);
  asym_080_N->SetMinimum(-1.0);
  asym_090_N->SetMaximum(1.0);
  asym_090_N->SetMinimum(-1.0);

  asym_180->SetTitle("A_{180} without radiative correction");
  asym_080->SetTitle("A_{80} without radiative correction");
  

  TCanvas *c1 = new TCanvas("c1","canvas1"); 
  c1->SetGrid();asym_090->Draw();

  TCanvas *c2 = new TCanvas("c2","canvas2"); 
  c2->SetGrid();asym_180->Draw();

  TCanvas *c3 = new TCanvas("c3","canvas3"); 
  c3->SetGrid();asym_080->Draw();

  TCanvas *c4 = new TCanvas("c4","canvas4"); 
  c4->SetGrid();asym_090_N->Draw();

  TCanvas *c5 = new TCanvas("c5","canvas5"); 
  c5->SetGrid();asym_180_N->Draw();

  TCanvas *c6 = new TCanvas("c6","canvas6"); 
  c6->SetGrid();asym_080_N->Draw();

  //Writing to text file for radiative correction
  //APAR_IN, EPAR_IN, APER_IN, EPER_IN, A1_FRW, E1_FRW, A2_FRW, E2_FRW, Q2_KIN_IN, XBJ_KIN_IN, EPS_KIN_IN, HSCA1PA_IN, HSCA1PE_IN, HSCA2PE_IN, HSNA1_IN, HSNA2_IN, R_FRW_IN, F1_FRW_IN, FL_FRW, FITA1_FRW, FITG1UNK, F2_FRW_IN, W_KIN_IN, HSTHETA_IN, G1_FRW, EG1_FRW, G2_FRW, EG2_FRW, D2_FRW, ED2_FRW
  
  TH1F *W_W_par = (TH1F*)A180->Get("W_W");
  TH1F *Ei_W_par = (TH1F*)A180->Get("Ei_W");
  TH1F *Ef_W_par = (TH1F*)A180->Get("Ef_W");
  TH1F *theta_W_par = (TH1F*)A180->Get("theta_W");
  TH1F *phi_W_par = (TH1F*)A180->Get("phi_W");

  TH1F *W_W_per = (TH1F*)A080->Get("W_W");
  TH1F *Ei_W_per = (TH1F*)A080->Get("Ei_W");
  TH1F *Ef_W_per = (TH1F*)A080->Get("Ef_W");
  TH1F *theta_W_per = (TH1F*)A080->Get("theta_W");
  TH1F *phi_W_per = (TH1F*)A080->Get("phi_W");

  ofstream rcfile;

  if(bin_option==1) rcfile.open("SANE_noRC_1.txt");
  if(bin_option==2) rcfile.open("SANE_noRC_2.txt");
  
  Kin_var *k = new Kin_var(0.0,0.0,0.0,0.0);
  Kin_var *l = new Kin_var(0.0,0.0,0.0,0.0);

  for(Int_t ibin = 1; ibin<=Wbin; ibin++) {
    Float_t W_mid = Wmin + (ibin-0.5)*Wbinsize;
    Float_t W_avg_par = W_W_par->GetBinContent(ibin);
    Float_t Ei_avg_par = Ei_W_par->GetBinContent(ibin);
    Float_t Ef_avg_par = Ef_W_par->GetBinContent(ibin);
    Float_t theta_avg_par = theta_W_par->GetBinContent(ibin);
    Float_t phi_avg_par = phi_W_par->GetBinContent(ibin);

    Float_t W_avg_per = W_W_per->GetBinContent(ibin);
    Float_t Ei_avg_per = Ei_W_per->GetBinContent(ibin);
    Float_t Ef_avg_per = Ef_W_per->GetBinContent(ibin);
    Float_t theta_avg_per = theta_W_per->GetBinContent(ibin);
    Float_t phi_avg_per = phi_W_per->GetBinContent(ibin);

    k->set(Ei_avg_par,Ef_avg_par,theta_avg_par,phi_avg_par);
    l->set(Ei_avg_per,Ef_avg_per,theta_avg_per,phi_avg_per);

    //Float_t W_avgp = W_Wp->GetBinContent(ibin);
    //Float_t Ei_avgp = Ei_Wp->GetBinContent(ibin);
    //Float_t Ef_avgp = Ef_Wp->GetBinContent(ibin);
    //Float_t theta_avgp = theta_Wp->GetBinContent(ibin);

    //l->set(Ei_avgp,Ef_avgp,theta_avgp);

    Float_t Apar = asym_180->GetBinContent(ibin);
    Float_t Epar = asym_180->GetBinError(ibin);
    Float_t Aper = asym_080->GetBinContent(ibin);
    Float_t Eper = asym_080->GetBinError(ibin);

    //Output in GeV and its square units, angle in degree
    rcfile << scientific;
    rcfile << setprecision(5);
    //Older file part
    rcfile << W_mid << "\t";
    rcfile << Apar << "\t";
    rcfile << Epar << "\t";
    rcfile << Aper << "\t";
    rcfile << Eper << "\t";
    rcfile << "0.0\t0.0\t0.0\t0.0\t";
    rcfile << k->Q2()/1e6 << "\t";
    rcfile << k->x_bj() << "\t";
    rcfile << k->epsilon() << "\t";
    rcfile << k->hsca1pa()/1e3 << "\t";
    rcfile << k->hsca1pe()/1e3 << "\t";
    rcfile << k->hsca2pe() << "\t";
    rcfile << k->hsna1()*1e3 << "\t";
    rcfile << k->hsna2() << "\t";
    rcfile << "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t";
    rcfile << W_avg_par << "\t";
    rcfile << theta_avg_par << "\t";
    rcfile << "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t";

    //New part in case of different kinematics, with MeV and degree
    rcfile << Ei_avg_par << "\t";
    rcfile << Ef_avg_par << "\t";
    rcfile << theta_avg_par << "\t";
    rcfile << phi_avg_par << "\t";

    rcfile << Ei_avg_per << "\t";
    rcfile << Ef_avg_per << "\t";
    rcfile << theta_avg_per << "\t";
    rcfile << phi_avg_per << "\t";
    
    rcfile << W_avg_par*1e3 << "\t";
    rcfile << k->Q2() << "\t";
    rcfile << k->x_bj() << "\t";
    rcfile << k->epsilon() << "\t";

    rcfile << W_avg_per*1e3 << "\t";
    rcfile << l->Q2() << "\t";
    rcfile << l->x_bj() << "\t";
    rcfile << l->epsilon() << "\n";
  }

  rcfile.close();

}


// Kinematic variables
class Kin_var {
private:
  Float_t sEi;
  Float_t sEf;
  Float_t stheta;
  Float_t stheta_degree;
  Float_t sphi;
  Float_t sphi_degree;
  Float_t sQ2;
  Float_t M;
public:
  Kin_var(Float_t, Float_t, Float_t, Float_t);
  ~Kin_var();
  void set(Float_t, Float_t, Float_t, Float_t);
  void summary();
  Float_t Ei() { return sEi; }
  Float_t Ef() { return sEf; }
  Float_t theta() { return stheta_degree; }
  Float_t phi() { return sphi_degree; }
  Float_t nu() { return sEi-sEf; }
  Float_t Q2() { return sQ2; }
  Float_t W2();
  Float_t x_bj();
  Float_t gamma();
  Float_t epsilon();
  Float_t hsca1pa();
  Float_t hsca1pe();
  Float_t hsca2pe();
  Float_t hsna1();
  Float_t hsna2();
};

Kin_var::Kin_var(Float_t Ei,Float_t Ef, Float_t theta, Float_t phi) {
  set(Ei,Ef,theta,phi);
}

Kin_var::~Kin_var() { }

void Kin_var::set(Float_t Ei,Float_t Ef, Float_t theta, Float_t phi) {
  // Energy in MeV, theta in degree, but stheta in radian
  sEi = Ei;
  sEf = Ef;
  // Ebeam loss and Ep loss are 2.425 MeV and 2.668 MeV for 73014
  // Ebeam loss and Ep loss are 2.425 MeV and 2.650 MeV for 72790
  stheta = theta*TMath::Pi()/180;
  stheta_degree = theta;
  sphi = phi*TMath::Pi()/180;
  sphi_degree = phi;
  sQ2 = 4*sEi*sEf*TMath::Sin(stheta/2)*TMath::Sin(stheta/2);
  M = 938.272;
}

void Kin_var::summary() {
  cout << "Summary of Kinematic variables" << endl;
  cout << "Ei = " << sEi << endl;
  cout << "Ef = " << sEf << endl;
  cout << "theta = " << stheta << endl;
  cout << "phi = " << sphi << endl;
  cout << "nu = " << nu() << endl;
  cout << "Q2 = " << Q2() << endl;
  cout << "W2 = " << W2() << endl;
  cout << "x_bj = " << x_bj() << endl;
  cout << "gamma = " << gamma() << endl;
}

Float_t Kin_var::W2() {
  return M*M + 2*M*nu() - sQ2;
}

Float_t Kin_var::x_bj() {
  return sQ2/(2*M*nu());
}

Float_t Kin_var::gamma() {
  return 2*M*x_bj()/sqrt(sQ2);
}

Float_t Kin_var::epsilon() {
  return 1/(1+2*(1+1/gamma()/gamma())*TMath::Tan(stheta/2)*TMath::Tan(stheta/2));
}

//# Col 55: HSCA1PA = E-EP*COS(TH)  // MeV              
//# Col 56: HSCA1PE = EP*SIN(TH)/COS(PHI)  // MeV         
//# Col 57: HSCA2PE = (E-EP*COS(TH) )/(EP*SIN(TH)COS(PHI))  // dimless
//# Col 58: HSNA1   = 1/(E+EP)  // MeV^-1
//# Col 59: HSNA2   = SQRT(Q2)/(2E)  // dimless

Float_t Kin_var::hsca1pa() {
  return sEi - sEf*TMath::Cos(stheta);
}

Float_t Kin_var::hsca1pe() {
  return sEf*TMath::Sin(stheta)/TMath::Cos(sphi);
}

Float_t Kin_var::hsca2pe() {
  return (sEi - sEf*TMath::Cos(stheta))/(sEf*TMath::Sin(stheta)*TMath::Cos(sphi));
}

Float_t Kin_var::hsna1() {
  return 1/(sEi + sEf);
}

Float_t Kin_var::hsna2() {
  return sqrt(sQ2)/(2*sEi);
}