Merge pull request #17251 from intrepid42/bb4l_71x

Powheg-RES integration

GeneratorInterface/Pythia8Interface/plugins/PowhegResHook.cc

@@ -0,0 +1,87 @@
+// Hook for setting shower scale in top and W resonances
+// for Powheg ttb_NLO_dec and b_bbar_4l processes
+// C++ port of algorithm by Jezo et. al. (arXiv:1607.04538, Appendix B.2)
+
+#include "Pythia8/Pythia.h"
+
+using namespace Pythia8;
+
+#include "GeneratorInterface/Pythia8Interface/plugins/PowhegResHook.h"
+
+double PowhegResHook::scaleResonance( const int iRes, const Event& event) {
+  calcScales_ = settingsPtr->flag("POWHEGres:calcScales");
+
+  double scale = 0.;
+
+  int nDau = event[iRes].daughterList().size();
+
+  if (!calcScales_ or nDau == 0) {
+    // No resonance found, set scale to high value
+    // Pythia will shower any MC generated resonance unrestricted
+    scale = 1e30;
+  }
+
+  else if (nDau < 3) {
+    // No radiating resonance found
+    scale = 0.8;
+  }
+
+  else if (abs(event[iRes].id()) == 6) {
+    // Find top daughters
+    int idw = -1, idb = -1, idg = -1;
+
+    for (int i = 0; i < nDau; i++) {
+      int iDau = event[iRes].daughterList()[i];
+      if (abs(event[iDau].id()) == 24) idw = iDau;
+      if (abs(event[iDau].id()) ==  5) idb = iDau;
+      if (abs(event[iDau].id()) == 21) idg = iDau;
+    }
+
+    // Get daughter 4-vectors in resonance frame
+    Vec4 pw(event[idw].p());
+    pw.bstback(event[iRes].p());
+
+    Vec4 pb(event[idb].p());
+    pb.bstback(event[iRes].p());
+
+    Vec4 pg(event[idg].p());
+    pg.bstback(event[iRes].p());
+
+    // Calculate scale
+    scale = sqrt(2*pg*pb*pg.e()/pb.e());
+  }
+
+  else if (abs(event[iRes].id()) == 24) {
+    // Find W daughters
+    int idq = -1, ida = -1, idg = -1;
+
+    for (int i = 0; i < nDau; i++) {
+      int iDau = event[iRes].daughterList()[i];
+      if      (event[iDau].id() == 21) idg = iDau;
+      else if (event[iDau].id()  >  0) idq = iDau;
+      else if (event[iDau].id()  <  0) ida = iDau;
+    }
+
+    // Get daughter 4-vectors in resonance frame
+    Vec4 pq(event[idq].p());
+    pq.bstback(event[iRes].p());
+
+    Vec4 pa(event[ida].p());
+    pa.bstback(event[iRes].p());
+
+    Vec4 pg(event[idg].p());
+    pg.bstback(event[iRes].p());
+
+    // Calculate scale
+    Vec4 pw = pq + pa + pg;
+    double q2 = pw*pw;
+    double csi = 2*pg.e()/sqrt(q2);
+    double yq = 1 - pg*pq/(pg.e()*pq.e());
+    double ya = 1 - pg*pa/(pg.e()*pa.e());
+
+    scale = sqrt(min(1-yq,1-ya)*pow2(csi)*q2/2);
+  }
+
+  return scale;
+}
+

GeneratorInterface/Pythia8Interface/plugins/PowhegResHook.h

@@ -0,0 +1,22 @@
+// Hook for setting shower scale in top and W resonances
+// for Powheg ttb_NLO_dec and b_bbar_4l processes
+// C++ port of algorithm by Jezo et. al. (arXiv:1607.04538, Appendix B.2)
+
+class PowhegResHook : public Pythia8::UserHooks {
+
+public:
+
+  // Constructor and destructor.
+  PowhegResHook() {}
+  ~PowhegResHook() {}
+
+  bool canSetResonanceScale() { return true; }
+
+  double scaleResonance( const int iRes, const Pythia8::Event& event);
+
+//--------------------------------------------------------------------------
+
+private:
+  bool calcScales_;
+
+};

GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc

@@ -30,6 +30,9 @@ using namespace Pythia8;
 #include "Pythia8Plugins/PowhegHooks.h"
 #include "GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.h"
 
