/* This file calculates the pion electromagnetic FF : 
   leading order PQCD + k_perp dependence + soft overlap. */

/**************************************************************************/
/*                                                                        */
/* IMPORTANT : all dimensional quantities in this program                 */
/*             are in units GeV                                           */  
/*                                                                        */
/**************************************************************************/

#include <pion_emff.h>

static int choice_coupling_s5;
static double Q_sqr_s5, x_s6, y_s7, sigma_s5;

/* The following function calculates the PQCD prediction for the 
   pion electromagnetic form factor 
   WITH inclusion of transverse momentum dependence */

double pion_emff(int choice_coupling, double Q_sqr)
{
     choice_coupling_s5 = choice_coupling; 
     Q_sqr_s5 = Q_sqr;
     
     return int_gauss5(int1_pion_emff, 0., 1., 30);
}

double int1_pion_emff(double x)
{
     x_s6 = x;
     
     return int_gauss6(int2_pion_emff, 0., 1., 30);
}

double int2_pion_emff(double y)
{
     y_s7 = y;

     return 1./pow(4. * PI, 2.) 
            * int_gauss7(int3_pion_emff, 0.01, 25., 30);
}

/*   b is in units GeV^-1 */

double int3_pion_emff(double b)
{
     int choice_sudakov;
     double Q, scale, scale1, scale2, wavef_x, wavef_y, sudakov_x, sudakov_y;

     Q = sqrt(Q_sqr_s5); 

     if(choice_coupling_s5 == 4)
     { 
       scale = pow(sqrt(x_s6 * y_s7)/sqrt(2.) * Q, 2.); 
     } 
     else
     {
       scale1 = x_s6 * y_s7 * Q_sqr_s5;
       scale2 = pow(0.366,2.); /* average transverse momentum squared */
       if(scale1 > scale2)
	 scale = scale1;
       else 
	 scale = scale2;
     }

     choice_sudakov = 1;

     switch(choice_sudakov)
     {
     case 1 : /* with k_perp dependence but without sudakov factor */
            {
             wavef_x = pion_wavef_impactpar(x_s6, b);
             wavef_y = pion_wavef_impactpar(y_s7, b);
             sudakov_x = 1.;
             sudakov_y = 1.;
             break;
            }
     }


     return b * wavef_x * sudakov_x * wavef_y * sudakov_y
            * 4./3. * (16. * PI * alpha_strong(choice_coupling_s5, scale)) 
            * bessK_0(sqrt(x_s6 * y_s7) * Q * b);
}


/*************************************************************************/
/*************************************************************************/


/* The function pion_emff_soft_overlap() calculates the soft overlap 
   contribution to the pion electromagnetic form factor */

double pion_emff_soft_overlap(double Q_sqr)
{
     double sigma, integ;
     
     sigma = 8. * pow(PI * F_PI,2.);  /* average squared transverse 
                                         quark momentum in pion 
                                         wavefunction in units GeV^2 */
     sigma_s5 = sigma;
     Q_sqr_s5 = Q_sqr;
     integ = int_gauss5(int_pion_emff_soft_overlap, 0.001, .999, 30);
     
     return pow(2. * PI,2.)/sigma * integ;
}

double int_pion_emff_soft_overlap(double x)
{     
     double distr;

     distr = pion_distribution_ampl(x);

     return pow(distr,2.)/(x * (1. - x)) 
            * exp(- (1. - x) * Q_sqr_s5 / (4. * sigma_s5 * x) );
}