Merge pull request #36509 from cericeci/SUEPpythia_dev_10_6_X

Backport of Factory for adding UserHooks to PythiaHadronizer and SUEP decay pythia Hook

GeneratorInterface/Pythia8Interface/interface/CustomHook.h

@@ -0,0 +1,8 @@
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "Pythia8/Pythia.h"
+#include "FWCore/PluginManager/interface/PluginFactory.h"
+
+// Automatic addition of user hooks to pythia without the need to edit the hadronizer/generator files
+typedef edmplugin::PluginFactory<Pythia8::UserHooks*(const edm::ParameterSet&)> CustomHookFactory;
+
+#define REGISTER_USERHOOK(type) DEFINE_EDM_PLUGIN(CustomHookFactory, type, #type)

GeneratorInterface/Pythia8Interface/interface/MultiUserHook.h

@@ -10,6 +10,7 @@ class MultiUserHook : public Pythia8::UserHooks {
   // Constructor and destructor.
   MultiUserHook() {}
 
+  std::vector<Pythia8::UserHooks*> hooks() const { return hooks_; }
   unsigned int nHooks() const { return hooks_.size(); }
   void addHook(Pythia8::UserHooks *hook) { hooks_.push_back(hook); }
 

GeneratorInterface/Pythia8Interface/interface/SuepShower.h

@@ -0,0 +1,61 @@
+#ifndef GeneratorInterface_Pythia8Interface_SuepShower_h
+#define GeneratorInterface_Pythia8Interface_SuepShower_h
+
+#include <vector>
+#include <utility>
+#include <cmath>
+#include <boost/math/tools/roots.hpp>
+#include <boost/bind/bind.hpp>
+#include "Pythia8/Pythia.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+class SuepShower {
+public:
+  // Constructor
+  SuepShower(double mass, double temperature, Pythia8::Rndm* rndmPtr);
+
+  // Empty destructor
+  ~SuepShower();
+
+  // Method to generate 4-momenta of dark mesons after the showering
+  std::vector<Pythia8::Vec4> generateShower(double energy);
+
+private:
+  // private variables
+  // Shower parameters
+  double darkmeson_mass_;
+  double mass_over_T_;
+
+  // Overall energy of the decaying particle
+  double mediator_energy_;
+
+  // For the numerical algorithm precision
+  boost::math::tools::eps_tolerance<double> tolerance_;
+
+  // Several auxiliar variables for generating the 4-momentum of showered particles. Following the naming of Appendix 1 of https://arxiv.org/pdf/1305.5226.pdf
+  // Median momentum in the M-B distribution
+  double p_m_;
+  // Two values of momentum at fMB(x)/fMB(p_m_) = e
+  double p_plus_, p_minus_;
+  // Auxiliars: fMB(p_plus_)/f'(p_plus_),  fMB(p_minus_)/f'(p_minus_)
+  double lambda_plus_, lambda_minus_;
+  // More auxiliars: lambda_plus_/(p_plus_ + p_minus_), lambda_minus_/(p_plus_ + p_minus_), 1-q_plus_-q_minus_
+  double q_plus_, q_minus_, q_m_;
+
+  // Pythia random number generator, to get the randomness into the shower
+  Pythia8::Rndm* fRndmPtr_;
+
+  // Methods
+  // Maxwell-Boltzmann distribution as a function of |p|, slightly massaged, Eq. 6 of https://arxiv.org/pdf/1305.5226.pdf
+  const double fMaxwellBoltzmann(double p);
+  // Maxwell-Boltzmann derivative as a function of |p|, slightly massaged
+  const double fMaxwellBoltzmannPrime(double p);
+  // Log(fMaxwellBoltzmann(x)/fMaxwellBoltzmann(xmedian))+1, as a function of |p|, to be solved to find p_plus_, p_minus_
+  const double logTestFunction(double p);
+  // Generate the four vector of a particle (dark meson) in the shower
+  const Pythia8::Vec4 generateFourVector();
+  // sum(scale*scale*p*p + m*m) - E*E, find roots in scale to reballance momenta and impose E conservation
+  const double reballanceFunction(double scale, const std::vector<Pythia8::Vec4>& shower);
+};
+
+#endif

