Support 3S states (only transverse polarization at the moment)

Add wave functions for upsilon(2S) and upsilon(3S)

gauss-boosted-upsilon-2s.dat

@@ -0,0 +1,29 @@
+# Parameters for upsilon and \gamma^* wave function overlap
+# 
+# Model: Boosted gaussian 
+#
+# Parameters from arXiv:0905.0102 (only transverse)
+#
+# Syntax of this file is the following:
+# Lines starting with # are comments
+# Parameters are given as a pair key:value
+# Values are in units of GeV^n and possible keys are:
+# e_f, m_f, M_V, N_T, N_L, R, S alpha and del
+#
+#
+#e_f: Effective charge: 
+e_f:0.333333
+#N_T, N_L, R_T, R_L: Normalization for transversial/lognitudial components
+#In 0905.0102 normalization factor 0.481 for the scalar part is given
+N_T:0.624
+N_L:0
+#sqrt(0.831)
+R:0.91159
+m_f:4.20
+M_V:10.023
+# delta: 0 or 1: one part of the wave functio is multiplied by delta
+del:1
+# State S=n, J/Psi is S=1
+S:2
+# alpha needed for higher S states, no effect for S=1
+alpha:-0.555

gauss-boosted-upsilon-3s.dat

@@ -0,0 +1,30 @@
+# Parameters for upsilon and \gamma^* wave function overlap
+# 
+# Model: Boosted gaussian 
+#
+# Parameters from arXiv:0905.0102 (only transverse)
+#
+# Syntax of this file is the following:
+# Lines starting with # are comments
+# Parameters are given as a pair key:value
+# Values are in units of GeV^n and possible keys are:
+# e_f, m_f, M_V, N_T, N_L, R, S alpha,alpha2 and del
+#
+#
+#e_f: Effective charge: 
+e_f:0.333333
+#N_T, N_L, R_T, R_L: Normalization for transversial/lognitudial components
+#In 0905.0102 normalization factor 0.481 for the scalar part is given
+N_T:0.668
+N_L:0
+#sqrt(1.028)
+R:1.013903348
+m_f:4.20
+M_V:10.3552
+# delta: 0 or 1: one part of the wave functio is multiplied by delta
+del:1
+# State S=n, J/Psi is S=1
+S:3
+# alpha needed for higher S states, no effect for S=1
+alpha:-1.219
+alpha2:0.217

src/gauss_boost.cpp

@@ -62,7 +62,8 @@ BoostedGauss::BoostedGauss(std::string file)
         return;
     }
     std::string line;
-
+    alpha2=0;
+
     while(!f.eof() )
     {
         std::getline(f, line);
@@ -84,8 +85,11 @@ BoostedGauss::BoostedGauss(std::string file)
             delta=StrToInt(line.substr(4));
         if (line.substr(0,2)=="S:")
             S=StrToInt(line.substr(2));
-        if (line.substr(0,5)=="alpha")
+        if (line.substr(0,6)=="alpha2")
+            alpha2=StrToReal(line.substr(7));      
+        else if (line.substr(0,5)=="alpha") // else if, as this would be true for alpha2 as well
             alpha=StrToReal(line.substr(6));
+
     }
     f.close();
 
