C++ Vec4D::Abs2方法代码示例

double Vrel(Vec4D p1,Vec4D& p2)
{
  double m21=(p1+p2).Abs2();
  double m1=p1.Abs2();
  double m2=p2.Abs2();
  return sqrt(Lambda(m21,m1,m2))/(m21-m1-m2);
}

Complex W_To_Lepton_Neutrino::InfraredSubtractedME_1_05(unsigned int i) {
  m_moms       = m_moms1[i];                // set to set of momenta to be used
  Vec4C epsW   = Polarization_Vector(m_moms[0])[m_spins[0]];
  Vec4C epsP   = conj(Polarization_Vector(m_moms[3])[m_spins[3]]);
  Vec4D q      = m_moms[1]+m_moms[3];       // fermion propagator momenta
  double q2    = q.Abs2();
  Vec4D Q      = m_moms[0]-m_moms[3];       // boson propagator momenta
  double Q2    = Q.Abs2();
  double m     = m_masses[1];               // fermion mass/propagator pole
  double m2    = sqr(m);
  double M     = m_masses[0];               // boson mass/propagator pole
  double M2    = sqr(M);
  m_moms[4]    = m_moms[5] = q;             // enter those into m_moms
  m_flavs[4]   = m_flavs[1];                // set to corresponding particle/antiparticle
  m_flavs[5]   = m_flavs[1].Bar();
  XYZFunc XYZ(6,m_moms,m_flavs,false);
  m_flavs[4] = m_flavs[5] = Flavour(kf_none);
  // two diagrams
  // M_1 = -ie^2/(2sqrt(2)sW) * 1/((pl+k)^2-m^2)
  //       * ubar(l)gamma^mu(-pl-k+m)gamma^nu P_L v(nu) eps_nu^W eps_mu^y*
  // M_2 = ie/(2sqrt(2)sW) * 1/(pW-k)^2-M^2)
  //       * ubar(l)gamma_rho P_L v(nu)
  //       * [-2g^{rho,nu}pW^mu + g^{rho,mu}pW^nu
  //          + g^{nu,mu}k^rho + 1/M^2(pW-k)^rho pW^nu pW^mu]
  //       * eps_nu^W eps_mu^y*
  Complex r1 = Complex(0.,0.);
  Complex r2 = Complex(0.,0.);
  Complex r3 = Complex(0.,0.);
  Lorentz_Ten3C ten31,ten32,ten33,ten34;
  for (unsigned int s=0; s<=1; s++) {
    r1 += XYZ.X(1,m_spins[1],epsP,4,s,1.,1.)
          *XYZ.X(4,s,epsW,2,m_spins[2],m_cR,m_cL);
    r2 += XYZ.X(1,m_spins[1],epsP,5,s,1.,1.)
          *XYZ.X(5,s,epsW,2,m_spins[2],m_cR,m_cL);
  }
  Vec4D p = m_moms[0];
  Vec4D k = m_moms[3];
  // index ordering rho(1),nu(2),mu(3)
  // -2g^{rho,nu}pW^mu
  ten31 = BuildTensor(MetricTensor(),-2.*p);
  // g^{rho,mu}pW^nu
  ten32 = BuildTensor(MetricTensor(),p).Transpose(2,3);
  // g^{nu,mu}k^rho
  ten33 = BuildTensor(MetricTensor(),k).Transpose(1,3);
  // 1/M^2(pW-k)^rho pW^nu pW^mu
  ten34 = -1./M2*BuildTensor(p-k,p,p);
  Lorentz_Ten3C ten = ten31+ten32+ten33+ten34;
  // v^\sigma = L^\sigma\mu\lambda epsW_\mu epsP_\lambda
  Vec4C v3 = Contraction(Contraction(ten,3,epsP),2,epsW);
  r3 = XYZ.X(1,m_spins[1],v3,2,m_spins[2],m_cR,m_cL);

  r1 *= (1.+m/sqrt(q2))/(q2-m2);
  r2 *= (1.-m/sqrt(q2))/(q2-m2);
  r3 *= -1./(Q2-M2);

  return (m_i*m_e*m_e)/(2.*m_sqrt2*m_sW)*(r1+r2/*+r3*/);
}

