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,524 @@ // -*- 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" 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; } Energy Qg = theQTildeSudakov->kinScale(); Energy2 pt2 = ZERO; 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); } if ( dipole()->bornEmitter() < 2 ) { pt2 = sqr((1.-z)*qtilde) - z*sqr(Qg); } 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(); 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())); } if ( dipole()->bornEmitter() < 2 ) { return sqr((1.-z)*qtilde) - z*sqr(Qg) >= 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(); 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); } if ( dipole()->bornEmitter() < 2 ) { pt2 = sqr((1.-z)*qtilde) - z*sqr(Qg); } 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 + 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/MatrixElement/Matchbox/Matching/QTildeMatching.h b/MatrixElement/Matchbox/Matching/QTildeMatching.h --- a/MatrixElement/Matchbox/Matching/QTildeMatching.h +++ b/MatrixElement/Matchbox/Matching/QTildeMatching.h @@ -1,201 +1,201 @@ // -*- C++ -*- // // QTildeMatching.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_QTildeMatching_H #define Herwig_QTildeMatching_H // // This is the declaration of the QTildeMatching class. // #include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h" #include "Herwig/Shower/ShowerHandler.h" #include "Herwig/Shower/QTilde/Base/PartnerFinder.h" -#include "Herwig/Shower/QTilde/Default/QTildeSudakov.h" +#include "Herwig/Shower/QTilde/Base/SudakovFormFactor.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief QTildeMatching implements NLO matching with the default shower. * */ class QTildeMatching: public Herwig::ShowerApproximation { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ QTildeMatching(); /** * The destructor. */ virtual ~QTildeMatching(); //@} public: /** * Return the shower approximation to the real emission cross * section for the given pair of Born and real emission * configurations. */ virtual CrossSection dSigHatDR() const; /** * Return the shower approximation splitting kernel for the given * pair of Born and real emission configurations in units of the * Born center of mass energy squared, and including a weight to * project onto the splitting given by the dipole used. */ virtual double me2() const; /** * Determine if the configuration is below or above the cutoff. */ virtual void checkCutoff(); /** * Determine all kinematic variables which are not provided by the * dipole kinematics; store all shower variables in the respective * dipole object for later use. */ virtual void getShowerVariables(); protected: /** * Return true, if the shower was able to generate an emission * leading from the given Born to the given real emission process. */ virtual bool isInShowerPhasespace() const; /** * Return true, if the shower emission leading from the given Born * to the given real emission process would have been generated * above the shower's infrared cutoff. */ virtual bool isAboveCutoff() const; /** * Calculate qtilde^2 and z for the splitting considered */ void calculateShowerVariables() const; /** * Return the splitting function as a function of the kinematic * variables */ double splitFn(const pair&) const; 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: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** 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; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ QTildeMatching & operator=(const QTildeMatching &); /** * The shower handler to be used */ Ptr::ptr theShowerHandler; /** * The qtilde partner finder for calculating the hard scales */ Ptr::ptr theQTildeFinder; /** * The qtilde Sudakov to access the cutoff */ - Ptr::ptr theQTildeSudakov; + Ptr::ptr theQTildeSudakov; /** * True, if PDF weight should be corrected for z/x mismatch at the * hard phase space boundary */ bool theCorrectForXZMismatch; }; } #endif /* Herwig_QTildeMatching_H */ diff --git a/Shower/QTilde/Base/SudakovFormFactor.cc b/Shower/QTilde/Base/SudakovFormFactor.cc --- a/Shower/QTilde/Base/SudakovFormFactor.cc +++ b/Shower/QTilde/Base/SudakovFormFactor.cc @@ -1,364 +1,1399 @@ // -*- C++ -*- // // SudakovFormFactor.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 SudakovFormFactor class. // #include "SudakovFormFactor.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ShowerKinematics.h" #include "ShowerParticle.h" #include "ThePEG/Utilities/DescribeClass.h" -#include "Herwig/Shower/ShowerHandler.h" +#include "Herwig/Shower/QTilde/QTildeShowerHandler.h" +#include "Herwig/Shower/QTilde/Default/FS_QTildeShowerKinematics1to2.h" +#include "Herwig/Shower/QTilde/Default/IS_QTildeShowerKinematics1to2.h" +#include "Herwig/Shower/QTilde/Default/Decay_QTildeShowerKinematics1to2.h" using namespace Herwig; -DescribeAbstractClass +DescribeClass describeSudakovFormFactor ("Herwig::SudakovFormFactor",""); void SudakovFormFactor::persistentOutput(PersistentOStream & os) const { os << splittingFn_ << alpha_ << pdfmax_ << particles_ << pdffactor_ << a_ << b_ << ounit(c_,GeV) << ounit(kinCutoffScale_,GeV) << cutOffOption_ << ounit(vgcut_,GeV) << ounit(vqcut_,GeV) << ounit(pTmin_,GeV) << ounit(pT2min_,GeV2); } void SudakovFormFactor::persistentInput(PersistentIStream & is, int) { is >> splittingFn_ >> alpha_ >> pdfmax_ >> particles_ >> pdffactor_ >> a_ >> b_ >> iunit(c_,GeV) >> iunit(kinCutoffScale_,GeV) >> cutOffOption_ >> iunit(vgcut_,GeV) >> iunit(vqcut_,GeV) >> iunit(pTmin_,GeV) >> iunit(pT2min_,GeV2); } void SudakovFormFactor::Init() { static ClassDocumentation documentation ("The SudakovFormFactor class is the base class for the implementation of Sudakov" " form factors in Herwig"); static Reference interfaceSplittingFunction("SplittingFunction", "A reference to the SplittingFunction object", &Herwig::SudakovFormFactor::splittingFn_, false, false, true, false); static Reference interfaceAlpha("Alpha", "A reference to the Alpha object", &Herwig::SudakovFormFactor::alpha_, false, false, true, false); static Parameter interfacePDFmax ("PDFmax", "Maximum value of PDF weight. ", &SudakovFormFactor::pdfmax_, 35.0, 1.0, 1000000.0, false, false, Interface::limited); static Switch interfacePDFFactor ("PDFFactor", "Include additional factors in the overestimate for the PDFs", &SudakovFormFactor::pdffactor_, 0, false, false); static SwitchOption interfacePDFFactorNo (interfacePDFFactor, "No", "Don't include any factors", 0); static SwitchOption interfacePDFFactorOverZ (interfacePDFFactor, "OverZ", "Include an additional factor of 1/z", 1); static SwitchOption interfacePDFFactorOverOneMinusZ (interfacePDFFactor, "OverOneMinusZ", "Include an additional factor of 1/(1-z)", 2); static SwitchOption interfacePDFFactorOverZOneMinusZ (interfacePDFFactor, "OverZOneMinusZ", "Include an additional factor of 1/z/(1-z)", 3); static SwitchOption interfacePDFFactorOverRootZ (interfacePDFFactor, "OverRootZ", "Include an additional factor of 1/sqrt(z)", 4); static SwitchOption interfacePDFFactorRootZ (interfacePDFFactor, "RootZ", "Include an additional factor of sqrt(z)", 5); static Switch interfaceCutOffOption ("CutOffOption", "The type of cut-off to use to end the shower", &SudakovFormFactor::cutOffOption_, 0, false, false); static SwitchOption interfaceCutOffOptionDefault (interfaceCutOffOption, "Default", "Use the standard Herwig cut-off on virtualities with the minimum" " virtuality depending on the mass of the branching particle", 0); static SwitchOption interfaceCutOffOptionFORTRAN (interfaceCutOffOption, "FORTRAN", "Use a FORTRAN-like cut-off on virtualities", 1); static SwitchOption interfaceCutOffOptionpT (interfaceCutOffOption, "pT", "Use a cut on the minimum allowed pT", 2); static Parameter interfaceaParameter ("aParameter", "The a parameter for the kinematic cut-off", &SudakovFormFactor::a_, 0.3, -10.0, 10.0, false, false, Interface::limited); static Parameter interfacebParameter ("bParameter", "The b parameter for the kinematic cut-off", &SudakovFormFactor::b_, 2.3, -10.0, 10.0, false, false, Interface::limited); static Parameter interfacecParameter ("cParameter", "The c parameter for the kinematic cut-off", &SudakovFormFactor::c_, GeV, 0.3*GeV, 0.1*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfaceKinScale ("cutoffKinScale", "kinematic cutoff scale for the parton shower phase" " space (unit [GeV])", &SudakovFormFactor::kinCutoffScale_, GeV, 2.3*GeV, 0.001*GeV, 10.0*GeV,false,false,false); static Parameter interfaceGluonVirtualityCut ("GluonVirtualityCut", "For the FORTRAN cut-off option the minimum virtuality of the gluon", &SudakovFormFactor::vgcut_, GeV, 0.85*GeV, 0.1*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfaceQuarkVirtualityCut ("QuarkVirtualityCut", "For the FORTRAN cut-off option the minimum virtuality added to" " the mass for particles other than the gluon", &SudakovFormFactor::vqcut_, GeV, 0.85*GeV, 0.1*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfacepTmin ("pTmin", "The minimum pT if using a cut-off on the pT", &SudakovFormFactor::pTmin_, GeV, 1.0*GeV, ZERO, 10.0*GeV, false, false, Interface::limited); } bool SudakovFormFactor::alphaSVeto(Energy2 pt2) const { double ratio=alphaSVetoRatio(pt2,1.); return UseRandom::rnd() > ratio; } double SudakovFormFactor::alphaSVetoRatio(Energy2 pt2, double factor) const { factor *= ShowerHandler::currentHandler()->renormalizationScaleFactor(); return ThePEG::Math::powi(alpha_->ratio(pt2, factor), splittingFn_->interactionOrder()); } bool SudakovFormFactor::PDFVeto(const Energy2 t, const double x, const tcPDPtr parton0, const tcPDPtr parton1, Ptr::transient_const_pointer beam) const { double ratio=PDFVetoRatio(t,x,parton0,parton1,beam,1.); return UseRandom::rnd() > ratio; } double SudakovFormFactor::PDFVetoRatio(const Energy2 t, const double x, const tcPDPtr parton0, const tcPDPtr parton1, Ptr::transient_const_pointer beam,double factor) const { assert(pdf_); Energy2 theScale = t * sqr(ShowerHandler::currentHandler()->factorizationScaleFactor()*factor); if (theScale < sqr(freeze_)) theScale = sqr(freeze_); double newpdf(0.0), oldpdf(0.0); newpdf=pdf_->xfx(beam,parton0,theScale,x/z()); oldpdf=pdf_->xfx(beam,parton1,theScale,x); if(newpdf<=0.) return 0.; if(oldpdf<=0.) return 1.; double ratio = newpdf/oldpdf; double maxpdf = pdfmax_; switch (pdffactor_) { case 0: break; case 1: maxpdf /= z(); break; case 2: maxpdf /= 1.-z(); break; case 3: maxpdf /= (z()*(1.-z())); break; case 4: maxpdf /= sqrt(z()); break; case 5: maxpdf *= sqrt(z()); break; default : throw Exception() << "SudakovFormFactor::PDFVetoRatio invalid PDFfactor = " << pdffactor_ << Exception::runerror; } if (ratio > maxpdf) { generator()->log() << "PDFVeto warning: Ratio > " << name() << ":PDFmax (by a factor of " << ratio/maxpdf <<") for " << parton0->PDGName() << " to " << parton1->PDGName() << "\n"; } return ratio/maxpdf ; } void SudakovFormFactor::addSplitting(const IdList & in) { bool add=true; for(unsigned int ix=0;ix::iterator it=particles_.begin(); it!=particles_.end();++it) { if(it->size()==in.size()) { bool match=true; for(unsigned int iy=0;iy::iterator itemp=it; --itemp; particles_.erase(it); it = itemp; } } } } Energy2 SudakovFormFactor::guesst(Energy2 t1,unsigned int iopt, const IdList &ids, double enhance,bool ident, double detune) const { unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_; double c = 1./((splittingFn_->integOverP(zlimits_.second,ids,pdfopt) - splittingFn_->integOverP(zlimits_.first ,ids,pdfopt))* alpha_->overestimateValue()/Constants::twopi*enhance*detune); assert(iopt<=2); if(iopt==1) { c/=pdfmax_; //symmetry of FS gluon splitting if(ident) c*=0.5; } else if(iopt==2) c*=-1.; if(splittingFn_->interactionOrder()==1) { double r = UseRandom::rnd(); if(iopt!=2 || c*log(r)interactionOrder()-1); c/=Math::powi(alpha_->overestimateValue()/Constants::twopi,nm); return t1 / pow (1. - nm*c*log(UseRandom::rnd()) * Math::powi(t1*UnitRemoval::InvE2,nm) ,1./double(nm)); } } double SudakovFormFactor::guessz (unsigned int iopt, const IdList &ids) const { unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_; double lower = splittingFn_->integOverP(zlimits_.first,ids,pdfopt); return splittingFn_->invIntegOverP (lower + UseRandom::rnd()*(splittingFn_->integOverP(zlimits_.second,ids,pdfopt) - lower),ids,pdfopt); } void SudakovFormFactor::doinit() { Interfaced::doinit(); pT2min_ = cutOffOption()==2 ? sqr(pTmin_) : ZERO; } const vector & SudakovFormFactor::virtualMasses(const IdList & ids) { static vector output; output.clear(); if(cutOffOption() == 0) { for(unsigned int ix=0;ixmass()); Energy kinCutoff= kinematicCutOff(kinScale(),*std::max_element(output.begin(),output.end())); for(unsigned int ix=0;ixmass()); output.back() += ids[ix]->id()==ParticleID::g ? vgCut() : vqCut(); } } else if(cutOffOption() == 2) { for(unsigned int ix=0;ixmass()); } else { throw Exception() << "Unknown option for the cut-off" << " in SudakovFormFactor::virtualMasses()" << Exception::runerror; } return output; } + + + + + + + + + + + + + + + +bool SudakovFormFactor::guessTimeLike(Energy2 &t,Energy2 tmin,double enhance, + double detune) { + Energy2 told = t; + // calculate limits on z and if lower>upper return + if(!computeTimeLikeLimits(t)) return false; + // guess values of t and z + t = guesst(told,0,ids_,enhance,ids_[1]==ids_[2],detune); + z(guessz(0,ids_)); + // actual values for z-limits + if(!computeTimeLikeLimits(t)) return false; + if(tupper return + if(!computeSpaceLikeLimits(t,x)) return false; + // guess values of t and z + t = guesst(told,1,ids_,enhance,ids_[1]==ids_[2],detune); + z(guessz(1,ids_)); + // actual values for z-limits + if(!computeSpaceLikeLimits(t,x)) return false; + if(t zLimits().second) return true; + Energy2 q2 = z()*(1.-z())*t; + if(ids_[0]->id()!=ParticleID::g && + ids_[0]->id()!=ParticleID::gamma ) q2 += masssquared_[0]; + if(q2>maxQ2) return true; + // compute the pts + Energy2 pt2 = z()*(1.-z())*q2 - masssquared_[1]*(1.-z()) - masssquared_[2]*z(); + // if pt2<0 veto + if(pt2 min + if(tmax<=tmin) return ShoKinPtr(); + // calculate next value of t using veto algorithm + Energy2 t(tmax); + // no shower variations to calculate + if(ShowerHandler::currentHandler()->showerVariations().empty()){ + // Without variations do the usual Veto algorithm + // No need for more if-statements in this loop. + do { + if(!guessTimeLike(t,tmin,enhance,detuning)) break; + } + while(PSVeto(t,maxQ2) || + SplittingFnVeto(z()*(1.-z())*t,ids,true,rho,detuning) || + alphaSVeto(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t)); + } + else { + bool alphaRew(true),PSRew(true),SplitRew(true); + do { + if(!guessTimeLike(t,tmin,enhance,detuning)) break; + PSRew=PSVeto(t,maxQ2); + if (PSRew) continue; + SplitRew=SplittingFnVeto(z()*(1.-z())*t,ids,true,rho,detuning); + alphaRew=alphaSVeto(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t); + double factor=alphaSVetoRatio(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t,1.)* + SplittingFnVetoRatio(z()*(1.-z())*t,ids,true,rho,detuning); + + tShowerHandlerPtr ch = ShowerHandler::currentHandler(); + + if( !(SplitRew || alphaRew) ) { + //Emission + q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; + if (q_ <= ZERO) break; + } + + for ( map::const_iterator var = + ch->showerVariations().begin(); + var != ch->showerVariations().end(); ++var ) { + if ( ( ch->firstInteraction() && var->second.firstInteraction ) || + ( !ch->firstInteraction() && var->second.secondaryInteractions ) ) { + + double newfactor = alphaSVetoRatio(splittingFn()->pTScale() ? + sqr(z()*(1.-z()))*t : + z()*(1.-z())*t,var->second.renormalizationScaleFactor) + * SplittingFnVetoRatio(z()*(1.-z())*t,ids,true,rho,detuning); + + double varied; + if ( SplitRew || alphaRew ) { + // No Emission + varied = (1. - newfactor) / (1. - factor); + } else { + // Emission + varied = newfactor / factor; + } + + map::iterator wi = ch->currentWeights().find(var->first); + if ( wi != ch->currentWeights().end() ) + wi->second *= varied; + else { + assert(false); + //ch->currentWeights()[var->first] = varied; + } + } + } + + } + while(PSRew || SplitRew || alphaRew); + } + q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; + if(q_ < ZERO) return ShoKinPtr(); + + // return the ShowerKinematics object + return createFinalStateBranching(q_,z(),phi(),pT()); +} + +ShoKinPtr SudakovFormFactor:: +generateNextSpaceBranching(const Energy startingQ, + const IdList &ids, + double x, + const RhoDMatrix & rho, + double enhance, + Ptr::transient_const_pointer beam, + double detuning) { + // First reset the internal kinematics variables that can + // have been eventually set in the previous call to the method. + q_ = ZERO; + z(0.); + phi(0.); + // perform the initialization + Energy2 tmax(sqr(startingQ)),tmin; + initialize(ids,tmin); + // check max > min + if(tmax<=tmin) return ShoKinPtr(); + // calculate next value of t using veto algorithm + Energy2 t(tmax),pt2(ZERO); + // no shower variations + if(ShowerHandler::currentHandler()->showerVariations().empty()){ + // Without variations do the usual Veto algorithm + // No need for more if-statements in this loop. + do { + if(!guessSpaceLike(t,tmin,x,enhance,detuning)) break; + pt2=sqr(1.-z())*t-z()*masssquared_[2]; + } + while(pt2 < pT2min()|| + z() > zLimits().second|| + SplittingFnVeto((1.-z())*t/z(),ids,false,rho,detuning)|| + alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t)|| + PDFVeto(t,x,ids[0],ids[1],beam)); + } + // shower variations + else + { + bool alphaRew(true),PDFRew(true),ptRew(true),zRew(true),SplitRew(true); + do { + if(!guessSpaceLike(t,tmin,x,enhance,detuning)) break; + pt2=sqr(1.-z())*t-z()*masssquared_[2]; + ptRew=pt2 < pT2min(); + zRew=z() > zLimits().second; + if (ptRew||zRew) continue; + SplitRew=SplittingFnVeto((1.-z())*t/z(),ids,false,rho,detuning); + alphaRew=alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t); + PDFRew=PDFVeto(t,x,ids[0],ids[1],beam); + double factor=PDFVetoRatio(t,x,ids[0],ids[1],beam,1.)* + alphaSVetoRatio(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t,1.)* + SplittingFnVetoRatio((1.-z())*t/z(),ids,false,rho,detuning); + + tShowerHandlerPtr ch = ShowerHandler::currentHandler(); + + if( !(PDFRew || SplitRew || alphaRew) ) { + //Emission + q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; + if (q_ <= ZERO) break; + } + + for ( map::const_iterator var = + ch->showerVariations().begin(); + var != ch->showerVariations().end(); ++var ) { + if ( ( ch->firstInteraction() && var->second.firstInteraction ) || + ( !ch->firstInteraction() && var->second.secondaryInteractions ) ) { + + + + double newfactor = PDFVetoRatio(t,x,ids[0],ids[1],beam,var->second.factorizationScaleFactor)* + alphaSVetoRatio(splittingFn()->pTScale() ? + sqr(1.-z())*t : (1.-z())*t,var->second.renormalizationScaleFactor) + *SplittingFnVetoRatio((1.-z())*t/z(),ids,false,rho,detuning); + + double varied; + if( PDFRew || SplitRew || alphaRew) { + // No Emission + varied = (1. - newfactor) / (1. - factor); + } else { + // Emission + varied = newfactor / factor; + } + + + map::iterator wi = ch->currentWeights().find(var->first); + if ( wi != ch->currentWeights().end() ) + wi->second *= varied; + else { + assert(false); + //ch->currentWeights()[var->first] = varied; + } + } + } + + } + while( PDFRew || SplitRew || alphaRew); + } + if(t > ZERO && zLimits().first < zLimits().second) q_ = sqrt(t); + else return ShoKinPtr(); + + pT(sqrt(pt2)); + // create the ShowerKinematics and return it + return createInitialStateBranching(q_,z(),phi(),pT()); +} + +void SudakovFormFactor::initialize(const IdList & ids, Energy2 & tmin) { + ids_=ids; + tmin = cutOffOption() != 2 ? ZERO : 4.*pT2min(); + masses_ = virtualMasses(ids); + masssquared_.clear(); + for(unsigned int ix=0;ix0) tmin=max(masssquared_[ix],tmin); + } +} + +ShoKinPtr SudakovFormFactor::generateNextDecayBranching(const Energy startingScale, + const Energy stoppingScale, + const Energy minmass, + const IdList &ids, + const RhoDMatrix & rho, + double enhance, + double detuning) { + // First reset the internal kinematics variables that can + // have been eventually set in the previous call to this method. + q_ = Constants::MaxEnergy; + z(0.); + phi(0.); + // perform initialisation + Energy2 tmax(sqr(stoppingScale)),tmin; + initialize(ids,tmin); + tmin=sqr(startingScale); + // check some branching possible + if(tmax<=tmin) return ShoKinPtr(); + // perform the evolution + Energy2 t(tmin),pt2(-MeV2); + do { + if(!guessDecay(t,tmax,minmass,enhance,detuning)) break; + pt2 = sqr(1.-z())*(t-masssquared_[0])-z()*masssquared_[2]; + } + while(SplittingFnVeto((1.-z())*t/z(),ids,true,rho,detuning)|| + alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t ) || + pt2masssquared_[0]-sqr(minmass)); + if(t > ZERO) { + q_ = sqrt(t); + pT(sqrt(pt2)); + } + else return ShoKinPtr(); + phi(0.); + // create the ShowerKinematics object + return createDecayBranching(q_,z(),phi(),pT()); +} + +bool SudakovFormFactor::guessDecay(Energy2 &t,Energy2 tmax, Energy minmass, + double enhance, double detune) { + // previous scale + Energy2 told = t; + // overestimated limits on z + if(tmax limits=make_pair(sqr(minmass/masses_[0]), + 1.-sqrt(masssquared_[2]+pT2min()+ + 0.25*sqr(masssquared_[2])/tm2)/tm + +0.5*masssquared_[2]/tm2); + zLimits(limits); + if(zLimits().secondtmax||zLimits().second limits; + if(ids_[0]->id()==ParticleID::g||ids_[0]->id()==ParticleID::gamma) { + // no emission possible + if(t<16.*(masssquared_[1]+pT2min())) { + t=-1.*GeV2; + return false; + } + // overestimate of the limits + limits.first = 0.5*(1.-sqrt(1.-4.*sqrt((masssquared_[1]+pT2min())/t))); + limits.second = 1.-limits.first; + } + // special case for radiated particle is gluon + else if(ids_[2]->id()==ParticleID::g||ids_[2]->id()==ParticleID::gamma) { + limits.first = sqrt((masssquared_[1]+pT2min())/t); + limits.second = 1.-sqrt((masssquared_[2]+pT2min())/t); + } + else if(ids_[1]->id()==ParticleID::g||ids_[1]->id()==ParticleID::gamma) { + limits.second = sqrt((masssquared_[2]+pT2min())/t); + limits.first = 1.-sqrt((masssquared_[1]+pT2min())/t); + } + else { + limits.first = (masssquared_[1]+pT2min())/t; + limits.second = 1.-(masssquared_[2]+pT2min())/t; + } + if(limits.first>=limits.second) { + t=-1.*GeV2; + return false; + } + zLimits(limits); + return true; +} + +bool SudakovFormFactor::computeSpaceLikeLimits(Energy2 & t, double x) { + if (t < 1e-20 * GeV2) { + t=-1.*GeV2; + return false; + } + pair limits; + // compute the limits + limits.first = x; + double yy = 1.+0.5*masssquared_[2]/t; + limits.second = yy - sqrt(sqr(yy)-1.+pT2min()/t); + // return false if lower>upper + zLimits(limits); + if(limits.second(particle.parents()[0]) : tShowerParticlePtr(); + } + else { + mother = particle.children().size()==2 ? + dynamic_ptr_cast(&particle) : tShowerParticlePtr(); + } + tShowerParticlePtr partner; + while(mother) { + tPPtr otherChild; + if(forward) { + for (unsigned int ix=0;ixchildren().size();++ix) { + if(mother->children()[ix]!=child) { + otherChild = mother->children()[ix]; + break; + } + } + } + else { + otherChild = mother->children()[1]; + } + tShowerParticlePtr other = dynamic_ptr_cast(otherChild); + if((inter==ShowerInteraction::QCD && otherChild->dataPtr()->coloured()) || + (inter==ShowerInteraction::QED && otherChild->dataPtr()->charged())) { + partner = other; + break; + } + if(forward && !other->isFinalState()) { + partner = dynamic_ptr_cast(mother); + break; + } + child = mother; + if(forward) { + mother = ! mother->parents().empty() ? + dynamic_ptr_cast(mother->parents()[0]) : tShowerParticlePtr(); + } + else { + if(mother->children()[0]->children().size()!=2) + break; + tShowerParticlePtr mtemp = + dynamic_ptr_cast(mother->children()[0]); + if(!mtemp) + break; + else + mother=mtemp; + } + } + if(!partner) { + if(forward) { + partner = dynamic_ptr_cast( child)->partner(); + } + else { + if(mother) { + tShowerParticlePtr parent; + if(!mother->children().empty()) { + parent = dynamic_ptr_cast(mother->children()[0]); + } + if(!parent) { + parent = dynamic_ptr_cast(mother); + } + partner = parent->partner(); + } + else { + partner = dynamic_ptr_cast(&particle)->partner(); + } + } + } + return partner; +} + +pair softPhiMin(double phi0, double phi1, double A, double B, double C, double D) { + double c01 = cos(phi0 - phi1); + double s01 = sin(phi0 - phi1); + double s012(sqr(s01)), c012(sqr(c01)); + double A2(A*A), B2(B*B), C2(C*C), D2(D*D); + if(abs(B/A)<1e-10 && abs(D/C)<1e-10) return make_pair(phi0,phi0+Constants::pi); + double root = sqr(B2)*C2*D2*sqr(s012) + 2.*A*B2*B*C2*C*D*c01*s012 + 2.*A*B2*B*C*D2*D*c01*s012 + + 4.*A2*B2*C2*D2*c012 - A2*B2*C2*D2*s012 - A2*B2*sqr(D2)*s012 - sqr(B2)*sqr(C2)*s012 + - sqr(B2)*C2*D2*s012 - 4.*A2*A*B*C*D2*D*c01 - 4.*A*B2*B*C2*C*D*c01 + sqr(A2)*sqr(D2) + + 2.*A2*B2*C2*D2 + sqr(B2)*sqr(C2); + if(root<0.) return make_pair(phi0,phi0+Constants::pi); + root = sqrt(root); + double denom = (-2.*A*B*C*D*c01 + A2*D2 + B2*C2); + double denom2 = (-B*C*c01 + A*D); + if(denom==ZERO || denom2==0) + return make_pair(phi0,phi0+Constants::pi); + double num = B2*C*D*s012; + return make_pair(atan2(B*s01*(-C*(num + root) / denom + D) / denom2, -(num + root ) / denom) + phi0, + atan2(B*s01*(-C*(num - root) / denom + D) / denom2, -(num - root ) / denom) + phi0); +} + +} + +double SudakovFormFactor::generatePhiForward(ShowerParticle & particle, + const IdList & ids, + ShoKinPtr kinematics, + const RhoDMatrix & rho) { + // no correlations, return flat phi + if(! dynamic_ptr_cast(ShowerHandler::currentHandler())->correlations()) + return Constants::twopi*UseRandom::rnd(); + // get the kinematic variables + double z = kinematics->z(); + Energy2 t = z*(1.-z)*sqr(kinematics->scale()); + Energy pT = kinematics->pT(); + // if soft correlations + Energy2 pipj,pik; + bool canBeSoft[2] = {ids[1]->id()==ParticleID::g || ids[1]->id()==ParticleID::gamma, + ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma }; + vector pjk(3,ZERO); + vector Ek(3,ZERO); + Energy Ei,Ej; + Energy2 m12(ZERO),m22(ZERO); + InvEnergy2 aziMax(ZERO); + bool softAllowed = dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()&& + (canBeSoft[0] || canBeSoft[1]); + if(softAllowed) { + // find the partner for the soft correlations + tShowerParticlePtr partner=findCorrelationPartner(particle,true,splittingFn()->interactionType()); + // remember we want the softer gluon + bool swapOrder = !canBeSoft[1] || (canBeSoft[0] && canBeSoft[1] && z < 0.5); + double zFact = !swapOrder ? (1.-z) : z; + // compute the transforms to the shower reference frame + // first the boost + Lorentz5Momentum pVect = particle.showerBasis()->pVector(); + Lorentz5Momentum nVect = particle.showerBasis()->nVector(); + Boost beta_bb; + if(particle.showerBasis()->frame()==ShowerBasis::BackToBack) { + beta_bb = -(pVect + nVect).boostVector(); + } + else if(particle.showerBasis()->frame()==ShowerBasis::Rest) { + beta_bb = -pVect.boostVector(); + } + else + assert(false); + pVect.boost(beta_bb); + nVect.boost(beta_bb); + Axis axis; + if(particle.showerBasis()->frame()==ShowerBasis::BackToBack) { + axis = pVect.vect().unit(); + } + else if(particle.showerBasis()->frame()==ShowerBasis::Rest) { + axis = nVect.vect().unit(); + } + else + assert(false); + // and then the rotation + LorentzRotation rot; + if(axis.perp2()>0.) { + double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); + rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); + } + else if(axis.z()<0.) { + rot.rotate(Constants::pi,Axis(1.,0.,0.)); + } + rot.invert(); + pVect *= rot; + nVect *= rot; + // shower parameters + Energy2 pn = pVect*nVect, m2 = pVect.m2(); + double alpha0 = particle.showerParameters().alpha; + double beta0 = 0.5/alpha0/pn* + (sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt)); + Lorentz5Momentum qperp0(particle.showerParameters().ptx, + particle.showerParameters().pty,ZERO,ZERO); + assert(partner); + Lorentz5Momentum pj = partner->momentum(); + pj.boost(beta_bb); + pj *= rot; + // compute the two phi independent dot products + pik = 0.5*zFact*(sqr(alpha0)*m2 - sqr(particle.showerParameters().pt) + 2.*alpha0*beta0*pn ) + +0.5*sqr(pT)/zFact; + Energy2 dot1 = pj*pVect; + Energy2 dot2 = pj*nVect; + Energy2 dot3 = pj*qperp0; + pipj = alpha0*dot1+beta0*dot2+dot3; + // compute the constants for the phi dependent dot product + pjk[0] = zFact*(alpha0*dot1+dot3-0.5*dot2/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) + +0.5*sqr(pT)*dot2/pn/zFact/alpha0; + pjk[1] = (pj.x() - dot2/alpha0/pn*qperp0.x())*pT; + pjk[2] = (pj.y() - dot2/alpha0/pn*qperp0.y())*pT; + m12 = sqr(particle.dataPtr()->mass()); + m22 = sqr(partner->dataPtr()->mass()); + if(swapOrder) { + pjk[1] *= -1.; + pjk[2] *= -1.; + } + Ek[0] = zFact*(alpha0*pVect.t()-0.5*nVect.t()/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) + +0.5*sqr(pT)*nVect.t()/pn/zFact/alpha0; + Ek[1] = -nVect.t()/alpha0/pn*qperp0.x()*pT; + Ek[2] = -nVect.t()/alpha0/pn*qperp0.y()*pT; + if(swapOrder) { + Ek[1] *= -1.; + Ek[2] *= -1.; + } + Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2])); + Ei = alpha0*pVect.t()+beta0*nVect.t(); + Ej = pj.t(); + double phi0 = atan2(-pjk[2],-pjk[1]); + if(phi0<0.) phi0 += Constants::twopi; + double phi1 = atan2(-Ek[2],-Ek[1]); + if(phi1<0.) phi1 += Constants::twopi; + double xi_min = pik/Ei/(Ek[0]+mag2), xi_max = pik/Ei/(Ek[0]-mag2), xi_ij = pipj/Ei/Ej; + if(xi_min>xi_max) swap(xi_min,xi_max); + if(xi_min>xi_ij) softAllowed = false; + Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); + if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { + aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); + } + else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { + double A = (pipj*Ek[0]- Ej*pik)/Ej/sqr(Ej); + double B = -sqrt(sqr(pipj)*(sqr(Ek[1])+sqr(Ek[2])))/Ej/sqr(Ej); + double C = pjk[0]/sqr(Ej); + double D = -sqrt(sqr(pjk[1])+sqr(pjk[2]))/sqr(Ej); + pair minima = softPhiMin(phi0,phi1,A,B,C,D); + aziMax = 0.5/pik/(Ek[0]-mag2)*(Ei-m12*(Ek[0]-mag2)/pik + max(Ej*(A+B*cos(minima.first -phi1))/(C+D*cos(minima.first -phi0)), + Ej*(A+B*cos(minima.second-phi1))/(C+D*cos(minima.second-phi0)))); + } + else + assert(false); + } + // if spin correlations + vector > wgts; + if(dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations()) { + // calculate the weights + wgts = splittingFn()->generatePhiForward(z,t,ids,rho); + } + else { + wgts = vector >(1,make_pair(0,1.)); + } + // generate the azimuthal angle + double phi,wgt; + static const Complex ii(0.,1.); + unsigned int ntry(0); + double phiMax(0.),wgtMax(0.); + do { + phi = Constants::twopi*UseRandom::rnd(); + // first the spin correlations bit (gives 1 if correlations off) + Complex spinWgt = 0.; + for(unsigned int ix=0;ix1e-10) { + generator()->log() << "Forward spin weight problem " << wgt << " " << wgt-1. + << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; + generator()->log() << "Weights \n"; + for(unsigned int ix=0;ixlog() << wgts[ix].first << " " << wgts[ix].second << "\n"; + } + // soft correlations bit + double aziWgt = 1.; + if(softAllowed) { + Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); + Energy Eg = Ek[0]+Ek[1]*cos(phi)+Ek[2]*sin(phi); + if(pipj*Eg>pik*Ej) { + if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { + aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; + } + else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { + aziWgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax); + } + if(aziWgt-1.>1e-10||aziWgt<-1e-10) { + generator()->log() << "Forward soft weight problem " << aziWgt << " " << aziWgt-1. + << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; + } + } + else { + aziWgt = 0.; + } + } + wgt *= aziWgt; + if(wgt>wgtMax) { + phiMax = phi; + wgtMax = wgt; + } + ++ntry; + } + while(wgtlog() << "Too many tries to generate phi in forward evolution\n"; + phi = phiMax; + } + // return the azimuthal angle + return phi; +} + +double SudakovFormFactor::generatePhiBackward(ShowerParticle & particle, + const IdList & ids, + ShoKinPtr kinematics, + const RhoDMatrix & rho) { + // no correlations, return flat phi + if(! dynamic_ptr_cast(ShowerHandler::currentHandler())->correlations()) + return Constants::twopi*UseRandom::rnd(); + // get the kinematic variables + double z = kinematics->z(); + Energy2 t = (1.-z)*sqr(kinematics->scale())/z; + Energy pT = kinematics->pT(); + // if soft correlations + bool softAllowed = dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations() && + (ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma); + Energy2 pipj,pik,m12(ZERO),m22(ZERO); + vector pjk(3,ZERO); + Energy Ei,Ej,Ek; + InvEnergy2 aziMax(ZERO); + if(softAllowed) { + // find the partner for the soft correlations + tShowerParticlePtr partner=findCorrelationPartner(particle,false,splittingFn()->interactionType()); + double zFact = (1.-z); + // compute the transforms to the shower reference frame + // first the boost + Lorentz5Momentum pVect = particle.showerBasis()->pVector(); + Lorentz5Momentum nVect = particle.showerBasis()->nVector(); + assert(particle.showerBasis()->frame()==ShowerBasis::BackToBack); + Boost beta_bb = -(pVect + nVect).boostVector(); + pVect.boost(beta_bb); + nVect.boost(beta_bb); + Axis axis = pVect.vect().unit(); + // and then the rotation + LorentzRotation rot; + if(axis.perp2()>0.) { + double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); + rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); + } + else if(axis.z()<0.) { + rot.rotate(Constants::pi,Axis(1.,0.,0.)); + } + rot.invert(); + pVect *= rot; + nVect *= rot; + // shower parameters + Energy2 pn = pVect*nVect; + Energy2 m2 = pVect.m2(); + double alpha0 = particle.x(); + double beta0 = -0.5/alpha0/pn*sqr(alpha0)*m2; + Lorentz5Momentum pj = partner->momentum(); + pj.boost(beta_bb); + pj *= rot; + double beta2 = 0.5*(1.-zFact)*(sqr(alpha0*zFact/(1.-zFact))*m2+sqr(pT))/alpha0/zFact/pn; + // compute the two phi independent dot products + Energy2 dot1 = pj*pVect; + Energy2 dot2 = pj*nVect; + pipj = alpha0*dot1+beta0*dot2; + pik = alpha0*(alpha0*zFact/(1.-zFact)*m2+pn*(beta2+zFact/(1.-zFact)*beta0)); + // compute the constants for the phi dependent dot product + pjk[0] = alpha0*zFact/(1.-zFact)*dot1+beta2*dot2; + pjk[1] = pj.x()*pT; + pjk[2] = pj.y()*pT; + m12 = ZERO; + m22 = sqr(partner->dataPtr()->mass()); + Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); + if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { + aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); + } + else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { + Ek = alpha0*zFact/(1.-zFact)*pVect.t()+beta2*nVect.t(); + Ei = alpha0*pVect.t()+beta0*nVect.t(); + Ej = pj.t(); + if(pipj*Ek> Ej*pik) { + aziMax = 0.5/pik/Ek*(Ei-m12*Ek/pik + (pipj*Ek- Ej*pik)/(pjk[0]-mag)); + } + else { + aziMax = 0.5/pik/Ek*(Ei-m12*Ek/pik); + } + } + else { + assert(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==0); + } + } + // if spin correlations + vector > wgts; + if(dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations()) { + // get the weights + wgts = splittingFn()->generatePhiBackward(z,t,ids,rho); + } + else { + wgts = vector >(1,make_pair(0,1.)); + } + // generate the azimuthal angle + double phi,wgt; + static const Complex ii(0.,1.); + unsigned int ntry(0); + double phiMax(0.),wgtMax(0.); + do { + phi = Constants::twopi*UseRandom::rnd(); + Complex spinWgt = 0.; + for(unsigned int ix=0;ix1e-10) { + generator()->log() << "Backward weight problem " << wgt << " " << wgt-1. + << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << z << " " << phi << "\n"; + generator()->log() << "Weights \n"; + for(unsigned int ix=0;ixlog() << wgts[ix].first << " " << wgts[ix].second << "\n"; + } + // soft correlations bit + double aziWgt = 1.; + if(softAllowed) { + Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); + if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { + aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; + } + else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { + aziWgt = max(ZERO,0.5/pik/Ek*(Ei-m12*Ek/pik + pipj*Ek/dot - Ej*pik/dot)/aziMax); + } + if(aziWgt-1.>1e-10||aziWgt<-1e-10) { + generator()->log() << "Backward soft weight problem " << aziWgt << " " << aziWgt-1. + << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; + } + } + wgt *= aziWgt; + if(wgt>wgtMax) { + phiMax = phi; + wgtMax = wgt; + } + ++ntry; + } + while(wgtlog() << "Too many tries to generate phi in backward evolution\n"; + phi = phiMax; + } + // return the azimuthal angle + return phi; +} + +double SudakovFormFactor::generatePhiDecay(ShowerParticle & particle, + const IdList & ids, + ShoKinPtr kinematics, + const RhoDMatrix &) { + // only soft correlations in this case + // no correlations, return flat phi + if( !(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations() && + (ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma ))) + return Constants::twopi*UseRandom::rnd(); + // get the kinematic variables + double z = kinematics->z(); + Energy pT = kinematics->pT(); + // if soft correlations + // find the partner for the soft correlations + tShowerParticlePtr partner = findCorrelationPartner(particle,true,splittingFn()->interactionType()); + double zFact(1.-z); + // compute the transforms to the shower reference frame + // first the boost + Lorentz5Momentum pVect = particle.showerBasis()->pVector(); + Lorentz5Momentum nVect = particle.showerBasis()->nVector(); + assert(particle.showerBasis()->frame()==ShowerBasis::Rest); + Boost beta_bb = -pVect.boostVector(); + pVect.boost(beta_bb); + nVect.boost(beta_bb); + Axis axis = nVect.vect().unit(); + // and then the rotation + LorentzRotation rot; + if(axis.perp2()>0.) { + double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); + rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); + } + else if(axis.z()<0.) { + rot.rotate(Constants::pi,Axis(1.,0.,0.)); + } + rot.invert(); + pVect *= rot; + nVect *= rot; + // shower parameters + Energy2 pn = pVect*nVect; + Energy2 m2 = pVect.m2(); + double alpha0 = particle.showerParameters().alpha; + double beta0 = 0.5/alpha0/pn* + (sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt)); + Lorentz5Momentum qperp0(particle.showerParameters().ptx, + particle.showerParameters().pty,ZERO,ZERO); + Lorentz5Momentum pj = partner->momentum(); + pj.boost(beta_bb); + pj *= rot; + // compute the two phi independent dot products + Energy2 pik = 0.5*zFact*(sqr(alpha0)*m2 - sqr(particle.showerParameters().pt) + 2.*alpha0*beta0*pn ) + +0.5*sqr(pT)/zFact; + Energy2 dot1 = pj*pVect; + Energy2 dot2 = pj*nVect; + Energy2 dot3 = pj*qperp0; + Energy2 pipj = alpha0*dot1+beta0*dot2+dot3; + // compute the constants for the phi dependent dot product + vector pjk(3,ZERO); + pjk[0] = zFact*(alpha0*dot1+dot3-0.5*dot2/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) + +0.5*sqr(pT)*dot2/pn/zFact/alpha0; + pjk[1] = (pj.x() - dot2/alpha0/pn*qperp0.x())*pT; + pjk[2] = (pj.y() - dot2/alpha0/pn*qperp0.y())*pT; + Energy2 m12 = sqr(particle.dataPtr()->mass()); + Energy2 m22 = sqr(partner->dataPtr()->mass()); + Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); + InvEnergy2 aziMax; + vector Ek(3,ZERO); + Energy Ei,Ej; + if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { + aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); + } + else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { + Ek[0] = zFact*(alpha0*pVect.t()+-0.5*nVect.t()/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) + +0.5*sqr(pT)*nVect.t()/pn/zFact/alpha0; + Ek[1] = -nVect.t()/alpha0/pn*qperp0.x()*pT; + Ek[2] = -nVect.t()/alpha0/pn*qperp0.y()*pT; + Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2])); + Ei = alpha0*pVect.t()+beta0*nVect.t(); + Ej = pj.t(); + aziMax = 0.5/pik/(Ek[0]-mag2)*(Ei-m12*(Ek[0]-mag2)/pik + pipj*(Ek[0]+mag2)/(pjk[0]-mag) - Ej*pik/(pjk[0]-mag) ); + } + else + assert(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==0); + // generate the azimuthal angle + double phi,wgt(0.); + unsigned int ntry(0); + double phiMax(0.),wgtMax(0.); + do { + phi = Constants::twopi*UseRandom::rnd(); + Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); + if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { + wgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; + } + else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { + if(qperp0.m2()==ZERO) { + wgt = 1.; + } + else { + Energy Eg = Ek[0]+Ek[1]*cos(phi)+Ek[2]*sin(phi); + wgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax); + } + } + if(wgt-1.>1e-10||wgt<-1e-10) { + generator()->log() << "Decay soft weight problem " << wgt << " " << wgt-1. + << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; + } + if(wgt>wgtMax) { + phiMax = phi; + wgtMax = wgt; + } + ++ntry; + } + while(wgtlog() << "Too many tries to generate phi\n"; + } + // return the azimuthal angle + return phi; +} + + +Energy SudakovFormFactor::calculateScale(double zin, Energy pt, IdList ids, + unsigned int iopt) { + Energy2 tmin; + initialize(ids,tmin); + // final-state branching + if(iopt==0) { + Energy2 scale=(sqr(pt)+masssquared_[1]*(1.-zin)+masssquared_[2]*zin); + if(ids[0]->id()!=ParticleID::g) scale -= zin*(1.-zin)*masssquared_[0]; + scale /= sqr(zin*(1-zin)); + return scale<=ZERO ? sqrt(tmin) : sqrt(scale); + } + else if(iopt==1) { + Energy2 scale=(sqr(pt)+zin*masssquared_[2])/sqr(1.-zin); + return scale<=ZERO ? sqrt(tmin) : sqrt(scale); + } + else if(iopt==2) { + Energy2 scale = (sqr(pt)+zin*masssquared_[2])/sqr(1.-zin)+masssquared_[0]; + return scale<=ZERO ? sqrt(tmin) : sqrt(scale); + } + else { + throw Exception() << "Unknown option in SudakovFormFactor::calculateScale() " + << "iopt = " << iopt << Exception::runerror; + } +} + +ShoKinPtr SudakovFormFactor::createFinalStateBranching(Energy scale,double z, + double phi, Energy pt) { + ShoKinPtr showerKin = new_ptr(FS_QTildeShowerKinematics1to2()); + showerKin->scale(scale); + showerKin->z(z); + showerKin->phi(phi); + showerKin->pT(pt); + showerKin->SudakovFormFactor(this); + return showerKin; +} + +ShoKinPtr SudakovFormFactor::createInitialStateBranching(Energy scale,double z, + double phi, Energy pt) { + ShoKinPtr showerKin = new_ptr(IS_QTildeShowerKinematics1to2()); + showerKin->scale(scale); + showerKin->z(z); + showerKin->phi(phi); + showerKin->pT(pt); + showerKin->SudakovFormFactor(this); + return showerKin; +} + +ShoKinPtr SudakovFormFactor::createDecayBranching(Energy scale,double z, + double phi, Energy pt) { + ShoKinPtr showerKin = new_ptr(Decay_QTildeShowerKinematics1to2()); + showerKin->scale(scale); + showerKin->z(z); + showerKin->phi(phi); + showerKin->pT(pt); + showerKin->SudakovFormFactor(this); + return showerKin; +} diff --git a/Shower/QTilde/Base/SudakovFormFactor.h b/Shower/QTilde/Base/SudakovFormFactor.h --- a/Shower/QTilde/Base/SudakovFormFactor.h +++ b/Shower/QTilde/Base/SudakovFormFactor.h @@ -1,688 +1,793 @@ // -*- 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 "ShowerKinematics.fh" #include "SudakovFormFactor.fh" 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), cutOffOption_(0), a_(0.3), b_(2.3), c_(0.3*GeV), kinCutoffScale_( 2.3*GeV ), vgcut_(0.85*GeV), vqcut_(0.85*GeV), pTmin_(1.*GeV), pT2min_(ZERO), 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)=0; + 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)=0; + 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)=0; + 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)=0; + 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)=0; + 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)=0; + 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)=0; + 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)=0; + 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)=0; + 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)=0; + 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: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} 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; //@} /** * Methods to set the kinematic variables for the branching */ //@{ /** * The energy fraction */ void z(double in) { z_=in; } /** * The azimuthal angle */ void phi(double in) { phi_=in; } /** * The transverse momentum */ void pT(Energy in) { pT_=in; } //@} /** * Set/Get the limits on the energy fraction for the splitting */ //@{ /** * Get the limits */ pair zLimits() const { return zlimits_;} /** * Set the limits */ void zLimits(pair in) { zlimits_=in; } //@} /** * 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: /** * @name Methods for the cut-off */ //@{ /** * The option being used */ unsigned int cutOffOption() const { return cutOffOption_; } /** * The kinematic scale */ Energy kinScale() const {return kinCutoffScale_;} /** * The virtuality cut-off on the gluon \f$Q_g=\frac{\delta-am_q}{b}\f$ * @param scale The scale \f$\delta\f$ * @param mq The quark mass \f$m_q\f$. */ Energy kinematicCutOff(Energy scale, Energy mq) const {return max((scale -a_*mq)/b_,c_);} /** * The virtualilty cut-off for gluons */ Energy vgCut() const { return vgcut_; } /** * The virtuality cut-off for everything else */ Energy vqCut() const { return vqcut_; } /** * The minimum \f$p_T\f$ for the branching */ Energy pTmin() const { return pTmin_; } /** * The square of the minimum \f$p_T\f$ */ Energy2 pT2min() const { return pT2min_; } /** * Calculate the virtual masses for a branchings */ const vector & virtualMasses(const IdList & ids); //@} /** * Set the PDF */ void setPDF(tcPDFPtr pdf, Energy scale) { pdf_ = pdf; freeze_ = scale; } +protected: + + /** @name Clone Methods. */ + //@{ + /** + * Make a simple clone of this object. + * @return a pointer to the new object. + */ + inline 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. + */ + inline 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 &); + 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_; /** * 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: /** * Option for the type of cut-off to be applied */ unsigned int cutOffOption_; /** * Parameters for the default Herwig cut-off option, i.e. the parameters for * the \f$Q_g=\max(\frac{\delta-am_q}{b},c)\f$ kinematic cut-off */ //@{ /** * The \f$a\f$ parameter */ double a_; /** * The \f$b\f$ parameter */ double b_; /** * The \f$c\f$ parameter */ Energy c_; /** * Kinematic cutoff used in the parton shower phase space. */ Energy kinCutoffScale_; //@} /** * Parameters for the FORTRAN-like cut-off */ //@{ /** * The virtualilty cut-off for gluons */ Energy vgcut_; /** * The virtuality cut-off for everything else */ Energy vqcut_; //@} /** * Parameters for the \f$p_T\f$ cut-off */ //@{ /** * The minimum \f$p_T\f$ for the branching */ Energy pTmin_; /** * The square of the minimum \f$p_T\f$ */ Energy2 pT2min_; //@} 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 */ diff --git a/Shower/QTilde/Default/QTildeSudakov.cc b/Shower/QTilde/Default/QTildeSudakov.cc deleted file mode 100644 --- a/Shower/QTilde/Default/QTildeSudakov.cc +++ /dev/null @@ -1,1061 +0,0 @@ -// -*- C++ -*- -// -// QTildeSudakov.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 QTildeSudakov class. -// - -#include "QTildeSudakov.h" -#include "ThePEG/Interface/ClassDocumentation.h" -#include "ThePEG/Interface/Parameter.h" -#include "ThePEG/Interface/Switch.h" -#include "ThePEG/PDT/ParticleData.h" -#include "ThePEG/EventRecord/Event.h" -#include "ThePEG/Repository/EventGenerator.h" -#include "ThePEG/Repository/CurrentGenerator.h" -#include "ThePEG/PDT/EnumParticles.h" -#include "Herwig/Shower/QTilde/Default/FS_QTildeShowerKinematics1to2.h" -#include "Herwig/Shower/QTilde/Default/IS_QTildeShowerKinematics1to2.h" -#include "Herwig/Shower/QTilde/Default/Decay_QTildeShowerKinematics1to2.h" -#include "ThePEG/Utilities/DescribeClass.h" -#include "Herwig/Shower/QTilde/Base/ShowerVertex.h" -#include "Herwig/Shower/QTilde/Base/ShowerParticle.h" -#include "Herwig/Shower/QTilde/QTildeShowerHandler.h" -#include "Herwig/Shower/QTilde/Base/PartnerFinder.h" -#include "Herwig/Shower/QTilde/Base/KinematicsReconstructor.h" - -using namespace Herwig; - -DescribeNoPIOClass -describeQTildeSudakov ("Herwig::QTildeSudakov","HwShower.so"); - -void QTildeSudakov::Init() { - - static ClassDocumentation documentation - ("The QTildeSudakov class implements the Sudakov form factor for ordering it" - " qtilde"); -} - -bool QTildeSudakov::guessTimeLike(Energy2 &t,Energy2 tmin,double enhance, - double detune) { - Energy2 told = t; - // calculate limits on z and if lower>upper return - if(!computeTimeLikeLimits(t)) return false; - // guess values of t and z - t = guesst(told,0,ids_,enhance,ids_[1]==ids_[2],detune); - z(guessz(0,ids_)); - // actual values for z-limits - if(!computeTimeLikeLimits(t)) return false; - if(tupper return - if(!computeSpaceLikeLimits(t,x)) return false; - // guess values of t and z - t = guesst(told,1,ids_,enhance,ids_[1]==ids_[2],detune); - z(guessz(1,ids_)); - // actual values for z-limits - if(!computeSpaceLikeLimits(t,x)) return false; - if(t zLimits().second) return true; - Energy2 q2 = z()*(1.-z())*t; - if(ids_[0]->id()!=ParticleID::g && - ids_[0]->id()!=ParticleID::gamma ) q2 += masssquared_[0]; - if(q2>maxQ2) return true; - // compute the pts - Energy2 pt2 = z()*(1.-z())*q2 - masssquared_[1]*(1.-z()) - masssquared_[2]*z(); - // if pt2<0 veto - if(pt2 min - if(tmax<=tmin) return ShoKinPtr(); - // calculate next value of t using veto algorithm - Energy2 t(tmax); - // no shower variations to calculate - if(ShowerHandler::currentHandler()->showerVariations().empty()){ - // Without variations do the usual Veto algorithm - // No need for more if-statements in this loop. - do { - if(!guessTimeLike(t,tmin,enhance,detuning)) break; - } - while(PSVeto(t,maxQ2) || - SplittingFnVeto(z()*(1.-z())*t,ids,true,rho,detuning) || - alphaSVeto(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t)); - } - else { - bool alphaRew(true),PSRew(true),SplitRew(true); - do { - if(!guessTimeLike(t,tmin,enhance,detuning)) break; - PSRew=PSVeto(t,maxQ2); - if (PSRew) continue; - SplitRew=SplittingFnVeto(z()*(1.-z())*t,ids,true,rho,detuning); - alphaRew=alphaSVeto(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t); - double factor=alphaSVetoRatio(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t,1.)* - SplittingFnVetoRatio(z()*(1.-z())*t,ids,true,rho,detuning); - - tShowerHandlerPtr ch = ShowerHandler::currentHandler(); - - if( !(SplitRew || alphaRew) ) { - //Emission - q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; - if (q_ <= ZERO) break; - } - - for ( map::const_iterator var = - ch->showerVariations().begin(); - var != ch->showerVariations().end(); ++var ) { - if ( ( ch->firstInteraction() && var->second.firstInteraction ) || - ( !ch->firstInteraction() && var->second.secondaryInteractions ) ) { - - double newfactor = alphaSVetoRatio(splittingFn()->pTScale() ? - sqr(z()*(1.-z()))*t : - z()*(1.-z())*t,var->second.renormalizationScaleFactor) - * SplittingFnVetoRatio(z()*(1.-z())*t,ids,true,rho,detuning); - - double varied; - if ( SplitRew || alphaRew ) { - // No Emission - varied = (1. - newfactor) / (1. - factor); - } else { - // Emission - varied = newfactor / factor; - } - - map::iterator wi = ch->currentWeights().find(var->first); - if ( wi != ch->currentWeights().end() ) - wi->second *= varied; - else { - assert(false); - //ch->currentWeights()[var->first] = varied; - } - } - } - - } - while(PSRew || SplitRew || alphaRew); - } - q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; - if(q_ < ZERO) return ShoKinPtr(); - - // return the ShowerKinematics object - return createFinalStateBranching(q_,z(),phi(),pT()); -} - -ShoKinPtr QTildeSudakov:: -generateNextSpaceBranching(const Energy startingQ, - const IdList &ids, - double x, - const RhoDMatrix & rho, - double enhance, - Ptr::transient_const_pointer beam, - double detuning) { - // First reset the internal kinematics variables that can - // have been eventually set in the previous call to the method. - q_ = ZERO; - z(0.); - phi(0.); - // perform the initialization - Energy2 tmax(sqr(startingQ)),tmin; - initialize(ids,tmin); - // check max > min - if(tmax<=tmin) return ShoKinPtr(); - // calculate next value of t using veto algorithm - Energy2 t(tmax),pt2(ZERO); - // no shower variations - if(ShowerHandler::currentHandler()->showerVariations().empty()){ - // Without variations do the usual Veto algorithm - // No need for more if-statements in this loop. - do { - if(!guessSpaceLike(t,tmin,x,enhance,detuning)) break; - pt2=sqr(1.-z())*t-z()*masssquared_[2]; - } - while(pt2 < pT2min()|| - z() > zLimits().second|| - SplittingFnVeto((1.-z())*t/z(),ids,false,rho,detuning)|| - alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t)|| - PDFVeto(t,x,ids[0],ids[1],beam)); - } - // shower variations - else - { - bool alphaRew(true),PDFRew(true),ptRew(true),zRew(true),SplitRew(true); - do { - if(!guessSpaceLike(t,tmin,x,enhance,detuning)) break; - pt2=sqr(1.-z())*t-z()*masssquared_[2]; - ptRew=pt2 < pT2min(); - zRew=z() > zLimits().second; - if (ptRew||zRew) continue; - SplitRew=SplittingFnVeto((1.-z())*t/z(),ids,false,rho,detuning); - alphaRew=alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t); - PDFRew=PDFVeto(t,x,ids[0],ids[1],beam); - double factor=PDFVetoRatio(t,x,ids[0],ids[1],beam,1.)* - alphaSVetoRatio(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t,1.)* - SplittingFnVetoRatio((1.-z())*t/z(),ids,false,rho,detuning); - - tShowerHandlerPtr ch = ShowerHandler::currentHandler(); - - if( !(PDFRew || SplitRew || alphaRew) ) { - //Emission - q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; - if (q_ <= ZERO) break; - } - - for ( map::const_iterator var = - ch->showerVariations().begin(); - var != ch->showerVariations().end(); ++var ) { - if ( ( ch->firstInteraction() && var->second.firstInteraction ) || - ( !ch->firstInteraction() && var->second.secondaryInteractions ) ) { - - - - double newfactor = PDFVetoRatio(t,x,ids[0],ids[1],beam,var->second.factorizationScaleFactor)* - alphaSVetoRatio(splittingFn()->pTScale() ? - sqr(1.-z())*t : (1.-z())*t,var->second.renormalizationScaleFactor) - *SplittingFnVetoRatio((1.-z())*t/z(),ids,false,rho,detuning); - - double varied; - if( PDFRew || SplitRew || alphaRew) { - // No Emission - varied = (1. - newfactor) / (1. - factor); - } else { - // Emission - varied = newfactor / factor; - } - - - map::iterator wi = ch->currentWeights().find(var->first); - if ( wi != ch->currentWeights().end() ) - wi->second *= varied; - else { - assert(false); - //ch->currentWeights()[var->first] = varied; - } - } - } - - } - while( PDFRew || SplitRew || alphaRew); - } - if(t > ZERO && zLimits().first < zLimits().second) q_ = sqrt(t); - else return ShoKinPtr(); - - pT(sqrt(pt2)); - // create the ShowerKinematics and return it - return createInitialStateBranching(q_,z(),phi(),pT()); -} - -void QTildeSudakov::initialize(const IdList & ids, Energy2 & tmin) { - ids_=ids; - tmin = cutOffOption() != 2 ? ZERO : 4.*pT2min(); - masses_ = virtualMasses(ids); - masssquared_.clear(); - for(unsigned int ix=0;ix0) tmin=max(masssquared_[ix],tmin); - } -} - -ShoKinPtr QTildeSudakov::generateNextDecayBranching(const Energy startingScale, - const Energy stoppingScale, - const Energy minmass, - const IdList &ids, - const RhoDMatrix & rho, - double enhance, - double detuning) { - // First reset the internal kinematics variables that can - // have been eventually set in the previous call to this method. - q_ = Constants::MaxEnergy; - z(0.); - phi(0.); - // perform initialisation - Energy2 tmax(sqr(stoppingScale)),tmin; - initialize(ids,tmin); - tmin=sqr(startingScale); - // check some branching possible - if(tmax<=tmin) return ShoKinPtr(); - // perform the evolution - Energy2 t(tmin),pt2(-MeV2); - do { - if(!guessDecay(t,tmax,minmass,enhance,detuning)) break; - pt2 = sqr(1.-z())*(t-masssquared_[0])-z()*masssquared_[2]; - } - while(SplittingFnVeto((1.-z())*t/z(),ids,true,rho,detuning)|| - alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t ) || - pt2masssquared_[0]-sqr(minmass)); - if(t > ZERO) { - q_ = sqrt(t); - pT(sqrt(pt2)); - } - else return ShoKinPtr(); - phi(0.); - // create the ShowerKinematics object - return createDecayBranching(q_,z(),phi(),pT()); -} - -bool QTildeSudakov::guessDecay(Energy2 &t,Energy2 tmax, Energy minmass, - double enhance, double detune) { - // previous scale - Energy2 told = t; - // overestimated limits on z - if(tmax limits=make_pair(sqr(minmass/masses_[0]), - 1.-sqrt(masssquared_[2]+pT2min()+ - 0.25*sqr(masssquared_[2])/tm2)/tm - +0.5*masssquared_[2]/tm2); - zLimits(limits); - if(zLimits().secondtmax||zLimits().second limits; - if(ids_[0]->id()==ParticleID::g||ids_[0]->id()==ParticleID::gamma) { - // no emission possible - if(t<16.*(masssquared_[1]+pT2min())) { - t=-1.*GeV2; - return false; - } - // overestimate of the limits - limits.first = 0.5*(1.-sqrt(1.-4.*sqrt((masssquared_[1]+pT2min())/t))); - limits.second = 1.-limits.first; - } - // special case for radiated particle is gluon - else if(ids_[2]->id()==ParticleID::g||ids_[2]->id()==ParticleID::gamma) { - limits.first = sqrt((masssquared_[1]+pT2min())/t); - limits.second = 1.-sqrt((masssquared_[2]+pT2min())/t); - } - else if(ids_[1]->id()==ParticleID::g||ids_[1]->id()==ParticleID::gamma) { - limits.second = sqrt((masssquared_[2]+pT2min())/t); - limits.first = 1.-sqrt((masssquared_[1]+pT2min())/t); - } - else { - limits.first = (masssquared_[1]+pT2min())/t; - limits.second = 1.-(masssquared_[2]+pT2min())/t; - } - if(limits.first>=limits.second) { - t=-1.*GeV2; - return false; - } - zLimits(limits); - return true; -} - -bool QTildeSudakov::computeSpaceLikeLimits(Energy2 & t, double x) { - if (t < 1e-20 * GeV2) { - t=-1.*GeV2; - return false; - } - pair limits; - // compute the limits - limits.first = x; - double yy = 1.+0.5*masssquared_[2]/t; - limits.second = yy - sqrt(sqr(yy)-1.+pT2min()/t); - // return false if lower>upper - zLimits(limits); - if(limits.second(particle.parents()[0]) : tShowerParticlePtr(); - } - else { - mother = particle.children().size()==2 ? - dynamic_ptr_cast(&particle) : tShowerParticlePtr(); - } - tShowerParticlePtr partner; - while(mother) { - tPPtr otherChild; - if(forward) { - for (unsigned int ix=0;ixchildren().size();++ix) { - if(mother->children()[ix]!=child) { - otherChild = mother->children()[ix]; - break; - } - } - } - else { - otherChild = mother->children()[1]; - } - tShowerParticlePtr other = dynamic_ptr_cast(otherChild); - if((inter==ShowerInteraction::QCD && otherChild->dataPtr()->coloured()) || - (inter==ShowerInteraction::QED && otherChild->dataPtr()->charged())) { - partner = other; - break; - } - if(forward && !other->isFinalState()) { - partner = dynamic_ptr_cast(mother); - break; - } - child = mother; - if(forward) { - mother = ! mother->parents().empty() ? - dynamic_ptr_cast(mother->parents()[0]) : tShowerParticlePtr(); - } - else { - if(mother->children()[0]->children().size()!=2) - break; - tShowerParticlePtr mtemp = - dynamic_ptr_cast(mother->children()[0]); - if(!mtemp) - break; - else - mother=mtemp; - } - } - if(!partner) { - if(forward) { - partner = dynamic_ptr_cast( child)->partner(); - } - else { - if(mother) { - tShowerParticlePtr parent; - if(!mother->children().empty()) { - parent = dynamic_ptr_cast(mother->children()[0]); - } - if(!parent) { - parent = dynamic_ptr_cast(mother); - } - partner = parent->partner(); - } - else { - partner = dynamic_ptr_cast(&particle)->partner(); - } - } - } - return partner; -} - -pair softPhiMin(double phi0, double phi1, double A, double B, double C, double D) { - double c01 = cos(phi0 - phi1); - double s01 = sin(phi0 - phi1); - double s012(sqr(s01)), c012(sqr(c01)); - double A2(A*A), B2(B*B), C2(C*C), D2(D*D); - if(abs(B/A)<1e-10 && abs(D/C)<1e-10) return make_pair(phi0,phi0+Constants::pi); - double root = sqr(B2)*C2*D2*sqr(s012) + 2.*A*B2*B*C2*C*D*c01*s012 + 2.*A*B2*B*C*D2*D*c01*s012 - + 4.*A2*B2*C2*D2*c012 - A2*B2*C2*D2*s012 - A2*B2*sqr(D2)*s012 - sqr(B2)*sqr(C2)*s012 - - sqr(B2)*C2*D2*s012 - 4.*A2*A*B*C*D2*D*c01 - 4.*A*B2*B*C2*C*D*c01 + sqr(A2)*sqr(D2) - + 2.*A2*B2*C2*D2 + sqr(B2)*sqr(C2); - if(root<0.) return make_pair(phi0,phi0+Constants::pi); - root = sqrt(root); - double denom = (-2.*A*B*C*D*c01 + A2*D2 + B2*C2); - double denom2 = (-B*C*c01 + A*D); - if(denom==ZERO || denom2==0) - return make_pair(phi0,phi0+Constants::pi); - double num = B2*C*D*s012; - return make_pair(atan2(B*s01*(-C*(num + root) / denom + D) / denom2, -(num + root ) / denom) + phi0, - atan2(B*s01*(-C*(num - root) / denom + D) / denom2, -(num - root ) / denom) + phi0); -} - -} - -double QTildeSudakov::generatePhiForward(ShowerParticle & particle, - const IdList & ids, - ShoKinPtr kinematics, - const RhoDMatrix & rho) { - // no correlations, return flat phi - if(! dynamic_ptr_cast(ShowerHandler::currentHandler())->correlations()) - return Constants::twopi*UseRandom::rnd(); - // get the kinematic variables - double z = kinematics->z(); - Energy2 t = z*(1.-z)*sqr(kinematics->scale()); - Energy pT = kinematics->pT(); - // if soft correlations - Energy2 pipj,pik; - bool canBeSoft[2] = {ids[1]->id()==ParticleID::g || ids[1]->id()==ParticleID::gamma, - ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma }; - vector pjk(3,ZERO); - vector Ek(3,ZERO); - Energy Ei,Ej; - Energy2 m12(ZERO),m22(ZERO); - InvEnergy2 aziMax(ZERO); - bool softAllowed = dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()&& - (canBeSoft[0] || canBeSoft[1]); - if(softAllowed) { - // find the partner for the soft correlations - tShowerParticlePtr partner=findCorrelationPartner(particle,true,splittingFn()->interactionType()); - // remember we want the softer gluon - bool swapOrder = !canBeSoft[1] || (canBeSoft[0] && canBeSoft[1] && z < 0.5); - double zFact = !swapOrder ? (1.-z) : z; - // compute the transforms to the shower reference frame - // first the boost - Lorentz5Momentum pVect = particle.showerBasis()->pVector(); - Lorentz5Momentum nVect = particle.showerBasis()->nVector(); - Boost beta_bb; - if(particle.showerBasis()->frame()==ShowerBasis::BackToBack) { - beta_bb = -(pVect + nVect).boostVector(); - } - else if(particle.showerBasis()->frame()==ShowerBasis::Rest) { - beta_bb = -pVect.boostVector(); - } - else - assert(false); - pVect.boost(beta_bb); - nVect.boost(beta_bb); - Axis axis; - if(particle.showerBasis()->frame()==ShowerBasis::BackToBack) { - axis = pVect.vect().unit(); - } - else if(particle.showerBasis()->frame()==ShowerBasis::Rest) { - axis = nVect.vect().unit(); - } - else - assert(false); - // and then the rotation - LorentzRotation rot; - if(axis.perp2()>0.) { - double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); - rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); - } - else if(axis.z()<0.) { - rot.rotate(Constants::pi,Axis(1.,0.,0.)); - } - rot.invert(); - pVect *= rot; - nVect *= rot; - // shower parameters - Energy2 pn = pVect*nVect, m2 = pVect.m2(); - double alpha0 = particle.showerParameters().alpha; - double beta0 = 0.5/alpha0/pn* - (sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt)); - Lorentz5Momentum qperp0(particle.showerParameters().ptx, - particle.showerParameters().pty,ZERO,ZERO); - assert(partner); - Lorentz5Momentum pj = partner->momentum(); - pj.boost(beta_bb); - pj *= rot; - // compute the two phi independent dot products - pik = 0.5*zFact*(sqr(alpha0)*m2 - sqr(particle.showerParameters().pt) + 2.*alpha0*beta0*pn ) - +0.5*sqr(pT)/zFact; - Energy2 dot1 = pj*pVect; - Energy2 dot2 = pj*nVect; - Energy2 dot3 = pj*qperp0; - pipj = alpha0*dot1+beta0*dot2+dot3; - // compute the constants for the phi dependent dot product - pjk[0] = zFact*(alpha0*dot1+dot3-0.5*dot2/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) - +0.5*sqr(pT)*dot2/pn/zFact/alpha0; - pjk[1] = (pj.x() - dot2/alpha0/pn*qperp0.x())*pT; - pjk[2] = (pj.y() - dot2/alpha0/pn*qperp0.y())*pT; - m12 = sqr(particle.dataPtr()->mass()); - m22 = sqr(partner->dataPtr()->mass()); - if(swapOrder) { - pjk[1] *= -1.; - pjk[2] *= -1.; - } - Ek[0] = zFact*(alpha0*pVect.t()-0.5*nVect.t()/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) - +0.5*sqr(pT)*nVect.t()/pn/zFact/alpha0; - Ek[1] = -nVect.t()/alpha0/pn*qperp0.x()*pT; - Ek[2] = -nVect.t()/alpha0/pn*qperp0.y()*pT; - if(swapOrder) { - Ek[1] *= -1.; - Ek[2] *= -1.; - } - Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2])); - Ei = alpha0*pVect.t()+beta0*nVect.t(); - Ej = pj.t(); - double phi0 = atan2(-pjk[2],-pjk[1]); - if(phi0<0.) phi0 += Constants::twopi; - double phi1 = atan2(-Ek[2],-Ek[1]); - if(phi1<0.) phi1 += Constants::twopi; - double xi_min = pik/Ei/(Ek[0]+mag2), xi_max = pik/Ei/(Ek[0]-mag2), xi_ij = pipj/Ei/Ej; - if(xi_min>xi_max) swap(xi_min,xi_max); - if(xi_min>xi_ij) softAllowed = false; - Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); - if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { - aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); - } - else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { - double A = (pipj*Ek[0]- Ej*pik)/Ej/sqr(Ej); - double B = -sqrt(sqr(pipj)*(sqr(Ek[1])+sqr(Ek[2])))/Ej/sqr(Ej); - double C = pjk[0]/sqr(Ej); - double D = -sqrt(sqr(pjk[1])+sqr(pjk[2]))/sqr(Ej); - pair minima = softPhiMin(phi0,phi1,A,B,C,D); - aziMax = 0.5/pik/(Ek[0]-mag2)*(Ei-m12*(Ek[0]-mag2)/pik + max(Ej*(A+B*cos(minima.first -phi1))/(C+D*cos(minima.first -phi0)), - Ej*(A+B*cos(minima.second-phi1))/(C+D*cos(minima.second-phi0)))); - } - else - assert(false); - } - // if spin correlations - vector > wgts; - if(dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations()) { - // calculate the weights - wgts = splittingFn()->generatePhiForward(z,t,ids,rho); - } - else { - wgts = vector >(1,make_pair(0,1.)); - } - // generate the azimuthal angle - double phi,wgt; - static const Complex ii(0.,1.); - unsigned int ntry(0); - double phiMax(0.),wgtMax(0.); - do { - phi = Constants::twopi*UseRandom::rnd(); - // first the spin correlations bit (gives 1 if correlations off) - Complex spinWgt = 0.; - for(unsigned int ix=0;ix1e-10) { - generator()->log() << "Forward spin weight problem " << wgt << " " << wgt-1. - << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; - generator()->log() << "Weights \n"; - for(unsigned int ix=0;ixlog() << wgts[ix].first << " " << wgts[ix].second << "\n"; - } - // soft correlations bit - double aziWgt = 1.; - if(softAllowed) { - Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); - Energy Eg = Ek[0]+Ek[1]*cos(phi)+Ek[2]*sin(phi); - if(pipj*Eg>pik*Ej) { - if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { - aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; - } - else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { - aziWgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax); - } - if(aziWgt-1.>1e-10||aziWgt<-1e-10) { - generator()->log() << "Forward soft weight problem " << aziWgt << " " << aziWgt-1. - << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; - } - } - else { - aziWgt = 0.; - } - } - wgt *= aziWgt; - if(wgt>wgtMax) { - phiMax = phi; - wgtMax = wgt; - } - ++ntry; - } - while(wgtlog() << "Too many tries to generate phi in forward evolution\n"; - phi = phiMax; - } - // return the azimuthal angle - return phi; -} - -double QTildeSudakov::generatePhiBackward(ShowerParticle & particle, - const IdList & ids, - ShoKinPtr kinematics, - const RhoDMatrix & rho) { - // no correlations, return flat phi - if(! dynamic_ptr_cast(ShowerHandler::currentHandler())->correlations()) - return Constants::twopi*UseRandom::rnd(); - // get the kinematic variables - double z = kinematics->z(); - Energy2 t = (1.-z)*sqr(kinematics->scale())/z; - Energy pT = kinematics->pT(); - // if soft correlations - bool softAllowed = dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations() && - (ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma); - Energy2 pipj,pik,m12(ZERO),m22(ZERO); - vector pjk(3,ZERO); - Energy Ei,Ej,Ek; - InvEnergy2 aziMax(ZERO); - if(softAllowed) { - // find the partner for the soft correlations - tShowerParticlePtr partner=findCorrelationPartner(particle,false,splittingFn()->interactionType()); - double zFact = (1.-z); - // compute the transforms to the shower reference frame - // first the boost - Lorentz5Momentum pVect = particle.showerBasis()->pVector(); - Lorentz5Momentum nVect = particle.showerBasis()->nVector(); - assert(particle.showerBasis()->frame()==ShowerBasis::BackToBack); - Boost beta_bb = -(pVect + nVect).boostVector(); - pVect.boost(beta_bb); - nVect.boost(beta_bb); - Axis axis = pVect.vect().unit(); - // and then the rotation - LorentzRotation rot; - if(axis.perp2()>0.) { - double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); - rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); - } - else if(axis.z()<0.) { - rot.rotate(Constants::pi,Axis(1.,0.,0.)); - } - rot.invert(); - pVect *= rot; - nVect *= rot; - // shower parameters - Energy2 pn = pVect*nVect; - Energy2 m2 = pVect.m2(); - double alpha0 = particle.x(); - double beta0 = -0.5/alpha0/pn*sqr(alpha0)*m2; - Lorentz5Momentum pj = partner->momentum(); - pj.boost(beta_bb); - pj *= rot; - double beta2 = 0.5*(1.-zFact)*(sqr(alpha0*zFact/(1.-zFact))*m2+sqr(pT))/alpha0/zFact/pn; - // compute the two phi independent dot products - Energy2 dot1 = pj*pVect; - Energy2 dot2 = pj*nVect; - pipj = alpha0*dot1+beta0*dot2; - pik = alpha0*(alpha0*zFact/(1.-zFact)*m2+pn*(beta2+zFact/(1.-zFact)*beta0)); - // compute the constants for the phi dependent dot product - pjk[0] = alpha0*zFact/(1.-zFact)*dot1+beta2*dot2; - pjk[1] = pj.x()*pT; - pjk[2] = pj.y()*pT; - m12 = ZERO; - m22 = sqr(partner->dataPtr()->mass()); - Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); - if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { - aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); - } - else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { - Ek = alpha0*zFact/(1.-zFact)*pVect.t()+beta2*nVect.t(); - Ei = alpha0*pVect.t()+beta0*nVect.t(); - Ej = pj.t(); - if(pipj*Ek> Ej*pik) { - aziMax = 0.5/pik/Ek*(Ei-m12*Ek/pik + (pipj*Ek- Ej*pik)/(pjk[0]-mag)); - } - else { - aziMax = 0.5/pik/Ek*(Ei-m12*Ek/pik); - } - } - else { - assert(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==0); - } - } - // if spin correlations - vector > wgts; - if(dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations()) { - // get the weights - wgts = splittingFn()->generatePhiBackward(z,t,ids,rho); - } - else { - wgts = vector >(1,make_pair(0,1.)); - } - // generate the azimuthal angle - double phi,wgt; - static const Complex ii(0.,1.); - unsigned int ntry(0); - double phiMax(0.),wgtMax(0.); - do { - phi = Constants::twopi*UseRandom::rnd(); - Complex spinWgt = 0.; - for(unsigned int ix=0;ix1e-10) { - generator()->log() << "Backward weight problem " << wgt << " " << wgt-1. - << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << z << " " << phi << "\n"; - generator()->log() << "Weights \n"; - for(unsigned int ix=0;ixlog() << wgts[ix].first << " " << wgts[ix].second << "\n"; - } - // soft correlations bit - double aziWgt = 1.; - if(softAllowed) { - Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); - if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { - aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; - } - else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { - aziWgt = max(ZERO,0.5/pik/Ek*(Ei-m12*Ek/pik + pipj*Ek/dot - Ej*pik/dot)/aziMax); - } - if(aziWgt-1.>1e-10||aziWgt<-1e-10) { - generator()->log() << "Backward soft weight problem " << aziWgt << " " << aziWgt-1. - << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; - } - } - wgt *= aziWgt; - if(wgt>wgtMax) { - phiMax = phi; - wgtMax = wgt; - } - ++ntry; - } - while(wgtlog() << "Too many tries to generate phi in backward evolution\n"; - phi = phiMax; - } - // return the azimuthal angle - return phi; -} - -double QTildeSudakov::generatePhiDecay(ShowerParticle & particle, - const IdList & ids, - ShoKinPtr kinematics, - const RhoDMatrix &) { - // only soft correlations in this case - // no correlations, return flat phi - if( !(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations() && - (ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma ))) - return Constants::twopi*UseRandom::rnd(); - // get the kinematic variables - double z = kinematics->z(); - Energy pT = kinematics->pT(); - // if soft correlations - // find the partner for the soft correlations - tShowerParticlePtr partner = findCorrelationPartner(particle,true,splittingFn()->interactionType()); - double zFact(1.-z); - // compute the transforms to the shower reference frame - // first the boost - Lorentz5Momentum pVect = particle.showerBasis()->pVector(); - Lorentz5Momentum nVect = particle.showerBasis()->nVector(); - assert(particle.showerBasis()->frame()==ShowerBasis::Rest); - Boost beta_bb = -pVect.boostVector(); - pVect.boost(beta_bb); - nVect.boost(beta_bb); - Axis axis = nVect.vect().unit(); - // and then the rotation - LorentzRotation rot; - if(axis.perp2()>0.) { - double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); - rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); - } - else if(axis.z()<0.) { - rot.rotate(Constants::pi,Axis(1.,0.,0.)); - } - rot.invert(); - pVect *= rot; - nVect *= rot; - // shower parameters - Energy2 pn = pVect*nVect; - Energy2 m2 = pVect.m2(); - double alpha0 = particle.showerParameters().alpha; - double beta0 = 0.5/alpha0/pn* - (sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt)); - Lorentz5Momentum qperp0(particle.showerParameters().ptx, - particle.showerParameters().pty,ZERO,ZERO); - Lorentz5Momentum pj = partner->momentum(); - pj.boost(beta_bb); - pj *= rot; - // compute the two phi independent dot products - Energy2 pik = 0.5*zFact*(sqr(alpha0)*m2 - sqr(particle.showerParameters().pt) + 2.*alpha0*beta0*pn ) - +0.5*sqr(pT)/zFact; - Energy2 dot1 = pj*pVect; - Energy2 dot2 = pj*nVect; - Energy2 dot3 = pj*qperp0; - Energy2 pipj = alpha0*dot1+beta0*dot2+dot3; - // compute the constants for the phi dependent dot product - vector pjk(3,ZERO); - pjk[0] = zFact*(alpha0*dot1+dot3-0.5*dot2/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) - +0.5*sqr(pT)*dot2/pn/zFact/alpha0; - pjk[1] = (pj.x() - dot2/alpha0/pn*qperp0.x())*pT; - pjk[2] = (pj.y() - dot2/alpha0/pn*qperp0.y())*pT; - Energy2 m12 = sqr(particle.dataPtr()->mass()); - Energy2 m22 = sqr(partner->dataPtr()->mass()); - Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); - InvEnergy2 aziMax; - vector Ek(3,ZERO); - Energy Ei,Ej; - if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { - aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); - } - else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { - Ek[0] = zFact*(alpha0*pVect.t()+-0.5*nVect.t()/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) - +0.5*sqr(pT)*nVect.t()/pn/zFact/alpha0; - Ek[1] = -nVect.t()/alpha0/pn*qperp0.x()*pT; - Ek[2] = -nVect.t()/alpha0/pn*qperp0.y()*pT; - Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2])); - Ei = alpha0*pVect.t()+beta0*nVect.t(); - Ej = pj.t(); - aziMax = 0.5/pik/(Ek[0]-mag2)*(Ei-m12*(Ek[0]-mag2)/pik + pipj*(Ek[0]+mag2)/(pjk[0]-mag) - Ej*pik/(pjk[0]-mag) ); - } - else - assert(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==0); - // generate the azimuthal angle - double phi,wgt(0.); - unsigned int ntry(0); - double phiMax(0.),wgtMax(0.); - do { - phi = Constants::twopi*UseRandom::rnd(); - Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); - if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { - wgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; - } - else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { - if(qperp0.m2()==ZERO) { - wgt = 1.; - } - else { - Energy Eg = Ek[0]+Ek[1]*cos(phi)+Ek[2]*sin(phi); - wgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax); - } - } - if(wgt-1.>1e-10||wgt<-1e-10) { - generator()->log() << "Decay soft weight problem " << wgt << " " << wgt-1. - << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; - } - if(wgt>wgtMax) { - phiMax = phi; - wgtMax = wgt; - } - ++ntry; - } - while(wgtlog() << "Too many tries to generate phi\n"; - } - // return the azimuthal angle - return phi; -} - - -Energy QTildeSudakov::calculateScale(double zin, Energy pt, IdList ids, - unsigned int iopt) { - Energy2 tmin; - initialize(ids,tmin); - // final-state branching - if(iopt==0) { - Energy2 scale=(sqr(pt)+masssquared_[1]*(1.-zin)+masssquared_[2]*zin); - if(ids[0]->id()!=ParticleID::g) scale -= zin*(1.-zin)*masssquared_[0]; - scale /= sqr(zin*(1-zin)); - return scale<=ZERO ? sqrt(tmin) : sqrt(scale); - } - else if(iopt==1) { - Energy2 scale=(sqr(pt)+zin*masssquared_[2])/sqr(1.-zin); - return scale<=ZERO ? sqrt(tmin) : sqrt(scale); - } - else if(iopt==2) { - Energy2 scale = (sqr(pt)+zin*masssquared_[2])/sqr(1.-zin)+masssquared_[0]; - return scale<=ZERO ? sqrt(tmin) : sqrt(scale); - } - else { - throw Exception() << "Unknown option in QTildeSudakov::calculateScale() " - << "iopt = " << iopt << Exception::runerror; - } -} - -ShoKinPtr QTildeSudakov::createFinalStateBranching(Energy scale,double z, - double phi, Energy pt) { - ShoKinPtr showerKin = new_ptr(FS_QTildeShowerKinematics1to2()); - showerKin->scale(scale); - showerKin->z(z); - showerKin->phi(phi); - showerKin->pT(pt); - showerKin->SudakovFormFactor(this); - return showerKin; -} - -ShoKinPtr QTildeSudakov::createInitialStateBranching(Energy scale,double z, - double phi, Energy pt) { - ShoKinPtr showerKin = new_ptr(IS_QTildeShowerKinematics1to2()); - showerKin->scale(scale); - showerKin->z(z); - showerKin->phi(phi); - showerKin->pT(pt); - showerKin->SudakovFormFactor(this); - return showerKin; -} - -ShoKinPtr QTildeSudakov::createDecayBranching(Energy scale,double z, - double phi, Energy pt) { - ShoKinPtr showerKin = new_ptr(Decay_QTildeShowerKinematics1to2()); - showerKin->scale(scale); - showerKin->z(z); - showerKin->phi(phi); - showerKin->pT(pt); - showerKin->SudakovFormFactor(this); - return showerKin; -} diff --git a/Shower/QTilde/Default/QTildeSudakov.h b/Shower/QTilde/Default/QTildeSudakov.h deleted file mode 100644 --- a/Shower/QTilde/Default/QTildeSudakov.h +++ /dev/null @@ -1,291 +0,0 @@ -// -*- C++ -*- -// -// QTildeSudakov.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_QTildeSudakov_H -#define HERWIG_QTildeSudakov_H -// -// This is the declaration of the QTildeSudakov class. -// - -#include "Herwig/Shower/QTilde/Base/SudakovFormFactor.h" - -namespace Herwig { - -using namespace ThePEG; - -/** \ingroup Shower - * - * The QTildeSudakov class implements the Sudakov form factor for evolution in - * \f$\tilde{q}^2\f$ using the veto algorithm. - * - * @see \ref QTildeSudakovInterfaces "The interfaces" - * defined for QTildeSudakov. - */ -class QTildeSudakov: public SudakovFormFactor { - -public: - - /** - * The default constructor. - */ - inline QTildeSudakov() {} - - /** - * 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 - * @param maxQ2 The maximum \f$Q^2\f$ for the emission - */ - 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 enhance The radiation enhancement factor - * @param beam The beam particle - */ - 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 branching - * @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 branching - * @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); - - /** - * 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: - - /** @name Clone Methods. */ - //@{ - /** - * Make a simple clone of this object. - * @return a pointer to the new object. - */ - inline 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. - */ - inline 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. - */ - QTildeSudakov & operator=(const QTildeSudakov &); - -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_QTildeSudakov_H */ diff --git a/Shower/QTilde/Makefile.am b/Shower/QTilde/Makefile.am --- a/Shower/QTilde/Makefile.am +++ b/Shower/QTilde/Makefile.am @@ -1,42 +1,41 @@ SUBDIRS = Matching pkglib_LTLIBRARIES = HwShower.la HwShower_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 25:0:0 HwShower_la_SOURCES = \ Couplings/ShowerAlphaQCD.h Couplings/ShowerAlphaQCD.cc \ Couplings/ShowerAlphaQED.h Couplings/ShowerAlphaQED.cc\ QTildeShowerHandler.h QTildeShowerHandler.fh QTildeShowerHandler.cc \ SplittingFunctions/HalfHalfOneSplitFn.h SplittingFunctions/HalfHalfOneSplitFn.cc\ SplittingFunctions/OneOneOneSplitFn.h SplittingFunctions/OneOneOneSplitFn.cc\ SplittingFunctions/OneOneOneMassiveSplitFn.h SplittingFunctions/OneOneOneMassiveSplitFn.cc\ SplittingFunctions/ZeroZeroOneSplitFn.h SplittingFunctions/ZeroZeroOneSplitFn.cc\ SplittingFunctions/OneHalfHalfSplitFn.h SplittingFunctions/OneHalfHalfSplitFn.cc\ SplittingFunctions/HalfOneHalfSplitFn.h SplittingFunctions/HalfOneHalfSplitFn.cc\ SplittingFunctions/CMWOneOneOneSplitFn.h SplittingFunctions/CMWOneOneOneSplitFn.cc\ SplittingFunctions/CMWHalfHalfOneSplitFn.h SplittingFunctions/CMWHalfHalfOneSplitFn.cc\ -Default/QTildeSudakov.cc Default/QTildeSudakov.h\ Default/Decay_QTildeShowerKinematics1to2.cc \ Default/Decay_QTildeShowerKinematics1to2.h \ Default/IS_QTildeShowerKinematics1to2.cc Default/IS_QTildeShowerKinematics1to2.h \ Default/FS_QTildeShowerKinematics1to2.cc Default/FS_QTildeShowerKinematics1to2.h \ Base/KinematicsReconstructor.cc \ Base/KinematicsReconstructor.tcc \ Base/KinematicsReconstructor.h \ Base/KinematicsReconstructor.fh \ Base/HardTree.cc Base/HardTree.h Base/HardTree.fh \ Base/HardBranching.h Base/HardBranching.fh Base/HardBranching.cc\ Base/PartnerFinder.h Base/PartnerFinder.fh Base/PartnerFinder.cc \ Base/ShowerVeto.h Base/ShowerVeto.fh Base/ShowerVeto.cc \ Base/FullShowerVeto.h Base/FullShowerVeto.fh Base/FullShowerVeto.cc \ SplittingFunctions/SplittingGenerator.cc SplittingFunctions/SplittingGenerator.h\ SplittingFunctions/SplittingGenerator.fh \ Base/ShowerTree.h Base/ShowerTree.fh Base/ShowerTree.cc \ ShowerConfig.h ShowerConfig.cc \ Base/Branching.h \ Base/ShowerParticle.cc Base/ShowerParticle.fh Base/ShowerParticle.h \ Base/ShowerKinematics.fh Base/ShowerKinematics.h Base/ShowerKinematics.cc \ Base/ShowerBasis.fh Base/ShowerBasis.h Base/ShowerBasis.cc \ Base/ShowerProgenitor.fh Base/ShowerProgenitor.h \ Base/SudakovFormFactor.cc Base/SudakovFormFactor.h Base/SudakovFormFactor.fh \ SplittingFunctions/SplittingFunction.h SplittingFunctions/SplittingFunction.fh \ SplittingFunctions/SplittingFunction.cc \ Base/ShowerVertex.cc Base/ShowerVertex.fh Base/ShowerVertex.h diff --git a/src/defaults/Shower.in b/src/defaults/Shower.in --- a/src/defaults/Shower.in +++ b/src/defaults/Shower.in @@ -1,310 +1,310 @@ # -*- ThePEG-repository -*- ############################################################ # Setup of default parton shower # # Useful switches for users are marked near the top of # this file. # # Don't edit this file directly, but reset the switches # in your own input files! ############################################################ library HwMPI.so library HwShower.so library HwMatching.so mkdir /Herwig/Shower cd /Herwig/Shower create Herwig::QTildeShowerHandler ShowerHandler newdef ShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler newdef ShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer # use LO PDFs for Shower, can be changed later newdef ShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF newdef ShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF newdef ShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF newdef ShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF ##################################### # initial setup, don't change these! ##################################### create Herwig::SplittingGenerator SplittingGenerator create Herwig::ShowerAlphaQCD AlphaQCD create Herwig::ShowerAlphaQED AlphaQED create Herwig::PartnerFinder PartnerFinder newdef PartnerFinder:PartnerMethod 1 newdef PartnerFinder:ScaleChoice 1 create Herwig::KinematicsReconstructor KinematicsReconstructor newdef KinematicsReconstructor:ReconstructionOption Colour3 newdef KinematicsReconstructor:InitialStateReconOption SofterFraction newdef KinematicsReconstructor:InitialInitialBoostOption LongTransBoost newdef /Herwig/Partons/RemnantDecayer:AlphaS AlphaQCD newdef /Herwig/Partons/RemnantDecayer:AlphaEM AlphaQED newdef ShowerHandler:PartnerFinder PartnerFinder newdef ShowerHandler:KinematicsReconstructor KinematicsReconstructor newdef ShowerHandler:SplittingGenerator SplittingGenerator newdef ShowerHandler:SpinCorrelations Yes newdef ShowerHandler:SoftCorrelations Singular ################################################################## # Intrinsic pT # # Recommended: # 1.9 GeV for Tevatron W/Z production. # 2.1 GeV for LHC W/Z production at 10 TeV # 2.2 GeV for LHC W/Z production at 14 TeV # # Set all parameters to 0 to disable ################################################################## newdef ShowerHandler:IntrinsicPtGaussian 1.3*GeV newdef ShowerHandler:IntrinsicPtBeta 0 newdef ShowerHandler:IntrinsicPtGamma 0*GeV newdef ShowerHandler:IntrinsicPtIptmax 0*GeV ############################################################# # Set up truncated shower handler. ############################################################# create Herwig::PowhegShowerHandler PowhegShowerHandler set PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler set PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler newdef PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:PartnerFinder PartnerFinder newdef PowhegShowerHandler:KinematicsReconstructor KinematicsReconstructor newdef PowhegShowerHandler:SplittingGenerator SplittingGenerator newdef PowhegShowerHandler:Interactions QCDandQED newdef PowhegShowerHandler:SpinCorrelations Yes newdef PowhegShowerHandler:SoftCorrelations Singular newdef PowhegShowerHandler:IntrinsicPtGaussian 1.3*GeV newdef PowhegShowerHandler:IntrinsicPtBeta 0 newdef PowhegShowerHandler:IntrinsicPtGamma 0*GeV newdef PowhegShowerHandler:IntrinsicPtIptmax 0*GeV newdef PowhegShowerHandler:ReconstructionOption OffShell5 ############################################################# # End of interesting user servicable section. # # Anything that follows below should only be touched if you # know what you're doing. # # Really. ############################################################# # # a few default values newdef ShowerHandler:MECorrMode 1 newdef ShowerHandler:ReconstructionOption OffShell5 newdef AlphaQCD:ScaleFactor 1.0 newdef AlphaQCD:NPAlphaS 2 newdef AlphaQCD:Qmin 0.935 newdef AlphaQCD:NumberOfLoops 2 newdef AlphaQCD:InputOption 1 newdef AlphaQCD:AlphaMZ 0.126234 # # # Lets set up all the splittings create Herwig::HalfHalfOneSplitFn QtoQGammaSplitFn set QtoQGammaSplitFn:InteractionType QED set QtoQGammaSplitFn:ColourStructure ChargedChargedNeutral set QtoQGammaSplitFn:AngularOrdered Yes create Herwig::HalfHalfOneSplitFn QtoQGSplitFn newdef QtoQGSplitFn:InteractionType QCD newdef QtoQGSplitFn:ColourStructure TripletTripletOctet set QtoQGSplitFn:AngularOrdered Yes create Herwig::OneOneOneSplitFn GtoGGSplitFn newdef GtoGGSplitFn:InteractionType QCD newdef GtoGGSplitFn:ColourStructure OctetOctetOctet set GtoGGSplitFn:AngularOrdered Yes create Herwig::OneOneOneMassiveSplitFn WtoWGammaSplitFn newdef WtoWGammaSplitFn:InteractionType QED newdef WtoWGammaSplitFn:ColourStructure ChargedChargedNeutral set WtoWGammaSplitFn:AngularOrdered Yes create Herwig::OneHalfHalfSplitFn GtoQQbarSplitFn newdef GtoQQbarSplitFn:InteractionType QCD newdef GtoQQbarSplitFn:ColourStructure OctetTripletTriplet set GtoQQbarSplitFn:AngularOrdered Yes create Herwig::OneHalfHalfSplitFn GammatoQQbarSplitFn newdef GammatoQQbarSplitFn:InteractionType QED newdef GammatoQQbarSplitFn:ColourStructure NeutralChargedCharged set GammatoQQbarSplitFn:AngularOrdered Yes create Herwig::HalfOneHalfSplitFn QtoGQSplitFn newdef QtoGQSplitFn:InteractionType QCD newdef QtoGQSplitFn:ColourStructure TripletOctetTriplet set QtoGQSplitFn:AngularOrdered Yes create Herwig::HalfOneHalfSplitFn QtoGammaQSplitFn newdef QtoGammaQSplitFn:InteractionType QED newdef QtoGammaQSplitFn:ColourStructure ChargedNeutralCharged set QtoGammaQSplitFn:AngularOrdered Yes # # Now the Sudakovs -create Herwig::QTildeSudakov SudakovCommon +create Herwig::SudakovFormFactor SudakovCommon newdef SudakovCommon:Alpha AlphaQCD newdef SudakovCommon:cutoffKinScale 0.0*GeV newdef SudakovCommon:PDFmax 1.0 newdef SudakovCommon:CutOffOption pT newdef SudakovCommon:pTmin 1.222798*GeV cp SudakovCommon QtoQGSudakov newdef QtoQGSudakov:SplittingFunction QtoQGSplitFn newdef QtoQGSudakov:PDFmax 1.9 cp SudakovCommon QtoQGammaSudakov set QtoQGammaSudakov:SplittingFunction QtoQGammaSplitFn set QtoQGammaSudakov:Alpha AlphaQED set QtoQGammaSudakov:PDFmax 1.9 cp QtoQGammaSudakov LtoLGammaSudakov # Technical parameter to stop evolution. set LtoLGammaSudakov:pTmin 0.000001 cp SudakovCommon GtoGGSudakov newdef GtoGGSudakov:SplittingFunction GtoGGSplitFn newdef GtoGGSudakov:PDFmax 2.0 cp SudakovCommon WtoWGammaSudakov newdef WtoWGammaSudakov:SplittingFunction WtoWGammaSplitFn set WtoWGammaSudakov:Alpha AlphaQED cp SudakovCommon GtoQQbarSudakov newdef GtoQQbarSudakov:SplittingFunction GtoQQbarSplitFn newdef GtoQQbarSudakov:PDFmax 120.0 cp SudakovCommon GammatoQQbarSudakov newdef GammatoQQbarSudakov:SplittingFunction GammatoQQbarSplitFn set GammatoQQbarSudakov:Alpha AlphaQED newdef GammatoQQbarSudakov:PDFmax 120.0 cp SudakovCommon GtobbbarSudakov newdef GtobbbarSudakov:SplittingFunction GtoQQbarSplitFn newdef GtobbbarSudakov:PDFmax 40000.0 cp SudakovCommon GtoccbarSudakov newdef GtoccbarSudakov:SplittingFunction GtoQQbarSplitFn newdef GtoccbarSudakov:PDFmax 2000.0 cp SudakovCommon QtoGQSudakov newdef QtoGQSudakov:SplittingFunction QtoGQSplitFn cp SudakovCommon QtoGammaQSudakov newdef QtoGammaQSudakov:SplittingFunction QtoGammaQSplitFn set QtoGammaQSudakov:Alpha AlphaQED cp SudakovCommon utoGuSudakov newdef utoGuSudakov:SplittingFunction QtoGQSplitFn newdef utoGuSudakov:PDFFactor OverOneMinusZ newdef utoGuSudakov:PDFmax 5.0 cp SudakovCommon dtoGdSudakov newdef dtoGdSudakov:SplittingFunction QtoGQSplitFn newdef dtoGdSudakov:PDFFactor OverOneMinusZ # # Now add the final splittings # do SplittingGenerator:AddFinalSplitting u->u,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting d->d,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting s->s,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting c->c,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting b->b,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting t->t,g; QtoQGSudakov # do SplittingGenerator:AddFinalSplitting g->g,g; GtoGGSudakov # do SplittingGenerator:AddFinalSplitting g->u,ubar; GtoQQbarSudakov do SplittingGenerator:AddFinalSplitting g->d,dbar; GtoQQbarSudakov do SplittingGenerator:AddFinalSplitting g->s,sbar; GtoQQbarSudakov do SplittingGenerator:AddFinalSplitting g->c,cbar; GtoccbarSudakov do SplittingGenerator:AddFinalSplitting g->b,bbar; GtobbbarSudakov do SplittingGenerator:AddFinalSplitting g->t,tbar; GtoQQbarSudakov # do SplittingGenerator:AddFinalSplitting gamma->u,ubar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->d,dbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->s,sbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->c,cbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->b,bbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->t,tbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->e-,e+; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->mu-,mu+; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->tau-,tau+; GammatoQQbarSudakov # do SplittingGenerator:AddFinalSplitting u->u,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting d->d,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting s->s,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting c->c,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting b->b,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting t->t,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting e-->e-,gamma; LtoLGammaSudakov do SplittingGenerator:AddFinalSplitting mu-->mu-,gamma; LtoLGammaSudakov do SplittingGenerator:AddFinalSplitting tau-->tau-,gamma; LtoLGammaSudakov do SplittingGenerator:AddFinalSplitting W+->W+,gamma; WtoWGammaSudakov # # Now lets add the initial splittings. Remember the form a->b,c; means # that the current particle b is given and we backward branch to new # particle a which is initial state and new particle c which is final state # do SplittingGenerator:AddInitialSplitting u->u,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting d->d,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting s->s,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting c->c,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting b->b,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting u->u,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting d->d,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting s->s,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting c->c,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting b->b,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting t->t,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting g->g,g; GtoGGSudakov # do SplittingGenerator:AddInitialSplitting g->d,dbar; GtoQQbarSudakov do SplittingGenerator:AddInitialSplitting g->u,ubar; GtoQQbarSudakov do SplittingGenerator:AddInitialSplitting g->s,sbar; GtoQQbarSudakov do SplittingGenerator:AddInitialSplitting g->c,cbar; GtoccbarSudakov do SplittingGenerator:AddInitialSplitting g->b,bbar; GtobbbarSudakov # do SplittingGenerator:AddInitialSplitting gamma->d,dbar; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->u,ubar; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->s,sbar; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->c,cbar; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->b,bbar; GammatoQQbarSudakov # do SplittingGenerator:AddInitialSplitting d->g,d; dtoGdSudakov do SplittingGenerator:AddInitialSplitting u->g,u; utoGuSudakov do SplittingGenerator:AddInitialSplitting s->g,s; QtoGQSudakov do SplittingGenerator:AddInitialSplitting c->g,c; QtoGQSudakov do SplittingGenerator:AddInitialSplitting b->g,b; QtoGQSudakov do SplittingGenerator:AddInitialSplitting dbar->g,dbar; dtoGdSudakov do SplittingGenerator:AddInitialSplitting ubar->g,ubar; utoGuSudakov do SplittingGenerator:AddInitialSplitting sbar->g,sbar; QtoGQSudakov do SplittingGenerator:AddInitialSplitting cbar->g,cbar; QtoGQSudakov do SplittingGenerator:AddInitialSplitting bbar->g,bbar; QtoGQSudakov # do SplittingGenerator:AddInitialSplitting d->gamma,d; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting u->gamma,u; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting s->gamma,s; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting c->gamma,c; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting b->gamma,b; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting dbar->gamma,dbar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting ubar->gamma,ubar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting sbar->gamma,sbar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting cbar->gamma,cbar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting bbar->gamma,bbar; QtoGammaQSudakov