diff --git a/MatrixElement/Matchbox/Matching/QTildeMatching.cc b/MatrixElement/Matchbox/Matching/QTildeMatching.cc --- a/MatrixElement/Matchbox/Matching/QTildeMatching.cc +++ b/MatrixElement/Matchbox/Matching/QTildeMatching.cc @@ -1,524 +1,538 @@ // -*- C++ -*- // // QTildeMatching.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the QTildeMatching class. // #include "QTildeMatching.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h" #include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h" +#include "Herwig/Shower/QTilde/Kinematics/KinematicHelpers.h" using namespace Herwig; QTildeMatching::QTildeMatching() : theCorrectForXZMismatch(true) {} QTildeMatching::~QTildeMatching() {} IBPtr QTildeMatching::clone() const { return new_ptr(*this); } IBPtr QTildeMatching::fullclone() const { return new_ptr(*this); } void QTildeMatching::checkCutoff() { if ( showerTildeKinematics() ) { showerTildeKinematics()-> prepare(realCXComb(),bornCXComb()); showerTildeKinematics()->dipole(dipole()); showerTildeKinematics()->getShowerVariables(); } } void QTildeMatching::getShowerVariables() { // already filled from checkCutoff in this case if ( showerTildeKinematics() ) return; // get the shower variables calculateShowerVariables(); // check for the cutoff dipole()->isAboveCutoff(isAboveCutoff()); // get the hard scale dipole()->showerHardScale(hardScale()); // check for phase space dipole()->isInShowerPhasespace(isInShowerPhasespace()); } bool QTildeMatching::isInShowerPhasespace() const { - assert((theQTildeSudakov->cutOffOption() == 0 || theQTildeSudakov->cutOffOption() == 2) && - "implementation only provided for default and pt cutoff"); - Energy qtildeHard = ZERO; Energy qtilde = dipole()->showerScale(); assert(!dipole()->showerParameters().empty()); double z = dipole()->showerParameters()[0]; // FF if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() > 1 ) { qtildeHard = theQTildeFinder-> calculateFinalFinalScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()], bornCXComb()->meMomenta()[dipole()->bornSpectator()]).first; } // FI if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() < 2 ) { qtildeHard = theQTildeFinder-> calculateInitialFinalScales(bornCXComb()->meMomenta()[dipole()->bornSpectator()], bornCXComb()->meMomenta()[dipole()->bornEmitter()],false).second; } // IF if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() > 1 ) { qtildeHard = theQTildeFinder-> calculateInitialFinalScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()], bornCXComb()->meMomenta()[dipole()->bornSpectator()],false).first; if ( z < (dipole()->bornEmitter() == 0 ? bornCXComb()->lastX1() : bornCXComb()->lastX2()) ) return false; } // II if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() < 2 ) { qtildeHard = theQTildeFinder-> calculateInitialInitialScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()], bornCXComb()->meMomenta()[dipole()->bornSpectator()]).first; if ( z < (dipole()->bornEmitter() == 0 ? bornCXComb()->lastX1() : bornCXComb()->lastX2()) ) return false; } + Energy2 pt2 = ZERO; - Energy Qg = theQTildeSudakov->kinScale(); - Energy2 pt2 = ZERO; + const vector & masses = theQTildeSudakov->virtualMasses({{ + bornCXComb()->mePartonData()[dipole()->bornEmitter() ], + realCXComb()->mePartonData()[dipole()->realEmitter() ], + realCXComb()->mePartonData()[dipole()->realEmission()] + }}); + + const Energy2 m22 = sqr(masses[2]); + if ( dipole()->bornEmitter() > 1 ) { - Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass()); - if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g ) - pt2 = sqr(z*(1.-z)*qtilde) - sqr(mu); - else - pt2 = sqr(z*(1.-z)*qtilde) - sqr((1.-z)*mu) - z*sqr(Qg); + const Energy2 m02 = sqr(masses[0]); + const Energy2 m12 = sqr(masses[1]); + + pt2 = QTildeKinematics::pT2_FSR(sqr(qtilde),z,m02,m12,m22); } - if ( dipole()->bornEmitter() < 2 ) { - pt2 = sqr((1.-z)*qtilde) - z*sqr(Qg); + else { + pt2 = QTildeKinematics::pT2_ISR(sqr(qtilde),z,m22); } if ( pt2 < max(theQTildeSudakov->pT2min(),sqr(safeCut()) )) return false; bool hardVeto = restrictPhasespace() && sqrt(pt2) >= dipole()->showerHardScale(); return qtilde <= qtildeHard && !hardVeto; } bool QTildeMatching::isAboveCutoff() const { - - assert((theQTildeSudakov->cutOffOption() == 0 || theQTildeSudakov->cutOffOption() == 2) && - "implementation only provided for default and pt cutoff"); Energy qtilde = dipole()->showerScale(); assert(!dipole()->showerParameters().empty()); double z = dipole()->showerParameters()[0]; - Energy Qg = theQTildeSudakov->kinScale(); + + + const vector & masses = theQTildeSudakov->virtualMasses({{ + bornCXComb()->mePartonData()[dipole()->bornEmitter() ], + realCXComb()->mePartonData()[dipole()->realEmitter() ], + realCXComb()->mePartonData()[dipole()->realEmission()] + }}); + + const Energy2 m22 = sqr(masses[2]); + if ( dipole()->bornEmitter() > 1 ) { - Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass()); - if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g ) - return sqr(z*(1.-z)*qtilde) - sqr(mu) >= - max(theQTildeSudakov->pT2min(),sqr(safeCut())); - else - return sqr(z*(1.-z)*qtilde) - sqr((1.-z)*mu) - z*sqr(Qg) >= - max(theQTildeSudakov->pT2min(),sqr(safeCut())); + const Energy2 m02 = sqr(masses[0]); + const Energy2 m12 = sqr(masses[1]); + + const Energy2 pt2 = QTildeKinematics::pT2_FSR(sqr(qtilde),z,m02,m12,m22); + + return pt2 >= max(theQTildeSudakov->pT2min(),sqr(safeCut())); } - if ( dipole()->bornEmitter() < 2 ) { - return - sqr((1.-z)*qtilde) - z*sqr(Qg) >= - max(theQTildeSudakov->pT2min(),sqr(safeCut())); + else { + const Energy2 pt2 = QTildeKinematics::pT2_ISR(sqr(qtilde),z,m22); + return pt2 >= max(theQTildeSudakov->pT2min(),sqr(safeCut())); } return false; } CrossSection QTildeMatching::dSigHatDR() const { assert(!dipole()->showerParameters().empty()); pair vars = make_pair(sqr(dipole()->showerScale()), dipole()->showerParameters()[0]); pair ij(dipole()->bornEmitter(), dipole()->bornSpectator()); double ccme2 = dipole()->underlyingBornME()->largeNColourCorrelatedME2(ij,theLargeNBasis); if(ccme2==0.)return 0.*nanobarn; double lnme2=dipole()->underlyingBornME()->largeNME2(theLargeNBasis); if(lnme2==0){ generator()->log() <<"\nQTildeMatching: "; generator()->log() <<"\n largeNME2 is ZERO, while largeNColourCorrelatedME2 is not ZERO." ; generator()->log() <<"\n This is too seriuos.\n" ; generator()->log() << Exception::runerror; } ccme2 *= dipole()->underlyingBornME()->me2() /lnme2; Energy2 prop = ZERO; if ( dipole()->bornEmitter() > 1 ) { prop = (realCXComb()->meMomenta()[dipole()->realEmitter()] + realCXComb()->meMomenta()[dipole()->realEmission()]).m2() - bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2(); } else { prop = 2.*vars.second*(realCXComb()->meMomenta()[dipole()->realEmitter()]* realCXComb()->meMomenta()[dipole()->realEmission()]); } // note alphas included downstream from subtractionScaleWeight() double xme2 = -8.*Constants::pi*ccme2*splitFn(vars)*realXComb()->lastSHat()/prop; xme2 *= pow(realCXComb()->lastSHat() / bornCXComb()->lastSHat(), bornCXComb()->mePartonData().size()-4.); double bornPDF = bornPDFWeight(dipole()->underlyingBornME()->lastScale()); if ( bornPDF == 0.0 ) return ZERO; xme2 *= bornPDF; xme2 *= dipole()->realEmissionME()->finalStateSymmetry() / dipole()->underlyingBornME()->finalStateSymmetry(); // take care of mismatch between z and x as we are approaching the // hard phase space boundary // TODO get rid of this useless scale option business and simplify PDF handling in here if ( dipole()->bornEmitter() < 2 && theCorrectForXZMismatch ) { Energy2 emissionScale = ZERO; if ( emissionScaleInSubtraction() == showerScale ) { emissionScale = showerFactorizationScale(); } else if ( emissionScaleInSubtraction() == realScale ) { emissionScale = dipole()->realEmissionME()->lastScale(); } else if ( emissionScaleInSubtraction() == bornScale ) { emissionScale = dipole()->underlyingBornME()->lastScale(); } double xzMismatch = dipole()->subtractionParameters()[0] / dipole()->showerParameters()[0]; double realCorrectedPDF = dipole()->bornEmitter() == 0 ? dipole()->realEmissionME()->pdf1(emissionScale,theExtrapolationX, xzMismatch) : dipole()->realEmissionME()->pdf2(emissionScale,theExtrapolationX, xzMismatch); double realPDF = dipole()->bornEmitter() == 0 ? dipole()->realEmissionME()->pdf1(emissionScale,theExtrapolationX,1.0) : dipole()->realEmissionME()->pdf2(emissionScale,theExtrapolationX,1.0); if ( realPDF == 0.0 || realCorrectedPDF == 0.0 ) return ZERO; xme2 *= realCorrectedPDF / realPDF; } Energy qtilde = sqrt(vars.first); double z = vars.second; Energy2 pt2 = ZERO; - Energy Qg = theQTildeSudakov->kinScale(); + + const vector & masses = theQTildeSudakov->virtualMasses({{ + bornCXComb()->mePartonData()[dipole()->bornEmitter() ], + realCXComb()->mePartonData()[dipole()->realEmitter() ], + realCXComb()->mePartonData()[dipole()->realEmission()] + }}); + + const Energy2 m22 = sqr(masses[2]); if ( dipole()->bornEmitter() > 1 ) { - Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass()); - if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g ) - pt2 = sqr(z*(1.-z)*qtilde) - sqr(mu); - else - pt2 = sqr(z*(1.-z)*qtilde) - sqr((1.-z)*mu) - z*sqr(Qg); + const Energy2 m02 = sqr(masses[0]); + const Energy2 m12 = sqr(masses[1]); + + pt2 = QTildeKinematics::pT2_FSR(sqr(qtilde),z,m02,m12,m22); } - if ( dipole()->bornEmitter() < 2 ) { - pt2 = sqr((1.-z)*qtilde) - z*sqr(Qg); + else { + pt2 = QTildeKinematics::pT2_ISR(sqr(qtilde),z,m22); } assert(pt2 >= ZERO); if ( profileScales() ) xme2 *= profileScales()->hardScaleProfile(dipole()->showerHardScale(),sqrt(pt2)); CrossSection res = sqr(hbarc) * realXComb()->jacobian() * subtractionScaleWeight() * xme2 / (2. * realXComb()->lastSHat()); return res; } double QTildeMatching::me2() const { throw Exception() << "QTildeMatching::me2(): Not intented to use. Disable the ShowerApproximationGenerator." << Exception::runerror; return 0.; } void QTildeMatching::calculateShowerVariables() const { Lorentz5Momentum n; Energy2 Q2 = ZERO; const Lorentz5Momentum& pb = bornCXComb()->meMomenta()[dipole()->bornEmitter()]; const Lorentz5Momentum& pc = bornCXComb()->meMomenta()[dipole()->bornSpectator()]; if ( dipole()->bornEmitter() > 1 ) { Q2 = (pb+pc).m2(); } else { Q2 = -(pb-pc).m2(); } if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() > 1 ) { double b = sqr(bornCXComb()->meMomenta()[dipole()->bornEmitter()].m())/Q2; double c = sqr(bornCXComb()->meMomenta()[dipole()->bornSpectator()].m())/Q2; double lambda = sqrt(1.+sqr(b)+sqr(c)-2.*b-2.*c-2.*b*c); n = (1.-0.5*(1.-b+c-lambda))*pc - 0.5*(1.-b+c-lambda)*pb; } if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() < 2 ) { n = bornCXComb()->meMomenta()[dipole()->bornSpectator()]; } if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() > 1 ) { double c = sqr(bornCXComb()->meMomenta()[dipole()->bornSpectator()].m())/Q2; n = (1.+c)*pc - c*pb; } if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() < 2 ) { n = bornCXComb()->meMomenta()[dipole()->bornSpectator()]; } // the light-cone condition is numerically not very stable, so we // explicitly push it on the light-cone here n.setMass(ZERO); n.rescaleEnergy(); double z = 0.0; if ( dipole()->bornEmitter() > 1 ) { z = 1. - (n*realCXComb()->meMomenta()[dipole()->realEmission()])/ (n*bornCXComb()->meMomenta()[dipole()->bornEmitter()]); } else { z = 1. - (n*realCXComb()->meMomenta()[dipole()->realEmission()])/ (n*realCXComb()->meMomenta()[dipole()->realEmitter()]); } // allow small violations (numerical inaccuracies) if ( z <= 0 && z >= -1e-6 ) { z = std::numeric_limits::epsilon(); } else if ( z >= 1 && z <= 1+1e-6 ) { z = 1-std::numeric_limits::epsilon(); } Energy2 qtilde2 = ZERO; Energy2 q2 = ZERO; if ( dipole()->bornEmitter() > 1 ) { q2 = (realCXComb()->meMomenta()[dipole()->realEmitter()] + realCXComb()->meMomenta()[dipole()->realEmission()]).m2(); qtilde2 = (q2 - bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2())/(z*(1.-z)); } else { q2 = -(realCXComb()->meMomenta()[dipole()->realEmitter()] - realCXComb()->meMomenta()[dipole()->realEmission()]).m2(); qtilde2 = (q2 + bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2())/(1.-z); } if ( qtilde2 < ZERO ) { qtilde2 = ZERO; } assert(qtilde2 >= ZERO && z > 0.0 && z < 1.0); dipole()->showerScale(sqrt(qtilde2)); dipole()->showerParameters().resize(1); dipole()->showerParameters()[0] = z; } double QTildeMatching::splitFn(const pair& vars) const { const Energy2& qtilde2 = vars.first; const double z = vars.second; double Nc = SM().Nc(); // final state branching if ( dipole()->bornEmitter() > 1 ) { // final state quark quark branching if ( abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 7 ) { Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass(); return ((sqr(Nc)-1.)/(2.*Nc))*(1+sqr(z)-2.*sqr(m)/(z*qtilde2))/(1.-z); } // final state gluon branching if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g ) { if ( realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) { // ATTENTION the factor 2 here is intentional as it cancels to the 1/2 // stemming from the large-N colour correlator return 2.*Nc*(z/(1.-z)+(1.-z)/z+z*(1.-z)); } if ( abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) { Energy m = realCXComb()->mePartonData()[dipole()->realEmission()]->hardProcessMass(); return (1./2.)*(1.-2.*z*(1.-z)+2.*sqr(m)/(z*(1.-z)*qtilde2)); } } // final state squark branching if ((abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) > 1000000 && abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 1000007) || (abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) > 2000000 && abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 2000007)){ Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass(); return ((sqr(Nc)-1.)/Nc)*(z-sqr(m)/(z*qtilde2))/(1.-z); } // final state gluino branching if (bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == 1000021){ Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass(); return Nc*(1.+sqr(z)-2.*sqr(m)/(z*qtilde2))/(1.-z); } } // initial state branching if ( dipole()->bornEmitter() < 2 ) { // g/g if ( realCXComb()->mePartonData()[dipole()->realEmitter()]->id() == ParticleID::g && realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) { // see above for factor of 2 return 2.*Nc*(z/(1.-z)+(1.-z)/z+z*(1.-z)); } // q/q if ( abs(realCXComb()->mePartonData()[dipole()->realEmitter()]->id()) < 7 && realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) { return ((sqr(Nc)-1.)/(2.*Nc))*(1+sqr(z))/(1.-z); } // g/q if ( realCXComb()->mePartonData()[dipole()->realEmitter()]->id() == ParticleID::g && abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) { return (1./2.)*(1.-2.*z*(1.-z)); } // q/g if ( abs(realCXComb()->mePartonData()[dipole()->realEmitter()]->id()) < 7 && abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) { return ((sqr(Nc)-1.)/(2.*Nc))*(1+sqr(1.-z))/z; } } return 0.0; } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void QTildeMatching::doinit() { assert(theShowerHandler && theQTildeFinder && theQTildeSudakov); theShowerHandler->init(); theQTildeFinder->init(); theQTildeSudakov->init(); hardScaleFactor(theShowerHandler->hardScaleFactor()); factorizationScaleFactor(theShowerHandler->factorizationScaleFactor()); renormalizationScaleFactor(theShowerHandler->renormalizationScaleFactor()); profileScales(theShowerHandler->profileScales()); restrictPhasespace(theShowerHandler->restrictPhasespace()); hardScaleIsMuF(theShowerHandler->hardScaleIsMuF()); ShowerApproximation::doinit(); } void QTildeMatching::doinitrun() { assert(theShowerHandler && theQTildeFinder && theQTildeSudakov); theShowerHandler->initrun(); theQTildeFinder->initrun(); theQTildeSudakov->initrun(); ShowerApproximation::doinitrun(); } void QTildeMatching::persistentOutput(PersistentOStream & os) const { os << theQTildeFinder << theQTildeSudakov << theShowerHandler << theCorrectForXZMismatch; } void QTildeMatching::persistentInput(PersistentIStream & is, int) { is >> theQTildeFinder >> theQTildeSudakov >> theShowerHandler >> theCorrectForXZMismatch; } // *** Attention *** The following static variable is needed for the type // description system in ThePEG. Please check that the template arguments // are correct (the class and its base class), and that the constructor // arguments are correct (the class name and the name of the dynamically // loadable library where the class implementation can be found). DescribeClass describeHerwigQTildeMatching("Herwig::QTildeMatching", "HwShower.so HwQTildeMatching.so"); void QTildeMatching::Init() { static ClassDocumentation documentation ("QTildeMatching implements NLO matching with the default shower."); static Reference interfaceQTildeFinder ("QTildeFinder", "Set the partner finder to calculate hard scales.", &QTildeMatching::theQTildeFinder, false, false, true, false, false); interfaceQTildeFinder.rank(-1); static Reference interfaceQTildeSudakov ("QTildeSudakov", "Set the partner finder to calculate hard scales.", &QTildeMatching::theQTildeSudakov, false, false, true, false, false); interfaceQTildeSudakov.rank(-1); static Reference interfaceShowerHandler ("ShowerHandler", "The QTilde shower handler to use.", &QTildeMatching::theShowerHandler, false, false, true, true, false); interfaceShowerHandler.rank(-1); static Switch interfaceCorrectForXZMismatch ("CorrectForXZMismatch", "Correct for x/z mismatch near hard phase space boundary.", &QTildeMatching::theCorrectForXZMismatch, true, false, false); static SwitchOption interfaceCorrectForXZMismatchYes (interfaceCorrectForXZMismatch, "Yes", "Include the correction factor.", true); static SwitchOption interfaceCorrectForXZMismatchNo (interfaceCorrectForXZMismatch, "No", "Do not include the correction factor.", false); interfaceCorrectForXZMismatch.rank(-1); } diff --git a/Shower/QTilde/SplittingFunctions/SudakovFormFactor.h b/Shower/QTilde/SplittingFunctions/SudakovFormFactor.h --- a/Shower/QTilde/SplittingFunctions/SudakovFormFactor.h +++ b/Shower/QTilde/SplittingFunctions/SudakovFormFactor.h @@ -1,647 +1,652 @@ // -*- C++ -*- // // SudakovFormFactor.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_SudakovFormFactor_H #define HERWIG_SudakovFormFactor_H // // This is the declaration of the SudakovFormFactor class. // #include "ThePEG/Interface/Interfaced.h" #include "Herwig/Shower/QTilde/SplittingFunctions/SplittingFunction.h" #include "Herwig/Shower/ShowerAlpha.h" #include "Herwig/Shower/QTilde/SplittingFunctions/SplittingGenerator.fh" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/PDF/BeamParticleData.h" #include "ThePEG/EventRecord/RhoDMatrix.h" #include "ThePEG/EventRecord/SpinInfo.h" #include "Herwig/Shower/QTilde/Kinematics/ShowerKinematics.fh" #include "SudakovFormFactor.fh" #include "SudakovCutOff.h" namespace Herwig { using namespace ThePEG; /** * A typedef for the BeamParticleData */ typedef Ptr::transient_const_pointer tcBeamPtr; /** \ingroup Shower * * This is the definition of the Sudakov form factor class. In general this * is the base class for the implementation of Sudakov form factors in Herwig. * The methods generateNextTimeBranching(), generateNextDecayBranching() and * generateNextSpaceBranching need to be implemented in classes inheriting from this * one. * * In addition a number of methods are implemented to assist with the calculation * of the form factor using the veto algorithm in classes inheriting from this one. * * In general the Sudakov form-factor, for final-state radiation, is given * by * \f[\Delta_{ba}(\tilde{q}_{i+1},\tilde{q}_i)= * \exp\left\{ * -\int^{\tilde{q}^2_i}_{\tilde{q}^2_{i+1}} * \frac{{\rm d}\tilde{q}^2}{\tilde{q}^2} * \int\frac{\alpha_S(z,\tilde{q})}{2\pi} * P_{ba}(z,\tilde{q})\Theta(p_T) * \right\}. * \f] * We can solve this to obtain the next value of the scale \f$\tilde{q}_{i+1}\f$ * given the previous value \f$\tilde{q}_i\f$ * in the following way. First we obtain a simplified form of the integrand * which is greater than or equal to the true integrand for all values of * \f$\tilde{q}\f$. * * In practice it is easiest to obtain this over estimate in pieces. The ShowerAlpha * object contains an over estimate for \f$\alpha_S\f$, the splitting function * contains both an over estimate of the spltting function and its integral * which is needed to compute the over estimate of the \f$\tilde{q}\f$ integrand, * together with an over estimate of the limit of the \f$z\f$ integral. * * This gives an overestimate of the integrand * \f[g(\tilde{q}^2) = \frac{c}{\tilde{q}^2}, \f] * where because the over estimates are chosen to be independent of \f$\tilde{q}\f$ the * parameter * \f[c = \frac{\alpha_{\rm over}}{2\pi}\int^{z_1}_{z_0}P_{\rm over}(z),\f] * is a constant independent of \f$\tilde{q}\f$. * * The guesst() member can then be used to generate generate the value of * \f$\tilde{q}^2\f$ according to this result. This is done by solving the Sudakov * form factor, with the over estimates, is equal to a random number * \f$r\f$ in the interval \f$[0,1]\f$. This gives * \f[\tilde{q}^2_{i+1}=G^{-1}\left[G(\tilde{q}^2_i)+\ln r\right],\f] * where \f$G(\tilde{q}^2)=c\ln(\tilde{q}^2)\f$ is the infinite integral * of \f$g(\tilde{q}^2)\f$ and \f$G^{-1}(x)=\exp\left(\frac{x}c\right)\f$ * is its inverse. * It this case we therefore obtain * \f[\tilde{q}^2_{i+1}=\tilde{q}^2_ir^{\frac1c}.\f] * The value of \f$z\f$ can then be calculated in a similar way * \f[z = I^{-1}\left[I(z_0)+r\left(I(z_1)-I(z_0)\right)\right],\f] * using the guessz() member, * where \f$I=\int P(z){\rm d}z\f$ and \f$I^{-1}\f$ is its inverse. * * The veto algorithm then uses rejection using the ratio of the * true value to the overestimated one to obtain the original distribution. * This is accomplished using the * - alphaSVeto() member for the \f$\alpha_S\f$ veto * - SplittingFnVeto() member for the veto on the value of the splitting function. * in general there must also be a chech that the emission is in the allowed * phase space but this is left to the inheriting classes as it will depend * on the ordering variable. * * The Sudakov form factor for the initial-scale shower is different because * it must include the PDF which guides the backward evolution. * It is given by * \f[\Delta_{ba}(\tilde{q}_{i+1},\tilde{q}_i)= * \exp\left\{ * -\int^{\tilde{q}^2_i}_{\tilde{q}^2_{i+1}} * \frac{{\rm d}\tilde{q}^2}{\tilde{q}^2} * \int\frac{\alpha_S(z,\tilde{q})}{2\pi} * P_{ba}(z,\tilde{q})\frac{x'f_a(\frac{x}z,\tilde{q}^2)}{xf_b(x,\tilde{q^2})} * \right\}, * \f] * where \f$x\f$ is the fraction of the beam momentum the parton \f$b\f$ had before * the backward evolution. * This can be solve in the same way as for the final-state branching but the constant * becomes * \f[c = \frac{\alpha_{\rm over}}{2\pi}\int^{z_1}_{z_0}P_{\rm over}(z)PDF_{\rm max},\f] * where * \f[PDF_{\rm max}=\max\frac{x'f_a(\frac{x}z,\tilde{q}^2)}{xf_b(x,\tilde{q^2})},\f] * which can be set using an interface. * In addition the PDFVeto() member then is needed to implement the relevant veto. * * @see SplittingGenerator * @see SplittingFunction * @see ShowerAlpha * @see \ref SudakovFormFactorInterfaces "The interfaces" * defined for SudakovFormFactor. */ class SudakovFormFactor: public Interfaced { /** * The SplittingGenerator is a friend to insert the particles in the * branchings at initialisation */ friend class SplittingGenerator; public: /** * The default constructor. */ SudakovFormFactor() : pdfmax_(35.0), pdffactor_(0), z_( 0.0 ),phi_(0.0), pT_(){} /** * Members to generate the scale of the next branching */ //@{ /** * Return the scale of the next time-like branching. If there is no * branching then it returns ZERO. * @param startingScale starting scale for the evolution * @param ids The PDG codes of the particles in the splitting * @param enhance The radiation enhancement factor * defined. */ virtual ShoKinPtr generateNextTimeBranching(const Energy startingScale, const IdList &ids, const RhoDMatrix & rho, double enhance, double detuning, Energy2 maxQ2); /** * Return the scale of the next space-like decay branching. If there is no * branching then it returns ZERO. * @param startingScale starting scale for the evolution * @param stoppingScale stopping scale for the evolution * @param minmass The minimum mass allowed for the spake-like particle. * @param ids The PDG codes of the particles in the splitting * defined. * @param enhance The radiation enhancement factor */ virtual ShoKinPtr generateNextDecayBranching(const Energy startingScale, const Energy stoppingScale, const Energy minmass, const IdList &ids, const RhoDMatrix & rho, double enhance, double detuning); /** * Return the scale of the next space-like branching. If there is no * branching then it returns ZERO. * @param startingScale starting scale for the evolution * @param ids The PDG codes of the particles in the splitting * @param x The fraction of the beam momentum * defined. * @param beam The beam particle * @param enhance The radiation enhancement factor */ virtual ShoKinPtr generateNextSpaceBranching(const Energy startingScale, const IdList &ids,double x, const RhoDMatrix & rho, double enhance, tcBeamPtr beam, double detuning); //@} /** * Generate the azimuthal angle of the branching for forward evolution * @param particle The branching particle * @param ids The PDG codes of the particles in the branchings * @param The Shower kinematics */ virtual double generatePhiForward(ShowerParticle & particle,const IdList & ids, ShoKinPtr kinematics, const RhoDMatrix & rho); /** * Generate the azimuthal angle of the branching for backward evolution * @param particle The branching particle * @param ids The PDG codes of the particles in the branchings * @param The Shower kinematics */ virtual double generatePhiBackward(ShowerParticle & particle,const IdList & ids, ShoKinPtr kinematics, const RhoDMatrix & rho); /** * Generate the azimuthal angle of the branching for ISR in decays * @param particle The branching particle * @param ids The PDG codes of the particles in the branchings * @param The Shower kinematics */ virtual double generatePhiDecay(ShowerParticle & particle,const IdList & ids, ShoKinPtr kinematics, const RhoDMatrix & rho); /** * Methods to provide public access to the private member variables */ //@{ /** * Return the pointer to the SplittingFunction object. */ tSplittingFnPtr splittingFn() const { return splittingFn_; } /** * Return the pointer to the ShowerAlpha object. */ tShowerAlphaPtr alpha() const { return alpha_; } /** * The type of interaction */ inline ShowerInteraction interactionType() const {return splittingFn_->interactionType();} //@} public: /** * Methods to access the kinematic variables for the branching */ //@{ /** * The energy fraction */ double z() const { return z_; } /** * The azimuthal angle */ double phi() const { return phi_; } /** * The transverse momentum */ Energy pT() const { return pT_; } //@} /** * Access the maximum weight for the PDF veto */ double pdfMax() const { return pdfmax_;} /** * Method to return the evolution scale given the * transverse momentum, \f$p_T\f$ and \f$z\f$. */ virtual Energy calculateScale(double z, Energy pt, IdList ids,unsigned int iopt); /** * Method to create the ShowerKinematics object for a final-state branching */ virtual ShoKinPtr createFinalStateBranching(Energy scale,double z, double phi, Energy pt); /** * Method to create the ShowerKinematics object for an initial-state branching */ virtual ShoKinPtr createInitialStateBranching(Energy scale,double z, double phi, Energy pt); /** * Method to create the ShowerKinematics object for a decay branching */ virtual ShoKinPtr createDecayBranching(Energy scale,double z, double phi, Energy pt); public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** * Methods to provide the next value of the scale before the vetos * are applied. */ //@{ /** * Value of the energy fraction and scale for time-like branching * @param t The scale * @param tmin The minimum scale * @param enhance The radiation enhancement factor * @return False if scale less than minimum, true otherwise */ bool guessTimeLike(Energy2 &t, Energy2 tmin, double enhance, double detune); /** * Value of the energy fraction and scale for time-like branching * @param t The scale * @param tmax The maximum scale * @param minmass The minimum mass of the particle after the branching * @param enhance The radiation enhancement factor */ bool guessDecay(Energy2 &t, Energy2 tmax,Energy minmass, double enhance, double detune); /** * Value of the energy fraction and scale for space-like branching * @param t The scale * @param tmin The minimum scale * @param x Fraction of the beam momentum. * @param enhance The radiation enhancement factor */ bool guessSpaceLike(Energy2 &t, Energy2 tmin, const double x, double enhance, double detune); //@} /** * Initialize the values of the cut-offs and scales * @param tmin The minimum scale * @param ids The ids of the partics in the branching */ void initialize(const IdList & ids,Energy2 &tmin); /** * Phase Space veto member to implement the \f$\Theta\f$ function as a veto * so that the emission is within the allowed phase space. * @param t The scale * @param maxQ2 The maximum virtuality * @return true if vetoed */ bool PSVeto(const Energy2 t,const Energy2 maxQ2); /** * Compute the limits on \f$z\f$ for time-like branching * @param scale The scale of the particle * @return True if lower limit less than upper, otherwise false */ bool computeTimeLikeLimits(Energy2 & scale); /** * Compute the limits on \f$z\f$ for space-like branching * @param scale The scale of the particle * @param x The energy fraction of the parton * @return True if lower limit less than upper, otherwise false */ bool computeSpaceLikeLimits(Energy2 & scale, double x); protected: /** * Methods to implement the veto algorithm to generate the scale of * the next branching */ //@{ /** * Value of the energy fraction for the veto algorithm * @param iopt The option for calculating z * @param ids The PDG codes of the particles in the splitting * - 0 is final-state * - 1 is initial-state for the hard process * - 2 is initial-state for particle decays */ double guessz (unsigned int iopt, const IdList &ids) const; /** * Value of the scale for the veto algorithm * @param t1 The starting valoe of the scale * @param iopt The option for calculating t * @param ids The PDG codes of the particles in the splitting * - 0 is final-state * - 1 is initial-state for the hard process * - 2 is initial-state for particle decays * @param enhance The radiation enhancement factor * @param identical Whether or not the outgoing particles are identical */ Energy2 guesst (Energy2 t1,unsigned int iopt, const IdList &ids, double enhance, bool identical, double detune) const; /** * Veto on the PDF for the initial-state shower * @param t The scale * @param x The fraction of the beam momentum * @param parton0 Pointer to the particleData for the * new parent (this is the particle we evolved back to) * @param parton1 Pointer to the particleData for the * original particle * @param beam The BeamParticleData object */ bool PDFVeto(const Energy2 t, const double x, const tcPDPtr parton0, const tcPDPtr parton1, tcBeamPtr beam) const; /** * The PDF veto ratio */ double PDFVetoRatio(const Energy2 t, const double x, const tcPDPtr parton0, const tcPDPtr parton1, tcBeamPtr beam,double factor) const; /** * The veto on the splitting function. * @param t The scale * @param ids The PDG codes of the particles in the splitting * @param mass Whether or not to use the massive splitting functions * @return true if vetoed */ bool SplittingFnVeto(const Energy2 t, const IdList &ids, const bool mass, const RhoDMatrix & rho, double detune) const { return UseRandom::rnd()>SplittingFnVetoRatio(t,ids,mass,rho,detune); } /** * The Splitting function veto ratio */ double SplittingFnVetoRatio(const Energy2 t, const IdList &ids, const bool mass, const RhoDMatrix & rho, double detune) const { return splittingFn_->ratioP(z_, t, ids,mass,rho)/detune; } /** * The veto on the coupling constant * @param pt2 The value of ther transverse momentum squared, \f$p_T^2\f$. * @return true if vetoed */ bool alphaSVeto(Energy2 pt2) const; /** * The alpha S veto ratio */ double alphaSVetoRatio(Energy2 pt2,double factor) const; //@} /** * Set the particles in the splittings */ void addSplitting(const IdList &); /** * Delete the particles in the splittings */ void removeSplitting(const IdList &); /** * Access the potential branchings */ const vector & particles() const { return particles_; } public: /** * Set the PDF */ void setPDF(tcPDFPtr pdf, Energy scale) { pdf_ = pdf; freeze_ = scale; } public: /** * Calculate the virtual masses for a branchings */ const vector & virtualMasses(const IdList & ids) { return cutoff_->virtualMasses(ids); } + /** + * The minimum pT2 + */ + Energy2 pT2min() { return cutoff_->pT2min(); } + protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const {return new_ptr(*this);} /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const {return new_ptr(*this);} //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ SudakovFormFactor & operator=(const SudakovFormFactor &) = delete; private: /** * Pointer to the splitting function for this Sudakov form factor */ SplittingFnPtr splittingFn_; /** * Pointer to the coupling for this Sudakov form factor */ ShowerAlphaPtr alpha_; /** * Pointer to the coupling for this Sudakov form factor */ SudakovCutOffPtr cutoff_; /** * Maximum value of the PDF weight */ double pdfmax_; /** * List of the particles this Sudakov is used for to aid in setting up * interpolation tables if needed */ vector particles_; /** * Option for the inclusion of a factor \f$1/(1-z)\f$ in the PDF estimate */ unsigned pdffactor_; private: /** * Member variables to keep the shower kinematics information * generated by a call to generateNextTimeBranching or generateNextSpaceBranching */ //@{ /** * The energy fraction */ double z_; /** * The azimuthal angle */ double phi_; /** * The transverse momentum */ Energy pT_; //@} /** * The limits of \f$z\f$ in the splitting */ pair zlimits_; /** * Stuff for the PDFs */ //@{ /** * PDf */ tcPDFPtr pdf_; /** * Freezing scale */ Energy freeze_; //@} private: /** * The evolution scale, \f$\tilde{q}\f$. */ Energy q_; /** * The Ids of the particles in the current branching */ IdList ids_; /** * The masses of the particles in the current branching */ vector masses_; /** * The mass squared of the particles in the current branching */ vector masssquared_; }; } #endif /* HERWIG_SudakovFormFactor_H */