GeneratorInterface/Pythia8Interface/plugins/BuildFile.xml

@@ -10,6 +10,8 @@
 </library>
 
 <library file="SLHA*.cc" name="GeneratorInterfacePythia8SLHAReader">
-  <use name="GeneratorInterface/Pythia8Interface"/>
   <use name="root"/>
 </library>
+
+<library file="Suep*.cc" name="GeneratorInterfacePythia8SUEPHook">
+</library>

GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc

@@ -19,6 +19,7 @@ using namespace Pythia8;
 #include "GeneratorInterface/Pythia8Interface/interface/Py8InterfaceBase.h"
 
 #include "GeneratorInterface/Pythia8Interface/plugins/ReweightUserHooks.h"
+#include "GeneratorInterface/Pythia8Interface/interface/CustomHook.h"
 
 // PS matchning prototype
 //
@@ -151,6 +152,9 @@ class Pythia8Hadronizer : public Py8InterfaceBase {
   //PT filter hook
   std::unique_ptr<PTFilterHook> fPTFilterHook;
 
+  //Generic customized hooks vector
+  std::unique_ptr<MultiUserHook> fCustomHooksVector;
+
   int EV1_nFinal;
   bool EV1_vetoOn;
   int EV1_maxVetoCount;
@@ -313,6 +317,16 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet &params)
                                                    0));
   }
 
