PowhegResHook.cc

Go to the documentation of this file.

// 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;
}