From acae3254eec8b7ed4662bc82f0dbfa9f6e77681b Mon Sep 17 00:00:00 2001
From: Pyry Runko <pyjorunk@hytkl2011-38.ws.jyu.fi>
Date: Mon, 17 Mar 2025 16:08:11 +0200
Subject: [PATCH] Using wilsonline reader as library
---
CMakeLists.txt | 6 +
src/CMakeLists.txt | 6 +-
src/diffractivecrosssection.cpp | 65 ++--
src/diffractivecrosssection.hpp | 31 +-
src/ipglasma.cpp | 666 --------------------------------
src/ipglasma.hpp | 90 -----
src/main.cpp | 13 +-
src/wilsonline.cpp | 403 -------------------
src/wilsonline.hpp | 79 ----
9 files changed, 72 insertions(+), 1287 deletions(-)
delete mode 100644 src/ipglasma.cpp
delete mode 100644 src/ipglasma.hpp
delete mode 100644 src/wilsonline.cpp
delete mode 100644 src/wilsonline.hpp
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4056a44..33617fe 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -16,6 +16,12 @@ include_directories(${GSL_INCLUDE_DIRS})
find_package(OpenMP REQUIRED)
+include_directories("../wilsonline_reader/src/")
+link_directories("../wilsonline_reader/build/src")
+
+
+
+
add_subdirectory(src)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}")
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 82452ea..e7fbb2b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -2,8 +2,9 @@
add_executable(
wline
main.cpp
- wilsonline.cpp
- ipglasma.cpp
+ #wilsonline.cpp
+ # ipglasma.cpp
+ diffractivecrosssection.cpp
)
if(OpenMP_CXX_FOUND)
@@ -14,6 +15,7 @@ endif()
target_link_libraries(
wline
PRIVATE
+ wilsonlinereader
GSL::gsl
GSL::gslcblas
OpenMP::OpenMP_CXX
diff --git a/src/diffractivecrosssection.cpp b/src/diffractivecrosssection.cpp
index 161f475..adc4978 100644
--- a/src/diffractivecrosssection.cpp
+++ b/src/diffractivecrosssection.cpp
@@ -6,13 +6,16 @@
*/
-#include "diffractiveCrossSection.hpp"
+#include "diffractivecrosssection.hpp"
+#include "ipglasma.hpp"
+#include "wilsonline.hpp"
#include <string>
#include <gsl/gsl_randist.h>
#include <fstream>
#include <sstream>
-#include <CrossSectiontdlib>
+#include <cstdlib>
#include <cmath>
+#include <string>
#include <complex>
#include <numeric>
#include <omp.h>
@@ -27,6 +30,15 @@ const int NC=3;
//int FindIndex(double val, std::vector<double> &vec);
+DiffractiveCrossSection::DiffractiveCrossSection(std::string file)
+{
+ ipglasma = new IPGlasma(file);
+}
+
+DiffractiveCrossSection::~DiffractiveCrossSection()
+{
+ delete ipglasma;
+}
// compute complex amplitude with one wilson line and two indices
std::complex<double> DiffractiveCrossSection::ComplexAmplitude(double xpom, WilsonLine quark, int q2[2])
@@ -44,7 +56,7 @@ std::complex<double> DiffractiveCrossSection::ComplexAmplitude(double xpom, Wils
// ?TODO?: (double?) check if dipole is less than lattice spacing
// First find corresponding grid indeces
- WilsonLine antiquark = GetWilsonLine(WilsonLineCoordinate(q2[0], q2[1]));
+ WilsonLine antiquark = ipglasma->GetWilsonLine(ipglasma->WilsonLineCoordinate(q2[0], q2[1]));
WilsonLine prod;
try {
@@ -61,7 +73,7 @@ std::complex<double> DiffractiveCrossSection::ComplexAmplitude(double xpom, Wils
double result = amp.real();
//if (result < 0) return 0;
//if (result > 1) return 1;
- if (isnan(result))
+ if (std::isnan(result))
{
cerr << "Wilson line trance NaN, antiquark coords " << q2[0] << ", " << q2[1] << "and quark wline" << quark << endl;
exit(1);
@@ -75,8 +87,8 @@ std::complex<double> DiffractiveCrossSection::ComplexAmplitude(double xpom, Wils
//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]);
+ double x = ipglasma->X(q[0]);
+ double y = ipglasma->Y(q[1]);
// calculate length of the coordinate vector
return sqrt(x*x+y*y);
@@ -87,8 +99,8 @@ double DiffractiveCrossSection::DistanceToOrigin(int q[2]){
//TODO: check
std::complex<double> DiffractiveCrossSection::SumAmplitudes()
{
- int L = xcoords.size();
- cout << "# Summing with lattice size L = " << L << endl;
+ int N = ipglasma->GetN(); //xcoords.size();
+ cout << "# Summing with lattice size N = " << N << endl;
double z0 = 0.5; //Momentum fractions
double z1 = 1-z0;
@@ -96,7 +108,7 @@ std::complex<double> DiffractiveCrossSection::SumAmplitudes()
double Mx = 1; //[GeV] Invarint mass of system X
// in GeV, see LoadData
- double lattice_spacing = xcoords[1]-xcoords[0];
+ double lattice_spacing = ipglasma->GetLatticeSpacing(); //xcoords[1]-xcoords[0];
cout << "# Lattice spacing is " << lattice_spacing << " GeV^-1." << endl;
@@ -123,9 +135,9 @@ std::complex<double> DiffractiveCrossSection::SumAmplitudes()
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;
+ quark_flavor s_quark;
+ s_quark.mass = 0.14; //GeV
+ s_quark.e_charge = u_quark.e_charge;
@@ -145,16 +157,16 @@ std::complex<double> DiffractiveCrossSection::SumAmplitudes()
#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 i =0; i < N; i += 1){
- for(int j =0; j < L; j += 1){
+ for(int j =0; j < N; j += 1){
int q1[2] = {i,j};
- const WilsonLine quark_wilsonline = GetWilsonLine(WilsonLineCoordinate(i, j));
+ const WilsonLine quark_wilsonline = ipglasma->GetWilsonLine(ipglasma->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 v =0; v < N; v +=1){
+ for(int w =0; w < N; w +=1){
// check that the dipole size is not zero
if (i == v && j == w) { continue;}
@@ -163,21 +175,22 @@ std::complex<double> DiffractiveCrossSection::SumAmplitudes()
// 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)};
+ double difference_of_x_vectors[2] = {ipglasma->X(i) - ipglasma->X(v), ipglasma->Y(j) - ipglasma->Y(w)};
+ double sum_of_z_times_x[2] = {z0*ipglasma->X(i) + z1*ipglasma->X(v), z0*ipglasma->Y(j) + z1*ipglasma->Y(w)};
// dipole size for wave function
- double r = sqrt(pow(X(i)-X(v),2) + pow(Y(j)-Y(w),2));
+ double r = sqrt(pow(ipglasma->X(i)-ipglasma->X(v),2) + pow(ipglasma->Y(j)-ipglasma->Y(w),2));
// wave functions
+ // TODO: make quark loop and sum wavefunctions together, check if mass same..., remember e-charge for iq = 0; iq < quarks.size(); iq++
+ //
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);
+ //double Psi_d = 0;
+ //if (u_quark.mass == d_quark.mass){ Psi_d = Psi_u;}
+ //else: Psi_d = PsiGammaqqbarLongitudinal(r, z0, Q, d_quark);
+ //double Psi_s = PsiGammaqqbarLongitudinal(r, z0, Q, s_quark);
- sum += ( Psi_u + Psi_d + Psi_c ) * ComplexAmplitude(0, quark_wilsonline, antiquark_indices)
+ sum += Psi_u * 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) ) );
diff --git a/src/diffractivecrosssection.hpp b/src/diffractivecrosssection.hpp
index 739ee1d..b3026b9 100644
--- a/src/diffractivecrosssection.hpp
+++ b/src/diffractivecrosssection.hpp
@@ -7,14 +7,16 @@
#ifndef diffractivecrosssection_hpp
#define diffractivecrosssection_hpp
-//#include <string>
-//#include <vector>
-//#include "wilsonline.hpp"
-//#include "ipglasma.hpp"
+#include <string>
+#include <vector>
+#include <ipglasma.hpp>
+#include <wilsonline.hpp>
class DiffractiveCrossSection {
public:
-
+ DiffractiveCrossSection(std::string file);
+ ~DiffractiveCrossSection();
+
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]);
@@ -27,23 +29,20 @@ public:
double PsiGammaqqbarLongitudinal(double r, double z, double Q, const quark_flavor &quark);
+ IPGlasma* GetIPGlasma() { return ipglasma; }
+
private:
- class IPGlasma;
+ IPGlasma *ipglasma;
+ //WilsonLine& IPGlasma::GetWilsonLine(double x, double y);
+ //WilsonLine& IPGlasma::GetWilsonLine(int i);
+ //int IPGlasma::WilsonLineCoordinate(int indx, int indy);
+ //double IPGlasma::X(int ix);
+ //double IPGlasma::Y(int iy);
//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 */
diff --git a/src/ipglasma.cpp b/src/ipglasma.cpp
deleted file mode 100644
index faad217..0000000
--- a/src/ipglasma.cpp
+++ /dev/null
@@ -1,666 +0,0 @@
-/*
- * Diffraction at sub-nucleon scale
- * Dipole amplitude for a IPglasma nucleus
- * Heikki Mäntysaari <mantysaari@bnl.gov>, 2015
- */
-
-
-#include "ipglasma.hpp"
-#include <string>
-#include <gsl/gsl_randist.h>
-#include <fstream>
-#include <sstream>
-#include <cstdlib>
-#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);
-
-
-/*
- * Calculate dipole amplitude form Wilson lines
- * Search the closest grid point that corresponds to the given quark/antiquark coordinate
- * Then calcualte 1 - 1/Nc Tr U(quark) U^dagger(antiquark)
- */
-double IPGlasma::Amplitude(double xpom, double q1[2], double q2[2] )
-{
- ApplyPeriodicBoundaryConditions(q1);
- ApplyPeriodicBoundaryConditions(q2);
-
- // 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]
- or q2[0] < xcoords[0] or q2[0] > xcoords[xcoords.size()-1]
- or q2[1] < ycoords[0] or q2[1] > ycoords[ycoords.size()-1])
- {
- if (periodic_boundary_conditions)
- cerr << "WTF, I'm here..." << endl;
-
- return 0;
-
- }
-
- double r = sqrt( pow(q1[0]-q2[0],2) + pow(q1[1]-q2[1],2));
-
- // Smaller dipole than the grid
- if (r < std::abs( xcoords[1] - xcoords[0]))
- return 0;
-
-
- // First find corresponding grid indeces
- WilsonLine quark = GetWilsonLine(q1[0], q1[1]);
- WilsonLine antiquark = GetWilsonLine(q2[0], q2[1]);
-
- //antiquark = antiquark.HermitianConjugate();
-
- WilsonLine prod;
- try {
- prod = quark.MultiplyByHermitianConjugate(antiquark);
- //prod = quark*antiquark;
- } catch (...) {
- cerr << "Matrix multiplication failed!" << endl;
- cout << "Quark: " << q1[0] << ", " << q1[1] << 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;
- return result;
- if (result > 1)
- return 1;
- if (result < 0)
- return 0;
-
- if (isnan(result))
- {
- cerr << "Wilson line trance NaN, quark coords " << q1[0] << ", " << q1[1] << " and " << q2[0] << ", " << q2[1] << endl;
- exit(1);
- }
-
- return result;
-}
-// Stupid copypaste
-double IPGlasma::AmplitudeImaginaryPart(double xpom, double q1[2], double q2[2] )
-{
-
- ApplyPeriodicBoundaryConditions(q1);
- ApplyPeriodicBoundaryConditions(q2);
-
- // Out of grid? Return 1 (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]
- or q2[0] < xcoords[0] or q2[0] > xcoords[xcoords.size()-1]
- or q2[1] < ycoords[0] or q2[1] > ycoords[ycoords.size()-1])
- return 0;
-double r = sqrt( pow(q1[0]-q2[0],2) + pow(q1[1]-q2[1],2));
- if (r < std::abs( xcoords[1] - xcoords[0]))
- return 0;
-
- // First find corresponding grid indeces
- WilsonLine quark = GetWilsonLine(q1[0], q1[1]);
- WilsonLine antiquark = GetWilsonLine(q2[0], q2[1]);
-// antiquark = antiquark.HermitianConjugate();
-
- WilsonLine prod;
- try {
- //prod = quark*antiquark;
- prod = quark.MultiplyByHermitianConjugate(antiquark);
-
- } catch (...) {
- cerr << "Matrix multiplication failed!" << endl;
- cout << "Quark: " << q1[0] << ", " << q1[1] << 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.imag();
- return result;
- if (result > 1)
- return 1;
- if (result < 0)
- return 0;
-
- return result;
-}
-
-WilsonLine& IPGlasma::GetWilsonLine(double x, double y)
-{
- double q[2]={x,y};
- ApplyPeriodicBoundaryConditions(q);
- x=q[0];
- y=q[1];
- std::vector<int> coords = LatticeCoordinates(x,y);
-
- // Handle edges
- if (coords[0] < 0)
- coords[0] = 0;
- if (coords[0] >= xcoords.size())
- coords[0] = xcoords.size()-1;
-
- if (coords[1] < 0)
- coords[1] = 0;
- if (coords[1] >= ycoords.size())
- coords[1] = ycoords.size()-1;
-
- return GetWilsonLine( WilsonLineCoordinate(coords[0],coords[1]));
-
-}
-
-IPGlasma::IPGlasma(std::string file)
-{
- int load = LoadData(file, 5.12/512.0); // By default assume lattice spacing 0.01fm, or 5.12fm 512^2 lattice
-
- if (load<0)
- exit(1);
-
- periodic_boundary_conditions=false;
-}
-
-IPGlasma::IPGlasma(std::string file, double step, WilsonLineDataFileType type)
-{
- int load =LoadData(file, step, type);
-
- if (load<0)
- exit(1);
-
- periodic_boundary_conditions=false;
-}
-
-
-// Note that here step is in fm!
-int IPGlasma::LoadData(std::string fname, double step, WilsonLineDataFileType type)
-{
- // Load data
- datafile=fname;
-
- if (type == BINARY)
- return LoadBinaryData(fname, step);
-
- // Syntax: x y [fm/indeces] matrix elements Re Im for elements (0,0), (0,1), (0,2), (1,0), ...
- std::ifstream f(fname.c_str());
-
- if (!f.is_open())
- {
- std::cerr << "Could not open file " << fname << " IPGlasma::LoadData " << std::endl;;
- return -1;
- }
- std::string line;
-
- while(!f.eof() )
- {
- std::getline(f, line);
- if (line[0]=='#' or line.length() < 10) // Comment line, or empty
- {
- continue;
- }
-
- // Parse
- std::stringstream ss(line);
- // Get coordinates
- double x,y;
- ss >> x;
- ss >> y;
-
- // Datafile is in fm, but we want to use GeVs in this code
- // Once we have loaded all points, we will sift all coordinates such that 0 is at the center
-
-
- x = step*x*FMGEV;
- y = step*y*FMGEV;
-
-
- // Read rows and columns
- std::vector< std::vector< std::complex<double> > > matrix;
- for (int row=0; row < 3; row++)
- {
- std::vector< std::complex<double> > tmprow;
- for (int col=0; col<3; col++)
- {
- double real, imag;
- ss >> real;
- ss >> imag;
- std::complex<double> element (real, imag);
- tmprow.push_back(element);
- }
- matrix.push_back(tmprow);
- }
-
- // Save Wilson line
- WilsonLine w(matrix);
- w = w.Transpose();
- wilsonlines.push_back(w);
-
- if (w.Size() != 3)
- {
- cerr << "Matrix at coordinate " << x << " " << y << " has size " << w.Size() << endl;
- exit(1);
- }
-
- // We assume that grid is symmetric, so don't save same valeus multiple times
- // In the datafile the coordinates are increasing
- if (xcoords.size() == 0 or x > xcoords[xcoords.size()-1])
- xcoords.push_back(x);
- if (ycoords.size() == 0 or y > ycoords[ycoords.size()-1])
- ycoords.push_back(y);
-
- }
-
- // Shift coordinates s.t. 0 fm is at the center
- double center = xcoords[xcoords.size()/2];
- for (int i=0; i<xcoords.size(); i++)
- {
- xcoords[i] -= center;
- ycoords[i] -= center;
- }
- f.close();
-
- if (xcoords.size() != ycoords.size())
- {
- cerr << "xcoords.size() != ycoords.size(), probably uncomplete input data? Datafile " << fname << " - IPGlasma::LoadData" << endl;
- return -1;
- }
-
- if (xcoords.size() < 10 or ycoords.size()<10)
- {
- cerr << "Grid size is " << xcoords.size() << " x " << ycoords.size() << ", this makes no sense! File " << fname << " - IPGlasma::LoadData" << endl;
- return -1;
- }
-
-
-
- // Now, given that we have a point (x,y), we can find the index xind such that
- // xcoords[xind] is closest to x, and similarly ind
- // Then, the corresponding Wilson line is
- // wilsonlines[ xcoords.size()*xind + yind]
- // Of course this is symmetric and we could just as well swap xind and yind
-
- SetSchwinger(false);
- // std::cout <<"# Loaded " << wilsonlines.size() << " Wilson lines from file " << datafile << ", grid size " << xcoords.size() << " x " << ycoords.size() << " grid range [" << xcoords[0] << ", " << xcoords[xcoords.size()-1] << "]" << " step size " << xcoords[1]-xcoords[0] << " GeV^-1" << std::endl;
-
-
- return 0;
-}
-
-
-/*
- * Load data from binary file
- */
-int IPGlasma::LoadBinaryData(std::string fname, double step)
-{
- std::ifstream InStream;
- InStream.precision(10);
- InStream.open(fname.c_str(), std::ios::in | std::ios::binary);
- int N;
- int Nc;
- double L,a;
-
- if(InStream.is_open())
- {
- // READING IN PARAMETERS -----------------------------------------------------------------//
- double temp;
- InStream.read(reinterpret_cast<char*>(&N), sizeof(int));
- InStream.read(reinterpret_cast<char*>(&Nc), sizeof(int));
- InStream.read(reinterpret_cast<char*>(&L), sizeof(double));
- InStream.read(reinterpret_cast<char*>(&a), sizeof(double));
- InStream.read(reinterpret_cast<char*>(&temp), sizeof(double));
-
- std::cout << "# Reading Wilson line: binary format, size is " << N << ", Nc " << Nc << ", length is [fm] " << L << ", a is [fm]" << a << std::endl;
-
- // Init Wilson lines - fill in later
- wilsonlines.clear();
- for (unsigned int i=0; i<N*N; i++)
- {
- WilsonLine w;
- wilsonlines.push_back(w);
- }
-
-
- // READING ACTUAL DATA --------------------------------------------------------------------//
- double ValueBuffer;
- int INPUT_CTR=0;
- double re,im;
-
- while( InStream.read(reinterpret_cast<char*>(&ValueBuffer), sizeof(double)))
- {
- if(INPUT_CTR%2==0) //this is the real part
- {
- re=ValueBuffer;
- }
- else // this is the imaginary part, write then to variable //
- {
- im=ValueBuffer;
-
- int TEMPINDX=((INPUT_CTR-1)/2);
- int PositionIndx = TEMPINDX / 9;
- int ix = PositionIndx / N;
- int iy = PositionIndx - N*ix;
-
-
- int MatrixIndx=TEMPINDX - PositionIndx*9;
- int j=MatrixIndx/3;
- int k=MatrixIndx-j*3;
-
- int indx = N*iy + ix;
- wilsonlines[indx].Set(j,k, std::complex<double> (re,im));
- }
- INPUT_CTR++;
- }
- }
- else
- {
- std::cerr << "ERROR COULD NOT OPEN FILE " << fname << std::endl;
- }
-
- InStream.close();
-
- // Fill x and y coordinate values in GeV^-1
- for (unsigned int i=0; i<N; i++)
- {
- xcoords.push_back(i*a*5.068);
- ycoords.push_back(i*a*5.068);
- }
-
- // Shift coordinates s.t. 0 fm is at the center
- double center = xcoords[xcoords.size()/2];
- for (int i=0; i<xcoords.size(); i++)
- {
- xcoords[i] -= center;
- ycoords[i] -= center;
- }
-
- SetSchwinger(false);
- return 0;
-}
-
-double IPGlasma::MinX()
-{
- return xcoords[0];
-}
-
-double IPGlasma::MaxX()
-{
- return xcoords[xcoords.size()-1];
-}
-
-
-double IPGlasma::XStep()
-{
- return xcoords[1] - xcoords[0];
-}
-
-
-void IPGlasma::ApplyPeriodicBoundaryConditions(double q[2])
-{
-
- if (periodic_boundary_conditions == false) return;
-
- double Lx = xcoords[xcoords.size()-1]-xcoords[0];
- double Ly = ycoords[ycoords.size()-1]-ycoords[0];
- while (q[0] < xcoords[0])
- q[0] = q[0] + Lx;
-
- while (q[0] > xcoords[xcoords.size()-1])
- q[0] = q[0] - Lx;
-
- while (q[1] < ycoords[0])
- q[1] = q[1] + Ly;
-
- while (q[1] > ycoords[ycoords.size()-1])
- q[1] = q[1] - Ly;
-
-
-}
-
-std::string IPGlasma::InfoStr()
-{
- std::stringstream ss;
- ss << "# IPGlasma loaded from file " << datafile << " lattice " << xcoords.size() << "^2 range [" << xcoords[0]/5.068 << ", " << xcoords[xcoords.size()-1]/5.068 << "] fm" << endl ;
- if (periodic_boundary_conditions) ss << "# Periodic boundary conditions" << endl;
- if (schwinger) ss << "# schwinger mechanism included, rc=" << schwinger_rc << " GeV^-1" << endl;
- return ss.str();
-}
-
-std::vector<double> &IPGlasma::GetXCoordinates()
-{
- return xcoords;
-}
-
-void IPGlasma::SetSchwinger(bool s, double rc)
-{
- schwinger = s;
- schwinger_rc = rc;
-
- if (s)
- {
- cerr << "NOTE: Schwinger Mechanism is not implemented in IPGlasma!" << endl;
- }
-}
-
-
-std::vector<int> IPGlasma::LatticeCoordinates(double x, double y)
-{
- std::vector<int> ret;
-
- // Note: My lattice is from -L/2 to L/2, so I need to shift the coordinates
- x = x + xcoords[xcoords.size()-1];
- y = y + ycoords[ycoords.size()-1];
- double lattice_spacing = xcoords[1]-xcoords[0];
-
- int ix = x/lattice_spacing;
- int iy = y/lattice_spacing;
-
- return std::vector<int> {ix, iy};
-}
-
-int IPGlasma::WilsonLineCoordinate(int xind, int yind)
-{
- return xcoords.size()*xind + yind;
-}
-
-
-
-// compute complex amplitude with one wilson line and two indices
-std::complex<double> IPGlasma::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 IPGlasma::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> IPGlasma::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: 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};
- 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)};
-
- //TODO: is this correct?
- // 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 = 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
- // Finish amplitude
- // 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 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,
- //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);
-}
-
diff --git a/src/ipglasma.hpp b/src/ipglasma.hpp
deleted file mode 100644
index a649095..0000000
--- a/src/ipglasma.hpp
+++ /dev/null
@@ -1,90 +0,0 @@
-/*
- * Diffraction at sub-nucleon scale
- * Dipole amplitude for a IPglasma nucleus
- * Heikki Mäntysaari <mantysaari@bnl.gov>, 2015
- */
-
-#ifndef ipglasma_hpp
-#define ipglasma_hpp
-
-#include <string>
-#include <vector>
-#include "wilsonline.hpp"
-
-enum WilsonLineDataFileType
-{
- BINARY,
- TEXT
-};
-
-class IPGlasma {
-public:
-
- IPGlasma(std::string fname);
- // Specify the step size, note that here step is in fm!
- IPGlasma(std::string fname, double step, WilsonLineDataFileType=TEXT);
-
- int LoadData(std::string fname, double step, WilsonLineDataFileType type=TEXT);
- int LoadBinaryData(std::string fname, double step);
-
- // Evaluate dipole ampltitude, qaurks at coordinates x1 and x2
- // Array points are x and y coordinates
- double Amplitude(double xpom, double q1[2], double q2[2]);
- double AmplitudeImaginaryPart(double xpom, double q1[2], double q2[2] );
-
- WilsonLine& GetWilsonLine( double x, double y); // Find Wilson line at the given point (x,y) [GeV^-1]
-
- std::string InfoStr();
-
- double MinX();
- double MaxX();
-
- double XStep(); // Grid spacing in x
-
- void SetPeriodicBoundaryConditions(bool s) { periodic_boundary_conditions = s; }
-
- std::vector<double> &GetXCoordinates();
-
- void SetSchwinger(bool s, double rc=0);
- void ApplyPeriodicBoundaryConditions(double q[2]);
-
- std::vector<int> LatticeCoordinates(double x, double y);
- int WilsonLineCoordinate(int xind, int yind);
- WilsonLine& GetWilsonLine(int i) { return wilsonlines[i]; }
-
- double X(int ix) { return xcoords[ix]; }
- double Y(int iy) { return ycoords[iy]; }
-
- //Pyry added:
- 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:
-
-
- 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 /* ipglasma_hpp */
diff --git a/src/main.cpp b/src/main.cpp
index 190455d..0e2f033 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -7,9 +7,11 @@
#include <iomanip>
#include <omp.h>
+#include "diffractivecrosssection.hpp"
+
using namespace std;
-int printDipoleAmplitudes(string fname, IPGlasma ipglasma)
+int printDipoleAmplitudes(string fname, DiffractiveCrossSection & diffxs)
{
//ofstream output0;
//Initialize params
@@ -43,7 +45,7 @@ int printDipoleAmplitudes(string fname, IPGlasma ipglasma)
q2[0] = -q1[0];
q2[1] = -q1[1];
- dipoleAmp += ipglasma.Amplitude(0, q1, q2);
+ dipoleAmp += diffxs.GetIPGlasma()->Amplitude(0, q1, q2);
}
dipoleAmp = dipoleAmp/N;
@@ -65,15 +67,16 @@ int main(int argc, char* argv[])
}
string fname = argv[1];
- IPGlasma ipglasma(fname, 0, BINARY);
+ //IPGlasma ipglasma(fname, 0, BINARY);
+ DiffractiveCrossSection cross_section(argv[1]);
if (argv[2]){
- printDipoleAmplitudes(fname, ipglasma);
+ printDipoleAmplitudes(fname, cross_section);
}
else{
- std::complex<double> sum = ipglasma.SumAmplitudes();
+ std::complex<double> sum = cross_section.SumAmplitudes();
// Save results
//std::filesystem::create_directory("output_"+fname);
ofstream outfile("da_sum_ +" + fname + ".txt");
diff --git a/src/wilsonline.cpp b/src/wilsonline.cpp
deleted file mode 100644
index f07f9ef..0000000
--- a/src/wilsonline.cpp
+++ /dev/null
@@ -1,403 +0,0 @@
-/*
- * Diffraction at sub-nucleon scale
- * Wilson line / SU(3) matrix handling
- * Heikki Mäntysaari <mantysaari@bnl.gov>, 2015
- * Standalone class, no external dependences, please!
- */
-
-#include "wilsonline.hpp"
-#include <vector>
-#include <complex>
-#include <iostream>
-#include <gsl/gsl_blas.h>
-#include <gsl/gsl_complex.h>
-#include <gsl/gsl_complex_math.h>
-
-using std::cerr;
-using std::cout;
-using std::endl;
-
-
-WilsonLine::WilsonLine() // Initialize everything to 0
-{
- wline_gsl_matrix = gsl_matrix_complex_alloc(size,size);
-
- // By default 3x3 zero matrix
- gsl_matrix_complex_set_zero(wline_gsl_matrix);
-
-
-}
-WilsonLine::WilsonLine(std::vector< std::vector< std::complex<double> > > &d)
-{
- wline_gsl_matrix = gsl_matrix_complex_alloc(size,size);
-
- for (int i=0; i<size; i++)
- {
- for (int j=0; j<size; j++)
- {
- gsl_matrix_complex_set(wline_gsl_matrix,i,j,gsl_complex_rect(d[i][j].real(),d[i][j].imag()));
- }
- }
-}
-
-WilsonLine::WilsonLine(gsl_matrix_complex* m)
-{
- wline_gsl_matrix = m;
-}
-
-WilsonLine::WilsonLine(const WilsonLine &m)
-{
- // Copy constructor: allocate memory for this, and copy data from m
- wline_gsl_matrix = gsl_matrix_complex_alloc(size,size);
- for (int i=0; i<size; i++)
- {
- for (int j=0; j<size; j++)
- {
- std::complex<double> tmp = m.Element(i,j);
- gsl_matrix_complex_set(wline_gsl_matrix, i,j,gsl_complex_rect(tmp.real(), tmp.imag()));
- }
- }
-}
-
-WilsonLine& WilsonLine::operator=(const WilsonLine &m)
-{
- for (int i=0; i<size; i++)
- {
- for (int j=0; j<size; j++)
- {
- std::complex<double> tmp = m.Element(i,j);
- gsl_matrix_complex_set(wline_gsl_matrix, i,j,gsl_complex_rect(tmp.real(), tmp.imag()));
- }
- }
- return *this;
-}
-
-WilsonLine::~WilsonLine()
-{
- gsl_matrix_complex_free(wline_gsl_matrix);
-}
-
-void WilsonLine::Set(int row, int column, std::complex<double> value)
-{
- if (row >= size)
- {
- cerr << "Invalid row index " << row << " num of rows in matrix " << size << endl;
- return;
- }
- if (column >= size )
- {
- cerr << "Invalid column index " << column << " num of rows in matrix " << size << endl;
- return;
- }
-
- gsl_matrix_complex_set(wline_gsl_matrix,row,column,
- gsl_complex_rect(value.real(),value.imag()));
-
-
-}
-
-WilsonLine WilsonLine::operator*(WilsonLine& w)
-{
- gsl_matrix_complex *result = gsl_matrix_complex_alloc(size,size);
- gsl_matrix_complex_set_zero(result);
-
-
- gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, gsl_complex_rect(1,0), wline_gsl_matrix, w.GetGslMatrix(), gsl_complex_rect(0,0), result);
-
- return WilsonLine(result);
-}
-
-// Multiplies this by w^\dagger, returns the product
-// Fast in BLAS
-WilsonLine WilsonLine::MultiplyByHermitianConjugate(WilsonLine& w)
-{
- gsl_matrix_complex *result = gsl_matrix_complex_alloc(size,size);
-
- gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, gsl_complex_rect(1,0), wline_gsl_matrix, w.GetGslMatrix(), gsl_complex_rect(0,0), result);
-
- return result;
-
-}
-
-WilsonLine WilsonLine::operator*(std::complex<double> t)
-{
- std::cerr << "Multiplication by a number not implemented" << endl;
- exit(1);
-}
-
-
-WilsonLine WilsonLine::operator+(WilsonLine& w)
-{
-
- std::cerr << "WilsonLine::operator+ not yet implemented for GSL matrix" << endl;
- exit(1);
-
-
-
-}
-
-
-
-
-WilsonLine WilsonLine::ComplexConjugate()
-{
-
- std::cerr << "WilsonLine:ComplexConjugate not yet implemented for GSL matrix" << endl;
- exit(1);
-
-}
-
-WilsonLine WilsonLine::Transpose()
-{
-
- std::cerr << "WilsonLine:Transpose not yet implemented for GSL matrix" << endl;
- exit(1);
-
-}
-
-WilsonLine WilsonLine::HermitianConjugate()
-{
-
- WilsonLine res = ComplexConjugate().Transpose();
- return res;
-}
-
-std::complex<double> WilsonLine::Trace()
-{
- gsl_complex sum;
- GSL_SET_COMPLEX(&sum, 0,0);
-
- for (unsigned int i=0; i< size; i++)
- {
- gsl_complex c = gsl_matrix_complex_get(wline_gsl_matrix, i, i);
- sum = gsl_complex_add(sum,c);
- }
-
- return std::complex<double>(GSL_REAL(sum), GSL_IMAG(sum));
-}
-
-int WilsonLine::Size()
-{
- return size;
-}
-
-std::complex<double> WilsonLine::Element(int row , int col) const
-{
- gsl_complex c = gsl_matrix_complex_get(wline_gsl_matrix, row, col);
- std::complex<double> cc(GSL_REAL(c), GSL_IMAG(c));
- return cc;
-
-}
-std::ostream& operator<<(std::ostream& os, WilsonLine& wl)
-{
- for (int r=0; r<wl.Size(); r++)
- {
- for (int i=0; i<wl.Size(); i++)
- os << wl.Element(r,i) ;
- os << endl;
- }
- return os;
-}
-
-/* Initialize as SU(3) generator
- * a=1...8
- */
-void WilsonLine::InitializeAsGenerator(int a)
-{
- if (a < 1 or a >8)
- {
- std::cerr << "SU(3) generator id must be within 1...8, asked " << a << std::endl;
- return;
- }
-
- /*
- std::complex<double> imag(0,1.0);
- std::vector< std::complex<double> > row;
- switch(a)
- {
- case 1:
- row.clear();
- row.push_back(0); row.push_back(0.5); row.push_back(0);
- data.push_back(row);
- row.clear();
- row.push_back(0.5); row.push_back(0); row.push_back(0);
- data.push_back(row);
- row.clear();
- row.push_back(0); row.push_back(0); row.push_back(0);
- data.push_back(row);
- break;
- case 2:
- row.clear();
- row.push_back(0); row.push_back(-0.5*imag); row.push_back(0);
- data.push_back(row);
- row.clear();
- row.push_back(0.5*imag); row.push_back(0); row.push_back(0);
- data.push_back(row);
- row.clear();
- row.push_back(0); row.push_back(0); row.push_back(0);
- data.push_back(row);
- break;
- case 3:
- row.clear();
- row.push_back(0.5); row.push_back(0); row.push_back(0);
- data.push_back(row);
- row.clear();
- row.push_back(0); row.push_back(-0.5); row.push_back(0);
- data.push_back(row);
- row.clear();
- row.push_back(0); row.push_back(0); row.push_back(0);
- data.push_back(row);
- break;
- case 4:
- row.clear();
- row.push_back(0); row.push_back(0); row.push_back(0.5);
- data.push_back(row);
- row.clear();
- row.push_back(0); row.push_back(0); row.push_back(0);
- data.push_back(row);
- row.clear();
- row.push_back(0.5); row.push_back(0); row.push_back(0);
- data.push_back(row);
- break;
- case 5:
- row.clear();
- row.push_back(0); row.push_back(0); row.push_back(-0.5*imag);
- data.push_back(row);
- row.clear();
- row.push_back(0); row.push_back(0); row.push_back(0);
- data.push_back(row);
- row.clear();
- row.push_back(0.5*imag); row.push_back(0); row.push_back(0);
- data.push_back(row);
- break;
- case 6:
- row.clear();
- row.push_back(0); row.push_back(0); row.push_back(0);
- data.push_back(row);
- row.clear();
- row.push_back(0); row.push_back(0); row.push_back(0.5);
- data.push_back(row);
- row.clear();
- row.push_back(0); row.push_back(0.5); row.push_back(0);
- data.push_back(row);
- break;
- case 7:
- row.clear();
- row.push_back(0); row.push_back(0); row.push_back(0);
- data.push_back(row);
- row.clear();
- row.push_back(0); row.push_back(0); row.push_back(-0.5*imag);
- data.push_back(row);
- row.clear();
- row.push_back(0); row.push_back(0.5*imag); row.push_back(0);
- data.push_back(row);
- break;
- case 8:
- row.clear();
- row.push_back(1.0/(2.0*std::sqrt(3))); row.push_back(0); row.push_back(0);
- data.push_back(row);
- row.clear();
- row.push_back(0); row.push_back(1.0/(2.0*std::sqrt(3))); row.push_back(0);
- data.push_back(row);
- row.clear();
- row.push_back(0); row.push_back(0); row.push_back(-2.0/(2.0*std::sqrt(3)));
- data.push_back(row);
- break;
- default:
- std::cerr << "Incorrect a!" << std::endl;
- };
- */
- cerr << "Generator to gsl_matrix not yet implemented! " << endl;
- exit(1);
-}
-
-void WilsonLine::InitializeAsIdentity()
-{
-
- for (int i=0; i<size; i++)
- {
- for (int j=0; j<size; j++)
- {
- if (i==j)
- gsl_matrix_complex_set(wline_gsl_matrix,i,j,gsl_complex_rect(1,0));
- else
- gsl_matrix_complex_set(wline_gsl_matrix,i,j,gsl_complex_rect(0,0));
- }
- }
-}
-
-
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_complex.h>
-#include <gsl/gsl_complex_math.h>
-
-/*
- * Calculate exp of matrix using GSL
- * Way too much memory allocation, so one should optimize if doing some
- * heavy numerics
- */
-WilsonLine WilsonLine::Exp()
-{
- std::cerr << "Test GSL matrix exponential before using this! WilsonLine WilsonLine::Exp() " << std::endl;
- gsl_matrix_complex* m = GetGslMatrix();
- gsl_matrix_complex* exp = gsl_matrix_complex_alloc(3,3);
- my_gsl_complex_matrix_exponential(exp,m,3);
- WilsonLine result;
- result.InitializeAsGslMatrix(exp);
-
- gsl_matrix_complex_free(m);
- gsl_matrix_complex_free(exp);
-
- return result;
-}
-
-
-gsl_matrix_complex* WilsonLine::GetGslMatrix()
-{
- return wline_gsl_matrix;
-}
-
-
-
-void my_gsl_complex_matrix_exponential(gsl_matrix_complex *eA, gsl_matrix_complex *A, int dimx)
-{
- int j,k=0;
- gsl_complex temp;
- gsl_matrix *matreal =gsl_matrix_alloc(2*dimx,2*dimx);
- gsl_matrix *expmatreal =gsl_matrix_alloc(2*dimx,2*dimx);
- //Converting the complex matrix into real one using A=[Areal, Aimag;-Aimag,Areal]
- for (j = 0; j < dimx;j++)
- for (k = 0; k < dimx;k++)
- {
- temp=gsl_matrix_complex_get(A,j,k);
- gsl_matrix_set(matreal,j,k,GSL_REAL(temp));
- gsl_matrix_set(matreal,dimx+j,dimx+k,GSL_REAL(temp));
- gsl_matrix_set(matreal,j,dimx+k,GSL_IMAG(temp));
- gsl_matrix_set(matreal,dimx+j,k,-GSL_IMAG(temp));
- }
-
- gsl_linalg_exponential_ss(matreal,expmatreal,.01);
-
- double realp;
- double imagp;
- for (j = 0; j < dimx;j++)
- for (k = 0; k < dimx;k++)
- {
- realp=gsl_matrix_get(expmatreal,j,k);
- imagp=gsl_matrix_get(expmatreal,j,dimx+k);
- gsl_matrix_complex_set(eA,j,k,gsl_complex_rect(realp,imagp));
- }
- gsl_matrix_free(matreal);
- gsl_matrix_free(expmatreal);
-}
-
-void WilsonLine::InitializeAsGslMatrix(gsl_matrix_complex* m)
-{
- gsl_matrix_complex_free(wline_gsl_matrix);
- wline_gsl_matrix = m;
-
-}
-
-
-
-
diff --git a/src/wilsonline.hpp b/src/wilsonline.hpp
deleted file mode 100644
index ad58728..0000000
--- a/src/wilsonline.hpp
+++ /dev/null
@@ -1,79 +0,0 @@
-/*
- * Diffraction at sub-nucleon scale
- * Wilson line / SU(3) matrix handling
- * Heikki Mäntysaari <mantysaari@bnl.gov>, 2015
- * Standalone class, no external dependences, please!
- */
-
-#ifndef wilsonline_hpp
-#define wilsonline_hpp
-
-#include <iostream>
-#include <complex>
-#include <vector>
-
-
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_linalg.h>
-
-class WilsonLine
-{
-public:
- WilsonLine(); // Initialize everything to 0
- WilsonLine(std::vector< std::vector< std::complex<double> > > &d);
- WilsonLine(gsl_matrix_complex *m);
- WilsonLine(const WilsonLine &m);
-
- ~WilsonLine();
-
- WilsonLine operator*(WilsonLine& w);
- WilsonLine operator*(std::complex<double> t);
- WilsonLine operator+(WilsonLine& w);
- WilsonLine& operator=(const WilsonLine& w);
-
-
-
- // Multiplies this by w^\dagger, returns the product
- // Fast in BLAS
- WilsonLine MultiplyByHermitianConjugate(WilsonLine& w);
-
-
- WilsonLine ComplexConjugate();
- WilsonLine Transpose();
- WilsonLine HermitianConjugate();
-
- void Set(int row, int column, std::complex<double> value);
-
- int Size(); // Size of NxN matrix
- std::complex<double> Element(int row, int col) const;
-
- std::complex<double> Trace();
-
- void InitializeAsGenerator(int a); // Initialize as Color matrix t^a
- void InitializeAsIdentity();
-
- WilsonLine Exp(); // Calculate exponential
-
- gsl_matrix_complex* GetGslMatrix();
-
- void InitializeAsGslMatrix(gsl_matrix_complex* m);
-
-private:
- gsl_matrix_complex *wline_gsl_matrix;
-
- static const int size = 3;
-};
-
-
-std::ostream& operator<<(std::ostream& os, WilsonLine& wl);
-
-
-// Matrix exponential, calculates as writing the matrix as 2Nx2N real matrix.
-// Downloaded from
-void my_gsl_complex_matrix_exponential(gsl_matrix_complex *eA, gsl_matrix_complex *A, int dimx);
-
-
-
-
-
-#endif /* wilsonline_hpp */
--
GitLab