diff --git a/src/ipglasma.cpp b/src/ipglasma.cpp
index 77781e6e350882e17c2925195f94c96af625389c..44b6c79c352be3174e1e2c728f83d37c8b7da5a8 100644
--- a/src/ipglasma.cpp
+++ b/src/ipglasma.cpp
@@ -14,6 +14,7 @@
#include <cmath>
#include <complex>
#include <numeric>
+#include <omp.h>
typedef unsigned int uint;
@@ -509,7 +510,7 @@ std::complex<double> IPGlasma::ComplexAmplitude(double xpom, WilsonLine quark, i
double result = amp.real();
//if (result < 0) return 0;
//if (result > 1) return 1;
- //if (isnan(result))
+ if (isnan(result))
{
cerr << "Wilson line trance NaN, antiquark coords " << q2[0] << ", " << q2[1] << "and quark wline" << quark << endl;
exit(1);
@@ -545,6 +546,8 @@ std::complex<double> IPGlasma::SumAmplitudes()
// in GeV, see LoadData
double lattice_spacing = xcoords[1]-xcoords[0];
+ cout << "# Lattice spacing is " << lattice_spacing << " GeV^-1." << endl;
+
//TODO: loop over all values
double l = 1./(100*lattice_spacing); //[GeV] This is the magnitude of the vector ... whose direction is chosen to be in the x direction (1,0) For starters << a (lattice const)
@@ -557,19 +560,19 @@ std::complex<double> IPGlasma::SumAmplitudes()
// quark flavors
- const double elementary_charge = sqrt(4.*M_PI*1./137.036);
+ const double elementary_charge = sqrt(4.*M_PI/137.036);
//TODO: check masses
quark_flavor u_quark;
- u_quark.mass = 2e-3; //GeV
+ u_quark.mass = 0.14; //GeV
u_quark.e_charge = 2./3.*elementary_charge;
quark_flavor d_quark;
- d_quark.mass = 2e-3;
+ d_quark.mass = 0.14; //GeV
d_quark.e_charge = -1./3.*elementary_charge;
quark_flavor c_quark;
- c_quark.mass = 1.;
+ c_quark.mass = 0.14; //GeV
c_quark.e_charge = u_quark.e_charge;
@@ -577,40 +580,66 @@ std::complex<double> IPGlasma::SumAmplitudes()
std::complex<double> sum = 0.0;
//double sum = 0.0;
-
+ cout << "# r (lattice units) WF_u" << endl << std::flush;
+ // bool is_parallel = false;
#pragma omp declare reduction(+:std::complex<double>:omp_out+=omp_in)
//TODO: check that this sums correctly
- #pragma omp parallel for reduction(+:sum)
- //error: user defined reduction not found for ‘sum’
-
+ //#pragma omp parallel for reduction(+:sum)
for(int i =0; i < L; i += 1){
+
+ //if (i ==0){
+ // if (omp_get_num_threads() > 1) {
+ // #pragma omp critical
+ // { is_parallel = true;}
+
+ //}
+ //if (is_parallel) { std::cout << "# OpenMP is running with " << omp_get_num_threads() << " threads!" << std::endl;}
+ //else { std::cout << "# OpenMP is not enabled, or running with 1 thread." << std::endl;}
+ // }
+
+
for(int j =0; j < L; j += 1){
int q1[2] = {i,j};
- WilsonLine quark = GetWilsonLine(WilsonLineCoordinate(i, j));
- double r1 = DistanceToOrigin(q1);
+ // WilsonLine quark = GetWilsonLine(WilsonLineCoordinate(i, j));
+ // double r1 = DistanceToOrigin(q1);
for(int v =0; v < L; v +=1){
- for(int w =0; w < L; w +=1){
-
+ for(int w =0; w < 100; w +=1){
+ double r = w*lattice_spacing;
+ if (v>0 || j>0 || i>0) {break;}
+ // check that the dipole size is not zero
+ if (i == v && j == w) { continue;}
+
+
int q2[2] = {v,w};
- double r2 = DistanceToOrigin(q2);
+ // double r2 = DistanceToOrigin(q2);
// phase factor vectors
double diffX[2] = {X(i) - X(v), Y(j) - Y(w)};
+ //cout << "# diffX is " << diffX << endl;
+
double sumZX[2] = {z0*X(i) + z1*X(v), z0*Y(j) + z1*Y(w)};
-
- // position for wave function
- double r = sqrt(pow(X(i)-X(v),2) + pow(Y(j)-Y(w),2));
+ //cout << "# sumZX is " << sumZX << endl;
+
+ //TODO: is this correct?
+ //TODO: remove case where r=0
+ // dipole size for wave function
+ //double r = sqrt(pow(X(i)-X(v),2) + pow(Y(j)-Y(w),2));
+ //cout << "# r in WF is " << r << endl;
// wave functions
double WF_u = WFGammaqqbarLongitudinal(r, z0, Q, u_quark);
+ //cout << "# WF_u is " << WF_u << endl;
+
+ cout << r << " " << WF_u << endl << std::flush;
+ continue;
double WF_d = WFGammaqqbarLongitudinal(r, z0, Q, d_quark);
double WF_c = WFGammaqqbarLongitudinal(r, z0, Q, c_quark);
- sum += (WF_u + WF_d + WF_c) * ComplexAmplitude(0, quark, q2)
- * exp( -i * ( std::inner_product(std::begin(vecL), std::end(vecL), std::begin(diffX), 0)
- + std::inner_product(std::begin(vecDelta), std::end(vecDelta), std::begin(sumZX), 0) ) );
+ // sum += (WF_u + WF_d + WF_c) * ComplexAmplitude(0, quark, q2)
+ // * exp( -i * ( std::inner_product(std::begin(vecL), std::end(vecL), std::begin(diffX), 0)
+ // + std::inner_product(std::begin(vecDelta), std::end(vecDelta), std::begin(sumZX), 0) ) );
//TODO: Finish eq (12) from notes
// Finish amplitude
@@ -624,27 +653,27 @@ std::complex<double> IPGlasma::SumAmplitudes()
}
}
}
- cout << "# total amplitude = " << sum << endl;
+
+ cout << "# total amplitude = " << endl;
+ cout << sum << endl;
return sum;
}
-// return longitudinal photon polarization wave function
+// longitudinal photon polarization wave function
// source arxiv.org/pdf/hep-ph/0606272
-// quark flavor f: 0=u, 1=d
double IPGlasma::WFGammaqqbarLongitudinal(double r, double z, double Q, const quark_flavor &quark)
{
- //Assuming r is the transverse size of the dipole
+ //r is the transverse size of the dipole
//constants: e_f (quark charge), e, N_C,
//variables: r, Q, z, epsilon
//epsilon = sqrt(z(1-z)Q^2-m^2_f
- double lattice_spacing = xcoords[1]-xcoords[0];
// fine structure constant set to 1/137
- double e = sqrt(4.*M_PI*1./137.036);
+ double e = sqrt(4.*M_PI/137.036);
- double epsilon = sqrt(z*(1-z)*Q*Q * pow(quark.mass, 2));
-
+ //double epsilon = sqrt(z*(1-z)*Q*Q + pow(quark.mass, 2));
+ double epsilon = 1;
return (quark.e_charge * e * sqrt(NC) * 2*Q*Q*z*(1-z) * std::cyl_bessel_k(0, epsilon*r)*0.5/M_PI);
}