+// Resonance scale hook
+#include "GeneratorInterface/Pythia8Interface/plugins/PowhegResHook.h"
+
 //decay filter hook
 #include "GeneratorInterface/Pythia8Interface/interface/ResonanceDecayFilterHook.h"
 
@@ -124,6 +127,9 @@ class Pythia8Hadronizer : public BaseHadronizer, public Py8InterfaceBase {
     PowhegHooks* fEmissionVetoHook;
     EmissionVetoHook1* fEmissionVetoHook1;
 
+    // Resonance scale hook
+    PowhegResHook* fPowhegResHook;
+
     //resonance decay filter hook
     ResonanceDecayFilterHook *fResonanceDecayFilterHook;
 
@@ -258,6 +264,9 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet &params) :
   fMasterGen->settings.addParm("PTFilter:scaleToFilter", 0.4,true,true,0.0, 10.);
   fMasterGen->settings.addParm("PTFilter:quarkRapidity",10.0,true,true,0.0, 10.);
   fMasterGen->settings.addParm("PTFilter:quarkPt",       -.1,true,true,-.1,100.);
+
+  //add settings for powheg resonance scale calculation
+  fMasterGen->settings.addFlag("POWHEGres:calcScales",false);
 
   // Reweight user hook
   //
@@ -401,6 +410,26 @@ bool Pythia8Hadronizer::initializeForInternalPartons()
       <<" UNKNOWN INITIAL STATE. \n The allowed initial states are: PP, PPbar, ElectronPositron \n";
   }
 
+  if((fMasterGen->settings.mode("POWHEG:veto") > 0 || fMasterGen->settings.mode("POWHEG:MPIveto") > 0) && !fEmissionVetoHook) {
+
+    if(fJetMatchingHook || fEmissionVetoHook1)
+      throw edm::Exception(edm::errors::Configuration,"Pythia8Interface")
+      <<" Attempt to turn on PowhegHooks by pythia8 settings but there are incompatible hooks on \n Incompatible are : jetMatching, emissionVeto1 \n";
+
+    fEmissionVetoHook = new PowhegHooks();
+
+    edm::LogInfo("Pythia8Interface") << "Turning on Emission Veto Hook from pythia8 code";
+    fMultiUserHook->addHook(fEmissionVetoHook);
+
+  }
+
+  bool PowhegRes = fMasterGen->settings.flag("POWHEGres:calcScales");
+  if (PowhegRes) {
+    edm::LogInfo("Pythia8Interface") << "Turning on resonance scale setting from CMSSW Pythia8Interface";
+    fPowhegResHook = new PowhegResHook();
+    fMultiUserHook->addHook(fPowhegResHook);
+  }
+
   bool resonanceDecayFilter = fMasterGen->settings.flag("ResonanceDecayFilter:filter");
   if (resonanceDecayFilter && !fResonanceDecayFilterHook) {
     fResonanceDecayFilterHook = new ResonanceDecayFilterHook;
@@ -465,6 +494,13 @@ bool Pythia8Hadronizer::initializeForExternalPartons()
     fMultiUserHook->addHook(fEmissionVetoHook);
 
   }
+
+  bool PowhegRes = fMasterGen->settings.flag("POWHEGres:calcScales");
+  if (PowhegRes) {
+    edm::LogInfo("Pythia8Interface") << "Turning on resonance scale setting from CMSSW Pythia8Interface";
+    fPowhegResHook = new PowhegResHook();
+    fMultiUserHook->addHook(fPowhegResHook);
+  }
 
   //adapted from main89.cc in pythia8 examples
   bool internalMatching = fMasterGen->settings.flag("JetMatching:merge");
@@ -589,7 +625,21 @@ bool Pythia8Hadronizer::generatePartonsAndHadronize()
   if (evtgenDecays) evtgenDecays->decay();
 
   event().reset(new HepMC::GenEvent);
-  return toHepMC.fill_next_event( *(fMasterGen.get()), event().get());
+  bool py8hepmc =  toHepMC.fill_next_event( *(fMasterGen.get()), event().get());
+
+  if (!py8hepmc) {
+    return false;
+  }
+
+  //fill compresse lhe weights for systematic uncertainties
+  if (fMasterGen->info.getWeightsCompressedSize() > 0) {
+    for (unsigned int i = 0; i < fMasterGen->info.getWeightsCompressedSize(); i++) {
+      double wgt = fMasterGen->info.getWeightsCompressedValue(i);
+      event()->weights().push_back(wgt);
+    }
+  }
+
+  return true;
 
 }