LN_Pair PHASIC::GetLN
(const Vec4D &pi,const Vec4D &pk,const int mode)
{
  double mi2(pi.Abs2()), mk2(pk.Abs2());
  double eps(pi*pk), kap(eps*eps-mi2*mk2);
  if (kap<0.0) return LN_Pair();
  kap=Sign(eps)*sqrt(kap);
  Vec4D l(((eps+kap)*pi-mi2*pk)/(2.0*kap));
  Vec4D n(((eps+kap)*pk-mk2*pi)/(2.0*kap));
  return LN_Pair(l,n,mode);
}

void Kinematics_Base::SetFixVec(Parton *const p,Vec4D mom,
				const Kin_Args &lt,const int mode) const
{
  if (p->GetNext()) SetFixVec(p->GetNext(),mom,lt,mode|4);
  if (p->FixSpec()==Vec4D()) return;
  Vec4D oldp(p->OldMomentum()), ref(p->FixSpec());
  if ((mode&3)==3 || ((mode&1)==1 && lt.m_mode==0)) {
    if (mode&4) {
      Poincare_Sequence lam(lt.m_lam);
      lam.Invert();
      mom=lam*mom;
    }
    else {
      oldp=lt.m_lam*oldp;
      ref=lt.m_lam*ref;
    }
  }
  if (IsEqual(oldp,mom,rpa->gen.SqrtAccu())) {
    p->SetFixSpec(ref);
    p->SetOldMomentum(oldp);
    return;
  }
  Vec4D np(0.0,cross(Vec3D(oldp),Vec3D(mom)));
  if (np.PSpat2()<=1.0e-6) {
    msg_Debugging()<<"Set fixed n_perp\n";
    np=Vec4D(0.0,cross(Vec3D(oldp),Vec3D(1.0,1.0,0.0)));
  }
  np*=1.0/np.PSpat();
  Vec4D lp(0.0,cross(Vec3D(oldp),Vec3D(np)));
  lp*=1.0/lp.PSpat();
  Vec4D pl(0.0,(Vec3D(ref)*Vec3D(lp))*lp);
  Vec4D pn(0.0,(Vec3D(ref)*Vec3D(np))*np);
  double D(oldp*ref), T(oldp.PSpat()), F(ref[0]);
  double Q(mom[0]), P(mom.PSpat()), S(mom.Abs2());
  Poincare rot(oldp,mom);
  if (oldp.Abs2()>1.0e-3 && mom.Abs2()>1.0e-3) {
    Poincare oldcms(oldp), newcms(mom);
    oldcms.Boost(ref);
    rot.Rotate(ref);
    newcms.BoostBack(ref);
  }
  else {
    double E((Q*D+P/T*(F*S-Q*D))/S);
    ref=Vec4D(E,Vec3D(mom)*(Q*E-D)/(P*P));
    ref+=pn+rot*pl;
  }
  p->SetFixSpec(ref);
  p->SetOldMomentum(mom);
}

 double pT2pythia(Cluster_Amplitude* ampl, const Cluster_Leg& RadAfterBranch,
             const Cluster_Leg& EmtAfterBranch,
             const Cluster_Leg& RecAfterBranch, int ShowerType){
   // Save type: 1 = FSR pT definition, else ISR definition
   int Type   = ShowerType;
   // Calculate virtuality of splitting
   int sign = (Type==1) ? 1 : -1;
   Vec4D Q(RadAfterBranch.Mom() + sign*EmtAfterBranch.Mom());
   double Qsq = sign * Q.Abs2();
   // Mass term of radiator
   DEBUG_VAR(ampl->MS());
   double m2Rad = ( abs(RadAfterBranch.Flav().Kfcode()) >= 4
                && abs(RadAfterBranch.Flav().Kfcode()) < 7)
                ? ampl->MS()->Mass2(RadAfterBranch.Flav())
                : 0.;
   // Construct 2->3 variables for FSR
   Vec4D  sum     = RadAfterBranch.Mom() + RecAfterBranch.Mom()
                  + EmtAfterBranch.Mom();
   double m2Dip = sum.Abs2();
   double x1 = 2. * (sum * RadAfterBranch.Mom()) / m2Dip;
   double x3 = 2. * (sum * EmtAfterBranch.Mom()) / m2Dip;
   // Construct momenta of dipole before/after splitting for ISR 
   Vec4D qBR(RadAfterBranch.Mom() - EmtAfterBranch.Mom() + RecAfterBranch.Mom());
   Vec4D qAR(RadAfterBranch.Mom() + RecAfterBranch.Mom()); 
   // Calculate z of splitting, different for FSR and ISR
   double z = (Type==1) ? x1/(x1+x3)
                      : (qBR.Abs2())/( qAR.Abs2());
   // Separation of splitting, different for FSR and ISR
   double pTpyth = (Type==1) ? z*(1.-z) : (1.-z);
   // pT^2 = separation*virtuality
   pTpyth *= (Qsq - sign*m2Rad);
   if(pTpyth < 0.) pTpyth = 0.;
   // Return pT
   return pTpyth;
 }

