Commit fb801f17 authored by Stefan Hoeche's avatar Stefan Hoeche

svn merge -r28982:29147 ^/branches/dire

parent f5c7d5fa
......@@ -38,13 +38,13 @@ namespace DIRE {
double Integral(const Splitting &s) const
{
double I=20.0/9.0*0.5*0.5*log((s.m_Q2+s.m_t0)/(s.m_Q2*sqr(s.m_eta)+s.m_t0));
return I*p_sk->GF()->CplMax(s)/(2.0*M_PI)*m_jmax;
return I*p_sk->GF()->CplMax(s)/(2.0*M_PI)*m_jmax*PDFEstimate(s);
}
double Estimate(const Splitting &s) const
{
double E=20.0/9.0*0.5*s.m_z/(sqr(s.m_z)+s.m_t0/s.m_Q2);
return E*p_sk->GF()->CplMax(s)/(2.0*M_PI)*m_jmax;
return E*p_sk->GF()->CplMax(s)/(2.0*M_PI)*m_jmax*PDFEstimate(s);
}
bool GeneratePoint(Splitting &s) const
......@@ -52,29 +52,37 @@ namespace DIRE {
double k2(s.m_t0/s.m_Q2);
s.m_z=sqrt(pow((1.0+k2)/(sqr(s.m_eta)+k2),-ran->Get())*(1.0+k2)-k2);
s.m_phi=2.0*M_PI*ran->Get();
do s.m_z2=pow(s.m_z,ran->Get());
while (1.0+sqr(1.0-s.m_z2)<2.0*ran->Get());
do s.m_z2=1.0/(1.0+ran->Get()*(1.0/s.m_z-1.0));
while (1.0-s.m_z2+sqr(s.m_z2)/2.0<ran->Get());
s.m_phi2=2.*M_PI*ran->Get();
return true;
}
bool Compute(Splitting &s,const int mode) const
{
s.m_x=s.m_z;
double xa(s.m_z2), za(s.m_z), Q2(s.m_Q2);
s.m_y=s.m_s=0.0;
s.m_x=xa+s.m_t*za/(Q2*xa);
s.m_mk2=Q2*(xa/za-1.0)+s.m_t/xa-s.m_mi2-s.m_ml2;
if (!Lorentz_IF::Compute(s,0)) return false;
Splitting c(s);
double amax(p_sk->GF()->CplMax(s));
double I(amax/(2.0*M_PI)*4.0/3.0*
(2.0*log(1.0/s.m_z)-(3.0-s.m_z)*(1.0-s.m_z)/2.0));
(2.0/s.m_z-s.m_z+2.0*log(s.m_z)-1.0));
for (s.m_s=s.m_Q2*pow(ran->Get(),Max(1.0/I,1.0e-3));
s.m_s>s.m_t0;s.m_s*=pow(ran->Get(),Max(1.0/I,1.0e-3))) {
s.m_y=s.m_s*s.m_x/s.m_Q2;
s.m_y=(c.m_t=s.m_s)*za/Q2;
s.m_x=xa+s.m_t*za/(Q2*xa)+s.m_s*za/Q2;
s.m_mk2=Q2*(xa/za-1.0)+s.m_s+s.m_t/xa-s.m_mi2-s.m_ml2;
if (p_sk->GF()->Coupling(c)/amax<ran->Get()) continue;
if (Lorentz_IF::Compute(s,0)) return true;
if (Lorentz_IF::Compute(s,0)) {
s.m_mk2=p_ms->Mass2(s.p_s->Flav());
return true;
}
}
s.m_mk2=p_ms->Mass2(s.p_s->Flav());
s.m_y=s.m_s=0.0;
return Lorentz_IF::Compute(s,0);
return true;
}
};// end of class FFFF_IF
......
......@@ -38,13 +38,13 @@ namespace DIRE {
double Integral(const Splitting &s) const
{
double I=20.0/9.0*0.5*0.5*log((s.m_Q2+s.m_t0)/(s.m_Q2*sqr(s.m_eta)+s.m_t0));
return I*p_sk->GF()->CplMax(s)/(2.0*M_PI)*m_jmax;
return I*p_sk->GF()->CplMax(s)/(2.0*M_PI)*m_jmax*PDFEstimate(s);
}
double Estimate(const Splitting &s) const
{
double E=20.0/9.0*0.5*s.m_z/(sqr(s.m_z)+s.m_t0/s.m_Q2);
return E*p_sk->GF()->CplMax(s)/(2.0*M_PI)*m_jmax;
return E*p_sk->GF()->CplMax(s)/(2.0*M_PI)*m_jmax*PDFEstimate(s);
}
bool GeneratePoint(Splitting &s) const
......@@ -52,29 +52,37 @@ namespace DIRE {
double k2(s.m_t0/s.m_Q2);
s.m_z=sqrt(pow((1.0+k2)/(sqr(s.m_eta)+k2),-ran->Get())*(1.0+k2)-k2);
s.m_phi=2.0*M_PI*ran->Get();
do s.m_z2=pow(s.m_z,ran->Get());
while (1.0+sqr(1.0-s.m_z2)<2.0*ran->Get());
do s.m_z2=1.0/(1.0+ran->Get()*(1.0/s.m_z-1.0));
while (1.0-s.m_z2+sqr(s.m_z2)/2.0<ran->Get());
s.m_phi2=2.*M_PI*ran->Get();
return true;
}
bool Compute(Splitting &s,const int mode) const
{
s.m_x=s.m_z;
s.m_y=s.m_t/s.m_Q2*sqr(s.m_z2/s.m_x);
s.m_s=s.m_z2*(p_ms->Mass(m_fl[1])-p_ms->Mass(m_fl[3])/(1.0-s.m_z2));
if (!Lorentz_II::Compute(s,0)) return false;
double xa(s.m_z2), za(s.m_z);
double Q2((s.p_c->Mom()+s.p_s->Mom()).Abs2()-s.m_mij2-s.m_mk2);
s.m_s=xa*(p_ms->Mass(m_fl[1])-p_ms->Mass(m_fl[3])/(1.0-xa));
s.m_x=(Q2+s.m_mij2-s.m_mj2+s.m_s)/(xa/za*Q2);
s.m_y=s.m_t/Q2*sqr(xa/za);
Splitting c(s);
c.m_x=xa-(c.m_y=(s.m_s+s.m_mi2+s.m_ml2)/(Q2/za));
if (!(Lorentz_II::Compute(s,0) &&
Lorentz_II::Compute(c,0))) return false;
double amax(p_sk->GF()->CplMax(s));
double I(amax/(2.0*M_PI)*4.0/3.0*
(2.0*log(1.0/s.m_z)-(3.0-s.m_z)*(1.0-s.m_z)/2.0));
for (s.m_s=s.m_Q2*pow(ran->Get(),Max(1.0/I,1.0e-3));
s.m_s>s.m_t0;s.m_s*=pow(ran->Get(),Max(1.0/I,1.0e-3))) {
(2.0/s.m_z-s.m_z+2.0*log(s.m_z)-1.0));
for (c.m_t=s.m_s=s.m_Q2*pow(ran->Get(),Max(1.0/I,1.0e-3));
s.m_s>s.m_t0;c.m_t=s.m_s*=pow(ran->Get(),Max(1.0/I,1.0e-3))) {
if (p_sk->GF()->Coupling(c)/amax<ran->Get()) continue;
if (Lorentz_II::Compute(s,0)) return true;
s.m_x=(Q2+s.m_mij2-s.m_mj2+s.m_s)/(xa/za*Q2);
c.m_x=xa-(c.m_y=(s.m_s+s.m_mi2+s.m_ml2)/(Q2/za));
if (Lorentz_II::Compute(s,0) &&
Lorentz_II::Compute(c,0)) return true;
}
s.m_s=s.m_z2*(p_ms->Mass(m_fl[1])-p_ms->Mass(m_fl[3])/(1.0-s.m_z2));
return Lorentz_II::Compute(s,0);
s.m_s=xa*(p_ms->Mass(m_fl[1])-p_ms->Mass(m_fl[3])/(1.0-xa));
s.m_x=(Q2+s.m_mij2-s.m_mj2+s.m_s)/(xa/za*Q2);
return true;
}
};// end of class FFFF_II
......
......@@ -111,7 +111,7 @@ bool Lorentz::SetLimits(Splitting &s) const
if (m_fl.size()>3) s.m_ml2=p_ms->Mass2(m_fl[3]);
s.m_mk2=p_ms->Mass2(s.p_s->Flav());
s.m_Q2=dabs((s.p_c->Mom()+s.p_s->Mom()).Abs2()
-s.m_mi2-s.m_mj2-s.m_mk2);
-s.m_mi2-s.m_ml2-s.m_mj2-s.m_mk2);
s.m_eta=s.p_c->GetXB();
return true;
}
......@@ -27,7 +27,7 @@ double Lorentz_FI::Jacobian(const Splitting &s) const
double Lorentz_FI::PDFEstimate(const Splitting &s) const
{
return pow(1.0e6,s.p_s->GetXB()*s.m_t0/Min(s.m_t1,s.m_Q2));
return 1.0;
}
int Lorentz_FI::Construct(Splitting &s,const int mode) const
......
......@@ -29,14 +29,15 @@ double Lorentz_IF::PDFEstimate(const Splitting &s) const
(s.m_eta,Min(s.m_t1,s.m_Q2),m_fl[0],s.p_c->Beam()-1);
double fn=p_sk->PS()->GetXPDF
(s.m_eta,Min(s.m_t1,s.m_Q2),m_fl[1],s.p_c->Beam()-1);
if (m_fl[1]==Flavour(kf_u).Bar()||m_fl[1]==Flavour(kf_d).Bar()) {
double fo0=p_sk->PS()->GetXPDF(s.m_eta,s.m_t0,m_fl[0],s.p_c->Beam()-1);
double fn0=p_sk->PS()->GetXPDF(0.2,s.m_t0,m_fl[1],s.p_c->Beam()-1);
if (m_fl[1].Mass(true)<1.0 && m_fl[0].Mass(true)>=1.0) {
double tcut(Max(s.m_t0,sqr(2.0*m_fl[0].Mass(true))));
double fo0=p_sk->PS()->GetXPDF(s.m_eta,tcut,m_fl[0],s.p_c->Beam()-1);
double fn0=p_sk->PS()->GetXPDF(0.2,tcut,m_fl[1],s.p_c->Beam()-1);
if (fo0 && dabs(fo0)<dabs(fo)) fo=fo0;
if (dabs(fn0)>dabs(fn)) fn=fn0;
}
if (fo==0.0) return 0.0;
return dabs(fn/fo)*pow(1.0e6,s.m_eta*s.m_t0/Min(s.m_t1,s.m_Q2));
if (fo<p_sk->PS()->PDFMin()) return 0.0;
return dabs(fn/fo);
}
int Lorentz_IF::Construct(Splitting &s,const int mode) const
......@@ -72,13 +73,15 @@ bool Lorentz_IF::Cluster(Splitting &s,const int mode) const
bool Lorentz_IF::Compute(Splitting &s,const int mode) const
{
s.m_y=s.m_t/s.m_Q2/(1.0-s.m_z);
s.m_x=s.m_z;
if (mode) {
s.m_y=s.m_t/s.m_Q2/(1.0-s.m_z);
s.m_x=s.m_z;
}
if (s.m_mk2==0.0)
return s.m_y>0.0 && s.m_y<1.0;
double muk2(s.m_mk2*s.m_z/s.m_Q2);
double zp((1.0-s.m_z)/(1.0-s.m_z+muk2));
return s.m_y>0.0 && s.m_y<zp;
return s.m_y>=0.0 && s.m_y<1.0 && s.m_x<1.0;
double muk2(s.m_mk2*s.m_x/(s.m_Q2+s.m_mk2));
double zp((1.0-s.m_x)/(1.0-s.m_x+muk2));
return s.m_y>=0.0 && s.m_y<zp && s.m_x<1.0;
}
Lorentz_IF_123::Lorentz_IF_123(const Kernel_Key &k):
......@@ -86,18 +89,25 @@ Lorentz_IF_123::Lorentz_IF_123(const Kernel_Key &k):
{
}
double Lorentz_IF_123::Jacobian(const Splitting &s) const
double Lorentz_IF_123::J(const Splitting &s) const
{
double za(s.m_x), xa(s.m_z2);
double Q2(s.m_mij2+s.m_mk2-(s.p_c->Mom()+s.p_s->Mom()).Abs2());
double za(s.m_z), xa(s.m_z2), Q2(s.m_Q2);
double saijk(Q2*xa/za+s.m_t/xa), sjk(saijk-Q2+s.m_s);
sjk-=s.m_mi2+s.m_ml2+s.m_mj2+s.m_mk2;
double ej((sjk+s.m_mj2-s.m_mk2)/(2.0*sjk));
double bai(sqrt(1.0+4.0*s.m_s*sjk/sqr(saijk)));
double bj(sqrt(1.0-4.0*s.m_mj2*sjk/sqr(s.m_mj2+sjk)));
double J(Q2*s.m_t/ej+za*sjk/sqr(bai)*(1.0+2.0*s.m_s/saijk)
double tjk(sjk+s.m_mj2+s.m_mk2), tai(-s.m_s+s.m_mi2+s.m_ml2);
double J(Q2*s.m_t/ej+za*tjk/sqr(bai)*(1.0-2.0*tai/saijk)
*(1.0-s.m_t/xa/saijk/ej)*(Q2*xa/za-s.m_t/xa));
J/=bai*bj*xa*za*sqr(saijk);
return J*Lorentz_IF::Jacobian(s);
return J/(bai*bj*xa*za*sqr(saijk));
}
double Lorentz_IF_123::Jacobian(const Splitting &s) const
{
Splitting c(s);
c.m_x=s.m_z;
return s.m_z2*J(s)*Lorentz_IF::Jacobian(c);
}
int Lorentz_IF_123::Construct(Splitting &s,const int mode) const
......@@ -110,12 +120,11 @@ int Lorentz_IF_123::Construct(Splitting &s,const int mode) const
break;
}
if (s.m_s<rpa->gen.SqrtAccu()) s.m_s=0.0;
double xa(s.m_z2), za(s.m_x);
double Q2(s.m_mij2+s.m_mk2-(s.p_c->Mom()+s.p_s->Mom()).Abs2());
double xa(s.m_z2), za(s.m_z), Q2(s.m_Q2);
Kin_Args ff(s.m_s*za/Q2,xa+s.m_s*za/Q2+s.m_t*za/(Q2*xa),s.m_phi,s.m_kin);
if (ff.m_z>1.0) return -1.0;
ff.m_mk2=Q2*(xa/za-1.0)+s.m_s+s.m_t/xa;
if (ff.m_mk2<sqr(sqrt(s.m_mj2)+sqrt(s.m_mk2))) return -1.0;
if (ff.m_z>1.0) return -1;
ff.m_mk2=Q2*(xa/za-1.0)+s.m_s+s.m_t/xa-s.m_mi2-s.m_ml2;
if (ff.m_mk2<sqr(sqrt(s.m_mj2)+sqrt(s.m_mk2))) return -1;
if (ConstructIFDipole
(s.m_mi2,s.m_ml2,s.m_mij2,s.m_mk2,b?p_ms->Mass2(b->Flav()):0.0,
-s.p_c->Mom(),s.p_s->Mom(),b?-b->Mom():Vec4D(),ff)<0)
......
......@@ -27,6 +27,8 @@ namespace DIRE {
class Lorentz_IF_123: public Lorentz_IF {
protected:
double J(const Splitting &s) const;
double Jacobian(const Splitting &s) const;
public:
......
......@@ -29,14 +29,15 @@ double Lorentz_II::PDFEstimate(const Splitting &s) const
(s.m_eta,Min(s.m_t1,s.m_Q2),m_fl[0],s.p_c->Beam()-1);
double fn=p_sk->PS()->GetXPDF
(s.m_eta,Min(s.m_t1,s.m_Q2),m_fl[1],s.p_c->Beam()-1);
if (m_fl[1]==Flavour(kf_u).Bar()||m_fl[1]==Flavour(kf_d).Bar()) {
double fo0=p_sk->PS()->GetXPDF(s.m_eta,s.m_t0,m_fl[0],s.p_c->Beam()-1);
double fn0=p_sk->PS()->GetXPDF(0.2,s.m_t0,m_fl[1],s.p_c->Beam()-1);
if (m_fl[1].Mass(true)<1.0 && m_fl[0].Mass(true)>=1.0) {
double tcut(Max(s.m_t0,sqr(2.0*m_fl[0].Mass(true))));
double fo0=p_sk->PS()->GetXPDF(s.m_eta,tcut,m_fl[0],s.p_c->Beam()-1);
double fn0=p_sk->PS()->GetXPDF(0.2,tcut,m_fl[1],s.p_c->Beam()-1);
if (fo0 && dabs(fo0)<dabs(fo)) fo=fo0;
if (dabs(fn0)>dabs(fn)) fn=fn0;
}
if (fo==0.0) return 0.0;
return dabs(fn/fo)*pow(1.0e6,s.m_eta*s.m_t0/Min(s.m_t1,s.m_Q2));
if (fo<p_sk->PS()->PDFMin()) return 0.0;
return dabs(fn/fo);
}
int Lorentz_II::Construct(Splitting &s,const int mode) const
......@@ -79,14 +80,16 @@ Lorentz_II_123::Lorentz_II_123(const Kernel_Key &k):
double Lorentz_II_123::Jacobian(const Splitting &s) const
{
return Lorentz_II::Jacobian(s);
Splitting c(s);
c.m_x=s.m_z;
return Lorentz_II::Jacobian(c);
}
int Lorentz_II_123::Construct(Splitting &s,const int mode) const
{
if (s.m_s<rpa->gen.SqrtAccu()) s.m_s=0.0;
double Q2((s.p_c->Mom()+s.p_s->Mom()).Abs2()-s.m_mij2-s.m_mk2);
double q2(Q2+s.m_mij2-s.m_mj2+s.m_s), xa(s.m_z2), za(s.m_x);
double q2(Q2+s.m_mij2-s.m_mj2+s.m_s), xa(s.m_z2), za(s.m_z);
Kin_Args ff(s.m_t/Q2*sqr(xa/za),q2/(xa/za*Q2),s.m_phi,s.m_kin);
if (ConstructIIDipole
(-s.m_s,s.m_mj2,s.m_mij2,s.m_mk2,
......
......@@ -232,7 +232,7 @@ Splitting Shower::GeneratePoint(Parton &p,const double &t)
if (sum==0.0) return Splitting();
win=Splitting(&p,NULL,ct);
}
win.m_t*=exp(log(ran->Get())*Max(2.0*M_PI/sum,1.0e-3));
win.m_t*=exp(log(ran->Get())*2.0*M_PI/sum);
if (win.m_t<m_tmin[p.Beam()?1:0]) return win;
double disc(sum*ran->Get()), csum(0.0);
for (size_t j(0);j<splits.size();++j)
......
......@@ -441,13 +441,18 @@ bool HepMC2_Interface::Sherpa2HepMC(ATOOLS::Blob_List *const blobs,
m_blob2genvertex.clear();
m_particle2genparticle.clear();
HepMC::GenVertex * vertex;
HepMC::GenVertex * vertex, * psvertex(NULL);
std::vector<HepMC::GenParticle*> beamparticles;
for (ATOOLS::Blob_List::iterator blit=blobs->begin();
blit!=blobs->end();++blit) {
// Call the Blob to vertex code, changes vertex pointer above
if (Sherpa2HepMC(*(blit),vertex)) {
event.add_vertex(vertex);
for (size_t i(0);i<(*blit)->NInP();++i)
if ((*blit)->InParticle(i)->ProductionBlob()==NULL) {
psvertex=vertex;
break;
}
if ((*blit)->Type()==ATOOLS::btp::Signal_Process) {
if ((**blit)["NLO_subeventlist"]) {
THROW(fatal_error,"Events containing correlated subtraction events"
......@@ -470,6 +475,17 @@ bool HepMC2_Interface::Sherpa2HepMC(ATOOLS::Blob_List *const blobs,
}
}
} // End Blob_List loop
if (beamparticles.empty()) {
Vec4D pbeam[2]={rpa->gen.PBeam(0),rpa->gen.PBeam(1)};
HepMC::FourVector pa(pbeam[0][1],pbeam[0][2],pbeam[0][3],pbeam[0][0]);
HepMC::FourVector pb(pbeam[1][1],pbeam[1][2],pbeam[1][3],pbeam[1][0]);
HepMC::GenParticle *inpart[2] = {
new HepMC::GenParticle(pa,(long int)rpa->gen.Beam1(),2),
new HepMC::GenParticle(pb,(long int)rpa->gen.Beam2(),2)};
psvertex->add_particle_in(inpart[0]);
psvertex->add_particle_in(inpart[1]);
event.set_beam_particles(inpart[0],inpart[1]);
}
if (beamparticles.size()==2) {
event.set_beam_particles(beamparticles[0],beamparticles[1]);
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment