Implementing more function in doubleobject.cc

doc/src/SIDISTMD.tex

@@ -49,21 +49,27 @@ \section{Structure of the observable}
 of semi-inclusive deep-inelastic scattering (SIDIS) multiplicities
 under the assumption that the (negative) virtuality of the $Q^2$ of
 the exchanged vector boson is much smaller than the $Z$ mass. This
-allows us to neglect weak contributions and write the multiplicity
-as:
+allows us to neglect weak contributions and write the cross section
+in TMD factorisation as:
 \begin{equation}
-\frac{d\sigma}{dxdQ^2 dz d p_T} = \frac{2p_T\pi \alpha^2
-  }{z^2 x Q^4}\left[1+(1-y)^2\right] H(Q,\mu) \sum_q e_q^2
-  \int_0^\infty db \,b J_0\left(\frac{bp_T}{z}\right)\overline{F}_q(x,b;\mu,\zeta) \overline{D}_{q}(z,b;\mu,\zeta)\,,
+  \frac{d\sigma}{dxdQdz d q_T} = \frac{4\pi \alpha^2q_T}{z x Q^3}Y_+ H(Q,\mu) \sum_q e_q^2
+  \int_0^\infty db \,b J_0\left(bq_T\right)\overline{F}_q(x,b;\mu,\zeta_1) \overline{D}_{q}(z,b;\mu,\zeta_2)\,,
 \end{equation}
-where:
+with:
+\begin{equation}
+Y_+=1+(1-y)^2=1+\left(1-\frac{Q^2}{xs}\right)^2\,,
+\end{equation}
+and where:
 \begin{equation}
 \overline{F}_i(x,b;\mu,\zeta) =xF_i(x,b;\mu,\zeta) = R_q(\mu_0,\zeta_0\rightarrow \mu,\zeta;b) \sum_{j}\int_x^1dy\,\mathcal{C}_{ij}(y;\mu_0,\zeta_0)\left[\frac{x}{y}f_j\left(\frac{x}{y},\mu_0\right)\right]\,,
 \end{equation}
 and:
 \begin{equation}
-\overline{D}_{i}(z,b;\mu,\zeta) =z^3D_{i}(z,b;\mu,\zeta) = R_q(\mu_0,\zeta_0\rightarrow \mu,\zeta;b) \sum_{j}\int_z^1dy\,\left[y^2\mathbb{C}_{ij}(y;\mu_0,\zeta_0)\right]\left[\frac{z}{y}d_j\left(\frac{z}{y},\mu_0\right)\right]\,.
+\overline{D}_{i}(z,b;\mu,\zeta) =z^3D_{i}(z,b;\mu,\zeta) = R_q(\mu_0,\zeta_0\rightarrow \mu,\zeta;b) \sum_{j}\int_z^1dy\,\left[y^2\mathbb{C}_{ij}(y;\mu_0,\zeta_0)\right]\left[\frac{z}{y}d_j\left(\frac{z}{y},\mu_0\right)\right]\,,
 \end{equation}
+are the distributions directly computed by APFEL++.
+
+
 
 \begin{thebibliography}{alp}
 

inc/apfel/doubleobject.h

@@ -65,15 +65,25 @@ namespace apfel
     std::vector<term<T>> GetTerms() const { return _terms; };
 
     /**
-     * @brief Funtion that evaluates the double distribution.
+     * @brief Function that evaluates the double distribution.
      * @param x: value of the first variable
      * @param z: value of the second variable
      * @return The value of the double distribution in (x, z)
      */
     double Evaluate(double const& x, double const& z) const;
 
     /**
-     * @brief Funtion that evaluates the derivative of the double
+     * @brief Function that evaluates the double distribution in one of
+     * the two variables leaving the other undetermined.
+     * @param iv: variable index. It can be either 1 or 2.
+     * @param y: value of the iv-th variable
+     * @return The value of the double object in y for any other value
+     * of the other variable
+     */
+    T Evaluate(int const& iv, double const& y) const;
+
+    /**
+     * @brief Function that evaluates the derivative of the double
      * distribution.
      * @param x: value of the first variable
      * @param z: value of the second variable
@@ -83,7 +93,19 @@ namespace apfel
     double Derive(double const& x, double const& z) const;
 
     /**
-     * @brief Funtion that evaluates the integral of the double
+     * @brief Function that evaluates the derivative of a double
+     * distribution in one of the two variables leaving the other
+     * undetermined.
+     * @param iv: variable index. It can be either 1 or 2.
+     * @param y: value of the iv-th variable in which the derivative
+     * has to be computed
+     * @return The value of the derivative of the double object in y
+     * for any other value of the other variable
+     */
+    T Derive(int const& iv, double const& y) const;
+
+    /**
+     * @brief Function that evaluates the integral of the double
      * distribution.
      * @param xl: value of the lower bound of the of the first variable
      * @param xu: value of the upper bound of the of the first variable
@@ -93,6 +115,18 @@ namespace apfel
      */
     double Integrate(double const& xl, double const& xu, double const& zl, double const& zu) const;
 
+    /**
+     * @brief Function that evaluates the integral of the double
+     * distribution in one of the two variables leaving the other
+     * undetermined.
+     * @param iv: variable index. It can be either 1 or 2.
+     * @param yl: value of the lower bound of the of the iv-th variable
+     * @param yu: value of the upper bound of the of the iv-th variable
+     * @return The value of the integral of the double distribution in
+     * [yl: yu] for any other value of the other variable
+     */
+    T Integrate(int const& iv, double const& yl, double const& yu) const;
+
     /**
      * @brief This function multiplies the distributions single terms
      * of the DoubleObject by a respective function.

src/kernel/doubleobject.cc

@@ -133,6 +133,29 @@ namespace apfel
     return result;
   }
 
+  //_________________________________________________________________________________
+  template<>
+  Distribution DoubleObject<Distribution>::Evaluate(int const& iv, double const& y) const
+  {
+    if (iv < 1 || iv >2)
+      throw std::runtime_error(error("Evaluate", "Function index out of range: it can be either 1 or 2."));
+
+    if (iv == 1)
+      {
+        Distribution result = _terms[0].coefficient * _terms[0].object1.Evaluate(y) * _terms[0].object2;
+        for (int i = 1; i < (int) _terms.size(); i++)
+          result += _terms[i].coefficient * _terms[i].object1.Evaluate(y) * _terms[i].object2;
+        return result;
+      }
+    else
+      {
+        Distribution result = _terms[0].coefficient * _terms[0].object2.Evaluate(y) * _terms[0].object1;
+        for (int i = 1; i < (int) _terms.size(); i++)
+          result += _terms[i].coefficient * _terms[i].object2.Evaluate(y) * _terms[i].object1;
+        return result;
+      }
+  }
+
   //_________________________________________________________________________________
   template<>
   double DoubleObject<Distribution>::Derive(double const& x, double const& z) const
@@ -146,6 +169,29 @@ namespace apfel
     return result;
   }
 
+  //_________________________________________________________________________________
+  template<>
+  Distribution DoubleObject<Distribution>::Derive(int const& iv, double const& y) const
+  {
+    if (iv < 1 || iv >2)
+      throw std::runtime_error(error("Derive", "Function index out of range: it can be either 1 or 2."));
+
+    if (iv == 1)
+      {
+        Distribution result = _terms[0].coefficient * _terms[0].object1.Derive(y) * _terms[0].object2;
+        for (int i = 1; i < (int) _terms.size(); i++)
+          result += _terms[i].coefficient * _terms[i].object1.Derive(y) * _terms[i].object2;
+        return result;
+      }
+    else
+      {
+        Distribution result = _terms[0].coefficient * _terms[0].object2.Derive(y) * _terms[0].object1;
+        for (int i = 1; i < (int) _terms.size(); i++)
+          result += _terms[i].coefficient * _terms[i].object2.Derive(y) * _terms[i].object1;
+        return result;
+      }
+  }
+
   //_________________________________________________________________________________
   template<>
   double DoubleObject<Distribution>::Integrate(double const& xl, double const& xu, double const& zl, double const& zu) const
@@ -159,6 +205,29 @@ namespace apfel
     return result;
   }
 
+  //_________________________________________________________________________________
+  template<>
+  Distribution DoubleObject<Distribution>::Integrate(int const& iv, double const& yl, double const& yu) const
+  {
+    if (iv < 1 || iv >2)
+      throw std::runtime_error(error("Integrate", "Function index out of range: it can be either 1 or 2."));
+
+    if (iv == 1)
+      {
+        Distribution result = _terms[0].coefficient * _terms[0].object1.Integrate(yl, yu) * _terms[0].object2;
+        for (int i = 1; i < (int) _terms.size(); i++)
+          result += _terms[i].coefficient * _terms[i].object1.Integrate(yl, yu) * _terms[i].object2;
+        return result;
+      }
+    else
+      {
+        Distribution result = _terms[0].coefficient * _terms[0].object2.Integrate(yl, yu) * _terms[0].object1;
+        for (int i = 1; i < (int) _terms.size(); i++)
+          result += _terms[i].coefficient * _terms[i].object2.Integrate(yl, yu) * _terms[i].object1;
+        return result;
+      }
+  }
+
   //_________________________________________________________________________
   template<>
   DoubleObject<Operator>& DoubleObject<Operator>::operator *= (DoubleObject<Operator> const& o)