Commit 679a3dc595b381909ffc6fa9d394b7709bf8e71c

  • avatar
  • Heikki Mäntysaari <heikki.mantysaari @j…u.fi> (Committer)
  • Tue Jul 05 13:55:06 EEST 2011
  • avatar
  • Heikki Mäntysaari <heikki.mantysaari @j…u.fi> (Author)
  • Tue Jul 05 13:55:06 EEST 2011
* FTIPSAT ic follows HERA parametrization
* Some small changes to Hankel transformation code
src/amplitude.cpp
(20 / 3)
  
77#include "datafile.hpp"
88//#include <bci.h>
99#include <vector>
10#include <sstream>
1011#include <gsl/gsl_sf_gamma.h>
1112#include <gsl/gsl_integration.h>
1213#include <gsl/gsl_sf_bessel.h>
511511
512512REAL Amplitude::InitialCondition(REAL ktsqr)
513513{
514 int status;
514515 switch(ic)
515516 {
516517 case INVPOWER:
517518 return pow(ktsqr+1, -1); // BK in Full Momentum Space, hep-ph/0504080
518519 break;
519520 case FTIPSAT:
520 if (ktsqr > 2500) return 0.0;
521 return gsl_sf_gamma_inc(0.0,ktsqr/4.0); // Ft of 2(1-exp(-r^2))
521 // Ft of 1-exp(-r^2*Q_s0^2 / 4) = 1/2*Gamma[0, k^2/Q_s0^2]
522 // Fitted to HERA data in arXiv:0902.1112: Q_s0^2 = 0.24GeV^2
523 // Q_s0^2 defined in amplitude.hpp
524 // NB: In ref. \alpha(s) depends on log(4C^2/r^2lambdaqcd^2), but
525 // tools.cpp:Alpha_s depends on log(Q^2/lambdaqcd^2). Fitted value
526 // for C^2 is 5.3
527
528 gsl_sf_result res;
529 status = gsl_sf_gamma_inc_e(0.0, ktsqr/Q0SQR, &res);
530 if (status==15) return 0.0; // Overflow
531 if (status)
532 cerr << "gsl_sf_gamma_inc_e failed at " << LINEINFO << " with "
533 << "code " << status << " res " << 0.5*res.val << endl;
534 return 0.5*res.val;
522535 break;
523536 case INVPOWER4:
524537 return pow(SQR(ktsqr)+1, -1);
557557 {
558558 return "Data is read from a file, don't know what initial condition was used";
559559 }
560 std::stringstream s;
560561 switch (ic)
561562 {
562563 case INVPOWER:
563564 return "(kt^2 + 1)^(-1), BK in full momentum space, hep-ph/0504080";
564565 break;
565566 case FTIPSAT:
566 return "Gamma[0, kt^2/4], FT of 2(1-exp(-r^2))";
567 s << "0.5*Gamma[0, kt^2/Q_s0^2], FT of 1-exp(-r^2 Q_s0^2/4),"
568 << " Q_s0^2 = " << Q0SQR << " GeV^2";
569 return s.str();
567570 break;
568571 case INVPOWER4:
569572 return "(kt^4+1)^(-1), arbitrary";
src/amplitude.hpp
(2 / 0)
  
3030const int INTERPOLATION_POINTS = 20;
3131const int INTERPOLATION_POINTS_DER=50; // 50 good if 2000 ktsqrpoints, 100 for 5000
3232
33const REAL Q0SQR = 0.24; // GeV^2 arXiv:0902.1112
34
3335enum INITIAL_CONDITION
3436{
3537 FTIPSAT, // \int d^2r 1/(2\pi r^2) exp(ik.r) 2(1-exp(-r^2))
src/config.hpp
(2 / 1)
  
2525
2626// Reqularization of the running coupling
2727const REAL MAXALPHA = 0.3; // Berger&Stasto 1010.0671
28const REAL ALPHAS = 0.2; // \alpha_s if RC=constant
2829
2930// Other constants
3031
31const REAL eps=0.000001; // Epsilon for runge kutta
32const REAL eps=0.000001;
3233
3334// Inline functions
3435
src/hankel.cpp
(18 / 6)
  
2424REAL Inthelperf_hankel(REAL k, void* p)
2525{
2626 Inthelper_hankel* par = (Inthelper_hankel*)p;
27 return k*gsl_sf_bessel_j0(k* par->r) * par->N->N(SQR(k), par->y, true);
27 return k*gsl_sf_bessel_j0(k* par->r) * par->N->N(SQR(k), par->y, false);
2828}
2929
3030void Hankel::PrintRAmplitude()
3131{
3232 cout << "# r N" << endl;
33/*
34 for (int i=0; i<points; i++)
35 {
36 REAL tmpr = gsl_dht_k_sample(dht, i);
37 cout << tmpr << " " << SQR(tmpr)*transformed[i] << endl;
3338
39 }
40 return;
41 */
42
43
3444 REAL minr = 1e-6; REAL maxr=1e1;
3545 REAL mult = std::pow(maxr/minr, 1.0/((REAL)points) );
3646
47 #pragma omp parallel for
3748 for (int i=0; i<points; i++)
3849 {
3950 REAL r = minr*std::pow(mult, i);
8686
8787 sample = new REAL[points];
8888 transformed = new REAL[points];
89 /*
90 dht = gsl_dht_new(points, 0, N->Ktsqrval(N->KtsqrPoints()-2) );
89/*
90 REAL maxk = std::sqrt( N->KtsqrPoints()-2 );
91 dht = gsl_dht_new(points, 0, maxk );
9192
9293 for (int i=0; i<points; i++)
9394 {
94 REAL tmpktsqr = gsl_dht_x_sample(dht, i);
95 sample[i]=N->N(tmpktsqr, y);
95 REAL tmpkt = gsl_dht_x_sample(dht, i);
96 sample[i]=N->N(SQR(tmpkt), y);
9697 }
9798
9899 gsl_dht_apply( dht, sample, transformed );
99 */
100*/
100101}
101102
102103Hankel::~Hankel()
src/main.cpp
(1 / 1)
  
562562void SinglePlotR()
563563{
564564 cout << "# y=" << y << ", ic=" << N->InitialConditionStr() << endl;
565 Hankel transformed(N, y, 100);
565 Hankel transformed(N, y, 300);
566566 transformed.PrintRAmplitude();
567567
568568}
src/solver_force.cpp
(8 / 3)
  
181181 abserr << " relerror: " << abserr/result << endl;
182182
183183 // Nonlinear term
184 result -= SQR(N(ktsqr, y));
184 // if rc is {MIN,MAX}k, then result is already multiplied by alpha_s
185 if (RunningCoupling()==MAXK or RunningCoupling() == MINK)
186 result -= SQR(N(ktsqr, y))*Alphabar_s(ktsqr);
187 else
188 result -= SQR(N(ktsqr, y));
185189
186190 if (RunningCoupling()==CONSTANT)
187 result*=0.2;
191 result*=ALPHAS;
188192 if (RunningCoupling()==PARENT_DIPOLE)
189193 result *= Alphabar_s(ktsqr);
190194
282282 << h << endl;
283283 for (int ktsqrind = 0; ktsqrind < KtsqrPoints(); ktsqrind++)
284284 {
285 n[ktsqrind][yind] = amplitude[ktsqrind];
285 AddDataPoint(ktsqrind, yind, amplitude[ktsqrind], 0.0);
286 //n[ktsqrind][yind] = amplitude[ktsqrind];
286287 }
287288 continue;
288289 }
src/tools.cpp
(5 / 2)
  
3838
3939 // Errors related to convergence of integrals are handled when
4040 // gsl_integration functions are called, don't do anything with them here
41 // 14 = failed to reach tolerance
42 // 18 = roundoff error prevents tolerance from being achieved
41 // 14 = failed to reach tolerance
42 // 18 = roundoff error prevents tolerance from being achieved
4343 if (gsl_errno == 14 or gsl_errno == 18)
4444 return;
45
46 // 15 in expint.c:363: gsl_sf_gamma_inc underflows
47 if (gsl_errno == 15 and string(file)=="expint.c" and line==363) return;
4548
4649 errors++;
4750 std::cerr << file << ":"<< line <<": Error " << errors << ": " <<reason