Page MenuHomeHEPForge
Files F8309510

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/EW/ElectroWeakReweighter.cc b/MatrixElement/EW/ElectroWeakReweighter.cc
--- a/MatrixElement/EW/ElectroWeakReweighter.cc
+++ b/MatrixElement/EW/ElectroWeakReweighter.cc
@@ -1,200 +1,659 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ElectroWeakReweighter class.
//
#include "ElectroWeakReweighter.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 "boost/numeric/ublas/matrix.hpp"
#include "boost/numeric/ublas/operation.hpp"
#include "EWProcess.h"
#include "HighEnergyMatching.h"
#include "ElectroWeakMatching.h"
+#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
+#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
+#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
+#include "ThePEG/Helicity/epsilon.h"
using namespace Herwig;
tEWCouplingsPtr ElectroWeakReweighter::staticEWCouplings_ = tEWCouplingsPtr();
ElectroWeakReweighter::ElectroWeakReweighter() {}
ElectroWeakReweighter::~ElectroWeakReweighter() {}
IBPtr ElectroWeakReweighter::clone() const {
return new_ptr(*this);
}
IBPtr ElectroWeakReweighter::fullclone() const {
return new_ptr(*this);
}
void ElectroWeakReweighter::persistentOutput(PersistentOStream & os) const {
os << EWCouplings_ << collinearSudakov_ << softSudakov_;
}
void ElectroWeakReweighter::persistentInput(PersistentIStream & is, int) {
is >> EWCouplings_ >> collinearSudakov_ >> softSudakov_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<ElectroWeakReweighter,ReweightBase>
describeHerwigElectroWeakReweighter("Herwig::ElectroWeakReweighter", "HwMEEW.so");
void ElectroWeakReweighter::Init() {
static ClassDocumentation<ElectroWeakReweighter> documentation
("There is no documentation for the ElectroWeakReweighter class");
static Reference<ElectroWeakReweighter,EWCouplings> interfaceEWCouplings
("EWCouplings",
"The object to calculate the electroweak couplings",
&ElectroWeakReweighter::EWCouplings_, false, false, true, false, false);
static Reference<ElectroWeakReweighter,CollinearSudakov> interfaceCollinearSudakov
("CollinearSudakov",
"The collinear Sudakov",
&ElectroWeakReweighter::collinearSudakov_, false, false, true, false, false);
static Reference<ElectroWeakReweighter,SoftSudakov> interfaceSoftSudakov
("SoftSudakov",
"The soft Sudakov",
&ElectroWeakReweighter::softSudakov_, false, false, true, false, false);
}
namespace {
void axpy_prod_local(const boost::numeric::ublas::matrix<Complex> & A,
const boost::numeric::ublas::matrix<complex<InvEnergy2> > & B,
boost::numeric::ublas::matrix<complex<InvEnergy2> > & C) {
assert(A.size2()==B.size1());
C.resize(A.size1(),B.size2());
for(unsigned int ix=0;ix<A.size1();++ix) {
for(unsigned int iy=0;iy<B.size2();++iy) {
C(ix,iy) = ZERO;
for(unsigned int iz=0;iz<A.size2();++iz) {
C(ix,iy) += A(ix,iz)*B(iz,iy);
}
}
}
}
+void axpy_prod_local(const boost::numeric::ublas::matrix<complex<InvEnergy2> > & A,
+ const boost::numeric::ublas::vector<complex<Energy2> > & B,
+ boost::numeric::ublas::vector<Complex > & C) {
+ assert(A.size2()==B.size());
+ C.resize(A.size1());
+ for(unsigned int ix=0;ix<A.size1();++ix) {
+ C(ix) = ZERO;
+ for(unsigned int iz=0;iz<A.size2();++iz) {
+ C(ix) += A(ix,iz)*B(iz);
+ }
+ }
+}
+
void axpy_prod_local(const boost::numeric::ublas::matrix<complex<InvEnergy2> > & A,
const boost::numeric::ublas::matrix<Complex> & B,
boost::numeric::ublas::matrix<complex<InvEnergy2> > & C) {
assert(A.size2()==B.size1());
C.resize(A.size1(),B.size2());
for(unsigned int ix=0;ix<A.size1();++ix) {
for(unsigned int iy=0;iy<B.size2();++iy) {
C(ix,iy) = ZERO;
for(unsigned int iz=0;iz<A.size2();++iz) {
C(ix,iy) += A(ix,iz)*B(iz,iy);
}
}
}
}
}
+
double ElectroWeakReweighter::weight() const {
EWCouplings_->initialize();
staticEWCouplings_ = EWCouplings_;
// cerr << "aEM\n";
// for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.1) {
// cerr << scale/GeV << " "
// << EWCouplings_->aEM(scale) << "\n";
// }
// cerr << "aS\n";
// for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
// cerr << scale/GeV << " "
// << EWCouplings_->aS(scale) << "\n";
// }
// cerr << "y_t\n";
// for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
// cerr << scale/GeV << " "
// << EWCouplings_->y_t(scale) << "\n";
// }
// cerr << "lambda\n";
// for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
// cerr << scale/GeV << " "
// << EWCouplings_->lambda(scale) << "\n";
// }
// cerr << "vev\n";
// for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
// cerr << scale/GeV << " "
// << EWCouplings_->vev(scale)/GeV << "\n";
// }
- collinearSudakov_->makePlots();
- Energy2 s = sqr(5000.*GeV);
- Energy2 t = -0.25*s;
- Energy2 u = -0.75*s;
- testEvolution(s,t,u);
+ // collinearSudakov_->makePlots();
+ // Energy2 s = sqr(5000.*GeV);
+ // Energy2 t = -0.25*s;
+ // Energy2 u = -0.75*s;
+ // testEvolution(s,t,u);
- // cerr << subProcess() << "\n";
- // cerr << *subProcess() << "\n";
- // cerr << subProcess()->outgoing()[0] << *subProcess()->outgoing()[0] << "\n";
- // cerr << subProcess()->outgoing()[0]->spinInfo() << "\n";
- // cerr << subProcess()->outgoing()[0]->spinInfo()->productionVertex() << "\n";
+ cerr << subProcess() << "\n";
+ cerr << *subProcess() << "\n";
+ cerr << subProcess()->outgoing()[0] << *subProcess()->outgoing()[0] << "\n";
+ cerr << subProcess()->outgoing()[0]->spinInfo() << "\n";
+ cerr << subProcess()->outgoing()[0]->spinInfo()->productionVertex() << "\n";
+ if(subProcess()->outgoing().size()!=2)
+ return 1.;
+ if(subProcess()->incoming().first->id()==ParticleID::g &&
+ subProcess()->incoming().second->id()==ParticleID::g) {
+ if(subProcess()->outgoing()[0]->id()==ParticleID::g &&
+ subProcess()->outgoing()[1]->id()==ParticleID::g)
+ return 1.;
+ else if(abs(subProcess()->outgoing()[0]->id())<=5 &&
+ subProcess()->outgoing()[0]->id()==-subProcess()->outgoing()[1]->id()) {
+ return reweightggqqbar();
+ }
+ else
+ assert(false);
+ }
+ else if(abs(subProcess()->incoming().first->id())<=5 &&
+ subProcess()->incoming().first->id()==-subProcess()->incoming().second->id()) {
+ if(subProcess()->outgoing()[0]->id()==ParticleID::g &&
+ subProcess()->outgoing()[1]->id()==ParticleID::g)
+ return reweightqqbargg();
+ else
+ assert(false);
+ }
+ else
+ assert(false);
assert(false);
staticEWCouplings_ = tEWCouplingsPtr();
}
void ElectroWeakReweighter::testEvolution(Energy2 s,Energy2 t, Energy2 u) const {
Energy highScale = sqrt(s);
Energy ewScale = coupling()->mZ();
Energy lowScale = 50.0*GeV;
for (unsigned int i=0; i<45;++i) {
EWProcess::Process process = (EWProcess::Process)i;
cerr << "process " << process << "\n";
// result for all EW and QCD SCET contributions:
boost::numeric::ublas::matrix<complex<InvEnergy2> > highMatch_val
= HighEnergyMatching::highEnergyMatching(highScale,s,t,u,process,true,true);
boost::numeric::ublas::matrix<Complex> highRunning_val
= softSudakov_->highEnergyRunning(highScale,ewScale,s,t,u,process);
boost::numeric::ublas::matrix<Complex> ewMatch_val =
ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,true);
boost::numeric::ublas::matrix<Complex> lowRunning_val =
softSudakov_->lowEnergyRunning(ewScale,lowScale,s,t,u,process);
boost::numeric::ublas::matrix<Complex> collinearHighRunning_val =
collinearSudakov_->highEnergyRunning(highScale,ewScale,s,process,false);
boost::numeric::ublas::matrix<Complex> collinearEWMatch_val =
collinearSudakov_->electroWeakMatching(ewScale,s,process,true);
boost::numeric::ublas::matrix<Complex> collinearLowRunning_val =
collinearSudakov_->lowEnergyRunning(ewScale,lowScale,s,process);
boost::numeric::ublas::matrix<Complex> lowMatchTemp_val =
boost::numeric::ublas::zero_matrix<Complex>(ewMatch_val.size1(),ewMatch_val.size2());
for (unsigned int ii=0; ii<ewMatch_val.size1(); ++ii) {
for (unsigned int jj=0; jj<ewMatch_val.size2(); ++jj) {
lowMatchTemp_val(ii,jj) = collinearEWMatch_val(ii,jj)*ewMatch_val(ii,jj);
}
}
boost::numeric::ublas::matrix<Complex> temp(highRunning_val.size1(),collinearHighRunning_val.size2());
boost::numeric::ublas::axpy_prod(highRunning_val,collinearHighRunning_val,temp);
boost::numeric::ublas::matrix<Complex> temp2(collinearLowRunning_val.size1(),lowRunning_val.size2());
boost::numeric::ublas::axpy_prod(collinearLowRunning_val,lowRunning_val,temp2);
boost::numeric::ublas::matrix<Complex> temp3(temp2.size1(),lowMatchTemp_val.size2());
boost::numeric::ublas::axpy_prod(temp2,lowMatchTemp_val,temp3);
temp2.resize(temp3.size1(),temp.size2());
boost::numeric::ublas::axpy_prod(temp3,temp,temp2);
boost::numeric::ublas::matrix<complex<InvEnergy2> > result(temp2.size1(),highMatch_val.size2());
axpy_prod_local(temp2,highMatch_val,result);
for(unsigned int ix=0;ix<result.size1();++ix) {
for(unsigned int iy=0;iy<result.size2();++iy) {
cerr << s*result(ix,iy) << " ";
}
cerr << "\n";
}
}
}
+
+namespace {
+
+void SackGluonPolarizations(Lorentz5Momentum &p1,
+ Lorentz5Momentum &p2,
+ Lorentz5Momentum &p3,
+ Lorentz5Momentum &p4,
+ Energy2 s, Energy2 t, Energy2 u, Energy2 m2,
+ vector<LorentzVector<Complex> > & eps3,
+ vector<LorentzVector<Complex> > & eps4,
+ unsigned int iopt) {
+ static const Complex I(0.,1.);
+ // p1 is p-, p2 is p+
+ // p3 is k-, p4 is k+
+ // both final-state
+ if(iopt==0) {
+ // swap t and u due Aneesh's defn
+ Energy3 den1 = sqrt((u*t-sqr(m2))*(s-4.*m2));
+ Energy3 den2 = sqrt(s*(u*t-sqr(m2)));
+ LorentzVector<Complex> eps3Para = (m2+t)/den1*p1 -(m2+u)/den1*p2 +(u-t)/den1*p3;
+ LorentzVector<Complex> eps3Perp = 2./den2*epsilon(p1,p2,p3);
+ LorentzVector<Complex> eps4Para = (m2+t)/den1*p2 -(m2+u)/den1*p1 +(u-t)/den1*p4;
+ LorentzVector<Complex> eps4Perp = 2./den2*epsilon(p1,p2,p4);
+ eps3.push_back(sqrt(0.5)*(eps3Para+I*eps3Perp));
+ eps3.push_back(sqrt(0.5)*(eps3Para-I*eps3Perp));
+ eps4.push_back(sqrt(0.5)*(eps4Para+I*eps4Perp));
+ eps4.push_back(sqrt(0.5)*(eps4Para-I*eps4Perp));
+ if(m2!=ZERO) assert(false);
+ }
+ // both initial-state
+ else if(iopt==1) {
+ if(m2!=ZERO) assert(false);
+ LorentzVector<Complex> eps3Para( 1., 0.,0.,0.);
+ LorentzVector<Complex> eps3Perp( 0.,-1.,0.,0.);
+ LorentzVector<Complex> eps4Para(-1.,0.,0., 0.);
+ LorentzVector<Complex> eps4Perp( 0., 1.,0.,0.);
+ eps3.push_back(sqrt(0.5)*(eps3Para+I*eps3Perp));
+ eps3.push_back(sqrt(0.5)*(eps3Para-I*eps3Perp));
+ eps4.push_back(sqrt(0.5)*(eps4Para+I*eps4Perp));
+ eps4.push_back(sqrt(0.5)*(eps4Para-I*eps4Perp));
+ }
+}
+
+}
+
+double ElectroWeakReweighter::reweightqqbargg() const {
+ // momenta and invariants
+ Lorentz5Momentum p1 = subProcess()->incoming().first ->momentum();
+ tcPDPtr q = subProcess()->incoming().first ->dataPtr();
+ Lorentz5Momentum p2 = subProcess()->incoming().second->momentum();
+ tcPDPtr qbar = subProcess()->incoming().second->dataPtr();
+ if(subProcess()->incoming().first->id()<0) {
+ swap(p1,p2 );
+ swap(q ,qbar);
+ }
+ Lorentz5Momentum p3 = subProcess()->outgoing()[0]->momentum();
+ Lorentz5Momentum p4 = subProcess()->outgoing()[1]->momentum();
+ tcPDPtr g = subProcess()->outgoing()[1]->dataPtr();
+ Energy2 s = (p1+p2).m2();
+ Energy2 t = (p1-p4).m2();
+ Energy2 u = (p1-p3).m2();
+ // // boost to partonci rest frame
+ // Lorentz5Momentum psum=p1+p2;
+ // LorentzRotation boost(-psum.boostVector());
+ // p1 *= boost;
+ // p2 *= boost;
+ // p3 *= boost;
+ // p4 *= boost;
+ // cerr << "testing momenta in reweight A " << p1/GeV << "\n";
+ // cerr << "testing momenta in reweight B " << p2/GeV << "\n";
+ // cerr << "testing momenta in reweight C " << p3/GeV << "\n";
+ // cerr << "testing momenta in reweight D " << p4/GeV << "\n";
+ // LO matrix element coefficents
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ bornQQGGweights,bornRRGGweights;
+ // quark left doublet
+ if(q->id()!=5) {
+ bornQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,true);
+ }
+ else {
+ bornQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,true);
+ }
+ // quark right singlet
+ if(abs(subProcess()->incoming().first->id())%2==0)
+ bornRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,true);
+ else
+ bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true);
+ // EW corrected matrix element coefficients
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ EWQQGGweights,EWRRGGweights;
+ // quark left doublet
+ if(q->id()!=5) {
+ EWQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,false);
+ }
+ else {
+ EWQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,false);
+ }
+ // quakr right singlet
+ if(abs(subProcess()->incoming().first->id())%2==0)
+ EWRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,false);
+ else
+ EWRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,false);
+ // cerr << "testing matrices\n";
+ // for(unsigned int ix=0;ix<bornRRGGweights.size1();++ix) {
+ // for(unsigned int iy=0;iy<bornRRGGweights.size2();++iy) {
+ // cerr << bornRRGGweights(ix,iy)*GeV2 << " ";
+ // }
+ // cerr << "\n";
+ // }
+ // for(unsigned int ix=0;ix<bornQQGGweights.size1();++ix) {
+ // for(unsigned int iy=0;iy<bornQQGGweights.size2();++iy) {
+ // cerr << bornQQGGweights(ix,iy)*GeV2 << " ";
+ // }
+ // cerr << "\n";
+ // }
+
+ SpinorWaveFunction qw(p1,q ,incoming);
+ SpinorBarWaveFunction qbarw(p2,qbar,incoming);
+ // VectorWaveFunction g1w(p3,getParticleData(ParticleID::g),outgoing);
+ // VectorWaveFunction g2w(p4,getParticleData(ParticleID::g),outgoing);
+ vector<LorentzVector<Complex> > eps3,eps4;
+ SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps3,eps4,0);
+ boost::numeric::ublas::matrix<Complex>
+ bornME = boost::numeric::ublas::zero_matrix<Complex>(3,3),
+ EWME = boost::numeric::ublas::zero_matrix<Complex>(3,3);
+ // testME = boost::numeric::ublas::zero_matrix<Complex>(3,3);
+ for(unsigned int iq=0;iq<2;++iq) {
+ if(iq==0) {
+ qw.reset (0);
+ qbarw.reset(1);
+ }
+ else {
+ qw.reset (1);
+ qbarw.reset(0);
+ }
+ LorentzVector<complex<Energy> > current = iq==0 ?
+ qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) :
+ qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave());
+ for(unsigned int i1=0;i1<2;++i1) {
+ complex<Energy> d31 = eps3[i1].dot(p1);
+ for(unsigned int i2=0;i2<2;++i2) {
+ // g1w.reset(2*i1);
+ // g2w.reset(2*i2);
+ boost::numeric::ublas::vector<complex<Energy2> > M(5);
+ Complex d34 = eps3[i1].dot(eps4[i2]);
+ complex<Energy> d42 = eps4[i2].dot(p2);
+ // M0 in paper
+ M(0) = qw.dimensionedWave().slash(eps3[i1])
+ .slash(p4-p2).vectorCurrent(qbarw.dimensionedWave()).dot(eps4[i2]);
+
+ // // really t channel
+ // Complex MEU = -Complex(0.,1.)*qw.dimensionedWave().slash(g1w.wave())
+ // .slash(p4-p2).vectorCurrent(qbarw.dimensionedWave()).dot(g2w.wave())/u;
+ // // really u channel
+ // Complex MET = -Complex(0.,1.)*qw.dimensionedWave().slash(g2w.wave())
+ // .slash(p3-p2).vectorCurrent(qbarw.dimensionedWave()).dot(g1w.wave())/t;
+ // // s channel
+ // Complex MES = -Complex(0.,1.)*(current.dot(p3-p4)*g1w.wave().dot(g2w.wave())
+ // -2.*current.dot(g1w.wave())*(g2w.wave().dot(p3))
+ // -2.*current.dot(g2w.wave())*(g1w.wave().dot(p4)))/s;
+
+ // // really t channel
+ // Complex MEU = -Complex(0.,1.)*qw.dimensionedWave().slash(eps3[i1])
+ // .slash(p4-p2).vectorCurrent(qbarw.dimensionedWave()).dot(eps4[i2])/u;
+ // // really u channel
+ // Complex MET = -Complex(0.,1.)*qw.dimensionedWave().slash(eps4[i2])
+ // .slash(p3-p2).vectorCurrent(qbarw.dimensionedWave()).dot(eps3[i1])/t;
+ // // s channel
+ // Complex MES = -Complex(0.,1.)*(current.dot(p3-p4)*eps3[i1].dot(eps4[i2])
+ // -2.*current.dot(eps3[i1])*(eps4[i2].dot(p3))
+ // -2.*current.dot(eps4[i2])*(eps3[i1].dot(p4)))/s;
+ // cerr << "NEW flows " << MEU+MES << " " << MET-MES << "\n";
+
+
+
+
+ // cerr << "testing new U " << MET << "\n";
+ // M4 in paper
+ M(2) = current.dot(eps4[i2])*d31;
+ // M5 in paper
+ M(3) = -current.dot(eps3[i1])*d42;
+ // M1 in paper (missing factor)
+ M(1) = current.dot(p4);
+ // M6 in paper
+ M(4) = M(1)*d31*d42/GeV2;
+ // M1 final factor
+ M(1) *= d34;
+ // coefficient of different contributions
+ boost::numeric::ublas::vector<Complex> Cborn(3),CEW(3),Ctest(3);
+
+
+ // Ctest(0) = 1./6.*( MEU+MET);
+ // Ctest(1) = 0.5*( MEU+MET);
+ // Ctest(2) = 0.5*(MEU+MES-MET+MES);
+ if(iq==0) {
+ axpy_prod_local(bornQQGGweights,M,Cborn);
+ axpy_prod_local(EWQQGGweights ,M,CEW );
+ }
+ else {
+ axpy_prod_local(bornRRGGweights,M,Cborn);
+ axpy_prod_local(EWRRGGweights ,M,CEW );
+ }
+ for(unsigned int ix=0;ix<3;++ix) {
+ for(unsigned int iy=0;iy<3;++iy) {
+ bornME(ix,iy) += Cborn(ix)*conj(Cborn(iy));
+ EWME (ix,iy) += CEW (ix)*conj(CEW (iy));
+ // testME (ix,iy) += Ctest (ix)*conj(Ctest (iy));
+ }
+ }
+ }
+ }
+ }
+ double born = 24.*real(bornME(0,0))+20./3.*real(bornME(1,1))+12.*real(bornME(2,2));
+ double EW = 24.*real(EWME(0,0))+20./3.*real(EWME(1,1))+12.*real(EWME(2,2));
+ // double test = 24.*real(testME(0,0))+20./3.*real(testME(1,1))+12.*real(testME(2,2));
+
+ // double gs2 = 4.*Constants::pi*ElectroWeakReweighter::coupling()->a3(sqrt(s));
+
+ // cerr << "testing born A " << 0.125*born/sqr(gs2)/9. << "\n";
+ // cerr << "testing born B " << 0.125*test/9. << "\n";
+
+ return EW/born;
+}
+
+boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ElectroWeakReweighter::evaluateRunning(EWProcess::Process process, Energy2 s,
+ Energy2 t, Energy2 u, bool born) const {
+ using namespace boost::numeric::ublas;
+ bool SU3save = coupling()->SU3();
+ bool EWsave = coupling()-> EW();
+ Energy highScale = sqrt(s);
+ Energy ewScale = coupling()->mZ();
+ Energy lowScale = ewScale;
+ // result for all EW and QCD SCET contributions:
+ // high energy matching
+ matrix<complex<InvEnergy2> > highMatch_val
+ = HighEnergyMatching::highEnergyMatching(highScale,s,t,u,process,!born,false);
+ // low energy matching
+ matrix<Complex>
+ ewMatch_val = ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,!born);
+ matrix<Complex> collinearEWMatch_val =
+ collinearSudakov_->electroWeakMatching(ewScale,s,process,!born);
+ matrix<Complex> highRunning_val,lowRunning_val,
+ collinearHighRunning_val,collinearLowRunning_val;
+ if(born) {
+ highRunning_val = identity_matrix<Complex>(softSudakov_->numberGauge(process));
+ lowRunning_val = identity_matrix<Complex>(softSudakov_->numberBrokenGauge(process));
+ collinearHighRunning_val = identity_matrix<Complex>(softSudakov_->numberGauge(process));
+ collinearLowRunning_val = identity_matrix<Complex>(softSudakov_->numberBrokenGauge(process));
+ }
+ else {
+ coupling()->SU3(false);
+ coupling()-> EW( true);
+ highRunning_val = softSudakov_->highEnergyRunning(highScale,ewScale,s,t,u,process);
+ lowRunning_val = softSudakov_->lowEnergyRunning(ewScale,lowScale,s,t,u,process);
+ collinearHighRunning_val = collinearSudakov_->highEnergyRunning(highScale,ewScale,s,process,false);
+ collinearLowRunning_val = collinearSudakov_->lowEnergyRunning(ewScale,lowScale,s,process);
+ };
+ matrix<Complex> lowMatchTemp_val =
+ zero_matrix<Complex>(ewMatch_val.size1(),ewMatch_val.size2());
+ for (unsigned int ii=0; ii<ewMatch_val.size1(); ++ii) {
+ for (unsigned int jj=0; jj<ewMatch_val.size2(); ++jj) {
+ lowMatchTemp_val(ii,jj) = collinearEWMatch_val(ii,jj)*ewMatch_val(ii,jj);
+ }
+ }
+ // perform all the multiplications
+ matrix<Complex> temp(highRunning_val.size1(),collinearHighRunning_val.size2());
+ axpy_prod(highRunning_val,collinearHighRunning_val,temp);
+ matrix<Complex> temp2(collinearLowRunning_val.size1(),lowRunning_val.size2());
+ axpy_prod(collinearLowRunning_val,lowRunning_val,temp2);
+ matrix<Complex> temp3(temp2.size1(),lowMatchTemp_val.size2());
+ axpy_prod(temp2,lowMatchTemp_val,temp3);
+ temp2.resize(temp3.size1(),temp.size2());
+ axpy_prod(temp3,temp,temp2);
+ matrix<complex<InvEnergy2> > result(temp2.size1(),highMatch_val.size2());
+ axpy_prod_local(temp2,highMatch_val,result);
+ // for(unsigned int ix=0;ix<result.size1();++ix) {
+ // for(unsigned int iy=0;iy<result.size2();++iy) {
+ // cerr << s*result(ix,iy) << " ";
+ // }
+ // cerr << "\n";
+ // }
+ // reset the couplings
+ coupling()->SU3(SU3save);
+ coupling()-> EW( EWsave);
+ // return the answer
+ return result;
+}
+
+
+double ElectroWeakReweighter::reweightggqqbar() const {
+ // momenta and invariants
+ Lorentz5Momentum p1 = subProcess()->incoming().first ->momentum();
+ Lorentz5Momentum p2 = subProcess()->incoming().second->momentum();
+ Lorentz5Momentum p3 = subProcess()->outgoing()[0]->momentum();
+ Lorentz5Momentum p4 = subProcess()->outgoing()[1]->momentum();
+ tcPDPtr qbar = subProcess()->outgoing()[0]->dataPtr();
+ tcPDPtr q = subProcess()->outgoing()[1]->dataPtr();
+ if(q->id()<0) {
+ swap(p3,p4 );
+ swap(q ,qbar);
+ }
+ Energy2 s = (p1+p2).m2();
+ Energy2 t = (p1-p4).m2();
+ Energy2 u = (p1-p3).m2();
+ // boost to partonic rest frame and rescale momenta of outgoing
+ // so zero mass
+ Lorentz5Momentum psum=p1+p2;
+ LorentzRotation boost(-psum.boostVector());
+ p1 *= boost;
+ p2 *= boost;
+ p3 *= boost;
+ p4 *= boost;
+ p3.setMass(ZERO);
+ p3.rescaleRho();
+ p4.setMass(ZERO);
+ p4.rescaleRho();
+ // LO matrix element coefficents
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ bornQQGGweights,bornRRGGweights;
+ // quark left doublet
+ if(q->id()<5) {
+ bornQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,true);
+ }
+ else {
+ bornQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,true);
+ }
+ // quark right singlet
+ if(q->id()==0) {
+ if(q->id()==6)
+ bornRRGGweights = evaluateRunning(EWProcess::tRtRGG,s,t,u,true);
+ else
+ bornRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,true);
+ }
+ else
+ bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true);
+ // EW corrected matrix element coefficients
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ EWQQGGweights,EWRRGGweights;
+ // quark left doublet
+ if(q->id()<5) {
+ EWQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,false);
+ }
+ else {
+ EWQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,false);
+ }
+ // quark right singlet
+ if(q->id()%2==0) {
+ if(q->id()==6)
+ EWRRGGweights = evaluateRunning(EWProcess::tRtRGG,s,t,u,false);
+ else
+ EWRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,false);
+ }
+ else
+ EWRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,false);
+ SpinorWaveFunction qw(p4,qbar,incoming);
+ SpinorBarWaveFunction qbarw(p3,q ,incoming);
+ vector<LorentzVector<Complex> > eps1,eps2;
+ SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps1,eps2,1);
+ boost::numeric::ublas::matrix<Complex>
+ bornME = boost::numeric::ublas::zero_matrix<Complex>(3,3),
+ EWME = boost::numeric::ublas::zero_matrix<Complex>(3,3);
+ // helicities of outgoing quarks
+ for(unsigned int iq=0;iq<2;++iq) {
+ if(iq==0) {
+ qw.reset (0);
+ qbarw.reset(1);
+ }
+ else {
+ qw.reset (1);
+ qbarw.reset(0);
+ }
+ LorentzVector<complex<Energy> > current = iq==0 ?
+ qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) :
+ qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave());
+ for(unsigned int i1=0;i1<2;++i1) {
+ complex<Energy> d31 = eps1[i1].dot(p3);
+ for(unsigned int i2=0;i2<2;++i2) {
+ boost::numeric::ublas::vector<complex<Energy2> > M(5);
+ Complex d34 = eps1[i1].dot(eps2[i2]);
+ complex<Energy> d42 = eps2[i2].dot(p4);
+ // M0 in paper
+ M(0) = qw.dimensionedWave().slash(eps1[i1])
+ .slash(p2-p4).vectorCurrent(qbarw.dimensionedWave()).dot(eps2[i2]);
+ // M4 in paper
+ M(2) = current.dot(eps2[i2])*d31;
+ // M5 in paper
+ M(3) = -current.dot(eps1[i1])*d42;
+ // M1 in paper (missing factor)
+ M(1) = current.dot(p2);
+ // M6 in paper
+ M(4) = M(1)*d31*d42/GeV2;
+ // M1 final factor
+ M(1) *= d34;
+ // coefficient of different contributions
+ boost::numeric::ublas::vector<Complex> Cborn(3),CEW(3);
+ if(iq==0) {
+ axpy_prod_local(bornQQGGweights,M,Cborn);
+ axpy_prod_local(EWQQGGweights ,M,CEW );
+ }
+ else {
+ axpy_prod_local(bornRRGGweights,M,Cborn);
+ axpy_prod_local(EWRRGGweights ,M,CEW );
+ }
+ for(unsigned int ix=0;ix<3;++ix) {
+ for(unsigned int iy=0;iy<3;++iy) {
+ bornME(ix,iy) += Cborn(ix)*conj(Cborn(iy));
+ EWME (ix,iy) += CEW (ix)*conj(CEW (iy));
+ }
+ }
+ }
+ }
+ }
+ double born = 24.*real(bornME(0,0))+20./3.*real(bornME(1,1))+12.*real(bornME(2,2));
+ double EW = 24.*real(EWME(0,0))+20./3.*real(EWME(1,1))+12.*real(EWME(2,2));
+ return EW/born;
+}
diff --git a/MatrixElement/EW/ElectroWeakReweighter.h b/MatrixElement/EW/ElectroWeakReweighter.h
--- a/MatrixElement/EW/ElectroWeakReweighter.h
+++ b/MatrixElement/EW/ElectroWeakReweighter.h
@@ -1,140 +1,164 @@
// -*- C++ -*-
#ifndef Herwig_ElectroWeakReweighter_H
#define Herwig_ElectroWeakReweighter_H
//
// This is the declaration of the ElectroWeakReweighter class.
//
#include "ThePEG/MatrixElement/ReweightBase.h"
#include "EWCouplings.h"
#include "CollinearSudakov.h"
#include "SoftSudakov.h"
namespace Herwig {
using namespace ThePEG;
/**
* The ElectroWeakReweighter class.
*
* @see \ref ElectroWeakReweighterInterfaces "The interfaces"
* defined for ElectroWeakReweighter.
*/
class ElectroWeakReweighter: public ReweightBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
ElectroWeakReweighter();
/**
* The destructor.
*/
virtual ~ElectroWeakReweighter();
//@}
public:
/**
* Return the weight for the kinematical configuation provided by
* the assigned XComb object (in the LastXCombInfo base class).
*/
virtual double weight() const;
/**
*
*/
static tEWCouplingsPtr coupling() {
assert(staticEWCouplings_);
return staticEWCouplings_;
}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
+ * Functions to reweight specific processes
+ */
+ //@{
+ /**
+ * Reweight \f$g g\to q\bar{q}\f$
+ */
+ double reweightggqqbar() const;
+
+ /**
+ * Reweight \f$q\bar{q}\to g g\f$
+ */
+ double reweightqqbargg() const;
+ //@}
+
+protected:
+
+ /**
* Check the evolution for a fixed s,t,u
*/
void testEvolution(Energy2 s,Energy2 t, Energy2 u) const;
+ /**
+ * Evalaute the running
+ */
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ evaluateRunning(EWProcess::Process process, Energy2 s,
+ Energy2 t, Energy2 u, bool born) const;
+
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ElectroWeakReweighter & operator=(const ElectroWeakReweighter &);
private:
/**
* The Electroweak Couplings
*/
EWCouplingsPtr EWCouplings_;
/**
* The Collinear Sudakov
*/
CollinearSudakovPtr collinearSudakov_;
/**
* The Soft Sudakov
*/
SoftSudakovPtr softSudakov_;
/**
* The couplings to allow global access
*/
static tEWCouplingsPtr staticEWCouplings_;
};
}
#endif /* Herwig_ElectroWeakReweighter_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 3:35 PM (1 d, 6 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4017482
Default Alt Text
(29 KB)

Event Timeline