diff --git a/MC/config/PWGGAJE/external/generator/generator_pythia8_powheg.C b/MC/config/PWGGAJE/external/generator/generator_pythia8_powheg.C
new file mode 100644
index 000000000..989f33971
--- /dev/null
+++ b/MC/config/PWGGAJE/external/generator/generator_pythia8_powheg.C
@@ -0,0 +1,105 @@
+R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT)
+///#include "FairGenerator.h"
+//#include "Generators/GeneratorPythia8.h"
+#include "Pythia8/Pythia.h"
+#include "Generators/GeneratorPythia8Param.h"
+#include "CommonUtils/FileSystemUtils.h"
+// Pythia8 generator with POWHEG
+//
+// Author: Marco Giacalone (marco.giacalone@cern.ch)
+
+// o2-sim-dpl-eventgen --nEvents 10 --generator external\
+   --configKeyValues "GeneratorExternal.fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/\
+   generator/generator_pythia8_powheg.C;GeneratorExternal.funcName=\
+   getGeneratorJEPythia8POWHEG(\"powheg.input\",\"${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_powheg.cfg\")"
+// or with iniFile
+// o2-sim -g external --noGeant -n 2 -j 8 --configFile $O2DPG_MC_CONFIG_ROOT/MC/config/PWGGAJE/ini/GeneratorPythia8POWHEG.ini
+
+namespace o2
+{
+namespace eventgen
+{
+
+using namespace Pythia8;
+
+// Pythia8 generator using POWHEG data that are generated during the initialization
+// of the external generator. The POWHEG configuration file is copied to the current
+// directory with the right name and the POWHEG events are generated using the pwhg_main_hvq executable.
+class GeneratorJEPythia8POWHEG : public o2::eventgen::GeneratorPythia8
+{
+public:
+  /// default constructor
+  GeneratorJEPythia8POWHEG(std::string confpath = "pwgpath")
+  {
+    // Check if file exist and is not empty
+    if (std::filesystem::exists(confpath) && std::filesystem::file_size(confpath) > 0) {
+      // Copy the file to the current directory
+      ifstream src(confpath);
+      ofstream dst("powheg.input");
+      gRandom->SetSeed(0);
+      int seed = gRandom->Integer(900000000);
+      bool isseed = false;
+      bool isnumevts = false;
+      std::string line;
+      while (std::getline(src, line)) {
+        if (line.find("iseed") != std::string::npos)
+        {
+          // Set the seed to the random number
+          line = "iseed " + std::to_string(seed);
+          isseed = true;
+        }
+        if (line.find("numevts") != std::string::npos)
+        {
+          // Set the number of events to the number of events defined in the configuration
+          line = "numevts " + std::to_string(mSimConfig.getNEvents());
+          // replace it in the file
+          isnumevts = true;
+        }
+        dst << line << std::endl;
+      }
+      if (!isseed) {
+        dst << "iseed " << seed << std::endl;
+      }
+      if (!isnumevts) {
+        dst << "numevts " << mSimConfig.getNEvents() << std::endl;
+      }
+      src.close();
+      dst.close();
+    } else {
+      LOG(fatal) << "POWHEG configuration file not found or empty" << std::endl;
+      exit(1);
+    }
+    // Generate the POWHEG events
+    std::string cmd = "pwhg_main_hvq";
+    system(cmd.c_str());
+  };
+
+private:
+  o2::conf::SimConfig mSimConfig = o2::conf::SimConfig::Instance(); // local sim config object
+};
+
+} // namespace eventgen
+} // namespace o2
+
+/** generator instance and settings **/
+
+FairGenerator *getGeneratorJEPythia8POWHEG(std::string powhegconf = "pwgpath", std::string pythia8conf = "")
+{
+  using namespace o2::eventgen;
+  // Expand paths for the POWHEG configuration file
+  powhegconf = o2::utils::expandShellVarsInFileName(powhegconf);
+  LOG(info) << "Using POWHEG configuration file: " << powhegconf;
+  auto myGen = new GeneratorJEPythia8POWHEG(powhegconf);
+  if(GeneratorPythia8Param::Instance().config.empty() && pythia8conf.empty()) {
+    LOG(fatal) << "No configuration provided for Pythia8";
+  }
+  else if (!pythia8conf.empty())
+  {
+    // Force the configuration for Pythia8 in case it is provided.
+    // Useful for setting up the generator in the hybrid configuration
+    // making it more versatile and not relying entirely on the parameters provided
+    // by ini file or static parameters
+    myGen->setConfig(pythia8conf);
+  }
+  return myGen;
+}
\ No newline at end of file
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_beauty_F05.input b/MC/config/PWGGAJE/external/powheg/powheg_beauty_F05.input
new file mode 100644
index 000000000..7418cbf9a
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_beauty_F05.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 4.75      ! mass of heavy quark in GeV
+facscfact 0.5   ! factorization scale factor: mufact=muref*facscfact
+renscfact 1     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 10000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_beauty_F05_R05.input b/MC/config/PWGGAJE/external/powheg/powheg_beauty_F05_R05.input
new file mode 100644
index 000000000..4dc7895b2
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_beauty_F05_R05.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 4.75      ! mass of heavy quark in GeV
+facscfact 0.5   ! factorization scale factor: mufact=muref*facscfact
+renscfact 0.5   ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 10000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_beauty_F2.input b/MC/config/PWGGAJE/external/powheg/powheg_beauty_F2.input
new file mode 100644
index 000000000..753049c7d
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_beauty_F2.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 4.75      ! mass of heavy quark in GeV
+facscfact 2     ! factorization scale factor: mufact=muref*facscfact
+renscfact 1     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 10000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_beauty_F2_R2.input b/MC/config/PWGGAJE/external/powheg/powheg_beauty_F2_R2.input
new file mode 100644
index 000000000..f4ae5ef68
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_beauty_F2_R2.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 4.75      ! mass of heavy quark in GeV
+facscfact 2     ! factorization scale factor: mufact=muref*facscfact
+renscfact 2     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 10000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_beauty_MHigh.input b/MC/config/PWGGAJE/external/powheg/powheg_beauty_MHigh.input
new file mode 100644
index 000000000..a3a6fe3aa
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_beauty_MHigh.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 5.00      ! mass of heavy quark in GeV
+facscfact 1     ! factorization scale factor: mufact=muref*facscfact
+renscfact 1     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 10000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_beauty_MLow.input b/MC/config/PWGGAJE/external/powheg/powheg_beauty_MLow.input
new file mode 100644
index 000000000..8f689db88
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_beauty_MLow.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 4.50      ! mass of heavy quark in GeV
+facscfact 1     ! factorization scale factor: mufact=muref*facscfact
+renscfact 1     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 10000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_beauty_R05.input b/MC/config/PWGGAJE/external/powheg/powheg_beauty_R05.input
new file mode 100644
index 000000000..0a685eba5
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_beauty_R05.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 4.75      ! mass of heavy quark in GeV
+facscfact 1     ! factorization scale factor: mufact=muref*facscfact
+renscfact 0.5   ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 10000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_beauty_R2.input b/MC/config/PWGGAJE/external/powheg/powheg_beauty_R2.input
new file mode 100644
index 000000000..3209b3e50
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_beauty_R2.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 4.75      ! mass of heavy quark in GeV
+facscfact 1     ! factorization scale factor: mufact=muref*facscfact
+renscfact 2     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 10000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_beauty_default.input b/MC/config/PWGGAJE/external/powheg/powheg_beauty_default.input
new file mode 100644
index 000000000..ab3d3a176
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_beauty_default.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 4.75      ! mass of heavy quark in GeV
+facscfact 1     ! factorization scale factor: mufact=muref*facscfact
+renscfact 1     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 10000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_charm_F05.input b/MC/config/PWGGAJE/external/powheg/powheg_charm_F05.input
new file mode 100644
index 000000000..9031ef445
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_charm_F05.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 1.5       ! mass of heavy quark in GeV
+facscfact 0.5   ! factorization scale factor: mufact=muref*facscfact
+renscfact 1     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 50000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_charm_F05_R05.input b/MC/config/PWGGAJE/external/powheg/powheg_charm_F05_R05.input
new file mode 100644
index 000000000..2366cbbbb
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_charm_F05_R05.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 1.5       ! mass of heavy quark in GeV
+facscfact 0.5   ! factorization scale factor: mufact=muref*facscfact
+renscfact 0.5   ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 50000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_charm_F2.input b/MC/config/PWGGAJE/external/powheg/powheg_charm_F2.input
new file mode 100644
index 000000000..9b1a2466d
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_charm_F2.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 1.5       ! mass of heavy quark in GeV
+facscfact 2     ! factorization scale factor: mufact=muref*facscfact
+renscfact 1     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 50000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_charm_F2_R2.input b/MC/config/PWGGAJE/external/powheg/powheg_charm_F2_R2.input
new file mode 100644
index 000000000..e5f77b0bf
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_charm_F2_R2.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 1.5       ! mass of heavy quark in GeV
+facscfact 2     ! factorization scale factor: mufact=muref*facscfact
+renscfact 2     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 50000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_charm_MHigh.input b/MC/config/PWGGAJE/external/powheg/powheg_charm_MHigh.input
new file mode 100644
index 000000000..214aba28b
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_charm_MHigh.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 1.7       ! mass of heavy quark in GeV
+facscfact 1     ! factorization scale factor: mufact=muref*facscfact
+renscfact 1     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 50000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_charm_MLow.input b/MC/config/PWGGAJE/external/powheg/powheg_charm_MLow.input
new file mode 100644
index 000000000..16df9c118
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_charm_MLow.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 1.3       ! mass of heavy quark in GeV
+facscfact 1     ! factorization scale factor: mufact=muref*facscfact
+renscfact 1     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 50000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_charm_R05.input b/MC/config/PWGGAJE/external/powheg/powheg_charm_R05.input
new file mode 100644
index 000000000..e794c82be
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_charm_R05.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 1.5       ! mass of heavy quark in GeV
+facscfact 1     ! factorization scale factor: mufact=muref*facscfact
+renscfact 0.5   ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 50000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_charm_R2.input b/MC/config/PWGGAJE/external/powheg/powheg_charm_R2.input
new file mode 100644
index 000000000..dfc07b4cd
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_charm_R2.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 1.5       ! mass of heavy quark in GeV
+facscfact 1     ! factorization scale factor: mufact=muref*facscfact
+renscfact 2     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 50000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/external/powheg/powheg_charm_default.input b/MC/config/PWGGAJE/external/powheg/powheg_charm_default.input
new file mode 100644
index 000000000..5d1aaabc8
--- /dev/null
+++ b/MC/config/PWGGAJE/external/powheg/powheg_charm_default.input
@@ -0,0 +1,37 @@
+! Random number generator initializing parameters
+iseed 0         ! initialize random number sequence
+
+! Heavy flavour production parameters
+numevts 100000  ! number of events to be generated
+ih1 1           ! hadron 1 type (1: proton; -1: antiproton)
+ih2 1           ! hadron 2 type (1: proton; -1: antiproton)
+!ndns1 131      ! pdf for hadron 1 (when using the hvqpdf package, hvqpdf numbering)
+!ndns2 131      ! pdf for hadron 2 (when using the hvqpdf package, hvqpdf numbering)
+lhans1 11000    ! pdf set for hadron 1 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+lhans2 11000    ! pdf set for hadron 2 (LHAGLUE number) (11000 = CT10NLO, 10550 = CTEQ66)
+ebeam1 6800     ! energy of beam 1 in GeV
+ebeam2 6800     ! energy of beam 2 in GeV
+qmass 1.5       ! mass of heavy quark in GeV
+facscfact 1     ! factorization scale factor: mufact=muref*facscfact
+renscfact 1     ! renormalization scale factor: muren=muref*renscfact
+!fixedscale 1   ! use ref. scale=qmass (default 0, use running scale)
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1    ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
+use-old-ubound 1  ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
+
+! Parameters that control the grid for Born variables generation
+ncall1 50000    ! number of calls for initializing the integration grid
+itmx1 5         ! number of iterations for initializing the integration grid
+ncall2 100000   ! number of calls for computing the integral and finding upper bound
+itmx2 5         ! number of iterations for computing the integral and finding upper bound
+foldcsi 5       ! number of folds on x integration
+foldy 5         ! number of folds on y integration
+foldphi 1       ! number of folds on phi integration
+nubound 500000  ! number of bbarra calls to setup norm of upper bounding function
+iymax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+ixmax 1         ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
+xupbound 2      ! increase upper bound for radiation generation
+bornktmin 0     ! (default 0d0) Generation cut: minimum kt in underlying Born
+bornsuppfact 0  ! (default 0d0) Mass parameter for Born suppression factor. If < 0 suppfact = 1
+storemintupb 0  ! store function calls in binary file for more precise generation
diff --git a/MC/config/PWGGAJE/ini/GeneratorPythia8POWHEG.ini b/MC/config/PWGGAJE/ini/GeneratorPythia8POWHEG.ini
new file mode 100644
index 000000000..45b8a0bcb
--- /dev/null
+++ b/MC/config/PWGGAJE/ini/GeneratorPythia8POWHEG.ini
@@ -0,0 +1,9 @@
+### The external generator derives from GeneratorPythia8.
+## The generator allows to run Pythia8 with POWHEG
+[GeneratorExternal]
+fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/generator/generator_pythia8_powheg.C
+funcName=getGeneratorJEPythia8POWHEG("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/powheg/powheg_beauty_default.input")
+
+[GeneratorPythia8]
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_powheg.cfg
+
diff --git a/MC/config/PWGGAJE/ini/tests/GeneratorPythia8POWHEG.C b/MC/config/PWGGAJE/ini/tests/GeneratorPythia8POWHEG.C
new file mode 100644
index 000000000..8c83ce1e3
--- /dev/null
+++ b/MC/config/PWGGAJE/ini/tests/GeneratorPythia8POWHEG.C
@@ -0,0 +1,46 @@
+int External() {
+  std::string path{"o2sim_Kine.root"};
+  std::ifstream powhegconf("powheg.input");
+  if (!powhegconf) {
+    std::cerr << "POWHEG configuration file not found\n";
+    return 1;
+  }
+  std::ifstream powhegout("pwgevents.lhe");
+  if (!powhegout) {
+    std::cerr << "POWHEG output file not found\n";
+    return 1;
+  }
+  powhegout.close();
+
+  TFile file(path.c_str(), "READ");
+  if (file.IsZombie()) {
+    std::cerr << "Cannot open ROOT file " << path << "\n";
+    return 1;
+  }
+
+  auto tree = (TTree *)file.Get("o2sim");
+  auto nEvents = tree->GetEntries();
+  std::string line;
+  int nevpowheg = -1;
+  while (std::getline(powhegconf, line)) {
+    if (line.find("numevts") != std::string::npos) {
+      // Read the number right after numevts
+      auto pos = line.find("numevts");
+      nevpowheg = std::stoi(line.substr(pos + 7));
+      if (nevpowheg != nEvents) {
+        std::cerr << "Number of events in POWHEG configuration file " << nevpowheg
+                  << " does not match the simulated number of events "
+                  << nEvents << "\n";
+        return 1;
+      }
+    }
+  }
+  if (nevpowheg == -1) {
+    std::cerr << "Number of events not found in POWHEG configuration file\n";
+    return 1;
+  }
+  powhegconf.close();
+  file.Close();
+
+  return 0;
+}
\ No newline at end of file
diff --git a/MC/config/PWGGAJE/pythia8/generator/pythia8_powheg.cfg b/MC/config/PWGGAJE/pythia8/generator/pythia8_powheg.cfg
new file mode 100644
index 000000000..7a3e2830e
--- /dev/null
+++ b/MC/config/PWGGAJE/pythia8/generator/pythia8_powheg.cfg
@@ -0,0 +1,79 @@
+Beams:eCM 13000. 		# GeV
+### processes
+Next:numberShowLHA = 1
+Next:numberShowInfo = 1
+Next:numberShowProcess = 1
+Next:numberShowEvent = 1
+Main:timesAllowErrors = 10
+! Read LHE file from POWHEG
+Beams:frameType = 4
+Beams:LHEF = pwgevents.lhe
+
+! Number of outgoing particles of POWHEG Born level process
+! (i.e. not counting additional POWHEG radiation)
+POWHEG:nFinal = 2
+
+! How vetoing is performed:
+!  0 - No vetoing is performed (userhooks are not loaded)
+!  1 - Showers are started at the kinematical limit.
+!      Emissions are vetoed if pTemt > pThard.
+!      See also POWHEG:vetoCount below
+POWHEG:veto = 1
+
+! After 'vetoCount' accepted emissions in a row, no more emissions
+! are checked. 'vetoCount = 0' means that no emissions are checked.
+! Use a very large value, e.g. 10000, to have all emissions checked.
+POWHEG:vetoCount = 3
+
+! Selection of pThard (note, for events where there is no
+! radiation, pThard is always set to be SCALUP):
+!  0 - pThard = SCALUP (of the LHA/LHEF standard)
+!  1 - the pT of the POWHEG emission is tested against all other
+!      incoming and outgoing partons, with the minimal value chosen
+!  2 - the pT of all final-state partons is tested against all other
+!      incoming and outgoing partons, with the minimal value chosen
+POWHEG:pThard = 2
+
+! Selection of pTemt:
+!  0 - pTemt is pT of the emitted parton w.r.t. radiating parton
+!  1 - pT of the emission is checked against all incoming and outgoing
+!      partons. pTemt is set to the minimum of these values
+!  2 - the pT of all final-state partons is tested against all other
+!      incoming and outgoing partons, with the minimal value chosen
+! WARNING: the choice here can give significant variations in the final
+! distributions, notably in the tail to large pT values.
+POWHEG:pTemt = 0
+
+! Selection of emitted parton for FSR
+!  0 - Pythia definition of emitted
+!  1 - Pythia definition of radiator
+!  2 - Random selection of emitted or radiator
+!  3 - Both are emitted and radiator are tried
+POWHEG:emitted = 0
+
+! pT definitions
+!  0 - POWHEG ISR pT definition is used for both ISR and FSR
+!  1 - POWHEG ISR pT and FSR d_ij definitions
+!  2 - Pythia definitions
+POWHEG:pTdef = 1
+
+! MPI vetoing
+!  0 - No MPI vetoing is done
+!  1 - When there is no radiation, MPIs with a scale above pT_1 are vetoed,
+!      else MPIs with a scale above (pT_1 + pT_2 + pT_3) / 2 are vetoed
+POWHEG:MPIveto = 0
+! Note that POWHEG:MPIveto = 1 should be combined with
+! MultipartonInteractions:pTmaxMatch = 2
+! which here is taken care of in main31.cc.
+
+! QED vetoing
+!  0 - No QED vetoing is done for pTemt > 0.
+!  1 - QED vetoing is done for pTemt > 0.
+!  2 - QED vetoing is done for pTemt > 0.   If a photon is found
+!      with pT>pThard from the Born level process, the event is accepted
+!      and no further veto of this event is allowed (for any pTemt).
+POWHEG:QEDveto = 2
+
+### decays
+ParticleDecays:limitTau0 on	
+ParticleDecays:tau0Max 10.