Pythia8Interface/plugins/TopRecoilHook.h

0001 // TopRecoilHook.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2020 Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005
0006 // Includes a user hook that corrects emission in top decay for dipole
0007 // from gluon to W, to instead be from gluon to top.
0008
0009 // Important: the top mass shift analysis encoded here is very primitive,
0010 // does not perform well at all, and should not be taken seriously.
0011 // The important part is that you see how the different scenarios
0012 // should be set up to operate as intended.
0013
0014 #include "Pythia8/Pythia.h"
0015 namespace Pythia8 {
0016
0017   //==========================================================================
0018
0019   // Write own derived UserHooks class for modified emission in top decay.
0020
0021   class TopRecoilHook : public UserHooks {
0022   public:
0023     // Constructor.
0024     //  doTopRecoil : eikonal correction in GW dipole on/off when no MEC applied.
0025     //  useOldDipole  : in GW dipole, use partons before or after branching.
0026     //  doList        : diagnostic output; set false for production runs.
0027     TopRecoilHook(bool doTopRecoilIn = true, bool useOldDipoleIn = false, bool doListIn = false) {
0028       doTopRecoil = doTopRecoilIn;
0029       useOldDipole = useOldDipoleIn;
0030       doList = doListIn;
0031       // Constructor also creates some histograms for analysis inside User Hook.
0032       wtCorr = new Hist("corrective weight", 100, 0., 2.);
0033     }
0034
0035     // Destructor prints histogram.
0036     ~TopRecoilHook() override { delete wtCorr; }
0037
0038     // Initialise. Only use hook for simple showers with recoilToColoured = off.
0039     bool initAfterBeams() override {
0040       // Switch off if recoilToColoured = on.
0041       bool recoilToColoured = settingsPtr->flag("TimeShower:recoilToColoured");
0042       if (recoilToColoured)
0043         doTopRecoil = false;
0044       // Flag if W mass term is already accounted for (true) or not (false).
0045       recoilDeadCone = settingsPtr->flag("TimeShower:recoilDeadCone");
0046       // All ok.
0047       return true;
0048     }
0049
0050     // Allow a veto after an FSR emission
0051     bool canVetoFSREmission() override { return doTopRecoil; }
0052
0053     // Access the event after an FSR emission, specifically inside top decay.
0054     bool doVetoFSREmission(int sizeOld, const Event& event, int iSys, bool inResonance) override {
0055       // Check that we are inside a resonance decay.
0056       if (!inResonance)
0057         return false;
0058
0059       // Check that it is a top decay.
0060       int iTop = partonSystemsPtr->getInRes(iSys);
0061       if (iTop == 0 || event[iTop].idAbs() != 6)
0062         return false;
0063
0064       // Skip first emission, where ME corrections are already made.
0065       int sizeOut = partonSystemsPtr->sizeOut(iSys);
0066       if (sizeOut == 2)
0067         return false;
0068
0069       // Location of trial new particles: radiator, emitted, recoiler.
0070       int iRad = sizeOld;
0071       int iEmt = sizeOld + 1;
0072       int iRec = sizeOld + 2;
0073
0074       // The above partons are after emission;
0075       // alternatively use the ones before.
0076       if (useOldDipole) {
0077         iRad = event[iRad].mother1();
0078         iRec = event[iRec].mother1();
0079       }
0080
0081       // Check if newly emitted gluon matches (anti)top colour line.
0082       if (event[iEmt].id() != 21)
0083         return false;
0084       if (event[iTop].id() == 6) {
0085         if (event[iEmt].col() != event[iTop].col())
0086           return false;
0087       } else {
0088         if (event[iEmt].acol() != event[iTop].acol())
0089           return false;
0090       }
0091
0092       // Recoiler should now be a W, else something is wrong.
0093       if (event[iRec].idAbs() != 24) {
0094         return false;
0095       }
0096
0097       // Denominator: eikonal weight with W as recoiler.
0098       double pRadRec = event[iRad].p() * event[iRec].p();
0099       double pRadEmt = event[iRad].p() * event[iEmt].p();
0100       double pRecEmt = event[iRec].p() * event[iEmt].p();
0101       double wtW = 2. * pRadRec / (pRadEmt * pRecEmt) - pow2(event[iRad].m() / pRadEmt);
0102       // If recoilDeadCone = on, include W mass term in denominator.
0103       if (recoilDeadCone)
0104         wtW -= pow2(event[iRec].m() / pRecEmt);
0105
0106       // Numerator: eikonal weight with top as recoiler.
0107       double pRadTop = event[iRad].p() * event[iTop].p();
0108       double pTopEmt = event[iTop].p() * event[iEmt].p();
0109       double wtT =
0110           2. * pRadTop / (pRadEmt * pTopEmt) - pow2(event[iRad].m() / pRadEmt) - pow2(event[iTop].m() / pTopEmt);
0111
0112       // Histogram weight ratio.
0113       wtCorr->fill(wtT / wtW);
0114
0115       // List relevant properties.
0116       if (doList) {
0117         partonSystemsPtr->list();
0118         event.list();
0119       }
0120
0121       // Accept/reject emission. Smooth suppression or step function.
0122       return (wtT < wtW * rndmPtr->flat());
0123     }
0124
0125   private:
0126     // Options and Histograms.
0127     bool doTopRecoil, useOldDipole, doList, recoilDeadCone;
0128     Hist* wtCorr;
0129   };
0130
0131 }  // namespace Pythia8