diff --git a/MatrixElement/Lepton/MEee2Higgs2SM.cc b/MatrixElement/Lepton/MEee2Higgs2SM.cc --- a/MatrixElement/Lepton/MEee2Higgs2SM.cc +++ b/MatrixElement/Lepton/MEee2Higgs2SM.cc @@ -1,354 +1,422 @@ // -*- C++ -*- // // MEee2Higgs2SM.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 MEee2Higgs2SM class. // #include "MEee2Higgs2SM.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Handlers/StandardXComb.h" #include "Herwig/MatrixElement/HardVertex.h" #include "ThePEG/PDT/EnumParticles.h" using namespace Herwig; void MEee2Higgs2SM::doinit() { ME2to2Base::doinit(); h0_ = getParticleData(ThePEG::ParticleID::h0); tcHwSMPtr hwsm=ThePEG::dynamic_ptr_cast(standardModel()); // do the initialisation if(hwsm) { FFHVertex_ = hwsm->vertexFFH(); HGGVertex_ = hwsm->vertexHGG(); + HWWVertex_ = hwsm->vertexWWH(); } else { throw InitException() << "Must have Herwig StandardModel object in " << "MEee2Higgs2SM::doinit()" << Exception::runerror; } } Energy2 MEee2Higgs2SM::scale() const { return sHat(); } unsigned int MEee2Higgs2SM::orderInAlphaS() const { return 0; } unsigned int MEee2Higgs2SM::orderInAlphaEW() const { return 2; } void MEee2Higgs2SM::getDiagrams() const { // specific the diagrams tcPDPtr ep = getParticleData(ParticleID::eplus ); tcPDPtr em = getParticleData(ParticleID::eminus); // outgoing quarks for(int i=1;i<=5;++i) { if(allowed_==0 || allowed_==1 || allowed_==i+2) { tcPDPtr lm = getParticleData(i); tcPDPtr lp = lm->CC(); add(new_ptr((Tree2toNDiagram(2), em, ep, 1, h0_, 3, lm, 3, lp, -1))); } } // outgoing leptons for( int i =11;i<=16;i+=2) { if(allowed_==0 || allowed_==2 || allowed_==(i+7)/2) { tcPDPtr lm = getParticleData(i); tcPDPtr lp = lm->CC(); add(new_ptr((Tree2toNDiagram(2), em, ep, 1, h0_, 3, lm, 3, lp, -1))); } } if(allowed_==0 || allowed_==12) { tcPDPtr g = getParticleData(ParticleID::g); add(new_ptr((Tree2toNDiagram(2), em, ep, 1, h0_, 3, g, 3, g, -1))); } + if(allowed_==0 || allowed_==13) { + tcPDPtr Wp = getParticleData(ParticleID::Wplus); + add(new_ptr((Tree2toNDiagram(2), em, ep, 1, h0_, 3, Wp, 3, Wp->CC(), -1))); + } + if(allowed_==0 || allowed_==14) { + tcPDPtr Z0 = getParticleData(ParticleID::Z0); + add(new_ptr((Tree2toNDiagram(2), em, ep, 1, h0_, 3, Z0, 3, Z0, -1))); + } } double MEee2Higgs2SM::me2() const { double aver=0.; // get the order right int ielectron(0),ipositron(1),ilp(2),ilm(3); if(mePartonData()[0]->id()!=11) swap(ielectron,ipositron); if(mePartonData()[2]->id()>mePartonData()[3]->id()) swap(ilp,ilm); // the arrays for the wavefunction to be passed to the matrix element vector fin; vector ain; for(unsigned int ihel=0;ihel<2;++ihel) { fin .push_back(SpinorWaveFunction(meMomenta()[ielectron], mePartonData()[ielectron],ihel,incoming)); ain .push_back(SpinorBarWaveFunction(meMomenta()[ipositron], mePartonData()[ipositron],ihel,incoming)); } // H -> f fbar - if(mePartonData()[2]->id()!=ParticleID::g) { + if(abs(mePartonData()[2]->id())<20) { vector aout; vector fout; for(unsigned int ihel=0;ihel<2;++ihel) { fout.push_back(SpinorBarWaveFunction(meMomenta()[ilm], mePartonData()[ilm],ihel,outgoing)); aout.push_back(SpinorWaveFunction(meMomenta()[ilp], mePartonData()[ilp],ihel,outgoing)); } HelicityME(fin,ain,fout,aout,aver); if(mePartonData()[ilm]->id()<=6) aver*=3.; } - else { + else if(mePartonData()[2]->id()==21) { vector g1,g2; for(unsigned int ihel=0;ihel<2;++ihel) { g1.push_back(VectorWaveFunction(meMomenta()[2],mePartonData()[2], 2*ihel,outgoing)); g2.push_back(VectorWaveFunction(meMomenta()[3],mePartonData()[3], 2*ihel,outgoing)); } ggME(fin,ain,g1,g2,aver); aver *= 8.; } + else if(abs(mePartonData()[2]->id())==24 ||mePartonData()[2]->id()==23) { + vector g1,g2; + for(unsigned int ihel=0;ihel<3;++ihel) { + g1.push_back(VectorWaveFunction(meMomenta()[2],mePartonData()[2], + ihel,outgoing)); + g2.push_back(VectorWaveFunction(meMomenta()[3],mePartonData()[3], + ihel,outgoing)); + } + WWME(fin,ain,g1,g2,aver); + } + else { + assert(false); + } return aver; } // the helicity amplitude matrix element ProductionMatrixElement MEee2Higgs2SM::HelicityME(vector fin, vector ain, vector fout, vector aout, double & aver) const { // the particles should be in the order // for the incoming // 0 incoming fermion (u spinor) // 1 incoming antifermion (vbar spinor) // for the outgoing // 0 outgoing fermion (ubar spinor) // 1 outgoing antifermion (v spinor) // me to be returned ProductionMatrixElement output(PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1Half,PDT::Spin1Half); // wavefunctions for the intermediate particles ScalarWaveFunction interh; // temporary storage of the different diagrams Complex diag; aver=0.; // sum over helicities to get the matrix element unsigned int inhel1,inhel2,outhel1,outhel2; for(inhel1=0;inhel1<2;++inhel1) { for(inhel2=0;inhel2<2;++inhel2) { interh = FFHVertex_->evaluate(sHat(),1,h0_,fin[inhel1],ain[inhel2]); for(outhel1=0;outhel1<2;++outhel1) { for(outhel2=0;outhel2<2;++outhel2) { diag = FFHVertex_->evaluate(sHat(),aout[outhel2], fout[outhel1],interh); output(inhel1,inhel2,outhel1,outhel2)=diag; aver +=real(diag*conj(diag)); } } } } return output; } Selector MEee2Higgs2SM::diagrams(const DiagramVector &) const { Selector sel;sel.insert(1.0, 0); return sel; } Selector MEee2Higgs2SM::colourGeometries(tcDiagPtr diag) const { static ColourLines neutral ( " " ); static ColourLines quarks ( "-5 4"); static ColourLines gluons ( "-5 4, 5 -4"); Selector sel; int id = abs((diag->partons()[2])->id()); if (id<=6 ) sel.insert(1.0, &quarks); else if (id==21) sel.insert(1.0, &gluons); else sel.insert(1.0, &neutral); return sel; } void MEee2Higgs2SM::persistentOutput(PersistentOStream & os) const { os << FFHVertex_ << HGGVertex_ << h0_ << allowed_; } void MEee2Higgs2SM::persistentInput(PersistentIStream & is, int) { is >> FFHVertex_ >> HGGVertex_ >> h0_ >> allowed_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigMEee2Higgs2SM("Herwig::MEee2Higgs2SM", "HwMELepton.so"); void MEee2Higgs2SM::Init() { static ClassDocumentation documentation ("The MEee2Higgs2SM class implements the matrix element for e+e- to" " SM particle via Higgs exchnage and is designed to be a process to test" " things involving scalar particles. "); static Switch interfaceallowed ("Allowed", "Allowed outgoing particles", &MEee2Higgs2SM::allowed_, 0, false, false); static SwitchOption interfaceallowedAll (interfaceallowed, "All", "All SM particles allowed", 0); static SwitchOption interfaceallowed1 (interfaceallowed, "Quarks", "Only the quarks allowed", 1); static SwitchOption interfaceallowed2 (interfaceallowed, "Leptons", "Only the leptons allowed", 2); static SwitchOption interfacealloweddown (interfaceallowed, "Down", "Only d dbar allowed", 3); static SwitchOption interfaceallowedup (interfaceallowed, "Up", "Only u ubar allowed", 4); static SwitchOption interfaceallowedstrange (interfaceallowed, "Strange", "Only s sbar allowed", 5); static SwitchOption interfaceallowedcharm (interfaceallowed, "Charm", "Only c cbar allowed", 6); static SwitchOption interfaceallowedbottom (interfaceallowed, "Bottom", "Only b bbar allowed", 7); static SwitchOption interfaceallowedtop (interfaceallowed, "Top", "Only t tbar allowed", 8); static SwitchOption interfaceallowedelectron (interfaceallowed, "Electron", "Only e+e- allowed", 9); static SwitchOption interfaceallowedMuon (interfaceallowed, "Muon", "Only mu+mu- allowed", 10); static SwitchOption interfaceallowedTau (interfaceallowed, "Tau", "Only tau+tau- allowed", 11); static SwitchOption interfaceallowedGluon (interfaceallowed, "Gluon", "Only gg allowed", 12); + static SwitchOption interfaceallowedW + (interfaceallowed, + "W", + "Only W+W- allowed", + 13); + static SwitchOption interfaceallowedZ + (interfaceallowed, + "Z", + "Only ZZ allowed", + 14); } void MEee2Higgs2SM::constructVertex(tSubProPtr sub) { // extract the particles in the hard process ParticleVector hard; hard.push_back(sub->incoming().first);hard.push_back(sub->incoming().second); hard.push_back(sub->outgoing()[0]);hard.push_back(sub->outgoing()[1]); if(hard[0]->id()id()) swap(hard[0],hard[1]); if(hard[2]->id()id()) swap(hard[2],hard[3]); vector fin; vector ain; SpinorWaveFunction( fin ,hard[0],incoming,false,true); SpinorBarWaveFunction(ain ,hard[1],incoming,false,true); double me; ProductionMatrixElement prodme; - if(hard[2]->id()!=ParticleID::g) { + if(abs(hard[2]->id())<20) { vector aout; vector fout; SpinorBarWaveFunction(fout,hard[2],outgoing,true ,true); SpinorWaveFunction( aout,hard[3],outgoing,true ,true); // calculate the matrix element prodme=HelicityME(fin,ain,fout,aout,me); } - else { + else if(abs(hard[2]->id())==21) { vector g1,g2; VectorWaveFunction(g1,hard[2],outgoing,true,true); VectorWaveFunction(g2,hard[3],outgoing,true,true); g1[1]=g1[2]; g2[1]=g2[2]; prodme=ggME(fin,ain,g1,g2,me); } + else { + vector g1,g2; + VectorWaveFunction(g1,hard[2],outgoing,true,false); + VectorWaveFunction(g2,hard[3],outgoing,true,false); + prodme=WWME(fin,ain,g1,g2,me); + } // construct the vertex HardVertexPtr hardvertex=new_ptr(HardVertex()); // set the matrix element for the vertex hardvertex->ME(prodme); // set the pointers and to and from the vertex for(unsigned int ix=0;ix<4;++ix) { (hard[ix]->spinInfo())->productionVertex(hardvertex); } } void MEee2Higgs2SM::rebind(const TranslationMap & trans) { // dummy = trans.translate(dummy); FFHVertex_ = trans.translate(FFHVertex_); HGGVertex_ = trans.translate(HGGVertex_); h0_ = trans.translate(h0_); ME2to2Base::rebind(trans); } IVector MEee2Higgs2SM::getReferences() { IVector ret = ME2to2Base::getReferences(); ret.push_back(FFHVertex_); ret.push_back(HGGVertex_); ret.push_back(h0_); return ret; } // the helicity amplitude matrix element ProductionMatrixElement MEee2Higgs2SM::ggME(vector fin, vector ain, vector g1, vector g2, double & aver) const { ProductionMatrixElement output(PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1,PDT::Spin1); // wavefunctions for the intermediate particles ScalarWaveFunction interh; // temporary storage of the different diagrams Complex diag; aver=0.; // sum over helicities to get the matrix element unsigned int inhel1,inhel2,outhel1,outhel2; for(inhel1=0;inhel1<2;++inhel1) { for(inhel2=0;inhel2<2;++inhel2) { interh = FFHVertex_->evaluate(sHat(),1,h0_,fin[inhel1],ain[inhel2]); for(outhel1=0;outhel1<2;++outhel1) { for(outhel2=0;outhel2<2;++outhel2) { diag = HGGVertex_->evaluate(sHat(),g1[outhel1],g2[outhel2],interh); output(inhel1,inhel2,2*outhel1,2*outhel2)=diag; - aver +=real(diag*conj(diag)); + aver += norm(diag); } } } } return output; } + +// the helicity amplitude matrix element +ProductionMatrixElement MEee2Higgs2SM::WWME(vector fin, + vector ain, + vector g1, + vector g2, + double & aver) const { + ProductionMatrixElement output(PDT::Spin1Half,PDT::Spin1Half, + PDT::Spin1,PDT::Spin1); + // wavefunctions for the intermediate particles + ScalarWaveFunction interh; + // temporary storage of the different diagrams + Complex diag; + aver=0.; + // sum over helicities to get the matrix element + unsigned int inhel1,inhel2,outhel1,outhel2; + for(inhel1=0;inhel1<2;++inhel1) { + for(inhel2=0;inhel2<2;++inhel2) { + interh = FFHVertex_->evaluate(sHat(),1,h0_,fin[inhel1],ain[inhel2]); + for(outhel1=0;outhel1<3;++outhel1) { + for(outhel2=0;outhel2<3;++outhel2) { + diag = HWWVertex_->evaluate(sHat(),g1[outhel1],g2[outhel2],interh); + output(inhel1,inhel2,outhel1,outhel2)=diag; + aver += norm(diag); + } + } + } + } + return output; +} diff --git a/MatrixElement/Lepton/MEee2Higgs2SM.h b/MatrixElement/Lepton/MEee2Higgs2SM.h --- a/MatrixElement/Lepton/MEee2Higgs2SM.h +++ b/MatrixElement/Lepton/MEee2Higgs2SM.h @@ -1,236 +1,254 @@ // -*- C++ -*- // // MEee2Higgs2SM.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_MEee2Higgs2SM_H #define HERWIG_MEee2Higgs2SM_H // // This is the declaration of the MEee2Higgs2SM class. // #include "ThePEG/MatrixElement/ME2to2Base.h" #include "Herwig/Models/StandardModel/StandardModel.h" #include "Herwig/MatrixElement/ProductionMatrixElement.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" namespace Herwig { using namespace ThePEG; /** * The MEee2Higgs2SM class implements the production of an \f$s\f$-channel * Higgs in \f$e^+e^-\f$ collisions in order to allow easy tests of Higgs * decays. It should not be used for physics studies. * * @see \ref MEee2Higgs2SMInterfaces "The interfaces" * defined for MEee2Higgs2SM. */ class MEee2Higgs2SM: public ME2to2Base { public: /** * The default constructor. */ inline MEee2Higgs2SM() : allowed_(0) {} /** @name Virtual functions required by the MEBase class. */ //@{ /** * Return the order in \f$\alpha_S\f$ in which this matrix * element is given. */ virtual unsigned int orderInAlphaS() const; /** * Return the order in \f$\alpha_{EW}\f$ in which this matrix * element is given. */ virtual unsigned int orderInAlphaEW() const; /** * The matrix element for the kinematical configuration * previously provided by the last call to setKinematics(), suitably * scaled by sHat() to give a dimension-less number. * @return the matrix element scaled with sHat() to give a * dimensionless number. */ virtual double me2() const; /** * Return the scale associated with the last set phase space point. */ virtual Energy2 scale() const; /** * Add all possible diagrams with the add() function. */ virtual void getDiagrams() const; /** * Get diagram selector. With the information previously supplied with the * setKinematics method, a derived class may optionally * override this method to weight the given diagrams with their * (although certainly not physical) relative probabilities. * @param dv the diagrams to be weighted. * @return a Selector relating the given diagrams to their weights. */ virtual Selector diagrams(const DiagramVector & dv) const; /** * Return a Selector with possible colour geometries for the selected * diagram weighted by their relative probabilities. * @param diag the diagram chosen. * @return the possible colour geometries weighted by their * relative probabilities. */ virtual Selector colourGeometries(tcDiagPtr diag) const; /** * set up the spin correlations */ virtual void constructVertex(tSubProPtr sub); //@} 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. */ inline virtual IBPtr clone() const {return new_ptr(*this);} /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ inline virtual IBPtr fullclone() const {return new_ptr(*this);} //@} 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(); /** * Rebind pointer to other Interfaced objects. Called in the setup phase * after all objects used in an EventGenerator has been cloned so that * the pointers will refer to the cloned objects afterwards. * @param trans a TranslationMap relating the original objects to * their respective clones. * @throws RebindException if no cloned object was found for a given * pointer. */ virtual void rebind(const TranslationMap & trans) ; /** * Return a vector of all pointers to Interfaced objects used in this * object. * @return a vector of pointers. */ virtual IVector getReferences(); //@} private: /** * The matrix element * @param fin The incoming spinor wavefunction * @param ain The incoming spinorbar wavefunction * @param fout The outgoing spinor bar wavefunction * @param aout The outgoing spinor wavefunction * @param me The spin averaged matrix element */ ProductionMatrixElement HelicityME(vector fin, vector ain, vector fout, vector aout,double& me) const; /** * \f$H\to gg\f$ matrix element * @param fin The incoming spinor wavefunction * @param ain The incoming spinorbar wavefunction * @param g1 Outgoing gluon wavefunction * @param g2 Outgoing gluon wavefunction * @param me The spin averaged matrix element */ ProductionMatrixElement ggME(vector fin, vector ain, vector g1, vector g2, double & me) const; + /** + * \f$H\to WW\f$ matrix element + * @param fin The incoming spinor wavefunction + * @param ain The incoming spinorbar wavefunction + * @param g1 Outgoing W wavefunction + * @param g2 Outgoing W wavefunction + * @param me The spin averaged matrix element + */ + ProductionMatrixElement WWME(vector fin, + vector ain, + vector g1, + vector g2, + double & me) const; private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MEee2Higgs2SM & operator=(const MEee2Higgs2SM &) = delete; private: /** * Pointer to the Higgs fermion-antifermion vertex */ AbstractFFSVertexPtr FFHVertex_; /** * Pointer to Higgs-gluon-gluon vertex */ AbstractVVSVertexPtr HGGVertex_; /** + * Pointer to Higgs-WW vertex + */ + AbstractVVSVertexPtr HWWVertex_; + + /** * Allowed outgoing particles */ int allowed_; /** * Pointer to the Higgs ParticleData object */ PDPtr h0_; }; } #endif /* HERWIG_MEee2Higgs2SM_H */ 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,243 +1,249 @@ // -*- 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_ << _theSM; } void HalfHalfZeroEWSplitFn::persistentInput(PersistentIStream & is, int) { is >> ghqq_ >> _theSM; } // 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,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,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 1. - sqrt(1. - 2.*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.); Energy2 tC = t/(z*(1-z)); getCouplings(gH,ids,tC); double mqt = mq/sqrt(tC); double mHt = mH/sqrt(tC); 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; + cerr << "testing kernal"; + for(unsigned int ix=0;ix<2;++ix) { + for(unsigned int iy=0;iy<2;++iy) + cerr << (*kernal)(ix,iy,0) << " "; + cerr << "\n"; + } // return the answer return kernal; }