#include "SHRiMPS/Eikonals/Single_Channel_Eikonal.H" #include "ATOOLS/Math/Gauss_Integrator.H" #include "ATOOLS/Math/Random.H" #include "ATOOLS/Org/MyStrStream.H" using namespace SHRIMPS; using namespace ATOOLS; Single_Channel_Eikonal:: Single_Channel_Eikonal(const deqmode::code & deq,const int & test) : m_test(test), m_deqmode(deq), p_ff1(NULL), p_ff2(NULL), p_convolution2D(NULL), p_integrator(NULL) { } Single_Channel_Eikonal::~Single_Channel_Eikonal() { if (p_ff1) { delete p_ff1; p_ff1 = NULL; } if (p_ff2) { delete p_ff2; p_ff2 = NULL; } if (p_convolution2D) { delete p_convolution2D; p_convolution2D = NULL; } if (p_integrator) { delete p_integrator; p_integrator = NULL; } } void Single_Channel_Eikonal:: Initialise(Form_Factor * ff1,Form_Factor * ff2, const double & lambda,const double & alpha, const double & Y,const double & ycutoff) { p_ff1 = ff1; m_ff1max = p_ff1->Maximum(); m_ff1bins = p_ff1->Bbins(); m_deltaff1 = m_ff1max/m_ff1bins; m_b1max = p_ff1->Bmax(); p_ff2 = ff2; m_ff2max = p_ff2->Maximum(); m_ff2bins = p_ff2->Bbins(); m_deltaff2 = m_ff2max/m_ff2bins; m_b2max = p_ff2->Bmax(); m_Bmax = Max(m_b1max,m_b2max); m_Bbins = Max(m_ff1bins,m_ff2bins); m_deltaB = m_Bmax/m_Bbins; m_beta2 = p_ff1->Beta0()*p_ff2->Beta0(); m_lambda = lambda; m_alpha = alpha; m_expfactor = 1./2.; m_Y = Y; m_ycutoff = ycutoff; m_yshift = m_Y-m_ycutoff; m_ybins = 20; m_deltay = 2.*m_yshift/m_ybins; if (m_test==10) { m_ycutoff = 0.; m_yshift = m_Y; m_lambda = 0.; } m_accu = 1.e-2; msg_Tracking()< > > (m_ff1bins+1,std::vector > (m_ff2bins+1,std::vector(2,0.))); m_grid2 = std::vector > > (m_ff1bins+1,std::vector > (m_ff2bins+1,std::vector(2,0.))); double value1, value2; for (int i=0;i<=m_ff1bins;i++) { for (int j=0;j<=m_ff2bins;j++) { InitialiseBoundaries(i,j,value1,value2); if (i==0 && j==0) m_ybins = AdjustGrid(i,j,value1,value2); else SolveSystem(i,j,value1,value2,m_ybins); } } msg_Tracking()< comp1, comp2; do { comp1 = m_grid1[i][j]; comp2 = m_grid2[i][j]; ysteps *= 2; m_deltay /= 2.; SolveSystem(i,j,val1,val2,ysteps); } while (CheckAccuracy(i,j,ysteps,comp1,comp2)); return ysteps; } void Single_Channel_Eikonal:: SolveSystem(const int & i,const int & j,double & val1, double & val2,const int & steps) { switch (m_deqmode) { case deqmode::RungeKutta2: return RungeKutta2(i,j,val1,val2,steps); case deqmode::RungeKutta4: return RungeKutta4(i,j,val1,val2,steps); case deqmode::RungeKutta4Transformed: default: return RungeKutta4Transformed(i,j,val1,val2,steps); } } bool Single_Channel_Eikonal:: CheckAccuracy(const int & i,const int & j,const int & ysteps, const std::vector & comp1, const std::vector & comp2) { double diff1(0.), diff2(0.), ave1(0.), ave2(0.); diff1 = diff2 = 0.; for (int step=2;step diff1) { diff1 = dabs(test/x1-1.); } if (dabs(testp/xp1-1.)> diff2) { diff2 = dabs(testp/xp1-1.); } } if (diff1>m_accu || diff2>m_accu) { return true; } return false; } void Single_Channel_Eikonal:: RungeKutta4Transformed(const int & i,const int & j, double & val1,double & val2,const int & steps) { double norm_ik(val1),norm_ki(val2),Omega_ik(norm_ik),Omega_ki(norm_ki); double deltay(m_deltay),y(0.); double barDelta(m_alpha*exp(-m_lambda*m_expfactor*(norm_ik+norm_ki))); double x1_ik(1.),x1_ki(1.),exp1,x2_ik,x2_ki,exp2; double x3_ik,x3_ki,exp3,x4_ik,x4_ki,exp4; double f1_ik,f1_ki,f2_ik,f2_ki,f3_ik,f3_ki,f4_ik,f4_ki; m_grid1[i][j].clear(); m_grid2[i][j].clear(); m_grid1[i][j].push_back(Omega_ik); m_grid2[i][j].push_back(Omega_ki); for (int step=0;step0.) { if (valuevalmax) valmax = value; } m_gridB.push_back(value); B += m_deltaB; } run = false; if (valmin/valmax>m_accu) { m_Bbins *= 2; m_Bmax *= 2.; msg_Tracking()<0.01)<<") " <<"in "< " <<"delta_B = "<PrintErrors(); } double Single_Channel_Eikonal:: IntegrateOutImpactParameters(const double & B,const double & Y) { p_convolution2D->SetB(B); p_convolution2D->SetY(Y); return p_integrator->Integrate(0.,m_b1max,m_accu,1)/m_beta2; } bool Single_Channel_Eikonal:: GeneratePositions(const double & B,double & b1,double & theta1) { double b2, value, y(0.); do { b1 = ran->Get()*m_b1max; theta1 = 2.*M_PI*ran->Get(); b2 = sqrt(B*B+b1*b1-2.*B*b1*cos(theta1)); value = Omega12(b1,b2,y)*Omega21(b1,b2,y); } while (valueGet()); return true; } double Single_Channel_Eikonal:: Omega12(const double & b1,const double & b2,const double & y, const bool & plot) const { if (b1>m_b1max || b1<0. || b2>m_b2max || b2<0. || y>m_yshift || y<-m_yshift) return 0.; double ff1(p_ff1->FourierTransform(b1)), ff2(p_ff2->FourierTransform(b2)); int ff1bin(int((m_ff1max-ff1)/m_deltaff1)); int ff2bin(int((m_ff2max-ff2)/m_deltaff2)); double yactual(y+m_yshift); int ybin(int(yactual/m_deltay)); if (ff1bin<0 || ff1bin>m_ff1bins || ff2bin<0 || ff2bin>m_ff2bins || ybin<0 || ybin>m_ybins) { msg_Error()<<"Error in "< " <<"ff1 = "< ff1bin = "<m_b1max || b1<0. || b2>m_b2max || b2<0. || y>m_yshift || y<-m_yshift) return 0.; double ff1(p_ff1->FourierTransform(b1)), ff2(p_ff2->FourierTransform(b2)); int ff1bin(int((m_ff1max-ff1)/m_deltaff1)); int ff2bin(int((m_ff2max-ff2)/m_deltaff2)); double yactual(m_yshift-y); int ybin(int(yactual/m_deltay)); if (ff1bin<0 || ff1bin>m_ff1bins || ff2bin<0 || ff2bin>m_ff2bins || ybin<0 || ybin>m_ybins) { msg_Error()<<"Error in "< " <<"ff1 = "< ff1bin = "<=m_Bmax) return 0.; int Bbin(int(B/m_deltaB)); return (m_gridB[Bbin]*((Bbin+1)*m_deltaB-B)+ m_gridB[Bbin+1]*(B-Bbin*m_deltaB))/m_deltaB; } double Single_Channel_Eikonal::Convolution2D::operator()(double b1) { p_convolution1D->Setb1(b1); return 2.*b1*p_integrator->Integrate(0.,M_PI,m_accu,1); } double Single_Channel_Eikonal::Convolution2D::Convolution1D:: operator()(double theta1) { double b2((m_b==0.)?m_b1:sqrt(m_b*m_b+m_b1*m_b1-2.*m_b*m_b1*cos(theta1))); double eik12(p_eikonal->Omega12(m_b1,b2,m_y)); double eik21(p_eikonal->Omega21(m_b1,b2,m_y)); double value(eik12*eik21); if (value>p_eikonal->MaxConv()) p_eikonal->SetMaxConv(value); return value; } void Single_Channel_Eikonal::PrintOmega_ik() { double b1(0.), b2(0.); double deltab1(3.), deltab2(3.), deltay(0.1); while (b1<3.) { while (b2<3.) { std::cout<<"Omega_ik for b1 = "< b1 = "<<(i*deltab1)<<" b2 = "<<(j*deltab2)<<"." < b1 = "<<(i*deltab1)<<" b2 = "<<(j*deltab2)<<"." <Lambda2()), kappa(p_ff1->Kappa()); double pref(beta02*Lambda2/(4.*M_PI)*exp(-sqr(b1)*Lambda2/(4.*(1.+kappa)))); double ana, num, y; msg_Tracking()<<" prefactor = "<