C++ Vec4D类代码示例

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);
}

 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;
 }

void Six_Particle_Mass::Evaluate(const Vec4D& mom1,const Vec4D& mom2,
                                 const Vec4D& mom3,const Vec4D& mom4,
                                 const Vec4D& mom5,const Vec4D& mom6,
                                 double weight, double ncount)
{
  Vec4D  p = mom1+mom2+mom3+mom4+mom5+mom6;
  p_histo->Insert(p.Mass(),weight,ncount);
}

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);
}

void EtaTracksAsym::Evaluate(const Vec4D & mom,double weight,double ncount) 
{
  
  double eta = 0.;
  eta=mom.Eta();
  
  if (eta<0.) {
    p_histo->Insert(-eta,-weight,ncount); 
  }
  else {
    p_histo->Insert(eta,weight,ncount); 
  }
} 

double PHASIC::ComputePhi
(Vec4D pijt,Vec4D pkt,Vec4D pi,const int mode)
{
  LN_Pair ln(GetLN((mode&1)?-pijt:pijt,(mode&2)?-pkt:pkt,0));
  Vec4D ktt(0.0,cross(Vec3D(ln.m_l),Vec3D(ln.m_n)));
  Poincare cms(ln.SaneCMS());
  cms.Boost(ln.m_l);
  cms.Boost(pi);
  Poincare zax(ln.m_l,Vec4D::ZVEC);
  if (mode!=3 && ktt.PSpat2()>rpa->gen.SqrtAccu()) {
    zax.Rotate(ktt);
  }
  else {
    msg_Debugging()<<"Set fixed n_perp\n";
    if (IsZero(ln.m_l.PPerp2())) ln.m_l[1]=ln.m_l[2]=0.0;
    zax=Poincare(ln.m_l,Vec4D::ZVEC);
    ktt=Vec4D(0.0,1.0,1.0,0.0);
  }
  zax.Rotate(pi);
  Poincare xax(ktt,Vec4D::XVEC);
  xax.Rotate(pi);
  return pi.Phi();
}

double EV_SC_ET::Calc(const Vec4D &mom)
{
  return mom.EPerp();
}

void Three_Particle_Correlation::
Evaluate(const Vec4D & mom1,const Vec4D & mom2,const Vec4D & mom3,double weight, double ncount)
{ 
  double eta1 = mom1.Eta();
  double eta2 = mom2.Eta();
  double eta3 = mom3.Eta();
  if (ATOOLS::dabs(eta1)>2. || ATOOLS::dabs(eta2)>2. || ATOOLS::dabs(eta3)>2.) return;
  double phi1 = mom1.Phi();
  double phi2 = mom2.Phi();
  double phi3 = mom3.Phi();

  double delta12 = ATOOLS::dabs(phi1-phi2); delta12 = ATOOLS::Min(delta12,2.*M_PI-delta12);
  double delta23 = ATOOLS::dabs(phi2-phi3); delta23 = ATOOLS::Min(delta23,2.*M_PI-delta23); 
  double delta31 = ATOOLS::dabs(phi3-phi1); delta31 = ATOOLS::Min(delta31,2.*M_PI-delta31); 

  if (ATOOLS::dabs(eta1-eta2)>2.) p_histo->Insert(delta12,weight,ncount);
  if (ATOOLS::dabs(eta2-eta3)>2.) p_histo->Insert(delta23,weight,ncount);
  if (ATOOLS::dabs(eta3-eta1)>2.) p_histo->Insert(delta31,weight,ncount);

  if (mom1.PPerp()>10. || mom2.PPerp()>10. || mom3.PPerp()>10.) {
    if (ATOOLS::dabs(eta1-eta2)>2.) p_histo10->Insert(delta12,weight,ncount);
    if (ATOOLS::dabs(eta2-eta3)>2.) p_histo10->Insert(delta23,weight,ncount);
    if (ATOOLS::dabs(eta3-eta1)>2.) p_histo10->Insert(delta31,weight,ncount);
  }

  if (mom1.PPerp()>20. || mom2.PPerp()>20. || mom3.PPerp()>20.) {
    if (ATOOLS::dabs(eta1-eta2)>2.) p_histo20->Insert(delta12,weight,ncount);
    if (ATOOLS::dabs(eta2-eta3)>2.) p_histo20->Insert(delta23,weight,ncount);
    if (ATOOLS::dabs(eta3-eta1)>2.) p_histo20->Insert(delta31,weight,ncount);
  }

  if (mom1.PPerp()>50. || mom2.PPerp()>50. || mom3.PPerp()>50.) {
    if (ATOOLS::dabs(eta1-eta2)>2.) p_histo50->Insert(delta12,weight,ncount);
    if (ATOOLS::dabs(eta2-eta3)>2.) p_histo50->Insert(delta23,weight,ncount);
    if (ATOOLS::dabs(eta3-eta1)>2.) p_histo50->Insert(delta31,weight,ncount);
  }

  if (mom1.PPerp()>100. || mom2.PPerp()>100. || mom3.PPerp()>100.) {
    if (ATOOLS::dabs(eta1-eta2)>2.) p_histo100->Insert(delta12,weight,ncount);
    if (ATOOLS::dabs(eta2-eta3)>2.) p_histo100->Insert(delta23,weight,ncount);
    if (ATOOLS::dabs(eta3-eta1)>2.) p_histo100->Insert(delta31,weight,ncount);
  }
} 

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;
  }
}

 static bool OrderV(const Vec4D &a,const Vec4D &b) 
 { return dabs(a.Eta())>dabs(b.Eta()); }

	void CopyVector(const Vec4D & vector, GLfixed *params) {
		params[0] = vector.x();
		params[1] = vector.y();
		params[2] = vector.z();
		params[3] = vector.w();
	}

void One_Particle_BeamAngle::EvaluateNLOcontrib(const Vec4D & mom,double weight, double ncount) 
{
  double ct = mom.CosTheta();
  p_histo->InsertMCB(ct,weight,ncount); 
} 

void One_Particle_M::Evaluate(const Vec4D & mom,double weight, double ncount) 
{
  double mass = mom.Mass();
  p_histo->Insert(mass,weight,ncount); 
}