PDF::CParam Default_Core_Scale::Calculate(Cluster_Amplitude *const ampl)
{
  DEBUG_FUNC("");
  msg_Debugging()<<*ampl<<"\n";
  if (ampl->Legs().size()==3 && ampl->NIn()==2) {
    double kt2cmin(ampl->Leg(2)->Mom().Abs2());
    return PDF::CParam(kt2cmin,kt2cmin,0.0,kt2cmin,-1);
  }
  double muf2(0.0), mur2(0.0), muq2(0.0);
  Cluster_Amplitude *campl(Cluster(ampl->Copy()));
  if (campl->Legs().size()!=ampl->Legs().size())
    msg_Debugging()<<*campl<<"\n";
  if (campl->Legs().size()!=4) {
    double q2((campl->Leg(0)->Mom()+campl->Leg(1)->Mom()).Abs2());
    Vec4D ewsum;
    for (size_t i(0);i<campl->Legs().size();++i)
      if (!campl->Leg(i)->Flav().Strong()) ewsum+=campl->Leg(i)->Mom();
    if (ewsum==Vec4D()) ewsum=campl->Leg(0)->Mom()+campl->Leg(1)->Mom();
    if (campl->NIn()==2 &&
	campl->Leg(0)->Flav().Strong() &&
	campl->Leg(1)->Flav().Strong()) {// HThat'/2
      q2=ewsum.PPerp();
      for (size_t i(0);i<campl->Legs().size();++i)
	if (campl->Leg(i)->Flav().Strong())
	  q2+=campl->Leg(i)->Mom().PPerp();
      q2=sqr(ewsum.Mass()+q2/2.0);
    }
    campl->Delete();
    return PDF::CParam(q2,dabs(ewsum.Abs2()),0.0,q2,-1);
  }
  Flavour_Vector fl; fl.resize(4);
  fl[0]=campl->Leg(0)->Flav();
  fl[1]=campl->Leg(1)->Flav();
  fl[2]=campl->Leg(2)->Flav();
  fl[3]=campl->Leg(3)->Flav();
  if (fl[0].Strong() && fl[1].Strong()) {// hh collision
    if (fl[2].Strong() && fl[3].Strong()) {
      msg_Debugging()<<"pure QCD like\n";
      double s(2.0*campl->Leg(0)->Mom()*campl->Leg(1)->Mom());
      double t(2.0*campl->Leg(0)->Mom()*campl->Leg(2)->Mom());
      double u(2.0*campl->Leg(0)->Mom()*campl->Leg(3)->Mom());
      muq2=muf2=mur2=-1.0/(1.0/s+1.0/t+1.0/u)/4.0;
    }
    else if (!fl[2].Strong() && !fl[3].Strong()) {
      msg_Debugging()<<"DY like\n";
      muq2=muf2=mur2=(campl->Leg(0)->Mom()+campl->Leg(1)->Mom()).Abs2();
    }
    else if (fl[2].Strong() && !fl[3].Strong()) {
      msg_Debugging()<<"jV like\n";
      muq2=muf2=mur2=campl->Leg(3)->Mom().MPerp2()/4.0;
    }
    else if (!fl[2].Strong() && fl[3].Strong()) {
      msg_Debugging()<<"Vj like\n";
      muq2=muf2=mur2=campl->Leg(2)->Mom().MPerp2()/4.0;
    }
    else THROW(fatal_error,"Internal error.");
  }
  else if (!fl[0].Strong() && !fl[1].Strong()) {// ll collision
    if (fl[2].Strong() && fl[3].Strong()) {
      msg_Debugging()<<"jets like\n";
      muq2=muf2=mur2=(campl->Leg(0)->Mom()+campl->Leg(1)->Mom()).Abs2();
    }
  }
  else {// lh collision
    msg_Debugging()<<"DIS like\n";
    muq2=muf2=mur2=dabs((campl->Leg(fl[0].Strong()?1:0)->Mom()+
			 campl->Leg(fl[2].Strong()?3:2)->Mom()).Abs2());
  }
  campl->Delete();
  msg_Debugging()<<"\\mu_f = "<<sqrt(muf2)<<"\n"
		 <<"\\mu_r = "<<sqrt(mur2)<<"\n"
		 <<"\\mu_q = "<<sqrt(muq2)<<"\n";
  return PDF::CParam(muf2,muq2,0.0,mur2,-1);
}