@@ -214,43 +218,84 @@ REAL BoostedGauss::PsiSqr_L_intz(REAL Qsqr, REAL r)
  */
 REAL BoostedGauss::Psi_T(REAL r, REAL z)
 {
-    return N_T*z*(1.0-z)*exp(-SQR(m_f)*SQR(R)/(8.0*z*(1.0-z))
-        - 2.0*z*(1-z)*SQR(r)/SQR(R) + SQR(m_f)*SQR(R)/2.0 )
-        * (1.0 + alpha*
-                ( 2.0 - SQR(m_f*R) + SQR(m_f*R)/(4.0*z*(1.0-z))
-                    - 4.0*z*(1.0-z)*SQR(r)/SQR(R) ) );
-}
+    return N_T*Psi(r,z);
+}             
 
 REAL BoostedGauss::Psi_L(REAL r, REAL z)
 {
-        return N_L*z*(1.0-z)*exp(-SQR(m_f)*SQR(R)/(8.0*z*(1.0-z))
-        - 2.0*z*(1.0-z)*SQR(r)/SQR(R) + SQR(m_f)*SQR(R)/2.0 )
-            * (1.0 + alpha*
-                ( 2.0 - SQR(m_f*R) + SQR(m_f*R)/(4.0*z*(1.0-z))
-                    - 4.0*z*(1.0-z)*SQR(r)/SQR(R) ) );
+    return N_L*Psi(r,z);
+}
+
+// Scalar part of the wave function without normalizatin factor
+REAL BoostedGauss::Psi(REAL r, REAL z)
+{
+    double res = z*(1.0-z)*exp(-SQR(m_f)*SQR(R)/(8.0*z*(1.0-z))
+        - 2.0*z*(1-z)*SQR(r)/SQR(R) + SQR(m_f)*SQR(R)/2.0 );
+    if (S==1)
+        return res;
+
+    double g = 2.0 - SQR(m_f*R) 
+        + SQR(m_f*R)/(4.*z*(1.-z)) - 4.*z*(1.-z)*SQR(r)/SQR(R);
+
+    if (S ==2)
+        return res *(1 + alpha * g);
+    else if (S==3)
+        return res*(1.+alpha*g+alpha2*(SQR(g) + 4.*(1.-4.*z*(1.-z)*SQR(r)/SQR(R))));
+    else
+    {
+        std::cerr << "Unknown state " << S << " in BoostedGauss::Psi_T" << std::endl;
+        return 0;
+    }
+}
+
+// \partial_r Psi(r,z)
+REAL BoostedGauss::Psi_DR(REAL r, REAL z)
+{
+    if (S==1)
+        return -4.0*z*(1.0-z)*r/SQR(R)*Psi(r,z); ;
+
+    // We need the "S=1" wave function denoted by G_nS in 0905.0102
+    double Gn = z*(1.0-z)*exp(-SQR(m_f)*SQR(R)/(8.0*z*(1.0-z))
+        - 2.0*z*(1-z)*SQR(r)/SQR(R) + SQR(m_f)*SQR(R)/2.0 ); 
+    double Gn_der = -4.0*z*(1.0-z)*r/SQR(R)*Gn;
+
+    if (S>1)
+    {
+        double der1 = -8*r*(1.-z)*z*alpha*Gn/SQR(R) 
+            + (2.-SQR(m_f*R)+SQR(m_f*R)/(4.*z*(1.-z))-4.*SQR(r)*z*(1.-z)/SQR(R))*alpha*Gn_der;
+        if (S==2)
+            return Gn_der+der1;
+
+        double der2 = (-((32*r*(1-z)*z)/SQR(R)) 
+            - (16*r*(1-z)*z*(2-SQR(m_f)*SQR(R) + (SQR(m_f)*SQR(R))/(4*(1-z)*z) - (4*r*r*(1-z)*z)/SQR(R)))/SQR(R)) 
+                * alpha2 * Gn;
+        double der2_2 = (4 * (1 - (4 * SQR(r) * (1 - z) * z) / SQR(R)) + (2 - SQR(m_f) * SQR(R) + (SQR(m_f) * SQR(R)) / (4 * (1 - z) * z) - (4 * r * r * (1 - z) * z) / SQR(R)) * (2 - SQR(m_f) * SQR(R) + (SQR(m_f) * SQR(R)) / (4 * (1 - z) * z) - (4 * r * r * (1 - z) * z) / SQR(R))) * alpha2 * Gn_der;
+        return Gn_der + der1 + der2 + der2_2;
+    }
+    std::cerr << "Unknown state " << S << " in BoostedGauss::Psi_DR" << std::endl;
+    return 0;
 }
 
 // \partial_r Psi_T(r,z)
 REAL BoostedGauss::Psi_T_DR(REAL r, REAL z)
 {
-    return -4.0*z*(1.0-z)*r/SQR(R)*Psi_T(r,z)
-        -8.0*alpha*N_T*SQR(z*(1.0-z)) * r/SQR(R)
-            * std::exp(-SQR(m_f)*SQR(R)/(8.0*z*(1.0-z))
-                - 2.0*z*(1-z)*SQR(r)/SQR(R) + SQR(m_f)*SQR(R)/2.0 )  ;
+    return N_T*Psi_DR(r,z);
 }
 
 // \partial_r PSI_L(r,z)
 REAL BoostedGauss::Psi_L_DR(REAL r, REAL z)
 {
-    return -4.0*z*(1-z)*r/SQR(R)*Psi_L(r,z)
-        -8.0*alpha*N_L*SQR(z*(1.0-z)) * r/SQR(R)
-            * std::exp(-SQR(m_f)*SQR(R)/(8.0*z*(1.0-z))
-                - 2.0*z*(1.0-z)*SQR(r)/SQR(R) + SQR(m_f)*SQR(R)/2.0 )  ;
+    return N_L*Psi_DR(r,z);
 }
 
 // \partial^2_r PSI_L(r,z)
 REAL BoostedGauss::Psi_L_D2R(REAL r, REAL z)
 {
+    if (S>2)
+    {
+        std::cerr << "BoostedGauss::Psi_L_D2R not implemented for S>2" << std::endl;
+        exit(1);
+    }
     return -4.0*z*(1.0-z)/SQR(R)*Psi_L(r,z) 
         + SQR(4.0*z*(1.0-z)*r/SQR(R))*Psi_L(r,z)
             - 8.0*alpha*N_T*SQR(z*(1.0-z))/std::pow(R,4.0)
@@ -270,7 +315,8 @@ std::string BoostedGauss::GetParamString()
     str << "e_f = "
     << e_f << ", m_f = " << m_f << ", M_V = " << M_V << ", N_T = "
     << N_T << ", N_L = " << N_L << ", R = " << R 
-    << ", delta = " << delta << " " << " S = " << S << ", alpha = " << alpha ;
+    << ", delta = " << delta << " " << " S = " << S << ", alpha = " << alpha
+    << ", alpha_2=" << alpha2 ;
     return str.str();
 }
 

