#include "CSSHOWER++/Main/CS_Cluster_Definitions.H" #include "CSSHOWER++/Showers/Shower.H" #include "PHASIC++/Channels/CSS_Kinematics.H" #include "ATOOLS/Org/Exception.H" #include "ATOOLS/Org/Run_Parameter.H" #include "ATOOLS/Org/MyStrStream.H" #include "ATOOLS/Org/Shell_Tools.H" #include "ATOOLS/Org/My_Limits.H" using namespace CSSHOWER; using namespace PHASIC; using namespace PDF; using namespace ATOOLS; CS_Cluster_Definitions::CS_Cluster_Definitions (Shower *const shower,const int kmode,const int pdfcheck,const int kfmode): p_shower(shower), m_kmode(kmode), m_pdfcheck(pdfcheck), m_kfmode(kfmode) {} Cluster_Param CS_Cluster_Definitions::Cluster(const Cluster_Config &ca) { DEBUG_FUNC(ca); CS_Parameters cs(KT2(ca.p_ampl, ca.p_ampl->Leg(ca.m_i), ca.p_ampl->Leg(ca.m_j), ca.p_ampl->Leg(ca.m_k), ca.m_mo,ca.p_ms,ca.m_kin,ca.m_mode|m_kmode)); return Cluster_Param(this,cs.m_wk,cs.m_kt2,cs.m_mu2,cs.m_cpl, cs.m_kin,cs.m_kmode,cs.m_pijt,cs.m_pkt,cs.m_lt); } Splitting_Function_Base *CS_Cluster_Definitions::GetSF (const ATOOLS::Cluster_Leg *i,const ATOOLS::Cluster_Leg *j, const ATOOLS::Cluster_Leg *k,const ATOOLS::Flavour &mo,CS_Parameters &cs) const { const SF_EEE_Map *cmap(&p_shower->GetSudakov()->FFFMap()); if (cs.m_mode==2) cmap=&p_shower->GetSudakov()->FFIMap(); else if (cs.m_mode==1) cmap=cs.m_col>=0?&p_shower->GetSudakov()->IFFMap(): &p_shower->GetSudakov()->FIFMap(); else if (cs.m_mode==3) cmap=cs.m_col>=0?&p_shower->GetSudakov()->IFIMap(): &p_shower->GetSudakov()->FIIMap(); SF_EEE_Map::const_iterator eees(cmap->find(ProperFlav(i->Flav(),i->Id()&3))); if (eees==cmap->end()) { msg_Debugging()<<"No splitting function (i)\n"; return NULL; } SF_EE_Map::const_iterator ees(eees->second.find(ProperFlav(j->Flav(),j->Id()&3))); if (ees==eees->second.end()) { msg_Debugging()<<"No splitting function (j)\n"; return NULL; } SF_E_Map::const_iterator es(ees->second.find(ProperFlav(mo,(i->Id()&3)|(j->Id()&3)))); if (es==ees->second.end()) { msg_Debugging()<<"No splitting function (ij)\n"; return NULL; } return es->second; } CS_Parameters CS_Cluster_Definitions::KT2 (const ATOOLS::Cluster_Amplitude *ampl, const ATOOLS::Cluster_Leg *i,const ATOOLS::Cluster_Leg *j, const ATOOLS::Cluster_Leg *k,const ATOOLS::Flavour &mo, const ATOOLS::Mass_Selector *const ms,const int ikin, const int kmode,const int force) { p_ms=ms; int kin(ikin<0?p_shower->KinScheme():ikin), col(1); if ((i->Id()&3)<(j->Id()&3)) { std::swap(i,j); col=-1; } p_b=ampl->Leg(i==ampl->Leg(0)?1:0); const Vec4D pi(i->Mom()), pj(j->Mom()), pk(k->Mom()); const double Q2=(pi+pj+pk).Abs2(); double mb2=p_b->Mom().Abs2(), mfb2=p_ms->Mass2(p_b->Flav()); if (mfb2==0.0 || IsEqual(mb2,mfb2,1.0e-6)) mb2=mfb2; double mi2=pi.Abs2(), mfi2=p_ms->Mass2(i->Flav()); double mj2=pj.Abs2(), mfj2=p_ms->Mass2(j->Flav()); double mk2=pk.Abs2(), mfk2=p_ms->Mass2(k->Flav()); if ((mfi2==0.0 && IsZero(mi2,1.0e-6)) || IsEqual(mi2,mfi2,1.0e-6)) mi2=mfi2; if ((mfj2==0.0 && IsZero(mj2,1.0e-6)) || IsEqual(mj2,mfj2,1.0e-6)) mj2=mfj2; if ((mfk2==0.0 && IsZero(mk2,1.0e-6)) || IsEqual(mk2,mfk2,1.0e-6)) mk2=mfk2; double mij2=p_ms->Mass2(mo); if (kmode&1) { mij2=(pi+pj).Abs2(); kin=0; } Kin_Args lt; CS_Parameters cs(sqrt(std::numeric_limits::max()), 1.0,1.0,0.0,0.0,0.0, ((i->Id()&3)?1:0)|((k->Id()&3)?2:0),kin,kmode&1); cs.m_kt2=-1.0; cs.m_wk=0.0; cs.m_col=col; if ((i->Id()&3)==0) { if ((j->Id()&3)==0) { if ((k->Id()&3)==0) { lt=ClusterFFDipole(mi2,mj2,mij2,mk2,pi,pj,pk,1|(kin?4:0)); if (lt.m_stat!=1) if (!force) return cs; double kt2=p_shower->KinFF()->GetKT2(Q2,lt.m_y,lt.m_z,mi2,mj2,mk2,mo,j->Flav()); cs=CS_Parameters(kt2,lt.m_z,lt.m_y,lt.m_phi,1.0,Q2,0,kin,kmode&1); cs.m_pkt=cs.m_pk=lt.m_pk; cs.m_pijt=lt.m_pi; } else { lt=ClusterFIDipole(mi2,mj2,mij2,mk2,pi,pj,-pk,1|8|(kin?4:0)); if (lt.m_pk[3]*pk[3]>0.0 || lt.m_pk[0]<0.0 || lt.m_stat!=1) if (!force) return cs; double kt2=p_shower->KinFI()->GetKT2(Q2,1.0-lt.m_y,lt.m_z,mi2,mj2,mk2,mo,j->Flav()); cs=CS_Parameters(kt2,lt.m_z,lt.m_y,lt.m_phi,1.0-lt.m_y,Q2,2,kin,kmode&1); cs.m_pkt=-(cs.m_pk=lt.m_pk); cs.m_pijt=lt.m_pi; } } } else { if ((j->Id()&3)==0) { if ((k->Id()&3)==0) { lt=ClusterIFDipole(mi2,mj2,mij2,mk2,mb2,-pi,pj,pk,-p_b->Mom(),3|(kin?4:0)); if ((kmode&1) && lt.m_mode) lt.m_stat=-1; if (lt.m_pi[3]*pi[3]>0.0 || lt.m_pi[0]<0.0 || lt.m_z<0.0 || lt.m_stat!=1) if (!force) return cs; double kt2=p_shower->KinIF()->GetKT2(Q2,lt.m_y,lt.m_z,mi2,mj2,mk2,mo,j->Flav()); cs=CS_Parameters(kt2,lt.m_z,lt.m_y,lt.m_phi,lt.m_z,Q2,1,lt.m_mode,kmode&1); cs.m_pkt=cs.m_pk=lt.m_pk; cs.m_pijt=-lt.m_pi; cs.m_lt=lt.m_lam; } else { lt=ClusterIIDipole(mi2,mj2,mij2,mk2,-pi,pj,-pk,3|(kin?4:0)); if (lt.m_pi[3]*pi[3]>0.0 || lt.m_pi[0]<0.0 || lt.m_z<0.0 || lt.m_stat!=1) if (!force) return cs; double kt2=p_shower->KinII()->GetKT2(Q2,lt.m_y,lt.m_z,mi2,mj2,mk2,mo,j->Flav()); cs=CS_Parameters(kt2,lt.m_z,lt.m_y,lt.m_phi,lt.m_z,Q2,3,kin,kmode&1); cs.m_pkt=-(cs.m_pk=lt.m_pk); cs.m_pijt=-lt.m_pi; cs.m_lt=lt.m_lam; } } } cs.m_col=col; KernelWeight(i,j,k,mo,cs,kmode); return cs; } double CS_Cluster_Definitions::GetX(const Cluster_Leg* l, Splitting_Function_Base* const sf) const { if ((l->Id() & 3) == 0) THROW(fatal_error, "Invalid call for parton ID="+ToString(l->Id())); int beam((l->Id() & 1) ? 0 : 1); if (sf) sf->Lorentz()->SetBeam(beam); Vec4D mom( (l->Mom()[0] < 0.0) ? -l->Mom() : l->Mom() ); return p_shower->ISR()->CalcX(mom); } Flavour CS_Cluster_Definitions::ProperFlav(const Flavour &fl,const int anti) const { Flavour pfl(fl); switch (pfl.Kfcode()) { case kf_gluon_qgc: pfl=Flavour(kf_gluon); break; default: break; } return anti?pfl.Bar():pfl; } void CS_Cluster_Definitions::KernelWeight (const ATOOLS::Cluster_Leg *i,const ATOOLS::Cluster_Leg *j, const ATOOLS::Cluster_Leg *k,const ATOOLS::Flavour &mo, CS_Parameters &cs,const int kmode) const { Splitting_Function_Base *cdip(GetSF(i,j,k,mo,cs)); if (cdip==NULL || !cdip->On()) { cs.m_wk=0.0; return; } if (m_pdfcheck && (cs.m_mode&1)) { int beam=i->Id()&1?0:1; if (p_shower->ISR()->PDF(beam) && !p_shower->ISR()->PDF(beam)->Contains(mo)) { msg_Debugging()<<"Not in PDF: "<Flav(),k->Id()&3)); if (!(kmode&32) && !cdip->Coupling()->AllowSpec(fls,1)) { msg_Debugging()<<"Invalid spectator "<Mom()+j->Mom()+k->Mom()).Abs2()); if (Q2<(cs.m_mode&1? p_shower->GetSudakov()->ISPT2Min(): p_shower->GetSudakov()->FSPT2Min())) { msg_Debugging()<<"Small Q2 "<SetMS(p_ms); cdip->SetFlavourSpec(fls); cs.m_mu2=cdip->Lorentz()->Scale(cs.m_z,cs.m_y,cs.m_kt2,Q2); cs.m_mu2=Max(cs.m_mu2,cs.m_mode&1? p_shower->GetSudakov()->ISPT2Min(): p_shower->GetSudakov()->FSPT2Min()); if (m_kfmode && !(kmode&16)) cs.m_mu2*=cdip->Coupling()->CplFac(cs.m_mu2); cs.m_cpl=cdip->PureQCD()?0:-1; if (!(kmode&2)) return; else { double eta=1.0; if (cs.m_mode==1) eta=GetX(i,cdip)*cs.m_z; else if (cs.m_mode==2) eta=GetX(k,cdip)*(1.0-cs.m_y); else if (cs.m_mode==3) eta=GetX(i,cdip)*cs.m_z; cs.m_wk=(*cdip)(cs.m_z,cs.m_y,eta,-1.0,Q2)*Q2/cs.m_kt2; if (IsBad(cs.m_wk)) cs.m_wk=0.0; SF_Lorentz *lf = cdip->Lorentz(); SF_Coupling *cf = cdip->Coupling(); msg_Debugging()<<"Kernel weight (NLO="<<(kmode&16) <<") [m="< w = "<j) std::swap(i,j); Vec4D_Vector after(ampl.Legs().size()-1); double mb2(0.0); if (i<2) { mb2=ampl.Leg(1-i)->Mom().Abs2(); double mfb2(p_ms->Mass2(ampl.Leg(1-i)->Flav())); if (mfb2==0.0 || IsEqual(mb2,mfb2,1.0e-6)) mb2=mfb2; } Vec4D pi(ampl.Leg(i)->Mom()), pj(ampl.Leg(j)->Mom()); Vec4D pk(ampl.Leg(k)->Mom()), pb(i<2?ampl.Leg(1-i)->Mom():Vec4D()); double mi2=pi.Abs2(), mfi2=p_ms->Mass2(ampl.Leg(i)->Flav()); double mj2=pj.Abs2(), mfj2=p_ms->Mass2(ampl.Leg(j)->Flav()); double mk2=pk.Abs2(), mfk2=p_ms->Mass2(ampl.Leg(k)->Flav()); if ((mfi2==0.0 && IsZero(mi2,1.0e-6)) || IsEqual(mi2,mfi2,1.0e-6)) mi2=mfi2; if ((mfj2==0.0 && IsZero(mj2,1.0e-6)) || IsEqual(mj2,mfj2,1.0e-6)) mj2=mfj2; if ((mfk2==0.0 && IsZero(mk2,1.0e-6)) || IsEqual(mk2,mfk2,1.0e-6)) mk2=mfk2; double mij2=p_ms->Mass2(mo); bool sk(true); if (kmode&1) { mij2=(pi+pj).Abs2(); kin=0; } Kin_Args lt; if (i>1) { if (k>1) lt=ClusterFFDipole(mi2,mj2,mij2,mk2,pi,pj,pk,2|(kin?4:0)); else lt=ClusterFIDipole(mi2,mj2,mij2,mk2,pi,pj,-pk,2|(kin?4:0)); if (lt.m_pk[3]*(k>1?pk[3]:-pk[3])<0.0 || lt.m_pk[0]<0.0) return Vec4D_Vector(); } else { if (k>1) { lt=ClusterIFDipole(mi2,mj2,mij2,mk2,mb2,-pi,pj,pk,-pb,2|(kin?4:0)); if ((kmode&1) && lt.m_mode) lt.m_stat=-1; } else lt=ClusterIIDipole(mi2,mj2,mij2,mk2,-pi,pj,-pk,2|(kin?4:0)); if (lt.m_pi[3]*pi[3]>0.0 || lt.m_pi[0]<0.0) return Vec4D_Vector(); } if (lt.m_stat<0) return Vec4D_Vector(); for (size_t l(0), m(0);m1?lt.m_pi:-lt.m_pi; else if (m==(size_t)k && sk) after[l]=k>1?lt.m_pk:-lt.m_pk; else after[l]=lt.m_lam*ampl.Leg(m)->Mom(); ++l; } if (ampl.Next()) ampl.Next()->SetKin(lt.m_mode&1); return after; } namespace CSSHOWER { std::ostream &operator<<(std::ostream &str,const CS_Parameters &cs) { return str<<"CS{kt="<