Pythia8Interface/plugins/PowhegResHook.cc

0001 // Hook for setting shower scale in top and W resonances
0002 // for Powheg ttb_NLO_dec and b_bbar_4l processes
0003 // C++ port of algorithm by Jezo et. al. (arXiv:1607.04538, Appendix B.2)
0004
0005 #include "Pythia8/Pythia.h"
0006
0007 using namespace Pythia8;
0008
0009 #include "GeneratorInterface/Pythia8Interface/plugins/PowhegResHook.h"
0010
0011 double PowhegResHook::scaleResonance(const int iRes, const Event& event) {
0012   calcScales_ = settingsPtr->flag("POWHEGres:calcScales");
0013
0014   double scale = 0.;
0015
0016   int nDau = event[iRes].daughterList().size();
0017
0018   if (!calcScales_ or nDau == 0) {
0019     // No resonance found, set scale to high value
0020     // Pythia will shower any MC generated resonance unrestricted
0021     scale = 1e30;
0022   }
0023
0024   else if (nDau < 3) {
0025     // No radiating resonance found
0026     scale = 0.8;
0027   }
0028
0029   else if (abs(event[iRes].id()) == 6) {
0030     // Find top daughters
0031     int idw = -1, idb = -1, idg = -1;
0032
0033     for (int i = 0; i < nDau; i++) {
0034       int iDau = event[iRes].daughterList()[i];
0035       if (abs(event[iDau].id()) == 24)
0036         idw = iDau;
0037       if (abs(event[iDau].id()) == 5)
0038         idb = iDau;
0039       if (abs(event[iDau].id()) == 21)
0040         idg = iDau;
0041     }
0042
0043     // Get daughter 4-vectors in resonance frame
0044     Vec4 pw(event[idw].p());
0045     pw.bstback(event[iRes].p());
0046
0047     Vec4 pb(event[idb].p());
0048     pb.bstback(event[iRes].p());
0049
0050     Vec4 pg(event[idg].p());
0051     pg.bstback(event[iRes].p());
0052
0053     // Calculate scale
0054     scale = sqrt(2 * pg * pb * pg.e() / pb.e());
0055   }
0056
0057   else if (abs(event[iRes].id()) == 24) {
0058     // Find W daughters
0059     int idq = -1, ida = -1, idg = -1;
0060
0061     for (int i = 0; i < nDau; i++) {
0062       int iDau = event[iRes].daughterList()[i];
0063       if (event[iDau].id() == 21)
0064         idg = iDau;
0065       else if (event[iDau].id() > 0)
0066         idq = iDau;
0067       else if (event[iDau].id() < 0)
0068         ida = iDau;
0069     }
0070
0071     // Get daughter 4-vectors in resonance frame
0072     Vec4 pq(event[idq].p());
0073     pq.bstback(event[iRes].p());
0074
0075     Vec4 pa(event[ida].p());
0076     pa.bstback(event[iRes].p());
0077
0078     Vec4 pg(event[idg].p());
0079     pg.bstback(event[iRes].p());
0080
0081     // Calculate scale
0082     Vec4 pw = pq + pa + pg;
0083     double q2 = pw * pw;
0084     double csi = 2 * pg.e() / sqrt(q2);
0085     double yq = 1 - pg * pq / (pg.e() * pq.e());
0086     double ya = 1 - pg * pa / (pg.e() * pa.e());
0087
0088     scale = sqrt(min(1 - yq, 1 - ya) * pow2(csi) * q2 / 2);
0089   }
0090
0091   return scale;
0092 }