diff --git a/Decay/General/SSVDecayer.cc b/Decay/General/SSVDecayer.cc --- a/Decay/General/SSVDecayer.cc +++ b/Decay/General/SSVDecayer.cc @@ -1,372 +1,343 @@ // -*- C++ -*- // // SSVDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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 SSVDecayer class. // #include "SSVDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/DecayMode.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" using namespace Herwig; using namespace ThePEG::Helicity; IBPtr SSVDecayer::clone() const { return new_ptr(*this); } IBPtr SSVDecayer::fullclone() const { return new_ptr(*this); } void SSVDecayer::setDecayInfo(PDPtr incoming, PDPair outgoing, vector vertex, map & inV, const vector > & outV, map fourV) { decayInfo(incoming,outgoing); - for(auto vert : vertex) { + for(auto vert : vertex) vertex_ .push_back(dynamic_ptr_cast(vert)); - perturbativeVertex_.push_back(dynamic_ptr_cast (vert)); - } vector itemp={ShowerInteraction::QCD,ShowerInteraction::QED}; for(auto & inter : itemp) { incomingVertex_[inter] = dynamic_ptr_cast(inV.at(inter)); fourPointVertex_[inter] = dynamic_ptr_cast(fourV.at(inter)); outgoingVertexS_[inter] = AbstractVSSVertexPtr(); outgoingVertexV_[inter] = AbstractVVVVertexPtr(); if(outV[0].at(inter)) { if (outV[0].at(inter)->getName()==VertexType::VSS) outgoingVertexS_[inter] = dynamic_ptr_cast(outV[0].at(inter)); else outgoingVertexV_[inter] = dynamic_ptr_cast(outV[0].at(inter)); } if(outV[1].at(inter)) { if (outV[1].at(inter)->getName()==VertexType::VSS) outgoingVertexS_[inter] = dynamic_ptr_cast(outV[1].at(inter)); else outgoingVertexV_[inter] = dynamic_ptr_cast(outV[1].at(inter)); } } } void SSVDecayer::persistentOutput(PersistentOStream & os) const { - os << vertex_ << perturbativeVertex_ + os << vertex_ << incomingVertex_ << outgoingVertexS_ << outgoingVertexV_ << fourPointVertex_; } void SSVDecayer::persistentInput(PersistentIStream & is, int) { - is >> vertex_ >> perturbativeVertex_ + is >> vertex_ >> incomingVertex_ >> outgoingVertexS_ >> outgoingVertexV_ >> fourPointVertex_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigSSVDecayer("Herwig::SSVDecayer", "Herwig.so"); void SSVDecayer::Init() { static ClassDocumentation documentation ("This implements the decay of a scalar to a vector and a scalar"); } void SSVDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { unsigned int isc(0),ivec(1); if(decay[0]->dataPtr()->iSpin() != PDT::Spin0) swap(isc,ivec); ScalarWaveFunction:: constructSpinInfo(const_ptr_cast(&part),incoming,true); ScalarWaveFunction:: constructSpinInfo(decay[isc],outgoing,true); VectorWaveFunction:: constructSpinInfo(vector_,decay[ivec],outgoing,true,false); } double SSVDecayer::me2(const int,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { unsigned int isc(0),ivec(1); if(outgoing[0]->iSpin() != PDT::Spin0) swap(isc,ivec); if(!ME()) { if(ivec==1) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin1))); else ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin1,PDT::Spin0))); } if(meopt==Initialize) { ScalarWaveFunction:: calculateWaveFunctions(rho_,const_ptr_cast(&part),incoming); swave_ = ScalarWaveFunction(part.momentum(),part.dataPtr(),incoming); // fix rho if no correlations fixRho(rho_); } VectorWaveFunction:: calculateWaveFunctions(vector_,momenta[ivec],outgoing[ivec],Helicity::outgoing,false); ScalarWaveFunction sca(momenta[isc],outgoing[isc],Helicity::outgoing); Energy2 scale(sqr(part.mass())); //make sure decay matrix element is in the correct order double output(0.); if(ivec == 0) { for(unsigned int ix = 0; ix < 3; ++ix) { (*ME())(0, ix, 0) = 0.; for(auto vert : vertex_) (*ME())(0, ix, 0) += vert->evaluate(scale,vector_[ix],sca, swave_); } } else { for(unsigned int ix = 0; ix < 3; ++ix) { (*ME())(0, 0, ix) = 0.; for(auto vert : vertex_) (*ME())(0, 0, ix) += vert->evaluate(scale,vector_[ix],sca,swave_); } } output = (ME()->contract(rho_)).real()/scale*UnitRemoval::E2; // colour and identical particle factors output *= colourFactor(part.dataPtr(),outgoing[0],outgoing[1]); // return the answer return output; } Energy SSVDecayer:: partialWidth(PMPair inpart, PMPair outa, PMPair outb) const { if( inpart.second < outa.second + outb.second ) return ZERO; - if(perturbativeVertex_.size()==1 && - perturbativeVertex_[0]) { - double mu1sq(sqr(outa.second/inpart.second)), - mu2sq(sqr(outb.second/inpart.second)); - tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first; - if(outa.first->iSpin() == PDT::Spin0) { - perturbativeVertex_[0]->setCoupling(sqr(inpart.second), outb.first, outa.first,in); - } - else { - swap(mu1sq,mu2sq); - perturbativeVertex_[0]->setCoupling(sqr(inpart.second), outa.first, outb.first,in); - } - double me2(0.); - if(mu2sq == 0.) - me2 = -2.*mu1sq - 2.; - else - me2 = ( sqr(mu2sq - mu1sq) - 2.*(mu2sq + mu1sq) + 1. )/mu2sq; - Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second, outa.second, - outb.second); - Energy output = pcm*me2*norm(perturbativeVertex_[0]->norm())/8./Constants::pi; - // colour factor - output *= colourFactor(inpart.first,outa.first,outb.first); - // return the answer - return output; - } - else { - return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb); - } + return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb); } double SSVDecayer::threeBodyME(const int , const Particle & inpart, const ParticleVector & decay, ShowerInteraction inter, MEOption meopt) { int iscal (0), ivect (1), iglu (2); // get location of outgoing scalar/vector if(decay[1]->dataPtr()->iSpin()==PDT::Spin0) swap(iscal,ivect); if(meopt==Initialize) { // create scalar wavefunction for decaying particle ScalarWaveFunction::calculateWaveFunctions(rho3_,const_ptr_cast(&inpart),incoming); swave3_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(),incoming); } // setup spin information when needed if(meopt==Terminate) { ScalarWaveFunction:: constructSpinInfo(const_ptr_cast(&inpart),incoming,true); ScalarWaveFunction:: constructSpinInfo(decay[iscal],outgoing,true); VectorWaveFunction:: constructSpinInfo(vector3_,decay[ivect],outgoing,true,false); VectorWaveFunction:: constructSpinInfo(gluon_, decay[iglu ],outgoing,true,false); return 0.; } // calculate colour factors and number of colour flows unsigned int nflow; vector cfactors = getColourFactors(inpart, decay, nflow); vector ME(nflow,new_ptr(GeneralDecayMatrixElement(PDT::Spin0, PDT::Spin0, PDT::Spin1, PDT::Spin1))); // create wavefunctions ScalarWaveFunction scal(decay[iscal]->momentum(), decay[iscal]->dataPtr(),outgoing); VectorWaveFunction::calculateWaveFunctions(vector3_,decay[ivect],outgoing,false); VectorWaveFunction::calculateWaveFunctions(gluon_, decay[iglu ],outgoing,true ); // gauge invariance test #ifdef GAUGE_CHECK gluon_.clear(); for(unsigned int ix=0;ix<3;++ix) { if(ix==1) gluon_.push_back(VectorWaveFunction()); else { gluon_.push_back(VectorWaveFunction(decay[iglu ]->momentum(), decay[iglu ]->dataPtr(),10, outgoing)); } } #endif if (! ((incomingVertex_[inter] && (outgoingVertexS_[inter] || outgoingVertexV_[inter])) || (outgoingVertexS_[inter] && outgoingVertexV_[inter]))) throw Exception() << "Invalid vertices for radiation in SSV decay in SSVDecayer::threeBodyME" << Exception::runerror; // sort out colour flows int S(1), V(2); if (decay[iscal]->dataPtr()->iColour()==PDT::Colour3bar && decay[ivect]->dataPtr()->iColour()==PDT::Colour8) swap(S,V); else if (decay[ivect]->dataPtr()->iColour()==PDT::Colour3 && decay[iscal]->dataPtr()->iColour()==PDT::Colour8) swap(S,V); Energy2 scale(sqr(inpart.mass())); const GeneralTwoBodyDecayer::CFlow & colourFlow = colourFlows(inpart, decay); double gs(0.); bool couplingSet(false); #ifdef GAUGE_CHECK double total=0.; #endif for(unsigned int iv = 0; iv < 3; ++iv) { for(unsigned int ig = 0; ig < 2; ++ig) { // radiation from the incoming scalar if((inpart.dataPtr()->coloured() && inter==ShowerInteraction::QCD) || (inpart.dataPtr()->charged() && inter==ShowerInteraction::QED) ) { assert(incomingVertex_[inter]); ScalarWaveFunction scalarInter = incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(), gluon_[2*ig],swave3_,inpart.mass()); assert(swave3_.particle()->id()==scalarInter.particle()->id()); Complex diag = 0.; for(auto vertex : vertex_) diag += vertex->evaluate(scale,vector3_[iv],scal,scalarInter); if(!couplingSet) { gs = abs(incomingVertex_[inter]->norm()); couplingSet = true; } for(unsigned int ix=0;ixdataPtr()->coloured() && inter==ShowerInteraction::QCD) || (decay[iscal]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) { assert(outgoingVertexS_[inter]); // ensure you get correct outgoing particle from first vertex tcPDPtr off = decay[iscal]->dataPtr(); if(off->CC()) off = off->CC(); ScalarWaveFunction scalarInter = outgoingVertexS_[inter]->evaluate(scale,3,off,gluon_[2*ig],scal,decay[iscal]->mass()); assert(scal.particle()->id()==scalarInter.particle()->id()); if(!couplingSet) { gs = abs(outgoingVertexS_[inter]->norm()); couplingSet = true; } Complex diag = 0.; for(auto vertex : vertex_) diag += vertex->evaluate(scale,vector3_[iv],scalarInter,swave3_); for(unsigned int ix=0;ixdataPtr()->coloured() && inter==ShowerInteraction::QCD) || (decay[ivect]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) { assert(outgoingVertexV_[inter]); // ensure you get correct outgoing particle from first vertex tcPDPtr off = decay[ivect]->dataPtr(); if(off->CC()) off = off->CC(); VectorWaveFunction vectorInter = outgoingVertexV_[inter]->evaluate(scale,3,off,gluon_[2*ig], vector3_[iv],decay[ivect]->mass()); assert(vector3_[iv].particle()->id()==vectorInter.particle()->id()); if(!couplingSet) { gs = abs(outgoingVertexV_[inter]->norm()); couplingSet = true; } Complex diag = 0.; for(auto vertex : vertex_) diag += vertex->evaluate(scale,vectorInter,scal,swave3_); for(unsigned int ix=0;ixevaluate(scale, gluon_[2*ig], vector3_[iv], scal, swave3_); for(unsigned int ix=0;ixcontract(*ME[iy],rho3_)).real(); } } // divide by alpha_(S,EM) output*=(4.*Constants::pi)/sqr(gs); #ifdef GAUGE_CHECK double ratio = output/total; if(abs(ratio)>1e-20) { generator()->log() << "Test of gauge invariance in decay\n" << inpart << "\n"; for(unsigned int ix=0;ixlog() << *decay[ix] << "\n"; generator()->log() << "Test of gauge invariance " << ratio << "\n"; } #endif // return the answer return output; } diff --git a/Decay/General/SSVDecayer.h b/Decay/General/SSVDecayer.h --- a/Decay/General/SSVDecayer.h +++ b/Decay/General/SSVDecayer.h @@ -1,233 +1,228 @@ // -*- C++ -*- // // SSVDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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_SSVDecayer_H #define HERWIG_SSVDecayer_H // // This is the declaration of the SSVDecayer class. // #include "GeneralTwoBodyDecayer.h" #include "ThePEG/Helicity/Vertex/Scalar/VSSVertex.h" #include "ThePEG/Helicity/Vertex/Vector/VVVVertex.h" #include "ThePEG/Helicity/Vertex/Scalar/VVSSVertex.h" #include "ThePEG/Repository/EventGenerator.h" namespace Herwig { using namespace ThePEG; using Helicity::VSSVertexPtr; /** \ingroup Decay * The SSVDecayer class implements the decay of a scalar to a vector * and a scalar in a general model. It holds an VSSVertex pointer * that must be typecast from the VertexBase pointer held in * GeneralTwoBodyDecayer. It implents the virtual functions me2() and * partialWidth(). * * @see GeneralTwoBodyDecayer */ class SSVDecayer: public GeneralTwoBodyDecayer { public: /** * The default constructor. */ SSVDecayer() {} /** @name Virtual functions required by the Decayer class. */ //@{ /** * Return the matrix element squared for a given mode and phase-space channel. * @param ichan The channel we are calculating the matrix element for. * @param part The decaying Particle. * @param outgoing The particles produced in the decay * @param momenta The momenta of the particles produced in the decay * @param meopt Option for the calculation of the matrix element * @return The matrix element squared for the phase-space configuration. */ double me2(const int ichan,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const; /** * Construct the SpinInfos for the particles produced in the decay */ virtual void constructSpinInfo(const Particle & part, ParticleVector outgoing) const; /** * Function to return partial Width * @param inpart The decaying particle. * @param outa One of the decay products. * @param outb The other decay product. */ virtual Energy partialWidth(PMPair inpart, PMPair outa, PMPair outb) const; /** * Has a POWHEG style correction */ virtual POWHEGType hasPOWHEGCorrection() { POWHEGType output = FSR; for(auto vertex : vertex_) { if(vertex->orderInAllCouplings()!=1) { output = No; break; } } return output; } /** * Three-body matrix element including additional QCD radiation */ virtual double threeBodyME(const int , const Particle & inpart, const ParticleVector & decay, ShowerInteraction inter, MEOption meopt); /** * Set the information on the decay */ virtual void setDecayInfo(PDPtr incoming, PDPair outgoing, vector, map &, const vector > &, map); //@} 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; //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ SSVDecayer & operator=(const SSVDecayer &) = delete; private: /** * Abstract pointer to AbstractFFVVertex */ vector vertex_; - - /** - * Pointer to the perturbative vertex - */ - vector perturbativeVertex_; /** * Abstract pointer to AbstractVSSVertex for QCD radiation from incoming scalar */ map incomingVertex_; /** * Abstract pointer to AbstractVSSVertex for QCD radiation from outgoing scalar */ map outgoingVertexS_; /** * Abstract pointer to AbstractVVVVertex for QCD radiation from outgoing vector */ map outgoingVertexV_; /** * Abstract pointer to AbstractVVSSVertex for QCD radiation from 4 point vertex */ map fourPointVertex_; /** * Spinor density matrix */ mutable RhoDMatrix rho_; /** * Scalar wavefunction */ mutable Helicity::ScalarWaveFunction swave_; /** * Vector wavefunction */ mutable vector vector_; /** * Spin density matrix for 3 body decay */ mutable RhoDMatrix rho3_; /** * Scalar wavefunction for 3 body decay */ mutable Helicity::ScalarWaveFunction swave3_; /** * Scalar wavefunction for 3 body decay */ mutable Helicity::ScalarWaveFunction scal_; /** * Vector wavefunction for 3 body decay */ mutable vector vector3_; /** * Vector wavefunction for 3 body decay */ mutable vector gluon_; }; } #endif /* HERWIG_SSVDecayer_H */ diff --git a/MatrixElement/Matchbox/Matching/ShowerApproximation.cc b/MatrixElement/Matchbox/Matching/ShowerApproximation.cc --- a/MatrixElement/Matchbox/Matching/ShowerApproximation.cc +++ b/MatrixElement/Matchbox/Matching/ShowerApproximation.cc @@ -1,716 +1,718 @@ // -*- C++ -*- // // ShowerApproximation.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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 ShowerApproximation class. // #include "ShowerApproximation.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h" #include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h" #include "Herwig/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h" using namespace Herwig; ShowerApproximation::ShowerApproximation() : HandlerBase(), theExtrapolationX(1.0), theBelowCutoff(false), theFFPtCut(1.0*GeV), theFFScreeningScale(ZERO), theFIPtCut(1.0*GeV), theFIScreeningScale(ZERO), theIIPtCut(1.0*GeV), theIIScreeningScale(ZERO), theSafeCut(0.0*GeV), theRestrictPhasespace(true), theHardScaleFactor(1.0), theRenormalizationScaleFactor(1.0), theFactorizationScaleFactor(1.0), theRealEmissionScaleInSubtraction(showerScale), theBornScaleInSubtraction(showerScale), theEmissionScaleInSubtraction(showerScale), theRealEmissionScaleInSplitting(showerScale), theBornScaleInSplitting(showerScale), theEmissionScaleInSplitting(showerScale), theRenormalizationScaleFreeze(1.*GeV), theFactorizationScaleFreeze(1.*GeV), maxPtIsMuF(false), theOpenZ(true) {} ShowerApproximation::~ShowerApproximation() {} void ShowerApproximation::setLargeNBasis() { assert(dipole()->realEmissionME()->matchboxAmplitude()); if ( !dipole()->realEmissionME()->matchboxAmplitude()->treeAmplitudes() ) return; if ( !theLargeNBasis ) { if ( !dipole()->realEmissionME()->matchboxAmplitude()->colourBasis() ) throw Exception() << "ShowerApproximation::setLargeNBasis(): Expecting a colour basis object." << Exception::runerror; theLargeNBasis = dipole()->realEmissionME()->matchboxAmplitude()->colourBasis()->cloneMe(); theLargeNBasis->clear(); theLargeNBasis->doLargeN(); } } void ShowerApproximation::setDipole(Ptr::tptr dip) { theDipole = dip; setLargeNBasis(); } Ptr::tptr ShowerApproximation::dipole() const { return theDipole; } Ptr::tptr ShowerApproximation::showerTildeKinematics() const { return Ptr::tptr(); } Ptr::tptr ShowerApproximation::showerInvertedTildeKinematics() const { return Ptr::tptr(); } void ShowerApproximation::checkCutoff() { assert(!showerTildeKinematics()); } void ShowerApproximation::getShowerVariables() { // check for the cutoff dipole()->isAboveCutoff(isAboveCutoff()); // get the hard scale dipole()->showerHardScale(hardScale()); // set the shower scale and variables for completeness dipole()->showerScale(dipole()->lastPt()); dipole()->showerParameters().resize(1); dipole()->showerParameters()[0] = dipole()->lastZ(); // check for phase space dipole()->isInShowerPhasespace(isInShowerPhasespace()); } bool ShowerApproximation::isAboveCutoff() const { if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() > 1 ) { return dipole()->lastPt() >= max(ffPtCut(),safeCut()); } else if ( ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() < 2 ) || ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() > 1 ) ) { return dipole()->lastPt() >= max(fiPtCut(),safeCut()); } else { assert(dipole()->bornEmitter() < 2 && dipole()->bornSpectator() < 2); return dipole()->lastPt() >= max(iiPtCut(),safeCut()); } return true; } Energy ShowerApproximation::hardScale() const { if ( !maxPtIsMuF ) { if ( !bornCXComb()->mePartonData()[0]->coloured() && !bornCXComb()->mePartonData()[1]->coloured() ) { Energy maxPt = (bornCXComb()->meMomenta()[0] + bornCXComb()->meMomenta()[1]).m(); maxPt *= hardScaleFactor(); return maxPt; } Energy maxPt = generator()->maximumCMEnergy(); vector::const_iterator p = bornCXComb()->meMomenta().begin() + 2; cPDVector::const_iterator pp = bornCXComb()->mePartonData().begin() + 2; for ( ; p != bornCXComb()->meMomenta().end(); ++p, ++pp ) if ( (**pp).coloured() ) maxPt = min(maxPt,p->mt()); if ( maxPt == generator()->maximumCMEnergy() ) maxPt = (bornCXComb()->meMomenta()[0] + bornCXComb()->meMomenta()[1]).m(); maxPt *= hardScaleFactor(); return maxPt; } else { return hardScaleFactor()*sqrt(bornCXComb()->lastShowerScale()); } } bool ShowerApproximation::isInShowerPhasespace() const { if ( !dipole()->isAboveCutoff() ) return false; if ( !restrictPhasespace() ) return true; InvertedTildeKinematics& kinematics = const_cast(*dipole()->invertedTildeKinematics()); tcStdXCombPtr tmpreal = kinematics.realXComb(); tcStdXCombPtr tmpborn = kinematics.bornXComb(); Ptr::tptr tmpdip = kinematics.dipole(); Energy hard = dipole()->showerHardScale(); Energy pt = dipole()->lastPt(); double z = dipole()->lastZ(); pair zbounds(0.,1.); kinematics.dipole(const_ptr_cast::tptr>(theDipole)); kinematics.prepare(realCXComb(),bornCXComb()); if ( pt > hard ) { kinematics.dipole(tmpdip); kinematics.prepare(tmpreal,tmpborn); return false; } try { zbounds = kinematics.zBounds(pt,openZ() ? kinematics.ptMax() : hard); } catch(...) { kinematics.dipole(tmpdip); kinematics.prepare(tmpreal,tmpborn); throw; } kinematics.dipole(tmpdip); kinematics.prepare(tmpreal,tmpborn); return z > zbounds.first && z < zbounds.second; } Energy2 ShowerApproximation::showerEmissionScale() const { Energy2 mur = sqr(dipole()->lastPt()); if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() > 1 ) { return mur + sqr(ffScreeningScale()); } else if ( ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() < 2 ) || ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() > 1 ) ) { return mur + sqr(fiScreeningScale()); } else { assert(dipole()->bornEmitter() < 2 && dipole()->bornSpectator() < 2); return mur + sqr(iiScreeningScale()); } return mur; } Energy2 ShowerApproximation::bornRenormalizationScale() const { return sqr(dipole()->underlyingBornME()->renormalizationScaleFactor()) * dipole()->underlyingBornME()->renormalizationScale(); } Energy2 ShowerApproximation::bornFactorizationScale() const { return sqr(dipole()->underlyingBornME()->factorizationScaleFactor()) * dipole()->underlyingBornME()->factorizationScale(); } Energy2 ShowerApproximation::realRenormalizationScale() const { return sqr(dipole()->realEmissionME()->renormalizationScaleFactor()) * dipole()->realEmissionME()->renormalizationScale(); } Energy2 ShowerApproximation::realFactorizationScale() const { return sqr(dipole()->realEmissionME()->factorizationScaleFactor()) * dipole()->realEmissionME()->factorizationScale(); } double ShowerApproximation::bornPDFWeight(Energy2 muf) const { if ( !bornCXComb()->mePartonData()[0]->coloured() && !bornCXComb()->mePartonData()[1]->coloured() ) return 1.; if ( muf < sqr(theFactorizationScaleFreeze) ) muf = sqr(theFactorizationScaleFreeze); double pdfweight = 1.; if ( bornCXComb()->mePartonData()[0]->coloured() && dipole()->underlyingBornME()->havePDFWeight1() ) pdfweight *= dipole()->underlyingBornME()->pdf1(muf,theExtrapolationX); if ( bornCXComb()->mePartonData()[1]->coloured() && dipole()->underlyingBornME()->havePDFWeight2() ) pdfweight *= dipole()->underlyingBornME()->pdf2(muf,theExtrapolationX); return pdfweight; } double ShowerApproximation::realPDFWeight(Energy2 muf) const { if ( !realCXComb()->mePartonData()[0]->coloured() && !realCXComb()->mePartonData()[1]->coloured() ) return 1.; if ( muf < sqr(theFactorizationScaleFreeze) ) muf = sqr(theFactorizationScaleFreeze); double pdfweight = 1.; if ( realCXComb()->mePartonData()[0]->coloured() && dipole()->realEmissionME()->havePDFWeight1() ) pdfweight *= dipole()->realEmissionME()->pdf1(muf,theExtrapolationX); if ( realCXComb()->mePartonData()[1]->coloured() && dipole()->realEmissionME()->havePDFWeight2() ) pdfweight *= dipole()->realEmissionME()->pdf2(muf,theExtrapolationX); return pdfweight; } double ShowerApproximation::scaleWeight(int rScale, int bScale, int eScale) const { double emissionAlpha = 1.; Energy2 emissionScale = ZERO; Energy2 showerscale = ZERO; if ( eScale == showerScale || bScale == showerScale || eScale == showerScale ) { showerscale = showerRenormalizationScale(); if ( showerscale < sqr(theRenormalizationScaleFreeze) ) showerscale = sqr(theFactorizationScaleFreeze); } if ( eScale == showerScale ) { emissionAlpha = SM().alphaS(showerscale); emissionScale = showerFactorizationScale(); } else if ( eScale == realScale ) { emissionAlpha = dipole()->realEmissionME()->lastXComb().lastAlphaS(); emissionScale = dipole()->realEmissionME()->lastScale(); } else if ( eScale == bornScale ) { emissionAlpha = dipole()->underlyingBornME()->lastXComb().lastAlphaS(); emissionScale = dipole()->underlyingBornME()->lastScale(); } double emissionPDF = realPDFWeight(emissionScale); double couplingFactor = 1.; if ( bScale != rScale ) { double bornAlpha = 1.; if ( bScale == showerScale ) { bornAlpha = SM().alphaS(showerscale); } else if ( bScale == realScale ) { bornAlpha = dipole()->realEmissionME()->lastXComb().lastAlphaS(); } else if ( bScale == bornScale ) { bornAlpha = dipole()->underlyingBornME()->lastXComb().lastAlphaS(); } double realAlpha = 1.; if ( rScale == showerScale ) { realAlpha = SM().alphaS(showerscale); } else if ( rScale == realScale ) { realAlpha = dipole()->realEmissionME()->lastXComb().lastAlphaS(); } else if ( rScale == bornScale ) { realAlpha = dipole()->underlyingBornME()->lastXComb().lastAlphaS(); } couplingFactor *= pow(realAlpha/bornAlpha,(double)(dipole()->underlyingBornME()->orderInAlphaS())); } Energy2 hardScale = ZERO; if ( bScale == showerScale ) { hardScale = showerFactorizationScale(); } else if ( bScale == realScale ) { hardScale = dipole()->realEmissionME()->lastScale(); } else if ( bScale == bornScale ) { hardScale = dipole()->underlyingBornME()->lastScale(); } double bornPDF = bornPDFWeight(hardScale); - if ( bornPDF < 1e-8 ) + if ( abs(bornPDF) < 1e-8 ) bornPDF = 0.0; - if ( emissionPDF < 1e-8 ) + if ( abs(emissionPDF) < 1e-8 ) emissionPDF = 0.0; if ( emissionPDF == 0.0 || bornPDF == 0.0 ) return 0.0; + double pdfRatio = emissionPDF/bornPDF; + pdfRatio = min(abs(pdfRatio),100000.); + return - emissionAlpha * emissionPDF * - couplingFactor / bornPDF; + emissionAlpha * pdfRatio * couplingFactor; } double ShowerApproximation::channelWeight(int emitter, int emission, int spectator, int) const { double cfac = 1.; double Nc = generator()->standardModel()->Nc(); if (realCXComb()->mePartonData()[emitter]->iColour() == PDT::Colour8){ if (realCXComb()->mePartonData()[emission]->iColour() == PDT::Colour8) cfac = Nc; else if ( realCXComb()->mePartonData()[emission]->iColour() == PDT::Colour3 || realCXComb()->mePartonData()[emission]->iColour() == PDT::Colour3bar) cfac = 0.5; else assert(false); } else if ((realCXComb()->mePartonData()[emitter] ->iColour() == PDT::Colour3 || realCXComb()->mePartonData()[emitter] ->iColour() == PDT::Colour3bar)) cfac = (sqr(Nc)-1.)/(2.*Nc); else assert(false); // do the most simple thing for the time being; needs fixing later if ( realCXComb()->mePartonData()[emission]->id() == ParticleID::g ) { Energy2 pipk = realCXComb()->meMomenta()[emitter] * realCXComb()->meMomenta()[spectator]; Energy2 pipj = realCXComb()->meMomenta()[emitter] * realCXComb()->meMomenta()[emission]; Energy2 pjpk = realCXComb()->meMomenta()[emission] * realCXComb()->meMomenta()[spectator]; return cfac *GeV2 * pipk / ( pipj * ( pipj + pjpk ) ); } return cfac * GeV2 / (realCXComb()->meMomenta()[emitter] * realCXComb()->meMomenta()[emission]); } double ShowerApproximation::channelWeight() const { double currentChannel = channelWeight(dipole()->realEmitter(), dipole()->realEmission(), dipole()->realSpectator(), dipole()->bornEmitter()); if ( currentChannel == 0. ) return 0.; double sum = 0.; for ( vector::tptr>::const_iterator dip = dipole()->partnerDipoles().begin(); dip != dipole()->partnerDipoles().end(); ++dip ) sum += channelWeight((**dip).realEmitter(), (**dip).realEmission(), (**dip).realSpectator(), (**dip).bornEmitter()); assert(sum > 0.0); return currentChannel / sum; } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void ShowerApproximation::doinit() { if ( profileScales() ) { if ( profileScales()->unrestrictedPhasespace() && restrictPhasespace() ) { generator()->log() << "ShowerApproximation warning: The scale profile chosen requires an unrestricted phase space,\n" << "however, the phase space was set to be restricted. Will switch to unrestricted phase space.\n" << flush; restrictPhasespace(false); } } HandlerBase::doinit(); } void ShowerApproximation::persistentOutput(PersistentOStream & os) const { os << theLargeNBasis << theBornXComb << theRealXComb << theTildeXCombs << theDipole << theBelowCutoff << ounit(theFFPtCut,GeV) << ounit(theFFScreeningScale,GeV) << ounit(theFIPtCut,GeV) << ounit(theFIScreeningScale,GeV) << ounit(theIIPtCut,GeV) << ounit(theIIScreeningScale,GeV) << ounit(theSafeCut,GeV) << theRestrictPhasespace << theHardScaleFactor << theRenormalizationScaleFactor << theFactorizationScaleFactor << theExtrapolationX << theRealEmissionScaleInSubtraction << theBornScaleInSubtraction << theEmissionScaleInSubtraction << theRealEmissionScaleInSplitting << theBornScaleInSplitting << theEmissionScaleInSplitting << ounit(theRenormalizationScaleFreeze,GeV) << ounit(theFactorizationScaleFreeze,GeV) << maxPtIsMuF << theHardScaleProfile << theOpenZ; } void ShowerApproximation::persistentInput(PersistentIStream & is, int) { is >> theLargeNBasis >> theBornXComb >> theRealXComb >> theTildeXCombs >> theDipole >> theBelowCutoff >> iunit(theFFPtCut,GeV) >> iunit(theFFScreeningScale,GeV) >> iunit(theFIPtCut,GeV) >> iunit(theFIScreeningScale,GeV) >> iunit(theIIPtCut,GeV) >> iunit(theIIScreeningScale,GeV) >> iunit(theSafeCut,GeV) >> theRestrictPhasespace >> theHardScaleFactor >> theRenormalizationScaleFactor >> theFactorizationScaleFactor >> theExtrapolationX >> theRealEmissionScaleInSubtraction >> theBornScaleInSubtraction >> theEmissionScaleInSubtraction >> theRealEmissionScaleInSplitting >> theBornScaleInSplitting >> theEmissionScaleInSplitting >> iunit(theRenormalizationScaleFreeze,GeV) >> iunit(theFactorizationScaleFreeze,GeV) >> maxPtIsMuF >> theHardScaleProfile >> theOpenZ; } // *** 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). DescribeAbstractClass describeHerwigShowerApproximation("Herwig::ShowerApproximation", "Herwig.so"); void ShowerApproximation::Init() { static ClassDocumentation documentation ("ShowerApproximation describes the shower emission to be used " "in NLO matching."); static Parameter interfaceFFPtCut ("FFPtCut", "Set the pt infrared cutoff", &ShowerApproximation::theFFPtCut, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); static Parameter interfaceFIPtCut ("FIPtCut", "Set the pt infrared cutoff", &ShowerApproximation::theFIPtCut, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); static Parameter interfaceIIPtCut ("IIPtCut", "Set the pt infrared cutoff", &ShowerApproximation::theIIPtCut, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); static Parameter interfaceSafeCut ("SafeCut", "Set the enhanced infrared cutoff for the Matching.", &ShowerApproximation::theSafeCut, GeV, 0.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); static Parameter interfaceFFScreeningScale ("FFScreeningScale", "Set the screening scale", &ShowerApproximation::theFFScreeningScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); static Parameter interfaceFIScreeningScale ("FIScreeningScale", "Set the screening scale", &ShowerApproximation::theFIScreeningScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); static Parameter interfaceIIScreeningScale ("IIScreeningScale", "Set the screening scale", &ShowerApproximation::theIIScreeningScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); static Switch interfaceRestrictPhasespace ("RestrictPhasespace", "Switch on or off phasespace restrictions", &ShowerApproximation::theRestrictPhasespace, true, false, false); static SwitchOption interfaceRestrictPhasespaceYes (interfaceRestrictPhasespace, "Yes", "Perform phasespace restrictions", true); static SwitchOption interfaceRestrictPhasespaceNo (interfaceRestrictPhasespace, "No", "Do not perform phasespace restrictions", false); static Parameter interfaceHardScaleFactor ("HardScaleFactor", "The hard scale factor.", &ShowerApproximation::theHardScaleFactor, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceRenormalizationScaleFactor ("RenormalizationScaleFactor", "The hard scale factor.", &ShowerApproximation::theRenormalizationScaleFactor, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceFactorizationScaleFactor ("FactorizationScaleFactor", "The hard scale factor.", &ShowerApproximation::theFactorizationScaleFactor, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceExtrapolationX ("ExtrapolationX", "The x from which on extrapolation should be performed.", &ShowerApproximation::theExtrapolationX, 1.0, 0.0, 1.0, false, false, Interface::limited); static Switch interfaceRealEmissionScaleInSubtraction ("RealEmissionScaleInSubtraction", "Set the scale choice for the real emission cross section in the matching subtraction.", &ShowerApproximation::theRealEmissionScaleInSubtraction, showerScale, false, false); static SwitchOption interfaceRealEmissionScaleInSubtractionRealScale (interfaceRealEmissionScaleInSubtraction, "RealScale", "Use the real emission scale.", realScale); static SwitchOption interfaceRealEmissionScaleInSubtractionBornScale (interfaceRealEmissionScaleInSubtraction, "BornScale", "Use the Born scale.", bornScale); static SwitchOption interfaceRealEmissionScaleInSubtractionShowerScale (interfaceRealEmissionScaleInSubtraction, "ShowerScale", "Use the shower scale", showerScale); interfaceRealEmissionScaleInSubtraction.rank(-1); static Switch interfaceBornScaleInSubtraction ("BornScaleInSubtraction", "Set the scale choice for the Born cross section in the matching subtraction.", &ShowerApproximation::theBornScaleInSubtraction, showerScale, false, false); static SwitchOption interfaceBornScaleInSubtractionRealScale (interfaceBornScaleInSubtraction, "RealScale", "Use the real emission scale.", realScale); static SwitchOption interfaceBornScaleInSubtractionBornScale (interfaceBornScaleInSubtraction, "BornScale", "Use the Born scale.", bornScale); static SwitchOption interfaceBornScaleInSubtractionShowerScale (interfaceBornScaleInSubtraction, "ShowerScale", "Use the shower scale", showerScale); interfaceBornScaleInSubtraction.rank(-1); static Switch interfaceEmissionScaleInSubtraction ("EmissionScaleInSubtraction", "Set the scale choice for the emission in the matching subtraction.", &ShowerApproximation::theEmissionScaleInSubtraction, showerScale, false, false); static SwitchOption interfaceEmissionScaleInSubtractionRealScale (interfaceEmissionScaleInSubtraction, "RealScale", "Use the real emission scale.", realScale); static SwitchOption interfaceEmissionScaleInSubtractionEmissionScale (interfaceEmissionScaleInSubtraction, "BornScale", "Use the Born scale.", bornScale); static SwitchOption interfaceEmissionScaleInSubtractionShowerScale (interfaceEmissionScaleInSubtraction, "ShowerScale", "Use the shower scale", showerScale); interfaceEmissionScaleInSubtraction.rank(-1); static Switch interfaceRealEmissionScaleInSplitting ("RealEmissionScaleInSplitting", "Set the scale choice for the real emission cross section in the splitting.", &ShowerApproximation::theRealEmissionScaleInSplitting, showerScale, false, false); static SwitchOption interfaceRealEmissionScaleInSplittingRealScale (interfaceRealEmissionScaleInSplitting, "RealScale", "Use the real emission scale.", realScale); static SwitchOption interfaceRealEmissionScaleInSplittingBornScale (interfaceRealEmissionScaleInSplitting, "BornScale", "Use the Born scale.", bornScale); static SwitchOption interfaceRealEmissionScaleInSplittingShowerScale (interfaceRealEmissionScaleInSplitting, "ShowerScale", "Use the shower scale", showerScale); interfaceRealEmissionScaleInSplitting.rank(-1); static Switch interfaceBornScaleInSplitting ("BornScaleInSplitting", "Set the scale choice for the Born cross section in the splitting.", &ShowerApproximation::theBornScaleInSplitting, showerScale, false, false); static SwitchOption interfaceBornScaleInSplittingRealScale (interfaceBornScaleInSplitting, "RealScale", "Use the real emission scale.", realScale); static SwitchOption interfaceBornScaleInSplittingBornScale (interfaceBornScaleInSplitting, "BornScale", "Use the Born scale.", bornScale); static SwitchOption interfaceBornScaleInSplittingShowerScale (interfaceBornScaleInSplitting, "ShowerScale", "Use the shower scale", showerScale); interfaceBornScaleInSplitting.rank(-1); static Switch interfaceEmissionScaleInSplitting ("EmissionScaleInSplitting", "Set the scale choice for the emission in the splitting.", &ShowerApproximation::theEmissionScaleInSplitting, showerScale, false, false); static SwitchOption interfaceEmissionScaleInSplittingRealScale (interfaceEmissionScaleInSplitting, "RealScale", "Use the real emission scale.", realScale); static SwitchOption interfaceEmissionScaleInSplittingEmissionScale (interfaceEmissionScaleInSplitting, "BornScale", "Use the Born scale.", bornScale); static SwitchOption interfaceEmissionScaleInSplittingShowerScale (interfaceEmissionScaleInSplitting, "ShowerScale", "Use the shower scale", showerScale); interfaceEmissionScaleInSplitting.rank(-1); static Parameter interfaceRenormalizationScaleFreeze ("RenormalizationScaleFreeze", "The freezing scale for the renormalization scale.", &ShowerApproximation::theRenormalizationScaleFreeze, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); interfaceRenormalizationScaleFreeze.rank(-1); static Parameter interfaceFactorizationScaleFreeze ("FactorizationScaleFreeze", "The freezing scale for the factorization scale.", &ShowerApproximation::theFactorizationScaleFreeze, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); interfaceFactorizationScaleFreeze.rank(-1); static Reference interfaceHardScaleProfile ("HardScaleProfile", "The hard scale profile to use.", &ShowerApproximation::theHardScaleProfile, false, false, true, true, false); static Reference interfaceLargeNBasis ("LargeNBasis", "Set the large-N colour basis implementation.", &ShowerApproximation::theLargeNBasis, false, false, true, true, false); interfaceLargeNBasis.rank(-1); static Switch interfaceMaxPtIsMuF ("MaxPtIsMuF", "", &ShowerApproximation::maxPtIsMuF, false, false, false); static SwitchOption interfaceMaxPtIsMuFYes (interfaceMaxPtIsMuF, "Yes", "", true); static SwitchOption interfaceMaxPtIsMuFNo (interfaceMaxPtIsMuF, "No", "", false); }