void II_DipoleSplitting::SetMomenta(const Vec4D *mom)
{
  m_mom.clear();
  for(int i=0;i<=m_m;i++) m_mom.push_back(mom[i]);

  m_pi = mom[m_i];
  m_pj = mom[m_j];
  m_pk = mom[m_k];

  m_xijk = (m_pk*m_pi-m_pi*m_pj-m_pj*m_pk)/(m_pk*m_pi);

  m_ptk  = m_pk;
  m_ptij = m_xijk*m_pi;

  Vec4D K  = m_pi-m_pj+m_pk;
  Vec4D Kt = m_ptij+m_pk;
  Vec4D KKt = K+Kt;
  for(int i=2;i<=m_m;i++) m_mom[i]-=2.*(m_mom[i]*KKt/KKt.Abs2()*KKt-m_mom[i]*K/K.Abs2()*Kt);

  m_vi   = (m_pi*m_pj)/(m_pi*m_pk);
  m_a = m_vi;

  m_Q2 = (-m_pi+m_pj-m_pk).Abs2();
  if (m_es==0) {
    m_kt2 = m_Q2*(1.-m_xijk)/m_xijk*m_vi;
  }
  else {
  m_kt2 = m_Q2/m_xijk*m_vi;
  switch (m_ft) {
  case 1:
    m_kt2*=(1.-m_xijk);
    break;
  case 4:
    m_kt2*=(1.-m_xijk);
    break;
  }
  }

//   m_pt1  =    m_pj;
//   m_pt2  =-1.*m_vi*m_pk;
  m_pt1  =    m_pj-m_vi*m_pk;
  m_pt2  =    m_ptij;

  switch (m_ft) {
  case 1:
    m_sff = 2./(1.-m_xijk)-(1.+m_xijk);
    m_av  = m_sff;
    break;
  case 2:
    m_sff = 1.-2.*m_xijk*(1.-m_xijk);
    m_av  = m_sff;
    break;
  case 3:
    m_sff = m_xijk;
    m_av  = m_sff + 2.0*(1.0-m_xijk)/m_xijk;
    break;
  case 4:
    m_sff = m_xijk/(1.-m_xijk)+m_xijk*(1.-m_xijk);
    m_av  = m_sff + (1.0-m_xijk)/m_xijk;
  }
}