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

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

/* int SimpleSum(int fx, int fy) { */
/*   return TMath::Power((fx + fy),2); */
/* } */

double fPi = TMath::Pi();


double DgamLn(double z);

/* int fZACheck = SimpleSum(3,7); */

double fTerm = 0.0;
double fReg = 1.0;
double fXu = 0;
double fXl = 0;
double fXt = 0;
double myT = -4.54226372679858090e-002;
bool WeissPrint = false;


void PionPlusReg(double s, double t_gev, double qsq, double iterm, double ireg, double& xu, double& xl, double& xt){

  double fGamPi = 0;
  double fGamRh = 0;
  
  if(WeissPrint)
    cout << "Inside pipsrg" << endl;
  if(WeissPrint)
    cout << "s = " << s << "  t = " << t_gev << "  qq = " << qsq << "  iterm = " << iterm << "  ireg = " << ireg << endl;

  double fS0 = 1.0, fAn = 0.94, fApi = 0.14, fArh = 0.77, fGpnn = 1.0, fGrpg = 0.103;
  double fGrnn = 3.4, fGkr = 6.1, fAlem = (1.0/137.0), fAv = 0.77, fConv = 389.0;
  double fI1 = -1.0,  fI2 = -1.0,  fI3 = -1.0;
  double fPP1 = 0.0, fPP2 = 0.0, fPP3 = 0.0;
  if ( iterm == 0 ) { fI1 = 1.0;  fI2 = 1.0;  fI3 = 1.0;  }
  if ( iterm == 1 ) { fI1 = 1.0;  fI2 = 0.0;  fI3 = 0.0;  }
  if ( iterm == 2 ) { fI1 = 0.0;  fI2 = 1.0;  fI3 = 0.0;  }
  if ( iterm == 3 ) { fI1 = 0.0;  fI2 = 0.0;  fI3 = 1.0;  }

  
  if ( ireg == 0.0 ){
    // Point particle propagators
    fPP1 = 1.0/TMath::Power(( t_gev - fApi*fApi),2);
    fPP2 = 1.0/TMath::Power(( t_gev - fArh*fApi),2);
    fPP3 = 1.0/( ( t_gev - fApi*fApi ) * ( t_gev - fArh * fArh) );
  }
  
  if ( ireg == 1.0 ){
    // Regge propagators
    if(WeissPrint)
      cout << "Regge propagators" << endl;
    double fAlPi0 = -0.70 * ( 0.14 * 0.14 );
    double fAlPi1 = 0.70;
    double fAlRh0 = 0.55;
    double fAlRh1 = 0.80;
    double fAlPi = fAlPi0 + fAlPi1 * t_gev;
    double fAlRh = fAlRh0 + fAlRh1 * t_gev;

    /* cout << "alpi = " << fAlPi << "  alrh = " << fAlRh << endl; */

    fAlPi = fAlPi0 + fAlPi1 * t_gev;

    /* cout << setw(8) << "alpi" << setw(12) << fAlPi << endl; */

    if ( (1.0 - fAlPi) < 0 )  
      fGamPi = TMath::Exp(DgamLn(1.0 - fAlPi));

    if ( (1.0 - fAlPi) >= 0 ) {
      fGamPi = TMath::Exp(DgamLn(1.0 - fAlPi));
    }


    if ( (1.0 - fAlRh) < 0 )  
      fGamRh = TMath::Exp(DgamLn(1.0 - fAlRh));
    if ( (1.0 - fAlRh) >= 0 ) {
      fGamRh = TMath::Exp(DgamLn(1.0 - fAlRh));
    }

    /* cout << setw(8) << "alpi" << setw(12) << fAlPi << setw(8) << "alpi1" << setw(12) << fAlPi1  */
    /* 	 << setw(8) << "gampi" << setw(12) << fGamPi  */
    /* 	 << endl;  */

    double fSr  = s/fS0;
    fPP1 = TMath::Power(fSr,(2.0*fAlPi))     * TMath::Power((fAlPi1/fAlPi),2) * TMath::Power(fGamPi,2);
    fPP2 = TMath::Power(fSr,(2.0*fAlRh-2))   * TMath::Power(fAlRh1,2) * TMath::Power(fGamRh,2);

    /* if(WeissPrint) */
    /* cout << setw(8) << "sr"    << setw(12) << fSr    << setw(8) << "alpi"  << setw(12) << fAlPi  */
    /* 	 << setw(8) << "alrh"  << setw(12) << fAlRh  << setw(8) << "alpi1" << setw(12) << fAlPi1   */
    /* 	 << setw(8) << "alrh1" << setw(12) << fAlRh1 << setw(8) << "gampi" << setw(12) << fGamPi  */
    /* 	 << setw(8) << "gamrh" << setw(12) << fGamRh << endl; */

    if(WeissPrint)
      cout << "TMath::Power(fSr,(fAlPi+fAlRh-1)) = " << TMath::Power(fSr,(fAlPi+fAlRh-1)) << endl;
    if(WeissPrint)
      cout << "TMath::Cos((fPi/180)*(fAlPi - fAlRh)) = " << TMath::Cos(fPi*(fAlPi - fAlRh)) << endl;

    // fPP3 = TMath::Power(fSr,(fAlPi+fAlRh-1)) * (fAlPi1/fAlPi) * fAlRh1 * fGamPi * fGamRh * TMath::Cos((fPi/180)*(fAlPi - fAlRh));
    // Question: Why do we have a multiplication by pi in cos()?
    fPP3 = TMath::Power(fSr,(fAlPi+fAlRh-1)) * (fAlPi1/fAlPi) * fAlRh1 * fGamPi * fGamRh * TMath::Cos(fPi*(fAlPi - fAlRh));
  }

  // common auxiliary variable

  double fTT = t_gev - fApi*fApi;
  double fTTQQ = t_gev - fApi * fApi - qsq;

  /* cout << setw(8) << "grpg" << setw(12) << fGrpg << setw(8) << "api" << setw(12) << fApi  */
  /*      << setw(8) << "grnn" << setw(12) << fGrnn << setw(8) << "pp1" << setw(12) << fPP1 << endl; */

  // Pion exchange
  double fFac1 = ((2.0*fGpnn*fGpnn)/(fApi*fApi))*4.0*fAn*fAn*fPP1;
  double fwxx1 = fFac1 * ( 4.0* fApi*fApi*t_gev + fTT*fTT);
  double fwyy1 = fFac1 * fTT * fTT;
  double fwzz1 = fFac1 * qsq * (-1.0*t_gev);

  // Rho exchange


  double fFac2 = ( 2.0 * ( fGrpg * fGrpg ) / ( fApi * fApi ) * fGrnn * fGrnn ) * fPP2;

  if(WeissPrint)
    cout << "fac2 = " << fFac2 << endl;

  double fwxx2 = 0.0;

  if(WeissPrint)
    cout << "fac2 = " << fFac2 << "   gkr = " << fGkr << "   t = " << t_gev 
	 << "   an = " << fAn << "   t = " << t_gev << "   s^2 = " << s * s << endl;

  double fwyy2 = -fFac2 * (1 - ( ( fGkr * fGkr * t_gev ) / ( 4.0 * fAn * fAn ) ) ) * t_gev * s * s;

  double fwzz2 = fFac2 * (1 + fGkr) * (1 + fGkr) * t_gev * t_gev * qsq;

  // Pi-Rho interference

  /* if(WeissPrint) */
  /* cout << setw(8) << "gpnn" << setw(12) << fGpnn << setw(8) << "api"  << setw(12) << fApi */
  /*      << setw(8) << "grpg" << setw(12) << fGrpg << setw(8) << "grnn" << setw(12) << fGrnn */
  /*      << setw(8) << "pp3"  << setw(12) << fPP3 << endl; */

  double fFac3 = ( 2 * fGpnn / fApi ) * ( ( fGrpg / fApi ) * fGrnn ) * fPP3;
  double fwxx3 = 0.0;
  
  /* if(WeissPrint) */
  /* cout << setw(8) << "fac3" << setw(12) << fFac3 << setw(8) << "gkr" << setw(12) << fGkr  */
  /*      << setw(8) << "an"   << setw(12) << fAn   << setw(8) << "tt"  << setw(12) << fTT  */
  /*      << setw(8) << "t"    << setw(12) << t_gev << setw(8) << "s"   << setw(12) << s  */
  /*      << setw(8) << "qq"   << setw(12) << qsq  */
  /*      << endl; */

  double fwyy3 = -fFac3 * fGkr * 2 * fTT * t_gev * s;
  double fwzz3 = ( fFac3 * ( 1 + fGkr ) * 8 * fAn *fAn * fTT * t_gev ) / ( s * qsq);

  // Add contributions
  double fwxx = fI1 * fwxx1 + fI2 * fwxx2 + fI3 * fwxx3;

  /* if(WeissPrint) */
  /*   cout << "  i1 = " << fI1 << "   wyy1 = " << fwyy1  */
  /* 	 << "  i2 = " << fI2 << "   wyy2 = " << fwyy2  */
  /* 	 << "  i3 = " << fI3 << "   wyy3 = " << fwyy3  << endl;     */

  double fwyy = fI1 * fwyy1 + fI2 * fwyy2 + fI3 * fwyy3;

  if(WeissPrint)
    cout << "wyy = " << fwyy << endl;

  /* cout << "  i1 = " << fI1 << "   wzz1 = " << fwzz1  */
  /*      << "  i2 = " << fI2 << "   wzz2 = " << fwzz2  */
  /*      << "  i3 = " << fI3 << "   wzz3 = " << fwzz3  << endl; */

  double fwzz = fI1 * fwzz1 + fI2 * fwzz2 + fI3 * fwzz3;

  // Include vmd form factor
  double fVmd = ( fAv * fAv ) / ( qsq + fAv * fAv);
  fwxx =  fwxx * fVmd * fVmd;

  /* cout << "wyy = " << fwyy << "   fvmd = " << fVmd << endl; */

  fwyy =  fwyy * fVmd * fVmd;

  /* cout << "wzz = " << fwzz << "   fvmd = " << fVmd << endl;   */

  fwzz =  fwzz * fVmd * fVmd;
  
  // Differential cross sections. Keep nucleon mass in phton flux factor
  double fxFac = ( fAlem/(4.0 * TMath::Power(( s - fAn * fAn),2) ) ) * fConv;

  /* cout << "xfac = " << fxFac << "   wxx = " << fwxx << "   wyy = " << fwyy << endl; */

  xu = fxFac * 0.5 * ( fwxx + fwyy );
  xt = fxFac * 0.5 * ( fwxx - fwyy );
  xl = fxFac * fwzz;

  /* cout << setw(8) << "xu" << setw(12) << xu << setw(8) << "xl" << setw(12) << xl  */
  /*      << endl; */

}

