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,525 +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()], - bornCXComb()->mePartonData()[dipole()->bornEmitter()]->iColour() == PDT::Colour3).first; + 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 ("QTildeSudakov", "Set the partner finder to calculate hard scales.", &QTildeMatching::theQTildeSudakov, false, false, true, false, false); interfaceQTildeSudakov.rank(-1); static Reference interfaceShowerHandler ("ShowerHandler", "The QTilde shower handler to use.", &QTildeMatching::theShowerHandler, false, false, true, true, false); interfaceShowerHandler.rank(-1); static Switch interfaceCorrectForXZMismatch ("CorrectForXZMismatch", "Correct for x/z mismatch near hard phase space boundary.", &QTildeMatching::theCorrectForXZMismatch, true, false, false); static SwitchOption interfaceCorrectForXZMismatchYes (interfaceCorrectForXZMismatch, "Yes", "Include the correction factor.", true); static SwitchOption interfaceCorrectForXZMismatchNo (interfaceCorrectForXZMismatch, "No", "Do not include the correction factor.", false); interfaceCorrectForXZMismatch.rank(-1); } diff --git a/Shower/QTilde/Base/PartnerFinder.cc b/Shower/QTilde/Base/PartnerFinder.cc --- a/Shower/QTilde/Base/PartnerFinder.cc +++ b/Shower/QTilde/Base/PartnerFinder.cc @@ -1,757 +1,628 @@ // -*- C++ -*- // // PartnerFinder.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 PartnerFinder class. // #include "PartnerFinder.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Shower/QTilde/Base/ShowerParticle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Utilities/Debug.h" #include "ThePEG/Utilities/DescribeClass.h" using namespace Herwig; DescribeClass describePartnerFinder ("Herwig::PartnerFinder","HwShower.so"); // some useful functions to avoid using #define namespace { // return bool if final-state particle inline bool FS(const tShowerParticlePtr a) { return a->isFinalState(); } // return colour line pointer inline Ptr::transient_pointer CL(const tShowerParticlePtr a, unsigned int index=0) { return a->colourInfo()->colourLines().empty() ? ThePEG::tColinePtr() : const_ptr_cast(a->colourInfo()->colourLines()[index]); } // return colour line size inline unsigned int CLSIZE(const tShowerParticlePtr a) { return a->colourInfo()->colourLines().size(); } inline Ptr::transient_pointer ACL(const tShowerParticlePtr a, unsigned int index=0) { return a->colourInfo()->antiColourLines().empty() ? ThePEG::tColinePtr() : const_ptr_cast(a->colourInfo()->antiColourLines()[index]); } inline unsigned int ACLSIZE(const tShowerParticlePtr a) { return a->colourInfo()->antiColourLines().size(); } } void PartnerFinder::persistentOutput(PersistentOStream & os) const { - os << partnerMethod_ << QEDPartner_ << scaleChoice_ - << _finalFinalConditions << _initialFinalDecayConditions - << _initialInitialConditions; + os << partnerMethod_ << QEDPartner_ << scaleChoice_; } void PartnerFinder::persistentInput(PersistentIStream & is, int) { - is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_ - >> _finalFinalConditions >> _initialFinalDecayConditions - >>_initialInitialConditions; + is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_; } void PartnerFinder::Init() { static ClassDocumentation documentation ("This class is responsible for finding the partners for each interaction types ", "and within the evolution scale range specified by the ShowerVariables ", "then to determine the initial evolution scales for each pair of partners."); static Switch interfacePartnerMethod ("PartnerMethod", "Choice of partner finding method for gluon evolution.", &PartnerFinder::partnerMethod_, 0, false, false); static SwitchOption interfacePartnerMethodRandom (interfacePartnerMethod, "Random", "Choose partners of a gluon randomly.", 0); static SwitchOption interfacePartnerMethodMaximum (interfacePartnerMethod, "Maximum", "Choose partner of gluon with largest angle.", 1); static Switch interfaceQEDPartner ("QEDPartner", "Control of which particles to use as the partner for QED radiation", &PartnerFinder::QEDPartner_, 0, false, false); static SwitchOption interfaceQEDPartnerAll (interfaceQEDPartner, "All", "Consider all possible choices which give a positive contribution" " in the soft limit.", 0); static SwitchOption interfaceQEDPartnerIIandFF (interfaceQEDPartner, "IIandFF", "Only allow initial-initial or final-final combinations", 1); static SwitchOption interfaceQEDPartnerIF (interfaceQEDPartner, "IF", "Only allow initial-final combinations", 2); static Switch interfaceScaleChoice ("ScaleChoice", "The choice of the evolution scales", &PartnerFinder::scaleChoice_, 0, false, false); static SwitchOption interfaceScaleChoicePartner (interfaceScaleChoice, "Partner", "Scale of all interactions is that of the evolution partner", 0); static SwitchOption interfaceScaleChoiceDifferent (interfaceScaleChoice, "Different", "Allow each interaction to have different scales", 1); - - - static Switch interfaceFinalFinalConditions - ("FinalFinalConditions", - "The initial conditions for the shower of a final-final colour connection", - &PartnerFinder::_finalFinalConditions, 0, false, false); - static SwitchOption interfaceFinalFinalConditionsSymmetric - (interfaceFinalFinalConditions, - "Symmetric", - "The symmetric choice", - 0); - static SwitchOption interfaceFinalFinalConditionsColoured - (interfaceFinalFinalConditions, - "Coloured", - "Maximal radiation from the coloured particle", - 1); - static SwitchOption interfaceFinalFinalConditionsAntiColoured - (interfaceFinalFinalConditions, - "AntiColoured", - "Maximal emission from the anticoloured particle", - 2); - static SwitchOption interfaceFinalFinalConditionsRandom - (interfaceFinalFinalConditions, - "Random", - "Randomly selected maximal emission from one of the particles", - 3); - - static Switch interfaceInitialFinalDecayConditions - ("InitialFinalDecayConditions", - "The initial conditions for the shower of an initial-final" - " decay colour connection.", - &PartnerFinder::_initialFinalDecayConditions, 0, false, false); - static SwitchOption interfaceInitialFinalDecayConditionsSymmetric - (interfaceInitialFinalDecayConditions, - "Symmetric", - "The symmetric choice", - 0); - static SwitchOption interfaceInitialFinalDecayConditionsMaximal - (interfaceInitialFinalDecayConditions, - "Maximal", - "Maximal radiation from the decay product", - 1); - static SwitchOption interfaceInitialFinalDecayConditionsSmooth - (interfaceInitialFinalDecayConditions, - "Smooth", - "Smooth matching in the soft limit", - 2); - - static Switch interfaceInitialInitialConditions - ("InitialInitialConditions", - "The initial conditions for the shower of an initial-initial" - " colour connection.", - &PartnerFinder::_initialInitialConditions, 0, false, false); - static SwitchOption interfaceInitialInitialConditionsSymmetric - (interfaceInitialInitialConditions, - "Symmetric", - "The symmetric choice", - 0); - static SwitchOption interfaceInitialInitialConditionsMaximiseB - (interfaceInitialInitialConditions, - "MaximiseB", - "Maximal radiation from parton b", - 1); - static SwitchOption interfaceInitialInitialConditionsMaximiseC - (interfaceInitialInitialConditions, - "MaximiseC", - "Maximal radiation from parton c", - 2); } void PartnerFinder::setInitialEvolutionScales(const ShowerParticleVector &particles, const bool isDecayCase, ShowerInteraction type, const bool setPartners) { // clear the existing partners for(ShowerParticleVector::const_iterator cit = particles.begin(); cit != particles.end(); ++cit) (*cit)->clearPartners(); // set them if(type==ShowerInteraction::QCD) { setInitialQCDEvolutionScales(particles,isDecayCase,setPartners); } else if(type==ShowerInteraction::QED) { setInitialQEDEvolutionScales(particles,isDecayCase,setPartners); } else { setInitialQCDEvolutionScales(particles,isDecayCase,setPartners); setInitialQEDEvolutionScales(particles,isDecayCase,false); } // print out for debugging if(Debug::level>=10) { for(ShowerParticleVector::const_iterator cit = particles.begin(); cit != particles.end(); ++cit) { generator()->log() << "Particle: " << **cit << "\n"; if(!(**cit).partner()) continue; generator()->log() << "Primary partner: " << *(**cit).partner() << "\n"; for(vector::const_iterator it= (**cit).partners().begin(); it!=(**cit).partners().end();++it) { generator()->log() << static_cast(it->type) << " " << it->weight << " " << it->scale/GeV << " " << *(it->partner) << "\n"; } } generator()->log() << flush; } } void PartnerFinder::setInitialQCDEvolutionScales(const ShowerParticleVector &particles, const bool isDecayCase, const bool setPartners) { // set the partners and the scales if(setPartners) { // Loop over particles and consider only coloured particles which don't // have already their colour partner fixed and that don't have children // (the latter requirement is relaxed in the case isDecayCase is true). // Build a map which has as key one of these particles (i.e. a pointer // to a ShowerParticle object) and as a corresponding value the vector // of all its possible *normal* candidate colour partners, defined as follows: // --- have colour, and no children (this is not required in the case // isDecayCase is true); // --- if both are initial (incoming) state particles, then the (non-null) colourLine() // of one of them must match the (non-null) antiColourLine() of the other. // --- if one is an initial (incoming) state particle and the other is // a final (outgoing) state particle, then both must have the // same (non-null) colourLine() or the same (non-null) antiColourLine(); // Notice that this definition exclude the special case of baryon-violating // processes (as in R-parity Susy), which will show up as particles // without candidate colour partners, and that we will be treated a part later // (this means that no modifications of the following loop is needed!) ShowerParticleVector::const_iterator cit, cjt; for(cit = particles.begin(); cit != particles.end(); ++cit) { // Skip colourless particles if(!(*cit)->data().coloured()) continue; // find the partners vector< pair > partners = findQCDPartners(*cit,particles); // must have a partner if(partners.empty()) { throw Exception() << "`Failed to make colour connections in " << "PartnerFinder::setQCDInitialEvolutionScales" << (**cit) << Exception::eventerror; } // Calculate the evolution scales for all possible pairs of of particles vector > scales; for(unsigned int ix=0;ix< partners.size();++ix) { scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second), isDecayCase)); } // In the case of more than one candidate colour partners, // there are now two approaches to choosing the partner. The // first method is based on two assumptions: // 1) the choice of which is the colour partner is done // *randomly* between the available candidates. // 2) the choice of which is the colour partner is done // *independently* from each particle: in other words, // if for a particle "i" its selected colour partner is // the particle "j", then the colour partner of "j" // does not have to be necessarily "i". // The second method always chooses the furthest partner // for hard gluons and gluinos. // store the choice int position(-1); // random choice if( partnerMethod_ == 0 ) { // random choice of partner position = UseRandom::irnd(partners.size()); } // take the one with largest angle else if (partnerMethod_ == 1 ) { if ((*cit)->perturbative() == 1 && (*cit)->dataPtr()->iColour()==PDT::Colour8 ) { assert(partners.size()==2); // Determine largest angle double maxAngle(0.); for(unsigned int ix=0;ixmomentum().vect(). angle(partners[ix].second->momentum().vect()); if(angle>maxAngle) { maxAngle = angle; position = ix; } } } else position = UseRandom::irnd(partners.size()); } else assert(false); // set the evolution partner (*cit)->partner(partners[position].second); for(unsigned int ix=0;ixdata().coloured()) continue; // find the partners vector< pair > partners = findQCDPartners(*cit,particles); // must have a partner if(partners.empty()) { throw Exception() << "`Failed to make colour connections in " << "PartnerFinder::setQCDInitialEvolutionScales" << (**cit) << Exception::eventerror; } // Calculate the evolution scales for all possible pairs of of particles vector > scales; int position(-1); for(unsigned int ix=0;ix< partners.size();++ix) { if(partners[ix].second) position = ix; scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second), isDecayCase)); } assert(position>=0); for(unsigned int ix=0;ixcharged()) continue; // find the potential partners vector > partners = findQEDPartners(*cit,particles,isDecayCase); if(partners.empty()) { throw Exception() << "Failed to find partner in " << "PartnerFinder::setQEDInitialEvolutionScales" << (**cit) << Exception::eventerror; } // calculate the probabilities double prob(0.); for(unsigned int ix=0;ixpartner()) { for(unsigned int ix=0;ixpartner()==partners[ix].second) { position = ix; break; } } } // set the partner if(setPartners||!(*cit)->partner()||position<0) { prob = UseRandom::rnd(); for(unsigned int ix=0;ixprob) { position = ix; break; } prob -= partners[ix].first; } if(position>=0&&(setPartners||!(*cit)->partner())) { (*cit)->partner(partners[position].second); } } // must have a partner if(position<0) throw Exception() << "Failed to find partner in " << "PartnerFinder::setQEDInitialEvolutionScales" << (**cit) << Exception::eventerror; // Calculate the evolution scales for all possible pairs of of particles vector > scales; for(unsigned int ix=0;ix< partners.size();++ix) { scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second), isDecayCase)); } // store all the possible partners for(unsigned int ix=0;ix PartnerFinder:: calculateInitialEvolutionScales(const ShowerPPair &particlePair, const bool isDecayCase) { bool FS1=FS(particlePair.first),FS2= FS(particlePair.second); - if(FS1 && FS2) { - bool colouredFirst = - particlePair.first->colourLine()&& - particlePair.first->colourLine()==particlePair.second->antiColourLine(); - return calculateFinalFinalScales(particlePair.first->momentum(), - particlePair.second->momentum(), - colouredFirst); - } + if(FS1 && FS2) + return calculateFinalFinalScales(particlePair.first->momentum(), + particlePair.second->momentum()); else if(FS1 && !FS2) { - ShowerPPair a(particlePair.second, particlePair.first); - pair rval = calculateInitialFinalScales(a.first->momentum(), - a.second->momentum(), + pair rval = calculateInitialFinalScales(particlePair.second->momentum(), + particlePair.first->momentum(), isDecayCase); - return pair(rval.second,rval.first); + return { rval.second, rval.first }; } else if(!FS1 &&FS2) return calculateInitialFinalScales(particlePair.first->momentum(),particlePair.second->momentum(),isDecayCase); else return calculateInitialInitialScales(particlePair.first->momentum(),particlePair.second->momentum()); } vector< pair > PartnerFinder::findQCDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles) { vector< pair > partners; ShowerParticleVector::const_iterator cjt; for(cjt = particles.begin(); cjt != particles.end(); ++cjt) { if(!(*cjt)->data().coloured() || particle==*cjt) continue; // one initial-state and one final-state particle if(FS(particle) != FS(*cjt)) { // loop over all the colours of both particles for(unsigned int ix=0; ixsourceNeighbours().first) { tColinePair cpair = col->sourceNeighbours(); for(cjt=particles.begin();cjt!=particles.end();++cjt) { if(( FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second))|| (!FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt) == cpair.second ))) { partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt)); } } } else if(col&&col->sinkNeighbours().first) { tColinePair cpair = col->sinkNeighbours(); for(cjt=particles.begin();cjt!=particles.end();++cjt) { if(( FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt) == cpair.second))|| (!FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second))) { partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt)); } } } col = ACL(particle); if(FS(particle)&&col&&col->sinkNeighbours().first) { tColinePair cpair = col->sinkNeighbours(); for(cjt=particles.begin();cjt!=particles.end();++cjt) { if(( FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt) == cpair.second))|| (!FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second ))) { partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt)); } } } else if(col&&col->sourceNeighbours().first) { tColinePair cpair = col->sourceNeighbours(); for(cjt=particles.begin();cjt!=particles.end();++cjt) { if(( FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second))|| (!FS(*cjt) && (ACL(*cjt) == cpair.first ||ACL(*cjt) == cpair.second))) { partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt)); } } } } // return the partners return partners; } vector< pair > PartnerFinder::findQEDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles, const bool isDecayCase) { vector< pair > partners; ShowerParticleVector::const_iterator cjt; double pcharge = particle->id()==ParticleID::gamma ? 1 : double(particle->data().iCharge()); vector< pair > photons; for(cjt = particles.begin(); cjt != particles.end(); ++cjt) { if(particle == *cjt) continue; if((**cjt).id()==ParticleID::gamma) photons.push_back(make_pair(1.,*cjt)); if(!(*cjt)->data().charged() ) continue; double charge = pcharge*double((*cjt)->data().iCharge()); if( FS(particle) != FS(*cjt) ) charge *=-1.; if( QEDPartner_ != 0 && !isDecayCase ) { // only include II and FF as requested if( QEDPartner_ == 1 && FS(particle) != FS(*cjt) ) continue; // only include IF is requested else if(QEDPartner_ == 2 && FS(particle) == FS(*cjt) ) continue; } if(particle->id()==ParticleID::gamma) charge = -abs(charge); // only keep positive dipoles if(charge<0.) partners.push_back(make_pair(-charge,*cjt)); } if(particle->id()==ParticleID::gamma&& partners.empty()) { return photons; } return partners; } pair PartnerFinder::calculateFinalFinalScales( const Lorentz5Momentum & p1, - const Lorentz5Momentum & p2, - bool colouredFirst) + const Lorentz5Momentum & p2) { static const double eps=1e-7; // Using JHEP 12(2003)045 we find that we need ktilde = 1/2(1+b-c+lambda) // ktilde = qtilde^2/Q^2 therefore qtilde = sqrt(ktilde*Q^2) // find momenta in rest frame of system // calculate quantities for the scales Energy2 Q2 = (p1+p2).m2(); double b = p1.mass2()/Q2; double c = p2.mass2()/Q2; if(b<0.) { if(b<-eps) { throw Exception() << "Negative Mass squared b = " << b << "in PartnerFinder::calculateFinalFinalScales()" << Exception::eventerror; } b = 0.; } if(c<0.) { if(c<-eps) { throw Exception() << "Negative Mass squared c = " << c << "in PartnerFinder::calculateFinalFinalScales()" << Exception::eventerror; } c = 0.; } // KMH & PR - 16 May 2008 - swapped lambda calculation from // double lam=2.*p1.vect().mag()/Q; to sqrt(kallen(1,b,c)), // which should be identical for p1 & p2 onshell in their COM // but in the inverse construction for the Nason method, this // was not the case, leading to misuse. - double lam=sqrt((1.+sqrt(b)+sqrt(c))*(1.-sqrt(b)-sqrt(c)) - *(sqrt(b)-1.-sqrt(c))*(sqrt(c)-1.-sqrt(b))); + const double lam=sqrt((1.+sqrt(b)+sqrt(c))*(1.-sqrt(b)-sqrt(c)) + *(sqrt(b)-1.-sqrt(c))*(sqrt(c)-1.-sqrt(b))); // symmetric case - unsigned int iopt=finalFinalConditions(); - Energy firstQ,secondQ; - if(iopt==0) { - firstQ = sqrt(0.5*Q2*(1.+b-c+lam)); - secondQ = sqrt(0.5*Q2*(1.-b+c+lam)); - } - // assymetric choice - else { - double kappab,kappac; - // calculate kappa with coloured line getting maximum - if((iopt==1&&colouredFirst)|| // first particle coloured+maximal for coloured - (iopt==2&&!colouredFirst)|| // first particle anticoloured+maximal for acoloured - (iopt==3&&UseRandom::rndbool(0.5))) { // random choice - kappab=4.*(1.-2.*sqrt(c)-b+c); - kappac=c+0.25*sqr(1.-b-c+lam)/(kappab-b); - } - else { - kappac=4.*(1.-2.*sqrt(b)-c+b); - kappab=b+0.25*sqr(1.-b-c+lam)/(kappac-c); - } - // calculate the scales - firstQ = sqrt(Q2*kappab); - secondQ = sqrt(Q2*kappac); - } - return pair(firstQ, secondQ); + return { sqrt(0.5*Q2*(1.+b-c+lam)), sqrt(0.5*Q2*(1.-b+c+lam)) }; } pair PartnerFinder::calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc, const bool isDecayCase) { if(!isDecayCase) { // In this case from JHEP 12(2003)045 we find the conditions // ktilde_b = (1+c) and ktilde_c = (1+2c) // We also find that c = m_c^2/Q^2. The process is a+b->c where // particle a is not colour connected (considered as a colour singlet). // Therefore we simply find that q_b = sqrt(Q^2+m_c^2) and // q_c = sqrt(Q^2+2 m_c^2) // We also assume that the first particle in the pair is the initial // state particle and the second is the final state one c - Energy2 mc2 = sqr(pc.mass()); - Energy2 Q2 = -(pb-pc).m2(); - return pair(sqrt(Q2+mc2), sqrt(Q2+2*mc2)); + const Energy2 mc2 = sqr(pc.mass()); + const Energy2 Q2 = -(pb-pc).m2(); + return { sqrt(Q2+mc2), sqrt(Q2+2*mc2) }; } else { // In this case from JHEP 12(2003)045 we find, for the decay // process b->c+a(neutral), the condition // (ktilde_b-1)*(ktilde_c-c)=(1/4)*sqr(1-a+c+lambda). // We also assume that the first particle in the pair is the initial // state particle (b) and the second is the final state one (c). // - We find maximal phase space coverage through emissions from // c if we set ktilde_c = 4.*(sqr(1.-sqrt(a))-c) // - We find the most 'symmetric' way to populate the phase space // occurs for (ktilde_b-1)=(ktilde_c-c)=(1/2)*(1-a+c+lambda) // - We find the most 'smooth' way to populate the phase space // occurs for... Energy2 mb2(sqr(pb.mass())); double a=(pb-pc).m2()/mb2; double c=sqr(pc.mass())/mb2; double lambda = 1. + a*a + c*c - 2.*a - 2.*c - 2.*a*c; lambda = sqrt(max(lambda,0.)); - double PROD = 0.25*sqr(1. - a + c + lambda); - double ktilde_b, ktilde_c,cosi(0.); - switch(initialFinalDecayConditions()) { - case 0: // the 'symmetric' choice - ktilde_c = 0.5*(1-a+c+lambda) + c ; - ktilde_b = 1.+PROD/(ktilde_c-c) ; - break; - case 1: // the 'maximal' choice - ktilde_c = 4.0*(sqr(1.-sqrt(a))-c); - ktilde_b = 1.+PROD/(ktilde_c-c) ; - break; - case 2: // the 'smooth' choice - // c is a problem if very small here use 1GeV as minimum - c = max(c,1.*GeV2/mb2); - cosi = (sqr(1-sqrt(c))-a)/lambda; - ktilde_b = 2.0/(1.0-cosi); - ktilde_c = (1.0-a+c+lambda)*(1.0+c-a-lambda*cosi)/(2.0*(1.0+cosi)); - break; - default: - throw Exception() << "Invalid option for decay shower's phase space" - << " PartnerFinder::calculateInitialFinalScales" - << Exception::abortnow; - } - return pair(sqrt(mb2*ktilde_b),sqrt(mb2*ktilde_c)); + const double PROD = 0.25*sqr(1. - a + c + lambda); + const double ktilde_c = 0.5*(1-a+c+lambda) + c ; + const double ktilde_b = 1.+PROD/(ktilde_c-c) ; + return { sqrt(mb2*ktilde_b), sqrt(mb2*ktilde_c) }; } } pair PartnerFinder::calculateInitialInitialScales(const Lorentz5Momentum& p1, const Lorentz5Momentum& p2) { // This case is quite simple. From JHEP 12(2003)045 we find the condition // that ktilde_b = ktilde_c = 1. In this case we have the process // b+c->a so we need merely boost to the CM frame of the two incoming // particles and then qtilde is equal to the energy in that frame - Energy Q = sqrt((p1+p2).m2()); - if(_initialInitialConditions==1) { - return pair(sqrt(2.0)*Q,sqrt(0.5)*Q); - } else if(_initialInitialConditions==2) { - return pair(sqrt(0.5)*Q,sqrt(2.0)*Q); - } else { - return pair(Q,Q); - } + const Energy Q = sqrt((p1+p2).m2()); + return {Q,Q}; } diff --git a/Shower/QTilde/Base/PartnerFinder.h b/Shower/QTilde/Base/PartnerFinder.h --- a/Shower/QTilde/Base/PartnerFinder.h +++ b/Shower/QTilde/Base/PartnerFinder.h @@ -1,297 +1,230 @@ // -*- C++ -*- // // PartnerFinder.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_PartnerFinder_H #define HERWIG_PartnerFinder_H // // This is the declaration of the PartnerFinder class. // #include "Herwig/Shower/QTilde/ShowerConfig.h" #include "ThePEG/Interface/Interfaced.h" #include "PartnerFinder.fh" namespace Herwig { using namespace ThePEG; /** * typedef of a pair of particle for calculating the evolution scales */ typedef pair ShowerPPair; /** \ingroup Shower * * This class is responsible of two related tasks: * - it finds the partners * - for each pair of partners (and interaction therefore) * it sets the initial evolution scales of both of them. * * In general the finding of the partners is performed by this class but * the calculation of the initial evolution scales should be implemented * for different shower evolution models in classes inheriting from this one. * Notice that a given particle has, in general, a different partner * for each different interaction; however, given a partner, its * initial evolution scale, Q, is purely a kinematical relationship * between the pair, without dependence on the dynamics (i.e. type of interaction). * * @see \ref PartnerFinderInterfaces "The interfaces" * defined for PartnerFinder. */ class PartnerFinder: public Interfaced { public: /** * The default constructor. */ - PartnerFinder() : partnerMethod_(0), QEDPartner_(0), scaleChoice_(0), - _finalFinalConditions(0), - _initialFinalDecayConditions(0), - _initialInitialConditions(0) {} + PartnerFinder() : partnerMethod_(0), QEDPartner_(0), scaleChoice_(0) {} /** * Given in input a collection of particles (ShowerParticle objects), * each of these methods set the initial evolution scales of those particles, * between the ones given in input, that do not have yet their * evolution scale set. * The input collection of particles can be either the full collection of * showering particles (kept in the main class ShowerHandler, * in the case isDecayCase is false, or simply, in the case isDecayCase * is true, the decaying particle and its decay products. * The methods returns true, unless something wrong (inconsistencies, * or undefined values) happens. * * These methods are virtual but in most cases inheriting classes should not * need to overide them as they simply find the relevant partner and call * one of the calculateScale members to calculate the scale. */ //@{ /** * Set the initial scales * @param particles The particles to be considered * @param isDecayCase Whether or not this is a decay * @param setPartners Whether to set the colour partners or just the scales */ virtual void setInitialEvolutionScales(const ShowerParticleVector &particles, const bool isDecayCase, ShowerInteraction, const bool setPartners=true); //@} protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const {return new_ptr(*this);} /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const {return new_ptr(*this);} //@} 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: /** - * Access function for the initial conditions for the shower - */ - //@{ - /** - * Initial conditions for the shower of a final-final colour connection - * - 0 is the symmetric choice - * - 1 is maximal emmision from the coloured particle - * - 2 is maximal emmision from the anticoloured particle - * - 3 is randomly selected maximal emmision - */ - unsigned int finalFinalConditions() const - {return _finalFinalConditions;} - - /** - * Initial conditions for the shower of an initial-final decay colour connection - * - 0 is the symmetric choice - * - 1 is maximal emission from the decay product - * - 2 is the smooth choice - */ - unsigned int initialFinalDecayConditions() const - {return _initialFinalDecayConditions;} - //@} - -protected: - - /** * Members to set the scales for different interactions */ //@{ /** * Set initial scales for a QCD interaction */ virtual void setInitialQCDEvolutionScales(const ShowerParticleVector &particles, const bool isDecayCase, const bool setPartners=true); /** * Set initial scales for a QED interaction */ virtual void setInitialQEDEvolutionScales(const ShowerParticleVector &particles, const bool isDecayCase, const bool setPartners=true); //@} /** * Find the QCD partners * @param particle The particle to find the partners for * @param particles The full set of particles to search */ vector< pair > findQCDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles); /** * Find the QED partners * @param particle The particle to find the partners for * @param particles The full set of particles to search */ vector< pair > findQEDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles, const bool isDecayCase); public: /** * Given a pair of particles, supposedly partners w.r.t. an interaction, * this method returns their initial evolution scales as a pair. * If something wrong happens, it returns the null (ZERO,ZERO) pair. * This method is used by the above setXXXInitialEvolutionScales * methods. * These methods must be overiden in inheriting classes */ //@{ /** * General method to calculate the initial evolution scales */ pair calculateInitialEvolutionScales(const ShowerPPair &, const bool isDecayCase); /** * Calculate the initial evolution scales given momenta */ pair calculateFinalFinalScales(const Lorentz5Momentum & p1, - const Lorentz5Momentum & p2, - bool colouredfirst); + const Lorentz5Momentum & p2); /** * Calculate the initial evolution scales given momenta */ pair calculateInitialInitialScales(const Lorentz5Momentum& p1, const Lorentz5Momentum& p2); /** * Calculate the initial evolution scales given momenta */ pair calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc, const bool isDecayCase); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ PartnerFinder & operator=(const PartnerFinder &); private: /** * Method for choosing colour partner */ int partnerMethod_; /** * Choice for the QED radiation partner */ int QEDPartner_; /** * Choice of the scale */ int scaleChoice_; - /** - * Flags controlling the initial conditions for the shower - */ - //@{ - /** - * Initial conditions for the shower with a final-final colour - * connection - */ - unsigned int _finalFinalConditions; - - /** - * Initial conditions for the shower with an initial-final decay colour - * connection. This is done according to the top decay colour - * connection calculation in JHEP12(2003)_045. The options act as follows: - * 0: This is the default 'symmetric' choice which more or less divides - * the phase space evenly between the parent and its charged child. - * 1: This 'maximal' choice maximises the phase space available for - * gluons emitted from the charged child. - * 2: This (experimental) 'smooth' choice does not suffer from - * a discontinuity at the boundary between the region populated by - * emissions from the charged child and the region populated by emissions - * from the parent. This does, however, mean that the phase space - * available for emissions from the charged child is fairly minimal. - */ - unsigned int _initialFinalDecayConditions; - - /** - * Initial conditions for the shower with an initial-initial colour - * connection. This is done according to the colour connection - * calculation in JHEP12(2003)_045. The options act as follows: - * 0: This is the default 'symmetric' choice which more or less divides - * the phase space evenly between the two incoming partons. - * 1: This increases the phase space for emission from "parton b". - * 2: This increases the phase space for emission from "parton c". - */ - unsigned int _initialInitialConditions; - //@} }; } #endif /* HERWIG_PartnerFinder_H */