ipglasma.cpp

/*
 * Diffraction at sub-nucleon scale
 * Dipole amplitude for a IPglasma nucleus
 * Heikki Mäntysaari <mantysaari@bnl.gov>, 2015
 */


#include "ipglasma.hpp"
#include <string>
#include <gsl/gsl_randist.h>
#include <fstream>
#include <sstream>
#include <cstdlib>
#include <cmath>
#include <complex>
#include <numeric>

typedef unsigned int uint;

using std::cout;
using std::cerr;
using std::endl;

const int NC=3;

int FindIndex(double val, std::vector<double> &vec);


/*
 * Calculate dipole amplitude form Wilson lines
 * Search the closest grid point that corresponds to the given quark/antiquark coordinate
 * Then calcualte 1 - 1/Nc Tr U(quark) U^dagger(antiquark)
 */
double IPGlasma::Amplitude(double xpom, double q1[2], double q2[2] )
{
    ApplyPeriodicBoundaryConditions(q1);
    ApplyPeriodicBoundaryConditions(q2);
    
    // Out of grid? Return 0 (probably very large dipole)

    if (q1[0] < xcoords[0] or q1[0] > xcoords[xcoords.size()-1]
        or q1[1] < ycoords[0] or q1[1] > ycoords[ycoords.size()-1]
        or q2[0] < xcoords[0] or q2[0] > xcoords[xcoords.size()-1]
        or q2[1] < ycoords[0] or q2[1] > ycoords[ycoords.size()-1])
    {
        if (periodic_boundary_conditions)
            cerr << "WTF, I'm here..." << endl;
            
        return 0;
        
    } 

    double  r = sqrt( pow(q1[0]-q2[0],2) + pow(q1[1]-q2[1],2));

	// Smaller dipole than the grid
	if (r < std::abs( xcoords[1] - xcoords[0])) 
			return 0;	
	
		
    // First find corresponding grid indeces
    WilsonLine quark = GetWilsonLine(q1[0], q1[1]);
    WilsonLine antiquark = GetWilsonLine(q2[0], q2[1]);
    
    //antiquark = antiquark.HermitianConjugate();
    
    WilsonLine prod;
    try {
        prod = quark.MultiplyByHermitianConjugate(antiquark);
        //prod =  quark*antiquark;
    } catch (...) {