File indexing completed on 2025-04-19 09:10:15
0001 #ifndef SHRIMPS_Eikonals_Single_Channel_Eikonal_H
0002 #define SHRIMPS_Eikonals_Single_Channel_Eikonal_H
0003
0004 #include "SHRiMPS/Eikonals/Form_Factors.H"
0005 #include "ATOOLS/Math/Gauss_Integrator.H"
0006 #include "ATOOLS/Math/MathTools.H"
0007 #include "ATOOLS/Org/Message.H"
0008
0009 namespace SHRIMPS {
0010 struct deqmode {
0011 enum code {
0012 RungeKutta4Transformed = -4,
0013 RungeKutta4 = 4,
0014 RungeKutta2 = 2
0015 };
0016 };
0017
0018 class Single_Channel_Eikonal : public ATOOLS::Function_Base {
0019 class Convolution2D : public ATOOLS::Function_Base {
0020 class Convolution1D : public ATOOLS::Function_Base {
0021 private:
0022 Single_Channel_Eikonal * p_eikonal;
0023 double m_b,m_b1,m_Y,m_y;
0024 double m_errmax12, m_errmax21;
0025 int m_test;
0026 public:
0027 Convolution1D(Single_Channel_Eikonal * eikonal,
0028 const double & Y,const int & test) :
0029 p_eikonal(eikonal), m_Y(Y),
0030 m_errmax12(0.), m_errmax21(0.), m_test(test) {}
0031 void SetY(const double & y) { m_y = y; }
0032 void SetB(const double & b) { m_b = b; }
0033 void Setb1(const double & b1) { m_b1 = b1; }
0034 double operator()(double theta1);
0035 void PrintErrors() {
0036 msg_Out()<<"Maximal errors in evaluating product of single terms: "<<std::endl<<" "
0037 <<"delta_max{Omega_12} = "<<m_errmax12<<", "
0038 <<"delta_max{Omega_21} = "<<m_errmax12<<"."<<std::endl;
0039 }
0040 };
0041 private:
0042 Single_Channel_Eikonal * p_eikonal;
0043 double m_b,m_y,m_Y,m_accu;
0044 Convolution1D * p_convolution1D;
0045 ATOOLS::Gauss_Integrator * p_integrator;
0046 int m_test;
0047 public:
0048 Convolution2D(Single_Channel_Eikonal * eikonal,const double & Y,const int & test) :
0049 p_eikonal(eikonal), m_Y(Y), m_accu(1.e-3), m_test(test)
0050 {
0051 p_convolution1D = new Convolution1D(p_eikonal,m_Y,m_test);
0052 p_integrator = new ATOOLS::Gauss_Integrator(p_convolution1D);
0053 }
0054 ~Convolution2D() {
0055 if (p_convolution1D) { delete p_convolution1D; p_convolution1D = NULL; }
0056 if (p_integrator) { delete p_integrator; p_integrator = NULL; }
0057 }
0058 void SetY(const double & y) {
0059 m_y = y;
0060 p_convolution1D->SetY(m_y);
0061 }
0062 void SetB(const double & b) {
0063 m_b = b;
0064 p_convolution1D->SetB(m_b);
0065 }
0066 double operator()(double b1);
0067 void PrintErrors() { p_convolution1D->PrintErrors(); }
0068 };
0069 private:
0070 int m_test;
0071 deqmode::code m_deqmode;
0072 Form_Factor * p_ff1, * p_ff2;
0073
0074 double m_Y, m_ycutoff, m_yshift;
0075 int m_ybins;
0076 double m_deltay;
0077 double m_beta2, m_lambda, m_alpha, m_expfactor;
0078 double m_b1max, m_b2max, m_ff1max, m_ff2max, m_Bmax;
0079 int m_ff1bins, m_ff2bins, m_Bbins;
0080 double m_deltaff1, m_deltaff2, m_deltaB;
0081 double m_accu, m_maxconv;
0082 bool m_inelastic;
0083
0084 std::vector<std::vector<std::vector<double> > > m_grid1, m_grid2;
0085 std::vector<double> m_gridB;
0086
0087
0088 Convolution2D * p_convolution2D;
0089 ATOOLS::Gauss_Integrator * p_integrator;
0090
0091
0092 void ProduceInitialGrids();
0093 void InitialiseBoundaries(const int & i,const int & j,double & val1,double & val2);
0094 void SolveSystem(const int & i,const int & j,double & val1,double & val2,const int & steps);
0095 void RungeKutta2(const int & i, const int & j,double & val1,double & val2,const int & steps);
0096 void RungeKutta4(const int & i, const int & j,double & val1,double & val2,const int & steps);
0097 void RungeKutta4Transformed(const int & i, const int & j,
0098 double & val1,double & val2,const int & steps);
0099 int AdjustGrid(const int & i, const int & j,double & val1,double & val2);
0100 bool CheckAccuracy(const int & i,const int & j,const int & ysteps,
0101 const std::vector<double> & comp1,const std::vector<double> & comp2);
0102
0103 void PrintOmega_ik();
0104 void TestSingleEikonal(const double & b1=0.,const double & b2=0.);
0105 void TestDEQSolvers();
0106 public:
0107 Single_Channel_Eikonal(const deqmode::code & deq=deqmode::RungeKutta4,const int & m_test=0);
0108 ~Single_Channel_Eikonal();
0109
0110 void Initialise(Form_Factor * ff1,Form_Factor * ff2,
0111 const double & lambda,const double & alpha,
0112 const double & Y,const double & ycutoff=0.);
0113
0114 bool GeneratePositions(const double & B,double & b1,double & theta1);
0115
0116 double Omega12(const double & b1,const double & b2,const double & y,const bool & plot=false) const;
0117 double Omega21(const double & b1,const double & b2,const double & y,const bool & plot=false) const;
0118 double operator()(double B);
0119
0120 void ProduceImpactParameterGrid(const double & y);
0121 double IntegrateOutImpactParameters(const double & B,const double & Y=0.);
0122
0123 void SetInelastic(const bool & inelastic) { m_inelastic = inelastic; }
0124 void SetMaxConv(const double & maxconv) { m_maxconv = maxconv; }
0125
0126 const bool & Inelastic() const { return m_inelastic; }
0127 const double & MaxConv() const { return m_maxconv; }
0128 double Prefactor() const { return ATOOLS::sqr(p_ff1->Prefactor()*p_ff2->Prefactor()); }
0129
0130 const double & Beta2() const { return m_beta2; }
0131 const double & Lambda2() const { return p_ff1->Lambda2(); }
0132 const double & Delta() const { return m_alpha; }
0133 const double & Kappa_i() const { return p_ff1->Kappa(); }
0134 const double & Kappa_k() const { return p_ff2->Kappa(); }
0135
0136 const Form_Factor * FF1() const { return p_ff1; }
0137 const Form_Factor * FF2() const { return p_ff1; }
0138 };
0139 }
0140
0141 #endif
