separated new stuff to their own class, needs testing

src/diffractivecrosssection.cpp

+/*
+ * Diffraction at sub-nucleon scale
+ * Dipole amplitude for a IPglasma nucleus
+ * Originally by Heikki Mäntysaari <mantysaari@bnl.gov>, 2015
+ * Edits and additions by Pyry Runko <pyry.j.runko@jyu.fi>, 2025
+ */
+
+
+#include "diffractiveCrossSection.hpp"
+#include <string>
+#include <gsl/gsl_randist.h>
+#include <fstream>
+#include <sstream>
+#include <CrossSectiontdlib>
+#include <cmath>
+#include <complex>
+#include <numeric>
+#include <omp.h>
+
+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);
+
+
+// compute complex amplitude with one wilson line and two indices
+std::complex<double> DiffractiveCrossSection::ComplexAmplitude(double xpom, WilsonLine quark, int q2[2])
+{
+    //TODO: make work with int?
+    //ApplyPeriodicBoundaryConditions(q1);
+    //ApplyPeriodicBoundaryConditions(q2);
+
+    // TODO: make work with indices (xMax?
+    // 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]
+    //    etc...
+    // ?TODO?: (double?) check if dipole is less than lattice spacing   
+
+    // First find corresponding grid indeces
+    WilsonLine antiquark = GetWilsonLine(WilsonLineCoordinate(q2[0], q2[1]));
+
+    WilsonLine prod;
+    try {
+        prod = quark.MultiplyByHermitianConjugate(antiquark);
+    } catch (...) {
+        cerr << "Matrix multiplication failed!" << endl;
+        cout << quark << endl;
+        cout << "Antiquark: " << q2[0] << ", " << q2[1] << endl;
+        cout << antiquark << endl;
+        exit(1);
+    }
+    std::complex<double> amp =  1.0 - 1.0/NC * prod.Trace();
+
+    double result = amp.real();
+    //if (result < 0) return 0;
+    //if (result > 1) return 1;
+    if (isnan(result))
+    {
+	cerr << "Wilson line trance NaN, antiquark coords " << q2[0] << ", " << q2[1] << "and quark wline" << quark << endl;
+        exit(1);
+    }
+
+    return amp;
+}
+
+
+
+//TODO: check
+double DiffractiveCrossSection::DistanceToOrigin(int q[2]){
+    // Get x and y coordinates corresponding to these indices
+    double x = X(q[0]);
+    double y = Y(q[1]);
+    // calculate length of the coordinate vector
+    return sqrt(x*x+y*y);
+
+}
+
+
+
+//TODO: check
+std::complex<double> DiffractiveCrossSection::SumAmplitudes()
+{
+    int L = xcoords.size();
+    cout << "# Summing with lattice size L = " << L << endl;
+    
+    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];
+    cout << "# Lattice spacing is " << lattice_spacing << " GeV^-1." << endl;
+
+
+    //TODO: loop over all values
+    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 +
+    // TODO: this is nonsense, 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); 
+
+    // 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;
+
+    quark_flavor d_quark;
+    d_quark.mass = 0.14; //GeV
+    d_quark.e_charge = -1./3.*elementary_charge;
+
+    quark_flavor c_quark;
+    c_quark.mass = 0.14; //GeV
+    c_quark.e_charge = u_quark.e_charge;
+
+
+
+    std::complex<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;
+
+    #pragma omp declare reduction(+:std::complex<double>:omp_out+=omp_in)
+    //TODO: check that this sums correctly
+    #pragma omp parallel for reduction(+:sum)
+    for(int i =0; i < L; i += 1){
+    
+        for(int j =0; j < L; j += 1){
+    
+            int q1[2] = {i,j};
+            const WilsonLine quark_wilsonline = GetWilsonLine(WilsonLineCoordinate(i, j));
+        //    double r1 = DistanceToOrigin(q1);
+    
+            for(int v =0; v < L; v +=1){
+                for(int w =0; w < L; w +=1){
+                    // check that the dipole size is not zero
+                    if (i == v && j == w) { continue;}
+
+
+                    int antiquark_indices[2] = {v,w};
+          //          double r2 = DistanceToOrigin(q2);
+
+                    // phase factor vectors
+                    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)};
+
+                    // dipole size for wave function
+                    double r = sqrt(pow(X(i)-X(v),2) + pow(Y(j)-Y(w),2));
+
+                    // wave functions
+                    double Psi_u = PsiGammaqqbarLongitudinal(r, z0, Q, u_quark);
+                    double Psi_d = 0;
+                    //TODO: finish mass checks
+                    if (u_quark.mass == d_quark.mass):{ Psi_d = Psi_u;}
+                    double Psi_d = PsiGammaqqbarLongitudinal(r, z0, Q, d_quark); 
+                    double Psi_c = PsiGammaqqbarLongitudinal(r, z0, Q, c_quark); 
+
+                    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
+                    //          Add z0 z1 integrations
+                    //          Add l Delta integrations
+                    //          Add delta functions
+                    //          Add prefactor and integration measure
+                    //          Compute modulus squared after averaging
+                }
+            }
+        }
+    }
+    
+    cout << "# total amplitude = " << endl;
+    cout << sum << endl;
+    return sum;
+}
+
+
+// normalized longitudinal photon polarization wave function
+// source arxiv.org/pdf/hep-ph/0606272
+// TODO: check factors of 2pi
+double DiffractiveCrossSection::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, 
+    //variables: r, Q, z, epsilon
+    //epsilon = sqrt(z(1-z)Q^2-m^2_f
+
+    // fine structure constant set to 1/137
+    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*z*(1-z) * std::cyl_bessel_k(0, epsilon*r)*0.5/M_PI);
+}
+

src/diffractivecrosssection.hpp

+/*
+ * Diffraction at sub-nucleon scale
+ * Dipole amplitude for a IPglasma nucleus
+ * Heikki Mäntysaari <mantysaari@bnl.gov>, 2015
+ */
+
+#ifndef diffractivecrosssection_hpp
+#define diffractivecrosssection_hpp
+
+//#include <string>
+//#include <vector>
+//#include "wilsonline.hpp"
+//#include "ipglasma.hpp"
+
+class DiffractiveCrossSection {
+public:
+    
+    double DistanceToOrigin(int q[2]);
+    std::complex<double> ComplexAmplitude(double xpom, int q1[2], int q2[2]);
+    std::complex<double> ComplexAmplitude(double xpom, WilsonLine quark, int q2[2]);
+    std::complex<double> SumAmplitudes();
+
+    struct quark_flavor{
+        double mass;
+        double e_charge;
+    };
+
+    double PsiGammaqqbarLongitudinal(double r, double z, double Q, const quark_flavor &quark);
+
+private:
+    
+    class IPGlasma;
+    //std::vector< double > xcoords;
+    //std::vector< double > ycoords;
+    //std::vector< WilsonLine  > wilsonlines;
+
+    //bool schwinger; // true if we use schwinger mechanism
+    //double schwinger_rc; // Use Schwinger for dipoles larger than schwinger_rc
+
+    //bool periodic_boundary_conditions;
+    
+    //std::string datafile;
+    
+    
+};
+
+const double FMGEV = 5.068;
+
+#endif /* diffractivecrosssection_hpp */