double DgamLn(double z){

  /*
**** a double precision routine ****
dgamln computes the natural log of the gamma function for z.gt.0.  the asymptotic expansion is 
used to generate values greater than zmin which are adjusted by the recursion g(z+1)=z*g(z) for z.le.zmin.  
the function was made as portable as possible by computing zmin from the number of base 10 digits in a word, 
rln=max(-alog10(r1mach(4)),0.5e-18) limited to 18 digits of (relative) accuracy. since integer arguments 
are common, a table look up on 100 values is used for speed of execution.
description of arguments
input      z is d0uble precision
z      - argument, z.gt.0.0d0
output      dgamln is double precision
dgamln  - natural log of the gamma function at z.ne.0.0d0
ierr    - error flag
ierr=0, normal return, computation completed
ierr=1, z.le.0.0d0,    no computation
  */

  /* cout << setw(8) << "z" << setw(12) << z << endl; */

  int nz, i;
  double con,   fln,   fz,    rln,  s,     tlg,     trm,  tst,     t1,      wdtol;
  double zdmy,  zinc,  zm,   zmin,  zp,      zsq,  d1mach,  dgamln,  l1m;
  double ierr,  i1m,   k,     mz,   i1mach;
  double fGLn[100] = {  0.00000000000000000e+00,     0.00000000000000000e+00,  6.93147180559945309e-01,     1.79175946922805500e+00,
			3.17805383034794562e+00,     4.78749174278204599e+00,  6.57925121201010100e+00,     8.52516136106541430e+00,
			1.06046029027452502e+01,     1.28018274800814696e+01,  1.51044125730755153e+01,     1.75023078458738858e+01,
			1.99872144956618861e+01,     2.25521638531234229e+01,  2.51912211827386815e+01,     2.78992713838408916e+01,
			3.06718601060806728e+01,     3.35050734501368889e+01,  3.63954452080330536e+01,     3.93398841871994940e+01,
			4.23356164607534850e+01,     4.53801388984769080e+01,  4.84711813518352239e+01,     5.16066755677643736e+01,
			5.47847293981123192e+01,     5.80036052229805199e+01,  6.12617017610020020e+01,     6.45575386270063311e+01,
			6.78897431371815350e+01,     7.12570389671680090e+01,  7.46582363488301644e+01,     7.80922235533153106e+01,
			8.15579594561150372e+01,     8.50544670175815174e+01,  8.85808275421976788e+01,     9.21361756036870925e+01,
			9.57196945421432025e+01,     9.93306124547874269e+01,  1.02968198614513813e+02,     1.06631760260643459e+02,
			1.10320639714757395e+02,     1.14034211781461703e+02,  1.17771881399745072e+02,     1.21533081515438634e+02,
			1.25317271149356895e+02,     1.29123933639127215e+02,  1.32952575035616310e+02,     1.36802722637326368e+02,
			1.40673923648234259e+02,     1.44565743946344886e+02,  1.48477766951773032e+02,     1.52409592584497358e+02,
			1.56360836303078785e+02,     1.60331128216630907e+02,  1.64320112263195181e+02,     1.68327445448427652e+02,
			1.72352797139162802e+02,     1.76395848406997352e+02,  1.80456291417543771e+02,     1.84533828861449491e+02,
			1.88628173423671591e+02,     1.92739047287844902e+02,  1.96866181672889994e+02,     2.01009316399281527e+02,
			2.05168199482641199e+02,     2.09342586752536836e+02,  2.13532241494563261e+02,     2.17736934113954227e+02,
			2.21956441819130334e+02,     2.26190548323727593e+02,  2.30439043565776952e+02,     2.34701723442818268e+02,
			2.38978389561834323e+02,     2.43268849002982714e+02,  2.47572914096186884e+02,     2.51890402209723194e+02,
			2.56221135550009525e+02,     2.60564940971863209e+02,  2.64921649798552801e+02,     2.69291097651019823e+02,
			2.73673124285693704e+02,     2.78067573440366143e+02,  2.82474292687630396e+02,     2.86893133295426994e+02,
			2.91323950094270308e+02,     2.95766601350760624e+02,  3.00220948647014132e+02,     3.04686856765668715e+02,
			3.09164193580146922e+02,     3.13652829949879062e+02,  3.18152639620209327e+02,     3.22663499126726177e+02,
			3.27185287703775217e+02,     3.31717887196928473e+02,  3.36261181979198477e+02,     3.40815058870799018e+02,
			3.45379407062266854e+02,     3.49954118040770237e+02,  3.54539085519440809e+02,     3.59134205369575399e+02 };
  double fCf[22] = {    8.33333333333333333e-02,    -2.77777777777777778e-03,  7.93650793650793651e-04,    -5.95238095238095238e-04,
			8.41750841750841751e-04,    -1.91752691752691753e-03,  6.41025641025641026e-03,    -2.95506535947712418e-02,
			1.79644372368830573e-01,    -1.39243221690590112e+00,  1.34028640441683920e+01,    -1.56848284626002017e+02,
			2.19310333333333333e+03,    -3.61087712537249894e+04,  6.91472268851313067e+05,    -1.52382215394074162e+07,
			3.82900751391414141e+08,    -1.08822660357843911e+10,  3.47320283765002252e+11,    -1.23696021422692745e+13,
			4.88788064793079335e+14,    -2.13203339609193739e+16};
  
  con =  1.83787706640934548;
  double f1Dmach[5] = {exp(-307), exp(307),2.22*exp(-16),4.44*exp(-16),0.30103};
  double fI1Mach[16] = {5,6,5,0,32,4,2,31,2147483647,2,24,-127,127,52,-1023,1023};

  // cout << "z = " << z << endl;
  if ( z < 0 ){
    dgamln = f1Dmach[1];
    ierr=1;
    return dgamln;
  }

  if ( z > 0 && z < 101 ){
    nz = z;
    fz = z - nz;
    // cout << "fz = " << fz << endl;
    if ( fz < 0 && nz < 100 ){
      dgamln = fGLn[nz-1];      
      // cout << "check 1: dgamln = " << dgamln << endl;
      return dgamln;
    }
  }
  
  if ( z > 101 || fz > 0 || nz > 100 ){

    wdtol = f1Dmach[3];
    wdtol = TMath::Max(wdtol,0.5*TMath::Exp(-18));
    l1m = fI1Mach[13];
    rln = f1Dmach[4]*l1m;

    /* cout << setw(8) << "rln" << setw(12) << rln << endl; */

    fln = TMath::Min(rln,20.00);
    fln = TMath::Max(fln,3.0);

    /* cout << setw(8) << "fln" << setw(12) << fln << endl; */

    fln = fln - 3.0;
    zm = 1.8000 + 0.3875*fln;

    /* cout << setw(8) << "zm" << setw(12) << zm << endl; */

    mz = floor(zm + 1);

    /* cout << setw(8) << "mz" << setw(12) << mz << endl; */

    zmin = mz;
    zdmy = z;
    zinc = 0.0;
    if (z < zmin){

      /* cout << setw(8) << "zmin" << setw(12) << zmin << setw(8) << "nz" << setw(12) << nz  */
      /* 	   << setw(8) << "zinc" << setw(12) << zinc  */
      /* 	   << endl;       */

      zinc = zmin - nz;
      zdmy = z + zinc;
    }

    /* cout << setw(8) << "zdmy" << setw(12) << zdmy << endl; */

    zp = 1.0/zdmy;
    t1 = fCf[0]*zp;

    // cout << "cf = " << fCf[0] << "  zp = " << zp << endl;

    s = t1;
    if ( zp >= wdtol ){
      zsq = zp*zp;
      tst = t1*wdtol;
      for (int l = 2; l < 23; l++ ){
	zp = zp*zsq;
        trm = fCf[l]*zp;
	if ( TMath::Abs(trm) <= tst )
	  s = s + trm;
      }
    }

    if ( zinc == 0 ){
      tlg = TMath::Log(z);
      dgamln = z*(tlg-1.0) + 0.50*(con-tlg) + s;
      // cout << "check: dgamln = " << dgamln << endl;
      // cout << "check 2: dgamln = " << dgamln << endl;
      return dgamln;
    }

    zp = 1.0;
    nz =  zinc;
    
    for (int m = 1; m < nz+1; m++ ){
      zp = zp*(z+(m-1));
    }
    tlg = TMath::Log(zdmy);

    /* cout << setw(8) << "zdmy" << setw(12) << zdmy << setw(8) << "tlg"     << setw(12) << tlg */
    /* 	 << setw(8) << "zp"   << setw(12) << zp   << setw(8) << "log(zp)" << setw(12) << log(zp) */
    /* 	 << setw(8) << "con"  << setw(12) << con  << setw(8) << "s"       << setw(12) << s  */
    /* 	 << endl; */

    dgamln = zdmy*(tlg-1.0) - log(zp) + 0.5*(con-tlg) + s;
    // cout << "m = " << m 
    // 	 << "   zdmy = " << zdmy
    // 	 << "   tlg = " << tlg 
    // 	 << "   zp = " << zp 
    // 	 << "   log(zp) = " << log(zp)
    // 	 << "   con = " << con
    // 	 << "   s = " << s << endl;

    return dgamln;
  }

  return 0;
}