Commit 5dad317438857f356a7a0dfeda472d04346cd085

  • avatar
  • Heikki Mäntysaari <heikki.mantysaari @j…u.fi> (Committer)
  • Thu Jul 14 08:46:22 EEST 2011
  • avatar
  • Heikki Mäntysaari <heikki.mantysaari @j…u.fi> (Author)
  • Thu Jul 14 08:46:22 EEST 2011
New feature to calculate saturation scale in r-space via definition
N(r=1/Q_s)=const
src/amplitude.cpp
(63 / 0)
  
2020#include <gsl/gsl_min.h>
2121#include <gsl/gsl_roots.h>
2222#include "interpolation.hpp"
23#include "hankel.hpp"
2324
2425
2526#include <cmath>
485485
486486 return std::sqrt(res);
487487
488}
489
490/*
491 * Solve saturation scale in r space
492 * Defined as N(r=1/Q_s) = const
493 */
494const REAL SATSCALE_R = 0.2;
495struct SatscaleRHelper
496{
497 Hankel *transf;
498 REAL y;
499};
500
501REAL SatscaleRHelperf(REAL lnr, void* p)
502{
503 SatscaleRHelper* par = (SatscaleRHelper*)p;
504 return par->transf->Amplitude_r(std::exp(lnr), par->y) - SATSCALE_R;
505}
506
507REAL Amplitude::SaturationScaleR(REAL y)
508{
509 const int MAX_ITER = 30;
510 const REAL ROOTFINDACCURACY = 0.05;
511
512 Hankel transform(this);
513 SatscaleRHelper help; help.transf = &transform;
514 help.y=y;
515 gsl_function f; f.params=&help;
516 f.function=SatscaleRHelperf;
517
518 const gsl_root_fsolver_type *T = gsl_root_fsolver_bisection;
519 gsl_root_fsolver *s = gsl_root_fsolver_alloc(T);
520 REAL minlnr = std::log(1e-5);
521 REAL maxlnr = std::log(10);
522 gsl_root_fsolver_set(s, &f, minlnr, maxlnr);
523
524 int iter=0; int status; REAL min,max;
525 do
526 {
527 iter++;
528 gsl_root_fsolver_iterate(s);
529 min = gsl_root_fsolver_x_lower(s);
530 max = gsl_root_fsolver_x_upper(s);
531 status = gsl_root_test_interval(min, max, 0, ROOTFINDACCURACY);
532 //cout << "y=" << y <<", interval " << std::exp(min) << " - " << std::exp(max) << endl;
533
534 } while (status == GSL_CONTINUE and iter < MAX_ITER);
535
536 REAL res = gsl_root_fsolver_root(s);
537
538 if (iter>=MAX_ITER)
539 {
540 cerr << "Solving failed at y=" << y << " at " << LINEINFO <<
541 " result " << std::exp(res) << " relaccuracy "
542 << (std::exp(max)-std::exp(min))/std::exp(res) << endl;
543 }
544
545
546 gsl_root_fsolver_free(s);
547
548 return std::exp(-res);
549
488550}
489551
490552REAL Amplitude::InitialCondition(REAL ktsqr)
src/amplitude.hpp
(1 / 0)
  
5454 void IntializeBSpline(int ktsqrind, REAL rapidity);
5555
5656 REAL SaturationScale(REAL y);
57 REAL SaturationScaleR(REAL y);
5758 REAL SolveKtsqr(REAL y, REAL amp);
5859
5960 void AddDataPoint(int ktsqrindex, int yindex, REAL value, REAL der);
src/interpolation.cpp
(1 / 1)
  
120120
121121REAL Interpolator::Derivative(REAL x)
122122{
123 REAL res, err; int status;
123 REAL res; int status;
124124 switch(method)
125125 {
126126 case INTERPOLATE_SPLINE:
src/main.cpp
(21 / 1)
  
4141 GENERATE_SINGLE_RPLOT, // Print amplitude in r-space at a given rapidity
4242 LOGLOG_DERIVATIVE, // Calculate d ln N(k^2) / d ln(k^2)
4343 SATURATION_SCALE, // Calculate Q_s as a function of y up to maxy
44 SATURATION_SCALE_R, // Calculate N(r=1/Q_s) = const saturation scale
4445 PT_SPECTRUM, // Print dN_{ch} / dyd^2p_T as a function of p_T
4546 PSEUDOY_SPECTRUM, // dN_{ch} / d\eta, pseudorapidity spectrum
4647 UGD // Plot uninteg. gluon density as a function of k_T
105105void SinglePlotR(); // FT to r-space
106106void Clear();
107107void SaturationScale();
108void SaturationScaleR();
108109void ParticleSpectrum_pt(); //dN_ch / dyd2pt as a function of sqrt(s)
109110void ParticleSpectrum_pseudoy();
110111void UnintegratedGluonDistribution(); // Plot unintegrated gluon density as a function of k_T
127127 {
128128 cout << "Usage: " << endl;
129129 cout << "-mode [MODE]: what to do, modes: GENERATE_DATA, GENERATE_PLOTS, SINGLE_PLOT, SINGLE_RPLOT, LOGLOG_DERIVATIVE " << endl;
130 cout << " SATURATION_SCALE PT_SPECTRUM PSEUDOY_SPECTRUM UGD" << endl;
130 cout << " SATURATION_SCALE[_R] PT_SPECTRUM PSEUDOY_SPECTRUM UGD" << endl;
131131 cout << "-method [METHOD]: what method is used to solve BK, methods: BRUTEFORCE, CHEBYSHEV" << endl;
132132 cout << "-output [prefix]: set output file prefix, filenames are prefix_y[rapidity].dat" << endl;
133133 cout << "-miny, -maxy: rapidity values to solve" << endl;
258258 mode=GENERATE_SINGLE_RPLOT;
259259 else if (string(argv[i+1])=="SATURATION_SCALE")
260260 mode = SATURATION_SCALE;
261 else if (string(argv[i+1])=="SATURATION_SCALE_R")
262 mode = SATURATION_SCALE_R;
261263 else if (string(argv[i+1])=="PT_SPECTRUM")
262264 mode = PT_SPECTRUM;
263265 else if (string(argv[i+1])=="PSEUDOY_SPECTRUM")
451451 SinglePlotR();
452452 else if (mode==SATURATION_SCALE)
453453 SaturationScale();
454 else if (mode==SATURATION_SCALE_R)
455 SaturationScaleR();
454456 else if (mode==PT_SPECTRUM)
455457 ParticleSpectrum_pt();
456458 else if (mode==PSEUDOY_SPECTRUM)
663663 delete[] yarray;
664664 delete[] qsarray;
665665 }
666}
667
668/*
669 * Calculate Saturation scale N(r=1/Q_s) = const
670 */
671void SaturationScaleR()
672{
673 cout << "# Saturation scale N(r=1/Q_s) = const" << endl;
674 cout << "# y Q_s [1/GeV]" << endl;
675 for (REAL y=0.5; y<maxy; y+=0.5)
676 {
677 cout << y << " " << N->SaturationScaleR(y) << endl;
678 }
679
666680}
667681
668682// dN_{ch} / dydp^2