src/gauss_boost.hpp

@@ -51,12 +51,16 @@ class BoostedGauss : public WaveFunction {
 
 
     private:
+        REAL Psi(REAL r, REAL z);
+        REAL Psi_DR(REAL r, REAL z);
+
         // Parameters
         REAL e_f;
         REAL N_T, N_L; // Constants for Psi_T and Psi_L
         REAL R;
         REAL m_f,M_V;
         REAL alpha; // needed for higher exited states like 2s
+        REAL alpha2; // needed for higher 3s
         int S;      // S=1 is J/Psi etc.
         int delta;  // There are two different longitudial wave functions,
                     // delta=0 or delta=1

src/test.cpp

@@ -107,6 +107,9 @@ TEST(structure_function_ipsat)
 
 }
 
+
+
+// Wave function normalizations
 struct inthelperf_boosted_gauss
 {
     double z;
@@ -141,29 +144,88 @@ double integrand_z(double z, void* p) {
     return result;
 }
 
-TEST(gauss_boosted_normalization)
+TEST(gauss_boosted_normalization_jpsi)
 {
     BoostedGauss wf("gauss-boosted_mzsat.dat");
     inthelperf_boosted_gauss p; p.wf=&wf;
 
+    double result, error;
 
-    double result;
-    double error;
+    gsl_function F;
+    F.function = integrand_z;
+    F.params = &p;
+
+    gsl_integration_workspace* w = gsl_integration_workspace_alloc(1000);
+    double eps=1e-5;
+    gsl_integration_qags(&F, eps, 1-eps, 0, 1e-6, 1000, w, &result, &error);
+    ASSERT_ALMOST_EQUAL(result, 1, 1e-3);
+
+    gsl_integration_workspace_free(w); 
+}
+
+TEST(gauss_boosted_normalization_upsilon_1s)
+{
+    BoostedGauss wf("gauss-boosted-upsilon.dat");
+    inthelperf_boosted_gauss p; p.wf=&wf;
+
+    double result, error;
 
     gsl_function F;
     F.function = integrand_z;
     F.params = &p;
 
     gsl_integration_workspace* w = gsl_integration_workspace_alloc(1000);
-    double eps=1e-4;
+    double eps=1e-5;
     gsl_integration_qags(&F, eps, 1-eps, 0, 1e-6, 1000, w, &result, &error);
+    ASSERT_ALMOST_EQUAL(result, 1, 1e-2);
 
-    gsl_integration_workspace_free(w);
+    gsl_integration_workspace_free(w); 
+}
 
-    ASSERT_ALMOST_EQUAL(result, 1, 1e-3);
+TEST(gauss_boosted_normalization_upsilon_2s)
+{
+    BoostedGauss wf("gauss-boosted-upsilon-2s.dat");
+    inthelperf_boosted_gauss p; p.wf=&wf;
 
+    ASSERT_ALMOST_EQUAL(wf.Psi_T(1,0.3),-0.0215006,1e-5);
+    ASSERT_ALMOST_EQUAL(wf.Psi_T_DR(1,0.5),0.056788,1e-5);
+
+    double result, error;
+
+    gsl_function F;
+    F.function = integrand_z;
+    F.params = &p;
+
+    gsl_integration_workspace* w = gsl_integration_workspace_alloc(1000);
+    double eps=1e-5;
+    gsl_integration_qags(&F, eps, 1-eps, 0, 1e-6, 1000, w, &result, &error);
+    ASSERT_ALMOST_EQUAL(result, 1, 1e-2);
+
+    gsl_integration_workspace_free(w); 
 }
 
+TEST(gauss_boosted_normalization_upsilon_3s)
+{
+    BoostedGauss wf("gauss-boosted-upsilon-3s.dat");
+    inthelperf_boosted_gauss p; p.wf=&wf;
+
+    // Values computed for the wave functions presented in 0905.0102
+     ASSERT_ALMOST_EQUAL(wf.Psi_T(2,0.3),-0.012636,1e-5);
+     ASSERT_ALMOST_EQUAL(wf.Psi_T_DR(2,0.3),0.0111469,1e-5);
+
+    double result, error;
+
+    gsl_function F;
+    F.function = integrand_z;
+    F.params = &p;
+
+    gsl_integration_workspace* w = gsl_integration_workspace_alloc(1000);
+    double eps=1e-5;
+    gsl_integration_qags(&F, eps, 1-eps, 0, 1e-6, 1000, w, &result, &error);
+    ASSERT_ALMOST_EQUAL(result, 1, 1e-2);
+
+    gsl_integration_workspace_free(w); 
+}
 
 // DO NOT REMOVE
 // Generates a main() function that runs all of your tests.