diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 13bdcabdd9fa02332d4cd5244b24881b552d3566..82452eab4493445dcb805570def7885d001141e0 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -6,6 +6,9 @@ add_executable(
ipglasma.cpp
)
+if(OpenMP_CXX_FOUND)
+ target_link_libraries(wline PUBLIC OpenMP::OpenMP_CXX)
+endif()
target_link_libraries(
@@ -13,4 +16,5 @@ target_link_libraries(
PRIVATE
GSL::gsl
GSL::gslcblas
+ OpenMP::OpenMP_CXX
)
diff --git a/src/ipglasma.cpp b/src/ipglasma.cpp
index ae7bf5af9093a9675d505f5394801e47bba58122..3a3d232572652fb83cc3d78cfd1792e3f42309d5 100644
--- a/src/ipglasma.cpp
+++ b/src/ipglasma.cpp
@@ -12,6 +12,8 @@
#include <sstream>
#include <cstdlib>
#include <cmath>
+#include <complex>
+#include <numeric>
typedef unsigned int uint;
@@ -584,7 +586,7 @@ double IPGlasma::DistanceToOrigin(int q[2]){
//TODO: check
-double IPGlasma::SumAmplitudes(int l)
+double IPGlasma::SumAmplitudes(int w)
{
double sum = 0.0;
int L = xcoords.size();
@@ -595,8 +597,8 @@ double IPGlasma::SumAmplitudes(int l)
int q1[2] = {i,j};
//TODO: pointer instead of copy?
quark = GetWilsonLine(WilsonLineCoordinate(i, j));
- for(int k =0; k < L; k +=1){
- int q2[2] = {k,l};
+ for(int v =0; v < L; v +=1){
+ int q2[2] = {v,w};
double r1 = DistanceToOrigin(q1);
double r2 = DistanceToOrigin(q2);
if (r1 > 1. && r2 > 1.){
@@ -607,9 +609,6 @@ double IPGlasma::SumAmplitudes(int l)
continue;
}
- // Haw w2 = WilsonLine(k,l)
- // Laske w1.MultiplbyHermitianConjugate(w2);
- // sum += amplitude(0, q1, q2, ipglasma);
sum += Amplitude(0, quark, q2);
}
}
@@ -625,31 +624,79 @@ double IPGlasma::SumAmplitudes(int l)
double IPGlasma::SumAmplitudes()
{
int L = xcoords.size();
- double z = 0.5;
- double Q = 0;
-
-
cout << "# Summing with lattice size L = " << L << endl;
- double sum = 0.0;
- WilsonLine quark;
- #pragma omp parallel for reduction(+:sum)
+
+ double z0 = 0.5; //Momentum fractions
+ double z1 = 1-z0;
+ double Q = sqrt(10); //[GeV] Sqrt of photon virtuality
+ double Mx = 1; //[GeV] Invarint mass of system X
+
+ // in GeV, see LoadData
+ double lattice_spacing = xcoords[1]-xcoords[0];
+
+ //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)
+ double vecL[2] = {l, 0.}; //this is the vector l
+ double thetaLD = M_PI;
+ double Delta = 0.5; //[GeV] This is the magnitude of the vector ... For starters <1 Gev
+ // given the first vector's (l) components and the inner product, find the second vector's (Delta) components
+ // sign of y component unknown, choose +
+ double vecDelta[2] = {Delta*l*cos(thetaLD), pow(Delta,2) - pow(Delta*l*cos(thetaLD),2)};
+
+
+ // quark flavors
+ const double elementary_charge = sqrt(4.*M_PI*1./137.036);
+
+ //TODO: check masses
+ quark_flavor u_quark;
+ u_quark.mass = 2e-3; //GeV
+ u_quark.e_charge = 2./3.*elementary_charge;
+
+ quark_flavor d_quark;
+ d_quark.mass = 2e-3;
+ d_quark.e_charge = -1./3.*elementary_charge;
+
+ quark_flavor c_quark;
+ c_quark.mass = 1.;
+ c_quark.e_charge = u_quark.e_charge;
+
+
+
+ std::complex<double> sum = 0.0;
+
+ //TODO: does pragma omp work?
+ //#pragma omp parallel for reduction(+:sum)
+ //error: user defined reduction not found for ‘sum’
+
for(int i =0; i < L; i += 1){
for(int j =0; j < L; j += 1){
int q1[2] = {i,j};
- quark = GetWilsonLine(WilsonLineCoordinate(i, j));
+ WilsonLine quark = GetWilsonLine(WilsonLineCoordinate(i, j));
double r1 = DistanceToOrigin(q1);
- for(int k =0; k < L; k +=1){
- for(int l =0; l < L; l +=1){
+ for(int v =0; v < L; v +=1){
+ for(int w =0; w < L; w +=1){
- int q2[2] = {k,l};
+ int q2[2] = {v,w};
double r2 = DistanceToOrigin(q2);
- double r = sqrt(pow((X(i)-X(k)),2)+pow(Y(j)-Y(l)),2);
-
- double WF = WFGammaqqbarLongitudinal(double r, double z, double Q, struct* quark_flavor);
- sum += Amplitude(0, quark, q2)*WF;
+ // phase factor vectors
+ double diffX[2] = {X(i) - X(v), Y(j) - Y(w)};
+ 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));
+
+ // wave functions
+ double WF_u = WFGammaqqbarLongitudinal(r, z0, Q, u_quark);
+ 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) * Amplitude(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
// Add phase factor
@@ -662,15 +709,15 @@ double IPGlasma::SumAmplitudes()
}
}
}
- cout << "# total sum = " << sum << endl;
- return sum;
+ cout << "# total amplitude squared = " << std::norm(sum) << endl;
+ return std::norm(sum);
}
// return 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, struct* quark_flavor)
+double IPGlasma::WFGammaqqbarLongitudinal(double r, double z, double Q, const quark_flavor &quark)
{
//Assuming r is the transverse size of the dipole
//constants: e_f (quark charge), e, N_C,
@@ -679,25 +726,10 @@ double IPGlasma::WFGammaqqbarLongitudinal(double r, double z, double Q, struct*
double lattice_spacing = xcoords[1]-xcoords[0];
// fine structure constant set to 1/137
- double e = sqrt(4*M_PI*1/137);
- double m_f = 0.16e-9; //GeV
- double e_f;
- if (f){ e_f = -1/3;}
- else { e_f = 2/3;}
-
- double z = 0.5; //GeV^0
- double Q = sqrt(10); //GeV
- double Delta = 0.5; //GeV
- double l = lattice_spacing/100; // << a^-1 (a = lattice const)
- double eps = sqrt(z*(1-z)*Q*Q*m_f*m_f);
+ double e = sqrt(4.*M_PI*1./137.036);
+
+ double epsilon = sqrt(z*(1-z)*Q*Q * pow(quark.mass, 2));
- return (e_f*e*sqrt(NC)*2*Q*Q*z*(1-z)*std::cyl_bessel_k(0, eps*r)*0.5/M_PI);
+ return (quark.e_charge * e * sqrt(NC) * 2*Q*Q*z*(1-z) * std::cyl_bessel_k(0, epsilon*r)*0.5/M_PI);
}
-// TODO: This has to be done inside the loop in SumAmplitudes
-double IPGlasma::DiffractiveCrossSection()
-{
-
- SumAmplitudes
- return 0;
-}
diff --git a/src/ipglasma.hpp b/src/ipglasma.hpp
index 81d30501e5078d03db9748b2af6e0f76e4901515..354240c0e505b8cebe230abf63fd1eb1619f66e2 100644
--- a/src/ipglasma.hpp
+++ b/src/ipglasma.hpp
@@ -60,16 +60,15 @@ public:
double Amplitude(double xpom, int q1[2], int q2[2]);
double Amplitude(double xpom, WilsonLine quark, int q2[2]);
double SumAmplitudes();
- double SumAmplitudes(int l);
- double WFGammaqqbarLongitudinal(double r, double z, double Q, struct* quark_flavor);
- struct q_flavor(double e_charge, double mass);
+ double SumAmplitudes(int w);
- //TODO: missä paras asettaa arvot?
- struct q_flavor{
- double e_charge;
+ struct quark_flavor{
double mass;
- }
-
+ double e_charge;
+ };
+
+ double WFGammaqqbarLongitudinal(double r, double z, double Q, const quark_flavor &quark);
+
private: