[Not intended to be Merged] x-space results

dSigmadptX.cpp

@@ -0,0 +1,153 @@
+/*
+ * =====================================================================================
+ *
+ *       Filename:  main.cpp
+ *
+ *    Description:  Main file that computes either the full hadronic or the
+ * Mellin partonic cross section.
+ *
+ *        Version:  1.0
+ *        Created:  17/02/2021 23:16:58
+ *       Revision:  none
+ *       Compiler:  gcc
+ *
+ *         Author:  Tanjona R. Rabemananjara, Roy Stegeman
+ *   Organization:  N3PDF
+ *
+ * =====================================================================================
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <exception>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#include "higgs-fo/hadronic.h"
+#include "higgs-fo/higgspt.h"
+#include "higgs-fo/higgsptpartonic.h"
+#include "higgs-fo/luminosity.h"
+#include "higgs-fo/params.h"
+#include "higgs-fo/partonic.h"
+#include "include/CombinedRes.h"
+#include "include/SmallptExp.h"
+#include "include/ThresExp.h"
+#include "include/HighEnergyExp.h"
+#include "yaml-cpp/yaml.h"
+
+// Exception for wrong inputs
+struct err_message : public std::exception {
+  const char *what() const throw() { return "Wrong Parameters!!"; }
+};
+
+int main(int argc, char *argv[]) {
+  LHAPDF::setVerbosity(0);
+  YAML::Node node = YAML::LoadFile(argv[1]);
+
+  int inorm = node["inorm"].as<int>();
+  int order = node["order"].as<int>();
+  int _nf = node["nf"].as<int>();
+  int channel = node["channel"].as<int>();
+  int scheme = node["scheme"].as<int>();
+
+  double _mh = node["mh"].as<double>();
+  double _mur = node["mur"].as<double>();
+  double _muf = node["muf"].as<double>();
+
+  double _sroot = node["sroot"].as<double>();
+  long double xmin = node["xmin"].as<long double>();
+  long double xmax = node["xmax"].as<long double>();
+  long double xbin = node["xbin"].as<long double>();
+  long double pt = node["pt"].as<long double>();
+
+  std::string pdfname = node["pdfname"].as<std::string>();
+  std::string filename = node["outfile"].as<std::string>();
+
+  std::string ord_fixod[2] = {"_LO", "_NLO"};
+  std::string par_chanl[5] = {"_gg_channel", "_gq_channel", "_qq_channel",
+                              "_qqb_channel", "_all_channels"};
+  std::string matsch[4] = {"_smallpt_asN.dat", "_threshold_asN.dat",
+                           "_combined_asN.dat", "_highenergy_asN.dat"};
+
+  try {
+    if (order < 0 || order > 1) throw err_message();
+    if (channel < 0 || channel > 5) throw err_message();
+    filename += ord_fixod[order];
+    filename += par_chanl[channel];
+    filename += matsch[scheme];
+  } catch (err_message &err) {
+    std::cout << err.what() << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  // Factors for Born cross-section
+  double factor;
+  double gf = 1.16637e-5;      // Fermi Constant
+  double gevpb = 3.8937966e8;  // GeV to pb
+
+  if (inorm == 1) {
+    std::cout << "ERROR in inorm!" << std::endl;
+    exit(EXIT_FAILURE);  // TODO: complete implementation!
+  } else if (inorm == 0) {
+    factor = gf / 288. / M_PI / sqrt(2.);  // large-top mass limit
+  } else {
+    std::cout << "ERROR in inorm!" << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  // Initialize PDF and extract alpha_s
+  LHAPDF::initPDFSetByName(pdfname);
+  double _as = LHAPDF::alphasPDF(_mur);
+  double _sigma0 = factor * gevpb * std::pow(_as, 2);
+
+  // Define parameters
+  PhysParams physparam;
+  physparam.nc = 3;
+  physparam.nf = _nf;
+  physparam.mh = _mh;
+  physparam.mur = _mur;
+  physparam.muf = _muf;
+  physparam.alphas = _as;
+  physparam.sroot = _sroot;
+  physparam.sigma0 = _sigma0;
+
+  // Init. combined resummation class
+  CombinedRes combres(order, channel, pdfname, &physparam);
+
+  // Construct output fie
+  std::ofstream output_file(filename);
+  output_file << "# PDF set name         : " << pdfname << "\n"
+              << "# Matching scheme      : " << scheme << "\n"
+              << "# Fixed Order          : " << order << "\n"
+              << "# Partonic channel     : " << channel << "\n"
+              << "# Center of M.E. (GeV) : " << _sroot << "\n"
+              << "# Higgs mass (GeV)     : " << _mh << "\n"
+              << "# Renorm. scale (GeV)  : " << _mur << "\n"
+              << "# Fact. scale (GeV)    : " << _muf << "\n";
+  const int space = 16;
+  output_file << "# [x values]" << std::setw(space) << "[dHpt (pb)]"
+              << "\n";
+
+  long double results;
+  long double xx = xmin;
+  while (xx <= xmax) {
+    results = combres.CombinedResExprX(xx, pt, scheme);
+
+    // Generate some output logs & write to output file
+    printf("x=%Le: dHdpt = %Le \n", xx, results);
+    output_file.setf(std::ios_base::scientific);
+    output_file << xx << std::setw(space) << results << "\n";
+    output_file.flush();
+
+    xx += xbin;
+  }
+
+  output_file.close();
+
+  return 0;
+}

include/CombinedRes.h

@@ -32,6 +32,7 @@
 #include "./IntMellin.h"
 #include "./SmallptExp.h"
 #include "./ThresExp.h"
+#include "./ThresXspace.h"
 #include "higgs-fo/partonic.h"
 
 class CombinedRes {
@@ -40,6 +41,7 @@ class CombinedRes {
   virtual ~CombinedRes();
 
   // Attribute that compute the expanded results
+  long double CombinedResExprX(long double x, long double pt, int scheme);
   std::complex<long double> CombinedResExpr(std::complex<long double> N,
                                             long double pt, int scheme);
 
@@ -52,6 +54,7 @@ class CombinedRes {
   // Init. resummation classes
   SmallptExp *SMALLPT;
   ThresExp *THRESHOLD;
+  ThresXspace *xTHRESHOLD;
   HighEnergyExp *HIGHENERGY;
   MellinTrans *MELLIN;
 

include/HighEnergyExp.h

@@ -38,71 +38,21 @@ class HighEnergyExp {
   virtual ~HighEnergyExp();
 
   // Attribute that compute the expanded results
-  /* std::vector<std::complex<long double>> HighEnergyExpExpr( */
-  /*     std::complex<long double> N, long double xp */
-  /* ); */
   std::complex<long double> HighEnergyExpExpr(std::complex<long double> N,
                                               long double pt);
+  long double HighEnergyExpExprX(long double x, long double pt);
 
-  // Matching Coefficients
-  long double Hth1gggH(long double xp);
-  long double Hth1gqqH(long double xp);
-  long double Hth1qqgH(long double xp);
+  // Coefficients in x space
+  long double C1ggx(long double x, long double xp);
+  long double C2ggx(long double x, long double xp);
 
+  // Coefficients in N space
   std::complex<long double> C1gg(std::complex<long double> NN, long double xp);
   std::complex<long double> C2gg(std::complex<long double> NN, long double xp);
-  // LO pt-distribution
-  std::complex<long double> LOgggH(std::complex<long double> NN,
-                                   long double xp);
-  std::complex<long double> LOgqqH(std::complex<long double> NN,
-                                   long double xp);
-  std::complex<long double> LOqqgH(std::complex<long double> NN,
-                                   long double xp);
-
-  // Matching function G)
-  long double GOgggH(long double);
-  long double GOgqqH(long double);
-  long double GOqqgH(long double);
 
  private:
   int NC, NF, CA, ORD, CHANNEL;
   long double LF, LR, LQ;
   long double MH2, MUF2, MUR2, SROOT;
   long double CF, aass, SIGMA0;
-
-  // Beta functions
-  long double Beta0;
-  long double Beta1;
-  long double Beta2;
-
-  // Cusp Anomalous Dimensions
-  long double Ath1g;
-  long double Ath2g;
-  long double Ath1q;
-  long double Ath2q;
-  long double Bth1g;
-  long double Bth1q;
-
-  // Sigma functions
-  // gg->g
-  long double Sigma22ggg(long double xp);
-  long double Sigma21ggg(long double xp);
-  long double Sigma20ggg(long double xp);
-  // gq->g
-  long double Sigma22gqg(long double xp);
-  long double Sigma21gqg(long double xp);
-  long double Sigma20gqg(long double xp);
-  // qq->g
-  long double Sigma22qqg(long double xp);
-  long double Sigma21qqg(long double xp);
-  long double Sigma20qqg(long double xp);
-
-  // Zeta functions
-  const long double zeta2 = gsl_sf_zeta_int(2);
-  const long double zeta3 = gsl_sf_zeta_int(3);
-  const long double zeta4 = gsl_sf_zeta_int(4);
-
-  // EulerGamma
-  // TODO: Move this to global
-  const long double EulerGamma = 0.5772156649;
 };

meson.build

@@ -48,6 +48,15 @@ higgsexen = executable(
   c_args : extra_args
 )
 
+# create executable for cross section as a function of x
+higgsexex = executable(
+  'higgs-x', 'dSigmadptX.cpp',
+  include_directories : inc,
+  dependencies : [glib_dep, higgsfo_lib, yaml_lib],
+  link_with : higgslib,
+  c_args : extra_args
+)
+
 pkg_mod = import('pkgconfig')
 pkg_mod.generate(
   libraries : higgslib,

src/CombinedRes.cpp

@@ -34,6 +34,7 @@ CombinedRes::CombinedRes(int order, int channel, std::string pdfname,
   EXACT_ORD = order == 0 ? 0 : order - 1;
 
   SMALLPT = new SmallptExp(order, channel, params);
+  xTHRESHOLD = new ThresXspace(order, channel, params);
   THRESHOLD = new ThresExp(order, channel, params);
   HIGHENERGY = new HighEnergyExp(order, channel, params);
   MELLIN = new MellinTrans(order, channel, pdfname, params);
@@ -95,13 +96,52 @@ std::complex<long double> CombinedRes::CombinedResExpr(
     return HighEnergyMellin;
   } else {
     std::complex<long double> SptMellin = SMALLPT->SmallptExpExpr(N, pt);
-    /* std::complex<long double> ThresMellin = THRESHOLD->ThresExpExpr(N, pt);
-     */
-    std::complex<long double> xThresMellin = MELLIN->xSpaceThres(N, pt);
+    std::complex<long double> ThresMellin = THRESHOLD->ThresExpExpr(N, pt);
+    // std::complex<long double> xThresMellin = MELLIN->xSpaceThres(N, pt);
 
     /* std::complex<long double> ExactMellinCmpx(ExactMellin[0], 0.); */
     mres = (1. - Matching(N, pt, scheme)) * SptMellin +
-           Matching(N, pt, scheme) * xThresMellin;
+           Matching(N, pt, scheme) * ThresMellin;
+    /* Matching(N, pt, scheme) * ThresMellin; */
+
+    /* return ExactMellinCmpx + mres; */
+    return mres;
+  }
+}
+
+long double CombinedRes::CombinedResExprX(long double x, long double pt,
+                                          int scheme) {
+  /* double pp = static_cast<double>(pt); */
+  // take only real part. Does not work for complex
+  /* double nn = static_cast<double>(N.real()); */
+  // std::complex<long double> mres;
+  /* std::vector<double> ResultsMellin; */
+  /* std::vector<double> zero(2, 0.0); */
+
+  // Compute exact FO from HpT-MON
+  /* if (ORD == 0) { */
+  /*   ResultsMellin = zero; */
+  /* } else { */
+  /*   ResultsMellin = MELLINPARTONIC->partonichiggsdpt(pp, nn); */
+  /* } */
+  /* std::vector<long double> ExactMellin(ResultsMellin.begin(), */
+  /*                                      ResultsMellin.end()); */
+
+  // Compute approximation from resummations
+  if (scheme == 3)  // High Energy
+  {
+    /* std::complex<long double> ExactMellinCmpx(ExactMellin[0], 0.); */
+    long double HighEnergyMellin = HIGHENERGY->HighEnergyExpExprX(x, pt);
+    /* return ExactMellinCmpx + HighEnergyMellin; */
+    return HighEnergyMellin;
+  } else {
+    // std::complex<long double> SptMellin = SMALLPT->SmallptExpExpr(N, pt);
+    // std::complex<long double> ThresMellin = THRESHOLD->ThresExpExpr(N, pt);
+    std::complex<long double> xThresMellin =
+        xTHRESHOLD->ThresXspaceExpr(x, 1, pt);
+
+    /* std::complex<long double> ExactMellinCmpx(ExactMellin[0], 0.); */
+    long double mres = static_cast<long double>(xThresMellin.real());
     /* Matching(N, pt, scheme) * ThresMellin; */
 
     /* return ExactMellinCmpx + mres; */