#include <iostream>
#include <fstream>
#include <cstdlib>
#include <stdlib.h>
#include <iomanip>
#include <stdio.h>
#include <math.h>
#include <cmath>
#include <string>
#include <sstream>

#include "SigWeiss.h"

using std::setw;
using std::setprecision;
using std::cout;
using std::cin;
using std::endl;
using namespace std;

int weiss( double w, double t, double qsq ) {

  double fSig_Weiss = 0;

  double fT_GeV   = -t;
  double fW_GeV   = w;
  double fQsq_GeV = qsq;

  Double_t fBjorkenX = 0;
  Double_t fAn = 0.94;
  Double_t fS0 = 4.0;
  Double_t fQQ0 = 1.0;
  Double_t fS = fW_GeV * fW_GeV;

  Double_t fXu = 0.0;
  Double_t fXl = 0.0;
  Double_t fXt = 0.0;

  if ( fQsq_GeV < fQQ0 ){ 
    Double_t fTerm = 0.0;
    Double_t fReg = 1.0;
    Double_t fXu0 = 0;
    Double_t fXl0 = 0;
    Double_t fXt0 = 0;

    PionPlusReg(fS, fT_GeV, fQsq_GeV, fTerm, fReg,fXu, fXl, fXt);

  }
  if ( fQsq_GeV >= fQQ0 ){ 
    Double_t fTerm = 0.0;
    Double_t fReg = 1.0;
    Double_t fXu0 = 0;
    Double_t fXl0 = 0;
    Double_t fXt0 = 0;

    fBjorkenX = fQsq_GeV / ( fS - fAn * fAn + fQsq_GeV );
    Double_t fS1 = ( ( 1.0 - fBjorkenX ) / fBjorkenX ) * fQQ0 + fAn * fAn;

    if ( fS1 > fS0 ){
	    	      
      PionPlusReg(fS1, fT_GeV, fQQ0, fTerm, fReg,fXu0, fXl0, fXt0);
      fXu = TMath::Power((fQQ0/fQsq_GeV), 4 ) * fXu0;
      fXl = TMath::Power((fQQ0/fQsq_GeV), 3 ) * fXl0;

    }

    if ( fS1 <= fS0 ){
	    
      Double_t fQQ1 = ( fBjorkenX / ( 1.0 - fBjorkenX ) ) * ( fS0 - fAn * fAn );
      PionPlusReg(fS0, fT_GeV, fQQ1, fTerm, fReg, fXu0, fXl0, fXt0);
      fXu = TMath::Power((fQQ1/fQsq_GeV), 4 ) * fXu0;
      fXl = TMath::Power((fQQ1/fQsq_GeV), 3 ) * fXl0;
	
    }	    
  }

  fT_GeV = -1.0 * fT_GeV;
  // fSig_Weiss = fXu + fEpsilon * fXl;
  fSig_Weiss = fXu + 0.95 * fXl;
  fSig_Weiss = fSig_Weiss / 2.0;

  return fSig_Weiss;
}


int eic() {

  double sigma = 0;


  
}