From 23d963772f25a9d1b75f1a74d11ba9a2b11b631a Mon Sep 17 00:00:00 2001 From: Laurids Jeppe Date: Thu, 22 Feb 2024 17:02:37 +0100 Subject: [PATCH 1/6] Replaced old bb4l hook by a version provided by the authors --- .../plugins/PowhegHooksBB4L.h | 985 +++++++++--------- .../Pythia8Interface/src/Py8InterfaceBase.cc | 6 +- 2 files changed, 481 insertions(+), 510 deletions(-) diff --git a/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h b/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h index 411273ed0d612..854b7077a408a 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h +++ b/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h @@ -1,505 +1,480 @@ -// PowhegHooksBB4L.h -// Copyright (C) 2017 Silvia Ferrario Ravasio, Tomas Jezo, Paolo Nason, Markus Seidel -// inspired by PowhegHooks.h by Richard Corke -// adjusted to work with EmissionVetoHook1 in CMSSW by Alexander Grohsjean - -#ifndef Pythia8_PowhegHooksBB4L_H -#define Pythia8_PowhegHooksBB4L_H - -// Includes -#include "Pythia8/Pythia.h" -#include -struct { - int radtype; -} radtype_; - -namespace Pythia8 { - - class PowhegHooksBB4L : public UserHooks { - public: - //--- Constructor and destructor ------------------------------------------- - PowhegHooksBB4L() : nFSRvetoBB4l(0) {} - ~PowhegHooksBB4L() override { std::cout << "Number of FSR vetoed in BB4l = " << nFSRvetoBB4l << std::endl; } - - //--- Initialization ----------------------------------------------------------------------- - bool initAfterBeams() override { - // settings of this class - vetoFSREmission = settingsPtr->flag("POWHEG:bb4l:FSREmission:veto"); - onlyDistance1 = settingsPtr->flag("POWHEG:bb4l:FSREmission:onlyDistance1"); - dryRunFSR = settingsPtr->flag("POWHEG:bb4l:FSREmission:dryRun"); - vetoAtPL = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoAtPL"); - vetoQED = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoQED"); - vetoPartonLevel = settingsPtr->flag("POWHEG:bb4l:PartonLevel:veto"); - excludeFSRConflicting = settingsPtr->flag("POWHEG:bb4l:PartonLevel:excludeFSRConflicting"); - debug = settingsPtr->flag("POWHEG:bb4l:DEBUG"); - scaleResonanceVeto = settingsPtr->flag("POWHEG:bb4l:ScaleResonance:veto"); - vetoDipoleFrame = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoDipoleFrame"); - pTpythiaVeto = settingsPtr->flag("POWHEG:bb4l:FSREmission:pTpythiaVeto"); - //vetoProduction = (settingsPtr->mode("POWHEG:veto")==1); - pTmin = settingsPtr->parm("POWHEG:bb4l:pTminVeto"); - return true; - } - - //--- PROCESS LEVEL HOOK --------------------------------------------------- - - // called at the LHE level - inline bool canVetoProcessLevel() override { return true; } - inline bool doVetoProcessLevel(Event &e) override { - // extract the radtype from the event comment - stringstream ss; - // use eventattribute as comments not filled when using edm input - //ss << infoPtr->getEventComments(); - ss << infoPtr->getEventAttribute("#rwgt"); - string temp; - ss >> temp >> radtype_.radtype; - assert(temp == "#rwgt"); - - // find last top and the last anti-top in the record - int i_top = -1, i_atop = -1; - for (int i = 0; i < e.size(); i++) { - if (e[i].id() == 6) - i_top = i; - if (e[i].id() == -6) - i_atop = i; - } - if (i_top != -1) - topresscale = findresscale(i_top, e); - else - topresscale = 1e30; - if (i_top != -1) - atopresscale = findresscale(i_atop, e); - else - atopresscale = 1e30; - // initialize stuff - doVetoFSRInit(); - // do not veto, ever - return false; - } - - //--- PARTON LEVEL HOOK ---------------------------------------------------- - - // called after shower - bool retryPartonLevel() override { return vetoPartonLevel || vetoAtPL; } - inline bool canVetoPartonLevel() override { return vetoPartonLevel || vetoAtPL; } - inline bool doVetoPartonLevel(const Event &e) override { - if (radtype_.radtype == 2) - return false; - if (debug) { - if (dryRunFSR && wouldVetoFsr) { - double scale = getdechardness(vetoTopCharge, e); - cout << "FSRdecScale = " << vetoDecScale << ", PLdecScale = " << scale << ", ratio = " << vetoDecScale / scale - << endl; - } - } - if (vetoPartonLevel) { - double topdecscale = getdechardness(1, e); - double atopdecscale = getdechardness(-1, e); - if ((topdecscale > topresscale) || (atopdecscale > atopresscale)) { - //if(dryRunFSR && ! wouldVetoFsr) mydatacontainer_.excludeEvent = excludeFSRConflicting?1:0; - return true; - } else - //if(dryRunFSR && wouldVetoFsr) mydatacontainer_.excludeEvent = excludeFSRConflicting?1:0; - return false; - } - if (vetoAtPL) { - if (dryRunFSR && wouldVetoFsr) - return true; - else - return false; - } - return false; - } - - //--- FSR EMISSION LEVEL HOOK ---------------------------------------------- - - // FSR veto: this should be true if we want PowhegHooksBB4l veto in decay - // OR PowhegHooks veto in production. (The virtual method - // PowhegHooks::canVetoFSREmission has been replaced by - // PowhegHooksBB4L::canVetoFSREmission, so FSR veto in production - // must be handled here. ISR and MPI veto are instead still - // handled by PowhegHooks.) - inline bool canVetoFSREmission() override { return vetoFSREmission; } // || vetoProduction; } - inline bool doVetoFSREmission(int sizeOld, const Event &e, int iSys, bool inResonance) override { - ////////////////////////////// - //VETO INSIDE THE RESONANCE // - ////////////////////////////// - if (inResonance && vetoFSREmission) { - int iRecAft = e.size() - 1; - int iEmt = e.size() - 2; - int iRadAft = e.size() - 3; - int iRadBef = e[iEmt].mother1(); - - // find the top resonance the radiator originates from - int iTop = e[iRadBef].mother1(); - int distance = 1; - while (abs(e[iTop].id()) != 6 && iTop > 0) { - iTop = e[iTop].mother1(); - distance++; - } - if (iTop == 0) { - infoPtr->errorMsg( - "Warning in PowhegHooksBB4L::doVetoFSREmission: emission in resonance not from top quark, not vetoing"); - return doVetoFSR(false, 0, 0); - //return false; - } - int iTopCharge = (e[iTop].id() > 0) ? 1 : -1; - - // calculate the scale of the emission - double scale; - //using pythia pT definition ... - if (pTpythiaVeto) - scale = pTpythia(e, iRadAft, iEmt, iRecAft); - //.. or using POWHEG pT definition - else { - Vec4 pr(e[iRadAft].p()), pe(e[iEmt].p()), pt(e[iTop].p()), prec(e[iRecAft].p()), psystem; - // The computation of the POWHEG pT can be done in the top rest frame or in the diple one. - // pdipole = pemt +prec +prad (after the emission) - // For the first emission off the top resonance pdipole = pw +pb (before the emission) = ptop - if (vetoDipoleFrame) - psystem = pr + pe + prec; - else - psystem = pt; - - // gluon splitting into two partons - if (e[iRadBef].id() == 21) - scale = gSplittingScale(psystem, pr, pe); - // quark emitting a gluon (or a photon) - else if (abs(e[iRadBef].id()) == 5 && ((e[iEmt].id() == 21) && !vetoQED)) - scale = qSplittingScale(psystem, pr, pe); - // other stuff (which we should not veto) - else - scale = 0; - } - - if (iTopCharge > 0) { - if (onlyDistance1) { - if (debug && (distance == 1) && scale > topresscale && !wouldVetoFsr) - cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() - << "; " << scale << endl; - return doVetoFSR((distance == 1) && scale > topresscale, scale, iTopCharge); - } else { - if (debug && scale > topresscale && !wouldVetoFsr) - cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() - << "; " << scale << endl; - return doVetoFSR(scale > topresscale, scale, iTopCharge); - } - } else if (iTopCharge < 0) { - if (onlyDistance1) { - if (debug && (distance == 1) && scale > atopresscale && !wouldVetoFsr) - cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() - << "; " << scale << endl; - return doVetoFSR((distance == 1) && scale > atopresscale, scale, iTopCharge); - } else { - if (debug && scale > topresscale && !wouldVetoFsr) - cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() - << "; " << scale << endl; - return doVetoFSR(scale > atopresscale, scale, iTopCharge); - } - } else { - cout << "Bug in PohwgeHooksBB4l" << endl; - } - } - ///////////////////////////////// - // VETO THE PRODUCTION PROCESS // - ///////////////////////////////// - // covered by multiuserhook, i.e. need to turn on EV1 - // else if(!inResonance && vetoProduction){ - // return EmissionVetoHook1::doVetoFSREmission(sizeOld, e, iSys, inResonance); - // } - - return false; - } - - inline bool doVetoFSR(bool condition, double scale, int iTopCharge) { - if (radtype_.radtype == 2) - return false; - if (condition) { - if (!wouldVetoFsr) { - wouldVetoFsr = true; - vetoDecScale = scale; - vetoTopCharge = iTopCharge; - } - if (dryRunFSR) - return false; - else { - nFSRvetoBB4l++; - return true; - } - } else - return false; - } - - inline void doVetoFSRInit() { - wouldVetoFsr = false; - vetoDecScale = -1; - vetoTopCharge = 0; - } - - //--- SCALE RESONANCE HOOK ------------------------------------------------- - // called before each resonance decay shower - inline bool canSetResonanceScale() override { return scaleResonanceVeto; } - // if the resonance is the (anti)top set the scale to: - // ---> (anti)top virtuality if radtype=2 - // ---> (a)topresscale otherwise - // if is not the top, set it to a big number - inline double scaleResonance(int iRes, const Event &e) override { - if (e[iRes].id() == 6) { - if (radtype_.radtype == 2) - return sqrt(e[iRes].m2Calc()); - else - return topresscale; - } else if (e[iRes].id() == -6) { - if (radtype_.radtype == 2) - return sqrt(e[iRes].m2Calc()); - else - return atopresscale; - } else - return pow(10.0, 30.); - } - - //--- Internal helper functions -------------------------------------------- - - // Calculates the scale of the hardest emission from within the resonance system - // translated by Markus Seidel modified by Tomas Jezo - inline double findresscale(const int iRes, const Event &event) { - double scale = 0.; - - int nDau = event[iRes].daughterList().size(); - - if (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 = pTmin; - } 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 { - scale = 1e30; - } - - return scale; - } - - // The following routine will match daughters of particle `e[iparticle]` to an expected pattern specified via the list of expected particle PDG ID's `ids`, - // id wildcard is specified as 0 if match is obtained, the positions and the momenta of these particles are returned in vectors `positions` and `momenta` - // respectively - // if exitOnExtraLegs==true, it will exit if the decay has more particles than expected, but not less - inline bool match_decay(int iparticle, - const Event &e, - const vector &ids, - vector &positions, - vector &momenta, - bool exitOnExtraLegs = true) { - // compare sizes - if (e[iparticle].daughterList().size() != ids.size()) { - if (exitOnExtraLegs && e[iparticle].daughterList().size() > ids.size()) { - cout << "extra leg" << endl; - exit(-1); - } - return false; - } - // compare content - for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) { - int di = e[iparticle].daughterList()[i]; - if (ids[i] != 0 && e[di].id() != ids[i]) - return false; - } - // reset the positions and momenta vectors (because they may be reused) - positions.clear(); - momenta.clear(); - // construct the array of momenta - for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) { - int di = e[iparticle].daughterList()[i]; - positions.push_back(di); - momenta.push_back(e[di].p()); - } - return true; - } - - inline double qSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2) { - p1.bstback(pt); - p2.bstback(pt); - return sqrt(2 * p1 * p2 * p2.e() / p1.e()); - } - - inline double gSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2) { - p1.bstback(pt); - p2.bstback(pt); - return sqrt(2 * p1 * p2 * p1.e() * p2.e() / (pow(p1.e() + p2.e(), 2))); - } - - // Routines to calculate the pT (according to pTdefMode) in a FS splitting: - // i (radiator before) -> j (emitted after) k (radiator after) - // For the Pythia pT definition, a recoiler (after) must be specified. - // (INSPIRED BY pythia8F77_31.cc double pTpythia) - inline double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch) { - // Convenient shorthands for later - Vec4 radVec = e[RadAfterBranch].p(); - Vec4 emtVec = e[EmtAfterBranch].p(); - Vec4 recVec = e[RecAfterBranch].p(); - int radID = e[RadAfterBranch].id(); - - // Calculate virtuality of splitting - Vec4 Q(radVec + emtVec); - double Qsq = Q.m2Calc(); - - // Mass term of radiator - double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ? pow2(particleDataPtr->m0(radID)) : 0.; - - // z values for FSR - double z, pTnow; - // Construct 2 -> 3 variables - Vec4 sum = radVec + recVec + emtVec; - double m2Dip = sum.m2Calc(); - - double x1 = 2. * (sum * radVec) / m2Dip; - double x3 = 2. * (sum * emtVec) / m2Dip; - z = x1 / (x1 + x3); - pTnow = z * (1. - z); - - // Virtuality - pTnow *= (Qsq - m2Rad); - - if (pTnow < 0.) { - cout << "Warning: pTpythia was negative" << endl; - return -1.; - } else - return (sqrt(pTnow)); - } - - inline double getdechardness(int topcharge, const Event &e) { - int tid = 6 * topcharge, wid = 24 * topcharge, bid = 5 * topcharge, gid = 21, wildcard = 0; - // find last top in the record - int i_top = -1; - Vec4 p_top, p_b, p_g, p_g1, p_g2; - for (int i = 0; i < e.size(); i++) - if (e[i].id() == tid) { - i_top = i; - p_top = e[i].p(); - } - if (i_top == -1) - return -1.0; - - // summary of cases - // 1.) t > W b - // a.) b > 3 ... error - // b.) b > b g ... h = sqrt(2*p_g*p_b*p_g.e()/p_b.e()) - // c.) b > other ... h = -1 - // return h - // 2.) t > W b g - // a.) b > 3 ... error - // b.) b > b g ... h1 = sqrt(2*p_g*p_b*p_g.e()/p_b.e()) - // c.) b > other ... h1 = -1 - // i.) g > 3 ... error - // ii.) g > 2 ... h2 = sqrt(2*p_g1*p_g2*p_g1.e()*p_g2.e()/(pow(p_g1.e(),2)+pow(p_g2.e(),2))) ); - // iii.) g > other ... h2 = -1 - // return max(h1,h2) - // 3.) else ... error - - vector momenta; - vector positions; - - // 1.) t > b W - if (match_decay(i_top, e, vector{wid, bid}, positions, momenta, false)) { - double h; - int i_b = positions[1]; - // a.+b.) b > 3 or b > b g - if (match_decay(i_b, e, vector{bid, gid}, positions, momenta)) - h = qSplittingScale(e[i_top].p(), momenta[0], momenta[1]); - // c.) b > other - else - h = -1; - return h; - } - // 2.) t > b W g - else if (match_decay(i_top, e, vector{wid, bid, gid}, positions, momenta, false)) { - double h1, h2; - int i_b = positions[1], i_g = positions[2]; - // a.+b.) b > 3 or b > b g - if (match_decay(i_b, e, vector{bid, gid}, positions, momenta)) - h1 = qSplittingScale(e[i_top].p(), momenta[0], momenta[1]); - // c.) b > other - else - h1 = -1; - // i.+ii.) g > 3 or g > 2 - if (match_decay(i_g, e, vector{wildcard, wildcard}, positions, momenta)) - h2 = gSplittingScale(e[i_top].p(), momenta[0], momenta[1]); - // c.) b > other - else - h2 = -1; - return max(h1, h2); - } - // 3.) else - else { - cout << "getdechardness" << endl; - cout << "top at position " << i_top << endl; - cout << "with " << e[i_top].daughterList().size() << " daughters " << endl; - for (unsigned i = 0; i < e[i_top].daughterList().size(); i++) { - int di = e[i_top].daughterList()[i]; - cout << "with daughter " << di << ": " << e[di].id() << endl; - } - exit(-1); - } - } - - //-------------------------------------------------------------------------- - - // Functions to return information - - // inline int getNFSRveto() { return nFSRveto; } - - //-------------------------------------------------------------------------- - - private: - // FSR emission veto flags - bool vetoFSREmission, dryRunFSR, wouldVetoFsr, onlyDistance1, vetoAtPL, vetoQED; - // Parton Level veto flags - bool vetoPartonLevel, excludeFSRConflicting; - // Scale Resonance veto flags - double scaleResonanceVeto; - // other flags - bool debug; - // internal: resonance scales - double topresscale, atopresscale; - // internal: inter veto communication - double vetoDecScale; - int vetoTopCharge; - bool vetoDipoleFrame; - bool pTpythiaVeto; - //bool vetoProduction; - double pTmin; - // Statistics on vetos - unsigned long int nFSRvetoBB4l; - }; - - //========================================================================== - -} // end namespace Pythia8 - -#endif // end Pythia8_PowhegHooksBB4L_H +// PowhegHooksBB4L.h +// Rewritten by T. Jezo in 2021. With various contributions from S. Ferrario +// Ravasio, B. Nachman, P. Nason and M. Seidel. Inspired by +// ttb_NLO_dec/main-PYTHIA8.f by P. Nason and E. Re and by PowhegHooks.h by R. +// Corke. +// +// Adapted for CMSSW by Laurids Jeppe. + +// # Introduction +// +// This hook is intended for use together with POWHEG-BOX-RES/b_bbar_4l NLO LHE +// events. This also includes events in which one of the W bosons was +// re-decayed hadronically. (Note that LHE format version larger than 3 may not +// require this hook). +// +// The hook inherits from PowhegHooks and as such it indirectly implements the +// following: +// - doVetoMPIStep +// - doVetoISREmission +// - doVetoMPIEmission +// and it overloads: +// - doVetoFSREmission, which works as follows (if POWHEG:veto==1): +// - if inResonance==true it vetoes all the emission that is harder than +// the scale of its parent (anti-)top quark or W^+(-) +// - if inResonance==false, it calls PowhegHooks::doVetoISREmission +// and it also implements: +// - doVetoProcessLevel, which is never used for vetoing (i.e. it always +// returns false). Instead it is used for the calculation of reconance scales +// using LHE kinematics. +// +// This version of the hooks is only suitable for use with fully compatible +// POWHEG-BOX Les Houches readers (such the one in main-PYTHIA82-lhef but +// not the one in main31.cc.) +// +// +// # Basic use +// +// In order to use this hook you must replace all the declarations and +// constructor calls of PowhegHooks to PowhegHooksBB4L: +// +// PowhegHooks *powhegHooks; -> PowhegHooksBB4L *powhegHooks; +// *powhegHooks = new PowhegHooks(); -> *powhegHooks = new PowhegHooksBB4L(); +// +// In order to switch it on set POWHEG:veto = 1 and +// POWHEG:bb4l:FSREmission:veto = 1. This will lead to a veto in ISR, FSR and +// MPI steps of pythia as expected using PowhegHooks in all the cases other than +// the case of FSR emission in resonance decay. Within resonance decays +// PowhegHooksBB4L takes over the control. +// +// Furthermore, this hook can also be used standalone without PowhegHooks, i.e. +// only the FSR emission from resonances will be vetoed (the default use in +// 1801.03944 and 1906.09166). In order to do that set +// POWHEG:bb4l:FSREmission:veto = 1 and POWHEG:veto = 0. Note that the this is +// not recommended for comparing against data but because it is easier to +// interpret it is often done in theoretical studies. +// +// Note that this version of the hook relies on the event "radtype" (1 for +// btilde, 2 for remnant) to be set by an external program, such as +// main-PYTHIA82-lhef in the radtype_ common block. +// There also exists a version of this hook in which the event "radtype" is +// read in from the .lhe file using pythia built in functionality. You need +// that version if you want to use this hook with main31.cc. +// +// +// # Expert use +// +// This hook also implements an alternative veto procedure which allows to +// assign a "SCALUP" type of scale to a resonance using the scaleResonance +// method. This is a much simpler veto but it is also clearly inferior as +// compared to the one implemented using the doVetoFSREmission method because +// the definition of the scale of the emission does not match the +// POWHEG-BOX-RES definition. Nevertheless, it can be activated using +// POWHEG:bb4l:ScaleResonance:veto = 1. Additionally one MUST switch off the +// other veto by calling on POWHEG:bb4l:FSREmission:veto = 0. +// +// The following settings are at the disposal of the user to control the +// behaviour of the hook +// - On/off switches for the veto: +// - POWHEG:bb4l:FSREmission:veto +// on/off switch for the default veto based on doFSREmission +// - POWHEG:bb4l:ScaleResonance:veto +// on/off switch for the alternative veto based on scaleResonance (only +// for expert users) +// - Important settings: +// - POWHEG:bb4l:ptMinVeto: MUST be set to the same value as the +// corresponding flag in POWHEG-BOX-RES +// - Alternatives for scale calculations +// - default: emission scale is calculated using POWHEG definitions and in +// the resonance rest frame +// - POWHEG:bb4l:FSREmission:vetoDipoleFrame: emission scale is calculated +// using POWHEG definitions in the dipole frame +// - POWHEG:bb4l:FSREmission:pTpythiaVeto: emission scale is calculated +// using Pythia definitions +// - Other flags: +// - POWHEG:bb4l:FSREmission:vetoQED: decides whether or not QED emission +// off quarks should also be vetoed (not implemented in the case of +// the ScaleResonance:veto) +// - POWHEG:bb4l:DEBUG: enables debug printouts on standard output + +#ifndef Pythia8_PowhegHooksBB4L_H +#define Pythia8_PowhegHooksBB4L_H + +#include "Pythia8/Pythia.h" +#include + +namespace Pythia8 { + +class PowhegHooksBB4L : public UserHooks { + + public: + + PowhegHooksBB4L() {} + ~PowhegHooksBB4L() { + std::cout << "Number of FSR vetoed in BB4l = " << nInResonanceFSRveto << std::endl; + } + + //--- Initialization ------------------------------------------------------- + bool initAfterBeams() { + // initialize settings of this class + vetoFSREmission = settingsPtr->flag("POWHEG:bb4l:FSREmission:veto"); + vetoDipoleFrame = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoDipoleFrame"); + pTpythiaVeto = settingsPtr->flag("POWHEG:bb4l:FSREmission:pTpythiaVeto"); + vetoQED = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoQED"); + scaleResonanceVeto = settingsPtr->flag("POWHEG:bb4l:ScaleResonance:veto"); + debug = settingsPtr->flag("POWHEG:bb4l:DEBUG"); + pTmin = settingsPtr->parm("POWHEG:bb4l:pTminVeto"); + nInResonanceFSRveto = 0; + return true; + } + + //--- PROCESS LEVEL HOOK --------------------------------------------------- + // This hook gets triggered for each event before the shower starts, i.e. at + // the LHE level. We use it to calculate the scales of resonances. + inline bool canVetoProcessLevel() { return true; } + inline bool doVetoProcessLevel(Event &e) { + // extract the radtype from the event comment + stringstream ss; + // use eventattribute as comments not filled when using edm input + ss << infoPtr->getEventComments(); + string temp; + ss >> temp >> radtype; + assert(temp == "#rwgt"); + // we only calculate resonance scales for btilde events (radtype==1) + // remnant events are not vetoed + if (radtype==2) return false; + // find last top and the last anti-top in the record + int i_top = -1, i_atop = -1, i_wp = -1, i_wm = -1; + for (int i = 0; i < e.size(); i++) { + if (e[i].id() == 6) i_top = i; + if (e[i].id() == -6) i_atop = i; + if (e[i].id() == 24) i_wp = i; + if (e[i].id() == -24) i_wm = i; + } + // if found calculate the resonance scale + topresscale = findresscale(i_top, e); + // similarly for anti-top + atopresscale = findresscale(i_atop, e); + // and for W^+ and W^- + wpresscale = findresscale(i_wp, e); + wmresscale = findresscale(i_wm, e); + + // do not veto, ever + return false; + } + + //--- FSR EMISSION LEVEL HOOK ---------------------------------------------- + // This hook gets triggered everytime the parton shower attempts to attach + // a FSR emission. + inline bool canVetoFSREmission() { return vetoFSREmission; } + inline bool doVetoFSREmission(int sizeOld, const Event &e, int iSys, bool inResonance) { + + // FSR VETO INSIDE THE RESONANCE (if it is switched on) + if (inResonance && vetoFSREmission) { + + // get the participants of the splitting: the recoiler, the radiator and the emitted + int iRecAft = e.size() - 1; + int iEmt = e.size() - 2; + int iRadAft = e.size() - 3; + int iRadBef = e[iEmt].mother1(); + + // find the resonance the radiator originates from + int iRes = e[iRadBef].mother1(); + int distance = 1; + while ( iRes > 0 && (abs(e[iRes].id()) !=6 && abs(e[iRes].id()) != 24) ) { + iRes = e[iRes].mother1(); + distance ++; + } + if (iRes == 0) { + infoPtr->errorMsg("Warning in PowhegHooksBB4L::doVetoFSREmission: emission in resonance not from the top quark or from the W boson, not vetoing"); + return doVetoFSR(false,0); + } + int iResId = e[iRes].id(); + + // calculate the scale of the emission + double scale; + //using pythia pT definition ... + if(pTpythiaVeto) + scale = pTpythia(e, iRadAft, iEmt, iRecAft); + //.. or using POWHEG pT definition + else{ + Vec4 pr(e[iRadAft].p()), pe(e[iEmt].p()), pres(e[iRes].p()), prec(e[iRecAft].p()), psystem; + // The computation of the POWHEG pT can be done in the top rest frame or in the diple one. + // pdipole = pemt +prec +prad (after the emission) + // For the first emission off the top resonance pdipole = pw +pb (before the emission) = ptop + if(vetoDipoleFrame) + psystem = pr+pe+prec; + else + psystem = pres; + + // gluon splitting into two partons + if (e[iRadBef].id() == 21) + scale = gSplittingScale(psystem, pr, pe); + // quark emitting a gluon (or a photon) + else if (abs(e[iRadBef].id()) <= 5 && ((e[iEmt].id() == 21) && ! vetoQED) ) + scale = qSplittingScale(psystem, pr, pe); + // other stuff (which we should not veto) + else { + scale = 0; + } + } + + // compare the current splitting scale to the correct resonance scale + if (iResId == 6) { + if ( debug && scale > topresscale) + cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; + return doVetoFSR(scale > topresscale, scale); + } + else if (iResId == -6){ + if ( debug && scale > atopresscale) + cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; + return doVetoFSR(scale > atopresscale, scale); + } + else if (iResId == 24) { + if ( debug && scale > wpresscale) + cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; + return doVetoFSR(scale > wpresscale, scale); + } + else if (iResId == -24){ + if ( debug && scale > wmresscale) + cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; + return doVetoFSR(scale > wmresscale, scale); + } + else { + infoPtr->errorMsg("Error in PowhegHooksBB4L::doVetoFSREmission: unimplemented case"); + exit(-1); + } + } + // In CMSSW, the production process veto is done in EmissionVetoHook1.cc + // so for events outside resonance, nothing needs to be done here + else { + return false; + } + } + + inline bool doVetoFSR(bool condition, double scale) { + if (radtype==2) return false; + if (condition) { + nInResonanceFSRveto++; + return true; + } + return false; + } + + //--- SCALE RESONANCE HOOK ------------------------------------------------- + // called before each resonance decay shower + inline bool canSetResonanceScale() { return scaleResonanceVeto; } + // if the resonance is the (anti)top or W+/W- set the scale to: + // - if radtype=2 (remnant): resonance virtuality + // - if radtype=1 (btilde): + // - (a)topresscale/wp(m)resscale for tops and Ws + // - a large number otherwise + // if is not the top, set it to a big number + inline double scaleResonance(int iRes, const Event &e) { + if(radtype == 2) + return sqrt(e[iRes].m2Calc()); + else { + if (e[iRes].id() == 6) + return topresscale; + else if (e[iRes].id() == -6) + return atopresscale; + else if (e[iRes].id() == 24) + return wpresscale; + else if (e[iRes].id() == 24) + return wmresscale; + else + return 1e30; + } + } + + //--- Internal helper functions -------------------------------------------- + // Calculates the scale of the hardest emission from within the resonance system + // translated by Markus Seidel modified by Tomas Jezo + inline double findresscale( const int iRes, const Event& event) { + + // return large scale if the resonance position is ill defined + if (iRes < 0) return 1e30; + + // get number of resonance decay products + int nDau = event[iRes].daughterList().size(); + + // iRes is not decayed, return high scale equivalent to + // unrestricted shower + if (nDau == 0) { + return 1e30; + } + // iRes did not radiate, this means that POWHEG pt scale has + // evolved all the way down to pTmin + else if (nDau < 3) { + return pTmin; + } + // iRes is a (anti-)top quark + 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 and return it + return sqrt(2*pg*pb*pg.e()/pb.e()); + } + // iRes is a W+(-) boson + 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()); + // and return it + return sqrt(min(1-yq,1-ya)*pow2(csi)*q2/2); + } + // in any other case just return a high scale equivalent to + // unrestricted shower + return 1e30; + } + + // The following routine will match daughters of particle `e[iparticle]` to an expected pattern specified via the list of expected particle PDG ID's `ids`, + // id wildcard is specified as 0 if match is obtained, the positions and the momenta of these particles are returned in vectors `positions` and `momenta` + // respectively + // if exitOnExtraLegs==true, it will exit if the decay has more particles than expected, but not less + inline bool match_decay(int iparticle, const Event &e, const vector &ids, vector &positions, vector &momenta, bool exitOnExtraLegs = true){ + // compare sizes + if (e[iparticle].daughterList().size() != ids.size()) { + if (exitOnExtraLegs && e[iparticle].daughterList().size() > ids.size()) { + cout << "extra leg" << endl; + exit(-1); + } + return false; + } + // compare content + for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) { + int di = e[iparticle].daughterList()[i]; + if (ids[i] != 0 && e[di].id() != ids[i]) + return false; + } + // reset the positions and momenta vectors (because they may be reused) + positions.clear(); + momenta.clear(); + // construct the array of momenta + for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) { + int di = e[iparticle].daughterList()[i]; + positions.push_back(di); + momenta.push_back(e[di].p()); + } + return true; + } + + inline double qSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){ + p1.bstback(pt); + p2.bstback(pt); + return sqrt( 2*p1*p2*p2.e()/p1.e() ); + } + + inline double gSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){ + p1.bstback(pt); + p2.bstback(pt); + return sqrt( 2*p1*p2*p1.e()*p2.e()/(pow(p1.e()+p2.e(),2)) ); + } + + // Routines to calculate the pT (according to pTdefMode) in a FS splitting: + // i (radiator before) -> j (emitted after) k (radiator after) + // For the Pythia pT definition, a recoiler (after) must be specified. + // (INSPIRED BY pythia8F77_31.cc double pTpythia) + inline double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch, + int RecAfterBranch) + { + + // Convenient shorthands for later + Vec4 radVec = e[RadAfterBranch].p(); + Vec4 emtVec = e[EmtAfterBranch].p(); + Vec4 recVec = e[RecAfterBranch].p(); + int radID = e[RadAfterBranch].id(); + + // Calculate virtuality of splitting + Vec4 Q(radVec + emtVec); + double Qsq = Q.m2Calc(); + + // Mass term of radiator + double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ? + pow2(particleDataPtr->m0(radID)) : 0.; + + // z values for FSR + double z, pTnow; + // Construct 2 -> 3 variables + Vec4 sum = radVec + recVec + emtVec; + double m2Dip = sum.m2Calc(); + + double x1 = 2. * (sum * radVec) / m2Dip; + double x3 = 2. * (sum * emtVec) / m2Dip; + z = x1 / (x1 + x3); + pTnow = z * (1. - z); + + + // Virtuality + pTnow *= (Qsq - m2Rad); + + if (pTnow < 0.) { + cout << "Warning: pTpythia was negative" << endl; + return -1.; + } + else + return(sqrt(pTnow)); + } + + // Functions to return statistics about the veto + inline int getNInResonanceFSRVeto() { return nInResonanceFSRveto; } + + //-------------------------------------------------------------------------- + + private: + // FSR emission veto flags + bool vetoFSREmission, vetoQED; + // scale Resonance veto flags + double scaleResonanceVeto; + // other flags + bool debug; + bool vetoDipoleFrame; + bool pTpythiaVeto; + double pTmin; + // veto counter + int nInResonanceFSRveto; + // internal: resonance scales + double topresscale, atopresscale, wpresscale, wmresscale; + int radtype; +}; + +//========================================================================== + +} // end namespace Pythia8 + +#endif // end Pythia8_PowhegHooksBB4L_H diff --git a/GeneratorInterface/Pythia8Interface/src/Py8InterfaceBase.cc b/GeneratorInterface/Pythia8Interface/src/Py8InterfaceBase.cc index 99f8ec98dfcd1..885a69a427c33 100644 --- a/GeneratorInterface/Pythia8Interface/src/Py8InterfaceBase.cc +++ b/GeneratorInterface/Pythia8Interface/src/Py8InterfaceBase.cc @@ -112,18 +112,14 @@ namespace gen { //add settings for powheg resonance scale calculation fMasterGen->settings.addFlag("POWHEGres:calcScales", false); fMasterGen->settings.addFlag("POWHEG:bb4l", false); - fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:onlyDistance1", false); fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:veto", false); - fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:dryRun", false); - fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:vetoAtPL", false); fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:vetoQED", false); - fMasterGen->settings.addFlag("POWHEG:bb4l:PartonLevel:veto", false); - fMasterGen->settings.addFlag("POWHEG:bb4l:PartonLevel:excludeFSRConflicting", false); fMasterGen->settings.addFlag("POWHEG:bb4l:DEBUG", false); fMasterGen->settings.addFlag("POWHEG:bb4l:ScaleResonance:veto", false); fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:vetoDipoleFrame", false); fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:pTpythiaVeto", false); fMasterGen->settings.addParm("POWHEG:bb4l:pTminVeto", 10.0, true, true, 0.0, 10.); + fMasterGen->settings.addFlag("POWHEG:bb4l:vetoAllRadtypes", false); fMasterGen->setRndmEnginePtr(&p8RndmEngine_); fDecayer->setRndmEnginePtr(&p8RndmEngine_); From d14557a4158da5791964c3cd6d29cf74056b0088 Mon Sep 17 00:00:00 2001 From: Laurids Jeppe Date: Thu, 22 Feb 2024 17:20:17 +0100 Subject: [PATCH 2/6] Fix reading radtype --- GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h b/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h index 854b7077a408a..11e2129b11a09 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h +++ b/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h @@ -136,7 +136,7 @@ class PowhegHooksBB4L : public UserHooks { // extract the radtype from the event comment stringstream ss; // use eventattribute as comments not filled when using edm input - ss << infoPtr->getEventComments(); + ss << infoPtr->getEventAttribute("#rwgt"); string temp; ss >> temp >> radtype; assert(temp == "#rwgt"); From 4599a1e4f7cb31f525c1b4b612652ef5e5c34a2e Mon Sep 17 00:00:00 2001 From: Laurids Jeppe Date: Thu, 22 Feb 2024 17:29:46 +0100 Subject: [PATCH 3/6] Added vetoAllRadtypes flag for ttb_NLO_dec --- .../Pythia8Interface/plugins/PowhegHooksBB4L.h | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h b/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h index 11e2129b11a09..2c579493457ae 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h +++ b/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h @@ -124,6 +124,7 @@ class PowhegHooksBB4L : public UserHooks { scaleResonanceVeto = settingsPtr->flag("POWHEG:bb4l:ScaleResonance:veto"); debug = settingsPtr->flag("POWHEG:bb4l:DEBUG"); pTmin = settingsPtr->parm("POWHEG:bb4l:pTminVeto"); + vetoAllRadtypes = settingsPtr->flag("POWHEG:bb4l:vetoAllRadtypes"); nInResonanceFSRveto = 0; return true; } @@ -142,7 +143,7 @@ class PowhegHooksBB4L : public UserHooks { assert(temp == "#rwgt"); // we only calculate resonance scales for btilde events (radtype==1) // remnant events are not vetoed - if (radtype==2) return false; + if (!vetoAllRadtypes && radtype==2) return false; // find last top and the last anti-top in the record int i_top = -1, i_atop = -1, i_wp = -1, i_wm = -1; for (int i = 0; i < e.size(); i++) { @@ -253,7 +254,7 @@ class PowhegHooksBB4L : public UserHooks { } inline bool doVetoFSR(bool condition, double scale) { - if (radtype==2) return false; + if (!vetoAllRadtypes && radtype==2) return false; if (condition) { nInResonanceFSRveto++; return true; @@ -271,7 +272,7 @@ class PowhegHooksBB4L : public UserHooks { // - a large number otherwise // if is not the top, set it to a big number inline double scaleResonance(int iRes, const Event &e) { - if(radtype == 2) + if(!vetoAllRadtypes && radtype == 2) return sqrt(e[iRes].m2Calc()); else { if (e[iRes].id() == 6) @@ -466,6 +467,7 @@ class PowhegHooksBB4L : public UserHooks { bool vetoDipoleFrame; bool pTpythiaVeto; double pTmin; + bool vetoAllRadtypes; // veto counter int nInResonanceFSRveto; // internal: resonance scales From 08e80c9dbcd4deb0daecac4855ff57553424056c Mon Sep 17 00:00:00 2001 From: Laurids Jeppe Date: Thu, 22 Feb 2024 17:38:13 +0100 Subject: [PATCH 4/6] Change indentation for comparibility --- .../plugins/PowhegHooksBB4L.h | 964 +++++++++--------- 1 file changed, 482 insertions(+), 482 deletions(-) diff --git a/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h b/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h index 2c579493457ae..75e6fc7111dac 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h +++ b/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h @@ -1,482 +1,482 @@ -// PowhegHooksBB4L.h -// Rewritten by T. Jezo in 2021. With various contributions from S. Ferrario -// Ravasio, B. Nachman, P. Nason and M. Seidel. Inspired by -// ttb_NLO_dec/main-PYTHIA8.f by P. Nason and E. Re and by PowhegHooks.h by R. -// Corke. -// -// Adapted for CMSSW by Laurids Jeppe. - -// # Introduction -// -// This hook is intended for use together with POWHEG-BOX-RES/b_bbar_4l NLO LHE -// events. This also includes events in which one of the W bosons was -// re-decayed hadronically. (Note that LHE format version larger than 3 may not -// require this hook). -// -// The hook inherits from PowhegHooks and as such it indirectly implements the -// following: -// - doVetoMPIStep -// - doVetoISREmission -// - doVetoMPIEmission -// and it overloads: -// - doVetoFSREmission, which works as follows (if POWHEG:veto==1): -// - if inResonance==true it vetoes all the emission that is harder than -// the scale of its parent (anti-)top quark or W^+(-) -// - if inResonance==false, it calls PowhegHooks::doVetoISREmission -// and it also implements: -// - doVetoProcessLevel, which is never used for vetoing (i.e. it always -// returns false). Instead it is used for the calculation of reconance scales -// using LHE kinematics. -// -// This version of the hooks is only suitable for use with fully compatible -// POWHEG-BOX Les Houches readers (such the one in main-PYTHIA82-lhef but -// not the one in main31.cc.) -// -// -// # Basic use -// -// In order to use this hook you must replace all the declarations and -// constructor calls of PowhegHooks to PowhegHooksBB4L: -// -// PowhegHooks *powhegHooks; -> PowhegHooksBB4L *powhegHooks; -// *powhegHooks = new PowhegHooks(); -> *powhegHooks = new PowhegHooksBB4L(); -// -// In order to switch it on set POWHEG:veto = 1 and -// POWHEG:bb4l:FSREmission:veto = 1. This will lead to a veto in ISR, FSR and -// MPI steps of pythia as expected using PowhegHooks in all the cases other than -// the case of FSR emission in resonance decay. Within resonance decays -// PowhegHooksBB4L takes over the control. -// -// Furthermore, this hook can also be used standalone without PowhegHooks, i.e. -// only the FSR emission from resonances will be vetoed (the default use in -// 1801.03944 and 1906.09166). In order to do that set -// POWHEG:bb4l:FSREmission:veto = 1 and POWHEG:veto = 0. Note that the this is -// not recommended for comparing against data but because it is easier to -// interpret it is often done in theoretical studies. -// -// Note that this version of the hook relies on the event "radtype" (1 for -// btilde, 2 for remnant) to be set by an external program, such as -// main-PYTHIA82-lhef in the radtype_ common block. -// There also exists a version of this hook in which the event "radtype" is -// read in from the .lhe file using pythia built in functionality. You need -// that version if you want to use this hook with main31.cc. -// -// -// # Expert use -// -// This hook also implements an alternative veto procedure which allows to -// assign a "SCALUP" type of scale to a resonance using the scaleResonance -// method. This is a much simpler veto but it is also clearly inferior as -// compared to the one implemented using the doVetoFSREmission method because -// the definition of the scale of the emission does not match the -// POWHEG-BOX-RES definition. Nevertheless, it can be activated using -// POWHEG:bb4l:ScaleResonance:veto = 1. Additionally one MUST switch off the -// other veto by calling on POWHEG:bb4l:FSREmission:veto = 0. -// -// The following settings are at the disposal of the user to control the -// behaviour of the hook -// - On/off switches for the veto: -// - POWHEG:bb4l:FSREmission:veto -// on/off switch for the default veto based on doFSREmission -// - POWHEG:bb4l:ScaleResonance:veto -// on/off switch for the alternative veto based on scaleResonance (only -// for expert users) -// - Important settings: -// - POWHEG:bb4l:ptMinVeto: MUST be set to the same value as the -// corresponding flag in POWHEG-BOX-RES -// - Alternatives for scale calculations -// - default: emission scale is calculated using POWHEG definitions and in -// the resonance rest frame -// - POWHEG:bb4l:FSREmission:vetoDipoleFrame: emission scale is calculated -// using POWHEG definitions in the dipole frame -// - POWHEG:bb4l:FSREmission:pTpythiaVeto: emission scale is calculated -// using Pythia definitions -// - Other flags: -// - POWHEG:bb4l:FSREmission:vetoQED: decides whether or not QED emission -// off quarks should also be vetoed (not implemented in the case of -// the ScaleResonance:veto) -// - POWHEG:bb4l:DEBUG: enables debug printouts on standard output - -#ifndef Pythia8_PowhegHooksBB4L_H -#define Pythia8_PowhegHooksBB4L_H - -#include "Pythia8/Pythia.h" -#include - -namespace Pythia8 { - -class PowhegHooksBB4L : public UserHooks { - - public: - - PowhegHooksBB4L() {} - ~PowhegHooksBB4L() { - std::cout << "Number of FSR vetoed in BB4l = " << nInResonanceFSRveto << std::endl; - } - - //--- Initialization ------------------------------------------------------- - bool initAfterBeams() { - // initialize settings of this class - vetoFSREmission = settingsPtr->flag("POWHEG:bb4l:FSREmission:veto"); - vetoDipoleFrame = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoDipoleFrame"); - pTpythiaVeto = settingsPtr->flag("POWHEG:bb4l:FSREmission:pTpythiaVeto"); - vetoQED = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoQED"); - scaleResonanceVeto = settingsPtr->flag("POWHEG:bb4l:ScaleResonance:veto"); - debug = settingsPtr->flag("POWHEG:bb4l:DEBUG"); - pTmin = settingsPtr->parm("POWHEG:bb4l:pTminVeto"); - vetoAllRadtypes = settingsPtr->flag("POWHEG:bb4l:vetoAllRadtypes"); - nInResonanceFSRveto = 0; - return true; - } - - //--- PROCESS LEVEL HOOK --------------------------------------------------- - // This hook gets triggered for each event before the shower starts, i.e. at - // the LHE level. We use it to calculate the scales of resonances. - inline bool canVetoProcessLevel() { return true; } - inline bool doVetoProcessLevel(Event &e) { - // extract the radtype from the event comment - stringstream ss; - // use eventattribute as comments not filled when using edm input - ss << infoPtr->getEventAttribute("#rwgt"); - string temp; - ss >> temp >> radtype; - assert(temp == "#rwgt"); - // we only calculate resonance scales for btilde events (radtype==1) - // remnant events are not vetoed - if (!vetoAllRadtypes && radtype==2) return false; - // find last top and the last anti-top in the record - int i_top = -1, i_atop = -1, i_wp = -1, i_wm = -1; - for (int i = 0; i < e.size(); i++) { - if (e[i].id() == 6) i_top = i; - if (e[i].id() == -6) i_atop = i; - if (e[i].id() == 24) i_wp = i; - if (e[i].id() == -24) i_wm = i; - } - // if found calculate the resonance scale - topresscale = findresscale(i_top, e); - // similarly for anti-top - atopresscale = findresscale(i_atop, e); - // and for W^+ and W^- - wpresscale = findresscale(i_wp, e); - wmresscale = findresscale(i_wm, e); - - // do not veto, ever - return false; - } - - //--- FSR EMISSION LEVEL HOOK ---------------------------------------------- - // This hook gets triggered everytime the parton shower attempts to attach - // a FSR emission. - inline bool canVetoFSREmission() { return vetoFSREmission; } - inline bool doVetoFSREmission(int sizeOld, const Event &e, int iSys, bool inResonance) { - - // FSR VETO INSIDE THE RESONANCE (if it is switched on) - if (inResonance && vetoFSREmission) { - - // get the participants of the splitting: the recoiler, the radiator and the emitted - int iRecAft = e.size() - 1; - int iEmt = e.size() - 2; - int iRadAft = e.size() - 3; - int iRadBef = e[iEmt].mother1(); - - // find the resonance the radiator originates from - int iRes = e[iRadBef].mother1(); - int distance = 1; - while ( iRes > 0 && (abs(e[iRes].id()) !=6 && abs(e[iRes].id()) != 24) ) { - iRes = e[iRes].mother1(); - distance ++; - } - if (iRes == 0) { - infoPtr->errorMsg("Warning in PowhegHooksBB4L::doVetoFSREmission: emission in resonance not from the top quark or from the W boson, not vetoing"); - return doVetoFSR(false,0); - } - int iResId = e[iRes].id(); - - // calculate the scale of the emission - double scale; - //using pythia pT definition ... - if(pTpythiaVeto) - scale = pTpythia(e, iRadAft, iEmt, iRecAft); - //.. or using POWHEG pT definition - else{ - Vec4 pr(e[iRadAft].p()), pe(e[iEmt].p()), pres(e[iRes].p()), prec(e[iRecAft].p()), psystem; - // The computation of the POWHEG pT can be done in the top rest frame or in the diple one. - // pdipole = pemt +prec +prad (after the emission) - // For the first emission off the top resonance pdipole = pw +pb (before the emission) = ptop - if(vetoDipoleFrame) - psystem = pr+pe+prec; - else - psystem = pres; - - // gluon splitting into two partons - if (e[iRadBef].id() == 21) - scale = gSplittingScale(psystem, pr, pe); - // quark emitting a gluon (or a photon) - else if (abs(e[iRadBef].id()) <= 5 && ((e[iEmt].id() == 21) && ! vetoQED) ) - scale = qSplittingScale(psystem, pr, pe); - // other stuff (which we should not veto) - else { - scale = 0; - } - } - - // compare the current splitting scale to the correct resonance scale - if (iResId == 6) { - if ( debug && scale > topresscale) - cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; - return doVetoFSR(scale > topresscale, scale); - } - else if (iResId == -6){ - if ( debug && scale > atopresscale) - cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; - return doVetoFSR(scale > atopresscale, scale); - } - else if (iResId == 24) { - if ( debug && scale > wpresscale) - cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; - return doVetoFSR(scale > wpresscale, scale); - } - else if (iResId == -24){ - if ( debug && scale > wmresscale) - cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; - return doVetoFSR(scale > wmresscale, scale); - } - else { - infoPtr->errorMsg("Error in PowhegHooksBB4L::doVetoFSREmission: unimplemented case"); - exit(-1); - } - } - // In CMSSW, the production process veto is done in EmissionVetoHook1.cc - // so for events outside resonance, nothing needs to be done here - else { - return false; - } - } - - inline bool doVetoFSR(bool condition, double scale) { - if (!vetoAllRadtypes && radtype==2) return false; - if (condition) { - nInResonanceFSRveto++; - return true; - } - return false; - } - - //--- SCALE RESONANCE HOOK ------------------------------------------------- - // called before each resonance decay shower - inline bool canSetResonanceScale() { return scaleResonanceVeto; } - // if the resonance is the (anti)top or W+/W- set the scale to: - // - if radtype=2 (remnant): resonance virtuality - // - if radtype=1 (btilde): - // - (a)topresscale/wp(m)resscale for tops and Ws - // - a large number otherwise - // if is not the top, set it to a big number - inline double scaleResonance(int iRes, const Event &e) { - if(!vetoAllRadtypes && radtype == 2) - return sqrt(e[iRes].m2Calc()); - else { - if (e[iRes].id() == 6) - return topresscale; - else if (e[iRes].id() == -6) - return atopresscale; - else if (e[iRes].id() == 24) - return wpresscale; - else if (e[iRes].id() == 24) - return wmresscale; - else - return 1e30; - } - } - - //--- Internal helper functions -------------------------------------------- - // Calculates the scale of the hardest emission from within the resonance system - // translated by Markus Seidel modified by Tomas Jezo - inline double findresscale( const int iRes, const Event& event) { - - // return large scale if the resonance position is ill defined - if (iRes < 0) return 1e30; - - // get number of resonance decay products - int nDau = event[iRes].daughterList().size(); - - // iRes is not decayed, return high scale equivalent to - // unrestricted shower - if (nDau == 0) { - return 1e30; - } - // iRes did not radiate, this means that POWHEG pt scale has - // evolved all the way down to pTmin - else if (nDau < 3) { - return pTmin; - } - // iRes is a (anti-)top quark - 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 and return it - return sqrt(2*pg*pb*pg.e()/pb.e()); - } - // iRes is a W+(-) boson - 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()); - // and return it - return sqrt(min(1-yq,1-ya)*pow2(csi)*q2/2); - } - // in any other case just return a high scale equivalent to - // unrestricted shower - return 1e30; - } - - // The following routine will match daughters of particle `e[iparticle]` to an expected pattern specified via the list of expected particle PDG ID's `ids`, - // id wildcard is specified as 0 if match is obtained, the positions and the momenta of these particles are returned in vectors `positions` and `momenta` - // respectively - // if exitOnExtraLegs==true, it will exit if the decay has more particles than expected, but not less - inline bool match_decay(int iparticle, const Event &e, const vector &ids, vector &positions, vector &momenta, bool exitOnExtraLegs = true){ - // compare sizes - if (e[iparticle].daughterList().size() != ids.size()) { - if (exitOnExtraLegs && e[iparticle].daughterList().size() > ids.size()) { - cout << "extra leg" << endl; - exit(-1); - } - return false; - } - // compare content - for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) { - int di = e[iparticle].daughterList()[i]; - if (ids[i] != 0 && e[di].id() != ids[i]) - return false; - } - // reset the positions and momenta vectors (because they may be reused) - positions.clear(); - momenta.clear(); - // construct the array of momenta - for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) { - int di = e[iparticle].daughterList()[i]; - positions.push_back(di); - momenta.push_back(e[di].p()); - } - return true; - } - - inline double qSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){ - p1.bstback(pt); - p2.bstback(pt); - return sqrt( 2*p1*p2*p2.e()/p1.e() ); - } - - inline double gSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){ - p1.bstback(pt); - p2.bstback(pt); - return sqrt( 2*p1*p2*p1.e()*p2.e()/(pow(p1.e()+p2.e(),2)) ); - } - - // Routines to calculate the pT (according to pTdefMode) in a FS splitting: - // i (radiator before) -> j (emitted after) k (radiator after) - // For the Pythia pT definition, a recoiler (after) must be specified. - // (INSPIRED BY pythia8F77_31.cc double pTpythia) - inline double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch, - int RecAfterBranch) - { - - // Convenient shorthands for later - Vec4 radVec = e[RadAfterBranch].p(); - Vec4 emtVec = e[EmtAfterBranch].p(); - Vec4 recVec = e[RecAfterBranch].p(); - int radID = e[RadAfterBranch].id(); - - // Calculate virtuality of splitting - Vec4 Q(radVec + emtVec); - double Qsq = Q.m2Calc(); - - // Mass term of radiator - double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ? - pow2(particleDataPtr->m0(radID)) : 0.; - - // z values for FSR - double z, pTnow; - // Construct 2 -> 3 variables - Vec4 sum = radVec + recVec + emtVec; - double m2Dip = sum.m2Calc(); - - double x1 = 2. * (sum * radVec) / m2Dip; - double x3 = 2. * (sum * emtVec) / m2Dip; - z = x1 / (x1 + x3); - pTnow = z * (1. - z); - - - // Virtuality - pTnow *= (Qsq - m2Rad); - - if (pTnow < 0.) { - cout << "Warning: pTpythia was negative" << endl; - return -1.; - } - else - return(sqrt(pTnow)); - } - - // Functions to return statistics about the veto - inline int getNInResonanceFSRVeto() { return nInResonanceFSRveto; } - - //-------------------------------------------------------------------------- - - private: - // FSR emission veto flags - bool vetoFSREmission, vetoQED; - // scale Resonance veto flags - double scaleResonanceVeto; - // other flags - bool debug; - bool vetoDipoleFrame; - bool pTpythiaVeto; - double pTmin; - bool vetoAllRadtypes; - // veto counter - int nInResonanceFSRveto; - // internal: resonance scales - double topresscale, atopresscale, wpresscale, wmresscale; - int radtype; -}; - -//========================================================================== - -} // end namespace Pythia8 - -#endif // end Pythia8_PowhegHooksBB4L_H +// PowhegHooksBB4L.h +// Rewritten by T. Jezo in 2021. With various contributions from S. Ferrario +// Ravasio, B. Nachman, P. Nason and M. Seidel. Inspired by +// ttb_NLO_dec/main-PYTHIA8.f by P. Nason and E. Re and by PowhegHooks.h by R. +// Corke. +// +// Adapted for CMSSW by Laurids Jeppe. + +// # Introduction +// +// This hook is intended for use together with POWHEG-BOX-RES/b_bbar_4l NLO LHE +// events. This also includes events in which one of the W bosons was +// re-decayed hadronically. (Note that LHE format version larger than 3 may not +// require this hook). +// +// The hook inherits from PowhegHooks and as such it indirectly implements the +// following: +// - doVetoMPIStep +// - doVetoISREmission +// - doVetoMPIEmission +// and it overloads: +// - doVetoFSREmission, which works as follows (if POWHEG:veto==1): +// - if inResonance==true it vetoes all the emission that is harder than +// the scale of its parent (anti-)top quark or W^+(-) +// - if inResonance==false, it calls PowhegHooks::doVetoISREmission +// and it also implements: +// - doVetoProcessLevel, which is never used for vetoing (i.e. it always +// returns false). Instead it is used for the calculation of reconance scales +// using LHE kinematics. +// +// This version of the hooks is only suitable for use with fully compatible +// POWHEG-BOX Les Houches readers (such the one in main-PYTHIA82-lhef but +// not the one in main31.cc.) +// +// +// # Basic use +// +// In order to use this hook you must replace all the declarations and +// constructor calls of PowhegHooks to PowhegHooksBB4L: +// +// PowhegHooks *powhegHooks; -> PowhegHooksBB4L *powhegHooks; +// *powhegHooks = new PowhegHooks(); -> *powhegHooks = new PowhegHooksBB4L(); +// +// In order to switch it on set POWHEG:veto = 1 and +// POWHEG:bb4l:FSREmission:veto = 1. This will lead to a veto in ISR, FSR and +// MPI steps of pythia as expected using PowhegHooks in all the cases other than +// the case of FSR emission in resonance decay. Within resonance decays +// PowhegHooksBB4L takes over the control. +// +// Furthermore, this hook can also be used standalone without PowhegHooks, i.e. +// only the FSR emission from resonances will be vetoed (the default use in +// 1801.03944 and 1906.09166). In order to do that set +// POWHEG:bb4l:FSREmission:veto = 1 and POWHEG:veto = 0. Note that the this is +// not recommended for comparing against data but because it is easier to +// interpret it is often done in theoretical studies. +// +// Note that this version of the hook relies on the event "radtype" (1 for +// btilde, 2 for remnant) to be set by an external program, such as +// main-PYTHIA82-lhef in the radtype_ common block. +// There also exists a version of this hook in which the event "radtype" is +// read in from the .lhe file using pythia built in functionality. You need +// that version if you want to use this hook with main31.cc. +// +// +// # Expert use +// +// This hook also implements an alternative veto procedure which allows to +// assign a "SCALUP" type of scale to a resonance using the scaleResonance +// method. This is a much simpler veto but it is also clearly inferior as +// compared to the one implemented using the doVetoFSREmission method because +// the definition of the scale of the emission does not match the +// POWHEG-BOX-RES definition. Nevertheless, it can be activated using +// POWHEG:bb4l:ScaleResonance:veto = 1. Additionally one MUST switch off the +// other veto by calling on POWHEG:bb4l:FSREmission:veto = 0. +// +// The following settings are at the disposal of the user to control the +// behaviour of the hook +// - On/off switches for the veto: +// - POWHEG:bb4l:FSREmission:veto +// on/off switch for the default veto based on doFSREmission +// - POWHEG:bb4l:ScaleResonance:veto +// on/off switch for the alternative veto based on scaleResonance (only +// for expert users) +// - Important settings: +// - POWHEG:bb4l:ptMinVeto: MUST be set to the same value as the +// corresponding flag in POWHEG-BOX-RES +// - Alternatives for scale calculations +// - default: emission scale is calculated using POWHEG definitions and in +// the resonance rest frame +// - POWHEG:bb4l:FSREmission:vetoDipoleFrame: emission scale is calculated +// using POWHEG definitions in the dipole frame +// - POWHEG:bb4l:FSREmission:pTpythiaVeto: emission scale is calculated +// using Pythia definitions +// - Other flags: +// - POWHEG:bb4l:FSREmission:vetoQED: decides whether or not QED emission +// off quarks should also be vetoed (not implemented in the case of +// the ScaleResonance:veto) +// - POWHEG:bb4l:DEBUG: enables debug printouts on standard output + +#ifndef Pythia8_PowhegHooksBB4L_H +#define Pythia8_PowhegHooksBB4L_H + +#include "Pythia8/Pythia.h" +#include + +namespace Pythia8 { + +class PowhegHooksBB4L : public UserHooks { + + public: + + PowhegHooksBB4L() {} + ~PowhegHooksBB4L() { + std::cout << "Number of FSR vetoed in BB4l = " << nInResonanceFSRveto << std::endl; + } + + //--- Initialization ------------------------------------------------------- + bool initAfterBeams() { + // initialize settings of this class + vetoFSREmission = settingsPtr->flag("POWHEG:bb4l:FSREmission:veto"); + vetoDipoleFrame = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoDipoleFrame"); + pTpythiaVeto = settingsPtr->flag("POWHEG:bb4l:FSREmission:pTpythiaVeto"); + vetoQED = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoQED"); + scaleResonanceVeto = settingsPtr->flag("POWHEG:bb4l:ScaleResonance:veto"); + debug = settingsPtr->flag("POWHEG:bb4l:DEBUG"); + pTmin = settingsPtr->parm("POWHEG:bb4l:pTminVeto"); + vetoAllRadtypes = settingsPtr->flag("POWHEG:bb4l:vetoAllRadtypes"); + nInResonanceFSRveto = 0; + return true; + } + + //--- PROCESS LEVEL HOOK --------------------------------------------------- + // This hook gets triggered for each event before the shower starts, i.e. at + // the LHE level. We use it to calculate the scales of resonances. + inline bool canVetoProcessLevel() { return true; } + inline bool doVetoProcessLevel(Event &e) { + // extract the radtype from the event comment + stringstream ss; + // use eventattribute as comments not filled when using edm input + ss << infoPtr->getEventAttribute("#rwgt"); + string temp; + ss >> temp >> radtype; + assert(temp == "#rwgt"); + // we only calculate resonance scales for btilde events (radtype==1) + // remnant events are not vetoed + if (!vetoAllRadtypes && radtype==2) return false; + // find last top and the last anti-top in the record + int i_top = -1, i_atop = -1, i_wp = -1, i_wm = -1; + for (int i = 0; i < e.size(); i++) { + if (e[i].id() == 6) i_top = i; + if (e[i].id() == -6) i_atop = i; + if (e[i].id() == 24) i_wp = i; + if (e[i].id() == -24) i_wm = i; + } + // if found calculate the resonance scale + topresscale = findresscale(i_top, e); + // similarly for anti-top + atopresscale = findresscale(i_atop, e); + // and for W^+ and W^- + wpresscale = findresscale(i_wp, e); + wmresscale = findresscale(i_wm, e); + + // do not veto, ever + return false; + } + + //--- FSR EMISSION LEVEL HOOK ---------------------------------------------- + // This hook gets triggered everytime the parton shower attempts to attach + // a FSR emission. + inline bool canVetoFSREmission() { return vetoFSREmission; } + inline bool doVetoFSREmission(int sizeOld, const Event &e, int iSys, bool inResonance) { + + // FSR VETO INSIDE THE RESONANCE (if it is switched on) + if (inResonance && vetoFSREmission) { + + // get the participants of the splitting: the recoiler, the radiator and the emitted + int iRecAft = e.size() - 1; + int iEmt = e.size() - 2; + int iRadAft = e.size() - 3; + int iRadBef = e[iEmt].mother1(); + + // find the resonance the radiator originates from + int iRes = e[iRadBef].mother1(); + int distance = 1; + while ( iRes > 0 && (abs(e[iRes].id()) !=6 && abs(e[iRes].id()) != 24) ) { + iRes = e[iRes].mother1(); + distance ++; + } + if (iRes == 0) { + infoPtr->errorMsg("Warning in PowhegHooksBB4L::doVetoFSREmission: emission in resonance not from the top quark or from the W boson, not vetoing"); + return doVetoFSR(false,0); + } + int iResId = e[iRes].id(); + + // calculate the scale of the emission + double scale; + //using pythia pT definition ... + if(pTpythiaVeto) + scale = pTpythia(e, iRadAft, iEmt, iRecAft); + //.. or using POWHEG pT definition + else{ + Vec4 pr(e[iRadAft].p()), pe(e[iEmt].p()), pres(e[iRes].p()), prec(e[iRecAft].p()), psystem; + // The computation of the POWHEG pT can be done in the top rest frame or in the diple one. + // pdipole = pemt +prec +prad (after the emission) + // For the first emission off the top resonance pdipole = pw +pb (before the emission) = ptop + if(vetoDipoleFrame) + psystem = pr+pe+prec; + else + psystem = pres; + + // gluon splitting into two partons + if (e[iRadBef].id() == 21) + scale = gSplittingScale(psystem, pr, pe); + // quark emitting a gluon (or a photon) + else if (abs(e[iRadBef].id()) <= 5 && ((e[iEmt].id() == 21) && ! vetoQED) ) + scale = qSplittingScale(psystem, pr, pe); + // other stuff (which we should not veto) + else { + scale = 0; + } + } + + // compare the current splitting scale to the correct resonance scale + if (iResId == 6) { + if ( debug && scale > topresscale) + cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; + return doVetoFSR(scale > topresscale, scale); + } + else if (iResId == -6){ + if ( debug && scale > atopresscale) + cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; + return doVetoFSR(scale > atopresscale, scale); + } + else if (iResId == 24) { + if ( debug && scale > wpresscale) + cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; + return doVetoFSR(scale > wpresscale, scale); + } + else if (iResId == -24){ + if ( debug && scale > wmresscale) + cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; + return doVetoFSR(scale > wmresscale, scale); + } + else { + infoPtr->errorMsg("Error in PowhegHooksBB4L::doVetoFSREmission: unimplemented case"); + exit(-1); + } + } + // In CMSSW, the production process veto is done in EmissionVetoHook1.cc + // so for events outside resonance, nothing needs to be done here + else { + return false; + } + } + + inline bool doVetoFSR(bool condition, double scale) { + if (!vetoAllRadtypes && radtype==2) return false; + if (condition) { + nInResonanceFSRveto++; + return true; + } + return false; + } + + //--- SCALE RESONANCE HOOK ------------------------------------------------- + // called before each resonance decay shower + inline bool canSetResonanceScale() { return scaleResonanceVeto; } + // if the resonance is the (anti)top or W+/W- set the scale to: + // - if radtype=2 (remnant): resonance virtuality + // - if radtype=1 (btilde): + // - (a)topresscale/wp(m)resscale for tops and Ws + // - a large number otherwise + // if is not the top, set it to a big number + inline double scaleResonance(int iRes, const Event &e) { + if(!vetoAllRadtypes && radtype == 2) + return sqrt(e[iRes].m2Calc()); + else { + if (e[iRes].id() == 6) + return topresscale; + else if (e[iRes].id() == -6) + return atopresscale; + else if (e[iRes].id() == 24) + return wpresscale; + else if (e[iRes].id() == 24) + return wmresscale; + else + return 1e30; + } + } + + //--- Internal helper functions -------------------------------------------- + // Calculates the scale of the hardest emission from within the resonance system + // translated by Markus Seidel modified by Tomas Jezo + inline double findresscale( const int iRes, const Event& event) { + + // return large scale if the resonance position is ill defined + if (iRes < 0) return 1e30; + + // get number of resonance decay products + int nDau = event[iRes].daughterList().size(); + + // iRes is not decayed, return high scale equivalent to + // unrestricted shower + if (nDau == 0) { + return 1e30; + } + // iRes did not radiate, this means that POWHEG pt scale has + // evolved all the way down to pTmin + else if (nDau < 3) { + return pTmin; + } + // iRes is a (anti-)top quark + 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 and return it + return sqrt(2*pg*pb*pg.e()/pb.e()); + } + // iRes is a W+(-) boson + 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()); + // and return it + return sqrt(min(1-yq,1-ya)*pow2(csi)*q2/2); + } + // in any other case just return a high scale equivalent to + // unrestricted shower + return 1e30; + } + + // The following routine will match daughters of particle `e[iparticle]` to an expected pattern specified via the list of expected particle PDG ID's `ids`, + // id wildcard is specified as 0 if match is obtained, the positions and the momenta of these particles are returned in vectors `positions` and `momenta` + // respectively + // if exitOnExtraLegs==true, it will exit if the decay has more particles than expected, but not less + inline bool match_decay(int iparticle, const Event &e, const vector &ids, vector &positions, vector &momenta, bool exitOnExtraLegs = true){ + // compare sizes + if (e[iparticle].daughterList().size() != ids.size()) { + if (exitOnExtraLegs && e[iparticle].daughterList().size() > ids.size()) { + cout << "extra leg" << endl; + exit(-1); + } + return false; + } + // compare content + for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) { + int di = e[iparticle].daughterList()[i]; + if (ids[i] != 0 && e[di].id() != ids[i]) + return false; + } + // reset the positions and momenta vectors (because they may be reused) + positions.clear(); + momenta.clear(); + // construct the array of momenta + for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) { + int di = e[iparticle].daughterList()[i]; + positions.push_back(di); + momenta.push_back(e[di].p()); + } + return true; + } + + inline double qSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){ + p1.bstback(pt); + p2.bstback(pt); + return sqrt( 2*p1*p2*p2.e()/p1.e() ); + } + + inline double gSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){ + p1.bstback(pt); + p2.bstback(pt); + return sqrt( 2*p1*p2*p1.e()*p2.e()/(pow(p1.e()+p2.e(),2)) ); + } + + // Routines to calculate the pT (according to pTdefMode) in a FS splitting: + // i (radiator before) -> j (emitted after) k (radiator after) + // For the Pythia pT definition, a recoiler (after) must be specified. + // (INSPIRED BY pythia8F77_31.cc double pTpythia) + inline double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch, + int RecAfterBranch) + { + + // Convenient shorthands for later + Vec4 radVec = e[RadAfterBranch].p(); + Vec4 emtVec = e[EmtAfterBranch].p(); + Vec4 recVec = e[RecAfterBranch].p(); + int radID = e[RadAfterBranch].id(); + + // Calculate virtuality of splitting + Vec4 Q(radVec + emtVec); + double Qsq = Q.m2Calc(); + + // Mass term of radiator + double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ? + pow2(particleDataPtr->m0(radID)) : 0.; + + // z values for FSR + double z, pTnow; + // Construct 2 -> 3 variables + Vec4 sum = radVec + recVec + emtVec; + double m2Dip = sum.m2Calc(); + + double x1 = 2. * (sum * radVec) / m2Dip; + double x3 = 2. * (sum * emtVec) / m2Dip; + z = x1 / (x1 + x3); + pTnow = z * (1. - z); + + + // Virtuality + pTnow *= (Qsq - m2Rad); + + if (pTnow < 0.) { + cout << "Warning: pTpythia was negative" << endl; + return -1.; + } + else + return(sqrt(pTnow)); + } + + // Functions to return statistics about the veto + inline int getNInResonanceFSRVeto() { return nInResonanceFSRveto; } + + //-------------------------------------------------------------------------- + + private: + // FSR emission veto flags + bool vetoFSREmission, vetoQED; + // scale Resonance veto flags + double scaleResonanceVeto; + // other flags + bool debug; + bool vetoDipoleFrame; + bool pTpythiaVeto; + double pTmin; + bool vetoAllRadtypes; + // veto counter + int nInResonanceFSRveto; + // internal: resonance scales + double topresscale, atopresscale, wpresscale, wmresscale; + int radtype; +}; + +//========================================================================== + +} // end namespace Pythia8 + +#endif // end Pythia8_PowhegHooksBB4L_H From 166e3b01bf2307e59b384b7df8db325d0250c689 Mon Sep 17 00:00:00 2001 From: Laurids Jeppe Date: Thu, 22 Feb 2024 18:19:16 +0100 Subject: [PATCH 5/6] Code format --- .../plugins/PowhegHooksBB4L.h | 203 +++++++++--------- 1 file changed, 105 insertions(+), 98 deletions(-) diff --git a/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h b/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h index 75e6fc7111dac..f0c5274f42778 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h +++ b/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h @@ -72,7 +72,7 @@ // POWHEG-BOX-RES definition. Nevertheless, it can be activated using // POWHEG:bb4l:ScaleResonance:veto = 1. Additionally one MUST switch off the // other veto by calling on POWHEG:bb4l:FSREmission:veto = 0. -// +// // The following settings are at the disposal of the user to control the // behaviour of the hook // - On/off switches for the veto: @@ -89,11 +89,11 @@ // the resonance rest frame // - POWHEG:bb4l:FSREmission:vetoDipoleFrame: emission scale is calculated // using POWHEG definitions in the dipole frame -// - POWHEG:bb4l:FSREmission:pTpythiaVeto: emission scale is calculated +// - POWHEG:bb4l:FSREmission:pTpythiaVeto: emission scale is calculated // using Pythia definitions // - Other flags: // - POWHEG:bb4l:FSREmission:vetoQED: decides whether or not QED emission -// off quarks should also be vetoed (not implemented in the case of +// off quarks should also be vetoed (not implemented in the case of // the ScaleResonance:veto) // - POWHEG:bb4l:DEBUG: enables debug printouts on standard output @@ -105,14 +105,10 @@ namespace Pythia8 { -class PowhegHooksBB4L : public UserHooks { - + class PowhegHooksBB4L : public UserHooks { public: - PowhegHooksBB4L() {} - ~PowhegHooksBB4L() { - std::cout << "Number of FSR vetoed in BB4l = " << nInResonanceFSRveto << std::endl; - } + ~PowhegHooksBB4L() { std::cout << "Number of FSR vetoed in BB4l = " << nInResonanceFSRveto << std::endl; } //--- Initialization ------------------------------------------------------- bool initAfterBeams() { @@ -130,7 +126,7 @@ class PowhegHooksBB4L : public UserHooks { } //--- PROCESS LEVEL HOOK --------------------------------------------------- - // This hook gets triggered for each event before the shower starts, i.e. at + // This hook gets triggered for each event before the shower starts, i.e. at // the LHE level. We use it to calculate the scales of resonances. inline bool canVetoProcessLevel() { return true; } inline bool doVetoProcessLevel(Event &e) { @@ -143,14 +139,19 @@ class PowhegHooksBB4L : public UserHooks { assert(temp == "#rwgt"); // we only calculate resonance scales for btilde events (radtype==1) // remnant events are not vetoed - if (!vetoAllRadtypes && radtype==2) return false; + if (!vetoAllRadtypes && radtype == 2) + return false; // find last top and the last anti-top in the record int i_top = -1, i_atop = -1, i_wp = -1, i_wm = -1; for (int i = 0; i < e.size(); i++) { - if (e[i].id() == 6) i_top = i; - if (e[i].id() == -6) i_atop = i; - if (e[i].id() == 24) i_wp = i; - if (e[i].id() == -24) i_wm = i; + if (e[i].id() == 6) + i_top = i; + if (e[i].id() == -6) + i_atop = i; + if (e[i].id() == 24) + i_wp = i; + if (e[i].id() == -24) + i_wm = i; } // if found calculate the resonance scale topresscale = findresscale(i_top, e); @@ -169,11 +170,9 @@ class PowhegHooksBB4L : public UserHooks { // a FSR emission. inline bool canVetoFSREmission() { return vetoFSREmission; } inline bool doVetoFSREmission(int sizeOld, const Event &e, int iSys, bool inResonance) { - // FSR VETO INSIDE THE RESONANCE (if it is switched on) if (inResonance && vetoFSREmission) { - - // get the participants of the splitting: the recoiler, the radiator and the emitted + // get the participants of the splitting: the recoiler, the radiator and the emitted int iRecAft = e.size() - 1; int iEmt = e.size() - 2; int iRadAft = e.size() - 3; @@ -182,66 +181,68 @@ class PowhegHooksBB4L : public UserHooks { // find the resonance the radiator originates from int iRes = e[iRadBef].mother1(); int distance = 1; - while ( iRes > 0 && (abs(e[iRes].id()) !=6 && abs(e[iRes].id()) != 24) ) { + while (iRes > 0 && (abs(e[iRes].id()) != 6 && abs(e[iRes].id()) != 24)) { iRes = e[iRes].mother1(); - distance ++; + distance++; } if (iRes == 0) { - infoPtr->errorMsg("Warning in PowhegHooksBB4L::doVetoFSREmission: emission in resonance not from the top quark or from the W boson, not vetoing"); - return doVetoFSR(false,0); + infoPtr->errorMsg( + "Warning in PowhegHooksBB4L::doVetoFSREmission: emission in resonance not from the top quark or from the " + "W boson, not vetoing"); + return doVetoFSR(false, 0); } int iResId = e[iRes].id(); // calculate the scale of the emission double scale; //using pythia pT definition ... - if(pTpythiaVeto) + if (pTpythiaVeto) scale = pTpythia(e, iRadAft, iEmt, iRecAft); //.. or using POWHEG pT definition - else{ + else { Vec4 pr(e[iRadAft].p()), pe(e[iEmt].p()), pres(e[iRes].p()), prec(e[iRecAft].p()), psystem; // The computation of the POWHEG pT can be done in the top rest frame or in the diple one. // pdipole = pemt +prec +prad (after the emission) - // For the first emission off the top resonance pdipole = pw +pb (before the emission) = ptop - if(vetoDipoleFrame) - psystem = pr+pe+prec; + // For the first emission off the top resonance pdipole = pw +pb (before the emission) = ptop + if (vetoDipoleFrame) + psystem = pr + pe + prec; else psystem = pres; - + // gluon splitting into two partons if (e[iRadBef].id() == 21) - scale = gSplittingScale(psystem, pr, pe); + scale = gSplittingScale(psystem, pr, pe); // quark emitting a gluon (or a photon) - else if (abs(e[iRadBef].id()) <= 5 && ((e[iEmt].id() == 21) && ! vetoQED) ) - scale = qSplittingScale(psystem, pr, pe); + else if (abs(e[iRadBef].id()) <= 5 && ((e[iEmt].id() == 21) && !vetoQED)) + scale = qSplittingScale(psystem, pr, pe); // other stuff (which we should not veto) else { - scale = 0; + scale = 0; } } - + // compare the current splitting scale to the correct resonance scale if (iResId == 6) { - if ( debug && scale > topresscale) - cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; + if (debug && scale > topresscale) + cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " + << scale << endl; return doVetoFSR(scale > topresscale, scale); - } - else if (iResId == -6){ - if ( debug && scale > atopresscale) - cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; + } else if (iResId == -6) { + if (debug && scale > atopresscale) + cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " + << scale << endl; return doVetoFSR(scale > atopresscale, scale); - } - else if (iResId == 24) { - if ( debug && scale > wpresscale) - cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; + } else if (iResId == 24) { + if (debug && scale > wpresscale) + cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " + << scale << endl; return doVetoFSR(scale > wpresscale, scale); - } - else if (iResId == -24){ - if ( debug && scale > wmresscale) - cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl; + } else if (iResId == -24) { + if (debug && scale > wmresscale) + cout << iResId << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " + << scale << endl; return doVetoFSR(scale > wmresscale, scale); - } - else { + } else { infoPtr->errorMsg("Error in PowhegHooksBB4L::doVetoFSREmission: unimplemented case"); exit(-1); } @@ -253,8 +254,9 @@ class PowhegHooksBB4L : public UserHooks { } } - inline bool doVetoFSR(bool condition, double scale) { - if (!vetoAllRadtypes && radtype==2) return false; + inline bool doVetoFSR(bool condition, double scale) { + if (!vetoAllRadtypes && radtype == 2) + return false; if (condition) { nInResonanceFSRveto++; return true; @@ -266,13 +268,13 @@ class PowhegHooksBB4L : public UserHooks { // called before each resonance decay shower inline bool canSetResonanceScale() { return scaleResonanceVeto; } // if the resonance is the (anti)top or W+/W- set the scale to: - // - if radtype=2 (remnant): resonance virtuality - // - if radtype=1 (btilde): + // - if radtype=2 (remnant): resonance virtuality + // - if radtype=1 (btilde): // - (a)topresscale/wp(m)resscale for tops and Ws // - a large number otherwise // if is not the top, set it to a big number inline double scaleResonance(int iRes, const Event &e) { - if(!vetoAllRadtypes && radtype == 2) + if (!vetoAllRadtypes && radtype == 2) return sqrt(e[iRes].m2Calc()); else { if (e[iRes].id() == 6) @@ -291,10 +293,10 @@ class PowhegHooksBB4L : public UserHooks { //--- Internal helper functions -------------------------------------------- // Calculates the scale of the hardest emission from within the resonance system // translated by Markus Seidel modified by Tomas Jezo - inline double findresscale( const int iRes, const Event& event) { - + inline double findresscale(const int iRes, const Event &event) { // return large scale if the resonance position is ill defined - if (iRes < 0) return 1e30; + if (iRes < 0) + return 1e30; // get number of resonance decay products int nDau = event[iRes].daughterList().size(); @@ -308,16 +310,19 @@ class PowhegHooksBB4L : public UserHooks { // evolved all the way down to pTmin else if (nDau < 3) { return pTmin; - } + } // iRes is a (anti-)top quark 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; + 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 @@ -329,7 +334,7 @@ class PowhegHooksBB4L : public UserHooks { pg.bstback(event[iRes].p()); // Calculate scale and return it - return sqrt(2*pg*pb*pg.e()/pb.e()); + return sqrt(2 * pg * pb * pg.e() / pb.e()); } // iRes is a W+(-) boson else if (abs(event[iRes].id()) == 24) { @@ -337,11 +342,14 @@ class PowhegHooksBB4L : public UserHooks { 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; + 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()); @@ -349,15 +357,15 @@ class PowhegHooksBB4L : public UserHooks { 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()); - // and return it - return sqrt(min(1-yq,1-ya)*pow2(csi)*q2/2); + 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()); + // and return it + return sqrt(min(1 - yq, 1 - ya) * pow2(csi) * q2 / 2); } // in any other case just return a high scale equivalent to // unrestricted shower @@ -368,19 +376,24 @@ class PowhegHooksBB4L : public UserHooks { // id wildcard is specified as 0 if match is obtained, the positions and the momenta of these particles are returned in vectors `positions` and `momenta` // respectively // if exitOnExtraLegs==true, it will exit if the decay has more particles than expected, but not less - inline bool match_decay(int iparticle, const Event &e, const vector &ids, vector &positions, vector &momenta, bool exitOnExtraLegs = true){ + inline bool match_decay(int iparticle, + const Event &e, + const vector &ids, + vector &positions, + vector &momenta, + bool exitOnExtraLegs = true) { // compare sizes if (e[iparticle].daughterList().size() != ids.size()) { if (exitOnExtraLegs && e[iparticle].daughterList().size() > ids.size()) { cout << "extra leg" << endl; exit(-1); } - return false; + return false; } // compare content for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) { int di = e[iparticle].daughterList()[i]; - if (ids[i] != 0 && e[di].id() != ids[i]) + if (ids[i] != 0 && e[di].id() != ids[i]) return false; } // reset the positions and momenta vectors (because they may be reused) @@ -395,41 +408,37 @@ class PowhegHooksBB4L : public UserHooks { return true; } - inline double qSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){ + inline double qSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2) { p1.bstback(pt); p2.bstback(pt); - return sqrt( 2*p1*p2*p2.e()/p1.e() ); + return sqrt(2 * p1 * p2 * p2.e() / p1.e()); } - inline double gSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){ + inline double gSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2) { p1.bstback(pt); - p2.bstback(pt); - return sqrt( 2*p1*p2*p1.e()*p2.e()/(pow(p1.e()+p2.e(),2)) ); + p2.bstback(pt); + return sqrt(2 * p1 * p2 * p1.e() * p2.e() / (pow(p1.e() + p2.e(), 2))); } // Routines to calculate the pT (according to pTdefMode) in a FS splitting: // i (radiator before) -> j (emitted after) k (radiator after) // For the Pythia pT definition, a recoiler (after) must be specified. // (INSPIRED BY pythia8F77_31.cc double pTpythia) - inline double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch, - int RecAfterBranch) - { - + inline double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch) { // Convenient shorthands for later Vec4 radVec = e[RadAfterBranch].p(); Vec4 emtVec = e[EmtAfterBranch].p(); Vec4 recVec = e[RecAfterBranch].p(); - int radID = e[RadAfterBranch].id(); + int radID = e[RadAfterBranch].id(); // Calculate virtuality of splitting Vec4 Q(radVec + emtVec); double Qsq = Q.m2Calc(); // Mass term of radiator - double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ? - pow2(particleDataPtr->m0(radID)) : 0.; + double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ? pow2(particleDataPtr->m0(radID)) : 0.; - // z values for FSR + // z values for FSR double z, pTnow; // Construct 2 -> 3 variables Vec4 sum = radVec + recVec + emtVec; @@ -437,19 +446,17 @@ class PowhegHooksBB4L : public UserHooks { double x1 = 2. * (sum * radVec) / m2Dip; double x3 = 2. * (sum * emtVec) / m2Dip; - z = x1 / (x1 + x3); + z = x1 / (x1 + x3); pTnow = z * (1. - z); - // Virtuality pTnow *= (Qsq - m2Rad); if (pTnow < 0.) { cout << "Warning: pTpythia was negative" << endl; return -1.; - } - else - return(sqrt(pTnow)); + } else + return (sqrt(pTnow)); } // Functions to return statistics about the veto @@ -473,10 +480,10 @@ class PowhegHooksBB4L : public UserHooks { // internal: resonance scales double topresscale, atopresscale, wpresscale, wmresscale; int radtype; -}; + }; -//========================================================================== + //========================================================================== -} // end namespace Pythia8 +} // end namespace Pythia8 -#endif // end Pythia8_PowhegHooksBB4L_H +#endif // end Pythia8_PowhegHooksBB4L_H From f39d218f0f59f642b3fec14d586abc9b0a2bbb0f Mon Sep 17 00:00:00 2001 From: Laurids Jeppe Date: Tue, 12 Mar 2024 18:16:06 +0100 Subject: [PATCH 6/6] Remove unused variable --- GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h b/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h index f0c5274f42778..781d36c87209a 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h +++ b/GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h @@ -180,10 +180,8 @@ namespace Pythia8 { // find the resonance the radiator originates from int iRes = e[iRadBef].mother1(); - int distance = 1; while (iRes > 0 && (abs(e[iRes].id()) != 6 && abs(e[iRes].id()) != 24)) { iRes = e[iRes].mother1(); - distance++; } if (iRes == 0) { infoPtr->errorMsg(