diff --git a/MatrixElement/General/GeneralQQHiggs.cc b/MatrixElement/General/GeneralQQHiggs.cc --- a/MatrixElement/General/GeneralQQHiggs.cc +++ b/MatrixElement/General/GeneralQQHiggs.cc @@ -1,644 +1,644 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEPP2QQH class. // #include "GeneralQQHiggs.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/MatrixElement/Tree2toNDiagram.h" #include "Herwig/Models/StandardModel/StandardModel.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/Cuts/Cuts.h" #include "Herwig/MatrixElement/HardVertex.h" using namespace Herwig; GeneralQQHiggs::GeneralQQHiggs() : quarkFlavour_(6), process_(0), shapeOpt_(2), mh_(), wh_(), alpha_(1.1) {} // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigGeneralQQHiggs("Herwig::GeneralQQHiggs", "Herwig.so"); void GeneralQQHiggs::Init() { static ClassDocumentation documentation ("The GeneralQQHiggs class implements the matrix elements for the " "production of the Higgs boson in association with a heavy quark-antiquark pair"); static Switch interfaceQuarkType ("QuarkType", "The type of quark", &GeneralQQHiggs::quarkFlavour_, 6, false, false); static SwitchOption interfaceQuarkTypeBottom (interfaceQuarkType, "Bottom", "Produce bottom-antibottom", 5); static SwitchOption interfaceQuarkTypeTop (interfaceQuarkType, "Top", "Produce top-antitop", 6); static Switch interfaceProcess ("Process", "Which subprocesses to include", &GeneralQQHiggs::process_, 0, false, false); static SwitchOption interfaceProcessAll (interfaceProcess, "All", "Include all subprocesses", 0); static SwitchOption interfaceProcess1 (interfaceProcess, "gg", "Include only gg -> QQbarH processes", 1); static SwitchOption interfaceProcessqbarqbarqbarqbar (interfaceProcess, "qqbar", "Include only qbar qbar -> QQbarH processes", 2); static Switch interfaceShapeOption ("ShapeScheme", "Option for the treatment of the Higgs resonance shape", &GeneralQQHiggs::shapeOpt_, 2, false, false); static SwitchOption interfaceStandardShapeFixed (interfaceShapeOption, "FixedBreitWigner", "Breit-Wigner s-channel resonance", 1); static SwitchOption interfaceStandardShapeRunning (interfaceShapeOption, "MassGenerator", "Use the mass generator to give the shape", 2); static SwitchOption interfaceStandardShapeOnShell (interfaceShapeOption, "OnShell", "Produce an on-shell Higgs boson", 0); static Parameter interfaceAlpha ("Alpha", "Power for the generation of the tranverse mass in the pT mapping", &GeneralQQHiggs::alpha_, 1.1, 0.0, 10.0, false, false, Interface::limited); } Energy2 GeneralQQHiggs::scale() const { return sHat(); } int GeneralQQHiggs::nDim() const { return 4 + int(shapeOpt_>0); } unsigned int GeneralQQHiggs::orderInAlphaS() const { return 2; } unsigned int GeneralQQHiggs::orderInAlphaEW() const { return 1; } IBPtr GeneralQQHiggs::clone() const { return new_ptr(*this); } IBPtr GeneralQQHiggs::fullclone() const { return new_ptr(*this); } void GeneralQQHiggs::setKinematics() { HwMEBase::setKinematics(); } void GeneralQQHiggs::persistentOutput(PersistentOStream & os) const { os << quarkFlavour_ << process_ << shapeOpt_ << ounit(mh_,GeV) << ounit(wh_,GeV) << hmass_ << GGGVertex_ << QQGVertex_ << QQHVertex_ << gluon_ << higgs_ << quark_ << antiquark_ << alpha_; } void GeneralQQHiggs::persistentInput(PersistentIStream & is, int) { is >> quarkFlavour_ >> process_ >> shapeOpt_ >> iunit(mh_,GeV) >> iunit(wh_,GeV) >> hmass_ >> GGGVertex_ >> QQGVertex_ >> QQHVertex_ >> gluon_ >> higgs_ >> quark_ >> antiquark_ >> alpha_; } void GeneralQQHiggs::doinit() { HwMEBase::doinit(); // stuff for the higgs mass mh_ = higgs_->mass(); wh_ = higgs_->width(); if(higgs_->massGenerator()) { hmass_=dynamic_ptr_cast(higgs_->massGenerator()); } if(shapeOpt_==2&&!hmass_) throw InitException() << "If using the mass generator for the line shape in GeneralQQHiggs::doinit()" << "the mass generator must be an instance of the GenericMassGenerator class" << Exception::runerror; // get the vertex pointers from the SM object tcHwSMPtr hwsm= dynamic_ptr_cast(standardModel()); // do the initialisation if(!hwsm) throw InitException() << "Wrong type of StandardModel object in " << "GeneralQQHiggs::doinit() the Herwig" << " version must be used" << Exception::runerror; GGGVertex_ = hwsm->vertexGGG(); QQGVertex_ = hwsm->vertexFFG(); // get the particle data objects gluon_=getParticleData(ParticleID::g); for(int ix=1;ix<=6;++ix) { quark_.push_back( getParticleData( ix)); antiquark_.push_back(getParticleData(-ix)); } } bool GeneralQQHiggs::generateKinematics(const double * r) { jacobian(1.); // CMS energy Energy rs = sqrt(sHat()); // quark mass Energy mq(quark_[quarkFlavour_-1]->mass()); // generate the higgs mass Energy mh(mh_); if(shapeOpt_!=0) { Energy mhmax = min(rs-2.*mq,higgs_->massMax()); Energy mhmin = max(ZERO ,higgs_->massMin()); if(mhmax<=mhmin) return false; double rhomin = atan2((sqr(mhmin)-sqr(mh_)), mh_*wh_); double rhomax = atan2((sqr(mhmax)-sqr(mh_)), mh_*wh_); mh = sqrt(mh_*wh_*tan(rhomin+r[4]*(rhomax-rhomin))+sqr(mh_)); jacobian(jacobian()*(rhomax-rhomin)); } if(rs ZERO ) { double ctm = 1.0 - sqr(ptmin/p1); if ( ctm <= 0.0 ) return false; ctmin = max(ctmin, -sqrt(ctm)); ctmax = min(ctmax, sqrt(ctm)); } double cos1 = getCosTheta(ctmin,ctmax,r[0]); double sin1(sqrt(1.-sqr(cos1))); double phi1 = Constants::twopi*UseRandom::rnd(); Lorentz5Momentum p13(sin1*p1*cos(phi1),sin1*p1*sin(phi1),cos1*p1, sqrt(sqr(p1)+m132),m13); meMomenta()[3].setVect(Momentum3(-sin1*p1*cos(phi1),-sin1*p1*sin(phi1),-cos1*p1)); meMomenta()[3].setMass(mq); meMomenta()[3].rescaleEnergy(); bool test=Kinematics::twoBodyDecay(p13,mq,mh,-1.+2*r[2],r[3]*Constants::twopi, meMomenta()[2],meMomenta()[4]); if(!test) return false; m232 = (meMomenta()[3]+meMomenta()[4]).m2(); double D = 2./(pow(sqr(mq)/(m132-sqr(mq)),alpha_)+ pow(sqr(mq)/(m232-sqr(mq)),alpha_)); jacobian(0.5*jacobian()*rs/m13*sqr(mq)*D*(rhomax-rhomin)/sHat()); } // second branch else { double rtemp = 2.*(r[1]-0.5); double rho = rhomin+rtemp*(rhomax-rhomin); if(alpha_==0) m232 = sqr(mq)*rho; else if(alpha_==1) m232 = sqr(mq)*(exp(rho)+1.); else m232 = sqr(mq)*(pow(rho,1./(1.-alpha_))+1.); Energy m23 = sqrt(m232); try { p1 = SimplePhaseSpace::getMagnitude(sHat(), m23, mq); p2 = SimplePhaseSpace::getMagnitude(m232,mq,mh); - } catch ( ImpossibleKinematics ) { + } catch ( ImpossibleKinematics & e ) { return false; } Energy ptmin = lastCuts().minKT(mePartonData()[2]); double ctmin = -1.0, ctmax = 1.0; if ( ptmin > ZERO ) { double ctm = 1.0 - sqr(ptmin/p1); if ( ctm <= 0.0 ) return false; ctmin = max(ctmin, -sqrt(ctm)); ctmax = min(ctmax, sqrt(ctm)); } double cos1 = getCosTheta(ctmin,ctmax,r[0]); double sin1(sqrt(1.-sqr(cos1))); double phi1 = Constants::twopi*UseRandom::rnd(); Lorentz5Momentum p23(-sin1*p1*cos(phi1),-sin1*p1*sin(phi1),-cos1*p1, sqrt(sqr(p1)+m232),m23); meMomenta()[2].setVect(Momentum3(sin1*p1*cos(phi1),sin1*p1*sin(phi1),cos1*p1)); meMomenta()[2].setMass(mq); meMomenta()[2].rescaleEnergy(); bool test=Kinematics::twoBodyDecay(p23,mq,mh,-1.+2*r[2],r[3]*Constants::twopi, meMomenta()[3],meMomenta()[4]); if(!test) return false; m132 = (meMomenta()[2]+meMomenta()[4]).m2(); double D = 2./(pow(sqr(mq)/(m132-sqr(mq)),alpha_)+ pow(sqr(mq)/(m232-sqr(mq)),alpha_)); jacobian(0.5*jacobian()*rs/m23*sqr(mq)*D*(rhomax-rhomin)/sHat()); } // calculate jacobian jacobian(0.125*jacobian()*p1*p2/sHat()); // check cuts vector out; tcPDVector tout; for(unsigned int ix=2;ix<5;++ix) { out .push_back(meMomenta()[ix]); tout.push_back(mePartonData()[ix]); } return lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]); } CrossSection GeneralQQHiggs::dSigHatDR() const { using Constants::pi; // jacobian factor for the higgs InvEnergy2 bwfact(ZERO); Energy moff = meMomenta()[4].mass(); if(shapeOpt_==1) { bwfact = mePartonData()[4]->generateWidth(moff)*moff/pi/ (sqr(sqr(moff)-sqr(mh_))+sqr(mh_*wh_)); } else if(shapeOpt_==2) { bwfact = hmass_->BreitWignerWeight(moff); } double jac1 = shapeOpt_==0 ? 1. : double(bwfact*(sqr(sqr(moff)-sqr(mh_))+sqr(mh_*wh_))/(mh_*wh_)); return sqr(hbarc)*me2()*jacobian()*jac1/sHat()/pow(Constants::twopi,3); } void GeneralQQHiggs::getDiagrams() const { tPDPtr Q = quark_[quarkFlavour_-1]; tPDPtr QB = antiquark_[quarkFlavour_-1]; // gg -> q qbar h0 subprocesses if(process_==0||process_==1) { // first t-channel add(new_ptr((Tree2toNDiagram(3), gluon_, QB, gluon_, 1, Q, 4, Q , 2, QB, 4, higgs_, -1))); add(new_ptr((Tree2toNDiagram(4), gluon_, QB, QB, gluon_, 1, Q, 3, QB, 2, higgs_, -2))); add(new_ptr((Tree2toNDiagram(3),gluon_,QB,gluon_, 1, Q, 2, QB, 5, QB, 5, higgs_, -3))); // interchange add(new_ptr((Tree2toNDiagram(3),gluon_,Q,gluon_, 2, Q, 4, Q , 1, QB, 4, higgs_, -4))); add(new_ptr((Tree2toNDiagram(4),gluon_,Q,Q,gluon_, 3, Q, 1, QB, 2, higgs_, -5))); add(new_ptr((Tree2toNDiagram(3),gluon_,Q,gluon_, 2, Q, 1, QB, 5, QB, 5, higgs_, -6))); // s-channel add(new_ptr((Tree2toNDiagram(2),gluon_,gluon_, 1, gluon_, 3, Q, 4, Q, 3, QB, 4, higgs_, -7))); add(new_ptr((Tree2toNDiagram(2),gluon_,gluon_, 1, gluon_, 3,Q, 3, QB, 5, QB, 5, higgs_, -8))); } // q qbar -> q qbar if(process_==0||process_==2) { for(unsigned int ix=1;ix<5;++ix) { // gluon s-channel add(new_ptr((Tree2toNDiagram(2),quark_[ix-1], antiquark_[ix-1], 1, gluon_, 3, Q, 4, Q , 3, QB, 4, higgs_, -9))); add(new_ptr((Tree2toNDiagram(2),quark_[ix-1], antiquark_[ix-1], 1, gluon_, 3, Q, 3, QB, 5, QB, 5, higgs_, -10))); } } } double GeneralQQHiggs::me2() const { // total matrix element double me(0.); // gg initiated processes if(mePartonData()[0]->id()==ParticleID::g) { VectorWaveFunction g1w(meMomenta()[0],mePartonData()[0],incoming); VectorWaveFunction g2w(meMomenta()[1],mePartonData()[1],incoming); SpinorBarWaveFunction qw(meMomenta()[2],mePartonData()[2],outgoing); SpinorWaveFunction qbarw(meMomenta()[3],mePartonData()[3],outgoing); ScalarWaveFunction higgs(meMomenta()[4],mePartonData()[4],1.,outgoing); vector g1,g2; vector q; vector qbar; for(unsigned int ix=0;ix<2;++ix) { g1w.reset(2*ix);g1.push_back(g1w); g2w.reset(2*ix);g2.push_back(g2w); qw.reset(ix);q.push_back(qw); qbarw.reset(ix);qbar.push_back(qbarw); } // calculate the matrix element me=ggME(g1,g2,q,qbar,higgs,0); } // q qbar initiated else { SpinorWaveFunction q1w(meMomenta()[0],mePartonData()[0],incoming); SpinorBarWaveFunction q2w(meMomenta()[1],mePartonData()[1],incoming); SpinorBarWaveFunction q3w(meMomenta()[2],mePartonData()[2],outgoing); SpinorWaveFunction q4w(meMomenta()[3],mePartonData()[3],outgoing); ScalarWaveFunction higgs(meMomenta()[4],mePartonData()[4],1.,outgoing); vector q1,q4; vector q2,q3; for(unsigned int ix=0;ix<2;++ix) { q1w.reset(ix);q1.push_back(q1w); q2w.reset(ix);q2.push_back(q2w); q3w.reset(ix);q3.push_back(q3w); q4w.reset(ix);q4.push_back(q4w); } // calculate the matrix element me = qqME(q1,q2,q3,q4,higgs,0); } return me*sHat()*UnitRemoval::InvE2; } double GeneralQQHiggs::ggME(vector &g1, vector &g2, vector & q, vector & qbar, ScalarWaveFunction & hwave, unsigned int iflow) const { // scale Energy2 mt(scale()); Energy mass = q[0].mass(); // matrix element to be stored if(iflow!=0) me_.reset(ProductionMatrixElement(PDT::Spin1,PDT::Spin1, PDT::Spin1Half,PDT::Spin1Half, PDT::Spin0)); // calculate the matrix element double output(0.),sumflow[2]={0.,0.}; double sumdiag[8]={0.,0.,0.,0.,0.,0.,0.,0.}; Complex diag[8],flow[2]; VectorWaveFunction interv; SpinorWaveFunction inters,QBoff; SpinorBarWaveFunction intersb,Qoff; for(unsigned int ihel1=0;ihel1<2;++ihel1) { for(unsigned int ihel2=0;ihel2<2;++ihel2) { interv = GGGVertex_->evaluate(mt,5,gluon_,g1[ihel1],g2[ihel2]); for(unsigned int ohel1=0;ohel1<2;++ohel1) { Qoff = QQHVertex_->evaluate(mt,3,q[ohel1].particle(), q[ohel1],hwave,mass); for(unsigned int ohel2=0;ohel2<2;++ohel2) { QBoff = QQHVertex_->evaluate(mt,3,qbar[ohel2].particle(), qbar[ohel2],hwave,mass); // 1st diagram inters = QQGVertex_->evaluate(mt,1,qbar[ohel2].particle(), qbar[ohel2],g2[ihel2],mass); diag[0] = QQGVertex_->evaluate(mt,inters,Qoff,g1[ihel1]); // 2nd diagram intersb = QQGVertex_->evaluate(mt,1,q[ohel1].particle(), q[ohel1],g1[ihel1],mass); diag[1] = QQHVertex_->evaluate(mt,inters,intersb,hwave); // 3rd diagram diag[2] = QQGVertex_->evaluate(mt,QBoff,intersb,g2[ihel2]); // 4th diagram inters = QQGVertex_->evaluate(mt,1,qbar[ohel2].particle(), qbar[ohel2],g1[ihel1],mass); diag[3] = QQGVertex_->evaluate(mt,inters,Qoff,g2[ihel2]); // 5th diagram intersb = QQGVertex_->evaluate(mt,1,q[ohel1].particle(), q[ohel1],g2[ihel2],mass); diag[4] = QQHVertex_->evaluate(mt,inters,intersb,hwave); // 6th diagram diag[5] = QQGVertex_->evaluate(mt,QBoff,intersb,g1[ihel1]); // 7th diagram diag[6] = QQGVertex_->evaluate(mt,qbar[ohel2],Qoff ,interv); // 8th diagram diag[7] = QQGVertex_->evaluate(mt,QBoff ,q[ohel1],interv); // colour flows flow[0]=diag[0]+diag[1]+diag[2]+(diag[6]+diag[7]); flow[1]=diag[3]+diag[4]+diag[5]-(diag[6]+diag[7]); // sums for(unsigned int ix=0;ix<8;++ix) sumdiag[ix] += norm(diag[ix]); for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]); // total output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1]) -0.25*flow[0]*conj(flow[1])); // store the me if needed if(iflow!=0) me_(2*ihel1,2*ihel2,ohel1,ohel2,0)=flow[iflow-1]; } } } } // select a colour flow flow_ = 1 + UseRandom::rnd2(sumflow[0],sumflow[1]); if(flow_==1) sumdiag[0]=sumdiag[1]=sumdiag[2]=0.; else sumdiag[3]=sumdiag[4]=sumdiag[5]=0.; // select a diagram from that flow double prob = UseRandom::rnd(); for(unsigned int ix=0;ix<8;++ix) { if(prob<=sumdiag[ix]) { diagram_=1+ix; break; } prob -= sumdiag[ix]; } // final part of colour and spin factors return output/48.; } double GeneralQQHiggs::qqME(vector & q1, vector & q2, vector & q3, vector & q4, ScalarWaveFunction & hwave, unsigned int iflow) const { // scale Energy2 mt(scale()); Energy mass = q3[0].mass(); // matrix element to be stored if(iflow!=0) me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1Half,PDT::Spin1Half, PDT::Spin0)); // calculate the matrix element double output(0.),sumdiag[2]={0.,0.}; Complex diag[2]; VectorWaveFunction interv; SpinorWaveFunction QBoff; SpinorBarWaveFunction Qoff; for(unsigned int ihel1=0;ihel1<2;++ihel1) { for(unsigned int ihel2=0;ihel2<2;++ihel2) { interv = QQGVertex_->evaluate(mt,5,gluon_,q1[ihel1],q2[ihel2]); for(unsigned int ohel1=0;ohel1<2;++ohel1) { Qoff = QQHVertex_->evaluate(mt,3,q3[ohel1].particle(), q3[ohel1],hwave,mass); for(unsigned int ohel2=0;ohel2<2;++ohel2) { QBoff = QQHVertex_->evaluate(mt,3,q4[ohel2].particle(), q4[ohel2],hwave,mass); // 1st diagram diag[0] = QQGVertex_->evaluate(mt,q4[ohel2],Qoff,interv); // 2nd diagram diag[1] = QQGVertex_->evaluate(mt,QBoff,q3[ohel1],interv); // sum of diagrams for(unsigned int ix=0;ix<2;++ix) sumdiag[ix] += norm(diag[ix]); diag[0] += diag[1]; output += norm(diag[0]); if(iflow!=0) me_(ihel1,ihel2,ohel1,ohel2,0) = diag[0]; } } } } // only 1 colour flow flow_=1; // select a diagram diagram_ = 9+UseRandom::rnd2(sumdiag[0],sumdiag[1]); // final part of colour and spin factors return output/18.; } Selector GeneralQQHiggs::colourGeometries(tcDiagPtr diag) const { // colour lines for gg -> Q Qbar H static const ColourLines cgg[10]= {ColourLines("1 4 5, -1 -2 3 , -3 -6 "), ColourLines("1 5 , -1 -2 -3 4, -4 -6 "), ColourLines("1 4 , -1 -2 3 , -3 -5 -6"), ColourLines("3 4 5, 1 2 -3 , -1 -6 "), ColourLines("4 5 , 1 2 3 -4, -1 -6"), ColourLines("3 4 , 1 2 -3 , -1 -5 -6"), ColourLines("1 3 4 5, -1 2, -2 -3 -6"), ColourLines("2 3 4 5, 1 -2, -1 -3 -6"), ColourLines("1 3 4, -1 2, -2 -3 -5 -6"), ColourLines("2 3 4, 1 -2, -1 -3 -5 -6")}; // colour lines for q qbar -> Q Qbar H static const ColourLines cqq[2]= {ColourLines("1 3 4 5, -2 -3 -6"), ColourLines("1 3 4 , -2 -3 -5 -6")}; // select the colour flow (as all ready picked just insert answer) Selector sel; switch(abs(diag->id())) { // gg -> q qbar subprocess case 1: case 2: case 3: case 4: case 5: case 6: sel.insert(1.0, &cgg[abs(diag->id())-1]); break; case 7: sel.insert(1.0, &cgg[5 + flow_]); break; case 8: sel.insert(1.0, &cgg[7 + flow_]); break; // q qbar -> q qbar subprocess case 9: case 10: sel.insert(1.0, &cqq[abs(diag->id())-9]); break; } return sel; } Selector GeneralQQHiggs::diagrams(const DiagramVector & diags) const { // select the diagram, this is easy for us as we have already done it Selector sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ) { if(diags[i]->id()==-int(diagram_)) sel.insert(1.0, i); else sel.insert(0., i); } return sel; } void GeneralQQHiggs::constructVertex(tSubProPtr sub) { // extract the particles in the hard process ParticleVector hard; hard.push_back(sub->incoming().first); hard.push_back(sub->incoming().second); for(unsigned int ix=0;ix<3;++ix) hard.push_back(sub->outgoing()[ix]); // identify the process and calculate the matrix element if(hard[0]->id()<0) swap(hard[0],hard[1]); if(hard[2]->id()==higgs_->id()) swap(hard[2],hard[4]); if(hard[3]->id()==higgs_->id()) swap(hard[3],hard[4]); if(hard[2]->id()<0) swap(hard[2],hard[3]); if(hard[2]->id()<0) swap(hard[2],hard[3]); if(hard[0]->id()==ParticleID::g) { vector g1,g2; vector q; vector qbar; // off-shell wavefunctions for the spin correlations VectorWaveFunction( g1,hard[0],incoming,false,true,true); VectorWaveFunction( g2,hard[1],incoming,false,true,true); SpinorBarWaveFunction(q ,hard[2],outgoing,true ,true); SpinorWaveFunction( qbar,hard[3],outgoing,true ,true); ScalarWaveFunction hwave( hard[4],outgoing,true); g1[1]=g1[2];g2[1]=g2[2]; ggME(g1,g2,q,qbar,hwave,flow_); } // q qbar -> Q Qbar Higgs else { vector q1,q4; vector q2,q3; // off-shell for spin correlations SpinorWaveFunction( q1,hard[0],incoming,false,true); SpinorBarWaveFunction(q2,hard[1],incoming,false,true); SpinorBarWaveFunction(q3,hard[2],outgoing,true ,true); SpinorWaveFunction( q4,hard[3],outgoing,true ,true); ScalarWaveFunction hwave( hard[4],outgoing,true); qqME(q1,q2,q3,q4,hwave,flow_); } // construct the vertex HardVertexPtr hardvertex=new_ptr(HardVertex()); // set the matrix element for the vertex hardvertex->ME(me_); // set the pointers and to and from the vertex for(unsigned int ix=0;ix<5;++ix) hard[ix]->spinInfo()->productionVertex(hardvertex); } void GeneralQQHiggs::setProcessInfo(unsigned int quark, PDPtr hin, AbstractFFSVertexPtr vertex, unsigned int shapeOpt, unsigned int proc) { quarkFlavour_ = quark; higgs_ = hin; QQHVertex_ = vertex; process_ = proc; shapeOpt_ = shapeOpt; } diff --git a/MatrixElement/Hadron/MEMinBias.cc b/MatrixElement/Hadron/MEMinBias.cc --- a/MatrixElement/Hadron/MEMinBias.cc +++ b/MatrixElement/Hadron/MEMinBias.cc @@ -1,164 +1,164 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEMinBias class. // #include "MEMinBias.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" //#include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/MatrixElement/Tree2toNDiagram.h" void MEMinBias::getDiagrams() const { int maxflav(2); // Pomeron data tcPDPtr pom = getParticleData(990); for ( int i = 1; i <= maxflav; ++i ) { for( int j=1; j <= i; ++j){ tcPDPtr q1 = getParticleData(i); tcPDPtr q1b = q1->CC(); tcPDPtr q2 = getParticleData(j); tcPDPtr q2b = q2->CC(); // For each flavour we add: //qq -> qq add(new_ptr((Tree2toNDiagram(3), q1, pom, q2, 1, q1, 2, q2, -1))); //qqb -> qqb add(new_ptr((Tree2toNDiagram(3), q1, pom, q2b, 1, q1, 2, q2b, -2))); //qbqb -> qbqb add(new_ptr((Tree2toNDiagram(3), q1b, pom, q2b, 1, q1b, 2, q2b, -3))); } } } Energy2 MEMinBias::scale() const { return sqr(10*GeV); } int MEMinBias::nDim() const { return 0; } void MEMinBias::setKinematics() { HwMEBase::setKinematics(); // Always call the base class method first. } bool MEMinBias::generateKinematics(const double *) { // generate the masses of the particles for ( int i = 2, N = meMomenta().size(); i < N; ++i ) { meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->generateMass()); } Energy q = ZERO; try { q = SimplePhaseSpace:: getMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass()); - } catch ( ImpossibleKinematics ) { + } catch ( ImpossibleKinematics & e ) { return false; } Energy pt = ZERO; meMomenta()[2].setVect(Momentum3( pt, pt, q)); meMomenta()[3].setVect(Momentum3(-pt, -pt, -q)); meMomenta()[2].rescaleEnergy(); meMomenta()[3].rescaleEnergy(); jacobian(1.0); return true; } double MEMinBias::me2() const { //tuned so it gives the correct normalization for xmin = 0.11 return csNorm_*(sqr(generator()->maximumCMEnergy())/GeV2); } CrossSection MEMinBias::dSigHatDR() const { return me2()*jacobian()/sHat()*sqr(hbarc); } unsigned int MEMinBias::orderInAlphaS() const { return 2; } unsigned int MEMinBias::orderInAlphaEW() const { return 0; } Selector MEMinBias::diagrams(const DiagramVector & diags) const { Selector sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ) sel.insert(1.0, i); return sel; } Selector MEMinBias::colourGeometries(tcDiagPtr diag) const { static ColourLines qq("1 4, 3 5"); static ColourLines qqb("1 4, -3 -5"); static ColourLines qbqb("-1 -4, -3 -5"); Selector sel; switch(diag->id()){ case -1: sel.insert(1.0, &qq); break; case -2: sel.insert(1.0, &qqb); break; case -3: sel.insert(1.0, &qbqb); break; } return sel; } IBPtr MEMinBias::clone() const { return new_ptr(*this); } IBPtr MEMinBias::fullclone() const { return new_ptr(*this); } ClassDescription MEMinBias::initMEMinBias; // Definition of the static class description member. void MEMinBias::persistentOutput(PersistentOStream & os) const { os << csNorm_; } void MEMinBias::persistentInput(PersistentIStream & is, int) { is >> csNorm_; } void MEMinBias::Init() { static ClassDocumentation documentation ("There is no documentation for the MEMinBias class"); static Parameter interfacecsNorm ("csNorm", "Normalization of the min-bias cross section.", &MEMinBias::csNorm_, 1.0, 0.0, 100.0, false, false, Interface::limited); } diff --git a/MatrixElement/Hadron/MEPP2HiggsJet.cc b/MatrixElement/Hadron/MEPP2HiggsJet.cc --- a/MatrixElement/Hadron/MEPP2HiggsJet.cc +++ b/MatrixElement/Hadron/MEPP2HiggsJet.cc @@ -1,875 +1,875 @@ // -*- C++ -*- // // MEPP2HiggsJet.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 MEPP2HiggsJet class. // #include "MEPP2HiggsJet.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/Cuts/Cuts.h" #include "ThePEG/MatrixElement/Tree2toNDiagram.h" #include "Herwig/MatrixElement/HardVertex.h" using namespace Herwig; const Complex MEPP2HiggsJet::_epsi = Complex(0.,-1.e-20); IBPtr MEPP2HiggsJet::clone() const { return new_ptr(*this); } IBPtr MEPP2HiggsJet::fullclone() const { return new_ptr(*this); } unsigned int MEPP2HiggsJet::orderInAlphaS() const { return 3; } unsigned int MEPP2HiggsJet::orderInAlphaEW() const { return 1; } void MEPP2HiggsJet::persistentOutput(PersistentOStream & os) const { os << _shapeopt << _maxflavour << _process << _minloop << _maxloop << _massopt << ounit(_mh,GeV) << ounit(_wh,GeV) << _hmass; } void MEPP2HiggsJet::persistentInput(PersistentIStream & is, int) { is >> _shapeopt >> _maxflavour >> _process >> _minloop >> _maxloop >> _massopt >> iunit(_mh,GeV) >> iunit(_wh,GeV) >> _hmass; } ClassDescription MEPP2HiggsJet::initMEPP2HiggsJet; // Definition of the static class description member. void MEPP2HiggsJet::Init() { static ClassDocumentation documentation ("The MEPP2HiggsJet class implements the matrix elements for" " Higgs+Jet production in hadron-hadron collisions.", "The theoretical calculations of \\cite{Baur:1989cm} and \\cite{Ellis:1987xu}" " were used for the Higgs+jet matrix element in hadron-hadron collisions.", "\\bibitem{Baur:1989cm} U.~Baur and E.~W.~N.~Glover," "Nucl.\\ Phys.\\ B {\\bf 339} (1990) 38.\n" "\\bibitem{Ellis:1987xu} R.~K.~Ellis, I.~Hinchliffe, M.~Soldate and " "J.~J.~van der Bij, Nucl.\\ Phys.\\ B {\\bf 297} (1988) 221."); static Parameter interfaceMaximumFlavour ("MaximumFlavour", "The maximum flavour of the quarks in the process", &MEPP2HiggsJet::_maxflavour, 5, 1, 5, false, false, Interface::limited); static Switch interfaceShapeOption ("ShapeScheme", "Option for the treatment of the Higgs resonance shape", &MEPP2HiggsJet::_shapeopt, 1, false, false); static SwitchOption interfaceStandardShapeFixed (interfaceShapeOption, "FixedBreitWigner", "Breit-Wigner s-channel resonanse", 1); static SwitchOption interfaceStandardShapeRunning (interfaceShapeOption, "MassGenerator", "Use the mass generator to give the shape", 2); static Switch interfaceProcess ("Process", "Which subprocesses to include", &MEPP2HiggsJet::_process, 0, false, false); static SwitchOption interfaceProcessAll (interfaceProcess, "All", "Include all subprocesses", 0); static SwitchOption interfaceProcess1 (interfaceProcess, "qqbar", "Only include the incoming q qbar subprocess", 1); static SwitchOption interfaceProcessqg (interfaceProcess, "qg", "Only include the incoming qg subprocess", 2); static SwitchOption interfaceProcessqbarg (interfaceProcess, "qbarg", "Only include the incoming qbar g subprocess", 3); static SwitchOption interfaceProcessgg (interfaceProcess, "gg", "Only include the incoming gg subprocess", 4); static Parameter interfaceMinimumInLoop ("MinimumInLoop", "The minimum flavour of the quarks to include in the loops", &MEPP2HiggsJet::_minloop, 6, 4, 6, false, false, Interface::limited); static Parameter interfaceMaximumInLoop ("MaximumInLoop", "The maximum flavour of the quarks to include in the loops", &MEPP2HiggsJet::_maxloop, 6, 4, 6, false, false, Interface::limited); static Switch interfaceMassOption ("MassOption", "Option for the treatment of the masses in the loop diagrams", &MEPP2HiggsJet::_massopt, 0, false, false); static SwitchOption interfaceMassOptionFull (interfaceMassOption, "Full", "Include the full mass dependence", 0); static SwitchOption interfaceMassOptionLarge (interfaceMassOption, "Large", "Use the heavy mass limit", 1); } bool MEPP2HiggsJet::generateKinematics(const double * r) { Energy ptmin = max(lastCuts().minKT(mePartonData()[2]), lastCuts().minKT(mePartonData()[3])); Energy e = sqrt(sHat())/2.0; // generate the mass of the higgs boson Energy2 mhmax2 = sHat()-4.*ptmin*e; Energy2 mhmin2 =ZERO; if(mhmax2<=mhmin2) return false; double rhomin = atan2((mhmin2-sqr(_mh)), _mh*_wh); double rhomax = atan2((mhmax2-sqr(_mh)), _mh*_wh); Energy mh = sqrt(_mh*_wh*tan(rhomin+r[1]*(rhomax-rhomin))+sqr(_mh)); // assign masses if(mePartonData()[2]->id()!=ParticleID::h0) { meMomenta()[2].setMass(ZERO); meMomenta()[3].setMass(mh); } else { meMomenta()[3].setMass(ZERO); meMomenta()[2].setMass(mh); } Energy q = ZERO; try { q = SimplePhaseSpace:: getMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass()); } - catch ( ImpossibleKinematics ) { + catch ( ImpossibleKinematics & e ) { return false; } Energy2 m22 = meMomenta()[2].mass2(); Energy2 m32 = meMomenta()[3].mass2(); Energy2 e0e2 = 2.0*e*sqrt(sqr(q) + m22); Energy2 e1e2 = 2.0*e*sqrt(sqr(q) + m22); Energy2 e0e3 = 2.0*e*sqrt(sqr(q) + m32); Energy2 e1e3 = 2.0*e*sqrt(sqr(q) + m32); Energy2 pq = 2.0*e*q; double ctmin = -1.0,ctmax = 1.0; Energy2 thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[2]); if ( thmin > ZERO ) ctmax = min(ctmax, (e0e2 - m22 - thmin)/pq); thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[2]); if ( thmin > ZERO ) ctmin = max(ctmin, (thmin + m22 - e1e2)/pq); thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[3]); if ( thmin > ZERO ) ctmax = min(ctmax, (e1e3 - m32 - thmin)/pq); thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[3]); if ( thmin > ZERO ) ctmin = max(ctmin, (thmin + m32 - e0e3)/pq); if ( ptmin > ZERO ) { double ctm = 1.0 - sqr(ptmin/q); if ( ctm <= 0.0 ) return false; ctmin = max(ctmin, -sqrt(ctm)); ctmax = min(ctmax, sqrt(ctm)); } if ( ctmin >= ctmax ) return false; double cth = getCosTheta(ctmin, ctmax, r); Energy pt = q*sqrt(1.0-sqr(cth)); phi(rnd(2.0*Constants::pi)); meMomenta()[2].setX(pt*sin(phi())); meMomenta()[2].setY(pt*cos(phi())); meMomenta()[2].setZ(q*cth); meMomenta()[3].setX(-pt*sin(phi())); meMomenta()[3].setY(-pt*cos(phi())); meMomenta()[3].setZ(-q*cth); meMomenta()[2].rescaleEnergy(); meMomenta()[3].rescaleEnergy(); vector out(2); out[0] = meMomenta()[2]; out[1] = meMomenta()[3]; tcPDVector tout(2); tout[0] = mePartonData()[2]; tout[1] = mePartonData()[3]; if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) ) return false; tHat(pq*cth + m22 - e0e2); uHat(m22 + m32 - sHat() - tHat()); // main piece jacobian((pq/sHat())*Constants::pi*jacobian()); // mass piece jacobian((rhomax-rhomin)*jacobian()); return true; } void MEPP2HiggsJet::getDiagrams() const { tcPDPtr h0=getParticleData(ParticleID::h0); tcPDPtr g =getParticleData(ParticleID::g); tcPDPtr q[6],qb[6]; for(int ix=0;ix H g if(_process==0||_process==1) {for(unsigned int ix=0;ix<_maxflavour;++ix) {add(new_ptr((Tree2toNDiagram(2), q[ix], qb[ix], 1, g , 3, h0, 3, g, -1)));}} // q g -> H g if(_process==0||_process==2) {for(unsigned int ix=0;ix<_maxflavour;++ix) {add(new_ptr((Tree2toNDiagram(3), q[ix], g, g, 2, h0, 1, q[ix], -2)));}} // qbar g -> H qbar if(_process==0||_process==3) {for(unsigned int ix=0;ix<_maxflavour;++ix) {add(new_ptr((Tree2toNDiagram(3), qb[ix], g, g, 2, h0, 1, qb[ix], -3)));}} // g g -> H g if(_process==0||_process==4) { // t channel add(new_ptr((Tree2toNDiagram(3), g, g, g, 1, h0, 2, g, -4))); // u channel add(new_ptr((Tree2toNDiagram(3), g, g, g, 2, h0, 1, g, -5))); // s channel add(new_ptr((Tree2toNDiagram(2), g, g, 1, g , 3, h0, 3, g, -6))); } } Energy2 MEPP2HiggsJet::scale() const { return meMomenta()[2].perp2()+ meMomenta()[2].m2(); } double MEPP2HiggsJet::me2() const { useMe(); double output(0.); // g g to H g if(mePartonData()[0]->id()==ParticleID::g&&mePartonData()[1]->id()==ParticleID::g) { // order of the particles unsigned int ih(2),ig(3); if(mePartonData()[3]->id()==ParticleID::h0){ig=2;ih=3;} VectorWaveFunction glin1(meMomenta()[ 0],mePartonData()[ 0],incoming); VectorWaveFunction glin2(meMomenta()[ 1],mePartonData()[ 1],incoming); ScalarWaveFunction hout(meMomenta()[ih],mePartonData()[ih],outgoing); VectorWaveFunction glout(meMomenta()[ig],mePartonData()[ig],outgoing); vector g1,g2,g4; for(unsigned int ix=0;ix<2;++ix) { glin1.reset(2*ix);g1.push_back(glin1); glin2.reset(2*ix);g2.push_back(glin2); glout.reset(2*ix);g4.push_back(glout); } // calculate the matrix element output = ggME(g1,g2,hout,g4,false); } // qg -> H q else if(mePartonData()[0]->id()>0&&mePartonData()[1]->id()==ParticleID::g) { // order of the particles unsigned int iq(0),iqb(3),ih(2),ig(1); if(mePartonData()[0]->id()==ParticleID::g){iq=1;ig=0;} if(mePartonData()[3]->id()==ParticleID::h0){iqb=2;ih=3;} // calculate the spinors and polarization vectors vector fin; vector fout; vector gin; SpinorWaveFunction qin (meMomenta()[iq ],mePartonData()[iq ],incoming); VectorWaveFunction glin(meMomenta()[ig ],mePartonData()[ig ],incoming); ScalarWaveFunction hout(meMomenta()[ih ],mePartonData()[ih ],outgoing); SpinorBarWaveFunction qout(meMomenta()[iqb],mePartonData()[iqb],outgoing); for(unsigned int ix=0;ix<2;++ix) { qin.reset(ix) ; fin.push_back( qin); qout.reset(ix) ;fout.push_back(qout); glin.reset(2*ix); gin.push_back(glin); } // calculate the matrix element output = qgME(fin,gin,hout,fout,false); } // qbar g -> H q else if(mePartonData()[0]->id()<0&&mePartonData()[1]->id()==ParticleID::g) { // order of the particles unsigned int iq(0),iqb(3),ih(2),ig(1); if(mePartonData()[0]->id()==ParticleID::g){iq=1;ig=0;} if(mePartonData()[3]->id()==ParticleID::h0){iqb=2;ih=3;} // calculate the spinors and polarization vectors vector fin; vector fout; vector gin; SpinorBarWaveFunction qin (meMomenta()[iq ],mePartonData()[iq ],incoming); VectorWaveFunction glin(meMomenta()[ig ],mePartonData()[ig ],incoming); ScalarWaveFunction hout(meMomenta()[ih ],mePartonData()[ih ],outgoing); SpinorWaveFunction qout(meMomenta()[iqb],mePartonData()[iqb],outgoing); for(unsigned int ix=0;ix<2;++ix) { qin.reset(ix) ; fin.push_back( qin); qout.reset(ix) ;fout.push_back(qout); glin.reset(2*ix); gin.push_back(glin); } // calculate the matrix element output = qbargME(fin,gin,hout,fout,false); } // q qbar to H g else if(mePartonData()[0]->id()==-mePartonData()[1]->id()) { // order of the particles unsigned int iq(0),iqb(1),ih(2),ig(3); if(mePartonData()[0]->id()<0){iq=1;iqb=0;} if(mePartonData()[2]->id()==ParticleID::g){ig=2;ih=3;} // calculate the spinors and polarization vectors vector fin; vector ain; vector gout; SpinorWaveFunction qin (meMomenta()[iq ],mePartonData()[iq ],incoming); SpinorBarWaveFunction qbin(meMomenta()[iqb],mePartonData()[iqb],incoming); ScalarWaveFunction hout(meMomenta()[ih ],mePartonData()[ih ],outgoing); VectorWaveFunction glout(meMomenta()[ig ],mePartonData()[ig ],outgoing); for(unsigned int ix=0;ix<2;++ix) { qin.reset(ix) ; fin.push_back( qin); qbin.reset(ix) ; ain.push_back( qbin); glout.reset(2*ix);gout.push_back(glout); } // calculate the matrix element output = qqbarME(fin,ain,hout,gout,false); } else throw Exception() << "Unknown subprocess in MEPP2HiggsJet::me2()" << Exception::runerror; // return the answer return output; } double MEPP2HiggsJet::qqbarME(vector & fin, vector & ain, ScalarWaveFunction & hout, vector & gout, bool calc) 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 higgs // 1 outgoing gluon // me to be returned ProductionMatrixElement newme(PDT::Spin1Half,PDT::Spin1Half, PDT::Spin0,PDT::Spin1); // get the kinematic invariants Energy2 s(sHat()),u(uHat()),t(tHat()),mh2(hout.m2()),et(scale()); // calculate the loop function complex A5 = Energy2(); for ( int ix=_minloop; ix<=_maxloop; ++ix ) { // full mass dependance if(_massopt==0) { Energy2 mf2=sqr(getParticleData(ix)->mass()); A5+= mf2*(4.+4.*double(s/(u+t))*(W1(s,mf2)-W1(mh2,mf2)) +(1.-4.*double(mf2/(u+t)))*(W2(s,mf2)-W2(mh2,mf2))); } // infinite mass limit else { A5+=2.*(s-mh2)/3.; } } // multiply by the rest of the form factors using Constants::pi; double g(sqrt(4.*pi*SM().alphaEM(mh2)/SM().sin2ThetaW())); double gs(sqrt(4.*pi*SM().alphaS(et))); Energy mw(getParticleData(ParticleID::Wplus)->mass()); complex A5c = A5 * Complex(0.,1.)*g*sqr(gs)*gs/(32.*s*sqr(pi)*mw); // compute the matrix element LorentzPolarizationVectorE fcurrent; complex fdotp; complex epsdot[2]; Complex diag; Lorentz5Momentum ps(fin[0].momentum()+ain[0].momentum()); ps.rescaleMass(); for(unsigned int ix=0;ix<2;++ix){epsdot[ix]=gout[ix].wave().dot(ps);} Energy2 denom(-ps*gout[0].momentum()); LorentzSpinorBar atemp; double output(0.); for(unsigned int ihel1=0;ihel1<2;++ihel1) { for(unsigned int ihel2=0;ihel2<2;++ihel2) { // compute the fermion current atemp=ain[ihel2].wave(); fcurrent=UnitRemoval::E*fin[ihel1].wave().vectorCurrent(atemp); fdotp = -(fcurrent.dot(gout[0].momentum())); for(unsigned int ghel=0;ghel<2;++ghel) { // calculate the matrix element diag=A5c*(fcurrent.dot(gout[ghel].wave()) -fdotp*epsdot[ghel]/denom); // calculate the matrix element output+=real(diag*conj(diag)); if(calc) newme(ihel1,ihel2,0,2*ghel)=diag; } } } // test with glover form // final colour/spin factors if(calc) _me.reset(newme); return output/9.; } double MEPP2HiggsJet::qgME(vector & fin, vector & gin, ScalarWaveFunction & hout, vector & fout,bool calc) const { // the particles should be in the order // for the incoming // 0 incoming fermion (u spinor) // 1 incoming gluon // for the outgoing // 0 outgoing higgs // 1 outgoing fermion (ubar spinor) // me to be returned ProductionMatrixElement newme(PDT::Spin1Half,PDT::Spin1, PDT::Spin0,PDT::Spin1Half); // get the kinematic invariants Energy2 s(sHat()),u(uHat()),t(tHat()),mh2(hout.m2()),et(scale()); // calculate the loop function complex A5 = Energy2(); for(int ix=_minloop;ix<=_maxloop;++ix) { if(_massopt==0) { Energy2 mf2=sqr(getParticleData(ix)->mass()); A5+= mf2*(4.+4.*double(u/(s+t))*(W1(u,mf2)-W1(mh2,mf2)) +(1.-4.*double(mf2/(s+t)))*(W2(u,mf2)-W2(mh2,mf2))); } else { A5+=2.*(u-mh2)/3.; } } // multiply by the rest of the form factors using Constants::pi; double g(sqrt(4.*pi*SM().alphaEM(mh2)/SM().sin2ThetaW())); double gs(sqrt(4.*pi*SM().alphaS(et))); Energy mw(getParticleData(ParticleID::Wplus)->mass()); complex A5c =A5*Complex(0.,1.)*g*sqr(gs)*gs/(32.*u*sqr(pi)*mw); // compute the matrix element LorentzPolarizationVectorE fcurrent; complex fdotp; complex epsdot[2]; Complex diag; Lorentz5Momentum pu(fin[0].momentum()+fout[0].momentum()); pu.rescaleMass(); for(unsigned int ix=0;ix<2;++ix){epsdot[ix]=gin[ix].wave().dot(pu);} Energy2 denom(pu*gin[0].momentum()); LorentzSpinorBar atemp; double output(0.); for(unsigned int ihel=0;ihel<2;++ihel) { for(unsigned int ohel=0;ohel<2;++ohel) { // compute the fermion current atemp=fout[ohel].wave(); fcurrent=UnitRemoval::E*fin[ihel].wave().vectorCurrent(atemp); fdotp=fcurrent.dot(gin[0].momentum()); for(unsigned int ghel=0;ghel<2;++ghel) { // calculate the matrix element diag=A5c*(fcurrent.dot(gin[ghel].wave())-fdotp*epsdot[ghel]/denom); // calculate the matrix element output+=real(diag*conj(diag)); if(calc) newme(ihel,2*ghel,0,ohel)=diag; } } } // final colour/spin factors if(calc) _me.reset(newme); return output/24.; } double MEPP2HiggsJet::qbargME(vector & fin, vector & gin, ScalarWaveFunction & hout, vector & fout,bool calc) const { // the particles should be in the order // for the incoming // 0 incoming antifermion (vbar spinor) // 1 incoming gluon // for the outgoing // 0 outgoing higgs // 1 outgoing antifermion (v spinor) // me to be returned ProductionMatrixElement newme(PDT::Spin1Half,PDT::Spin1, PDT::Spin0,PDT::Spin1Half); // get the kinematic invariants Energy2 s(sHat()),u(uHat()),t(tHat()),mh2(hout.m2()),et(scale()); // calculate the loop function complex A5 = Energy2(); for(int ix=_minloop;ix<=_maxloop;++ix) { if(_massopt==0) { Energy2 mf2=sqr(getParticleData(ix)->mass()); A5+= mf2*(4.+4.*double(u/(s+t))*(W1(u,mf2)-W1(mh2,mf2)) +(1.-4.*double(mf2/(s+t)))*(W2(u,mf2)-W2(mh2,mf2))); } else { A5+=2.*(u-mh2)/3.; } } // multiply by the rest of the form factors using Constants::pi; double g(sqrt(4.*pi*SM().alphaEM(mh2)/SM().sin2ThetaW())); double gs(sqrt(4.*pi*SM().alphaS(et))); Energy mw(getParticleData(ParticleID::Wplus)->mass()); complex A5c = A5*Complex(0.,1.)*g*sqr(gs)*gs/(32.*u*sqr(pi)*mw); // compute the matrix element LorentzPolarizationVectorE fcurrent; complex fdotp; complex epsdot[2]; Complex diag; Lorentz5Momentum pu(fin[0].momentum()+fout[0].momentum()); pu.rescaleMass(); for(unsigned int ix=0;ix<2;++ix){epsdot[ix]=gin[ix].wave().dot(pu);} Energy2 denom(pu*gin[0].momentum()); LorentzSpinorBar atemp; double output(0.); for(unsigned int ihel=0;ihel<2;++ihel) { for(unsigned int ohel=0;ohel<2;++ohel) { // compute the fermion current atemp=fin[ihel].wave(); fcurrent=UnitRemoval::E*fout[ohel].wave().vectorCurrent(atemp); fdotp=fcurrent.dot(gin[0].momentum()); for(unsigned int ghel=0;ghel<2;++ghel) { // calculate the matrix element diag=A5c*(fcurrent.dot(gin[ghel].wave())-fdotp*epsdot[ghel]/denom); // calculate the matrix element output+=real(diag*conj(diag)); if(calc) newme(ihel,2*ghel,0,ohel)=diag; } } } // final colour/spin factors if(calc) _me.reset(newme); return output/24.; } Selector MEPP2HiggsJet::diagrams(const DiagramVector & diags) const { Selector sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ) { if(abs(diags[i]->id())<4) sel.insert(1.0, i); else sel.insert(_diagwgt[abs(diags[i]->id())-4], i); } return sel; } Selector MEPP2HiggsJet::colourGeometries(tcDiagPtr diag) const { // colour lines for q qbar -> h0 g static const ColourLines cqqbar("1 3 5,-2 -3 -5"); // colour lines for q g -> h0 q static const ColourLines cqg("1 2 -3, 3 -2 5"); // colour lines for qbar q -> h0 qbar static const ColourLines cqbarg("-1 -2 3, -3 2 -5"); // colour lines for g g -> h0 g static const ColourLines cgg[6]={ColourLines("1 2 5, -3 -5, 3 -2 -1"), ColourLines("-1 -2 -5, 3 5, -3 2 1"), ColourLines("1 5, -1 -2 3, -3 2 -5"), ColourLines("-1 -5, 1 2 -3, 3 -2 5"), ColourLines("1 3 5, -5 -3 -2, 2 -1"), ColourLines("-1 -3 -5, 5 3 2 ,-2 1")}; // select the colour flow Selector sel; if ( diag->id() == -1) sel.insert(1.0, &cqqbar); else if ( diag->id() == -2) sel.insert(1.0, &cqg); else if ( diag->id() == -3) sel.insert(1.0, &cqbarg); else { sel.insert(0.5, &cgg[2*(abs(diag->id())-4) ]); sel.insert(0.5, &cgg[2*(abs(diag->id())-4)+1]); } // return the answer return sel; } double MEPP2HiggsJet::ggME(vector g1, vector g2, ScalarWaveFunction & hout, vector g4, bool calc) const { // the particles should be in the order // for the incoming // 0 first incoming gluon // 1 second incoming gluon // for the outgoing // 0 outgoing higgs // 1 outgoing gluon // me to be returned ProductionMatrixElement newme(PDT::Spin1,PDT::Spin1, PDT::Spin0,PDT::Spin1); // get the kinematic invariants Energy2 s(sHat()),u(uHat()),t(tHat()),mh2(hout.m2()),et(scale()); // calculate the loop functions Complex A4stu(0.),A2stu(0.),A2tsu(0.),A2ust(0.); Complex A5s(0.),A5t(0.),A5u(0.); for(int ix=_minloop;ix<=_maxloop;++ix) { Energy2 mf2=sqr(getParticleData(ix)->mass()); // loop functions if(_massopt==0) { A4stu+=A4(s,t,u,mf2); A2stu+=A2(s,t,u,mf2); A2tsu+=A2(u,s,t,mf2); A2ust+=A2(t,s,u,mf2); A5s+= mf2/s*(4.+4.*double(s/(u+t))*(W1(s,mf2)-W1(mh2,mf2)) +(1.-4.*double(mf2/(u+t)))*(W2(s,mf2)-W2(mh2,mf2))); A5t+= mf2/t*(4.+4.*double(t/(s+u))*(W1(t,mf2)-W1(mh2,mf2)) +(1.-4.*double(mf2/(s+u)))*(W2(t,mf2)-W2(mh2,mf2))); A5u+= mf2/u*(4.+4.*double(u/(s+t))*(W1(u,mf2)-W1(mh2,mf2)) +(1.-4.*double(mf2/(s+t)))*(W2(u,mf2)-W2(mh2,mf2))); } else { A4stu=-1./3.; A2stu=-sqr(s/mh2)/3.; A2tsu=-sqr(t/mh2)/3.; A2ust=-sqr(u/mh2)/3.; A5s+=2.*(s-mh2)/3./s; A5t+=2.*(t-mh2)/3./t; A5u+=2.*(u-mh2)/3./u; } } Complex A3stu=0.5*(A2stu+A2ust+A2tsu-A4stu); // compute the dot products for the matrix element complex eps[3][4][2]; Energy2 pdot[4][4]; pdot[0][0]=ZERO; pdot[0][1]= g1[0].momentum()*g2[0].momentum(); pdot[0][2]=-1.*g1[0].momentum()*g4[0].momentum(); pdot[0][3]=-1.*g1[0].momentum()*hout.momentum(); pdot[1][0]= pdot[0][1]; pdot[1][1]= ZERO; pdot[1][2]=-1.*g2[0].momentum()*g4[0].momentum(); pdot[1][3]=-1.*g2[0].momentum()*hout.momentum(); pdot[2][0]= pdot[0][2]; pdot[2][1]= pdot[1][2]; pdot[2][2]= ZERO; pdot[2][3]= g4[0].momentum()*hout.momentum(); pdot[3][0]=pdot[0][3]; pdot[3][1]=pdot[1][3]; pdot[3][2]=pdot[2][3]; pdot[3][3]=mh2; for(unsigned int ix=0;ix<2;++ix) { eps[0][0][ix]=InvEnergy(); eps[0][1][ix]=g1[ix].wave().dot(g2[0].momentum())/pdot[0][1]; eps[0][2][ix]=-1.*g1[ix].wave().dot(g4[0].momentum())/pdot[0][2]; eps[0][3][ix]=-1.*g1[ix].wave().dot(hout.momentum())/ pdot[0][3]; eps[1][0][ix]=g2[ix].wave().dot(g1[0].momentum())/ pdot[1][0]; eps[1][1][ix]=InvEnergy(); eps[1][2][ix]=-1.*g2[ix].wave().dot(g4[0].momentum())/pdot[1][2]; eps[1][3][ix]=-1.*g2[ix].wave().dot(hout.momentum())/ pdot[1][3]; eps[2][0][ix]=g4[ix].wave().dot(g1[0].momentum())/ pdot[2][0]; eps[2][1][ix]=g4[ix].wave().dot(g2[0].momentum())/ pdot[2][1]; eps[2][2][ix]=InvEnergy(); eps[2][3][ix]=-1.*g4[ix].wave().dot(hout.momentum())/ pdot[2][3]; } // prefactors using Constants::pi; double g(sqrt(4.*pi*SM().alphaEM(mh2)/SM().sin2ThetaW())); double gs(sqrt(4.*pi*SM().alphaS(et))); Energy mw(getParticleData(ParticleID::Wplus)->mass()); Energy3 pre=g*sqr(mh2)*gs*sqr(gs)/(32.*sqr(pi)*mw); // compute the matrix element double output(0.); Complex diag[4],wdot[3][3]; _diagwgt[0]=0.; _diagwgt[1]=0.; _diagwgt[2]=0.; for(unsigned int ihel1=0;ihel1<2;++ihel1) { for(unsigned int ihel2=0;ihel2<2;++ihel2) { for(unsigned int ohel=0;ohel<2;++ohel) { wdot[0][1]=g1[ihel1].wave().dot(g2[ihel2].wave()); wdot[0][2]=g1[ihel1].wave().dot(g4[ohel ].wave()); wdot[1][0]=wdot[0][1]; wdot[1][2]=g2[ihel2].wave().dot(g4[ohel ].wave()); wdot[2][0]=wdot[0][2]; wdot[2][1]=wdot[1][2]; // last piece diag[3]=pre*A3stu*(eps[0][2][ihel1]*eps[1][0][ihel2]*eps[2][1][ohel]- eps[0][1][ihel1]*eps[1][2][ihel2]*eps[2][0][ohel]+ (eps[2][0][ohel ]-eps[2][1][ohel ])*wdot[0][1]/pdot[0][1]+ (eps[1][2][ihel2]-eps[1][0][ihel2])*wdot[0][2]/pdot[0][2]+ (eps[0][1][ihel1]-eps[0][2][ihel1])*wdot[1][2]/pdot[1][2]); // first piece diag[3]+=pre*( +A2stu*(eps[0][1][ihel1]*eps[1][0][ihel2]-wdot[0][1]/pdot[0][1])* (eps[2][0][ohel ]-eps[2][1][ohel ]) +A2ust*(eps[0][2][ihel1]*eps[2][0][ohel ]-wdot[0][2]/pdot[0][2])* (eps[1][2][ihel2]-eps[1][0][ihel2]) +A2tsu*(eps[1][2][ihel2]*eps[2][1][ohel ]-wdot[1][2]/pdot[1][2])* (eps[0][1][ihel1]-eps[0][2][ihel1]) ); output+=real(diag[3]*conj(diag[3])); // matrix element if needed if(calc) newme(2*ihel1,2*ihel2,0,2*ohel)=diag[3]; // different diagrams diag[0] = A5t*UnitRemoval::InvE*(-eps[0][3][ihel1]* (-2.*eps[2][1][ohel ]*eps[1][0][ihel2]*pdot[2][1]*pdot[1][0] -2.*eps[1][2][ihel2]*eps[2][0][ohel ]*pdot[1][2]*pdot[2][0] +wdot[1][2]*(pdot[0][1]+pdot[0][2])) -2.*eps[2][1][ohel ]*pdot[2][1]*wdot[0][1] -2.*eps[1][2][ihel2]*pdot[1][2]*wdot[0][2] +wdot[1][2]*(eps[0][1][ihel1]*pdot[0][1]+ eps[0][2][ihel1]*pdot[0][2])); diag[1] = A5u*UnitRemoval::InvE*(-eps[1][3][ihel2]* (+2.*eps[0][1][ihel1]*eps[2][0][ohel ]*pdot[0][1]*pdot[2][0] +2.*eps[0][2][ihel1]*eps[2][1][ohel ]*pdot[0][2]*pdot[2][1] -wdot[0][2]*(pdot[1][0]+pdot[1][2])) +2.*eps[2][0][ohel ]*pdot[2][0]*wdot[0][1] +2.*eps[0][2][ihel1]*pdot[0][2]*wdot[2][1] -wdot[0][2]*(eps[1][0][ihel2]*pdot[1][0]+ eps[1][2][ihel2]*pdot[1][2])); diag[2] = A5s*UnitRemoval::InvE*(-eps[2][3][ohel ]* (+2.*eps[0][1][ihel1]*eps[1][2][ihel2]*pdot[0][1]*pdot[1][2] -2.*eps[1][0][ihel2]*eps[0][2][ihel1]*pdot[1][0]*pdot[1][3] +wdot[0][1]*(pdot[2][0]-pdot[2][1])) +2.*eps[0][1][ihel1]*pdot[0][1]*wdot[1][2] -2.*eps[1][0][ihel2]*pdot[1][0]*wdot[0][2] +wdot[0][1]*(eps[2][0][ohel]*pdot[2][0]- eps[2][1][ohel]*pdot[2][1])); _diagwgt[0]+=real(diag[0]*conj(diag[0])); _diagwgt[1]+=real(diag[1]*conj(diag[1])); _diagwgt[2]+=real(diag[2]*conj(diag[2])); } } } // final colour and spin factors if(calc){_me.reset(newme);} return 3.*output/32.; } void MEPP2HiggsJet::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]); // ensure correct order or particles if((hard[0]->id()==ParticleID::g&&hard[1]->id()!=ParticleID::g)|| (hard[0]->id()<0&&hard[1]->id()<6)) swap(hard[0],hard[1]); if(hard[2]->id()!=ParticleID::h0) swap(hard[2],hard[3]); // different processes // g g to H g if(hard[0]->id()==ParticleID::g) { vector g1,g2,g4; VectorWaveFunction(g1,hard[0],incoming,false,true,true); VectorWaveFunction(g2,hard[1],incoming,false,true,true); VectorWaveFunction(g4,hard[3],outgoing,true ,true,true); ScalarWaveFunction hout(hard[2],outgoing,true); g1[1]=g1[2];g2[1]=g2[2];g4[1]=g4[2]; ggME(g1,g2,hout,g4,true); } // qg -> H q else if(hard[0]->id()>0&&hard[1]->id()==ParticleID::g) { vector g2; vector qin; vector qout; SpinorWaveFunction( qin,hard[0],incoming,false,true); VectorWaveFunction( g2,hard[1],incoming,false,true,true); SpinorBarWaveFunction(qout,hard[3],outgoing,true ,true); ScalarWaveFunction hout(hard[2],outgoing,true); g2[1]=g2[2]; qgME(qin,g2,hout,qout,true); } // qbar g -> H q else if(hard[0]->id()<0&&hard[1]->id()==ParticleID::g) { vector g2; vector qin; vector qout; SpinorBarWaveFunction( qin,hard[0],incoming,false,true); VectorWaveFunction( g2,hard[1],incoming,false,true,true); SpinorWaveFunction( qout,hard[3],outgoing,true ,true); ScalarWaveFunction hout(hard[2],outgoing,true); g2[1]=g2[2]; qbargME(qin,g2,hout,qout,true); } // q qbar to H g else if(hard[0]->id()==-hard[1]->id()) { vector qbar; vector q; vector g4; SpinorWaveFunction( q ,hard[0],incoming,false,true); SpinorBarWaveFunction(qbar,hard[1],incoming,false,true); VectorWaveFunction( g4,hard[3],outgoing,true ,true,true); ScalarWaveFunction hout(hard[2],outgoing,true); g4[1]=g4[2]; qqbarME(q,qbar,hout,g4,true); } else throw Exception() << "Unknown subprocess in MEPP2HiggsJet::constructVertex()" << Exception::runerror; // construct the vertex HardVertexPtr hardvertex=new_ptr(HardVertex()); // set the matrix element for the vertex hardvertex->ME(_me); // set the pointers and to and from the vertex for(unsigned int ix=0;ix<4;++ix) { (hard[ix]->spinInfo())-> productionVertex(hardvertex); } } int MEPP2HiggsJet::nDim() const { return 2; } void MEPP2HiggsJet::doinit() { ME2to2Base::doinit(); tcPDPtr h0=getParticleData(ParticleID::h0); _mh = h0->mass(); _wh = h0->generateWidth(_mh); if(h0->massGenerator()) { _hmass=dynamic_ptr_cast(h0->massGenerator()); } if(_shapeopt==2&&!_hmass) throw InitException() << "If using the mass generator for the line shape in MEPP2HiggsJet::doinit()" << "the mass generator must be an instance of the GenericMassGenerator class" << Exception::runerror; } CrossSection MEPP2HiggsJet::dSigHatDR() const { using Constants::pi; InvEnergy2 bwfact; Energy moff = mePartonData()[2]->id()==ParticleID::h0 ? meMomenta()[2].mass() : meMomenta()[3].mass(); if(_shapeopt==1) { tcPDPtr h0 = mePartonData()[2]->id()==ParticleID::h0 ? mePartonData()[2] : mePartonData()[3]; bwfact = h0->generateWidth(moff)*moff/pi/ (sqr(sqr(moff)-sqr(_mh))+sqr(_mh*_wh)); } else { bwfact = _hmass->BreitWignerWeight(moff); } return me2()*jacobian()/(16.0*sqr(Constants::pi)*sHat())*sqr(hbarc)* (sqr(sqr(moff)-sqr(_mh))+sqr(_mh*_wh))/(_mh*_wh)*bwfact; } diff --git a/MatrixElement/Hadron/MEPP2QQHiggs.cc b/MatrixElement/Hadron/MEPP2QQHiggs.cc --- a/MatrixElement/Hadron/MEPP2QQHiggs.cc +++ b/MatrixElement/Hadron/MEPP2QQHiggs.cc @@ -1,633 +1,633 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEPP2QQH class. // #include "MEPP2QQHiggs.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/MatrixElement/Tree2toNDiagram.h" #include "Herwig/Models/StandardModel/StandardModel.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/Cuts/Cuts.h" #include "Herwig/MatrixElement/HardVertex.h" using namespace Herwig; MEPP2QQHiggs::MEPP2QQHiggs() : quarkFlavour_(6), process_(0), shapeOpt_(2), mh_(), wh_(), alpha_(1.1) {} ClassDescription MEPP2QQHiggs::initMEPP2QQHiggs; // Definition of the static class description member. void MEPP2QQHiggs::Init() { static ClassDocumentation documentation ("The MEPP2QQHiggs class implements the matrix elements for the " "production of the Higgs boson in association with a heavy quark-antiquark pair"); static Switch interfaceQuarkType ("QuarkType", "The type of quark", &MEPP2QQHiggs::quarkFlavour_, 6, false, false); static SwitchOption interfaceQuarkTypeBottom (interfaceQuarkType, "Bottom", "Produce bottom-antibottom", 5); static SwitchOption interfaceQuarkTypeTop (interfaceQuarkType, "Top", "Produce top-antitop", 6); static Switch interfaceProcess ("Process", "Which subprocesses to include", &MEPP2QQHiggs::process_, 0, false, false); static SwitchOption interfaceProcessAll (interfaceProcess, "All", "Include all subprocesses", 0); static SwitchOption interfaceProcess1 (interfaceProcess, "gg", "Include only gg -> QQbarH processes", 1); static SwitchOption interfaceProcessqbarqbarqbarqbar (interfaceProcess, "qqbar", "Include only qbar qbar -> QQbarH processes", 2); static Switch interfaceShapeOption ("ShapeScheme", "Option for the treatment of the Higgs resonance shape", &MEPP2QQHiggs::shapeOpt_, 2, false, false); static SwitchOption interfaceStandardShapeFixed (interfaceShapeOption, "FixedBreitWigner", "Breit-Wigner s-channel resonance", 1); static SwitchOption interfaceStandardShapeRunning (interfaceShapeOption, "MassGenerator", "Use the mass generator to give the shape", 2); static SwitchOption interfaceStandardShapeOnShell (interfaceShapeOption, "OnShell", "Produce an on-shell Higgs boson", 0); static Parameter interfaceAlpha ("Alpha", "Power for the generation of the tranverse mass in the pT mapping", &MEPP2QQHiggs::alpha_, 1.1, 0.0, 10.0, false, false, Interface::limited); } Energy2 MEPP2QQHiggs::scale() const { return sHat(); // return sqr(mePartonData()[2]->mass()+mePartonData()[3]->mass()+ // mePartonData()[4]->mass()); } int MEPP2QQHiggs::nDim() const { return 4 + int(shapeOpt_>0); } unsigned int MEPP2QQHiggs::orderInAlphaS() const { return 2; } unsigned int MEPP2QQHiggs::orderInAlphaEW() const { return 1; } IBPtr MEPP2QQHiggs::clone() const { return new_ptr(*this); } IBPtr MEPP2QQHiggs::fullclone() const { return new_ptr(*this); } void MEPP2QQHiggs::setKinematics() { HwMEBase::setKinematics(); } void MEPP2QQHiggs::persistentOutput(PersistentOStream & os) const { os << quarkFlavour_ << process_ << shapeOpt_ << ounit(mh_,GeV) << ounit(wh_,GeV) << hmass_ << GGGVertex_ << QQGVertex_ << QQHVertex_ << gluon_ << higgs_ << quark_ << antiquark_ << alpha_; } void MEPP2QQHiggs::persistentInput(PersistentIStream & is, int) { is >> quarkFlavour_ >> process_ >> shapeOpt_ >> iunit(mh_,GeV) >> iunit(wh_,GeV) >> hmass_ >> GGGVertex_ >> QQGVertex_ >> QQHVertex_ >> gluon_ >> higgs_ >> quark_ >> antiquark_ >> alpha_; } void MEPP2QQHiggs::doinit() { HwMEBase::doinit(); // stuff for the higgs mass higgs_=getParticleData(ParticleID::h0); mh_ = higgs_->mass(); wh_ = higgs_->width(); if(higgs_->massGenerator()) { hmass_=dynamic_ptr_cast(higgs_->massGenerator()); } if(shapeOpt_==2&&!hmass_) throw InitException() << "If using the mass generator for the line shape in MEPP2QQHiggs::doinit()" << "the mass generator must be an instance of the GenericMassGenerator class" << Exception::runerror; // get the vertex pointers from the SM object tcHwSMPtr hwsm= dynamic_ptr_cast(standardModel()); // do the initialisation if(!hwsm) throw InitException() << "Wrong type of StandardModel object in " << "MEPP2QQHiggs::doinit() the Herwig" << " version must be used" << Exception::runerror; GGGVertex_ = hwsm->vertexGGG(); QQGVertex_ = hwsm->vertexFFG(); QQHVertex_ = hwsm->vertexFFH(); // get the particle data objects gluon_=getParticleData(ParticleID::g); for(int ix=1;ix<=6;++ix) { quark_.push_back( getParticleData( ix)); antiquark_.push_back(getParticleData(-ix)); } } bool MEPP2QQHiggs::generateKinematics(const double * r) { jacobian(1.); // CMS energy Energy rs = sqrt(sHat()); // quark mass Energy mq(quark_[quarkFlavour_-1]->mass()); // generate the higgs mass Energy mh(mh_); if(shapeOpt_!=0) { Energy mhmax = min(rs-2.*mq,higgs_->massMax()); Energy mhmin = max(ZERO ,higgs_->massMin()); if(mhmax<=mhmin) return false; double rhomin = atan2((sqr(mhmin)-sqr(mh_)), mh_*wh_); double rhomax = atan2((sqr(mhmax)-sqr(mh_)), mh_*wh_); mh = sqrt(mh_*wh_*tan(rhomin+r[4]*(rhomax-rhomin))+sqr(mh_)); jacobian(jacobian()*(rhomax-rhomin)); } if(rs ZERO ) { double ctm = 1.0 - sqr(ptmin/p1); if ( ctm <= 0.0 ) return false; ctmin = max(ctmin, -sqrt(ctm)); ctmax = min(ctmax, sqrt(ctm)); } double cos1 = getCosTheta(ctmin,ctmax,r[0]); double sin1(sqrt(1.-sqr(cos1))); double phi1 = Constants::twopi*UseRandom::rnd(); Lorentz5Momentum p13(sin1*p1*cos(phi1),sin1*p1*sin(phi1),cos1*p1, sqrt(sqr(p1)+m132),m13); meMomenta()[3].setVect(Momentum3(-sin1*p1*cos(phi1),-sin1*p1*sin(phi1),-cos1*p1)); meMomenta()[3].setMass(mq); meMomenta()[3].rescaleEnergy(); bool test=Kinematics::twoBodyDecay(p13,mq,mh,-1.+2*r[2],r[3]*Constants::twopi, meMomenta()[2],meMomenta()[4]); if(!test) return false; m232 = (meMomenta()[3]+meMomenta()[4]).m2(); double D = 2./(pow(sqr(mq)/(m132-sqr(mq)),alpha_)+ pow(sqr(mq)/(m232-sqr(mq)),alpha_)); jacobian(0.5*jacobian()*rs/m13*sqr(mq)*D*(rhomax-rhomin)/sHat()); } // second branch else { double rtemp = 2.*(r[1]-0.5); double rho = rhomin+rtemp*(rhomax-rhomin); if(alpha_==0) m232 = sqr(mq)*rho; else if(alpha_==1) m232 = sqr(mq)*(exp(rho)+1.); else m232 = sqr(mq)*(pow(rho,1./(1.-alpha_))+1.); Energy m23 = sqrt(m232); try { p1 = SimplePhaseSpace::getMagnitude(sHat(), m23, mq); p2 = SimplePhaseSpace::getMagnitude(m232,mq,mh); - } catch ( ImpossibleKinematics ) { + } catch ( ImpossibleKinematics & e ) { return false; } Energy ptmin = lastCuts().minKT(mePartonData()[2]); double ctmin = -1.0, ctmax = 1.0; if ( ptmin > ZERO ) { double ctm = 1.0 - sqr(ptmin/p1); if ( ctm <= 0.0 ) return false; ctmin = max(ctmin, -sqrt(ctm)); ctmax = min(ctmax, sqrt(ctm)); } double cos1 = getCosTheta(ctmin,ctmax,r[0]); double sin1(sqrt(1.-sqr(cos1))); double phi1 = Constants::twopi*UseRandom::rnd(); Lorentz5Momentum p23(-sin1*p1*cos(phi1),-sin1*p1*sin(phi1),-cos1*p1, sqrt(sqr(p1)+m232),m23); meMomenta()[2].setVect(Momentum3(sin1*p1*cos(phi1),sin1*p1*sin(phi1),cos1*p1)); meMomenta()[2].setMass(mq); meMomenta()[2].rescaleEnergy(); bool test=Kinematics::twoBodyDecay(p23,mq,mh,-1.+2*r[2],r[3]*Constants::twopi, meMomenta()[3],meMomenta()[4]); if(!test) return false; m132 = (meMomenta()[2]+meMomenta()[4]).m2(); double D = 2./(pow(sqr(mq)/(m132-sqr(mq)),alpha_)+ pow(sqr(mq)/(m232-sqr(mq)),alpha_)); jacobian(0.5*jacobian()*rs/m23*sqr(mq)*D*(rhomax-rhomin)/sHat()); } // calculate jacobian jacobian(0.125*jacobian()*p1*p2/sHat()); // check cuts vector out; tcPDVector tout; for(unsigned int ix=2;ix<5;++ix) { out .push_back(meMomenta()[ix]); tout.push_back(mePartonData()[ix]); } return lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]); } CrossSection MEPP2QQHiggs::dSigHatDR() const { using Constants::pi; // jacobian factor for the higgs InvEnergy2 bwfact; Energy moff = meMomenta()[4].mass(); if(shapeOpt_==1) { bwfact = mePartonData()[4]->generateWidth(moff)*moff/pi/ (sqr(sqr(moff)-sqr(mh_))+sqr(mh_*wh_)); } else if(shapeOpt_==2) { bwfact = hmass_->BreitWignerWeight(moff); } double jac1 = shapeOpt_==0 ? 1. : double(bwfact*(sqr(sqr(moff)-sqr(mh_))+sqr(mh_*wh_))/(mh_*wh_)); return sqr(hbarc)*me2()*jacobian()*jac1/sHat()/pow(Constants::twopi,3); } void MEPP2QQHiggs::getDiagrams() const { tPDPtr Q = quark_[quarkFlavour_-1]; tPDPtr QB = antiquark_[quarkFlavour_-1]; // gg -> q qbar h0 subprocesses if(process_==0||process_==1) { // first t-channel add(new_ptr((Tree2toNDiagram(3), gluon_, QB, gluon_, 1, Q, 4, Q , 2, QB, 4, higgs_, -1))); add(new_ptr((Tree2toNDiagram(4), gluon_, QB, QB, gluon_, 1, Q, 3, QB, 2, higgs_, -2))); add(new_ptr((Tree2toNDiagram(3),gluon_,QB,gluon_, 1, Q, 2, QB, 5, QB, 5, higgs_, -3))); // interchange add(new_ptr((Tree2toNDiagram(3),gluon_,Q,gluon_, 2, Q, 4, Q , 1, QB, 4, higgs_, -4))); add(new_ptr((Tree2toNDiagram(4),gluon_,Q,Q,gluon_, 3, Q, 1, QB, 2, higgs_, -5))); add(new_ptr((Tree2toNDiagram(3),gluon_,Q,gluon_, 2, Q, 1, QB, 5, QB, 5, higgs_, -6))); // s-channel add(new_ptr((Tree2toNDiagram(2),gluon_,gluon_, 1, gluon_, 3, Q, 4, Q, 3, QB, 4, higgs_, -7))); add(new_ptr((Tree2toNDiagram(2),gluon_,gluon_, 1, gluon_, 3,Q, 3, QB, 5, QB, 5, higgs_, -8))); } // q qbar -> q qbar if(process_==0||process_==2) { for(unsigned int ix=1;ix<5;++ix) { // gluon s-channel add(new_ptr((Tree2toNDiagram(2),quark_[ix-1], antiquark_[ix-1], 1, gluon_, 3, Q, 4, Q , 3, QB, 4, higgs_, -9))); add(new_ptr((Tree2toNDiagram(2),quark_[ix-1], antiquark_[ix-1], 1, gluon_, 3, Q, 3, QB, 5, QB, 5, higgs_, -10))); } } } double MEPP2QQHiggs::me2() const { // total matrix element double me(0.); // gg initiated processes if(mePartonData()[0]->id()==ParticleID::g) { VectorWaveFunction g1w(meMomenta()[0],mePartonData()[0],incoming); VectorWaveFunction g2w(meMomenta()[1],mePartonData()[1],incoming); SpinorBarWaveFunction qw(meMomenta()[2],mePartonData()[2],outgoing); SpinorWaveFunction qbarw(meMomenta()[3],mePartonData()[3],outgoing); ScalarWaveFunction higgs(meMomenta()[4],mePartonData()[4],1.,outgoing); vector g1,g2; vector q; vector qbar; for(unsigned int ix=0;ix<2;++ix) { g1w.reset(2*ix);g1.push_back(g1w); g2w.reset(2*ix);g2.push_back(g2w); qw.reset(ix);q.push_back(qw); qbarw.reset(ix);qbar.push_back(qbarw); } // calculate the matrix element me=ggME(g1,g2,q,qbar,higgs,0); } // q qbar initiated else { SpinorWaveFunction q1w(meMomenta()[0],mePartonData()[0],incoming); SpinorBarWaveFunction q2w(meMomenta()[1],mePartonData()[1],incoming); SpinorBarWaveFunction q3w(meMomenta()[2],mePartonData()[2],outgoing); SpinorWaveFunction q4w(meMomenta()[3],mePartonData()[3],outgoing); ScalarWaveFunction higgs(meMomenta()[4],mePartonData()[4],1.,outgoing); vector q1,q4; vector q2,q3; for(unsigned int ix=0;ix<2;++ix) { q1w.reset(ix);q1.push_back(q1w); q2w.reset(ix);q2.push_back(q2w); q3w.reset(ix);q3.push_back(q3w); q4w.reset(ix);q4.push_back(q4w); } // calculate the matrix element me = qqME(q1,q2,q3,q4,higgs,0); } return me*sHat()*UnitRemoval::InvE2; } double MEPP2QQHiggs::ggME(vector &g1, vector &g2, vector & q, vector & qbar, ScalarWaveFunction & hwave, unsigned int iflow) const { // scale Energy2 mt(scale()); Energy mass = q[0].mass(); // matrix element to be stored if(iflow!=0) me_.reset(ProductionMatrixElement(PDT::Spin1,PDT::Spin1, PDT::Spin1Half,PDT::Spin1Half, PDT::Spin0)); // calculate the matrix element double output(0.),sumflow[2]={0.,0.}; double sumdiag[8]={0.,0.,0.,0.,0.,0.,0.,0.}; Complex diag[8],flow[2]; VectorWaveFunction interv; SpinorWaveFunction inters,QBoff; SpinorBarWaveFunction intersb,Qoff; for(unsigned int ihel1=0;ihel1<2;++ihel1) { for(unsigned int ihel2=0;ihel2<2;++ihel2) { interv = GGGVertex_->evaluate(mt,5,gluon_,g1[ihel1],g2[ihel2]); for(unsigned int ohel1=0;ohel1<2;++ohel1) { Qoff = QQHVertex_->evaluate(mt,3,q[ohel1].particle(), q[ohel1],hwave,mass); for(unsigned int ohel2=0;ohel2<2;++ohel2) { QBoff = QQHVertex_->evaluate(mt,3,qbar[ohel2].particle(), qbar[ohel2],hwave,mass); // 1st diagram inters = QQGVertex_->evaluate(mt,1,qbar[ohel2].particle(), qbar[ohel2],g2[ihel2],mass); diag[0] = QQGVertex_->evaluate(mt,inters,Qoff,g1[ihel1]); // 2nd diagram intersb = QQGVertex_->evaluate(mt,1,q[ohel1].particle(), q[ohel1],g1[ihel1],mass); diag[1] = QQHVertex_->evaluate(mt,inters,intersb,hwave); // 3rd diagram diag[2] = QQGVertex_->evaluate(mt,QBoff,intersb,g2[ihel2]); // 4th diagram inters = QQGVertex_->evaluate(mt,1,qbar[ohel2].particle(), qbar[ohel2],g1[ihel1],mass); diag[3] = QQGVertex_->evaluate(mt,inters,Qoff,g2[ihel2]); // 5th diagram intersb = QQGVertex_->evaluate(mt,1,q[ohel1].particle(), q[ohel1],g2[ihel2],mass); diag[4] = QQHVertex_->evaluate(mt,inters,intersb,hwave); // 6th diagram diag[5] = QQGVertex_->evaluate(mt,QBoff,intersb,g1[ihel1]); // 7th diagram diag[6] = QQGVertex_->evaluate(mt,qbar[ohel2],Qoff ,interv); // 8th diagram diag[7] = QQGVertex_->evaluate(mt,QBoff ,q[ohel1],interv); // colour flows flow[0]=diag[0]+diag[1]+diag[2]+(diag[6]+diag[7]); flow[1]=diag[3]+diag[4]+diag[5]-(diag[6]+diag[7]); // sums for(unsigned int ix=0;ix<8;++ix) sumdiag[ix] += norm(diag[ix]); for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]); // total output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1]) -0.25*flow[0]*conj(flow[1])); // store the me if needed if(iflow!=0) me_(2*ihel1,2*ihel2,ohel1,ohel2,0)=flow[iflow-1]; } } } } // select a colour flow flow_ = 1 + UseRandom::rnd2(sumflow[0],sumflow[1]); if(flow_==1) sumdiag[0]=sumdiag[1]=sumdiag[2]=0.; else sumdiag[3]=sumdiag[4]=sumdiag[5]=0.; // select a diagram from that flow double prob = UseRandom::rnd(); for(unsigned int ix=0;ix<8;++ix) { if(prob<=sumdiag[ix]) { diagram_=1+ix; break; } prob -= sumdiag[ix]; } // final part of colour and spin factors return output/48.; } double MEPP2QQHiggs::qqME(vector & q1, vector & q2, vector & q3, vector & q4, ScalarWaveFunction & hwave, unsigned int iflow) const { // scale Energy2 mt(scale()); Energy mass = q3[0].mass(); // matrix element to be stored if(iflow!=0) me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1Half,PDT::Spin1Half, PDT::Spin0)); // calculate the matrix element double output(0.),sumdiag[2]={0.,0.}; Complex diag[2]; VectorWaveFunction interv; SpinorWaveFunction QBoff; SpinorBarWaveFunction Qoff; for(unsigned int ihel1=0;ihel1<2;++ihel1) { for(unsigned int ihel2=0;ihel2<2;++ihel2) { interv = QQGVertex_->evaluate(mt,5,gluon_,q1[ihel1],q2[ihel2]); for(unsigned int ohel1=0;ohel1<2;++ohel1) { Qoff = QQHVertex_->evaluate(mt,3,q3[ohel1].particle(), q3[ohel1],hwave,mass); for(unsigned int ohel2=0;ohel2<2;++ohel2) { QBoff = QQHVertex_->evaluate(mt,3,q4[ohel2].particle(), q4[ohel2],hwave,mass); // 1st diagram diag[0] = QQGVertex_->evaluate(mt,q4[ohel2],Qoff,interv); // 2nd diagram diag[1] = QQGVertex_->evaluate(mt,QBoff,q3[ohel1],interv); // sum of diagrams for(unsigned int ix=0;ix<2;++ix) sumdiag[ix] += norm(diag[ix]); diag[0] += diag[1]; output += norm(diag[0]); if(iflow!=0) me_(ihel1,ihel2,ohel1,ohel2,0) = diag[0]; } } } } // only 1 colour flow flow_=1; // select a diagram diagram_ = 9+UseRandom::rnd2(sumdiag[0],sumdiag[1]); // final part of colour and spin factors return output/18.; } Selector MEPP2QQHiggs::colourGeometries(tcDiagPtr diag) const { // colour lines for gg -> Q Qbar H static const ColourLines cgg[10]= {ColourLines("1 4 5, -1 -2 3 , -3 -6 "), ColourLines("1 5 , -1 -2 -3 4, -4 -6 "), ColourLines("1 4 , -1 -2 3 , -3 -5 -6"), ColourLines("3 4 5, 1 2 -3 , -1 -6 "), ColourLines("4 5 , 1 2 3 -4, -1 -6"), ColourLines("3 4 , 1 2 -3 , -1 -5 -6"), ColourLines("1 3 4 5, -1 2, -2 -3 -6"), ColourLines("2 3 4 5, 1 -2, -1 -3 -6"), ColourLines("1 3 4, -1 2, -2 -3 -5 -6"), ColourLines("2 3 4, 1 -2, -1 -3 -5 -6")}; // colour lines for q qbar -> Q Qbar H static const ColourLines cqq[2]= {ColourLines("1 3 4 5, -2 -3 -6"), ColourLines("1 3 4 , -2 -3 -5 -6")}; // select the colour flow (as all ready picked just insert answer) Selector sel; switch(abs(diag->id())) { // gg -> q qbar subprocess case 1: case 2: case 3: case 4: case 5: case 6: sel.insert(1.0, &cgg[abs(diag->id())-1]); break; case 7: sel.insert(1.0, &cgg[5 + flow_]); break; case 8: sel.insert(1.0, &cgg[7 + flow_]); break; // q qbar -> q qbar subprocess case 9: case 10: sel.insert(1.0, &cqq[abs(diag->id())-9]); break; } return sel; } Selector MEPP2QQHiggs::diagrams(const DiagramVector & diags) const { // select the diagram, this is easy for us as we have already done it Selector sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ) { if(diags[i]->id()==-int(diagram_)) sel.insert(1.0, i); else sel.insert(0., i); } return sel; } void MEPP2QQHiggs::constructVertex(tSubProPtr sub) { // extract the particles in the hard process ParticleVector hard; hard.push_back(sub->incoming().first); hard.push_back(sub->incoming().second); for(unsigned int ix=0;ix<3;++ix) hard.push_back(sub->outgoing()[ix]); // identify the process and calculate the matrix element if(hard[0]->id()<0) swap(hard[0],hard[1]); if(hard[2]->id()==ParticleID::h0) swap(hard[2],hard[4]); if(hard[3]->id()==ParticleID::h0) swap(hard[3],hard[4]); if(hard[2]->id()<0) swap(hard[2],hard[3]); if(hard[0]->id()==ParticleID::g) { vector g1,g2; vector q; vector qbar; // off-shell wavefunctions for the spin correlations VectorWaveFunction( g1,hard[0],incoming,false,true,true); VectorWaveFunction( g2,hard[1],incoming,false,true,true); SpinorBarWaveFunction(q ,hard[2],outgoing,true ,true); SpinorWaveFunction( qbar,hard[3],outgoing,true ,true); ScalarWaveFunction hwave( hard[4],outgoing,true); g1[1]=g1[2];g2[1]=g2[2]; ggME(g1,g2,q,qbar,hwave,flow_); } // q qbar -> Q Qbar Higgs else { vector q1,q4; vector q2,q3; // off-shell for spin correlations SpinorWaveFunction( q1,hard[0],incoming,false,true); SpinorBarWaveFunction(q2,hard[1],incoming,false,true); SpinorBarWaveFunction(q3,hard[2],outgoing,true ,true); SpinorWaveFunction( q4,hard[3],outgoing,true ,true); ScalarWaveFunction hwave( hard[4],outgoing,true); qqME(q1,q2,q3,q4,hwave,flow_); } // construct the vertex HardVertexPtr hardvertex=new_ptr(HardVertex()); // set the matrix element for the vertex hardvertex->ME(me_); // set the pointers and to and from the vertex for(unsigned int ix=0;ix<5;++ix) hard[ix]->spinInfo()->productionVertex(hardvertex); } diff --git a/MatrixElement/Hadron/MEPP2WJet.cc b/MatrixElement/Hadron/MEPP2WJet.cc --- a/MatrixElement/Hadron/MEPP2WJet.cc +++ b/MatrixElement/Hadron/MEPP2WJet.cc @@ -1,839 +1,839 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEPP2WJet class. // #include "MEPP2WJet.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/MatrixElement/Tree2toNDiagram.h" #include "Herwig/Models/StandardModel/StandardModel.h" #include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "ThePEG/Cuts/Cuts.h" #include "Herwig/MatrixElement/HardVertex.h" using namespace Herwig; MEPP2WJet::MEPP2WJet() : _process(0), _maxflavour(5), _plusminus(0), _wdec(0), _widthopt(1) {} void MEPP2WJet::doinit() { HwMEBase::doinit(); _wplus = getParticleData(ThePEG::ParticleID::Wplus ); _wminus = getParticleData(ThePEG::ParticleID::Wminus); // cast the SM pointer to the Herwig SM pointer ThePEG::Ptr::transient_const_pointer hwsm=ThePEG::dynamic_ptr_cast< ThePEG::Ptr ::transient_const_pointer>(standardModel()); // do the initialisation if(!hwsm) throw InitException() << "Must be Herwig::StandardModel in MEPP2WJet::doinit()" << Exception::runerror; // set the vertex pointers _theFFWVertex = hwsm->vertexFFW(); _theQQGVertex = hwsm->vertexFFG(); } ClassDescription MEPP2WJet::initMEPP2WJet; // Definition of the static class description member. void MEPP2WJet::Init() { static ClassDocumentation documentation ("The MEPP2WJet class implements the matrix element for W + jet production"); static Parameter interfaceMaxFlavour ( "MaxFlavour", "The heaviest incoming quark flavour this matrix element is allowed to handle " "(if applicable).", &MEPP2WJet::_maxflavour, 5, 0, 8, false, false, true); static Switch interfaceProcess ("Process", "Which subprocesses to include", &MEPP2WJet::_process, 0, false, false); static SwitchOption interfaceProcessAll (interfaceProcess, "All", "Include all subprocesses", 0); static SwitchOption interfaceProcessqqbar (interfaceProcess, "qqbar", "Only include q qbar -> W g process", 1); static SwitchOption interfaceProcessqg (interfaceProcess, "qg", "Only include the q g -> W q process", 2); static SwitchOption interfaceProcessqbarg (interfaceProcess, "qbarg", "Only include the qbar g -> W qbar process", 3); static Switch interfacePlusMinus ("Wcharge", "Which intermediate W bosons to include", &MEPP2WJet::_plusminus, 0, false, false); static SwitchOption interfacePlusMinusAll (interfacePlusMinus, "Both", "Include W+ and W-", 0); static SwitchOption interfacePlusMinusPlus (interfacePlusMinus, "Plus", "Only include W+", 1); static SwitchOption interfacePlusMinusMinus (interfacePlusMinus, "Minus", "Only include W-", 2); static Switch interfaceWDecay ("WDecay", "Which processes to include", &MEPP2WJet::_wdec, 0, false, false); static SwitchOption interfaceWDecayAll (interfaceWDecay, "All", "Include all SM fermions as outgoing particles", 0); static SwitchOption interfaceWDecayQuarks (interfaceWDecay, "Quarks", "Only include outgoing quarks", 1); static SwitchOption interfaceWDecayLeptons (interfaceWDecay, "Leptons", "All include outgoing leptons", 2); static SwitchOption interfaceWDecayElectron (interfaceWDecay, "Electron", "Only include outgoing e nu_e", 3); static SwitchOption interfaceWDecayMuon (interfaceWDecay, "Muon", "Only include outgoing mu nu_mu", 4); static SwitchOption interfaceWDecayTau (interfaceWDecay, "Tau", "Only include outgoing tauu nu_tau", 5); static SwitchOption interfaceWDecayUpDown (interfaceWDecay, "UpDown", "Only include outgoing u dbar/ d ubar", 6); static SwitchOption interfaceWDecayUpStrange (interfaceWDecay, "UpStrange", "Only include outgoing u sbar/ s ubar", 7); static SwitchOption interfaceWDecayUpBottom (interfaceWDecay, "UpBottom", "Only include outgoing u bbar/ b ubar", 8); static SwitchOption interfaceWDecayCharmDown (interfaceWDecay, "CharmDown", "Only include outgoing c dbar/ d cbar", 9); static SwitchOption interfaceWDecayCharmStrange (interfaceWDecay, "CharmStrange", "Only include outgoing c sbar/ s cbar", 10); static SwitchOption interfaceWDecayCharmBottom (interfaceWDecay, "CharmBottom", "Only include outgoing c bbar/ b cbar", 11); static Switch interfaceWidthOption ("WidthOption", "The option for handling the width of the off-shell W boson", &MEPP2WJet::_widthopt, 1, false, false); static SwitchOption interfaceWidthOptionFixedDenominator (interfaceWidthOption, "FixedDenominator", "Use a fixed with in the W propagator but the full matrix element" " in the numerator", 1); static SwitchOption interfaceWidthOptionAllRunning (interfaceWidthOption, "AllRunning", "Use a running width in the W propagator and the full matrix " "element in the numerator", 2); } void MEPP2WJet::getDiagrams() const { // which intgermediates to include bool wplus = _plusminus==0 || _plusminus==1; bool wminus = _plusminus==0 || _plusminus==2; // possible incoming and outgoing particles typedef std::vector > Pairvector; // possible parents Pairvector parentpair; parentpair.reserve(6); // don't even think of putting 'break' in here! switch(_maxflavour) { case 5: parentpair.push_back(make_pair(ParticleID::b, ParticleID::cbar)); parentpair.push_back(make_pair(ParticleID::b, ParticleID::ubar)); [[fallthrough]]; case 4: parentpair.push_back(make_pair(ParticleID::s, ParticleID::cbar)); parentpair.push_back(make_pair(ParticleID::d, ParticleID::cbar)); [[fallthrough]]; case 3: parentpair.push_back(make_pair(ParticleID::s, ParticleID::ubar)); [[fallthrough]]; case 2: parentpair.push_back(make_pair(ParticleID::d, ParticleID::ubar)); [[fallthrough]]; default: ; } // possible children Pairvector childpair; childpair.reserve(9); childpair.push_back(make_pair(ParticleID::eminus, ParticleID::nu_ebar)); childpair.push_back(make_pair(ParticleID::muminus, ParticleID::nu_mubar)); childpair.push_back(make_pair(ParticleID::tauminus, ParticleID::nu_taubar)); childpair.push_back(make_pair(ParticleID::d, ParticleID::ubar)); childpair.push_back(make_pair(ParticleID::s, ParticleID::ubar)); childpair.push_back(make_pair(ParticleID::b, ParticleID::ubar)); childpair.push_back(make_pair(ParticleID::d, ParticleID::cbar)); childpair.push_back(make_pair(ParticleID::s, ParticleID::cbar)); childpair.push_back(make_pair(ParticleID::b, ParticleID::cbar)); // gluon for diagrams tcPDPtr g = getParticleData(ParticleID::g); // loop over the children bool lepton,quark; Pairvector::const_iterator child = childpair.begin(); for (; child != childpair.end(); ++child) { // allowed leptonic decay lepton=child->first>10&& (_wdec==0||_wdec==2|| (abs(child->first)-5)/2==int(_wdec)); // allowed quark decay quark =abs(child->second)<10&& (_wdec==0||_wdec==1|| (abs(child->second)==2&&(abs(child->first)+11)/2==int(_wdec))|| (abs(child->second)==4&&(abs(child->first)+17)/2==int(_wdec))); // if decay not allowed skip if(!(quark||lepton)) continue; // decay products tcPDPtr lNeg1 = getParticleData(child->first); tcPDPtr lNeg2 = getParticleData(child->second); tcPDPtr lPos1 = lNeg2->CC(); tcPDPtr lPos2 = lNeg1->CC(); Pairvector::const_iterator parent = parentpair.begin(); for (; parent != parentpair.end(); ++parent) { // parents tcPDPtr qNeg1 = getParticleData(parent->first); tcPDPtr qNeg2 = getParticleData(parent->second); tcPDPtr qPos1 = qNeg2->CC(); tcPDPtr qPos2 = qNeg1->CC(); // diagrams // q qbar annhilation processes if(_process==0||_process==1) { // q qbar -> W- g if(wminus) { add(new_ptr((Tree2toNDiagram(3), qNeg1, qNeg2, qNeg2, 1, _wminus, 2, g, 4, lNeg1, 4, lNeg2, -1))); add(new_ptr((Tree2toNDiagram(3), qNeg1, qNeg1, qNeg2, 2, _wminus, 1, g, 4, lNeg1, 4, lNeg2, -2))); } // q qbar -> W+ g if(wplus) { add(new_ptr((Tree2toNDiagram(3), qPos1, qPos2, qPos2, 1, _wplus, 2, g, 4, lPos1, 4, lPos2, -3))); add(new_ptr((Tree2toNDiagram(3), qPos1, qPos1, qPos2, 2, _wplus, 1, g, 4, lPos1, 4, lPos2, -4))); } } // q g compton if(_process==0||_process==2) { if(wminus) { add(new_ptr((Tree2toNDiagram(3), qNeg1, qPos1, g , 1, _wminus, 2, qPos1, 4, lNeg1, 4, lNeg2, -5))); add(new_ptr((Tree2toNDiagram(2), qNeg1, g, 1, qNeg1, 3, _wminus, 3, qPos1, 4, lNeg1, 4, lNeg2, -6))); } if(wplus) { add(new_ptr((Tree2toNDiagram(3), qPos1, qNeg1, g, 1, _wplus, 2, qNeg1, 4, lPos1, 4, lPos2, -7))); add(new_ptr((Tree2toNDiagram(2), qPos1, g, 1, qNeg1, 3, _wplus, 3, qNeg1, 4, lPos1, 4, lPos2, -8))); } } // qbar g compton if(_process==0||_process==3) { if(wminus) { add(new_ptr((Tree2toNDiagram(3), qNeg2, qPos2, g, 1, _wminus, 2, qPos2, 4, lNeg1, 4, lNeg2, -9 ))); add(new_ptr((Tree2toNDiagram(2), qNeg2, g, 1, qNeg2, 3, _wminus, 3, qPos2, 4, lNeg1, 4, lNeg2, -10))); } if(wplus) { add(new_ptr((Tree2toNDiagram(3), qPos2, qNeg2, g, 1, _wplus, 2, qNeg2, 4, lPos1, 4, lPos2, -11))); add(new_ptr((Tree2toNDiagram(2), qPos2, g, 1, qPos2, 3, _wplus, 3, qNeg2, 4, lPos1, 4, lPos2, -12))); } } } } } unsigned int MEPP2WJet::orderInAlphaS() const { return 1; } unsigned int MEPP2WJet::orderInAlphaEW() const { return 2; } void MEPP2WJet::persistentOutput(PersistentOStream & os) const { os << _theFFWVertex << _theQQGVertex << _wplus << _widthopt << _wminus << _process << _maxflavour << _plusminus << _wdec; } void MEPP2WJet::persistentInput(PersistentIStream & is, int) { is >> _theFFWVertex >> _theQQGVertex >> _wplus >> _widthopt >> _wminus >> _process >> _maxflavour >> _plusminus >> _wdec; } int MEPP2WJet::nDim() const { return 5; } Selector MEPP2WJet::colourGeometries(tcDiagPtr diag) const { // colour lines for q qbar -> W g static const ColourLines cqqbar[4]={ColourLines("1 -2 5,-3 -5"), ColourLines("1 5, -5 2 -3"), ColourLines("1 -2 5,-3 -5,6 -7"), ColourLines("1 5, -5 2 -3,6 -7")}; // colour lines for q g -> W q static const ColourLines cqg [4]={ColourLines("1 2 -3,3 5"), ColourLines("1 -2,2 3 5"), ColourLines("1 2 -3,3 5,6 -7"), ColourLines("1 -2,2 3 5,6 -7")}; // colour lines for qbar q -> W qbar static const ColourLines cqbarg[4]={ColourLines("-1 -2 3,-3 -5"), ColourLines("-1 2,-2 -3 -5"), ColourLines("-1 -2 3,-3 -5,6 -7"), ColourLines("-1 2,-2 -3 -5,6 -7")}; // select the correct line unsigned int icol = mePartonData()[3]->coloured() ? 2 : 0; Selector sel; switch(abs(diag->id())) { case 1 : case 3: sel.insert(1.0, &cqqbar[icol]); break; case 2 : case 4: sel.insert(1.0, &cqqbar[icol+1]); break; case 5 : case 7: sel.insert(1.0, &cqg[icol]); break; case 6 : case 8: sel.insert(1.0, &cqg[icol+1]); break; case 9 : case 11: sel.insert(1.0, &cqbarg[icol]); break; case 10 : case 12: sel.insert(1.0, &cqbarg[icol+1]); break; } return sel; } Selector MEPP2WJet::diagrams(const DiagramVector & diags) const { Selector sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ) { int id=abs(diags[i]->id()); if (id <= 2 ) sel.insert(meInfo()[id- 1],i); else if(id <= 4 ) sel.insert(meInfo()[id- 3],i); else if(id <= 6 ) sel.insert(meInfo()[id- 5],i); else if(id <= 8 ) sel.insert(meInfo()[id- 7],i); else if(id <= 10) sel.insert(meInfo()[id- 9],i); else if(id <= 12) sel.insert(meInfo()[id-11],i); } return sel; } Energy2 MEPP2WJet::scale() const { return _scale; } CrossSection MEPP2WJet::dSigHatDR() const { return me2()*jacobian()/(16.0*sqr(Constants::pi)*sHat())*sqr(hbarc); } bool MEPP2WJet::generateKinematics(const double * r) { // initialize jacobian jacobian(1.); // cms energy Energy ecm=sqrt(sHat()); // find the right W pointer tcPDPtr wdata = mePartonData()[3]->iCharge()+mePartonData()[4]->iCharge() > 0 ? _wplus :_wminus; // first generate the mass of the off-shell gauge boson // minimum mass of the tcPDVector ptemp; ptemp.push_back(mePartonData()[3]); ptemp.push_back(mePartonData()[4]); Energy2 minMass2 = max(lastCuts().minSij(mePartonData()[3],mePartonData()[4]), lastCuts().minS(ptemp)); // minimum pt of the jet Energy ptmin = max(lastCuts().minKT(mePartonData()[2]), lastCuts().minKT(wdata)); // maximum mass of the gauge boson so pt is possible Energy2 maxMass2 = min(ecm*(ecm-2.*ptmin),lastCuts().maxS(ptemp)); if(maxMass2<=ZERO||minMass2massMin())); maxMass2 = min(maxMass2,sqr(wdata->massMax())); // return if not kinematically possible if(minMass2>maxMass2) return false; // generation of the mass Energy M(wdata->mass()),Gamma(wdata->width()); Energy2 M2(sqr(M)),MG(M*Gamma); double rhomin = atan2((minMass2-M2),MG); double rhomax = atan2((maxMass2-M2),MG); _mw2=M2+MG*tan(rhomin+r[1]*(rhomax-rhomin)); Energy mw=sqrt(_mw2); // jacobian jacobian(jacobian()*(sqr(_mw2-M2)+sqr(MG))/MG*(rhomax-rhomin)/sHat()); // set the masses of the outgoing particles in the 2-2 scattering meMomenta()[2].setMass(ZERO); Lorentz5Momentum pw(mw); // generate the polar angle of the hard scattering double ctmin(-1.0), ctmax(1.0); Energy q(ZERO); try { q = SimplePhaseSpace::getMagnitude(sHat(), meMomenta()[2].mass(),mw); } - catch ( ImpossibleKinematics ) { + catch ( ImpossibleKinematics & e ) { return false; } Energy2 pq = sqrt(sHat())*q; if ( ptmin > ZERO ) { double ctm = 1.0 - sqr(ptmin/q); if ( ctm <= 0.0 ) return false; ctmin = max(ctmin, -sqrt(ctm)); ctmax = min(ctmax, sqrt(ctm)); } if ( ctmin >= ctmax ) return false; double cth = getCosTheta(ctmin, ctmax, r[0]); // momenta of particle in hard scattering Energy pt = q*sqrt(1.0-sqr(cth)); double phi=2.0*Constants::pi*r[2]; meMomenta()[2].setVect(Momentum3( pt*sin(phi), pt*cos(phi), q*cth)); pw.setVect( Momentum3(-pt*sin(phi),-pt*cos(phi),-q*cth)); meMomenta()[2].rescaleEnergy(); pw.rescaleEnergy(); // set the scale _scale = _mw2+sqr(pt); // generate the momenta of the W decay products meMomenta()[3].setMass(mePartonData()[3]->mass()); meMomenta()[4].setMass(mePartonData()[4]->mass()); Energy q2 = ZERO; try { q2 = SimplePhaseSpace::getMagnitude(_mw2, meMomenta()[3].mass(), meMomenta()[4].mass()); - } catch ( ImpossibleKinematics ) { + } catch ( ImpossibleKinematics & e ) { return false; } double cth2 =-1.+2.*r[3]; double phi2=Constants::twopi*r[4]; Energy pt2 =q2*sqrt(1.-sqr(cth2)); Lorentz5Momentum pl[2]={Lorentz5Momentum( pt2*cos(phi2), pt2*sin(phi2), q2*cth2,ZERO, meMomenta()[3].mass()), Lorentz5Momentum(-pt2*cos(phi2),-pt2*sin(phi2),-q2*cth2,ZERO, meMomenta()[4].mass())}; pl[0].rescaleEnergy(); pl[1].rescaleEnergy(); Boost boostv(pw.boostVector()); pl[0].boost(boostv); pl[1].boost(boostv); meMomenta()[3] = pl[0]; meMomenta()[4] = pl[1]; // check passes all the cuts vector out(3); out[0] = meMomenta()[2]; out[1] = meMomenta()[3]; out[2] = meMomenta()[4]; tcPDVector tout(3); tout[0] = mePartonData()[2]; tout[1] = mePartonData()[3]; tout[2] = mePartonData()[4]; if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) ) return false; // jacobian jacobian((pq/sHat())*Constants::pi*jacobian()/8./sqr(Constants::pi)*q2/mw); return true; } double MEPP2WJet::me2() const { InvEnergy2 output(ZERO); // construct spinors for the leptons (always the same) vector lm; vector lp; SpinorBarWaveFunction lmout(meMomenta()[3],mePartonData()[3],outgoing); SpinorWaveFunction lpout(meMomenta()[4],mePartonData()[4],outgoing); for(unsigned int ix=0;ix<2;++ix) { lmout.reset(ix);lm.push_back(lmout); lpout.reset(ix);lp.push_back(lpout); } // q g to q W if(mePartonData()[0]->id()<=6&&mePartonData()[0]->id()>0&& mePartonData()[1]->id()==ParticleID::g) { // polarization states for the particles vector fin; vector gin; vector fout; SpinorWaveFunction qin (meMomenta()[0],mePartonData()[0],incoming); VectorWaveFunction glin(meMomenta()[1],mePartonData()[1],incoming); SpinorBarWaveFunction qout(meMomenta()[2],mePartonData()[2],outgoing); for(unsigned int ix=0;ix<2;++ix) { qin.reset(ix) ; fin.push_back(qin); glin.reset(2*ix); gin.push_back(glin); qout.reset(ix);fout.push_back(qout); } output=qgME(fin,gin,fout,lm,lp); } // qbar g to qbar W else if(mePartonData()[0]->id()>=-6&&mePartonData()[0]->id()<0&& mePartonData()[1]->id()==ParticleID::g) { vector ain; vector gin; vector aout; SpinorBarWaveFunction qbin (meMomenta()[0],mePartonData()[0],incoming); VectorWaveFunction glin (meMomenta()[1],mePartonData()[1],incoming); SpinorWaveFunction qbout(meMomenta()[2],mePartonData()[2],outgoing); for(unsigned int ix=0;ix<2;++ix) { qbin.reset(ix) ; ain.push_back(qbin); glin.reset(2*ix) ; gin.push_back(glin); qbout.reset(ix);aout.push_back(qbout); } output=qbargME(ain,gin,aout,lm,lp); } // q qbar to g W else { vector fin; vector ain; vector gout; SpinorWaveFunction qin (meMomenta()[0],mePartonData()[0],incoming); SpinorBarWaveFunction qbin(meMomenta()[1],mePartonData()[1],incoming); VectorWaveFunction glout(meMomenta()[2],mePartonData()[2],outgoing); for(unsigned int ix=0;ix<2;++ix) { qin.reset(ix) ; fin.push_back(qin); qbin.reset(ix) ; ain.push_back(qbin); glout.reset(2*ix); gout.push_back(glout); } output=qqbarME(fin,ain,gout,lm,lp); } return output*sHat(); } InvEnergy2 MEPP2WJet::qqbarME(vector & fin, vector & ain, vector & gout, vector & lm, vector & lp, bool calc) const { // if calculation spin corrections construct the me if(calc) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1,PDT::Spin1Half, PDT::Spin1Half)); // some integers unsigned int ihel1,ihel2,ohel1,ohel2,ohel3; // find the right W pointer tcPDPtr wdata = mePartonData()[3]->iCharge()+mePartonData()[4]->iCharge() > 0 ? _wplus :_wminus; // compute the W current for speed VectorWaveFunction bcurr[2][2]; for(ohel2=0;ohel2<2;++ohel2) { for(ohel3=0;ohel3<2;++ohel3) { bcurr[ohel2][ohel3] = _theFFWVertex->evaluate(_mw2,_widthopt,wdata, lp[ohel3],lm[ohel2]); } } double me[3]={0.,0.,0.}; Complex diag[2]; SpinorWaveFunction inters; SpinorBarWaveFunction interb; for(ihel1=0;ihel1<2;++ihel1) { for(ihel2=0;ihel2<2;++ihel2) { for(ohel1=0;ohel1<2;++ohel1) { // intermediates for the diagrams inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[0], fin[ihel1],gout[ohel1]); interb=_theQQGVertex->evaluate(_scale,5,mePartonData()[1], ain[ihel2],gout[ohel1]); for(ohel2=0;ohel2<2;++ohel2) { for(ohel3=0;ohel3<2;++ohel3) { diag[0] = _theFFWVertex->evaluate(_mw2,fin[ihel1],interb, bcurr[ohel2][ohel3]); diag[1] = _theFFWVertex->evaluate(_mw2,inters,ain[ihel2], bcurr[ohel2][ohel3]); // diagram contributions me[1] += norm(diag[0]); me[2] += norm(diag[1]); // total diag[0] += diag[1]; me[0] += norm(diag[0]); if(calc) _me(ihel1,ihel2,2*ohel1,ohel2,ohel3) = diag[0]; } } } } } // results // initial state spin and colour average double colspin=1./9./4.; // and C_F N_c from matrix element colspin *= 4.; // colour factor for the W decay if(mePartonData()[3]->coloured()) colspin*=3.; DVector save; for(unsigned int ix=0;ix<3;++ix) { me[ix]*=colspin; if(ix>0) save.push_back(me[ix]); } meInfo(save); return me[0] * UnitRemoval::InvE2; } InvEnergy2 MEPP2WJet::qgME(vector & fin, vector & gin, vector & fout, vector & lm, vector & lp, bool calc) const { // if calculation spin corrections construct the me if(calc) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1, PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1Half)); // find the right W pointer tcPDPtr wdata = mePartonData()[3]->iCharge()+mePartonData()[4]->iCharge() > 0 ? _wplus :_wminus; // some integers unsigned int ihel1,ihel2,ohel1,ohel2,ohel3; // compute the leptonic W current for speed VectorWaveFunction bcurr[2][2]; for(ohel2=0;ohel2<2;++ohel2) { for(ohel3=0;ohel3<2;++ohel3) { bcurr[ohel2][ohel3] = _theFFWVertex->evaluate(_mw2,_widthopt,wdata, lp[ohel3],lm[ohel2]); } } // compute the matrix elements double me[3]={0.,0.,0.}; Complex diag[2]; SpinorWaveFunction inters; SpinorBarWaveFunction interb; for(ihel1=0;ihel1<2;++ihel1) { for(ihel2=0;ihel2<2;++ihel2) { for(ohel1=0;ohel1<2;++ohel1) { // intermediates for the diagrams interb=_theQQGVertex->evaluate(_scale,5,mePartonData()[2], fout[ohel1],gin[ihel2]); inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[0], fin[ihel1],gin[ihel2]); for(ohel2=0;ohel2<2;++ohel2) { for(ohel3=0;ohel3<2;++ohel3) { diag[0]=_theFFWVertex->evaluate(_mw2,fin[ihel1],interb, bcurr[ohel2][ohel3]); diag[1]=_theFFWVertex->evaluate(_mw2,inters,fout[ohel1], bcurr[ohel2][ohel3]); // diagram contributions me[1] += norm(diag[0]); me[2] += norm(diag[1]); // total diag[0] += diag[1]; me[0] += norm(diag[0]); if(calc) _me(ihel1,2*ihel2,ohel1,ohel2,ohel3) = diag[0]; } } } } } // results // initial state spin and colour average double colspin=1./24./4.; // and C_F N_c from matrix element colspin *=4.; // colour factor for the W decay if(mePartonData()[3]->coloured()) colspin*=3.; DVector save; for(unsigned int ix=0;ix<3;++ix) { me[ix]*=colspin; if(ix>0) save.push_back(me[ix]); } meInfo(save); return me[0] * UnitRemoval::InvE2; } InvEnergy2 MEPP2WJet::qbargME(vector & fin, vector & gin, vector & fout, vector & lm, vector & lp, bool calc) const { // if calculation spin corrections construct the me if(calc) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1, PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1Half)); // find the right W pointer tcPDPtr wdata = mePartonData()[3]->iCharge()+mePartonData()[4]->iCharge() > 0 ? _wplus :_wminus; // some integers unsigned int ihel1,ihel2,ohel1,ohel2,ohel3; // compute the leptonic W current for speed VectorWaveFunction bcurr[2][2]; for(ohel2=0;ohel2<2;++ohel2) { for(ohel3=0;ohel3<2;++ohel3) { bcurr[ohel2][ohel3] = _theFFWVertex->evaluate(_mw2,_widthopt,wdata, lp[ohel3],lm[ohel2]); } } // compute the matrix elements double me[3]={0.,0.,0.}; Complex diag[2]; SpinorWaveFunction inters; SpinorBarWaveFunction interb; for(ihel1=0;ihel1<2;++ihel1) { for(ihel2=0;ihel2<2;++ihel2) { for(ohel1=0;ohel1<2;++ohel1) { // intermediates for the diagrams inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[2], fout[ohel1],gin[ihel2]); interb=_theQQGVertex->evaluate(_scale,5,mePartonData()[0], fin[ihel1],gin[ihel2]); for(ohel2=0;ohel2<2;++ohel2) { for(ohel3=0;ohel3<2;++ohel3) { diag[0]= _theFFWVertex->evaluate(_mw2,inters,fin[ihel1], bcurr[ohel2][ohel3]); diag[1]= _theFFWVertex->evaluate(_mw2,fout[ohel1],interb, bcurr[ohel2][ohel3]); // diagram contributions me[1] += norm(diag[0]); me[2] += norm(diag[1]); // total diag[0] += diag[1]; me[0] += norm(diag[0]); if(calc) _me(ihel1,2*ihel2,ohel1,ohel2,ohel3) = diag[0]; } } } } } // results // initial state spin and colour average double colspin=1./24./4.; // and C_F N_c from matrix element colspin *= 4.; // colour factor for the W decay if(mePartonData()[3]->coloured()) colspin*=3.; DVector save; for(unsigned int ix=0;ix<3;++ix) { me[ix]*=colspin; if(ix>0) save.push_back(me[ix]); } meInfo(save); return me[0] * UnitRemoval::InvE2; } void MEPP2WJet::constructVertex(tSubProPtr sub) { // extract the particles in the hard process ParticleVector hard(5); // incoming hard[0]=sub->incoming().first; hard[1]=sub->incoming().second; if((hard[0]->id()<0&&hard[1]->id()<=6)|| hard[0]->id()==ParticleID::g) swap(hard[0],hard[1]); // outgoing for(unsigned int ix=0;ix<3;++ix) { unsigned int iloc; PPtr mother=sub->outgoing()[ix]->parents()[0]; if(mother&&abs(mother->id())==ParticleID::Wplus) { if(sub->outgoing()[ix]->id()>0) iloc=3; else iloc=4; } else iloc=2; hard[iloc]=sub->outgoing()[ix]; } // wavefunctions for the W decay products vector lm; vector lp; SpinorBarWaveFunction(lm,hard[3],outgoing,true,true); SpinorWaveFunction (lp,hard[4],outgoing,true,true); // identify hard process and calculate matrix element // q g to q W if(hard[0]->id()<=6&&hard[0]->id()>0&&hard[1]->id()==ParticleID::g) { vector fin; vector gin; vector fout; SpinorWaveFunction (fin ,hard[0],incoming,false,true); VectorWaveFunction (gin ,hard[1],incoming,false,true,true); SpinorBarWaveFunction (fout,hard[2],outgoing,true ,true); gin[1]=gin[2]; qgME(fin,gin,fout,lm,lp,true); } // qbar g to qbar W else if(hard[0]->id()>=-6&&hard[0]->id()<0&&hard[1]->id()==ParticleID::g) { vector ain; vector gin; vector aout; SpinorBarWaveFunction(ain ,hard[0],incoming,false,true); VectorWaveFunction (gin ,hard[1],incoming,false,true,true); SpinorWaveFunction (aout,hard[2],outgoing,true ,true); gin[1]=gin[2]; qbargME(ain,gin,aout,lm,lp,true); } // q qbar to g W else { vector fin; vector ain; vector gout; SpinorWaveFunction (fin ,hard[0],incoming,false,true); SpinorBarWaveFunction(ain ,hard[1],incoming,false,true); VectorWaveFunction (gout,hard[2],outgoing,true ,true,true); gout[1]=gout[2]; qqbarME(fin,ain,gout,lm,lp,true); } // construct the vertex HardVertexPtr hardvertex=new_ptr(HardVertex()); // set the matrix element for the vertex hardvertex->ME(_me); // set the pointers and to and from the vertex for(unsigned int ix=0;ix<5;++ix) (hard[ix]->spinInfo())->productionVertex(hardvertex); } diff --git a/MatrixElement/Hadron/MEPP2ZJet.cc b/MatrixElement/Hadron/MEPP2ZJet.cc --- a/MatrixElement/Hadron/MEPP2ZJet.cc +++ b/MatrixElement/Hadron/MEPP2ZJet.cc @@ -1,858 +1,858 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEPP2ZJet class. // #include "MEPP2ZJet.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/MatrixElement/Tree2toNDiagram.h" #include "Herwig/Models/StandardModel/StandardModel.h" #include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "ThePEG/Cuts/Cuts.h" #include "Herwig/MatrixElement/HardVertex.h" using namespace Herwig; MEPP2ZJet::MEPP2ZJet() : _process(0), _maxflavour(5), _zdec(0), _gammaZ(0), _widthopt(1), _pprob(0.5) {} void MEPP2ZJet::doinit() { HwMEBase::doinit(); _z0 = getParticleData(ThePEG::ParticleID::Z0 ); _gamma = getParticleData(ThePEG::ParticleID::gamma); // cast the SM pointer to the Herwig SM pointer ThePEG::Ptr::transient_const_pointer hwsm=ThePEG::dynamic_ptr_cast< ThePEG::Ptr ::transient_const_pointer>(standardModel()); // do the initialisation if(!hwsm) throw InitException() << "Must be Herwig::StandardModel in MEPP2ZJet::doinit()" << Exception::runerror; // set the vertex pointers _theFFZVertex = hwsm->vertexFFZ(); _theFFPVertex = hwsm->vertexFFP(); _theQQGVertex = hwsm->vertexFFG(); } ClassDescription MEPP2ZJet::initMEPP2ZJet; // Definition of the static class description member. void MEPP2ZJet::Init() { static ClassDocumentation documentation ("The MEPP2ZJet class implements the matrix element for Z/gamma+ jet production"); static Parameter interfaceMaxFlavour ( "MaxFlavour", "The heaviest incoming quark flavour this matrix element is allowed to handle " "(if applicable).", &MEPP2ZJet::_maxflavour, 5, 0, 8, false, false, true); static Switch interfaceZDecay ("ZDecay", "Which process to included", &MEPP2ZJet::_zdec, 0, false, false); static SwitchOption interfaceZDecayAll (interfaceZDecay, "All", "Include all SM fermions as outgoing particles", 0); static SwitchOption interfaceZDecayQuarks (interfaceZDecay, "Quarks", "All include the quarks as outgoing particles", 1); static SwitchOption interfaceZDecayLeptons (interfaceZDecay, "Leptons", "Only include the leptons as outgoing particles", 2); static SwitchOption interfaceZDecayChargedLeptons (interfaceZDecay, "ChargedLeptons", "Only include the charged leptons as outgoing particles", 3); static SwitchOption interfaceZDecayNeutrinos (interfaceZDecay, "Neutrinos", "Only include the neutrinos as outgoing particles", 4); static SwitchOption interfaceZDecayElectron (interfaceZDecay, "Electron", "Only include e+e- as outgoing particles", 5); static SwitchOption interfaceZDecayMuon (interfaceZDecay, "Muon", "Only include mu+mu- as outgoing particles", 6); static SwitchOption interfaceZDecayTau (interfaceZDecay, "Tau", "Only include tau+tau- as outgoing particles", 7); static SwitchOption interfaceZDecayNu_e (interfaceZDecay, "Nu_e", "Only include nu_e ne_ebar as outgoing particles", 8); static SwitchOption interfaceZDecaynu_mu (interfaceZDecay, "Nu_mu", "Only include nu_mu nu_mubar as outgoing particles", 9); static SwitchOption interfaceZDecaynu_tau (interfaceZDecay, "Nu_tau", "Only include nu_tau nu_taubar as outgoing particles", 10); static SwitchOption interfaceZDecayDown (interfaceZDecay, "Down", "Only include d dbar as outgoing particles", 11); static SwitchOption interfaceZDecayUp (interfaceZDecay, "Up", "Only include u ubar as outgoing particles", 12); static SwitchOption interfaceZDecayStrange (interfaceZDecay, "Strange", "Only include s sbar as outgoing particles", 13); static SwitchOption interfaceZDecayCharm (interfaceZDecay, "Charm", "Only include c cbar as outgoing particles", 14); static SwitchOption interfaceZDecayBottom (interfaceZDecay, "Bottom", "Only include b bbar as outgoing particles", 15); static SwitchOption interfaceZDecayTop (interfaceZDecay, "Top", "Only include t tbar as outgoing particles", 16); static Switch interfaceProcess ("Process", "Which subprocesses to include", &MEPP2ZJet::_process, 0, false, false); static SwitchOption interfaceProcessAll (interfaceProcess, "All", "Include all subprocesses", 0); static SwitchOption interfaceProcessqqbar (interfaceProcess, "qqbar", "Only include q qbar -> Z/gamma g process", 1); static SwitchOption interfaceProcessqg (interfaceProcess, "qg", "Only include the q g -> Z/gamma q process", 2); static SwitchOption interfaceProcessqbarg (interfaceProcess, "qbarg", "Only include the qbar g -> Z/gamma qbar process", 3); static Parameter interfacePhotonProbablity ("PhotonProbablity", "Probability for using the \\f$1/s^2\\f$ piece for the" " generation of the gauge boson mass", &MEPP2ZJet::_pprob, 0.5, 0.0, 1.0, false, false, Interface::limited); static Switch interfaceGammaZ ("GammaZ", "Which terms to include", &MEPP2ZJet::_gammaZ, 0, false, false); static SwitchOption interfaceGammaZAll (interfaceGammaZ, "All", "Include both gamma and Z terms", 0); static SwitchOption interfaceGammaZGamma (interfaceGammaZ, "Gamma", "Only include the photon", 1); static SwitchOption interfaceGammaZZ (interfaceGammaZ, "Z", "Only include the Z", 2); static Switch interfaceWidthOption ("WidthOption", "The option for handling the width of the off-shell W boson", &MEPP2ZJet::_widthopt, 1, false, false); static SwitchOption interfaceWidthOptionFixedDenominator (interfaceWidthOption, "FixedDenominator", "Use a fxied with in the W propagator but the full matrix element" " in the numerator", 1); static SwitchOption interfaceWidthOptionAllRunning (interfaceWidthOption, "AllRunning", "Use a running width in the W propagator and the full matrix " "element in the numerator", 2); } void MEPP2ZJet::getDiagrams() const { // which intermediates to include bool gamma = _gammaZ==0 || _gammaZ==1; bool Z0 = _gammaZ==0 || _gammaZ==2; // pointer for gluon tcPDPtr g = getParticleData(ParticleID::g); bool quark,lepton; for ( int ix=1; ix<17; ++ix ) { // increment counter to switch between quarks and leptons if(ix==7) ix+=4; // is it a valid quark process quark=ix<=6&&(_zdec==0||_zdec==1||_zdec-10==ix); // is it a valid lepton process lepton=ix>=11&&ix<=16&& (_zdec==0||_zdec==2||(_zdec==3&&ix%2==1)|| (_zdec==4&&ix%2==0)||(ix%2==0&&(ix-10)/2==_zdec-7)|| (ix%2==1&&(ix-9)/2==_zdec-4)); // if not a validf process continue if(!(quark||lepton)) continue; // pointer for Z decay products tcPDPtr lm = getParticleData(ix); tcPDPtr lp = lm->CC(); for (int i = 1; i <= _maxflavour; ++i ) { tcPDPtr q = getParticleData(i); tcPDPtr qb = q->CC(); // q qbar -> Z g -> l+l- g if(_process==0||_process==1) { if(gamma) add(new_ptr((Tree2toNDiagram(3), q, q, qb, 1, _gamma, 2, g, 4, lm, 4, lp, -1))); if(Z0) add(new_ptr((Tree2toNDiagram(3), q, q, qb, 1, _z0, 2, g, 4, lm, 4, lp, -2))); if(gamma) add(new_ptr((Tree2toNDiagram(3), q, q, qb, 2, _gamma, 1, g, 4, lm, 4, lp, -3))); if(Z0) add(new_ptr((Tree2toNDiagram(3), q, q, qb, 2, _z0, 1, g, 4, lm, 4, lp, -4))); } // q g -> Z q -> l+l- qbar if(_process==0||_process==2) { if(gamma) add(new_ptr((Tree2toNDiagram(3), q, q, g, 1, _gamma, 2, q, 4, lm, 4, lp, -5))); if(Z0) add(new_ptr((Tree2toNDiagram(3), q, q, g, 1, _z0, 2, q, 4, lm, 4, lp, -6))); if(gamma) add(new_ptr((Tree2toNDiagram(2), q, g, 1, q, 3, _gamma, 3, q, 4, lm, 4, lp, -7))); if(Z0) add(new_ptr((Tree2toNDiagram(2), q, g, 1, q, 3, _z0, 3, q, 4, lm, 4, lp, -8))); } // qbar g -> Z qbar -> l+l- qbar if(_process==0||_process==3) { if(gamma) add(new_ptr((Tree2toNDiagram(3), qb, qb, g, 1, _gamma, 2, qb, 4, lm, 4, lp, -9 ))); if(Z0) add(new_ptr((Tree2toNDiagram(3), qb, qb, g, 1, _z0, 2, qb, 4, lm, 4, lp, -10))); if(gamma) add(new_ptr((Tree2toNDiagram(2), qb, g, 1, qb, 3, _gamma, 3, qb, 4, lm, 4, lp, -11))); if(Z0) add(new_ptr((Tree2toNDiagram(2), qb, g, 1, qb, 3, _z0, 3, qb, 4, lm, 4, lp, -12))); } } } } unsigned int MEPP2ZJet::orderInAlphaS() const { return 1; } unsigned int MEPP2ZJet::orderInAlphaEW() const { return 2; } void MEPP2ZJet::persistentOutput(PersistentOStream & os) const { os << _theFFZVertex << _theFFPVertex << _theQQGVertex << _z0 << _widthopt << _gamma << _process << _maxflavour << _zdec << _pprob << _gammaZ; } void MEPP2ZJet::persistentInput(PersistentIStream & is, int) { is >> _theFFZVertex >> _theFFPVertex >> _theQQGVertex >> _z0 >> _widthopt >> _gamma >> _process >> _maxflavour >> _zdec >> _pprob >> _gammaZ; } int MEPP2ZJet::nDim() const { return 5; } Selector MEPP2ZJet::colourGeometries(tcDiagPtr diag) const { // colour lines for q qbar -> Z g static const ColourLines cqqbar[4]={ColourLines("1 2 5,-3 -5"), ColourLines("1 5,-5 2 -3"), ColourLines("1 2 5,-3 -5, 6 -7"), ColourLines("1 5,-5 2 -3, 6 -7")}; // colour lines for q g -> Z q static const ColourLines cqg [4]={ColourLines("1 2 -3,3 5"), ColourLines("1 -2,2 3 5"), ColourLines("1 2 -3,3 5, 6 -7"), ColourLines("1 -2,2 3 5, 6 -7")}; // colour lines for qbar q -> Z qbar static const ColourLines cqbarg[4]={ColourLines("-1 -2 3,-3 -5"), ColourLines("-1 2,-2 -3 -5"), ColourLines("-1 -2 3,-3 -5, 6 -7"), ColourLines("-1 2,-2 -3 -5, 6 -7")}; // select the correct line unsigned int icol = mePartonData()[3]->coloured() ? 2 : 0; Selector sel; switch(abs(diag->id())) { case 1 : case 2: sel.insert(1.0, &cqqbar[icol]); break; case 3 : case 4: sel.insert(1.0, &cqqbar[icol+1]); break; case 5 : case 6: sel.insert(1.0, &cqg[icol]); break; case 7 : case 8: sel.insert(1.0, &cqg[icol+1]); break; case 9 : case 10: sel.insert(1.0, &cqbarg[icol]); break; case 11 : case 12: sel.insert(1.0, &cqbarg[icol+1]); break; } return sel; } Selector MEPP2ZJet::diagrams(const DiagramVector & diags) const { Selector sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ) { int id=abs(diags[i]->id()); if (id <= 4 ) sel.insert(meInfo()[id-1],i); else if(id <= 8 ) sel.insert(meInfo()[id-5],i); else if(id <= 12) sel.insert(meInfo()[id-9],i); } return sel; } Energy2 MEPP2ZJet::scale() const { return _scale; } CrossSection MEPP2ZJet::dSigHatDR() const { return me2()*jacobian()/(16.0*sqr(Constants::pi)*sHat())*sqr(hbarc); } bool MEPP2ZJet::generateKinematics(const double * r) { // initialize jacobian jacobian(1.); // cms energy Energy ecm=sqrt(sHat()); // first generate the mass of the off-shell gauge boson // minimum mass of the tcPDVector ptemp; ptemp.push_back(mePartonData()[3]); ptemp.push_back(mePartonData()[4]); Energy2 minMass2 = max(lastCuts().minSij(mePartonData()[3],mePartonData()[4]), lastCuts().minS(ptemp)); // minimum pt of the jet Energy ptmin = max(lastCuts().minKT(mePartonData()[2]), lastCuts().minKT(_z0)); // maximum mass of the gauge boson so pt is possible Energy2 maxMass2 = min(ecm*(ecm-2.*ptmin),lastCuts().maxS(ptemp)); if(maxMass2<=ZERO||minMass2massMin())); maxMass2 = min(maxMass2,sqr(_z0->massMax())); // also impose the limits from the ParticleData object if(maxMass2mass()),Gamma(_z0->width()); Energy2 M2(sqr(M)),MG(M*Gamma); double rhomin = atan2((minMass2-M2),MG); double rhomax = atan2((maxMass2-M2),MG); if(r[1]<_pprob) { double rand=r[1]/_pprob; _mz2=minMass2*maxMass2/(minMass2+rand*(maxMass2-minMass2)); } else { double rand=(r[1]-_pprob)/(1.-_pprob); _mz2=M2+MG*tan(rhomin+rand*(rhomax-rhomin)); } Energy mz=sqrt(_mz2); InvEnergy2 emjac1 = _pprob*minMass2*maxMass2/(maxMass2-minMass2)/sqr(_mz2); InvEnergy2 emjac2 = (1.-_pprob)*MG/(rhomax-rhomin)/(sqr(_mz2-M2)+sqr(MG)); // jacobian jacobian(jacobian()/sHat()/(emjac1+emjac2)); // set the masses of the outgoing particles to 2-2 scattering meMomenta()[2].setMass(ZERO); Lorentz5Momentum pz(mz); // generate the polar angle of the hard scattering double ctmin(-1.0), ctmax(1.0); Energy q(ZERO); try { q = SimplePhaseSpace::getMagnitude(sHat(), meMomenta()[2].mass(),mz); } - catch ( ImpossibleKinematics ) { + catch ( ImpossibleKinematics & e ) { return false; } Energy2 pq = sqrt(sHat())*q; if ( ptmin > ZERO ) { double ctm = 1.0 - sqr(ptmin/q); if ( ctm <= 0.0 ) return false; ctmin = max(ctmin, -sqrt(ctm)); ctmax = min(ctmax, sqrt(ctm)); } if ( ctmin >= ctmax ) return false; double cth = getCosTheta(ctmin, ctmax, r[0]); Energy pt = q*sqrt(1.0-sqr(cth)); double phi = 2.0*Constants::pi*r[2]; meMomenta()[2].setVect(Momentum3( pt*sin(phi), pt*cos(phi), q*cth)); pz.setVect( Momentum3(-pt*sin(phi),-pt*cos(phi),-q*cth)); meMomenta()[2].rescaleEnergy(); pz.rescaleEnergy(); // set the scale _scale = _mz2+sqr(pt); // generate the momenta of the Z decay products meMomenta()[3].setMass(mePartonData()[3]->mass()); meMomenta()[4].setMass(mePartonData()[4]->mass()); Energy q2 = ZERO; try { q2 = SimplePhaseSpace::getMagnitude(_mz2, meMomenta()[3].mass(), meMomenta()[4].mass()); - } catch ( ImpossibleKinematics ) { + } catch ( ImpossibleKinematics & e ) { return false; } double cth2 =-1.+2.*r[3]; double phi2=Constants::twopi*r[4]; Energy pt2 =q2*sqrt(1.-sqr(cth2)); Lorentz5Momentum pl[2]={Lorentz5Momentum( pt2*cos(phi2), pt2*sin(phi2), q2*cth2,ZERO, meMomenta()[3].mass()), Lorentz5Momentum(-pt2*cos(phi2),-pt2*sin(phi2),-q2*cth2,ZERO, meMomenta()[4].mass())}; pl[0].rescaleEnergy(); pl[1].rescaleEnergy(); Boost boostv(pz.boostVector()); pl[0].boost(boostv); pl[1].boost(boostv); meMomenta()[3] = pl[0]; meMomenta()[4] = pl[1]; // check passes all the cuts vector out(3); tcPDVector tout(3); for(unsigned int ix=0;ix<3;++ix) { out[ ix] = meMomenta()[ix+2]; tout[ix] = mePartonData()[ix+2]; } if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) ) return false; // jacobian jacobian((pq/sHat())*Constants::pi*jacobian()/8./sqr(Constants::pi)*q2/mz); return true; } double MEPP2ZJet::me2() const { InvEnergy2 output(ZERO); // construct spinors for the leptons (always the same) vector lm; vector lp; SpinorBarWaveFunction lmout(meMomenta()[3],mePartonData()[3],outgoing); SpinorWaveFunction lpout(meMomenta()[4],mePartonData()[4],outgoing); for(unsigned int ix=0;ix<2;++ix) { lmout.reset(ix);lm.push_back(lmout); lpout.reset(ix);lp.push_back(lpout); } // q g to q Z if(mePartonData()[0]->id()<=6&&mePartonData()[0]->id()>0&& mePartonData()[1]->id()==ParticleID::g) { // polarization states for the particles vector fin; vector gin; vector fout; SpinorWaveFunction qin (meMomenta()[0],mePartonData()[0],incoming); VectorWaveFunction glin(meMomenta()[1],mePartonData()[1],incoming); SpinorBarWaveFunction qout(meMomenta()[2],mePartonData()[2],outgoing); for(unsigned int ix=0;ix<2;++ix) { qin.reset(ix) ; fin.push_back(qin); glin.reset(2*ix); gin.push_back(glin); qout.reset(ix);fout.push_back(qout); } output=qgME(fin,gin,fout,lm,lp); } // qbar g to qbar Z else if(mePartonData()[0]->id()>=-6&&mePartonData()[0]->id()<0&& mePartonData()[1]->id()==ParticleID::g) { vector ain; vector gin; vector aout; SpinorBarWaveFunction qbin (meMomenta()[0],mePartonData()[0],incoming); VectorWaveFunction glin (meMomenta()[1],mePartonData()[1],incoming); SpinorWaveFunction qbout(meMomenta()[2],mePartonData()[2],outgoing); for(unsigned int ix=0;ix<2;++ix) { qbin .reset(ix ); ain .push_back(qbin ); glin .reset(2*ix); gin .push_back(glin ); qbout.reset(ix ); aout.push_back(qbout); } output=qbargME(ain,gin,aout,lm,lp); } // q qbar to g Z else { vector fin; vector ain; vector gout; SpinorWaveFunction qin (meMomenta()[0],mePartonData()[0],incoming); SpinorBarWaveFunction qbin(meMomenta()[1],mePartonData()[1],incoming); VectorWaveFunction glout(meMomenta()[2],mePartonData()[2],outgoing); for(unsigned int ix=0;ix<2;++ix) { qin.reset(ix) ; fin.push_back(qin); qbin.reset(ix) ; ain.push_back(qbin); glout.reset(2*ix); gout.push_back(glout); } output=qqbarME(fin,ain,gout,lm,lp); } return output*sHat(); } InvEnergy2 MEPP2ZJet::qqbarME(vector & fin, vector & ain, vector & gout, vector & lm, vector & lp, bool calc) const { // if calculation spin corrections construct the me if(calc) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1,PDT::Spin1Half, PDT::Spin1Half)); // diagrams to include bool gamma = _gammaZ==0 || _gammaZ==1; bool Z0 = _gammaZ==0 || _gammaZ==2; // some integers unsigned int ihel1,ihel2,ohel1,ohel2,ohel3; // compute the leptonic photon and Z currents for speed VectorWaveFunction bcurr[2][2][2]; for(ohel2=0;ohel2<2;++ohel2) { for(ohel3=0;ohel3<2;++ohel3) { // photon current if(gamma) bcurr[0][ohel2][ohel3]= _theFFPVertex->evaluate(_mz2,1,_gamma,lp[ohel3],lm[ohel2]); // Z current if(Z0) bcurr[1][ohel2][ohel3]= _theFFZVertex->evaluate(_mz2,_widthopt,_z0,lp[ohel3],lm[ohel2]); } } // compute the matrix elements double me[5]={0.,0.,0.,0.,0.}; Complex diag[4]; SpinorWaveFunction inters; SpinorBarWaveFunction interb; for(ihel1=0;ihel1<2;++ihel1) { for(ihel2=0;ihel2<2;++ihel2) { for(ohel1=0;ohel1<2;++ohel1) { // intermediates for the diagrams inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[0], fin[ihel1],gout[ohel1]); interb=_theQQGVertex->evaluate(_scale,5,mePartonData()[1], ain[ihel2],gout[ohel1]); for(ohel2=0;ohel2<2;++ohel2) { for(ohel3=0;ohel3<2;++ohel3) { diag[0] = gamma ? _theFFPVertex->evaluate(_mz2,fin[ihel1],interb, bcurr[0][ohel2][ohel3]) : 0.; diag[1]= Z0 ? _theFFZVertex->evaluate(_mz2,fin[ihel1],interb, bcurr[1][ohel2][ohel3]) : 0.; diag[2]= gamma ? _theFFPVertex->evaluate(_mz2,inters,ain[ihel2], bcurr[0][ohel2][ohel3]) : 0.; diag[3]= Z0 ? _theFFZVertex->evaluate(_mz2,inters,ain[ihel2], bcurr[1][ohel2][ohel3]) : 0.; // diagram contributions me[1] += norm(diag[0]); me[2] += norm(diag[1]); me[3] += norm(diag[2]); me[4] += norm(diag[3]); // total diag[0] += diag[1] + diag[2] + diag[3]; me[0] += norm(diag[0]); if(calc) _me(ihel1,ihel2,2*ohel1,ohel2,ohel3) = diag[0]; } } } } } // results // initial state spin and colour average double colspin = 1./9./4.; // and C_F N_c from matrix element colspin *= 4.; // and for Z decay products if(mePartonData()[3]->coloured()) colspin *= 3.; DVector save; for(unsigned int ix=0;ix<5;++ix) { me[ix] *= colspin; if(ix>0) save.push_back(me[ix]); } meInfo(save); return me[0]*UnitRemoval::InvE2; } InvEnergy2 MEPP2ZJet::qgME(vector & fin, vector & gin, vector & fout, vector & lm, vector & lp, bool calc) const { // diagrams to include bool gamma = _gammaZ==0 || _gammaZ==1; bool Z0 = _gammaZ==0 || _gammaZ==2; // if calculation spin corrections construct the me if(calc) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1, PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1Half)); // some integers unsigned int ihel1,ihel2,ohel1,ohel2,ohel3; // compute the leptonic photon and Z currents for speed VectorWaveFunction bcurr[2][2][2]; for(ohel2=0;ohel2<2;++ohel2) { for(ohel3=0;ohel3<2;++ohel3) { // photon current if(gamma) bcurr[0][ohel2][ohel3]= _theFFPVertex->evaluate(_mz2,1,_gamma,lp[ohel3],lm[ohel2]); // Z current if(Z0) bcurr[1][ohel2][ohel3]= _theFFZVertex->evaluate(_mz2,_widthopt,_z0,lp[ohel3],lm[ohel2]); } } // compute the matrix elements double me[5]={0.,0.,0.,0.,0.}; Complex diag[4]; SpinorWaveFunction inters; SpinorBarWaveFunction interb; Energy2 _scale=scale(); for(ihel1=0;ihel1<2;++ihel1) { for(ihel2=0;ihel2<2;++ihel2) { for(ohel1=0;ohel1<2;++ohel1) { // intermediates for the diagrams interb=_theQQGVertex->evaluate(_scale,5,mePartonData()[2], fout[ohel1],gin[ihel2]); inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[0], fin[ihel1],gin[ihel2]); for(ohel2=0;ohel2<2;++ohel2) { for(ohel3=0;ohel3<2;++ohel3) { diag[0]=gamma ? _theFFPVertex->evaluate(_mz2,fin[ihel1],interb, bcurr[0][ohel2][ohel3]) : 0.; diag[1]=Z0 ? _theFFZVertex->evaluate(_mz2,fin[ihel1],interb, bcurr[1][ohel2][ohel3]) : 0.; diag[2]=gamma ? _theFFPVertex->evaluate(_mz2,inters,fout[ohel1], bcurr[0][ohel2][ohel3]) : 0.; diag[3]=Z0 ? _theFFZVertex->evaluate(_mz2,inters,fout[ohel1], bcurr[1][ohel2][ohel3]) : 0.; // diagram contributions me[1] += norm(diag[0]); me[2] += norm(diag[1]); me[3] += norm(diag[2]); me[4] += norm(diag[3]); // total diag[0] += diag[1] + diag[2] + diag[3]; me[0] += norm(diag[0]); if(calc) _me(ihel1,2*ihel2,ohel1,ohel2,ohel3) = diag[0]; } } } } } // results // initial state spin and colour average double colspin = 1./24./4.; // and C_F N_c from matrix element colspin *= 4.; // and for Z decay products if(mePartonData()[3]->coloured()) colspin *= 3.; DVector save; for(unsigned int ix=0;ix<5;++ix) { me[ix] *= colspin; if(ix>0) save.push_back(me[ix]); } meInfo(save); return me[0]*UnitRemoval::InvE2; } InvEnergy2 MEPP2ZJet::qbargME(vector & fin, vector & gin, vector & fout, vector & lm, vector & lp, bool calc) const { // diagrams to include bool gamma = _gammaZ==0 || _gammaZ==1; bool Z0 = _gammaZ==0 || _gammaZ==2; // if calculation spin corrections construct the me if(calc) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1, PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1Half)); // some integers unsigned int ihel1,ihel2,ohel1,ohel2,ohel3; // compute the leptonic photon and Z currents for speed VectorWaveFunction bcurr[2][2][2]; for(ohel2=0;ohel2<2;++ohel2) { for(ohel3=0;ohel3<2;++ohel3) { // photon current if(gamma) bcurr[0][ohel2][ohel3]= _theFFPVertex->evaluate(_mz2,1,_gamma,lp[ohel3],lm[ohel2]); // Z current if(Z0) bcurr[1][ohel2][ohel3]= _theFFZVertex->evaluate(_mz2,_widthopt,_z0,lp[ohel3],lm[ohel2]); } } // compute the matrix elements double me[5]={0.,0.,0.,0.,0.}; Complex diag[4]; SpinorWaveFunction inters; SpinorBarWaveFunction interb; Energy2 _scale=scale(); for(ihel1=0;ihel1<2;++ihel1) { for(ihel2=0;ihel2<2;++ihel2) { for(ohel1=0;ohel1<2;++ohel1) { // intermediates for the diagrams inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[2], fout[ohel1],gin[ihel2]); interb=_theQQGVertex->evaluate(_scale,5,mePartonData()[0], fin[ihel1],gin[ihel2]); for(ohel2=0;ohel2<2;++ohel2) { for(ohel3=0;ohel3<2;++ohel3) { diag[0]= gamma ? _theFFPVertex->evaluate(_mz2,inters,fin[ihel1], bcurr[0][ohel2][ohel3]) : 0.; diag[1]= Z0 ? _theFFZVertex->evaluate(_mz2,inters,fin[ihel1], bcurr[1][ohel2][ohel3]) : 0.; diag[2]= gamma ? _theFFPVertex->evaluate(_mz2,fout[ohel1],interb, bcurr[0][ohel2][ohel3]) : 0.; diag[3]= Z0 ? _theFFZVertex->evaluate(_mz2,fout[ohel1],interb, bcurr[1][ohel2][ohel3]) : 0.; // diagram contributions me[1] += norm(diag[0]); me[2] += norm(diag[1]); me[3] += norm(diag[2]); me[4] += norm(diag[3]); // total diag[0] += diag[1] + diag[2] + diag[3]; me[0] += norm(diag[0]); if(calc) _me(ihel1,2*ihel2,ohel1,ohel2,ohel3) = diag[0]; } } } } } // results // initial state spin and colour average double colspin = 1./24./4.; // and C_F N_c from matrix element colspin *= 4.; // and for Z decay products if(mePartonData()[3]->coloured()) colspin*=3.; DVector save; for(unsigned int ix=0;ix<5;++ix) { me[ix] *= colspin; if(ix>0) save.push_back(me[ix]); } meInfo(save); return me[0]*UnitRemoval::InvE2; } void MEPP2ZJet::constructVertex(tSubProPtr sub) { // extract the particles in the hard process ParticleVector hard(5); // incoming hard[0]=sub->incoming().first; hard[1]=sub->incoming().second; if((hard[0]->id()<0&&hard[1]->id()<=6)|| hard[0]->id()==ParticleID::g) swap(hard[0],hard[1]); // outgoing for(unsigned int ix=0;ix<3;++ix) { unsigned int iloc; PPtr mother=sub->outgoing()[ix]->parents()[0]; if(mother&&(mother->id()==ParticleID::gamma||mother->id()==ParticleID::Z0)) { if(sub->outgoing()[ix]->id()>0) iloc=3; else iloc=4; } else iloc=2; hard[iloc]=sub->outgoing()[ix]; } // wavefunctions for the Z decay products vector lm; vector lp; SpinorBarWaveFunction(lm,hard[3],outgoing,true,true); SpinorWaveFunction (lp,hard[4],outgoing,true,true); // identify hard process and calculate matrix element // q g to q Z if(hard[0]->id()<=6&&hard[0]->id()>0&&hard[1]->id()==ParticleID::g) { vector fin; vector gin; vector fout; SpinorWaveFunction (fin ,hard[0],incoming,false,true); VectorWaveFunction (gin ,hard[1],incoming,false,true,true); SpinorBarWaveFunction (fout,hard[2],outgoing,true ,true); gin[1]=gin[2]; qgME(fin,gin,fout,lm,lp,true); } // qbar g to qbar Z else if(hard[0]->id()>=-6&&hard[0]->id()<0&&hard[1]->id()==ParticleID::g) { vector ain; vector gin; vector aout; SpinorBarWaveFunction(ain ,hard[0],incoming,false,true); VectorWaveFunction (gin ,hard[1],incoming,false,true,true); SpinorWaveFunction (aout,hard[2],outgoing,true ,true); gin[1]=gin[2]; qbargME(ain,gin,aout,lm,lp,true); } // q qbar to g Z else { vector fin; vector ain; vector gout; SpinorWaveFunction (fin ,hard[0],incoming,false,true); SpinorBarWaveFunction(ain ,hard[1],incoming,false,true); VectorWaveFunction (gout,hard[2],outgoing,true ,true,true); gout[1]=gout[2]; qqbarME(fin,ain,gout,lm,lp,true); } // construct the vertex HardVertexPtr hardvertex=new_ptr(HardVertex()); // set the matrix element for the vertex hardvertex->ME(_me); // set the pointers and to and from the vertex for(unsigned int ix=0;ix<5;++ix) (hard[ix]->spinInfo())->productionVertex(hardvertex); } diff --git a/MatrixElement/HwMEBase.cc b/MatrixElement/HwMEBase.cc --- a/MatrixElement/HwMEBase.cc +++ b/MatrixElement/HwMEBase.cc @@ -1,301 +1,301 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the HwMEBase class. // #include "HwMEBase.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "ThePEG/PDT/EnumParticles.h" #include "Herwig/PDT/GenericMassGenerator.h" #include "ThePEG/Cuts/Cuts.h" #include "Herwig/Shower/Core/Base/Branching.h" #include "Herwig/Shower/RealEmissionProcess.h" using namespace Herwig; void HwMEBase::persistentOutput(PersistentOStream & os) const { os << massOption_ << rescaleOption_; } void HwMEBase::persistentInput(PersistentIStream & is, int) { is >> massOption_ >> rescaleOption_; } AbstractClassDescription HwMEBase::initHwMEBase; // Definition of the static class description member. void HwMEBase::Init() { static ClassDocumentation documentation ("The HwMEBase class is the base class for matrix elements in Herwig" " and provides the virtual members for hard radiation corrections in the" " shower."); } int HwMEBase::nDim() const { unsigned ndim = 1; for(unsigned int ix=0;ix & masses, double & mjac, const double *r) { assert(massOption_.size()+2==mePartonData().size()); mjac = 1.; masses.clear(); masses.resize(massOption_.size(),ZERO); Energy ecm = sqrt(sHat()); Energy emin(ZERO); int noff(0); for(unsigned int ix=0;ixhardProcessMass(); emin += masses[ix]; } else if (massOption_[ix]==2) { emin += mePartonData()[ix+2]->massMin(); ++noff; } } // check allowed if(emin>ecm) return false; // if nothing off-shell return if(noff==0) return true; int iloc = nDim()-noff; emin = ecm - emin; // generate the masses for(unsigned int ix=0;ixmassMin(); emin += mmin; Energy mmax = min(mePartonData()[ix+2]->massMax(),emin); if(mmin>mmax) return false; tGenericMassGeneratorPtr gen = mePartonData()[ix+2]->massGenerator() ? dynamic_ptr_cast(mePartonData()[ix+2]->massGenerator()) : tGenericMassGeneratorPtr(); if(gen) { double jtemp(0.); masses[ix] = gen->mass(jtemp,*mePartonData()[ix+2],mmin,mmax,r[iloc]); mjac *= jtemp; } else { Energy mon(mePartonData()[ix+2]->hardProcessMass()); Energy width(mePartonData()[ix+2]->width()); double rhomin = atan2((sqr(mmin)-sqr(mon)), mon*width); double rhomax = atan2((sqr(mmax)-sqr(mon)), mon*width); masses[ix] = sqrt(mon*width*tan(rhomin+r[iloc]*(rhomax-rhomin))+sqr(mon)); mjac *= (rhomax-rhomin)/Constants::pi; } emin -= masses[ix]; if(emin masses; double mjac(0.); if(!generateMasses(masses,mjac,r)) return false; // set up the momenta for ( int i = 2, N = meMomenta().size(); i < N; ++i ) { meMomenta()[i] = Lorentz5Momentum(masses[i-2]); } double ctmin = -1.0, ctmax = 1.0; Energy q = ZERO; try { q = SimplePhaseSpace:: getMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass()); } - catch ( ImpossibleKinematics ) { + catch ( ImpossibleKinematics & e) { return false; } Energy e = sqrt(sHat())/2.0; Energy2 m22 = meMomenta()[2].mass2(); Energy2 m32 = meMomenta()[3].mass2(); Energy2 e0e2 = 2.0*e*sqrt(sqr(q) + m22); Energy2 e1e2 = 2.0*e*sqrt(sqr(q) + m22); Energy2 e0e3 = 2.0*e*sqrt(sqr(q) + m32); Energy2 e1e3 = 2.0*e*sqrt(sqr(q) + m32); Energy2 pq = 2.0*e*q; Energy2 thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[2]); if ( thmin > ZERO ) ctmax = min(ctmax, (e0e2 - m22 - thmin)/pq); thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[2]); if ( thmin > ZERO ) ctmin = max(ctmin, (thmin + m22 - e1e2)/pq); thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[3]); if ( thmin > ZERO ) ctmax = min(ctmax, (e1e3 - m32 - thmin)/pq); thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[3]); if ( thmin > ZERO ) ctmin = max(ctmin, (thmin + m32 - e0e3)/pq); Energy ptmin = max(lastCuts().minKT(mePartonData()[2]), lastCuts().minKT(mePartonData()[3])); if ( ptmin > ZERO ) { double ctm = 1.0 - sqr(ptmin/q); if ( ctm <= 0.0 ) return false; ctmin = max(ctmin, -sqrt(ctm)); ctmax = min(ctmax, sqrt(ctm)); } double ymin2 = lastCuts().minYStar(mePartonData()[2]); double ymax2 = lastCuts().maxYStar(mePartonData()[2]); double ymin3 = lastCuts().minYStar(mePartonData()[3]); double ymax3 = lastCuts().maxYStar(mePartonData()[3]); double ytot = lastCuts().Y() + lastCuts().currentYHat(); if ( ymin2 + ytot > -0.9*Constants::MaxRapidity ) ctmin = max(ctmin, sqrt(sqr(q) + m22)*tanh(ymin2)/q); if ( ymax2 + ytot < 0.9*Constants::MaxRapidity ) ctmax = min(ctmax, sqrt(sqr(q) + m22)*tanh(ymax2)/q); if ( ymin3 + ytot > -0.9*Constants::MaxRapidity ) ctmax = min(ctmax, sqrt(sqr(q) + m32)*tanh(-ymin3)/q); if ( ymax3 + ytot < 0.9*Constants::MaxRapidity ) ctmin = max(ctmin, sqrt(sqr(q) + m32)*tanh(-ymax3)/q); if ( ctmin >= ctmax ) return false; double cth = getCosTheta(ctmin, ctmax, r[0]); Energy pt = q*sqrt(1.0-sqr(cth)); phi(rnd(2.0*Constants::pi)); meMomenta()[2].setVect(Momentum3( pt*sin(phi()), pt*cos(phi()), q*cth)); meMomenta()[3].setVect(Momentum3(-pt*sin(phi()), -pt*cos(phi()), -q*cth)); meMomenta()[2].rescaleEnergy(); meMomenta()[3].rescaleEnergy(); vector out(2); out[0] = meMomenta()[2]; out[1] = meMomenta()[3]; tcPDVector tout(2); tout[0] = mePartonData()[2]; tout[1] = mePartonData()[3]; if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) ) return false; tHat(pq*cth + m22 - e0e2); uHat(m22 + m32 - sHat() - tHat()); jacobian((pq/sHat())*Constants::pi*jacobian()*mjac); // compute the rescaled momenta return rescaleMomenta(meMomenta(),mePartonData()); } bool HwMEBase::rescaleMomenta(const vector & momenta, const cPDVector & data) { assert(momenta.size()==4&&data.size()==4); // default just use the ones we generated rescaledMomenta_=momenta; if(rescaleOption_==1) return true; Energy mnew[2] = {0*MeV, ZERO}; if(rescaleOption_==0) { mnew[0] = ZERO; mnew[1] = ZERO; } else if(rescaleOption_==2) { mnew[0] = data[2]->hardProcessMass(); mnew[1] = data[3]->hardProcessMass(); } else if(rescaleOption_==3) { if(abs(data[2]->id())!=abs(data[3]->id())) return true; mnew[0] = 0.5*(momenta[2].mass()+momenta[3].mass()); mnew[1] = mnew[0]; } else { assert(false); } Lorentz5Momentum pcm(momenta[2]+momenta[3]); Energy m0=pcm.m(); if(m01e-10*(rescaledMomenta_[2].t()+rescaledMomenta_[2].mass())) rescaledMomenta_[2].rescaleRho(); else { rescaledMomenta_[2].setX(ZERO); rescaledMomenta_[2].setY(ZERO); rescaledMomenta_[2].setZ(ZERO); } rescaledMomenta_[2].boost(-bv); rescaledMomenta_[3].boost(bv); rescaledMomenta_[3].setMass(mnew[1]); rescaledMomenta_[3].setE(0.5*(sqr(m0)-sqr(mnew[0])+sqr(mnew[1]))/m0); if(rescaledMomenta_[3].t()-rescaledMomenta_[3].mass()>1e-10*(rescaledMomenta_[3].t()+rescaledMomenta_[3].mass())) rescaledMomenta_[3].rescaleRho(); else { rescaledMomenta_[3].setX(ZERO); rescaledMomenta_[3].setY(ZERO); rescaledMomenta_[3].setZ(ZERO); } rescaledMomenta_[3].boost(-bv); return true; } double HwMEBase::getCosTheta(double ctmin, double ctmax, const double r) { double cth = 0.0; static const double eps = 1.0e-6; if ( 1.0 + ctmin <= eps && 1.0 - ctmax <= eps ) { jacobian(jacobian()*(ctmax - ctmin)); cth = ctmin + r*(ctmax - ctmin); } else if ( 1.0 + ctmin <= eps ) { cth = 1.0 - (1.0 - ctmax)*pow((1.0 - ctmin)/(1.0 - ctmax), r); jacobian(jacobian()*log((1.0 - ctmin)/(1.0 - ctmax))*(1.0 - cth)); } else if ( 1.0 - ctmax <= eps ) { cth = -1.0 + (1.0 + ctmin)*pow((1.0 + ctmax)/(1.0 + ctmin), r); jacobian(jacobian()*log((1.0 + ctmax)/(1.0 + ctmin))*(1.0 + cth)); } else { double zmin = 0.5*(1.0 - ctmax); double zmax = 0.5*(1.0 - ctmin); double A1 = -ctmin/(zmax*(1.0-zmax)); double A0 = -ctmax/(zmin*(1.0-zmin)); double A = r*(A1 - A0) + A0; double z = A < 2.0? 2.0/(sqrt(sqr(A) + 4.0) + 2 - A): 0.5*(A - 2.0 + sqrt(sqr(A) + 4.0))/A; cth = 1.0 - 2.0*z; jacobian(jacobian()*2.0*(A1 - A0)*sqr(z)*sqr(1.0 - z)/(sqr(z) + sqr(1.0 - z))); } return cth; } bool HwMEBase::softMatrixElementVeto(ShowerProgenitorPtr, ShowerParticlePtr,Branching) { assert(false); return false; } RealEmissionProcessPtr HwMEBase::generateHardest(RealEmissionProcessPtr,ShowerInteraction) { assert(false); return RealEmissionProcessPtr(); } RealEmissionProcessPtr HwMEBase::applyHardMatrixElementCorrection(RealEmissionProcessPtr) { assert(false); return RealEmissionProcessPtr(); } void HwMEBase::initializeMECorrection(RealEmissionProcessPtr , double & , double & ) { assert(false); } diff --git a/MatrixElement/MEfftoVH.cc b/MatrixElement/MEfftoVH.cc --- a/MatrixElement/MEfftoVH.cc +++ b/MatrixElement/MEfftoVH.cc @@ -1,421 +1,421 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEfftoVH class. // #include "MEfftoVH.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/MatrixElement/Tree2toNDiagram.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/Cuts/Cuts.h" #include "Herwig/Models/StandardModel/StandardModel.h" #include "Herwig/MatrixElement/HardVertex.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/PDF/PolarizedBeamParticleData.h" using namespace Herwig; void MEfftoVH::persistentOutput(PersistentOStream & os) const { os << _shapeopt << _wplus << _wminus << _z0 << _higgs << _vertexFFW << _vertexFFZ << _vertexWWH << _maxflavour << ounit(_mh,GeV) << ounit(_wh,GeV) << _hmass; } void MEfftoVH::persistentInput(PersistentIStream & is, int) { is >> _shapeopt >> _wplus >> _wminus >> _z0 >> _higgs >> _vertexFFW >> _vertexFFZ >> _vertexWWH >> _maxflavour >> iunit(_mh,GeV) >> iunit(_wh,GeV) >> _hmass; } AbstractClassDescription MEfftoVH::initMEfftoVH; // Definition of the static class description member. void MEfftoVH::Init() { static ClassDocumentation documentation ("The MEfftoVH class is the base class for the Bjirken process f fbar -> V H"); static Switch interfaceShapeOption ("ShapeScheme", "Option for the treatment of the Higgs resonance shape", &MEfftoVH::_shapeopt, 2, false, false); static SwitchOption interfaceStandardShapeFixed (interfaceShapeOption, "FixedBreitWigner", "Breit-Wigner s-channel resonanse", 1); static SwitchOption interfaceStandardShapeRunning (interfaceShapeOption, "MassGenerator", "Use the mass generator to give the shape", 2); static SwitchOption interfaceStandardShapeOnShell (interfaceShapeOption, "OnShell", "Produce an on-shell Higgs boson", 0); static Parameter interfaceMaxFlavour ( "MaxFlavour", "The heaviest incoming quark flavour this matrix element is allowed to handle " "(if applicable).", &MEfftoVH::_maxflavour, 5, 1, 5, false, false, true); } unsigned int MEfftoVH::orderInAlphaS() const { return 0; } unsigned int MEfftoVH::orderInAlphaEW() const { return 3; } Energy2 MEfftoVH::scale() const { return sHat(); } int MEfftoVH::nDim() const { return 4 + int(_shapeopt>0); } void MEfftoVH::setKinematics() { DrellYanBase::setKinematics(); } Selector MEfftoVH::diagrams(const DiagramVector & diags) const { Selector sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ) sel.insert(1.0, i); return sel; } Selector MEfftoVH::colourGeometries(tcDiagPtr ) const { static ColourLines c1(""); static ColourLines c2("6 -7"); static ColourLines c3("1 -2"); static ColourLines c4("1 -2, 6 -7"); Selector sel; if(mePartonData()[0]->coloured()) { if(mePartonData()[4]->coloured()) sel.insert(1.0, &c4); else sel.insert(1.0, &c3); } else { if(mePartonData()[4]->coloured()) sel.insert(1.0, &c2); else sel.insert(1.0, &c1); } return sel; } void MEfftoVH::doinit() { DrellYanBase::doinit(); // get the vedrtex pointers from the SM object tcHwSMPtr hwsm= dynamic_ptr_cast(standardModel()); // do the initialisation if(hwsm) { _vertexFFW = hwsm->vertexFFW(); _vertexFFZ = hwsm->vertexFFZ(); } else throw InitException() << "Wrong type of StandardModel object in " << "MEfftoVH::doinit() the Herwig" << " version must be used" << Exception::runerror; // get the particle data objects for the intermediates _wplus = getParticleData(ParticleID::Wplus ); _wminus = getParticleData(ParticleID::Wminus); _z0 = getParticleData(ParticleID::Z0); // higgs stuff _mh = _higgs->mass(); _wh = _higgs->width(); if(_higgs->massGenerator()) { _hmass=dynamic_ptr_cast(_higgs->massGenerator()); } if(_shapeopt==2&&!_hmass) throw InitException() << "If using the mass generator for the line shape in MEfftoVH::doinit()" << "the mass generator must be an instance of the GenericMassGenerator class" << Exception::runerror; } double MEfftoVH::me2() const { vector fin,aout; vector ain,fout; SpinorWaveFunction q(meMomenta()[0],mePartonData()[0],incoming); SpinorBarWaveFunction qbar(meMomenta()[1],mePartonData()[1],incoming); SpinorBarWaveFunction f(meMomenta()[3],mePartonData()[3],outgoing); SpinorWaveFunction fbar(meMomenta()[4],mePartonData()[4],outgoing); for(unsigned int ix=0;ix<2;++ix) { q.reset(ix) ; fin.push_back(q); qbar.reset(ix); ain.push_back(qbar); f.reset(ix) ;fout.push_back(f); fbar.reset(ix);aout.push_back(fbar); } return helicityME(fin,ain,fout,aout,false)*sHat()*UnitRemoval::InvE2; } double MEfftoVH::helicityME(vector & fin , vector & ain , vector & fout, vector & aout, bool calc) const { // scale Energy2 mb2(scale()); // matrix element to be stored ProductionMatrixElement menew(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin0, PDT::Spin1Half,PDT::Spin1Half); // work out the id of the vector boson int incharge = mePartonData()[0]->iCharge()+mePartonData()[1]->iCharge(); tcPDPtr vec; if(incharge==0) vec = _z0; else if(incharge>0) vec = _wplus; else vec = _wminus; // vertex for vector boson AbstractFFVVertexPtr vertex = vec==_z0 ? _vertexFFZ : _vertexFFW; // wavefunction for the scalar ScalarWaveFunction higgs(meMomenta()[2],mePartonData()[2],1.,outgoing); // calculate the matrix element VectorWaveFunction inter[2]; unsigned int ihel1,ihel2,ohel1,ohel2; Complex diag; double me(0.); for(ihel1=0;ihel1<2;++ihel1) { for(ihel2=0;ihel2<2;++ihel2) { // wavefunction for the intermediate 1st vector inter[0] = vertex->evaluate(mb2,1,vec,fin[ihel1],ain[ihel2]); // boson decay piece for(ohel1=0;ohel1<2;++ohel1) { for(ohel2=0;ohel2<2;++ohel2) { inter[1] = vertex->evaluate(sqr(vec->mass()),1,vec, aout[ohel2],fout[ohel1]); diag = _vertexWWH->evaluate(mb2,inter[1],inter[0],higgs); me += norm(diag); menew(ihel1,ihel2,0,ohel1,ohel2) = diag; } } } } // spin factor me *=0.25; tcPolarizedBeamPDPtr beam[2] = {dynamic_ptr_cast(mePartonData()[0]), dynamic_ptr_cast(mePartonData()[1])}; if( beam[0] || beam[1] ) { RhoDMatrix rho[2] = {beam[0] ? beam[0]->rhoMatrix() : RhoDMatrix(mePartonData()[0]->iSpin()), beam[1] ? beam[1]->rhoMatrix() : RhoDMatrix(mePartonData()[1]->iSpin())}; me = menew.average(rho[0],rho[1]); } // incoming colour factor if(mePartonData()[0]->coloured()) me /= 3.; // outgoing colour factor if(mePartonData()[3]->coloured()) me *= 3.; if(calc) _me.reset(menew); return me; } void MEfftoVH::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]); hard.push_back(sub->outgoing()[2]); // ensure right order if(hard[0]->id()<0) swap(hard[0],hard[1]); if(hard[3]->dataPtr()->iSpin()==PDT::Spin0) swap(hard[2],hard[3]); if(hard[4]->dataPtr()->iSpin()==PDT::Spin0) swap(hard[2],hard[4]); if(hard[3]->id()<0) swap(hard[3],hard[4]); vector fin,aout; vector ain,fout; SpinorWaveFunction( fin ,hard[0],incoming,false,true); SpinorBarWaveFunction(ain ,hard[1],incoming,false,true); ScalarWaveFunction( hard[2],outgoing,true); SpinorBarWaveFunction(fout,hard[3],outgoing,true ,true); SpinorWaveFunction( aout,hard[4],outgoing,true ,true); helicityME(fin,ain,fout,aout,true); // construct the vertex HardVertexPtr hardvertex=new_ptr(HardVertex()); // set the matrix element for the vertex hardvertex->ME(_me); // set the pointers and to and from the vertex for(unsigned int ix=0;ix<5;++ix) { tcSpinPtr spin = hard[ix]->spinInfo(); if(ix<2) { tcPolarizedBeamPDPtr beam = dynamic_ptr_cast(hard[ix]->dataPtr()); if(beam) spin->rhoMatrix() = beam->rhoMatrix(); } spin->productionVertex(hardvertex); } } bool MEfftoVH::generateKinematics(const double * r) { using Constants::pi; // workout the ID of the vector boson tcPDPtr vec = mePartonData()[0]->iCharge()+mePartonData()[1]->iCharge()!=0 ? _wplus : _z0; // order determined randomly Energy e = sqrt(sHat())/2.0; Energy mh(_mh),mv; double jac(1.); if(UseRandom::rndbool()) { double rhomax,rhomin; // generate the mass of the Higgs if(_shapeopt!=0) { Energy mhmax = min(2.*e-vec->massMin(),mePartonData()[2]->massMax()); Energy mhmin = max(ZERO ,mePartonData()[2]->massMin()); if(mhmax<=mhmin) return false; rhomin = atan2((sqr(mhmin)-sqr(_mh)), _mh*_wh); rhomax = atan2((sqr(mhmax)-sqr(_mh)), _mh*_wh); mh = sqrt(_mh*_wh*tan(rhomin+r[3]*(rhomax-rhomin))+sqr(_mh)); jac *= rhomax-rhomin; } // generate the mass of the vector boson Energy2 mvmax2 = sqr(min(2.*e-mh,vec->massMax())); Energy2 mvmin2 = sqr(vec->massMin()); if(mvmax2<=mvmin2) return false; rhomin = atan2((mvmin2-sqr(vec->mass())), vec->mass()*vec->width()); rhomax = atan2((mvmax2-sqr(vec->mass())), vec->mass()*vec->width()); mv = sqrt(vec->mass()*vec->width()*tan(rhomin+r[1]*(rhomax-rhomin)) +sqr(vec->mass())); jac *= rhomax-rhomin; } else { // generate the mass of the vector boson Energy2 mvmax2 = sqr(min(2.*e,vec->massMax())); Energy2 mvmin2 = sqr(vec->massMin()); if(mvmax2<=mvmin2) return false; double rhomin = atan2((mvmin2-sqr(vec->mass())), vec->mass()*vec->width()); double rhomax = atan2((mvmax2-sqr(vec->mass())), vec->mass()*vec->width()); mv = sqrt(vec->mass()*vec->width()*tan(rhomin+r[1]*(rhomax-rhomin)) +sqr(vec->mass())); jac *= rhomax-rhomin; // generate the mass of the Higgs if(_shapeopt!=0) { Energy mhmax = min(2.*e-mv,mePartonData()[2]->massMax()); Energy mhmin = max(ZERO ,mePartonData()[2]->massMin()); if(mhmax<=mhmin) return false; rhomin = atan2((sqr(mhmin)-sqr(_mh)), _mh*_wh); rhomax = atan2((sqr(mhmax)-sqr(_mh)), _mh*_wh); mh = sqrt(_mh*_wh*tan(rhomin+r[3]*(rhomax-rhomin))+sqr(_mh)); jac *= rhomax-rhomin; } } if(mh+mv>2.*e) return false; // assign masses meMomenta()[2].setMass(mh); for(unsigned int ix=3;ix<5;++ix) meMomenta()[ix] = Lorentz5Momentum(mePartonData()[ix]->generateMass()); Energy q = ZERO; Lorentz5Momentum pvec(mv); try { q = SimplePhaseSpace:: getMagnitude(sHat(), meMomenta()[2].mass(), mv); } - catch ( ImpossibleKinematics ) { + catch ( ImpossibleKinematics & e) { return false; } Energy ptmin = max(lastCuts().minKT(mePartonData()[2]), lastCuts().minKT(vec)); Energy2 m22 = meMomenta()[2].mass2(); Energy2 m32 = pvec .mass2(); Energy2 e0e2 = 2.0*e*sqrt(sqr(q) + m22); Energy2 e1e2 = 2.0*e*sqrt(sqr(q) + m22); Energy2 e0e3 = 2.0*e*sqrt(sqr(q) + m32); Energy2 e1e3 = 2.0*e*sqrt(sqr(q) + m32); Energy2 pq = 2.0*e*q; double ctmin = -1.0,ctmax = 1.0; Energy2 thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[2]); if ( thmin > ZERO ) ctmax = min(ctmax, (e0e2 - m22 - thmin)/pq); thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[2]); if ( thmin > ZERO ) ctmin = max(ctmin, (thmin + m22 - e1e2)/pq); thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[3]); if ( thmin > ZERO ) ctmax = min(ctmax, (e1e3 - m32 - thmin)/pq); thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[3]); if ( thmin > ZERO ) ctmin = max(ctmin, (thmin + m32 - e0e3)/pq); if ( ptmin > ZERO ) { double ctm = 1.0 - sqr(ptmin/q); if ( ctm <= 0.0 ) return false; ctmin = max(ctmin, -sqrt(ctm)); ctmax = min(ctmax, sqrt(ctm)); } if ( ctmin >= ctmax ) return false; jacobian(1.); double cth = getCosTheta(ctmin, ctmax, r[0]); Energy pt = q*sqrt(1.0-sqr(cth)); double phi = UseRandom::rnd()*2.0*pi; meMomenta()[2].setX(pt*sin(phi)); meMomenta()[2].setY(pt*cos(phi)); meMomenta()[2].setZ(q*cth); pvec .setX(-pt*sin(phi)); pvec .setY(-pt*cos(phi)); pvec .setZ(-q*cth); meMomenta()[2].rescaleEnergy(); pvec .rescaleEnergy(); // decay of the vector boson bool test=Kinematics::twoBodyDecay(pvec,meMomenta()[3].mass(), meMomenta()[4].mass(), -1.+2*r[2],r[3]*2.*pi, meMomenta()[3],meMomenta()[4]); if(!test) return false; // check cuts vector out; tcPDVector tout; for(unsigned int ix=2;ix<5;++ix) { out .push_back(meMomenta()[ix]); tout.push_back(mePartonData()[ix]); } if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) ) return false; // jacobian factors // main piece jacobian((pq/sHat())*pi*jacobian()); // mass piece jacobian(jac*jacobian()); // decay piece Energy p2 = Kinematics::pstarTwoBodyDecay(mv,meMomenta()[3].mass(), meMomenta()[4].mass()); jacobian(p2/mv/8./sqr(pi)*jacobian()); // jacobian factor for the gauge boson jacobian((sqr(sqr(mv)-sqr(vec->mass()))+sqr(vec->mass()*vec->width()))/ (vec->mass()*vec->width())*jacobian()/sHat()); return true; } CrossSection MEfftoVH::dSigHatDR() const { using Constants::pi; // jacobian factor for the higgs InvEnergy2 bwfact(ZERO); Energy moff =meMomenta()[2].mass(); if(_shapeopt==1) { tcPDPtr h0 = mePartonData()[2]->iSpin()==PDT::Spin0 ? mePartonData()[2] : mePartonData()[3]; bwfact = h0->generateWidth(moff)*moff/pi/ (sqr(sqr(moff)-sqr(_mh))+sqr(_mh*_wh)); } else if(_shapeopt==2) { bwfact = _hmass->BreitWignerWeight(moff); } double jac1 = _shapeopt!=0 ? double(bwfact*(sqr(sqr(moff)-sqr(_mh))+sqr(_mh*_wh))/(_mh*_wh)) : 1.; // answer return jac1*me2()*jacobian()/(16.0*sqr(pi)*sHat())*sqr(hbarc); } diff --git a/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.cc --- a/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.cc +++ b/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.cc @@ -1,392 +1,392 @@ // -*- C++ -*- // // FFMassiveInvertedTildeKinematics.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 FFMassiveInvertedTildeKinematics class. // #include "FFMassiveInvertedTildeKinematics.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/MatrixElement/Matchbox/Phasespace/RandomHelpers.h" using namespace Herwig; FFMassiveInvertedTildeKinematics::FFMassiveInvertedTildeKinematics() : theFullJacobian(true) {} FFMassiveInvertedTildeKinematics::~FFMassiveInvertedTildeKinematics() {} IBPtr FFMassiveInvertedTildeKinematics::clone() const { return new_ptr(*this); } IBPtr FFMassiveInvertedTildeKinematics::fullclone() const { return new_ptr(*this); } bool FFMassiveInvertedTildeKinematics::doMap(const double * r) { // todo - SW: Sort out all of the notation in the matchbox // kinematics to match the manual, before manual release. if ( ptMax() < ptCut() ) { jacobian(0.0); return false; } Lorentz5Momentum emitter = bornEmitterMomentum(); Lorentz5Momentum spectator = bornSpectatorMomentum(); double mapping = 1.0; vector values(6); pair ptz = generatePtZ(mapping, r, &values); if ( mapping == 0.0 ) { jacobian(0.0); return false; } // pt and zPrime = qi.nk / (qi+qj).nk are the generated variables Energy pt = ptz.first; Energy2 pt2 = sqr(pt); double zPrime = ptz.second; // Define the scale Energy2 Qijk = sqr(lastScale()); // Most of the required values have been calculated in ptzAllowed double y = values[0]; double z = values[1]; double A = values[2]; double xk = values[3]; double xij = values[4]; double suijk = values[5]; double suijk2 = sqr(suijk); // masses double mui2 = sqr( realEmitterData()->hardProcessMass() / lastScale() ); double mu2 = sqr( realEmissionData()->hardProcessMass() / lastScale() ); - double muj2 = sqr( realSpectatorData()->hardProcessMass() / lastScale() ); + //double muj2 = sqr( realSpectatorData()->hardProcessMass() / lastScale() ); double Mui2 = sqr( bornEmitterData()->hardProcessMass() / lastScale() ); double Muj2 = sqr( bornSpectatorData()->hardProcessMass() / lastScale() ); // Construct reference momenta nk, nij, nt Lorentz5Momentum nij = ( suijk2 / (suijk2-Mui2*Muj2) ) * (emitter - (Mui2/suijk)*spectator); Lorentz5Momentum nk = ( suijk2 / (suijk2-Mui2*Muj2) ) * (spectator - (Muj2/suijk)*emitter); // Following notation in notes, qt = sqrt(wt)*nt double phi = 2.*Constants::pi*r[2]; Lorentz5Momentum qt = getKt(nij,nk,pt,phi); // Construct qij, qk, qi and qj Lorentz5Momentum qij = xij*nij + (Mui2/(xij*suijk))*nk; Lorentz5Momentum spe = xk*nk + (Muj2/(xk*suijk))*nij; Lorentz5Momentum em = zPrime*qij + ((pt2/Qijk + mui2 - zPrime*zPrime*Mui2)/(xij*suijk*zPrime))*nk + qt; Lorentz5Momentum emm = (1.-zPrime)*qij + ((pt2/Qijk + mu2 - sqr(1.-zPrime)*Mui2)/(xij*suijk*(1.-zPrime)))*nk - qt; em.setMass(realEmitterData()->hardProcessMass()); em.rescaleEnergy(); emm.setMass(realEmissionData()->hardProcessMass()); emm.rescaleEnergy(); spe.setMass(realSpectatorData()->hardProcessMass()); spe.rescaleEnergy(); // book realEmitterMomentum() = em; realEmissionMomentum() = emm; realSpectatorMomentum() = spe; // Compute and store the jacobian // The jacobian here corresponds to dpt2 / sqr(lastscale) NOT dpt2 / pt2. // This jacobian is the one-particle phase space // jac s.t.: dy dz = jac* dpt2/sqr(lastScale()) dz double jac = 0.0; // SW - Change in notation here that needs to be fixed (j<->nothing<->k) Energy2 mi2 = sqr(realEmitterData()->hardProcessMass()); Energy2 mj2 = sqr(realEmissionData()->hardProcessMass()); Energy2 mk2 = sqr(realSpectatorData()->hardProcessMass()); Energy2 mij2 = sqr(bornEmitterData()->hardProcessMass()); Energy2 sbar = Qijk - mi2 - mj2 -mk2; if ( theFullJacobian && bornSpectatorData()->hardProcessMass() != ZERO ) { Energy2 sijk = Qijk*suijk; double lambdaK = 1. + (mk2/sijk); double lambdaIJ = 1. + (mij2/sijk); // Compute dy/dzPrime and pt2* dy/dpt2 double dyBydzPrime = (1./sbar) * ( -pt2*(1.-2.*zPrime)/sqr(zPrime*(1.-zPrime)) - mi2/sqr(zPrime) + mj2/sqr(1.-zPrime) ); InvEnergy2 dyBydpt2 = 1./(sbar*zPrime*(1.-zPrime)); // Compute dA/dzPrime and dA/dpt2 double dABydzPrime = (sbar/sijk) * dyBydzPrime; InvEnergy2 dABydpt2 = (sbar/sijk) * dyBydpt2; // Compute dxk/dzPrime, dxk/dpt2, dxij/dzPrime and dxij/dpt2 double factor = (0.5/lambdaK) * (-1. - (1./sqrt( sqr(lambdaK + (mk2/sijk)*lambdaIJ - A) - 4.*lambdaK*lambdaIJ*mk2/sijk)) * (lambdaK + (mk2/sijk)*lambdaIJ - A) ); double dxkBydzPrime = factor * dABydzPrime; InvEnergy2 dxkBydpt2 = factor * dABydpt2; double dxijBydzPrime = (mk2/sijk) * (1./sqr(xk)) * dxkBydzPrime; InvEnergy2 dxijBydpt2 = (mk2/sijk) * (1./sqr(xk)) * dxkBydpt2; Energy2 dqiDotqkBydzPrime = xij*xk*0.5*sijk + zPrime*dxijBydzPrime*xk*0.5*sijk + zPrime*xij*dxkBydzPrime*0.5*sijk + 0.5*(mk2/sijk)*(pt2 + mi2) * (-1./(xk*xij*sqr(zPrime)) - dxkBydzPrime/(zPrime*xij*sqr(xk)) - dxijBydzPrime/(zPrime*xk*sqr(xij))); double dqiDotqkBydpt2 = dxijBydpt2*zPrime*xk*0.5*sijk + zPrime*xij*dxkBydpt2*0.5*sijk + (0.5*mk2/sijk) * (1./(zPrime*xk*xij)) * (1. + (pt2+mi2)*(-dxkBydpt2/xk - dxijBydpt2/xij) ); // Compute dzBydzPrime and dzBydpt2 Energy2 qiDotqk = (zPrime*xij*xk*sijk*0.5) + (mk2/ ( 2.*xk*xij*sijk*zPrime))*(pt2 + mi2); double dzBydzPrime = (1./sbar) * ( 2.*qiDotqk*dyBydzPrime/sqr(1.-y) + (1./(1.-y))*2.*dqiDotqkBydzPrime ); InvEnergy2 dzBydpt2 = (1./sbar) * ( 2.*qiDotqk*dyBydpt2/sqr(1.-y) + (1./(1.-y))*2.*dqiDotqkBydpt2 ); // Compute the jacobian jac = Qijk * abs(dzBydpt2*dyBydzPrime - dzBydzPrime*dyBydpt2); } // This is correct in the massless spectator case. else { jac = Qijk / (sbar*zPrime*(1.-zPrime)); } // Mapping includes only the variable changes/jacobians mapping *= jac; jacobian( (sqr(sbar)/rootOfKallen(Qijk,mij2,mk2)) * mapping*(1.-y)/(16.*sqr(Constants::pi)) / sHat() ); // Store the parameters subtractionParameters().resize(3); subtractionParameters()[0] = y; subtractionParameters()[1] = z; subtractionParameters()[2] = zPrime; return true; } Energy FFMassiveInvertedTildeKinematics::lastPt() const { Energy scale = (bornEmitterMomentum()+bornSpectatorMomentum()).m(); // masses double mui2 = sqr( realEmitterData()->hardProcessMass() / scale ); double mu2 = sqr( realEmissionData()->hardProcessMass() / scale ); double muj2 = sqr( realSpectatorData()->hardProcessMass() / scale ); double y = subtractionParameters()[0]; double zPrime = subtractionParameters()[2]; Energy ret = scale * sqrt( y * (1.-mui2-mu2-muj2) * zPrime*(1.-zPrime) - sqr(1.-zPrime)*mui2 - sqr(zPrime)*mu2 ); return ret; } double FFMassiveInvertedTildeKinematics::lastZ() const { return subtractionParameters()[2]; } Energy FFMassiveInvertedTildeKinematics::ptMax() const { Energy scale = (bornEmitterMomentum()+bornSpectatorMomentum()).m(); // masses double mui2 = sqr( realEmitterData()->hardProcessMass() / scale ); double mu2 = sqr( realEmissionData()->hardProcessMass() / scale ); double muj2 = sqr( realSpectatorData()->hardProcessMass() / scale ); Energy ptmax = rootOfKallen( mui2, mu2, sqr(1.-sqrt(muj2)) ) / ( 2.-2.*sqrt(muj2) ) * scale; return ptmax > 0.*GeV ? ptmax : 0.*GeV; } // NOTE: bounds calculated at this step may be too loose // These apply to zPrime, which is stored in lastZ() pair FFMassiveInvertedTildeKinematics::zBounds(Energy pt, Energy hardPt) const { hardPt = hardPt == ZERO ? ptMax() : min(hardPt,ptMax()); if(pt>hardPt) return make_pair(0.5,0.5); Energy scale = (bornEmitterMomentum()+bornSpectatorMomentum()).m(); // masses double mui2 = sqr( realEmitterData()->hardProcessMass() / scale ); double mu2 = sqr( realEmissionData()->hardProcessMass() / scale ); double muj2 = sqr( realSpectatorData()->hardProcessMass() / scale ); double zp = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) + rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) * sqrt( 1.-sqr(pt/hardPt) ) ) / ( 2.*sqr(1.-sqrt(muj2)) ); double zm = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) - rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) * sqrt( 1.-sqr(pt/hardPt) ) ) / ( 2.*sqr(1.-sqrt(muj2)) ); return make_pair(zm,zp); } bool FFMassiveInvertedTildeKinematics::ptzAllowed(pair ptz, vector* values) const { Energy pt = ptz.first; Energy2 pt2 = sqr(pt); double zPrime = ptz.second; // masses double mui2 = sqr( realEmitterData()->hardProcessMass() / lastScale() ); double mu2 = sqr( realEmissionData()->hardProcessMass() / lastScale() ); double muj2 = sqr( realSpectatorData()->hardProcessMass() / lastScale() ); double Mui2 = sqr( bornEmitterData()->hardProcessMass() / lastScale() ); double Muj2 = sqr( bornSpectatorData()->hardProcessMass() / lastScale() ); // Calculate the scales that we need Energy2 Qijk = sqr(lastScale()); double suijk = 0.5*( 1. - Mui2 - Muj2 + sqrt( sqr(1.-Mui2-Muj2) - 4.*Mui2*Muj2 ) ); double bar = 1.-mui2-mu2-muj2; double y = ( pt2/Qijk + sqr(1.-zPrime)*mui2 + zPrime*zPrime*mu2 ) / (zPrime*(1.-zPrime)*bar); // Calculate A:=xij*w double A = (1./(suijk*zPrime*(1.-zPrime))) * ( pt2/Qijk + zPrime*mu2 + (1.-zPrime)*mui2 - zPrime*(1.-zPrime)*Mui2 ); // Calculate the scaling factors, xk and xij double lambdaK = 1. + (Muj2/suijk); double lambdaIJ = 1. + (Mui2/suijk); double xk = (1./(2.*lambdaK)) * ( (lambdaK + (Muj2/suijk)*lambdaIJ - A) + sqrt( sqr(lambdaK + (Muj2/suijk)*lambdaIJ - A) - 4.*lambdaK*lambdaIJ*Muj2/suijk) ); double xij = 1. - ( (Muj2/suijk) * (1.-xk) / xk ); // Calculate z double z = ( (zPrime*xij*xk*suijk/2.) + (Muj2/ ( 2.*xk*xij*suijk*zPrime))*(pt2/Qijk + mui2) ) / ( (xij*xk*suijk/2.) + (Muj2/(2.*xk*xij))*(Mui2/suijk + A) ); // check (y,z) phase space boundary // TODO: is y boundary necessary? double ym = 2.*sqrt(mui2)*sqrt(mu2)/bar; double yp = 1. - 2.*sqrt(muj2)*(1.-sqrt(muj2))/bar; // These limits apply to z, not zPrime double zm = ( (2.*mui2+bar*y)*(1.-y) - sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) / ( 2.*(1.-y)*(mui2+mu2+bar*y) ); double zp = ( (2.*mui2+bar*y)*(1.-y) + sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) / ( 2.*(1.-y)*(mui2+mu2+bar*y) ); if ( yyp || zzp ) return false; assert( (*values).size() == 6); (*values)[0] = y; (*values)[1] = z; (*values)[2] = A; (*values)[3] = xk; (*values)[4] = xij; (*values)[5] = suijk; return true; } // This is used to generate pt and zPrime pair FFMassiveInvertedTildeKinematics::generatePtZ(double& jac, const double * r, vector * values) const { double kappaMin = ptCut() != ZERO ? sqr(ptCut()/ptMax()) : sqr(0.1*GeV/GeV); double kappa; using namespace RandomHelpers; if ( ptCut() > ZERO ) { pair kw = generate(inverse(0.,kappaMin,1.),r[0]); kappa = kw.first; jac *= kw.second; } else { pair kw = generate((piecewise(), flat(1e-4,kappaMin), match(inverse(0.,kappaMin,1.))),r[0]); kappa = kw.first; jac *= kw.second; } Energy pt = sqrt(kappa)*ptMax(); pair zLims = zBounds(pt); pair zw = generate(inverse(0.,zLims.first,zLims.second)+ inverse(1.,zLims.first,zLims.second),r[1]); double z = zw.first; jac *= zw.second; jac *= sqr(ptMax()/lastScale()); if( !ptzAllowed(make_pair(pt,z), values )) jac = 0.; return make_pair(pt,z); } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void FFMassiveInvertedTildeKinematics::persistentOutput(PersistentOStream &) const { } void FFMassiveInvertedTildeKinematics::persistentInput(PersistentIStream &, int) { } void FFMassiveInvertedTildeKinematics::Init() { static ClassDocumentation documentation ("FFMassiveInvertedTildeKinematics inverts the final-final tilde " "kinematics involving a massive particle."); } // *** Attention *** The following static variable is needed for the type // description system in ThePEG. Please check that the template arguments // are correct (the class and its base class), and that the constructor // arguments are correct (the class name and the name of the dynamically // loadable library where the class implementation can be found). DescribeClass describeHerwigFFMassiveInvertedTildeKinematics("Herwig::FFMassiveInvertedTildeKinematics", "Herwig.so"); diff --git a/MatrixElement/Powheg/MEPP2GammaGammaPowheg.cc b/MatrixElement/Powheg/MEPP2GammaGammaPowheg.cc --- a/MatrixElement/Powheg/MEPP2GammaGammaPowheg.cc +++ b/MatrixElement/Powheg/MEPP2GammaGammaPowheg.cc @@ -1,2047 +1,2047 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEPP2GammaGammaPowheg class. // #include "MEPP2GammaGammaPowheg.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/MatrixElement/Tree2toNDiagram.h" #include "ThePEG/Cuts/Cuts.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "Herwig/Models/StandardModel/StandardModel.h" #include "Herwig/Utilities/Maths.h" #include "Herwig/Shower/RealEmissionProcess.h" #include "ThePEG/Utilities/DescribeClass.h" using namespace Herwig; // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeMEPP2GammaGammaPowheg("Herwig::MEPP2GammaGammaPowheg", "HwMEHadron.so HwPowhegMEHadron.so"); unsigned int MEPP2GammaGammaPowheg::orderInAlphaS() const { return 0; } unsigned int MEPP2GammaGammaPowheg::orderInAlphaEW() const { return 2; } IBPtr MEPP2GammaGammaPowheg::clone() const { return new_ptr(*this); } IBPtr MEPP2GammaGammaPowheg::fullclone() const { return new_ptr(*this); } MEPP2GammaGammaPowheg::MEPP2GammaGammaPowheg() : contrib_(1), power_(0.1), process_(0), threeBodyProcess_(0), maxflavour_(5), alphaS_(0.), fixedAlphaS_(false), supressionFunction_(0), supressionScale_(0), lambda_(20.*GeV), preQCDqqbarq_(5.), preQCDqqbarqbar_(0.5), preQCDqg_(50.), preQCDgqbar_(50.), preQEDqqbarq_(40.), preQEDqqbarqbar_(0.5), preQEDqgq_(1.), preQEDgqbarqbar_(1.), minpT_(2.*GeV), scaleChoice_(0), scalePreFactor_(1.) {} void MEPP2GammaGammaPowheg::getDiagrams() const { tcPDPtr gamma = getParticleData(ParticleID::gamma); tcPDPtr g = getParticleData(ParticleID::g); for(int ix=1;ix<=maxflavour_;++ix) { tcPDPtr qk = getParticleData(ix); tcPDPtr qb = qk->CC(); // gamma gamma if(process_==0 || process_ == 1) { add(new_ptr((Tree2toNDiagram(3), qk, qk, qb, 1, gamma, 2, gamma, -1))); add(new_ptr((Tree2toNDiagram(3), qk, qk, qb, 2, gamma, 1, gamma, -2))); } // gamma +jet if(process_==0 || process_ == 2) { add(new_ptr((Tree2toNDiagram(3), qk, qb, qb, 1, gamma, 2, g, -4))); add(new_ptr((Tree2toNDiagram(3), qk, qk, qb, 2, gamma, 1, g, -5))); add(new_ptr((Tree2toNDiagram(3), qk, qk, g, 1, gamma, 2, qk, -6))); add(new_ptr((Tree2toNDiagram(2), qk, g, 1, qk, 3, gamma, 3, qk, -7))); add(new_ptr((Tree2toNDiagram(3), g, qb, qb, 2, gamma, 1, qb, -8))); add(new_ptr((Tree2toNDiagram(2), g, qb, 1, qb, 3, gamma, 3, qb, -9))); } // gamma + jet + gamma if((process_==0 && contrib_==1) || process_ == 3) { // gamma + g + gamma if(threeBodyProcess_==0 || threeBodyProcess_==1) { add(new_ptr((Tree2toNDiagram(4), qk, qk, qk, qb, 1, gamma, 2, gamma, 3, g, -10))); add(new_ptr((Tree2toNDiagram(4), qk, qk, qk, qb, 3, gamma, 2, gamma, 1, g, -12))); } // Z + q + gamma if(threeBodyProcess_==0 || threeBodyProcess_==2) { add(new_ptr((Tree2toNDiagram(4),qk,qk,qk,g,1,gamma,2,gamma,3,qk, -20))); add(new_ptr((Tree2toNDiagram(4),qk,qk,qk,g,2,gamma,1,gamma,3,qk, -21))); add(new_ptr((Tree2toNDiagram(3),qk,qk,g,1,gamma,2,qk,5,gamma,5,qk,-22))); } // Z + qbar + gamma if(threeBodyProcess_==0 || threeBodyProcess_==3) { add(new_ptr((Tree2toNDiagram(4),g,qb,qb,qb,3,gamma,2,gamma,1,qb ,-30))); add(new_ptr((Tree2toNDiagram(4),g,qb,qb,qb,2,gamma,3,gamma,1,qb ,-31))); add(new_ptr((Tree2toNDiagram(3),g,qb,qb ,2,gamma,1,qb,5,gamma,5,qb,-32))); } } } } Energy2 MEPP2GammaGammaPowheg::scale() const { Energy2 scale; if(scaleChoice_==0) { Energy pt; if(meMomenta()[2].perp(meMomenta()[0].vect())>= meMomenta()[3].perp(meMomenta()[0].vect())){ pt = meMomenta()[2].perp(meMomenta()[0].vect()); } else { pt = meMomenta()[3].perp(meMomenta()[0].vect()); } scale = sqr(pt); } else if(scaleChoice_==1) { scale = sHat(); } return scalePreFactor_*scale; } int MEPP2GammaGammaPowheg::nDim() const { return HwMEBase::nDim() + ( contrib_>=1 ? 3 : 0 ); } bool MEPP2GammaGammaPowheg::generateKinematics(const double * r) { // radiative variables if(contrib_>=1) { zTilde_ = r[nDim()-1]; vTilde_ = r[nDim()-2]; phi_ = Constants::twopi*r[nDim()-3]; } // set the jacobian jacobian(1.0); // set up the momenta for ( int i = 2, N = meMomenta().size(); i < N; ++i ) meMomenta()[i] = Lorentz5Momentum(ZERO); // generate sHat Energy2 shat(sHat()); if(mePartonData().size()==5) { double eps = sqr(meMomenta()[2].mass())/shat; jacobian(jacobian()*(1.-eps)); shat *= eps+zTilde_*(1.-eps); } // momenta of the core process double ctmin = -1.0, ctmax = 1.0; Energy q = ZERO; try { q = SimplePhaseSpace:: getMagnitude(shat, meMomenta()[2].mass(), ZERO); } - catch ( ImpossibleKinematics ) { + catch ( ImpossibleKinematics & e ) { return false; } Energy e = 0.5*sqrt(shat); Energy2 m22 = meMomenta()[2].mass2(); Energy2 e0e2 = 2.0*e*sqrt(sqr(q) + m22); Energy2 e1e2 = 2.0*e*sqrt(sqr(q) + m22); Energy2 e0e3 = 2.0*e*sqrt(sqr(q)); Energy2 e1e3 = 2.0*e*sqrt(sqr(q)); Energy2 pq = 2.0*e*q; if(mePartonData().size()==4) { Energy2 thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[2]); if ( thmin > ZERO ) ctmax = min(ctmax, (e0e2 - m22 - thmin)/pq); thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[2]); if ( thmin > ZERO ) ctmin = max(ctmin, (thmin + m22 - e1e2)/pq); thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[3]); if ( thmin > ZERO ) ctmax = min(ctmax, (e1e3 - thmin)/pq); thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[3]); if ( thmin > ZERO ) ctmin = max(ctmin, (thmin - e0e3)/pq); Energy ptmin = max(lastCuts().minKT(mePartonData()[2]), lastCuts().minKT(mePartonData()[3])); if ( ptmin > ZERO ) { double ctm = 1.0 - sqr(ptmin/q); if ( ctm <= 0.0 ) return false; ctmin = max(ctmin, -sqrt(ctm)); ctmax = min(ctmax, sqrt(ctm)); } double ymin2 = lastCuts().minYStar(mePartonData()[2]); double ymax2 = lastCuts().maxYStar(mePartonData()[2]); double ymin3 = lastCuts().minYStar(mePartonData()[3]); double ymax3 = lastCuts().maxYStar(mePartonData()[3]); double ytot = lastCuts().Y() + lastCuts().currentYHat(); if ( ymin2 + ytot > -0.9*Constants::MaxRapidity ) ctmin = max(ctmin, sqrt(sqr(q) + m22)*tanh(ymin2)/q); if ( ymax2 + ytot < 0.9*Constants::MaxRapidity ) ctmax = min(ctmax, sqrt(sqr(q) + m22)*tanh(ymax2)/q); if ( ymin3 + ytot > -0.9*Constants::MaxRapidity ) ctmax = min(ctmax, tanh(-ymin3)); if ( ymax3 + ytot < 0.9*Constants::MaxRapidity ) ctmin = max(ctmin, tanh(-ymax3)); if ( ctmin >= ctmax ) return false; } double cth = getCosTheta(ctmin, ctmax, r[0]); Energy pt = q*sqrt(1.0-sqr(cth)); phi(rnd(2.0*Constants::pi)); meMomenta()[2].setVect(Momentum3( pt*sin(phi()), pt*cos(phi()), q*cth)); meMomenta()[3].setVect(Momentum3(-pt*sin(phi()), -pt*cos(phi()), -q*cth)); meMomenta()[2].rescaleEnergy(); meMomenta()[3].rescaleEnergy(); // jacobian tHat(pq*cth + m22 - e0e2); uHat(m22 - shat - tHat()); jacobian(pq/shat*Constants::pi*jacobian()); // end for 2->2 processes if(mePartonData().size()==4) { vector out(2); out[0] = meMomenta()[2]; out[1] = meMomenta()[3]; tcPDVector tout(2); tout[0] = mePartonData()[2]; tout[1] = mePartonData()[3]; if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) ) return false; return true; } // special for 2-3 processes pair x = make_pair(lastX1(),lastX2()); // partons pair partons = make_pair(mePartonData()[0],mePartonData()[1]); // If necessary swap the particle data objects so that // first beam gives the incoming quark if(lastPartons().first ->dataPtr()!=partons.first) { swap(x.first,x.second); } // use vTilde to select the dipole for emission // gamma gamma g processes if(mePartonData()[4]->id()==ParticleID::g) { if(vTilde_<=0.5) { dipole_ = IIQCD1; vTilde_ = 4.*vTilde_; } else { dipole_ = IIQCD2; vTilde_ = 4.*(vTilde_-0.25); } jacobian(2.*jacobian()); } // gamma gamma q processes else if(mePartonData()[4]->id()>0&&mePartonData()[4]->id()<6) { if(vTilde_<=1./3.) { dipole_ = IIQCD2; vTilde_ = 3.*vTilde_; } else if(vTilde_<=2./3.) { dipole_ = IFQED1; vTilde_ = 3.*vTilde_-1.; } else { dipole_ = FIQED1; vTilde_ = 3.*vTilde_-2.; } jacobian(3.*jacobian()); } // gamma gamma qbar processes else if(mePartonData()[4]->id()<0&&mePartonData()[4]->id()>-6) { if(vTilde_<=1./3.) { dipole_ = IIQCD1; vTilde_ = 3.*vTilde_; } else if(vTilde_<=2./3.) { dipole_ = IFQED2; vTilde_ = 3.*vTilde_-1.; } else { dipole_ = FIQED2; vTilde_ = 3.*vTilde_-2.; } jacobian(3.*jacobian()); } else { assert(false); } // initial-initial dipoles if(dipole_<=4) { double z = shat/sHat(); double vt = vTilde_*(1.-z); double vJac = 1.-z; Energy pT = sqrt(shat*vt*(1.-vt-z)/z); if(pT pnew(5); pnew [0] = Lorentz5Momentum(ZERO,ZERO,0.5*rs*x.first, 0.5*rs*x.first,ZERO); pnew [1] = Lorentz5Momentum(ZERO,ZERO,-0.5*rs*x.second, 0.5*rs*x.second,ZERO) ; pnew [2] = meMomenta()[2]; pnew [3] = meMomenta()[3]; pnew [4] = Lorentz5Momentum(pT*cos(phi_),pT*sin(phi_), pT*sinh(rapidity), pT*cosh(rapidity), ZERO); pnew[4].rescaleEnergy(); Lorentz5Momentum K = pnew [0]+pnew [1]-pnew [4]; Lorentz5Momentum Kt = pcmf; Lorentz5Momentum Ksum = K+Kt; Energy2 K2 = K.m2(); Energy2 Ksum2 = Ksum.m2(); for(unsigned int ix=2;ix<4;++ix) { pnew [ix].boost(blab); pnew [ix] = pnew [ix] - 2.*Ksum*(Ksum*pnew [ix])/Ksum2 +2*K*(Kt*pnew [ix])/K2; pnew[ix].rescaleEnergy(); } pcmf = Lorentz5Momentum(ZERO,ZERO, 0.5*rs*(x.first-x.second), 0.5*rs*(x.first+x.second)); pcmf.rescaleMass(); blab = pcmf.boostVector(); for(unsigned int ix=0;ix1e-20) { rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); rot.rotateX(Constants::pi); } if(abs(1.-q.e()/q.vect().mag())>1e-6) rot.boostZ(q.e()/q.vect().mag()); pin *= rot; if(pin.perp2()/GeV2>1e-20) { Boost trans = -1./pin.e()*pin.vect(); trans.setZ(0.); rot.boost(trans); } rot.invert(); Energy Q = sqrt(-q.m2()); meMomenta()[4] = rot*Lorentz5Momentum( 0.5*Q*xT*cos(phi_), 0.5*Q*xT*sin(phi_), -0.5*Q*x2,0.5*Q*sqrt(sqr(x2)+sqr(xT))); meMomenta()[3] = rot*Lorentz5Momentum(-0.5*Q*xT*cos(phi_),-0.5*Q*xT*sin(phi_), -0.5*Q*x3,0.5*Q*sqrt(sqr(x3)+sqr(xT))); double ratio; if(dipole_<=6) { ratio = 2.*((meMomenta()[3]+meMomenta()[4])*meMomenta()[0])/sHat(); } else { ratio = 2.*((meMomenta()[3]+meMomenta()[4])*meMomenta()[1])/sHat(); } jacobian(jacobian()*ratio); } else { assert(false); } vector out(3); tcPDVector tout(3); for(unsigned int ix=0;ix<3;++ix) { out[ix] = meMomenta() [2+ix]; tout[ix] = mePartonData()[2+ix]; } return lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]); } double MEPP2GammaGammaPowheg::me2() const { // Born configurations if(mePartonData().size()==4) { // gamma gamma core process if(mePartonData()[3]->id()==ParticleID::gamma) { return 2.*Constants::twopi*alphaEM_* loGammaGammaME(mePartonData(),meMomenta(),true); } // V jet core process else if(mePartonData()[3]->id()==ParticleID::g) { return 2.*Constants::twopi*alphaS_* loGammagME(mePartonData(),meMomenta(),true); } else if(mePartonData()[3]->id()>0) { return 2.*Constants::twopi*alphaS_* loGammaqME(mePartonData(),meMomenta(),true); } else if(mePartonData()[3]->id()<0) { return 2.*Constants::twopi*alphaS_* loGammaqbarME(mePartonData(),meMomenta(),true); } else assert(false); } // hard emission configurations else { if(mePartonData()[4]->id()==ParticleID::g) return sHat()*realGammaGammagME (mePartonData(),meMomenta(),dipole_,Hard,true); else if(mePartonData()[4]->id()>0&&mePartonData()[4]->id()<6) return sHat()*realGammaGammaqME (mePartonData(),meMomenta(),dipole_,Hard,true); else if(mePartonData()[4]->id()<0&&mePartonData()[4]->id()>-6) return sHat()*realGammaGammaqbarME(mePartonData(),meMomenta(),dipole_,Hard,true); else assert(false); } } CrossSection MEPP2GammaGammaPowheg::dSigHatDR() const { // couplings if(!fixedAlphaS_) alphaS_ = SM().alphaS(scale()); alphaEM_ = SM().alphaEM(); // cross section CrossSection preFactor = jacobian()/(16.0*sqr(Constants::pi)*sHat())*sqr(hbarc); loME_ = me2(); if( contrib_== 0 || mePartonData().size()==5 || (mePartonData().size()==4&& mePartonData()[3]->coloured())) return loME_*preFactor; else return NLOWeight()*preFactor; } Selector MEPP2GammaGammaPowheg::diagrams(const DiagramVector & diags) const { if(mePartonData().size()==4) { if(mePartonData()[3]->id()==ParticleID::gamma) { Selector sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ){ sel.insert(meInfo()[abs(diags[i]->id())], i); } return sel; } else { Selector sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ){ sel.insert(meInfo()[abs(diags[i]->id())%2], i); } return sel; } } else { Selector sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ) { if(abs(diags[i]->id()) == 10 && dipole_ == IIQCD2 ) sel.insert(1., i); else if(abs(diags[i]->id()) == 12 && dipole_ == IIQCD1 ) sel.insert(1., i); else if(abs(diags[i]->id()) == 20 && dipole_ == IIQCD2 ) sel.insert(1., i); else if(abs(diags[i]->id()) == 21 && dipole_ == IFQED1 ) sel.insert(1., i); else if(abs(diags[i]->id()) == 22 && dipole_ == FIQED1 ) sel.insert(1., i); else sel.insert(0., i); } return sel; } } Selector MEPP2GammaGammaPowheg::colourGeometries(tcDiagPtr diag) const { // colour lines for V gamma static ColourLines cs("1 -2"); static ColourLines ct("1 2 -3"); // colour lines for q qbar -> V g static const ColourLines cqqbar[2]={ColourLines("1 -2 5,-3 -5"), ColourLines("1 5,-5 2 -3")}; // colour lines for q g -> V q static const ColourLines cqg [2]={ColourLines("1 2 -3,3 5"), ColourLines("1 -2,2 3 5")}; // colour lines for g qbar -> V qbar static const ColourLines cgqbar[2]={ColourLines("-3 -2 1,-1 -5"), ColourLines("-2 1,-1 -3 -5")}; // colour lines for q qbar -> V gamma g static const ColourLines cqqbarg[4]={ColourLines("1 2 3 7,-4 -7"), ColourLines("1 2 7,-4 3 -7"), ColourLines("1 7,-4 3 2 -7"), ColourLines("1 2 7,-4 3 -7")}; // colour lines for q g -> V gamma q static const ColourLines cqgq [3]={ColourLines("1 2 3 -4,4 7"), ColourLines("1 2 3 -4,4 7"), ColourLines("1 2 -3,3 5 7")}; // colour lines for gbar -> V gamma qbar static const ColourLines cqbargqbar[3]={ColourLines("1 -2 -3 -4,-1 -7"), ColourLines("1 -2 -3 -4,-1 -7"), ColourLines("1 -2 -3,-1 -5 -7")}; Selector sel; switch(abs(diag->id())) { case 1 :case 2 : sel.insert(1.0, &ct); break; case 3 : sel.insert(1.0, &cs); break; case 4 : sel.insert(1.0, &cqqbar[0]); break; case 5: sel.insert(1.0, &cqqbar[1]); break; case 6: sel.insert(1.0, &cqg[0]); break; case 7: sel.insert(1.0, &cqg[1]); break; case 8: sel.insert(1.0, &cgqbar[0]); break; case 9: sel.insert(1.0, &cgqbar[1]); break; case 10: case 11: case 12: case 13: sel.insert(1.0, &cqqbarg[abs(diag->id())-10]); break; case 20: case 21: case 22: sel.insert(1.0, &cqgq[abs(diag->id())-20]); break; case 30: case 31: case 32: sel.insert(1.0, &cqbargqbar[abs(diag->id())-30]); break; default: assert(false); } return sel; } void MEPP2GammaGammaPowheg::persistentOutput(PersistentOStream & os) const { os << FFPvertex_ << FFGvertex_ << contrib_ << power_ << gluon_ << prefactor_ << process_ << threeBodyProcess_<< maxflavour_ << alphaS_ << fixedAlphaS_ << supressionFunction_ << supressionScale_ << ounit(lambda_,GeV) << alphaQCD_ << alphaQED_ << ounit(minpT_,GeV) << preQCDqqbarq_ << preQCDqqbarqbar_ << preQCDqg_ << preQCDgqbar_ << preQEDqqbarq_ << preQEDqqbarqbar_ << preQEDqgq_ << preQEDgqbarqbar_ << scaleChoice_ << scalePreFactor_; } void MEPP2GammaGammaPowheg::persistentInput(PersistentIStream & is, int) { is >> FFPvertex_ >> FFGvertex_ >> contrib_ >> power_ >> gluon_ >> prefactor_ >> process_ >> threeBodyProcess_ >> maxflavour_ >> alphaS_ >> fixedAlphaS_ >> supressionFunction_ >> supressionScale_ >> iunit(lambda_,GeV) >> alphaQCD_ >> alphaQED_ >> iunit(minpT_,GeV) >> preQCDqqbarq_ >> preQCDqqbarqbar_ >> preQCDqg_ >> preQCDgqbar_ >> preQEDqqbarq_ >> preQEDqqbarqbar_ >> preQEDqgq_ >> preQEDgqbarqbar_ >> scaleChoice_ >> scalePreFactor_; } void MEPP2GammaGammaPowheg::Init() { static ClassDocumentation documentation ("TheMEPP2GammaGammaPowheg class implements gamma gamma production at NLO"); static Switch interfaceProcess ("Process", "Which processes to include", &MEPP2GammaGammaPowheg::process_, 0, false, false); static SwitchOption interfaceProcessAll (interfaceProcess, "All", "Include all the processes", 0); static SwitchOption interfaceProcessGammaGamma (interfaceProcess, "GammaGamma", "Only include gamma gamma", 1); static SwitchOption interfaceProcessVJet (interfaceProcess, "VJet", "Only include gamma + jet", 2); static SwitchOption interfaceProcessHard (interfaceProcess, "Hard", "Only include hard radiation contributions", 3); static Switch interfaceThreeBodyProcess ("ThreeBodyProcess", "The possible three body processes to include", &MEPP2GammaGammaPowheg::threeBodyProcess_, 0, false, false); static SwitchOption interfaceThreeBodyProcessAll (interfaceThreeBodyProcess, "All", "Include all processes", 0); static SwitchOption interfaceThreeBodyProcessqqbar (interfaceThreeBodyProcess, "qqbar", "Only include q qbar -> gamma gamma g processes", 1); static SwitchOption interfaceThreeBodyProcessqg (interfaceThreeBodyProcess, "qg", "Only include q g -> gamma gamma q processes", 2); static SwitchOption interfaceThreeBodyProcessgqbar (interfaceThreeBodyProcess, "gqbar", "Only include g qbar -> gamma gamma qbar processes", 3); static Switch interfaceContribution ("Contribution", "Which contributions to the cross section to include", &MEPP2GammaGammaPowheg::contrib_, 1, false, false); static SwitchOption interfaceContributionLeadingOrder (interfaceContribution, "LeadingOrder", "Just generate the leading order cross section", 0); static SwitchOption interfaceContributionPositiveNLO (interfaceContribution, "PositiveNLO", "Generate the positive contribution to the full NLO cross section", 1); static SwitchOption interfaceContributionNegativeNLO (interfaceContribution, "NegativeNLO", "Generate the negative contribution to the full NLO cross section", 2); static Parameter interfaceMaximumFlavour ("MaximumFlavour", "The maximum flavour allowed for the incoming quarks", &MEPP2GammaGammaPowheg::maxflavour_, 5, 1, 5, false, false, Interface::limited); static Parameter interfaceAlphaS ("AlphaS", "The value of alphaS to use if using a fixed alphaS", &MEPP2GammaGammaPowheg::alphaS_, 0.118, 0.0, 0.2, false, false, Interface::limited); static Switch interfaceFixedAlphaS ("FixedAlphaS", "Use a fixed value of alphaS", &MEPP2GammaGammaPowheg::fixedAlphaS_, false, false, false); static SwitchOption interfaceFixedAlphaSYes (interfaceFixedAlphaS, "Yes", "Use a fixed alphaS", true); static SwitchOption interfaceFixedAlphaSNo (interfaceFixedAlphaS, "No", "Use a running alphaS", false); static Switch interfaceSupressionFunction ("SupressionFunction", "Choice of the supression function", &MEPP2GammaGammaPowheg::supressionFunction_, 0, false, false); static SwitchOption interfaceSupressionFunctionNone (interfaceSupressionFunction, "None", "Default POWHEG approach", 0); static SwitchOption interfaceSupressionFunctionThetaFunction (interfaceSupressionFunction, "ThetaFunction", "Use theta functions at scale Lambda", 1); static SwitchOption interfaceSupressionFunctionSmooth (interfaceSupressionFunction, "Smooth", "Supress high pT by pt^2/(pt^2+lambda^2)", 2); static Parameter interfaceSupressionScale ("SupressionScale", "The square of the scale for the supression function", &MEPP2GammaGammaPowheg::lambda_, GeV, 20.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); static Switch interfaceSupressionScaleChoice ("SupressionScaleChoice", "Choice of the supression scale", &MEPP2GammaGammaPowheg::supressionScale_, 0, false, false); static SwitchOption interfaceSupressionScaleChoiceFixed (interfaceSupressionScaleChoice, "Fixed", "Use a fixed scale", 0); static SwitchOption interfaceSupressionScaleChoiceVariable (interfaceSupressionScaleChoice, "Variable", "Use the pT of the hard process as the scale", 1); static Reference interfaceShowerAlphaQCD ("ShowerAlphaQCD", "Reference to the object calculating the QCD coupling for the shower", &MEPP2GammaGammaPowheg::alphaQCD_, false, false, true, false, false); static Reference interfaceShowerAlphaQED ("ShowerAlphaQED", "Reference to the object calculating the QED coupling for the shower", &MEPP2GammaGammaPowheg::alphaQED_, false, false, true, false, false); static Parameter interfacepreQCDqqbarq ("preQCDqqbarq", "The constant for the Sudakov overestimate for the " "q qbar -> V Gamma +g with emission from the q", &MEPP2GammaGammaPowheg::preQCDqqbarq_, 23.0, 0.0, 1000.0, false, false, Interface::limited); static Parameter interfacepreQCDqqbarqbar ("preQCDqqbarqbar", "The constant for the Sudakov overestimate for the " "q qbar -> V Gamma +g with emission from the qbar", &MEPP2GammaGammaPowheg::preQCDqqbarqbar_, 23.0, 0.0, 1000.0, false, false, Interface::limited); static Switch interfaceScaleChoice ("ScaleChoice", "The scale choice to use", &MEPP2GammaGammaPowheg::scaleChoice_, 0, false, false); static SwitchOption interfaceScaleChoicepT (interfaceScaleChoice, "pT", "Use the pT of the photons", 0); static SwitchOption interfaceScaleChoiceMGammaGamma (interfaceScaleChoice, "MGammaGamma", "Use the mass of the photon pair", 1); static Parameter interfaceScalePreFactor ("ScalePreFactor", "Prefactor to change factorization/renormalisation scale", &MEPP2GammaGammaPowheg::scalePreFactor_, 1.0, 0.1, 10.0, false, false, Interface::limited); // prefactor_.push_back(preQCDqg_); // prefactor_.push_back(preQCDgqbar_); // prefactor_.push_back(preQEDqqbarq_); // prefactor_.push_back(preQEDqqbarqbar_); // prefactor_.push_back(preQEDqgq_); // prefactor_.push_back(preQEDgqbarqbar_); } double MEPP2GammaGammaPowheg::NLOWeight() const { // if leading-order return if(contrib_==0) return loME_; // prefactors CFfact_ = 4./3.*alphaS_/Constants::twopi; TRfact_ = 1./2.*alphaS_/Constants::twopi; // scale Energy2 mu2 = scale(); // virtual pieces double virt = CFfact_*subtractedVirtual(); // extract the partons and stuff for the real emission // and collinear counter terms // hadrons pair hadrons= make_pair(dynamic_ptr_cast(lastParticles().first->dataPtr() ), dynamic_ptr_cast(lastParticles().second->dataPtr())); // momentum fractions pair x = make_pair(lastX1(),lastX2()); // partons pair partons = make_pair(mePartonData()[0],mePartonData()[1]); // If necessary swap the particle data objects so that // first beam gives the incoming quark if(lastPartons().first ->dataPtr()!=partons.first) { swap(x.first,x.second); swap(hadrons.first,hadrons.second); } // convert the values of z tilde to z pair z; pair zJac; double rhomax(pow(1.-x.first,1.-power_)); double rho = zTilde_*rhomax; z.first = 1.-pow(rho,1./(1.-power_)); zJac.first = rhomax*pow(1.-z.first,power_)/(1.-power_); rhomax = pow(1.-x.second,1.-power_); rho = zTilde_*rhomax; z.second = 1.-pow(rho,1./(1.-power_)); zJac.second = rhomax*pow(1.-z.second,power_)/(1.-power_); // calculate the PDFs pair oldqPDF = make_pair(hadrons.first ->pdf()->xfx(hadrons.first ,partons.first ,scale(), x.first )/x.first , hadrons.second->pdf()->xfx(hadrons.second,partons.second,scale(), x.second)/x.second); // real/coll q/qbar pair newqPDF = make_pair(hadrons.first ->pdf()->xfx(hadrons.first ,partons.first ,scale(), x.first /z.first )*z.first /x.first , hadrons.second->pdf()->xfx(hadrons.second,partons.second,scale(), x.second/z.second)*z.second/x.second); // real/coll gluon pair newgPDF = make_pair(hadrons.first ->pdf()->xfx(hadrons.first ,gluon_,scale(), x.first /z.first )*z.first /x.first , hadrons.second->pdf()->xfx(hadrons.second,gluon_,scale(), x.second/z.second)*z.second/x.second); // coll terms // g -> q double collGQ = collinearGluon(mu2,zJac.first,z.first, oldqPDF.first,newgPDF.first); // g -> qbar double collGQbar = collinearGluon(mu2,zJac.second,z.second, oldqPDF.second,newgPDF.second); // q -> q double collQQ = collinearQuark(x.first ,mu2,zJac.first ,z.first , oldqPDF.first ,newqPDF.first ); // qbar -> qbar double collQbarQbar = collinearQuark(x.second,mu2,zJac.second,z.second, oldqPDF.second,newqPDF.second); // collinear remnants double coll = collQQ+collQbarQbar+collGQ+collGQbar; // real emission contribution double real1 = subtractedReal(x,z. first,zJac. first,oldqPDF. first, newqPDF. first,newgPDF. first, true); double real2 = subtractedReal(x,z.second,zJac.second,oldqPDF.second, newqPDF.second,newgPDF.second,false); // the total weight double wgt = loME_ + loME_*virt + loME_*coll + real1 + real2; return contrib_ == 1 ? max(0.,wgt) : max(0.,-wgt); } double MEPP2GammaGammaPowheg::loGammaGammaME(const cPDVector & particles, const vector & momenta, bool first) const { double output(0.); // analytic formula for speed if(!first) { Energy2 th = (momenta[0]-momenta[2]).m2(); Energy2 uh = (momenta[0]-momenta[3]).m2(); output = 4./3.*Constants::pi*SM().alphaEM(ZERO)*(th/uh+uh/th)* pow(double(particles[0]->iCharge())/3.,4); } // HE code result else { // wavefunctions for the incoming fermions SpinorWaveFunction em_in( momenta[0],particles[0],incoming); SpinorBarWaveFunction ep_in( momenta[1],particles[1],incoming); // wavefunctions for the outgoing bosons VectorWaveFunction v1_out(momenta[2],particles[2],outgoing); VectorWaveFunction v2_out(momenta[3],particles[3],outgoing); vector f1; vector a1; vector v1,v2; // calculate the wavefunctions for(unsigned int ix=0;ix<2;++ix) { em_in.reset(ix); f1.push_back(em_in); ep_in.reset(ix); a1.push_back(ep_in); v1_out.reset(2*ix); v1.push_back(v1_out); v2_out.reset(2*ix); v2.push_back(v2_out); } vector me(4,0.0); me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1,PDT::Spin1)); vector diag(2,0.0); SpinorWaveFunction inter; for(unsigned int ihel1=0;ihel1<2;++ihel1) { for(unsigned int ihel2=0;ihel2<2;++ihel2) { for(unsigned int ohel1=0;ohel1<2;++ohel1) { for(unsigned int ohel2=0;ohel2<2;++ohel2) { inter = FFPvertex_->evaluate(ZERO,5,f1[ihel1].particle(), f1[ihel1],v1[ohel1]); diag[0] = FFPvertex_->evaluate(ZERO,inter,a1[ihel2],v2[ohel2]); inter = FFPvertex_->evaluate(ZERO,5,f1[ihel1].particle(), f1[ihel1] ,v2[ohel2]); diag[1] = FFPvertex_->evaluate(ZERO,inter,a1[ihel2],v1[ohel1]); // individual diagrams for (size_t ii=0; ii<2; ++ii) me[ii] += std::norm(diag[ii]); // full matrix element diag[0] += diag[1]; output += std::norm(diag[0]); // storage of the matrix element for spin correlations me_(ihel1,ihel2,2*ohel1,2*ohel2) = diag[0]; } } } } // store diagram info, etc. DVector save(3); for (size_t i = 0; i < 3; ++i) save[i] = 0.25 * me[i]; meInfo(save); // spin and colour factors output *= 0.125/3./norm(FFPvertex_->norm()); } return output; } double MEPP2GammaGammaPowheg::loGammaqME(const cPDVector & particles, const vector & momenta, bool first) const { double output(0.); // analytic formula for speed if(!first) { Energy2 sh = (momenta[0]+momenta[1]).m2(); Energy2 th = (momenta[0]-momenta[2]).m2(); Energy2 uh = (momenta[0]-momenta[3]).m2(); output = -1./3./sh/th*(sh*sh+th*th+2.*uh*(sh+th+uh))* 4.*Constants::pi*SM().alphaEM(ZERO)* sqr(particles[0]->iCharge()/3.); } // HE result else { vector fin; vector gin; vector fout; vector vout; SpinorWaveFunction qin (momenta[0],particles[0],incoming); VectorWaveFunction glin(momenta[1],particles[1],incoming); VectorWaveFunction wout(momenta[2],particles[2],outgoing); SpinorBarWaveFunction qout(momenta[3],particles[3],outgoing); // polarization states for the particles for(unsigned int ix=0;ix<2;++ix) { qin.reset(ix) ; fin.push_back(qin); qout.reset(ix); fout.push_back(qout); glin.reset(2*ix); gin.push_back(glin); wout.reset(2*ix); vout.push_back(wout); } me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1, PDT::Spin1,PDT::Spin1Half)); // compute the matrix elements double me[3]={0.,0.,0.}; Complex diag[2]; SpinorWaveFunction inters; SpinorBarWaveFunction interb; for(unsigned int ihel1=0;ihel1<2;++ihel1) { for(unsigned int ihel2=0;ihel2<2;++ihel2) { for(unsigned int ohel1=0;ohel1<2;++ohel1) { // intermediates for the diagrams interb= FFGvertex_->evaluate(scale(),5,particles[3], fout[ohel1],gin[ihel2]); inters= FFGvertex_->evaluate(scale(),5,particles[0], fin[ihel1],gin[ihel2]); for(unsigned int vhel=0;vhel<2;++vhel) { diag[0] = FFPvertex_->evaluate(ZERO,fin[ihel1],interb,vout[vhel]); diag[1] = FFPvertex_->evaluate(ZERO,inters,fout[ohel1],vout[vhel]); // diagram contributions me[1] += norm(diag[0]); me[2] += norm(diag[1]); // total diag[0] += diag[1]; me[0] += norm(diag[0]); me_(ihel1,2*ihel2,2*vhel,ohel1) = diag[0]; } } } } // results // initial state spin and colour average double colspin = 1./24./4.; // and C_F N_c from matrix element colspin *= 4.; DVector save; for(unsigned int ix=0;ix<3;++ix) { me[ix] *= colspin; if(ix>0) save.push_back(me[ix]); } meInfo(save); output = me[0]/norm(FFGvertex_->norm()); } return output; } double MEPP2GammaGammaPowheg::loGammaqbarME(const cPDVector & particles, const vector & momenta, bool first) const { double output(0.); // analytic formula for speed if(!first) { Energy2 sh = (momenta[0]+momenta[1]).m2(); Energy2 uh = (momenta[0]-momenta[2]).m2(); Energy2 th = (momenta[0]-momenta[3]).m2(); output = -1./3./sh/th*(sh*sh+th*th+2.*uh*(sh+th+uh))* 4.*Constants::pi*SM().alphaEM()* sqr(particles[1]->iCharge()/3.); } // HE result else { vector ain; vector gin; vector aout; vector vout; VectorWaveFunction glin (momenta[0],particles[0],incoming); SpinorBarWaveFunction qbin (momenta[1],particles[1],incoming); VectorWaveFunction wout (momenta[2],particles[2],outgoing); SpinorWaveFunction qbout(momenta[3],particles[3],outgoing); // polarization states for the particles for(unsigned int ix=0;ix<2;++ix) { qbin .reset(ix ); ain .push_back(qbin ); qbout.reset(ix ); aout.push_back(qbout); glin.reset(2*ix); gin.push_back(glin); wout.reset(2*ix); vout.push_back(wout); } // if calculation spin corrections construct the me me_.reset(ProductionMatrixElement(PDT::Spin1,PDT::Spin1Half, PDT::Spin1,PDT::Spin1Half)); // compute the matrix elements double me[3]={0.,0.,0.}; Complex diag[2]; SpinorWaveFunction inters; SpinorBarWaveFunction interb; for(unsigned int ihel1=0;ihel1<2;++ihel1) { for(unsigned int ihel2=0;ihel2<2;++ihel2) { for(unsigned int ohel1=0;ohel1<2;++ohel1) { // intermediates for the diagrams inters= FFGvertex_->evaluate(scale(),5,particles[3], aout[ohel1],gin[ihel1]); interb= FFGvertex_->evaluate(scale(),5,particles[1], ain[ihel2],gin[ihel1]); for(unsigned int vhel=0;vhel<2;++vhel) { diag[0]= FFPvertex_->evaluate(ZERO,inters,ain[ihel2],vout[vhel]); diag[1]= FFPvertex_->evaluate(ZERO,aout[ohel1],interb,vout[vhel]); // diagram contributions me[1] += norm(diag[0]); me[2] += norm(diag[1]); // total diag[0] += diag[1]; me[0] += norm(diag[0]); me_(2*ihel1,ihel2,2*vhel,ohel1) = diag[0]; } } } } // results // initial state spin and colour average double colspin = 1./24./4.; // and C_F N_c from matrix element colspin *= 4.; DVector save; for(unsigned int ix=0;ix<3;++ix) { me[ix] *= colspin; if(ix>0) save.push_back(me[ix]); } meInfo(save); output = me[0]/norm(FFGvertex_->norm()); } return output; } double MEPP2GammaGammaPowheg::loGammagME(const cPDVector & particles, const vector & momenta, bool first) const { double output(0.); // analytic formula for speed if(!first) { Energy2 uh = (momenta[0]-momenta[2]).m2(); Energy2 th = (momenta[0]-momenta[3]).m2(); output = 8./9.*double((th*th+uh*uh)/uh/th)* 4.*Constants::pi*SM().alphaEM(ZERO)* sqr(particles[0]->iCharge()/3.); } else { vector fin; vector ain; vector gout; vector vout; SpinorWaveFunction qin (momenta[0],particles[0],incoming); SpinorBarWaveFunction qbin(momenta[1],particles[1],incoming); VectorWaveFunction wout(momenta[2],particles[2],outgoing); VectorWaveFunction glout(momenta[3],particles[3],outgoing); // polarization states for the particles for(unsigned int ix=0;ix<2;++ix) { qin.reset(ix) ; fin.push_back(qin); qbin.reset(ix) ; ain.push_back(qbin); glout.reset(2*ix); gout.push_back(glout); wout.reset(2*ix); vout.push_back(wout); } // if calculation spin corrections construct the me if(first) me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1,PDT::Spin1)); // compute the matrix elements double me[3]={0.,0.,0.}; Complex diag[2]; SpinorWaveFunction inters; SpinorBarWaveFunction interb; for(unsigned int ihel1=0;ihel1<2;++ihel1) { for(unsigned int ihel2=0;ihel2<2;++ihel2) { for(unsigned int ohel1=0;ohel1<2;++ohel1) { // intermediates for the diagrams inters= FFGvertex_->evaluate(scale(),5,particles[0], fin[ihel1],gout[ohel1]); interb= FFGvertex_->evaluate(scale(),5,particles[1], ain[ihel2],gout[ohel1]); for(unsigned int vhel=0;vhel<2;++vhel) { diag[0]= FFPvertex_->evaluate(ZERO,fin[ihel1],interb,vout[vhel]); diag[1]= FFPvertex_->evaluate(ZERO,inters,ain[ihel2],vout[vhel]); // diagram contributions me[1] += norm(diag[0]); me[2] += norm(diag[1]); // total diag[0] += diag[1]; me[0] += norm(diag[0]); if(first) me_(ihel1,ihel2,vhel,2*ohel1) = diag[0]; } } } } // results // initial state spin and colour average double colspin = 1./9./4.; // and C_F N_c from matrix element colspin *= 4.; DVector save; for(unsigned int ix=0;ix<3;++ix) { me[ix] *= colspin; if(ix>0) save.push_back(me[ix]); } meInfo(save); output = me[0]/norm(FFGvertex_->norm()); } return output; } InvEnergy2 MEPP2GammaGammaPowheg:: realGammaGammagME(const cPDVector & particles, const vector & momenta, DipoleType dipole, RadiationType rad, bool ) const { // matrix element double sum = realME(particles,momenta); // loop over the QCD and QCD dipoles InvEnergy2 dipoles[2]; pair supress[2]; // compute the two dipole terms unsigned int iemit = 4, ihard = 3; double x = (momenta[0]*momenta[1]-momenta[iemit]*momenta[1]- momenta[iemit]*momenta[0])/(momenta[0]*momenta[1]); Lorentz5Momentum Kt = momenta[0]+momenta[1]-momenta[iemit]; vector pa(4),pb(4); // momenta for q -> q g/gamma emission pa[0] = x*momenta[0]; pa[1] = momenta[1]; Lorentz5Momentum K = pa[0]+pa[1]; Lorentz5Momentum Ksum = K+Kt; Energy2 K2 = K.m2(); Energy2 Ksum2 = Ksum.m2(); pa[2] = momenta[2]-2.*Ksum*(Ksum*momenta[2])/Ksum2+2*K*(Kt*momenta[2])/K2; pa[2].setMass(momenta[2].mass()); pa[3] = momenta[ihard] -2.*Ksum*(Ksum*momenta[ihard])/Ksum2+2*K*(Kt*momenta[ihard])/K2; pa[3].setMass(ZERO); cPDVector part(particles.begin(),--particles.end()); part[3] = particles[ihard]; // first leading-order matrix element double lo1 = loGammaGammaME(part,pa); // first dipole dipoles[0] = 1./(momenta[0]*momenta[iemit])/x*(2./(1.-x)-(1.+x))*lo1; supress[0] = supressionFunction(momenta[iemit].perp(),pa[3].perp()); // momenta for qbar -> qbar g/gamma emission pb[0] = momenta[0]; pb[1] = x*momenta[1]; K = pb[0]+pb[1]; Ksum = K+Kt; K2 = K.m2(); Ksum2 = Ksum.m2(); pb[2] = momenta[2]-2.*Ksum*(Ksum*momenta[2])/Ksum2+2*K*(Kt*momenta[2])/K2; pb[2].setMass(momenta[2].mass()); pb[3] = momenta[ihard] -2.*Ksum*(Ksum*momenta[ihard])/Ksum2+2*K*(Kt*momenta[ihard])/K2; pb[3].setMass(ZERO); // second LO matrix element double lo2 = loGammaGammaME(part,pb); // second dipole dipoles[1] = 1./(momenta[1]*momenta[iemit])/x*(2./(1.-x)-(1.+x))*lo2; supress[1] = supressionFunction(momenta[iemit].perp(),pb[3].perp()); for(unsigned int ix=0;ix<2;++ix) dipoles[ix] *= 4./3.; // denominator for the matrix element InvEnergy2 denom = abs(dipoles[0]) + abs(dipoles[1]); // contribution if( denom==ZERO || dipoles[(dipole-1)/2]==ZERO ) return ZERO; sum *= abs(dipoles[(dipole-1)/2])/denom; // final coupling factors InvEnergy2 output; if(rad==Subtraction) { output = alphaS_*alphaEM_* (sum*UnitRemoval::InvE2*supress[(dipole-1)/2].first - dipoles[(dipole-1)/2]); } else { output = alphaS_*alphaEM_*sum*UnitRemoval::InvE2; if(rad==Hard) output *=supress[(dipole-1)/2].second; else if(rad==Shower) output *=supress[(dipole-1)/2].first ; } return output; } InvEnergy2 MEPP2GammaGammaPowheg::realGammaGammaqME(const cPDVector & particles, const vector & momenta, DipoleType dipole, RadiationType rad, bool ) const { double sum = realME(particles,momenta); // initial-state QCD dipole double x = (momenta[0]*momenta[1]-momenta[4]*momenta[1]- momenta[4]*momenta[0])/(momenta[0]*momenta[1]); Lorentz5Momentum Kt = momenta[0]+momenta[1]-momenta[4]; vector pa(4); pa[0] = momenta[0]; pa[1] = x*momenta[1]; Lorentz5Momentum K = pa[0]+pa[1]; Lorentz5Momentum Ksum = K+Kt; Energy2 K2 = K.m2(); Energy2 Ksum2 = Ksum.m2(); pa[2] = momenta[2]-2.*Ksum*(Ksum*momenta[2])/Ksum2+2*K*(Kt*momenta[2])/K2; pa[2].setMass(momenta[2].mass()); pa[3] = momenta[3] -2.*Ksum*(Ksum*momenta[3])/Ksum2+2*K*(Kt*momenta[3])/K2; pa[3].setMass(ZERO); cPDVector part(particles.begin(),--particles.end()); part[1] = particles[4]->CC(); double lo1 = loGammaGammaME(part,pa); InvEnergy2 D1 = 0.5/(momenta[1]*momenta[4])/x*(1.-2.*x*(1.-x))*lo1; // initial-final QED dipole vector pb(4); x = 1.-(momenta[3]*momenta[4])/(momenta[4]*momenta[0]+momenta[0]*momenta[3]); pb[3] = momenta[4]+momenta[3]-(1.-x)*momenta[0]; pb[0] = x*momenta[0]; pb[1] = momenta[1]; pb[2] = momenta[2]; double z = momenta[0]*momenta[3]/(momenta[0]*momenta[3]+momenta[0]*momenta[4]); part[1] = particles[1]; part[3] = particles[4]; double lo2 = loGammaqME(part,pb); Energy pT = sqrt(-(pb[0]-pb[3]).m2()*(1.-x)*(1.-z)*z/x); InvEnergy2 DF = 1./(momenta[4]*momenta[3])/x*(1./(1.-x+z)-2.+z)*lo2; InvEnergy2 DI = 1./(momenta[0]*momenta[3])/x*(1./(1.-x+z)-1.-x)*lo2; DI *= sqr(double(particles[0]->iCharge())/3.); DF *= sqr(double(particles[0]->iCharge())/3.); InvEnergy2 denom = abs(D1)+abs(DI)+abs(DF); pair supress; InvEnergy2 term; if ( dipole == IFQED1 ) { term = DI; supress = supressionFunction(pT,pb[3].perp()); } else if( dipole == FIQED1 ) { term = DF; supress = supressionFunction(pT,pb[3].perp()); } else { term = D1; supress = supressionFunction(momenta[4].perp(),pa[3].perp()); } if( denom==ZERO || term == ZERO ) return ZERO; sum *= abs(term)/denom; // final coupling factors InvEnergy2 output; if(rad==Subtraction) { output = alphaS_*alphaEM_* (sum*UnitRemoval::InvE2*supress.first - term); } else { output = alphaS_*alphaEM_*sum*UnitRemoval::InvE2; if(rad==Hard) output *= supress.second; else if(rad==Shower) output *= supress.first ; } // final coupling factors return output; } InvEnergy2 MEPP2GammaGammaPowheg:: realGammaGammaqbarME(const cPDVector & particles, const vector & momenta, DipoleType dipole, RadiationType rad, bool) const { double sum = realME(particles,momenta); // initial-state QCD dipole double x = (momenta[0]*momenta[1]-momenta[4]*momenta[1]-momenta[4]*momenta[0])/ (momenta[0]*momenta[1]); Lorentz5Momentum Kt = momenta[0]+momenta[1]-momenta[4]; vector pa(4); pa[0] = x*momenta[0]; pa[1] = momenta[1]; Lorentz5Momentum K = pa[0]+pa[1]; Lorentz5Momentum Ksum = K+Kt; Energy2 K2 = K.m2(); Energy2 Ksum2 = Ksum.m2(); pa[2] = momenta[2]-2.*Ksum*(Ksum*momenta[2])/Ksum2+2*K*(Kt*momenta[2])/K2; pa[2].setMass(momenta[2].mass()); pa[3] = momenta[3] -2.*Ksum*(Ksum*momenta[3])/Ksum2+2*K*(Kt*momenta[3])/K2; pa[3].setMass(ZERO); cPDVector part(particles.begin(),--particles.end()); part[0] = particles[4]->CC(); double lo1 = loGammaGammaME(part,pa); InvEnergy2 D1 = 0.5/(momenta[0]*momenta[4])/x*(1.-2.*x*(1.-x))*lo1; // initial-final QED dipole vector pb(4); x = 1.-(momenta[3]*momenta[4])/(momenta[4]*momenta[1]+momenta[1]*momenta[3]); pb[3] = momenta[4]+momenta[3]-(1.-x)*momenta[1]; pb[0] = momenta[0]; pb[1] = x*momenta[1]; pb[2] = momenta[2]; double z = momenta[1]*momenta[3]/(momenta[1]*momenta[3]+momenta[1]*momenta[4]); part[0] = particles[0]; part[3] = particles[4]; double lo2 = loGammaqbarME(part,pb); Energy pT = sqrt(-(pb[1]-pb[3]).m2()*(1.-x)*(1.-z)*z/x); InvEnergy2 DF = 1./(momenta[4]*momenta[3])/x*(2./(1.-x+z)-2.+z)*lo2; InvEnergy2 DI = 1./(momenta[0]*momenta[3])/x*(2./(1.-x+z)-1.-x)*lo2; InvEnergy2 term; DI *= sqr(double(particles[1]->iCharge())/3.); DF *= sqr(double(particles[1]->iCharge())/3.); InvEnergy2 denom = abs(D1)+abs(DI)+abs(DF); pair supress; if ( dipole == IFQED2 ) { term = DI; supress = supressionFunction(pT,pb[3].perp()); } else if( dipole == FIQED2 ) { term = DF; supress = supressionFunction(pT,pb[3].perp()); } else { term = D1; supress = supressionFunction(momenta[4].perp(),pa[3].perp()); } if( denom==ZERO || dipole==ZERO ) return ZERO; sum *= abs(term)/denom; // final coupling factors InvEnergy2 output; if(rad==Subtraction) { output = alphaS_*alphaEM_* (sum*UnitRemoval::InvE2*supress.first - term); } else { output = alphaS_*alphaEM_*sum*UnitRemoval::InvE2; if(rad==Hard) output *= supress.second; else if(rad==Shower) output *= supress.first ; } // final coupling factors return output; } double MEPP2GammaGammaPowheg:: realME(const cPDVector & particles, const vector & momenta) const { vector qin; vector qbarin; vector wout,pout,gout; SpinorWaveFunction q_in; SpinorBarWaveFunction qbar_in; VectorWaveFunction g_out; VectorWaveFunction v_out (momenta[2],particles[2],outgoing); VectorWaveFunction p_out (momenta[3],particles[3],outgoing); // q qbar -> gamma gamma g if(particles[4]->id()==ParticleID::g) { q_in = SpinorWaveFunction (momenta[0],particles[0],incoming); qbar_in = SpinorBarWaveFunction (momenta[1],particles[1],incoming); g_out = VectorWaveFunction (momenta[4],particles[4],outgoing); } // q g -> gamma gamma q else if(particles[4]->id()>0) { q_in = SpinorWaveFunction (momenta[0],particles[0],incoming); qbar_in = SpinorBarWaveFunction (momenta[4],particles[4],outgoing); g_out = VectorWaveFunction (momenta[1],particles[1],incoming); } else if(particles[4]->id()<0) { q_in = SpinorWaveFunction (momenta[4],particles[4],outgoing); qbar_in = SpinorBarWaveFunction (momenta[1],particles[1],incoming); g_out = VectorWaveFunction (momenta[0],particles[0],incoming); } else assert(false); for(unsigned int ix=0;ix<2;++ix) { q_in.reset(ix); qin.push_back(q_in); qbar_in.reset(ix); qbarin.push_back(qbar_in); g_out.reset(2*ix); gout.push_back(g_out); p_out.reset(2*ix); pout.push_back(p_out); v_out.reset(2*ix); wout.push_back(v_out); } vector diag(6 , 0.); Energy2 mu2 = scale(); double sum(0.); for(unsigned int ihel1=0;ihel1<2;++ihel1) { for(unsigned int ihel2=0;ihel2<2;++ihel2) { for(unsigned int whel=0;whel<2;++whel) { for(unsigned int phel=0;phel<2;++phel) { for(unsigned int ghel=0;ghel<2;++ghel) { // first diagram SpinorWaveFunction inters1 = FFPvertex_->evaluate(ZERO,5,qin[ihel1].particle(),qin[ihel1],pout[phel]); SpinorBarWaveFunction inters2 = FFPvertex_->evaluate(ZERO,5,qin[ihel1].particle()->CC(), qbarin[ihel2],wout[whel]); diag[0] = FFGvertex_->evaluate(mu2,inters1,inters2,gout[ghel]); // second diagram SpinorWaveFunction inters3 = FFGvertex_->evaluate(mu2,5,qin[ihel1].particle(),qin[ihel1],gout[ghel]); SpinorBarWaveFunction inters4 = FFPvertex_->evaluate(ZERO,5,qbarin[ihel2].particle(), qbarin[ihel2],pout[phel]); diag[1] = FFPvertex_->evaluate(ZERO,inters3,inters4,wout[whel]); // fourth diagram diag[2] = FFPvertex_->evaluate(ZERO,inters3,inters2,pout[phel]); // fifth diagram SpinorBarWaveFunction inters5 = FFGvertex_->evaluate(mu2,5,qbarin[ihel2].particle(), qbarin[ihel2],gout[ghel]); diag[3] = FFPvertex_->evaluate(ZERO,inters1,inters5,wout[whel]); // sixth diagram SpinorWaveFunction inters6 = FFPvertex_->evaluate(ZERO,5,qbarin[ihel2].particle()->CC(), qin[ihel1],wout[whel]); diag[4] = FFGvertex_->evaluate(mu2,inters6,inters4,gout[ghel]); // eighth diagram diag[5] = FFPvertex_->evaluate(ZERO,inters6,inters5,pout[phel]); // sum Complex dsum = std::accumulate(diag.begin(),diag.end(),Complex(0.)); sum += norm(dsum); } } } } } // divide out the em and strong couplings sum /= norm(FFGvertex_->norm()*FFPvertex_->norm()); // final spin and colour factors spin = 1/4 colour = 4/9 if(particles[4]->id()==ParticleID::g) sum /= 9.; // final spin and colour factors spin = 1/4 colour = 4/(3*8) else sum /= 24.; // finally identical particle factor return 0.5*sum; } double MEPP2GammaGammaPowheg::subtractedVirtual() const { double v = 1+tHat()/sHat(); double born = (1-v)/v+v/(1-v); double finite_term = born* (2./3.*sqr(Constants::pi)-3.+sqr(log(v))+sqr(log(1-v))+3.*log(1-v))+ 2.+2.*log(v)+2.*log(1-v)+3.*(1-v)/v*(log(v)-log(1-v))+ (2.+v/(1-v))*sqr(log(v))+(2.+(1-v)/v)*sqr(log(1-v)); double virt = ((6.-(2./3.)*sqr(Constants::pi))* born-2.+finite_term); return virt/born; } double MEPP2GammaGammaPowheg::subtractedReal(pair x, double z, double zJac, double oldqPDF, double newqPDF, double newgPDF,bool order) const { double vt = vTilde_*(1.-z); double vJac = 1.-z; Energy pT = sqrt(sHat()*vt*(1.-vt-z)/z); // rapidities double rapidity; if(order) { rapidity = -log(x.second*sqrt(lastS())/pT*vt); } else { rapidity = log(x.first *sqrt(lastS())/pT*vt); } // CMS system Energy rs=sqrt(lastS()); Lorentz5Momentum pcmf = Lorentz5Momentum(ZERO,ZERO,0.5*rs*(x.first-x.second), 0.5*rs*(x.first+x.second)); pcmf.rescaleMass(); Boost blab(pcmf.boostVector()); // emission from the quark radiation vector pnew(5); if(order) { pnew [0] = Lorentz5Momentum(ZERO,ZERO,0.5*rs*x.first/z, 0.5*rs*x.first/z,ZERO); pnew [1] = Lorentz5Momentum(ZERO,ZERO,-0.5*rs*x.second, 0.5*rs*x.second,ZERO) ; } else { pnew[0] = Lorentz5Momentum(ZERO,ZERO,0.5*rs*x.first, 0.5*rs*x.first,ZERO); pnew[1] = Lorentz5Momentum(ZERO,ZERO,-0.5*rs*x.second/z, 0.5*rs*x.second/z,ZERO) ; } pnew [2] = meMomenta()[2]; pnew [3] = meMomenta()[3]; pnew [4] = Lorentz5Momentum(pT*cos(phi_),pT*sin(phi_), pT*sinh(rapidity), pT*cosh(rapidity), ZERO); Lorentz5Momentum K = pnew [0]+pnew [1]-pnew [4]; Lorentz5Momentum Kt = pcmf; Lorentz5Momentum Ksum = K+Kt; Energy2 K2 = K.m2(); Energy2 Ksum2 = Ksum.m2(); for(unsigned int ix=2;ix<4;++ix) { pnew [ix].boost(blab); pnew [ix] = pnew [ix] - 2.*Ksum*(Ksum*pnew [ix])/Ksum2 +2*K*(Kt*pnew [ix])/K2; } // phase-space prefactors double phase = zJac*vJac/z; // real emission q qbar vector output(4,0.); double realQQ(0.),realGQ(0.); if(!(zTilde_<1e-7 || vt<1e-7 || 1.-z-vt < 1e-7 )) { cPDVector particles(mePartonData()); particles.push_back(gluon_); // calculate the full 2->3 matrix element realQQ = sHat()*phase*newqPDF/oldqPDF* realGammaGammagME(particles,pnew,order ? IIQCD1 : IIQCD2,Subtraction,false); if(order) { particles[0] = gluon_; particles[4] = mePartonData()[0]->CC(); realGQ = sHat()*phase*newgPDF/oldqPDF* realGammaGammaqbarME(particles,pnew,IIQCD2,Subtraction,false); } else { particles[1] = gluon_; particles[4] = mePartonData()[1]->CC(); realGQ = sHat()*phase*newgPDF/oldqPDF* realGammaGammaqME (particles,pnew,IIQCD1,Subtraction,false); } } // return the answer return realQQ+realGQ; } double MEPP2GammaGammaPowheg::collinearQuark(double x, Energy2 mu2, double jac, double z, double oldPDF, double newPDF) const { if(1.-z < 1.e-8) return 0.; return CFfact_*( // this bit is multiplied by LO PDF sqr(Constants::pi)/3.-5.+2.*sqr(log(1.-x )) +(1.5+2.*log(1.-x ))*log(sHat()/mu2) // NLO PDF bit +jac /z * newPDF /oldPDF * (1.-z -(1.+z )*log(sqr(1.-z )/z ) -(1.+z )*log(sHat()/mu2)-2.*log(z )/(1.-z )) // + function bit +jac /z *(newPDF /oldPDF -z )* 2./(1.-z )*log(sHat()*sqr(1.-z )/mu2)); } double MEPP2GammaGammaPowheg::collinearGluon(Energy2 mu2, double jac, double z, double oldPDF, double newPDF) const { if(1.-z < 1.e-8) return 0.; return TRfact_*jac/z*newPDF/oldPDF* ((sqr(z)+sqr(1.-z))*log(sqr(1.-z)*sHat()/z/mu2) +2.*z*(1.-z)); } void MEPP2GammaGammaPowheg::doinit() { HwMEBase::doinit(); vector mopt(2,1); massOption(mopt); // get the vertices we need // get a pointer to the standard model object in the run static const tcHwSMPtr hwsm = dynamic_ptr_cast(standardModel()); if (!hwsm) throw InitException() << "hwsm pointer is null in" << " MEPP2GammaGamma::doinit()" << Exception::abortnow; // get pointers to all required Vertex objects FFPvertex_ = hwsm->vertexFFP(); FFGvertex_ = hwsm->vertexFFG(); gluon_ = getParticleData(ParticleID::g); // sampling factors prefactor_.push_back(preQCDqqbarq_); prefactor_.push_back(preQCDqqbarqbar_); prefactor_.push_back(preQCDqg_); prefactor_.push_back(preQCDgqbar_); prefactor_.push_back(preQEDqqbarq_); prefactor_.push_back(preQEDqqbarqbar_); prefactor_.push_back(preQEDqgq_); prefactor_.push_back(preQEDgqbarqbar_); } RealEmissionProcessPtr MEPP2GammaGammaPowheg:: generateHardest(RealEmissionProcessPtr born, ShowerInteraction inter) { beams_.clear(); partons_.clear(); bool QCDAllowed = inter !=ShowerInteraction::QED; bool QEDAllowed = inter !=ShowerInteraction::QCD; // find the incoming particles // and get the particles to be showered ParticleVector incoming,particlesToShower; pair x; for(unsigned int ix=0;ixbornIncoming().size();++ix) { incoming.push_back( born->bornIncoming()[ix] ); beams_.push_back(dynamic_ptr_cast(born->hadrons()[ix]->dataPtr())); partons_.push_back( born->bornIncoming()[ix]->dataPtr() ); particlesToShower.push_back( born->bornIncoming()[ix] ); if(ix==0) x.first = incoming.back()->momentum().rho()/born->hadrons()[ix]->momentum().rho(); else x.second = incoming.back()->momentum().rho()/born->hadrons()[ix]->momentum().rho(); } // find the parton which should be first if( ( particlesToShower[1]->id() > 0 && particlesToShower[0]->id() < 0 ) || ( particlesToShower[0]->id() == ParticleID::g && particlesToShower[1]->id() < 6 && particlesToShower[1]->id() > 0 ) ) { swap(particlesToShower[0],particlesToShower[1]); swap(partons_[0],partons_[1]); swap(beams_ [0],beams_ [1]); swap(x.first ,x.second ); } // check that quark is along +ve z direction quarkplus_ = particlesToShower[0]->momentum().z() > ZERO; // outgoing partons for(unsigned int ix=0;ixbornOutgoing().size();++ix) { particlesToShower.push_back( born->bornOutgoing()[ix] ); } if(particlesToShower.size()!=4) return RealEmissionProcessPtr(); if(particlesToShower[2]->id()!=ParticleID::gamma) swap(particlesToShower[2],particlesToShower[3]); if(particlesToShower[3]->id()==ParticleID::gamma) { if(QCDAllowed) return hardQCDEmission(born,particlesToShower,x); } else { if(QEDAllowed) return hardQEDEmission(born,particlesToShower,x); } return born; } RealEmissionProcessPtr MEPP2GammaGammaPowheg:: hardQCDEmission(RealEmissionProcessPtr born, ParticleVector particlesToShower, pair x) { Energy rootS = sqrt(lastS()); // limits on the rapidity of the jet double minyj = -8.0,maxyj = 8.0; // generate the hard emission vector pT; Energy pTmax(-GeV); cPDVector selectedParticles; vector selectedMomenta; int iemit(-1); for(unsigned int ix=0;ix<4;++ix) { pT.push_back(0.5*generator()->maximumCMEnergy()); double a = alphaQCD_->overestimateValue()/Constants::twopi* prefactor_[ix]*(maxyj-minyj); cPDVector particles; for(unsigned int iy=0;iydataPtr()); if(ix<2) particles.push_back(gluon_); else if(ix==2) { particles.push_back(particles[0]->CC()); particles[0] = gluon_; } else { particles.push_back(particles[1]->CC()); particles[1] = gluon_; } vector momenta(5); do { pT[ix] *= pow(UseRandom::rnd(),1./a); double y = UseRandom::rnd()*(maxyj-minyj)+ minyj; double vt,z; if(ix%2==0) { vt = pT[ix]*exp(-y)/rootS/x.second; z = (1.-pT[ix]*exp(-y)/rootS/x.second)/(1.+pT[ix]*exp( y)/rootS/x.first ); if(z>1.||z1.||z1.-z || vt<0.) continue; if(ix%2==0) { momenta[0] = particlesToShower[0]->momentum()/z; momenta[1] = particlesToShower[1]->momentum(); } else { momenta[0] = particlesToShower[0]->momentum(); momenta[1] = particlesToShower[1]->momentum()/z; } double phi = Constants::twopi*UseRandom::rnd(); momenta[2] = particlesToShower[2]->momentum(); momenta[3] = particlesToShower[3]->momentum(); if(!quarkplus_) y *= -1.; momenta[4] = Lorentz5Momentum(pT[ix]*cos(phi),pT[ix]*sin(phi), pT[ix]*sinh(y),pT[ix]*cosh(y), ZERO); Lorentz5Momentum K = momenta[0] + momenta[1] - momenta[4]; Lorentz5Momentum Kt = momenta[2]+momenta[3]; Lorentz5Momentum Ksum = K+Kt; Energy2 K2 = K.m2(), Ksum2 = Ksum.m2(); for(unsigned int iy=2;iy<4;++iy) { momenta [iy] = momenta [iy] - 2.*Ksum*(Ksum*momenta [iy])/Ksum2 +2*K*(Kt*momenta [iy])/K2; } // matrix element piece double wgt = alphaQCD_->ratio(sqr(pT[ix]))*z/(1.-vt)/prefactor_[ix]/loME_; if(ix==0) wgt *= sqr(pT[ix])*realGammaGammagME(particles,momenta,IIQCD1,Shower,false); else if(ix==1) wgt *= sqr(pT[ix])*realGammaGammagME(particles,momenta,IIQCD2,Shower,false); else if(ix==2) wgt *= sqr(pT[ix])*realGammaGammaqbarME(particles,momenta,IIQCD1,Shower,false); else if(ix==3) wgt *= sqr(pT[ix])*realGammaGammaqME(particles,momenta,IIQCD2,Shower,false); wgt *= 4.*Constants::pi/alphaS_; // pdf piece double pdf[2]; if(ix%2==0) { pdf[0] = beams_[0]->pdf()->xfx(beams_[0],partons_ [0], scale(), x.first ) /x.first; pdf[1] = beams_[0]->pdf()->xfx(beams_[0],particles[0], scale()+sqr(pT[ix]),x.first /z)*z/x.first; } else { pdf[0] = beams_[1]->pdf()->xfx(beams_[1],partons_ [1], scale() ,x.second ) /x.second; pdf[1] = beams_[1]->pdf()->xfx(beams_[1],particles[1], scale()+sqr(pT[ix]),x.second/z)*z/x.second; } if(pdf[0]<=0.||pdf[1]<=0.) continue; wgt *= pdf[1]/pdf[0]; if(wgt>1.) generator()->log() << "Weight greater than one in " << "MEPP2GammaGammaPowheg::hardQCDEmission() " << "for channel " << ix << " Weight = " << wgt << "\n"; if(UseRandom::rnd()minpT_); if(pT[ix]>minpT_ && pT[ix]>pTmax) { pTmax = pT[ix]; selectedParticles = particles; selectedMomenta = momenta; iemit=ix; } } // if no emission if(pTmaxpT()[ShowerInteraction::QCD] = minpT_; return born; } // construct the HardTree object needed to perform the showers // create the partons ParticleVector newparticles; newparticles.push_back(selectedParticles[0]->produceParticle(selectedMomenta[0])); newparticles.push_back(selectedParticles[1]->produceParticle(selectedMomenta[1])); for(unsigned int ix=2;ixdataPtr()-> produceParticle(selectedMomenta[ix])); } newparticles.push_back(selectedParticles[4]->produceParticle(selectedMomenta[4])); // identify the type of process // gluon emission if(newparticles.back()->id()==ParticleID::g) { newparticles[4]->incomingColour(newparticles[0]); newparticles[4]->incomingColour(newparticles[1],true); } // quark else if(newparticles.back()->id()>0) { iemit=1; newparticles[4]->incomingColour(newparticles[1]); newparticles[1]-> colourConnect(newparticles[0]); } // antiquark else { iemit=0; newparticles[4]->incomingColour(newparticles[0],true); newparticles[1]-> colourConnect(newparticles[0]); } // add incoming int ispect = iemit==0 ? 1 : 0; if(particlesToShower[0]==born->bornIncoming()[0]) { born->incoming().push_back(newparticles[0]); born->incoming().push_back(newparticles[1]); } else { born->incoming().push_back(newparticles[1]); born->incoming().push_back(newparticles[0]); swap(iemit,ispect); } // add the outgoing for(unsigned int ix=2;ixoutgoing().push_back(newparticles[ix]); // emitter spectator etc born->emitter (iemit ); born->spectator(ispect); born->emitted(4); // x values pair xnew; for(unsigned int ix=0;ix<2;++ix) { double x = born->incoming()[ix]->momentum().rho()/born->hadrons()[ix]->momentum().rho(); if(ix==0) xnew.first = x; else xnew.second = x; } born->x(xnew); // max pT born->pT()[ShowerInteraction::QCD] = pTmax; born->interaction(ShowerInteraction::QCD); // return the process return born; } RealEmissionProcessPtr MEPP2GammaGammaPowheg:: hardQEDEmission(RealEmissionProcessPtr born, ParticleVector particlesToShower, pair x) { // return if not emission from quark if(particlesToShower[0]->id()!=ParticleID::g && particlesToShower[1]->id()!=ParticleID::g ) return RealEmissionProcessPtr(); // generate the hard emission vector pT; Energy pTmax(-GeV); cPDVector selectedParticles; vector selectedMomenta; int iemit(-1); pair mewgt(make_pair(0.,0.)); for(unsigned int ix=0;ixdataPtr()); selectedMomenta.push_back(particlesToShower[ix]->momentum()); } selectedParticles.push_back(getParticleData(ParticleID::gamma)); swap(selectedParticles[3],selectedParticles[4]); selectedMomenta.push_back(Lorentz5Momentum()); swap(selectedMomenta[3],selectedMomenta[4]); Lorentz5Momentum pin,pout; double xB; unsigned int iloc; if(particlesToShower[0]->dataPtr()->charged()) { pin = particlesToShower[0]->momentum(); xB = x.first; iloc = 6; } else { pin = particlesToShower[1]->momentum(); xB = x.second; iloc = 7; } pout = particlesToShower[3]->momentum(); Lorentz5Momentum q = pout-pin; Axis axis(q.vect().unit()); LorentzRotation rot; double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); rot = LorentzRotation(); if(axis.perp2()>1e-20) { rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); rot.rotateX(Constants::pi); } if(abs(1.-q.e()/q.vect().mag())>1e-6) rot.boostZ(q.e()/q.vect().mag()); Lorentz5Momentum ptemp = rot*pin; if(ptemp.perp2()/GeV2>1e-20) { Boost trans = -1./ptemp.e()*ptemp.vect(); trans.setZ(0.); rot.boost(trans); } rot.invert(); Energy Q = sqrt(-q.m2()); double xT = sqrt((1.-xB)/xB); double xTMin = 2.*minpT_/Q; double wgt(0.); double a = alphaQED_->overestimateValue()*prefactor_[iloc]/Constants::twopi; Lorentz5Momentum p1,p2,p3; do { wgt = 0.; // intergration variables dxT/xT^3 xT *= 1./sqrt(1.-2.*log(UseRandom::rnd())/a*sqr(xT)); // dz double zp = UseRandom::rnd(); double xp = 1./(1.+0.25*sqr(xT)/zp/(1.-zp)); if(xT1.) continue; // phase-space piece of the weight wgt = 4.*sqr(1.-xp)*(1.-zp)*zp/prefactor_[iloc]/loME_; // coupling Energy2 pT2 = 0.25*sqr(Q*xT); wgt *= alphaQED_->ratio(pT2); // matrix element wgt *= 4.*Constants::pi/alphaEM_; // PDF double pdf[2]; if(iloc==6) { pdf[0] = beams_[0]->pdf()-> xfx(beams_[0],partons_[0],scale() ,x.first ); pdf[1] = beams_[0]->pdf()-> xfx(beams_[0],partons_[0],scale()+pT2,x.first /xp); } else { pdf[0] = beams_[1]->pdf()-> xfx(beams_[1],partons_[1],scale() ,x.second ); pdf[1] = beams_[1]->pdf()-> xfx(beams_[1],partons_[1],scale()+pT2,x.second/xp); } if(pdf[0]<=0.||pdf[1]<=0.) { wgt = 0.; continue; } wgt *= pdf[1]/pdf[0]; // matrix element piece double phi = Constants::twopi*UseRandom::rnd(); double x2 = 1.-(1.-zp)/xp; double x3 = 2.-1./xp-x2; p1=Lorentz5Momentum(ZERO,ZERO,0.5*Q/xp,0.5*Q/xp,ZERO); p2=Lorentz5Momentum( 0.5*Q*xT*cos(phi), 0.5*Q*xT*sin(phi), -0.5*Q*x2,0.5*Q*sqrt(sqr(xT)+sqr(x2))); p3=Lorentz5Momentum(-0.5*Q*xT*cos(phi),-0.5*Q*xT*sin(phi), -0.5*Q*x3,0.5*Q*sqrt(sqr(xT)+sqr(x3))); selectedMomenta[iloc-6] = rot*p1; selectedMomenta[3] = rot*p3; selectedMomenta[4] = rot*p2; if(iloc==6) { mewgt.first = sqr(Q)*realGammaGammaqME(selectedParticles,selectedMomenta,IFQED1,Shower,false); mewgt.second = sqr(Q)*realGammaGammaqME(selectedParticles,selectedMomenta,FIQED1,Shower,false); wgt *= mewgt.first+mewgt.second; } else { mewgt.first = sqr(Q)*realGammaGammaqbarME(selectedParticles,selectedMomenta,IFQED2,Shower,false); mewgt.second = sqr(Q)*realGammaGammaqbarME(selectedParticles,selectedMomenta,FIQED2,Shower,false); wgt *= mewgt.first+mewgt.second; } if(wgt>1.) generator()->log() << "Weight greater than one in " << "MEPP2GammaGammaPowheg::hardQEDEmission() " << "for IF channel " << " Weight = " << wgt << "\n"; } while(xT>xTMin&&UseRandom::rnd()>wgt); // if no emission if(xTpT()[ShowerInteraction::QED] = minpT_; return born; } pTmax = 0.5*xT*Q; iemit = mewgt.first>mewgt.second ? 2 : 3; // construct the object needed to perform the showers // create the partons ParticleVector newparticles; newparticles.push_back(selectedParticles[0]->produceParticle(selectedMomenta[0])); newparticles.push_back(selectedParticles[1]->produceParticle(selectedMomenta[1])); for(unsigned int ix=2;ix dataPtr()->produceParticle(selectedMomenta[ix==2 ? 2 : 4 ])); } newparticles.push_back(selectedParticles[3]->produceParticle(selectedMomenta[3])); // make the colour connections bool col = newparticles[3]->id()<0; if(particlesToShower[0]->dataPtr()->charged()) { newparticles[3]->incomingColour(newparticles[1],col); newparticles[1]->colourConnect (newparticles[0],col); } else { newparticles[3]->incomingColour(newparticles[0],col); newparticles[0]->colourConnect (newparticles[1],col); } bool FSR = iemit==3; // add incoming particles if(particlesToShower[0]==born->bornIncoming()[0]) { born->incoming().push_back(newparticles[0]); born->incoming().push_back(newparticles[1]); } else { born->incoming().push_back(newparticles[1]); born->incoming().push_back(newparticles[0]); } // IS radiatng particle unsigned int iemitter = born->incoming()[0]->dataPtr()->charged() ? 0 : 1; // add outgoing particles if(particlesToShower[2]==born->bornOutgoing()[0]) { born->outgoing().push_back(newparticles[2]); born->outgoing().push_back(newparticles[3]); } else { born->outgoing().push_back(newparticles[3]); born->outgoing().push_back(newparticles[2]); } born->outgoing().push_back(newparticles[4]); // outgoing radiating particle unsigned int ispectator = born->outgoing()[0]->dataPtr()->charged() ? 2 : 3; // get emitter and spectator right if(FSR) swap(iemitter,ispectator); born->emitter (iemitter ); born->spectator(ispectator); born->emitted(4); // x values pair xnew; for(unsigned int ix=0;ix<2;++ix) { double x = born->incoming()[ix]->momentum().rho()/born->hadrons()[ix]->momentum().rho(); if(ix==0) xnew.first = x; else xnew.second = x; } born->x(xnew); // max pT born->pT()[ShowerInteraction::QED] = pTmax; // return the process born->interaction(ShowerInteraction::QED); return born; } diff --git a/Sampling/exsample/cell.icc b/Sampling/exsample/cell.icc --- a/Sampling/exsample/cell.icc +++ b/Sampling/exsample/cell.icc @@ -1,490 +1,490 @@ // -*- C++ -*- // // cell.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions // // Copyright (C) 2008-2017 Simon Platzer -- simon.plaetzer@desy.de, The Herwig Collaboration // // ExSample is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // namespace exsample { template void cell_info::select (Random& rnd_gen, std::vector& p) { std::transform(lower_left_.begin(),lower_left_.end(), upper_right_.begin(),p.begin(), rnd_gen); ++attempted_; } template void cell_info::select (Random& rnd_gen, std::vector& p, const std::vector& sample) { conditional_transform(lower_left_.begin(),lower_left_.end(), upper_right_.begin(),sample.begin(), p.begin(),rnd_gen); ++attempted_; } template void cell_info::explore(Random& rnd_gen, const adaption_info& ainfo, Function* function, statistics* stats, SlaveStatistics& opt, double detuning) { function->start_presampling(); unsigned long n_sampled = 0; std::vector ll = lower_left_; std::vector ur = upper_right_; double val = 0.; std::vector pos (ll.size()); std::vector< std::pair > > vals; unsigned long ivalnonzero = 0; while (n_sampled < ainfo.presampling_points) { std::transform(ll.begin(),ll.end(), ur.begin(),pos.begin(), rnd_gen); val = function->evaluate(pos) * detuning; vals.push_back( std::pair > (val,pos) ); if ( val != 0 ) ivalnonzero++; ++n_sampled; } while ( ivalnonzero > 0 ) { double avg = 0; double err = 0; double maxval = 0; std::vector maxpos (ll.size()); - unsigned long imax; + unsigned long imax(0); for ( unsigned long ival=0; ival < vals.size(); ival++ ) { val = std::abs(vals[ival].first); if ( val == 0 ) continue; avg += val; err += sqr(val); if ( val > maxval ) { maxval = val; maxpos = vals[ival].second; imax = ival; } } avg /= ivalnonzero; err /= ivalnonzero; err = sqrt(err-sqr(avg)); if ( maxval <= avg+sqrt(ivalnonzero/2.)*err ) { overestimate_ = maxval; last_max_position_ = maxpos; break; } vals.erase(vals.begin()+imax); ivalnonzero--; } for ( unsigned long ival=0; ival < vals.size(); ival++ ) { val = vals[ival].first; stats->presampled(val); opt.select(val); selected(pos,std::abs(val),ainfo); } function->stop_presampling(); } template void cell_info::explore (Random& rnd_gen, const adaption_info& ainfo, Function* function, double detuning) { function->start_presampling(); unsigned long n_sampled = 0; std::vector ll = lower_left_; std::vector ur = upper_right_; double val = 0.; std::vector pos (ll.size()); std::vector< std::pair > > vals; while (n_sampled < ainfo.presampling_points) { std::transform(ll.begin(),ll.end(), ur.begin(),pos.begin(), rnd_gen); val = function->evaluate(pos) * detuning; if ( std::abs(val) > 0 ) vals.push_back( std::pair > (std::abs(val),pos) ); ++n_sampled; } while ( vals.size() > 0 ) { double avg = 0; double err = 0; double maxval = 0; std::vector maxpos (ll.size()); - unsigned long imax; + unsigned long imax(0); for ( unsigned long ival=0; ival < vals.size(); ival++ ) { double thisval = vals[ival].first; avg += thisval; err += sqr(thisval); if ( thisval > maxval ) { maxval = thisval; maxpos = vals[ival].second; imax = ival; } } avg /= vals.size(); err /= vals.size(); err = sqrt(err-sqr(avg)); if ( maxval <= avg+sqrt(vals.size()/2.)*err ) { overestimate_ = maxval; last_max_position_ = maxpos; break; } vals.erase(vals.begin()+imax); } function->stop_presampling(); } template void cell_info::put (OStream& os) const { os << overestimate_; ostream_traits::separator(os); os << volume_; ostream_traits::separator(os); os << lower_left_.size(); ostream_traits::separator(os); for (std::size_t k = 0; k < lower_left_.size(); ++k) { os << lower_left_[k]; ostream_traits::separator(os); } for (std::size_t k = 0; k < upper_right_.size(); ++k) { os << upper_right_[k]; ostream_traits::separator(os); } for (std::size_t k = 0; k < mid_point_.size(); ++k) { os << mid_point_[k]; ostream_traits::separator(os); } for (std::size_t k = 0; k < last_max_position_.size(); ++k) { os << last_max_position_[k]; ostream_traits::separator(os); } for (std::size_t k = 0; k < avg_weight_.size(); ++k) { os << avg_weight_[k].first; ostream_traits::separator(os); os << avg_weight_[k].second; ostream_traits::separator(os); } os << attempted_; ostream_traits::separator(os); os << accepted_; ostream_traits::separator(os); os << parametric_missing_map_.size(); ostream_traits::separator(os); for ( std::map,int>::const_iterator p = parametric_missing_map_.begin(); p != parametric_missing_map_.end(); ++p ) { p->first.put(os); os << p->second; ostream_traits::separator(os); } } template void cell_info::get (IStream& is) { std::size_t dim; is >> overestimate_ >> volume_ >> dim; lower_left_.resize(dim); for (std::size_t k = 0; k < lower_left_.size(); ++k) { is >> lower_left_[k]; } upper_right_.resize(dim); for (std::size_t k = 0; k < upper_right_.size(); ++k) { is >> upper_right_[k]; } mid_point_.resize(dim); for (std::size_t k = 0; k < mid_point_.size(); ++k) { is >> mid_point_[k]; } last_max_position_.resize(dim); for (std::size_t k = 0; k < last_max_position_.size(); ++k) { is >> last_max_position_[k]; } avg_weight_.resize(dim); for (std::size_t k = 0; k < avg_weight_.size(); ++k) { is >> avg_weight_[k].first >> avg_weight_[k].second; } is >> attempted_ >> accepted_ >> dim; for ( size_t k = 0; k < dim; ++k ) { bit_container in; in.get(is); is >> parametric_missing_map_[in]; } } template std::pair cell::split (std::pair split_d, Random& rnd_gen, Function* function, const adaption_info& ainfo, const std::vector& sampled, double detuning) { assert(!missing_events() && !info().parametric_compensating()); split_dimension_ = split_d.first; split_point_ = split_d.second; std::vector lower_left1 = info().lower_left(); std::vector upper_right1 = info().upper_right(); std::vector lower_left2 = info().lower_left(); std::vector upper_right2 = info().upper_right(); upper_right1[split_dimension_] = split_point_; lower_left2[split_dimension_] = split_point_; std::pair children; if (sampled.empty()) children = std::pair(cell(lower_left1,upper_right1,ainfo), cell(lower_left2,upper_right2,ainfo)); else children = std::pair (cell(lower_left1,upper_right1,sampled,ainfo), cell(lower_left2,upper_right2,sampled,ainfo)); if (info().last_max_position()[split_dimension_] <= split_point_) { children.first.info().overestimate(info().overestimate(),info().last_max_position(),1.0); children.second.info().explore(rnd_gen,ainfo,function,detuning); } else { children.second.info().overestimate(info().overestimate(),info().last_max_position(),1.0); children.first.info().explore(rnd_gen,ainfo,function,detuning); } cell_info_.reset(0); children.first.integral(children.first.info().volume() * children.first.info().overestimate()); children.second.integral(children.second.info().volume() * children.second.info().overestimate()); return children; } template void cell::put (OStream& os) const { os << split_dimension_; ostream_traits::separator(os); os << split_point_; ostream_traits::separator(os); os << integral_; ostream_traits::separator(os); os << missing_events_; ostream_traits::separator(os); if (cell_info_) { os << "has_cell_info"; ostream_traits::separator(os); cell_info_->put(os); } else { os << "has_no_cell_info"; ostream_traits::separator(os); } } template void cell::get (IStream& is) { std::string info_tag; is >> split_dimension_ >> split_point_ >> integral_ >> missing_events_ >> info_tag; if (info_tag == "has_cell_info") { cell_info_.reset(new cell_info()); cell_info_->get(is); } } inline cell_info::cell_info() : overestimate_(0.), volume_(0.), lower_left_(), upper_right_(), mid_point_(), last_max_position_(), avg_weight_(), attempted_(0), accepted_(0) {} inline cell_info::cell_info(const std::vector& ll, const std::vector& ur, const adaption_info& ainfo) : overestimate_(0.), volume_(), lower_left_(ll), upper_right_(ur), mid_point_(), last_max_position_(), avg_weight_(std::vector > (ainfo.dimension,std::make_pair(0.,0.))), attempted_(0), accepted_(0) { std::vector delta; std::transform(ur.begin(),ur.end(), ll.begin(),std::back_inserter(delta), std::minus()); volume_ = std::accumulate(delta.begin(),delta.end(),1.,std::multiplies()); std::transform(ur.begin(),ur.end(), ll.begin(),std::back_inserter(mid_point_), std::plus()); for (std::size_t k = 0; k < ainfo.dimension; ++k) mid_point_[k] /= 2.; } inline cell_info::cell_info(const std::vector& ll, const std::vector& ur, const std::vector& sampled_variables, const adaption_info& ainfo) : overestimate_(0.), volume_(), lower_left_(ll), upper_right_(ur), mid_point_(), last_max_position_(), avg_weight_(std::vector > (ainfo.dimension,std::make_pair(0.,0.))), attempted_(0), accepted_(0) { std::vector delta; conditional_transform(ur.begin(),ur.end(), ll.begin(),sampled_variables.begin(), std::back_inserter(delta), std::minus()); volume_ = std::accumulate(delta.begin(),delta.end(),1.,std::multiplies()); std::transform(ur.begin(),ur.end(), ll.begin(),std::back_inserter(mid_point_), std::plus()); for (std::size_t k = 0; k < ainfo.dimension; ++k) mid_point_[k] /= 2.; } inline int cell_info::parametric_missing(const bit_container& id) const { std::map,int>::const_iterator mit = parametric_missing_map_.find(id); if (mit == parametric_missing_map_.end()) return 0; return mit->second; } inline void cell_info::parametric_missing(const bit_container& id, int n) { if (n == 0) { std::map,int>::iterator mit = parametric_missing_map_.find(id); if (mit != parametric_missing_map_.end()) parametric_missing_map_.erase(mit); return; } parametric_missing_map_[id] = n; } inline void cell_info::increase_parametric_missing(const bit_container& id) { std::map,int>::iterator mit = parametric_missing_map_.find(id); if (mit != parametric_missing_map_.end()) { mit->second += 1; if (mit->second == 0) parametric_missing_map_.erase(mit); } else parametric_missing_map_[id] = 1; } inline void cell_info::decrease_parametric_missing(const bit_container& id) { std::map,int>::iterator mit = parametric_missing_map_.find(id); if (mit != parametric_missing_map_.end()) { mit->second -= 1; if (mit->second == 0) parametric_missing_map_.erase(mit); } else assert(false); } inline void cell_info::selected(const std::vector& p, double weight, const adaption_info& ainfo) { for (std::size_t k = 0; k < p.size(); ++k) { if (ainfo.adapt[k]) { if (p[k] < mid_point_[k]) avg_weight_[k].first += weight; else avg_weight_[k].second += weight; } } } inline std::pair cell_info::get_split (const adaption_info& ainfo, bool& worth) const { std::size_t split_d = 0; double gain = 0.; for (std::size_t k = 0; k < ainfo.dimension; ++k) { double xgain = 0.; double left = avg_weight_[k].first; double right = avg_weight_[k].second; if (left+right > 0.) { xgain = std::abs(left-right)/(left+right); } if (xgain > gain) { gain = xgain; split_d = k; } } worth = (gain >= ainfo.gain_threshold); return std::make_pair(split_d,mid_point_[split_d]); } inline bool cell_info::contains_parameter (const std::vector& point, const std::vector& sampled) const { std::vector::const_iterator p = point.begin(); std::vector::const_iterator l = lower_left_.begin(); std::vector::const_iterator u = upper_right_.begin(); std::vector::const_iterator f = sampled.begin(); for (; p < point.end(); ++p, ++f, ++l, ++u) if (!(*f)) { if (((*l) > (*p)) || ((*u) < (*p))) return false; } return true; } inline cell::cell() : split_dimension_(0), split_point_(0.), integral_(0.), missing_events_(0), cell_info_(nullptr) {} inline cell::cell(const std::vector& ll, const std::vector& ur, const adaption_info& ainfo) : split_dimension_(0), split_point_(0.), integral_(0.), missing_events_(0), cell_info_(new cell_info(ll,ur,ainfo)) {} inline cell::cell(const std::vector& ll, const std::vector& ur, const std::vector& sampled_variables, const adaption_info& ainfo) : split_dimension_(0), split_point_(0.), integral_(0.), missing_events_(0), cell_info_(new cell_info(ll,ur,sampled_variables,ainfo)) {} inline cell::cell(const cell& x) : split_dimension_(x.split_dimension_), split_point_(x.split_point_), integral_(x.integral_), missing_events_(x.missing_events_), cell_info_(nullptr) { if (x.cell_info_) cell_info_.reset(new cell_info(*x.cell_info_)); } inline cell& cell::operator=(const cell& x) { if (this == &x) return *this; split_dimension_ = x.split_dimension_; split_point_ = x.split_point_; integral_ = x.integral_; missing_events_ = x.missing_events_; if (x.cell_info_) cell_info_.reset(new cell_info(*x.cell_info_)); return *this; } }