diff --git a/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_HFLepton_OO.C b/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_HFLepton_OO.C
new file mode 100644
index 000000000..55a23a437
--- /dev/null
+++ b/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_HFLepton_OO.C
@@ -0,0 +1,234 @@
+#include "Pythia8/Pythia.h"
+#include "Pythia8/HeavyIons.h"
+#include "FairGenerator.h"
+#include "FairPrimaryGenerator.h"
+#include "Generators/GeneratorPythia8.h"
+#include "TRandom3.h"
+#include "TParticlePDG.h"
+#include "TDatabasePDG.h"
+
+#include <map>
+#include <unordered_set>
+//#include <utility>	// for std::pair
+
+using namespace Pythia8;
+
+class GeneratorPythia8GapTriggeredHFLeptonOO : public o2::eventgen::GeneratorPythia8
+{
+public:
+  /// default constructor
+  GeneratorPythia8GapTriggeredHFLeptonOO() = default;
+
+  /// constructor
+  GeneratorPythia8GapTriggeredHFLeptonOO(TString configsignal, int quarkPdg = 4, int lInputTriggerRatio = 5, int lInputExternalID = 0)
+  {
+
+    lGeneratedEvents = 0;
+    lInverseTriggerRatio = lInputTriggerRatio;
+    lExternalID = lInputExternalID;
+    mQuarkPdg = quarkPdg;
+
+    auto seed = (gRandom->TRandom::GetSeed() % 900000000);
+
+    int offset = (int)(gRandom->Uniform(lInverseTriggerRatio)); // create offset to mitigate edge effects due to small number of events per job
+    lGeneratedEvents += offset;
+    
+    cout << "Initalizing extra PYTHIA object used to generate min-bias events..." << endl;
+    TString pathconfigMB = gSystem->ExpandPathName("${O2DPG_MC_CONFIG_ROOT}//MC/config/common/pythia8/generator/pythia8_OO_536.cfg");
+    pythiaObjectMinimumBias.readFile(pathconfigMB.Data());
+    pythiaObjectMinimumBias.readString("Random:setSeed on");
+    pythiaObjectMinimumBias.readString("Random:seed " + std::to_string(seed));
+    pythiaObjectMinimumBias.init();
+    cout << "Initalization of gap event is complete" << endl;
+
+    cout << "Initalizing extra PYTHIA object used to generate signal events..." << endl;
+    TString pathconfigSignal = gSystem->ExpandPathName(configsignal.Data());
+    pythiaObjectSignal.readFile(pathconfigSignal.Data());
+    pythiaObjectSignal.readString("Random:setSeed on");
+    pythiaObjectSignal.readString("Random:seed " + std::to_string(seed));
+    pythiaObjectSignal.readString("Beams:idA = 1000080160");
+    pythiaObjectSignal.readString("Beams:idB = 1000080160");
+    pythiaObjectSignal.readString("Beams:eCM = 5360.0");
+    pythiaObjectSignal.readString("Beams:frameType = 1");
+    pythiaObjectSignal.readString("ParticleDecays:limitTau0 = on");
+    pythiaObjectSignal.readString("ParticleDecays:tau0Max = 10.");
+    pythiaObjectSignal.readString("HeavyIon:SigFitNGen = 0");
+    pythiaObjectSignal.readString("HeavyIon:SigFitDefPar = 2.15,18.42,0.33");
+    pythiaObjectSignal.init();
+    cout << "Initalization of signal event is complete" << endl;
+
+    // flag the generators using type
+    // addCocktailConstituent(type, "interesting");
+    // addCocktailConstitent(0, "minbias");
+    // Add Sub generators
+    addSubGenerator(0, "default generator");
+    addSubGenerator(1, "charm lepton");
+    addSubGenerator(2, "beauty forced decay");
+    addSubGenerator(3, "beauty no foced decay");
+  }
+
+  ///  Destructor
+  ~GeneratorPythia8GapTriggeredHFLeptonOO() = default;
+
+  void addTriggerOnDaughter(int nb, int pdg)
+  {
+    mNbDaughter = nb;
+    mPdgDaughter = pdg;
+  };
+  void setQuarkRapidity(float yMin, float yMax)
+  {
+    mQuarkRapidityMin = yMin;
+    mQuarkRapidityMax = yMax;
+  };
+  void setDaughterRapidity(float yMin, float yMax)
+  {
+    mDaughterRapidityMin = yMin;
+    mDaughterRapidityMax = yMax;
+  };
+
+protected:
+  //__________________________________________________________________
+  Bool_t generateEvent() override
+  {
+    /// reset event
+    mPythia.event.reset();
+
+    // Simple straightforward check to alternate generators
+    if (lGeneratedEvents % lInverseTriggerRatio == 0) {
+      // Generate event of interest
+      Bool_t lGenerationOK = kFALSE;
+      while (!lGenerationOK) {
+        if (pythiaObjectSignal.next()) {
+          lGenerationOK = selectEvent(pythiaObjectSignal.event);
+        }
+      }
+      mPythia.event = pythiaObjectSignal.event;
+      notifySubGenerator(lExternalID);
+    } else {
+      // Generate minimum-bias event
+      Bool_t lGenerationOK = kFALSE;
+      while (!lGenerationOK) {
+        lGenerationOK = pythiaObjectMinimumBias.next();
+      }
+      mPythia.event = pythiaObjectMinimumBias.event;
+      notifySubGenerator(0);
+    }
+
+    lGeneratedEvents++;
+    // mPythia.next();
+
+    return true;
+  }
+
+  bool selectEvent(const Pythia8::Event& event)
+  {
+    bool isGoodAtPartonLevel = false, isGoodAtDaughterLevel = (mPdgDaughter != 0) ? false : true;
+    int nbDaughter = 0;
+    for (auto iPart{0}; iPart < event.size(); ++iPart) {
+      // search for Q-Qbar mother with at least one Q in rapidity window
+      if (!isGoodAtPartonLevel) {
+        auto daughterList = event[iPart].daughterList();
+        bool hasQ = false, hasQbar = false, atSelectedY = false;
+        for (auto iDau : daughterList) {
+          if (event[iDau].id() == mQuarkPdg) {
+            hasQ = true;
+          }
+          if (event[iDau].id() == -mQuarkPdg) {
+            hasQbar = true;
+          }
+          if ((std::abs(event[iDau].id()) == mQuarkPdg) && (event[iDau].y() > mQuarkRapidityMin) && (event[iDau].y() < mQuarkRapidityMax))
+            atSelectedY = true;
+        }
+        if (hasQ && hasQbar && atSelectedY) {
+          isGoodAtPartonLevel = true;
+        }
+      }
+      // search for mNbDaughter daughters of type mPdgDaughter in rapidity window
+      if (!isGoodAtDaughterLevel) {
+        int id = std::abs(event[iPart].id());
+        float rap = event[iPart].y();
+        if (id == mPdgDaughter) {
+          int motherindexa = event[iPart].mother1();
+          if (motherindexa > 0) {
+            int idmother = std::abs(event[motherindexa].id());
+            if (int(std::abs(idmother) / 100.) == 4 || int(std::abs(idmother) / 1000.) == 4 || int(std::abs(idmother) / 100.) == 5 || int(std::abs(idmother) / 1000.) == 5) {
+              if (rap > mDaughterRapidityMin && rap < mDaughterRapidityMax) {
+                nbDaughter++;
+                if (nbDaughter >= mNbDaughter) isGoodAtDaughterLevel = true;
+              }
+            }
+          }
+        }
+      }
+      // we send the trigger
+      if (isGoodAtPartonLevel && isGoodAtDaughterLevel) {
+        return true;
+      }
+    }
+    return false;
+  };
+
+private:
+  // Interface to override import particles
+  Pythia8::Event mOutputEvent;
+
+  // Properties of selection
+  int mQuarkPdg;
+  float mQuarkRapidityMin;
+  float mQuarkRapidityMax;
+  int mPdgDaughter;
+  int mNbDaughter;
+  float mDaughterRapidityMin;
+  float mDaughterRapidityMax;
+
+  // Control gap-triggering
+  Long64_t lGeneratedEvents;
+  int lInverseTriggerRatio;
+  // ID for different generators
+  int lExternalID;
+
+  // Base event generators
+  Pythia8::Pythia pythiaObjectMinimumBias; ///Minimum bias collision generator
+  Pythia8::Pythia pythiaObjectSignal;      ///Signal collision generator
+};
+
+// Predefined generators:
+
+// Charm-enriched forced decay
+FairGenerator* GeneratorPythia8GapTriggeredCharmLepton(int inputTriggerRatio, int inputExternalID, int pdgLepton, float yMinQ = -1.5, float yMaxQ = 1.5, float yMinL = -1, float yMaxL = 1)
+{
+  auto myGen = new GeneratorPythia8GapTriggeredHFLeptonOO("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_pp_cr2_forceddecayscharm.cfg", 4, inputTriggerRatio, inputExternalID);
+  auto seed = (gRandom->TRandom::GetSeed() % 900000000);
+  myGen->readString("Random:setSeed on");
+  myGen->readString("Random:seed " + std::to_string(seed));
+  myGen->setQuarkRapidity(yMinQ, yMaxQ);
+  myGen->addTriggerOnDaughter(2, pdgLepton);
+  myGen->setDaughterRapidity(yMinL, yMaxL);
+  return myGen;
+}
+
+// Beauty-enriched forced decay
+FairGenerator* GeneratorPythia8GapTriggeredBeautyForcedDecays(int inputTriggerRatio, int inputExternalID, int pdgLepton, float yMinQ = -1.5, float yMaxQ = 1.5, float yMinL = -1, float yMaxL = 1)
+{
+  auto myGen = new GeneratorPythia8GapTriggeredHFLeptonOO("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_bbbar_forceddecayscharmbeauty.cfg", 5, inputTriggerRatio, inputExternalID);
+  auto seed = (gRandom->TRandom::GetSeed() % 900000000);
+  myGen->readString("Random:setSeed on");
+  myGen->readString("Random:seed " + std::to_string(seed));
+  myGen->setQuarkRapidity(yMinQ, yMaxQ);
+  myGen->addTriggerOnDaughter(2, pdgLepton);
+  myGen->setDaughterRapidity(yMinL, yMaxL);
+  return myGen;
+}
+
+// Beauty-enriched no forced decay
+FairGenerator* GeneratorPythia8GapTriggeredBeautyNoForcedDecays(int inputTriggerRatio, int inputExternalID, int pdgLepton, float yMinQ = -1.5, float yMaxQ = 1.5, float yMinL = -1, float yMaxL = 1)
+{
+  auto myGen = new GeneratorPythia8GapTriggeredHFLeptonOO("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_bbbar.cfg", 5, inputTriggerRatio, inputExternalID);
+  auto seed = (gRandom->TRandom::GetSeed() % 900000000);
+  myGen->readString("Random:setSeed on");
+  myGen->readString("Random:seed " + std::to_string(seed));
+  myGen->setQuarkRapidity(yMinQ, yMaxQ);
+  myGen->addTriggerOnDaughter(2, pdgLepton);
+  myGen->setDaughterRapidity(yMinL, yMaxL);
+  return myGen;
+}
diff --git a/MC/config/PWGEM/ini/GeneratorHFGapTriggered_BeautyForcedDecay_Gap2_OO_electron.ini b/MC/config/PWGEM/ini/GeneratorHFGapTriggered_BeautyForcedDecay_Gap2_OO_electron.ini
new file mode 100644
index 000000000..f25401e61
--- /dev/null
+++ b/MC/config/PWGEM/ini/GeneratorHFGapTriggered_BeautyForcedDecay_Gap2_OO_electron.ini
@@ -0,0 +1,8 @@
+#NEV_TEST> 5
+[GeneratorExternal]
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_HFLepton_OO.C
+funcName = GeneratorPythia8GapTriggeredBeautyForcedDecays(2, 2, 11, -1.5, +1.5, -1.2, +1.2)
+
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg
+includePartonEvent = true
diff --git a/MC/config/PWGEM/ini/GeneratorHFGapTriggered_BeautyForcedDecay_Gap2_OO_muon.ini b/MC/config/PWGEM/ini/GeneratorHFGapTriggered_BeautyForcedDecay_Gap2_OO_muon.ini
new file mode 100644
index 000000000..b129f9fe8
--- /dev/null
+++ b/MC/config/PWGEM/ini/GeneratorHFGapTriggered_BeautyForcedDecay_Gap2_OO_muon.ini
@@ -0,0 +1,8 @@
+#NEV_TEST> 5
+[GeneratorExternal]
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_HFLepton_OO.C
+funcName = GeneratorPythia8GapTriggeredBeautyForcedDecays(2, 2, 13, -6, 0, -5, -1)
+
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg
+includePartonEvent = true
diff --git a/MC/config/PWGEM/ini/GeneratorHFGapTriggered_BeautyNoForcedDecay_Gap2_OO_electron.ini b/MC/config/PWGEM/ini/GeneratorHFGapTriggered_BeautyNoForcedDecay_Gap2_OO_electron.ini
new file mode 100644
index 000000000..8cd7e5ab6
--- /dev/null
+++ b/MC/config/PWGEM/ini/GeneratorHFGapTriggered_BeautyNoForcedDecay_Gap2_OO_electron.ini
@@ -0,0 +1,8 @@
+#NEV_TEST> 5
+[GeneratorExternal]
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_HFLepton_OO.C
+funcName = GeneratorPythia8GapTriggeredBeautyNoForcedDecays(2, 3, 11, -1.5, +1.5, -1.2, +1.2)
+
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg
+includePartonEvent = true
diff --git a/MC/config/PWGEM/ini/GeneratorHFGapTriggered_BeautyNoForcedDecay_Gap2_OO_muon.ini b/MC/config/PWGEM/ini/GeneratorHFGapTriggered_BeautyNoForcedDecay_Gap2_OO_muon.ini
new file mode 100644
index 000000000..97741b31b
--- /dev/null
+++ b/MC/config/PWGEM/ini/GeneratorHFGapTriggered_BeautyNoForcedDecay_Gap2_OO_muon.ini
@@ -0,0 +1,8 @@
+#NEV_TEST> 5
+[GeneratorExternal]
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_HFLepton_OO.C
+funcName = GeneratorPythia8GapTriggeredBeautyNoForcedDecays(2, 3, 13, -6, 0, -5, -1)
+
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg
+includePartonEvent = true
diff --git a/MC/config/PWGEM/ini/GeneratorHFGapTriggered_Charm_Gap2_OO_electron.ini b/MC/config/PWGEM/ini/GeneratorHFGapTriggered_Charm_Gap2_OO_electron.ini
new file mode 100644
index 000000000..9691afc04
--- /dev/null
+++ b/MC/config/PWGEM/ini/GeneratorHFGapTriggered_Charm_Gap2_OO_electron.ini
@@ -0,0 +1,8 @@
+#NEV_TEST> 5
+[GeneratorExternal]
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_HFLepton_OO.C
+funcName = GeneratorPythia8GapTriggeredCharmLepton(2, 1, 11, -1.5, +1.5, -1.2, +1.2)
+
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg
+includePartonEvent = true
diff --git a/MC/config/PWGEM/ini/GeneratorHFGapTriggered_Charm_Gap2_OO_muon.ini b/MC/config/PWGEM/ini/GeneratorHFGapTriggered_Charm_Gap2_OO_muon.ini
new file mode 100644
index 000000000..59d5babf4
--- /dev/null
+++ b/MC/config/PWGEM/ini/GeneratorHFGapTriggered_Charm_Gap2_OO_muon.ini
@@ -0,0 +1,8 @@
+#NEV_TEST> 5
+[GeneratorExternal]
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_HFLepton_OO.C
+funcName = GeneratorPythia8GapTriggeredCharmLepton(2, 1, 13, -6, 0, -5, -1)
+
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg
+includePartonEvent = true
diff --git a/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_BeautyForcedDecay_Gap2_OO_electron.C b/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_BeautyForcedDecay_Gap2_OO_electron.C
new file mode 100644
index 000000000..d4c7056e4
--- /dev/null
+++ b/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_BeautyForcedDecay_Gap2_OO_electron.C
@@ -0,0 +1,117 @@
+int External()
+{
+  int checkPdgDecay = 11;
+  std::string path{"o2sim_Kine.root"};
+  TFile file(path.c_str(), "READ");
+  if (file.IsZombie()) {
+    std::cerr << "Cannot open ROOT file " << path << "\n";
+    return 1;
+  }
+
+  float ratioTrigger = 1./3; // one event triggered out of 3
+  auto tree = (TTree*)file.Get("o2sim");
+  std::vector<o2::MCTrack>* tracks{};
+  tree->SetBranchAddress("MCTrack", &tracks);
+
+  int nLeptonsInAcceptance{};
+  int nLeptons{};
+  int nAntileptons{};
+  int nLeptonsToBeDone{};
+  int nAntileptonsToBeDone{};
+  int nSignalPairs{};
+  int nLeptonPairs{};
+  int nLeptonPairsToBeDone{};
+  auto nEvents = tree->GetEntries();
+
+  for (int i = 0; i < nEvents; i++) {
+    printf("Event %d\n",i);
+    tree->GetEntry(i);
+    int nElectrons = 0;
+    int nPositrons = 0;
+    int nElectronsToBeDone = 0;
+    int nPositronsToBeDone = 0;
+    int nOpenBeautyPos = 0;
+    int nOpenBeautyNeg = 0;
+    int nPositronsElectronsInAcceptance = 0;
+    for (auto& track : *tracks) {
+      auto pdg = track.GetPdgCode();
+      auto y = track.GetRapidity();
+      if (pdg == checkPdgDecay) {
+        int igmother = track.getMotherTrackId();
+        if (igmother > 0) {
+          auto gmTrack = (*tracks)[igmother];
+          int gmpdg = gmTrack.GetPdgCode();
+          if (int(std::abs(gmpdg)/100.) == 5 || int(std::abs(gmpdg)/1000.) == 5 || int(std::abs(gmpdg)/100.) == 4 || int(std::abs(gmpdg)/1000.) == 4) {
+            nLeptons++;
+            nElectrons++;
+            if (-1 < y && y < 1) nPositronsElectronsInAcceptance++;
+            if (track.getToBeDone()) {
+              nLeptonsToBeDone++;
+              nElectronsToBeDone++;
+            }
+          }
+        }    
+      } else if (pdg == -checkPdgDecay) {
+        int igmother = track.getMotherTrackId();
+        if (igmother > 0) {
+          auto gmTrack = (*tracks)[igmother];
+          int gmpdg = gmTrack.GetPdgCode();
+          if (int(TMath::Abs(gmpdg)/100.) == 4 || int(TMath::Abs(gmpdg)/1000.) == 4 || int(std::abs(gmpdg)/100.) == 5 || int(std::abs(gmpdg)/1000.) == 5) {
+            nAntileptons++;
+            nPositrons++;
+            if (-1 < y && y < 1) nPositronsElectronsInAcceptance++;
+            if (track.getToBeDone()) {
+              nAntileptonsToBeDone++;
+              nPositronsToBeDone++;
+            }
+          }
+        }
+      } else if (pdg == 511 || pdg == 521 || pdg == 531 || pdg == 5122 || pdg == 5132 || pdg == 5232 || pdg == 5332) {
+        nOpenBeautyPos++;
+      }  else if (pdg == -511 || pdg == -521 || pdg == -531 || pdg == -5122 || pdg == -5132 || pdg == -5232 || pdg == -5332) {
+        nOpenBeautyNeg++;
+      }
+    }
+    if (nOpenBeautyPos > 0 && nOpenBeautyNeg > 0) {
+      nSignalPairs++;
+    }
+    if (nPositronsElectronsInAcceptance > 1) {
+      nLeptonsInAcceptance++;
+    }
+    if (nElectrons > 0 && nPositrons > 0) {
+      nLeptonPairs++;
+    }
+    if (nElectronsToBeDone > 0 && nPositronsToBeDone > 0) nLeptonPairsToBeDone++;
+  }
+  std::cout << "#events: " << nEvents << "\n"
+    << "#leptons: " << nLeptons << "\n"
+    << "#antileptons: " << nAntileptons << "\n"
+    << "#leptons to be done: " << nLeptonsToBeDone << "\n"
+    << "#antileptons to be done: " << nAntileptonsToBeDone << "\n"
+    << "#Open-beauty hadron pairs: " << nSignalPairs << "\n"
+    << "#leptons in acceptance: " << nLeptonsInAcceptance << "\n"
+    << "#Electron-positron pairs " << nLeptonPairs << "\n"
+    << "#Electron-positron pairs to be done: " << nLeptonPairsToBeDone << "\n";
+  if (nLeptons == 0 && nAntileptons == 0) {
+    std::cerr << "Number of leptons, number of anti-leptons should all be greater than 1.\n";
+    return 1;
+  }
+  if (nLeptonPairs != nLeptonPairsToBeDone) {
+    std::cerr << "The number of electron-positron pairs should be the same as the number of electron-positron pairs which should be transported.\n";
+    return 1;
+  }
+  if (nLeptons != nLeptonsToBeDone) {
+    std::cerr << "The number of leptons should be the same as the number of leptons which should be transported.\n";
+    return 1;
+  }
+  if (nLeptonsInAcceptance ==  (nEvents/ratioTrigger)) {
+    std::cerr << "The number of leptons in acceptance should be at least equaled to the number of events.\n";
+    return 1;
+  }
+  if (nLeptonPairs <  nLeptonsInAcceptance) {
+    std::cerr << "The number of positron-electron pairs should be at least equaled to the number of leptons in acceptance.\n";
+    return 1;
+  }
+
+  return 0;
+}
diff --git a/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_BeautyForcedDecay_Gap2_OO_muon.C b/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_BeautyForcedDecay_Gap2_OO_muon.C
new file mode 100644
index 000000000..dea67792e
--- /dev/null
+++ b/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_BeautyForcedDecay_Gap2_OO_muon.C
@@ -0,0 +1,117 @@
+int External()
+{
+  int checkPdgDecay = 13;
+  std::string path{"o2sim_Kine.root"};
+  TFile file(path.c_str(), "READ");
+  if (file.IsZombie()) {
+    std::cerr << "Cannot open ROOT file " << path << "\n";
+    return 1;
+  }
+
+  float ratioTrigger = 1./3; // one event triggered out of 3
+  auto tree = (TTree*)file.Get("o2sim");
+  std::vector<o2::MCTrack>* tracks{};
+  tree->SetBranchAddress("MCTrack", &tracks);
+
+  int nLeptonsInAcceptance{};
+  int nLeptons{};
+  int nAntileptons{};
+  int nLeptonsToBeDone{};
+  int nAntileptonsToBeDone{};
+  int nSignalPairs{};
+  int nLeptonPairs{};
+  int nLeptonPairsToBeDone{};
+  auto nEvents = tree->GetEntries();
+
+  for (int i = 0; i < nEvents; i++) {
+    printf("Event %d\n",i);
+    tree->GetEntry(i);
+    int nElectrons = 0;
+    int nPositrons = 0;
+    int nElectronsToBeDone = 0;
+    int nPositronsToBeDone = 0;
+    int nOpenBeautyPos = 0;
+    int nOpenBeautyNeg = 0;
+    int nPositronsElectronsInAcceptance = 0;
+    for (auto& track : *tracks) {
+      auto pdg = track.GetPdgCode();
+      auto y = track.GetRapidity();
+      if (pdg == checkPdgDecay) {
+        int igmother = track.getMotherTrackId();
+        if (igmother > 0) {
+          auto gmTrack = (*tracks)[igmother];
+          int gmpdg = gmTrack.GetPdgCode();
+          if (int(std::abs(gmpdg)/100.) == 5 || int(std::abs(gmpdg)/1000.) == 5 || int(std::abs(gmpdg)/100.) == 4 || int(std::abs(gmpdg)/1000.) == 4) {
+            nLeptons++;
+            nElectrons++;
+            if (-1 < y && y < 1) nPositronsElectronsInAcceptance++;
+            if (track.getToBeDone()) {
+              nLeptonsToBeDone++;
+              nElectronsToBeDone++;
+            }
+          }
+        }    
+      } else if (pdg == -checkPdgDecay) {
+        int igmother = track.getMotherTrackId();
+        if (igmother > 0) {
+          auto gmTrack = (*tracks)[igmother];
+          int gmpdg = gmTrack.GetPdgCode();
+          if (int(TMath::Abs(gmpdg)/100.) == 4 || int(TMath::Abs(gmpdg)/1000.) == 4 || int(std::abs(gmpdg)/100.) == 5 || int(std::abs(gmpdg)/1000.) == 5) {
+            nAntileptons++;
+            nPositrons++;
+            if (-1 < y && y < 1) nPositronsElectronsInAcceptance++;
+            if (track.getToBeDone()) {
+              nAntileptonsToBeDone++;
+              nPositronsToBeDone++;
+            }
+          }
+        }
+      } else if (pdg == 511 || pdg == 521 || pdg == 531 || pdg == 5122 || pdg == 5132 || pdg == 5232 || pdg == 5332) {
+        nOpenBeautyPos++;
+      }  else if (pdg == -511 || pdg == -521 || pdg == -531 || pdg == -5122 || pdg == -5132 || pdg == -5232 || pdg == -5332) {
+        nOpenBeautyNeg++;
+      }
+    }
+    if (nOpenBeautyPos > 0 && nOpenBeautyNeg > 0) {
+      nSignalPairs++;
+    }
+    if (nPositronsElectronsInAcceptance > 1) {
+      nLeptonsInAcceptance++;
+    }
+    if (nElectrons > 0 && nPositrons > 0) {
+      nLeptonPairs++;
+    }
+    if (nElectronsToBeDone > 0 && nPositronsToBeDone > 0) nLeptonPairsToBeDone++;
+  }
+  std::cout << "#events: " << nEvents << "\n"
+    << "#leptons: " << nLeptons << "\n"
+    << "#antileptons: " << nAntileptons << "\n"
+    << "#leptons to be done: " << nLeptonsToBeDone << "\n"
+    << "#antileptons to be done: " << nAntileptonsToBeDone << "\n"
+    << "#Open-beauty hadron pairs: " << nSignalPairs << "\n"
+    << "#leptons in acceptance: " << nLeptonsInAcceptance << "\n"
+    << "#Electron-positron pairs " << nLeptonPairs << "\n"
+    << "#Electron-positron pairs to be done: " << nLeptonPairsToBeDone << "\n";
+  if (nLeptons == 0 && nAntileptons == 0) {
+    std::cerr << "Number of leptons, number of anti-leptons should all be greater than 1.\n";
+    return 1;
+  }
+  if (nLeptonPairs != nLeptonPairsToBeDone) {
+    std::cerr << "The number of electron-positron pairs should be the same as the number of electron-positron pairs which should be transported.\n";
+    return 1;
+  }
+  if (nLeptons != nLeptonsToBeDone) {
+    std::cerr << "The number of leptons should be the same as the number of leptons which should be transported.\n";
+    return 1;
+  }
+  if (nLeptonsInAcceptance ==  (nEvents/ratioTrigger)) {
+    std::cerr << "The number of leptons in acceptance should be at least equaled to the number of events.\n";
+    return 1;
+  }
+  if (nLeptonPairs <  nLeptonsInAcceptance) {
+    std::cerr << "The number of positron-electron pairs should be at least equaled to the number of leptons in acceptance.\n";
+    return 1;
+  }
+
+  return 0;
+}
diff --git a/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_BeautyNoForcedDecay_Gap2_OO_electron.C b/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_BeautyNoForcedDecay_Gap2_OO_electron.C
new file mode 100644
index 000000000..ac2a627cf
--- /dev/null
+++ b/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_BeautyNoForcedDecay_Gap2_OO_electron.C
@@ -0,0 +1,112 @@
+int External()
+{
+  int checkPdgDecay = 11;
+  std::string path{"o2sim_Kine.root"};
+  TFile file(path.c_str(), "READ");
+  if (file.IsZombie()) {
+    std::cerr << "Cannot open ROOT file " << path << "\n";
+    return 1;
+  }
+
+  float ratioTrigger = 1./3; // one event triggered out of 3
+  auto tree = (TTree*)file.Get("o2sim");
+  std::vector<o2::MCTrack>* tracks{};
+  tree->SetBranchAddress("MCTrack", &tracks);
+
+  int nLeptonsInAcceptance{};
+  int nLeptons{};
+  int nAntileptons{};
+  int nLeptonsToBeDone{};
+  int nAntileptonsToBeDone{};
+  int nSignalPairs{};
+  int nLeptonPairs{};
+  int nLeptonPairsToBeDone{};
+  auto nEvents = tree->GetEntries();
+
+  for (int i = 0; i < nEvents; i++) {
+    tree->GetEntry(i);
+    int nElectrons = 0;
+    int nPositrons = 0;
+    int nElectronsToBeDone = 0;
+    int nPositronsToBeDone = 0;
+    int nOpenBeautyPos = 0;
+    int nOpenBeautyNeg = 0;
+    int nPositronsElectronsInAcceptance = 0;
+    for (auto& track : *tracks) {
+      auto pdg = track.GetPdgCode();
+      auto y = track.GetRapidity();
+      if (pdg == checkPdgDecay) {
+        int igmother = track.getMotherTrackId();
+        if (igmother > 0) {
+          auto gmTrack = (*tracks)[igmother];
+          int gmpdg = gmTrack.GetPdgCode();
+          if (int(std::abs(gmpdg)/100.) == 5 || int(std::abs(gmpdg)/1000.) == 5 || int(std::abs(gmpdg)/100.) == 4 || int(std::abs(gmpdg)/1000.) == 4) {
+            nLeptons++;
+            nElectrons++;
+            if (-1 < y && y < 1) nPositronsElectronsInAcceptance++;
+            if (track.getToBeDone()) {
+              nLeptonsToBeDone++;
+              nElectronsToBeDone++;
+            }
+          }
+        }    
+      } else if (pdg == -checkPdgDecay) {
+        int igmother = track.getMotherTrackId();
+        if (igmother > 0) {
+          auto gmTrack = (*tracks)[igmother];
+          int gmpdg = gmTrack.GetPdgCode();
+          if (int(TMath::Abs(gmpdg)/100.) == 4 || int(TMath::Abs(gmpdg)/1000.) == 4 || int(std::abs(gmpdg)/100.) == 5 || int(std::abs(gmpdg)/1000.) == 5) {
+            nAntileptons++;
+            nPositrons++;
+            if (-1 < y && y < 1) nPositronsElectronsInAcceptance++;
+            if (track.getToBeDone()) {
+              nAntileptonsToBeDone++;
+              nPositronsToBeDone++;
+            }
+          }
+        }
+      } else if (pdg == 511 || pdg == 521 || pdg == 531 || pdg == 5122 || pdg == 5132 || pdg == 5232 || pdg == 5332) {
+        nOpenBeautyPos++;
+      }  else if (pdg == -511 || pdg == -521 || pdg == -531 || pdg == -5122 || pdg == -5132 || pdg == -5232 || pdg == -5332) {
+        nOpenBeautyNeg++;
+      }
+    }
+    if (nOpenBeautyPos > 0 && nOpenBeautyNeg > 0) {
+      nSignalPairs++;
+    }
+    if (nPositronsElectronsInAcceptance > 1) {
+      nLeptonsInAcceptance++;
+    }
+    if (nElectrons > 0 && nPositrons > 0) {
+      nLeptonPairs++;
+    }
+    if (nElectronsToBeDone > 0 && nPositronsToBeDone > 0) nLeptonPairsToBeDone++;
+  }
+  std::cout << "#events: " << nEvents << "\n"
+    << "#leptons: " << nLeptons << "\n"
+    << "#antileptons: " << nAntileptons << "\n"
+    << "#leptons to be done: " << nLeptonsToBeDone << "\n"
+    << "#antileptons to be done: " << nAntileptonsToBeDone << "\n"
+    << "#Open-beauty hadron pairs: " << nSignalPairs << "\n"
+    << "#leptons in acceptance: " << nLeptonsInAcceptance << "\n"
+    << "#Electron-positron pairs: " << nLeptonPairs << "\n"
+    << "#Electron-positron pairs to be done: " << nLeptonPairsToBeDone << "\n";
+  if (nLeptons == 0 && nAntileptons == 0) {
+    std::cerr << "Number of leptons, number of anti-leptons 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;
+  }
+  if (nLeptons != nLeptonsToBeDone) {
+    std::cerr << "The number of leptons should be the same as the number of leptons which should be transported.\n";
+    return 1;
+  }
+  if (nLeptonsInAcceptance ==  (nEvents/ratioTrigger)) {
+    std::cerr << "The number of leptons in acceptance should be at least equaled to the number of events.\n";
+    return 1;
+  }
+
+  return 0;
+}
diff --git a/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_BeautyNoForcedDecay_Gap2_OO_muon.C b/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_BeautyNoForcedDecay_Gap2_OO_muon.C
new file mode 100644
index 000000000..ac2a627cf
--- /dev/null
+++ b/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_BeautyNoForcedDecay_Gap2_OO_muon.C
@@ -0,0 +1,112 @@
+int External()
+{
+  int checkPdgDecay = 11;
+  std::string path{"o2sim_Kine.root"};
+  TFile file(path.c_str(), "READ");
+  if (file.IsZombie()) {
+    std::cerr << "Cannot open ROOT file " << path << "\n";
+    return 1;
+  }
+
+  float ratioTrigger = 1./3; // one event triggered out of 3
+  auto tree = (TTree*)file.Get("o2sim");
+  std::vector<o2::MCTrack>* tracks{};
+  tree->SetBranchAddress("MCTrack", &tracks);
+
+  int nLeptonsInAcceptance{};
+  int nLeptons{};
+  int nAntileptons{};
+  int nLeptonsToBeDone{};
+  int nAntileptonsToBeDone{};
+  int nSignalPairs{};
+  int nLeptonPairs{};
+  int nLeptonPairsToBeDone{};
+  auto nEvents = tree->GetEntries();
+
+  for (int i = 0; i < nEvents; i++) {
+    tree->GetEntry(i);
+    int nElectrons = 0;
+    int nPositrons = 0;
+    int nElectronsToBeDone = 0;
+    int nPositronsToBeDone = 0;
+    int nOpenBeautyPos = 0;
+    int nOpenBeautyNeg = 0;
+    int nPositronsElectronsInAcceptance = 0;
+    for (auto& track : *tracks) {
+      auto pdg = track.GetPdgCode();
+      auto y = track.GetRapidity();
+      if (pdg == checkPdgDecay) {
+        int igmother = track.getMotherTrackId();
+        if (igmother > 0) {
+          auto gmTrack = (*tracks)[igmother];
+          int gmpdg = gmTrack.GetPdgCode();
+          if (int(std::abs(gmpdg)/100.) == 5 || int(std::abs(gmpdg)/1000.) == 5 || int(std::abs(gmpdg)/100.) == 4 || int(std::abs(gmpdg)/1000.) == 4) {
+            nLeptons++;
+            nElectrons++;
+            if (-1 < y && y < 1) nPositronsElectronsInAcceptance++;
+            if (track.getToBeDone()) {
+              nLeptonsToBeDone++;
+              nElectronsToBeDone++;
+            }
+          }
+        }    
+      } else if (pdg == -checkPdgDecay) {
+        int igmother = track.getMotherTrackId();
+        if (igmother > 0) {
+          auto gmTrack = (*tracks)[igmother];
+          int gmpdg = gmTrack.GetPdgCode();
+          if (int(TMath::Abs(gmpdg)/100.) == 4 || int(TMath::Abs(gmpdg)/1000.) == 4 || int(std::abs(gmpdg)/100.) == 5 || int(std::abs(gmpdg)/1000.) == 5) {
+            nAntileptons++;
+            nPositrons++;
+            if (-1 < y && y < 1) nPositronsElectronsInAcceptance++;
+            if (track.getToBeDone()) {
+              nAntileptonsToBeDone++;
+              nPositronsToBeDone++;
+            }
+          }
+        }
+      } else if (pdg == 511 || pdg == 521 || pdg == 531 || pdg == 5122 || pdg == 5132 || pdg == 5232 || pdg == 5332) {
+        nOpenBeautyPos++;
+      }  else if (pdg == -511 || pdg == -521 || pdg == -531 || pdg == -5122 || pdg == -5132 || pdg == -5232 || pdg == -5332) {
+        nOpenBeautyNeg++;
+      }
+    }
+    if (nOpenBeautyPos > 0 && nOpenBeautyNeg > 0) {
+      nSignalPairs++;
+    }
+    if (nPositronsElectronsInAcceptance > 1) {
+      nLeptonsInAcceptance++;
+    }
+    if (nElectrons > 0 && nPositrons > 0) {
+      nLeptonPairs++;
+    }
+    if (nElectronsToBeDone > 0 && nPositronsToBeDone > 0) nLeptonPairsToBeDone++;
+  }
+  std::cout << "#events: " << nEvents << "\n"
+    << "#leptons: " << nLeptons << "\n"
+    << "#antileptons: " << nAntileptons << "\n"
+    << "#leptons to be done: " << nLeptonsToBeDone << "\n"
+    << "#antileptons to be done: " << nAntileptonsToBeDone << "\n"
+    << "#Open-beauty hadron pairs: " << nSignalPairs << "\n"
+    << "#leptons in acceptance: " << nLeptonsInAcceptance << "\n"
+    << "#Electron-positron pairs: " << nLeptonPairs << "\n"
+    << "#Electron-positron pairs to be done: " << nLeptonPairsToBeDone << "\n";
+  if (nLeptons == 0 && nAntileptons == 0) {
+    std::cerr << "Number of leptons, number of anti-leptons 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;
+  }
+  if (nLeptons != nLeptonsToBeDone) {
+    std::cerr << "The number of leptons should be the same as the number of leptons which should be transported.\n";
+    return 1;
+  }
+  if (nLeptonsInAcceptance ==  (nEvents/ratioTrigger)) {
+    std::cerr << "The number of leptons in acceptance should be at least equaled to the number of events.\n";
+    return 1;
+  }
+
+  return 0;
+}
diff --git a/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_Charm_Gap2_OO_electron.C b/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_Charm_Gap2_OO_electron.C
new file mode 100644
index 000000000..ed19c336e
--- /dev/null
+++ b/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_Charm_Gap2_OO_electron.C
@@ -0,0 +1,117 @@
+int External()
+{
+
+  int checkPdgDecay = 11;
+  std::string path{"o2sim_Kine.root"};
+  TFile file(path.c_str(), "READ");
+  if (file.IsZombie()) {
+    std::cerr << "Cannot open ROOT file " << path << "\n";
+    return 1;
+  }
+
+  float ratioTrigger = 1./3.; // one event triggered out of 3
+  auto tree = (TTree*)file.Get("o2sim");
+  std::vector<o2::MCTrack>* tracks{};
+  tree->SetBranchAddress("MCTrack", &tracks);
+
+  int nLeptonsInAcceptance{};
+  int nLeptons{};
+  int nAntileptons{};
+  int nLeptonsToBeDone{};
+  int nAntileptonsToBeDone{};
+  int nSignalPairs{};
+  int nLeptonPairs{};
+  int nLeptonPairsToBeDone{};
+  auto nEvents = tree->GetEntries();
+
+  for (int i = 0; i < nEvents; i++) {
+    tree->GetEntry(i);
+    int nElectrons = 0;
+    int nPositrons = 0;
+    int nElectronsToBeDone = 0;
+    int nPositronsToBeDone = 0;
+    int nOpenCharmPos = 0;
+    int nOpenCharmNeg = 0;
+    int nPositronsElectronsInAcceptance = 0;
+    for (auto& track : *tracks) {
+      auto pdg = track.GetPdgCode();
+      auto y = track.GetRapidity();
+      if (pdg == checkPdgDecay) {
+        int igmother = track.getMotherTrackId();
+        if (igmother > 0) {
+          auto gmTrack = (*tracks)[igmother];
+          int gmpdg = gmTrack.GetPdgCode();
+          if (int(std::abs(gmpdg)/100.) == 4 || int(std::abs(gmpdg)/1000.) == 4) {
+            nLeptons++;
+            nElectrons++;
+            if (-1 < y && y < 1) nPositronsElectronsInAcceptance++;
+            if (track.getToBeDone()) {
+              nLeptonsToBeDone++;
+              nElectronsToBeDone++;
+            }
+          }
+        }    
+      } else if (pdg == -checkPdgDecay) {
+        int igmother = track.getMotherTrackId();
+        if (igmother > 0) {
+          auto gmTrack = (*tracks)[igmother];
+          int gmpdg = gmTrack.GetPdgCode();
+          if (int(TMath::Abs(gmpdg)/100.) == 4 || int(TMath::Abs(gmpdg)/1000.) == 4) {
+            nAntileptons++;
+            nPositrons++;
+            if (-1 < y && y < 1) nPositronsElectronsInAcceptance++;
+            if (track.getToBeDone()) {
+              nAntileptonsToBeDone++;
+              nPositronsToBeDone++;
+            }
+          }
+        }
+      } else if (pdg == 411 || pdg == 421 || pdg == 431 || pdg == 4122 || pdg == 4132 || pdg == 4232 || pdg == 4332) {
+        nOpenCharmPos++;
+      }  else if (pdg == -411 || pdg == -421 || pdg == -431 || pdg == -4122 || pdg == -4132 || pdg == -4232 || pdg == -4332) {
+        nOpenCharmNeg++;
+      }
+    }
+    if (nOpenCharmPos > 0 && nOpenCharmNeg > 0) {
+      nSignalPairs++;
+    }
+    if (nPositronsElectronsInAcceptance > 1) {
+      nLeptonsInAcceptance++;
+    }
+    if (nElectrons > 0 && nPositrons > 0) {
+      nLeptonPairs++;
+    }
+    if (nElectronsToBeDone > 0 && nPositronsToBeDone > 0) nLeptonPairsToBeDone++;
+  }
+  std::cout << "#events: " << nEvents << "\n"
+    << "#leptons: " << nLeptons << "\n"
+    << "#antileptons: " << nAntileptons << "\n"
+    << "#leptons to be done: " << nLeptonsToBeDone << "\n"
+    << "#antileptons to be done: " << nAntileptonsToBeDone << "\n"
+    << "#open-charm hadron pairs: " << nSignalPairs << "\n"
+    << "#leptons in acceptance: " << nLeptonsInAcceptance << "\n"
+    << "#Electron-positron pairs: " << nLeptonPairs << "\n"
+    << "#Electron-positron pairs to be done: " << nLeptonPairsToBeDone << "\n";
+  if (nLeptons == 0 && nAntileptons == 0) {
+    std::cerr << "Number of leptons, number of anti-leptons should all be greater than 1.\n";
+    return 1;
+  }
+  if (nLeptonPairs != nLeptonPairsToBeDone) {
+    std::cerr << "The number of electron-positron pairs should be the same as the number of electron-positron pairs which should be transported.\n";
+    return 1;
+  }
+  if (nLeptons != nLeptonsToBeDone) {
+    std::cerr << "The number of leptons should be the same as the number of leptons which should be transported.\n";
+    return 1;
+  }
+  if (nLeptonsInAcceptance ==  (nEvents/ratioTrigger)) {
+    std::cerr << "The number of leptons in acceptance should be at least equaled to the number of events.\n";
+    return 1;
+  }
+  if (nLeptonPairs <  nLeptonsInAcceptance) {
+    std::cerr << "The number of positron-electron pairs should be at least equaled to the number of leptons in acceptance.\n";
+    return 1;
+  }
+
+  return 0;
+}
diff --git a/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_Charm_Gap2_OO_muon.C b/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_Charm_Gap2_OO_muon.C
new file mode 100644
index 000000000..0b837b413
--- /dev/null
+++ b/MC/config/PWGEM/ini/tests/GeneratorHFGapTriggered_Charm_Gap2_OO_muon.C
@@ -0,0 +1,116 @@
+int External()
+{
+  int checkPdgDecay = 13;
+  std::string path{"o2sim_Kine.root"};
+  TFile file(path.c_str(), "READ");
+  if (file.IsZombie()) {
+    std::cerr << "Cannot open ROOT file " << path << "\n";
+    return 1;
+  }
+
+  float ratioTrigger = 1./3.; // one event triggered out of 3
+  auto tree = (TTree*)file.Get("o2sim");
+  std::vector<o2::MCTrack>* tracks{};
+  tree->SetBranchAddress("MCTrack", &tracks);
+
+  int nLeptonsInAcceptance{};
+  int nLeptons{};
+  int nAntileptons{};
+  int nLeptonsToBeDone{};
+  int nAntileptonsToBeDone{};
+  int nSignalPairs{};
+  int nLeptonPairs{};
+  int nLeptonPairsToBeDone{};
+  auto nEvents = tree->GetEntries();
+
+  for (int i = 0; i < nEvents; i++) {
+    tree->GetEntry(i);
+    int nElectrons = 0;
+    int nPositrons = 0;
+    int nElectronsToBeDone = 0;
+    int nPositronsToBeDone = 0;
+    int nOpenCharmPos = 0;
+    int nOpenCharmNeg = 0;
+    int nPositronsElectronsInAcceptance = 0;
+    for (auto& track : *tracks) {
+      auto pdg = track.GetPdgCode();
+      auto y = track.GetRapidity();
+      if (pdg == checkPdgDecay) {
+        int igmother = track.getMotherTrackId();
+        if (igmother > 0) {
+          auto gmTrack = (*tracks)[igmother];
+          int gmpdg = gmTrack.GetPdgCode();
+          if (int(std::abs(gmpdg)/100.) == 4 || int(std::abs(gmpdg)/1000.) == 4) {
+            nLeptons++;
+            nElectrons++;
+            if (-1 < y && y < 1) nPositronsElectronsInAcceptance++;
+            if (track.getToBeDone()) {
+              nLeptonsToBeDone++;
+              nElectronsToBeDone++;
+            }
+          }
+        }    
+      } else if (pdg == -checkPdgDecay) {
+        int igmother = track.getMotherTrackId();
+        if (igmother > 0) {
+          auto gmTrack = (*tracks)[igmother];
+          int gmpdg = gmTrack.GetPdgCode();
+          if (int(TMath::Abs(gmpdg)/100.) == 4 || int(TMath::Abs(gmpdg)/1000.) == 4) {
+            nAntileptons++;
+            nPositrons++;
+            if (-1 < y && y < 1) nPositronsElectronsInAcceptance++;
+            if (track.getToBeDone()) {
+              nAntileptonsToBeDone++;
+              nPositronsToBeDone++;
+            }
+          }
+        }
+      } else if (pdg == 411 || pdg == 421 || pdg == 431 || pdg == 4122 || pdg == 4132 || pdg == 4232 || pdg == 4332) {
+        nOpenCharmPos++;
+      }  else if (pdg == -411 || pdg == -421 || pdg == -431 || pdg == -4122 || pdg == -4132 || pdg == -4232 || pdg == -4332) {
+        nOpenCharmNeg++;
+      }
+    }
+    if (nOpenCharmPos > 0 && nOpenCharmNeg > 0) {
+      nSignalPairs++;
+    }
+    if (nPositronsElectronsInAcceptance > 1) {
+      nLeptonsInAcceptance++;
+    }
+    if (nElectrons > 0 && nPositrons > 0) {
+      nLeptonPairs++;
+    }
+    if (nElectronsToBeDone > 0 && nPositronsToBeDone > 0) nLeptonPairsToBeDone++;
+  }
+  std::cout << "#events: " << nEvents << "\n"
+    << "#leptons: " << nLeptons << "\n"
+    << "#antileptons: " << nAntileptons << "\n"
+    << "#leptons to be done: " << nLeptonsToBeDone << "\n"
+    << "#antileptons to be done: " << nAntileptonsToBeDone << "\n"
+    << "#open-charm hadron pairs: " << nSignalPairs << "\n"
+    << "#leptons in acceptance: " << nLeptonsInAcceptance << "\n"
+    << "#Electron-positron pairs: " << nLeptonPairs << "\n"
+    << "#Electron-positron pairs to be done: " << nLeptonPairsToBeDone << "\n";
+  if (nLeptons == 0 && nAntileptons == 0) {
+    std::cerr << "Number of leptons, number of anti-leptons should all be greater than 1.\n";
+    return 1;
+  }
+  if (nLeptonPairs != nLeptonPairsToBeDone) {
+    std::cerr << "The number of electron-positron pairs should be the same as the number of electron-positron pairs which should be transported.\n";
+    return 1;
+  }
+  if (nLeptons != nLeptonsToBeDone) {
+    std::cerr << "The number of leptons should be the same as the number of leptons which should be transported.\n";
+    return 1;
+  }
+  if (nLeptonsInAcceptance ==  (nEvents/ratioTrigger)) {
+    std::cerr << "The number of leptons in acceptance should be at least equaled to the number of events.\n";
+    return 1;
+  }
+  if (nLeptonPairs <  nLeptonsInAcceptance) {
+    std::cerr << "The number of positron-electron pairs should be at least equaled to the number of leptons in acceptance.\n";
+    return 1;
+  }
+
+  return 0;
+}