Merge pull request #29807 from intrepid42/BiasedTauDecayer_94X

Backport of #29567 (biased tau decayer) to CMSSW_9_4_X

GeneratorInterface/Pythia8Interface/interface/BiasedTauDecayer.h

@@ -0,0 +1,30 @@
+#include "Pythia8/ParticleDecays.h"
+#include "Pythia8/Pythia.h"
+
+//==========================================================================
+
+// Specialized decayer for resonance decays to taus to allowing biasing to
+// leptonic decays
+//
+class BiasedTauDecayer : public Pythia8::DecayHandler {
+public:
+  BiasedTauDecayer(Pythia8::Info* infoPtr,
+                   Pythia8::Settings* settingsPtr,
+                   Pythia8::ParticleData* particleDataPtr,
+                   Pythia8::Rndm* rndmPtr,
+                   Pythia8::Couplings* couplingsPtr);
+
+  bool decay(std::vector<int>& idProd,
+             std::vector<double>& mProd,
+             std::vector<Pythia8::Vec4>& pProd,
+             int iDec,
+             const Pythia8::Event& event) override;
+
+private:
+  Pythia8::TauDecays decayer;
+  bool filter_;
+  bool eMuDecays_;
+  std::vector<int> idProdSave;
+  std::vector<double> mProdSave;
+  std::vector<Pythia8::Vec4> pProdSave;
+};

GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc

@@ -36,6 +36,9 @@ using namespace Pythia8;
 #include "GeneratorInterface/Pythia8Interface/plugins/PowhegResHook.h"
 #include "GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h"
 
+//biased tau decayer
+#include "GeneratorInterface/Pythia8Interface/interface/BiasedTauDecayer.h"
+
 //decay filter hook
 #include "GeneratorInterface/Pythia8Interface/interface/ResonanceDecayFilterHook.h"
 
@@ -138,6 +141,9 @@ class Pythia8Hadronizer : public Py8InterfaceBase {
     std::auto_ptr<PowhegResHook> fPowhegResHook;
     std::auto_ptr<PowhegHooksBB4L> fPowhegHooksBB4L;
 
+    // biased tau decayer
+    std::unique_ptr<BiasedTauDecayer> fBiasedTauDecayer;
+
     //resonance decay filter hook
     std::auto_ptr<ResonanceDecayFilterHook> fResonanceDecayFilterHook;
 
@@ -422,6 +428,15 @@ bool Pythia8Hadronizer::initializeForInternalPartons()
     fMultiUserHook->addHook(fMergingHook.get());
   }
 
+  bool biasedTauDecayer = fMasterGen->settings.flag("BiasedTauDecayer:filter");
+  if (biasedTauDecayer) {
+    fBiasedTauDecayer.reset(new BiasedTauDecayer(&(fMasterGen->info), &(fMasterGen->settings),
+	&(fMasterGen->particleData), &(fMasterGen->rndm), &(fMasterGen->couplings) ) );
+    std::vector<int> handledParticles;
+    handledParticles.push_back(15);
+    fMasterGen->setDecayPtr(fBiasedTauDecayer.get(),handledParticles);
+  }
+
   bool resonanceDecayFilter = fMasterGen->settings.flag("ResonanceDecayFilter:filter");
   if (resonanceDecayFilter) {
     fResonanceDecayFilterHook.reset(new ResonanceDecayFilterHook);
@@ -553,6 +568,15 @@ bool Pythia8Hadronizer::initializeForExternalPartons()
     fMultiUserHook->addHook(fMergingHook.get());
   }
 
