From 7fb6be529073e5057cd08c4c5de6f381a48c46ff Mon Sep 17 00:00:00 2001
From: "Pyry Runko (Puck)" <pyjorunk@student.jyu.fi>
Date: Wed, 5 Mar 2025 16:19:44 +0200
Subject: [PATCH] Amplitude calculation with fixed z, l and Delta
---
src/ipglasma.cpp | 78 ++++++++++++++++++++++--------------------------
src/ipglasma.hpp | 2 +-
src/main.cpp | 12 ++------
3 files changed, 39 insertions(+), 53 deletions(-)
diff --git a/src/ipglasma.cpp b/src/ipglasma.cpp
index 07e201b..faad217 100644
--- a/src/ipglasma.cpp
+++ b/src/ipglasma.cpp
@@ -550,19 +550,20 @@ std::complex<double> IPGlasma::SumAmplitudes()
//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
+ double length_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 vector_l[2] = {length_l, 0.}; //this is the vector l
+ double theta_lDelta = M_PI;
+ double length_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)};
+ // TODO: fix
+ double vector_Delta[2] = {length_Delta*length_l*cos(theta_lDelta), pow(length_Delta,2) - pow(length_Delta*length_l*cos(theta_lDelta),2)};
// quark flavors
const double elementary_charge = sqrt(4.*M_PI/137.036);
- //TODO: check masses
+ // masses from Kowalski Myotka Watt p. 13 "The BGBK analysis found..."
quark_flavor u_quark;
u_quark.mass = 0.14; //GeV
u_quark.e_charge = 2./3.*elementary_charge;
@@ -578,73 +579,63 @@ std::complex<double> IPGlasma::SumAmplitudes()
std::complex<double> sum = 0.0;
- //double sum = 0.0;
+ std::complex<double> imaginary_unit(0.,1.);
+ std::complex<double> complex_zero = (0.,0.);
+
+ int n_threads;
+ #pragma omp parallel
+ {
+ if (omp_get_thread_num()==0){
+ n_threads = omp_get_num_threads();
+ }
+ }
+ cout << "# Number of OMP threads is " << n_threads << endl;
- 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)
+ #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));
+ WilsonLine quark_wilsonline = GetWilsonLine(WilsonLineCoordinate(i, j));
// double r1 = DistanceToOrigin(q1);
for(int v =0; v < L; v +=1){
- for(int w =0; w < 100; w +=1){
+ for(int w =0; w < L; w +=1){
// check that the dipole size is not zero
if (i == v && j == w) { continue;}
- int q2[2] = {v,w};
+ int antiquark_indices[2] = {v,w};
// 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)};
- //cout << "# sumZX is " << sumZX << endl;
+ double difference_of_x_vectors[2] = {X(i) - X(v), Y(j) - Y(w)};
+ double sum_of_z_times_x[2] = {z0*X(i) + z1*X(v), z0*Y(j) + z1*Y(w)};
//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;
+ double Psi_u = PsiGammaqqbarLongitudinal(r, z0, Q, u_quark);
- double WF_d = WFGammaqqbarLongitudinal(r, z0, Q, d_quark);
- double WF_c = WFGammaqqbarLongitudinal(r, z0, Q, c_quark);
+ double Psi_d = PsiGammaqqbarLongitudinal(r, z0, Q, d_quark);
+ double Psi_c = PsiGammaqqbarLongitudinal(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 += ( Psi_u + Psi_d + Psi_c ) * ComplexAmplitude(0, quark_wilsonline, antiquark_indices)
+ * exp( - imaginary_unit * ( std::inner_product(std::begin(vector_l), std::end(vector_l), std::begin(difference_of_x_vectors), complex_zero)
+ + std::inner_product(std::begin(vector_Delta), std::end(vector_Delta), std::begin(sum_of_z_times_x), complex_zero) ) );
//TODO: Finish eq (12) from notes
// Finish amplitude
- // Add phase factor
// Add z0 z1 integrations
// Add l Delta integrations
// Add delta functions
// Add prefactor and integration measure
- //
+ // Compute modulus squared after averaging
}
}
}
@@ -656,9 +647,10 @@ std::complex<double> IPGlasma::SumAmplitudes()
}
-// longitudinal photon polarization wave function
+// normalized longitudinal photon polarization wave function
// source arxiv.org/pdf/hep-ph/0606272
-double IPGlasma::WFGammaqqbarLongitudinal(double r, double z, double Q, const quark_flavor &quark)
+// TODO: check factors of 2pi
+double IPGlasma::PsiGammaqqbarLongitudinal(double r, double z, double Q, const quark_flavor &quark)
{
//r is the transverse size of the dipole
//constants: e_f (quark charge), e, N_C,
@@ -669,6 +661,6 @@ double IPGlasma::WFGammaqqbarLongitudinal(double r, double z, double Q, const qu
double e = sqrt(4.*M_PI/137.036);
double epsilon = sqrt(z*(1-z)*Q*Q + pow(quark.mass, 2));
- return (quark.e_charge * e * sqrt(NC) * 2*Q*Q*z*(1-z) * std::cyl_bessel_k(0, epsilon*r)*0.5/M_PI);
+ return (quark.e_charge * e * sqrt(NC) * 2*Q*z*(1-z) * std::cyl_bessel_k(0, epsilon*r)*0.5/M_PI);
}
diff --git a/src/ipglasma.hpp b/src/ipglasma.hpp
index 049caad..a649095 100644
--- a/src/ipglasma.hpp
+++ b/src/ipglasma.hpp
@@ -66,7 +66,7 @@ public:
double e_charge;
};
- double WFGammaqqbarLongitudinal(double r, double z, double Q, const quark_flavor &quark);
+ double PsiGammaqqbarLongitudinal(double r, double z, double Q, const quark_flavor &quark);
private:
diff --git a/src/main.cpp b/src/main.cpp
index 9c7da0e..190455d 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -60,20 +60,14 @@ int main(int argc, char* argv[])
{
//TODO if num_threads > 1 cout << "using omp";
- if (omp_get_max_threads() > 1){
- cout << "Running with OpenMP." << " Using " << omp_get_max_threads() << " threads." << std::endl;
- }
- else { cout << "Running without OpenMP." << endl;}
-
-
- if (argc < 3 || argc > 4){
- cerr << "Usage: " << argv[0] << " filename 1(integrate)/0(print amplitudes) [task_id]" << endl;
+ if (argc != 2){
+ cerr << "# Usage: " << argv[0] << " filename 0(integrate)/1(print dipole amplitudes)" << endl;
}
string fname = argv[1];
IPGlasma ipglasma(fname, 0, BINARY);
- if (!argv[2]){
+ if (argv[2]){
printDipoleAmplitudes(fname, ipglasma);
}
--
GitLab