#include "SHRiMPS/Cross_Sections/Sigma_Elastic.H" #include "SHRiMPS/Tools/Special_Functions.H" #include "ATOOLS/Math/Gauss_Integrator.H" #include "ATOOLS/Math/Random.H" #include "ATOOLS/Org/MyStrStream.H" #include "ATOOLS/Org/Run_Parameter.H" #include "ATOOLS/Org/Exception.H" using namespace SHRIMPS; using namespace ATOOLS; double Sigma_Elastic::dSigma_dt::operator()(double B) { // factors: // * B from B integration in polar coordinates // * 2 pi J0 from integral of angle between Q and B // * (arg) usual argument of elastic scattering return B * 2.*M_PI*SF.Jn(0,B*m_Q) * p_sigma_el->GetDiffArgument(B); } Sigma_Elastic::Sigma_Elastic() : m_tmin(0), m_tmax(1.), m_summed(0.), m_steps(1000), m_delta((m_tmax-m_tmin)/m_steps) { } double Sigma_Elastic::GetValue(const double & B) { return ATOOLS::sqr(p_eikonal->Prefactor()*(1.-exp(-(*p_eikonal)(B)/2.))); } double Sigma_Elastic::GetCombinedValue(const double & B) { return ATOOLS::sqr(GetDiffArgument(B)); } double Sigma_Elastic::GetDiffArgument(const double & B) { double value(0.); for (size_t i=0;isize();i++) { for (size_t j=0;j<(*p_eikonals)[i].size();j++) { Omega_ik * eikonal = (*p_eikonals)[i][j]; value += eikonal->Prefactor()*(1.-exp(-(*eikonal)(B)/2.)); } } return value; } void Sigma_Elastic::FillGrids() { m_diffgrid.clear(); m_intgrid.clear(); FillDiffQGrid(); FillIntQGridAndNormalize(); } void Sigma_Elastic::FillDiffQGrid() { msg_Out()<Picobarn()/(4.*M_PI) * ATOOLS::sqr(integrator.Integrate(0.,MBpars.GetEikonalParameters().bmax, MBpars.GetEikonalParameters().accu,1.)); if (dabs(value<0.)) value = 0.; m_diffgrid.push_back(value); } } void Sigma_Elastic::FillIntQGridAndNormalize() { m_intgrid.push_back(0.); m_summed = 0.; double average, binvalue; for (size_t i=1;iGet()); unsigned int i(0); while (random-m_intgrid[i]>=0) i++; return m_tmin+(i-1)*m_delta + m_delta *(random-m_intgrid[i-1])/(m_intgrid[i]-m_intgrid[i-1]); } double Sigma_Elastic::Test() { const Eikonal_Parameters & eikparams(MBpars.GetEikonalParameters()); const FormFactor_Parameters & ffparams(MBpars.GetFFParameters()); const double EulerGamma= 0.577215664901532860606512090082 ; double a(ffparams.Lambda2/(8.*(1.+ffparams.kappa))); double c(eikparams.beta02*ffparams.Lambda2*(1.+ffparams.kappa)* exp(2.*eikparams.Delta*eikparams.Ymax)/(8.*M_PI)); double alpha(2.*M_PI*ffparams.norm); ExpInt expint; double ei(expint.GetExpInt(-c)), ei2(expint.GetExpInt(-c/2.)); return alpha*(EulerGamma+ei-ei2+log(c/4.))/(2.*a)*rpa->Picobarn(); }