Commit c38f7db53e5f24d4703387f6bb4346d238f980ed

  • avatar
  • Heikki Mäntysaari <heikki.mantysaari @j…u.fi> (Committer)
  • Wed Jun 22 09:03:42 EEST 2011
  • avatar
  • Heikki Mäntysaari <heikki.mantysaari @j…u.fi> (Author)
  • Wed Jun 22 09:03:42 EEST 2011
* Code to solve saturation scale
* Support for multiple running coupling prescriptions
src/amplitude.cpp
(5 / 5)
  
4040 interpolation_rapidity=-1.0;
4141 interpolation_points = INTERPOLATION_POINTS;
4242
43 running_coupling = false;
43 running_coupling = CONSTANT;
4444}
4545
4646/*
324324{
325325 Sathelper helper;
326326 helper.y=y; helper.N=this;
327 //helper.gammac = 0.6275;
328 helper.gammac = 0.5;
327 helper.gammac = 0.6275;
328 //helper.gammac = 0.5;
329329
330330 gsl_function f; f.function=&SaturationHelperf;
331331 f.params=&helper;
548548
549549}
550550
551void Amplitude::SetRunningCoupling(bool rc)
551void Amplitude::SetRunningCoupling(RUNNING_COUPLING rc)
552552{
553553 running_coupling=rc;
554554}
555555
556bool Amplitude::RunningCoupling()
556RUNNING_COUPLING Amplitude::RunningCoupling()
557557{
558558 return running_coupling;
559559}
src/amplitude.hpp
(11 / 4)
  
3535 INVPOWER4, // 1/(k^4+1), very arbitrary
3636 GAUSS // Exp[-(Log[k^2] + 2)^2/5], hep-ph/0110325
3737};
38
3938
39enum RUNNING_COUPLING
40{
41 PARENT_DIPOLE, // scale of the parent dipole
42 MAXK, // max of {k,k'}, in NL term just k
43 CONSTANT // no running
44};
45
46
4047class Amplitude
4148{
4249 public:
8686
8787 void SetInterpolationPoints(int p);
8888
89 void SetRunningCoupling(bool rc);
90 bool RunningCoupling();
89 void SetRunningCoupling(RUNNING_COUPLING rc);
90 RUNNING_COUPLING RunningCoupling();
9191
9292 unsigned int YPoints();
9393 unsigned int KtsqrPoints();
122122 REAL maxy;
123123 REAL delta_y;
124124
125 bool running_coupling;
125 RUNNING_COUPLING running_coupling;
126126
127127 bool datafile; // True if data is read from external file
128128
src/main.cpp
(57 / 4)
  
1313#include "tools.hpp"
1414#include "chebyshev.hpp"
1515#include "hankel.hpp"
16#include "interpolation.hpp"
1617#include <gsl/gsl_errno.h>
1718#include <cmath>
1819#include <fstream>
8181string datafile="output";
8282int avg=0;
8383string file_prefix="output";
84bool rungekutta = false;
8485
8586// Parameters for ChebyshevAmplitudeSolver
8687int chebyshev_degree=100;
123123 //cout << "-minktsqr, -maxktsqr: range of k_T^2 to plot, doesn't affect to limits when solving BK" << endl;
124124 cout << "-ic [initial condition]: set initial condition, possible ones are FTIPSAT, INVPOWER, INVPOWER4, GAUSS " << endl;
125125 cout << "-kc: apply kinematical constraint" << endl;
126 cout << "-rungekutta: use 2nd order runge kutta (only with method BRUTEFORCE)" << endl;
126127 cout << "-avg [avgs]: number or averagements" << endl;
127128 cout << "-data [datafile]: read data from datafiles from path datafile_y[rapdity].dat" << endl;
128129 cout << " -maxdatay [yval]: set maximum y value for datafiles, -delta_datay [yval] difference of yvals for datafiles" << endl;
154154 kc=true;
155155 else if (string(argv[i])=="-rc")
156156 running_coupling=true;
157 else if (string(argv[i])=="-rungekutta")
158 rungekutta = true;
157159 else if (string(argv[i])=="-scale_sat")
158160 scale_sat=true;
159161 else if (string(argv[i])=="-avg")
258258 else if (method==CHEBYSHEV_SERIES) N = new ChebyshevAmplitudeSolver;
259259
260260 N->SetInitialCondition(ic);
261 N->SetRunningCoupling(running_coupling);
261 N->SetRunningCoupling(PARENT_DIPOLE);
262262 N->SetKinematicConstraint(kc);
263263 N->SetMaxY(maxy);
264264 N->SetMaxKtsqr(maxktsqr);
334334 }
335335 else
336336 {
337 if (rungekutta)
338 {
339 ((BruteForceSolver*)N)->SetRungeKutta(true);
340 infostr << "# 2nd order Runge Kutta is applied" << endl;
341 cout << "# 2nd order Runge Kutta is applied" << endl;
342 }
337343 cout << "# Generating data..." << endl;
338344 ((BruteForceSolver*)N)->Solve(maxy);
339345 }
555555void SaturationScale()
556556{
557557 cout << "# Saturation scale Q_s as a function of y" << endl;
558 cout << "# y Q_s " << endl;
559 for (REAL tmpy=0; tmpy<maxy; tmpy += 0.1)
558 cout << "# y Bspline-Q_s Q_s " << endl;
559
560 // Bspline interpolation
561 int points = static_cast<int>(maxy/0.1);
562 REAL *yarray = new REAL[points];
563 REAL *qsarray = new REAL[points];
564
565 const int interpolation_points=80;
566 int interpolation_start, interpolation_end;
567
568 for (int i=0; i<points; i++)
560569 {
561 cout << std::scientific << std::setprecision(15) << tmpy << " " << N->SaturationScale(tmpy) << endl;
570 if (i-interpolation_points/2<0)
571 {
572 interpolation_start=0; interpolation_end=interpolation_points-1;
573 }
574 else if (i+interpolation_points/2 > points)
575 {
576 interpolation_end = points-1;
577 interpolation_start = points-1-interpolation_points;
578 }
579 else
580 {
581 interpolation_start = i-interpolation_points/2;
582 interpolation_end = i+interpolation_points/2;
583 }
584 int interpo_points = interpolation_end - interpolation_start + 1;
585
586 REAL *yarray = new REAL[interpo_points];
587 REAL *qsarray = new REAL[interpo_points];
588 for (int j=interpolation_start; j<=interpolation_end; j++)
589 {
590 yarray[j-interpolation_start] = j*0.1;
591 qsarray[j-interpolation_start] = N->SaturationScale(j*0.1);
592 }
593
594 Interpolator inter(yarray, qsarray, interpo_points);
595 inter.SetMethod(INTERPOLATE_BSPLINE);
596 inter.Initialize();
597
598 REAL tmpy = 0.1*i;
599 cout << std::scientific << std::setprecision(15) << tmpy << " "
600 << inter.Evaluate(tmpy) << " " << N->SaturationScale(tmpy) << endl;
601
602 delete[] yarray;
603 delete[] qsarray;
562604 }
563605}
src/solver_force.cpp
(38 / 6)
  
55
66#include "amplitude.hpp"
77#include "datafile.hpp"
8#include <algorithm>
89#include "solver_force.hpp"
910#include <vector>
1011#include <gsl/gsl_sf_gamma.h>
2929 REAL ktsqr;
3030 Amplitude* N;
3131 REAL y;
32 REAL offset;
3233};
3334
3435
5050 * sgn(k'^2 - k^2)/k^2 [N(k^2) + k'^2 N'(k^2) ] + N(k^2)/( Sqrt[5] k^2 )
5151 */
5252
53 REAL kn = par->ktsqr * (par->N->N(par->ktsqr, par->y) - par->offset );
54
5355 if (std::abs(ktsqr - par->ktsqr) < 1e-15)
5456 {
5557 //cerr << "ktsqr \\approx par->ktsqr and we can't handle this! y=" << par->y
6565 */
6666 }
6767 else
68 result += (ktsqr*par->N->N(ktsqr, par->y) - par->ktsqr*par->N->N(par->ktsqr, par->y))
68 result += (ktsqr*par->N->N(ktsqr, par->y) - kn)
6969 / std::abs(ktsqr - par->ktsqr);
7070
71 result += par->ktsqr * par->N->N(par->ktsqr, par->y) / sqrt( 4.0*SQR(ktsqr) + SQR(par->ktsqr));
71 result += kn / sqrt( 4.0*SQR(ktsqr) + SQR(par->ktsqr));
7272
7373 result /= ktsqr;
74
7475
76
7577 //cout << ktsqr << " " << result << endl;
7678 return result;
7779}
117117 return result;
118118}
119119
120REAL BruteForceSolver::RapidityDerivative(REAL ktsqr, REAL y)
120/*
121 * Calculate \partial_Y N
122 * offset is used in 2nd order Runge Kutta where the derivative
123 * \partial_Y N(Y) = F(Y, N(Y)) is evaluated at point
124 * F(Y-\delta Y, N(Y) - \delta Y F(Y, N(Y)) )
125 * It's default value is 0
126 * NB: It can't be used with the kinematical constraint as in that case the
127 * equation is not local in Y
128 *
129 */
130REAL BruteForceSolver::RapidityDerivative(REAL ktsqr, REAL y, REAL offset)
121131{
122 REAL alphabar=0.2;
123 if (RunningCoupling())
132 REAL alphabar=1.0;
133 if (RunningCoupling()==CONSTANT)
134 alphabar = 0.2;
135 if (RunningCoupling()==PARENT_DIPOLE)
124136 alphabar = Alpha_s(ktsqr)*Nc/M_PI;
125137
126138 inthelper_bkmom inthelp;
127139 inthelp.N=this;
128140 inthelp.y=y;
129141 inthelp.ktsqr = ktsqr;
142 inthelp.offset=offset;
130143 gsl_function int_helper;
131144
132145 if (kinematic_constraint==false)
212212 newn = n[ktsqrind][yind-1] + dy*tmpder
213213 + 1.0/2.0*dy*( tmpder - old_der);
214214 }
215
216 if (rungekutta and yind>1)
217 {
218 REAL der2 = RapidityDerivative(tmpkt, yvals[yind-1] - dy,
219 - dy*tmpder);
220 newn = n[ktsqrind][yind-1]
221 + 3.0/2.0*dy*tmpder - 1.0/2.0*dy*der2;
222 }
215223
216224
217225 AddDataPoint(ktsqrind, yind, newn, tmpder );
271271
272272BruteForceSolver::BruteForceSolver()
273273{
274
274 rungekutta = false;
275275}
276276BruteForceSolver::~BruteForceSolver()
277277{
278278
279}
280
281void BruteForceSolver::SetRungeKutta(bool rk)
282{
283 rungekutta = rk;
279284}
src/solver_force.hpp
(3 / 2)
  
1616{
1717 public:
1818 void Solve(REAL maxy);
19 REAL RapidityDerivative(REAL ktsqr, REAL y);
19 REAL RapidityDerivative(REAL ktsqr, REAL y, REAL offset=0);
2020
21 void SetRungeKutta(bool rk);
2122
2223 BruteForceSolver();
2324 ~BruteForceSolver();
2425 private:
25
26 bool rungekutta;
2627
2728
2829};
src/solver_force2.cpp
(1 / 1)
  
112112REAL BruteForceSolver2::RapidityDerivative(REAL u, REAL y)
113113{
114114 REAL alphabar=0.2;
115 if (RunningCoupling())
115 if (RunningCoupling()==PARENT_DIPOLE)
116116 alphabar = Alpha_s(Ktsqr(u))*Nc/M_PI;
117117
118118 inthelper_bkmom2 inthelp;