diff --git a/MatrixElement/Lepton/MEee2Mesons2.cc b/MatrixElement/Lepton/MEee2Mesons2.cc --- a/MatrixElement/Lepton/MEee2Mesons2.cc +++ b/MatrixElement/Lepton/MEee2Mesons2.cc @@ -1,352 +1,411 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEee2Mesons2 class. // #include "MEee2Mesons2.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/MatrixElement/Tree2toNDiagram.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "ThePEG/Cuts/Cuts.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "Herwig/Decay/DecayIntegrator.h" #include "ThePEG/PDF/PolarizedBeamParticleData.h" using namespace Herwig; typedef LorentzVector<complex<InvEnergy> > LorentzPolarizationVectorInvE; void MEee2Mesons2::getDiagrams() const { // make sure the current got initialised current_->init(); // max energy Energy Emax = generator()->maximumCMEnergy(); // incoming particles tPDPtr em = getParticleData(ParticleID::eminus); tPDPtr ep = getParticleData(ParticleID::eplus); // loop over the modes int nmode=0; int ndiag=0; for(unsigned int imode=0;imode<current_->numberOfModes();++imode) { // get the external particles for this mode int iq(0),ia(0); tPDVector out = current_->particles(0,imode,iq,ia); current_->decayModeInfo(imode,iq,ia); if(iq==2&&ia==-2) continue; PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(em,out,1.,ep,Emax)); PhaseSpaceChannel channel(mode); if(!current_->createMode(0,tcPDPtr(), IsoSpin::IUnknown, IsoSpin::I3Unknown, imode,mode,0,-1,channel,Emax)) continue; for(unsigned int ix=0;ix<mode->channels().size();++ix) { ThePEG::Ptr<ThePEG::Tree2toNDiagram>::pointer diag = mode->channels()[ix].createDiagram(); ndiag+=1; diag = new_ptr((*diag,-ndiag)); add(diag); } + nmult_= mode->outgoing().size(); } - nmult_=3; } Energy2 MEee2Mesons2::scale() const { return sHat(); } int MEee2Mesons2::nDim() const { return 3*nmult_-5; } void MEee2Mesons2::setKinematics() { HwMEBase::setKinematics(); } bool MEee2Mesons2::generateKinematics(const double * r) { using Constants::pi; // Save the jacobian dPS/dr for later use. jacobian(1.0); // set the masses of the outgoing particles for ( int i = 2, N = meMomenta().size(); i < N; ++i ) { meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->mass()); } double ctmin = -1.0, ctmax = 1.0; double cth = getCosTheta(ctmin, ctmax, r[0]); phi(rnd(2.0*Constants::pi)); unsigned int i1(2),i2(3),i3(4); Energy q = ZERO; Energy e = sqrt(sHat())/2.0; Energy2 pq; if(nDim()==1) { try { q = SimplePhaseSpace:: getMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass()); } catch ( ImpossibleKinematics ) { return false; } pq = 2.0*e*q; Energy pt = q*sqrt(1.0-sqr(cth)); 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(); Energy2 m22 = meMomenta()[2].mass2(); Energy2 m32 = meMomenta()[3].mass2(); Energy2 e0e2 = 2.0*e*sqrt(sqr(q) + m22); tHat(pq*cth + m22 - e0e2); uHat(m22 + m32 - sHat() - tHat()); } else if(nDim()==4) { double rm=r[1]; if(rm<1./3.) { rm *=3.; } else if(rm<2./3.) { rm = 3.*rm-1.; swap(i2,i3); } else { rm = 3.*rm-2.; swap(i1,i2); } tcPDPtr res = getParticleData(213); Energy mass = res->mass(), width = res->width(); Energy2 m2max = sqr(sqrt(sHat())-meMomenta()[i3].mass()); Energy2 m2min = sqr(meMomenta()[i1].mass()+meMomenta()[i2].mass()); double rhomin = atan((m2min-sqr(mass))/mass/width); double rhomax = atan((m2max-sqr(mass))/mass/width); double rho = rhomin+rm*(rhomax-rhomin); Energy2 m2 = mass*(width*tan(rho)+mass); Energy mv = sqrt(m2); try { q = SimplePhaseSpace:: getMagnitude(sHat(), mv, meMomenta()[i3].mass()); } catch ( ImpossibleKinematics ) { return false; } pq = 2.0*e*q; Energy pt = q*sqrt(1.0-sqr(cth)); Lorentz5Momentum poff; poff.setMass(mv); poff.setVect(Momentum3( pt*sin(phi()), pt*cos(phi()), q*cth)); meMomenta()[i3].setVect(Momentum3(-pt*sin(phi()), -pt*cos(phi()), -q*cth)); poff.rescaleEnergy(); meMomenta()[i3].rescaleEnergy(); // decay of the intermediate bool test=Kinematics::twoBodyDecay(poff,meMomenta()[i1].mass(), meMomenta()[i2].mass(), -1.+2*r[2],r[3]*2.*pi, meMomenta()[i1],meMomenta()[i2]); if(!test) return false; // decay piece of the jacobian Energy p2 = Kinematics::pstarTwoBodyDecay(mv,meMomenta()[i1].mass(), meMomenta()[i2].mass()); jacobian(p2/mv/8./sqr(pi)*jacobian()); // mass piece jacobian((rhomax-rhomin)*( sqr(m2-sqr(mass))+sqr(mass*width)) /mass/width*jacobian()/sHat()); } + else if(nDim()==7) { + tcPDPtr res = getParticleData(213); + Energy mass = res->mass(), width = res->width(); + Energy m12min = meMomenta()[2].mass()+meMomenta()[3].mass(); + Energy m34min = meMomenta()[4].mass()+meMomenta()[5].mass(); + + Energy m12max = sqrt(sHat())-m34min; + double rhomin = atan((sqr(m12min)-sqr(mass))/mass/width); + double rhomax = atan((sqr(m12max)-sqr(mass))/mass/width); + double rho = rhomin+r[1]*(rhomax-rhomin); + jacobian(jacobian()*(rhomax-rhomin)); + Energy2 m122 = mass*(width*tan(rho)+mass); + Energy m12 = sqrt(m122); + + Energy m34max = sqrt(sHat())-m12; + rhomin = atan((sqr(m34min)-sqr(mass))/mass/width); + rhomax = atan((sqr(m34max)-sqr(mass))/mass/width); + rho = rhomin+r[2]*(rhomax-rhomin); + jacobian(jacobian()*(rhomax-rhomin)); + Energy2 m342 = mass*(width*tan(rho)+mass); + Energy m34 = sqrt(m342); + try { + q = SimplePhaseSpace:: + getMagnitude(sHat(), m12, m34); + } + catch ( ImpossibleKinematics ) { + return false; + } + pq = 2.0*e*q; + Energy pt = q*sqrt(1.0-sqr(cth)); + Lorentz5Momentum p12,p34; + p12.setMass(m12); + p34.setMass(m34); + p12.setVect(Momentum3( pt*sin(phi()), pt*cos(phi()), q*cth)); + p34.setVect(Momentum3(-pt*sin(phi()), -pt*cos(phi()), -q*cth)); + p12.rescaleEnergy(); + p34.rescaleEnergy(); + // decay of the intermediate + bool test=Kinematics::twoBodyDecay(p12,meMomenta()[2].mass(), + meMomenta()[3].mass(), + -1.+2*r[3],r[4]*2.*pi, + meMomenta()[2],meMomenta()[3]); + if(!test) return false; + test=Kinematics::twoBodyDecay(p34,meMomenta()[4].mass(), + meMomenta()[5].mass(), + -1.+2*r[5],r[6]*2.*pi, + meMomenta()[4],meMomenta()[5]); + if(!test) return false; + // decay piece of the jacobian + Energy p2 = Kinematics::pstarTwoBodyDecay(m12,meMomenta()[2].mass(), + meMomenta()[3].mass()); + jacobian(p2/m12/8./sqr(pi)*jacobian()); + p2 = Kinematics::pstarTwoBodyDecay(m34,meMomenta()[4].mass(), + meMomenta()[5].mass()); + jacobian(p2/m34/8./sqr(pi)*jacobian()); + // mass piece + jacobian((sqr(m122-sqr(mass))+sqr(mass*width))/mass/width/sHat()*jacobian()); + jacobian((sqr(m342-sqr(mass))+sqr(mass*width))/mass/width/sHat()*jacobian()); + } else { assert(false); } vector<LorentzMomentum> out(meMomenta().size()-2); tcPDVector tout(meMomenta().size()); for(unsigned int ix=2;ix<meMomenta().size();++ix) { out[ix-2] = meMomenta()[ix]; tout[ix-2] = mePartonData()[ix]; } if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) ) return false; jacobian((pq/sHat())*Constants::pi*jacobian()); return true; } double MEee2Mesons2::me2() const { // cerr << "testing in me2\n"; // for(unsigned int ix=0;ix<mePartonData().size();++ix) // cerr << mePartonData()[ix]->PDGName() << " " << meMomenta()[ix]/GeV << "\n"; // wavefunctions for the incoming leptons SpinorWaveFunction em_in( meMomenta()[0],mePartonData()[0],incoming); SpinorBarWaveFunction ep_in( meMomenta()[1],mePartonData()[1],incoming); vector<SpinorWaveFunction> f1; vector<SpinorBarWaveFunction> a1; 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); } // compute the leptonic current LorentzPolarizationVectorInvE lepton[2][2]; InvEnergy2 pre = SM().alphaEM(sHat())*4.*Constants::pi/sHat(); for(unsigned ix=0;ix<2;++ix) { for(unsigned iy=0;iy<2;++iy) { lepton[ix][iy]= pre*f1[ix].dimensionedWave().vectorCurrent(a1[iy].dimensionedWave()); } } // work out the mapping for the hadron vector vector<unsigned int> constants(meMomenta().size()+1); vector<PDT::Spin > ispin(meMomenta().size()-2); vector<int> hadrons(meMomenta().size()-2); int itemp(1); unsigned int ix(meMomenta().size()); do { --ix; ispin[ix-2] = mePartonData()[ix]->iSpin(); hadrons[ix-2] = mePartonData()[ix]->id(); itemp *= ispin[ix-2]; constants[ix] = itemp; } while(ix>2); constants[meMomenta().size()] = 1; constants[0] = constants[1] = constants[2]; // cerr << "testing constants\n"; // for(unsigned int ix=0;ix<constants.size();++ix) // cerr << constants[ix] << " "; // cerr << "\n"; // cerr << "testing ispn\n"; // for(unsigned int ix=0;ix<ispin.size();++ix) // cerr << ispin[ix] << " "; // cerr << "\n"; // calculate the matrix element me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,ispin)); // calculate the hadron current unsigned int imode = current_->decayMode(hadrons); Energy q = sqrt(sHat()); vector<Lorentz5Momentum> momenta; tPDVector out; for(unsigned int ix=2;ix<meMomenta().size();++ix) { momenta.push_back(meMomenta()[ix]); out.push_back(const_ptr_cast<tPDPtr>(mePartonData()[ix])); } vector<LorentzPolarizationVectorE> hadron(current_->current(tcPDPtr(), IsoSpin::IUnknown, IsoSpin::I3Unknown, imode,-1, q,out,momenta,DecayIntegrator2::Calculate)); // for(unsigned int ix=0;ix<hadron.size();++ix) // cerr << hadron[ix].x()/GeV << " " << hadron[ix].y()/GeV << " " << hadron[ix].z()/GeV << " " << hadron[ix].t()/GeV << "\n"; // compute the matrix element vector<unsigned int> ihel(meMomenta().size()); double output(0.); for(unsigned int hhel=0;hhel<hadron.size();++hhel) { // map the index for the hadrons to a helicity state for(unsigned int ix=meMomenta().size()-1;ix>1;--ix) { ihel[ix]=(hhel%constants[ix-1])/constants[ix]; } // loop over the helicities of the incoming leptons for(ihel[1]=0;ihel[1]<2;++ihel[1]){ for(ihel[0]=0;ihel[0]<2;++ihel[0]) { Complex amp = lepton[ihel[0]][ihel[1]].dot(hadron[hhel]); me_(ihel)= amp; output += std::norm(amp); // cerr << "testing ME "; // for(unsigned int ix=0;ix<ihel.size();++ix) cerr << ihel[ix] << " "; // cerr << me_(ihel) << "\n"; } } } output *= 0.25*sqr(pow(sqrt(sHat())/q,int(momenta.size()-2))); // polarization stuff tcPolarizedBeamPDPtr beam[2] = {dynamic_ptr_cast<tcPolarizedBeamPDPtr>(mePartonData()[0]), dynamic_ptr_cast<tcPolarizedBeamPDPtr>(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())}; output = me_.average(rho[0],rho[1]); } return output; } CrossSection MEee2Mesons2::dSigHatDR() const { return me2()*jacobian()/(16.0*sqr(Constants::pi)*sHat())*sqr(hbarc); } unsigned int MEee2Mesons2::orderInAlphaS() const { return 0; } unsigned int MEee2Mesons2::orderInAlphaEW() const { return 0; } Selector<MEBase::DiagramIndex> MEee2Mesons2::diagrams(const DiagramVector & diags) const { // This example corresponds to the diagrams specified in the example // in the getDiagrams() function. Selector<DiagramIndex> sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ) if ( diags[i]->id() == -1 ) sel.insert(1.0, i); else if ( diags[i]->id() == -2 ) sel.insert(1.0, i); else if ( diags[i]->id() == -3 ) sel.insert(1.0, i); // You probably do not want equal weights here... return sel; // If there is only one possible diagram you can override the // MEBase::diagram function instead. } Selector<const ColourLines *> MEee2Mesons2::colourGeometries(tcDiagPtr) const { static ColourLines none(""); Selector<const ColourLines *> sel; sel.insert(1.0, &none); return sel; } IBPtr MEee2Mesons2::clone() const { return new_ptr(*this); } IBPtr MEee2Mesons2::fullclone() const { return new_ptr(*this); } void MEee2Mesons2::doinit() { HwMEBase::doinit(); } void MEee2Mesons2::persistentOutput(PersistentOStream & os) const { os << current_ << nmult_; } void MEee2Mesons2::persistentInput(PersistentIStream & is, int) { is >> current_ >> nmult_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass<MEee2Mesons2,HwMEBase> describeHerwigMEee2Mesons2("Herwig::MEee2Mesons2", "HwMELeptonLowEnergy.so"); void MEee2Mesons2::Init() { static ClassDocumentation<MEee2Mesons2> documentation ("The MEee2Mesons2 class simluation the production of low multiplicity" " events via the weak current"); static Reference<MEee2Mesons2,WeakCurrent> interfaceWeakCurrent ("WeakCurrent", "The reference for the decay current to be used.", &MEee2Mesons2::current_, false, false, true, false, false); } void MEee2Mesons2::constructVertex(tSubProPtr sub) { }