diff --git a/Decay/Perturbative/SMTopDecayer.cc b/Decay/Perturbative/SMTopDecayer.cc --- a/Decay/Perturbative/SMTopDecayer.cc +++ b/Decay/Perturbative/SMTopDecayer.cc @@ -1,1141 +1,1142 @@ // -*- C++ -*- // // SMTopDecayer.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 SMTopDecayer class. // #include "SMTopDecayer.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/DecayMode.h" #include "Herwig/Decay/DecayVertex.h" #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" #include "Herwig/PDT/ThreeBodyAllOn1IntegralCalculator.h" #include "Herwig/Shower/RealEmissionProcess.h" #include "Herwig/Shower/Core/Base/ShowerProgenitor.h" #include "Herwig/Shower/Core/Base/ShowerParticle.h" #include "Herwig/Shower/Core/Base/Branching.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" using namespace Herwig; using namespace ThePEG::Helicity; SMTopDecayer::SMTopDecayer() : _wquarkwgt(6,0.),_wleptonwgt(3,0.), _xg_sampling(1.5), _initialenhance(1.), _finalenhance(2.3), _useMEforT2(true) { _wleptonwgt[0] = 0.302583; _wleptonwgt[1] = 0.301024; _wleptonwgt[2] = 0.299548; _wquarkwgt[0] = 0.851719; _wquarkwgt[1] = 0.0450162; _wquarkwgt[2] = 0.0456962; _wquarkwgt[3] = 0.859839; _wquarkwgt[4] = 3.9704e-06; _wquarkwgt[5] = 0.000489657; generateIntermediates(true); } bool SMTopDecayer::accept(tcPDPtr parent, const tPDVector & children) const { if(abs(parent->id()) != ParticleID::t) return false; int id0(0),id1(0),id2(0); for(tPDVector::const_iterator it = children.begin(); it != children.end();++it) { int id=(**it).id(),absid(abs(id)); if(absid==ParticleID::b&&double(id)/double(parent->id())>0) { id0=id; } else { switch (absid) { case ParticleID::nu_e: case ParticleID::nu_mu: case ParticleID::nu_tau: id1 = id; break; case ParticleID::eminus: case ParticleID::muminus: case ParticleID::tauminus: id2 = id; break; case ParticleID::b: case ParticleID::d: case ParticleID::s: id1 = id; break; case ParticleID::u: case ParticleID::c: id2=id; break; default : break; } } } if(id0==0||id1==0||id2==0) return false; if(double(id1)/double(id2)>0) return false; return true; } ParticleVector SMTopDecayer::decay(const Particle & parent, const tPDVector & children) const { int id1(0),id2(0); for(tPDVector::const_iterator it = children.begin(); it != children.end();++it) { int id=(**it).id(),absid=abs(id); if(absid == ParticleID::b && double(id)/double(parent.id())>0) continue; //leptons if(absid > 10 && absid%2==0) id1=absid; if(absid > 10 && absid%2==1) id2=absid; //quarks if(absid < 10 && absid%2==0) id2=absid; if(absid < 10 && absid%2==1) id1=absid; } unsigned int imode(0); if(id2 >=11 && id2<=16) imode = (id1-12)/2; else imode = id1+1+id2/2; bool cc = parent.id() == ParticleID::tbar; ParticleVector out(generate(true,cc,imode,parent)); //arrange colour flow PPtr pparent=const_ptr_cast(&parent); out[1]->incomingColour(pparent,out[1]->id()<0); ParticleVector products = out[0]->children(); if(products[0]->hasColour()) products[0]->colourNeighbour(products[1],true); else if(products[0]->hasAntiColour()) products[0]->colourNeighbour(products[1],false); return out; } void SMTopDecayer::persistentOutput(PersistentOStream & os) const { os << _wvertex << _wquarkwgt << _wleptonwgt << _wplus << _alpha << _initialenhance << _finalenhance << _xg_sampling << _useMEforT2; } void SMTopDecayer::persistentInput(PersistentIStream & is, int) { is >> _wvertex >> _wquarkwgt >> _wleptonwgt >> _wplus >> _alpha >> _initialenhance >> _finalenhance >> _xg_sampling >> _useMEforT2; } ClassDescription SMTopDecayer::initSMTopDecayer; // Definition of the static class description member. void SMTopDecayer::Init() { static ClassDocumentation documentation ("This is the implementation of the SMTopDecayer which " "decays top quarks into bottom quarks and either leptons " "or quark-antiquark pairs including the matrix element for top decay", "The matrix element correction for top decay \\cite{Hamilton:2006ms}.", "%\\cite{Hamilton:2006ms}\n" "\\bibitem{Hamilton:2006ms}\n" " K.~Hamilton and P.~Richardson,\n" " ``A simulation of QCD radiation in top quark decays,''\n" " JHEP {\\bf 0702}, 069 (2007)\n" " [arXiv:hep-ph/0612236].\n" " %%CITATION = JHEPA,0702,069;%%\n"); static ParVector interfaceQuarkWeights ("QuarkWeights", "Maximum weights for the hadronic decays", &SMTopDecayer::_wquarkwgt, 6, 1.0, 0.0, 10.0, false, false, Interface::limited); static ParVector interfaceLeptonWeights ("LeptonWeights", "Maximum weights for the semi-leptonic decays", &SMTopDecayer::_wleptonwgt, 3, 1.0, 0.0, 10.0, false, false, Interface::limited); static Parameter interfaceEnhancementFactor ("InitialEnhancementFactor", "The enhancement factor for initial-state radiation in the shower to ensure" " the weight for the matrix element correction is less than one.", &SMTopDecayer::_initialenhance, 1.0, 1.0, 10000.0, false, false, Interface::limited); static Parameter interfaceFinalEnhancementFactor ("FinalEnhancementFactor", "The enhancement factor for final-state radiation in the shower to ensure" " the weight for the matrix element correction is less than one", &SMTopDecayer::_finalenhance, 1.6, 1.0, 1000.0, false, false, Interface::limited); static Parameter interfaceSamplingTopHardMEC ("SamplingTopHardMEC", "The importance sampling power for choosing an initial xg, " "to sample xg according to xg^-_xg_sampling", &SMTopDecayer::_xg_sampling, 1.5, 1.2, 2.0, false, false, Interface::limited); static Switch interfaceUseMEForT2 ("UseMEForT2", "Use the matrix element correction, if available to fill the T2" " region for the decay shower and don't fill using the shower", &SMTopDecayer::_useMEforT2, true, false, false); static SwitchOption interfaceUseMEForT2Shower (interfaceUseMEForT2, "Shower", "Use the shower to fill the T2 region", false); static SwitchOption interfaceUseMEForT2ME (interfaceUseMEForT2, "ME", "Use the Matrix element to fill the T2 region", true); static Reference interfaceCoupling ("Coupling", "Pointer to the object to calculate the coupling for the correction", &SMTopDecayer::_alpha, false, false, true, false, false); } double SMTopDecayer::me2(const int, const Particle & inpart, const ParticleVector & decay, MEOption meopt) const { if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1Half,PDT::Spin1Half))); // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(inpart.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&inpart), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&inpart), incoming); } // setup spin info when needed if(meopt==Terminate) { // for the decaying particle if(inpart.id()>0) { SpinorWaveFunction:: constructSpinInfo(_inHalf,const_ptr_cast(&inpart),incoming,true); SpinorBarWaveFunction::constructSpinInfo(_inHalfBar,decay[0],outgoing,true); SpinorWaveFunction ::constructSpinInfo(_outHalf ,decay[1],outgoing,true); SpinorBarWaveFunction::constructSpinInfo(_outHalfBar,decay[2],outgoing,true); } else { SpinorBarWaveFunction:: constructSpinInfo(_inHalfBar,const_ptr_cast(&inpart),incoming,true); SpinorWaveFunction::constructSpinInfo(_inHalf,decay[0],outgoing,true); SpinorBarWaveFunction::constructSpinInfo(_outHalfBar,decay[1],outgoing,true); SpinorWaveFunction ::constructSpinInfo(_outHalf ,decay[2],outgoing,true); } } if ( ( decay[1]->momentum() + decay[2]->momentum() ).m() < decay[1]->data().constituentMass() + decay[2]->data().constituentMass() ) return 0.0; // spinors for the decay product if(inpart.id()>0) { SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar ,decay[0],outgoing); SpinorWaveFunction ::calculateWaveFunctions(_outHalf ,decay[1],outgoing); SpinorBarWaveFunction::calculateWaveFunctions(_outHalfBar,decay[2],outgoing); } else { SpinorWaveFunction ::calculateWaveFunctions(_inHalf ,decay[0],outgoing); SpinorBarWaveFunction::calculateWaveFunctions(_outHalfBar,decay[1],outgoing); SpinorWaveFunction ::calculateWaveFunctions(_outHalf ,decay[2],outgoing); } Energy2 scale(sqr(inpart.mass())); if(inpart.id() == ParticleID::t) { //Define intermediate vector wave-function for Wplus tcPDPtr Wplus(getParticleData(ParticleID::Wplus)); VectorWaveFunction inter; unsigned int thel,bhel,fhel,afhel; for(thel = 0;thel<2;++thel){ for(bhel = 0;bhel<2;++bhel){ inter = _wvertex->evaluate(scale,1,Wplus,_inHalf[thel], _inHalfBar[bhel]); for(afhel=0;afhel<2;++afhel){ for(fhel=0;fhel<2;++fhel){ (*ME())(thel,bhel,afhel,fhel) = _wvertex->evaluate(scale,_outHalf[afhel], _outHalfBar[fhel],inter); } } } } } else if(inpart.id() == ParticleID::tbar) { VectorWaveFunction inter; tcPDPtr Wminus(getParticleData(ParticleID::Wminus)); unsigned int tbhel,bbhel,afhel,fhel; for(tbhel = 0;tbhel<2;++tbhel){ for(bbhel = 0;bbhel<2;++bbhel){ inter = _wvertex-> evaluate(scale,1,Wminus,_inHalf[bbhel],_inHalfBar[tbhel]); for(afhel=0;afhel<2;++afhel){ for(fhel=0;fhel<2;++fhel){ (*ME())(tbhel,bbhel,fhel,afhel) = _wvertex->evaluate(scale,_outHalf[afhel], _outHalfBar[fhel],inter); } } } } } double output = (ME()->contract(_rho)).real(); if(abs(decay[1]->id())<=6) output *=3.; return output; } void SMTopDecayer::doinit() { DecayIntegrator::doinit(); //get vertices from SM object tcHwSMPtr hwsm = dynamic_ptr_cast(standardModel()); if(!hwsm) throw InitException() << "Must have Herwig::StandardModel in " << "SMTopDecayer::doinit()"; _wvertex = hwsm->vertexFFW(); //initialise _wvertex->init(); //set up decay modes _wplus = getParticleData(ParticleID::Wplus); DecayPhaseSpaceModePtr mode; DecayPhaseSpaceChannelPtr Wchannel; tPDVector extpart(4); vector wgt(1,1.0); extpart[0] = getParticleData(ParticleID::t); extpart[1] = getParticleData(ParticleID::b); //lepton modes for(int i=11; i<17;i+=2) { extpart[2] = getParticleData(-i); extpart[3] = getParticleData(i+1); mode = new_ptr(DecayPhaseSpaceMode(extpart,this)); Wchannel = new_ptr(DecayPhaseSpaceChannel(mode)); Wchannel->addIntermediate(extpart[0],0,0.0,-1,1); Wchannel->addIntermediate(_wplus,0,0.0,2,3); Wchannel->init(); mode->addChannel(Wchannel); addMode(mode,_wleptonwgt[(i-11)/2],wgt); } //quark modes unsigned int iz=0; for(int ix=1;ix<6;ix+=2) { for(int iy=2;iy<6;iy+=2) { // check that the combination of particles is allowed if(_wvertex->allowed(-ix,iy,ParticleID::Wminus)) { extpart[2] = getParticleData(-ix); extpart[3] = getParticleData( iy); mode = new_ptr(DecayPhaseSpaceMode(extpart,this)); Wchannel = new_ptr(DecayPhaseSpaceChannel(mode)); Wchannel->addIntermediate(extpart[0],0,0.0,-1,1); Wchannel->addIntermediate(_wplus,0,0.0,2,3); Wchannel->init(); mode->addChannel(Wchannel); addMode(mode,_wquarkwgt[iz],wgt); ++iz; } else { throw InitException() << "SMTopDecayer::doinit() the W vertex" << "cannot handle all the quark modes" << Exception::abortnow; } } } } void SMTopDecayer::dataBaseOutput(ofstream & os,bool header) const { if(header) os << "update decayers set parameters=\""; // parameters for the DecayIntegrator base class for(unsigned int ix=0;ix<_wquarkwgt.size();++ix) { os << "newdef " << name() << ":QuarkWeights " << ix << " " << _wquarkwgt[ix] << "\n"; } for(unsigned int ix=0;ix<_wleptonwgt.size();++ix) { os << "newdef " << name() << ":LeptonWeights " << ix << " " << _wleptonwgt[ix] << "\n"; } DecayIntegrator::dataBaseOutput(os,false); if(header) os << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } void SMTopDecayer::doinitrun() { DecayIntegrator::doinitrun(); if(initialize()) { for(unsigned int ix=0;ixmaxWeight(); else _wquarkwgt [ix-3] = mode(ix)->maxWeight(); } } } WidthCalculatorBasePtr SMTopDecayer::threeBodyMEIntegrator(const DecayMode & dm) const { // identify W decay products int sign = dm.parent()->id() > 0 ? 1 : -1; int iferm(0),ianti(0); for(ParticleMSet::const_iterator pit=dm.products().begin(); pit!=dm.products().end();++pit) { int id = (**pit).id(); if(id*sign != ParticleID::b) { if (id*sign > 0 ) iferm = id*sign; else ianti = id*sign; } } assert(iferm!=0&&ianti!=0); // work out which mode we are doing int imode(-1); for(unsigned int ix=0;ixexternalParticles(2)->id() == ianti && mode(ix)->externalParticles(3)->id() == iferm ) { imode = ix; break; } } assert(imode>=0); // get the masses we need Energy m[3] = {mode(imode)->externalParticles(1)->mass(), mode(imode)->externalParticles(3)->mass(), mode(imode)->externalParticles(2)->mass()}; return new_ptr(ThreeBodyAllOn1IntegralCalculator (3,_wplus->mass(),_wplus->width(),0.0,*this,imode,m[0],m[1],m[2])); } InvEnergy SMTopDecayer::threeBodydGammads(const int imode, const Energy2 mt2, const Energy2 mffb2, const Energy mb, const Energy mf, const Energy mfb) const { Energy mffb(sqrt(mffb2)); Energy mw(_wplus->mass()); Energy2 mw2(sqr(mw)),gw2(sqr(_wplus->width())); Energy mt(sqrt(mt2)); Energy Eb = 0.5*(mt2-mffb2-sqr(mb))/mffb; Energy Ef = 0.5*(mffb2-sqr(mfb)+sqr(mf))/mffb; Energy Ebm = sqrt(sqr(Eb)-sqr(mb)); Energy Efm = sqrt(sqr(Ef)-sqr(mf)); Energy2 upp = sqr(Eb+Ef)-sqr(Ebm-Efm); Energy2 low = sqr(Eb+Ef)-sqr(Ebm+Efm); InvEnergy width=(dGammaIntegrand(mffb2,upp,mt,mb,mf,mfb,mw)- dGammaIntegrand(mffb2,low,mt,mb,mf,mfb,mw)) /32./mt2/mt/8/pow(Constants::pi,3)/(sqr(mffb2-mw2)+mw2*gw2); // couplings width *= 0.25*sqr(4.*Constants::pi*generator()->standardModel()->alphaEM(mt2)/ generator()->standardModel()->sin2ThetaW()); width *= generator()->standardModel()->CKM(*mode(imode)->externalParticles(0), *mode(imode)->externalParticles(1)); if(abs(mode(imode)->externalParticles(2)->id())<=6) { width *=3.; if(abs(mode(imode)->externalParticles(2)->id())%2==0) width *=generator()->standardModel()->CKM(*mode(imode)->externalParticles(2), *mode(imode)->externalParticles(3)); else width *=generator()->standardModel()->CKM(*mode(imode)->externalParticles(3), *mode(imode)->externalParticles(2)); } // final spin average assert(!std::isnan(double(width*MeV))); return 0.5*width; } Energy6 SMTopDecayer::dGammaIntegrand(Energy2 mffb2, Energy2 mbf2, Energy mt, Energy mb, Energy mf, Energy mfb, Energy mw) const { Energy2 mt2(sqr(mt)) ,mb2(sqr(mb)) ,mf2(sqr(mf )),mfb2(sqr(mfb )),mw2(sqr(mw )); Energy4 mt4(sqr(mt2)),mb4(sqr(mb2)),mf4(sqr(mf2)),mfb4(sqr(mfb2)),mw4(sqr(mw2)); return -mbf2 * ( + 6 * mb2 * mf2 * mfb2 * mffb2 + 6 * mb2 * mt2 * mfb2 * mffb2 + 6 * mb2 * mt2 * mf2 * mffb2 + 12 * mb2 * mt2 * mf2 * mfb2 - 3 * mb2 * mfb4 * mffb2 + 3 * mb2 * mf2 * mffb2 * mffb2 - 3 * mb2 * mf4 * mffb2 - 6 * mb2 * mt2 * mfb4 - 6 * mb2 * mt2 * mf4 - 3 * mb4 * mfb2 * mffb2 - 3 * mb4 * mf2 * mffb2 - 6 * mb4 * mf2 * mfb2 + 3 * mt4 * mf4 + 3 * mb4 * mfb4 + 3 * mb4 * mf4 + 3 * mt4 * mfb4 + 3 * mb2 * mfb2 * mffb2 * mffb2 + 3 * mt2 * mfb2 * mffb2 * mffb2 - 3 * mt2 * mfb4 * mffb2 + 3 * mt2 * mf2 * mffb2 * mffb2 - 3 * mt2 * mf4 * mffb2 - 3 * mt4 * mfb2 * mffb2 - 3 * mt4 * mf2 * mffb2 - 6 * mt4 * mf2 * mfb2 + 6 * mt2 * mf2 * mfb2 * mffb2 + 12 * mt2 * mf2 * mw4 + 12 * mb2 * mfb2 * mw4 + 12 * mb2 * mt2 * mw4 + 6 * mw2 * mt2 * mfb2 * mbf2 - 12 * mw2 * mt2 * mf2 * mffb2 - 6 * mw2 * mt2 * mf2 * mbf2 - 12 * mw2 * mt2 * mf2 * mfb2 - 12 * mw2 * mb2 * mfb2 * mffb2 - 6 * mw2 * mb2 * mfb2 * mbf2 + 6 * mw2 * mb2 * mf2 * mbf2 - 12 * mw2 * mb2 * mf2 * mfb2 - 12 * mw2 * mb2 * mt2 * mfb2 - 12 * mw2 * mb2 * mt2 * mf2 + 12 * mf2 * mfb2 * mw4 + 4 * mbf2 * mbf2 * mw4 - 6 * mfb2 * mbf2 * mw4 - 6 * mf2 * mbf2 * mw4 - 6 * mt2 * mbf2 * mw4 - 6 * mb2 * mbf2 * mw4 + 12 * mw2 * mt2 * mf4 + 12 * mw2 * mt4 * mf2 + 12 * mw2 * mb2 * mfb4 + 12 * mw2 * mb4 * mfb2) /mw4 / 3.; } void SMTopDecayer::initializeMECorrection(RealEmissionProcessPtr born, double & initial, double & final) { + if(born->bornOutgoing().size()!=2) return; // check the outgoing particles PPtr part[2]; for(unsigned int ix=0;ixbornOutgoing().size();++ix) { part[ix]= born->bornOutgoing()[ix]; } // check the final-state particles and get the masses if(abs(part[0]->id())==ParticleID::Wplus&&abs(part[1]->id())==ParticleID::b) { _ma=part[0]->mass(); _mc=part[1]->mass(); } else if(abs(part[1]->id())==ParticleID::Wplus&&abs(part[0]->id())==ParticleID::b) { _ma=part[1]->mass(); _mc=part[0]->mass(); } else { return; } // set the top mass _mt=born->bornIncoming()[0]->mass(); // set the gluon mass _mg=getParticleData(ParticleID::g)->constituentMass(); // set the radiation enhancement factors initial = _initialenhance; final = _finalenhance; // reduced mass parameters _a=sqr(_ma/_mt); _g=sqr(_mg/_mt); _c=sqr(_mc/_mt); double lambda = sqrt(1.+sqr(_a)+sqr(_c)-2.*_a-2.*_c-2.*_a*_c); _ktb = 0.5*(3.-_a+_c+lambda); _ktc = 0.5*(1.-_a+3.*_c+lambda); useMe(); } RealEmissionProcessPtr SMTopDecayer::applyHardMatrixElementCorrection(RealEmissionProcessPtr born) { // Get b and a and put them in particle vector ba in that order... ParticleVector ba; for(unsigned int ix=0;ixbornOutgoing().size();++ix) ba.push_back(born->bornOutgoing()[ix]); if(abs(ba[0]->id())!=5) swap(ba[0],ba[1]); assert(born->bornIncoming().size()==1); // Now decide if we get an emission into the dead region. // If there is an emission newfs stores momenta for a,c,g // according to NLO decay matrix element. vector newfs = applyHard(ba,_ktb,_ktc); // If there was no gluon emitted return. if(newfs.size()!=3) return RealEmissionProcessPtr(); // Sanity checks to ensure energy greater than mass etc :) bool check = true; tcPDPtr gluondata=getParticleData(ParticleID::g); if (newfs[0].e()data().constituentMass()) check = false; if (newfs[1].e()mass()) check = false; if (newfs[2].e()constituentMass()) check = false; // Return if insane: if (!check) return RealEmissionProcessPtr(); // // Set masses in 5-vectors: newfs[0].setMass(ba[0]->mass()); newfs[1].setMass(ba[1]->mass()); newfs[2].setMass(ZERO); // The next part of this routine sets the colour structure. // To do this for decays we assume that the gluon comes from c! // First create new particle objects for c, a and gluon: PPtr newg = gluondata->produceParticle(newfs[2]); PPtr newc = ba[0]->data().produceParticle(newfs[0]); PPtr newa = ba[1]->data().produceParticle(newfs[1]); born->spectator(0); born->emitted(3); // decaying particle born->incoming().push_back(born->bornIncoming()[0]->dataPtr()-> produceParticle(born->bornIncoming()[0]->momentum())); // colour flow newg->incomingColour(born->incoming()[0],ba[0]->id()<0); newg->colourConnect(newc ,ba[0]->id()<0); if(born->bornOutgoing()[0]->id()==newc->id()) { born->outgoing().push_back(newc); born->outgoing().push_back(newa); born->emitter(1); } else { born->outgoing().push_back(newa); born->outgoing().push_back(newc); born->emitter(2); } born->outgoing().push_back(newg); // boost for the W LorentzRotation trans(ba[1]->momentum().findBoostToCM()); trans.boost(newfs[1].boostVector()); born->transformation(trans); if(!inTheDeadRegion(_xg,_xa,_ktb,_ktc)) { generator()->log() << "SMTopDecayer::applyHardMatrixElementCorrection()\n" << "Just found a point that escaped from the dead region!\n" << " _xg: " << _xg << " _xa: " << _xa << " newfs.size(): " << newfs.size() << endl; } born->interaction(ShowerInteraction::QCD); return born; } vector SMTopDecayer:: applyHard(const ParticleVector &p,double ktb, double ktc) { // ********************************* // // First we see if we get a dead // // region event: _xa,_xg // // ********************************* // vector fs; // Return if there is no (NLO) gluon emission: double weight = getHard(ktb,ktc); if(weight>1.) { generator()->log() << "Weight greater than 1 for hard emission in " << "SMTopDecayer::applyHard xg = " << _xg << " xa = " << _xa << "\n"; weight=1.; } // Accept/Reject if (weightmomentum() + p[1]->momentum(); // Define/assign momenta of c,a and the gluon: if(abs(p[0]->id())==5) { pc_lab = p[0]->momentum(); pa_lab = p[1]->momentum(); } else { pc_lab = p[1]->momentum(); pa_lab = p[0]->momentum(); } // Calculate the boost to the b rest frame: SpinOneLorentzRotation rot0(pb_lab.findBoostToCM()); // Calculate the rotation matrix to position a along the +z direction // in the rest frame of b and does a random rotation about z: SpinOneLorentzRotation rot1 = rotateToZ(rot0*pa_lab); // Calculate the boost from the b rest frame back to the lab: // and the inverse of the random rotation about the z-axis and the // rotation required to align a with +z: SpinOneLorentzRotation invrot = rot0.inverse()*rot1.inverse(); // ************************************ // // Now we construct the momenta in the // // b rest frame using _xa,_xg. // // First we construct b, then c and g, // // finally we generate a by momentum // // conservation. // // ************************************ // Lorentz5Momentum pa_brf, pb_brf(_mt), pc_brf, pg_brf; // First we set the top quark to being on-shell and at rest. // Second we set the energies of c and g, pc_brf.setE(0.5*_mt*(2.-_xa-_xg)); pg_brf.setE(0.5*_mt*_xg); // then their masses, pc_brf.setMass(_mc); pg_brf.setMass(ZERO); // Now set the z-component of c and g. For pg we simply start from // _xa and _xg, while for pc we assume it is equal to minus the sum // of the z-components of a (assumed to point in the +z direction) and g. double root=sqrt(_xa*_xa-4.*_a); pg_brf.setZ(_mt*(1.-_xa-_xg+0.5*_xa*_xg-_c+_a)/root); pc_brf.setZ(-1.*( pg_brf.z()+_mt*0.5*root)); // Now set the y-component of c and g's momenta pc_brf.setY(ZERO); pg_brf.setY(ZERO); // Now set the x-component of c and g's momenta pg_brf.setX(sqrt(sqr(pg_brf.t())-sqr(pg_brf.z()))); pc_brf.setX(-pg_brf.x()); // Momenta b,c,g are now set. Now we obtain a from momentum conservation, pa_brf = pb_brf-pc_brf-pg_brf; pa_brf.setMass(pa_brf.m()); pa_brf.rescaleEnergy(); // ************************************ // // Now we orient the momenta and boost // // them back to the original lab frame. // // ************************************ // // As in herwig6507 we assume that, in the rest frame // of b, we have aligned the W boson momentum in the // +Z direction by rot1*rot0*pa_lab, therefore // we obtain the new pa_lab by applying: // invrot*pa_brf. pa_lab = invrot*pa_brf; pb_lab = invrot*pb_brf; pc_lab = invrot*pc_brf; pg_lab = invrot*pg_brf; fs.push_back(pc_lab); fs.push_back(pa_lab); fs.push_back(pg_lab); return fs; } double SMTopDecayer::getHard(double ktb, double ktc) { // zero the variables _xg = 0.; _xa = 0.; _xc = 0.; // Get a phase space point in the dead region: double volume_factor = deadRegionxgxa(ktb,ktc); // if outside region return -1 if(volume_factor<0) return volume_factor; // Compute the weight for this phase space point: double weight = volume_factor*me(_xa,_xg)*(1.+_a-_c-_xa); // Alpha_S and colour factors - this hard wired Alpha_S needs removing. weight *= (4./3.)/Constants::pi *(_alpha->value(_mt*_mt*_xg*(1.-_xa+_a-_c) /(2.-_xg-_xa-_c))); return weight; } bool SMTopDecayer::softMatrixElementVeto(ShowerProgenitorPtr initial, ShowerParticlePtr parent,Branching br) { // check if we need to apply the full correction long id[2]={abs(initial->progenitor()->id()),abs(parent->id())}; // the initial-state correction if(id[0]==ParticleID::t&&id[1]==ParticleID::t) { Energy pt=br.kinematics->pT(); // check if hardest so far // if not just need to remove effect of enhancement bool veto(false); // if not hardest so far if(pthighestpT()) veto=!UseRandom::rndbool(1./_initialenhance); // if hardest so far do calculation else { // values of kappa and z double z(br.kinematics->z()),kappa(sqr(br.kinematics->scale()/_mt)); // parameters for the translation double w(1.-(1.-z)*(kappa-1.)),u(1.+_a-_c-(1.-z)*kappa),v(sqr(u)-4.*_a*w*z); // veto if outside phase space if(v<0.) veto=true; // otherwise calculate the weight else { v = sqrt(v); double xa((0.5*(u+v)/w+0.5*(u-v)/z)),xg((1.-z)*kappa); double f(me(xa,xg)), J(0.5*(u+v)/sqr(w)-0.5*(u-v)/sqr(z)+_a*sqr(w-z)/(v*w*z)); double wgt(f*J*2./kappa/(1.+sqr(z)-2.*z/kappa)/_initialenhance); // This next `if' prevents the hardest emission from the // top shower ever entering the so-called T2 region of the // phase space if that region is to be populated by the hard MEC. if(_useMEforT2&&xg>xgbcut(_ktb)) wgt = 0.; if(wgt>1.) { generator()->log() << "Violation of maximum for initial-state " << " soft veto in " << "SMTopDecayer::softMatrixElementVeto" << "xg = " << xg << " xa = " << xa << "weight = " << wgt << "\n"; wgt=1.; } // compute veto from weight veto = !UseRandom::rndbool(wgt); } // if not vetoed reset max if(!veto) initial->highestpT(pt); } // if vetoing reset the scale if(veto) parent->vetoEmission(br.type,br.kinematics->scale()); // return the veto return veto; } // final-state correction else if(id[0]==ParticleID::b&&id[1]==ParticleID::b) { Energy pt=br.kinematics->pT(); // check if hardest so far // if not just need to remove effect of enhancement bool veto(false); // if not hardest so far if(pthighestpT()) return !UseRandom::rndbool(1./_finalenhance); // if hardest so far do calculation // values of kappa and z double z(br.kinematics->z()),kappa(sqr(br.kinematics->scale()/_mt)); // momentum fractions double xa(1.+_a-_c-z*(1.-z)*kappa),r(0.5*(1.+_c/(1.+_a-xa))),root(sqr(xa)-4.*_a); if(root<0.) { generator()->log() << "Imaginary root for final-state veto in " << "SMTopDecayer::softMatrixElementVeto" << "\nz = " << z << "\nkappa = " << kappa << "\nxa = " << xa << "\nroot^2= " << root; parent->vetoEmission(br.type,br.kinematics->scale()); return true; } root=sqrt(root); double xg((2.-xa)*(1.-r)-(z-r)*root); // xfact (below) is supposed to equal xg/(1-z). double xfact(z*kappa/2./(z*(1.-z)*kappa+_c)*(2.-xa-root)+root); // calculate the full result double f(me(xa,xg)); // jacobian double J(z*root); double wgt(f*J*2.*kappa/(1.+sqr(z)-2.*_c/kappa/z)/sqr(xfact)/_finalenhance); if(wgt>1.) { generator()->log() << "Violation of maximum for final-state soft veto in " << "SMTopDecayer::softMatrixElementVeto" << "xg = " << xg << " xa = " << xa << "weight = " << wgt << "\n"; wgt=1.; } // compute veto from weight veto = !UseRandom::rndbool(wgt); // if vetoing reset the scale if(veto) parent->vetoEmission(br.type,br.kinematics->scale()); // return the veto return veto; } // otherwise don't veto else return !UseRandom::rndbool(1./_finalenhance); } double SMTopDecayer::me(double xw,double xg) { double prop(1.+_a-_c-xw),xg2(sqr(xg)); double lambda=sqrt(1.+_a*_a+_c*_c-2.*_a-2.*_c-2.*_a*_c); double denom=(1.-2*_a*_a+_a+_c*_a+_c*_c-2.*_c); double wgt=-_c*xg2/prop+(1.-_a+_c)*xg-(prop*(1 - xg)+xg2) +(0.5*(1.+2.*_a+_c)*sqr(prop-xg)*xg+2.*_a*prop*xg2)/denom; return wgt/(lambda*prop); } // This function is auxiliary to the xab function. double SMTopDecayer::xgbr(int toggle) { return 1.+toggle*sqrt(_a)-_c*(1.-toggle*sqrt(_a))/(1.-_a); } // This function is auxiliary to the xab function. double SMTopDecayer::ktr(double xgb, int toggle) { return 2.*xgb/ (xgb+toggle*sqrt((1.-1./_a) *(xgb-xgbr( 1)) *(xgb-xgbr(-1)))); } // Function xab determines xa (2*W energy fraction) for a given value // of xg (2*gluon energy fraction) and kappa tilde (q tilde squared over // m_top squared). Hence this function allows you to draw 1: the total // phase space volume in the xa vs xg plane 2: for a given value of // kappa tilde (i.e. starting evolution scale) the associated contour // in the xa vs xg plane (and hence the regions that either shower can // populate). This calculation is done assuming the emission came from // the top quark i.e. kappa tilde here is the q tilde squared of the TOP // quark divided by m_top squared. double SMTopDecayer::xab(double xgb, double kt, int toggle) { double xab; if(toggle==2) { // This applies for g==0.&&kt==ktr(a,c,0.,xgb,1). xab = -2.*_a*(xgb-2.)/(1.+_a-_c-xgb); } else if(toggle==1) { // This applies for kt==1&&g==0. double lambda = sqrt(sqr(xgb-1.+_a+_c)-4.*_a*_c); xab = (0.5/(kt-xgb))*(kt*(1.+_a-_c-xgb)-lambda) + (0.5/(kt+xgb*(1.-kt)))*(kt*(1.+_a-_c-xgb)+lambda); } else { // This is the form of xab FOR _g=0. double ktmktrpktmktrm = kt*kt - 4.*_a*(kt-1.)*xgb*xgb / (sqr(1.-_a-_c-xgb)-4.*_a*_c); if(fabs(kt-(2.*xgb-2.*_g)/(xgb-sqrt(xgb*xgb-4.*_g)))/kt>1.e-6) { double lambda = sqrt((sqr(1.-_a-_c-xgb)-4.*_a*_c)*ktmktrpktmktrm); xab = (0.5/(kt-xgb))*(kt*(1.+_a-_c-xgb)-lambda) + (0.5/(kt+xgb*(1.-kt)))*(kt*(1.+_a-_c-xgb)+lambda); } else { // This is the value of xa as a function of xb when kt->infinity. // Where we take any kt > (2.*xgb-2.*_g)/(xgb-sqrt(xgb*xgb-4.*_g)) // as being effectively infinite. This kt value is actually the // maximum allowed value kt can have if the phase space is calculated // without the approximation of _g=0 (massless gluon). This formula // for xab below is then valid for _g=0 AND kt=infinity only. xab = ( 2.*_c+_a*(xgb-2.) + 3.*xgb - xgb*(_c+xgb+sqrt(_a*_a-2.*(_c-xgb+1.)*_a+sqr(_c+xgb-1.))) - 2. )/2./(xgb-1.); } } if(std::isnan(xab)) { double ktmktrpktmktrm = ( sqr(xgb*kt-2.*(xgb-_g)) -kt*kt*(1.-1./_a)*(xgb-xgbr( 1)-_g/(1.+sqrt(_a))) *(xgb-xgbr(-1)-_g/(1.-sqrt(_a))) )/ (xgb*xgb-(1.-1./_a)*(xgb-xgbr( 1)-_g/(1.+sqrt(_a))) *(xgb-xgbr(-1)-_g/(1.-sqrt(_a))) ); double lambda = sqrt((xgb-1.+sqr(sqrt(_a)+sqrt(_c-_g))) *(xgb-1.+sqr(sqrt(_a)-sqrt(_c-_g)))* ktmktrpktmktrm); xab = (0.5/(kt-xgb+_g))*(kt*(1.+_a-_c+_g-xgb)-lambda) + (0.5/(kt+xgb*(1.-kt)-_g))*(kt*(1.+_a-_c+_g-xgb)+lambda); if(std::isnan(xab)) throw Exception() << "TopMECorrection::xab complex x_a value.\n" << " xgb = " << xgb << "\n" << " xab = " << xab << "\n" << " toggle = " << toggle << "\n" << " ktmktrpktmktrm = " << ktmktrpktmktrm << Exception::eventerror; } return xab; } // xgbcut is the point along the xg axis where the upper bound on the // top quark (i.e. b) emission phase space goes back on itself in the // xa vs xg plane i.e. roughly mid-way along the xg axis in // the xa vs xg Dalitz plot. double SMTopDecayer::xgbcut(double kt) { double lambda2 = 1.+_a*_a+_c*_c-2.*_a-2.*_c-2.*_a*_c; double num1 = kt*kt*(1.-_a-_c); double num2 = 2.*kt*sqrt(_a*(kt*kt*_c+lambda2*(kt-1.))); return (num1-num2)/(kt*kt-4.*_a*(kt-1.)); } double SMTopDecayer::xaccut(double kt) { return 1.+_a-_c-0.25*kt; } double SMTopDecayer::z(double xac, double kt, int toggle1, int toggle2) { double z = -1.0; if(toggle2==0) { z = (kt+toggle1*sqrt(kt*(kt-4.*(1.+_a-_c-xac))))/(2.*kt); } else if(toggle2==1) { z = ((1.+_a+_c-xac)+toggle1*(1.+_a-_c-xac)) /(2.*(1.+_a-xac)); } else if(toggle2==2) { z = 0.5; } else { throw Exception() << "Cannot determine z in SMTopDecayer::z()" << Exception::eventerror; } return z; } double SMTopDecayer::xgc(double xac, double kt, int toggle1, int toggle2) { double tiny(1.e-6); double xaToMinBoundary(xac*xac-4.*_a); if(xaToMinBoundary<0) { if(fabs(xaToMinBoundary/(1.-_a)/(1.-_a))0.&&(4.*v3+u2)<0.) output = cos( atan(sqrt(-4.*v3-u2)/u)/3.); else if(u>0.&&(4.*v3+u2)>0.) output = cosh(atanh(sqrt( 4.*v3+u2)/u)/3.); else output = cos(( atan(sqrt(-4.*v3-u2)/u) +Constants::pi)/3.); output *= 2.*sqrt(-v); } else { output = sinh(log((u+sqrt(4.*v3+u2))/(2.*sqrt(v3)))/3.); output *= 2.*sqrt(v); } if(!isfinite(output)) { throw Exception() << "TopMECorrection::xginvc0:\n" << "possible numerical instability detected.\n" << "\n v = " << v << " u = " << u << "\n4.*v3+u2 = " << 4.*v3+u2 << "\n_a = " << _a << " ma = " << sqrt(_a*_mt*_mt/GeV2) << "\n_c = " << _c << " mc = " << sqrt(_c*_mt*_mt/GeV2) << "\n_g = " << _g << " mg = " << sqrt(_g*_mt*_mt/GeV2) << Exception::eventerror; } return ( 1.+_a +output)/3.; } double SMTopDecayer::approxDeadMaxxa(double xg,double ktb,double ktc) { double maxxa(0.); double x = min(xginvc0(xg,ktc), xab(xg,(2.*xg-2.*_g)/(xg-sqrt(xg*xg-4.*_g)),0)); double y(-9999999999.); if(xg>2.*sqrt(_g)&&xg<=xgbcut(ktb)) { y = max(xab(xg,ktb,0),xab(xg,1.,1)); } else if(xg>=xgbcut(ktb)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) { y = max(xab(xg,ktr(xg,1),2),xab(xg,1.,1)); } if(xg>2.*sqrt(_g)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) { if(x>=y) { maxxa = x ; } else { maxxa = -9999999.; } } else { maxxa = -9999999.; } return maxxa; } double SMTopDecayer::approxDeadMinxa(double xg,double ktb,double ktc) { double minxa(0.); double x = min(xginvc0(xg,ktc), xab(xg,(2.*xg-2.*_g)/(xg-sqrt(xg*xg-4.*_g)),0)); double y(-9999999999.); if(xg>2.*sqrt(_g)&&xg<=xgbcut(ktb)) { y = max(xab(xg,ktb,0),xab(xg,1.,1)); } else if(xg>=xgbcut(ktb)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) { if(_useMEforT2) y = xab(xg,1.,1); else y = max(xab(xg,ktr(xg,1),2),xab(xg,1.,1)); } if(xg>2.*sqrt(_g)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) { if(x>=y) { minxa = y ; } else { minxa = 9999999.; } } else { minxa = 9999999.; } return minxa; } // This function returns true if the phase space point (xg,xa) is in the // kinematically allowed phase space. bool SMTopDecayer::inTheAllowedRegion(double xg , double xa) { bool output(true); if(xg<2.*sqrt(_g)||xg>1.-sqr(sqrt(_a)+sqrt(_c))) output = false; if(xaxab(xg,(2.*xg-2.*_g)/(xg-sqrt(xg*xg-4.*_g)),0)) output = false; return output; } // This function returns true if the phase space point (xg,xa) is in the // approximate (overestimated) dead region. bool SMTopDecayer::inTheApproxDeadRegion(double xg , double xa, double ktb, double ktc) { bool output(true); if(!inTheAllowedRegion(xg,xa)) output = false; if(xaapproxDeadMaxxa(xg,ktb,ktc)) output = false; return output; } // This function returns true if the phase space point (xg,xa) is in the // dead region. bool SMTopDecayer::inTheDeadRegion(double xg , double xa, double ktb, double ktc) { bool output(true); if(!inTheApproxDeadRegion(xg,xa,ktb,ktc)) output = false; if(xa>xaccut(ktc)) { if(xgxgc(xa,ktc, 1,0)) { output = false; } if(xg>xgc(max(xaccut(ktc),2.*sqrt(_a)),ktc,-1,2)&& xgxgmax) throw Exception() << "TopMECorrection::deadRegionxgxa:\n" << "upper xg bound is less than the lower xg bound.\n" << "\n_xg = " << _xg << "\n2.*sqrt(_g) = " << 2.*sqrt(_g) << "\n_a = " << _a << " ma = " << sqrt(_a*_mt*_mt/GeV2) << "\n_c = " << _c << " mc = " << sqrt(_c*_mt*_mt/GeV2) << "\n_g = " << _g << " mg = " << sqrt(_g*_mt*_mt/GeV2) << Exception::eventerror; double xamin(approxDeadMinxa(_xg,ktb,ktc)); double xamax(approxDeadMaxxa(_xg,ktb,ktc)); // Are the bounds sensible? If not return. if(xamax<=xamin) return -1.; _xa=1.+_a-(1.+_a-xamax)*pow((1.+_a-xamin)/(1.+_a-xamax),UseRandom::rnd()); // If outside the allowed region return -1. if(!inTheDeadRegion(_xg,_xa,ktb,ktc)) return -1.; // The integration volume for the weight double xg_vol,xa_vol; if(_xg_sampling==2.) { xg_vol = (xgmax-xgmin) / (xgmax*xgmin); } else { xg_vol = (pow(xgmax,_xg_sampling-1.)-pow(xgmin,_xg_sampling-1.)) / ((_xg_sampling-1.)*pow(xgmax*xgmin,_xg_sampling-1.)); } xa_vol = log((1.+_a-xamin)/(1.+_a-xamax)); // Here we return the integral volume factor multiplied by the part of the // weight left over which is not included in the BRACES function, i.e. // the part of _xg^-2 which is not absorbed in the integration measure. return xg_vol*xa_vol*pow(_xg,_xg_sampling-2.); } LorentzRotation SMTopDecayer::rotateToZ(Lorentz5Momentum v) { // compute the rotation matrix LorentzRotation trans; // rotate so in z-y plane trans.rotateZ(-atan2(v.y(),v.x())); // rotate so along Z trans.rotateY(-acos(v.z()/v.vect().mag())); // generate random rotation double c,s,cs; do { c = 2.*UseRandom::rnd()-1.; s = 2.*UseRandom::rnd()-1.; cs = c*c+s*s; } while(cs>1.||cs==0.); double cost=(c*c-s*s)/cs,sint=2.*c*s/cs; // apply random azimuthal rotation trans.rotateZ(atan2(sint,cost)); return trans; }