#ifndef SHRIMPS_Eikonals_Single_Channel_Eikonal_H #define SHRIMPS_Eikonals_Single_Channel_Eikonal_H #include "SHRiMPS/Eikonals/Form_Factors.H" #include "ATOOLS/Math/Gauss_Integrator.H" #include "ATOOLS/Math/MathTools.H" #include "ATOOLS/Org/Message.H" namespace SHRIMPS { struct deqmode { enum code { RungeKutta4Transformed = -4, RungeKutta4 = 4, RungeKutta2 = 2 }; }; class Single_Channel_Eikonal : public ATOOLS::Function_Base { class Convolution2D : public ATOOLS::Function_Base { class Convolution1D : public ATOOLS::Function_Base { private: Single_Channel_Eikonal * p_eikonal; double m_b,m_b1,m_Y,m_y; double m_errmax12, m_errmax21; int m_test; public: Convolution1D(Single_Channel_Eikonal * eikonal, const double & Y,const int & test) : p_eikonal(eikonal), m_Y(Y), m_errmax12(0.), m_errmax21(0.), m_test(test) {} void SetY(const double & y) { m_y = y; } void SetB(const double & b) { m_b = b; } void Setb1(const double & b1) { m_b1 = b1; } double operator()(double theta1); void PrintErrors() { msg_Out()<<"Maximal errors in evaluating product of single terms: "<SetY(m_y); } void SetB(const double & b) { m_b = b; p_convolution1D->SetB(m_b); } double operator()(double b1); void PrintErrors() { p_convolution1D->PrintErrors(); } }; private: int m_test; deqmode::code m_deqmode; Form_Factor * p_ff1, * p_ff2; double m_Y, m_ycutoff, m_yshift; int m_ybins; double m_deltay; double m_beta2, m_lambda, m_alpha, m_expfactor; double m_b1max, m_b2max, m_ff1max, m_ff2max, m_Bmax; int m_ff1bins, m_ff2bins, m_Bbins; double m_deltaff1, m_deltaff2, m_deltaB; double m_accu, m_maxconv; bool m_inelastic; std::vector > > m_grid1, m_grid2; std::vector m_gridB; Convolution2D * p_convolution2D; ATOOLS::Gauss_Integrator * p_integrator; void ProduceInitialGrids(); void InitialiseBoundaries(const int & i,const int & j,double & val1,double & val2); void SolveSystem(const int & i,const int & j,double & val1,double & val2,const int & steps); void RungeKutta2(const int & i, const int & j,double & val1,double & val2,const int & steps); void RungeKutta4(const int & i, const int & j,double & val1,double & val2,const int & steps); void RungeKutta4Transformed(const int & i, const int & j, double & val1,double & val2,const int & steps); int AdjustGrid(const int & i, const int & j,double & val1,double & val2); bool CheckAccuracy(const int & i,const int & j,const int & ysteps, const std::vector & comp1,const std::vector & comp2); void PrintOmega_ik(); void TestSingleEikonal(const double & b1=0.,const double & b2=0.); void TestDEQSolvers(); public: Single_Channel_Eikonal(const deqmode::code & deq=deqmode::RungeKutta4,const int & m_test=0); ~Single_Channel_Eikonal(); void Initialise(Form_Factor * ff1,Form_Factor * ff2, const double & lambda,const double & alpha, const double & Y,const double & ycutoff=0.); bool GeneratePositions(const double & B,double & b1,double & theta1); double Omega12(const double & b1,const double & b2,const double & y,const bool & plot=false) const; double Omega21(const double & b1,const double & b2,const double & y,const bool & plot=false) const; double operator()(double B); void ProduceImpactParameterGrid(const double & y); double IntegrateOutImpactParameters(const double & B,const double & Y=0.); void SetInelastic(const bool & inelastic) { m_inelastic = inelastic; } void SetMaxConv(const double & maxconv) { m_maxconv = maxconv; } const bool & Inelastic() const { return m_inelastic; } const double & MaxConv() const { return m_maxconv; } double Prefactor() const { return ATOOLS::sqr(p_ff1->Prefactor()*p_ff2->Prefactor()); } const double & Beta2() const { return m_beta2; } const double & Lambda2() const { return p_ff1->Lambda2(); } const double & Delta() const { return m_alpha; } const double & Kappa_i() const { return p_ff1->Kappa(); } const double & Kappa_k() const { return p_ff2->Kappa(); } const Form_Factor * FF1() const { return p_ff1; } const Form_Factor * FF2() const { return p_ff1; } }; } #endif