+  fCustomHooksVector.reset(new MultiUserHook);
+  if (params.exists("UserCustomization")) {
+    const std::vector<edm::ParameterSet> userParams =
+        params.getParameter<std::vector<edm::ParameterSet>>("UserCustomization");
+    for (const auto &pluginParams : userParams) {
+      fCustomHooksVector->addHook(
+              CustomHookFactory::get()->create(pluginParams.getParameter<std::string>("pluginName"), pluginParams));
+    }
+  }
+
   if (params.exists("VinciaPlugin")) {
     fMasterGen.reset(new Pythia);
     fvincia.reset(new Vincia::VinciaPlugin(fMasterGen.get()));
@@ -446,6 +460,13 @@ bool Pythia8Hadronizer::initializeForInternalPartons() {
     fMultiUserHook->addHook(fPTFilterHook.get());
   }
 
+  if (fCustomHooksVector->nHooks() > 0) {
+    edm::LogInfo("Pythia8Interface") << "Adding customized user hooks";
+    for (const auto &fUserHook : fCustomHooksVector.get()->hooks()) {
+      fMultiUserHook->addHook(fUserHook);
+    }
+  }
+
   if (fMultiUserHook->nHooks() > 0) {
     fMasterGen->setUserHooksPtr(fMultiUserHook.get());
   }
@@ -514,6 +535,13 @@ bool Pythia8Hadronizer::initializeForExternalPartons() {
     fMultiUserHook->addHook(fEmissionVetoHook1.get());
   }
 
+  if (fCustomHooksVector->nHooks() > 0) {
+    edm::LogInfo("Pythia8Interface") << "Adding customized user hooks";
+    for (const auto &fUserHook : fCustomHooksVector.get()->hooks()) {
+      fMultiUserHook->addHook(fUserHook);
+    }
+  }
+
   if (fMasterGen->settings.mode("POWHEG:veto") > 0 || fMasterGen->settings.mode("POWHEG:MPIveto") > 0) {
     if (fJetMatchingHook.get() || fEmissionVetoHook1.get())
       throw edm::Exception(edm::errors::Configuration, "Pythia8Interface")

GeneratorInterface/Pythia8Interface/plugins/SuepDecay.cc

@@ -0,0 +1,57 @@
+#include "SuepDecay.h"
+
+SuepDecay::SuepDecay(const edm::ParameterSet& iConfig)
+    : idMediator_(iConfig.getParameter<int>("idMediator")),
+      idDark_(iConfig.getParameter<int>("idDark")),
+      temperature_(iConfig.getParameter<double>("temperature")) {}
+
+bool SuepDecay::initAfterBeams() {
+  mDark_ = particleDataPtr->m0(idDark_);
+  bool medDecay = particleDataPtr->mayDecay(idMediator_);
+  if (!medDecay) {
+    infoPtr->errorMsg("Error in SuepDecay::initAfterBeams: mediator decay should be enabled");
+    return false;
+  }
+
+  //construct the shower helper
+  suep_shower_ = std::make_unique<SuepShower>(mDark_, temperature_, rndmPtr);
+
+  return true;
+}
+
+//based on https://gitlab.com/simonknapen/suep_generator/-/blob/master/suep_main.cc:AttachSuepShower
+bool SuepDecay::doVetoProcessLevel(Pythia8::Event& event) {
+  Pythia8::Vec4 pMediator, pDark;
+
+  // Find the mediator in the event
+  for (int i = 0; i < event.size(); ++i) {
+    //mediator w/ distinct daughters = last copy (decayed)
+    if (event[i].id() == idMediator_ && event[i].daughter1() != event[i].daughter2() && event[i].daughter1() > 0 &&
+        event[i].daughter2() > 0) {
+      pMediator = event[i].p();
+      // undo mediator decay
+      event[i].undoDecay();
+
+      // Generate the shower, output are 4 vectors in the rest frame of the shower, adding energy here avoids issues if scalar is off-shell
+      std::vector<Pythia8::Vec4> suep_shower_fourmomenta = suep_shower_->generateShower(pMediator.mCalc());
+      // Loop over hidden sector mesons and append to the event
+      int firstDaughter = event.size();
+      for (auto& pDark : suep_shower_fourmomenta) {
+        // Boost to the lab frame, i.e. apply the mediator boost
+        pDark.bst(pMediator);
+        // Append particle to the event w/ hidden meson pdg code. Magic number 91 means it is produced as a normal decay product
+        event.append(idDark_, 91, i, 0, 0, 0, 0, 0, pDark.px(), pDark.py(), pDark.pz(), pDark.e(), mDark_);
+      }
+
+      // Change the status code of the mediator to reflect that it has decayed.
+      event[i].statusNeg();
+
+      //set daughters of the mediator: daughter1 < daughter2 > 0 -> the particle has a range of decay products from daughter1 to daughter2
+      event[i].daughters(firstDaughter, event.size() - 1);
+      break;
+    }
+  }
+
+  //allow event to continue
+  return false;
+}

GeneratorInterface/Pythia8Interface/plugins/SuepDecay.h

@@ -0,0 +1,29 @@
+#ifndef GeneratorInterface_Pythia8Interface_SuepDecay_h
+#define GeneratorInterface_Pythia8Interface_SuepDecay_h
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/PluginManager/interface/PluginFactory.h"
+#include "Pythia8/Pythia.h"
+#include "GeneratorInterface/Pythia8Interface/interface/SuepShower.h"
+#include <memory>
+#include "GeneratorInterface/Pythia8Interface/interface/CustomHook.h"
+
+// Adapted by Kevin Pedro to run on cmssw as a user hook
+class SuepDecay : public Pythia8::UserHooks {
+public:
+  SuepDecay(const edm::ParameterSet& iConfig);
+  ~SuepDecay() override {}
+
+  bool initAfterBeams() override;
+
+  bool canVetoProcessLevel() override { return true; }
+  bool doVetoProcessLevel(Pythia8::Event& event) override;
+
+protected:
+  int idMediator_, idDark_;
+  float temperature_, mDark_;
+  std::unique_ptr<SuepShower> suep_shower_;
+};
+
+REGISTER_USERHOOK(SuepDecay);
+#endif

GeneratorInterface/Pythia8Interface/src/CustomHookFactory.cc

@@ -0,0 +1,3 @@
+#include "GeneratorInterface/Pythia8Interface/interface/CustomHook.h"
+
+EDM_REGISTER_PLUGINFACTORY(CustomHookFactory, "CustomHookFactory");