Skip to content
Snippets Groups Projects
Commit 44b667e5 authored by Pyry Runko (Puck)'s avatar Pyry Runko (Puck)
Browse files

building the integral, not finished!

parent 54bab6d1
Branches feature-idea-support
No related tags found
No related merge requests found
......@@ -624,11 +624,14 @@ double IPGlasma::SumAmplitudes(int l)
//TODO: check
double IPGlasma::SumAmplitudes()
{
double sum = 0.0;
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;
//int thrownAway = 0;
#pragma omp parallel for reduction(+:sum)
for(int i =0; i < L; i += 1){
for(int j =0; j < L; j += 1){
......@@ -636,35 +639,22 @@ double IPGlasma::SumAmplitudes()
int q1[2] = {i,j};
quark = GetWilsonLine(WilsonLineCoordinate(i, j));
double r1 = DistanceToOrigin(q1);
//#pragma omp parallel for reduction(+:sum) collapse(2)
for(int k =0; k < L; k +=1){
for(int l =0; l < L; l +=1){
int q2[2] = {k,l};
double r2 = DistanceToOrigin(q2);
if (r1 > 1. && r2 > 1.){
// cout << "#!!Throwing away large dipole separation!" << endl;
//cout << "#Quark coordinate indices " << q1[0] << "," << q1[1] << " and " << q2[0] << "," << q2[1] << endl;
//cout << "#Quark coordinates " << X(q1[0]) << "," << Y(q1[1]) << " and " << X(q2[0]) << "," << Y(q2[1]) << endl;
//cout << "#Dipole size is " << r1-r2 << endl;
//cout << "#Dipole amplitude is " << Amplitude(0, q1, q2) << endl;
// thrownAway += 1;
continue;
}
r = sqrt(pow((X(i)-X(k)),2)+pow(Y(j)-Y(l)),2);
// Haw w2 = WilsonLine(k,l)
// Laske w1.MultiplbyHermitianConjugate(w2);
sum += Amplitude(0, q1, q2);
//sum += Amplitude(0, quark, q2);
double WF = WFGammaqqbarLongitudinal(double r, double z, double Q, struct* quark_flavor);
sum += Amplitude(0, quark, q2)*WF*WF;
}
}
}
}
cout << "# total sum = " << sum << endl;
// cout << "# i counted " << iWasHere << " values" << endl;
// cout << "# threw away " << thrownAway << " values (" << 100.*thrownAway*pow(L,-4) << "%)" << endl;
return sum;
}
......@@ -672,25 +662,34 @@ double IPGlasma::SumAmplitudes()
// return longitudinal photon polarization wave function
// source arxiv.org/pdf/hep-ph/0606272
// quark flavor f: 0=u, 1=d
double IPGlasma::WFGammaqqbar(double r, int f)
double IPGlasma::WFGammaqqbarLongitudinal(double r, double z, double Q, struct* quark_flavor)
{
//Assuming 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
double lattice_spacing = xcoords[1]-xcoords[0];
// fine structure constant set to 1/137
double e = cmath.sqrt(4*cmath.pi*1/137)
double m = 0.16e-9 //GeV
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};
if (f){ e_f = -1/3;}
else { e_f = 2/3;}
double z = 0.5 //GeV^0
double QQ = 10 //GeV^2
double Delta = 0.5 //GeV
double l = a/100 // << a^-1 (a = lattice const)
double eps = cmath.sqrt(z*(1-z)*QQ*m_f*m_f)
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);
return e_f*e*cmath.sqrt(NC)*2*Q*z*(1-z)*std::cyl_bessel_k(0, eps*r)0.5/cmath.pi
return (e_f*e*sqrt(NC)*2*Q*Q*z*(1-z)*std::cyl_bessel_k(0, eps*r)*0.5/M_PI);
}
// TODO: This has to be done inside the loop in SumAmplitudes
double IPGlasma::DiffractiveCrossSection()
{
SumAmplitudes
return 0;
}
......@@ -61,7 +61,14 @@ public:
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);
//TODO: missä paras asettaa arvot?
struct q_flavor{
double e_charge;
double mass;
}
private:
......
......@@ -55,66 +55,6 @@ int printDipoleAmplitudes(string fname, IPGlasma ipglasma)
return 0;
}
int thread_id = omp_get_thread_num();
int num_threads = omp_get_num_threads();
// Print the number of threads (only from the master thread to avoid too many prints)
if (thread_id == 0) {
cout << "# Number of threads: " << num_threads << endl;
}
// Print the thread id for debugging (can be limited to a few threads for easier reading)
if (thread_id < 5) { // Change 5 to a larger number if needed
cout << "# Thread " << thread_id << " is starting" << endl;
}
}
cout << "# Hello from main.cpp sumAmplitudes!" << endl;
double sum = 0.0;
int L = 512; //ipglasma->xcoords.size();
cout << "# Summing with lattice size L = " << L << endl;
//WilsonLine quark;
//int thrownAway = 0;
//int iWasHere = 0;
#pragma omp parallel for reduction(+:sum)//,thrownAway,iWasHere) //collapse(4)
for(int i =0; i < L; i += 1){
for(int j =0; j < L; j += 1){
//#pragma omp parallel for reduction(+:sum) collapse(2)
for(int k =0; k < L; k +=1){
// cout << "# Thread " << omp_get_thread_num() << " is processing (" << i << ", " << j << ", " << k <<")" << endl;
for(int l =0; l < L; l +=1){
int q1[2] = {i,j};
double r1 = ipglasma->DistanceToOrigin(q1);
int q2[2] = {k,l};
//double r1 = ipglasma->DistanceToOrigin(q1);
double r2 = ipglasma->DistanceToOrigin(q2);
if (r1 > 1. && r2 > 1.){
// cout << "#!!Throwing away large dipole separation!" << endl;
//cout << "#Quark coordinate indices " << q1[0] << "," << q1[1] << " and " << q2[0] << "," << q2[1] << endl;
//cout << "#Quark coordinates " << X(q1[0]) << "," << Y(q1[1]) << " and " << X(q2[0]) << "," << Y(q2[1]) << endl;
//cout << "#Dipole size is " << r1-r2 << endl;
//cout << "#Dipole amplitude is " << Amplitude(0, q1, q2) << endl;
//thrownAway += 1;
continue;
}
// Haw w2 = WilsonLine(k,l)
// Laske w1.MultiplbyHermitianConjugate(w2);
sum += ipglasma->Amplitude(0, q1, q2);
//sum += Amplitude(0, quark, q2);
//iWasHere += 1;
}
}
}
}
cout << "# total sum = " << sum << endl;
//cout << "# i counted " << iWasHere << " values" << endl;
//cout << "# threw away " << thrownAway << " values (" << 100.*thrownAway*pow(L,-4) << "%)" << endl;
return sum;
}
int main(int argc, char* argv[])
{
......@@ -125,20 +65,20 @@ int main(int argc, char* argv[])
string fname = argv[1];
IPGlasma ipglasma(fname, 0, BINARY);
int task_id = 0;
if (!argv[2]){
printDipoleAmplitudes(fname, ipglasma);
}
int task_id = 0;
if (argc > 3){
task_id= atoi(argv[3]);
cout << "#Integrating array job #" << task_id << endl;
ipglasma->SumAmplitudes(task_id);
ipglasma.SumAmplitudes(task_id);
}
else if (!argv[2]){
printDipoleAmplitudes(fname, ipglasma);
}
else{
double sum = ipglasma->SumAmplitudes();
double sum = ipglasma.SumAmplitudes();
// Save results
//std::filesystem::create_directory("output_"+fname);
ofstream outfile("da_sum_ +" + fname + ".txt");
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment