needs testing: integral for fixed l, Delta, z, Q, Mx

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
 )

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;
-}

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: