diff --git a/MC/config/PWGDQ/external/generator/GeneratorBottomonia.C b/MC/config/PWGDQ/external/generator/GeneratorBottomonia.C
new file mode 100644
index 000000000..25665ede2
--- /dev/null
+++ b/MC/config/PWGDQ/external/generator/GeneratorBottomonia.C
@@ -0,0 +1,280 @@
+//
+// generators for bottomonia considering at midrapidity and forward rapidity
+//
+
+R__ADD_INCLUDE_PATH($O2DPG_ROOT/MC/config/PWGDQ/EvtGen)
+R__LOAD_LIBRARY(libpythia6)
+#include "GeneratorCocktail.C"
+#include "GeneratorEvtGen.C"
+
+namespace o2
+{
+namespace eventgen
+{
+
+/////////////////////////////////////////////////////////////////////////////
+class O2_GeneratorParamUpsilon1SFwdY : public GeneratorTGenerator
+{
+
+ public:
+  O2_GeneratorParamUpsilon1SFwdY() : GeneratorTGenerator("ParamUpsilon1S")
+  {
+    paramUpsilon1S = new GeneratorParam(1, -1, PtUpsilon1Spp13TeV, YUpsilon1Spp13TeV, V2Upsilon1Spp13TeV, IpUpsilon1Spp13TeV);
+    paramUpsilon1S->SetMomentumRange(0., 1.e6);
+    paramUpsilon1S->SetPtRange(0, 999.);
+    paramUpsilon1S->SetYRange(-4.2, -2.3);
+    paramUpsilon1S->SetPhiRange(0., 360.);
+    paramUpsilon1S->SetDecayer(new TPythia6Decayer());
+    paramUpsilon1S->SetForceDecay(kNoDecay); // particle left undecayed
+    // - - paramUpsilon1S->SetTrackingFlag(1);  // (from AliGenParam) -> check this
+    setTGenerator(paramUpsilon1S);
+  };
+
+  ~O2_GeneratorParamUpsilon1SFwdY()
+  {
+    delete paramUpsilon1S;
+  };
+
+  Bool_t Init() override
+  {
+    GeneratorTGenerator::Init();
+    paramUpsilon1S->Init();
+    return true;
+  }
+
+  void SetNSignalPerEvent(Int_t nsig) { paramUpsilon1S->SetNumberParticles(nsig); }
+
+  //-------------------------------------------------------------------------//
+  static Double_t PtUpsilon1Spp13TeV(const Double_t* px, const Double_t* /*dummy*/)
+  {
+    // Upsilon1S pt shape from LHCb pp@13TeV arXiv:1804.09214
+    Double_t x = *px;
+    Float_t p0, p1, p2, p3;
+
+    p0 = 4.11195e+02;
+    p1 = 1.03097e+01;
+    p2 = 1.62309e+00;
+    p3 = 4.84709e+00;
+
+    return (p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t YUpsilon1Spp13TeV(const Double_t* py, const Double_t* /*dummy*/)
+  {
+    // Upsilon1S y shape from LHCb pp@13TeV arXiv:1804.09214
+    Double_t x = *py;
+    Float_t p0, p1;
+
+    p0 =  3.07931e+03;
+    p1 = -3.53102e-02;
+
+    return (p0 * (1. + p1 * x * x));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t V2Upsilon1Spp13TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
+  {
+    // Upsilon(1S) v2
+    return 0.;
+  }
+
+  //-------------------------------------------------------------------------//
+  static Int_t IpUpsilon1Spp13TeV(TRandom*)
+  {
+    return 553;
+  }
+
+ private:
+  GeneratorParam* paramUpsilon1S = nullptr;
+};
+
+/////////////////////////////////////////////////////////////////////////////
+class O2_GeneratorParamUpsilon2SFwdY : public GeneratorTGenerator
+{
+
+ public:
+  O2_GeneratorParamUpsilon2SFwdY() : GeneratorTGenerator("ParamUpsilon2S")
+  {
+    paramUpsilon2S = new GeneratorParam(1, -1, PtUpsilon2Spp13TeV, YUpsilon2Spp13TeV, V2Upsilon2Spp13TeV, IpUpsilon2Spp13TeV);
+    paramUpsilon2S->SetMomentumRange(0., 1.e6);
+    paramUpsilon2S->SetPtRange(0, 999.);
+    paramUpsilon2S->SetYRange(-4.2, -2.3);
+    paramUpsilon2S->SetPhiRange(0., 360.);
+    paramUpsilon2S->SetDecayer(new TPythia6Decayer());
+    paramUpsilon2S->SetForceDecay(kNoDecay); // particle left undecayed
+    // - - paramUpsilon2S->SetTrackingFlag(1);  // (from AliGenParam) -> check this
+    setTGenerator(paramUpsilon2S);
+  };
+
+  ~O2_GeneratorParamUpsilon2SFwdY()
+  {
+    delete paramUpsilon2S;
+  };
+
+  Bool_t Init() override
+  {
+    GeneratorTGenerator::Init();
+    paramUpsilon2S->Init();
+    return true;
+  }
+
+  void SetNSignalPerEvent(Int_t nsig) { paramUpsilon2S->SetNumberParticles(nsig); }
+
+  //-------------------------------------------------------------------------//
+  static Double_t PtUpsilon2Spp13TeV(const Double_t* px, const Double_t* /*dummy*/)
+  {
+    // Upsilon2S pt shape from LHCb pp@13TeV arXiv:1804.09214
+    Double_t x = *px;
+    Float_t p0, p1, p2, p3;
+
+    p0 = 8.15699e+01;
+    p1 = 1.48060e+01;
+    p2 = 1.50018e+00;
+    p3 = 6.34208e+00;
+
+    return (p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t YUpsilon2Spp13TeV(const Double_t* py, const Double_t* /*dummy*/)
+  {
+    // Upsilon2s y shape from LHCb pp@13TeV arXiv:1804.09214
+    Double_t x = *py;
+    Float_t p0, p1;
+
+    p0 =  7.50409e+02;
+    p1 = -3.57039e-02;
+
+    return (p0 * (1. + p1 * x * x));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t V2Upsilon2Spp13TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
+  {
+    // Upsilon(2S) v2
+    return 0.;
+  }
+
+  //-------------------------------------------------------------------------//
+  static Int_t IpUpsilon2Spp13TeV(TRandom*)
+  {
+    return 100553;
+  }
+
+ private:
+  GeneratorParam* paramUpsilon2S = nullptr;
+};
+
+/////////////////////////////////////////////////////////////////////////////
+class O2_GeneratorParamUpsilon3SFwdY : public GeneratorTGenerator
+{
+
+ public:
+  O2_GeneratorParamUpsilon3SFwdY() : GeneratorTGenerator("ParamUpsilon3S")
+  {
+    paramUpsilon3S = new GeneratorParam(1, -1, PtUpsilon3Spp13TeV, YUpsilon3Spp13TeV, V2Upsilon3Spp13TeV, IpUpsilon3Spp13TeV);
+    paramUpsilon3S->SetMomentumRange(0., 1.e6);
+    paramUpsilon3S->SetPtRange(0, 999.);
+    paramUpsilon3S->SetYRange(-4.2, -2.3);
+    paramUpsilon3S->SetPhiRange(0., 360.);
+    paramUpsilon3S->SetDecayer(new TPythia6Decayer());
+    paramUpsilon3S->SetForceDecay(kNoDecay); // particle left undecayed
+    // - - paramUpsilon3S->SetTrackingFlag(1);  // (from AliGenParam) -> check this
+    setTGenerator(paramUpsilon3S);
+  };
+
+  ~O2_GeneratorParamUpsilon3SFwdY()
+  {
+    delete paramUpsilon3S;
+  };
+
+  Bool_t Init() override
+  {
+    GeneratorTGenerator::Init();
+    paramUpsilon3S->Init();
+    return true;
+  }
+
+  void SetNSignalPerEvent(Int_t nsig) { paramUpsilon3S->SetNumberParticles(nsig); }
+
+  //-------------------------------------------------------------------------//
+  static Double_t PtUpsilon3Spp13TeV(const Double_t* px, const Double_t* /*dummy*/)
+  {
+    // Upsilon3S pt shape from LHCb pp@13TeV arXiv:1804.09214
+    Double_t x = *px;
+    Float_t p0, p1, p2, p3;
+
+    p0 = 3.51590e+01;
+    p1 = 2.30813e+01;
+    p2 = 1.40822e+00;
+    p3 = 9.38026e+00;
+
+    return (p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t YUpsilon3Spp13TeV(const Double_t* py, const Double_t* /*dummy*/)
+  {
+    // Upsilon3s y shape from LHCb pp@13TeV arXiv:1804.09214
+    Double_t x = *py;
+    Float_t p0, p1;
+
+    p0 =  3.69961e+02;
+    p1 = -3.54650e-02;
+
+    return (p0 * (1. + p1 * x * x));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t V2Upsilon3Spp13TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
+  {
+    // Upsilon(3S) v2
+    return 0.;
+  }
+
+  //-------------------------------------------------------------------------//
+  static Int_t IpUpsilon3Spp13TeV(TRandom*)
+  {
+    return 200553;
+  }
+
+ private:
+  GeneratorParam* paramUpsilon3S = nullptr;
+};
+
+
+} // namespace eventgen
+} // namespace o2
+
+FairGenerator* GeneratorCocktailBottomoniaToMuonEvtGen_pp13TeV()
+{
+
+  auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen<GeneratorCocktail>();
+
+  auto genUpsilon1S = new o2::eventgen::O2_GeneratorParamUpsilon1SFwdY;
+  genUpsilon1S->SetNSignalPerEvent(1); // 1 Upsilon(1S) generated per event by GeneratorParam
+
+  auto genUpsilon2S = new o2::eventgen::O2_GeneratorParamUpsilon2SFwdY;
+  genUpsilon2S->SetNSignalPerEvent(1); // 1 Upsilon(2S) generated per event by GeneratorParam
+
+  auto genUpsilon3S = new o2::eventgen::O2_GeneratorParamUpsilon3SFwdY;
+  genUpsilon3S->SetNSignalPerEvent(1); // 1 Upsilon(3S) generated per event by GeneratorParam
+
+  genCocktailEvtGen->AddGenerator(genUpsilon1S, 1); // add Upsilon(1S) generator
+  genCocktailEvtGen->AddGenerator(genUpsilon2S, 1); // add Upsilon(2S) generator
+  genCocktailEvtGen->AddGenerator(genUpsilon3S, 1); // add Upsilon(3S) generator
+
+  TString pdgs = "553;100553;200553";
+  std::string spdg;
+  TObjArray* obj = pdgs.Tokenize(";");
+  genCocktailEvtGen->SetSizePdg(obj->GetEntriesFast());
+  for (int i = 0; i < obj->GetEntriesFast(); i++) {
+    spdg = obj->At(i)->GetName();
+    genCocktailEvtGen->AddPdg(std::stoi(spdg), i);
+    printf("PDG %d \n", std::stoi(spdg));
+  }
+  genCocktailEvtGen->SetForceDecay(kEvtDiMuon);
+
+  return genCocktailEvtGen;
+}
\ No newline at end of file
diff --git a/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedBottomoniaSignals_gaptriggered_dq.C b/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedBottomoniaSignals_gaptriggered_dq.C
new file mode 100644
index 000000000..de4ab2b7a
--- /dev/null
+++ b/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedBottomoniaSignals_gaptriggered_dq.C
@@ -0,0 +1,76 @@
+#include "FairGenerator.h"
+#include "Generators/GeneratorPythia8.h"
+#include "Pythia8/Pythia.h"
+#include "TRandom.h"
+#include "GeneratorBottomonia.C"
+#include <string>
+
+using namespace o2::eventgen;
+using namespace Pythia8;
+
+class GeneratorPythia8BottomoniaInjectedGapTriggeredDQ : public o2::eventgen::GeneratorPythia8 {
+public:
+
+  /// default constructor
+  GeneratorPythia8BottomoniaInjectedGapTriggeredDQ() = default;
+
+  /// constructor
+  GeneratorPythia8BottomoniaInjectedGapTriggeredDQ(int inputTriggerRatio = 5, int gentype = 0) {
+
+    mGeneratedEvents = 0;
+    mGeneratorParam = 0x0;
+    mInverseTriggerRatio = inputTriggerRatio; 
+    switch (gentype) {
+      case 0: // generate bottomonia cocktail at forward rapidity
+        mGeneratorParam = (Generator*)GeneratorCocktailBottomoniaToMuonEvtGen_pp13TeV(); 
+        break;
+    }
+    mGeneratorParam->Init();  
+  }
+
+  ///  Destructor
+  ~GeneratorPythia8BottomoniaInjectedGapTriggeredDQ() = default;
+
+protected:
+
+Bool_t importParticles() override
+  {
+    GeneratorPythia8::importParticles();
+    bool genOk = false;
+    if (mGeneratedEvents % mInverseTriggerRatio == 0) { // add injected prompt signals to the stack
+      bool genOk = false;
+      while (!genOk) {    
+        genOk =  (mGeneratorParam->generateEvent() && mGeneratorParam->importParticles()) ? true : false ;
+      }
+      int originalSize = mParticles.size();
+      for (int ipart=0; ipart < mGeneratorParam->getParticles().size(); ipart++) {
+        TParticle part = TParticle(mGeneratorParam->getParticles().at(ipart));
+        if(part.GetFirstMother() >= 0) part.SetFirstMother(part.GetFirstMother() + originalSize);
+        if(part.GetFirstDaughter() >= 0) part.SetFirstDaughter(part.GetFirstDaughter() + originalSize);
+        if(part.GetLastDaughter() >= 0) part.SetLastDaughter(part.GetLastDaughter() + originalSize);
+        mParticles.push_back(part); 
+        // encodeParticleStatusAndTracking method already called in GeneratorEvtGen.C 
+      }	   
+      mGeneratorParam->clearParticles(); 
+    }
+
+    mGeneratedEvents++;
+    return true;
+  }
+
+
+private:
+  Generator* mGeneratorParam; 
+  // Control gap-triggering
+  unsigned long long mGeneratedEvents;
+  int mInverseTriggerRatio;
+};
+
+// Predefined generators:
+FairGenerator *GeneratorPythia8InjectedBottomoniaGapTriggered(int inputTriggerRatio, int gentype) {
+  auto myGen = new GeneratorPythia8BottomoniaInjectedGapTriggeredDQ(inputTriggerRatio,gentype);
+  auto seed = (gRandom->TRandom::GetSeed() % 900000000);
+  myGen->readString("Random:setSeed on");
+  myGen->readString("Random:seed " + std::to_string(seed));
+  return myGen;
+}
diff --git a/MC/config/PWGDQ/ini/Generator_InjectedBottomoniaFwdy_TriggerGap.ini b/MC/config/PWGDQ/ini/Generator_InjectedBottomoniaFwdy_TriggerGap.ini
new file mode 100755
index 000000000..029230bd0
--- /dev/null
+++ b/MC/config/PWGDQ/ini/Generator_InjectedBottomoniaFwdy_TriggerGap.ini
@@ -0,0 +1,7 @@
+### The external generator derives from GeneratorPythia8.
+[GeneratorExternal]
+fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedBottomoniaSignals_gaptriggered_dq.C
+funcName=GeneratorPythia8InjectedBottomoniaGapTriggered(0,3)
+
+[GeneratorPythia8]
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg
diff --git a/MC/config/PWGDQ/ini/tests/Generator_InjectedBottomoniaFwdy_TriggerGap.C b/MC/config/PWGDQ/ini/tests/Generator_InjectedBottomoniaFwdy_TriggerGap.C
new file mode 100644
index 000000000..9641dc2d8
--- /dev/null
+++ b/MC/config/PWGDQ/ini/tests/Generator_InjectedBottomoniaFwdy_TriggerGap.C
@@ -0,0 +1,101 @@
+int External()
+{
+    int checkPdgSignal[] = {553,100553,200553};
+    int checkPdgDecay = 13;
+    double rapiditymin = -4.3; double rapiditymax = -2.3;
+    std::string path{"o2sim_Kine.root"};
+    std::cout << "Check for\nsignal PDG " << checkPdgSignal << "\ndecay PDG " << checkPdgDecay << "\n";
+    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");
+    std::vector<o2::MCTrack>* tracks{};
+    tree->SetBranchAddress("MCTrack", &tracks);
+
+    int nLeptons{};
+    int nAntileptons{};
+    int nLeptonPairs{};
+    int nLeptonPairsToBeDone{};
+    int nSignalUpsilon1S{};
+    int nSignalUpsilon2S{};
+    int nSignalUpsilon1SWithinAcc{};
+    int nSignalUpsilon2SWithinAcc{};
+    auto nEvents = tree->GetEntries();
+    o2::steer::MCKinematicsReader mcreader("o2sim", o2::steer::MCKinematicsReader::Mode::kMCKine);
+    Bool_t isInjected = kFALSE;
+
+    for (int i = 0; i < nEvents; i++) {
+        tree->GetEntry(i);
+        for (auto& track : *tracks) {
+            auto pdg = track.GetPdgCode();
+	        auto rapidity =  track.GetRapidity();
+	        auto idMoth = track.getMotherTrackId();
+            if (pdg == checkPdgDecay) {
+                // count leptons
+                nLeptons++;
+            } else if(pdg == -checkPdgDecay) {
+                // count anti-leptons
+                nAntileptons++;
+            } else if (pdg == checkPdgSignal[0] || pdg == checkPdgSignal[1] || pdg == checkPdgSignal[2]) {
+                if(idMoth < 0){
+                    // count signal PDG 
+                    if (pdg == checkPdgSignal[0]) {
+                        nSignalUpsilon1S++;
+                    } else if (pdg == checkPdgSignal[1]) {
+                        nSignalUpsilon2S++;
+                    } else {
+                        nSignalUpsilon3S++;
+                    }
+                    
+                    // count signal PDG within acceptance 
+                    if (rapidity > rapiditymin && rapidity < rapiditymax) {
+                        if (pdg == checkPdgSignal[0]) {
+                            nSignalUpsilon1SWithinAcc++;
+                        } else if (pdg == checkPdgSignal[1]) {
+                            nSignalUpsilon2SWithinAcc++;
+                        } else {
+                            nSignalUpsilon3SWithinAcc++;
+                        }
+                    }
+		        }
+		        auto child0 = o2::mcutils::MCTrackNavigator::getDaughter0(track, *tracks);
+                auto child1 = o2::mcutils::MCTrackNavigator::getDaughter1(track, *tracks);
+                if (child0 != nullptr && child1 != nullptr) {
+                    // check for parent-child relations
+                    auto pdg0 = child0->GetPdgCode();
+                    auto pdg1 = child1->GetPdgCode();
+                    std::cout << "First and last children of parent " << checkPdgSignal << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
+                    if (std::abs(pdg0) == checkPdgDecay && std::abs(pdg1) == checkPdgDecay && pdg0 == -pdg1) {
+                        nLeptonPairs++;
+                        if (child0->getToBeDone() && child1->getToBeDone()) {
+                            nLeptonPairsToBeDone++;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    std::cout << "#events: " << nEvents << "\n"
+              << "#leptons: " << nLeptons << "\n"
+              << "#antileptons: " << nAntileptons << "\n"
+              << "#signal (Upsilon(1S)): " << nSignalUpsilon1S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalUpsilon1SWithinAcc << "\n"
+              << "#signal (Upsilon(2S)): " << nSignalUpsilon2S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalUpsilon2SWithinAcc << "\n"
+              << "#signal (Upsilon(2S)): " << nSignalUpsilon3S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalUpsilon3SWithinAcc << "\n"
+              << "#lepton pairs: " << nLeptonPairs << "\n"
+              << "#lepton pairs to be done: " << nLeptonPairs << "\n";
+
+
+    if (nLeptonPairs == 0 || nLeptons == 0 || nAntileptons == 0) {
+        std::cerr << "Number of leptons, number of anti-leptons as well as number of lepton pairs should all be greater than 1.\n";
+        return 1;
+    }
+    if (nLeptonPairs != nLeptonPairsToBeDone) {
+        std::cerr << "The number of lepton pairs should be the same as the number of lepton pairs which should be transported.\n";
+        return 1;
+    }
+
+    return 0;
+}
\ No newline at end of file