diff --git a/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc b/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc --- a/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc +++ b/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc @@ -1,222 +1,243 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the HalfHalfZeroEWSplitFn class. // #include "HalfHalfZeroEWSplitFn.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/ParticleData.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" +#include "Herwig/Models/StandardModel/SMFFHVertex.h" using namespace Herwig; IBPtr HalfHalfZeroEWSplitFn::clone() const { return new_ptr(*this); } IBPtr HalfHalfZeroEWSplitFn::fullclone() const { return new_ptr(*this); } void HalfHalfZeroEWSplitFn::persistentOutput(PersistentOStream & os) const { os << ghqq_; } void HalfHalfZeroEWSplitFn::persistentInput(PersistentIStream & is, int) { is >> ghqq_; } // The following static variable is needed for the type description system in ThePEG. DescribeClass describeHerwigHalfHalfZeroEWSplitFn("Herwig::HalfHalfZeroEWSplitFn", "HwShower.so"); void HalfHalfZeroEWSplitFn::Init() { static ClassDocumentation documentation ("The HalfHalfZeroEWSplitFn class implements the splitting q->qH"); } void HalfHalfZeroEWSplitFn::doinit() { SplittingFunction::doinit(); tcSMPtr sm = generator()->standardModel(); double sw2 = sm->sin2ThetaW(); ghqq_ = 1./sqrt(4.*sw2); + + _theSM = dynamic_ptr_cast(generator()->standardModel()); } void HalfHalfZeroEWSplitFn::getCouplings(double & gH, const IdList & ids) const { if(abs(ids[2]->id())==ParticleID::h0) { //get quark masses Energy mq; if(abs(ids[0]->id())==ParticleID::c) mq = getParticleData(ParticleID::c)->mass(); else if(abs(ids[0]->id())==ParticleID::b) mq = getParticleData(ParticleID::b)->mass(); else if(abs(ids[0]->id())==ParticleID::t) mq = getParticleData(ParticleID::t)->mass(); Energy mW = getParticleData(ParticleID::Wplus)->mass(); gH = ghqq_*(mq/mW); } else assert(false); } +void HalfHalfZeroEWSplitFn::getCouplings(double & gH, const IdList & ids, const Energy2 t) const { + if(abs(ids[2]->id())==ParticleID::h0) { + //get quark masses + Energy mq; + if(abs(ids[0]->id())==ParticleID::c) + mq = _theSM->mass(t,getParticleData(ParticleID::c)); + else if(abs(ids[0]->id())==ParticleID::b) + mq = _theSM->mass(t,getParticleData(ParticleID::b)); + else if(abs(ids[0]->id())==ParticleID::t) + mq = _theSM->mass(t,getParticleData(ParticleID::t)); + Energy mW = getParticleData(ParticleID::Wplus)->mass(); + //Energy mW = _theSM->mass(t,getParticleData(ParticleID::Wplus)); + gH = ghqq_*(mq/mW); + } + else + assert(false); +} + double HalfHalfZeroEWSplitFn::P(const double z, const Energy2 t, const IdList &ids, const bool mass, const RhoDMatrix & rho) const { double gH(0.); - getCouplings(gH,ids); + getCouplings(gH,ids,t); double val = (1.-z); Energy mq, mH; //get masses if(mass) { mq = ids[0]->mass(); mH = ids[2]->mass(); } else { // to assure the particle mass in non-zero if(abs(ids[0]->id())==ParticleID::c) mq = getParticleData(ParticleID::c)->mass(); else if(abs(ids[0]->id())==ParticleID::b) mq = getParticleData(ParticleID::b)->mass(); else if(abs(ids[0]->id())==ParticleID::t) mq = getParticleData(ParticleID::t)->mass(); mH = getParticleData(ParticleID::h0)->mass(); } val += (4.*sqr(mq) - sqr(mH))/(t*(1. - z)*z); val *= sqr(gH); return colourFactor(ids)*val; } double HalfHalfZeroEWSplitFn::overestimateP(const double z, const IdList & ids) const { double gH(0.); getCouplings(gH,ids); return sqr(gH)*colourFactor(ids)*(1.-z); } double HalfHalfZeroEWSplitFn::ratioP(const double z, const Energy2 t, const IdList & ids, const bool mass, const RhoDMatrix & rho) const { double gH(0.); - getCouplings(gH,ids); + getCouplings(gH,ids,t); double val = 1.; Energy mq, mH; if(mass) { mq = ids[0]->mass(); mH = ids[2]->mass(); } else { // to assure the particle mass in non-zero if(abs(ids[0]->id())==ParticleID::c) mq = getParticleData(ParticleID::c)->mass(); else if(abs(ids[0]->id())==ParticleID::b) mq = getParticleData(ParticleID::b)->mass(); else if(abs(ids[0]->id())==ParticleID::t) mq = getParticleData(ParticleID::t)->mass(); mH = getParticleData(ParticleID::h0)->mass(); } val += (4.*sqr(mq) - sqr(mH))/(t*(1. - z)*z); return val; } double HalfHalfZeroEWSplitFn::integOverP(const double z, const IdList & ids, unsigned int PDFfactor) const { double gH(0.); getCouplings(gH,ids); double pre = colourFactor(ids)*sqr(gH); switch (PDFfactor) { case 0: //OverP return pre*(z-sqr(z)/2.); case 1: //OverP/z return pre*(log(z)-z); case 2: //OverP/(1-z) return pre*z; case 3: //OverP/[z(1-z)] return pre*log(z); default: throw Exception() << "HalfHalfZeroEWSplitFn::integOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } double HalfHalfZeroEWSplitFn::invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor) const { double gH(0.); getCouplings(gH,ids); double pre = colourFactor(ids)*sqr(gH); switch (PDFfactor) { case 0: return max((-pre+sqrt(sqr(pre)-2.*pre*r))/pre, (-pre-sqrt(sqr(pre)-2.*pre*r))/pre); case 1: //OverP/z case 2: //OverP/(1-z) return r/pre; case 3: //OverP/[z(1-z)] return exp(r/pre); default: throw Exception() << "HalfHalfZeroEWSplitFn::invIntegOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } bool HalfHalfZeroEWSplitFn::accept(const IdList &ids) const { if(ids.size()!=3) return false; if(ids[2]->id()==ParticleID::h0) { if(ids[0]->id()==ids[1]->id() && (ids[0]->id()==4 || ids[0]->id()==5 || ids[0]->id()==6)) return true; } return false; } vector > HalfHalfZeroEWSplitFn::generatePhiForward(const double, const Energy2, const IdList & , const RhoDMatrix &) { // no dependence on the spin density matrix, dependence on off-diagonal terms cancels // and rest = splitting function for Tr(rho)=1 as required by defn return vector >(1,make_pair(0,1.)); } vector > HalfHalfZeroEWSplitFn::generatePhiBackward(const double, const Energy2, const IdList & , const RhoDMatrix &) { // no dependence on the spin density matrix, dependence on off-diagonal terms cancels // and rest = splitting function for Tr(rho)=1 as required by defn return vector >(1,make_pair(0,1.)); } DecayMEPtr HalfHalfZeroEWSplitFn::matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool) { // calculate the kernal DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin0))); //get masses Energy mq, mH; if(abs(ids[0]->id())==ParticleID::c) mq = getParticleData(ParticleID::c)->mass(); else if(abs(ids[0]->id())==ParticleID::b) mq = getParticleData(ParticleID::b)->mass(); else if(abs(ids[0]->id())==ParticleID::t) mq = getParticleData(ParticleID::t)->mass(); mH = getParticleData(ParticleID::h0)->mass(); double gH(0.); - getCouplings(gH,ids); + getCouplings(gH,ids,t); double mqt = mq/sqrt(t); double mHt = mH/sqrt(t); double num1 = gH*(1.+z)*mqt; double num2 = gH*sqrt(-sqr(mqt)*(1.-z) - sqr(mHt)*z + z*(1.-z)*(sqr(mqt)+z*(1.-z))); //watch this double dnum = sqrt(2.)*sqrt((1.-z)*sqr(z)); Complex phase = exp(Complex(0.,1.)*phi); Complex cphase = conj(phase); (*kernal)(0,0,0) = num1/dnum; (*kernal)(0,1,0) = cphase*num2/dnum; (*kernal)(1,0,0) = -phase*num2/dnum; (*kernal)(1,1,0) = num1/dnum; // return the answer return kernal; } diff --git a/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.h b/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.h --- a/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.h +++ b/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.h @@ -1,207 +1,219 @@ // -*- C++ -*- #ifndef Herwig_HalfHalfZeroEWSplitFn_H #define Herwig_HalfHalfZeroEWSplitFn_H // // This is the declaration of the HalfHalfZeroEWSplitFn class. // #include "SplittingFunction.h" +#include "Herwig/Models/StandardModel/StandardModel.h" namespace Herwig { using namespace ThePEG; /** * The HalfHalfZeroEWSplitFn class implements the splitting function for * \f$\frac12\to q\frac12 h\f$ where the spin-0 higgs particle is a massive scalar boson. * * @see \ref HalfHalfZeroEWSplitFnInterfaces "The interfaces" * defined for HalfHalfZeroEWSplitFn. */ class HalfHalfZeroEWSplitFn: public SplittingFunction { public: /** * Concrete implementation of the method to determine whether this splitting * function can be used for a given set of particles. * @param ids The PDG codes for the particles in the splitting. */ virtual bool accept(const IdList & ids) const; /** * Methods to return the splitting function. */ //@{ /** * The concrete implementation of the splitting function, \f$P(z,t)\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms * @param rho The spin density matrix */ virtual double P(const double z, const Energy2 t, const IdList & ids, const bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the overestimate of the splitting function, * \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. */ virtual double overestimateP(const double z, const IdList & ids) const; /** * The concrete implementation of the * the ratio of the splitting function to the overestimate, i.e. * \f$P(z,t)/P_{\rm over}(z)\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms * @param rho The spin density matrix */ virtual double ratioP(const double z, const Energy2 t, const IdList & ids, const bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the indefinite integral of the * overestimated splitting function, \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double integOverP(const double z, const IdList & ids, unsigned int PDFfactor=0) const; /** * The concrete implementation of the inverse of the indefinite integral. * @param r Value of the splitting function to be inverted * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor=0) const; //@} /** * Method to calculate the azimuthal angle * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiForward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Method to calculate the azimuthal angle for backward evolution * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiBackward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Calculate the matrix element for the splitting * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. */ virtual DecayMEPtr matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool timeLike); protected: /** - * Get the couplings + * Get the couplings without running masses */ void getCouplings(double & gH, const IdList & ids) const; + /** + * Get the couplings with running masses + */ + void getCouplings(double & gH, const IdList & ids, const Energy2 t) 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; //@} 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(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ HalfHalfZeroEWSplitFn & operator=(const HalfHalfZeroEWSplitFn &) = delete; private: /** * Higgs couplings */ double ghqq_; + + /** + * Pointer to the SM object. + */ + tcHwSMPtr _theSM; + }; } #endif /* Herwig_HalfHalfZeroEWSplitFn_H */