+  bool biasedTauDecayer = fMasterGen->settings.flag("BiasedTauDecayer:filter");
+  if (biasedTauDecayer) {
+    fBiasedTauDecayer.reset(new BiasedTauDecayer(&(fMasterGen->info), &(fMasterGen->settings),
+	&(fMasterGen->particleData), &(fMasterGen->rndm), &(fMasterGen->couplings) ) );
+    std::vector<int> handledParticles;
+    handledParticles.push_back(15);
+    fMasterGen->setDecayPtr(fBiasedTauDecayer.get(),handledParticles);
+  }
+
   bool resonanceDecayFilter = fMasterGen->settings.flag("ResonanceDecayFilter:filter");
   if (resonanceDecayFilter) {
     fResonanceDecayFilterHook.reset(new ResonanceDecayFilterHook);

GeneratorInterface/Pythia8Interface/src/BiasedTauDecayer.cc

@@ -0,0 +1,97 @@
+#include "Pythia8/Pythia.h"
+#include "GeneratorInterface/Pythia8Interface/interface/BiasedTauDecayer.h"
+using namespace Pythia8;
+
+//==========================================================================
+
+// Specialized decayer for resonance decays to taus to allowing biasing to
+// leptonic decays
+//
+
+BiasedTauDecayer::BiasedTauDecayer(
+    Info* infoPtr, Settings* settingsPtr, ParticleData* particleDataPtr, Rndm* rndmPtr, Couplings* couplingsPtr) {
+  decayer = TauDecays();
+  decayer.init(infoPtr, settingsPtr, particleDataPtr, rndmPtr, couplingsPtr);
+  filter_ = settingsPtr->flag("BiasedTauDecayer:filter");
+  eMuDecays_ = settingsPtr->flag("BiasedTauDecayer:eMuDecays");
+}
+
+bool BiasedTauDecayer::decay(
+    std::vector<int>& idProd, std::vector<double>& mProd, std::vector<Vec4>& pProd, int iDec, const Event& event) {
+  if (!filter_)
+    return false;
+  if (idProd[0] != 15 && idProd[0] != -15)
+    return false;
+  int iStart = event[iDec].iTopCopyId();
+  int iMom = event[iStart].mother1();
+  int idMom = event[iMom].idAbs();
+  if (idMom != 23 && idMom != 24 && idMom != 25)
+    return false;
+  int iDau1 = event[iMom].daughter1();
+  int iDau2 = event[iMom].daughter2();
+  int iBot1 = event[iDau1].iBotCopyId();
+  int iBot2 = event[iDau2].iBotCopyId();
+  int iDecSis = (iDec == iBot1) ? iBot2 : iBot1;
+  // Check if sister has been decayed
+  // Since taus decays are correlated, use one decay, store the other
+  bool notDecayed = event[iDecSis].status() > 0 ? true : false;
+  if (notDecayed) {
+    // bias for leptonic decays
+    bool hasLepton = (eMuDecays_) ? false : true;
+    Event decay;
+    int i1 = -1;
+    int i2 = -1;
+    while (!hasLepton) {
+      decay = event;
+      decayer.decay(iDec, decay);
+      // check for lepton in first decay
+      i1 = decay[iDec].daughter1();
+      i2 = decay[iDec].daughter2();
+      for (int i = i1; i < i2 + 1; ++i) {
+        if (decay[i].isLepton() && decay[i].isCharged()) {
+          hasLepton = true;
+          break;
+        }
+      }
+      if (hasLepton)
+        break;
+      // check for lepton in second decay
+      i1 = decay[iDecSis].daughter1();
+      i2 = decay[iDecSis].daughter2();
+      for (int i = i1; i < i2 + 1; ++i) {
+        if (decay[i].isLepton() && decay[i].isCharged()) {
+          hasLepton = true;
+          break;
+        }
+      }
+    }
+    // Return decay products
+    i1 = decay[iDec].daughter1();
+    i2 = decay[iDec].daughter2();
+    for (int i = i1; i < i2 + 1; ++i) {
+      idProd.push_back(decay[i].id());
+      mProd.push_back(decay[i].m());
+      pProd.push_back(decay[i].p());
+    }
+    // Store correlated decay products
+    i1 = decay[iDecSis].daughter1();
+    i2 = decay[iDecSis].daughter2();
+    idProdSave.clear();
+    mProdSave.clear();
+    pProdSave.clear();
+    for (int i = i1; i < i2 + 1; ++i) {
+      idProdSave.push_back(decay[i].id());
+      mProdSave.push_back(decay[i].m());
+      pProdSave.push_back(decay[i].p());
+    }
+  } else {
+    // Return stored decay products
+    for (size_t i = 0; i < idProdSave.size(); ++i) {
+      idProd.push_back(idProdSave[i]);
+      mProd.push_back(mProdSave[i]);
+      pProd.push_back(pProdSave[i]);
+    }
+  }
+
+  return true;
+}

GeneratorInterface/Pythia8Interface/src/Py8InterfaceBase.cc

@@ -52,6 +52,10 @@ bool Py8InterfaceBase::readSettings( int )
    if(!fMasterGen.get()) fMasterGen.reset(new Pythia);
    fDecayer.reset(new Pythia);
 
+   //add settings for resonance decay filter
+   fMasterGen->settings.addFlag("BiasedTauDecayer:filter", false);
+   fMasterGen->settings.addFlag("BiasedTauDecayer:eMuDecays", true);
+
    //add settings for resonance decay filter
    fMasterGen->settings.addFlag("ResonanceDecayFilter:filter",false);
    fMasterGen->settings.addFlag("ResonanceDecayFilter:exclusive",false);