/*Include{{{*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <cmath>
#include <string.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
#include <map>
#include <errno.h>
#include <sstream>
//#include <iterator>
//from root
/*#include <TMath.h>*/
//#include <TRandom3.h>
//#include <TLorentzVector.h>
//#include <TMatrix.h>
//#include <TFile.h>
//#include <TTree.h>
//#include <TString.h>
//#include <TStopwatch.h>
/*#include <TRandom2.h>*/
/*}}}*/

/*Constants{{{*/
using namespace::std;
static const double PI=3.14159265358979323846;
static const double ALPHA=1/137.036; //fine-structure constant
static const double Deg2Rad = PI/180.0;
static const double MeV2GeV = 0.001;
static const double AP=ALPHA/PI;
static const double ELECTRON_MASS=0.510998918*MeV2GeV; //GeV
static const double NQe = 1.6*1e-19;
static const double Avogadro = 6.02*1e+23;
static const double M2=ELECTRON_MASS*ELECTRON_MASS*MeV2GeV; //GeV
static const double AMU=931.494043*MeV2GeV; //GeV
static const double PROTON_MASS=938.272029*MeV2GeV; //GeV
static const double DELTA_MASS=1232.*MeV2GeV; //GeV
static const double PI_MASS=139.57*MeV2GeV; //GeV
static const double HBARC=197.326968*MeV2GeV; //hbar*c GeV fm
static const double NBARN_TO_GEV2= 1./3.894 * 1e-5; //nbarn-->GeV^-2
static const double GEV2_TO_NBARN= 3.894 * 1e5; //GeV^-2-->nbarn
static const double NBARN_TO_CM2=1e-33; //nbarn-->CM^2
/*}}}*/

#include "christy_bosted_inelastic.h"
#include "DIS.h"

int main(){
    double E0, Ep, Theta;
    cout <<"---- E0 = (GeV)"<<endl; cin >> E0;
    cout <<"---- Ep (GeV) = "; cin >> Ep;
    cout <<"---- Theta (Degree) = "; cin >> Theta;

    double nu = E0-Ep;
    double Q2 = 4.0*E0*Ep*pow(sin(Theta/2.0),2);
    double xb = Q2 / (2.0*PROTON_MASS*nu);
    double W2 = PROTON_MASS*PROTON_MASS + 2.0*PROTON_MASS*nu - Q2;

    string model = ""; 
    int mod = 0; cout<<"-- Which model (1->F2ALLM, 2->F1F2IN09, 3->CJ or 4->PB) "; cin >> mod;
    if(mod==1) model="F2ALLM";
    else if(mod==2) model="F1F2IN09";
    else if(mod==3) model="CJ";
    else if(mod==4) model="PB";
    else {cerr<<"*** Wrong model#!!!!"<<endl; return -111;}

    DIS* dis = new DIS();
    dis->SetModel(Model.c_str());

    dis->SetKin(E0, Ep, Theta*Deg2Rad);
    //dis->SetKin(Model.c_str(), E0, Ep, Theta*Deg2Rad);

    double F1_p = dis->GetF1("Proton");
    //double F1_p = dis->GetF1(Model.c_str(),"Proton");
    double F2_p = dis->GetF2("Proton");
    //double F2_p = dis->GetF2(Model.c_str(),"Proton");

    double xs_p = dis->Sigma("Proton");
    //double xs_p = dis->Sigma(Model.c_str(),"Proton");
    double xs_n = dis->Sigma("Neutron");
    //double xs_n = dis->Sigma(Model.c_str(),"Neutron");
    double xs_d = dis->Sigma("Deutron");
    //double xs_d = dis->Sigma(Model.c_str(),"Deutron");
    double xs_al = dis->Sigma(27, 13);
    //double xs_al = dis->Sigma(Model.c_str(),27, 13);
}